Thermodynamics and Steady State of Quantum Motors and Pumps Far from Equilibrium

Abstract

In this article, we briefly review the 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.

1. Introduction

In recent years, there has been a sustained growth in the interest in different forms of nanomachines. This was boosted by seminal experiments [1,2,3,4,5,6,7,8], the blooming of new theoretical proposals [9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30], and the latest developments towards the understanding of the fundamental physics underlying such systems [31,32,33,34,35,36,37,38,39,40,41,42,43,44,45]. 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” [13,46,47,48,49,50,51,52,53]. They typically consist of an electromechanical device connected to electronic reservoirs and controlled by nonequilibrium sources; see Figure 1. 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 [13,39,54], Thouless quantum pumps and motors [13,22,55], or Anderson quantum motors [56], 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 [26,42,57,58,59]. The effect of decoherence has also been addressed [22,39,54,60], as well as the influence of the friction forces and the system-lead coupling in the dynamics of quantum motors and pumps [22]. Indeed, the thermodynamics of those systems has proven to be a key aspect to study. In the last few years, different individual efforts have coalesced to give rise to a new field dubbed “quantum thermodynamics” [48,49,50,51,52,53,61], 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 [48,49,50,51,52,53,61]. 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 [52]. 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 [52] and the references therein. However, such parameterizations usually hide the physics behind the nonlinearities and require the optimization of a number of variables that grow 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 [48] and the 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 connections 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 the 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 [10,11,22,24,26,34,35,37,39,42,44,49,50,51,62,63,64,65,66,67,68,69,70,71]. Because of the delayed response of the electronic degrees of freedom to the mechanical motion, one should include the so-called nonadiabatic corrections with 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 it is adequate or not to assume a constant terminal velocity [22,39,42]. 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, especially 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 in 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 2, 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 3, 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 4, 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 an 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 5, 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 6, we analyze and give physical interpretation to some of the relations obtained in Section 4. Finally, in Section 7, 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.

2. 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 a starting point, we consider as the local system the region where electronic and mechanical degrees of freedom are present and coupled to each other, like the examples shown in Figure 1. Such a local system is generically modeled by the following Hamiltonian:

$${\widehat{H}}_{\mathrm{local}}={\widehat{H}}_{\mathrm{s}}(\widehat{\mathit{X}})+\frac{{\widehat{\mathit{P}}}^2}{2m_{\mathrm{eff}}}+U(\widehat{\mathit{X}},t),$$

(1)

where $\widehat{\mathit{X}}=({\widehat{X}}_1,\dots ,{\widehat{X}}_N)$ is the vector of mechanical coordinates and $\widehat{\mathit{P}}=({\widehat{P}}_1,\dots ,{\widehat{P}}_N)$ collects their associated momenta, $m_{\mathrm{eff}}$ is the effective mass related to $\widehat{\mathit{X}}$, and $U(\widehat{\mathit{X}},t)$ represents some external potential, of a 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 ${\widehat{H}}_{\mathrm{s}}$ includes both the electronic degrees of freedom and their coupling to the mechanical ones through:

$${\widehat{H}}_{\mathrm{s}}(\widehat{\mathit{X}})=\sum _i{E}_i(\widehat{\mathit{X}})|i\rangle \langle i|,$$

(2)

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

$${\widehat{H}}_{\mathrm{total}}={\widehat{H}}_{\mathrm{local}}+\sum _r{\widehat{H}}_r+\sum _r{\widehat{H}}_{\mathrm{s},r}.$$

(3)

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

$${\widehat{H}}_r=\sum _{k\sigma }{ϵ}_{rk}{\widehat{c}}_{rk\sigma }^{†}{\widehat{c}}_{rk\sigma },$$

(4)

where ${\widehat{c}}_{rk\sigma }^{†}$ (${\widehat{c}}_{rk\sigma }$) creates (annihilates) an electron in the r-reservoir with state-index k and spin projection $\sigma$. As usual, the reservoirs are assumed to be always in equilibrium, characterized by a temperature $T_r$ and electrochemical potential ${\mu }_r$. The coupling between the local system and the r-lead is determined by the tunnel Hamiltonian:

$${\widehat{H}}_{\mathrm{s},r}=\sum _{k\sigma \ell }\left(t_{r\ell }{\widehat{d}}_{\ell \sigma }^{†}{\widehat{c}}_{rk\sigma }+\mathrm{h}.\mathrm{c}.\right),$$

(5)

