Linear response theory for quantum Gaussian processes

TL;DR

This paper develops a linear response theory for quantum Gaussian systems, establishing a fluctuation dissipation theorem for their covariance matrices that applies to non-equilibrium states and broad physical platforms.

Contribution

It introduces a fluctuation dissipation theorem for quantum Gaussian systems' covariance matrices, extending linear response theory to non-equilibrium and infinite-dimensional regimes.

Findings

Derived a fluctuation dissipation theorem for Gaussian covariance matrices.

Simplified analysis of Gaussian systems under time-dependent Lindbladian dynamics.

Applicable to various physical platforms like opto-mechanical systems and quantum heat devices.

Abstract

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.

Figures

Equations

$$\displaystyle\Omega=i\hbar\left(\begin{array}[]{cc}{\mathbb{0}}_{N}&{\mathbbm{I}}_{N}\\ -{\mathbbm{I}}_{N}&{\mathbb{0}}_{N}\end{array}\right).$$

$$\chi (\eta )$$

$$\chi (\eta )={\mathrm{e}}^{\frac{1}{2}{\eta }^T\Omega \sigma \Omega \eta -{\mathrm{d}}^T\Omega \eta }.$$

$$\mathrm{d}$$

$${\sigma }_{ij}$$

$$\mathcal{M}:\mathrm{d}\mapsto X\mathrm{d}+\mathrm{f},$$

$$\mathcal{M}:\sigma \mapsto X\sigma X^T+Y.$$

$$Y+\frac{\Omega }{2}-\frac{1}{2}X\Omega X^T\ge 0.$$

$${\mathcal{M}}_{\lambda }{\sigma }_{\lambda }={\sigma }_{\lambda },$$

$${\sigma }_{\lambda }=X_{\lambda }{\sigma }_{\lambda }{X}_{\lambda }^T+Y_{\lambda }.$$

$$\sigma (t)={\mathcal{M}}_{\lambda (t)}\bullet {\mathcal{M}}_{\lambda (t-1)}\cdots \bullet {\mathcal{M}}_{\lambda (1)}{\sigma }_0,$$

$$\sigma (t)={\sigma }_0+\sum _{s=1}^t\lambda (t-s)\Phi (s),$$

$$\displaystyle\Phi(t)=-\Delta_{t}~{}\big{(}X_{0}^{t}~{}\varsigma~{}{X_{0}^{T}}^{t}\big{)}.$$

$$\sigma (t)$$

$$\Phi (t)$$

$$\sigma (t)-{\sigma }_0$$

$$\displaystyle=-\lambda~{}\big{(}X_{0}^{s}~{}\varsigma~{}{X_{0}^{T}}^{s}\big{)}_{s=0}^{s=t}+{\cal O}(\lambda^{2}).$$

$$\underset{s\to \infty }{\lim }{X}_0^s\varsigma X_0^{Ts}=0.$$

$${\sigma }_{\lambda }-{\sigma }_0=\lambda \varsigma +\mathcal{O}({\lambda }^2).$$

$$\Phi (t)=\frac{i}{\hslash }{S}_G^t\tilde{\sigma }{S_G^T}^t,$$

$$\Phi (t)=-{\partial }_t\mathrm{Corr}(\Sigma (t),{\Lambda }_0{)}_0,$$

$${\Lambda }_0=\sum _{i,j}{C}_{ij}\left(R_i\circ R_j-{⟨R_i\circ R_j⟩}_0\right),$$

$$\varsigma =4{\sigma }_{0}C{\sigma }_0+\Omega C\Omega .$$

$$\sigma (t)$$

$$={\mathcal{M}}_0^t{\sigma }_0+\sum _{s=0}^t\lambda (s){\mathcal{M}}_0^{t-s}\mathrm{M}{\mathcal{M}}_0^{s-1}{\sigma }_0$$

$$={\sigma }_0+\sum _{s=1}^t\lambda (s){\mathcal{M}}_0^{t-s}\mathrm{M}{\sigma }_0$$

$$={\sigma }_0+\sum _{s=1}^t\lambda (s)\Phi (t-s).$$

