Two-Membrane Cavity Optomechanics: Linear and Non-Linear Dynamics

Abstract

In this paper, we review the linear and non-linear dynamics of an optomechanical system made of a two-membrane etalon in a high-finesse Fabry–Pérot cavity. This two-membrane setup has the capacity to modify on demand the single-photon optomechanical coupling, and in the linearized interaction regime to cool simultaneously two mechanical oscillators. It is a promising platform for realizing cavity optomechanics with multiple resonators. In the non-linear regime, an analytical approach based on slowly varying amplitude equations allows us to derive a consistent and full characterization of the non-linear displacement detection, enabling a truthful detection of membrane displacements much above the usual linear sensing limited by the cavity linewidth. Such a high quality system also shows a pre-synchronization regime.

1. Introduction

Cooperative effects, enhanced interactions and nontrivial dynamics occur when multiple mechanical resonators are placed within an optical cavity [1,2,3,4,5,6,7,8,9,10]; for example, one can induce and control the coherent exchange of excitations [11,12,13,14], or study self-oscillations and their synchronization in the case of two or more mechanical resonators [12,15,16,17,18,19,20,21,22,23,24]. In this paper, we review the linear and non-linear dynamics of an optomechanical system made of a two-membrane etalon in a high-finesse Fabry–Pérot cavity [8,10,23]. When the optical cavity is driven on the red sideband, the linear dynamics of such a system is explored: the optomechanical coupling can be controlled on demand [8,9,25] by a local control of the membrane position along the cavity axis, and multiple oscillators can be simultaneously cooled [8,26], or exploited for photon-mediated coherent interaction and heat transfer between separate resonators [13,14]. However, the radiation pressure interaction is proportional to the photon number and it may have non-linear effects on both the mechanical and optical degrees of freedom which become evident when the mechanical motion is excited [27] by means of laser driving on the blue sideband of the optical cavity. Optical backaction in this case counteracts the internal mechanical friction, and when the total effective damping becomes equal to zero, a Hopf bifurcation into a regime of self-induced mechanical oscillations takes place [23,28,29,30,31,32,33,34]. A fixed amplitude limit cycle is established, with a free running oscillation phase, which may lock to external forces or to other optomechanical oscillators [35], leading to synchronization (see references [12,15,16,22,36,37,38,39,40,41,42] for theoretical characterizations, and references [17,18,19,20,21,24,43,44,45,46] for experimental demonstrations in optomechanical and electromechanical devices). In the specific case of the two-membrane-in-the-middle setup of interest here, self-organized synchronization, phase-locking, and the transition between in-phase and antiphase regimes have been qualitatively demonstrated [24], without calibration of the mechanical displacements. An alternative way of describing the non-linear optomechanical dynamics in this mechanical parametric oscillation regime is that the modulation of the radiation pressure interaction induced by the mechanical motion causes a cavity frequency shift comparable to or larger than the optical linewidth, yielding a nontrivial modification of the cavity response.

In this review article we first focus on studying the dynamics of the linearized fluctuations around the stable stationary state of the system, showing the tunability of the single-photon optomechanical coupling rate, and also that both membranes of the etalon can be simultaneously cooled by means of a red-detuned cavity driving. Later on we explore a regime, which can be called pre-synchronization regime, obtained when the blue-detuned driving is weak and slightly below the synchronization threshold. In this regime only one of the two membrane resonators, “master”, is driven to a limit cycle through the Hopf bifurcation, while the other oscillator is only partially synchronized with the “master” resonator because the amplitude of the synchronized component does not prevail over the thermal motion. We show that when multiple mechanical resonators are detected by the same single probe field simultaneously interacting with all of them (such as, for example, in references [17,18,19,43,45]), and at least one resonator enters a limit cycle, one has a non-trivial, non-linear dynamics of the system, which has to be properly considered, yielding a highly non-linear calibration of the displacement measurement obtained by means of the output probe readout.

This review article is organized as follows: in Section 2 we briefly describe the multimode optomechanical system under study, exploring both the linearized dynamics of the fluctuation within the stability conditions, and the non-linear regime, which is obtained when the Hopf bifurcation is crossed, and one has mechanical self-oscillations. In Section 3 we describe the experimental setup, and in Section 4 we present the linear dynamics for the membrane resonators, in particular their controllable coupling and the possibility to cool both of them. In Section 5 we analyze the non-linear dynamics of the mechanical modes at the onset of synchronization, and we provide the analytical recipe able to describe in a quantitative way both the full numerical solution of the Langevin equations, and the experimental results for the optical probe spectrum. Section 6 is for concluding remarks.

2. Theoretical Description of the System Dynamics

The following Hamiltonian describes an optomechanical system where two optical modes interact with two mechanical modes via radiation pressure:

$$\begin{array}{c}\hfill H=H_{\mathrm{pump}}+H_{\mathrm{probe}}+H_{\mathrm{mech}}+H_{\mathrm{int}},\end{array}$$

(1)

decomposed as follows:

$$\begin{array}{cc}\hfill H_{\mathrm{pump}}& =\hslash {\omega }_{c1}{a}_1^{†}{a}_1+\mathrm{i}\hslash E_1\left(a_1^{†}{\mathrm{e}}^{-\mathrm{i}{\omega }_{L1}t}-a_1{\mathrm{e}}^{\mathrm{i}{\omega }_{L1}t}\right),\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill H_{\mathrm{probe}}& =\hslash {\omega }_{c2}{a}_2^{†}{a}_2+\mathrm{i}\hslash E_2\left(a_2^{†}{\mathrm{e}}^{-\mathrm{i}{\omega }_{L2}t}-a_2{\mathrm{e}}^{\mathrm{i}{\omega }_{L2}t}\right),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill H_{\mathrm{mech}}& =\sum _{j=1,2}\hslash {\omega }_j{b}_j^{†}{b}_j,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill H_{\mathrm{int}}& =-\sum _{i,j=1,2}\hslash g_{ij}{a}_i^{†}{a}_i(b_j+b_j^{†}).\hfill \end{array}$$

(5)

Two cavity modes with resonance frequencies, ${\omega }_{c1}$ and ${\omega }_{c2}$, described by two bosonic annihilation operators $a_1$, and $a_2$, are driven at frequencies ${\omega }_{L1}$, and ${\omega }_{L2}$, respectively. The mode $a_1$, called PUMP, is used for engineering the optomechanical interaction, while $a_2$, denominated PROBE, performs a continuous detection of the mechanical motion. To prevent any interference between the two drivings they usually have a different cavity mode from each other (different polarization or frequency). $E_i=\sqrt{2{\kappa }_{in,i}{P}_i/\hslash {\omega }_{Li}}$ includes the driving dissipation rates, with ${\kappa }_{in,i}$ the decay rate of the i-th cavity mode through the input mirror, and $P_i$ the corresponding laser power.

The Hamiltonian $H_{\mathrm{mech}}$ describes two membrane oscillators with bosonic annihilation mechanical operator $b_j$ ($j=1,2$), each with mass $m_j$ and resonance frequency ${\omega }_j$, and displacement $q_j=(b_j^{†}+b_j)x_{\mathrm{zpf},j}$, where $x_{\mathrm{zpf},j}={(\hslash /2{\omega }_j{m}_j)}^{1/2}$ is the zero point motion of the j-th oscillator. The radiation pressure dispersive interaction couples the mechanical modes with the probe and pump optical modes. The single-photon coupling rates $g_{ij}=-(d{\omega }_{ci}/d{x}_j)x_{\mathrm{zpf},j}$ quantify the optomechanical interaction (see Figure 1).

