Quantum energy exchange and refrigeration: a full-counting statistics approach

The large-scale heat engine was instrumental in the development of classical thermodynamics in the 19th century. In order to establish the theory of thermodynamics from quantum principles, studies of the nanoscale analogue, the quantum heat engine (QHE), are ongoing [1–3]. What is 'quantum' in QHEs [4, 5]? The construction of the QHE differs it from the classical heat engine by having a quantum system, comprising a set of discrete-quantized states, as the working fluid analogue in many models of such machines [6]. Also, the degrees of freedom of the reservoirs are described by quantum statistical distributions such as the Bose–Einstein or Fermi–Dirac functions. The quantum system is either periodically driven by a classical field, as in the four-stroke Otto engine [7–9], or operated continuously, to realize e.g. an absorption refrigerator [10]. It is necessary to consider the ramifications of non-trivial quantum effects such as quantum coherences [11–14], quantum correlations [15] and reservoir engineering e.g. by squeezing [16–20], on the performance of heat engines to potentially overcome classical bounds.

A central aspect of nanoscale QHEs is that they are not restricted to the weak system–bath coupling limit. Classical-macroscopic thermodynamics is a weak-coupling theory; the effect of the boundary between the system and its environment is small relative to the bulk behavior. In contrast, small systems may strongly couple to their environment, in the sense that the interaction energy between the system and the bath is comparable to the frequencies of the isolated system. While quantum thermodynamical machines were traditionally analyzed under a strict weak-coupling assumption, it is now recognized that to properly characterize the performance of quantum engines, one must develop methods that are not limited in this respect [21–32].

What might be the impact of a strong system–bath interaction energy? Consider, for example, a heat conducting two-terminal nanojunction [33]. A weak-coupling (Born–Markov) treatment of the thermal current colossally fails beyond the strict weak-coupling regime [34]: while a weak-coupling theory predicts a linear enhancement of the thermal current with the increase of the system–bath interaction energy [35–37], a strong-coupling treatment administers a turnover behavior [35, 36, 38–40]. This crossover, between first enhancement of the current and subsequent suppression with increasing interaction strength, appears once the system–bath interaction energy becomes comparable to the system's natural frequencies.

Physically, strong system–bath interactions are responsible for three-phonon (and more) scattering contributions to thermal transport problems, beyond the weak-coupling resonant term. Alternatively, rather than focusing on the distinction between strong and weak-coupling, one may classify system–bath interaction operators based on whether they are additive or non-additive in the different reservoirs. These two classes of models, additive and non-additive, realize distinctive energy transport characteristics and refrigeration as we show in this work.

An autonomous absorption refrigerator transfers thermal energy from a cold bath to a hot bath without an input power by using thermal energy provided from a so-called work reservoir. In a recent study [41], we demonstrated that the smallest system, a qubit, is incapable of operating as a quantum absorption refrigerator (QAR) when the baths are coupled weakly-additively to the system. However, we demonstrated that by coupling the system in a non-additive manner to three heat baths, which are spectrally structured, a cooling function was achieved. Moreover, we showed that the system reached the Carnot bound when the reservoirs were characterized by a single frequency component. This study thus clearly illustrates that non-additive models can bring in a new function, which is missing in additive scenarios. On the other hand, as was pointed out in [41] and illustrated earlier on in [35, 36, 38], the dynamics of non-additive models can be treated with standard kinetic quantum master equations (QMEs), albeit with a non-additive, baths-cooperative dissipator.

The objective of this paper is to present a rigorous, thermodynamically consistent formalism for the calculation of energy exchange (current and cumulants) in additive and non-additive interaction models, thus develop the groundwork of the qubit-QAR presented in [41]. Our formalism utilizes a full-counting statistics (FCS) approach that provides cumulants of energy exchange to all orders [26, 42–45]. In this method, rather than focusing directly on the averaged heat current, the first cumulant, we analyze the probability distribution of exchanged energy with each bath, within a certain interval of time t. More specifically, we study the Fourier transform of this probability distribution function, the so-called characteristic function. This function can be derived from an equation of motion for the counting-field dependent reduced density operator, which we organize in the form of a Redfield equation. The formalism is applied to treat both additive and non-additive interaction models. Under the secular approximation, our derivation is thermodynamically consistent: we confirm that the derived cumulant generating function (CGF) satisfies the entropy production fluctuation theorem, which is a microscopic statement of the second law [42, 43]. We exemplify our work on the two-state system, culminating this study in the exploration of the cooling performance of a qubit refrigerator.

The paper is organized as follows. We motivate the study of non-additive interaction models in section 2. Section 3 contains our method development. We define the model, derive the counting-field dependent Redfield-type equation and use it to calculate the CGF. Closed-form expressions for the energy current in two-terminal setups, assuming either additive or non-additive system–bath interaction operators are also presented. In section 4, we exemplify our analytical results on a two-state system and present simulation results of the two-terminal nonequilibrium spin-boson (NESB) model. The application of our formalism to study a QAR is given in section 5. We summarize our work in section 6. Throughout the paper we work with units where ℏ = 1, k_B = 1.

In standard models of open quantum systems, under the weak-coupling approximation, the different baths dictate additive decoherence and dissipative dynamics [46, 47]. As a result, heat currents flow independently between the quantum system and each individual bath [35, 36, 48]. Let us clarify this point. The commonly used open quantum-system Hamiltonian, with a system ${\hat{H}}_{S}$ coupled to N reservoirs is

$$\begin{eqnarray}&&\hat{H}={\hat{H}}_{S}+\displaystyle \sum _{\mu =1}^{N}{\hat{H}}_{\mu }+{\hat{H}}_{{\rm{int}}},\,\ {\hat{H}}_{{\rm{int}}}=\displaystyle \sum _{\mu =1}^{N}{\hat{S}}_{\mu }\otimes {\hat{V}}_{\mu }.\end{eqnarray} \tag{ 1 }$$

Here, ${\hat{S}}_{\mu }$ is a system operator that is coupled to the μ bath's degrees of freedom and ${\hat{V}}_{\mu }$ is an operator of the μ bath. Note that there is an inherent assumption of additivity of the interaction of different reservoirs with the system. Assuming weak system–bath coupling, an equation of motion can be derived for the reduced density matrix of the system, σ, valid to second order in the interaction Hamiltonian. Under the Markovian approximation, we receive the QME (ℏ ≡ 1)

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\sigma }{{\rm{d}}{t}}=-{\rm{i}}[{\hat{H}}_{S},\sigma ]+\displaystyle \sum _{\mu =1}^{N}{D}_{\mu }[\sigma (t)],\end{eqnarray} \tag{ 2 }$$

with the dissipators ${D}_{\mu }[\sigma (t)]$ that contain information about the μth bath's effect on the system. Decoherence and relaxation dynamics is therefore additive in the different baths, i.e., the total relaxation rate of the system is the sum of relaxation rates induced by each bath. From equation (2), the rate of energy relaxation from the system can be written as

$$\begin{eqnarray}&&\langle {\dot{\hat{H}}}_{S}\rangle =\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\mathrm{Tr}[{\hat{H}}_{S}\sigma (t)]=\displaystyle \sum _{\mu =1}^{N}\mathrm{Tr}[{\hat{H}}_{S}{D}_{\mu }[\sigma (t)]].\end{eqnarray} \tag{ 3 }$$

We can readily identify the currents flowing towards the system from the different baths as $\langle {J}_{\mu }(t)\rangle \equiv \mathrm{Tr}[{\hat{H}}_{S}{D}_{\mu }[\sigma (t)]]$. This analysis is obviously limited to the additive interaction Hamiltonian—and when the dynamics can be recast into (2).

In this manuscript, we present a thermodynamically consistent formalism for the calculation of energy exchange in situations that do not necessarily follow the evolution equation (2). In particular, we consider two classes of system–bath interaction models, shown as schemes in figure 1, ${\hat{H}}_{{\rm{int}}}=\gamma {\sum }_{p}{\hat{S}}^{p}\otimes {\hat{V}}_{{\rm{B}}}^{p}$, with either ${\hat{V}}_{{\rm{B}}}^{p}={\sum }_{\mu =1}^{N}{\hat{B}}_{\mu }^{p}$ (panel a) or ${\hat{V}}_{{\rm{B}}}^{p}={{\rm{\Pi }}}_{\mu =1}^{N}{\hat{B}}_{\mu }^{p}$ (panel b). Here, ${\hat{S}}^{p}$ is a system operator, ${\hat{B}}_{\mu }^{p}$ is a μ-bath operator, and γ is introduced to keep track of the perturbative expansion. The baths may couple to multiple system's operators, counted by the index p. The first model is additive (ADD) and separable in the reservoirs' interaction operators ${\hat{B}}_{\mu }^{p}$, as in equation (1). The second model is non-additive (NADD), or multiplicative, in ${\hat{B}}_{\mu }^{p}$ allowing for a cooperative effect of the baths. In both cases, however, we can study the system's dynamics and its energy transport characteristics by using a perturbation theory treatment up to second order in γ and calculating the CGF.

Zoom In Zoom Out Reset image size

We distinguish here between additive and non-additive system–bath interaction models while one usually makes a division between weak- and strong-coupling limits—for additive Hamiltonians of the form (1). In the context of quantum transport, additive models were examined in the weak-coupling approximation, then improved by including higher order interaction effects, e.g. by exercising perturbation theory to the fourth order [49], using dressing [35, 36, 38] or mapping approaches [50, 51]. Recently developed numerically-exact simulation methods [22, 52–56] interrogate strong coupling effects, though they often have limitations in system size. Here, we analyze a general interaction model for multiple baths, and use it to study additive and non-additive interaction models. In both cases, we demonstrate how to derive the energy currents and noise flowing out of each bath, and our analytical results can provide physical insights to the underlying mechanisms for energy exchange.

Intuitively, the ADD interaction model is applicable when the reservoirs are physically separated from each other and therefore do not interact. The NADD model can show up in different scenarios. For example, in the treatment of the nonequilibrium (two-bath) spin-boson model one can perform a unitary (polaron) transformation of the Hamiltonian resulting in a product of displacement operators for the system–bath coupling, then proceed to study the dynamics under the Born–Markov approximation [35, 36, 38]. The bare system–bath interaction energy is incorporated inside the displacement operator, which is an exponential function, thus its effect is included to high orders. NADD models can also be accomplished by engineering many-body Hamiltonians based on e.g. resonant conditions and selection rules [10, 57]. In [48], the NADD Hamiltonian represents a chemical engine where reactants are destroyed and a chemical product is created, along with an excitation.

The NADD model may realize behaviors that cannot be captured by an ADD Hamiltonian. We exemplify this fact by studying a two-level system to realize the smallest QAR, with a qubit as its working fluid analogue [41]. A cooperative bath behavior is what allows for refrigeration.

We present a perturbative formalism for deriving the characteristic function for quantum energy exchange between a system and multiple attached thermal reservoirs. In steady state, we ultimately obtain the CGF for quantum energy exchange. This function provides the energy current and any higher order cumulant of the current such as the noise power.

Assuming weak system–bath coupling and memoryless reservoirs, the procedure is formulated in the language of a counting-field dependent Markovian QME—the Redfield equation—generalized to describe the FCS of energy exchange in a multi-terminal geometry. Our formalism is valid in the weak-coupling limit but allows for system transitions induced cooperatively by N baths.

3.1. Hamiltonian

Our discussion begins with a general open quantum system Hamiltonian, a system ${\hat{H}}_{S}$ coupled to N reservoirs,

$$\begin{eqnarray}&&\hat{H}={\hat{H}}_{S}+\displaystyle \sum _{\mu =1}^{N}{\hat{H}}_{\mu }+\gamma \displaystyle \sum _{p}{\hat{S}}^{p}\otimes {\hat{V}}_{{\rm{B}}}^{p}.\end{eqnarray} \tag{ 4 }$$

