Linear response theory for quantum Gaussian processes

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 σ_λ as its unique fixed point⁵ . For any arbitrary time t the linear response simplifies to:

$$\begin{eqnarray}&&\sigma (t)-{\sigma }_{0}=-\lambda {\int }_{0}^{t}{\rm{d}}{s}\,{\partial }_{s}\left({X}_{0}^{s}\,\varsigma \,{{X}_{0}^{{T}^{s}}}^{}\right)+{\mathscr{O}}({\lambda }^{2})=-\lambda {\left({X}_{0}^{s}\,\varsigma \,{{X}_{0}^{T}}^{s}\right)}_{s=0}^{s=t}+{\mathscr{O}}({\lambda }^{2}).\end{eqnarray} \tag{ 16 }$$

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

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{s\to \infty }{X}_{0}^{s}\,\varsigma \ {X}_{0}^{{Ts}}=0.\end{eqnarray} \tag{ 17 }$$

We immediately revive the static linear response:

$$\begin{eqnarray}&&{\sigma }_{\lambda }-{\sigma }_{0}=\lambda \varsigma +{\mathscr{O}}({\lambda }^{2}).\end{eqnarray} \tag{ 18 }$$

Consider a system at thermal equilibrium with the density matrix ρ₀ = exp(−β H₀)/Z₀. Here, the quadratic Hamiltonian is given by ${H}_{0}=\tfrac{1}{2}{R}^{T}GR$, and β is the inverse temperature. The system is disturbed by adding a perturbative time-dependent term to its Hamiltonian: H(t) = H₀ − λ(t) h with $h=\tfrac{1}{2}{R}^{T}gR$. Notice that both G and g are symmetric 2N by 2N real matrices. The response function of the covariance matrix reads as:

$$\begin{eqnarray}&&{\rm{\Phi }}(t)=\displaystyle \frac{{\rm{i}}}{{\hslash }}{S}_{G}^{t}\tilde{\sigma }{{S}_{G}^{{T}^{t}}}^{},\end{eqnarray} \tag{ 19 }$$

with $\tilde{\sigma }:= ({\rm{\Omega }}g{\sigma }_{0}-{\sigma }_{0}g{\rm{\Omega }})$ and ${S}_{G}^{t}:= \exp (-{\rm{i}}{t}{\rm{\Omega }}G)$. Therefore, the linear response of Hamiltonian drivings can be identified only by having σ₀, G and g.

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:

$$\begin{eqnarray}&&{\rm{\Phi }}(t)=-{\partial }_{t}\mathrm{Corr}{\left({\rm{\Sigma }}(t),{{\rm{\Lambda }}}_{0}\right)}_{0},\end{eqnarray} \tag{ 20 }$$

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 ${{\mathscr{M}}}_{0}$, such that ${{\rm{\Sigma }}}_{m,n}(t)={R}_{m}(t)\circ {R}_{n}(t)$. In addition, Λ₀ is the symmetric logarithmic derivative (SLD). The SLD is a Hermitian operator with a vanishing expectation value: ${\left\langle {{\rm{\Lambda }}}_{0}\right\rangle }_{0}=0$—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:

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{0}=\displaystyle \sum _{i,j}{C}_{{ij}}\left({R}_{i}\circ {R}_{j}-{\left\langle {R}_{i}\circ {R}_{j}\right\rangle }_{0}\right),\end{eqnarray} \tag{ 21 }$$

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

$$\begin{eqnarray}&&\varsigma =4{\sigma }_{0}C{\sigma }_{0}+{\rm{\Omega }}C{\rm{\Omega }}.\end{eqnarray} \tag{ 22 }$$

We start by expanding equation (11) and keeping the terms up to the first order in λ. This yields:

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}\sigma (t) & = & ({{\mathscr{M}}}_{0}+\lambda (t)M)\,\cdot \cdot \cdot \,\bullet \,({{\mathscr{M}}}_{0}+\lambda (1)M){\sigma }_{0}={{\mathscr{M}}}_{0}^{t}{\sigma }_{0}+\displaystyle \sum _{s=0}^{t}\lambda (s){{\mathscr{M}}}_{0}^{t-s}M{{\mathscr{M}}}_{0}^{s-1}{\sigma }_{0}\\ & = & {\sigma }_{0}+\displaystyle \sum _{s=1}^{t}\lambda (s){{\mathscr{M}}}_{0}^{t-s}M{\sigma }_{0}={\sigma }_{0}+\displaystyle \sum _{s=1}^{t}\lambda (s){\rm{\Phi }}(t-s).\end{array}\end{array}\end{eqnarray} \tag{ 23 }$$

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

$$\begin{eqnarray}&&({{\mathscr{M}}}_{0}+\lambda {M})({\sigma }_{0}+\lambda \varsigma )={\sigma }_{0}+\lambda \varsigma +{\mathscr{O}}({\lambda }^{2})\Rightarrow {M}{\sigma }_{0}=(1-{{\mathscr{M}}}_{0})\varsigma .\end{eqnarray} \tag{ 24 }$$

By substituting in the response function, we have:

$$\begin{eqnarray*}&&{\rm{\Phi }}(t)=({{\mathscr{M}}}_{0}^{t}-{{\mathscr{M}}}_{0}^{t+1})\varsigma =-{{\rm{\Delta }}}_{t}\varsigma (t),\end{eqnarray*}$$

where we define $\varsigma (t)\equiv {{\mathscr{M}}}_{0}^{t}\varsigma$. In section 2 we explained how the map ${{\mathscr{M}}}_{0}$ 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 σ_λ:

$$\begin{eqnarray}\begin{array}{rcl}{{\mathscr{M}}}_{0}^{t}{\sigma }_{\lambda } & = & {X}_{0}^{t}\,{\sigma }_{\lambda }\,{{X}_{0}^{{T}^{t}}}^{}+\displaystyle \sum _{s=0}^{t-1}{X}_{0}^{s}\,{Y}_{0}\,{{X}_{0}^{T}}^{s}={X}_{0}^{t}({\sigma }_{0}+\lambda \varsigma ){{X}_{0}^{{T}^{t}}}^{}+{Y}_{0}(t)+{\mathscr{O}}({\lambda }^{2})\\ & = & {X}_{0}^{t}{\sigma }_{0}{{X}_{0}^{{T}^{t}}}^{}+{Y}_{0}(t)+\lambda {X}_{0}^{t}\varsigma {{X}_{0}^{{T}^{t}}}^{}+{\mathscr{O}}({\lambda }^{2}),\end{array}\end{eqnarray} \tag{ 25 }$$

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 +{\mathscr{O}}({\lambda }^{2})$ on the left-hand side of (25), we have:

$$\begin{eqnarray}&&{{\mathscr{M}}}_{0}^{t}{\sigma }_{0}+\lambda {{\mathscr{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}}}^{}+{\mathscr{O}}({\lambda }^{2}).\end{eqnarray} \tag{ 26 }$$

Therefore, we identify the first two terms of the right-hand side as ${{\mathscr{M}}}_{0}^{t}{\sigma }_{0}$, and $\varsigma (t)={X}_{0}^{t}\varsigma {{X}_{0}^{{T}^{t}}}^{}$. Plugging this into (25) completes our proof of equation (13).

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

$$\begin{eqnarray}&&{\sigma }_{0}=\mathop{\mathrm{lim}}\limits_{t\to \infty }{{\mathscr{M}}}_{0}^{t}\sigma =\mathop{\mathrm{lim}}\limits_{t\to \infty }\left({X}_{0}^{t}\sigma {{X}_{0}^{{T}^{t}}}^{}+{\int }_{0}^{t}{\rm{d}}{s}\,{X}_{0}^{s}\,{Y}_{0}\,{{X}_{0}^{{T}^{s}}}^{}\right).\end{eqnarray} \tag{ 27 }$$

Particularizing to σ_λ = σ₀ + λς, yields:

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{0} & = & \mathop{\mathrm{lim}}\limits_{t\to \infty }\left({X}_{0}^{t}{\sigma }_{\lambda }{{X}_{0}^{{T}^{t}}}^{}+{\displaystyle \int }_{0}^{t}{\rm{d}}{s}\,{X}_{0}^{s}\,{Y}_{0}\,{{X}_{0}^{{T}^{s}}}^{}\right)=\mathop{\mathrm{lim}}\limits_{t\to \infty }\left({X}_{0}^{t}({\sigma }_{0}+\lambda {\varsigma }_{\lambda }){{X}_{0}^{{T}^{t}}}^{}+{\displaystyle \int }_{0}^{t}{{\rm{d}}{s}{X}}_{0}^{s}\,{Y}_{0}\,{{X}_{0}^{{T}^{s}}}^{}\right)\\ & = & \lambda \,\mathop{\mathrm{lim}}\limits_{t\to \infty }{X}_{0}^{t}\varsigma \ {{X}_{0}^{{T}^{t}}}^{}+\mathop{\mathrm{lim}}\limits_{t\to \infty }\left({X}_{0}^{t}{\sigma }_{0}{{X}_{0}^{{T}^{t}}}^{}+{\displaystyle \int }_{0}^{t}{\rm{d}}{s}\,{X}_{0}^{s}\,{Y}_{0}\,{{X}_{0}^{{T}^{s}}}^{}\right).\end{array}\end{eqnarray} \tag{ 28 }$$

Since the equality holds for any value of the parameter, the term proportional to λ should vanish, which proves equation (17).

Consider a Gaussian system with H(t) = H₀ − λ(t) h, where ${H}_{0}=\tfrac{1}{2}{R}^{T}GR$, and $h=\tfrac{1}{2}{R}^{T}gR$. The thermal state ${\rho }_{\lambda }=\exp (-\beta {H}_{\lambda })/{{Z}}_{\lambda }$—being Z_λ the partition function—is a fixed point of the unitary dynamics produced by H_λ ≡ H₀ − λ h. 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 ${X}_{0}=\exp (-{\rm{i}}{\rm{\Omega }}G)=: {S}_{G}$. By initially preparing the system at the thermal state ρ₀ (corresponding to λ = 0) the linear response reads as:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Phi }}(t) & = & -{\partial }_{t}\left({S}_{G}^{t}\,\varsigma \,{{S}_{G}^{{T}^{t}}}^{}\right)=-{\partial }_{t}\left({S}_{G}^{t}\,{\partial }_{\lambda }\mathrm{tr}{\left[{\rm{\Sigma }}{\rho }_{\lambda }\right]}_{\lambda =0}\,{{S}_{G}^{{T}^{t}}}^{}\right)\\ & = & -{\partial }_{t}\left({S}_{G}^{t}\,\mathrm{tr}\left[{\rm{\Sigma }}\,{\left({\partial }_{\lambda }{\rho }_{\lambda }\right)}_{\lambda =0}\right]{{S}_{G}^{{T}^{t}}}^{}\right)=-\mathrm{tr}\left[{\partial }_{t}{\rm{\Sigma }}(t)\,{\left({\partial }_{\lambda }{\rho }_{\lambda }\right)}_{\lambda =0}\right],\end{array}\end{eqnarray} \tag{ 29 }$$

