Statistics of heat transport across capacitively coupled double quantum dot circuit

TL;DR

This paper investigates heat transport and fluctuations in a double quantum dot system, revealing how interactions and coupling influence heat current behavior and confirming fluctuation theorems through advanced theoretical methods.

Contribution

It introduces a comprehensive analysis of heat fluctuations using Lindblad and Schwinger-Keldysh approaches, highlighting the effects of Coulomb interaction and system-reservoir coupling.

Findings

Heat current and fluctuations are non-monotonic functions of interaction and coupling strength.

Large Coulomb interaction and coupling suppress heat transport and fluctuations.

Steady-state fluctuation theorem is verified for the quantum dot system.

Abstract

We study heat current and the full statistics of heat fluctuations in a capacitively-coupled double quantum dot system. This work is motivated by recent theoretical studies and experimental works on heat currents in quantum dot circuits. As expected intuitively, within the (static) mean-field approximation, the system at steady-state decouples into two single-dot equilibrium systems with renormalized dot energies, leading to zero average heat flux and fluctuations. This reveals that dynamic correlations induced between electrons on the dots is solely responsible for the heat transport between the two reservoirs. To study heat current fluctuations, we compute steady-state cumulant generating function for heat exchanged between reservoirs using two approaches : Lindblad quantum master equation approach, which is valid for arbitrary coulomb interaction strength but weak system-reservoir coupling strength, and the saddle point approximation for Schwinger-Keldysh coherent state path integral, which is valid for arbitrary system-reservoir coupling strength but weak coulomb interaction strength. Using thus obtained generating functions, we verify steady-state fluctuation theorem for stochastic heat flux and study the average heat current and its fluctuations. We find that the heat current and its fluctuations change non-monotonically with the coulomb interaction strength ($U$) and system-reservoir coupling strength ($\Gamma$) and are suppressed for large values of $U$ and $\Gamma$.

Equations

$$\hat{H}$$

$$P[\Delta Q;T-T_0]$$

$$\mathcal{Z}[\chi ;T-T_0]$$

$$H_{\chi }=\sum _{\alpha =L,R}{ϵ}_{\alpha }{c}_{\alpha }^{†}{c}_{\alpha }+U{c}_L^{†}{c}_L{c}_R^{†}{c}_R+{\sum _k}_{\alpha =L,R}{ϵ}_{\alpha ,k}{d}_{\alpha k}^{†}{d}_{\alpha k}$$

$$\displaystyle+\mathop{\sum_{k}}_{\alpha=L,R}\big{[}g_{\alpha k}e^{-i(\epsilon_{\alpha k}-\mu_{\alpha})\chi_{\alpha}}d_{\alpha k}^{\dagger}c_{\alpha}+g_{\alpha k}^{*}e^{i(\epsilon_{\alpha k}-\mu_{\alpha})\chi_{\alpha}}c_{\alpha}^{\dagger}d_{\alpha k}\big{]},$$

$$\rho (T_0)$$

$$\mathcal{Z}[\chi ;T-T_0]$$

$$\frac{d}{dT}{\rho }_s^{\chi }(T)$$

$$\frac{d}{dT}|P_s^{\chi }(T)⟩$$

$$|P_s^{\chi }(T)⟩$$

$${\mathcal{L}}^{\chi }$$

$$\displaystyle\mathcal{Z}[\chi;T-T_{0}]=e^{-\frac{\left(\Gamma_{L}+\Gamma_{R}\right)}{2}(T-T_{0})}\Bigg{\{}\cosh\left[\Lambda[\chi]\left(T-T_{0}\right)\right]+\frac{\sinh\left[\Lambda[\chi]\left(T-T_{0}\right)\right]}{\Lambda[\chi]}\Bigg{[}\underset{\alpha=L,R}{\sum}\Gamma_{\alpha}\Big{(}\frac{1}{2}$$

$$\displaystyle+\Big{[}f_{\alpha}\left[f_{SL}\left(1-f_{SR}\right)+f_{SR}\left(1-f_{SL}\right)\right]\left(e^{iU\chi_{\alpha}/2}-1\right)+\left(1-f_{\alpha}\right)\left[1-f_{SL}\left(1-f_{SR}\right)-f_{SR}\left(1-f_{SL}\right)\right]\left(e^{-iU\chi_{\alpha}/2}-1\right)\Big{]}\Big{)}\Bigg{]}\Bigg{\}},$$

$$\Lambda [\chi ]=\sqrt{{\left(\frac{{\Gamma }_L+{\Gamma }_R}{2}\right)}^2+{\Gamma }_L{\Gamma }_R\left[{\bar{f}}_L[1-{\bar{f}}_R]\left(e^{iU\chi /2}-1\right)+{\bar{f}}_R[1-{\bar{f}}_L]\left(e^{-iU\chi /2}-1\right)\right]}.$$

$$|P_s^0(T_0)⟩$$

$$\mathcal{F}[\chi ]$$

$$C_1$$

$$C_2$$

$$C_3$$

$$C_4$$

$$\mathcal{Z}[\chi ;T-T_0]$$

$$\mathcal{Z}[\chi ;T-T_0]=\frac{1}{\mathcal{N}}\int \mathcal{D}[\{{\psi }_{\alpha }^{†}(\tau ),{\psi }_{\alpha }(\tau )\}]$$

$$\int \mathcal{D}[\{{\psi }_{\alpha k}^{†}(\tau ),{\psi }_{\alpha k}(\tau )\}]e^{\frac{i}{\hslash }{S}^{\chi }[\{{\psi }_{\alpha }^{†}(\tau ),{\psi }_{\alpha }(\tau )\},\{{\psi }_{\alpha k}^{†}(\tau ),{\psi }_{\alpha k}(\tau )\}]},$$

$$S^{\chi }[\{{\psi }_{\alpha }^{†}(\tau ),{\psi }_{\alpha }(\tau )\},\{{\psi }_{\alpha k}^{†}(\tau ),{\psi }_{\alpha k}(\tau )\}]$$