where $t_{r\ell }$ denotes the tunnel amplitude, assumed to be k and $\sigma$ independent for simplicity, and the fermion operator ${\widehat{d}}_{\ell \sigma }^{†}$ (${\widehat{d}}_{\ell \sigma }$) creates (annihilates) one electron with spin $\sigma$ in the ℓ-orbital of the local system. The tunnel-coupling strengths ${\Gamma }_{r\ell }=2\pi {\rho }_r{\left|t_{r\ell }\right|}^2$ then characterize the rate at which the electrons enter/leave the local system from/to the r-reservoir, where ${\rho }_r$ is the density of states in the r-lead. Note that ${\widehat{H}}_{\mathrm{s}}$ is defined in the eigenstate basis, while ${\widehat{H}}_{\mathrm{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 [72].

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 $\widehat{\mathit{P}}$ operator, which yields:

$$m_{\mathrm{eff}}\frac{\mathrm{d}\dot{\widehat{\mathit{X}}}}{\mathrm{d}t}+\widehat{\nabla U}(\widehat{\mathit{X}},t)=-{\widehat{\nabla H}}_{\mathrm{s}}(\widehat{\mathit{X}}).$$

(6)

The measured value $A_{\mathrm{measured}}$ of an observable described by an operator $\widehat{A}$ can always be taken as its mean value $A=\langle \widehat{A}\rangle$ plus some fluctuation ${\xi }_A$ around it, i.e., $A_{\mathrm{measured}}=A+{\xi }_A$. We will work under the nonequilibrium Born–Oppenheimer approximation [13,22,34,35,37,39,49,73,74] (or Ehrenfest approximation [33,63,75,76,77]), 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 (6) and describe the mechanical motion only through the mean value $\mathit{X}$, which is reasonable for large or massive objects. With this in mind, we obtain the following Langevin equation of motion:

$$m_{\mathrm{eff}}\frac{\mathrm{d}\dot{\mathit{X}}}{\mathrm{d}t}+{\mathit{F}}_{\mathrm{ext}}=\mathit{F}+\mathit{\xi },$$

(7)

where $\mathit{F}=-〈{\widehat{\nabla H}}_{\mathrm{s}}〉=i〈[{\widehat{H}}_{\mathrm{s}}(\widehat{\mathit{X}}),\widehat{\mathit{P}}]〉$ and $\mathit{\xi }$ account for the mean value and the fluctuation of the CIF, respectively (throughout this manuscript, we take $\hslash =1$ for simplicity). As we shall see later on, the external force applied to the mechanical part of the local system, ${\mathit{F}}_{\mathrm{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 define ${\mathit{F}}_{\mathrm{ext}}$ with a minus sign for better clarity in future discussions. The main task, therefore, relies on the calculation of the CIFs from appropriate formalisms capable of describing the dynamics of the electronic part of the system. Once these forces are calculated, we can use Equation (7) to integrate the classical equations of motion and obtain the effective dynamics of the complete electromechanical system.

In most previous works, $\mathit{F}$ is expanded up to first order in $\dot{\mathit{X}}$, i.e., $\mathit{F}\approx {\mathit{F}}^{(0)}-\mathit{\gamma }·\dot{\mathit{X}}$. The resulting CIF is then the sum of an adiabatic contribution ${\mathit{F}}^{(0)}$ and its first nonadiabatic correction ${\mathit{F}}^{(1)}=-\gamma \dot{\mathit{X}}$, respectively. Under this approximation, Equation (7) turns into:

$$m_{\mathrm{eff}}\frac{\mathrm{d}\dot{\mathit{X}}}{\mathrm{d}t}+{\mathit{F}}_{\mathrm{ext}}={\mathit{F}}^{(0)}-\mathit{\gamma }·\dot{\mathit{X}}+\mathit{\xi }.$$

(8)

Explicit formulas for the calculation of ${\mathit{F}}^{(0)}$, $\mathit{\gamma }$, and $\mathit{\xi }$ in terms of Green functions and scattering matrices were derived in [10,11,34,35,37] and extended in [22,39] 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 [62,64]. In [49,50,51], 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 [44]; modeling the system as a Luttinger liquid [24]; and using a time-dependent slave-boson approximation [26]. In [42], explicit expressions for the CIFs within the Coulomb blockade regime of transport were obtained using a real-time diagrammatic approach [78], which we present in more detail in Section 7 when considering the example of Figure 1c.

3. Mechanical Steady State

In this paper, we will restrict ourselves to systems that perform overall cyclic motions. Immediate examples are shown in Figure 1a,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 1b,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 1b, 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 1d, 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., [79,80,81]. 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 (7) on a closed trajectory defined in the space of $\mathit{X}$. By assuming a circular trajectory, the dynamics can be described by an angle $\theta$, its associated angular velocity $\dot{\theta }$, the moment of inertia $\mathcal{I}$, and the torques $\mathcal{F}$, ${\mathcal{F}}_{\mathrm{ext}}$, and ${\xi }_{\theta }$. Using this, we obtain an effective angular Langevin equation equivalent to Equation (7),

$$\frac{\mathrm{d}\dot{\theta }}{\mathrm{d}t}=\frac{1}{\mathcal{I}}\left[\mathcal{F}-{\mathcal{F}}_{\mathrm{ext}}+{\xi }_{\theta }\right].$$

(9)

We assume that after a long waiting time, the system arrives at the steady-state regime where the mechanical motion becomes periodic, and it is then characterized by a time period $\tau$ such that $\theta (t+\tau )=\theta (t)$ and $\dot{\theta }(t+\tau )=\dot{\theta }(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 $\theta$ and $\dot{\theta }$, i.e., the mean trajectories with or without the stochastic force approximately coinciding. This occurs, for example, at low temperatures or in mechanical systems with a large moment of inertia [22,39,42]. In the following, we will just ignore ${\xi }_{\theta }$ for practical purposes (this is the opposite regime of another type of nanomotors, the Brownian motors [1,82]). Under the above assumptions, we integrate both sides of Equation (9) from an initial position ${\theta }_i$ to a final one ${\theta }_f$ and obtain:

$$\begin{array}{ccc}\hfill \frac{\mathcal{I}}{2}\left[{\dot{\theta }}_f^2-{\dot{\theta }}_i^2\right]& =& {\int }_{{\theta }_i}^{{\theta }_f}\left[\mathcal{F}-{\mathcal{F}}_{\mathrm{ext}}\right]\mathrm{d}\theta .\hfill \end{array}$$

(10)

The torques in this equation are, in general, intricate functions of both $\theta$ and $\dot{\theta }$.

Therefore, the calculation of the $\theta$-dependent angular velocity usually requires the resolution of a transcendental equation (see, e.g., [42]). Alternatively, one can obtain $\theta (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 of taking the terminal velocity as constant during the whole cycle [13,22,24,26,39,49,50,51,56,74]. 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 $\dot{\theta }$ in one period. We now address this interesting issue in more detail. First, we take the integral in Equation (10) over the whole period. This gives:

$${\mathcal{W}}_{\mathrm{ext}}={\mathcal{W}}_F,\mathrm{where}{\mathcal{W}}_F={\int }_0^{\tau }\mathcal{F}\dot{\theta }\mathrm{d}t,\mathrm{and}{\mathcal{W}}_{\mathrm{ext}}={\int }_0^{\tau }{\mathcal{F}}_{\mathrm{ext}}\dot{\theta }\mathrm{d}t.$$

(11)

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 that 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 $\dot{\theta }$ for later when treating a concrete example in Section 7.1). If we now expand $\mathcal{F}$ in terms of $\dot{\theta }$, Equation (11) yields:

$$\begin{array}{ccc}\hfill {\mathcal{W}}_{\mathrm{ext}}& =& \sum _k\left({\int }_0^{2\pi }{\left.\frac{\partial {}^k\mathcal{F}}{\partial {\dot{\theta }}^k}\right|}_{\dot{\theta }=0}\frac{\mathrm{d}\theta }{k!}\right){\dot{\theta }}^k.\hfill \end{array}$$

(12)

Two important conclusions can be extracted from the above formal solution. First, there may be conditions where some roots of Equation (12) are complex numbers, meaning that the assumption $\dot{\theta }=\mathrm{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 [22,42] that the moment of inertia $\mathcal{I}$ not only affects the time that it takes the mechanical system to reach the stationary regime, but also the internal range in the angular velocity, i.e., the difference $\Delta \dot{\theta }={\dot{\theta }}_{\max }-{\dot{\theta }}_{\min }$ in one period. According to Equation (12), $\dot{\theta }$ is independent of $\mathcal{I}$, while the variation of $\dot{\theta }$ scales with ${\mathcal{I}}^{-1}$, cf. Equation (10). Then, the ratio $\Delta \dot{\theta }/\dot{\theta }$, which is the relevant quantity in our analysis, should vanish for large $\mathcal{I}$ values, justifying the constant velocity assumption for large or massive mechanical systems.

The value of ${\mathcal{F}}_{\mathrm{ext}}$ is supposed to be controllable externally, as well as the voltage and temperature biases, which, in turn, affect the current-induced torque $\mathcal{F}$. Therefore, under the above discussed conditions, $\dot{\theta }$ 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 voltage $eV={\mu }_L-{\mu }_R$. By considering the current-induced torque up to its first nonadiabatic correction, i.e., $\mathcal{F}\approx {\mathcal{F}}^{(0)}-\gamma \dot{\theta }$, and assuming that ${\mathcal{F}}_{\mathrm{ext}}$ is independent of $\dot{\theta }$, the following relation must hold, according to the above discussion,

$$\dot{\theta }\approx \frac{Q_{I_R}}{2\pi \overline{\gamma }}\left(V-\frac{{\mathcal{W}}_{\mathrm{ext}}}{Q_{I_R}}\right)=\frac{Q_{I_R}{V}_{\mathrm{eff}}}{2\pi \overline{\gamma }},$$

(13)

where $\overline{\gamma }$ is the average electronic friction coefficient along the cycle, $Q_{I_R}$ is the pumped charge to the right lead, and we used Onsager’s reciprocal relation between $\mathcal{F}$ and the charge current $I_R$ in the absence of magnetic fields [22,39,42,49]. Alternatively, if we assume that the external torque is of the form ${\mathcal{F}}_{\mathrm{ext}}={\gamma }_{\mathrm{ext}}\dot{\theta }$, one finds:

$$\dot{\theta }\approx \frac{Q_{I_R}V}{2\pi }\frac{1}{(\overline{\gamma }+{\overline{\gamma }}_{\mathrm{ext}})}=\frac{Q_{I_R}V}{2\pi {\gamma }_{\mathrm{eff}}}.$$

(14)

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 voltage V or the electronic friction coefficient $\gamma$. Numerical simulations in [22,42] showed that the 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 cannot be compensated by the work done by the CIF, and thus, $\dot{\theta }=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 ($\dot{\theta }=0$ and those given by the above equations). Although in Section 7.1, we will take $\dot{\theta }$ as constant when discussing a specific example, the reader should keep in mind that this approximation does not always hold, especially at very small voltages or $\mathcal{I}$.

4. 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 allowed 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 on systematic 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:

$$\widehat{H}=\sum _r{\widehat{H}}_r+{\widehat{H}}_{\mathrm{s}}+\sum _r{\widehat{H}}_{\mathrm{s},r}.$$

(15)

The mechanical part of the local system, when treated classically in Equation (7), introduces an explicit time dependence into ${\widehat{H}}_{\mathrm{s}}$ that, in turn, makes $\mathrm{d}\langle \widehat{H}\rangle /\mathrm{d}t\ne 0$. According to the above Hamiltonian, the total internal energy of the electronic system, $U=\langle \widehat{H}\rangle$, 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 is:

$${\dot{U}}_{\beta }=J_{\beta }^E+〈\frac{\partial {\widehat{H}}_{\beta }}{\partial t}〉,\mathrm{where}{J}_{\beta }^E=i\langle \left[\widehat{H},{\widehat{H}}_{\beta }\right]\rangle$$

(16)

which is the mean value of the energy flux entering in the subsystem $\beta =\{r,(\mathrm{s},r),\mathrm{s}\}$. The value of $\langle {\partial }_t{\widehat{H}}_{\beta }\rangle$ is zero when we evaluate the energy flux in the reservoirs and the tunnel couplings, but equals $-{\dot{\mathcal{W}}}_F$ when $\beta$ is evaluated in the local system. Although, strictly speaking, ${\mathcal{W}}_F$ is time independent when considering a closed trajectory, we used the symbol ${\dot{\mathcal{W}}}_F$ to denote the power delivered by the CIF, i.e., ${\dot{\mathcal{W}}}_F=-\langle {\partial }_t{\widehat{H}}_{\mathrm{s}}\rangle =\mathit{F}·\dot{\mathit{X}}$. The latter comes from the definition of the CIF given in Equation (7). Note that the energy fluxes fulfill the condition ${\sum }_r(J_r^E+J_{\mathrm{s},r}^E)+J_{\mathrm{s}}^E=0$. Therefore, the variation of the total internal energy of the electrons yields the following conservation rule:

$$\dot{U}=\sum _r{\dot{U}}_r+\sum _r{\dot{U}}_{\mathrm{s},r}+{\dot{U}}_{\mathrm{s}}=\sum _r\left(J_r^E+J_{\mathrm{s},r}^E\right)+J_{\mathrm{s}}^E-{\dot{\mathcal{W}}}_F=-{\dot{\mathcal{W}}}_F.$$

(17)

Conservation of the total number of particles implies a relation between the particle currents of the reservoirs ${\dot{N}}_r$ and that of the local system ${\dot{N}}_{\mathrm{s}}$:

$${\dot{N}}_{\mathrm{total}}=\frac{\mathrm{d}}{\mathrm{d}t}\langle {\widehat{N}}_{\mathrm{total}}\rangle =\sum _r\frac{\mathrm{d}}{\mathrm{d}t}\langle {\widehat{N}}_r\rangle +\frac{\mathrm{d}}{\mathrm{d}t}\langle {\widehat{N}}_{\mathrm{s}}\rangle =0\Rightarrow \sum _r{\dot{N}}_r=-{\dot{N}}_{\mathrm{s}}.$$

(18)

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 (17) over a time period $\tau$ of the cyclic motion. Under this condition, the above-defined local quantities only depend periodically on time, i.e., $U_{\mathrm{s}}(t+\tau )=U_{\mathrm{s}}(t)$, $U_{\mathrm{s},r}(t+\tau )=U_{\mathrm{s},r}(t)$, and ${\dot{N}}_{\mathrm{s}}(t+\tau )={\dot{N}}_{\mathrm{s}}(t)$. Therefore, the following quantities should evaluate to zero, i.e.,

$${\int }_0^{\tau }{\dot{U}}_{\mathrm{s},r}\mathrm{d}t=0,{\int }_0^{\tau }{\dot{U}}_{\mathrm{s}}\mathrm{d}t=0,\mathrm{and}{\int }_0^{\tau }{\dot{N}}_{\mathrm{s}}\mathrm{d}t=0,$$

(19)

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 (19) can be used together with Equation (18) to prove the charge current conservation between reservoirs, ${\sum }_r{\int }_0^{\tau }{I}_r\mathrm{d}t=0$, where the charge current of the r reservoir is defined as $I_r=e{\dot{N}}_r$, with $e>0$ being minus the electron’s charge.

We will take the following definition for the heat current $J_r$ in the reservoir r:

$$J_r=J_r^E-{\mu }_r{\dot{N}}_r.$$

(20)

In [43,50,83], the authors proposed a different definition for the heat current of the reservoirs, $J_r=J_r^E+\left(J_{\mathrm{s},r}^E/2\right)-{\mu }_r{\dot{N}}_r$. However, the inclusion or not of half the heat current of the coupling region, $\left(J_{\mathrm{s},r}^E/2\right)$, does not make any difference in the present paper, as this quantity integrates to zero over a cycle, Equation (19), and does not contribute to the rate of entropy production, as we will see in the next section, Equation (30). Replacing Equation (20) in Equation (17) and integrating over a period result in:

$$\sum _r\left(Q_{J_r}+Q_{I_r}\delta V_r\right)=-{\mathcal{W}}_F,\mathrm{where}{Q}_{J_r}={\int }_0^{\tau }{J}_r\mathrm{d}t,\mathrm{and}{Q}_{I_r}={\int }_0^{\tau }{I}_r\mathrm{d}t,$$

(21)

where $\delta V_r=({\mu }_r-{\mu }_0)/e$ and ${\mu }_0$ is an arbitrary reference’s potential. We, in addition, defined the quantities $Q_{J_r}$ and $Q_{I_r}$, which are, respectively, the total heat and charge pumped to reservoir r in a cycle. Note that in the above equation, we used conservation of the total charge.

Before we 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 $\chi$, which parameterizes the closed mechanical trajectory. This, for example, can be accounted for by taking $\chi =\Omega \mathrm{mod}(t,\tau )$, such that $\dot{\chi }=\Omega =2\pi /\tau$. In other words, one can always find a natural scale, $\Omega$ 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 reads ${\mathcal{W}}_F={\int }_0^{\tau }{F}_{\Omega }\Omega \mathrm{d}t$, where $F_{\Omega }=\mathit{F}·{\partial }_{\chi }\mathit{X}$. Of course, with this $\chi$-parameterization, the velocity is constant by definition, but at the expense that now the calculation of $F_{\Omega }$ requires the knowledge of the mechanical trajectory. For circular trajectories and the conditions discussed in Section 3, the coordinate $\chi$ would be $\left|\theta \right|$. For other cases, it would not be that simple to find $\chi$, and one should first solve the system’s equation of motion.

In addition to the mechanical velocity $\Omega$, the CIF and the currents can also be expanded in terms of the remaining nonequilibrium sources, corresponding to voltage and temperature deviations $\delta \mathit{V}={(\delta V_1,\delta V_2,\dots )}^{\mathrm{T}}$ and $\delta \mathit{T}={(\delta T_1,\delta T_2,\dots )}^{\mathrm{T}}$, respectively, from their equilibrium values $V_{\mathrm{eq}}$ and $T_{\mathrm{eq}}$. For an arbitrary observable R, this general expansion takes the form:

$$R=\sum _{\left|\mathit{\alpha }\right|\ge 0}{R}^{\mathit{\alpha }}=\sum _{\left|\mathit{\alpha }\right|\ge 0}\frac{{(\Omega ,\delta \mathit{V},\delta \mathit{T})}^{\mathit{\alpha }}}{\mathit{\alpha }!}{\left({\partial }^{\mathit{\alpha }}R\right)}_{\mathrm{eq}},$$

(22)

where we used the following multi-index notation: $\mathit{\alpha }=\left(n_{\Omega },n_{V_1},\dots ,n_{T_1},\dots \right)$; $\mathit{\alpha }!=n_{\Omega }!n_{V_1}!\dots n_{T_1}!\dots$; ${(\Omega ,\delta \mathit{V},\delta \mathit{T})}^{\mathit{\alpha }}={\Omega }^{n_{\Omega }}\delta V_1^{n_{V_1}}\dots \delta T_1^{n_{T_1}}\dots$; and:

$${\partial }^{\mathit{\alpha }}=\frac{{\partial }^{\left(n_{\Omega }+n_{V_1}+...n_{T_1}+...\right)}}{{\partial \Omega }^{n_{\Omega }}\partial V_1^{n_{V_1}}\dots \partial T_1^{n_{T_1}}\dots }.$$

(23)

This, in turn, allowed us to recognize in the integral quantities of Equation (21) the following expansion coefficients:

$${\mathcal{W}}_F^{\mathit{\alpha }}={\int }_0^{\tau }{F}_{\Omega }^{\mathit{\alpha }}\Omega \mathrm{d}t,Q_{J_r}^{\mathit{\alpha }}={\int }_0^{\tau }{J}_r^{\mathit{\alpha }}\mathrm{d}t,\mathrm{and}{Q}_{I_r}^{\mathit{\alpha }}={\int }_0^{\tau }{I}_r^{\mathit{\alpha }}\mathrm{d}t,$$

(24)

where $F_{\Omega }^{\mathit{\alpha }}$, $J_r^{\mathit{\alpha }}$, and $I_r^{\mathit{\alpha }}$ all have the form given by Equation (22). To find a consistent order-by-order conservation relation from Equation (21), we should first note that the involved terms enter in different orders: ${\mathrm{W}}_F^{\mathit{\alpha }}$ contains an additional $\Omega$ from the time-integral as compared with the pumped currents, while the term $Q_{I_r}^{\mathit{\alpha }}\delta V_r$ is one order higher in $\delta V_r$ than the other two. Therefore, the following relation must hold for every order $\mathit{\alpha }$ of the expansion:

$${\mathit{Q}}_J^{\mathit{\alpha }}·\mathbf{1}+{\mathit{Q}}_I^{\mathit{\alpha }(\mathit{V})}·\delta \mathit{V}=-{\mathcal{W}}_F^{\mathit{\alpha }(\Omega )},$$

(25)

where $\mathbf{1}={(1,1,\dots )}^{\mathrm{T}}$, ${\mathit{Q}}_J^{\mathit{\alpha }}=(Q_{J_1}^{\mathit{\alpha }},Q_{J_2}^{\mathit{\alpha }},\dots )$, ${\mathit{Q}}_I^{\mathit{\alpha }(\mathit{V})}=(Q_{I_1}^{\mathit{\alpha }(V_1)},Q_{I_2}^{\mathit{\alpha }(V_2)},\dots )$, and we used the shorthand $\mathit{\alpha }(\Omega )=(n_{\Omega }-1,n_{V_1},\dots ,n_{T_1},\dots )$ and $\mathit{\alpha }(V_r)=(n_{\Omega },n_{V_1},\dots ,n_{V_r}-1,\dots ,n_{T_1},\dots )$. Equation (25) results in being 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 $\Omega$, $\delta \mathit{V}$, and $\delta \mathit{T}$. Importantly, the above conservation rule can be extended to other nonequilibrium sources like, e.g., spin polarization in ferromagnetic leads, such that other types of currents could also be considered in the energy transfer process between the electronic and mechanical parts of the system.

5. Entropy Production And Efficiency

In 1865, Rudolf Clausius proposed a new state function, the thermodynamic entropy S, 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 $\delta Q_J$ that is transferred in a reversible thermodynamic process, $\delta S=\delta Q_{\mathrm{rev}}/T$. Here, we are assuming that each reservoir r is in its own equilibrium at a constant temperature $T_r$ and chemical potential ${\mu }_r$ (we are not going to treat time-dependent $T_r$ or ${\mu }_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., $\delta Q_{J_r}=T_r\delta S_r$. Therefore, from Equation (16) and the definition of the heat current given in Equation (20), we can write:

$${\dot{U}}_r=T_r{\dot{S}}_r+{\dot{N}}_r{\mu }_r.$$

(26)

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., [84,85,86]. 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.

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

$$\sum _r{\dot{U}}_r=\sum _r{T}_r{\dot{S}}_r+\sum _r{\dot{N}}_r{\mu }_r.$$

(27)

If we take $T_r=\delta T_r+T_0$ and ${\mu }_r={\mu }_0+\delta {\mu }_r$, add and subtract the change of the internal energy of the local system ${\dot{U}}_s$ and that of the couplings between the local system and the reservoirs ${\sum }_r{\dot{U}}_{\mathrm{s},r}$, one can rearrange the above equation to the following:

$$T_0\sum _r{\dot{S}}_r=\left({\dot{U}}_{\mathrm{s}}+\sum _r{\dot{U}}_r+\sum _r{\dot{U}}_{\mathrm{s},r}\right)-\sum _r{\dot{N}}_r\delta {\mu }_r-\sum _r{\dot{S}}_r\delta T_r-\left({\dot{U}}_{\mathrm{s}}+\sum _r{\dot{U}}_{\mathrm{s},r}\right)-{\mu }_0\left(\sum _r{\dot{N}}_r\right).$$

(28)

Note that the values of ${\mu }_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.

Replacing ${\dot{S}}_r$ by $=J_r/T_r$ in the right-hand side of the above equation, using energy conservation (17) and particle number conservation (18), allows one to rewrite Equation (28) as:

$$\begin{array}{ccc}\hfill T_0{\dot{S}}_{\mathrm{res}}& =& -F_{\Omega }\Omega -\sum _r{I}_r\delta V_r-\sum _r{J}_r\left(\frac{\delta T_r}{T_r}\right)-\left({\dot{U}}_{\mathrm{s}}+\sum _r{\dot{U}}_{\mathrm{s},r}\right)+{\mu }_0{\dot{N}}_{\mathrm{s}},\hfill \end{array}$$

(29)

where ${\dot{S}}_{\mathrm{res}}={\sum }_r{\dot{S}}_r$ is the variation of the entropy of the electrons of all reservoirs, and we used $e{\dot{N}}_r=I_r$. The CIF can be split into “equilibrium” and “nonequilibrium” terms, $F^{(\mathrm{eq})}$ and $F^{(\mathrm{ne})}$, respectively, where one can prove that $F^{(\mathrm{eq})}$ is always conservative [13,22,42,63]. We are interested in the steady-state situation of our local system. As discussed around Equation (19), 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:

$${\dot{S}}_{\mathrm{res}}^{\left(\mathrm{rev}\right)}=-\frac{1}{T_0}\left(F_{\Omega }^{(\mathrm{eq})}\Omega +{\dot{U}}_{\mathrm{s}}+\sum _r{\dot{U}}_{\mathrm{s},r}-{\mu }_0{\dot{N}}_{\mathrm{s}}\right).$$

(30)

Obviously, this quantity will not contribute to the total entropy production. Therefore, the rate of entropy production ${\dot{S}}_{\mathrm{res}}^{\left(\mathrm{irrev}\right)}$ yields:

$${\dot{S}}_{\mathrm{res}}^{(\mathrm{irrev})}=-\frac{1}{T_0}\left(F^{(}{\mathrm{ne})}_{\Omega }\Omega -\mathit{I}·\delta \mathit{V}-\mathit{J}·\delta \mathcal{T}\right),$$

(31)

where $\delta \mathcal{T}={(\delta T_1/T_1,\delta T_2/T_2,\dots )}^{\mathrm{T}}$, and the currents are defined through $\mathit{I}=(I_1,I_2,\dots )$ and $\mathit{J}=(J_1,J_2,\dots )$. Integrating Equation (31) over a cycle and taking into account the second law of thermodynamics, one finds:

$$0\ge {\mathcal{W}}_F+{\mathit{Q}}_I·\delta \mathit{V}+{\mathit{Q}}_J·\delta \mathcal{T}.$$

(32)

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 $\delta \mathcal{T}=0$, ${\mathit{Q}}_I·\delta \mathit{V}<0$, and ${\mathcal{W}}_F>0$, then the system should operate as a nanomotor driven by electric currents, and the following relation holds:

$$1\ge -\frac{{\mathcal{W}}_F}{{\mathit{Q}}_I·\delta \mathit{V}},$$

(33)

while for ${\mathit{Q}}_I·\delta \mathit{V}>0$ and ${\mathcal{W}}_F<0$, the system operates as a charge pump, and Equation (32) implies:

$$1\ge -\frac{{\mathit{Q}}_I·\delta \mathit{V}}{{\mathcal{W}}_F}.$$

(34)

Notice that, because of the steady-state condition, ${\mathcal{W}}_F$ equals ${\mathcal{W}}_{\mathrm{ext}}$, where ${\mathcal{W}}_{\mathrm{ext}}$ can be taken as the output or the input energy, depending on the considered type of process. Therefore, the above formulas describe the efficiency $\eta$ 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 $\delta \mathit{V}=0$ and $\delta \mathcal{T}\ne 0$, where Equation (32) yields:

$$1\ge -\frac{{\mathcal{W}}_F}{{\mathit{Q}}_J·\delta \mathcal{T}},\mathrm{and}1\ge -\frac{{\mathit{Q}}_J·\delta \mathcal{T}}{{\mathcal{W}}_F}.$$

(35)

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 $\delta \mathcal{T}$, the above formulas differ from what is expected from energy conservation solely. This is clear in a two-lead system, where $\eta$ is limited by 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:

$${\mathit{Q}}_J·\delta \mathcal{T}=Q_{J_{\mathrm{hot}}}\left(1-\frac{T_{\mathrm{cold}}}{T_{\mathrm{hot}}}\right)<0\Rightarrow \left(1-\frac{T_{\mathrm{cold}}}{T_{\mathrm{hot}}}\right)\ge -\frac{{\mathcal{W}}_F}{Q_{J_{\mathrm{hot}}}},$$

(36)

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 (32) to set the bounds of their associated efficiencies.

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 $\left|\mathit{\alpha }\right|=0$. Importantly, Equation (25) with $\left|\mathit{\alpha }\right|=0$ seems to impose the evaluation of terms with negative coefficients. For example, in the exponent of ${\mathit{Q}}_I$, one would be tempted to evaluate $\mathit{\alpha }(V_1)=(0,-1,0,0\dots )$, but this term is not defined in the expansions of Equation (24), and then, it should be taken as zero. Therefore, Equation (25) yields,

$${\mathit{Q}}_J^{\mathbf{0}}·\mathbf{1}=\sum _r{\int }_0^{\tau }{\left(J_r\right)}_{\mathrm{eq}}\mathrm{d}t=0.$$

(37)

This simply reflects the fact that no net heat current occurs at equilibrium. For $\left|\mathit{\alpha }\right|=1$, all relations are summarized in the following three cases:

$$\sum _r{\int }_0^{\tau }{\left(\frac{\partial J_r}{\partial \Omega }\right)}_{\mathrm{eq}}\Omega \mathrm{d}t=0,\sum _r{\int }_0^{\tau }{\left(\frac{\partial J_r}{\partial V_i}\right)}_{\mathrm{eq}}\delta V_i\mathrm{d}t=0,\mathrm{and}\sum _r{\int }_0^{\tau }{\left(\frac{\partial J_r}{\partial T_i}\right)}_{\mathrm{eq}}\delta T_i\mathrm{d}t=0,$$

(38)

where we used the fact that equilibrium forces are conservative [13,22,42,50,63], i.e., ${\int }_0^{\tau }{F}_{\Omega }^{(}\mathrm{eq})\Omega \mathrm{d}t=0$, and that there are no net charge currents in equilibrium, thus ${\int }_0^{\tau }{(I_r)}_{\mathrm{eq}}\mathrm{d}t=0$. The above relations mean that, at first order, there is a conservation of pumped heat between reservoirs. For $\left|\mathit{\alpha }\right|=2$, there are many relations and cases, but we restrict ourselves to only a few of them. For $n_{\Omega }=2$, $n_{V_i}=0$, and $n_{T_j}=0$, Equation (25) gives:

$$\begin{array}{ccc}\hfill \sum _r{0}^{\tau }\frac{{\Omega }^2}{2}{\left(\frac{{\partial }^2{J}_r}{\partial {\Omega }^2}\right)}_{\mathrm{eq}}\mathrm{d}t& =& {-}_0^{\tau }{\left(\frac{\partial F_{\Omega }}{\partial \Omega }\right)}_{\mathrm{eq}}{\Omega }^2\mathrm{d}t.\hfill \end{array}$$

(39)

Now, the quantity ${({\partial }_{\Omega }{F}_{\Omega })}_{\mathrm{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_{\Omega }=0$, $n_{V_i}=1$, $n_{V_j}=1$, and $n_{T_k}=0$, Equation (25) yields:

$$\begin{array}{ccc}\hfill \sum _r{0}^{\tau }\delta V_i\delta V_j{\left(\frac{{\partial }^2{J}_r}{\partial V_i\partial V_j}\right)}_{\mathrm{eq}}\mathrm{d}t& =& -2_0^{\tau }\delta V_i\delta V_j{\left(\frac{\partial I_i}{\partial V_j}\right)}_{\mathrm{eq}}\mathrm{d}t,\hfill \end{array}$$

(40)

where we used Onsager’s reciprocity relation ${({\partial }_{V_j}{I}_i)}_{\mathrm{eq}}={({\partial }_{V_i}{I}_j)}_{\mathrm{eq}}$ [42,49]. The quantities ${({\partial }_{V_j}{I}_i)}_{\mathrm{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_{\Omega }=1$, $n_{V_i}=1$, and $n_{T_k}=0$, Equation (25) results in:

$$\begin{array}{ccc}\hfill \sum _r{0}^{\tau }\Omega \delta V_i{\left(\frac{{\partial }^2{J}_r}{\partial \Omega \partial V_i}\right)}_{\mathrm{eq}}\mathrm{d}t{+}_0^{\tau }\Omega {\left(\frac{\partial I_i}{\partial \Omega }\right)}_{\mathrm{eq}}\delta V_i\mathrm{d}t& =& {-}_0^{\tau }\delta V_i{\left(\frac{\partial F_{\Omega }}{\partial V_i}\right)}_{\mathrm{eq}}\Omega \mathrm{d}t.\hfill \end{array}$$

(41)

Now, using Onsager’s reciprocity relation ${(-{\partial }_{V_i}{F}_{\Omega })}_{\mathrm{eq}}={({\partial }_{\Omega }{I}_i)}_{\mathrm{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 $\left|\mathit{\alpha }\right|=2$ and beyond, but the number of relations and cases grows very fast with $\left|\mathit{\alpha }\right|$, and each relation may have its own physical interpretation. The study of higher-order terms in $\left|\mathit{\alpha }\right|$ 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 $\Omega$. 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 $\Omega$ 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 $\Omega$ expansion provided by the theory. The zeroth order terms in $\Omega$ reads:

$${\mathit{Q}}_J^{(0)}·\mathbf{1}+{\mathit{Q}}_I^{(0)}·\delta \mathit{V}=0.$$

(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:

$${\mathit{Q}}_J^{(1)}·\mathbf{1}+{\mathit{Q}}_I^{(1)}·\delta \mathit{V}=-{\mathcal{W}}_F^{(0)},$$

(43)

and can be understood as a generalization of the motor–pump relation, ${\mathit{Q}}_I^{(1)}·\delta \mathit{V}\approx -{\mathcal{W}}_F^{(0)}$, 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 $\delta \mathit{V}=0$ and $\delta \mathit{T}=0$, Equation (38) implies ${\mathcal{W}}_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 $\Omega$ gives:

$${\mathit{Q}}_J^{(2)}·\mathbf{1}+{\mathit{Q}}_I^{(2)}·\delta \mathit{V}=-{\mathcal{W}}_F^{(1)},$$

(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, ${\mathcal{W}}_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 (25), we take $\left|\mathit{\alpha }\right|$ up to some truncation value ${\mathit{\alpha }}_{\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 }_{\left|\mathit{\alpha }\right|=0}^{{\mathit{\alpha }}_{\max }}{\mathcal{W}}_F^{\mathit{\alpha }(\Omega )}}{{\sum }_{\left|\mathit{\alpha }\right|=0}^{{\mathit{\alpha }}_{\max }}{\mathit{Q}}_I^{\mathit{\alpha }(\mathit{V})}·\delta \mathit{V}}.$$

(45)

To illustrate this, let us take the case of a local system coupled to left and right reservoirs at voltages $V_L$ and $V_R$, respectively. We assume the leads are at the same temperature, and we set $\delta V_L=-\delta V_R=\delta 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 $\delta V=V_L-V_R$ and $\Omega$ up to ${\mathit{\alpha }}_{\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},$$

(46)

where the superscripts indicate the order in the expansions in $\Omega$ and $\delta V$, respectively, and we defined $I=(I_L-I_R)/2$. As we already mentioned, the zeroth order contributions ${\mathcal{W}}_F^{(0,0)}$ and $Q_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 ${\mathit{\alpha }}_{\max }=2$ but only in terms of $\Omega$, 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},$$

(47)

where now ${\mathcal{W}}_F^{(0)}$ and $Q_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 (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,58,59,87,88,89,90,91,92,93,94] and quantum motors [42].

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

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathit{p}=\mathit{W}\mathit{p},$$

(48)

where the vector $\mathit{p}=\left\{p_i(t)\right\}$ describes the dot’s occupation probabilities and $\mathit{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 $\dot{\mathit{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 $\mathit{p}={\sum }_{k\ge 0}{\mathit{p}}^{(k)}$, with ${\mathit{p}}^{(k)}$ of order ${(\Omega /\Gamma )}^k$. Here, $\Omega$ and $\Gamma$ 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 $\Omega <\Gamma$. The above master equation, in turn, takes the following hierarchical structure [12,78]:

$$\mathit{W}{\mathit{p}}^{(0)}=\mathbf{0},\mathrm{and}\mathit{W}{\mathit{p}}^{(k)}=\frac{\mathrm{d}}{\mathrm{d}t}{\mathit{p}}^{(k-1)}.$$

(49)

In the first equation, ${\mathit{p}}^{(0)}(\mathit{X})$ represents the steady-state solution the electronic system arrives at when the mechanical system is “frozen” at the point $\mathit{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 ${\mathit{p}}^{(k)}(\mathit{X},\dot{\mathit{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., $\mathit{W}=\mathit{W}(\mathit{X})$. The normalization condition on the dot’s density matrix implies ${\mathit{e}}^{\mathrm{T}}{\mathit{p}}^{(k)}={\delta }_{k0}$, with ${\mathit{e}}^{\mathrm{T}}$ the trace over the dot’s Hilbert space. This allows for the definition of an invertible pseudo-kernel $\tilde{\mathit{W}}$, whose matrix elements are ${\tilde{W}}_{ij}=W_{ij}-W_{ii}$, such that we can write:

$${\mathit{p}}^{(k)}={\left[{\tilde{\mathit{W}}}^{-1}\frac{\mathrm{d}}{\mathrm{d}t}\right]}^k{\mathit{p}}^{(0)}.$$

(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 ${\mathit{W}}^R$. This implies the same expansion as before, and it reads:

$$R^{(k)}={\mathit{e}}^{\mathrm{T}}{\mathit{W}}^R{\mathit{p}}^{(k)}.$$

(51)

The observables we are going to address here are the charge and heat currents $I_r$ and $J_r$ 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 $n_i$ and energy $E_i$ associated with the dot state $|i\rangle$ [87,95]:

$$W_{ij}^{I_r}=-e(n_i-n_j){\left[W_r\right]}_{ij},W_{ij}^{J_r}=-\left[E_i-E_j-{\mu }_r(n_i-n_j)\right]{\left[W_r\right]}_{ij}.$$

(52)

The CIF was derived in Section 2 under the assumption that the mechanical part of the system, characterized by coordinates $\mathit{X}$ and associated momenta $\mathit{P}$, only interacts with local parameters of the quantum dot through their many-body eigenenergies, cf. Equation (2). This implies that the $\nu$-component of the CIF operator, defined as ${\widehat{F}}_{\nu }=-\partial {\widehat{H}}_{\mathrm{s}}/\partial X_{\nu }$, is local in the quantum dot basis. For this local operator, then, we can simply define a diagonal kernel of the form:

$$W_{ij}^{F_{\nu }}=-\frac{\partial E_i}{\partial X_{\nu }}{\delta }_{ij},$$

(53)

such that Equation (51) gives the k order in the $\Omega$ expansion for any observable. Importantly, when adding up all contributions from the leads, the above kernel definitions, together with the sum rule ${\sum }_i{W}_{ij}=0$, lead to the following conservation rule for all orders in $\Omega$ [87]:

$$\sum _r{J}_r^{(k)}=-\frac{\mathrm{d}}{\mathrm{d}t}{\langle {\widehat{H}}_{\mathrm{s}}\rangle }^{(k-1)}-{\mathit{F}}^{(k-1)}·\dot{\mathit{X}}-\sum _r\frac{{\mu }_r}{e}{I}_r^{(k)}.$$

(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 $\tau$, and taking the time integral of the above equation, we recover the frequency expansion of Equation (25):

$$Q_J^{(k)}+Q_I^{(k)}\delta V=-{\mathcal{W}}_F^{(k-1)},$$

(55)

where the bias voltage $\delta V$ is defined through ${\mu }_L=-{\mu }_R=e\delta V/2$, $Q_J=Q_{J_L}+Q_{J_R}$ and $Q_I=(Q_{I_L}-Q_{I_R})/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],

$${\widehat{H}}_{\mathrm{s}}=\sum _{\ell =L,R}{ϵ}_{\ell }{\widehat{n}}_{\ell }+U{\widehat{n}}_L{\widehat{n}}_R+\frac{U^{\prime }}{2}\sum _{\ell =L,R}{\widehat{n}}_{\ell }({\widehat{n}}_{\ell }-1)-\frac{t_c}{2}\left({\widehat{d}}_{L\sigma }^{†}{\widehat{d}}_{R\sigma }+\mathrm{h}.\mathrm{c}.\right),$$

(56)

where ${ϵ}_{\ell }$ is the on-site energy and ${\widehat{n}}_{\ell }={\sum }_{\sigma }{\widehat{d}}_{\ell \sigma }^{†}{\widehat{d}}_{\ell \sigma }$ the number operator in the ℓ-dot, U and $U^{\prime }$ the inter- and intra-dot charging energies, respectively, and $t_c$ the interdot hopping amplitude. We will work in the strong coupling regime $t_c\gg \Gamma$, 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^{\prime }\to \infty$, 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 $\mathit{p}={(p_0,p_{\mathrm{b}\uparrow },p_{\mathrm{b}\downarrow },p_{\mathrm{a}\uparrow },p_{\mathrm{a}\downarrow },p_{\uparrow \uparrow },p_{\uparrow \downarrow },p_{\downarrow \uparrow },p_{\downarrow \downarrow })}^{\mathrm{T}}$. These elements thus denote the probability for the DQD to be either empty ($p_0$), occupied by one electron in the $\ell =\mathrm{b}$ (or $\mathrm{a}$) orbital with spin $\sigma$ ($p_{\ell \sigma }$), or by two electrons ($p_{\sigma {\sigma }^{\prime }}$), one of them in the left dot and with spin $\sigma$ and the other electron in the right dot and with spin ${\sigma }^{\prime }$. The many-body eigenenergies are therefore $E_0=0$ for the empty DQD,

$$E_{\mathrm{b}/\mathrm{a}}=\frac{{ϵ}_L+{ϵ}_R}{2}\mp \sqrt{{\left(\frac{{ϵ}_L-{ϵ}_R}{2}\right)}^2+{\left(\frac{t_c}{2}\right)}^2}$$

(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 $\theta$ describing the orientation of the rotor axis. We assume the following dependence through the on-site energies:

$${ϵ}_L(\theta )={\overline{ϵ}}_L+\lambda x_L={\overline{ϵ}}_L+{\delta }_{ϵ}cos(\theta ),\mathrm{and}{ϵ}_R(\theta )={\overline{ϵ}}_R+\lambda x_R={\overline{ϵ}}_R+{\delta }_{ϵ}sin(\theta ),$$

(58)

which defines a circular trajectory in energy space of radius ${\delta }_{ϵ}=\lambda r_0$ ($r_0$ measures the actual radius of the rotor) around the working point $({\overline{ϵ}}_L,{\overline{ϵ}}_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 $\dot{\theta }$, where: $\mathcal{I}=m_{\mathrm{eff}}{r}_0^2$ is the rotor’s moment of inertia, $\mathcal{F}={\sum }_i(-{\partial }_{\theta }{E}_i)p_i=r_0\mathit{F}·\widehat{\theta }$ is the current-induced torque, ${\mathcal{F}}_{\mathrm{ext}}$ is the torque produced by the external force, which is assumed to be constant along the whole trajectory, and ${\xi }_{\theta }$ accounts for the force fluctuations. The stationary regime is reached once the rotor performs a periodic motion (characterized by a time period $\tau$) 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 $\Delta \dot{\theta }$ over one period is small as compared to the average value of $\dot{\theta }$ 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 3, when we take $\mathcal{I}$ sufficiently large, $\Delta \dot{\theta }\to 0$, and it is possible to approximate $\dot{\theta }$ 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 $\dot{\theta }$ is not a priori known and, with it, the direction of the trajectory for the line integrals defining ${\mathcal{W}}_F^{(k)}$. To become independent of this issue, we denote by s the sign of $\dot{\theta }$, and notice that ${\mathcal{F}}^{(k)}=f_{\theta }^{(k)}{\dot{\theta }}^k$, where $f_{\theta }^{(k)}={\partial }_{\dot{\theta }}^k{\mathcal{F}|}_{\dot{\theta }=0}/k!$ only depends parametrically on $\theta$, so we have:

$${\mathcal{W}}_F^{(k)}={\oint }_{C(s)}{\mathit{F}}^{(k)}·\mathrm{d}\mathit{X}={\int }_0^{2\pi s}{\mathcal{F}}^{(k)}\mathrm{d}\theta =s\left[{\int }_0^{2\pi }{f}_{\theta }^{(k)}\mathrm{d}\theta \right]{\dot{\theta }}^k=s{\mathcal{C}}_F^{(k)}{\dot{\theta }}^k,$$

(59)

where, importantly, the coefficients ${\mathcal{C}}_F^{(k)}$ are independent of s, so the ${\dot{\theta }}^k$ and the extra sign s accompanying ${\mathcal{C}}_F^{(k)}$ in the right-hand side of the equation give the correct sign for ${\mathcal{W}}_F^{(k)}$. 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}}_{\mathrm{ext}},$$

(60)

where we defined ${\mathcal{C}}_{\mathrm{ext}}=2\pi {\mathcal{F}}_{\mathrm{ext}}$. The stability of the final solution is inherited from Equation (9), and in our case with $\dot{\theta }$ constant, it is given by:

$$\sum _k{\mathcal{C}}_F^{(k)}k{\dot{\theta }}^{k-1}<0.$$

(61)

Notice the specific case where ${\mathcal{C}}_{\mathrm{ext}}={\mathcal{C}}_F^{(0)}$, 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 $\dot{\theta }=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 $\dot{\theta }$, such that the solution is:

$$\dot{\theta }=\frac{{\mathcal{C}}_{\mathrm{ext}}-{\mathcal{C}}_F^{(0)}}{{\mathcal{C}}_F^{(1)}},$$

(62)

together with the condition ${\mathcal{C}}_F^{(1)}<0$, in agreement 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 ${\mathcal{C}}_F^{(1)}$ is fixed to this order, the sign in $\dot{\theta }$ only depends on the relation between ${\mathcal{C}}_{\mathrm{ext}}$ and ${\mathcal{C}}_F^{(0)}$. This, together with the stationary condition ${\mathcal{W}}_{\mathrm{ext}}={\mathcal{W}}_F$, establishes the operation mode of the device, in the sense that with this information, we can deduce the sign of ${\mathcal{W}}_{\mathrm{ext}}=s{\mathcal{C}}_{\mathrm{ext}}$ and, hence, the direction of the energy flow. If, for example, ${\mathcal{W}}_{\mathrm{ext}}<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 ${\mathcal{W}}_{\mathrm{ext}}>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 ${\mathcal{C}}_{\mathrm{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 in ${\mathcal{C}}_{\mathrm{ext}}$ where we can expect some kind of pumping mechanism. Depending then on the specific type of pumping we have in mind, the ranges in ${\mathcal{C}}_{\mathrm{ext}}$ 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 ${\mathcal{C}}_{\mathrm{ext}}$ 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 ${\mathcal{C}}_F^{(2)}$ in Equation (60), where we obtain:

$${\dot{\theta }}_{\pm }=-\frac{{\mathcal{C}}_F^{(1)}}{2{\mathcal{C}}_F^{(2)}}\pm \sqrt{{\left(\frac{{\mathcal{C}}_F^{(1)}}{2{\mathcal{C}}_F^{(2)}}\right)}^2+\frac{{\mathcal{C}}_{\mathrm{ext}}-{\mathcal{C}}_F^{(0)}}{{\mathcal{C}}_F^{(2)}}},$$

(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}}_{\mathrm{ext}}}{{\mathcal{C}}_F^{(2)}},\pm {\mathcal{C}}_F^{(2)}<0.$$

(64)

The first inequality ensures a positive argument in the square root of Equation (63), so $\dot{\theta }$ 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 get the stable solution: if ${\mathcal{C}}_F^{(2)}>0$, then ${\dot{\theta }}_{-}$, and if ${\mathcal{C}}_F^{(2)}<0$, then ${\dot{\theta }}_{+}$. Importantly, far from equilibrium (e.g., $\delta V\gg k_{\mathrm{B}}T$), it is possible to arrive at the odd situation where ${\mathcal{C}}_F^{(1)}>0$, and the “dissipated” energy ${\mathcal{W}}_F^{(1)}=s{\mathcal{C}}_F^{(1)}\dot{\theta }={\mathcal{C}}_F^{(1)}\left|\dot{\theta }\right|$ is positive. This, however, does not prevent the rotor from reaching the stationary regime (to this order in $\dot{\theta }$), 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 $\dot{\theta }$ now depends on the relation between the first term and the square root, together with the sign of ${\mathcal{C}}_F^{(2)}$, so this analysis is not that simple as in the linear case. However, once we have the ${\mathcal{C}}_F$-coefficients, we can always determine if the rotor is able to reach the stationary regime and infer whether ${\mathcal{W}}_{\mathrm{ext}}$ is positive or negative and, with it, the operation mode of the device.

For higher orders in $\dot{\theta }$, 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 $\dot{\theta }$. Interestingly, by including higher order terms in Equation (60), it could happen that for a fixed choice of the parameters (${\mathcal{F}}_{\mathrm{ext}}$, $\delta V$, $\delta T$, etc.), the system presents more than one stable solution, and even with different signs in ${\mathcal{W}}_{\mathrm{ext}}$, so the initial condition on $\dot{\theta }$ 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 $\mathcal{F}$ up to third order in the $\dot{\theta }$ expansion, where we find two stable solutions ($\dot{\theta }\sim -7.4\times {10}^{-3}{k}_{\mathrm{B}}T$ and $6.7\times {10}^{-3}{k}_{\mathrm{B}}T$) and one unstable solution in between ($\dot{\theta }=2\times {10}^{-3}{k}_{\mathrm{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 $\dot{\theta }=0$, the rotor gets trapped in a local minimum.

In Figure 3, we show how the ${\mathcal{C}}_F$-coefficients contribute to the current-induced work as a function of the bias voltage for different values of the orbit radius ${\delta }_{ϵ}$. According to Equation (59), these coefficients can only be compared upon multiplication with the k-power of some frequency of reference ${\Omega }_{\mathrm{x}}$. We assess this frequency of reference by performing the following analysis: from Equation (50), we can estimate the maximum allowed value for $\dot{\theta }$ compatible with the frequency expansion. Since ${\tilde{\mathit{W}}}^{-1}\propto {\Gamma }^{-1}$, with $\Gamma ={\Gamma }_L+{\Gamma }_R$, and:

$$\frac{\mathrm{d}}{\mathrm{d}t}=\dot{\theta }\frac{\partial }{\partial \theta }=\dot{\theta }\sum _i\frac{\partial {ϵ}_i}{\partial \theta }\frac{\partial }{\partial {ϵ}_i},$$

(65)

we obtain:

$${\mathit{p}}^{(k)}={\left[{\tilde{\mathit{W}}}^{-1}\dot{\theta }\sum _i\frac{\partial {ϵ}_i}{\partial \theta }\frac{\partial }{\partial {ϵ}_i}\right]}^k{\mathit{p}}^{(0)}\propto {\left[\frac{\Omega }{\Gamma }\frac{{\delta }_{ϵ}}{k_{\mathrm{B}}T}\right]}^k,$$

(66)

where we use the fact that ${\partial }_{\theta }{ϵ}_i\propto {\delta }_{ϵ}$ and the energy derivative of both ${\mathit{p}}^{(0)}$ and ${\tilde{\mathit{W}}}^{-1}$ are proportional to $1/k_{\mathrm{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_{\mathrm{B}}T}{{\delta }_{ϵ}},$$

(67)

from which we can estimate the maximum allowed frequency ${\Omega }_{\max }$ as:

$${\Omega }_{\max }=\frac{k_{\mathrm{B}}T}{{\delta }_{ϵ}}\Gamma .$$

(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 ${\Omega }_{\mathrm{x}}={\Omega }_{\max }/2$. As we can see, for a long range of the bias voltage ($V\lesssim 10k_{\mathrm{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 ${\mathcal{C}}_{\mathrm{ext}}\sim {\mathcal{C}}_F^{(0)}$. For $V\gtrsim 10k_{\mathrm{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 (${\mathcal{C}}_F^{(4)}$) has a negligible impact on the third-order solutions.

To clarify the role of ${\mathcal{C}}_F$-coefficients in the operation regime of the device, let us take for example ${\delta }_{ϵ}=100k_{\mathrm{B}}T$ (see Figure 3d) and $\delta V\sim 2k_{\mathrm{B}}T$. As the sign in ${\mathcal{C}}_F^{(1)}$ is negative, the linear-order solution given by Equation (62) is stable. If we start with ${\mathcal{C}}_{\mathrm{ext}}=0$, then $\dot{\theta }$ is positive, so for $0<{\mathcal{C}}_{\mathrm{ext}}<{\mathcal{C}}_F^{(0)}$, the device acts as a motor (see Figure 2a). When ${\mathcal{C}}_{\mathrm{ext}}>{\mathcal{C}}_F^{(0)}>0$, $\dot{\theta }$ becomes negative, and hence, the device operates as a pump (see Figure 2b). Additionally, when ${\mathcal{C}}_{\mathrm{ext}}<0$, the angular velocity is still positive, but as we changed the sign in ${\mathcal{C}}_{\mathrm{ext}}$, we have that ${\mathcal{W}}_{\mathrm{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 ${\mathcal{W}}_{\mathrm{ext}}=s{\mathcal{C}}_{\mathrm{ext}}$ changes with s and the sign of ${\mathcal{C}}_{\mathrm{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:

$${\mathcal{C}}_F^{(k)}(-\delta V)={(-1)}^{k+1}{\mathcal{C}}_F^{(k)}(\delta V),$$

(69)

since ${\Omega }_{\mathrm{x}}$ remains constant for a fixed orbit radius. For the chosen parameters ${\overline{ϵ}}_L={\overline{ϵ}}_R$, ${\Gamma }_L={\Gamma }_R$, and $\delta T=0$, the transformation $\delta V\to -\delta V$ can be regarded as the inversion operation $(L,R)\to (R,L)$, which changes the sign of $\dot{\theta }$, i.e., $\dot{\theta }(-\delta V)=-\dot{\theta }(\delta V)$. In this sense, if there is some finite temperature gradient, this operation should also involve the sign inversion of $\delta 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 $\dot{\theta }(-\delta V)=-\dot{\theta }(\delta V)$ if we also change the sign of the external force, such that ${\mathcal{C}}_{\mathrm{ext}}(-\delta V)=-{\mathcal{C}}_{\mathrm{ext}}(\delta 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 ${\overline{ϵ}}_L\ne {\overline{ϵ}}_R$ or ${\Gamma }_L\ne {\Gamma }_R$, such that the change in the sign of $\delta V$ (or $\delta 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 {\dot{\theta }}^k}\right|}_{\dot{\theta }=0}\frac{\mathrm{d}\theta }{k!},{\mathcal{C}}_J^{(k)}={\int }_0^{2\pi }{\left.\frac{\partial {}^k{J}_{\theta }}{\partial {\dot{\theta }}^k}\right|}_{\dot{\theta }=0}\frac{\mathrm{d}\theta }{k!},$$

(70)

such that $Q_I^{(k)}=s{\mathcal{C}}_I^{(k)}{\dot{\theta }}^{k-1}$, and $Q_J^{(k)}=s{\mathcal{C}}_J^{(k)}{\dot{\theta }}^{k-1}$. These contributions can be evaluated independently from the ${\mathcal{C}}_F$-coefficients, and in Figure 3, these are shown in circles, which gives numerical agreement for the energy conservation principle:

$${\mathcal{C}}_J^{(k)}+{\mathcal{C}}_I^{(k)}\delta V=-{\mathcal{C}}_F^{(k-1)},$$

(71)

in the considered orders of $\dot{\theta }$. For the considered example in Figure 3, the current coefficients ${\mathcal{C}}_I$ and ${\mathcal{C}}_J$ also show (independently) a definite parity with respect to the bias voltage. In fact, as Equation (71) suggests, the heat current coefficients ${\mathcal{C}}_J^{(k)}$ present the same parity as ${\mathcal{C}}_F^{(k-1)}$, while the charge current coefficients ${\mathcal{C}}_I^{(k)}$ need to have the opposite parity since they are multiplied by $\delta 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, determines 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 $\eta =$ (output power)/(input power), is always equal to one. 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 [87] 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 by ${\mathcal{W}}_{\mathrm{ext}}/\tau$, under the condition ${\mathcal{W}}_{\mathrm{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 $-{\mathit{Q}}_I·\delta \mathit{V}/\tau$, and hence, the efficiency of this type of motor is:

$$\eta =-\frac{{\mathcal{W}}_{\mathrm{ext}}}{{\mathit{Q}}_I·\delta \mathit{V}}=1+\frac{Q_J}{Q_I\delta V}.$$

(72)

Equation (33) establishes that the maximum efficiency for this device is $\eta \le 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 $\dot{\theta }$, we consider Equation (60) up to first order in ${\mathcal{C}}_F^{(k)}$, then, as discussed in Section 6, 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:

$$\eta =-\frac{{\mathcal{C}}_{\mathrm{ext}}/\delta V}{{\mathcal{C}}_I^{(0)}{\dot{\theta }}^{-1}+{\mathcal{C}}_I^{(1)}+{\mathcal{C}}_I^{(2)}\dot{\theta }}.$$

(73)

In the limits of the motor’s operation regime, given by ${\mathcal{C}}_{\mathrm{ext}}=0$ and ${\mathcal{C}}_{\mathrm{ext}}={\mathcal{C}}_F^{(0)}$, it is easy to see that the efficiency goes to zero, since for ${\mathcal{C}}_{\mathrm{ext}}=0$, the numerator in the above expression is zero, and for ${\mathcal{C}}_{\mathrm{ext}}={\mathcal{C}}_F^{(0)}$ the rotor’s velocity goes to zero, so the denominator grows to infinity due to the contribution ${\mathcal{C}}_I^{(0)}/\dot{\theta }$ from the leakage current. The same happens if we consistently include higher order terms in this expression. For example, if we use Equation (63) for the rotor’s angular velocity, then we should add ${\mathcal{C}}_I^{(3)}{\dot{\theta }}^2$ in the denominator.

Away from this region, we enter in the “pumping domain” characterized by a charge current, which opposes the “natural” direction dictated by the bias voltage and ${\mathcal{W}}_{\mathrm{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:

$$\eta =-\frac{{\mathcal{C}}_I^{(0)}{\dot{\theta }}^{-1}+{\mathcal{C}}_I^{(1)}+{\mathcal{C}}_I^{(2)}\dot{\theta }}{{\mathcal{C}}_{\mathrm{ext}}/\delta V}.$$

(74)

In this regime there is, however, an additional condition to be fulfilled, which is $Q_I\delta V>0$. Regarding the different orders in $Q_I$, it usually happens that close to the transition point ${\mathcal{C}}_{\mathrm{ext}}={\mathcal{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 to reverse the direction of the charge current. Going away from this region, the angular velocity acquires some finite value, reducing the zeroth-order contribution ${\mathcal{C}}_I^{(0)}/\dot{\theta }$ to the point where it is equal to the higher order contributions ${\mathcal{C}}_I^{(1)}+{\mathcal{C}}_I^{(2)}\dot{\theta }$, thus marking the activation point of the charge pump. In Figure 4a, we show the efficiency of the device as a function of the external force for different orbit sizes and fixed bias $\delta V=2k_{\mathrm{B}}T$. In all cases, the device starts from ${\mathcal{C}}_{\mathrm{ext}}=0$ as 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 at ${\mathcal{C}}_{\mathrm{ext}}={\mathcal{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 charge $Q_I^{(1)}$ is smaller than its quantized limit, and thereby, it takes a larger value of ${\mathcal{C}}_{\mathrm{ext}}$ to compensate the amount of pumped leakage current $Q_I^{(0)}$ in a cycle.

Up to this point, we have discussed only the effect of the lowest order terms of the expansion in $\dot{\theta }$, given by Equations (62), (73), and (74). For the smallest orbit, in addition, we show in the dotted red line the next-order efficiency, obtained from Equation (63) and adding the $Q_I^{(3)}$ term to Equations (73) and (74). We can see that for the motor regime, there are no significant changes, but for the pump regime, important differences appear. Firstly, in the ${\mathcal{C}}_{\mathrm{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 Figure 5a,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 ${\mathcal{C}}_{\mathrm{ext}}<0$ region. There, the efficiency is nonzero for ${\mathcal{C}}_{\mathrm{ext}}/{\mathcal{C}}_F^{(0)}<-8$ (see the inset in Figure 4a), 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 ${\mathcal{F}}_{\mathrm{ext}}$ in Equation (7) and the chosen parameters, in this region, the sign of $\dot{\theta }$ remains the same as that when ${\mathcal{F}}_{\mathrm{ext}}=0$. There, the zeroth- and first-order contributions to the charge current flow in the same direction. In the analyzed case, again the third-order term is the one that reverses the direction of the total charge current; see Figure 5b. It is important to mention that the purpose of the present discussion is only to highlight deviations from the linear solution of Equation (62), not to analyze the convergence of the total pumped current. For the larger values of ${\delta }_{ϵ}$ used in Figure 4a, we do not show the next-order corrections, as they are negligible in the shown range of ${\mathcal{C}}_{\mathrm{ext}}$.

An analysis similar to the above one can be carried out for the device driven by a temperature gradient $\delta T$, defined through $T_L=T+\delta T/2$ and $T_R=T-\delta T/2$ and no bias voltage applied, such that for $\delta T>0$, we have $T_{\mathrm{hot}}=T_L$ and $T_{\mathrm{cold}}=T_R$. When ${\mathcal{W}}_{\mathrm{ext}}>0$, we 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_{\mathrm{hot}}}/\tau$, and Equation (35) implies:

$$\eta =-\frac{{\mathcal{W}}_{\mathrm{ext}}}{Q_{J_{\mathrm{hot}}}}=1+\frac{Q_{J_{\mathrm{cold}}}}{Q_{J_{\mathrm{hot}}}}\le 1-\frac{T_{\mathrm{cold}}}{T_{\mathrm{hot}}}={\eta }_{\mathrm{carnot}}.$$

(75)

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

$$\eta =-\frac{{\mathcal{C}}_{\mathrm{ext}}}{{\mathcal{C}}_{J_{\mathrm{hot}}}^{(0)}{\dot{\theta }}^{-1}+{\mathcal{C}}_{J_{\mathrm{hot}}}^{(1)}+{\mathcal{C}}_{J_{\mathrm{hot}}}^{(2)}\dot{\theta }},$$

(76)

is defined in the same range for ${\mathcal{C}}_{\mathrm{ext}}$ as in the electric motor. Regarding Figure 4b, the engine’s normalized efficiency $\eta /{\eta }_{\mathrm{carnot}}$ looks similar to that of the electric motor. Perhaps the only difference here is that for the smallest orbit ${\delta }_{ϵ}=10k_{\mathrm{B}}T$, the efficiency maximum is very low, such that it cannot 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 $\delta T$. Therefore, in addition to the ${\mathcal{W}}_{\mathrm{ext}}<0$ condition, the overall amount of pumped heat in the cold reservoir should be negative, i.e., $Q_{J_{\mathrm{cold}}}<0$. The efficiency of the refrigerator, or coefficient of performance (COP), is then given by:

$$\mathrm{COP}=\frac{Q_{J_{\mathrm{cold}}}}{{\mathcal{W}}_{\mathrm{ext}}}=\frac{Q_{J_{\mathrm{hot}}}}{Q_J}-1\le \frac{T_{\mathrm{cold}}}{T_{\mathrm{hot}}-T_{\mathrm{cold}}}={\mathrm{COP}}_{\mathrm{carnot}},$$

(77)

where we can consistently expand $Q_{J_{\mathrm{cold}}}$ in terms of $\dot{\theta }$. In Figure 4b, we show the lowest order contribution from Equation (62), as the next-order calculation does not change significantly the efficiencies in the considered regimes of the parameters. Again, we can observe in Figure 4b a gap region where the device is frustrated since the work delivered by the rotor is not enough 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 4a are very similar. Additionally, for the chosen value $\delta T=0.5T$ and small orbit radius, the device can only work as a heat engine, and the refrigerator cannot be activated even if the external force is large, as happens for ${\delta }_{ϵ}=20k_{\mathrm{B}}T$ (solid blue line). The reason for this relies on the competition between the different orders in the pumped heat $Q_{J_{\mathrm{cold}}}$: the reduction of the leakage pumped current, $Q_{J_{\mathrm{cold}}}^{(0)}$, by increasing the magnitude of $\dot{\theta }$, is accompanied by an increase of the second-order contribution, $Q_{J_{\mathrm{cold}}}^{(2)}$, such that the first-order term, $Q_{J_{\mathrm{cold}}}^{(1)}$, may not be able to compensate these two and the requirement $Q_{J_{\mathrm{cold}}}<0$ cannot be fulfilled; see Figure 6.

8. Summary

Throughout this manuscript, we revisited some fundamental aspects related to 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, in Section 3. 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 took part in the quantum transport problem. This enabled us to derive an order-by-order scheme for the energy conservation law, Equation (25). 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 7, 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 at which the steady-state velocity changes its sign and, with it, the direction of the energy flow. When considering a specific type of pumped current (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 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 confirm numerically 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 participated in the pumped heat of the cold reservoir, highlighting the importance of the order-by-order conservation laws.