An all-optical feedback assisted steady state of an optomechanical array

The drive toward exploring quantum features at the mesoscopic and macroscopic scale has intensified in recent times with marked advances in understanding of the quantum–classical transition [1]. The prime motivations have been to test the limits of the quantum theory and to search for various quantum communication protocols [2] while discovering novel applications of mesoscopic systems as components of quantum memories and repeaters required in the implementation of various quantum information processing tasks [3, 4] and new sensor technology [5, 6].

In this respect vastly different physical systems such as ultra-cold atoms coupled to nanomechanical systems [7–10], superconducting qubits coupled to microwave resonators [11] and opto- and electro-mechanical cavities [12–34] have recently attracted a lot of attention in probing quantum signatures at the mesoscopic scale. By and large the goal has been to achieve coherent control over the motion of the mechanical element via coupling to electromagnetic and optical modes in the quantum regime [35]. Such a coupling inevitably has the potential to engineer a hybrid quantum device composed of otherwise incompatible degrees of freedom of different physical systems. Among the plethora of such proposals, special attention has been devoted to explore further the potentials of optomechanics, in particular. Optomechanics makes use of radiation pressure of light on a mechanical element resulting in reversible coupling between optics and mechanics, two vastly different physical systems at the extremes of the quantum world. Recent progress in the field not only promises to strongly couple these very disparate modes identifying a whole new class of engineered quantum systems but also to control the dynamics of such systems [36].

In an optomechanical (OM) set-up, the resonance frequency of the optical cavity under consideration can be altered by the displacement of the mechanical element which may for example form one of the mirrors in a Fabry–Pérot cavity. As a result of radiation pressure of light, reversible coupling is generated between the optical and mechanical modes. The strength of the OM coupling is however derived from the magnitude of the radiation pressure force which can be enhanced by strong coherent driving in the considered OM system. As has been extensively shown in the resolved sideband regime [37–43] (see Hartmann and Plenio [24]), OM interaction can generate entanglement between the optical and mechanical modes within an OM system. Further for two or more coupled OM systems forming an OM array, intra cavity OM entanglement can be distributed over inter cavity modes in the steady state [34, 44, 45]. In this work we are interested in devising a scheme whereby entanglement can not only be generated but can also be controlled as well. The ability to control entanglement would be a useful resource in the design and implementation of quantum repeaters and quantum memories crucial in the implementation of quantum information processing. It may also enable new measurement protocols based on mechanical sensors, for example, weak force sensors [6, 46]. With this objective in mind we apply an all-optical feedback scheme to earlier proposals, [45] considering coupled OM arrays. All-optical feedback was considered extensively in the context of cascaded cavities [47] and associated amplitude and quadrature feedback.

We apply similar coupling configurations to a system of OM cavities in our analysis and compare the steady state dynamics of various OM networks both in the presence and absence of an all-optical feedback loop. Specifically we address two different topologies: firstly we consider a pair of independent but reversibly coupled OM systems. Secondly we consider an OM arrangement where two different optical modes are coupled via a common mechanical element. In fact our results also apply to a hybrid system in which one optical mode and one microwave mode interact via a common mechanical element as in [48, 49]. In each case the all-optical feedback loop is implemented by an additional irreversible coupling between the two optical modes. Our results show that while a reversibly coupled OM network can result in steady state entanglement between the inter cavity mechanical modes, implementing an all optical feedback loop can enhance the generated entanglement. Hence an all-optical feedback can be used as a resource to control entanglement in the steady state of an OM network. Furthermore in the last section of this work, we extend our analysis to a large array comprising of reversibly coupled OM ports. Each OM port in the array is composed of a pair of OM systems coupled via an all-optical feedback loop. Our results show bipartite steady state entanglement only exists between mechanical resonators from different OM ports with different parity, and can be enhanced in the presence of feedback.

2.1. Two coupled OM cavities

We begin by first considering the set-up in figure 1 which comprises of two OM cavities (OM 1 and 2), coupled reversibly via their optical modes. Additionally the optical mode a₁ of OM 1 is also coupled irreversibly to mode a₂ of OM 2. Simultaneous presence of unidirectional and bidirectional coupling between the optical modes constitute an all-optical feedback loop [47]. The system OM 1 acts as a source driving the state of the OM 2 cavity which in turn influences the state of source cavity through reversible coupling of arbitrary strength κ.

Zoom In Zoom Out Reset image size

In the usual convention, radiation pressure is induced by strong external coherent driving on each OM system hence displacing the cavity field by its steady state amplitude. The resulting OM interaction is then linearized such that a quadratic coupling is induced between the optical and mechanical modes in each system [45]. In this linearized regime and in the rotating frame of the external drive, the OM interaction between each individual cavity takes the form (in units ℏ = 1) [45]

$$\begin{equation} {H}_{\mathrm{I}}= \sum_{j=1}^{2} \Delta_{j}\hat{a}_{j}^{\dagger}\hat{a}_{j}+\omega_{m_j}\hat{b}_{j}^{\dagger}\hat{b}_{j}+ g_{j}(\hat{a}_{j}+ \hat{a}_{j}^{\dagger})(\hat{b}_{j}+ \hat{b}_{j}^{\dagger}), \end{equation} \tag{ 1 }$$

where Δ_j = Ω_j − ω_Lj, where ω_Lj is the frequency of the external laser drive pumping the jth cavity with resonance frequency Ω_j, ω_mj is the fundamental resonance frequency of the jth mechanical mode, g_j is the effective OM coupling strength proportional to the steady state amplitude of the cavity field due to linearization of the radiation pressure force and $\hat {a}_{j}(\hat {b}_{j})$ and $\hat {a}_{j}^{\dagger }(\hat {b}_{j}^{\dagger })$ are the respective annihilation and creation operators for the optical (mechanical) modes for each cavity. In the presence of the OM interaction of the form (1), the optical and mechanical modes of each individual OM network will become quantum entangled. In this work we are interested in achieving long-distance entanglement in coupled OM cavity arrays. Moreover, we want to explore the influence of an all-optical feedback on inter-cavity entanglement. In the presence of an all-optical feedback and under the Born–Markov approximation [50], the open system dynamics of the coupled OM cavities is given by the following master equation:

$$\begin{eqnarray} \frac{{\rm d}\rho}{{\rm d}t}&=& -\mathrm{i} [{H}_{\mathrm{I}},\rho]-\mathrm{i}\kappa [\hat{a}_1 \hat{a}_2^\dagger \,{\rm e}^{{\rm i}(\omega_{L_2} -\omega_{L_1})t}+\hat{a}_2 \hat{a}_1^\dagger\,{\rm e}^{{\rm i}(\omega_{L_1} -\omega_{L_2})t}, \rho] \nonumber \\ & &+\sqrt{\Gamma_{1}\Gamma_{2}}\left([\hat{a}_{1}\rho,\hat{a}_{2}^{\dagger}]\,{\rm e}^{{\rm i}(\omega_{L_2} -\omega_{L_1})t}+[\hat{a}_{2},\rho \hat{a}_{1}^{\dagger}]\,{\rm e}^{{\rm i}(\omega_{L_1} -\omega_{L_2})t}\right) \nonumber \\ & & +\sum_{j=1}^2\Gamma_j{\cal D}[\hat{a}_j]\rho +\gamma_{j} (\bar{n}_{j}+1){\cal D}[\hat{b}_{j}]\rho+\gamma_{j} \bar{n}_{j}{\cal D}[\hat{b}_{j}^\dagger]\rho, \end{eqnarray} \tag{ 2 }$$

where the linewidth of each OM cavity is given by Γ_j while the mechanical resonators are allowed to decay to a thermal bath with average thermal occupancy $\bar {n}_{j}$ at a rate γ_j. The damping superoperator in equation (2) ${\cal D}[A]$ is defined in the usual way by

$$\begin{equation} {\cal D}[A]\rho=A\rho A^\dagger -\textstyle\frac{1}{2} (A^\dagger A\rho+\rho A^\dagger A). \end{equation} \tag{ 3 }$$

