Linear response theory for quantum Gaussian processes

Fluctuation dissipation theorems connect the linear response of a physical system to a perturbation to the steady-state correlation functions. Until now, most of these theorems have been derived for finite-dimensional systems. However, many relevant physical processes are described by systems of infinite dimension in the Gaussian regime. In this work, we find a linear response theory for quantum Gaussian systems subject to time dependent Gaussian channels. In particular, we establish a fluctuation dissipation theorem for the covariance matrix that connects its linear response at any time to the steady state two-time correlations. The theorem covers non-equilibrium scenarios as it does not require the steady state to be at thermal equilibrium. We further show how our results simplify the study of Gaussian systems subject to a time dependent Lindbladian master equation. Finally, we illustrate the usage of our new scheme through some examples. Due to broad generality of the Gaussian formalism, we expect our results to find an application in many physical platforms, such as opto-mechanical systems in the presence of external noise or driven quantum heat devices.


Introduction
Fluctuation dissipation theorems (FDTs) provide very powerful tools to study the linear response of physical systems close to their steady state. The aim of such theorems is to establish and quantify a connection between (i) the linear response of the system under study to a (time-dependent) perturbation, and (ii) the steady-state correlation functions. Different versions of FDT appear, depending on whether the system under study is classical [1][2][3][4][5][6] or quantum [7][8][9][10][11][12][13], or whether the steady state is thermal [7], or a generic non-equilibrium steady state [5,6,8]. Response functions have been used to estimate noise [14,15], to study topological insulators [16] or to witness and quantify non-Markovianity of quantum systems [17]. For a thermal system, or a thermal system subject to a quench, the FDT is connected to the quantum Fisher information [18,19]. On the account of the fact that the quantum Fisher information is a witness of multipartite entanglement [20][21][22], one can benefit its connection to the FDT in order to detect multipartite entanglement close to thermal equilibrium. Some recent works report violation of FDTs under certain circumstances [23,24].
To date, the majority of theoretical works in the quantum domain have been focused on finite dimensional systems and close to thermal equilibrium. Many physical processes of interest, however, are described by continuous variable systems in the Gaussian regime with infinite-dimensional Hilbert space. These systems also find an application for quantum information technologies and, in fact, have been successfully used for quantum teleportation [25], crafting cluster states with enormous number of entangled states [26] or secure quantum key distribution [27]. It is therefore relevant and timely to establish FDT for quantum Gaussian systems. These theorems should be phrased in terms of the natural tools used to describe Gaussian systems, based on first and second moments rather than density matrices.
In this work we address all these issues and provide a linear response theory for Gaussian continuous variable quantum systems. More specifically, we consider processes described by Gaussian quantum channels and derive the linear response of the covariance matrix. The formalism can find an application in many different scenarios, since it covers the case of time-dependent fluctuations and non-equilibrium scenarios, as it does not assume the initial state to be thermal.
The structure of the article is as follows: in section 2 we review quantum Gaussian channels. The aim of this section is to provide the minimal necessary tools for this study; for more about Gaussian quantum channels see [28,29] and the references therein. In section 3 we set our framework and present the main results. We prove these results in section 4. In section 5, we discuss the application of our theorem for those cases in which the channel is described by a Lindbladian master equation. We show how to use our results in section 6 through some examples. Finally, in section 7 we conclude and discuss future directions.