$$\displaystyle=\sum_{\alpha,\alpha^{\prime}=L,R}\int_{c}d\tau\int_{c}d\tau^{\prime}\Big{[}\psi_{\alpha}^{\dagger}(\tau)[G_{sys}^{0}]^{-1}_{\alpha\alpha^{\prime}}(\tau,\tau^{\prime})\psi_{\alpha^{\prime}}(\tau^{\prime})+\sum_{k,k^{\prime}}\psi_{\alpha k}^{\dagger}(\tau)[G_{res}^{0}]^{-1}_{\alpha k\alpha^{\prime}k^{\prime}}(\tau,\tau^{\prime})\psi_{\alpha^{\prime}k^{\prime}}(\tau^{\prime})\Big{]}$$

$$\displaystyle-\sum_{\alpha=L,R}\sum_{k}\int_{c}d\tau\Big{[}g_{\alpha k}e^{-i(\epsilon_{\alpha k}-\mu_{\alpha})\chi_{\alpha}(\tau)}\psi_{\alpha k}^{\dagger}(\tau)\psi_{\alpha}(\tau)+g_{\alpha k}^{*}e^{i(\epsilon_{\alpha k}-\mu_{\alpha})\chi_{\alpha}(\tau)}\psi_{\alpha}^{\dagger}(\tau)\psi_{\alpha k}(\tau)\big{]}-\int_{c}d\tau U\psi_{L}^{\dagger}(\tau)\psi_{L}(\tau)\psi_{R}^{\dagger}(\tau)\psi_{R}(\tau),$$

$$\displaystyle\sum_{\alpha_{1}=L,R}\int_{c}d\tau_{1}\Big{[}\big{(}i\hbar\frac{\partial}{\partial\tau}-\epsilon_{\alpha}\big{)}\delta_{\alpha\alpha_{1}}\delta^{c}(\tau,\tau_{1})\Big{]}[G_{sys}^{0}]_{\alpha_{1}\alpha^{\prime}}(\tau_{1},\tau^{\prime})=\delta_{\alpha\alpha^{\prime}}\delta^{c}(\tau,\tau^{\prime})$$

$$\displaystyle\sum_{\alpha_{1}=L,R}\int_{c}d\tau_{1}\Big{[}\big{(}-i\hbar\frac{\partial}{\partial\tau^{\prime}}-\epsilon_{\alpha^{\prime}}\big{)}\delta_{\alpha_{1}\alpha^{\prime}}\delta^{c}(\tau_{1},\tau^{\prime})\Big{]}[G_{sys}^{0}]_{\alpha\alpha_{1}}(\tau,\tau_{1})=\delta_{\alpha\alpha^{\prime}}\delta^{c}(\tau,\tau^{\prime})$$

$$\displaystyle\sum_{\alpha_{1}=L,R}\sum_{k_{1}}\int_{c}d\tau_{1}\Big{[}\big{(}i\hbar\frac{\partial}{\partial\tau}-\epsilon_{\alpha k}\big{)}\delta_{\alpha\alpha_{1}}\delta_{kk_{1}}\delta^{c}(\tau,\tau_{1})\Big{]}[G_{res}^{0}]_{\alpha_{1}k_{1}\alpha^{\prime}k^{\prime}}(\tau_{1},\tau^{\prime})=\delta_{\alpha\alpha^{\prime}}\delta_{kk^{\prime}}\delta^{c}(\tau,\tau^{\prime})$$

$$\displaystyle\sum_{\alpha_{1}=L,R}\sum_{k_{1}}\int_{c}d\tau_{1}\Big{[}\big{(}-i\hbar\frac{\partial}{\partial\tau^{\prime}}-\epsilon_{\alpha^{\prime}k^{\prime}}\big{)}\delta_{\alpha_{1}\alpha^{\prime}}\delta_{k_{1}k^{\prime}}\delta^{c}(\tau_{1},\tau^{\prime})\Big{]}[G_{res}^{0}]_{\alpha k\alpha_{1}k_{1}}(\tau,\tau_{1})=\delta_{\alpha\alpha^{\prime}}\delta_{kk^{\prime}}\delta^{c}(\tau,\tau^{\prime})$$

$$[G_{sys}^0{]}_{\alpha {\alpha }^{\prime }}(T_0^{-},{\tau }^{\prime })=-e^{{\beta }_S\left({ϵ}_{\alpha }-{\mu }_S\right)}[G_{sys}^0{]}_{\alpha {\alpha }^{\prime }}(T_0^{+},{\tau }^{\prime })$$

$$[G_{sys}^0{]}_{\alpha {\alpha }^{\prime }}(\tau ,T_0^{-})=-e^{-{\beta }_S\left({ϵ}_{{\alpha }^{\prime }}-{\mu }_S\right)}[G_{sys}^0{]}_{\alpha {\alpha }^{\prime }}(\tau ,T_0^{+})$$

$$[G_{res}^0{]}_{\alpha k{\alpha }^{\prime }{k}^{\prime }}(T_0^{-},{\tau }^{\prime })=-e^{{\beta }_{\alpha }\left({ϵ}_{\alpha k}-{\mu }_{\alpha }\right)}[G_{res}^0{]}_{\alpha k{\alpha }^{\prime }{k}^{\prime }}(T_0^{+},{\tau }^{\prime })$$

$$[G_{res}^0{]}_{\alpha k{\alpha }^{\prime }{k}^{\prime }}(\tau ,T_0^{-})=-e^{-{\beta }_{{\alpha }^{\prime }}\left({ϵ}_{{\alpha }^{\prime }{k}^{\prime }}-{\mu }_{{\alpha }^{\prime }}\right)}[G_{res}^0{]}_{\alpha k{\alpha }^{\prime }{k}^{\prime }}(\tau ,T_0^{+})$$

$$[G_{sys}^0{]}_{\alpha {\alpha }^{\prime }}(\tau ,{\tau }^{\prime })=-\frac{i}{\hslash }{e}^{-\frac{i}{\hslash }{ϵ}_{\alpha }\left(\tau -{\tau }^{\prime }\right)}$$

$${\delta }_{\alpha {\alpha }^{\prime }}\left[\Theta (\tau ,{\tau }^{\prime })\left[1-f_S({ϵ}_{\alpha })\right]-\Theta ({\tau }^{\prime },\tau )f_S({ϵ}_{\alpha })\right]$$

$$[G_{res}^0{]}_{\alpha k{\alpha }^{\prime }{k}^{\prime }}(\tau ,{\tau }^{\prime })=-\frac{i}{\hslash }{e}^{-\frac{i}{\hslash }{ϵ}_{\alpha k}\left(\tau -{\tau }^{\prime }\right)}$$

Full text

I Introduction