We consider for both probe and pump modes, the frame rotating at the corresponding laser driving frequency. Fluctuation–dissipation processes couple the cavity modes and the resonators to their related thermal reservoir at temperature T, and they are included in the Heisenberg picture adding noise and dissipative terms. The quantum Langevin equations [27,47,48], for $i,j=1,2$, are:

$$\begin{array}{cc}\hfill \dot{a_i}& =(-\mathrm{i}{\Delta }_i^{\left(0\right)}-{\kappa }_i)a_i+E_i+\mathrm{i}\sum _{j=1,2}{g}_{ij}(b_j+b_j^{†})a_i+\sqrt{2{\kappa }_{in,i}}{a}_{in,i}+\sqrt{2{\kappa }_{ex,i}}{a}_{ex,i},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \dot{b_j}& =(-{\gamma }_j-\mathrm{i}{\omega }_j)b_j+\mathrm{i}\sum _{i=1,2}{g}_{ij}{a}_i^{†}{a}_i+\sqrt{2{\gamma }_j}{b}_j^{in},\hfill \end{array}$$

(7)

where ${\Delta }_i^{\left(0\right)}={\omega }_{ci}-{\omega }_{Li}$, ${\kappa }_i={\kappa }_{in,i}+{\kappa }_{ex,i}$ represents the amplitude decay rate of the cavity for the pump and the probe modes, ${\kappa }_{ex,i}$ is the optical loss rate through all the ports different from the input one, and ${\gamma }_j$ represents the amplitude decay rate of the j-th oscillator. $a_{in,i}\left(t\right)$, $a_{ex,i}\left(t\right)$ and $b_j^{in}$ are the corresponding zero-mean noise reservoir operators. They are all Gaussian, and white noises uncorrelated from each other with correlation functions $\langle f{\left(t\right)}^{†}f\left(t^{\prime }\right)\rangle =\overline{n}\delta (t-t^{\prime })$ and $\langle f\left(t\right)f{\left(t^{\prime }\right)}^{†}\rangle =(\overline{n}+1)\delta (t-t^{\prime })$ where $f\left(t\right)$ is either $a_{in,i}\left(t\right)$, $a_{ex,i}\left(t\right)$ or $b_j^{in}$, and $\overline{n}={[exp(\hslash {\omega }_f/k_{b}T)-1]}^{-1}$ is the mean thermal occupation number for the related mode.

We shall restrict, in this review, to study the dynamics at temperature $T\simeq 300$ K only, when a classical treatment of the above quantum Langevin equations is reasonable. At this temperature a different treatment of mechanical and optical noise terms should be considered: the optical frequencies ${\omega }_{c,i}/2\pi \simeq {10}^{14}$ Hz, so that ${\overline{n}}_{c,i}\simeq 0$, are dominated by photon shot noise, but for large enough driving powers fluctuations of the intracavity field might occur, due either to vacuum fluctuations or technical laser noise. At mechanical frequencies ${\omega }_1/2\pi \simeq {\omega }_2/2\pi \simeq {10}^6$ Hz, implying ${\overline{n}}_{m,j}\simeq k_{b}T/\hslash {\omega }_j\gg 1$, the thermal noise is dominant for the mechanical mode. For this reason, we consider classical complex random noises, ${\beta }_j^{in}\left(t\right)$, $j=1,2$, which replace the mechanical quantum thermal noise $b_j^{in}\left(t\right)$, and ${\alpha }_i^{opt}\left(t\right)$, which replace the sum of optical vacuum noises $\sqrt{{\kappa }_{in,i}/{\kappa }_i}{a}_{in,i}\left(t\right)+\sqrt{{\kappa }_{ex,i}/{\kappa }_i}{a}_{ex,i}\left(t\right)$, with correlation functions [39,42]

$$\begin{array}{cc}\hfill \langle {\beta }_j^{in}\left(t\right){\beta }_{j^{\prime }}^{in}\left(t^{\prime }\right)\rangle & =\langle {\alpha }_i^{opt}\left(t\right){\alpha }_{i^{\prime }}^{opt}\left(t^{\prime }\right)\rangle =0,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \langle {\beta }_{j^{\prime }}^{in}\left(t^{\prime }\right){\beta }_j^{in,*}\left(t\right)\rangle & =\langle {\beta }_j^{in,*}\left(t\right){\beta }_{j^{\prime }}^{in}\left(t^{\prime }\right)\rangle =({\overline{n}}_m+1/2){\delta }_{j{j}^{\prime }}\delta (t-t^{\prime }),\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \langle {\alpha }_i^{opt}\left(t^{\prime }\right){\alpha }_{i^{\prime }}^{opt,*}\left(t\right)\rangle & =\langle {\alpha }_i^{opt,*}\left(t\right){\alpha }_{i^{\prime }}^{opt}\left(t^{\prime }\right)\rangle =(1/2){\delta }_{i{i}^{\prime }}\delta (t-t^{\prime }).\hfill \end{array}$$

(10)

Therefore, the following coupled classical Langevin equations well approximates the set of coupled quantum Langevin equations for the optical and mechanical complex amplitudes ${\alpha }_i\left(t\right)$ and ${\beta }_j\left(t\right)$ [38,39,42],

$$\begin{array}{cc}\hfill {\dot{\alpha }}_i\left(t\right)& =\left(-\mathrm{i}{\Delta }_i^{\left(0\right)}-{\kappa }_i\right){\alpha }_i\left(t\right)+E_i+\sum _{j=1,2}2\mathrm{i}{g}_{ij}\mathrm{Re}\left[{\beta }_j\left(t\right)\right]{\alpha }_i\left(t\right)+\sqrt{2{\kappa }_i}{\alpha }_i^{opt}\left(t\right),\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\dot{\beta }}_j\left(t\right)& =(-\mathrm{i}{\omega }_j-{\gamma }_j){\beta }_j\left(t\right)+\mathrm{i}\sum _{i=1,2}{g}_{ij}{\left|{\alpha }_i\left(t\right)\right|}^2+\sqrt{2{\gamma }_j}{\beta }_j^{in}\left(t\right).\hfill \end{array}$$

(12)

In our case, the quasi-resonant weak probe beam is used only for the detection of the mechanical motion of the oscillators. We are far from the regime where the drivings relative phase can be used to control nonreciprocal effects [49,50].

2.1. Linear Dynamics

It is often useful to introduce the linearized approximate description of Equations (11) and (12). In the following we remove the index 1 when referring to the cavity pump field, that is ${\alpha }_1\to \alpha$, ${\kappa }_1\to \kappa$. We also ignore the Langevin equation of the weak probe beam, which is in resonance with the optomechanical cavity, and therefore realizes a perfectly non-invasive detection of the mechanical mode, without any backaction, as it occurs in a Michelson interferometer readout [34]. The cavity field and the mechanical mode are split into an average coherent amplitude and a fluctuating term $\alpha \left(t\right)={\alpha }_0$ + $\delta \alpha \left(t\right)$ and ${\beta }_j\left(t\right)={\beta }_{j,0}$ + $\delta {\beta }_j\left(t\right)$. In the steady state case, the average amplitudes are:

$$\begin{array}{cc}\hfill {\alpha }_0& =\frac{E}{\kappa +\mathrm{i}\left[{\Delta }^{\left(0\right)}-{\displaystyle \sum _{j=1,2}}{g}_j\left({\beta }_{0,j}+{\beta }_{0,j}^{*}\right)\right]},\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\beta }_{0,j}& =\frac{\mathrm{i}{g}_j{\left|{\alpha }_0\right|}^2}{{\gamma }_j+\mathrm{i}{\omega }_j}.\hfill \end{array}$$

(14)

Neglecting the second order terms, the fluctuation terms $\delta \alpha \left(t\right)$ and $\delta {\beta }_j\left(t\right)$ give the following linearized Langevin equations [27,51]:

$$\begin{array}{cc}\hfill \delta \dot{\alpha }\left(t\right)& =-\left(\mathrm{i}\Delta +\kappa \right)\delta \alpha \left(t\right)+\mathrm{i}\sum _{j=1,2}{g}_j[\delta {\beta }_j\left(t\right)+\delta {\beta }_j^{*}\left(t\right)]{\alpha }_0+\sqrt{2\kappa }\delta {\alpha }^{opt}\left(t\right),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \delta {\dot{\beta }}_j\left(t\right)& =-(\mathrm{i}{\omega }_j+{\gamma }_j)\delta {\beta }_j\left(t\right)+\mathrm{i}{g}_j[{\alpha }_0^{*}\delta \alpha \left(t\right)+{\alpha }_0\delta {\alpha }^{*}\left(t\right)]+\sqrt{2{\gamma }_j}{\beta }_j^{in}\left(t\right),\hfill \end{array}$$

(16)

where $\Delta ={\Delta }^{\left(0\right)}-{\sum }_{j=1,2}{g}_j({\beta }_{0,j}+{\beta }_{0,j}^{*})$ corresponds to a modification of the cavity detuning, which is the parameter controlled in the experiments. For our purpose, it is easier to solve the linearized Langevin equations in the frequency space:

$$\begin{array}{cc}\hfill -\mathrm{i}\omega \delta \tilde{\alpha }\left(\omega \right)& =-\left(\mathrm{i}\Delta +\kappa \right)\delta \tilde{\alpha }\left(\omega \right)+\mathrm{i}\sum _{j=1,2}{g}_j{\alpha }_0\left[\delta {\tilde{\beta }}_j\left(\omega \right)+\delta {\tilde{\beta }}_j^{*}\left(\omega \right)\right]+\sqrt{2\kappa }\delta {\tilde{\alpha }}^{opt}\left(\omega \right),\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill -\mathrm{i}\omega \delta {\tilde{\beta }}_j\left(\omega \right)& =-(\mathrm{i}{\omega }_j+{\gamma }_j)\delta {\tilde{\beta }}_j\left(\omega \right)+\mathrm{i}{g}_j\left[{\alpha }_0^{*}\delta \tilde{\alpha }\left(\omega \right)+{\alpha }_0\delta {\tilde{\alpha }}^{*}\left(\omega \right)\right]+\sqrt{2{\gamma }_j}\delta {\tilde{\beta }}_j^{in}\left(\omega \right),\hfill \end{array}$$

(18)

where $\delta {\tilde{\beta }}_j\left(\omega \right)={\int }_{-\infty }^{+\infty }dt{\mathrm{e}}^{\mathrm{i}\omega t}\delta {\beta }_j\left(t\right)$ and $\delta \tilde{\alpha }\left(\omega \right)={\int }_{-\infty }^{+\infty }dt{\mathrm{e}}^{\mathrm{i}\omega t}\delta \alpha \left(t\right)$ are the Fourier transform of $\delta {\beta }_j\left(t\right)$ and $\delta \alpha \left(t\right)$. Performing the calculations, and defining the optical and mechanical susceptibility ${\tilde{\chi }}_c\left(\omega \right)={\left[\kappa +\mathrm{i}\left(\Delta -\omega \right)\right]}^{-1}$ and ${\tilde{\chi }}_{m,j}\left(\omega \right)={\left[{\gamma }_j+\mathrm{i}\left({\omega }_j-\omega \right)\right]}^{-1}$, respectively, we obtain:

$$\begin{array}{cc}\hfill {\tilde{\chi }}_c^{-1}\left(\omega \right)\delta \tilde{\alpha }\left(\omega \right)& =\mathrm{i}\sum _{j=1,2}{g}_j{\alpha }_0\left[\delta {\tilde{\beta }}_j\left(\omega \right)+\delta {\tilde{\beta }}_j^{*}\left(\omega \right)\right]+\sqrt{2\kappa }\delta {\tilde{\alpha }}^{opt}\left(\omega \right),\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\tilde{\chi }}_{m,j}^{-1}\left(\omega \right)\delta {\tilde{\beta }}_j\left(\omega \right)& =\mathrm{i}{g}_j\left[{\alpha }_0^{*}\delta \tilde{\alpha }\left(\omega \right)+{\alpha }_0\delta {\tilde{\alpha }}^{*}\left(\omega \right)\right]+\sqrt{2{\gamma }_j}\delta {\tilde{\beta }}_j^{in}\left(\omega \right).\hfill \end{array}$$

(20)

To obtain the modified susceptibility for the mechanical oscillator we solve the previous system of equations inserting Equation (19) into (20). Defining $\Sigma \left(\omega \right)={\mathrm{i}\left|\alpha \right|}^2\left[{\tilde{\chi }}_c^{*}\left(\omega \right)-{\tilde{\chi }}_c\left(\omega \right)\right]$, and $\delta {\tilde{u}}_j^{opt}=g_j\left[{\alpha }_0^{*}\delta {\tilde{\alpha }}^{opt}{\tilde{\chi }}_c\left(\omega \right)+{\alpha }_0\delta {\tilde{\alpha }}^{opt{}^{*}}{\tilde{\chi }}_c^{*}\left(\omega \right)\right]$, we get the linearized Langevin equations for the two mechanical oscillators in the rotating wave approximation, that is neglecting the counter-rotating terms $\delta {\tilde{\beta }}_j^{*}$ [52], which is valid in the red-detuned and resolved-sideband regime (even moderate):

$$\begin{array}{cc}\hfill \left[{\tilde{\chi }}_{m,1}^{-1}+\mathrm{i}{g}_1^2\Sigma \left(\omega \right)\right]\delta {\tilde{\beta }}_1\left(\omega \right)& =\mathrm{i}{g}_1{g}_2\Sigma \left(\omega \right)\delta {\tilde{\beta }}_2\left(\omega \right)+\mathrm{i}\sqrt{2\kappa }\delta {\tilde{u}}_1^{opt}+\sqrt{2{\gamma }_1}\delta {\tilde{\beta }}_1^{in}\left(\omega \right),\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \left[{\tilde{\chi }}_{m,2}^{-1}+\mathrm{i}{g}_2^2\Sigma \left(\omega \right)\right]\delta {\tilde{\beta }}_2\left(\omega \right)& =\mathrm{i}{g}_1{g}_2\Sigma \left(\omega \right)\delta {\tilde{\beta }}_1\left(\omega \right)+\mathrm{i}\sqrt{2\kappa }\delta {\tilde{u}}_2^{opt}+\sqrt{2{\gamma }_2}\delta {\tilde{\beta }}_2^{in}\left(\omega \right).\hfill \end{array}$$

