Current-induced forces in single-resonance systems

As it is nowadays clear, electric currents can induce mechanical forces in nanodevices. Such current-induced forces (CIFs) [1–6], sometimes dubbed 'electron wind forces', have attracted increasing attention in many areas of condensed matter physics, including nanoelectromechanics [7–13], electromigration [14, 15], molecular wires and nanojuctions [16–25], the development of current-driven nanodevices [26–31], and quantum thermodynamics [32–38]. CIFs also play an important role in the understanding and development of molecular Brownian motors [39–41]. Among the most interesting examples are those where the forces are nonconservative. These types of forces may result in phenomena such as cooling, heating, or amplification of mechanical motions [8, 11, 20, 22]. Furthermore, it is also the basis of exciting proposals such as adiabatic quantum motors [27]. The most simplified description of a quantum motor consists of a system connected to leads and where a current of quantum particles drives some mechanical degrees of freedom. This is the reverse of what happens in a quantum pump, where the driving Hamiltonian's parameters, which may arise from a mechanical or an external electrical manipulation, induces a current [1, 27, 29, 33, 35, 36, 42]. Understanding the rules that dictate the interplay between the electronic and the mechanical degrees of freedom is crucial for the design of optimal electromechanical devices, and it can even contribute to the emergence of devices with novel features.

Aside from the potential applications, the type of devices discussed above may result particularly attractive from the theoretical point of view. This is so as quantum motors, quantum pumps, and some nanoelectromechanical devices can be interpreted as macroscopic manifestations of quantum behavior. Thus, these systems may provide new insights into the transition between the classical and quantum worlds, both from a dynamical [2, 29] and a thermodynamical [35–38] perspective. Moreover, the effective Hamiltonian of open quantum systems turns out to be non-Hermitian, which adds an extra richness to the problem. For example, the non-Hermiticity makes the dynamics of electrons to be affected in a nontrivial way by the change of some parameters, especially close to the so-called quantum dynamical phase transitions (QDPTs) [43–52].

In this work, we perform an in-depth study of the CIFs and the electronic-friction tensor of devices where a single energy level is relevant and, hence, they can be described by one of the most common Hamiltonian models in quantum transport [53–55], see figure 1. The main motivation behind that is to find the general conditions that lead to the improvement of the performance of different forms of nanoelectromechanical devices and quantum machines. Moreover, since CIFs may lead to the mechanical failure of conducting devices such as nanowires, the present work may also contribute to a better understanding of the circumstances that increase the chances of such failures.

Zoom In Zoom Out Reset image size

The work is organized as follows. In section 2 we discuss the main theoretical aspects of the studied system. These include, a brief introduction to the QDPTs and CIFs relevant to our problem, sections 2.1 and 2.2 respectively. In section 3 we present the results of this work and in section 4 we highlight the main findings and discuss their importance.

To avoid overwhelming the readers with the derivations of the many analytic expressions, we put all the nonessential comments and mathematical details in the appendices from A to H. We also include in appendix I an example of how to apply our results to concrete physical systems.

As we will see, the understanding of CIFs and electronic friction requires the comprehension of the different quantum dynamical regimes of the electrons. These regimes are separated by QDPTs, abrupt change in the decay dynamics of particles. As we will also see, this is intimately related to the maximization of the work done by CIFs. For that reason, in the first part of this section, we briefly review the QDPTs that our model Hamiltonian may undergo, before going to the theoretical description of CIFs.

2.1. Quantum dynamical phase transitions

The intuitive idea of dynamical phase transitions is better understood by considering the classical problem of a damped harmonic oscillator in the absence of external forces. There, the oscillator presents two well defined dynamical regimes, damped and overdamped motions, which typically can be reached by moving a single parameter [43, 56]. In this case, there is an analytical discontinuity in the plot of some dynamical observables, like the oscillation frequency, versus the control parameter. Due to its similarity with thermodynamical phase transitions, the phenomenon is known as dynamical phase transitions [43, 56].

In quantum mechanics, the same kind of phenomenon appears. There, the movement of a single parameter of the Hamiltonian may produce abrupt changes (non-analytical) in the decay dynamics of quantum particles or in their spectra [44, 45, 47–49, 52, 57].

For time-independent Hamiltonians, QDPTs are usually analyzed through the poles of the retarded (or advanced) Green function ${G}^{r}\left(\varepsilon \right)$ (or ${G}^{a}\left(\varepsilon \right)$) [44],

$$\begin{equation}{\boldsymbol{G}}^{r/a}\left(\varepsilon \right)=\underset{\eta \to {0}^{+}}{\mathrm{lim}}{\left(\left(\varepsilon {\pm}\mathrm{i}\eta \right)\boldsymbol{I}-\boldsymbol{H}\right)}^{-1}.\end{equation} \tag{ 1 }$$

This is natural as all spectral and dynamical properties of the system can be obtained from the elements of G ^r . For instance, the local density of states (LDOS) at a site n of a tight-binging Hamiltonian is given by:

$$\begin{equation}\begin{matrix}\hfill {N}_{n}\left(\varepsilon \right)=\underset{\eta \to {0}^{+}}{\mathrm{lim}}-\enspace \frac{1}{\pi }\mathrm{I}\mathrm{m}\left\{{G}_{nn}^{r}\left(\varepsilon +\mathrm{i}\eta \right)\right\}\hfill \end{matrix}.\end{equation} \tag{ 2 }$$

On the other hand, the survival probability of a particle in the nth site, ${P}_{n}\left(t\right)={\left\vert \left\langle {\psi }_{n}\left(t\right)\right\vert \left.{\psi }_{n}\left(0\right)\right\rangle \right\vert }^{2}$, can be expressed as [44]

$$\begin{equation}\begin{matrix}\hfill {P}_{n}\left(t\right)={\left\vert {\Theta}\left(t\right)\underset{-\infty }{\overset{\infty }{\int }}\mathrm{d}\varepsilon \enspace {N}_{n}\left(\varepsilon \right){\text{e}}^{-\text{i}\frac{\varepsilon t}{\hslash }}\right\vert }^{2}\hfill \end{matrix},\end{equation} \tag{ 3 }$$

where Θ(t) is the Heaviside step function.

2.1.1. Tight-binding model

In this paper, we will consider the tight-binding system depicted in figure 1. It consists of a single energy level with energy E_d coupled, with hopping constants V_L and V_R , to two semi-infinite tight-binding chains, with site energy E₀ and hopping V₀. This Hamiltonian represents a minimum model for a quantum dot coupled to two reservoirs and is widely used in the context of quantum transport. The total Hamiltonian reads

$$\begin{align}\hfill \hat{H}& =\sum\limits _{n=-\infty }^{\infty }{E}^{\left(n\right)}\left\vert n\right\rangle \left\langle n\right\vert \hfill \\ \hfill & \quad -{V}^{\left(n\right)}\left(\left\vert n\right\rangle \left\langle n+1\right\vert +\left\vert n+1\right\rangle \left\langle n\right\vert \right),\hfill \end{align} \tag{ 4 }$$

where E⁽⁰⁾ = E_d and E⁽ⁿ⁾ = E₀ for n ≠ 0, and the hopping constants V⁽⁻¹⁾ = V_L , V⁽⁰⁾ = V_R , and V⁽ⁿ⁾ = V₀ for $n\ne \left\{-1,0\right\}$. Note that we have disregarded the spin degree of freedom, which implies that our model is only valid in the absence of magnetic fields. We will keep this assumption throughout the manuscript.

For convenience, we will define ${\hat{H}}_{S}$ as the Hamiltonian of the system which includes the first sites of each semi-infinite chain representing the leads. Then, the three by three matrix H _S contains the site energies E⁽⁻¹⁾ = E₀, E⁽⁰⁾ = E_d , and E⁽¹⁾ = E₀, as well the couplings V_L and V_R .

The matrix representation of the eigenvalue equation $\hat{H}\left\vert \psi \right\rangle =\varepsilon \left\vert \psi \right\rangle$, which in principle is of infinite dimension, can be reduced to an effective system of finite dimension by using the decimation technique [54, 58]. There, the sites connected to leads are corrected by a self-energy associated to the G ^r (ɛ)

$$\begin{equation}{{\Sigma}}_{0}^{r}\left(\varepsilon \right)={\Delta}\left(\varepsilon \right)-\mathrm{i}{\Gamma}\left(\varepsilon \right).\end{equation} \tag{ 5 }$$

The real (Δ) and imaginary (Γ) parts of Σ are

$$\begin{equation}\begin{matrix}\hfill {\Delta}\left(\varepsilon \right)\hfill & \hfill =\hfill & \hfill \begin{cases}\hfill \frac{\varepsilon -{E}_{0}}{2}-\sqrt{{\left(\frac{\varepsilon -{E}_{0}}{2}\right)}^{2}-{V}_{0}^{2}}\quad \hfill & \hfill \left(\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}\enspace \mathrm{A}\right)\hfill \\ \hfill \frac{\varepsilon -{E}_{0}}{2}\quad \hfill & \hfill \left(\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}\enspace \mathrm{B}\right)\hfill \\ \hfill \frac{\varepsilon -{E}_{0}}{2}+\sqrt{{\left(\frac{\varepsilon -{E}_{0}}{2}\right)}^{2}-{V}_{0}^{2}}\quad \hfill & \hfill \left(\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}\enspace \mathrm{C}\right)\hfill \end{cases}\hfill \end{matrix}\;\;,\end{equation} \tag{ 6 }$$

$$\begin{equation}\begin{matrix}\hfill {\Gamma}\left(\varepsilon \right)\hfill & \hfill =\hfill & \hfill \begin{cases}\hfill 0\quad \hfill & \hfill \left(\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}\enspace \mathrm{A}\right)\hfill \\ \hfill \sqrt{{V}_{0}^{2}-{\left(\frac{\varepsilon -{E}_{0}}{2}\right)}^{2}}\quad \hfill & \hfill \left(\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}\enspace \mathrm{B}\right)\hfill \\ \hfill 0\quad \hfill & \hfill \left(\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}\enspace \mathrm{C}\right)\hfill \end{cases}\hfill \end{matrix}\;\;,\end{equation} \tag{ 7 }$$

where ɛ − E₀ ⩾ 2V₀ for case A, ɛ − E₀ ∈ [−2V₀, 2V₀] for case B, and ɛ − E₀ ⩽ −2V₀ for case C. The region where Γ ≠ 0 (case B), the band, corresponds to the energy regions where excitations propagate freely inside the leads. Using the above self-energy we can define the effective Hamiltonian ${\hat{H}}_{\mathrm{e}\mathrm{ff}}^{r}\left(\varepsilon \right)={\hat{H}}_{S}+{\hat{{\Sigma}}}^{r}\left(\varepsilon \right)$, where ${\hat{{\Sigma}}}^{r}\left(\varepsilon \right)$ is given by ${\hat{{\Sigma}}}^{r}\left(\varepsilon \right)={{\Sigma}}_{0}^{r}\left(\varepsilon \right)\left(\left\vert -1\right\rangle \left\langle -1\right\vert +\left\vert 1\right\rangle \left\langle 1\right\vert \right)$. In matrix form ${\hat{H}}_{\mathrm{e}\mathrm{ff}}^{r}$ is

$$\begin{equation}\begin{matrix}\hfill {\boldsymbol{H}}_{\mathrm{e}\mathrm{ff}}^{r}\left(\varepsilon \right)\hfill & \hfill =\hfill & \hfill \left(\begin{matrix}\hfill {E}_{0}+{{\Sigma}}_{0}^{r}\left(\varepsilon \right)\hfill & \hfill -{V}_{L}\hfill & \hfill 0\hfill \\ \hfill -{V}_{L}\hfill & \hfill {E}_{d}-\mathrm{i}{{\Gamma}}_{\eta }\hfill & \hfill -{V}_{R}\hfill \\ \hfill 0\hfill & \hfill -{V}_{R}\hfill & \hfill {E}_{0}+{{\Sigma}}_{0}^{r}\left(\varepsilon \right)\hfill \end{matrix}\right)\hfill \end{matrix}.\end{equation} \tag{ 8 }$$

By means of equation (1), ${\boldsymbol{H}}_{\mathrm{e}\mathrm{ff}}^{r}\left(\varepsilon \right)$ can be used to calculate G ^r , while ${\boldsymbol{G}}^{a}={\left[{\boldsymbol{G}}^{r}\right]}^{{\dagger}}$. The resulting eigenvalue equation $\left(\varepsilon \boldsymbol{I}-{\boldsymbol{H}}_{\mathrm{e}\mathrm{ff}}^{r}\left(\varepsilon \right)\right)\left\vert \psi \right\rangle =0$ is now finite and exact. The price to be paid is that now the coefficients of ${\boldsymbol{H}}_{\mathrm{e}\mathrm{ff}}^{r}$ are energy-dependent. Notice that, because ${{\Sigma}}^{r}\in \mathbb{C}$, the effective Hamiltonian becomes non-Hermitian, and thus, its eigenvalues, or equivalently, the poles of the associated Green's function G ^r are not necessarily real. Importantly, these poles can change abruptly with the parameters of the Hamiltonian and, in consequence, also the LDOS [see panels from (c) to (e) in figure 2] and the dynamics of the system [see panels from (f) to (h) in figure 2].

Zoom In Zoom Out Reset image size

Finally, one last detail. Due to technical reasons, we perturbatively couple a third lead to the system, see Γ_η in equation (8). This reservoir is modeled by a self-energy in the wideband approximation Σ_η = −iΓ_η with Γ_η positive. The purpose of this extra reservoir is to add a small broadening to the LDOS and to guarantee an occupation to the system. Both effects are particularly important when the eigenenergies of the effective Hamiltonian lay outside the band or when V_L = V_R = 0. Despite this, in all calculations, we take the limit Γ_η → 0, and thus its effects on the poles structure is negligible.

2.1.2. Poles and dynamical phase transitions

The poles of the system correspond to the singular points of the Green's function. The specific tight-binding model described in the previous section has two poles $\left\{{\varepsilon }_{p,+},{\varepsilon }_{p,-}\right\}$ given by the solutions of the secular equation $\mathrm{det}\left(\varepsilon \boldsymbol{I}-{\boldsymbol{H}}_{\mathrm{e}\mathrm{ff}}\right)=0$. Only three variables are really necessary to describe such poles, E_d , V_L , and V_R . The effect of E₀ is only to shift the poles' position while V₀ gives the energy scale. For that reason, we define the dimensionless parameters

$$\begin{align}\hfill {{\epsilon}}_{p,{\pm}}=\frac{{\varepsilon }_{p,{\pm}}-{E}_{0}}{{V}_{0}},\qquad {{\epsilon}}_{d}=\frac{{E}_{d}-{E}_{0}}{{V}_{0}}\\ \hfill \mathrm{a}\mathrm{n}\mathrm{d}\quad {v}^{2}={v}_{L}^{2}+{v}_{R}^{2}.\end{align} \tag{ 9 }$$

where ${v}_{L}=\frac{{V}_{L}}{{V}_{0}}$ and ${v}_{R}=\frac{{V}_{R}}{{V}_{0}}$. The real and imaginary part of the dimensionless poles are, respectively [46, 47],

$$\begin{align}\hfill \mathrm{R}\mathrm{e}\left({{\epsilon}}_{p,{\pm}}\right)& =\begin{cases}{{\epsilon}}_{d}+\frac{{v}^{2}}{\left(1-{v}^{2}\right)}\left[\frac{{{\epsilon}}_{d}}{2}{\pm}\sqrt{{\left(\frac{{{\epsilon}}_{d}}{2}\right)}^{2}+{v}^{2}-1}\right]\quad \hfill \\ {{\epsilon}}_{d}+\frac{{v}^{2}}{\left(1-{v}^{2}\right)}\frac{{{\epsilon}}_{d}}{2}\quad \hfill \end{cases}\hfill \\ \hfill \mathrm{I}\mathrm{m}\left({{\epsilon}}_{p,{\pm}}\right)& =\begin{cases}0\quad \hfill \\ {\pm}\frac{{v}^{2}}{\left(1-{v}^{2}\right)}\sqrt{1-\left[{\left(\frac{{{\epsilon}}_{d}}{2}\right)}^{2}+{v}^{2}\right]}\hfill \end{cases}\hspace{-2.5pt},\left.\right.\hfill \end{align} \tag{ 10 }$$

