Nonequilibrium quantum fluctuation relations for harmonic systems in nonthermal environments

We formulate exact generalized nonequilibrium fluctuation relations for the quantum mechanical harmonic oscillator coupled to multiple harmonic baths. Each of the different baths is prepared in its own individual (in general nonthermal) state. Starting from the exact solution for the oscillator dynamics we study fluctuations of the oscillator position as well as of the energy current through the oscillator under general nonequilibrium conditions. In particular, we formulate a fluctuation-dissipation relation for the oscillator position autocorrelation function that generalizes the standard result for the case of a single bath at thermal equilibrium. Moreover, we show that the generating function for the position operator fullfills a generalized Gallavotti-Cohen-like relation. For the energy transfer through the oscillator, we determine the average energy current together with the current fluctuations. Finally, we discuss the generalization of the cumulant generating function for the energy transfer to nonthermal bath preparations.


Introduction
Fluctuation relations [1][2][3][4][5][6][7] build on the fundamental connection between the response of a physical system to a weak externally applied force and the fluctuations in the system without the external force. This connection was first observed for thermal equilibrium by William Sutherland [8,9] and Albert Einstein [10][11][12]. They established the relation between the mobility of a Brownian particle, which is a quantity that measures the response to an external electric field, and the diffusion constant, which is a quantity that characterizes the fluctuating forces at equilibrium. The famous Johnson-Nyquist relation [13,14] gives the corresponding connection between the electrical resistance of a circuit and charge fluctuations in the resistor. A more general relation has been derived by Callen and Welton [15] in form of the quantum fluctuation-dissipation theorem (FDT) which relates the Fourier transform Ψ(ω) of the symmetric equilibrium correlation function of an observable to the Fourier transform Φ(ω) of the (antisymmetric) response function of this observable in thermal equilibrium at temperature T = (k B β) −1 . It was recognized by Green [16,17] and Kubo [18] that the FDT in Eq. (1) is a particular case of the more general linear response theory which is an invaluable tool to model and understand experimental data in all fields of physics. However, often situations are encountered where the assumption of thermal equilibrium is invalid, for example, for systems strongly driven by external fields, charge currents in systems with large differences in the electric potential, heat currents in systems with strong temperature gradients, or systems in solvents and disordered media which themselves are in metastable quasi-equilibria only. It has been a longstanding task in statistical physics to generalize linear response theory and FDTs to such nonequilibrium situations and, by this, to build a unifying theoretical framework of the spectral characteristics of environmental noise. Generalized nonequilibrium fluctuation theorems have been formulated for classical nonstationary Markov processes [19] and for stationary Markov processes far away from thermal equilibrium [20,21]. They relate the higher-order nonlinear response to higherorder correlation functions of stationary nonequilibrium fluctuations. A fully nonlinear, exact and universal classical fluctuation relation has been provided by Bochkov and Kuzovlev [22]. It gives the fluctuation relation at any order for systems that are in a thermal state in absence of external forces. It solely builds on the time-reversal invariance of the equations of motion and the assumption of a thermally equilibrated initial state. The quantum version was provided by Andrieux and Gaspard [23] and lead to fundamental insights [1] into the fact that work injected to or extracted from a system is not a quantum mechanical operator or observable, because it characterizes a process rather than a state of the system [24].
Recently, growing interest in nonequilibrium fluctuation relations arose from alternative formulations by Evans et al. [25] and by Gallavotti and Cohen [26] for the statistics of nonequilibrium fluctuations in steady states and by Jarzynski [27] and Crooks [28] on the statistics of work performed by a transient time-dependent perturbation [1]. The reviews [1][2][3][4][5][6][7] summarize the actual progress in this field.
Most studies so far consider systems initially in thermal equilibrium, described by the canonical distribution with the system Hamiltonian H 0 and the partition function Z 0 = Tr[e −βH 0 ]. In this work, we want to give up this assumption and formulate generalized nonequilibrium fluctuation relations for nonthermal initial states. To do so, we consider the dissipative quantum mechanical harmonic oscillator [29][30][31][32][33][34][35][36][37]. Building on our previous work in Ref. [38] we study a central oscillator coupled to an arbitrary number of harmonic baths each of which can be prepared in its own individual initial state. The fluctuations of the baths are thus still Gaussian, but not necessarily thermally distributed. Because the exact solution for the system dynamics is known, we can analytically calculate all observables and correlation functions of interest, and thus investigate the validity of nonthermal nonequilibrium fluctuation relations for this admittedly restricted model situation.
The structure of the paper is as follows. We introduce the model, its classical equation of motion and the basic notions in Sec. 2. Then, in Sec. 3, we calculate the symmetric and antisymmetric correlation functions of the oscillator position for the case of general nonthermal bath states. In Sec. 3.3, we formulate the generalized nonequilibrium fluctuation relation for the oscillator position correlation functions. This constitutes one major result of this work. In Sec. 4, we calculate the generating function for the position operator of the oscillator and show that it fullfills a generalized Gallavotti-Cohen relation under nonequilibrium conditions at arbitrary times. Sec. 5 is devoted to energy transfer and we present the derivation of the average energy current. In Sec. 6, we calculate the energy current fluctuations and generalize the well-known cumulant generating function of the heat transfer for thermal baths to general bath preparations, before we summarize in Sec. 7.

