Two-membrane cavity optomechanics: non-linear dynamics

Multimode optomechanical systems [1] are attracting an increasing interest for the study of collective dynamical effects, both at quantum and classical level. Two different situations are mainly considered from both the theoretical and experimental point of view: (i) a group of mechanical oscillators interacting via radiation pressure with the same optical mode [2–17] (e.g. multiple membranes within the same optical cavity); (ii) an array of mechanical oscillators each interacting locally with a single optical mode, and coupled by the tunnelling of photons and phonons between neighbouring sites [18–23], (e.g. optomechanical crystals in one and two dimensions [24]).

Several features of multimode optomechanical systems have already been investigated in the literature, such as long-range collective interactions [4, 7, 9] yielding an effective increase of the optomechanical coupling, slowing and stopping light [18], correlated quantum many-body states [20], reservoir engineering and dynamical phase transitions [6], graphene-like Dirac physics [22], topological phases of sound and light [23], transport in a one-dimensional chain [21, 25, 26], superradiance and collective gain [27], and nonreciprocal routing of electromagnetic signals [28, 29].

The radiation pressure interaction is inherently non-linear and the effects of such non-linearity on the mechanical motion are easily manifested when the optical cavity is driven on the blue sideband, when optical backaction is responsible for mechanical antidamping [1]. When the latter overcomes the internal mechanical friction, a Hopf bifurcation towards a regime of self-induced mechanical oscillations takes place [30–36], with a fixed amplitude, and a free running oscillation phase, which may lock to external forces or to other optomechanical oscillators [37]. This mutual phase-locking of self-oscillating resonators is at the basis of optomechanical synchronization, which has been thoroughly investigated both theoretically [5, 19, 20, 38–45], and experimentally [46–55] under different configurations. The non-linear effects of radiation pressure manifest themselves whenever the mechanical motion produces a cavity frequency shift comparable or larger than the optical linewidth, resulting in a nontrivial modification of the cavity response to the external driving. This is responsible for a variety of non-linear phenomena beyond synchronization, such as phonon lasing [56], mode competition [57], and chaos [58–60]. This radiation-pressure-induced non-linear behaviour may occur not only when the mechanical resonators are driven to large amplitude via the parametric amplification provided by blue-sideband driving, but also in the strong optomechanical coupling regime [61] where even intrinsic Brownian motion induces cavity frequency fluctuations larger than the optical linewidth [62, 63]. In both situations, the optomechanical non-linearity plays a fundamental role, affecting optomechanical displacement measurement and transduction, and this role can be exploited for extending in a nontrivial way the dynamic range of optomechanical sensors beyond the cavity linewidth regime [64].

Here we experimentally explore the non-linear dynamics of the multimode optomechanical setup first demonstrated in reference [14], realized by placing a membrane cavity within a high-finesse Fabry–Pérot cavity. Reference [14] reported a ∼2.47 gain in the optomechanical coupling strength of the membrane relative motion with respect to the single membrane case, and showed the capability to tune the single-photon optomechanical coupling on demand. Reference [55] recently demonstrated synchronization of this two-membrane cavity optomechanical system, by operating with a low-finesse cavity in the strongly unresolved sideband regime. Here instead we focus onto the pre-synchronization regime of weak blue-detuned driving, where only one of the two membrane resonators enters into a limit cycle through the Hopf bifurcation, while the other resonator remains in a mixed condition where the modulation of the radiation pressure force induced by the excited oscillator does not prevail over the thermal motion. We provide a detailed, quantitative analysis of the dynamics in this regime, with a significant agreement between the experimental data, the numerical simulation, and the analytical treatment based on amplitude equations of reference [45]. This quantitative analysis is based on a detailed treatment of the optical detection apparatus including the probe and calibration tones, and provides an accurate, reliable, measurement of the displacement of both membranes, even in the non-linear regime where the frequency modulation caused by the two membranes' motion is significantly larger than the cavity linewidth. A remarkable result of this analysis is that, in the presence of a self-oscillating resonator in a limit cycle, non-linear corrections to the displacement measurement by the probe cavity output must be applied not only to the excited resonator but also to the small-amplitude, unexcited one. This implies that in multimode optomechanical systems, whenever multiple mechanical resonators are detected by the same single probe field (such as for example in references [46, 47, 50, 51, 53]), and at least one resonator enters a limit cycle, one has to properly include the full non-linear dynamics of the system in order to extract the correct displacement measurement from the output probe spectrum.

The paper is organized as follows: in section 2 we provide the basic theoretical description of the multimode optomechanical system under study. In section 3 we describe the experimental setup, and in section 4 we derive in detail the probe beam power spectral density, including all the non-linear effects. In section 5 we analyse the non-linear dynamics of the mechanical modes at the onset of synchronization and we provide an analytical description in very good agreement with the numerical and experimental results. Section 6 is for concluding remarks.

We study the non-linear dynamics of a multimode optomechanical system, formed by two electromagnetic and two mechanical modes, at room temperature, which justifies a treatment in terms of classical amplitudes, and implies that thermal noise will be dominant for the mechanical modes and treated as classical complex random noises. The two optical modes with frequencies ω_ci , i = (1, 2), total cavity amplitude decay rates κ_i = κ_in,i + κ_ex,i with κ_ex,i optical loss rates through all the ports different from the input one κ_in,i , and driven at frequencies ω_Li , interact via radiation-pressure with two mechanical modes with resonance frequencies ω_j , j = (1, 2), mass m_j , and amplitude decay rates γ_j . Their dynamics is described by the set of coupled classical Langevin equations for the corresponding optical and mechanical complex amplitudes α_i (t) and β_j (t) [40, 41, 44], respectively,

$$\begin{equation}{\dot {\alpha }}_{i}(t)=\left(\mathrm{i}{{\Delta}}_{i}^{(0)}-{\kappa }_{i}\right){\alpha }_{i}(t) + {E}_{i} + \sum\limits _{j=1,2}2\mathrm{i}{g}_{ij}\enspace \text{Re}[{\beta }_{j}(t)]{\alpha }_{i}(t) + \sqrt{2{\kappa }_{i}}{\alpha }_{i}^{\mathrm{o}\mathrm{p}\mathrm{t}}(t),\end{equation} \tag{ 1 }$$

$$\begin{equation}{\dot {\beta }}_{j}(t)=(-\mathrm{i}{\omega }_{j}-{\gamma }_{j}){\beta }_{j}(t)+\mathrm{i}\sum\limits _{i=1,2}{g}_{ij}\vert {\alpha }_{i}(t){\vert }^{2}+\sqrt{2{\gamma }_{j}}{\beta }_{j}^{\mathrm{i}\mathrm{n}}(t),\end{equation} \tag{ 2 }$$

where ${{\Delta}}_{i}^{(0)}={\omega }_{Li}-{\omega }_{ci}$ are the detunings, ${E}_{i}=\sqrt{2{\kappa }_{\mathrm{i}\mathrm{n},i}{P}_{i}/\hslash {\omega }_{Li}}$ the driving rates with P_i the associated laser input powers, g_ij = −(dω_ci /dx_j )x_zpf,j the single-photon optomechanical coupling rates, ${x}_{\mathrm{z}\mathrm{p}\mathrm{f},\mathrm{j}}=\sqrt{\hslash /2{m}_{j}{\omega }_{j}}$ the spatial width of the jth oscillator zero point motion, and ${\beta }_{j}^{\mathrm{i}\mathrm{n}}(t)$, and ${\alpha }_{i}^{\mathrm{o}\mathrm{p}\mathrm{t}}(t)$, are the mechanical and optical noise terms, respectively. These noises are uncorrelated from each other and the only nonzero correlation functions are $\langle {\beta }_{j}^{\mathrm{i}\mathrm{n},{\ast}}(t){\beta }_{{j}^{\prime }}^{\mathrm{i}\mathrm{n}}({t}^{\prime })\rangle =({\bar{n}}_{j}+1/2){\delta }_{j{j}^{\prime }}\delta (t-{t}^{\prime })$, where ${\bar{n}}_{j}={\left[\mathrm{exp}\left(\hslash {\omega }_{j}/{k}_{\mathrm{B}}T\right)-1\right]}^{-1}\simeq {k}_{\mathrm{B}}T/\hslash {\omega }_{j}\gg 1$ is the mean thermal occupation number. Multimode optomechanical systems formed by two electromagnetic modes and two mechanical modes have been proposed and demonstrated [28, 29] for the nonreciprocal routing of signals controlled by the relative phase of multiple external and off-resonant drives. In the present case the weak, quasi-resonant driving probe beam is used only for detecting the mechanical motion and we are far from the regime where one can use and control the relative phase of the drivings for nonreciprocal effects. In this system, under appropriate parameter regimes, the pump cavity mode (i = 1) may drive the oscillators into a self-sustained limit cycle [30–35], which may eventually become synchronized. Synchronization may occur on a long timescale, determined by the inverse of the typically small parameters Δω = ω₂ − ω₁ (typically never larger than few kHz), and γ_j (order of Hz). Therefore it is physically useful to derive from the full dynamics of the classical Langevin equations (1) and (2), approximate equations able to correctly describe the slow, long-time dynamics of the two mechanical resonators, leading eventually to synchronization.