for v² ≠ 1. The first case of $\mathrm{R}\mathrm{e}\left({{\epsilon}}_{p,{\pm}}\right)$ and $\mathrm{I}\mathrm{m}\left({{\epsilon}}_{p,{\pm}}\right)$ corresponds to the condition ${\left(\frac{{{\epsilon}}_{d}}{2}\right)}^{2}+{v}^{2}{\geqslant}1$, while the other one corresponds to the condition ${\left(\frac{{{\epsilon}}_{d}}{2}\right)}^{2}+{v}^{2}{\leqslant}1$.

In figures 2(a) and (b) we show the real and imaginary part of the poles, _p,±, as a function of the dimensionless parameter _d . In the lower panels of the figure, see panels (c)–(h), we show the typical spectral and dynamical characteristic of three types of poles defining the resonant state, the virtual state, and the localized state.

Resonant states. They correspond to poles with a nonzero imaginary part. See the green dot marked as (1) in figure 2(a). When the system presents such a pole, the LDOS is typically a Lorentzian peak centered at approximately Re(_p ) and whose width is given by Im(_p ). This description is particularly accurate when Re(_p ) is close to the center of the band. However, when Re(_p ) approaches one of the band edges the maximum of the LDOS is shifted toward the closest band edge and its shape deviated from a perfect Lorentzian function. After a short quadratic decay, the dynamics of such systems typically shows an exponential decay [59].

Localized states. They correspond to poles with Im(_p ) = 0 and where the Re(_p ) lays outside the band (|Re(_p )| > 2). For instance, see the green dot marked as (3) in figure 2(b). When the system presents such a pole, its LDOS typically shows a very narrow peak outside the band (a Dirac's delta in the limit Γ_η → 0). The survival probability may show a small decay at short times but then it remains constant with time.

Virtual states. They also correspond to poles with Im(_p ) = 0 and |Re(_p )| > 2, see the green dot marked as (2) in figure 2(b). However, they are poles of the nonphysical Riemann sheet of Σ^r , see references [44, 46] for more details. Due to this, there is not a peak in the LDOS at Re(_p ). However, the effect of these poles on the LDOS within the band is the same as that of localized states. There is an accumulation of states near the closest band edge, which grows as the pole approaches this band edge. See the examples shown in panels (d) and (e).

The three type of poles described above determine two QDPT namely the resonant-virtual and the virtual-localized transitions. The resonant-virtual QDPT is given by [46, 47]

$$\begin{equation}{\left(\frac{{{\epsilon}}_{d}}{2}\right)}^{2}+{v}^{2}=1,\end{equation} \tag{ 11 }$$

while the virtual-localized QDPT is given by [46, 47]

$$\begin{equation}{v}^{2}=2\mp {{\epsilon}}_{d}.\end{equation} \tag{ 12 }$$

A more detailed classification may lead to another QDPT as resonant states can be subdivided into two additional type of poles, see references [44, 46]. However, for the purpose of this work this subclassification is irrelevant and will be ignored.

One final remark is still in order. In more general scenarios of a quantum dot coupled to m identical semi-infinite chains with coupling V_i , the corresponding poles are still given by equation (10), but the quadratic coupling v² should be replaced by

$$\begin{equation}\begin{matrix}\hfill {v}^{2}={\left(\frac{{V}_{1}}{{V}_{0}}\right)}^{2}+{\left(\frac{{V}_{2}}{{V}_{0}}\right)}^{2}+\cdots +{\left(\frac{{V}_{m}}{{V}_{0}}\right)}^{2}\hfill \end{matrix}.\end{equation} \tag{ 13 }$$

2.2. Current-induced forces

When a current of quantum particles is coupled to a mesoscopic mechanical device, the classical dynamics of the mechanical part can be described by an effective Langevin equation [1, 36]:

$$\begin{equation}\begin{matrix}\hfill {M}_{\nu }{\ddot {X}}_{\nu }+\frac{\partial U}{\partial {X}_{\nu }}={\mathcal{F}}_{\nu }+{\xi }_{\nu }\hfill \end{matrix}.\end{equation} \tag{ 14 }$$

Here, on the left-hand side we have the terms corresponding to the classical degrees of freedom, where the coordinates $\left\{{X}_{\nu }\right\}$ are associated with masses M_ν , and $\frac{\partial U}{\partial {X}_{\nu }}$ accounts for any external forces. The terms on the right-hand side of equation (14) give the CIFs, or the forces arising from the interaction between the mechanical device and the quantum particles (electrons in our case). The term ξ_ν describes the quantum fluctuations of the CIFs while ${\mathcal{F}}_{\nu }$ is the mean value of the force operator

$$\begin{equation}\begin{matrix}\hfill {\mathcal{F}}_{\nu }=\left\langle -\frac{\partial \hat{H}}{\partial {X}_{\nu }}\right\rangle \hfill \end{matrix},\end{equation} \tag{ 15 }$$

where $\hat{H}$ is the electronic Hamiltonian.

2.2.1. CIFs in the Keldysh formalism

The mean value of the force operator can be calculated in terms of Green's functions by

$$\begin{equation}\begin{matrix}\hfill {\mathcal{F}}_{\nu }=\int \frac{\mathrm{d}\varepsilon }{2\pi \mathrm{i}}\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }{\boldsymbol{\mathcal{G}}}^{{< }}\right]\hfill \end{matrix},\end{equation} \tag{ 16 }$$

where ${\mathbf{\Lambda }}_{\nu }=-\frac{\partial \boldsymbol{H}}{\partial {X}_{\nu }}$ and ${\boldsymbol{\mathcal{G}}}^{{< }}$ is the lesser Green's function. This function, which for shortness we omitted its energy and time dependence (${\boldsymbol{\mathcal{G}}}^{{< }}\equiv {\boldsymbol{\mathcal{G}}}^{{< }}\left(\varepsilon ,\tau \right)$), is indeed the Wigner transform of what generalizes the density matrix and is called the particle propagator ${\mathcal{G}}^{{< }}\left(x,t;{x}^{\prime },{t}^{\prime }\right)=\left(\mathrm{i}/\hslash \right)\left\langle {\psi }^{{\dagger}}\left({x}^{\prime },{t}^{\prime }\right)\psi \left(x,t\right)\right\rangle$ [1, 53]. The closely related Green function ${\mathcal{G}}^{{ >}}\left(x,t;{x}^{\prime },{t}^{\prime }\right)=-\left(\mathrm{i}/\hslash \right)\left\langle \psi \left(x,t\right){\psi }^{{\dagger}}\left({x}^{\prime },{t}^{\prime }\right)\right\rangle$ is dubbed the hole propagator, since the order of the creation and annihilation operators are reversed. Both functions are directly linked to observables and kinetic properties, such as particle densities and currents [53, 60].

In complex scenarios like the one treated here, where the time-dependent Hamiltonians are slowly varying functions of some parameters, the function ${\boldsymbol{\mathcal{G}}}^{{< }}$ can be calculated by resorting to an adiabatic expansion in terms of the time-dependent parameters [53, 60]. There, ${\boldsymbol{\mathcal{G}}}^{{< }}$ is given, up to first order, by [1],

$$\begin{align}\hfill {\boldsymbol{\mathcal{G}}}^{{< }}& \simeq {\boldsymbol{G}}^{{< }}-\frac{\mathrm{i}\hslash }{2}\sum\limits _{\nu }\dot {{X}_{\nu }}\left[\left({\partial }_{\varepsilon }{\boldsymbol{G}}^{{< }}\right){\mathbf{\Lambda }}_{\nu }{G}^{a}-{G}^{r}{\mathbf{\Lambda }}_{\nu }\left({\partial }_{\varepsilon }{\boldsymbol{G}}^{{< }}\right)\right.\hfill \\ \hfill & \quad \left.+\left({\partial }_{\varepsilon }{G}^{r}\right){\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{{< }}-{\boldsymbol{G}}^{{< }}{\mathbf{\Lambda }}_{\nu }\left({\partial }_{\varepsilon }{G}^{a}\right)\right],\hfill \end{align} \tag{ 17 }$$

where G ^< is the frozen lesser Green's function, given by

$$\begin{equation}{\boldsymbol{G}}^{{< }}={\boldsymbol{G}}^{r}{\mathbf{\Sigma }}^{{< }}{\boldsymbol{G}}^{a}.\end{equation} \tag{ 18 }$$

Here, G ^r/a are given by equation (1) and Σ^< is the lesser self-energy, ${\mathbf{\Sigma }}^{{< }}=2\mathrm{i}\sum _{\alpha }{f}_{\alpha }{\mathbf{\Gamma }}_{\alpha }$. The latter contains the information of the reservoirs, where f_α is the Fermi–Dirac distribution function of reservoir α with chemical potential μ_α and temperature T_α . The function Γ_α is the imaginary part of the self-energy of the reservoir α, and describes the escape rate from the region of the contact in units of ℏ. Importantly, in deriving equation (17) it was assumed that the self-energies Σ^< and Σ^r/a are time-independent. However, as we included the first sites of the leads into the effective Hamiltonian, equation (8), a variation of the couplings V_L and V_R does not change the self-energies.

Then, under the adiabatic approximation, ${\mathcal{F}}_{\nu }$ is

$$\begin{equation}\begin{matrix}\hfill {\mathcal{F}}_{\nu }\simeq {F}_{\nu }-\sum\limits _{{\nu }^{\prime }}{\gamma }_{\nu {\nu }^{\prime }}{\dot {X}}_{{\nu }^{\prime }}\hfill \end{matrix},\end{equation} \tag{ 19 }$$

where the 'frozen' force or, from now on, simply the force, is

$$\begin{equation}\begin{matrix}\hfill {F}_{\nu }=\int \frac{\mathrm{d}\varepsilon }{2\pi \mathrm{i}}\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{{< }}\right]\hfill \end{matrix}.\end{equation} \tag{ 20 }$$

The terms γ_νν' are the components of the electronic-friction tensor γ . The diagonal elements of γ provide the mechanical dissipation while the nondiagonal components are Lorentz-like terms which allow energy transfer between modes. This tensor can be linearly decomposed into a symmetric and an antisymmetric contribution γ = γ ^s + γ ^a . The antisymmetric component is finite only for nonzero voltage biases or temperature gradients between the reservoirs. Thus, close to equilibrium conditions only the symmetric contribution of the electronic-friction tensor needs to be taken into account.

2.2.2. Expansion of the CIFs

In this work, we are interested only on the leading orders of the expansion of $\mathcal{F}$ in terms of the different nonequilibrium sources δμ, δT, and $\dot { \overrightarrow {X}}$,

$$\begin{align}\hfill {\mathcal{F}}_{\nu }& \approx {\left.{\mathcal{F}}_{\nu }\right\vert }_{\text{eq}}+\sum\limits _{{\nu }^{\prime }}{\left.\frac{\partial {\mathcal{F}}_{\nu }}{\partial {\dot {X}}_{{\nu }^{\prime }}}\right\vert }_{\text{eq}}{\dot {X}}_{{\nu }^{\prime }}\hfill \\ \hfill & \quad +\sum\limits _{\alpha }{\left.\frac{\partial {\mathcal{F}}_{\nu }}{\partial {\mu }_{\alpha }}\right\vert }_{\text{eq}}\delta {\mu }_{\alpha }+\sum\limits _{\alpha }{\left.\frac{\partial {\mathcal{F}}_{\nu }}{\partial {T}_{\alpha }}\right\vert }_{\text{eq}}\delta {T}_{\alpha }.\hfill \end{align} \tag{ 21 }$$

The equilibrium contribution of the force,

$$\begin{equation}{F}_{\nu }^{\text{eq}}={\left.{\mathcal{F}}_{\nu }\right\vert }_{\text{eq}},\end{equation} \tag{ 22 }$$

is conservative and thus can be written as the gradient of a potential, ${ \overrightarrow {F}}^{\text{eq}}=-\nabla {U}^{\text{eq}}$. As we show in appendix A, the equilibrium potential U^eq is

$$\begin{equation}\begin{matrix}\hfill {U}^{\text{eq}}=-\int \frac{\mathrm{d}\varepsilon }{2\pi \mathrm{i}}{f}_{0}\left(\varepsilon \right)\mathrm{ln}\left[\frac{\mathrm{det}\left({\boldsymbol{G}}^{a}\right)}{\mathrm{det}\left({\boldsymbol{G}}^{r}\right)}\right]\hfill \end{matrix},\end{equation} \tag{ 23 }$$

where f₀ is the equilibrium Fermi–Dirac distribution function given by μ₀ and T₀, the equilibrium chemical potential and temperature respectively.

Comparing equations (21) and (19), it is obvious that ${\left.\frac{\partial {\mathcal{F}}_{\nu }}{\partial {\dot {X}}_{{\nu }^{\prime }}}\right\vert }_{\text{eq}}={\gamma }_{\nu ,{\nu }^{\prime }}{\vert }_{\text{eq}}$. Therefore, only the equilibrium contribution of the symmetric component of γ is relevant for the present work, where ${\gamma }_{\nu {\nu }^{\prime }}^{s,\mathrm{e}\mathrm{q}}\equiv {\gamma }_{\nu {\nu }^{\prime }}^{\text{eq}}$ and

$$\begin{equation}\begin{matrix}\hfill {\gamma }_{\nu {\nu }^{\prime }}^{\text{eq}}=\int \frac{\hslash \enspace \mathrm{d}\varepsilon }{4\pi }\enspace \mathrm{Tr}\left[\left({\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}_{\text{eq}}^{{< }}{\mathbf{\Lambda }}_{{\nu }^{\prime }}+{\mathbf{\Lambda }}_{{\nu }^{\prime }}{\boldsymbol{G}}_{\text{eq}\;}^{{< }}{\mathbf{\Lambda }}_{\nu }\right){\partial }_{\varepsilon }{\boldsymbol{G}}_{\text{eq}}^{{ >}}\right]\hfill \end{matrix}.\end{equation} \tag{ 24 }$$

Here, ${\boldsymbol{G}}_{\text{eq}}^{\lessgtr }$ are the equilibrium greater (>) and lesser (<) Green's functions. The above equation can be rewritten, see appendix B, in a simpler form,

$$\begin{equation}{\gamma }_{\nu {\nu }^{\prime }}^{\text{eq}}=\int \frac{\hslash \enspace \mathrm{d}\varepsilon }{4\pi }\left(-{\partial }_{\varepsilon }{f}_{0}\right)\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right],\end{equation} \tag{ 25 }$$

where A is the spectral function, $\boldsymbol{A}=\mathrm{i}\left({\boldsymbol{G}}^{r}-{\boldsymbol{G}}^{a}\right)$. This equation shows that only states with energies close to the Fermi energy contribute to the equilibrium electronic-friction tensor. Those states are also responsible for the electric current [54] and the nonequilibrium contribution of the CIFs, as we will see afterward. Therefore, equation (25) highlights a fundamental physical connection between these quantities. In the low temperature limit, equation (25) becomes

$$\begin{equation}\underset{{k}_{\mathrm{B}}{T}_{0}\to {0}^{+}}{\mathrm{lim}}\enspace {\gamma }_{\nu {\nu }^{\prime }}^{\text{eq}}=\frac{\hslash }{4\pi }\enspace \mathrm{Tr}{\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right]}_{\varepsilon ={\mu }_{0}}.\end{equation} \tag{ 26 }$$

Stochastic forces, in equation (14), can be obtained by using the above equations and the fluctuation–dissipation relation [1], given by

$$\begin{equation}{D}_{\nu {\nu }^{\prime }}^{\text{eq}}=2{k}_{\mathrm{B}}{T}_{0}{\gamma }_{\nu {\nu }^{\prime }}^{\text{eq}},\end{equation} \tag{ 27 }$$

where k_B is the Boltzmann constant and it is assumed that stochastic forces are described by $\left\langle {\xi }_{\nu }\left(t\right){\xi }_{{\nu }^{\prime }}\left({t}^{\prime }\right)\right\rangle ={D}_{\nu {\nu }^{\prime }}^{\text{eq}}\delta \left(t-{t}^{\prime }\right)$.

At low voltage biases and temperature gradients, the nonequilibrium forces are given by

$$\begin{equation}{F}_{\nu }^{\text{ne}}=\sum\limits _{\alpha }{\left.\frac{\partial {\mathcal{F}}_{\nu }}{\partial {T}_{\alpha }}\right\vert }_{\text{eq}}\delta {T}_{\alpha }+{\left.\frac{\partial {\mathcal{F}}_{\nu }}{\partial {\mu }_{\alpha }}\right\vert }_{\text{eq}}\delta {\mu }_{\alpha },\end{equation} \tag{ 28 }$$

