Fluctuation-dissipation theorem for non-equilibrium quantum systems

The fluctuation-dissipation theorem (FDT) is a central result in statistical physics, both for classical and quantum systems. It establishes a relationship between the linear response of a system under a time-dependent perturbation and time correlations of certain observables in equilibrium. Here we derive a generalization of the theorem which can be applied to any Markov quantum system and makes use of the symmetric logarithmic derivative (SLD). There are several important benefits from our approach. First, such a formulation clarifies the relation between classical and quantum versions of the equilibrium FDT. Second, and more important, it facilitates the extension of the FDT to arbitrary quantum Markovian evolutions, as given by quantum maps. Third, it brings out the full connection between the FDT and the Quantum Fisher information, the figure of merit in quantum metrology.

looks for correction terms in the original equilibrium FDT [20,21], whereas the second keeps the very mathematical structure of the theorem by redefining the magnitude conjugated to an external parameter [22,23].
Here, we adopt the second strategy to prove a FDT for generic quantum Markovian systems. The key point in our derivation is the use of the symmetric logarithmic derivative (SLD), Λ λ , of a density matrix ρ λ depending on a real parameter λ, defined as: The SLD is an observable with zero average, Λ λ λ = Tr[Λ λ ρ λ ] = 0, as can be easily proved by taking the trace of the above equation. It is intimately related to the quantum Fisher information (QFI), F λ = Tr Λ 2 λ ρ λ , which plays a prominent role in metrology, since the uncertainty of any unbiased observable A (i.e. with A λ = λ), satisfies the Crámer-Rao bound Var(A) λ ≥ 1/F λ [24][25][26][27][28][29][30]. In this paper we show that the SLD provides a novel definition of an observable conjugated to an external parameter, which is extremely useful to derive a completely general FDT for quantum Markov systems and to relate previous versions of the FDT for classical and quantum systems.
We start by applying the SLD to the simplest case of a fluctuation-dissipation relation for the static susceptibility. Consider a quantum system whose density matrix ρ λ depends on an external parameter λ. Taking ρ 0 as a reference state, we are interested on the change of the expected value of a generic observable B under a small change in λ. More precisely, for small λ, the expected value of a generic observable B can be written as: where is the static susceptibility of observable B. Using the SLD, the derivation of a fluctuation-dissipation relation is straightforward: which is the symmetrized correlation between observables B and Λ 0 , since Λ 0 0 = 0. In the Appendix A we show that Eq. (4), when particularized to a thermal state ρ λ = e −β(H0−λA) /Z(λ), with β = 1/(kT ) and Z(λ) ≡ Tr e −β(H0−λA) and expressed in the eigenbasis of H 0 , yields the standard fluctuation-dissipation relation. We now turn to the case of a generic Markov evolution given by the composition of completely positive and trace preserving (CPTP) maps. Let ξ λ (ρ) be a CPTP map that depends on a parameter λ. We assume that each map ξ λ has an invariant state π λ , i.e., ξ λ (π λ ) = π λ .
We study a small time-dependent perturbation λ(t) affecting the invariant state π 0 . More precisely, we consider the evolution, in discrete time steps t = 1, 2, . . . , The linear response of an observable B can be written as: where is the response function of the observable B under the perturbation λ.
One can extend the above definition to the case of maps acting for a short time ∆t. In the continuous limit, ∆t → 0, the sum in (6) is replaced by the integral The generalized susceptibility is defined as the Fourier transform of the response function (φ B (t) is assumed to vanish for t < 0 due to causality): The generalized susceptibility has interesting properties such as the Kramers-Kronig relation between the real and imaginary parts χ

