Dynamical and quantum effects of collective dissipation in optomechanical systems

Cavity optomechanical systems (OMs) consist of a set of cavity light modes coupled to one or more mechanical elements typically by radiation pressure forces. OMs encompass many physical implementations ranging from the canonical Fabry–Perot resonator with a moving end mirror, to intracavity membranes, to the co-localized photonic and phononic modes of optomechanical crystals, to mention some [1]. The nonlinear character of the optomechanical interaction enables OMs to work in various dynamical regimes and to exhibit a rich dynamical behavior. OMs can settle both in fixed points and display phenomena such as optical bistability [2, 3], as well as into limit cycles or self-sustained oscillations [4, 5], in which dynamical multistability is found [6–8]. Moreover, arrays of coupled OMs can synchronize their self-sustained oscillations, as it has been shown both theoretically [9–11] and experimentally [12–14]. Spontaneous (or mutual) quantum synchronization (recently overviewed in [15]) has also been considered in several systems [15–26], also for common bath (CB) dissipation [16–18, 20–23], being OMs a promising platform to study this phenomenon [19–21, 24–26]. In the quantum regime optomechanical cooling allows efficient ground state cooling of the mechanical modes deep in the resolved sideband regime [27, 28], and mechanical occupation numbers below one have been experimentally achieved [29, 30]. A variety of non-classical states and correlations have been predicted for OMs [1], some of them being recently observed [31–33]. In particular entanglement in the asymptotic state has been considered between light and mirror [34, 35], and membranes [36, 37], even at finite temperatures.

The ubiquitous interaction of a system with its surroundings introduces noise and damping, and might lead to the complete erasure of quantum coherence [38]. In the case of a spatially extended multipartite system—e.g. (optomechanical) array—the spatial structure of the environment—like a field, a lattice or a photonic crystal—becomes crucial in determining different dissipation scenarios. Two prototypical models are often considered: the separate bath model (SB) in which the units of the system dissipate into different uncorrelated environments, and the CB model in which they dissipate collectively into the same environment. The particular dissipation scenario influences deeply the system. While SB dissipation usually destroys quantum correlations [38], dissipation in a CB leads to different results and enables phenomena such as decoherence free/noiseless subspaces [39], dark states [40], superradiance [38, 41], dissipation-induced synchronization of linear networks of quantum harmonic oscillators [16, 18] and of non-interacting spins [22], or no sudden-death of entanglement in systems of decoupled oscillators [42, 43], that would not be present in these systems for SB dissipation. CB dissipation was first considered in the context of superradiance of atoms at distances smaller than the emitted radiation wavelength [41] and often assumed to arise when the spatial extension of a system is smaller than the correlation length of the structured bath. A recent analysis of the CB/SB crossover in a lattice environment reveals the failure of this simple criteria and that actually collective dissipation can emerge even at large distances between the system's units (also for 2D and 3D environments [17, 43]).

A detailed analysis of the CB/SB crossover in optomechanical arrays has not yet been reported and most of the works generally assume independent dissipation affecting optical and mechanical units. On the other hand, collective mechanical dissipation has been recently reported in some experimental platforms [44, 45], such as OMs composed by two coupled nanobeams in a photonic crystal platform (environment). Photonic crystals are indeed also a known tool to suppress mechanical dissipation [46, 47]. In the peculiar case of devices [44, 45] the main dissipation mechanism of the nanobeams is found to be the emission of elastic radiation by the center of mass coordinate of the beams. Furthermore, an analogous collective dissipation mechanism has been experimentally observed in piezoelectric resonators with a similar geometry [48, 49]. In general collective dissipation can lead to quantitative and qualitative differences and we show how it can be beneficial for collective phenomena and quantum correlations. In this work we investigate the effects of CB dissipation in OMs. We focus on a particular OM consisting of two coupled mechanical elements each optomechanically coupled to a different optical mode (figure 1). We consider two dissipation schemes: independent SB for all optical and mechanical modes (figure 1(a)), or CB for the mechanical modes and SB for the optical ones (figure 1(b)). We analyze the effects of collective dissipation both in the classical and in the quantum regime, focusing on two phenomena already addressed in the literature in presence of SB, namely synchronization and entanglement [9, 36]. CB dissipation will be shown to have beneficial effects when compared to the SB case, both on classical synchronization (section 3) and on asymptotic entanglement and optomechanical cooling (section 4).

Zoom In Zoom Out Reset image size

2.1. Hamiltonian model

A simple configuration to assess collective or independent forms of dissipation consists of two optomechanical coupled units, as sketched in figure 1. We allow for a finite mechanical–mechanical coupling to be compared with dissipative coupling effects, while mixing between optical modes is not addressed here for simplicity. The total Hamiltonian is:

$$\begin{eqnarray}&&{\hat{{ \mathcal H }}}_{S}={\hat{{ \mathcal H }}}_{1}+{\hat{{ \mathcal H }}}_{2}+\displaystyle \frac{k}{2}{({\hat{x}}_{1}-{\hat{x}}_{2})}^{2},\end{eqnarray} \tag{ 1 }$$

where ${\hat{{ \mathcal H }}}_{\mathrm{1,2}}$ are the Hamiltonians for each optomechanical oscillator and the last term is the mechanical coupling between them. This kind of mechanical coupling can be found for instance in nanoresonators clamped to the same substrate as in [44, 45, 50], as well as it can be induced optically using off-resonant lasers [51]. Alternatively, interaction with a common (and single) optical field can also induce coupling between mechanical modes of OMs, as in [11, 14]. The Hamiltonian for each unit of the optomechanical system reads as [52]:

$$\begin{eqnarray}&&{\hat{{ \mathcal H }}}_{j}=\displaystyle \frac{1}{2}m{\omega }_{m,j}^{2}{\hat{x}}_{j}^{2}+\displaystyle \frac{1}{2m}{\hat{p}}_{j}^{2}-{\hslash }{{\rm{\Delta }}}_{0}{\hat{a}}_{j}^{\dagger }{\hat{a}}_{j}+{\rm{i}}{\hslash }\sqrt{\displaystyle \frac{\kappa {P}_{{\rm{in}}}}{{\hslash }{\omega }_{L}}}({\hat{a}}_{j}^{\dagger }-{\hat{a}}_{j})-{\hslash }\displaystyle \frac{{\omega }_{c}}{{L}_{{\rm{om}}}}{\hat{x}}_{j}{\hat{a}}_{j}^{\dagger }{\hat{a}}_{j},\ j=1,2\end{eqnarray} \tag{ 2 }$$

in the frame rotating with the input laser frequencies [1]. The first two terms of the Hamiltonian correspond to the mechanical elements, described as harmonic oscillators of frequency ${\omega }_{m,j}$, effective mass m, and position and momentum operators ${\hat{x}}_{j}$, and ${\hat{p}}_{j}$, where $[{\hat{x}}_{j},{\hat{p}}_{j}]={\rm{i}}{\hslash }$. The third and fourth terms describe the driven optical cavities in the frame rotating at the laser frequency, where ${{\rm{\Delta }}}_{0}={\omega }_{L}-{\omega }_{c}$, ${\omega }_{c}$ is the cavity resonance frequency, ${\omega }_{L}$ is the laser frequency, P_in is the laser input power, κ is the cavity energy decay rate, and ${\hat{a}}_{j}^{\dagger }$, ${\hat{a}}_{j}$ are the creation and annihilation operators of the light modes, with $[{\hat{a}}_{j},{\hat{a}}_{j}^{\dagger }]=1$. Notice that we make the simplifying assumption of identical optomechanical oscillators except for the mechanical frequencies, in which allowing disparity is essential for the study of synchronization. Finally the last term in equation (2) describes the optomechanical coupling with a strength parametrized by the effective optomechanical length L_om, whose precise definition is system dependent and might involve the particular properties of the interacting optical and mechanical modes (see for instance its definition for optomechanical crystals [53]). Finally we remark that this model is appropriate for many different OMs as long as the parameters of the system, such as m or L_om, appropriately refer to specific devices (like e.g. examples as diverse as in [53, 54]).

