Thermodynamics and steady state of quantum motors and pumps far from equilibrium
Raúl A. Bustos-Marún Hernan L. Calvo





[Electronic  address: ]rbustos@famaf.unc.edu.ar
Instituto de Física Enrique Gaviola (CONICET) and FaMAF, Universidad Nacional de Córdoba, Argentina
Facultad de Ciencias Químicas, Universidad Nacional de Córdoba, Argentina


Instituto de Física Enrique Gaviola (CONICET) and FaMAF, Universidad Nacional de Córdoba, Argentina
Departamento de Física, Universidad Nacional de Río Cuarto, Ruta 36, Km 601, 5800 Río Cuarto, Argentina


In this article, we briefly review dynamical and thermodynamical aspects of different forms of quantum motors and quantum pumps. We then extend previous results to provide new theoretical tools for a systematic study of those phenomena at far-from-equilibrium conditions.
We mainly focus on two key topics: (1) The steady-state regime of quantum motors and pumps, paying particular attention to the role of higher-order terms in the nonadiabatic expansion of the current-induced forces. (2) The thermodynamical properties of such systems, emphasizing systematic ways of studying the relationship between different energy fluxes (charge and heat currents, and mechanical power) passing through the system when beyond-first-order expansions are required. We derive a general order-by-order scheme based on energy conservation to rationalize how every order of the expansion of one form of energy flux is connected with the others. We use this approach to give a physical interpretation of the leading terms of the expansion. Finally, we illustrate the above-discussed topics in a double quantum dot within the Coulomb-blockade regime and capacitively coupled to a mechanical rotor. We find many exciting features of this system for arbitrary nonequilibrium conditions: A definite parity of the expansion coefficients with respect to the voltage or temperature biases; negative friction coefficients; and the fact that, under fixed parameters, the device can exhibit multiple steady states where it may operate as a quantum motor or as a quantum pump depending on the initial conditions.




Thermodynamics and steady state of quantum motors and pumps far from equilibrium
    Hernán L. Calvo
    December 30, 2023
================================================================================





§ INTRODUCTION

In recent years, there has been a sustained growth in the interest in different forms of nanomachines. This was boosted by seminal experiments <cit.>,
the blooming of new theoretical proposals <cit.>,
 and the latest developments towards the understanding of the fundamental physics underlying such systems <cit.>.
Quantum mechanics has proven to be crucial in the description of a broad family of nanomachines which can be put together under the generic name of “quantum motors” and “quantum pumps” <cit.>. They typically consist in an electromechanical device connected to electronic reservoirs and controlled by nonequilibrium sources, see Figure <ref>. These nonequilibrium sources may include temperature gradients and bias voltages among the reservoirs, or even an external driving of the internal parameters of the system. The dimensions of the electronic component of these devices are normally within the characteristic coherence length of the electrons flowing through them, hence the essential role of quantum mechanics in their description.

Various aspects of quantum motors and pumps have been extensively studied in the literature. For example, it has been shown that quantum interferences can be exploited to boost the performance of these devices. Remarkably, some systems operate solely due to quantum interference, e.g., quantum pumps and motors based on chaotic quantum dots <cit.>, Thouless quantum pumps and motors <cit.>, or Anderson quantum motors <cit.>, among others.
On the other hand, the strong Coulomb repulsion between electrons in quantum-dot-based pumps and motors has shown to enhance the performance (or even induce the activation) of these nanodevices <cit.>. 
The effect of decoherence has also been addressed <cit.> as well as the influence of the friction forces and the system-lead coupling in the dynamics of quantum motors and pumps <cit.>.
Indeed, the thermodynamics of those systems has proven to be a key aspect to study.
In the last years, different individual efforts have coalesced to give rise to a new field dubbed “quantum thermodynamics” <cit.>, which studies the relations among the different energy fluxes that drive the motion of those machines where quantum mechanics plays a fundamental role.



Despite the progress in the theoretical description of quantum motors and pumps, most of the research has focused on parameter conditions that lie close to the thermodynamic equilibrium, i.e., small bias voltages, temperature gradients, or frequencies of the external driving <cit.>.
This is reasonable since under such conditions the linear response regime of the nonequilibrium sources gives an accurate description of the problem, greatly simplifying its general treatment.
For example, in this regime it is common to define dimensionless figures of merit made by some combination of linear response coefficients which give a measurement of the efficiency or the maximum power that quantum devices can achieve.
It is also known that such figures of merit fail in nonlinear regime conditions <cit.>.
Although efforts have been made in this direction, currently there is not a nonlinear version for the figures of merit and the performance must be calculated from the microscopic details of the system's dynamics.
One strategy to deal with such situations is to use phenomenological models where the linear response coefficients are parameterized with respect to the voltage biases, the temperature gradients, or to other relevant parameters of the system, see <cit.> and references therein. However, such parameterizations usually hide the physics behind the nonlinearities and require the optimization of a number of variables that grows very fast with the complexity of the model.
The weakly nonlinear regime of transport has also been explored within the scattering matrix formalism. Under these conditions, it is enough to expand the response coefficients up to second order of the voltage biases and temperature gradients, which can be done by using standard quantum transport techniques. This approach has been applied to a variety of situations, where it proved to be a valuable strategy, see <cit.> and references therein. However, it would be also important to extend this method to more general situations without hindering the description of the physical processes that take part in the nonlinear effects, while keeping the deep connexions between the response coefficients.

Regarding the dynamics of quantum motors and pumps, one can notice that most of the works in the literature assume a constant terminal velocity of the driving parameters without a concrete model for them.
A typical problem is that when these devices are coupled to nonequilibrium sources, nonconservative current-induced forces (CIFs) appear. These CIFs come, in first place, by assuming a type of Born-Oppenheimer approximation where the electronic and mechanical degrees of freedom can be treated separately, and, secondly, by evaluating the mean value of the force operator <cit.>. Because of the delayed response of the electronic degrees of freedom to the mechanical motion, one should include the so-called nonadiabatic corrections to the CIFs. This phenomenon is translated into a possible complex dependency of the CIFs on the velocity of the mechanical degrees of freedom. When the effect of the mechanical velocities on CIFs can be treated in linear response, it is clear whether is adequate or not to assume a constant terminal velocity <cit.>.
However, in far-from-equilibrium conditions this subject has not been fully addressed.

In this article, we discuss two key aspects of far-from-equilibrium quantum motors and pumps: Their steady-state dynamics, specially when CIFs present nontrivial dependencies on the terminal velocities; and their nonequilibrium thermodynamical properties, when a linear response description is not enough. We provide a systematic expansion to study the 
relations between the different energy fluxes that drive the quantum device. These aspects are illustrated in a concrete example where nonlinear effects due to nonequilibrium sources play a major role on the steady-state properties of the system. We show that these nonlinearities may result in, e.g., negative friction coefficients or motor/pump coexistence regimes.

Our work is organized as follows: In Section <ref> 
we present the general model that describes the considered type of systems and we derive an effective Langevin equation that characterizes the dynamics of the mechanical degrees of freedom, treated classically in the present context.
In Section <ref>, we discuss in general terms the steady-state dynamics of quantum motors and pumps, highlighting some key aspects that differentiate close-to and far-from equilibrium conditions.
In Section <ref> we derive, on general grounds, the first law of thermodynamics for the kind of systems treated. Then we expand the different energy fluxes passing through the system in terms of the nonequilibrium sources (temperature gradients, bias voltages, and velocities) for arbitrary number of reservoirs.
In this way, we obtain an order-by-order relation between the different energy fluxes entering and leaving the device.
In Section <ref>, we perform a derivation of the rate of entropy production from first principles. Then, based on the second law of thermodynamics, we discuss the limits of the efficiency for different forms of quantum motors and pumps in general nonequilibrium conditions.
In Section <ref>, we analyze and give physical interpretation to some of the relations obtained in Section <ref>.
Finally, in Section <ref> we consider the CIFs for strongly interacting electrons in a particular example based on a double quantum dot system coupled to a mechanical rotor. We then analyze in detail the effects of higher-order terms in the CIFs on the final steady state of the electromechanical system. 



§ CURRENT-INDUCED FORCES AND LANGEVIN EQUATION
 

In this section we introduce the generic model for the treatment of CIFs and the standard method employed
in the description of the dynamics of the mechanical degrees of freedom. As starting point, we consider as local system the region where  electronic and mechanical degrees of freedom are present and coupled to each other, like the examples shown in Figure <ref>.
Such a local system is generically modeled by the following Hamiltonian

    Ĥ_local = Ĥ_s(X̂)+P̂^2/2m_eff+U(X̂,t),

where X̂=(X̂_1,...,X̂_N) is the vector of mechanical coordinates and P̂=(P̂_1,...,P̂_N) collects their associated momenta, m_eff is the effective mass related to X̂, and U(X̂,t) represents some external potential, of mechanical nature, that may be acting on the local system.
The explicit time-dependence on this potential thus emphasizes the fact that an external agent can exert some effective work on the local system.
The Hamiltonian Ĥ_s includes both the electronic degrees of freedom and their coupling to the mechanical ones through

    Ĥ_s(X̂) =  ∑_i E_i(X̂)|i⟩⟨i|,

where the sum runs over all possible electronic many-body eigenstates |i⟩.
The local system is then coupled to macroscopic reservoirs and the total Hamiltonian, including the mechanical degrees of freedom, reads

    Ĥ_total = Ĥ_local+∑_r Ĥ_r+∑_r Ĥ_s,r.

Each lead r is described as a reservoir of noninteracting electrons through the Hamiltonian

    Ĥ_r = ∑_kσϵ_rkĉ_rkσ^†ĉ_rkσ,

where ĉ_rkσ^† (ĉ_rkσ) creates (annihilates) an electron in the r-reservoir with state-index k and spin projection σ.
As usual, the reservoirs are assumed to be always in equilibrium, characterized by a temperature T_r and electrochemical potential μ_r.
The coupling between the local system and the r-lead is determined by the tunnel Hamiltonian

    Ĥ_s,r = ∑_kσℓ( t_rℓd̂_ℓσ^†ĉ_rkσ+h.c.),

where t_rℓ denotes the tunnel amplitude, assumed to be k and σ independent for simplicity, and the fermion operator
d̂_ℓσ^† (d̂_ℓσ) creates (annihilates) one
electron with spin σ in the ℓ-orbital of the local system. The tunnel-coupling strengths Γ_rℓ = 2πρ_r |t_rℓ|^2 then characterize the rate at which the electrons enter/leave the local system from/to the r-reservoir, where ρ_r is the density of states in the r-lead.
Note that Ĥ_s was defined in the eigenstate basis while Ĥ_s,r is written in terms of single-particle field operators.
The tunnel matrix elements accounting for transitions between different eigenstates can then be obtained from linear superpositions of the above tunnel amplitudes <cit.>.

To obtain an effective description of the dynamics of the mechanical degrees of freedom through a Langevin equation, we start from the Heisenberg equation of motion for the P̂ operator, which yields

    m_effẊ̂̇t+∇̂ ̂Û(X̂,t) = -∇̂ ̂Ĥ_s(X̂).

The measured value A_measured of an observable described by an operator Â can always be taken as its mean value A= ⟨Â|$⟩ plus some fluctuationξ_Aaround it, i.e.,A_measured = A +ξ_A.
We will work under the nonequilibrium Born-Oppenheimer approximation
 <cit.>
(or Ehrenfest approximation <cit.>),
where the dynamics of the electronic and mechanical degrees of freedom
can be separated and the latter is treated classically. This
allows us to neglect the fluctuations of the terms appearing in the left-hand
side of Equation (<ref>), and describe the mechanical motion only through the mean valueX, which is reasonable for large or massive objects. With this in mind, we obtain the following Langevin equation of motion:

    m_effẊt+F_ext = F+ξ,

whereF = -⟨∇̂Ĥ_s⟩= i ⟨[Ĥ_s(X̂),P̂]⟩andξaccount for the mean value and the fluctuation of
the CIF, respectively.[Throughout this manuscript we take ħ=1 for simplicity.]
As we shall see later on, the external force applied to the mechanical part of the local system,F_ext, plays the role of an eventual “load” force for a quantum motor or a “driving” force for a quantum pump.
As this force will be typically opposed to the CIF, we defineF_extwith a minus sign for better clarity in future discussions.
The main task, therefore, relies on the calculation of the CIFs from appropriate formalisms capable to describe the dynamics of the electronic part of the system.
Once these forces are calculated, we can use Equation (<ref>) to integrate the classical equations of motion and 
obtain the effective dynamics of the complete electromechanical system.

In most previous worksFis expanded up to first order inẊ, i.e.,F ≈F^(0) - γ ·Ẋ.
The resulting CIF is then the sum of an adiabatic contributionF^(0)and its first nonadiabatic correctionF^(1) = -γẊ, respectively. Under this approximation, Equation (<ref>) turns into

    m_effẊt+F_ext = F^(0)-γ·Ẋ+ξ.

Explicit formulas for the calculation ofF^(0),γ, andξin terms of Green functions and scattering matrices were derived
in <cit.>, and extended in <cit.> to account for decoherent events.
Although these expressions were obtained in the context of noninteracting particles, they can be used in effective Hamiltonians derived from first principles calculations <cit.>.
In <cit.> the CIFs were obtained from the Floquet Green's function formalism.
The role of Coulomb interactions was addressed through different formalisms and methods like, e.g., many-body perturbation theory based on nonequilibrium Green's functions <cit.>; modeling the system as a Luttinger liquid <cit.>; and using a time-dependent slave-boson approximation <cit.>.
In <cit.>, explicit expressions for the CIFs within the Coulomb blockade regime of transport were obtained using a real-time diagrammatic approach <cit.>, which we present in more detail in Section <ref> when considering the example of Figure <ref>c.



§ MECHANICAL STEADY-STATE 


In this paper, we will restrict ourselves to systems that perform overall cyclic motions.
Immediate examples are shown in Figure <ref>a and c, where the rotation angle of the rotor can be assigned as the natural mechanical coordinate. 
On the other hand, the examples shown in Figure <ref>b and d may also, under certain circumstances, sustain cyclic motion, though the general coordinate could be not so obvious.  
As a possibility for the quantum shuttle of Figure <ref>b, the cyclic motion would involve a cyclic reversal of the bias voltage (ac-driven). This ac-driven case, though intriguing, goes beyond the scope of the present manuscript as we are not considering here time-dependent biases.
Another scenario would be that of Figure <ref>d, where the cyclic motion is in principle attainable by periodically hitting the device.
Note that we are not dealing with the steady-state of sets of interacting nanomotors, such as those described in, e.g., Refs. <cit.>. Instead, here we are interested in the steady-state of the mechanical part of isolated quantum motors and pumps that interact solely with the electrons of a set of reservoirs, and where, at most, Coulomb interactions are only taken into account within the local system.

To discuss the dynamics of cyclic motions in simple terms, we start by projecting Equation (<ref>)
on a closed trajectory defined in the space ofX.
By assuming a circular trajectory, the dynamics can be described by an angleθ, its associated angular velocityθ̇, the moment of inertiaI, and the torquesF,F_ext, andξ_θ.
Using this, we obtain an effective angular Langevin equation equivalent to Equation (<ref>),

    θ̇t = 1/I[ F - F_ext +ξ_θ] .

We assume that after a long waiting time the system arrives to the steady-state regime where the mechanical motion becomes periodic, and it is then characterized by a time periodτsuch thatθ(t+τ)=θ(t)andθ̇(t+τ)=θ̇(t).
Moreover, we will assume that the stochastic force plays a minor role in the above equation, such that it does not affect the mean values of the dynamical variablesθandθ̇, i.e., the mean trajectories with or without the stochastic force approximately coincide.
This occurs, for example, at low temperatures or in mechanical systems with a large moment of inertia <cit.>.
In the following we will just ignoreξ_θfor practical purposes.[
This is the opposite regime of another type of nanomotors: The Brownian motors <cit.>.
]
Under the above assumptions, we integrate both sides of Equation (<ref>) from an initial positionθ_ito a final oneθ_fand obtain:

    I/2[θ̇_f^2-θ̇_i^2]     =    ∫_θ_i^θ_f [ F  - F_ext ] dθ.

The torques in this equation are, in general, intricate functions of bothθandθ̇.
Therefore, the calculation of theθ-dependent angular velocity usually requires the resolution of a transcendental equation (see, e.g., <cit.>).
Alternatively, one can obtainθ(t)from the numerical integration of the equation of motion by standard techniques like, e.g., the Runge-Kutta method.
All this greatly complicates the study of quantum motors and pumps, to the point where it becomes almost impossible to draw any general conclusion.
For this reason, one common simplification consists in taking the terminal velocity as constant during the whole cycle
 <cit.>.
Indeed, this description is exact if the external agent compels the constant velocity condition to be fulfilled in a controllable manner, as is often conceived in quantum pumping protocols. However, this is not the case in general, and typically one expects internal variations forθ̇in one period. We now address this interesting issue in more detail. First, we take the integral in Equation (<ref>) over the whole period. This gives

    W_ext = W_F,   where  W_F= ∫_0^τF θ̇ dt,   and  W_ext = ∫_0^τF_ext θ̇ dt.

The above stationary state condition thus establishes that the work originated from the CIF is always compensated by the external mechanical work in the case this regime can be reached.
Now, let us assume for a moment that the terminal velocity of a nanodevice is constant and positive (we leave the discussion of the effect of the sign ofθ̇for later when treating a concrete example in Section <ref>).
If we now expandFin terms ofθ̇, Equation (<ref>) yields

    𝒲_ext    =    ∑_k(∫_0^2π. ^k Fθ̇^k|_θ̇=0dθ/k!) θ̇^k .

Two important conclusions can be extracted from the above formal solution.
First, there may be conditions where some roots of Equation (<ref>) are complex numbers, meaning that
the assumptionθ̇ = const.is nonsense, as the periodicity condition required for the steady-state regime would not be fulfilled. 
Second, for real solutions, it was shown in <cit.> that the moment of inertiaInot only affects the time that takes to the mechanical system to reach the stationary regime, but also the internal range in the angular velocity, i.e., the differenceΔθ̇= θ̇_max - θ̇_minin one period.
According to Equation (<ref>),θ̇is independent ofI, while the variation ofθ̇scales withI^-1, cf. Equation (<ref>).
Then the ratioΔθ̇/θ̇, which is the relevant quantity in our analysis, should vanish for largeIvalues, justifying the constant velocity assumption for large or massive mechanical systems.

The value ofF_extis supposed to be controllable externally as well as the voltage and temperature biases which, in turn, affect the current-induced torqueF.
Therefore, under the above discussed conditions,θ̇can be thought as a parameter that surely depends on the internal details of the system, but it is also tunable by external “knobs”. Let us analyze a concrete example:
Consider a local system connected to two leads at the same temperature and with a small bias voltageeV=μ_L-μ_R. By considering the current induced torque up to its first nonadiabatic correction, i.e.,F ≈F^(0) - γθ̇,
and assuming thatF_extis independent ofθ̇, the following relation must hold, according to the above discussion,

    θ̇≈Q_I_R/2πγ̅( V- W_ext/Q_I_R) = Q_I_R 
    V_eff/2πγ̅,