which in the low-temperature limit (see appendix C) gives

$$\begin{align}\hfill {F}_{\nu }^{\text{ne}}& =\sum\limits _{\alpha }\frac{\pi }{3}{\left({k}_{\mathrm{B}}{T}_{0}\right)}^{2}\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }{\left.\frac{\partial \left({\boldsymbol{G}}^{r}{{\Gamma}}_{\alpha }{\boldsymbol{G}}^{a}\right)}{\partial \varepsilon }\right\vert }_{\varepsilon ={\mu }_{0}}\right]\frac{\delta {T}_{\alpha }}{{T}_{0}}\hfill \\ \hfill & \quad +\frac{1}{\pi }\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }{\left({\boldsymbol{G}}^{r}{{\Gamma}}_{\alpha }{\boldsymbol{G}}^{a}\right)}_{\varepsilon ={\mu }_{0}}\right]\delta {\mu }_{\alpha }.\hfill \end{align} \tag{ 29 }$$

Here, δμ_α = μ_α − μ₀, δT_α = T_α − T₀, and the expression is only valid for δT_α ⩽ T₀.

2.2.3. Work

One of the main motivations of this work is to study the general conditions that maximize the performance of different nanoelectromechanical devices and quantum machines. In this respect, one of the central quantities to evaluate is the work W done by CIFs during cyclic motions. If one is dealing with systems where only two or three parameters are required to describe their motion, then it is useful to study the curl of the CIFs, defined as

$$\begin{equation}{\left(\nabla {\times} \overrightarrow {F}\right)}_{\rho }={\partial }_{\nu }{F}_{{\nu }^{\prime }}-{\partial }_{{\nu }^{\prime }}{F}_{\nu },\end{equation} \tag{ 30 }$$

where ρ, ν, and ν' are the indexes of the coordinates and ν < ν'. Now, the work can be obtained from

$$\begin{equation}\begin{matrix}\hfill W=\underset{S}{\iint }\left(\nabla {\times} \overrightarrow {F}\right){\cdot}\mathrm{d} \overrightarrow {S}\hfill \end{matrix},\end{equation} \tag{ 31 }$$

where the integration is done on the surface S enclosed by the closed curve C that describes the cyclic motion.

Now inserting equation (29) into equation (30), gives

$$\begin{align}\hfill {\left(\nabla {\times} \overrightarrow {F}\right)}_{\rho }& =\sum\limits _{\alpha }\frac{\delta {\mu }_{\alpha }}{\pi \mathrm{i}}\enspace \mathrm{Tr}{\left[{\left\{{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}{\boldsymbol{G}}^{r}{\mathbf{\Gamma }}_{\alpha }{\boldsymbol{G}}^{a}\right\}}_{a}\right]}_{\varepsilon ={\mu }_{0}}\hfill \\ \hfill & \quad +\frac{\pi }{3i}{\left({k}_{\mathrm{B}}{T}_{0}\right)}^{2}\frac{\delta {T}_{\alpha }}{{T}_{0}}\hfill \\ \hfill & \quad {\times}\frac{\partial }{\partial \varepsilon }\enspace \mathrm{Tr}{\left[{\left\{{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}{\boldsymbol{G}}^{r}{\mathbf{\Gamma }}_{\alpha }{\boldsymbol{G}}^{a}\right\}}_{a}\right]}_{\varepsilon ={\mu }_{0}},\hfill \end{align} \tag{ 32 }$$

where we used ${\boldsymbol{G}}^{r/a}{\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{r/a}={\partial }_{\nu }\left({\boldsymbol{G}}^{r/a}\right)$, ∂_ν' Λ_ν = ∂_ν Λ_ν', and ${\left\{{\mathbf{\Lambda }}_{\nu }..{\mathbf{\Lambda }}_{{\nu }^{\prime }}..\right\}}_{a}=\left({\mathbf{\Lambda }}_{\nu }..{\mathbf{\Lambda }}_{{\nu }^{\prime }}..-{\mathbf{\Lambda }}_{{\nu }^{\prime }}..{\mathbf{\Lambda }}_{\nu }..\right)$. Equation (32) gives the work done per unit area in the low temperatures limit and for small bias voltages or temperature gradients.

2.2.4. CIFs in our model Hamiltonian

Let us return to the tight-binding model described in section 2.1.1. For the analysis of the forces, especially for the equilibrium forces, it is crucial to have a definite occupation of the system, even for localized states or when V_L = V_R = 0. As discussed in section 2.1.2, when the system is at the localized state, the poles of G ^r , or the eigenenergies of H _eff, lay outside the band. At this condition, the imaginary part of the self-energies, Γ, is zero. The same happens for V_L = V_R = 0. This implies the nonphysical condition that localized states do not contribute to the forces, or even worse, that isolated systems do not have forces at all. To solve this problem, we added the third reservoir (η) to the model, assuming its occupation is given by the equilibrium Fermi function f₀ = (f_L + f_R )/2, and taking the limit Γ_η → 0 in all calculations.

Given the Hamiltonian of equation (8) and the definition of G ^r/a (equation (1)) and G ^< (equation (18)), we obtain a closed-form expression for the Green functions, (see appendix D). As can be seen from the expressions of the appendix, the third reservoir does not affect the calculations for energies within the band, provided Γ_η is small. However, it allows the calculation of equilibrium G ^< (and thus equilibrium forces) at all conditions. Then, all occupied states, independently if they are inside or outside the band, contribute to the equilibrium forces. One important point to emphasize is that the nonequilibrium contribution of G ^< (and thus ${ \overrightarrow {F}}^{\left(\mathrm{n}\mathrm{e}\right)}$) is always zero for states outside the band (localized states), or for V_L = V_R = 0. This will have important consequences when studying the maximization of the curl of the forces.

According to equation (16), the forces depend on how the elements of the electronic Hamiltonian are affected by the movement of the mechanical degrees of freedom, but this is system-dependent. Therefore, to gain generality, we take as 'mechanical' variables the same parameters of the Hamiltonian (E_d , V_L , V_R ), which leads to dimensionless forces F_ν , see appendix E. One can readily calculate the forces in terms of 'physical' variables q_i by

$$\begin{equation}{F}_{i}=\sum\limits _{\nu }{F}_{\nu }\frac{\partial {X}_{\nu }}{\partial {q}_{i}},\end{equation} \tag{ 33 }$$

where ${X}_{\nu }=\left\{{E}_{d},{V}_{L},{V}_{R}\right\}$. Notice that while forces may depend on the chosen variables, the work in equation (31) is independent of them, as long as the curve C remains the same, see appendix E. Therefore, using our atypical choice of 'mechanical' variables in equation (32) provides a useful tool to analyze in a general way the performance of nanoelectromechanical devices or quantum machines without having to specify the details of their physical implementation.

3.1. Equilibrium potential

Using equation (23) and the expression for G^r/a obtained in appendix D, we calculate the equilibrium potential

$$\begin{equation}{U}^{\text{eq}}=-\int \frac{\mathrm{d}\varepsilon }{2\pi \mathrm{i}}{f}_{0}\enspace \mathrm{ln}\left[\frac{\mathrm{det}\left({\boldsymbol{G}}^{a}\right)}{\mathrm{det}\left({\boldsymbol{G}}^{r}\right)}\right],\end{equation} \tag{ 34 }$$

where

$$\begin{equation}\frac{\mathrm{det}\left({\boldsymbol{G}}^{a}\right)}{\mathrm{det}\left({\boldsymbol{G}}^{r}\right)}=\frac{{\left({\epsilon}-{\tilde {{\Sigma}}}_{0}^{r}\right)}^{2}\left[{\epsilon}-{{\epsilon}}_{d}+\mathrm{i}{\tilde {{\Gamma}}}_{\eta }-{v}^{2}{\tilde {{\Sigma}}}_{0}^{r}\right]}{{\left({\epsilon}-{\tilde {{\Sigma}}}_{0}^{a}\right)}^{2}\left[{\epsilon}-{{\epsilon}}_{d}-\mathrm{i}{\tilde {{\Gamma}}}_{\eta }-{v}^{2}{\tilde {{\Sigma}}}_{0}^{a}\right]}.\end{equation} \tag{ 35 }$$

Here, we used $\mathrm{det}{\left({\boldsymbol{G}}^{r}\right)}^{-1}=\mathrm{det}\left(\varepsilon \boldsymbol{I}-{\boldsymbol{H}}_{\text{eff}}\right)$ and = (ɛ − E₀)/V₀. It is worth mentioning that equation (34) is also valid for systems coupled to more than two reservoirs where v² is given by equation (13).

In figure 3 we show the potential energy for some representative values of the normalized Fermi energy ${{\epsilon}}_{F}=\frac{\left({\mu }_{0}-{E}_{0}\right)}{{V}_{0}}$. The green and cyan lines superimposed to the plot show the different QDPTs. As can be seen, there is not a clear correlation between them and the equilibrium potential. However, the figures show a strong dependence of U^eq with the Fermi energy. This is important since then, e.g., gate voltages can be used to control the equilibrium position of the system or, even up to some extent, its dynamics. In this regard, having a simple but general expression for U^eq may result particularly useful.

Zoom In Zoom Out Reset image size

3.2. Curl and electronic friction

The matrices Λ_ν , with our choice of mechanical coordinates, see section 2.2.4, are

$$\begin{align}\hfill {\mathbf{\Lambda }}_{{V}_{L}}& =\left(\begin{matrix}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),\qquad {\mathbf{\Lambda }}_{{V}_{R}}=\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \end{matrix}\right)\hfill \\ \hfill \mathrm{a}\mathrm{n}\mathrm{d}& \quad {\mathbf{\Lambda }}_{{E}_{p}}=-\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right).\hfill \end{align} \tag{ 36 }$$

The spectral function can be calculated from $\boldsymbol{A}=2\enspace \mathrm{I}\mathrm{m}\left({\boldsymbol{G}}^{r}\right)$. Using this in equation (32), it is possible to obtain a simple expression for the elements of the curl of the forces (see appendices F and G)

$$\begin{align}\hfill \left(\nabla \right.{\times} \overrightarrow {F}\left.\right)& \simeq \overrightarrow {g}\left({{\epsilon}}_{F}\right)\tilde {T}\left({{\epsilon}}_{F}\right){N}_{d}\left({{\epsilon}}_{F}\right)\delta \mu \hfill \\ \hfill & \quad +\frac{{\pi }^{2}}{3}{\left.\frac{\partial \left( \overrightarrow {g}\tilde {T}{N}_{d}\right)}{\partial {\epsilon}}\right\vert }_{{\epsilon}={{\epsilon}}_{F}}{\left({k}_{\mathrm{B}}{T}_{0}\right)}^{2}\frac{\delta T}{{T}_{0}}.\hfill \end{align} \tag{ 37 }$$

Here we take μ₀ = (μ_L + μ_R )/2, T₀ = (T_L + T_R )/2, δμ = μ_L − μ_R , and δT = T_L − T_R . The transmittance between reservoirs L and R is $\tilde {T}\left(\varepsilon \right)$, ${N}_{d}\left(\varepsilon \right)$ is the LDOS of the dot, and

$$\begin{equation} \overrightarrow {g}=\frac{1}{{V}_{0}}\left(2\frac{\left({\epsilon}-{{\epsilon}}_{d}\right)}{{v}_{L}{v}_{R}},\frac{1}{{v}_{R}},\frac{1}{{v}_{L}}\right).\end{equation} \tag{ 38 }$$

The curl of the force represents the work per unit area where, in our case, the unit area is given in units of energy and therefore the curl has units of one over energy.

As there are three components of the curl (equation (37)) and two coefficients of the expansion (one for δT and one for δμ), there are a total of six components of the expansion of the curls. However, due to the symmetry of the problem, the coefficients of the expansion of ${\left(\nabla {\times}F\right)}_{{V}_{L}}$ are equivalent to the coefficients of the expansion of ${\left(\nabla {\times}F\right)}_{{V}_{R}}$ with v_L and v_R interchanged. Therefore, we are going to consider only four coefficients namely ${W}_{{E}_{d}}^{\left(\delta \mu \right)}$, ${W}_{{E}_{d}}^{\left(\delta T\right)}$, ${W}_{{V}_{L}}^{\left(\delta \mu \right)}$, and ${W}_{{V}_{L}}^{\left(\delta T\right)}$, where

$$\begin{align}\hfill {\left(\nabla {\times}F\right)}_{{V}_{L}}& =\left(\delta \mu \right){W}_{{V}_{L}}^{\left(\delta \mu \right)}+{\left({k}_{\mathrm{B}}{T}_{0}\right)}^{2}\frac{\delta T}{{T}_{0}}{W}_{{V}_{L}}^{\left(\delta T\right)}\hfill \\ \hfill {\left(\nabla {\times}F\right)}_{{E}_{d}}& =\left(\delta \mu \right){W}_{{E}_{d}}^{\left(\delta \mu \right)}+{\left({k}_{\mathrm{B}}{T}_{0}\right)}^{2}\frac{\delta T}{{T}_{0}}{W}_{{E}_{d}}^{\left(\delta T\right)}.\hfill \end{align} \tag{ 39 }$$

When one multiplies the coefficients ${W}_{\rho }^{\left(\delta \mu \right)}$ or ${W}_{\rho }^{\left(\delta T\right)}$ (where ρ = E_d , V_L , or V_R ) by δμ or ${\left({k}_{\mathrm{B}}{T}_{0}\right)}^{2}\frac{\delta T}{{T}_{0}}$ respectively, the result gives the low temperature limit of the work done per unit area for a cyclic movement of the other two variables.

Having a simple expression for the coefficients ${W}_{\rho }^{\left(\delta \mu \right)}$ or ${W}_{\rho }^{\left(\delta T\right)}$ may be important when studying different forms of nanoelectromechanical devices, adiabatic quantum motors, or quantum heat engines. Let us recall that due to Onsager's reciprocal relations, the coefficients of the expansion of the work in terms of δμ and δT gives, respectively, the charge pumped by adiabatic quantum pumps and the heat pumped by adiabatic quantum heat pumps [27, 32]. Therefore, the expressions found for the coefficients and the conclusions regarding their general behavior are also valid for such systems.

In figure 4 we show some examples of the behavior of the coefficients of the expansion of the curl. In general, their behavior is highly dependent on the Fermi energy, just as with the case of U^eq. This feature may be useful, e.g., to control the dynamics of a system by using a gate voltage as a knob. Indeed, it is possible to change, in a given region of the parameter's space, the sign of the coefficients, and thus the preferred direction of motion of the system by changing _F .

Zoom In Zoom Out Reset image size

The green solid lines and the cyan dashed lines indicate the QDPTs in figure 4. Although not obvious at first glance, we find that the behavior of the W coefficients is related to the type of poles of the Green's function describing the effective Hamiltonian, equation (37). This is so as the shape of the LDOS is determined by the type of pole, which, in turn, depends on the parameters of the Hamiltonian.

The ${W}_{{V}_{L}}$ coefficients depend on the square of the LDOS (or its derivative). Note that, in our system the transmittance $\tilde {T}\left(\varepsilon \right)$ is proportional to the LDOS within a wide band approximation [53]. Therefore it is expected a coincidence between the regions with the maximum value of the coefficients and the regions with a maximum value of the LDOS, see equation (37). According to the discussions of section 2.1, for resonant states (see the region enclosed by the cyan dashed line in figure 4), the maximum of the LDOS approximately coincides (for poles not too close to the band edges) with the real part of the pole, see figure 2(c). In such a case it is expected a maximum of the ${W}_{{V}_{L}}$ coefficients for _F = Re(_p ), or

$$\begin{equation}{v}_{R}^{2}=\frac{1}{1+\frac{1}{2\left(\frac{{{\epsilon}}_{F}}{{{\epsilon}}_{d}}-1\right)}}-{v}_{L}^{2},\end{equation} \tag{ 40 }$$

where we used equation (10). The above equation is plotted in figures 4(b) and (d) as violet dashed lines. As can be seen, the regions with a maximum value of the ${W}_{{V}_{L}}$ coefficients are close to the line given by equation (40).

For the ${W}_{{E}_{d}}$ coefficients the same analysis can be done except that they depend on the square of the LDOS multiplied by the term ( − _d ), which shifts the maximum of the function from Re(_p ). In this case, the maximum value of the coefficients should be close to the regions where the following function is maximum