Studying transport processes in nano sized electronic quantum-dot junctions has been an active research area for last two decades Van2002 ; Joachim2005 ; Cuniberti2006 ; Tao2006 ; Cuevas2010 ; Baldea2016 ; Xiang2016 . The motivation being two fold : desire to design more efficient electronic devices and heat engines Sothmann2014 , and also as a platform for testing fundamental principles. From the technological perspective, useful devices have been proposed theoretically and few are tested experimentally. For example, nano diodes Aviram1974 , transistors Cardamone2006 , switches and other electronic elements relevant for device applications have been proposed. Understanding of charge and heat transport in nano systems are relevant for these applications. However, due to the small size, fluctuations of fluxes flowing through these systems are not negligible. These fluctuations are not arbitrary but follow universal relations called fluctuation theorems which generalize second law of thermodynamics to small scale Lebowitz1999 ; Jarzynski1997 ; Crooks1998 ; Crooks1999 ; Kurchan2000 ; Tasaki2000 ; Jarzynski2004 ; Monnai2005 ; Esposito2009 ; Campisi2011 ; Seifert2012 . These beautiful identities relate number of microscopic realizations of transport processes which produce certain amount of entropy to those which annihilate the same amount of entropy. Nano-electronic devices have served as useful platforms for testing these identities Utsumi2010 ; Kung2012 ; Pekola2015 ; Pekola2018 . These theorems are not only aesthetically appealing, but also are used to gain insights into transport processes. For example they have been used to characterize efficiency fluctuations Verley2014unlikely ; Verley2014universal ; Esposito2015efficiency of nano heat engines, which is an important fundamental generalization of Carnot’s analysis Callen2013 of macroscopic heat engines to micro scale. Thermoelectric engines, which are of current theoretical and experimental interest, constitutes one such class of nano heat engines Dubi2011 ; Sothmann2014 ; Perroni2016 ; Benenti2017 ; Cui2017 that convert heat to electrical work.

Although heat flow plays a central role in determining the efficiencies of these engines, heat currents at nano electronic junctions are not as well explored as the charge currents. Recently there has been some interest in exploring the effects of various many-body interactions in thermoelectric heat engines. For example, effects of electron-phonon and electron-electron interactions on efficiencies of two-terminal and three-terminal thermoelectric engines are studied Azema2012 ; Zimbovskaya2014 ; Agarwalla2015 ; Perroni2015 ; Sierra2016 ; Thierschmann2016 ; Erdman2017 ; Dare2017 ; Friedman2017 ; Walldorf2017 ; Svilans2018 ; Daroca2018 . Furthermore, important experimental advancements in measuring heat currents in nano-electric junctions have been achieved recently Jezouin2013 ; Dutta2017 ; Cui2017 .

Motivated by these works, we study steady-state heat flux and fluctuations across capacitively (coulomb) coupled double quantum dot system. This system is known to act as a heat rectifier in some parameter regime Ruokola2011 . Capacitive coupling has been used to probe charge fluctuations in nano-junctions Cuetara2011 ; Golubev2011 ; Cuetara2015 , to understand coulomb drag effects, where electric charge flux in a circuit induces a charge flux in another capacitively coupled circuit Moldoveanu2009 ; Aita2013 ; Kaasbjerg2016 ; Keller2016 ; Narozhny2016 ; Zhou2019 . In this work we are interested in the study of heat-flux that is induced due to coulomb interactions between electrons in a capacitively coupled double-quantum dot system. To study heat flux and its fluctuations, we calculate cumulant generating function, defined using two point measurement scheme, using two different approaches valid in different parameter regimes (far above the Kondo-temperature Hewson1997 ). Lindblad quantum master equation Breuer2002 ; Bagrets2003 ; Harbola2006 ; Harbola2007 valid at high temperatures, weak system reservoir coupling strength and arbitrary coulomb interaction strength, and Schwinger-Keldysh Rammer2007 ; Altland2010 ; Kamenev2011 ; Stefanucci2013 saddle-point method (random-phase approximation with mean-field dressed propagators) Hamann1969 ; Hamann1970 ; Morandi1974 ; Altland2010 ; Kamenev2011 valid for weak coulomb interaction strength and arbitrary system reservoir coupling strength. We verify steady-state heat fluctuation theorem and calculate heat flux and its fluctuations. Heat flux flowing through the same model system Wang2018 and fluctuations of heat flow in a variant of this model have been studied recently to understand near-field radiative heat transfer within bare random-phase approximation Tang2018 . Here we present results that are valid beyond bare random-phase approximation (random phase approximation with mean-field dressed propagators). We find that the steady-state scaled cumulant generating functions obtained using both the approximation schemes satisfy Gallavotti-Cohen symmetry and hence the steady-state fluctuation theorem for the heat fluctuations. Heat flux and its fluctuations are non-monotonic functions of coulomb interaction strength and decay exponentially for asymptotically large coulomb interaction strength. Similar non-monotonic behavior is seen with respect to system-reservoir coupling strength. The flux and its fluctuations are suppressed as a power law ($\Gamma^{-4}$) for large coupling strength.

In section II we introduce the model system. In section III we define moment generating function for heat fluctuations using two-point measurement scheme and calculate it using two approximation schemes and discuss heat flux and fluctuations in two subsections. We conclude in section IV.

II Model system

Schematic of the model system considered in this work is shown in fig. (1). It consists of two capacitively (Coulomb) coupled quantum dots (each having single orbital) individually coupled to two different fermionic reservoirs. The whole system is described by the following Hamiltonian,

[TABLE]

Here $c_{\alpha}^{{\dagger}}$ ($c_{\alpha}$) and $d_{\alpha k}^{\dagger}$ ($d_{\alpha k}$) stand for fermionic creation (annihilation) operator for creating (annihilating) an electron in the $\alpha^{th}$ ($\alpha=L,R$) quantum dot and in the state labeled by ’$k$’ in the $\alpha^{th}$ fermionic reservoir respectively. The first term in Eq. (1) represents Hamiltonian of two isolated quantum dots each having a single orbital with energies $\epsilon_{\alpha}$, the second term represents coulomb interaction between electrons on the two quantum dots, the third term is the Hamiltonian for the free electrons in the reservoirs, and the last term stands for hybridization between electrons on quantum dots and the reservoirs. Throughout this work we assume wide-band approximation i.e., we assume that $g_{\alpha k}$ is independent of $k$ and density of states of reservoirs are constant functions of energy. We note that the Hamiltonian given in Eq. (1) is a variant of the famous Anderson Hamiltonian Anderson1961 ; Anderson1978 ; Hewson1997 .

