Heat-exchange statistics in driven open quantum systems

In small systems, fluctuations of thermodynamic quantities become important. The study of fluctuations in small systems led to the development of stochastic thermodynamics [1, 2], whose predictions, such as the Jarzinski equality [3] and the Crooks fluctuation relations [4], have been experimentally verified in mechanical [5], biological [6], molecular [7, 8], optical [9], and electronic systems [10]. Thermal fluctuations can also be exploited to extract work by using nonequilibrium feedback, allowing for physical realizations of the celebrated Maxwellʼs demon [11–15].

In all the mentioned experiments, the underlying physics is essentially classical. Down to atomic size and low temperatures, however, quantum mechanics is expected to come into play. Thermodynamics and statistical mechanics of quantum systems, as well as different types of quantum fluctuation relations and their observability in different settings, have been the subject of much attention lately [16–22]. However, while experimental demonstrations are still lacking, the very definitions of basic concepts, such as work, are still under debate [16, 23–25]. This difficulty is related to the measurement process: on one hand, work is not by itself an 'observable' in quantum-mechanical sense; on the other hand, its definition in terms of a double-projective measurement, which prevails in literature [16], involves coherence loss between different-energy eigenstates after the first measurement. An interferometric measurement of the work distribution using an ancilla was considered in [20, 21]. Measurement of the environment, rather than the system, has also been proposed, in order to keep track of the energy exchanges between the two [26–28].

Here, we provide a general framework for understanding the statistics of energy exchanges between a periodically driven quantum system and a heat bath. The problem of dissipation in driven systems has a long story [29–32]. However, the statistics of the dissipated heat were considered only recently and either in very specific cases and/or in the weak-drive limit [25, 27, 33]. We set the stage by deriving from microscopic principles a generalized master equation describing the dynamics of the energy exchanges. We focus on a driven two-level system; under a limited set of assumptions, we provide analytical expressions for the mean heat power and for all the cumulants of the energy distribution after a given time. In the long-time limit, we also write down the full probability distribution in a general form. By highlighting the contribution of each individual energy exchange, the full distribution provides us with a deeper insight into the thermalization process.

Within the stated assumptions, our results apply to any type of drive and coupling mechanism. Among different regimes, that of strong driving is particularly interesting, as it cannot be addressed with perturbative techniques—the system and the drive 'hybridize' and thermalization takes place by exchanging photons at dressed-state energies. We illustrate these concepts by discussing two analytically solvable cases in which the same driven quantum system is coupled to the environment in two different ways. As we show, different coupling mechanisms can result in very different energy-exchange processes. Our formalism can be straightforwardly generalized to consider multi-level systems and multiple heat baths; as such, it could be used to study the performance of quantum heat pumps and engines [34–46].

We consider a generic quantum system subject to a periodic drive and weakly coupled to its environment. The environment is a macroscopic object, which we assume to be thermalized at all times. Our quantity of interest is the probability density function (PDF) $P(Q,t)$ for the environment to exchange the amount of energy Q between times 0 and t. The moments of Q are conveniently expressed in terms of the characteristic function (CF) $G(\nu ,t)=\int _{-\infty }^{\infty }dQP(Q,t){{e}^{iQ\nu }}\$, as follows: $\langle {{Q}^{n}}(t)\rangle ={{\left( -i \right)}^{n}}{{\left. \frac{{{d}^{n}}G(\nu ,t)}{d{{\nu }^{n}}} \right|}_{\nu =0}}$.

We write the total Hamiltonian as ${{H}_{T}}(t)=H(t)+{{H}_{R}}+V$, where H(t) and H_R are the system and reservoir Hamiltonians, respectively, and $V=S\otimes B$ is a coupling term consisting of operators S and B acting in the Hilbert space of system and reservoir, respectively. In the following, we will take $\hbar =1$ and denote with ${{\rho }_{T}},$ ρ and ${{\rho }_{R}}$ the total (system + environment), system and environment density operators, respectively.

The PDF for the heat amount Q to be transferred to the reservoir between times ${{t}_{0}}=0$ and t is given by [47, 48]

$$\begin{eqnarray}&&P(Q,t)=\mathop{\sum }\limits_{{{e}_{1}},{{e}_{2}}}\delta ({{e}_{2}}-{{e}_{1}}-Q)p[{{e}_{2}};{{e}_{1}}]p[{{e}_{1}}],\end{eqnarray} \tag{ 1 }$$

where $p[{{e}_{2}};{{e}_{1}}]$ is the conditional probability that a measurement of H_R gives e₂ at time t when it gave e₁ at time t₀ and $p[{{e}_{1}}]$ is the probability to measure e₁ at time t₀. Introducing the projector ${{P}_{{{e}_{j}}}}$ on the j-th state of the reservoir of energy e_j and using the property $P_{{{e}_{j}}}^{2}={{P}_{{{e}_{j}}}}$, we have

$$\begin{eqnarray}\begin{array}{rcl} p[{{e}_{2}},{{e}_{1}}]p[{{e}_{1}}] & = & {\rm Tr}[{{P}_{{{e}_{2}}}}U(t){{P}_{{{e}_{1}}}}{{\rho }_{T}}(0){{P}_{{{e}_{1}}}}{{U}^{\dagger }}(t){{P}_{{{e}_{2}}}}] \\ {} & = & {\rm Tr}[{{U}^{\dagger }}(t){{P}_{{{e}_{2}}}}U(t){{P}_{{{e}_{1}}}}{{\rho }_{T}}(0){{P}_{{{e}_{1}}}}]. \\ \end{array}\end{eqnarray} \tag{ 2 }$$

where U(t) is the evolution operator generated by H_T(t). The CF is given by [48]

$$\begin{eqnarray}&&G(\nu ,t)\equiv \int _{-\infty }^{\infty }dQP(Q,t){{e}^{iQ\nu }}=\mathop{\sum }\limits_{{{e}_{1}},{{e}_{2}}}p[{{e}_{2}};{{e}_{1}}]p[{{e}_{1}}]{{e}^{i\nu ({{e}_{2}}-{{e}_{1}})}}.\end{eqnarray} \tag{ 3 }$$

We assume the initial total density matrix ${{\rho }_{T}}(0)=\rho (0)\otimes {{\rho }_{R}}(0)$ to be factorized into the system density matrix $\rho (0)$ and the thermalized environment density matrix ${{\rho }_{R}}(0)={{e}^{-\beta {{H}_{R}}}}/{{Z}_{R}}$, where Z_R is the partition function of the environment. In particular, this assumption implies that all projectors ${{P}_{{{e}_{j}}}}$ commute with $\rho (0)$, so that the initial measurement of the observable of interest (here, H_R) does not change the subsequent dynamics [25, 48].

After noticing that ${{\sum }_{{{e}_{j}}}}{{P}_{{{e}_{j}}}}{{e}^{\pm i\nu {{e}_{j}}}}={{e}^{\pm i\nu {{H}_{R}}}}$, we write (3) as

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

and finally as

$$\begin{eqnarray}&&G(\nu ,t)={\rm Tr}[\rho _{T}^{\nu }(t)],\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&\rho _{T}^{\nu }(t)={{U}_{\nu /2}}(t){{\rho }_{T}}(0)U_{-\nu /2}^{\dagger }(t)\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&{{U}_{\nu }}(t)={{e}^{i\nu {{H}_{R}}}}U(t){{e}^{-i\nu {{H}_{R}}}}\ \end{eqnarray} \tag{ 7 }$$

satisfies the equation of motion $id{{U}_{\nu }}(t)/dt={{H}_{\nu }}(t){{U}_{\nu }}(t)$, with ${{H}_{\nu }}(t)={{e}^{i\nu {{H}_{R}}}}{{H}_{T}}(t){{e}^{-i\nu {{H}_{R}}}}$. The equation for $\rho _{T}^{\nu }(t)$ can be compared to similar expressions obtained from the full-counting statistics theory of charge transport [49], with ν playing the role of the counting field. In general, the evolution of $\rho _{T}^{\nu }$ is not unitary for $\nu \ne 0$.

For a periodic drive, $H(t+\tau )=H(t)$, where τ is the drive period and $\Omega =2\pi /\tau$ the corresponding angular frequency. We assume that the unitary dynamics generated by H are known; this knowledge is captured in our formalism by the use of Floquet states and quasi-energies [31]. Floquet states are particular solutions of the Schrödinger equation of the form $|{{\Psi }_{\alpha }}(t)\rangle ={{e}^{-i{{\epsilon }_{\alpha }}t}}|{{\phi }_{\alpha }}(t)\rangle$ where the Floquet mode $|{{\phi }_{\alpha }}(t)\rangle$ satisfies $|{{\phi }_{\alpha }}(t+\tau )\rangle =|{{\phi }_{\alpha }}(t)\rangle$ and ${{\epsilon }_{\alpha }}$ is its corresponding quasi-energy.

The derivation of the GQME for the generalized density matrix of the system ${{\rho }^{\nu }}={\rm T}{{{\rm r}}_{R}}[\rho _{T}^{\nu }]$ is done in complete analogy to the standard master equation [32] and is presented in appendix A. We assume weak coupling between the system and the environment and a fast autocorrelation time of the environment (Born–Markov approximation). We model the environment as a collection of harmonic oscillators with bosonic occupation ${{n}_{B}}(\omega )$ and Ohmic spectral density $J(\omega )=\eta \omega$ with η as a dimensionless coupling coefficient. The so-obtained GQME can be simplified in different ways, depending on the time scales of the problem. In the fast-driving regime [31, 50–52], one can safely neglect oscillating terms of the form ${{e}^{ik\Omega t}}$, with $k\ne 0$, and, in most instances, also of the form ${{e}^{i({{\epsilon }_{\alpha }}-{{\epsilon }_{\beta }})t}}$, with $\alpha \ne \beta$. This amounts to a full secular approximation [52], which admits an analytical solution.