In the interaction picture of the unitary operator $\hat {U}=\mathrm {e}^{\sum _{j=1}^{2}\mathrm {i}(\Delta _{j}\hat {a}^{\dagger }_{j}\hat {a}_{j}+\omega _{m_j}\hat {b}^{\dagger }_{j}\hat {b}_{j})t}$, the master equation (2) takes the form

$$\begin{eqnarray}&&\fl \frac{{\rm d}\rho}{{\rm d}t}= -\mathrm{i} \sum_{j=1}^{2}g_{j}[(\hat{a}_{j}\,{\rm e}^{-{\rm i} \Delta_{j} t}+\hat{a}^{\dagger}_{j}\,{\rm e}^{{\rm i} \Delta_{j} t})(\hat{b}_{j}\,{\rm e}^{-{\rm i} \omega_{m_j} t}+\hat{b}^{\dagger}_{j}\,{\rm e}^{{\rm i} \omega_{m_j}t}),\rho]-\mathrm{i}\kappa [\hat{a}_1 \hat{a}_2^\dagger \,{\rm e}^{{\rm i}(\omega_{L_2} -\omega_{L_1} + \Delta_{2}- \Delta_{1})t}\nonumber \\ &&+\ \hat{a}_2 \hat{a}_1^\dagger \,{\rm e}^{{\rm i}(\omega_{L_1} -\omega_{L_2}+ \Delta_{1}- \Delta_{2})t}, \rho] +\sqrt{\Gamma_{1}\Gamma_{2}}([\hat{a}_{1}\rho,\hat{a}_{2}^{\dagger}]\,{\rm e}^{{\rm i}(\omega_{L_2} -\omega_{L_1}+ \Delta_{2}- \Delta_{1})t}+ [\hat{a}_{2},\rho \hat{a}_{1}^{\dagger}]\nonumber \\ &&\times\,{\rm e}^{{\rm i}(\omega_{L_1} -\omega_{L_2}+ \Delta_{1}- \Delta_{2})t}) +\sum_{j=1}^2\Gamma_j{\cal D}[\hat{a}_j]\rho +\gamma (\bar{n}+1){\cal D}[\hat{b}_{j}]\rho+\gamma \bar{n}{\cal D}[\hat{b}_{j}^\dagger]\rho. \end{eqnarray} \tag{ 4 }$$

For simplicity, we take both OM systems to be identical such that Ω₁ = Ω₂, both mechanical resonators have the same frequency, ω_mj = ω_m, mechanical damping rate γ_j = γ and thermal occupation number $\bar {n}_{j}=\bar {n}$. Each OM cavity is however driven on a different sideband. In particular we drive OM 1 exclusively on the blue sideband, Δ₁ = −ω_m and OM 2 on the red sideband, Δ₂ = ω_m. Under the rotating-wave approximation (RWA) [45, 50], the OM interaction simplifies to the following form:

$$\begin{equation} {\skew3\tilde{H}}_{\mathrm{I}}= g_{1}(\hat{a}_{1}\hat{b}_{1} + \hat{a}_{1}^{\dagger}\hat{b}_{1}^{\dagger}) + g_{2}(\hat{a}_{2}^{\dagger}\hat{b}_{2}+\hat{a}_{2}\hat{b}_{2}^{\dagger}) \end{equation} \tag{ 5 }$$

and the master equation describing the open system dynamics of the coupled OM network takes the following form:

$$\begin{eqnarray} \frac{{\rm d}\rho}{{\rm d}t}&=& -\mathrm{i} [{\skew3\tilde{H}}_{\mathrm{I}},\rho]-\mathrm{i}\kappa [\hat{a}_1 \hat{a}_2^\dagger +\hat{a}_2 \hat{a}_1^\dagger, \rho] +\sqrt{\Gamma_{1}\Gamma_{2}}([\hat{a}_{1}\rho,\hat{a}_{2}^{\dagger}]+[\hat{a}_{2},\rho \hat{a}_{1}^{\dagger}]) \nonumber \\ & & +\sum_{j=1}^2\Gamma_j{\cal D}[\hat{a}_j]\rho +\gamma (\bar{n}+1){\cal D}[\hat{b}_{j}]\rho+\gamma \bar{n}{\cal D}[\hat{b}_{j}^\dagger]\rho. \end{eqnarray} \tag{ 6 }$$

The first two commutators on the right hand side of the above equation denotes the unitary evolution of the coupled OM network, the successive two commutators represents the irreversible coupling between the optical modes of the first and the second cavity and the last six terms takes into account the damping of optical and mechanical modes. Reversible coupling between the optical modes can be generated either through evanescent coupling of the adjacent optical modes or through an optical fiber [34]. Irreversible coupling between the source and the driven cavity can be established using an optical circulator [45, 51]. A non-reciprocal optical device such as Faraday rotator can be used to establish irreversible coupling between the optical modes of two adjacent OM cavities.

This particular choice of detuning configurations in such a coupled array induces a two mode squeezing interaction between the optical and mechanical modes in OM 1, and a beam splitter interaction resulting in state transfer between optical and mechanical modes in OM 2. This configuration of coupled dynamics has been analyzed in detail in [45], where two independent coupling configurations between individual OM systems were analyzed, i.e. the OM systems were considered to be either reversibly or irreversibly coupled via their optical ports. In the present work, however we extend the consideration to an added coupling channel between coupled OM systems 1 and 2. This additional coupling is engineered as an unidirectional coupling emanating from the optical mode of OM cavity 1 to the optical mode of OM cavity 2, while also maintaining a reversible coupling between the two inter cavity optical modes: hence implementing an all-optical feedback loop [47] in the OM array. As has been established previously, coupling between such driven OM systems results in distribution of entanglement generated between optical and mechanical modes in OM cavity 1 over inter cavity modes in the steady state [45]. However in the present work we are interested in exploring the possibility of controlling/altering the steady state inter cavity phonon–phonon entanglement via an all-optical feedback loop.

Alternatively such a network of coupled cascaded quantum systems shown in figure 1 can also be conveniently modeled using the SLH framework [52, 53]. H is the Hamiltonian operator for the system, L is the coupling vector and S is the scattering matrix. The SLH formalism is useful in forming reducible networks of coupled open quantum systems with the objective of implementing quantum control analysis and design. In the SLH formalism, an open quantum system [3] is described by a set of three parameters, G = (S,L,H) where S denotes the scattering matrix, L is the coupling vector and H is the Hamiltonian operator for the system. Hence the set of coupled quantum OM systems summarized in figure 1 forms a series product connection under the SLH formalism [52]

$$\begin{eqnarray} {\bf G_{2}} \triangleleft {\bf G_{1}}&=&({\bf S}_{\rm eff},{\bf L}_{\rm eff},{\bf H_{\rm eff}})\nonumber \\ &=&\left({\bf S_{2}}{\bf S_{1}},{\bf L_{2}}+{\bf S_{2}}{\bf L_{1}},H_{1}+H_{2}+\frac{1}{2\mathrm{i}}({\bf L_{2}}^{\dagger}{\bf S_{2}}{\bf L_{1}}-{\bf L_{1}}^{\dagger}{\bf S_{2}}^{\dagger}{\bf L_{2}})\right), \end{eqnarray} \tag{ 7 }$$

where the specific variables of the series product are S_j =  I, ${\bf{ L_{j}}}=\sqrt {\Gamma _{j}}\hat {a}_{j}$ corresponding to each OM system and the combined coupling vector is ${\bf{ L}}_{\mathrm { eff}}=\sqrt {\Gamma _{1}}\hat {a}_{1} +\sqrt {\Gamma _{2}}\hat {a}_{2}$. The Hamiltonian operator of the combined system is given as

$$\begin{equation} {\bf H_{eff}}={\skew3\tilde{H}}_{\mathrm{I}}+\kappa(\hat{a}_1 \hat{a}_2^\dagger +\hat{a}_2 \hat{a}_1^\dagger)+\frac{1}{2\mathrm{i}}\left({\bf L_{2}}^{\dagger}{\bf S_{2}}{\bf L_{1}}- {\bf L_{1}^{\dagger}}{\bf S_{2}}{\bf L_{2}}\right), \end{equation} \tag{ 8 }$$