We adapt here the slowly varying amplitude equations approach of reference [5] to the case with noise studied here, as discussed in detail in reference [45], discarding here the limiting case of chaotic motion of the two resonators, which however occurs only at extremely large driving powers, and are not physically meaningful for the Fabry–Pérot cavity system considered here, it is known that each mechanical resonator, after an initial transient regime, sets itself into a dynamics of the following form

$$\begin{equation}{\beta }_{j}(t)={\beta }_{0,j}+{A}_{j}(t){\text{e}}^{-\text{i}{\omega }_{\mathrm{r}\mathrm{e}\mathrm{f}}t},\end{equation} \tag{ 3 }$$

where β_0,j is the approximately constant, static shift of the jth resonator, A_j (t) is the corresponding slowly-varying complex amplitudes, and ω_ref ≫ Δω is a reference mechanical frequency, of the order of ω_j .

From equations (1) and (2) one gets the set of coupled amplitude equations [45]

$$\begin{align}{\dot {A}}_{1}(t)=\left[-{\gamma }_{1}-\mathrm{i}{\Delta}{\omega }_{1}\right]{A}_{1}(t)& +\mathrm{i}{d}_{1}{A}_{1}(t)+\mathrm{i}{d}_{12}{A}_{2}(t)+\mathrm{i}\sum\limits _{i=1,2}{g}_{i1}{\eta }_{i}^{\mathrm{o}\mathrm{p}\mathrm{t}}(t)+\sqrt{2{\gamma }_{1}}{\beta }_{1}^{\mathrm{i}\mathrm{n}}(t),\end{align} \tag{ 4 }$$

$$\begin{align}{\dot {A}}_{2}(t)=\left[-{\gamma }_{2}-\mathrm{i}{\Delta}{\omega }_{2}\right]{A}_{2}(t)& +\mathrm{i}{d}_{2}{A}_{2}(t)+\mathrm{i}{d}_{12}{A}_{1}(t)+\mathrm{i}\sum\limits _{i=1,2}{g}_{i2}{\eta }_{i}^{\mathrm{o}\mathrm{p}\mathrm{t}}(t)+\sqrt{2{\gamma }_{2}}{\beta }_{2}^{\mathrm{i}\mathrm{n}}(t),\end{align} \tag{ 5 }$$

with Δω_j = ω_j − ω_ref, and where

$$\begin{equation}{d}_{1}=\left(\frac{{g}_{11}^{2}{\mathcal{F}}_{1}}{{g}_{1}^{b}}+\frac{{g}_{21}^{2}{\mathcal{F}}_{2}}{{g}_{2}^{b}}\right),\end{equation} \tag{ 6 }$$

$$\begin{equation}{d}_{2}=\left(\frac{{g}_{12}^{2}{\mathcal{F}}_{1}}{{g}_{1}^{b}}+\frac{{g}_{22}^{2}{\mathcal{F}}_{2}}{{g}_{2}^{b}}\right),\end{equation} \tag{ 7 }$$

$$\begin{equation}{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),\end{equation} \tag{ 8 }$$

are non-linear coefficients because of their dependence upon the regular dimensionless auxiliary functions ${\mathcal{F}}_{i}$, which are given by

$$\begin{equation}{\mathcal{F}}_{i}=\frac{{E}_{i}^{2}}{\vert {A}_{i}^{b}\vert }\sum\limits _{n=-\infty }^{\infty }\frac{{J}_{n}\left(-{\xi }_{i}\right){J}_{n+1}\left(-{\xi }_{i}\right)}{[\mathrm{i}n{\omega }_{\mathrm{r}\mathrm{e}\mathrm{f}}-{\mathcal{W}}_{i}][-\mathrm{i}(n+1){\omega }_{\mathrm{r}\mathrm{e}\mathrm{f}}-{\mathcal{W}}_{i}^{{\ast}}]},\end{equation} \tag{ 9 }$$

and which can be easily shown to be a function of even powers of $\vert {A}_{i}^{b}\vert$ only. We have defined the bright complex amplitudes ${A}_{i}^{b}(t)=\vert {A}_{i}^{b}(t)\vert {\text{e}}^{\text{i}{\theta }_{i}(t)}={\sum }_{j}{g}_{ij}{A}_{j}(t)/{g}_{i}^{b}$, ${g}_{i}^{b}=\sqrt{{g}_{i1}^{2}+{g}_{i2}^{2}}$, ${\xi }_{i}=2{g}_{i}^{b}\vert {A}_{i}^{b}\vert /{\omega }_{\mathrm{r}\mathrm{e}\mathrm{f}}$, ${\mathcal{W}}_{i}=\text{i}\left[{{\Delta}}_{i}^{(0)}+\sum {g}_{ij}({\beta }_{0,j}+{\beta }_{0,j}^{{\ast}})\right]-{\kappa }_{i}$, and J_n is the nth Bessel function of the first kind. ${\eta }_{i}^{\mathrm{o}\mathrm{p}\mathrm{t}}$ is the term describing the noise of optical origin [45]. As already shown in references [5, 45], equations (4) and (5) provide a general and very accurate description of the dynamics of the two mechanical resonators.

The experimental setup for studying the non-linear dynamics in an optomechanical system constituted of a two-membrane sandwich within a cavity, is shown in figure 1. A laser beam at wavelength λ₀ = 1064 nm is split in a probe beam with intensity P_probe = 5.9 μW, modulated by an electro-optical modulator (EOM), and a pump beam, detuned by Δ₁ from the cavity resonance by means of an acousto-optic modulator (AOM). The reflected probe beam is locked to the optical cavity by means of a Pound–Drever–Hall (PDH) technique, and the thermal voltage spectral noise (VSN) is measured by homodyne detection of the light reflected by the optical cavity. The pump beam is used for engineering the optomechanical interaction, and in particular to realise laser driving of the mechanical modes.

Zoom In Zoom Out Reset image size

The optical and mechanical properties of the optomechanical system were investigated in reference [14]. The membrane-cavity length, realised with two equal membranes (Norcada), was measured to be L_c = 53.571(9) μm, and the membrane thickness is L_m = 106(1) nm that is found assuming the index of refraction of Si₃N₄ given in reference [65]. Assuming rectangular membranes, and the nominal values provided by the manufacturer for the stress, σ = 0.825 GPa, and for the density ρ = 3100 kg m⁻³, the side lengths were estimated to be ${L}_{x}^{(1)}=$ 1.519(6) mm, ${L}_{y}^{(1)}=$ 1.536(6) mm, and ${L}_{x}^{(2)}=$ 1.522(6) mm, ${L}_{y}^{(2)}=$ 1.525(6) mm. We studied the dynamics of the lower frequency mode of the two membranes: for the first membrane we measured ω₁ ≃ 2π × 230.795 kHz, γ₁ ≃ 2π × 1.64 Hz, while for the second ω₂ ≃ 2π × 233.759 kHz, γ₂ ≃ 2π × 9.37 Hz. The membrane-cavity is placed in the middle of an optical cavity with empty cavity finesse ${\mathcal{F}}_{0}=$ 50125(25), which reduces to $\mathcal{F}=$ 12463(13) when the membrane-sandwich is placed in. Such finesse corresponds to a cavity intensity decay rate $2\kappa ={\tau }^{-1}=\mathrm{F}\mathrm{S}\mathrm{R}/\mathcal{F}\simeq 2\pi {\times}$ 134 kHz, with FSR ≃ 2π × 1.67 GHz.

