Functional Integral approach to time-dependent heat exchange in open quantum systems: general method and applications

We establish the path integral approach for the time-dependent heat exchange of an externally driven quantum system coupled to a thermal reservoir. We derive the relevant influence functional and present an exact formal expression for the moment generating functional which carries all statistical properties of the heat exchange process for general linear dissipation. The general method is applied to the time-dependent average heat transfer in the dissipative two-state system. We show that the heat can be written as a convolution integral which involves the population and coherence correlation functions of the two-state system and additional correlations due to a polarization of the reservoir. The corresponding expression can be solved in the weak-damping limit both for white noise and for quantum mechanical coloured noise. The implications of pure quantum effects are discussed. Altogether a complete description of the dynamics of the average heat transfer ranging from the classical regime down to zero temperature is achieved.


Introduction
In recent years enormous advancements have been accomplished in the fast growing field of quantum thermodynamics [1,2,3,4]. This area of research is spurred both for fundamental and practical reasons. On the fundamental side, one of the main question is how classical thermodynamics emerges from quantum physics [5,6]. On the practical side, there is an urgent need for a thorough understanding: how quantum systems exchange energy and heat as a function of time? This is a prerequisite for building quantum devices and/or quantum heat engines with potentially enormous technological impact [6,7]. Regardless of the new developments, quantum thermodynamics can still be considered as a "young" field of research. This is because there is an intrinsic difficulty in considering the concept of classical thermodynamics at quantum level since one has to face with many problems. A striking example is the attempt to generalize the definition of work and heat done and dissipated in a quantum system, which is still heavily debated [1,8]. The most common approach is to measure the energy of the system twice and then iterate the protocol to reconstruct statistical information [9,10,11,1,3,12,13,14]. Despite its simplicity and usefulness, this approach is limited by the fact that a standard projective measurement would lead to wave function collapse and thus would seriously perturb the system dynamics. Alternative approaches have been proposed but the main issues are still open [15,16]. A complementary approach is to measure the energy, or heat, dissipated because of interactions between the system and an environment. In this case, it is essential to study the degrees of freedom of the environment itself. In this direction, a promising experimental proposal was put forward in Ref. [17,18], where the setup consists in an engineered resistor acting as an environment with a sufficiently small heat capacity, allowing the measurement of the change in temperature when a quantum of energy is emitted or absorbed. A second, and more technical, problem lies in the approximations generally made in describing the dynamics of a quantum system which interacts with an environment. In fact a rigorous microscopic description of energy exchange and dissipation is still lacking. In the following, we assume the standard definition of heat [11,1,3,19,20,21] and focus on the second issue by presenting a microscopic derivation of the process of heat exchange of a quantum system with the surroundings. We will approach the issue with functional integral techniques. In such a formalism, all quantum effects of the environment on a microscopic system, which are related to energy transfer, can be included in a specific influence functional of the coordinates of the system only. Within the requirements of a linear environment all conceivable quantum noise sources may be accounted for by appropriate spectral properties of the influence kernels. The central problem of this study is to develop a general formalism for finding all statistical features of the energy transfer between the system and the environment. In order to investigate the properties and advantages of this approach, and to map out potential of the method we will apply it to a specific quantum systems. As in the context of dissipative quantum dynamics [22,23], the influence functional does not depend on the particular quantum system under consideration. In addition, this formalism allows the study of the full dynamics of the system and supports insights into the physical processes of emission, absorption, and dissipation of energy. Moreover, one can take advantage of many numerical and non-perturbative methods [24,25,26,27,28], as well as different well-established analytical techniques [29,30,22,28,31,32,33], all of them available on the "market place" of dissipative quantum mechanics. We will apply the general formalism to the simple but non-trivial dissipative two-state or spin-boson model [22,34,35,36,37,38]. We will derive an exact formal expression for the time evolution of the average heat and heat power in the presence of an external driving field. We will show that for a constant bias the heat power can be expressed in the form of a convolution which involves the population and coherence correlation functions of the model system superimposed by polarization correlations of the heat bath accounting for heat exchange. We will present results in analytic form for weak Ohmic friction which cover both the white-noise and the coloured quantum noise regime. In Sect. 2 we will present the microscopic model and introduce the characteristic function which includes the entire statistics of the energy transfer process. Specifically, we will define the characteristic function as the trace of a generalized density matrix of the system-plus-bath entity. This will enable us perform the reduction of the generating functional to the coordinates of the system alone by a suitable generalization of the Feynman-Vernon method. The reduction will be performed for general coloured quantum noise in Sect. 3. The explicit expression for the influence functional generalizes the Feynman-Vernon functional. While the latter accounts for decoherence and relaxation of the reduced system, the generalized version captures also the quantum effects of the environment upon the time-dependent energy transfer process. To investigate potential of the method, we will consider in Sect. 4 the average heat in the dissipative two-state system and we will present an exact formal expression suitable as a source for performing dynamical simulations in the presence of external driving. In addition, we will give exact expressions in analytic form for a constant bias and weak Ohmic friction which cover the entire temperature domain ranging from the classical white noise regime down to zero temperature where coloured quantum noise prevails. We point out that these effects could be observed using a setup similar to the one proposed in Ref. [17,18], provided that the detector has sensitivity to resolve in real time the dynamics of the system at sufficiently low temperatures. Finally, in Sect 5 we will present our conclusions.