The optical and mechanical environments are described as an infinite collection of harmonic oscillators in thermal equilibrium. The standard model for the environment of an optical cavity, and the master/Langevin equations describing its dynamics, can be found in many textbooks on open quantum systems [38, 55, 56]. We describe the mechanical environments using the model of [57], in which the Born–Markov master equations of two coupled nonidentical harmonic oscillators for the CB/SB cases are derived, starting from the following Hamiltonians. In the SB case each mechanical oscillator is coupled to its own reservoir, so there are two different baths described by:

$$\begin{eqnarray}&&{\hat{{ \mathcal H }}}_{{\rm{bath}},j}^{{\rm{SB}}}=\displaystyle \sum _{\alpha =1}^{\infty }{\hslash }\,{\omega }_{\alpha ,j}{\hat{r}}_{\alpha ,j}^{\dagger }{\hat{r}}_{\alpha ,j},\quad {\hat{{ \mathcal H }}}_{{\rm{int}},j}^{{\rm{SB}}}=\displaystyle \sum _{\alpha =1}^{\infty }\lambda \,{\hat{q}}_{\alpha ,j}{\hat{x}}_{j},\quad j=1,2,\end{eqnarray} \tag{ 3 }$$

where ${\hat{r}}_{\alpha ,j}^{\dagger }$, ${\hat{r}}_{\alpha ,j}$ are the creation and annihilation operators of the bath modes, and ${\hat{q}}_{\alpha ,j}$ their position operators. On the other hand, in the CB case there is only one collective environment for both system units:

$$\begin{eqnarray}&&{\hat{{ \mathcal H }}}_{{\rm{bath}}}^{{\rm{CB}}}=\displaystyle \sum _{\alpha =1}^{\infty }{\hslash }\,{\omega }_{\alpha }{\hat{r}}_{\alpha }^{\dagger }{\hat{r}}_{\alpha },\quad {\hat{{ \mathcal H }}}_{{\rm{int}}}^{{\rm{CB}}}=\displaystyle \sum _{\alpha =1}^{\infty }2\,\lambda \,{\hat{q}}_{\alpha }{\hat{x}}_{+},\end{eqnarray} \tag{ 4 }$$

resulting in coupling through the center of mass coordinate ${\hat{x}}_{+}=({\hat{x}}_{1}+{\hat{x}}_{2})/\sqrt{2}$ [16, 17, 38, 39, 41–44]. Notice that we have included a factor two in equation (4) to enforce an equal energy damping rate for both models: in SB case two independent channels dissipate, in CB case only one coordinate (center of mass) dissipates, but it is coupled to the bath with double strength $2\lambda$. This allows for a quantitative, not only qualitative, comparison of both dissipation regimes.

2.2. Dimensionless parameters and observables

A set of dimensionless observables and parameters can be introduced similarly to [6]. First we define the quantities ${x}_{{\rm{om}}}=\kappa {L}_{{\rm{om}}}/{\omega }_{c}$ and ${n}_{{\rm{\max }}}=4{P}_{{\rm{in}}}/{\hslash }{\omega }_{L}\kappa$, which we use to achieve a dimensionless time and dimensionless operators of the system:

$$\begin{eqnarray}&&{t}^{{\prime} }=\kappa t,\quad {\hat{x}}_{j}^{{\prime} }=\displaystyle \frac{{\hat{x}}_{j}}{{x}_{{\rm{om}}}},\quad {\hat{p}}_{j}^{{\prime} }=\displaystyle \frac{{\hat{p}}_{j}}{m\,\kappa \,{x}_{{\rm{om}}}},\quad {\hat{a}}_{j}^{{\prime} }=\displaystyle \frac{{\hat{a}}_{j}}{\sqrt{{n}_{{\rm{\max }}}}},\end{eqnarray} \tag{ 5 }$$

where the dimensionless quantities are denoted by a prime. Then the following dimensionless parameters appear in the equations of motion:

$$\begin{eqnarray}&&{\omega }_{m,j}^{{\prime} }=\displaystyle \frac{{\omega }_{m,j}}{\kappa },\ {{\rm{\Gamma }}}^{{\prime} }=\displaystyle \frac{\gamma }{\kappa },\ {K}_{c}^{{\prime} }=\displaystyle \frac{k}{m{\kappa }^{2}},\ {{\rm{\Delta }}}_{0}^{{\prime} }=\displaystyle \frac{{{\rm{\Delta }}}_{0}}{\kappa },\ {{ \mathcal P }}^{{\prime} }=\displaystyle \frac{4{P}_{{\rm{in}}}{\omega }_{c}}{{{mL}}_{{\rm{om}}}^{2}{\kappa }^{4}}.\end{eqnarray} \tag{ 6 }$$

Unless stated otherwise we will work with these dimensionless parameters in the following, and thus we can drop the primes.

2.3. Nonlinear classical equations

After defining the set of dimensionless parameters and observables used in this work, we can write down the dimensionless equations describing the classical dynamics of this system [9], both for SB and CB. The classical dynamics can be obtained by the quantum equations of motion taking the first moments and making the approximation of factorizing expectation values of nonlinear terms, as usual. The quantum dynamics can be obtained either from a master equation approach as well as from the quantum Langevin equations of the system (see next section) [1, 56, 57]:

$$\begin{eqnarray}&&{\dot{a}}_{j}=\left[{\rm{i}}({{\rm{\Delta }}}_{0}+{x}_{j})-\displaystyle \frac{1}{2}\right]{a}_{j}+\displaystyle \frac{1}{2},\quad j=1,2,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\dot{{x}}}_{{j}}={{p}}_{{j}},\,{\dot{{p}}}_{{j}}=-{{\omega }}_{{m},{j}}^{2}{{x}}_{{j}}+{(-1)}^{{j}}{{K}}_{{c}}({{x}}_{1}-{{x}}_{2})-{{ \mathcal D }}_{{\rm{SB}}/{\rm{CB}}}+{ \mathcal P }| {{a}}_{{j}}{| }^{2},\end{eqnarray} \tag{ 8 }$$

where variables without hat denote expectation values, e.g., $x=\langle \hat{x}\rangle$. The mechanical dissipation terms are defined as ${{ \mathcal D }}_{{\rm{SB}}}={\rm{\Gamma }}{p}_{j}$ for the SB case [9], and ${{ \mathcal D }}_{{\rm{CB}}}={\rm{\Gamma }}({p}_{1}+{p}_{2})$ for the CB case [57].