The model
In a system-bath model approach, we consider the one-dimensional harmonic oscillator bilinearly coupled to a finite number N B of different and mutually uncoupled baths of harmonic oscillators. The total Hamiltonian is H = H S + H B + H SB , where ( = 1, k B = 1 throughout the work) is the contribution of the central oscillator with frequency Ω, is the contribution of the bath oscillators with frequencies ω α ν , and is the coupling part. In these expressions, the position and momentum operators Q α ν and P α ν fullfill the canonical commutation relation [Q α ν , P α ′ µ ] = iδ νµ δ αα ′ . The labels α, α ′ = 1, . . . , N B are used to identify a particular bath, while the indices ν, µ = 1, . . . , N α identify a single oscillator from bath α.
The coupling term contains the counter term which serves to eliminate the potential renormalization due to the coupling of the oscillator to the baths [39,40]. Throughout this work, we assume factorizing initial states ρ(0) = ρ S (0) N B α=1 ρ α B (0) corresponding to the choice of isolated systems that are brought into contact at t = 0 + . Notice, however, that we keep the initial distributions ρ α B (0) of the baths arbitrary and do not necessarily assume thermal equilibrium.

The exact solution for the operator dynamics
Starting from the Heisenberg equation of motion for the system and bath operators, one inserts the formal solution for the bath operator dynamics into the equation of motion of the central oscillator to obtain the quantum Langevin equation with the damping or friction kernel and the noise term The noise term η(t) together with the initial slip term K(t)Q(0) appears as a fluctuating force in Eq. (6). Due to our choice of factorizing initial states, the noise terms of different baths are uncorrelated, i. e., η α (t)η β (s) = η α (t) η β (s) for α = β. Nevertheless, the fluctuating forces ξ α (t) = η α (t) + K α (t)Q(0) including the initial slip term are correlated because of the coupling to the central oscillator [41]. These correlations vanish if the expectation values are calculated with respect to the non-factorizing initial state that is obtained out of ρ(0) through the unitary transformation with the displacement operator exp[iQ Nα ν=1 λ α ν P α ν /(ω α ν ) 2 ]. At this point, we note that explicit expressions for the correlation functions of the fluctuating forces depend on the choice of the initial distributions ρ α B (0) and thus, the fluctuations are in general associated with a nonstationary Gaussian operator noise. Only in the limit of long times, these fluctuations become stationary again (see Appendix A).
As is well established in the literature [29][30][31][32][33][34][35][36][37], the full solution for the central oscillator dynamics can be constructed from the solution u(t) ∈ R of the corresponding classical equation of motion, The relevant solution u(t) is specified by u(t) = 0 for t < 0 and by the initial conditions u(0) = 0 andu(0) = 0. It is given by the Fourier transform of the function Given u(t), the solution for the dynamics of the central oscillator operators can be obtained from the matrix equation We here introduced the matrices and denote by the respective index R or I the real or imaginary part of the partial Fourier transforms of the classical solution u(t),