We write the resulting GQME in the Schrödinger picture and in the Floquet basis. Under the stated assumptions, the GQME is Markovian and time-independent. The populations $\rho _{\alpha \alpha }^{\nu }=\langle {{\Psi }_{\alpha }}|{{\rho }^{\nu }}|{{\Psi }_{\alpha }}\rangle$ are decoupled from the coherences and satisfy a vector equation of the form $\mathop{\overrightarrow{{{\rho }^{\nu }}}}\limits^{.}=\mathcal{A}\cdot \overrightarrow{{{\rho }^{\nu }}}$. For a two-level system, $\overrightarrow{{{\rho }^{\nu }}}=\{\rho _{11}^{\nu },\rho _{22}^{\nu }\}$ and

$$\begin{eqnarray}\mathcal{A}=\left( \begin{array}{ccccccccccccccc} A_{11}^{\nu }-{{A}_{11}}-{{A}_{21}} & A_{12}^{\nu } \\ A_{21}^{\nu } & A_{22}^{\nu }-{{A}_{22}}-{{A}_{12}} \\ \end{array} \right).\end{eqnarray} \tag{ 8 }$$

A derivation of (8) can be found in appendix B. To define the coefficients in (8), let us first write the matrix elements of the coupling operator S in the Floquet basis [31] as $\langle {{\Psi }_{\alpha }}(t)|S|{{\Psi }_{\beta }}(t)\rangle ={{\sum }_{k}}{{S}_{\alpha \beta ,k}}{{e}^{i{{\Delta }_{\alpha \beta ,k}}t}}$ where ${{S}_{\alpha \beta ,k}}=\frac{1}{\tau }\int _{0}^{\tau }dt\langle {{\phi }_{\alpha }}(t)|S|{{\phi }_{\beta }}(t)\rangle {{e}^{-ik\Omega t}}$ and ${{\Delta }_{\alpha \beta ,k}}={{\epsilon }_{\alpha }}-{{\epsilon }_{\beta }}+k\Omega$. The quantities ${{S}_{\alpha \beta ,k}}$ and ${{\Delta }_{\alpha \beta ,k}}$ can be thought of as coupling amplitudes and transition energies; the corresponding (partial) transition rates are ${{a}_{\alpha \beta ,k}}=|{{S}_{\alpha \beta ,k}}{{|}^{2}}\;s({{\Delta }_{\alpha \beta ,k}})$, where $s(\epsilon )=\theta (\epsilon )J(\epsilon ){{n}_{B}}(\epsilon )+\theta (-\epsilon )J(-\epsilon )\left[ {{n}_{B}}(-\epsilon )+1 \right]$ and $\theta (\epsilon )$ is the Heaviside step function. The rates appearing in (8) are obtained by summing the partial rates over all k contributions; namely, $A_{\alpha \beta }^{\nu }={{\sum }_{k}}{{e}^{-i\nu {{\Delta }_{\alpha \beta ,k}}}}{{a}_{\alpha \beta ,k}}$ and ${{A}_{\alpha \beta }}\equiv A_{\alpha \beta }^{\nu =0}$. Taking the limit $\nu \to 0$ in (8), one recovers a master equation found in the literature, whose dynamic steady state (DSS, also known as quasi-stationary state) is characterized by ${{z}^{{\rm st}}}\equiv \rho _{22}^{{\rm st}}-\rho _{11}^{{\rm st}}=({{A}_{21}}-{{A}_{12}})/({{A}_{12}}+{{A}_{21}})$ [31].

Even before solving the GQME, we notice that the knowledge of ${{\mathop{{\vec{\rho }}}\limits^{.}}^{\nu }}$ can be used to evaluate the mean heat power $\langle \dot{Q}\rangle =d\langle Q\rangle /dt$ transferred to the reservoir at any time (see appendix C):

$$\begin{eqnarray}&&\langle \dot{Q}\rangle =-\mathop{\sum }\limits_{\alpha ,\beta ,k}{{\Delta }_{\alpha \beta ,k}}\ {{a}_{\alpha \beta ,k}}\ {{\rho }_{\beta \beta }}.\end{eqnarray} \tag{ 9 }$$

An equivalent result was recently derived using the standard master equation approach and associating the relaxation processes to the dissipated heat [44, 53, 54]. Equation (9) can be used to classify open driven quantum systems into two categories: those that exchange heat at steady state, and those that do not. In general and in contrast to the undriven case, a driven system will keep exchanging energy with the environment unless its dressed (Floquet) states are decoupled (for instance, by symmetry) from the noise.

The GQME defined by (8) has an analytical solution. The corresponding CF can be written as

$$\begin{eqnarray}&&G(\nu ,t)=c_{-}^{\nu }{{e}^{{{\xi }_{-}}(\nu )t}}+c_{+}^{\nu }{{e}^{{{\xi }_{+}}(\nu )t}},\end{eqnarray} \tag{ 10 }$$

where ${{\xi }_{\pm }}(\nu )$ are the eigenvalues of $\mathcal{A}$ and $c_{\pm }^{\nu }$ the projection of the initial density matrix $\rho (0)$ onto the corresponding eigenvectors $v_{\pm }^{\nu }$, normalized so that ${\rm Tr}v_{\pm }^{\nu }=1$. Many properties of the heat distribution can be obtained from (10) ; see appendix D. We first notice that the first few moments of the distribution can be straightforwardly calculated; in a rigorously defined long-time limit, the cumulants ${{\langle \langle Q\rangle \rangle }_{n}}$ of the PDF [55] are all linear in time and solely determined by the dominant eigenvalue ${{\xi }_{+}}(0)$:

$$\begin{eqnarray}&&{{\langle \langle Q\rangle \rangle }_{n}}={{\left( -i \right)}^{n}}\xi _{+}^{(n)}(0)t,\end{eqnarray} \tag{ 11 }$$

where $\xi _{+}^{(n)}$ is the n-th derivative of ξ₊ with respect to ν. In the following, we will consider the first three cumulants, which coincide with the central moments, namely, mean value, variance, and skewness of the PDF.

Furthermore, the explicit knowledge of the CF allows us to retrieve the full PDF at any given time. This requires inverting the Fourier transform, a task which, in general, must be performed numerically. However, some analytical insight can be gained in the long-time limit. We first notice that $\Re {\rm e}\left[ {{\xi }_{-}}(\nu ) \right]$ is always negative, while $\Re {\rm e}\left[ {{\xi }_{+}}(\nu ) \right]$ is negative everywhere except at ${{\nu }_{n}}=n\tau$, $n\in \mathcal{Z}$, where it vanishes. From (10), we thus see that ${{{\rm lim} }_{t\to \infty }}G(\nu ,t)=0$ for $\nu \ne {{\nu }_{n}}$. This property suggests performing a series expansion around each ν_n. After the expansion, the Fourier transform can be inverted analytically. In the following, we expand ${{\xi }_{+}}(\nu )$ up to the second order. This expansion results in a Gaussian approximation for the PDF; accordingly, it correctly reproduces the first two moments of the original distribution. The PDF reads

$$\begin{eqnarray}\begin{array}{rcl} P(Q,t) & = & w(Q,t)\Omega \mathop{\sum }\limits_{k}[{{p}_{\downarrow }}\delta (Q+{{\epsilon }_{1}}-{{\epsilon }_{2}}+k\Omega ) \\ {} & {} & +\;{{p}_{\uparrow }}\delta (Q+{{\epsilon }_{2}}-{{\epsilon }_{1}}+k\Omega ) \\ {} & {} & +(1-{{p}_{\uparrow }}-{{p}_{\downarrow }})\delta (Q+k\Omega )], \\ \end{array}\end{eqnarray} \tag{ 12 }$$

where ${{p}_{\uparrow }}={{A}_{21}}/({{A}_{21}}+{{A}_{12}}){{\rho }_{11}}(0)$, ${{p}_{\downarrow }}={{A}_{12}}/({{A}_{21}}+{{A}_{12}}){{\rho }_{22}}(0)$ and

$$\begin{eqnarray}&&w(Q,t)=\frac{1}{\sqrt{2\pi }}\frac{1}{\sqrt{2bt}\Omega }{\rm exp} \left[ -\frac{{{\left( Q-at\Omega \right)}^{2}}}{4bt{{\Omega }^{2}}} \right],\end{eqnarray} \tag{ 13 }$$

with real coefficients $a=-i\xi _{+}^{(1)}(0)\geqslant 0$ and $b=-1/2\xi _{+}^{(2)}(0)\geqslant 0$.

Equations (12) and (13) allow for an insightful interpretation of the energy-exchange process. The delta functions account for the fact that the exchanges take place only in multiples of ${{\Delta }_{\alpha \beta ,k}}$. This discretization results from energy being available in photons of energy Ω, the drive frequency, and ${{\epsilon }_{2}}-{{\epsilon }_{1}}$, the dressed energy gap of the driven system. Furthermore, they tell us that the system is either found in the same Floquet state or has undergone a transition upwards (downwards), with probability $1-{{p}_{\uparrow }}-{{p}_{\downarrow }}$ and ${{p}_{\uparrow }}$ $({{p}_{\downarrow }})$, respectively. On the other hand, the total probability that a certain amount of energy has been exchanged is dictated by the weight function $w(Q,t)$. An appealing feature of (12) is that it clearly shows how the heat-exchange distribution 'builds up' on individual exchanges of well-defined energy. One should keep in mind, however, that it was derived under the assumption that $G(\nu ,t)$ is localized at the nodes $\nu =n\tau$. This assumption holds true only in the limit where many energy quanta have been exchanged. For very long times, $t\gg 1/b$, one can neglect the discretization due to the Dirac combs in (12) and find that $P(Q,t)\approx w(Q,t)$.

