# Mechanical overtone frequency combs

Matthijs H. J. de Jong  Department of Precision and Microsystems Engineering,
Delft University of Technology, Mekelweg 2, 2628CD Delft, The Netherlands
Kavli Institute of Nanoscience, Department of Quantum Nanoscience, Delft
University of Technology, Lorentzweg 1, 2628CJ Delft, The Netherlands  Adarsh
Ganesan  Ahmedabad University, Ahmedabad, Gujarat 380009, India  National
Institute of Standards and Technology, Gaithersburg, Maryland 20899, USA
Andrea Cupertino  Department of Precision and Microsystems Engineering, Delft
University of Technology, Mekelweg 2, 2628CD Delft, The Netherlands  Simon
Gröblacher  Kavli Institute of Nanoscience, Department of Quantum Nanoscience,
Delft University of Technology, Lorentzweg 1, 2628CJ Delft, The Netherlands
Richard A. Norte  Department of Precision and Microsystems Engineering, Delft
University of Technology, Mekelweg 2, 2628CD Delft, The Netherlands  Kavli
Institute of Nanoscience, Department of Quantum Nanoscience, Delft University
of Technology, Lorentzweg 1, 2628CJ Delft, The Netherlands
<EMAIL_ADDRESS>
###### Abstract

Mechanical frequency combs are poised to bring the applications and utility of
optical frequency combs into the mechanical domain. So far, their main
challenge has been strict requirements on drive frequencies and power, which
complicate operation. We demonstrate a straightforward mechanism to create a
frequency comb consisting of mechanical overtones (integer multiples) of a
single eigenfrequency, by monolithically integrating a suspended dielectric
membrane with a counter-propagating optical trap. The periodic optical field
modulates the dielectrophoretic force on the membrane at the overtones of a
membrane’s motion. These overtones share a fixed frequency and phase relation,
and constitute a mechanical frequency comb. The periodic optical field also
creates an optothermal parametric drive that requires no additional power or
external frequency reference. This combination of effects results in an easy-
to-use mechanical frequency comb platform that requires no precise alignment,
no additional feedback or control electronics, and only uses a single,
$\mathrm{mW}$ continuous wave laser beam. This highlights the overtone
frequency comb as the straightforward future for applications in sensing,
metrology and quantum acoustics.

##  I Introduction

Over the last quarter century, optical frequency combs have become key tools
for metrology, timing and spectroscopy [1, 2], and are indispensable in many
laboratories around the world. The fixed frequency and phase relations between
the many different tones of a comb have revolutionized fields as astronomy [3]
or cosmology [4], and allowed tests of fundamental physics with atomic clocks
[5, 4]. Recent developments in optomechanics have expanded methods to create
optical frequency combs by interactions with mechanical resonators [6, 7, 8,
9, 10]. But in the last few years, a new paradigm of frequency combs has
appeared that are completely mechanical in nature: Phononic frequency combs
[11, 12, 13], also called acoustic or mechanical frequency combs. The ideas
behind this field began in nonlinear dynamics, where it was realized that
mixing in coupled oscillators may lead to a series of sidebands [14], that can
be regarded as a frequency comb if there exists a fixed phase relation [12]
between these sidebands. Experimental demonstrations have shown mechanical
frequency combs exist in different mechanical systems [15, 16, 17, 18, 19, 20,
21, 22, 23, 10], and have explored connections to well-known concepts in
nonlinear dynamics such as bifurcations [15, 18, 24], 3- or 4-wave mixing [14,
17, 25, 26], and symmetry-breaking [27]. The fixed frequency and phase
relation between the comb teeth allows the application of techniques known
from optical combs. These can e.g. improve position sensing accuracy in
optically opaque materials [28] such as underwater or medical imaging.
Mechanical frequency combs can further be used to track and stabilize
mechanical resonances [29, 30], and enhance Brillouin microscopy [31, 32].
Recent proposals from the field of quantum acoustics [33] foresee an important
role for devices with multiple mechanical modes with equal frequency spacing.
These may also be useful in scaling up transduction between optical and
microwave frequencies [34], and engineering the interference between
subsequent modes may aid in coupling phonons to multiple qubits [35]. Until
now, mechanical frequency combs have been hamstrung by the (generally) non-
linear phononic dispersion relation, which demanded high drive powers,
carefully designed mode frequencies or engineered mechanical non-linearities
to obtain an evenly spaced set of modes.

In a different regime of breakthrough physics, optical trapping has allowed us
to manipulate and control small particles ranging from single atoms [36, 37]
to micrometers [38] in size. These particles are confined by the potential
created by a strongly focused laser, which has enabled exploration of cutting-
edge fields in both fundamental physics and biology. Using optical traps
(tweezers), biologists can precisely manipulate anything from single strands
of DNA [39] to whole living cells [40]. Optically trapped ultra-cold atoms and
levitated nanoparticles are perfect testbeds for fundamental physics involving
gravity and mesoscopic quantum mechanics [41, 42, 43]. Although optical
trapping and frequency combs are widespread techniques, there is little direct
overlap between these regimes of physics.

In this work, we uniquely interface optical traps with frequency combs via a
mechanism that enables mechanical frequency combs without requiring feedback
control, external drives and frequency references, or precision optics. We
observe that due to a weak counterpropagating optical trap [44, 45], strongly
driven silicon nitride (Si3N4) membranes vibrate not only at their mechanical
eigenfrequencies, but at perfect integer multiples of a single frequency,
which can form a mechanical frequency comb. The standing-wave optical field
exerts a dielectrophoretic force on the membrane, which is modulated by the
membrane’s motion as it crosses extrema of the optical field. If the
displacement is small, this interaction can suppress mechanical dissipation
[46], but if the displacement becomes of the order of a quarter wavelength, it
creates integer multiple copies of the original membrane motion (”overtones”
[47], see the supplementary information (SI) Sec. A). These overtones share a
fixed frequency and phase relation, and thus form a frequency comb while
avoiding all difficulties of engineering a linear mechanical dispersion
relation. Strikingly, the overtones are solely dependent on the amplitude of
motion and the optical field, and thus independent of drive mechanism. We
utilize an optothermal parametric drive based on the same single, unmodulated
optical field to bring the membrane to self-oscillation. This allows us to
build a mechanical comb without any additional pump tone or frequency
reference, which makes overtone combs uniquely simple to generate. We will
first describe the mechanism that creates the overtones and study their
behavior in the frequency domain. Then, in the the time domain we will show
the fixed phase relation that makes the overtones act as a frequency comb.

##  II Results

###  II.1 Overtone and driving mechanism

Figure 1: Schematic of overtone frequency comb. a: The suspended Si3N4
membrane is subject to an out-of-plane laser, focused through a microscope
objective. Part of the light is reflected from the Si backplane and interferes
with the incident light, creating a counterpropagating wave optical trap such
that the field intensity (green) is periodic. The membrane is clamped to the
substrate, and motion is predominantly out-of-plane (purple, red dashed
lines). b: The optical field intensity causes a (spatially) periodic
modulation of the elastic potential through the dielectrophoretic force (red,
blue lines). For small motion of a membrane at frequency $\omega_{0}$ (cyan),
the restoring force component from the optical field switches sign twice per
oscillation (arrow pairs), same as the elastic potential. However, for larger
membrane motion (grey, orange), more optical extrema are crossed so the
optical field component switches sign multiple times per oscillation,
efficiently generating frequency components at $n\omega_{0}$ ($n=2,3,4,...$)
which form an overtone frequency comb.

In this section, we first describe the physical system and the mechanism that
creates the overtone frequency combs, and then introduce a more quantitative
model for their dynamics. The system consists of a suspended $t=$
$80\text{\,}\mathrm{nm}$ thick Si3N4 trampoline membrane [48] (see Methods),
shown schematically in Fig.1a. It rests $\sim 10$ $\mathrm{\SIUnitSymbolMicro
m}$ above a backplane formed by the silicon (Si) substrate, and a laser
($\lambda=633$ $\mathrm{nm}$, $p\lesssim 3$ $\mathrm{mW}$) from a commercial
Polytec MSA400 laser Doppler vibrometer is incident on the membrane. At this
wavelength, the Si3N4 reflects $\lesssim 30\%$ of the light, and Si reflects
about $35\%$. Part of the light thus forms a standing wave, Fig. 1a, periodic
in the direction of the mechanical motion.

The dielectric Si3N4 experiences a dielectrophoretic force proportional to the
gradient of the optical intensity, similar to a particle in a
counterpropagating-wave optical trap [44, 45]. The trap also exerts a
radiation pressure force, but this is negligible in our system ([6] and SI
Sec. C). If the dielectric moves (e.g. by driving a mechanical eigenmode), it
will experience a restoring force from its own elastic potential (blue solid
line in Fig. 1b), with an additional component from the optical field (red
solid line). For small motion ($|x|\ll\lambda/4$, cyan line in Fig. 1b), the
optical potential functions as an additional spring [46]. However, if the
motion of the membrane is of the order of the optical potential period
($\lambda/2$), the modulated potential generates the overtones.

This can be seen as follows: The restoring force of the mechanical potential
switches sign twice per oscillation (arrow pairs in Fig. 1b). If the motion is
large enough (grey, orange lines) for the resonator to cross multiple extrema
of the optical field, the optical component of the restoring force switches
direction $2n$ times per oscillation ($n$ integer number of optical extrema).
If the original motion was at mechanical eigenfrequency $\omega_{0}$, this
effect generates motional components at $n\omega_{0}$, which are the overtones
of the original eigenmode. We will show in Sec. II.3 that these overtones have
a fixed phase relation and thus form a mechanical frequency comb. Because the
overtones originate from the combination of the optical and mechanical
potential, they completely avoid the difficulties of engineering the
mechanical dispersion while still resulting in perfectly evenly spaced tones.
In contrast to other combs, the tones do not exist around some carrier
frequency.

Figure 2: Overtone frequency combs. a: Measured spectrum of a mechanical
frequency comb of the fundamental mode ($\omega_{0}=2\pi\times 118.049$
$\mathrm{kHz}$) of a suspended Si3N4 membrane, consisting of 35 peaks spaced
by $\omega_{0}$. The overtones are spectrally flat until $\sim 2200$
$\mathrm{kHz}$, after which their amplitude drops exponentially. Inset shows
the location of the laser spot to generate and read out the comb. The
additional peak around $2083\text{\,}\mathrm{kHz}$ are from a different
mechanical mode, see SI Sec. F. b: Simulated overtone comb spectrum, by
integrating Eq. (1), matching the measurement of a. c Simultaneous measurement
of fundamental-mode frequency comb (first four teeth) and higher-order
mechanical modes of the membrane. Insets: Laser spot position (left) and mode
shapes. d,e Mechanical frequency combs of the second mode
($257\text{\,}\mathrm{kHz}$) and third modes ($359\text{\,}\mathrm{kHz}$).
Difference in the noise floor is due to different decoders used due to
overloading. The $257\text{\,}\mathrm{kHz}$ mode has a near-degeneracy, so the
comb spectrum is less clean. Insets show the laser position to drive and read
out the different modes.

We model the overtone frequency comb using a resonator described by
displacement $x(t)$ ($x=0$ mechanical equilibrium), with resonance frequency
$\omega_{0}$ and decay rate $\gamma$, where we have divided by the simulated
effective mass $m_{\mathrm{eff}}\simeq 12\times 10^{-12}$ $\mathrm{kg}$. We
write the dielectrophoretic force as a proportionality constant
$F_{\mathrm{o}}$ (units of force) times the gradient of the periodic part of
the optical intensity, $\nabla
E^{2}\propto\sin\left(\frac{4\pi}{\lambda}(x-x_{\mathrm{off}})\right)$. This
way we can move finite-size effects of the membrane into $F_{\mathrm{o}}$,
which we numerically evaluate in SI Sec. C. We obtain the equation of motion

$\ddot{x}+\gamma\dot{x}+\omega_{0}^{2}x=\frac{F_{\mathrm{o}}}{m_{\mathrm{eff}}}\sin\left(\frac{4\pi}{\lambda}(x-x_{\mathrm{off}})\right),$
(1)

which requires only a suitable initial condition $|x|_{t=0}\gtrsim\lambda/4$
to demonstrate the creation of overtones. This condition is much larger than
typical interferometric position measurements, which is why we use a laser
Doppler vibrometer that is capable of resolving such large displacements (see
Methods). Eq (1) also shows that the overtones are independent of the choice
of drive (e.g. piezoelectric, electrostatic, thermal). By utilizing the
spatially-periodic optical field through optothermal effects [49], we can
create a parametric drive powerful enough to bring the membrane to self-
oscillation. That is, when the resonator moves through the field, the optical
intensity it experiences is modulated at twice the frequency of the original
motion. This modulates the resonator frequency through absorption (thermal
expansion), thus creating an optothermal parametric drive (see SI Sec. D). Our
membranes are patterned with a specific photonic crystal that enhances the
absorption of $633\text{\,}\mathrm{nm}$ laser light. This allows a single,
continuous-wave laser beam of $\mathrm{mW}$ power to bring the membrane into
self-oscillation. While we verify that the overtones can also be generated via
inertial (piezo) driving (see SI Sec. E), we leverage the optothermal drive to
avoid any external pump tone or frequency reference. This property of overtone
combs is unique within the mechanical frequency combs, and allows for
significantly simpler setups.

###  II.2 Overtone frequency comb

We measure the overtone frequency comb for the fundamental mode
$\omega_{0}=2\pi\times 118.049$ $\mathrm{kHz}$ of our membrane in Fig. 2a. In
the velocity power spectrum, we see a series of 35 peaks spanning the full
detection bandwidth of the setup ($4.2\text{\,}\mathrm{MHz}$), spaced by
$\omega_{0}$. This spacing is perfectly uniform, limited by the spectral
resolution of the setup to $4.7\times 10^{-8}$ relative spacing difference. A
comparison of this comb to other mechanical combs is included in the SI Sec.
A. Through integration of the velocity signal, we verify the displacement
$x\gg\lambda/4$. By numerically integrating Eq. (1), we can reproduce this
overtone comb, shown in Fig. 2b, using only the optical potential strength
$F_{\mathrm{o}}$, offset $x_{\mathrm{off}}$ and initial conditions
$[x_{0},v_{0}]$ as fit parameters (assuming steady state, so
$\gamma=F_{\mathrm{d}}=0$). This highlights that the non-linearity comes from
the optical field, without introducing mechanical nonlinearities previously
used to explain this behavior [47].

Figure 3: Comb dynamics. a: Measured displacement of membrane, showing clear
growth before reaching a plateau around $t=24$ $\mathrm{s}$. b: Measured
spectrum of membrane motion close to fundamental mode and two of its overtones
(5th and 15th). A thermal shift of $9\text{\,}\mathrm{Hz}$ of the fundamental
mode is visible as a $5\times 9$ and $15\times 9$ $\mathrm{Hz}$ shift in the
overtones (white dashed lines are horizontal to guide the eye). c: First six
overtones increasing in power as the membrane displacement increases (markers)
in a, with simulated dynamics (solid lines) based on integration of Eq. (1)
matching quantitatively to the highest four. d: Membrane displacement of a
different device (different chip) than a, showing growth from close to thermal
regime to steady state. e: Overtone amplitudes extracted from the spectrum of
measurement of d, showing all 35 overtones within the detection bandwidth
(black, dotted lines are overtones with numbers between the labeled ones). f:
Spectrum showing measured fundamental mode and some selected overtones in the
comb regime (blue) and thermal (green), along with Lorentzian fit with center
frequencies $n\omega_{0}$ ($n$ integer) and identical linewidths $\gamma_{0}$.

The overtones can be distinguished from the other mechanical eigenmodes of the
membrane. The mechanical eigenmodes can be detected and identified (Fig. 2c).
When the overtone comb is generated, additional peaks appear in the spectrum
(dotted black lines). This demonstrates that the overtones are not affected by
the mechanical dispersion relation and do not require engineered non-linear
resonances [12, 17]. Furthermore, it is possible to generate combs from the
second and third mechanical eigenmodes, by selecting the right laser position
on the membrane, Fig. 2d,e. This makes the frequency spacing variable, limited
by our ability to drive a particular eigenmode. Finally, for some laser
positions the overtone comb at $\omega_{0}$ interacts with a different
mechanical eigenmode ($\omega_{\mathrm{h}}$), which could allow extension of
the comb bandwidth (SI Sec. F).