(22)

Neglecting the optical noise terms (high temperature regime), we finally obtain:

$$\begin{array}{cc}\hfill \left[{\tilde{\chi }}_{m,1}^{eff-1}+\frac{{\Sigma }_1\left(\omega \right){\Sigma }_2\left(\omega \right)}{{\tilde{\chi }}_{m,2}^{eff-1}}\right]\delta {\tilde{\beta }}_1\left(\omega \right)& \approx \sqrt{2{\gamma }_1}\delta {\tilde{\beta }}_1^{in}\left(\omega \right),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \left[{\tilde{\chi }}_{m,2}^{eff-1}+\frac{{\Sigma }_1\left(\omega \right){\Sigma }_2\left(\omega \right)}{{\tilde{\chi }}_{m,1}^{eff-1}}\right]\delta {\tilde{\beta }}_2\left(\omega \right)& \approx \sqrt{2{\gamma }_2}\delta {\tilde{\beta }}_2^{in}\left(\omega \right),\hfill \end{array}$$

(24)

where ${\tilde{\chi }}_{m,j}^{eff-1}={\tilde{\chi }}_{m,j}^{-1}+\mathrm{i}{\Sigma }_j\left(\omega \right)$, and ${\Sigma }_j\left(\omega \right)=g_j^2\Sigma \left(\omega \right)$. From the previous equations, we can derive the effective susceptibility of the two mechanical oscillators driven by a red detuned pump beam [48]:

$$\begin{array}{c}\hfill {\tilde{\chi }}_{m,j}^{eff}=\frac{{\tilde{\chi }}_{m,j}^{eff}}{1+{\Sigma }_1\left(\omega \right){\tilde{\chi }}_{m,1}^{eff}{\Sigma }_2\left(\omega \right){\tilde{\chi }}_{m,2}^{eff}}.\end{array}$$

(25)

2.2. Non–Linear Dynamics

In the non-linear regime, the chaotic motion of the two oscillators is not considered, as it occurs at extremely large driving powers which are not physically meaningful for the optomechanical system discussed in this review. We might find the dynamics of the system by considering the slowly varying amplitude equations approach of references [16,22]. After an initial transient regime, each mechanical oscillator sets itself into the following dynamics:

$$\begin{array}{c}\hfill {\beta }_j\left(t\right)={\beta }_{0,j}+A_j\left(t\right){\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathrm{ref}}t},\end{array}$$

(26)

where ${\beta }_{0,j}$ is the approximately constant, static shift of the j-th resonator, $A_j\left(t\right)$ is the corresponding slowly–varying complex amplitudes, and ${\omega }_{\mathrm{ref}}\gg \Delta \omega$ is a reference mechanical frequency, of the order of ${\omega }_j$. Equation (26) implies that we will study the long–time dynamics of the two mechanical resonators in the frame rotating at the fast reference frequency ${\omega }_{\mathrm{ref}}$. Inserting the Equation (26) into the Equation (11), and performing explicitly the integrals one gets:

$$\begin{array}{cc}\hfill {\alpha }_i\left(t\right)& ={\mathrm{e}}^{\mathrm{i}{\psi }_i\left(t\right)}{\int }_0^{t}d{t}^{\prime }{\mathrm{e}}^{{\mathcal{W}}_i(t-t^{\prime })}[E_i+\sqrt{2{\kappa }_i}{\alpha }_i^{opt}\left(t^{\prime }\right)]{\mathrm{e}}^{-\mathrm{i}{\psi }_i\left(t^{\prime }\right)}\equiv \hfill \\ \hfill & \equiv {\alpha }_i^E\left(t\right)+\delta {\alpha }_i\left(t\right),\hfill \end{array}$$

(27)

where the intracavity field ${\alpha }_i^E\left(t\right)$ is proportional to the driving rate $E_i$, and $\delta {\alpha }_i\left(t\right)$ is related to the input noise ${\alpha }_i^{opt}\left(t\right)$. We defined the bright complex amplitudes $A_i^b\left(t\right)=\left|A_i^b\left(t\right)\right|{\mathrm{e}}^{\mathrm{i}{\theta }_i\left(t\right)}=\sum g_{ij}{A}_j\left(t\right)/g_i^b$ with $g_i^b={(g_{i1}^2+g_{i2}^2)}^{1/2}$, and ${\mathcal{W}}_i=\mathrm{i}[{\Delta }_i^{\left(0\right)}-\sum g_{ij}({\beta }_{0,j}+{\beta }_{0,j}^{*})]-{\kappa }_i$, and ${\psi }_i\left(t\right)={\xi }_{i}sin({\omega }_{\mathrm{ref}}t-{\theta }_i)$ with ${\xi }_i=2g_i^b\left|A_i^b\right|/{\omega }_{\mathrm{ref}}$.

By performing the Jacobi–Anger expansion for the ${\mathrm{e}}^{-\mathrm{i}{\psi }_i\left(t^{\prime }\right)}$ factor within the integral [16,30], that is, ${\mathrm{e}}^{-\mathrm{i}\xi sin\varphi }={\sum }_n{J}_n(-\xi ){\mathrm{e}}^{\mathrm{i}n\varphi }$, ($\varphi ={\omega }_{\mathrm{ref}}{t}^{\prime }-{\theta }_i$, and $J_n$ is the n-th Bessel function of the first kind), and neglecting a transient decay term, the expression for the intracavity field amplitude ${\alpha }_i^E\left(t\right)$ is:

$$\begin{array}{c}\hfill {\alpha }_i^E\left(t\right)=E_i{\mathrm{e}}^{\mathrm{i}{\psi }_i\left(t\right)}\sum _{n=-\infty }^{\infty }\frac{J_n\left(-{\xi }_i\right){\mathrm{e}}^{\mathrm{i}n({\omega }_{\mathrm{ref}}t-{\theta }_i)}}{\mathrm{i}n{\omega }_{\mathrm{ref}}-{\mathcal{W}}_i}.\end{array}$$

(28)

For the fluctuation term we have:

$$\begin{array}{c}\hfill \delta {\alpha }_i\left(t\right)=\sqrt{2{\kappa }_i}{\mathrm{e}}^{\mathrm{i}{\psi }_i\left(t\right)}{\int }_0^{t}d{t}^{\prime }{\mathrm{e}}^{{\mathcal{W}}_i(t-t^{\prime })}{\alpha }_i^{opt}\left(t^{\prime }\right).\end{array}$$

(29)

We might derive an equation for the unknown quantities $A_j\left(t\right)$ and ${\beta }_{0,j}$ by inserting these relations into the radiation pressure force term within the Equation (12) for the mechanical motion. Since the intracavity optical fluctuations are negligible, one can approximate at first order the radiation pressure term in $\delta {\alpha }_i\left(t\right)$,

$$\begin{array}{c}\hfill \mathrm{i}{g}_{ij}|{\alpha }_i{\left(t\right)|}^2\simeq \mathrm{i}{g}_{ij}{\left|{\alpha }_i^E\left(t\right)\right|}^2+\mathrm{i}{g}_{ij}{\eta }_i^{opt}\left(t\right),\end{array}$$

