Quantum heat statistics with time-evolving matrix product operators

We present a numerically exact method to compute the full counting statistics of heat transfer in non-Markovian open quantum systems, which is based on the time-evolving matrix product operator (TEMPO) algorithm. This approach is applied to the paradigmatic spin-boson model in order to calculate the mean and fluctuations of the heat transferred to the environment during thermal equilibration. We show that system-reservoir correlations make a significant contribution to the heat statistics at low temperature and present a variational theory that quantitatively explains our numerical results. We also demonstrate a fluctuation-dissipation relation connecting the mean and variance of the heat distribution at high temperature. Our results reveal that system-bath interactions make a significant contribution to heat transfer even when the dynamics of the open system is effectively Markovian. The method presented here provides a flexible and general tool to predict the fluctuations of heat transfer in open quantum systems in non-perturbative regimes.


I. INTRODUCTION
The importance of heat management at the nanoscale has grown in tandem with advances in the fabrication and control of small devices, motivating increasing interest in the non-equilibrium thermodynamics of open quantum systems [1][2][3][4]. For example, quantum thermal machines have been studied in such diverse experimental platforms as single-electron transistors [5][6][7], trapped ions [8][9][10], superconducting circuits [11], and spin ensembles [12,13]. Numerous technologically or biologically important systems are also naturally described as quantum heat engines, including lasers [14], light-emitting diodes [15], and light-harvesting complexes [16][17][18][19]. These minuscule machines all operate far from equilibrium and are significantly affected by quantum and thermal noise. Strong coupling may blur the boundary between system and environment [20,21], potentially leading to non-Markovian effects [22,23] with interesting thermodynamic consequences [24][25][26][27][28]. In addition, the importance of fluctuations at small scales means that the statistical character of thermodynamic quantities such as work and heat cannot be ignored [29,30]. These features together give rise to a rich and varied phenomenology with important ramifications for emerging quantum technologies.
A crucial limiting factor for the performance of quantum devices is the transfer of heat to and from their surroundings. A detailed understanding of heat transfer is therefore essential to optimise control protocols while minimising wasteful dissipation [31][32][33]. More generally, heat flux is a fundamental source of irreversibility and entropy production in open quantum systems [34,35]. Entropy production limits the efficiency of heat engines and refrigerators [36], determines the energy cost of information erasure [37] and feedback control [38], constrains current fluctuations far from equilibrium [39][40][41][42][43], and can be directly measured in well controlled quantum settings [44][45][46]. However, modelling heat transfer in strongly coupled systems is a difficult theoretical problem because it requires access to the energetics of the bath.
On the contrary, the majority of techniques for describing open quantum systems either neglect the environment's dynamics completely or treat it via an effective or approximate description [47]. An accurate, tractable method to predict the fluctuations of heat transfer in generic open quantum systems is still lacking.
Here, we fill this gap by developing an efficient numerical method to compute heat statistics using the path-integral formulation of dissipative quantum mechanics [48]. Previous research has shown that the probability distributions of heat and work can be formally derived within this framework [49][50][51]. However, a direct evaluation of the corresponding path integral is only possible for a few exactly solvable models, while numerical approaches based on the quasi-adiabatic path integral (QUAPI) method [52,53] require careful finetuning to avoid error accumulation [54,55]. We solve this problem by generalising the TEMPO algorithm [56] to calculate the characteristic function of energy changes in the bath, equivalent to the Fourier transform of the heat probability distribution. This algorithm exploits a tensor-network representation of the QUAPI propagator to describe complicated non-Markovian evolutions efficiently [57]. As a result, we obtain a flexible and accurate tool to describe fluctuating heat transfer in generic, strongly coupled open quantum systems, which can be extended to deal with time-dependent Hamiltonians [58] or multiple baths [55].
The canonical open quantum system comprises a small, few-state system coupled to a bosonic bath. This general setting is known to be amenable to efficient tensornetwork descriptions [59,60]. For the sake of concreteness, in this work we focus on the paradigmatic spinboson model, which describes quantum dots [61], ultracold atomic impurities [62] and superconducting circuits [63], to name just a few examples. We demonstrate our approach by applying it to the non-equilibrium quantum thermodynamics of this important model. We first verify the accuracy of our method by comparison with the exact solution in the limit of the independent boson arXiv:2008.06491v3 [quant-ph] 10 Jun 2021 model. Then we compute the time-dependent heat transfer and its fluctuations across a range of parameters in the unbiased spin-boson model, including the challenging low-temperature and strong-coupling regimes. We interpret our results using the notion of generalised equilibration in strong-coupling thermodynamics [21], and develop analytical models that quantitatively explain the mean heat exchange in the high-temperature and lowtemperature limits. We also show numerically that the heat distribution obeys a fluctuation-dissipation relation (FDR) in the high-temperature limit, which is similar to the well-known FDR of the work distribution [64]. Interestingly, our results show that the system-bath interaction energy makes a considerable contribution to the heat statistics, even in the weak-coupling and hightemperature regime where a Markovian description of the system dynamics alone is accurate. This underlines the need to interpret with great care the standard Markovian description of quantum thermodynamics [65], which is based on properties of the open system alone.
A brief outline of the paper is as follows. In the next section, we introduce the spin-boson model and define the thermodynamic quantities of interest. Details of our numerical method are provided in Sec. III. We then present results for the independent-boson and spin-boson models in Sec. IV, before concluding in Sec. V. Units where = 1 = k B are used throughout.