The results discussed above are general and can be applied to a two-level system with any drive and coupling operator to the environment. As an illustrative, analytically tractable example, we discuss a two-level system interacting with a monochromatic drive, described by the Rabi Hamiltonian $H(t)=\omega /2(1+{{\sigma }_{z}})+g[{{e}^{i(\Omega t-\varphi )}}{{\sigma }_{+}}+{{e}^{-i(\Omega t-\varphi )}}{{\sigma }_{-}}]$ where σ_i are the usual Pauli operators, ω is the bare energy gap, g the amplitude and φ the phase of the driving field. The Floquet states can be written in analytic form (see appendix E) and the quasi-energy gap ${{\epsilon }_{2}}-{{\epsilon }_{1}}$ equals the Rabi frequency ${{\Omega }_{R}}=\sqrt{{{\Delta }^{2}}+4{{g}^{2}}}$, where $\Delta =\Omega -\omega$ is the drive detuning. We consider two types of coupling operators, $S={{\sigma }_{x}}$ (transverse coupling) and $S={{\sigma }_{z}}$ (longitudinal coupling). Both can be realized, e. g., in solid-state devices [56, 57].

When $s={{\sigma }_{x}}$, (9) tells that the DSS heat current is, in general, nonvanishing. In physical terms, the system is continuously 'pumped' by the drive and emits photons to the environment, resulting in a net heat flow to the environment. In this case, the dependence of the PDF on the initial state is negligible at long times. Therefore, we directly take the DSS as the initial state, for which ${{p}_{\downarrow }}={{p}_{\uparrow }}$. The full PDF is shown in the main panel of figure 1 at times ${{t}_{1}}/\tau =80$ (purple) and ${{t}_{2}}/\tau =700$ (blue). It was obtained from the analytical solution (10) by numerically inverting the Fourier transform. The structure of the PDF is the same as in (12): a series of Dirac combs modulated by an envelope which moves in time. As best seen in the Inset, the spectrum of possible energies is composed by symmetric triplets centered at integer multiples of Ω and spaced by an amount ${{\Omega }_{R}}$. There is a strong analogy between this result and the Mollow triplet [58] observed in quantum-optics experiments [59, 60]. As for the envelope function of the PDF, it is manifestly non-Gaussian at early times. In the long-time limit, the PDF tends to a Gaussian as the relative importance of higher-order cumulants decreases [see equation (11)].

Zoom In Zoom Out Reset image size

In figures 1(b) and (c), we show how the PDF is affected by changes in the drive frequency and the environmental temperature. We take a long time ($t/\tau =700$) and use the analytical results (13) and (11). In panel (b), we plot the Gaussian envelope function $w(Q,t)$ for three representative cases. In passing from low to high temperature, the mean value of the distribution is barely affected, while the variance strongly increases. The fact that temperature has little influence on the average dissipated heat is related to the fact that—differently from an undriven system—the transition rates ${{A}_{\alpha \beta }}$ do not satisfy the detailed balance. On the contrary, the variance of the distribution depends on temperature, due to the fact that a higher-temperature environment leads to stronger noise effects. In passing from a resonant to a red-detuned drive, the PDF shows reduced average and variance because of the reduction of the energy injected in the system that can then be dissipated. This trend is confirmed by the behavior of the central moments of the PDF as a function of Δ, shown in panel (c). The maxima of the moments are reached at resonance, i.e., when $\Delta =0$.

We next consider the noise operator $S={{\sigma }_{z}}$. In this case, (9) tells that the mean heat power vanishes at DSS. Therefore, the heat depends critically on the initial state since it is dissipated before the reaching of the DSS. The absence of a net heat flow at DSS is due to the symmetry of the problem, which forbids transitions with energy exchange $n\Omega$ and thereby causes the drive and the environment to be effectively decoupled [54]. The resulting GQME is formally equivalent to that of an undriven system; its derivation can be found in appendix F. The corresponding PDF reads

$$\begin{eqnarray}&&P(Q,t)={{p}_{\downarrow }}(t)\delta (Q-{{\Omega }_{R}})+{{p}_{\uparrow }}(t)\delta (Q+{{\Omega }_{R}})+{{p}_{0}}(t)\delta (Q),\end{eqnarray} \tag{ 14 }$$

where ${{p}_{\downarrow }}(t)={{\Gamma }_{-}}/\Gamma \left( 1-{{e}^{-\Gamma t}} \right){{\rho }_{22}}(0)$, ${{p}_{\uparrow }}(t)={{\Gamma }_{+}}/\Gamma \left( 1-{{e}^{-\Gamma t}} \right){{\rho }_{11}}(0)$ and ${{p}_{0}}(t)=1-{{p}_{\downarrow }}(t)-{{p}_{\uparrow }}(t)$, ${{\Gamma }_{\pm }}=2\pi S_{12,0}^{2}\;s(\pm {{\Omega }_{R}})$ and $\Gamma ={{\Gamma }_{+}}+{{\Gamma }_{-}}$. Notice that the structure of this PDF is the same as the more general one (12), with weight function $w(E,t)=1$ for $-\frac{\Omega }{2}\lt E\lt \frac{\Omega }{2}$ and 0 otherwise. Figure 2 (a) shows a typical time evolution of the PDF as the system reaches its DSS. Figure 2 (b) shows the dependence of the first three central moments on the detuning Δ, starting from the ground state of the undriven system. If, by contrast, the initial state is in a coherent superposition of the undriven system eigenstates $|0\rangle$ and $|1\rangle$, then the drive phase φ also affects the heat distribution, depending on the projection of the initial state onto the $\widehat{xy}$ plane. An example is shown in figure 2 (c) for an equal superposition of $|0\rangle$ and $|1\rangle$. Notably, there exists a range of values where the mean exchanged heat becomes negative, meaning that the environment is more likely to emit energy and the system to absorb it.

Zoom In Zoom Out Reset image size

We have developed a general formalism to study the distribution of heat exchanges between a periodically driven quantum system and a heat bath. Our approach, based on a combination of Floquet theory and a generalized quantum master equation, fully takes into account the effects of the drive. Its power is epitomized by the possibility of obtaining the full heat-exchange distribution, as in equations (12) and (13).

Some of the assumptions made in this work could be further relaxed. The full secular approximation may be replaced by a partial secular approximation, which gives more accurate results in the vicinity of quasi-energy crossings [31, 52]. Several techniques developed in the context of full-counting statistics may be readily adapted to the current setting, including weak-coupling Markovian approaches [61, 62], non-Markovian perturbative expansion [63], recursive schemes [64, 65] and waiting-time analysis [66–68]. These techniques, originally developed for undriven systems, can find application here as the Floquet picture turns a time-dependent problem into a time-independent one.

Our predictions could be tested in a variety of physical systems; for example, superconducting quantum bits embedded in a resistive environment [25, 26] and/or integrated in a circuit-quantum-electrodynamics architecture [69, 70]. Their experimental confirmation would improve our current understanding of the thermodynamics of driven quantum systems and may as well open new avenues for the realization of thermal machines based on quantum mechanics.

We would like to thank J. Pekola and R. Fazio for a careful reading of the manuscript. This work has received funding from the European Union FP7/2007–2013 under REA grant agreement no 630925—COHEAT, from MIUR-FIRB2012—Project HybridNanoDev (Grant No. RBFR1236VV) and from MIUR-FIRB2013—Project Coca (Grant No. RBFR1379UX). S.G. acknowledges financial support from the Finnish National Graduate School in Nanoscience (NGS-NANO) and the Network in Condensed Matter and Materials Physics (CMMP).

Starting from the definition of $G(\nu ,t)$ discussed in the main text, we can write a generalized quantum master equation (GQME) for the system generalized density matrix ${{\rho }^{\nu }}={\rm T}{{{\rm r}}_{R}}[\rho _{T}^{\nu }]$, in analogy with the usual quantum master equation [32]. The latter is recovered in the limit $\nu =0$.

We start from the equation of motion for the total generalized density operator $\rho _{T}^{\nu }(t)$ in (6)

$$\begin{eqnarray}&&{{\mathop{{{\rho }_{T}}}\limits^{.}}^{\nu }}(t)=-i[H(t)+{{H}_{R}},\rho _{T}^{\nu }]-i({{V}_{\frac{\nu }{2}}}\rho _{T}^{\nu }-\rho _{T}^{\nu }{{V}_{-\frac{\nu }{2}}}),\end{eqnarray} \tag{ A.1 }$$

where ${{V}_{\nu }}={{e}^{i\nu {{H}_{R}}}}V{{e}^{-i\nu {{H}_{R}}}}=S\otimes {{B}_{\nu }}=S\otimes {{e}^{i\nu {{H}_{R}}}}B{{e}^{-i\nu {{H}_{R}}}}$. We work in the interaction picture defined by $H(t)+{{H}_{R}}$, so that a generic operator O is transformed as ${{O}^{(I)}}={{e}^{i{{H}_{R}}t}}U_{0}^{\dagger }(t)O{{U}_{0}}(t){{e}^{-i{{H}_{R}}t}}$, where ${{U}_{0}}(t)$ is the evolution operator generated by H(t). After a simplification, we get

$$\begin{eqnarray}&&\dot{\rho }_{T}^{\nu ,(I)}(t)=-i(V_{\frac{\nu }{2}}^{(I)}\rho _{T}^{\nu ,(I)}-\rho _{T}^{\nu ,(I)}V_{-\frac{\nu }{2}}^{(I)}).\end{eqnarray} \tag{ A.2 }$$