The dynamical behavior of OMs can be very diverse depending on the parameter choice [1]. When the mechanical resonator is 'too slow' with respect to the cavity lifetime (low mechanical quality factor, and ${\omega }_{m}$ smaller than κ, in dimensional units), the intracavity power follows adiabatically changes in the mechanical displacement, and a fixed point in which radiation pressure equilibrates with the mechanical restoring force is reached [2]. Otherwise, when the mechanical resonator is able to follow the fast cavity dynamics (high quality factors, and ${\omega }_{m}$ comparable to κ), dynamical backaction enables striking phenomena such as damping (cooling) or anti-damping (heating) of the mechanical motion by the optical cavity [1, 58]. In fact optomechanical cooling is enhanced in the red-sideband regime (${{\rm{\Delta }}}_{0}\lt 0$) whereas anti-damping in the blue-sideband regime (${{\rm{\Delta }}}_{0}\gt 0$). For enough laser power, ${ \mathcal P }$, anti-damping can overcome intrinsic mechanical damping resulting in a dynamical instability and self-sustained oscillations [1, 58]. In these conditions the OM experiences a supercritical Hopf bifurcation when varying continuously ${{\rm{\Delta }}}_{0}$ from negative to positive values [6, 59]. Moreover, for strong driving a dynamical multistability is found, and self-sustained oscillations are possible even in the red-sideband regime [6]. The regime of self-sustained oscillations has been explored in the context of spontaneous synchronization in presence of independent dissipation for two coupled mechanical units [9] and in section 3 we are going to address the effect of dissipative coupling (CB).

2.4. Linear quantum Langevin equations

From the above Hamiltonian model and the input–output formalism we can obtain a set of Langevin equations describing the noise driven damped dynamics of the operators [1] (see also the appendix), that will be used in section 4 to explore a possible improvement of non-classical effects in presence of collective dissipation. We are going to consider the system in a stable and stationary state where a linearized treatment is a suitable approximation for the fluctuations dynamics (for more details see for instance [1, 34, 60]). Low temperatures and high mechanical and optical quality factors [1] also contribute to maintain low levels of noise (fluctuations). Indeed setting the OMs in a stable fixed point [34, 60] the fluctuation operators can be defined as: $\delta \hat{O}=\langle O{\rangle }_{{\rm{st}}}-\hat{O}$, where $\langle O{\rangle }_{{\rm{st}}}$ is the constant solution of equations (7) and (8), and the linear equations for the fluctuation operators are:

$$\begin{eqnarray}&&\delta {\dot{\hat{Q}}}_{j}=-(\langle {x}_{j}{\rangle }_{{\rm{st}}}+{{\rm{\Delta }}}_{0})\delta {\hat{P}}_{j}-\langle {P}_{j}{\rangle }_{{\rm{st}}}\,\delta {\hat{x}}_{j}-\displaystyle \frac{1}{2}\delta {\hat{Q}}_{j}+\sqrt{\displaystyle \frac{{\omega }_{m}}{{ \mathcal P }}}\,{\hat{Q}}_{{\rm{in}},j},\quad j=1,2,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\delta {\dot{\hat{P}}}_{j}=(\langle {x}_{j}{\rangle }_{{\rm{st}}}+{{\rm{\Delta }}}_{0})\delta {\hat{Q}}_{j}+\langle {Q}_{j}{\rangle }_{{\rm{st}}}\,\delta {\hat{x}}_{j}-\displaystyle \frac{1}{2}\delta {\hat{P}}_{j}+\sqrt{\displaystyle \frac{{\omega }_{m}}{{ \mathcal P }}}\,{\hat{P}}_{{\rm{in}},j},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\delta {\dot{\hat{x}}}_{j}=\delta {\hat{p}}_{j}-{d}_{{\rm{SB}}/{\rm{CB}}}^{x}+{\hat{x}}_{{\rm{in}},j},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\delta {\dot{\hat{p}}}_{j}=-{\omega }_{m}^{2}\delta {\hat{x}}_{j}+{(-1)}^{j}{K}_{c}(\delta {\hat{x}}_{1}-\delta {\hat{x}}_{2})-{d}_{{\rm{SB}}/{\rm{CB}}}^{p}+{ \mathcal P }(\langle {Q}_{j}{\rangle }_{{\rm{st}}}\,\delta {\hat{Q}}_{j}+\langle {P}_{j}{\rangle }_{{\rm{st}}}\,\delta {\hat{P}}_{j})+{\omega }_{m}\,{\hat{p}}_{{\rm{in}},j}.\end{eqnarray} \tag{ 12 }$$

In the following we are going to consider the light quadratures ${\hat{Q}}_{j}=({\hat{a}}_{j}+{\hat{a}}_{j}^{\dagger })/\sqrt{2}$ and ${\hat{P}}_{j}={\rm{i}}({\hat{a}}_{j}^{\dagger }-{\hat{a}}_{j})/\sqrt{2}$. We also notice that we are considering identical optomechanical oscillators in the study of quantum correlations (section 4). As usual, the fluctuation operators are rescaled by the parameter $\zeta ={\tilde{x}}_{{zpf}}/{x}_{{\rm{o}}{\rm{m}}}$ with ${\tilde{x}}_{{zpf}}=\sqrt{{\hslash }/m{\omega }_{m}}$, so that $\delta {\hat{O}}_{{\rm{new}}}=\delta {\hat{O}}_{{\rm{old}}}/\zeta$.

The mechanical dissipation terms are:

$$\begin{eqnarray}&&{d}_{{\rm{SB}}}^{x}={\rm{\Gamma }}\delta {\hat{x}}_{j},\quad {d}_{{\rm{SB}}}^{p}={\rm{\Gamma }}\delta {\hat{p}}_{j},\quad {d}_{{\rm{CB}}}^{x}={\rm{\Gamma }}(\delta {\hat{x}}_{1}+\delta {\hat{x}}_{2}),\quad {d}_{{\rm{CB}}}^{p}={\rm{\Gamma }}(\delta {\hat{p}}_{1}+\delta {\hat{p}}_{2}),\end{eqnarray} \tag{ 13 }$$

where the known equivalence in position and momentum damping results from the rotating wave approximation in the system-bath coupling, a valid approximation in the limit ${\rm{\Gamma }}\ll {\omega }_{m}$ (i.e. high mechanical quality factors) [38]. In the Markovian limit, the zero mean Gaussian noise terms of equations (9)–(12) are characterized by the following symmetrized correlations [34]:

$$\begin{eqnarray}&&\langle \{{\hat{Q}}_{{\rm{in}},j}(t),{\hat{Q}}_{{\rm{in}},j}(t^{\prime} )\}\rangle =\langle \{{\hat{P}}_{{\rm{in}},j}(t),{\hat{P}}_{{\rm{in}},j}(t^{\prime} )\}\rangle =\delta (t-t^{\prime} ),\quad j=1,2,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\langle \{{\hat{x}}_{{\rm{in}},j}(t),{\hat{x}}_{{\rm{in}},j}(t^{\prime} )\}\rangle =\displaystyle \frac{1}{2}\langle \{{\hat{p}}_{{\rm{in}},j}(t),{\hat{p}}_{{\rm{in}},j}(t^{\prime} )\}\rangle ={\rm{\Gamma }}(2{n}_{{\rm{th}}}+1)\delta (t-t^{\prime} ),\end{eqnarray} \tag{ 15 }$$

where $\{...,\ldots \}$ denotes anticommutator, $\langle ...\rangle$ denotes an ensemble average, $\delta (t-t^{\prime} )$ is the Dirac delta, and ${n}_{{\rm{th}}}={(\exp [{\hslash }{\omega }_{m}/{k}_{B}T]-1)}^{-1}$ is the phonon occupancy number of the mechanical baths, assumed to be at the same temperature T. The light noise correlations for the corresponding quadratures are obtained by considering optical environments at zero temperature for which the only nonvanishing correlations are $\langle {\hat{a}}_{{\rm{in}},j}(t){\hat{a}}_{{\rm{in}},j}^{\dagger }(t^{\prime} )\rangle =\delta (t-t^{\prime} )$ for j = 1, 2 [34]. In the CB case additional cross-correlations appear:

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\langle \{{\hat{x}}_{{\rm{in}},1}(t),{\hat{x}}_{{\rm{in}},2}(t^{\prime} )\}\rangle =\displaystyle \frac{1}{2}\langle \{{\hat{p}}_{{\rm{in}},1}(t),{\hat{p}}_{{\rm{in}},2}(t^{\prime} )\}\rangle ={\rm{\Gamma }}(2{n}_{{\rm{th}}}+1)\delta (t-t^{\prime} ).\end{eqnarray} \tag{ 16 }$$

Finally we remark that as the fluctuations dynamics is linear and the noise is Gaussian, states initially Gaussian will remain Gaussian at all times [61, 62]. Notably Gaussian states are completely characterized by their first and second moments, which are all encoded in the covariance matrix of the system [61, 62]. Rewriting the Langevin equations (9)–(12) as $\dot{\vec{R}}={ \mathcal M }\vec{R}+\vec{D}$, where $\vec{R}$ is the vector of the fluctuations of position, momentum, and light quadratures of the system, $\vec{D}$ is the vector containing the zero mean Gaussian noise terms, and ${ \mathcal M }$ is the matrix generating the dynamics of the system, the following equation for the covariance matrix can be written down:

$$\begin{eqnarray}&&\dot{{ \mathcal C }}={ \mathcal M }\,{ \mathcal C }+{ \mathcal C }\,{{ \mathcal M }}^{T}+{ \mathcal N },\end{eqnarray} \tag{ 17 }$$

with the covariance matrix and the noise covariance matrix defined respectively as:

$$\begin{eqnarray}&&{{ \mathcal C }}_{{ij}}(t)=\displaystyle \frac{\langle {R}_{i}(t){R}_{j}(t)+{R}_{j}(t){R}_{i}(t)\rangle }{2},\quad {{ \mathcal N }}_{{ij}}(t)=\displaystyle \frac{\langle {D}_{i}(t){D}_{j}(t)+{D}_{j}(t){D}_{i}(t)\rangle }{2},\end{eqnarray} \tag{ 18 }$$

with $i,j=1,\,\ldots ,\,8$. By solving equation (17) at any time we can completely characterize the quantum correlations present in the system [63, 64].

In this section we analyze the effects of collective dissipation on the synchronization of the first moments of the system, or classical synchronization. In the SB case classical synchronization has been already studied in [9], where it has been shown that the OMs described by equations (7) and (8) can synchronize with a locked phase difference either tending to 0, in-phase synchronization, or to π, anti-phase synchronization. To study the presence of synchronization we integrate numerically equations (7) and (8) in both CB/SB cases and in a parameter region in which the system is self-oscillating. Once the system is in the steady state we compute a synchronization measure. Synchronization is characterized by a Pearson correlation function defined as:

$$\begin{eqnarray}&&{C}_{{x}_{1},{x}_{2}}(t,{\rm{\Delta }}t)=\displaystyle \frac{\overline{\delta {x}_{1}\delta {x}_{2}}}{\sqrt{\overline{\delta {x}_{1}^{2}}\,\overline{\delta {x}_{2}^{2}}}},\end{eqnarray} \tag{ 19 }$$

where the bar represents a time average with time window ${\rm{\Delta }}t$, i.e. $\overline{x}=\tfrac{1}{{\rm{\Delta }}t}{\int }_{t}^{t+{\rm{\Delta }}t}{\rm{d}}{s}\,x(s)$, and $\delta x=x(t)-\overline{x}$. This correlation function has proven to be a useful indicator of synchronization both in classical and quantum systems [15, 65]. In particular this measure provides an absolute scale for synchronization strength as it takes values between –1 and 1, where –1 means perfect anti-phase synchronization and 1 perfect in-phase synchronization. Furthermore this measure can be generalized also in presence of time delays accounting for regimes where OMs synchronize with a phase difference depending on the parameter values. In this case two delayed functions, ${x}_{1}(t)$ and ${x}_{2}(t+\tau )$, need to be considered in equation (19). The best figure of merit accounting for delay is found maximizing the synchronization measure for different values of τ: in particular we take an interval of time ${\rm{\Delta }}t$ such that it includes many periodic oscillations, and we compute the correlation function (19) between ${x}_{1}(t)$ and ${x}_{2}(t+\tau )$, where τ is varied between $[0,{ \mathcal T }]$, being ${ \mathcal T }$ the period of the oscillation. By keeping the maximum correlation C_max, and the temporal shift at which it occurs ${\tau }_{{\rm{\max }}}$, we obtain a measure of the degree or quality of synchronization, C_max, and the locked phase difference, ${\tau }_{{\rm{\max }}}/{ \mathcal T }$, in units of $2\pi$.

The degree of synchronization accounting for delay, C_max, as a function of the frequency detuning, ${\rm{\Delta }}{\omega }_{m}={\omega }_{m,2}-{\omega }_{m,1}$, and the mechanical coupling strength, K_c, is shown in figure 2, where synchronization is found for coupling overcoming detuning between the mechanical oscillators. The case of dissipation in a CB is represented also including a comparison with the SB case studied in [9] in which a similar synchronization diagram was already observed. In the synchronized regions, ${C}_{{\rm{\max }}}\approx 1$, the mechanical elements oscillate at the same frequency describing sinusoidal trajectories with constant amplitude. Conversely, when there is no synchronization ($0\lesssim {C}_{{\rm{\max }}}\lesssim 0.6$ in the represented case), the mechanical oscillations have different frequencies and display strong amplitude modulations.

Zoom In Zoom Out Reset image size

In presence of synchronization two further regimes are recognized depending on whether the locked phase difference tends to 0 or to π. At threshold the system first anti-synchronize and by further increasing the coupling for a given detuning there is a transition to in-phase synchronization (respectively, dashed and dashed-dotted lines of figure 2, for both CB (salmon lines) and SB (blue lines)). The comparison allows to identify how collective dissipation favors spontaneous synchronization: indeed, the synchronized area for CB is much larger than for SB, as the first anti-phase synchronization threshold coupling is lowered (salmon and blue dashed 'π-sync.' lines). The anti-phase synchronization regime is also enlarged for CB case as the second threshold to in-phase synchronization increases for CB with respect to SB. The fact the anti-phase synchronization is favored is consistent with the particular dissipation form in the CB case: it is the common coordinate ${x}_{1}+{x}_{2}$ that is damped by the environment, see equation (8), leading to a 'preferred' amplitude of motion in ${x}_{1}-{x}_{2}$, and thus an anti-phase locking.

The effect of dissipation on the synchronization threshold is quantified for CB and SB in figure 3, for fixed detuning, looking at the mechanical coupling needed for the system to synchronize when increasing the damping rate Γ. This mechanical coupling threshold, K_c^thres, is defined as the minimum amount of mechanical coupling, K_c, needed for the oscillators to synchronize, and can be easily obtained as the synchronization of the system is signaled by a sharp transition of the synchronization indicator from small values (for our parameters below ${C}_{{\rm{\max }}}\approx 0.6$) to the values around unity that characterize synchronized oscillations. The damping strength does not influence significantly the synchronization threshold in SB case, in stark contrast with the CB case, which is very sensitive to the damping strength and improving significantly when increasing Γ. The enlargement in the (anti-)synchronization area for the CB case is also displayed in the full parameter plots (figure 3(b)). This clearly shows how the presence of an additional (dissipative) coupling term between the oscillators in the CB case triggers the emergence of spontaneous synchronization. In other words, both reactive and dissipative couplings contribute similarly in overcoming detuning effects in the context of synchronization. Furthermore CB specially favors anti-phase synchronization being the relative motion of the oscillators shielded by damping effects.

Zoom In Zoom Out Reset image size