Expectation values
Equation (12) allows us to express central oscillator expectation values for t ≥ 0 in terms of the initial ones at t = 0. The linear expectation values are given by the equation where depends on the initial bath expectation values For the quadratic expectation values we define the correlator of two operators A and B by and write Σ AB (t) ≡ Σ A(t)B(t) for better readability. For correlators of operators related to bath oscillators at initial time, we define and write We then obtain with Eq. (12) the relation where

The thermodynamic limit
In the thermodynamic limit N α → ∞ for all α = 1, . . . , N B we can replace summations by introducing the densities of states of the baths, that converge to continuous functions. Since the coupling constants λ α ν enter Eqs. (7) and (11) as (λ α ν ) 2 , they have to scale as 1/ √ N α to obtain finite results for the sum over N α terms. We thus introduce continuous functions λ α (ω) according to and define the bath spectral functions Note that we here use the definition of the bath spectral function of Ref. [42] without the factor π/2, which corresponds to the definition of Ref. [30] with an additional 1/ω factor. The linear expectation values X α ν have to scale as 1/ √ N α , because they appear in Eq. (18) with the prefactors λ α ν . We introduce continuous functions X α,Q (ω) and X α,P (ω) according to Moreover, we have to separate the N α diagonal terms Σ α νν from the N 2 α off-diagonal terms Σ α νµ with ν = µ that require an additional 1/N α prefactor for convergence in the thermodynamic limit. Hence, we define with continuous functions σ (1) α,XY (ω) and σ (2) α,XY (ω 1 , ω 2 ) (X, Y = Q, P ) as the matrix entries of Σ The function F (z) in the thermodynamic limit can be obtained via contour integration with the result for Im z > 0. The complex functions Γ α (z) follow from analytic continuation of γ α (ω) = ∓(2/π) Im Γ α (±ω + i0 + ) into the upper half of the complex plane. If the function F (z) has no poles for Im z > 0, the classical function u(t) from Eq. (10) is the inverse Fourier transform of a continuous function. We can use the Riemann-Lebesgue lemma valid for any integrable function f (ω) and conclude, that u(t) → 0 in the long-time limit t → ∞. In turn, poles of F (z) correspond to undamped oscillations in u(t), such that the central oscillator will approach a stationary state only if isolated modes do not exist. The possibility of lim t→∞ u(t) = 0, i. e. the existence of isolated poles in F (z), is closely connected to a breaking of ergodicity in the sense of the mean-square of a stochastic observable [43][44][45][46]. Precise conditions for lim t→∞ u(t) = 0, as well as a general discussion of equilibration and thermalization of the central oscillator, can be found in Ref. [38]. Throughout this work, we assume that F (z) has no isolated poles, such that the classical solutions for t → ∞ approach zero, i. e., U(t) → 0. Then, the central oscillator equilibrates and the asymptotic state is Gaussian with the expectation values in the long-time limit lim t→∞ X(t) = 0 and lim t→∞ Σ(t) = Σ ∞ [38].

Nonequilibrium fluctuation relation for the oscillator position
The results from the previous section allow us to derive a generalized nonequilibrium fluctuation relation of the form of Eq. (1). For this, we determine the symmetric and the antisymmetric correlation functions of the central oscillator position Q. Their Fourier transforms are then shown to obey a generalized nonequilibrium fluctuation relation in form of a characteristic proportionality relation.

The symmetric correlation function
We define the symmetric correlation function of the central oscillator position as Inserting the solution for Q(t) from Eq. (12) and performing the thermodynamic limit N α → ∞, we obtain with the two functions and In the long-time limit t → ∞ the terms involving u(t),u(t) and Q(t) vanish according to our assumption of continuity of F (z). For the remaining terms Ψ (1) (t, s) and Ψ (2) (t, s) we rewrite the partial Fourier transform of Eq. (15) Since u(t) vanishes at long times, the partial Fourier transform u(t + s, ω) behaves asymptotically as where is the full Fourier transform of the function u(t) ‡. Using this asymptotic behaviour in the expressions for Ψ (1) (t, s) and Ψ (2) (t, s) we see that the off-diagonal term Ψ (2) (t, s) contains only oscillatory terms in the two frequencies ω 1 and ω 2 . If we recall the Riemann-Lebesgue lemma, Eq. (30), we conclude, that Ψ (2) (t, s) vanishes in the longtime limit. Following the same line of reasoning we find that the only non-zero term in the limit t → ∞ comes from Ψ (1) (t, s) and involves |u(ω)| 2 while the arising oscillating terms vanish. In particular, where denotes the frequency-resolved energy distribution functions of the initial bath states. We next Fourier transform Eq. (37) and obtain (39) ‡ We use the same symbol u for the function and its Fourier transform for ease of readability. Time arguments are denoted as t, τ or s, while frequency arguments are denoted by ω.