Model and general settings
We consider a quantum mechanical particle of mass M , position q and momentum p moving in a timedependent potential V (q, t). We assume that the particle is coupled to a thermal bath, modelled as a set of N independent harmonic oscillators [39,22,23,40]. The Hamiltonians of the system S and reservoir R are where x α and p α are position and momentum of the harmonic oscillator α, respectively. The system-reservoir coupling is chosen to be bi-linear in the positions variables of system and reservoir. It reads with c α the interaction strength. The second term in H I has been added to compensate renormalization of the potential V (q) caused by the bi-linear coupling term in H I [22,23]. All effects of the reservoir on the system induced by the coupling H I are captured by the spectral density of the coupling Being interested in the continuum limit N → ∞, we consider smooth spectral densities of the coupling with power-law behaviour at low frequencies and spectral cut-off, e.g., exponential cutoff, at high frequencies, Here, γ s is the coupling constant,ω a reference frequency and ω c the high-frequency cut-off. The important case of Ohmic damping is represented by s = 1 with M γ 1 ≡ η being the viscosity coefficient. The regime s > 1 (0 < s < 1) describes super-(sub-) Ohmic dissipation [41,22,34,37]. In the sequel we consistently use units in which k B = = 1.
Our goal is to determine the amount and the statistics of the exchanged energy between the system and the reservoir in a time interval t (hereafter we put the initial time t 0 = 0). This information may be obtained by performing a projective measurement of the energy of the environment both at the beginning and at the end of the time interval. Different kinds of the preparation of the initial state are feasible. In the following, we consider the density matrix of the composite system at time zero in factorized form [22], We assume that ρ R (0) is the canonical density of the unperturbed reservoir, where Z R is its partition function. Then H R commutes with ρ(0), and thus the initial measurement of the heat of the reservoir does not affect the dynamics. With this assumption, we can formulate the conditional probability P (Q, t) for the output of the two measures and, therefore, for the dissipated heat Q [19]. The central quantity of interest, which includes the entire statistics of the heat transfer process, is the characteristic function G ν (t) defined as Fourier transform of the probability distribution P (Q, t), Because of the second equality, the function G ν (t) is referred to as the moment generating function (MGF). The n-th derivative of the MGF at ν = 0 gives the nth moment, The first moment, n = 1, yields the average heat transferred to the reservoir [19], A major advantage of using the MGF is that it can be written as the trace of a generalized density operator. Thus, known techniques can be applied to solve the dynamics. Following Ref. [19], we may write where the trace is over all degrees of freedom and U (t, 0) is the unitary evolution operator of the composite system, i.e., the corresponding generator is the total Hamiltonian H(t) = H S (t) + H R + H I . In Ref. [19] the MGF has been studied within the weak-coupling Born-Markov approximation by solving a master equation associated with a generalized density matrix of the total system. In this work we will advance significantly by investigating the MGF G ν (t) from a quite different perspective. In essence, we shall derive a path integral formulation for the characteristic function G ν (t) which may serve as a basis for the study of all statistical moments of the heat transfer process, as specified in Eq. (9).
We do not specify a priori the particular form of the initial state ρ S (0). The only assumption we make for convenience is that the system-bath complex is initially in a product state and that the system-bath interaction is switched on at time t = 0 + . However, the approach envisaged here can easily be generalized to other scenarios of the initial preparation.
3. Path integral representation of the characteristic generating function G ν (t) In this section we take the route of obtaining a closed formal expression for the characteristic function G ν (t). It is advantageous to carry out the reduction of the dynamics of the system-plus-reservoir entity to the dynamics of the system alone within the functional integral description [22,23], a technique introduced by Feynman and Vernon already in 1963 [42]. The present problem, however, is substantially different because of the additional operator terms e iνHR and e −iνHR in the original expression (11). The coordinate representation of Eq. (11) reads where the N -component vector x ′ f stands for the bath coordinates (x ′ 1f , ...., x ′ N f ), and q f represents a given position of the particle. For convenience, the factor e −iνHR in Eq. (11) is included in the initial density matrix of the reservoir with the shifted canonical density operator of the reservoir ρ (ν) Using completeness of the position eigenstates we can rewrite G ν (t) as Here, is the probability amplitude for the transition from the initial state i to the final statef according to evolution with U (t, 0). The path integral representation of the propagator is [22,42,23] K Here, S S [q] is the action of the system, S R [x] the action of the reservoir, and S I [q, x] is the action associated with the interaction H I [q, x]. With the factorizing form (13) of ρ (ν) (0) one then obtains The propagating function J (ν) describes the time evolution of the RDM resulting from the generalized density operator given in Eq. (11) upon integrating out the reservoir. It reads For F (ν) [q, q ′ ] = 1, J (ν) (q f , t; q i , q ′ i , 0) describes propagation of the RDM in the absence of the systemreservoir interaction H I . This illustrates that the functional F (ν) [q, q ′ ] is the influence functional for the characteristic function G ν (t). This functional carries all effects of the system-reservoir coupling both on the dynamics of the reduced system and on the heat exchange process with the environment. The still formal expression reads Here the amplitude is the propagator of the bath in presence of the system-bath coupling H I given in Eq. (3). In the limit The expression in Eq. (18) is formally similar to the one obtained for the RDM ρ S (t) reported in [22] in the context of quantum dissipative system. However, the substantial difference is the extra factor R which entail independent integrations of the end/initial point x ′ f /i and x f /i . It is convenient to express the influence functional (20) in terms of the linear combinations Here, q 0 is a characteristic length of the system S, which is introduced to render the paths η(τ ) and ξ(τ ) dimensionless. The path η(τ ), usually referred to as quasi-classical, measures propagation of the system along the diagonal of the density matrix, whereas the path ξ(τ ) measures off-diagonal excursions, i.e., quantum fluctuations. In order to evaluate Eq. (20) we follow the standard approach [22,34] of gaussian path integral integrations. The details of the calculation are presented in Appendix A. With the resulting expressions (A.1) and (A.14) -(A.17) , the influence functional is found to read The influence actions of the first exponential factor are The local action Φ (0) loc [η, ξ] originates from the potential counterterm in the interaction (3). The nonlocal term Φ (0) cor [η, ξ] is the standard Feynman-Vernon influence action functional which governs decoherence and relaxation of the reduced system [22]. The real part of Φ The imaginary part L ′′ (τ ) is proportional to the time derivative of the friction kernel M γ(τ ) in the associated classical equation of motion (M is the inertial mass of the system), The second exponential factor in Eq. (23) controls the statistics of the heat transport process between the system and the environment, The influence action ∆Φ (ν) [η, ξ] carries in addition to the time-correlations of the aforementioned types also time-correlations between quasiclassical propagation at different time and between quantum fluctuations and subsequent quasiclassical propagation. The spectral representations of the kernels are As ν → 0, the functional e i∆Φ (ν) [η,ξ] approaches unity. With the influence functional (23) and the shortform Here we have assumed that the system is initially in a diagonal state of its density matrix. The time-dependent average heat transfer between the system and the reservoir is now obtained according to the relation (10) as where . The resulting expression for Φ (1) in whichL ′ (τ ) andL ′′ (τ ) are real and imaginary part of the time-derivative of the correlation function L(τ ) given in Eq. (25). The explicit expressions (29) and (30) for the MGF G ν (t) and the average heat transfer Q(t) with the expressions (27) and (31) for ∆Φ (ν) [η, ξ] and Φ (1) [η, ξ], respectively, are the main results of this section. They hold for general linear dissipation with any type of memory-friction. We conclude this section with a short additional consideration. Firstly, we see that the term in (31) which describes time correlations between quasiclassical paths yields a contribution to Q(t) even if the system is frozen in its initial position. The constant path q(τ ) = q ′ (τ ) = 1 2 η i q 0 yields the contribution Here p(η i ) = η i | ρ S (0) |η i is the occupation probability of the initial state η i , upon using the relation (26), we can relate the heat portion δQ i (t) to the friction kernel M γ(t), Remarkably, the heat portion (32) generally occurs when the system dynamics is slow on the scale 1/ω c . Evidently, for the thermal initial state (7) of the uncoupled reservoir, the heat contribution δQ i (t) depends on the particular initial state η i of the system. In fact, when the coupling H I [η i q 0 /2, x] is switched on at time zero, the thermal state of the uncoupled reservoir relaxes to a shifted canonical state. The initial slip (32) on the time scale 1/ω c represents the respective heat transfer. If we had calculated the functional Φ (1) [η, ξ] for a canonical initial state with shifted oscillator positions, x α −ηq 0 c α /(2m α ω 2 α ), and had subtracted the irrelevant constant polarization energy (ηq 0 ) 2 µ/8, the heat contribution would be With the choiceη = η i the initial slip term (34) would be absent. The relevant shifted thermal initial state can be arranged, e.g., by trapping the system at times t < 0 in the particular state η(t) = η i and releasing the constraint at time zero. Keeping this in mind, we disregard the initial slip contribution in the sequel.

