Second-order adiabatic expansions of heat and charge currents with nonequilibrium Green's functions
Sebastián E. Deghi Raúl A. Bustos-Marún





^1Instituto de Física Enrique Gaviola and Facultad de
Matemática, Astronomía, Física y Computacin, Universidad Nacional de
Córdoba, Ciudad Universitaria, Córdoba, 5000, Argentina

^2Facultad de Ciencias Químicas, Universidad Nacional de Córdoba, Ciudad
Universitaria, Córdoba, 5000, Argentina



Due to technological needs, nanoscale heat management, energy conversion and quantum thermodynamics have become key areas of research, putting heat pumps and nanomotors center stage.
The treatment of these particular systems often requires the use of adiabatic expansions in terms of the frequency of the external driving or the velocity of some classical degree of freedom.
However, due to the difficulty of getting the expressions, most works have only explored first-order terms. Despite this, adiabatic expansions have allowed the study of intriguing phenomena such as adiabatic quantum pumps and motors, or electronic friction.
Here, we use nonequilibrium Green's functions, within a Schwinger-Keldysh approach, to develop second-order expressions for the energy, heat, and charge currents.
We illustrate, through two simple models, how the obtained formulas produce physically consistent results, and allow for the thermodynamic study of unexplored phenomena, such as second-order monoparametric pumping.

Second-order adiabatic expansions of heat and charge currents with
nonequilibrium Green's functions.

    Sebastián E. Deghi,^1 and Raúl A. Bustos-Marún^1,2
    April 29, 2024
======================================================================================================





§ INTRODUCTION


The conversion of heat into work (either mechanical or electrical)
has been at the center of technological and scientific interest since
the first studies of steam engines <cit.>. In the past
decades, miniaturization has guided this interest toward micro- and
nanoscale heat-to-work conversion <cit.>.
This was not only due to the fundamental questions posed by the new
size scale <cit.>, but also to practical
reasons such as the need for efficient heat management of nanodevices.
However, the nanoscale presents new challenges. On the one hand, the
regions of interest now involved a finite number of particles, as
opposed to classical thermodynamics. On the other hand, the explicit
treatment of quantum mechanical effects often becomes unavoidable.

Importantly, nanoscale quantum machines (devices with strong quantum
effects) require appropriate and efficient computational treatments.
This is so, since processes with vastly different time scales may
coexist there. Therefore, without any type of approximation, the calculations
must then be carried out with a temporal resolution given by the fastest
degrees of freedom (DOFs) but with total times determined by the slowest
DOFs. This clearly makes numerical calculations highly inefficient.
Even worse, sometimes contrasting theoretical approaches may be needed
for the different DOFs.

A successful strategy followed in the past consisted of treating some
slow DOFs classically (typically nuclear, mechanical, or external
ones). In contrast, the fast DOFs are treated fully quantum (typically
electrons or optical phonons). For nanoelectronic or nanoelectromechanical
devices, the time scale separation also drives the establishment of
steady-state currents which parametrically depend on the position
and velocity of the slow DOFs. This situation is amenable to some
kind of adiabatic approximation where the observable of interest is
expanded in terms of the frequency of the external driving or the
velocity of the classical DOFs. In this context, different theoretical
approaches have been used including non-equilibrium Green’s
functions (NEGF) <cit.>,
real-time diagrammatic approaches <cit.>
(or similar approaches based on the adiabatic expansions of the system's
reduced density matrix <cit.>), scattering formalism
<cit.>, DFT-based calculations
<cit.>, and hierarchical equation of motion approaches <cit.>.

Within the context of quantum transport, adiabatic approximations
up to the first order have allowed the treatment of fascinating phenomena
such as adiabatic quantum pumping <cit.>, adiabatic
quantum motors<cit.>, negative friction coefficients
<cit.>, reciprocity breaking <cit.>,
hysteresis <cit.>, current-noise induced
by thermal oscillations <cit.>, etc. Regardless of the
notable advances made, one may wonder about the unexplored phenomena
that may await beyond first-order adiabatic treatments. In this sense,
Kershaw et al. <cit.> pioneered the second-order
treatment of electric currents within a NEGF approach. However, a
thermodynamically consistent second-order treatment of quantum machines
necessarily requires addressing heat currents. Moreover, explicit
expressions of the observables are desirable for numerical calculations,
instead of implicit expressions that need to be developed for the
cases of interest. In the present manuscript, we present the complete
and explicit expressions, up to the second order of the adiabatic
expansions of heat and charge currents. Our results are based on NEGFs
within a Schwinger-Keldysh approach. The obtained formulas are general
and restricted only by the assumption of time-independent self-energies.[It is important to highlight that this quite common assumption is
not a limitation for a wide class of systems. For example, in a tight-binding
approach, one can redefine the local system by adding the sites of
the leads with time dependence <cit.>.]

This manuscript is organized as follows. Sec. <ref>
presents the general theory based on the NEGF formalism. Having developed
the theoretical tools, in Sec. <ref>, we provide
the second-order adiabatic expansion for charge, energy, and heat
currents. Sec. <ref> contains two different models illustrating
how the corrections for each order work and the usefulness of the
expressions for thermodynamic analysis. In particular, the last example
treats second-order quantum pumping, a phenomenon not previously described
up to our knowledge. The examples also include a numerical verification
of our formulas through the first law of thermodynamics. Finally,
Sec. <ref> comprises a summary and a brief discussion.
In an effort not to overload the readers, all the nonessential comments
and mathematical details lay in the Appendices from <ref>
to <ref>.



§ GENERAL THEORY




 §.§ Generic Hamiltonian


To define the different parts of the system and their associated creation
(annihilation) operators, we introduce the total Hamiltonian H
that describes the total system. It consists of a core region
S connected to several macroscopic reservoirs α (the leads).
The core region, which typically has nanometric dimensions, will be
referred to as the local system. The total Hamiltonian takes
the form

    H     =α∑H_α+H_S+α∑V_α,

where H_S, in general, corresponds to a time-dependent local
system, H_α describes the time-independent α-lead,
and the term V_α contains the time-independent interaction
between the local system and the α-lead. The Hamiltonian of
the local system depends parametrically on time through a multidimensional
vector X, which describes the mechanical degrees
of freedom. The local system's Hamiltonian assumes the following form
when the second quantization is applied

    H_S    =l,s∑h_l,s(X)d_l^†d_s.

In the notation we have used, the operators d_l^† and
d_l create or annihilate, respectively, an electron within the
local system. The Hamiltonians of the leads can be written as

    H_α    =k∑ϵ_α kc_α k^†c_α k,

where ϵ_α k is the energy of the α-lead in
the state k, c_α k^† creates an electron in the
α-lead with the state k, whereas the c_α k annihilates
it. For simplicity, the index k also includes the electron's spin.
Finally, the tunneling interaction describes the coupling between
the local system and the α-leads

    V_α    =k,l∑(t_α k,lc_α k^†d_l+t_α k,l^*d_l^†c_α k),

where t_α k,l are the tunneling amplitudes between the leads
and the local system. 



 §.§ Dynamics of non-equilibrium open quantum systems


