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

In recent years enormous advancements have been accomplished in the field of quantum thermodynamics both on the fundamental and practical side [1, 2]. From a fundamental perspective, a key question is how classical thermodynamics emerges from quantum dynamics [3, 4]. From a practical point of view, there is urgent need to understand how quantum systems exchange energy and heat over time [5, 6]. This is an essential step towards building new quantum devices with potentially enormous technological impact [4, 7]. To this end, the statistics of the energy supplied to and removed from a closed quantum system by a time-dependent force has been largely investigated [1, 2, 8–15].

The case in which a driven quantum system exchanges energy with a thermal reservoir is of particular interest in view of potential applications. Here the challenge is to render an adequate description of the underlying non-unitary dynamics. Because of the similarities between exchange of energy and exchange of particles, many correspondences exist between heat transfer and electron transport. The steady-state heat current between reservoirs has been studied with a Landauer-type transport formula in a nonequilibrium spin-boson model [16, 17] and in molecular junction models [18, 19]. Generalized Landauer formulae of the Meir-Wingreen type, which take also inelastic processes into account, have been applied to determine the steady-state heat current across a molecular bridge [20] and through a nano-circuit [21], and to study the Kondo effect in heat transport via a local two-state system (TSS) [22].

Heat pumping across a two-state system has been studied in weak-coupling Redfield theory [23], in the strong-coupling noninteracting-blip approximation [24], and in the Markov approximation [25]. The heat flux in a driven resonant level model has been discussed in the adiabatic regime within a scattering approach [26]. The effects of strong driving and non-Markovian dynamics on the heat flux have been studied numerically within a stochastic Liouville–von Neumann scheme [27].

By now, the full counting statistics in steady-state electron transport is well understood on the basis of Levitov–Lesovik-type formulae [28] and master equation approaches [29–31]. Much less attention has been given to heat transport. A study of the counting statistics for heat transport through a 1D harmonic chain based on the Levitov–Lesovik formula is reported in [18]. The transient behavior of the counting statistics has been investigated within nonequilibrium Green's function techniques [32]. The counting statistics of heat transfer has also been studied with master equation approaches [33, 34] and with Lindblad theory [35].

By starting out from the standard definition of heat [1, 10, 33–36], we develop a new microscopic approach to heat exchange of a general quantum system with a thermal reservoir. Our method is based on influence functional integral techniques in which all effects related to heat exchange in the reduced system are included in a specific influence functional. For any linear environment, all conceivable quantum noise sources are captured by appropriate spectral properties of the influence kernels. As in the context of dissipative quantum dynamics [37–39], the influence functional does not depend on the particular quantum system under consideration. The formalism yields, upon combining quantum coherence with random emission and absorption of energy, the full dynamics of the heat transfer process. Importantly, one can take advantage of many numerical and non-perturbative methods [40–44], as well as different well-established analytical techniques [37, 44–50], all of them available on the 'market place' of dissipative quantum mechanics.

The general formalism will be applied to the simple but non-trivial dissipative two-state or spin-boson model (SBM) [37, 51–56]. 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 covering both the white-noise regime at high temperatures and the coloured quantum noise regime at low temperatures. The findings can potentially be tested upon using a setup similar to the one proposed in [57, 58], provided that the detector has sensitivity to resolve in real time the dynamics of the system at sufficiently low temperatures.

In section 2 we will present the microscopic model and introduce the characteristic function as the trace of a generalized density matrix of the system-plus-bath entity. In section 3 we perform the reduction of the characteristic function to the coordinates of the system alone for general coloured quantum noise and determine the respective influence functional. To investigate potential of the method, we will consider in section 4 the average heat in the dissipative TSS and we will present an exact formal expression suitable as a source for performing dynamical simulations in the presence of external driving. We will give an exact expression 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. Finally, in section 5 we will present our conclusions.

We consider aquantum mechanical particle of mass M, position q and momentum p moving in a time-dependent potential $V(q,t)$. We assume that the particle is coupled to a thermal bath, modelled as a set of N independent harmonic oscillators [37, 38, 51, 59]. We consider a single reservoir, but the extension to multiple reservoirs is straightforward. The Hamiltonians of the system S and reservoir R are

$$\begin{eqnarray}&&{{H}_{{\rm S}}}=\frac{{{p}^{2}}}{2M}+V(q,t),\qquad {\rm and}\qquad {{H}_{{\rm R}}}=\mathop{\sum }\limits_{\alpha =1}^{N}\left[ \frac{p_{\alpha }^{2}}{2{{m}_{\alpha }}}+\frac{1}{2}{{m}_{\alpha }}\omega _{\alpha }^{2}x_{\alpha }^{2} \right].\end{eqnarray} \tag{ 1 }$$

Here ${{x}_{\alpha }}$ and ${{p}_{\alpha }}$ 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

$$\begin{eqnarray}&&{{H}_{{\rm I}}}=-q\mathop{\sum }\limits_{\alpha =1}^{N}{{c}_{\alpha }}{{x}_{\alpha }}+\frac{1}{2}\mathop{\sum }\limits_{\alpha =1}^{N}\frac{c_{\alpha }^{2}}{{{m}_{\alpha }}\omega _{\alpha }^{2}}{{q}^{2}}\;,\end{eqnarray} \tag{ 2 }$$

with ${{c}_{\alpha }}$ the interaction strength. The second term in ${{H}_{{\rm I}}}$ has been added to compensate renormalization of the potential $V(q,t)$ caused by the bi-linear coupling term in ${{H}_{{\rm I}}}$ [37, 38]. All effects of the reservoir on the system induced by the coupling ${{H}_{{\rm I}}}$ are captured by the spectral density of the coupling

$$\begin{eqnarray}&&J(\omega )=\frac{\pi }{2}\mathop{\sum }\limits_{\alpha =1}^{N}\frac{c_{\alpha }^{2}}{{{m}_{\alpha }}{{\omega }_{\alpha }}}\delta (\omega -{{\omega }_{\alpha }}).\end{eqnarray} \tag{ 3 }$$

Being interested in the continuum limit $N\to \infty$, 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

$$\begin{eqnarray}&&J(\omega )=M{{\gamma }_{s}}\omega {{(\omega /\tilde{\omega })}^{s-1}}{{{\rm e}}^{-\omega /{{\omega }_{{\rm c}}}}}.\end{eqnarray} \tag{ 4 }$$

Here, ${{\gamma }_{s}}$ is the coupling constant, $\tilde{\omega }$ a reference frequency and ${{\omega }_{{\rm c}}}$ the high-frequency cut-off. The case of Ohmic damping is represented by s = 1. The regime $s\gt 1$ ($0\lt s\lt 1$) describes super- (sub-) Ohmic dissipation [37, 52, 55, 60]. In the sequel we consistently use units in which $k_{{\rm B}}^{{}}=\hbar =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. 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. Here we suppose that at initial time t = 0 the system and the bath are uncoupled, and the bath is in thermal equilibrium

$$\begin{eqnarray}&&\rho (0)={{\rho }_{{\rm S}}}(0)\;\otimes \;{{\rho }_{{\rm R}}}(0)={{\rho }_{{\rm S}}}(0)\;\otimes \;Z_{{\rm R}}^{-1}{{{\rm e}}^{-{{H}_{{\rm R}}}/T}},\end{eqnarray} \tag{ 5 }$$

where ${{Z}_{{\rm R}}}$ is the partition function of the reservoir. We do not specify a priori the particular form of the initial state ${{\rho }_{{\rm S}}}(0)$. Importantly, different kinds of the initial preparation are feasible [37].

With the initial state (5) the Hamiltonian ${{H}_{{\rm R}}}$ commutes with $\rho (0)$. Hence the measurement of the heat of the reservoir at t = 0 does not affect the dynamics. As a result, we can formulate the conditional probability $P(Q,t)$ for the output of the two measures and, therefore, for the dissipated heat Q [33].

The central quantity of interest, which includes the entire statistics of the heat transfer process, is the characteristic function ${{G}_{\nu }}(t)$ defined as Fourier transform of the probability distribution $P(Q,t)$,

$$\begin{eqnarray}&&{{G}_{\nu }}(t)=\int _{-\infty }^{+\infty }{\rm d}Q\;P(Q,t)\;{{{\rm e}}^{{\rm i}\nu Q}}=\mathop{\sum }\limits_{n=0}^{\infty }\frac{{{({\rm i}\nu )}^{n}}}{n!}\langle {{Q}^{n}}(t)\rangle .\end{eqnarray} \tag{ 6 }$$

Because of the second equality, the function ${{G}_{\nu }}(t)$ is referred to as the moment generating function (MGF).

The nth derivative of the MGF at $\nu =0$ gives the nth moment

$$\begin{eqnarray}&&\langle {{Q}^{n}}(t)\rangle ={{\left. {{(-{\rm i})}^{n}}\frac{{{{\rm d}}^{n}}}{{\rm d}{{\nu }^{n}}}{{G}_{\nu }}(t) \right|}_{\nu =0}}.\end{eqnarray} \tag{ 7 }$$

The first moment, n = 1, yields the average heat $\langle Q(t)\rangle$ transferred to the reservoir [33]. In addition, the nth cumulant of $P(Q,t)$ is obtained by taking the nth derivative of ${\rm ln} {{G}_{\nu }}(t)$ with respect to ν at $\nu =0$.

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 [33], we may write

$$\begin{eqnarray}&&{{G}_{\nu }}(t)={\rm Tr}\;[{{{\rm e}}^{{\rm i}\nu {{H}_{{\rm R}}}}}U(t,0){{{\rm e}}^{-{\rm i}\nu {{H}_{{\rm R}}}}}\rho (0){{U}^{\dagger }}(t,0)],\end{eqnarray} \tag{ 8 }$$

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}_{{\rm S}}}(t)+{{H}_{{\rm R}}}+{{H}_{{\rm I}}}$.

In [33] 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}_{\nu }}(t)$ from a quite different perspective. In essence, we shall derive the path integral representation of the characteristic function ${{G}_{\nu }}(t)$. This may serve as a firm basis for the study of all statistical moments of the heat transfer process, as specified in equation (7).