Here, ${\hat{S}}^{p}$ is a system operator and ${\hat{V}}_{{\rm{B}}}^{p}$ is an operator of the environment (bath) that depends on the nature and number of reservoirs that are included in the model. Note that the summation over p allows for many possibilities for system coupling to the reservoirs in different configurations with operators ${\hat{S}}^{p}$. The Hamiltonian is time-independent. Therefore, in the steady-state limit, energy exchange between the system and the reservoirs corresponds to heat transfer.

To study energy transfer over the time interval [0, t], between the system and the different reservoirs, we write down the energy current from a specific reservoir ν as a change in the ν bath energy ${J}_{\nu }(t)=-\tfrac{{\rm{d}}{\hat{H}}_{\nu }(t)}{{\rm{d}}t}$. Operators are written in the Heisenberg representation, $\hat{A}(t)={\hat{U}}^{\dagger }(t)\hat{A}\hat{U}(t)$, where $\hat{U}(t)={{\rm{e}}}^{-{\rm{i}}\hat{H}t}$. Therefore, the total energy exchange at the ν terminal is given by the integrated current

$$\begin{eqnarray}&&{Q}_{\nu }(t,{t}_{0}=0)={\int }_{0}^{t}{J}_{\nu }({t}^{{\rm{{\prime} }}}){\rm{d}}{t}^{{\rm{{\prime} }}}={\hat{H}}_{\nu }(0)-{\hat{H}}_{\nu }(t).\end{eqnarray} \tag{ 5 }$$

Employing the two-time measurement protocol [42, 44], the characteristic function for energy transfer at multiple baths takes the form

$$\begin{eqnarray}&&{ \mathcal Z }(\lambda =\{{\lambda }_{\mu }\},t)\equiv \mathrm{Tr}\left[{{\rm{e}}}^{{\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }(0)}{{\rm{e}}}^{-{\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }(t)}\rho (0)\right],\end{eqnarray} \tag{ 6 }$$

where we have introduced λ_μ as a real-valued parameter for counting energy at the μ reservoir, and ρ(0) is the total density matrix at the initial time t = 0. In the long-time limit, the CGF is defined as

$$\begin{eqnarray}&&G(\lambda )\equiv \mathop{\mathrm{lim}}\limits_{t\to \infty }\displaystyle \frac{1}{t}\mathrm{ln}{ \mathcal Z }(\lambda ).\end{eqnarray} \tag{ 7 }$$

Differentiating G(λ) with respect to (iλ) yields the steady-state energy current cumulants. We now recast equation (6) in the structure of the Liouville equation by explicitly writing down the time evolution operators, factorizing the exponentials, using the cyclic property of the trace operation and the assumption that ρ(0) commutes with ${\hat{H}}_{\mu }$. We compactly write

$$\begin{eqnarray}{ \mathcal Z }(\{{\lambda }_{\mu }\},t) & = & \mathrm{Tr}[{\hat{U}}^{-\lambda }(t)\rho (0){\hat{U}}^{\lambda \dagger }(t)]\\ & = & \mathrm{Tr}[{\rho }^{\lambda }(t)].\end{eqnarray} \tag{ 8 }$$

Here, we define the time evolution operators ${\hat{U}}^{\lambda \dagger }(t)\equiv {{\rm{e}}}^{{\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }/2}{\hat{U}}^{\dagger }(t){{\rm{e}}}^{-{\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }/2}$, and ${\hat{U}}^{-\lambda }(t)\equiv {{\rm{e}}}^{-{\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }/2}\hat{U}(t){{\rm{e}}}^{{\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }/2}$. In our convention, the sign of the subscript determines the sign of the counting terms in the exponent. We also define the counting-field dependent time evolution operator,

$$\begin{eqnarray}&&{\hat{U}}^{-\lambda }(t)\equiv {{\rm{e}}}^{-{\rm{i}}{\hat{H}}^{-\lambda }t},\end{eqnarray} \tag{ 9 }$$

with ${\hat{H}}^{-\lambda }\equiv {{\rm{e}}}^{-{\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }/2}\hat{H}{{\rm{e}}}^{{\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }/2}$. Since ${\hat{H}}_{\mu }$ commutes with all operators on the right hand side of equation (4), except ${\hat{V}}_{{\rm{B}}}^{p}$, we obtain a counting-dependent total Hamiltonian as, for example,

$$\begin{eqnarray}&&{\hat{H}}^{-\lambda }={\hat{H}}_{S}+\displaystyle \sum _{\mu =1}^{N}{\hat{H}}_{\mu }+\gamma \displaystyle \sum _{p}{\hat{S}}^{p}\otimes {\hat{V}}_{{\rm{B}}}^{p,-\lambda },\end{eqnarray} \tag{ 10 }$$

where ${\hat{V}}_{{\rm{B}}}^{p,-\lambda }={{\rm{e}}}^{-{\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }/2}\,{\hat{V}}_{{\rm{B}}}^{p}\,{{\rm{e}}}^{{\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }/2}$. In order to ultimately reach the CGF, we present in section 3.2 a QME to obtain the counting-field dependent density matrix in equation (8).

3.2. Counting-field dependent second-order QME

Continuing from equation (8), we switch to the interaction representation and define the counting-field density matrix,

$$\begin{eqnarray}&&{\rho }_{I}^{\lambda }(t)\equiv {\hat{U}}_{I}^{-\lambda }(t)\rho (0){\hat{U}}_{I}^{\lambda \,\dagger }(t),\end{eqnarray} \tag{ 11 }$$

with ${\hat{U}}_{I}^{-\lambda }(t)={\rm{T}}\,{{\rm{e}}}^{-{\rm{i}}{\int }_{0}^{t}{\hat{H}}_{{\rm{int}}}^{-\lambda }(\tau ){\rm{d}}\tau }$, ${\hat{H}}_{{\rm{int}}}^{-\lambda }(\tau )=\gamma {\sum }_{p}{\hat{S}}^{p}(\tau )\otimes {\hat{V}}_{{\rm{B}}}^{p,-\lambda }(\tau )$. Operators are now written in the interaction representation, $\hat{A}(t)={\hat{U}}_{0}^{\dagger }(t)\hat{A}{\hat{U}}_{0}(t)$ where ${\hat{U}}_{0}(t)={{\rm{e}}}^{-{\rm{i}}{\hat{H}}_{0}t}$, ${\hat{H}}_{0}={\hat{H}}_{S}+{\sum }_{\mu =1}^{N}{\hat{H}}_{\mu }$. Equation (11) can be written as a differential equation,

$$\begin{eqnarray}&&{\dot{\rho }}_{I}^{\lambda }(t)={\rm{i}}\,{\rho }_{I}^{\lambda }(t)\,{\hat{H}}_{{\rm{int}}}^{\lambda }(t)-{\rm{i}}\,{\hat{H}}_{{\rm{int}}}^{-\lambda }(t)\,{\rho }_{I}^{\lambda }(t),\end{eqnarray} \tag{ 12 }$$

which reduces to the quantum Liouville equation when λ = 0. For convenience, we omit below the interaction representation descriptor 'I' moving forward.

We follow the standard derivation of the second order Markovian master equation [46]—while taking into account the counting information. An important assumption is that the thermal baths are prepared in a canonical thermal state. For details, see appendix A. The result of this treatment is the counting-field dependent Redfield-type equation [47] for ${\sigma }^{\lambda }(t)\equiv {\mathrm{Tr}}_{{\rm{B}}}[{\rho }^{\lambda }(t)]$, satisfying (Schrödinger representation),

$$\begin{eqnarray}{\dot{\sigma }}_{{nm}}^{\lambda }(t) & = & -{\rm{i}}{[{\hat{H}}_{S},\sigma (t)]}_{{nm}}+\displaystyle \sum _{p,q,j,k}[-{\sigma }_{{km}}^{\lambda }(t){R}_{{\rm{B}};{nj},{jk}}^{p,{q}^{(+)}}({E}_{k,j})-{\sigma }_{{nj}}^{\lambda }(t){R}_{{\rm{B}};{jk},{km}}^{p,{q}^{(-)}}({E}_{j,k})\\ & & +{\sigma }_{{jk}}^{\lambda }(t)({R}_{{\rm{B}};{km},{nj}}^{q,\lambda ,p,-{\lambda }^{(-)}}({E}_{k,m})+{R}_{{\rm{B}};{km},{nj}}^{q,\lambda ,p,-{\lambda }^{(+)}}({E}_{j,n}))].\end{eqnarray} \tag{ 13 }$$

Here, ${E}_{m,n}={E}_{m}-{E}_{n}$ and E_n are the energy eigenstates of the system Hamiltonian ${\hat{H}}_{S}$. p and q sum over the different operators that couple to the system. The other indices, j, k, count eigenstates of the system. The different terms in this equation are defined in appendix A. The characteristic function is obtained by tracing ${\sigma }^{\lambda }(t)$ over the system states.

Equation (13) can be further simplified by performing the secular approximation, so as to decouple population and coherence dynamics. The population dynamics, ${p}_{n}(t)\equiv {\sigma }_{{nn}}(t)$, then follows

$$\begin{eqnarray}&&{\dot{p}}_{n}^{\lambda }(t)=-{p}_{n}^{\lambda }(t)\displaystyle \sum _{p,q,j}{M}_{{\rm{B}};{nj},{jn}}^{p;q}({E}_{n,j})+\displaystyle \sum _{p,q,j}{p}_{j}^{\lambda }(t)\,{M}_{{\rm{B}};{nj},{jn}}^{q,\lambda ;p,-\lambda }({E}_{j,n}),\end{eqnarray} \tag{ 14 }$$

with

$$\begin{eqnarray}&&{M}_{{\rm{B}};{ab},{cd}}^{p,\lambda ;{p}^{{\prime} },{\lambda }^{{\prime} }}(s)\equiv {\gamma }^{2}\,{S}_{{ab}}^{p}{S}_{{cd}}^{{p}^{{\prime} }}\,\langle {\hat{V}}_{{\rm{B}}}^{p,\lambda }(s){\hat{V}}_{{\rm{B}}}^{{p}^{{\prime} },{\lambda }^{{\prime} }}(0)\rangle .\end{eqnarray} \tag{ 15 }$$

Here, the subscript B reminds us that the correlation function involves multiple thermal baths, a, b, c, and d are eigenstates of the system, p and p' index the different system–bath coupling terms. The averages are performed with respect to the initial condition. The rate constants are Fourier transforms of the correlation functions, ${M}_{{\rm{B}};{ab},{cd}}^{p,\lambda ;{p}^{{\prime} },{\lambda }^{{\prime} }}(\omega )\equiv {\int }_{-\infty }^{\infty }{{\rm{e}}}^{{\rm{i}}\omega s}\,{M}_{{\rm{B}};{ab},{cd}}^{p,\lambda ;{p}^{{\prime} },{\lambda }^{{\prime} }}(s)\,{\rm{d}}{s}$. Equation (14) is a Markovian-secular counting-field dependent QME for the system population dynamics, with the system coupled to N baths. We can write it down in a compact matrix form, with the counting-dependent Liouvillian ${{ \mathcal L }}^{\lambda }$, as

$$\begin{eqnarray}&&| {\dot{p}}^{\lambda }(t)\rangle ={{ \mathcal L }}^{\lambda }| {p}^{\lambda }(t)\rangle .\end{eqnarray} \tag{ 16 }$$

Analogous derivations of the counting-field dependent reduced density operator, for charge transfer problems, were performed in e.g. [58–62]. Recall that equation (16) gives us the characteristic function ${ \mathcal Z }$ from equation (8) and ultimately the CGF via equation (7).

While our derivation is standard, our objective here has been to highlight that equations (13) and (14) are valid for both ADD and NADD interaction Hamiltonians. These two equations are central to this paper and can be used to calculate energy transfer cumulants for a general open quantum system setup—given by the Hamiltonian in equation (4).

3.3. Energy current