$$\begin{equation}\left({{\epsilon}}_{F}-{{\epsilon}}_{d}\right){\left(\frac{\mathrm{I}\mathrm{m}\left({{\epsilon}}_{p}\right)}{{\left[{{\epsilon}}_{F}-\mathrm{R}\mathrm{e}\left({{\epsilon}}_{p}\right)\right]}^{2}+\mathrm{I}\mathrm{m}{\left({{\epsilon}}_{p}\right)}^{2}}\right)}^{2}.\end{equation} \tag{ 41 }$$

Here, we used the fact that in our system the transmittance is proportional to the LDOS. We assumed the pole is close to the center of the band, and thus Γ can be taken in a wide band approximation and the LDOS has a Lorentzian shape centered around Re(_p ) and with a width given by Im(_p ). And we applied all this to equation (37). Note that _p is a function of _d , v_L , and v_R , see equation (10). The above condition, numerically found, is plotted in figures 4(a) and (c) as violet dashed lines. As can be seen, the regions with a maximum value of the ${W}_{{E}_{d}}$ coefficients are close to this line.

All the above discussions are valid as long as |_F | ≪ 2. Let us recall that only when _p is far from the band edges, the LDOS can be described by a Lorentzian function centered around Re(_p ). When Re(_p ) approaches one of the band edges (|Re(_p )| ≈ 2), the LDOS is largely distorted exhibiting a peak with a maximum shifted toward the closest band edge. If we move the pole even further, it will correspond to a virtual state (see figure 2(a)). For virtual states (|Re(_p )| ⩾ 2 and Im(_p ) = 0), the LDOS shows a maximum almost at the band edges ≈ ±2, see figure 2(d). The height of this maximum grows until the virtual-localized QDPT is reached (where the maximum is exactly at = ±2) and then it decreases. Therefore, in all these cases the maximum of the LDOS is expected to approximately coincide with one of the band edges. This is the reason why for |_F | ≈ 2 the maximum value of the curl coincide approximately with the virtual-localized QDPT, see figure 5. For localized states, obviously there is also a maximum outside the band, see figure 2(e). However, as discussed in section 2.2.4, states with energies outside the band do not contribute to nonequilibrium forces. Note also that the transmittance $\tilde {T}\left(\varepsilon \right)$ is zero outside the band and thus $\left(\nabla {\times} \overrightarrow {F}\right)=0$ according to equation (37).

Zoom In Zoom Out Reset image size

In summary, as a rule of thumb when |_F | ≪ 2, the regions of the parameter space with maximum values of the ${W}_{{V}_{L}}$ coefficients are given by equation (40), while for the ${W}_{{E}_{d}}$ coefficients, these regions are shifted and given by the maximum of equation (41). On the other hand, when |_F | ≈ 2, the regions of the parameter space with the maximum values of both coefficients (${W}_{{V}_{L}}$ and ${W}_{{E}_{d}}$) approximately coincide with the virtual-localized QDPT, equation (12).

Here, we also obtained closed-form expressions for the elements of the electronic-friction tensor γ , see appendix H. Figure 6 shows some elements of γ for the same conditions used in figure 4. As can be seen, in general, the regions in the parameter space that maximize the W coefficients roughly coincide with the regions that maximize the elements of γ . This is reasonable since γ and W show a similar, although not exactly the same, dependence with the LDOS and the transmittance, see appendix H. Therefore, one can use much the same arguments to explain why some regions of the parameter space present a maximum.

Zoom In Zoom Out Reset image size

One may wonder, what are the regions in the parameter space that maximizes the W coefficients independently of _F , i.e., if one could move _F at will. This is shown in figures 7 and 8. We do not show the W^(δT) coefficients as their behavior is approximately the same as that of the W^(δμ) coefficients. As can be seen, the regions that maximize the W coefficients coincide in general with the virtual-localized QDPT. Note however that in figure 7 the region v_L ≈ v_R ≈ 0 presents values of ${W}_{{E}_{d}}^{\left(\delta \mu \right)}$ similar to those of the virtual-localized QDPT. Regarding the elements of the electronic-friction tensor, they roughly follow the behavior of the W coefficients.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, to help better understand how the above results can be applied to real physical systems, we added a discussion about a concrete tight-binding example in appendix I.

When designing or studying a nanoelectromechanical device or a quantum machine, it is crucial to understand the effect of the system's parameters on their performance. Here, we carry out a thorough study of the CIF and the electronic-friction tensor within one of the most common Hamiltonian models in quantum transport [53–55]. It is important to point out that any single-particle Hamiltonian, where there is only one relevant state, can be reduced to the model used here by means of a decimation procedure, see, e.g. reference [54].

We derived analytical formulas for the equilibrium potential (equations (23) and (35)), the low-temperature limit of the electronic-friction tensor (equation (26)), and the low-temperature limit of the curl of the CIFs up to leading order in δT and δμ (equations (37) and (32)). We showed the connection between QDPTs, the curl of the CIFs, and the electronic-friction tensor. Using this as a guide, we found some rule of thumbs that can be useful to quickly identify the regions of the parameters space that should be enclosed to maximize the work done by the device. Moreover, because of the strategy we used to carry out the analysis, our results are independent of the details of how the mechanical and electronic degrees of freedom are coupled.

Again, we want to emphasize that, due to Onsager's reciprocal relations, our analysis of the coefficients of the expansion of the curl (${W}_{\rho }^{\left(\delta \mu \right)}$ and ${W}_{\rho }^{\left(\delta T\right)}$) and the formulas shown in equations (32) and (37), are also valid for adiabatic quantum pumps and quantum heat pumps [27, 32].

The present work may also help to connect the configuration of conducting devices with current-induced structural-failures. Basically, our simple Hamiltonian also represents a minimal model of an impurity in a one-dimensional conductor. Therefore, according to our results, some particular values of the parameters of the system should dramatically increase nonconservative forces and the curl of the force, which should ultimately lead to a mechanical failure. That reasoning, of course, depends on the vibrational modes involved and the dynamical feedback between them and the electronic currents [16]. However, it is an interesting direction for further research.

Other appealing directions of extending this work include analyzing other common Hamiltonian models, such as double quantum dots, and studying systems within the Coulomb blockade regime of quantum transport.

We acknowledge financial support by Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET); Secretaría de Ciencia y Tecnología de la Universidad Nacional de Córdoba (SECYT-UNC); and Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT, PICT-2018-03587).

No new data were created or analysed in this study.

By writing the Fermi functions as f_L/R = f₀ + Δf_L/R , one can split the zero order term of the adiabatic expansion of the CIFs F_ν into the equilibrium (${F}_{\nu }^{\text{eq}}$) and nonequilibrium (${F}_{\nu }^{\text{ne}}$) contributions

$$\begin{align*}\hfill {F}_{\nu }& =\int \frac{\mathrm{d}\varepsilon }{\pi }{f}_{0}\sum\limits _{\alpha }\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{r}{\mathbf{\Gamma }}_{\alpha }{\boldsymbol{G}}^{a}\right]\hfill \\ \hfill & \quad +\int \frac{\mathrm{d}\varepsilon }{\pi }\sum\limits _{\alpha }{\Delta}{f}_{\alpha }\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{r}{\mathbf{\Gamma }}_{\alpha }{\boldsymbol{G}}^{a}\right]\hfill \\ \hfill & ={F}_{\nu }^{\text{eq}}+{F}_{\nu }^{\text{ne}}.\hfill \end{align*}$$

Given their definition, the following relation holds for the Green functions

$$\begin{equation*}{\left({\boldsymbol{G}}^{r}\right)}^{-1}-{\left({\boldsymbol{G}}^{a}\right)}^{-1}=2\mathrm{i}\sum\limits _{\alpha }{\mathbf{\Gamma }}_{\alpha }.\end{equation*}$$

Therefore, equilibrium force can be written as

$$\begin{equation*}{F}_{\nu }^{\text{eq}}=\int \frac{\mathrm{d}\varepsilon }{2\pi \mathrm{i}}{f}_{0}\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\left({\boldsymbol{G}}^{a}-{\boldsymbol{G}}^{a}\right)\right].\end{equation*}$$

Now, let us take a matrix U which diagonalize G ^r ,

$$\begin{equation*}{\boldsymbol{U}}^{-1}{\boldsymbol{G}}^{r}\boldsymbol{U}={\boldsymbol{D}}^{r}.\end{equation*}$$

Obviously, U also diagonalize ${\left({\boldsymbol{G}}^{r}\right)}^{-1}$. Then, using $\frac{\partial \boldsymbol{H}}{\partial {X}_{\nu }}=\frac{\partial {\left({\boldsymbol{G}}^{r}\right)}^{-1}}{\partial {X}_{\nu }}={\left({\boldsymbol{G}}^{r}\right)}^{-1}\frac{\partial {\boldsymbol{G}}^{r}}{\partial {X}_{\nu }}{\left({\boldsymbol{G}}^{r}\right)}^{-1}$, one finds

$$\begin{equation*}\mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{r}\right]=\mathrm{Tr}\left[-\frac{\partial {\left({\boldsymbol{D}}^{r}\right)}^{-1}}{\partial {X}_{\nu }}{\boldsymbol{D}}^{r}\right].\end{equation*}$$

The elements of D ^r are the eigenvalues of G ^r , thus

$$\begin{equation*}\mathrm{Tr}\left[-\frac{\partial {\left({\boldsymbol{D}}^{r}\right)}^{-1}}{\partial {X}_{\nu }}{\boldsymbol{D}}^{r}\right]=-\frac{\partial \mathrm{ln}\left(\mathrm{det}\left({\boldsymbol{G}}^{r}\right)\right)}{\partial {X}_{\nu }}.\end{equation*}$$

A similar argument can be used for G ^a . The equilibrium force results in

$$\begin{equation*}{F}_{\nu }^{\text{eq}}=\frac{\partial }{\partial {X}_{\nu }}\left[\int \frac{\mathrm{d}\varepsilon }{2\pi \mathrm{i}}{f}_{0}\enspace \mathrm{ln}\left(\frac{\mathrm{det}\left({\boldsymbol{G}}^{a}\right)}{\mathrm{det}\left({\boldsymbol{G}}^{r}\right)}\right)\right],\end{equation*}$$

where one can identify the equilibrium potential U^eq as

$$\begin{equation*}{U}^{\text{eq}}=-\int \frac{\mathrm{d}\varepsilon }{2\pi \mathrm{i}}{f}_{0}\enspace \mathrm{ln}\left(\frac{\mathrm{det}\left({\boldsymbol{G}}^{a}\right)}{\mathrm{det}\left({\boldsymbol{G}}^{r}\right)}\right).\end{equation*}$$

Note that $\left\vert \frac{\mathrm{det}\left({\boldsymbol{G}}^{a}\right)}{\mathrm{det}\left({\boldsymbol{G}}^{r}\right)}\right\vert =1$, and then one can write

$$\begin{equation*}\frac{\mathrm{det}\left({\boldsymbol{G}}^{a}\right)}{\mathrm{det}\left({\boldsymbol{G}}^{r}\right)}=\mathrm{exp}\left(-\mathrm{i}\chi \right).\end{equation*}$$

Taking $z=\mathrm{det}\left({\boldsymbol{G}}^{r}\right)$ and considering χ as the principal value of the logarithm, gives

$$\begin{equation*}\chi =\mathrm{arctan}\left(\frac{2\enspace \mathrm{R}\mathrm{e}\left(z\right)\enspace \mathrm{I}\mathrm{m}\left(z\right)}{{\left(\mathrm{R}\mathrm{e}\left(z\right)\right)}^{2}-{\left(\mathrm{I}\mathrm{m}\left(z\right)\right)}^{2}}\right).\end{equation*}$$

Using this, the equilibrium potential can be written as

$$\begin{equation*}{U}^{\text{eq}}=\int \frac{\mathrm{d}\varepsilon }{2\pi }{f}_{0}\chi \left(\varepsilon \right).\end{equation*}$$

Note that it is not necessary to integrate from −∞ to ∞ but only over energies where Im( G ^r ) ≠ 0 (where χ(ɛ) ≠ 0). As the function 'arctan' is bounded and monotonically increasing, the following inequalities hold

$$\begin{align*}\hfill \pi & {\geqslant}\chi {\geqslant}-\pi \hfill \\ \hfill \frac{1}{2}{\int }_{{\ast}}\mathrm{d}\varepsilon \enspace {f}_{0}& {\geqslant}{U}^{\text{eq}}{\geqslant}-\frac{1}{2}{\int }_{{\ast}}\mathrm{d}\varepsilon \enspace {f}_{0},\hfill \end{align*}$$

where the symbol ∗ means that the integrals only cover the energy ranges where Im( G ^r ) ≠ 0, or, in other words, the integrals only run over energy ranges where the leads support propagating states, or the system presents localized states. Then, the above equation shows that the equilibrium potential is bounded by the occupation of the reservoirs and the system.

As mentioned, here we are interested only in the symmetric component of the electronic-friction tensor at equilibrium, equation (24). Under this condition, the lesser and greater Green's functions are given by [1, 53]:

$$\begin{align*}\hfill {\boldsymbol{G}}_{\text{eq}}^{{< }}& =\mathrm{i}{f}_{0}\boldsymbol{A}\hfill \\ \hfill {\boldsymbol{G}}_{\text{eq}}^{{ >}}& =-\mathrm{i}\left[1-{f}_{0}\right]\boldsymbol{A},\hfill \end{align*}$$

where A is the spectral function, $\boldsymbol{A}=\mathrm{i}\left({\boldsymbol{G}}^{r}-{\boldsymbol{G}}^{a}\right)$. Using this we can write equation (24) as

$$\begin{align}\hfill {\gamma }_{\nu {\nu }^{\prime }}^{\text{eq}}=\int \frac{\hslash \enspace \mathrm{d}\varepsilon }{4\pi }{f}_{0}\left({\partial }_{\varepsilon }\left(1-{f}_{0}\right)\mathrm{Tr}\left[{\left\{{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right\}}_{s}\right]\right.& \hfill & \hfill \\ \hfill \left.+\left(1-{f}_{0}\right)\mathrm{Tr}\left[{\left\{{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right\}}_{s}{\partial }_{\varepsilon }\boldsymbol{A}\right]\right),\end{align} \tag{ B1 }$$

where we used the notation ${\left\{{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right\}}_{s}=\left({\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}+{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}\right)$.

Employing the cyclic property of the trace of matrix multiplications, one can readily show

$$\begin{equation}\mathrm{Tr}\left[{\left\{{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right\}}_{s}\right]=2\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right].\end{equation} \tag{ B2 }$$

The above equation can be used to further prove

$$\begin{equation}{\partial }_{\varepsilon }\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right]=\mathrm{Tr}\left[{\left\{{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\right\}}_{s}{\partial }_{\varepsilon }\boldsymbol{A}\right].\end{equation} \tag{ B3 }$$

Combining equations (B2) and (B3) we find

$$\begin{align}\hfill & 2{\partial }_{\varepsilon }\left(\left(1-{f}_{0}\right)\mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right]\right)=\hfill \\ \hfill & \quad \left(1-{f}_{0}\right){\partial }_{\varepsilon }\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right]\hfill \\ \hfill & \quad +\mathrm{Tr}\left[{\left\{{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\right\}}_{s}{\partial }_{\varepsilon }\left(\left(1-{f}_{0}\right)\boldsymbol{A}\right)\right].\hfill \end{align} \tag{ B4 }$$

Now, after some algebra and using equations (B2) and (B4) the following relation is obtained

$$\begin{align}\hfill {f}_{0}\enspace & \mathrm{Tr}\left[{\left\{{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\right\}}_{s}{\partial }_{\varepsilon }\left(\left(1-{f}_{0}\right)\boldsymbol{A}\right)\right]\hfill & \hfill \\ \hfill & =2{\partial }_{\varepsilon }\left({f}_{0}\left(1-{f}_{0}\right)\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right]\right)\hfill \\ \hfill & \quad -2\left({\partial }_{\varepsilon }{f}_{0}\right)\left(1-{f}_{0}\right)\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right]\hfill \\ \hfill & \quad -{f}_{0}\left(1-{f}_{0}\right){\partial }_{\varepsilon }\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right].\hfill \end{align} \tag{ B5 }$$

Inserting equation (B5) into equation (B1) yields

