Study of the energy variation in many-body open quantum systems: role of interactions in the weak and strong coupling regimes

TL;DR

This paper derives an expression for energy variation in many-body open quantum systems, highlighting the role of interactions and virtual processes in energy exchange across leads, with implications for thermal conductance calculations.

Contribution

It introduces a new framework to analyze energy exchange in interacting quantum systems connected to multiple leads, emphasizing the impact of virtual processes.

Findings

Energy variation contributions from each lead can differ from the lead's energy change.

Discrepancies are due to virtual energy exchange among leads via the central region.

Implications for calculating thermal conductance in interacting systems.

Abstract

We derive an expression for the rate of change of the energy of an interacting many-body system connected to macroscopic leads. We show that the energy variation is the sum of contributions from each different lead. Unlike the charge current each of this contribution can differ from the rate of change of the energy of the lead. We demonstrate that the discrepancy between the two is due to the direct exchange of energy among the considered lead and all other leads. We conclude that the microscopic mechanism behind it are virtual processes via the interacting central region. We also speculate on what are the implications of our findings in the calculation of the thermal conductance of an interacting system.

Equations

$${\mathcal{J}}_C^{(\alpha )}(t)=2\text{Re}{\left\{\int d{x}_{1}d\bar{1}d\bar{2}H(1,\bar{1})G(\bar{1},\bar{2}){\Sigma }_{\alpha }(\bar{2},{\bar{1}}^{+})\right\}}_{z_1=t}$$

$$\frac{d}{dt}⟨{\hat{V}}_C^{(\alpha )}(t)⟩=-{\mathcal{J}}_{\alpha }(t)-{\mathcal{J}}_C^{(\alpha )}(t)-\Delta {\mathcal{J}}_{\alpha }(t)+\frac{\partial }{\partial t}⟨{\hat{V}}_C^{(\alpha )}(t)⟩,$$

$${\dot{\mathcal{E}}}_C\equiv \frac{d}{dt}⟨{\hat{H}}_C(t)⟩=i⟨[\hat{H}(t),{\hat{H}}_C(t)]⟩$$

$${\dot{\mathcal{E}}}_C(z_1)=i⟨[\hat{H}(z_1),{\hat{H}}_D(z_1)]⟩=i⟨[{\hat{H}}_T(z_1),{\hat{H}}_D(z_1)]⟩$$

$$\displaystyle=\sum\limits_{\alpha}\biggl{(}i\left\langle\left[\int d\mathbb{x}_{1}\left(\hat{\psi}^{\dagger}(1)T(1)\hat{\psi}_{\alpha}(1)+\hat{\psi}_{\alpha}^{\dagger}(1)T^{*}(1)\hat{\psi}(1)\right),\int d\mathbb{x}_{2}\hat{\psi}^{\dagger}(2)h(2)\hat{\psi}(2)\right]\right\rangle+$$

$$\displaystyle+i\biggl{\langle}\biggl{[}\int d\mathbb{x}_{1}\left(\hat{\psi}^{\dagger}(1)T(1)\hat{\psi}_{\alpha}(1)+\hat{\psi}_{\alpha}^{\dagger}(1)T^{*}(1)\hat{\psi}(1)\right),\frac{1}{2}\int_{\gamma}dz_{1}^{\prime}\int\int d\mathbb{x}_{2}d\mathbb{x}_{1}^{\prime}\hat{\psi}^{\dagger}(2)\hat{\psi}^{\dagger}(1^{\prime})v(1^{\prime},2)\hat{\psi}(1^{\prime})\hat{\psi}(2)\biggr{]}\biggr{\rangle}\biggr{)}$$

$$\displaystyle=\sum\limits_{\alpha}\biggl{(}-i\int d\mathbb{x}_{1}\;T(1)h(1)\left\langle\hat{\psi}^{\dagger}(1)\hat{\psi}_{\alpha}(1)\right\rangle+i\int d\mathbb{x}_{1}\;T^{*}(1)h(1)\left\langle\hat{\psi}_{\alpha}^{\dagger}(1)\hat{\psi}(1)\right\rangle+$$

$$-\frac{i}{2}{\int }_{\gamma }d{z}_1^{\prime }\int \int d{x}_{1}d{x}_1^{\prime }T(1^{\prime })v(1,1^{\prime })(-i{)}^2⟨{\hat{\psi }}^{†}(1){\hat{\psi }}^{†}(1^{\prime })\hat{\psi }(1){\hat{\psi }}_{\alpha }(1^{\prime })⟩+$$

$$+\frac{i}{2}{\int }_{\gamma }d{z}_1^{\prime }\int \int d{x}_{1}d{x}_1^{\prime }T(1)v(1,1^{\prime })(-i{)}^2⟨{\hat{\psi }}^{†}(1){\hat{\psi }}^{†}(1^{\prime })\hat{\psi }(1^{\prime }){\hat{\psi }}_{\alpha }(1)⟩+$$

$$-\frac{i}{2}{\int }_{\gamma }d{z}_1^{\prime }\int \int d{x}_{1}d{x}_1^{\prime }{T}^{\ast }(1)v(1,1^{\prime })(-i{)}^2⟨{\hat{\psi }}_{\alpha }^{†}(1){\hat{\psi }}^{†}(1^{\prime })\hat{\psi }(1^{\prime })\hat{\psi }(1)⟩+$$

$$\displaystyle+\frac{i}{2}\int_{\gamma}dz_{1}^{\prime}\int\int d\mathbb{x}_{1}d\mathbb{x}_{1}^{\prime}\;T^{*}(1^{\prime})v(1,1^{\prime})(-i)^{2}\left\langle\hat{\psi}^{\dagger}_{\alpha}(1^{\prime})\hat{\psi}^{\dagger}(1)\hat{\psi}(1^{\prime})\hat{\psi}(1)\right\rangle\biggr{)}$$

$$\int d\mathbb{x}_{1}\;h(1)\left[T^{*}(1)i\left\langle\hat{\psi}_{\alpha}^{\dagger}(1)\hat{\psi}(1)\right\rangle-T(1)i\left\langle\hat{\psi}^{\dagger}(1)\hat{\psi}_{\alpha}(1)\right\rangle\right]=2Re\biggl{\{}\int d\mathbb{x}_{1}\;h(1)T^{*}(1)G_{C\alpha}(1;1^{+})\biggr{\}}$$

$$i{\int }_{\gamma }d{z}_1^{\prime }\int \int d{x}_{1}d{x}_1^{\prime }{T}_{\alpha }^{\ast }(1^{\prime })v(1,1^{\prime })(-i{)}^2⟨{\hat{\psi }}_{\alpha }^{†}(1^{\prime }){\hat{\psi }}^{†}(1)\hat{\psi }(1^{\prime })\hat{\psi }(1)⟩+h.c.$$

$$\displaystyle=2Re\biggl{\{}i\int_{\gamma}dz_{1}^{\prime}\int\int d\mathbb{x}_{1}d\mathbb{x}_{1}^{\prime}\;v(1,1^{\prime})T_{\alpha}^{*}(1^{\prime})G^{(2)}_{CCC\alpha}(1^{\prime},1;1^{+},1^{\prime+})\biggr{\}}$$

$$\mathcal{\dot{E}}_{C}=\sum\limits_{\alpha}2Re\biggl{\{}\int d\mathbb{x}\;h(1)T^{*}(1)G_{C\alpha}(1^{+})+i\int_{\gamma}dz^{\prime}\int\int d\mathbb{x}d\mathbb{x^{\prime}}\;v(1;1^{\prime})T_{\alpha}^{*}(1^{\prime})G^{(2)}_{CCC\alpha}(1^{\prime},1;1^{+},1^{\prime+})\biggr{\}}.$$