Definition of the Gaussian scenario
We consider N-mode bosonic systems with the quadrature vector R q q p p , ... ; , ...
Here and throughout the article, T stands for transpose, and the elements q i and p i represent the position and momentum of the ith mode, respectively. The quadratures respect the bosonic algebra: In our notation, N 0 is an N×N matrix of zeros, while N  is the identity matrix of size N. In the rest of this work we set 1  = unless otherwise mentioned. We recall that, by definition, Gaussian systems are those with a Gaussian characteristic function [28,29]. In turn, the characteristic function, denoted by c h ( ), reads as: -W being the Weyl operator, and the phase space vector η belongs to N 2  . Therefore, the characteristic function of a Gaussian system has the following shape: Here, we use the displacement vector, denoted by d, and the covariance matrix, denoted by σ, which are given by: , in order to lighten our notation. The covariance matrix σ obeys the uncertainty principle: σ+Ω/20 [28,30].
As already advanced, Gaussian systems are fully described by their first and second moments. Thus, Gaussian channels can be completely identified by their action on the displacement vector and the covariance matrix. Denoting an arbitrary Gaussian quantum channel by M, in the most generic case it operates on the quadratures vector and the covariance matrix as follows [28,29,31]:   Î´are real matrices. The complete positivity of the map dictates that [28]: Hereafter, without loss of generality, we restrict to zero-mean Gaussian states, and focus on Gaussian channels which map zero-mean states to zero-mean states. This is to say: d=f=0. On this account, the map M could be alternatively characterized by the set {X, Y}. In section 5 we review how to bring the particular case of (quadratic) Lindbladian master equations into the standard form of Gaussian channels.

Framework and main results
We work with a one-parameter family of Gaussian quantum channels M l , where λ is a real parameter, and can represent the strength of an external magnetic field or temperature, to name a few. See section 6 for some examples. Let σ λ be a fixed point covariance matrix of the Gaussian channel M l . This implies that: 1 where å stands for summation (not to be mistaken with Σ, the matrix of second order operators), and Φ(t) is the response function which reads: Here, Δ t stands for time differentiation, i.e. Δ t ( f (t))≡f (t+1)−f (t) and X t 0 represents t times the application of X 0 l l= | . If each map is applied for an infinitesimal time δt→0, we have the continuous version of the response function: Thus, in order to find the response function we simply need to find: (i) the static linear response ς, and (ii) the time evolution of ς under the unperturbed channel X 0 . Three further comments/results are in order:

Static linear response
Our result is fully consistent with the static linear response. Consider a scenario in which (i) the perturbation is constant in time, i.e. t l l = ( ) , and (ii) the map has σ λ as its unique fixed point 5 . For any arbitrary time t the linear response simplifies to: In the limit of t  ¥ the upper value vanishes (see section 4.2), that is: We immediately revive the static linear response:

Kubo's response function
Consider a system at thermal equilibrium with the density matrix ρ 0 =exp(−β H 0 )/Z 0 . Here, the quadratic Hamiltonian is given by H R GR T 0 1 2

=
, and β is the inverse temperature. The system is disturbed by adding a perturbative time-dependent term to its Hamiltonian: . Notice that both G and g are symmetric 2N by 2N real matrices. The response function of the covariance matrix reads as: 5 The static linear response holds true only for a scenario in which the steady state is unique. Indeed for the unitary dynamics and thermal states this is not the case, nonetheless our main result-equation (15)-still holds, from which we obtaine the Kubo response. with ). Therefore, the linear response of Hamiltonian drivings can be identified only by having σ 0 , G and g.

Alternative expression for the linear response
Practically, equations (13) and (15) provide a very useful formalization to find the linear response, for them relying on two elements that are easy to calculate, but they do not seem to follow the usual structure of FDT. However, one can write down an alternative FDT that looks more similar to the traditional one, in particular as presented in [32]. This reads as: where the correlation function is evaluated according to the fixed point of the unperturbed map, hence the index '0'. Here the elements of the matrix Σ(t) are the time evolution of the quadratures under 0 In addition, Λ 0 is the symmetric logarithmic derivative (SLD). The SLD is a Hermitian operator with a vanishing expectation value: 0 0 0 áL ñ = -again the index '0' indicates that the trace is evaluated over the fixed point of the unperturbed map. In our case it can be written down as a linear combination of second order quadratures: with the matrix of coefficients C being the solution of the following equation [33,34]:

Proof of main results
Here we present the proof of equations (13), (17), and (19). The proof of equation (20) is presented in the appendix B.