$$\begin{align}\hfill {\gamma }_{\nu {\nu }^{\prime }}^{\text{eq}}& ={\left.\frac{\hslash }{2\pi }{f}_{0}\left(1-{f}_{0}\right)\mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right]\right\vert }_{-\infty }^{\infty }\hfill \\ \hfill & \quad -\int \frac{\hslash \enspace \mathrm{d}\varepsilon }{2\pi }\left({\partial }_{\varepsilon }{f}_{0}\right)\left(1-{f}_{0}\right)\mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right]\hfill \\ \hfill & \quad -\int \frac{\hslash \enspace \mathrm{d}\varepsilon }{4\pi }{f}_{0}\left(1-{f}_{0}\right){\partial }_{\varepsilon }\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right].\hfill \end{align} \tag{ B6 }$$

As $\mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right]$ is bounded, and ${\mathrm{lim}}_{\varepsilon \to {\pm}\infty }\enspace {f}_{0}\left(1-{f}_{0}\right)=0$, then

$$\begin{equation}{\left.\frac{1}{2\pi }{f}_{0}\left(1-{f}_{0}\right)\mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right]\right\vert }_{-\infty }^{\infty }=0.\end{equation} \tag{ B7 }$$

Using this in equation (B6) and after some algebra we finally find

$$\begin{equation*}{\gamma }_{\nu {\nu }^{\prime }}^{\text{eq}}=\int \frac{\hslash \enspace \mathrm{d}\varepsilon }{4\pi }\left(-{\partial }_{\varepsilon }{f}_{0}\right)\mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }\boldsymbol{A}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\boldsymbol{A}\right].\end{equation*}$$

Let us consider the contribution of the α lead to the CIF (F_α ),

$$\begin{equation*}F=\sum\limits _{\alpha }{F}_{\alpha },\end{equation*}$$

and define the compact support function ${\varphi }_{\alpha }\left(\varepsilon \right)$ such that

$$\begin{equation*}{F}_{\alpha }=\underset{\mathbb{R}}{\int }{\varphi }_{\alpha }\left(\varepsilon \right){f}_{\alpha }\enspace \mathrm{d}\varepsilon ,\end{equation*}$$

where f_α is the Fermi distribution function of reservoir α and

$$\begin{equation}{\varphi }_{\alpha }=\frac{1}{\pi }\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\alpha }{\boldsymbol{G}}^{r}{\mathbf{\Gamma }}_{\alpha }{\boldsymbol{G}}^{a}\right].\end{equation} \tag{ C1 }$$

Expanding F_α up to second order in μ_α and T_α and taking the small temperature limit of T_α yields

$$\begin{align}\hfill & {F}_{\alpha }=\underset{\infty }{\overset{{\mu }_{0}}{\int }}{\varphi }_{\alpha }\left(\varepsilon \right)\mathrm{d}\varepsilon +\frac{{\pi }^{2}}{6}{\left.\frac{\partial {\varphi }_{\alpha }}{\partial \varepsilon }\right\vert }_{\varepsilon ={\mu }_{\alpha }}{\left({k}_{\mathrm{B}}{T}_{0}\right)}^{2}\hfill \\ \hfill & \quad \quad +{\left.{\varphi }_{\alpha }\right\vert }_{\varepsilon ={\mu }_{0}}\delta {\mu }_{\alpha }+\frac{{\pi }^{2}}{3}{\left.\frac{\partial {\varphi }_{\alpha }}{\partial \varepsilon }\right\vert }_{\varepsilon ={\mu }_{0}}{\left({k}_{\mathrm{B}}{T}_{0}\right)}^{2}\frac{\delta {T}_{\alpha }}{{T}_{0}}\hfill \\ \hfill & \quad \quad +\frac{1}{2}{\left.\frac{\partial {\varphi }_{\alpha }}{\partial \varepsilon }\right\vert }_{\varepsilon ={\mu }_{0}}\delta {\mu }_{\alpha }^{2}+\frac{{\pi }^{2}}{6}{\left.\frac{\partial {\varphi }_{\alpha }}{\partial \varepsilon }\right\vert }_{\varepsilon ={\mu }_{0}}{k}_{\mathrm{B}}^{2}\delta {T}_{\alpha }^{2}.\hfill \end{align} \tag{ C2 }$$

As in the main text, here we used δT_α = T_α − T₀ and δμ_α = μ_α − μ₀, where T₀ and μ₀ are the average temperature and chemical potential of the reservoirs connected to the system. Note that the first two terms of the right-hand side of equation (C2) will contribute to the equilibrium force while the last two terms are second order in an expansion in terms of δμ_α and δT_α . Therefore, by inserting equation (C1) into the above equation, taking only the first-order terms and summing up all F_α contributions, one readily arrives at equation (29).

The determinant of the effective Hamiltonian given in equation (8) is

$$\begin{align}\hfill \mathrm{det}\left(\varepsilon \boldsymbol{I}-{\boldsymbol{H}}_{\text{eff}}\right)& ={V}_{0}^{3}{\left({\epsilon}-\frac{{{\Sigma}}_{0}^{r}}{{V}_{0}}\right)}^{2}\hfill \\ \hfill & \quad {\times}\left[{\epsilon}-{{\epsilon}}_{d}+\mathrm{i}\frac{{{\Gamma}}_{\eta }^{r}}{{V}_{0}}-\left({v}_{R}^{2}+{v}_{L}^{2}\right)\frac{{{\Sigma}}_{0}^{r}}{{V}_{0}}\right].\hfill \end{align} \tag{ D1 }$$

The retarded Green function is

$$\begin{equation*}{\boldsymbol{G}}^{r}\left(\varepsilon \right)=\left(\begin{matrix}\hfill {G}_{11}^{r}\hfill & \hfill {G}_{12}^{r}\hfill & \hfill {G}_{13}^{r}\hfill \\ \hfill {G}_{21}^{r}\hfill & \hfill {G}_{22}^{r}\hfill & \hfill {G}_{23}^{r}\hfill \\ \hfill {G}_{31}^{r}\hfill & \hfill {G}_{32}^{r}\hfill & \hfill {G}_{33}^{r}\hfill \end{matrix}\right)={\left[{\epsilon}\boldsymbol{I}-{\boldsymbol{H}}_{\text{eff}}\right]}^{-1}.\end{equation*}$$

For convenience, we changed the notation for the indexes of the elements of G , with respect to that of equation (4) in the main text. Thus, e.g., $-\left(1/\pi \right)\mathrm{I}\mathrm{m}\left({G}_{2,2}^{r}\right)$ is the LDOS of the dot (site '0' according to equation (4)), while $-\left(1/\pi \right)\mathrm{I}\mathrm{m}\left({G}_{1,1}^{r}\right)$ is the LDOS of the first site of the left chain (site '−1' according to equation (4)). In all the appendixes, we followed the present convention.

Using equation (D1) one can readily calculate the elements of ${\boldsymbol{G}}^{r}\left(\varepsilon \right)$. Written them in terms of the dimensionless quantities (, _d , v_L/R , ${\tilde {{\Gamma}}}_{\eta }=\frac{{{\Gamma}}_{\eta }}{{V}_{0}}$, and ${\tilde {{\Sigma}}}_{0}^{r}=\frac{{{\Sigma}}_{0}^{r}\left(\varepsilon \right)}{{V}_{0}}$) they are

$$\begin{align*}\hfill {G}_{11}^{r}& =\frac{\left({\epsilon}-{{\epsilon}}_{d}+\mathrm{i}{\tilde {{\Gamma}}}_{\eta }-{v}_{R}^{2}{\tilde {{\Sigma}}}_{0}^{r}\right)}{{V}_{0}\left[{\epsilon}-{{\epsilon}}_{d}+\mathrm{i}{\tilde {{\Gamma}}}_{\eta }-\left({v}_{R}^{2}+{v}_{L}^{2}\right){\tilde {{\Sigma}}}_{0}^{r}\right]\left({\epsilon}-{\tilde {{\Sigma}}}_{0}^{r}\right)}\hfill \\ \hfill {G}_{22}^{r}& =\frac{1}{{V}_{0}\left[{\epsilon}-{{\epsilon}}_{d}+\mathrm{i}{\tilde {{\Gamma}}}_{\eta }-\left({v}_{R}^{2}+{v}_{L}^{2}\right){\tilde {{\Sigma}}}_{0}^{r}\right]}\hfill \\ \hfill {G}_{33}^{r}& =\frac{\left({\epsilon}-{{\epsilon}}_{d}+\mathrm{i}{\tilde {{\Gamma}}}_{\eta }-{v}_{L}^{2}{\tilde {{\Sigma}}}_{0}^{r}\right)}{{V}_{0}\left[{\epsilon}-{{\epsilon}}_{d}+{\tilde {{\Gamma}}}_{\eta }-\left({v}_{R}^{2}+{v}_{L}^{2}\right){\tilde {{\Sigma}}}_{0}^{r}\right]\left({\epsilon}-{\tilde {{\Sigma}}}_{0}^{r}\right)}\hfill \\ \hfill {G}_{21}^{r}& ={G}_{12}^{r}=\frac{-{v}_{L}{\tilde {{\Sigma}}}_{0}^{r}}{{V}_{0}\left[{\epsilon}-{{\epsilon}}_{d}+\mathrm{i}{\tilde {{\Gamma}}}_{\eta }-\left({v}_{R}^{2}+{v}_{L}^{2}\right){\tilde {{\Sigma}}}_{0}^{r}\right]}\hfill \\ \hfill {G}_{32}^{r}& ={G}_{23}^{r}=\frac{-{v}_{R}{\tilde {{\Sigma}}}_{0}^{r}}{{V}_{0}\left[{\epsilon}-{{\epsilon}}_{d}+\mathrm{i}{\tilde {{\Gamma}}}_{\eta }-\left({v}_{R}^{2}+{v}_{L}^{2}\right){\tilde {{\Sigma}}}_{0}^{r}\right]}\hfill \\ \hfill {G}_{31}^{r}& ={G}_{13}^{r}=\frac{{v}_{L}{v}_{R}{\left({\tilde {{\Sigma}}}_{0}^{r}\right)}^{2}}{{V}_{0}\left[{\epsilon}-{{\epsilon}}_{d}+\mathrm{i}{\tilde {{\Gamma}}}_{\eta }-\left({v}_{R}^{2}+{v}_{L}^{2}\right){\tilde {{\Sigma}}}_{0}^{r}\right]}.\hfill \end{align*}$$

The advanced Green function is calculated from ${\boldsymbol{G}}^{a}={\left[{\boldsymbol{G}}^{r}\right]}^{{\dagger}}$, while the lesser Green function is given by

$$\begin{equation*}{\boldsymbol{G}}^{{< }}={\boldsymbol{G}}^{r}{\mathbf{\Sigma }}^{{< }}{\boldsymbol{G}}^{a},\end{equation*}$$

where ${\mathbf{\Sigma }}^{{< }}={\mathbf{\Sigma }}_{L}^{{< }}+{\mathbf{\Sigma }}_{R}^{{< }}+{\mathbf{\Sigma }}_{\eta }^{{< }}$. We will split G ^< into two contributions (${\boldsymbol{G}}_{\left(L+R\right)}^{{< }}$ and G _η )

$$\begin{equation*}{\boldsymbol{G}}^{{< }}=\underset{{\boldsymbol{G}}_{\left(L+R\right)}^{{< }}}{\underbrace{{\boldsymbol{G}}^{r}\left({\mathbf{\Sigma }}_{L}^{{< }}+{\mathbf{\Sigma }}_{R}^{{< }}\right){\boldsymbol{G}}^{a}}}+\underset{{\boldsymbol{G}}_{\eta }^{{< }}}{\underbrace{{\boldsymbol{G}}^{r}{\mathbf{\Sigma }}_{\eta }^{{< }}{\boldsymbol{G}}^{a}}},\end{equation*}$$

where

$$\begin{align*}\hfill \left({\mathbf{\Sigma }}_{L}^{{< }}+{\mathbf{\Sigma }}_{R}^{{< }}\right)& =2\mathrm{i}\left(\begin{matrix}\hfill \left({f}_{0}+\frac{{\Delta}f}{2}\right){{\Gamma}}_{L}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \left({f}_{0}-\frac{{\Delta}f}{2}\right){{\Gamma}}_{R}\hfill \end{matrix}\right)\hfill \\ \hfill {\mathbf{\Sigma }}_{\eta }^{{< }}& =2\mathrm{i}\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {f}_{0}{{\Gamma}}_{\eta }^{\left(0\right)}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right).\hfill \end{align*}$$

Here, f_L/R = f₀ ± Δf with f₀ = (f_L + f_R )/2 and Δf = f_L − f_R . We also assumed the third lead is at equilibrium with f_η = f₀. Then, the elements of ${\boldsymbol{G}}_{\left(L+R\right)}^{{< }}$ and ${\boldsymbol{G}}_{\eta }^{{< }}$ are

$$\begin{align*}\hfill {\left[{G}_{\left(L+R\right)}^{{< }}\right]}_{ij}& =\underset{{\left[{G}_{\left(L+R\right)}^{{< }}\right]}_{ij}^{\left(\mathrm{e}\mathrm{q}\right)}}{\underbrace{2\mathrm{i}{f}_{0}\left({G}_{i1}^{r}{G}_{j1}^{r{\ast}}{{\Gamma}}_{L}+{G}_{i3}^{r}{G}_{j3}^{r{\ast}}{{\Gamma}}_{R}\right)}}\hfill \\ \hfill & \quad +\underset{{\left[{G}_{\left(L+R\right)}^{{< }}\right]}_{ij}^{\left(\mathrm{n}\mathrm{e}\right)}}{\underbrace{2\mathrm{i}\frac{{\Delta}f}{2}\left({G}_{i1}^{r}{G}_{j1}^{r{\ast}}{{\Gamma}}_{L}-{G}_{i3}^{r}{G}_{j3}^{r{\ast}}{{\Gamma}}_{R}\right)}}\hfill \end{align*}$$

and

$$\begin{equation*}{\left[{G}_{\eta }^{{< }}\right]}_{ij}=2\mathrm{i}{f}_{0}\left({G}_{i2}^{r}{G}_{j2}^{r{\ast}}{{\Gamma}}_{\eta }^{\left(0\right)}\right),\end{equation*}$$

where we identified the equilibrium (eq) and nonequilibrium (ne) contributions to ${\left[{\boldsymbol{G}}_{\left(L+R\right)}^{{< }}\right]}_{ij}$. The equilibrium contributions only contains f₀ terms while nonequilibrium contributions only contains Δf terms.

Therefore, the nonequilibrium G ^< is given by

$$\begin{equation*}{\boldsymbol{G}}_{\left(\mathrm{n}\mathrm{e}\right)}^{{< }}={\left[{\boldsymbol{G}}_{\left(L+R\right)}^{{< }}\right]}_{ij}^{\left(\mathrm{n}\mathrm{e}\right)},\end{equation*}$$

while the equilibrium contribution is

$$\begin{equation*}{\boldsymbol{G}}_{\left(\mathrm{e}\mathrm{q}\right)}^{{< }}=\begin{cases}{\left[{\boldsymbol{G}}_{\left(L+R\right)}^{{< }}\right]}^{\left(\mathrm{e}\mathrm{q}\right)}\quad \hfill & \left\vert {\epsilon}\right\vert {\leqslant}2\hfill \\ \underset{{{\Gamma}}_{\eta }\to 0}{\mathrm{lim}}{\boldsymbol{G}}_{\eta }^{{< }}\quad \hfill & \left\vert {\epsilon}\right\vert { >}2\hfill \end{cases}.\end{equation*}$$

Let $\left\{{q}_{j}\right\}$ and $\left\{{X}_{\nu }\right\}$ be two sets of different mechanical variables related through the expression ${q}_{j}\left({X}_{\nu }\right)$, which we assume continuous but not necessarily lineal. Then, CIFs can be written in terms of $\left\{{X}_{\nu }\right\}$ as $F\left({X}_{\nu }\right)$ or in terms of $\left\{{q}_{j}\right\}$ as $F\left({q}_{j}\right)$. Both descriptions are related through

$$\begin{align*}\hfill F\left({X}_{\nu }\right)& =\int \frac{\mathrm{d}\varepsilon }{2\pi \mathrm{i}}\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{{< }}\right]\hfill \\ \hfill & =\sum\limits _{j}\underset{F\left({q}_{j}\right)}{\underbrace{\left(\int \frac{\mathrm{d}\varepsilon }{2\pi \mathrm{i}}\enspace \mathrm{Tr}\left[{\mathbf{\Lambda }}_{j}{\boldsymbol{G}}^{{< }}\right]\right)}}\frac{\partial {q}_{j}}{\partial {X}_{\nu }},\hfill \end{align*}$$

where $\frac{\partial {q}_{j}}{\partial {X}_{\nu }}$ is the covariant derivative.

The work done by the forces over a closed loop C in the space of the variables is independent of the choice of the sets of mechanical variables used to describe the system, as can be readily shown

$$\begin{align*}\hfill W\left( \overrightarrow {X}\right)& =\underset{C\left( \overrightarrow {X}\right)}{\int }\sum\limits _{\nu }F\left({X}_{\nu }\right)\mathrm{d}{X}_{\nu }\hfill \\ \hfill & =\sum\limits _{\nu }\sum\limits _{j}\underset{C\left( \overrightarrow {q}\right)}{\int }F\left({q}_{j}\right)\frac{\partial {q}_{j}}{\partial {X}_{\nu }}\mathrm{d}{X}_{\nu }\hfill \\ \hfill & =\underset{C\left( \overrightarrow {q}\right)}{\int }\sum\limits _{j}F\left({q}_{j}\right)\mathrm{d}{q}_{j}=W\left( \overrightarrow {q}\right).\hfill \end{align*}$$

In the main text, we used this property to choose a set of convenient variables to become independent of the particular way in which mechanical variables affect the electronic Hamiltonian.

Our goal here is to obtained a closed-form expression for ${\left(\nabla {\times}F\right)}_{\rho }$. We start by evaluating the derivative of the force ∂_ν' F_ν = ∂F_ν /∂X_ν'

$$\begin{equation*}{\partial }_{{\nu }^{\prime }}{F}_{\nu }=\frac{1}{2\pi \mathrm{i}}\underset{-\infty }{\overset{\infty }{\int }}{\partial }_{{\nu }^{\prime }}\left\{\mathrm{Tr}\left({\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{{< }}\right)\right\}\mathrm{d}\varepsilon .\end{equation*}$$

The derivative of the lesser Green function is

$$\begin{equation*}{\partial }_{{\nu }^{\prime }}{\boldsymbol{G}}^{{< }}={\boldsymbol{G}}^{r}{\mathbf{\Lambda }}_{{\nu }^{\prime }}{\boldsymbol{G}}^{{< }}+{\boldsymbol{G}}^{{< }}{\mathbf{\Lambda }}_{{\nu }^{\prime }}{G}^{a},\end{equation*}$$

where we assume ∂_ν' Σ^< = 0 and used ∂_ν' G ^r/a = G ^r/a Λ_ν' G ^r/a . Then, the derivative of the trace is

$$\begin{equation*}\begin{matrix}\hfill {\partial }_{{\nu }^{\prime }}\enspace \mathrm{Tr}\left({\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{{< }}\right)=\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \mathrm{Tr}\left(\left[{\partial }_{{\nu }^{\prime }}{\mathbf{\Lambda }}_{\nu }-{\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{r}{\mathbf{\Lambda }}_{{\nu }^{\prime }}-{\left({\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{r}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\right)}^{{\dagger}}\right]{\boldsymbol{G}}^{{< }}\right),\hfill \end{matrix}\end{equation*}$$

where we used ${\mathbf{\Lambda }}_{{\nu }^{\prime }}^{{\dagger}}={\mathbf{\Lambda }}_{{\nu }^{\prime }}$ and ${\left[{\boldsymbol{G}}^{r}\right]}^{{\dagger}}={\boldsymbol{G}}^{a}$ to make ${\mathbf{\Lambda }}_{{\nu }^{\prime }}{G}^{a}{\mathbf{\Lambda }}_{\nu }={\left[{\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{r}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\right]}^{{\dagger}}$. Using this, the derivative of the force can be written as

$$\begin{align*}\hfill {\partial }_{{\nu }^{\prime }}{F}_{\nu }& =\frac{1}{2\pi \mathrm{i}}\underset{-\infty }{\overset{\infty }{\int }}\enspace \mathrm{Tr}\left(\left[{\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{r}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\right.\right.\hfill \\ \hfill & \quad \left.\left.+{\left({\mathbf{\Lambda }}_{\nu }{\boldsymbol{G}}^{r}{\mathbf{\Lambda }}_{{\nu }^{\prime }}\right)}^{{\dagger}}-{\partial }_{{\nu }^{\prime }}{\mathbf{\Lambda }}_{\nu }\right]{\boldsymbol{G}}^{{< }}\right)\mathrm{d}\varepsilon .\hfill \end{align*}$$

Now as ∂_ν' Λ_ν = ∂_ν Λ_ν', the following holds

$$\begin{equation*}{\partial }_{{\nu }^{\prime }}{F}_{\nu }-{\partial }_{\nu }{F}_{{\nu }^{\prime }}=\frac{1}{2\pi \mathrm{i}}\underset{-\infty }{\overset{\infty }{\int }}\enspace \mathrm{Tr}\left(\left[{\mathbf{\Pi }}_{{\nu }^{\prime }\nu }-{\mathbf{\Pi }}_{\nu {\nu }^{\prime }}\right]{\boldsymbol{G}}^{{< }}\right)\mathrm{d}\varepsilon \end{equation*}$$

where, to compact the notation, we define

$$\begin{equation*}{\mathbf{\Pi }}_{{\nu }^{\prime }\nu }={\mathbf{\Lambda }}_{{\nu }^{\prime }}{\boldsymbol{G}}^{r}{\mathbf{\Lambda }}_{\nu }+{\left({\mathbf{\Lambda }}_{{\nu }^{\prime }}{\boldsymbol{G}}^{r}{\mathbf{\Lambda }}_{\nu }\right)}^{{\dagger}}.\end{equation*}$$

Finally, the curl of the force yields

$$\begin{equation*}{\left(\nabla {\times} \overrightarrow {F}\right)}_{\rho }=\frac{1}{2\pi \mathrm{i}}\underset{-\infty }{\overset{\infty }{\int }}\enspace \mathrm{Tr}\left({{\epsilon}}^{\rho {\nu }^{\prime }\nu }{\mathbf{\Pi }}_{{\nu }^{\prime }\nu }{\boldsymbol{G}}^{{< }}\right)\mathrm{d}\varepsilon ,\end{equation*}$$

where ^ρμν is the Levi-Civita symbol.

We start by defining the index order of the variables

$$\begin{align*}\hfill {E}_{d}& \to \#1\hfill \\ \hfill {V}_{L}& \to \#2\hfill \\ \hfill {V}_{R}& \to \#3.\hfill \end{align*}$$

With this order, the following relations hold

$$\begin{align*}\hfill {{\epsilon}}^{{E}_{d}\mu \nu }{\mathbf{\Pi }}_{\mu \nu }& ={\left\{{\mathbf{\Lambda }}_{{V}_{R}}\left({\boldsymbol{G}}^{r}-{\boldsymbol{G}}^{a}\right){\mathbf{\Lambda }}_{{V}_{L}}\right\}}_{a}\hfill \\ \hfill {{\epsilon}}^{{V}_{L}\mu \nu }{\mathbf{\Pi }}_{\mu \nu }& ={\left\{{\mathbf{\Lambda }}_{{E}_{d}}\left({\boldsymbol{G}}^{r}-{\boldsymbol{G}}^{a}\right){\mathbf{\Lambda }}_{{V}_{R}}\right\}}_{a}\hfill \\ \hfill {{\epsilon}}^{{V}_{R}\mu \nu }{\mathbf{\Pi }}_{\mu \nu }& ={\left\{{\mathbf{\Lambda }}_{{V}_{L}}\left({\boldsymbol{G}}^{r}-{\boldsymbol{G}}^{a}\right){\mathbf{\Lambda }}_{{E}_{d}}\right\}}_{a},\hfill \end{align*}$$

where ${\left\{{\boldsymbol{\lambda }}_{\nu }\dots {\boldsymbol{\lambda }}_{{\nu }^{\prime }}\right\}}_{a}=\left({\boldsymbol{\lambda }}_{\nu }\dots {\boldsymbol{\lambda }}_{{\nu }^{\prime }}-{\boldsymbol{\lambda }}_{{\nu }^{\prime }}\dots {\boldsymbol{\lambda }}_{\nu }\right)$.

Defining ${{\Sigma}}_{L/R}={v}_{L/R}^{2}{{\Sigma}}_{0}^{r}$ and Σ_L/R = Δ_L/R − iΓ_L/R , one can deduce the following

$$\begin{equation}\frac{1}{{V}_{L/R}^{2}}{\left\vert {{\Sigma}}_{L/R}\right\vert }^{2}{{\Gamma}}_{0}=\frac{{V}_{L/R}^{2}{{\Gamma}}_{0}}{{\left\vert \varepsilon -\left({E}_{0}+{{\Sigma}}_{0}^{r}\right)\right\vert }^{2}}\end{equation} \tag{ G1 }$$

$$\begin{equation}={{\Gamma}}_{L/R}.\end{equation} \tag{ G2 }$$

This relation will be useful afterward. Another formula that will be useful is the expression for the elements of G^<. For simplicity, we will assume Γ_η = 0 from the beginning, as we are interested only in the nonequilibrium contribution of the force, see appendix D. Then, the lesser self-energy yields

$$\begin{equation*}{\mathbf{\Sigma }}^{{< }}=2\mathrm{i}\left(\begin{matrix}\hfill {f}_{L}{{\Gamma}}_{L}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {f}_{R}{{\Gamma}}_{R}\hfill \end{matrix}\right).\end{equation*}$$

It will be helpful to define ${\boldsymbol{G}}^{\lesssim }=\frac{1}{2\mathrm{i}}{\boldsymbol{G}}^{{< }}$. In our case, the elements of G ^≲ are

$$\begin{equation}{G}_{ij}^{\lesssim }=\left\{{G}_{i1}^{r}{G}_{j1}^{r{\ast}}{f}_{L}{{\Gamma}}_{L}+{G}_{i3}^{r}{G}_{j3}^{r{\ast}}{f}_{R}{{\Gamma}}_{R}\right\}.\end{equation} \tag{ G3 }$$

Note that ${G}_{ij}^{\lesssim {\ast}}={G}_{ji}^{\lesssim }$. In the following, we will use all the above properties and definitions to derive the expression for each of the elements of the curl.

G.1. Element ρ = E_d

This element of the curl is given by

$$\begin{equation*}{\left(\nabla {\times} \overrightarrow {F}\right)}_{{E}_{d}}=\frac{1}{2\pi \mathrm{i}}\underset{-\infty }{\overset{\infty }{\int }}\enspace \mathrm{Tr}\left[\left({{\epsilon}}^{{E}_{d}\mu \nu }{\mathbf{\Pi }}_{\mu \nu }\right){\boldsymbol{G}}^{{< }}\right],\end{equation*}$$

where

$$\begin{equation*}{{\epsilon}}^{{E}_{d}\mu \nu }{\mathbf{\Pi }}_{\mu \nu }={\left\{{\mathbf{\Lambda }}_{{V}_{R}}\left({\boldsymbol{G}}^{r}-{\boldsymbol{G}}^{a}\right){\mathbf{\Lambda }}_{{V}_{L}}\right\}}_{a}.\end{equation*}$$

Now, the terms of the right-hand side of the above equation are

$$\begin{equation*}{\mathbf{\Lambda }}_{{V}_{R}}\left({\boldsymbol{G}}^{r}-{\boldsymbol{G}}^{a}\right){\mathbf{\Lambda }}_{{V}_{L}}=\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {a}_{32}\hfill & \hfill {a}_{31}\hfill & \hfill 0\hfill \\ \hfill {a}_{22}\hfill & \hfill {a}_{21}\hfill & \hfill 0\hfill \end{matrix}\right)\end{equation*}$$

$$\begin{equation*}{\mathbf{\Lambda }}_{{V}_{L}}\left({\boldsymbol{G}}^{r}-{\boldsymbol{G}}^{a}\right){\mathbf{\Lambda }}_{{V}_{R}}=\left(\begin{matrix}\hfill 0\hfill & \hfill {a}_{23}\hfill & \hfill {a}_{22}\hfill \\ \hfill 0\hfill & \hfill {a}_{13}\hfill & \hfill {a}_{12}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),\end{equation*}$$

and therefore

$$\begin{equation*}{{\epsilon}}^{{E}_{d}\mu \nu }{\mathbf{\Pi }}_{\mu \nu }=\left(\begin{matrix}\hfill 0\hfill & \hfill -{a}_{23}\hfill & \hfill -{a}_{22}\hfill \\ \hfill {a}_{32}\hfill & \hfill {a}_{31}-{a}_{13}\hfill & \hfill -{a}_{12}\hfill \\ \hfill {a}_{22}\hfill & \hfill {a}_{21}\hfill & \hfill 0\hfill \end{matrix}\right).\end{equation*}$$

Here, we define ${a}_{ij}=\left({G}_{ij}^{r}-{G}_{ji}^{r{\ast}}\right)$ to compact the notation. As the retarded Green functions is symmetric in our case, the following relations holds ${a}_{ij}=2\mathrm{i}\enspace \mathrm{I}\mathrm{m}\left\{{G}_{ij}^{r}\right\}$, a₃₁ = a₁₃, a₃₂ = a₂₃, and a₂₁ = a₁₂. Using this, we obtain

$$\begin{align}\hfill \mathrm{Tr}\left[\left({{\epsilon}}^{{E}_{d}\mu \nu }{\mathbf{\Pi }}_{\mu \nu }\right){\boldsymbol{G}}^{{< }}\right]=& \hfill & \hfill \\ \hfill 2\mathrm{i}\left\{{a}_{23}{b}_{12}+{a}_{22}{b}_{13}+{a}_{12}{b}_{23}\right\}.\end{align} \tag{ G4 }$$

where we defined ${b}_{ij}=\left({G}_{ij}^{\lesssim }-{G}_{ij}^{{\lesssim }^{{\ast}}}\right)=2\mathrm{i}\enspace \mathrm{I}\mathrm{m}\left\{{G}_{ij}^{\lesssim }\right\}$ with ${G}_{ij}^{\lesssim }$ given by equation (G3). After some algebra, the right-hand side terms of the above equation can be written as

$$\begin{align*}\hfill {a}_{23}{b}_{12}& =-4\enspace \mathrm{I}\mathrm{m}\left\{{G}_{23}^{r}\right\}\hfill \\ \hfill & \quad {\times}\left(\mathrm{I}\mathrm{m}\left\{{G}_{11}^{r}{G}_{21}^{r{\ast}}\right\}{f}_{L}{{\Gamma}}_{L}+\mathrm{I}\mathrm{m}\left\{{G}_{13}^{r}{G}_{23}^{r{\ast}}\right\}{f}_{R}{{\Gamma}}_{R}\right)\hfill \end{align*}$$

$$\begin{align*}\hfill {a}_{22}{b}_{13}& =-4\enspace \mathrm{I}\mathrm{m}\left\{{G}_{22}^{r}\right\}\hfill \\ \hfill & \quad {\times}\left(\mathrm{I}\mathrm{m}\left\{{G}_{11}^{r}{G}_{31}^{r{\ast}}\right\}{f}_{L}{{\Gamma}}_{L}+\mathrm{I}\mathrm{m}\left\{{G}_{13}^{r}{G}_{33}^{r{\ast}}\right\}{f}_{R}{{\Gamma}}_{R}\right)\hfill \end{align*}$$

$$\begin{align*}\hfill {a}_{12}{b}_{23}& =-4\enspace \mathrm{I}\mathrm{m}\left\{{G}_{12}^{r}\right\}\hfill \\ \hfill & \quad {\times}\left(\mathrm{I}\mathrm{m}\left\{{G}_{21}^{r}{G}_{31}^{{r}^{{\ast}}}\right\}{f}_{L}{{\Gamma}}_{L}+\mathrm{I}\mathrm{m}\left\{{G}_{23}^{r}{G}_{33}^{{r}^{{\ast}}}\right\}{f}_{R}{{\Gamma}}_{R}\right).\hfill \end{align*}$$

Using the expression for G ^r given in appendix D and equation (G1) we finds

$$\begin{align}\hfill \mathrm{I}\mathrm{m}\left\{{G}_{11}^{r}{G}_{21}^{r{\ast}}\right\}{f}_{L}{{\Gamma}}_{L}+\mathrm{I}\mathrm{m}\left\{{G}_{13}^{r}{G}_{23}^{r{\ast}}\right\}{f}_{R}{{\Gamma}}_{R}& \hfill \\ \hfill =\frac{{{\Gamma}}_{L}{{\Gamma}}_{R}}{{\left\vert \varepsilon -{E}_{d}-{{\Sigma}}_{L}-{{\Sigma}}_{R}\right\vert }^{2}}\frac{\left[{f}_{R}-{f}_{L}\right]}{{V}_{L}},\end{align} \tag{ G5 }$$

$$\begin{align}\hfill \mathrm{I}\mathrm{m}\left\{{G}_{11}^{r}{G}_{31}^{r{\ast}}\right\}{f}_{L}{{\Gamma}}_{L}+\mathrm{I}\mathrm{m}\left\{{G}_{13}^{r}{G}_{33}^{r{\ast}}\right\}{f}_{R}{{\Gamma}}_{R}& \hfill \\ \hfill =-\frac{{{\Gamma}}_{L}{{\Gamma}}_{R}}{{\left\vert \varepsilon -{E}_{d}-{{\Sigma}}_{L}-{{\Sigma}}_{R}\right\vert }^{2}}\frac{\left(\varepsilon -{E}_{d}\right)}{{V}_{L}{V}_{R}}\left[{f}_{R}-{f}_{L}\right]\end{align} \tag{ G6 }$$

and

$$\begin{align}\hfill \mathrm{I}\mathrm{m}\left\{{G}_{21}^{r}{G}_{31}^{r{\ast}}\right\}{f}_{L}{{\Gamma}}_{L}+\mathrm{I}\mathrm{m}\left\{{G}_{23}^{r}{G}_{33}^{r{\ast}}\right\}{f}_{R}{{\Gamma}}_{R}& =\hfill \\ \hfill \frac{{{\Gamma}}_{L}{{\Gamma}}_{R}}{{\left\vert \varepsilon -{E}_{d}-{{\Sigma}}_{L}-{{\Sigma}}_{R}\right\vert }^{2}}\frac{\left[{f}_{R}-{f}_{L}\right]}{{V}_{R}}.\end{align} \tag{ G7 }$$

Now introducing the transmittance $\tilde {T}\left(\varepsilon \right)$ [53, 54]

$$\begin{equation}\begin{matrix}\hfill \tilde {T}\left(\varepsilon \right)\hfill & \hfill =\hfill & \hfill 2{{\Gamma}}_{L}{\left\vert {G}_{22}^{r}\right\vert }^{2}2{{\Gamma}}_{R}\hfill \\ \hfill \hfill & \hfill =\hfill & \hfill 4{{\Gamma}}_{L}\frac{1}{{\left\vert \varepsilon -{E}_{d}-{{\Sigma}}_{L}-{{\Sigma}}_{R}\right\vert }^{2}}{{\Gamma}}_{R}\hfill \end{matrix}\end{equation} \tag{ G8 }$$

and using equations (G5)–(G7) into equation (G4) gives

$$\begin{align*}\hfill \mathrm{Tr}\left[\left({{\epsilon}}^{{E}_{d}\mu \nu }{\mathbf{\Pi }}_{\mu \nu }\right){\boldsymbol{G}}^{{< }}\right]=-2\mathrm{i}\left\{\frac{\mathrm{I}\mathrm{m}\left\{{G}_{23}^{r}\right\}}{{V}_{L}}{\times}\right.\\ \hfill \left.-\frac{\mathrm{I}\mathrm{m}\left\{{G}_{22}^{r}\right\}}{{V}_{L}{V}_{R}}\left(\varepsilon -{E}_{d}\right)+\frac{\mathrm{I}\mathrm{m}\left\{{G}_{12}^{r}\right\}}{{V}_{R}}\right\}\tilde {T}\left(\varepsilon \right)\left[{f}_{R}-{f}_{L}\right].\end{align*}$$

Using the expressions for the elements of G ^r , it is not difficult to prove:

$$\begin{equation*}\mathrm{I}\mathrm{m}\left\{{G}_{23}^{r}\right\}\frac{1}{{V}_{L}}+\mathrm{I}\mathrm{m}\left\{{G}_{12}^{r}\right\}\frac{1}{{V}_{R}}=-\frac{\left(\varepsilon -{E}_{d}\right)}{{V}_{L}{V}_{R}}\mathrm{I}\mathrm{m}\left\{{G}_{22}^{r}\right\}.\end{equation*}$$

Therefore, the element ρ = E_d of the curl can be written as

$$\begin{equation*}{\left(\nabla {\times} \overrightarrow {F}\right)}_{{E}_{d}}=-2\underset{-\infty }{\overset{\infty }{\int }}\frac{\left(\varepsilon -{E}_{d}\right)}{{V}_{L}{V}_{R}}{N}_{d}\tilde {T}\left(\varepsilon \right)\left[{f}_{R}-{f}_{L}\right]\mathrm{d}\varepsilon ,\end{equation*}$$

where N_d is the LDOS of the dot

$$\begin{equation}\begin{matrix}\hfill {N}_{d}=-\frac{1}{\pi }\underset{\eta \to {0}^{+}}{\mathrm{lim}}\enspace \mathrm{I}\mathrm{m}\left\{{G}_{22}^{r}\left(\varepsilon +\mathrm{i}\eta \right)\right\}\hfill \end{matrix}.\end{equation} \tag{ G9 }$$

G.2. Element ρ = V_L

This element of the curl is given by

$$\begin{equation*}{\left({\nabla }_{ \overrightarrow {X}}{\times} \overrightarrow {F}\right)}_{{V}_{L}}=\frac{1}{2\pi \mathrm{i}}\underset{-\infty }{\overset{\infty }{\int }}\enspace \mathrm{Tr}\left[\left({{\epsilon}}^{{V}_{L}\mu \nu }{\mathbf{\Pi }}_{\mu \nu }\right){\boldsymbol{G}}^{{< }}\right],\end{equation*}$$

where

$$\begin{equation*}{{\epsilon}}^{{V}_{L}\mu \nu }{\mathbf{\Pi }}_{\mu \nu }={\left\{{\mathbf{\Lambda }}_{{E}_{d}}\left({\boldsymbol{G}}^{r}-{\boldsymbol{G}}^{a}\right){\mathbf{\Lambda }}_{{V}_{R}}\right\}}_{a}.\end{equation*}$$

Now, the terms of the right-hand side of the above equation are

$$\begin{equation*}{\mathbf{\Lambda }}_{{E}_{d}}\left({\boldsymbol{G}}^{r}-{\boldsymbol{G}}^{a}\right){\mathbf{\Lambda }}_{{V}_{R}}=-\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {a}_{23}\hfill & \hfill {a}_{22}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right)\end{equation*}$$

$$\begin{equation*}{\mathbf{\Lambda }}_{{V}_{R}}\left({\boldsymbol{G}}^{r}-{\boldsymbol{G}}^{a}\right){\mathbf{\Lambda }}_{{E}_{d}}=-\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {a}_{32}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {a}_{22}\hfill & \hfill 0\hfill \end{matrix}\right).\end{equation*}$$

Therefore

$$\begin{equation*}{{\epsilon}}^{{V}_{L}\mu \nu }{\mathbf{\Pi }}_{\mu \nu }=\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {a}_{32}-{a}_{23}\hfill & \hfill -{a}_{22}\hfill \\ \hfill 0\hfill & \hfill {a}_{22}\hfill & \hfill 0\hfill \end{matrix}\right).\end{equation*}$$

Using this we obtain

$$\begin{equation*}\mathrm{Tr}\left[\left({{\epsilon}}^{{V}_{L}\mu \nu }{\mathbf{\Pi }}_{\mu \nu }\right){\boldsymbol{G}}^{{< }}\right]=2\mathrm{i}{a}_{22}{b}_{23},\end{equation*}$$

where the elements of the matrix G^≲ are given in equation (G3). After some algebra it is not difficult to show

$$\begin{equation*}\begin{aligned}\hfill {a}_{22}{b}_{23}& =-4\enspace \mathrm{I}\mathrm{m}\left\{{G}_{22}^{r}\right\}\hfill \\ \hfill & \quad {\times}\left(\mathrm{I}\mathrm{m}\left\{{G}_{21}^{r}{G}_{31}^{r{\ast}}\right\}{f}_{L}{{\Gamma}}_{L}+\mathrm{I}\mathrm{m}\left\{{G}_{23}^{r}{G}_{33}^{r{\ast}}\right\}{f}_{R}{{\Gamma}}_{R}\right)\hfill \end{aligned}.\end{equation*}$$

Using the expression for G ^r given in appendix D and equation (G1) one finds

$$\begin{align}\hfill \mathrm{I}\mathrm{m}\left\{{G}_{21}^{r}{G}_{31}^{r{\ast}}\right\}{f}_{L}{{\Gamma}}_{L}+\mathrm{I}\mathrm{m}\left\{{G}_{23}^{r}{G}_{33}^{r{\ast}}\right\}{f}_{R}{{\Gamma}}_{R}& =\hfill \\ \hfill \frac{{{\Gamma}}_{L}{{\Gamma}}_{R}}{{\left\vert \varepsilon -{E}_{d}-{{\Sigma}}_{L}-{{\Sigma}}_{R}\right\vert }^{2}}\frac{1}{{V}_{R}}\left[{f}_{R}-{f}_{L}\right]\end{align} \tag{ G10 }$$

Using equations (G8)–(G10) into equation (G10), we obtain

$$\begin{equation*}\mathrm{Tr}\left[\left({{\epsilon}}^{{V}_{L}\mu \nu }{{\Pi}}_{\mu \nu }\right){G}^{{< }}\right]=-2\pi \mathrm{i}\frac{1}{{V}_{R}}{N}_{d}\tilde {T}\left(\varepsilon \right)\left[{f}_{R}-{f}_{L}\right].\end{equation*}$$

Therefore, the element ρ = V_L of the curl can be written as

$$\begin{equation*}{\left(\nabla {\times} \overrightarrow {F}\right)}_{{V}_{L}}=-\frac{1}{{V}_{R}}\underset{-\infty }{\overset{\infty }{\int }}{N}_{d}\tilde {T}\left(\varepsilon \right)\left[{f}_{R}-{f}_{L}\right]\mathrm{d}\varepsilon .\end{equation*}$$

G.3. Element ρ = V_R

The expression for the element ρ = V_R of the curl is obtained by following a similar procedure to that of ${\left(\nabla {\times} \overrightarrow {F}\right)}_{{V}_{L}}$ but replacing V_L by V_R . The result is

$$\begin{equation*}{\left(\nabla {\times} \overrightarrow {F}\right)}_{{V}_{R}}=-\frac{1}{{V}_{L}}\underset{-\infty }{\overset{\infty }{\int }}{N}_{d}\tilde {T}\left(\varepsilon \right)\left[{f}_{R}-{f}_{L}\right]\mathrm{d}\varepsilon .\end{equation*}$$

G.4. Final expression of the curl

The results of appendices G.1–G.3, can be written in the compact form

$$\begin{equation*}\left(\nabla {\times} \overrightarrow {F}\right)=\underset{-\infty }{\overset{\infty }{\int }} \overrightarrow {g}\left({E}_{d},{V}_{L},{V}_{R}\right){N}_{d}\tilde {T}\left(\varepsilon \right)\left[{f}_{L}-{f}_{R}\right]\mathrm{d}\varepsilon ,\end{equation*}$$

where

$$\begin{equation*} \overrightarrow {g}\left({E}_{d},{V}_{L},{V}_{R}\right)=\left(2\frac{\varepsilon -{E}_{d}}{{V}_{L}{V}_{R}},\frac{1}{{V}_{R}},\frac{1}{{V}_{L}}\right).\end{equation*}$$

The three by three equilibrium friction tensor is symmetric, thus only six elements are necessary to describe it. To derive the expression for the elements of γ we used the following relations

$$\begin{equation*}\begin{matrix}\hfill {\mathbf{\Lambda }}_{{V}_{L}}\boldsymbol{A}\hfill & \hfill =\hfill & \hfill -2\left(\begin{matrix}\hfill \mathrm{I}\mathrm{m}\enspace {G}_{12}^{r}\hfill & \hfill \mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\hfill & \hfill \mathrm{I}\mathrm{m}\enspace {G}_{23}^{r}\hfill \\ \hfill \mathrm{I}\mathrm{m}\enspace {G}_{11}^{r}\hfill & \hfill \mathrm{I}\mathrm{m}\enspace {G}_{12}^{r}\hfill & \hfill \mathrm{I}\mathrm{m}\enspace {G}_{13}^{r}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right)\hfill \end{matrix},\end{equation*}$$

$$\begin{equation*}\begin{matrix}\hfill {\mathbf{\Lambda }}_{{V}_{R}}\boldsymbol{A}\hfill & \hfill =\hfill & \hfill -2\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \mathrm{I}\mathrm{m}\enspace {G}_{13}^{r}\hfill & \hfill \mathrm{I}\mathrm{m}\enspace {G}_{23}^{r}\hfill & \hfill \mathrm{I}\mathrm{m}\enspace {G}_{33}^{r}\hfill \\ \hfill \mathrm{I}\mathrm{m}\enspace {G}_{12}^{r}\hfill & \hfill \mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\hfill & \hfill \mathrm{I}\mathrm{m}\enspace {G}_{23}^{r}\hfill \end{matrix}\right)\hfill \end{matrix}\end{equation*}$$

and

$$\begin{equation*}\begin{matrix}\hfill {\mathbf{\Lambda }}_{{E}_{d}}\boldsymbol{A}\hfill & \hfill =\hfill & \hfill 2\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \mathrm{I}\mathrm{m}\enspace {G}_{12}^{r}\hfill & \hfill \mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\hfill & \hfill \mathrm{I}\mathrm{m}\enspace {G}_{23}^{r}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right)\hfill \end{matrix}.\end{equation*}$$