whereγ̅is the average electronic friction coefficient along the cycle,Q_I_Ris the pumped charge to the right lead, and we used the Onsager's reciprocal relation betweenFand the charge currentI_Rin the absence of magnetic fields <cit.>.
Alternatively, if we assume that the external torque is of the formF_ext = γ_extθ̇, one finds

    θ̇≈Q_I_R V/2π1/(γ̅+γ̅_ext) = Q_I_R V/2πγ_eff.

Note in the above equations that, at least in the present order, the effect of the external forces can be described as a renormalization of the bias voltageVor the electronic friction coefficientγ.
Numerical simulations in <cit.> show that above equations agree well in general with the steady-state velocities found by integrating the equation of motion.
However, at very small voltages, essential differences may appear.
There is a critical voltage below which the dissipated energy per cycle can not be compensated by the work done by the CIF and thusθ̇= 0. We dubbed this the “nonoperational” regime of the motor.
Moreover, when increasing the bias voltage, there is an intermediate region where a hysteresis cycle appears, and two values of the velocity are possible (θ̇= 0and those given
by the above equations).
Although in Section <ref>, we will takeθ̇as constant when discussing a specific example, the reader should keep in mind that this approximation not always holds, especially at very small voltages orI.



§ ORDER-BY-ORDER ENERGY CONSERVATION


In the previous sections we introduced and discussed the role of the mechanical degrees of freedom, emphasizing certain parameter restrictions which allow for the simplification in their dynamics. In this section we are going to derive, on general grounds, essential relations between the electronic and mechanical degrees of freedom from the point of view of energy conservation. Importantly, we will focus in sistematic expansions beyond the standard linear regime of nonequilibrium sources like, e.g., the mechanical velocity, the bias voltage and the temperature gradient in systems composed by an arbitrary number of reservoirs.

As already pointed out, we are treating the mechanical degrees of freedom classically, such that their effect on the electronic degrees of freedom enters as a parametric dependence in the electronic part of the total Hamiltonian, which now reads

    Ĥ = ∑_r Ĥ_r + Ĥ_s + ∑_r Ĥ_s,r.

The mechanical part of the local system, when treated classically in Equation (<ref>), introduces an explicit time dependence intoĤ_sthat, in turn, it makesd⟨Ĥ|/⟩dt ≠0. According to the above Hamiltonian, the total internal energy of the electronic system,U=⟨Ĥ|$⟩, can be split into energy contributions from the reservoirs, the local system, and the tunnel couplings.
The time-variation of the internal energies associated with the different partitions of the system are:

    U̇_β	= J_β^E + ⟨Ĥ_βt⟩,   where   
    J_β^E = i ⟨[Ĥ,Ĥ_β]|⟩eq:dotU_l

is the mean value of the energy flux entering in the subsystem β= {r , (s,r), s}.
The value of ⟨∂_t Ĥ_β|$⟩ is zero when we evaluate the energy flux in the reservoirs and the tunnel couplings, but equals-Ẇ_Fwhenβis evaluated in the local system.[Although strictly speaking W_F is time-independent when considering a closed trajectory, we use the symbol Ẇ_F to denote the power delivered by the CIF, i.e., Ẇ_F = - ⟨∂_t Ĥ_s|=⟩F·Ẋ.] 
The latter comes from the definition of the CIF given in Equation (<ref>).
Note that the energy fluxes fulfill the condition∑_r (J_r^E+J_s,r^E)+J_s^E = 0. Therefore, the variation of the total internal energy of the electrons yields the following conservation rule

    U̇ = ∑_r U̇_r+∑_r U̇_s,r+U̇_s  = ∑_r ( J_r^E + J_s,r^E ) + J_s^E -Ẇ_F = -Ẇ_F.

Conservation of the total number of particles implies a relation between the particle currents of the reservoirsṄ_rand that of the local systemṄ_sṄ_total = t⟨N̂_total|=⟩∑_r t⟨N̂_r|+⟩t⟨N̂_s|=⟩ 0   ⇒  ∑_r Ṅ_r = -Ṅ_s.

This is so since no particle can be assigned to the coupling region.
Now, let us assume the system is in the steady-state regime and integrate Equation (<ref>) 
over a time periodτof the cyclic motion. Under this condition, the above defined local quantities only depend periodically on time, i.e.,U_s(t+τ) = U_s(t),U_s,r(t+τ) = U_s,r(t), andṄ_s(t+τ) = Ṅ_s(t). Therefore, the following quantities should evaluate to zero, i.e.,

    ∫_0^τU̇_s,r dt = 0 ,   ∫_0^τU̇_s dt = 0,   and  ∫_0^τṄ_s dt = 0,

as these are integrals of a total derivative of some periodic function. This means that no energy (or particles) is accumulated/extracted indefinitely within the finite regions defined by the local system or its coupling to the leads. 
Equation (<ref>) can be used together with Equation (<ref>)
to prove charge current conservation between reservoirs,∑_r ∫_0^τI_r   dt = 0, where the charge current of ther-reservoir is defined asI_r = e Ṅ_r, withe>0being minus the electron's charge.

We will take the following definition for the heat currentJ_rin the reservoirr,[
In <cit.> the authors proposed a different definition for the heat current of the reservoirs, J_r=J_r^E+(J_s,r^E/2)-μ_rṄ_r.
However, the inclusion or not of half the heat current of the coupling region, (J_s,r^E/2), does not make any difference in the present paper as this quantity integrates to zero over a cycle, Equation (<ref>), and does not contribute to the rate of entropy production as we will see in the next section, Equation (<ref>).
]
    J_r = J_r^E - μ_r Ṅ_r.

Replacing this expression in Equation (<ref>), and integrating over a period, results in

    ∑_r (Q_J_r+Q_I_rδ V_r) = -𝒲_F,   where   Q_J_r=∫_0^τ J_r  dt, 
      and   Q_I_r=∫_0^τ I_r  dt,

whereδV_r = (μ_r - μ_0)/eandμ_0is an arbitrary reference's potential.
We, in addition, defined the quantitiesQ_J_randQ_I_rwhich are, respectively, the total heat and charge pumped to reservoirrin a cycle.
Note that in the above equation we have used conservation of the total charge.

Before to continue, we would like to emphasize that the periodic motion imposed by the steady-state regime can be further exploited to reduce the number of mechanical coordinates to an effective description of the CIF. In principle, one can always recognize a generalized coordinateχ, that parameterizes the closed mechanical trajectory. This, for example, can be accounted for by takingχ= Ω  mod(t,τ), such thatχ̇= Ω= 2π/τ. In other words, one can always find a natural scale,Ωin the present case, for all mechanical velocities along the trajectory. We will take this natural scale as the expansion variable for the CIF.
The current-induced work then readsW_F = ∫_0^τF_Ω  Ω  dt, whereF_Ω= F ·∂_χX. Of course, with thisχ-parameterization, the velocity is constant by definition, but at the expense that now the calculation ofF_Ωrequires the knowledge of the mechanical trajectory.
[
For circular trajectories and the conditions discussed in Section <ref>, the coordinate χ would be 
|θ|. For other cases, it could not be that simple to find χ and one should first solve the system's equation of motion.
]

In addition to the mechanical velocityΩ, the CIF and the currents can also be expanded in terms of the remaining nonequilibrium sources, 
corresponding to voltage and temperature deviationsδV = (δV_1, δV_2,...)^TandδT = (δT_1, δT_2,...)^T, respectively, from their equilibrium valuesV_eqandT_eq. For an arbitrary observableR, this general expansion takes the form:

    R = ∑_|α| ≥ 0 R^α = ∑_|α| ≥ 0(Ω,δ V, δT)^α/α!( ∂^α R )_eq,

where we used the following multi-index notation:α=(n_Ω,n_V_1,...,n_T_1,...);α!=n_Ω!n_V_1!...n_T_1!...;(Ω,δV,δT)^α = Ω^n_ΩδV_1^n_V_1...δT_1^n_T_1...;
and

    ∂^α = ∂^(n_Ω+n_V_1+...n_T_1+...)/∂Ω^n_Ω∂ V_1^n_V_1...∂ T_1^n_T_1....

This, in turn, allows us to recognize in the integral quantities of Equation (<ref>) the following expansion coefficients

    W_F^α = ∫_0^τ F_Ω^αΩdt,    Q_J_r^α = ∫_0^τ J_r^αdt,
      and   Q_I_r^α = ∫_0^τ I_r^αdt,

whereF_Ω^α,J_r^α, andI_r^αall have the form given by Equation (<ref>). To find a consistent order-by-order conservation relation from Equation (<ref>), we should first note that the involved terms enter in different 
orders:W_F^αcontains an additionalΩfrom the time-integral as compared with the pumped currents, while the termQ_I_r^αδV_ris one order higher inδV_rthan the other two. Therefore, the following relation must hold for every orderαof the expansion

    Q_J^α·1 + Q_I^α(V)·δV = - W_F^α(Ω),

where1=(1,1,...)^T,Q_J^α= (Q_J_1^α,Q_J_2^α,...),Q_I^α(V) = (Q_I_1^α(V_1),Q_I_2^α(V_2),...),
and we used the shorthandα(Ω) = (n_Ω-1,n_V_1,...,n_T_1,...)andα(V_r) = (n_Ω,n_V_1,...,n_V_r-1,...,n_T_1,...). Equation (<ref>) results important for a systematic study of far-from-equilibrium systems as it helps to rationalize how every order of the expansion of an energy flux is connected with the others. The equation is completely general and valid for any value of the relevant quantitiesΩ,δV, andδT. Importantly, the above conservation rule can be extended to other nonequilibrium sources like, e.g., spin polarization in ferromagnetic leads, such that other type of currents could also be considered in the energy transfer process between electronic and mechanical parts of the system.