III Moment generating function

When the quantum dots are brought together and are coupled to two reservoirs, energy and particles are exchanged. In this work, we are interested in calculating statistics of steady-state fluxes flowing through the double quantum dot system. In the long-time limit, only the heat flows between the left and the right reservoirs. The heat fluxes at the left and the right interface are balanced at the steady-state. The particle flux between system and reservoirs vanishes at steady-state. Physical reason for this is : (i) Since the coupling between the two quantum dots does not change number of particles on the dots (i.e., the particle exchange between the two dots is not allowed by microscopic dynamics), the net number of particles exchanged between the left(right) dot and the left(right) reservoir is constrained to [math] and $1$, hence the particle flux and fluctuations are suppressed at the steady-state. (ii) Similarly energy cannot indefinitely accumulate on the system due to the boundedness of the systems energy spectrum, energy flux at the left and the right interfaces balance out at long-times.

Distribution function, $P[\Delta Q;T-T_{0}]$, for heat ($\Delta Q$) flowing from the right reservoir to the left reservoir within a time $T-T_{0}$ can be obtained using two-point measurement protocol Esposito2009 ; Campisi2011 for the observable corresponding to the operator $\frac{1}{2}[(H_{L}-\mu_{L}N_{L})-(H_{R}-\mu_{R}N_{R})]$ as,

[TABLE]

where $\mathcal{Z}[\chi;T-T_{0}]$ is the moment generating function for $\Delta Q$ and is given as,

[TABLE]

where $\mathcal{U}_{\chi}(T_{1},T_{2})=e^{-\frac{i}{\hbar}(T_{1}-T_{2})H_{\chi}}$ with

[TABLE]

and $\chi_{L}=-\chi_{R}=\frac{\chi}{2}$. The trace in Eq. (III) is over the combined Fock space of the system and reservoirs and $\rho(T_{0})$ is the density matrix of the whole system at initial time $T_{0}$ taken here as

[TABLE]

i.e., system and reservoirs initial states are non-interacting equilibrium states with different temperatures and chemical potentials.

Below we calculate $\mathcal{Z}[\chi;T-T_{0}]$ approximately using two approaches, (i) Lindblad quantum master equation approach where coupling between system and reservoirs is assumed to be weak and, (ii) Schwinger-Keldysh path-integral approach where coulomb interaction strength is assumed to be weak but system reservoir coupling can be arbitrary.

III.1 Lindblad quantum master equation approach

$\mathcal{Z}[\chi;T-T_{0}]$ defined in Eq. (III) can be expressed as,

[TABLE]

where \rho_{s}^{\chi}(T)=\mbox{\boldmathTr}_{\text{Res}}\big{[}\mathcal{U}_{\chi}(T,T_{0})\rho(T_{0})\mathcal{U}_{\{0,0\}}(T_{0},T)\big{]} is the counting field dependent reduced system density matrix at time $T$ obtained by tracing out the two reservoirs. Using the standard Born-Markov-Secular approximations (neglecting Lamb shifts), (counting field dependent) Lindblad quantum master equation can be derived for $\rho_{s}^{\chi}(T)$ Breuer2002 ; Bagrets2003 ; Harbola2006 ; Harbola2007 , which is given as,

[TABLE]

here $\alpha\neq\bar{\alpha}=L,R$; $\Gamma_{\alpha}=\frac{2\pi}{\hbar}|g_{\alpha}|^{2}\rho_{\alpha}$ and $f_{\alpha}(\hat{o})=\left(e^{\beta_{\alpha}(\hat{o}-\mu_{\alpha})}+1\right)^{-1}$ with $\beta_{\alpha}$ and $\mu_{\alpha}$ being the temperature and chemical potential of the $\alpha^{\text{th}}$ reservoir. It is important to note that if (static) mean-field approximation is made here, i.e., replacing $c_{\bar{\alpha}}^{\dagger}c_{\bar{\alpha}}$ by \langle c_{\bar{\alpha}}^{\dagger}c_{\bar{\alpha}}\rangle_{S}=\lim_{\left(T-T_{0}\right)\to\infty}\frac{\mbox{\boldmathTr}_{s}\big{[}c_{\bar{\alpha}}^{\dagger}c_{\bar{\alpha}}\rho_{s}^{\chi}(T)\big{]}}{\mbox{\boldmathTr}_{s}\big{[}\rho_{s}^{\chi}(T)\big{]}}, the right hand side of Eq. (III.1) can be separated into two terms which depend only on the dynamics of individual dots whose energies are renormalized by coupling to the other dot. This results in two decoupled quantum dots which equilibrate with their own reservoirs at long-time. Thus within this approximation, heat flux and fluctuations through the system vanish at steady-state. Hence mean-field approximation leads to no steady-state heat flux and fluctuations. And one needs to go beyond the mean-field approximation for having non-zero flux and fluctuations at steady state.

By taking matrix elements in the occupation number basis of the two dots, $|N_{L},N_{R}\rangle$ (with $N_{\alpha}=0,1$), it can be seen that the populations ($\langle N_{L},N_{R}|\rho_{s}^{\chi}(T)|N_{L},N_{R}\rangle$) are decoupled from the coherences ($\langle N_{L},N_{R}|\rho_{s}^{\chi}(T)|N_{L}^{\prime},N_{R}^{\prime}\rangle$), which die out exponentially fast with time. Further, we restrict ourselves to a parameter regime : $\epsilon_{\alpha}=\left(\mu_{\alpha}-\frac{U}{2}\right)$, which simplifies the analysis. In this regime, we only need to solve the following $2\times 2$ matrix equation,

[TABLE]

where

[TABLE]

and the Liouvillian $\mathcal{L}^{\chi}$ is given as,

[TABLE]

here $\bar{f}_{\alpha}=\left(e^{\beta_{\alpha}U/2}+1\right)^{-1}$. Note that the structure of the Liouvillian given in Eq. (8) is very similar to the case of charge transport through a resonant level system Goswami2015 when the two many-body states of the level are identified with the singly occupied and doubly the (un)occupied states of the double quantum dot system. Using solution of Eq. (8) in Eq. (6) (equivalent to $\mathcal{Z}[\chi;T-T_{0}]=\langle\mathcal{I}|e^{\mathcal{L}^{\chi}(T-T_{0})}|P_{s}^{0}(T_{0})\rangle$, with $\langle\mathcal{I}|=\begin{bmatrix}1&1\end{bmatrix}$), $\mathcal{Z}[\chi;T-T_{0}]$ is obtained as,