An interesting effect of collective dissipation is that actually spontaneous synchronization can emerge even in absence of direct mechanical coupling (see figure 2 for ${\rm{\Delta }}{\omega }_{m}\lesssim 0.01$ and figure 3(a) for ${\rm{\Gamma }}\gtrsim 0.05$). This is a specific signature of CB dissipation, being indeed not possible between uncoupled OMs with SB, and shows the constructive role played by the dissipative coupling due to CB even in absence of any other optical or mechanical coherent coupling. The threshold for synchronization between decoupled oscillators (K_c = 0) in the CB case is indeed set by the damping strength occurring for Γ overcoming detuning effects (${\rm{\Delta }}{\omega }_{m}$) as nicely shown in the inset of figure 3(a). We also notice that, for Γ small enough, synchronization is always in-phase (red points) while when increasing it different phase (delays) can be favored at threshold (figure 3). In particular, in the CB case, the phase-locked difference does not always correspond to in-phase or anti-phase synchronization, as it also takes intermediate values $\sim \pi /2$. This occurs both with (figure 3(a)) and without (inset in figures 3(a) and (b)) direct mechanical coupling, being more evident for larger dissipation rates (figures 3(a) and (b)). The exact mechanism behind these $\pi /2$ values of the locked phase-difference is an open question.

In this section we study the entanglement between mechanical modes in the asymptotic state of the dynamics. As explained in section 2 we set the system in a stable fixed point and we study the linear fluctuations around the mean state. The stationary state of the fluctuations is independent of the initial conditions and we characterize it obtaining numerically the stationary covariance matrix from equation (17). From this covariance matrix we can compute entanglement between the different degrees of freedom [61, 62], for which we use the logarithmic negativity, a valid entanglement monotone for mixed Gaussian states [61, 62]. Its definition for modes described by quadratures with commutation relations $[{q}_{\mathrm{1,2}},{p}_{\mathrm{1,2}}]=i$ is ${E}_{N}=\max \{0,-{\rm{ln}}2{\tilde{\nu }}_{-}\}$, where ${\tilde{\nu }}_{-}$ is the smallest symplectic eigenvalue of the partially transposed covariance matrix of the subsystem [62].

One of the main results discussed below is the tight relation between the effective temperature of the mechanical modes and the amount of asymptotic entanglement that is preserved between them: the linear regime in which the system is operated basically consists of four coupled units which in the stationary state reach a given effective temperature depending on the effectivity of the optomechanical cooling. It seems that this effective temperature regulates the entanglement behavior, which is that of a coupled thermal system. Furthermore we find that the CB case has higher mechanical entanglement and is easier to cool down than the SB case, specially when increasing the relative mechanical coupling strength ${K}_{c}/{\omega }_{m}^{2}$. We analyze next the dependence of cooling and entanglement on the different system parameters, showing that this is a puzzle with several pieces. At this point we remark that as in the following we are dealing with identical optomechanical oscillators, quantities such as the occupation number of each of the mechanical oscillators (${n}_{\mathrm{eff}}$), or the optical-mechanical entanglement between a mechanical mode and its optical mode, are the same for the two optomechanical oscillators.

Mechanical coupling K_c. We first consider asymptotic entanglement when varying the mechanical coupling between the oscillators. A minimum coupling is needed to attain non-vanishing entanglement between the mechanical modes (solid lines in figure 4(a)), both for CB and SB. Indeed, only above some minimum mechanical coupling entanglement overcomes heating effects as shown (through the mechanical occupancy number ${n}_{\mathrm{eff}}$) in figure 4(b). This is known to occur also in the simplest configuration of a pair of coupled harmonic oscillators at finite temperature [66, 67]. We mention that actually for similar (effective) temperatures and similar mechanical coupling strengths we find that both harmonic and OMs display similar amounts of entanglement. Slightly above this threshold the build-up of entanglement is similar in presence of CB or SB (figure 4(a)), as also the heating of the mechanical components (figure 4(b)), but further increasing the coupling differences arise and entanglement worsen for SB. This effect is accompanied by a stronger heating in the SB than in the CB case (figure 4(b)). Overall, larger entanglement, as well as better cooling (smaller ${n}_{\mathrm{eff}}$), are found for the CB case.

Zoom In Zoom Out Reset image size

From these observations we infer that the mechanical coupling between the oscillators influences their entanglement through two competing mechanisms: on the one hand it is the coupling that entangles the two mechanical oscillators, on the other K_c modifies the effective temperature of the oscillators. This is because the normal mode frequencies change with K_c reducing their resonance with the optical modes and therefore hindering their cooling (see discussion below and figure 5). As the heating produced increasing K_c is more pronounced for independent dissipation, best entanglement values are attained in the CB case. This picture is confirmed when looking at light-mirror entanglement: K_c is not expected to contribute significantly to the light-mirror correlations in each OM oscillator, while the effects of heating will still be present. As we see from figure 4(a) (dashed-dotted lines) the optical-mechanical entanglement is present for all couplings (no threshold) and diminishes when increasing K_c. Both for CB and SB cases the loss of entanglement is similar to the effective heating (increasing effective temperature) in accordance with the above considerations.

Zoom In Zoom Out Reset image size

Cavity-laser detuning ${{\rm{\Delta }}}_{0}$. The detuning ${{\rm{\Delta }}}_{0}$, set to values comparable to the mechanical frequencies, allows for quanta exchange between the very disparate optical and mechanical modes. In this way the optical reservoir, effectively at T = 0, acts as a heat sink for the mechanical degrees of freedom when the system operates in the red-sideband regime [27, 28]. Thus, resonance between optical and mechanical degrees of freedom controls the effectiveness of the cooling process. Looking at figures 4(c) and (d) we see that the entanglement between oscillators and their effective temperature display a strong inverse behavior, i.e. as the effective temperature drops the entanglement grows and vice versa, and the maximum of entanglement is very close to the minimum of temperature for both CB/SB cases. This behavior is found also for different temperatures and mechanical couplings. Despite the fact that the optimal ${{\rm{\Delta }}}_{0}$ values for cooling and entanglement are slightly (less than 2%) shifted one with respect to the other, the strong anticorrelations between the two curves and the small size of these shifts, make us conclude that the laser detuning ${{\rm{\Delta }}}_{0}$ modifies the asymptotic mechanical entanglement mainly through the degree of optomechanical cooling. We also note that the detuning for optimal cooling in both CB/SB cases is displaced with respect to known values for single OM ${{\rm{\Delta }}}_{0}\approx -{\omega }_{m}$ [27, 28] (see discussion on optimal detuning). Finally we note that entanglement and cooling are again enhanced in the CB case, and that the optimal detunings for maximum entanglement are quite different in the CB and SB case (the same applies to the detunings for optimal cooling).

Laser power ${ \mathcal P }$. Similarly to the case of ${{\rm{\Delta }}}_{0}$, ${ \mathcal P }$ modifies the mechanical entanglement mainly by modifying the degree of optomechanical cooling as it follows from the strong inverse behavior of figures 4(e) and (f). Recall from the linear quantum Langevin equations, that ${ \mathcal P }$ has the role of a coupling strength between the optical and mechanical quantum fluctuations. It is then to be expected that optomechanical cooling improves with ${ \mathcal P }$. However, when ${ \mathcal P }$ becomes very high, it is known [28, 68] that optical and mechanical modes begin to hybridize and heating from quantum backaction noise enters the scene. Thus at first higher ${ \mathcal P }$ implies an enhancement of the optomechanical cooling effect, since this coupling parametrizes the strength of the optomechanical cooling rate [27, 28]. As ${ \mathcal P }$ is further increased, the regime of strong coupling is reached and heating of the mechanical modes by radiation pressure becomes important. The competition between these two effects leads to a minimum of ${n}_{\mathrm{eff}}$ as a function of the optomechanical coupling strength, as derived for a single OM in [28, 68] and observed in figure 4(f). Finally, the CB case requires less power to achieve the same cooling as compared to the SB case.