A. Quantum thermodynamics of relaxation processes
We are interested in the non-equilibrium thermodynamics of an open quantum system coupled to a large heat bath. The Hamiltonian of such a system can be written generically aŝ whereĤ S is the free Hamiltonian of the quantum system,Ĥ B is the free Hamiltonian of the environment, and H I is the Hamiltonian that describes the interaction between these two components. Following the standard approach [66], the bath is modelled by an infinite collection of harmonic oscillators coupled linearly to the system, so thatĤ Here,â j is the annihilation operator for mode j of the bath, ω j is the corresponding mode frequency, g j is a coupling constant, andŜ is an arbitrary system operator.
The bath is characterized by its spectral density function [52] Non-equilibrium processes at the nanoscale feature significant and measurable fluctuations. Therefore, thermodynamic process variables such as work, W , and heat, Q, must be promoted to stochastic quantities described by the corresponding probability distributions, P (W ) and P (Q). Thermodynamic work is associated with changes in the external conditions defining the Hamiltonian, while heat is defined here to be the change in energy of the bath. Operationally, each of these quantities can be extracted from a two-point measurement ofĤ (work) orĤ B (heat) at the beginning and end of the evolution, either with direct projective measurements [67] or via ancillary probes [68][69][70][71]. Therefore, under strong-coupling conditions where the commutator [Ĥ B ,Ĥ I ] is non-negligible, work and heat are simultaneously measurable only if the system-bath interaction vanishes at the beginning and end of the evolution [21]. This is the relevant scenario for cyclic thermal machines, for example, and also the one that we assume here.
Following the above reasoning, we consider the relaxation dynamics starting from a product statê is the reduced state of the open quantum system andρ B (t) is the state of the bath with a thermal initial conditionρ B (0) = e −βĤ B /Tr[e −βĤ B ] at inverse temperature β = 1/T . The system evolves in time according toρ(t) =Û (t)ρ(0)Û † (t), whereÛ (t) = e −iĤt is the time evolution operator. The energy and entropy change of the system are given respectively by where we denote time-dependent expectation values by is the von Neumann entropy of the stateρ. Note that, unless the initial and final states of the system are in thermal equilibrium, neither ∆U nor ∆S as defined above necessarily correspond to variations of thermodynamic potentials. Since the Hamiltonian is time-independent during the relaxation process, all energy transferred during the evolution is in the form of heat exchanged with the bath. The mean heat absorbed by the bath is given by The first law of thermodynamics states that where W is the average work performed on the entire system by switching the system-bath interaction on and off at the endpoints of the evolution. Assuming that this switching is instantaneous, we have that W = − Ĥ I t (i.e., Ĥ I t is the mean interaction energy just before it is switched off), which follows from energy conservation, Ĥ t = Ĥ 0 , and the fact that Ĥ I 0 = 0 for an interaction of the form of Eq. (41). The average heat dissipated into the bath therefore comprises two contributions: the change in the system's internal energy and the system-bath interaction energy developed throughout the relaxation process. This dissipation is associated with an average entropy production which obeys Σ ≥ 0 in accordance with the second law [34,35], where equality holds for reversible processes.