Application to the spin-boson model: the average heat transfer
Consider a particle moving in a double well potential and in contact with its environment and assume that the thermal energy is small compared to the energy splitting in the individual wells. Then only the two lowest states become relevant to the dynamics. This scenario is conveniently described by the spinboson model. It is of great interest in different fields ranging from nuclear magnetic systems to quantum optics [34,24,25,35]. We now apply the method developed in the previous section to the driven spin-boson model [22,34,35].
Here our interest is focused on the transferred average heat and on the heat power. The Hamiltonian of the two-state system (TSS) in pseudo-spin representation is The state basis is formed by the two localized states |R and |L of the double well, which are eigenstates of σ z with eigenvalues +1 and −1, respectively. The position operator is q = q 0 σ z /2 and has eigenvalues ±q 0 /2 being associated with the minima of the double-well potential V (q). The tunneling operator σ x = |R L| − |L R| transfers the particle between the two wells with tunneling amplitude ∆, while describes an externally applied bias with a constant part ǫ 0 , and a driving term ǫ 1 (t). The MGF expression in Eq. (29) takes for the TSS in the {ξ, η}-representation the form The corresponding expression for the transfer of heat reads Here we have utilized that the local action Φ

Exact formal solution for the average heat
The double path sum in Eq. (37) can be visualized as a path sum for a single path that successively visits the four states of the RDM, i.e., the two diagonal or sojourn states labelled by η = ±1, and the two off-diagonal or blip states labelled by ξ = ±1. Starting out from a diagonal state of the RDM, the path dwells in sojourn j during time interval s j = t 2j+1 − t 2j , and in a blip during time interval τ j = t 2j − t 2j−1 , where t j is the flip time. The piecewise constant paths η(τ )-and ξ(τ ) with 2n flips are given by where t 0 = 0 and t 2n+1 = t. Here the initial state of the RDM is η i = η 0 and the final state is η f = η n . The RDM at initial time t = 0 is given by the initial populations p L = 1 2 (1 + p 0 ) and p R = 1 2 (1 − p 0 ) of the left and right well. We have −1 ≤ p 0 ≤ 1 and The path sum in the expressions (37) results in a series in the number of flip pairs, and in each term the summation of all possible arrangements in the time-ordered alternating series of sojourns and blips, Here the integration symbol D 2n {t j } is a short form of the 2n time-ordered integrations, Consider first the weight factors of the paths resulting from the internal dynamics of the TSS. The weight to switch per unit time from the diagonal state η j to the directly following off-diagonal state ξ j+1 is −i η j ξ j+1 ∆/2, and the weight to switch per unit time from the offdiagonal state ξ j+1 to the directly following diagonal state η j+1 is −i ξ j+1 η j+1 ∆/2. Thus, the aggregated weight of 2n flips depends only on the initial and final sojourn state, The weight to stay in a sojourn is unity, whereas the weight to stay in the blip ξ j during time interval ]. Accordingly, the contributions of n blips are accumulated in the bias weight factor For the piecewise constant path (39) in the influence action Φ The resulting correlation factor for 2n flips is [34,43,44] The first exponential factor of G n represents the intrablip correlations, and the second exponential factor the interblip correlations of the n blips. With the short form W j,k = W (t j − t k ) the interactions between all sojourns and all subsequent blips are in the phase factor H n . In the second form, H n is written in terms of the entire blip correlations φ n,k with sojourn k. The interactions between blip k and subsequent blip j, and the interactions between sojourn k and subsequent blip j are written as With these individual terms the average heat takes the form Here the double prime in {η j = ±1} ′′ indicates summation over the internal sojourn states η 1 , · · · , η n−1 .
For the piecewise constant paths (39) the functional (31) describes time correlations between sojourns and blips with any time-ordering. By analogy with the relation of the kernels in the actions (31) and (24), the time correlations between the {η i }-and {ξ j }-charges in the action Φ (1) [η (n) , ξ (n) ] are specified by the derivative of the function W (τ ) defined in Eq. (45). The imaginary partẆ ′′ (τ ) is connected with the damping kernel according tȯ Interestingly enough, the summation over the final sojourn states η n = ±1 in the expression (49) cancels, by reason of the extra factor η n , all correlations of blips and sojourns in Φ (1) except the correlations with the two final sojourns states. Summation of the leftover correlations yields Here U n,k (t) represents the correlations of sojourn k with sojourn n and V n,j (t) the correlations of blip j with sojourn n. Explicit dependence on the end time t is indicated for later purpose. We obtain analogous to the forms (48) the expressions We readily find where and where r = 1, 2. The time correlations carried by the kernelẆ (τ ) are in the functions The expression (53) with (54) and (55) holds for linear dissipation with arbitrary memory-friction.