[TABLE]

where $f_{S\alpha}=\left(e^{\beta_{S}\left(\epsilon_{\alpha}-\mu_{S}\right)}+1\right)^{-1}$ with $\beta_{S}$ and $\mu_{S}$ are temperature and chemical potential of the uncoupled quantum dots, and

[TABLE]

To arrive at explicit expression for $\mathcal{Z}[\chi;T-T_{0}]$, we have used the initial condition i.e.,

[TABLE]

which is equivalent to $\rho_{s}^{\chi}(T_{0})=\frac{e^{-\beta_{S}\left[H_{S}-\mu_{S}N_{S}\right]}}{Tr[e^{-\beta_{S}\left[H_{S}-\mu_{S}N_{S}\right]}]}$.

In the long-time limit (i.e., $T-T_{0}\to\infty$), the scaled cumulant generating function defined as,

[TABLE]

is given by

[TABLE]

This scaled cumulant generating function has the same form as that of for charge transport through a resonant level model Goswami2015 . This is due to the mapping between the two models as discussed earlier. It is straight forward to see that the cumulant generating function, $\mathcal{F}[\chi]$, satisfies Gallavotti-Cohen symmetry : $\mathcal{F}[-\chi-i(\beta_{L}-\beta_{R})]=\mathcal{F}[\chi]$. This symmetry leads to the detailed steady-state fluctuation theorem for the distribution function for heat flow : $\lim_{(T-T_{0})\to\infty}\frac{P[+\Delta Q;T-T_{0}]}{P[-\Delta Q;T-T_{0}]}=e^{(\beta_{L}-\beta_{R})\Delta Q}$.

Further, using the above long-time limit scaled cumulant generating function, $\mathcal{F}[\chi]$, cumulants of heat flux can be obtained as $C_{n}=i^{n}\frac{d}{d\chi^{n}}\mathcal{F}[\chi]$. Analytical expression for first four scaled cumulants, i.e, heat flux ($C_{1}$), heat noise ($C_{2}$), third cumulant ($C_{3}$) and fourth cumulant ($C_{4}$) are given as,

[TABLE]

Figure (2) shows the four cumulants as a function of coulomb interaction strength. It is clear from this figure that heat flux and fluctuations are suppressed exponentially for large $U$. This is due to the exponential dependence of Fermi functions on $U$. Physically, the transition between the singly occupied states to doubly (un)occupied state become less probable as $U$ is increased Ruokola2011 . Further, we note that for intermediate values of $U$, fluctuations of heat are enhanced.

Since, $\mathcal{F}[\chi]$ is a periodic function of $\chi$ with period $\frac{4\pi}{U}$, this implies that $\underset{\left(T-T_{0}\right)\to\infty}{\lim}P[\Delta Q;T-T_{0}]$ has the Dirac comb structure : $\underset{\left(T-T_{0}\right)\to\infty}{\lim}P[\Delta Q;T-T_{0}]=\sum_{n=-\infty}^{+\infty}p[n;T-T_{0}]\delta[\Delta Q-\frac{nU}{2}]$ with $p[n;T-T_{0}]=\int_{0}^{2\pi}\frac{d\chi}{2\pi}e^{\mathcal{F}[\frac{2\chi}{U}]\left(T-T_{0}\right)}e^{i\chi n}$.

$p[n;T-T_{0}]$ is computed numerically and is shown in fig. (3) along with $\ln\frac{p[n,T-T_{0}]}{p[-n,T-T_{0}]}$ in the inset, demonstrating the validity of steady-state Gallavotti-Cohen fluctuation theorem for the stochastic heat flow.

In the next sub-section we present results obtained within saddle-point approximation for path-integral formulation of Schwinger-Keldysh technique.

III.2 Schwinger-Keldysh path integral approach

We compute the moment generating function ($\mathcal{Z}[\chi;T-T_{0}]$ ) using path-integral on Schwinger-Keldysh contour. The results obtained are valid for arbitrary dot-reservoir coupling strength. However the effect of the coulomb interaction is incorporated approximately.

$\mathcal{Z}[\chi;T-T_{0}]$, defined in Eq. (III), can be expressed as

[TABLE]

where $\mathcal{T}_{c}e^{-\frac{i}{\hbar}\int_{c}d\tau H_{\chi(\tau)}(\tau)}$ is the evolution operator defined on the Schwinger-Keldysh contour Rammer2007 ; Haug2008 ; Kita2010 ; Kamenev2011 ; Stefanucci2013 shown in Fig. (4), going from $T_{0}$ to $T$ and back to $T_{0}$. Here $H_{\chi(\tau)}(\tau)=H_{\chi}$ on forward contour and $H_{\chi(\tau)}(\tau)=H_{\chi=0}$ on the backward contour.

$\mathcal{Z}[\chi;T-T_{0}]$ can be expressed as a functional integral using Grassman field variables Negele1988 ; Rammer2007 ; Altland2010 ; Kamenev2011 , $\{\psi_{\alpha}^{\dagger}(\tau),\psi_{\alpha}(\tau)\}$ for the system and $\{\psi_{\alpha k}^{\dagger}(\tau),\psi_{\alpha k}(\tau)\}$ for the reservoirs. This gives,

[TABLE]

where $\mathcal{N}$ is the normalization constant (independent of $\chi$) such that $\mathcal{Z}[\chi;T-T_{0}]|_{\chi=0}=1$. Here, we do not compute $\mathcal{N}$ explicitly, and modify it at intermediate steps by absorbing all constants (’$\chi$’ independent). Its value is determined finally by imposing $\mathcal{Z}[\chi;T-T_{0}]|_{\chi=0}=1$. $S^{\chi}[\{\psi_{\alpha}^{\dagger}(\tau),\psi_{\alpha}(\tau)\},\{\psi_{\alpha k}^{\dagger}(\tau),\psi_{\alpha k}(\tau)\}]$ is the action of the whole system, given as,

[TABLE]