The antisymmetric correlation function
The antisymmetric correlation function of the oscillator position Q is given by Inserting the solution Q(t) of Eq. (12), using the property [Q α ν (0), P α µ (0)] = iδ νµ and performing the thermodynamic limit N α → ∞, we obtain which is independent of the initial bath preparation as expected [39]. Similar to the calculation of the symmetric correlation function, we obtain for the antisymmetric response function in the long-time limit Its Fourier transform readily follows as

The generalized nonequilibrium fluctuation relation
To formulate the general nonequilibrium fluctuation relation, we compare Eqs. (39) and (43) and obtain for general initial preparations and an arbitrary number N B of independent harmonic baths the relation This is one major result of the present work and illustrates that the relation is crucially determined by the frequency-resolved energy distributions E α (ω) of the initial bath states defined in Eq. (38) and the bath spectral functions γ α (ω) given in Eq. (26). A comparison with the thermal fluctuation-dissipation theorem in Eq. (1) shows that in the considered nonthermal situation we have to exchange the thermal energy distribution with the average of the individual energy distributions of the baths weighted with their spectral functions.
In the case when all baths are initially distributed thermally at the same temperature T according the thermal equilibrium Bose-Einstein distribution function, we have E α (ω) = E th (ω, T ) for all α = 1, . . . , N B . This reproduces the equilibrium fluctuation-dissipation theorem Eq. (1) [39].
A natural question then is under which initial bath preparations the central oscillator thermalizes, i. e., reaches a stationary state which is thermally distributed with a given temperature T . By comparing Eqs. (44) and (1), we obtain the condition for which the fluctuations of the central oscillator for t → ∞ are thermal. This condition certainly is satisfied whenever all baths are thermal and have equal temperature, but can also be satisfied for other nonthermal initial bath preparations. In turn, if this condition is satisfied, the quantity is a constant, i. e., independent of ω. It is then tempting to understand this quantity as an "effective" temperature characterizing the general initial bath preparation. However, the above condition does not guarantee true thermalization of the central oscillator, which is essential for a meaningful notion of temperature. For a more detailed discussion of this question, see Ref. [38].

Generating function for the position operator of the oscillator
In this section we show that the dissipative oscillator model allows us to study the connection between transient and steady state fluctuation relations. We calculate the generating function for the central oscillator position operator and show that it fullfills a Gallavotti-Cohen symmetry relation [26] valid for arbitrary times and a Gaussian initial state of the central oscillator. This additional Gaussian assumption is not necessary in the long-time limit and we obtain an exact result for the steady state fluctuation relation.
We define the generating function for the position operator according to With that, all the cumulants Q n (t) of the position operator follow by performing the respective derivative, For instance, we have Q(t) = Q(t) and Q 2 (t) = Σ QQ (t).
It is convenient to represent the generating function in terms of the Wigner function of the central oscillator where we write W S (x, t) = W S (q, p, t) with x = (q, p) T and dx = dq dp for abbreviation. The Wigner function W S (x, t) at time t ≥ 0 can be obtained from the propagating function J W (x,x, t) = J W (q, p,q,p, t) in Wigner representation, that is defined by the relation and can be evaluated to [38] with U(t), I(t), and C(t) from Eqs. (13), (18), and (23). Performing the Gaussian integral over x we obtain where e 1 = (1, 0) T . In the long-time limit t → ∞, where U(t) → 0 according to our assumption of continuity of F (z), the integration in Eq. (54) evaluates to one because the initial Wigner function is normalized. We then obtain with Σ ∞ QQ = lim t→∞ Σ QQ (t). The results obeys the symmetry Z ∞ Q (ξ) = Z ∞ Q (−ξ). For finite times, we can restrict ourselves to Gaussian initial states of the central oscillator, and obtain In order to see when the Gallavotti-Cohen relation is fullfilled, we calculate Hence, the relation Z Q (−ξ + iA, t) = Z Q (ξ, t) is fullfilled at any arbitrary time t, if This implies that the oscillator fluctuates symmetrically around its momentary position average Q(t) since Z Q− Q (−ξ, t) = Z Q− Q (ξ, t). On the other hand, however, the symmetry point for the generating function of the position operator, which in the stationary state is ξ = 0, is shifted by the momentary position expectation value scaled by the momentary position variance, i. e. Z Q (−ξ + iA/2, t) = Z Q (ξ + iA/2, t). Note that this relation holds in general and also when the central oscillator has not yet reached its equilibrium state. This transient fluctuation relation is linked with the steady state fluctuation relation from above by realizing that lim t→∞ A(t) = 0.