The Ohmic case
Of particular importance is the case of Ohmic dissipation, which is the limit s → 1 in the spectral density of Eq. (5). In the Ohmic universality limit ω c → ∞ we have where K is a dimensionless damping strength, K = M γ 1 q 2 0 /(2π). Because the kernelẆ ′′ (τ ) is memoryless, the sojourn-sojourn correlation U n,ℓ (t) is restricted to the nearestneighbour term ℓ = n − 1 in which the intermediate sojourn has length zero, Observing that a blip of zero length does not interact with other blips, the ξ-summation in Eq. (54) cancels all sojourn-sojourn correlations except those of the first with the second. Thus we have U n,n−1 = δ n,1 U 1,0 . In addition, the odd term R (1,a) 1 (t) is zero because the single blip left over has zero length, thus readily Here δ = ∆ e −W adia /2 is a renormalized transition amplitude which takes into account the polarization cloud of the bath modes in the adiabatic limit. We have W adia = 2K ln(ω c /∆ r ) and hence δ = ∆ r ≡ ∆(∆/ω c ) K/(1−K) in the regime T ≪ ∆, and W adia = 2K ln[ω c /(2πT )] and hence δ = ∆ r (2πT /∆ r ) K in the regime T ≫ ∆. Consider next the temperature-dependent term Q(t) 2 . We see from the imaginary part of the expression (56) that the blip-sojourn correlations X j,k are restricted to the nearest-neighbour term, X j,k = πK δ j,k+1 . Thus we have φ n,k = πKξ k+1 . With this the expression (54) with (55) takes the form At this point it is expedient to turn to the average heat power P (t) , which is the time derivative of the heat, P (t) = Q (t) . Since V n,ℓ (t = t 2n ) = 0, the derivative of Q(t) 2 is restricted to the derivative of V n,ℓ (t = t 2n ) in Eq. (59). Observing thatẄ ′ (τ ) coincides with L ′ (τ ) given in Eq. (25), we readily obtain This form is the exact formal path sum expression for the average heat power of the two-state system in the Ohmic scaling limit. Apart from the last factor, the expression represents the dynamics of the RDM of the TSS under all the time-correlations carried by the intra-and inter-blip correlation factor G n . The residual factor with the time correlations of all intermediate flip times with the final time t is specific to the mean of the heat power. The expression (60) may be used to calculate the time dependence of the heat power for any time-dependent driving of the system according to Eq. (36).
To sound out validity and potential of the presented method, we now study the expression (60) for a constant bias ǫ 0 . First, we see that in the term of ∆ 2n the blip ξ ℓ , which may take any position in the sequence of n blips, plays a particular role. In the weak damping limit, we may safely disregard the interblip correlations of the preceding and the subsequent blips with the blip ℓ, Λ j,ℓ = 0 and Λ ℓ,k = 0, as well as all interblip correlations across this blip, i.e., Λ j,k = 0, where j < ℓ and k > ℓ. Then the series (60) can be expressed in terms of a convolution which involves in the time segment the population correlation function σ z (τ ) [40] and the coherence correlation function σ x (τ ) of the TSS [44,43,22]. The resulting expression is Here the components which are symmetric under bias inversion ǫ 0 → −ǫ 0 have index s, and the respective antisymmetric components have index a. The initial conditions are σ z (0) a = σ x (0) s = σ x (0) a = 0 and σ z (0) s = 1.
In the further analysis of the average heat power P (t) we can thus rely on analytic results for the nonequilibrium correlation functions σ z (t) and σ x (t) reported in the recent literature [40,41,22]. To conclude this subsection, we consider as a preparatory work for the subsequent studies the total heat transferred from the TSS to the reservoir in the weak damping regime. Writing the localized states |L and |R as linear combinations of the ground state |g and excited state |e of H S , and assuming initial preparation according to Eq. (40) and thermal occupation of the states |g and |e at time infinity, we obtain where E ini and E eq indicate the initial and equilibrium energy of the system respectively. In the following subsections, we shall apply these results to study the heat transfer for weak damping in the Markov and in the quantum noise regime.
Observing that the function L ′ (τ ) is memory-less and the initial conditions for σ z (t) and σ x (t) in Eq. (61) are σ z (0) = 1 and σ x (0) = 0, the expression (61) for the heat power reduces to the concise form Since the pair interaction W ′ (τ ) ∝ τ , the interblip correlations Λ j,k are zero. As a result, the series for the coherence correlation function σ x (t) is easily summed in Laplace space, yielding [44,22] Importantly, the residuum of the pole at the origin, which is the equilibrium value σ x eq = δ 2 /(2T ∆) of σ x (t) reached at time infinity, together with the factor πKT ∆ cancels exactly the constant term πKδ 2 /2 in Eq. (65). This secures that the heat power in fact vanishes as t → ∞, and the dynamics of P (t) is determined by the simple zeros λ 1 , λ 2 and λ 3 of the denominator in Eq. (66). The resulting expression is where "cycl." denotes addition of the terms obtained by cyclic permutation of the indices. The behaviours of λ 1 , λ 2 and λ 3 reflect the characteristics of a cubic equation with real coefficients, which depend on δ, ǫ and the scaled temperature ϑ. Physically, it is natural to write them in terms of an oscillation frequency Ω, a decoherence rate γ and a relaxation rate γ r as λ 1,2 = −γ ∓ i Ω and λ 3 = −γ r . In the temperature regime δ 2 + ǫ 2 0 ≪ T ≪ δ 2 + ǫ 2 0 /(2πK) the bath coupling H I is a perturbation. In the leading one-boson exchange approximation, the rates increase linearly with the scaled temperature ϑ, . As temperature is raised, multi-boson exchange processes become increasingly important. In the temperature range well above δ 2 + ǫ 2 0 /(2πK), the decoherence and relaxation rates show drastically different behaviours, In this so-called Kondo regime, the relaxation rate decreases with increasing temperature as T 2K−1 , whereas the decoherence rate increases linearly with T . Integration of the expression (67) yields the average heat Q(t) transferred to the reservoir until time t.
Readily the average heat transferred until time infinity is found from Eq. (67) in the Markovian regime as It is straightforward to see by use of the Vieta relations for the frequencies λ 1 , λ 2 and λ 3 , that the expression (70) holds for all T in the Markov regime. Importantly, the result (70) is just the high-temperature limit of the expression (62).