where ${\rm{\Sigma }}(t)={S}_{G}^{t}\,{\rm{\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 equation (A.1)) we have:

$$\begin{eqnarray}&&{\partial }_{t}{\rm{\Sigma }}(t)=-\displaystyle \frac{{\rm{i}}}{2{\hslash }}[{R}^{T}G\,R,{\rm{\Sigma }}(t)].\end{eqnarray} \tag{ 30 }$$

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

$$\begin{eqnarray}&&{\rm{\Phi }}(t)=-\displaystyle \frac{{\rm{i}}}{2{\hslash }}\mathrm{tr}\left[\left[{\partial }_{\lambda }{\rho }_{\lambda }{| }_{\lambda =0},{R}^{T}G\,R\right]\,{\rm{\Sigma }}(t)\right].\end{eqnarray} \tag{ 31 }$$

To proceed further, we notice that [H₀ − λh, ρ_λ] = 0, and hence, by taking the derivative with respect to λ, and evaluating at λ = 0, we have:

$$\begin{eqnarray}&&[{\rho }_{0},{R}^{T}g\,R]=\left[{\partial }_{\lambda }{\rho }_{\lambda }{| }_{\lambda =0},{R}^{T}G\,R\right].\end{eqnarray} \tag{ 32 }$$

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

$$\begin{eqnarray}&&{\rm{\Phi }}(t)=\displaystyle \frac{{\rm{i}}}{2{\hslash }}\,\mathrm{tr}\left[\left[{\rm{\Sigma }}(t),{R}^{T}g\,R\right]\,{\rho }_{0}\right],\end{eqnarray} \tag{ 33 }$$

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

$$\begin{eqnarray}&&{\rm{\Phi }}(t)=\displaystyle \frac{{\rm{i}}}{{\hslash }}{\left\langle [{\rm{\Sigma }}(t),h]\right\rangle }_{0}.\end{eqnarray} \tag{ 34 }$$

The Kubo response function (33) can be brought into a more useful shape. To this end, let us look at an individual element of Σ, say ${{\rm{\Sigma }}}_{{lm}}=\tfrac{1}{2}({R}_{l}{R}_{m}+{R}_{m}{R}_{l})$. The response function of this object has two parts, i.e. ${{\rm{\Phi }}}_{{{\rm{\Sigma }}}_{{lm}}}(t)=\tfrac{1}{2}\left({{\rm{\Phi }}}_{{R}_{m}{R}_{l}}(t)+{{\rm{\Phi }}}_{{R}_{l}{R}_{m}}(t)\right)$. We have:

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{{R}_{l}{R}_{m}}(t) & = & \displaystyle \frac{{\rm{i}}}{{\hslash }}\,{\left\langle \left[{R}_{l}(t){R}_{m}(t),h\right]\right\rangle }_{{\rho }_{0}}\\ & = & \displaystyle \frac{{\rm{i}}}{2{\hslash }}\displaystyle \sum _{{{nn}}^{{\prime} }}{g}_{{{nn}}^{{\prime} }}{\left\langle \left[{R}_{l}(t){R}_{m}(t),{R}_{n}{R}_{{n}^{{\prime} }}\right]\right\rangle }_{0}\\ & = & \displaystyle \frac{{\rm{i}}}{2\,{\hslash }}\displaystyle \sum _{{l}^{{\prime} }{m}^{{\prime} }{{nn}}^{{\prime} }}{g}_{{{nn}}^{{\prime} }}{\left({S}_{G}^{t}\right)}_{{{ll}}^{{\prime} }}{\left({S}_{G}^{t}\right)}_{{{mm}}^{{\prime} }}{\left\langle \left[{R}_{{l}^{{\prime} }}{R}_{{m}^{{\prime} }},{R}_{n}{R}_{{n}^{{\prime} }}\right]\right\rangle }_{0}\\ & = & \displaystyle \frac{{\rm{i}}}{2\,{\hslash }}\displaystyle \sum _{{l}^{{\prime} }{m}^{{\prime} }{{nn}}^{{\prime} }}{g}_{{{nn}}^{{\prime} }}{\left({S}_{G}^{t}\right)}_{{{ll}}^{{\prime} }}{\left({S}_{G}^{t}\right)}_{{{mm}}^{{\prime} }}\left\{{{\rm{\Omega }}}_{{m}^{{\prime} }{n}^{{\prime} }}{\left\langle {R}_{{l}^{{\prime} }}{R}_{n}\right\rangle }_{0}\right.\\ & & \left.+{{\rm{\Omega }}}_{{m}^{{\prime} }n}{\left\langle {R}_{{l}^{{\prime} }}{R}_{{n}^{{\prime} }}\right\rangle }_{0}+{{\rm{\Omega }}}_{{l}^{{\prime} }{n}^{{\prime} }}{\left\langle {R}_{n}{R}_{{m}^{{\prime} }}\right\rangle }_{0}+{{\rm{\Omega }}}_{{l}^{{\prime} }n}{\left\langle {R}_{{n}^{{\prime} }}{R}_{{m}^{{\prime} }}\right\rangle }_{0}\right\}\\ & = & \displaystyle \frac{{\rm{i}}}{{\hslash }}{\left\{{S}_{G}^{t}{\left\langle {{RR}}^{T}\right\rangle }_{0}{g}^{T}{{\rm{\Omega }}}^{T}{{S}_{G}^{{T}^{t}}}^{}+{S}_{G}^{t}{\rm{\Omega }}g{\left\langle {{RR}}^{T}\right\rangle }_{0}{{S}_{G}^{{T}^{t}}}^{}\right\}}_{{lm}},\end{array}\end{eqnarray} \tag{ 35 }$$

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

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{{\rm{\Sigma }}}_{{lm}}}(t)=\displaystyle \frac{{\rm{i}}}{{\hslash }}{\left({S}_{G}^{t}{\sigma }_{0}{g}^{T}{{\rm{\Omega }}}^{T}{{S}_{G}^{{T}^{t}}}^{}+{S}_{G}^{t}{\rm{\Omega }}g{\sigma }_{0}{{S}_{G}^{{T}^{t}}}^{}\right)}_{{lm}}=\displaystyle \frac{{\rm{i}}}{{\hslash }}{\left({S}_{G}^{t}\tilde{\sigma }{{S}_{G}^{{T}^{t}}}^{}\right)}_{{lm}},\end{eqnarray} \tag{ 36 }$$

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

