Heat-exchange statistics in driven open quantum systems

As the dimensions of physical systems approach the nanoscale, the laws of thermodynamics must be reconsidered due to the increased importance of fluctuations and quantum effects. While the statistical mechanics of small classical systems is relatively well understood, the quantum case still poses challenges. Here we set up a formalism that allows to calculate the full probability distribution of energy exchanges between a periodically driven quantum system and a thermalized heat reservoir. The formalism combines Floquet theory with a generalized master equation approach. For a driven two-level system and in the long-time limit, we obtain a universal expression for the distribution, providing clear physical insight into the exchanged energy quanta. We illustrate our approach in two analytically solvable cases and discuss the differences in the corresponding distributions. Our predictions could be directly tested in a variety of systems, including optical cavities and solid-state devices.


Introduction
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,12,13,14,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,17,18,19,20,21,22].However, while experimental demonstrations are still lacking, the very definition of basic concepts, such as, work, is still under debate [23,24,16,25].This difficulty is related to the measurement process: on the 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 an ancilla was considered in [20,21].It has also been proposed to measure the environment, rather than the system, to keep track of the energy exchanges between the two [26,27,28].
Here, we provide a general framework to understand 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,30,31,32].However, the statistics of the dissipated heat was 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 pertubative 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,35,36,37,38,39,40,41,42,43,44,45,46].