$$({\mathcal{M}}_0+\lambda \mathrm{M})({\sigma }_0+\lambda \varsigma )={\sigma }_0+\lambda \varsigma +\mathcal{O}({\lambda }^2)$$

$$\Rightarrow \mathrm{M}{\sigma }_0=(1-{\mathcal{M}}_0)\varsigma .$$

$$\Phi (t)=({\mathcal{M}}_0^t-{\mathcal{M}}_0^{t+1})\varsigma =-{\Delta }_t\varsigma (t),$$

$${\mathcal{M}}_0^t{\sigma }_{\lambda }$$

$$=X_0^t({\sigma }_0+\lambda \varsigma ){X_0^T}^t+Y_0(t)+\mathcal{O}({\lambda }^2)$$

$$=X_0^t{\sigma }_0{X_0^T}^t+Y_0(t)+\lambda X_0^t\varsigma {X_0^T}^t+\mathcal{O}({\lambda }^2).$$

$${\mathcal{M}}_0^t{\sigma }_0+\lambda {\mathcal{M}}_0^t\varsigma =X_0^t{\sigma }_0{X_0^T}^t+Y_0(t)+\lambda X_0^t\varsigma {X_0^T}^t+\mathcal{O}({\lambda }^2).$$

$${\sigma }_0=\underset{t\to \infty }{\lim }{\mathcal{M}}_0^t\sigma =\underset{t\to \infty }{\lim }\left(X_0^t\sigma {X_0^T}^t+{\int }_0^{t}ds{X}_0^s{Y}_0{X_0^T}^s\right).$$

Full text

1 Introduction

Fluctuation dissipation theorems (FDT) 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, 8, 6]. 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.

2 Definition of the Gaussian scenario

We consider $N$-mode Bosonic systems with the quadrature vector $R\equiv(q_{1},~{}\dots q_{N};~{}p_{1},~{}\dots p_{N})^{T}$. Here and throughout the article, $T$ stands for transpose, and the elements $q_{i}$ and $p_{i}$ represent the position and momentum of the $i$th mode, respectively. The quadratures respect the Bosonic algebra: $[R_{i},R_{j}]=\Omega_{i,j}$, where $\Omega$ is the symplectic matrix

[TABLE]

In our notation, ${\mathbb{0}}_{N}$ is an $N\times N$ matrix of zeros, while ${\mathbbm{I}}_{N}$ is the identity matrix of size $N$. In the rest of this work we set $\hbar=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 $\chi(\eta)$, reads as:

[TABLE]

with $\rho$ being the density matrix, $W_{\eta}\equiv{\rm e}^{-\eta^{T}~{}\Omega~{}R}$ being the Weyl operator, and the phase space vector $\eta$ belongs to ${\mathbbm{R}}^{2N}$. Therefore, the characteristic function of a Gaussian system has the following shape:

[TABLE]

Here, we use the displacement vector, denoted by ${\rm d}$, and the covariance matrix, denoted by ${\sigma}$, which are given by:

[TABLE]

where we define $\Sigma_{ij}\equiv R_{i}\circ R_{j}=\frac{1}{2}(\bar{R_{i}}~{}\bar{R_{j}}+\bar{R_{j}}~{}\bar{R_{i}})$, and $\bar{R_{i}}=R_{i}-{\rm d}_{i}$, in order to lighten our notation. The covariance matrix $\sigma$ obeys the uncertainty principle: $\sigma+\Omega/2\geq 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 ${\cal M}$, in the most generic case it operates on the quadratures vector and the covariance matrix as follows [28, 29, 31]:

[TABLE]

Here, ${\rm f}\in{\mathbbm{R}}^{2N}$, while $X$ and $Y\in{\mathbbm{R}}^{2N}\times{\mathbbm{R}}^{2N}$ are real matrices. The complete positivity of the map dictates that [28]:

[TABLE]

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: ${\rm d}={\rm f}=0$. On this account, the map ${\cal 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.

3 Framework and main results

We work with a one-parameter family of Gaussian quantum channels ${\cal M}_{\lambda}$, where $\lambda$ 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 $\sigma_{\lambda}$ be a fixed point covariance matrix of the Gaussian channel ${\cal M}_{\lambda}$. This implies that:

[TABLE]

or alternatively

[TABLE]

Note that we allow both $X_{\lambda}$ and $Y_{\lambda}$ depend on the parameter $\lambda$. Since we are interested in the linear response, we work in a regime where the parameter $\lambda$ can be considered as a linear contribution to the channel and its fixed point covariance matrix. Therefore, we assume that ${\cal M}_{\lambda}={\cal M}_{0}+\lambda{\rm M}+{\cal O}(\lambda^{2})$, and $\sigma_{\lambda}=\sigma_{0}+\lambda\varsigma+{\cal O}(\lambda^{2})$, and we safely ignore the second and higher orders. In other words, if we normalize ${\cal M}_{0}$ and ${\rm M}$ (and doing the same for $\sigma_{0}$ and $\varsigma$) such that they have the same operator norm, then $|\lambda|\ll 1$. Notice that, neither ${\rm M}$ is a Gaussian quantum channel on its own, nor is $\varsigma$ a covariance matrix. We refer to $\varsigma$ as the static linear response of the covariance matrix, and we have: $\varsigma\equiv\partial_{\lambda}\sigma_{\lambda}|_{\lambda=0}$. In a Markovian scenario, the time evolution of the covariance matrix is described by the consecutive operations of the map. Choosing an initial covariance matrix $\sigma_{0}$, that is the fixed point of ${\cal M}_{0}$, at discrete time steps $t=1,2,3,\dots$ we have:

[TABLE]

where we allow for a time dependent parameter $\lambda(t)$. The aim is to characterize the linear response of $\sigma(t)$ in terms of steady state correlations, that is elements of steady state covariance matrix.

Our main result expresses the linear response of the covariance matrix as follows:

[TABLE]

where $\sum$ stands for summation (not to be mistaken with $\Sigma$, the matrix of second order operators), and $\Phi(t)$ is the response function which reads:

[TABLE]

Here, $\Delta_{t}$ stands for time differentiation, i.e., $\Delta_{t}(f(t))\equiv f(t+1)-f(t)$ and $X_{0}^{t}$ represents $t$ times the application of $X_{\lambda}|_{\lambda=0}$. If each map is applied for an infinitesimal time $\delta t\to 0$, we have the continuous version of the response function:

[TABLE]

Thus, in order to find the response function we simply need to find: (i) the static linear response $\varsigma$, and (ii) the time evolution of $\varsigma$ under the unperturbed channel $X_{0}$. Three further comments/results are in order:

3.1 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., $\lambda(t)=\lambda$, and (ii) the map has $\sigma_{\lambda}$ as its unique fixed point 111The 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—Eq. (17)—still holds, from which we obtaine the Kubo response.. For any arbitrary time $t$ the linear response simplifies to:

[TABLE]

In the limit of $t\to\infty$ the upper value vanishes [see section 4.2], that is:

[TABLE]

We immediately revive the static linear response:

[TABLE]

3.2 Kubo’s response function

Consider a system at thermal equilibrium with the density matrix $\rho_{0}={\rm exp}(-\beta H_{0})/Z_{0}$. Here, the quadratic Hamiltonian is given by $H_{0}=\frac{1}{2}R^{T}G~{}R$, and $\beta$ is the inverse temperature. The system is disturbed by adding a perturbative time-dependent term to its Hamiltonian: $H(t)=H_{0}-\lambda(t)h$ with $h=\frac{1}{2}R^{T}g~{}R$. Notice that both $G$ and $g$ are symmetric $2N$ by $2N$ real matrices. The response function of the covariance matrix reads as:

[TABLE]

with ${\tilde{\sigma}}\coloneqq(\Omega g\sigma_{0}-\sigma_{0}g\Omega)$ and $S_{G}^{~{}t}\coloneqq{\rm exp}(-it\Omega G)$. Therefore, the linear response of Hamiltonian drivings can be identified only by having $\sigma_{0}$, $G$ and $g$.

3.3 Alternative expression for the linear response

Practically, Eqs. (15) and (17) 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:

[TABLE]

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 $\Sigma(t)$ are the time evolution of the quadratures under ${\cal M}_{0}$, such that $\Sigma_{m,n}(t)=R_{m}(t){\circ}R_{n}(t)$. In addition, $\Lambda_{0}$ is the symmetric logarithmic derivative (SLD). The SLD is a Hermitian operator with a vanishing expectation value: $\left<\Lambda_{0}\right>_{0}=0$—again the index “[math]” 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:

[TABLE]

with the matrix of coefficients $C$ being the solution of the following equation [33, 34]:

[TABLE]

4 Proof of main results

Here we present the proof of Eqs. (15), (19), and (21). The proof of Eq. (22) is presented in the B.

4.1 The response function

We start by expanding Eq. (13) and keeping the terms up to the first order in $\lambda$. This yields:

[TABLE]

In the last line we define the response function $\Phi(t)={\cal M}_{0}^{t}{\rm M}\sigma_{0}$. To proceed further, we need to identify how ${\rm M}$ acts on $\sigma_{0}$. To this aim, we notice that criterion (11) implies that:

[TABLE]

By substituting in the response function, we have:

[TABLE]

where we define $\varsigma(t)\equiv{\cal M}_{0}^{t}\varsigma$. In section 2 we explained how the map ${\cal M}_{0}$ applies to covariance matrices, however, $\varsigma$ is not a covariance matrix. Thus, we have to identify how the map acts on the static linear response $\varsigma$. To this end, we focus on the time evolution of the covariance matrix $\sigma_{\lambda}$:

[TABLE]

Where we define $Y_{0}(t)\equiv\sum_{s=0}^{t-1}X_{0}^{s}~{}Y_{0}~{}{X_{0}^{T}}^{s}$. By substituting $\sigma_{\lambda}=\sigma_{0}+\lambda\varsigma+{\cal O}(\lambda^{2})$ on the left hand side of (4.1), we have:

[TABLE]

Therefore, we identify the first two terms of the right hand side as ${\cal M}_{0}^{t}\sigma_{0}$, and $\varsigma(t)=X_{0}^{t}\varsigma{X_{0}^{T}}^{t}$. Plugging this into (4.1) completes our proof of Eq. (15).

4.2 Proof of Equation (19)

The map ${\cal M}_{0}$ has a unique fixed point $\sigma_{0}$. This is to say, for any initial covariance matrix $\sigma$ we have:

[TABLE]

Particularizing to $\sigma_{\lambda}=\sigma_{0}+\lambda\varsigma$, yields:

[TABLE]

Since the equality holds for any value of the parameter, the term proportional to $\lambda$ should vanish, which proves Eq. (19).

4.3 FDT for thermal states and Hamiltonian evolutions (Kubo’s response function)

Consider a Gaussian system with $H(t)=H_{0}-\lambda(t)h$, where $H_{0}=\frac{1}{2}~{}R^{T}G~{}R$, and $h=\frac{1}{2}~{}R^{T}g~{}R$. The thermal state $\rho_{\lambda}=\exp(-\beta H_{\lambda})/{\rm Z_{\lambda}}$—being $Z_{\lambda}$ the partition function—is a fixed point of the unitary dynamics produced by $H_{\lambda}\equiv H_{0}-\lambda h$. Recall that, such unitary dynamics corresponds to a symplectic transformation that operates on the covariance matrix (see A for details). In particular we have $X_{0}={\exp}(-i\Omega G)\eqqcolon S_{G}$. By initially preparing the system at the thermal state $\rho_{0}$ (corresponding to $\lambda=0$) the linear response reads as:

[TABLE]

where $\Sigma(t)=S_{G}^{~{}t}~{}\Sigma~{}{S_{G}^{T}}^{t}$ is a $2N$ by $2N$ matrix that represents the Heisenberg picture evolution of all of the quadratures. From the Heisenberg equation (or by using Eq. (90)) we have:

[TABLE]

By plugging (32) into (4.3), and using the cyclic property of the trace we have:

[TABLE]

To proceed further, we notice that $[H_{0}-\lambda h~{},~{}\rho_{\lambda}]=0$, and hence, by taking the derivative with respect to $\lambda$, and evaluating at $\lambda=0$, we have:

[TABLE]

By inserting the above identity in (33), and using again the cyclic property of the trace, we have:

[TABLE]

which can be rewritten in the shape of the standard Kubo-response function:

[TABLE]

The Kubo response function (35) can be brought into a more useful shape. To this end, let us look at an individual element of $\Sigma$, say $\Sigma_{lm}=\frac{1}{2}(R_{l}R_{m}+R_{m}R_{l})$. The response function of this object has two parts, i.e., $\Phi_{\Sigma_{lm}}(t)=\frac{1}{2}~{}(\Phi_{R_{m}R_{l}}(t)+\Phi_{R_{l}R_{m}}(t))$. We have:

[TABLE]

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_{G}^{t}~{}R$, 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 $\Phi_{R_{m}R_{l}}(t)$, and adding it up to the above result, one obtains:

[TABLE]

with ${\tilde{\sigma}}\coloneqq(\Omega g\sigma_{0}-\sigma_{0}g\Omega)$. In a more compact form, for any second order moment we have:

[TABLE]

The equation (39) 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 $\varsigma$.

5 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:

[TABLE]

Since we are interested in Gaussian dynamics, the Hamiltonian is quadratic in the quadrature operators, and can be written as:

[TABLE]

while the Lindbladian operators can be written as:

[TABLE]

with $c_{m}\in\mathbbm{C}^{2N}$ 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:

[TABLE]

with the matrix $A=-i\Omega(G-{\rm Im}(CC^{\dagger}))$ and $D=\Omega{\rm Re}(CC^{\dagger})\Omega$ and, the rectangular matrix $C$ is defined as $C=(c_{1}^{T};c_{2}^{T};\dots;c_{m}^{T})^{T}\in{\mathbbm{C}}^{2N\times m}$ (See the C or [36, 37, 38] for the derivation.) From here, the application of the map in an infinitesimal time $\delta t$ can be identified as:

[TABLE]

Finally, the stationary state covariance matrix can be obtained by letting the left hand side of (43b) equal to zero. This leads to solving the Lyapunov equation that reads as [39]:

[TABLE]

The answer to this equation exists and is unique if all of the eigenvalues of $A$ have negative real parts and is given by:

[TABLE]

6 Examples

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

6.1 Driven harmonic oscillator

The thermalization and/or dynamics of quantum open systems is sometimes described in the collision model framework (see for instance [40, 41, 42, 43, 44]). Following [45] 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 $\delta t$ with individual environmental modes. Let $\sigma_{E}(t)$ denote the state of the environmental mode that interacts with the system at time $t$:

[TABLE]

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 $\sigma_{s}(t)$. After colliding with the $t/\delta t$th mode of the environment, the covariance matrix of the system maps to:

[TABLE]

with the index $E$ meaning that we trace out the environmental mode. Moreover, the symplectic matrix $S_{\eta}$ depends on the interaction time $\delta t$ but we have dropped this dependence for lightening our notation. It is given by [45]:

[TABLE]

Here, the parameter $\eta$ quantifies the thermalization rate—which again depends on the interaction length $\delta t$—in particular for $\eta=0$ the system thermalizes after one application of the map (in this case, the system and the environment exchange their states), whereas for $\eta=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 Eq. (48), we can describe the dynamics of the system covariance matrix by means of a quantum Gaussian channel—i.e., with the form of Eq. (8). In particular, one can easily check that $X=\sqrt{\eta}~{}{\mathbbm{I}}_{2}$, and $Y=(1-\eta)c(t){\mathbbm{I}}_{2}$. We notice that a consecutive application of this map, with fixed $c$ (i.e., $\sigma_{E}(t)=c_{0}\mathbbm{I}_{2}$, $\forall t$) brings any initial system covariance matrix to the steady state $\sigma_{s}(\infty)=c_{0}{\mathbbm{I}}_{2}$. In fact, even if the parameter $\eta$ 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}+\lambda(t)$, with $|\lambda(t)|\ll c_{0}$. Moreover, the initial covariance matrix of the system is $\sigma_{s}(0)=c_{0}\mathbbm{I}_{2}$. By using our FDT we aim at identifying the response function. First, notice that $\varsigma=\partial_{\lambda}\sigma_{\lambda}|_{\lambda=0}=\mathbbm{I}_{2}$. In addition, since $X=\sqrt{\eta}~{}\mathbbm{I}_{2}$, we have $\varsigma(t)=X^{t/\delta t}\varsigma{X^{T}}^{t/\delta t}~{}{\mathbbm{I}}_{2}=\eta^{t/\delta t}~{}{\mathbbm{I}}_{2}$. Therefore, we have $\Phi(t)=-\partial_{t}\eta^{t/\delta t}~{}{\mathbbm{I}}_{2}=-\log(\eta^{1/\delta t})\eta^{t/\delta t}~{}{\mathbbm{I}}_{2}$, which vanishes exponentially in time.

