Energy current and its statistics in the nonequilibrium spin-boson model: Majorana fermion representation

We study the statistics of thermal energy transfer in the nonequilibrium (two-bath) spin-boson model. This quantum many-body impurity system serves as a canonical model for quantum energy transport. Our method makes use of the Majorana fermion representation for the spin operators, in combination with the Keldysh nonequilibrium Green's function approach. We derive an analytical expression for the cumulant generating function of the model in the steady state limit, and show that it satisfies the Gallavotti-Cohen fluctuation symmetry. We obtain analytical expressions for the heat current and its noise, valid beyond the sequential and the co-tunnelling regimes. Our results satisfy the quantum mechanical bound for heat current in interacting nanojunctions. Results are compared with other approximate theories, as well as with a non-interacting model, a fully harmonic thermal junction.


I. INTRODUCTION
The spin-boson (SB) model comprises a two-state system (spin) interacting with a dissipative thermal environment, a collection of harmonic modes. It is one of the (conceptually) simplest, yet non-trivial models in the theory of open quantum systems [1,2]. The model has found diverse applications in condensed phases physics, chemical dynamics, and quantum optics. In particular, it offers a rich platform for studying complex physical processes such as dissipative spin-dynamics [1], charge and energy transfer phenomena in condensed phases [1,3], Kondo physics [2], and decoherence dynamics of superconducting qubits [2,4].
In such applications, the spin system can represent donor-acceptor charge states, a magnetic impurity [1], or a truncated harmonic spectrum, mimicking an anharmonic oscillator [5,6].
The bosonic bath may stand for a collection of lattice phonons, electromagnetic modes, bound electron-hole pairs, and other composite bosonic excitations [1,2].
From the theoretical perspective, the NESB model is an extremely rich platform for studying nonequilibrium quantum physics. One is interested in studying its transport characteristics, including transient dynamics and steady state properties, while covering different regimes: low-to-high temperatures, weak-to-strong system-bath coupling, adiabatic-tononadiabatic spin dynamics, with or without a (magnetic) spin biasing field, from linear response to the far from equilibrium regime. This challenge could be tackled by extending open quantum system methodologies, previously developed to treat the dissipative dynamics of the (traditional) SB model, to treat the more complicated, nonequilibrium, two-bath version.
Among the techniques developed to study the characteristics of the thermal heat current in the NESB model we recount perturbative quantum master equation tools: Redfield equation [5-7, 9, 25-28], the noninteracting blip approximation (NIBA) [6,7,28,29], as well as Keldysh nonequilibrium Green's function (NEGF) methods [30,31]. Computational studies had further established the non-monotonic behavior of the heat current with the spin-bath coupling energy, including studies based on the multi-layer multi-configuration Hartree approach [32], the iterative influence functional path integral technique (in the spin-fermion representation) [11,33], Monte Carlo simulations [10], and the hierarchical equation of motion [34,35].
Beyond the analysis of the thermal conductance or the energy current, in small systems the fluctuations of the current are expected to reveal plethora of information, such as current correlations to all orders [36,37]. In fact, rather than focusing on the thermal conductance, it is more demanding yet highly profitable to pursue the probability distribution P t (Q) of the transferred energy Q within a certain interval of time t. This measure is also known as full-counting statistics (FCS) in the context of electron transport. Obtaining the FCS for interacting systems is a highly desirable, yet formidable task. The FCS of the NESB model has been analyzed so far in two different limits: (i) in the sequential tunneling limit i.e., to the lowest order in the system-bath tunnelling strength, by employing the Redfield-type quantum master equation approach [27,29], and (ii) in the strong coupling and/or high temperature regime, following NIBA-type quantum master equations [29,38]. A theory interpolating these two limits was presented in Ref. [28]. However, these studies still miss the low temperature limit [31].
In this paper, we use the Schwinger-Keldysh NEGF approach [39,40] in combination with the Majorana fermion representation for the system-spin operators, and obtain the cumulant generating function (CGF) of the NESB model beyond the weak spin-bath coupling limit. The crucial impetus to introduce the Majorana representation is that in the fermion representation we are able to use Wick's theorem, thus obtain relevant nonequilibrium spinspin correlation functions-while including the counting parameter. Our results go beyond the sequential and co-tunnelling limits, and we are particularly able to capture to all orders the low temperature regime. As well, we observe deviations from the weak spin-bath coupling limit. On the other hand, while our result is valid beyond the strictly weak-coupling limit, it does miss the strong-coupling behavior as received in Refs. [10,11,28,32,41] using the NIBA approach. The main outcome of our study is an analytic expression for the FCS of the NESB model, capturing quantum effects, interactions, and far-from-equilibrium function.
The paper is organized as follows. We introduce the nonequilibrium spin-boson model and the Majorana fermion representation in Sec. II. In Sec. III, we present our main results for the CGF, followed by a discussion over different limits and numerical examples.
We further compare our expressions to previous theories on the NESB model, and to the harmonic oscillator-junction model. We conclude in Sec. V. The derivation of the CGF is explained in details in the Appendix.