We examine Eq. (1) to better understand the behavior of this mechanism
(details in SI Sec. A). Firstly, when the displacement $x$ increases, more
overtones appear due to the increasing number of optical extrema crossed. Each
overtone will again be modulated, so they have equal power up to a certain
cutoff, which is beneficial for many applications ($\sim 2200$ $\mathrm{kHz}$
in Fig. 2a,b). Secondly, the strength of the optical field ($F_{\mathrm{o}}$)
controls the power of each overtone relative to the original mode
$\omega_{0}$. Finally, the offset $x_{\mathrm{off}}$ between mechanical and
optical zeros determines the relative intensity of the odd and even number
overtones. The position offset $x_{\mathrm{off}}\simeq 40$ $\mathrm{nm}$ is
consistent between different membranes in this work, and is half the thickness
of the Si3N4 membranes. We can repeatably create the overtone comb in
different membranes, and study the effect of the optical beam itself on the
membranes and overtone comb in SI Sec. G.

###  II.3 Comb dynamics

The comb dynamics can be visualized by starting a measurement with the laser
spot positioned away from the membrane, as shown in Fig. 3a. We move the laser
to the membrane center at $t=$ $8.7\text{\,}\mathrm{s}$ such that the
optothermal parametric driving starts increasing displacement. We analyze the
dynamics by cutting the recorded time signal into intervals and performing a
Fourier transform on each. Then we concatenate the spectra such that we can
study its behavior over time (Fig. 3b) and monitor the power of individual
overtones by taking linecuts (Fig. 3c).

In Fig. 3b, the fundamental mode shows a $9\text{\,}\mathrm{Hz}$ upwards shift
in frequency, which is reproduced as an $n\times 9$ $\mathrm{Hz}$ upwards
shift for the $n$th overtone (e.g. $n=5,15$ in the figure) and likely
originates from slow thermalization (SI Sec. D). This illustrates thermal
tuning would be an effective mechanism to control and tune the comb spacing,
without compromising comb uniformity.

We plot the power of the overtones as the comb grows in Fig. 3c. We can
simulate and reproduce this growth quantitatively by integrating Eq. (1).
Comparison between the experimental data (markers) and simulation (solid
lines) in Fig. 3c shows good agreement with fit parameters
[$F_{\mathrm{o}}=3.8$ $\mathrm{pN}$, $x_{\mathrm{off}}=40$ $\mathrm{nm}$,
$x_{0}=5$ $\mathrm{nm}$, $v_{0}=1$ $\mathrm{nm}\text{\,}{\mathrm{s}}^{-1}$,
$\gamma/2\pi=0.8$ $\mathrm{Hz}$ and $F_{\mathrm{d}}=2.1$ $\mathrm{pN}$]. Most
important for the overtones is $F_{\mathrm{o}}$, which is in excellent
agreement with simulated force on the order of $\mathrm{pN}$ (SI Sec. C). All
extracted curves share a vertical offset to account for the total detection
efficiency. This shows Eq. (1) reproduces the dynamics of the resonator and
individual overtones.

In Figs. 3d,e, we observe the membrane motion from thermal state at $t<5$
$\mathrm{s}$ to the steady state of the comb at $t>33$ $\mathrm{s}$.
Initially, we detect only the fundamental mode $\omega_{0}$, until the
resonator starts crossing multiple optical extrema. Higher overtones appear as
the maximum displacement grows. In contrast to Fig. 2a, the fundamental mode
is the most powerful. The detection efficiency likely varies for each overtone
due to their shape [47]. The displacement stops growing at $t\simeq 26$
$\mathrm{s}$, and decreases slightly before the system reaches steady state
(see SI Sec. H). It is likely limited in amplitude by a mechanical non-
linearity, but the amplitude overshoot and the oscillation of amplitude of
overtones suggests that the interaction with the optical field play a role. In
the steady state, the fractional frequency stability of the 30th overtone is
$7.5\cdot 10^{-10}$ over a 6-hour period (SI Sec. H), limited by thermal
drifts. The frequency of the overtones determined solely by the mechanical
frequency, and is thus not affected by drifts in the laser frequency. Changes
in the laser power will affect the overtone amplitudes, mainly via the
optothermal parametric driving. The uniformity of the comb is not affected by
mechanical frequency shifts (nor by the laser), thus uniformity is constant.

We isolate and plot the spectrum around several of the overtones in Fig. 3f.
In the thermal regime (green markers), only the fundamental mode is visible
and we can fit a Lorentzian with linewidth $\gamma_{0}=2\pi\times 0.07$
$\mathrm{Hz}$ to the peak. For the comb in steady state, the fundamental mode
retains its linewidth $\gamma_{0}$, and we can derive that the comb mechanism
does not cause additional noise, SI Sec. I. All overtones possess the same
linewidth $\gamma_{0}$, which does not match the $n\gamma_{0}$ scaling
expected from a frequency-fluctuation limited linewidth, but matches a decay-
rate limited system. These properties combined should allow frequency combs
with single-$\mathrm{mHz}$ linewidths based on ultra-high-Q membrane
resonators [50].

To finally show that the overtones form a comb, i.e. that a fixed phase and
frequency relation exists, we show the time-domain signal in Fig. 4a. Unlike
combs centered around a carrier (such as soliton-based mechanical [17] and
optomechanical [8] frequency combs), there is no component with periodicity
longer than the mechanical period $2\pi/\omega_{0}$ (Figs. 4b,c). This
highlights the different physical and dynamical processes behind the overtone
comb. If the fundamental mode is the dominant component in the comb, we get a
sinusoid (Figs. 4b,c, blue line). When other components are dominant,
particularly $2\omega_{0}$ and $3\omega_{0}$, we get the peaked curves (Fig.
4c, green and orange lines), which retain the periodicity of the fundamental
mode. The shape of these curves proves that the entire overtone comb is phase-
coherent, as they are formed by a sum of cosines with the same phase offset
(see SI Sec. J). This behavior is retained not only in the steady-state, but
also during comb growth, Fig. 4d. There we plot the time signal at various
stages during the measurement of Fig. 3d,e, which shows a smooth transition
from sinusoidal to peaked behavior as the comb grows. Thus the overtones form
a mechanical frequency comb.

Figure 4: Time domain frequency comb. a: Measured velocity of resonator in
frequency comb regime, displaying no pattern with period longer than
$1/\omega_{0}$. b: Zoom-in of time-domain signal in a, showing the
$118\text{\,}\mathrm{kHz}$ fundamental mode dominates. c: Three different
measured time signals in steady state, with varying relative strengths of the
overtones (blue: $\omega_{0}$ dominates, orange: $2\omega_{0}$ and
$3\omega_{0}$ dominate, green: $1\omega_{0}$ to $6\omega_{0}$ similar in
strength). All traces can be reproduced only with phase-coherent addition of
the different overtones. d Time domain signal of comb during growth, same
measurement as Fig. 3e at the indicated times. This shows the smooth
transition from thermal regime (bottom) to the comb in steady state (top),
traces are offset vertically for clarity.

##  III Discussion

We have discovered a mechanism which uniquely interfaces two breakthrough
concepts: optical trapping and frequency combs. This allows for mechanical
frequency combs whose simplicity stands out from those based on previous
mechanisms. This is realized by integrating a suspended dielectric membrane in
a weak optical trap. The dielectrophoretic force from the optical field
modulates the mechanical potential of the membrane. This modulation creates
integer multiple copies (overtones) of the membrane’s motion, which forms a
frequency comb. We show combs of up to 35 overtones in a
$4.2\text{\,}\mathrm{MHz}$ bandwidth with control over the frequency spacing,
excellent uniformity, stability and no added mechanical noise. The integration
of the membrane in the optical trap allows us to combine the overtone comb
with optothermal parametric driving, which brings the membrane to self-
oscillate. We thus realize a frequency comb that requires no external drive or
control frequencies and which uses a minimal setup (laser, microscope
objective and vacuum chamber). This makes it more versatile and easier to use
than other ways of generating these combs. In summary, this mechanism unlocks
potential of mechanical frequency combs for sensing, timing and metrology
applications at the microscale, and provides native integration with phononic
circuits.

##  IV Methods

###  IV.1 Laser Doppler Vibrometer

Figure 5: Setup and membranes. a: Schematic of the laser Doppler vibrometer
(LDV) setup, where we access the chip containing the membranes using a
microscope objective directly outside a vacuum chamber. b: SEM image of the
suspended membrane, with nominal design parameters in white. c: Top-down view
of corner of the membrane, highlighting the interface of the photonic crystal
with the tether. d: Tilted view to show Si surface below c, to illustrate etch
roughness and an etched imprint of the photonic crystal pattern (bottom
right).

The setup used in this work is shown in Fig. 5a. It consists of a commercial
laser Doppler vibrometer (LDV), Polytec MSA400, which is depicted
schematically. Light from the LDV goes through a microscope objective, which
focuses it on the chip containing the membrane resonators. This chip is placed
in a vacuum chamber and pumped down until the pressure is $<5\times 10^{-6}$
$\mathrm{mbar}$, to reduce gas damping. There is a piezoelectric shaker
mounted to the sample holder, by which we can drive the membrane, though for
the majority of the measurements we use the thermal parametric driving
described in the SI Sec. D. The reflected light from the membranes is Doppler-
shifted due to their out-of-plane motion, which is then detected by the LDV.
To extend the time we can continuously measure, we add a Rohde & Schwarz
RTB2004 digital oscilloscope to readout the LDV decoder. Using the history
function of this oscilloscope, we can chain 16 measurements of 20 million data
points each, which allows $38\text{\,}\mathrm{s}$ of time signal at
$8.33\text{\,}\mathrm{MHz}$ sampling rate. However, this method for reading
out the velocity comes at a cost of the calibrated readout that the LDV
postprocessing offers.

The LDV outputs a voltage signal proportional to velocity or displacement,
depending on the LDV decoder used. If using the LDV postprocessing, it can
easily be Fourier-transformed and exported. If using the oscilloscope and
history function, we sequentially read out the measurements afterwards and
concatenate them in the correct order in post-processing. We can then further
extract information from this signal either by integrating (to obtain the
displacement), or by using Scipy’s short-time Fourier transform function to
obtain the time-varying behavior of the comb spectrum.

###  IV.2 Trampoline membranes

The membrane structures used in this work are fabricated out of
$100\text{\,}\mathrm{nm}$ thick stoichiometric low-pressure chemical vapor
deposition Si3N4 on top of a $1\text{\,}\mathrm{mm}$ Si chip. This was done by
writing the pattern using electron beam lithography and an inductively coupled
plasma (ICP) etch to transfer that pattern to the Si3N4. The membranes are
then released using a second ICP etch, at
$-120\text{\,}\mathrm{\SIUnitSymbolCelsius}$ for $30\text{\,}\mathrm{s}$,
resulting in about $10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ of undercut and
a final Si3N4 thickness of $80\text{\,}\mathrm{nm}$. A SEM image of a released
trampoline is shown in Fig. 5b, with the nominal design parameters added in
white. The resulting trampolines have well-characterized mechanical properties
[48, 51], with fundamental mode frequencies at $120\text{\,}\mathrm{kHz}$ and
Q-factors typically of 1 million.

The membranes are pattered with a periodic array of holes (Fig. 5c) that forms
a photonic crystal, though it is not designed to have a high reflectivity at
the operating wavelength of our LDV ($\lambda=633$ $\mathrm{nm}$). Instead,
the holes function as etch release holes to evenly release the membrane in the
final ICP etch, which is essential for the fabrication yield. For these
membranes, the periodic array of holes functions as a method to control the
mass and resonance frequencies [51], and possesses internal optical resonances
that facilitate absorption which we detail in SI Sec. D. Finally, the release
etch imprints the photonic crystal pattern on the Si backplane, which roughens
the surface. Fig. 5d shows the final Si surface roughness from the release
etch (top left half) and the imprint below the periodic hole pattern (bottom
right half).

## Data Availability

The data (raw data, analysis and calculation scripts, and finite element
simulations) that support the findings of this study are available at
https://doi.org/10.4121/19821016.v1.

## Acknowledgements

We would like to acknowledge Farbod Alijani, Peter Steeneken, Lior Michaeli,
and Albert Schliesser for interesting discussions, Dongil Shin for help with
the thermal simulations, and Wouter Westerveld and Gerard Verbiest for their
comments on the manuscript. R.N. would like to acknowledge support from the
Limitless Space Institute’s I2 Grant. M.J., A.C., S.G., and R.N. acknowledge
valuable support from the Kavli Nanolab Delft and from the Technical Support
Staff at PME, 3mE Delft, in particular from Gideon Emmaneel and Patrick van
Holst.

## Author contributions

M.J. performed the experiments, simulations and data analysis. A.G. and M.J.
developed the theory. A.C. developed the fabrication process and fabricated
the devices. M.J., A.G., S.G., and R.N. developed the interpretation of the
results. R.N. supervised the project. All authors contributed to writing and
editing the manuscript.

## Competing interests

The authors declare no competing interests.

##  V Supplementary information

This document contains the supplementary information.

###  V.1 Mechanical overtones

Figure S1: Mechanical overtones. a: Resonator modes of a string. b: Frequency
diagram of eigenmodes of a harmonic string and the anharmonic SiN trampoline
membranes used in this work. c: SiN trampoline modes, grey box denotes
degenerate modes. The color scale denotes normalized absolute displacement,
from minimum (blue) to maximum (red).

The motion of a mechanical resonator can typically be described by considering
a linear superposition of the normal modes of the structure. In some
structures, the normal modes are related by integer multiples of some
frequency, and they are said to be _harmonic_. This is the case for a model
string (Fig. S1a,b) and a desirable property for many musical instruments.
However for most of the resonators in this work, the normal modes are not
related by integer multiples, thus the spectrum is _anharmonic_. In Fig. S1c,
we have simulated and plotted the first six eigenmodes of the trampoline
membrane resonator of the main text, and plotted their frequencies in Fig.
S1b.

In this work, we introduce an optics-based non-linearity to the equation of
motion, which affects the description of the resonator motion in terms of the
normal modes given by its structure. Briefly summarized, we consider only the
lowest frequency normal mode given by the mechanical structure, as a simple
harmonic oscillator,

$\ddot{x}+\gamma\dot{x}+\omega_{0}^{2}x=0.$ (S1)

The solutions for this equation are well-known, $x\propto e^{i\omega_{0}t}$,
which oscillates at the frequency $\omega_{0}$ of the normal mode. By
introducing the gradient force term, $\propto\sin(x)$ (non-linear), we get
terms in our solution that oscillate at integer multiples of $\omega_{0}$,
without being related to the (other) normal modes of the mechanical resonator.
To distinguish them from the other normal modes, we will refer to the
components at $n\omega_{0}$ ($n$ integer) as overtones. Others have detected
these components of the mechanical motion before, and refer to them in a
similar manner [52, 47, 53]. They ascribe the origin of these overtones to be
mechanical in nature, while we propose an optical origin.

####  V.1.1 Small displacement

To expand the description of these overtones, we consider the equation of
motion as Eq. 1 in the main text, where we have added a resonant drive term
$F_{\mathrm{d}}$ and absorbed the effective mass $m_{\mathrm{eff}}$ into
$F_{\mathrm{o}}$ and $F_{\mathrm{d}}$ for convenience. The shorthand
$\Lambda=4\pi/\lambda$ is similarly used. We get

$\ddot{x}+\gamma\dot{x}+\omega_{0}^{2}x=F_{\mathrm{o}}\sin\left(\Lambda(x-x_{\mathrm{off}})\right)+F_{\mathrm{d}}e^{i\omega_{0}t}.$
(S2)

The periodic part of the optical potential is of the form
$E\propto\sin\left(\frac{2\pi}{\lambda}(x-x_{\mathrm{off}})\right)$. The
gradient of the optical intensity, $\nabla
E^{2}\propto\sin\left(\Lambda(x-x_{\mathrm{off}})\right)$ after some algebra
and absorbing the relevant constants into $F_{\mathrm{o}}$. The constant
offset between the potentials, $x_{\mathrm{off}}$, can be taken out of the
sine by