where ${\tilde {H}}_{\mathrm {I}}$ has been defined earlier in equation (5). The master equation of the reduced network can then be deduced from the three variables as [52]

$$\begin{equation} \frac{\rm d}{{\rm d}t}\rho=\mathrm{i}[\rho,{\bf H_{eff}}]+{\cal D}[{\bf L_{eff}}]\rho. \end{equation} \tag{ 9 }$$

Substituting the expressions for H_eff and ${\cal D}[{\bf{ L_{eff}}}]$ in the above equation one gets

$$\begin{eqnarray*} \fl \frac{\rm d}{{\rm d}t}\rho &=&-\mathrm{i} [{\skew3\tilde{H}}_{\mathrm{I}},\rho]-\mathrm{i}\kappa [\hat{a}_1 \hat{a}_2^\dagger +\hat{a}_2 \hat{a}_1^\dagger, \rho]+\Gamma_{1}\mathcal{D}[\hat{a}_{1}]\rho +\Gamma_{2} \mathcal{D}[\hat{a}_{2}]\rho -\frac{\sqrt{\Gamma_{1}\Gamma_{2}}}{2}[\hat{a}_1 \hat{a}_2^\dagger -\hat{a}_2 \hat{a}_1^\dagger,\rho ]\nonumber \\\fl &&\quad + \sqrt{\Gamma_{1}\Gamma_{2}}\left(\hat{a}_{1}\rho \hat{a}_{2}^{\dagger}-\frac{1}{2}\hat{a}_{1}^{\dagger}\hat{a}_{2}\rho-\frac{1}{2}\rho\hat{a}_{1}^{\dagger}\hat{a}_{2}\right)+\sqrt{\Gamma_{1}\Gamma_{2}}\left(\hat{a}_{2}\rho \hat{a}_{1}^{\dagger}-\frac{1}{2}\hat{a}_{2}^{\dagger}\hat{a}_{1}\rho-\frac{1}{2}\rho\hat{a}_{2}^{\dagger}\hat{a}_{1}\right). \end{eqnarray*}$$

Simplifying the above equation we get

$$\begin{eqnarray} \frac{\rm d}{{\rm d}t}\rho &=&-\mathrm{i} [{\skew3\tilde{H}}_{\mathrm{I}},\rho]-\mathrm{i}\kappa [\hat{a}_1 \hat{a}_2^\dagger +\hat{a}_2 \hat{a}_1^\dagger, \rho]+\Gamma_{1}\mathcal{D}[\hat{a}_{1}]\rho +\Gamma_{2} \mathcal{D}[\hat{a}_{2}]\rho \nonumber \\ & &+ \sqrt{\Gamma_{1}\Gamma_{2}}(\hat{a}_{1}\rho \hat{a}_{2}^{\dagger}+\hat{a}_{2}\rho \hat{a}_{1}^{\dagger}-\hat{a}_{1} \hat{a}_{2}^{\dagger} \rho-\rho\hat{a}_{2} \hat{a}_{1}^{\dagger} ). \end{eqnarray} \tag{ 10 }$$

Taking into account the damping of mechanical modes b₁ and b₂ by coupling them to identical independent thermal reservoirs, the dynamics of the coupled OM network is then given by

$$\begin{eqnarray}&&\fl \frac{\rm d}{{\rm d}t}\rho =-\mathrm{i} [{\skew3\tilde{H}}_{\mathrm{I}},\rho]-\mathrm{i}\kappa [\hat{a}_1 \hat{a}_2^\dagger +\hat{a}_2 \hat{a}_1^\dagger, \rho] + \sqrt{\Gamma_{1}\Gamma_{2}}(\hat{a}_{1}\rho \hat{a}_{2}^{\dagger}+\hat{a}_{2}\rho \hat{a}_{1}^{\dagger}-\hat{a}_{1} \hat{a}_{2}^{\dagger} \rho-\rho\hat{a}_{2} \hat{a}_{1}^{\dagger} )\nonumber \\ &&+\sum_{j=1}^2\Gamma_j{\cal D}[\hat{a}_j]\rho +\gamma (\bar{n}+1){\cal D}[\hat{b}_{j}]\rho+\gamma \bar{n}{\cal D}[\hat{b}_{j}^\dagger]\rho. \end{eqnarray} \tag{ 11 }$$

It is thus clear that the master equation (11) is identical to the master equation (6) of the OM network.

It is easy to see from equation (6) that the Hamiltonian, the Lindblad operators and the all-optical feedback loop are all bilinear in bosonic operators. Thus it is guaranteed that if the coupled OM array is initially prepared in a Gaussian state, then the evolution described by the master equation (6) maintains this Gaussian character. Thus to solve the master equation (6) for an initial Gaussian state, it is consistent to make a Gaussian ansatz for the normal-ordered quantum characteristic function of the following form $\chi (\epsilon ,\eta ,x,y,{t}) =\exp [-z^{\mathrm {T}}{{\bf{ A}}}(t)\,z+\mathrm {i} z^{\mathrm {T}}h(t)]$, where A(t) is a time-dependent 8 × 8 matrix and h(t) is a 1 × 8 time-dependent vector and z^T = (,*,η,η*,x,x*,y,y*). An explicit solution of the master equation (6) is outlined in appendix A. From the solution of the normal-ordered quantum characteristic function it is a straightforward exercise to calculate all the correlators between the desired modes.

To quantize the entanglement between any two modes, we utilize the logarithmic negativity as a measure of entanglement for Gaussian states [54]. For a two-mode Gaussian continuous variable state characterized by its covariance matrix V, the logarithmic negativity is computed as $\mathcal {N}= {\mathrm { Max}}[0,-{\mathrm { log}}(2 \nu _{-})]$ [54] where ν₋ is the smallest of the symplectic eigenvalues of the covariance matrix given by $\nu _{-}= \sqrt {\sigma /2-\sqrt {(\sigma ^{2}-4\,{\mathrm { Det}} {{\bf{ V}}})}/{\mathrm { 2}}}$. Here σ = Det  A₁ + Det  B₁ − 2 Det  C₁ and ${\bf{ V}}= \left(\!{{{\bf{ A}_{1}}} \atop {{\bf{ C}_{1}^{T}}} } \; {{{\bf{ C}_{1}}} \atop {{\bf{ B}_{1}}}}\!\right)$, where A₁ ( B₁) accounts for the local variances of the two modes and C₁ for the inter-mode correlations.

We compute logarithmic negativity measuring the quantum entanglement between the mechanical modes b₁ and b₂ in the steady state and in figures 2(a) and (b) plot it as a function of the reversible coupling strength between the optical modes κ and the thermal occupancy of each mechanical mode $\bar {n}$. We find finite steady state entanglement between the two inter cavity mechanical modes in the steady state. We observe the entanglement persists for a finite amount of non-zero thermal occupation number $\bar {n}$. We also note the entanglement between the mechanical modes is higher for smaller values of coupling strength κ between the optical modes.

Zoom In Zoom Out Reset image size

To explore the advantage of using an all-optical feedback, in figures 2(a) and (b) we also plot the entanglement between the mechanical modes b₁ and b₂ in absence of all-optical feedback. Clearly the presence of feedback increases the degree of entanglement between inter cavity phonons. It can be seen from figure 2(a) that presence of reversible and irreversible couplings between the optical modes does have an impact on the entanglement between the mechanical modes. The influence is clearly marked for small values of κ and under which case feedback generates stronger entanglement between the mechanical modes. From figure 2(b) it is clear that an all-optical feedback results in generating entangled state of mechanical modes which is more robust to thermal fluctuations.

It is also interesting to see how the inter cavity phonon–phonon entanglement $\mathcal {N}$, varies as a function of the different OM coupling strengths g₁ and g₂ in the OM array. These are plotted in figures 3(a) and (b) in the presence and absence of feedback, respectively. Hence the added coupling channel between the inter cavity optical modes forming the feedback loop, can be used to increase the steady state entanglement in such an OM array.