We then perform a series expansion up to the second order in the system-environment coupling and invoke the Born approximation, which implies that $\rho _{T}^{(I)}(t)={{\rho }^{(I)}}(t)\otimes {{\rho }_{R}}(0)$ [32]. We introduce a generalized correlation function

$$\begin{eqnarray}&&{{g}_{\nu }}(t^{\prime} -t^{\prime\prime} )={\rm T}{{{\rm r}}_{R}}\left[ B_{-\frac{\nu }{2}}^{(I)}(t^{\prime} )\rho _{R}^{(I)}B_{\frac{\nu }{2}}^{(I)}(t^{\prime\prime} ) \right].\end{eqnarray} \tag{ A.3 }$$

Its explicit evaluation gives

$$\begin{eqnarray}&&{{g}_{\nu }}(\tau )=\mathop{\sum }\limits_{n,n^{\prime} }{{e}^{i({{e}_{n}}-{{e}_{n^{\prime} }})(\tau +\nu )}}{{B}_{nn^{\prime} }}{{B}_{n^{\prime} n}}{{p}_{n}},\end{eqnarray} \tag{ A.4 }$$

where e_n are the eigenenergies of the bath and ${{p}_{n}}={{e}^{-\beta {{e}_{n}}}}/{{Z}_{R}}$ the corresponding occupation probabilities. From (A.4) we see that

$$\begin{eqnarray}&&{{g}_{\nu }}(\tau )={{g}_{0}}(\tau +\nu )\equiv g(\tau +\nu ),\end{eqnarray} \tag{ A.5 }$$

implying that the generalized correlation function ${{g}_{\nu }}$ is a time-shifted ordinary correlation function. For a bosonic bath, the correlation function takes the form

$$\begin{eqnarray}\begin{array}{rcl} g(\tau ) & = & \int _{0}^{\infty }d\omega J(\omega )\left[ {{e}^{i\omega \tau }}{{n}_{B}}(\omega )+{{e}^{-i\omega \tau }}\left( {{n}_{B}}(\omega )+1 \right) \right] \\ {} & = & \int _{0}^{\infty }d\omega J(\omega )\left[ {\rm cos} (\omega \tau ){\rm coth} (\beta \omega /2)-i{\rm sin} (\omega \tau ) \right]\ , \\ \end{array}\end{eqnarray} \tag{ A.6 }$$

where $J(\omega )$ is the spectral density associated to the coupling operator B and ${{n}_{B}}(\omega )$ the Bose distribution function. Performing the Born–Markov approximation and extending the integration times to infinity, we obtain

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial }{\partial t}{{\rho }^{\nu ,(I)}} & = & -\int _{0}^{\infty }d\tau \\ \end{array}\end{eqnarray} \tag{ A.7 }$$

Sometimes it is convenient to go back to the Schrödinger picture. In this case we have:

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial }{\partial t}{{\rho }^{\nu }} & = & -i\left[ H,{{\rho }^{\nu }} \right]-\int _{0}^{\infty }d\tau \\ {} & {} & g(\tau )S(0)S(-\tau ){{\rho }^{\nu }} \\ {} & {} & -g(-\nu +\tau )S(-\tau ){{\rho }^{\nu }}S(0) \\ {} & {} & -g(-\nu -\tau )S(0){{\rho }^{\nu }}S(-\tau ) \\ {} & {} & +g(-\tau ){{\rho }^{\nu }}S(-\tau )S(0). \\ \end{array}\end{eqnarray} \tag{ A.8 }$$

where $S(-\tau )$ stands for $U_{0}^{\dagger }(t-\tau ,t)S{{U}_{0}}(t-\tau ,t)$.

We consider the case of a periodic drive, so that $H(t+\tau )=H(t)$ with $\tau =2\pi /\Omega$. The dynamics of the driven system are captured by considering the Floquet states $|{{\Psi }_{\alpha }}(t)\rangle$ and their corresponding quasi-energies ${{\epsilon }_{\alpha }}$. We also introduce the Floquet modes $|{{\phi }_{\alpha }}(t)\rangle$, satisfying $|{{\phi }_{\alpha }}(t+\tau )\rangle =|{{\phi }_{\alpha }}(t)\rangle$ and are related to the Floquet states by $|{{\Psi }_{\alpha }}(t)\rangle ={{e}^{-i{{\epsilon }_{\alpha }}t}}|{{\phi }_{\alpha }}(t)\rangle$.

The evolution operator of the system in the Floquet basis reads ${{U}_{0}}={{\sum }_{\alpha }}|{{\Psi }_{\alpha }}(t)\rangle \langle {{\Psi }_{\alpha }}(0)|$. We work in the interaction picture defined by ${{U}_{0}}(t)$. Then,

$$\begin{eqnarray}\begin{array}{rcl} {{S}^{(I)}}(t) & = & \mathop{\sum }\limits_{\alpha ,\beta }|{{\Psi }_{\alpha }}(0)\rangle \langle {{\Psi }_{\alpha }}(t)|S|{{\Psi }_{\beta }}(t)\rangle \langle {{\Psi }_{\beta }}(0)| \\ {} & = & \mathop{\sum }\limits_{\alpha ,\beta ,k}{{S}_{\alpha \beta ,k}}{{e}^{i{{\Delta }_{\alpha \beta ,k}}t}}|{{\Psi }_{\alpha }}(0)\rangle \langle {{\Psi }_{\beta }}(0)|, \\ \end{array}\end{eqnarray} \tag{ B.1 }$$

with

$$\begin{eqnarray}&&{{S}_{\alpha \beta ,k}}=\frac{1}{\tau }\int _{0}^{\tau }dt\langle {{\phi }_{\alpha }}(t)|S|{{\phi }_{\beta }}(t)\rangle {{e}^{-ik\Omega t}},\end{eqnarray} \tag{ B.2 }$$

$$\begin{eqnarray}&&{{\Delta }_{\alpha \beta ,k}}={{\epsilon }_{\alpha }}-{{\epsilon }_{\beta }}+k\Omega \ .\end{eqnarray} \tag{ B.3 }$$

Notice the property ${{S}_{\alpha \beta ,k}}=S_{\beta \alpha ,-k}^{*}$.

By substituting (B.2) into (A.7), and denoting with $\rho _{\alpha \beta }^{\nu }=\langle {{\Psi }_{\alpha }}|{{\rho }^{\nu }}|{{\Psi }_{\beta }}\rangle$, we get

$$\begin{eqnarray}\begin{array}{rcl} \dot{\rho }_{\alpha \beta }^{\nu ,(I)} & = & -\mathop{\sum }\limits_{\gamma ,\delta ,k,k^{\prime} }\int _{0}^{\infty }d\tau \\ {} & {} & g(\tau ){{S}_{\alpha \gamma ,k}}{{e}^{i{{\Delta }_{\alpha \gamma ,k}}t}}{{S}_{\gamma \delta ,k^{\prime} }}{{e}^{i{{\Delta }_{\gamma \delta ,k^{\prime} }}(t-\tau )}}\rho _{\delta \beta }^{\nu ,(I)} \\ {} & {} & -g(-\nu +\tau ){{S}_{\alpha \gamma ,k}}{{e}^{i{{\Delta }_{\alpha \gamma ,k}}(t-\tau )}}\rho _{\gamma \delta }^{\nu ,(I)}{{S}_{\delta \beta ,k^{\prime} }}{{e}^{i{{\Delta }_{\delta \beta ,k^{\prime} }}t}} \\ \end{array}\end{eqnarray} \tag{ B.4 }$$

$$\begin{eqnarray*}\begin{array}{ll} {} & {} & -g(-\nu -\tau ){{S}_{\alpha \gamma ,k}}{{e}^{i{{\Delta }_{\alpha \gamma ,k}}t}}\rho _{\gamma \delta }^{\nu ,(I)}{{S}_{\delta \beta ,k^{\prime} }}{{e}^{i{{\Delta }_{\delta \beta ,k^{\prime} }}(t-\tau )}} \\ {} & {} & +g(-\tau )\rho _{\alpha \gamma }^{\nu ,(I)}{{S}_{\gamma \delta ,k}}{{e}^{i{{\Delta }_{\gamma \delta ,k}}(t-\tau )}}{{S}_{\delta \beta ,k^{\prime} }}{{e}^{i{{\Delta }_{\delta \beta ,k^{\prime} }}t}}. \\ \end{array}\end{eqnarray*}$$

Equation (B.4) can be rewritten as:

$$\begin{eqnarray}\begin{array}{rcl} \dot{\rho }_{\alpha \beta }^{\nu ,(I)} & = & -\mathop{\sum }\limits_{\gamma ,\delta ,k,k^{\prime} }{{I}_{+}}(0,{{\Delta }_{\gamma \delta ,k^{\prime} }}){{S}_{\alpha \gamma ,k}}{{S}_{\gamma \delta ,k^{\prime} }}{{e}^{i({{\Delta }_{\alpha \gamma ,k}}+{{\Delta }_{\gamma \delta ,k^{\prime} }})t}}\rho _{\delta \beta }^{\nu ,(I)} \\ {} & {} & -{{I}_{+}}(-\nu ,{{\Delta }_{\alpha \gamma ,k}}){{S}_{\alpha \gamma ,k}}{{S}_{\delta \beta ,k^{\prime} }}{{e}^{i({{\Delta }_{\alpha \gamma ,k}}+{{\Delta }_{\delta \beta ,k^{\prime} }})t}}\rho _{\gamma \delta }^{\nu ,(I)} \\ {} & {} & -{{I}_{-}}(-\nu ,{{\Delta }_{\delta \beta ,k^{\prime} }}){{S}_{\alpha \gamma ,k}}{{S}_{\delta \beta ,k^{\prime} }}{{e}^{i({{\Delta }_{\alpha \gamma ,k}}+{{\Delta }_{\delta \beta ,k^{\prime} }})t}}\rho _{\gamma \delta }^{\nu ,(I)} \\ {} & {} & +{{I}_{-}}(0,{{\Delta }_{\gamma \delta ,k}}){{S}_{\gamma \delta ,k}}{{S}_{\delta \beta ,k^{\prime} }}{{e}^{i({{\Delta }_{\gamma \delta ,k}}+{{\Delta }_{\delta \beta ,k^{\prime} }})t}}\rho _{\alpha \gamma }^{\nu ,(I)}, \\ \end{array}\end{eqnarray} \tag{ B.5 }$$