$\sin\left(\Lambda(x-x_{\mathrm{off}})\right)=\sin\left(\Lambda
x\right)s_{x}+\cos\left(\Lambda x\right)c_{x},$ (S3)

with $s_{x}=\cos(\Lambda x_{\mathrm{off}})$ and $c_{x}=\sin(\Lambda
x_{\mathrm{off}})$. For a small displacement $x$, we can use the Taylor series
and truncate the higher order terms, so $s_{x}\sin(x)+c_{x}\cos(x)\simeq
c_{x}+s_{x}x-\frac{c_{x}}{2}x^{2}$. The equation of motion then becomes

$\ddot{x}+\gamma\dot{x}+\omega_{0}^{2}x=F_{\mathrm{o}}\left(c_{x}+s_{x}\Lambda
x-\frac{c_{x}}{2}\Lambda^{2}x^{2}\right)+F_{\mathrm{d}}e^{i\omega_{0}t}.$ (S4)

This equation admits solutions of the form

$x=\sum_{n=0}^{\infty}A_{n}e^{in\omega_{0}t},$ (S5)

which are the integer multiples of our original frequency $\omega_{0}$; the
overtones. The term $n=0$ corresponds to a static position offset from the
zero of the mechanical potential. Substituting the solution Eq. (S5) into Eq.
(S4) and gathering all terms by their frequency ($e^{in\omega_{0}t}$ for every
$n$ separately) allows us to extract the amplitudes of the individual
overtones in terms of the parameters of our system. We get

$\displaystyle A_{0}$ $\displaystyle\simeq c_{x}F_{\mathrm{o}}/\omega_{0}^{2}$
(S6) $\displaystyle A_{1}$
$\displaystyle=\frac{F_{\mathrm{d}}}{i\gamma\omega_{0}-\Lambda
F_{\mathrm{o}}s_{x}+F_{\mathrm{o}}\frac{c_{x}}{2}\Lambda^{2}A_{0}},$
$\displaystyle A_{2}$
$\displaystyle=\frac{F_{\mathrm{o}}\frac{c_{x}}{2}\Lambda^{2}A_{1}^{2}}{3\omega_{0}^{2}-2i\gamma\omega_{0}+F_{\mathrm{o}}s_{x}\Lambda-
F_{\mathrm{o}}\frac{c_{x}}{2}\Lambda^{2}A_{0}},$ $\displaystyle A_{n}$
$\displaystyle=\frac{F_{\mathrm{o}}\Lambda\sum_{j,k=0}^{\infty}A_{j}A_{k}}{(n^{2}-1)\omega_{0}^{2}-i\gamma
n\omega_{0}+F_{\mathrm{o}}s_{x}\Lambda-
F_{\mathrm{o}}\frac{c_{x}}{2}\Lambda^{2}A_{0}},$

where the summation contains only the terms where $j+k=n$. This sequence of
overtone amplitudes is monotonically decreasing (see Fig. S2): Every
$A_{n}<A_{n-1}$. In the small-displacement case, the overtones will be
negligible compared to $A_{1}$.

####  V.1.2 Large displacement

If instead the displacement is not small, we can still use the same method. We
need to keep all the terms of the Taylor expansion,

$\sin(x)=\sum_{n=0}^{\infty}\frac{(-1)^{n}x^{2n+1}}{(2n+1)!},\quad\cos(x)=\sum_{n=0}^{\infty}\frac{(-1)^{n}x^{2n}}{(2n)!}.$
(S7)

We can use the same ansatz of Eq. (S5), and extract the amplitudes by
collecting all terms of the same frequency. For $n=0$, the solution converges
to

$A_{0}\simeq c_{x}F_{\mathrm{o}}/\omega_{0}^{2},$ (S8)

as the contribution from the higher order terms of the expansion scale with
$1/\Lambda^{n}$. Similarly, the denominator for $n>0$ can be truncated to
obtain

$\displaystyle A_{1}$
$\displaystyle=\frac{F_{\mathrm{d}}}{i\gamma\omega_{0}-\Lambda
F_{\mathrm{o}}s_{x}+F_{\mathrm{o}}\frac{c_{x}}{2}\Lambda^{2}A_{0}}$ (S9)
$\displaystyle A_{n}$
$\displaystyle=\frac{F_{\mathrm{o}}\sum_{j}\left(s_{x}s_{j}\sum\limits_{k,\ell,...=0}^{\infty}A_{k,\ell,...}^{j}+c_{x}c_{j}\sum\limits_{k,\ell,...=0}^{\infty}A_{k,\ell,...}^{j}\right)}{(n^{2}-1)\omega_{0}^{2}-i\gamma
n\omega_{0}+\Lambda
F_{\mathrm{o}}s_{x}-F_{\mathrm{o}}\frac{c_{x}}{2}\Lambda^{2}A_{0}},$

with $s_{j}=\frac{(-1)^{(j-1)/2}\Lambda^{j}}{j!}$ for odd $j$ ($s_{j}=0$ for
even $j$) and $c_{j}=\frac{(-1)^{j/2}\Lambda^{j}}{j!}$ for even $j$ ($c_{j}=0$
for odd $j$). The summation only contains the terms of $j$ amplitudes
($A_{k},A_{\ell},...$) whose indices $k,\ell,...$ add up to $n$. To clarify,
for $n=3$, we sum over $A_{1}A_{1}A_{1}$, $A_{1}A_{2}$, $A_{2}A_{1}$,
$A_{3}A_{0}$, $A_{0}A_{3}$ as well as many terms with more $A_{0}$’s. For
every $n$, we get an infinite number of contributions to each term. The
denominator contains at most two terms that depend on $n$, but the rest is
constant and the whole denominator converges for each $n$; thus we can
truncate the denominator, as we have done in Eq. (S9). In contrast, to
evaluate the numerator we can take sequentially higher terms in the Taylor
expansion (i.e. increase $j$), find all combinations of $j$ integers that sum
to $n$, and add them to the term of frequency $n\omega_{0}$.

It is difficult to analytically express all the combinations of $j$ integers
that sum up to $n$, so we numerically evaluate the overtone amplitudes in Fig.
S2. Here, we take the expansion up to $13^{\mathrm{th}}$ order, which results
in close to $2.5\cdot 10^{6}$ terms contributing to $n=12$. We fix the terms
with $n=0,1$ to a constant value, in the experiment we do not directly apply a
resonant drive ($A_{1}$) and there are radiation pressure effects that could
result in a static position offset ($A_{0}$). From the numerical evaluations,
we consistently see that for small displacement, the power in each subsequent
overtone drops exponentially. However, for large displacement this is no
longer true, and the higher-order terms from the expansion cause strong
overtones that decrease much slower in amplitude with overtone number. This is
reproduced in the experiments shown in the main text.

Figure S2: Overtone amplitudes. Numerically evaluated overtone amplitudes
$A_{n}$ in the small- and large-displacement case. The terms $A_{0}=1\cdot
10^{-7}$ and $A_{1}=3\cdot 10^{-7}$ are taken constant, with $F_{\mathrm{o}}$
corresponding to $100\text{\,}\mathrm{pN}$. For large displacement, the
overtones have much higher amplitude than in the small-displacement case. The
amplitude also decays more slowly with mode number.

####  V.1.3 Numerical simulation of overtones

So far, we have derived analytically that the addition of an optical $\sin(x)$
non-linearity to a mechanical harmonic oscillator creates components of
mechanical motion at integer multiples of the original frequency $\omega_{0}$.
If the displacement is large enough, these overtones have significant
amplitude and form a frequency comb. We have made some simplifying assumptions
(e.g. $x_{\mathrm{off}}=0$), which do not necessarily hold in practice. The
remedy that, we perform numerical simulations. These allow us to better
understand the roles of resonator displacement $x$, the the optical
(dielectrophoretic) force $F_{\mathrm{o}}$, and the position offset
$x_{\mathrm{off}}$ of the optical potential with respect to the mechanical
potential.

Figure S3: a: Numerical simulation of overtones generated depending on initial
amplitude. Parameters $F_{\mathrm{o}}=100$ $\mathrm{pN}$,
$\omega_{0}=2\pi\times 118.049$ $\mathrm{kHz}$, $\gamma=2\pi\times 0.02$
$\mathrm{Hz}$, $F_{\mathrm{d}}=0$ $\mathrm{N}$, $x_{\mathrm{off}}=30$
$\mathrm{nm}$. b: Simulation of overtones generated depending on optical force
$F_{\mathrm{o}}$, for $x_{o}=200$ $\mathrm{nm}$. Other parameters identical to
a. Inset shows the frequency shift for stronger optical fields. c: Simulation
of overtones depending on offset $x_{\mathrm{off}}$. For $x_{\mathrm{off}}=0$
$\mathrm{nm}$, we predominantly drive the even overtones through the cosine-
expansion, while for $x_{\mathrm{off}}=\lambda/8=79$ $\mathrm{nm}$ we
predominantly drive the odd overtones through the sine-expansion. At
$x_{\mathrm{off}}=35$ $\mathrm{nm}$, odd and even modes are approximately
equal in power. Individual traces are offset vertically, $x_{0}=1000$
$\mathrm{nm}$ and other parameters are identical to a.

First, we investigate the dependence of the frequency comb on displacement
$x$. In Fig. S3a, we simulate the motion of the resonator starting from the
value initial amplitude $x_{0}$ ($v_{0}=0$) indicated in the legend. We
include dissipation but it is sufficiently small that the amplitude does not
significantly decrease within the simulation time. For small amplitude of
motion ($x_{0}<x_{\mathrm{off}}$), we cross no optical extrema and only the
fundamental mode is visible, no higher overtones appear. Once the motion is
large enough to cross one optical extremum ($100\text{\,}\mathrm{nm}$),
several overtones appear with exponentially decreasing power, as derived
previously. If the motion is large enough to cross multiple optical extrema
($x_{0}>\lambda/4\simeq 150$ $\mathrm{nm}$), more overtones appear. The lower
overtones have approximately equal amplitude (empirically until
$\omega=3x_{0}/(\lambda/4)\times\omega_{0}$), since these also interact with
the optical potential and get modulated to drive higher overtones. Note that
the lower overtones have the same amplitude regardless of the resonator
displacement. The higher overtones, above
$\omega=3x_{0}/(\lambda/4)\times\omega_{0}$, have an amplitude that decays
exponentially.

Secondly, we study the dependence of the frequency comb on dielectrophoretic
force $F_{\mathrm{o}}$. For identical initial position $x_{0}$, we plot the
frequency comb for various values of $F_{\mathrm{o}}$ in Fig. S3b. For
$F_{\mathrm{o}}=0$ $\mathrm{N}$, we see the dominant fundamental mode but also
some higher harmonics (3,5,7,…). These are from the numerical accuracy of our
simulation, they do not follow the exponentially decaying trend of the other
curves. For increasing values of $F_{\mathrm{o}}$, the amplitude of each of
the individual overtones increases linearly, but the exponential fall-off for
higher overtone numbers is unchanged. At some point, the non-linearities start
to shift the frequency of the modes, but this trap strength is well beyond the
regime of our setup ($F_{\mathrm{o}}=100$ $\mathrm{nN}$ corresponds to
approximately $6\text{\,}\mathrm{W}$ of incident laser power).

Finally, we report on the dependence of the frequency comb on position offset
$x_{\mathrm{off}}$. In the analytical case, we simplified using
$x_{\mathrm{off}}=0$ $\mathrm{nm}$, retaining only the cosine terms so the
even overtones are strong (dark blue). By shifting the optical potential with
respect to the mechanical zero by $\lambda/8=79$ $\mathrm{nm}$, we see mainly
the odd-numbered overtones appear. For this value of $x_{\mathrm{off}}$, we
can rewrite the cosine into a sine, such that the Taylor-series expansion only
contains odd terms. At other values of $x_{\mathrm{off}}$, we have a weighted
average of the two series expansions. In the middle, at $x_{\mathrm{off}}=35$
$\mathrm{nm}$, we see odd and even numbered modes with approximately equal
power.

To summarize: The addition of an optical $\sin(x)$ non-linearity to a
mechanical harmonic oscillator creates components of mechanical motion at
integer multiples of the original frequency $\omega_{0}$, which we call the
overtones of $\omega_{0}$. These overtones can have significant amplitude, if
the mechanical displacement is large enough. The displacement $x$ controls the
number of overtones visible, while the dielectrophoretic force
$F_{\mathrm{o}}$ determines their power relative to the fundamental mode of
the comb. Based on the offset between the optical and mechanical potential
$x_{\mathrm{off}}$, we can control the relative power in the even- or odd-
numbered overtones. We have thus derived an analytical model for the overtone
frequency comb, and expanded our qualitative understanding with numerical
simulations.

####  V.1.4 Exclusion of mechanical nonlinearity

We have introduced an optical nonlinearity to explain the mechanical frequency
combs, but we have not excluded conventional mechanical nonlinearities (e.g. a
Duffing term, $\propto x^{3}$). These nonlinearities can also lead to
mechanical frequency combs, even in the absence of other modes to couple to
[47, 23]. In Fig. S4a, we have simulated and plotted the spectrum of a
mechanical resonator with a Duffing term $c_{\mathrm{duff}}x^{3}$, such that
the equation of motion is

$\ddot{x}+\gamma\dot{x}+\omega_{0}^{2}x+c_{\mathrm{duff}}x^{3}=0.$ (S10)

We start from an initial condition of $x_{0}=500$ $\mathrm{nm}$ for all
traces. For a sufficiently strong nonlinearity, we see a frequency comb.
However, due to the nonlinearity, the frequency is shifted away from
$\omega_{0}$. This is expected behavior for nonlinear resonators, oscillating
at $\omega_{0}$ for small amplitude frequency shifting for larger amplitudes
where the nonlinearity becomes dominant. The sign of $c_{\mathrm{duff}}$
controls the direction of the frequency shift, and a hardening nonlinearity
(increasing frequency) is expected for most mechanical resonators.

In Fig. S4b, we have plotted the fundamental mode during the measurement shown
also in Fig. 3e. At the start of the measurement, the amplitude is small and
there is no frequency comb, whereas at the end the amplitude is large and the
comb is clear. There is no shift of the frequency of the fundamental mode
visible in this measurement, which excludes the explanation of a Duffing
nonlinear term being the origin of the comb.

It was proposed by others that strong, higher order nonlinearities (e.g.
$\propto x^{3}$, $x^{5}$, $x^{7}$, etc.) could be the origin of this comb
[47]. However, to reproduce the frequency comb without a similar frequency
shift as in the Duffing case, one would need strong nonlinear terms of very
high orders. Rather than assuming as many nonlinear terms as we have overtones
(up to $35$), we have introduced an optical nonlinearity that can explain all
overtones from a single effect.

Figure S4: Frequency comb by mechanical nonlinearity. a: Simulated
displacement spectrum of a mechanical resonator with a Duffing nonlinearity.
All traces share the same initial displacement, $x_{0}=500$ $\mathrm{nm}$,
with different Duffing coefficients. For a sufficiently large nonlinearity,
the frequency comb appears, but the frequency of the fundamental mode (and
comb spacing) is changed. b: Measured frequency of the fundamental mode of
Fig. 3e, showing no shift in frequency between low displacement (no comb) and
high displacement (comb).

####  V.1.5 Comparison to other combs

There have been demonstrations of mechanical frequency combs generated via
different mechanisms in literature. For applications, the comb properties are
more important than the generation mechanism, and we have summarized a large
part of the combs available from literature in Table S1. The entries are
sorted by material/geometry platform, and multiple publications with subtle
difference in the comb generation mechanism may be present within the
individual categories.

The overtone comb of this work has a bandwidth on par with the largest
available in literature (limited by the setup), but also has a relatively
large resolution. The bandwidth is given by a combination of the maximum
displacement and mechanical frequency, both can easily be engineered. We
estimate the optothermal parametric drive is limited to $\lesssim 500$
$\mathrm{kHz}$ in our membrane, which would constrain the mechanical
frequencies that can be operated without an external drive. The frequency
resolution could be enhanced by lowering the mechanical frequency (i.e. make
the membrane tethers longer and thinner), and can likely be tuned thermally. A
further benefit is that the overtone combs are spectrally flat (same amplitude
in different overtones), in contrast to other mechanical combs.