$$G_{C\alpha }(1;1^{\prime })=-i⟨{\mathcal{T}}_{\gamma }\hat{\psi }(1){\hat{\psi }}_{\alpha }^{†}(1^{\prime })⟩=-i⟨{\mathcal{T}}_{\gamma }\tilde{\psi }(1){\tilde{\psi }}_{\alpha }^{†}(1^{\prime })S⟩$$

$$S=\sum _{k=0}^{\infty }\frac{(-i{)}^k}{k!}{\int }_{\gamma }d{\bar{z}}_1\cdots {\int }_{\gamma }d{\bar{z}}_k{\tilde{H}}_T(\bar{1})\dots {\tilde{H}}_T(\bar{k})$$

$$-i⟨{\mathcal{T}}_{\gamma }\tilde{\psi }(1){\tilde{\psi }}_{\alpha }^{†}(1^{\prime })\sum _{k=0}^{\infty }\frac{(-i{)}^k}{k!}{\int }_{\gamma }\int d{\bar{z}}_{1}d{\bar{x}}_1\left({\tilde{\psi }}^{†}(\bar{1})T(\bar{1}){\tilde{\psi }}_{\alpha }(\bar{1})+h.c.\right)\times \cdots \times {\int }_{\gamma }d{\bar{z}}_k{\tilde{H}}_T(\bar{k})⟩=$$

$$=\int d\bar{1}\sum _{k=0}^{\infty }\frac{(-i{)}^k}{k!}(-i)⟨{\mathcal{T}}_{\gamma }{\tilde{\psi }}_{\alpha }(\bar{1}){\tilde{\psi }}_{\alpha }^{†}(1^{\prime })⟩T(\bar{1})⟨{\mathcal{T}}_{\gamma }\tilde{\psi }(1){\tilde{\psi }}^{†}(\bar{1})\times \cdots \times {\int }_{\gamma }d{\bar{z}}_k{\tilde{H}}_T(\bar{k})⟩+$$

$$\displaystyle+\text{$\big{(}$$

k-1$remainingterms$\big{)}

$$}$$

$

$$G_{C\alpha }(1;1^{\prime })=\int d\bar{1}(-i)⟨{\mathcal{T}}_{\gamma }\tilde{\psi }(1){\tilde{\psi }}^{†}(\bar{1})\sum _{k=0}^{\infty }\frac{(-i{)}^{k-1}}{(k-1)!}\prod _{p=1}^k{\int }_{\gamma }d{\bar{z}}_p{\tilde{H}}_T(\bar{p})⟩\times T_{\alpha }(\bar{1})(-i)⟨{\mathcal{T}}_{\gamma }{\tilde{\psi }}_{\alpha }(\bar{1}){\tilde{\psi }}_{\alpha }^{†}(1^{\prime })⟩$$

$$=\int d\bar{1}G(1;\bar{1})T_{\alpha }(\bar{1})g_{\alpha }(\bar{1};1^{\prime })$$

$$\displaystyle\sum\limits_{\alpha}2Re\biggl{\{}\int d\mathbb{x}_{1}h(1)T_{\alpha}^{*}(1)G_{C\alpha}(1;1^{+})\biggr{\}}=\sum\limits_{\alpha}2Re\biggl{\{}\int d\mathbb{x}_{1}d\bar{1}\;h(1)G(1;\bar{1})\Sigma_{\alpha}(\bar{1};1^{+})\biggr{\}}$$

$$G_{CCC\alpha }^{(2)}(1^{\prime },1;1^{+},1^{\prime +})=(-i{)}^2⟨{\mathcal{T}}_{\gamma }\hat{\psi }(1^{\prime })\hat{\psi }(1){\hat{\psi }}_{\alpha }^{†}(1^{\prime +}){\hat{\psi }}^{†}(1^{+})⟩=(-i{)}^2⟨{\mathcal{T}}_{\gamma }\tilde{\psi }(1^{\prime })\tilde{\psi }(1){\tilde{\psi }}_{\alpha }^{†}(1^{\prime +}){\tilde{\psi }}^{†}(1^{+})S⟩$$

$$\displaystyle=(-i)^{2}\biggl{\langle}\mathcal{T}_{\gamma}\tilde{\psi}(1^{\prime})\tilde{\psi}(1)\tilde{\psi}_{\alpha}^{\dagger}(1^{\prime+})\tilde{\psi}^{\dagger}(1^{+})\times\sum_{k=0}^{\infty}\frac{(-i)^{k}}{k!}\int_{\gamma}\int d\bar{1}\left(\tilde{\psi}^{\dagger}(\bar{1})T(\bar{1})\tilde{\psi}_{\alpha}(\bar{1})+h.c.\right)\times\dots\times\int_{\gamma}d\bar{z}_{k}\tilde{H}_{T}(\bar{k})\biggr{\rangle}$$

$$=\int d\bar{1}(-i{)}^2⟨{\mathcal{T}}_{\gamma }\tilde{\psi }(1^{\prime })\tilde{\psi }(1){\tilde{\psi }}^{†}(\bar{1}){\tilde{\psi }}^{†}(1^{+})\sum _{k=0}^{\infty }\frac{(-i{)}^k}{k!}\prod _{p=1}^k{\int }_{\gamma }d{z}_p{\tilde{H}}_T(\bar{p})⟩\times (-i)⟨{\mathcal{T}}_{\gamma }{\tilde{\psi }}_{\alpha }(\bar{1}){\tilde{\psi }}_{\alpha }^{†}(1^{\prime +})⟩T(\bar{1})$$

$$=\int d\bar{1}{G}^{(2)}(1^{\prime },1;1^{+},\bar{1})T(\bar{1})g_{\alpha }(\bar{1};1^{\prime +})$$

$$i\sum _{\alpha }\int d{x}^{\prime }d1v(1^{\prime },1)T^{\ast }(1^{\prime })G_{CCC\alpha }^{(2)}(1^{\prime },1;1^{+},1^{\prime +})=i\sum _{\alpha }\int d{x}^{\prime }d1d\bar{1}v(1^{\prime },1)T^{\ast }(1^{\prime })G^{(2)}(1^{\prime },1;1^{+},\bar{1})T(\bar{1})g_{\alpha \alpha }(\bar{1};1^{\prime +})$$

$$=-i\sum _{\alpha }\int d{x}^{\prime }d1d\bar{1}v(1^{\prime },1)G^{(2)}(1^{\prime },1;\bar{1},1^{+})T(\bar{1})g_{\alpha }(\bar{1};1^{\prime +})T^{\ast }(1^{\prime })=\sum _{\alpha }\int d{x}^{\prime }d1d\bar{1}\Sigma (1^{\prime };1)G(1;\bar{1}){\Sigma }_{\alpha }(\bar{1};1^{\prime +})$$

$$\int d\bar{1}v(1,\bar{1})G^{(2)}(1,\bar{1};1^{\prime },{\bar{1}}^{+})=i\int d\bar{1}\Sigma (1,\bar{1})G(\bar{1},1^{\prime }).$$

$$\frac{d}{dt}\langle\hat{H}_{C}(t)\rangle=\sum\limits_{\alpha}\mathcal{J}_{C}^{(\alpha)}(t)\equiv\sum\limits_{\alpha}2Re\biggl{\{}\int d\mathbb{x}_{1}d\bar{1}d\bar{2}\biggl{[}h(1)\delta(1,\bar{1})+\Sigma(1;\bar{1})\biggr{]}G(\bar{1};\bar{2})\Sigma_{\alpha}(\bar{2},1^{+})\biggr{\}}_{z_{1}=t}$$

$$\frac{d}{dt}\mathcal{E}_{C}(t)=\frac{d}{dz_{1}}\left[i\int d\mathbb{x}_{1}d\bar{1}\;\biggl{[}h(1)\delta(1,\bar{1})+\frac{1}{2}\Sigma(1;\bar{1})\biggr{]}G(\bar{1};1^{+})\right]\bigg{|}_{z_{1}=t}.$$