II. MODEL
The NESB model comprises a two-state (spin) system coupled to two bosonic reservoirs (ν = L, R), which are maintained at different temperatures. The generic form of the full Hamiltonian is Here, σ i (i = x, y, z) are different components of the Pauli matrix, ω 0 and ∆ represents level detuning and the hopping between the spin states, respectively. b † j,ν (b j,ν ) is the creation (annihilation) operator of the j-th phonon mode in the ν-th reservoir. The last term describes the system-bath coupling term with λ j,ν as the coupling strength. For simplicity, we focus here on the unbiased case with degenerate spin levels (ω 0 = 0). Performing a unitary transformation, given by U = 1 √ 2 (σ x + σ z ), the transformed Hamiltonian reads We are interested here in obtaining the steady state energy current and its statistics beyond the weak system-bath coupling limit. Unlike the Redfield master equation technique, which captures only resonant energy transfer processes due to its underlying weak coupling approximation [6], the Keldysh nonequilibrium Green's function (NEGF) method offers a well established procedure so as to treat the system-bath interaction in a systematic-perturbative way [39,40]. However, the validity of Wick's theorem is a crucial requirement for practicing the method. Due to the lack of standard bosonic or fermionic commutation relations for spin operators, the NEGF approach is in fact unsuitable to be used in the spin representation of the NESB model. However this problem can be avoided by mapping the impurity spin to fermions, using the Majorana-fermion representation [31,42,43].
Explicitly, the spin operators can be expressed as σ = − i 2 η × η, i.e., Majorana fermions satisfy the anti commutation relation, {η α , η β } = 0, for α = β, η 2 α = 1, and unlike the Dirac fermions, they are real η α = η † α . Therefore, these fermions can be constructed in terms of ordinary Dirac fermions (f , g) and their conjugates as In this context, it is important to introduce the so-called copy-switching operator in terms of which the Majorana fermions can be expressed as σ α = τ x η α . Note that τ x commutes with all Majorana fermion operators and therefore is a constant of motion. Also, With the help of this operator, the spin-spin correlator reduces to correlator involving two Majorana fermions In this mixed Majorana-Dirac representation, the full Hamiltonian reads where B ν ≡ j λ j,ν (b j,ν +b † j,ν ) is a ν bath operator coupled to the spin system. Note that in this representation, the system-bath coupling term is no longer given in a bilinear form. For later use, we also identify the components of the Hamiltonian asH =