Quantum mechanical energy transfer between nonequilibrium baths
We now study the quantum mechanical transfer of energy between nonequilibrium baths. To keep the discussion simple, we concentrate on the case of the energy transfer between two baths, i. e., N B = 2, and denote them as left (α = l) and right (α = r) reservoir. In particular, we are interested in the form of the expectation value of the energy current operator which can be defined for instance for the left junction according to [47][48][49][50][51][52] with the anticommutator defined as {A, B} = AB + BA.
For the calculation of the expectation value I(t) we need the solutions of the Heisenberg equations of motion for the left bath operators, We insert these equations and the solution Eq. (12) for Q(t) into Eq. (60) and perform the thermodynamic limit to obtain I(t) = I 1 (t) + I 2 (t) + I 3 (t) with In these equations, Ψ(t, τ ) is the symmetric position autocorrelation function given in Eq. (31) and the diagonal and non-diagonal contributions to I 1 (t) and I 2 (t) are and l,P P (ω 1 , ω 2 ) . (64b) To perform the long-time limit, we follow the line of reasoning of Sec. 3.1. The terms containing the linear expectation value Q(t) vanish. The off-diagonal terms I The remaining diagonal terms can be simplified algebraically using the asymptotic behaviours Eq. (35) of the partial Fourier transform u(t, ω) and Eq. (37) of the symmetric correlation function Ψ(t, τ ) and by applying the Riemann-Lebesgue lemma Eq. (30). We finally obtain the expectation value of the energy current from the left reservoir to the central oscillator in the long-time limit as We can rewrite this expression into the final form by using, that the Fourier transform u(ω) in Eq. (36) is the inverse of the Fourier transform in Eq. (10), such that u I (ω) = − Im F (ω+i0 + ) = −(π/2)[γ l (ω)+γ r (ω)]|u(ω)| 2 . Expression (66) generalizes Eq. (4.2) of Ref. [49] and reproduces it for the special case of thermal baths. Obviously, the asymptotic energy current vanishes exactly, if γ r (ω) = 0 for only one bath, or if E r (ω) = E l (ω) for equal bath preparations.
For two thermal baths with E α (ω) = E th (ω, T α ), and T l = T r + ∆T where ∆T ≪ T l , we can expand the energy distribution function as where sinh −1 x = 1/ sinh x. Thus, we obtain the linear response result growing linearly with the difference ∆T of the temperatures of the left and right bath.

Nonequilibrium fluctuations of the transferred energy
In this section, we consider the energy which is transferred from one bath (say, the left) to the central oscillator in presence of the second bath (say, the right). Moreover, we are interested in the fluctuations of the transferred energy. We note in passing that we use the more general term of "energy" instead of "heat" since the definition of heat in the strict sense requires purely thermal environments. The energy that is transferred from the left bath to the rest of the system until time t is obtained from the difference of the energy of the left bath between times t and 0. This involves the measurement of the observable H l B at two different times. Following the idea of two-time quantum measurements, the corresponding generating function can be written as [50][51][52] Z where the prime indicates that the expectation value has to be taken with respect to the projected density matrix Here |φ a is an eigenstate of the operator H l B , i. e., H l B |φ a = a|φ a . Writing the generating function as a series in powers of iξ, we obtain [50,52] ln where W n (t) denotes the nth order cumulant of the operator In the following, we calculate the moments of the energy transfer operator W (t) entering Eq. (71). In particular, we are interested in the long-time limit of these quantities.

