Generation of stable Gaussian cluster states in optomechanical systems with multifrequency drives

We show how to dissipatively stabilize the quantum state of N mechanical resonators in an optomechanical system, where the resonators interact by radiation pressure with N optical modes, which are driven by properly selected multifrequency drives. We analyze the performance of this approach for the stationary preparation of Gaussian cluster states.

Here, we employ reservoir engineering techniques, based on the use of multifrequency drives, to control the stationary quantum properties of many mechanical resonators.This approach builds upon several previous results regarding the entanglement of mechanical resonators [49,50,51,52,53,55,54,56,16,23,22].In the case of two mechanical modes, these approaches work by driving the system with laser fields resonant with both the red and blue mechanical sidebands, such that specific collective Bogoliubov modes can be cooled, resulting in two-mode mechanical squeezing, that is, entanglement.
In the present work, we study a multimode optomechanical system and show that in order to achieve full steady-state control of N mechanical resonators, one needs N optical modes, each driven on both the red and blue sidebands corresponding to all the mechanical resonators (see Fig. 1).This requires the use of drivings with 2 N 2 frequency components.Thereby, the mechanical resonators can be laser-cooled to any Gaussian state, depending on the corresponding driving amplitudes and phases.In particular, we analyze in detail the conditions for the stationary generation of Gaussian cluster states [61,62,63], which constitute the fundamental resource of measurement-based quantum computation in the continuous variable setting [64,65].Differently from the more common approach of gatebased quantum computation, where a sequence of unitary quantum operations is applied on a quantum processor, in measurement-based quantum computation, a quantum algorithm is realized by performing a sequence of measurements over certain highly entangled quantum states.These are the so-called cluster states, characterized by graphs that identify the structure of the entanglement between the components of the quantum processor.Therefore, in this approach, a fundamental part of the computation is the preparation of the cluster state.In the continuous variable setting, very large Gaussian cluster states of optical fields have been demonstrated [66].Here, we suggest a path to achieving similar results with mechanical systems, thereby providing a step forward towards the realization of quantum computation over mechanical degrees of freedom [59,60].
Our approach is somewhat similar to that of Ref. [55] where, however, a single cavity mode is considered, so that the protocol requires the sequential application of specific driving pulses.In our proposal, on the contrary, the state preparation is truly stationary.In this respect, our result is analogous to the scheme proposed in Ref. [67], which shows how to stabilize the state of many microwave modes.
The article is organized as follows.In section 2, we introduce the system, discuss the approximations at the basis of the analytical model, and analyze the conditions under which it is possible to have full control of the mechanical steady state.Then, in section 3, we specialize these results to the dissipative preparation of Gaussian cluster states.The performance of the scheme is analyzed in the next section 4. Finally, in section 5, we draw our conclusions.In the appendices, we report additional details about the derivation of the analytical model (Appendix A), and of the corresponding steady state (Appendix B).