§ ENTROPY PRODUCTION AND EFFICIENCY


In 1865, Rudolf Clausius proposed a new state function, the thermodynamic entropyS, that turned out to be crucial to study the limits and the efficiency of different physical processes.
The thermodynamic entropy is defined as the amount of heatδQ_Jwhich is transferred in a reversible thermodynamic process,δS = δQ_rev/T.
Here we are assuming that each reservoirris in its own equilibrium at a constant temperatureT_rand chemical potentialμ_r(we are not going to treat time-dependentT_rorμ_r). As the reservoirs are considered to be macroscopic, the aforementioned equilibrium state is not altered by the coupling to the local system. This assumption allows us to associate the reservoirs' heat flux with the variation in their thermodynamic entropy, i.e.,δQ_J_r = T_r δS_r.
[
The information theory entropy, known as the Shannon or von Neumann entropy, times the Boltzmann constant equals the thermodynamic entropy for equilibrium states.
There is a debate on whether this equality can be extended to nonequilibrium states, see, e.g. <cit.>.
However, for the purpose of this article, only the change of the entropy of the reservoirs is needed, not that of the system.
Besides, we will not need to evaluate the entropy from the density matrix.
Therefore, the thermodynamic definition of entropy suffices in our case.
]
Therefore, from Equation (<ref>) and the definition of the heat current given in Equation (<ref>), we can write

    U̇_r = T_r Ṡ_r+ Ṅ_r μ_r.

Now, let us consider the sum of the internal energy over the set of all reservoirs coupled to the local system,

    ∑_r U̇_r = ∑_ r T_r Ṡ_r + ∑_r Ṅ_r μ_r.

If we takeT_r=δT_r+T_0andμ_r=μ_0+δμ_r,
[The values of μ_0 and T_0 are completely arbitrary, and there is no need to identify them with the chemical potential and temperature of the central region, which can be ill-defined far from
equilibrium.] add and subtract the change of the internal energy of the local systemU̇_sand that of the couplings between the local system and the reservoirs∑_rU̇_s,r, one can rearrange the above equation to the following:

    T_0 ∑_r Ṡ_r = ( U̇_s+∑_r U̇_r+∑_r U̇_s,r) - ∑_r Ṅ_r δμ_r-∑_r Ṡ_r
    δ T_r - (U̇_s+∑_r U̇_s,r) - μ_0 (∑_r Ṅ_r ).

ReplacingṠ_rby=J_r/T_rin the right-hand side of the above equation, using energy
conservation (<ref>), and particle number conservation
(<ref>), allows one to rewrite Equation (<ref>)
as

    T_0Ṡ_res    =     -F_ΩΩ-∑_rI_rδ V_r-∑_rJ_r(δ T_r/T_r)-(U̇_s+∑_rU̇_s,r)+μ_0Ṅ_s,

whereṠ_res = ∑_rṠ_ris the variation of the entropy of the electrons of all reservoirs, and we usedeṄ_r=I_r.
The CIF can be split into “equilibrium” and “nonequilibrium” terms,F^(eq)andF^(ne)respectively, where one can prove thatF^(eq)is always conservative <cit.>.
We are interested in the steady-state situation of our local system. As discussed around Equation (<ref>), the change of the internal energy of the electronic part of the local system and that of the coupling region must be zero after a cycle, as energy cannot be accumulated indefinitely within a finite region.
The same argument is true for the number of particles accumulated in a cycle, which should be zero. At steady state, we therefore recognize the reversible component of the entropy variation as that given by

    Ṡ_res^(rev) =-1/T_0(F^(eq)_ΩΩ+U̇_s+∑_rU̇_s,r-μ_0Ṅ_s).

Obviously, this quantity will not contribute to the total entropy production. Therefore, the rate of entropy productionṠ_res^(irrev)yields

    Ṡ_res^(irrev) = - 1/T_0( F^(ne)_ΩΩ - I·δV-J·δT),

whereδT = (δT_1 / T_1, δT_2/T_2,...)^T, and the currents are defined throughI = (I_1,I_2,...), andJ = (J_1,J_2,...). Integrating Equation (<ref>) over a cycle and taking into account the second law of thermodynamics, one finds

    0 ≥W_F + Q_I ·δV + Q_J ·δT.

The above general formula, also valid far from equilibrium, can be used to set efficiency bounds for energy transfer processes between the electronic and mechanical degrees of freedom.
For example, if we takeδT = 0,Q_I ·δV < 0, andW_F>0, then the system should operate as a nanomotor driven by electric currents, and the following relation holds

    1 ≥ -W_F/Q_I ·δV,

while forQ_I ·δV>0andW_F<0, the system operates as a charge pump, and Equation (<ref>) implies

    1 ≥ -Q_I ·δV/W_F.

Notice that, because of the steady-state condition,W_FequalsW_ext, whereW_extcan be taken as the output or the input energy, depending on the considered type of process.
Therefore, the above formulas describe the efficiencyηof the device's process, defined as the ratio between the output and input energies per cycle. It is also interesting to note that the above equations reflect no more than energy conservation in this particular case.
A different situation occurs forδV=0andδT ≠0, where Equation (<ref>) yields

    1 ≥ -W_F/Q_J ·δT,   and  
    1 ≥ -Q_J ·δT/W_F.

The first equation thus corresponds to a quantum heat engine and the second one to a quantum heat pump, respectively. Now, because of the factorδT, the above formulas differ from what it is expected from energy conservation solely. This is clear in a two-lead system, whereηis limited by the Carnot's efficiency of heat engines and refrigerators, respectively. To illustrate this, let us consider a hot and a cold reservoir, and set the temperature of the cold reservoir as the reference. For the heat engine, this gives

    Q_J ·δT = Q_J_hot(1-T_cold/T_hot) < 0   ⇒  (1-T_cold/T_hot) ≥ 
    -W_F/Q_J_hot,

where the left-hand side of the second equation represents the Carnot limit for heat engines.
Other energy transfer processes mixing voltage and temperature biases can also be analyzed in the context of Equation (<ref>) to set the bounds of their associated efficiencies.



§ PUMP-MOTOR RELATIONS


It is clear from Equation (<ref>) that there is an infinite number of relations that can be used to connect the pumped heat or charge and the work done by the CIF.
In this section we give some physical interpretation to the leading orders in the general expansion, which highlights the utility of Equation (<ref>).

We start from the order-by-order relations by taking|α|=0. Importantly, Equation (<ref>) with|α|=0seems to impose the evaluation of terms with negative coefficients. For example, in the exponent ofQ_Ione would be tempted to evaluateα(V_1)=(0,-1,0,0...), but this term is not defined in the expansions of Equation (<ref>), and then it should be taken as zero. Therefore, 
Equation (<ref>) yields,

    Q_J^0·1 = ∑_r ∫_0^τ( J_r )_eq dt = 0.

This simply reflects the fact that no net heat current occurs at equilibrium.
For|α|=1all relations are summarized in the following three cases:

    ∑_r ∫_0^τ(J_rΩ)_eqΩ dt = 0,   ∑_r ∫_0^τ(J_rV_i)_eqδ V_i  dt = 0,   and  ∑_r ∫_0^τ(J_rT_i)_eqδ T_i  dt = 0,

where we used the fact that equilibrium forces are conservative <cit.>,
i.e,∫_0^τF_Ω^(eq) Ωdt=0,
and that there are no net charge currents in equilibrium, thus∫_0^τ(I_r)_eq dt=0. The above relations mean that, at first order, there is a conservation of pumped heat between reservoirs.
For|α|=2, there are many relations and cases, but we restrict ourselves to only a few of them. Forn_Ω=2,n_V_i=0, andn_T_j=0, Equation (<ref>) gives

    ∑_r_0^τΩ^2/2(∂^2J_r/∂Ω^2)_eqdt     =     -_0^τ(∂ F_Ω/∂Ω)_eqΩ^2dt.

Now, the quantity(∂_ΩF_Ω)_eqis minus the electronic friction coefficient at equilibrium. Therefore, this relation shows 
that the energy dissipated as friction in the motor is delivered as heat to the reservoirs. More precisely, as second order pumped heat. Forn_Ω=0,n_V_i=1,n_V_j=1, andn_T_k=0, Equation (<ref>) yields

    ∑_r_0^τδ V_iδ V_j(∂^2J_r/∂ V_i∂ V_j)_eqdt   =   -2_0^τδ V_iδ V_j(∂ I_i/∂ V_j)_eqdt,

where we used the Onsager's reciprocity relation(∂_V_j I_i)_eq = (∂_V_i I_j)_eq<cit.>.
The quantities(∂_V_j I_i)_eqare the linear conductances in the limit of small bias voltages.
Therefore this relation shows that these leakage currents, defined as those currents which can not be used to perform any useful work, are also dissipated as heat in the reservoirs, a phenomenon known as Joule heating or Joule law <cit.>.
Finally, forn_Ω=1,n_V_i=1, andn_T_k=0,
Equation (<ref>) results in

    ∑_r_0^τΩδ V_i(∂^2J_r/∂Ω∂ V_i)_eqdt+_0^τΩ(∂ I_i/∂Ω)_eqδ V_idt     =     -_0^τδ V_i(∂ F_Ω/∂ V_i)_eqΩdt.