Without any loss of generality let $\lambda(t)=\lambda_{0}\cos\nu t$, with $\nu$ being the [potentially tunable] modulation frequency. We shall choose $\nu\delta t\ll 1$, such that the consecutive interaction with the environment is smooth. Specifically, it would be interesting to find the amplitude of the response for different $\nu$. To this aim, it is useful to study the linear response of the covariance matrix to the strength of the perturbation:

[TABLE]

where we define $\tilde{\eta}=-\log(\eta^{1/\delta t})$ to lighten our notation. In Fig. 1 we depict the linear response, i.e., $\partial_{\lambda}\sigma_{s}(t)|_{\lambda=0}$ versus time, for different frequencies. 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 $\nu=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 $\tilde{\eta}t$, regardless of the modulation frequency $\nu$. As time increases, the last term in Eq. (6.1) vanishes exponentially. Thus, at long times the system will be oscillating around its initial state unless for the case with $\nu=0$, i.e., when we have a constant perturbation. The maximum of the response at long times is achieved at $t^{*}$, the solution of $\tan(\nu t^{*})=\nu/\tilde{\eta}$. The value of linear response at such maximum is $\tilde{\eta}/\sqrt{\nu^{2}+\tilde{\eta}^{2}}$. Clearly, for any value of $\tilde{\eta}$ the modulation frequency with the biggest linear response corresponds to $\nu=0$. Finally, if $\left|\nu/\tilde{\eta}\right|\to\infty$, the response goes to zero, hence the noise is canceled out.