where $\chi_{L}(\tau)=-\chi_{R}(\tau)=\frac{\chi}{2}$ on the forward contour and $\chi_{L}(\tau)=\chi_{R}(\tau)=0$ on the backward contour. Further, $[G_{sys}^{0}]^{-1}_{\alpha\alpha^{\prime}}(\tau,\tau^{\prime})$ and $[G_{res}^{0}]^{-1}_{\alpha k\alpha^{\prime}k^{\prime}}(\tau,\tau^{\prime})$ are matrix elements (with indices spanning state labels and contour times) of inverse of matrices with elements satisfying the following Schwinger-Dyson or Kadanoff-Baym equations on the Schwinger-Keldysh contour,

[TABLE]

with the following Kubo-Martin-Schwinger boundary conditions Martin1959 ; vanLeeuwen2013 ; Stefanucci2013 enforcing the information of the initial state of the system and reservoirs,

[TABLE]

Equation (III.2) is one of the ways to take care of the initial state information in the Schwinger-Keldysh path integral formalism Utsumi2002 ; Altland2010 ; Kamenev2011 ; Secchi2018 . The solution of Eqs. (III.2) along with the boundary conditions Eqs. (III.2) are,

[TABLE]

where $f_{X}(x)=\left(e^{\beta_{X}(x-\mu_{X})}+1\right)^{-1}$ ($X=L,R,S$). We integrate over the reservoir Grassman fields in Eq. (III.2) to get $\mathcal{Z}[\chi;T-T_{0}]$ as a path integral only over system Grassman fields as,

[TABLE]

with

[TABLE]

where

[TABLE]

with the self-energy acquired by the system due to the coupling with the reservoirs given as,

[TABLE]

The path integral given in Eq. (III.2) over system Grassman fields cannot be evaluated exactly due to the presence of the quartic term (with coupling constant $U$). Hence, we proceed to evaluate it approximately. To that end, we decouple the above quartic term by introducing auxiliary real fields (aka Hubbard-Stratonovich decoupling) which can be interpreted as fluctuating external potentials Stratonovich1957 ; Hubbard1959 ; Kleinert1978 ; Negele1988 ; Altland2010 ; Kamenev2011 . This gives,

[TABLE]

where

[TABLE]

with $[G_{\phi}^{c}]^{-1}(\tau,\tau^{\prime})$ being the inverse of $[G_{\phi}^{c}](\tau,\tau^{\prime})$ which satisfy the following equation,

[TABLE]

Since the Grassman path integral in Eq. (III.2) is quadratic in terms of system fields $\{\psi_{\alpha}^{\dagger}(\tau),\psi_{\alpha}(\tau)\}$, it can be performed exactly (here we have used the identity \ln\det[\mathbb{A}]=\mbox{\boldmathTr}[\ln\mathbb{A}]) to get,

[TABLE]

with

[TABLE]

where $Tr$ stands for trace over contour time and orbital indices. The above algebraic gymnastics doesn’t solve the problem as the final path integral, Eq. (III.2), has an action, Eq. (III.2), which is highly non-linear, nevertheless it is a bosonic path integral, which can be approximately evaluated using saddle-point/stationary-phase method. Within saddle-point approximation, action, $S^{\chi}[\{\phi_{\alpha}(\tau)\}]$ is functional Taylor expanded around the path $\{\phi_{\alpha}^{0}(\tau)\}$ which makes action stationary, i.e.,

[TABLE]

Further, the action is approximated by retaining terms in the functional Taylor expansions up to quadratic order, making the action functional, a quadratic form. This quadratic functional integral can be analytically evaluated to get a functional Fredholm determinant multiplied by the exponential of the action evaluated at the stationary path. With this hand-wavy description of saddle-point/stationary-phase approximation, we move ahead.

The saddle point equations for the action given in Eq. (III.2) are obtained as,

[TABLE]

where an infinitesimal forward shift ($\tau^{+}$) of the second argument compared to the first argument of $[G_{\phi_{x}^{0}}^{c}]_{xx}(\tau,\tau^{+})$ along the Schwinger-Keldysh contour can be deduced by consistently decoupling the fermionic quartic term using Hubbard-Stratonovich fields in discretized path integral. Otherwise there will be an ambiguity, as $[G_{\phi_{x}^{0}}^{c}]_{xx}(\tau,\tau^{\prime})$ is discontinuous at $\tau=\tau^{\prime}$ with jump discontinuity of magnitude $\frac{i}{\hbar}$. Equations (III.2) for $\{\phi_{\alpha}^{0}(\tau)\}$ together with Eqs. (29) for $[G_{\phi^{0}}^{c}]_{\alpha\alpha^{\prime}}(\tau,\tau^{\prime})$, constitute self-consistent system of equations, which may possess more than one solution. When more than one stationary solution exist, then the functional integral is approximated by summing over the result obtained by Gaussian approximating the action around each of the stationary solutions. Expanding the action given in Eq. (III.2) around the stationary path and retaining only quadratic term, we get approximate expression for $\mathcal{Z}[\chi;T-T_{0}]$ as (assuming that there is a unique stationary path),

[TABLE]

with $S^{\chi}_{app}[\{\phi_{\alpha}(\tau)\}]$ representing approximate action given as,

[TABLE]

where $P_{\alpha\alpha}^{0}(\tau,\tau^{\prime})=[G_{\phi^{0}}^{c}]_{\alpha\alpha}(\tau,\tau^{\prime})[G_{\phi^{0}}^{c}]_{\alpha\alpha}(\tau^{\prime},\tau)$ for $\alpha=L,R$ are the contour-ordered polarization propagators within random-phase approximation expressed in terms of mean-field system fermion propagators (solutions of Eqs. (29) and Eqs. (III.2)). Note that the polarization dependent terms in Eq. (III.2) represent the leading order correction to the mean-field (saddle-point) contribution. After a change of variables (shift transformation $\{\phi_{\alpha}(\tau)\}\to\{\phi_{\alpha}(\tau)+\phi_{\alpha}^{0}(\tau)\}$), followed by performing the final path integral over $\phi_{\alpha}(\tau)$ and using the identity \ln\det[\mathbb{A}]=\mbox{\boldmathTr}[\ln\mathbb{A}], we get,

[TABLE]