A. Working expressions for the FCS
The complete information over the energy transport statistics can be obtained from the so-called cumulant generating function, G(ξ), for heat exchange. We begin by defining the energy current operator as the rate of change of energy in one of the reservoirs, say L, and write down the heat current as dt . The operators are written in the Heisenberg picture, and they evolve with respect to the total HamiltonianH in Eq. (7). Therefore, the total energy change in the L solid within the time interval t 0 = 0 to t, where t 0 (t) is the initial (final) observation time, is given by the integrated current Following this definition, we write down the characteristic function Z(ξ) based on the twotime measurement protocol [36,37], Here, ξ is the "counting-field", keeping track of the net amount of energy transferred from the solid L to the spin. ... represents an average with respect to the total density matrix at the initial time, ρ T (0). We assume a factorized initial state, ρ with reservoirs prepared at a canonical state with inverse temperature β ν = T −1 ν , ρ ν (0) = e −βν Hν /Tr ν [e −βν Hν ], and an arbitrary state for the spin system ρ S (0). We also use the definition, for the counting field-dependent unitary evolution. Here, p = ±ξ/2 corresponds to the forward and backward evolution branches. Note that due to the measurement protocol, the modified HamiltonianH p acquires a phase in the system-bath coupling term, modifying only the left-bath operators, Here, is the a bath operator, dressed by the counting field. In the second line of Eq. (10), the operators are written in the interaction picture with respect to the non-interacting part of the Hamiltonian H S + H L + H R . T c is the contourordered operator which orders operators according to their contour time; earlier contour-time operators are placed to the right of later-time terms. In the long time limit, the CGF is defined as Here, Q n represent cumulants. Specifically, the second cumulant is Taking derivatives of the CGF with respect to ξ immediately hands over the current and its higher order fluctuations, or cumulants. However, instead of working with the CGF directly, one can manipulate the so-called generalized current, defined as by following the nonequilibrium version of Feynman-Hellman theorem first introduced by Gogolin et al. [44]-in the context of counting statistics for charge transport. The key advantage in treating the generalized current, rather than the CGF, lies in the fact that the problem can be treated with the diagrammatic NEGF technique, as developed originallywithout the counting field [45][46][47][48].
Using the NEGF with counting fields as developed in [49], an expression for the generalized energy current can be formally organized as When ξ = 0, this expression reduces to the standard Meir-Wingreen formula [50] for heat current [30]. The symbol tilde represents that operators within the Green's functions evolve with the dressed (counting field-dependent) HamiltonianH p given in Eq. (12).Π <,> xx (ω) are the Fourier transformed lesser and greater components of the spin-spin correlators, namely, Σ <,> ν (ω) are the self-energy components emerging due to the coupling of the spin to the solids, responsible for transferring energy in and out of the system, Here,n ν (ω) ≡ [1 + n ν (ω)] with n ν (ω) = (e βν ω − 1) −1 as the Bose-Einstein distribution function and β ν = 1/T ν is the inverse temperature. Γ ν (ω) = 2π j λ 2 j,ν δ(ω − ω j ) is the spectral function for the ν reservoir.

B. Main results
To receive the generalized current, our primary objective is to obtain the components Π </> xx (ω). These terms are obtained using the NEGF method following a first order perturbation expansion with respect to the interaction of the bath with the spin. We summarize here the central results; details are given in the Appendix.
The lesser and greater components are obtained to the lowest non-zero order in the nonlinear self-energy. They are given as with Here M(ω, ξ) = C 2 (ω) + 4 A(ω, ξ) includes the two terms, If we eliminate the counting parameter, ξ = 0, Π </> xx (ω) provides the imaginary components of the response function Π R xx (ω), matching the results of Ref. (31).
Using these expressions, the CGF for the NESB model, with the temperature-dependent transmission function This expression is valid with an arbitrary form for the spectral function Γ ν (ω). The CGF further satisfies the steady state Gallavotti-Cohen fluctuation symmetry, [51]. Eq. (23) constitutes the main result of our work.
The cumulants of the energy flux can be readily obtained by taking derivatives of the CGF with respect to the counting field ξ. For example, the heat current and its noise are given by The result for the current agrees with the derivation in Ref.
[31]-once we organize our expressions, In the next subsection, we discuss interesting limits of the general results.

C. Special limits
Incoherent sequential tunnelling. When the system-bath coupling is weak and the reservoirs' temperatures are high, Γ L,R ≪ ∆ ≤ T L,R , the above generating function reduces to the result obtained from the Redfield quantum master equation approach [29], when directly employing the Born-Markov approximation. We now derive this result. Following Eq. (23), the generalized current can be simplified to where M ′ (ω, ξ) = ∂M (ω,ξ) ∂(iξ) . To the lowest order O(Γ 2 L,R ), working in the limit Γ L,R ≪ ∆ ≤ T L,R , the poles in the integrand can be approximated by By employing the residue theorem, the integration in Eq. (27) and the generating function reduces to This expression matches the result obtained in Ref. [29]. This CGF also respects the fluctuation symmetry. It immediately yields the heat current in the weak coupling limit [5] I weak Co-tunnelling. At low temperatures, Γ ν ≪ T ν ≤ ∆, the process of sequential tunnelling is exponentially suppressed since incoming phonons are off-resonance-with frequencies below the spin energy gap, ω ≪ ∆. The dominant contribution to the current and higher order fluctuations thus comes from coherent two-phonon co-tunnelling processes. In this limit, the transmission function of Eq. (24) is given by T co with fluctuation symmetry being satisfied. Here, ω h , the upper limit in the integral should be determined by the smaller energy scale, temperature of the cutoff frequency of the baths.
The co-tunneling (co) heat current then becomes This expression was previously achieved in two ways: (i) By using a systematic perturbative treatment [25], and (ii) working with the so-called Born-Oppenheimer approach for heat exchange [52], by assuming slow bath and a fast (high frequency) impurity. In the case of an Ohmic bath, Γ ν (ω) ∝ ω s with s = 1, the heat current scales as I co SB ∝ T 4 L − T 4 R , thus the thermal conductance scales with T 3 , in agreement with numerically exact simulations on the NESB model [10]. As well, in this low temperature limit the NESB junction behaves similarly to a fully harmonic junction, as we discuss in Sec. III D.
Note that in contrast to the CGF received in Eq. (23) and Eq. (29), the CGF in the co-tunnelling limit is symmetric with respect to Γ L,R (ω). Therefore, in this limit the system does not support the thermal rectification effect. Moreover, in this limit the cumulants C n = ∂ n G SB (ξ) ∂(iξ) n ξ=0 scale as C n ∝ 1/∆ 2 , whereas in the sequential tunneling limit cumulants grow as C n ∝ ∆ n .