Optimal detuning. As we have observed in the discussion of figures 4(c) and (d), the asymptotic mechanical entanglement changes with ${{\rm{\Delta }}}_{0}$ mainly because of the effectiveness of the optomechanical cooling process. The detuning optimizing each process indeed takes very similar values as mentioned above. In figure 5 we analyze in more detail the detunings ${{\rm{\Delta }}}_{0}$ which are optimal for optomechanical cooling and mechanical entanglement as a function of the relative mechanical coupling strength ${K}_{c}/{\omega }_{m}^{2}$, and address some of the differences between the CB/SB cases that have been pointed out in the above discussions. In figure 5(a) we can see that as ${K}_{c}/{\omega }_{m}^{2}$ increases, the ${{\rm{\Delta }}}_{0}$ maximizing entanglement is shifted towards more negative values. In figure 5(b) we have plotted the maximum entanglement for the optimal value of ${{\rm{\Delta }}}_{0}$ as a function of ${K}_{c}/{\omega }_{m}^{2}$, and we see that it increases in both cases and it is larger for the CB case. We note that the increase with ${K}_{c}/{\omega }_{m}^{2}$ is similar to that in 4(a), but that here we compute the value at the optimal ${{\rm{\Delta }}}_{0}$ for each ${K}_{c}/{\omega }_{m}^{2}$, ${{\rm{\max }}}_{{{\rm{\Delta }}}_{0}}{E}_{N}({K}_{c}/{\omega }_{m}^{2})$, thus increasing a bit the level of achieved entanglement. Analogous results are obtained when studying the detuning for minimum effective temperature (not shown here).

In figure 5(c) we address the question of how the detunings that maximize entanglement or minimize the effective temperature vary with the relative mechanical coupling strength. In particular we plot the ${{\rm{\Delta }}}_{0}$ that optimizes cooling (solid lines) and the one that optimizes entanglement (dashed lines) as a function of ${K}_{c}/{\omega }_{m}^{2}$. The solid black line corresponds to ${{\rm{\Delta }}}_{0}=-{{\rm{\Omega }}}_{+}$, the dashed-dotted line to ${{\rm{\Delta }}}_{0}=-{{\rm{\Omega }}}_{-}$, and the cross-line to ${{\rm{\Delta }}}_{0}=-\overline{{\rm{\Omega }}}=-({{\rm{\Omega }}}_{+}+{{\rm{\Omega }}}_{-})/2$. Where ${{\rm{\Omega }}}_{\pm }$ are the frequencies of the normal modes of the isolated mechanical system ${x}_{\pm }=({x}_{1}\pm {x}_{2})/\sqrt{2}$ (i.e. without optomechanical coupling), and which take the values ${{\rm{\Omega }}}_{+}={\omega }_{m}$ and ${{\rm{\Omega }}}_{-}=\sqrt{{\omega }_{m}^{2}+2{K}_{c}}$. From this figure it is appreciated that the shift between the optimal detunings for cooling and entanglement increases slightly with the mechanical coupling, but it remains small in all the range both for CB and SB. We also observe that in the SB case, the optimal strategy for cooling is to set (${{\rm{\Delta }}}_{0}\approx -\overline{{\rm{\Omega }}}$), and thus the minimum mechanical effective temperature is achieved when both mechanical normal modes (of the isolated system) are cooled at the same rate. Notice that given the small difference between the detunings optimizing each quantity, the same strategy applies, to a good approximation, to obtain the maximum entanglement. On the other hand, in the CB case the optimal ${{\rm{\Delta }}}_{0}$ for cooling and entanglement is shifted towards $-{{\rm{\Omega }}}_{+}$, consistent with the fact that dissipation enters through this coordinate. Why is it then not the case that the optimal detuning coincides with ${{\rm{\Omega }}}_{+}$? Well, the common coordinate ${x}_{+}$ is a normal mode of the isolated mechanical system, but it is not of the full OM, so ${x}_{+}$ and ${x}_{-}$ are coupled through the optical degrees of freedom. This is why it is necessary to cool both coordinates, although it is more important to cool the one that dissipates most (${x}_{+}$). In fact, as we increase ${K}_{c}/{\omega }_{m}^{2}$, the role of the optomechanical coupling is diminished, moving the optimal detuning towards ${{\rm{\Omega }}}_{+}$. This also explains why the SB case is harder to cool than the CB one: in SB we need to address both normal modes ${x}_{\pm }$ with a single laser frequency, whereas for CB it is better to address only ${x}_{+}$.

Resolved sideband. In the following we explore mechanical entanglement (figure 6(a)) and optomechanical cooling (figure 6(b)) as we go deeper into the resolved sideband regime ${\omega }_{m}\gt 1$, where the parameter ${\omega }_{m}$ measures how well the linewidth of the optical modes resolves the frequency of the mechanical ones. It is known that the efficiency of the cooling process of a single OM increases initially as ${\omega }_{m}$ increases [27, 28, 68] (recall that ${\omega }_{m}/\kappa \mapsto {\omega }_{m}$ in our dimensionless units, see section 2.2), and thus lower phonon numbers can be achieved when the intracavity power is adjusted optimally [68]. It is then interesting to analyze whether increasing ${\omega }_{m}$ can enhance cooling and entanglement in the CB/SB cases too. The maximal mechanical entanglement and the minimal mechanical occupancy number optimized in terms of ${ \mathcal P }$ for each value of ${\omega }_{m}$ are displayed in figure 6.

Zoom In Zoom Out Reset image size

In figures 6(a), (b), we can see how for the CB case the maximal entanglement (minimal occupancy number) increases (decreases) as ${\omega }_{m}$ is increased. The SB case at first behaves as the CB case, but for larger ${\omega }_{m}$ maximal entanglement begins to diminish and the phonon number increases. Notice that again the asymptotic mechanical entanglement and the effective mechanical temperature are strongly anticorrelated following an inverse behavior in both CB/SB cases. The increase of the minimum effective temperature in the SB case is more evident in the high temperature case (dashed lines figures 6(a) and (b)), and occurs even at optimal detuning ${{\rm{\Delta }}}_{0}=-\overline{{\rm{\Omega }}}$ (not shown here). The difference between CB and SB cases can be explained as follows: for CB, cooling the coordinate ${x}_{+}$ (by setting ${{\rm{\Delta }}}_{0}=-{\omega }_{m}=-{{\rm{\Omega }}}_{+}$) is almost optimal, and the more we resolve this frequency with our optical modes the better. For SB however, we need to cool both ${x}_{+}$ and ${x}_{-}$, but this has to be done with a single optical frequency: the result is that the better we address ${{\rm{\Omega }}}_{+}$, the worse we address ${{\rm{\Omega }}}_{-}$ which gets heated up. Thus in a SB situation it will be more appropriate to have one optical mode targeted at ${{\rm{\Delta }}}_{1}=-{{\rm{\Omega }}}_{+}$, and another at ${{\rm{\Delta }}}_{2}=-{{\rm{\Omega }}}_{-}$, provided we can resolve the difference between ${{\rm{\Omega }}}_{+}$ and ${{\rm{\Omega }}}_{-}$ which requires high mechanical couplings or very low κ.