Zoom In Zoom Out Reset image size

To gain some more insight into the constructive role played by an all-optical feedback in coupled OM array we consider a pragmatic regime where the linewidth of each cavity mode is very large compared to all other system parameters. Under this regime both the cavity modes can be adiabatically eliminated and stochastic differential equations for the mechanical modes can be arrived at. The detailed calculation is outlined in appendix B. The effect of feedback is clearly imprinted on the steady state correlators between the mechanical modes given by equations (B.5) and (B.6). The steady state correlator $|\langle \hat {b}_{1} \hat {b}_{2} \rangle |$ in presence and absence of feedback is plotted in figure 4. As can be clearly seen from figure 4, feedback has a clear influence on quantum correlations between the two mechanical modes. In particular for low values of κ, feedback helps in building stronger quantum correlations between the two mechanical modes. For values of κ ∼ Γ, there is no evident difference in the strength of quantum correlations between the mechanical modes in the presence and absence of all-optical feedback in OM network array. This observation is also imprinted in figure 2(a).

Zoom In Zoom Out Reset image size

2.2. Mechanics assisted optical coupling

In the present section we shall discuss the effect of all-optical feedback in a different setting of coupled optical and mechanical modes. We consider a scenario where the two optical modes are coupled to a common mechanical mode with a bilinear interaction as shown in figure 5. The Hamiltonian describing the interaction between the three coupled bosonic modes thus takes the following form:

$$\begin{equation} \skew3\tilde{H}=\sum_{j=1,2}\Delta_{j}\hat{a}_{j}^{\dagger} \hat{a}_{j} +\omega_{m} \hat{b}^{\dagger} \hat{b}+(g_{1}(\hat{a}_{1}^{\dagger}+ \hat{a}_{1} )+g_{2}(\hat{a}_{2}+\hat{a}_{2}^{\dagger} ))(\hat{b}+\hat{b}^{\dagger}), \end{equation} \tag{ 12 }$$

where Δ_j = Ω_j − ω_Lj, where ω_Lj is the frequency of the external laser drive pumping the jth cavity whose resonance frequency is Ω_j. The above Hamiltonian, for instance, can also describe the coherent coupling between three bilinearly coupled harmonic oscillators arranged in an open chain. In pursuit of main aim of this work we consider a scenario where in addition to an indirect interaction between the two optical modes, the two optical cavity modes also interact with a forward feed coupling. Under the Born–Markov approximation, the open system dynamics of the coupled OM cavities then takes the following form:

$$\begin{eqnarray} &&\fl \frac{{\rm d} \rho}{{\rm d} t}=-\mathrm{i}[\skew3\tilde{H},\rho]+\sqrt{\Gamma_{1} \Gamma_{2}} [\hat{a}_{1} \rho, \hat{a}_{2}^{\dagger} ]\,{\rm e} ^{{\rm i} (\omega_{L_{2}}-\omega_{L_{1}})t} +\sqrt{\Gamma_{1} \Gamma_{2}} [\hat{a}_{2} , \rho \hat{a}_{1}^{\dagger}]\,{\rm e} ^{{\rm i} (\omega_{L_{1}}-\omega_{L_{2}})t}\nonumber\\ &&+\sum_{j=1,2} \frac{\Gamma_{j}}{2}\mathcal{L}_{\hat{a}_{j}} \rho + \frac{\gamma_{1}}{2}(\bar{n}+1)\mathcal{L}_{\hat{b}} \rho+\frac{\gamma_{1}}{2}(\bar{n})\mathcal{L}_{\hat{b}^{\dagger}} \rho. \end{eqnarray} \tag{ 13 }$$

Working in the interaction picture of the bare frequency of the optical and the mechanical modes and further choosing the physical parameters such that Δ₁ = −ω_m and Δ₂ = ω_m, then under the RWA the coupled system dynamics is governed by the following master equation:

$$\begin{eqnarray}&&\fl \frac{{\rm d} \rho}{{\rm d} t}=-\mathrm{i} g_{1}[\hat{a}_{1}^{\dagger}\hat{b}^{\dagger}+\hat{a}_{1}\hat{b},\rho]-\mathrm{i} g_{2}[\hat{a}_{2}^{\dagger}\hat{b}+\hat{b}^{\dagger} \hat{a}_{2},\rho] +\sqrt{\Gamma_{1} \Gamma_{2}} [\hat{a}_{1} \rho, \hat{a}_{2}^{\dagger} ] +\sqrt{\Gamma_{1} \Gamma_{2}} [\hat{a}_{2} , \rho \hat{a}_{1}^{\dagger}]\nonumber\\ &&+\sum_{j=1,2} \frac{\Gamma_{j}}{2}\mathcal{L}_{\hat{a}_{j}} \rho + \frac{\gamma_{1}}{2}(\bar{n}+1)\mathcal{L}_{\hat{b}} \rho+\frac{\gamma_{1}}{2}(\bar{n})\mathcal{L}_{\hat{b}^{\dagger}} \rho, \end{eqnarray} \tag{ 14 }$$

where an explicit time independent form of the master equation has been arrived at by choosing the physical parameters such that Ω₁ = Ω₂ = (ω_L1 + ω_L2)/2 and ω_m = (ω_L1 − ω_L2)/2. Taking a similar approach to the previous section, it is possible to solve the master equation (14) and the solution is outlined in appendix C.

Zoom In Zoom Out Reset image size

From the steady state solution of the master equation (14) we compute the logarithmic negativity between the optical and mechanical modes. We plot the logarithmic negativity as an indicator of steady state entanglement between the optical and mechanical modes as a function of OM coupling strength g₁ in figure 6. The entanglement between the modes a₁ and b is plotted in figure 6(a), while figure 6(b) shows the entanglement between the modes a₂ and b. Each figure compares how the steady state entanglement between the relevant mode varies in the presence as well as absence of all-optical feedback. We note that feedback does not bring about qualitative change in the steady state of optical and mechanical modes a₁,b. On the other hand we observe that the presence of feedback generates entanglement between otherwise unentangled modes a₂ and b. The degree of entanglement between the modes a₂ and a₁ (not shown here) remains relatively small both in the presence and absence of the all-optical feedback.

Zoom In Zoom Out Reset image size

Some more insight into understanding the role of all-optical feedback in generating quantum entanglement between modes a₂ and b can be obtained by writing Langevin's equations of motion for the modes a₂ and b. In the absence of an all-optical feedback and adiabatically eliminating the cavity mode a₁ results in

$$\begin{equation} \frac{\rm d}{{\rm d}t}\hat{a}_{2}=-\mathrm{i}g_{2}\hat{b}-\frac{\Gamma}{2}\hat{a}_{2}+\sqrt{\Gamma}\hat{a}_{in}(2,t), \end{equation} \tag{ 15 }$$

$$\begin{equation} \frac{\rm d}{{\rm d}t}\hat{b}=-\mathrm{i}g_{2}\hat{a}_{2}-\left(\frac{\gamma}{2}-2\frac{g_{1}^{2}}{\Gamma}\right)\hat{b}-\frac{2 \mathrm{i} g_{1}}{\sqrt{\Gamma}}\hat{a}_{in}^{\dagger}(1,t)+\sqrt{\gamma}\hat{b}_{in}(1,t). \end{equation} \tag{ 16 }$$

It is a straightforward exercise to arrive at steady state expressions for the modes a₂ and b. In the steady state one obtains