We organize a closed expression for the energy current from equation (16) by differentiating G with respect to the counting-field λ. The characteristic function is given by a trace over the system states. From equation (8), ${ \mathcal Z }(\lambda ,t)=\langle 1| {p}^{\lambda }(t)\rangle$ with $\langle 1| ={(1,1,\ldots 1)}^{T}$ as the identity vector. As mentioned before, differentiating G(λ) with respect to (iλ_ν) returns the steady-state current between the system and the ν bath, $\langle {J}_{\nu } \rangle ={{\rm{lim}}}_{t\to \infty }\tfrac{1}{t}\tfrac{{\rm{\partial }}{\rm{ln}}Z(\lambda ,t)}{{\rm{\partial }}({\rm{i}}{\lambda }_{\nu })}{| }_{\lambda =0}$,

$$\begin{eqnarray} &&\langle {J}_{\nu } \rangle = \langle 1| \displaystyle \frac{{\rm{\partial }}{L}^{\lambda }}{{\rm{\partial }}({\rm{i}}{\lambda }_{\nu })}{\left.{| }_{\lambda =0}| {p}^{{ss}} \rangle =\displaystyle \sum _{n,j}{\left(\displaystyle \frac{{\rm{\partial }}{L}^{\lambda }}{{\rm{\partial }}({\rm{i}}{\lambda }_{\nu })}\right|}_{\lambda =0}\right)}_{{nj}}{p}_{j}^{{ss}},\end{eqnarray} \tag{ 17 }$$

with the steady state populations p_j^ss, which we get by solving equation (16) in the long-time limit. The FCS approach thus provides a rigorous working expression of the current for ADD and NADD models on the same footing. Higher order cumulants can also be calculated by taking higher order derivatives of G, such as the second cumulant, the noise power, calculated in section 4.2. Beyond the Markov approximation, it is useful to note that non-Markovian effects do not enter the steady state current, but they do lead to corrections to higher order cumulants. Such non-Markovian effects can be evaluated with techniques developed for the study of full-counting statistics in charge transport [63].

So far, we considered a system coupled to N thermal baths, with N counting parameters, λ₁, λ₂, ..., λ_ν, ..., λ_N. For simplicity, in what follows we count energy at a single-bath only, the ν bath, and denote λ = λ_ν. We study two different types of system–bath couplings; additive ${\hat{V}}_{{\rm{B}}}={\sum }_{\mu =1}^{N}{\hat{B}}_{\mu }$ and non-additive ${\hat{V}}_{{\rm{B}}}={\prod }_{\mu =1}^{N}{\hat{B}}_{\mu }$. The latter case may be realized by building up a compound interaction operator, e.g., through a unitary transformation of operators. Further, for simplicity, we treat a single system–bath operator, i.e., we ignore the p, q summation in equation (14). Finally, we henceforth ignore the γ² factor.

The counting-dependent Liouvillian ${{ \mathcal L }}^{\lambda }$ is made of correlation functions for the ν bath (15),

$$\begin{eqnarray}&&{M}_{\nu ;{nj},{jn}}^{\lambda ;-\lambda }(s)= \langle {{\rm{e}}}^{{\rm{i}}{\hat{H}}_{\nu }\lambda }{{\rm{e}}}^{{\rm{i}}{\hat{H}}_{\nu }s}{\hat{B}}_{\nu }{{\rm{e}}}^{-{\rm{i}}{\hat{H}}_{\nu }s}{{\rm{e}}}^{-{\rm{i}}{\hat{H}}_{\nu }\lambda }{\hat{B}}_{\nu } \rangle {S}_{{nj}}{S}_{{jn}}.\end{eqnarray} \tag{ 18 }$$

Its Fourier transform in frequency domain gives the relation,

$$\begin{eqnarray}&&{M}_{\nu ;{nj},{jn}}^{\lambda ;-\lambda }(\omega )={{\rm{e}}}^{-{\rm{i}}\omega \lambda }{M}_{\nu ;{nj},{jn}}(\omega ),\end{eqnarray} \tag{ 19 }$$

which is useful for calculating derivatives in equation (17). It is also useful to note the detailed balance relation with counting parameters for a specific bath ν at inverse temperature β_ν,

$$\begin{eqnarray}&&{M}_{\nu ;{nj},{jn}}^{\lambda ;-\lambda }(\omega )={{\rm{e}}}^{({\beta }_{\nu }-{\rm{i}}\lambda )\omega }{M}_{\nu ;{nj},{jn}}(-\omega ).\end{eqnarray} \tag{ 20 }$$

For details, see appendix B.

3.3.1. Additive bath interaction model: ${\hat{V}}_{{\rm{B}}}={\sum }_{\mu =1}^{N}{\hat{B}}_{\mu }$

In the ADD case, ${\hat{H}}_{{\rm{int}}}^{\lambda }=\gamma \,\hat{S}\otimes {\hat{V}}_{{\rm{B}}}^{\lambda }$ with ${\hat{V}}_{{\rm{B}}}^{\lambda }(t)={\hat{B}}_{1}(t)+{\hat{B}}_{2}(t)\,+...+\,{\hat{B}}_{\nu }^{\lambda }(t)\,+...+\,{\hat{B}}_{N}(t)$. Recall that baths' operators are written in the interaction representation, the baths are initially uncorrelated and we assume that $\langle {\hat{V}}_{{\rm{B}}}\rangle =0$. This results in the rate constants in equation (14) being additive in the N reservoirs,

$$\begin{eqnarray}&&{M}_{{\rm{B}};{nj},{jn}}^{\lambda ;-\lambda }({E}_{j,n})={M}_{1;{nj},{jn}}({E}_{j,n})\,+...+\,{M}_{\nu ;{nj},{jn}}^{\lambda ;-\lambda }({E}_{j,n})\,+...+\,{M}_{N;{nj},{jn}}({E}_{j,n}).\end{eqnarray} \tag{ 21 }$$

The Liouvillian is then given by adding up all contributions,

$$\begin{eqnarray}&&| {\dot{p}}^{\lambda }(t) \rangle =({L}_{1}\,+...+\,{L}_{\nu }^{\lambda }\,+...+\,{L}_{N})| {p}^{\lambda }(t) \rangle .\end{eqnarray} \tag{ 22 }$$

Using equations (14), (17) and (19) we receive for the energy current

$$\begin{eqnarray} \langle {J}_{\nu } \rangle & = & \displaystyle \sum _{n,j}{p}_{j}^{{ss}}{E}_{n,j}| {S}_{{nj}}{| }^{2}{\int }_{-\infty }^{\infty }{\rm{d}}s\,{{\rm{e}}}^{{{\rm{i}}E}_{j,n}s} \langle {\hat{B}}_{\nu }(s){\hat{B}}_{\nu }(0) \rangle \\ & = & {\rm{Tr}}[{\hat{{\rm{H}}}}_{{\rm{S}}}{L}_{\nu }\,{{\rm{p}}}^{{\rm{ss}}}],\end{eqnarray} \tag{ 23 }$$

with ${L}_{\nu }\,{p}^{{ss}}$ as the dissipator of the dynamics due to the ν reservoir. It is worth commenting that this relation can be immediately derived from the energy balance equation for the system energy operator, $\tfrac{{\rm{d}}\mathrm{Tr}[{\hat{H}}_{S}\sigma (t)]}{{\rm{d}}{t}}={\sum }_{\mu =1}^{N}\langle {J}_{\mu }\rangle$, as discussed in section 2. Nevertheless, the formalism presented here can be used to feasibly provide higher order cumulants.

3.3.2. Non-additive bath interaction: ${\hat{V}}_{{\rm{B}}}={\prod }_{\mu =1}^{N}{\hat{B}}_{\mu }$

We now assume that the system is coupled to the environment according to the following form, ${\hat{V}}_{{\rm{B}}}^{\lambda }(t)={\hat{B}}_{1}(t)\otimes {\hat{B}}_{2}\,\otimes ...\otimes \,{\hat{B}}_{\nu }^{\lambda }(t)\,\otimes ...\otimes \,{\hat{B}}_{N}(t)$. This structure arises e.g. in the study of the spin-boson model after a polaron transformation [38]. We insert this product form into the definition of the correlation function and arrive at

$$\begin{eqnarray}&&{M}_{{\rm{B}};{ab},{cd}}^{\lambda ;{\lambda }^{{\prime} }}(s)= \langle {\hat{B}}_{1}(s){\hat{B}}_{1}(0) \rangle \langle {\hat{B}}_{2}(s){\hat{B}}_{2}(0) \rangle ... \langle {\hat{B}}_{\nu }^{\lambda }(s){\hat{B}}_{\nu }^{{\lambda }^{{\prime} }}(0) \rangle ... \langle {\hat{B}}_{N}(s){\hat{B}}_{N}(0) \rangle \,{S}_{{ab}}{S}_{{cd}},\end{eqnarray} \tag{ 24 }$$

since the bath initial state is factorized, ${\rho }_{{\rm{B}}}={\prod }_{\mu =1}^{N}{\rho }_{\mu }$. In frequency domain, the correlation function turns into a convolution,

$$\begin{eqnarray}{M}_{{\rm{B}};{nj},{jn}}^{\lambda ;-\lambda }({E}_{j,n}) & = & \displaystyle \frac{1}{{(2\pi )}^{N-1}}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{2}\,\ldots \,{\rm{d}}{\omega }_{\nu }\,\ldots \,{\rm{d}}{\omega }_{N}[{M}_{1;{nj};{jn}}({E}_{j,n}-{\omega }_{2}\,-...-\,{\omega }_{\nu }\,-...-{\omega }_{N})\,{M}_{2;{nj},{jn}}({\omega }_{2})\\ & & ...\,{M}_{\nu ;{nj},{jn}}^{\lambda ;-\lambda }({\omega }_{\nu })\ldots \,{M}_{N;{nj},{jn}}({\omega }_{N})].\end{eqnarray} \tag{ 25 }$$

This rate constant describes a cooperative effect: the system makes a transition between the n and j states, and energy is absorbed or released simultaneously-partially to the N baths. From equation (14) and the relation in (19), we receive the energy current from equation (17),

$$\begin{eqnarray} \langle {J}_{\nu } \rangle & = & -\displaystyle \frac{1}{{(2\pi )}^{N-1}}\displaystyle \sum _{n,j}{p}_{j}^{{ss}}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{2}\,\ldots \,{\rm{d}}{\omega }_{\nu }\,\ldots \,{\rm{d}}{\omega }_{N}\,({\omega }_{\nu })[{M}_{1;{nj};{jn}}({E}_{j,n}-{\omega }_{2}-...-{\omega }_{\nu }-...-{\omega }_{N})\,\\ & & {M}_{2;{nj},{jn}}({\omega }_{2})...{M}_{\nu ;{nj},{jn}}({\omega }_{\nu })\ldots \,{M}_{N;{nj},{jn}}({\omega }_{N})].\end{eqnarray} \tag{ 26 }$$

We emphasize that this result holds without specifying the system, bath or the system–bath interaction Hamiltonian, aside from it being a product form.

4.1. Cumulant generating function

Spin-bath models are central to the theory of open quantum systems [64, 65]. In particular, the spin-boson model has found immense applications in condensed phase physics, chemical dynamics, quantum optics and quantum technologies [64]. It serves to develop approximation schemes perturbative and non-perturbative in the system–bath coupling strength, see for example recent studies [66, 67]. Beyond questions over quantum decoherence, dissipation, and thermalization, which can be addressed within the single-bath spin-boson model, the two-bath, NESB model has been put forward as a minimal model for exploring the fundamentals of thermal energy transfer in anharmonic nano-junctions [35]. When the two reservoirs are maintained at different temperatures away from linear response, nonlinear functionality such as the diode effect can develop in the junction [35–38]. From the theoretical perspective, the NESB model is a rich platform for developing methodologies in nonequilibrium quantum dynamics. For recent comprehensive studies, see e.g. [34, 68, 69].

In this section, we consider a nonequilibrium spin-bath model: a qubit coupled via a ${\hat{\sigma }}_{x}$ operator to N reservoirs. Our objective is to work with the formalism of section 3 and derive expressions for the current and noise in ADD and NADD spin-bath models. For simplicity, we continue with a single counting parameter λ = λ_ν, which counts energy at the ν contact. The system Hamiltonian is written in its energy eigenbasis, and the interaction operator with the bath is off-diagonal. We do not specify the reservoirs and their interaction with the qubit,