Within the context of the Schwinger-Keldysh approach for the NEGF
formalism <cit.>
we define the elements of the retarded and advanced Green’s
functions, 𝒢^R and 𝒢^A
respectively. For a system that evolves, not necessarily in an equilibrium
process, from t' to t times, they
are

    𝒢_l,s^R(t,t')     =-i/ħΘ(t-t')⟨{ d_s(t),d_l^†(t')}⟩ ,

and

    𝒢_s,l^A(t',t)     =[𝒢_l,s^R(t,t')]^*,

where {∙,∙} is the anticommutator and
⟨∙⟩ is the quantum expectation value.

For the states of the local system, the lesser Green’s
function elements are given by

    𝒢_l,s^<(t,t')     =i/ħ⟨ d_l^†(t')d_s(t)⟩ ,

whereas for the states that propagate between the local system and
leads, the lesser Green’s function elements are

    𝒢_α k,l^<(t,t')     =i/ħ⟨ c_α k^†(t')d_l(t)⟩ .


Furthermore, the retarded and advanced self-energies are assumed to
be stationary.[For our purposes, stationary implies that the self-energies only depend
on the time difference t-t', which means that both V_α and H_α are independent of the mechanical
degrees of freedom.] The elements of the retarded self-energy for the α-lead take
the form

    Σ_α,l,s^R(t,t')     =k∑t_l,α kg_α k^R(t,t')t_α k,s,

where

    g_α k^R(t,t')     =-i/ħΘ(t-t')e^-i/ħϵ_α k(t-t').

Hence, the elements of total retarded self-energy comprise the sum
of all leads

    Σ_l,s^R(t,t')     =α∑Σ_α,l,s^R(t,t').

To obtain the advanced self-energy, we calculate the adjoint of the
retarded one

    Σ_α,s,l^A(t',t)     =[Σ_α,l,s^R(t,t')]^*.

Likewise, the elements of the lesser self-energy of the α-lead
are given by

    Σ_α,l,s^<(t,t')     =k∑t_l,α kg_α k^<(t,t')t_α k,s,

where

    g_α k^<(t,t')     =2π i/ħf_α(ϵ_α k)e^-i/ħϵ_α k(t-t').

Here, f_α is the α-lead's Fermi-Dirac distribution.
To conclude, the elements of the total lesser self-energy are

    Σ_l,s^<(t,t')     =α∑Σ_α,l,s^<(t,t').


Given the essential operators, the next step is to delve into the
quantum dynamics, which in Schwinger-Keldysh formalism is expressed
by the integrodifferential Dyson equation in terms of the retarded
Green’s function 

    -iħ∂_t'𝒢^R(t,t')     =δ(t-t')+𝒢^R(t,t')H_S(t')
        +∫𝒢^R(t,t_1)Σ^R(t_1,t')dt_1.

Here, the Hamiltonian of the local system, retarded Green’s
function and self-energy are given by Eqs. (<ref>),
(<ref>) and (<ref>),
respectively.

At this point, we can apply many techniques to solve Eq. (<ref>).
One way is to use numerical methods <cit.>. Another one
is to develop an adiabatic expansion of the observables. To carry
out the latter technique, we must decompose the defined operators
and Dyson’s equation into two-time scales, which is
done through the use of a Wigner transform.



 §.§ Adiabatic expansion for retarded Green’s function


The idea is to separate the fast microscopic dynamics, followed by
the electrons, from the slow macroscopic changes, determined by the
classical mechanical degrees of freedom. Once the time scales have
been distinguished, a gradient expansion will be carried out on the
slow variables. First of all, we have to define the following transform
of time variables

    T     =t+t'/2,
    τ    =t-t'.

The time scale T is called the slow variable, while τ is
the fast variable. The Wigner transform for the retarded Green’s
function is defined as

    𝒢^R(T,ε)     =∫𝒢^R(t,t')e^i/ħετdτ,

whereas its inverse transformation is

    𝒢^R(t,t')     =1/2πħ∫𝒢^R(T,ε)e^-i/ħετdε.

For the elements of the retarded self-energy described in Eq. (<ref>),
Σ^R(t,t')≡Σ^R(τ), the
Wigner transform becomes a Fourier transform, where


    Σ^R(ε)     =    ∫Σ^R(τ)e^i/ħετdτ,

and


    Σ_α,l,s^R(ε)     =Δ_α,l,s(ε)-iΓ_α,l,s(ε).

Here, the level-width functions Γ_α,l,s are given by[The definitions of the level-width functions may differ from that
of other authors <cit.>. Here, we followed
the convention used in <cit.>.]

    Γ_α,l,s(ε)     =πk∑δ(ε-ϵ_α k)t_l,α kt_α k,s.

The level-shift functions Δ_α,l,s can be calculated
from Γ_α,l,s using the Kramers-Kroning relation <cit.>.
After applying the adjoint operator to Eq. (<ref>)
we get the advanced self-energy Σ_α,l,s^A. For the
elements of the lesser self-energy given in Eq. (<ref>),
their Wigner transform give

    Σ_α,l,s^<(ε)     =2if_α(ε)Γ_α,l,s(ε).

Given the operators in the energy domain, we must apply the Moyal
product to Eq. (<ref>),
yielding (see Appendix <ref>)

    1     =-1/2iħ∂_T𝒢^R(T,ε)+ε𝒢^R(T,ε)
        -𝒢^R(T,ε)e^iħ/2(∂_ε∂_T-∂_T∂_ε)Σ^R(ε)
        -𝒢^R(T,ε)e^iħ/2(∂_ε∂_T-∂_T∂_ε)Ĥ_S(T),

The above expression implies that the unknown retarded green's function
𝒢^R can be expressed as an infinite sum
over its derivatives, which are also unknown in principle. To solve
this, we have to implement an iterative approach (see Appendix <ref>)
which starts at zero-order by approximating 𝒢^R
by the adiabatic Green's function G^R, defined as

    G^R    =[ε I-H_S(T)-Σ^R(ε)]^-1.

The above adiabatic (or frozen) Green's function G^R is formally
equal to 𝒢^R(T,ε)
in the limit of infinitely slow mechanical degrees of freedom, which
constitute a form of the Born-Oppenheimer approximation.

By the chain rule, the slow-time T derivatives of the Hamiltonian
of the local system, up to the second order, are

    ∂_TH_S    =[ν=1]M∑_νẊ_ν,
    ∂_T^2H_S    =[μ,ν=1]M∑_μνẊ_μẊ_ν+[ν=1]M∑_νẌ_ν.

Here, Ẋ_ν and Ẍ_ν are the time derivatives
of the mechanical degree of freedom ν, ∂_TX_ν
and ∂_T^2X_ν respectively. Additionally, the elements
of the matrices _ν and _μν are
given by

    [_ν]_l,s    =∂ h_l,s/∂ X_ν,    [_μν]_l,s    =∂^2h_l,s/∂ X_μ∂ X_ν.

Then, applying the iterative method twice to second-order terms, the
resulting retarded Green's function takes the form

    𝒢^R   ≃ G^R(T,ε)+iħ/2[ν=1]M∑Ω_1,ν^R(T,ε)Ẋ_ν
        +(iħ/2)^2{[μ,ν=1]M∑Ω_11,μν^RẊ_μẊ_ν+[ν=1]M∑Ω_2,ν^RẌ_ν} .

To make the notation more compact, we have defined the following operators

    Ω_1,ν^R    =Ξ_1[G^R,_ν,G^R],
    Ω_11,μν^R    =Ξ_11[G^R,K_μν^R,G^R],
    Ω_2,ν^R    =1/2Ξ_2[G^R,_ν,G^R].

where we have also introduced the definitions

    K_μν^R    =_μG^R(T,ε)_ν+1/2_μν,

and

    Ξ_1[A,B,C]     =(∂_εA)BC-AB(∂_εC),
    Ξ_2[A,B,C]     =(∂_ε^2A)BC-2(∂_εA)B(∂_εC)
        +AB(∂_ε^2C),
    Ξ_11[A,B,C]     =Ξ_2[A,B,C]-A(∂_ε^2B)C.


Take into account that the order to which the Planck's constant is
elevated is consistent with the order of the adiabatic expansion.
In this sense, it can be used as a book-keeping
parameter to keep track of the expansion order for a given term.

Finally, to calculate the advanced Green function, we use the relation

    𝒢^A    =[𝒢^R(T,ε)]^†.





 §.§ Adiabatic expansion for lesser Green’s function


To find the Wigner’s transform of the lesser Green’s
functions, Eqs. (<ref>)
and (<ref>),
we start with the relation <cit.>

    𝒢^<(t,t')     =∬𝒢^R(t,t_1)Σ^<(t_1,t_2)𝒢^A(t_2,t')dt_1dt_2.

After applying the Wigner transform to both sides of the above formula
and, afterwards, the Moyal product, we arrive at the expression

    𝒢^<    =𝒢^R(T,ε)e^iħ/2(∂_ε∂_T-∂_T∂_ε)Σ^<(ε)×
        e^iħ/2(∂_ε∂_T-∂_T∂_ε)𝒢^A(T,ε).

As we discussed for the retarded Green’s function,
𝒢^< takes the form of a set of infinite
sums (see Appendix <ref>),
where the elements of Σ^< are given in Eq. (<ref>).
To solve this, we first insert in Eq. (<ref>)
the second-order expansion of the operators 𝒢^R
and 𝒢^A derived before, and then we truncate
the series at second order. This yields

    𝒢^<   ≃ G^<(T,ε)+iħ/2[ν=1]M∑Ω_1,ν^<(T,ε)Ẋ_ν
        +(iħ/2)^2{[μ,ν=1]M∑Ω_11,μν^<Ẋ_μẊ_ν+[ν=1]M∑Ω_2,ν^<Ẍ_ν} ,

Where the adiabatic lesser Green function is

    G^<    =G^R(T,ε)Σ^<(ε)G^A(T,ε),

and the remaining operators are defined by

    Ω_1,ν^<    =Ξ_1[G^R,_ν,G^<]+Ξ_1[G^<,_ν,G^A],
    Ω_11,μν^<    =Ξ_11[G^R,K_μν^R,G^<]+Ξ_11[G^<,K_μν^A,G^A]
        +Ξ_11[G^R,_μG^<_ν,G^A],
    Ω_2,ν^<    =1/2{Ξ_2[G^<,_ν,G^A]+Ξ_2[G^R,_ν,G^<]} .

In the above expressions we have introduced the following term to
compact the notation

    K_μν^A    =_μG^A(T,ε)_ν+1/2_μν.





§ OBSERVABLES


We have already developed explicit formulas for evaluating up-to-second-order
corrections to the adiabatic Green's functions. In the following section,
we will apply these results to obtain close expressions for the second-order
adiabatic corrections of three observables of interest within quantum
transport: charge, energy, and heat currents.



 §.§ Charge current


For an α-lead, the average charge current through it can be
written as the mean value of the time derivative of the number operator
N_α

    I_α    =-e⟨Ṅ_α(t)⟩ ,

where N_α=k∑c_α k^†c_α k, and
e is the modulus of the electron charge. Then, applying Ehrenfest’s
theorem, the commutation relations for fermions, and using Eqs. (<ref>)
and (<ref>)
in Eq. (<ref>), we get a formula
for the charge current in terms of the lesser Green's function

    I_α    =2eRe{k,l∑t_α k,l𝒢_l,α k^<(t,t)} .


To evaluate the time-dependent current I_α, we must proceed
in accord with the Langreth’s rules <cit.>,
yielding

    I_α    =2e∫Re{tr[𝒢^R(t,t')_α^<(t',t)..
       ..+𝒢^<(t,t')_α^A(t',t)]} dt'.


The next step is to insert the Wigner transform of the Green's functions
into Eq. (<ref>)
and make a gradient expansion, the resulting expression reads (see
Appendix <ref>)

    I_α    =2e/h[N=0]∞∑(-1)^N/N!∫Re{tr{ C_α,N^I}} dε,

where

    C_α,N^I    =(iħ/2)^N(∂_T^N𝒢^R∂_ε^N_α^<+∂_T^N𝒢^<∂_ε^N_α^A).


It is important to highlight that, although this expression is exact,
its explicit evaluation necessarily requires approximations. For this
purpose, we will use the up-to-second-order adiabatic expansion
of the Green's functions found previously. After some algebra, one
can find the up-to-second-order adiabatic expansion of the charge
current, giving

    I_α   ≃ I_α^(0)+I_α^(1)+I_α^(1,1)+I_α^(2).

The first two terms of the above expression are well known. The former,
I_α^(0), is equivalent to Landauer's formula,
while the second one, I_α^(1), is the quantum
pumping contribution to the charge current <cit.>

    I_α^(0)    =-e𝒩_α^(0),
    
    I_α^(1)    =-eν∑𝒩_α,ν^(1)Ẋ_ν.

Here, we have introduced the terms 

    𝒩_α^(0)    =∫𝒦_α^(0)(T,ε)dε,    𝒩_α,ν^(1)    =∫𝒦_α,ν^(1)(T,ε)dε.

We acquire the preceding definition for practical purposes, which
will be useful later, being the kernel’s integrals
defined as

    𝒦_α^(0)    =-2/hRe{tr(Φ_α^(0))} ,
    𝒦_α,ν^(1)    =1/πIm{tr(∂_εG^R_νΦ_α^(0)+∂_εG^<_νΦ_α^(1))} ,

where

    Φ_α^(0)    =G^r_α^<(ε)+G^<_α^A(ε),    Φ_α^(1)    =G^A_α^A(ε).


The remaining terms in Eq. (<ref>)
constitute the second-order adiabatic terms of the charge current.
These formulas can be written as

    I_α^(2)    =-eν∑𝒩_α,ν^(2)Ẍ_ν,
    
    I_α^(1,1)    =-eμ,ν∑𝒩_α,μν^(1,1)Ẋ_μẊ_ν.

The current derivatives 𝒩_α,ν^(2)=-1/e∂ I_α/∂Ẍ_ν
and 𝒩_α,μν^(1,1)=-1/e∂^2I_α/∂Ẋ_μ∂Ẋ_ν
are

    𝒩_α,ν^(2)    =∫𝒦_α,ν^(2)(T,ε)dε,    𝒩_α,μν^(1,1)    =∫𝒦_α,μν^(1,1)(T,ε)dε,

The above kernels of integrals are

    𝒦_α,ν^(2)    =ħ/2πRe{tr{∂_ε^2G^R_νΦ_α^(0)+∂_ε^2G^<_νΦ_α^(1)}} ,
    𝒦_α,μν^(1,1)    =ħ/πRe{tr{Φ_μν^(2)Φ_α^(0)+Φ_μν^(3)Φ_α^(1)}} ,

where

    Φ_μν^(2)    =∂_ε(∂_εG^RK_μν^R),
    Φ_μν^(3)    =∂_ε(∂_εG^R_μG^<_ν+∂_εG^<K_μν^A).





 §.§ Energy current


Analogous to the charge current, the energy current flowing in the
α-lead Ė_α is defined as the mean value of
the time derivative of the Hamiltonian of the α-lead H_α,

    Ė_α    =⟨Ḣ_α(t)⟩ .

As we proceeded in the previous section for the charge current, we
will apply Ehrenfest's theorem, the commutation relations for fermions,
and Eqs. (<ref>)
and (<ref>)
into Eq. (<ref>). This gives

    Ė_α    =-2Re{k,l∑ϵ_α kt_α k,l𝒢_l,α k^<(t,t)} .

Note that the essential difference between Eqs. (<ref>)
and (<ref>)
are the energy weights. This suggests employing a method akin to the
charge currents strategy. However, despite using the same assumption,
there are some subtle deviations from the reasoning for load currents,
as illustrated in Appendix <ref>.
The final formula is

    Ė_α    =-2Re{ iħ∫tr[𝒢^R(t,t')∂_t'_α^<(t',t)..
       ..+𝒢^<(t,t')∂_t'_α^A(t',t)]} dt'.

Unlike Eq. (<ref>),
the energy current formula shown above contains time derivatives over
self-energies. Up to this point, the expression for energy current
seems to be similar to the charge current. By applying the Wigner
transform and the gradient expansion to Eq. (<ref>)
we arrive at (see Appendix <ref>)

    Ė_α    =-2/h[N=0]∞∑(-1)^N/N!∫εRe{tr{ C_α,N^I}} dε
        +2/h[N=0]∞∑(-1)^N/N!∫Re{tr{ C_α,N^E}} dε.


The additional term C_α,N^E is given by[Each of the C_α,N^I functions are dimensionless, whereas
the C_α,N^E functions have energy units.]

    C_α,N^E    =(iħ/2)^N+1(∂_T^N+1𝒢^R∂_ε^N_α^<+∂_T^N+1𝒢^<∂_ε^N_α^A).


After introducing the second-order adiabatic expansion of the retarded
and lesser Green’s functions into Eq. (<ref>),
we get the adiabatic expansion up to the second-order of the energy current

    Ė_α   ≃Ė_α^(0)+Ė_α^(1)+Ė_α^(1,1)+Ė_α^(2).


The adiabatic energy current is 

    Ė_α^(0)    =ℰ̇_α^(0)=∫ε𝒦_α^(0)(T,ε)dε,

and its first correction Ė_α^(1) takes
the form

    Ė_α^(1)    =ν∑ℰ̇_α,ν^(1)Ẋ_ν,

where

    ℰ̇_α,ν^(1)    =∫ε𝒦_α,ν^(1)(T,ε)dε.


These terms are equivalent to those found previously by other authors,
see for example ludovico2016Feb. However, the last two terms
of Eq. (<ref>)
represent original expressions, to our knowledge. They read as

    Ė_α^(2)    =ν∑ℰ̇_α,ν^(2)Ẍ_ν,
    Ė_α^(1,1)    =μ,ν∑ℰ̇_α,μν^(1,1)Ẋ_μẊ_ν,

where the energy fluxes take the form

    ℰ̇_α,ν^(2)    =∫ε𝒦_α,ν^(2)(T,ε)dε,
    ℰ̇_α,μν^(1,1)    =∫ε𝒦_α,μν^(1,1)(T,ε)dε.


Note that, by comparing the adiabatic expansion of the charge current
Eq. (<ref>)
with the energy one (Eq. <ref>),
we recover, order-by-order, the straightforward interpretation of
the energy current. It comes from the integration of all particles
entering and leaving a given lead, but multiplied by the energy each
particle carries. Mathematically this is a consequence of a subtle
cancellation of the second term of Eq. (<ref>).



 §.§ Heat current


The heat current flowing through the α-lead is the last observable
we will study up to the second order. We can determine the expectation
value of heat current based on the first law of thermodynamics <cit.>.
The total energy gained or lost by a lead α is the sum of
the heat and the work done by the particles being exchanged with the
local system. In this way, the heat current J_α is given
by

    J_α    =Ė_α-μ_α⟨Ṅ_α(t)⟩ .


The formal definition of Eq. (<ref>)
contains two mean values that we developed in subsections <ref>
and <ref>. Using the previous results, we can
straightforwardly write the adiabatic expansion, up to the second order,
of the heat current, just as we have done before

    J_α   ≃ J_α^(0)+J_α^(1)+J_α^(1,1)+J_α^(2).


Employing Eqs. (<ref>-<ref>)
and (<ref>-<ref>),
we found

    J_α^(0)    =𝒥_α^(0),     J_α^(2)    =ν∑𝒥_α,ν^(2)Ẍ_ν,
    
    J_α^(1)    =ν∑𝒥_α,ν^(1)Ẋ_ν,     J_α^(1,1)    =μ,ν∑𝒥_α,μν^(1,1)Ẋ_μẊ_ν,

where

    𝒥_α^(0)    =ℰ̇_α^(0)-μ_α𝒩_α^(0),    𝒥_α,ν^(2)    =ℰ̇_α,ν^(2)-μ_α𝒩_α,ν^(2),
    𝒥_α,ν^(1)    =ℰ̇_α,ν^(1)-μ_α𝒩_α,ν^(1),    𝒥_α,μν^(1,1)    =ℰ̇_α,μν^(1,1)-μ_α𝒩_α,μν^(1,1).


Likewise, the above results can be written as integrals, giving

    𝒥_α^(0)    =∫(ε-μ_α)𝒦_α^(0)(T,ε)dε,
    𝒥_α,ν^(1)    =∫(ε-μ_α)𝒦_α,ν^(1)(T,ε)dε,
    𝒥_α,ν^(2)    =∫(ε-μ_α)𝒦_α,ν^(2)(T,ε)dε,
    𝒥_α,μν^(1,1)    =∫(ε-μ_α)𝒦_α,μν^(1,1)(T,ε)dε.


Thus far, we have offered a complete set of quantum transport observables
to understand how the second-order corrections work. These results
became a starting point to explore thermodynamic properties from a
quantum point of view.



§ MODELS


In the previous sections, we formulated the theory of adiabatic expansions
up to the second order for a set of fundamental observables. Now, we will
evaluate the charge, energy, and heat currents for two different time-dependent
devices. In studying these examples, we have three goals in mind.
First, we want to illustrate the usage of the developed formulas and
their utility for the thermodynamic analysis of nanomachines. Second,
we want to test their validity by showing that the results they produce
are physically reasonable. In particular, we will show that energy
and particle number are conserved at all times and for all parameters
in the order-by-order expansions <cit.>. This involves
the comparison of proven first-order formulas with the above-derived
second-order terms. The first two points are shown in the first analyzed
device, an atomic rotor in contact with two fixed leads (see Fig.
<ref>). The
second example consists of an oscillating quantum point contact (see
Fig. <ref>).
There we want to illustrate the kind of phenomena that can arise from
the second-order terms of the observables. In this sense, we should
recall that traditional adiabatic quantum pumping <cit.>,
requires the movement of at least two out-of-phase parameters. Instead,
the studied device shows monoparametric pumping thanks to the second-order
terms of charge currents, without requiring anharmonicities <cit.>
or quantum interference effects mediated by magnetic fields <cit.>.



 §.§ Driven atomic rotor


To model the atomic rotor we assume a minimal tight-binding model
consisting of a one-site system attached to two identical conduction
channels (the leads) through which an electrical current can flow.
We will name each of them as L and R. The isolated site is linked
to the rotor and represents a quantum dot with a site energy ϵ_d.
Two semi-infinite tight-binding chains set up the leads, where t_L
and t_R are the tunneling between them and the rotor <cit.>,
as sketched in Fig. <ref>.



To keep matters simple, we assume only the hopping terms change with
time, having a linear dependence on the distance between the rotor
and the respective leads. Furthermore, we force the atomic rotor to
follow a circular trajectory with a fixed radius R_0 and a given
angular frequency ω. Based on this outline, we arrive at (see
Appendix <ref>)

    t_L    =t_max(1+a-a√(1+β(1-sin(ω t)))),
    
    t_R    =t_max(1+a-a√(1+β(1-cos(ω t)))),

where t_max is the maximum (in absolute value) tunneling
amplitude, a sets the decay of the hopping terms with the distant
between the dot and the leads, and

    β    =2R_0L/(L-R_0)^2.


The description of the local system requires at least a three-site
Hamiltonian, which includes the site energy of the central dot and
the first sites of the leads. In this way, we circumvent the limitation
mentioned earlier of time-independent self-energies Σ^r,a,<(ε).
The resulting adiabatic Green's functions G^r,a,<, whose explicit
forms can be found in <cit.>, do not commute, in general,
with each other or Λ, as would be the case for a single site
local system; see the example of Ref. kershaw2017. Thus,
despite its simplicity, the analyzed example provides a challenging
test for the developed formulas.

Once we set up the model, we apply the charge current formulas. For
this purpose, we rewrite Eq. (<ref>)
(see Appendix <ref>)
for the current over the L-lead as

    I_L    =_I_L^(0)+_I_L^(1)ω+(_I_L^(1,1)+_I_L^(2))ω^2.

The term _I_L^(0) is the adiabatic contribution
given by Eq. (<ref>),
_I_L^(1) is the first-order correction
dictated by Eq. (<ref>),
and _I_L^(1,1) and _I_L^(2)
stand for the second-order corrections of the charge current, given
by Eqs. (<ref>)
and (<ref>)
respectively. 

The procedure used for the charge current can be extended to the remaining
studied observables. For example, the heat current formula takes the
form

    J_L    =_J_L^(0)+_J_L^(1)ω+(_J_L^(1,1)+_J_L^(2))ω^2.


Despite their complexity, our goal is to ensure the validity of the
formulas. To achieve this point, we will use the order-by-order energy
conservation formulas for quantum transport <cit.>.
For our model, it reads

    Q_J_L^(n)-δμ_L/eQ_I_L^(n)+Q_J_R^(n)-δμ_R/eQ_I_R^(n)+W^(n-1)    =0.

 The above formula holds for any integer n, where n represents
the order of the adiabatic expansion of the quantity of interest,
W^(n-1) is the n-1 order of the work done by the
rotor in a single cycle (see Appendix <ref>),
Q_J_α^(n) is the heat pumped per revolution
of the rotor to the α-lead, and Q_I_α^(n)
is the pumped charge in a cycle. These quantities are evaluated as:

    W^(n-1)    =[0]τ∫F^(n-1)·Xdt
    
    Q_J_α^(n)    =[0]τ∫J_α^(n)dt,
    
    Q_I_α^(n)    =[0]τ∫I_α^(n)dt.

Here, J_α^(n) and I_α^(n)
are defined in Section <ref>, τ is the period
of the movement, and F^(n-1) is the
n-1 order of the adiabatic expansion of the electronic force F.
The latter is a vector whose components are defined by

    F_ν(t)     =     -iħTr[Λ_ν𝒢^<(t,t)].


It is important to highlight that the fulfillment of Eq. (<ref>)
with second-order corrections of the heat and charge currents involves
the comparison of developed formulas with the well-known expressions
for the zero and first-order electronic forces, F^(0)
and F^(1), which can be found in Refs. <cit.>.

A further magnitude of significance is the total energy pumped to
the α-lead, which reads:

    Q_Ė_α^(n)    =Q_J_α^(n)-δμ_L/eQ_I_α^(n).

Rewriting Eq. (<ref>)
using the above formula, we obtain a version of the first law of thermodynamics
which relates the total energy of electrons flowing between leads
and the mechanical work. However, in this version electron energy
and mechanical work are related through the different orders of their
adiabatic expansion.

Finally, charge conservation is another property that must be satisfied
for every order and can easily be tested in the present model

    Q_I_L^(n)+Q_I_R^(n)    =0.


Fig. (<ref>)-(a)
illustrates how the different orders of the charge current add up
over one period of the rotor. Figs. <ref>-(b)
and <ref>-(c)
verify that the studied model satisfies charge conservation at every
order for different values of the dot's energy ϵ_d. Figures
<ref>-(d)
and <ref>-(e)
show, for different values of the dot's energy ϵ_d, order-by-order
energy conservation for the first and second adiabatic corrections
of the heat and charge currents.

We can also use the present example to highlight the importance of
the developed formulas for the thermodynamic analysis of systems.
For example, in Fig. <ref>
we can see that there is charge pumping in both orders [see (b)
and (c)] but energy pumping only occurs at first order. Note
that for n=1 [Fig. (d)] and at ϵ_d close to
0, energy is on average going out of the left lead (blue line below
zero) and entering into the right lead (red line above zero). On the
contrary, at second-order of the currents [n=2, see (e)],
energy coming from the external driving just dissipates through both
leads (red and blue curves are always positive). Here, W^(1)
[geen line of Fig. <ref>
(e)] is the mechanical energy being dissipated by the electronic
friction, see section <ref>
and Ref. bustosmarun2019. The negative value of W^(0)
in Fig. (d) (the work per cycle done by CIFs) is due to the sign
choice of Ẋ(t), which in the present case implies
that the system is being forced to act as a pump, not as a motor,
see Ref. bustosmarun2019. As a final remark, we note that
for the present example second-order pumping points in a different
direction than the first-order one. This implies that even small deviations
from the adiabaticity should diminish the efficiency of adiabatic
quantum pumps.






 §.§ Oscillating QPC or STM's tip


The system studied in this section consists of two conductors (not
necessarily of the same materials) close enough to let charge particles
tunnel from one to the other. Furthermore, the distance between the
conductors oscillates at a given frequency ω. This may represent
different physical situations, which are in principle within experimental
possibilities <cit.>. Some of them are depicted in Fig.
<ref>.



We modeled the different leads L and R using two distinct semi-infinite
homogeneous tight-binding chains, as in the previous example, each
characterized by a site energy and a hopping. The site energy of the
lead L (R) is ϵ_0L (ϵ_0R) and its hopping
is t_0L (t_0R). The hopping between the leads t_LR
depends on their distance, exhibiting an exponential behavior (see
Appendix <ref>)

    t_LR    =t_maxe^a(1-x/x_min),

where t_max is the maximum tunneling amplitude, a
is the decay factor, x is the distance between the leads, and x_min
is the minimum value of x. Assuming a forced oscillatory movement
and taking δ x+x_min as the maximum length, the
distance between leads yields

    x     =x_min+δ x/2(1-cos(ω t)).


For the calculation of charge and heat currents, we rewrite the corresponding
formulas just as we did in Eqs. (<ref>)
and (<ref>).
Moreover, for evaluating the fulfillment of the order-by-order energy
conservation in the present model, we will also resort to the same
definitions given in Eqs. (<ref>),
(<ref>),
and (<ref>). 



In Fig. <ref>
we evaluate the current for a system in equilibrium, i.e. (μ_L,T_L)=(μ_R,T_R),
with nonsymmetrical leads, which in this case means different energy
sites ϵ_0L≠ϵ_0R. There, the zero-order heat
J^(0) and charge currents I^(0)(adiabatic contributions)
are zero, as can be seen for I^(0) in Fig. <ref>-(a)
(blue solid line). Unlike the former example, the present one has
only one independent mechanical degree of freedom. This implies a
null contribution per cycle of the first-order heat and charge currents,
J^(1) and I^(1), see Ref. brouwer1998. Note in
Fig. <ref>-(a)
that the cyan curve, I^(0)+I^(1), is antisymmetric with respect
to I=0. However, the second-order term may allow particles and
heat to pump through the leads, even for a single-parameter system.
Note that, the dotted black line in Fig. <ref>-(a),
I^(0)+I^(1)+I^(2)+I^(1,1), is not antisymmetric for I=0
but has a shift towards negative values. In this example, it is then
clear that second-order terms comprise the leading order for charge
and heat pumping, exhibiting an unexplored form of quantum pumping.

In Figs. <ref>-(b)
and <ref>-(c),
we showed that the second-order term preserves, respectively, the
total charge and energy of the complete system (local system plus
the electron's leads). From this figure, it is also interesting to
note that, despite the charge current being pumped due to the second-order
terms [see (b)], there is no heat pumping. In Figs. <ref>
(c) we can see that the heat per cycle Q_J^(2) is positive
for left and right leads, i.e., the external work (W) is dissipated
as heat into the two leads. This mean that particles pumping from
the left lead,[According to Eq. (<ref>),
the number of particles pumped per cycle is -Q_I_L^(2)/e.] does not compensate for the external work (W) being dissipated
as heat to this lead. This highlight, as in the previous example,
the usefulness of the developed formulas for thermodynamic analysis.

For a symmetrical system (identical leads), we found (not shown in
the manuscript) that no net particles move from one lead to another
in the former configuration, even for the second-order correction.
This behavior is expected since, due to inversion symmetry, there
is no reason for the current to go in a preferential direction. In
this sense, this also validates the second-order formulas.



§ CONCLUSION


Throughout this work, we have presented expressions, based on NEGFs
within a Schwinger-Keldysh approach, for the second-order adiabatic
expansion of three different observables of quantum transport (energy,
heat, and charge currents). To our knowledge, NEGF second-order expansions
of heat and energy currents have not been treated before. Moreover,
the explicit formulas for the observables are general and ready to
be used in a wide range of systems. We only assumed that self-energies
are time-independent, which is naturally accomplished for most quantum
transport problems.

We have illustrated how the developed formulas produce physically
consistent results (including order-by-order energy and particle conservation
laws) using two simple models. We must emphasize that order-by-order
energy conservation strictly involves the comparison of the here-developed
expressions (second-order heat and charge currents) with a well-established
first-order expression (first-order current-induced force). Moreover,
we have shown in both examples how the developed formulas provide
useful thermodynamical information about the devices. Additionally,
in the second model, we have analyzed a phenomenon where second-order
monoparametric pumping takes place. This provides a glimpse of the
type of phenomenon that our expressions allow one to describe.

As for possible extensions of the present work, it would be interesting
a deeper exploration of second-order phenomena and their potential
use, as well as the development of similar formulas for other observables
like current-induced forces (also known as electron wind forces <cit.>)
or spin-transfer torque <cit.>.



§ ACKNOWLEDGMENTS


This work was financially supported by Consejo Nacional de Investigaciones
Cientficas y Tcnicas (CONICET, PIP-2022-59241); Secretara de Ciencia
y Tecnologa de la Universidad Nacional de Crdoba (SECyT-UNC, Proyecto
Formar 2020); and Agencia Nacional de Promocin Cientfica y Tecnolgica
(ANPCyT, PICT-2018-03587).





§ WIGNER TRANSFORM OF DYSON’S EQUATION


Eq. (<ref>)
is the equation of motion of the retarded Green's function, usually
called Dyson's equation of motion. We can apply a Wigner transform
to this expression. The first term on the right-hand side of Eq. (<ref>)
is the Dirac delta, whose Wigner transform is the identity. Before
applying the Wigner transform to the second term, we need to rewrite
it as

    A(t,t')     =𝒢^R(t,t')H_s(t')
        =∫𝒢^R(t,t_2)h(t_2,t')dt_2,

where

    h(t_2,t')     =H_S(t_2)δ(t_2-t').

Now, we can apply the Wigner transform to A(t,t'), by
using the Moyal product <cit.>, yielding

    A    =𝒢^R(T,ε)e^iħ/2(∂_ε∂_T-∂_T∂_ε)Ĥ_S(T).

where we used h(T,ε)=H_s(T).
The exponential operator is interpreted as usual in terms of the Maclaurin
series of the exponential function, i.e., e^ô=∑_nô^n/n!
The arrows in the exponential indicate the direction in which the
derivatives should be applied. For example, the symbol ∂_ε
means: derivative with respect to ε of the function to
the right, Ĥ_S in the above case, while ∂_ε
is the same but the derivative should be applied to the function to
the left, 𝒢^R(T,ε)
in the above case.

The third term of Eq. (<ref>),

    B(t,t')     =∫𝒢^R(t,t_1)Σ^R(t_1,t')dt_1,

can also be Wigner transformed by using the Moyal product, giving

    B    =𝒢^R(T,ε)e^iħ/2(∂_ε∂_T-∂_T∂_ε)Σ^R(ε).


To deal with the right-hand side of Eq. (<ref>)
we need to first apply the chain rule,


    ∂_t'𝒢^R    =1/2∂_T𝒢^R-∂_τ𝒢^R,

where T and τ are defined in Eqs. (<ref>)
and (<ref>) respectively.
Then, after Wigner transform it we obtained

    (∂_t'𝒢^r)    =1/2∂_T𝒢^R(T,ε)+iε/ħ𝒢^R(T,ε),

where we used the fact that Green’s functions have
compact support in the energy domain. The final result of all the
above is Eq. (<ref>)
of the main text.

In Eqs. (<ref>) and
(<ref>) we used a compact
notation, based on a exponential operator. However, the same equations
can also be written as two infinite sums

    A    =[N=0]∞∑1/N!(iħ/2)^NC_N^A(T,ε),

and

    B    =[N=0]∞∑1/N!(iħ/2)^NC_N^B(T,ε),

where

    C_N^A    =[j=0]N∑(-1)^N-jNj∂_T^N-j∂_ε^j𝒢^R∂_T^j∂_ε^N-jH_s,
    
    C_N^B    =[j=0]N∑(-1)^N-jNj∂_T^N-j∂_ε^j𝒢^R∂_T^j∂_ε^N-jΣ^R.

On the one hand, the slow time T only affects the local system,
according to our assumption. Therefore, the dynamics of the leads
only depends on the energy ε. On the other hand, the
Hamiltonian of the local system H_s only depends on T and
not on ε. Then

    ∂_ε^N-jH_s    =0 ∀ N-1>j≥0.
    ∂_T^jΣ^R    =0 ∀ N≥ j>1.

This can be used in Eq. (<ref>)
of the main text to obtain a simpler expression. The results is

    1     =-iħ/2∂_T𝒢^R+ε𝒢^R-[N=0]∞∑1/N!(iħ/2)^NC_N.

where we have introduced Eqs. (<ref>)
and (<ref>)
into Eqs. (<ref>)
and (<ref>)
respectively, and we have defined

    C_N    =(-1)^N∂_T^N𝒢^R∂_ε^NΣ^R+∂_ε^N𝒢^R∂_T^NH_s.

The zero and first order terms of last expression can be compared
with Eq. 20 of Ref. bode2012Jan.



§ ADIABATIC EXPANSION FOR RETARDED GREEN'S FUNCTION


The adiabatic retarded Green function, given in Eq. (<ref>),
satisfies the condition

    [G^R(T,ε)]^-1    =ε I-H_s(T)-Σ^R(T,ε).

Then, the derivative of G^R with respect to ε gives

    ∂_εG^R    =-G^r(1-∂_εΣ^R)G^r.

where we used ∂ A=-A∂[A^-1]A.

Next, we have to change the sequence of the terms in the Eq. (<ref>)
to introduce the retarded Green's function and its derivative, given
in Eqs. (<ref>) and (<ref>).
Then, using the slow-time derivative, Eqs. (<ref>),
and multiplying the adiabatic retarded Green’s function
on both sides of the obtained equation, the resulting expression takes
the form

    𝒢^R    =G^R+iħ/2[ν=1]M∑∂_εΣ^R_νG^RẊ_ν
        -iħ/2∂_T𝒢^R[G^R]^-1∂_εG^R+[N=2]∞∑1/N!(iħ/2)^NC_N.

Importantly, no approximations have been made in Eq. (<ref>).
Therefore, as long as the adiabatic expansion is valid (the series
is convergent), the equation provides the exact evolution of 𝒢^R
along a path driven by the mechanical degrees of freedom X_μ.
However, note that 𝒢^R (the exact retarded
Green's function) depends in turns on the time derivative of 𝒢^R,
which is not known. In principle, the only retarded Green's function
known in advance is the adiabatic one, G^r. To solve this, one
applies the iterative method twice. This consists of replacing 𝒢^R
on the right hand side by a first approximation (G^r), then using
the resulting Green function to improve the approximation of 𝒢^R
on the right hand side, and so on. During this iterative process one
identify terms of different orders given by the order of the derivatives
with respect to X_μ and cut the infinite series at a given
order (second order in our case). Note that the order also coincide
with the exponent of ħ or the exponent of the driving frequency
ω for the case X_μ∝cos(ω t).
By using this method we arrive, after some algebra, at Eq. (<ref>)
of the main text.



§ ADIABATIC EXPANSION FOR LESSER GREEN'S FUNCTION


Starting from Eq. (<ref>),
which gives the lesser Green’s function Wigner transforms
in close notation, we can write its gradient expansion as a double
infinite sum

    𝒢^<    =[N,K=0]∞∑1/N!K!(iħ/2)^N+KW_N,K^<(T,ε),

where we have introduced the nested definition

    W_N,K^<    =[j=0]N∑[l=0]K∑(-1)^NNjKl∂_T^N-j∂_ε^j𝒢^RC_j,l^N,K

and

    C_j,l^N,K    =(-1)^j+l∂_T^j∂_ε^N-j{∂_T^K-l∂_ε^lΣ^<∂_T^l∂_ε^K-l𝒢^A} .

The aforementioned expression is obtained by straightforwardly making
use of the Moyal product twice.

The lesser self-energy is not time-dependent, since we have assumed
that the retarded self-energy is not controlled by the slow time.
This implies

    ∂_T^K-lΣ^<    =0 ∀ N-l≠0,
    ∂_T^j-mΣ^<    =0 ∀ j-m≠0.

Putting Eqs. (<ref>)
and (<ref>)
in Eq. (<ref>),
result in the simplified expression

    𝒢^<    =[Q=0]∞∑1/Q!(iħ/2)^QW_Q^<(T,ε),

where

    W_Q^<    =[q=0]Q∑[s=0]q∑[j=0]Q-q∑(-1)^qQq+jq+jjqsC_q,j^Q,

and

    C_q,j^Q    =∂_T^q∂_ε^j𝒢^R∂_ε^Q-(j+s)Σ^<∂_ε^s∂_T^Q-q𝒢^A.

The zero and first-order terms of Eq. (<ref>)
coincide with Eq. 25 of Ref. bode2012Jan and Eq. 3 of Ref.
bode2011, once the iterative method explained in the previous
section is applied. The result up to the second order is shown in
Eq. (<ref>) of
the main text.



§ ADIABATIC EXPANSION FOR CHARGE CURRENTS


We start by rewriting Eq. (<ref>)
as


    I_α    =2eRe{tr[F_α(t)]} ,

where the auxiliary operator F_α is defined by

    F_α    =∫[𝒢^R(t,t')_α^<(t',t).
       .+𝒢^<(t,t')_α^A(t',t)]dt'

Now, our goal is to transform the F_α operator to the energy
domain and then apply the inverse Wigner transform. This allows us
to recover the charge current in the time domain written as an integral
with a kernel in the energy domain. To carry out his technique, we
start with the relationship between F_α and its Wigner transform
F_α, whose formula takes the form

    F_α    =1/2πħ∫F_α(T,ε)dε.


The next step is to put the gradient expansion in effect over Eq.
(<ref>),
leading to

    F_α    =[N=0]∞∑[j=0]N∑1/N!(-iħ/2)^NNjC_α,N,j^F(T,ε),

where

    C_α,N,j^F    =(-1)^j∂_T^N-j∂_ε^j𝒢^R∂_T^j∂_ε^N-j_α^<
        +(-1)^j∂_T^N-j∂_ε^j𝒢^<∂_T^j∂_ε^N-j_α^A.


As we have mentioned, we assumed that the leads are unaffected by
the slow time variation. For such reason. Therefore, the following
rules hold

    ∂_T^j_α^A    =0 ∀ N≥ j≥1,
    ∂_T^j_α^<    =0 ∀ N≥ j≥1.

Putting Eqs. (<ref>)
and (<ref>)
in Eq. (<ref>),
the Wigner transform of the F_α operator reads

    F_α    =[N=0]∞∑1/N!(-iħ/2)^NC_α,N^F(T,ε),

where

    C_α,N^F    =∂_T^N𝒢^R∂_ε^N_α^<+∂_T^N𝒢^<∂_ε^N_α^A.

To get the desired result, we must first plug the last equation into
Eq. (<ref>)
and then into Eq. (<ref>),
leading to the formula outlined in Eq. (<ref>).



§ ENERGY CURRENT MAIN FORMULA


Taking Eq. (<ref>)
as a starting point, we intend to get an expression of the energy
current in terms of self-energies, which lets us identify the leads.
The quoted formula can be read as the product of the hooping between
the local system and the leads, with those elements of the lesser
Green’s function that connect the local system and
leads, given by Eq. (<ref>).
However, we want to express the energy current in terms of the Green's
functions and self-energies of the local system only, Eq. (<ref>).
For this purpose, we will use the Keldysh technique and the Langreth
theorem of analytic continuation to achieve this goal <cit.>.
The mean result is

    J_α    =-2Re{∫𝒥_α(t,t')dt'} ,

where

    𝒥_α    =l,s∑(𝒢_l,s^R(t,t')f_α,s,l^<(t',t)+𝒢_l,s^<(t,t')f_α,s,l^A(t',t)).

Here, 𝒢_l,s^A are the elements of the retarded Green’s
function given in Eq. (<ref>),
and 𝒢_l,s^< are the lesser Green’s
function specified in Eq. (<ref>).
In addition, we have introduced the following auxiliary functions

    f_α,s,l^<    =k∑ϵ_α kt_s,α kg_α k^<(t',t)t_α k,l,
    
    f_α,s,l^A    =k∑ϵ_α kt_s,α kg_α k^A(t',t)t_α k,l.


As you can notice, these propagators satisfy the differential equations:

    (iħ∂/∂ t'-ϵ_α k)g_α k^<(t',t)     =0,
    (iħ∂/∂ t'-ϵ_α k)g_α k^A(t',t)     =δ(t'-t).

Moreover, and following our assumptions, we also have:

    ∂_t'(t_s,α kt_α k,l)     =0.

Then, putting Eqs. (<ref>)
and (<ref>)
into Eqs. (<ref>)
and (<ref>),
respectively, and using the previous condition, we arrive at

    f_α,s,l^<    =∂_t'_α,s,l^<(t',t),
    
    f_α,s,l^A    =∂_t'_α,s,l^A(t',t)+k∑t_s,α kδ(t'-t)t_α k,l.

Note that we have used the definitions of the self-energies given
in Eqs. (<ref>)
and (<ref>). The
last equations allow us to write the following

    𝒥_α    =𝒥_α^(0)(t,t')+𝒥_α^(1)(t,t'),

where

    𝒥_α^(A)    =tr{𝒢^R(t,t')_α^0}δ(t'-t),
    𝒥_α^(B)    =tr{𝒢^R(t,t')∂_t'_α^<(t',t)+𝒢^<(t,t')∂_t'_α^A(t',t)} .

For compactness, we have defined above the operator

    _α,s,l^0    =k∑t_s,α kt_α k,l.


With this, the energy current takes the form

    J_α    =J_α^(A)+J_α^(B),

where

    J_α^(A)    =-2Re{∫𝒥_α^(A)(t,t')dt'} ,
    
    J_α^(B)    =-2Re{∫𝒥_α^(B)(t,t')dt'} .

However, with the aid of the properties (𝒢^<)^=-𝒢^<
and (_α^0)^†=_α^0,
the following can be proved

    J_α^(A)    =0.

Therefore, the total energy current is set only by J_α^(A),
resulting in Eq. (<ref>).



§ ADIABATIC EXPANSION FOR ENERGY CURRENTS


In this section, we apply a method analogous to the one used for the
charge current. We start by defining the auxiliary operator

    M_α    =∫[𝒢^R(t,t')∂_t'_α^<(t',t).
       .+𝒢^<(t,t')∂_t'_α^A(t',t)]dt'

This enables us to write the energy current as

    J_α    =-2Re{ iħtr[M_α(t)]} .


Once again, the procedure is to apply the Wigner transform to the
M_α operator and afterward the inverse transform. Then,
the gradient expansion of Eq. (<ref>)
takes the form

    M_α    =[N=0]∞∑[j=0]N∑1/N!(-iħ/2)^NNjC_α,N,j^M(T,ε),

where

    C_α,N,j^M    =(-1)^j∂_T^N-j∂_ε^j𝒢^R∂_T^j∂_ε^N-j∂_t'_α^<
        +(-1)^j∂_T^N-j∂_ε^j𝒢^<∂_T^j∂_ε^N-j∂_t'_α^A.

where ∂_t'_α is the Wigner
transform of ∂_t'_α(t',t).

The charge current formula, shown in Eq. (<ref>),
can be compared with the above transform. The main difference between
both lies in the self-energies, where the one given in Eq. (<ref>)
implies the Wigner transform of the self-energy time derivative. To
find the expressions for these terms, we just need to use the definitions
of the Wigner coordinates, given in Eqs. (<ref>)
and (<ref>), followed by
the chain rule. We then apply the Wigner transform to the resulting
expressions and assume that the self-energies, which are solely energy-dependent
(_α(ε)≡_α),
vanish at high and low energies. The result is

    ∂_t'_α^<    =1/2∂_T_α^<-iε/ħ_α^<,
    ∂_t'_α^A    =1/2∂_T_α^A-iε/ħ_α^A.


Then, plugging the Eqs. (<ref>)
and (<ref>)
into Eq. (<ref>),
together with the conditions of Eqs. (<ref>)
and (<ref>),
the Wigner transform of M_α can be written as

    M_α    =-1/2[N=0]∞∑1/N!(-iħ/2)^N-1C_α,N^M(T,ε),

where

    C_α,N^M    =∂_T^N𝒢^R∂_ε^N(ε_α^<)+∂_T^N𝒢^<∂_ε^N(ε_α^A).

We can apply the chain rule to ∂_ε^N(ε_α^<)
giving

    ∂_ε^N(ε_α^<)     =ε∂_ε^N_α^<+N∂_ε^N-1_α^<,
    ∂_ε^N(ε_α^A)     =ε∂_ε^N_α^A+N∂_ε^N-1_α^A.

These results allow us to rewrite the M_α of
Eq. (<ref>)
in a closed form

    M_α    =-iε/ħF_α-1/2[N=0]∞∑1/N!(-iħ/2)^NC_α,N(T,ε),

where the term F_α is that defined for the charge
current (see Eq. (<ref>)),
and

    C_α,N    =∂_T^N+1𝒢^R∂_ε^N_α^<+∂_T^N+1𝒢^<∂_ε^N_α^A.


Finally, we need to apply the inverse Wigner transform to Eq. (<ref>)
and then put this into the energy current definition [Eq. (<ref>)].
The result is Eq. (<ref>).



§ DRIVEN ATOMIC ROTOR MODEL


In subsection <ref>, we provide a brief
account of the fundamental constituents of the atomic rotor model.
We will use two semi-infinite tight-binding chains to represent the
leads attached to the quantum dot and include the first site of each
chain as part of the local system. Then, the Hamiltonian of the local
system reads:

    𝐇_S    =[  ϵ_0 -t_L    0; -t_L  ϵ_d -t_R;    0 -t_R  ϵ_0 ].


The adiabatic retarded Green's function, given by Eq. (<ref>),
takes the form:

    𝐆^r    =(ε𝐈-𝐇_S-Σ_L^R-Σ_R^R)^-1,

where the self-energies coming from the decimation of the left and
right leads (Σ_L^R and Σ_R^R
respectively) are given by

    Σ_L^R    =[ Σ_0^R(ε)        0        0;        0        0        0;        0        0        0 ],    Σ_R^R    =[        0        0        0;        0        0        0;        0        0 Σ_0^R(ε) ].

Here, Σ_0^R is the self-energy (in the energy domain)
of a semi-infite tight-binding chain. It is often expressed as

    Σ_0^R(ε)     =Δ_0(ε)-iΓ_0(ε),

where

    Δ_0    =t_0ε-ϵ_0/2t_0-√((ε-ϵ_0/2t_0)^2-1)   ε-ϵ_0/2t_0≥1,
    ε-ϵ_0/2t_0    1≥ε-ϵ_0/2t_0≥-1,
    ε-ϵ_0/2t_0+√((ε-ϵ_0/2t_0)^2-1)    -1≥ε-ϵ_0/2t_0,
    
    Γ_0    =t_0
    0    ε-ϵ_0/2V_0≥1,
    √(1-(ε-ϵ_0/2t_0)^2)    1≥ε-ϵ_0/2t_0≥-1,
    
    0     -1≥ε-ϵ_0/2t_0.

Above, ϵ_0 and t_0 are, respectively, the site energy
and the hopping between neighboring sites of the tight-binding chain.

In our model, the mechanical DOF affects only t_L and t_R,
which can be taken as generalized coordinates related to the position
of the rotor θ by some function. Then, the matrix operator
Λ_θ can be written as


    _θ    =    ∂ t_L/∂θ_t_L+∂ t_R/∂θ_t_R.

Here, the matrix operators _t_L
and _t_R, necessary to calculate
the charge and heat currents (see sections <ref>
and <ref>), as well as the electronic force
(see Section <ref>), are:


    _t_L    =-[ 0 1 0; 1 0 0; 0 0 0 ]  and   _t_R    =-[ 0 0 0; 0 0 1; 0 1 0 ].


So far, we have provided a general tight-binding model for a quantum
dot connected with two leads. The next stage requires the specification
of the geometric parameters of the atomic rotor. Typically the dependence
of hopping parameters with the distance are modeled by exponential
functions. In our case, that would mean taking

    t_L    =t_me^a(1-Δ r_L/r_0L),
    
    t_R    =t_me^a(1-Δ r_R/r_0R).

where Δ r_L and Δ r_R are the spatial distances
between the quantum dot and the closest chain sites to the L and
R leads, respectively, while r_0L and r_0R are the smallest
values of Δ r_L and Δ r_R along the trajectory.
Using polar coordinates, we have

    Δ r_L    =√(L^2+R_0^2-2R_0Lsin(θ)),
    Δ r_R    =√(L^2+R_0^2-2R_0Lcos(θ)).

where θ, L, and R_0 are defined in Fig. <ref>,
and r_0L and r_0R satisfy the condition

    r_0L    =r_0R=L-R_0.


Plugging the Eqs. (<ref>),
(<ref>)
and (<ref>)
into Eqs. (<ref>)
and (<ref>)
make the model available to evaluate charge and heat currents. However,
since our purpose is to outline typical behaviors of the currents
and not to focus on the particularities of the used models, we decided
to use a linearized version of Eqs. (<ref>)
and (<ref>).
In this way, we arrive at Eqs. (<ref>)
and (<ref>).

Taking into account the model, we are now able to put the charging
current calculation into action. For a fixed radius and angular position
given by θ=ω t, we rewrite Eqs. (<ref>),
(<ref>), (<ref>)
and (<ref>)
as

    I_L^(0)    =_I_L^(0),     I_L^(1,1)    =_I_L^(1,1)ω^2,
    
    I_L^(1)    =_I_L^(1)ω,     I_L^(2)    =_I_L^(2)ω^2,

where

    _I_L^(1)    =-e(𝒩_L,t_L^(1)∂ t_L/∂θ+𝒩_L,t_R^(1)∂ t_R/∂θ),
    _I_L^(2)    =-e(𝒩_L,t_L^(2)∂^2t_L/∂θ^2+𝒩_L,t_R^(2)∂^2t_R^2/∂θ^2),
    _I_L^(1,1)    =-e(𝒩_L,t_Lt_L^(1,1)(∂ t_L/∂θ)^2+𝒩_L,t_Rt_R^(1,1)(∂ t_R/∂θ)^2)
        -e(𝒩_L,t_Lt_R^(1,1)+𝒩_L,t_Rt_L^(1,1))∂ t_L/∂θ∂ t_R/∂θ.


Similar expressions can be obtained for the heat currents, see Sections
<ref> and <ref>.



§ DRIVEN QUANTUM POINT CONTACT MODEL


In subsection <ref>, we discussed
the physical behavior of this kind of device and provided a brief
introduction to its modeling via two tight-binding semi-infinite chains
with tunneling between them. From this starting point, our next step
is to give the Hamiltonian of the local system, which reads

    𝐇_S    =[  ϵ_0L -t_LR; -t_LR  ϵ_0R ].


The adiabatic retarded Green's function assumes the same form as Eq.
(<ref>),
The self-energies in this case are

    Σ_L^R    =[ Σ_0L^R(ε)         0;         0         0 ],    Σ_R^R    =[         0         0;         0 Σ_0R^R(ε) ].

where Σ_0L^r (Σ_0R^r) is given by Eq. (<ref>)
but replacing the site energy ϵ_0 with ϵ_0L
(ϵ_0R) and the hopping t_0 by t_0L (t_0R).

Since in this system only one parameter moves with time, we have a
single _ν matrix given by

    =-[ 0 1; 1 0 ].

Finally, the dependence of the hopping with the separation between
leads is given in Eqs. (<ref>)
and (<ref>).



§ ELECTRONIC FORCES


To verify the order-by-order energy conservation, the evaluation of
the work done by the electronic forces at a given order n is necessary.
The work W at a given order n is

    W^(n)    =[0]τ∫(ν∑F_ν^(n)Ẋ_ν)dt,

where Ẋ_ν is the slow time derivative of the classical
mechanical degree of freedom ν, F_ν^(n) is
the n-th order of the adiabatic expansion of the electronic force
F_ν. The expressions for F_ν^(n) up to
first-order are well known (see for example Refs. bode2012Jan,deghi2021).
They read

    F_μ^(0)    =-1/2π i∫tr{_μG^<} dε,
    
    F_μ^(1)    =-ν∑γ_μνẊ_ν,

where

    γ_μν    =ħ/2π∫tr{ G^<_μ∂_εG^R_ν-G^<_ν∂_εG^A_μ} dε.