Experimentally we have studied the weak driving regime of the onset of synchronization where only one of the two membrane resonators, say first oscillator, enters into a limit cycle through the Hopf bifurcation associated with the parametric instability [32]. In this regime the other resonator remains in a mixed condition where the modulation of the radiation pressure force induced by the first oscillator is not yet able to prevail over the thermal noise contribution [5, 45].

Zoom In Zoom Out Reset image size

In figure 2 are compared simulated and experimental results for a set of parameters in this weak driving regime. On panel (b) is reported the VSN of the homodyne signal as a function of time, where the frequencies are counted from the frequency of the fundamental mode of the first oscillator, ω₁, (marked by an orange square symbol, while the second mode is marked by a dark-green square symbol). During the first 10 s, the pump beam is turned off, and the VSN shows the thermal displacement of the fundamental modes of the two membranes. The magenta square symbol marks the external tone used for determining the single-photon optomechanical coupling g₁ ≃ 2π × 0.43 Hz, and g₂ ≃ 2π × 0.70 Hz, and for calibrating the VSN in displacement spectral noise (DSN) [66]. Finally, after 10 s the pump beam is turned on at a blue detuning Δ₁/2π = 259.350 kHz with a power of 4.25 μW for studying the dynamics of the optomechanical system, and after 25 s increased to 6.0 μW. Panel (c) of figure 2 shows the variances corresponding to the integral over the frequency range denoted by the squares symbols in panel (b). Panel (a) of figure 2 shows the numerical simulation of the non-linear dynamics of the system provided by integration of equations (1) and (2). We note the remarkable agreement between simulated and experimental data. Results in figure 2 shows that when the first membrane reaches a limit cycle, sidebands around its central frequency appears on the output field. In particular the experimental data and simulations demonstrate the appearance of the sideband due to the second mode (light-green square symbol) at 2ω₁ − ω₂.

For a detailed quantitative description of these results, we need to reconsider the non-linear dynamics of the mean cavity-field amplitude α₂ of the probe beam described by equation (1), and replace the input noise operator with the input field with a sinusoidal modulation at ω_b for calibration. Moreover, for simplicity we drop everywhere the index 2 when referring to the cavity probe field, that is α₂ → α, κ₂ → κ and so on. We get

$$\begin{equation}\dot {\alpha }=-\kappa \alpha +\mathrm{i}\left[{{\Delta}}^{(0)}+2{g}_{1}\enspace \text{Re}[{\beta }_{1}(t)]+2{g}_{2}\enspace \text{Re}[{\beta }_{2}(t)]\right]\alpha +{E}_{\mathrm{i}\mathrm{n}},\end{equation} \tag{ 10 }$$

where ${E}_{\mathrm{i}\mathrm{n}}=\sqrt{2{\kappa }_{\mathrm{i}\mathrm{n}}}\enspace {\text{e}}_{\mathrm{i}\mathrm{n}}\mathrm{exp}[-\mathrm{i}\beta \enspace \mathrm{sin}({\omega }_{b}t)]$, and ${\text{e}}_{\mathrm{i}\mathrm{n}}=\sqrt{{\mathcal{P}}_{\mathrm{i}\mathrm{n}}/\hslash {\omega }_{L}}$. We assume that the behaviour of the two oscillators is also sinusoidal as in equation (3), i.e., β_j (t) = β_0,j + |A_j |exp[i(ω_mj t + ϕ_j )], and that the first oscillator reaches a limit cycle with an amplitude |A₁| for which g₁|A₁| is much larger that g₂|A₂|. The solution of equation (10) can be found considering an expansion in terms of = g₂|A₂|/g₁|A₁|, that is α = ∑_j ^j α_j . The zero-order solution, α₀, satisfies

$$\begin{equation}{\dot {\alpha }}_{0}=-\kappa {\alpha }_{0}+\mathrm{i}\left[{\Delta}+2{g}_{1}\vert {A}_{1}\vert \mathrm{cos}({\omega }_{1}t+{\phi }_{1})\right]{\alpha }_{0}+\sqrt{2{\kappa }_{\mathrm{i}\mathrm{n}}}\enspace {\text{e}}_{\mathrm{i}\mathrm{n}}\enspace {\text{e}}^{-\text{i}\beta \mathrm{sin}({\omega }_{b}t)},\end{equation} \tag{ 11 }$$

where Δ = Δ⁽⁰⁾ + 2g₁ Re[β_0,1] + 2g₂ Re[β_0,2], and a first-order perturbation solution α₁, driven by the amplitude of the second oscillator,

$$\begin{equation}{\dot {\alpha }}_{1}=-\kappa {\alpha }_{1} + \mathrm{i}\left[{\Delta}+2{g}_{1}\vert {A}_{1}\vert \mathrm{cos}({\omega }_{1}t + {\phi }_{1})\right]{\alpha }_{1} + \mathrm{i}2{g}_{2}\vert {A}_{2}\vert \mathrm{cos}({\omega }_{2}t + {\phi }_{2})\frac{{\alpha }_{0}}{{\epsilon}}.\end{equation} \tag{ 12 }$$

Solution of equation (11) provides the leading-order contribution to the cavity response function $\mathcal{C}={\sum }_{j}{{\epsilon}}^{j}{\mathcal{C}}_{j}$, with$\mathcal{C}=\alpha /{\text{e}}_{\mathrm{i}\mathrm{n}}$,

$$\begin{equation}{\mathcal{C}}_{0}=\sqrt{2{\kappa }_{\mathrm{i}\mathrm{n}}}\sum\limits _{b}{J}_{b}(-\beta ){\text{e}}^{\text{i}b{\omega }_{b}t}\enspace {C}_{b}(\xi ),\end{equation} \tag{ 13 }$$

where ξ = 2g₁|A₁|/ω₁, and

$$\begin{equation}{C}_{b}(\xi )=\sum\limits _{m,n}\frac{{J}_{m-n}(-\xi ){J}_{m}(-\xi )}{\mathrm{i}(m{\omega }_{1}+b{\omega }_{b})-\mathcal{W}}\enspace {\text{e}}^{\text{i}n({\omega }_{1}t+{\phi }_{1})},\end{equation} \tag{ 14 }$$

with $\mathcal{W}=\mathrm{i}{\Delta}-\kappa$. The first-order solution can be found to be

$$\begin{equation}{\epsilon}{\mathcal{C}}_{1}=\sqrt{2{\kappa }_{\mathrm{i}\mathrm{n}}}\sum\limits _{b}{J}_{b}(-\beta ){\text{e}}^{\text{i}b{\omega }_{b}t}\left[\delta {C}_{b}^{+}(\xi ){\text{e}}^{\text{i}({\omega }_{2}t+{\phi }_{2})}+\delta {C}_{b}^{-}(\xi ){\text{e}}^{-\text{i}({\omega }_{2}t+{\phi }_{2})}\right],\end{equation} \tag{ 15 }$$

where

$$\begin{equation}\delta {C}_{b}^{{\pm}}(\xi )=\frac{i{g}_{2}\vert {A}_{2}\vert }{\sqrt{2}}\sum\limits _{m,n}\frac{{J}_{m-n}(-\xi ){J}_{m}(-\xi )}{\mathrm{i}(m{\omega }_{1}+b{\omega }_{b})-\mathcal{W}}\frac{{\text{e}}^{\text{i}n({\omega }_{1}t+{\phi }_{1})}}{\mathrm{i}(m{\omega }_{1}+b{\omega }_{b}{\pm}{\omega }_{2})-\mathcal{W}}.\end{equation} \tag{ 16 }$$

