# Supplementary Information:
Optomechanical Synchronization across Multi-Octave Frequency Spans

Caique C. Rodrigues1,2<EMAIL_ADDRESS>1Applied Physics Department, Gleb
Wataghin Physics Institute, University of Campinas, Campinas, SP, Brazil
2Photonics Research Center, University of Campinas, Campinas, SP, Brazil
3Department of Electrical Engineering, Columbia University, New York, New York
10027, USA
4Department of Applied Physics and Applied Mathematics, Columbia University,
New York, New York 10027, USA  Cauê M. Kersul1,2 1Applied Physics Department,
Gleb Wataghin Physics Institute, University of Campinas, Campinas, SP, Brazil
2Photonics Research Center, University of Campinas, Campinas, SP, Brazil
3Department of Electrical Engineering, Columbia University, New York, New York
10027, USA
4Department of Applied Physics and Applied Mathematics, Columbia University,
New York, New York 10027, USA  André G. Primo1,2 1Applied Physics Department,
Gleb Wataghin Physics Institute, University of Campinas, Campinas, SP, Brazil
2Photonics Research Center, University of Campinas, Campinas, SP, Brazil
3Department of Electrical Engineering, Columbia University, New York, New York
10027, USA
4Department of Applied Physics and Applied Mathematics, Columbia University,
New York, New York 10027, USA  Michal Lipson3,4 1Applied Physics Department,
Gleb Wataghin Physics Institute, University of Campinas, Campinas, SP, Brazil
2Photonics Research Center, University of Campinas, Campinas, SP, Brazil
3Department of Electrical Engineering, Columbia University, New York, New York
10027, USA
4Department of Applied Physics and Applied Mathematics, Columbia University,
New York, New York 10027, USA  Thiago P. M. Alegre1,2 1Applied Physics
Department, Gleb Wataghin Physics Institute, University of Campinas, Campinas,
SP, Brazil
2Photonics Research Center, University of Campinas, Campinas, SP, Brazil
3Department of Electrical Engineering, Columbia University, New York, New York
10027, USA
4Department of Applied Physics and Applied Mathematics, Columbia University,
New York, New York 10027, USA  Gustavo S. Wiederhecker1,2<EMAIL_ADDRESS>1Applied Physics Department, Gleb Wataghin Physics Institute, University of
Campinas, Campinas, SP, Brazil
2Photonics Research Center, University of Campinas, Campinas, SP, Brazil
3Department of Electrical Engineering, Columbia University, New York, New York
10027, USA
4Department of Applied Physics and Applied Mathematics, Columbia University,
New York, New York 10027, USA

††preprint: APS/123-QED

## Supplementary Note 1

Cavity Characterization. The whole experimental setup and the cavity geometry
are shown in Supplementary Figure 1.

Supplementary Figure 1: a) Experimental setup used in the article. A tunable
laser (TL) goes into a beam splitter (BS), in which one of the arms goes to a
HCN cell wavelength reference, and the other arm goes to a electro-optical
modulator (EOM) that is controlled by a radio frequency generator (RFG,
Agilent PSG E8251). The modulated field after interacting with the sample
inside a vacuum chamber of $\approx$ 0.1 mbar goes to another beam splitter
which we finally obtain our results. The output signal then reaches a fast
photodetector (FPD), which measures both the temporal trace using an
oscilloscope (OSC, DSO9254A) and also the spectral content at a electrical
spectrum analyzer (ESA, Keysight N9030), but also a slow photodetector (PD)
which gives the Lorentzian shape optical transmission. The final part of the
setup is a feedback loop (LaseLock) that goes back into the tunable laser that
makes the laser wavelength stable by self referenciation, avoiding unwanted
drifts during the data acquisition; b) Illustration of the nitride double disk
cavity geometry used in the experiment.

The optical and the mechanical modes used in this experiment are shown in
Supplementary Figure 2, with their best fits in red. The model of these curves
are given by well-known Equation 1 and Equation 2

$T(\Delta)=\left|\frac{s_{\text{out}}}{s_{\text{in}}}\right|^{2}=\frac{\left(1-2\eta\right)^{2}+\frac{4\Delta^{2}}{\kappa^{2}}}{1+\frac{4\Delta^{2}}{\kappa^{2}}}\quad\quad\quad\quad\text{(Optical
Transmission Spectrum)}$ (1)

$\mathcal{S}_{PP}[\Omega]=\mathcal{S}^{\text{min}}_{PP}+\frac{\left(\mathcal{S}^{\text{max}}_{PP}-\mathcal{S}^{\text{min}}_{PP}\right)\left(\Gamma_{m}\Omega_{m}\right)^{2}}{\left(\Omega^{2}-\Omega_{m}^{2}\right)^{2}+\left(\Gamma_{m}\Omega\right)^{2}}\quad\quad\quad\quad\text{(Power
Spectral Density)}$ (2)

Supplementary Figure 2: a) Experimental optical transmission spectrum of the
cavity; b) Experimental power spectrum density (PSD). The best fits of both
curves are shown in red.

The value measured for the vacuum optomechanical coupling rate was
$g_{0}/2\pi=16.2$ kHz, where we have followed the M. L. Gorodetksy et al.
article [3]. The function $s_{\text{in}}^{2}$ can be interpreted as the power
reaching the cavity, i.e., $s_{\text{in}}^{2}=P_{\text{in}}$, which is also
valid for the output field $s_{\text{out}}^{2}=P_{\text{out}}$. The power
spectral density $\mathcal{S}_{PP}$ (or just PSD) shown in Equation 2 is in
dBm units. The parameter $\eta=\kappa_{e}/\kappa$ is the coupling between the
optical fiber taper and the cavity. For now on we are omitting the sub-index
of $s_{\text{in}}(t)\rightarrow s(t)$, because we are not using
$s_{\text{out}}$ in any future calculations, so there will be no ambiguity in
just writing $s(t)$ for the input field.

## Supplementary Note 2

Fractional Synchronization. As mentioned in the article we also observed
several fractional order synchronizations, i.e., $\rho=p/q$ not an integer. We
have shown in the main article, however, only the threshold to observe the tip
of the Arnold tongues; here, we present in Supplementary Figure 3 the whole
experimental map obtained.

Supplementary Figure 3: Experimental fractional order Arnold tongues maps. The
$p:q$ order is shown in white inside each map, which are arranged from lower
frequency to higher frequencies from the left to the right.

These maps, however, requires a really strong modulation depth $\varepsilon$,
impossible to study them in the weak perturbation regime of our semi-
analytical model, as we are drastically changing the dynamic of the system.
The importance of these data is to prove the existence of this kind of
injection locking in optomechanics and to motivate the study of such
phenomenon in future works, maybe the possibility of achieving such regimes
using weak perturbations with different experimental parameters or even in
some other cavity design to enhance these effects.

## Supplementary Note 3

Phase Noise Analysis and Routes to Frequency Division Optimization. In the
main text we have shown only the phase noise measurements for the $4:1$
injection, here we will show the PN density spectra for the high-harmonics, 2
: 1 to 4 : 1, as a function of the RF power, as shown in Supplementary Figure
4. For the smallest modulation depth (-19 dBm $\approx 1.5\%$) the $2:1$ PN is
flat around -70 dBc/Hz, in contrast with the $3:1$ and $4:1$ cases. For small
modulations, both the $3:1$ and $4:1$ PN spectra appear to be transitioning
from the OMO free-running spectrum to the injection-locked regime
characterized by the flat plateau. This is expected as the farther we are from
$\Omega_{0}$ because smaller the interaction, in such a way that higher
modulation depths are needed to achieve the same low PN levels.

Supplementary Figure 4: Phase noise spectral densities for the cases $2:1$, at
a), $3:1$, at b), and $4:1$, c), for various modulation depth. The black curve
is the OMO free-running PN for reference.

A compact way to visualize the evolution of the phase noise as a function of
the RF power is averaging the phase noise around low frequencies, as done in
Supplementary Figure 5.

Supplementary Figure 5: Average phase noise around 100 Hz and 1 kHz for high
harmonic injections as function of the RF power.

We address possible paths to optimize the injection locked phase noise. The
model used for superharmonic injection phase noise is

$\mathcal{L}_{\text{out}}(\Omega)=\frac{(\Delta\Omega_{n}/n)^{2}\mathcal{L}_{\text{inj}}(\Omega)\cos^{2}\theta+\Omega^{2}\mathcal{L}_{\text{free}}(\Omega)}{\Delta\Omega_{n}^{2}\cos^{2}\theta+\Omega^{2}},$
(3)

which is a weighted average of $\mathcal{L}_{\text{inj}}(\Omega)$ and
$\mathcal{L}_{\text{free}}(\Omega)$. If any one of these is much smaller than
the other, i.e.,
$\mathcal{L}_{\text{inj}}(\Omega)\ll\mathcal{L}_{\text{free}}(\Omega)$, we
cannot expect
$\mathcal{L}_{\text{out}}(\Omega)\approx\mathcal{L}_{\text{inj}}(\Omega)$,
which would be desirable to a frequency divider. That being said, we must find
ways to improve $\mathcal{L}_{\text{free}}(\Omega)$ such that
$\mathcal{L}_{\text{free}}(\Omega)$ approaches
$\mathcal{L}_{\text{inj}}(\Omega)$. According to [2] and [5], an
optomechanical cavity in the unresolved sideband regime is dominated by
thermomechanical noise and its phase noise spectral density is given by

$\mathcal{L}_{\text{free}}(\Omega)=\left(\frac{2\Gamma_{m}}{n_{x}}\right)\left(\bar{n}_{\text{th}}+\frac{1}{2}\right)\left(\frac{1}{\Omega^{2}}+\frac{\nu_{\text{om}}^{\prime
2}}{\Omega^{2}}\frac{1}{\gamma_{\text{om}}^{\prime
2}+\Omega^{2}}+\frac{\eta_{I}^{2}}{\gamma_{\text{om}}^{\prime
2}+\Omega^{2}}\right),$ (4)

in which three “general rules” towards phase noise reduction are noticed:
reducing $\Gamma_{m}$ (mechanical oscillator’s linewidth); reducing
$\bar{n}_{\text{th}}$ (thermal phonon number); and increasing $n_{x}$
(coherent phonon number).

Quality factors up to $Q_{m}\approx 10000$, 10 times larger than in our
experiment, are reported in silicon nitride double disks [9]. Such large
mechanical quality factors were obtained just by manipulating the thickness of
the nitride film. Also, increasing the optical pump power $P_{0}$ should
readily increase the coherent phonon occupation $n_{x}$. To obtain an estimate
on a feasible amplitude enhancement, we simulate the oscillation amplitude $x$
(in units of its zero point fluctuation, $x_{\text{zpf}}$) as a function of
the optical detuning $\Delta_{x}$, where the subindex indicates that we
already counted the static optomechanical shift $x_{0}$, i.e.,
$\Delta_{x}=\omega_{l}-\omega_{0}-Gx_{0}$. Our simulations were performed for
a set of optical powers, as shown in Supplementary Figure 6, to find out the
precise scaling of $n_{x}$ with $P_{0}$. The numerical simulation details and
the parameters used are covered in the next section.

Supplementary Figure 6: a) Mechanical oscillator’s amplitude, in units of
$x_{\text{zpf}}$, as a function of the optical detuning $\Delta_{x}$, in units
of $\kappa$; b) Maximum oscillator’s amplitude for a given optical power
$P_{0}$, in mW. The floating numbers indicates how many times we must multiply
the actual pump power used in the experiment to reach that level. The
$\times$1 case is the actual numerical simulation using the data obtained in
the laboratory, which is $P_{0}=480\mu$W.

The terms $\nu_{\text{om}}^{\prime}$, $\gamma_{\text{om}}^{\prime}$ and
$\eta_{I}$ are the frequency shift due to amplitude change at the limit cycle,
the damping rate of the amplitude fluctuation and transfers from displacement
amplitude noise to photon number phase noise, respectively [2]. As Equation 4
show us, all these terms also contribute to the phase noise level, however, as
we are using an optomechanical cavity in an unresolved sideband regime, we are
neglecting their contributions because the term $1/\Omega^{2}$ dominates.
Therefore, combining higher quality factors with larger oscillation
amplitudes, we could achieve a net 30 dBc/Hz improvement, paving the way
towards future improvements in optomechanical frequency dividers.

## Supplementary Note 4

Numerical Simulation. The numerical simulations of this section are not
straightforward to perform, and it is worth discussing how they were carried
out carefully. One issue faced when solving the following coupled nonlinear
ODE

$\dot{a}=i\Delta(t)a-\frac{\kappa}{2}a-iGxa+\sqrt{\kappa_{e}}s_{0}\sqrt{1+\varepsilon(t)\sin{\Theta_{d}(t)}}\quad\quad\text{and}\quad\quad\ddot{x}+\Gamma_{m}\dot{x}+\Omega_{m}^{2}x=-\frac{\hslash
G}{m_{\text{eff}}}|a|^{2}$ (5)

is the stiff nature of the system, characterized by the need of very small
discretization steps despite the relative smoothness of the solutions. To
tackle this system, we used well known numerical packages
DifferentialEquations.jl, FFTW.jl, Sundials.jl and DSP.jl available in Julia
language which implements robust methods for such systems. The simulation was
done as follows: we first set an optical detuning function $\Delta(t)$ to
sweep linearly from $\Delta_{i}$ to $\Delta_{f}$, where the sub-indexes $i$
and $f$ means initial and final, respectively. We have chosen $\Delta_{i}>0$
because we want to access the blue side of the optical mode, where the self-
sustained dynamic is naturally accessible. After reaching $\Delta_{f}$, we
wait a few cycles of the mechanical oscillator to make sure the system is in a
stationary regime and then turned on the modulation depth $\varepsilon(t)$, in
which we modeled as a Heaviside step function. With the modulation depth
online we, once again, waited a few microseconds to stabilize the energy
inside the cavity, and then finally turned on the RF frequency sweep. In the
laboratory our RF frequency sweep was linear between $\Omega_{d}^{i}$ and
$\Omega_{d}^{f}$ with constant velocity $d\Omega_{d}/dt=\dot{\Omega}_{d}$, so
we modeled $\Theta_{d}(t)$ as a parabola, i.e.,
$d\Theta_{d}(t)/dt=\Omega_{d}(t)=\Omega_{d}^{i}+\dot{\Omega}_{d}t$. The value
chosen for $\dot{\Omega}_{d}$ need to be small to guarantee adiabaticity,
which clearly is the case in the laboratory. A good threshold for adiabaticity
is to sweep the RF tone over the mechanical resonance (of linewidth
$\Gamma_{m}$) within the mechanical lifetime, $\tau_{m}\approx
2\pi/\Gamma_{m}$, i.e.,
$\dot{\Omega}_{d}\approx\Gamma_{m}/\tau_{m}\approx\Gamma_{m}^{2}/2\pi$. For
our purposes, a RF frequency sweep velocity of $\dot{\Omega}_{d}\approx
0.1\Gamma_{m}^{2}$ was enough to ensure adiabaticity. A summary of all said is
shown in Supplementary Figure 7, highlighting the main aspects of the dynamic.

Supplementary Figure 7: a) Complete time domain simulation, showing important
aspects of the synchronization. The purple region is the transient region
which the mechanical oscillator gains amplitude. The two vertical dashed black
lines shows where exactly we turned on the modulation $\varepsilon$ and the RF
sweep $\dot{\Omega}_{d}$. The pink region is where injection locking is
happening; b) Transient region of Supplementary Figure 7(a) showing the regime
from non-oscillating cavity to self-sustained oscillation; c) Phase space of
Supplementary Figure 7(b); d) Temporal trace of the black part of
Supplementary Figure 7(b).

These are the raw data that we obtain from the simulation. To obtain from
these data the Arnold tongues we can take the length of the synchronized
region of Supplementary Figure 7(a) - the pink region of the plot - for each
modulation depth $\varepsilon$. However, we must clarify how we find this pink
region, i.e., the specific point where we say that synchronization occurs is a
bit blur in the time domain, and that’s why we construct a spectrogram, which
is the Fourier transform of our signal in function of time, as shown in
Supplementary Figure 8(a). However, instead of plotting the spectrogram as a
function of time, as we know the value of the driving RF frequency
$\Omega_{d}$ for each time $t$, we can plot the spectrogram already in
function of the RF frequency, and that is what was done in Supplementary
Figure 8(a). One way to obtain the synchronized region is to take the
horizontal slice of this spectrogram just above the mechanical oscillation
frequency $\Omega_{0}/2\pi$, which is the horizontal dashed red line, and its
plot is shown in Supplementary Figure 8(b). A second way, which is more well
known in the literature, is to plot the difference between the driving
frequency and oscillator’s frequency ($\Omega_{d}-\Omega_{0}$) as a function
of the drive frequency itself (or, in our case, the driving frequency minus a
constant natural mechanical frequency $\Omega_{m}$), as shown in Supplementary
Figure 8(c).

Supplementary Figure 8: a) Spectrogram of the transmission signal
(Supplementary Figure 7(a)) after the RF sweep is turned on. The dashed black
line is the value of the drive frequency. The vertical bold black line is the
value of the mechanical mode frequency; b) Horizontal dashed red slice of the
spectrogram shown in a); c) Typical plot of synchronization systems showing
the mismatch from driving frequency and bare oscillation frequency, making
clear where these became the same, defining a synchronized state. We have used
$\varepsilon=2\%$ for these simulations.

The Arnold tongues constructed using the explanation above are shown in
Supplementary Figure 9, which was already presented at the article as Fig.
3(a). As we can see, the simulation shows bigger synchronized region for the
case $p:q=2:1$ than $1:1$, and also a pretty wide $4:1$ AT, but a small $3:1$,
the same trend of the experimental data. Supplementary Table 1 shows the
parameters used in the simulations done and Supplementary Table 1 shows the
conversion from RF power, in dBm, to modulation depth $\varepsilon$, in $\%$,
which is based in experimental data. The actual formula for the RF power
$P_{RF}$ is shown inside the plot.

Supplementary Figure 9: Simulated Arnold tongues using injection frequency
$\Omega_{d}=p\Omega_{0}/q$ for the cases $p=\\{1,2,3,4\\}$ and $q=1$, in
order, from a) to d). To simulate these maps we used Equation 5, being these
the ones that we will always call “full model”. The dashed black line at a) is
the region that we analyzed in Supplementary Figure 8(c).

Parameters | Values
---|---
$P_{0}$ | 425 $\mu$W
$\lambda$ | 1560nm
$\Delta_{i}$ | 8$\kappa$
$\Delta_{f}$ | 0.35$\kappa$
$d\Delta/dt$ | 10${}^{2}\Gamma_{m}^{2}$
$d\Omega_{d}/dt$ | 0.075$\Gamma_{m}^{2}$
$\eta$ | 0.75
$\kappa/2\pi$ | 1.16 GHz
$Q_{\text{opt}}$ | 165000
$\Omega_{m}/2\pi$ | 31.86 MHz
$\Gamma_{m}/2\pi$ | 25.37 kHz
$Q_{m}$ | 1255
$g_{0}/2\pi$ | 16.2 kHz
Supplementary Table 1: Parameters values used in every simulations, unless
explicitly mentioned the opposite.

Supplementary Figure 10: Conversion from modulation depth, in percentage, to
RF power, in dBm. The inset shows the actual formula, in S.I. units, of this
graph. $V_{\pi}=$ 5.5 V and $R=$ 50 $\Omega$.

## Supplementary Note 5

Semi-Analytical Model. The Hamiltonian of our system, regardless dissipative
considerations, can be modeled by

$H=\hslash\left(\omega_{0}+Gx\right)a^{\dagger}a+\frac{p^{2}}{2m_{\text{eff}}}+\frac{m_{\text{eff}}\Omega_{m}x^{2}}{2}+i\hslash\sqrt{\kappa_{e}}s_{0}\left(e^{-i\omega_{l}t}a^{\dagger}-e^{i\omega_{l}t}a\right),$
(6)

in which $\hslash$ is the reduced Planck’s constant, $\omega_{0}$ is the
unperturbed angular frequency of the optical mode, $a^{\dagger}$ and $a$ are
the creation and annihilation operators for photons with energy
$\hslash\omega_{0}$, respectively, $G$ is the first order coefficient of the
Taylor expansion of $\omega(x)=\omega_{0}+Gx$ evaluated at mechanical
equilibrium position $\left(\text{i.e., }G=(d\omega/dx)|_{x=0}\right)$, $p$
and $x$ are the momentum and position operator of the mechanical oscillator,
respectively, $\Omega_{m}$ is the unperturbed angular frequency of the
mechanical mode, $m_{\text{eff}}$ is the mechanical oscillator effective mass,
$i$ is the complex unity, $\kappa_{e}$ is the external optical coupling rate,
$s_{0}^{2}$ is the input power and $\omega_{l}$ is the optical pump angular
frequency. As long as we are not interested in quantum phenomena we can study
our dynamical system just looking to the average value of these operators, and
we can also introduce the optical and mechanical loses $\kappa$ and
$\Gamma_{m}$ directly in the equations of motion [1] as

$\dot{a}=i\Delta
a-\frac{\kappa}{2}a-iGxa+\sqrt{\kappa_{e}}s_{0}\quad\quad\text{and}\quad\quad\ddot{x}+\Gamma_{m}\dot{x}+\Omega_{m}^{2}x=-\frac{\hslash
G}{m_{\text{eff}}}|a|^{2},$ (7)

where we already changed $a$ to the slow rotating frame of reference
$a\rightarrow ae^{-i\omega_{l}t}$ to let the equation autonomous. We define
the bare optical detuning $\Delta=\omega_{l}-\omega_{0}$ as the difference
between optical pump frequency and unperturbed optical mode frequency. To
introduce the amplitude modulation used in the experiment we can simply
multiply $s_{0}$ by a factor $\sqrt{1+\varepsilon\sin{\Theta_{d}(t)}}$ in
Equation 7, i.e., $s(t)=s_{0}\sqrt{1+\varepsilon\sin{\Theta_{d}(t)}}$, in
which $\varepsilon$ is the modulation depth and is correlated with the RF
power as shown in Supplementary Table 1. The term $\Theta_{d}(t)$ is the phase
of this modulation, in which most of the time will just be $\Omega_{d}t$.
Simulating Equation 7 as it is shown require us to know $G$ and
$m_{\text{eff}}$ but, because they are normalization dependent, we are going
to avoid this using $g_{0}=Gx_{\text{zpf}}$, the optomechanical single-photon
coupling strength, and also the
$x_{\text{zpf}}=\sqrt{\hslash/2m_{\text{eff}}\Omega_{m}}$, the zero points
fluctuation amplitude of $x$. New normalizations will be used to study self-
sustained oscillations, which are given by

$x(t)=x_{0}+\delta
x(t)=\left(\frac{\kappa}{2g_{0}}\right)\sqrt{\frac{\hbar}{2m_{\text{eff}}\Omega_{m}}}\left(\widetilde{x}_{0}+\widetilde{x}(t)\right)\quad\quad\text{and}\quad\quad
t=\frac{\widetilde{t}}{\Omega_{m}},$ (8)

where all tilde variables are now adimensional. The terms $x_{0}$ and $\delta
x(t)$ are the DC and AC components of $x(t)$, respectively, being
$\widetilde{x}_{0}$ and $\widetilde{x}(t)$ their adimensional version. We can
then rewrite Equation 7 as

$\frac{da}{d\widetilde{t}}=-\frac{\kappa}{2\Omega_{m}}a+i\left(\frac{\Delta_{x}}{\Omega_{m}}-\frac{\kappa}{2\Omega_{m}}\widetilde{x}\right)a+\sqrt{\frac{\kappa_{e}}{\Omega_{m}}}\frac{s}{\sqrt{\Omega_{m}}}\quad\quad\text{and}\quad\quad\frac{d^{2}\widetilde{x}}{d\widetilde{t}^{2}}+\frac{1}{Q_{m}}\frac{d\widetilde{x}}{d\widetilde{t}}+\widetilde{x}=-\widetilde{x}_{0}-\frac{\mathcal{C}_{0}}{Q_{m}}|a|^{2},$
(9)

where the new variables $\Delta_{x}$, $Q_{m}$ and $\mathcal{C}_{0}$ are the
optical detuning with the static optomechanical shift correction,
$\Delta_{x}=\omega_{l}-\omega_{0}-Gx_{0}$, the mechanical oscillator quality
factor, $Q_{m}=\Omega_{m}/\Gamma_{m}$, and the single-photon cooperativity
$\mathcal{C}_{0}=4g_{0}^{2}/\Gamma_{m}\kappa$, respectively. The derivatives
now will be all with respect to $\widetilde{t}$ unless explicitly said the
opposite. The motivations to construct a semi-analytical model in this article
are (i) to prove that each term $F_{n}$ from the power expansion of the
optical force is mainly responsible for the $p:q$ Arnold tongue width
$\Delta\Omega(p,q)$, (ii) to obtain a semi-analytical formula for these
$\Delta\Omega(p,q)$, (iii) to show that the influence of the optical detuning
greatly change the synchronization region, (iv) to prove that the symmetry
breaking term $F_{2}$, which is neglected in many articles, is actually
crucial for the dynamic and (v) to explain/predict those sidebands around the
synchronization region. We then start uncoupling Equation 9 using adiabatic
considerations: our optomechanical cavity has mechanical linewidth
$\Gamma_{m}$ much smaller than the optical linewidth $\kappa$, as well the
mechanical frequency $\Omega_{m}$ also much smaller than $\kappa$, a regime
called unresolved sidebands. We can then assume that $a(\widetilde{t})$ is
always in equilibrium with $\widetilde{x}(\widetilde{t}-\widetilde{\tau})$,
where $\widetilde{\tau}$ is some adimensional time delay that we will deduce
later in Equation 23, so we can write $a(\widetilde{t})$ as

$a(\widetilde{t})\approx\sqrt{\frac{\kappa_{e}}{\kappa}}\frac{s(\widetilde{t})}{\sqrt{\kappa}}\frac{2}{1-i\left[\frac{2\Delta_{x}}{\kappa}-\widetilde{x}(\widetilde{t}-\widetilde{\tau})\right]},$
(10)

and then we can analyze the whole system just looking to one equation

$\ddot{\widetilde{x}}(\widetilde{t})+\frac{\dot{\widetilde{x}}(\widetilde{t})}{Q_{m}}+\widetilde{x}(\widetilde{t})=-\widetilde{x}_{0}-\frac{\mathcal{C}_{0}}{Q_{m}}\left(\frac{2\kappa_{e}}{\kappa}\right)\left(\frac{2s_{0}^{2}}{\kappa}\right)\left(\frac{1+\varepsilon\sin{\Theta_{d}(\widetilde{t})}}{1+\left[\frac{2\Delta_{x}}{\kappa}-\widetilde{x}(\widetilde{t}-\widetilde{\tau})\right]^{2}}\right),$
(11)

where we define $f(\widetilde{t})$ and $f_{0}$ as

$f(\widetilde{t})=f_{0}\left[1+\varepsilon\sin{\Theta_{d}(\widetilde{t})}\right]\quad\text{and}\quad
f_{0}=\frac{\mathcal{C}_{0}}{Q_{m}}\left(\frac{2\kappa_{e}}{\kappa}\right)\left(\frac{2s_{0}^{2}}{\kappa}\right),$
(12)

which allows us to rewrite Equation 10 as

$\ddot{\widetilde{x}}(\widetilde{t})+\frac{\dot{\widetilde{x}}(\widetilde{t})}{Q_{m}}+\widetilde{x}(\widetilde{t})=-\widetilde{x}_{0}-\frac{f(\widetilde{t})}{1+\left[\frac{2\Delta_{x}}{\kappa}-\widetilde{x}(\widetilde{t}-\widetilde{\tau})\right]^{2}},$
(13)

However, Equation 13 is still very complicated because it is a non-autonomous
delay differential equation, so we will expand the RHS in a power series of
$\widetilde{x}(\widetilde{t}-\widetilde{\tau})$ as

$\frac{1}{1+\left[\frac{2\Delta_{x}}{\kappa}-\widetilde{x}(\widetilde{t}-\widetilde{\tau})\right]^{2}}=F_{0}+F_{1}\widetilde{x}(\widetilde{t}-\widetilde{\tau})+F_{2}\widetilde{x}^{2}(\widetilde{t}-\widetilde{\tau})+F_{3}\widetilde{x}^{3}(\widetilde{t}-\widetilde{\tau})+...$
(14)

in which the actual form of these first coefficients (which were already shown
in Fig. 1(b) of the article) are

$F_{0}=\frac{1}{1+\frac{4\Delta_{x}^{2}}{\kappa^{2}}}\quad,\quad
F_{1}=\frac{2\left(\frac{2\Delta_{x}}{\kappa}\right)}{\left(1+\frac{4\Delta_{x}^{2}}{\kappa^{2}}\right)^{2}}\quad,\quad
F_{2}=\frac{\left(\frac{12\Delta_{x}^{2}}{\kappa^{2}}-1\right)}{\left(1+\frac{4\Delta_{x}^{2}}{\kappa^{2}}\right)^{3}}\quad\text{and}\quad
F_{3}=\frac{4\left(\frac{2\Delta_{x}}{\kappa}\right)\left(\frac{4\Delta_{x}^{2}}{\kappa^{2}}-1\right)}{\left(1+\frac{4\Delta_{x}^{2}}{\kappa^{2}}\right)^{4}},$
(15)

which shows that large normalized detuning ($\Delta_{x}/\kappa\gg 1$) leads to
negligible values, as each $F_{n+1}$ term decreases faster than $F_{n}$ as a
function of $\Delta_{x}/\kappa$, i.e.,

$\frac{F_{n}}{F_{n+1}}\sim\left(\frac{\Delta_{x}}{\kappa}\right).$ (16)

The value of each $F_{n}$ in our experiment was found to be $F_{0}=0.6711$,
$F_{1}=0.6306$, $F_{2}=0.1421$ and $F_{3}=-0.2897$. Substituting Equation 14
in Equation 13 reveals the nonlinear nature of the optical feedback into the
mechanical oscillator,

$\ddot{\widetilde{x}}(\widetilde{t})+\frac{\dot{\widetilde{x}}(\widetilde{t})}{Q_{m}}+\widetilde{x}(\widetilde{t})=-\widetilde{x}_{0}-f(\widetilde{t})\left[F_{0}+F_{1}\widetilde{x}(\widetilde{t}-\widetilde{\tau})+F_{2}\widetilde{x}^{2}(\widetilde{t}-\widetilde{\tau})+F_{3}\widetilde{x}^{3}(\widetilde{t}-\widetilde{\tau})\right].$
(17)

To remove the delay dependence we can expand
$\widetilde{x}(\widetilde{t}-\widetilde{\tau})$ in powers of
$\widetilde{\tau}$ as

$\widetilde{x}^{n}(\widetilde{t}-\widetilde{\tau})=\widetilde{x}^{n}(\widetilde{t})-n\widetilde{\tau}\widetilde{x}^{n-1}(\widetilde{t})\dot{\widetilde{x}}(\widetilde{t})+O(\widetilde{\tau}^{2}),$
(18)

where we are neglecting $O(\widetilde{\tau}^{2})$ since we have that
$\widetilde{\tau}^{2}$ is of order $O(\Omega_{m}^{2}/\kappa^{2})$, as we will
verify soon. We can then group all these terms in an arrangement that
highlights how far from the ideal harmonic oscillator this system is, as shown
in Equation 19

$\ddot{\widetilde{x}}+\left[\frac{1}{Q_{m}}-\widetilde{\tau}f(\widetilde{t})\left(F_{1}+2F_{2}\widetilde{x}+3F_{3}\widetilde{x}^{2}\right)\right]\dot{\widetilde{x}}+\left[1+f(\widetilde{t})\left(F_{1}+F_{2}\widetilde{x}+F_{3}\widetilde{x}^{2}\right)\right]\widetilde{x}=-f_{1}(\widetilde{t})F_{0}-\left(f_{0}F_{0}+\widetilde{x}_{0}\right),$
(19)

where $f_{1}(\widetilde{t})$ is the AC component of $f(\widetilde{t})$, i.e.,
$f(\widetilde{t})=f_{0}+f_{1}(\widetilde{t})$ with
$f_{1}(\widetilde{t})=\varepsilon f_{0}\sin{\Theta_{d}(\widetilde{t})}$. The
last term of Equation 19 in parentheses must be zero because is the only DC
component of the whole equation and, as one can verify, solving
$f_{0}F_{0}+\widetilde{x}_{0}=0$ returns the same static correction for the
mechanical displacement $x(t)$ from the linearized optomechanics. Equation 19
is, finally, a well looking shape equation, very similar with Equation 8 from
reference [7] where their $\zeta$ and $f_{\text{e}}$ are related with our
$F_{1}$ and $F_{0}$, besides we are not using temperature dynamics here. The
generalization of Equation 19 until $O(\widetilde{\tau}^{2})$ is given by

$\ddot{\widetilde{x}}+\left[\frac{1}{Q_{m}}-\widetilde{\tau}f(\widetilde{t})\sum_{n=1}^{n_{\text{max}}}nF_{n}\widetilde{x}^{n-1}\right]\dot{\widetilde{x}}+\left[1+f(\widetilde{t})\sum_{n=1}^{n_{\text{max}}}F_{n}\widetilde{x}^{n-1}\right]\widetilde{x}=-f_{1}(\widetilde{t})F_{0},$
(20)

but we are not going to analyze this system, we will stick with the case
$n_{\text{max}}=3$. We can do a final approximation which is to neglect terms
of order $O(\widetilde{\tau}\varepsilon)$ as far $\varepsilon$ is kept small,
knowing a priori that $\widetilde{\tau}$ is already small. We then have
Equation 21

$\ddot{\widetilde{x}}+\left[\frac{1}{Q_{m}}-\widetilde{\tau}f_{0}\left(F_{1}+2F_{2}\widetilde{x}+3F_{3}\widetilde{x}^{2}\right)\right]\dot{\widetilde{x}}+\left[1+\left(f_{0}+f_{1}\right)\left(F_{1}+F_{2}\widetilde{x}+F_{3}\widetilde{x}^{2}\right)\right]\widetilde{x}=-f_{1}F_{0}$
(21)

and finally it is an ODE that all the terms proportional to $\widetilde{x}$
have a parametric excitation $f_{1}$, but terms proportional to
$\dot{\widetilde{x}}$ do not. The only formula that is missing is
$\widetilde{\tau}=\widetilde{\tau}(\Delta_{x})$ for us to start studying
Equation 21. To find this missing expression note (i) that the term
$f_{0}F_{1}$ acts like a constant shift in the frequency, so we can associate
with it the optical spring effect and (ii) that the term
$-\widetilde{\tau}f_{0}F_{1}$ acts like a constant change in the mechanical
linewidth, so we can associate it with the optical cooling/heating. Doing that
interconnection with the linearized optomechanical equations [1] we can
identify an analytic expression for $\widetilde{\tau}$ because from our model
we have that $\delta\Gamma_{m}^{\text{linear}}=-\widetilde{\tau}f_{0}F_{1}$
and, from the linearized optomechanical equations,

$\delta\Gamma_{m}^{\text{linear}}=\frac{\mathcal{C}_{0}}{Q_{m}}\left(\frac{2\kappa_{e}}{\kappa}\right)\left(\frac{2s_{0}^{2}}{\kappa}\right)\left(\frac{1}{1+\left(\frac{2\Delta_{x}}{\kappa}\right)^{2}}\right)\left(\frac{1}{1+\left(\frac{2\Delta_{x}}{\kappa}+\frac{2\Omega_{m}}{\kappa}\right)^{2}}-\frac{1}{1+\left(\frac{2\Delta_{x}}{\kappa}-\frac{2\Omega_{m}}{\kappa}\right)^{2}}\right),$
(22)

remembering that $\delta\Gamma_{m}^{\text{linear}}$ is in units of
$\Omega_{m}$. Solving then for $\widetilde{\tau}$

$\widetilde{\tau}=\frac{1}{2\left(\frac{2\Delta_{x}}{\kappa}\right)}\left[\frac{1+\left(\frac{2\Delta_{x}}{\kappa}\right)^{2}}{1+\left(\frac{2\Delta_{x}}{\kappa}-\frac{2\Omega_{m}}{\kappa}\right)^{2}}-\frac{1+\left(\frac{2\Delta_{x}}{\kappa}\right)^{2}}{1+\left(\frac{2\Delta_{x}}{\kappa}+\frac{2\Omega_{m}}{\kappa}\right)^{2}}\right].$
(23)

We must emphasize here that our $\widetilde{\tau}$ should, rigorously, be
rewritten as $\widetilde{\tau}_{\text{linear}}$, because that is just the
first order correction of $\widetilde{\tau}$. Nevertheless, the function
$\widetilde{\tau}=\widetilde{\tau}(\Delta_{x})$ has every property that we
expect: (i) is positive for every value of $\Delta_{x}$, (ii) it is
approximately $\Omega_{m}/\kappa$ and (iii) is also consistent with the fact
that far from the resonance there is no mechanical response, i.e.,
$\widetilde{\tau}(|\Delta_{x}|\gg\kappa)=0$, as shown in Supplementary Figure
11 using our experimental parameters and also for three other cases.

Supplementary Figure 11: Adimensional linear mechanical relaxation time
$\widetilde{\tau}$ (in units of $\Omega_{m}/\kappa$) as function of the
optical detuning $\Delta_{x}$ (in units of $\kappa$).

The Arnold tongues simulations using this semi-analytical model from Equation
21 are shown in Supplementary Figure 12, and the comparison between this and
Supplementary Figure 9 is striking. Before finishing this section we want to
prove that the connection between $F_{n}$’s and the entire AT width
$\Delta\Omega(n+1,1)$ is really strong and, for that, we simulated again
Equation 21 but here we want to show not the AT map, but some particular
spectrograms while we considered only one $F_{n}$ each time we simulate the
system, as shown in Supplementary Figure 13.

Supplementary Figure 12: Simulated Arnold tongues using injection frequency
$\Omega_{d}=p\Omega_{0}/q$ for the cases $p=\\{1,2,3,4\\}$ and $q=1$, in
order, from a) to d). To simulate these maps we used Equation 21.
Supplementary Figure 13: These simulations were done considering only one
$\varepsilon F_{n}$ term in Equation 21, e.g., the first row of this image-
matrix-like consist of considering $\varepsilon F_{1}=\varepsilon
F_{2}=\varepsilon F_{3}=0$ while $\varepsilon F_{0}\neq 0$ and doing the
injection locking for all $p=1-4$. The same logic is valid for the others
rows, as indicated in the right side of the figure. In these simulations we
used $\varepsilon=5\%$.

And as we can see, simulating Equation 21 proves that the major dependence of
each AT is indeed the parametric term $\varepsilon F_{n}$, as each one of
these terms alone almost reproduces the whole dynamic of the system in a
specific region. This hierarchical dependence is explicitly calculated, as
shown in Equation 48.

## Supplementary Note 6

The averaging method of Krylov-Bogoliubov-Mitropolsky (KBM). While the
numerical simulation does predict many features and give us many insights
about the observed data, it does not provide a direct prediction of the
synchronization behavior. To pursue further analytical insight, we resort to
the KBM method to derive amplitude and phase equations describing the coupling
optomechanical oscillator enslaved by the driving modulated RF signal. We
start here introducing new adimensional time $T$ and displacement $y$ given by

$\frac{\widetilde{x}}{y}=L_{\widetilde{x}}=-\frac{F_{2}}{3F_{3}}+\sqrt{\left(\frac{F_{2}}{3F_{3}}\right)^{2}+\left(\frac{1-f_{0}F_{1}Q_{m}\widetilde{\tau}}{3f_{0}F_{3}Q_{m}\widetilde{\tau}}\right)}\quad\quad\text{and}\quad\quad\frac{\widetilde{t}}{T}=L_{\widetilde{t}}=\frac{1}{\sqrt{1+f_{0}F_{1}}}.$
(24)

At first glance, it seems like an awkward choice of normalization, but as we
discuss below, they have a clear physical interpretation. Equation 24 is the
positive root of the coefficient of $\dot{\widetilde{x}}(\widetilde{t})$ in
Equation 21, which makes the amplitude of $y$ near the value of the limit
circle of a van-der Pol oscillator $\approx O(1)$, i.e., we are renormalizing
$\widetilde{x}$ by the positive solution of

$\frac{1}{Q_{m}}-\widetilde{\tau}f_{0}\left(F_{1}+2F_{2}\widetilde{x}+3F_{3}\widetilde{x}^{2}\right)=0.$
(25)

The choice of the new time scale makes the oscillation frequency of the
oscillator, already accounted by the optical spring effect, about $\approx
O(1)$. After these normalizations we have that Equation 21 becomes

$\frac{d^{2}y}{dT^{2}}-\mu(1-y)(1+\sigma
y)\frac{dy}{dT}+\left[1+\varepsilon\alpha\sin{\left(\omega
T\right)}\right]y+\\\ +\beta\left[1+\varepsilon\sin{\left(\omega
T\right)}\right]y^{2}+\gamma\left[1+\varepsilon\sin{\left(\omega
T\right)}\right]y^{3}=F\varepsilon\sin{\left(\omega T\right)},$ (26)

with new adimensional parameters defined as

$\omega=\frac{\Omega_{d}}{\Omega_{m}}L_{\widetilde{t}},\quad\quad\quad\mu=\left(\widetilde{\tau}f_{0}F_{1}Q_{m}-1\right)L_{\widetilde{t}},\quad\quad\quad\sigma=1+\frac{2\widetilde{\tau}f_{0}F_{2}Q_{m}}{\widetilde{\tau}f_{0}F_{1}Q_{m}-1}L_{\widetilde{x}},$
(27)

$\alpha=f_{0}F_{1}L_{\widetilde{t}}^{2},\quad\quad\quad\beta=f_{0}F_{2}L_{\widetilde{x}}L_{\widetilde{t}}^{2},\quad\quad\quad\gamma=f_{0}F_{3}L^{2}_{\widetilde{x}}L_{\widetilde{t}}^{2},\quad\quad\text{and}\quad\quad
F=-f_{0}F_{0}L_{\widetilde{t}}^{2}/L_{\widetilde{x}}.$

It is evident that every parametric term $\alpha$, $\beta$, $\gamma$ and $F$
is proportional to $f_{0}$, regardless of the convoluted terms
$L_{\widetilde{t}}$ and $L_{\widetilde{x}}$, meaning that higher optical pump
intensity enhance these terms. Also, each of these terms are proportional to
one $F_{n}$, making clear distinction where each nonlinearity really is. Such
model returns us the same used by Shreyas Y. Shah at [6] if $\mu=\gamma=F=0$
and also neglecting the autonomous quadratic term $\beta y^{2}$ (which is the
term that comes from a odd power potential, making the problem parity
asymmetric), then

$\frac{d^{2}y}{dT^{2}}+\left[1+\alpha\varepsilon\sin{\left(\omega
T\right)}\right]y+\beta\varepsilon\sin{\left(\omega T\right)}y^{2}=0,$ (28)

but to leave in the exact shape of the one used there we should change the
independent variable $\omega T\rightarrow U+\frac{\pi}{2}$ and then

$\frac{d^{2}y}{dU^{2}}+\left[\frac{1}{\omega^{2}}+\frac{\alpha}{\omega^{2}}\varepsilon\cos{\left(U\right)}\right]y+\frac{\beta}{\omega^{2}}\varepsilon\cos{\left(U\right)}y^{2}=0,$
(29)

from which we interpret the parameters as

$\delta^{\text{Shah}}=\frac{1}{\omega^{2}}\quad,\quad
D_{1}^{\text{Shah}}=\frac{\alpha}{\omega^{2}}\quad,\quad
D_{2}^{\text{Shah}}=\frac{\beta}{\omega^{2}}\quad,\quad\gamma^{\text{Shah}}=\varepsilon.$
(30)

However, we will not use the multi-scale method to study synchronization
neither to find bifurcations, our analysis will be based in the KBM method of
averaging. The value of the parameters obtained from simulations are
$\mu=9.813\times 10^{-4}$, $\sigma=1.665\times 10^{0}$, $\alpha=2.383\times
10^{-2}$, $\beta=4.396\times 10^{-3}$, $\gamma=-7.340\times 10^{-3}$ and
$F=-3.098\times 10^{-2}$.

We begin our KBM analyzes from Equation 26, which is a nonlinear oscillators
of the form

$\frac{d^{2}y}{dT^{2}}+y=K\left(T,y,\frac{dy}{dT}\right),$ (31)

where $K$ is small compared to $y$. If $K(T,y,\frac{dy}{dT})=0$, we would have
the ideal harmonic oscillator with solution $y=A\sin{\left(T+\Phi\right)}$ for
any choice of constants $A$ and $\Phi$. If we now try solving Equation 31 with
slowly varying amplitude and phase $(A(T),\Phi(T))$ as ansatz, i.e.,

$y=A(T)\sin{\left[T+\Phi(T)\right]}\quad\quad\text{and}\quad\quad\frac{dy}{dT}=A(T)\cos{\left[T+\Phi(T)\right]},$
(32)

we can show [4] that this system has general solution given by Equation 33

$\begin{cases}\frac{dA}{dT}=\cos{\left(\phi\right)}K\left(T,A\sin{\phi},A\cos{\phi}\right)\\\
\\\ \phi(T)=T+\Phi(T)\\\ \\\
\frac{d\Phi}{dT}=-\frac{\sin{\left(\phi\right)}}{A}K\left(T,A\sin{\phi},A\cos{\phi}\right)\end{cases}$
(33)

The KBM method take its place here, where we average these equations over one
period, however, the integral in $T$ is replaced over a integral in $\phi$
considering that $d\phi\approx dT$, which is correct to zero order in
$\Phi(T)$, so

$\left<\frac{dA}{dT}\right>_{T}\approx\left<\frac{dA}{dT}\right>_{\phi}=\frac{1}{2\pi}\int_{0}^{2\pi}\cos{\left(\phi\right)}K\left(\phi-\Phi,A\sin{\phi},A\cos{\phi}\right)d\phi$
(34)

$\left<\frac{d\Phi}{dT}\right>_{T}\approx\left<\frac{d\Phi}{dT}\right>_{\phi}=-\frac{1}{2\pi
A}\int_{0}^{2\pi}\sin{\left(\phi\right)}K\left(\phi-\Phi,A\sin{\phi},A\cos{\phi}\right)d\phi$
(35)

and if our system were autonomous the integrals from Equation 34 and Equation
35 would be relatively easy to proceed, however, we have an external drive and
this makes our system non-autonomous. To proceed the integral, we need to deal
with $\Phi$, which we will just let constant during integration, arguing that
$\Phi$ is a slow varying function of $T$. The general form of
$K\left(T,y,\frac{dy}{dT}\right)$ for our system can be splitted in two
contribution: one autonomous and another one non-autonomous, i.e.,

$K\left(T,y,\frac{dy}{dT}\right)=K_{\text{auto}}\left(y,\frac{dy}{dT}\right)+K_{\text{non-
auto}}\left(T,y,\frac{dy}{dT}\right)$ (36)

in which each part is written as

$K_{\text{auto}}\left(y,\frac{dy}{dT}\right)=\mu\left(1-y\right)(1+\sigma
y)\frac{dy}{dT}-\beta y^{2}-\gamma y^{3}$ (37)

$\text{and}\quad\quad K_{\text{non-
auto}}\left(T,y,\frac{dy}{dT}\right)=\varepsilon\sin(\omega T)\left(F-\alpha
y-\beta y^{2}-\gamma y^{3}\right).$ (38)

Substituting $K\left(T,y,\frac{dy}{dT}\right)$ into Equation 34 and Equation
35 and then performing the integration over one period of $\phi$ we obtain our
amplitude and phase differential equations given by

$\left<\frac{dA}{dT}\right>_{\phi}=\frac{\mu A}{2}\left(1-\frac{\sigma
A^{2}}{4}\right)+\\\
+\frac{\varepsilon\sin(\pi\omega)}{\pi}\left(\frac{\omega\left[F\left(\omega^{2}-9\right)+2\beta
A^{2}\right]\sin{\left[(\pi-\Phi)\omega\right]}}{\left(\omega^{2}-9\right)\left(\omega^{2}-1\right)}-\frac{A\left[\alpha\left(\omega^{2}-16\right)-6\gamma
A^{2}\right]\cos{\left[(\pi-\Phi)\omega\right]}}{\left(\omega^{2}-16\right)\left(\omega^{2}-4\right)}\right),$
(39)

$\left<\frac{d\Phi}{dT}\right>_{\phi}=\frac{3\gamma A^{2}}{8}+\\\
-\frac{\varepsilon\sin(\pi\omega)}{\pi}\left(\frac{\left[F\left(\omega^{2}-9\right)+6\beta
A^{2}\right]\cos{\left[(\pi-\Phi)\omega\right]}}{\left(\omega^{2}-9\right)\left(\omega^{2}-1\right)A}-\frac{2\left[\alpha\left(\omega^{2}-16\right)-12A^{2}\gamma\right]\sin{\left[(\pi-\Phi)\omega\right]}}{\omega\left(\omega^{2}-16\right)\left(\omega^{2}-4\right)}\right).$
(40)

This averaging technique is the essence of the KBM method to obtain amplitude
and phase equations of nonlinear oscillators. Before we proceed, we can study
Equation 39 and Equation 40 for the case $\varepsilon=0$, which would give us
exact solutions for both $A(T)$ and $\Phi(T)$ as

$\begin{cases}A(T)=\frac{\pm
2}{\sqrt{\sigma+\left(\frac{4}{A_{0}}-\sigma\right)e^{-\mu
T}}}\quad\Rightarrow\quad\lim_{T\rightarrow\infty}A(T)=A_{\infty}=\pm\frac{2}{\sqrt{\sigma}}\\\
\\\ \Phi(T)=\Phi_{0}+\frac{3\gamma}{2\mu\sigma}\ln{\left[1+\frac{\sigma
A_{0}^{2}}{4}\left(e^{\mu
T}-1\right)\right]}\quad\Rightarrow\quad\lim_{T\rightarrow\infty}\Phi(T)=\Phi_{\infty}=\frac{3\gamma}{2\sigma}T\end{cases}$
(41)

for constants $A_{0}$ and $\Phi_{0}$. Here we conclude that even with zero
modulation depth there exist a frequency shift contribution that comes from
the Duffing term $\gamma$. The steady oscillation frequency $\Omega_{0}$ is
given by

$\Omega_{0}(\varepsilon=0)=\lim_{t\rightarrow\infty}\frac{d\phi}{dt}=\lim_{t\rightarrow\infty}\frac{d\phi}{dT}\frac{dT}{dt}=\frac{\Omega_{m}}{L_{\widetilde{t}}}\frac{d}{dT}\left(T+\Phi_{\infty}\right)=\Omega_{m}\sqrt{1+f_{0}F_{1}}\left(1+\frac{3\gamma}{2\sigma}\right)$
(42)

which can be used to estimate the Duffing term from the measured oscillation
frequency. Remembering the definitions of $\omega$ and $\rho$ we can obtain a
really nice relation between them as

$\omega=\frac{\Omega_{d}}{\Omega_{m}\sqrt{1+f_{0}F_{1}}}\quad\quad\&\quad\quad\rho=\frac{p}{q}=\frac{\Omega_{d}}{\Omega_{0}}\quad\quad\Rightarrow\quad\quad\frac{\omega}{\rho}=1+\frac{3\gamma}{2\sigma}+O(\varepsilon),$
(43)

and its clear now that if $\rho$ is an integer, it does not mean that $\omega$
is an integer. In our case (and most of the cases) the Duffing correction in
the frequency is really small, as a matter of fact we have

$\frac{3\gamma}{2\sigma}\approx-6.6125\times 10^{-3}=O(10^{-3})\ll O(1)$ (44)

so we will considerer $\omega\approx\rho$ for now on, because that is an
excellent approach to obtain simple, but good, analytical results. The phase
equation can be expanded in the vicinity of some integer $\rho=\\{1,2,3,4\\}$
to give us insight into higher harmonic synchronization. We have then the
following cases

$\begin{cases}\lim_{\rho\rightarrow
1}\left<\frac{d\Phi}{dT}\right>_{\phi}\approx\frac{3\gamma
A^{2}}{8}-\varepsilon\left(\frac{F}{2A}-\frac{3\beta
A}{8}\right)\cos{\left(\Phi\right)},\\\ \\\ \lim_{\rho\rightarrow
2}\left<\frac{d\Phi}{dT}\right>_{\phi}\approx\frac{3\gamma
A^{2}}{8}+\varepsilon\left(\frac{\alpha}{4}+\frac{\gamma
A^{2}}{4}\right)\sin{\left(2\Phi\right)}.\\\ \\\ \lim_{\rho\rightarrow
3}\left<\frac{d\Phi}{dT}\right>_{\phi}\approx\frac{3\gamma
A^{2}}{8}-\varepsilon\frac{\beta A}{8}\cos{\left(3\Phi\right)}.\\\ \\\
\lim_{\rho\rightarrow
4}\left<\frac{d\Phi}{dT}\right>_{\phi}\approx\frac{3\gamma
A^{2}}{8}-\varepsilon\frac{\gamma
A^{2}}{16}\sin{\left(4\Phi\right)}.\end{cases}$ (45)

When the system is locked to the driving signal, we know that the amplitude
$A(T)$ of the oscillator is almost constant (this is the Kuramoto
approximation) and that the phase $\phi(T)=T+\Phi(T)$ is a linear function of
time $T$ because, otherwise, the oscillation frequency $\Omega_{0}$ would not
be static, i.e., it would fluctuate around some mean frequency. In other
words, we are imposing that the derivative of $\phi(T)$ to be constant during
the locking, so we can say that $\Phi(T)=\delta_{\rho}T$ for some frequency
mismatch $\delta_{\rho}$ of our bare oscillator, i.e.,

$\phi(T)=(1+\delta_{\rho})T\quad\Rightarrow\quad\frac{d\phi}{dT}=1+\delta_{\rho},$
(46)

and then we can solve for $\varepsilon=\varepsilon(\delta_{\rho})$ for various
$\delta_{\rho}$. The cases $\rho=1$ and $\rho=2$ are

$\begin{cases}\delta_{1}=\frac{3\gamma
A^{2}}{8}-\varepsilon\left(\frac{F}{2A}-\frac{3\beta
A}{8}\right)\cos{\left(\delta_{1}T\right)},\\\ \\\ \delta_{2}=\frac{3\gamma
A^{2}}{8}+\varepsilon\left(\frac{\alpha}{4}+\frac{\gamma
A^{2}}{4}\right)\sin{\left(2\delta_{2}T\right)},\end{cases}\quad\Rightarrow\quad\begin{cases}\varepsilon(\delta_{1})>\left|\left(\delta_{1}-\frac{3\gamma
A^{2}}{8}\right)/\left(\frac{F}{2A}-\frac{3\beta A}{8}\right)\right|,\\\ \\\
\varepsilon(\delta_{2})>\left|\left(\delta_{2}-\frac{3\gamma
A^{2}}{8}\right)/\left(\frac{\alpha}{4}+\frac{\gamma
A^{2}}{4}\right)\right|,\end{cases}$ (47)

which defines a region in a $\varepsilon-\delta_{\rho}$ space which is, as we
would guess, the Arnold tongues. The AT maps using this approach are shown in
Supplementary Figure 14 for three different oscillation amplitudes
$A=\\{2,3,4\\}$, which has the same effects of ones obtained from the
experiment even after lots of approximations.

Supplementary Figure 14: Simulated Arnold tongues using Equation 45 for three
different oscillations amplitude. a) $A=2$; b) $A=3$; c) $A=4$. The softening
Duffing effect is enhanced as we increase the amplitude $A$, as we can verify
as the whole region is getting away from $\delta_{\rho}=0$ to the left, i.e.,
negative values of $\delta_{\rho}$.

The reduction of the oscillation frequency $\Omega_{0}$ as we increase the
amplitude $A$ is called softening Duffing, which is expected because we have
$\gamma<0$ (a consequence of the chosen optical detuning
$0<\Delta_{x}<\kappa/2$). If we had chosen other detuning $\Delta_{x}$, for
example, $\Delta_{x}>\kappa/2$ or $\Delta_{x}<0$, we would have $\gamma>0$ and
we would see the hardening Duffing effect, which is the shift of $\Omega_{0}$
to higher frequencies. Not only this but for high enough amplitudes we can
obtain larger $2:1$ AT than $1:1$, which is the case for $A=3$ and $A=4$,
showing that our model still has features about high-harmonic synchronization,
the same features of the experiment. To finalize this section, we conclude
that each of the terms $F,\alpha,\beta$ and $\gamma$ are directly proportional
to the tongue width $\Delta\Omega(p,q)$ with $p=\\{1,2,3,4\\}$ and $q=1$,
respectively, as we can see at the denumerator of Equation 47, i.e.,

$\begin{cases}\Delta\Omega(1,1)=2\times\left(\frac{F}{2A}-\frac{3\beta
A}{8}\right)\approx\frac{F}{A}\propto F\propto F_{0}\\\ \\\
\Delta\Omega(2,1)=2\times\left(\frac{\alpha}{4}+\frac{\gamma
A^{2}}{4}\right)\approx\frac{\alpha}{2}\propto\alpha\propto F_{1}\\\ \\\
\Delta\Omega(3,1)=2\times\frac{\beta A}{8}\propto\beta\propto F_{2}\\\ \\\
\Delta\Omega(4,1)=2\times\frac{\gamma A^{2}}{16}\propto\gamma\propto
F_{3}\end{cases}$ (48)

and these are the semi analytical expressions for the tongue width of each
harmonic, which we could engineer it to achieve wider Arnold tongues for
different harmonics choosing different $F$, $\alpha$, $\beta$ and $\gamma$ as
we design the geometry and the materials of our optomechanical cavity.

## Supplementary Note 7

Sidebands around the carrier at the synchronized region. To explain the
sidebands around the synchronization region we will first linearize Equation
39 and Equation 40 expanding $A(T)$ and $\Phi(T)$ as

$A(T)=\overline{A}+\delta
A(T)\quad\quad\text{and}\quad\quad\Phi(T)=\overline{\Phi}+\delta\Phi(T),$ (49)

and then diagonalize the linear part of the system

$\frac{d}{dT}\begin{pmatrix}\delta A\\\
\delta\Phi\end{pmatrix}=\begin{pmatrix}H_{AA}&H_{A\Phi}\\\ H_{\Phi
A}&H_{\Phi\Phi}\end{pmatrix}\begin{pmatrix}\delta A\\\
\delta\Phi\end{pmatrix},$ (50)

in which the actual form of $H_{AA}$, $H_{A\Phi}$, $H_{\Phi A}$ and
$H_{\Phi\Phi}$ are too big to be shown and not important for our present
analysis. The eigenvalues of this system give us the first order correction in
frequency and damping of our oscillator. The evolution of this system is

$\begin{pmatrix}\delta A(T)\\\
\delta\Phi(T)\end{pmatrix}=\begin{pmatrix}\delta A_{+}\\\
\delta\Phi_{+}\end{pmatrix}e^{\left(\lambda_{\text{Re}}^{+}+i\lambda_{\text{Im}}^{+}\right)T}+\begin{pmatrix}\delta
A_{-}\\\
\delta\Phi_{-}\end{pmatrix}e^{\left(\lambda_{\text{Re}}^{-}+i\lambda_{\text{Im}}^{-}\right)T},$
(51)

where $\lambda^{\pm}=\lambda_{\text{Re}}^{\pm}+i\lambda_{\text{Im}}^{\pm}$ are
the eigenvalues. Using these in $y(T)$ we can now search for sidebands, i.e.,

$y(T)=A(T)\sin{\left[T+\Phi(T)\right]}\approx\frac{1}{2i}\left(\overline{A}+\delta
A\right)\left[\left(1+i\delta\Phi\right)e^{iT}e^{i\overline{\Phi}}-\left(1-i\delta\Phi\right)e^{-iT}e^{-i\overline{\Phi}}\right]\approx\\\
\approx\overline{A}\sin{\left(T+\overline{\Phi}\right)}+\delta{A}(T)\sin{\left(T+\overline{\Phi}\right)}+\overline{A}\delta\Phi(T)\cos{\left(T+\overline{\Phi}\right)}.$
(52)

It is evident now that our oscillator has more than one single frequency
because of the product $\delta A(T)$ with the sine function and also because
of the product $\delta\Phi(T)$ with the cosine. The frequencies
$\Omega_{\text{SB}}^{\pm}$ and the linewidths $\Gamma_{\text{SB}}^{\pm}$ of
these new sidebands are given by

$\Omega_{\text{SB}}^{\pm}=\frac{\lambda_{\text{Im}}^{\pm}}{L_{\widetilde{t}}}\quad\quad\text{and}\quad\quad\Gamma_{\text{SB}}^{\pm}=\frac{\lambda_{\text{Re}}^{\pm}}{L_{\widetilde{t}}},$
(53)

and the graph of these using our semi-analytical model for the case $\rho=1$
is shown in Supplementary Figure 15, which has an excellent agreement with the
experimental data in both frequency and linewidth, showing us that these
sidebands are indeed a coupling between phase and amplitude of the mechanical
oscillator.

Supplementary Figure 15: a) Detuning from $\Omega_{0}/2\pi$ showing both
sidebands emerging as we increase the modulation depth; b) Linewidth of these
sidebands.

## Supplementary Note 8

1:1 Arnold Tongue Cusp Deviation. The purple cusp simulation in Fig.2 (d3)
from the main text (also shown again in Supplementary Figure 16(a)), which is
not present in the experimental data, was found to be a bifurcation in the
amplitude dynamics. To show this we must investigate Equation 39 in more
detail:

$\left<\frac{dA}{dT}\right>_{\phi}=\frac{\mu A}{2}\left(1-\frac{\sigma
A^{2}}{4}\right)+\\\
+\frac{\varepsilon\sin(\pi\omega)}{\pi}\left(\frac{\omega\left[F\left(\omega^{2}-9\right)+2\beta
A^{2}\right]\sin{\left[(\pi-\Phi)\omega\right]}}{\left(\omega^{2}-9\right)\left(\omega^{2}-1\right)}-\frac{A\left[\alpha\left(\omega^{2}-16\right)-6\gamma
A^{2}\right]\cos{\left[(\pi-\Phi)\omega\right]}}{\left(\omega^{2}-16\right)\left(\omega^{2}-4\right)}\right).$
(54)

We are searching for bifurcations around the $1:1$ Arnold tongue, so it’s
natural to expand such equation around $\omega=1$, identically as we’ve done
with the phase $\Phi$ in Equation 45, i.e.,

$\left<\frac{dA}{dT}\right>_{\phi}=\frac{\mu A}{2}\left(1-\frac{\sigma
A^{2}}{4}\right)-\varepsilon\left(\frac{F}{2}-\frac{\beta
A^{2}}{8}\right)\sin{\left[\Phi(T)\right]}.$ (55)

The fixed points of this equation are such that $\dot{A}=\dot{\Phi}=0$. As the
bifurcation occurs at the boundary of the Arnold tongue, we will assume a
fixed phase $\Phi$ of $\pi/2$ and $3\pi/2$ (which, respectively, minimize and
maximize $\sin{\Phi}$), such that the equation that we must solve for $A$ is

$0=\frac{\mu A}{2}\left(1-\frac{\sigma
A^{2}}{4}\right)\pm\varepsilon\left(\frac{F}{2}-\frac{\beta A^{2}}{8}\right).$
(56)

We know that every $n$-degree polynomial has $n$ complex roots, so we can plot
the absolute value of each of these roots as a function of $\varepsilon$ (we
will use the same $\mu$, $\sigma$, $F$, and $\beta$ parameters from the
simulation), and the result is shown in Supplementary Figure 16(b)

Supplementary Figure 16: a) Experimental $1:1$ Arnold tongue with its
simulated boundary in purple. The cusp occurs around -16 dBm, which
corresponds to $\varepsilon=1.9\%$; b) Amplitude evolution of each $A_{n}$
root from Equation 56 as the modulation depth $\varepsilon$ increases.

As we see, two of the roots ($A_{1}$ and $A_{2}$) of Equation 56 become
degenerate (the orange and blue curves of Supplementary Figure 16(b)) for a
modulation depth around $\varepsilon=1.83\%$. The cusp that we have at the
Supplementary Figure 16(a) happens around -16 dBm, which corresponds to a
modulation depth of $1.9\%$. We must not forget that Equation 56 is one of the
final steps of our semi-analytical analysis, and many simplifications were
involved. Nevertheless, the purple curve in Supplementary Figure 16 was
obtained using the full optomechanical model, and a difference of only
$0.07\%$ in $\varepsilon$ is observed. Similar cusps were observed in other
injection-locking experiments [8]. A more detailed analysis and more
experimental data should be gathered to understand why our experiment does not
exhibit this feature.

## Supplementary references

## References

  * [1] Markus Aspelmeyer, Tobias J. Kippenberg, and Florian Marquardt.  Cavity optomechanics.  Reviews of Modern Physics, 86(4):1391–1452, 12 2014.
  * [2] King Yan Fong, Menno Poot, Xu Han, and Hong X. Tang.  Phase noise of self-sustained optomechanical oscillators.  Phys. Rev. A, 90:023825, Aug 2014.
  * [3] M. L. Gorodetksy, A. Schliesser, G. Anetsberger, S. Deleglise, and T. J. Kippenberg.  Determination of the vacuum optomechanical coupling rate using frequency noise calibration.  Opt. Express, 18(22):23236–23246, Oct 2010.
  * [4] E.A. Jackson.  Perspectives of Nonlinear Dynamics: Volume 1.  Perspectives of Nonlinear Dynamics. Cambridge University Press, 1989.
  * [5] Cijy Mathai, Sunil A. Bhave, and Siddharth Tallur.  Modeling the colors of phase noise in optomechanical oscillators.  OSA Continuum, 2(7):2253–2259, Jul 2019.
  * [6] Shreyas Y. Shah, Mian Zhang, Richard Rand, and Michal Lipson.  Master-slave locking of optomechanical oscillators over a long distance.  Phys. Rev. Lett., 114:113602, Mar 2015.
  * [7] Keren Shlomi, D. Yuvaraj, Ilya Baskin, Oren Suchoi, Roni Winik, and Eyal Buks.  Synchronization in an optomechanical cavity.  Phys. Rev. E, 91:032910, Mar 2015.
  * [8] J. Simonet, M. Warden, and E. Brun.  Locking and arnold tongues in an infinite-dimensional system: The nuclear magnetic resonance laser with delayed feedback.  Phys. Rev. E, 50:3383–3391, Nov 1994.
  * [9] Mian Zhang, Gustavo Luiz, Shreyas Shah, Gustavo Wiederhecker, and Michal Lipson.  Eliminating anchor loss in optomechanical resonators using elastic wave interference.  Applied Physics Letters, 105(5):051904, Aug 2014.