$$\begin{eqnarray}&&{\rm{\Phi }}(t)=\displaystyle \frac{{\rm{i}}}{{\hslash }}{S}_{G}^{t}\tilde{\sigma }{{S}_{G}^{{T}^{t}}}^{}.\end{eqnarray} \tag{ 37 }$$

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 ς.

The thermalization and/or dynamics of quantum open systems is sometimes described in the collision model framework (see for instance [41–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 $\delta t$ with individual environmental modes. Let σ_E(t) denote the state of the environmental mode that interacts with the system at time t:

$$\begin{eqnarray}&&{\sigma }_{E}(t)=c(t){{\mathbb{I}}}_{2}.\end{eqnarray} \tag{ 45 }$$

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:

$$\begin{eqnarray}&&{\sigma }_{s}(t)\mapsto {\sigma }_{s}(t+\delta t)={\left[{S}_{\eta }({\sigma }_{s}(t)\oplus c(t){{\mathbb{I}}}_{2}){S}_{\eta }^{T}\right]}_{E},\end{eqnarray} \tag{ 46 }$$

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

$$\begin{eqnarray}{S}_{\eta }=\left(\begin{array}{cc}\sqrt{\eta }\,{{\mathbb{I}}}_{2} & \sqrt{1-\eta }\,{{\mathbb{I}}}_{2}\\ -\sqrt{1-\eta }\,{{\mathbb{I}}}_{2} & \sqrt{\eta }\,{{\mathbb{I}}}_{2}\end{array}\right),\ \ \eta \in [0,1].\end{eqnarray} \tag{ 47 }$$

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 $X=\sqrt{\eta }\,{{\mathbb{I}}}_{2}$, and $Y=(1-\eta )c(t){{\mathbb{I}}}_{2}$. We notice that a consecutive application of this map, with fixed c (i.e. ${\sigma }_{E}(t)={c}_{0}{{\mathbb{I}}}_{2}$, ∀t) brings any initial system covariance matrix to the steady state ${\sigma }_{s}(\infty )={c}_{0}{{\mathbb{I}}}_{2}$. 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₀ + λ(t), with $| \lambda (t)| \ll {c}_{0}$. Moreover, the initial covariance matrix of the system is ${\sigma }_{s}(0)={c}_{0}{{\mathbb{I}}}_{2}$. By using our FDT we aim at identifying the response function. First, notice that $\varsigma ={\partial }_{\lambda }{\sigma }_{\lambda }{| }_{\lambda =0}={{\mathbb{I}}}_{2}$. In addition, since $X=\sqrt{\eta }\,{{\mathbb{I}}}_{2}$, we have $\varsigma (t)={X}^{t/\delta t}\varsigma {{X}^{T}}^{t/\delta t}\,{{\mathbb{I}}}_{2}={\eta }^{t/\delta t}\,{{\mathbb{I}}}_{2}$. Therefore, we have ${\rm{\Phi }}(t)=-{\partial }_{t}{\eta }^{t/\delta t}\,{{\mathbb{I}}}_{2}=-\mathrm{log}({\eta }^{1/\delta t}){\eta }^{t/\delta t}\,{{\mathbb{I}}}_{2}$, which vanishes exponentially in time.

Without any loss of generality let $\lambda (t)={\lambda }_{0}\cos \nu t$, with ν 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 ν. To this aim, it is useful to study the linear response of the covariance matrix to the strength of the perturbation:

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{\lambda }{\sigma }_{s}(t){| }_{\lambda =0} & = & -{\displaystyle \int }_{0}^{t}{\partial }_{s}{\eta }^{s/\delta t}\cos \nu (t-s){\rm{d}}{s}\,{{\mathbb{I}}}_{2}\\ & = & \displaystyle \frac{\tilde{\eta }}{{\nu }^{2}+{\tilde{\eta }}^{2}}\,\left[\tilde{\eta }\,\cos \nu t+\nu \sin \nu t-\tilde{\eta }{{\rm{e}}}^{-\tilde{\eta }t}\right]{{\mathbb{I}}}_{2},\end{array}\end{eqnarray} \tag{ 48 }$$

where we define $\tilde{\eta }=-\mathrm{log}({\eta }^{1/\delta t})$ to lighten our notation. In figure 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 ν = 0, specifically illustrates the relaxation of the system to a new thermal state after a quench.

Zoom In Zoom Out Reset image size

Indeed, the linear response is initially zero. For small t, it grows with $\tilde{\eta }t$, 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 $\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 ν = 0. Finally, if $\left|\nu /\tilde{\eta }\right|\to \infty$, the response goes to zero, hence the noise is canceled out.

An OPO coupled to a vacuum field is described by the Hamiltonian $H={\rm{i}}\epsilon ({a}^{\dagger 2}-{a}^{2})/4$, with $\epsilon \geqslant 0$ denoting the effective pump intensity. Here, we use the standard definition of the annihilation and creation operators, that read as $a=(x+{\rm{i}}\,p)/\sqrt{2}$ and ${a}^{\dagger }=(x-{\rm{i}}\,p)/\sqrt{2}$, respectively. The coupling to the vacuum is described by the Lindbladian operator $L=\sqrt{\kappa }\,a$, 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:

$$\begin{eqnarray}G=\displaystyle \frac{\epsilon }{2}\left(\begin{array}{cc}0 & 1\\ 1 & 0\end{array}\right),\ \ C=\sqrt{\displaystyle \frac{\kappa }{2}}{\left(1{i}\right)}^{T}.\end{eqnarray} \tag{ 49 }$$

From here, one can obtain the matrices A and D that appear in equation (43):

$$\begin{eqnarray}A=\left(\begin{array}{cc}\displaystyle \frac{1}{2}(\epsilon -\kappa ) & 0\\ 0 & -\displaystyle \frac{1}{2}(\epsilon +\kappa )\end{array}\right),\ \ D=\displaystyle \frac{\kappa }{2}\,{{\mathbb{I}}}_{2}.\end{eqnarray} \tag{ 50 }$$

Thus, the steady state covariance matrix exists if κ > , and reads as $\sigma =\tfrac{1}{2}\mathrm{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}+{\rm{i}}({L}_{1}^{\dagger }{L}_{2}-{L}_{1}{L}_{2}^{\dagger })/2$, with ${H}_{j}={\rm{i}}{\epsilon }_{j}({a}_{j}^{\dagger 2}-{a}_{j}^{2})/4$, and ${L}_{j}=\sqrt{{\hslash }\kappa }{a}_{j}$. In turn the dissipation is given by a single operator L = L₁ + L₂ For convenience, in what follows we work in a representation, where the quadrature vector is represented by R = (x₁ ... x_N p₁ ... p_N)^T. The matrices G and C read as:

$$\begin{eqnarray}\begin{array}{rcl}G & = & \displaystyle \frac{1}{2}\left(\begin{array}{cccc}0 & 0 & {\epsilon }_{1} & -\kappa \\ 0 & 0 & \kappa & {\epsilon }_{2}\\ {\epsilon }_{1} & \kappa & 0 & 0\\ -\kappa & {\epsilon }_{2} & 0 & 0\end{array}\right),\\ C & = & \sqrt{\displaystyle \frac{\kappa }{2}}{\left(11ii\right)}^{T}.\end{array}\end{eqnarray} \tag{ 51 }$$

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

$$\begin{eqnarray}\begin{array}{rcl}A & = & \left(\begin{array}{cc}\displaystyle \frac{{\epsilon }_{1}-\kappa }{2} & 0\\ -\kappa & \displaystyle \frac{{\epsilon }_{2}-\kappa }{2}\end{array}\right)\,\oplus \,\left(\begin{array}{cc}-\displaystyle \frac{{\epsilon }_{1}+\kappa }{2} & 0\\ -\kappa & -\displaystyle \frac{{\epsilon }_{2}+\kappa }{2}\end{array}\right),\\ D & = & \displaystyle \frac{1}{2}\left(\begin{array}{cc}\kappa & \kappa \\ \kappa & \kappa \end{array}\right)\,\oplus \,\displaystyle \frac{1}{2}\left(\begin{array}{cc}\kappa & \kappa \\ \kappa & \kappa \end{array}\right).\end{array}\end{eqnarray} \tag{ 52 }$$

Again, the criterion for having a steady state is $\kappa \gt \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}_{\alpha }^{-1}$. After doing some straightforward algebra, one finds:

$$\begin{eqnarray}{{\rm{e}}}^{{A}_{x}}=\left(\begin{array}{cc}{{\rm{e}}}^{\tfrac{{\epsilon }_{1}-\kappa }{2}} & 0\\ \displaystyle \frac{-2\kappa }{{\epsilon }_{1}-{\epsilon }_{2}}\left({{\rm{e}}}^{\tfrac{{\epsilon }_{1}-\kappa }{2}}-{{\rm{e}}}^{\tfrac{{\epsilon }_{2}-\kappa }{2}}\right) & {{\rm{e}}}^{\tfrac{{\epsilon }_{2}-\kappa }{2}}\end{array}\right),\end{eqnarray} \tag{ 53 }$$

and

$$\begin{eqnarray}{{\rm{e}}}^{{A}_{p}}=\left(\begin{array}{cc}{{\rm{e}}}^{-\tfrac{{\epsilon }_{1}+\kappa }{2}} & 0\\ \displaystyle \frac{2\kappa }{{\epsilon }_{1}-{\epsilon }_{2}}\left({{\rm{e}}}^{-\tfrac{{\epsilon }_{1}+\kappa }{2}}-{{\rm{e}}}^{-\tfrac{{\epsilon }_{2}+\kappa }{2}}\right) & {{\rm{e}}}^{-\tfrac{{\epsilon }_{2}+\kappa }{2}}\end{array}\right).\end{eqnarray} \tag{ 54 }$$

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

$$\begin{eqnarray}\sigma =\displaystyle \frac{1}{2}\left(\begin{array}{cc}\displaystyle \frac{\kappa }{\kappa -{\epsilon }_{1}} & \displaystyle \frac{-2\kappa {\epsilon }_{1}}{{g}_{-}}\\ \displaystyle \frac{-2\kappa {\epsilon }_{1}}{{g}_{-}} & \displaystyle \frac{-\kappa {h}_{+}}{{g}_{-}}\end{array}\right)\oplus \displaystyle \frac{1}{2}\left(\begin{array}{cc}\displaystyle \frac{\kappa }{\kappa +{\epsilon }_{1}} & \displaystyle \frac{2\kappa {\epsilon }_{1}}{{g}_{+}}\\ \displaystyle \frac{2\kappa {\epsilon }_{1}}{{g}_{+}} & \displaystyle \frac{\kappa {h}_{-}}{{g}_{+}}\end{array}\right),\end{eqnarray} \tag{ 55 }$$

where we define g_± = (₁ + ₂  ±  2κ)(₁  ± κ), and ${h}_{\pm }=({\epsilon }_{1}^{2}+{\epsilon }_{1}\,{\epsilon }_{2}\,\pm \,{\epsilon }_{1}\,\kappa +2{\kappa }^{2}\,\mp \kappa \,{\epsilon }_{2}){\left({\epsilon }_{2}\mp k\right)}^{-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:

$$\begin{eqnarray}&&{\rm{\Phi }}(t)=-{\partial }_{t}\left[{X}_{0}^{t}\,\varsigma \,{{X}_{0}^{{T}^{t}}}^{}\right].\end{eqnarray} \tag{ 56 }$$

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, Φ_x(t) read as:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{x}(t)=-{\partial }_{t}\left[{{\rm{e}}}^{{A}_{x}t}{| }_{\lambda =0}{\left({\partial }_{\lambda }{\sigma }_{x}\right)}_{\lambda =0}{{\rm{e}}}^{{A}_{x}^{T}t}{| }_{\lambda =0}\right],\end{eqnarray} \tag{ 57 }$$

with σ_x being the first matrix in the rhs of equation (55). Setting κ 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 figure 2. Our first observation is that the position quadratures are more responsive to the perturbation, with $\left\langle {x}_{2}^{2}\right\rangle$ having the biggest amplitude of oscillations. From a sensing (estimation) point of view, this means that $\left\langle {x}_{2}^{2}\right\rangle$ 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 $\left\langle {x}_{2}^{2}\right\rangle$). We further notice that the response to a perturbation with bigger ν is smaller—except for the $\left\langle {p}_{1}{p}_{2}\right\rangle$. 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:

$$\begin{eqnarray}&&\chi (\omega )={\int }_{-\infty }^{\infty }{\rm{\Phi }}(t){{\rm{e}}}^{{\rm{i}}\omega t}{\rm{d}}{t}={\int }_{0}^{\infty }{\rm{\Phi }}(t){{\rm{e}}}^{{\rm{i}}\omega t}{\rm{d}}{t},\end{eqnarray} \tag{ 58 }$$

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 $\lambda (t)=\lambda \cos \nu t$, the linear response at long times reads as:

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{\lambda }{\sigma }_{t}{| }_{\lambda =0} & \approx & {\displaystyle \int }_{0}^{\infty }{\rm{\Phi }}(\tau )\cos \nu (t-\tau ){\rm{d}}\tau \\ & = & \cos (\nu \,t){\displaystyle \int }_{0}^{\infty }{\rm{\Phi }}(\tau )\cos (\nu \tau ){\rm{d}}\tau +\sin (\nu \,t){\displaystyle \int }_{0}^{\infty }{\rm{\Phi }}(\tau )\sin (\nu \tau ){\rm{d}}\tau \\ & = & \cos (\nu t)\mathrm{Re}\chi (\nu )+\sin (\nu t)\mathrm{Im}\chi (\nu )=| \chi (\nu )| \cos (\nu t-\alpha ),\end{array}\end{eqnarray} \tag{ 59 }$$

where we define $\alpha ={\cot }^{-1}\left(\tfrac{\mathrm{Im}\chi (\omega )}{\mathrm{Re}\chi (\omega )}\right)$. 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We start by the first moments. Namely, for the R_j element, we can calculate the terms that appear above, individually. For the Hamiltonian part we have

$$\begin{eqnarray}&&{\rm{i}}[H,{R}_{j}]=\displaystyle \frac{{\rm{i}}}{2}{G}_{l,m}[{R}_{l}{R}_{m},{R}_{j}]=\displaystyle \frac{{\rm{i}}}{2}{G}_{l,m}({{\rm{\Omega }}}_{{lj}}{R}_{m}+{{\rm{\Omega }}}_{{mj}}{R}_{l})=-{\rm{i}}{\left({\rm{\Omega }}{GR}\right)}_{j}.\end{eqnarray} \tag{ C.3 }$$

The dissipation part of the dynamics has two contributions. The first part can be written as (we drop the index of the Lindbladian operators to avoid confusion, but we shall reconsider them later on):

$$\begin{eqnarray}\begin{array}{rcl}{L}^{\dagger }{R}_{j}L & = & {R}^{T}{c}^{* }{R}_{j}{c}^{T}R={c}_{l}^{* }{c}_{m}{R}_{l}{R}_{j}{R}_{m}={c}_{l}^{* }{c}_{m}{R}_{l}{R}_{m}{R}_{j}+{c}_{l}^{* }{c}_{m}{{\rm{\Omega }}}_{{jm}}{R}_{l}\\ & = & {c}_{l}^{* }{c}_{m}{R}_{l}{R}_{m}{R}_{j}+{\left({\rm{\Omega }}{{cc}}^{\dagger }R\right)}_{j}.\end{array}\end{eqnarray} \tag{ C.4 }$$

In the same manner, we can deal with the second part:

$$\begin{eqnarray}\begin{array}{rcl}-\displaystyle \frac{1}{2}\left\{{L}^{\dagger }L,{R}_{j}\right\} & = & -\displaystyle \frac{1}{2}\left({R}^{T}{c}^{* }{c}^{T}{{RR}}_{j}+{R}_{j}{R}^{T}{c}^{* }{c}^{T}R\right)=-\displaystyle \frac{1}{2}\left({c}_{l}^{* }{c}_{m}{R}_{l}{R}_{m}{R}_{j}+{c}_{l}^{* }{c}_{m}{R}_{j}{R}_{l}{R}_{m}\right)\\ & = & -\displaystyle \frac{1}{2}\left({c}_{l}^{* }{c}_{m}{R}_{l}{R}_{m}{R}_{j}+{c}_{l}^{* }{c}_{m}{R}_{l}{R}_{j}{R}_{m}+{c}_{l}^{* }{c}_{m}{{\rm{\Omega }}}_{{jl}}{R}_{m}\right)\\ & = & -\displaystyle \frac{1}{2}\left({c}_{l}^{* }{c}_{m}{R}_{l}{R}_{m}{R}_{j}+{c}_{l}^{* }{c}_{m}{R}_{l}{R}_{m}{R}_{j}+{c}_{l}^{* }{c}_{m}{{\rm{\Omega }}}_{{jm}}{R}_{l}+{c}_{l}^{* }{c}_{m}{{\rm{\Omega }}}_{{jl}}{R}_{m}\right)\\ & = & -{c}_{l}^{* }{c}_{m}{R}_{l}{R}_{m}{R}_{j}-\displaystyle \frac{1}{2}{\left({\rm{\Omega }}{c}^{* }{c}^{T}R+{\rm{\Omega }}{{cc}}^{\dagger }R\right)}_{j}.\end{array}\end{eqnarray} \tag{ C.5 }$$

Putting the two parts together, we have:

$$\begin{eqnarray}&&{L}^{\dagger }{R}_{j}L-\displaystyle \frac{1}{2}\left\{{{LL}}^{\dagger },{R}_{j}\right\}=\displaystyle \frac{1}{2}{\left({\rm{\Omega }}{{cc}}^{\dagger }R-{\rm{\Omega }}{c}^{* }{c}^{T}R\right)}_{j}={\left[i{\rm{\Omega }}\mathrm{Im}({{cc}}^{\dagger })R\right]}_{j}.\end{eqnarray} \tag{ C.6 }$$

By considering all of the Lindbladian operators, in a more compact form we have:

$$\begin{eqnarray}&&\dot{R}=\left(-{\rm{i}}{\rm{\Omega }}G\,+{\rm{i}}{\rm{\Omega }}\displaystyle \sum _{k=1}^{m}\,\mathrm{Im}({c}_{k}{c}_{k}^{\dagger })\right)R={AR},\end{eqnarray} \tag{ C.7 }$$

with $A:= -{\rm{i}}{\rm{\Omega }}\left(G\,-\mathrm{Im}({{CC}}^{\dagger })\right)$, and $C:= {\left({c}_{1}^{T};{c}_{2}^{T};...{c}_{m}^{T}\right)}^{T}$ is a 2N × m matrix containing the dissipation coefficients.

We choose the same procedure to deal with the second moments. To begin with, we explore the Hamiltonian term:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{i}}[H,{R}_{j}{R}_{k}] & = & {\rm{i}}\left([H,{R}_{j}]{R}_{k}+{R}_{j}[H,{R}_{k}]\right)=-{\rm{i}}{\left({\rm{\Omega }}{GR}\right)}_{j}{R}_{k}-{{\rm{i}}{R}}_{j}{\left({\rm{\Omega }}{GR}\right)}_{k}\\ & = & -{\rm{i}}{\left({\rm{\Omega }}{{GRR}}^{T}\right)}_{{jk}}-{\rm{i}}{\left({{RR}}^{T}{G}^{T}{{\rm{\Omega }}}^{T}\right)}_{{jk}}.\end{array}\end{eqnarray} \tag{ C.8 }$$