B. Heat statistics
By definition, the heat transfer is the energy change that would be registered by projective energy measurements on the bath at the beginning and end of the process. We denote byΠ n = |E n E n | the projector onto the eigenstate |E n ofĤ B with eigenvalue E n . The heat distribution is then defined by where p n = Tr[(1⊗Π n )ρ(0)] is the probability of measuring initial energy E n , and p n→m = Tr[Π mÛ (t)(ρ S (0) ⊗ Π n )Û † (t)] is the conditional probability for the transition E n → E m [49]. The fluctuating heat exchange can be characterised by the statistical moments Here, we have introduced the characteristic function where u is the counting field parameter. Using Eq. (11), one easily obtains [29] It is convenient to define a modified time evolution operator asV This allows us to rewrite Eq. (15) as χ(u) = Tr [ρ(t, u)], with the modified density matrix Definingρ S (t, u) = Tr B [ρ(t, u)] as the reduced modified system density matrix, we have The form in Eq. (18) facilitates the calculation of the heat statistics by means of path-integral techniques.

III. PATH INTEGRAL METHODS
A. Influence functional for the modified density matrix The dynamics of the modified reduced density matrixρ S (t, u) can be formulated as a path integral [48], in which the effects of the environment on the open quantum system are captured by an influence functional that is non-local in time. This Feynman-Vernon influence functional is well suited to numerically discretised approaches such as QUAPI [52,53], upon which the TEMPO algorithm is built [56]. Here we describe how to obtain the influence functional modified by the counting field u.
To derive a discretised path integral for the modified density matrix,ρ The path integral forρ S (t, u) is constructed by inserting resolutions of the identity in the eigenbasis of the system coordinateŜ at each time step and then tracing over the bath. We use the notation s ± k for the eigenstates ofŜ, where the index k indicates the time t k = k∆ and the superscript + (−) is used to label eigenvectors inserted on the left (right) of the density matrix. Given our product initial condition in Eq. (5), we find is a product of free propagators for the system, with where ∆ k = ∆ for k < N and ∆ N = ∆/2, while the modified influence functional is Above, we have definedĤ env u (s) = s|Ĥ env u |s as the environment Hamiltonian conditioned on a particular eigenvalue s of the system coordinate.
To evaluate the influence functional explicitly for the spin-boson model, we introduce a compact superoperator notation [51,54]. Let us move to an interaction picture with respect to the free Hamiltonian From Eq. (17), we derive the differential equation withH The solution for the modified reduced density matrix isρ where the influence superoperator is given by and we introduced the time-ordering symbol ← − T , which reorders superoperators such that time increases from right to left, and the reservoir average, for any superoperator X , Since the interaction HamiltonianĤ I is linear and the reservoir thermal state is Gaussian, we may express Eq. (25) exactly using a time-ordered cumulant expansion up to second order [72]: The exponent is evaluated using well-known properties of bosonic thermal states; see Appendix A for details. The result is expressed in terms of three correlation functions: where J(ω) is the spectral density function of the bath defined in Eq. (4). Following Ref. [54], we recover the path integral representation from Eq. (24) by simply discretising time into N intervals and inserting resolutions of the identity at each time step,1 = s ± k s ± k s ± k , as in Eq. (20). In the interaction picture, the free propagators F ({s ± k }) do not appear and we obtain the influence functional in the form Here, ∆k = k − k and η qq k−k (u) are the discretised correlation functions Our expression for I s ± k , u matches the one recently derived in Ref. [55] and it is straightforward to verify that, for u = 0, it reduces to the original influence functional described in Ref. [52].
The form of Eq. (31) emphasises that the environment introduces memory into the evolution by coupling the system coordinate to itself at different times. Crucially, however, the correlation functions η α (t, u) decay to zero for sufficiently large t and therefore the memory time of the environment is finite. This insight forms the basis of the TEMPO algorithm described in the following section.
B. TEMPO algorithm TEMPO [56] is an efficient algorithm to compute path sums of the form of Eq. (20), given an influence functional of the form of Eq. (31). The standard TEMPO algorithm can be applied directly to our problem, with the only novelty being that here the influence functional is parametrised by the counting field u. We therefore provide only a brief summary of TEMPO here, directing the interested reader to Ref. [56] for a detailed description.
The key assumption of both the QUAPI and TEMPO methods is that the non-local time correlations encoded in the influence functional have a finite range, i.e. η α (t, u) ≈ 0 for t > τ C , where τ C is the bath memory time. Therefore, in the discretised form of the modified influence functional (31) one can introduce a maximum value of |k − k | beyond which the coefficients η α k−k (u) are negligible for all u. As a result, we may approximate The assumption of finite memory depth allows for an efficient description of the quantum dynamics through an iterative tensor propagation scheme, which forms the basis of QUAPI [52]. To see this, note that the summand in Eq. (20) can be viewed as an (N + 1)-index object called the augmented density tensor (ADT), denoted with d the dimension of the system S). The modified density matrix is found by summing over all but the final index, i.e.
where the remaining index σ N = {s + N , s − N } is determined by the values of s ± N on the left-hand side. The ADT is built iteratively starting from the initial condition with δ σ µ the Kronecker delta symbol, the ADT at the nth time step is given by the contraction with the Einstein summation convention assumed. Due to the finite memory depth K, the propagator (37) acts non-trivially on at most K indices of the ADT, since I ∆k (s ± n , s ± n−∆k ) = 1 for ∆k > K. At the nth time step, therefore, when n > K one needs only to store the object A σn···σ n−K , with the remaining indices summed over. (For the first K time steps one stores the full ADT.) The limiting factor for QUAPI is the computational resources needed to store and perform contractions on Kindex tensors. The TEMPO approach circumvents this limitation by representing the ADT and the propagators as tensor networks, which can be stored efficiently using truncated singular-value decompositions, enabling very large values of K to be reached. Comparing with previously developed methods, which are able to perform at values up to K ∼ 10 [55], in the calculations we show in this paper the TEMPO approach performs at a value of K = 500. The tensor-network representation is efficient due to the finite range of temporal correlations contained in the ADT. This is analogous to the well-known ability of tensor networks to represent many-body quantum states exhibiting short-ranged spatial correlations [73]. In the present case, the bond dimension, i.e. the number of singular values retained during the construction of the tensor network, quantifies correlations between different time points induced by the non-Markovian environment. The bond dimension is controlled by retaining only those singular values λ greater than a cutoff λ C . We define the cutoff as λ C = λ max 10 −p/10 , with λ max the highest singular value. The accuracy of the algorithm is therefore controlled by the exponent p as well as the memory depth K and the numerical time step ∆.