Finally we observe that $\mathcal{C}$ is composed by a series of sidebands at frequencies nω₁ + bω_b ± ω₂, which are reproduced on the reflected output field determined by the input–output relation ${\text{e}}_{\mathrm{o}\mathrm{u}\mathrm{t}}=-{\text{e}}_{\mathrm{i}\mathrm{n}}\enspace \mathrm{exp}[-\mathrm{i}\beta \enspace \mathrm{sin}({\omega }_{b}t)]+\sqrt{2{\kappa }_{\mathrm{i}\mathrm{n}}}\alpha$.

The cavity reflection function, $\mathcal{R}={\text{e}}_{\mathrm{o}\mathrm{u}\mathrm{t}}/{\text{e}}_{\mathrm{i}\mathrm{n}}={\sum }_{j}{{\epsilon}}^{j}{\mathcal{R}}_{j}$, has a leading-order contribution

$$\begin{equation}{\mathcal{R}}_{0}=\sum\limits _{b}{J}_{b}(-\beta )\enspace {\text{e}}^{\text{i}b\enspace {\omega }_{b}t}\enspace {R}_{b}(\xi ).\end{equation} \tag{ 17 }$$

We observe that ${\mathcal{R}}_{0}$ is the superposition of the cavity reflection function R_b (ξ) at each sideband frequencies bω_b of the input field

$$\begin{equation}{R}_{b}(\xi )=-1+2{\kappa }_{\mathrm{i}\mathrm{n}}{C}_{b}(\xi ).\end{equation} \tag{ 18 }$$

The first-order contribution is

$$\begin{equation}{\mathcal{R}}_{1}=\sqrt{2{\kappa }_{\mathrm{i}\mathrm{n}}}{\mathcal{C}}_{1}.\end{equation} \tag{ 19 }$$

We consider now the modulation amplitude β, and the amplitude of the second mode, as small perturbations, that is β ≪ 1, and $\mathcal{R}\simeq {\mathcal{R}}_{0}+{\epsilon}{\mathcal{R}}_{1}$. In particular we focus on the contributions at DC, and at the five frequencies: ω₁, ω₂, ω_sm = 2ω₁ − ω₂, ω_b , and ω_sb = 2ω₁ − ω_b .

In this case we have

$$\begin{align}\mathcal{R}\simeq {\mathcal{R}}_{\mathrm{D}\mathrm{C}}& +{\mathcal{R}}_{+}({\omega }_{1}){\text{e}}^{\text{i}({\omega }_{1}t+{\phi }_{1})}+{\mathcal{R}}_{-}({\omega }_{1}){\text{e}}^{-\text{i}({\omega }_{1}t+{\phi }_{1})}\\ & +{\mathcal{R}}_{+}({\omega }_{2}){\text{e}}^{\text{i}({\omega }_{2}t+{\phi }_{2})}+{\mathcal{R}}_{-}({\omega }_{2}){\text{e}}^{-\text{i}({\omega }_{2}t+{\phi }_{2})}\\ & +{\mathcal{R}}_{+}({\omega }_{\mathrm{s}\mathrm{m}}){\text{e}}^{\text{i}({\omega }_{\mathrm{s}\mathrm{m}}t+{\phi }_{\mathrm{s}\mathrm{m}})}+{\mathcal{R}}_{-}({\omega }_{\mathrm{s}\mathrm{m}}){\text{e}}^{-\text{i}({\omega }_{\mathrm{s}\mathrm{m}}t+{\phi }_{\mathrm{s}\mathrm{m}})}\\ & +{\mathcal{R}}_{+}({\omega }_{b}){\text{e}}^{\text{i}{\omega }_{b}t}+{\mathcal{R}}_{-}({\omega }_{\mathrm{s}\mathrm{b}}){\text{e}}^{-\text{i}{\omega }_{b}t}\\ & +{\mathcal{R}}_{+}({\omega }_{\mathrm{s}\mathrm{b}}){\text{e}}^{\text{i}({\omega }_{\mathrm{s}\mathrm{b}}t+{\phi }_{\mathrm{s}\mathrm{b}})}+{\mathcal{R}}_{-}({\omega }_{\mathrm{s}\mathrm{b}}){\text{e}}^{-\text{i}({\omega }_{\mathrm{s}\mathrm{b}}t+{\phi }_{\mathrm{s}\mathrm{b}})},\end{align} \tag{ 20 }$$

The contribution at DC is provided by equation (17) for b = 0, and the time-independent term (n = 0) of equation (14)

$$\begin{equation}{\mathcal{R}}_{\mathrm{D}\mathrm{C}}={J}_{0}(-\beta )\left[-1+2{\kappa }_{\mathrm{i}\mathrm{n}}\sum\limits _{m}\frac{{J}_{m}^{2}(-\xi )}{\mathrm{i}m{\omega }_{1}-\mathcal{W}}\right];\end{equation} \tag{ 21 }$$

at ω₁ is given by equation (17) with b = 0, and in equation (14) the term with n = ±1

$$\begin{equation}{\mathcal{R}}_{{\pm}}({\omega }_{1})={J}_{0}(-\beta )2{\kappa }_{\mathrm{i}\mathrm{n}}\sum\limits _{m}\frac{{J}_{m}(-\xi ){J}_{m\mp 1}(-\xi )}{\mathrm{i}m{\omega }_{1}-\mathcal{W}};\end{equation} \tag{ 22 }$$

at ω₂ is provided by equation (15) with b = 0, and in equation (16) the term with n = 0

$$\begin{equation}{\mathcal{R}}_{{\pm}}({\omega }_{2})={J}_{0}(-\beta )2{\kappa }_{\mathrm{i}\mathrm{n}}\frac{\mathrm{i}{g}_{2}\vert {A}_{2}\vert }{\sqrt{2}}\sum\limits _{m}\frac{{J}_{m}^{2}(-\xi )}{\mathrm{i}m{\omega }_{1}-\mathcal{W}}\frac{1}{\mathrm{i}(m{\omega }_{1}{\pm}{\omega }_{2})-\mathcal{W}};\end{equation} \tag{ 23 }$$

for ω_sm = 2ω₁ − ω₂, b = 0 in equation (15), and n = ∓2 in equation (16)

$$\begin{equation}{\mathcal{R}}_{{\pm}}({\omega }_{\mathrm{s}\mathrm{m}}) = {J}_{0}(-\beta )2{\kappa }_{\mathrm{i}\mathrm{n}}\frac{\mathrm{i}{g}_{2}\vert {A}_{2}\vert }{\sqrt{2}}\sum\limits _{m}\frac{{J}_{m}(-\xi ){J}_{m\mp 2}(-\xi )}{\mathrm{i}m{\omega }_{1}-\mathcal{W}}\frac{1}{\mathrm{i}(m{\omega }_{1}\mp {\omega }_{2})-\mathcal{W}};\end{equation} \tag{ 24 }$$

at ω_b is provided by equation (17) with b = ±1, and in equation (14) the term with n = 0

$$\begin{equation}{\mathcal{R}}_{{\pm}}({\omega }_{b})={J}_{{\pm}1}(-\beta )\left[-1+2{\kappa }_{\mathrm{i}\mathrm{n}}\sum\limits _{m}\frac{{J}_{m}^{2}(-\xi )}{\mathrm{i}(m{\omega }_{1}{\pm}{\omega }_{b})-\mathcal{W}}\right];\end{equation} \tag{ 25 }$$

and the sideband at ω_sb = 2ω₁ − ω_b for b = ±1, and n = 0

$$\begin{equation}{\mathcal{R}}_{{\pm}}({\omega }_{\mathrm{s}\mathrm{b}})={J}_{\mp 1}(-\beta )2{\kappa }_{\mathrm{i}\mathrm{n}}\sum\limits _{m}\frac{{J}_{m}(-\xi ){J}_{m\mp 2}(-\xi )}{\mathrm{i}(m{\omega }_{1}\mp {\omega }_{b})-\mathcal{W}}.\end{equation} \tag{ 26 }$$