Now, using the Onsager's reciprocity relation(-∂_V_i F_Ω)_eq = (∂_ΩI_i)_eq<cit.>, one finds that the first term in the left-hand side of the above equation vanishes, which is an unexpected conservation relation for this second order pumped heat.
We remark that the utility of Equation (<ref>) relies on the fact that it provides a physical interpretation for the connection between different order contributions that participate in the energy conservation rule.
One can continue analyzing the other relations for|α|=2and beyond, but the number of relations and cases grow very fast with|α|and each relation may have its own physical interpretation.
The study of higher-order terms in|α|may result useful when addressing particular nonlinear effects in the involved energy currents.
However, for the purpose of the present article we believe that the above analysis is enough as to illustrate the approach proposed by Equation (<ref>) regarding multi-index expansions.

As the full dynamics of the complete system typically involves a formidable task, many methods in quantum transport treat this problem through a perturbative expansion inΩ.
This is the case of the real-time diagrammatic theory we use in Section <ref>, where the effective dynamics of the electronic part of the system is described through a perturbative expansion in the characteristic frequency of the driving parameters (which in our case is modeled by the mechanical system), while the voltage and temperature biases are treated exactly.
In this case, the expansion inΩcomes naturally from the theory itself.
In situations like this, it may result more useful to simplify the expansions by restricting ourselves to those nonequilibrium sources whose perturbative treatment is inherent to the used formalism, as in in <cit.>.

Regarding the above discussion, we now use Equation (<ref>) to describe how the different energy contributions are linked in anΩexpansion provided by the theory. The zeroth order terms inΩreads

    Q_J^(0)·1 + Q_I^(0)·δV = 0.

This equation shows that the total amount of heat delivered by the leads comes from the bias voltage maintained between them. Its interpretation
is similar to that of Equation (<ref>). The heat term can be associated with the leakage energy current, i.e., the energy flowing
from source to drain leads without being transferred to the local system. The next order in this expansion yields

    Q_J^(1)·1 + Q_I^(1)·δV = -W_F^(0),

and can be understood as a generalization of the motor-pump relation,Q_I^(1) ·δV ≈- W_F^(0), discussed in <cit.> for arbitrary bias voltages and temperature gradients. Therefore, according to Equation (<ref>), deviations of the mentioned relation at finite voltages are due to the pumped heat induced by the mechanical motion of the local system.
When we evaluate the currents in equilibrium by settingδV = 0andδT=0, Equation (<ref>) impliesW_F^(0)=0, meaning that no external work is done in a cycle, and the pumped energy from the leads can again be considered as a leakage current since no net effect on the mechanical system is performed.
If, on the other hand, some bias is present (either thermal or electric), the energy transfer from the leads to the local system imprints a mechanical motion which, in turn, produces some useful work.
The second order term inΩgives

    Q_J^(2)·1 + Q_I^(2)·δV = -W_F^(1),

and generalizes Equation (<ref>) to finite voltage and temperature biases. The right-hand side of Equation (<ref>) can be interpreted as the dissipated energy of the mecanical system, which is delivered to the electronic reservoirs. However, in the above equationW_F^(1)is not guaranteed to be always negative and, from the point of view of the mechanical system, this can be interpreted as a negative friction coefficient.

Now, we return to the multi-variable expansion of the energy currents to establish which orders should be considered in a consistent calculation of the efficiency of quantum motors and pumps. Assuming that in Equation (<ref>) we take|α|up to some truncation valueα_max, the order-by-order scheme implies that, for example, the efficiency of the electrically driven quantum motor should be given by

    η = - ∑_|α| = 0^α_maxW_F^α(Ω)/∑_|α| = 0^α_maxQ_I^α(V)·δV.

To illustrate this, let us take the case of a local system coupled to left and right reservoirs at voltagesV_LandV_R, respectively.
We assume the leads are at the same temperature and we setδV_L = - δV_R = δV/2as the voltage biases. Depending on which kind of expansion we take, one can obtain different expressions for the efficiencies. On the one hand, when expanding in terms ofδV = V_L-V_RandΩup toα_max = 2, the above general expression yields

    η = - W_F^(0,0)+W_F^(1,0)+W_F^(0,1)/(Q_I^(0,0)+Q_I^(1,0)+Q_I^(0,1)) δ V = 
    -W_F^(1,0)+W_F^(0,1)/(Q_I^(1,0)+Q_I^(0,1)) δ V,

where the superscripts indicate the order in the expansions inΩandδV, respectively, and we definedI = (I_L - I_R)/2.
As we already mentioned, the zeroth order contributionsW_F^(0,0)andQ_I^(0,0)are simply zero as they correspond to the work done by
the conservative part of the CIF and the equilibrium charge current, respectively. On the other hand, when performing an expansion up toα_max = 2but only in terms ofΩ, Equation (<ref>) turns into

    η = - W_F^(0)+W_F^(1)/(Q_I^(0)+Q_I^(1)+Q_I^(2)) δ V,

where nowW_F^(0)andQ_I^(0)are nonzero in general, since they are not necessarily evaluated at equilibrium. Although we here restricted ourselves to the efficiency of a nanomotor driven by electric currents, its extension to other operational modes of the device can be obtained from Equation (<ref>) in a similar way to that of Equation (<ref>). This procedure then allows us to obtain efficiency expressions for arbitrary expansions in the nonequilibrium sources, which could be useful in the evaluation of the device's performance far from equilibrium.



§ QUANTUM MOTORS AND PUMPS IN THE COULOMB BLOCKADE REGIME


In this section we consider the CIFs in the so-called Coulomb blockade regime of transport. In this regime, the strong electrostatic repulsion that takes place inside a small quantum dot (usually taken as the local system) highly impacts the device's transport properties, as for small bias voltages no additional charges can flow through the dot and the current gets completely blocked. The full system dynamics in this strongly interacting regime can not be described by, e.g., the scattering matrix approach, and one needs to move to some other theoretical framework. A suitable methodology is given by the real-time diagrammatic theory <cit.>, which allows for an effective treatment of the quantum dot dynamics by performing a double expansion in both the tunnel coupling between the dot and the leads and the frequency associated with the external driving parameters. Since then, many extensions and application examples appeared in the context of quantum pumps <cit.> and quantum motors <cit.>.

To lowest order in the tunnel coupling, the dot's reduced density matrix obeys the following master equation:

    tp = Wp,

where the vectorp = { p_i(t) }describes the dot's occupation probabilities andWis the evolution kernel matrix accounting for the 
transition rates between the quantum dot states, due to its coupling to the leads. In the context of CIFs, we assume that the time scale of the mechanical motion, characterized byẊ, is large as compared to the typical dwell time of the electrons in the local system. This allows for an expansion of the reduced density matrix asp = ∑_k ≥0 p^(k), withp^(k)of order(Ω/Γ)^k. HereΩandΓdenote the characteristic scales for the velocity of the mechanical degrees of freedom and the tunnel rate of the electronic system, respectively, and we always 
assumeΩ<Γ. The above master equation, in turn, takes the following hierarchical structure <cit.>:

    Wp^(0) = 0,   and  Wp^(k) = d/dtp^(k-1).

In the first equation,p^(0)(X)represents the steady-state solution the electronic system arrives when the mechanical system is “frozen” at the pointX. As the mechanical coordinate moves in time, this order corresponds to the adiabatic electronic response to the mechanical motion. This has not to be confused with the steady-state regime of the mechanical system, which obviously takes a much longer time to be reached. The second equation contains higher-order nonadiabatic correctionsp^(k)(X,Ẋ)due to retardation effects in the electronic response to the mechanical motion.
In all these equations, the matrix elements of the kernel are of zeroth order in the mechanical velocities, i.e.,W = W(X). 
Normalization condition on the dot's density matrix impliese^T p^(k) = δ_k0, withe^Tthe trace over the dot's Hilbert space. This allows for the definition of an
invertible pseudo-kernelW̃, whose matrix elements areW̃_ij = W_ij-W_ii, such that we can write:

    p^(k) = [ W̃^-1d/dt]^k p^(0).

Once we know the different orders of the reduced density matrix, it is possible to compute the expectation value of any observableRafter 
calculation of its corresponding kernelW^R. This implies the same expansion as before, and it reads

    R^(k) = e^TW^Rp^(k).

The observables we are going to address here are the charge and heat currentsI_randJ_rassociated with ther-lead, and defined in Section <ref>.
In lowest-order in tunneling, and under the assumption that coherences are completely decoupled from the occupations, it is possible to obtain simple expressions for 
the current kernels in terms of the number of particlesn_iand energyE_iassociated to the dot state|i⟩<cit.>:

    W^I_r_ij = -e(n_i-n_j) [W_r]_ij,      W^J_r_ij = -[ E_i-E_j-μ_r (n_i-n_j) ] [W_r]_ij.

The CIF was derived in Section <ref> under the assumption that the mechanical part of the system, characterized by 
coordinatesXand associated momentaP, only interacts with local parameters of the quantum dot through their many-body eigenenergies, 
cf. Equation (<ref>). This implies that theν-component of the CIF operator, defined asF̂_ν= -∂Ĥ_s/∂X_ν, is local in the 
quantum dot basis. For this local operator, then, we can simply define a diagonal kernel of the form

    W^F_ν_ij = -E_iX_νδ_ij,

such that Equation (<ref>) gives thek-order in theΩexpansion for any observable.
Importantly, when adding up all contributions from the leads, the above kernel definitions, together with the sum rule∑_i W_ij = 0, leads to the following conservation rule for all orders inΩ<cit.>:

    ∑_r J_r^(k) =  - d/dt⟨Ĥ_s|^⟩(k-1) - F^(k-1)·Ẋ - ∑_r μ_r/e I_r^(k).