The response function
We start by expanding equation (11) and keeping the terms up to the first order in λ. This yields: To proceed further, we need to identify how M acts on σ 0 . To this aim, we notice that criterion (9) implies that: By substituting in the response function, we have: . In section 2 we explained how the map 0 M applies to covariance matrices, however, ς is not a covariance matrix. Thus, we have to identify how the map acts on the static linear response ς. To this end, we focus on the time evolution of the covariance matrix σ λ : Therefore, we identify the first two terms of the right-hand side as t . Plugging this into(25) completes our proof of equation (13). (17) The map 0 M has a unique fixed point σ 0 . This is to say, for any initial covariance matrix σ we have:

Proof of equation
Particularizing to σ λ =σ 0 +λς, yields: Since the equality holds for any value of the parameter, the term proportional to λ should vanish, which proves equation (17).

FDT for thermal states and Hamiltonian evolutions
, and h R gR Recall that, such unitary dynamics corresponds to a symplectic transformation that operates on the covariance matrix (see appendix A for details). In particular we have By initially preparing the system at the thermal state ρ 0 (corresponding to λ=0) the linear response reads as: is a 2N by 2N matrix that represents the Heisenberg picture evolution of all of the quadratures. From the Heisenberg equation (or by using equation (A.1)) we have: By plugging (30) into (29), and using the cyclic property of the trace we have: To proceed further, we notice that [H 0 −λh, ρ λ ]=0, and hence, by taking the derivative with respect to λ, and evaluating at λ=0, we have: By inserting the above identity in (31), and using again the cyclic property of the trace, we have: which can be rewritten in the shape of the standard Kubo-response function: The Kubo response function (33) can be brought into a more useful shape. To this end, let us look at an ). The response function of this object has two parts, i.e.
where from first to the second equation we use the definition of h, from second to the third we use the fact that for a unitary dynamics R t S R G t = ( ) , from the third to the fourth we expand the commutators and benefit from the canonical commutation relation, and in the last line we reorder everything to show them as product of matrices. By writing the same expression for t R R m l F ( ), and adding it up to the above result, one obtains: ). In a more compact form, for any second order moment we have: The equation (37) can be considered as the Kubo response function for the covariance matrix of a Gaussian system. Since it is purely defined in terms of the original Hamiltonian (G) and the driving force (g), it does not require finding ς.

FDT for Lindbladian master equations
The stationary state-if it exists-and the elements of the Gaussian channel equivalent to the Lindbladian master equation are found routinely. Let us have a quick reminder about how to formalize this-one could also see for instance [35]. Consider the following master equation: Since we are interested in Gaussian dynamics, the Hamiltonian is quadratic in the quadrature operators, and can be written as: while the Lindbladian operators can be written as: being a vector of size 2N. We have ignored some constants in both expressions above, and also a linear dependence of H on the quadratures, since they shall not affect our results significantly. With these definitions, one can write down the master equation for the covariance matrix and the quadratures vector as follows: and, the rectangular matrix C is defined as C c c c ; ; ...; [36][37][38][39] for the derivation.) From here, the application of the map in an infinitesimal time δt can be identified as: Finally, the stationary state covariance matrix can be obtained by letting the left-hand side of (41b) equal to zero. This leads to solving the Lyapunov equation that reads as [40]: ( ) The answer to this equation exists and is unique if all of the eigenvalues of A have negative real parts and is given by:

Examples
We conclude our study by applying our formalism to two physically relevant examples: a driven harmonic oscillator and a cascaded optimal parametric oscillator.