$$\begin{eqnarray}{\hat{H}}_{S} & = & \displaystyle \frac{{\omega }_{0}}{2}{\hat{\sigma }}_{z}=\displaystyle \frac{{\omega }_{0}}{2}(| 1\rangle \langle 1| -| 0\rangle \langle 0| ),\\ {\hat{H}}_{{\rm{int}}} & = & \hat{S}\otimes {\hat{V}}_{{\rm{B}}},\,\hat{S}={\hat{\sigma }}_{x}=| 0\rangle \langle 1| +| 1\rangle \langle 0| .\end{eqnarray} \tag{ 27 }$$

In the language of the general discussion, S₀₀ = S₁₁ = 0, S₀₁ = S₁₀ = 1, and E_1,0 = ω₀, E_0,1 = −ω₀. The QME (14) thus reduces to

$$\begin{eqnarray*}&&\left(\begin{array}{c}{\dot{p}}_{0}^{\lambda }(t)\\ {\dot{p}}_{1}^{\lambda }(t)\end{array}\right)=\left(\begin{array}{cc}-{M}_{{\rm{B}}}(-{\omega }_{0}) & {M}_{{\rm{B}}}^{\lambda ;-\lambda }({\omega }_{0})\\ {M}_{{\rm{B}}}^{\lambda ;-\lambda }(-{\omega }_{0}) & -{M}_{{\rm{B}}}({\omega }_{0})\end{array}\right)\left(\begin{array}{c}{p}_{0}^{\lambda }(t)\\ {p}_{1}^{\lambda }(t)\end{array}\right).\end{eqnarray*}$$

Note that we have dropped the subscripts ab, cd on M since S₀₁S₁₀ = 1. It is convenient to identify the following four rates,

$$\begin{eqnarray}&&\left(\begin{array}{c}{\dot{p}}_{0}^{\lambda }(t)\\ {\dot{p}}_{1}^{\lambda }(t)\end{array}\right)=\left(\begin{array}{cc}-{k}_{0\to 1} & {k}_{1\to 0}^{\lambda }\\ {k}_{0\to 1}^{\lambda } & -{k}_{1\to 0}\end{array}\right)\left(\begin{array}{c}{p}_{0}^{\lambda }(t)\\ {p}_{1}^{\lambda }(t)\end{array}\right).\end{eqnarray} \tag{ 28 }$$

The two eigenvalues of this matrix are

$$\begin{eqnarray*}&&{\mu }_{\pm }^{\lambda }=-\displaystyle \frac{1}{2}({k}_{0\to 1}+{k}_{1\to 0})\pm \displaystyle \frac{1}{2}\sqrt{{({k}_{0\to 1}+{k}_{1\to 0})}^{2}-4({k}_{0\to 1}{k}_{1\to 0}-{k}_{0\to 1}^{\lambda }{k}_{1\to 0}^{\lambda })}.\end{eqnarray*}$$

The rate matrix in equation (28) is the Liouvillian ${ \mathcal L }$ and is therefore used to calculate the CGF, energy current and noise power.

Recall that the characteristic function ${ \mathcal Z }(\lambda ,t)$ is the trace over the counting-field dependent reduced density matrix, and that it can be used to obtain the CGF,

$$\begin{eqnarray}&&G(\lambda )=\mathop{\mathrm{lim}}\limits_{t\to \infty }\displaystyle \frac{1}{t}\mathrm{ln}\langle I| {p}^{\lambda }(t)\rangle ,\end{eqnarray} \tag{ 29 }$$

where $\langle 1| \,=\,{(1,1)}^{T}$ is the identity vector. In the long-time limit, only the smallest-magnitude eigenvalue of the Liouvillian survives, and the CGF is given by

$$\begin{eqnarray}&&G(\lambda )={\mu }_{+}^{\lambda }.\end{eqnarray} \tag{ 30 }$$

Note that ${\mu }_{+}^{\lambda =0}=0$, which is consistent with the normalization condition of the probabilities. We show simulations of G for both ADD and NADD models in figures 2–3 and 4–5, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Let us now consider the expansion

$$\begin{eqnarray}&&G(\lambda )={a}_{1}({\rm{i}}\lambda )+{a}_{2}{({\rm{i}}\lambda )}^{2}+{a}_{3}{({\rm{i}}\lambda )}^{3}+...\end{eqnarray} \tag{ 31 }$$

The imaginary part of G(λ) is an odd function in λ. The slope of Im G(λ) around λ = 0 corresponds to the averaged energy current. Similarly, the real part of G(λ) is an even function in λ. The coefficients correspond to the cumulants of energy exchange: energy current (a₁), noise power (2a₂), skewness, etc.

4.2. Energy current and noise power

The energy currents and the corresponding zero-frequency noise powers can be readily obtained by taking derivatives of the CGF with respect to the counting-field. The first cumulant (current) and second cumulant (noise power) at the ν bath are given, respectively, by

$$\begin{eqnarray} \langle {J}_{\nu } \rangle & \equiv & \mathop{{\rm{lim}}}\limits_{t\to \infty }\displaystyle \frac{ \langle {Q}_{\nu } \rangle }{t}={\left.\displaystyle \frac{{\rm{\partial }}G(\lambda )}{{\rm{\partial }}({\rm{i}}\lambda )}\right|}_{\lambda =0}{}_{},\\ \langle {S}_{\nu } \rangle & \equiv & \mathop{{\rm{lim}}}\limits_{t\to \infty }\displaystyle \frac{ \langle {Q}_{\nu }^{2} \rangle - \langle {Q}_{\nu }{ \rangle }^{2}}{t}={\left.\displaystyle \frac{{{\rm{\partial }}}^{2}G(\lambda )}{{\rm{\partial }}{({\rm{i}}\lambda )}^{2}}\right|}_{\lambda =0},\end{eqnarray} \tag{ 32 }$$

In the long-time limit, we solve equation (28) for λ = 0 and find the steady state population ${p}_{0}^{{ss}}=\tfrac{{k}_{1\to 0}}{{k}_{1\to 0}+{k}_{0\to 1}}$, with p₀ + p₁ = 1. We then reach the following expressions,

$$\begin{eqnarray} \langle {J}_{\nu } \rangle & = & {p}_{0}^{{ss}}{\left.\displaystyle \frac{{\rm{\partial }}{k}_{0\to 1}^{\lambda }}{{\rm{\partial }}({\rm{i}}\lambda )}\right|}_{\lambda =0}{\left.{}_{}+{p}_{1}^{{ss}}\displaystyle \frac{{\rm{\partial }}{k}_{1\to 0}^{\lambda }}{{\rm{\partial }}({\rm{i}}\lambda )}\right|}_{\lambda =0}\\ \langle {S}_{\nu } \rangle & = & {p}_{0}^{{ss}}{\left.\displaystyle \frac{{{\rm{\partial }}}^{2}{k}_{0\to 1}^{\lambda }}{{\rm{\partial }}{({\rm{i}}\lambda )}^{2}}\right|}_{\lambda =0}{\left.{}_{}+{p}_{1}^{{ss}}\displaystyle \frac{{{\rm{\partial }}}^{2}{k}_{1\to 0}^{\lambda }}{{\rm{\partial }}{({\rm{i}}\lambda )}^{2}}\right|}_{\lambda =0}+\displaystyle \frac{2}{{k}_{1\to 0}+{k}_{0\to 1}}\left[{\left.\left.\left(\displaystyle \frac{{\rm{\partial }}{k}_{0\to 1}^{\lambda }}{{\rm{\partial }}({\rm{i}}\lambda )}\displaystyle \frac{{\rm{\partial }}{k}_{1\to 0}^{\lambda }}{{\rm{\partial }}({\rm{i}}\lambda )}\right.\right)\right|}_{\lambda =0}- \langle {J}_{\nu }{ \rangle }^{2}\right].\end{eqnarray} \tag{ 33 }$$

It is significant to note that these results are general, valid for arbitrary bath coupling operator ${\hat{V}}_{{\rm{B}}}$. In what follows we discuss the so-called 'thermodynamic uncertainty relation', then calculate the current and noise in the ADD and NADD models.

4.2.1. Thermodynamics uncertainty relation

In recent studies, an interesting, thermodynamic uncertainty relation was discovered for Markov processes in steady state. It connects the steady state averaged current, its fluctuations, and the entropy production rate in the nonequilibrium process Σ_Q as [70–72],

$$\begin{eqnarray}&&\displaystyle \frac{\langle S\rangle }{\langle J{\rangle }^{2}}\displaystyle \frac{{{\rm{\Sigma }}}_{Q}}{{k}_{B}}\geqslant 2.\end{eqnarray} \tag{ 34 }$$

This relation points to a crucial trade off between precision and dissipation: a precise process with little noise requires high thermodynamic-entropic cost. For a two-terminal junction with T_R > T_L, the entropy production rate is directly related to the current, ${{\rm{\Sigma }}}_{Q}=\langle J\rangle {\rm{\Delta }}\beta$ with ${\rm{\Delta }}\beta =\left(\tfrac{1}{{T}_{L}}-\tfrac{1}{{T}_{R}}\right)$. The thermodynamic uncertainty relation then reduces to

$$\begin{eqnarray}&&\displaystyle \frac{\langle {S}_{R}\rangle }{\langle {J}_{R}\rangle }{\rm{\Delta }}\beta \geqslant 2{k}_{B}.\end{eqnarray} \tag{ 35 }$$

In linear response, $\langle {J}_{R}\rangle =\kappa {\rm{\Delta }}T$, ${\rm{\Delta }}\beta ={\rm{\Delta }}T/{T}_{a}^{2}$, with $2{T}_{a}={T}_{L}+{T}_{R}$ twice the averaged temperature. The inequality then reaches the Green–Kubo relation, $\langle {S}_{R}\rangle =2{k}_{B}{T}_{a}^{2}\kappa$, with κ as the thermal conductance.

In sections 4.2.2 and 4.2.3 we show that the relation holds in additive, non-additive, symmetric, asymmetric, classical and quantum regimes. This follows from the fact that we use assume Markovianity of the dynamics throughout these different cases. Beyond that, the inequality can be violated depending on the behavior of high order transport coefficients, as we recently proved in [73].

4.2.2. Additive bath interaction model: ${\hat{V}}_{{\rm{B}}}={\hat{B}}_{L}+{\hat{B}}_{R}$

In the case of an additive system–bath coupling with two reservoirs, ${\hat{V}}_{{\rm{B}}}={\hat{B}}_{L}+{\hat{B}}_{R}$, it can be shown that the rate constants are given by (we count energy at the R terminal)

$$\begin{eqnarray*}{k}_{0\to 1}^{\lambda } & = & {\int }_{-\infty }^{\infty }{\rm{d}}{s}\,{{\rm{e}}}^{-{\rm{i}}{\omega }_{0}s}[\langle {\hat{B}}_{L}(s){\hat{B}}_{L}(0)\rangle +\langle {\hat{B}}_{R}^{\lambda }(s){\hat{B}}_{R}^{-\lambda }(0)\rangle ]\\ & \equiv & {k}_{0\to 1}^{L}+{k}_{0\to 1}^{R,\lambda }.\end{eqnarray*}$$

Explicitly, one can readily prove that the counting-field rate constant relates to the λ = 0 expression as ${k}_{0\to 1}^{R,\lambda }={{\rm{e}}}^{{\rm{i}}{\omega }_{0}\lambda }{k}_{0\to 1}^{R}$, shown previously in equation (19). Similarly, ${k}_{1\to 0}^{R,\lambda }={{\rm{e}}}^{-{\rm{i}}{\omega }_{0}\lambda }{k}_{1\to 0}^{R}$. Using these two relations in (33), the energy current and noise power reduce to