This equation, equivalent to (<ref>), thus expresses the first principle of thermodynamics, which in this case relates the total heat flowing from/into the leads with the variation of the internal energy of the dot and the power contributions due to both mechanical and electrochemical external sources.
By considering a system coupled to two reservoirsLandR, with periodic motion characterized by a time periodτ, and taking the time integral of the above equation, we recover the frequency expansion of Equation (<ref>):

    Q_J^(k)+Q_I^(k)δ V = -W_F^(k-1),

where the bias voltageδVis defined throughμ_L = - μ_R = e δV/2,Q_J= Q_J_L+Q_J_R, andQ_I= (Q_I_L-Q_I_R)/2.
The sign convention employed here implies, for example, that if the left hand side of Equation (<ref>) is 
positive, then there is some amount of energy entering into the leads in one cycle, and work is being extracted from the local system.



 §.§ Example: Double quantum dot coupled to a rotor 


In this section we illustrate the discussions of the previous sections in a concrete example based on a double quantum dot (DQD) system locally coupled to a mechanical rotor (see Figure <ref>c). We 
assume a capacitive coupling between the dots and the fixed charges in the rotor such that no charge flows between the two subsystems. The DQD system
is described as in Refs. <cit.>, 

    Ĥ_s = ∑_ℓ = L, Rϵ_ℓn̂_ℓ + U n̂_L n̂_R + U'/2∑_ℓ = L,Rn̂_ℓ (n̂_ℓ-1)
    - t_c/2( d̂_Lσ^†d̂_Rσ + h.c.),

whereϵ_ℓis the onsite energy andn̂_ℓ= ∑_σd̂_ℓσ^†d̂_ℓσthe number operator 
in theℓ-dot,UandU'the inter and intradot charging energies, respectively, andt_cthe interdot hopping amplitude. We will work in the 
strong coupling regimet_c ≫Γ, such that nondiagonal elements (coherences) in the reduced density matrix are decoupled from the diagonal 
ones (occupations), and can be disregarded in first order in tunneling. To simplify the analysis (by reducing the number of states in the two-charge 
block), we work in the limitU'→∞, such that double occupation in a single dot is energetically forbidden.