Driven harmonic oscillator
The thermalization and/or dynamics of quantum open systems is sometimes described in the collision model framework (see for instance [41][42][43][44][45]). Following [46] we use the collision model to address the dynamics of a system in presence of thermal noise. The system consists of a single bosonic mode, while the environment consists of an infinite number of bosonic modes. The system mode consecutively interacts for some time t d with individual environmental modes. Let σ E (t) denote the state of the environmental mode that interacts with the system at time t: Here (c(t)−2)/2 represents the mean photon number of the environment mode. Furthermore, we denote the covariance matrix of the system at time t with σ s (t) . After colliding with the t/δtth mode of the environment, the covariance matrix of the system maps to: with the index E meaning that we trace out the environmental mode. Moreover, the symplectic matrix S η depends on the interaction time δt but we have dropped this dependence for lightening our notation. It is given by [46]: Here, the parameter η quantifies the thermalization rate-which again depends on the interaction length δt-in particular for η=0 the system thermalizes after one application of the map (in this case, the system and the environment exchange their states), whereas for η=1 it never thermalizes (in this case, the system and the environment do not interact, and their states remain unchanged). By plugging this symplectic transformation into equation (46), we can describe the dynamics of the system covariance matrix by means of a quantum Gaussian channel-i.e. with the form of equation (6). In particular, one can easily check that , ∀t) brings any initial system covariance matrix to the steady state c s 0 2  s ¥ = ( ) . In fact, even if the parameter η is time dependent, the steady state will remain the same, therefore, we chose it to be constant. Now let us assume that the time dependence of c(t) is a linear correction, i.e, at any time we have c(t)=c 0 +λ(t), with t c 0 l  | ( )| . Moreover, the initial covariance matrix of the system is c 0 . By using our FDT we aim at identifying the response function. First, notice that , which vanishes exponentially in time.
Without any loss of generality let t t cos 0 l l n = ( ) , with ν being the (potentially tunable) modulation frequency. We shall choose t 1 nd  , such that the consecutive interaction with the environment is smooth.
Specifically, it would be interesting to find the amplitude of the response for different ν. To this aim, it is useful to study the linear response of the covariance matrix to the strength of the perturbation: These graphs show how the system responds to a perturbation imposed by manipulating the environment degrees of freedom. Such perturbation can be realized by e.g. changing the temperature or the frequency of the modes in the environment. The case with ν=0, specifically illustrates the relaxation of the system to a new thermal state after a quench.
Indeed, the linear response is initially zero. For small t, it grows with t h , regardless of the modulation frequency ν. As time increases, the last term in equation (48) vanishes exponentially. Thus, at long times the system will be oscillating around its initial state unless for the case with ν=0, i.e. when we have a constant perturbation. The maximum of the response at long times is achieved at t * , the solution of t tan * n n h = ( )˜. The value of linear response at such maximum is 2 2 h n h +˜. Clearly, for any value of h the modulation frequency with the biggest linear response corresponds to ν=0. Finally, if n h  ¥ |˜| , the response goes to zero, hence the noise is canceled out.