In this section we take the route of obtaining a closed formal expression for the characteristic function ${{G}_{\nu }}(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 [37, 38], a technique introduced by Feynman and Vernon already in 1963 [61]. The present problem, however, is rather intricate because of the additional operator terms ${{{\rm e}}^{{\rm i}\nu {{H}_{{\rm R}}}}}$ and ${{{\rm e}}^{-{\rm i}\nu {{H}_{{\rm R}}}}}$ in the original expression (8).

Using completeness of the position eigenstates the coordinate representation of equation (8) may be written as

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \qquad {{G}_{\nu }}(t) & = & \int {\rm d}{{{\boldsymbol{x}} }_{f}}{\rm d}{\boldsymbol{x}} _{f}^{\prime }{\rm d}{{{\boldsymbol{x}} }_{i}}{\rm d}{\boldsymbol{x}} _{i}^{\prime }\int {\rm d}{{q}_{f}}{\rm d}{{q}_{i}}{\rm d}q_{i}^{\prime } \\ {} & {} & \times \ \langle {\boldsymbol{x}} _{f}^{\prime }|{{{\rm e}}^{{\rm i}\nu {{H}_{{\rm R}}}}}|{{{\boldsymbol{x}} }_{f}}\rangle \langle {{{\boldsymbol{x}} }_{f}},{{q}_{f}}|U(t,0)|{{{\boldsymbol{x}} }_{i}},{{q}_{i}}\rangle \langle {{{\boldsymbol{x}} }_{i}},{{q}_{i}}|{{\rho }^{(\nu )}}(0)|{\boldsymbol{x}} _{i}^{\prime },q_{i}^{\prime }\rangle \\ {} & {} & \times \langle {\boldsymbol{x}} _{i}^{\prime },q_{i}^{\prime }|{{U}^{\dagger }}(t,0)|{\boldsymbol{x}} _{f}^{\prime },{{q}_{f}}\rangle . \\ \end{array}\end{eqnarray} \tag{ 9 }$$

The N-component vector ${\boldsymbol{x}} _{f}^{\prime }$ represents the bath coordinates $(x_{1f}^{\prime },....,x_{Nf}^{\prime })$, and q_f is the position of the particle. For convenience, the factor ${{{\rm e}}^{-{\rm i}\nu {{H}_{{\rm R}}}}}$ in equation (8) is included in the initial density matrix of the bath

$$\begin{eqnarray}&&{{\rho }^{(\nu )}}(0)={{\rho }_{{\rm S}}}(0)\;\otimes \;Z_{{\rm R}}^{-1}{{{\rm e}}^{-{{H}_{{\rm R}}}/T-{\rm i}\,\nu {{H}_{{\rm R}}}}}.\end{eqnarray} \tag{ 10 }$$

The expression (9) is formally similar to the one obtained for the reduced density matrix (RDM) ${{\rho }_{{\rm S}}}(t)$ reported in [37] in the context of dissipative quantum dynamics. However, the substantial difference is the shifted canonical density operator $\rho _{{\rm R}}^{(\nu )}$ and the extra factor $\langle {\boldsymbol{x}} _{f}^{\prime }|{{{\rm e}}^{{\rm i}\nu {{H}_{{\rm R}}}}}|{{{\boldsymbol{x}} }_{f}}\rangle$. The latter entails independent integrations of the end/initial point ${\boldsymbol{x}} _{f/i}^{\prime }$ and ${{{\boldsymbol{x}} }_{f/i}}$. The functional integral representation of the expression (9) reads

$$\begin{eqnarray}&&{{G}_{\nu }}(t)=\int {\rm d}{{q}_{f}}{\rm d}{{q}_{i}}{\rm d}q_{i}^{\prime }\;\langle {{q}_{i}}|{{\rho }_{{\rm S}}}(0)|q_{i}^{\prime }\rangle \;\int _{q(0)={{q}_{i}}}^{q(t)={{q}_{f}}}\mathcal{D}q\;{{{\rm e}}^{{\rm i}{{S}_{{\rm S}}}[q]}}\int _{q^{\prime} (0)=q_{i}^{\prime }}^{q^{\prime} (t)={{q}_{f}}}\mathcal{D}q^{\prime} \;{{{\rm e}}^{-{\rm i}{{S}_{{\rm S}}}[q^{\prime} ]}}{{\mathcal{F}}^{(\nu )}}[q,q^{\prime} ],\end{eqnarray} \tag{ 11 }$$

where ${{S}_{{\rm S}}}[q]$ is the action of the system. The double path sum represents the propagation of the RDM, while the influence functional ${{\mathcal{F}}^{(\nu )}}[q,q^{\prime} ]$ accounts for all effects of the environmental coupling relevant to the generating function ${{G}_{\nu }}(t)$. The influence functional ${{\mathcal{F}}^{(\nu )}}[q,q^{\prime} ]$ is found upon integrating out the bath variables in the expression (9). Details of the calculation are presented in appendix A.

It is convenient to rewrite the expression (11) in terms of the linear combinations

$$\begin{eqnarray}&&\eta (\tau )=[q(\tau )+q^{\prime} (\tau )]/{{q}_{0}},\qquad {\rm and}\qquad \xi (\tau )=[q(\tau )-q^{\prime} (\tau )]/{{q}_{0}}.\end{eqnarray} \tag{ 12 }$$

Here, q₀ is a characteristic length of the system S, introduced to render the paths $\eta (\tau )$ and $\xi (\tau )$ dimensionless. The path $\eta (\tau )$, usually referred to as quasi-classical, measures propagation of the system along the diagonal of the density matrix, whereas the path $\xi (\tau )$ measures off-diagonal excursions, i.e., quantum fluctuations. Assuming that the system is initially in a diagonal state of its density matrix, and employing the shortform ${{S}_{{\rm S}}}[\eta ,\xi ]={{S}_{{\rm S}}}[q]-{{S}_{{\rm S}}}[q^{\prime} ]$ the MGF then takes the form

$$\begin{eqnarray}&&\qquad {{G}_{\nu }}(t)=\int {\rm d}{{\eta }_{i}}\;\langle {{\eta }_{i}}|{{\rho }_{{\rm S}}}(0)|{{\eta }_{i}}\rangle \int {\rm d}{{\eta }_{f}}\int _{\eta (0)={{\eta }_{i}}}^{\eta (t)={{\eta }_{f}}}\mathcal{D}\eta \int _{\xi (0)=0}^{\xi (t)=0}\mathcal{D}\xi \ {{{\rm e}}^{{\rm i}\,{{S}_{{\rm S}}}[\eta ,\xi ]}}\ {{\mathcal{F}}^{(\nu )}}[\eta ,\xi ].\end{eqnarray} \tag{ 13 }$$

With the results of appendix A, equations (A.3) and (A.14)–(A.17) , the influence functional is found to read

$$\begin{eqnarray}&&\qquad {{\mathcal{F}}^{(\nu )}}[\eta ,\xi ]={{\mathcal{F}}_{{\rm FV}}}[\eta ,\xi ]\;{{{\rm e}}^{{\rm i}\Delta {{\Phi }^{(\nu )}}[\eta ,\xi ]}}\qquad {\rm with}\qquad {{\mathcal{F}}_{{\rm FV}}}[\eta ,\xi ]={{{\rm e}}^{-(\Phi _{{\rm loc}}^{(0)}[\eta ,\xi ]+\Phi _{{\rm cor}}^{(0)}[\eta ,\xi ])}}.\end{eqnarray} \tag{ 14 }$$

The first factor is the familiar Feynman–Vernon influence functional represented by the actions

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \Phi _{{\rm loc}}^{(0)}[\eta ,\xi ] & = & \ {\rm i}\;\frac{\mu q_{0}^{2}}{2}\int _{0}^{t}{\rm d}t^{\prime} \;\xi (t^{\prime} )\eta (t^{\prime} )\quad \ \ {\rm with}\quad \ \ \mu =\frac{2}{\pi }\int _{0}^{\infty }{\rm d}\omega \frac{J(\omega )}{\omega }, \\ \Phi _{{\rm cor}}^{(0)}[\eta ,\xi ] & = & \int _{0}^{t}{\rm d}t^{\prime} \int _{0}^{{{t}^{\prime }}}{\rm d}t^{\prime\prime} \;\xi (t^{\prime} )[L^{\prime} (t^{\prime} -t^{\prime\prime} )\xi (t^{\prime\prime} )+{\rm i}\;L^{\prime\prime} (t^{\prime} -t^{\prime\prime} )\eta (t^{\prime\prime} )]. \\ \end{array}\end{eqnarray} \tag{ 15 }$$

The local action $\Phi _{{\rm loc}}^{(0)}[\eta ,\xi ]$ originates from the potential counterterm in the interaction (2). The time-nonlocal action functional $\Phi _{{\rm cor}}^{(0)}[\eta ,\xi ]$ governs decoherence and relaxation of the reduced system [37]. The real part of $\Phi _{{\rm cor}}^{(0)}[\eta ,\xi ]$ is the noise action, and the imaginary part is the friction action. The former describes time-correlations of the quantum fluctuations, whereas the latter describes time-correlations between quasiclassical propagation and subsequent off-diagonal excursions. The kernel $L(\tau )$ is the equilibrium time-correlation function of the collective bath mode $\mathcal{E}(t)={{q}_{0}}{{\sum }_{\alpha }}{{c}_{\alpha }}{{x}_{\alpha }}(t)$,

$$\begin{eqnarray}&&L(\tau )=L^{\prime} (\tau )+{\rm i}L^{\prime\prime} (\tau )=\frac{q_{0}^{2}}{\pi }\int _{0}^{\infty }{\rm d}\omega \;J(\omega )\left[ \;{\rm coth} \left( \frac{\omega }{2T} \right){\rm cos} \omega \tau -{\rm i}{\rm sin} \omega \tau \; \right].\end{eqnarray} \tag{ 16 }$$

The second influence functional factor ${{{\rm e}}^{{\rm i}\Delta {{\Phi }^{(\nu )}}[\eta ,\xi ]}}$ in equation (14) controls the entire statistics of the heat transport process between the system and the environment

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \Delta {{\Phi }^{(\nu )}}[\eta ,\xi ] & = & \int _{0}^{t}{\rm d}t^{\prime} \int _{0}^{{{t}^{\prime }}}{\rm d}t^{\prime\prime} \left\{ [\eta (t^{\prime} )\eta (t^{\prime\prime} )-\xi (t^{\prime} )\xi (t^{\prime\prime} )]L_{1}^{(\nu )}(t^{\prime} -t^{\prime\prime} ) \right. \\ {} & {} & \left. +[\eta (t^{\prime} )\xi (t^{\prime\prime} )-\xi (t^{\prime} )\eta (t^{\prime\prime} )]L_{2}^{(\nu )}(t^{\prime} -t^{\prime\prime} ) \right\}. \\ \end{array}\end{eqnarray} \tag{ 17 }$$

The influence action $\Delta {{\Phi }^{(\nu )}}[\eta ,\xi ]$ 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

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} L_{1}^{(\nu )}(\tau ) & = & \ \frac{q_{0}^{2}}{\pi }\int _{0}^{\infty }{\rm d}\omega \;\frac{{\rm sin} (\nu \;\omega /2){\rm sinh} [(\omega /T+{\rm i}\;\nu \;\omega )/2]}{{\rm sinh} [\omega /(2T)]}J(\omega ){\rm cos} (\omega \tau ), \\ L_{2}^{(\nu )}(\tau ) & = & -{\rm i}\frac{q_{0}^{2}}{\pi }\int _{0}^{\infty }{\rm d}\omega \;\frac{{\rm sin} (\nu \;\omega /2){\rm cosh} [(\omega /T+{\rm i}\;\nu \;\omega )/2]}{{\rm sinh} [\omega /(2T)]}J(\omega ){\rm sin} (\omega \tau ). \\ \end{array}\end{eqnarray} \tag{ 18 }$$

Evidently, ${{\mathcal{F}}^{(\nu )}}[q,q^{\prime} ]$ coincides with the Feynman–Vernon functional at $\nu =0$, ${{\mathcal{F}}^{(0)}}[q,q^{\prime} ]={{\mathcal{F}}_{{\rm FV}}}[q,q^{\prime} ]$.

The nth moment of the probability distribution $P(Q,t)$ is readily found from the relation (7) as

$$\begin{eqnarray}&&\langle {{Q}^{n}}(t)\rangle =\int {\rm d}{{\eta }_{i}}\;\langle {{\eta }_{i}}|{{\rho }_{{\rm S}}}(0)|{{\eta }_{i}}\rangle \int {\rm d}{{\eta }_{f}}\int _{\eta (0)={{\eta }_{i}}}^{\eta (t)={{\eta }_{f}}}\mathcal{D}\eta \int _{\xi (0)=0}^{\xi (t)=0}\mathcal{D}\xi \ {{{\rm e}}^{{\rm i}\,{{S}_{{\rm S}}}[\eta ,\xi ]}}\ {{\mathcal{F}}_{{\rm FV}}}[\eta ,\xi ]\;{{\Phi }^{(n)}}[\eta ,\xi ]\end{eqnarray} \tag{ 19 }$$

with

$$\begin{eqnarray}&&{{\Phi }^{(n)}}[\eta ,\xi ]={{(-{\rm i})}^{n}}{{\left. \frac{{{{\rm d}}^{n}}}{{\rm d}{{\nu }^{n}}}{{{\rm e}}^{{\rm i}\Delta {{\phi }^{(\nu )}}[\eta ,\xi ]}} \right|}_{\nu =0}}.\end{eqnarray} \tag{ 20 }$$

The energy functional ${{\Phi }^{(1)}}[\eta ,\xi ]$ in the average heat transfer $\langle Q(t)\rangle$ is expressed in terms of the time-derivative $\dot{L}(\tau )$ of the equilibrium correlation function of the collective bath mode (16) as

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\Phi }^{(1)}}[\eta ,\xi ] & = & \frac{1}{2}\int _{0}^{t}{\rm d}t^{\prime} \int _{0}^{{{t}^{\prime }}}{\rm d}t^{\prime\prime} \left\{ [\xi (t^{\prime} )\xi (t^{\prime\prime} )-\eta (t^{\prime} )\eta (t^{\prime\prime} )]{{{\dot{L}}}^{\prime\prime }}(t^{\prime} -t^{\prime\prime} ) \right. \\ {} & {} & \qquad \qquad \quad \left. +\;{\rm i}[\eta (t^{\prime} )\xi (t^{\prime\prime} )-\xi (t^{\prime} )\eta (t^{\prime\prime} )]{{{\dot{L}}}^{\prime }}(t^{\prime} -t^{\prime\prime} ) \right\}. \\ \end{array}\end{eqnarray} \tag{ 21 }$$