Theoretical framework
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 . 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 ⊗ 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 = 1 and denote with ρ T , ρ and ρ 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] where p[e 2 ; e 1 ] is the conditional probability that a measurement of H R gives e 2 at time t when it gave e 1 at time t 0 and p[e 1 ] is the probability to measure e 1 at time t 0 .
Introducing the projector P e j on the j-th state of the reservoir of energy e j and using the property P where U(t) is the evolution operator generated by H T (t).The CF is given by [48] p[e 2 ; e 1 ]p[e 1 ]e iν(e 2 −e 1 ) . ( We assume the initial total density matrix ρ T (0) = ρ(0) ⊗ ρ R (0) to be factorized into the system density matrix ρ(0) and the thermalized environment density matrix ρ R (0) = e −β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 ρ(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 e j P e j e ±iνe j = e ±iνH R , we write (3) as and finally as where and satisfies the equation of motion The equation for ρ ν T (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 ρ ν T is not unitary for ν = 0.For a periodic drive, H(t + τ ) = H(t), where τ is the drive period and Ω = 2π/τ the corresponding angular frequency.We assume that the unitary dynamics generated by H is known; this knowledge is captured in our formalism by the use of Floquet states and quasienergies [31].Floquet states are particular solutions of the Schrödinger equation of the form |Ψ α (t) = e −iǫαt |φ α (t) where the Floquet mode |φ α (t) satisfies |φ α (t + τ ) = |φ α (t) and ǫ α is its corresponding quasienergy.
The derivation of the GQME for the generalized density matrix of the system T ] is done in complete analogy to the standard master equation [32] and it is presented in Appendix A. We assume weak coupling between system and environment and fast autocorrelation time of the environment (Born-Markov approximation).We model the environment as a collection of harmonic oscillators with bosonic occupation n B (ω) and Ohmic spectral density J(ω) = ηω with η a dimensionless coupling coefficient.The so-obtained GQME can be simplified in different manners, depending on the time scales of the problem.In the fast-driving regime [31,50,51,52], one can safely neglect oscillating terms of the form e ikΩt , with k = 0, and, in most instances, also of the form e i(ǫα−ǫ β )t , with α = β.This amounts to a full secular approximation [52], which admits an analytical solution.

General results
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 ρ ν αα = Ψ α |ρ ν |Ψ α are decoupled from the coherences and satisfy a vector equation of the form ˙ ρ ν = A • ρ ν .For a two-level system, ρ ν = {ρ ν 11 , ρ ν 22 } and 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 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 got an analytical solution.The corresponding CF can be written as where ξ ± (ν) are the eigenvalues of A and c ν ± the projection of the initial density matrix ρ(0) onto the corresponding eigenvectors v ν ± , normalized so that Tr v ν ± = 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 Q n of the PDF [55] are all linear in time and solely determined by the dominant eigenvalue ξ + (0): where ξ + 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 ℜe [ξ − (ν)] is always negative, while ℜe [ξ + (ν)] is negative everywhere except at ν n = nτ , n ∈ Z, where it vanishes.From (10), we thus see that lim t→∞ G(ν, t) = 0 for ν = ν 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 ξ + (ν) 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 where with real coefficients a = −iξ + (0) ≥ 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 ∆ αβ,k .This discretization results from energy being available in photons of energy Ω, the drive frequency, and ǫ 2 − ǫ 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 ↑ − p ↓ and p ↑ (p ↓ ), 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(ν, t) is localized at the nodes ν = nτ .This assumption holds true only in the limit where many energy quanta have been exchanged.For very long times, t ≫ 1/b, one can neglect the discretization due to the Dirac combs in (12) and find that P (Q, t) ≈ w(Q, t).

Examples
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 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 quasienergy gap ǫ 2 −ǫ 1 equals the Rabi frequency Ω R = √ ∆ 2 + 4g 2 , where ∆ = Ω − ω is the drive detuning.We consider two types of coupling operators, S = σ x (transverse coupling) and S = σ z (longitudinal coupling).Both can be realized, e. g., in solid-state devices [56,57].
When S = σ 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 ↓ = p ↑ .The full PDF is shown in the main panel of figure 1 at times t 1 /τ = 80 (purple) and t 2 /τ = 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 Ω 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)].
In figure 1 (b,c), we show how the PDF is affected by changes in the drive frequency and the temperature of the environment.We take a long time (t/τ = 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 αβ 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 behaviour 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 ∆ = 0.
We next consider the noise operator S = σ 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 the symmetry of the problem, which forbids transitions with energy exchange nΩ and thereby causes the drive and the environment to be effectively decoupled [54].The resulting GQME is formally equivalent to that for an undriven system; its derivation can be found in [?].The corresponding PDF reads where Notice that the structure of this PDF is the same as the more general (12), with weight function w(E, t) = 1 for − Ω 2 < E < Ω 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 and |1 , then the drive phase ϕ also affects the heat distribution, depending on the projection of the initial state onto the xy plane.An example is shown in figure 2 (c) for an equal superposition of |0 and |1 .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.

Conclusions
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 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 Eqs. ( 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 quasienergy crossings [31,52].Several techniques developed in the context of the full-counting statistics may be readily adapted to the present setting, including weak-coupling Markovian approaches [61,62], non-Markovian perturbative expansion [63], recursive schemes [64,65] and waiting-time analysis [66,67,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.implying that the generalized correlation function g ν is a time-shifted ordinary correlation function.For a bosonic bath, the correlation function takes the form where J(ω) is the spectral density associated to the coupling operator B and n B (ω) the Bose distribution function.Performing the Born-Markov approximation and extending the integration times to infinity, we obtain Sometimes it is convenient to go back to the Schrödinger picture.In this case we have: where S(−τ ) stands for U † 0 (t − τ, t)SU 0 (t − τ, t).

Appendix B. GQME for a periodic drive
We consider the case of a periodic drive, so that H(t + τ ) = H(t) with τ = 2π/Ω.The dynamics of the driven system is captured by considering the Floquet states |Ψ α (t) and their corresponding quasienergies ǫ α .We also introduce the Floquet modes |φ α (t) , satisfying |φ α (t + τ ) = |φ α (t) and related to the Floquet states by The evolution operator of the system in the Floquet basis reads U 0 = α |Ψ α (t) Ψ α (0)|.We work in the interaction picture defined by U 0 (t).Then Equation (B.4) can be rewritten as: ρν,(I) where we have introduced the integrals of the correlation function The integrals (B.6) can be explicitly evaluated: Using the identity ∞ 0 e iωt dt = πδ(ω) + P( i ω ) and neglecting principal-value integrals, we obtain where and θ(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.

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Ωt , k = 0.A further class of terms, oscillating as e i(ǫα−ǫ β )t , α = β, 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(ν, t).
In order to determine G(ν, t), it is enough to solve (B.5) for the diagonal elements ρ ν ii .We define At this point we focus on a two-level system, which simplifies the analytical treatment.The GQME (B.5) can be written as: (B.11) We can also write (B.11) in matrix form by introducing a vector ρ ν = {ρ ν 11 , ρ ν 22 } and a matrix A so that ˙ ρ ν = A • ρ ν .(B.12) In the limit ν → 0, one recovers the standard master equation

Appendix C. Average heat power
Let us calculate the time variation of the average heat, which we refer to as the mean heat power.It is given by From the definition of the CF we have that The evolution of the generating function is given by Ġ(ν, t) = ρν 11 + ρν 22 .Using (B.11) we obtain 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 non-dissipative at steady-state.

Appendix D. General results in the long-time limit
We define The eigenvalues of A are 12 (ν) , (D.3) The corresponding eigenvectors are 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 ν ± so that Tr[v ν ± ] = 1.
The initial density matrix can be characterized by the quantity z = ρ 22 (0) − ρ 11 (0).We can write it as Putting things together, we finally write the CF as 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.

Analytical expressions for the first moments
The moments dν n ν=0 can be calculated directly from (D.9).From a physical point of view, the central moments (δQ) n , with δQ = Q − Q , are the most meaningful quantities.For n = 1, 2, 3, these coincide with the cumulants Q n [55].In the long-time limit, the cumulants can be calculated directly from the dominant eigenvalues ξ + (ν), as This is possible because at ν = 0, ξ + (0) = 0 and ξ − (0) = −Σ < 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 Q = Q , the variance σ 2 = (δQ) 2 , and the skewness κ = (δQ) 3 .In the long time limit, all these quantities increase linearly in time and are proportional to the derivatives of the dominant eigenvalue ξ + at ν = 0: + (0)t , (D.12) κ = iξ so that also ℜe[Υ

Solution in the long-time limit
In the long-time limit, G(ν, t) ≈ c ν + e ξ + (ν)t vanishes everywhere except in a neighborhood of ν n = nτ .We recall that ξ + is periodic in τ .In a neighborhood of ν n , we perform a second-order (Gaussian) approximation and write where a = −iξ + (0) with a ≤ 0 and b > 0. Then we approximate the CF as where G n = c ν=nτ + .Explicitly: The Fourier transform in (3) can be inverted analytically starting from (D.18).The result is where we have introduced the weight function Making use of (D. 19) and (D.20), we obtain 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 ≫ 1/b, for which the CF is localized at the nodes ν = nτ .This is also the limit where many photons have been exchanged, as the standard deviation of the Gaussian function 2Ω √ bt 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) ≈ w(Q, t).The coefficient a and b are then directly linked to the average exchanged heat Q = Ωat and variance σ 2 = 4Ω 2 bt, as expected from the analysis of the central moments.

Appendix E. Two-level system with a monochromatic drive
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 where σ i are the usual Pauli operators, defined with respect to the fixed basis {|0 , |1 }.The Floquet modes read |φ 1 = cos θ|0 − e −i(Ωt+ϕ) sin θ|1 , and the corresponding quasienergies are where we have introduced the quantities ∆ 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(ω) = ηω (we assume that Ω is much smaller than the cutoff frequency of the environment).Notice from (B.9) that for E > 0 we have s(E)/s(−E) = e −βE .

Transverse coupling (S = σ x )
For S = σ x only the terms with k = ±1 in (B.2) contribute and we have (E.5) In the case, the GQME reads The GQME (E.6) can be solved analytically.By using the notation Γ ± = 2πS 2  12,0 s(±Ω R ) + s(Ω R ) and Γ = Γ + + Γ − , the CF reads (E.7) From this CF we see that the moments have a particular symmetry.As where (E.9) The PDF is immediately obtained from the inverse Fourier transform of the CF.It reads where p ↓ (t) is the probability for the environment to absorb a photon (the system emits it), p ↑ (t) to emit a photon (the system absorbs it), and p 0 (t) = 1 − p ↓ (t) − p ↑ (t) corresponds to no emissions or absorptions.

Figure 2 . 1 √ 2 (
Figure 2. Dissipated energy for longitudinal coupling.(a) PDF at different times.(b, c) First three central moments of the dissipated heat as a function of the detuning ∆ (b) and the drive phase ϕ (c).The initial state is |0 in (a, b) and 1 √ 2 (|0 + |1 ) in (c).A high temperature k B T /ω = 0.5 is considered in (c).All other parameters are as in figure 1.