Cascaded optical parametric oscillator (OPO)
An OPO coupled to a vacuum field is described by the Hamiltonian H a a i 4 , with 0   denoting the effective pump intensity. Here, we use the standard definition of the annihilation and creation operators, that read as a , respectively. The coupling to the vacuum is described by the Lindbladian operator L a k = , with κ>0 being the damping cavity rate. For this system, it is not difficult to see that the operators G and C read as follow: From here, one can obtain the matrices A and D that appear in equation (43): Thus, the steady state covariance matrix exists if κ>ò, and reads as diag , Again, the criterion for having a steady state is max , [35]. Before trying to solve the Lyapunov equation, we note that the matrices A and G are block diagonal, and can be solved in the corresponding blocks. Therefore, the resulting covariance matrix will be a direct sum of two different terms [35] (one is completely in position subspace, the other one in the momentum subspace, while there are no correlations between the two). In particular, we need to find the exponential of non-Hermitian matrices A x and A p , the position and momentum subspaces of the matrix A, respectively. To this end, we use the Jordan canonical form of the matrices, that is A V J V 1 = a aa a -. After doing some straightforward algebra, one finds: where we define g ± =(ò 1 +ò 2 ± 2κ)(ò 1 ±κ), and h k 2 ) . This state is always entangled [35]. By making use of the above covariance matrix we can find the other necessary element for our FDT, namely 0 V s = ¶ l l l= | . In turn, the response reads as: Notice that, since all the matrices involved in the above equation are block-diagonals, the response will be blockdiagonal as well. Therefore, we deal with response function in the position and momentum blocks separately. For instance, Φ x (t) read as: with σ x being the first matrix in the rhs of equation (55). Setting κ as the driving parameter, that is by choosing t t cos k k l n = + ( ) , the linear response for different values of modulation frequency is depicted in figure 2. Our first observation is that the position quadratures are more responsive to the perturbation, with x 2 2 á ñ having the biggest amplitude of oscillations. From a sensing (estimation) point of view, this means that x 2 2 á ñ is the most sensitive quadrature measurement in estimation of λ (however, one can design measurements which are superposition of different quadratures and perform even better than x 2 2 á ñ). We further notice that the response to a perturbation with bigger ν is smaller-except for the p p 1 2 á ñ. In fact we can prove that, after long enough time, the amplitude of oscillations of the linear response monotonically decreases with ν. To see this, we recall that the amplitude of oscillations at long times is given by the magnitude of the dynamical susceptibility. The latter is defined as the Fourier transform of the response function: where the integration over negative times is ignored because Φ(t<0)=0 (due to causality). On the other hand, for any perturbation of the form t t cos l l n = ( ) , the linear response at long times reads as: Also in the first line, we use the fact that the response function vanishes exponentially at large enough t, so that we can replace the upper bound of the integral with infinity. Thus, the dynamical susceptibility characterizes the amplitude of oscillations at long times. In figure 3 we depict the dynamical susceptibility of position and momentum quadratures. The monotonic decrease in χ x, p (ω) with ω makes it clear why in figure 2 we see the amplitude of oscillations decrease with increasing ω. Thus, in an estimation scenario it is more suitable to modulate the perturbation with a small or vanishing frequency. Whereas, if we treat the perturbation as a noise, we shall modulate it with a higher frequency in order to cancel it ot. Comparing the two panels of figure 3 also clarifies why the position quadratures have a considerably larger response than the momentum quadratures. Here we have set the parameters ò 1 =1, ò 2 =1.1, and the coupling κ is derived with the profile t 1.5 cos k l n = + ( ). Overall, the position quadratures have a bigger response than the momentum ones. In particular, the position of the second oscillator is the most sensitive one. Moreover, it is seen that the bigger frequency ν=0.5 leads to a smaller response. In fact, for all observables except p p 1 2 á ñ, the linear response monotonically decreases with increasing the frequency. See figure 3.

Conclusions
We have derived a linear response theory for the covariance matrix of Gaussian systems subjected to timedependent Gaussian quantum channels. Our method establishes a connection between the linear response to a time dependent perturbation on the one hand, and on the other hand (i) the static linear response of the system and (ii) the building blocks of the Gaussian channel itself. When dealing with thermal states evolving under unitary dynamics, we revive Kubo's linear response theory. We further present an alternative expression for Kubo's response theory, that is more suitable for Gaussian dynamics. We have then showcased how for any arbitrary (Gaussian) Lindbladian master equation, the two ingredients (i) and (ii) can be identified straightforwardly.
Through the examples of thermalization of a harmonic oscillator, and the cascaded parametric oscillator, we have illustrated how to use our formalism. Since Lindbladian master equations appear often in the description of open quantum systems [47][48][49][50], we expect our results to find an application in many different setups. In particular, they can be used to improve our understanding of opto-mechanical systems [51,52], quantum heat devices [53][54][55][56][57][58], or for the study of the open dynamics of quantum systems in the vicinity of non-thermal steady states.
where we use the fact that R GR R GR , Appendix B. Proof of the alternative shape of the response function On the one hand, by expanding ς(t) with the help of equation (22) we have: where we have defined the non-symmetric two-time correlation matrix σ t0 with the elements represents the time evolution of the quadratures vector, under the unperturbed map. On the other hand, with the help of Wick's theorem, one can expand the fourth order moments that appear in the right-hand side of equation (20), in terms of second order moments. Specifically-by breaking the elements of Σ(t) into two parts as in t R t R t R t R t m n m n n m ,