In order to point out further differences between the CB/SB cases it is also interesting to study the relative linear optomechanical coupling, $g/{\omega }_{m}$ (with $g=({\omega }_{c}{\tilde{x}}_{{zpf}}/{L}_{{\rm{om}}})\sqrt{0.5\langle {a}^{\dagger }a{\rangle }_{{\rm{st}}}}$ [1], in dimensional units), at which optimal cooling is achieved (inset figure 6(b)). From this figure we can observe that as we operate the system deeper in the resolved sideband regime (larger ${\omega }_{m}$) this optomechanical coupling for optimal cooling diminishes in the CB case. This is so because the exchange between optical and mechanical modes is more efficient, and it is the case both for one OM (yellow solid lines and blue cross-lines) [68], and for the CB case (where addressing ${x}_{+}$ matters most). In fact, if we increase it too much, we end up heating the mechanical modes due to radiation pressure as shown in figure 4(f) and in [28, 68]. Peculiarly, in the SB case, the optimal $g/{\omega }_{m}$ increases monotonically with ${\omega }_{m}$. This can be interpreted as follows: as before, the more we resolve, the worse we cool one of the coordinates ${x}_{\pm };$ the interesting thing is that this effect is initially counteracted by increasing $g/{\omega }_{m}$ (see small values of ${\omega }_{m}$, figure 6(b)), which increases the range of frequency that the mechanical mode is able to exchange energy with the optical modes. This is just the well-known fact that stronger couplings allow for resonance with further detuned frequencies, which here counteracts the decreasing optical linewidth. However, as a side effect, radiation pressure heating also grows with the optomechanical coupling, and thus the maximum achievable entanglement (degree of cooling) eventually decreases as shown in figures 6(a) and (b).

For a more quantitative discussion of cooling we simulate the case of a single OM (yellow solid lines and blue cross-lines, figure 6(b) and inset), modeled by equations (9)–(12) for one of the two OMs in the SB case and with K_c = 0. The case of a single OM is studied in [68], where the optimization of the cooling process in terms of the linear optomechanical coupling is addressed. This is equivalent to what it is shown in figure 6(b) and in the inset, and we plot their results (see last equations of [68]) in blue cross-lines, using the parameters of the low temperature case. We have also obtained the same quantities numerically (yellow solid lines) finding good agreement between the analytical expressions of [68] and the exact results. Comparing the single OM case and the CB case, we find that they follow very similar behaviors both for the minimal occupancy number as well as for the optimal $g/{\omega }_{m}$ for cooling. In particular, lower phonon numbers are achieved in the single OM case, while the optimal $g/{\omega }_{m}$ for cooling takes very close values in both cases. On the contrary, the behavior displayed in the SB case is qualitatively different as discussed above. The similitudes between the single OM case and the CB case reinforce the the overall idea that in the CB case we mainly need to address only one dissipative mode, i.e. ${x}_{+}$.

Miniaturization and dense packing of optomechanical devices can give rise to collective effects in dissipation [17], what we have called here CB, in contrast to the widely used SB dissipation. We have explored the consequences of this form of dissipation in two interesting facets of these kind of nonlinear systems: classical synchronization and steady-state quantum correlations. We have found that the parameter range where synchronization occurs is notably enlarged for CB, and also that in the CB case more dissipation can be beneficial, leading to synchronization even in the absence of mechanical coupling as both the dissipative coupling and the coherent one lead to synchronization between detuned mechanical elements.

With respect to steady-state entanglement of mechanical motion, we have found that the CB case requires less input laser power, because it is more effectively used for cooling, yielding higher levels of entanglement. This can be understood when considering that in this case it is possible to address only the ${x}_{+}$ mechanical eigenmode, which is the only one dissipating, whereas for SB both modes dissipate equally. This easy picture is quantitatively shown to hold exploring different parameter regions. Furthermore we have shown that going to the resolved sideband regime, the CB case gets monotonically better entanglement and cooling, whereas in SB it improves only up to a point, after which it begins to worsen. This is also due to the inability to cool both eigenmodes with a perfectly resolved sideband at only one of the eigenfrequencies, for which we also propose a two-tone cooling scheme as solution.

As stated in the introduction, a collective dissipation channel is found to be the dominant one in the optomechanical setup of references [44, 45], consisting on the emission of elastic radiation by the common mechanical mode ${x}_{+}$, and thus being analogous to our CB model. Moreover the fact that this same mechanism is observed in other mechanical oscillators with a similar geometry [48, 49], indicates that collective dissipation is relevant in closely clamped mechanical resonators. Despite the mechanical coupling strength in this device [44, 45] is an order of magnitude smaller than the regime studied here, it is possible to increase it up to the analyzed values by reducing the separation of the beams' clamping points and increasing the size of the overhang or common substrate, as shown in [69]. Interestingly this last point suggests that strong mechanical coupling and collective dissipation might go hand in hand. On the other hand, a complementary approach consists on inducing the mechanical coupling by strongly driving the system with off-resonant lasers [51].

In this work we have chosen to study the particular collective dissipation channel ${x}_{+}$ motivated by the numerous previous works studying this particular type of dissipation emerging from microscopic physical models in different contexts, starting from works on superradiance [41]. Furthermore this collective dissipation has actually been reported in some mechanical platforms as discussed in the introduction [44, 45, 48, 49]. However, other forms of collective dissipation are also possible: ${x}_{-}$ dissipation can arise in microscopic lattice models [17], and it has been reported in certain systems of trapped ions [70], while more complex forms of collective jumps are considered in [40]. Moreover, quantum correlations and quantum synchronization for ${x}_{-}$ CB dissipation have been studied both for harmonic oscillators [16, 57], for two level systems [43], and quantum Van der Pol oscillators [20, 21]. In general, it would be interesting to study the robustness of our results also when both collective and individual dissipation channels are present, as this is probably the most common case in experiments.

Finally we discuss the experimental feasibility of the chosen parameters values, noticing that most of them have been reached in experimental platforms (see [1] section IV for a review on the subject). Specifically, the dimensionless laser power is fixed to ${ \mathcal P }=12$, which in terms of the more commonly quoted g, is $g/{\omega }_{m}\approx 0.079$, and when varied takes values at most of $g/{\omega }_{m}\approx 0.19$ in figure 6, or $g/{\omega }_{m}\approx 0.37$ when the dynamical instability is reached in figures 4(e), (f). The values ranges in figure 4 have been achieved for instance in the experimental setups of references [71–73], with the exception of the highest range values explored for ${ \mathcal P }$ and ${K}_{c}/{\omega }_{m}^{2}$, for which previous considerations hold. The resolved sideband regime has been reached in numerous optomechanical devices, from the GHz mechanical frequencies of optomechanical crystals [53], to MHz frequencies in microtoroidal and LC resonators where values such as ${\omega }_{m}/\kappa \approx 11$ [72] or ${\omega }_{m}/\kappa \approx 63$ [73] are reported. Moreover in optimized optomechanical crystals [29, 74] bath phononic occupancies between 10 and 100 are reported for GHz mechanical modes, together with mechanical quality factors of the order of ${Q}_{m}\sim {10}^{4}\mbox{--}{10}^{6}$. Therefore, based on an analysis of actual state-of-the-art systems [1], this study represents a first step in establishing dynamical and quantum effects of collective dissipation expected to play a key role towards miniaturization of optomechanical devices.

This work has been supported by the EU through the H2020 Project QuProCS (Grant Agreement 641277), by MINECO/AEI/FEDER through projects NoMaQ FIS2014-60343-P, QuStruct FIS2015-66860-P and EPheQuCS FIS2016-78010-P.

