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

We formulate a full-counting statistics description to study energy exchange in multi-terminal junctions. Our approach applies to quantum systems that are coupled either additively or non-additively (cooperatively) to multiple reservoirs. We derive a Markovian Redfield-type equation for the counting-field dependent reduced density operator. Under the secular approximation, we confirm that the cumulant generating function satisfies the heat exchange fluctuation theorem. Our treatment thus respects the second law of thermodynamics. We exemplify our formalism on a multi-terminal two-level quantum system, and apply it to realize the smallest quantum absorption refrigerator, operating through engineered reservoirs, and achievable only through a cooperative bath interaction model.


I. INTRODUCTION
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][2][3]. What is "quantum" in quantum heat engines [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][8][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][12][13][14], quantum correlations [15] and reservoir engineering e.g. by squeezing [16][17][18][19][20], on the performance of heat engines to potentially overcome classical bounds.
A central aspect of nanoscale quantum heat engines 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 weakcoupling 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][22][23][24][25][26][27][28][29][30][31][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][36][37], a strong-coupling treatment administers a turnover behavior [35,36,[38][39][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 Ref. [41] and illustrated earlier on in Refs. [35,36,38], the dynamics of non-additive models can be treated with standard kinetic quantum master equations, 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 Ref. [41]. Our formalism utilizes a full-counting statistics (FCS) approach that provides cumulants of energy exchange to all orders [26,[42][43][44][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 Sec. II. Sec. III 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 Sec. IV, we exemplify our analytical results on a two-state system and present simulation results of the two-terminal nonequilibrium spin-boson model. The application of our formalism to study a quantum absorption refrigerator is given in Sec. V. We summarize our work in Sec. VI. Throughout the paper we work with units where = 1, k B = 1.

II. ADDITIVE AND NON-ADDITIVE INTERACTION MODELS
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-quantumsystem Hamiltonian, with a systemĤ S coupled to N reservoirs iŝ Here,Ŝ ν is a system operator that is coupled to the ν bath's degrees of freedom andV ν 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 quantum master equation (QME) ( ≡ 1) with the dissipators D ν [σ(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 Eq. (2), the rate of energy relaxation from the system can be written as We can readily identify the currents flowing towards the system from the different baths as 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 systembath interaction models, shown as schemes in Fig. ( Here,Ŝ p is a system operator,B 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 operatorsB p ν , as in Eq. (1). The second model is non-additive (NADD), or multiplicative, inB 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.
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][53][54][55][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 Ref. [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 quantum absorption refrigerator, with a qubit as its working fluid analogue [41]. A cooperative bath behavior is what allows for refrigeration.

III. METHOD: FULL-COUNTING STATISTICS FOR ENERGY EXCHANGE
We present a projection operator formalism [58,59] for deriving the characteristic function for quantum energy exchange between a system and multiple attached thermal reservoirs. In steady state, we ultimately obtain the cumulant generating function 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 quantum master equation-the Redfield equation-generalized to describe the full-counting statistics 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.

A. Hamiltonian
Our discussion begins with a general open quantum system Hamiltonian, a systemĤ S coupled to N reservoirs, Here,Ŝ p is a system operator andV p N 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Ŝ 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 reservoir ν as a change in the ν bath energy J ν (t) = − dĤν (t) dt . Operators are written in the Heisenberg representation,Â(t) =Û † (t)ÂÛ (t), whereÛ (t) = e −iĤt . Therefore, the total energy exchange at the ν terminal is given by the integrated current Employing the two-time measurement protocol [42,44], the characteristic function for energy transfer at multiple baths takes the form 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 Derivatives of G(λ) with respect to (iλ) yield the steady-state energy current cumulants. We now recast Eq. (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Ĥ ν . We compactly write Here, we define the time evolution 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, withĤ −λ ≡ e −i ν λνĤν /2Ĥ e i ν λνĤν /2 . SinceĤ ν commutes with all operators on the right hand side of Eq. (4), exceptV p N , we obtain a counting-dependent total Hamiltonian aŝ whereV p,−λ N = e −i ν λνĤν /2V p N e i ν λνĤν /2 . In order to ultimately reach the CGF, we describe the counting-field dependent Nakajima-Zwanzig equation in Appendix A, and present in Sec. III B a second order quantum master equation to obtain the counting-field dependent density matrix in Eq. (8).
B. Counting-field dependent second-order quantum master equation Continuing from Eq. (8), we switch to the interaction representation and define the counting-field density matrix, Operators are now written in the interaction (11) can be written as a differential equation,ρ λ 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 σ λ (t) ≡ Tr B ρ λ (t) , satisfying (Schrödinger representation), Here, E m,n = E m − E n and E n are the energy eigenstates of the system HamiltonianĤ 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. Eq. (13) can be further simplified by performing the secular approximation, so as to decouple population and coherence dynamics. The population dynamics, p n (t) ≡ σ nn (t), then followṡ with The averages are performed with respect to the initial condition. The rate constants are Fourier transforms of the correlation functions, M p,λ;p ,λ N ; ab,cd (ω) ≡ ∞ −∞ e iωs M p,λ;p ,λ N ; ab,cd (s) ds. Eq. (14) is a Markovian-secular counting-field dependent quantum master equation 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 L λ as Analogous derivations of the counting-field dependent reduced density operator, for charge transfer problems, were performed in e.g. Refs. [60][61][62][63][64]. Recall that Eq. (16) gives us the characteristic function Z from Eq. (8) and ultimately the CGF via Eq. (7). While our derivation is standard, our objective here has been to highlight that Eqs. (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 Eq. (4).

C. Energy current
We organize a closed expression for the energy current from Eq. (16) by differentiating G with respect to the counting-field λ. The characteristic function is given by a trace over the system states. From Eq. (8), Z(λ, t) = 1|p λ (t) with 1| = (1, 1, ...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, J ν = lim t→∞ with the steady state populations p ss j , which we get by solving Eq. (16) in the long-time limit. The full-counting statistics approach thus provides a rigorous working expression of the current for ADD and NADD models, at 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 Sec. IV B. 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 [65].
So far, we considered a system coupled to N thermal baths, with N counting parameters, λ 1 , λ 2 , ..., λ ν , ..., λ 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; additiveV N = N ν=1B ν and non-additiveV N = N ν=1B ν . 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 Eq. (14). Finally, we henceforth ignore the γ 2 factor.
The counting-dependent Liouvillian L λ is made of correlation functions for the ν bath (15), Its Fourier transform in frequency domain gives the relation, which is useful for calculating derivatives in Eq. (17). It is also useful to note the detailed balance relation with counting parameters for a specific bath ν at inverse temperature β ν , For details, see Appendix B.

Additive bath interaction model:VN
Recall that baths' operators are written in the interaction representation, the baths are initially uncorrelated and we assume that V N = 0. This results in the rate constants in Eq. (14) being additive in the N reservoirs, The Liouvillian is then given by adding up all contributions, Using Eqs. (14), (17) and (19) we receive for the energy current with L ν 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, dTr[Ĥ S σ(t)] dt = N ν=1 J ν , as discussed in Sec. II. Nevertheless, the formalism presented here can be used to feasibly provide higher order cumulants. We now assume that the system is coupled to the environment according to the following form,V λ . 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 since the bath initial state is factorized, ρ B = N ν=1 ρ ν . In frequency domain, the correlation function turns into a convolution, 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 Eq. (14) and the relation in (19), we receive the energy current from Eq. (17), We emphasize that this result holds without specifying the system, bath or the system-bath interaction Hamiltonian, aside from it being a product form.

A. Cumulant generating function
Spin-bath models are central to the theory of open quantum systems [66,67]. In particular, the spin-boson model has found immense applications in condensed phases physics, chemical dynamics, quantum optics and quantum technologies [66]. It serves to develop approximation schemes perturbative and non-perturbative in the system-bath coupling strength, see for example recent studies [68,69]. Beyond questions over quantum decoherence, dissipation, and thermalization, which can be addressed within the single-bath spin-boson model, the two-bath, nonequilibrium spin-boson (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][36][37][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. Refs. [34,70,71].
In this section, we consider a nonequilibrium spin-bath model: a qubit coupled via aσ x operator to N reservoirs. Our objective is to work with the formalism of Sec. III 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, In the language of the general discussion, S 00 = S 11 = 0, S 01 = S 10 = 1, and E 1,0 = ω 0 , E 0,1 = −ω 0 . The QME (14) thus reduces to .
Note that we have dropped the subscripts ab, cd on M since S 01 S 10 = 1. It is convenient to identify the following four rates, The two eigenvalues of this matrix are The rate matrix in Eq. (28) is the Liouvillian L and is therefore used to calculate the CGF, energy current and noise power.
Recall that the characteristic function Z(λ, t) is the trace over the counting-field dependent reduced density matrix, and that it can be used to obtain the CGF, where 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 Note that µ λ=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.
Let us now consider the expansion 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 1 ), noise power (2a 2 ), skewness, etc.

B. 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 In the long-time limit, we solve Eq. (28) for λ = 0 and find the steady state population p ss 0 = k1→0 k1→0+k0→1 , with p 0 + p 1 = 1. We then reach the following expressions, It is significant to note that these results are general, valid for arbitrary bath coupling operatorV N . In what follows we discuss the so-called "thermodynamic uncertainty relation", then calculate the current and noise in the ADD and NADD models.

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 [72][73][74], 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, Σ Q = J ∆β with ∆β = 1 T L − 1 T R . The thermodynamic uncertainty relation then reduces to In linear response, J R = κ∆T , ∆β = ∆T /T 2 a , with 2T a = T L + T R twice the averaged temperature. The inequality then reaches the Green-Kubo relation, S R = 2k B T 2 a κ, with κ as the thermal conductance. In Secs. IV B 2 and IV B 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 Ref. [75].

Additive bath interaction model:VN =BL +BR
In the case of an additive system-bath coupling,V N =B L +B R , it can be shown that the rate constants are given by (we count energy at the R terminal) Explicitly, one can readily prove that the counting field rate constant relates to the λ = 0 expression as k R,λ 0→1 = e iω0λ k R 0→1 , shown previously in Eq. (19). Similarly, k R,λ 1→0 = e −iω0λ k R 1→0 . Using these two relations in (33), the energy current and noise power reduce to 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 below Eq. (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,76]. The noise expression agrees with the zero-frequency limit of Ref. [77]. Let us consider the nonequilibrium spin-boson model, The system is given by a spin, and the thermal baths are modeled by a collection of noninteracting bosons, Here,σ x andσ z are the Pauli matrices as in Eq. (27), ∆ is the energy gap between the spin levels, and ω 0 is the tunnelling energy. The two reservoirs, L and R, include uncoupled harmonic oscillators, withb † ν,j (b ν,j ) as the bosonic creation (annihilation) operator of the jth mode in the ν reservoir. We solve the transport behavior of the model in the separable bath interactionV N limit,V N = ν,j γ ν,j (b † ν,j +b ν,j ), γ ν,j is the system-bath interaction energy. For simplicity, let us take ∆ = 0 in Eq. (37). We then perform a rotation,Ŵ †σ zŴ =σ x ,Ŵ †σ xŴ =σ z , witĥ W = 1 √ 2 (σ x +σ z ), and receive the transformed HamiltonianĤ ADD =Ŵ †ĤŴ , This Hamiltonian offers a convenient starting point for a perturbative treatment. Using the results of Sec. IV A, one can readily write down the CGF of the model (30) with k ν 1→0 = Γ ν (ω 0 )(1 + n ν (ω 0 )), and k ν 0→1 = Γ ν (ω 0 )n ν (ω 0 ), the results for which are shown in Fig. 2. Here, n ν (ω) = (e βν ω − 1) −1 is the Bose-Einstein distribution function, Γ ν (ω) = 2π j∈ν γ 2 ν,j δ(ω − ω ν,j ) is the spectral density function. For details, see [38]. The energy current and noise power at the R terminal are Other approaches for treating this model in a perturbative manner are formulated in e.g. Refs. [49,70,78]. We simulate both symmetric and asymmetric junctions in Fig. 2, and demonstrate the thermal rectification effect, an asymmetry of the current and noise with respect to ±∆T . In Fig. 3 we present the current and its noise for the weak-coupling additive model (39). We further demonstrate that the thermodynamic uncertainty relation, Eq. (35), is satisfied in both the classical T a ω 0 and quantum T a ω 0 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 Γ ν (ω) = α ν ωe −|ω|/ων with α ν as a dimensionless coefficient and ω ν as the cutoff frequency.

Non-additive bath interaction:VN =BL ⊗BR
We now discuss the NADD modelV N =B L ⊗B R . In this case, the Liouvillian cannot be separated to left and right-reservoir assisted processes. For a two-level system, the rates, Eq. (25), simplify to From Eq. (19) we immediately build the energy current and noise power expressions (33) This energy current can be deduced from Eq. (26). It agrees with the expression used ad hoc in Ref. [35] and with the result of Ref. [38], where counting was done by resolving the population dynamics. It can be rationalized by interpreting ω M L (ω 0 − ω)M R (ω) p ss 1 as a relaxation process with the energy ω emitted to the R reservoir, and the amount of ω 0 − ω 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 Eq. (37), we setV N = ν,j γ ν,j (b † ν,j +b ν,j ) and interrogate the strong-system bath coupling limit by performing the polaron transformation [79],Ĥ P =P †ĤP withP = e iσzΩ/2 , whereσ ± = 1 2 (σ x ± iσ y ),Ω = νΩ ν , andΩ ν = 2i j γν,j ων,j (b † ν,j −b ν,j ). The tunneling splitting ω 0 is thus dressed by a product of shift operators, Π ν=L,R exp −2 j γν,j ων,j (b † ν,j −b ν,j ) , a non-additive bath coupling model. The CGF of the model is given by Eq. (29), with the correlation functions where the real and imaginary parts, The energy current of the model is given by Eq. (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 Refs. [39,40], going beyond the assumption of V N = 0. Other studies had employed the polaron picture as a starting point for higher order perturbative treatments [80,81]. We display the CGF in Fig. 4, and exemplify the current and its noise in Fig.  5, exposing a thermal diode effect. Panel (a) in Fig. 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 .  (42). Parameters are the same as in Fig. 4. In (c), the lower bound of 2 is shown as a dashed line.

V. QUANTUM ABSORPTION REFRIGERATOR
Based on the formalism organized in this paper, we are now ready to derive the working equations of the qubitrefrigerator, which was described in Ref. [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,82,83]. 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 Sec. IV B 2, resulting in Using the detailed balance relation and the fact that (e −β H,W ω0 − e −β C ω0 ) > 0, we conclude that J C < 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 Ref. [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 Ref. [41] is presented in this paper. The design is based on a three-bath non-additive interaction form,V N =B H ⊗B C ⊗B W withŜ =σ x . The three baths are characterized by different spectral properties and maintained at different temperatures. Following Sec. IV, 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 The energy current, from the cold bath to the system, is given by Eq. (33), Analogous expressions can be written for J H and J W . We also calculate the noise power using Eq. (33) and arrive at, Based on these expressions, we study the operation of an absorption refrigerator. As was demonstrated in Ref. [41], the design is quite robust, and a qubit-chiller can be realized with the bath correlation function M ν (ω) taking a variety of forms, such as a Heaviside box function. Another convenient form is a bimodal Gaussian function [41], which by construction satisfies the detailed balance relation. Here, θ ν is a central frequency that characterizes the spectrum, σ ν 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 defined by 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,82,83]. 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 Ref. [41]. In Fig. 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 J C is presented in Fig. 7 for these two cases. Both setups achieve cooling J C > 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 .
It is useful to note that according to Eq. (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.

VI. SUMMARY
We provided the theoretical groundwork for the bath-cooperative qubit-QAR of Ref. [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 full-counting statistics approach, we derived a counting-field dependent Redfield-type equation. After the secular approximation, we obtained the cumulant generating function of the model, particularly, the averaged energy current and its noise power for both additive and non-additive models, on the same footing. Our work illustrated that by studying non-additive models within a weak coupling method one could capture strong-coupling (e.g. multi-phonon) effects that are not realized with an additive weak-coupling treatment.
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 quantum absorption refrigerator 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 quantum heat engines. Analysis of composite heat engine models, e.g., made of several qubits [84,85], 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 full-counting statistics analysis of heat exchange is paramount for various reasons. Fundamentally, testing whether the cumulant generating function 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 cumulant generating function. 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 [86]. 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. The formally exact Nakajima-Zwanzig generalized quantum master equation (GQME) [47] can be extended to include counting information, as we now explain. We consider an open quantum system Hamiltonian,Ĥ =Ĥ S + H B +Ĥ int , with the total density operator ρ(t). In the absence of counting parameters, the quantum Liouville equation can be written as (Schrödinger picture),ρ where L mnm n = H mm δ nn − δ mm H n n . Following Eq. (12), we generalize this standard definition to our needs, where L λ mnm n = H −λ mm δ nn − δ mm H λ n n .
The counting-field dependent HamiltonianĤ ±λ is defined below Eq. (9). Using projection operators, we proceed and write down the exact GQME [46,47], with the kernel Here the system Liouvillian is L S (·) = [Ĥ S , ·] and L λ int is defined as in Eq. (A3). To arrive at this GQME we assumed a factorized initial condition ρ(0) = σ(0) ⊗ ρ B (0), and that Ĥ ±λ int ≡ Tr B Ĥ ±λ int (0)ρ(0) = 0. The projection operatorŝ Q = 1−P are chosen such that we project out the degrees of freedom of the baths,P(·) = ρ B (0)Tr B (·). Eqs. (A4)-(A5) are simple generalizations of the exact Nakajima-Zwanzig quantum master equation to include counting information. Approximate expressions, developed to solve this formal result, see e.g. Ref. [87], can be readily generalized to the present case. In particular, in the weak coupling limit the exact kernel reduces to Once further performing the Markovian approximation, Eq. (A4) with (A6) precisely recovers the Redfield equation with counting parameters (13).
In what follows, we develop the counting-field dependent Redfield equation in more 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 Eq. (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 Eq. (12). We trace out the bath degrees of freedom, and obtain the reduced density matrix σ λ (t) ≡ Tr B ρ λ (t) , satisfyinġ Note the critical assumption made here, namely, Ĥ ±λ int ≡ Tr B Ĥ ±λ int (0)ρ(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 nonequilibrium spin-boson model under the polaron transformation, see Refs. [34,38,88], Ĥ ±λ int 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 aŝ H =Ĥ S + Ĥ int + N ν=1Ĥ ν +Ĥ int − Ĥ int , and derive a quantum master equation in the eigenbasis of the new system Hamiltonian, Ĥ S + Ĥ int , see e.g. Refs. [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, ρ(0) = σ(0) ⊗ ρ B (0), ρ B (0) = N ν=1 ρ ν with the reservoirs prepared in a canonical state at inverse temperature β ν , ρ ν (0) = exp[−β νĤν ]/Z ν , 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 [89]. We proceed with the Born approximation, ρ λ (t − s) ∼ σ λ (t − s) ⊗ ρ 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, ρ ν ≈ ρ ν (0).
To progress from Eq. (A7), we recall thatĤ λ int = γ pŜ p ⊗V p,λ N and define the two-time bath correlation functions as V p,λ Note that when operators appear with the same sign for the counting-field, it cancels out, The second equality arises due to the cyclic properties of the trace operation and the fact thatĤ ν and ρ B commute. The last equality relies on the fact that the bath is stationary. We define this coupling correlation function as, M p,λ;p ,λ N ; ab,cd (s) ≡ γ 2 S p ab S p cd V p,λ N (s)V p ,λ The terms S p ab and S p cd are the matrix elements of theŜ p andŜ p operators respectively, in the eigenbasis of the system HamiltonianĤ S .
We now study Eq. (A7) in the Markovian limit, assuming the dynamics of the bath is much faster than those of the system. Firstly, we replace σ λ (t − s) → σ λ (t) since in the weak coupling limit, the interaction parameter γ is small. Recall that Eq. (A7) 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,90] 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], σ λ nm (t) = −i[Ĥ S , σ(t)] nm + p,q,j,k − σ λ km (t)R p,q (+) N ; nj,jk (E k,j ) − σ λ nj (t)R p,q (−) N ; jk,km (E j,k ) + σ λ jk (t) R q,λ,p,−λ (−) N ; km,nj (E k,m ) + R q,λ,p,−λ (+) N ; km,nj (E j,n ) , where E m,n = E m − E n and E n are the energy eigenstates of the system HamiltonianĤ 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. We continue and simplify Eq. (A11) by performing the secular approximation so as to decouple population and coherence dynamics. If the characteristic timescale of the subsystem τ S ≈ E −1 n,m is much shorter than the system's relaxation time τ R , the function e i En,mt 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) ≡ σ nn (t), is given in Eq. (14) in the main text.
Appendix C: Fluctuation relation for two level system under two-terminal transport 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, with ∆β = β L − β R . This symmetry, in fact, holds for both eigenvalues in Eq. (29). We prove it by showing that for both additive and non-additive interaction models. First, we recall that for an additive system-bath interaction model, the rate constants satisfy k λ The L-reservoir induced rates take the same form but with λ = 0. Using the detailed balance relation, k ν 0→1 = e −βν ω0 k ν 1→0 , we can readily show that which proves Eq. (C2). It is easy to generalize this proof for systems including more than two states. In the non-additive case, the rate constants satisfy Using the detailed balance relation, which holds for each component separately, M ν (ω) = e ωβν M ν (−ω), we confirm that which immediately proves Eq. (C2). A more general proof, for systems with more than two states, was given in Ref. [88]. It is straightforward to generalize these results, to describe quantum transport in three-terminal setups. In such models, the system can act as a quantum absorption refrigerator [10], as discussed in Sec. V.