Quantum noise regime
For temperatures of the order δ 2 + ǫ 2 0 or lower, the power spectrum S(ω) is coloured because of quantum mechanical noise and the functions L ′ (τ ) and W ′ (τ ) carry memory. Nevertheless, it is possible to calculate exactly the one-boson contribution to the rates γ and γ r [41,22]. The resulting expressions are where Ω = δ 2 + ǫ 2 0 and δ = ∆ r . The correlation functions σ z (τ ) and σ x (τ ) going down into the expression of Eq. (61) take the form Next, the convolution in Eq. (61) is conveniently evaluated by using for L ′ (τ ) its spectral representation in Eq. (25) and writing the coth( ω 2T )-function in terms of the sum representation where ω n = 2πT n is a bosonic Matsubara frequency. Then the τ -integration is straightforward. The remaining ω-integral picks up (i) a contribution from the singularity at ω = Ω formed in the limit K → 0 and (ii) a contribution from the residua of the infinite sequence of poles stringed along the imaginary axis at ω = i ω n . This split-up into a system and a Matsubara part is a general feature of equilibrium correlation functions in the quantum regime. The first contribution has a part which cancels the first term in Eq. (61), and the remaining part is The contribution from the Matsubara poles is odd in the bias and is The time-dependent coefficients u 1 (t) and u 2 (t) are These functions can be written as linear combinations of hypergeometric 2 F 1 -functions which reduce to sine and cosine integral functions at T = 0. The function u 1 (t) depends logarithmically on the cutoff ω c . As the Markov regime at T well above Ω is approached, the expression (74) smoothly matches on the expression (67) with the one-boson rates (68). The heat power term P (t) 2 is a pure quantum contribution, which is negligibly small in this regime. Conversely, as the temperature is lowered, quantum effects appear in P (t) 1 through the coth-function in the amplitude and in the rate expressions (71). This function indicates emission and absorption of quanta with energy ω in thermal equilibrium. Moreover, the term P (t) 2 becomes increasingly important as T is decreased, and may even dominate the time-dependence of the heat power at low T . The combined expression can be conveniently written in terms of the amplitudes With these forms, the antisymmetric part of the total average heat power P (t = P (t) 1 + P (t) 2 is conveniently written in terms of the amplitude A rel (t) of the relaxation contribution, the amplitude A osc (t) of the oscillatory contribution, and a phase shift ϕ(t), The resulting expression for the total mean heat power is The oscillatory part of P (t) , and equivalently of the integrated power Q(t) , describes seesaw transport of heat between system and reservoir. The damped oscillatory parts are superimposed by relaxation contributions. Considering Eq. (79) and comparing it with Eq. (74), one can deduce that the degree of deviation from unity of the amplitude factors A rel (t) and A osc (t) indicate the relative strength of the pure quantum contribution P (t) 2 . Indeed, in the case A rel (t) = A osc (t) = 1 and ϕ(t) = 0 one recovers the form (74) in which the pure quantum contribution P (t) 2 = 0 is absent. The amplitudes A rel (t) and A osc (t), and the phase ϕ(t) are plotted in Figs. 1 and 2 as a functions of time for three different temperatures. The amplitude A rel (t) approaches unity in the regime t > 1/Ω for T ≪ Ω and in the regime t > 1/T for T ≫ Ω. For fixed T , the amplitude A rel (t) is monotonically decreasing as time t is lowered and eventually changes sign at time t 0 . In the temperature regime T ≫ Ω, the instant t 0 depends inversely on T as t 0 ≈ ln(2)/(2πT ) and eventually enters the core region, t 0 = O(1/ω c ), at temperatures of the order ω c . Reversely, the instant t 0 increases with decreasing temperature and reaches at T = 0 the maximal value t 0 = 1/(πΩ). The amplitude A osc (t) has a minimum at the time t 0 , at which A rel (t) is zero, whereas the phase ϕ(t) undergoes a sudden jump from −π/2 to π/2 as t passes t 0 from below. At times below t 0 , the pure quantum contribution P (t) 2 dominates the behaviour of the heat power P (t) . Both the amplitudesA osc (t) and the phase ϕ(t) approach constant values A osc,∞ and ϕ ∞ at times t ≫ 1/Ω. These values get larger as T is decreased, and at low T they are notably different from the values A osc = 1 and ϕ = 0 holding in the absence of the pure quantum contribution P (t) 2 . The asymptotic values are largest at T = 0 and are where B ∞ = arctan[2 ln(ω c /Ω) − C E ], and where C E is Euler's constant. This shows that the direct time dependence of P (t) and Q(t) at low T strongly depends on the pure quantum contribution P (t) 2 for all times. These non-Markovian effects could be measured at temperature T < Ω, provided that the detector used can resolve the signal with the sufficient accuracy for time t < 1/γ, 1/γ r . On the bottom line, the contribution of P (t) 2 to the overall heat Q ∞ is negligibly small. In the asymptotic weak-damping limit, the total heat transferred to the bath is found from the expression (79) as With the form (71) for γ r , Q ∞ coincides with the expression (62). Plots of P (t) and Q(t) are shown for different temperatures and different initial conditions in Figs. 3 and 4. The oscillations of Q(t) are clearly visible for the case T = 0.1 δ both for initial condition p 0 = 1 in panel a) and for p 0 = −1 in panel b). For p 0 = 1, the three curves stay together until time t ≈ 15 δ and then diverge to reach the respective asymptotic values indicated by the arrows on the right border. We see from the curves in Fig. 4 b), which correspond to the initial condition p 0 = −1, that the heat can be positive or negative, dependent on temperature.

Conclusions
In this paper we established the functional integral description for the time-dependent heat exchange of a quantum system coupled to a thermal reservoir. We presented the exact formal solution for the moment generating functional which carries all statistical features of the heat exchange process for general linear dissipation. We derived an exact formal expression for the transferred heat and applied this formalism to the dissipative two-state system. We showed that the difference between the dynamics of the heat transfer and the dynamics of the reduced density matrix (RDM) is an additional time-nonlocal correlation function which correlates intermediate states of the RDM with the final state.
To investigate the potential of the presented method, we calculated the dynamics of the average heat power and average heat in analytic form for weak Ohmic dissipation both in the Markovian regime relevant at high temperatures and in the non-Markovian quantum noise regime holding when temperature is of the order of the level splitting Ω or lower. In the latter regime, the heat is represented by a convolution integral which involves the population and coherence correlation functions of the dissipative TSS and the polarization correlation function of the reservoir. We find that the heat transfer receives contributions both from singularities related to the system dynamics and from Matsubara singularities resulting from the system-reservoir correlations. The latter yield significant contributions in the quantum noise regime while they are absent in the Markovian regime.
Altogether we have achieved a complete description of the dynamics of the heat transfer for weak damping ranging from the classical regime down to zero temperature

Appendix A. Gaussian integration
In this Appendix we outline the evaluation of the multiple integral expression (20). For simplicity, we denote quantities related to the reservoir oscillator α by an index α, and drop the redundant label R. Further, we use the inverse temperature β = 1/T , but return to T when we employ the expressions (A.17) in Section 3. At first, we see that the influence functional F (ν) [q, q ′ ] can be written in the product form in which the term F (ν) α [q, q ′ ] is the contribution of the bath oscillator α. We have The matrices of the density operators ρ (ν) α (0) = e −(β+i ν)Hα /Z α and e i νHα depend solely on coordinates of the bath oscillator α and may be written as The prefactor A 1α (τ ) and the Euclidean action B 1α (x, y; τ ) for an imaginary-time interval τ are given by A 1α (τ ) = m α ω α 2π sinh(ω α τ ) The phase is determined by the real-time action in the presence of the external force c α q(τ ) [42,22] φ α [q, where B 2α (x, y; t) = m α ω α 2 sin(ω α t) [(x 2 + y 2 ) cos(ω α t) − 2xy] , C α (t ′ ) = c α sin(ω α t) sin(ω α t ′ ) . (A.9) The first term in Eq. (A.8) is the phase resulting from the internal oscillator dynamics, the second and third terms depend linearly on the initial and final oscillator position and on the history of the system. Finally, ψ α [q] is a global phase which does not depend on the reservoir mode, where µ α = c 2 α /(m α ω 2 α ). Here the first term originates from the counter term in the interaction H I given in Eq. (3). We can thus write the generalized influence functional F (ν) [q, q ′ ] in the form where K α captures the integrations of the bath mode α, iα -integrals can be done by using again twice the relation (A.13). The various pre-exponential factors that accrued multiply to unity. In addition, the arising exponent is conveniently expressed in terms of the paths η(τ ) and ξ(τ ) given in Eq. (22). After all, the resulting expression for the influence functional (A.2) may be written in the form (A.14) with 2m α ω α sin(νω α /2) cosh[(β + i ν)ω α /2] sinh(βω α /2) sin(ω α τ ) .