Moreover, for the dissipators we have:

$$\begin{eqnarray}&&{L}^{\dagger }{R}_{j}{R}_{k}L={c}_{m}^{* }{c}_{n}{R}_{m}{R}_{j}{R}_{k}{R}_{n}={c}_{m}^{* }{c}_{n}{R}_{m}{R}_{n}{R}_{j}{R}_{k}+{c}_{m}^{* }{c}_{n}{R}_{m}{R}_{k}{{\rm{\Omega }}}_{{jn}}+{c}_{m}^{* }{c}_{n}{R}_{m}{R}_{j}{{\rm{\Omega }}}_{{kn}}.\end{eqnarray} \tag{ C.9 }$$

Let us deal with the two terms appearing in the anticommutator separately. The first one reads:

$$\begin{eqnarray}&&-\displaystyle \frac{1}{2}{L}^{\dagger }{{LR}}_{j}{R}_{k}=-\displaystyle \frac{1}{2}{c}_{m}^{* }{c}_{n}{R}_{m}{R}_{n}{R}_{j}{R}_{k},\end{eqnarray} \tag{ C.10 }$$

and finally the second one reads as:

$$\begin{eqnarray}\begin{array}{rcl}-\displaystyle \frac{1}{2}{R}_{j}{R}_{k}{L}^{\dagger }L & = & -\displaystyle \frac{1}{2}\,{c}_{m}^{* }{c}_{n}{R}_{j}{R}_{k}{R}_{m}{R}_{n}=-\displaystyle \frac{1}{2}\left\{{c}_{m}^{* }{c}_{n}{R}_{j}{R}_{m}{R}_{n}{R}_{k}+{c}_{m}^{* }{c}_{n}{R}_{j}{R}_{m}{{\rm{\Omega }}}_{{kn}}+{c}_{m}^{* }{c}_{n}{R}_{j}{R}_{n}{{\rm{\Omega }}}_{{km}}\right\}\\ & = & -\displaystyle \frac{1}{2}\left\{{c}_{m}^{* }{c}_{n}{R}_{m}{R}_{n}{R}_{j}{R}_{k}+{c}_{m}^{* }{c}_{n}{R}_{m}{R}_{k}{{\rm{\Omega }}}_{{jn}}+{c}_{m}^{* }{c}_{n}{R}_{n}{R}_{k}{{\rm{\Omega }}}_{{jm}}+{c}_{m}^{* }{c}_{n}{R}_{j}{R}_{m}{{\rm{\Omega }}}_{{kn}}\right.\\ & & \left.+{c}_{m}^{* }{c}_{n}{R}_{j}{R}_{n}{{\rm{\Omega }}}_{{km}}\right\}.\end{array}\end{eqnarray} \tag{ C.11 }$$

