Heat and charge transport in interacting nanoconductors driven by time-modulated temperatures

We investigate the quantum transport of the heat and the charge through a quantum dot coupled to fermionic contacts under the influence of time modulation of temperatures


Introduction
Nanoelectronics [1] and the emerging field of Thermotronics [2][3][4][5] are at the forefront of manipulating electron charge and energy fluxes through electrostatic and thermal gradients, respectively.The interplay between these technologies, known as Thermoelectricity [6] explores how a thermal gradient influences charge current and vice versa, expanding the functionality of quantum conductors.Investigating time-dependent transport in nanostructures unlocks unique insights not achievable with static fields [7][8][9].The ability to control these nanoscale systems with electrical or thermal drivings paves the way for new possibilities in quantum technologies.Electrically modulated quantum conductors lead to innovative devices like electron pumps [10][11][12][13][14][15], dynamical Coulomb Blockade quantum systems [16,17], and AC-driven nano-electromechanical systems for sensors [18].This setup can function as a quantized emitter, operating as a quantum capacitor in the adiabatic regime and behaving like an RC circuit in the GHz range, with a unique charge relaxation resistance quantized as R 0 = h/2e 2 [7,19,20].They have demonstrated to serve as single electron sources for quantum computing applications and metrology [12,[21][22][23][24][25][26].
Similarly, thermally modulated nanoconductors, within the framework of quantum thermodynamics [27], are rapidly advancing.Examples include quantum thermal machines [26,[28][29][30][31][32], thermal diodes [33,34], thermal transistors [35][36][37][38], thermal memristors [39] and thermal capacitors [2][3][4].Despite significant progress, time-dependent transport in nanostructures, especially in thermally driven systems, remains a challenging and vibrant field.Thus far, researchers have primarily tackled this issue in two scenarios: the incoherent transport regime [40] and the adiabatic driving [41][42][43][44], characterized by a slow time-modulation leading to a net transport of charge or heat.Adiabatic pumping for the design of thermal machines [45,46] is rooted in geometrical concepts like the Berry connection and it has very broad applications in electronics, thermotronics [5], and quantum information processing.Recent studies on temperature-driven dynamics in two-level systems engineer an effective temperature through oscillator frequency manipulation [47].However, most of the existing temperature-driven studies heavily rely on scattering theory, lacking the consideration of Coulomb interaction beyond mean-field treatment and being unable to capture dynamical excitations induced by nonadiabatic temperature driving in the presence of nontrivial interactions.
Figure 1: Lateral quantum dot system coupled to a left reservoir (l) and a right reservoir (r) that are described by H C,l and H C,r .Each reservoir is under the influence of a modulated temperature in time with T l (t) and T r (t).Left and right tunneling barriers are described by H T,l and H T,r , respectively, as indicated.The central part corresponds a spinful quantum dot.Plunger gates V l and V r control the barrier transparency, an additional gate V g is applied to the quantum dot region to tune the dot level position denoted by ε σ .
New and exciting functionalities are expected to arise in time-dependent thermally-driven interacting nanoconductors away from the adiabatic regime.This issue is the central objective of our work.The primary challenge in introducing a thermal bias lies in achieving a microscopic formulation.In this context, reconciling the macroscopic nature of temperature gradients and thermal forces, which arise from statistical averaging, appears to be incompatible.Then, representing these effects in a microscopic quantum mechanical Hamiltonian is not straightforward.However, in 1964, Luttinger [48,49] presented a clever and ingenious solution.He introduced a scalar potential, Ψ, referred as a gravitational field, which couples to the energy density of the system, J E .This scalar field can be viewed as a mechanical conjugate to the energy density of the system defined as follows Luttinger justified his trick by the request for this scalar field to satisfy the Einstein relation, i.e., the potential adjusts itself to balance the thermal force, resulting in an identity ∇Ψ = ∇T /T in the thermal equilibrium.In this respect, the gravitational field or thermomechanical potential serves as a local proxy for local temperature variations.This theoretical treatment finds application in determining linear responses to static gradients by introducing the scalar field [50,51] or equivalently, a vector potential [52].It proves effective in various scenarios, including classical systems [53], thermoelectrical transport in quantum systems using densityfunctional theory framework [50] and transient current calculations [51,54] in the stationary regime.In the latter, the adiabatic limit for time-dependent temperature is analyzed using the nonequilibrium Keldysh Green function (NEGF) formalism [45,46,49,55], employing Luttinger's trick.
Our work applies Luttinger's idea to a quantum conductor tunnel-coupled to two electronic reservoirs subjected to time-dependent temperature modulation [see Fig. 1].We focus on a quantum dot (QD), the most generic quantum conductor, featuring multiple local levels with diverse interactions such as Coulomb interaction and spin-orbit coupling.Our goal is to derive generic expressions of the electric and heat currents flowing through the QD junctions in the framework of the NEGF formalism.For that purpose we adopt the gravitational field description to address the temperature modulation in time.Our formalism to be proposed is quite widely applicable in that it can consider diverse and complex quantum conductors and, in particular, it allows one to deal with the effect of strong Coulomb interaction in a systematic way, in contrast to previous methods.Also, it can be easily extended to involve geometric configurations such as quantum mechanical interferometers.
Our formalism, in the course of its development, addresses a very important technical aspect regarding the correct definition of the dynamical contact energy and its coupling to the gravitational field.We adopt a tight-binding model to describe the system, which divides the total Hamiltonian into pieces, namely, the contact (reservoir) regions, the quantum dot, and the tunneling barriers connecting the contacts to the QD.The contact energy responsible for the heat current in the corresponding contact should then include not only the energy stored in the contact ℓ = l, r (described by the Hamiltonian H Cℓ , see Eq. ( 3)) but also the half of the energy stored in the tunneling barrier coupled to that contact (described by the Hamiltonian H Tℓ , see Eq. ( 5)) [56,57].Together with additional physical requirements (to be explained in Sec. 3) it is natural to couple the Luttinger field to the contact energy Q ℓ (which therefore includes H Tℓ ).This is the most important ingredient for a correct application of the Luttinger's trick to the calculation of the dynamical heat current.
Accommodating the Luttinger's trick and employing the NEGF technique, we formulate expressions for charge and heat currents with a particular emphasis on the linear response regime.Notably, our formalism yields comprehensive expressions for the currents in relation to the QD NEGFs [see Eqs.(34) and (37)].In special conditions, by a help of the charge conservation and a sum rule on the energy change rates, the currents can be obtained solely in terms of the retarded and advanced components of the QD NEGFs [see Eq. ( 41)], which streamlines the computation of the currents.It should be noted that our expressions for currents are the analogue to those for electrical current induced by a time-dependent electrical modulation originally presented in the pioneering work by Jauho, Wingreen and Meir [8], but in our case the time-dependent modulation is done thermally instead.
In demonstrating our Luttinger formalism, we initially apply it to the noninteracting case in which the fulfillment of the Onsager reciprocity relation validates our approach and moreover our results align with those obtained from the scattering theory.Extending our formalism to the interacting case within the Hartree approximation reveals that the Coulomb interaction can alter the responses for charging and energy relaxations with distinct temperature dependences.The success of our formalism in these applications suggests its potential for studying dynamical heat transport in the presence of the Coulomb blockade effect or many-body correlations.
Our paper is organized as follows: In Sec. 2 we introduce our model Hamiltonian and implement the Luttinger's trick onto it.We also define the charge and heat currents and find out the relevant sum rules.In Sec. 3 we express the charge and heat currents in terms of QD NEGFs by solving the corresponding Dyson's equations for the NEGF.Later, we restrict ourselves to the linear response regime and further simplify the expressions for the charge and heat currents with a help of sum rules, which are the main results of our work.Sections 4 and 5 demonstrate the applications of our Luttinger formalism, considering the noninteracting and interacting cases, respectively.In Sec. 6, we summarize our work and discuss the possible applications and extensions of our formalism and the experimental proposals.

Thermal fields and currents
For our study, we consider a quantum dot coupled to two fermionic reservoirs.To represent this configuration, we establish the Hamiltonian in the context of a tight-binding model.To construct this Hamiltonian, we first consider the distinct components that comprise our nanoconductor, and then provide an explanation for how to add the term corresponding to the gravitational field.The tight-binding Hamiltonian comprises three contributions.First, the Hamiltonian for the electrodes, denoted as ℓ that takes the values of left (l) and right (r), is given by where c † ℓkσ is the creation operator for an electron in the lead ℓ with wavevector k, spin σ(=↑, ↓), and energy ε ℓk measured with respect to the Fermi level.Second, the central conductor considered as an interacting localized multi-level quantum dot is described as Here d m annihilates an electron on the localized level m where the index m denotes both the orbital level and the spin.H D can describe any kind of interactions inside the quantum dot such as the Coulomb interactions and the spin-orbit couplings.Thirdly, the tunneling Hamiltonian that connects each of the two electrodes with the central conductor reads where the tunneling amplitude between the lead ℓ and the dot level m is denoted by t ℓkσ,m .Hereafter, for simplicity, we assume that the tunneling amplitude is momentumindependent t ℓkσ,m ≡ t ℓσ,m .Accordingly, we define the escaping tunneling rates Γ ℓ,mm ′ ≡ σ πρ 0 t ℓσ,m t * ℓσ,m ′ /ħ h for level m and m ′ due to the coupling to the contact ℓ, where ρ 0 is the density of states in the contacts at the Fermi energy.Later it is convenient to express the rates in the form of a N d × N d matrix Γ ℓ whose matrix elements are given by Γ ℓ,mm ′ .

Luttinger's trick and Hamiltonian
The idea about the Luttinger's trick consists of introducing new fields, dubbed as gravitational fields Ψ ℓ (t) which are coupled to the contact energies.We determine the precise form of this coupling based on the following two arguments.First, as mentioned before, the contact energies in the framework of a tight-binding formulation should be redefined to account for the reactance or energy stored in the barrier [56], leading to Q ℓ as specified in Eq. ( 2).Therefore, in our setup, the term responsible for the dynamical thermal driving is introduced as where each of the Ψ ℓ (t) field is coupled to the excess energy for the contact ℓ with respect to its equilibrium value (at Ψ ℓ = 0); 〈• • • 〉 0 denotes the expectation value at equilibrium.The coupling to the excess energy is reasonable because for sufficiently weak driving or at low temperature only the excitations around the Fermi level will contribute to the electric and thermal transports.Furthermore, mathematically, it introduces only the additional time-varying number, − ℓ Ψ ℓ (t) 〈Q ℓ 〉 0 which does not affect the dynamics of states.The constant λ in Q ℓ , the coupling coefficient of the gravitational field to the tunneling barrier energy, is set to be 1/2, but for a time being we keep this symbol as it is in order to track down how the precise value of this coefficient influences the results.Second, we require that the coupling to the gravitational field should not affect the effective coupling between the dot and the leads [49,55] since the dot-lead coupling should be immune to the temperature in the leads.Interestingly, we find that the coupling in Eq. ( 6) automatically satisfies this requirement, at least, in the linear response regime, as we will see later [see The complete Hamiltonian for our setup is then described by where we have introduced the time-dependent contact and tunneling Hamiltonians defined as for α = C, T and with λ C = 1 and λ T = λ.We assume that the thermal drivings on both the contacts are periodic with the same frequency Ω = 2π/τ: Ψ ℓ (t +τ) = Ψ ℓ (t).More specifically, we take a sinusoidal time dependence: Here the factors Ψ ℓ are real and can be zero.Since the dynamical variation of the temperature is taken into account via the gravitaional field, both the contacts are assumed to have the same chemical potential and the same base temperature T (or the inverse temperature β) so that the thermal populations of both the uncoupled contacts are specified by the same Fermi distribution function f (ε).

Charge currents and charge conservation
We focus in our study on the charge and heat currents transversing the contacts, driven by the dynamical change of the temperature.The charge currents, or the change rates of the charges in the contacts (ℓ = l, r) are described in terms of the charge current operators where N ℓ ≡ kσ c † ℓkσ c ℓkσ is the charge number operator for the contact ℓ.Under this consideration, the charge currents, the expectation values of the charge current operators for the contacts and the quantum dot can be obtained by evaluating where is the charge number operator for the quantum dot.It should be noted that the charge conservation condition, [H, ℓ N ℓ + N D ] = 0, guarantees that their sum vanishes at all time t: Later we will use this equality to simplify our final results.

Heat currents and sum rule
The heat currents for the electrodes are derived from the time derivative of the energies stored at the contacts.Then, according to our choice of the contact energy Q ℓ which incorporates the contribution from the neighboring tunneling barrier, the heat current for contact ℓ is given by [56, 57] It should be noted 1 that multiple choices for the definition of the heat currents for contacts are a priori possible in the Luttinger formalism, while we find that Eq. ( 13) is the most suitable one.The power supplied by the time-dependent thermal source or the power dissipated is defined as P(t) = 〈∂ H/∂ t〉 = 〈∂ H G /∂ t〉 and explicitly expressed as which contains the source contributions from the contacts and the tunneling barriers.Under the thermal driving, the energy conservation does not hold: P(t) ̸ = 0.However, one can still find an equality similar to Eq. (12).To this purpose, we define the energy change rates as Then, from Eq. ( 7), the obvious commutation relation [H(t), H(t)] = 0 leads to a useful equality which is to be exploited later.

Charge and heat currents in terms of NEGF's
In the seminal research [8], the charge current flowing through a quantum dot under a dynamical electric drive was formulated in a closed form in terms of the interacting QD Green functions.In a similar way, we formulate here the charge and heat currents in terms of solely the QD NEGFs when the system is driven by oscillating temperatures.For such purpose we adopt the nonequilibrium Keldysh formalism and employ the equation-of-motion technique 1 A potential argument could state that the appropriate selection of the contact energy operator is Q ℓ = H Cℓ,Ψ (t) + λH Tℓ,Ψ (t), since it incorporates the sources, that is, the contribution to the energy from the coupling terms of the gravitional fields.Evaluating this alongside Eq. ( 2), we observe their difference being second order in Ψ ℓ .As we are operating in a linear response regime, both definitions offer equivalent outcomes for the heat current to this order.The determination of the correct contact energy operator definition requires a study into the nonlinear regime, which currently surpasses our research scope.Besides, the power definition is of second order in Ψ ℓ , which means that the source contribution to the heat fluxes is of second order and it is consistent with the indistinguishability between the two candidate definitions for the contact energy.
together with the Langreth rules [58].Throughout our paper, we use the retarded/advanced and lesser QD NEGFs defined as respectively, and the contact and QD-contact NEGFs are defined in a similar way: For examples, Based on the definition of the Green's functions, the charge currents (11) and the expectation values for the contact energy operators (2) are readily expressed in terms of the NEGF´s: and Note that above we have introduced a time-varying effective tunneling amplitudes, for convenience.Knowing the time dependence of two expectation values, 〈H Cℓ 〉 and 〈H Tℓ 〉, one can evaluate not only the heat currents (13) and the power (14) but also the energy change rates for contact ℓ and the tunneling parts via Now, by employing the equation-of-motion technique [8], the lesser Green's functions, G < m,ℓkσ (t, t) and G < ℓkσ,ℓkσ (t, t), can be cast in terms of solely the QD NEGFs [see Appendix A for details].After some algebraic manipulations on G < m,ℓkσ (t, t), the contact charge current (19a) and the expectation value of energy stored in the tunneling barrier (20b) have compact forms in terms of self energies and the QD NEGFs: Here G a and Σ a (a = R, A, <) are the N d ×N d matrix representations of the QD Green's functions and the self energies, respectively.The traces are done over the QD orbital degrees of freedoms.The matrix elements of the self energies are defined as where g a ℓkσ (t, t ′ ) with a = R, A, < correspond to the uncoupled contact Green's functions governed by the time-dependent Hamiltonian, ε ℓk (t)c † ℓkσ c ℓkσ , with the time-varying energy, It should be noted that the gravitational field Ψ ℓ (t) enters into the self energy through the timevarying tunneling amplitude t ℓkσ,m (t) as well as the dynamically driven contact energy ε ℓk (t).
Similarly, by expressing G < ℓkσ,ℓkσ (t, t) in terms of the QD NEGFs, the expectation value of the contact Hamiltonian is written as where E Cℓ0 ≡ −i kσ ε ℓk g < ℓkσ (t, t) are the equilibrium energies contained in the unperturbed contacts.Note that this constant is irrelevant in our study of heat flow because it does not contribute to the heat current once the time derivative is taken.The matrix elements of the self-energy-like term Ξ a b ℓ in Eq. ( 26) are defined as with a, b = R, A, <.
In the presence of the time-dependent terms in the Hamiltonian, the Green's functions G(t, t ′ ) depend not on the time difference t − t ′ but on t and t ′ separately.However, since the driving is periodic with period τ, the Green's functions and the self energies are also periodic with G(t +τ, t ′ +τ) = G(t, t ′ ) so that it is more convenient to apply the Fourier transformation and to move from the time domain to the frequency domain.We adopt the mixed time-energy representation for the Fourier transformation: Then, the integrals in Eqs. ( 23) and ( 26) can be expressed in terms of the Fourier components: with Importantly, the charge and heat currents [see Eqs. ( 23) and ( 26) together with Eq. ( 13)] are now expressed in the frequency domain employing Eqs.(29a) and (29b) solely in terms of the QD NEGF.This constitutes an exact form for the charge and heat currents of an interacting conductor coupled to fermionic reservoirs.
One may want to interpret the above ac-driven current in terms of the photo-assisted tunneling as done in the Tien-Gordon approach [59].In this approach, the effect of the driving is transformed to the appearance of the quasi-energy states (ω + nΩ) due to the absorption/emission processes of photons and the resulting current is then given by the sum of all possible processes, each of which is weighted by proper Bessel functions with argument proportional to Ψ ℓ .This simple interpretation does not work in our case though.The main reason is that, while in the Tien-Gordon approach the time-dependent field is coupled to the charge numbers (c † ℓkσ c ℓkσ ) only, in the temperature-driven case it is coupled to the excitation energies (ε ℓk ) as well.Then, the photon-assisted processes with different n are no longer independent of each other but are instead intermingled with each other and have very complicated energy dependencies [see Eq. ( 84)].Hence, it is not practical to interpret the temperature-driven current in terms of the photon-assisted processes.
Hereafter we adopt an approximation scheme using a linear expansion in the Luttinger field Ψ ℓ (t).The reasons are three-fold: (1) The Luttinger scheme was originally proposed to work only in the linear response regime, (2) the second requirement of our Luttinger setup, making the effective dot-contact coupling independent of Ψ ℓ , is satisfied only in the linear response regime, and (3) in a practical manner the use of the above exact form is limited because it is quite difficult to obtain manageable analytical expressions for Σ a ℓ (n, ω) and Ξ ab ℓ (n, n ′ , ω) working for arbitrary magnitude of Ψ ℓ .

Linear response regime
After deriving the general forms for the charge and heat currents in the system, our intention is to supply a manageable formulation of such currents in terms of the equilibrium or dynamical QD Green's functions.This objective is reachable within the linear response regime, in which we treat the amplitudes Ψ ℓ as the smallest parameters and keep up to their linear order in the charge and heat fluxes.
Considering that our periodic driving ( 9) is of the form, Ψ ℓ (t) = (Ψ ℓ /2)(e iΩt +e −iΩt ), in the linear expansion with respect to Ψ ℓ , only the Fourier components with n, n ′ = 0, ±1 of Σ a ℓ (n, ω) and Ξ a b ℓ (n, n ′ , ω) should remain finite.While their explicit expressions and derivations can be found in Appendix B, one should note that from Eq. (83) the linear expansion of Σ so that Σ R/A ℓ (t, t ′ ) becomes independent of the gravitational field only when λ = 1/2.That is, only at this choice of λ, the requirement that the dot-contact hybridization should be immune to the temperature change is met.
In the same spirit of the linear expansion, only the n = 0, ±1 Fourier components of the QD Green's functions should survive: G a (0, ω) is in the zeroth order of Ψ ℓ , while G a (±1, ω) are linear in Ψ ℓ : up to the linear order in Ψ ℓ , In particular, the n = 0 components of the self energies and the QD Green's functions should exactly correspond to their equilibrium values at Ψ ℓ = 0.
In the linear response regime, the charge and heat currents are described by their first Fourier components only, with ℓ/D (Ω)] * since no current flows at Ψ ℓ = 0.In the following sections, we are going to express the Fourier components I c/h ℓ (Ω) of the charge and heat currents in terms of the equilibrium QD Green's functions, G a (0, ω) and the nonequilibrium linear components, G a (±1, ω).

Charge/heat currents and power
In order to obtain the charge current in the linear response regime, we express the contact charge current (23a) in terms of the Fourier components G a (n, ω) and Σ a ℓ (n, ω) by using Eq.(29a) and then keep only the linear-order terms in Ψ ℓ [refer to Appendix B].The derivation is quite straightforward and we obtain Tr and where we have defined [refer to Appendix C for details].As expected, the charge current depends not only on the equilibrium QD Green's functions but also on the dynamical ones, G R/A/< (1, ω), even though the linear response (Ψ ℓ → 0) is taken.It is obviously because our perturbations are dynamical and the dynamical excitations of the system, even though being small, cannot be described solely in terms of the equilibrium Green's functions.Now we turn to the heat transport.In order to find the expressions for the heat currents, the power, and the energy change rates in the linear response regime, we write down the energies (23b) and (26) in terms of the Fourier components of G a (n, ω), Σ a ℓ (n, ω), and Ξ ab ℓ (n, n ′ , ω) by using Eqs.(29a) and (29b) and then keep only the linear-order terms in Ψ ℓ [see Appendix B for details].Following the explicit derivation in Appendix C, the contact heat current can be explicitly obtained as where iΩ 4 E Tℓ0 is an additional unphysical term.Equations ( 34) and (37) are our main results.
The expression given in (37) for the contact heat current needs a discussion.The additional term I h ℓT (Ω) which is proportional to E Tℓ0 is an artefact of our Luttinger's trick.In our setup, we dynamically drive the contact by the field Ψ ℓ (t) and the tunneling barrier by the field λΨ ℓ (t) so that an effective energy capacitor which is dynamically driven by the field difference (1−λ)Ψ ℓ (t) is formed.However, this effect is not contained in the original system and is solely due to the Luttinger's setup itself.This artificial setup then gives rise to an additional heat transfer (1 − λ)Ψ ℓ (t)E Tℓ0 (up to the linear order) between the contact ℓ and the tunneling barrier, resulting in I h ℓT (Ω) in the contact heat current.In fact, in the next section for the noninteracting system, we compare the results from our Luttinger's trick and those obtained from the scattering theory based on the dissipation-fluctuation theorem and find that two results are identical only when this artefact is not taken into account.Therefore, we will drop out the term I h ℓT (Ω) from now on.Finally, we find the expression for the power of dissipation (14) in the linear response regime [see Eq. ( 101)].In particular, we focus on the time average of the power which is simplified to As expected, it reflects that the time-averaged power is directly related to the real part of the Fourier component of the contact heat currents.It is because the dissipation happens at the contacts and the time average picks up only the dissipative effect.

Application of conservation law and sum rule
Note that our expressions for the charge and heat currents, Eqs. ( 34) and ( 37) necessitate the knowledge of the nonequilibrium components G < (1, ω) as wells as G R/A (1, ω).Technically, it is much harder to theoretically obtain the lesser Green's functions than the retarded/advanced Green's functions because the lesser ones reflect not only the energy excitations of the system but also their nonequilibrium distribution.Therefore, our expressions for the currents would be more useful if the knowledge of the lesser Green's functions is avoided, especially when the quantum dot is interacting.
We have found that it is possible as long as the hybridization matrix Γ ℓ is proportional to the identity matrix: This condition cannot accommodate the general situations in real experiments, but most of qualitative features can be still captured within this condition.So this simplifying condition is acceptable at least for theoretical studies.First, the expression for the contact charge current, Eq. ( 34) can be simplified further if the charge conservation is taken into account.If we apply the charge conservation (12) up to the linear order, then we have which enables one to write down the integral dω Tr[G < (1, ω)] in terms of the other components of the QD Green's functions [see Eq. ( 93)].Then, we can get a nice expression for the interacting charge current at the contact ℓ in the linear response regime which writes solely in terms of the retarded/advanced QD Green's functions: with The expression of the contact heat current, Eq. ( 37) can be also further simplified in a similar way.Recall that we have the sum rule (16) for the energy change rates, which gives rise to in the linear response regime.This sum rule enables us to replace the integral dω ω Tr[G < (1, ω)] with ones of other QD Green's functions.However, this sum rule cannot be constructed without knowing explicitly the form of the QD Hamiltonian H D since the change rate W D (Ω) depends on H D .Therefore, our strategy is as follows: Once H D is known, we find out the expression for W D (Ω) in terms of the QD NEGFs in the linear response regime and then write down the integral dω ω Tr[G < (1, ω)] in terms of other QD NEGFs by using the sum rule (43) and the explicit expression for the partial sum, ℓ (W In the following sections, we apply our formalism to two specific examples: the noninteracting case and the interacting case with the Hartree approximation.Especially, in the noninteracting case, we demonstrate the justification of the choice λ = 1/2 in more details.

Noninteracting case
As a first application of our formalism derived in the previous sections, we consider the noninteracting case in which the QD Hamiltonian (4) is given by We assume that the dot-lead coupling is spin-independent so that Γ, now being 2 × 2 matrix, is diagonal and given by Γ = Γ 1.The orbital index m is now replaced by the spin index σ, and the trace over m is now the summation over σ.It is then quite straightforward to derive and solve the Dyson's equations for the QD NEGFs, G R/A σ (t, t ′ ) and to find out their Fourier components in the linear response regime [refer to Appendix D for detailed derivations].The n = 0 (equilibrium) components of the retarded/advanced QD Green's functions are found to be and their n = 1 components exactly vanish at ℓσ (±1, ω) = 0. On the other hand, G < σ (±1, ω) does not vanish at λ = 1/2, so the dynamical components of the QD Green's functions are still relevant in time-dependent charge and heat transports.However, as explained in the previous section, we do not seek out the explicit expression for G < σ (±1, ω), but instead resort to the charge conservation (40) and the sum rule (43) to derive the integrals of G < σ (±1, ω).Starting from the general formulas of the charge and heat currents, Eqs. ( 41) and (37), one can derive the explicit expressions for the charge and heat currents [see Eqs.(107) and where the diagonal components are decomposed into The self admittances are then found to be and the cross admittances are given by where we have defined These linear-response admittances can be expressed in terms of an equivalent RC circuit, as shown in Fig. 2.Under this equivalence, the cross admittances L l r and K l r represent the thermoelectric (Y = L) and thermal conductance (Y = K) between the two contacts, coming from the parallel configuration of a resistor R Y,l r and a capacitor C Y,l r so that On the other hand, the self admittances L ℓ and K ℓ represent the electrical/thermal conductances for charging/heating and relaxing in the contact ℓ, coming from the serial configurations of a resistor R Y,ℓ and a capacitor C Y,ℓ so that Generally, the above resistances and capacitances are functions of the frequency Ω.However, in the low-frequency (ħ hΩ ≪ k B T ) limit, they can be approximated to constants.We also assume the low-temperature (k B T ≪ Γ ℓ ) condition to apply the Sommerfeld approximation.
The self resistances and capacitances are then approximated to and where h is the thermal relaxation resistance for a single mode mesoscopic capacitor.On the other hand, the low-frequency cross resistances and capacitances are found to and The detailed analysis of the resistances and capacitances will be present in Sec.4.2.

Why the choice of λ = 1/2 ?
In the original works [56,57] which explored the heat transport in time domain, it is found that the meaningful definition of the contact energy should be determined to be Q ℓ (t) with λ = 1/2 [see Eq. ( 2)] because only this choice is consistent with the first and second laws of thermodynamics.While this argument alone justifies the choice of λ = 1/2, in this section we intend to find more evidence which supports the choice of λ = 1/2 by comparing our results for I c ℓ (Ω) and I h ℓ (Ω) with the previous ones obtained by different methods for the similar systems.By using the equation-of-motion method, Rosselló, López, and Lim [60] have investigated the dynamical heat current through the similar setup as ours but with a single contact which is driven by an ac electric voltage.In calculating the heat current, they also took into account the energy barrier contribution with λ = 1/2, as proposed in Ref. [56].As expected, our self thermoelectric admittance L ℓ (Ω) is exactly equal to their self electrothermal admittance M ℓ (Ω)/T [see Eq. (36) in Ref. [60]] which measures the heat current through the contact with respect to the ac electric driving in the contact.It justifies that the dynamical gravitational field in our setup should be coupled to Q ℓ with λ = 1/2 in order to get the correct dynamical charge current.Note that this agreement, M ℓ (Ω) = T L ℓ (Ω) (where T is the background temperature, the common temperature in the contacts) also reflects the fact that the reciprocal relation, or the so-called Onsager's relation [61] should hold in the thermal transport.However, in the point of view of the fluctuation-dissipation theorem, both admittances L ℓ (Ω) and M ℓ (Ω) in the linear response regime are related to the same fluctuation In our setup, the perturbative (gravitational) field is coupled to Q ℓ and 〈dN ℓ /d t〉 is measured, while in Rosselló's work the external (electric) field is coupled to N ℓ and 〈dQ ℓ /d t〉 is measured.So, apparently, the agreement might be mathematically trivial because both used the same The second previous work to be compared with is the work done by Lim, López, and Sánchez [62] which has applied the scattering theory to the single-contact QD setup to obtain the dynamical charge and heat current in the linear response regime when the contact is driven either by an ac voltage or by an ac temperature.Note that in the scattering theory approach the barrier plays the role of the scatterer only, so no energy is stored in it and the contact energy is defined with respect to H ℓ only.They calculated the low-frequency responses of the currents with respect to the ac voltage and temperature, via the fluctuation-dissipation relation which they assumes to hold.We have found that our low-frequency expansion of the self admittances, Eqs. ( 52) and ( 53) are in good agreement with the scattering-theory predictions [see Eqs. ( 7), (8), and (9) in Ref. [62]].This agreement strongly justifies our use of λ = 1/2, especially because they are derived from two different approaches: We have directly calculated the dissipative part by adopting the Luttinger's trick, while in Ref. [62] the admittances were obtained by calculating the fluctuations based on the scattering theory.Also, this agreement implies that the fluctuation-dissipation theorem holds for thermal transport, at least in the non-interacting and single-contact case.

More analysis on charge/heat currents
Since the physical discussion on the low-frequency self admittances, R L/K,ℓ and C L/K,ℓ has already been done in Ref. [62], we focus on the cross admittances here.First, the low-frequency and low-temperature expansions of the cross resistances R L/K,l r are physically reasonable and are in agreement with the scattering-theory prediction.For example, the low-temperature cross thermal conductance 1/R K,l r [see the black dotted line in Fig. 3 (a)] was predicted to be proportional to the electric conductance (e 2 /h)(2Γ l Γ r /Γ ) σ ρ σ (0) [63], which is well reflected in Eq. ( 55), while the low-temperature cross thermoelectric conductance 1/R L,l r [see the black dotted line in Fig. 3 (c)] is proportional to σ ρ ′ σ (0) as expected.It should be noted that the fluctuation-dissipation theorem applied to the heat transport through two-contact systems is no longer valid because scattering events that connect two different terminals induce a nonvanishing term for the equilibrium heat-heat correlation function at the low temperature limit, which is incompatible with the expected behavior of K l r (Ω) [63,64].Therefore, one cannot exploit the fluctuation-dissipation theorem to study the dynamic heat transport through the quantum-dot systems described in terms of the tight-binding model.It signifies that our Luttinger formalism is the promising candidate for the systematic study of dynamical temperature driving.
Figure 3 displays the dot-level dependencies of the thermal/thermoelectric conductances and capacitances for a wide range of temperatures.The thermoelectric conductance 1/R L,l r and capacitances C L,l r share similar dependence on the dot level, whose qualitative feature does not change much as the temperature increases [see Figs. 3 (c) and (d)].While the thermal conductance 1/R K,l r and capacitances C K,l r also share similar dependence on the dot level, they change from single-peak shape to double-peak one as the temperature increases [see

(d)
Figure 3: (a) The cross thermal conductance 1/R K,l r and (b) the cross thermal capacitance C K,l r , and (c) the cross thermoelectric conductance 1/R L,l r and (d) the cross thermoelectric capacitance C L,l r in the low-frequency limit as functions of the spin-degenerate dot level ε d in the noninteracting case with the symmetric coupling, Γ l = Γ r .Each curve corresponds to the different temperature k B T /ħ hΓ whose value is annotated.The resistances and capacitors for the cross admittances are defined via Eq.( 50) and the corresponding cross thermal and thermoelectric admittances are evaluated via Eq.( 48) by using given values of temperatures and the value of frequency chosen numerically as small as possible.The black dotted lines correspond to the low-temperature limits as given by Eqs. ( 54) and (55).
Figs. 3 (a) and (b)].Recall that the heat current depends not only on the carrier occupation but also the carrier energy.At high temperatures, high-energy carrier can make more contribution to the heat current, which is the reason why the heat current can be larger at the off-resonant condition, as demonstrated in Figs. 3 (a) and (b).
Very interesting property peculiar to the noninteracting condition can be found in the RC times defined as for Y = L, K.In the noninteracting case, the self and cross RC times are always equal to each other, that is, for both Y = L, K.It is because the self and cross admittances share the same frequency dependence: From Eqs. ( 47) and ( 48), one can find that L ℓ (Ω) and L l r (Ω) are proportional to σ P 1σ (Ω), while K ℓ (Ω) and K l r (Ω) are proportional to σ P 2σ (Ω).As we will see in the next section, this is not the case for interacting QD systems.That is, the comparison between the self and cross response times can be used to measure the effect of the interaction.

(b)
Figure 4: (a) The RC times τ K,l r = τ K,ℓ from thermal admittance and (b) the RC time τ L,l r = τ L,ℓ from thermoelectric admittance in the low-frequency limit as functions of the spin-degenerate dot level ε d in the noninteracting case with the same parameters as in Fig. 3.The black dotted lines correspond to the low-temperature limits as given by Eq. ( 59).
shows the dot-level dependence of the RC times.The black dotted lines correspond to the low-temperature limit which is given by While τ L,l r shows the behavior similar to the QD density of states which broadens as the temperature increases, τ K,l r features the double-peak structure at high temperatures.
Finally, the time-averaged power (38) for the noninteracting case is obtained as We then found that in the low-frequency limit only the second term remains finite so that only the cross thermal resistance R K,l r is responsible for the energy dissipation.Interestingly, this second term for the dissipation is identical to its electric counterpart, V 2 /2R where V is the electric voltage drop and R is the electric resistance between the two contacts.We again would like to stress that Eq. ( 60) is a natural outcome obtained by following the procedure based on our Luttinger formalism, without resorting to some heuristic arguments.Therefore our formalism is proven to provide a systematic way to investigate the dynamical heat transport in the linear response regime.

Interacting case: Hartree approximation
Our formalism is not limited to the noninteracting case.The charge and heat currents, Eqs. ( 41) and ( 37) can be calculated as long as the equilibrium and n = 1 components of the retarded/advanced QD Green's functions are provided.In this section we take into account the Coulomb interaction in the quantum dot which is now described by Unfortunately, in the presence of finite Coulomb interaction, it is impossible to obtain any analytical form of the QD NEGFs without a proper approximation.Here we adopt the simplest (d) approximation: the Hartree approximation which reduces the two-particle correlations to oneparticle ones as Then, the Dyson's equations for the QD NEGFs are found to be basically similar to those for the noninteracting case but with the retarded/advanced self energies being now replaced by Note that lesser self energy Σ < σ (t, t ′ ) remains unchanged compared to the noninteracting case.This additional term in Σ R/A σ,HF (t, t ′ ) induces two changes compared to the noninteracting case: (1) The effective dot level is shifted from the unperturbed one, where is the equilibrium QD occupation which should be determined in a self-consistent way.Note that the equilibrium QD Green's functions now depend on ε σ,HF : (2) The n = 1 Fourier component of the QD Green's functions are now finite: where is the n = 1 Fourier component of the QD occupation.By using the charge conservation, one can obtain the explicit expression of n σ (1, Ω) which is found to be with Then, following the recipe proposed in our formalism [see Appendix E for further explanations], the charge and heat currents can be obtained [see Eqs.(121a) and (121b)] and the corresponding self/cross thermoelectric and thermal admittances are found to be and These admittances clearly display the corrections due to the Coulomb interaction, as shown in Fig. 5: The resonance is shifted and the curves are slightly deformed compared to the noninteracting case, but no qualitative changes are observed.For examples, the low-frequency and low-temperature cross resistances R L,l r and R K,l r are found to be identical to those in the noninteracting case except ε σ being replaced by ε σ,HF .Therefore, it is not convenient to find a solid evidence on the effect of the Coulomb interaction from the dot-level dependence of the resistances and capacitances.
Instead, we focus on the RC times.In the presence of the Coulomb interaction, the self and cross RC times are not equal to each other any longer:  71) with Eq. ( 72)] and make their frequency dependence different from each other.As matter of fact, the difference between them is proportional to U for weak Coulomb interaction.By performing the low-frequency and low-temperature expansions   similar to those done in the noninteracting case, one can find that, for spin-degenerate case with ρ ↑ (0 These expansions show that (1) both τ L,l r −τ L,ℓ and τ K,l r −τ K,ℓ are finite in the presence of the Coulomb interaction and (2) τ L,l r − τ L,ℓ saturates in the T = 0 limit, while τ K,l r − τ K,ℓ scales as T 2 .These asymtotic behaviors for low temperatures are clearly manifested for a rather wide range of the temperature, as demonstrated in Figs. 6 (c) and (d).We expect that these temperature dependencies of the differences between the RC times can be used to identify the effect of the Coulomb interaction as long as U is sufficiently small.It should be noted that these temperature dependences due to the Coulomb corrections originate from the dynamical excitations reflected in the n = 1 components of the QD Green's functions, G R/A< σ (1, ω), and not from the effective level shift in ε σ,HF .Our Luttinger formalism takes into account the dynamical excitations systematically and correctly, even in the presence of the Coulomb interaction.Therefore, the utility of our formalism becomes more evident in the interacting systems.
In order to investigate the effect of the strong Coulomb interaction on the thermoelectric and thermal admittances, one should go beyond the Hartree approximation, taking into account the higher-order terms in the equation-of-motion method.For example, as long as the Coulomb blockade is concerned, one can apply the Meir-Wingreen-Lee approximations [65,66] to our formalism, which will be our future work.

Conclusions
We have formulated a general formalism, based on the Luttinger's trick, to calculate the dynamic charge and heat currents through a multi-level quantum dot which is driven by timedependent temperatures.Our Luttinger formalism is built on the correct definition of the dynamical contact energy and the requirement for the effective dot-contact coupling.It provides the general expressions for the linear-response charge and heat currents, given by Eqs. ( 34) and (37), respectively, which can be calculated as long as the n = 0 (equilibrium) and n = 1 Fourier components of the retarded/advanced and lesser QD Green's functions are known.Furthermore, with the help of charge conservation and the vanishing sum of the energy change rates, the knowledge of the lesser QD Green's functions can be avoided as long as some conditions are met.
The physically important point of our formalism is that it goes beyond the static and adiabatic condition for temperature modulation and it can capture naturally and systematically the effect due to dynamical excitations driven by the nonadiabatic temperature driving.Our formalism is then expected to be very adequate and useful, considering the state-of-the-art technology which enables very fast modulation of temperatures.Furthermore, our formalism works even in the presence of the Coulomb interaction and can reveal the role of the interaction in nonequilibrium thermal responses.Its application to the interacting case studied in Sec. 5 clearly demonstrates the success of our formalism in this regard, even though it has been done only in the Hartree approximation.
While the applications of our formalism present in this paper [see Secs. 4 and 5] are limited to the single-level quantum-dot system, it can be applied to any multi-level quantum-dot junctions so as to study the effect of multiple levels and spin-orbit interactions on the temperaturedriven transport.Since our formalism allows the quantum dot to be any finite systems, it can also consider the case of multiple quantum dots.Furthermore, some interferometer-like geometry such as the junction with a quantum-embedded ring can be also incorporated into our formalism.As long as one ignores the Coulomb interactions in these systems, the equation-ofmotion method can be readily exploited to derive the required retarded/advanced QD Green's functions and subsequently the dynamical charge and heat currents.
Real challenge is to take into account the Coulomb interaction in a non-perturbative way.There is no analytical solution in this case without proper approximations.Usually, it is very hard to obtain the exact form of equilibrium QD Green's functions let alone their dynamical (n = 1) Fourier components.One promising method to deal with the Coulomb interaction in a non-perturbative way is to use the numerical renormalization group [67] in order to obtain the QD Green's functions in a numerical way.While originally the numerical renormalization group is restricted to the equilibrium case, recently its extension, so-called the time-dependent numerical renormalization group, is improved to deal with the nonequilibrium case with a periodic driving [68,69], so that the two-time QD Green's functions are successfully calculated.While it has some issues about the accuracy, we believe that this method can be safely used to study the linear response regime in which the driving is sufficiently weak.
Another interesting theoretical challenge is to extend the Luttinger formalism beyond the linear regime.The original proposal of the Luttinger's trick [48] was based on the linear response regime so that the gradient of the gravitational field is identified with the gradient of the temperature.In fact, there is no physically reliable justification of the use of the Luttinger's trick beyond the linear regime.However, recently some tried to extend the Luttinger's idea into the nonlinear-regime study of the nanostructures: the steady-state and transient behaviors of charge and heat currents after sudden quench of the gravitational field [50,54] and temperature-driven adiabatic pumping [49,55].The predictions made by these works are quite interesting and nontrivial, while their validity of the prediction is uncertain and to be confirmed in experiments.So, knowing that there is no systematic way to deal with the timedependent temperature, it may be very physically interesting to extend our formalism beyond the linear response regime so that a new physics is explored.
Finally, we would like to address briefly the experimental realization of our scheme.In order to test our predictions, a time-dependent modulation and control of temperature in the reservoirs should be experimentally implemented.We think, for example, that our theory may be tested in spin qubits [40] with detunnings of order of few meV that corresponds to ac frequencies for the temperature modulation of about hundreds of MHz.These are frequencies that are experimentally accessible nowadays.Alternatively, for higher frequencies (in the GHz and THz range), a design of time-dependent temperature signals by employing a collection of quantum harmonic oscillators that mediate the interactions between the quantum system and a thermal bath has been recently proposed in Ref. [47].It is worth noting that in order to realize our prediction, the driving frequency should be not too low.It is because our scheme is based on the quantum-mechanically coherent state during the driving so the period of the driving should be shorter than the decoherence time, which in turn requires frequencies in the GHz regime in quantum-dot systems [19].In short, the main hurdle to deal with is to maintain the quantum coherence long enough during the temperature modulation.
technique with respect to the Hamiltonian (7).For convenience, we introduce the time-varying tunneling amplitudes (21) and contact energy (25) so that the contact and tunneling Hamiltonians can be written as apart from the time-dependent numbers which do not affect the dynamics of the Green's functions.
Via the equation-of-motion method, it is straightforward to obtain the following Dyson's equations: where τ 3 is the third Pauli matrix in the Keldysh space and g ℓkσ (t, t ′ ) denotes the uncoupled contact Green's function matrix as introduced in the main text.Explicitly, the uncoupled contact Green's functions are given by where f (ε) is the Fermi function at the temperature T and the phase factor is evaluated as with In order to evaluate the charge and heat currents, we need to get the expressions for ), (20a), and (20b)] from the Dyson's equations (77a) and (77b): We insert the above expressions into Eqs.(19a), (20a), and (20b) and obtain Eqs. ( 23) and ( 26) in terms of the relevant self energies Σ a ℓ defined as Eq. ( 24) and the self-energy-like terms Ξ a b ℓ defined as Eq. ( 27).

B Explicit expressions and linear expansions of Σ
R/A/< ℓ (t, t ′ ) and To evaluate the charge and heat currents via Eqs.( 23) and ( 26), one needs to know the explicit expressions for the self energies Σ a ℓ (t, t ′ ) [see Eq. ( 24)] and the self-energy-like forms Ξ a b ℓ (t, t ′ , t ′′ ) [see Eq. ( 27)].For simplicity, we take the wide-band limit with a constant density of states ρ 0 for both the contacts, which allows us to replace the sum k F (ε ℓk ) by the integral ρ 0 ∞ −∞ dε F (ε).Using the explicit expressions for g R/A/< ℓkσ (t, t ′ ), one can express the self energies where ω = ε/ħ h and we have used f (ω) instead of f (ħ hω), for simplicity.The integration over ω can be done for Σ R/A ℓσ , giving rise to where we have used the fact that |Ψ ℓ | ≤ 1: The temperature oscillation amplitude cannot be larger than the base temperature T itself.The corresponding Fourier components, Σ R/A ℓ,n (ω) can be also obtained in an analytical form (we do not write down the detailed expression here) which involves the regularized generalized hypergeometric functions.However, unfortunately, no simple analytical expressions for Σ < ℓ (t, t ′ ) nor Σ < ℓ,n (ω) are available.Specifically, by using the identity e iβ sin Ωt = ∞ n=−∞ e iΩt J n (β) where J n (x) are the first kind Bessel functions and with a help of recurrence relations for J n (x), one can obtain with ω m ≡ ω − mΩ.For general value of Ψ ℓ , this expression is not adequate for analytical nor numerical analysis since it requires the summation over m from −∞ to ∞.
The situation becomes worse for Ξ ab ℓ (t, t ′ , t ′′ ) whose explicit expressions involve more complicated integration: Unfortunately, for general values of Ψ ℓ , the integrations over ω and the Fourier transformation do not yield any manageable analytical expressions.On the other hand, we have found that the linear expansion of the self energies with respect to Ψ ℓ yields reliable and analytical expressions.Therefore, in our study we focus on the linear response regime.Before presenting the linear expansion of the self energies, the reliability of the linear expansion should be examined.The linear expansion in Ψ ℓ approximates the exponential function in the integrals, Eqs. ( 82) and ( 85) as before the integration over ω is done.While Ψ ℓ is assumed to be small enough, in fact, ωΨ ℓ is not small for large values of ω which definitely happen during the integration over ω, which may disqualify the use of the linear expansion.However, we have confirmed that this expansion produces correct results.Our justifications are two-fold.First, in the linear response regime, only the contact excitations close to the Fermi level are relevant: Note that ε = ħ hω is the contact excitation energy.Hence, the correctness of the approximation at higher energies does not matter.Secondly, we have explicitly adopted a regularization function η(ω) to the gravitational field so that the thermal driving is really effective only to the low-energy states: Ψ ℓ → η(ω)Ψ ℓ .Specifically, we have chosen a Gaussian regularization η(ω) = exp[−(ω/ω 0 ) 2 ] with a constant ω 0 which determines the range of energies to be meaningfully coupled to the gravitational field.Then, the linear expansion is well justified because ωη(ω) Ψ(t) is small for all values of ω: Note that η(ω) decreases exponentially with ω.Then, at the final stage of the calculations, we take ω 0 → ∞ to restore the original coupling of the gravitational field.We have confirmed that the result obtained from the regularization and taking the limit of ω 0 → ∞ is identical to that obtained by using the linear expansion, Eq. ( 86) from the beginning.
Applying the linear expansion (86) and performing the Fourier transformation (28) to and where the definition of ∆ f (ω, ω ′ ) is given by Eq. ( 36).As one can see, in the linear response regime, the n = 0 Fourier components are nothing but the equilibrium values at Ψ ℓ = 0, and the n = ±1 components are linear in Ψ ℓ while all the higher components (|n| ≥ 2) vanish up to the linear order in Ψ ℓ .One can note that the n = ±1 components become a lot simpler at λ = 1/2: In particular, Σ and and with Note that owing to the constant I ∞ the equilibrium contributions, Ξ RA/R</<A ℓ (0, 0, ω) are divergingly large.However, we have found that these diverging contributions cancel out each other exactly in the zeroth-order term of the linear expansion of Eq. ( 26): where in the last line we have used the equality, which holds generally in equilibrium for any interacting quantum-dot junctions.

C Linear regime: Charge and heat currents and sum rules
The charge currents in the linear regime can be obtained by expressing (23a) in terms of the Fourier components G(n, ω) and Σ ℓ (n, ω) via (29a) and then by keeping terms up to the linear order in Ψ ℓ .It is then quite straightforward to obtain the explicit expression, giving rise to Eq. ( 34).In the derivation, we have used the expressions for Σ R/A/< ℓ (n, ω) [see Eqs.(87) and ( 88)] with λ = 1/2.
Our system conserves the total charge so that the charge conservation ( 40 By inserting Eq. ( 93) into the contact charge current (34), one gets Eq. ( 41).
For obtaining the heat currents we first take the average energies that are then expanded into with and and the tunneling barrier energy at equilibrium Then, from Eqs. ( 13) and ( 22), the contact heat currents and the energy change rates in the linear regime are expressed in terms of E C/Tℓ (Ω): where the Fourier components of the energy change rates in the linear regime are defined as Employing these expressions and for λ = 1/2 the Fourier component of the heat current in the linear regime is obtained as (37).
We have another similar sum rule for the energy change rates, Eq. ( 16), which is written in the linear regime as Eq. ( 43).While W D (t) = d 〈H D 〉 /d t = (i/ħ h)〈[H, H D ]〉 requires the specification of H D , the other terms in the sum rule can be written as This expression can be used to write the integral dω ω Tr[G < (1, ω)] in terms of other QD Green's functions, via the sum rule (43).Finally, we find the expression for the power of dissipation ( 14) in the linear response regime.Up to the lowest order in Ψ ℓ , the power reads which is of the second order in Ψ ℓ .

D Noninteracting Case: QD Green Functions and Charge/Heat Currents
By applying the equation-of-motion method to the noninteracting Hamiltonian (44) By using the linear expansions (87) of Σ R/A σ (n, ω) and by keeping up to the linear order in Ψ ℓ , the n = 0 (equilibrium) and n = 1 components of the retarded/advanced QD Green's functions are found to be (1, ω).Note that at λ = 1/2, Σ R/A ℓ,σ (1, ω) is zero up to the linear order in Ψ ℓ so that G R/A σ (1, ω) also vanishes.Then, by exploiting the properties of the equilibrium noninteracting QD Green's functions (45) and the vanishingness of Σ R/A ℓ,σ (1, ω), one can get the explicit expression for the charge current from the general formula of the charge current (41): with with Ψ l = Ψ r and vice versa.
In order to find the explicit expression for the heat current, one should know the integrals of G < σ (1, ω), that is, dω G < σ (1, ω) and dω ωG < σ (1, ω).In the noninteracting case, one can easily derive the linear expansions of the lesser QD Green's functions G < σ (t, ω) from the Dyson's equation (103).However, as proposed in the main text, we instead exploit the charge conservation (93) and the sum rule (43) for the energy change rates in order to express the integrals of G < σ (1, ω), that is, dω G < σ (1, ω) and dω ωG < σ (1, ω) in terms of the equilibrium retarded/advanced QD Green's functions.Below we derive their explicit expressions.From the charge conservation (93), by using the properties of the equilibrium QD Green's functions (45) and G R/A σ (1, ω) = 0, one can get Refer the definition of P 1σ (Ω) to Eq. (49).For the noninteracting QD Hamiltonian, the energy change rate W D (t) is identified as and its Fourier component in the linear regime is found to be Now, by inserting the expressions for W D (Ω) and the partial sum (100) into the sum rule (43), one can find that By inserting this integral into the general expression for the heat current, Eq. ( 37), we obtain the explicit expression for the heat current for the noninteracting case.

E Interacting case: The Hartree approximation
In the presence of the Coulomb interaction described by the interacting Hamiltonian (61), the Dyson's equation for the QD NEGFs has an additional term proportional to U: where n σ(0) and n σ(1, Ω) are the n = 0 (equilibrium) and n = 1 Fourier components of the QD occupation 〈n σ (t)〉, given by Eqs. ( 65) and ( 68), respectively.Hence, the equilibrium retarded/advanced QD Green's functions are modified accordingly [see Eq. ( 66)] and, according to Eq. (106b), the n = 1 components of the retarded/advanced QD Green's functions are now finite [see Eq. ( 67)].Since G R/A σ (±1, ω) are now finite, the integrals of dω G < σ (1, ω) and dω ωG < σ (1, ω) will have additional terms.First, from the charge conservation (93), by using the properties of the equilibrium QD Green's functions (66) and the explicit expressions (67) By inserting this integral into the general expression for the heat current, Eq. ( 37), we can get the explicit expression for the heat current.Finally, the charge and heat currents in the Hartree approximation are obtained as

Figure 2 :
Figure 2: Equivalent RC circuit in the low temperature and low frequency limits.The RC circuit describes both the thermoelectrical (with Y = L) and the thermal admittance (with Y = K) in terms of resistors R Y,ℓ , capacitors C Y,ℓ and the cross resistance and capacitances R Y,l r , and C Y,l r .Refer to the main text for their specifications for either the thermoelectrical or the thermal transports.

Figure 5 :
Figure5: (a) The cross thermal conductance 1/R K,l r and (b) the cross thermal capacitance C K,l r , and (c) the cross thermoelectric conductance 1/R L,l r and (d) the cross thermoelectric capacitance C L,l r in the low-frequency limit as functions of the spin-degenerate dot level ε d in the interacting case with U/ħ hΓ = 0.9 and the symmetric coupling, Γ l = Γ r .Each curve corresponds to the different temperature k B T /ħ hΓ whose value is annotated.
) as demonstrated in Figs. 6 (a) and (b).The Coulomb corrections in Y ℓ and Y l r (Y = L, K) are different [compare Eq. (

Figure 6 :
Figure 6: (a) The RC times τ K,l r (solid lines) and τ K,l (dotted lines) from thermal admittance and (b) the RC time τ L,l r (solid lines) and τ L,l (dotted lines) from thermoelectric admittance in the low-frequency limit as functions of the spin-degenerate dot level ε d in the interacting case with the same parameters as in Fig. 5. (c) max ε d [τ K,l r − τ K,l ] and (d) max ε d [τ L,l r − τ L,l ] as functions of the temperature.The dotted lines are asymptotes in the T → 0 limit.
, ω) vanish at λ = 1/2.The linear expansion is now applied to Ξ ab ℓ (n, n ′ , ω) [see Eq. (30)].Up to the linear order in Ψ ℓ , only the Fourier components with |n| ≤ 1, |n ′ | ≤ 1 and |n + n ′ | ≤ 1 are relevant: ) holds.As long as the condition (39) holds, the charge conservation condition can be used to remove dω Tr[G < (1, ω)] from the expression for the charge current.By inserting the expressions of the QD charge current (35) and the contact charge current (34) into the charge conservation (40), one can solve the integral dω Tr[G < (1, ω)] in terms of other QD Green's functions: with λ = 1/2, recover the non-interacting Dyson's equations (103) but with the self energy being now replaced by the Hartree one defined as Eq.(63).The Fourier components of the Hartree self energies in the linear response regime are then Σ R/A σ,HF (0, ω) = ±iΓ +