**6. Pump-Motor Relations**

It is clear from Equation (25) 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 (25).

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

$$\mathbf{Q}\_{l}^{\mathbf{0}} \cdot \mathbf{1} = \sum\_{r} \int\_{0}^{\tau} \left( \mathbf{J}\_{r} \right)\_{\mathbf{eq}} \, \mathrm{d}t = 0. \tag{37}$$

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

$$\sum\_{r} \int\_{0}^{\tau} \left(\frac{\partial f\_{r}}{\partial \Omega}\right)\_{\text{eq}} \Omega \,\text{d}t = 0, \quad \sum\_{r} \int\_{0}^{\tau} \left(\frac{\partial f\_{r}}{\partial V\_{i}}\right)\_{\text{eq}} \delta V\_{i} \,\text{d}t = 0, \quad \text{and} \quad \sum\_{r} \int\_{0}^{\tau} \left(\frac{\partial f\_{r}}{\partial T\_{i}}\right)\_{\text{eq}} \delta T\_{i} \,\text{d}t = 0, \tag{38}$$

where we used the fact that equilibrium forces are conservative [13,22,42,50,63], i.e., *τ* 0 *F*(eq) Ω Ωd*t* = 0, and that there are no net charge currents in equilibrium, thus *τ* 0 (*Ir*)eqd*<sup>t</sup>* = 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. For *n*Ω = 2, *nVi* = 0, and *nTj* = 0, Equation (25) gives:

$$\sum\_{r} \int\_{0}^{\tau} \frac{\Omega^{2}}{2} \left(\frac{\partial^{2} I\_{r}}{\partial \Omega^{2}}\right)\_{\text{eq}} \text{d}t \quad = \quad - \int\_{0}^{\tau} \left(\frac{\partial F\_{\Omega}}{\partial \Omega}\right)\_{\text{eq}} \Omega^{2} \text{d}t. \tag{39}$$

Now, the quantity (*∂*Ω*F*Ω)eq is 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 a second-order pumped heat. For *n*Ω = 0, *nVi* = 1, *nVj* = 1, and *nTk*= 0, Equation (25) yields:

$$\sum\_{r} \int\_{0}^{\overline{\tau}} \delta V\_{\overline{i}} \delta V\_{\overline{j}} \left( \frac{\partial^{2} I\_{r}}{\partial V\_{\overline{i}} \partial V\_{\overline{j}}} \right)\_{\text{eq}} \text{d}t \quad = \quad -2 \int\_{0}^{\overline{\tau}} \delta V\_{\overline{i}} \delta V\_{\overline{j}} \left( \frac{\partial I\_{\overline{i}}}{\partial V\_{\overline{j}}} \right)\_{\text{eq}} \text{d}t,\tag{40}$$

where we used Onsager's reciprocity relation (*∂Vj Ii*)eq = (*∂Vi Ij*)eq [42,49]. The quantities (*∂Vj Ii*)eq are the linear conductances in the limit of small bias voltages. Therefore, this relation shows that these leakage currents, defined as those currents that cannot be used to perform any useful work, are also dissipated as heat in the reservoirs, a phenomenon known as Joule heating or the Joule law [27,43,50,83]. Finally, for *n*Ω = 1, *nVi* = 1, and *nTk* = 0, Equation (25) results in:

$$\sum\_{r} \int\_{0}^{\overline{\tau}} \Omega \delta V\_{i} \left(\frac{\partial^{2} l\_{r}}{\partial \Omega \partial V\_{i}}\right)\_{\text{eq}} \text{d}t + \int\_{0}^{\overline{\tau}} \Omega \left(\frac{\partial I\_{i}}{\partial \Omega}\right)\_{\text{eq}} \delta V\_{i} \text{d}t = -\int\_{0}^{\overline{\tau}} \delta V\_{i} \left(\frac{\partial F\_{\Omega}}{\partial V\_{i}}\right)\_{\text{eq}} \Omega \text{d}t.\tag{41}$$