IV. SPIN-BOSON MODEL RESULTS
Although our method is general, in the following we specialise to the spin-boson model describing a single spin one-half interacting with a bosonic bath of harmonic oscillators [66]. In this case, the terms in Eq. (1) take the formĤ Above,Ŝ z andŜ x are the spin operators for the system. We focus on an Ohmic spectral density function of the form where α is a dimensionless coupling constant and ω C is a large cutoff frequency. In the following, we consider two different limits of the spin-boson model: the independent boson model with Ω = 0, and the unbiased spin-boson model with ω 0 = 0 and Ω = 0. The independent boson model is exactly solvable, allowing us to verify the accuracy of our numerical method. We then turn to the unbiased spin-boson model, an archetypal example of a non-integrable open quantum system.

A. Independent boson model
The independent boson (IB) model is described by Eqs. (39)- (41) with Ω = 0. The Hamiltonian can be diagonalised by a polaron transformation, which takes the general form This describes a spin-dependent displacement of each bath oscillator by an amount proportional to f j . The choice f j = g j diagonalises the IB Hamiltonian aŝ P †ĤP =Ĥ 0 − 1 2 E r , whereĤ 0 =Ĥ S +Ĥ B is the free Hamiltonian and we have defined the reorganisation energy which determines the shift in ground-state energy due to the system-bath interaction. In the IB model, [Ĥ,Ĥ S ] = 0, meaning that the local energy of the spin is conserved and ∆U = 0. Therefore, the heat dissipated into the bath is associated purely with the system-bath interaction, as discussed in Sec. II A. In particular, we show in Appendix B that the heat characteristic function is independent of the state of the spin and given explicitly by Differentiation of this quantity yields closed-form expressions for arbitrary cumulants of the heat distribution, given in Appendix B. Specifically, the mean heat is found to be which is strictly positive and independent of temperature. Interestingly, these properties are shared by all odd cumulants of the heat distribution in the IB model. For an Ohmic spectral density, we have Q = αω 3 C t 2 /(1 + ω 2 C t 2 ), which monotonically approaches the reorganisation energy in the long-time limit: For an Ohmic spectral density function, Eq. (46) depends on only two parameters, the coupling strength and the frequency cutoff. While ω C sets the timescale of the heat transfer process, the mean exchanged heat scales linearly with α. At first glance, it is not obvious that for strong coupling our method will be able to give the correct prediction, as this regime is in general difficult to model. It is therefore of interest to demonstrate the validity of the numerical method for different values of α.
The mean heat is plotted as a function of evolution time for several different coupling strengths in Fig. 1. We use these results to validate the numerical algorithm, whose results are shown in the same plot. We find excellent agreement between our simulations and the exact solution for each value of α considered. A simple estimate of the accuracy of our approach is obtained by comparing the asymptotic heat values to the exact result in Eq. (47). For the convergence parameters we have used, we find a relative discrepancy of δQ/Q = 0.04% in the case of α = 0.1, which increases to δQ/Q = 0.67% in the case of α = 1.5. These discrepancies could be further reduced by increasing the accuracy of TEMPO through changing the convergence parameters ∆, p and K. For an in depth discussion on the accuracy of the mean heat calculations with respect to the convergence parameters and value of the counting field, see Appendix C.
To quantify the fluctuations of the exchanged heat, we consider the variance Q 2 = Q 2 − Q 2 , which is given by Unlike the mean heat in Eq. (46), which is independent of temperature, the variance in Eq. (48) depends on the inverse temperature of the bath β. We show that our method is accurate for both a lower and a comparable temperature k B T with respect to the energy scale of the system ω 0 . Fig. 2 shows the variance as a function of time, for different values of temperature and coupling strength. The numerical predictions again match the analytical solutions given by Eq. (48). Note that in order to get a better match between the solutions for high coupling, α = 1.5, the value of the counting field at which the numerical derivative of χ(u) is taken has been set to u = 0.005, compared to the value of u = 0.01 in the case of α = 0.1. This suggests that high coupling strength cases require in general more computational precision than low coupling cases, although not higher precision in the singular-value decomposition cutoff or time-step. The relative discrepancy in the asymptotic values between analytical and numerical solutions for Q 2 in the case of T = 1 are found to be δQ 2 /Q 2 = 0.12% for α = 0.1, and δQ 2 /Q 2 = 0.06% for α = 1.5. In the case of T = 0.1, α = 0.1, the relative discrepancy is δQ 2 /Q 2 = 0.13%.