The explicit expressions (13) and (19) for the MGF ${{G}_{\nu }}(t)$ and the nth moment $\langle {{Q}^{n}}(t)\rangle$ with the expressions (17) and (20) for $\Delta {{\Phi }^{(\nu )}}[\eta ,\xi ]$ and ${{\Phi }^{(n)}}[\eta ,\xi ]$, respectively, are the main results of this section. They hold for general linear dissipation with any type of memory-friction. Evidently, the energy functional (21) carries correlations of the reduced system between different times. Increasingly more convoluted correlations appear, if we proceed to higher moments of the heat transfer process. Considering the already respectable complexity inherent in the first moment, we focus in the sequel on the discussion of the average heat and postpone discussion of statistical properties with the method proposed here to subsequent work.

We conclude this section with a short additional consideration. Firstly, we see that the term in (21) which describes time correlations between quasiclassical paths yields a contribution to $\langle Q(t)\rangle$ even if the system is frozen in its initial position. The constant path $q(\tau )=q^{\prime} (\tau )=\frac{1}{2}{{\eta }_{i}}{{q}_{0}}$ yields the contribution

$$\begin{eqnarray}&&\langle \delta Q(t)\rangle =\int {\rm d}{{\eta }_{i}}\;p({{\eta }_{i}})\langle \delta {{Q}_{i}}(t)\rangle \qquad {\rm with}\qquad \langle \delta {{Q}_{i}}(t)\rangle ={{\Phi }^{(1)}}[{{\eta }_{i}},0].\end{eqnarray} \tag{ 22 }$$

Here $p({{\eta }_{i}})=\langle {{\eta }_{i}}|{{\rho }_{{\rm S}}}(0)|{{\eta }_{i}}\rangle$ is the occupation probability of the initial state ${{\eta }_{i}}$. Upon using the relation $L^{\prime\prime} (\tau )=Mq_{0}^{2}\;\dot{\gamma }(\tau )/2$, where M is the inertial mass of the system, the heat portion $\langle \delta {{Q}_{i}}(t)\rangle$ can be directly related to the friction kernel $M\gamma (t)$,

$$\begin{eqnarray}&&\langle \delta {{Q}_{i}}(t)\rangle ={{({{\eta }_{i}}{{q}_{0}}/2)}^{2}}M\left[ \;\gamma ({{0}^{+}})-\gamma (t) \right].\end{eqnarray} \tag{ 23 }$$

Remarkably, the heat portion (22) generally occurs when the system dynamics is slow on the scale $1/{{\omega }_{{\rm c}}}$.

Apparently, for the thermal initial state (5) of the uncoupled reservoir, the heat contribution $\langle \delta {{Q}_{i}}(t)\rangle$ depends on the particular initial state ${{\eta }_{i}}$ of the system. In fact, when the coupling ${{H}_{{\rm I}}}[{{\eta }_{i}}{{q}_{0}}/2,{\boldsymbol{x}} ]$ is switched on at time zero, the thermal state of the uncoupled reservoir relaxes to a shifted canonical state. The initial slip (22) on the time scale $1/{{\omega }_{{\rm c}}}$ represents the respective heat transfer. If we had calculated the functional ${{\Phi }^{(1)}}[\eta ,\xi ]$ for a canonical initial state with shifted oscillator positions, ${{x}_{\alpha }}-\bar{\eta }{{q}_{0}}{{c}_{\alpha }}/(2{{m}_{\alpha }}\omega _{\alpha }^{2})$, and had subtracted the irrelevant constant polarization energy ${{(\bar{\eta }{{q}_{0}})}^{2}}\;\mu /8$, the heat contribution would be

$$\begin{eqnarray}&&\langle \delta {{Q}_{i}}(t)\rangle ={{({{q}_{0}}/2)}^{2}}\;{{\eta }_{i}}({{\eta }_{i}}-\bar{\eta })M\left[ \;\gamma ({{0}^{+}})-\gamma (t)\; \right].\end{eqnarray} \tag{ 24 }$$

With the choice $\bar{\eta }={{\eta }_{i}}$ the initial slip term (24) would be absent. The relevant shifted thermal initial state can be arranged, e.g., by trapping the system at times $t\lt 0$ in the particular state $\eta (t)={{\eta }_{i}}$ and releasing the constraint at time zero. Keeping this in mind, we disregard the initial slip contribution in the sequel.

Consider a particle moving in a double well potential and being 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 dissipative two-state or SBM. The SBM is of great interest in different fields ranging from nuclear magnetic systems to quantum optics [37, 40, 41, 51–53].

In this section we apply the method developed in the previous section to the driven SBM [37, 52, 53]. Here our interest is focused on the transferred average heat and on the heat power. The Hamiltonian of the TSS in pseudo-spin representation is

$$\begin{eqnarray}&&{{H}_{{\rm S}}}=-\frac{\Delta }{2}{{\sigma }_{x}}-\frac{\epsilon (t)}{2}{{\sigma }_{z}}.\end{eqnarray} \tag{ 25 }$$

The state basis is formed by the two localized states $|R\rangle$ and $|L\rangle$ of the double well, which are eigenstates of ${{\sigma }_{z}}$ with eigenvalues +1 and −1, respectively. The position operator is $q={{q}_{0}}{{\sigma }_{z}}/2$ and has eigenvalues $\pm {{q}_{0}}/2$ being associated with the minima of the double-well potential V(q). The tunneling operator ${{\sigma }_{x}}=|R\rangle \langle L|-|L\rangle \langle R|$ transfers the particle between the two wells with tunneling amplitude Δ, while

$$\begin{eqnarray}&&\epsilon (t)={{\epsilon }_{0}}+{{\epsilon }_{1}}(t)\end{eqnarray} \tag{ 26 }$$

describes an externally applied bias with a constant part ${{\epsilon }_{0}}$, and a driving term ${{\epsilon }_{1}}(t)$.

4.1. Exact formal solution for the average heat

For the TSS the double path sum (13) 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 $\eta =\pm 1$, and the two off-diagonal or blip states labelled by $\xi =\pm 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 ${{\tau }_{j}}={{t}_{2j}}-{{t}_{2j-1}}$, where t_j is the flip time. For a path with $2m$ flips the initial state of the RDM is ${{\eta }_{i}}={{\eta }_{0}}$ and the final state is ${{\eta }_{f}}={{\eta }_{m}}$. We choose that the RDM at initial time t = 0 is given by the populations $p_{L}^{{}}=\frac{1}{2}(1+{{p}_{0}})$ and $p_{{\rm R}}^{{}}=\frac{1}{2}(1-{{p}_{0}})$ of the left and right well. We have $-1\leqslant {{p}_{0}}\leqslant 1$ and

$$\begin{eqnarray}&&\langle {{\eta }_{0}}|{{\rho }_{{\rm S}}}(0)|{{\eta }_{0}}\rangle =\frac{1}{2}(1-{{p}_{0}}){{\delta }_{{{\eta }_{0}},1}}+\frac{1}{2}(1+{{p}_{0}}){{\delta }_{{{\eta }_{0}},-1}}.\end{eqnarray} \tag{ 27 }$$

The analysis of the path sum (19) in the SBM is reported in appendix B. The resulting exact formal series expression for the nth moment $\langle {{Q}^{n}}(t)\rangle$ is found to read

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \langle {{Q}^{n}}(t)\rangle & = & \mathop{\sum }\limits_{{{\eta }_{0}}=\pm 1}\langle {{\eta }_{0}}|{{\rho }_{{\rm S}}}(0)|{{\eta }_{0}}\rangle \;{{\eta }_{0}}\;\mathop{\sum }\limits_{m=1}^{\infty }{{\left( \frac{-{{\Delta }^{2}}}{4} \right)}^{m}}\int _{0}^{t}{{\mathcal{D}}_{2m}}\{{{t}_{j}}\} \\ {} & {} & \times \mathop{\sum }\limits_{\{{{\xi }_{j}}=\pm 1\}}{{B}_{m}}\;{{G}_{m}}\mathop{\sum }\limits_{{{\{{{\eta }_{j}}=\pm 1\}}^{\prime\prime }}}{{H}_{m}}\mathop{\sum }\limits_{{{\eta }_{m}}=\pm 1}{{\eta }_{m}}\;{{\Phi }^{(n)}}[{{\eta }^{(m)}},{{\xi }^{(m)}}]. \\ \end{array}\end{eqnarray} \tag{ 28 }$$

The integration symbol $\int _{0}^{t}{{\mathcal{D}}_{2m}}\{{{t}_{j}}\}$ is a short form for the $2m$ time-ordered integrations, equation (B.3). The bias weight of the m blips is in the factor B_m defined in equation (B.4). The factors G_m and H_m result from the noise action and from the friction action in equation (15), respectively, and they are defined in equation (B.6) with (B.7). The weight factor G_m carries the intrablip and interblip correlations of the m blips, while the weight factor H_m includes the correlations of the m sojourns with the respective subsequent blips. Finally, the double prime in $\{{{\eta }_{j}}=\pm 1\}^{\prime\prime}$ indicates summation over the internal sojourn states ${{\eta }_{1}},\cdots ,{{\eta }_{m-1}}$.

From now on we restrict ourselves to the discussion of the average heat $\langle Q(t)\rangle$. For the piecewise constant paths (B.1) the functional (21) describes time correlations between sojourns and blips with any time-ordering. By analogy with the relation of the kernels in the actions (21) and (15), the time correlations between the $\{{{\eta }_{i}}\}$- and $\{{{\xi }_{j}}\}$-charges in the action ${{\Phi }^{(1)}}[{{\eta }^{(m)}},{{\xi }^{(m)}}]$ are specified by the derivative of the function $W(\tau )$ defined in equation (B.5). The imaginary part ${{\dot{W}}^{\prime\prime }}(\tau )$ is connected with the damping kernel by ${{\dot{W}}^{\prime\prime }}(\tau )=(Mq_{0}^{2}/2)\gamma (\tau )$.

Interestingly enough, the summation over the final sojourn states ${{\eta }_{m}}=\pm 1$ in the expression (28) for the first moment n = 1, by reason of the extra factor ${{\eta }_{m}}$, cancels all correlations of blips and sojourns in ${{\Phi }^{(1)}}$ except the correlations with the two final sojourns states. Summation of the leftover correlations yields

$$\begin{eqnarray}&&\mathop{\sum }\limits_{{{\eta }_{m}}=\pm 1}{{\eta }_{m}}{{\Phi }^{(1)}}[{{\eta }^{(m)}},{{\xi }^{(m)}}]=\mathop{\sum }\limits_{k=0}^{m-1}{{\mathcal{U}}_{m,k}}(t)\;{{\eta }_{k}}+{\rm i}\mathop{\sum }\limits_{j=1}^{m}{{\mathcal{V}}_{m,j}}(t){{\xi }_{j}}.\end{eqnarray} \tag{ 29 }$$

Here ${{\mathcal{U}}_{m,k}}(t)$ represents the correlations of sojourn k with the final sojourn m and ${{\mathcal{V}}_{m,j}}(t)$ the correlations of blip j with the final sojourn m. Explicit dependence on the end time t is indicated for later purpose. We obtain analogous to the forms (B.7) the expressions

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\mathcal{U}}_{m,k}}(t) & = & {{{\dot{W}}}^{\prime\prime }}(t-{{t}_{2k+1}})+{{{\dot{W}}}^{\prime\prime }}({{t}_{2m}}-{{t}_{2k}})-{{{\dot{W}}}^{\prime\prime }}(t-{{t}_{2k}})-{{{\dot{W}}}^{\prime\prime }}({{t}_{2m}}-{{t}_{2k+1}}), \\ {{\mathcal{V}}_{m,k}}(t) & = & {{{\dot{W}}}^{\prime }}(t-{{t}_{2k-1}})+{{{\dot{W}}}^{\prime }}({{t}_{2m}}-{{t}_{2k}})-{{{\dot{W}}}^{\prime }}(t-{{t}_{2k}})-{{{\dot{W}}}^{\prime }}({{t}_{2m}}-{{t}_{2k-1}}). \\ \end{array}\end{eqnarray} \tag{ 30 }$$