The main benefit of the overtone comb is that it does not require an external
drive, which is unique for a mechanical comb. This greatly simplifies the
necessary setup to operate the comb. In terms of uniformity of the frequency
spacing, overtone combs come close to the level of electromechanical combs,
where the uniformity is given by the electric drives directly. However, our
measured uniformity is limited by the detection setup. Overtone combs are also
among the most stable, especially comparing our 6-hour frequency stability
(Sec. V.8) to the $10\text{\,}\mathrm{s}$ Allan deviation of Ref. [54].

Type | Bandwidth ($\mathrm{Hz}$) | Resolution ($\mathrm{Hz}$) | Ext. drive power (dBm) | Uniformity (-) | Stability (-) | Footprint (${\mathrm{\SIUnitSymbolMicro m}}^{2}$)
---|---|---|---|---|---|---
Electromechanical [11] | $20-100$ |  | $0.005-10$ |  |  $-15$ [15] | $1\cdot 10^{-9}$ |  | $3\cdot 1-150\cdot 50$
Nanostrings [16] | $10-25$ | $(\times 10^{3})$ | $0.5-10$ | $(\times 10^{3})$ |  $-2.5$ [23] |  |  | $55\cdot 0.27$
Free-beam [17] | $60-150$ | $(\times 10^{3})$ | $2-10$ |  $(\times 10^{3})$ [26] | $5$ |  $5\cdot 10^{-6}$ [54] |  $1\cdot 10^{-8}$ [54] | $1100\cdot 350$
Coupled beams [18] | $10-200$ |  |  $1-30$ [55] |  | $-26$ |  |  | $500\cdot 50$
Bulk acoustic [20] | $20$ |  | $0.7-2$ |  | $-67$ |  |  | $23000\cdot 1000$
Circular 2D [22, 27] | $2-11$ | $(\times 10^{6})$ | $75-400$ | $(\times 10^{3})$ | $10$ |  |  |  $5-8$ (diameter)
This work | $4.2$ | $(\times 10^{6})$ | $118-359$ | $(\times 10^{3})$ | None | $4.7\cdot 10^{-8}$ | $7.5\cdot 10^{-10}$ | $750\cdot 750$
Table S1: Mechanical frequency comb performance. Overview of reported
mechanical frequency comb properties from literature, organized by platform
(irrespective of comb mechanism). The overtone comb has a bandwidth (frequency
span) on par with the largest combs, but also a relatively large resolution
(frequency spacing). We report the absolute frequency stability over a 6-hour
period (see SI Sec. V.8), whereas [54] computes the Allan deviation over a
$10\text{\,}\mathrm{s}$ period. Our work is the only mechanical frequency comb
that does not require any external drive.

###  V.2 Higher harmonics in LDVs

There are several mechanisms by which spurious higher harmonics can appear in
measurements from a LDV [56]. In the following paragraphs, we will discuss
these mechanisms and show that we can exclude them as the source of the
observed frequency comb. We will treat all parts of the setup shown in Fig. 5:
the optical parts, the photodetector, the decoder and the data acquisition.

Multi-path interferences may happen when the LDV laser beam is reflected by
more than one surface [57]. The geometry of our system facilitates having
multiple reflection sources, since both the trampoline membrane and the
substrate contribute to the reflection. We assume an ideal LDV, with an
incident electrical field with amplitude $E_{0}$, a moving membrane, and a
stationary substrate, as shown in Fig. S5a. The directly reflected (desired)
beam from the membrane has power $R_{\mathrm{SiN}}|E_{0}|^{2}$
($R_{\mathrm{SiN}}=0.3$). Assuming no scattering ($T+R=1$, for transmission
$T$ and reflection $R$), the substrate contributes
$T_{\mathrm{SiN}}^{2}R_{\mathrm{Si}}|E_{0}|^{2}\simeq 0.17|E_{0}|^{2}$
($R_{\mathrm{Si}}=0.35$). Only the membrane contribution has a Doppler shift
due to its motion, $+v$. Neglecting effects from beam divergence, multiply-
reflected beams contribute
$T_{\mathrm{SiN}}^{2}R_{\mathrm{Si}}^{2}R_{\mathrm{SiN}}$,
$T_{\mathrm{SiN}}^{2}R_{\mathrm{Si}}^{3}R_{\mathrm{SiN}}^{2}$, and so on.
These terms have Doppler shifts with the opposite sign to the directly
reflected beam, as they interact with the membrane from the opposite
direction. The first term has a Doppler shift equal in magnitude to the
desired signal ($-v$), while subsequent terms have multiple times the Doppler
shift ($-2v,-3v,...$). These terms thus directly affect the magnitude of the
velocity that the LDV senses. For harmonic membrane motion
$v(t)=v_{0}\cos(\omega_{0}t)$, the multiply-reflected beams lead to an error
in $v_{0}$, but their Doppler shift does not contribute directly to higher
harmonics of the observed membrane motion.

Figure S5: Sources of higher harmonics. a: Schematic of reflected components
of a single optical beam on our device (red arrows). For each component, the
Doppler shift due to membrane velocity $v$ is shown, as well as the relative
power. b: Frequency comb measured using the VD-09 decoder (maximum bandwidth
set to $1.5\text{\,}\mathrm{MHz}$ to allow maximum detectable amplitude).
Inset shows microscope image with laser position highlighted by the red
circle. c: Frequency comb measured using the DD-300 decoder, measured on the
same device at a slightly different position.

In a non-ideal laser Doppler vibrometer, multi-path interferences can be
identified by presence of ripples and spikes in the velocity signal. We can
compare the measured velocity signals with simulated signals using the model
of [57]. For a harmonic oscillator with $\omega_{0}=118$ $\mathrm{kHz}$,
maximum displacement amplitude $x_{\mathrm{max}}$, the demodulated velocity
signal has the form

$v(t)=-\frac{1}{2}\omega_{0}x_{\mathrm{max}}\left(\frac{\theta_{2}-1}{2\theta\cos\left(\frac{4\pi}{\lambda}\cos(\omega
t)+\Delta\phi\right)+\theta^{2}+1}\right)\sin(\omega t).$ (S11)

Here, $\theta$ is the electric field ratio of the correct and the unwanted
reflected signals ($\theta=\infty$ for an ideal vibrometer), $\Delta\phi$ is
the phase offset between the two beams (constant). In the worst-case scenario
($R_{\mathrm{SiN}}=0.3$, $R_{\mathrm{Si}}=0.35$, no scattering losses),
$\theta\simeq 1.33$. The unwanted beam likely suffers from scattering losses
more than the correct beam, since its path crosses the Si3N4 domain twice and
reflects of the etch-roughened Si surface (see Fig. 5d).

We simulate the velocity signal with multi-path interferences, and compare it
with the observed velocity signal in Fig. S6. The simulated signal of Fig. S6b
($\theta=10$) has the characteristic ripples of multi-path interference. By
comparing the simulated signal to the measured signal (Fig. S6a), one sees
that the characteristic ripples are absent. While it is unlikely that the
multi-path interference effect is completely absent in our setup, it can be
corrected by a combination of amplitude- and phase-locked loops [58]. Either
way, the observed velocity signal does not correspond to multi-path
interference effects.

We can further exclude multi-path interferences as the source of the comb by
measuring outside the membrane, on the Si chip. At the position indicated on
the microscope image of Fig. S6d, the Si and Si3N4 layers are touching and
thus have the same motion. The multi-path interference as described in Ref.
[57] does not happen. We use a second beam (’drive beam’ in Fig. S6d) to
generate a comb of the fundamental mode of the membrane. At the measurement
position outside the membrane, we can detect this comb, red curve in Fig. S6c.
For verification, without the drive beam present (blue curve), no comb is
detected while the fundamental mode is visible. Thus we can exclude multi-path
interference effects as the source of the observed frequency combs.

Figure S6: Exluding multi-path interferences. a: Measured velocity signal for
different times in Fig. 4b. The data has been smoothed with a first order
Savitzky-Golay filter of width $3$, to emphasize the presence of ripples. b:
Simulated velocity signal of a resonator with frequency
$118\text{\,}\mathrm{kHz}$ with multi-path interferences present. The
characteristic ripples of the interference are visible at the extrema of the
velocity curve in b, but absent in a. c: Measured velocity spectrum outside
the membrane, on chip. A second (drive) beam generates a frequency comb of the
membrane’s fundamental mode, which we detect at the position indicated on the
microscope image d. Without this beam on the membrane, no comb is observed.

High-aperture effects can occur when the LDV is operated with a microscope
objective, based on the Guoy phase delay of the optical beam [56, 59]. The
majority of the measurements in this work were done using a Mitutoyo Plan APO
5x objective with a numerical aperture of $0.14$. With this numerical
aperture, the amplitude error and harmonic distortion should be limited to
$<1\%$ [56]. Additionally, we compare the frequency combs obtained from the
same device using different lenses (different numerical apertures) in Sec.
V.7, which show no relation between the numerical aperture and frequency comb.
Thus we can exclude high-aperture effects as the source of the observed
frequency comb.

The electronic parts of the LDV (sketched in Fig. 5a) may contribute to a
nonlinear response that would generate harmonics [56]. We record frequencies
well below the $20\text{\,}\mathrm{MHz}$ specified maximum frequency of the
detector, thus this component should not contribute higher harmonics.

In a conventional heterodyne interferometer, harmonics of a signal would
naturally arise for sufficiently large displacement. The interferometer phase
is linear only in a limited regime, typically $x\ll\lambda$. However, a
Doppler vibrometer measures the frequency shift of the reflected light rather
than its phase. The frequency shift can be extracted with a frequency-to-
voltage converter, in the LDV decoder. This allows an LDV to measure
velocities associated with displacements much larger than the optical
wavelength, without suffering from the harmonics expected in a displacement
interferometer.

There are two LDV decoders available in our setup, the VD-09 velocity decoder
and the DD-300 displacement decoder. The VD-09 has a lower maximum operating
frequency ($1.5-2.5$ $\mathrm{MHz}$), while the DD-300 has a lower detection
maximum. In Fig. S5b,c, we show the frequency combs measured from the same
device with the different decoders. These measurements show that the frequency
comb is not a decoder artifact. Thus we can exclude the electronic parts of
the measurement setup as a source of the observed frequency comb.

In summary, we have described common mechanisms that yield higher harmonics in
laser Doppler vibrometer measurements. All these mechanisms can be excluded as
the cause of the frequency combs, by qualitative arguments and the
measurements shown in Figs. S5, S6, and S17. In combination with the fact that
there is a visible change in the membrane when the frequency comb occurs (Fig.
S15), we conclude that the observed mechanical frequency comb is a real,
physical phenomenon and not an artifact of the measurement setup.

###  V.3 Membrane in an optical trap

Figure S7: Simulation of optical trap. a: Electric field norm of incident
$633\text{\,}\mathrm{nm}$ light on $100\text{\,}\mathrm{nm}$ Si3N4 membrane
suspended $\sim$$6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ above a Si
backplane. The reflected light forms a standing wave pattern. Zoom shows
center of the Si3N4 membrane. Bottom inset shows the simulation geometry
schematically. b: Force exerted by the optical field on the Si3N4 domain as
the position of the membrane is varied. The colors represent different
membrane thicknesses. c: Normalized strength of the trap for different Si3N4
thicknesses.

Optical traps can be used to confine small particles, as the light exerts a
force proportional to the gradient of the optical intensity [44, 60]. For
uncharged dielectric particles, the electrical field polarizes the particles,
and any change in field (from a gradient) will result in a dielectrophoretic
force acting on that particle [61]. In this section, we will simulate a laser
incident on a dielectric (Si3N4) membrane, and evaluate the dielectrophoretic
and radiation pressure forces. We will show that the optical gradient is
sufficient for the dielectrophoretic force to affect the mechanical motion,
while the radiation pressure force does not affect the dynamics significantly.

####  V.3.1 Dielectrophoretic force

We simulate the physical system consisting of a $t\simeq$
$80\text{\,}\mathrm{nm}$ thin Si3N4 slab suspended in air above a much thicker
Si slab, shown in Fig. S7a. We study a 2D-axisymmetric domain centered around
our laser spot, which is incident from the top of the figure. We assume that
our laser Doppler vibrometer emits a Gaussian beam for convenience. For a
regular optical trap, the optical gradient and thus the force exerted is
weakest along the laser propagation direction, while it is stronger in the
plane normal to the propagation direction. However, the Si backplane in our
geometry reflects part of the light and produces a standing wave (Fig. S7a),
producing a periodic optical intensity. Combined with the fact that the
membranes are much larger than the optical gradient laterally ($\sim 750\times
750$ ${\mathrm{\SIUnitSymbolMicro m}}^{2}$), we have a one-dimensional
counterpropagating-wave optical trap in the out-of-plane direction. A simple
(unpatterned) Si backplane has been shown before to be sufficient to create an
interference pattern for such an optical trap [45].

We consider only motion in the out-of-plane direction of the membrane, so we
are interested in the gradient of the electric field norm in the $y$-direction
(Fig. S7a). The total force $F_{\mathrm{de}}$ induced by the optics on the
dielectric Si3N4 structure can then be found by integrating this gradient over
the membrane domain,

$F_{\mathrm{de}}=\int_{V_{SiN}}\alpha\frac{\partial|E|^{2}}{\partial
y}\mathrm{d}V,$ (S12)

with $\alpha$ the polarizability. For a linear isotropic material, the
constant $\alpha=P_{x}/E_{x}$ can be found from the simulation by dividing the
polarization in any direction $P_{x}$ by the electric field in that direction
$E_{x}$.

We need relate the equation for $F_{\mathrm{de}}$ (units of force) to the
coefficient $F_{\mathrm{o}}$ (units of force) used in Eq. 1. The expression
for $F_{\mathrm{de}}$ depends on the position of the resonator with respect to
the Si backplane, Fig. S7b, whereas $F_{\mathrm{o}}$ is a constant.
$F_{\mathrm{de}}(x)$ in Fig. S7b shows the periodic behavior expected from an
optical trap (i.e. periodic sign change), as well as a slight asymmetry due to
radiation pressure. It shares the periodicity that we have separately
introduced in Eq. 1 in the main text,
$\sin\left(\frac{4\pi}{\lambda}(x-x_{\mathrm{off}})\right)$, but does not
exactly follow the same curve due to the finite size (thickness) of the
membrane with respect to the wavelength. To simplify the dynamical
calculations, we have thus chosen to separate out the periodic behavior and
leave it proportional to the optical intensity,
$\sin\left(\frac{4\pi}{\lambda}(x-x_{\mathrm{off}})\right)$, and keep
$F_{\mathrm{o}}$ as a constant equal to the maximum of $|F_{\mathrm{de}}(x)|$.
Thus for a $t=$ $80\text{\,}\mathrm{nm}$ membrane subject to our
$3.6\text{\,}\mathrm{mW}$ laser, we expect a maximum out-of-plane force in the
order of $60\text{\,}\mathrm{pN}$. This confirms the values of
$F_{\mathrm{o}}$ used in the numerical simulations in the main text and in
Sec. V.1, which show that the dynamics of our membrane can be significantly
altered by the presence of the optical field.

The different curves of Fig. S7b illustrate that the total optical force
exerted on the membrane does not decrease monotonically with membrane
thickness. When we calculate the effective strength of the trap (peak-to-peak
force difference divided by the thickness of the membrane, in Fig. S7c), we
see that it decreases for increasing Si3N4 thickness. There are minima around
$t=\lambda/(2n_{\mathrm{SiN}})$ where $n_{\mathrm{SiN}}$ the index of
refraction; the half-integer number of optical periods leads to a total
cancellation of the net force on the membrane $F_{\mathrm{de}}$.

####  V.3.2 Radiation pressure force

In an optical cavity, the interaction between radiation pressure and a
mechanically compliant element can lead to optomechanical frequency combs [6].
The radiation pressure of an optical cavity coupled to a mechanical resonator
can result in an effective cubic non-linearity, and by sufficiently driving
the optical cavity one can bring the system in the self-oscillation regime.
This creates an optical frequency comb with a spacing equal to the mechanical
frequency. We must thus distinguish the frequency combs from the
dielectrophoretic force from those originating from the radiation pressure.

We follow the analysis of Miri et al. [6] as they analyze the stability of a
mechanically compliant end-mirror optomechanically coupled to an optical
cavity. Our Si3N4 membrane functions as the compliant end mirror, while the Si
substrate is the other (fixed) end mirror. For a given set of parameters, we
evaluate whether the frequency comb will appear through equations (9a,b) of
[6], which amounts to checking the Routh-Hurwitz stability criterion. For
completion, the conditions for stability are

$\displaystyle 3\mathcal{G}^{2}+4\Delta\mathcal{G}+\Delta^{2}+\kappa^{2}/4$
$\displaystyle>0,$ (S13)
$\displaystyle\mathcal{G}^{4}+d_{1}\mathcal{G}^{3}+d_{2}\mathcal{G}^{2}+d_{3}\mathcal{G}+d_{4}$
$\displaystyle>0,$

with coefficients

$\displaystyle d_{1}=$ $\displaystyle 4\Delta,$ (S14) $\displaystyle d_{2}=$
$\displaystyle\gamma\kappa+6\Delta^{2}+\gamma^{2}-6\omega_{0}^{2}+\frac{\kappa^{2}}{2}-\frac{2\gamma\omega_{0}^{2}}{\kappa}-\frac{2\omega_{0}^{2}\kappa}{\gamma},$
$\displaystyle d_{3}=$
$\displaystyle\Delta\left(2\gamma\kappa+4\Delta^{2}+2\gamma^{2}-8\omega_{0}^{2}+\kappa^{2}-\frac{2\gamma\omega_{0}^{2}}{\kappa}-\frac{2\omega_{0}^{2}\kappa}{\gamma}\right),$
$\displaystyle d_{4}=$
$\displaystyle\Delta^{2}\left(\Delta^{2}+\gamma^{2}+\gamma\kappa-2\omega_{0}^{2}+\frac{\kappa^{2}}{2}\right)+\frac{\gamma^{2}\kappa^{2}}{4}$
$\displaystyle+\gamma\omega_{0}^{2}\kappa+\frac{\gamma\kappa^{3}}{4}+\left(\omega_{0}^{2}+\frac{\kappa^{2}}{4}\right)^{2}.$

The key quantity is the static optomechanical frequency shift
$\mathcal{G}=\frac{2g_{0}^{2}}{\omega_{0}}|a|^{2}$, with $g_{0}=2\pi\times 1$
$\mathrm{Hz}$ the vacuum optomechanical coupling strength typical for these
membranes [62]. Additionally, $a=\left(\frac{\kappa
P_{\ell}/4\hbar\omega_{\ell}}{(\kappa/2+i\Delta)}\right)^{1/2}$ is the average
optical field amplitude in the cavity, determined by the cavity decay rate
$\kappa$, laser power and frequency $P_{\ell},\omega_{\ell}$ and detuning
$\Delta$. Finally, the mechanical frequency $\omega_{0}=2\pi\times 120$
$\mathrm{kHz}$ and linewidth $\gamma=2\pi\times 0.1$ $\mathrm{Hz}$ are rounded
from the values measured and reported in the main text.

To estimate $\kappa$, we calculate the finesse of the Fabry-Pérot cavity
formed by a $10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ distance between an
$R=0.3$ Si3N4 and an $R=0.35$ Si backplane. This is a low estimate for
$\kappa$, as the Si surface is rough from fabrication (Fig. 5d), and we do not
consider additional losses from that. We obtain $\kappa\sim 2\pi\times 430$
$\mathrm{GHz}$. We assume the optimal detuning to reach the frequency comb
regime, which is around $\Delta=-\kappa/4$. Combining these values, we can
evaluate the expressions for the stability and determine a minimum input laser
power to reach the frequency comb regime. This is at approximately
$1.5\text{\,}\mathrm{kW}$, several orders of magnitude larger than the
$3.6\text{\,}\mathrm{mW}$ measured laser power output of our laser Doppler
vibrometer. Thus the optomechanical (radiation pressure) instability is not
the origin of our frequency combs.

In conclusion, we have simulated the behavior of a Si3N4 membrane subject to
an incident laser. Due to the Si backplane, reflected light forms a standing
wave which leads to a counterpropagating optical trap. We have calculated the
force exerted on the membrane by the optical trap, its dielectrophoretic part
is sufficient to generate the overtone frequency comb as described in the main
text. We have also verified that the optical cavity formed between the Si3N4
membrane and Si backplane is not good enough to lead to an optomechanical
frequency comb based on the radiation pressure.

###  V.4 Optothermal parametric driving

Figure S8: Optical behavior of photonic crystal. a: Schematic of simulation
setup, consisting of a Si3N4 photonic crystal unit cell in air, with a plane
wave incident in the out-of-plane direction (red arrows). There are periodic
boundary conditions on the outer surfaces, with the brown markings indicating
one of the surface pairs. b: Norm of the E-field along the mid-surface of the
Si3N4 domain, showing a standing wave for $\lambda=630.8$ $\mathrm{nm}$. c:
Reflection and transmission of the photonic crystal as a function of
wavelength, and that of a flat Si3N4 surface for comparison (dotted lines).
The feature in the reflection/transmission is correlated with the mode of the
photonic crystal shown in b, indicated by the total E-field integrated over
the Si3N4 domain (green). This feature is also absent for an unpatterned
surface (green, dotted). d: SEM image of photonic crystal at the membrane
edge, with inset showing single hole.

Driving the trampoline resonators to sufficiently large amplitudes to see the
frequency comb from the optical trapping potential requires a strong driving
mechanism. We have previously used a piezo shaker to drive these membranes
[51] (also described in Sec. V.5), but for the majority of this work we use a
optothermal parametric drive mechanism where the membrane reaches the self-
oscillation regime [49]. That means that we can use a continuous-wave laser
without power or frequency modulation, and do not require any additional
signal or connection to the chip containing the resonator.

In the following subsections, we study the optothermal parametric driving.
Since it is based on absorption, which is normally not too large for
$\lambda=633$ $\mathrm{nm}$ in Si3N4, we first show the absorption is enhanced
by our photonic crystal structure. Then, we simulate the slow ($>0.1$
$\mathrm{s}$) thermal dynamics from steady state absorption and the fast ($\ll
0.1$ $\mathrm{s}$) thermal dynamics due to modulated intensity from the
standing wave optical field. Finally, we compare the simulated changes in
temperature to the optothermal self-oscillation theory [49] and other effects
from literature.

####  V.4.1 Absorption in the photonic crystal