D. Comparison between the NESB model and the harmonic oscillator junction
In the harmonic oscillator (HO) junction, a single harmonic oscillator of frequency ω 0 , replaces the spin impurity of the NESB model, Eq. (1). The resulting Hamiltonian is fully harmonic, and it can be readily solved exactly to yield the CGF [53,54] Further, (ii) there is a crucial sign difference in this CGF as compared to G SB (ξ) in Eq. (23).
This sign difference reflects on the nonlinear nature of the spin. A similar sign-difference between harmonic and spin impurity nanojunctions has been observed in vibrationallyassisted electron conducting junctions [21,24]. The above expression immediately provides the Landauer expression for the heat current, and the noise In the weak coupling limit, the CGF of the HO model reduces to the standard result obtained by a low order QME [27,29] G weak HO (ξ) = with C HO (ω 0 ) = Γ L (ω 0 )+Γ R (ω 0 ), The heat current then reduces to the familiar result, The co-tunnelling limit is more subtle, and we exemplify it now when calculating the current.
We break the transmission function (34) into two contributions (leaving for a moment the numerator) T HO (ω) = T o (ω) + T e (ω), Assuming the hierarchy of energies Γ ν ≪ T ν ≤ ω 0 , we note that the function ω[n L (ω)−n R (ω)] changes slowly at the vicinity of ω 0 , in the regime where the functions T e,o (ω) have significant weight. Therefore, the integral (35) over the odd component (approximately) cancels out, and the current is solely determined by the even term, T e (ω) ∼ 1/ω 2 0 , to yield This result reproduces exactly the behavior of the NESB model in the corresponding limit, Eq. (32). This correspondence is not surprising: At low temperatures (smaller than the energy spacing in the quantum impurity) and at weak system-bath coupling, the NESB and the HO junctions should behave rather similarly. For a comprehensive analysis of the harmonic-mode thermal junction, see Ref. [55].

E. Steady state population and a bound on heat current
Besides transport properties, we use the Majorana formalism and calculate the steady state population of the ground and excited states in the eigenbasis of the spin. This can be obtained by calculating σ z , given as, The function C(ω) is defined in Eq. (21). In the weak coupling limit, we receive the same result as obtained in Ref. [6], .
The population of the states are p g = 1 2 (1 − σ z ) and p e = 1 2 (1 + σ z ). Recently, a rigorous quantum mechanical bound for the heat current in interacting systems has been derived, valid at the high temperature-yet in the quantum regime [56]. We now confirm that the heat current derived in our work, Eq. (25), does not violate the bound.
This further affirms the validity and usefulness of our result.
In the following analysis we make use of the inequality 0 ≤ [n L (ω) − n R (ω)] ≤ (T L − T R )/( ω) for ω > 0 and T L > T R . As well, we recall on the positivity of the transmission function T SB (ω) > 0. Furthermore, we assume an Ohmic spectral density function for the reservoirs, Γ ν (ω) = γ ν ω, ν = L, R (see Ref. [56] for a detailed discussion over different spectral functions). Putting these pieces together, we conclude that the heat current of Eq.
(25) satisfies the following inequality which precisely matches with the bound organized in Ref. [56] for the NESB model. We conclude that our expression for the current thus does not violate a fundamental bound, unlike the prediction of the Redfield QME, see Ref. [56].