The above come from the definition of the spectral function, $\boldsymbol{A}=\mathrm{i}\left({\boldsymbol{G}}^{r}-{\boldsymbol{G}}^{a}\right)$, and the matrices Λ_ν , equation (36). To evaluate the elements of the above matrices one can use the explicit expressions for the elements of the retarded Green's function given in appendix D. After some algebra we find

$$\begin{align}\hfill \mathrm{I}\mathrm{m}\enspace {G}_{11}^{r}& =-\frac{{{\Gamma}}_{R}}{{V}_{0}^{2}}\frac{1}{{\left\vert \varepsilon -{E}_{d}-{{\Sigma}}_{R}-{{\Sigma}}_{L}\right\vert }^{2}}\hfill \\ \hfill & \quad {\times}\left\{{\left(\varepsilon -{E}_{d}\right)}^{2}{\left(\frac{{V}_{R}}{{V}_{0}}\right)}^{-2}-2\left(\varepsilon -{E}_{d}\right){{\Delta}}_{0}\right.\hfill \\ \hfill & \quad \left.+\left[{\left(\frac{{V}_{R}}{{V}_{0}}\right)}^{2}+{\left(\frac{{V}_{L}}{{V}_{0}}\right)}^{2}\right]{\left\vert {{\Sigma}}_{0}^{r}\right\vert }^{2}\right\}\hfill \end{align} \tag{ H1 }$$

$$\begin{equation}\mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}=-\frac{1}{{\left\vert \varepsilon -{E}_{d}-{{\Sigma}}_{R}-{{\Sigma}}_{L}\right\vert }^{2}}\left[{{\Gamma}}_{L}+{{\Gamma}}_{R}\right]\end{equation} \tag{ H2 }$$

$$\begin{equation}\mathrm{I}\mathrm{m}\enspace {G}_{12}^{r}=\frac{\left(\varepsilon -{E}_{d}\right)}{{V}_{L}}\frac{1}{{\left\vert \varepsilon -{E}_{d}-{{\Sigma}}_{R}-{{\Sigma}}_{L}\right\vert }^{2}}{{\Gamma}}_{L}\end{equation} \tag{ H3 }$$

$$\begin{align}\hfill \mathrm{I}\mathrm{m}\enspace {G}_{13}^{r}& =\frac{1}{{V}_{L}{V}_{R}}\frac{{\left(\frac{{V}_{L}}{{V}_{0}}\right)}^{2}{\left(\frac{{V}_{R}}{{V}_{0}}\right)}^{2}{{\Gamma}}_{0}}{{\left\vert \varepsilon -{E}_{d}-{{\Sigma}}_{R}-{{\Sigma}}_{L}\right\vert }^{2}}\left\{2\left({E}_{d}-\varepsilon \right){{\Delta}}_{0}\right.\hfill \\ \hfill & \quad \left.+{\left\vert {{\Sigma}}_{0}^{r}\right\vert }^{2}\left[{\left(\frac{{V}_{R}}{{V}_{0}}\right)}^{2}+{\left(\frac{{V}_{L}}{{V}_{0}}\right)}^{2}\right]\right\}\hfill \end{align} \tag{ H4 }$$