$$\begin{eqnarray}&& \left( \begin{array}{@{}c@{}} \hat{a}_{2}(\infty)\\ \hat{b}(\infty) \end{array} \right)= \lim_{t \to \infty} \int_{0}^{\infty}\left( \begin{array}{@{}cc@{}} x_{1}(t') & x_{2}(t') \\ x_{3}(t') & x_{4}(t') \end{array} \right) A(t-t')\,{\rm d}t', \nonumber\\&& A(t-t') =\left( \begin{array}{@{}c@{}} \sqrt{\Gamma} \hat{a}_{in}(2,t-t') \\ (-2\mathrm{i}g_{1}/\sqrt{\Gamma})\hat{a}_{in}^{\dagger}(1,t-t')+\sqrt{\gamma} \hat{b}_{in}(1,t-t') \end{array} \right),\\ &&\left( \begin{array}{@{}cc@{}} x_{1}(t') & x_{2}(t') \\ x_{3}(t') & x_{4}(t') \end{array} \right)=\mathrm{e}^{Bt'}, B=\left( \begin{array}{@{}cc@{}} -\Gamma/2 & -\mathrm{i}g_{2} \\ -\mathrm{i}g_{2} & -(\gamma/2-2g_{1}^{2}/\Gamma) \end{array} \right).\nonumber \end{eqnarray} \tag{ 17 }$$

It is easy to verify that presence of completely uncorrelated noise operators $\hat {a}_{in}(1,t),\hat {a}_{in}(2,t)$ results in $\langle \hat {a}_{2}(\infty )\hat {b}(\infty )\rangle=0$. Thus in absence of all-optical feedback between the optical modes results in an uncorrelated steady state of modes a₂ and b.

On the other hand, in the presence of all-optical feedback between the optical modes a₁ and a₂ adiabatically eliminating the cavity mode a₁ results in

$$\begin{eqnarray}&&\fl \frac{\rm d}{{\rm d}t}\left( \begin{array}{@{}c@{}} \hat{a}_{2} \\ \hat{a}_{2}^{\dagger} \\ \hat{b} \\ \hat{b}^{\dagger} \end{array} \right)=\left( \begin{array}{@{}cccc@{}} -\Gamma/2 & 0 & -\mathrm{i}g_{2} & (2\mathrm{i}g_{1}/\Gamma)\sqrt{\Gamma^{2}} \\ 0 & -\Gamma/2 & -(2\mathrm{i}g_{1}/\Gamma)\sqrt{\Gamma^{2}} & \mathrm{i}g_{2} \\ -\mathrm{i}g_{2} & 0 & -(\gamma/2 -2g_{1}^{2}/\Gamma) & 0 \\ 0 & \mathrm{i}g_{2} & 0 & -(\gamma/2 -2g_{1}^{2}/\Gamma) \end{array} \right)\left( \begin{array}{@{}c@{}} \hat{a}_{2} \\ \hat{a}_{2}^{\dagger} \\ \hat{b} \\ \hat{b}^{\dagger} \end{array} \right)\nonumber \\ &&+\left( \begin{array}{@{}c@{}} -\sqrt{\Gamma}\hat{a}_{in}(1,t) \\ -\sqrt{\Gamma}\hat{a}_{in}^{\dagger}(1,t) \\ \sqrt{\gamma}\hat{b}_{in}(1,t)-(2\mathrm{i}g_{1}/\sqrt{\Gamma})\hat{a}_{in}^{\dagger}(1,t) \\ \sqrt{\gamma}\hat{b}_{in}^{\dagger}(1,t)+(2\mathrm{i}g_{1}/\sqrt{\Gamma})\hat{a}_{in}(1,t) \end{array} \right). \end{eqnarray} \tag{ 18 }$$

It is easy to obtain expressions for steady state correlations between modes a₂ and b. Assuming that the modes a₂ and b are in contact with a zero-temperature reservoir one arrives at

$$\begin{eqnarray}\fl &&\langle \hat{a}_{2}(\infty)\hat{b}(\infty)\rangle=\int_{0}^{\infty}(\Gamma y_{1,1}y_{3,2}+2\mathrm{i}g_{1}y_{1,1}y_{3,3}+\gamma y_{1,3}y_{3,4}-2\mathrm{i}g_{1}y_{1,4}y_{3,2}+\frac{4g_{1}^{2}}{\Gamma}y_{1,4}y_{3,3})\,{\rm d}t',~~\nonumber\\ &&\left( \begin{array}{@{}cccc@{}} y_{1,1} & y_{1,2} & y_{1,3} & y_{1,4} \\ y_{2,1} & y_{2,2} & y_{2,3} & y_{2,4} \\ y_{3,1} & y_{3,2} & y_{3,3} & y_{3,4} \\ y_{4,1} & y_{4,2} & y_{4,3} & y_{4,4} \\ \end{array} \right)=\mathrm{e}^{C t'},\nonumber \\ && C =\left( \begin{array}{@{}cccc@{}} -\Gamma/2 & 0 & -\mathrm{i}g_{2} & (2\mathrm{i}g_{1}/\Gamma)\sqrt{\Gamma^{2}} \\ 0 & -\Gamma/2 & -(2\mathrm{i}g_{1}/\Gamma)\sqrt{\Gamma^{2}} & \mathrm{i}g_{2} \\ -\mathrm{i}g_{2} & 0 & -(\gamma/2 -2g_{1}^{2}/\Gamma) & 0 \\ 0 & \mathrm{i}g_{2} & 0 & -(\gamma/2 -2g_{1}^{2}/\Gamma) \end{array} \right). \end{eqnarray} \tag{ 19 }$$

Thus in the presence of all-optical feedback between the optical modes a₁ and a₂, we clearly have a non-vanishing steady state correlation function between the modes a₂ and b. This feature is corroborated in figure 6(b).

We now extend our analysis to study the influence of all-optical feedback on a large network of coupled OM cavities. We consider an array of N-OM ports. Each 'port' itself comprises of a unit formed by two OM cavities coupled through an all-optical feedback, as described in section 2.1. To invoke a genuine OM array with non-trivial quantum features in the steady state, we consider a scenario where there is a finite inter-mode coupling between different ports. In the absence of an inter-port coupling, quantum mechanical evolution of each individual port is described by the master equation (6). Denoting ρ_n for the state of nth port in a network with N-coupled OM cavities, the joint state of the OM array is thus simply given by ρ = ⊗^N_n=1ρ_n. The optical and mechanical modes associated with an ith 'port' are labeled a_i,1,a_i,2 and b_i,1,b_i,2, respectively. The inter-port coupling between ith and i + 1th ports then takes the form

$$\begin{equation} V_{\rm int}=\chi (a^{\dagger}_{i,2}a_{i+1,1}+ a^{\dagger}_{i+1,1} a_{i,2}), \end{equation} \tag{ 20 }$$

where χ is the coupling strength. A schematic of such a model is shown in figure 7. This sets the stage to understand the influence of all-optical feedback on a network of coupled OM cavities. We now numerically solve the master equation for the network of coupled OM cavities with inter-port coupling of the form (20), where each individual port is evolving under the master equation (6). From the solution of the master equation we compute the steady state entanglement between distant mechanical modes and the result is plotted in figure 8. Our simulations show that entanglement can only be generated between mechanical modes of different parity, i.e. b_i,1 and b_i+1,2. This result arises because entanglement can only be generated between any two oscillators when one OM system driven on the blue sideband is coupled to the next OM system driven on the red sideband, and vice versa. However driving any two oscillators on the same sideband, does not result in a distribution or generation of entanglement [45]. Therefore mechanical modes such as b_i,1 and b_i+1,1 remain separable for all parameters regardless of the absence or presence of the feedback loop.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As illustrated in figure 8, we find once again that an all-optical feedback loop helps in generating stronger quantum entanglement between mechanical oscillators belonging to different ports. This observation follows from the results displayed in figure 2 which show that the presence of feedback enhances the entanglement between the inter cavity mechanical modes in each port. Hence in the OM array composed of N-ports, where each port has a similar topology to that discussed in section 2.1, we notice stronger entanglement between mechanical modes b_i,1 and b_i+1,2. Another noteworthy feature of figure 8 is that for smaller values of reversible coupling κ between the optical modes, more entanglement is generated between distant mechanical modes. However, the spatial range of entanglement is short. On the other hand, for larger values of inter-optical coupling κ (∝χ) spatial range of entanglement between the mechanical modes grows. This, however, happens at the cost of reduced magnitude of entanglement between the mechanical modes.