Putting the last three expressions together in the dissipation part of the master equation yields:

$$\begin{eqnarray}\begin{array}{rcl}{L}^{\dagger }{R}_{j}{R}_{k}L-\displaystyle \frac{1}{2}\left\{{L}^{\dagger }L,{R}_{j}{R}_{k}\right\} & = & {{\rm{\Omega }}}_{{jn}}{c}_{m}^{* }{c}_{n}{R}_{m}{R}_{k}+{{\rm{\Omega }}}_{{kn}}{c}_{n}{c}_{m}^{* }{R}_{m}{R}_{j}-\displaystyle \frac{1}{2}\left\{{{\rm{\Omega }}}_{{jm}}{c}_{m}^{* }{c}_{n}{R}_{n}{R}_{k}\right.\\ & & \left.+{{\rm{\Omega }}}_{{jn}}{c}_{n}{c}_{m}^{* }{R}_{m}{R}_{k}+{{\rm{\Omega }}}_{{km}}{c}_{m}^{* }{c}_{n}{R}_{j}{R}_{n}+{{\rm{\Omega }}}_{{kn}}{c}_{n}{c}_{m}^{* }{R}_{j}{R}_{m}\right\}\\ & = & {\left[{\rm{\Omega }}{{cc}}^{\dagger }{{RR}}^{T}\right]}_{{jk}}-{\left[{{RR}}^{T}{c}^{* }{c}^{T}{\rm{\Omega }}\right]}_{{jk}}+{\left[{\rm{\Omega }}{c}^{* }{c}^{T}{\rm{\Omega }}\right]}_{{jk}}\\ & & -\displaystyle \frac{1}{2}{\left[{\rm{\Omega }}{c}^{* }{c}^{T}{{RR}}^{T}\right]}_{{jk}}-\displaystyle \frac{1}{2}{\left[{\rm{\Omega }}{{cc}}^{\dagger }{{RR}}^{T}\right]}_{{jk}}\\ & & +\displaystyle \frac{1}{2}{\left[{{RR}}^{T}{{cc}}^{\dagger }{\rm{\Omega }}\right]}_{{jk}}+\displaystyle \frac{1}{2}{\left[{{RR}}^{T}{c}^{* }{c}^{T}{\rm{\Omega }}\right]}_{{jk}}\\ & = & \displaystyle \frac{1}{2}{\left[{\rm{\Omega }}{{cc}}^{\dagger }{{RR}}^{T}-{\rm{\Omega }}{c}^{* }{c}^{T}{{RR}}^{T}\right]}_{{jk}}+\displaystyle \frac{1}{2}{\left[{{RR}}^{T}{{cc}}^{\dagger }{\rm{\Omega }}-{{RR}}^{T}{c}^{* }{c}^{T}{\rm{\Omega }}\right]}_{{jk}}\\ & & +{\left[{\rm{\Omega }}{c}^{* }{c}^{T}{\rm{\Omega }}\right]}_{{jk}}\\ & = & {\rm{i}}{\left[{\rm{\Omega }}\mathrm{Im}({{cc}}^{\dagger }){{RR}}^{T}\right]}_{{jk}}+{\rm{i}}{\left[{{RR}}^{T}\mathrm{Im}({{cc}}^{\dagger }){\rm{\Omega }}\right]}_{{jk}}+{\left[{\rm{\Omega }}{c}^{* }{c}^{T}{\rm{\Omega }}\right]}_{{jk}}.\end{array}\end{eqnarray} \tag{ C.12 }$$