where we have introduced the integrals of the correlation function

$$\begin{eqnarray}&&{{I}_{\pm }}(\nu ,E)=\int _{0}^{\infty }d\tau g(\nu \pm \tau ){{e}^{-iE\tau }}.\end{eqnarray} \tag{ B.6 }$$

The integrals (B.6) can be explicitly evaluated:

$$\begin{eqnarray}\begin{array}{rcl} {{I}_{\pm }}(\nu ,E) & = & \int _{0}^{\infty }d\tau {{e}^{-iE\tau }}\int _{0}^{\infty }d\omega \\ {} & {} & J(\omega )\left[ {{e}^{i\omega \nu }}{{e}^{\pm i\omega \tau }}\;{{n}_{B}}(\omega )+{{e}^{-i\omega \nu }}{{e}^{\mp i\omega \tau }}\left( {{n}_{B}}(\omega )+1 \right) \right]. \\ \end{array}\end{eqnarray} \tag{ B.7 }$$

Using the identity $\int _{0}^{\infty }{{e}^{i\omega t}}dt=\pi \delta (\omega )+\mathcal{P}(\frac{i}{\omega })$ and neglecting principal-value integrals, we obtain

$$\begin{eqnarray}&&{{I}_{\pm }}(\nu ,E)=\pi {{e}^{\pm }}^{iE\nu }s(\pm E),\end{eqnarray} \tag{ B.8 }$$

where

$$\begin{eqnarray}&&s(E)=\theta (E)J(E){{n}_{B}}(E)+\theta (-E)J(-E)\left[ {{n}_{B}}(-E)+1 \right]\end{eqnarray} \tag{ B.9 }$$

and $\theta (E)$ is the Heaviside step function.

Our results so far apply to any periodically driven system. In the following, we will take the fast-drive limit.

Appendix. GQME in the full secular approximation

We assume Ω to be faster than the decoherence rates that appear in the GQME. This allows us to neglect fast-oscillating terms of the form ${{e}^{ik\Omega t}}$, $k\ne 0$. A further class of terms, oscillating as ${{e}^{i({{\epsilon }_{\alpha }}-{{\epsilon }_{\beta }})t}}$, $\alpha \ne \beta$, can be safely neglected in most [31, 51], but not all [52], instances. Altogether, these two approximations amount to a full secular approximation, which results in a decoupling of the equations for the diagonal and off-diagonal elements in (B.5). We note in passing that in this approximation, the principal-value integrals that we neglected in deriving (B.8) would not contribute to $G(\nu ,t)$.

In order to determine $G(\nu ,t),$ it is enough to solve (B.5) for the diagonal elements $\rho _{ii}^{\nu }$. We define

$$\begin{eqnarray}\begin{array}{rcl} {{a}_{\alpha \beta ,k}} & = & 2\pi s({{\Delta }_{\alpha \beta ,k}})|{{S}_{\alpha \beta ,k}}{{|}^{2}}, \\ {{A}_{\alpha \beta }} & = & \mathop{\sum }\limits_{k}{{a}_{\alpha \beta ,k}}, \\ A_{\alpha \beta }^{\nu } & = & \mathop{\sum }\limits_{k}{{e}^{-i\nu {{\Delta }_{\alpha \beta ,k}}}}{{a}_{\alpha \beta ,k}}. \\ \end{array}\end{eqnarray} \tag{ B.10 }$$

At this point, we focus on a two-level system, which simplifies the analytical treatment. The GQME (B.5) can be written as:

$$\begin{eqnarray}\begin{array}{rcl} \dot{\rho }_{11}^{\nu } & = & [A_{11}^{\nu }-{{A}_{11}}-{{A}_{21}}]\rho _{11}^{\nu }+A_{12}^{\nu }\rho _{22}^{\nu }, \\ \dot{\rho }_{22}^{\nu } & = & A_{21}^{\nu }\rho _{11}^{\nu }+[A_{22}^{\nu }-{{A}_{22}}-{{A}_{12}}]\rho _{22}^{\nu }. \\ \end{array}\end{eqnarray} \tag{ B.11 }$$

We can also write (B.11) in matrix form by introducing a vector $\overrightarrow{{{\rho }^{\nu }}}=\{\rho _{11}^{\nu },\rho _{22}^{\nu }\}$ and a matrix $\mathcal{A}$ so that

$$\begin{eqnarray}&&\mathop{\overrightarrow{{{\rho }^{\nu }}}}\limits^{.}=\mathcal{A}\cdot \overrightarrow{{{\rho }^{\nu }}}.\end{eqnarray} \tag{ B.12 }$$

In the limit $\nu \to 0$, one recovers the standard master equation

$$\begin{eqnarray}\begin{array}{rcl} {{{\dot{\rho }}}_{11}} & = & -{{A}_{21}}{{\rho }_{11}}+{{A}_{12}}{{\rho }_{22}}, \\ {{{\dot{\rho }}}_{22}} & = & {{A}_{21}}{{\rho }_{11}}-{{A}_{12}}{{\rho }_{22}}, \\ \end{array}\end{eqnarray} \tag{ B.13 }$$

whose stationary solutions or dynamical steady state (DSS) are the well-known [31]

$$\begin{eqnarray}\begin{array}{rcl} \rho _{11}^{{\rm st}} & = & \frac{{{A}_{12}}}{{{A}_{12}}+{{A}_{21}}}, \\ \rho _{22}^{{\rm st}} & = & \frac{{{A}_{21}}}{{{A}_{12}}+{{A}_{21}}}. \\ \end{array}\end{eqnarray} \tag{ B.14 }$$

Notice that for a two-level system any noise operator S is traceless, so that ${{S}_{22}}=-{{S}_{11}}$. As a result, the coefficients appearing in (B.10) and (B.11) satisfy the relation ${{a}_{22,k}}={{a}_{11,k}}$, which in turn implies ${{A}_{22}}={{A}_{11}}$ and $A_{22}^{\nu }=A_{11}^{\nu }$.

Let us calculate the time variation of the average heat, which we refer to as the mean heat power. It is given by

$$\begin{eqnarray}&&\langle \dot{Q}\rangle =\frac{d}{dt}\langle Q\rangle =-i{{\left. \frac{d}{dt}\frac{d}{d\nu }G(\nu ,t) \right|}_{\nu =0}}=-i{{\left. \frac{d}{d\nu }\dot{G}(\nu ,t) \right|}_{\nu =0}}.\end{eqnarray} \tag{ C.1 }$$

From the definition of the CF we have that $G(\nu ,t)={\rm T}{{{\rm r}}_{S}}[{{\rho }^{\nu }}]=\rho _{11}^{\nu }+\rho _{22}^{\nu }$. The evolution of the generating function is given by $\dot{G}(\nu ,t)=\dot{\rho }_{11}^{\nu }+\dot{\rho }_{22}^{\nu }$. Using (B.11) we obtain

$$\begin{eqnarray}&&\langle \dot{Q}\rangle =-\mathop{\sum }\limits_{\alpha ,\beta ,k}{{\Delta }_{\alpha \beta ,k}}\ {{a}_{\alpha \beta ,k}}\ {{\rho }_{\beta \beta }}.\end{eqnarray} \tag{ C.2 }$$

This equation directly relates the dissipated heat power to the dynamics of the system. Given a solution ρ of the standard master equation, which is in general easier to obtain than the one for the GQME, we can immediately predict the dissipated power using (C.2). The most interesting application is in determining if and how much heat is dissipated when the system reaches the DDS (B.14). As discussed in the main text and in the following examples, this allows us to classify the driven systems in dissipative or nondissipative at steady state.

We define

$$\begin{eqnarray}\begin{array}{rcl} \Sigma & = & {{A}_{12}}+{{A}_{21}}, \\ {{\Upsilon }_{\alpha \beta }}(\nu ) & = & A_{\alpha \beta }^{\nu }-{{A}_{\alpha \beta }}, \\ \Upsilon _{\alpha \beta }^{(2)}(\nu ) & = & A_{\alpha \beta }^{\nu }A_{\beta \alpha }^{\nu }-{{A}_{\alpha \beta }}{{A}_{\beta \alpha }}. \\ \end{array}\end{eqnarray} \tag{ D.1 }$$

The eigenvalues of $\mathcal{A}$ are

$$\begin{eqnarray}&&{{\xi }_{\pm }}(\nu )=\frac{1}{2}\left( r\pm h \right)\end{eqnarray} \tag{ D.2 }$$

with

$$\begin{eqnarray}&&h=\sqrt{{{\Sigma }^{2}}+4\Upsilon _{12}^{(2)}(\nu )}\ ,\end{eqnarray} \tag{ D.3 }$$

$$\begin{eqnarray}&&r=-\Sigma +{{\Upsilon }_{11}}(\nu )+{{\Upsilon }_{22}}(\nu )\ .\end{eqnarray} \tag{ D.4 }$$

The corresponding eigenvectors are

$$\begin{eqnarray}&&v_{\pm }^{\nu }=\frac{1}{{{A}_{12}}-{{A}_{21}}+2A_{21}^{\nu }\pm h}\left( \frac{{{A}_{12}}-{{A}_{21}}\pm h}{2A_{21}^{\nu }} \right).\end{eqnarray} \tag{ D.5 }$$