The homodyne technique is implemented by mixing on a beam-splitter the reflected field e_out with an intense local oscillator ${\text{e}}_{\mathrm{l}\mathrm{o}}\enspace {\text{e}}^{\text{i}{\phi }_{\mathrm{l}\mathrm{o}}}$, and detecting the fields at the output of the beam-splitter, i.e., ${\text{e}}_{j}=[{\text{e}}_{\mathrm{l}\mathrm{o}}\enspace {\text{e}}^{\text{i}{\phi }_{lo}}+{(-1)}^{j}{\text{e}}_{\mathrm{o}\mathrm{u}\mathrm{t}}]/\sqrt{2}$. The output power are ${\mathcal{P}}_{{\pm}}=\vert {\text{e}}_{\mathrm{l}\mathrm{o}}\enspace {\text{e}}^{\text{i}{\phi }_{\mathrm{l}\mathrm{o}}}{\pm}{\text{e}}_{\mathrm{o}\mathrm{u}\mathrm{t}}{\vert }^{2}/2$, and the differential current

$$\begin{equation}{I}_{\mathrm{D}}=S[{\mathcal{P}}_{+}-{\mathcal{P}}_{-}]=2S\sqrt{{\mathcal{P}}_{\mathrm{l}\mathrm{o}}{\mathcal{P}}_{\mathrm{i}\mathrm{n}}}\enspace \text{Re}\left[\mathcal{R}\enspace {\text{e}}^{-\text{i}{\phi }_{\mathrm{l}\mathrm{o}}}\right],\end{equation} \tag{ 27 }$$

where S is the sensitivity of the photodiodes, and the phase ϕ_lo represents the controllable phase difference between the local oscillators and the reflected field. The local oscillator phase is locked to have zero DC signal, which turns out to be ${\phi }_{\mathrm{l}\mathrm{o}}=\mathrm{arctan}\left[\text{Re}({\mathcal{R}}_{\mathrm{D}\mathrm{C}})/\text{Im}({\mathcal{R}}_{\mathrm{D}\mathrm{C}})\right]$. This phase is plotted in figure 3 for our experimental parameters, showing that it does not deviate from the optimal value of π/2. The error signal, assuming optimal detection for the local oscillator phase ϕ_lo ≃ π/2, is

$$\begin{equation}{V}_{\mathrm{H}}(t)={g}_{\mathrm{T}}S2\sqrt{{\mathcal{P}}_{\mathrm{l}\mathrm{o}}{\mathcal{P}}_{\mathrm{i}\mathrm{n}}}\enspace \text{Im}\left[\mathcal{R}\right],\end{equation} \tag{ 28 }$$

where g_T is the transimpedance gain. This voltage contains the signature of any modulation frequency of the reflected field, provided that it falls within the bandwidth of the electronic system. The single-sided power spectral density ${\mathcal{S}}_{\mathrm{W}}(\omega )=\int \mathrm{d}\tau \enspace {\text{e}}^{\text{i}\omega \tau }{\langle {V}_{\mathrm{H}}(t+\tau ){V}_{\mathrm{H}}(t)\rangle }_{t}/{R}_{0}$ on a termination resistor R₀, provides a normalized amplitude spectral noise for each well separated frequency, that is

$$\begin{equation}{\mathcal{V}}_{\mathrm{H}}({\Omega})=\sqrt{\frac{{\mathcal{S}}_{\mathrm{W}}({\Omega})}{{\mathcal{S}}_{0}}}=\frac{1}{2}\left\vert {\mathcal{R}}_{+}({\Omega})-{\mathcal{R}}_{-}^{{\ast}}({\Omega})\right\vert ,\end{equation} \tag{ 29 }$$

with ${\mathcal{S}}_{0}={({g}_{\mathrm{T}}2S)}^{2}{\mathcal{P}}_{\mathrm{l}\mathrm{o}}{\mathcal{P}}_{\mathrm{i}\mathrm{n}}/{R}_{0}$. In figure 3 are reported the amplitude spectral noise at the frequencies of interest for our experiment.

Zoom In Zoom Out Reset image size

We notice that for probe detuning Δ = 0, and when the cavity field is weakly modulated at frequency ω₁, that is ξ ≪ 1, the output signal is linear in the displacement

$$\begin{equation}{\mathcal{R}}_{+}({\omega }_{1})=-{[{\mathcal{R}}_{-}({\omega }_{1})]}^{{\ast}}\sim -\frac{\xi }{2}\frac{2{\kappa }_{\mathrm{i}\mathrm{n}}}{\kappa }\frac{\mathrm{i}{\omega }_{1}}{\kappa +\mathrm{i}{\omega }_{1}},\hspace{42.67912pt}\end{equation} \tag{ 30 }$$

and the amplitude spectral noise is

$$\begin{equation}\sqrt{{\mathcal{S}}_{\mathrm{W}}({\omega }_{1})}\left[\mathrm{W}/\sqrt{\mathrm{H}\mathrm{z}}\right]=\frac{{g}_{\mathrm{T}}S}{\sqrt{{R}_{0}}}\cdot 2\sqrt{{\mathcal{P}}_{\mathrm{l}\mathrm{o}}{\mathcal{P}}_{\mathrm{i}\mathrm{n}}}\cdot \frac{\mathcal{F}}{{\lambda }_{0}}\frac{\delta \tilde {x}({\omega }_{1})\left[\mathrm{m}/\sqrt{\mathrm{H}\mathrm{z}}\right]}{\sqrt{1+{\omega }_{1}^{2}/{\kappa }^{2}}} \cdot \eta ,\end{equation} \tag{ 31 }$$

where $\mathcal{F}=\mathrm{F}\mathrm{S}\mathrm{R}/2\kappa$, $\xi =\delta \tilde {x}({\omega }_{1}){g}_{1}/{\omega }_{1}{x}_{\mathrm{z}\mathrm{p}\mathrm{f}}$, and η = [2κ_in/κ] × [g₁ λ₀/2 FSR x_zpf] ∼ 0.25 × 0.24 ∼ 0.06 for our setup. The average power falling on each photodiode is approximately ${\mathcal{P}}_{\mathrm{l}\mathrm{o}}/2$. The shot noise in the differential signal has a flat spectrum with spectral density ${S}_{\mathcal{P}\mathcal{P}}^{\mathrm{s}\mathrm{n}}\sim 2{\times}2\hslash {\omega }_{L}{\mathcal{P}}_{\mathrm{l}\mathrm{o}}/2$, which sets a limit to the sensitivity of the detection [67, 68]

$$\begin{equation}\delta \tilde {x}({\omega }_{1})\left[\mathrm{m}/\sqrt{\mathrm{H}\mathrm{z}}\right]=\frac{1}{\sqrt{2{\mathcal{P}}_{\mathrm{i}\mathrm{n}}/\hslash {\omega }_{L}}}\frac{{\lambda }_{0}}{\mathcal{F}}\frac{1}{\eta }\sqrt{1+\frac{{\omega }_{1}^{2}}{{\kappa }^{2}}}\sim 4{\times}1{0}^{-16}\enspace \mathrm{m}/\sqrt{\mathrm{H}\mathrm{z}},\end{equation} \tag{ 32 }$$