Optomechanics is rapidly emerging as an exciting branch of quantum technology where it is possible to engineer unique hybrid quantum devices composed of seemingly incompatible physical systems. In the present work we have focused on a regime where each individual OM cavity is very strongly pumped and thus allowing us to linearize the inherently nonlinear interaction between the optical and the mechanical modes. It should, however, be noted that there has also been interest in exploring the nonlinear OM interaction by going beyond this linearization technique and thus allowing one to explore the non-Gaussian steady states of OM cavities [55, 56].

Expanding from previous work on distribution of entanglement in OM arrays, our scheme illustrates an all-optical feedback can be implemented to control the distributed entanglement in coupled OM arrays as well as generate otherwise non-existent quantum correlations as in the mechanics assisted OM interaction scheme. It should also be noted that our control protocols apply equally well to hybrid optical–microwave schemes as in the recent experiments [48]. Our scheme involving an all-optical feedback has a clear advantage in generating spatially separated macroscopic entangled states. In this work we have shown that using a cascaded coupling between the optical modes, the degree of quantum correlations in an OM network can be further enhanced. As shown in figures 2 and 3, using an all-optical feedback results in generating stronger entanglement between distant mechanical modes and is more robust to thermal fluctuations of the mechanical resonators. The advantage of using an all-optical feedback is clearly imprinted in the physical regime when the reversible coupling between the optical modes is small. We believe this is not a handicap for our scheme and in fact provides more flexibility in choosing the separation between distant nodes of an OM network. For an alternate topology considered and as shown in figure 6, we again find that an all-optical feedback may help in generating quantum correlations between distant optical and mechanical modes. We, however, do note that for the range of parameters considered in figure 6 steady state logarithmic negativity is one order of magnitude smaller as compared to figures 2 and 3. This brings us closer to make a short remark on the physical meaning of non-zero value of logarithmic negativity and a possible way to enhance entanglement between the mechanical modes.

As is well established now quantum entanglement is a vital requirement for the realization of various tasks including quantum cryptography, quantum metrology and quantum computing. However, entanglement is not very robust to environment induced decoherence and thus generally requires entanglement distillation as an auxiliary tool to counteract the degradation of coherence. Entanglement distillation refers to extracting a small number of maximally entangled states from a larger ensemble of weakly entangled states [57]. Logarithmic negativity, on the other hand, provides an upper-bound on the amount of distillable entanglement [58]. Thus one can make use of continuous variable entanglement distillation protocol to increase the entanglement of the shared state [59]. It should be mentioned that for distilling Gaussian continuous variable states, non-Gaussian operations including photon addition/photon subtraction needs to be employed [60]. An alternative approach is based on distilling Gaussian entanglement using quantum memory [61].

The ability to entangle a distributed array of mechanical elements will provide a path to new kinds of quantum enabled sensors. The Gaussian entanglement between harmonic oscillators that we have discussed in this paper is completely analogous to Gaussian entanglement for multi mode light fields. The enhanced metrology schemes that have been proposed in that case, for example [62], can carry over to the case of mechanical modes.

Since it is comparatively easier to distill and detect quantum correlations between optical modes, as compared to directly detecting quantum entanglement between mechanical modes, it has been suggested to swap the non-local correlations from the mechanical modes back to the optical modes [34, 37, 44]. As shown in figure 9 this can, for instance, be implemented using two auxiliary light modes, each initially prepared in classical uncorrelated states. These auxiliary modes can be two modes of distant cavities, and the geometry so arranged that each entangled mirror couples independently to the two modes. The non-local correlations may then be transferred from the movable mirrors to the initially uncorrelated auxiliary modes, which may eventually become entangled. Quantum correlated optical modes can then be eventually distilled. Thus, using standard homodyne measurement techniques, the entire correlation matrix of the two optical auxiliary modes can be reconstructed. A presence of non-zero quantum correlations between these optical modes will be an indirect signature of non-zero quantum correlations between the mechanical modes. The setup shown in figure 9 thus can also be used to ascertain the entanglement between the mechanical modes.

Zoom In Zoom Out Reset image size

In this work we have outlined a possibility of generating remote macroscopic entanglement. The possibility of generating such macroscopic superpositions could shed light on understanding the mysteries of nature [1] and well be a useful resource in the construction of long-distance quantum communication networks [3]. The model of OM network discussed in this work is general and we believe the ideas presented here could potentially be tested in number of physical settings. The parameters we have considered to describe the all-optical feedback scheme in this work do not place any stringent requirements on the degree of OM coupling strength nor the mechanical damping. Hence the ideas proposed in this paper would be suitable to a number of existing OM systems [63–65]. Comparing specific parameters from Chan et al [65] where mechanical damping γ = 7.5 kHz for cavity damping rate Γ/2π = 214 MHz and OM coupling strength g = 1.1 MHz fit in the regime of parameters we have considered in our calculations. A suitable system to realize our scheme would be the implementation of cavity optomechanics using ultra-cold atoms as discussed in [64]. Mechanics assisted coupling between the optical modes as required in the topology shown in figure 5 can be achieved using collective vibrational degree of freedom of an ensemble of ultra-cold atoms [64]. Another interesting test-bed for our results could be provided by optomechanical crystal (OMC) array which is a periodic structure that comprises both a photonic and a phononic crystal [66]. Coherent coupling between the photons and phonons in the crystal provides the ingredients for the OM interaction discussed in this work.

Other than creating massive quantum mechanical superposition of optical and mechanical modes, OM interaction provides other interesting avenues for implementing various tasks of quantum information processing. Storing and on-demand retrieval of quantum light, without compromising with its quantum character, has attracted wide attention both in the fields of quantum optics and quantum information processing [67, 68]. When it comes to implementing quantum light memory, electromagnetically induced transparency (EIT) is a much favored technique [69], where one can retain the large and highly desirable nonlinear optical properties associated with the resonant response of a material system. In analogy to the EIT in the atomic medium, optomechanically induced transparency has been experimentally demonstrated in [70]. Moreover as illustrated in [66], an optical waveguide can be coupled to an OMC array and quantum state of the light from the waveguide can be coherently transferred onto long-lived mechanical vibrations of the OMC array. OM transducers provides another promising application where coherent interaction between optics and mechanics can be exploited. As proposed in [3], an OM network can be employed to mediate interactions between distant nodes of a quantum network. Solid state qubits or electronic spin/charge degrees of freedom could form nodes of the network. Coupling can be achieved between (a) the evanescent field of the microcavity and motion of the mechanical resonator and (b) also between local nodes and vibrational degree of freedom of the mechanical resonator. Thus, indirect interaction mediated via mechanical resonator can then be set between optical mode and a distant node of the quantum network.

To summarize we have analyzed the influence of an all-optical feedback loop on the steady state dynamics of OM arrays. We considered different topologies forming OM networks and found that entanglement can be distributed over the chain of modes in the steady state and at finite temperatures. This opens up an interesting possibility to study spatially separated massive Schrödinger cat states. In each case our results demonstrate that all-optical feedback can be used to enhance and thus control entanglement between inter cavity mechanical and optical modes. For the OM network resulting from two optical modes with a common mechanical element, we note that the presence of feedback in fact generates entanglement between the optical mode a₂ and mechanics b which otherwise remain unentangled. We have also extended the analysis to the case of an OM network comprised of N-OM ports. We have shown that an all-optical feedback also helps in distributing stronger entanglement between inter cavity mechanical cavity modes. To conclude, we believe that the experimental demonstration of the quantum nature of macroscopic mechanical objects would help us to test long-standing questions about macroscopic quantum coherence and long-lived mechanical states.

CJ acknowledges support from the ORS scheme and kind hospitality offered by the Centre for Engineered Quantum Systems, University of Queensland where this work was initiated. UA acknowledges support from the University of Queensland postdoctoral fellowship and grant. This work was supported by the Australian Research Council grant CE110001013.