Next, we add up the contributions of the $2m-1$ intermediate sojourns. Further, we utilize that in equation (28) the contributions of the blip states which are odd under the substitutions $\{{{\xi }_{j}}\}\to \{-{{\xi }_{j}}\}$ cancel each other out. The resulting expression for $\langle Q(t)\rangle$ is conveniently divided into the $\mathcal{U}$- and $\mathcal{V}$-contributions, ${{\langle Q(t)\rangle }_{1}}$ and ${{\langle Q(t)\rangle }_{2}}$, and these into those which are symmetric (s) and antisymmetric (a) under bias inversion [$\epsilon (t)\to -\epsilon (t)$]. We readily find

$$\begin{eqnarray}&&\langle Q(t)\rangle ={{\langle Q(t)\rangle }_{1}}+{{\langle Q(t)\rangle }_{2}},\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray}&&\qquad {{\langle Q(t)\rangle }_{{\rm R}}}=\frac{1}{2}\;\mathop{\sum }\limits_{m=1}^{\infty }{{\left( \frac{-{{\Delta }^{2}}}{2} \right)}^{m}}\int _{0}^{t}{{\mathcal{D}}_{2m}}\{{{t}_{j}}\}\mathop{\sum }\limits_{\{{{\xi }_{j}}=\pm 1\}}\left( {{\mathop{\prod }\limits_{k=0}^{m-1}}^{{}}}{\rm cos} ({{\phi }_{m,k}}) \right)\;{{G}_{m}}\left[ \mathcal{R}_{m}^{(r,{\rm s})}(t)+R_{m}^{(r,{\rm a})}(t)\; \right],\end{eqnarray} \tag{ 32 }$$

and where $r=1,2$. The time correlations carried by the kernel $\dot{W}(\tau )$ are in the functions

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \mathcal{R}_{m}^{(1,{\rm s})}(t) & = & {\rm cos} ({{\varphi }_{m}})\left[ \;{{\mathcal{U}}_{m,0}}(t)-{\rm tan} ({{\phi }_{m,0}})\mathop{\sum }\limits_{\ell =1}^{m-1}{\rm tan} ({{\phi }_{m,\ell }}){{\mathcal{U}}_{m,\ell }}(t)\; \right], \\ \mathcal{R}_{m}^{(1,{\rm a})}(t) & = & {{p}_{0}}{\rm sin} ({{\varphi }_{m}})\mathop{\sum }\limits_{\ell =0}^{m-1}{\rm tan} ({{\phi }_{m,\ell }})\;{{\mathcal{U}}_{m,\ell }}(t), \\ \mathcal{R}_{m}^{(2,{\rm s})}(t) & = & -{\rm cos} ({{\varphi }_{m}}){\rm tan} ({{\phi }_{m,0}})\mathop{\sum }\limits_{\ell =1}^{m}{{\mathcal{V}}_{m,\ell }}(t){{\xi }_{\ell }}, \\ \mathcal{R}_{m}^{(2,{\rm a})}(t) & = & {{p}_{0}}{\rm sin} ({{\varphi }_{m}})\mathop{\sum }\limits_{\ell =1}^{m}{{\mathcal{V}}_{m,\ell }}(t)\;{{\xi }_{\ell }}. \\ \end{array}\end{eqnarray} \tag{ 33 }$$

The expression (31) with (32) and (33) holds for linear dissipation with arbitrary memory-friction.

4.2. The Ohmic case

Of particular importance is the case of Ohmic dissipation, which is the limit $s\to 1$ in the spectral density of equation (4). In the Ohmic universality limit ${{\omega }_{{\rm c}}}\to \infty$ we have

$$\begin{eqnarray}&&W(\tau )=2K{\rm ln} \left[ \frac{{{\omega }_{{\rm c}}}}{\pi T}{\rm sinh} (\pi T|\tau |) \right]+{\rm i}\;\pi K{\rm sign}(\tau ),\end{eqnarray} \tag{ 34 }$$

where K is a dimensionless damping strength, $K=M{{\gamma }_{1}}q_{0}^{2}/(2\pi )$. Because the kernel ${{\dot{W}}^{\prime\prime }}(\tau )$ is memoryless, ${{\dot{W}}^{\prime\prime }}(\tau )=2\pi K\delta (\tau )$, the sojourn–sojourn correlation ${{\mathcal{U}}_{m,\ell }}(t)$ is restricted to the nearest-neighbour term $\ell =m-1$ in which the intermediate blip has length zero

$$\begin{eqnarray}&&{{\mathcal{U}}_{m,m-1}}=-2\pi K\;\delta ({{t}_{2m}}-{{t}_{2m-1}}).\end{eqnarray} \tag{ 35 }$$

Observing that a blip of zero length does not interact with other blips, the ξ-summation in equation (32) cancels all sojourn–sojourn correlations except those of the first with the second. Thus we have ${{\mathcal{U}}_{m,m-1}}={{\delta }_{m,1}}\;{{\mathcal{U}}_{1,0}}$. In addition, the odd term $\mathcal{R}_{1}^{(1,{\rm a})}(t)$ is zero because the single blip left over has zero length, thus readily

$$\begin{eqnarray}&&{{\langle Q(t)\rangle }_{1}}=(\pi K/2){{\delta }^{2}}\;t.\end{eqnarray} \tag{ 36 }$$

Here $\delta =\Delta \;{{{\rm e}}^{-{{W}_{{\rm adia}}}/2}}$ is a dressed transition amplitude taking into account the polarization cloud of the bath modes in the adiabatic limit. We have ${{W}_{{\rm adia}}}=2K{\rm ln} ({{\omega }_{{\rm c}}}/{{\Delta }_{{\rm r}}})$ and hence $\delta ={{\Delta }_{{\rm r}}}\equiv \Delta {{(\Delta /{{\omega }_{{\rm c}}})}^{K/(1-K)}}$ in the regime $T\ll \Delta$, and ${{W}_{{\rm adia}}}=2K{\rm ln} [{{\omega }_{{\rm c}}}/(2\pi T)]$ and hence $\delta ={{\Delta }_{{\rm r}}}{{(2\pi T/{{\Delta }_{{\rm r}}})}^{K}}$ in the regime $T\gg \Delta$.

Consider next the temperature-dependent term ${{\langle Q(t)\rangle }_{2}}$. We see from the imaginary part of the expression (34) that the blip-sojourn correlations ${{X}_{j,k}}$ are restricted to the nearest-neighbour term, ${{X}_{j,k}}=\pi K\;{{\delta }_{j,k+1}}$. Thus we have ${{\phi }_{m,k}}=\pi K{{\xi }_{k+1}}$. With this the expression (32) with (33) takes the form

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\langle Q(t)\rangle }_{2}} & = & \frac{1}{2}\;\mathop{\sum }\limits_{m=1}^{\infty }{{\left( \frac{-{{\Delta }^{2}}}{2} \right)}^{m}}{\rm cos} {{(\pi K)}^{m}}\int _{0}^{t}{{\mathcal{D}}_{2m}}\{{{t}_{j}}\}\mathop{\sum }\limits_{\{{{\xi }_{j}}=\pm 1\}}{{G}_{m}} \\ {} & {} & \quad \times \ \left( -{\rm cos} ({{\varphi }_{m}}){\rm tan} (\pi K){{\xi }_{1}}+{{p}_{0}}{\rm sin} ({{\varphi }_{m}}) \right)\mathop{\sum }\limits_{\ell =1}^{m}{{\mathcal{V}}_{m,\ell }}(t)\;{{\xi }_{\ell }}. \\ \end{array}\end{eqnarray} \tag{ 37 }$$

At this point it is expedient to turn to the average heat power $\langle P(t)\rangle$, which is the time derivative of the heat, $\langle P(t)\rangle =\langle \dot{Q}(t)\rangle$. Since ${{\mathcal{V}}_{m,\ell }}(t={{t}_{2m}})=0$, the derivative of ${{\langle Q(t)\rangle }_{2}}$ is restricted to the derivative of ${{\mathcal{V}}_{m,\ell }}(t={{t}_{2m}})$ in equation (37). Observing that ${{\ddot{W}}^{\prime }}(\tau )$ coincides with $L^{\prime} (\tau )$ given in equation (16), we readily obtain

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \langle P(t)\rangle & = & \frac{\pi K}{2}{{\delta }^{2}}\;-\;\frac{1}{2}\;\mathop{\sum }\limits_{m=1}^{\infty }{{\left( \frac{-{{\Delta }^{2}}}{2} \right)}^{m}}{\rm cos} {{(\pi K)}^{m}}\int _{0}^{t}{{\mathcal{D}}_{2m}}\{{{t}_{j}}\}\mathop{\sum }\limits_{\{{{\xi }_{j}}=\pm 1\}}{{G}_{m}} \\ {} & {} & \times \left( {\rm cos} ({{\varphi }_{m}}){\rm tan} (\pi K){{\xi }_{1}}-{{p}_{0}}{\rm sin} ({{\varphi }_{m}}) \right)\mathop{\sum }\limits_{\ell =1}^{m}\left[ \;L^{\prime} (t-{{t}_{2\ell -1}})-L^{\prime} (t-{{t}_{2\ell }}) \right]\;{{\xi }_{\ell }}. \\ \end{array}\end{eqnarray} \tag{ 38 }$$

Equation (38) is the exact formal path sum expression for the average heat power of the TSS in the Ohmic scaling limit for any coupling strength K. Except for the last factor, the expression represents the dynamics of the RDM of the TSS under all the quantum noise-induced time-correlations carried by the functional G_m. The residual factor correlates the flip times of the TSS with the final time t. This factor is specific to the mean of the heat power. The expression (38) may be used to calculate the time dependence of the heat power for any time-dependent driving of the system according to equation (26).

To sound out validity and potential of the presented method, we now study the expression (38) for a constant bias $\epsilon _{0}^{{}}$. First, we see that in the term of ${{\Delta }^{2m}}$ the blip ${{\xi }_{\ell }}$, which may take any position in the sequence of m 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 ℓ, ${{\Lambda }_{j,\ell }}=0$ and ${{\Lambda }_{\ell ,k}}=0$, as well as all interblip correlations across this blip, i.e., ${{\Lambda }_{j,k}}=0$, where $j\lt \ell$ and $k\gt \ell$. Then the series (38) can be expressed in terms of a convolution which involves in the time segment the population correlation function $\langle {{\sigma }_{z}}(\tau )\rangle$ [51] and the coherence correlation function $\langle {{\sigma }_{x}}(\tau )\rangle$ of the TSS [37, 62, 63]. The resulting expression is

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \langle P(t)\rangle & = & \frac{\pi K}{2}{{\delta }^{2}}-\frac{\Delta }{2}\int _{0}^{t}{\rm d}\tau \left( \;{{\langle {{\sigma }_{x}}(t-\tau )\rangle }_{{\rm s}}}\;{{\langle {{\sigma }_{z}}(\tau )\rangle }_{{\rm s}}}-{{\langle {{\sigma }_{z}}(t-\tau )\rangle }_{{\rm a}}}\;{{\langle {{\sigma }_{x}}(\tau )\rangle }_{{\rm a}}} \right. \\ {} & {} & -\left. \;{{p}_{0}}\left[ {{\langle {{\sigma }_{x}}(t-\tau )\rangle }_{{\rm a}}}\;{{\langle {{\sigma }_{z}}(\tau )\rangle }_{{\rm s}}}-{{\langle {{\sigma }_{z}}(t-\tau )\rangle }_{{\rm s}}}\;{{\langle {{\sigma }_{x}}(\tau )\rangle }_{{\rm a}}}\; \right]\; \right)\;L^{\prime} (\tau ) \\ {} & {} & -\frac{1}{2}{\rm tan} (\pi K)\int _{0}^{t}{\rm d}\tau \;L^{\prime} (\tau )\;\frac{{\rm d}{{\langle {{\sigma }_{z}}(\tau )\rangle }_{{\rm s}}}}{{\rm d}\tau }. \\ \end{array}\end{eqnarray} \tag{ 39 }$$