and reproduces the shot-noise limited displacement detection filtered by the cavity response, for which the shot-noise limited sensitivity is not flat in the spectrum (Mizuno's sensitivity theorem [69]). The cavity response length in this linear detection regime might be estimated as ${\lambda }_{0}/2\mathcal{F}\simeq$ 43 pm, which corresponds to ${\xi }_{\mathrm{c}\mathrm{a}\mathrm{v}}=2{g}_{1}/{\omega }_{1}\cdot {\lambda }_{0}/2\mathcal{F}{x}_{\mathrm{z}\mathrm{p}\mathrm{f}}\simeq 0.284$.

On the contrary, when the excited mechanical mode reaches a limit cycle with large amplitude, that is for ξ approaching unity, the reflected signal, due to the non-linear response of the cavity, presents a reduction of the signal of the unexcited mode and the appearance of a sideband at ω_sm and ω_sb (see figure 3). In this case any attempt to determine the mechanical displacement from the measured phase of the output field requires careful attention, also because the calibration tone, which is implemented by modulation of the input field, is essentially unaffected (see magenta line in figure 3). In general, two main ratios of the normalized spectral amplitudes ${\mathcal{V}}_{\mathrm{H}}({\Omega})$, which are independent from β, κ_in, and from the coupling term with the second mechanical mode, g₂|A₂|, describe the non-linear dynamics of the optomechanical system in terms of the normalized mechanical amplitude of the first mechanical oscillator ξ: (i) the ratio derived by equation (22) and the product of β and equation (25); (ii) the ratio derived by equations (23) and (24)

$$\begin{equation}{\mathcal{T}}_{1}=\frac{{\mathcal{V}}_{\mathrm{H}}({\omega }_{1})\cdot \beta }{{\mathcal{V}}_{\mathrm{H}}({\omega }_{b})},\hspace{25.0pt}{\mathcal{T}}_{\mathrm{s}\mathrm{b}}=\frac{{\mathcal{V}}_{\mathrm{H}}({\omega }_{\mathrm{s}\mathrm{b}})}{{\mathcal{V}}_{\mathrm{H}}({\omega }_{2})}.\end{equation} \tag{ 33 }$$

In panel (a) of figure 4 are reported these ratios that allow us to estimate ξ for our experimental realisations. In panel (b) the ratios ${\mathcal{T}}_{1}(\xi ){\xi }^{-1}$, and ${\mathcal{T}}_{2}(\xi )={\mathcal{V}}_{\mathrm{H}}({\omega }_{2})\cdot \beta /{\mathcal{V}}_{\mathrm{H}}({\omega }_{b})$ (note that ${\mathcal{T}}_{2}$ does depend on g₂|A₂|), normalized to their maximum values, which are reached for ξ → 0, are shown as ${\mathcal{N}}_{1}$, and ${\mathcal{N}}_{2}$, respectively. By definition, they represent the correction factors that relate the displacement amplitudes detected via the reflected probe spectrum, to the effective displacements amplitudes. ${\mathcal{N}}_{1}$, and${\mathcal{N}}_{2}$ are equal to 1, corresponding to the usual linear detection regime, for ξ ≪ 1. In the present case instead, for the results reported in figures 2 and 7 we estimate ${\xi }_{1}^{\mathrm{s}\mathrm{t}}\simeq 1.05$, ${\mathcal{N}}_{1}\simeq 0.70$ and ${\mathcal{N}}_{2}\simeq 0.42$. Such analysis allows us to deduce the effective limit cycle displacement amplitude of the first oscillator to be ${q}_{1}^{\mathrm{s}\mathrm{t}}={\xi }_{1}^{\mathrm{s}\mathrm{t}}{\omega }_{1}{x}_{\mathrm{z}\mathrm{p}\mathrm{f}}/{g}_{1}\simeq$ 263 pm, and the observed limit cycle displacement amplitude ${q}_{1}^{\mathrm{o}\mathrm{b}}={q}_{1}^{\mathrm{s}\mathrm{t}}{\mathcal{N}}_{1}\simeq$ 183 pm. This analysis does not allow us, for the moment, to draw any quantitative conclusion for what concerns the second oscillator. We will be able to do that in the following section, when we will analyse the experimental time traces. We anticipate here that also the optical detection of the unexcited oscillator is affected by the nonlinear response of the cavity, and the complementary analysis in time domain will allow us to determine the proper correction to the calibration readout ${\mathcal{N}}_{2}$.

Zoom In Zoom Out Reset image size

In contrast to the systems implemented in [55, 57], our high finesse optical system, used for detection, does not allow us to reveal the effective motion of each membrane independently. However, by using the above analysis of the reflected spectrum, and the numerical simulations, it is possible to unambiguously infer from the homodyne detection of the probe beam shown in figure 2 that the dynamics of the two membranes is characterized in this parameter region by a pre-synchronisation regime. In fact the numerical integration of equations (1) and (2) with the parameters reported in table 1, which refers to the results of figure 2, also shows that when the output probe beam exhibits sidebands around the frequency of the fundamental mode of the first oscillator [see figure 2(a))], the second oscillator has a nonzero amplitude of oscillation at the frequency of the first oscillator, that is, it starts to synchronize with the first mode. This is emphasised by the magenta box in panel (b) of figure 5, which shows the numerical simulations of the spectra of the fundamental mechanical mode of the two oscillators for the experimental parameters of table 1. A more quantitative description of such pre-synchronization process of the second resonator with the excited first one, acting as the 'master' oscillator, is provided by the synchronization measure [45]

$$\begin{equation}{P}_{{\theta }_{-}}(t)=\frac{1}{{\Delta}t}{\int }_{t-{\Delta}t/2}^{t+{\Delta}t/2}\mathrm{cos}\left[{\theta }_{1}(t)-{\theta }_{2}(t)\right]\mathrm{d}t,\end{equation} \tag{ 34 }$$

where θ_j (t) = arg[β_j (t)], reported in figure 5(c), which shows an increase of the phase anti-correlation between the two oscillators. The effect is small due to the weak driving regime, but it is nonetheless unambiguously present.

Zoom In Zoom Out Reset image size

We are also able to provide a consistent analytical description of this pre-synchronization dynamics of the two membrane modes starting from the slowly varying amplitude equations (4) and (5). To study the regime when the first oscillator reaches a limit cycle while the second is not excited, it is convenient to take ω_ref = ω₁ as a reference in equations (4) and (5), so that Δω₁ = 0 and Δω₂ = Δω. One can make quantitative predictions on such a regime assuming that $\vert {A}_{1}\vert \gg \vert {A}_{2}\vert ,\sqrt{2{\bar{n}}_{1}}$. Moreover, in our experiment the optical noise is negligible and we will not consider the terms associated with ${\eta }_{i}^{\mathrm{o}\mathrm{p}\mathrm{t}}(t)$. With the above approximations, one can neglect both thermal noise and A₂ contributions in equation (4), which becomes

$$\begin{equation}{\dot {A}}_{1}(t) = \left[-{\gamma }_{1}+\mathrm{i}{d}_{1}(\vert {A}_{1}\vert )\right]{A}_{1}(t) = -\left[{\gamma }_{1}^{\mathrm{e}\mathrm{ff}}(\vert {A}_{1}\vert )-\mathrm{i}{\Delta}{\omega }_{1}^{\mathrm{e}\mathrm{ff}}(\vert {A}_{1}\vert )\right]{A}_{1}(t),\end{equation} \tag{ 35 }$$

where we have made explicit the dependence of d₁ on |A₁|, and defined ${\Delta}{\omega }_{1}^{\mathrm{e}\mathrm{ff}}(\vert {A}_{1}\vert )=-\text{Re}[{d}_{1}(\vert {A}_{1}\vert )]$, and ${\gamma }_{1}^{\mathrm{e}\mathrm{ff}}(\vert {A}_{1}\vert )={\gamma }_{1}+\text{Im}[{d}_{1}(\vert {A}_{1}\vert )]$. The effective mechanical damping can be cast as

$$\begin{equation}{\gamma }_{1}^{\mathrm{e}\mathrm{ff}}(\vert {A}_{1}\vert )={\gamma }_{1}\left[1+\frac{{g}_{1}}{{\gamma }_{1}\vert {A}_{1}\vert }\enspace \text{Im}\left[{E}_{1}^{2}{{\Sigma}}_{1}+{E}_{2}^{2}{{\Sigma}}_{2}\right]\right],\end{equation} \tag{ 36 }$$

where we have used the fact that in the considered regime $\vert {A}_{i}^{b}\vert \simeq \vert {A}_{1}\vert ({g}_{i1}/{g}_{i}^{b})$, ξ_j ≃ ξ₁ = 2g₁|A₁|/ω₁, assumed g_i1 ≃ g₁, and

$$\begin{equation}{{\Sigma}}_{j}\equiv {\Sigma}({\xi }_{1},{\kappa }_{i},{{\Delta}}_{i})=\sum\limits _{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}^{{\ast}}]}.\end{equation} \tag{ 37 }$$