B. Unbiased spin boson model
We now turn to the spin-boson model with Ω = 0, focussing on the unbiased case where ω 0 = 0. In this context, TEMPO has previously been used to pinpoint the localisation phase transition [56], which occurs when T = 0 and at a critical value of the coupling α [74,75], and to study non-Markovian dynamics induced by a spatially correlated environments [57]. Here we use it to investigate the non-equilibrium thermodynamics of relaxation over a range of temperatures and coupling strengths. In the following, we take Ω = 1, which defines our unit of energy.

High temperature and weak coupling
We begin by studying the regime of weak coupling and relatively high temperature, with α = 0.1 and T = 5. The mean heat transfer is plotted in Fig. 3 as a function of time, starting from a pure initial state,ρ S (0) = |Ψ 0 Ψ 0 |. Specifically, we consider three different initial conditions: We also consider two values of the cutoff, ω C = 5 and ω C = 50. Inspection of these results suggests that the heat transfer, Q is a sum of two contributions. The first contribution is the heat transferred directly from the system as it relaxes to a thermal stateρ eq S ∝ e −βĤ S . The corresponding change in internal energy will be The second contribution to the mean heat transfer is associated with switching on the system-bath interaction, and is equivalent to the work done in a cyclic process as discussed in Sec. II A. If we assume that this contribution is the reorganisation energy, as in the independent boson model, we expect This approximation shows near-perfect agreement with the long-time limit of the numerical results, as demonstrated by the dashed lines in Fig. 3. Notice that Eq. (49) is independent of the details of the bath spectral density (i.e. α and ω c ), while E r does not depend in any way on the spin degrees of freedom. This indicates that, at high temperature and weak coupling, the displacement of the bath modes is not affected by the thermalisation of the spin. Instead, these two processes give rise to independent and additive contributions to the mean heat transfer. These distinct modes of heat transfer take place on different time scales. This is illustrated by the blue lines in both the ω C = 5 and ω C = 50 case of Fig. 3, corresponding to the low-energy initial state |Ψ 0 = |← . First, heat is transferred to the environment as the systembath interaction forces the bath modes to rapidly adjust to their new equilibrium. This takes place over a time set by the inverse cutoff, ω −1 C ≈ 0.2 for ω C = 5 and ω −1 C ≈ 0.02 for ω C = 50. Then, the direction of heat flow reverses as the bath gives up energy in order to bring the spin to thermal equilibrium, which occurs on a slower timescale fixed by the inverse of the thermalisation rate, which can be estimated as γ ≈ (π/4)J(Ω) coth(βΩ/2) from standard weak-coupling theories, e.g. the secular Born-Markov master equation [76], giving γ −1 ≈ 0.8 for ω C = 5 and γ −1 ≈ 0.65 for ω C = 50. A comparison between the two different values of ω C in Fig. 3 shows how a larger frequency cut-off determines a shorter timescale for the heat transfer process, for fixed T and α. (E r is ten times larger in the ω C = 50 case, so that the energy due to the displacement of the bath modes dominantes over that due to the spin thermalisation.) It is worth emphasising that the system-bath interaction energy gives a significant contribution to the heat transfer, even though the system dynamics is very well captured by a Markovian, weak-coupling description. Indeed, for the parameters considered in Fig. 3 and ω C = 5, the reorganisation energy is comparable to the natural energy scale of the spin, since E r = Ω/2. Nevertheless, Fig. 4 shows that in this regime the calculated spin dynamics (solid curves) matches the corresponding Born-Markov and weak-coupling approximated problem (dashdotted curves), within the limits of such an approximation, the coupling strength being set to α = 0.1. The discrepancy shown in Fig. 4 is 10%.