When the evolution is unitary under the Hamiltonian
is called absorptive part of the susceptibility [2][3][4][5], since the energy absorbed by the system due to the perturbation is proportional to χ B (ω). The static susceptibility χ s B can be related to the response function by considering a constant perturbation λ(t) = λ for t ≥ 0 [2][3][4][5]. In this case B(t) → B λ when t → ∞ and Eq. (7) implies that the static susceptibility is the integral of the response function or, equivalently, the generalized susceptibility at zero frequency: χ s B = χ B (ω = 0). To obtain the FDT we calculate ρ(t) up to linear terms in λ. For that, we write ξ λ = ξ 0 + λξ 1 + . . . (notice that ξ 1 is not a CPTP map) and the invariant state as π λ = π 0 + λπ 1 + . . . The invariance of π λ under the map ξ λ implies: ξ 1 (π 0 ) + ξ 0 (π 1 ) = π 1 (9) and the SLD of π λ with respect to λ at λ = 0 obeys Expanding the evolution equation (5) up to linear terms, we obtain where we have used the invariance of π 0 under ξ 0 . The expected value of B is Comparing (12) with (6), we immediately get Using (9) and (10): where ∆B(t) = B(t + 1) − B(t) and B(t) =ξ t 0 (B) is the evolution of the observable B in the generalized Heisenberg picture for quantum maps [31][32][33]. Hereξ 0 (·) is the adjoint map (not necessarily trace preserving) with respect to the scalar product between operators given by the trace, i.e., Tr[Aξ 0 (B)] = Tr[ξ 0 (A)B], for all pair of operators A and B.
The fluctuation-dissipation relation for the static case (4) is recovered from (14) if the correlations between B(t) and Λ vanish in the limit t → ∞: Finally, the continuous-time version of theorem (14) is Eqs. (14) and (16) are our main result. These results are the quantum counterparts of the non-equilibrium classical FDT derived by Agarwal in [13] and revived recently in [22,23]. Notice that, in the classical scenario, the conjugate variable is defined as the derivative of the logarithm of the steady state probability distribution. On the quantum scenario, the non-commutativity of observables does not allow to simply replace the classical conjugate variable with the derivative of the logarithm of the density matrix. Nevertheless, the choice of the SLD observable solves this non-commutativity issue.
In the Appendix B we prove that the latter expression (16), when particularized to Hamiltonian dynamics and equilibrium states of the form π λ = e −β(H0−λA) /Z(λ), is completely equivalent to the standard Kubo formula [2][3][4][5]: It is worth it to point out that, in the quantum case, the Kubo formula does not yield a simple FDT; more precisely, the response function cannot be expressed in terms of the time derivative of a two-time correlation in equilibrium. Such a relationship can only be derived in the frequency domain for the absorptive part of the generalized susceptibility and the Fourier transform of the correlation (see [3][4][5] and the Appendix B). On the contrary, our generalized FDT, namely Eqs. (14) and (16), expresses the response function in terms of correlations in the time domain and can be equally applied to classical and quantum systems. This uniform formulation is possible due to the introduction of the SLD. In the classical case, the SLD coincides with the normal derivative and consequently, for a thermal state with Hamiltonian , whereas for a quantum system with [H 0 , A] = 0, the SLD yields a nontrivial conjugated variable as shown below. We highlight the usefulness of our FDT for quantum metrology explicitely through the examples that follow. However, let us remark that such a link holds for any map that fits our framework. We illustrate the new generalized FDT with a simple example consisting of two harmonic oscillators with a modulated interaction. The Hamiltonian reads: where a i and a † i are the ladder operators of i-th oscillator and δ denotes the frequency detuning between the oscillators. We assume |J(t)| ω so that linear response theory holds at any time. First, we consider the two harmonic oscillators placed in a bath at temperature T , and examine the response of the system to a perturbation J(t) = J 0 around the thermal equilibrium state ρ 0 = exp(−βH 0 )/Z. Classically, the susceptibility is defined as being the conjugate variable. Notice that if the detuning is zero, then [H 0 , A] = 0, and the SLD reads [34]. On the contrary, any δ = 0 forces that [H 0 , A] = 0, and the SLD cannot be anymore identified trivially as the conjugate variable A [12]. For a finite detuning, the SLD takes the form [35,36]). This additional coefficient arises from the non-commutativity between H 0 and A, and reads: The use of the SLD allows us to show that the QFI is: The additional coefficient is bounded, 0 < C(δ) ≤ 1. It has a maximum at δ = 0, and then it decreases monotonically with the detuning, therefore, the precision on the estimation of J decreases as the detune increases. We now address the problem of a time-dependent modulation of J(t) in a non-equilibrium environment induced by two thermal baths at different temperatures T 1 and T 2 . A master equation that has been widely used to describe the reduced state of the oscillators consists of a Linblad equation [31,[37][38][39] with two independent dissipators D i Here is the mean occupation number of the i-th oscillator, and γ is the dissipation rate, which is assumed to be equal for both oscillators. The equation is only valid for small J(t)-compared to ω, i.e., the system's energy scale-and J(t) < δ [40]. In particular, it does not predict the thermalization of the full two-oscillator system for finite J when T 1 = T 2 . Obviously, for J(t) = 0, both oscillators evolve independently reaching a stationary state in which each oscillator is at thermal equilibrium with its own bath: where H i denotes the free Hamiltonian of the i-th oscillator. Since the model under study is quadratic in creation and annihilation operators, one may equivalently describe it by means of its covariance matrix (CM) which contains only first and second moments (i.e., Gaussian), notably simplifying the calculations. The latter is described with the help of the quadratures which satisfy the standard commutation relation [x j , p k ] = iδ jk . In turn, the matrix elements of the CM are defined as follows: The CM corresponding to the stationary state of the master equation (21) is provided in the Appendix C, where we also find (i) the SLD, and (ii) the time evolution of all the quadratic observables. Note that according to Eq. (16), (i) and (ii) are the two key elements required for evaluation of φ B (t). Specifically, the SLD writes as: with c 1 and c 2 being real numbers. Therefore, the SLD is a non-local operator, hence the response of any local observable vanishes.
To proceed further, let us consider J(t) = J 0 (1 − cos ν t). It is convenient to define the time dependent counter part of the static susceptibility as follows: which quantifies the deviations of the observable from its initial value, and does not depend on J 0 . Note that for a constant perturbation (i.e., ν = 0) we recover χ B (t → ∞) = χ s B . In Fig. 1 (a) we depict χ B (t) versus time for B = x 1 x 2 evaluated at three different modulation frequencies ν = 0, ν = δ/2, and ν = δ. For the constant perturbation, we observe that at short times χ B (t) is oscillating but, as time passes, the system relaxes into a new steady state such that B(t → ∞) = B J . On the contrary, for ν = 0, the system never relaxes to a stationary state. In this case, the time dependent susceptibility can be approximated at sufficiently large times by: where α = arctan Re χ B (ν)/Im χ B (ν) . Effectively, the generalized suceptibility |χ B (ν)|, maps onto the amplitude of the oscillations of the time dependent susceptibility, and its dependence on ν is illustrated in Fig. 1 (b). One can see that the dynamical susceptibility is a flat function for ν δ, increases sharply as ν turns out to be resonant with the detuning, and strongly decreases afterwards. A similar behaviour is observed for any other nonlocal observable B = {x 1 p 2 ; p 1 x 2 ; p 1 p 2 }. Futhermore, the metrological character of the FDT when expressed through the SLD, can be seen by identifying which is the best measurement to detect the smallest interaction coupling J(t) = J 0 . If we have no prior knowledge about J 0 , linear response theory ensures that T r[B(ρ(t) − ρ 0 )] = J 0 χ B (t). Unavoidable, an expectation value bears a statistical error, Var(B) 0 /n, by repeating the measurement n times 1 . Thus, in order to infer the value of the perturbation, the linear response of the system must be larger than the statistical error, i.e., |J 0 χ B (t)| > Var(B) 0 /n, setting a lower bound on the smallest perturbation that one may detect by measuring B. The inverse of this lower bound defines the sensitivity F (B) t of the observable B to the perturbation J 0 : Notice that the above definition is not restricted to the model discussed here. It can be applied for estimation of a generic parameter which is parametrized through a divisible map as in (5). For a constant perturbation, as t → ∞, χ B (t → ∞) ≈ Corr(B, Λ 0 ) 0 . Substituting this in Eq. (28), and using the Cauchy-Schwartz inequality, the upper bound on the sensitivity of B at t → ∞ is given by: with F 0 being the QFI evaluated at J 0 = 0. The QFI has been previously used to quantify the ultimate precision of parameter estimation in non-equilibrium steady states of spin models [41]. The bound (29) is saturated by performing a measurement in the Λ 0 basis. In Fig. 2, we depict F (B) t for three observables B = {Λ 0 ; x 1 x 2 ; x 1 p 2 }. Although no observable can overperform the sensitivity of the SLD at t → ∞, at shorter times this fails to be the case as shown in the inset of Fig. 2. More interestingly, the maximum value of F (B) t is not necessarily achieved at t → ∞. A similar behavior is observed for a time-dependent perturbation J(t) = J 0 (1 − cos νt). Although in this case inequality (29) does not apply, the FDT as expressed in Eq. (16) links Λ 0 with the response function. From the oscillatory character of the time dependent susceptibility, the optimal times to best estimate J 0 are given by Eq. (27).
In summary, we have presented a novel formulation of the Fluctuation-Dissipation Theorem (FDT) in terms of the symmetric logarithmic derivative (SLD), which is completely general to describe the effect of a perturbation on a quantum systems, within the linear re-  sponse. Such a formulation presents some clear advantages. First, it unifies FDT for classical and quantum systems. Second, it can be straightforwardly applied to any Markovian evolution by means of quantum maps. This permits the extension of the FDT to nonequilibrium dynamics. Third, it provides an explicit connection between the susceptibility of an observable to an external perturbation and the figure of merit in quantum metrology.
Our FDT can be used to generalize the detection of multipartite entanglement to the non-equilibrium steady states (NESS). Although at thermal equilibrium, or after quenching a thermal state, such relations are known [7,8], there are less studies for a general NESS state. In this regard, the sensitivity measure Eq. (28)which is a lower bound on the QFI [28,42]-can be used to detect multipartite entanglement. Finally, our results prompt an interesting open question on the relationship between the quantum versions of FDT and fluctuation theorems (FT). Needless to mention that, the quantum versions of FT have been subject of extensive studies, that have produced a rich literature on the subject [43][44][45]. The relationship between FDT and FT has been made clear for classical systems. On the other hand, the FT for quantum CPTP maps requires a condition that is not necessary in our derivation of the FDT, namely, the Kraus operators of the map must be ladder operators in a relevant basis: the eigenbasis of the instantaneous stationary state π 0 for generic CPTP maps [46] or, for periodically driven systems in contact with equilibrium reservoirs, Floquet eigenstates [47] or displaced energy eigenstates [48]. Our FDT suggests that one could obtain a more general FT by using the SLD.

A Quantum systems in thermal equilibrium
A relevant particularization of the static fluctuation dissipation relation Eq. (4) is the case of a quantum system with Hamiltonian H 0 − λA at thermal equilibrium. In such a case the density matrix is ρ λ = e −β(H0−λA) /Z(λ), where β = 1/(k B T ) is the inverse temperature and Z(λ) ≡ Tr e −β(H0−λA) is the partition function of the system. The equilibrium static susceptibility of an arbitrary observable B under the perturbation λA is denoted as χ s B and obeys Eq. (4). Furthermore, when the SLD is expressed in the eigenbasis of the unperturbed Hamiltonian H 0 |n = E n |n , one can obtain some interesting relationships for the equilibrium static susceptibility. It is convenient to rewrite the SLD using the Feynman's formula: In the eigenbasis of H 0 , the formula reads: if E n = E m . Otherwise, the matrix element is β n| A |m e −βEn . Therefore: Here p n = e −βEn /Z(0) is the population of level E n at equilibrium. Using Eq. (30): Eq. (1) can be written as: for E n = E m . For E n = E m by choosing the eigenbasis |n such that A is diagonal in the eigen-subspaces of H 0 we have 2 : Therefore: Using this expression, the susceptibility of A reduces to: where the last sum runs over all n and m with E n = E m . One can distinguish the Curie and van Vleck terms in the FDT [4,5]. It is interesting to consider the observableÃ = Λ 0 , i.e., the SLD of ρ λ at λ = 0. Since [Ã, Λ 0 ] = 0 and Λ 0 0 = 0, the susceptibility ofÃ is equivalent to the QFI, i.e., we have χ sÃ = Var(Λ 0 ) 0 = F 0 . Therefore,Ã identifies the most sensitive observable to the perturbation (rather than A), and its sensitivity saturates the Cramér-Rao bound, that is given by:

B Kubo relations
Our main results Eqs. (14) and (16) are FDTs for generic quantum Markov systems. One can recover the familiar Kubo quantum FDT for states π λ = e −β(H0−λA) /Z(λ) and Hamiltonian evolution. In this case, the Eq. (16) reads (as before, we denote by φ B (t) the response function of observable B under the perturbation −λ(t)A): To proceed further, we need an additional formula for the SLD. For any real function f (x): Differentiating this equation with respect to λ and setting λ = 0, yields: Particularizing to f (x) = e −βx and using the definition of the SLD given by Eq. (1) the above equation reduces to: Introducing this last result into Eq. (38) we finally reach at: which is the standard Kubo formula. It is interesting to recall that the Kubo formula allows one to express the absorptive part of the susceptibility in terms of a two time correlation in the frequency domain. This is the so-called quantum FDT, which reads [4]: which is equivalent to: Here we see the clear advantage of the introduction of the SLD conjugated variable Λ 0 : it allows us to express the response function in the time domain as a correlation, both in the quantum and in the classical case.

C The response function of coupled harmonic oscillators
Here with the help of the stationary state of the coupled harmonic oscillators, we find the response function of any observable B to the perturbation. The covariance matrix (CM) corresponding to the density matrix of Eq. (23), i.e., for a vanishing interaction, is described by the following diagonal matrix: For a non-zero coupling, the CM is not diagonal anymore, as the coupling establishes correlations amongst the two harmonic oscillators. The corresponding CM is given by the following matrix [37], with ζ = γ 2 +δ 2 4J 2 +γ 2 +δ 2 , D = 2J 2 (N1+N2+1) γ 2 +δ 2 , and C = J(N1−N2) γ 2 +δ 2 . Next, we need to identify (i) the SLD associated to J, and (ii) the time evolution of the desired observable B(t), under the unperturbed map. One can find (i) with the help of Eq. (47). We know that for such a Gaussian state, Λ 0 can be expressed as a linear combination of all the second order moments of the quadratures [35,36], namely with d j s and c j s being coefficients which are to be determined. To this end, we make benefit of Eq. (4) of the main text, which states that for any observable B we have: Next, we imply this relation to all of the quadratic observables appearing in (48), i.e., we choose B ∈ {x 2 1 , x 2 2 , x 1 p 1 + p 1 x 1 , . . . }. With such choices of B, the left hand side of Eq. (49) can be evaluated by taking derivative from the covariance matrix σ ∞ . Moreover, by using the Wick's theorem for Gaussian distributions [49], one can easily simplify the right hand side and write it down in terms of the elements of σ ∞ as well. We shall start by focusing on the local terms. As an example, for B = x 2 1 we have: where we used the fact that σ ∞ J=0 , has no off-diagonal terms, as stated by Eq. (46). For B = p 2 1 one finds a similar equation, with the change of coefficients d 1 ↔ d 2 . This implies that, d 1 = d 2 = 0. In the same manner, one can find that d 4 = d 5 = 0. For the other two local terms i.e., B ∈ { 1 2 (x 1 p 1 + p 1 x 1 ), 1 2 (x 2 p 2 + p 2 x 2 )}, Eq. (49) leads to: whence, d 3 = d 6 = 0 as well. This confirms that the coefficients associated to local observables are zero, i.e., d i = 0 ∀i. However, the non-local coefficients are non-zero. For B ∈ {x 1 x 2 , x 1 p 2 , p 1 x 2 , p 1 p 2 }, using Eq. (49) yields: By replacing Eqs. (58) and (54) into the definition of the response function, Eq. (16), we obtain: φ x1x2 (t) = Corr Λ 0 , x 1 x 2 (t) 0 = e −γt c 1 (σ ∞ 11 ) J=0 (σ ∞ 22 ) J=0 [cos(ω 1 t) cos(ω 2 t) + sin(ω 1 t) sin(ω 2 t)] + c 2 (σ ∞ 11 ) J=0 (σ ∞ 22 ) J=0 [cos(ω 1 t) sin(ω 2 t) − sin(ω 1 t) cos(ω 2 t)] = N 2 − N 1 γ 2 + δ 2 e −γt δ cos δt + γ sin δt where we define θ = arctan γ/δ. By using the same strategy, it is easy to verify that the response function of the other non-local observables read as: