# Frequency stabilization of self-sustained oscillations in a sideband-driven
electromechanical resonator

B. Zhang  Department of Physics, the Hong Kong University of Science and
Technology, Clear Water Bay, Kowloon, Hong Kong, China  William Mong Institute
of Nano Science and Technology, the Hong Kong University of Science and
Technology, Clear Water Bay, Kowloon, Hong Kong, China  Yingming Yan
Department of Physics, the Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  William Mong Institute of Nano
Science and Technology, the Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  X. Dong  Department of Physics,
the Hong Kong University of Science and Technology, Clear Water Bay, Kowloon,
Hong Kong, China  William Mong Institute of Nano Science and Technology, the
Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong
Kong, China  M. I. Dykman<EMAIL_ADDRESS>Department of Physics and
Astronomy, Michigan State University, East Lansing, Michigan 48824, USA  H. B.
Chan<EMAIL_ADDRESS>Department of Physics, the Hong Kong University of Science
and Technology, Clear Water Bay, Kowloon, Hong Kong, China  William Mong
Institute of Nano Science and Technology, the Hong Kong University of Science
and Technology, Clear Water Bay, Kowloon, Hong Kong, China

###### Abstract

We present a method to stabilize the frequency of self-sustained vibrations in
micro- and nanomechanical resonators. The method refers to a two-mode system
with the vibrations at significantly different frequencies. The signal from
one mode is used to control the other mode. In the experiment, self-sustained
oscillations of micromechanical modes are excited by pumping at the blue-
detuned sideband of the higher-frequency mode. Phase fluctuations of the two
modes show near perfect anti-correlation. They can be compensated in either
one of the modes by a stepwise change of the pump phase. The phase change of
the controlled mode is proportional to the pump phase change, with the
proportionality constant independent of the pump amplitude and frequency. This
finding allows us to stabilize the phase of one mode against phase diffusion
using the measured phase of the other mode. We demonstrate that phase
fluctuations of either the high or low frequency mode can be significantly
reduced. The results open new opportunities in generating stable vibrations in
a broad frequency range via parametric downconversion in nonlinear resonators.

Subject Areas

Mechanics, Nonlinear Dynamics, Nanophysics

##  I Introduction

Self-sustained oscillations play a major role in diverse fields ranging from
biological systems to lasers and clocks [1]. One of the most important classes
of systems that can display self-sustained oscillations is mechanical
resonators. Examples include quartz oscillators that have long served as time
references in ordinary clocks. In recent years, silicon micromechanical
devices [2, 3, 4, 5], because of their small size, have become viable
alternatives of quartz oscillators in numerous applications, from
telecommunications to smart watches. Self-oscillating micro- and
nanomechanical devices are also used in high-precision measurements. Shifts in
the vibration frequency enable measuring force [6], charge [7], spin [8] and
mass [9]. A well-known example is frequency modulation atomic force microscopy
[10]. The decrease of the size of mechanical resonators is generally expected
to lead to an increase of fluctuations [11]. Suppressing fluctuations is a
central problem in micro- and nanomechanical systems, and in particular in the
systems that display self-sustained oscillations.

In micromechanical devices, the feedback needed to excite self-sustained
oscillations is often implemented via an external circuit in which the signal
generated by vibrations is phase-shifted and amplified to serve as the
periodic drive [12, 13, 14, 15, 16, 17, 18, 19]. Alternatively, the feedback
mechanism can also be intrinsically built into the driven system, such as by
photothermal [20, 21, 22, 23, 24] or other effects [25], or by measurement
backaction [26].

Phase fluctuations in self-sustained oscillations are induced by different
sources, from the thermal noise that arise from the coupling of the resonator
to the environment described by the fluctuation-dissipation theorem to noise
in various parts of the feedback circuit [27, 28, 29]. There have been many
efforts to improve the phase stability of self-sustained oscillations in
micro- and nano-mechanical systems. In the linear regime, phase fluctuations
induced by thermal noise decreases with the oscillation amplitude [30, 31,
32]. It is therefore beneficial to maximize the oscillation amplitude by
operating the feedback circuit with large external drive. However, if the
oscillation amplitude is increased beyond the linear regime, the
eigenfrequency of the resonator becomes dependent on the amplitude. As a
result, amplitude fluctuations are converted into frequency fluctuations. It
was shown that, by designing the dependence of the eigenfrequency on amplitude
to be non-monotonic, phase noise can be minimized at the extremum where the
eigenfrequency is locally independent of amplitude [33, 34].

The nonlinear coupling to high order modes has also been exploited to reduce
phase noise. In these systems, the two modes are tuned into internal
resonance, with the ratio of the eigenfreqeuncies of the modes being an
integer. With energy exchanged between the modes, the phase fluctuations of
the self-sustained oscillations of the lower mode become much weaker [35, 36].

Other approaches focus on mitigating the effects of noise generated by the
feedback loop. For example, noise in the phase of the periodic drive can
produce phase noise in the oscillations. This effect can be significantly
reduced by driving a resonator into nonlinearity and picking an operation
point where the vibration frequency is insensitive to the phase of the drive
[37, 38]. More recent works extended these concepts [39, 40] and explored the
possibility of surpassing the limit imposed by the thermomechanical noise
[41]. Similar approaches were also used in minimizing the phase noise in self-
sustained oscillations of two vibrational modes under nondegenerate parametric
drive [42].

An important distinction of self-sustained oscillations from forced
oscillations is that the oscillation phase is arbitrary. In the presence of
noise, fluctuations of the phase accumulate in time, resulting in phase
diffusion [43, 30, 44, 45, 46, 47, 48]. Phase diffusion leads to decoherence,
determining the spectral width of self-sustained oscillations. Minimizing
phase diffusion is therefore paramount in applications of self-sustained
oscillations ranging from high sensitivity detection to high stability
frequency standards. To our knowledge, the possibility of significant
reduction of phase diffusion in self-sustained oscillations has neither been
explored theoretically nor demonstrated experimentally in micro- and nano-
mechanical oscillators.

In this paper, we present a scheme to stabilize the frequency of an important
type of self-sustained oscillations induced by a drive with a frequency that
is orders of magnitude different from the frequency of the vibrations of
interest. We consider the vibrations induced by dynamical backaction in a
sideband-driven micromechanical resonator with two nonlinearly coupled
vibrational modes. The two modes have vastly different eigenfrequencies and
decay rates, with the higher mode serving as an analog of a photon cavity mode
[47, 49] in the context of cavity optomechanics [50, 51]. Sufficiently strong
pump applied at the blue-detuned sideband of the upper mode excites self-
sustained vibrations [52, 53, 54, 55, 56, 57, 58, 59] of both modes. In the
stationary regime, the phase diffusion of the modes is anti-correlated. The
sum of the phases remains essentially constant [47].

We show that the sum of the phases can be adjusted by the phase of the
sideband pump and study the transient process in which such adjustment occurs
when the pump phase is changed. For a step change of the pump phase, the
phases of both modes settle to new values after a transient, the duration of
which is determined by the smallest eigenvalue of the linearized equations of
motion about the stable state. The phase change of each mode is proportional
to the pump phase change, which is assumed to be small. Importantly, the two
proportionality constants add up to unity. This finding, together with the
phase anti-correlation of the two modes, allow us to stabilize the phase of
one mode by measuring the phase of the other mode and then compensating for
the phase diffusion by adjusting the phase of the pump. We demonstrate that
phase fluctuations of either the high or low frequency mode can be
significantly reduced. This results in a much narrower spectral linewidth. Our
scheme does not require a frequency reference near the mode to be stabilized,
in contrast to direct feedback that stabilizes a particular mode by measuring
its phase. By exploiting the nonlinear coupling between two modes, the
frequency reference could be orders of magnitude different from the
eigenfrequency of the mode that is stabilized.

Figure 1: (a) A sketch of the spectra of the mode 1 (mechanical mode) and mode
2 (cavity mode). Energy of the pump that leads to self-sustained vibrations is
supplied to the upper sideband of mode 2. (b) Scanning electron micrograph of
electromechanical resonator including the measurement scheme. The white scale
bar at the lower left corner measures 30 $\mu$m. (c) Vibration profiles of
mode 1 (top) and mode 2 (bottom). (d) Dependence of the square of vibration
amplitude of mode 1 on the detuning of the pump frequency from the upper
sideband for mode 1 at pump current amplitude of 189 $\mu$A. (e) Similar plot
for mode 2.

##  II TWO-MODE ELECTROMECHANICAL SYSTEM

###  II.1 Excitation and detection of self-sustained vibrations

In our experiment, the micromechanical resonator consists of a square
polysilicon plate ($100\,\mu$m$\,\times\,100\,\mu$m$\,\times\,3.5\,\mu$m) that
is suspended by two beams of the same cross-section area ($1.3\,\mu$m$\
\mathrm{by}\ 2\,\mu$m) but different lengths of $80\,\mu$m and $75\,\mu$m
respectively [Fig. 1(b)]. The beams are covered by 30 nm of gold to reduce the
electrical resistance. We utilize two vibrational modes in our experiment. For
the low frequency mode (mode 1), the plate undergoes translational motion
normal to the substrate as the two beams are deformed [top of Fig. 1(c)]. The
resonance frequency $\omega_{1}/2\pi$ is 47030.7 Hz and the damping constant
$\Gamma_{1}/2\pi$ is 0.48 Hz. The high frequency mode (mode 2) is the lowest
in-plane mode of the shorter beam [bottom of Fig. 1(c)], with resonant
frequency $\omega_{2}/2\pi=1867195.4$ Hz and damping constant
$\Gamma_{2}/2\pi=36.2$ Hz. Because of its significantly larger decay rate, in
the presence of a weak external drive mode 2 can play a role similar to
optical/microwave modes in cavity optomechanics [50, 51]. The change of the
tension caused by the vibrations leads to a nonlinear coupling between the
modes.

Forced vibrations of each mode can be excited by periodic electrostatic forces
when small probe ac voltages are applied to the electrodes nearby. For mode 2,
vibrations can also be excited by the periodic Lorentz force when an ac
current is passed through the beams in a magnetic field perpendicular to the
substrate. Motion of the plate in mode 1 is detected by measuring the
capacitance change between the plate and the underlying electrode. Vibrations
in mode 2 are detected by the change in the current due to the back
electromotive force as the beam vibrates in the magnetic field. All
measurements are carried out at a temperature of 4 K and pressure of
$<10^{-5}$ torr. These measurements allowed us to determine the values of the
eigenfrequencies and decay rates of the modes.

To induce self-sustained oscillations, a periodic pump current at frequency
$\omega_{F}$ close to the frequency sum ($\omega_{1}+\omega_{2}$) is applied
through the beam to generate a Lorentz force as shown in Fig. 1(a). This force
modulates the nonlinear coupling between the two modes. The equations of
motion that describe the effect of such resonant modulation on the coupled-
mode system are:

$\displaystyle{\ddot{q}}_{1}+\omega_{1}^{2}q_{1}+2\mathrm{\Gamma}_{1}{\dot{q}}_{1}+\frac{\gamma_{1}}{m_{1}}q_{1}^{3}+\frac{\gamma}{m_{1}}q_{1}q_{2}^{2}$
$\displaystyle=\frac{F}{m_{1}}q_{2}cos{\left(\omega_{F}t+\theta_{F}\right),}$
$\displaystyle{\ddot{q}}_{2}+\omega_{2}^{2}q_{2}+2\mathrm{\Gamma}_{2}{\dot{q}}_{2}+\frac{\gamma_{2}}{m_{2}}q_{2}^{3}+\frac{\gamma}{m_{2}}q_{2}q_{1}^{2}$
$\displaystyle=\frac{F}{m_{2}}q_{1}cos{\left(\omega_{F}t+\theta_{F}\right),}$
(1)

where $q_{1,2}$ are the displacements of the plate and the beam center,
respectively, $m_{1,2}$ are the effective masses, and $\gamma_{1,2}$ are the
Duffing nonlinearities of the two modes. The two modes are found to be
dispersively coupled, where the vibration of one mechanical mode creates
tension in the beams and in turn modifies the resonance frequency of the other
mode (Appendix A). In Eq. (II.1), the dispersive coupling is represented by
the terms with coupling constant $\gamma$ that correspond to coupling energy
$\frac{1}{2}\gamma q_{1}^{2}q_{2}^{2}$. Parameter $F$ gives the amplitude of
the parametric sideband pumping that is proportional to the amplitude of the
applied pump current $I_{pump}$. The phase of the pump is $\theta_{F}$. In the
Hamiltonian of the coupled modes, the term corresponding to the pump is
$-Fq_{1}q_{2}\cos{\left(\omega_{F}t+\theta_{F}\right)}$.

We assume that the pump frequency is close to $\omega_{1}+\omega_{2}$, with
the pump detuning $\Delta_{F}=\omega_{F}-(\omega_{1}+\omega_{2})$ much smaller
than $\omega_{F}$. By making the substitution
$u_{1}=\frac{1}{2}(q_{1}-i\omega_{1}^{-1}{\dot{q}}_{1})e^{-i\omega_{1}t}$ to
convert $q_{1}$ and $q_{2}$ to slowly varying complex vibration amplitudes
$u_{1}$ and $u_{2}$, and by applying the rotating wave approximation, we
obtain equations for $u_{1,2}$:

$\displaystyle{\dot{u}}_{1}+\Gamma_{1}u_{1}-i\frac{3\gamma_{11}}{2\omega_{1}}u_{1}|u_{1}|^{2}-i\frac{\gamma_{12}}{\omega_{1}}u_{1}|u_{2}|^{2}$
$\displaystyle=\frac{F_{1}e^{i\theta_{F}}}{4i\omega_{1}}u_{2}^{\ast},$
$\displaystyle{\dot{u}}_{2}+(\Gamma_{2}+i\Delta_{F})u_{2}-i\frac{3\gamma_{22}}{2\omega_{2}}u_{2}|u_{2}|^{2}-i\frac{\gamma_{21}}{\omega_{2}}u_{2}|u_{1}|^{2}$
$\displaystyle=\frac{F_{2}e^{i\theta_{F}}}{4i\omega_{2}}u_{1}^{\ast},$ (2)

where parameters $\gamma_{11}$ ($1.91\times
10^{22}\mathrm{rad^{2}s^{-2}m^{-2}}$), $\gamma_{22}$ ($7.38\times
10^{25}\mathrm{rad^{2}s^{-2}m^{-2}}$), $\gamma_{12}$ ($8.41\times
10^{22}\mathrm{rad^{2}s^{-2}m^{-2}}$), $\gamma_{21}$ ($1.26\times
10^{25}\mathrm{rad^{2}s^{-2}m^{-2}}$) and $F_{1,2}$ contain not only direct
contributions from $\gamma_{1}$, $\gamma_{2}$, $\gamma$ and $F$ respectively,
but also account for the renormalization that arises from the nonlinear
coupling of the modes [47] and from modulating not just the mode coupling, but
applying force to one or the other mode [47]; we note that
$F_{1}\textfractionsolidus
F_{2}=\gamma_{12}\textfractionsolidus\gamma_{21}=m_{2}/m_{1}$.

We first consider the case in which vibration amplitudes are small so that the
terms involving $\gamma_{11}$, $\gamma_{22}$, $\gamma_{12}$ and $\gamma_{21}$,
can be neglected. With $\Gamma_{2}>>\Gamma_{1}$, we apply the adiabatic
approximation by setting ${\dot{u}}_{2}=0$ in the second equation in Eq.
(II.1) because after a short transient, both $|{\dot{u}}_{1}/{u_{1}}|$ and
$|{\dot{u}}_{2}/{u_{2}}|$ become much smaller than
$\left({\mathrm{\Gamma}_{2}}^{2}+\Delta_{F}^{2}\right)^{1/2}$. It follows from
Eq. (II.1) that for $\Delta_{F}=0$,
$u_{2}\approx\frac{F_{2}e^{i\theta_{F}}}{4i\omega_{2}\mathrm{\Gamma}_{2}}u_{1}^{\ast}$.
Putting $u_{2}$ back into the first equation yields an effective damping
constant of mode 1 that decreases with the pumping power as ${\Gamma_{1}-\
F}_{1}F_{2}/16\omega_{1}\omega_{2}\Gamma_{2}$.

For sufficiently large pumping power, the effective damping becomes negative
and self-sustained vibrations are excited. In fact, both modes start
oscillating at the same time, although with different amplitudes. Figures 1(d)
and 1(e) shows these amplitudes for modes 1 and 2 respectively as a function
of $\Delta_{F}$ at pump current amplitude of 189$\ \mu$A. The vibrations are
excited when $\Delta_{F}$ is increased beyond a threshold value $-\omega_{c}$.
This value can be found from the condition of the loss of stability of the
stationary state $u_{1}=u_{2}=0$, giving [47]

$\displaystyle\omega_{c}=\left(\Gamma_{1}+\
\Gamma_{2}\right)\left[\Xi^{2}-1\right]^{\frac{1}{2}\ },$
$\displaystyle\Xi=\left(F_{1}F_{2}/16\omega_{1}\omega_{2}\Gamma_{1}\Gamma_{2}\right)^{1/2}.$
(3)

The vibration amplitude increases with $\Delta_{F}$. At
$\Delta_{F}=\omega_{c}$, the zero-amplitude state becomes stable again.
Previous experiments [47] have determined the dependence of $\omega_{c}$ on
the pump power and studied the bistability and hysteretic behavior as
$\Delta_{F}$ is varied.

###  II.2 Anti-correlated phase diffusion of the two modes

If the pump is not too strong, self-sustained vibrations of the two modes are
sinusoidal, but their frequency deviates from $\omega_{1}$ and $\omega_{2}$.
The deviation depends on the amplitude and frequency of the pump. To analyze
the self-sustained vibrations, we write the variables $u_{1},u_{2}$ in the
form of the amplitudes $r_{1},r_{2}$ and phases $\phi_{1},\phi_{2}$ of the
vibrations and take into account that the vibrations in the rotating frame
have frequency $\delta\omega$,

$\displaystyle u_{1}=r_{1}e^{i\phi_{1}(t)}\exp(i\delta\omega\,t),\quad
u_{2}=r_{2}e^{i\phi_{2}(t)}\exp(-i\delta\omega\,t).$ (4)

In the steady state, $r_{1,2}$ and $\phi_{1,2}$ are independent of time. In
the laboratory frame, self-sustained oscillations in modes 1 and 2 occur at
frequencies of $\omega_{\mathrm{self1}}=\omega_{1}+\delta\omega$ and
$\omega_{\mathrm{self2}}=\omega_{2}+\Delta_{F}-\delta\omega$, respectively.

It is convenient to write the equations of motion in the form

$\displaystyle\dot{r}_{i}=-\Gamma_{i}r_{i}+\frac{F_{i}}{4\omega_{i}}\frac{r_{3-i}}{r_{i}}\sin\Theta\quad(i=1,2),$
$\displaystyle\dot{\phi}_{+}=-\Delta_{F}+\sum_{i=1,2}\left[\frac{\gamma^{(+)}_{i\,3-i}}{\omega_{i}}r_{3-i}^{2}-\frac{F_{i}}{4\omega_{i}}\frac{r_{3-i}}{r_{i}}\cos\Theta\right]$
(5)

and

$\displaystyle\dot{\phi}_{-}=\Delta_{F}-2\delta\omega-\sum_{i=1,2}(-1)^{i}\left[\frac{\gamma^{(-)}_{i\,3-i}}{\omega_{i}}r_{3-i}^{2}-\frac{F_{i}}{4\omega_{i}}\frac{r_{3-i}}{r_{i}}\cos\Theta\right],$
(6)

Here

$\phi_{\pm}=\phi_{1}\pm\phi_{2},\quad\Theta=\theta_{F}-\phi_{+}.$

The explicit form of the coefficients $\gamma^{(\pm)}_{ij}$ is given in
Appendix B.

The stationary solution of Eqs. (II.2) and (6) gives the amplitudes
$r_{1,2}^{(0)}$ and the phases $\phi_{1,2}^{(0)}$ of the self-sustained
vibrations of the coupled modes, as well as the vibration frequency in the
rotating frame $\delta\omega$. They are also listed in Appendix B. In
particular, the total phase is

$\displaystyle\phi_{+}^{(0)}=\theta_{F}-\Theta^{(0)},\quad\Theta^{(0)}=\arcsin(1/\Xi),$
(7)

where $\Xi$ is defined in Eq. (II.1).

An important feature of Eq. (II.2) is that the evolution of the amplitudes
$r_{1,2}$ and the total phase $\phi_{+}$ is independent of the phase
difference $\phi_{-}$. At the same time, the evolution of $\phi_{-}$ depends
on $r_{1,2}$ and $\phi_{+}$ as well as the phase $\theta_{F}$ of the driving
field but not on $\phi_{-}$ itself. It is this feature that enables an
efficient control of the mode phase.

While fluctuations in $\phi_{1}+\phi_{2}$ do not accumulate, fluctuations in
$\phi_{1}-\phi_{2}$ do. Figure 2(a) plots the measured phase change
${\delta\phi}_{1,2}(t)=\phi_{1}(t)-\phi_{1}(0)$ as a function of time as
$\theta_{F}$ is kept constant. ${\delta\phi}_{1}(t)$ is almost perfectly anti-
correlated with ${\delta\phi}_{2}(t)$, in agreement with Eq. (7). Previous
works have demonstrated that the system undergoes phase diffusion, with the
mean square phase change increasing with time [47]. Depending on the origin of
phase fluctuations, the diffusion was found to be either normal or anomalous.
This phase diffusion leads to decoherence, determining the linewidth of the
self-sustained vibrations. Figures 2(b) and (c) plot the results of Fig. 2(a)
in frames rotating at the frequencies of self-sustained vibrations of the two
modes, $\omega_{\mathrm{self1}}$ and $\omega_{\mathrm{self2}}$ respectively.
The motion of two modes in the tangential direction is not confined, but is
correlated.

Figure 2:  (a) Change of the vibration phase ${\delta\phi}_{1}$ of mode 1
(green) and ${\delta\phi}_{2}$ of mode 2 (yellow), measured as a function of
time. The black line represents the sum ${\delta\phi}_{1}+{\delta\phi}_{2}$.
The pump current amplitude is 189 $\mu$A corresponding to $F$ = 2.3 mN$m^{-1}$
and frequency detuning $\Delta_{F}$ is 300 Hz. The same results are plotted in
panels (b) and (c) in the frames rotating at the frequencies
$\omega_{\mathrm{self1}}$ and $\omega_{\mathrm{self2}}$ of the self-sustained
vibrations of the two modes. The variables are $X_{i}=r_{i}\cos(\phi_{i})$ and
$Y_{i}=r_{i}\sin(\phi_{i})$, $i=1,2$.

##  III FREQUENCY STABILIZATION

With the phase anti-correlation of the two modes, one can measure the phase
diffusion of one mode and infer the phase of the other mode without measuring
it. This property enables the design of a scheme that uses the measured phase
of one mode to stabilize the phase of the other mode, provided that there is a
mechanism to alter the latter. We show below that either $\phi_{1}$ or
$\phi_{2}$ can be stabilized by adjusting the pump phase $\theta_{F}$ and
using the measured $\phi_{2}$ or $\phi_{1}$ as the error signal. An antecedent
less efficient scheme was implemented in a previous experiment [47] to prove
the concept for stabilizing mode 2. Here, we analyze in detail,
quantitatively, the effect on the two modes when the pump phase $\theta_{F}$
is adjusted. Then we develop and optimize the stabilization scheme in our
experimental system. Furthermore, we demonstrate that mode 1 can also be
stabilized by measuring mode 2, despite requiring much larger changes in
$\theta_{F}$.

###  III.1 Effect of a change of pump phase

For short durations of phase accumulation the noise-induced changes of the
phases are small and can be compensated by small changes of $\theta_{F}$. The
effect can be analyzed by linearizing Eqs. (II.2) about the stationary
solution $r_{1,2}^{(0)},\phi_{+}^{(0)}$. It is convenient to write the
linearized equations for the increments

$\displaystyle x_{1}=\zeta-\zeta^{(0)},\quad
x_{2}=\theta_{F}-\phi_{+}-\Theta^{(0)},\quad
x_{3}=\frac{r_{1}-r_{1}^{(0)}}{r_{1}^{(0)}},$

where $\zeta\equiv r_{2}/r_{1}$. The corresponding equations have the form

$\displaystyle\dot{x}_{n}=\Lambda_{nm}x_{m}\quad(n,m=1,2,3),$ (8)

where we imply summation over repeated indices. Appendix B gives the explicit
form of the matrix $\hat{\Lambda}$. The initial conditions that corresponds to
incrementing $\theta_{F}$ by $\delta\theta_{F}$ at $t=0$ are
$x_{1}(0)=x_{3}(0)=0,x_{2}(0)=\delta\theta_{F}$.

For stable self-sustained vibrations described by Eq. (II.2) the eigenvalues
$\lambda_{\nu}$ ($\nu=1,2,3$) of the matrix $\hat{\Lambda}$ have negative real
parts. Therefore $x_{1,2,3}(t)$ become exponentially small for $t\gg t_{r}$,
where the relaxation time is $t_{r}=1/\min|\mathrm{Re}\,\lambda_{\nu}|$. This
means that, over the relaxation time, $\phi_{+}$ is incremented by
$\delta\theta_{F}$, i.e., $\delta\phi_{+}(t)\to\delta\theta_{F}$. For the rest
of the paper, we use $\delta\phi_{\pm}(t)$ and $\Delta\phi_{\pm}$ to denote
the transient and settled values of $\phi_{\pm}$ changed from the initial
values respectively.

The eigenvalues $\lambda_{k}$ can be found analytically near the bifurcation
point $\Delta_{F}=-\omega_{c}$ where self-sustained vibrations emerge, see
Appendix B. Here the eigenvalues $\lambda_{1,2}$ are complex conjugate, with
the real part $-(\Gamma_{1}+\Gamma_{2})$. The eigenvalue $\lambda_{3}$ is real
and small in the absolute value, $\lambda_{3}\propto\Delta_{F}+\omega_{c}$.
Our numerical analysis shows that, for the parameters of our system, the
eigenvalue $\lambda_{3}$ remains small in the absolute value throughout the
range we have explored, see Fig. 3 (e) and (f). Therefore the relaxation time
is $t_{r}=-1/\lambda_{3}$.

The increment $\delta\phi_{-}$(t) of $\phi_{-}$ is described by the linearized
equation (6), with $d\delta\phi_{-}/dt$ being a sum of the appropriately
weighted $x_{1,2,3}(t)$. Hence, $\delta\phi_{-}$(t) is given by the integrals
over time of $x_{1,2,3}(t)$. Since all variables $x_{n}(t)$ are proportional
to $\delta\theta_{F}$ and all of them decay for $t\gg t_{r}$, we see that
$\delta\phi_{-}(t)\to C\delta\theta_{F}$ for $t\gg t_{r}$, where $C$ is
determined by the parameters of the system, see Appendix C. Therefore
incrementing $\theta_{F}$ leads, after the transient time $\gtrsim t_{r}$, to

$\displaystyle\Delta\phi_{2}=g\delta\theta_{F},\quad\Delta\phi_{1}=(1-g)\delta\theta_{F},\quad
g=(1-C)/2.$ (9)

####  III.1.1 Implementing phase increment of the driving field

The increment of the phase $\theta_{F}$ of the driving field is obtained by
changing the driving frequency $\omega_{F}$ for a short time and then bringing
it back to the original value. One scenario is where this time is of order
$1/\omega_{F}$. Equations (II.1), (II.2), and (6) do not apply on this time
scale. However, the change of the “slow” variables over time $\sim
1/\omega_{F}$ associated with the change of $\theta_{F}$ is small and can be
disregarded. It applies to the case implemented in the experiment where
$\theta_{F}$ is controlled by a digital function generator.

Alternatively, the phase can be incremented by changing the frequency
$\omega_{F}$ over time much longer than $1/\omega_{F}$ but yet much shorter
than $1/\max|\lambda_{\nu}|$. The dynamics in this case is described by Eqs.
(II.2), and (6). It can be modeled by assuming that $\Delta_{F}$ is changed by
$\delta\theta_{F}\delta(t)$. As seen from Eq. (II.2) this is equivalent to the
initial condition $x_{2}(t\to+0)=\delta\theta_{F}$. However, $\phi_{-}$ is
also incremented by $\delta\theta_{F}$.

Yet another way of changing the phases is to apply a pulse of $\Delta_{F}(t)$
with $\Delta_{F}(t)$ varying slowly on the scale $t_{r}$. The change of the
variables $x_{n},\phi_{-}$ can be then described in the adiabatic
approximation. It leads to simple explicit expressions for this change, see
Appendix C. Unexpectedly, the result allows one also to express the parameter
$C$ in Eq. (9) in a simple explicit form, which does not require diagonalizing
the matrix $\hat{\Lambda}$.

We note that the frequency of self-sustained vibrations can be obtained in an
alternative way. One could omit $\delta\omega$ in defining the variables
$r_{1,2},\phi_{1,2}$ in Eq. (4). This omission will not affect the equations
of motion (II.2) for $r_{1,2},\phi_{+}$, that is, they will have the same
form. The stationary solution of these equations gives
$r_{1,2}^{(0)},\phi_{+}^{(0)}$, and this solution is stable for
$\Delta_{F}+\omega_{c}>0$, since Re $\lambda_{1,2,3}<0$. On the other hand,
Eq. (6) will have a solution $\dot{\phi}_{-}$ equals constant. The value of
this constant gives $-2\delta\omega$, see Eq. (B). In the stationary state,
the complex amplitudes $u_{1,2}(t)$ oscillate at frequencies $\delta\omega$
and $-\delta\omega$, respectively.

####  III.1.2 Measurement of the effect of a step change of the pump phase

Figures 3(a) and 3(b) show that both the amplitude and phase are perturbed by
a step change in the pump phase, for modes 1 and 2 respectively. Details of
the calculations of the transient are described in Appendix C. There are two
stages in the transient response. The system first relaxes quickly following
the spirals, the decay rate and frequency of which are determined by the real
and imaginary parts of the eigenvalues $\lambda_{1,2}$ respectively. As
described earlier in this section, $\lambda_{1,2}$ are complex conjugates of
each other. The second part of the transient is marked by the arrow pointing
to the left in Fig. 3(a) for mode 1. On this part of the phase trajectory the
system approaches the new stable state at a much slower rate given by
$\lambda_{3}$ and the fast variables follow the slow variable adiabatically.
For mode 2, the second stage of relaxation occurs over a much smaller range in
Fig. 3(b) and is difficult to identify.

After the transient, the amplitudes $r_{1,2}$ of self-sustained vibrations go
back to the original values. The phases change by ${\Delta\phi_{1,2}}$ as
given by Eq. (9), assuming that phase diffusion induced by noise is
negligible. In particular, ${\Delta\phi}_{1}+{\Delta\phi}_{2}$ equals
${\delta\theta}_{F}$ after the transient, according to Eq. (7).

As described later, the algorithms to stabilize the phases $\phi_{1,2}$
requires knowing the ratio $g={\Delta\phi}_{2}/\delta\theta_{F}$ that denotes
the fraction of the change in pump phase taken up by mode 2. The fraction of
the change in pump phase taken up by mode 1 is
$1-g={\Delta\phi}_{1}/\delta\theta_{F}$. The value of the parameter $g$ is
determined by Eqs. (9) and (28). We find that $g$ depends only on system
parameters including the damping constants, effective masses and nonlinear
coupling coefficients, but is independent of the pump detune frequency and
amplitude (Appendix C). For our system parameters, $g=0.94$. A change in the
phase of the pump therefore affects both mode 1 and mode 2. However, the
effect on mode 2 far exceeds that on mode 1.

We measure the response of the two modes for a duration of 0.5 s after
$\theta_{F}$ is abruptly changed by $1^{\circ}$. $\theta_{F}$ is controlled by
a digital function generator. It changes over a time scale
$\sim\omega_{F}^{-1}$, as indicated above, much shorter than
$|\lambda_{1,2,3}^{-1}|$. The measurement is repeated 7200 times and the
averaged phase changes scaled by the pump step,
${\delta\phi}_{1,2}(t)/\delta\theta_{F}$, are plotted in Figs. 3(c) and 3(d),
respectively, as green and yellow dots at detuning $\Delta_{F}=$ 1000 Hz and
amplitude of 379 $\mu$A of the pump current that corresponds to $F$ = 4.5
mN$m^{-1}$. Averaging is necessary to resolve the change from detection noise
and the noise introduced by phase diffusion. The agreement between measurement
and theory is good on the slow-relaxation time scale, which in our case was
0.3 s, with no fitting parameters. However, it is not feasible for us to
measure the fast relaxation represented by the solid lines in the insets, as
it requires a much larger bandwidth for the lock-in amplifier that
significantly increases the detection noise.

Apart from the transient response, the measured values of ${\Delta\phi}_{1,2}$
in the steady state also agrees with calculations: mode 2 accounts for most of
the phase change (94%) in response to $\delta\theta_{F}$ while the remaining
6% occurs in mode 1. As plotted in Fig. 3(e), the green dashed line represents
the eigenvalue $\lambda_{3}$ that determines the rate of the slow relaxation.
$|\lambda_{3}|$ first increases with increasing pump detuning, then remains
largely constant when pump detuning exceeds 2 kHz. The red and blue lines plot
the real parts of $\lambda_{1,2}$ that determine the rate of the fast
relaxation. Figure 3(f) plots the imaginary parts of eigenvalues
$\lambda_{1,2}$ that are associated with the frequency of spiraling in Figs.
3(a) and 3(b). The relaxation rates change for other pump detuning values.

Figure 3: (a) Numerical calculation of the dynamics of mode 1 in the frame
rotating at frequency $\omega_{\mathrm{self1}}$ in response to a sudden change
of $\theta_{F}$ by $1^{\circ}$; $X_{1}=r_{1}\cos\phi_{1}$ and
$Y_{1}=r_{1}\sin\phi_{1}$. The pump current amplitude is 379 $\mu$A and
frequency detuning $\Delta_{F}$ is 1000 Hz. Initially, $X_{1}$ is set to zero.
The black line is part of the circles
$X_{1}^{2}+Y_{1}^{2}=r_{1}^{(0)}{}^{2}$.(b) Corresponding plot for mode 2. (c)
Measured dependence of $\phi_{1}$ on time (green circles). The solid line
represents numerical simulations with no fitting parameters. Inset: zoom in at
shorter times. (d) Corresponding plot for mode 2. (e) Real parts of the
eigenvalues for matrix $\mathbf{\Lambda}$ in the linearized equation of motion
plotted as a function of the pump detuned frequency. $\lambda_{1,2,3}$ are
plotted as red, blue and green dashed lines respectively. Near the
supercritical Hopf bifurcation point, analytical expressions (Appendix B) give
the green solid line for $\lambda_{3}$ and the red square for $\lambda_{1}$.
(f) Similar plot for the imaginary parts of $\lambda_{1,2,3}$.

The ability to control the phase of the modes by adjusting $\delta\theta_{F}$,
together with the near perfect anti-correlation in the phase diffusion of
modes 1 and 2, allows us to implement a scheme to stabilize the phase of mode
2/1 over time scales longer than $|\lambda_{3}^{-1}|$, by measuring the phase
of mode 1/2 that vibrates at a different frequency.

##  IV FREQUENCY STABILIZATION

###  IV.1 Stabilization Algorithm

We first discuss how the phase of mode 2 can be stabilized by measuring mode
1. Figure 4(a) shows the stabilization algorithm. The procedure starts at time
t = 0, at which we define $\phi_{1}\left(t=0\right)=\phi_{2}\left(t=0\right)=\
\ \theta_{F}\left(t=0\right)=0$ for convenience. In each cycle, $\phi_{1}$ is
measured. Then the pump phase $\theta_{F}$ is adjusted accordingly in the
subsequent cycle to bring $\phi_{2}$ close to zero. The duration of each cycle
$\Delta t$ is chosen to exceed the relaxation time $|\lambda_{3}^{-1}|$. The
$p^{th}$ cycle lies in the time interval $t_{p}<t<t_{p+1}$, where
$t_{p}=p\Delta t$ and $p$ is a positive integer. $\theta_{F}$ is modified only
at the beginning of a cycle at $t_{p}$, to a value $\theta_{F,p}$ determined
by the algorithm described below. It remains fixed at other times.

Figure 4: (a) Algorithm for stabilizing $\phi_{2}$ by measuring $\phi_{1}$ and
adjusting $\theta_{F}$. (b) The duration of one cycle equals
$t_{\mathrm{wait}}+t_{\mathrm{measure}}$. The $p^{th}$ cycle lasts from time
$t_{p}$ to $t_{p+1}$. $\theta_{F}$ is changed only at the beginning of each
cycle. (c) $\phi_{2}$ is measured as a function of time when the stabilization
routine is turned off for 20 min (yellow) and on for 20 min (brown). The pump
current amplitude is 189 $\mu$A and frequency detuning is 300 Hz. Inset: The
spectra of the self-sustained vibrations of mode 2 when the stabilization
algorithm in (b) is implemented. The spectral width is significantly reduced
with stabilization turned on. (d) Dependence of the standard deviation of
$\phi_{2}$ on $t_{\mathrm{measure}}$ at a fixed $t_{\mathrm{wait}}$ of 300 ms.

Let $\phi_{1,\mathrm{diffuse}}\left(t\right)$ and
$\phi_{2,\mathrm{diffuse}}\left(t\right)$ represent the phase change of modes
1 and 2 respectively that arises solely due to phase diffusion. In the
$p^{th}$ cycle, after the system has settled in response to changes in
$\theta_{F}$ that occurred at $t_{p}$, the phases of the two modes are given
by:

$\displaystyle\phi_{1}\left(t\right)=\phi_{1,\mathrm{diffuse}}\left(t\right)+\left(1-g\right)\theta_{F,p},$
$\displaystyle\phi_{2}\left(t\right)=\phi_{2,\mathrm{diffuse}}\left(t\right)+g\theta_{F,p}.$
(10)

Here, the second terms on the right-hand side account for the effect of the
change in pump phase $\theta_{F,p}$ that is split between modes 1 and 2
according to a fixed ratio $g$ that was discussed in section III.1. Our goal
is to choose $\theta_{F,p+1}$ in the next cycle to maintain
$\phi_{2}\left(t\right)$ close to zero. As discussed in section II.2, there is
near perfect anti-correlation in the phase diffusion of modes 1 and 2 for
fixed $\theta_{F,p}$. In other words,
$\phi_{2,\mathrm{diffuse}}\left(t\right)$ can be inferred by measuring mode 1,
without measuring mode 2 itself. Rearranging the first equation in Eq. (IV.1)
yields:

$\displaystyle\phi_{2,\mathrm{diffuse}}\left(t\right)\
=-\phi_{1,\mathrm{diffuse}}\left(t\right)\
=\left(1-g\right)\theta_{F,p}-\phi_{1}\left(t\right).$ (11)

The measurement of $\phi_{1}\left(t\right)$, however, takes a finite amount of
time of data averaging to reduce the noise introduced by the detection
circuit. As shown in Fig. 4(b), at the beginning of each cycle, the system is
first allowed to settle in response to the change in $\theta_{F}$. Then the
measurement of $\phi_{1}$ is performed over time $t_{\mathrm{measure}}$ to
yield an averaged value ${\widetilde{\phi}}_{1,p}$. The pump phase
$\theta_{F,p+1}$ for the next cycle is chosen by replacing
$\phi_{1}\left(t\right)$ in Eq. (11) by ${\widetilde{\phi}}_{1,p}$ and then
substituting into the second equation in Eq. (IV.1) with
$\phi_{2}\left(t\right)$ set to zero, yielding:

$\displaystyle\theta_{F,p+1}=\frac{1}{g}[{\widetilde{\phi}}_{1,p}-\left(1-g\right)\theta_{F,p}].$
(12)

For stabilization of the phase of mode 2, the pump phase for the $(p+1)^{th}$
cycle is thus determined by both the measured phase of mode 1 and the pump
phase chosen for the $p^{th}$ cycle.

To keep $\phi_{2}(t)$ small we choose the time interval $\Delta t$ to be
short. Consequently, the increment
$\theta_{F,p+1}-\theta_{F,p}=(\tilde{\phi}_{1,p}-\theta_{F,p})/g$ is small. To
understand the evolution of the increments of $\theta_{F,p}$ we assume that
right before the $(p+1)$th correction, i.e., for $t=t_{p+1}-\delta$ with small
$\delta>0$, there are small deviations of $\phi_{2}$ from the target value of
zero, so that $\phi_{2}(t)=\varepsilon_{p}$, with $|\varepsilon_{p}|\ll 1$
being a random number. If we set
$\phi_{1}(t_{p+1}-\delta)=\tilde{\phi}_{1,p}$, i.e., if we disregard the
change of $\phi_{1}(t)$ over time $t_{\mathrm{measure}}$ and ignore the
uncertainty in measuring $\phi_{1}$, we have from Eq. (IV.1)
$\tilde{\phi}_{1,p}=\theta_{F,p}-\varepsilon_{p}$. The increment in pump phase
is then given by:

$\displaystyle\theta_{F,p+1}-\theta_{F,p}=-\varepsilon_{p}/g.$ (13)

With $g\sim 0.94$ in our experiment, $|\theta_{F,p+1}-\theta_{F,p}|$ is indeed
small. In a way one can think of this equation as a discrete analog of a
Langevin equation. The values of $\varepsilon_{p}$ are determined by the
imprecision of the measurements of the phase and of the control, and it is
reasonable to assume that $\varepsilon_{p}$ with different $p$ are
incorrelated. Then $\braket{\theta_{F,p}^{2}}$ increases with the number of
cycles as $p\braket{\varepsilon_{p}^{2}}/g^{2}$.

As discussed in section III.1, a step change in $\theta_{F}$ leads to changes
in both $\phi_{1}$ and $\phi_{2}$. Therefore the phase of mode 1 can also be
stabilized by measuring the phase of mode 2 and adjusting $\theta_{F}$ to
compensate for the phase diffusion. However, since
$\frac{{\Delta\phi}_{1}}{\delta\theta_{F}}=1-g<<1$, the required
$\delta\theta_{F}$ is much larger when compared to what is needed for
stabilization of mode 2, for comparable phase diffusion. The algorithm for
stabilizing mode 1 is shown in Fig. 5(a). It is similar to that for
stabilizing mode 2, except that mode 2 is measured instead and the factor $g$
in Eq. (12) and Eq. (13) is replaced with $1-g$:

$\displaystyle\theta_{F,p+1}=\frac{1}{1-g}\left[{\widetilde{\phi}}_{2,p}-g\theta_{F,p}\right]$
(14) $\displaystyle\theta_{F,p+1}-\theta_{F,p}=-\varepsilon_{p}/(1-g),$ (15)

where ${\widetilde{\phi}}_{2,p}$ is the measured value of $\phi_{2}$ at the
end of the $p^{th}$ cycle. In Eq. (14), $\theta_{F,p}$ is multiplied by $\sim$
16\. It is therefore essential to take into account the effect of the pump
phase change on mode 2 for the stabilization scheme to work properly for mode
1. By comparing Eq. (15) to Eq. (13), we find that the increment in pump phase
required for stabilizing mode 1 is much larger than that for stabilizing mode
2, which is consistent with the fact that mode 1 only takes up a small
fraction ($\sim 6\%$) of a step change of the pump phase.

Figure 5: (a) Algorithm for stabilizing $\phi_{1}$ by measuring $\phi_{2}$ and
adjusting $\theta_{F}$. (b) $\phi_{1}$ is measured as a function of time when
the stabilization routine is turned off (green) and on (dark green). The pump
current amplitude is 189 $\mu$A and frequency detune is 300 Hz. (c) The
spectra of the self-sustained vibrations of mode 1 with and without
stabilization. (d) Dependence of the standard deviation of $\phi_{1}$ on
$t_{\mathrm{measure}}$ at a fixed $t_{\mathrm{wait}}$ of 300 ms.

###  IV.2 Implementation

We demonstrate stabilization of the vibration frequencies for a pumping
current of 189 $\mu$A at detune frequency $\Delta_{F}=300$ Hz. Figure 4(c)
shows how the phase of mode 2 is stabilized by adjusting $\theta_{F}$ based on
measurement of $\phi_{1}.\ \phi_{1,2}$ are measured with two lockin amplifiers
referenced to $\omega_{\mathrm{self1,2}}$ respectively, with time constant of
1 ms. Measurements are recorded every 4.4 ms. $t_{\mathrm{wait}}$ is chosen as
300 ms that exceeds the settling time in response to a step change in
$\theta_{F}$. $t_{\mathrm{measure}}$ is chosen as 44 ms. Without
stabilization, $\phi_{2}$ undergoes diffusion (yellow). The phase fluctuations
determine the width of the spectral peak in the inset of Fig. 4b. When the
stabilization is turned on (brown), fluctuations in $\phi_{2}$ at time scales
exceeding $\sim 1$ s are significantly reduced. Spectral components close to
$\omega_{\mathrm{self2}}$ are suppressed to levels given by the noise in the
detection circuit. For example, the single side-band phase noise at offset of
0.1 Hz is reduced by 29 dB with phase stabilization turned on.

We optimize the stabilization algorithm with respect to
$t_{\mathrm{measure}}$. Figure 4(d) plots the standard deviation
$\sigma_{\phi_{2}}$ i.e., the root mean squared value of $\phi_{2}$, measured
as a function of $t_{\mathrm{measure}}$ for constant $t_{\mathrm{wait}}$ of
300 ms. As $t_{\mathrm{measure}}$ increases from the minimum value of 4.4 ms,
$\sigma_{\phi_{2}}$ initially decreases because a longer
$t_{\mathrm{measure}}$ reduces the uncertainty introduced by the detection
noise in measuring $\phi_{1}$ that is used to determine $\theta_{F}$ in the
subsequent stabilization cycle. However, with further increase of
$t_{\mathrm{measure}}$, $\sigma_{\phi_{2}}$ increases. This rise occurs
because the value $\tilde{\phi}_{1p}$, which is used to determine the control
signal $\theta_{F,p+1}$ and is obtained by averaging $\phi_{1}(t)$ over time
$t_{\mathrm{measure}}$, is increasingly different from the actual value
$\phi_{1}(t_{p+1})$ required for efficient control. This difference is due to
the phase diffusion. The optimal $t_{\mathrm{measure}}$ depends on the
detection noise and both the intensity and the correlation time of the noise
that leads to the phase diffusion. It should be shorter than this correlation
time. We find that in our experiment, when $t_{\mathrm{measure}}>200$ ms,
phase diffusion becomes the dominant effect and leads to the rise in
$\sigma_{\phi_{2}}$. The optimal measured $\sigma_{\phi_{2}}$ is $\sim 12$
mrad attained at $t_{\mathrm{measure}}$ $\sim 44$ ms.

Figure 5(b) shows the stabilization of the phase of mode 1 by adjusting
$\theta_{F}$ based on measurement of $\phi_{2}$ using Eq. (14). The spectral
width of mode 1 is also reduced significantly, as shown in Fig. 5(c). In Fig.
5(d), the optimal value of $\sigma_{\phi_{1}}$ is $\sim 9.5$ mrad attained at
$t_{\mathrm{measure}}\sim 88$ ms. A larger $t_{\mathrm{measure}}$ is needed
compared to the case of stabilizing mode 2, consistent with the fact that the
uncertainty in detecting $\phi_{2}$ is larger than that of $\phi_{1}$ in our
experiment.

##  V DISCUSSION AND OUTLOOK

Efficient stabilization of the phase of one mode relies on accurate
determination of the phase of the other mode. The effect of the noise
introduced by the detection circuits on the latter can be reduced by
increasing the amplitude of the self-sustained oscillations. One way to
achieve high amplitudes for the self-sustained vibrations is to choose large
pump detuning, as shown in Fig. 1(d) and 1(e). When the pump detuning is
chosen to exceed $\omega_{c}$, the zero-amplitude state becomes stable. For
parameters chosen in our experiment, we have not encountered problems with
stabilizing mode 2. However, in stabilizing mode 1, the change in the pump
phase required to compensate for the phase diffusion is much larger than the
phase diffusion itself [Eq. (14)]. We find that occasionally the pump phase
changes of up to 20∘ are required in our experiment. Such large change of the
pump phase could perturb the system to switch to the zero-amplitude state if
the pump detuning is large. To avoid this problem, it is necessary to limit
the pump detuning frequency to values smaller than $\omega_{c}$ so that the
zero amplitude state does not become stable.

In the stabilization algorithm, the pump phase is adjusted once per cycle. As
a result, phase diffusion at time scales shorter than the cycle period are not
compensated for. The duration of one stabilization cycle $\Delta t$ is the sum
of $t_{\mathrm{wait}}$ and $t_{\mathrm{measure}}$. In our experiment,
$t_{\mathrm{wait}}$ is chosen to barely exceed the settling time of the system
in response to in step change in the pump phase. As shown in Figs. 4(d) and
5(d), $t_{\mathrm{measure}}$ that gives optimal $\sigma_{\phi_{1,2}}$ is much
shorter than $t_{\mathrm{wait}}$. Further improvements to extend the
stabilization to shorter time scales would rely on tuning system parameters
(Appendix D) or design new resonators to achieve smaller settling time.

To summarize, we present an algorithm to stabilize the frequency of self-
sustained oscillations that are induced by dynamical backaction in a sideband-
driven micromechanical resonator. It is based on periodically shifting the
phase of the pump to compensate for phase diffusion. When the system undergoes
phase diffusion in the presence of weak noise, the phases of modes 1 and 2 add
up to a constant (within our detection limit) that can be adjusted by the
phase of the sideband pump. For a step change of the pump phase, the phases of
both modes settle to new values after a transient. The duration of the
transient is determined by the smallest real part of the eigenvalues of the
linearized equations of motion about the vibrational state. If the pump change
is small, the phase change of each mode is proportional to the pump phase
change. The two proportionality constants add up to one. This finding,
together with the phase anti-correlation of the two modes, allow us to
stabilize the phase of one mode by measuring the phase of the other mode and
then compensating for the phase diffusion by adjusting the phase of the pump.
We demonstrate that phase fluctuations of either the high or low frequency
mode can be significantly reduced, resulting in a much narrower spectral
linewidth. The control algorithm requires knowledge of the phase shifts of the
modes in response to a shift in the pump phase. After the transient, the
proportionality constant between the mode phase shifts and the pump phase
shift is independent of the parameters of the pump. It is fully determined by
the parameters of the resonators.

The results open new opportunities in generating stable mechanical vibrations
via parametric down-conversion in nonlinear micromechanical resonators, with
controlled amplitude and phase. We reiterate that our scheme is distinct from
direct feedback that involves stabilizing a particular mode by measuring its
phase because a frequency reference near this mode is not required. As long as
a frequency reference is available to determine the phase diffusion of one of
the modes, the phase of the other mode can be stabilized even if its
eigenfrequency is orders of magnitude different from that of the former. In
principle, one can design resonators with the upper mode reaching microwave
frequencies in which the construction of signal sources with low phase noise
and small long term drift requires complex circuitry. The analysis also can be
extended to mechanical resonators coupled to optical or microwave cavities
[50, 51].

Acknowledgements: This work is supported by the Research Grants Council of
Hong Kong SAR, China (Project No. 16304219). M.I.D. acknowledges partial
support from the US Defense Advanced Research Projects Agency (Grant No.
HR0011-23-2-004) and from the Moore Foundation (Grant No. 12214)

##  Appendix A DISPERSIVE COUPLING BETWEEN MODES 1 AND 2

Modes 1 and 2 in the resonator are dispersively coupled, with interaction
energy of the form $\frac{1}{2}\gamma q_{1}^{2}q_{2}^{2}$. Figure 6 shows that
the decrease in the resonant frequency of one mode is proportional to the
square vibration amplitude of the other mode. In this measurement, the
sideband pumping is removed. Instead, ac voltages at frequencies close to
$\omega_{1,2}$ are applied to the side gates of modes 1 and 2 respectively to
generate periodic electrostatic forces to excite forced vibrations at
amplitudes $A_{1,2}$. From the slopes of the linear fits, $\gamma$ is
determined to be $6.41\times 10^{12}\ kg\ rad^{2}s^{-2}m^{-2}$.

Figure 6: (a) Scaled shifts of the frequency of mode 1 versus the square of
the amplitude $A_{2}$ of forced vibrations in mode 2. The line is a linear fit
through the origin. Inset: Spectra of mode 1 under a periodic electrostatic
force. Each resonance peak corresponds to a data point in the main panel. (b)
A similar plot for mode 2.

##  Appendix B EIGENVALUES FOR THE LINEARIZED EQUATIONS OF MOTION NEAR THE
SUPERCRITICAL BIFURCATION POINT

The nonlinearity parameters in the equation for the phases of self-sustained
vibrations (II.2) have the form

$\displaystyle\gamma^{(\pm)}_{ij}=\gamma_{ij}\pm\frac{3\gamma_{jj}\omega_{i}}{2\omega_{j}}\quad(i,j=1,2;i\neq
j)$ (16)

By setting $\dot{r}_{1,2}=\dot{\phi}_{\pm}=0$ in Eqs. (II.2), we obtain the
ratio of the stationary amplitudes $r_{1,2}^{(0)}$ and the phase
$\phi_{+}^{(0)}\equiv\phi_{1}^{(0)}+\phi_{2}^{(0)}$, cf. Ref. [47]:

$\displaystyle\sin\Theta^{(0)}=1/\Xi,\quad\Theta^{(0)}=\theta_{F}-\phi_{+}^{(0)},$
$\displaystyle\zeta^{(0)}\equiv
r_{2}^{(0)}/r_{1}^{(0)}=(\Gamma_{1}F_{2}\omega_{1}/\Gamma_{2}F_{1}\omega_{2})^{1/2},$
(17)

where $\Xi$ is given by Eq. (II.1).

The vibration frequency in the rotating frame $\delta\omega$ and the squared
radius of the vibrations of mode 1 are given by the expressions

$\displaystyle\delta\omega=\left(\frac{\gamma_{12}}{\omega_{1}}\zeta^{(0)}{}^{2}+\frac{3\gamma_{11}}{2\omega_{1}}\right){r_{1}^{(0)}}^{2}-\Gamma_{1}\left(\Xi^{2}-1\right)^{\frac{1}{2}}$
$\displaystyle{r_{1}^{(0)}}^{2}=\frac{\Gamma_{2}\omega_{2}F_{1}}{\gamma^{(+)}_{21}\Gamma_{2}F_{1}+\gamma^{(+)}_{12}\Gamma_{1}F_{2}}\left(\Delta_{F}+\omega_{c}\right),$
(18)

where the critical detuning $\omega_{c}$ is given by Eq. (II.1),
$\omega_{c}=(\Gamma_{1}+\Gamma_{2})(\Xi^{2}-1)^{1/2}$.

The linearized equations of motion (8) for the increments

$x_{1}=\delta\zeta,\quad x_{2}=\theta_{F}-\delta\phi_{+}-\Theta^{(0)},\quad
x_{3}=\delta r_{1}/r_{1}^{(0)}$

of the dynamical variables $\zeta=r_{2}/r_{1},\phi_{+}=\phi_{1}+\phi_{2}$, and
$r_{1}$ follow from Eq. (II.2). We will write them in a vector form for the
vector ${\bf x}$ with components $x_{1},x_{2},x_{3}$. We emphasize that the
unknown quantity is $\delta\phi_{+}$. It is a dynamical variable, whereas
$\theta_{F}$ is the control parameter. Therefore
$\dot{x}_{2}\equiv-d\delta\phi_{+}/dt$. The equations read

$\displaystyle\dot{\bf x}=\hat{\Lambda}{\bf
x};\quad\Lambda_{11}=-\left(\Gamma_{1}+\Gamma_{2}\right),\quad\Lambda_{12}=\left(\Gamma_{2}-\Gamma_{1}\right)\zeta^{\left(0\right)}\left(\Xi^{2}-1\right)^{1/2},\quad\Lambda_{13}=0$
$\displaystyle\Lambda_{21}=-\left(\Gamma_{2}-\Gamma_{1}\right)\left(\Xi^{2}-1\right)^{1/2}/\zeta^{(0)}-\left(2\gamma_{12}^{+}/\omega_{1}\right){r_{1}^{\left(0\right)}}^{2}\zeta^{\left(0\right)},\quad\Lambda_{22}=-(\Gamma_{1}+\Gamma_{2}),\quad\Lambda_{23}=-2[\left(\gamma_{12}^{+}/\omega_{1}\right){\zeta^{\left(0\right)}}^{2}+(\gamma_{21}^{+}/\omega_{2})]r_{1}^{\left(0\right)}{}^{2},$
$\displaystyle\Lambda_{31}=\left(\Gamma_{1}\Gamma_{2}F_{1}\omega_{2}/F_{2}\omega_{1}\right)^{1/2},\quad\Lambda_{32}=\Gamma_{1}\left(\Xi^{2}-1\right)^{1/2},\quad\Lambda_{33}=0.$
(19)

The analysis of the dynamics is simplified near the supercritical bifurcation
point $\Delta_{F}=-\omega_{c}$ where the state of self-sustained vibrations
emerges. Here $r_{1,2}^{(0)}<<1$, as seen from Eq. (B), whereas the ratio
$\zeta^{(0)}$ is not small. Therefore, as a first step, we can set
$r_{1}^{(0)}=0$ in Eqs. (B). This leads to a system of two coupled equations
for $\delta\zeta$ and $\delta\phi_{+}$. The corresponding eigenvalues are

$\displaystyle\lambda_{1,2}\approx-\left(\Gamma_{1}+\Gamma_{2}\right)\pm
i\omega_{c}(\Gamma_{2}-\Gamma_{1})/(\Gamma_{2}+\ \Gamma_{1}).$ (20)

Both the real and the imaginary parts of $\lambda_{1,2}$ are finite at the
bifurcation point, as shown in Figs. 3(e) and (f).

The third eigenvalue of the matrix $\hat{\Lambda}$ is small,
$|\lambda_{3}|\ll\Gamma_{1,2}$ near the bifurcation point. It describes the
slow decay of $x_{3}\equiv\delta r_{1}$. The variables $x_{1},x_{2}$ follow
this decay adiabatically. They can be found from Eq. (B) by setting
$\dot{x}_{1}=\dot{x}_{2}=0$. Then both $x_{1}$ and $x_{2}$ are proportional to
$\Lambda_{23}\propto r_{1}^{(0)}{}^{2}x_{3}$. The equation for $x_{3}$ then
takes the form $\dot{x}_{3}=\lambda_{3}x_{3}$ with $\lambda_{3}\propto-
r_{1}^{(0)}{}^{2}\propto\Delta_{F}+\omega_{c}$.

##  Appendix C INCREMENT OF THE PHASE DIFFERENCE

In this section, we find the values of the phases $\phi_{1},\phi_{2}$ after
the transient process that follows incrementing the driving field phase
$\theta_{F}$ by $\delta\theta_{F}$. This can be done by diagonalizing the
matrix $\hat{\Lambda}$. Respectively, we introduce right and left eigenvectors
of this matrix, ${\bf y}_{\nu}$ and $\bar{\bf y}_{\mu}$,

$\hat{\Lambda}{\bf y}_{\nu}=\lambda_{\nu}{\bf y}_{\nu},\quad\bar{\bf
y}_{\nu}\hat{\Lambda}=\lambda_{\nu}\bar{\bf y}_{\nu},\quad\bar{\bf
y}_{\nu^{\prime}}\cdot{\bf y}_{\nu}=\delta_{\nu\nu^{\prime}}.$

We then expand ${\bf x}(t)$ in ${\bf y}$ and use Eq. (B) to obtain

$\displaystyle{\bf x}(t)=\sum_{\nu}C_{\nu}(t){\bf y}_{\nu},\quad
C_{\nu}(t)=\exp(\lambda_{\nu}t)C_{\nu}(0),$ $\displaystyle C_{\nu}(0)=\bar{\bf
y}_{\nu}\cdot{\bf x}(0)$ (21)

To fully describe time evolution of the coupled oscillators in response to a
perturbation, and in particular in response to the increment of the phase of
the drive $\theta_{F}$, this equation has to be complemented by the expression
for the increment of the phase $\delta\phi_{-}$. The equation for
$\delta\phi_{-}$ can be obtained by linearizing Eq. (6), which gives

$\displaystyle\frac{d}{dt}\delta\phi_{-}={\bm{\Phi}}\cdot{\bf
x},\quad\Phi_{1}=-(\Gamma_{1}+\Gamma_{2})(\Xi^{2}-1)^{1/2}/\zeta^{(0)}$
$\displaystyle+2(\gamma_{12}^{-}/\omega_{1})r_{1}^{(0)}{}^{2}\zeta^{(0)},\quad\Phi_{2}=\Gamma_{1}-\Gamma_{2},$
$\displaystyle\Phi_{3}=2[(\gamma_{12}^{-}/\omega_{1})\zeta^{(0)}{}^{2}-(\gamma_{21}^{-}/\omega_{2})]r_{1}^{(0)}{}^{2}.$
(22)

The solution of Eqs. (C) and (C) was used in Fig. 3 to illustrate the
evolution of the parameters in time.

For $t\gg t_{r}$, where $t_{r}=\max|\mathrm{Re}\,\lambda_{1,2,3}|^{-1}$, all
components of the vector ${\bf x}(t)$ go to zero, and therefore
$\delta\phi_{+}(t)\to\delta\theta_{F}$. To find the change of $\phi_{-}$ we
integrate Eq. (C) over time, which gives in the limit of large time

$\displaystyle\delta\phi_{-}(t)\to-\sum_{\nu}\bigl{(}\bar{\bf
y}_{\nu}\cdot{\bf x}(0)\bigr{)}\bigl{(}{\bm{\Phi}}\cdot{\bf
y}_{\nu}\bigr{)}/\lambda_{\nu}.$ (23)

For ${\bf x}(0)=(0,\delta\theta_{F},0)$ this gives $\delta\phi_{-}(t)\to
C\delta\theta_{F}$, with

$\displaystyle C=-\sum_{\nu}(\bar{\bf
y}_{\nu})_{2}\bigl{(}{\bm{\Phi}}\cdot{\bf y}_{\nu}\bigr{)}/\lambda_{\nu}.$
(24)

This expression has been used in Eq. (9) to find the values approached by the
phases $\phi_{1}$ and $\phi_{2}$ after a change of the phase of the drive. We
note that, even though two of the three eigenvalues $\lambda_{\nu}$ are
complex, it is clear that $C$ is real.

###  C.1 Adiabatic change of the drive frequency

An insight into the value of the parameter $C$ can be gained by studying the
dynamics of the system for slowly varying frequency of the driving field,
where $\Delta_{F}$ slowly increases and then decreases back to its value,
$\Delta_{F}\to\Delta_{F}+\delta\Delta_{F}(t)$. We assume that
$\delta\Delta_{F}(t)$ is small and
$|d\delta\Delta_{F}(t)/dt|\ll|\delta\Delta_{F}(t)|/t_{r}$. In this case, as
seen from Eq. (II.2), the equation of motion for the increments ${\bf x}$
takes the form

$\dot{x}_{i}=\sum_{j}\Lambda_{ij}x_{j}-\delta\Delta_{F}(t)\,\delta_{i,2},\quad{\bf
x}(0)={\bf 0}.$

Changing to variables ${\bf y}$, as in Eq. (C), we obtain for the expansion
coefficients $C_{\nu}$

$\displaystyle\dot{C}_{\nu}=\lambda_{\nu}C_{\nu}-\delta\Delta_{F}(t)(\bar{\bf
y}_{\nu})_{2}.$ (25)

where $(\bar{\bf y}_{\nu})_{2}$ is the second component of the vector
$\bar{\bf y}_{\nu}$.

Assuming that $\delta\Delta_{F}(t)$ varies slowly, we can disregard
$\dot{C}_{\nu}$ in this equation, which gives

$C_{\nu}=\lambda_{\nu}^{-1}(\bar{\bf y}_{\nu})_{2}\delta\Delta_{F}(t).$

After a slow pulse $\delta\Delta_{F}(t)$ the coefficients $C_{\nu}$ become
equal to zero.

We now consider Eq. (6) for $\delta\phi_{-}$. It has the form

$\displaystyle\frac{d}{dt}\delta\phi_{-}=\delta\Delta_{F}(t)\sum_{\nu}(\bar{\bf
y}_{\nu})_{2}\bigl{(}{\bm{\Phi}}\cdot{\bf y}_{\nu}\bigr{)}/\lambda_{\nu}$
$\displaystyle+\delta\Delta_{F}(t)\left[1-2\frac{d}{d\Delta_{F}}\delta\omega\right].$
(26)