Notice that in this representation the trace of a density operator is obtained by summing the two vector components. In the above equation, we have chosen the normalization constant for $v_{\pm }^{\nu }$ so that ${\rm Tr}[v_{\pm }^{\nu }]=1$.

The initial density matrix can be characterized by the quantity $z={{\rho }_{22}}(0)-{{\rho }_{11}}(0)$. We can write it as

$$\begin{eqnarray}&&\rho (0)=\left( \frac{{{\rho }_{11}}(0)}{{{\rho }_{22}}(0)} \right)=\left( \frac{\frac{1-z}{2}}{\frac{1+z}{2}} \right)=c_{+}^{\nu }v_{+}^{\nu }+c_{-}^{\nu }v_{-}^{\nu }\end{eqnarray} \tag{ D.6 }$$

with

$$\begin{eqnarray}&&c_{-}^{\nu }=\frac{1}{2\;h}\left[ (h-A_{12}^{\nu }-A_{21}^{\nu })+({{\Upsilon }_{21}}(\nu )-{{\Upsilon }_{12}}(\nu ))z \right],\end{eqnarray} \tag{ D.7 }$$

$$\begin{eqnarray}&&c_{+}^{\nu }=\frac{1}{2\;h}\left[ (h+A_{12}^{\nu }+A_{21}^{\nu })+({{\Upsilon }_{12}}(\nu )-{{\Upsilon }_{21}}(\nu ))z \right].\end{eqnarray} \tag{ D.8 }$$

Putting things together, we finally write the CF as

$$\begin{eqnarray}&&G(\nu ,t)=c_{-}^{\nu }{{e}^{{{\xi }_{-}}(\nu )t}}+c_{+}^{\nu }{{e}^{{{\xi }_{+}}(\nu )t}}\ .\end{eqnarray} \tag{ D.9 }$$

Equation (D.9) is the exact solution of the GQME (B.11). It can be used to calculate all the moments and cumulants of the PDF. In general, in order to obtain the full PDF, we must numerically invert the Fourier transform in (3). In the following, however, we show that some analytical insight can be gained in the long-time limit.

Appendix. Analytical expressions for the first moments

The moments $\langle {{Q}^{n}}(t)\rangle ={{\left( -i \right)}^{n}}{{\left. \frac{{{d}^{n}}G(\nu ,t)}{d{{\nu }^{n}}} \right|}_{\nu =0}}$ can be calculated directly from (D.9). From a physical point of view, the central moments $\langle {{(\delta Q)}^{n}}\rangle$, with $\delta Q=Q-\langle Q\rangle$, are the most meaningful quantities. For $n=1,2,3$, these coincide with the cumulants ${{\langle \langle Q\rangle \rangle }_{n}}$ [55]. In the long-time limit, the cumulants can be calculated directly from the dominant eigenvalues ${{\xi }_{+}}(\nu )$, as

$$\begin{eqnarray}&&{{\langle \langle Q\rangle \rangle }_{n}}={{\left( -i \right)}^{n}}t\frac{{{\partial }^{n}}}{\partial {{\nu }^{n}}}{{\xi }_{+}}(\nu )\left| _{\nu =0}={{\left( -i \right)}^{n}}\xi _{+}^{(n)}(0)t\ \right..\end{eqnarray} \tag{ D.10 }$$

This is possible because at $\nu =0$, ${{\xi }_{+}}(0)=0$ and ${{\xi }_{-}}(0)=-\Sigma \lt 0$. As a result, all the terms involving derivatives of ξ₋ undergo exponential decay.

A PDF is often characterized by considering its first three central moments, namely, the mean value $\bar{Q}=\langle Q\rangle$, the variance ${{\sigma }^{2}}=\langle {{(\delta Q)}^{2}}\rangle$, and the skewness $\kappa =\langle {{(\delta Q)}^{3}}\rangle$. In the long time limit, all these quantities increase linearly in time and are proportional to the derivatives of the dominant eigenvalue ξ₊ at $\nu =0$:

$$\begin{eqnarray}&&\bar{Q}=-i\xi _{+}^{(1)}(0)t,\end{eqnarray} \tag{ D.11 }$$

$$\begin{eqnarray}&&{{\sigma }^{2}}=-\xi _{+}^{(2)}(0)t,\end{eqnarray} \tag{ D.12 }$$

$$\begin{eqnarray}&&\kappa =i\xi _{+}^{(3)}(0)t\ .\end{eqnarray} \tag{ D.13 }$$

Appendix. Eigenvalue analysis

Using (B.10) and (D.1), ${{\Upsilon }_{\alpha \alpha }}(\nu )$ can be rewritten as (for $\nu \in \mathbb{R}$)

$$\begin{eqnarray}&&{{\Upsilon }_{\alpha \alpha }}(\nu )=\mathop{\sum }\limits_{k}({{e}^{-i\nu k\Omega }}-1){{a}_{\alpha \alpha ,k}}\ ,\end{eqnarray} \tag{ D.14 }$$

so that $\Re {\rm e}[{{\Upsilon }_{\alpha \alpha }}(\nu )]\leqslant 0$. Likewise,

$$\begin{eqnarray}&&\Upsilon _{\alpha \beta }^{(2)}(\nu )=\mathop{\sum }\limits_{k,k^{\prime} }({{e}^{-i\nu (k+{{k}^{\prime }})\Omega }}-1){{a}_{\alpha \beta ,k}}{{a}_{\beta \alpha ,k^{\prime} }},\end{eqnarray} \tag{ D.15 }$$

so that also $\Re {\rm e}[\Upsilon _{\alpha \beta }^{(2)}(\nu )]\leqslant 0$.

Altogether, this implies that $\Re {\rm e}\ r\leqslant -({{A}_{12}}+{{A}_{21}})$ and $\Re {\rm e}\ h\leqslant ({{A}_{12}}+{{A}_{21}})$. This implies that $\Re {\rm e}\ {{\xi }_{-}}\lt 0$, $\Re {\rm e}\ {{\xi }_{+}}\leqslant 0$ and

$$\begin{eqnarray}&&\Re {\rm e}\ {{\xi }_{+}}=0\;{\rm for}\;\nu =n\tau ,\ n\in \mathcal{Z}\ .\end{eqnarray} \tag{ D.16 }$$

Notice that ${{\xi }_{\pm }}(\nu )$ are periodic with period τ, i.e. ${{\xi }_{\pm }}(\nu +n\tau )={{\xi }_{\pm }}(\nu )$ and that for the subdominant eigenvalues ${\rm max} [\Re {\rm e}\ {{\xi }_{-}}]=-\Sigma$ and the maximum is reached for $\nu =n\tau$. This means that the terms associated to the eigenvalue ξ₋ vanish exponentially as soon as $t\gg 1/{\rm min} [\Re {\rm e}\ {{\xi }_{-}}]=1/\Sigma$ and the behavior of $G(\nu ,t)$ is determined only by ξ₊. Furthermore, due to the fact that $\Re {\rm e}\ {{\xi }_{+}}\leqslant 0$ for $\nu \ne n\tau$, $G(\nu ,t)$ is different for zero only in a neighborhood of the resonances $\nu =n\tau$ and vanishes exponentially away from them.

D.1. Solution in the long-time limit

In the long-time limit, $G(\nu ,t)\approx c_{+}^{\nu }{{e}^{{{\xi }_{+}}(\nu )t}}$ vanishes everywhere except in a neighborhood of ${{\nu }_{n}}=n\tau$. We recall that ξ₊ is periodic in τ. In a neighborhood of ν_n, we perform a second-order (Gaussian) approximation and write

$$\begin{eqnarray}&&{{\xi }_{+}}(\nu )\approx i(\nu \Omega -2n\pi )a-{{\left( \nu \Omega -2n\pi \right)}^{2}}b,\end{eqnarray} \tag{ D.17 }$$

where $a=-i\xi _{+}^{(1)}(0)$ and $b=-1/2\xi _{+}^{(2)}(0)$ with $a\leqslant 0$ and $b\gt 0$. Then we approximate the CF as

$$\begin{eqnarray}&&G(\nu ,t)=\mathop{\sum }\limits_{n=-\infty }^{\infty }{{G}_{n}}{{e}^{i(\nu \Omega -2n\pi )at-{{\left( \nu \Omega -2n\pi \right)}^{2}}bt}},\end{eqnarray} \tag{ D.18 }$$

where ${{G}_{n}}=c_{+}^{\nu =n\tau }$. Explicitly:

$$\begin{eqnarray}\begin{array}{rcl} {{G}_{n}} & = & 1+\frac{{{\Upsilon }_{12}}(n\tau )+{{\Upsilon }_{21}}(n\tau )+[{{\Upsilon }_{12}}(n\tau )-{{\Upsilon }_{21}}(n\tau )]z}{2({{A}_{12}}+{{A}_{21}})} \\ {} & = & 1+\frac{{{\rho }_{11}}(0){{\Upsilon }_{21}}(n\tau )+{{\rho }_{22}}(0){{\Upsilon }_{12}}(n\tau )}{{{A}_{12}}+{{A}_{21}}}. \\ \end{array}\end{eqnarray} \tag{ D.19 }$$

using (B.10) and (D.1), we can write

$$\begin{eqnarray}&&{{\Upsilon }_{\alpha \beta }}(n\tau )=\left[ {{e}^{-in({{\epsilon }_{\alpha }}-{{\epsilon }_{\beta }})\tau }}-1 \right]{{A}_{\alpha \beta }}.\end{eqnarray} \tag{ D.20 }$$

The Fourier transform in (3) can be inverted analytically starting from (D.18). The result is

