Numerically exact counting statistics of energy current in the Kondo regime

TL;DR

This paper employs the inchworm Quantum Monte Carlo method to precisely analyze the full counting statistics of energy and particle currents in a strongly correlated quantum dot, revealing the crossover from Coulomb blockade to Kondo physics.

Contribution

It introduces a numerically exact approach to study energy current fluctuations in the Kondo regime, surpassing traditional master equation methods.

Findings

Identification of heat fluctuations and entropy production in quantum thermoelectric devices.

Observation of the crossover from Coulomb blockade to Kondo physics in energy current fluctuations.

Demonstration of the failure of conventional master equations in capturing this crossover.

Abstract

We use the inchworm Quantum Monte Carlo method to investigate the full counting statistics of particle and energy currents in a strongly correlated quantum dot. Our method is used to extract the heat fluctuations and entropy production of a quantum thermoelectric device, as well as cumulants of the particle and energy currents. The energy--particle current cross correlations reveal information on the preparation of the system and the interplay of thermal and electric currents. We furthermore demonstrate the signature of a crossover from Coulomb blockade to Kondo physics in the energy current fluctuations, and show how the conventional master equation approach to full counting statistics systematically fails to capture this crossover.

Equations

$$\hat{H}={\hat{H}}_D+{\hat{H}}_B+{\hat{H}}_V.$$

$${\hat{H}}_D$$

$${\hat{H}}_B$$

$${\hat{H}}_V$$

$${\Gamma }_{\ell }\left(\omega \right)=\pi \sum _{k\sigma }{\left|V_{k\sigma ,\ell }\right|}^2\delta \left(\omega -{\epsilon }_{k\sigma ,\ell }\right).$$

$${\Gamma }_{\ell }\left(\omega \right)=\frac{{\Gamma }_{\ell }}{\left(1+e^{\nu \left(\omega -{\Omega }_c\right)}\right)\left(1+e^{-\nu \left(\omega +{\Omega }_c\right)}\right)}.$$

$$Z\left(t;\lambda ,\chi \right)\equiv \sum _{△n_{\ell },△{\epsilon }_{\ell }}P\left(△n_{\ell },△{\epsilon }_{\ell },t\right)e^{i\lambda △n_{\ell }}{e}^{i\chi △{\epsilon }_{\ell }}.$$

$$Z\left(t;\lambda ,\chi \right)$$

$$\equiv \text{Tr}\left[T_{\mathcal{C}}{e}^{-i{\int }_{\mathcal{C}}dz{\hat{H}}_{\left(\lambda ,\chi \right)}(z)}{\hat{\rho }}_0\right],$$

$${\hat{U}}_{\left(\lambda ,\chi \right)}\left(t,t_0\right)=e^{i\frac{\chi }{2}{\hat{H}}_{\ell }}{e}^{i\frac{\lambda }{2}{\hat{N}}_{\ell }}\hat{U}\left(t,t_0\right)e^{-i\frac{\lambda }{2}{\hat{N}}_{\ell }}{e}^{-i\frac{\chi }{2}{\hat{H}}_{\ell }}$$

$${\hat{H}}_{\left(\lambda ,\chi \right)}\left(z\right)={\hat{H}}_D\left(z\right)+{\hat{H}}_B\left(z\right)+{\hat{H}}_{V\left(s_{\mathcal{C}}\lambda ,s_{\mathcal{C}}\chi \right)}\left(z\right).$$

$${\hat{H}}_{V\left(\lambda ,\chi \right)}=\sum _{k\sigma ,\ell }\left(V_{k\sigma ,\ell }{\hat{a}}_{k\sigma ,\ell }^{†}{\hat{d}}_{\sigma }{e}^{(\lambda +{\epsilon }_{k\sigma ,\ell }\chi )/2}+\text{H.c.}\right).$$

$$p\left(z_1,z_2;\lambda ,\chi \right)\equiv$$

$$p_{\varphi {\varphi }^{\prime }}\left(z_1,z_2;\lambda ,\chi \right)\equiv$$

$$Z\left(t;\lambda ,\chi \right)=\sum _{\varphi }{p}_{\varphi \varphi }\left(t^{+},t^{-};\lambda ,\chi \right).$$

$$p_{\varphi {\varphi }^{\prime }}\left(z_1,z_2;\lambda ,\chi \right)=\sum _{n=0}^{\infty }\sum _{\mathcal{C}}{w}_{\mathrm{loc}}^{\left(n\right)}\left(\mathcal{C}\right)w_{\mathrm{hyb}}^{\left(n\right)}\left(\mathcal{C};\lambda ,\chi \right).$$

$$w_{\mathrm{loc}}^{\left(n\right)}\left(\mathcal{C}\right)={\left(-i\right)}^{n_{+}-n_{-}}\prod _{i=0}^{2n-1}{p}_{{\varphi }_i{\varphi }_{i+1}}^{\left(0\right)}\left(z_i,z_{i+1}\right).$$

$$w_{\mathrm{hyb}}^{\left(n\right)}\left(\mathcal{C};\lambda ,\chi \right)=\sum _{\{\mathbf{X}\}}\mathrm{sign}\left(\mathbf{X}\right)\prod _{i=0}^n{\tilde{\Delta }}^{\lambda ,\chi }\left(z_i,z_{X_i}\right).$$

$${\widetilde{△}}_{\ell }^{\lambda ,\chi }\left(z_1,z_2\right)$$

$$+e^{-i{\lambda }_{\ell }\left(1-{\delta }_{\nu {\nu }^{\prime }}\right)}\theta \left(z_2-z_1\right){△}_{\ell }^{<}\left(z_1-z_2;\chi \right),$$

$${\Delta }_{\ell }^{>}\left(t_1-t_2;\chi \right)$$

$$\times {\Gamma }_{\ell }\left(\omega \right)\left[1-f\left(\omega -{\mu }_{\ell }\right)\right],$$

$${\Delta }_{\ell }^{<}\left(t_1-t_2;\chi \right)$$

$$\times {\Gamma }_{\ell }\left(\omega \right)f\left(\omega -{\mu }_{\ell }\right).$$

$${\Delta }_{\ell }^{\lessgtr }\left(t_1-t_2;\chi \right)={\Delta }_{\ell }^{\lessgtr }\left(t_1-t_2\mp {\chi }_{\ell }\left(1-{\delta }_{\nu {\nu }^{\prime }}\right);0\right)$$

$$Z\left(t;\lambda ,\chi \right)={\mathrm{Tr}}_S\left[\hat{\sigma }\left(t;\lambda ,\chi \right)\right]=\sum _{\varphi }{p}_{\varphi }\left(t;\lambda ,\chi \right).$$