Here the index s (a) indicates that the respective term is symmetric (antisymmetric) under bias inversion, $\epsilon _{0}^{{}}\to -\epsilon _{0}^{{}}$ . The initial conditions are ${{\langle {{\sigma }_{z}}(0)\rangle }_{{\rm a}}}={{\langle {{\sigma }_{x}}(0)\rangle }_{{\rm s}}}={{\langle {{\sigma }_{x}}(0)\rangle }_{{\rm a}}}=0$ and ${{\langle {{\sigma }_{z}}(0)\rangle }_{{\rm s}}}=1$.

In the further analysis of the average heat power $\langle P(t)\rangle$ we can thus employ known analytic results for the nonequilibrium correlation functions $\langle {{\sigma }_{z}}(t)\rangle$ and $\langle {{\sigma }_{x}}(t)\rangle$ reported in the recent literature [37, 51, 60].

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\rangle$ and $|R\rangle$ as linear combinations of the ground state $|g\rangle$ and excited state $|e\rangle$ of ${{H}_{{\rm S}}}$, and assuming initial preparation according to equation (27) and thermal occupation of the states $|g\rangle$ and $|e\rangle$ at time infinity, we obtain

$$\begin{eqnarray}&&{{Q}_{\infty }}\equiv \langle Q(t\to \infty )\rangle ={{\langle E\rangle }_{{\rm ini}}}-{{\langle E\rangle }_{{\rm eq}}}=\frac{{{p}_{0}}\;\epsilon _{0}^{{}}}{2}+\frac{\sqrt{{{\delta }^{2}}+\epsilon _{0}^{2}}}{2}{\rm tanh} \left( \frac{\sqrt{{{\delta }^{2}}+\epsilon _{0}^{2}}}{2T} \right),\end{eqnarray} \tag{ 40 }$$