To solve the master equation (6) we make a Gaussian ansatz for the quantum characteristic function of the following form $\chi (\epsilon ,\eta ,x,y,{ t}) =\exp [-z^{\mathrm {T}}{{\bf{ A}}}(t)\,z+\mathrm {i} z^{\mathrm {T}}h(t) ]$, where A(t) is a time-dependent 8 × 8 matrix and h(t) is a 1 × 8 time-dependent vector and z^T = (,*,η,η*,x,x*,y,y*). The corresponding partial differential equation for χ(,η,x,y,t) then becomes

$$\begin{equation} \frac{\partial }{\partial t} \chi(\epsilon,\eta,x,y,{t}) =z^{\mathrm{T}} {{\bf N}} z\chi(\epsilon,\eta,x,y,{t}) +z^{\mathrm{T}} {{\bf M}} \nabla \chi(\epsilon,\eta,x,y,{t}), \end{equation} \tag{ A.1 }$$

where $\nabla = \left (\frac {\partial }{ \partial \epsilon },\frac {\partial }{ \partial \epsilon ^{*}},\frac {\partial }{ \partial \eta },\frac {\partial }{ \partial \eta ^{*}},\frac {\partial }{\partial x},\frac {\partial }{ \partial x^{*}},\frac {\partial }{ \partial y},\frac {\partial }{ \partial y^{*}}^{\mathrm {T}}\right )$

$$\begin{eqnarray}\fl {{\bf N}}&=&\left( \begin{array}{@{}cccccccc@{}} 0 & 0 & \mathrm{i}g_{1}/2 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & -\mathrm{i}g_{1}/2 & 0 & 0 & 0 & 0\\ \mathrm{i}g_{1}/2 & 0 & 0 & -\gamma \bar{n}/2 & 0 & 0 & 0 & 0\\ 0 & -\mathrm{i}g_{1}/2 & -\gamma \bar{n}/2 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\gamma \bar{n}/2\\ 0 & 0 & 0 & 0 & 0 & 0 & -\gamma \bar{n}/2 & 0\\ \end{array} \right ),\\ \fl {{\bf M}}&=&\left( \begin{array}{@{}cccccccc@{}} -\Gamma_{1}/2 & 0 & 0 & -\mathrm{i}g_{1} & \mathrm{i} \kappa & 0 & 0 & 0\\ 0 & -\Gamma_{1}/2 & \mathrm{i}g_{1} & 0 & 0 & -\mathrm{i} \kappa & 0 & 0\\ 0 & -\mathrm{i}g_{1} & -\gamma /2 & 0 & 0 & 0 & 0 & 0\\ \mathrm{i}g_{1} & 0 & 0 & -\gamma /2 & 0 & 0 & 0 & 0\\ \mathrm{i} \kappa -\sqrt{\Gamma_{1} \Gamma_{2}} & 0 & 0 & 0 & -\Gamma_{2}/2 & 0 & \mathrm{i} g_{2} & 0\\ 0 & -\mathrm{i} \kappa -\sqrt{\Gamma_{1} \Gamma_{2}} & 0 & 0 & 0 & -\Gamma_{2}/2 & 0 & -\mathrm{i}g_{2}\\ 0 & 0 & 0 & 0 & \mathrm{i}g_{2} & 0 & -\gamma/2 & 0\\ 0 & 0 & 0 & 0 & 0 & -\mathrm{i}g_{2} & 0 & -\gamma/2\\ \end{array} \right ).\nonumber \end{eqnarray} \tag{ A.2 }$$

Using the Gaussian ansatz for the quantum characteristic function χ(,η,x,y,t), it easily follows that

$$\begin{equation} \frac{\partial \chi }{ \partial t}=-z^{\rm T} \frac{{\rm d} {{\bf A}}}{{\rm d} t} z \chi+\mathrm{i} z^{\rm T} \frac{{\rm d} h}{{\rm d} t} \chi, \nabla \chi =-{\rm 2} {{\bf A}} z \chi +\mathrm{i }h \chi. \end{equation} \tag{ A.3 }$$

Using (A.3), the partial differential equation (A.1) for χ becomes

$$\begin{equation} -z^{\rm T} \frac{{\rm d} {{\bf A}}}{{\rm d} t} z\chi+\mathrm{i} z^{\rm T} \frac{{\rm d} h}{{\rm d} t} \chi = z^{\rm T} {{\bf N}} z \chi-{\rm 2}z^{\rm T}{{\bf M}}{{{\bf A}}} z \chi +\mathrm{i} z^{\rm T} {{\bf M}} h \chi. \end{equation} \tag{ A.4 }$$

Recalling that A(t) is symmetric, the symmetric part of the (A.4) results in following two matrix differential equations:

$$\begin{equation} -\frac{{\rm d} {{\bf A}}(t)}{{\rm d} t}=-{{\bf M}}{{\bf A}}-{{\bf A}} {{\bf M}}^{\rm T}+{{\bf N}};\, \frac{{\rm d} h}{{\rm d} t}= {\bf M} h. \end{equation} \tag{ A.5 }$$

The above coupled matrix differential equations can now be numerically solved to get the time evolved quantum characteristic function for the coupled OM cavities. For an initial Gaussian state evolving according to master equation (A.5), it is sufficient to fully characterize the quantum correlations between various optical and mechanical modes in terms of their Wigner covariance matrix. The covariance matrix V is a 8 × 8 real symmetric matrix V_i,j = (〈R_iR_j + R_jR_i〉)/2 where i,j{a,b} and $R^{\mathrm {T}}=(\hat {q}_{a}, \hat {p}_{a},\hat {q}_{b},\hat {p}_{b})$. Here $\hat {q}_{i}$ and $\hat {p}_{i}$ are the position and momentum quadratures of the ith mode. From the expression of the quantum characteristic function it is straightforward to extract the Wigner covariance matrix as follows, $\langle {\hat {b}_{1}}^{\dagger m} {\hat {b}_{2}}^{n} \rangle = (\frac {\partial }{\partial \eta })^{m}(-\frac {\partial }{\partial y^{*}})^{n}\chi (\epsilon =0,x=0,{t})|_{\eta =0, y=0}$.

In this section we study in more detail how feedback affects the steady state correlations between the mechanical modes b₁ and b₂. We focus our attention to a regime where each cavity mode has a large line width and can thus be adiabatically eliminated. An equivalent way of describing the dynamics of coupled OM network array is to write the following coupled Langevin equations:

$$\begin{equation} \frac{\rm d}{{\rm d}t}\hat{a}_{1}=-\mathrm{i}g_{1}\hat{b}_{1}^{\dagger}-\mathrm{i}\kappa\hat{a}_{2}-\frac{\Gamma}{2}\hat{a}_{1}+\sqrt{\Gamma}\hat{a}_{in}(1,t), \end{equation} \tag{ B.1 }$$

$$\begin{equation} \frac{\rm d}{{\rm d}t}\hat{a}_{2}=-\mathrm{i}g_{2}\hat{b}_{2}-\mathrm{i}\kappa\hat{a}_{1}-\frac{\Gamma}{2}\hat{a}_{2}-\sqrt{\Gamma^{2}}\hat{a}_{1}+\sqrt{\Gamma}\hat{a}_{in}(1,t-\tau),\end{equation} \tag{ B.2 }$$

$$\begin{equation} \frac{\rm d}{{\rm d}t}\hat{b}_{1}=-\mathrm{i}g_{1}\hat{a}_{1}^{\dagger}-\frac{\gamma}{2}\hat{b}_{1}+\sqrt{\gamma}\hat{b}_{in}(1,t),\end{equation} \tag{ B.3 }$$

$$\begin{equation} \frac{\rm d}{{\rm d}t}\hat{b}_{2}=-\mathrm{i}g_{2}\hat{a}_{2}-\frac{\gamma}{2}\hat{b}_{2}+\sqrt{\gamma}\hat{b}_{in}(2,t), \end{equation} \tag{ B.4 }$$

where for simplifying the calculations we have assume that γ₁ = γ₂ = γ and Γ₁ = Γ₂ = Γ. Neglecting the time delay between the source and the driven cavity τ and adiabatically eliminating both the cavity modes $\hat {a}_{1}$ and $\hat {a}_{2}$, Langevin equations for the mechanical modes become