$$\frac{d}{dt}⟨{\hat{V}}_C^{(\alpha )}(t)⟩=2Re{\left\{-i\int d{x}_1{T}_{\alpha }^{\ast }(1)\frac{d}{d{z}_1}{G}_{C\alpha }(1;1^{+})\right\}}_{z_1=t}=$$

$$\displaystyle 2Re\left\{\int d\mathbb{x}_{1}d\bar{1}\;G(1;\bar{1})\Sigma_{\alpha}(\bar{1};1^{+})h_{\alpha}(1)\right\}_{z_{1}=t}-2Re\biggl{\{}\int d\mathbb{x}_{1}d\bar{1}d\bar{2}\biggl{[}h(1)\delta(1,\bar{1})+\Sigma(1;\bar{1})\biggr{]}G(\bar{1};\bar{2})\Sigma_{\alpha}(\bar{2},1^{+})\biggr{\}}_{z_{1}=t}$$

$$-$$

$$i\frac{d}{dz}G(1;1^{\prime })=\delta (1,1^{\prime })+h(1)G(1;1^{\prime })+\int d\bar{1}\left(\Sigma (1;\bar{1})+{\Sigma }_{emb}(1;\bar{1})\right)G(\bar{1};1^{\prime })$$

$$i\frac{d}{d{z}^{\prime }}{g}_{\alpha }(1;1^{\prime })=-\delta (1,1^{\prime })-g_{\alpha }(1;1^{\prime })h(1^{\prime }),$$

$$G_{\beta \alpha }(1;1^{\prime })\equiv \int d\bar{1}d\bar{2}(-i)⟨{\mathcal{T}}_{\gamma }{\tilde{\psi }}_{\beta }(1){\tilde{\psi }}_{\beta }^{†}(\bar{1})⟩T_{\beta }^{\ast }(\bar{1})T_{\alpha }(\bar{2})(-i)⟨{\mathcal{T}}_{\gamma }{\tilde{\psi }}_{\alpha }(\bar{2}){\tilde{\psi }}_{\alpha }^{†}(1^{\prime })⟩$$

$$\times (-i)⟨{\mathcal{T}}_{\gamma }\tilde{\psi }(\bar{1}){\tilde{\psi }}^{†}(\bar{2})\sum _{k=0}^{\infty }\frac{(-i{)}^k}{k!}\prod _{p=1}^k{\int }_{\gamma }d{z}_p{\tilde{H}}_T(\bar{p})⟩$$

$$=\int d\bar{1}d\bar{2}{g}_{\beta }(1;\bar{1})T_{\beta }^{\ast }(\bar{1})G(\bar{1};\bar{2})T_{\alpha }(\bar{2})g_{\alpha }(\bar{2};1^{\prime })$$

Full text

I Introduction

Understanding how quantum systems conduct energy is currently at the centre of investigation in different branches of physics. Recent advances in the manipulation of nanoscale systems jukka2009 have given a new boost to theoretical investigations in that direction; nevertheless how energy and most intriguingly heat propagates through a medium has always been at the centre of theoretical debate. A better understanding of the energy exchange mechanism is also crucial to design devices which can act as (local) refrigerators mikko2018 , energy harvesters jaliel2019 or can be used as a heat valve jukka2018 . The study of energy currents at nanoscale also helps in the characterization of certain physical effects such as Seebeck effect svilans2018 ; clemens2019 , violation of the Wiedemann-Franz law bivas2017 and rectification of currents which are useful in applications at nanoscales giazotto2015 . Another interesting problem is understanding the connection between the diffusion properties in a closed system and its conduction properties when it is connected to external leads (open system) dhar2018 ; dhar2019 . This question also arises in the case of transport in classical systems li2003 because it is known that the transport properties are modified by the nature of the scattering and the type of reservoirs chosen.

The Landauer-Büttiker formula for the particle current through a noninteracting system is a milestone in transport theory landauer1957 ; buttiker1986 . It was readily extended to the case of interacting systems by Meir and Wingreen Meir1992 ; Jauho1994 giving the possibility of studying the effects of interactions on the transport properties of many-body quantum systems. It relates the current to the spectral features of the system, therefore giving a microscopic picture of the mesoscopic transport properties. Specifically, in the original work the transport features of the single-impurity Anderson model (SIAM) are linked to the spectral function of the system coupled to the leads in both the Coulomb blockade and Kondo regimes. The Meir-Wingreen expression, which holds in the stationary state, has also been extended to transient times allowing us to study transport in driven systems riku2013 ; Haugbook ; stefleebook ; riku2016 ; riku2018 . The energy current is defined as the variation of the Hamiltonian of the lead; in this case one obtains an expression analogous to the Meir-Wingreen one for particle current, which holds for both interacting and non-interacting systems provided that the leads are assumed to be non-interacting. In the case of transport through an interacting central region, these expressions correctly capture the many-body interactions through the transmission function or alternatively through the spectral density of the central correlated region wang2014 . A strictly related issue is the definition of heat current, namely the contribution to energy current which is irreversibly dissipated into the lead and that is consistent with the laws of thermodynamics. The bone of contention is the coupling energy term which needs to be assigned either to the system, or to the terminal, or be split among the two esposito2015 ; misha2014 ; misha2016 ; misha2018 ; eich2014 ; eich2016 ; lopez2015 . The debate about the definition of heat current, although interesting and very important, goes beyond the purpose of this work.

In our work, we are interested in the energy variation of the central correlated region due to the reservoirs and we study the difference of this quantities with the energy variation of the reservoirs itself. We use the the non-equilibrium Green’s function formalism which naturally allows us to account for both many-body interactions and coupling to external reservoirs through the self-energy functional. We show that the two energy variation need not to be equal (and opposite in sign) as it is for particle currents; we discuss this discrepancy by looking at a specific example: the single-impurity Anderson model.

II Non-equilibrium Green’s function approach

We consider a correlated central region where particles can interact through a two-body potential connected to $N$ leads where particles are assumed to be non-interacting. This is the only assumption we make and it is the same needed to derive the Meir-Wingreen expression (and its time dependent version). A schematic representation of the system we consider is shown in Fig. 1 panel $a)$. The Hamiltonian of the total system is given by $\hat{H}(z)=\hat{H}_{C}(z)+\sum_{\alpha=1}^{N}\hat{H}_{\alpha}(z)+\sum_{\alpha=1}^{N}\hat{V}_{C}^{(\alpha)}(z)$, where $z$ is a generic complex time variable on the Keldysh contour $\gamma$ Haugbook ; stefleebook ; bonitzbook ; Keldysh1965 (Fig. 1 panel $b)$). In the central region the Hamiltonian reads $\hat{H}_{C}(z)=\int d\mathbb{x}_{1}\;\hat{\psi}^{\dagger}(1)h(1)\hat{\psi}(1)+\frac{1}{2}\int d\mathbb{x}_{1}d1^{\prime}\;\hat{\psi}^{\dagger}(1)\hat{\psi}^{\dagger}(1^{\prime})v(1,1^{\prime})\hat{\psi}(1^{\prime})\hat{\psi}(1)$.