We note that such approximation implies a regime where the second mode is still dominated by a thermal dynamics, i.e., pre-synchronized regime; on the contrary, if also the second mode would have reached a limit cycle, synchronized with the first, the amplitude g₂ A₂ would have not been negligible anymore with respect to g₁ A₁, and the dynamics would have been governed by the more general equations (4) and (5) [55, 57]. Equation (35) can be solved by rewriting it in terms of modulus and phase, ${A}_{1}={I}_{1}\enspace {\text{e}}^{\text{i}{\phi }_{1}}$,

$$\begin{equation}{\dot {I}}_{1}(t)=-{\gamma }_{1}^{\mathrm{e}\mathrm{ff}}({I}_{1})\enspace {I}_{1}(t),\end{equation} \tag{ 38 }$$

$$\begin{equation}{\dot {\phi }}_{1}(t)={\Delta}{\omega }_{1}^{\mathrm{e}\mathrm{ff}}({I}_{1}).\end{equation} \tag{ 39 }$$

After a transient these equations yield a steady state with a constant radius of the limit cycle of the first oscillator, ${I}_{1}^{\mathrm{s}\mathrm{t}}$, corresponding in our case of not too strong driving, to the smallest positive root of the implicit equation ${\gamma }_{1}^{\mathrm{e}\mathrm{ff}}({I}_{1}^{\mathrm{s}\mathrm{t}})=0$, which can be cast as

$$\begin{equation}a\vert {\xi }_{1}^{\mathrm{s}\mathrm{t}}\vert =-\text{Im}\left[{E}_{1}^{2}{\Sigma}({\xi }_{1}^{\mathrm{s}\mathrm{t}},{\kappa }_{1},{{\Delta}}_{1})+{E}_{2}^{2}{\Sigma}({\xi }_{1}^{\mathrm{s}\mathrm{t}},{\kappa }_{2},{{\Delta}}_{2})\right],\end{equation} \tag{ 40 }$$

with ${\xi }_{1}^{\mathrm{s}\mathrm{t}}=2{g}_{1}{I}_{1}^{\mathrm{s}\mathrm{t}}/{\omega }_{1}$, and $a={\omega }_{1}{\gamma }_{1}/2{g}_{1}^{2}$. As a consequence, at long times,${\phi }_{1}^{\mathrm{s}\mathrm{t}}(t)\simeq t{\Delta}{\omega }_{1}^{\mathrm{s}\mathrm{t}}$ with ${\Delta}{\omega }_{1}^{\mathrm{s}\mathrm{t}}={\Delta}{\omega }_{1}^{\mathrm{e}\mathrm{ff}}({I}_{1}^{\mathrm{s}\mathrm{t}})$ so that ${A}_{1}(t)\simeq {A}_{1}^{\mathrm{s}\mathrm{t}}(t)={I}_{1}^{\mathrm{s}\mathrm{t}}\enspace \mathrm{exp}[\mathrm{i}t{\Delta}{\omega }_{1}^{\mathrm{s}\mathrm{t}}]$.

In figure 6 we show the left and right side of equation (40) for the experimental parameters of table 1, which provides the optomechanical parameters for the results reported in figures 2 and 7. We infer from the intersection point, which corresponds to find the smallest root of${\gamma }_{1}^{\mathrm{e}\mathrm{ff}}({I}_{1}^{\mathrm{s}\mathrm{t}})=0$, a value ${\xi }_{1}^{\mathrm{s}\mathrm{t}}=1.054$, a steady-state displacement amplitude ${q}_{1}^{\mathrm{s}\mathrm{t}}=2\vert {A}_{1}\vert {x}_{\mathrm{z}\mathrm{p}\mathrm{f}}=$ 263.0 pm, and ${\Delta}{\omega }_{1}^{\mathrm{e}\mathrm{ff}}=-2\pi \cdot$0.04 Hz, confirming the value of ${q}_{1}^{\mathrm{s}\mathrm{t}}$ deduced in figure 4. We emphasise once more that, as shown in figure 4, for ξ ≪ 1 the spectral amplitude of the sideband of the output field is linear with ξ and provides a direct measurement of the mechanical position coordinate q₁; on the contrary for ξ ≥ 1 linearity is no more valid and a proper correction factor should be considered. In our case, as obtained in figure 4, the theoretical correction factor is ${\mathcal{N}}_{1}\simeq 0.70$, corresponding to an expected observable stationary limit cycle amplitude of ${q}_{1}^{\mathrm{o}\mathrm{b}}\simeq$ 183 pm fully consistent with the analysis described in the previous section. It is worth noting that, due to the oscillating behaviour of the Bessel functions, equation (40) may have more than one solution at sufficiently large power (see below the oblique black dashed line in figure 6), corresponding to the multistability phenomenon analysed in reference [32] and experimentally verified in reference [34].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now consider the dynamics of the second oscillator inserting the steady-state solution for A₁(t) into equation (5), which becomes

$$\begin{equation}{\dot {A}}_{2}(t)=\left[-{\gamma }_{2}-\mathrm{i}{\Delta}\omega +\mathrm{i}{d}_{2}({I}_{1}^{\mathrm{s}\mathrm{t}})\right]{A}_{2}(t)+\mathrm{i}{d}_{12}({I}_{1}^{\mathrm{s}\mathrm{t}}){A}_{1}^{\mathrm{s}\mathrm{t}}(t)+\sqrt{2{\gamma }_{2}}{\beta }_{2}^{\mathrm{i}\mathrm{n}}(t).\end{equation} \tag{ 41 }$$

The stationary solution can be obtained via Fourier transform and it can be written as

$$\begin{equation}{A}_{2}(t)=\frac{\mathrm{i}{d}_{12}({I}_{1}^{\mathrm{s}\mathrm{t}})}{{\gamma }_{2}^{\mathrm{e}\mathrm{ff}}+\mathrm{i}{\Delta}{\bar{\omega }}_{2}^{\mathrm{e}\mathrm{ff}}}{A}_{1}^{\mathrm{s}\mathrm{t}}(t)+\sqrt{2{\gamma }_{2}}{\int }_{0}^{t}\mathrm{d}s\enspace {\text{e}}^{-\left({\gamma }_{2}^{\mathrm{e}\mathrm{f}\mathrm{f}}+\mathrm{i}{\Delta}{\omega }_{2}^{\mathrm{e}\mathrm{f}\mathrm{f}}\right)s}{\beta }_{2}^{\mathrm{i}\mathrm{n}}(t-s),\end{equation} \tag{ 42 }$$

where ${\gamma }_{2}^{\mathrm{e}\mathrm{ff}}={\gamma }_{2}+\text{Im}[{d}_{2}({I}_{1}^{\mathrm{s}\mathrm{t}})]$ is positive, i.e., the second resonator is still damped despite the pump driving and it is not driven into a limit cycle,${\Delta}{\bar{\omega }}_{2}^{\mathrm{e}\mathrm{ff}}={\Delta}\omega +\text{Re}[{d}_{1}({I}_{1}^{\mathrm{s}\mathrm{t}})-{d}_{2}({I}_{1}^{\mathrm{s}\mathrm{t}})]$, and ${\Delta}{\omega }_{2}^{\mathrm{e}\mathrm{ff}}={\Delta}\omega -\text{Re}[{d}_{2}({I}_{1}^{\mathrm{s}\mathrm{t}})]$. Therefore the first term on the right-hand side of equation (42) is the synchronized component oscillating at the same frequency of the first master oscillator [see figure 5(b))], while the second term is the thermal noise component at its natural frequency. This equation describes how the second resonator is driven towards synchronization with the first resonator, and full synchronization and phase locking is achieved when the thermal contribution is negligible, i.e., when $\vert {d}_{12}({I}_{1}^{\mathrm{s}\mathrm{t}}){\vert }^{2}{I}_{1}^{\mathrm{s}\mathrm{t},\mathrm{2}}\gg {\gamma }_{2}{\gamma }_{2}^{\prime }{\bar{n}}_{2}$ (where we have exploited the fact that${\Delta}{\bar{\omega }}_{2}^{\mathrm{e}\mathrm{ff}}\simeq {\Delta}{\omega }_{2}^{\mathrm{e}\mathrm{ff}}\simeq {\Delta}\omega$). This transition to synchronization is consistent with the theoretical analysis made in references [5, 45], which, in the regime of not large driving power studied here, predicts an onset of synchronization with very different limit cycle amplitudes, even in the presence of thermal noise. In the four-dimensional phase space of the mechanical oscillators it manifests itself via a Neimark–Sacker bifurcation corresponding to the birth of a stable torus around the existing limit cycle [36].