(30)

where

$$\begin{array}{c}\hfill |{\alpha }_i^E{\left(t\right)|}^2=E_i^2\sum _{n,m=-\infty }^{\infty }\frac{J_n\left(-{\xi }_i\right)J_m\left(-{\xi }_i\right){\mathrm{e}}^{\mathrm{i}(n-m)({\omega }_{\mathrm{ref}}t-{\theta }_i)}}{(\mathrm{i}n{\omega }_{\mathrm{ref}}-{\mathcal{W}}_i)(-\mathrm{i}m{\omega }_{\mathrm{ref}}-{\mathcal{W}}_i^{*})},\end{array}$$

(31)

and

$$\begin{array}{c}\hfill {\eta }_i^{opt}\left(t\right)={\alpha }_i^E\left(t\right)\delta {\alpha }_i^{*}\left(t\right)+{\alpha }_i^{E,*}\left(t\right)\delta {\alpha }_i\left(t\right).\end{array}$$

(32)

Assuming ${\beta }_{0,j}$ approximately constant, and neglecting all the terms that oscillate faster than ${\omega }_{\mathrm{ref}}$, i.e., keeping only the resonant terms in Equation (31) [$n-m=0$ for ${\beta }_{0,j}$ and $n-m=-1$ for the amplitudes $A_j\left(t\right)$], we get

$$\begin{array}{c}\hfill ({\gamma }_j+\mathrm{i}{\omega }_j){\beta }_{0,j}=\mathrm{i}\sum _{i=1,2}\sum _{n=-\infty }^{\infty }\frac{g_{ij}{E}_i^2{J}_n{\left(-{\xi }_i\right)}^2}{(\mathrm{i}n{\omega }_{\mathrm{ref}}-{\mathcal{W}}_i)(-\mathrm{i}n{\omega }_{\mathrm{ref}}-{\mathcal{W}}_i^{*})},\end{array}$$

(33)

which is an implicit equation for ${\beta }_{0,j}$ because ${\mathcal{W}}_i$ depends upon ${\beta }_{0,j}$, and the value of ${\beta }_{0,j}$ can be obtained numerically. For the slowly varying amplitudes $A_j\left(t\right)$ we get instead:

$$\begin{array}{cc}\hfill {\dot{A}}_j\left(t\right)& =\left[-{\gamma }_j-\mathrm{i}\Delta {\omega }_j\right]A_j\left(t\right)+\sqrt{2{\gamma }_j}{\beta }_j^{in}\left(t\right)+\mathrm{i}\sum _{i=1,2}{g}_{ij}{\eta }_i^{opt}\left(t\right)+\hfill \\ \hfill & +\mathrm{i}\sum _{i=1,2}\sum _{n=-\infty }^{\infty }\frac{g_{ij}{\mathrm{e}}^{\mathrm{i}{\theta }_i}{E}_i^2{J}_n\left(-{\xi }_i\right)J_{n+1}\left(-{\xi }_i\right)}{[\mathrm{i}n{\omega }_{\mathrm{ref}}-{\mathcal{W}}_i][-\mathrm{i}(n+1){\omega }_{\mathrm{ref}}-{\mathcal{W}}_i^{*}]},\hfill \end{array}$$

(34)

where $\Delta {\omega }_j={\omega }_j-{\omega }_{\mathrm{ref}}$, and

$$\begin{array}{c}\hfill {\Delta }_i^{\left(0\right)}\to {\Delta }_i={\Delta }_i^{\left(0\right)}-\sum _{j=1,2}{g}_{ij}({\beta }_{0,j}+{\beta }_{0,j}^{*}),\end{array}$$

(35)

correspond to a modification of the cavity detunings, as already explained in the previous Section, so that ${\mathcal{W}}_i=-\mathrm{i}{\Delta }_i-{\kappa }_i$ can be regarded as a given parameter when we calculate the amplitude.

Finally, the Equation (34) can be rewritten in a more transparent form by making explicit the equation for each amplitude, and by defining the following regular dimensionless auxiliary functions ${\mathcal{F}}_i\left(\right|A_i^b|,{\kappa }_i,{\Delta }_i)$, $i=1,2$, as:

$$\begin{array}{c}\hfill {\mathcal{F}}_i=\frac{E_i^2}{|A_i^b|}\sum _{n=-\infty }^{\infty }\frac{J_n\left(-{\xi }_i\right)J_{n+1}\left(-{\xi }_i\right)}{[\mathrm{i}n{\omega }_{\mathrm{ref}}-{\mathcal{W}}_i][-\mathrm{i}(n+1){\omega }_{\mathrm{ref}}-{\mathcal{W}}_i^{*}]}\prime \end{array}$$

(36)

which can be easily shown to be a function of even powers of $|A_i^b|$ only. One finally gets the set of coupled slowly-varying complex amplitude equations:

$$\begin{array}{cc}\hfill {\dot{A}}_1\left(t\right)& =\left[-{\gamma }_1-\mathrm{i}(\Delta {\omega }_1-d_1)\right]A_1\left(t\right)+\mathrm{i}{d}_{12}{A}_2\left(t\right)+\mathrm{i}\sum _{i=1,2}{g}_{i1}{\eta }_i^{opt}\left(t\right)+\sqrt{2{\gamma }_1}{\beta }_1^{in}\left(t\right),\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill {\dot{A}}_2\left(t\right)& =\left[-{\gamma }_2-\mathrm{i}(\Delta {\omega }_2-d_2)\right]A_2\left(t\right)+\mathrm{i}{d}_{12}{A}_1\left(t\right)+\mathrm{i}\sum _{i=1,2}{g}_{i2}{\eta }_i^{opt}\left(t\right)+\sqrt{2{\gamma }_2}{\beta }_2^{in}\left(t\right),\hfill \end{array}$$

(38)

where

$$\begin{array}{cc}\hfill d_1& =\left(\frac{g_{11}^2{\mathcal{F}}_1}{g_1^b}+\frac{g_{21}^2{\mathcal{F}}_2}{g_2^b}\right),\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill d_2& =\left(\frac{g_{12}^2{\mathcal{F}}_1}{g_1^b}+\frac{g_{22}^2{\mathcal{F}}_2}{g_2^b}\right),\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill d_{12}& =\left(\frac{g_{11}{g}_{12}{\mathcal{F}}_1}{g_1^b}+\frac{g_{21}{g}_{22}{\mathcal{F}}_2}{g_2^b}\right),\hfill \end{array}$$

(41)

are non-linear coefficients because of their dependence upon the regular dimensionless auxiliary functions ${\mathcal{F}}_i$, which, in turn, depend upon the corresponding variable $|A_i^b{|}^2={|g_{i1}{A}_1+g_{i2}{A}_2|}^2/(g_{i1}^2+g_{i2}^2)$. As already shown in references [16,22], Equations (37) and (38) give an accurate and general description of two mechanical resonators dynamics.

3. Experimental Setup