where ${{\langle E\rangle }_{{\rm ini}}}$ and ${{\langle E\rangle }_{{\rm 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.

4.3. Markov regime

In the temperature regime $T\gg \sqrt{{{\delta }^{2}}+\epsilon _{0}^{2}}$, the power spectrum $S(\omega )=K\pi \;\omega {\rm coth} (\frac{\omega }{2T})$ of the collective bath mode $\mathcal{E}(t)={{q}_{0}}{{\sum }_{\alpha }}{{c}_{\alpha }}{{x}_{\alpha }}(t)$ is essentially white, $S(\omega )=\vartheta$, where $\vartheta \equiv 2\pi K\;T$. Then the real parts of the polarization correlation function $L(\tau )$ and the flip-pair correlation function $W(\tau )$ of equation (34) take the form $L^{\prime} (\tau )=\vartheta \;\delta (\tau )$ and $W^{\prime} (\tau )={{W}_{{\rm adia}}}+\vartheta \;|\tau |$. Here ${{W}_{{\rm adia}}}=2K{\rm ln} [{{\omega }_{{\rm c}}}/(2\pi T)]$ is the adiabatic contribution which dresses the transfer matrix element according to $\delta =\Delta \;{{{\rm e}}^{-{{W}_{{\rm adia}}}/2}}={{\Delta }_{{\rm r}}}{{(2\pi T/{{\Delta }_{{\rm r}}})}^{K}}$.

Observing that the function $L^{\prime} (\tau )$ is memory-less and the initial conditions for $\langle {{\sigma }_{z}}(t)\rangle$ and $\langle {{\sigma }_{x}}(t)\rangle$ in equation (39) are $\langle {{\sigma }_{z}}(0)\rangle =1$ and $\langle {{\sigma }_{x}}(0)\rangle =0$, the expression (39) for the heat power reduces to the concise form

$$\begin{eqnarray}&&\langle P(t)\rangle =\frac{\pi }{2}K{{\delta }^{2}}-\pi KT\Delta \;[{{\langle {{\sigma }_{x}}(t)\rangle }_{{\rm s}}}-{{p}_{0}}{{\langle {{\sigma }_{x}}(t)\rangle }_{{\rm a}}}\;]\ .\end{eqnarray} \tag{ 41 }$$

Since the pair interaction $W^{\prime} (\tau )\propto \tau$, the interblip correlations ${{\Lambda }_{j,k}}$ are zero. As a result, the series for the coherence correlation function $\langle {{\sigma }_{x}}(t)\rangle$ is easily summed in Laplace space, yielding [37, 63]

$$\begin{eqnarray}&&\langle {{\hat{\sigma }}_{x}}(\lambda )\rangle =\frac{{{\delta }^{2}}}{\Delta \;\lambda }\;\frac{K\pi [{{\delta }^{2}}+\lambda (\vartheta +\lambda )]+\;\epsilon _{0}^{{}}\lambda }{{{\delta }^{2}}(\vartheta +\lambda )+\lambda [\;\epsilon _{0}^{2}+{{(\vartheta +\lambda )}^{2}}\;]}\qquad {\rm with}\qquad \vartheta =2\pi KT.\end{eqnarray} \tag{ 42 }$$

Importantly, the residuum of the pole at the origin, which is the equilibrium value ${{\langle {{\sigma }_{x}}\rangle }_{{\rm eq}}}={{\delta }^{2}}/(2T\Delta )$ of $\langle {{\sigma }_{x}}(t)\rangle$ reached at time infinity, together with the factor $\pi KT\Delta$ cancels exactly the constant term $\pi K{{\delta }^{2}}/2$ in equation (41). This secures that the heat power in fact vanishes as $t\to \infty$, and the dynamics of $\langle P(t)\rangle$ is determined by the simple zeros ${{\lambda }_{1}}$, ${{\lambda }_{2}}$ and ${{\lambda }_{3}}$ of the denominator in equation (42). The resulting expression is

$$\begin{eqnarray}&&\langle P(t)\rangle =\pi KT{{\delta }^{2}}\left( \frac{{{p}_{0}}\;{{\epsilon }_{0}}{{\lambda }_{1}}-K\pi [\;{{\delta }^{2}}+(\vartheta +{{\lambda }_{1}}){{\lambda }_{1}}\;]}{{{\lambda }_{1}}({{\lambda }_{1}}-{{\lambda }_{2}})({{\lambda }_{1}}-{{\lambda }_{3}})}{{{\rm e}}^{{{\lambda }_{1}}t}}\;+\;{\rm cycl}{\rm .} \right),\end{eqnarray} \tag{ 43 }$$

where 'cycl.' denotes addition of the terms obtained by cyclic permutation of the indices. The behaviours of ${{\lambda }_{1}}$, ${{\lambda }_{2}}$ and ${{\lambda }_{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 ${{\gamma }_{{\rm r}}}$ as ${{\lambda }_{1,2}}=-\gamma \mp {\rm i}\;\Omega$ and ${{\lambda }_{3}}=-{{\gamma }_{{\rm r}}}$.

In the temperature regime $\sqrt{{{\delta }^{2}}+\epsilon _{0}^{2}}\ll T\ll \sqrt{{{\delta }^{2}}+\epsilon _{0}^{2}}/(2\pi K)$ the bath coupling ${{H}_{{\rm I}}}$ is a perturbation. In the leading one-boson exchange approximation, the rates increase linearly with the scaled temperature ϑ,

$$\begin{eqnarray}&&\qquad {{\gamma }_{{\rm r}}}=\frac{{{\delta }^{2}}}{{{\delta }^{2}}+\epsilon _{0}^{2}}\;\vartheta +\mathcal{O}({{\vartheta }^{3}}),\qquad \gamma =\frac{{{\delta }^{2}}+2\epsilon _{0}^{2}}{2({{\delta }^{2}}+\epsilon _{0}^{2})}\;\vartheta +\mathcal{O}({{\vartheta }^{3}})\;,\qquad \Omega =\sqrt{{{\delta }^{2}}+\epsilon _{0}^{2}}+\mathcal{O}({{\vartheta }^{2}}).\end{eqnarray} \tag{ 44 }$$

As temperature is raised, multi-boson exchange processes become increasingly important. In the temperature range well above $\sqrt{{{\delta }^{2}}+\epsilon _{0}^{2}}/(2\pi K)$, the decoherence and relaxation rates show drastically different behaviours

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \qquad {{\gamma }_{{\rm r}}}=\frac{{{\delta }^{2}}}{\vartheta }+\mathcal{O}(1/{{\vartheta }^{3}}),\qquad \gamma =\vartheta -\frac{{{\delta }^{2}}}{\vartheta }+\mathcal{O}(1/{{\vartheta }^{3}}), \\ \Omega ={{\epsilon }_{0}}+\frac{{{\delta }^{2}}(4\epsilon _{0}^{2}-{{\delta }^{2}})}{8\epsilon _{0}^{3}}\frac{1}{{{\vartheta }^{2}}}+\mathcal{O}(1/{{\vartheta }^{4}}). \\ \end{array}\end{eqnarray} \tag{ 45 }$$

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 (43) yields the average heat $\langle Q(t)\rangle$ transferred to the reservoir until time t. Readily the average heat transferred until time infinity is found from equation (43) in the Markovian regime as

$$\begin{eqnarray}&&{{Q}_{\infty }}=\frac{{{p}_{0}}\;\epsilon _{0}^{{}}}{2}+\frac{{{\delta }^{2}}+\epsilon _{0}^{2}}{4T}.\end{eqnarray} \tag{ 46 }$$

It is straightforward to see by use of the Vieta relations for the frequencies ${{\lambda }_{1}}$, ${{\lambda }_{2}}$ and ${{\lambda }_{3}}$, that the expression (46) holds for all T in the Markov regime. Importantly, the result (46) is just the high-temperature limit of the expression (40).

4.4. Quantum noise regime

For temperatures of the order $\sqrt{{{\delta }^{2}}+\epsilon _{0}^{2}}$ or lower, the power spectrum $S(\omega )$ is coloured because of quantum mechanical noise and the functions $L^{\prime} (\tau )$ and $W^{\prime} (\tau )$ carry memory. Nevertheless, it is possible to calculate exactly the one-boson contribution to the rates γ and ${{\gamma }_{{\rm r}}}$ [37, 60]. The resulting expressions are

$$\begin{eqnarray}&&{{\gamma }_{{\rm r}}}=\pi K\frac{{{\delta }^{2}}}{\Omega }{\rm coth} \left( \frac{\Omega }{2T} \right),\qquad {\rm and}\qquad \gamma =\frac{1}{2}{{\gamma }_{{\rm r}}}+2\pi KT\frac{\epsilon _{0}^{2}}{{{\Omega }^{2}}},\end{eqnarray} \tag{ 47 }$$

where $\Omega =\sqrt{{{\delta }^{2}}+\epsilon _{0}^{2}}$ and $\delta ={{\Delta }_{{\rm r}}}$. The correlation functions $\langle {{\sigma }_{z}}(\tau )\rangle$ and $\langle {{\sigma }_{x}}(\tau )\rangle$ going down into the expression of equation (39) take the form

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \langle {{\sigma }_{z}}(\tau )\rangle & = & \frac{\epsilon _{0}^{{}}}{\Omega }{\rm tanh} \left( \frac{\Omega }{2T} \right)\left[ 1-\;{{{\rm e}}^{-{{\gamma }_{{\rm r}}}\tau }} \right]+\frac{\epsilon _{0}^{2}}{{{\Omega }^{2}}}\;{{{\rm e}}^{-{{\gamma }_{{\rm r}}}\tau }}+\frac{{{\delta }^{2}}}{{{\Omega }^{2}}}{\rm cos} (\Omega \tau ){{{\rm e}}^{-\gamma \tau }}\; \\ \langle {{\sigma }_{x}}(\tau )\rangle & = & \frac{{{\delta }^{2}}}{\Delta \Omega }{\rm tanh} \left( \frac{\Omega }{2T} \right)\left[ 1-\;{{{\rm e}}^{-{{\gamma }_{{\rm r}}}\tau }} \right]+\frac{\epsilon _{0}^{{}}{{\delta }^{2}}}{\Delta {{\Omega }^{2}}}\left[ \;{{{\rm e}}^{-{{\gamma }_{{\rm r}}}\tau }}-{\rm cos} (\Omega \tau ){{{\rm e}}^{-\gamma \tau }}\; \right]. \\ \end{array}\end{eqnarray} \tag{ 48 }$$

Next, the convolution in equation (39) is conveniently evaluated by using for $L^{\prime} (\tau )$ its spectral representation in equation (16) and writing the ${\rm coth} (\frac{\omega }{2T})$-function in terms of the sum representation

$$\begin{eqnarray}&&{\rm coth} \left( \frac{\omega }{2T} \right)=\mathop{\sum }\limits_{n=-\infty }^{\infty }\frac{2\omega T}{{{\omega }^{2}}+\omega _{n}^{2}},\end{eqnarray} \tag{ 49 }$$

where ${{\omega }_{n}}=2\pi T\;n$ is a bosonic Matsubara frequency. Then the τ-integration is straightforward. The remaining ω-integral picks up (i) a contribution from the singularity at $\omega =\Omega$ formed in the limit $K\to 0$ and (ii) a contribution from the residua of the infinite sequence of poles stringed along the imaginary axis at $\omega ={\rm i}\;{{\omega }_{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 equation (39), and the remaining part is

$$\begin{eqnarray}&&\qquad {{\langle P(t)\rangle }_{1}}=\frac{\pi }{2}K{{\delta }^{2}}\left( {{{\rm e}}^{-{{\gamma }_{{\rm r}}}t}}+\frac{{{p}_{0}}\;\epsilon _{0}^{{}}}{\Omega }{\rm coth} \left( \frac{\Omega }{2T} \right)\left[ \;{{{\rm e}}^{-{{\gamma }_{{\rm r}}}t}}-{\rm cos} (\Omega t){{{\rm e}}^{-\gamma t}} \right] \right).\end{eqnarray} \tag{ 50 }$$

The contribution from the Matsubara poles is odd in the bias and is

$$\begin{eqnarray}&&\qquad {{\langle P(t)\rangle }_{2}}=K{{\delta }^{2}}\;\frac{{{p}_{0}}\;\epsilon _{0}^{{}}}{\Omega }\left( \left. {{u}_{2}}(t)\left[ {\rm cos} (\Omega t){{{\rm e}}^{-\gamma t}}-{{{\rm e}}^{-{{\gamma }_{{\rm r}}}t}}\; \right]+[{{u}_{1}}(0)-{{u}_{1}}(t)]\;{\rm sin} (\Omega t)\;{{{\rm e}}^{-\gamma t}} \right] \right),\end{eqnarray} \tag{ 51 }$$

where

$$\begin{eqnarray}&&\qquad {{u}_{k}}(t)=\frac{2\pi T}{\Omega }\mathop{\sum }\limits_{n=1}^{\infty }\;\frac{\omega _{n}^{k}}{\omega _{n}^{2}+{{\Omega }^{2}}}{{{\rm e}}^{-{{\omega }_{n}}(t+\frac{1}{{{\omega }_{{\rm c}}}})}},\qquad k=1,2.\end{eqnarray} \tag{ 52 }$$

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 ${{\omega }_{{\rm c}}}$.

As the Markov regime at T well above Ω is approached, the expression (50) smoothly matches on the expression (43) with the one-boson rates (44). The heat power term ${{\langle P(t)\rangle }_{2}}$ is a pure quantum contribution, which is negligibly small in this regime. Conversely, as the temperature is lowered, quantum effects appear in ${{\langle P(t)\rangle }_{1}}$ through the coth-function in the amplitude and in the rate expressions (47). This function indicates emission and absorption of quanta with energy ω in thermal equilibrium. Moreover, the term ${{\langle P(t)\rangle }_{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

$$\begin{eqnarray}&&\qquad {{A}_{{\rm rel}}}(t)=1-\frac{2}{\pi }{\rm tanh} \left( \frac{\Omega }{2T} \right)\;{{u}_{2}}(t),\qquad B(t)=\frac{2}{\pi }{\rm tanh} \left( \frac{\Omega }{2T} \right)\;[\;{{u}_{1}}(0)-{{u}_{1}}(t)].\end{eqnarray} \tag{ 53 }$$

With these forms, the antisymmetric part of the total average heat power $\langle P(t\rangle ={{\langle P(t)\rangle }_{1}}+{{\langle P(t)\rangle }_{2}}$ is conveniently written in terms of the amplitude ${{A}_{{\rm rel}}}(t)$ of the relaxation contribution, the amplitude ${{A}_{{\rm osc}}}(t)$ of the oscillatory contribution, and a phase shift $\varphi (t)$,

$$\begin{eqnarray}&&{{A}_{{\rm osc}}}(t)=\sqrt{{{A}_{{\rm rel}}}{{(t)}^{2}}+B{{(t)}^{2}}},\qquad \quad \varphi (t)={\rm arctan} [B(t)/{{A}_{{\rm rel}}}(t)].\end{eqnarray} \tag{ 54 }$$

The resulting expression for the total mean heat power $\langle P(t)\rangle ={{\langle P(t)\rangle }_{1}}+{{\langle P(t)\rangle }_{2}}$ is

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \langle P(t)\rangle & = & \frac{\pi }{2}K{{\delta }^{2}}\{{\rm ,}{{{\rm e}}^{-{{\gamma }_{{\rm r}}}t}}+\;\frac{{{p}_{0}}\;\epsilon _{0}^{{}}}{\Omega }{\rm coth} \left( \frac{\Omega }{2T} \right)\left( {{A}_{{\rm rel}}}(t)\;{{{\rm e}}^{-{{\gamma }_{{\rm r}}}t}}-{{A}_{{\rm osc}}}(t){\rm sign}[\varphi (t)]\; \right. \\ {} & {} & \left. \times {\rm cos} [\;\Omega t+\varphi (t)]{{{\rm e}}^{-\gamma t}} \right)\}. \\ \end{array}\end{eqnarray} \tag{ 55 }$$

The oscillatory part of $\langle P(t)\rangle$ describes seesaw transport of heat between system and reservoir. The damped oscillatory parts are superimposed by relaxation contributions.

Considering equation (55) and comparing it with equation (50), one can deduce that the degree of deviation from unity of the amplitude factors ${{A}_{{\rm rel}}}(t)$ and ${{A}_{{\rm osc}}}(t)$ indicate the relative strength of the pure quantum contribution ${{\langle P(t)\rangle }_{2}}$. Indeed, in the case ${{A}_{{\rm rel}}}(t)={{A}_{{\rm osc}}}(t)=1$ and $\varphi (t)=0$ one recovers the form (50) in which the pure quantum contribution ${{\langle P(t)\rangle }_{2}}=0$ is absent. The amplitudes ${{A}_{{\rm rel}}}(t)$ and ${{A}_{{\rm osc}}}(t)$, and the phase $\varphi (t)$ are plotted in figures 1 and 2 as a functions of time for three different temperatures. The amplitude ${{A}_{{\rm rel}}}(t)$ approaches unity in the regime $t\gt 1/\Omega$ for $T\ll \Omega$ and in the regime $t\gt 1/T$ for $T\gg \Omega$. For fixed T, the amplitude ${{A}_{{\rm rel}}}(t)$ is monotonically decreasing as time t is lowered and eventually changes sign at time t₀. In the temperature regime $T\gg \Omega$, the instant t₀ depends inversely on T as ${{t}_{0}}\approx {\rm ln} (2)/(2\pi T)$ and eventually enters the core region, ${{t}_{0}}=\mathcal{O}(1/{{\omega }_{{\rm c}}}),$ at temperatures of the order ${{\omega }_{{\rm c}}}$. Reversely, the instant t₀ increases with decreasing temperature and reaches at T = 0 the maximal value ${{t}_{0}}=1/(\pi \Omega )$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The amplitude ${{A}_{{\rm osc}}}(t)$ has a minimum at the time t₀, at which ${{A}_{{\rm rel}}}(t)$ is zero, whereas the phase $\varphi (t)$ undergoes a sudden jump from $-\pi /2$ to $\pi /2$ as t passes t₀ from below. At times below t₀, the pure quantum contribution ${{\langle P(t)\rangle }_{2}}$ dominates the behaviour of the heat power $\langle P(t)\rangle$. Both the amplitudes ${{A}_{{\rm osc}}}(t)$ and the phase $\varphi (t)$ approach constant values ${{A}_{{\rm osc},\infty }}$ and ${{\varphi }_{\infty }}$ at times $t\gg 1/\Omega$. These values get larger as T is decreased, and at low T they are notably different from the values ${{A}_{{\rm osc}}}=1$ and $\varphi =0$ holding in the absence of the pure quantum contribution ${{\langle P(t)\rangle }_{2}}$. The asymptotic values are largest at T = 0 and are

$$\begin{eqnarray}&&{{A}_{{\rm osc},\infty }}=\sqrt{1+B_{\infty }^{2}},\qquad {\rm and}\qquad {{\varphi }_{\infty }}={\rm arctan} ({{B}_{\infty }}),\end{eqnarray} \tag{ 56 }$$

where ${{B}_{\infty }}=(2/\pi )[\;{\rm ln} ({{\omega }_{{\rm c}}}/\Omega )-{{C}_{{\rm E}}}]$, and where ${{C}_{{\rm E}}}$ is Euler's constant. This shows that the direct time dependence of $\langle P(t)\rangle$ and $\langle Q(t)\rangle$ at low T strongly depends on the pure quantum contribution ${{\langle P(t)\rangle }_{2}}$ for all times. These non-Markovian effects could be measured at temperature $T\lt \Omega$, provided that the detector used can resolve the signal with the sufficient accuracy for time $t\lt 1/\gamma ,1/{{\gamma }_{{\rm r}}}$.

On the bottom line, the contribution of ${{\langle P(t)\rangle }_{2}}$ to the overall heat ${{Q}_{\infty }}$ is negligibly small. In the asymptotic weak-damping limit, the total heat transferred to the bath is found from the expression (55) as

$$\begin{eqnarray}&&{{Q}_{\infty }}=\frac{\pi }{2}K{{\delta }^{2}}\left[ \;1+{{p}_{0}}\frac{\epsilon _{0}^{{}}}{\Omega }{\rm coth} \left( \frac{\Omega }{2T} \right)\; \right]\;\frac{1}{{{\gamma }_{{\rm r}}}}.\end{eqnarray} \tag{ 57 }$$

With the form (47) for ${{\gamma }_{{\rm r}}}$, ${{Q}_{\infty }}$ coincides with the expression (40).

Plots of $\langle P(t)\rangle$ and $\langle Q(t)\rangle$ are shown for different temperatures and different initial conditions in figures 3 and 4. The oscillations of $\langle Q(t)\rangle$ are clearly visible for the case $T=0.1\ \delta$ 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\approx 15\ \delta$ and then diverge to reach the respective asymptotic values indicated by the arrows on the right border. We see from the curves in figure 4 (b), which correspond to the initial condition ${{p}_{0}}=-1$, that the heat can be positive or negative, dependent on temperature.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 TSS. We showed that the difference between the dynamics of the heat transfer and the dynamics of the 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.

We acknowledge the support of the MIUR-FIRB2012—Project HybridNanoDev (Grant No.RBFR1236VV), EU FP7/2007—2013 under REA grant agreement no 630925—COHEAT, MIUR-FIRB2013—Project Coca (Grant No. RBFR1379UX), and the COST Action MP1209. UW acknowleges support from the Deutsche Forschungsgemeinschaft through SFB/TRR21.

In this appendix we outline the evaluation of the influence functional ${{\mathcal{F}}^{(\nu )}}[q,q^{\prime} ]$. Here we use the inverse temperature $\beta =1/T$, but return to T when we employ the expressions (A.17) in section 3.

Upon comparing the expressions (9) and (11) the influence functional ${{\mathcal{F}}^{(\nu )}}[q,q^{\prime} ]$ is found to be defined as

$$\begin{eqnarray}&&\qquad {{\mathcal{F}}^{(\nu )}}[q,q^{\prime} ]=\int {\rm d}{{{\boldsymbol{x}} }_{f}}{\rm d}{\boldsymbol{x}} _{f}^{\prime }{\rm d}{{{\boldsymbol{x}} }_{i}}{\rm d}{\boldsymbol{x}} _{i}^{\prime }\langle {{{\boldsymbol{x}} }_{i}}|\rho _{{\rm R}}^{(\nu )}(0)|{\boldsymbol{x}} _{i}^{\prime }\rangle \langle {\boldsymbol{x}} _{f}^{\prime }|{{{\rm e}}^{{\rm i}\nu {{H}_{{\rm R}}}}}|{{{\boldsymbol{x}} }_{f}}\rangle F[q,{{{\boldsymbol{x}} }_{f}},{{{\boldsymbol{x}} }_{i}}]{{F}^{*}}[q^{\prime} ,{\boldsymbol{x}} _{f}^{\prime },{\boldsymbol{x}} _{i}^{\prime }].\end{eqnarray} \tag{ A.1 }$$

Here the amplitude

$$\begin{eqnarray}&&F[q,{{{\boldsymbol{x}} }_{f}},{{{\boldsymbol{x}} }_{i}}]=\int _{{\boldsymbol{x}} (0)={{{\boldsymbol{x}} }_{i}}}^{{\boldsymbol{x}} (t)={{{\boldsymbol{x}} }_{f}}}\mathcal{D}{\boldsymbol{x}} \;{{{\rm e}}^{{\rm i}\,({{S}_{{\rm R}}}[{\boldsymbol{x}} ]+{{S}_{{\rm I}}}[q,{\boldsymbol{x}} ])}}\end{eqnarray} \tag{ A.2 }$$

is the propagator of the bath exposed to the system-bath coupling, ${{S}_{{\rm I}}}[q,{\boldsymbol{x}} ]=-\int {\rm d}\tau {{H}_{{\rm I}}}[q,{\boldsymbol{x}} ]$.

At first, we see that the influence functional ${{\mathcal{F}}^{(\nu )}}[q,q^{\prime} ]$ can be written in the product form

$$\begin{eqnarray}&&{{\mathcal{F}}^{(\nu )}}[q,q^{\prime} ]=\mathop{\prod }\limits_{\alpha =1}^{N}\mathcal{F}_{\alpha }^{(\nu )}[q,q^{\prime} ].\end{eqnarray} \tag{ A.3 }$$

The term $\mathcal{F}_{\alpha }^{(\nu )}[q,q^{\prime} ]$ is the contribution of the bath oscillator α. Dropping the redundant label R we have

$$\begin{eqnarray}&&\mathcal{F}_{\alpha }^{(\nu )}[q,q^{\prime} ]=\int {\rm d}{{x}_{\alpha f}}{\rm d}x_{\alpha f}^{\prime }{\rm d}{{x}_{\alpha i}}{\rm d}x_{\alpha i}^{\prime }\langle {{x}_{\alpha i}}|{{\rho }^{(\nu )}}_{\alpha }(0)|x_{\alpha i}^{\prime }\rangle \langle x_{\alpha f}^{\prime }|{{{\rm e}}^{{\rm i}\nu {{H}_{\alpha }}}}|{{x}_{\alpha f}}\rangle F[q,{{x}_{\alpha f}},{{x}_{\alpha i}}]{{F}^{*}}[q^{\prime} ,x_{\alpha f}^{\prime },x_{\alpha i}^{\prime }].\end{eqnarray} \tag{ A.4 }$$

The matrices of the density operators $\rho _{\alpha }^{(\nu )}(0)={{{\rm e}}^{-(\beta +{\rm i}\,\nu ){{H}_{\alpha }}}}/{{Z}_{\alpha }}$ and ${{{\rm e}}^{{\rm i}\,\nu {{H}_{\alpha }}}}$ depend solely on coordinates of the bath oscillator α and may be written as

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \langle {{x}_{i\alpha }}|\rho _{\alpha }^{(\nu )}(0)|x_{i\alpha }^{\prime }\rangle & = & \frac{1}{{{Z}_{\alpha }}}{{\mathcal{A}}_{1\alpha }}(\beta +{\rm i}\nu ){\rm exp} \left\{ -{{\mathcal{B}}_{1\alpha }}({{x}_{i\alpha }},x_{i\alpha }^{\prime };\beta +{\rm i}\;\nu ) \right\}, \\ \langle x_{f\alpha }^{\prime }|{{{\rm e}}^{{\rm i}\nu {{H}_{\alpha }}}}|{{x}_{f\alpha }}\rangle & = & {{\mathcal{A}}_{1\alpha }}(-{\rm i}\nu ){\rm exp} \left\{ -{{\mathcal{B}}_{1\alpha }}({{x}_{f\alpha }},x_{f\alpha }^{\prime };-{\rm i}\;\nu ) \right\}. \\ \end{array}\end{eqnarray} \tag{ A.5 }$$

The prefactor ${{\mathcal{A}}_{1\alpha }}(\tau )$ and the Euclidean action ${{\mathcal{B}}_{1\alpha }}(x,y;\tau )$ for an imaginary-time interval τ are given by

$$\begin{eqnarray}&&\qquad {{\mathcal{A}}_{1\alpha }}(\tau )={{\left[ \frac{{{m}_{\alpha }}{{\omega }_{\alpha }}}{2\pi {\rm sinh} ({{\omega }_{\alpha }}\tau )} \right]}^{1/2}},\qquad {{\mathcal{B}}_{1\alpha }}(x,y;\tau )=\frac{{{m}_{\alpha }}{{\omega }_{\alpha }}}{2{\rm sinh} ({{\omega }_{\alpha }}\tau )}[({{x}^{2}}+{{y}^{2}}){\rm cosh} ({{\omega }_{\alpha }}\tau )-2xy],\end{eqnarray} \tag{ A.6 }$$

and ${{Z}_{\alpha }}=1/[2{\rm sinh} (\beta {{\omega }_{\alpha }}/2)]$ is the partition function of the reservoir mode α.

Correspondingly, the propagator of the driven oscillator α takes the form

$$\begin{eqnarray}&&\qquad {{F}_{\alpha }}[q,{{x}_{f\alpha }},{{x}_{i\alpha }}]={{\mathcal{A}}_{2\alpha }}(t){\rm exp} \left\{ {\rm i}{{\phi }_{\alpha }}[q,{{x}_{f\alpha }},{{x}_{i\alpha }};t] \right\}\qquad {\rm with}\qquad {{\mathcal{A}}_{2\alpha }}(t)={{\left[ \frac{{{m}_{\alpha }}{{\omega }_{\alpha }}}{2\pi {\rm sin} ({{\omega }_{\alpha }}t)} \right]}^{1/2}}.\end{eqnarray} \tag{ A.7 }$$

The phase is determined by the real-time action in the presence of the external force ${{c}_{\alpha }}q(\tau )$ [37, 61]

$$\begin{eqnarray}&&{{\phi }_{\alpha }}[q,{{x}_{f\alpha }},{{x}_{i\alpha }};t]={{\mathcal{B}}_{2\alpha }}({{x}_{i\alpha }},{{x}_{f\alpha }},t)+{{x}_{i\alpha }}\int _{0}^{t}{\rm d}t^{\prime} \;{{\mathcal{C}}_{\alpha }}(t-t^{\prime} )q(t^{\prime} )+{{x}_{f\alpha }}\int _{0}^{t}{\rm d}t^{\prime} \;{{\mathcal{C}}_{\alpha }}(t^{\prime} )q(t^{\prime} )+{{\psi }_{\alpha }}[q],\end{eqnarray} \tag{ A.8 }$$

where

$$\begin{eqnarray}&&\qquad {{\mathcal{B}}_{2\alpha }}(x,y;t)=\frac{{{m}_{\alpha }}{{\omega }_{\alpha }}}{2{\rm sin} ({{\omega }_{\alpha }}t)}[({{x}^{2}}+{{y}^{2}}){\rm cos} ({{\omega }_{\alpha }}t)-2xy],\qquad {{\mathcal{C}}_{\alpha }}(t^{\prime} )=\frac{{{c}_{\alpha }}}{{\rm sin} ({{\omega }_{\alpha }}t)}{\rm sin} ({{\omega }_{\alpha }}t^{\prime} ).\end{eqnarray} \tag{ A.9 }$$

The first term in equation (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, ${{\psi }_{\alpha }}[q]$ is a global phase which does not depend on the reservoir mode

$$\begin{eqnarray}&&{{\psi }_{\alpha }}[q]=-\frac{{{\mu }_{\alpha }}}{2}\int _{0}^{t}{\rm d}t^{\prime} {{q}^{2}}(t^{\prime} )-\frac{{\rm sin} ({{\omega }_{\alpha }}t)}{{{m}_{\alpha }}{{\omega }_{\alpha }}}\int _{0}^{t}{\rm d}t^{\prime} \int _{0}^{{{t}^{\prime }}}{\rm d}t^{\prime\prime} {{\mathcal{C}}_{\alpha }}(t-t^{\prime} ){{\mathcal{C}}_{\alpha }}(t^{\prime\prime} )q(t^{\prime} )q(t^{\prime\prime} ),\end{eqnarray} \tag{ A.10 }$$

where ${{\mu }_{\alpha }}=c_{\alpha }^{2}/({{m}_{\alpha }}\omega _{\alpha }^{2})$. Here the first term originates from the counter term in the interaction ${{H}_{{\rm I}}}$ given in equation (2). We can thus write the generalized influence functional ${{\mathcal{F}}^{(\nu )}}[q,q^{\prime} ]$ in the form

$$\begin{eqnarray}&&\mathcal{F}_{\alpha }^{(\nu )}[q,q^{\prime} ]=\frac{1}{{{Z}_{\alpha }}}{{\mathcal{A}}_{1\alpha }}(\beta +{\rm i}\;\nu ){{\mathcal{A}}_{1\alpha }}(-{\rm i}\;\nu ){{\mathcal{A}}_{2\alpha }}(t)\mathcal{A}_{2\alpha }^{*}(t)\;{{{\rm e}}^{{\rm i}{{\psi }_{\alpha }}[q]}}\;{{\mathcal{K}}_{\alpha }}[q,q^{\prime} ]\;,\end{eqnarray} \tag{ A.11 }$$

where ${{\mathcal{K}}_{\alpha }}$ captures the integrations of the bath mode α,

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\mathcal{K}}_{\alpha }} & = & \int {\rm d}{{x}_{f\alpha }}{\rm d}x_{f\alpha }^{\prime }{\rm d}{{x}_{i\alpha }}{\rm d}x_{{\rm i}\alpha }^{\prime }{\rm exp} \left\{ i\;[{{\mathcal{B}}_{2\alpha }}({{x}_{i\alpha }},{{x}_{f\alpha }};t)-{{\mathcal{B}}_{2\alpha }}(x_{i\alpha }^{\prime },x_{f\alpha }^{\prime };t)] \right\} \\ {} & {} & \times {\rm exp} \left\{ -[{{\mathcal{B}}_{1\alpha }}({{x}_{i\alpha }},x_{i\alpha }^{\prime };\beta +{\rm i}\;\nu )+{{\mathcal{B}}_{1\alpha }}({{x}_{f\alpha }},x_{f\alpha }^{\prime };-{\rm i}\;\nu )] \right\} \\ {} & {} & \times {\rm exp} \left\{ {\rm i}\int _{0}^{t}{\rm d}t^{\prime} \left[ {{x}_{i\alpha }}{{\mathcal{C}}_{\alpha }}(t-t^{\prime} )+{{x}_{f\alpha }}{{\mathcal{C}}_{\alpha }}(t^{\prime} ) \right]q(t^{\prime} )-[x_{i\alpha }^{\prime }{{\mathcal{C}}_{\alpha }}(t-t^{\prime} )-x_{f\alpha }^{\prime }{{\mathcal{C}}_{\alpha }}(t^{\prime} )]q^{\prime} (t^{\prime} )] \right\}. \\ \end{array}\end{eqnarray} \tag{ A.12 }$$

The exponent of the integrand of the fourfold integral (A.12) is a Gaussian form. Written in terms of the symmetric and antisymmetric cordinates $x_{i\alpha }^{(\pm )}=({{x}_{i,\alpha }}\pm x_{i\alpha }^{\prime })/\sqrt{2}$ and $x_{f\alpha }^{(\pm )}=({{x}_{f\alpha }}\pm x_{f\alpha }^{\prime })/\sqrt{2}$, the Gaussian form is free of mixed terms $x_{i\alpha }^{(+)}x_{f\alpha }^{(+)}$. Hence both the $x_{i\alpha }^{(+)}$- and the $x_{f\alpha }^{(+)}$-integrations may be performed by completing the square

$$\begin{eqnarray}&&\int _{-\infty }^{\infty }{\rm d}x\;{{{\rm e}}^{-a{{x}^{2}}+bx}}=\sqrt{\pi /a}\ {{{\rm e}}^{\frac{{{b}^{2}}}{4a}}}.\end{eqnarray} \tag{ A.13 }$$

With the shift $x_{f\alpha }^{(-)}=y_{f\alpha }^{(-)}+x_{i\alpha }^{(-)}{\rm cos} ({{\omega }_{\alpha }}t)$ the arising quadratic form of the exponent is free of the mixed term $y_{f\alpha }^{(-)}x_{i\alpha }^{(-)}$. Thus the remaining $y_{f\alpha }^{(-)}$- and $x_{i\alpha }^{(-)}$-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 $\eta (\tau )$ and $\xi (\tau )$ given in equation (12). After all, the resulting expression for the influence functional (A.4) may be written in the form

$$\begin{eqnarray}&&\mathcal{F}_{\alpha }^{(\nu )}[\eta ,\xi ]={{{\rm e}}^{-\Phi _{\alpha }^{(0)}[\eta ,\xi ]}}\ {{{\rm e}}^{{\rm i}\,\Delta \Phi _{\alpha }^{(\nu )}[\eta ,\xi ]}}\end{eqnarray} \tag{ A.14 }$$

with

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \Phi _{\alpha }^{(0)}[\eta ,\xi ] & = & {\rm i}\;\frac{{{\mu }_{\alpha }}q_{0}^{2}}{2}\int _{0}^{t}{\rm d}t^{\prime} \eta (t^{\prime} )\xi (t^{\prime} )+\int _{0}^{t}{\rm d}t^{\prime} \\ {} & {} & \times \int _{0}^{{{t}^{\prime }}}{\rm d}t^{\prime\prime} \xi (t^{\prime} )[L_{\alpha }^{\prime }(t^{\prime} -t^{\prime\prime} )\xi (t^{\prime\prime} )+{\rm i}L_{\alpha }^{\prime\prime }(t^{\prime} -t^{\prime\prime} )\eta (t^{\prime\prime} )], \\ \end{array}\end{eqnarray} \tag{ A.15 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \Delta \Phi _{\alpha }^{(\nu )}[\eta ,\xi ] & = & \int _{0}^{t}{\rm d}t^{\prime} \int _{0}^{{{t}^{\prime }}}{\rm d}t^{\prime\prime} \left\{ [\;\eta (t^{\prime} )\eta (t^{\prime\prime} )-\xi (t^{\prime} )\xi (t^{\prime\prime} )]L_{\alpha ,1}^{(\nu )}(t^{\prime} -t^{\prime\prime} ) \right. \\ {} & {} & +\left. [\;\eta (t^{\prime} )\xi (t^{\prime\prime} )-\xi (t^{\prime} )\eta (t^{\prime\prime} )\;]L_{\alpha ,2}^{(\nu )}(t^{\prime} -t^{\prime\prime} ) \right\}. \\ \end{array}\end{eqnarray} \tag{ A.16 }$$

The kernels are

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{L}_{\alpha }}(\tau ) & = & L_{\alpha }^{\prime }(\tau )+{\rm i}L_{\alpha }^{\prime\prime }(\tau )=\frac{c_{\alpha }^{2}q_{0}^{2}}{2{{m}_{\alpha }}{{\omega }_{\alpha }}}[{\rm coth} \frac{\beta {{\omega }_{\alpha }}}{2}{\rm cos} ({{\omega }_{\alpha }}\tau )-{\rm i}\;{\rm sin} ({{\omega }_{\alpha }}\tau )], \\ L_{\alpha ,1}^{(\nu )}(\tau ) & = & \frac{c_{\alpha }^{2}q_{0}^{2}}{2{{m}_{\alpha }}{{\omega }_{\alpha }}}\frac{{\rm sin} (\nu {{\omega }_{\alpha }}/2){\rm sinh} [(\beta +{\rm i}\;\nu ){{\omega }_{\alpha }}/2]}{{\rm sinh} (\beta {{\omega }_{\alpha }}/2)}{\rm cos} ({{\omega }_{\alpha }}\tau ), \\ L_{\alpha ,2}^{(\nu )}(\tau ) & = & \frac{c_{\alpha }^{2}q_{0}^{2}}{2{{m}_{\alpha }}{{\omega }_{\alpha }}}\frac{{\rm sin} (\nu {{\omega }_{\alpha }}/2){\rm cosh} [(\beta +{\rm i}\;\nu ){{\omega }_{\alpha }}/2]}{{\rm sinh} (\beta {{\omega }_{\alpha }}/2)}{\rm sin} ({{\omega }_{\alpha }}\tau ). \\ \end{array}\end{eqnarray} \tag{ A.17 }$$

Here we sketch the derivation of the expression (28). First, the piecewise constant paths with $2m$ flips are

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\eta }^{(m)}}(\tau ) & = & \mathop{\sum }\limits_{j=0}^{m}{{\eta }_{j}}\left[ \;\Theta (\tau -{{t}_{2j}})-\Theta (\tau -{{t}_{2j+1}}) \right], \\ {{\xi }^{(m)}}(\tau ) & = & \mathop{\sum }\limits_{j=1}^{m}{{\xi }_{j}}\left[ \;\Theta (\tau -{{t}_{2j-1}})-\Theta (\tau -{{t}_{2j}})\; \right], \\ \end{array}\end{eqnarray} \tag{ B.1 }$$

where ${{t}_{0}}=0$ and ${{t}_{2m+1}}=t$. Here the initial state of the RDM is ${{\eta }_{i}}={{\eta }_{0}}$ and the final state is ${{\eta }_{f}}={{\eta }_{m}}$.

The path sums in equations (13) and (19) result 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

$$\begin{eqnarray}&&\int \mathcal{D}\eta \int \mathcal{D}\xi \;\cdots \to \mathop{\sum }\limits_{m=0}^{\infty }\int _{0}^{t}{{\mathcal{D}}_{2m}}\{{{t}_{j}}\}\mathop{\sum }\limits_{\{{{\xi }_{j}}=\pm 1\}}\mathop{\sum }\limits_{\{{{\eta }_{j}}=\pm 1\}}\cdots ,\end{eqnarray} \tag{ B.2 }$$

$$\begin{eqnarray}&&\int _{0}^{t}{{\mathcal{D}}_{2m}}\{{{t}_{j}}\}=\int _{0}^{t}{\rm d}{{t}_{2m}}\cdots \int _{0}^{{{t}_{3}}}{\rm d}{{t}_{2}}\int _{0}^{{{t}_{2}}}{\rm d}{{t}_{1}}.\end{eqnarray} \tag{ B.3 }$$

The weight to switch per unit time from the diagonal state ${{\eta }_{j}}$ to the directly following off-diagonal state ${{\xi }_{j+1}}$ is $-{\rm i}\;{{\eta }_{j}}\;{{\xi }_{j+1}}\;\Delta /2$, and the weight to switch per unit time from the offdiagonal state ${{\xi }_{j+1}}$ to the directly following diagonal state ${{\eta }_{j+1}}$ is $-{\rm i}\;{{\xi }_{j+1}}\;{{\eta }_{j+1}}\;\Delta /2$. Thus, the aggregated weight of $2m$ flips depends only on the initial and final sojourn state and is ${{(-1)}^{m}}{{({{\Delta }^{2}}/4)}^{m}}\;{{\eta }_{0}}\;{{\eta }_{m}}$. The weight to stay in a sojourn is unity, whereas the weight to stay in the blip ${{\xi }_{j}}$ during time interval ${{t}_{2j}}-{{t}_{2j-1}}$ is ${\rm exp} [\;{\rm i}\;{{\xi }_{j}}\int _{{{t}_{2j-1}}}^{{{t}_{2j}}}{\rm d}t^{\prime} \epsilon (t^{\prime} )]$. Accordingly, the contributions of m blips are accumulated in the bias weight factor

$$\begin{eqnarray}&&{{B}_{m}}={{{\rm e}}^{{\rm i}\,{{\varphi }_{m}}}},\qquad {{\varphi }_{m}}=\mathop{\sum }\limits_{j=1}^{m}{{\xi }_{j}}\left[ {{\epsilon }_{0}}{{\tau }_{j}}+\int _{{{t}_{2j-1}}}^{{{t}_{2j}}}{\rm d}\tau \;{{\epsilon }_{1}}(\tau ) \right].\end{eqnarray} \tag{ B.4 }$$

For the piecewise constant path (B.1) in the influence action $\Phi _{{\rm cor}}^{(0)}[\eta ,\xi ]$ the time correlations among the $\{{{\eta }_{i}}\}$- and $\{{{\xi }_{j}}\}$- labels or charges are expressed in terms of the second integral of the function $L(\tau )$,

$$\begin{eqnarray}&&W(\tau )=W^{\prime} (\tau )+{\rm i}\;W^{\prime\prime} (\tau )=\frac{q_{0}^{2}}{\pi }\int _{0}^{\infty }{\rm d}\omega \frac{J(\omega )}{{{\omega }^{2}}}\left[ {\rm coth} \left( \frac{\omega }{2T} \right)[\;1-{\rm cos} (\omega \tau )]+{\rm i}\;{\rm sin} (\omega \tau )\; \right].\end{eqnarray} \tag{ B.5 }$$

The resulting correlation factor for $2m$ flips is ${{({{{\rm e}}^{-\Phi _{{\rm cor}}^{(0)}[\eta ,\xi ]}})}_{m}}={{G}_{m}}\;{{H}_{m}}$, where [52, 62, 63]

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{G}_{m}} & = & {\rm exp} \left\{ -\mathop{\sum }\limits_{j=1}^{m}W_{2j,2j-1}^{\prime } \right\}\ {\rm exp} \left\{ -\mathop{\sum }\limits_{j=2}^{m}\mathop{\sum }\limits_{k=1}^{j-1}{{\xi }_{j}}{{\Lambda }_{j,k}}{{\xi }_{k}} \right\}, \\ {{H}_{m}} & = & {\rm exp} \left\{ {\rm i}\;\mathop{\sum }\limits_{k=0}^{m-1}{{\phi }_{m,k}}{{\eta }_{k}} \right\}\qquad {\rm with}\qquad {{\phi }_{m,k}}=\mathop{\sum }\limits_{j=k+1}^{m}{{\xi }_{j}}{{X}_{j,k}}, \\ \end{array}\end{eqnarray} \tag{ B.6 }$$

and where

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\Lambda }_{j,k}} & = & W_{2j,2k-1}^{\prime }+W_{2j-1,2k}^{\prime }-W_{2j,2k}^{\prime }-W_{2j-1,2k-1}^{\prime }, \\ {{X}_{j,k}} & = & W_{2j,2k+1}^{\prime\prime }+W_{2j-1,2k}^{\prime\prime }-W_{2j,2k}^{\prime\prime }-W_{2j-1,2k+1}^{\prime\prime }. \\ \end{array}\end{eqnarray} \tag{ B.7 }$$

Here we used the short form ${{W}_{j,k}}=W({{t}_{j}}-{{t}_{k}})$. The functions ${{\Lambda }_{j,k}}$ and ${{X}_{j,k}}$ describe the interactions between blip k and subsequent blip j, and the interactions between sojourn k and subsequent blip j, respectively. With these individual terms the moment $\langle {{Q}^{n}}(t)\rangle$ takes the form given in equation (28).