Lower temperature and stronger coupling
We now consider the heat transfer at intermediate and low temperatures. In Fig. 5 we show the mean heat transfer for temperatures T = 1 and T = 0.1, starting from the state |Ψ 0 = |↑ . We see the same monotonic relaxation behaviour as was observed at high temperature (cf. the orange curve in Fig. 3), albeit proceeding on a slower timescale as the temperature is reduced. Outside of the high-temperature limit, the asymptotic value of Q can no longer be well approximated by Eq. (50), shown by the dashed lines in Fig. 5. We find that the the spin's internal energy change and the total heat transfer are smaller in magnitude than Eqs. (49) and (50)  as Fig. 5 and Fig. 6 both show. This demonstrates that the tendency of the spin to minimise its local free energy defined byĤ S competes with the displacing effect ofĤ I on the bath modes. As a consequence of this interplay, both ∆U and W depend non-trivially on system-bath correlations generated during the relaxation process. The effect of the correlations with the bath is indeed to decrease the magnitude of ∆U with respect to the value − Ω 2 tanh Ω 2T predicted by Eq. (49), and represented in Fig. 6 by the dashed lines. Such discrepancy is starkly greater for stronger coupling.
In order to understand this, we note that at strong system-bath coupling the equilibrium state must be generalised to [21] i.e. the reduction of a global thermal state. This takes into account correlations with the bath and reduces to the standard formρ eq S ∝ e −βĤ S in the weak-coupling limit. Assuming that the open quantum system couples to the bath locally in space, the interaction Hamiltonian is a local degree of freedom that is also expected to thermalise, in the sense that We emphasise that these thermalisation conditions hold only for local subsystems: they do not imply that the system as a whole attains thermal equilibrium in the longtime limit.
We estimate the effect of system-bath correlations on heat transfer using the variational approach pioneered by Silbey and Harris [77], which has been successfully applied to understand various static and dynamic properties of the spin-boson model [78][79][80]. The method is briefly summarised here with further details given in Appendix D. The basic idea is to express the Hamiltonian in a different basis by applying a unitary transformation that mixes the system and bath degrees of freedom. A judicious choice of transformation -determined in this case by a variational principle -leads to a weak effective interaction termĤ I in the new basis, even though the bare interactionĤ I may be strong.
Specifically, the Hamiltonian is diagonalised approximately using the polaron transformation in Eq. (43), with the displacements {f j } interpreted as variational parameters. After the transformation, the Hamiltonian is written asP †ĤP =Ĥ 0 +Ĥ I ≈Ĥ 0 . Here,Ĥ 0 is the free Hamiltonian in the variational frame, which is given up to a constant byĤ where Ω is a renormalised tunnelling matrix element to be defined below. The neglected interaction term,Ĥ I , describes residual transitions between dressed states of the system and environment and is proportional to the spin tunnelling amplitude Ω. The effect ofĤ I is made as small as possible by choosing the variational parameters to minimise the Feynman-Bogoliubov upper bound on the free energy; see Appendix D for details. This is achieved by taking f j = g j φ(ω j ) with which must be solved self-consistently for Ω . The heat transfer is then found by approximating e −βĤ ≈ P e −βĤ 0P † in Eqs. (51) and (52), yielding This has the same form as Eq. (50) but with both the tunnelling matrix element Ω and reorganisation energy E r = 1 2 dωJ(ω)φ(ω)/ω renormalised. The variational theory predicts that both the spin tunnelling matrix element and the reorganisation energy are reduced relative to their bare values, since Ω /Ω ≤ 1 and φ(ω j ) ≤ 1. Physically, this occurs because the tunnelling between spin states |↑ ↔ |↓ induced byĤ S is suppressed by the spin-dependent mode displacements generated byĤ I , which reduce the effective overlap between the two spin states. The equilibrium state emerges from a balance of these two competing effects, which explains why both ∆U and Q are reduced at low temperature relative to Eqs. (49) and (50).
We show in Fig. 7 that the variational theory gives a good quantitative approximation to the mean heat transfer at low temperature, T = 0.1, with the best agreement at weak coupling. At higher temperatures on the order of T = 1 and above, we find that the approximation breaks down completely because the renormalisation of the tunnelling amplitude is overestimated, leading to values Ω Ω. This failure is presumably due to the neglect of thermally activated transitions generated byĤ I , which become relevant at temperatures βΩ 1. On the other hand Fig. 7 shows that the additive ansatz given by Eq. (50) performs worse than the variational theory across all the coupling range.
At very strong coupling, the variational theory performs well at all temperatures. In this regime, strong correlations with the bath lead to an almost maximally mixed equilibrium state of the spin, corresponding to a vanishing tunnelling rate in the variational frame, Ω → 0. As a result, the heat transfer for an initial state |Ψ 0 = |↑ reduces to the bare reorganisation energy, E r . This behaviour is shown in the inset of Fig. 5, where the solid curve converges to Q ≈ E r to a good approximation. The dynamics of the heat transfer is correspondingly fast in this regime since it depends only on the bath cutoff scale, ω C .