To calculate the derivative of $\delta\omega$ in the adiabatic limit we use
Eq. (B). As seen from this equation, $\delta\omega$ is linear in $\Delta_{F}$.
Since $\zeta^{(0}$ is independent of $\Delta_{F}$, we obtain

$\displaystyle\frac{d}{d\Delta_{F}}\delta\omega=\left(\frac{\gamma_{12}}{\omega_{1}}\zeta^{(0)}{}^{2}+\frac{3\gamma_{11}}{2\omega_{1}}\right)\frac{\Gamma_{2}\omega_{2}F_{1}}{\gamma^{(+)}_{21}\Gamma_{2}F_{1}+\gamma^{(+)}_{12}\Gamma_{1}F_{2}}.$
(27)

In the adiabatic limit the value of $\phi_{-}$ does not change. Therefore the
coefficient of $\delta\Delta_{F}(t)$ in Eq. (C.1) is equal to zero. Since the
first term in the right-hand side of this equation is $-C\delta\Delta_{F}(t)$,
the condition $\dot{\phi}_{-}=0$ provides a compact form of the parameter $C$
in terms of $d\delta\omega/d\Delta_{F}$,

$\displaystyle C=1-2\frac{d\delta\omega}{d\Delta_{F}}.$ (28)

An important consequence is that $C$ is independent of the value of
$\Delta_{F}$ itself. Moreover, given that $F_{1}\propto F_{2}$, we see that
$C$ is independent of the drive amplitude as well.

##  Appendix D DEPENDENCE OF THE SETTLING TIME ON PUMP CURRENT

The time required for the two modes to settle in response to a step change in
$\theta_{F}$ is determined by relaxation time $|\lambda_{3}^{-1}|$. In terms
of stabilizing the frequencies of self-sustained vibrations, it is desirable
to have large $|\lambda_{3}|$ so that the duration of each cycle in the
stabilization algorithm can be reduced and $\theta_{F}$ can be updated more
frequently to compensate for phase diffusion. Figure 7 plots the dependence of
$\lambda_{3}$ on the amplitude of the pump current when the pump detune
frequency is fixed at 0 Hz. Initially, $\lambda_{3}$ becomes more negative
rapidly as the pump current amplitude increases and then levels off. For the
data reported in the main text, the pump current amplitude is chosen to be 189
$\mu$A. Further increase in the pump current amplitude only marginally shorten
the relaxation time, while significantly increase the risk of damaging the
thin gold layer on the suspended beams.

Figure 7: Dependence of eigenvalue $\lambda_{3}$ of matrix $\mathbf{K}$ in the
equation of motion linearized about a stable self-sustained vibrations state
on the amplitude of the pump current. The pump current frequency detune is 0
Hz.

## References

  * Jenkins [2013] A. Jenkins, Self-oscillation, Physics Reports 525, 167 (2013).
  * Kim  _et al._ [2007] B. Kim, R. N. Candler, M. A. Hopcroft, M. Agarwal, W.-T. Park, and T. W. Kenny, Frequency stability of wafer-scale film encapsulated silicon based MEMS resonators, Sensors and Actuators A: Physical 136, 125 (2007).
  * van Beek and Puers [2012] J. T. M. van Beek and R. Puers, A review of MEMS oscillators for frequency reference and timing applications, Journal of Micromechanics and Microengineering 22, 013001 (2012).
  * Zaliasl  _et al._ [2015] S. Zaliasl, J. C. Salvia, G. C. Hill, L. W. Chen, K. Joo, R. Palwai, N. Arumugam, M. Phadke, S. Mukherjee, H.-C. Lee, C. Grosjean, P. M. Hagelin, S. Pamarti, T. S. Fiez, K. A. A. Makinwa, A. Partridge, and V. Menon, A 3 ppm 1.5 × 0.8 mm 2 1.0 µA 32.768 kHz MEMS-Based Oscillator, IEEE Journal of Solid-State Circuits 50, 291 (2015).
  * Roshan  _et al._ [2016] M. H. Roshan, S. Zaliasl, K. Joo, K. Souri, R. Palwai, W. Chen, S. Pamarti, J. C. Doll, N. Miller, C. Arft, S. Tabatabaei, C. Sechen, A. Partridge, and V. Menon, 11.1 Dual-MEMS-resonator temperature-to-digital converter with 40 K resolution and FOM of 0.12pJK2, in  _2016 IEEE International Solid-State Circuits Conference (ISSCC)_ (IEEE, San Francisco, CA, USA, 2016) pp. 200–201.
  * Eichler  _et al._ [2013] A. Eichler, J. Moser, M. Dykman, and A. Bachtold, Symmetry breaking in a mechanical resonator made from a carbon nanotube, Nature Communications 4, 2843 (2013).
  * Cleland and Roukes [1998] A. N. Cleland and M. L. Roukes, A nanometre-scale mechanical electrometer, Nature 392, 160 (1998).
  * Rugar  _et al._ [2004] D. Rugar, R. Budakian, H. J. Mamin, and B. W. Chui, Single spin detection by magnetic resonance force microscopy, Nature 430, 329 (2004).
  * Hanay  _et al._ [2012] M. S. Hanay, S. Kelber, A. K. Naik, D. Chi, S. Hentz, E. C. Bullard, E. Colinet, L. Duraffourg, and M. L. Roukes, Single-protein nanomechanical mass spectrometry in real time, Nature Nanotechnology 7, 602 (2012).
  * Albrecht  _et al._ [1991] T. R. Albrecht, P. Grütter, D. Horne, and D. Rugar, Frequency Modulation Detection Using High Q Cantilevers for Enhanced Force Microscope Sensitivity, J. Appl. Phys. 69, 668 (1991).
  * Vig and Kim [1999] J. R. Vig and Y. Kim, Noise in Microelectromechanical System Resonators, IEEE Trans Ultrason Ferroelectr Freq Control 46, 1558 (1999).
  * Yurke  _et al._ [1995a] B. Yurke, D. S. Greywall, A. N. Pargellis, and P. A. Busch, Theory of amplifier-noise evasion in an oscillator employing a nonlinear resonator, Phys. Rev. A 51, 4211 (1995a).
  * Feng  _et al._ [2008] X. L. Feng, C. J. White, A. Hajimiri, and M. L. Roukes, A self-sustaining ultrahigh-frequency nanoelectromechanical oscillator, Nature Nanotechnology 3, 342 (2008).
  * Kenig  _et al._ [2012a] E. Kenig, M. C. Cross, L. G. Villanueva, R. B. Karabalin, M. H. Matheny, R. Lifshitz, and M. L. Roukes, Optimal operating points of oscillators using nonlinear resonators, Physical Review E 86, 056207 (2012a).
  * Villanueva  _et al._ [2011] L. G. Villanueva, R. B. Karabalin, M. H. Matheny, E. Kenig, M. C. Cross, and M. L. Roukes, A Nanoscale Parametric Feedback Oscillator, Nano Letters 11, 5054 (2011).
  * Antonio  _et al._ [2012a] D. Antonio, D. H. Zanette, and D. López, Frequency stabilization in nonlinear micromechanical oscillators, Nature Communications 3, 10.1038/ncomms1813 (2012a).
  * Chen  _et al._ [2016] C. Chen, D. H. Zanette, J. R. Guest, D. A. Czaplewski, and D. López, Self-Sustained Micromechanical Oscillator with Linear Feedback, Physical Review Letters 117, 017203 (2016).
  * Huang  _et al._ [2019a] L. Huang, S. M. Soskin, I. A. Khovanov, R. Mannella, K. Ninios, and H. B. Chan, Frequency stabilization and noise-induced spectral narrowing in resonators with zero dispersion, Nature Communications 10, 10.1038/s41467-019-11946-8 (2019a).
  * Defoort  _et al._ [2021] M. Defoort, S. Hentz, S. W. Shaw, and O. Shoshani, Amplitude stabilization in a synchronized nonlinear nanomechanical oscillator, ArXiv211208788 Cond-Mat Physicsphysics (2021), arxiv:2112.08788 [cond-mat, physics:physics] .
  * Velikovich  _et al._ [1991] A. Velikovich, G. Golubev, V. Golubchenko, and D. Luchinsky, Tristability and self-oscillations in a double resonator system, Optics Communications 80, 444 (1991).
  * Zalalutdinov  _et al._ [2001] M. Zalalutdinov, A. Zehnder, A. Olkhovets, S. Turner, L. Sekaric, B. Ilic, D. Czaplewski, J. M. Parpia, and H. G. Craighead, Autoparametric optical drive for micromechanical oscillators, Applied Physics Letters 79, 695 (2001).
  * Barton  _et al._ [2012] R. A. Barton, I. R. Storch, V. P. Adiga, R. Sakakibara, B. R. Cipriany, B. Ilic, S. P. Wang, P. Ong, P. L. McEuen, J. M. Parpia, and H. G. Craighead, Photothermal Self-Oscillation and Laser Cooling of Graphene Optomechanical Systems, Nano Lett. 12, 4681 (2012).
  * De Alba  _et al._ [2017] R. De Alba, T. S. Abhilash, R. H. Rand, H. G. Craighead, and J. M. Parpia, Low-Power Photothermal Self-Oscillation of Bimetallic Nanowires, Nano Letters 17, 3995 (2017).
  * Roxworthy and Aksyuk [2018] B. J. Roxworthy and V. A. Aksyuk, Electrically tunable plasmomechanical oscillators for localized modulation, transduction, and amplification, Optica 5, 71 (2018).
  * Ochs  _et al._ [2022] J. S. Ochs, D. K. J. Boneß, G. Rastelli, M. Seitner, W. Belzig, M. I. Dykman, and E. M. Weig, Frequency Comb from a Single Driven Nonlinear Nanomechanical Mode, Physical Review X 12, 041019 (2022).
  * Etaki  _et al._ [2013] S. Etaki, F. Konschelle, Ya. M. Blanter, H. Yamaguchi, and H. S. J. van der Zant, Self-sustained oscillations of a torsional SQUID resonator induced by Lorentz-force back-action, Nature Communications 4, 1803 (2013).
  * Schmid  _et al._ [2016] S. Schmid, L. G. Villanueva, and M. L. Roukes,  _Fundamentals of Nanomechanical Resonators_ (Springer, Switzerland, 2016).
  * Demir and Hanay [2020] A. Demir and M. S. Hanay, Fundamental Sensitivity Limitations of Nanomechanical Resonant Sensors Due to Thermomechanical Noise, IEEE Sens. J. 20, 1947 (2020).
  * Bachtold  _et al._ [2022] A. Bachtold, J. Moser, and M. I. Dykman, Mesoscopic physics of nanomechanical systems, Rev. Mod. Phys. 94, 045005 (2022).
  * Rytov [1956a] S. M. Rytov, Fluctuations in Oscillating Systems of the Thompson Type .1., Sov. Phys. JETP-USSR 2, 217 (1956a).
  * Rytov [1956b] S. M. Rytov, Fluctuations in Oscillating Systems of the Thompson Type .2., Sov. Phys. JETP-USSR 2, 225 (1956b).
  * Rubiola [2009] E. Rubiola,  _Phase Noise and Frequency Stability in Oscillators_ (Cambridge University Press, Cambridge, 2009).
  * Huang  _et al._ [2019b] L. Huang, S. M. Soskin, I. A. Khovanov, R. Mannella, K. Ninios, and H. B. Chan, Frequency stabilization and noise-induced spectral narrowing in resonators with zero dispersion, Nat. Commun. 10, 3930 (2019b).
  * Miller  _et al._ [2021] N. J. Miller, S. W. Shaw, and M. I. Dykman, Suppressing Frequency Fluctuations of Self-Sustained Vibrations in Underdamped Nonlinear Resonators, Phys. Rev. Applied 15, 014024 (2021).
  * Antonio  _et al._ [2012b] D. Antonio, D. H. Zanette, and D. Lopez, Frequency Stabilization in Nonlinear Micromechanical Oscillators, Nat. Commun. 3, 806 (2012b).
  * Zhao  _et al._ [2017] C. Zhao, G. Sobreviela, M. Pandit, S. Du, X. Zou, and A. Seshia, Experimental Observation of Noise Reduction in Weakly Coupled Nonlinear MEMS Resonators, J. Microelectromechanical Syst. 26, 1196 (2017).
  * Greywall  _et al._ [1994] D. S. Greywall, B. Yurke, P. A. Busch, A. N. Pargellis, and R. L. Willett, Evading Amplifier Noise in Nonlinear Oscillators, Phys. Rev. Lett. 72, 2992 (1994).
  * Yurke  _et al._ [1995b] B. Yurke, D. S. Greywall, A. N. Pargellis, and P. A. Busch, Theory of Amplifier-Noise Evasion in an Oscillator Employing a Nonlinear Resonator, Phys. Rev. A 51, 4211 (1995b).
  * Kenig  _et al._ [2012b] E. Kenig, M. C. Cross, L. G. Villanueva, R. B. Karabalin, M. H. Matheny, R. Lifshitz, and M. L. Roukes, Optimal Operating Points of Oscillators Using Nonlinear Resonators, Phys. Rev. E 86, 056207 (2012b).
  * Kenig  _et al._ [2013] E. Kenig, M. C. Cross, J. Moehlis, and K. Wiesenfeld, Phase Noise of Oscillators with Unsaturated Amplifiers, Phys. Rev. E 88, 062922 (2013).
  * Villanueva  _et al._ [2013] L. G. Villanueva, E. Kenig, R. B. Karabalin, M. H. Matheny, R. Lifshitz, M. C. Cross, and M. L. Roukes, Surpassing Fundamental Limits of Oscillators Using Nonlinear Resonators, Phys. Rev. Lett. 110, 177208 (2013).
  * Kenig  _et al._ [2012c] E. Kenig, M. C. Cross, R. Lifshitz, R. B. Karabalin, L. G. Villanueva, M. H. Matheny, and M. L. Roukes, Passive Phase Noise Cancellation Scheme, Phys. Rev. Lett. 108, 264102 (2012c).
  * Berstein [1938] I. L. Berstein, On fluctuations in the neighborhood of periodic motion of an auto-oscillating system, Sov Acad Sci Dokl. 20, 11 (1938).
  * Haken [1966] H. Haken, Theory of Intensity and Phase Fluctuations of a Homogeneously Broadened Laser, Z. Phys. 190, 327 (1966).
  * Lax [1967] M. Lax, Classical Noise .v. Noise in Self Sustained Oscillators, Phys. Rev. 160, 290 (1967).
  * Agrawal and Seshia [2014] D. K. Agrawal and A. A. Seshia, An analytical formulation for phase noise in MEMS oscillators, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 61, 1938 (2014).
  * Sun  _et al._ [2016] F. Sun, X. Dong, J. Zou, M. I. Dykman, and H. B. Chan, Correlated anomalous phase diffusion of coupled phononic modes in a sideband-driven resonator, Nature Communications 7, 10.1038/ncomms12694 (2016).
  * Demir  _et al._ [2000] A. Demir, A. Mehrotra, and J. Roychowdhury, Phase noise in oscillators: A unifying theory and numerical methods for characterization, IEEE Trans. Circuits Syst. Fundam. Theory Appl. 47, 655 (2000).
  * Mahboob  _et al._ [2012] I. Mahboob, K. Nishiguchi, H. Okamoto, and H. Yamaguchi, Phonon-cavity electromechanics, Nature Physics 8, 387 (2012).
  * Kippenberg and Vahala [2008] T. J. Kippenberg and K. J. Vahala, Cavity Optomechanics: Back-Action at the Mesoscale, Science 321, 1172 (2008).
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Reviews of Modern Physics 86, 1391 (2014).
  * Dykman [1978] M. Dykman, Heating and cooling of local and quasilocal vibrations by a nonresonance field, Sov. Phys. Solid State 20, 1306 (1978).
  * Marquardt  _et al._ [2006] F. Marquardt, J. G. E. Harris, and S. M. Girvin, Dynamical Multistability Induced by Radiation Pressure in High-Finesse Micromechanical Optical Cavities, Physical Review Letters 96, 103901 (2006).
  * Kippenberg  _et al._ [2005] T. J. Kippenberg, H. Rokhsari, T. Carmon, A. Scherer, and K. J. Vahala, Analysis of Radiation-Pressure Induced Mechanical Oscillation of an Optical Microcavity, Physical Review Letters 95, 033901 (2005).
  * Carmon  _et al._ [2005] T. Carmon, H. Rokhsari, L. Yang, T. J. Kippenberg, and K. J. Vahala, Temporal Behavior of Radiation-Pressure-Induced Vibrations of an Optical Microcavity Phonon Mode, Physical Review Letters 94, 223902 (2005).
  * Metzger  _et al._ [2008] C. Metzger, M. Ludwig, C. Neuenhahn, A. Ortlieb, I. Favero, K. Karrai, and F. Marquardt, Self-Induced Oscillations in an Optomechanical System Driven by Bolometric Backaction, Physical Review Letters 101, 133903 (2008).
  * Bagheri  _et al._ [2011] M. Bagheri, M. Poot, M. Li, W. P. H. Pernice, and H. X. Tang, Dynamic manipulation of nanomechanical resonators in the high-amplitude regime and non-volatile mechanical memory operation, Nature Nanotechnology 6, 726 (2011).
  * Faust  _et al._ [2012] T. Faust, P. Krenn, S. Manus, J. Kotthaus, and E. Weig, Microwave cavity-enhanced transduction for plug and play nanomechanics at room temperature, Nature Communications 3, 728 (2012).
  * Buters  _et al._ [2015] F. M. Buters, H. J. Eerkens, K. Heeck, M. J. Weaver, B. Pepper, S. de Man, and D. Bouwmeester, Experimental exploration of the optomechanical attractor diagram and its dynamics, Physical Review A 92, 013811 (2015).