From now onwards we set $\hbar=1$. Further extracting $-\frac{1}{2}\ln\det\begin{pmatrix}0&U\delta^{c}(\tau,\tau^{\prime})\\ U\delta^{c}(\tau,\tau^{\prime})&0\end{pmatrix}$ from $\ln\mathcal{N}$ and combining with $-\frac{1}{2}\ln\det\begin{pmatrix}iP_{LL}^{0}(\tau,\tau^{\prime})&\frac{\delta^{c}(\tau,\tau^{\prime})}{U}\\ \frac{\delta^{c}(\tau,\tau^{\prime})}{U}&iP_{RR}^{0}(\tau,\tau^{\prime})\end{pmatrix}$ and using the identities, $\det\begin{pmatrix}\mathbb{I}&\mathbb{A}\\ \mathbb{B}&\mathbb{I}\end{pmatrix}=\det\left[\mathbb{I}-\mathbb{A}\mathbb{B}\right]$ and \ln\det[\mathbb{A}]=\mbox{\boldmathTr}[\ln\mathbb{A}], we get,

[TABLE]

$Tr$ in the above equation now stands only for trace over the contour time. The set of approximations made till now can be termed as mean-field dressed random-phase approximation based on the Feynman diagram representation of the final expression.

Approximate expression for $\ln\mathcal{Z}[\chi;T-T_{0}]$ given in Eq. (III.2) is valid for arbitrary measurement times $(T-T_{0})$. But, in this work we are only interested in steady-state, hence we take $(T-T_{0})\to\infty$ and neglect information contained in the initial state of the system. For solving self consistent system of equations given in Eqs. (29) and Eqs. (III.2), we approximate, $\{\phi_{\alpha}^{0}(\tau)\}$ as independent of contour time ($\{\phi_{\alpha}^{0}(\tau)\}=\{\phi_{\alpha}\}$) meaning, we assume that the stationary paths, $\{\phi_{\alpha}^{0}(\tau)\}$ as independent of time and are same on the forward and the backward branches of the contour. At this level, neglecting fluctuations of the Hubbard field, or stated equivalently, approximating the path integral within self-consistent Hartree-Fock/mean-field approximation leads to no heat flux and fluctuations at steady-state. This is because within this approximation only the second and the third terms (apart from normalization factor) which are independent of counting field are retained in Eq. (III.2). Hence fluctuations of Hubbard fields around their mean-field values are necessary to have finite heat flux and fluctuations. Within this approximation, equation for $[G_{\phi}^{c}]_{\alpha\alpha^{\prime}}(\tau,\tau^{\prime})$, Eqs. (29), is solved in the frequency domain by first projecting it onto real times (which gives four Keldysh components for each $\alpha,\alpha^{\prime}$) (notice that $[G_{\phi}^{c}]_{\alpha\alpha^{\prime}}(\tau,\tau^{\prime})\propto\delta_{\alpha\alpha^{\prime}}$ is block diagonal in orbital space) and sending all temporal integrals from $-\infty$ to $+\infty$, followed by Fourier transforming to frequency domain Rammer2007 ; Haug2008 . The solution of Eq. (29) is then given as,

[TABLE]

where $\Gamma_{\alpha}=2\pi|g_{\alpha}|^{2}\rho_{\alpha}$. With this, $[G_{\phi^{0}}]_{\alpha\alpha}(\omega)$ is a function of constant stationary paths $\{\phi^{0}\}$, which are determined self-consistently using Eq. (III.2), which in real-time read as,

[TABLE]

Expressing these equations in frequency domain using Eq. (III.2), we get,

[TABLE]

The $\omega$ integrals in the above equations can be analytically performed to get,

[TABLE]

where $\textbf{Im}\Psi[z]$ is the imaginary part of digamma function evaluated at ’$z$’ Milne1972 . Eq. (III.2) are coupled non-linear self-consistent equations for $\{\phi_{\alpha}^{0}\}$ which are difficult to solve analytically. However if we specialize to a special parameter regime $\epsilon_{\alpha}=\mu_{\alpha}-\frac{U}{2}$ and noting that $\textbf{Im}\Psi[z]=0$ for real $z$, it is clear that $\phi_{L}^{0}=\phi_{R}^{0}=\frac{U}{2}$ is always a stable solution for Eqs. (III.2), if

[TABLE]

Here onwards we confine ourselves to this regime.

We simplify the expression for $\ln\mathcal{Z}[\chi;T-T_{0}]$ given in Eq. (III.2) using the assumption that $\{\phi_{\alpha}^{0}(\tau)\}$ are independent of contour time, hence $\int_{c}d\tau\frac{\phi_{L}^{0}(\tau)\phi_{R}^{0}(\tau)}{U}=0$ and absorbing \mbox{\boldmathTr}\ln[[G_{\phi^{0}}^{c}]^{-1}], which is independent of $\chi$ into $\ln\mathcal{N}$. Expanding the logarithmic term in Taylor series, projecting on to real times and sending intermediate time integrals to $-\infty$ to $+\infty$ and going over to the frequency domain, we get the long-time expression for scaled cumulant generating function as,

[TABLE]

Here

[TABLE]

We notice that,

[TABLE]

where $\tilde{P}=P|_{\chi=0}$ and $\Lambda_{\alpha}(\omega)=e^{i\frac{\chi_{\alpha}}{2}\omega\sigma_{z}}$ (where $\sigma_{z}$ is the Pauli matrix). Further, $\tilde{P}_{\alpha\alpha}^{0}(\omega)$ can be expressed in terms of retarded, advanced and Keldysh projections of counting-field independent polarization propagators. After performing $\omega^{\prime}$ integral in Eq. (44), we get,

[TABLE]

with $\mathcal{U}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&-1\\ 1&1\end{pmatrix}$. Explicit expressions for counting field independent Keldysh rotated polarization propagators are given as : $[\tilde{P}_{\alpha\alpha}^{0}]^{A}(\omega)=-\left([\tilde{P}_{\alpha\alpha}^{0}]^{R}(\omega)\right)^{*}$ and $[\tilde{P}_{\alpha\alpha}^{0}]^{K}(\omega)=\left[1+2n_{\alpha}(\omega)\right]\left([\tilde{P}_{\alpha\alpha}^{0}]^{R}(\omega)-[\tilde{P}_{\alpha\alpha}^{0}]^{A}(\omega)\right)$ with $n_{\alpha}(\omega)=\frac{1}{e^{\beta_{\alpha}\omega}-1}$ and

[TABLE]