If we add this to the same dissipation part, but with the elements $j\leftrightarrow k$, and divide by two, we have:

$$\begin{eqnarray}&&{L}^{\dagger }{{\rm{\Sigma }}}_{{jk}}L-\displaystyle \frac{1}{2}\left\{{L}^{\dagger }L,{{\rm{\Sigma }}}_{{jk}}\right\}={\rm{i}}{\left[{\rm{\Omega }}\mathrm{Im}({{cc}}^{\dagger }){\rm{\Sigma }}\right]}_{{jk}}+{\rm{i}}{\left[{\rm{\Sigma }}\mathrm{Im}({{cc}}^{\dagger }){\rm{\Omega }}\right]}_{{jk}}+{\left[{\rm{\Omega }}\mathrm{Re}({{cc}}^{\dagger }){\rm{\Omega }}\right]}_{{jk}}.\end{eqnarray} \tag{ C.13 }$$

Thus, for the second moments, after adding up all dissipators, we will have the following master equation:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\rm{\Sigma }}}{{\rm{d}}{t}}=A{\rm{\Sigma }}+{\rm{\Sigma }}{A}^{T}+D,\end{eqnarray} \tag{ C.14 }$$

with $A=-{\rm{i}}{\rm{\Omega }}(G-\mathrm{Im}({{CC}}^{\dagger }))$, and $D={\rm{\Omega }}\,\mathrm{Re}({{CC}}^{\dagger }){\rm{\Omega }}$, with the same definition of C that we have for first moments, i.e. $C:= {\left({c}_{1}^{T};{c}_{2}^{T};...{c}_{m}^{T}\right)}^{T}$.