$$\begin{eqnarray}&&P(Q,t)=w(Q,t)\mathop{\sum }\limits_{n=-\infty }^{\infty }{{G}_{n}}{{e}^{-inQ\tau }}\ ,\end{eqnarray} \tag{ D.21 }$$

where we have introduced the weight function

$$\begin{eqnarray}&&w(Q,t)=\frac{1}{\sqrt{2\pi }}\frac{1}{\sqrt{2bt}\Omega }{\rm exp} \left[ -\frac{{{\left( Q-at\Omega \right)}^{2}}}{4bt{{\Omega }^{2}}} \right].\end{eqnarray} \tag{ D.22 }$$

Making use of (D.19) and (D.20), we obtain

$$\begin{eqnarray}&&P(Q,t)=w(Q,t)\mathop{\sum }\limits_{n=-\infty }^{\infty }\left[ {{p}_{\downarrow }}{{e}^{-in(Q+{{\epsilon }_{1}}-{{\epsilon }_{2}})\tau }}+{{p}_{\uparrow }}{{e}^{-in(Q+{{\epsilon }_{2}}-{{\epsilon }_{1}})\tau }}+(1-{{p}_{\uparrow }}-{{p}_{\downarrow }}){{e}^{-inQ\tau }} \right],\end{eqnarray} \tag{ D.23 }$$

where

$$\begin{eqnarray}\begin{array}{rcl} {{p}_{\uparrow }} & = & \frac{{{A}_{21}}{{\rho }_{11}}(0)}{{{A}_{12}}+{{A}_{21}}}, \\ {{p}_{\downarrow }} & = & \frac{{{A}_{12}}{{\rho }_{22}}(0)}{{{A}_{12}}+{{A}_{21}}}. \\ \end{array}\end{eqnarray} \tag{ D.24 }$$

Now we recall the identity (Dirac comb)

$$\begin{eqnarray}&&\mathop{\sum }\limits_{n}{{e}^{-inQ\tau }}=\Omega \mathop{\sum }\limits_{n}\delta (Q+n\Omega ),\end{eqnarray} \tag{ D.25 }$$

with $\tau =2\pi /\Omega$. Substituting (D.25) into (D.23), we finally get

$$\begin{eqnarray}\begin{array}{rcl} P(Q,t) & = & w(Q,t)\Omega \mathop{\sum }\limits_{n=-\infty }^{\infty }\left[ {{p}_{\downarrow }}\delta (Q+{{\epsilon }_{1}}-{{\epsilon }_{2}}+n\Omega )+{{p}_{\uparrow }}\delta (Q+{{\epsilon }_{2}}-{{\epsilon }_{1}}+n\Omega ) \right. \\ {} & {} & \left. +(1-{{p}_{\uparrow }}-{{p}_{\downarrow }})\delta (Q+n\Omega ) \right]. \\ \end{array}\end{eqnarray} \tag{ D.26 }$$

The analytical result (D.26) stems from the Gaussian expansion (D.17) and is valid in the long-time limit. This limit is defined by the condition $t\gg 1/b$, for which the CF is localized at the nodes $\nu =n\tau$. This is also the limit where many photons have been exchanged, as the standard deviation of the Gaussian function $\Omega \sqrt{2bt}$ is much larger than the photon energy Ω. In this limit, the delta functions in (D.26) tend to a continuum and, as a result, we can write $P(Q,t)\approx w(Q,t)$. The coefficient a and b are then directly linked to the average exchanged heat $\bar{Q}=\Omega at$ and variance ${{\sigma }^{2}}=2{{\Omega }^{2}}bt$, as expected from the analysis of the central moments.

Let us consider a two-level system driven by a transverse monochromatic drive. This is also referred to as the Rabi model, described by the Hamiltonian

$$\begin{eqnarray}&&H(t)=\frac{\omega }{2}(1+{{\sigma }_{z}})+g\left[ {{e}^{i(\Omega t-\varphi )}}{{\sigma }_{+}}+{{e}^{-i(\Omega t-\varphi )}}{{\sigma }_{-}} \right],\end{eqnarray} \tag{ E.1 }$$

where σ_i are the usual Pauli operators, defined with respect to the fixed basis $\{|0\rangle ,|1\rangle \}$. The Floquet modes read

$$\begin{eqnarray}\begin{array}{rcl} |{{\phi }_{1}}\rangle & = & {\rm cos} \theta |0\rangle -{{e}^{-i(\Omega t+\varphi )}}{\rm sin} \theta |1\rangle , \\ |{{\phi }_{2}}\rangle & = & {\rm sin} \theta |0\rangle {{e}^{-i\varphi }}+{{e}^{-i\Omega t}}{\rm cos} \theta |1\rangle \\ \end{array}\end{eqnarray} \tag{ E.2 }$$

and the corresponding quasi-energies are

$$\begin{eqnarray}&&{{\epsilon }_{1,2}}=\frac{1}{2}(\Delta \pm {{\Omega }_{R}}),\end{eqnarray} \tag{ E.3 }$$

where we have introduced the quantities $\Delta =\omega -\Omega$, ${\rm tan} 2\theta =2\;g/\Delta$ and ${{\Omega }_{R}}=\sqrt{{{\Delta }^{2}}+4\;{{g}^{2}}}$.

In order to write the GQME (B.11), we need to calculate the matrix elements of the noise operator S in (B.2). We consider an Ohmic spectral density $J(\omega )=\eta \omega$ (we assume that Ω is much smaller than the cutoff frequency of the environment). Notice from (B.9) that for $E\gt 0$ we have $s(E)/s(-E)={{e}^{-\beta E}}$.

E.1. Transverse coupling ($S={{\sigma }_{x}}$)

For $S={{\sigma }_{x}}$ only the terms with $k=\pm 1$ in (B.2) contribute and we have

$$\begin{eqnarray}\begin{array}{rcl} {{S}_{11,-1}} & = & {{S}_{11,1}}=-{{S}_{22,-1}}=-{{S}_{22,1}}=-\frac{1}{2}{\rm sin} 2\theta , \\ {{S}_{12,-1}} & = & {{S}_{21,1}}={{{\rm cos} }^{2}}\theta , \\ {{S}_{21,-1}} & = & {{S}_{12,1}}=-{\rm sin} 2\theta . \\ \end{array}\end{eqnarray} \tag{ E.4 }$$

From these expressions, the calculation follows the one leading to (D.26) with the definitions B.10.

E.2. Longitudinal coupling ($S={{\sigma }_{z}}$)

For $S={{\sigma }_{z}}$, only the terms with k = 0 contribute and we have

$$\begin{eqnarray}\begin{array}{rcl} {{S}_{11,0}} & = & -{{S}_{22,0}}={\rm cos} 2\theta , \\ {{S}_{12,0}} & = & {{S}_{21,0}}={\rm sin} 2\theta . \\ \end{array}\end{eqnarray} \tag{ E.5 }$$

In the case, the GQME reads

$$\begin{eqnarray}\begin{array}{rcl} \dot{\rho }_{11}^{\nu } & = & 2\pi S_{12,0}^{2}\left[ -\rho _{11}^{\nu }s({{\Omega }_{R}})+\rho _{22}^{\nu }{{e}^{i{{\Omega }_{R}}\nu }}s(-{{\Omega }_{R}}) \right], \\ \dot{\rho }_{22}^{\nu } & = & 2\pi S_{12,0}^{2}\left[ -\rho _{22}^{\nu }s(-{{\Omega }_{R}})+\rho _{11}^{\nu }{{e}^{-i{{\Omega }_{R}}\nu }}s({{\Omega }_{R}}) \right]. \\ \end{array}\end{eqnarray} \tag{ E.6 }$$

The GQME (E.6) can be solved analytically. By using the notation ${{\Gamma }_{\pm }}=2\pi S_{12,0}^{2}\;s(\pm {{\Omega }_{R}})+s({{\Omega }_{R}})$ and $\Gamma ={{\Gamma }_{+}}+{{\Gamma }_{-}}$, the CF reads

$$\begin{eqnarray}\begin{array}{rcl} G(\nu ,t) & = & \frac{{{\Gamma }_{+}}}{\Gamma }{{\rho }_{11}}(0)\left( 1-{{e}^{-\Gamma t}} \right){{e}^{-i{{\Omega }_{R}}\nu }} \\ {} & {} & +\frac{{{\Gamma }_{-}}}{\Gamma }{{\rho }_{22}}(0)\left( 1-{{e}^{-\Gamma t}} \right){{e}^{i{{\Omega }_{R}}\nu }} \\ {} & {} & +\frac{{{\Gamma }_{+}}}{\Gamma }\left[ {{\rho }_{22}}(0)+{{\rho }_{11}}(0){{e}^{-\Gamma t}} \right] \\ {} & {} & +\frac{{{\Gamma }_{-}}}{\Gamma }\left[ {{\rho }_{11}}(0)+{{\rho }_{22}}(0){{e}^{-\Gamma t}} \right]. \\ \end{array}\end{eqnarray} \tag{ E.7 }$$

From this CF, we see that the moments have a particular symmetry. As ${{\partial }^{k}}/\partial {{\nu }^{k}}\left( {{e}^{\pm i{{\Omega }_{R}}\nu }} \right)={{\left( \pm i{{\Omega }_{R}} \right)}^{k}}{{e}^{\pm i{{\Omega }_{R}}\nu }}$, the kth moment reads

$$\begin{eqnarray}&&\langle {{Q}^{k}}\rangle \equiv \frac{1}{{{i}^{k}}}\frac{{{\partial }^{k}}G(\nu ,t)}{\partial {{\nu }^{k}}}{{|}_{\nu =0}}=\Omega _{R}^{k}\left[ {{p}_{\downarrow }}(t)+{{\left( -1 \right)}^{k}}{{p}_{\uparrow }}(t) \right]\end{eqnarray} \tag{ E.8 }$$