IV. NUMERICAL RESULTS
In Figs. 1-3, we present simulations demonstrating the behavior of the heat current I and the second cumulant S , based on Eq. (23), as a function of the system-bath coupling, averaged temperature, and temperature difference. We focus on the following questions regarding the operation of the NESB nanojunction: (i) How are the current and noise influenced by the system-bath coupling strength? (Fig.   1 and 3). (ii) What are the signatures of operation far from equilibrium, as opposed to the linear response regime? (Fig. 1 and 3) (iii) What is the temperature dependence of the heat current? (Fig. 2) (iv) Thermal diode effect: Can we enhance this effect if we go beyond the weak spin-bath coupling? (3) (v) What is the relation between the Majorana-based treatment and other techniques? (Figs. 1-3). Fig. 1 displays the current and the noise as obtained from Eqs. (26), as well as the weak coupling (Redfield) limit [11,29], and the NIBA approximation [11,41]. We use an Ohmic spectral function for the baths with an exponential cutoff, Γ ν (ω) = πα ν ωe −ω/ωc . In accord with previous results (for the heat current [11]), we find that Redfield equation dramatically overestimates the current and the noise in comparison to the (more accurate) Majorana and NIBA results. Majorana treatment shows a saturation of the current and its noise at large α, while under NIBA these quantities quickly decay beyond α ∼ 0.15. Since the temperature is rather high, ∆ = T a , with T a = (T L + T R )/2, we expect the NIBA to be rather accurate here [10,11,41]. We also confirm in panel (a) that in linear response (LR), the conductance, I LR /∆T , is proportional to the thermal noise in the junction, in accord with the Green-Kubo relation, Far from equilibrium [see panel (b)], we obviously observe violations of the above relation.
However, it is interesting to note that the current and noise still follow a similar functional form within the three different methods. show up in anharmonic nanojunctions [12].
Next, we discuss the operation of the NESB as a heat diode, as suggested in Ref. [5].
To materialize this effect, it is necessary to (i) include anharmonic interactions, and (ii) introduce a spatial asymmetry [8]. The NESB model naturally includes an anharmonic potential. We break here the left-right symmetry by using different coupling strengths at the contacts, α L = α R . In Fig. 3, we analyze the ratio between the forward and backward currents as we switch the temperatures of the two baths, We set α L =0.01, 0.2, and modify α R over a broad range of values.
Based on Eq. (30), we can readily confirm that under the Redfield formalism the rectification ratio R does not depend on the absolute value of α (given the linearity of the current with α), only on the ratio α R /α L . In contrast, the Majorana treatment, which goes beyond weak coupling, reveals that the diode effect is enhanced as we increase the coupling strength itself. This result points out to the crucial role of many-body interactions in realizing the diode function.

V. CONCLUSIONS
We have studied the statistics of energy transfer in the nonequilibrium spin-boson model.

By combining Majorana fermion representation for the spin operators with the Schwinger-
Keldysh Green's function approach, we were able to derive an analytical expression for the CGF of the model. This function, which we confirmed here to satisfy the fluctuation symmetry for heat exchange, hands over the complete information over the energy statistics in the steady state limit. Our approach goes beyond the weak-coupling (Redfield) and the co- tunnelling limits. Surprisingly, the CGF of the NESB model has a similar structure as in the harmonic oscillator junction, besides sign differences and the appearance of a temperaturedependent transmission function-in the NESB model. These differences reflect on the nonlinear nature of the spin-boson system.
We have presented numerical examples for the heat current and its noise, and compared our results to previously-developed quantum master equation approaches, namely Redfield and the NIBA. We have further demonstrated that a heat diode becomes more effective as we increase the system-bath coupling. Additional improvements to the Majorana formulation presented here could be made, e.g., by developing a polaron-transformed Majorana fermion-NEGF approach [58]. Future work will be focused on simulating counting statistics in the NESB model beyond perturbative approaches [57]. Our goal is to evaluate the generalized current, Eq. (15). It is given in terms of the Keeping in mind the nonequilibrium setup, we introduce the ξ-dependent contourordered Green's function for the σ x component, Recall that · · · ξ means that operators are evolving with the dressed Hamiltonian of Eq. (12). Here τ, τ ′ are the contour times. When projecting to real time (t, t ′ ), we re-ceive four different terms, namely, time-ordered (t), anti-time ordered (t), lesser (<) and greater (>) Green's functions.
In Eq. (A6), G η (τ, τ ′ ) = −i T c η z (τ )η z (τ ′ ) is the Green's function involving the z-th component of the Majorana fermion,λ is the Nambu matrix andΣ L , Σ R are the bare Green's functions for the Bosonic baths, Recall that the operators of the left reservoirs are dressed by the additional ξ dependence, i.e., , when τ is on the upper (lower) branch. Given the perturbative nature of our treatment, the self-energy contribution from the baths is additive.