The indices $1=(\mathbb{x}_{1},z_{1})$, $1^{\prime}=(\mathbb{x}_{1}^{\prime},z_{1}^{\prime})$ are collective indices for the position-spin coordinate $\mathbb{x}=(\mathbb{r},\sigma)$ and complex time $z$, $h(1)$ the single particle Hamiltonian in the central region and $v(1,2)=\delta_{\gamma}(z_{1}-z_{2})v(\mathbb{x}_{1},z_{1};\mathbb{x}_{2},z_{2})$ a generic two-body interaction. $\hat{H}_{\alpha}(z)=\int d\mathbb{x}_{1}\hat{\psi}^{\dagger}_{\alpha}(1)h_{\alpha}(1)\hat{\psi}_{\alpha}(1)$ describes particles in the region of the $\alpha$-th lead with $h_{\alpha}(1)$ the single particle Hamiltonian. The Hamiltonian accounting for the coupling between the interacting region and the leads is chosen to be tunnel-like and given by $\hat{V}_{C}^{(\alpha)}(z)=\int d\mathbb{x}_{1}\left(\hat{\psi}^{\dagger}(1)T_{\alpha}(1)\hat{\psi}_{\alpha}(1)+\text{h.c.}\right)$, where $T(1)$ are the tunneling energies between particles in the interacting region and those in the $\alpha$-th lead. The key object in the NEGF framework is the single-particle Green’s function (SPGF) $G(1,1^{\prime})=-i\left\langle\mathcal{T}_{\gamma}\hat{\psi}(1)\hat{\psi}^{\dagger}(1^{\prime})\right\rangle_{0}$ with $\mathcal{T}_{\gamma}$ the time-ordering operator over the Keldysh contour $\gamma$ and $\left\langle\dots\right\rangle_{0}$ the average over the initial many-body state. The knowledge of $G$ gives access to all single particle quantities (e.g. density, momentum distribution, density-of-states). For an interacting system the SPGF can be found by solving either the Kadanoff-Baym equations stefleebook ; bonitzbook ; robert2009 ; robert2016 , the Dyson equation logullo2016 ; talarico2019 , or other equivalent techniques spicka1986 ; bonitz2012 ; eckstein2014 ; daniel2018 . A crucial point is the choice of a self-energy, which is typically a functional of the SPGF itself. The reason to introduce the self-energy is to truncate the hierarchy of equations which couples the equation-of-motion for the SPGF to those of higher-order Green’s functions stefleebook . Henceforth, we shall assume that a self-energy, accounting for both the many-body interactions and the coupling to the leads, has been chosen and that a solution for the SPGF has been found. Under these assumptions we shall compute the energy variation in the interacting region in a way which is consistent with the approximations embodied in the chosen self-energy.

III Energy variation in the central region

We start considering the rate of variation of the energy in the central (interacting) region $d\langle\hat{H}_{C}(t)\rangle/dt=i\langle[\hat{H}(z),\hat{H}_{C}(z)]\rangle_{z=t}+\partial\langle\hat{H}_{C}(z)\rangle/\partial z|_{z=t}\equiv\mathcal{J}(t)+P(t)$ where we defined the energy current $\mathcal{J}(t)$ and the power $P(t)$ due to external drive; in the following we will disregard this last term as it is not relevant for the purposes of our discussion. By computing the commutator of $\hat{H}_{C}(z)$ with the total Hamiltonian (see Appendix A for details), we obtain an expression in terms of $G_{C\alpha}(1,1^{\prime})=-i\langle\mathcal{T}_{\gamma}\hat{\psi}(1)\hat{\psi}_{\alpha}^{\dagger}(1^{\prime})\rangle$ and $G^{(2)}_{CCC\alpha}(1,2;3,4)=(-i)^{2}\langle\mathcal{T}_{\gamma}\hat{\psi}(1)\hat{\psi}(2)\hat{\psi}_{\alpha}^{\dagger}(4)\hat{\psi}^{\dagger}(3)\rangle$.

This expression can be manipulated by means of the S-matrix expansion in the interaction picture with respect to $\hat{H}_{C}(t)+\sum\limits_{\alpha}\hat{H}_{\alpha}(t)$ to obtain $\mathcal{J}(t)=\sum\limits_{\alpha}\mathcal{J}_{C}^{(\alpha)}(t)$ (see Appendix B for details) with:

[TABLE]

where $H(1,1^{\prime})=\left[h(1)\delta(1,1^{\prime})+\Sigma(1,1^{\prime})\right]$ and $1^{+}=(\mathbb{x}_{1},z_{1}^{+})$, with $z_{1}^{+}$ a time infinitesimally later than $z_{1}$ on the Keldysh-contour.

Eq.  (1) expresses the variation of the energy in the interacting region as sum of different contributions each coming from the individual leads. The single contribution from the lead $\alpha$ is in turn made of two terms; the first one is similar to the expression of the particles’ current due to the $\alpha$-th lead: $\mathcal{I}_{C}^{(\alpha)}(t)=2\text{Re}\left\{\int d\mathbb{x}_{1}d\bar{1}d\bar{2}\;G(\bar{1},\bar{2})\Sigma_{\alpha}(\bar{2},\bar{1}^{+})\right\}_{z_{1}=t}$, with the significant difference that it contains the single-particle Hamiltonian $h(1)$ and it accounts for the energy carried by the flowing particles; the effect of the interactions is present only in the SPGF $G$ through a change in the density of states of the particles in the central region. This term is analogous to the energy variation of the lead where the single particle Hamiltonian of the central region is replaced by the single particle Hamiltonian of the $\alpha$-th lead. The second term contains the many-body self-energy explicitly and therefore accounts for the transport of the particle-particle interaction energy; it can be seen as a mechanism of energy redistribution in the central region due to particle-particle scattering: non-interacting electrons coming from the leads scatter in the central region releasing part of their energy before tunneling into a new non-interacting state of another lead.

IV Inter-leads coupling

The expression in Eq. (1) has to be compared with the Meir-Wingreen like expression $J_{\alpha}(t)=i\langle[\hat{H}(z),\hat{H}_{\alpha}(z)]\rangle_{z=t}$, that describes the energy variation of the lead $\alpha$ due to the coupling with the central correlated region. For the particle currents, it is true that the current into the lead is equal and opposite to the charge variation in the central system due to that lead, namely $\mathcal{I}_{\alpha}=-\mathcal{I}_{C}^{(\alpha)}$. Here we defined the current in the $\alpha$ lead as $\mathcal{I}_{\alpha}=i\langle[\hat{H}(z),\hat{N}_{\alpha}(z)]\rangle_{z=t}$ and the currents into the central region due to the lead $\alpha$ as $\sum\limits_{\alpha}\mathcal{I}_{C}^{(\alpha)}=i\langle[\hat{H}(z),\hat{N}_{C}(z)]\rangle_{z=t}$. This relations translates into formal mathematical expression the intuitive concept of the locality of the number operator. Basically, they state that particles flowing out of a lead necessarily enter the central region causing an equally and opposite variation. In analogy with the particle current we would expect that the following relation for the energy current $\mathcal{J}_{\alpha}(t)=-\mathcal{J}_{C}^{(\alpha)}(t)$ should always be true as well. However, as we prove in the following, this is not always the case. It is sufficient to look at the time derivative $d\langle\hat{V}_{C}^{(\alpha)}(t)\rangle/dt$ of the coupling Hamiltonian between the terminal $\alpha$ and the central region. It is easy to prove (see Appendix C for details) that

[TABLE]

where $\Delta J_{\alpha}(t)=\sum\limits_{\beta\neq\alpha}\int d\mathbb{x}_{1}T_{\beta}(1)G_{\beta\alpha}(1;1^{+})T_{\alpha}^{*}(1)$ and $G_{\beta\alpha}(1,1^{\prime})=-i\langle\mathcal{T}_{\gamma}\hat{\psi}_{\beta}(1)\hat{\psi}_{\alpha}^{\dagger}(1^{\prime})\rangle$ is the lead $\beta$- lead $\alpha$ Green’s function. It has been shown misha2016 that the DC component of $\frac{d}{dt}\langle\hat{V}_{C}^{(\alpha)}(t)\rangle=0$ and therefore $\mathcal{J}_{\alpha}(t)=-\mathcal{J}_{C}^{(\alpha)}(t)-\Delta\mathcal{J}_{\alpha}(t)$. Therefore, we have shown that a variation of the energy of the lead is not necessarily accompanied by an equal and opposite change in the energy of the central region as in the case for the particle variation. The difference between the two is given by the term $\Delta\mathcal{J}_{\alpha}(t)$ which arises from the direct propagation of particles from any other lead $\beta(\neq\alpha)$ and the lead $\alpha$ via virtual scattering through the central region. This interpretation is supported by the physical meaning of $G_{\beta\alpha}(1,1^{\prime})=\int d\bar{1}d\bar{2}\;g_{\beta}(1;\bar{1})T_{\beta}^{*}(\bar{1})G(\bar{1};\bar{2})T_{\alpha}(\bar{2})g_{\alpha}(\bar{2};1^{\prime})$ (see Appendix C for details): it represents the propagation of a particle from one lead to another through the central region. The expression $\Sigma_{\beta\alpha}(1;1^{\prime})=T_{\beta}^{*}(1)G(1;1^{\prime})T_{\alpha}(1^{\prime})$ is also called the inbedding self-energy stefleebook and accounts for the back-action of the central region onto the leads. This term is usually either negligible or zero altogether. In the weak coupling limit this term is vanishingly small due to the fact that it is of fourth order in the coupling between the central region and the leads. Moreover it also vanishes when the two integrals in the expression for $G_{\beta\alpha}(1;1^{\prime})$ have different spatial supports, namely when $T_{\alpha}(1)$ and $T_{\beta}(1^{\prime})$ for $\alpha\neq\beta$ are non-zero on different spatial regions and therefore their product vanishes.