where

$$\begin{eqnarray}\begin{array}{rcl} {{p}_{\downarrow }}(t) & = & \frac{{{\Gamma }_{-}}}{\Gamma }{{\rho }_{22}}(0)\left( 1-{{e}^{-\Gamma t}} \right), \\ {{p}_{\uparrow }}(t) & = & \frac{{{\Gamma }_{+}}}{\Gamma }{{\rho }_{11}}(0)\left( 1-{{e}^{-\Gamma t}} \right). \\ \end{array}\end{eqnarray} \tag{ E.9 }$$

The PDF is immediately obtained from the inverse Fourier transform of the CF. It reads

$$\begin{eqnarray}&&P(Q)={{p}_{\downarrow }}(t)\delta (Q-{{\Omega }_{R}})+{{p}_{\uparrow }}(t)\delta (Q+{{\Omega }_{R}})+{{p}_{0}}(t)\delta (Q),\end{eqnarray} \tag{ E.10 }$$

where ${{p}_{\downarrow }}(t)$ is the probability for the environment to absorb a photon (the system emits it), ${{p}_{\uparrow }}(t)$ to emit a photon (the system absorbs it), and ${{p}_{0}}(t)=1-{{p}_{\downarrow }}(t)-{{p}_{\uparrow }}(t)$ corresponds to no emissions or absorptions.

Appendix. Dependence on the initial state

The moments $\langle {{Q}^{k}}\rangle$ in (E.8) depend on the initial state. If the system is initially thermalized in the Floquet basis, i.e. ${{\rho }_{22}}(0)={{e}^{-\beta {{\Omega }_{R}}}}{{\rho }_{11}}(0)$, by using the relation $s({{\Omega }_{R}})=s(-{{\Omega }_{R}}){{e}^{-\beta {{\Omega }_{R}}}}$, we verify that

$$\begin{eqnarray}\langle {{Q}^{k}}\rangle =\left\{ \begin{array}{ccccccccccccccc} 0 & \;{\rm for}\;k\;{\rm odd}, \\ \Omega _{R}^{k}\frac{2\left( 1-{{e}^{-\Gamma t}} \right){{e}^{\beta \epsilon }}}{{{\left( 1+{{e}^{\beta \epsilon }} \right)}^{2}}} & \;{\rm for}\;k\;{\rm even}. \\ \end{array} \right.\end{eqnarray} \tag{ E.11 }$$

That means that $\langle Q\rangle =0$ and, in general, $\langle {{Q}^{2k+1}}\rangle =0$. From this result [or directly from (E.8)], it is immediate to see that the average heat power $\langle \dot{Q}\rangle =0$ at DSS, as discussed in the main text.

If the system is initially in the state $|{{\psi }_{0}}\rangle ={\rm cos} \delta |0\rangle -{{e}^{i\gamma }}{\rm sin} \delta |1\rangle$, the corresponding populations in the Floquet basis are

$$\begin{eqnarray}&&{{\rho }_{11}}(0)=\frac{1}{2}[1+{\rm cos} 2\delta {\rm cos} 2\theta -{\rm sin} 2\delta {\rm sin} 2\theta {\rm cos} (\gamma +\varphi )]\end{eqnarray} \tag{ E.12 }$$

and ${{\rho }_{22}}(0)=1-{{\rho }_{11}}(0)$. For a resonant drive, i.e. $\theta =\pi /4$, $\delta =\pi /4$ and $\gamma =0$, we have ${{\rho }_{11}}(0)=1/2(1-{\rm cos} \varphi )$. In the long time limit, we have

$$\begin{eqnarray}&&\frac{\langle Q\rangle }{{{\Omega }_{R}}}=\frac{s(-{{\Omega }_{R}})}{s({{\Omega }_{R}})+s(-{{\Omega }_{R}})}-\frac{1}{2}(1-{\rm cos} \varphi )\ .\end{eqnarray} \tag{ E.13 }$$

In the low temperature regime, i.e., $s({{\Omega }_{R}})/s(-{{\Omega }_{R}})\ll 1$, we have $\langle Q\rangle /{{\Omega }_{R}}=1/2(1+{\rm cos} \varphi )$, which vanishes for $\varphi =\pi$. In the high temperature regime, i.e., $s({{\Omega }_{R}})/s(-{{\Omega }_{R}})\approx 1$, we have $\langle Q\rangle /{{\Omega }_{R}}=1/2{\rm cos} \varphi$. which has no definite sign and attains its minimum $\langle Q\rangle /{{\Omega }_{R}}=-1/2$ at $\varphi =\pi$. A negative $\langle Q\rangle /{{\Omega }_{R}}$ implies a net heat flow from the environment into the system.

For completeness, we also present a derivation of the GQME for an undriven two-level system, starting from (A.8). We take $H=\frac{\omega }{2}{{\sigma }_{z}}$, $\omega \gt 0$, and $S={{\sigma }_{x}}$. This gives $S(-\tau )={{e}^{i\omega \tau }}{{\sigma }_{-}}+{{e}^{-i\omega \tau }}{{\sigma }_{+}}$. Substituting in (A.8), we obtain

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial }{\partial t}{{\rho }^{\nu }}(t) & = & -i\left[ H,{{\rho }^{\nu }} \right]-\int _{0}^{\infty }d\tau \\ {} & {} & g(\tau ){{\sigma }_{x}}[{{e}^{i\omega \tau }}{{\sigma }_{-}}+{{e}^{-i\omega \tau }}{{\sigma }_{+}}]{{\rho }^{\nu }} \\ {} & {} & -g(-\nu +\tau )[{{e}^{i\omega \tau }}{{\sigma }_{-}}+{{e}^{-i\omega \tau }}{{\sigma }_{+}}]{{\rho }^{\nu }}{{\sigma }_{x}} \\ {} & {} & -g(-\nu -\tau ){{\sigma }_{x}}{{\rho }^{\nu }}[{{e}^{i\omega \tau }}{{\sigma }_{-}}+{{e}^{-i\omega \tau }}{{\sigma }_{+}}] \\ {} & {} & +g(-\tau ){{\rho }^{\nu }}[{{e}^{i\omega \tau }}{{\sigma }_{-}}+{{e}^{-i\omega \tau }}{{\sigma }_{+}}]{{\sigma }_{x}}. \\ \end{array}\end{eqnarray} \tag{ F.1 }$$

using (B.8), we obtain

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial }{\partial t}{{\rho }^{\nu }}(t) & = & -i\left[ H,{{\rho }^{\nu }} \right] \\ {} & {} & -{{\sigma }_{x}}[{{I}_{+}}(0,-\omega ){{\sigma }_{-}}+{{I}_{+}}(0,\omega ){{\sigma }_{+}}]{{\rho }^{\nu }} \\ {} & {} & +[{{I}_{+}}(-\nu ,-\omega ){{\sigma }_{-}}+{{I}_{+}}(-\nu ,\omega ){{\sigma }_{+}}]{{\rho }^{\nu }}{{\sigma }_{x}} \\ {} & {} & +{{\sigma }_{x}}{{\rho }^{\nu }}[{{I}_{-}}(-\nu ,-\omega ){{\sigma }_{-}}+{{I}_{-}}(-\nu ,\omega ){{\sigma }_{+}}] \\ {} & {} & -{{\rho }^{\nu }}[{{I}_{-}}(0,-\omega ){{\sigma }_{-}}+{{I}_{-}}(0,\omega ){{\sigma }_{+}}]{{\sigma }_{x}}. \\ \end{array}\end{eqnarray} \tag{ F.2 }$$

We define the rates $d=\pi J(\omega )\left[ {{n}_{B}}(\omega )+1 \right]$ and $u=\pi J(\omega ){{n}_{B}}(\omega )$. Then

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial }{\partial t}{{\rho }^{\nu }}(t) & = & -i\left[ H,{{\rho }^{\nu }} \right] \\ {} & {} & -{{\sigma }_{x}}[d{{\sigma }_{-}}+u{{\sigma }_{+}}]{{\rho }^{\nu }} \\ {} & {} & +[{{e}^{-i\omega \nu }}d{{\sigma }_{-}}+{{e}^{i\omega \nu }}u{{\sigma }_{+}}]{{\rho }^{\nu }}{{\sigma }_{x}} \\ {} & {} & +{{\sigma }_{x}}{{\rho }^{\nu }}[{{e}^{i\omega \nu }}u{{\sigma }_{-}}+{{e}^{-i\omega \nu }}d{{\sigma }_{+}}] \\ {} & {} & -{{\rho }^{\nu }}[u{{\sigma }_{-}}+d{{\sigma }_{+}}]{{\sigma }_{x}}. \\ \end{array}\end{eqnarray} \tag{ F.3 }$$

The time-independent differential equations for the relevant matrix elements read

$$\begin{eqnarray}\begin{array}{rcl} \dot{\rho }_{11}^{\nu } & = & -2u\rho _{11}^{\nu }+2d{{e}^{i\omega \nu }}\rho _{22}^{\nu }, \\ \dot{\rho }_{22}^{\nu } & = & 2u{{e}^{-i\omega \nu }}\rho _{11}^{\nu }-2d\rho _{22}^{\nu }. \\ \end{array}\end{eqnarray} \tag{ F.4 }$$

The equations (F.4) have the same structure as the (E.6). We thus obtain the same CF and PDF, with the substitutions $\omega \to {{\Omega }_{R}}$, $d\to \pi S_{12,0}^{2}\;s(-{{\Omega }_{R}})$ and $u\to \pi S_{12,0}^{2}\;s({{\Omega }_{R}})$.