# Self-organized synchronization of mechanically coupled resonators based on optomechanics gain-loss balance

P. Djorwe, Yan Pennec, Bahram Djafari-Rouhani

► To cite this version:

P. Djorwe, Yan Pennec, Bahram Djafari-Rouhani. Self-organized synchronization of mechanically coupled resonators based on optomechanics gain-loss balance. Physical Review B, 2020, 102 (15), 155410, 10 p. 10.1103/PhysRevB.102.155410. hal-03321540

HAL Id: hal-03321540

https://hal.science/hal-03321540

Submitted on 31 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Self-organized synchronization of mechanically coupled resonators based on optomechanics gain-loss balance

P. Djorwé,¹,²,³,* Y. Pennec,¹,† and B. Djafari-Rouhani¹,‡

¹Institut d'Electronique, de Microélectronique et Nanotechnologie, UMR CNRS 8520 Université de Lille, Sciences et technologies, Villeneuve d'Ascq 59652, France

²Catalan Institute of Nanoscience and Nanotechnology (ICN2), BIST, Campus UAB, Bellaterra, 08193 Barcelona, Spain

³Department of Physics, Faculty of Science, The University of Ngaoundere, P.O. Box 454, Ngaoundere, Cameroon

<figure><img src="image_2.png" /></figure>

(Received 28 July 2020; revised 31 August 2020; accepted 2 September 2020; published 12 October 2020)

We investigate self-organized synchronization in a blue-detuned optomechanical cavity that is mechanically coupled to an undriven mechanical resonator. By controlling the strength of the driving field, we engineer a mechanical gain that balances the losses of the undriven resonator. This gain-loss balance corresponds to the threshold where both coupled mechanical resonators enter simultaneously into self-sustained limit cycle oscillations regime. This leads to rich sets of collective dynamics such as in-phase and out-of-phase synchro- nizations, depending on the mechanical coupling rate, the frequency mismatch between the resonators, and the external driving strength through the mechanical gain and the optical spring effect. Moreover, we show that the introduction of a quadratic coupling, which results from a quadratically coupling of the optical cavity mode to the mechanical displacement, enhances the in-phase synchronization. This work shows how phonon transfer can optomechanically induce synchronization in a coupled mechanical resonator array and opens up new avenues for phonon processing and novel memories concepts.

DOI: 10.1103/PhysRevB.102.155410

I. INTRODUCTION

Recent progress in nanoengineering has led to exotic op-
tomechanical crystals [1–3], which are able to confine several
optical and mechanical modes on a single chip. A lot of
attention has been paid to these optomechanical arrays owing
to their capabilities to promote new phenomena and appli-
cations. These include collective nonlinear dynamics [4,5],
quantum many-body dynamics of photons and phonons [6],
long-range coupling of phonon modes [7,8], photons and
phonons transport [9,10], Anderson localization [11], as well
as topological phases of sound and light [12,13].

Owing to their practical applications in rf communica-
tion [14], signal-processing [15], and clock synchronization
[16], collective phenomena such as synchronization [17–21]
and frequency locking [22,23] were recently realized in op-
tomechanical systems. These schemes are different one from
another, and each of them having its own specificity. Synchro-
nization of two oscillators in contact and sharing a common
optical mode was investigated in Ref. [17], while the same
configuration was later on extended to multiple resonators in
Ref. [18]. Another synchronized system, consisting of two
spatially separated oscillators integrated in a common optical
racetrack cavity was reported in Ref. [19]. Synchronization of
the mechanical dynamics of a pair of optomechanical crys-
tal cavities, which are intercoupled with a mechanical link

and support independent optical modes, has been reported in
Ref. [20]. The observation of antiphase synchronous oscil-
lations has been observed in this system, which is the first
reported work on mechanically coupled resonators in optome-
chanics. More recently, self-organized synchronization of two
optically coupled membranes placed inside an optomechan-
ical cavity has been demonstrated in Ref. [21]. It has been
shown that the system enters into the synchronization regime
via a torus birth bifurcation line. Also, the phase-locking phe-
nomenon and the transition between in-phase and antiphase
regimes were directly observed. Besides these optomechani-
cal synchronization investigations, frequency locking of two
and multiple optomechanical systems were reported, respec-
tively, in Refs. [22] and [23]. In all these works, light is
the key element used to couple the optomechanical systems
involved except in Ref. [20]. It follows that there is a lack
of work on optomechanical synchronization where phonons
are mediating the coupling between the oscillators. Theoret-
ical investigation of quantum many-body dynamics has been
investigated along this line in Ref. [6], where both photons
and phonons can hop to nearest neighbor sites, and each
optomechanical cell in the array is driven.

In the present work, we propose a blue-detuned optome-
chanical system that is mechanically coupled to an undriven
mechanical resonator. By driving the optomechanical system,
we engineer a mechanical gain that balances the losses of the
undriven oscillator. This gain-loss balance threshold induces
rich sets of collective nonlinear dynamics in the system. Our
proposal is different from those listed above in at least three
points. Indeed, (i) the collective dynamics are induced through
mechanical gain-loss balance, (ii) only phonons are mediating

<footnote>* philippe.djorwe@icn2.cat</footnote>

<footnote>†yan.pennec@univ-lille.fr</footnote>

<footnote>‡bahram.djafari-rouhani@univ-lille.fr</footnote>
the energy transfer between the driven and the undriven res-
onator, and (iii) only one optomechanical cell is excited and
the energy is transferred to the rest of the system. The figure
of merit of our proposal is to use only one driving laser to
drive a set of mechanically coupled mechanical resonators,
and our system can also be used to perform phonon transfer
in an array of mechanical resonators. Moreover, by adopting
the general case of nondegenerated mechanical resonators, we
have identified different sets of synchronized states emerging
in our proposal. These synchronized states depend on the
mechanical coupling, the frequency mismatch of the two res-
onators, and the external driving that controls the mechanical
gain and the optical spring effect. Both in-phase and out-of-
phase synchronization are present, and the phase transition
between these two phenomena happens at a phase flip which
has been shown to be related to a high enough value of
the mechanical gain. Moreover, we have used the quadratic
coupling to enhance the in-phase synchronization process. As
these rich collective phenomena can be controlled by only
tuning the driving field, our proposal appears as an efficient
platform to perform applications such as phonons transfer
and phonons processing. Our proposal can be implemented
in the experimental setup used in Ref. [20] providing that
one driving laser is switched off during the experiment. A
similar idea to our proposal has been recently implemented
to demonstrate an optomechanical transducer based on two
elastically coupled cantilevers [24]. In this experiment, one
of the cantilevers is trapped by a harmonically oscillating
optical field inside a fiber-based cavity, while the other one
has been left undriven. The elastic coupling between the two
cantilevers results from their connection to the same over-
hang. By switching on the optical driving, the two cantilevers
strongly hybridize enabling exchange and coherent transfer of
the phononic information between the two cantilevers. There-
fore, the features investigated in our proposal can be figured
out in Ref. [24] as well. The rest of the work is organized
as follows. In Sec. II, the model and its dynamical equations
are described. The emerged collective dynamics are presented
in Sec. III, without the quadratic coupling. The frequency
mismatch effect together with qualitative explanations of the
transitional phases are figured out in Sec. IV. Section V is de-
voted to the effect of quadratic coupling on the enhancement
of synchronization, and Sec. VI concludes the work.

II. MODELLING AND DYNAMICAL EQUATIONS

Our benchmark system is the one sketched in Fig. 1(a),
which consists of an optomechanical cavity whose vibrating
element is mechanically coupled to an undriven auxiliary me-
chanical resonator. In the rotating frame of the driving fields
(ω<sub>p</sub>), the Hamiltonian (ħ = 1) describing this system is

$$
H = H_{\text{OM}} + H_{\text{int}} + H_{\text{drive}}, \quad (1)
$$

with

$$
\begin{align}
H_{\text{OM}} &= -\Delta a^{\dagger} a - g a^{\dagger} a (b_{1}^{\dagger} + b_{1}) + \sum_{j=1,2} \omega_{j} b_{j}^{\dagger} b_{j}, \nonumber \\
H_{\text{int}} &= -J(b_{1}b_{2}^{\dagger} + b_{1}^{\dagger}b_{2}), \nonumber \\
H_{\text{drive}} &= E(a^{\dagger} + a). \tag{2}
\end{align}
$$

In this Hamiltonian, *a* and *b*<sub>*j*</sub> are the annihilation bosonic field operators describing the optical and mechanical resonators, respectively. The mechanical displacements *x*<sub>*j*</sub> are connected to *b*<sub>*j*</sub> as *x*<sub>*j*</sub> = *x*<sub>*ZPF*</sub>(*b*<sub>*j*</sub> + *b*<sub>*j*</sub><sup>†</sup>), where *x*<sub>*ZPF*</sub> is the zero-point fluctuation amplitude of the mechanical resonator. The mechanical frequency of the *j*th resonator is *ω*<sub>*j*</sub> and Δ = *ω*<sub>*p*</sub> − *ω*<sub>cav</sub> is the detuning between the optical drive (*ω*<sub>*p*</sub>) and the cavity (*ω*<sub>cav</sub>) frequencies. The mechanical coupling strength between the two mechanical resonators is *J*, and the optomechanical coupling is *g*. The amplitude of the driving pump is *E*. For a large photon number in the system, the quantum operators in Eq. (2) can be treated as complex scalar fields, which are defined as the mean values of the operators: 〈*a*〉 = *α* and 〈*b*<sub>*j*</sub>〉 = *β*<sub>*j*</sub>. This leads to the following set of classical nonlinear equations of motion for our system,

$$
\begin{align}
\dot{\alpha} &= \left[ i(\Delta + g(\beta_1^* + \beta_1)) - \frac{\kappa}{2} \right] \alpha - i\sqrt{\kappa} \alpha^{\text{in}}, \nonumber \\
\dot{\beta}_1 &= -\left( i\omega_1 + \frac{\gamma_1}{2} \right) \beta_1 + iJ\beta_2 + ig\alpha^*\alpha, \tag{3} \\
\dot{\beta}_2 &= -\left( i\omega_2 + \frac{\gamma_2}{2} \right) \beta_2 + iJ\beta_1, \nonumber
\end{align}
$$

where optical (κ) and mechanical (γ<sub>j</sub>) dissipations have been
added, and the amplitude of the driving pump has been sub-
stituted as E = √κα<sup>in</sup> in order to account for losses. In this
form, the input laser power P<sub>in</sub> acts through α<sup>in</sup> = √(P<sub>in</sub>/ħω<sub>p</sub>).
Throughout the work, we assume the hierarchy of parame-
ters γ<sub>j</sub>, g ≪ κ ≪ ω<sub>j</sub>, similar to those encountered in resolved
sideband experiment [25]. Moreover, our simulated results are
based on numerical integration of Eq. (3) by using fourth-
order Runge Kutta method on Matlab software. As we are
interested in the steady state solutions of the mechanical res-
onators, only the temporal traces in which all the transient
behaviors have died out are kept in our analysis.

Starting from two nondegenerated mechanical resonators,
ω₁ ≠ ω₂ and γ₁ ≠ γ₂, we aim to synchronize their motions.
Our strategy is based on engineering and control of the me-
chanical gain, by driving the optomechanical cavity with a
blue detuned electromagnetic field. This strategy is sketched
on Fig. 1(b) where the effective damping (γ<sub>eff</sub> = γ₁ + γ<sub>opt</sub>) is
plotted versus the driving field. As we increase the driving
strength, this generates optical damping or gain (γ<sub>opt</sub>) that
balances the intrinsic damping (γ₁) of the driven mechanical
resonator [see green dot in Fig. 1(b)] which starts emitting
phonons that are mechanically transferred to the second me-
chanical resonator. This phonon transfer process is revealed
in Fig. 1(c) (see the inset), which depicts the normalized
mechanical amplitudes of both resonators (see Appendix A)
versus the driving field. The blue (red) curve is the amplitude
of the driven (undriven) mechanical resonator. It can be seen
that the driven resonator (blue curve) is transferring phonons
to the undriven one (red), and this transfer process is more
pronounced when the optical damping balances the intrinsic
damping [compare green dots in Fig. 1(b) and Fig. 1(c)].
By further increasing the strength of the driving field, the
engineered gain increases too and balances the losses of the
second mechanical resonator [see magenta dot in Fig. 1(b)].
From this point, the whole system reaches the phonon lasing
<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Benchmark system consisting of a blue-detuned driven optomechanical cavity, which is mechanically coupled to an undriven mechanical resonator. (b) Effective damping $\gamma_{\text{eff}}$ versus the driving field $\alpha^{\text{in}}$, showing both phonon amplification thresholds for the driven resonator at $\gamma_1 = |\gamma_{\text{opt}}|$ (green dot) and the one for the whole system at $\gamma_2 = |\gamma_{\text{eff}}|$ (pink dot). (c) Normalized mechanical amplitudes of the driven (blue) and the undriven (red) mechanical resonator versus the driving field. The green dashed box is zoomed out in the inset to show phonon lasing of the driven resonator. The vertical double arrow is reminiscent of a phase-flip transition that will be explained in Secs. III and IV. (d) Diagram displaying collective dynamics in the $(\alpha^{\text{in}}, J)$ parameter's space. Once the losses are balanced by the gain, both mechanical resonators carry out different types of dynamics depending on the values of $\alpha^{\text{in}}$ and $J$. The coupling has been fixed at $J = 2.2 \times 10^{-2}\omega_m$ for (b) and (c). The other used parameters are, $\omega_1 = \omega_m$, $\omega_2 = 1.002\omega_m$, $\gamma_1 = 10^{-3}\omega_m$, $\gamma_2 = \frac{\gamma_1}{3}$, $\kappa = 10^{-1}\omega_m$, $g = 2.5 \times 10^{-4}\omega_m$, $\Delta = \omega_m$.</figcaption></figure>

threshold and both mechanical resonators simultaneously start emitting phonons [see magenta dot in Fig. 1(c)]. After this amplification phase, the system settles into a self-sustained mechanical oscillations regime, above which complex nonlinear behaviors such as period doubling and chaos could emerge for strong enough driving strength [26]. As we are looking for collective dynamics such as synchronization and frequency locking phenomena, we limit ourselves in this work to the self-sustained oscillations regime. In order to optimize the input power needed to balance gain and losses in the whole system, we require that the second resonator has a higher quality factor than the first one ($\gamma_1 > \gamma_2$). This requirement ensures that our system can reach low-power phonon lasing, which is a useful prerequisite for experimental test of our proposal. Indeed, the sooner $\gamma_1 = |\gamma_{\text{opt}}|$ the faster the overall losses are balanced in the whole system by supplying less driven power [compare the green and pink dots in Fig. 1(b) and Fig. 1(c)].

### III. COLLECTIVE DYNAMICAL STATES

Beyond the phonon lasing threshold, both mechanical resonators exhibit limit cycle oscillations, which are correlated

through their phases and/or their vibrational amplitudes. These correlations lead to different sets of collective nonlinear behaviors, ranging from in-phase to out-of-phase synchro- nizations with equal or mismatched amplitudes. We have characterized these behaviors from steady state solutions of the mechanical resonators, where all the transient behavior has died out. The overall dynamics is shown in Fig. 1(d). The blue area is the linear regime, which is characterized by fixed points and where there is an energy flow from the driven resonator ($\beta_1$) into the undriven one ($\beta_2$) as depicted in the inset of Fig. 1(c). At the gain and loss balance, the paramet- ric instability threshold is reached, and phonon amplification process happens (gray area) until the mechanical amplitudes settle into the self-sustained oscillations regime due to the intrinsic nonlinearity in the system. Within the self-sustained oscillations regime, and depending on the strength of both mechanical coupling $J$ and driving field $\alpha^{\text{in}}$, the resonators display several sets of collective dynamical state as shown by different colors beyond the gray area in Fig. 1(d). Two quantities have been defined to characterize these dynamical states, the phase difference between the resonators and their mismatched amplitudes that is termed error synchronization later on. The instantaneous phase $\phi_i$ of a given resonator is
<figure><img src="image_4.png" /><figcaption>FIG. 2. (a), (b) 0, π-synchronization states. (c), (d) Phase, antiphase synchronizations. We highlight that (a)–(d) are related to the yellow dots A, B, C, and D indicated in Fig. 1(d). In (a) and (b), $J = 2.2 \times 10^{-2}\omega_m$ [see upper horizontal dashed line in Fig. 1(d)] and $\alpha^{\text{in}} = (1.2 \times 10^2, 2.2 \times 10^2)\sqrt{\omega_m}$, respectively. In (c) and (d), $J = 5 \times 10^{-3}\omega_m$ [see lower horizontal dashed line in Fig. 1(d)] and $\alpha^{\text{in}} = (10^2, 3 \times 10^2)\sqrt{\omega_m}$, respectively.</figcaption></figure>

defined as $\phi_i = \text{atan}(\frac{\text{Im}(\beta_i)}{\text{Re}(\beta_i)})$, where $\beta$ is a mechanical state variable. The averaged phase difference ($\Delta\phi$) between both resonator is

$$ \Delta\phi = \langle |\phi_i - \phi_j| \rangle \quad i, j = 1, 2, \quad (4) $$

where $\langle \cdot \rangle$ denotes an average over time. Similarly, the steady state amplitude of a given resonator is captured by $\epsilon \equiv \text{rms}(\beta_i)$, where “rms” is the root-mean-square value. We define the synchronization error here as a difference between the amplitudes of the resonators,

$$ \text{Err} = |\epsilon_i - \epsilon_j| \quad i, j = 1, 2. \quad (5) $$

Depending on these two quantities, the system exhibits four different dynamical states as depicted beyond the gray area in Fig. 1(d). In the green area, one has both $\Delta\phi \sim 0$ and $\text{Err} \sim 0$, resulting in a 0-synchronized state of the mechanical resonators as dynamically depicted in Fig. 2(a). It can be clearly seen that both resonators exhibit a similar behavior in this regime. In the red area however, the resonators are out of phase ($\Delta\phi \sim \pi$), where they oscillate with identical amplitudes ($\text{Err} \sim 0$) with a phase difference of $\pi$. Such a dynamic is depicted in Fig. 2(b), which reveals $\pi$ synchronization of both mechanical resonators over time. In the areas corresponding to a weaker mechanical coupling $J$, the resonators carry out only phase synchronizations. Indeed, phase synchronization is observed in the light-green area where both resonators are in phase ($\Delta\phi \sim 0$), but they have different amplitudes $\text{Err} \neq 0$. Such a synchronized state is dynamically shown in Fig. 2(c). Moreover, the resonators exhibit an antiphase synchronization in the magenta area, where they evolve out of phase from each other with dissimilar vibrational amplitudes ($\text{Err} \neq 0$) as it can be seen in Fig. 2(d). These sets of dynamical states enabled by our proposal reveal its

<figure><img src="image_5.png" /><figcaption>FIG. 3. Phase difference $\Delta\phi$ (dash-dotted curve) and error synchronization $\text{Err}$ (full curve) between the mechanical resonators for $J = 2.2 \times 10^{-2}\omega_m$. The other parameters are the same as in Fig. 1.</figcaption></figure>

performance in carrying out rich collective dynamics compared to those known in the state of the art optomechanical systems. The key point of this performance here is the presence of gain and loss, which induces the different phase transitions observed in Fig. 1(d) as we will explain later on. From these displayed collective dynamics, one can deduce two features depending on either the mechanical coupling $J$ or the driven strength $\alpha^{\text{in}}$ is adjusted. For weak mechanical coupling for instance, the two mechanical resonators start exhibiting phase synchronization, which ends up to complete synchronization as $J$ is increasing. This is the case for the transition from the light green to the green areas in Fig. 1(d). Similar transition happens when the resonators are out of phase as well, where the increase of $J$ adjusts their vibrational amplitudes [see switching from magenta to red areas in Fig. 1(d)]. One can point out also the fact that the driving strength mostly induces switching related to the phase of the resonators. Indeed, as $\alpha^{\text{in}}$ is increasing, the resonators switch from being in phase ($\Delta\phi \sim 0$) to completely be out of phase ($\Delta\phi \sim \pi$). This feature can be seen for instance from the transitions between the green and the red areas in Fig. 1(d). This switching that is termed later on as phase flip is a well-known nonlinear phenomenon characterized by a sudden jump of the phase difference roughly from 0 to $\pi$ [27,28]. We have figured out this phenomenon in Fig. 3, where the phase difference given in Eq. (4) (dash-dotted curve) and the synchronization error defined in Eq. (5) (full curve) are represented for a given $J = 2.2 \times 10^{-2}\omega_m$. It can be clearly seen that after the amplification regime, both mechanical resonators adjust their amplitude and phase, and synchronize for a while within the range ($100 \lesssim \alpha^{\text{in}} \lesssim 135)\sqrt{\omega_m}$. Beyond the upper limit of this interval, the phase difference suddenly jumps from roughly 0 to $\pi$ while the synchronization error remains almost zero ($\text{Err} \sim 0$). This phase-flip phenomenon marks the threshold of the out of phase synchronization, which is induced both by gain-loss competition and optical spring effect through an adjustment of the driving strength. We have found that the phase-flip transition is likewise revealed through the amplitude jump phenomenon, pointed out by the green double arrow indicated in Fig. 1(c). This can be qualitatively explained by an analytical approach based on the mechanical eigenmodes of our system.
IV. ANALYTICAL APPROXIMATIONS AND FREQUENCY MISMATCH EFFECT

Based on well-known analytical approximations (see Appendix A), the mechanical eigenmodes of our system can be obtained by integrating the optical intracavity field $\alpha(t)$ out of the set of Eq. (3). The resulting eigenmodes yield

$$ \lambda_{\pm} = \frac{\omega_{\text{eff}} + \omega_2}{2} - \frac{i}{4}(\gamma_{\text{eff}} + \gamma_2) \pm \frac{\sigma}{4}, \quad (6) $$

where $\omega_{\text{eff}} = \omega_1 + \delta\omega_{\text{opt}}$ and $\gamma_{\text{eff}} = \gamma_1 + \gamma_{\text{opt}}$ are the effective frequency and damping of the driven mechanical resonator, respectively. The quantities $\delta\omega_{\text{opt}}$ and $\gamma_{\text{opt}}$ are, respectively, the optical spring effect and damping induced by the driving field, and they are given by

$$ \delta\omega_{\text{opt}} = -\frac{2\kappa(g\alpha^{\text{in}})^2}{\omega_{\text{lock}}\epsilon_1} \text{Re}\left(\sum_n \frac{J_{n+1}(-\epsilon_1)J_n(-\epsilon_1)}{h_{n+1}^*h_n}\right), \quad (7) $$

and

$$ \gamma_{\text{opt}} = \frac{2(g\kappa\alpha^{\text{in}})^2}{\epsilon_1} \sum_n \frac{J_{n+1}(-\epsilon_1)J_n(-\epsilon_1)}{|h_{n+1}^*h_n|^2}, \quad (8) $$

where $\epsilon_1 = \frac{2g\text{Re}(A_1)}{\omega_{\text{lock}}}$ is a normalized amplitude, with $A_1$ and $\omega_{\text{lock}}$ the mechanical amplitude of the first resonator and its locked frequency, respectively. $J_n$ is the Bessel function of the first kind of order $n$, $h_n = i(n\omega_{\text{lock}} - \bar{\Delta}) + \frac{\kappa}{2}$, and $\bar{\Delta} = \Delta_1 + 2g(\bar{\beta}_1)$ the effective detuning. In Eq. (6), $\sigma$ is a complex quantity defined as

$$ \sigma = \sqrt{16J^2 + [2(\omega_{\text{eff}} - \omega_2) + i(\gamma_{\text{eff}} - \gamma_2)]^2}, \quad (9) $$

whose real and imaginary part read

$$ \text{Re}(\sigma) = \sqrt{\frac{\chi^2 + (4\Delta\omega_{\text{eff}}\Delta\gamma_{\text{eff}})^2 + \chi}{2}}, \quad (10) $$

and

$$ \text{Im}(\sigma) = \sqrt{\frac{\chi^2 + (4\Delta\omega_{\text{eff}}\Delta\gamma_{\text{eff}})^2 - \chi}{2}}, \quad (11) $$

where $\chi = 16J^2 + 4\Delta\omega_{\text{eff}}^2 - \Delta\gamma_{\text{eff}}^2$ with $\Delta\omega_{\text{eff}} = \omega_{\text{eff}} - \omega_2$ and $\Delta\gamma_{\text{eff}} = \gamma_2 - \gamma_{\text{eff}}$. From Eqs. (6)–(11), it appears that the mechanical supermodes can be characterized by their frequencies $[\text{Re}(\lambda_{\pm}) = \omega_{\pm}]$ and dampings $[\text{Im}(\lambda_{\pm}) = \gamma_{\pm}]$ given by

$$ \omega_{\pm} = \frac{\omega_{\text{eff}} + \omega_2}{2} \pm \frac{\text{Re}(\sigma)}{4}, \quad (12) $$

and

$$ \gamma_{\pm} = -\frac{(\gamma_{\text{eff}} + \gamma_2)}{4} \pm \frac{\text{Im}(\sigma)}{4}. \quad (13) $$

In order to reveal the mechanical supermodes involved in our study, we fast Fourier transform (FFT) the steady state solutions of the temporal traces of each mechanical resonator for different $\alpha^{\text{in}}$, where all the transient behavior has died out. From this FFT data, we collect the dominant peaks together with their corresponding frequencies, which represent the amplitudes and the frequencies of the supermodes, respectively. Applying this procedure to the data of Fig. 1(c) for instance, has led to the supermodes depicted in Fig. 4.

<figure><img src="image_6.png" /><figcaption>FIG. 4. Mechanical supermodes obtained by fast Fourier transforming the steady state solutions of the temporal traces of Fig. 1(c). (a) Normalized mechanical amplitudes ($\epsilon_{1,2}$) of the resonators. These amplitudes can be viewed as two upper ($\epsilon_+$) and two lower supermode amplitudes with the corresponding frequencies $\omega_{\pm}$. (b) Supermode frequencies ($\omega_{\pm}$). The black color corresponds to the supermode oscillating with the highest amplitude, and the gray color is related to the supermode oscillating with the lower amplitude (see further details in Fig. 5). The vertical lines correspond to: (i) the phonon lasing of the driven resonator (green line), (ii) the phonon lasing in the whole system (magenta line), and (iii) the phase-flip transition (cyan line). These curves are plotted versus $\alpha^{\text{in}}$ for a fixed $J = 2.2 \times 10^{-2}\omega_m$, and the other parameters are as in Fig. 1.</figcaption></figure>

In Fig. 4(a), it can be clearly seen that each mechanical resonator is characterized by two supermodes ($\epsilon_{\pm}$) which are well captured through their supermode frequencies ($\omega_{\pm}$) in Fig. 4(b). More importantly, the lower supermode amplitudes ($\epsilon_-$) died out exactly where the phase-flip and amplitude jump phenomena happen [see green rectangle in Fig. 4(a)]. Furthermore, the supermode frequencies merge to a single frequency ($\omega_{\text{lock}}$) at the same transitional point [Fig. 4(b)], which clearly reveals a frequency locking effect and can be explained by the mode competition phenomenon as well [21].

From Eqs. (10)–(13) together with the frequency dynamics figured out in Fig. 4(b), one can analytically predict the transitional phases involved in our system. These transitions
are indicated by the vertical lines in Fig. 4(b) and they refer to the phonon lasing threshold phenomena (green and magenta lines) and the phase-flip effect (cyan line). Since the first transition (green line) is not directly related to the synchronization feature, the specific cases of the global phonon lasing (magenta line) and phase flip are detailed in the following. For degenerated mechanical resonators for instance ($\omega_1 = \omega_2$), one gets $\Delta\omega_{eff} \sim 0$ for a weak driving field. At the balanced gain and loss, it results that $\text{Im}(\sigma) = 0$ and $\text{Re}(\sigma) = \sqrt{16J^2 - \Delta\gamma_{eff}^2}$. This leads to a well-known exceptional point (EP) at $J = \frac{\Delta\gamma_{eff}}{4}$, which can be used for mass sensing [29] and it also marks the threshold of phonon lasing where the mechanical resonators exhibit self-sustained limit cycle oscillations [26].

In our proposal however, there is a frequency mismatch $\delta\omega = \omega_2 - \omega_1$, which acts as a perturbation on the EP. At the balanced gain-loss, the requirement of EP ($\text{Re}(\sigma) = \text{Im}(\sigma) = 0$) is fulfilled in our system providing that $\chi = 0$ and $\Delta\omega_{eff} \Delta\gamma_{eff} = 0$. The former condition is satisfied for $4J^2 + \Delta\omega_{eff}^2 = \gamma_2^2$, while the latter condition cannot be satisfied at this weak driving field regime, since $\Delta\gamma_{eff} \simeq 2\gamma_2$ and $\Delta\omega_{eff} = \delta\omega_{opt} - \delta\omega$. Therefore, one obtains

$$ \omega_{\pm} = \frac{\omega_{\text{eff}} + \omega_2}{2} \pm \sqrt{\frac{\Delta\omega_{\text{eff}}\Delta\gamma_{\text{eff}}}{8}} \\ \approx \frac{1}{2}(\omega_{\text{eff}} + \omega_2 \pm \sqrt{\gamma_2\Delta\omega_{\text{eff}}}), \quad (14) $$

and

$$ \gamma_{\pm} = -\frac{(\gamma_{\text{eff}} + \gamma_2)}{4} \pm \sqrt{\frac{\Delta\omega_{\text{eff}}\Delta\gamma_{\text{eff}}}{8}} \\ \approx \frac{1}{2}\left(-\frac{(\gamma_{\text{eff}} + \gamma_2)}{2} \pm \sqrt{\gamma_2\Delta\omega_{\text{eff}}}\right) \quad (15) $$

which clearly reveals the effect of the frequency mismatch ($\delta\omega$) on getting $\text{Re}(\sigma) = \text{Im}(\sigma) = 0$. From Eq. (14) and Eq. (15), it results that two supermodes are still possible although the losses are balanced by the gain in the whole system as shown in Fig. 4. In order to reach the transition point where the phase flip happens, one needs to further increase the driving strength $\alpha^{\text{in}}$ (see Fig. 4). However, increasing $\alpha^{\text{in}}$ leads to the increase of both $\delta\omega_{opt}$ and $\gamma_{opt}$ as it appears in Eq. (7) and Eq. (8). For enough driving strength, one reaches $\delta\omega_{opt} = \delta\omega$ leading to $\Delta\omega_{eff} = 0$. Therefore, the previous unfulfilled condition ($\Delta\omega_{eff} \Delta\gamma_{eff} = 0$) is now satisfied and both Eq. (12) and Eq. (13) yield

$$ \omega_{\pm} = \frac{\omega_{\text{eff}} + \omega_2}{2}, \quad (16) $$

and

$$ \gamma_{\pm} = -\frac{(\gamma_{\text{eff}} + \gamma_2)}{4} \pm \sqrt{\chi_0}, \quad (17) $$

since $J < \frac{\Delta\gamma_{eff}}{4}$ at this strong driving regime and we have defined $\chi_0 = \Delta\gamma_{eff}^2 - 16J^2$. These equations reveal that one supermode has vanished [$\text{Re}(\sigma) = 0$] and this means that each mechanical resonator is now oscillating with a single frequency given by Eq. (16). From this phase-flip transition,

<figure><img src="image_7.png" /><figcaption>FIG. 5. Different spectra corresponding to the states within the areas delimited by vertical lines in Fig. 4. (a) Spectrum corresponding to $\alpha^{\text{in}} = 54\sqrt{\omega_m}$ which is an example of the situation found before the green vertical line. (b) Spectrum at $\alpha^{\text{in}} = 65\sqrt{\omega_m}$ which corresponds to a case encountered between the green and magenta vertical lines. (c) Spectrum for $\alpha^{\text{in}} = 90\sqrt{\omega_m}$ which reveals a situation found between the magenta and cyan vertical lines. (d) Spectrum which corresponds to $\alpha^{\text{in}} = 150\sqrt{\omega_m}$ and this shows the cases beyond the cyan vertical line. Despite the dynamical switching of the peaks between the frequencies $\omega_{\pm}$, attention should be also paid on the growing intensity from (a) to (d), which reveals the action of $\alpha^{\text{in}}$.</figcaption></figure>

each resonator carries out oscillations that decay according to one of the dissipations given in Eq. (17).

From the above discussion, it follows that both mechanical gain and optical spring are the key points for the collective phenomena arising in our proposal. Owing to the fact that these quantities are tuned through the driving field, our proposal appears as a nice platform to generate collective phenomena, where the system is wholly controlled externally. To bear in mind the dynamics of the supermodes, we have provided the spectra given in Fig. 5, which result from a strong hybridization of the mechanical resonators once the optical driving is switched on. These spectra capture the energy transfer between the supermodes for the different situations presented in Fig. 4(b). The spectrum in Fig. 5(a) depicts the dynamics before the green line, Fig. 5(b) shows the situation between the green and magenta lines, the case within the area bounded by the magenta and cyan line is given in Fig. 5(c), while the case beyond the cyan line is shown by the spectrum in Fig. 5(d). It can be seen that the energy is first transferred from the higher frequency supermode ($\omega_+$) to the lower one ($\omega_-$), and from the phase-flip transition, the whole energy is swapped back to the supermode $\omega_+$. This feature can be further exploited to perform phonon information processing between the mechanical resonators and to improve novel memory concepts. The rough idea is that the more synchronized the resonators are, the more information they can exchange. Thus, the study of synchronization and desynchronization of mechanical wave may be helpful to understand further the phonon information processing in mechanical structures. Furthermore, the classical idea of synchronization states that one can understand two synchronized oscillators by studying only one of them, since they present
<figure><img src="image_8.png" /><figcaption>FIG. 6. Numerical diagram displaying collective dynamics in the ($\alpha^{\text{in}}$, $\delta\omega$) parameter's space for $J = 2.2 \times 10^{-2}\omega_m$. The other parameters are as in Fig. 1.</figcaption></figure>

the same dynamical behavior. Therefore, synchronization enables us to store and to replicate a synchronized state in several samples for later usage, which is the idea behind the memory concepts.

Based on the fact that the frequency mismatch has been needed to explain the frequency locking effect through Re($\sigma$) = 0 or $\delta\omega_{\text{opt}} = \delta\omega$, a particular attention must be paid on it. Therefore, the effect of $\delta\omega$ on the organization of the collective dynamics is shown in Fig. 6, for a fixed coupling strength $J = 2.2 \times 10^{-2}\omega_m$. It can be seen that the frequency mismatch very weakly shifts the lasing threshold (transition between blue and gray areas), compared to the coupling strength $J$ [as in Fig. 1(d)]. This feature simply reveals the robustness of the amplification process at the vicinity of the EP, against moderate frequency mismatches as aforementioned. From Fig. 6, this robustness is roughly preserved up to 1% of frequency mismatch. However, just above the lasing threshold, the mismatch induces new reorganization of collective dynamics. For a weak mismatch ($\delta\omega < 2 \times 10^{-3}\omega_m$), $\pi$ synchronization is the dominant dynamics. This is expected and is explained by the fact that weak driving strength is required for the mechanical gain and optical spring effect to reverse the sign of $\chi$ and to compensate the mismatch, respectively. For large mismatch however, strong driving strength is required to generate large gain and optical spring effect. Therefore, it is hard to reach $\pi$ synchronization in such a configuration as it can be seen for $\delta\omega \gtrsim 5 \times 10^{-3}\omega_m$, where only phase synchronization emerges. For intermediate values of the mismatch ($2 \times 10^{-3}\omega_m \lesssim \delta\omega \lesssim 4 \times 10^{-3}\omega_m$), both resonators can synchronize for a while before switching to out of phase synchronization. Owing to the difficulties of engineering identical mechanical resonators in micro/nanofabrication technologies, it could be interesting to seek an enhancement strategy of the observed collective phenomena as presented in the next section.

## V. QUADRATIC COUPLING ENHANCES IN-PHASE DYNAMICS

In order to assist the role played by the optical spring effect and to enhance the effect of the mechanical gain, we

<figure><img src="image_9.png" /><figcaption>FIG. 7. Enhancement effect of the quadratic coupling on in-phase synchronization. (a) Phase difference $\Delta\phi$ (dash-dotted curve) and error synchronization Err (full curve) between the mechanical resonators for $J = 2.2 \times 10^{-2}\omega_m$. (b) Numerical diagram displaying collective dynamics in the ($\alpha^{\text{in}}$, $J$) parameter's space. Quadratic coupling of $g_{ck}/g = -10^{-3}$ has been accounted, and the other parameters are the same as in Fig. 1.</figcaption></figure>

have introduced the quadratic coupling ($g_{ck}$) in the system. This well-known nonlinearity can be generated by inserting the driven mechanical resonator inside the optical cavity for instance [30]. The resulted system can be thought of as a membrane-in-the-middle setup where the moving element is mechanically coupled to the undriven mechanical resonator. In such a system, the frequency of the driven mechanical system yields $\tilde{\omega}_1 = \omega_1 - g_{ck}\alpha^*\alpha$, which can be tuned through the driving field (see Appendix B). This frequency control is useful to bring closer the frequencies of the nondegenerate mechanical resonators involved in our system, even at the vicinity of an EP [29]. This can be seen by comparing Fig. 7(a) to Fig. 3. Indeed, 0 synchronization is established before $\alpha^{\text{in}} = 10^2\sqrt{\tilde{\omega}_m}$ in Fig. 7(a) (see dash-dotted curve) where the quadratic term is $g_{ck}/g = -10^{-3}$, whereas this happens above such a value in Fig. 3 where $g_{ck}/g = 0$. It can also be seen that 0 synchronization is wider in Fig. 7(a) than in Fig. 3, revealing the effect of the quadratic coupling in enhancing in-phase synchronization. This enhancement effect of the quadratic term is depicted in Fig. 7(b), where in-phase synchronization (green area) is clearly improved. Moreover, by comparing the diagrams depicted in Fig. 7(b) and Fig. 1(d), it results that quadratic coupling also reorganizes localization of the different dynamical states. We have checked the dynamics for the labeled points in Fig. 7(b),
<figure><img src="image_10.png" /><figcaption>FIG. 8. Dynamical states. (a) to (h) correspond to the points labeled A to H in Fig. 7(b).</figcaption></figure>

and the corresponding dynamical states satisfy the expected behaviors (see Fig. 8). Furthermore, the dynamical state of the gray area intertwined in the magenta zone has been figured out, and it corresponds to a strong coupling regime where both resonators exchange energy through Rabi oscillations (Appendix B).

## VI. CONCLUSION

We have investigated collective dynamics in a blue-detuned optomechanical cavity that is mechanically coupled to an undriven mechanical resonator. When the optomechanically engineered mechanical gain balances the losses of the undriven resonator, phonon lasing threshold happens and both resonators simultaneously exhibit self-sustained limit cycles, leading to interesting sets of collective dynamics. Depending on the mechanical coupling and the external driving field, we have observed in-phase and out-of-phase synchronizations that mainly result from the process of driving that induces both mechanical gain and optical spring. Qualitative explanations of the main phase transitions arising in our proposal have been provided, based on mechanical eigenvalues evaluated through analytical approximations. Furthermore, we have used quadratic optomechanical coupling to enhance in-phase synchronization between the nondegenerated mechanical resonators involved in the system. Our work opens new avenues towards collective dynamics in an array of mechanically coupled resonators by only driving one of them. This scheme can be extended to related systems including electromechanical and superconducting microwave setups.

## ACKNOWLEDGMENTS

This work was supported by the European Commission FET OPEN H2020 project PHENOMEN-Grant Agreement No. 713450. P.D. acknowledges the funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 754510, and the support from Severo Ochoa Program (MINECO, Grant No. SEV-2017-0706) and funding from the CERCA Programme/Generalitat de Catalunya.

## APPENDIX A: MECHANICAL SUPERMODES

Our aim here is to give few details about the steps leading to the mechanical supermodes mentioned in Sec. IV. To avoid reproducing some expressions that are already explained above, we start from Eq. (3) in the main text which is the classical nonlinear set of equations describing our proposal. Only one mechanical resonator ($\beta_1$) is driven by a blue-detuned electromagnetic field while the second ($\beta_2$) is just mechanically coupled to it. Therefore, the following optomechanical transformation will be only applied to $\beta_1$. However, we underline that both resonators can oscillate for enough driving strength, owing to the mechanical coupling $J$ as explained in the main text. Once the limit cycles regime is reached, this driven resonator carries out oscillations whose amplitudes change only slowly over time. Thus, we solve the equation for $\alpha$ assuming a fixed amplitude for the mechanical oscillations and then use the result to solve the equation for $\beta_1$. Under this assumption, the mechanical oscillation can be
described by the ansatz [26],

$$ \beta_j(t) = \bar{\beta}_j + A_j \exp(-i\omega_{\text{lock}}t), \quad (A1) $$

where $\bar{\beta}_j$ is a constant shift in the origin of the resonator and the amplitude $A_j$ is taken to be a slowly varying function of time. In such a regime of single frequency ansatz [see Fig. 4(b) for instance], we have denoted the locked frequency by $\omega_{\text{lock}}$. We substitute this ansatz into the equation for $\alpha$ and use the assumption of a slowly evolving amplitude to solve it, first neglecting the time dependence of $A_j$. We then obtain the intracavity field in the form,

$$ \alpha(t) = e^{-i\theta(t)} \sum_n \alpha_n e^{in\omega_{\text{lock}}t}. \quad (A2) $$

The phase is $\theta(t) = -\epsilon_1 \sin \omega_{\text{lock}}t$ and the amplitudes of the different harmonics of the optical field are

$$ \alpha_n = -i\sqrt{\kappa}\alpha^{\text{in}} \frac{J_n(-\epsilon_1)}{h_n}, \quad (A3) $$

where $\epsilon_1 = \frac{2g\text{Re}(A_1)}{\omega_{\text{lock}}}$, $\tilde{\Delta} = \Delta + 2g\text{Re}(\bar{\beta}_1)$, $h_n = i(n\omega_{\text{lock}} - \tilde{\Delta}) + \frac{\kappa}{2}$, and $J_n$ is the Bessel function of the first kind of order $n$.

As we are interested in the regime of limit cycles of the resonators, a rotating wave approximation can be made in which we drop all the terms (in the mechanical dynamics) except the constant one and the term oscillating at $\omega_{\text{lock}}$. Hence, we substitute Eq. (A2) in the equation for $\beta_1$ [see Eq. (3)] which, by equating constant terms, leads to the zero-frequency components,

$$ \bar{\beta}_1 = \frac{1}{\omega_1 - i \frac{\gamma_m}{2}} \left( g\kappa \sum_n \frac{(\alpha^{\text{in}} J_n(-\epsilon))^2}{|h_n|^2} + J\bar{\beta}_2 \right). \quad (A4) $$

that induce a shifts of the cavity frequencies,

$$ \delta = 2g\text{Re}(\bar{\beta}_1). \quad (A5) $$

The equations of motion for the oscillating part of $\beta_j$ are deduced from $\beta_r^j(t) = \beta_j(t) - \bar{\beta}_j \equiv A_j \exp(-i\omega_{\text{lock}}t)$ and read

$$ \begin{aligned} \dot{\beta}_r^1(t) &= -i(\omega_1 + \delta\omega_{\text{opt}})\beta_r^1 - \frac{\gamma_1 + \gamma_{\text{opt}}}{2}\beta_r^1 + iJ\beta_r^2 \\ \dot{\beta}_r^2(t) &= -i\omega_2\beta_r^2 - \frac{\gamma_2}{2}\beta_r^2 + iJ\beta_r^1. \end{aligned} \quad (A6) $$

Here the optical spring effect $\delta\omega_{\text{opt}}$ and the optical damping $\gamma_{\text{opt}}$ coming both from the average dynamics of the cavity are given by Eq. (7) and Eq. (8) in the main text, respectively. From Eq. (A6), it is possible to define an tive Hamiltonian in order to figure out supermodes involved in the system. Such supermodes will be deduced from the eigenmodes of the effective model, describing the mechanical resonators. Indeed, the real parts of the eigenmodes give the eigenfrequencies of the coupled system while their imaginary parts stand for the dissipations rate of the system. In the limit cycles regime, the constant shift $\bar{\beta}_j$ is weak compared to the amplitude of the mechanical resonator ($\bar{\beta}_j \ll A_j$). This means that $\beta_j(t) \cong \beta_r^j(t)$, and Eq. (A6) can be assumed as a set of equations describing the effective system that reads

$$ \begin{aligned} \dot{\beta}_1 &= -\left(i\omega_{\text{eff}} + \frac{\gamma_{\text{eff}}}{2}\right)\beta_1 + iJ\beta_2, \\ \dot{\beta}_2 &= -\left(i\omega_2 + \frac{\gamma_2}{2}\right)\beta_2 + iJ\beta_1, \end{aligned} \quad (A7) $$

where $\omega_{\text{eff}} = \omega_1 + \delta\omega_{\text{opt}}$ and $\gamma_{\text{eff}} = \gamma_1 + \gamma_{\text{opt}}$ define the effective frequencies and the effective damping, respectively.

Furthermore, Eq. (A7) can be rewritten in the compact form,

$$ \partial t \Psi = -i H_{\text{eff}} \Psi \quad (A8) $$

with the effective Hamiltonian,

$$ H_{\text{eff}} = \begin{bmatrix} \omega_{\text{eff}} - i \frac{\gamma_{\text{eff}}}{2} & -J \\ -J & \omega_2 - i \frac{\gamma_2}{2} \end{bmatrix} \quad (A9) $$

and the state vector $\Psi = (\beta_1, \beta_2)^T$. The eigenvalues of the Hamiltonian given in Eq. (A9) are obtained by solving the equation,

$$ \det(H_{\text{eff}} - \lambda I) = 0, \quad (A10) $$

and that yields to the following eigenvalues $\lambda_-$ and $\lambda_+$,

$$ \lambda_{\pm} \approx \frac{\omega_{\text{eff}} + \omega_2}{2} - \frac{i}{4}(\gamma_{\text{eff}} + \gamma_2) \pm \frac{\sigma}{2}, \quad (A11) $$

as given by Eq. (6) in the main text with the complex quantity $\sigma$ defined by Eq. (9) as well. The frequencies and the dissipations of the supermodes discussed in the main text are given by the real and imaginary parts of $\lambda_\pm$, respectively,

$$ \omega_{\pm} = \text{Re}(\lambda_{\pm}) \text{ and } \gamma_{\pm} = \text{Im}(\lambda_{\pm}). $$

## APPENDIX B: QUADRATIC COUPLING AND DYNAMICAL STATES

By taking into account the quadratic coupling (or cross-Kerr) in the model, the Hamiltonian becomes

$$ H = H_{\text{OM,ck}} + H_{\text{int}} + H_{\text{drive}}, \quad (B1) $$

with

$$ \begin{aligned} H_{\text{OM,ck}} &= -\Delta a^\dagger a + \sum_{j=1,2} \omega_j b_j^\dagger b_j \\ &\quad -ga^\dagger a(b_1^\dagger + b_1) - g_{ck}a^\dagger ab^\dagger b, \\ H_{\text{int}} &= -J(b_1 b_2^\dagger + b_1^\dagger b_2), \\ H_{\text{drive}} &= E(a^\dagger + a), \end{aligned} \quad (B2) $$

where ck stands for “cross-Kerr” and accounts for the quadratic coupling. This leads to the following classical set of nonlinear equations,

$$ \begin{aligned} \dot{\alpha} &= \left[ i(\Delta + g(\beta_1^* + \beta_1) + g_{ck}\beta_1^*\beta_1) - \frac{\kappa}{2} \right] \alpha - i\sqrt{\kappa} \alpha^{\text{in}}, \\ \dot{\beta}_1 &= -\left(i\tilde{\omega}_1 + \frac{\gamma_1}{2}\right)\beta_1 + iJ\beta_2 + ig\alpha^*\alpha, \\ \dot{\beta}_2 &= -\left(i\omega_2 + \frac{\gamma_2}{2}\right)\beta_2 + iJ\beta_1, \end{aligned} \quad (B3) $$

where $\tilde{\omega}_1 = \omega_1 - g_{ck}\alpha^*\alpha$ is the optically tunable mechanical frequency mentioned in the main text.

The cross-Kerr effect has been pointed out in the main text through Fig. 7, and it has been shown that it enhances in-phase synchronization. This enhancement is related to the control of the frequency mismatch through $\tilde{\omega}_1 = \omega_1 - g_{ck}\alpha^*\alpha$. As we
have the frequency hierarchy of $\omega_1 < \omega_2$, we have conveniently used a negative quadratic term ($g_{ck} < 0$) in order to minimize the effect of $\delta\omega$. The dynamical states of the points labeled in Fig. 7(b) are shown in Fig. 8. It can be clearly seen that these dynamics agree well with the collective behaviors

displayed in Fig. 7(b). Furthermore, the dynamical state carried out in the gray area intertwined in the magenta zone in Fig. 7(b) is revealed. It results that the mechanical resonators exhibit Rabi oscillations within this regime, which means that they are strongly coupled and can exchange energy.

[1] M. Eichenfield, J. Chan, R. M. Camacho, K. J. Vahala, and O. Painter, Nature (London) **462**, 78 (2009).

[2] M. Eichenfield, R. Camacho, J. Chan, K. J. Vahala, and O. Painter, Nature (London) **459**, 550 (2009).

[3] D. E. Chang, A. H. Safavi-Naeini, M. Hafezi, and O. Painter, New J. Phys. **13**, 023003 (2011).

[4] G. Heinrich, M. Ludwig, J. Qian, B. Kubala, and F. Marquardt, Phys. Rev. Lett. **107**, 043603 (2011).

[5] C. B. Wallin, R. De Alba, D. Westly, G. Holland, S. Grutzik, R. H. Rand, A. T. Zehnder, V. A. Aksyuk, S. Krylov, and B. R. Ilic, Phys. Rev. Lett. **121**, 264301 (2018).

[6] M. Ludwig and F. Marquardt, Phys. Rev. Lett. **111**, 073603 (2013).

[7] A. Xuereb, C. Genes, and A. Dantan, Phys. Rev. Lett. **109**, 223601 (2012).

[8] A. Xuereb, C. Genes, G. Pupillo, M. Paternostro, and A. Dantan, Phys. Rev. Lett. **112**, 133604 (2014).

[9] W. Chen and A. A. Clerk, Phys. Rev. A **89**, 033854 (2014).

[10] M. Schmidt, V. Peano, and F. Marquardt, New J. Phys. **17**, 023025 (2015).

[11] T. F. Roque, V. Peano, O. M. Yevtushenko, and F. Marquardt, New J. Phys. **19**, 013006 (2017).

[12] V. Peano, C. Brendel, M. Schmidt, and F. Marquardt, Phys. Rev. X **5**, 031011 (2015).

[13] C. Brendel, V. Peano, O. Painter, and F. Marquardt, Phys. Rev. B **97**, 020102(R) (2018).

[14] S. Bregni, Synchronization of Digital Telecommunications Networks (Wiley, Chichester, England, 2002).

[15] S. Strogatz, *Sync: The Emerging Science of Spontaneous Order* (Hyperion, New York, 2003).

[16] T. Bahder, *Clock Synchronization and Navigation in the Vicinity of the Earth* (Nova Science, New York, 2009).

[17] M. Zhang, G. S. Wiederhecker, S. Manipatruni, A. Barnard, P. McEuen, and M. Lipson, Phys. Rev. Lett. **109**, 233906 (2012).

[18] M. Zhang, S. Shah, J. Cardenas, and M. Lipson, Phys. Rev. Lett. **115**, 163902 (2015).

[19] M. Bagheri, M. Poot, L. Fan, F. Marquardt, and H. X. Tang, Phys. Rev. Lett. **111**, 213902 (2013).

[20] M. F. Colombano, G. Arregui, N. E. Capuj, A. Pitanti, J. Maire, A. Griol, B. Garrido, A. Martinez, C. M. Sotomayor-Torres, and D. Navarro-Urrios, Phys. Rev. Lett. **123**, 017402 (2019).

[21] J. Sheng, X. Wei, C. Yang, and H. Wu, Phys. Rev. Lett. **124**, 053604 (2020).

[22] S. Y. Shah, M. Zhang, R. Rand, and M. Lipson, Phys. Rev. Lett. **114**, 113602 (2015).

[23] E. Gil-Santos, M. Labousse, C. Baker, A. Goetschy, W. Hease, C. Gomez, A. Lemaître, G. Leo, C. Ciuti, and I. Favero, Phys. Rev. Lett. **118**, 063605 (2017).

[24] H. Fu, Z.-c. Gong, L.-p. Yang, T.-h. Mao, C.-p. Sun, S. Yi, Y. Li, and G.-y. Cao, Phys. Rev. Appl. **9**, 054024 (2018).

[25] S. Hong, R. Riedinger, I. Marinković, A. Wallucks, S. G. Hofer, R. A. Norte, M. Aspelmeyer, and S. Gröblacher, Science **358**, 203 (2017).

[26] P. Djorwe, Y. Pennec, and B. Djafari-Rouhani, Phys. Rev. E **98**, 032201 (2018).

[27] R. Karnatak, N. Punetha, A. Prasad, and R. Ramaswamy, Phys. Rev. E **82**, 046219 (2010).

[28] A. Sharma, M. D. Shrimali, and S. K. Dana, Chaos: Interdiscip. J. Nonlinear Sci. **22**, 023147 (2012).

[29] P. Djorwe, Y. Pennec, and B. Djafari-Rouhani, Phys. Rev. Appl. **12**, 024002 (2019).

[30] T. K. Paraïso, M. Kalaee, L. Zang, H. Pfeifer, F. Marquardt, and O. Painter, Phys. Rev. X **5**, 041024 (2015).