The first moment
Using Eqs. (61a) and (61b) the linear expectation value of the energy transfer operator follows as In the long-time limit t → ∞, we expect from the definition in Eq. (72) and from the result I(t) → I ∞ of the last section that W (t) grows linearly with time. It is thus useful to consider W (t) /t instead of W (t) .
The explicit calculation of W (t) /t in the long-time limit is achieved by inserting the solution Q(t) from Eq. (12), performing the thermodynamic limit according to Sec. 2.3, and analytically carrying out the remaining time integrations. The result is As expected, this expression coincides with the expectation value of the energy current operator given in Eq. (66).

The second moment
The second moment of the energy transfer operator is Expanding the square yields a sum of terms containing expectation values of products of four operators. We may reorder the operator products using the commutators A general expectation value of a product of four operators can be ascribed to a sum of products of expectation values of one or two operators for Gaussian states. The assumption of a Gaussian bath state is justified in the thermodynamic and long-time limit on general grounds [53,54]. In Ref. [38] it is shown that the state of the central oscillator becomes Gaussian for t → ∞, independent of its initial preparation, if the classical solution u(t) vanishes asymptotically-the situation of interest here.
For an explicit result, we insert the solution Q(t) from Eq. (12), perform the thermodynamic limit, use the results for the position correlation functions from Sec. 3, and carry out the remaining time integrals in the long-time limit. The result for the second order cumulant then reads This expression generalizes the result for the second moment given in Eq. (9) in Ref. [50,51]. Eq. (77) reduces to this equation for the special case of thermal baths with E α (ω) = E th (ω, T α ) = ωf α (ω) + ω/2 which then lead to the expressions f α (±ω) = 1/[exp(±ω/T α ) − 1] in Ref. [50,51].

The generating function for the energy transfer
We have seen that the well-known results [50,52] for the first and second moment of the heat transfer operator for the special case of thermal baths are well reproduced by our more general results. The generalization follows by the corresponding replacements of the thermal distribution functions of the baths by the general initial distributions. Hence we can now follow the same line of reasoning and generalize the steady state expression of the cumulant generating function for the heat transfer given in Eq. (8) in Ref. [50,51] with the result We observe that G(ξ) fullfills the symmetry relation where A = β r − β l with Since the constants β α should be independent of ω the existence of the symmetry (79) implies a condition on the initial bath preparation. In particular, the energy distribution functions should be thermal, i. e. E α (ω) = E th (ω, T α ). Note that this is a condition on the combination E α (ω) of the initial bath variances σ (1) α,QQ (ω) and σ (1) α,P P (ω), not on the individual functions [see Eq. (38)]. It can be fullfilled for nonthermal bath preparations as well [38].
From the relation (79) it follows that the probability distribution of the transferred energy, fullfills the steady state fluctuation theorem We remark that the exchange fluctuation relation (82) can only be proven rigorously when the initial preparation is indeed free of correlations and also the interaction of the system and the bath is switched off at some final time [1,55,56]. The role

Summary
Most studies related to fluctuation relations so far consider the special case when the systems are initially in thermal equilibrium, but do not restrict their analyses to a specific model. In the present study, we give up the assumption of initial thermal states and allow for nonthermal bath preparations. The price we have to pay for this generalization is the restriction to an analytically solvable model for which we obtain exact generalized nonequilibrium fluctuation relations. On the one hand, we can give the explicit expressions for the symmetric and antisymmetric autocorrelation functions of the central oscillator position. Then, a generalized nonequilibrium fluctuation relation follows which only involves the bath spectral functions and the frequency-resolved energy distribution of the initial bath states. The general expression also contains the special case of a single thermal bath and coincides with the well-known equilibrium fluctuation-dissipation theorem. Moreover, we discuss the conditions under which the generating function of the oscillator position fullfills a Gallavotti-Cohen relation at arbitrary times. This relation reflects the fact that the oscillator position fluctuates symmetrically around its momentary average position. On the other hand, we have elucidated the quantum mechanical energy transfer through the central oscillator by calculating the time-dependent energy current and the second moment of the current fluctuations. Based on this result we generalize the cumulant generating function for energy transfer, which is well-known for thermal baths, to the nonthermal situation.