In this Section we describe the experimental setup used to study the linear and non-linear dynamics of an optomechanical system constituted of a two-membrane etalon inserted in a high-finesse Fabry–Pérot cavity. The membrane-etalon is realised with two equal $1.5\mathrm{m}\mathrm{m}\times 1.5\mathrm{m}\mathrm{m}$ membranes (Norcada), with a etalon length $L_{\mathrm{c}}=53.571\left(9\right) \mathsf{\mu }\mathrm{m}$. The membranes have a thickness $L_{\mathrm{m}}=106\left(1\right) \mathrm{n}\mathrm{m}$, considering the index of refraction of Si${}_3$N${}_4$ given in reference [53]. The membrane-sandwich is placed at the center of a high-finesse Fabry–Pérot cavity. Inserting the membrane-etalon, the empty cavity finesse ${\mathcal{F}}_0$ = 50,125 (25) decreases to $\mathcal{F}$ = 12,463 (13), corresponding to an intensity decay rate of the cavity $2\kappa ={\tau }^{-1}=\mathrm{FSR}/\mathcal{F}\simeq 2\pi \times 134 \mathrm{kHz}$, with $\mathrm{FSR}\simeq 2\pi \times 1.67 \mathrm{GHz}$. The mechanical and optical properties of this system were explored in references [8].

Figure 2 shows the experimental setup. A laser beam at wavelength ${\lambda }_0=1064 \mathrm{n}\mathrm{m}$ is split into a pump beam and a probe beam. An acousto-optic modulator (AOM) shifts the frequency of the pump beam by ${\Delta }_1$ from the cavity resonance. The probe beam is in resonance with the optical cavity and has an intensity $P_{probe}=5.9 \mathsf{\mu }\mathrm{W}$. The light reflected by the optical cavity is homodyne detected and the thermal voltage spectral noise (VSN) is measured. The pump beam is used to laser drive the mechanical modes. In the following Sections, we focus on studying the linear and non-linear dynamics with the fundamental modes of the two membranes.

4. Experiments in the Linear Regime

In the linear regime, we first show that shifting the position of the membranes along the cavity axis with the piezo controllers the optomechanical interaction can be tuned and controlled. Figure 3a shows a simulation of the shift of the resonance frequencies of the optomechanical cavity as a function of the position of the two oscillators, $q_j$, along the standing wave resonating in the cavity. The corresponding gradient field is represented as superimposed vector plot, and shows how the two optomechanical couplings $G_j=\partial {\omega }_{cav}/\partial q_j=g_j/x_{\mathrm{zpf},j}$ change by displacing the membranes [8]. Note that there are positions $q_j$ for which only one of the membranes is coupled (red and magenta dots), or both (green dot).

To prove tunability of the optomechanical coupling rate, the probe beam was frequency locked to the optical cavity using the PDH technique, and the thermal VSN of the two membranes is measured with homodyne detection of the light reflected by the optical cavity (see Figure 2). Figure 3b shows the detected thermal VSN, which clearly demonstrates the capacity to switch off and on the optomechanical interaction in a controlled way by shifting the position of each membrane (see Figure 3b top and middle spectra, and also the magenta and red dots in Figure 3a) where only one oscillator is in a position in which it interacts with the cavity field. In the bottom spectrum of Figure 3b both membranes are instead coupled to the optical cavity [see also the green dot in Figure 3a]. We measured ${\omega }_{m1}=2\pi \times 235.810 \mathrm{kHz}$, ${\gamma }_{m1}=2\pi \times 1.64 \mathrm{Hz}$ for the lower frequency fundamental mode on the left (red–dashed line), while for the fundamental mode of the second resonator we measured ${\omega }_{m2}=2\pi \times 236.580 \mathrm{kHz}$, ${\gamma }_{m2}=2\pi \times 9.37 \mathrm{Hz}$ (orange–dashed line).

To avoid any optomechanical effect with the probe beam, such as cooling or optical spring effect, we used a very low power probe field as resonant as possible (${\Delta }_2\simeq 0$) with a cavity mode. The corresponding measured single photon optomechanical coupling rates $g_j=G_j{x}_{\mathrm{zpf},j}{\Theta }_j$, and ${\Theta }_j$, the dimensionless transverse overlap between the optical cavity mode and the j-th mechanical mode, are $g_1=2\pi \times 0.30\mathrm{Hz}$ and $g_2=2\pi \times 0.28\mathrm{Hz}$. These couplings are comparable to the ones obtained in a similar setup with a single membrane-in-the-middle [54,55]. With this setup, we can tune the optomechanical interaction of both resonators with the optical mode, switching on the “cooling” pump beam with a variable detuning ${\Delta }_1\equiv \Delta$ with respect the cavity resonance. In this Section we focus on the case of a red-detuned drive which enhances the beam-splitter interaction between the mechanical modes and the cavity mode allowing the cooling of the former, and whose dynamics is described by the treatment of Section 2.1. We observe the simultaneous cooling of the fundamental modes of the two distinct membranes [26], demonstrated by the significant suppression of the area below the measured displacement spectral noise (DSN) when the driving laser is resonant with the red sideband and the pump power is not negligible (see Figure 4). The behaviour of the DSN is determined by the effective susceptibility of each resonator, modified by their common interaction with the driven pump mode, and explicitly given by Equation (25), obtained in Section 2.1. In more detail, in Figure 4a,b we report the measured DSN as a function of the detuning $\Delta$ normalized to the mean mechanical frequency ${\overline{\omega }}_m$. In Figure 4a we use a lower power of the cooling beam with respect to that used in Figure 4b. We note the lower optomechanical coupling of the left mode for the results in Figure 4b, which corresponds to an ineffective optical cooling. In Figure 4c instead the DSN is shown as a function of the pump beam power $P_{pump}$, at a fixed detuning $\Delta \sim {\overline{\omega }}_m$.

5. Experiments in the Non-Linear Regime

In Figure 5 the VSN measured with homodyne signal is reported as a function of time. The frequencies are normalized with respect to the frequency of the fundamental mode of the first membrane, ${\omega }_1$, (underlined by an orange-dashed line, while the second membrane mode is marked by a green-dashed line). The parameters were set to drive the mechanical oscillators in a weak regime.

In the first 10 $\mathrm{s}$, we switched off the pump beam, and measured the thermal noise displacement of the fundamental modes of the two resonators. The magenta-dashed line highlights the external signal injected to estimate the single-photon optomechanical coupling $g_1\simeq 2\pi \times 0.43 \mathrm{Hz}$, and $g_2\simeq 2\pi \times 0.70 \mathrm{Hz}$, and to calibrate the VSN in displacement spectral noise (DSN) [34,56]. After 10 $\mathrm{s}$, the blue detuned (${\Delta }_1/2\pi =259.350 \mathrm{kHz}$) pump beam is switched on with a power of $4.25$ $\mathsf{\mu }$$\mathrm{W}$, to study the non-linear dynamics of the optomechanical system. After 25 $\mathrm{s}$ the power is further increased to $6.0$ $\mathsf{\mu }$$\mathrm{W}$. In the experimental data, we observe the appearance of a sideband due to the second mechanical mode at frequency $2{\omega }_1-{\omega }_2$.