Si3N4 is a often-used optical material because of its low absorption at
telecom wavelengths, with typical absorption coefficients $\mu=1-2$
${\mathrm{cm}}^{-1}$ at $1550\text{\,}\mathrm{nm}$ [63]. However at
$633\text{\,}\mathrm{nm}$, absorption is higher and sensitive to deposition
parameters (reports vary from $\mu=2$ ${\mathrm{cm}}^{-1}$ [64] to $\mu=5000$
${\mathrm{cm}}^{-1}$ [65]. Works using a similar setup but low-stress Si3N4
measured $0.5\text{\,}\mathrm{\char 37\relax}$ absorption for a
$50\text{\,}\mathrm{nm}$ film ($\mu=100$ ${\mathrm{cm}}^{-1}$) [66], so for
simulation purposes we assume $1\text{\,}\mathrm{\char 37\relax}$ for our
$80\text{\,}\mathrm{nm}$ film. It is difficult to verify the accuracy of this
assumption, especially since our trampoline resonators are patterned with an
array of circular holes that form a photonic crystal [67] and previously cited
absorption coefficients cover unpatterned films. These works suggest
absorption should be dependent mainly on material properties (Si3N4 deposition
parameters) thus be the same for all membranes on a chip. However, we
experience quite a variance in the response of our membranes to our laser beam
(SI Sec. V.7), and also depending on the exact positioning of the laser on the
membrane, so we propose a different effect.

Photonic crystal patterns can be designed to yield a certain reflectivity at a
target wavelength [67], but for the resonators in this work the photonic
crystal functions as a mechanism to accurately control the resonance frequency
[51], as well as allowing their smooth undercut. Due to this, the photonic
crystal parameters (lattice spacing $a=1356$ $\mathrm{nm}$, hole radius
$r=481$ $\mathrm{nm}$) are larger than the wavelength, which allows standing
waves in the Si3N4 domain that couple to waves incident on the plane of the
membrane. These modes enhance the field inside the Si3N4 and thus could
increase the absorption while being subject to fabrication imperfections such
that not all membranes or locations on membranes absorb equally, leading to
the variance in the observed optothermal self-oscillation.

We simulate the behavior of the photonic crystal in COMSOL by taking a unit
cell with the photonic crystal parameters ($a$, $r$, and Si3N4 thickness
$t=100$ $\mathrm{nm}$), as shown in Fig. S8a,d. We use $n=2.016$ for the Si3N4
domain and $n=1$ for vacuum. Parallel to the plane of the membrane we define
two periodic ports, for input (reflection) and output (transmission), which
are located sufficiently far away from the Si3N4 domain. These ports border a
perfectly-matched layer that absorbs scattered light that is not absorbed by
the periodic ports. The incident wave is an plane wave with normal incidence
to the slab, with a polarization in the $x+y$ direction, and the entire
simulation domain has periodic boundary conditions in $x$ and $y$ directions
(brown lines in Fig. S8a.

From this simulation, we find that a standing wave pattern forms in the plane
of the Si3N4 membrane, shown in Fig. S8b. When we evaluate the reflection and
transmission of this photonic crystal, we see only a small feature around
$\lambda=630.8$ $\mathrm{nm}$, Fig. S8c. However, the field integrated over
the Si3N4 domain is much enhanced. These modes are sensitive to material
thickness, $a$ and $r$, and not unique to the parameters chosen for our
simulation. We propose that standing waves of this form cause strong(er)
absorption of light in some of our structures, depending on fabrication
imperfections and precise laser position.

####  V.4.2 Slow thermal behavior

Figure S9: Thermal behavior of photonic crystal. a: Schematic of simulation of
the laser heating the suspended membrane. Boundary conditions and geometry are
described in the text. Results are extracted along two linecuts, horizontal
(left column) and diagonal (right column) in b, c, d. Starting from room
temperature, an incident heat flux representing absorption from the laser
light heats up the membrane over time. Colors indicate time, solid lines
assume thermal conductivity $k=$
$20\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$,
dashed lines use $k=$
$2\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$. b:
Temperature along the linecuts shows heating in the center, which settles
within $0.3\text{\,}\mathrm{s}$. c Displacement of the geometry away from
(positve) or towards (negative) the center of the domain due to thermal
expansion. The maximum displacement is $\sim 45$ $\mathrm{nm}$. d: Von Mises
stress (smoothed, only first and last timestep shown). Left and right panel
have different y-axis scale.

To build our understanding of the effects of the absorbed light, we simulate
the effect of a heat source on the photonic crystal membrane (Fig. S9a). The
membrane structure is symmetric, so we can reduce the simulation domain to
only a quarter of a membrane. The outer edges of the domain connect the
suspended Si3N4 connects to the (much more massive) Si substrate, so we apply
a fixed temperature condition (purple). The heating due to the laser light is
considered as an incident heat flux in the center of the membrane (orange),
defined via the total deposited heat as $P\mu_{\mathrm{Si_{3}N_{4}}}$ with
laser power $P=3.6$ $\mathrm{mW}$ and absorption coefficient
$\mu_{\mathrm{Si_{3}N_{4}}}=0.01$. The heat can be lost through the fixed
temperature boundary condition at the outer edges of the domain, but also
through radiative transfer in the out-of-plane direction, where we assume the
environment is at room temperature. We use a surface emissivity of $0.1$,
which is valid for room temperature and a film thickness of
$100\text{\,}\mathrm{nm}$ [68], and noticeably different than the $0.6-0.9$
value used for bulk Si3N4.

The choice of material parameters of Si3N4 is important for our simulation.
The Young’s Modulus $E=250$ $\mathrm{GPa}$, Poisson’s ratio $\nu=0.23$ and
density $\rho=3100$ $\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ are relatively
well-known, and the fabrication pre-stress is known to be
$1.0\text{\,}\mathrm{GPa}$. The thermal properties reported in literature vary
(see [66] and references therein), we take the specific heat $C_{p}=700$
$\mathrm{J}\text{\,}{\mathrm{kg}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ and thermal
expansion coefficient $\Upsilon=2.3\times 10^{-6}$ ${\mathrm{K}}^{-1}$.
Literature values for the thermal conductivity vary between $k=0.34$
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ [69] and
$k=20$ $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ [70].
We simulate using both values and plot the results as solid ($k=20$
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$) or dashed
lines ($k=2$ $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$)
in Fig. S9b,c,d.

The simulation results show that for both thermal conductivities, the
temperature (Fig. S9b) reaches a steady state after $\lesssim
0.3$$\mathrm{s}$, though the temperature is significantly higher for the low
thermal conductivity case ($470\text{\,}\mathrm{K}$) than for the high thermal
conductivity case ($349\text{\,}\mathrm{K}$). Due to this, the material
expands by 20 (45) $\mathrm{nm}$ for the high (low) thermal conductivity case
((Fig. S9c). On a lateral dimension of $150\text{\,}\mathrm{\SIUnitSymbolMicro
m}$ of the membrane this is negligible as it corresponds to approximately a
0.2 (0.4) $\mathrm{nm}$ increase in the lattice spacing $a$. Because of the
width of the pad (middle of membrane) with respect to the tether, the stress
in the pad is strongly reduced from the initial film stress
($1\text{\,}\mathrm{GPa}$), to $<20$ $\mathrm{MPa}$. By the heating, this can
be reduced to close to zero stress (compare yellow-dashed line in Fig. S9d to
the blue line). Overall, these changes are not so significant and unlikely to
change the dynamics of the membrane by themselves.

####  V.4.3 Fast thermal behavior

Figure S10: Fast thermal dynamics. a: Incident laser power (grey) and average
temperature difference with respect to steady state ($\simeq 348$ $\mathrm{K}$
for $120\text{\,}\mathrm{kHz}$ motion through standing wave optical trap
field. Colors represent integration over domains depicted in the insets. b:
Resulting stress change with respect to steady state.

After simulating the steady state in temperature, material deformation and
stress described in the previous subsection, we will now include the time-
modulated optical intensity. We assume that the resonator moves initially in
its fundamental mode at $120\text{\,}\mathrm{kHz}$, so the optical field
intensity is modulated by the $\sin$ term of Eq. 1. We make sure to keep the
time-averaged power the same as in our steady-state simulation. We simulate a
full period, and plot the incident heat flux (grey) and resulting temperature
averaged over circular domains at various distances from the center (colors)
in Fig. S10a.

The incident heat flux modulates the temperature of the membrane with an
amplitude of about $0.3\text{\,}\mathrm{K}$ and a delay on the order of a
microsecond. This delay is smallest in the domain where the heat flux is
directly incident (blue), and increases as we move further away from the
center (orange through red). The change in temperature also decreases as we
move further away, which indicates that the effect is local to the laser beam
spot. We can also extract the average von Mises stress, Fig. S10b, according
to
$\sigma_{\mathrm{avg}}=1/\sqrt{2}\sqrt{(\sigma_{xx}-\sigma_{yy})^{2}+6\sigma_{xy}}$,
with $\sigma_{xx,yy,xy}$ representing the stress components in the respective
coordinates. This shows a clear inverse correlation with the temperature, and
follows the trend of decreasing in amplitude and increasing in delay as we
move further away from the incident laser beam spot.

The resulting stress difference is relatively small compared to the total
local steady-state stress ($\sim 24$ $\mathrm{MPa}$). However, it can forms a
parametric driving mechanism since stress is modulated at double the frequency
of the original motion, guaranteed by the optical standing wave. This
optothermal parametric driving is thus automatically frequency-matched to the
dominant motion. The delay between heat flux and stress suggests that this
mechanism of driving is limited in our case to $\lesssim 500$ $\mathrm{kHz}$,
which corroborates with the fact that only observe combs based on the first
few eigenmodes.

####  V.4.4 Self-oscillation

Figure S11: Simulated optical power necessary to create a limit-cycle of
amplitude $R$, using the absorption expected for a bare Si3N4 membrane, and
with one enhanced by absorption through internal standing waves due to the
photonic crystal. Black dashed lines indicate maximum laser power used in
experiments, $3.6\text{\,}\mathrm{mW}$.

We have shown so far that the motion of the membrane through the optical
standing wave results in a modulated stress via the absorption. However, this
alone is not sufficient to drive the resonator to larger amplitudes since
dissipation is also present (though small). To estimate whether this effect is
strong enough to bring our resonator into limit-cycle oscillations under a
$P\lesssim 3$ $\mathrm{mW}$ drive, we follow the model of [49]. There, the
authors derive the minimum incident laser power to bring a resonator into the
self-oscillation regime (Hopf bifurcation),

$P_{\mathrm{Hopf}}=\frac{3(B^{2}+1)(B^{2}+4)}{2\pi^{2}AQ\psi K},$ (S15)

where

$\displaystyle K$
$\displaystyle=\eta_{1}\lambda^{4}R^{4}+\eta_{2}\lambda^{2}R^{2}+\eta_{3},$
(S16) $\displaystyle\eta_{1}$ $\displaystyle=-4\pi^{2}(B^{2}+1)C,$
$\displaystyle\eta_{2}$
$\displaystyle=-24\pi^{2}(B^{2}+1)x_{\mathrm{off}}^{2}C/\lambda^{2}-24\pi^{2}(B^{2}+4)x_{\mathrm{off}}D/\lambda$
$\displaystyle\quad+3(B^{2}+1)C,$ $\displaystyle\eta_{3}$
$\displaystyle=-4D(B^{2}+4)(8\pi^{2}x_{\mathrm{off}}^{2}/\lambda^{2}-3)x_{\mathrm{off}}/\lambda.$

Here, $A$ describes the heating due to the laser (units
$\mathrm{K}\text{\,}{\mathrm{W}}^{-1}$), $B$ is a dimensionless constant, $C$
represents the stiffness change due to temperature (units
${\mathrm{K}}^{-1}$), $D$ is the optothermal forcing term (units
${\mathrm{K}}^{-1}$), $Q$ is the mechanical Q-factor of our resonator, and
$\psi$ is the fraction of optical field that forms a standing wave. These
equations give the laser power necessary to create a limit-cycle of amplitude
$R$.

The authors of [49] give helpful guide on how to estimate the parameters of
Eq. (S15) from the simulations of Secs. V.4.3 and V.4.2. To estimate the
thermal parameters $A$ and $B$, we require the steady state temperature at the
point of illumination (center) of Fig. S9, $T_{\mathrm{dc}}$ (difference with
respect to the room temperature environment), and the amplitude of the
temperature change due to the standing wave optical field from Fig. S10,
$T_{\mathrm{ac}}$. Using

$\displaystyle A$ $\displaystyle=BT_{\mathrm{dc}}$ (S17) $\displaystyle B$
$\displaystyle=\frac{T_{\mathrm{ac}}\omega}{\sqrt{T_{\mathrm{dc}}^{2}-T_{\mathrm{ac}}^{2}}}$

and the values of $T_{\mathrm{dc}}$ and $T_{\mathrm{ac}}$ reported in the
previous subsections, we get $A=222425$ $\mathrm{K}\text{\,}{\mathrm{W}}^{-1}$
(independent of thermal conductivity $k$), and $B=1257$ ($k=2$
$\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$) or $B=3971$
($k=20$ $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$).

The forcing parameters $C$ (parametric) and $D$ (direct) can be estimated from
other simulations. By calculating the difference between eigenfrequencies of
our structure at the minimum and maximum temperatures of Fig. S10,
$\Delta\omega_{0}$ for the fundamental membrane mode, we can use
$C=2\Delta\omega_{0}/(\omega_{0}T_{\mathrm{ac}})$. We extract
$\Delta\omega=2\pi\times 200$ $\mathrm{Hz}$ from a simulation where we impose
only the stress change due to the thermal behavior, which gives $C\simeq 0.01$
${\mathrm{K}}^{-1}$ The parameter $D$ should be zero, because the tensile
stress means that a temperature change does not lead to a direct out-of-plane
displacement. This is true if the out-of-plane heat gradient is small enough,
the tensile stress does not reach zero, and as long as our entire structure
has the same thermal expansion coefficient (i.e. is from the same material).
The latter is what distinguishes us from previous work on the optothermal
excitation (bolometric backaction) driving a cantilever into self-oscillation
[71], where the gold layer induces deflection under temperature change such
that $D\neq 0$.

Finally, the we take mechanical quality factor $Q=1\times 10^{6}$ and
$\psi=1$. The latter we implicitly assumed already in Sec. V.4.3, where the
standing wave field cancels out fully (modulated power goes to zero). This is
likely an overestimation due to losses and transmission through the Si
backplane. The assumptions behind the derivation of [49] leading to Eq. (S15)
are somewhat different than for our system, since they use small displacement
$x\ll\lambda$.

We evaluate Eq. (S15) and plot the power necessary for a limit-cycle of radius
$R$ in Fig. S11. We distinguish between the absorption that we expect for an
unpatterned Si3N4 film (1% light absorption, blue), and if it would be
enhanced a factor 3 by the standing wave (orange). For non-zero limit cycle
amplitude $R$, the orange curve lies below the (maximum) power we send in,
$3.6\text{\,}\mathrm{mW}$ indicated by the black dashed line, which shows that
our membrane can be driven to self-oscillation by the optothermal parametric
effect. This requires a slightly higher absorption than would be expected
based on the material parameters, which can be achieved due to the internal
standing optical waves shown in the previous section. Fig. S11 shows that
limit cycles with amplitudes larger than $R\simeq 30$ $\mathrm{nm}$ can exist.
For $R\rightarrow 0$, no limit cycle appears which is due to $D=0$ in our
system. By evaluating the second derivative of $P$ with respect to $R$,
$\mathrm{d}^{2}P/\mathrm{d}R^{2}$, we find that the Hopf bifurcation is
supercritical and thus our limit-cycle is stable [49].

The optothermal parametric drive forms an effective periodic modulation, which
could also lead to a frequency comb, similar to the dielectrophoretic force.
By including a periodic variation in the frequency $\omega_{0}$, we simulate
this effect. The equation of motion is

$\ddot{x}+\gamma\dot{x}+\omega_{0}^{2}(1+F_{\mathrm{th}}\sin{2\omega_{0}t})x=0,$
(S18)

such that the parametric term oscillates at twice the frequency of motion.
With $F_{\mathrm{th}}=1\cdot 10^{-5}$, we are close to the simulated frequency
shift $\Delta\omega$. Numerically integrating Eq. (S18) yields an increase of
velocity over time (even in the presence of damping $\gamma/2\pi=0.2$
$\mathrm{Hz}$). However, the higher harmonics generated by the inclusion of
this parametric drive are negligible, $60\text{\,}\mathrm{dB}$ lower than the
fundamental mode. Higher harmonics (above $2\omega_{0}$) do not appear visibly
in the spectrum. While the optothermal parametric drive is sufficiently strong
to drive the membranes to self-oscillation, it is not strong enough to result
in a frequency comb without the dielectrophoretic force.

To conclude this section, we have simulated the optical and thermal behavior
of our resonators and analyzed the possibility of self-oscillation. We find
that the motion of the membrane through the standing wave optical field
modulates the absorbed light intensity, which in turn modulates the tensile
stress in the material. This is sufficiently strong to drive our resonator to
self-oscillation, if the absorption is slightly higher than what is expected
based on purely the material parameters.

###  V.5 Overtone comb with piezo driving

Figure S12: Comb with sine drive. Observed velocity power spectra for
different drive powers of a $104.524\text{\,}\mathrm{kHz}$ sine wave. The
overtones appear when the resonator is driven with high powers on resonance,
such that the displacement amplitude is large. The individual traces are
offset by $5\text{\,}\mathrm{dB}$ vertically.

In this section, we demonstrate that it is possible to generate the overtone
frequency comb through inertial driving with a piezo shaker. We do this by
using several membranes that do not show any sign of overtones when the laser
is positioned anywhere on the membrane for at least five minutes. This way, we
exclude thermal parametric driving. The piezo shaker is mounted on the back-
side of the sample holder and driven via an Agilent 33220A arbitrary waveform
generator and a a FLC A400 voltage amplifier with $20\times$ gain.

In Fig. S12, we show the velocity power spectrum of a device that is driven on
resonance by a sine wave of varying power. Without any drive, the fundamental
mode ($104.524\text{\,}\mathrm{kHz}$) can be observed, together with the
second and third modes above $250\text{\,}\mathrm{kHz}$. At higher powers, the
power in the first mode grows and overtones appear. For the highest drive
powers, we retrieve a comb similar to the one reported in the main text. This
demonstrates that using an external drive tone, we can utilize the overtone
mechanism to generate a frequency comb.

If the comb can be generated with the piezo by applying a white noise drive,
it would be possible to avoid an external frequency reference that must be
matched to the device. In Fig. S13a, we apply a white noise drive of varying
peak-to-peak amplitude. Again, the overtone comb is absent for lower drive
powers, but clearly present once the drive is sufficiently strong. Fig. S13b
shows the fundamental mode for the lowest and highest drive cases.

Figure S13: Overtone frequency comb through piezo driving. a: Displacement
power spectrum for various white noise driving powers. The overtone frequency
comb only appears for the highest driving powers. b: Fundamental mechanical
mode without and with driving. c: Displacement of the resonator as the drive
power is slowly increased from $0$ to $3\text{\,}\mathrm{V}$. The displacement
shows a less smooth signal than those obtained through thermal parametric
driving in the main text. d: Extracted strengths of the individual comb teeth
over time from the signal shown in c. Jumps around $25\text{\,}\mathrm{s}$
likely originate from the waveform generator internal switches. e: spectral
maps of fundamental and first two overtones. The jumps observed in d cause the
vertical lines.

We perform a similar analysis of the white-noise driven comb dynamics as for
the optothermal parametrically driven comb in the main text. The displacement
of the mechanical resonator over time as we increase the drive power is shown
in Fig. S13c, starting from $0\text{\,}\mathrm{V}$ at $t\leq 2.5$ $\mathrm{s}$
to $3\text{\,}\mathrm{V}$ at $t\geq 25$ $\mathrm{s}$. Compared to those
obtained via thermal parametric driving, the white noise drive results in a
less smooth displacement signal. The amplitude is also noticeably smaller,
which may be related to using a different displacement encoder, which is not
impedance-matched to the oscilloscope used to record the time signal. By
extracting the different overtones, Fig. S13d, we see the fundamental mode
increase in power first, the first and second overtone grow some seconds later
and the higher overtones only appear once the amplitude has grown
sufficiently. This matches the behavior of the comb shown in the main text.
The bumps between $t=20$ $\mathrm{s}$ and $t=28$ $\mathrm{s}$ are due to the
waveform generator internal switches.

We have shown that we can generate the overtone frequency comb using different
driving mechanisms. Using inertial driving via a piezoelectric shaker, we
observe qualitatively the same comb behavior as when using the thermal
parametric drive. However, the white noise drive results in a noisier
displacement signal and comb teeth powers. To highlight this, we show the
displacement power spectrum around the first three tones in Fig. S13e, where
the vertical lines indicate fast changes of comb power. This, combined with
the added complexity of requiring a piezoelectric shaker, voltage source and
amplifier, motivates the choice of using thermal parametric driving for the
majority of the measurements in this work.

###  V.6 Extension via comb interactions

Figure S14: Comb interactions. a: Spectrum around a higher-order mechanical
mode, showing a symmetrized copy of frequency comb of the fundamental mode,
$\omega_{0}$. The teeth are offset from the integer multiples of the original
mode (dashed lines), but still equally spaced by $\omega_{0}$. Inset shows
mode shape of the (5,5) mode expected to be the center peak at
$2083.3\text{\,}\mathrm{kHz}$. Readout position is identical to Fig. 2c of the
main text, where the (5,5) mode has considerable amplitude. b: The teeth of
the secondary comb follow the $9\text{\,}\mathrm{Hz}$ shift of the fundamental
comb mode as in Fig. 3, and shift symmetrically away from the center peak.

The bandwidth (span) of a frequency comb is an important property for many
applications. While the overtone comb already performs on-par with the largest
bandwidth mechanical combs (Table S1), we illustrate how the bandwidth can be
extended further. This can be done by letting the frequency comb interact with
a higher-order mechanical eigenmode $\omega_{\mathrm{h}}$ of the membrane, to
generate a comb with frequency spacing $\omega_{0}$ centered around
$\omega_{\mathrm{h}}$. This comb around $\omega_{\mathrm{h}}$ has harmonics
appearing on both sides of $\omega_{\mathrm{h}}$. In contrast, the comb of
$\omega_{0}$ has harmonics only on one side. This means that the interaction
with $\omega_{\mathrm{h}}$ could double the span.

We observe a copy of the fundamental mode overtone comb ($\omega_{0}$ spacing)
symmetrically around the higher-order mode of the system. In Fig. S14a, we
plot the spectrum around a higher order at $\omega_{\mathrm{h}}=$
$2083.3\text{\,}\mathrm{kHz}$ which we designate as the (5,5) mode shown in
the inset. The simulated frequency of the (5,5) mode is at
$2046\text{\,}\mathrm{kHz}$, and this is the closest mode with significant
mode amplitude at the readout position (same as in Fig. 2c of the main text).
The symmetric pattern of $\omega_{0}$-spaced comb teeth is offset from the
overtones of $\omega_{0}$ (black dashed lines).

To corroborate this interaction between the frequency comb as a whole and a
higher-order mode of the membrane, we closely study the frequency of two of
the comb teeth over time. Originating from the same measurement as Fig. 3 of
the main text, the fundamental mode at $\omega_{0}$ has a
$9\text{\,}\mathrm{Hz}$ upwards frequency shift over the duration of the
measurement. In Fig. S14b, we plot the spectrum of two of the comb teeth (the
first ones symmetrically around the central peak). These peaks show the same
shift in frequency, but one shifts upwards while the other shifts downwards in
frequency; they are clearly spaced around the central (5,5) mode peak by
$\omega_{0}$. The teeth further away from the center follow the $9n$
$\mathrm{Hz}$ scaling of the fundamental mode comb as described in the main
text. This interaction thus unlocks the possibility of shifting the overtone
comb upwards in frequency and thereby extend its span.

###  V.7 Effect of optics on comb and membrane

In this section, we demonstrate the overtone frequency comb on multiple
devices to corroborate the measurements shown in the main text. Additionally,
we also describe the visible changes in the membranes under an optical
microscope when driven into the overtone comb regime, and we study the effect
of the optics (power and focus) on the frequency comb.

####  V.7.1 Comb prevalence and visible deformation

Figure S15: Visual changes in membrane structure. a: Schematic sample maps of
two nominally identical chips containing $25$ membrane devices each, using $5$
different photonic crystal designs. Green checkmarks indicate which devices
show frequency comb behavior when optically addressed, blank squares show
membranes collapsed during fabrication. b: Microscope images of three
different devices before (left) and during (right) frequency comb behavior,
showing visible changes. The color pattern along the outer edge of the
membrane either disappears fully (top), partially (middle) or shows as a
slight discoloration centered around the laser spot (bottom), indicating
induced change in the photonic crystal lattice spacing and hole radius.

We observe the overtone frequency comb behavior on multiple devices, shown in
Fig. S15a. We fabricated two nominally identical chips containing $25$
membranes each. The membranes are identical with the exception of the photonic
crystal parameters, as described in earlier work [51]. All of the designs
feature the same internal modes that allow for strong absorption of
$\lambda=633$ $\mathrm{nm}$ light as in Sec. V.4, albeit at slightly different
wavelengths. Each chip contains the same set of $5$ photonic crystal designs,
oriented as marked in Fig. S15a.

Of the $45$ surviving devices on the two chips, $27$ demonstrated the
frequency comb behavior based on optothermal parametric driving (Sec. V.4),
and the remainder could be driven into the overtone comb regime by sufficient
piezoelectric shaking (Sec. V.5). This confirms the generality of the overtone
mechanism.

We can corroborate the relation of the overtone comb and (powerful) heating
through absorption by studying the membranes using the microscope of the laser
Doppler vibrometer. The presence of the optothermal parametrically driven
frequency comb behavior is correlated with visible changes in the membrane
structure. We show this by taking a microscope image immediately after the
laser spot is moved to the center of the membrane (Fig. S15b, left column),
and $>30$ $\mathrm{s}$ later when the system has reached steady state (Fig.
S15b, right column). There is an interference pattern at the edges of the
membrane, which is likely from the change in photonic crystal pattern due to a
stress-gradient from the edge of the Si3N4. The laser beam heats the membrane
and via the thermal expansion of the Si3N4, the stress changes such that this
edge deforms and the interference pattern changes.

Figure S16: Comb dependence on laser focus. a,c,e: Microscope images of the
same device with the laser focus before (a), at (c) and after (e) the plane of
the membrane structure. b,d,f: Frequency comb spectra measured with the
respective laser foci, showing no major change. g: Simulated trapping force as
a function of beam waist for constant total beam power.

We show three devices to illustrate different changes observable in this
interference pattern in Fig. S15b. In some devices (top row), the pattern
disappears completely, which suggests these have the strongest absorption: The
deformation is larger than in other devices and the steady-state temperature
is thus high. Simultaneously, these devices show strong frequency combs. In
other devices (Fig. S15b, middle row), the interference pattern disappears
only partially. This is typically along one diagonal, which suggests that the
stress along the other diagonal is not fully removed; we can switch between
the diagonals by positioning our laser spot. The steady-state temperature is
likely less than in the case where the pattern disappears completely, and also
the frequency comb appears less pronounced. Finally, a few devices show a
small discoloration around the laser spot, which suggests that the
interference pattern is modified only slightly. These devices display only
weak frequency combs. These behaviors corroborate that the thermal parametric
driving through absorption of the light is the origin of the frequency combs
that we observe.

####  V.7.2 Comb dependence on optical focus

We study the effect of laser focus on the frequency comb behavior. By
deliberately defocusing the laser beam, and using a membrane that has a weak
frequency comb behavior, we expect the focus to change the number of overtones
appearing in the spectrum.

In Fig. S16a,c,e, we show microscope images of the laser focus before, at and
after the plane of the membrane. In b,d,f, we plot the velocity power spectrum
measured from the devices with their respective foci, which are nearly
identical. This shows that the formation of the frequency comb does not depend
significantly on the laser focusing. Simulations of the optical trapping force
as a function of beam waist (Fig. S16g) confirm this behavior: The trapping
force does not change for beam waists $>3$ $\mathrm{\SIUnitSymbolMicro m}$.
Simultaneously, the thermal parametric force is limited to the area of the
beam spot, so as long as the total beam power is constant, the changes in
power density and beam area will cancel and the thermal parametric driving
will be constant. This shows that the generation of an overtone frequency comb
does not require a tightly focused laser beam.

####  V.7.3 Comb dependence on optical power

Figure S17: Comb dependence on optical power. Relative powers in the
individual overtones, normalized the overtone powers without filter and to the
power of the fundamental mode for each filter setting. With less optical
power, the overtones decrease in power with respect to the fundamental mode,
as simulated in Fig. S3. The labeled optical densities correspond to a
transmission of $0.1\mathrm{D}:79\%$, $0.2\mathrm{D}:63\%$,
$0.3\mathrm{D}:50\%$, $0.4\mathrm{D}:32\%$, $0.5\mathrm{D}:32\%$,
$0.6\mathrm{D}:25\%$, $1.0\mathrm{D}:10\%$.

We also study the effect of laser power on the frequency comb. We can add an
absorptive neutral density filter to reduce the power of the optical beam.
This decreased $F_{\mathrm{o}}$ and should decrease the power of the overtones
with respect to the fundamental mode, as simulated in Fig. S3. In Fig. S17, we
show the relative overtone powers measured on a single device, for various
filter strengths. The power in every overtone is normalized to the overtone
power for the measurement without optical filter. Subsequently, the powers of
the individual filter traces are then normalized for the power in their
respective fundamental mode (overtone number $0$). $F_{\mathrm{o}}$ affects
the power relative to the fundamental mode, not the absolute power.

From Fig. S17, we see a trend where the higher optical density filters lead to
lower relative overtone powers. This matches the simulations shown in Fig. S3.
For higher density filters and higher overtone numbers ($n\geq 4$), the signal
quality is too poor extract a meaningful relative overtone power. Nonetheless,
it is a clear validation of the model of Eq. (S2) as the correct description
of the source of the frequency combs.

###  V.8 Comb stability

Figure S18: Comb stability. a: Extracted power of different overtones
($\omega_{0}$, $5\omega_{0}$, $15\omega_{0}$ and $30\omega_{0}$, measured on
the same device for several timescales. There is no correlation in the
amplitudes on timescales longer than a minute, but there appears to be a
beating pattern on short timescales. b: Frequency of different overtones,
measured for several timescales. The frequency is stable on short timescales,
but drifts due to temperature on timescales longer than a minute. The drift is
the same (scaled by overtone number) for all overtones.

An important property of a frequency comb is the stability of the comb teeth
over time, both in amplitude and in frequency. To gain insight in how our
overtone frequency comb performs, we park our laser on a single device and
measure the spectrum on various timescales. From each spectrum, extract the
frequency and amplitude of some of the overtones, and plot them in Fig. S18.

In Fig. S18a, we plot the power of the overtones at $\omega_{0}$,
$5\omega_{0}$, $15\omega 0$ and $30\omega_{0}$, offset vertically for clarity.
On the timescales $>1$ $\mathrm{min}$, the amplitude of the individual
overtones fluctuates around a constant mean with a 95% confidence interval
between $2.0\text{\,}\mathrm{dB}$ and $2.5\text{\,}\mathrm{dB}$. On shorter
timescales, a correlation between the amplitudes is visible which can be
reproduced by simulations if decay is included.

In Fig. S18b, we show the frequency of the same set of overtones, where the
width of the displayed part of the spectrum is proportional to the overtone
number. The color indicates power and is the same for the bottom three rows,
but multiplied by a factor 10 for the top row for clarity. All comb teeth
display the same frequency shift over time, scaled by overtone number. This
correspondence is _exact_ , to within the resolution of our measurement. On
short timescales ($<16$ $\mathrm{s}$), the frequency of the fundamental mode
is stable to within the resolution of our measurement
($0.78\text{\,}\mathrm{Hz}$). However, we can use the exact frequency
correspondence and observe a slight shift of the overtone at $30\omega_{0}$,
$4.68\text{\,}\mathrm{Hz}$ over the measurement time, or
$7.9\text{\times}{10}^{-8}\text{\,}{\mathrm{s}}^{-1}$ relative frequency
stability. This is likely limited by the environment temperature causing a
frequency shift of the mechanical mode. That shift is clearly visible as the
slow drift on longer timescales, as the right-most column was measured
overnight and the frequency seems to stabilize. There are also oscillations on
a $10\text{\,}\mathrm{min}$ timescale, which we attribute to the room
airconditioning. If we select the most stable 6-hour period of our overnight
measurement, the fractional stability reaches
$7.5\text{\times}{10}^{-10}\text{\,}{\mathrm{s}}^{-1}$.

###  V.9 Noise in overtone frequency combs

We study fluctuations in the frequency comb as a response to (thermal) noise.
By starting from Eq. 1 and including a Langevin force noise term $\xi(t)$, we
obtain

$\ddot{x}+\gamma\dot{x}+\omega_{0}^{2}x=F_{\mathrm{o}}\sin\left(\frac{4\pi}{\lambda}(x-x_{\mathrm{off}})\right)+F_{\mathrm{d}}e^{i\omega_{d}t}+\xi(t),$
(S19)

Where the noise has zero mean $\langle\xi(t)\rangle=0$ and is correlated
$\langle\xi(t)\xi(t^{\prime})\rangle=2\gamma T\delta(t-t^{\prime})$ such that
we have assumed Markovian noise at temperature $T$, valid since resonator
$Q=\frac{\omega_{0}}{\gamma}\gg 1$ [72]. We expand the sine-squared nonlinear
term in the same way as in Sec. V.1, where we keep all the terms since the
displacement is not small. We can rewrite this to

$\ddot{x}+\gamma\dot{x}+\omega_{0}^{2}x=\sum_{n}^{\infty}\alpha_{n}x^{n}+F_{\mathrm{d}}e^{i\omega_{d}t}+\xi(t),$
(S20)

such that if $x_{\mathrm{off}}=0$, all odd coefficients
($\alpha_{1},~{}\alpha_{3},~{}\alpha_{5},$…) are zero. Let us first consider
the oscillation of only the first harmonic term,
$x=\frac{1}{2}\left(Ae^{i\omega_{0}t}+A^{*}e^{-i\omega_{0}t}\right)$. Then
expand the series of $\alpha_{n}x^{n}$ and truncate all the tones that are not
resonant. The first few terms read

$\displaystyle\alpha_{1}x$
$\displaystyle=\frac{\alpha_{1}}{2}\left(Ae^{i\omega_{d}t}+A^{*}e^{-i\omega_{d}t}\right)\Rightarrow\alpha_{1}x$
(S21) $\displaystyle\alpha_{2}x$
$\displaystyle=\frac{\alpha_{2}}{4}\left(A^{2}e^{2i\omega_{d}t}+2|A|^{2}+A^{*2}e^{-2i\omega_{d}t}\right)\Rightarrow
0$ $\displaystyle\alpha_{3}x$
$\displaystyle=\frac{\alpha_{3}}{8}\left(\begin{multlined}A^{3}e^{3i\omega_{d}t}+3|A|^{2}Ae^{i\omega_{d}t}\\\
+3|A|^{2}A^{*}e^{-i\omega_{d}t}+A^{*3}e^{-3i\omega_{d}t}\end{multlined}A^{3}e^{3i\omega_{d}t}+3|A|^{2}Ae^{i\omega_{d}t}\\\
+3|A|^{2}A^{*}e^{-i\omega_{d}t}+A^{*3}e^{-3i\omega_{d}t}\right)\Rightarrow\frac{3\alpha_{3}}{4}|A|^{2}x$
$\displaystyle\alpha_{4}x$
$\displaystyle=\frac{\alpha_{4}}{16}\left(\begin{multlined}A^{4}e^{4i\omega_{d}t}+4|A|^{2}A^{2}e^{2i\omega_{d}t}+6|A|^{4}\\\
+4|A|^{2}A^{*2}e^{-2i\omega_{d}t}+A^{*4}e^{-4i\omega_{d}t}\end{multlined}A^{4}e^{4i\omega_{d}t}+4|A|^{2}A^{2}e^{2i\omega_{d}t}+6|A|^{4}\\\
+4|A|^{2}A^{*2}e^{-2i\omega_{d}t}+A^{*4}e^{-4i\omega_{d}t}\right)\Rightarrow
0.$

We can thus rewrite Eq. (S20) as

$\ddot{x}+\gamma\dot{x}+\omega_{0}^{2}x=\left[\alpha_{1}x+0+\frac{3\alpha_{3}}{4}|A|^{2}x+0+\frac{5\alpha_{5}}{8}|A|^{4}x+...\right]\\\
+F_{\mathrm{d}}e^{i\omega_{d}t}+\xi(t).$ (S22)

Now we can introduce a new set of coefficients $\beta_{n}$ such that
$\beta_{1}=1,~{}\beta_{2}=0,~{}\beta_{3}=\frac{3}{4},~{}\beta_{4}=0,~{}\beta_{5}=\frac{5}{8}$.
Rearranging the terms now gives us

$\displaystyle\ddot{x}+\gamma\dot{x}+\left[1-\sum_{n}\alpha_{n}\beta_{n}|A|^{n-1}\right]\omega_{0}^{2}x$
$\displaystyle=F_{\mathrm{d}}e^{i\omega_{d}t}+\xi(t)$ (S23)
$\displaystyle\ddot{x}+\gamma\dot{x}+\Omega_{0}^{2}x$
$\displaystyle=F_{\mathrm{d}}e^{i\omega_{d}t}+\xi(t).$

This way, we work with the frequency of the first harmonic, $\Omega_{0}$,
taking into account the shift due to the non-linearity. In the case where
$x_{\mathrm{off}}=0$, all odd $\alpha_{n}=0$, and all even $\beta_{n}=0$ in
general, which implies that the frequency of the first harmonic is not shifted
by the non-linearity. Intuitively, that makes sense since the electric field
is symmetric around the resonator rest position for $x_{\mathrm{off}}=0$.

We obtain the amplitude of the first harmonic by switching off the noise
$\xi(t)$, and get

$|A|^{2}=\frac{F_{d}^{2}}{\left(\left[1-\sum_{n}\alpha_{n}\beta_{n}|A|^{n-1}\right]\omega_{0}^{2}-\omega_{d}^{2}\right)^{2}+\gamma^{2}\omega_{d}^{2}}.$
(S24)

Similarly, we can find the amplitude response to fluctuations $\xi(t)$ by
turning off the coherent drive,

$\ddot{x}+\gamma\dot{x}+\Omega_{0}^{2}x=\xi(t).$ (S25)

In Fourier space, we then get

$X(\omega)=\frac{\xi(\omega)}{\Omega_{0}^{2}-\omega^{2}+i\gamma\omega}$ (S26)

with autocorrelation

$\left\langle X(\omega)X(\omega^{\prime})\right\rangle=\frac{2\gamma
k_{\mathrm{B}}T\delta(\omega+\omega^{\prime})}{(\Omega_{0}^{2}-\omega^{2}+i\gamma\omega)(\Omega_{0}^{2}-\omega^{\prime
2}+i\gamma\omega^{\prime})}.$ (S27)

Then we can find the mean square amplitude response to fluctuations as

$\displaystyle\overline{\delta x^{2}}$ $\displaystyle=\left\langle
X(t)^{2}\right\rangle=\int\frac{\mathrm{d}\omega}{2\pi}\frac{\mathrm{d}\omega^{\prime}}{2\pi}\left\langle
X(\omega)X(\omega^{\prime})\right\rangle e^{i(\omega+\omega^{\prime})t}$ (S28)
$\displaystyle=\frac{2k_{\mathrm{B}}T}{\Omega_{0}^{2}}.$

Thus the mean square amplitude of the first harmonic due to thermal
fluctuations depends on the optical non-linearity via $\Omega_{0}^{2}$.

Now to obtain the net motion we consider both driving and noise terms (use Eq.
(S20)) and make the ansatz that our solution consists of a set of harmonics
with a perturbation $\delta x$

$x=\sum_{n=1}^{\infty}A_{n}e^{in\omega_{0}t}+\sum_{n=1}^{\infty}A_{n}^{*}e^{-in\omega_{0}t}+\delta
x=B+\delta x.$ (S29)

The nonlinear term $\sum_{n}\alpha_{n}x^{n}$ from Eq. (S20) can then be
expanded and simplified by neglecting higher-order noise terms and keeping
only the linear one. If $\delta x\ll B$, $(B+\delta x)^{n}\simeq
B^{n}+nB^{n-1}\delta x$. Now if we let the sum in the nonlinear term go to a
reasonable finite number $N_{\mathrm{max}}$, we split the sum and shift the
index of one of the two,

$\displaystyle\sum_{n=1}^{\infty}\alpha_{n}(B^{n}+nB^{n-1}\delta
x)=\sum_{n=1}^{N_{\mathrm{max}}}\alpha_{n}B^{n}+\sum_{n=0}^{N_{\mathrm{max}}-1}\alpha_{n+1}(n+1)B^{n}\delta
x$ (S30) $\displaystyle=\alpha_{1}\delta
x+\sum_{n=1}^{N_{\mathrm{max}}-1}B^{n}(\alpha_{n}+\alpha_{n+1}(n+1)\delta
x)+\alpha_{N_{\mathrm{max}}}B^{N_{\mathrm{max}}}$ .

If we then have $\alpha_{n}\simeq\alpha_{n+1}$ for all reasonable $n$, and
small noise such that $\delta(n+1)\ll 1$, we can then absorb the term at
$N_{\mathrm{max}}$ back into the sum,

$\sum_{n=1}^{\infty}\alpha_{n}(B^{n}+nB^{n-1}\delta x)\simeq\alpha_{1}\delta
x+\sum_{n=1}^{N_{\mathrm{max}}}\alpha_{n}B^{n}.$ (S31)

Thus we obtain the same solution in terms of the harmonics as we would do
without considering fluctuations $\delta x$, with only the additional term
$\alpha_{1}\delta x$. The equation of motion thus becomes

$\sum_{n=1}^{\infty}A_{n}\left(\omega_{0}^{2}-n^{2}\omega_{d}+in\gamma\omega_{d}\right)e^{in\omega_{d}t}+\ddot{\delta
x}+\gamma\dot{\delta x}+\omega_{0}^{2}\delta x=\alpha_{1}\delta
x+\sum_{n=1}^{N_{\mathrm{max}}}\alpha_{n}B^{n}+F_{d}e^{i\omega_{d}t}+\xi(t)$
(S32)

Based on the previous analysis, we can linearize this to

$\sum_{n=1}^{\infty}A_{n}\left(\omega_{0}^{2}-n^{2}\omega_{d}+in\gamma\omega_{d}\right)e^{in\omega_{d}t}+\ddot{\delta
x}+\gamma\dot{\delta x}+\Omega_{0}^{2}\delta x+F_{d}e^{i\omega_{d}t}+\xi(t)$
(S33)

where
$\Omega_{m}^{2}=\left[1-\sum_{n}\alpha_{n}\beta_{n}|A|^{n-1}\right]\omega_{0}^{2}$.
While the amplitudes $A_{n}$ are obtained in Sec. V.1, this analysis shows
that the root-mean square fluctuations stay constant at
$\frac{2k_{\mathrm{B}}T}{\Omega_{m}^{2}}$ for all values of $n$. This
indicates that the comb generation does not cause additional fluctuations or
noise.

Thus we have confirmed our observations of Fig. 3f in the main text: The
Lorentzian linewidth of our resonator fundamental mode is the same in the
thermal regime as it is when driven into the overtone comb, since the optical
non-linearity does not add noise to the system.

###  V.10 Phase-coherence

To show the phase-coherence of all tones within the comb, we take a closer
look at the time domain signal. When the overtones at
$2,3,4,...\times\omega_{0}$ have considerable amplitude, we can extract the
phases of each component separately from the shape of a single period in the
time domain. To do so, fit the sum of the first $n$ cosine terms with
amplitude coefficients $a_{1},...,a_{n}$ and phase offsets
$\phi_{1},...,\phi_{n}$ using

$\displaystyle A(t)$
$\displaystyle=a_{1}\cos(2\pi\omega_{0}t+\phi_{1})+a_{2}\cos(4\pi\omega_{0}t+\phi_{2})+...$
(S34) $\displaystyle+a_{16}\cos(32\pi\omega_{0}t+\phi_{n}).$

When all phase offsets are the same, $\phi_{1}=\phi_{2}=...=\phi_{n}$, the
overtone comb is phase-coherent.

To fit Eq.(S34), we extract the amplitudes of the first $n=16$ overtones from
the measurement of the whole time signal ($\sim 36$ $\mathrm{s}$). This we use
as an initial guess for a fit that optimizes them to the final amplitudes that
best describe a single period ($\sim 8$ $\mathrm{\SIUnitSymbolMicro s}$). We
calculate two curves, one where all phase offsets $\phi_{1},...,\phi_{16}$ are
the same (orange in Fig. S19a,b) and one where all the phase offsets are
random (black). From this, it is clear that all overtones have the same phase
offet and thus the overtone comb is phase-coherent.

Figure S19: Phase coherence of overtone comb. a,b: Single period (blue) of the
measured time signals of Fig. 4c, together with fits containing the sum of the
first 16 cosines (overtones, at $\omega_{0}$, …, $16\omega_{0}$). In orange,
the cosines all have the same phase offset, while in black the phase offset is
random. This shows the phase-coherence of our overtone comb.

## References

  * Cundiff and Ye [2003] S. T. Cundiff and J. Ye, Colloquium: Femtosecond optical frequency combs, Reviews of Modern Physics 75, 325 (2003).
  * Fortier and Baumann [2019] T. Fortier and E. Baumann, 20 years of developsments in optical frequency comb technology and applications, Communications Physics 2, 153 (2019).
  * Metcalf  _et al._ [2019] A. J. Metcalf, T. Anderson, C. F. Bender, S. Blakeslee, W. Brand, D. R. Carlson, W. D. Cochran, S. A. Diddams, M. Endl, C. Fredrick, S. Halverston, D. D. Hickstein, F. Hearty, J. Jennings, S. Kanodia, K. F. Kaplan, E. Levi, E. Lubar, S. Mahadevan, A. Monson, J. P. Ninan, C. Nitroy, S. Osterman, S. B. Papp, F. Quinlan, L. Ramsey, P. Robertson, A. Roy, C. Schwab, S. Sigurdsson, K. Srinivasan, G. Stefansson, D. A. Sterner, R. Terrien, A. Wolszczan, J. T. Wright, and G. Ycas, Stellar spectroscopy in the near-infrared with a laser frequency comb, Optica 6, 233 (2019).
  * Boulder Atomic Clock Optical Network Collaboration [2021] Boulder Atomic Clock Optical Network Collaboration, Frequency ratio measurements at 18-digit accuracy using an optical clock network, Nature 591, 564 (2021).
  * Rosenband  _et al._ [2007] T. Rosenband, D. B. Hume, P. O. Schmidt, C. W. Chou, A. Brusch, L. Lorini, W. H. Oskay, R. E. Drullinger, T. M. Fortier, J. E. Stalnaker, S. A. Diddams, W. C. Swann, N. R. Newbury, W. M. Itano, D. M. Wineland, and J. C. Bergquist, Frequency ratio of Al+ and Hg+ single-ion optical clocks; metrology at the 17th decimal place, Science 319, 1808 (2007).
  * Miri  _et al._ [2018] M.-A. Miri, G. D’Aguanno, and A. Alù, Optomechanical frequency combs, New Journal of Physics 20, 043013 (2018).
  * Mercadé  _et al._ [2020] L. Mercadé, L. L. Martín, A. Griol, D. NAvarro-Urrios, and A. Martínez, Microwave oscillator and frequency comb in a silicon optomechanical cavity with a full phononic bandgap, Nanophotonics 9, 3535 (2020).
  * Zhang  _et al._ [2021] J. Zhang, B. Peng, S. Kim, F. Monifi, X. Jiang, Y. Li, P. Yu, L. Liu, A. Alù, and L. Yang, Optomechanical dissipative solitons, Nature 600, 75 (2021).
  * Wu  _et al._ [2022] S. Wu, Y. Liu, Q. Liu, S.-P. Wang, Z. Chen, and T. Li, Hybridized frequency combs in multimode cavity electromechanical system, Physical Review Letters 128, 153901 (2022).
  * Han  _et al._ [2022] X. Han, C.-L. Zou, W. Fu, M. Xu, Y. Xu, and H. X. Tang, Superconducting cavity electromechanics: The realization of an acoustic frequency comb at microwave frequencies, Physical Review Letters 129, 107701 (2022).
  * Mahboob  _et al._ [2012] I. Mahboob, Q. Wilmart, K. Nishiguchi, A. Fujiwara, and H. Yamaguchi, Tuneable electromechanical comb generation, Applied Physics Letters 100, 113109 (2012).
  * Cao  _et al._ [2014] L. S. Cao, D. X. Qi, R. W. Peng, M. Wang, and P. Schmelcher, Phononic frequency combs through nonlinear resonances, Physical Review Letters 112, 075505 (2014).
  * Maksymov  _et al._ [2022] I. S. Maksymov, B. Q. H. Nguyen, A. Pototsky, and S. A. Suslov, Acoustic, phononic, Brillouin light scattering and Faraday wave based frequency combs: Physical foundations and applications, Sensors 22, 3921 (2022).
  * Erbe  _et al._ [2000] A. Erbe, H. Krömmer, A. Kraus, R. H. Blick, G. Corso, and K. Richter, Mechanical mixing in nonlinear nanomechanical resonators, Applied Physics Letters 77, 3102 (2000).
  * Mahboob  _et al._ [2016] I. Mahboob, R. Dupuy, K. Nishiguchi, A. Fujiwara, and H. Yamaguchi, Hopf and period-doubling bifurcations in an electromechanical resonator, Applied Physics Letters 109, 073101 (2016).
  * Seitner  _et al._ [2017] M. J. Seitner, M. Abdi, A. Ridolfo, M. J. Hartmann, and E. M. Weig, Parametric oscillation, frequency mixing, and injection locking of strongly coupled nanomechanical resonator modes, Physical Review Letters 118, 254301 (2017).
  * Ganesan  _et al._ [2017a] A. Ganesan, C. Do, and A. Seshia, Phononic frequency comb via intrinsic three-wave mixing, Physical Review Letters 118, 033903 (2017a).
  * Czaplewski  _et al._ [2018] D. A. Czaplewski, C. Chen, D. Lopez, O. Shoshani, A. M. Eriksson, S. Strachan, and S. M. Shaw, Bifurcation generated mechanical frequency comb, Physical Review Letters 121, 244302 (2018).
  * Park and Ansari [2019] M. Park and A. Ansari, Formation, evolution, and tuning of frequency combs in microelectromechanical resonators, Journal of Microelectromechanical Systems 28, 429 (2019).
  * Goryachev  _et al._ [2020] M. Goryachev, S. Galliou, and M. E. Tobar, Generation of ultralow power phononic combs, Physical Review Research 2, 023035 (2020).
  * Singh  _et al._ [2020] R. Singh, A. Sarkar, C. Guria, R. J. T. Nicholl, S. Chakraborty, K. I. Bolotin, and S. Ghosh, Giant tunable mechanical nonlinearity in graphene-silicon nitride hybrid resonator, Nano Letters 20, 4659 (2020).
  * Chiout  _et al._ [2021] A. Chiout, F. Correia, M.-Q. Zhao, A. T. C. Johnson, D. Pierucci, F. Oehler, A. Ouerghi, and J. Chaste, Multi-order phononic frequency comb generation within a MoS2 electromechanical resonators, Applied Physics Letters 119, 173102 (2021).
  * 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 152, 041019 (2022).
  * Batista and Lisboa de Souza [2020] A. A. Batista and A. A. Lisboa de Souza, Frequency-comb response of a parametrically driven Duffing oscillator to a small added ac excitation, Journal of Applied Physics 128, 244901 (2020).
  * Ganesan  _et al._ [2017b] A. Ganesan, C. Do, and A. Seshia, Frequency transitions in phononic four-wave mixing, Applied Physics Letters 111, 064101 (2017b).
  * Ganesan  _et al._ [2018] A. Ganesan, C. Do, and A. Seshia, Phononic frequency comb via three-mode parametric resonance, Applied Physics Letters 100, 021906 (2018).
  * Keşkekler  _et al._ [2022] A. Keşkekler, H. Arjmandi, P. G. Steeneken, and F. Alijani, Symmetry-breaking induced frequency combs in graphene resonators, Nano Letters 22, 6048 (2022).
  * Wu  _et al._ [2019] H. Wu, Z. Qian, H. Zhang, X. Xu, B. Xue, and J. Zhai, Precise underwater distance measurement by dual acoustic frequency combs, Annalen der Physik 531, 1900283 (2019).
  * Ganesan and Seshia [2019] A. Ganesan and A. Seshia, Resonance tracking in a micromechanical device using phononic frequency combs, Scientific Reports 9, 9452 (2019).
  * Wall  _et al._ [2020] W. S. Wall, R. L. Kubena, Y.-K. Yook, J. Koehl, and R. J. Joyce, Phase noise transfer in high-Q quartz phononic frequency combs, in  _2020 IEEE International Ultrasonics Symposium (IUS)_ (2020).
  * Aleman  _et al._ [2020] A. Aleman, S. Muralidhar, A. A. Awad, J. Åkerman, and D. Hanstorp, Frequency comb enhanced Brillouin microscopy, Optics Express 28, 29540 (2020).
  * Muralidhar  _et al._ [2020] S. Muralidhar, A. A. Awad, A. Alemán, R. Khymyn, M. Dvornik, D. Hanstorp, and J. Åkerman, Sustained coherent spin wave emission using frequency combs, Physical Review B 101, 224423 (2020).
  * Chu  _et al._ [2017] Y. Chu, P. Kharel, W. H. Renninger, L. D. Burkhart, L. Frunzio, P. T. Rakich, and R. J. Schoelkopf, Quantum acoustics with superconducting qubits, Science 358, 199 (2017).
  * Bochmann  _et al._ [2013] J. Bochmann, A. Vainsencher, D. D. Awschalom, and A. N. Cleland, Nanomechanical coupling between microwave and optical photons, Nature Physics 9, 712 (2013).
  * Sletten  _et al._ [2019] L. R. Sletten, B. A. Morres, J. J. Viennot, and K. W. Lehnert, Resolving phonon Fock states in a multimode cavity with a double-slit qubit, Physical Review X 9, 021056 (2019).
  * Hu and Kimble [1994] Z. Hu and H. J. Kimble, Observation of a single atom in a magneto-optical trap, Optics Letters 19, 1888 (1994).
  * Endres  _et al._ [2016] M. Endres, H. Bernien, A. Keesling, H. Levine, E. R. Anschuetz, A. Krajenbrink, C. Senko, V. Vuletic, M. Greiner, and M. D. Lukin, Atom-by-atom assembly of defect-free one-dimensional cold atom arrays, Science 354, 1024 (2016).
  * Dholakia and Čižmar [2011] K. Dholakia and T. Čižmar, Shaping the future of manipulation, Nature Photonics 5, 335 (2011).
  * Wang  _et al._ [1997] M. D. Wang, H. Yin, R. Landick, J. Gelles, and S. M. Block, Stretching DNA with optical tweezers, Biophysical Journal 72, 1335 (1997).
  * Ashkin and Dziedzic [1989] A. Ashkin and J. M. Dziedzic, Optical trapping and manipulation of single living cells using infra-red laser beams, Berichte der Bunsengesellschaft für physicaklisch Chemie 93, 254 (1989).
  * Guccione  _et al._ [2013] G. Guccione, M. Hosseini, S. Adlong, M. T. Johnsson, J. Hope, B. C. Buchler, and P. K. Lam, Scattering-free optical levitation of a cavity mirror, Physical Review Letters 111, 183001 (2013).
  * Bose  _et al._ [2017] S. Bose, A. Mazumdar, G. W. Morley, H. Ulbricht, M. Toroš, M. Paternostro, A. A. Geraci, P. F. Barker, M. S. Kim, and G. Milburn, Spin entanglement witness for quantum gravity, Physical Review Letters 119, 240401 (2017).
  * Delić  _et al._ [2020] U. Delić, M. Reisenbauer, K. Dare, D. Grass, V. Vuletić, N. Kiesel, and M. Aspelmeyer, Cooling of a levitated nanoparticle to the motional quantum ground state, Science 367, 892 (2020).
  * Ashkin [1970] A. Ashkin, Acceleration and trapping of particles by radiation pressure, Physical Review Letters 24, 156 (1970).
  * Zemánek  _et al._ [1999] P. Zemánek, A. Jonáš, L. Šrámek, and M. Liška, Optical trapping of nanoparticles and microparticles by a Gaussian standing wave, Optics Letters 24, 1448 (1999).
  * Ni  _et al._ [2012] K.-K. Ni, R. Norte, D. J. Wilson, J. D. Hood, D. E. Chang, O. Painter, and H. J. Kimble, Enhancement of mechanical Q factors by optical trapping, Physical Review Letters 108, 214302 (2012).
  * Yang  _et al._ [2019] F. Yang, F. Rochau, J. S. Huber, A. Brieussel, G. Rastelli, E. M. Weig, and E. Scheer, Spatial modulation of nonlinear flexural vibrations of membrane resonators, Physical Review Letters 122, 154301 (2019).
  * Norte  _et al._ [2016] R. A. Norte, J. P. Moura, and S. Gröblacher, Mechanical resonators for quantum optomechanics experiments at room temperature, Physical Review Letters 116, 147202 (2016).
  * Aubin  _et al._ [2004] K. Aubin, M. Zalalutdinov, T. Alan, R. B. Reichenbach, R. Rand, A. Zehnder, J. Parpia, and H. Craighead, Limit cycle oscillations in CW-laser driven MEMS, Journal of Microelectromechanical systems 13, 1018 (2004).
  * Høj  _et al._ [2021] D. Høj, F. Wang, W. Gao, U. B. Hoff, O. Sigmund, and U. L. Andersen, Ultra-coherent nanomechanical resonators based on inverse design, Nature Communications 12, 5766 (2021).