$$\begin{equation}\mathrm{I}\mathrm{m}\enspace {G}_{23}^{r}=\frac{\left(\varepsilon -{E}_{d}\right)}{{V}_{R}}\frac{1}{{\left\vert \varepsilon -{E}_{d}-{{\Sigma}}_{R}-{{\Sigma}}_{L}\right\vert }^{2}}{{\Gamma}}_{R}.\end{equation} \tag{ H5 }$$

Using the above formulas, and the expressions for the transmittance $\tilde {T}\left(\varepsilon \right)$ and the LDOS ${N}_{d}\left(\varepsilon \right)$ (equations (G8) and (G9) respectively), the elements of the electronic-friction tensor yield:

H.1. Element ${\gamma }_{{E}_{d}{E}_{d}}^{\text{eq}}$

$$\begin{align}\hfill \underset{{k}_{\mathrm{B}}{T}_{0}\to {0}^{+}}{\mathrm{lim}}\enspace {\gamma }_{{E}_{d}{E}_{d}}^{\text{eq}}& =\hfill & \frac{1}{4\pi }\enspace \mathrm{Tr}{\left[{{\Lambda}}_{{E}_{d}}A{{\Lambda}}_{{E}_{d}}A\right]}_{\varepsilon ={\mu }_{0}}\hfill \\ \hfill & =\hfill & \pi {N}_{d}\left({\mu }_{0}\right){N}_{d}\left({\mu }_{0}\right)\hfill \end{align},$$

where we used $\begin{matrix}\hfill \mathrm{Tr}\left[{{\Lambda}}_{{E}_{d}}A{{\Lambda}}_{{E}_{d}}A\right]\hfill & \hfill =\hfill & \hfill 4\enspace \mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\enspace \mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\hfill \end{matrix}$.

H.2. Element ${\gamma }_{{V}_{L}{V}_{L}}^{\text{eq}}$

$$\begin{align}\hfill \underset{{k}_{\mathrm{B}}{T}_{0}\to {0}^{+}}{\mathrm{lim}}\enspace {\gamma }_{{V}_{L}{V}_{L}}^{\text{eq}}& =\hfill & \frac{1}{4\pi }\enspace \mathrm{Tr}{\left[{{\Lambda}}_{{V}_{L}}A{{\Lambda}}_{{V}_{L}}A\right]}_{\varepsilon ={\mu }_{0}}\hfill \\ \hfill & =\hfill & \hfill 4\pi \frac{{\left({\mu }_{0}-{E}_{d}\right)}^{2}{V}_{L}^{2}}{{\left\{{V}_{R}^{2}+{V}_{L}^{2}\right\}}^{2}}{N}_{d}\left({\mu }_{0}\right){N}_{d}\left({\mu }_{0}\right)\hfill \\ \hfill & \hfill & +\frac{1}{2\pi }\frac{1}{{V}_{L}^{2}}\tilde {T}\left({\mu }_{0}\right).\hfill \end{align}$$

Here we used

$$\begin{equation*}\begin{matrix}\hfill \mathrm{Tr}\left[{{\Lambda}}_{{V}_{L}}A{{\Lambda}}_{{V}_{L}}A\right]=8\left(\mathrm{I}\mathrm{m}\enspace {G}_{12}^{r}\enspace \mathrm{I}\mathrm{m}\enspace {G}_{12}^{r}+\mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\enspace \mathrm{I}\mathrm{m}\enspace {G}_{11}^{r}\right)\hfill \end{matrix}\end{equation*}$$

and

$$\begin{align*}\hfill \mathrm{I}\mathrm{m}\enspace {G}_{12}^{r}\enspace \mathrm{I}\mathrm{m}\enspace {G}_{12}^{r}+\mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\enspace \mathrm{I}\mathrm{m}\enspace {G}_{11}^{r}& \hfill & \hfill \\ \hfill =2{\pi }^{2}\frac{{\left(\varepsilon -{E}_{d}\right)}^{2}{V}_{L}^{2}}{{\left\{{V}_{R}^{2}+{V}_{L}^{2}\right\}}^{2}}\left(-\frac{1}{\pi }\mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\right)\left(-\frac{1}{\pi }\mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\right)& \hfill & \hfill \\ \hfill +\frac{1}{4{V}_{L}^{2}}4\frac{1}{{\left\vert \varepsilon -{E}_{d}-{{\Sigma}}_{R}-{{\Sigma}}_{L}\right\vert }^{2}}{{\Gamma}}_{R}{{\Gamma}}_{L}.\end{align*}$$

The latter expression requires some algebra but comes from equations (H1)–(H5).

H.3. Element ${\gamma }_{{V}_{R}{V}_{R}}^{\text{eq}}$

$$\begin{align}\hfill \underset{{k}_{\mathrm{B}}{T}_{0}\to {0}^{+}}{\mathrm{lim}}\enspace {\gamma }_{{V}_{R}{V}_{R}}^{\text{eq}}& =\hfill & \frac{1}{4\pi }\enspace \mathrm{Tr}{\left[{{\Lambda}}_{{V}_{R}}A{{\Lambda}}_{{V}_{R}}A\right]}_{\varepsilon ={\mu }_{0}}\hfill \\ \hfill & =\hfill & \hfill 4\pi \frac{{\left({\mu }_{0}-{E}_{d}\right)}^{2}{V}_{R}^{2}}{{\left\{{V}_{R}^{2}+{V}_{L}^{2}\right\}}^{2}}{N}_{d}\left({\mu }_{0}\right){N}_{d}\left({\mu }_{0}\right)\hfill \\ \hfill & \hfill & +\frac{1}{2\pi }\frac{1}{{V}_{R}^{2}}\tilde {T}\left({\mu }_{0}\right).\hfill \end{align}$$

Due to the symmetry of the problem, this element is equal to ${\gamma }_{{V}_{L}{V}_{L}}^{\text{eq}}$ but replacing V_L by V_R in the formulas.

H.4. Element ${\gamma }_{{E}_{p}{V}_{L}}^{\text{eq}}$

$$\begin{align}\hfill \underset{{k}_{\mathrm{B}}{T}_{0}\to {0}^{+}}{\mathrm{lim}}\enspace {\gamma }_{{E}_{d}{V}_{L}}^{\text{eq}}& =\hfill & \frac{1}{4\pi }\enspace \mathrm{Tr}{\left[{{\Lambda}}_{{E}_{d}}A{{\Lambda}}_{{V}_{L}}A\right]}_{\varepsilon ={\mu }_{0}}\hfill \\ \hfill & =\hfill & \hfill 2\pi \frac{\left(\mu -{E}_{d}\right){V}_{L}}{{V}_{L}^{2}+{V}_{R}^{2}}{N}_{d}\left({\mu }_{0}\right){N}_{d}\left({\mu }_{0}\right),\hfill \end{align}$$

where we used $\begin{matrix}\hfill \mathrm{Tr}\left[{{\Lambda}}_{{E}_{d}}A{{\Lambda}}_{{V}_{L}}A\right]\hfill & \hfill =\hfill & \hfill -8\enspace \mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\enspace \mathrm{I}\mathrm{m}\enspace {G}_{12}^{r}\hfill \end{matrix}$.

H.5. Element ${\gamma }_{{E}_{d}{V}_{R}}^{\text{eq}}$

$$\begin{align}\hfill \underset{{k}_{\mathrm{B}}{T}_{0}\to {0}^{+}}{\mathrm{lim}}\enspace {\gamma }_{{E}_{d}{V}_{R}}^{\text{eq}}& =\hfill & \hfill \frac{1}{4\pi }\enspace \mathrm{Tr}{\left[{{\Lambda}}_{{E}_{P}}A{{\Lambda}}_{{V}_{R}}A\right]}_{\varepsilon ={\mu }_{0}}\hfill \\ \hfill & =\hfill & \hfill 2\pi \frac{\left(\mu -{E}_{d}\right){V}_{R}}{{V}_{L}^{2}+{V}_{R}^{2}}{N}_{d}\left({\mu }_{0}\right){N}_{d}\left({\mu }_{0}\right).\hfill \end{align}$$

Due to the symmetry of the problem, this element is equal to ${\gamma }_{{E}_{p}{V}_{L}}^{\text{eq}}$ but replacing V_L by V_R in the formulas.

H.6. Element ${\gamma }_{{V}_{L}{V}_{R}}^{\text{eq}}$

$$\begin{align}\hfill \underset{{k}_{\mathrm{B}}{T}_{0}\to {0}^{+}}{\mathrm{lim}}\enspace {\gamma }_{{V}_{L}{V}_{R}}^{\text{eq}}& =\hfill & \frac{1}{4\pi }\enspace \mathrm{Tr}{\left[{{\Lambda}}_{{V}_{L}}A{{\Lambda}}_{{V}_{R}}A\right]}_{\varepsilon ={\mu }_{0}}\hfill \\ \hfill \hfill & \hfill =\hfill & 4\pi \frac{{\left({\mu }_{0}-{E}_{d}\right)}^{2}{V}_{L}{V}_{R}}{{\left\{{V}_{R}^{2}+{V}_{L}^{2}\right\}}^{2}}{N}_{d}\left({\mu }_{0}\right){N}_{d}\left({\mu }_{0}\right)\hfill \\ \hfill & \hfill & -\frac{1}{2\pi {V}_{L}}\frac{1}{{V}_{R}}\tilde {T}\left({\mu }_{0}\right).\hfill \end{align}$$

Here we used

$$\begin{equation*}\begin{matrix}\hfill \mathrm{Tr}\left[{{\Lambda}}_{{V}_{L}}A{{\Lambda}}_{{V}_{R}}A\right]=8\left(\mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\enspace \mathrm{I}\mathrm{m}\enspace {G}_{13}^{r}+\mathrm{I}\mathrm{m}\enspace {G}_{23}^{r}\enspace \mathrm{I}\mathrm{m}\enspace {G}_{12}^{r}\right)\hfill \end{matrix}\end{equation*}$$

and

$$\begin{align*}\hfill \mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\enspace \mathrm{I}\mathrm{m}\enspace {G}_{13}^{r}+\mathrm{I}\mathrm{m}\enspace {G}_{23}^{r}\enspace \mathrm{I}\mathrm{m}\enspace {G}_{12}^{r}& \hfill & \hfill \\ \hfill =2{\pi }^{2}\frac{{\left(\varepsilon -{E}_{d}\right)}^{2}{V}_{L}{V}_{R}}{{\left\{{V}_{R}^{2}+{V}_{L}^{2}\right\}}^{2}}\left(-\frac{1}{\pi }\mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\right)\left(-\frac{1}{\pi }\mathrm{I}\mathrm{m}\enspace {G}_{22}^{r}\right)& \hfill & \hfill \\ \hfill -\frac{1}{4{V}_{L}}\frac{1}{{V}_{R}}4\frac{1}{{\left\vert \varepsilon -{E}_{d}-{{\Sigma}}_{R}-{{\Sigma}}_{L}\right\vert }^{2}}{{\Gamma}}_{L}{{\Gamma}}_{R}.\end{align*}$$

The latter expression requires some algebra but comes from equations (H1)–(H5).

In this section, we will analyze a concrete tight-binding example to help the readers understand how to use our results in physical systems. Our example consists of an atomic wire with a bend where the corner atom is free to move in the x–y plane. This example was already studied in reference [26] where the authors found that nonconservative CIFs drive the movement of the corner atom in a periodic orbit. This kind of system can be considered as an example of an adiabatic quantum motor [27]. We will not enter into the experimental details of how to realize such a proposal as this is beyond the scope of the present work. However, a good candidate to further study the proposal could be a quantum dot mounted on a nanopillar and coupled to two leads forming an angle. Another possibility would be to replace the nanopillar with some appropriate organic molecule.

The electronic Hamiltonian of the example is the same as that shown in equation (4) and figure 1. The fact that leads form an angle does not change the electronic Hamiltonian. For simplicity, we will assume the angle α between the leads is 90°, where the L-lead goes along the 'y' axis while the R-lead goes along the 'x' axis. Typically the hoppings vary exponentially with the distance between atoms, but here we will take the linearized version of this, such as

$$\begin{align}\hfill {V}_{L}& ={V}_{0}+{V}^{\prime }\delta y\hfill \\ \hfill {V}_{R}& ={V}_{0}+{V}^{\prime }\delta x,\hfill \end{align} \tag{ I1 }$$

where δx and δy are the deviations from the equilibrium position of the dot along the 'x' and 'y' axis respectively. Note that the above equation can be readily modify to take into account α different from 90°. Just to put some numbers, we will take the same tight-binding parameters as in reference [26]: E₀ = 0 eV, E_d = 2 eV, V₀ = −3.88 eV, ${V}^{\prime }=4{\times}3.88/2.5\enspace \mathrm{e}\mathrm{V}$ Å⁻¹, and ɛ_F = E₀ (the leads' bands are half-filled). This corresponds to the renormalized parameters _d = 2/3.88, _F = 0, and ${v}_{L}^{\left(\delta y=0\right)}={v}_{R}^{\left(\delta x=0\right)}=1$. According to equation (33), the CIFs F_x and F_y acting on the x and y directions are related to the dimensionless forces ${F}_{{V}_{L}}$ and ${F}_{{V}_{R}}$ by

$$\begin{align}\hfill {F}_{x}& ={F}_{{V}_{R}}{V}^{\prime }\hfill \\ \hfill {F}_{y}& ={F}_{{V}_{L}}{V}^{\prime }.\hfill \end{align} \tag{ I2 }$$

For small displacement from the equilibrium position, the work done by CIFs in a periodic orbit of the dot's position can be written as

$$\begin{equation}W={\left(\nabla {\times}\mathbf{F}\right)}_{z}\iint \mathrm{d}x\enspace \mathrm{d}y,\end{equation} \tag{ I3 }$$

where ∬dx dy is the area A enclosed by the movement of x and y (a figure-of-eight orbit gives a zero area). As highlighted in appendix E, W is independent of how the variables of the Hamiltonian (V_L and V_R ) depends on the movement of mechanical variables (δx and δy). Thus, assuming E_d is constant along the trajectory and taking the small temperature limit and for small bias voltages, the work yields

$$\begin{align}\hfill W& ={\left(\nabla {\times}\mathbf{F}\right)}_{{E}_{d}}\iint \mathrm{d}{V}_{L}\enspace \mathrm{d}{V}_{R}\hfill \\ \hfill & ={\left(\nabla {\times}\mathbf{F}\right)}_{{E}_{d}}{\left({V}^{\prime }\right)}^{2}A\hfill \\ \hfill & =\delta \mu {W}_{{E}_{d}}^{\left(\delta \mu \right)}{\left({V}^{\prime }\right)}^{2}A\hfill \\ \hfill & =\tilde {T}\left({{\epsilon}}_{F}\right){N}_{d}\left({{\epsilon}}_{F}\right)2\frac{\left({{\epsilon}}_{F}-{{\epsilon}}_{d}\right)}{{v}_{L}{v}_{R}}\frac{\delta \mu }{{V}_{0}}{\left({V}^{\prime }\right)}^{2}A.\hfill \end{align} \tag{ I4 }$$

For the chosen parameters, we numerically found the optimal working point (values of v_L and v_R that maximizes ${W}_{{E}_{d}}^{\left(\delta \mu \right)}$ for _F = 0). These values are approximately v_L = v_R = 0.51, which corresponds to $\delta x=\delta y=\left(0.51-1\right){V}_{0}/{V}^{\prime }=-0.306\;\mathring{A}$.

Taking _F ≈ |2|, which dramatically increases ${W}_{{E}_{d}}^{\left(\delta \mu \right)}$ according to the discussion of section 3, the optimal working point is near the virtual-localized QDPT, see equation (12). This occurs for

$$\begin{align}\hfill {v}_{L}^{2}+{v}_{R}^{2}& =2\mp {{\epsilon}}_{d}\hfill \\ \hfill {\left(1+1.6\delta x\right)}^{2}+{\left(1+1.6\delta y\right)}^{2}& =2\mp 0.515,\hfill \end{align} \tag{ I5 }$$

where the minus sign corresponds to _F ≈ −2 (Fermi energy close to the lower band edge) while the plus sign corresponds to _F ≈ 2 (Fermi energy close to the upper band edge).