$$\begin{eqnarray}\langle {J}_{R}\rangle & = & {\omega }_{0}\,[{p}_{0}^{{ss}}\,{k}_{0\to 1}^{R}-{p}_{1}^{{ss}}\,{k}_{1\to 0}^{R}],\\ \langle {S}_{R}\rangle & = & {\omega }_{0}^{2}({p}_{0}^{{ss}}\,{k}_{0\to 1}^{R}+{p}_{1}^{{ss}}\,{k}_{1\to 0}^{R})-\displaystyle \frac{2}{{k}_{1\to 0}+{k}_{0\to 1}}(\langle {J}_{R}{\rangle }^{2}+{\omega }_{0}^{2}{k}_{0\to 1}^{R}{k}_{1\to 0}^{R}).\end{eqnarray} \tag{ 36 }$$

The expression for the current agrees with previous studies in which the energy current was defined 'by hand' in the weak-coupling limit (see the discussion in section 2 below equation (3)) [35, 36] by constructing an energy current operator [37], or by performing a classical counting statistics analysis (by resolving the Markovian master equation) [38, 74]. The noise expression agrees with the zero-frequency limit of [75].

Let us consider the NESB model, The system is given by a spin, and the thermal baths are modeled by a collection of noninteracting bosons,

$$\begin{eqnarray}&&\hat{H}=\displaystyle \frac{{\rm{\Delta }}}{2}{\hat{\sigma }}_{z}+\displaystyle \frac{{\omega }_{0}}{2}{\hat{\sigma }}_{x}+{\hat{\sigma }}_{z}\otimes {\hat{V}}_{{\rm{B}}}+\displaystyle \sum _{\mu ,j}{\omega }_{\mu ,j}{\hat{b}}_{\mu ,j}^{\dagger }{\hat{b}}_{\mu ,j}.\end{eqnarray} \tag{ 37 }$$

Here, ${\hat{\sigma }}_{x}$ and ${\hat{\sigma }}_{z}$ are the Pauli matrices as in equation (27), Δ is the energy gap between the spin levels, and ω₀ is the tunnelling energy. The two reservoirs, L and R, include uncoupled harmonic oscillators, with ${\hat{b}}_{\mu ,j}^{\dagger }$ (${\hat{b}}_{\mu ,j}$) as the bosonic creation (annihilation) operator of the jth mode in the μ = L, R reservoir. We solve the transport behavior of the model in the separable bath interaction ${\hat{V}}_{{\rm{B}}}$ limit, ${\hat{V}}_{{\rm{B}}}={\sum }_{\mu ,j}{\gamma }_{\mu ,j}({\hat{b}}_{\mu ,j}^{\dagger }+{\hat{b}}_{\mu ,j})$, ${\gamma }_{\mu ,j}$ is the system–bath interaction energy. For simplicity, let us take Δ = 0 in equation (37). We then perform a rotation, ${\hat{W}}^{\dagger }{\hat{\sigma }}_{z}\hat{W}={\hat{\sigma }}_{x}$, ${\hat{W}}^{\dagger }{\hat{\sigma }}_{x}\hat{W}={\hat{\sigma }}_{z},$ with $\hat{W}=\tfrac{1}{\sqrt{2}}({\hat{\sigma }}_{x}+{\hat{\sigma }}_{z})$, and receive the transformed Hamiltonian ${\hat{H}}_{{\rm{ADD}}}={\hat{W}}^{\dagger }\hat{H}\hat{W}$,

$$\begin{eqnarray}&&{\hat{H}}_{\mathrm{ADD}}=\displaystyle \frac{{\omega }_{0}}{2}{\hat{\sigma }}_{z}+{\hat{\sigma }}_{x}\displaystyle \sum _{\mu ,j}{\gamma }_{\mu ,j}({\hat{b}}_{\mu ,j}^{\dagger }+{\hat{b}}_{\mu ,j})+\displaystyle \sum _{\mu ,j}{\omega }_{\mu ,j}{\hat{b}}_{\mu ,j}^{\dagger }{\hat{b}}_{\mu ,j}.\end{eqnarray} \tag{ 38 }$$

This Hamiltonian offers a convenient starting point for a perturbative treatment. Using the results of section 4.1, one can readily write down the CGF of the model (30) with ${k}_{1\to 0}^{\mu }={{\rm{\Gamma }}}_{\mu }({\omega }_{0})(1+{n}_{\mu }({\omega }_{0}))$, and ${k}_{0\to 1}^{\mu }={{\rm{\Gamma }}}_{\mu }({\omega }_{0}){n}_{\mu }({\omega }_{0})$, the results for which are shown in figure 2. Here, ${n}_{\mu }(\omega )={({{\rm{e}}}^{{\beta }_{\mu }\omega }-1)}^{-1}$ is the Bose–Einstein distribution function, ${{\rm{\Gamma }}}_{\mu }(\omega )=2\pi {\sum }_{j\in \mu }{\gamma }_{\mu ,j}^{2}\delta (\omega -{\omega }_{\mu ,j})$ is the spectral density function. For details, see [38]. The energy current and noise power at the R terminal are

$$\begin{eqnarray*}&&\langle {J}_{R}\rangle =\displaystyle \frac{{\omega }_{0}\,{{\rm{\Gamma }}}_{L}({\omega }_{0}){{\rm{\Gamma }}}_{R}({\omega }_{0})[{n}_{R}({\omega }_{0})-{n}_{L}({\omega }_{0})]}{{{\rm{\Gamma }}}_{L}({\omega }_{0})[1+2{n}_{L}({\omega }_{0})]+{{\rm{\Gamma }}}_{R}({\omega }_{0})[1+2{n}_{R}({\omega }_{0})]},\end{eqnarray*}$$

$$\begin{eqnarray}&&\langle {S}_{R}\rangle =\displaystyle \frac{{\omega }_{0}^{2}\,{{\rm{\Gamma }}}_{L}({\omega }_{0}){{\rm{\Gamma }}}_{R}({\omega }_{0})[{n}_{R}({\omega }_{0})(1+{n}_{L}({\omega }_{0}))+{n}_{L}({\omega }_{0})(1+{n}_{R}({\omega }_{0}))]-2\langle {J}_{R}{\rangle }^{2}}{{{\rm{\Gamma }}}_{L}({\omega }_{0})[1+2{n}_{L}({\omega }_{0})]+{{\rm{\Gamma }}}_{R}({\omega }_{0})[1+2{n}_{R}({\omega }_{0})]}.\end{eqnarray} \tag{ 39 }$$

Other approaches for treating this model in a perturbative manner are formulated in e.g. [49, 68, 76].

We simulate both symmetric and asymmetric junctions in figure 2, and demonstrate the thermal rectification effect, an asymmetry of the current and noise with respect to ±ΔT. In figure 3 we present the current and its noise for the weak-coupling additive model (39). We further demonstrate that the thermodynamic uncertainty relation, equation (35), is satisfied in both the classical T_a ≫ ω₀ and quantum T_a ≪ ω₀ regimes. It is significant to note that in the high temperature limit the noise decreases when increasing ΔT, yet the uncertainty relation holds. Throughout our two-terminal simulations we use an ohmic spectral density function ${{\rm{\Gamma }}}_{\mu }(\omega )={\alpha }_{\mu }\omega {{\rm{e}}}^{-| \omega | /{\omega }_{\mu }}$ with α_μ as a dimensionless coefficient and ω_μ as the cutoff frequency.

4.2.3. Non-additive bath interaction: ${\hat{V}}_{{\rm{B}}}={\hat{B}}_{L}\otimes {\hat{B}}_{R}$

We now discuss the NADD model ${\hat{V}}_{{\rm{B}}}={\hat{B}}_{L}\otimes {\hat{B}}_{R}$. In this case, the Liouvillian cannot be separated to left and right-reservoir assisted processes. For a two-level system, the rates, equation (25), simplify to

$$\begin{eqnarray}{k}_{1\to 0}^{\lambda } & = & {\int }_{-\infty }^{\infty }{\rm{d}}{s}\,{{\rm{e}}}^{{\rm{i}}{\omega }_{0}s}\langle {\hat{B}}_{L}(s){\hat{B}}_{L}(0)\rangle \langle {\hat{B}}_{R}^{\lambda }(s){\hat{B}}_{R}^{-\lambda }(0)\rangle =\displaystyle \frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{M}_{L}({\omega }_{0}-\omega ){M}_{R}^{\lambda ;-\lambda }(\omega ),\\ {k}_{0\to 1}^{\lambda } & = & {\int }_{-\infty }^{\infty }{\rm{d}}{s}\,{{\rm{e}}}^{-{\rm{i}}{\omega }_{0}s}\langle {\hat{B}}_{L}(s){\hat{B}}_{L}(0)\rangle \langle {\hat{B}}_{R}^{\lambda }(s){\hat{B}}_{R}^{-\lambda }(0)\rangle =\displaystyle \frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{M}_{L}(-{\omega }_{0}-\omega ){M}_{R}^{\lambda ;-\lambda }(\omega ).\end{eqnarray} \tag{ 40 }$$

From equation (19) we immediately build the energy current and noise power expressions (33)

$$\begin{eqnarray*}&&\langle {J}_{R}\rangle =\displaystyle \frac{1}{2\pi }\left[{p}_{0}^{{ss}}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,\omega \,{M}_{L}(-{\omega }_{0}+\omega ){M}_{R}(-\omega )-{p}_{1}^{{ss}}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,\omega \,{M}_{L}({\omega }_{0}-\omega ){M}_{R}(\omega )\right],\end{eqnarray*}$$