Nonetheless, it is possible to envisage situations in which this term gives a non-vanishing contribution: short quantum wires (mean free path comparable with the size of the wire), more than one lead coupled to the same spatial region or spatially extended coupling between the leads and the central region. As a consequence, the energy variation of the lead is not only related to the one of the central region, as one would expect, but also accounts for contribution coming from the direct exchange of energy with other leads. Therefore, as one can easily deduce, it is crucial to consider this term in situations where one wants to infer the thermal transport properties of a system (the central region) by measuring the properties of the reservoirs (e.g. particle distribution, temperatures).

V SIAM in the Kondo regime

In order to show the consequences of the results presented so far we look at the energy-transport in the single-impurity Anderson model anderson ; glazman2001 described by the Hamiltonian $\hat{H}=\varepsilon\sum_{\sigma}\hat{d}_{\sigma}^{\dagger}\hat{d}_{\sigma}+U\hat{d}_{\uparrow}^{\dagger}\hat{d}_{\uparrow}\hat{d}_{\downarrow}^{\dagger}\hat{d}_{\downarrow}+\sum_{\alpha,k\sigma}(\varepsilon_{\alpha k\sigma}+\mu_{\alpha})\hat{c}_{\alpha k\sigma}^{\dagger}\hat{c}_{\alpha k\sigma}-g\sum_{\alpha,k\sigma}(\hat{c}_{\alpha k\sigma}\hat{d}_{\sigma}^{\dagger}+h.c.)$. Here, $\hat{d}_{\sigma}^{\dagger}\left(\hat{d}_{\sigma}\right)$ corresponds to creation (annihilation) of electron on the impurity level with spin $\sigma$, $\varepsilon$ denotes the single-particle energies, and $U$ is the electronic charging energy. The operator $\hat{c}_{\alpha k\sigma}^{\dagger}\left(\hat{c}_{\alpha k\sigma}\right)$ creates (annihilates) electron with state $k$ and spin $\sigma$ in the lead $\alpha=L,R$ with chemical potential $\mu_{\alpha}$. Finally, $g$ denotes the tunneling amplitude between the terminal and the impurity level. In the wide-band limit approximation we are left with a frequency-independent coupling $\Gamma=|g|^{2}$. In what follow we consider unbiased voltage ($\mu=\mu_{L}=\mu_{R}=0$), a finite and symmetric thermal bias $T_{L}=T+\Delta T/2$, $T_{R}=T-\Delta T/2$, with $T=(T_{L}+T_{R})/2=0.125$ and $\Delta T=T_{L}-T_{R}=0.15$ and we study the dynamical quantities as function of the impurity level position $\varepsilon$. We present our results as a function of the shifted single-particle energies or gate-voltage $v_{g}=\varepsilon+U/2$, such that the particle-hole symmetric point $\varepsilon=-U/2$ correspond to $v_{g}=0$.

Under certain conditions the SIAM may reveal Kondo effects characterized by the formation of a correlated asymmetric resonance in the spectral structure of the system haldane ; glazman2001 . The conditions to attain such regime are a large charging energy and tunnel coupling at low temperature, $U\gg\Gamma\gg T$ and $T\ll T_{K}$. The expression for the threshold temperature or Kondo temperature is given by $T_{K}=\frac{1}{2}\sqrt{\Gamma U}exp\left(\pi\frac{(v_{g}^{2}-U^{2}/4)}{\Gamma U}\right)$ hewson and is strictly valid in the Kondo regime where $-U/2+\Gamma/2<v_{g}<U/2-\Gamma/2$.

The interactions in the correlated impurity are treated in the self-consistent $GW$ approximation (see Appendix D for details), which on the one hand is guaranteed to satisfy macroscopic conservation laws stefleebook , and on the other hand allows us to explore certain features of the Kondo regime rubio2008 . We solve the Dyson equation numerically using the method described in Ref. logullo2016 ; talarico2019 . We checked that the obtained expression for the total current through the interaction region $\mathcal{J}(t)$ in Eq. (1) is consistent with the derivative of the total energy in the corresponding region $d\mathcal{E}(t)/dt$ where \mathcal{E}(t)=\langle\hat{H}_{C}(t)\rangle=\left[i\int d\mathbb{x}_{1}d2\;\biggl{[}h(1)\delta(1,2)+\frac{1}{2}\Sigma(1,2)\biggr{]}G(2,1^{+})\right]. We find, within numerical errors, a remarkable agreement between these two quantities at any time of the dynamics and for any set of the parameters considered (Appendix D). This demonstrates that the approximations done on the two-particle Green’s function, while working out an expression for $\mathcal{J}(t)$, are indeed consistent with those used to solve the Dyson equation for the single-particle Green’s function and embodied into the chosen self-energy.

Before moving to the systematic study of the currents, we first characterize the regimes of parameters where the Kondo effects may occur. Fig. 2 panel $a)$ shows a map in the gate-voltage/frequency plane of the non-equilibrium spectral function of the system $\mathcal{A}(\omega)=i[G^{R}(\omega)-G^{A}(\omega)]$ obtained from the solution of the Dyson equation in a strong coupling regime with the leads $\Gamma=1.3$ and for different charging energies $U=0,2,4$. As the interaction increases, one can identify a range in the gate-voltage where the spectral function is centered at the chemical-potential of the leads, that correspond exactly to the range where Kondo correlations are expected $-U/2+\Gamma/2<v_{g}<U/2-\Gamma/2$ (white dashed lines). Fig. 2 panel $b)$ shows instead how the Kondo correlations reshape the spectral function at the charge degeneracy point $v_{g}=0$. Going from zero to large interactions, the Lorentzian shape is lost to leave a more asymmetric and pinned structure.

Fig. 2 panel $c)$ shows the Kondo temperature $T_{K}$ for the charging energy $U=4$ as a function of $v_{g}$ in the range of gate-voltage where Kondo correlations take place. As we can see the average temperature $T$ (red solid line) lies slightly above the Kondo temperature $T_{K}$ in a region close to the degeneracy point $v_{g}=0$. Nonetheless, the spectral function is still pinned at the chemical potential of the leads in the whole interval and its shape is very different from a Lorentzian, a clear signature of the Kondo regime. This may be explained from the highly out-of-equilibrium scenario ($T\sim\Delta T$) that we consider; it is quite reasonable that the system does not equilibrate to the average temperature $T$ of the leads and instead reaches a lower temperature that allows the Kondo correlations to persist.

VI Particle and energy currents