$$\begin{eqnarray*} &&\fl \frac{\rm d}{{\rm d}t}\hat{b}_{1}^{\dagger}=\left(-\frac{\gamma}{2}+\frac{\Gamma}{2}\frac{g_{1}^{2}}{\sigma^{*}}\right)\hat{b}_{1}^{\dagger}+\sqrt{\gamma}\hat{b}_{ in}^{\dagger}(1,t)+\left(\mathrm{i}\frac{g_{1}}{2 \sigma^{*}}\Gamma^{3/2}+g_{1}\frac{\kappa \sqrt{\Gamma}}{\sigma^{*}}\right)\hat{a}_{in}-\mathrm{i} \hat{g}_{1}\hat{g}_{2}\frac{\kappa}{\sigma^{*}}\hat{b}_{2},\\&&\fl \frac{\rm d}{{\rm d}t}\hat{b}_{2}=\left(-\frac{\gamma}{2}-\frac{\Gamma}{2}\frac{g_{2}^{2}}{\sigma^{*}}\right)\hat{b}_{2}+\sqrt{\gamma}\hat{b}_{in}(2,t)+g_{1}g_{2}\frac{(\sqrt{\Gamma^{2}}+\mathrm{i}\kappa)}{\sigma^{*}}\hat{b}_{1}^{\dagger} \\&&+\left(\mathrm{i}g_{2}\sqrt{\Gamma}\frac{(\sqrt{\Gamma^{2}}+\mathrm{i}\kappa)}{\sigma^{*}}-\mathrm{i}g_{2}\frac{\Gamma^{3/2}}{2\sigma^{*}}\right)\hat{a}_{in}, \end{eqnarray*}$$

where $\sigma =\Gamma ^{2}/4+\mathrm {i}\sqrt {\Gamma ^{2}}\kappa +\kappa ^{2}$.

After a lengthy but otherwise straightforward exercise it is possible to obtain the steady state correlators between the two mechanical modes. Assuming that both the cavity modes are in contact with a zero temperature reservoir then the only non-zero correlator between the mechanical modes takes the form

$$\begin{equation} \langle \hat{b}_{1} \hat{b}_{2} \rangle_{\rm feedback=0} =\frac{4\mathrm{ i} \kappa g_{1}g_{2}}{\alpha^2}\gamma\left(\Gamma (g_{1}^2+g_{2}^2) \frac{2 \beta_{2}}{\beta_{1}}+\alpha\right)\frac{\beta_{2}}{\beta_{1}^{2}-4\beta_{2}^{2}}, \end{equation} \tag{ B.5 }$$

$$\begin{eqnarray}&&\fl \langle \hat{b}_{1} \hat{b}_{2} \rangle_{\rm feedback \ne 0} =\frac{4 (\sqrt{\Gamma^{2}}+\mathrm{i} \kappa) g_{1}g_{2}}{\alpha^2-\mathrm{i}16 \kappa \sqrt{\Gamma^{2}} g_{1}^{2}g_{2}^{2}}\gamma\left(\!\Gamma (g_{1}^2+g_{2}^2) \frac{2 \tilde{\beta}_{2}}{\tilde{\beta}_{1}} +\sqrt{\alpha^2-\mathrm{i}16 \kappa \sqrt{\Gamma^{2}} g_{1}^{2}g_{2}^{2}}\right)\frac{\tilde{\beta}_{2}}{\tilde{\beta}_{1}^{2}-4\tilde{\beta}_{2}^{2}},\nonumber\\ \end{eqnarray} \tag{ B.6 }$$

where

$$\begin{equation} \beta_{1}=\gamma+2\Gamma\frac{(g_{2}^{2}-g_{1}^{2})}{\Gamma^{2}+4\kappa^{2}},\end{equation} \tag{ B.7 }$$

$$\begin{equation} \beta_{2}=\frac{\alpha}{\Gamma^{2}+4\kappa^{2}}, \end{equation} \tag{ B.8 }$$

$$\begin{equation} \tilde{\beta}_{1}=\gamma+2\Gamma\frac{(g_{2}^{2}-g_{1}^{2})}{\Gamma^{2}+4\kappa^{2}-\mathrm{i}4\kappa\sqrt{\Gamma^{2}}},\end{equation} \tag{ B.9 }$$

$$\begin{equation} \tilde{\beta}_{2}=\frac{\sqrt{\alpha^2-\mathrm{i}16 \kappa \sqrt{\Gamma^{2}} g_{1}^{2}g_{2}^{2}}}{\Gamma^{2}+4\kappa^{2}-\mathrm{i}4\kappa\sqrt{\Gamma^{2}}},\end{equation} \tag{ B.10 }$$

$$\begin{equation} \alpha=\sqrt{\Gamma^{2}g_{1}^{4}+2(\Gamma^{2}+8\kappa^{2})g_{1}^{2}g_{2}^{2}+\Gamma^{2}g_{2}^{4}}. \end{equation} \tag{ B.11 }$$

As done before, to solve the master equation (14) we again convert it into a partial differential equation for the quantum characteristic function. Following the same steps as outlined in the previous section we again define a normal-ordered characteristic function $\chi (\epsilon ,\eta ,x,{t})=\langle {\mathrm { e}}^{\epsilon \hat {a}_{1}^{\dagger }}\,{\mathrm { e}}^{-\epsilon ^{*}\hat {a}_{1}}\,{\mathrm { e}}^{\eta \hat {a}_{2}^{\dagger }}\,{\mathrm { e}}^{-\eta ^{*}\hat {a}_{2}}\,{\mathrm { e}}^{x \hat {b}^{\dagger }}\,{\mathrm { e}}^{-x^{*}\hat {b}}\rangle$. Now getting an expression for the time evolved quantum characteristic function reduces to solving the following coupled matrix differential equations

$$\begin{equation} \frac{\partial }{\partial t} \chi(\epsilon,\eta,x,{t}) =z^{\mathrm{T}} {{\bf N}} z\chi(\epsilon,\eta,x,{t}) +z^{\mathrm{T}} {{\bf M}} \nabla \chi(\epsilon,\eta,x,{t}), \end{equation} \tag{ C.1 }$$

where $\nabla = (\frac {\partial }{ \partial \epsilon },\frac {\partial }{ \partial \epsilon ^{*}},\frac {\partial }{\partial x},\frac {\partial }{ \partial x^{*}},\frac {\partial }{ \partial \eta },\frac {\partial }{ \partial \eta ^{*}})^{\mathrm {T}}$

$$\begin{eqnarray} {{\bf N}}&=&\left( \begin{array}{@{}cccccc@{}} 0 & 0 & \mathrm{i} g_{1}/2 & 0 & 0 & 0\\ 0 & 0 & 0 & -\mathrm{i} g_{1}/2 & 0 & 0\\ \mathrm{i} g_{1}/2 & 0 & 0 & -\gamma_{1}\bar{n}/2 & 0 & 0\\ 0 & -\mathrm{i} g_{1}/2 & -\gamma_{1}\bar{n}/2 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ \end{array} \right ),\\ \nonumber {{\bf M}}&=&\left( \begin{array}{@{}cccccc@{}} -\Gamma_{1}/2 & 0 & 0 & -\mathrm{i}g_{1} & 0 & 0\\ 0 & -\Gamma_{1}/2 & \mathrm{i}g_{1} & 0 & 0 & 0\\ 0 & -\mathrm{i}g_{1} & -\gamma_{1}/2 & 0 & \mathrm{i}g_{2} & 0\\ \mathrm{i}g_{1} & 0 & 0 & -\gamma_{1}/2 & 0 & -\mathrm{i} g_{2}\\ -\sqrt{\Gamma_{1} \Gamma_{2}} & 0 & \mathrm{i}g_{2} & 0 & -\Gamma_{2}/2 & 0\\ 0 & -\sqrt{\Gamma_{1} \Gamma_{2}} & 0 & -\mathrm{i}g_{2} & 0 & -\Gamma_{2}/2\\ \end{array} \right ). \end{eqnarray} \tag{ C.2 }$$