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

$$
\hat{H} = \sum\_{r} \hat{H}\_{r} + \hat{H}\_{\sf s} + \sum\_{r} \hat{H}\_{\sf s}.\tag{15}
$$

The mechanical part of the local system, when treated classically in Equation (7), introduces an explicit time dependence into *H* ˆ s that, in turn, makes d *H*ˆ /d*t* = 0. According to the above Hamiltonian, the total internal energy of the electronic system, *U* = *H*ˆ , 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{M}\_{\hat{\boldsymbol{\beta}}} = J\_{\hat{\boldsymbol{\beta}}}^{\mathrm{E}} + \left\langle \frac{\partial \hat{H}\_{\hat{\boldsymbol{\beta}}}}{\partial t} \right\rangle, \quad \text{where} \quad J\_{\hat{\boldsymbol{\beta}}}^{\mathrm{E}} = \dot{\boldsymbol{\epsilon}} \left\langle \left[ \hat{H}\_{\boldsymbol{\prime}} \hat{H}\_{\hat{\boldsymbol{\beta}}} \right] \right\rangle \tag{16}$$

which is the mean value of the energy flux entering in the subsystem *β* = {*<sup>r</sup>*,(s,*<sup>r</sup>*), <sup>s</sup>}. The value of *∂tH*<sup>ˆ</sup> *β* is zero when we evaluate the energy flux in the reservoirs and the tunnel couplings, but equals − W˙ *F* when *β* is evaluated in the local system. Although, strictly speaking, W*F* is time independent when considering a closed trajectory, we used the symbol W˙ *F* to denote the power delivered by the CIF, i.e., W˙ *F* = − *∂tH*<sup>ˆ</sup> s = *F* · *X*˙ . The latter comes from the definition of the CIF given in Equation (7). Note that the energy fluxes fulfill the condition ∑*r*(*J<sup>E</sup> r* + *J<sup>E</sup>*s,*r*) + *JE* s = 0. Therefore, the variation of the total internal energy of the electrons yields the following conservation rule:

$$
\dot{\mathbf{U}} = \sum\_{r} \dot{\mathbf{U}}\_{r} + \sum\_{r} \dot{\mathbf{U}}\_{\mathbf{s},r} + \dot{\mathbf{U}}\_{\mathbf{s}} = \sum\_{r} \left( f\_{r}^{E} + f\_{\mathbf{s},r}^{E} \right) + f\_{\mathbf{s}}^{E} - \dot{\mathcal{W}}\_{\mathbf{F}} = -\dot{\mathcal{W}}\_{\mathbf{F}}.\tag{17}
$$

Conservation of the total number of particles implies a relation between the particle currents of the reservoirs *N*˙ *r* and that of the local system *N*˙ s:

$$
\dot{N}\_{\text{total}} = \frac{\text{d}}{\text{d}t} \left< \dot{N}\_{\text{total}} \right> = \sum\_{r} \frac{\text{d}}{\text{d}t} \left< \dot{N}\_{r} \right> + \frac{\text{d}}{\text{d}t} \left< \dot{N}\_{\text{s}} \right> = 0 \quad \Rightarrow \quad \sum\_{r} \dot{N}\_{r} = -\dot{N}\_{\text{s}}.\tag{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 *τ* of the cyclic motion. Under this condition, the above-defined local quantities only depend periodically on time, i.e., *<sup>U</sup>*s(*<sup>t</sup>* + *τ*) = *<sup>U</sup>*s(*t*), *<sup>U</sup>*s,*r*(*<sup>t</sup>* + *τ*) = *<sup>U</sup>*s,*<sup>r</sup>*(*t*), and *N*˙ s(*<sup>t</sup>* + *τ*) = *N*˙ <sup>s</sup>(*t*). Therefore, the following quantities should evaluate to zero, i.e.,

$$
\int\_0^\tau \dot{\mathcal{U}}\_{\mathbf{s},r} \, \mathrm{d}t = 0, \quad \int\_0^\tau \dot{\mathcal{U}}\_{\mathbf{s}} \, \mathrm{d}t = 0, \quad \mathrm{and} \quad \int\_0^\tau \dot{\mathcal{N}}\_{\mathbf{s}} \, \mathrm{d}t = 0,\tag{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, ∑*r τ* 0 *Ir* d*t* = 0, where the charge current of the *r* reservoir is defined as *Ir* = *eN*˙ *r*, with *e* > 0 being minus the electron's charge.

We will take the following definition for the heat current *Jr* in the reservoir *r*:

$$J\_r = J\_r^E - \mu\_r \dot{N}\_r.\tag{20}$$

In [43,50,83], the authors proposed a different definition for the heat current of the reservoirs, *Jr* = *JE r* + *J<sup>E</sup>*s,*r*/2 − *<sup>μ</sup>rN*˙ *r*. However, the inclusion or not of half the heat current of the coupling region, *J<sup>E</sup>*s,*r*/2, 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\_{lr} + Q\_{lr} \delta V\_{r} \right) = -\mathcal{W}\_{\mathbf{F}'} \quad \text{where} \quad Q\_{lr} = \int\_{0}^{\tau} l\_{r} \, \mathrm{d}t, \quad \text{and} \quad Q\_{l\_{r}} = \int\_{0}^{\tau} l\_{r} \, \mathrm{d}t,\tag{21}$$

where *δVr* = (*μr* − *μ*0)/*<sup>e</sup>* and *μ*0 is an arbitrary reference's potential. We, in addition, defined the quantities *QJr* and *QIr*, 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 *χ*, which parameterizes the closed mechanical trajectory. This, for example, can be accounted for by taking *χ* = Ω mod(*<sup>t</sup>*, *<sup>τ</sup>*), such that *χ*˙ = Ω = 2*π*/*<sup>τ</sup>*. In other words, one can always find a natural scale, Ω in the present case, for all mechanical velocities along the trajectory. We will take this natural scale as the expansion variable for the CIF. The current-induced work then reads W*F* = *τ*0 *F*Ω Ω d*t*, where *F*Ω = *F* · *∂χ<sup>X</sup>*. Of course, with this *χ*-parameterization, the velocity is constant by definition, but at the expense that now the calculation of *F*Ω requires the knowledge of the mechanical trajectory. For circular trajectories and the conditions discussed in Section 3, the coordinate *χ* would be |*θ*|. For other cases, it would not be that simple to find *χ*, and one should first solve the system's equation of motion.

In addition to the mechanical velocity Ω, the CIF and the currents can also be expanded in terms of the remaining nonequilibrium sources, corresponding to voltage and temperature deviations *δV* = (*<sup>δ</sup>V*1, *δV*2, ...)<sup>T</sup> and *δT* = (*<sup>δ</sup>T*1, *δT*2, ...)T, respectively, from their equilibrium values *<sup>V</sup>*eq and *<sup>T</sup>*eq. For an arbitrary observable *R*, this general expansion takes the form:

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

where we used the following multi-index notation: *α* = *<sup>n</sup>*Ω, *nV*1 , ..., *nT*1 , ...; *α*! = *<sup>n</sup>*Ω!*nV*1 !...*nT*1 !...; (<sup>Ω</sup>, *δV*, *δT*)*α* = Ω*<sup>n</sup>*Ω *δVnV*<sup>1</sup> 1 ...*δTnT*<sup>1</sup> 1 ...; and:

$$\Theta^{\mathfrak{a}} = \frac{\mathfrak{d}^{\left(n\Omega + n\nu\_1 + \dots n\tau\_1 + \dots\right)}}{\mathfrak{d}\Omega^{\mathfrak{u}\Omega}\,\mathfrak{d}V\_1^{\mathfrak{u}V\_1}\dots\mathfrak{d}T\_1^{\mathfrak{u}\Gamma\_1}\dots}}.\tag{23}$$

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

$$\mathcal{D}\mathcal{W}\_{\mathbb{F}}^{\mathfrak{a}} = \int\_{0}^{\tau} F\_{\Omega}^{\mathfrak{a}} \Omega \mathrm{d}t, \quad \mathcal{Q}\_{l\_{r}}^{\mathfrak{a}} = \int\_{0}^{\tau} J\_{r}^{\mathfrak{a}} \mathrm{d}t, \quad \text{and} \quad \mathcal{Q}\_{l\_{r}}^{\mathfrak{a}} = \int\_{0}^{\tau} I\_{r}^{\mathfrak{a}} \mathrm{d}t,\tag{24}$$

where *<sup>F</sup><sup>α</sup>*Ω, *Jαr* , and *Iαr* 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: <sup>W</sup>*<sup>α</sup>F* contains an additional Ω from the time-integral as compared with the pumped currents, while the term *Q<sup>α</sup>Ir δVr* is one order higher in *δVr* than the other two. Therefore, the following relation must hold for every order *α* of the expansion:

$$
\mathbb{Q}\_I^{\mathfrak{a}} \cdot \mathbf{1} + \mathbb{Q}\_I^{\mathfrak{a}(V)} \cdot \delta V = -\mathcal{W}\_F^{\mathfrak{a}(\Omega)}\,\,\,\,\,\,\tag{25}
$$

where **1** = (1, 1, ...)T, *QαJ* = (*Q<sup>α</sup>J*1 , *Q<sup>α</sup>J*2 , ...), *Q<sup>α</sup>*(*V*) *I* = (*Q<sup>α</sup>*(*<sup>V</sup>*1) *I*1 , *Q<sup>α</sup>*(*<sup>V</sup>*2) *I*2 , ...), and we used the shorthand *<sup>α</sup>*(Ω)=(*<sup>n</sup>*Ω − 1, *nV*1 , ..., *nT*1 , ...) and *<sup>α</sup>*(*Vr*)=(*<sup>n</sup>*Ω, *nV*1 , ..., *nVr* − 1, ..., *nT*1 , ...). 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 Ω, *δV*, and *δT*. Importantly, the above conservation rule can be extended to other nonequilibrium sources like, e.g., spin polarization in ferromagnetic leads, such that other types of currents could also be considered in the energy transfer process between the electronic and mechanical parts of the system.