Diagonalization of the above DQD Hamiltonian yields the bonding (b) and antibonding (a) basis for the single-electron charge block, and the reduced density matrix 
readsp = (p_0, p_b↑, 
p_b↓, p_a↑, p_a↓, p_↑↑, p_↑↓, p_↓↑, 
p_↓↓)^T. These elements thus denote the probability for the DQD to be either empty (p_0), occupied by one electron in theℓ= b(ora) orbital with spinσ(p_ℓσ), or by two electrons (p_σσ'), one of them in the left dot 
and with spinσand the other electron in the right dot and with spinσ'. The many-body eigenenergies are thereforeE_0 = 0for the 
empty DQD,

    E_b/a = ϵ_L+ϵ_R/2∓√((ϵ_L-ϵ_R/2)^2+(t_c/2)^2)

for single-occupation in the bonding or the antibonding orbital, andE_2 = ϵ_L+ϵ_R + Ufor the doubly occupied DQD, respectively.

To account for the coupling between electronic and mechanical degrees of freedom we take as mechanical coordinate the angleθdescribing the orientation of the rotor axis. We assume the following dependence through the onsite energies:

    ϵ_L(θ) = ϵ̅_L + λ x_L = ϵ̅_L + δ_ϵcos(θ),   and  ϵ_R(θ) = ϵ̅_R + λ x_R = ϵ̅_R + δ_ϵsin(θ),

which defines a circular trajectory in energy space of radiusδ_ϵ= λr_0(r_0measures the actual radius of the rotor) around 
the working point(ϵ̅_L,ϵ̅_R). For simplicity in what follows, we will focus in the tangential component of the force only, 
by assuming that the radial component is always compensated by internal forces in the rotor (e.g., fixed charges along the rotor's axis). By 
applying this form in Equation (<ref>), we obtain the equation of motion (<ref>) for the angular velocityθ̇, where:I=m_eff r_0^2is the rotor's moment of inertia,F = ∑_i (-∂_θE_i) p_i = r_0 F·θ̂is the current-induced torque,F_extis the torque produced by the external force, which is assumed to be constant along the whole trajectory, andξ_θaccounts for the force fluctuations. The stationary regime is reached once the rotor performs a periodic motion (characterized by a time periodτ) with no overall acceleration, i.e., when Equation (<ref>) is fulfilled. We show in Figure <ref> some examples of the evolution of the rotor's angular velocity for different initial conditions. This is plotted as a function of the number of cycles performed by the rotor, instead of time, to unify scaling along the horizontal axis. Here we take a sufficiently large moment of inertia such that the variation Δθ̇ over one period is small as compared to the averaged value of θ̇ in one cycle. This, in turn, implies a large number of cycles for the rotor to reach the steady-state regime.

As discussed in Section <ref>, when we takeIsufficiently largeΔθ̇ →0, and it is possible to approximateθ̇as constant during the whole cycle. This allows us to take out this quantity from the integrals defining the different orders of the energy currents, as we did in Equation (<ref>). The problem is, however, that the sign inθ̇is not 
a priori known and, with it, the direction of the 
trajectory for the line integrals definingW_F^(k). To become independent of this issue, we denote bysthe sign ofθ̇, and notice thatF^(k) = f_θ^(k) θ̇^k, wheref_θ^(k) = ∂_θ̇^k F |_θ̇=0/k!only depends parametrically onθ, so we have

    W_F^(k) = ∮_C(s)F^(k)·dX
    = ∫_0^2π sF^(k)dθ
    = s [ ∫_0^2π f_θ^(k)dθ] θ̇^k
    = s C_F^(k)θ̇^k,

where, importantly, the coefficientsC_F^(k)are independent ons, so theθ̇^kand the extra signsaccompanyingC_F^(k)in the right hand side of the equation give the correct sign forW_F^(k). As discussed in Section <ref>, the steady-state condition in Equation (<ref>) can be viewed as the equation for the final velocity that the rotor acquires as a function of the external force, which in the present case turns into:

    0 = ∑_k C_F^(k)θ̇^k - C_ext,

where we definedC_ext = 2πF_ext.
The stability of the final solution is inherited from Equation (<ref>) 
and, in our case withθ̇constant, it is given by

    ∑_k C_F^(k) k θ̇^k-1 < 0.

Notice the specific case whereC_ext = C_F^(0), i.e., the external work equals the bias work coming 
from the first order currents, cf. Equation (<ref>) fork=1. In this case there is always a trivial solution, given byθ̇ = 0. Of course, this solution is useless as the system becomes frozen, and there is no energy transfer between the leads 
and the mechanical system. However, this situation marks one of the transition points between the operation modes of the device. If we now 
consider an infinitesimal difference between these two coefficients, we can faithfully treat Equation (<ref>) up to linear orden inθ̇, such that the solution is

    θ̇ = C_ext-C_F^(0)/C_F^(1),

together with the conditionC_F^(1)<0, in agreement with Equation (<ref>). This implies that the integral of the friction coefficient 
over a cycle needs to be positive, otherwise the rotor can not reach the stationary regime. As the sign inC_F^(1)is fixed to this order, the sign inθ̇only depends on the relation betweenC_extandC_F^(0). This, together with the stationary conditionW_ext = W_F, establishes the operation mode of the 
device, in the sense that with this information we can deduce the sign ofW_ext =  s C_ext, and hence the direction 
of the energy flow. If, for example,W_ext < 0, the external energy is delivered into the DQD, which in turn it goes as energy 
current to the leads, so the device acts as an energy pump. On the other hand, ifW_ext > 0the energy current coming from the leads goes into the DQD and it is transformed into mechanical work, so the device operates as a motor. These two scenarios are shown in Figure <ref>, where we take C_ext>0, such that the positive sign in the final velocity implies that the device acts as a motor (see panel a), while a negative sign implies that the device operates as a pump (see panel b).
Of course, when regarding the device as an energy pump, the above criterion only establishes those regions inC_extwhere we can expect some kind of pumping mechanism. Depending 
then on the specific type of pumping we have in mind, the ranges inC_extwill be subject to additional conditions.
For example, if we would like to have this device acting as a quantum charge pump, then we should check those regimes inC_extwhere the electrons flow against the bias voltage. This discussion will be reserved to the next section, and for now we will only define the operation mode from the direction of the energy flow.

Let us now consider second-order effects due toC_F^(2)in Equation (<ref>), where we obtain:

    θ̇_± = -C_F^(1)/2C_F^(2)±√(( C_F^(1)/2C_F^(2))^2+C_ext-C_F^(0)/C_F^(2)),

together with the conditions:

    ( C_F^(1)/2C_F^(2))^2 > C_F^(0)-C_ext/C_F^(2),     ±C_F^(2) < 0.

The first inequality ensures a positive argument in the square root of Equation (<ref>), soθ̇is a real number, while the second 
inequality comes from the stability condition given by Equation (<ref>), and tells us which branch one should choose in Equation (<ref>) to get the stable solution: IfC_F^(2)>0thenθ̇_-, and ifC_F^(2)<0thenθ̇_+. Importantly, far from 
equilibrium (e.g.,δV ≫k_B T) it is possible to arrive to the odd situation whereC_F^(1)>0, and the “dissipated” energyW_F^(1) = s C_F^(1) θ̇ = C_F^(1) |θ̇|is positive. This, however, does not prevent the 
rotor to reach the stationary regime (to this order inθ̇), since the lower-order terms compensate this energy gain. Therefore, one can end up in a situation where the second-order energy current comes from the leads and it is delivered into the mechanical system (thus favouring the motor regime), contrary to the standard situation where the “dissipated” energy flows to the leads. Regarding Equation (<ref>), the sign ofθ̇now depends on the relation between the first term and the square root, together with the sign ofC_F^(2), so this analysis is not that simple as in the linear case. However, once we have theC_F-coefficients we can always determine if the rotor is able to reach the stationary regime and infer whetherW_extis positive or negative and, with it, the operation mode of the device.

For higher orders inθ̇, the above analysis for the operation mode of the device is the same, but more ingredients 
may come into play due to the (order-dependent) specific solutions for the stationary value ofθ̇. Interestingly, by including 
higher-order terms in Equation (<ref>), it could happen that for a fixed choice of the parameters (F_ext,δV,δT, etc.) the system presents more than one stable solution, and even with different signs inW_ext, so the initial condition onθ̇decides the operation mode of the device once the rotor reaches the stationary regime. In Figure <ref>c we show this nolinear effect in the dynamics of the rotor's velocity. This was done by taking F up to third order in the θ̇ expansion, where we find two stable solutions (θ̇∼ -7.4 × 10^-3  k_BT and 6.7 × 10^-3  k_BT ) and one unstable solution in between (θ̇ = 2 × 10^-3  k_BT). Notice that, in this case, if the rotor starts with a positive initial condition below the unstable solution (see gray arrow), then it can not reach the negative stable solution, as once θ̇= 0 the rotor gets trapped in a local minimum.In Figure <ref> we show how the C_F-coefficients contribute to the current-induced work as a function of the bias voltage for different values of the orbit radius δ_ϵ. According to Equation (<ref>), these coefficients can only be compared upon multiplication with the k-power of some frequency of reference Ω_x. We assess this frequency of reference by performing the following analysis: From Equation (<ref>) we can estimate the maximum allowed value for θ̇ compatible with the frequency expansion. Since W̃^-1∝Γ^-1, with Γ = Γ_L + Γ_R, and

    t = θ̇θ = θ̇∑_i ϵ_iθϵ_i,

we obtain

    p^(k) = [ W̃^-1θ̇∑_i ϵ_iθϵ_i]^k p^(0)∝[ Ω/Γδ_ϵ/k_BT]^k,

where we use the fact that ∂_θϵ_i ∝δ_ϵ and the energy derivative of both p^(0) and W̃^-1 are proportional to 1/k_BT. As the consistence of the frequency expansion relies on the convergence of the occupations, this yields the adiabaticity condition <cit.>
    Ω/Γ≪k_BT/δ_ϵ,

from which we can estimate the maximum allowed frecuency Ω_max as

    Ω_max = k_BT/δ_ϵΓ.

Obviously, this extreme value sets the point where the expansion could diverge, so to illustrate the different order contributions to the current-induced work we take, in Figure <ref>, an intermediate referece frequency Ω_x = Ω_max/2.
As we can see, for a long range of the bias voltage (V ≲ 10   k_BT), both the zeroth and first order terms contribute, while higher-order terms are almost negligible. The zeroth order contributions (black) show a linear dependence whose slope increases with the orbit size, up to some saturation value, related with the quantization of the pumped charge. The first order terms (red) remain almost constant and negative in this bias regime. This implies that the linear order treatment given in Equation (<ref>) is enough, as long as C_ext∼C_F^(0). For V ≳ 10   k_BT, we can observe non-equilibrium effects like “inverse dissipation”, i.e. a positive contribution from the first order term (red) in Figure <ref>a. Additionally, the higher-order terms (blue and green) can become larger than the first two and they need to be included in the calculation of the final velocity through Equation (<ref>) or in the current-induced force appearing in Equation (<ref>). This, in turn, could lead to several stable solutions whose validity should be determined through a systematic convergence analysis that is beyond the scope of the present work. For the particular choice of parameters used in Figure <ref>c we checked that the next-order coefficient (C_F^(4)) has a negligible impact on the third order solutions.To clarify the role of C_F-coefficients in the operation regime of the device, let us take for example δ_ϵ = 100   k_B T (see Figure <ref>d) and δ V ∼ 2   k_BT. As the sign in C_F^(1) is negative, the linear order solution given by Equation (<ref>) is stable. If we start with 
C_ext = 0, then θ̇ is positive, so for 0 < C_ext < C_F^(0) the device acts as a motor (see Figure <ref>a). When C_ext > C_F^(0) > 0, θ̇ becomes negative and hence the device operates as a pump (see Figure <ref>b). Additionally, when C_ext < 0 the angular velocity is still positive, but as we changed the sign in C_ext we have that W_ext < 0, so the device operates as a pump. This is the other transition point between the two operation modes of the device, i.e., the sign in W_ext = s C_ext changes with s and the sign of C_ext. It is also interesting to notice how the operation regimes change with the sign of the bias voltage. For the force coefficients, we can see in the figure that they have definite parity with respect to the bias voltage

    C_F^(k)(-δ V) = (-1)^k+1C_F^(k)(δ V),

since Ω_x remains constant for a fixed orbit radius. For the chosen parameters ϵ̅_L = ϵ̅_R, Γ_L = Γ_R, and δ T = 0, the transformation δ V → - δ V can be regarded as the inversion operation (L,R) → (R,L), which changes the sign of θ̇, i.e., θ̇(-δ V) = - θ̇(δ V). In this sense, if there is some finite temperature gradient, this operation should also involve the sign inversion of δ T. To infer how this bias inversion affects the final velocity of the device, we can replace Equation (<ref>) in Equation (<ref>). We can therefore recognize the transformation θ̇(-δ V) = -θ̇(δ V) if we also change the sign of the external force, such that C_ext(-δ V) = - C_ext(δ V). The even/odd parity in the k-order coefficient thus implies that the current-induced work is invariant under such a transformation, cf. Equation (<ref>), and the ranges for the operation regimes of the device remain the same if we invert the sign of the external force. Of course, this analysis is no longer valid in more general situations where ϵ̅_L ≠ϵ̅_R or Γ_L ≠Γ_R, such that the change in the sign of δ V (or δ T) can not be related with the left-right inversion operation.

As an additional test for Equations (<ref>) and (<ref>), we define equivalent coefficients for the amount of transported charge and heat in a cycle

    C_I^(k) = ∫_0^2π.^k I_θθ̇^k|_θ̇=0dθ/k!,     C_J^(k) = ∫_0^2π.^k J_θθ̇^k|_θ̇=0dθ/k!,

such thatQ_I^(k) = s C_I^(k) θ̇^k-1, andQ_J^(k) = s C_J^(k) θ̇^k-1. These contributions can be 
evaluated independently from theC_F-coefficients and in Figure <ref> these are shown in circles, which gives numerical agreement for the 
energy conservation principle

    C_J^(k)+C_I^(k)δ V = -C_F^(k-1),

in the considered orders ofθ̇.
For the considered example in Figure <ref>, the current coefficientsC_IandC_Jalso show 
(independently) a definite parity with respect to the bias voltage. In fact, as Equation (<ref>) suggests, the heat current coefficientsC_J^(k)present the same parity asC_F^(k-1), while the charge current coefficientsC_I^(k)need to have the opposite 
parity since they are multiplied byδV.



  §.§.§ Motor-pump efficiencies


As we stated above, the sign of the external force, together with the rotor's stationary condition and the first law of thermodynamics, 
determine the direction of the energy flow and, with it, the operation mode of the device. Obviously, as no other power sources are involved in 
this example, the efficiency of this energy conversion, defined asη=(output power)/(input power), is always equal to 1.
This, however, only establishes those regions in the parameter space where we can expect the device operating as a quantum motor or a quantum energy pump. In this section we discuss the particular conditions that appear when the device operates through a specific type of current. In this sense, the motor regime corresponds to the situation in which a transport current (e.g., charge, heat, spin, etc.) flowing through the leads in response to a bias (voltage, temperature, spin polarization, etc.) delivers some amount of energy into the local system which can be used as mechanical work. The pump regime, on the other hand, corresponds to the inverse operation in which the external work is exploited to produce a current flowing against the imposed bias. This topic was also discussed in <cit.> for charge and heat currents in a DQD device, where limitations to the efficiency of the considered processes were attributed to the different orders appearing in the frequency expansion of the currents. We here provide a similar analysis in terms of our explicit model for the mechanical system. The inclusion of the external force in the description of the model, as we shall see next, appears as the key ingredient in bridging the motor and pump regimes for a given choice of the bias.

For the device acting as a motor, the output power should be given byW_ext/τ, under the conditionW_ext > 0, but we still need to specify the input power. If we consider that the mechanical rotor is driven by the electric current, i.e., due to some applied bias voltage and no thermal gradient applied, the input power is given by-Q_I·δV/τ, and hence the efficiency of this type of motor is

    η = -W_ext/Q_I ·δV = 1 + Q_J/Q_I δ V.

Equation (<ref>) establishes that the maximum efficiency for this device isη≤1, and the above equation tells us that 
the heat current produced by the bias voltage reduces the motor's performance. As in this work we calculate such currents through an expansion
in the angular velocity, the efficiency is also limited by this expansion.
If in the calculation ofθ̇we consider Equation (<ref>) up to first order inC_F^(k), then, as discussed in Section <ref>, order-by-order energy conservation demands that the currents are to be considered up to second order, and in terms of the current coefficients this takes the form:

    η = -C_ext/δ V/C_I^(0)θ̇^-1+ C_I^(1) + C_I^(2)θ̇.

In the limits of the motor's operation regime, given byC_ext = 0andC_ext = C_F^(0), it is easy to 
see that the efficiency goes to zero, since forC_ext = 0the numerator in the above expression is zero, and forC_ext = C_F^(0)the rotor's velocity goes to zero, so the denominator grows to infinity due to the contributionC_I^(0)/θ̇from the leakage current. The same happens if we consistently include higher-order terms in this expression. For example,
if we use Equation (<ref>) for the rotor's angular velocity, then we should addC_I^(3) θ̇^2in the denominator.



Away from this region, we enter in the “pumping domain” characterized by a charge current which opposes to the “natural” direction dictated by 
the bias voltage, andW_ext<0. In this sense, the input power and the output power are inverted with respect to the motor region and,
consequently, the efficiency of this “battery charger” device is given by

    η = -C_I^(0)θ̇^-1+ C_I^(1) + C_I^(2)θ̇/C_ext/ δ V.

In this regime there is, however, an additional condition to be fulfilled, which isQ_I δV > 0. Regarding the different orders inQ_I, it usually  happens that close to the transition pointC_ext = C_F^(0)the charge current still flows in the bias direction, since it is dominated by the leakage current. In this region we say that the pumping mechanism is “frustrated” as the energy delivered by the rotor is not 
enough as to reverse the direction of the charge current. Going away from this region the angular velocity acquires some finite value, reducing the 
zeroth-order contributionC_I^(0)/θ̇to the point where it is equal to the higher-order contributionsC_I^(1) + C_I^(2)θ̇, thus marking the activation point of the charge pump. In Figure <ref>a we show the efficiency of the device as a function of the external force for different orbit sizes and fixed biasδV = 2 k_BT. In all cases, the device starts fromC_ext = 0as a motor and its efficiency reaches a maximum which increases with the orbit size. Soon after this point the motor's efficiency decreases to zero due to the leakage current effect, which becomes dominant atC_ext = C_F^(0).
From this point we can observe the gapped region for the frustrated pump, which is more pronounced for small orbits, since the first-order pumped chargeQ_I^(1)is smaller than its quantized limit and thereby it takes a larger value ofC_extto compensate the amount of pumped leakage currentQ_I^(0)in a cycle.

Up to this point, we have discussed only the effect of the lowest-order terms of the expansion in θ̇, given by Equations (<ref>), (<ref>) and (<ref>). For the smallest orbit, in addition, we show in dotted red line the next-order efficiency, obtained from Equation (<ref>) and adding the Q_I^(3) term to Equations  (<ref>) and (<ref>). We can see that for the motor regime there are no significant changes, but for the pump regime important differences appear.
Firstly, in the C_ext>0 region, there is a cut-off for the external force in which the pumping mechanism is again frustrated, i.e., the charge current again points in the bias direction. As can be seen in Figures <ref>a and c, this decreasing in the efficiency is not  attributed to the extra heat dissipated to the reservoirs, as one may expect. Note in the figure that the extra contribution to the pumped heat, Q_J^(3), is negligible as compared to the lowest order terms. What happens here is that the third-order contribution to the pumped charge Q_I^(3) rapidly becomes dominant in the charge pump region, causing a sudden drop in the efficiency and, with it, the appearance of a second frustrated-pump region. Secondly, another higher-order effect appears in the C_ext<0 region. There, the efficiency is nonzero for C_ext/C_F^(0) < -8 (see inset in Figure <ref>a), meaning that the pumping mechanism can be activated even when the external force points in the same direction as that of the current-induced force. [Given the convention used for F_ext in Equation (<ref>) and the chosen parameters, in this region the sign of θ̇ remains the same as that when F_ext = 0. There, the zeroth and first order contributions to the charge current flow in the same direction.] 
Again, the third-order term is the one that reverses the direction of the total charge current, see Figure <ref>b. It is important to mention that the purpose of the present discussion is only to highlight deviations from the linear solution of Equation (<ref>), not to analyse the convergence of the total pumped current. For the larger values of δ_ϵ used in Figure <ref>a we do not show the next order corrections, as they are negligible in the shown range of C_ext.


An analysis similar to the above one can be carried out for the device driven by a temperature gradientδT, defined throughT_L = T+δT/2andT_R = T-δT/2and no bias voltage applied, such that forδT > 0we haveT_hot = T_LandT_cold = T_R. WhenW_ext>0we have a motor device driven by a heat current in response to a thermal gradient (heat engine), then the input power should be 
given by-Q_J_hot/ τ, and Equation (<ref>) implies

    η = -W_ext/Q_J_hot = 1 + Q_J_cold/Q_J_hot≤ 1-T_cold/T_hot = η_carnot.

As compared with Equation (<ref>), we can see that the efficiency of the quantum heat engine, now given by

    η = -C_ext/C_J_hot^(0)θ̇^-1 + C_J_hot^(1) + C_J_hot^(2)θ̇,

is defined in the same range forC_extas in the electric motor. Regarding Figure <ref>b, the engine's normalized efficiencyη/ η_carnotlooks similar to that of the electric motor. Perhaps the only difference here is that for the smallest orbitδ_ϵ= 10 k_B Tthe efficiency maximum is very low, such that it can not be appreciated on the employed scale of the plot. 



Now we move to the heat pump region where the device acts as a refrigerator, as we demand that the heat current flows against the direction dictated byδT. Therefore, in addition to theW_ext < 0condition, the overall amount of pumped heat in the cold reservoir should be 
negative, i.e.,Q_J_cold < 0. The efficiency of the refrigerator, or coefficient of performance (COP), is then given by

    COP = Q_J_cold/W_ext = Q_J_hot/Q_J-1 ≤T_cold/T_hot-T_cold = COP_carnot,

where we can consistently expandQ_J_coldin terms ofθ̇. In Figure <ref>b we show the lowest order contribution from 
Equation (<ref>) as the next-order calculation does not change significantly the efficiencies in the considered regimes of the parameters. Again, we can observe in Figure <ref>b a gap region where the device is frustrated since the work delivered by the rotor is not enough as to reverse the direction of the heat current.
One of the differences with the electric counterpart is that, for the refrigerator, the normalized COP develops a maximum that is always smaller than that of the quantum heat engine, while the obtained efficiency maxima (motor and pump) for a fixed orbit in Figure <ref>a are very similar. Additionally, for the chosen valueδT = 0.5   Tand small orbit radius, the device can only work as a heat engine, and the refrigerator can not be activated even if the external force is large, as happens forδ_ϵ= 20   k_BT(solid blue line). The reason for this relies on the competition between the different orders in the pumped heatQ_J_cold: The reduction of the leakage pumped current, Q_J_cold^(0), by increasing the magnitude ofθ̇is accompanied by an increase of the second-order contribution,Q_J_cold^(2), such that the first order term,Q_J_cold^(1), may not be able to compensate these two and the requirementQ_J_cold < 0can not be fulfilled, see Figure <ref>.



§ SUMMARY


Throughout this manuscript, we revisited some fundamental aspects related with the physics of quantum motors and pumps. Previous results based on the steady-state properties and energy conservation law were extended to deal with arbitrary nonequilibrium conditions in a systematic way. By considering the dynamics of the mechanical degrees of freedom through a Langevin equation, we were able to treat the motor/pump protocols on the same footing. This allowed us to describe the related energy transfer processes through a single parameter: the external force done on the local system.

In the steady-state regime, we treated in general terms the validity of the constant velocity assumption, Section <ref>. For arbitrary orders of the nonadiabatic expansion in the CIFs, this was linked to the separation between the electronic and mechanical dynamic scales through a large moment of inertia.
We then performed a general expansion (in terms of nonequilibrium sources) of the energy fluxes that take part in the quantum transport problem. This enabled us to derive an order-by-order scheme for the energy conservation law, Equation (<ref>). This Equation may be of help in recognizing the physical processes that enter at each order in the expansion, thereby providing a useful tool for the analysis of nonlinear effects. To illustrate this, we discussed the leading orders of the global expansion, and showed how different types of expansions of the energy fluxes change the expressions for the efficiency of quantum motors and pumps.

In Section <ref>, we introduced a specific example of a quantum motor/pump based on a double quantum dot. There, we discussed in more depth how higher-order terms of the CIFs affect the stationary state conditions. We found that multiple solutions for the device's terminal velocity could in principle be available for a fixed choice of parameters (voltage and temperature biases, and external force). In such a case, the stability of such solutions imposes an additional constraint on the force coefficients, and the final steady state strongly depends on initial conditions. Interestingly, it is possible to obtain more than one stable solution, each of them belonging to a different operation mode of the device.
The treated example is also appealing as it is possible to study the transition between different operational modes by continuously moving the external force.
This corresponds to the point in which the steady-state velocity changes its sign and, with it, the direction of the energy flow.
When considering a specific type of pumped currents (charge or heat), there is an intermediate region where the pumping mechanism is “frustrated”.
In this situation, the energy delivered by the external force is not enough as to reverse the natural direction of the charge or heat currents.
We found other interesting features of the studied example such as negative friction coefficients at finite voltages or a definite parity of the expansion coefficients with respect to the bias voltage and the temperature gradient, which is a manifestation of the inversion symmetry in the total energy flux.
We also used this example to numerically confirm the order-by-order energy conservation law up to third order in the final velocity.
Finally, for heat currents, we found parameter conditions under which the device can never work as a “refrigerator”, even for large values of the external force.
We explained this behavior in terms of the competition between the different orders that participate in the pumped heat of the cold reservoir, highlighting the importance of the order-by-order conservation laws.


Acknowlegdments.– This work was supported by Consejo Nacional de Investigaciones Científicas y 
Técnicas (CONICET), Secretaría de Ciencia y Tecnología – Universidad Nacional de Córdoba (SECYT–UNC). All authors are members of CONICET.

apsrev4-1_title