In this appendix we write down the main steps to derive the Langevin equations for the mechanical oscillators in the CB case. Notice that in this section we are not using the dimensionless variables and parameters introduced in section 2.2. As usual we assume the cavities and the oscillators to dissipate independently so that we can derive the Langevin equations separately, as if there was no optomechanical coupling [1, 75]. First we write down the mechanical Hamiltonian in the normal mode basis together with the bath Hamiltonians written in equation (4):

$$\begin{eqnarray}&&\hat{{ \mathcal H }}=\displaystyle \frac{1}{2}{m}({{\rm{\Omega }}}_{+}^{2}{\hat{{x}}}_{+}^{2}+{{\rm{\Omega }}}_{-}^{2}{\hat{{x}}}_{-}^{2})+\displaystyle \frac{1}{2{m}}({\hat{{p}}}_{+}^{2}+{\hat{{p}}}_{-}^{2})+2{\lambda }\displaystyle \sum _{{\alpha }=1}^{{\rm{\infty }}}{\hat{{x}}}_{+}{\hat{{q}}}_{{\alpha }}+{\hat{{ \mathcal H }}}_{{\rm{bath}}}^{{\rm{CB}}},\end{eqnarray} \tag{ A.1 }$$

where we recall that we are assuming identical oscillators, and thus ${{\rm{\Omega }}}_{+}={\omega }_{m}$, ${{\rm{\Omega }}}_{-}=\sqrt{{\omega }_{m}^{2}+2k/m}$. Making the rotating wave approximation in the system-bath coupling and writing the Hamiltonian using the creation and annihilation operators we arrive to:

$$\begin{eqnarray}&&\hat{{ \mathcal H }}={\hslash }{{\rm{\Omega }}}_{+}{\hat{{b}}}_{+}^{\dagger }{\hat{{b}}}_{+}+{\hslash }{{\rm{\Omega }}}_{-}{\hat{{b}}}_{-}^{\dagger }{\hat{{b}}}_{-}+2{\gamma }\displaystyle \sum _{{\alpha }=1}^{{\rm{\infty }}}{\hslash }({\hat{{b}}}_{+}^{\dagger }{\hat{{r}}}_{{\alpha }}+{\hat{{b}}}_{+}{\hat{{r}}}_{{\alpha }}^{\dagger })+{\hat{{ \mathcal H }}}_{{\rm{bath}}}^{{\rm{CB}}},\end{eqnarray} \tag{ A.2 }$$

where $2\gamma$ is the CB dissipation rate which is twice the SB one γ. Now we can obtain the Langevin equations of motion for the modes ${x}_{\pm }$ following the standard input–ouput formalism [76]:

$$\begin{eqnarray}&&{\dot{\hat{b}}}_{+}=-{\rm{i}}{{\rm{\Omega }}}_{+}{\hat{b}}_{+}-2\gamma {\hat{b}}_{+}+\sqrt{4\gamma }{\hat{b}}_{{\rm{in}}},\quad {\dot{\hat{b}}}_{-}=-{\rm{i}}{{\rm{\Omega }}}_{-}{\hat{b}}_{-},\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&{\dot{\hat{b}}}_{+}^{\dagger }={\rm{i}}{{\rm{\Omega }}}_{+}{\hat{b}}_{+}^{\dagger }-2\gamma {\hat{b}}_{+}^{\dagger }+\sqrt{4\gamma }{\hat{b}}_{{\rm{in}}}^{\dagger },\quad {\dot{\hat{b}}}_{-}^{\dagger }={\rm{i}}{{\rm{\Omega }}}_{-}{\hat{b}}_{-}^{\dagger },\end{eqnarray} \tag{ A.4 }$$

where the equations for the mode ${x}_{-}$ are just the Heisenberg equations of motion, since this mode is not coupled to a thermal bath. In the Markovian limit the zero mean Gaussian noise terms are characterized by the following correlations [34]:

$$\begin{eqnarray}&&\langle {\hat{b}}_{{\rm{in}}}^{\dagger }(t){\hat{b}}_{{\rm{in}}}(t^{\prime} )\rangle ={n}_{{\rm{th}}}\delta (t-t^{\prime} ),\quad \langle {\hat{b}}_{{\rm{in}}}(t){\hat{b}}_{{\rm{in}}}^{\dagger }(t^{\prime} )\rangle =({n}_{{\rm{th}}}+1)\delta (t-t^{\prime} ),\end{eqnarray} \tag{ A.5 }$$

where n_th is the thermal occupancy of the phonon bath defined previously. From these equations we can obtain the equations for the position and momentum operators:

$$\begin{eqnarray}&&{\dot{\hat{x}}}_{+}=\displaystyle \frac{{\hat{p}}_{+}}{m}-2\gamma {\hat{x}}_{+}+\sqrt{4\gamma }\,{\tilde{x}}_{{zpf}}\left[\displaystyle \frac{{\hat{b}}_{{\rm{in}}}+{\hat{b}}_{{\rm{in}}}^{\dagger }}{\sqrt{2}}\right],\quad {\dot{\hat{x}}}_{-}=\displaystyle \frac{{\hat{p}}_{-}}{m},\end{eqnarray} \tag{ A.6 }$$

$$\begin{eqnarray}&&{\dot{\hat{p}}}_{+}=-m{{\rm{\Omega }}}_{+}^{2}{\hat{x}}_{+}-2\gamma {\hat{p}}_{+}+\sqrt{4\gamma }\,m{{\rm{\Omega }}}_{+}{\tilde{x}}_{{zpf}}\left[{\rm{i}}\displaystyle \frac{{\hat{b}}_{{\rm{in}}}^{\dagger }-{\hat{b}}_{{\rm{in}}}}{\sqrt{2}}\right],\quad \quad {\dot{\hat{p}}}_{-}=-m{{\rm{\Omega }}}_{-}^{2}{\hat{x}}_{-}.\end{eqnarray} \tag{ A.7 }$$

Finally we move to the coupled oscillators picture, where ${\hat{x}}_{\pm }=({\hat{x}}_{1}\pm {\hat{x}}_{2})/\sqrt{2}$ and ${\hat{p}}_{\pm }=({\hat{p}}_{1}\pm {\hat{p}}_{2})/\sqrt{2}$, obtaining:

$$\begin{eqnarray}&&{\dot{\hat{x}}}_{j}=\displaystyle \frac{{\hat{p}}_{j}}{m}-\gamma ({\hat{x}}_{1}+{\hat{x}}_{2})+{\tilde{x}}_{{zpf}}{\hat{x}}_{{\rm{in}}},\quad j=1,2\end{eqnarray} \tag{ A.8 }$$

$$\begin{eqnarray}&&{\dot{\hat{p}}}_{j}=-m{\omega }_{m}^{2}{\hat{x}}_{j}-\gamma ({\hat{p}}_{1}+{\hat{p}}_{2})+{(-1)}^{j}k({\hat{x}}_{1}-{\hat{x}}_{2})+m{\omega }_{m}{\tilde{x}}_{{zpf}}{\hat{p}}_{{\rm{in}}},\end{eqnarray} \tag{ A.9 }$$

where we have defined ${\hat{x}}_{{\rm{in}}}=\sqrt{\gamma }({\hat{b}}_{{\rm{in}}}+{\hat{b}}_{{\rm{in}}}^{\dagger })$ and ${\hat{p}}_{{\rm{in}}}={\rm{i}}\sqrt{\gamma }({\hat{b}}_{{\rm{in}}}^{\dagger }-{\hat{b}}_{{\rm{in}}})$. From these equations we can obtain equations (11) and (12) by coupling the oscillators to the cavities and applying the nondimensionalization/rescaling procedure described in section 2.2. Furthermore, note that as the noise terms are the same for both oscillators there appear the cross-correlations defined in equation (16).