The onset of synchronization is observed in the weak driving regime. Only one of the two oscillators enters into a limit cycle through the Hopf bifurcation associated with the parametric instability [30]. Instead, the other membrane stays in a mixed condition where the radiation pressure force induced by the oscillations of the first membrane cannot prevail over the thermal noise [16,22]. For convenience we take ${\omega }_{\mathrm{ref}}={\omega }_1$ as a reference in Equations (37) and (38), so that $\Delta {\omega }_1=0$ and $\Delta {\omega }_2=\Delta \omega$. As one of the two membranes remains in a thermal state, we make the assumption that $|A_1|\gg |A_2|,\sqrt{2{\overline{n}}_{m,1}}$, where ${\overline{n}}_j={\left[exp\left(\hslash {\omega }_j/k_{B}T\right)-1\right]}^{-1}\simeq k_{B}T/\hslash {\omega }_j\gg 1$ is the mean thermal occupation number. Moreover, we will not consider the optical noise and the terms associated with ${\eta }_i^{opt}\left(t\right)$, as they are negligible in our experiment. With the above approximations, both $A_2$ and thermal noise contributions can be neglected, and Equation (37) becomes:

$$\begin{array}{c}\hfill {\dot{A}}_1\left(t\right)=\left[-{\gamma }_1+\mathrm{i}{d}_1\left(\right|A_1\left|\right)\right]A_1\left(t\right)=-\left[{\gamma }_1^{eff}\left(\right|A_1\left|\right)-\mathrm{i}\Delta {\omega }_1^{eff}\left(\right|A_1\left|\right)\right]A_1\left(t\right),\end{array}$$

(42)

where the dependence of $d_1$ on $|A_1|$ has been made explicit, and we defined $\Delta {\omega }_1^{eff}\left(\right|A_1\left|\right)=-\mathrm{Re}[d_1\left(\right|A_1\left|\right)]$, and ${\gamma }_1^{eff}\left(\right|A_1\left|\right)={\gamma }_1+\mathrm{Im}[d_1\left(\right|A_1\left|\right)]$. The mechanical effective damping can be cast as:

$$\begin{array}{c}\hfill {\gamma }_1^{eff}\left(\right|A_1\left|\right)={\gamma }_1\left[1+\frac{g_1}{{\gamma }_1\left|A_1\right|}\mathrm{Im}\left[E_1^2{\Sigma }_1+E_2^2{\Sigma }_2\right]\right],\end{array}$$

(43)

where assuming $g_{i1}\simeq g_1$, in the considered regime $|A_i^b|\simeq |A_1|(g_{i1}/g_i^b)$, ${\xi }_j\simeq {\xi }_1=2g_1\left|A_1\right|/{\omega }_1$, and

$$\begin{array}{c}\hfill {\Sigma }_j\equiv \Sigma ({\xi }_1,{\kappa }_i,{\Delta }_i)=\sum _n\frac{J_n\left(-{\xi }_1\right)J_{n+1}\left(-{\xi }_1\right)}{[\mathrm{i}n{\omega }_1-{\mathcal{W}}_j][-\mathrm{i}(n+1){\omega }_1-{\mathcal{W}}_j^{*}]}.\end{array}$$

(44)

With such approximation we imply we are in a pre-synchronization regime where, in the second oscillator, the thermal noise prevails; if the second mode is able to arrive to a limit cycle and synchronize with the first, the amplitude $g_2{A}_2$ would not be negligible anymore with respect to $g_1{A}_1$, and the dynamics would be described by the general Equations (37) and (38) [24,57]. We can solve Equation (42), rewriting it in terms of modulus and phase, $A_1=I_1{\mathrm{e}}^{\mathrm{i}{\varphi }_1}$,

$$\begin{array}{cc}\hfill {\dot{I}}_1\left(t\right)& =-{\gamma }_1^{eff}\left(I_1\right)I_1\left(t\right),\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill {\dot{\varphi }}_1\left(t\right)& =\Delta {\omega }_1^{eff}\left(I_1\right).\hfill \end{array}$$

(46)

After a transient these equations give, for the first oscillator, a steady state solution with a constant radius, $I_1^{\mathrm{st}}$, which in our case of weak driving, corresponds to the smallest positive root of the equation ${\gamma }_1^{eff}\left(I_1^{\mathrm{st}}\right)=0$, which can be cast as:

$$\begin{array}{c}\hfill a|{\xi }_1^{\mathrm{st}}| =-\mathrm{Im}\left[E_1^2\Sigma ({\xi }_1^{\mathrm{st}},{\kappa }_1,{\Delta }_1)+E_2^2\Sigma ({\xi }_1^{\mathrm{st}},{\kappa }_2,{\Delta }_2)\right],\end{array}$$

(47)

with ${\xi }_1^{\mathrm{st}}=2g_1{I}_1^{\mathrm{st}}/{\omega }_1$, and $a={\omega }_1{\gamma }_1/2g_1^2$. As a consequence, at long times, ${\varphi }_1^{\mathrm{st}}\left(t\right)\simeq t\Delta {\omega }_1^{\mathrm{st}}$ with $\Delta {\omega }_1^{\mathrm{st}}=\Delta {\omega }_1^{eff}\left(I_1^{\mathrm{st}}\right)$ so that $A_1\left(t\right)\simeq A_1^{\mathrm{st}}\left(t\right)=I_1^{\mathrm{st}}exp[\mathrm{i}t\Delta {\omega }_1^{\mathrm{st}}]$.

In Figure 6 we show the left and right side of Equation (47). We infer from the intersection point of Figure 6b, which corresponds to find the smallest root of ${\gamma }_1^{eff}\left(I_1^{\mathrm{st}}\right)=0$, a value ${\xi }_1^{\mathrm{st}}=1.054$, a steady–state displacement amplitude $q_1^{\mathrm{st}}=2\left|A_1\right|x_{\mathrm{zpf}}=263.0 \mathrm{p}\mathrm{m}$, and $\Delta {\omega }_1^{eff}=-2\pi ·0.04\mathrm{Hz}$. Such behaviour has been investigated for the determination of $g_0$ [34]. For $\xi \ll 1$, the sideband output field has a spectral amplitude linear with $\xi$, and it is possible to perform a direct measurement of the position coordinate $q_1$; for $\xi \ge 1$ we should consider a correction factor because linearity is no more valid. In our case, the theoretical correction factor is ${\mathcal{N}}_1\simeq 0.70$, which corresponds to an expected observable stationary limit cycle amplitude of $q_1^{\mathrm{ob}}\simeq 183 \mathrm{p}\mathrm{m}$ (see Figure 7) [23]. Because of the oscillating behaviour of the Bessel functions, Equation (47) considering sufficiently large pump power, may have more than one solution [oblique black–dashed line in Figure 6b], which corresponds to the multistability phenomenon theoretically analysed, and then verified in references [30,32].

We now focus on the dynamics of the second resonator inserting the steady–state solution for $A_1\left(t\right)$ into Equation (38), which becomes:

$$\begin{array}{c}\hfill {\dot{A}}_2\left(t\right)=\left[-{\gamma }_2-\mathrm{i}\Delta \omega +\mathrm{i}{d}_2\left(I_1^{\mathrm{st}}\right)\right]A_2\left(t\right)+\mathrm{i}{d}_{12}\left(I_1^{\mathrm{st}}\right)A_1^{\mathrm{st}}\left(t\right)+\sqrt{2{\gamma }_2}{\beta }_2^{in}\left(t\right).\end{array}$$

(48)

The stationary solution is obtained performing the Fourier transform, and it is written as:

$$\begin{array}{c}\hfill A_2\left(t\right)=\frac{\mathrm{i}{d}_{12}\left(I_1^{\mathrm{st}}\right)}{{\gamma }_2^{eff}+\mathrm{i}\Delta {\overline{\omega }}_2^{eff}}{A}_1^{\mathrm{st}}\left(t\right)+\sqrt{2{\gamma }_2}{\int }_0^{t}ds{\mathrm{e}}^{-\left({\gamma }_2^{eff}+\mathrm{i}\Delta {\omega }_2^{eff}\right)s}{\beta }_2^{in}(t-s),\end{array}$$

(49)

where ${\gamma }_2^{eff}={\gamma }_2+\mathrm{Im}\left[d_2\left(I_1^{\mathrm{st}}\right)\right]$ is positive, that is, despite the pump driving, the second resonator is still damped, and does not reach a limit cycle, $\Delta {\overline{\omega }}_2^{eff}=\Delta \omega +\mathrm{Re}[d_1\left(I_1^{\mathrm{st}}\right)-d_2\left(I_1^{\mathrm{st}}\right)]$, and $\Delta {\omega }_2^{eff}=\Delta \omega -\mathrm{Re}\left[d_2\left(I_1^{\mathrm{st}}\right)\right]$. Therefore on the right hand side of Equation (49) the first term is the synchronized component which oscillates at the same frequency of the first master oscillator, while the second term is the thermal noise component at its natural frequency. This equation describes the synchronization of the second oscillator with the first one. When $|d_{12}\left(I_1^{\mathrm{st}}\right){|}^2{I}_1^{\mathrm{st},2}\gg {\gamma }_2{\gamma }_2^{\prime }{\overline{n}}_{m,2}$ (where $\Delta {\overline{\omega }}_2^{eff}\simeq \Delta {\omega }_2^{eff}\simeq \Delta \omega$), i.e., the thermal noise is negligible, the two resonators achieve phase locking and full synchronization. The case of synchronization is consistent with the theoretical analysis of references [16,22]. In our case of small driving powers, an onset of synchronization with very different limit cycle amplitudes is predicted, even in the presence of thermal noise. Considering the four dimensional phase space of the mechanical resonators, this condition is described by a Neimark–Sacker bifurcation, which corresponds to the birth of a stable torus around the existing limit cycle [35].

This analysis is validated considering the experimental time traces shown in Figure 7, where the displacement amplitudes of the two membranes $q_1$ and $q_2$ are determined by means of the calibration signal [34,56]. For $t<4 \mathrm{s}$ (the pump beam is off), a value of the displacement amplitudes higher than the thermal ones is observed because of a slightly blue-detuned probe beam. Such detuning is evaluated observing that, for the second mode, [green curve in panel b)], the standard deviation of the calibrated measured position $\Delta q_2^{\Delta }\simeq 3.50 \mathrm{p}\mathrm{m}$, while the estimated standard deviation of the thermal displacement is $\Delta q_2^{\mathrm{th}}\simeq 3.32 \mathrm{p}\mathrm{m}$, so that

$$\begin{array}{c}\hfill \frac{\Delta q_2^{\Delta }}{\Delta q_2^{\mathrm{th}}}=\frac{1}{\sqrt{1+C(\Delta )}}\simeq 1.054,\end{array}$$

(50)

where, considering a quasi-resonant field [27],

$$\begin{array}{c}\hfill C(\Delta )\sim -\frac{2g_j^2{E}^2}{\gamma \kappa }\frac{4{\omega }_j}{{({\kappa }^2+{\omega }_j^2)}^2}\Delta .\end{array}$$

(51)

We estimate a small blue detuning of the probe ${\Delta }_2\simeq 2\pi ·3.9 \mathrm{kHz}$.

Finally, after $t>4 \mathrm{s}$, the pump beam is switched on and the measured mechanical amplitudes are observed: we notice that as the amplitude of the first oscillator increases, the measured amplitude of the second membrane decreases below the thermal value. The observed steady-state limit cycle displacement amplitude of the first resonator is $q_1^{\mathrm{ob}}\simeq 184 \mathrm{p}\mathrm{m}$ [orange curve in Figure 7d]. The effective steady-state displacement amplitude of the first oscillator, $q_1^{\mathrm{st}}=2\left|A_1\right|x_{\mathrm{zpf}}\simeq 262 \mathrm{p}\mathrm{m}$ (blue curve of panel d), is obtained by the 10 blue trajectories simulated with the parameters of reference [23]. Moreover, the slope of the trajectories follows the measured one, suggesting that our approach of using the slowly varying complex amplitudes of the two resonators is effective, and able to catch all the properties of the non-linear dynamics.

We notice that our model is also able to describe the dynamics of the second mode with very good agreement. When the pump beam is switched on, the measured $q_2$ [the green curve in Figure 7b] follows the dynamics of the effective mechanical displacement [blue trajectories in Figure 7b]. After 7 $\mathrm{s}$$q_1$ goes into the limit cycle, and then the effective displacement differs from the measured one. An estimation of the effective steady-state amplitude of the second contribution in Equation (49) is provided by $\Delta q_2^{\mathrm{th}}\sqrt{{\gamma }_2/{\gamma }_2^{eff}}$, with

$$\begin{array}{c}\hfill \frac{{\gamma }_2^{eff}}{{\gamma }_2}=1+\frac{\mathrm{Im}\left[d_2\left(I_1^{\mathrm{st}}\right)\right]}{{\gamma }_2}=1+\frac{g_2^2{E}_1^2}{{\gamma }_2{I}_1{g}_1}\mathrm{Im}\left[{\Sigma }_1+\eta {\Sigma }_2\right]=1-\frac{{\gamma }_1{g}_2^2}{{\gamma }_2{g}_1^2}.\end{array}$$

(52)

In our case, we notice that the effective amplitude is larger than the thermal amplitude by a factor $\sqrt{{\gamma }_2/{\gamma }_2^{eff}}\simeq 1.38$, because there is a small effective driving, even if is not enough for the appearance of a limit cycle. In conclusion, we notice that the non-linear dynamics of the system affects also the small effective amplitude displacement of the second resonator, which shows a fictitious cooling effect, that is instead only a manifestation of a non-linear detection of the displacement amplitude in such a regime. When the amplitude of one oscillator brings a modulation larger than the cavity linewidth, the motional amplitudes detected by the probe beam are non-linearly modified and an appropriate correction factor has to be considered. This happens also for the unexcited resonator whose motional amplitude gives a much smaller cavity frequency modulation than the cavity linewidth.

6. Conclusions

The optomechanical performance of a Fabry–Pérot cavity containing a two-membrane etalon and driven in the red- and blue-detuning has been reviewed. We studied the linear regime of such a system, and experimentally demonstrated the ability to tune on demand the optomechanical coupling of the two resonators with the optical field, and showed the simultaneous optical cooling of the two distinct oscillators. Hereafter, we studied the non-linear regime where a pump beam blue-detuned from the resonance frequency of the optical cavity, brings one of the two oscillators into a self-sustaining limit cycle. We find out that the system is in a pre-synchronized regime when driven by a weak pump beam, where the unexcited oscillator has a synchronized, small component at the frequency of the excited oscillator. We get a perfect agreement between the numerical simulations, the experimental results, and an analytical approach based on slowly varying amplitude equations. Both large and small amplitude resonator motions are transduced in a nontrivial way in this non-linear detection regime.