We now look at the steady-state properties of the time-dependent expressions for the particle and energy currents, namely the limit $t\to\infty$. In the case of two symmetrically coupled terminals the particle current is given by the Meir-Wingreen formula $\underset{t\rightarrow\infty}{\lim}\mathcal{I}_{\alpha}(t)\equiv\mathcal{I}_{\alpha}^{(S)}=\int d\omega\Gamma(\omega)[f_{\alpha}(\omega)-f_{\bar{\alpha}}(\omega)]\mathcal{A}(\omega)$. Here $\Gamma(\omega)=\Gamma_{L}(\omega)\Gamma_{R}(\omega)/[\Gamma_{L}(\omega)+\Gamma_{R}(\omega)]$, $\bar{\alpha}\neq\alpha$, $f_{\alpha}(\omega)=(1+e^{\beta_{\alpha}(\omega-\mu_{\alpha})})^{-1}$ is the Fermi-Dirac distribution function for electrons in the lead $\alpha$ and $\mathcal{A}(\omega)$ is the non-equilibrium spectral function of the central region. For the energy current we compare $\bar{\mathcal{J}}_{C}^{(L)}\equiv\underset{t\rightarrow\infty}{\lim}\mathcal{J}_{C}^{(L)}(t)$ with $\underset{t\rightarrow\infty}{\lim}\mathcal{J}_{L}(t)\equiv\mathcal{J}_{L}^{(S)}=\int d\omega\;\omega\Gamma(\omega)[f_{\alpha}(\omega)-f_{\bar{\alpha}}(\omega)]\mathcal{A}(\omega)$. The latest expression is the energy counterpart of the Meir-Wingreen formula for the particle current and properly describes how the energy runs across the non-interacting leads. It is widely used to characterize thermoelectric properties of correlated materials in linear response  kim2002 ; lei2002 ; costi2010 . Fig. 3 shows the agreement between the particle current in the central region due to the left lead $\mathcal{I}_{C}^{(L)}(t),(t\Gamma\gg 1)$ and the opposite of the current that flows into the same lead $-\mathcal{I}^{(S)}_{L}$ for both weak coupling $\Gamma\ll\Delta T$ panel $a)$ and strong coupling $\Gamma\gg\Delta T$ panel $b)$, for different values of the interaction strength $U$. The two currents are actually equal and opposite regardless the regime consider and for any range of parameters chosen as one would expected. We observe that for positive gate voltages particles flow from the central region to the left lead $\left(\mathcal{I}_{C}^{(L)}(t)=-\mathcal{I}^{(S)}_{L}<0\right)$ because the position of the spectral function of the system lies above the chemical potential of the lead ($\mu_{L}=0$). The opposite situation is realized for negative gate potentials where particles flow from the left lead to the central region $\left(\mathcal{I}_{C}^{(L)}(t)=-\mathcal{I}^{(S)}_{L}>0\right)$. In the weak-coupling case, panel Fig. 3 $a)$, the particle current is fully specified by the window function of the leads $\mathcal{I}_{L}\sim\Delta f(v_{g})=f_{L}(v_{g})-f_{R}(v_{g})$. This may be understood from the shape of the spectral function that, although being a Lorentzian function, is much sharper and narrow than the energy window of the leads ($\Gamma\ll\Delta T$) and it can be approximated with $\mathcal{A}(\omega)\sim\delta(\omega-v_{g})$. This picture is slightly different when the interaction is present since, in this case, the overall effect is to broaden the spectral function that cannot anymore be treated as a delta function. Nonetheless, the physical scenario is qualitatively the same and only quantitatively different. The strong-coupling case, Fig. 3 $b)$, does not present significant differences outside of the Kondo region. When the system enters the Kondo regime (green and gray shaded areas) the behavior of the current is drastically different from the window function $\Delta f(\omega)$ reflecting the change in the shape of the spectral function in this regime.

In contrast to the particle current, a perfect agreement between the energy current across the interacting region due to the left reservoir $\bar{\mathcal{J}}_{C}^{(L)}$ and the opposite current in the same lead $-\mathcal{J}_{L}^{(S)}$ is found only in the weak-coupling scenario as shown in Fig. 3 panel $c)$. At strong coupling we observe a qualitatively different behavior and not only a merely quantitative deviation between this two expressions, Fig. 3 panel $d)$, as it was for the particle currents. Interactions are not responsible for it as the effect is present, and actually more pronounced, in the non-interacting case $(U=0)$. The maximum difference between the two currents is at particle-hole symmetric point $(v_{g}=0)$ where $\bar{\mathcal{J}}_{C}^{(L)}$ is zero and $\mathcal{J}_{L}^{(S)}$ has its maximum. This can be explained with the following argument: the central region is at resonance, it is completely transparent and all the energy flows from one lead to the other without affecting the energy of the central region. The presence of Kondo correlations manifest themselves as a plateau in the energy currents in the range of gate-voltage where they are predicted, namely for $-U/2+\Gamma/2<v_{g}<U/2-\Gamma/2$. These regions are shown in Fig.3 as green and gray shaded areas and they corresponds to the regions in which the spectral function is pinned at the chemical potential of the leads, as shown by the horizontal solid lines in Fig.2 panel a). One interesting observation is that the extension of the plateaux of the lead energy current overestimates this regions, whereas the one in the energy variation of the central region shows a better agreement. As we noticed above the pinning of the spectral function is an hallmark of the emergence of the Kondo bound state and it can be probed through transport experiments martin2018 . In these experiments properties of the central region are inferred from measurements onto the leads. We have just shown that in the determination of the extent of the Kondo regime in gate voltage, it might be beneficial to consider this extra term to get a more precise estimation of the Kondo region for single level quantum dots in the strong coupling regime. Interestingly, in the case of a thermal bias in the leads the plateau induced by the Kondo cloud appears in the energy current, whereas in the case of small bias voltage it appears in the particle current as shown in Ref. talarico2019 . In general interactions increase the change of energy in the central region with respect to the non-interacting case. Nevertheless, there is a reduction between the energy current in the left lead $\mathcal{J}^{(S)}_{L}$ and the energy change in the central interacting region predicted by Eq. 1.

VII Virtual processes

In Sec.IV we have shown that the energy variation of one lead and the lead’s contribution to the variation of the the energy of the central region differ by $\Delta J_{\alpha}(t)$. This term has been shown to arise from the direct coupling of the lead to all other leads via second order processes through the central region. Moreover, we observe that in the strong coupling regime of the SIAM the two currents are indeed different and that, this effect does not come from the strong correlations present in this regime, but actually it is more pronounced in the non-interacting case. If we compute the term $-\mathcal{J}_{L}^{(S)}-\mathcal{J}_{C}^{(L)}$ and compare it with $\Delta\mathcal{J}_{L}$ computed separately we obtain a perfect agreement as shown in Fig. 3 $e)-f)$ for both weak and strong coupling respectively. Recently, this term has been shown josefsson2019 to be useful in computing the energy variation in the lead, that has been also measured in a setup analogous to the one we have considered. It is worth to notice that, whenever the interaction is present, the energy rate contribution $\mathcal{J}_{C}^{(L)}(t)$ is never zero when the density of states of the dot is in the energy window of the leads. Thus, the presence of the interaction redistributes the energy inside the dot and forbids from a perfect transparency as the one observed in the non-interacting case.