Model and basic dynamics
We consider an optomechanical system that involves the interaction of N resonant modes of a mechanical object (as, for example, in Refs.[3,4,5,7,6,29]) at frequencies Ω j and N resonant modes of an optical cavity at the frequencies ω k , for k, j ∈ {1, . . ., N}.In this system, each cavity mode interacts with all the mechanical modes with strengths g k j .The system dynamics are controlled by driving the cavity with laser fields constituted of 2N 2 spectral components at the frequencies λ km , where k ∈ {1, . . ., N} and m ∈ {1, . . ., 2N}.The frequencies where ϵ km is the complex amplitude of the spectral component at frequency λ km .
As is customary with optomechanical systems, we consider the fluctuations of the system variables around the corresponding average values [68].These fluctuations are described by the bosonic annihilation operators a k for the optical modes and b j for the mechanical modes.As detailed in Appendix A, the equations for the fluctuations can be linearized for large values of the average fields.In particular, simple equations can be obtained by assuming sufficiently small couplings g k j and expanding the average fields at the lowest relevant order in g k j .Thereby, including dissipation and noise of the optical and mechanical modes at rates κ k and γ j , respectively, the quantum Langevin equations for the slowly varying operators defined by a k = e −i ω k t a k [where ω k includes the frequency shift due to the optomechanical interaction, see Eq. (A.13)] and b j = e −i Ω j t b j , take the form (see Appendix A.1) where and where a in k and b in j are the input noise operators characterized by the correlation functions with n j the thermal phonon number at temperature T defined by The driving frequencies can be selected in order to make resonant only those specific terms that allow for full control of the mechanical dynamics, while the effects of the residual time-dependent non-resonant terms are negligible whenever their effective coupling strengths are much smaller than the corresponding frequencies.To be specific, we select, for m ≤ N, such that each cavity mode is driven both on the red and blue mechanical sidebands corresponding to each mechanical resonator (see Fig. 1).We also assume that the system is in the resolved sideband regime [68] and that all the mechanical frequencies are sufficiently different from each other so that (see also Appendix A.2) This relation justifies a rotating-wave approximation under which all the non-resonant terms are neglected.Consequently, the quantum Langevin equations can be further approximated and rewritten as follows where we have introduced the matrices of interaction coefficients and the collective coupling strengths When the strengths and phases of the driving fields, ϵ km , are properly selected such that [see Eq. ( 3) and ( 8)] then these matrices define a Bogoliubov transformation, and the corresponding collective mechanical Bogoliubov modes fulfill the equations where and the collective noise operators are In particular, we note that if all the mechanical dissipation rates are equal, i.e. γ j = γ for all j, then, according to Eq. ( 10), W = γ 2 1 and T = 0. Eq. ( 12) describes the laser cooling of all the mechanical Bogoliubov modes c j .If the values of γ j (i.e., the couplings to the thermal reservoir) are negligible, these modes are cooled to the vacuum of modes c j , characterized by the correlation functions The corresponding result for the original modes (when These expressions are equal to the correlations evaluated over a pure Gaussian state |ψ⟩ (such as ⟨ψ|b j b † j ′ |ψ⟩) defined in terms of a Gaussian unitary operator U applied to the vacuum of the b modes when U is characterized by the relation [see Eq. ( 11)] such that the inverse transformation is given by This means that, by properly tuning the interaction strengths, which constitute the matrices X and Y, it is possible to drive the system to pure stationary Gaussian states of the form of Eq. ( 17).In particular, since we have no limitations over the matrices X and Y, which determine the Bogoliubov transformation, this scheme allows for the dissipative preparation of any pure zero-average Gaussian state of the mechanical modes.In other terms, this scheme realizes a sort of laser cooling to any collective mechanical zero-average Gaussian state.
In the next sections, we apply this result to the dissipative preparation of Gaussian cluster states.

Stabilization of cluster states
Let us now study the performance of our scheme for the preparation of a Gaussian cluster state with an adjacency matrix A (symmetric and with A j, j = 0).It is defined as the state generated by applying the operator [61,62,63,64,65] with x j = b j + b † j , to a set of bosonic modes each in the zero eigenstate, |0⟩ p , of the momentum operators p j = −i b j + i b † j , i.e. |ψ cluster ⟩ = C |0⟩ p .These are infinitely squeezed states and, in practice, one can prepare approximate cluster states with finite squeezing.Here we study the preparation of states of the form where S r is the squeezing operator S r = e j , and |0⟩ is the vacuum.This state approaches an ideal cluster state in the limit of infinite squeezing, i.e. |ψ cluster ⟩ = lim r→∞ |ψ r ⟩.
The model of the previous section stabilizes a cluster state (21) when the unitary matrix that generates the transformation in Eq. ( 19) is Using this relation, together with Eqs. ( 3), ( 8), ( 11) and (18), one can determine the strength and phases of the driving fields for the stabilization of a cluster state.In detail, first, we compute how the unitary C S r transforms the operator b j .And we find with c r = cosh(r) and s r = sinh(r) (where we used the relations . Then, we compare this result to Eq. ( 19), observing that the two are equal, meaning the steady state of our system approximates a cluster state, when the following expressions are satisfied Finally, given any value of r and any matrix A (which determine the expected cluster state), and any set of real positive values of g k (which determines the effective collective linearized optomechanical interaction strengths [see Eq. ( 9)]), we can use Eq. ( 24), together with Eqs. ( 3), ( 5) and ( 8), to determine the strengths and phases, ϵ k,m , of all the driving fields needed for the preparation of the cluster state (21).They are given by It is worth noting that the values of r and g k can be considered arbitrary only up to a certain extent.In fact, our model is applicable only when the relation in Eq. ( 6) for the validity of the rotating wave approximation holds.In the present case, this condition reduces to for all k, m, m ′ with m m ′ (see also Appendix A.2).In particular, if we assume that all the optomechanical couplings are equal, i.e. g k, j = g k ′ , j ′ for all k, k ′ , j, j ′ , and g k = g k ′ ≡ g for all k, k ′ , and that the mechanical frequencies vary linearly with the mode index, i.e Ω j = j Ω, then the rotating wave approximation is valid when This highlights the fact that both g and r can not be taken too large.

Results
Let us now quantify the performance of our protocol for the preparation of cluster states of various dimensions and geometries.We characterize the performance of our protocol in terms of the fidelity to prepare the state (see Sec. 4.1 below) and in terms of the corresponding variance of the nullifiers (see Sec. 4.2 below).

Fidelity
We analyze the fidelity between the expected state (21) and the one obtained at the steady state of the quantum Langevin equations (12).The fidelity can be computed with the formula [69] Fidelity between the steady state of the engineered dynamics and the expected state, Eq. ( 28), as a function of the mechanical dissipation rate γ ≡ γ j for all j, for cluster states composed of (a) N = 4 modes and (b) N = 10 modes, and corresponding to linear (red, upper lines), rectangular (yellow, middle lines), and fully connected (blue, lower lines) graphs as depicted in the upper part of the figure.The nodes of the graphs represent the mechanical modes.The entry A j, j ′ of the adjacency matrices [introduced in Eq. ( 20)] corresponding to these graphs is equal to A j, j ′ = 1 if there is an edge connecting the nodes j and j ′ , and equal to A j, j ′ = 0 otherwise.The other parameters are r = 2, temperature = 10mK, where V (c)  r and V (c) st are the covariance matrices of the position x (c) j = c j + c † j and momentum p (c) j = −i c j + i c † j operators for, respectively, the expected states (21) and the steady state of the system dynamics.In particular, when Eq. ( 22) is true, the correlations of the collective mechanical operators c j , defined by Eq. ( 18), over the cluster state (21), are equal to Eq. (15).And this entails that Moreover, the matrix V (c) st can be computed by solving the Eq. ( 12) using standard techniques as discussed in Appendix B.

Variance of the nullifiers
A cluster state can be equivalently defined as the zero eigenstate of the nullifiers such that X j |ψ cluster ⟩ = 0 for all j.This entails that for realistic approximated cluster states, these operators are quantum squeezed (i.e., their variance is below the vacuum noise level).Consequently, the variance of the nullifiers constitutes a valid metric to characterize the quality of a cluster state.In particular, one finds that the variances of the nullifiers over the expected, approximated cluster state ( 21) are   28), as a function of the temperature.The other parameters and line styles are as in Fig. 2, and γ j = 5 × 10 −6 κ for all j.= e −2 r ∀ j .
We note that this quantity is equal to the squeezing value of the original modes over which one applies the operator C (20) to define a cluster state [see Eq. ( 21)].In other terms, it is equal to the variance of the momentum quadratures p j of the original modes over the state S r |0⟩, i.e., ⟨0|S † r p 2 j S r |0⟩ = e −2 r .This quantity is used to estimate the effectiveness of a cluster state as a resource for measurement-based quantum computation.Specifically, various estimates [70,71,72,73] indicate that a value of this quantity between -10 and -20 dB [that is, a value of 10 log 10 e −2 r between -10 and -20] is needed to achieve fault-tolerant measurement-based quantum computation.We also note that Ref. [72] suggests that the value of squeezing ⟨0|S † r p 2 j S r |0⟩ determines the potentiality of a cluster state regardless of the  Fidelity, Eq. ( 28), as a function of the number of modes N. Upper (red), middle (yellow) and lower (blue) marks are, respectively, for the linear, rectangular and fully connected geometry.The other parameters and line styles are as in Fig. 2, and γ j = 5 × 10 −6 κ for all j.purity of the state itself.Therefore, an approximate cluster state can be a useful resource for quantum computation even if the fidelity with a pure, ideal cluster state is particularly low, provided that the squeezing is sufficiently large.For this reason, the value of the variance of the nullifiers X 2 j st over the steady state represents a useful estimate of the quality of the steady state as a valid resource for quantum computation.It can be expressed as follows.We  .Fidelity, Eq. ( 28), as a function of g = g k , for all k.The other parameters and line styles are as in Fig. 2, and γ j = 5 × 10 −6 κ for all j.
consider the vector of operators x = (x 1 , . . .x N , p 1 , . . .p N ) T and the N × 2N matrix such that the nullifiers can be written as X j = {O x} j .Correspondingly, X 2 j st can be expressed, in terms of the matrix of correlations associated to the original modes

Numerical results
We compute numerically the covariance matrix V (c) st , of the collective modes and the corresponding matrix V (b)  st for the original modes, as detailed in (Appendix B).These matrices are used to evaluate the fidelity using Eqs.( 28) and ( 29) and the variance of the nullifiers using Eq.(33).We show the results for three different graph geometries, linear (red lines), rectangular (yellow lines) and fully connected (blue lines) and various number of modes, as depicted in the upper part of Fig. 2.
First, we examine the effect of noise.The system reaches the expected state (21) exactly, only if the mechanical dissipation rates γ j are negligible (see Figs. 2 and 3).At finite γ j the final state is closer to the expected result the smaller the temperature, as depicted in Figs. 4 and  5. It's important to note that in general, even at zero temperature, quantum noise prevents the exact generation of the expected state.However, sufficiently large mechanical quality factors, as those used in Figs.4-9, allow for very good cluster state preparation at low temperatures.Furthermore, the detrimental effect of noise becomes more pronounced with a larger number of modes and for cluster states described by graphs with larger number of connections.In particular, we analyze the scaling with system size in Fig. 6.We observe that the preparation is not very sensitive to the number of modes itself but depends more strongly on the number of connections of each mode.This is evident comparing the results of the linear and rectangular clusters with those of the fully connected one, which exhibit much lower fidelity at large N.
In general, we observe good state preparation when the light-induced mechanical dissipation rate, which is of the order of 4 g 2 k /κ k [68], is sufficiently large to overcome the effect of thermal noise, which affects the dynamics of the collective modes (11).Here, thermal noise is characterized by the correlations of the collective noise operators f † j (t) f j (t ′ ) = δ(t − t ′ ) Ξ j , where [see Eqs. ( 14), ( 24), (A.20)-(A.22)] and where N j , in the subscript of the summation, indicates the set of modes adjacent to the mode b j according to the matrix A. Therefore, the condition for a faithful preparation of a cluster state is 4 with Ξ j defined in Eq. ( 34).This quantity is the effective or quantum cooperativity [74] for the collective mechanical Bogoliubov mode c j (11).We remark that it depends not only on the natural mechanical thermal noise, but also on the cluster state we aim to achieve through the factor e 2r and the sets of adjacent modes N j .We also note that, in the case of a fully connected cluster, Ξ j is independent of j and is given by for all j.In Figs.2-5 and 9, we emphasize the significance of the condition in Eq. (35), by adding vertical lines corresponding to the values where the quantum cooperativity of a fully connected graph is one, that is where Ξ ⋆ = 4 g 2 k /κ k .For these values, the fidelity starts to increase to noticeable amounts, but, at the same time, the variance of the nullifiers is already very low, indicating strong squeezing.Specifically, this value is at the level of recent estimates [71,73] for the thresholds of squeezing required for fault-tolerant quantum computation.Therefore, although not exactly equal to the expected state, the steady state is already a good approximation of a cluster state.These results are evaluated for mechanical frequencies in the range 2π × 10 − 2π × 100MHz, in the resolved sideband regime, and at a temperature of 10mK (except for the results shown as a function of temperature itself in Figs.(4) and 5).The corresponding mechanical quality factors, in Figs.4-9, range between Ω 1 /γ 1 = 10 7 and Ω 20 /γ 20 = 2 × 10 8 .However, Fig. 3 show that the variance of the nullifiers remains quite low also for quality factors roughly one order of magnitude smaller, in the range 10 6 − 10 7 .These values are sufficient to achieve good preparation of cluster states with r = 2.For larger values of r the fidelity decreases (see Fig. 7) because, as discussed above, the effective noise of the Bogoliubov modes increases [see Eqs.(34) and (35)].Nevertheless, at the same time the variance of the nullifiers remains very low (see Fig. 8), indicating that the steady state is still a good approximation of a cluster state.
To achieve higher fidelities one should either reduce the value of squeezing r or select larger mechanical frequencies and/or larger values of g.However, increasing both g and r simultaneously may be problematic due to the requirements for the validity of the rotating wave approximation expressed in Eq. ( 27).This condition remains valid for all the values of r reported in Figs.7 and 8 (where for r = 4, we find g e r 2 Ω ∼ 0.09), and for all the values of g used in Fig. (9) (which is evaluated for r = 2).However, for larger value of r and g, the validity of the rotating wave approximation can be questionable unless one increases also the mechanical frequencies (or, more precisely, the value of min j j ′ ∈{1,...N} Ω j − Ω j ′ ).

Experimental implementation
We comment here on the possibility of observing and applying these results in real experiments.The parameters used for the numerical simulation, though challenging, are within present-day technological capabilities in both optical and microwave regimes.Mechanical resonator frequencies in optomechanics can vary widely, from tens of kHz to tens of GHz [68].Notably, various multimode optomechanical experiments have reported mechanical frequencies in the tens of MHz, similar to those used in our simulations [1,2,4,4,5,16,11,17,22,23,24].While the quality factors in our simulation are quite large, it's worth noting that the largest reported mechanical quality factor to date is > 10 9 [75,76], which is one to three orders of magnitude larger than those used here.
Additionally, our numerical results correspond to cryogenic temperatures.The temperature of 10 mK used in Figs. 2, 3, 6-9 is demanding, but we observe that in Fig. 5, a high level of squeezing (and thus good cluster state preparation) is achieved even at temperatures around 100 mK.Moreover, the temperature requirement can be further relaxed by employing higher mechanical frequencies in the GHz range [15,21,26,27,29].
The generated state can be probed by extending the approach discussed in Ref. [53].One should utilize N additional cavity modes and drive them with additional probe fields resonant with the red mechanical sidebands to transfer the mechanical state to the electromagnetic field.The reflected light can then be measured by homodyne detection, and the mechanical covariance matrix can be reconstructed.A similar detection strategy can be employed to perform the measurements needed for a quantum algorithm when using this state for measurement-based quantum computation.

Conclusions
In conclusions, we have studied a scheme for the dissipative stabilization of the quantum state of a multimode optomechanical system, constituted of N mechanical and N optical modes.We have identified the conditions that the driving fields must satisfy in order to achieve full control over the quantum properties of the collective steady state of N mechanical resonators.We have demonstrated that this can be attained for mechanical resonators of different frequencies in the resolved sideband regime, by employing driving fields with 2N 2 frequency components that are resonant with all the blue and red sidebands of all the optical modes, corresponding to all the mechanical modes.In particular, we have applied this result to the stabilization of mechanical cluster states and we have identified the corresponding conditions for the strengths of all the driving components.We have shown that we achieve good preparation, in terms of the fidelity and of the variance of the nullifiers, using system parameters within present day technological capabilities.with M given by the block matrix where K and G are diagonal matrices with elements K k,k = κ k and G k,k = g k , 0 indicates the null matrix, and W and T are defined in Eq. ( 13).Moreover f is the vector of input noise operators f = √ κ 1 a in 1 , . . ., √ κ N a in N , f 1 , . . ., f N , √ κ 1 a in 1 † , . . ., f † 1 , . . .

T
, with f k given in Eq. ( 14).The correlations of these noise operators can be expressed as where N is the correlation matrix We can also write the 2N × 2N matrix for the steady state correlations, solely for the mechanical degrees of freedom, as where Z is the 2N × 4N matrix Z = 0 1 0 0 0 0 0 1 .(B.8) The steady state covariance matrix V (c) st , the elements of which are expressed in terms of the vector of operators x (c) = x (c) 1 , . . .

Figure 1 .
Figure 1.Frequency scheme.N optical modes at frequencies ω k and linewidth κ k , for k ∈ {1, . . .N}, interact with N mechanical modes at frequencies Ω j , for j ∈ {1, . . .N}.The system is driven by multifrequency drives (upwards gray arrows) resonant with all the red and blue sidebands of all optical and mechanical modes, and with strengths ϵ k, j and ϵ k, j+N .
Here and in the following the symbol δ k,k ′ indicates a Kronecker delta while δ(t − t ′ ) a Dirac delta. b

Figure 2 .
Figure 2. Fidelity between the steady state of the engineered dynamics and the expected state, Eq. (28), as a function of the mechanical dissipation rate γ ≡ γ j for all j, for cluster states composed of (a) N = 4 modes and (b) N = 10 modes, and corresponding to linear (red, upper lines), rectangular (yellow, middle lines), and fully connected (blue, lower lines) graphs as depicted in the upper part of the figure.The nodes of the graphs represent the mechanical modes.The entry A j, j ′ of the adjacency matrices [introduced in Eq. (20)] corresponding to these graphs is equal to A j, j ′ = 1 if there is an edge connecting the nodes j and j ′ , and equal to A j, j ′ = 0 otherwise.The other parameters are r = 2, temperature = 10mK, Ω j = j Ω with Ω = 2π × 10MHz, κ ≡ κ k = 0.02Ω, and g ≡ g k = 0.16κ for all k.The vertical lines indicate the values for which Ξ ⋆ = 4 g 2 k /κ k (quantum cooperativity for the fully connected graph equal to one) [see Eqs.(35) and(36)].
Figure 2. Fidelity between the steady state of the engineered dynamics and the expected state, Eq. (28), as a function of the mechanical dissipation rate γ ≡ γ j for all j, for cluster states composed of (a) N = 4 modes and (b) N = 10 modes, and corresponding to linear (red, upper lines), rectangular (yellow, middle lines), and fully connected (blue, lower lines) graphs as depicted in the upper part of the figure.The nodes of the graphs represent the mechanical modes.The entry A j, j ′ of the adjacency matrices [introduced in Eq. (20)] corresponding to these graphs is equal to A j, j ′ = 1 if there is an edge connecting the nodes j and j ′ , and equal to A j, j ′ = 0 otherwise.The other parameters are r = 2, temperature = 10mK, Ω j = j Ω with Ω = 2π × 10MHz, κ ≡ κ k = 0.02Ω, and g ≡ g k = 0.16κ for all k.The vertical lines indicate the values for which Ξ ⋆ = 4 g 2 k /κ k (quantum cooperativity for the fully connected graph equal to one) [see Eqs.(35) and(36)].

10 Figure 3 .
Figure 3. Variance of the nullifiers over the steady state ⟨X 2 j ⟩ st (solid lines), and corresponding value for the expected state (dashed lines), ⟨X 2 j ⟩ r [see Eqs.(31) and (33)] evaluated in dB, as a function of γ and for the parameters of Fig. 2. Panels (a), (b), and (c) are for N = 4. Panels (d), (e) and (f) are for N = 10.Red [(a) and (d)], yellow [(b) and (e)], and blue [(c) and (f)] lines are respectively for the linear, rectangular, and fully connected graphs reported in the upper part of Fig. 2. Different shades of red and yellow indicate the results for different nullifiers [different j, see Eq. (30)] as reported in the insets.In the case of the fully connected graphs [blue lines, (c) and (f)], the result is the same for all the nullifiers.The vertical lines indicate the values for which Ξ ⋆ = 4 g 2 k /κ k (quantum cooperativity for the fully connected graph equal to one) [see Eqs.(35) and(36)].

Figure 4 .
Figure 4. Fidelity, Eq. (28), as a function of the temperature.The other parameters and line styles are as in Fig.2, and γ j = 5 × 10 −6 κ for all j.

10 Figure 5 .
Figure 5. Variance of the nullifiers as a function of the temperature.The other parameters and line styles are as in Fig. 3, and γ j = 5 × 10 −6 κ for all j.
Figure 6.Fidelity, Eq. (28), as a function of the number of modes N. Upper (red), middle (yellow) and lower (blue) marks are, respectively, for the linear, rectangular and fully connected geometry.The other parameters and line styles are as in Fig.2, and γ j = 5 × 10 −6 κ for all j.

Figure 7 .
Figure 7. Fidelity, Eq. (28), as a function of r.The other parameters and line styles are as in Fig.2, and γ j = 5 × 10 −6 κ for all j.

10 Figure 8 .
Figure 8. Variance of the nullifiers as a function of the squeezing parameter r.The other parameters and line styles are as in Fig.3, and γ j = 5 × 10 −6 κ for all j.

Figure 9
Figure 9. Fidelity, Eq. (28), as a function of g = g k , for all k.The other parameters and line styles are as in Fig.2, and γ j = 5 × 10 −6 κ for all j.

4 )
with the matrices Ξ µν defined in Eq. (A.21).The corresponding equation for the 4N × 4N correlation matrix, C, with elements C j,k = a j a k , takes the form Ċ = M C + C M T + N .(B.5) Introducing the linear operator LC ≡ M C + C M T , we can write the steady state solution as C st = −L −1 N .(B.6)

/ 2 ,
x (c) N , p (c) 1 , . . .p (c) N T as V (c) st ℓ,ℓ ′ = x (c) ℓ x (c) ℓ ′ st + x (c) ℓ ′ x (c)ℓ st can be written asV (c) st = R C (c) st + C (c) x (c) = R a.Finally the corresponding covariance matrix for the original modes V(b)  st can be found introducing the Bogoliubov matrix B which transforms the field operators c j and c † j into b j and b† j , i.e. b j = j ′ B j, j ′ c j ′ + B j, j ′ +N c † j ′ and b † j = j ′ B j+N, j ′ c j ′ + B j+N, j ′ +N c † j ′ , that is [see Eqs.(18) and (19)]B = X † −Y T −Y † X T .