Using $\epsilon_{\alpha}=\mu_{\alpha}-\frac{U}{2}$ and $\phi_{L}^{0}=\phi_{R}^{0}=\frac{U}{2}$ in Eq. (47) gives,

[TABLE]

Before proceeding further, we note that $P_{\alpha\alpha}^{0}(\omega)$ given in Eqs. (47) $\&$ (III.2) is a meromorphic function with simple poles in the lower complex plane. Fourier transform (which can easily be obtained) of Eq. (47) display oscillations at characteristic frequency $\left(\epsilon_{\alpha}+\phi_{\alpha}-\mu_{\alpha}\right)$ and decay in time with rates depending linearly on $\Gamma_{\alpha}$ and $\beta_{\alpha}^{-1}$, whereas Fourier transform of Eq. (III.2) displays pure decay behavior, as $\left(\epsilon_{\alpha}+\phi_{\alpha}-\mu_{\alpha}\right)=0$.

Finally on using $P_{\alpha\alpha}^{0}(\omega)$ expressed above Eq. (45) in terms of $\Lambda_{\alpha}(\omega)$ and Keldysh rotated quantities (Eq. (46)) in Eq. (III.2) (and fixing $\ln\mathcal{N}$ by imposing normalization condition : $\ln\mathcal{Z}[\chi=0;T-T_{0}]=0$), the final expression for long-time limit scaled cumulant generating function is obtained as,

[TABLE]

with the transmission function given by,

[TABLE]

where $[\tilde{P}_{\alpha\alpha}^{0}]^{R}(\omega)$ are given in Eq. (III.2) and $n_{\alpha}(\omega)$ is the bosonic distribution function. The algebraic form of the scaled cumulant generating function given in Eq. (III.2) is similar to that of heat transport across two-terminal bosonic harmonic junctions Saito2007 ; Agarwalla2012 ; Wang2014 . Unlike harmonic junctions, the transmission function given in Eq. (50) depends on temperatures of the reservoirs. We note that similar expression for scaled cumulant generating function (Eq. (III.2)) for a variant of the model considered here Tang2018 and transmission function (Eq. (50)) for the same model Wang2018 were obtained using bare random phase approximation recently. In contrast to these works, analytical expressions for polarization functions could be obtained here by invoking wide-band assumption. Using $P_{\alpha\alpha}^{0}(\omega)$ as defined in Eq. (47) with $\phi_{L}^{0}=\phi_{R}^{0}=0$ in Eq. (50) gives transmission function obtained in Ref. (Wang2018 ) for wide-band reservoirs case.

Expressions for first four long-time limit scaled cumulants ($C_{n}=i^{n}\frac{d^{n}}{d\chi^{n}}\mathcal{F}[\chi]$) are given as,

[TABLE]

Figure (5) shows the plot long-time limit of first four scaled cumulants of heat flowing from right reservoir to the left reservoir as a function of system reservoir coupling strength ($\Gamma_{L}=\Gamma_{R}=\Gamma$). It can be seen that heat flux and fluctuations exhibit non-monotonic behavior as a function of system-reservoir coupling strength ($\Gamma_{L}=\Gamma_{R}=\Gamma$). This can be understood as follows : as noted previously electron density fluctuations on dots are solely responsible for heat flux, which are exponentially suppressed in time with rate depending linearly on $\Gamma$, hence the heat flux is suppressed for large $\Gamma$. At low temperatures ($\beta_{L}^{-1}\approx 0,\beta_{R}^{-1}\approx 0$), only low frequency behavior of transmission function is important, for small $\omega$, $\mathbb{T}(\omega)\approx\frac{16U^{2}}{\pi^{2}\Gamma^{4}}\omega^{2}$ (here $\Gamma_{L}=\Gamma_{R}=\Gamma$). Hence at low temperatures and large system-reservoir coupling strength ($\Gamma$), heat current and fluctuations decay as a power law ($\approx\Gamma^{-4}$) with system-reservoir coupling strength. We note that similar non-monotonic behavior of particle flux through double-quantum dot system is observed recently Yadalam2018 . Further using this approximate expression for transmission function in the expression for heat flux ($C_{1}$ given in Eq. (III.2)) we get the result Milne1972 $C_{1}=\frac{16\pi^{2}U^{2}}{15\Gamma^{4}}(\beta_{R}^{-4}-\beta_{L}^{-4})$ (Stefan-Boltzmann law Landau1975 ) which is already noted in Ref. (Wang2018 ).

We note that, $\mathcal{F}[-\chi-i\left(\beta_{L}-\beta_{R}\right)]=\mathcal{F}[\chi]$, which is the steady-state Gallavotti-Cohen fluctuation symmetry, this symmetry leads to the standard steady-state fluctuation theorem for the stochastic heat flux flowing from the right reservoir to the left reservoir ($P[\Delta Q;T-T_{0}]$) i.e., $\underset{(T-T_{0})\to\infty}{\lim}\frac{P[+\Delta Q;T-T_{0}]}{P[-\Delta Q;T-T_{0}]}=e^{(\beta_{L}-\beta_{R})\Delta Q}$.

Figure (6) shows $P[\Delta Q;T-T_{0}]$ and $P[-\Delta Q;T-T_{0}]$ along with the inset plot of $\ln\frac{p[\Delta Q;T-T_{0}]}{p[-\Delta Q;T-T_{0}]}$ showing the validity of steady-state Gallavotti-Cohen fluctuation theorem.

IV Conclusion

In this work we studied heat flux and fluctuations of heat flowing across capacitively coupled double quantum dot circuit. We calculated moment generating function using two theoretical approaches valid in different parameter regimes. We found using Lindblad quantum master equation that heat flux and fluctuations exhibit non-monotonic behavior as a function of coulomb interaction strength and exponentially decay for strong coulomb interaction strength. Similarly using saddle point approximation scheme for Schwinger-Keldysh path integrals, heat flux and fluctuations are found to exhibit non-monotonic behavior as a function of system reservoir coupling strength and decay as inverse fourth power of system-reservoir coupling strength for large system-reservoir coupling strength. Further we have verified that the scaled cumulant generating function obtained using both the approximation schemes has Gallavotti-Cohen symmetry and hence the steady-state fluctuation theorem for the fluctuating heat flux is verified.

Acknowledgements

H. Y. and U. H. acknowledge the financial support from the Indian Institute of Science (India).