The lead-lead coupling has been already discussed in  wang2012 ; wang2014 with the remarkable difference that there a lead-lead interaction was present from the beginning in the Hamiltonian of the system. Moreover, the effect of such coupling on the energy transfer and its microscopic mechanism has not been discussed. As a concluding remark, we would like to highlight the importance and the implications of the term $\Delta\mathcal{J}_{\alpha}$ whenever it is not vanishing. When probing the conduction properties, either electrical or thermal, measurements are performed on the leads. In the case of particle currents the rate of change of particles in the lead equals the change of particles in the central region due to the lead itself. Therefore, measure of the charge current into the lead is equivalent to measure the change in the particle currents in the central region. Hence it is meaningful to infer the electrical conductivity of the central region from the measured current in the leads. Nevertheless, when it comes to compute the thermal conductivity of the system in the central correlated region $\kappa=\frac{\partial\dot{\mathcal{Q}}_{L}}{\partial\Delta T}\mid_{\mathcal{I}_{L}^{(S)}=0}$, with the heat current given by $\dot{\mathcal{Q}}_{L}=(\mathcal{J}_{L}^{(S)}-\mu_{L}\mathcal{I}_{L}^{(S)})$, one has to take extra care on what is the actual interpretation of what it is computed. Indeed, according to Eq.2, at stationarity and in the absence of external drive, one has $\kappa=-\frac{\partial}{\partial\Delta T}\left(\bar{\mathcal{J}}_{C}^{(L)}-\mu_{L}\mathcal{I}_{C}^{L}+\Delta\mathcal{J}_{L}\right)\mid_{\mathcal{I}_{C}^{(L)}=0}$, where we replace $\mathcal{I}_{L}^{(S)}\rightarrow-\mathcal{I}_{C}^{(L)}=-\underset{t\rightarrow\infty}{\lim}\mathcal{I}_{C}^{(L)}(t)$. The last equation shows that to compute the thermal conductance of the central region, the lead-lead term $\Delta\mathcal{J}_{L}$ needs to be considered. The term, as we discussed, accounts for the energy transfer through the contacts or terminals, thus in the strong coupling or for spatially extended couplings, what it is actually computed is the interface thermal conductance between tow leads, as if the two were directly coupled via a renormalized tunneling. We conclude that in the case of strong coupling and/or when the coupling of the leads are spatially extended, the usual definition of thermal conductance can be misleading as a figure of merit for the thermal properties of the central region. This is in contrast with the electrical conductance which is, instead, always consistent with the rate of change of the particle currents in the middle region. As a concluding remark, we would like to comment on the possibility to define, by analogy with the lead, the heat rate $\dot{\mathcal{Q}}_{C}^{(L)}=\bar{\mathcal{J}}_{C}^{(L)}-\mu_{L}\mathcal{I}_{C}^{L}$, namely the heat flowing from the central region to the left lead. However, this definition is in contrast with the thermodynamic formulation of heat as we will show in a forthcoming work.

VIII Conclusions

We have shown that the energy variation of an interacting system coupled to $N$ leads is equal to the sum of different contributions each related with the energy flowing from/to the individual leads. Each term is in turn made of two further contributions, a first one accounting for the energy transport due to the flow of particles, and a second one accounting for the particle-particle scattering which redistribute the energy in the interacting region. This expression has been compared to the energy variation into the corresponding lead and we have found that the two are not always equal and opposite in sign as it is in the case of particle current. The difference between these two currents is due to the direct exchange of energy among two of the leads and the microscopic mechanism behind it is the coupling of these leads through virtual processes involving the central interacting region.

Acknowledgements.

NLG akcnowledges P. Erdeman, and M. V. Moskalets and P. Burset Atienza for insightful comments on the expression for the energy current. The authors acknowledge financial support from the Academy of Finland Center of Excellence program (Project no. 312058) and the Academy of Finland (Project no. 287750). NLG acknowledges financial support from the Turku Collegium for Science and Medicine (TCSM). Numerical simulations were performed exploiting the Finnish CSC facilities under the Project no. 2000962 (”Thermoelectric effects in nanoscale devices”).

Appendix A ENERGY VARIATION OF THE CENTRAL REGION

We discuss here how to compute the expression of the energy variation in the central region:

[TABLE]

which can be recast in a form very similar to the time dependent Meir-Wingreen one for the particle current  Meir1992 ; Jauho1994 ; riku2013 in terms of the single-particle Green’s function, the embedding self-energies and, in this case, of the many-body self-energy. Here we show how the energy variation can be expressed in terms of single- and two-particle mixed Green’s functions. In the next subsection we demonstrate how to manipulate these expressions in order to recast them in terms of single-particle quantities of the central region. First, we need to compute explicitly the commutators entering into the definition of the energy current:

[TABLE]

where $z$ is the complex time variable on the path $\gamma$ in the complex-time plane and $1\equiv(\mathbb{x}_{1},z_{1})$, $2\equiv(\mathbb{x}_{2},z_{2})$ and $1^{\prime}\equiv(\mathbb{x}_{1}^{\prime},z_{1}^{\prime})$ are multi-indexes for position, spin and complex time. It is understood that at the end we take $z_{1}=t$ to project the complex time variable onto the real axis.

We can rewrite the first two terms in the following way

[TABLE]

where we use the definition of the mixed Green’s functions: $G_{C\alpha}(1,1^{\prime})=-i\left\langle\mathcal{T}_{\gamma}\hat{\psi}(1)\hat{\psi}_{\alpha}^{\dagger}(1^{\prime})\right\rangle$ and $G_{\alpha C}(1;1^{\prime})=-i\left\langle\mathcal{T}_{\gamma}\hat{\psi}_{\alpha}(1)\hat{\psi}^{\dagger}(1^{\prime})\right\rangle$, with $\mathcal{T}_{\gamma}$ the time-ordering operator over the Keldysh contour $\gamma$ in the complex-time plane and where we exploit the property $G_{\alpha C}(1;1^{\prime})=-[G_{C\alpha}(1^{\prime};1)]^{*}$. Moreover $1^{+}\equiv(\mathbb{x}_{1},z_{1}^{+})$ with $z_{1}^{+}$ being a time infinitesimally greater than $z_{1}$ on the Keldysh contour.

A similar manipulation can be done in the other four terms arising from the commutator of the interacting Hamiltonian. By a change of the integration variables and using the relation $\left\langle\hat{\psi}^{\dagger}(1)\hat{\psi}^{\dagger}(1^{\prime})\hat{\psi}(1^{\prime})\hat{\psi}_{\alpha}(1)\right\rangle=-\left\langle\hat{\psi}^{\dagger}(1^{\prime})\hat{\psi}^{\dagger}(1)\hat{\psi}(1^{\prime})\hat{\psi}_{\alpha}(1)\right\rangle$, the last four terms are equal in pair, this cancels out the factor $1/2$ and we are left with:

[TABLE]

where we use the definition of the mixed two-particle Green’s functions $G^{(2)}_{CCC\alpha}(1,2;3,4)=(-i)^{2}\left\langle\mathcal{T}_{\gamma}\hat{\psi}(1)\hat{\psi}(2)\hat{\psi}^{\dagger}_{\alpha}(4)\hat{\psi}^{\dagger}(3)\right\rangle$ and $G^{(2)}_{C\alpha CC}(1,2;3,4)=(-i)^{2}\left\langle\mathcal{T}_{\gamma}\hat{\psi}(1)\hat{\psi}_{\alpha}(2)\hat{\psi}^{\dagger}(4)\hat{\psi}^{\dagger}(3)\right\rangle$ together with the relation $G^{(2)}_{C\alpha CC}(1,1^{\prime};1^{\prime+},1^{+})=\left[G^{(2)}_{CCC\alpha}(1^{\prime},1;1^{+},1^{\prime+})\right]^{*}$.

By collecting all terms we obtain:

[TABLE]

Appendix B THE MIXED GREEN’S FUNCTIONS

We have seen that it is possible to express the energy current flowing through the central interacting region as the sum of two contributions. These contributions contain the mixed Green’s functions accounting for the propagation of both particle in the central region and in the leads. It is important to point of that such a situation occurs also in the case of the calculation of the particle (charge) current  riku2013 . The usual way to find an expression for the mixed Green’s functions $G^{C\alpha}(1;1^{\prime})$ is based on the equation-of-motion approach  Haugbook ; stefleebook . Since the leads are noninteracting, the method allows to write down a closed set of equations and then a general formula for $G^{C\alpha}(1;1^{\prime})$ in terms of the single particle Green’s function. Even though appealing for its simplicity, the equation-of-motion technique cannot straightforward apply to the mixed two-particle Green’s function $G_{2}^{CCC\alpha}(1,2;3,4)$ since the analysis gets quite complicated and a closed set of equations can be found only if one relies on some physical approximations. The latter must be chosen consistently with the approximations used for the single-particle Green’s function.