We corroborate such analysis by considering the experimental time traces shown in figure 7. The voltage V_xj and V_yj [panels (a) and (b)], are the slowly varying quadratures of the voltage signal, integrated over a bandwidth of 70 Hz–150 Hz around the mechanical frequencies ω₁, and ω₂, as a function of time. In panels (c) and (d) are reported the associated phase-space distributions. Panels (e) and (f) show the phase-space distributions of the calibration tone before and after the pump is turned on, confirming that the calibration tone is not appreciably affected.

By means of the calibration tone [66], firstly, we determine the displacement amplitudes of the two oscillators, shown in figures 8(b) and (c), which, before the pump beam is turned on (t < 4 s), show higher values than the thermal ones. This is ascribed to a slightly blue-detuning of the probe beam. We evaluate such detuning observing that, for the second mode, green curve in panel (b), the calibrated measured position standard deviation ${\Delta}{q}_{2}^{{\Delta}}\simeq$ 3.50 pm, while the estimated thermal position standard deviation is ${\Delta}{q}_{2}^{\mathrm{t}\mathrm{h}}\simeq$ 3.32 pm, so that

$$\begin{equation}\frac{{\Delta}{q}_{2}^{{\Delta}}}{{\Delta}{q}_{2}^{\mathrm{t}\mathrm{h}}}=\frac{1}{\sqrt{1+C({\Delta})}}\simeq 1.054,\end{equation} \tag{ 43 }$$

where, for almost resonant field [1],

$$\begin{equation}C({\Delta})\sim -\frac{2{g}_{j}^{2}{E}^{2}}{\gamma \kappa }\frac{4{\omega }_{j}}{{({\kappa }^{2}+{\omega }_{j}^{2})}^{2}}{\Delta}.\end{equation} \tag{ 44 }$$

This fact allows us to estimate a small blue detuning of the probe Δ₂ ≃ 2π × 3.9 kHz, which is the value provided in table 1.

Zoom In Zoom Out Reset image size

Finally we analyse the measured mechanical amplitudes after the pump beam has been turned on (t > 4 s): the amplitude of the first oscillator increases, while the measured amplitude of the second one reduces below the thermal value. For the first oscillator, the observed steady-state limit cycle displacement amplitude is ${q}_{1}^{\mathrm{o}\mathrm{b}}\simeq$ 184 pm, orange curve in figure 8(c). Such value agrees very well with the expected one, shown as dark-orange filled circle in panel (a) of figure 8. Panel (a) represents the effective steady-state mechanical amplitudes obtained as solution of equation (40) (blue curves), and the observed one (orange curves), that is, reduced by the correction factor reported in figure 4, and calculated in section 4. For a given set of parameters we observe that both the effective and observed steady-state mechanical amplitudes reach a maximum as a function of g, and then decrease. This result is confirmed by integration of equation (4) [and equation (5)], including the noise contributions, and finding the steady-state as mean amplitude after the oscillator has reached the limit cycle, values which are reported as diamond symbols in figure 8(a). The effective steady-state displacement amplitude of the first oscillator, ${q}_{1}^{\mathrm{s}\mathrm{t}}=2\vert {A}_{1}\vert {x}_{\mathrm{z}\mathrm{p}\mathrm{f}}\simeq$ 262 pm, is also confirmed by the ten blue trajectories simulated with the parameters of table 1, and reported in figure 8(c). We note that even the slope of the trajectories follows with accuracy the measured one, implying that our approach in terms of slowly-varying complex amplitudes of the two oscillators, is effective, and able to grasp all the features of the non-linear dynamics.

Lastly, we observe that even the dynamics of the second mode is very well described by our model. In fact, it is evident that after the pump is turned on, the behaviour of the observed q₂ [green curve in figure 8(b))] follows the dynamics of the effective mechanical displacement [blue trajectories in figure 8(b))] only until q₁ reaches the limit cycle (after 7 s), and since then the observed displacement differs from the effective one. An estimation of the effective steady-state amplitude of the second contribution in equation (42) is provided by ${\Delta}{q}_{2}^{\mathrm{t}\mathrm{h}}\sqrt{{\gamma }_{2}/{\gamma }_{2}^{\mathrm{e}\mathrm{ff}}}$, with

$$\begin{equation}\frac{{\gamma }_{2}^{\mathrm{e}\mathrm{ff}}}{{\gamma }_{2}}=1+\frac{\text{Im}[{d}_{2}({I}_{1}^{\mathrm{s}\mathrm{t}})]}{{\gamma }_{2}}=1+\frac{{g}_{2}^{2}{E}_{1}^{2}}{{\gamma }_{2}{I}_{1}{g}_{1}}\enspace \text{Im}\left[{{\Sigma}}_{1}+\eta \enspace {{\Sigma}}_{2}\right]=1-\frac{{\gamma }_{1}{g}_{2}^{2}}{{\gamma }_{2}{g}_{1}^{2}}.\end{equation} \tag{ 45 }$$

We note that, in our case, the effective amplitude is larger than the thermal one by a factor $\sqrt{{\gamma }_{2}/{\gamma }_{2}^{\mathrm{e}\mathrm{ff}}}\simeq 1.38$, that is, there is a small effective driving, although not enough for the appearance of a limit cycle. Also, from the correction factor ${\mathcal{N}}_{2}\simeq 0.42$, we estimate an observed displacement of ${q}_{2}^{\mathrm{o}\mathrm{b}}={\Delta}{q}_{2}^{{\Delta}}\cdot \sqrt{{\gamma }_{2}/{\gamma }_{2}^{\mathrm{e}\mathrm{ff}}}\cdot {\mathcal{N}}_{2}\simeq$ 3.50 pm × 1.38 × 0.42 ≃ 2.0 pm, in great agreement with the measured value of 2.1 pm. In conclusion, we observe that even the small effective amplitude displacement of the second oscillator is strongly affected by the non-linear cavity dynamics, exhibiting a fictitious cooling effect, which is instead only a manifestation of detection in this nonlinear regime. This is a somehow unexpected effect of the nonlinear regime in which our system operates. As soon as the amplitude of one of the resonators yields a frequency modulation larger than the cavity linewidth, all the optically detected motional amplitudes are nonlinearly modified and appropriate calibration factors ${\mathcal{N}}_{j}$ must be considered. This occurs also to the unexcited resonator whose motional amplitude corresponds to a cavity frequency modulation much smaller than the cavity linewidth.

We have presented a detailed experimental analysis of the dynamics of the multimode optomechanical system introduced in reference [14], formed by a sandwich of two membrane mechanical resonators placed within a high-finesse cavity, and interacting with a pump and a probe cavity mode. We have focussed onto the non-linear regime where a blue-detuned pump drives one of the two oscillators into a self-sustained limit cycle. In the weak driving regime studied here, the system is in a pre-synchronized situation where the unexcited oscillator has a small, synchronized component at the frequency of the excited (master) oscillator, which is however dominated by the fluctuating thermal noise component. We find perfect agreement between the experimental results, the numerical simulations, and an analytical approach based on slowly-varying amplitude equations. This analytical study allows to derive a full and consistent description of the displacement detection by the probe beam in this non-linear regime, enabling the faithful detection of membrane displacements well above the usual sensing limit corresponding to the cavity linewidth. In this non-linear detection regime, both large and small amplitude resonator motion are transduced in a nontrivial way by the non-linear response of the optical probe beam.

The data that support the findings of this study are available upon reasonable request from the authors.

We acknowledge the support of the European Union Horizon 2020 Programme for Research and Innovation through the Project No. 732894 (FET Proactive HOT) and the Project QuaSeRT funded by the QuantERA ERA-NET Cofund in Quantum Technologies. P Piergentili acknowledges support from the European Union's Horizon 2020 Programme for Research and Innovation under Grant Agreement No. 722923 (Marie Curie ETN—OMT).