Now, using Onsager's reciprocity relation (−*∂Vi <sup>F</sup>*Ω)eq = (*∂*Ω *Ii*)eq [22,39,42,49], 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 (25) 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 |*α*| = 2 and beyond, but the number of relations and cases grows very fast with |*α*|, and each relation may have its own physical interpretation. The study of higher-order terms in |*α*| may result in being 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 to illustrate the approach proposed by Equation (25) 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 7, 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 in being more useful to simplify the expansions by restricting ourselves to those nonequilibrium sources whose perturbative treatment is inherent to the formalism used, as in [42,87].

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

$$\mathbf{Q}\_{l}^{(0)} \cdot \mathbf{1} + \mathbf{Q}\_{l}^{(0)} \cdot \delta V = 0. \tag{42}$$

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 (40). 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:

$$
\mathbf{Q}\_I^{(1)} \cdot \mathbf{1} + \mathbf{Q}\_I^{(1)} \cdot \delta \mathbf{V} = -\mathcal{W}\_{\mathbf{F}}^{(0)} \, , \tag{43}
$$

and can be understood as a generalization of the motor–pump relation, *Q*(1) *I* · *δV* ≈ −W(0) *F* , discussed in [13] for arbitrary bias voltages and temperature gradients. Therefore, according to Equation (43), 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* = 0 and *δT* = 0, Equation (38) implies W(0) *F* = 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:

$$\mathbf{Q}\_{I}^{(2)} \cdot \mathbf{1} + \mathbf{Q}\_{I}^{(2)} \cdot \delta \mathbf{V} = -\mathcal{W}\_{F}^{(1)} \, , \tag{44}$$

and generalizes Equation (39) to finite voltage and temperature biases. The right-hand side of Equation (39) can be interpreted as the dissipated energy of the mechanical system, which is delivered to the electronic reservoirs. However, in the above equation, W(1) *F* 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 (25), 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:

$$\eta = -\frac{\sum\_{|\mathbf{a}|=0}^{\mathbf{a}\_{\text{max}}} \mathcal{W}\_F^{\mathbf{a}(\Omega)}}{\sum\_{|\mathbf{a}|=0}^{\mathbf{a}\_{\text{max}}} \mathcal{Q}\_I^{\mathbf{a}(V)} \cdot \delta V}. \tag{45}$$

To illustrate this, let us take the case of a local system coupled to left and right reservoirs at voltages *VL* and *VR*, respectively. We assume the leads are at the same temperature, and we set *δVL* = −*δVR* = *δV*/2 as 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* = *VL* − *VR* and Ω up to *α*max = 2, the above general expression yields:

$$\eta = -\frac{\mathcal{W}\_F^{(0,0)} + \mathcal{W}\_F^{(1,0)} + \mathcal{W}\_F^{(0,1)}}{\left(Q\_I^{(0,0)} + Q\_I^{(1,0)} + Q\_I^{(0,1)}\right)\delta V} = -\frac{\mathcal{W}\_F^{(1,0)} + \mathcal{W}\_F^{(0,1)}}{\left(Q\_I^{(1,0)} + Q\_I^{(0,1)}\right)\delta V} \tag{46}$$

where the superscripts indicate the order in the expansions in Ω and *δV*, respectively, and we defined *I* = (*IL* − *IR*)/2. As we already mentioned, the zeroth order contributions W(0,0) *F* and *Q*(0,0) *I* 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 = 2 but only in terms of Ω, Equation (45) turns into:

$$\eta = -\frac{\mathcal{W}\_F^{(0)} + \mathcal{W}\_F^{(1)}}{\left(Q\_I^{(0)} + Q\_I^{(1)} + Q\_I^{(2)}\right) \delta V},\tag{47}$$

where now W(0) *F* and *Q*(0) *I* 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 (32) in a similar way to that of Equation (45). 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.

### **7. 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 cannot 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 [78], 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 [12,57–59,87–94] and quantum motors [42].

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

$$\frac{d}{dt}p = \mathcal{W}p.\tag{48}$$

where the vector *p* = {*pi*(*t*)} describes the dot's occupation probabilities and *W* is 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 *X* ˙ , 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 as *p* = ∑*k*≥0 *<sup>p</sup>*(*k*), with *p*(*k*) of order (Ω/Γ)*<sup>k</sup>*. 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 [12,78]:

$$\mathcal{W}p^{(0)} = \mathbf{0}, \quad \text{and} \quad \mathcal{W}p^{(k)} = \frac{\mathbf{d}}{\mathbf{d}t}p^{(k-1)}.\tag{49}$$

In the first equation, *p*(0)(*X*) represents the steady-state solution the electronic system arrives at when the mechanical system is "frozen" at the point *X*. As the mechanical coordinate moves in time, this order corresponds to the adiabatic electronic response to the mechanical motion. This needs 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 corrections *<sup>p</sup>*(*k*)(*<sup>X</sup>*, *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* = *<sup>W</sup>*(*X*). The normalization condition on the dot's density matrix implies *e*T *p*(*k*) = *δk*0, with *e*T the trace over the dot's Hilbert space. This allows for the definition of an invertible pseudo-kernel *W* ˜ , whose matrix elements are *W* ˜ *ij* = *Wij* − *Wii*, such that we can write:

$$\mathcal{p}^{(k)} = \left[\vec{\mathbf{W}}^{-1} \frac{\mathbf{d}}{\mathbf{d}t}\right]^k \mathcal{p}^{(0)}.\tag{50}$$

Once we know the different orders of the reduced density matrix, it is possible to compute the expectation value of any observable *R* after calculation of its corresponding kernel *WR*. This implies the same expansion as before, and it reads:

$$\mathcal{R}^{(k)} = \mathfrak{e}^{\mathrm{T}} \mathsf{W}^{\mathbb{R}} \mathfrak{p}^{(k)}.\tag{51}$$

The observables we are going to address here are the charge and heat currents *Ir* and *Jr* associated with the *r* lead and defined in Section 4. 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 particles *ni* and energy *Ei* associated with the dot state |*i* [87,95]:

$$\mathcal{W}\_{\rm ij}^{\rm l} = -\varepsilon (n\_{\rm i} - n\_{\rm j}) [\mathcal{W}\_{\rm r}]\_{\rm ij}, \qquad \mathcal{W}\_{\rm ij}^{\rm l} = -\left[E\_{\rm i} - E\_{\rm j} - \mu\_{\rm r} (n\_{\rm i} - n\_{\rm j})\right] [\mathcal{W}\_{\rm r}]\_{\rm ij}.\tag{52}$$

The CIF was derived in Section 2 under the assumption that the mechanical part of the system, characterized by coordinates *X* and associated momenta *P*, only interacts with local parameters of the quantum dot through their many-body eigenenergies, cf. Equation (2). This implies that the *ν*-component of the CIF operator, defined as *F* ˆ *ν* = −*∂H*<sup>ˆ</sup> s/*∂X<sup>ν</sup>*, is local in the quantum dot basis. For this local operator, then, we can simply define a diagonal kernel of the form:

$$\mathcal{W}\_{ij}^{F\_{\mathcal{V}}} = -\frac{\partial E\_i}{\partial X\_{\mathcal{V}}} \delta\_{ij\prime} \tag{53}$$

such that Equation (51) gives the *k* 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 Wij* = 0, lead to the following conservation rule for all orders in Ω [87]:

$$\sum\_{r} I\_{r}^{(k)} = -\frac{\mathbf{d}}{\mathbf{d}t} \left< \mathcal{H}\_{\mathbf{s}} \right>^{(k-1)} - F^{(k-1)} \cdot \mathbf{X} - \sum\_{r} \frac{\mu\_{r}}{\varepsilon} I\_{r}^{(k)}.\tag{54}$$

This equation, equivalent to (17), 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 reservoirs *L* and *R*, with periodic motion characterized by a time period *τ*, and taking the time integral of the above equation, we recover the frequency expansion of Equation (25):

*Q*(*k*) *J* + *Q*(*k*) *I δV* = −W(*k*−<sup>1</sup>) *F* , (55)

where the bias voltage *δV* is defined through *μL* = −*μ<sup>R</sup>* = *eδV*/2, *QJ* = *QJL* + *QJR* and *QI* = (*QIL* − *QIR* )/2. The sign convention employed here implies, for example, that if the left-hand side of Equation (55) 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.

### *7.1. 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 1c). 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 [42,87,89],

$$\hat{H}\_{\mathfrak{k}} = \sum\_{\ell=L,R} \varepsilon\_{\ell} \hat{\mathfrak{n}}\_{\ell} + \mathsf{U} \hat{\mathfrak{n}}\_{L} \hat{\mathfrak{n}}\_{R} + \frac{\mathsf{U}^{\prime}}{2} \sum\_{\ell=L,R} \hat{\mathfrak{n}}\_{\ell} (\mathfrak{n}\_{\ell} - 1) - \frac{t\_{\mathfrak{c}}}{2} \left( \hat{d}\_{L\sigma}^{\dagger} \hat{d}\_{R\sigma} + \text{h.c.} \right), \tag{56}$$

where - is the on-site energy and *n*ˆ - = ∑*σ* ˆ *d*†*σ* ˆ*dσ* the number operator in the --dot, *U* and *U* the interand intra-dot charging energies, respectively, and *tc* the interdot hopping amplitude. We will work in the strong coupling regime *tc* Γ, such that non-diagonal 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 limit *U* → <sup>∞</sup>, 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 reads *p* = (*p*0, *p*b↑, *p*b↓, *p*a↑, *p*a↓, *p*↑↑, *p*↑↓, *p*↓↑, *<sup>p</sup>*↓↓)T. These elements thus denote the probability for the DQD to be either empty (*p*0), occupied by one electron in the - = b (or a) 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 *<sup>σ</sup>*. The many-body eigenenergies are therefore *E*0 = 0 for the empty DQD,

$$E\_{\mathbf{b}/\mathbf{a}} = \frac{\boldsymbol{\varepsilon}\_{L} + \boldsymbol{\varepsilon}\_{R}}{2} \mp \sqrt{\left(\frac{\boldsymbol{\varepsilon}\_{L} - \boldsymbol{\varepsilon}\_{R}}{2}\right)^{2} + \left(\frac{t\_{c}}{2}\right)^{2}}\tag{57}$$

for single occupation in the bonding or the antibonding orbital and *E*2 =  *L* +  *R* + *U* for the doubly-occupied DQD, respectively.

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

$$\varepsilon\_L(\theta) = \overline{\epsilon}\_L + \lambda \mathbf{x}\_L = \overline{\epsilon}\_L + \delta\_\mathbf{c} \cos(\theta), \quad \text{and} \quad \varepsilon\_R(\theta) = \overline{\epsilon}\_R + \lambda \mathbf{x}\_R = \overline{\epsilon}\_R + \delta\_\mathbf{c} \sin(\theta), \tag{58}$$

which defines a circular trajectory in energy space of radius *δ* = *λ<sup>r</sup>*0 (*<sup>r</sup>*0 measures the actual radius of the rotor) around the working point (¯*L*, ¯*R*). For simplicity, in what follows, we will focus on 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 (7), we obtain the equation of motion (9) for the angular velocity ˙ *θ*, where: I = *<sup>m</sup>*eff*<sup>r</sup>*20 is the rotor's moment of inertia, F = ∑*i*(−*∂θEi*)*pi* = *<sup>r</sup>*0*F* · *θ*ˆ is the current-induced torque, Fext is 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 (11) is fulfilled. We show in Figure 2 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 average value of ˙ *θ* in one cycle. This, in turn, implies a large number of cycles for the rotor to reach the steady-state regime.

**Figure 2.** Evolution of the rotor's angular velocity (average over cycle) as a function of the number of cycles for different initial conditions. In (**<sup>a</sup>**,**b**), we choose *δ* = 100 *k*B*T* and *V* = 2 *k*B*T*, and we take in (a) Cext = 0.5 C(0) *F* , so the device operates as a motor, while in (b), we use Cext = 1.5 C(0) *F* , so the device operates as a pump, cf. Equation (62). (**c**) A situation in which two stable velocities with opposite sign are possible, for *δ* = 10 *k*B*T*, *V* = 15.42 *k*B*T*, and Cext = C(0) *F* + 0.05 *k*B*T*. The other parameters are: *δT* = 0, *U* = 20 *k*B*T*, *tc* = 10 *k*B*T*, Γ*L* = Γ*R* = 0.25 *k*B*T*, while for the trajectories, we used as working point ¯*L* = ¯*R* = −6 *k*B*T* + *δ*/√2.

As discussed in Section 3, when we take I sufficiently 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 (12). 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 defining W(*k*) *F* . To become independent of this issue, we denote by *s* the sign of ˙*θ*, and notice that F(*k*) = *f*(*k*) *θ* ˙*θk*, where *f*(*k*) *θ* = *<sup>∂</sup>k*˙*θ*F|˙*θ*=<sup>0</sup>/*k*! only depends parametrically on *θ*, so we have:

$$\mathcal{W}\_{\rm F}^{(k)} = \oint\_{C(s)} \mathbf{F}^{(k)} \cdot \mathbf{d} \mathbf{X} = \int\_0^{2\pi s} \mathcal{F}^{(k)} \mathbf{d}\theta = s \left[ \int\_0^{2\pi} f\_{\theta}^{(k)} \mathbf{d}\theta \right] \theta^k = s \mathcal{C}\_{\rm F}^{(k)} \theta^k,\tag{59}$$

where, importantly, the coefficients C(*k*) *F* are independent of *s*, so the ˙*θk* and the extra sign *s* accompanying C(*k*) *F* in the right-hand side of the equation give the correct sign for W(*k*) *F* . As discussed in Section 3, the steady-state condition in Equation (11) 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 = \sum\_{k} \mathcal{C}\_{F}^{(k)} \dot{\theta}^{k} - \mathcal{C}\_{\text{ext}\prime} \tag{60}$$

where we defined Cext = 2*π*Fext. The stability of the final solution is inherited from Equation (9), and in our case with ˙ *θ* constant, it is given by:

$$\sum\_{k} \mathcal{C}\_{F}^{(k)} k \dot{\theta}^{k-1} < 0. \tag{61}$$

Notice the specific case where Cext = C(0) *F* , i.e., the external work equals the bias work coming from the first-order currents, cf. Equation (55) for *k* = 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 (60) up to linear order in ˙ *θ*, such that the solution is:

$$\dot{\theta} = \frac{\mathcal{C}\_{\text{ext}} - \mathcal{C}\_{F}^{(0)}}{\mathcal{C}\_{F}^{(1)}},\tag{62}$$

together with the condition C(1) *F* < 0, in agreemen<sup>t</sup> with Equation (61). This implies that the integral of the friction coefficient over a cycle needs to be positive; otherwise, the rotor cannot reach the stationary regime. As the sign in C(1) *F* is fixed to this order, the sign in ˙ *θ* only depends on the relation between Cext and C(0) *F* . This, together with the stationary condition Wext = W*<sup>F</sup>*, establishes the operation mode of the device, in the sense that with this information, we can deduce the sign of Wext = *<sup>s</sup>*Cext and, hence, the direction of the energy flow. If, for example, Wext < 0, the external energy is delivered into the DQD, which in turn goes as energy current to the leads, so the device acts as an energy pump. On the other hand, if Wext > 0, the 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 2, where we take Cext > 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 in Cext where we can expect some kind of pumping mechanism. Depending then on the specific type of pumping we have in mind, the ranges in Cext will 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 in Cext where 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 to C(2) *F*in Equation (60), where we obtain:

$$\theta\_{\pm} = -\frac{\mathcal{C}\_{\text{F}}^{(1)}}{2\mathcal{C}\_{\text{F}}^{(2)}} \pm \sqrt{\left(\frac{\mathcal{C}\_{\text{F}}^{(1)}}{2\mathcal{C}\_{\text{F}}^{(2)}}\right)^{2} + \frac{\mathcal{C}\_{\text{ext}} - \mathcal{C}\_{\text{F}}^{(0)}}{\mathcal{C}\_{\text{F}}^{(2)}}},\tag{63}$$

together with the conditions:

$$
\left(\frac{\mathcal{C}\_F^{(1)}}{2\mathcal{C}\_F^{(2)}}\right)^2 > \frac{\mathcal{C}\_F^{(0)} - \mathcal{C}\_{\text{ext}}}{\mathcal{C}\_F^{(2)}}, \qquad \pm \mathcal{C}\_F^{(2)} < 0. \tag{64}
$$

The first inequality ensures a positive argumen<sup>t</sup> in the square root of Equation (63), so ˙ *θ* is a real number, while the second inequality comes from the stability condition given by Equation (61) and tells us which branch one should choose in Equation (63) to ge<sup>t</sup> the stable solution: if C(2) *F* > 0, then ˙ *θ*−, and if C(2) *F* < 0, then ˙ *θ*+. Importantly, far from equilibrium (e.g., *δV k*B*T*), it is possible to arrive at the odd situation where C(1) *F* > 0, and the "dissipated" energy W(1) *F* = *s*C(1) *F* ˙ *θ* = C(1) *F* | ˙ *θ*| is positive. This, however, does not prevent the rotor from reaching 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 favoring the motor regime), contrary to the standard situation where the "dissipated" energy flows to the leads. Regarding Equation (63), the sign of ˙ *θ* now depends on the relation between the first term and the square root, together with the sign of C(2) *F* , so this analysis is not that simple as in the linear case. However, once we have the C*F*-coefficients, we can always determine if the rotor is able to reach the stationary regime and infer whether Wext is 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 (60), it could happen that for a fixed choice of the parameters ( Fext, *δV*, *δT*, etc.), the system presents more than one stable solution, and even with different signs in Wext, so the initial condition on ˙ *θ* decides the operation mode of the device once the rotor reaches the stationary regime. In Figure 2c, we show this nonlinear 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−<sup>3</sup> *k*B*T* and 6.7 × 10−<sup>3</sup> *k*B*T*) and one unstable solution in between ( ˙ *θ* = 2 × 10−<sup>3</sup> *k*B*T*). Notice that, in this case, if the rotor starts with a positive initial condition below the unstable solution (see gray arrow), then it cannot reach the negative stable solution, as once ˙ *θ* = 0, the rotor gets trapped in a local minimum.

In Figure 3, 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 (59), 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 (50), we can estimate the maximum allowed value for ˙ *θ* compatible with the frequency expansion. Since *W*˜ −1 ∝ <sup>Γ</sup>−1, with Γ = Γ*L* + Γ*R*, and:

$$\frac{\mathbf{d}}{\delta \mathbf{d}} = \dot{\theta} \frac{\partial}{\partial \theta} = \dot{\theta} \sum\_{i} \frac{\partial \varepsilon\_{i}}{\partial \theta} \frac{\partial}{\partial \varepsilon\_{i}} \,. \tag{65}$$

we obtain:

$$\mathcal{P}^{(k)} = \left[\tilde{\mathcal{W}}^{-1}\theta \sum\_{i} \frac{\partial \varepsilon\_{i}}{\partial \theta} \frac{\partial}{\partial \varepsilon\_{i}}\right]^{k} \mathcal{P}^{(0)} \propto \left[\frac{\Omega}{\Gamma} \frac{\delta\_{c}}{k\_{\text{B}}T}\right]^{k},\tag{66}$$

where we use the fact that *∂θi* ∝ *δ* and the energy derivative of both *p*(0) and *W*˜ −1 are proportional to 1/*k*B*T*. As the consistency of the frequency expansion relies on the convergence of the occupations, this yields the adiabaticity condition [42,87]:

$$\frac{\Omega}{\Gamma} \ll \frac{k\_{\text{B}}T}{\delta\_{\text{c}}} \,\tag{67}$$

from which we can estimate the maximum allowed frequency Ωmax as:

$$
\Omega\_{\text{max}} = \frac{k\_{\text{B}}T}{\delta\_{\text{c}}} \Gamma. \tag{68}
$$

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 3, an intermediate reference frequency Ωx = Ωmax/2. As we can see, for a long range of the bias voltage ( *V* 10 *k*B*T*), 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 (62) is enough, as long as Cext ∼ C(0) *F* . For *V* - 10 *k*B*T*, we can observe non-equilibrium effects like "inverse dissipation", i.e., a positive contribution from the first-order term (red) in Figure 3a. 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 (60) or in the current-induced force appearing in Equation (9). This, in turn, could lead to several stable solutions whose validity should be determined through a systematic convergence analysis, which is beyond the scope of the present work. For the particular choice of parameters used in Figure 2c, we check if the next-order coefficient (C(4) *F* ) has a negligible impact on the third-order solutions.

**Figure 3.** Order-by-order contributions to the current-induced work as a function of the bias voltage for different orbit radius: (**a**) *δ* = 10 and Ωx = 2.5 × <sup>10</sup>−2, (**b**) *δ* = 20 and Ωx = 1.25 × <sup>10</sup>−2, (**c**) *δ* = 50 and Ωx = 5 × <sup>10</sup>−3, and (**d**) *δ* = 100 and Ωx = 2.5 × <sup>10</sup>−3, in units of *k*B*T*. We show, in solid lines, the zeroth- (black), first- (red), second- (blue) and third- (green) orders in ˙ *θ*, respectively. The circles accompanying the lines correspond to (C(*k*) *J* + C(*k*) *I δV*)Ω(*k*−<sup>1</sup>) *x* , with C(*k*) *I* and C(*k*) *J* defined in Equation (70), and numerically confirm Equations (55) and (71). The remaining parameters coincide with those of Figure 2.

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 3d) and *δV* ∼ 2 *k*B*T*. As the sign in C(1) *F* is negative, the linear-order solution given by Equation (62) is stable. If we start with Cext = 0, then ˙ *θ* is positive, so for 0 < Cext < C(0) *F* , the device acts as a motor (see Figure 2a). When Cext > C(0) *F* > 0, ˙*θ* becomes negative, and hence, the device operates as a pump (see Figure 2b). Additionally, when Cext < 0, the angular velocity is still positive, but as we changed the sign in Cext, we have that Wext < 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 Wext = *<sup>s</sup>*Cext changes with *s* and the sign of Cext. 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:

$$\mathcal{C}\_{\mathcal{F}}^{(k)}(-\delta V) = (-1)^{k+1} \mathcal{C}\_{\mathcal{F}}^{(k)}(\delta V),\tag{69}$$

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 (*<sup>L</sup>*, *R*) → (*<sup>R</sup>*, *<sup>L</sup>*), which changes the sign of ˙ *θ*, i.e., ˙ *<sup>θ</sup>*(−*δV*) = − ˙*<sup>θ</sup>*(*<sup>δ</sup>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 (69) in Equation (60). We can therefore recognize the transformation ˙ *<sup>θ</sup>*(−*δV*) = − ˙*θ*(*δV*) if we also change the sign of the external force, such that Cext(−*δV*) = −Cext(*<sup>δ</sup>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 (59), 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*) cannot be related to the left-right inversion operation.

As an additional test for Equations (55) and (25), we define equivalent coefficients for the amount of transported charge and heat in a cycle:

$$\mathcal{C}\_{I}^{(k)} = \int\_{0}^{2\pi} \left. \frac{\partial^{k} I\_{\theta}}{\partial \theta^{k}} \right|\_{\theta=0} \frac{\mathbf{d}\theta}{k!} , \qquad \mathcal{C}\_{I}^{(k)} = \int\_{0}^{2\pi} \left. \frac{\partial^{k} I\_{\theta}}{\partial \theta^{k}} \right|\_{\theta=0} \frac{\mathbf{d}\theta}{k!} , \tag{70}$$

such that *Q*(*k*) *I* = *s*C(*k*) *I* ˙*<sup>θ</sup>k*−1, and *Q*(*k*) *J* = *s*C(*k*) *J* ˙*θk*−1. These contributions can be evaluated independently from the C*F*-coefficients, and in Figure 3, these are shown in circles, which gives numerical agreemen<sup>t</sup> for the energy conservation principle:

$$\mathcal{C}\_{I}^{(k)} + \mathcal{C}\_{I}^{(k)} \delta V = -\mathcal{C}\_{F}^{(k-1)},\tag{71}$$

in the considered orders of *θ*. For the considered example in Figure 3, the current coefficients C*I* and C*J* also show (independently) a definite parity with respect to the bias voltage. In fact, as Equation (71) suggests, the heat current coefficients C(*k*) *J* present the same parity as C(*k*−<sup>1</sup>) *F* , while the charge current coefficients C(*k*) *I*need to have the opposite parity since they are multiplied by *δV*.