$$\begin{eqnarray}\langle {S}_{R}\rangle & = & \displaystyle \frac{2}{{k}_{1\to 0}+{k}_{0\to 1}}\left(\displaystyle \frac{1}{4{\pi }^{2}}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,\omega \,{M}_{L}({\omega }_{0}-\omega ){M}_{R}(\omega ){\int }_{-\infty }^{\infty }{\rm{d}}{\omega }^{{\prime} }\,{\omega }^{{\prime} }{M}_{L}(-{\omega }_{0}-{\omega }^{{\prime} }){M}_{R}({\omega }^{{\prime} })-\langle {J}_{R}{\rangle }^{2}\right)\\ & & +\displaystyle \frac{{p}_{0}^{{ss}}}{2\pi }{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{\omega }^{2}{M}_{L}(-{\omega }_{0}-\omega ){M}_{R}(\omega )+\displaystyle \frac{{p}_{1}^{{ss}}}{2\pi }{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{\omega }^{2}{M}_{L}({\omega }_{0}-\omega ){M}_{R}(\omega ).\end{eqnarray} \tag{ 41 }$$

This energy current can be deduced from equation (26). It agrees with the expression used ad hoc in [35] and with the result of [38], where counting was done by resolving the population dynamics. It can be rationalized by interpreting $\omega \,{M}_{L}({\omega }_{0}-\omega ){M}_{R}(\omega ){p}_{1}^{{ss}}$ as a relaxation process with the energy ω emitted to the R reservoir, and the amount of ω₀ − ω disposed into the L-reservoir. The baths thus work cooperatively. A similar interpretation holds for the other term. It is significant to note that this expression has been achieved under relatively general conditions, without specifying the nature of the bath coupling operator, aside from assuming it is in a product form.

Back to equation (37), we set ${\hat{V}}_{{\rm{B}}}={\sum }_{\mu ,j}{\gamma }_{\mu ,j}({\hat{b}}_{\mu ,j}^{\dagger }+{\hat{b}}_{\mu ,j})$ and interrogate the strong-system–bath coupling limit by performing the polaron transformation [77], ${\hat{H}}_{P}={\hat{P}}^{\dagger }\hat{H}\hat{P}$ with $\hat{P}={{\rm{e}}}^{{\rm{i}}{\hat{\sigma }}_{z}\hat{{\rm{\Omega }}}/2}$,

$$\begin{eqnarray}&&{\hat{H}}_{{\rm{NADD}}}=\displaystyle \frac{{\rm{\Delta }}}{2}{\hat{\sigma }}_{z}+\displaystyle \frac{{\omega }_{0}}{2}({\hat{\sigma }}_{+}{{\rm{e}}}^{{\rm{i}}\hat{{\rm{\Omega }}}}+{\hat{\sigma }}_{-}{{\rm{e}}}^{-{\rm{i}}\hat{{\rm{\Omega }}}})+\displaystyle \sum _{\mu ,j}{\omega }_{\mu ,j}{\hat{b}}_{\mu ,j}^{\dagger }{\hat{b}}_{\mu ,j},\end{eqnarray} \tag{ 42 }$$

where ${\hat{\sigma }}_{\pm }=\tfrac{1}{2}({\hat{\sigma }}_{x}\pm {\rm{i}}{\hat{\sigma }}_{y})$, $\hat{{\rm{\Omega }}}={\sum }_{\mu }{\hat{{\rm{\Omega }}}}_{\mu }$, and ${\hat{{\rm{\Omega }}}}_{\mu }=2{\rm{i}}{\sum }_{j}\tfrac{{\gamma }_{\mu ,j}}{{\omega }_{\mu ,j}}({\hat{b}}_{\mu ,j}^{\dagger }-{\hat{b}}_{\mu ,j})$. The tunneling splitting ω₀ is thus dressed by a product of shift operators, ${{\rm{\Pi }}}_{\mu =L,R}\exp \left[-2{\sum }_{j}\tfrac{{\gamma }_{\mu ,j}}{{\omega }_{\mu ,j}}({\hat{b}}_{\mu ,j}^{\dagger }-{\hat{b}}_{\mu ,j})\right]$, a non-additive bath coupling model. The CGF of the model is given by equation (4.1), with the correlation functions

$$\begin{eqnarray}&&{M}_{\mu }(t)=\displaystyle \frac{{\omega }_{0}}{2}{{\rm{e}}}^{-{{ \mathcal Q }}_{\mu }(t)},\end{eqnarray} \tag{ 43 }$$

where the real and imaginary parts, ${{ \mathcal Q }}_{\mu }(t)={{{ \mathcal Q }}_{\mu }^{{\prime} }}_{}(t)+{\rm{i}}{{ \mathcal Q }}_{\mu }(t)^{\prime\prime}$ fulfill

$$\begin{eqnarray}{{ \mathcal Q }}_{\mu }^{{\prime\prime} }{}_{}(t) & = & 2{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{{\rm{\Gamma }}}_{\mu }(\omega )}{\pi {\omega }^{2}}\sin (\omega t){\rm{d}}\omega ,\\ {{ \mathcal Q }}_{\mu }^{{\prime} }{}_{}(t) & = & 2{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{{\rm{\Gamma }}}_{\mu }(\omega )}{\pi {\omega }^{2}}[1-\cos (\omega t)][1+2{n}_{\mu }(\omega )]\,{\rm{d}}\omega .\end{eqnarray} \tag{ 44 }$$

The energy current of the model is given by equation (41), which can be solved analytically in e.g. the so-called Marcus (strong coupling high temperature) limit [35, 38]. Extensions to this result were discussed in [39, 40], going beyond the assumption of $\langle {\hat{V}}_{{\rm{B}}}\rangle =0$. Other studies had employed the polaron picture as a starting point for higher order perturbative treatments [78, 79]. We display the CGF in figure 4, and exemplify the current and its noise in figure 5, exposing a thermal diode effect. Panel (a) in figure 5 also reveals the turnover effect of strong system–bath coupling, as a symmetric junction with α_μ = 0.4 shows a smaller current than with α_μ = 0.2. The behavior of the current noise is quite interesting as it may increase or decrease with ΔT.

Based on the formalism organized in this paper, we are now ready to derive the working equations of the qubit-refrigerator, which was described in [41]. Most interestingly, a non-additive bath coupling model can realize the qubit-QAR, but this function is missing in the additive case.

An autonomous absorption refrigerator transfers thermal energy from a cold (C) bath to a hot (H) bath without an input power by using thermal energy provided from a so-called work (W) reservoir, where T_W > T_H > T_C. A common design of a QAR consists of a three-level system, where each transition between a pair of levels is coupled to only one of the three thermal baths [10, 80, 81]. By tuning the level spacing of the three-level system, one can operate the system as a refrigerator, extract energy from the cold bath, and dump it into the hot one. The three-level QAR operates optimally when system–bath coupling is weak.

In a recent study [41], we demonstrated that the smallest system, a qubit, is incapable of operating as a refrigerator when the baths are coupled additively to the system. The energy current from the cold bath in this case can be derived from the expressions in section 4.2.2, resulting in

$$\begin{eqnarray}&&\langle {J}_{C}\rangle =-\displaystyle \frac{{\omega }_{0}}{{k}_{0\to 1}+{k}_{1\to 0}}\left[{k}_{1\to 0}^{C}{k}_{1\to 0}^{H}\left(\displaystyle \frac{{k}_{0\to 1}^{H}}{{k}_{1\to 0}^{H}}-\displaystyle \frac{{k}_{0\to 1}^{C}}{{k}_{1\to 0}^{C}}\right)+{k}_{1\to 0}^{C}{k}_{1\to 0}^{W}\left(\displaystyle \frac{{k}_{0\to 1}^{W}}{{k}_{1\to 0}^{W}}-\displaystyle \frac{{k}_{0\to 1}^{C}}{{k}_{1\to 0}^{C}}\right)\right].\end{eqnarray} \tag{ 45 }$$

Using the detailed balance relation and the fact that $({{\rm{e}}}^{-{\beta }_{H,W}{\omega }_{0}}-{{\rm{e}}}^{-{\beta }_{C}{\omega }_{0}})\gt 0$, we conclude that $\langle {J}_{C}\rangle \lt 0$ regardless of the details of the model. Equation (45) shows that under the additive model at weak-coupling, every two reservoirs exchange energy independently and thus thermal energy always flows to the colder bath, and the chiller performance is unattainable. In contrast, in [41] we showed that by coupling the system in a non-additive manner to three heat baths which are spectrally structured, a QAR could be achieved. Moreover, we demonstrated that the system could reach the Carnot bound if the reservoirs were characterized by a single frequency component [41].

The groundwork of the qubit-QAR of [41] is presented in this paper. The design is based on a three-bath non-additive interaction form, ${\hat{V}}_{{\rm{B}}}={\hat{B}}_{H}\otimes {\hat{B}}_{C}\otimes {\hat{B}}_{W}$ with $\hat{S}={\hat{\sigma }}_{x}$. The three baths are characterized by different spectral properties and maintained at different temperatures. Following section 4, we perform an FCS analysis for each bath, to obtain the current directed to the system from each terminal. For simplicity, we count energy only from C using the counting parameter λ. Transitions between the system's states are dictated by the rate constants

$$\begin{eqnarray}{k}_{1\to 0}^{\lambda } & = & \displaystyle \frac{1}{{(2\pi )}^{2}}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{\rm{d}}\omega ^{\prime} \,{M}_{C}(\omega ){M}_{H}(\omega ^{\prime} ){M}_{W}({\omega }_{0}-\omega -\omega ^{\prime} ),\\ {k}_{0\to 1}^{\lambda } & = & \displaystyle \frac{1}{{(2\pi )}^{2}}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{\rm{d}}\omega ^{\prime} \,{M}_{C}(\omega ){M}_{H}(\omega ^{\prime} ){M}_{W}(-{\omega }_{0}-\omega -\omega ^{\prime} ).\end{eqnarray} \tag{ 46 }$$

The energy current, from the cold bath to the system, is given by equation (33),

$$\begin{eqnarray}\langle {J}_{C}\rangle & = & \displaystyle \frac{1}{{(2\pi )}^{2}}\left[\,{p}_{0}^{{ss}}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{\rm{d}}\omega ^{\prime} \,\omega \,{M}_{C}(-\omega ){M}_{H}(-\omega ^{\prime} ){M}_{W}(-{\omega }_{0}+\omega +\omega ^{\prime} )\right.\\ & & \left.-{p}_{1}^{{ss}}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{\rm{d}}\omega ^{\prime} \,\omega \,{M}_{C}(\omega ){M}_{H}(\omega ^{\prime} ){M}_{W}({\omega }_{0}-\omega -\omega ^{\prime} )\right].\end{eqnarray} \tag{ 47 }$$

Analogous expressions can be written for $\langle {J}_{H}\rangle$ and $\langle {J}_{W}\rangle$. We also calculate the noise power using equation (33) and arrive at,

$$\begin{eqnarray}\langle {S}_{C}\rangle & = & \displaystyle \frac{2}{{k}_{1\to 0}+{k}_{0\to 1}}\left(\displaystyle \frac{1}{{(2\pi )}^{4}}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{\rm{d}}\omega ^{\prime} \omega \,{M}_{C}(\omega ){M}_{H}({\omega }^{{\prime} }){M}_{W}(-{\omega }_{0}-\omega -{\omega }^{{\prime} })\right.\\ & & \left.\times {\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{\rm{d}}\omega ^{\prime} \omega \,{M}_{C}(\omega ){M}_{H}(\omega ^{\prime} ){M}_{W}({\omega }_{0}-\omega -\omega ^{\prime} )-\langle {J}_{C}{\rangle }^{2}\right)\\ & & +\displaystyle \frac{{p}_{0}^{{ss}}}{{(2\pi )}^{2}}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{\rm{d}}\omega ^{\prime} {\omega }^{2}{M}_{C}(\omega ){M}_{H}({\omega }^{{\prime} }){M}_{W}(-{\omega }_{0}-\omega -{\omega }^{{\prime} })\\ & & +\displaystyle \frac{{p}_{1}^{{ss}}}{{(2\pi )}^{2}}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{\rm{d}}\omega ^{\prime} {\omega }^{2}{M}_{C}(\omega ){M}_{H}({\omega }^{{\prime} }){M}_{W}({\omega }_{0}-\omega -{\omega }^{{\prime} }).\end{eqnarray} \tag{ 48 }$$

Based on these expressions, we study the operation of an absorption refrigerator. As was demonstrated in [41], the design is quite robust, and a qubit-chiller can be realized with the bath correlation function M_μ(ω), μ = C, H, W, taking a variety of forms, such as a Heaviside box function. Another convenient form is a bimodal Gaussian function [41],

$$\begin{eqnarray}&&{M}_{\mu }(\omega )=\displaystyle \frac{1}{{\sigma }_{\mu }}({{\rm{e}}}^{-{(\omega -{\theta }_{\mu })}^{2}/2{\sigma }_{\mu }^{2}}+{{\rm{e}}}^{-{(\omega +{\theta }_{\mu })}^{2}/2{\sigma }_{\mu }^{2}}{{\rm{e}}}^{{\beta }_{\mu }\omega }),\end{eqnarray} \tag{ 49 }$$

which by construction satisfies the detailed balance relation. Here, ${\theta }_{\mu }$ is a central frequency that characterizes the spectrum of the μ bath, and σ_μ is a width parameter. When σ → 0, the function collapses to a Dirac delta function—at positive and negative frequencies. By further setting θ_C + θ_W = θ_H, one can analytically prove that the QAR can approach the Carnot bound [41]. In this special limit, the cooling window is given by

$$\begin{eqnarray}&&\left(\displaystyle \frac{{T}_{W}-{T}_{H}}{{T}_{W}-{T}_{C}}\right)\displaystyle \frac{{T}_{C}}{{T}_{H}}\geqslant \displaystyle \frac{{\theta }_{C}}{{\theta }_{H}},\end{eqnarray} \tag{ 50 }$$

which precisely corresponds to the cooling condition as obtained for a three-level or three-qubit QAR—analyzed with a Markovian master equation with additive dissipators [10, 80, 81].

We now demonstrate a cooling performance over a broad range of parameters, beyond the resonance condition and the Delta function (highly engineered-bath) limit analyzed in [41]. In figure 6, we present the bimodal functions M_μ(ω) for the three baths, and consider two cases: (a) we maintain the resonance condition θ_C + θ_W = θ_H, but depart from the standard setting by using an un-structured work reservoir, employing a very broad function M_W(ω). (b) We structure the reservoirs with a smaller width σ_W, but do not satisfy the resonance condition. The current $\langle {J}_{C}\rangle$ is presented in figure 7 for these two cases. Both setups achieve cooling $\langle {J}_{C}\rangle \gt 0$, demonstrating that the refrigerator is robust for a variety of setups. It survives even when the work reservoir is practically structureless. As well, it can operate beyond the strict resonance condition by tuning the width σ_W.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is useful to note that according to equation (50), we identify the optimal cooling windows for panels (a) and (b) as β_H ≥ 0.4 and β_H ≥ 0.5, respectively. Indeed these inequalities serve as good estimates for the cooling window when σ_W is small. Finally, we recall that as long as one of the bath correlation functions is structureless in frequency domain (decays fast in time domain), the overall correlation function, which is a product of the individual time correlation functions (24) dies quickly, justifying the Markov assumption that is underlying this analysis.

We provided the theoretical groundwork for the bath-cooperative qubit-QAR of [41]. More generally, we presented here a rigorous, thermodynamically consistent treatment of quantum energy transport in small systems while going beyond the standard, additive interaction form. Using a FCS approach, we derived a counting-field dependent Redfield-type equation. After the secular approximation, we obtained the CGF of the model, particularly, the averaged energy current and its noise power for both additive and non-additive models, on the same footing.

We exemplified our results on a qubit system: we studied the transport behavior of a two-terminal setup, and the cooling function of a three-terminal QAR model. In the latter case, we demonstrated that the cooling function was robust, surviving for a broad range of baths' parameters.

Applications explored in this work rely on the secular approximation (decoupling population and coherence dynamics). Nevertheless, the formalism could be used beyond that. The counting-field dependent Redfield equation could be employed to examine the role of quantum coherences in energy transport behavior within multi-level quantum systems, and in the operation of QHEs. Analysis of composite heat engine models, e.g., made of several qubits [82, 83], is left for future studies. It is furthermore essential to examine the behavior of the qubit absorption refrigerator with numerically exact approaches and confirm our predictions.

A FCS analysis of heat exchange is paramount for various reasons. Fundamentally, testing whether the CGF satisfies the fluctuation relation immediately reports on the thermodynamic consistency of the employed method. Practically, one can calculate the current and its noise far from equilibrium—directly from the CGF. While most studies in quantum transport and quantum thermodynamics are focused on the calculation of the averaged energy current within a device, evaluating the current noise is critical for estimating the precision of a process [84]. This nonequilibrium fluctuation-dissipation trade off is captured by the thermodynamic uncertainty relation, which is satisfied within our Markovian QME treatment. In future work we will interrogate this relation in strongly-coupled quantum heat machines, beyond the Markovian limit.

DS acknowledges support from an NSERC Discovery Grant and the Canada Research Chair program. The work of HMF was supported by NSERC CGS-M and NSERC PGS-D programs. BKA acknowledges support from the CQIQC at the University of Toronto.

In what follows, we develop the counting-field dependent Redfield equation in details by following the standard derivation of the second order Markovian master equation [46]. In doing so, we highlight on the different approximations involved and elaborate on their impact on the range of systems that can be simulated with this approach.

We begin with equation (12) and follow the standard derivation of the second order Markovian master equation [46]: we integrate this equation and insert the formal solution to ρ^λ(t) back into equation (12). We trace out the bath degrees of freedom, and obtain the reduced density matrix ${\sigma }^{\lambda }(t)\equiv {\mathrm{Tr}}_{{\rm{B}}}[{\rho }^{\lambda }(t)]$, satisfying

$$\begin{eqnarray}{\dot{\sigma }}^{\lambda }(t) & = & -{\displaystyle \int }_{0}^{t}{\mathrm{Tr}}_{{\rm{B}}}[{\hat{H}}_{{\rm{int}}}^{-\lambda }(t){\hat{H}}_{{\rm{int}}}^{-\lambda }(t-s){\rho }^{\lambda }(t-s)]\,{\rm{d}}{s}+{\displaystyle \int }_{0}^{t}{\mathrm{Tr}}_{{\rm{B}}}[{\hat{H}}_{{\rm{int}}}^{-\lambda }(t){\rho }^{\lambda }(t-s){\hat{H}}_{{\rm{int}}}^{\lambda }(t-s)]\,{\rm{d}}{s}\\ & & +{\displaystyle \int }_{0}^{t}{\mathrm{Tr}}_{{\rm{B}}}[{\hat{H}}_{{\rm{int}}}^{-\lambda }(t-s){\rho }^{\lambda }(t-s){\hat{H}}_{{\rm{int}}}^{\lambda }(t)]\,{\rm{d}}{s}-{\displaystyle \int }_{0}^{t}{\mathrm{Tr}}_{{\rm{B}}}[{\rho }^{\lambda }(t-s){\hat{H}}_{{\rm{int}}}^{\lambda }(t-s){\hat{H}}_{{\rm{int}}}^{\lambda }(t)]\,{\rm{d}}{s}.\end{eqnarray} \tag{ A1 }$$

Note the critical assumption made here, namely, $\langle {\hat{H}}_{{\rm{int}}}^{\pm \lambda }\rangle \equiv {\mathrm{Tr}}_{{\rm{B}}}[{\hat{H}}_{{\rm{int}}}^{\pm \lambda }(0)\rho (0)]=0$. This restricts the structure of the interaction operator, or the range of temperatures and interaction energies for which results are valid. For example, in the NESB model under the polaron transformation, see [34, 38, 85], $\langle {\hat{H}}_{{\rm{int}}}^{\pm \lambda }\rangle$ decays exponentially with temperature and system–bath coupling, thus the method as described here is valid in a high temperature and strong-coupling regime. More generally, one could re-organize the total Hamiltonian as $\hat{H}={\hat{H}}_{S}+\langle {\hat{H}}_{{\rm{int}}}\rangle +{\sum }_{\mu =1}^{N}{\hat{H}}_{\mu }+{\hat{H}}_{{\rm{int}}}-\langle {\hat{H}}_{{\rm{int}}}\rangle$, and derive a QME in the eigenbasis of the new system Hamiltonian, $[{\hat{H}}_{S}+\langle {\hat{H}}_{{\rm{int}}}\rangle ]$, see e.g. [39, 40]. The resulting QME correctly describes both the high and low temperature limits of the spin-boson model (within that order in perturbation theory). In this work, however, with the objective of keeping our presentation and results simple and transparent, we ignore the averaged interaction term.

We assume that the initial density matrix takes a factorized form, $\rho (0)=\sigma (0)\otimes {\rho }_{{\rm{B}}}(0)$, ${\rho }_{{\rm{B}}}(0)={\prod }_{\mu =1}^{N}{\rho }_{\mu }$ with the reservoirs prepared in a canonical state at inverse temperature ${\beta }_{\mu }$, ${\rho }_{\mu }(0)=\exp [-{\beta }_{\mu }{\hat{H}}_{\mu }]/{Z}_{\mu }$, and Z_μ as the partition function for the μ bath. Memory effects from an initially correlated unfactorized initial density matrix are relevant on short timescales in the weak-coupling approximation and do not affect long-time behavior once these correlations have died out [86]. We proceed with the Born approximation, ${\rho }^{\lambda }(t-s)\sim {\sigma }^{\lambda }(t-s)\otimes {\rho }_{{\rm{B}}}(t-s)$, which is justified because the system has little effect on the much larger baths. We assume the reservoirs do not change from their initial state over time, ${\rho }_{\mu }\approx {\rho }_{\mu }(0)$.

To progress from equation (A1), we recall that ${\hat{H}}_{{\rm{int}}}^{\lambda }=\gamma {\sum }_{p}{\hat{S}}^{p}\otimes {\hat{V}}_{{\rm{B}}}^{p,\lambda }$ and define the two-time bath correlation functions as $\langle {\hat{V}}_{{\rm{B}}}^{p,\lambda }(t-s){\hat{V}}_{{\rm{B}}}^{p^{\prime} ,\lambda ^{\prime} }(t)\rangle \equiv {\mathrm{Tr}}_{{\rm{B}}}[{\hat{V}}_{{\rm{B}}}^{p,\lambda }(t-s){\hat{V}}_{{\rm{B}}}^{p^{\prime} ,\lambda ^{\prime} }(t){\rho }_{{\rm{B}}}]$. Note that when operators appear with the same sign for the counting-field, it cancels out,

$$\begin{eqnarray}\langle {\hat{V}}_{{\rm{B}}}^{p,\pm \lambda }(t){\hat{V}}_{{\rm{B}}}^{p^{\prime} ,\pm \lambda }(t-s)\rangle & = & {\mathrm{Tr}}_{{\rm{B}}}\left[{{\rm{e}}}^{\pm {\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }/2}\,{\hat{V}}_{{\rm{B}}}^{p}(t)\,{{\rm{e}}}^{\pm {\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }/2}{{\rm{e}}}^{\pm {\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }/2}\,{\hat{V}}_{{\rm{B}}}^{p^{\prime} }(t-s)\,{{\rm{e}}}^{\pm {\rm{i}}{\sum }_{\mu }{\lambda }_{\mu }{\hat{H}}_{\mu }/2}{\rho }_{{\rm{B}}}\right]\\ & = & \langle {\hat{V}}_{{\rm{B}}}^{p}(t){\hat{V}}_{{\rm{B}}}^{p^{\prime} }(t-s)\rangle =\langle {\hat{V}}_{{\rm{B}}}^{p}(s){\hat{V}}_{{\rm{B}}}^{p^{\prime} }(0)\rangle .\end{eqnarray} \tag{ A2 }$$

The second equality arises due to the cyclic properties of the trace operation and the fact that ${\hat{H}}_{\mu }$ and ρ_B commute. The last equality relies on the fact that the bath is stationary. We define this coupling correlation function as,

$$\begin{eqnarray}&&{M}_{{\rm{B}};{ab},{cd}}^{p,\lambda ;{p}^{{\prime} },{\lambda }^{{\prime} }}(s)\equiv {\gamma }^{2}\,{S}_{{ab}}^{p}{S}_{{cd}}^{{p}^{{\prime} }}\,\langle {\hat{V}}_{{\rm{B}}}^{p,\lambda }(s){\hat{V}}_{{\rm{B}}}^{{p}^{{\prime} },{\lambda }^{{\prime} }}(0)\rangle ,\end{eqnarray} \tag{ A3 }$$

and its half Fourier transform,

$$\begin{eqnarray}&&{R}_{{\rm{B}};{ab},{cd}}^{p,\lambda ;{p}^{{\prime} },{\lambda }^{{\prime} (+)}}(\omega )\equiv {\int }_{0}^{\infty }{{\rm{e}}}^{{\rm{i}}\omega s}{M}_{{\rm{B}};{ab},{cd}}^{p,\lambda ;{p}^{{\prime} },{\lambda }^{{\prime} }}(s)\,{\rm{d}}{s},\,\,\,\,\,\,{R}_{{\rm{B}};{ab},{cd}}^{p,\lambda ;{p}^{{\prime} },{\lambda }^{{\prime} (-)}}(\omega )\equiv {\int }_{-\infty }^{0}{{\rm{e}}}^{{\rm{i}}\omega s}{M}_{{\rm{B}};{ab},{cd}}^{p,\lambda ;{p}^{{\prime} },{\lambda }^{{\prime} }}(s)\,{\rm{d}}{s}.\end{eqnarray} \tag{ A4 }$$

The terms ${S}_{{ab}}^{p}$ and ${S}_{{cd}}^{{p}^{{\prime} }}$ are the matrix elements of the ${\hat{S}}^{p}$ and ${\hat{S}}^{{p}^{{\prime} }}$ operators, respectively, in the eigenbasis of the system Hamiltonian ${\hat{H}}_{S}$.

We now study equation (A1) in the Markovian limit, assuming the dynamics of the bath is much faster than those of the system. Firstly, we replace ${\sigma }^{\lambda }(t-s)\to {\sigma }^{\lambda }(t)$ since in the weak-coupling limit, the interaction parameter γ is small. Recall that equation (A1) is written in the interaction representation. Second, we assume that bath correlation functions decay to small values much faster than any characteristic system timescale. This allows us to take the upper limit of the integrals to infinity, where beyond t, the contributions from the integrand due to the bath-bath correlation function, having decayed quickly, are negligible. The validity of the Markovian master equation has been extensively studied [46, 87] and is valid for our purposes as long as γ is small and the (broadband) bath is characterized by resonant modes with the system.

These approximations result in the Markovian QME, the counting-field dependent Redfield-type equation [47],

$$\begin{eqnarray}{\dot{\sigma }}_{{nm}}^{\lambda }(t) & = & -{\rm{i}}{[{\hat{H}}_{S},\sigma (t)]}_{{nm}}+\displaystyle \sum _{p,q,j,k}[-{\sigma }_{{km}}^{\lambda }(t){R}_{{\rm{B}};{nj},{jk}}^{p,{q}^{(+)}}({E}_{k,j})-{\sigma }_{{nj}}^{\lambda }(t){R}_{{\rm{B}};{jk},{km}}^{p,{q}^{(-)}}({E}_{j,k})\\ & & +{\sigma }_{{jk}}^{\lambda }(t)({R}_{{\rm{B}};{km},{nj}}^{q,\lambda ,p,-{\lambda }^{(-)}}({E}_{k,m})+{R}_{{\rm{B}};{km},{nj}}^{q,\lambda ,p,-{\lambda }^{(+)}}({E}_{j,n}))],\end{eqnarray} \tag{ A5 }$$

where ${E}_{m,n}={E}_{m}-{E}_{n}$ and E_n are the energy eigenstates of the system Hamiltonian ${\hat{H}}_{S}$. At this point, we switched back to the Schrödinger representation. Recall that p and q sum over the different operators that couple to the system. The other indices, j, k, count eigenstates of the system. This derivation can be extended to higher orders in the system–bath interaction within the projection operator formalism beyond the lowest order Redfield equation [46].

We continue and simplify equation (A5) by performing the secular approximation so as to decouple population and coherence dynamics. If the characteristic timescale of the subsystem ${\tau }_{{\rm{S}}}\approx {E}_{n,m}^{-1}$ is much shorter than the system's relaxation time τ_R, the function ${{\rm{e}}}^{{\rm{i}}{E}_{n,m}t}$ oscillates rapidly over τ_R and thus the contribution of the coherence terms in the reduced density matrix would be averaged out to zero. The resulting population dynamics, ${p}_{n}(t)\equiv {\sigma }_{{nn}}(t)$, is given in equation (14) in the main text.

To prove the detailed balance relation for a specific bath ν, we examine the correlation function $\langle {\hat{B}}_{\nu }^{\lambda }(s){\hat{B}}_{\nu }^{{\lambda }^{{\prime} }}(0) \rangle ;$ recall that the average is performed with respect to the initial, canonical state ${\rho }_{\nu }={{\rm{e}}}^{-{\beta }_{\nu }{\hat{H}}_{\nu }}/{Z}_{\nu }$. Using the cyclic property of the trace we get

$$\begin{eqnarray} &&\langle {\hat{B}}_{\nu }^{\lambda }(s){\hat{B}}_{\nu }^{{\lambda }^{{\prime} }}(0) \rangle =\displaystyle \frac{1}{{Z}_{\nu }}{{\rm{Tr}}}_{\nu }[{\hat{B}}_{\nu }^{{\lambda }^{{\prime} }}{{\rm{e}}}^{{\rm{i}}{\hat{H}}_{\nu }(s+{\rm{i}}{\beta }_{\nu })}{\hat{B}}_{\nu }^{\lambda }{{\rm{e}}}^{-{\rm{i}}{\hat{H}}_{\nu }(s+{\rm{i}}{\beta }_{\nu })}{{\rm{e}}}^{-{\beta }_{\nu }{\hat{H}}_{\nu }}].\end{eqnarray} \tag{ B1 }$$

From here, we readily organize as the detailed balance relation,

$$\begin{eqnarray} &&\langle {\hat{B}}_{\nu }^{\lambda }(s){\hat{B}}_{\nu }^{{\lambda }^{{\prime} }}(0) \rangle = \langle {\hat{B}}_{\nu }^{{\lambda }^{{\prime} }}(0){\hat{B}}_{\nu }^{\lambda }(s+{\rm{i}}{\beta }_{\nu }) \rangle = \langle {\hat{B}}_{\nu }^{{\lambda }^{{\prime} }}(-s-{\rm{i}}{\beta }_{\nu }){\hat{B}}_{\nu }^{\lambda }(0) \rangle .\end{eqnarray} \tag{ B2 }$$

Note that equation (B2) has the counting-field terms switched. Moreover, we can open up the counting-field dependent terms and achieve

$$\begin{eqnarray*} &&\langle {\hat{B}}_{\nu }^{\lambda }(s){\hat{B}}_{\nu }^{{\lambda }^{{\prime} }}(0) \rangle =\displaystyle \frac{1}{{Z}_{\nu }}{{\rm{Tr}}}_{\nu }[{\hat{B}}_{\nu }{{\rm{e}}}^{{\rm{i}}{\hat{H}}_{\nu }\left[s+{\rm{i}}{\beta }_{\nu }+\tfrac{1}{2}(\lambda -{\lambda }^{{\prime} }\right)]}{\hat{B}}_{\nu }{{\rm{e}}}^{-{\rm{i}}{\hat{H}}_{\nu }[s+{\rm{i}}{\beta }_{\nu }+\tfrac{1}{2}(\lambda -{\lambda }^{{\prime} })]}{{\rm{e}}}^{-{\beta }_{\nu }{\hat{H}}_{\nu }}],\end{eqnarray*}$$

which yields,

$$\begin{eqnarray} \langle {\hat{B}}_{\nu }^{\lambda }(s){\hat{B}}_{\nu }^{{\lambda }^{{\prime} }}(0) \rangle & = & \langle {\hat{B}}_{\nu }(0){\hat{B}}_{\nu }(s+{\rm{i}}{\beta }_{\nu }+(\lambda -{\lambda }^{{\prime} })/2) \rangle \\ & = & \langle {\hat{B}}_{\nu }(-s-{\rm{i}}{\beta }_{\nu }-(\lambda -{\lambda }^{{\prime} })/2){\hat{B}}_{\nu }(0) \rangle .\end{eqnarray} \tag{ B3 }$$

From here, we immediately deduce that

$$\begin{eqnarray} &&\langle {\hat{B}}_{\nu }^{\lambda }(s){\hat{B}}_{\nu }^{-\lambda }(0) \rangle = \langle {\hat{B}}_{\nu }(-s-{\rm{i}}{\beta }_{\nu }-\lambda ){\hat{B}}_{\nu }(0) \rangle ,\end{eqnarray} \tag{ B4 }$$

receiving the detailed balance result,

$$\begin{eqnarray}&&{M}_{\nu ;{nj},{jn}}^{\lambda ;-\lambda }(\omega )={{\rm{e}}}^{({\beta }_{\nu }-{\rm{i}}\lambda )\omega }{M}_{\nu ;{nj},{jn}}(-\omega ),\end{eqnarray} \tag{ B5 }$$

which is essential to much of our analysis.

The steady state fluctuation theorem for entropy production is a microscopic statement of the second law of thermodynamics [42, 43]. It is crucial to validate it here, so as to establish the thermodynamic consistency of our treatment. we prove it by verifying the following symmetry of the CGF,

$$\begin{eqnarray}&&{\mu }_{+}^{\lambda }={\mu }_{+}^{-\lambda +{\rm{i}}{\rm{\Delta }}\beta },\end{eqnarray} \tag{ C1 }$$

with ${\rm{\Delta }}\beta ={\beta }_{L}-{\beta }_{R}$. This symmetry, in fact, holds for both eigenvalues in equation (4.1). We prove it by showing that

$$\begin{eqnarray}&&{k}_{0\to 1}^{\lambda }{k}_{1\to 0}^{\lambda }={k}_{0\to 1}^{-\lambda +{\rm{i}}{\rm{\Delta }}\beta }{k}_{1\to 0}^{-\lambda +{\rm{i}}{\rm{\Delta }}\beta },\end{eqnarray} \tag{ C2 }$$

for both additive and non-additive interaction models. First, we recall that for an additive system–bath interaction model, the rate constants satisfy ${k}_{1\to 0}^{\lambda }={k}_{1\to 0}^{L}+{k}_{1\to 0}^{R,\lambda }$, where

$$\begin{eqnarray}{k}_{0\to 1}^{R,\lambda } & = & {\int }_{-\infty }^{\infty }{\rm{d}}{s}\,{{\rm{e}}}^{-{\rm{i}}{\omega }_{0}s}\langle {\hat{B}}_{R}(s){\hat{B}}_{R}(0)\rangle \,{{\rm{e}}}^{{\rm{i}}{\omega }_{0}\lambda },\\ {k}_{1\to 0}^{R,\lambda } & = & {\int }_{-\infty }^{\infty }{\rm{d}}{s}\,{{\rm{e}}}^{{\rm{i}}{\omega }_{0}s}\langle {\hat{B}}_{R}(s){\hat{B}}_{R}(0)\rangle \,{{\rm{e}}}^{-{\rm{i}}{\omega }_{0}\lambda }.\end{eqnarray} \tag{ C3 }$$

The L-reservoir induced rates take the same form but with λ = 0. Using the detailed balance relation, ${k}_{0\to 1}^{\mu }={{\rm{e}}}^{-{\beta }_{\mu }{\omega }_{0}}{k}_{1\to 0}^{\mu }$, we can readily show that

$$\begin{eqnarray}&&{k}_{0\to 1}^{-\lambda +{\rm{i}}{\rm{\Delta }}\beta }=({k}_{1\to 0}^{L}+{k}_{1\to 0}^{R,\lambda }){{\rm{e}}}^{-{\omega }_{0}{\beta }_{L}},\,\,\ {k}_{1\to 0}^{-\lambda +{\rm{i}}{\rm{\Delta }}\beta }=({k}_{0\to 1}^{L}+{k}_{0\to 1}^{R,\lambda }){{\rm{e}}}^{{\omega }_{0}{\beta }_{L}},\end{eqnarray} \tag{ C4 }$$

which proves equation (C2). It is easy to generalize this proof for systems including more than two states.

In the non-additive case, the rate constants satisfy

$$\begin{eqnarray}{k}_{0\to 1}^{\lambda } & = & {\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{M}_{L}(-{\omega }_{0}-\omega ){M}_{R}(\omega )\,{{\rm{e}}}^{-{\rm{i}}\omega \lambda },\\ {k}_{1\to 0}^{\lambda } & = & {\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{M}_{L}({\omega }_{0}-\omega ){M}_{R}(\omega )\,{{\rm{e}}}^{-{\rm{i}}\omega \lambda }.\end{eqnarray} \tag{ C5 }$$

Using the detailed balance relation, which holds for each component separately, ${M}_{\mu }(\omega )={{\rm{e}}}^{\omega {\beta }_{\mu }}{M}_{\mu }(-\omega )$, we confirm that

$$\begin{eqnarray}{k}_{0\to 1}^{-\lambda +{\rm{i}}{\rm{\Delta }}\beta } & = & {{\rm{e}}}^{-{\beta }_{L}{\omega }_{0}}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{M}_{L}({\omega }_{0}-\omega ){M}_{R}(\omega ){{\rm{e}}}^{-{\rm{i}}\omega \lambda }\\ {k}_{1\to 0}^{-\lambda +{\rm{i}}{\rm{\Delta }}\beta } & = & {{\rm{e}}}^{{\beta }_{L}{\omega }_{0}}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,{M}_{L}(-{\omega }_{0}-\omega ){M}_{R}(\omega ){{\rm{e}}}^{-{\rm{i}}\omega \lambda },\end{eqnarray} \tag{ C6 }$$

which immediately proves equation (C2). A more general proof, for systems with more than two states, was given in [85]. It is straightforward to generalize these results, to describe quantum transport in three-terminal setups. In such models, the system can act as a QAR [10], as discussed in section 5.