Two-membrane cavity optomechanics: non-linear dynamics

We study the non-linear dynamics of a multimode optomechanical system constituted of a driven high-finesse Fabry-P\'erot cavity containing two vibrating dielectric membranes. The analytical study allows to derive a full and consistent description of the displacement detection by a probe beam in the non-linear regime, enabling the faithful detection of membrane displacements well above the usual sensing limit corresponding to the cavity linewidth. In the weak driving regime where the system is in a pre-synchronized situation, the unexcited oscillator has a small, synchronized component at the frequency of the excited one; both large and small amplitude resonator motions are transduced in a nontrivial way by the non-linear response of the optical probe beam. We find perfect agreement between the experimental results, the numerical simulations, and an analytical approach based on slowly-varying amplitude equations.


Introduction
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,3,4,5,6,7,8,9,10,11,12,13,14,15,16,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 tunneling of photons and phonons between neighboring sites [18,19,20,21,22,23], (e.g. optomechanical crystals in one and two dimensions [24]).
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,31,32,33,34,35,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,39,40,41,42,43,44,45], and experimentally [46,47,48,49,50,51,52,53,54,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,59,60]. This radiation-pressure-induced nonlinear behavior 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 Ref. [14], realized by placing a membrane cavity within a high-finesse Fabry-Pérot cavity. Ref. [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. Ref. [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 Ref. [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 Refs. [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 Sec. 2 we provide the basic theoretical description of the multimode optomechanical system under study. In Sec. 3 we describe the experimental setup, and in Sec. 4 we derive in detail the probe beam power spectral density, including all the non-linear effects. In Sec. 5 we analyze 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. Sec. 6 is for concluding remarks.

Theoretical description of the system dynamics
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, where ∆ (0) i = ω Li − ω ci are the detunings, E i = 2κ in,i P i /hω 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 zpf,j = h/2m j ω j the spatial width of the j-th oscillator zero point motion, and β in j (t), and α opt i (t), are the mechanical and optical noise terms, respectively. These noises are uncorrelated from each other and the only nonzero correlation functions are β in, * j (t)β in j (t ) = (n j + 1/2)δ jj δ(t − t ), wherē n j = [exp (hω j /k B T ) − 1] −1 k B T /hω j 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,31,32,33,34,35], which may eventually become synchronized. Synchronization may occur on a long timescale, determined by the inverse of the typically small parameters ∆ω = ω 2 − ω 1 (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)-(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 Ref. [5] to the case with noise studied here, as discussed in detail in Ref. [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 where β 0,j is the approximately constant, static shift of the j−th resonator, A j (t) is the corresponding slowly-varying complex amplitudes, and ω ref ∆ω is a reference mechanical frequency, of the order of ω j .
From eqs. (1)-(2) one gets the set of coupled amplitude equations [45] with ∆ω j = ω j − ω ref , and where are non-linear coefficients because of their dependence upon the regular dimensionless auxiliary functions F i , which are given by and which can be easily shown to be a function of even powers of |A b i | only. We have defined the bright complex amplitudes and J n is the n-th Bessel function of the first kind. η opt i is the term describing the noise of optical origin [45]. As already shown in Refs. [5,45], eqs. (4)-(5) provide a general and very accurate description of the dynamics of the two mechanical resonators.

Experimental setup
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 Fig. 1. A laser beam at wavelength λ 0 = 1064 nm is split in a probe beam with intensity Experimental setup for studying the non-linear dynamics in an optomechanical system constituted of a two-membrane sandwich within a cavity. A probe beam, frequency modulated by an electro-optical modulator (EOM), impinges on the optical cavity. The reflected beam is split: one component is detected, demodulated and low-pass amplified for generating the Pound-Drever-Hall (PDH) error signal able to lock the laser to the cavity; the second component is analyzed by homodyne detection in order to detect the mechanical motion. A further beam, the pump beam, detuned by ∆ 1 from the cavity resonance by means of an acousto-optic modulator (AOM), is turned on for engineering the optomechanical interaction, and in particular to realize laser driving of the mechanical modes. HWP denotes a half-waveplate, QWP a quarter-waveplate, BS a beam-splitter, and PBS a polarizing beam-splitter. P probe = 5.9 µW, modulated by an electro-optical modulator (EOM), and a pump beam, detuned by ∆ 1 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.
The optical and mechanical properties of the optomechanical system were investigated in Ref. [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 3 N 4 given in Ref. [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 3 , the side lengths were estimated to be L (1) (6) mm. We studied the dynamics of the lower frequency mode of the two membranes: for the first membrane we measured ω 1 2π × 230.795 kHz, γ 1 2π × 1.64 Hz, while for the second ω 2 2π × 233.759 kHz, γ 2 2π × 9.37 Hz. The membrane-cavity is placed in the middle of an optical cavity with empty cavity finesse F 0 = 50 125 (25), which reduces to F = 12 463(13) when the membrane-sandwich is placed in. Such finesse corresponds to a cavity intensity decay rate 2κ = τ −1 = FSR/F 2π × 134 kHz, with FSR 2π × 1.67 GHz.

Spectral analysis of the probe beam
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].
In Fig. 2 are compared simulated and experimental results for a set of parameters in this weak driving regime. On panel b) is reported the voltage spectral noise (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, ω 1 , (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 1 2π × 0.43 Hz, and g 2 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 ∆ 1 /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 Variances corresponding to the frequency ranges highlighted by the squares symbols: orange curve for the first mode; dark-green curve for the second mode; light-green curve for the sideband at 2ω 1 −ω 2 ; magenta curve for the calibration tone at ω b = 2π ×226.000 kHz. Residual detuning beat tone at ω det = 2π ×259.350 kHz is outside the displayed frequency range. Areas A and A indicate the steady-state regime reached for the two power settings.  2). We note the remarkable agreement between simulated and experimental data. Results in Fig. 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ω 1 − ω 2 .
For a detailed quantitative description of these results, we need to reconsider the non-linear dynamics of the mean cavity-field amplitude α 2 of the probe beam described by eq. (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 α 2 → α, κ 2 → κ and so on. We geṫ where , and e in = P in /hω L . We assume that the behaviour of the two oscillators is also sinusoidal as in eq.
, and that the first oscillator reaches a limit cycle with an amplitude |A 1 | for which g 1 |A 1 | is much larger that g 2 |A 2 |. The solution of eq. (10) can be found considering an expansion in terms of = g 2 |A 2 |/g 1 |A 1 |, that is α = j j α j . The zero-order solution, α 0 , satisfieṡ , and a first-order perturbation solution α 1 , driven by the amplitude of the second oscillator, Solution of eq. (11) provides the leading-order contribution to the cavity response where ξ = 2g 1 |A 1 |/ω 1 , and with W = i∆ − κ. The first-order solution can be found to be where Finally we observe that C is composed by a series of sidebands at frequencies nω 1 +bω b ± ω 2 , which are reproduced on the reflected output field determined by the input-output relation e out = −e in exp[−iβ sin(ω b t)] + √ 2κ in α. The cavity reflection function, R = e out /e in = j j R j , has a leading-order contribution We observe that R 0 is the superposition of the cavity reflection function R b (ξ) at each sideband frequencies bω b of the input field The first-order contribution is We consider now the modulation amplitude β, and the amplitude of the second mode, as small perturbations, that is β 1, and R R 0 + R 1 . In particular we focus on the contributions at DC, and at the five frequencies: ω 1 , ω 2 , ω sm = 2ω 1 − ω 2 , ω b , and In this case we have The contribution at DC is provided by eq. (17) for b = 0, and the time-independent term (n = 0) of eq. (14) at ω 1 is given by eq. (17) with b = 0, and in eq. (14) the term with n = ±1 at ω 2 is provided by eq. (15) with b = 0, and in eq. (16) the term with n = 0 for ω sm = 2ω 1 − ω 2 , b = 0 in eq. (15), and n = ∓2 in eq. (16) at ω b is provided by eq. (17) with b = ±1, and in eq. (14) the term with n = 0 and the sideband at ω sb = 2ω 1 − ω b for b = ±1, and n = 0 The homodyne technique is implemented by mixing on a beam-splitter the reflected field e out with an intense local oscillator e lo e iφ lo , and detecting the fields at the output of the beam-splitter, i.e., e j = [e lo e iφ lo + (−1) j e out ]/ √ 2. The output power are P ± = |e lo e iφ lo ± e out | 2 /2, and the differential current 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 φ lo = arctan [Re(R DC )/Im(R DC )]. This phase is plotted in Fig. 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  29), and the experimental parameters measured for the results reported in Fig. 2, considering a small nonzero probe detuning ∆ ∼ 2 × 3.9 kHz, β = 2 × 10 −2 rad, and the dimensionless amplitude A 2 = k B T /mω 2 2 /2 x zpf ∼ 3634. Black curve is the DC contribution (21); orange curve at ω 1 , eq. (22); dark-green curve at ω 2 , eq. (23); light-green curve at ω sm = 2ω 1 − ω 2 , eq. (24); magenta curve at ω b , eq. (25); blue line curve at ω sb = 2ω 1 − ω b , eq. (26). Black dashed curve is the local oscillator phase, φ lo , which, for our experimental parameters, can be considered π/2.
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 S W (ω) = dτ e iωτ V H (t+ τ )V H (t) t /R 0 on a termination resistor R 0 , provides a normalized amplitude spectral noise for each well separated frequency, that is with S 0 = (g T 2S) 2 P lo P in /R 0 . In Fig. 3 are reported the amplitude spectral noise at the frequencies of interest for our experiment. We notice that for probe detuning ∆ = 0, and when the cavity field is weakly modulated at frequency ω 1 , that is ξ 1, the output signal is linear in the displacement and the amplitude spectral noise is where F = FSR/2κ, ξ = δx(ω 1 )g 1 /ω 1 x zpf , and η = [2κ in /κ] · [g 1 λ 0 /2 FSR x zpf ] ∼ 0.25 · 0.24 ∼ 0.06 for our setup. The average power falling on each photodiode is approximately P lo /2. The shot noise in the differential signal has a flat spectrum with spectral density S sn PP ∼ 2 × 2hω L P lo /2, which sets a limit to the sensitivity of the detection [67,68]  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 λ 0 /2F 43 pm, which corresponds to ξ cav = 2g 1 /ω 1 · λ 0 /2Fx zpf 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 Fig. 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 Fig. 3). In general, two main ratios of the normalized spectral amplitudes V H (Ω), which are independent from β, κ in , and from the coupling term with the second mechanical mode, g 2 |A 2 |, 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 eq. (22) and the product of β and eq. (25); ii) the ratio derived by eq. (23) and eq. (24) In panel a) of Fig. 4 are reported these ratios that allow us to estimate ξ for our experimental realisations. In panel b) the ratios T 1 (ξ)ξ −1 , and T 2 (ξ) = V H (ω 2 )·β/V H (ω b ) (note that T 2 does depend on g 2 |A 2 |), normalized to their maximum values, which are reached for ξ → 0, are shown as N 1 , and 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. N 1 , and 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 Fig. 2, and Fig. 7 we estimate ξ st 1 1.05, N 1 0.70 and N 2 0.42. Such analysis allows us to deduce the effective limit cycle displacement amplitude of the first oscillator to be q st 1 = ξ st 1 ω 1 x zpf /g 1 263 pm, and the observed limit cycle displacement amplitude q ob 1 = q st 1 N 1 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 N 2 .

Non-linear mechanical dynamics at the onset of synchronization
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 Fig. 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 eqs. (1)-(2) with the parameters reported in Table 1, which refers to the results of Fig. 2, also shows that when the output probe beam exhibits sidebands around the frequency of the fundamental mode of the first oscillator [see Fig. 2a)], 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 Fig. 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]  where θ j (t) = arg[β j (t)], reported in Fig. 5c), 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. We are also able to provide a consistent analytical description of this presynchronization dynamics of the two membrane modes starting from the slowly varying amplitude equations (4)- (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 = ω 1 as a reference in eqs. (4)-(5), so that ∆ω 1 = 0 and ∆ω 2 = ∆ω. One can make quantitative predictions on such a regime assuming that |A 1 | |A 2 |, √ 2n 1 . Moreover, in our experiment the optical noise is negligible and we will not consider the terms associated with η opt i (t). With the above approximations, one can neglect both thermal noise and A 2 contributions in Eq. (4), which becomeṡ where we have made explicit the dependence of d 1 on |A 1 |, and defined ∆ω ef f The effective mechanical damping can be cast as where we have used the fact that in the considered regime |A b i | |A 1 |(g i1 /g b i ), ξ j ξ 1 = 2g 1 |A 1 |/ω 1 , assumed g i1 g 1 , and 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 2 A 2 would have not been negligible anymore with respect to g 1 A 1 , and the dynamics would have been governed by the more general eqs. (4)-(5) [55,57]. Eq. (35) can be solved by rewriting it in terms of modulus and phase, A 1 = I 1 e iφ 1 , After a transient these equations yield a steady state with a constant radius of the limit cycle of the first oscillator, I st 1 , corresponding in our case of not too strong driving, to the smallest positive root of the implicit equation γ ef f 1 (I st 1 ) = 0, which can be cast as a |ξ st with ξ st 1 = 2g 1 I st 1 /ω 1 , and a = ω 1 γ 1 /2g 2 1 . As a consequence, at long times, φ st Fig. 6 we show the left and right side of eq. (40) for the experimental parameters of Table. 1, which provides the optomechanical parameters for the results reported in Fig. 2, and Fig. 7. We infer from the intersection point, which corresponds to find the smallest root of γ ef f 1 (I st 1 ) = 0, a value ξ st 1 = 1.054, a steady-state displacement Table 1. Optomechanical parameters for the results reported in Fig. 8.
2π · 3.9 kHz g 1 2π · 0.4225 Hz g 2 2π · 0.6965 Hz γ 1 2π · 1.64 Hz γ 2 2π · 9.37 Hz κ loss 2π · 50.35 kHz κ in 2π · 8.35 kHz κ ex 2π · 58.7 kHz FSR 2π · 1.67 GHz P pump 4.25 µW P probe 5.9 µW λ 0 Figure 6. Steady-state solution for the mechanical displacement amplitude that reaches a limit cycle. The intersection of right and left side of eq. (40) for the parameters of Table 1, which are reported as solid blue, and dashed orange curves, respectively, determines the steady-state value of ξ st 1 . We determine ξ st 1 = 1.054, and an effective steady-state amplitude q st 1 = 2|A 1 |x zpf = 263.0 pm, which confirms the results shown in Fig. 4. The vertical black dashed, and dot-dashed lines represent the values for the thermal displacement q th = k B T /m 1 ω 2 1 3.365 pm corresponding to ξ th 0.0112, and for the cavity response length λ 0 /2F 43 pm, corresponding to ξ cav 0.284, respectively. The oblique light-green line represents the left term in eq. (40) for the second less coupled mode, for which the equation is not satisfied for the parameters in Table 1. In fact, the threshold power, that is the minimum power for finding a root, equivalently for the optical damping to exceed the intrinsic one, for the first mode is ∼ 3 µW, while for the second 6.75 µW. For power larger than 6.75 µW both modes might establish a limit cycle [57]. The oblique black dashed line indicates the boundary between the region with only one solution and with multiple solutions, that is the multistability parameter region, which occurs for a pump power larger than ∼ 667 µW.
amplitude q st 1 = 2|A 1 |x zpf = 263.0 pm, and ∆ω ef f 1 = −2π · 0.04 Hz, confirming the value of q st 1 deduced in Fig. 4. We emphasise once more that, as shown in Fig. 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 1 ; on the contrary for ξ ≥ 1 linearity is no more valid and a proper correction factor should be considered. In our case, as obtained in Fig. 4, the theoretical correction factor is N 1 0.70, corresponding to an expected observable stationary limit cycle amplitude of q ob 1 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, eq. (40) may have more than one solution at sufficiently large power (see below the oblique black dashed line in Fig. 6), corresponding to the multistability phenomenon analysed in Ref. [32] and experimentally verified in Ref. [34].
We now consider the dynamics of the second oscillator inserting the steady-state solution for A 1 (t) into Eq. (5), which becomeṡ The stationary solution can be obtained via Fourier transform and it can be written as where γ ef f 2 = γ 2 +Im[d 2 (I st 1 )] is positive, i.e., the second resonator is still damped despite the pump driving and it is not driven into a limit cycle, ∆ω ef f . This transition to synchronization is consistent with the theoretical analysis made in Refs. [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 [37].
We corroborate such analysis by considering the experimental time traces shown in  Table. 1.
frequencies ω 1 , and ω 2 , as a function of time. In panel c) and d) are reported the associated phase-space distributions. Panel 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 Fig. 8b), and Fig. 8c), 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 ∆q ∆ where, for almost resonant field [1], This fact allows us to estimate a small blue detuning of the probe ∆ 2 2π · 3.9 kHz, which is the value provided in Table 1.
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 ob 1 184 pm, orange curve in Fig. 8c). Such value agrees very well with the expected one, shown as dark-orange filled circle in panel a) of Fig. 8. Panel a) represents the effective steadystate mechanical amplitudes obtained as solution of eq. (40) (blue curves), and the observed one (orange curves), that is, reduced by the correction factor reported in Fig. 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 eq. (4) [and eq. (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 Fig. 8a). The effective steady-state displacement amplitude of the first oscillator, q st 1 = 2|A 1 |x zpf 262 pm, is also confirmed by the 10 blue trajectories simulated with the parameters of Table 1, and reported in Fig. 8c). 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 2 [green curve in Fig. 8b)] follows the dynamics of the effective mechanical displacement [blue trajectories in Fig. 8b)] only until q 1 reaches the limit cycle (after 7 s), a) Observed (orange curve) and effective (blue curve), steady-state amplitudes of the first mode, q st sm , as a function of the single-photon optomechanical coupling rate of the first mode g ≡ g 1 . Solid, and dashed curves represent the solution of eq. (40) for two pump settings: P probe = 5.9 µW, P pump = 4.25 µW, and P probe = 16 µW, P pump = 18.7 µW, respectively. Filled diamonds, which confirm the results obtained by finding the first zero of eq. (40) for the second power setting, represent the effective steady-state amplitudes obtained as the mean of the amplitudes after the oscillator has reached the limit cycle, determined by integrating eq.
We note that, in our case, the effective amplitude is larger than the thermal one by a factor γ 2 /γ ef f 2 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 N 2 0.42, we estimate an observed displacement of q ob 2 = ∆q ∆ 2 · γ 2 /γ ef f 2 ·N 2 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 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.

Conclusion
We have presented a detailed experimental analysis of the dynamics of the multimode optomechanical system introduced in Ref. [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 focused onto the non-linear regime where a bluedetuned 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.