Rather than the equations-of-motion technique we find a general expression for both mixed Green’s functions by a direct expansion of the $S$-matrix in the interaction picture with respect to the coupling Hamiltonian  Haugbook . Despite the fact that the derivation is somehow more complex for the mixed single-particle Green’s function, it is very general and can be extended easily to the mixed n-particle Green’s functions.

Let’s look at the first term of the energy current, namely the one containing the the single-particle Hamiltonian $h(1)$. The key idea of this approach is to express the contour-ordered mixed single-particle Green’s function $G_{C\alpha}(1;1^{\prime})=-i\left\langle\mathcal{T}_{\gamma}\hat{\psi}(1)\hat{\psi}_{\alpha}^{\dagger}(1^{\prime})\right\rangle$ in terms of the contour-ordered single particle Green’s function of the particles in the central region $G_{CC}(1;1^{\prime})=-i\left\langle\mathcal{T}_{\gamma}\hat{\psi}(1)\hat{\psi}^{\dagger}(1^{\prime})\right\rangle$ and the one in the leads $g_{\alpha}(1;1^{\prime})=-i\left\langle\mathcal{T}_{\gamma}\hat{\psi}_{\alpha}(1)\hat{\psi}_{\alpha}^{\dagger}(1^{\prime})\right\rangle$.

The derivation follows by writing the Green’s function $G_{C\alpha}(1;1^{\prime})$ in terms of the interaction-picture operators (denoted by a tilde) with respect to the free Hamiltonians $\hat{H}_{C}$ and $\hat{H}_{\alpha}$ of both the central region and the terminal $\alpha$. Therefore the evolution operator will be expressed in terms of the coupling Hamiltonian $\tilde{H}_{T}$ in the interaction picture. Thus in this picture the mixed single-particle Green’s function can be written as:

[TABLE]

where we defined the S-matrix as:

[TABLE]

By inserting the explicit form of the S-matrix into Eq. 6 we obtain

[TABLE]

The key point in the second step is the assumption that the leads are noninteracting allowing to use the Wick’s theorem for the $\alpha$-operators. Besides, we use the fact that in the interaction picture the operators $\psi$ and $\psi_{\alpha}$ are independent and thus the expectations values can be factorized. Finally, by relabeling all integration variables in the remaining $k-1$ it turns out that all these terms are equal and hence we get a factor $k$, therefore:

[TABLE]

where in the last line we use $\hat{\psi}_{\alpha}(1)=\tilde{\psi}_{\alpha}(1)$, and restore the S-matrix expansion for the single particle Green’s function.

The first term in Eq.5 can be rewritten as:

[TABLE]

where we define the embedding self-energy $\Sigma_{\alpha}(1;2)=T_{\alpha}(1)g_{\alpha}(1;2)T_{\alpha}^{*}(2)$. It accounts for the presence of the $\alpha-$ lead. The second term of Eq.5 containing the mixed two-particle Green’s function $G_{CCC\alpha}^{(2)}$ can be manipulated in a similar way and rewrite it in terms of the contour-ordered two particle Green’s function of the central interacting region $G^{(2)}$:

[TABLE]

where in the last step we recognize the series expansion in the interaction picture of the two-particle Green’s function of the central region and the single-particle Green’s function of the lead $\alpha$.

Inserting this result into the second term of Eq. (5) we obtain:

[TABLE]

where we use $v(1,1^{\prime})=v(1^{\prime},1)$, the symmetry relations of the two-particle Green’s function $G^{(2)}(1,2;3,4)=-G^{(2)}(1,2;4,3)$ and the definition of the embedding self-energy $\Sigma_{\alpha}(1;2)$.

The crucial point in obtaining the above result is to use the relation linking the two-particle Green’s function in the interacting region with the many-body self-energy:

[TABLE]

which is nothing but the relation which defines the many-body self-energy itself.

Now we can write the final expression in terms of the currents, which in Eq. 1 gives:

[TABLE]

The expression for the variation of the energy has to be compared to the derivative of the total energy:

[TABLE]

Appendix C ENERGY VARIATION OF THE CONTACT REGION

Here we show how to compute the variation of energy in the contact region. Instead of considering directly the commutator, we exploit some useful identities that make the calculation faster and straightforward. The average variation of energy in the tunneling region then reads:

[TABLE]

where we use Eq. (8) together with the equation of motion for the single particle Green’s functions in the central region and in the leads:

[TABLE]

and where $\Sigma_{emb}(1;1^{\prime})=\sum\limits_{\beta}\Sigma_{\beta}(1;1^{\prime})$.

It is easy to recognize that the first term is nothing but the opposite of the variation of the energy of the lead $\alpha$: $\dot{\mathcal{E}}_{\alpha}=\frac{d}{dt}\langle\hat{H}_{\alpha}(t)\rangle$. The second term is $-\mathcal{J}_{C}^{(\alpha)}(t)$ ,namely the opposite of the variation of the energy on the central region due to the coupling with the lead $\alpha$.

The third term is the most interesting one and describes the direct coupling of the lead $\alpha$ with all the others. The coupling is mediated by the central region and therefore is of the fourth order in the coupling between the central region and the leads. To see it explicitly, we consider the mixed single particle Green’s function $G_{\beta\alpha}(1;1^{\prime})\equiv-i\left\langle\mathcal{T}_{\gamma}\hat{\psi}_{\beta}(1)\hat{\psi}_{\alpha}^{\dagger}(1^{\prime})\right\rangle$ and use once again the $S$-matrix expansion:

[TABLE]

Physically this term represents exactly the scattering of a particle (or a hole) from the lead $\alpha$ to the lead $\beta$ through the central region. It is now easy to see that the third term in Eq. (C) can be written as:

[TABLE]

It is simple to check that the term $\beta=\alpha$ does not contribute since it is purely imaginary.

Appendix D SELF CONSISTENT GW APPROXIMATION

In this section we recall the main equations to treat the interaction within the GW approximation. It is useful to split the full interaction self-energy into its Hartree and exchange-correlation parts

[TABLE]

The Hartree term, or first order polarization diagram, it is local in time and can be written as $\Sigma(1;2)_{H}=-i\delta(1;2)\int d3v(1;3)G(3;3^{+})$.

In the GW approximation the exchange-correlation part of the self-energy is given as a product of the Green function $G$ with a dynamically screened interaction $W$:

[TABLE]

The dressed interaction $W$ satisfies the self-consistent equation

[TABLE]

with $v$ the bare two-body interaction, and $P(1;2)=-iG(1;2)G(2;1)$ the irreducible polarization diagram. As it clear from the previous expressions the self-energy happen to be dependent on the Green’s function and thus it must be determined self-consistently in conjunction with the Dyson equation talarico2019 . Moreover, the GW self-energy can be express as the functional derivative of a so-called $\Phi$ functional, i.e., $\Sigma_{GW}[G]=\delta\Phi[G]/\delta G$ . This has been proven to be an effective way to guarantee macroscopic conservation laws stefleebook .

The reliability of the fully self-consistent GW approximation can be easily seen from Fig. 4 where we compare the expression obtained for the total current through the interaction region $\mathcal{J}(t)$ in Eq. (1) alongside the derivative of the total energy in the same region $d\mathcal{E}(t)/dt$, with \mathcal{E}(t)=\langle\hat{H}_{C}(t)\rangle=\left[i\int d\mathbb{x}_{1}d2\;\biggl{[}h(1)\delta(1,2)+\frac{1}{2}\Sigma(1,2)\biggr{]}G(2,1^{+})\right]. It is possible to appreciate the agreement within the numerical accuracy of the two quantities in the strong coupling and for large charging energy at different gate-voltage. The agreement holds also for the other set of parameters considered (not shown) and confirms the validity of the approximation.