$$\frac{\partial \mathbf{p}\left(t;\lambda ,\chi \right)}{\partial t}=M\left(\lambda ,\chi \right)\mathbf{p}\left(t;\lambda ,\chi \right).$$

$$M_{\varphi {\varphi }^{\prime }}=\sum _{\ell }{M}_{\varphi {\varphi }^{\prime }}^{\ell },$$

$$M_{\phi\phi^{\prime}}^{\ell}=\left\{\begin{array}[]{cc}\underset{\phi^{\prime\prime}}{\sum}\Gamma_{\ell}f_{\ell}\left(\Delta E_{\phi^{\prime}\phi^{\prime\prime}}\right)&\phi=\phi^{\prime},\\ -\Gamma_{\ell}f_{\ell}\left(\Delta E_{\phi^{\prime}\phi}\right)&n_{\phi^{\prime}}-n_{\phi}=1,\\ -\Gamma_{\ell}(1-f_{\ell}\left(\Delta E_{\phi\phi^{\prime}}\right))&n_{\phi^{\prime}}-n_{\phi}=-1,\\ 0&n_{\phi^{\prime}}-n_{\phi}=\pm 2,0.\end{array}\right.$$

$$M_{\varphi {\varphi }^{\prime }}^{\ell }\to M_{\varphi {\varphi }^{\prime }}^{\ell }{e}^{i\lambda {ϵ}_{\varphi {\varphi }^{\prime }}}{e}^{-i\chi \Delta E_{{\varphi }^{\prime }\varphi }}.$$

$${\mathcal{S}}_{\ell }\left(t;\lambda ,\chi \right)=\log Z\left(t;\lambda ,\chi \right),$$

$$C_{\ell ,m,k}\left(t\right)={\frac{{\partial }^{m+k}{\mathcal{S}}_{\ell }\left(t;\lambda ,\chi \right)}{\partial {\left(i\lambda \right)}^m\partial {\left(i\chi \right)}^k}|}_{\lambda ,\chi =0}.$$

$$C_{\ell ,n}^p\left(t\right)\equiv C_{\ell ,n,0}\left(t\right),$$

$$C_{\ell ,n}^E\left(t\right)\equiv C_{\ell ,0,n}\left(t\right).$$

$$C_{\ell }^{\left(\times \right)}\left(t\right)\equiv C_{\ell ,1,1}\left(t\right).$$

Full text

I Introduction

Energy transport through small junctions is a major paradigm at the heart of attempts to formulate thermodynamic principles applicable to nanoscale quantum systems. Allahverdyan and Nieuwenhuizen (2001); Martinez and Paz (2013); Ludovico et al. (2014); Esposito et al. (2015a); Bruch et al. (2016); Ochoa et al. (2016); Haughian et al. (2018) The latter are essential for understanding and improving operational efficiency in nanoelectronic devices,Pop (2010) where quantum entanglement effects are inseparably intertwined with performance fluctuations. Kosloff and Levy (2014); Uzdin et al. (2015); Esposito et al. (2015b); Niedenzu et al. (2018)

Fluctuations in the energy and particle currents Blanter and Büttiker (2000); Pietzonka and Seifert (2018) can be studied through consideration of the full counting statistics (FCS) of particle and energy transfer. Levitov et al. (1996); Nazarov (2007); Esposito et al. (2007); Bednorz and Belzig (2008) A particularly challenging problem in this field is the study of electronic correlations and their impact on the performance of quantum devices coupled to conducting reservoirs. A schematic view of a junction is shown in Fig. 1. This comprises a central molecule or quantum dot where confined electrons interact strongly, contacted by weakly interacting metallic left ($L$) and right ($R$) reservoirs. A voltage bias is applied across the junction, and an additional thermal or temperature bias may either not be applied (left panel) or applied (right panel).

Under certain conditions, the system may be characterized by a threshold known as the Kondo temperature $T_{K}$. Below this threshold, transport properties are dominated by the formation of a correlated resonance in the conductance spectrum.Haldane (1978); Hewson (1993) A bias voltage (or chemical potential difference between the leads) splits the Kondo resonance into a pair of peaks centered near the lead chemical potentials. Meir et al. (1993); Cronenwett (1998); Schmidt and Wölfle (2010) The resonance width is set by $T_{K}$, and it is typically much sharper than features related to noninteracting resonant tunneling, Ng and Lee (1988) as seen in experiments on strongly correlated single-molecule transistors. Cronenwett (1998); Liang et al. (2002)

A thermal bias (right panel of Fig. 1) may be applied to the junction by keeping the two reservoirs at different temperatures $T_{L}$ and $T_{R}$. The interplay between electric and thermal biases can be used to implement e.g. heat engines, heat pumps and refrigerators, depending on the direction of the thermal gradient and the relative direction of flow of particles $I_{p}$ and heat $I_{h}$. Kosloff and Levy (2014); Whitney (2015); Benenti et al. (2017); Whitney et al. (2018)

Recent experiments on molecular junctions show that both the Seebeck and Peltier coefficients can be measured at the nanoscale Reddy et al. (2007); Lee et al. (2013); Brechet et al. (2013); Kim et al. (2014); Rincón-García et al. (2016); Cui et al. (2017a, 2018) and the quantized character of thermal transport in nanojunctions was demonstrated. Schwab et al. (2000); Cui et al. (2017b) In addition, a recent measurement of thermal gradient induced particle noise was reported in the literature. Lumbroso et al. (2018)

In the linear regime, quantum coherence universally weakens heat engine performance, Brandner et al. (2017) and reductions in power and efficiency have also been reported in nonlinear quantum heat engines. Karimi and Pekola (2016) However, there is no comprehensive answer to the question of whether quantum effects improve or reduce the performance of nonlinear heat engines at the nanoscale. In certain specific cases, quantum coherence effects have been shown to increase the thermal efficiency, Chen et al. (2017a) to exceed the Carnot efficiency and approach a perfect efficiency of 1, Toyabe et al. (2010) to increase the power output Scully et al. (2011) and to reduce fluctuations in the power. Ptaszyński (2018)

Whereas it is possible to study nanoscale thermoelectricity within a single-particle approximation with Green’s function Crépieux et al. (2011); Ludovico et al. (2016); Covito et al. (2018) or DFT-based Pauly et al. (2008); Zotti et al. (2014); Eich et al. (2014) methods, the effect of Kondo physics on the thermoelectric performance of devices is not fully understood beyond linear response. An exception is recent work based on the matrix product states method, which gives access to the current generated by a temperature bias in the Kondo regime. Dorda et al. (2016) In related work, simulations of the current within the non-crossing approximation show a strong enhancement in the Seebeck coefficient by asymmetric dot–lead couplings in the Kondo regime. Pérez Daroca et al. (2018)

Going beyond the study of heat and charge currents, one may investigate higher-order correlations in the transferred charge and energy. These higher order statistical cumulants, and their underlying probability distributions, are obtained from the cumulant generating function for the FCS, and can reveal a whole range of information not present in the first moment. Landauer (1998) This includes particle traversal times, Ridley et al. (2017) waiting time distributions Tang et al. (2014); Rudge and Kosov (2016); Ptaszyński (2017); Ridley et al. (2018); Kosov (2018) shot noise Büttiker (1990); Levitov et al. (1996); Blanter and Büttiker (2000) and associated quasiparticle charges. Yamauchi et al. (2011) In addition, the steady state cumulant generating function can be related to fluctuations in the thermal efficiency. Esposito et al. (2015b)

To date, the majority of the exact FCS theory has focused on noninteracting models of nano junctions, beginning with the work of Levitov and Lesovik on the counting statistics of charge transport in the steady state regime, which reduces to the Landauer–Büttiker theory of shot noise, in which the propagating charge carriers are assumed to be coherent waves in the device region. Landauer (1998); Blanter and Büttiker (2000); OuYang et al. (2018) This program culminated with the path integral nonequilibrium Green’s function (PINEGF) approach pioneered by Tang et al.,  Tang and Wang (2014); Tang et al. (2014); Yu et al. (2016); Tang et al. (2017, 2018) valid for FCS calculations in the transient regime following a quench and for arbitrary time-dependent driving fields.

In strongly correlated systems, the quantum master equation (QME) approach has proven popular for the investigation of charge Flindt et al. (2009, 2010); Albert et al. (2012) and energy Esposito et al. (2007); Clerk (2011); Sánchez and Büttiker (2012); Silaev et al. (2014); Ptaszyński (2018) FCS in Coulomb-blockaded systems. It lends itself well to analytical investigation, so it has been used to derive exact relations for optimized thermal efficiencies in molecular heat engines Esposito et al. (2009a); Abah and Lutz (2014) and exact fluctuation relations in the counting statistics of heat exchange. Friedman et al. (2018) However, the QME method is generally limited to the weak impurity–bath coupling regime: $\Gamma$, the typical strength of hybridization between the lead states and those of the molecule or dot, must be substantially smaller than all other energy scales in the problem. In recent years, the QME method was extended beyond the sequential tunneling Coulomb-blockade regime to include spin-flip scattering cotunneling processes, so that noise and FCS can be studied at higher values of $\Gamma$. Thielmann et al. (2005); Leijnse and Wegewijs (2008); Kaasbjerg and Belzig (2015); Seja et al. (2016); Walldorf et al. (2017)

However, the Kondo regime and current noise resulting from the formation of a Kondo peak  Hewson (1993); Cohen et al. (2014a) are inaccessible to perturbative QME approaches. Miwa et al. (2017) Indeed, the study of current fluctuations in the Kondo regime has mostly been carried out by methods specialized to very specific parameter regimes. Sela et al. (2006); Gogolin and Komnik (2006); Carr et al. (2011) The situation has changed with the recent development of a numerically exact method for the computation of full particle counting statistics (FCS). Ridley et al. (2018) In this context, ‘numerically exact’ refers to the controlled evaluation of a quantity with known error bounds that can be made as small as desired with enough computer time. This approach was based on the inchworm quantum Monte Carlo (iQMC) algorithm, Cohen et al. (2015); Antipov et al. (2017); Chen et al. (2017b, c) a powerful method for the stochastic evaluation of propagators in strongly correlated impurity models. It circumvents the dynamical sign problem that plagues conventional QMC methodsGull et al. (2011) and scales polynomially in time. Cohen et al. (2015) In Ref. Ridley et al., 2018, numerically exact calculations of the current noise were presented, showing a crossover from interaction-induced suppression to enhancement as the temperature was lowered below $T_{K}$. This finding corroborates the experimental data showing Kondo-enhanced noise in strongly interacting impurity models,  Delattre et al. (2009) and also indicates the failure of noninteracting theory to properly account for noise in the Kondo regime.

In the present study we extend the iQMC FCS method to the full counting statistics of nonequilibrium energy transfer. We treat junctions with electron–electron interactions, thus taking steps towards a sorely-needed quantum theory of thermoelectricity. Whitney et al. (2018) Our methodology enables numerically exact calculations of cumulants of charge, energy and heat transfer in strongly correlated regimes, as well as cross-correlations between particle and energy transfer. Through these quantities, we obtain access to the entropy production.

The structure of the paper is as follows: in Section II we introduce the model of the quantum junction utilized in this study. Section III introduces the techniques we use (both numerically exact iQMC and approximate QME) to calculate particle and energy transfer statistics in the transient regime, and Section IV applies this to quantities of interest in the study of thermoelectric properties of junctions. In Section V we present our main results. This includes a comparison with QME data, a discussion of the effects of interaction and temperature on heat cumulants and associated heat current fluctuations, and calculations of entropy production. Finally, our conclusions can be found in Section VI.

II The Model

We consider a junction described by the nonequilibrium Anderson impurity model (AIM) Hamiltonian,

[TABLE]

Here, $\hat{H}_{D}$ describes a dot (or molecule) region with on-site interactions; $\hat{H}_{B}=\sum_{\ell}{\hat{H}_{\ell}}$ describes two large noninteracting reservoirs or leads, denoted by $\ell\in\{L,R\}$; and $\hat{H}_{V}$ is a hybridization Hamiltonian coupling the dot to both leads. These terms are given by

[TABLE]

The $\hat{d}^{\dagger}_{\sigma}$ ($\hat{d}_{\sigma}$) operators create (annihilate) an electron with spin $\sigma\in\left\{\ \uparrow,\downarrow\right\}$ on the dot, and the $\hat{a}_{k\sigma,\ell}^{\dagger}$ ($\hat{a}_{k\sigma,\ell}$) perform the analogous operation on the single particle level $k$ of lead $\ell$. $\varepsilon_{\sigma}$ and $\varepsilon_{k\sigma,\ell}$ are single particle energies on the dot and lead levels, respectively, and $U$ sets the strength of Coulomb repulsion on the dot. Throughout this work, we impose particle–hole symmetry by setting $\epsilon_{\sigma}=-\frac{U}{2}$. Finally, the $V_{k\sigma}$ control tunneling between the dot and leads.

We will simulate dynamics starting from an initial state $\rho_{0}=\rho_{L}\otimes\rho_{D}\otimes\rho_{R}$, where each lead $\ell$ is in a thermal state characterized by a temperature $T_{\ell}=1/\beta_{\ell}$ and a chemical potential $\mu_{L/R}=\pm V/2$. We assume that the dot is prepared in one of the four eigenstates of $\hat{H}_{D}$: the empty state $\left|0\right\rangle$, full state $\left|\uparrow\downarrow\right\rangle$, or one of the two magnetized half-occupied states $\left|\sigma\right\rangle$. In what follows these local or ‘atomic’ states are referred to collectively by the index $\phi$. At time $t=0$, the system begins to evolve with respect to $\hat{H}$, corresponding to instantaneously adding the hybridization term $\hat{H}_{V}$ to the Hamiltonian to model a sudden coupling of the leads to the dot. This setup is often referred to as either a coupling quench or a partitioned approach to the switch-on problem. Ridley and Tuovinen (2018)

The $\varepsilon_{k\sigma,\ell}$ and $V_{k\sigma,\ell}$ are effectively set by the coupling density $\Gamma_{\ell}\left(\omega\right)$ that describes the electron escape rate of lead $\ell$:

[TABLE]

We take this to be a flat band with a smooth cutoff,Werner et al. (2009)

[TABLE]

In what follows, we set $\Gamma_{\ell}=\frac{1}{2}$ such that $\Gamma\equiv\underset{\ell}{\sum}\Gamma_{\ell}=1$ determines our unit of energy. We take the leads’ band cutoff $\Omega_{c}$ to be $10\Gamma$, and their edge width $\frac{1}{\nu}$ to be $\Gamma/10$.

III FCS of energy and particles

A chemical potential bias and/or a temperature difference between the leads engenders charge and energy fluxes across the junction. For the interface between any lead $\ell$ and the dot, FCS yields a way to study the statistics of both electron transfer, $\Delta\hat{N}_{\ell}\left(t\right)\equiv\hat{N}_{\ell}\left(t\right)-\hat{N}_{\ell}\left(0\right)$; and energy transfer, $\Delta\hat{H}_{\ell}\left(t\right)\equiv\hat{H}_{\ell}\left(t\right)-\hat{H}_{\ell}\left(0\right)$. Here $\hat{N}_{\ell}\left(t\right)$ and $\hat{H}_{\ell}\left(t\right)$, respectively, are operators corresponding to the total number of electrons and total energy in lead $\ell$. In the following section, we briefly introduce some of the main concepts and definitions of FCS for the convenience of the reader.

III.1 FCS from iQMC

Within the FCS formalism, the statistics of currents at the interface with lead $\ell$ are obtained from a two-point measurement of the total number of particles or energy within reservoir $\ell$ at times [math] and $t$. Multiple realizations of the quantum process yield information on the probability $P\left(\Delta n_{\ell},\Delta\epsilon_{\ell},t\right)$ to observe transfer of $\Delta n_{\ell}$ electrons and $\Delta\epsilon_{\ell}$ energy across the interface. These probabilities are encoded within the generating function Esposito et al. (2009b); Tang and Wang (2014)

[TABLE]

The Fourier variables $\lambda$ and $\chi$ are called counting fields. In general there may be a separate counting field for every lead in a multiterminal counting experiment; the generalization of Eq. (7) to the multiterminal case is straightforward.

Assuming the initial state of the system to be an eigenstate of the decoupled dot Hamiltonian, Eq. (7) can be expressed in terms of counting-field-modified propagators

[TABLE]

where $\mathcal{C}$ is the Keldysh contour comprising a forward branch $\mathcal{C}_{+}$ (from $t_{0}$ to $t$) followed by a backward branch $\mathcal{C}_{-}$ (from $t$ to $t_{0}$), $z$ is a contour-time variable, and $T_{\mathcal{C}}$ is the time ordering operator on the contour. Time evolution is governed by the modified propagator and Hamiltonian, respectively given by

[TABLE]

and

[TABLE]

Here, $s_{\mathcal{C}}=\pm 1$ at the $\mathcal{C}_{\pm}$ branch of the contour and

[TABLE]

The object evaluated in the iQMC method is the auxiliary-field-modified propagator between two points $z_{1}$ and $z_{2}$ on the contour,

[TABLE]

where the second line describes projection to the atomic $\phi$ basis.

The propagator of Eq. (12) is analogous to propagators used in other real time, continuous time hybridization-expansion Monte Carlo methods: whereas earlier implementations only implicitly relied on such propagators,Mühlbacher and Rabani (2008); Werner et al. (2009); Schiró and Fabrizio (2009); Schiró (2010); Cohen and Rabani (2011) later work explicitly took advantage of their properties.Cohen et al. (2015); Antipov et al. (2017); Chen et al. (2017b, c); Dong et al. (2017); Boag et al. (2018); Krivenko et al. (2019); Gull et al. (2011); Cohen et al. (2013, 2014a, 2014b); Kubiczek et al. (2019) Here, the hybridization Eq. (11) is modified with respect to its physical counterpart $\hat{H}_{V\,\left(\lambda=0,\chi=0\right)}$ by $\lambda$- and $\chi$-dependent factors. The first of these modifications was already present in the previous iQMC FCS paper,Ridley et al. (2018) where it was used to address particle (but not energy) counting statistics.

The bare Monte Carlo process associated with evaluating Eq. (12) in the hybridization expansion results in a dynamical sign problem, such that the computational cost of the simulation increases exponentially with time. The inchworm algorithmCohen et al. (2015) overcomes this for at least some parameters, such that the FCS can be obtained efficiently.Ridley et al. (2018) In the present work, no new conceptual difficulties beyond the introduction of a second counting field (see below) arise. We therefore refer the reader to the literature for a thorough discussion of how iQMC can be applied to particle counting statistics.Ridley et al. (2018) In the rest of this chapter, we limit ourselves to a self-contained description of the the modifications to the theory underlying iQMC, that are needed in order to also access energy counting statistics.

The cumulant generating function at time $t$ can be extracted directly from the propagator between the contour times $t^{+}$ and $t^{-}$ corresponding to the physical time $t$ on the two branches of the Keldysh contour:

[TABLE]

By expanding Eq. (12) in powers of $\hat{H}_{V}$, the propagator can be represented as a sum over configurations suitable for stochastic Monte Carlo sampling,

[TABLE]

Since only even orders contribute in this combination of expansion, model and propagator, we use $n$ to denote the expansion order, but a term of order $n$ actually contains $2n$ hybridization vertices of the form in Eq. (11). These vertices occur at times $\mathcal{C}=\{\zeta_{1},\ldots,\zeta_{2n}\}$. The vertex times, all of which are located within the part of the contour between $z_{1}$ and $z_{2}$, are integrated over in the second summation.

The $w_{\mathrm{loc}}^{(n)}$ are products of interacting (but purely local) atomic state propagators $p_{\ell\kappa}^{\left(0\right)}\left(z_{1},z_{2}\right)$, defined as in Eq. (12) but with $\hat{H}_{\left(\lambda,\chi\right)}$ replaced by $\hat{H}_{0}\equiv\hat{H}_{D}+\hat{H}_{B}$ in the time-ordered integral. The states $\phi_{i}$ between each pair of vertices are fully determined by $\mathcal{C}$ and the initial condition $\phi_{0}$. Setting $\phi_{2n}=\phi_{0}$ allows us to write these objects in the following form, which is independent of both $\lambda$ and $\chi$:

[TABLE]

The $w_{\mathrm{hyb}}^{(n)}$, on the other hand, explicitly depend on the counting fields:

[TABLE]

In general, $\mathbf{X}$ are all permutations of the integers $(1,\ldots,n)$, and $\tilde{\Delta}^{\lambda,\chi}=\sum_{\ell}\tilde{\Delta}^{\lambda_{\ell},\chi_{\ell}}_{\ell}$, but below we always set $\lambda_{\ell}=\lambda\delta_{\ell L}$ and $\chi_{\ell}=\chi\delta_{\ell L}$. The contribution to the counting-field-modified hybridization from lead $\ell$ is given by

[TABLE]

where $\theta$ is the Heaviside function on the contour and $\nu,\nu^{\prime}\in\left\{\mathcal{C}_{+},\mathcal{C}_{-}\right\}$ are the branch indices of $z_{1}$ and $z_{2}$, respectively. Finally, the $\lambda$-unmodified lesser and greater hybridization components for lead $\ell$ are obtained from the level width function of Eq. (5):

[TABLE]

Here $t_{1}$ and $t_{2}$ are physical times corresponding to the contour times $z_{1}$ and $z_{2}$. For leads initially in equilibrium or steady state, the two-time hybridization functions depend only on the time difference $t_{1}-t_{2}$, and it is convenient to use the Fourier shift property

[TABLE]

to avoid numerical issues with strongly oscillating frequency integrals at finite values of $\chi$.

III.2 FCS from QME

The quantum master equation approach to FCS provides an approximate expression for the dynamics of the reduced density matrix modified by the counting field, $\hat{\sigma}\left(t;\lambda,\chi\right)\equiv\mathrm{Tr}_{B}\left\{\hat{\rho}\left(t;\lambda,\chi\right)\right\}$. Bagrets and Nazarov (2003)

Since the dynamics of the off-diagonal elements of $\hat{\sigma}$ is completely decoupled from that of the diagonal ones in the AIM, it is sufficient to consider the diagonal elements of $\hat{\sigma}\left(t;\lambda,\chi\right)$, the counting field-modified “populations” $p_{\phi}\left(t;\lambda,\chi\right)$. The generating function $Z\left(t;\lambda,\chi\right)$ is then given by the trace with respect to the dot subsystem Esposito et al. (2009b)

[TABLE]

We collect the $p_{\phi}$ into a vector $\mathbf{p}$, such that—within the QME approximation—$\mathbf{p}$ satisfies the following rate equation:

[TABLE]

Here, $M$ is a matrix with elements given at $\lambda,\chi=0$ by Datta (2005); Nitzan (2013); Levy et al. (2019)

[TABLE]

where

[TABLE]

The terms $M_{\phi\phi^{\prime}}^{\ell}$ in this expression correspond to transition rates from dot state $\phi^{\prime}$ to state $\phi$ as mediated by the lead $\ell$, and we have defined $n_{\phi}$ to be the occupation of the decoupled dot state $\phi$, $\Delta E_{\phi^{\prime}\phi}\equiv E_{\phi^{\prime}}-E_{\phi}$ and $f_{\ell}(\epsilon)=1/(1+e^{\beta_{\ell}(\epsilon-\mu_{\ell})})$. The level width $\Gamma_{\ell}$ is taken to be frequency independent, an additional approximation which could be easily removed, but is justified for the band width and shape used here. Covito et al. (2018); Ridley et al. (2019)

In the presence of the counting field, the transition rates undergo the modification

[TABLE]

Here, $\epsilon_{ij}=\pm 1$ for $n_{i}=n_{j}\pm 1$, and $\epsilon_{ij}=0$ otherwise. The sign of the exponential factor $-i\chi\Delta E_{\phi^{\prime}\phi}$ arises from the fact that an increase of energy in the dot corresponds to a decrease in lead $\ell$. Schaller et al. (2013); Agarwalla et al. (2015)

We note that more sophisticated QME approaches exist which take into account, e.g., higher-order tunneling processes. Kaasbjerg and Belzig (2015) However, Markovian master equations generally fail to address the tunneling-induced level broadening of dot states, and the accurate evaluation of the deeply non-Markovian memory kernels that characterize correlated regimes requires a numerical method on par with iQMC. Cohen and Rabani (2011); Cohen et al. (2013); Härtle et al. (2015)

IV Cumulants and thermoelectric quantities

The cumulants of the statistical distribution for particle and energy transfer at the interface of the dot with lead $\ell$, as well as their cross-correlations, may be obtained from derivatives of the cumulant generating function,

[TABLE]

with respect to the counting fields $\lambda$ and $\chi$ and with the convention introduced above that these fields are nonzero only for properties stemming from lead $\ell$. In general, we define

[TABLE]

In terms of these derivatives, we may define the $n$-th order cumulants of the separate particle and energy transfer processes

[TABLE]

We may also consider the cross-correlations measuring the statistical influence of particle change on a change in the energy:

[TABLE]

One may then evaluate the first and second cumulants of transferred heat from particle, energy and cross-correlation statistics Kilgour et al. (2019)

[TABLE]

In the long time regime, the first and second cumulants yield fluxes, $I_{\ell}^{p/E/h}$, and zero-frequency quantum noises, $S_{\ell}^{p/E/h}$, for particle, energy and heat transfer: Esposito et al. (2015b)

[TABLE]

Given the quantities above, the overall steady state entropy production rate can be evaluated from the expression Deffner and Lutz (2011); Whitney (2013); Pietzonka and Seifert (2018)

[TABLE]

We note that this relation is only valid for the steady state, Esposito et al. (2015a) and the second law of thermodynamics necessarily implies that $\sigma>0$.

V Results

We consider a biased junction, $\mu_{L}>\mu_{R}$, both with ($T_{L}<T_{R}$) and without ($T_{L}=T_{R}$) a temperature gradient. In the former case, at appropriately tuned parameters, the system may embody either a heat engine (electron flow against the voltage gradient is driven by the temperature bias) or a heat pump/refrigerator (heat flow against the thermal gradient is driven by the voltage bias).

V.1 $U=0$ benchmarks and failure of master equations

We begin with benchmarks of the iQMC method in the noninteracting limit $U=0$ against exact path-integral nonequilibrium Green’s function (PINEGF) results. Tang and Wang (2014); Yu et al. (2016) The iQMC method is numerically exact and one might trivially expect to simply find agreement. Nevertheless, such benchmarks are important whenever a numerical method is generalized, and the iQMC method has not previously been applied to the calculation of energy counting statistics. We furthermore examine the applicability of the QME approach to FCS in the noninteracting limit.

Fig. 2 displays the real (top panels) and imaginary (bottom panels) components of the moment generating function $Z\left(t;\lambda,\chi\right)$ on lead $L$, with counting fields set to $\lambda=0.5$ and $\chi=0.1$. We note that typically one is most interested in the derivatives of the generating function at the $\lambda,\chi\rightarrow 0$ limit, where they give the exact cumulants. However, logarithmic derivatives of $Z\left(t;\lambda,\chi\right)$ vary slowly with the counting fields, and the finite values of $\lambda$ and $\chi$ chosen here are therefore physically relevant and easier to work with numerically. Ridley et al. (2018) The bias voltage is $V=1\Gamma$ and the initial condition is the empty state $\left|0\right\rangle$. There is no thermal bias, i.e. $T_{L}/T_{R}=1.0$. The left panels are at a higher temperature $\beta_{L}=\beta_{R}=0.4/\Gamma$, whereas the right panels are at a lower temperature $\beta_{L}=\beta_{R}=50/\Gamma$. Each panel directly compares results from iQMC (black) with QME data (dashed blue) and exact PINEGF data (dotted red).

To within numerical errors, the iQMC and PINEGF methods agree at all times and both temperatures. Notably, the noninteracting limit is a nontrivial benchmark for the iQMC algorithm employed here, which is based on an expansion in the dot–lead hybridization rather than the interaction. As might be expected, the QME fails to correctly capture the time evolution in both cases, and fails qualitatively even at describing the steady state (constant time derivatives at long times) at the lower temperature. Braggio et al. (2006)

In Fig. 3, we focus on the higher temperature regime in which the QME is expected to work well, at $\beta_{L}=0.4/\Gamma$. The left panels repeat the calculation in the left panels of Fig. 2, but with an applied temperature bias $T_{L}/T_{R}=0.1$, so that heat is driven against the direction of the particle current. Note that this is set up so that the average temperature is higher than $\beta_{L}$, a supposedly even more favorable condition for QME. Interestingly, however, QME fails to capture the high temperature steady state not just quantitatively, but even in terms of the sign of the long-time derivative (slope) of $\operatorname{Im}{Z}$ (see lower left panel of Fig. 3).

The right-hand panels of Fig. 3 show that in this temperature regime the three methods are in excellent agreement for the case of zero energy counting field, $\chi=0$. We have verified that this remains true in the thermally unbiased case. Mixed particle–energy counting therefore appears to be more challenging to capture within the QME framework than particle counting alone.

V.2 Finite-$U$ FCS in QME and iQMC

Having established that the iQMC method provides reliable FCS in the presence of both the particle and energy counting fields, we continue to explore the interacting case, where no analytical results are available. We concentrate on energy and particle cumulants up to second order, in addition to energy–particle cross correlations $C_{1/2,\ell}^{p/E}\left(t\right)$ and $C_{pE,\ell}^{\left(\times\right)}\left(t\right)$. Using the symmetry properties $Z\left(t;\lambda,\chi\right)=Z\left(t;-\lambda,\chi\right)^{*}$ and $Z\left(t;\lambda,\chi\right)=Z\left(t;\lambda,-\chi\right)^{*}$, all these quantities can be extracted from finite-difference logarithmic derivatives based on just four simulations of the generating function in the $\left(\lambda,\chi\right)$ plane.

In the rest of this section, dashed lines in the figures will denote QME data and solid lines or symbols denote iQMC results. In Figs. 4–8, three different initial states of the dot are shown as color-coded curves: empty (black), half-filled (red) and fully-occupied (blue).

Figs. 4 and 5, respectively, show the first cumulants of particle and energy transfer on lead $\ell=L$ (left panels) and $\ell=R$ (right panels), where the counting fields are introduced identically for each lead.

We set the interaction strength to $U=8\Gamma$, the bias voltage to $V=1\Gamma$ and the temperatures to be equal on both leads, such that currents are driven only by the voltage. The approximate Kondo temperature of the system can be estimated as $k_{B}T_{K}\approx 0.1\Gamma$, in accordance with the Bethe ansatz formula $k_{B}T_{K}=\sqrt{U\Gamma/2}\exp\left(-\pi U/8\Gamma+\pi\Gamma/2U\right)$. Hewson (1993) We note in passing that recent work proposes nonequilibrium logarithmic corrections to this formula, Li et al. (2017) but these corrections are small enough that we are able to safely find temperatures above and below $T_{K}$ in what follows. In particular, the upper panels of Figs. 4–5 are at a high temperature $T=2.5\Gamma\approx 25T_{K}$, whereas the lower panels are at a low temperature $T=0.02\Gamma\approx T_{K}/5$ is shown.

Fig. 4 shows the total change in electron number with increasing time, from the quench time up until $t=8/\Gamma$. The upper panels shows a remarkable agreement between the iQMC and QME data at high temperatures, with the exception of the short-time regime in which higher-order scattering processes are dominant. Antipov et al. (2016) This plot shows significant initial condition dependence in the transferred charges, which may be understood in terms of Coulomb blockade physics: (i) in the initially unoccupied case (black line), one electron can tunnel onto the dot from the left lead, whose Fermi level is at $\mu_{L}=\Gamma/2$, before the Coulomb repulsion slows down the rate of charge transfer out of this lead, explaining the change to a less negative gradient in $C_{1,L}^{p}\left(t\right)$ after a relaxation time of approximately $\tau_{r}\approx 1/\Gamma$. (ii) When the dot is initially fully occupied (blue line), the particle number may only change on the left lead on the timescale of $\tau_{r}$ by tunneling up the voltage gradient to put the dot in a magnetized state, at which point the Coulomb blockade forces the current to take the same value as in case (i). This tunneling process is activated by the high temperature $k_{B}T>eV$ of the junction. (iii) When the dot is initially in the spin-up or spin-down state (red line), the steady state population is immediately established at the quench time, explaining why the gradient of this cumulant is identical to the gradient which only forms after time $\tau_{r}$ has elapsed in cases (i) and (ii).

In the lower panels of Fig. 4, the first cumulants are shown for the low temperature ($\beta=50/\Gamma$) case with an otherwise unchanged set of parameters. Here we see a clear temperature-dependent transition in the QME data, as all three dashed lines relax to a constant (zero current) value in a clear signature of Coulomb blockade preventing transport through the dot. However, the iQMC data shows an enhanced current magnitude. This is indicative of the additional transport channels enabled by the formation of a Kondo resonance in the nonequilibrium conductance, and illustrates the breakdown of our naive QME approach at low temperatures.

We now turn to the first cumulant of transferred energy for the same system parameters, which is displayed in Fig. 5. The QMC data clearly shows that the energy of the left lead undergoes an increase for all initial conditions. This energy gain is uniform across the entire junction, as can be seen by considering the cumulants in the right lead, shown in the right-hand panels of Fig. 5. It is therefore not associated with the direction of flow of particle current, which has a different sign in each lead. Instead, this energy increase can be interpreted as a charging effect, i.e. it arises from the change in the electrostatic energy driven by the addition or removal of an electron from the dot.

We point out that our data on energy cumulants (and later cross correlations) is significantly noisier than the corresponding results for particle cumulants. The main reason is that in order to obtain converged results for these cumulants from the logarithmic finite difference derivatives, relatively small values of $\chi$ were needed. Since the absolute value of the differences from the $\chi=0$ values is therefore small, this results in large relative errors when evaluating numerical logarithmic derivatives with respect to $\chi$. We do not include rigorous uncertainty estimates, which are expensive to obtain in this case.Cohen et al. (2015); Antipov et al. (2017) As a rough guideline, by comparing uncorrelated runs (not shown) we find that the data is reliable to within a small factor (say,  2–3) of the apparent stochastic noise.

Whereas the total transferred energy at low temperature is equal for the initially fully occupied and empty states, when the temperature is raised this symmetry is broken. This can be seen in the gap that opens up between the blue and black lines in Fig. 5 in both the iQMC and the QME data. Also for high temperatures, the symmetry between the energies transferred to the left and right leads is broken, such that the order of the blue and black lines is reversed between $L$ and $R$. To understand this, it is enough to consider the rates for dot–lead transfer processes in the $L$- and $R$-lead cases separately. In order to magnetize the dot when the initial condition is $\left|0\right\rangle$ or $\left|\uparrow\downarrow\right\rangle$, an energy of $-4\Gamma$ must be transferred to the dot from one of the lead states $k$. In the initially empty state, a particle is transferred out of one lead $\ell$ and the corresponding transfer rate $R_{\ell}$ is proportional to the Fermi function $f\left(-U/2-\mu_{\ell}\right)$. In the initially fully occupied state the rate is proportional to the hole occupation $f\left(U/2-\mu_{\ell}\right)$. At infinite $k_{B}T$, $R_{L}=R_{R}=1/2$ for both initial conditions, so there is no symmetry breaking in the system. In the opposite regime of $k_{B}T=0$, $R_{L}=R_{R}=1$ for the empty state and $R_{L}=R_{R}=0$ for the full state, once again symmetric. However, at intermediate temperatures $k_{B}T\sim\left|\triangle E\pm V\right|$, the left–right symmetry is broken. In this case, $R_{L}>R_{R}$ for the empty state and $R_{L}<R_{R}$ for the full state. If the system is initially empty, energy transfer to the left lead is favored in the short time regime. The converse is true for the full initial state.

The QME correctly captures the qualitative, though not the quantitative, aspects of the high and intermediate temperature energy transfer due to charge dynamics. Notably, however, it fails to predict even the direction of overall energy transfer in the initially magnetized case. This can be seen by comparing the dashed and solid red curves in the top panel of Fig. 5. Less surprisingly, QME also spuriously predicts no energy transfer in this case at low temperatures. This is because at the level of QME, the system essentially begins in its steady state under such conditions.

In Fig. 6 we present the time evolution of the statistical correlation between the changes in particle number and energy in each lead, $C_{\ell}^{\left(\times\right)}\left(t\right)=\left\langle\triangle\epsilon_{\ell}\left(t\right)\triangle n_{\ell}\left(t\right)\right\rangle$. Immediately one observes that the sign of this quantity depends on the initial condition, making it an excellent observable for experimental detection of the initial system preparation. In the case of an initially occupied (unoccupied) dot, the particle number in the leads increases (decreases) during the transient regime $t<\tau_{r}$. However, for either initial condition, the total energy in the leads undergoes the same increase, as shown in Fig. 5. This results in the dependence of sign on initial condition in Fig. 6. At long times and low temperature (bottom panels), $C_{\ell}^{\left(\times\right)}\left(t\right)$ increases, hinting at a buildup of correlations between charge and energy transport that is completely absent within the QME.

V.3 Heat transfer statistics and thermodynamics

We now examine the thermodynamic effects of simultaneously applying both a voltage bias and a temperature gradient. We will first focus on the flow of heat (see Eq. (29)) and its fluctuations through the system. In Fig. 7 the first cumulant of heat transfer, $C_{1,L}^{h}\left(t\right)$, is shown in the presence of a voltage bias of $V=1\Gamma$. The left lead is held at a high temperature $T_{L}=2.5\Gamma$, while $T_{R}$ is set so as to enforce a temperature ratio of either $T_{L}/T_{R}=1.0$ (upper panel) or $T_{L}/T_{R}=0.1$ (lower panel). Therefore, in the second case, a thermal bias is applied in a direction opposite to that of the voltage bias. We observe that while the left lead is heated either with or without a temperature bias, this heating is greatly enhanced when energy is driven into it from the hotter right lead. In addition, mild deviations between the iQMC and QME persist for the the gradients of the cumulants even at long times, even though both leads are at a high temperature. This occurs for the same reasons pointed out in the discussion of Fig. 3.

In Fig. 8 the second cumulants of heat transfer are shown for the same parameters as in Fig. 7. The thermal fluctuations are rather large even without a temperature bias, and their enhancement by this bias is not as strong. The QME is not able to capture these fluctuations as well as it captures the mean heat transfer, resulting in more significant deviations from the numerically exact data. These deviations are more obvious in the presence of a thermal bias, despite the larger average temperature of the system.

In our analysis so far, we have considered the entire time evolution following the coupling quench. We now investigate the steady state thermodynamics that eventually emerges. This is extracted from the long-time dynamics of the cumulants (see Eq. (30)). In Fig. 9 we present the particle current (upper panels) and the associated shot noise (lower panels), plotted on a logarithmic scale as functions of the left–right lead temperature ratio $T_{L}/T_{R}$ at constant $T_{L}$. $T_{L}$ takes on the same two values that were considered in previous plots, one of which (shown in left panels) is above the Kondo temperature at $U=8\Gamma$ and the other of which (right panels) is below it. The right lead temperature varies from $100T_{L}$ to $0.1T_{L}$. QME data is shown at interaction strength $U=8\Gamma$ (black dashed line) and $U=4\Gamma$ (green dashed line), and iQMC data at $U=8\Gamma$ (black crosses) and $U=4\Gamma$ (green crosses). Our estimate for the magnitude of numerical uncertainties Cohen et al. (2015); Antipov et al. (2017) is smaller than the symbol size.

The definition we have used for the particle current can be interpreted as the flow rate into the left lead. Its sign at the parameters considered here is always negative, i.e. particles are flowing from left to right, consistent with the directionality induced by the voltage gradient. Within the QME picture, the system is in the Coulomb blockade regime, such that reducing the interaction strength consistently increases the magnitude of the particle current and noise (cf. green and black dashed curves in Fig. 9).

As might be expected, the QME successfully reproduces the numerically exact iQMC results for both current and noise at high (left and right) temperatures. It therefore proves to be an excellent approximation in the left part of the left panels (high $T_{L}$ and $T_{R}$), and a reasonable one at the right edge of these panels (high $T_{L}$ and low $T_{R}$). In the low temperature case where $\beta_{L}=50/\Gamma$, the QME fails qualitatively. In particular, the QME predicts complete suppression of both current and noise when both leads are at low temperatures. In the physical regime explored here, as the temperature decreases, this may first be attributed to higher-order scattering processes, and later to the formation of a correlated Kondo transport channel. Nonmonotonic behavior can be observed in the iQMC results at high $T_{L}/T_{R}$. Here, the equilibrium system might be expected to be driven deeper into the Kondo regime by the lower average temperature, but this effect competes with the nonequilibrium fluxes that eventually break down the Kondo singlet.

In Fig. 10, we move on to examine the steady state values of heat current (upper panels), heat current fluctuations (middle panels) and the logarithm of the entropy production rate normalized by the left lead temperature (lower panels). In these plots the entropy is temperature-normalized so that the $T_{L}>T_{K}$ and $T_{L}<T_{K}$ plots may be compared on the same scale, and we consider the logarithm because the growth of entropy production with increasing average temperature is known to be exponential. As before, we consider temperatures $T_{L}=2.5\Gamma>T_{K}$ (left panels) and $T_{L}=0.02\Gamma<T_{K}$ (right panels).

At high temperature (left panels), QME provides a good estimate of the heat current, as it did for the particle current. QME consistently underestimates the fluctuations in heat current, but does predict the correct trend in that fluctuations decrease when the right lead is cooled down. However, it also predicts a strong and spurious dependence on $U$ at low temperature ratios, which is not observed in the iQMC results. A reversal in the direction of the steady state heat current as a function of the temperature ratio, as illustrated in Fig. 1, is observed in both methods (upper left panel). The temperature ratio at which this occurs is indistinguishable from 1 within our numerical resolution, indicating that the direction of heat flux is completely determined by that of the temperature gradient in this parameter regime.

The QME picture for the high temperature entropy generation $\sigma$ (lower panel of Fig. 10) is surprisingly accurate, though this is partially due to the logarithmic scale. A monotonic and approximately exponential growth with the temperature ratio is observed, with different exponents on the two sides of $T_{L}/T_{R}=1$.

For low $T_{L}$ (right panels of Fig. 10, QME fails qualitatively except where $T_{R}$ is large. The change in sign at $T_{L}/T_{R}=1$ observed at high temperature vanishes at low temperature. In the QME approximation, this occurs simply because the heat current vanishes entirely. The iQMC result does show a finite heat current which remains positive for $T_{R}<T_{L}$, such that the left lead continues to heat in this case. This effect and the wrong QME result can be understood by considering the two components of the heat current (see Eq. (29)) separately. At such low lead temperatures, the energy current (not shown) is essentially zero, a fact that is captured by QME; while the particle current continues to be carried by the Kondo channel, not captured by QME. The lack of sign reversal in the heat current is accompanied by an increase in its fluctuations, again in opposition to the high temperature result. This may be interpreted as a precursor to an eventual reversal that may manifest if the right lead is cooled further, or if the voltage bias is reduced.

Unlike at high temperatures, the QME estimate for the low temperature entropy production (lower right panel of Fig. 10) fails catastrophically. Here, it predicts a spurious structure with multiple peaks and troughs instead of the correct behavior, which is essentially exponential. At the lower temperature, the difference in exponents at the two sides of $T_{L}/T_{R}=1$ is either nonexistent or too small to detect within our numerical resolution. The detectable exponent is clearly larger than that appearing in the high temperature case, as is the overall entropy production rate. This result is most likely connected to the opening of the Kondo transport channel, and invites further theoretical analysis beyond the scope of the present work.

VI Conclusions

We presented a numerically exact method for evaluating full counting statistics (FCS) in nonequilibrium quantum junctions based on the inchworm quantum Monte Carlo (iQMC) approach. The method accounts for both particle and energy transport statistics, and is applicable in a wide range of parameters that includes the strongly correlated Kondo regime.

We benchmarked the method against the nonequilibrium Green’s function formalism in the noninteracting case. We also carried out an extensive comparison between the iQMC and a simple quantum master equation (QME) approach at different temperature regimes, showing clearly where the latter method fails. Surprisingly, we found that the QME approximation fails to produce the correct energy counting statistics even at temperatures significantly higher than $T_{K}$, although agreement between QME and iQMC improves as the temperature tends to infinity. We further found that the presence of a thermal gradient across the molecular junction makes the agreement worse even at high temperatures, and when the gradient increases the overall average temperature.

At temperatures above the Kondo threshold, we found that there is a left–right symmetry breaking in the system that can be observed in the first cumulant of energy transfer: energy is distributed among the leads in away that depends on the initial preparation of the system. When the temperature is lowered below $T_{K}$, the iQMC calculations predict finite values for the particle and heat current and noise, whereas the QME method predicts full suppression of both current and noise due to the Coulomb blockade effect. In general, the disagreement between the QME and iQMC approaches is most significant in the noise and the energy–particle cross-correlations, confirming that noise measurements offer more sensitive probes of higher-order scattering processes and many-body correlations than average currents.

Finally, we investigated the steady state entropy production rates at different interaction strengths using both the Monte Carlo and master equation methods. Among other things, we found that the entropy production rate from master equations spuriously predicts an opposite trend to that computed from iQMC as the average temperature of the system is reduced significantly below $T_{K}$.

This paper provides the basis for future investigations into several more specific questions, including the FCS of energy transport in periodically-driven systems and the properties of levitons in the Kondo regime, as in Ref. Suzuki, 2017. Another interesting direction is models with multiple orbitals, where QME may fail at any temperatures due to an inability to properly account for bath induced coherences in the system. In general, the tools presented here can to study a variety of fundamental questions in quantum thermodynamics, by obtaining the time-dependence of entropy production, testing thermodynamic uncertainty relations and validating fluctuation–dissipation relations in strongly correlated quantum many-body systems.

Acknowledgements.

G.C. acknowledges support by the Israel Science Foundation (Grant No. 1604/16). E.G. was supported by DOE ER 46932. M.G. acknowledges support by the National Science Foundation (Grant CHE-1565939). International exchange and collaboration was supported by Grant No. 2016087 from the United States-Israel Binational Science Foundation (BSF) and by University of California San Diego.