6.2 Cascaded optical parametric oscillator

An optical parametric oscillator (OPO) coupled to a vacuum field is described by the Hamiltonian $H=i\epsilon(a^{\dagger 2}-a^{2})/4$, with $\epsilon\geq 0$ denoting the effective pump intensity. Here, we use the standard definition of the annihilation and creation operators, that read as $a=(x+i~{}p)/\sqrt{2}$ and $a^{\dagger}=(x-i~{}p)/\sqrt{2}$, respectively. The coupling to the vacuum is described by the Lindbladian operator $L=\sqrt{\kappa}~{}a$, with $\kappa>0$ being the damping cavity rate. For this system, it is not difficult to see that the operators $G$ and $C$ read as follow:

[TABLE]

From here, one can obtain the matrices $A$ and $D$ that appear in Eq. (45):

[TABLE]

Thus, the steady state covariance matrix exists if $\kappa>\epsilon$, and reads as $\sigma=\frac{1}{2}\text{diag}(\kappa/(\kappa-\epsilon)~{},~{}\kappa/(\kappa+\epsilon))$.

The cascaded OPO consists of two interacting optical oscillators with local Lindbladian operators. The Hamiltonian reads as $H=H_{1}+H_{2}+i(L_{1}^{\dagger}L_{2}~{}-~{}L_{1}L_{2}^{\dagger})/2$, with $H_{j}=i\epsilon_{j}(a_{j}^{\dagger 2}-a_{j}^{2})/4$, and $L_{j}=\sqrt{\hbar\kappa}a_{j}$. In turn the dissipation is given by a single operator $L=L_{1}+L_{2}$  For convenience, in what follows we work in a representation, where the quadrature vector is represented by $R=(x_{1}~{}~{}\dots~{}~{}x_{N}~{}~{}p_{1}~{}~{}\dots~{}~{}p_{N})^{T}$. The matrices $G$ and $C$ read as:

[TABLE]

Therefore, one can find the matrices $A$ and $D$ as follow:

[TABLE]

Again, the criterion for having a steady state is $\kappa>\max\{\epsilon_{1}~{},~{}\epsilon_{2}\}$ [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_{\alpha}=V_{\alpha}J_{\alpha}V^{-1}_{\alpha}$. After doing some straightforward algebra, one finds:

[TABLE]

and

[TABLE]

Thus we have the dynamics element for our FDT. Moreover, by putting in the expression of the CM, one finds [35]:

[TABLE]

where we define $g_{\pm}=(\epsilon_{1}~{}+~{}\epsilon_{2}~{}\pm~{}2\kappa)(\epsilon_{1}~{}\pm\kappa)$, and $h_{\pm}=(\epsilon_{1}^{2}~{}+~{}\epsilon_{1}~{}\epsilon_{2}~{}\pm~{}\epsilon_{1}~{}\kappa~{}+~{}2\kappa^{2}~{}\mp\kappa~{}\epsilon_{2})(\epsilon_{2}~{}\mp~{}k)^{-1}$. 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 $\varsigma=\partial_{\lambda}\sigma_{\lambda}|_{\lambda=0}$. In turn, the response reads as:

[TABLE]

Notice that, since all the matrices involved in the above equation are block-diagonals, the response will be block-diagonal as well. Therefore, we deal with response function in the position and momentum blocks separately. For instance, $\Phi_{x}(t)$ read as:

[TABLE]

with $\sigma_{x}$ being the first matrix in the rhs of Eq. (85). Setting $\kappa$ as the driving parameter, that is by choosing $\kappa(t)=\kappa+\lambda\cos\nu t$, the linear response for different values of modulation frequency is depicted in Fig. 2. Our first observation is that the position quadratures are more responsive to the perturbation, with $\left<x_{2}^{2}\right>$ having the biggest amplitude of oscillations. From a sensing (estimation) point of view, this means that $\left<x_{2}^{2}\right>$ is the most sensitive quadrature measurement in estimation of $\lambda$ (however, one can design measurements which are superposition of different quadratures and perform even better than $\left<x_{2}^{2}\right>$). We further notice that the response to a perturbation with bigger $\nu$ is smaller—except for the $\left<p_{1}p_{2}\right>$. In fact we can prove that, after long enough time, the amplitude of oscillations of the linear response monotonically decreases with $\nu$.

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:

[TABLE]

where the integration over negative times is ignored because $\Phi(t<0)=0$ (due to causality). On the other hand, for any perturbation of the form $\lambda(t)=\lambda\cos\nu t$, the linear response at long times reads as:

[TABLE]

where we define $\alpha=\cot^{-1}(\frac{{\rm Im}\chi(\omega)}{{\rm Re}\chi(\omega)})$. 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 Fig. 3 we depict the dynamical susceptibility of position and momentum quadratures. The monotonic decrease in $\chi_{x,p}(\omega)$ with $\omega$ makes it clear why in Fig. 2 we see the amplitude of oscillations decrease with increasing $\omega$. 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 Fig. 3 also clarifies why the position quadratures have a considerably larger response than the momentum quadratures.

7 Conclusions

We have derived a linear response theory for the covariance matrix of Gaussian systems subjected to time-dependent 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 [46, 47, 48, 49], 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 [50, 51], quantum heat devices [52, 53, 54, 55, 56, 57], or for the study of the open dynamics of quantum systems in the vicinity of non-thermal steady states.

acknowledgments

We thank Anna Sanpera, Andreu Riera-Campeny, and Janek Kolodynski for fruitful discussions. Support from the Spanish MINECO (QIBEQI FIS2016-80773-P, ConTrAct FIS2017-83709- R, and Severo Ochoa SEV-2015-0522), the ERC CoG QITBOX, the AXA Chair in Quantum Information Science, Fundacio Privada Cellex, and the Generalitat de Catalunya (CERCA Program and SGR1381) is acknowledged. ————————-