C. Heat fluctuation-dissipation relation in the spin boson model
As a final demonstration of our method, we study the temperature dependence of the heat fluctuations in the spin-boson model. Fig. 8 shows the asymptotic variance of the heat distribution at long times, starting from the initial state |ψ(0) = |↑ . We see that the fluctuations increase with temperature, and grow approximately linearly with T at high temperature. This linear behaviour of Q 2 ∞ can be understood as a manifestation of the fluctuation-dissipation relation (FDR) that is well known in the context of nonequilibrium work distributions. If the distribution of work W is Gaussian, the Jarzynski equality directly implies the FDR [64] W − ∆F = β W 2 /2, where ∆F is the equilibrium free energy change. In the case of the independent-boson model, heat is identical to work sincê H S is a conserved quantity, while ∆F = 0 because the process is cyclic. It follows that we can write an equivalent FDR for the heat distribution: At high temperature, this relation holds at all times in the independent-boson model, as can be seen by comparing Eqs. (46) and (48) in the limit βω C 1. In the spin-boson model, we no longer have equality between work and heat since ∆U = 0. Nevertheless, we find numerically that the FDR (57) approximately holds at high temperatures, βω C 1, as shown in Fig. 9. This behaviour stems from the fact that the spin's contribution to the heat fluctuations is limited by its finite energy splitting Ω, whereas the contribution of the spin-boson interaction energy can grow arbitrarily large. The heat fluctuations are thus dominated by independent-boson physics at high temperature. For strong coupling, where the spin energy scale Ω is negligible compared to the reorganisation energy E r , we show in Fig. 9 that the heat fluctuations are essentially identical in the spin-boson and independent-boson models at all temperatures. One limitation encountered in the calculation of the data shown in Fig. 9, is that the TEMPO algorithm was not able to compute the variance up to equilibrium time for very low temperatures. Indeed, the lowest temperature shown is T = 0.4. Exploring the validity of the heat FDR in other scenarios, e.g. multipartite open quantum systems, is an interesting avenue for future work.

V. CONCLUSIONS
A better understanding of dissipation in open quantum systems is a fundamental goal of quantum thermodynamics as well as being crucial for quantum device engineering. We have shown that this goal can be successfully addressed by an extension of the TEMPO algorithm [81] to evaluate the characteristic function of the heat distribution. We have demonstrated the validity and flexibility of our approach by calculating the mean and variance of the heat transfer in the spin-boson model over a range of temperatures and system-bath coupling strengths. Our results clearly demonstrate the importance of systemenvironment correlations at low temperatures. Even at high temperature and weak coupling, we find significant contributions to the heat statistics from the systemenvironment interaction energy that are not captured by the standard weak-coupling master equation. This indicates that system-reservoir interactions are an important source of dissipation that must be accounted for when designing thermodynamic protocols [82][83][84][85][86][87], even in the weak-coupling regime.
Our approach to calculating heat statistics can be extended in several promising directions. It is straightforward to adapt the method to situations with a timedependent system Hamiltonian, which would enable the characterisation of heat statistics for driven open systems. This problem, which is theoretically challenging even for Markovian environments outside of the slowdriving regime, has numerous applications in quantum control, such as quantum information processing [88] and erasure [89], enhanced engine cycles through thermodynamic shortcuts [33,90], and tailored quantum light sources [31,91]. It is also possible to incorporate multiple baths within our framework by combining the corresponding influence functionals together. This would allow the study of the full counting statistics of quantum heat transport in non-equilibrium steady states [55], including highly non-Markovian regimes. In general, we expect that the method presented here will facilitate further research into the non-equilibrium quantum thermodynamics of strongly coupled open systems.
where we have used Eq. (A3) and Eq. (A4). The exponent of the modified influence functional in Eq. (27) can then be written as where the tilde denotes an operator in the interaction picture with respect toĤ 0 , i.e.
Using the Baker-Campbell-Hausdorff formula, e A e B = exp(A + B + 1 2 [A, B] + . . .), and neglecting an irrelevant phase factor, we obtainÛ (t) =Û 0 (t)Û I (t), whereÛ 0 (t) = e −iĤ0t is the free propagator and is the interaction-picture propagator, which describes a spin-dependent displacement for each mode of magnitude We find that the mean heat depends only on the imaginary part of the characteristic function and is given by the linear slope of the function depicted in Fig. 10 in an interval [0, u ], with an error of the order O (u ). The value u must be such that within the interval it defines, the real part of the characteristic function can still be approximated by a constant function, and the slope of the imaginary part is linear. u will depend on the model, as shown by comparing the two figures in Fig. 10, and on the physical parameters α, T and ω C . Indeed, Fig. 2 shows for example that while for α = 0.1 it is sufficient to take u = 0.01, for stronger coupling such as α = 1.5 it is necessary to set u = 0.005 to achieve the same precision.
We have found that in order to achieve a function Q that is constant in the long time limit, for the parameters considered in this work the value of u can't be greater than u = 0.01. In general, decreasing the value of u below u = 0.005 will increase the computational time but not improve significantly the precision of the result.

TEMPO memory depth
In Section III B we have discussed the finite memory depth K of the TEMPO algorithm that allows it to efficiently propagate the ADT. In this subsection of the Appendix we will show how the memory depth affects the convergence of the mean heat and the variance of the heat distribution studied throughout this work. The form of the correlation functions in Eq. (28)- (30) sets the minimum value of K∆ needed. Indeed, K∆ has to be large enough to so that the discretised correlation functions are zero. Preliminary calculations have shown that, for the values of temperature and coupling strength considered, this requirement is satisfied around the value K∆ = 5. Fig. 11 shows that in the independent-boson model, for a fixed value of ∆, both the mean heat and the variance of the heat distribution reach the predicted asymptotic value for K > 300, for all the values of T and α depicted. For values K < 100, however, the asymptotic TEMPO result diverges greatly from the predicted one. This clearly shows how our method, which is able to operate at high values of the memory depth, has a much greater accuracy than other methods which operate in the region K < 100.