\section{Introduction}
Over the last quarter century, optical frequency combs have become key tools for metrology, timing and spectroscopy~\cite{Cundiff2003,Fortier2019}, and are indispensable in many laboratories around the world. Their use took off after the invention of the mode-locked laser based frequency combs~\cite{Jones2000}, which improved the achievable clock accuracy to one picosecond per day and was awarded the 2005 Nobel prize in physics~\cite{Hall2006,Hansch2006}. Optical frequency combs have revolutionized fields such as astronomy~\cite{Metcalf2019} and cosmology~\cite{BACON2021}, allowed tests of fundamental physics with atomic clocks~\cite{Rosenband2007,McGrew2018,BACON2021} and their use spans applications ranging from navigation~\cite{Coddington2009} to medical diagnostics~\cite{Thorpe2008} and atmospheric monitoring~\cite{Bjork2016}. To keep with their widespread use, frequency combs have spurred remarkable laser development with the aim of operating in different spectral regimes and at varying frequency resolutions. However, a central challenge has been to enable systems of lower size, weight, complexity and power usage~\cite{Delhaye2007,Faist2016,Manurkar2018,Yu2018,Zhang2021,Hu2021}. Towards this goal, a frequency comb that does not require precisely aligned optics or feedback controls would open novel regimes of miniaturization and versatility. Recent investigations have discovered phonon-based frequency combs~\cite{Mahboob2012,Cao2014,Mahboob2016,Czaplewski2018,Maksymov2022} as a mechanical analogue to optical frequency combs. These small-scale, alignment free devices could provide a natively on-chip alternative to their optical counterparts. Furthermore, they could extend sensing capabilities to optically opaque materials and form an ideal testbed for the dynamics underlying frequency combs based on bistability or parametric driving~\cite{Seitner2017,Ganesan2017,Ganesan2018}. Until now, these mechanical frequency combs have been limited by requiring high drive powers, carefully designed mode frequencies and engineered mechanical non-linearities. We will address all of these issues in this work.

In a different regime of breakthrough physics, optical trapping has allowed us to manipulate and control small particles ranging from single atoms~\cite{Hu1994,Endres2016} to micrometers~\cite{Dholakia2011} 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~\cite{Wang1997} to whole living cells~\cite{Ashkin1989}. Optically trapped ultra-cold atoms and levitated nanoparticles are perfect testbeds for fundamental physics involving gravity and mesoscopic quantum mechanics~\cite{Guccione2013,Bose2017,Delic2020}. The unique capability to optically manipulate particles has transformed our experimental reach in both fundamental physics and biology, earning it the 2018 Nobel prize in physics~\cite{Ashkin1970,Ashkin1986}.

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 them via a novel mechanism that enables a mechanical frequency comb without requiring feedback control, a frequency reference or precision optics. We observe that due to a weak counterpropagating optical trap~\cite{Ashkin1970,Zemanek1999}, strongly driven silicon nitride (Si$_3$N$_4$) membranes vibrate not only at their mechanical eigenfrequencies, but at perfect integer multiples of a single frequency, thus forming 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~\cite{Ni2012,Chang2012,Muller2015,Page2016}, but if the displacement becomes of the order of a quarter wavelength, it creates integer multiple copies ('overtones', see the supplementary information (SI) Sec.~\ref{Overtones}) of the original membrane motion. Strikingly, these displacement amplitudes are achieved through a driving mechanism based on the same standing optical field. Through an optothermal parametric drive, we can generate overtone frequency combs using only a single, unmodulated mW (continuous wave) laser beam and no additional electronics. We first describe these mechanisms, then demonstrate the existence and our understanding of the frequency combs and finally study the dynamical behavior and performance. Our work paves the way towards mechanical frequency combs without additional feedback schemes, external frequency references or precisely aligned optics.

\section{Overtone and driving mechanism}
\begin{figure}
\includegraphics[width = 0.5\textwidth]{Figure1.pdf}
\caption{\textbf{Schematic of overtone frequency comb. a:} The suspended Si$_3$N$_4$ 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). \textbf{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.}
\label{Schematiccomb}
\end{figure}
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=$~\SI{80}{\nano\meter} thick Si$_3$N$_4$ trampoline membrane~\cite{Norte2016} (SI Sec.~\ref{PhCmembranes}), shown schematically in Fig.\ref{Schematiccomb}\textbf{a}. It rests $\sim 10$~\si{\micro\meter} above a backplane formed by the silicon (Si) substrate, and a laser ($\lambda = 633$~\si{\nano\meter}, $p \lesssim 3$~\si{\milli\watt}) from a commercial Polytec MSA400 laser Doppler vibrometer is incident on the membrane. At this wavelength, the Si$_3$N$_4$ reflects $\lesssim 30\%$ of the light, and Si reflects about $35\%$. Part of the light thus forms a standing wave, Fig.~\ref{Schematiccomb}\textbf{a}, periodic in the direction parallel to the mechanical motion.

The dielectric Si$_3$N$_4$ experiences a dielectrophoretic force proportional to the optical intensity, similar to a particle in a counterpropagating-wave optical trap~\cite{Ashkin1970,Zemanek1999}, which we simulate in SI Sec.~\ref{Opticaltrap}. 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.~\ref{Schematiccomb}\textbf{b}), with an additional component from the optical field (red solid line). For small motion ($|x| \ll \lambda/4$, cyan line in Fig.~\ref{Schematiccomb}\textbf{b}), the optical potential functions as an additional spring which can increase the mechanical Q-factor~\cite{Ni2012}.

If the motion of the membrane becomes of the order of the optical potential period ($\lambda/2$), the overtone frequency comb emerges. This can be seen as follows: The restoring force of the mechanical potential switches sign twice per oscillation (arrow pairs in Fig.~\ref{Schematiccomb}\textbf{b}). 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. These thus form a frequency comb based on the dielectrophoretic force from with the optical field, rather than the on the radiation pressure force (\cite{Miri2018} and SI Sec.~\ref{Opticaltrap}). The optical origin of this mechanism makes it independent of the phononic dispersion relation within the material, i.e. we can use a perfectly linear mechanical resonator since the source of the comb is optical.

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\times10^{-12}$~\si{\kilo\gram}. We write the dielectrophoretic force as a proportionality constant $F_\mathrm{o}$ (units of force) times a term periodic in $x$ with the same $\lambda/2$ period as the optical intensity $I = E^2 \propto \sin^2(\frac{2\pi}{\lambda} (x-x_\mathrm{off}))$. This way we can move finite-size effects of the membrane and the radiation pressure force into $F_\mathrm{o}$, which we numerically evaluate in SI Sec.~\ref{Opticaltrap}. Strictly, the dielectrophoretic force should be proportional to the gradient of $I$, but the key factor impacting the dynamics, the periodicity, is not affected by this simplification. Finally, we add some resonant drive $F_\mathrm{d}$, to obtain the equation of motion
\begin{equation}
\ddot{x} + \gamma \dot{x} + \omega_0^2 x = \frac{F_\mathrm{o}}{m_\mathrm{eff}} \sin^2\left(\frac{2\pi}{\lambda}(x-x_\mathrm{off})\right) + \frac{F_\mathrm{d}}{m_\mathrm{eff}} e^{i\omega_0 t}.
\label{EOM}
\end{equation}
\begin{figure*}
\includegraphics[width = 1.0\textwidth]{Figure2.png}
\caption{\textbf{Overtone frequency combs. a}: Measured spectrum of a mechanical frequency comb of the fundamental mode ($\omega_0 = 2\pi\times 118.049$~\si{\kilo\hertz}) of a suspended Si$_3$N$_4$ membrane, consisting of 35 peaks spaced by $\omega_0$. The overtones are spectrally flat until $\sim 2200$~\si{\kilo\hertz}, after which their amplitude drops exponentially. Inset shows the location of the laser spot to generate and read out the comb. \textbf{b}: Simulated overtone comb spectrum, by integrating Eq.~\eqref{EOM}, matching the measurement of \textbf{a}. \textbf{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. \textbf{d,e} Mechanical frequency combs of the second mode (257 kHz) and third modes (359 kHz). Difference in the noise floor is due to different encoders used due to overloading. Insets show the laser position to drive and read out the different modes.}
\label{Frequencycomb}
\end{figure*}
The overtone mechanism and the resulting frequency comb is independent of the choice of drive mechanism (e.g. piezoelectric, electrostatic, thermal), as long as the resonator displacement is large enough. We stress that this displacement is much larger than typical for these micromechanical resonators, and intrinsically in the optically non-linear regime for any interferometric position measurement~\cite{Dolleman2017}. That is why we use a laser Doppler vibrometer, which is capable of resolving displacements much larger than the wavelength of the light used and is not affected by same non-linear behavior.

The spatially-periodic optical field allows us to drive the resonator through optothermal self-oscillation~\cite{Aubin2004}. As 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. We study this in detail in Sec.~\ref{Thermalparametricdriving} and use it for the remainder of this work. It allows us to generate overtone frequency combs without externally modulating the laser intensity, and without any additional connection or power to the chip. For completeness, we also verify that we can create an overtone comb using a piezoelectric shaker (see Sec.~\ref{Combwithpiezodriving}) thus showing that the overtone mechanism is independent from the (optothermal) self-oscillation.

\section{Overtone frequency comb}

We measure the overtone frequency comb for the fundamental mode $\omega_0 = 118.049$~\si{\kilo\hertz} of our membrane in Fig.~\ref{Frequencycomb}\textbf{a}. In the velocity power spectrum, we see a series of 35 peaks spanning the full detection bandwidth of the setup, spaced by $\omega_0$ (spacing has $2.1\times 10^{-7}$ relative uniformity, limited by spectral resolution). Through integration of the velocity signal, we verify the displacement $x \gg \lambda/4$. By numerically integrating Eq.~\eqref{EOM} we can reproduce this overtone comb, shown in Fig.~\ref{Frequencycomb}\textbf{b}, 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 the strength of our model: We introduce the non-linearity via the optical field, thus we can avoid introducing mechanical nonlinearities previously used to explain this behavior~\cite{Zhou2017,Yang2019,Yang2021}.
\begin{figure*}
\includegraphics[width = 1.0\textwidth]{Figure3.pdf}
\caption{\textbf{Comb dynamics. a}: Measured displacement of membrane, showing clear growth before reaching a plateau around $t = 24$~\si{s}. \textbf{b}: Measured spectrum of membrane motion close to fundamental mode and two of its overtones (5th and 15th). A thermal shift of \SI{9}{\hertz} of the fundamental mode is visible as a $5\times9$ and $15\times9$~\si{\hertz} shift in the overtones (white dashed lines are horizontal to guide the eye). Color scale is power (dB), yellow the highest and dark blue the lowest.  \textbf{c}: First six overtones increasing in power as the membrane displacement increases (markers) in \textbf{a}, with simulated dynamics (solid lines) based on integration of Eq.~\eqref{EOM} matching quantitatively to the highest four. \textbf{d}: Membrane displacement, showing growth from close to thermal regime to steady state. \textbf{e}: Overtone amplitudes extracted from the spectrum of measurement of \textbf{d}, showing all 35 overtones within the detection bandwidth (black, dotted lines are overtones with numbers between the labeled ones).  \textbf{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$.}
\label{Slowshift}
\end{figure*}

The overtone comb can be distinguished from the other mechanical eigenmodes of the membrane. By positioning the laser on a corner of the membrane, the optothermal drive is less effective (Fig.~\ref{Frequencycomb}\textbf{c}), so we can simultaneously detect the first few comb teeth, spaced by $\omega_0$, and the higher order eigenmodes (insets show mode shapes, absolute displacement, simulated). This demonstrates that the overtones are not eigenmodes of the mechanical resonator, and thus they are not affected by the dispersion relation of phonons in Si$_3$N$_4$ and do not require engineered non-linear resonances~\cite{Cao2014,Ganesan2017}. Furthermore, it is possible to generate combs from the second and third mechanical eigenmodes, by selecting the right laser position on the membrane, Fig.~\ref{Frequencycomb}\textbf{d,e}. This allows us to select a frequency spacing, limited only 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.~\ref{Seccombinteractions}).

We examine Eq.~\eqref{EOM} to better understand the behavior of this mechanism (details in SI Sec.~\ref{Overtones}). Firstly, when the displacement $x$ increases, more overtones appear due to the increasing number of optical extrema crossed. Each overtone will again be modulated by the optical potential, transferring energy between them and generating overtones of higher orders. This results in the overtones having equal power up to a certain cutoff $n \propto \frac{3x}{\lambda/4}$, $\sim 2200$~\si{\kilo\hertz} in Fig.~\ref{Frequencycomb}\textbf{a,b}, after which their strength drops exponentially. Secondly, the strength of the optical field ($F_\mathrm{o}$) controls the power of each overtone relative to the original mode $\omega_0$. If $F_\mathrm{o}$ is too large, the mechanical frequency and all overtones will shift by the induced non-linearity. Finally, the offset $x_\mathrm{off}$ between mechanical and optical zeros controls the relative intensity of the odd and even number overtones. The position offset $x_\mathrm{off}\simeq 40$~\si{\nano\meter} is consistent between different membranes in this work, and is half the thickness of the Si$_3$N$_4$ membranes. We study the overtone comb in different membranes, and the effect of the optical beam itself on the membranes and overtone comb in SI Sec.~\ref{Comblaser}.

\section{Comb dynamics}
In this section we quantitatively study the dynamical behavior (growth) of overtone frequency combs. In Fig.~\ref{Slowshift}\textbf{a}, we show a measurement where the laser starts at a position away from the membrane, and we move it to the membrane center at $t=$~\SI{8.7}{\second} such that the optothermal parametric driving increases the displacement. It drives the overtone frequency comb, which we analyze 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 the frequency behavior over time (Fig.~\ref{Slowshift}\textbf{b}) and monitor the power of individual overtones by taking linecuts (Fig.~\ref{Slowshift}\textbf{c}).

In Fig.~\ref{Slowshift}\textbf{b}, the fundamental mode shows a \SI{9}{\hertz} upwards shift in frequency, which is reproduced as an $n\times9$~\si{\hertz} upwards shift for the $n$th overtone (e.g. $n=5,15$ in the figure) and likely originates from slow thermalization (SI Sec.~\ref{Slowthermal}). This points to thermal tuning as an effective mechanism to control and tune the resonance frequency and thus comb spacing. Based on our frequency comb mechanism, thermal tuning would not compromise the comb uniformity.

We plot the power of the overtones as the comb grows in Fig.~\ref{Slowshift}\textbf{c}. We can simulate and reproduce this growth quantitatively by integrating Eq.~\eqref{EOM}. Comparison between the experimental data (markers) and simulation (solid lines) in Fig.~\ref{Slowshift}\textbf{c} shows good agreement with fit parameters [$F_\mathrm{o} = 66.7$~\si{\pico\newton}, $x_\mathrm{off} = 40$~\si{\nano\meter}, $x_0 = 5$~\si{\nano\meter}, $v_0 = 1$~\si{\nano\meter\per\second} and $F_\mathrm{d} = 31.9$~\si{\pico\newton}]. Most important for the overtones is $F_\mathrm{o}$, which is in excellent agreement with simulated force of \SI{60}{\pico\newton} (SI Sec.~\ref{Opticaltrap}). All extracted curves share a vertical offset to account for the total detection efficiency. This shows Eq.~\eqref{EOM} correctly describes the dynamics of the resonator and individual overtones when the comb is growing.

In Figs.~\ref{Slowshift}\textbf{d,e}, we observe the membrane motion from thermal state at $t<5$~\si{\second} to the steady state of the comb at $t>33$~\si{\second}. Initially, we detect only the fundamental mode $\omega_0$, until the displacement is sufficient to cross an optical extremum. After this, higher overtones appear as it the membrane crosses more extrema every oscillation. However, the displacement stops growing at $t\simeq 26$~\si{\second}, and decreases slightly before the system reaches steady state. 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, we study and quantify several properties relevant to the performance of the frequency comb, see SI. Sec~\ref{Combstability}.

We isolate and plot the spectrum around several of the overtones in Fig.~\ref{Slowshift}\textbf{f}. 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$~\si{\hertz} 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.~\ref{Phasenoisederivation}. 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-\si{\milli\hertz} linewidths based on ultra-high-Q membrane resonators~\cite{Hoj2021}.

We show the time-domain signal of the overtone comb in steady state in Fig.~\ref{Timedomaincomb}\textbf{a}. There is no component with periodicity longer than the mechanical period $2\pi/\omega_0$, which is distinctly different from soliton-based mechanical~\cite{Ganesan2017} and optomechanical~\cite{Zhang2021} frequency combs. This highlights the different physical and dynamical processes behind the overtone comb. The time-domain signal is always periodic with $\omega_0$ as shown in Figs.~\ref{Timedomaincomb}\textbf{b,c}. If the fundamental mode is the dominant component in the comb, we get a sinusoid (blue). When other components are dominant, particularly $2\omega_0$ and $3\omega_0$, we get the peaked curves (green, orange), 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.~\ref{Secphasecoherence}). This behavior is retained not only in the steady-state, but also during comb growth, Fig.~\ref{Timedomaincomb}\textbf{d}. There we plot the time signal at various stages during the measurement of Fig.~\ref{Slowshift}\textbf{d,e}, which shows a smooth transition from sinusoidal to peaked behavior as the comb grows. The periodicity and phase-coherence are retained.

\begin{figure}
\includegraphics[width = 0.5\textwidth]{Figure4.pdf}
\caption{\textbf{Time domain frequency comb. a}: Measured velocity of resonator in frequency comb regime, displaying no pattern with period longer than $1/\omega_0$. \textbf{b}: Zoom-in of time-domain signal in \textbf{a}, showing the \SI{118}{\kilo\hertz} fundamental mode dominates. \textbf{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. \textbf{d} Time domain signal of comb during growth, same measurement as Fig.~\ref{Slowshift}\textbf{e} 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.}
\label{Timedomaincomb}
\end{figure}

\section{Conclusion}
We have discovered an optomechanical mechanism which uniquely interfaces two breakthrough concepts: optical trapping and frequency combs. This allows for frequency combs whose simplicity stands out from previous paradigms. This frequency comb is realized by integrating a strongly driven suspended dielectric membrane in a weak optical trap. The dielectrophoretic force from the periodic optical field is modulated by the membrane's motion. This creates integer multiple copies (overtones) of the membrane's frequency, which form a frequency comb. We show combs of up to 35 overtones in a \SI{4.2}{\mega\hertz} 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. 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 principle, the mechanism behind these frequency combs is agnostic to optical wavelength, and requires only a suspended dielectric material. This gives us complete freedom in designing the mechanical properties of our resonators. By integrating ultra-high-Q membranes~\cite{Hoj2021} or choosing tunable mechanical frequencies, the possibilities of these frequency combs can be greatly enhanced. 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.

\subsection*{Acknowledgements}
We would like to acknowledge Farbod Alijani, Simon Gröblacher and Peter Steeneken 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., 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.

\subsection*{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. R.N. supervised the project. All authors contributed to writing and editing the manuscript.

\subsection*{Data availability statement}
All measurement data, analysis and calculation scripts, and simulations are available at \href{https://figshare.com/s/3266314a46fbed6c04c4}{https://figshare.com/s/3266314a46fbed6c04c4}, to be replaced with a DOI upon acceptance.

\clearpage

\setcounter{figure}{0}
\renewcommand{\thefigure}{S\arabic{figure}}
\setcounter{equation}{0}
\renewcommand{\theequation}{S\arabic{equation}}

\section{Supplementary information}
This part of the document contains the supplementary information

\subsection{Mechanical overtones}\label{Overtones}
\begin{figure*}
\includegraphics[width = 1.0\textwidth]{FigureS1.png}
\caption{\textbf{Mechanical overtones. a}: Resonator modes of a string. \textbf{b}: Frequency diagram of eigenmodes of a harmonic string and the anharmonic SiN trampoline membranes used in this work. \textbf{c}: SiN trampoline modes, grey box denotes degenerate modes. The color scale denotes normalized absolute displacement, from minimum (blue) to maximum (red).}
\label{Schematicovertones}
\end{figure*}
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 \emph{harmonic}. This is the case for a model string (Fig.~\ref{Schematicovertones}\textbf{a,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 \emph{anharmonic}. In Fig.~\ref{Schematicovertones}\textbf{c}, we have simulated and plotted the first six eigenmodes of the trampoline membrane resonator of the main text, and plotted their frequencies in Fig.~\ref{Schematicovertones}\textbf{b}.

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,
\begin{equation}
\ddot{x} + \gamma \dot{x} + \omega_0^2 x = 0.
\end{equation}
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 $\sin^2(x)$ (non-linear) term, 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 \textbf{overtones}. Others have detected these components of the mechanical motion before, and refer to them in a similar manner~\cite{Zhou2017,Yang2019,Yang2021}. They ascribe the origin of these overtones to be mechanical in nature, while we propose an optical origin.

\subsubsection{Small displacement}
To expand the description of these overtones, we consider the equation of motion as Eq.~\eqref{EOM} in the main text, copied below for convenience,
\begin{equation}
\ddot{x} + \gamma \dot{x} + \omega_0^2 x = F_\mathrm{o} \sin^2\left(\frac{2\pi}{\lambda}(x-x_\mathrm{off})\right) + F_\mathrm{d} e^{i\omega_0 t}.
\label{SIEOM}
\end{equation}
We use the power reduction formula $\sin^2(y) = \frac{1}{2} - \frac{1}{2}\cos(2y)$, and assume a small displacement $x$ such that we can use the Taylor expansion for the cosine, $\cos(2y) \simeq 1 - 2y^2$. In the experiments, the displacement is distinctly not small, but it is instructive to consider the situation nonetheless. We get
\begin{equation}
\ddot{x} + \gamma \dot{x} + \omega_0^2 x = F_\mathrm{o} \frac{4\pi^2}{\lambda^2}\left( x - x_\mathrm{off} \right)^2 + F_\mathrm{d} e^{i\omega_0 t}.
\label{SIsmallamplitudeTaylor}
\end{equation}
For simplicity, we can take the offset between mechanical and optical potentials $x_\mathrm{off} = 0$, to get
\begin{equation}
\ddot{x} + \gamma \dot{x} + \omega_0^2 x = F_\mathrm{o} \frac{4\pi^2}{\lambda^2}x^2 + F_\mathrm{d} e^{i\omega_0 t}.
\label{SIsmallamplitudeEOM}
\end{equation}
This equation admits solutions of the form
\begin{equation}
x = \sum_{n=1}^\infty A_n e^{in\omega_0 t},
\label{Solutionforx}
\end{equation}
which are the integer multiples of our original frequency $\omega_0$; the overtones. Substituting the solution Eq.~\eqref{Solutionforx} into Eq.~\eqref{SIsmallamplitudeEOM} 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
\begin{equation}
\begin{aligned}
A_1 &= i \frac{F_\mathrm{d}}{\gamma \omega_0}, \\
A_2 &= F_\mathrm{o} \frac{4\pi^2}{\lambda^2} \frac{A_1^2}{(-3\omega_0^2 - 2i\gamma\omega_0)}, \\
A_3 &= F_\mathrm{o} \frac{4\pi^2}{\lambda^2} \frac{2A_1 A_2}{(-8\omega_0^2 - 3i\gamma\omega_0)}, \\
A_n &= F_\mathrm{o} \frac{4\pi^2}{\lambda^2} \frac{\sum_{j=1}^\infty \sum_{k=1}^\infty  A_j A_k}{[(1-n^2)\omega_0^2-i\gamma n\omega_0]}, \quad \mathrm{for} j+k = n .
\end{aligned}
\label{SIAforsmallamp}
\end{equation}
The sequence of overtone amplitudes is monotonically decreasing (see Fig.~\ref{Overtoneamplitudes}): Every $A_n <  A_{n-1}$. The driving strength, mechanical frequency and decay rate determine $A_1$, which will be small unless we drive the resonator very strongly (which would violate the assumption of small displacement). This means that in the small-displacement case, the overtones will be negligible.

\subsubsection{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 for the cosine,
\begin{equation}
\cos(y) = \sum_{n=0}^\infty \frac{(-1)^{n} y^{2n}}{(2n)!}.
\end{equation}
Instead of Eq.~\eqref{SIsmallamplitudeTaylor}, we get
\begin{equation}
\ddot{x} + \gamma \dot{x} + \omega_0^2 x = \frac{1}{2} F_\mathrm{o}\sum_{n=1}^\infty \frac{(-1)^{n-1} \left(\frac{4\pi}{\lambda} x \right)^{2n}}{(2n)!} + F_\mathrm{d} e^{i\omega_0 t},
\end{equation}
where we have made the same simplification $x_\mathrm{off} = 0$, and taken care of the $\frac{1}{2}F_\mathrm{o}$ term by letting the sum start from $n=1$. We can still use the same ansatz of Eq.~\eqref{Solutionforx}, and get
\begin{widetext}
\begin{equation}
\sum_{n=1}^\infty A_n \left((1-n^2)\omega_0^2 -in\gamma\omega_0 \right)e^{in\omega_0 t} = \frac{F_\mathrm{o}}{2} \sum_{n=1}^\infty \left(\frac{(-1)^{n-1}}{(2n)!} \left(\frac{4\pi}{\lambda}\right)^{2n} \left(\prod_k^{2n} A_k \right) e^{i \sum_{k} k \omega_0 t}\right) + F_\mathrm{d} e^{i\omega_0 t}.
\label{SIAforbigamp}
\end{equation}
\end{widetext}
Again, to extract the amplitudes of the individual overtones $A_n$ we collect the terms with the right frequency $e^{in\omega_0 t}$. In the small-amplitude case, we only deal with products of two terms, e.g. for $n=4$ we have $A_1A_3$. Here, we have retained the higher terms in the Taylor expansion, so we also obtain products of four terms, e.g. for $n=4$ we have $A_1^4$. The relevant terms for each $A_n$ are those where the indices $k$ sum up to $n$, such that $e^{i\sum_k k \omega_0 t} = e^{i n\omega_0 t}$. The first overtones ($n = 1,2,3$) are the same as in Eq.~\eqref{SIAforsmallamp}, but for $n\geq 4$, there is a difference. The first few terms read
\begin{equation}
\begin{aligned}
A_4 &= \frac{\frac{4F_\mathrm{o}\pi^2}{\lambda^2} (2A_1A_3 + A_2^2) -\frac{8F_\mathrm{o}\pi^4}{3\lambda^4} A_1^4}{(-15\omega_0^2 - 4i\gamma\omega_0)}, \\
A_5 &= \frac{\frac{4F_\mathrm{o}\pi^2}{\lambda^2}(2A_1A_4 + 2A_2A_3)-\frac{8F_\mathrm{o}\pi^4}{3\lambda^4}4A_1^3A_2}{(-24\omega_0^2 - 5i\gamma\omega_0)}, \\
A_6 &= \frac{\splitfrac{\frac{4F_\mathrm{o}\pi^2}{\lambda^2}(2A_1A_5 + 2A_2A_4 + A_3^2)}
{-\frac{8F_\mathrm{o}\pi^4}{3\lambda^4}(4A_1^3A_3 + 6A_1^2A_2^2) + \frac{64F_\mathrm{o}\pi^6}{45\lambda^6}A_1^6}}{(-35\omega_0^2 - 6i\gamma\omega_0)}.
\end{aligned}
\end{equation}
It is difficult to analytically express all the combinations of even numbers of integers that sum up to $n$, so we numerically evaluate the overtone amplitudes in Fig.~\ref{Overtoneamplitudes}. For small displacement, the power in each subsequent overtone drops exponentially. However, for large displacement this is no longer true under the right combinations of $F_\mathrm{o}$ and $F_\mathrm{d}$, and the higher-order terms from the expansion cause strong overtones that do not decrease exponentially in amplitude with overtone number, and can thus form a frequency comb.

\begin{figure}
\includegraphics[width = 0.5\textwidth]{FigureS2.pdf}
\caption{\textbf{Overtone amplitudes.} Numerically evaluated overtone amplitudes $A_n$ in the small- and large-displacement case. For large displacement, the overtones are much more pronounced than in the small-displacement case.}
\label{Overtoneamplitudes}
\end{figure}

\subsubsection{Numerical simulation of overtones}
So far, we have derived analytically that the addition of an optical $\sin^2(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, this 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.

\begin{figure}
\includegraphics[width = 0.5\textwidth]{FigureS3.png}
\caption{\textbf{a}: Numerical simulation of overtones generated depending on initial amplitude. Parameters $F_\mathrm{o} = 100$~\si{\pico\newton}, $\omega_0 = 2\pi \times 118.049$~\si{\kilo\hertz}, $\gamma = 2\pi \times 0.02$~\si{\hertz}, $F_\mathrm{d} = 0$~\si{\newton}, $x_\mathrm{off} = 30$~\si{\nano\meter}. \textbf{b}: Simulation of overtones generated depending on optical force $F_\mathrm{o}$, for $x_o = 200$~\si{\nano\meter}. Other parameters identical to \textbf{a}. \textbf{c}: Simulation of overtones depending on offset $x_\mathrm{off}$. For $x_\mathrm{off} = 0$~\si{\nano\meter}, we predominantly drive the even overtones through the cosine-expansion, while for $x_\mathrm{off} = \lambda/8 = 79$~\si{\nano\meter} we predominantly drive the odd overtones through the sine-expansion. At $x_\mathrm{off} = 35$, odd and even modes are approximately equal in power. Individual traces are offset vertically, $x_0 = 1000$~\si{\nano\meter} and other parameters are identical to \textbf{a}.}
\label{fignumericalsim}
\end{figure}

First, we investigate the dependence of the frequency comb on displacement $x$. In Fig.~\ref{fignumericalsim}\textbf{a}, 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 (\SI{100}{\nano\meter}), 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$~\si{\nano\meter}, more overtones appear. The lower overtones have approximately equal amplitude (empirically until $\omega = 3 x_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 = 3 x_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.~\ref{fignumericalsim}\textbf{b}. For $F_\mathrm{o} = 0$~\si{\newton}, 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$~\si{\nano\newton} corresponds to approximately \SI{6}{\watt} 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$~\si{\nano\meter}, 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$~\si{\nano\meter}, 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$~\si{\nano\meter}, we see odd and even numbered modes with approximately equal power.

To summarize: The addition of an optical $\sin^2(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.

\subsection{Trampoline membranes and setup}\label{PhCmembranes}
\begin{figure}
\includegraphics[width = 0.5\textwidth]{FigureS4.pdf}
\caption{\textbf{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. \textbf{b}: SEM image of the suspended membrane, with nominal design parameters in white. \textbf{c}: Top-down view of corner of the membrane, highlighting the interface of the photonic crystal with the tether. \textbf{d}: Tilted view to show Si surface below \textbf{c}, to illustrate etch roughness and an etched imprint of the photonic crystal pattern (bottom right).}
\label{Setupandmembrane}
\end{figure}
The setup used in this work is shown in Fig.~\ref{Setupandmembrane}\textbf{a}. 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}$~\si{\milli\bar}, 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 Sec.~\ref{Thermalparametricdriving}. 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 the readout in place of the LDV decoder. Using the history function of this oscilloscope, we can chain 16 measurements of 20 million data points each, which allows \SI{38}{\second} of time signal at \SI{8.33}{\mega\hertz} sampling rate. However, this method for reading out the velocity comes at a cost of the calibrated readout that the LDV decoder offers.

The LDV outputs a voltage signal proportional to velocity or displacement, depending on the LDV encoder used. If using the LDV decoder, 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.

The membrane structures used in this work are fabricated out of \SI{100}{\nano\meter} thick stoichiometric low-pressure chemical vapor deposition Si$_3$N$_4$ on top of a \SI{1}{\milli\meter} 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 Si$_3$N$_4$. The membranes are then released using a second ICP etch, at \SI{-120}{\celsius} for \SI{30}{\second}, resulting in about \SI{10}{\micro\meter} of undercut and a final Si$_3$N$_4$ thickness of \SI{80}{\nano\meter}. An SEM image of a released trampoline is shown in Fig.~\ref{Setupandmembrane}\textbf{b}, with the nominal design parameters added in white. The resulting trampolines have well-characterized mechanical properties~\cite{Norte2016,deJong2021}, with fundamental mode frequencies at \SI{120}{\kilo\hertz} and Q-factors typically of 1 million.

The membranes are pattered with a periodic array of holes (Fig.~\ref{Setupandmembrane}\textbf{c}) that forms a photonic crystal, though it is not designed to have a high reflectivity at the operating wavelength of our LDV ($\lambda = 633$~\si{\nano\meter}). 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~\cite{deJong2021}, and possesses internal optical resonances that facilitate absorption which we detail in Sec.~\ref{Thermalparametricdriving}. Finally, the release etch imprints the photonic crystal pattern on the Si backplane, which roughens the surface. Fig.~\ref{Setupandmembrane}\textbf{d} shows the final Si surface roughness from the release etch (top left half) and the imprint below the periodic hole pattern (bottom right half).

\subsection{Membrane in an optical trap}\label{Opticaltrap}
\begin{figure*}
\includegraphics[width = \textwidth]{FigureS5.png}
\caption{\textbf{Simulation of optical trap. a:} Electric field norm of incident \SI{633}{\nano\meter} light on \SI{100}{\nano\meter} Si$_3$N$_4$ membrane suspended $\sim$\SI{6}{\micro\meter} above a Si backplane. The reflected light forms a standing wave pattern. Zoom shows center of the Si$_3$N$_4$ membrane. Bottom inset shows the simulation geometry schematically. \textbf{b:} Force exerted by the optical field on the Si$_3$N$_4$ domain as the position of the membrane is varied. \textbf{c:} Normalized strength of the trap for different Si$_3$N$_4$ thicknesses.}
\label{figOpticaltrap}
\end{figure*}
Optical traps can be used to confine small particles, as the light exerts a force proportional to the gradient of the optical intensity~\cite{Ashkin1970,Ashkin1986}. 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~\cite{Pohl1951}. In this section, we will simulate a laser incident on a dielectric (Si$_3$N$_4$) 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.

\subsubsection{Dielectrophoretic force}
We simulate the physical system consisting of a $t\simeq$~\SI{80}{\nano\meter} thin Si$_3$N$_4$ slab suspended in air above a much thicker Si slab, shown in Fig.~\ref{figOpticaltrap}\textbf{a}. 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.~\ref{figOpticaltrap}\textbf{a}), producing a periodic optical intensity. Combined with the fact that the membranes are much larger than the optical gradient laterally ($\sim 750\times750$~\si{\micro\meter\squared}), 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~\cite{Zemanek1999}.

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.~\ref{figOpticaltrap}\textbf{a}). The total force $F_\mathrm{de}$ induced by the optics on the dielectric Si$_3$N$_4$ structure can then be found by integrating this gradient over the membrane domain,
\begin{equation}
F_\mathrm{de} = \int_{V_{SiN}} \alpha \frac{\partial |E|^2}{\partial y} \mathrm{d}V,
\label{Fototal}
\end{equation}
with $\alpha$ the induced dipole moment. 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.~\eqref{EOM}. The expression for $F_\mathrm{de}$ depends on the position of the resonator with respect to the Si backplane, Fig.~\ref{figOpticaltrap}\textbf{b}, whereas $F_\mathrm{o}$ is a constant. $F_\mathrm{de}(x)$ in Fig.~\ref{figOpticaltrap}\textbf{b} 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.~\eqref{EOM} in the main text, $\sin^2(\frac{2\pi}{\lambda}(x-x_\mathrm{off})$, 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^2(\frac{2\pi}{\lambda}(x-x_\mathrm{off})$, and keep $F_\mathrm{o}$ as a constant equal to the maximum of $|F_\mathrm{de}(x)|$. Thus for a $t = $~\SI{80}{\nano\meter} membrane subject to our \SI{3.6}{\milli\watt} laser, we expect a maximum out-of-plane force in the order of \SI{60}{\pico\newton}. This confirms the values of $F_\mathrm{o}$ used in the numerical simulations in the main text and in Sec.~\ref{Overtones}, which show that the dynamics of our membrane can be significantly altered by the presence of the optical field.

The different curves of Fig.~\ref{figOpticaltrap}\textbf{b} 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.~\ref{figOpticaltrap}\textbf{c}), we see that it decreases for increasing Si$_3$N$_4$ 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}$.

\subsubsection{Radiation pressure force}
In an optical cavity, the interaction between radiation pressure and a mechanically compliant element can lead to optomechanical frequency combs~\cite{Miri2018}. 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.~\cite{Miri2018} as they analyze the stability of a mechanically compliant end-mirror optomechanically coupled to an optical cavity. Our Si$_3$N$_4$ 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 \cite{Miri2018}, which amounts to checking the Routh-Hurwitz stability criterion. For completion, the conditions for stability are
\begin{equation}
\begin{aligned}
3\mathcal{G}^2 + 4\Delta\mathcal{G} + \Delta^2 + \kappa^2/4 &>0, \\
\mathcal{G}^4 + d_1 \mathcal{G}^3 + d_2 \mathcal{G}^2 + d_3 \mathcal{G} + d_4 &> 0,
\end{aligned}
\end{equation}
with coefficients
\begin{equation}
\begin{aligned}
d_1 = &4\Delta, \\
d_2 = &\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}, \\
d_3 = &\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), \\
d_4 = &\Delta^2 \left( \Delta^2 + \gamma^2 +\gamma \kappa - 2 \omega_0^2 + \frac{\kappa^2}{2}  \right) + \frac{\gamma^2 \kappa^2}{4} \\
& + \gamma\omega_0^2 \kappa + \frac{\gamma\kappa^3}{4} + \left( \omega_0^2 + \frac{\kappa^2}{4} \right)^2.
\end{aligned}
\end{equation}
The key quantity is the static optomechanical frequency shift $\mathcal{G} = \frac{2 g_0^2}{\omega_0} |a|^2$, with $g_0 = 2\pi \times 1$~\si{\hertz} the vacuum optomechanical coupling strength typical for these membranes~\cite{deJong2020}. 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$~\si{\kilo\hertz} and linewidth $\gamma = 2\pi \times 0.1$~\si{\hertz} 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 \SI{10}{\micro\meter} distance between an $R = 0.3$ Si$_3$N$_4$ and an $R = 0.35$ Si backplane. This is a low estimate for $\kappa$, as the Si surface is rough from fabrication (Fig.~\ref{Setupandmembrane}\textbf{d}), and we do not consider additional losses from that. We obtain $\kappa \sim 2\pi\times 430$~\si{\giga\hertz}. 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 \SI{1.5}{\kilo\watt}, several orders of magnitude larger than the \SI{3.6}{\milli\watt} 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 Si$_3$N$_4$ 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 Si$_3$N$_4$ membrane and Si backplane is not good enough to lead to an optomechanical frequency comb based on the radiation pressure.

\subsection{Optothermal parametric driving}\label{Thermalparametricdriving}
\begin{figure*}
\includegraphics[width = \textwidth]{FigureS6.pdf}
\caption{\textbf{Optical behavior of photonic crystal. a:} Schematic of simulation setup, consisting of a Si$_3$N$_4$ 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. \textbf{b:} Norm of the E-field along the mid-surface of the Si$_3$N$_4$ domain, showing a standing wave for $\lambda = 630.8$~\si{\nano\meter}. \textbf{c:} Reflection and transmission of the photonic crystal as a function of wavelength, and that of a flat Si$_3$N$_4$ surface for comparison (dotted lines). The feature in the reflection/transmission is correlated with the mode of the photonic crystal shown in \textbf{b}, indicated by the total E-field integrated over the Si$_3$N$_4$ domain (green). This feature is also absent for an unpatterned surface (green, dotted). \textbf{d}: SEM image of photonic crystal at the membrane edge, with inset showing single hole.}
\label{PhCOptical}
\end{figure*}

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~\cite{deJong2021} (also described in Sec.~\ref{Combwithpiezodriving}), but for the majority of this work we use a optothermal parametric drive mechanism where the membrane reaches the self-oscillation regime~\cite{Aubin2004}. 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$~\si{\nano\meter} in Si$_3$N$_4$, we first how the absorption is enhanced by our photonic crystal structure. Then, we simulate the slow ($>0.1$~\si{\second}) thermal dynamics from steady state absorption and the fast ($\ll 0.1$~\si{\second}) 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~\cite{Aubin2004} and other effects from literature.

\subsubsection{Absorption in the photonic crystal}
Si$_3$N$_4$ is a often-used optical material because of its low absorption at telecom wavelengths, with typical absorption coefficients $\alpha = 1-2$~\si{\per\centi\meter} at \SI{1550}{\nano\meter}~\cite{Steinlechner2017}. However at \SI{633}{\nano\meter}, absorption is higher and sensitive to deposition parameters (reports vary from $\alpha = 2$~\si{\per\centi\meter}~\cite{Frigg2019} to $\alpha = 5000$~\si{\per\centi\meter}~\cite{Gorin2008}. Works using a similar setup but low-stress Si$_3$N$_4$ measured \SI{0.5}{\percent} absorption for a \SI{50}{\nano\meter} film ($\alpha = 100$~\si{\per\centi\meter})~\cite{Chien2018}, so for simulation purposes we assume \SI{1}{\percent} for our \SI{80}{\nano\meter} 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~\cite{Gaertner2018} and previously cited absorption coefficients cover unpatterned films. These works suggest absorption should be dependent mainly on material properties (Si$_3$N$_4$ 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.~\ref{Comblaser}), 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~\cite{Gaertner2018}, but for the resonators in this work the photonic crystal functions as a mechanism to accurately control the resonance frequency~\cite{deJong2021}, as well as allowing their smooth undercut. Due to this, the photonic crystal parameters (lattice spacing $a = 1356$~\si{\nano\meter}, hole radius $r = 481$~\si{\nano\meter}) are larger than the wavelength, which allows standing waves in the Si$_3$N$_4$ domain that couple to waves incident on the plane of the membrane. These modes enhance the field inside the Si$_3$N$_4$ 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 Si$_3$N$_4$ thickness $d = 100$~\si{\nano\meter}), as shown in Fig.~\ref{PhCOptical}\textbf{a,d}. We use $n=2.016$ for the Si$_3$N$_4$ 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 Si$_3$N$_4$ 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.~\ref{PhCOptical}\textbf{a}.

From this simulation, we find that a standing wave pattern forms in the plane of the Si$_3$N$_4$ membrane, shown in Fig.~\ref{PhCOptical}\textbf{b}. When we evaluate the reflection and transmission of this photonic crystal, we see only a small feature around $\lambda = 630.8$~\si{\nano\meter}, Fig.~\ref{PhCOptical}\textbf{c}. However, the field integrated over the Si$_3$N$_4$ 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.

\subsubsection{Slow thermal behavior}\label{Slowthermal}
\begin{figure}
\includegraphics[width = 0.5\textwidth]{FigureS7.pdf}
\caption{\textbf{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 \textbf{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 =$~\SI{20}{\watt\per\meter\per\kelvin}, dashed lines use $k =$ ~\SI{2}{\watt\per\meter\per\kelvin}. \textbf{b:} Temperature along the linecuts shows heating in the center, which settles within \SI{0.3}{\second}. \textbf{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$~\si{\nano\meter}. \textbf{d:} Von Mises stress (smoothed, only first and last timestep shown). Left and right panel have different y-axis scale.}
\label{PhCthermal}
\end{figure}
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.~\ref{PhCthermal}\textbf{a}). 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 Si$_3$N$_4$ 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\alpha_\mathrm{Si_3N_4}$ with laser power $P = 3.6$~\si{\milli\watt} and absorption coefficient $\alpha_\mathrm{Si_3N_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 $\lambda = 633$~\si{\nano\meter} and a film thickness of \SI{100}{\nano\meter}~\cite{Zhang2020}, and noticeably different than the $0.6-0.9$ value used for bulk Si$_3$N$_4$.

The choice of material parameters of Si$_3$N$_4$ is important for our simulation. The Young's Modulus $E = 250$~\si{\giga\pascal}, Poisson's ratio $\nu = 0.23$ and density $\rho = 3100$ \si{\kilo\gram\per\meter\cubed} are relatively well-known, and the fabrication pre-stress is known to be \SI{1.0}{\giga\pascal}. The thermal properties reported in literature vary (see~\cite{Chien2018} and references therein), we take the specific heat $C_p = 700$~\si{\joule\per\kilo\gram\per\kelvin} and thermal expansion coefficient $\alpha = 2.3\times 10^{-6}$~\si{\per\kelvin}. Literature values for the thermal conductivity vary between $k = 0.34$~\si{\watt\per\meter\per\kelvin}~\cite{Alam2012} and $k = 20$~\si{\watt\per\meter\per\kelvin}~\cite{Larsen2013}. We simulate using both values and plot the results as solid ($k = 20$~\si{\watt\per\meter\per\kelvin}) or dashed lines ($k = 2$~\si{\watt\per\meter\per\kelvin}) in Fig.~\ref{PhCthermal}\textbf{b,c,d}.

The simulation results show that for both thermal conductivities, the temperature (Fig.~\ref{PhCthermal}\textbf{b}) reaches a steady state after $\lesssim 0.3$\si{\second}, though the temperature is significantly higher for the low thermal conductivity case (\SI{470}{\kelvin}) than for the high thermal conductivity case (\SI{349}{\kelvin}). Due to this, the material expands by 20 (45)~\si{\nano\meter} for the high (low) thermal conductivity case ((Fig.~\ref{PhCthermal}\textbf{c}). On a lateral dimension of \SI{150}{\micro\meter} of the membrane this is negligible as it corresponds to approximately a 0.2 (0.4) \si{\nano\meter} 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 (\SI{1}{\giga\pascal}), to $<20$~\si{\mega\pascal}. By the heating, this can be reduced to close to zero stress (compare yellow-dashed line in Fig.~\ref{PhCthermal}\textbf{d} to the blue line). Overall, these changes are not so significant and unlikely to change the dynamics of the membrane by themselves.

\subsubsection{Fast thermal behavior}\label{Fastthermal}
\begin{figure}
\includegraphics[width = 0.5\textwidth]{FigureS8.pdf}
\caption{\textbf{Fast thermal dynamics. a:} Incident laser power (grey) and average temperature difference with respect to steady state ($\simeq 348$~\si{\kelvin} for \SI{120}{\kilo\hertz} motion through standing wave optical trap field. Colors represent integration over domains depicted in the insets. \textbf{b:} Resulting stress change with respect to steady state.}
\label{PhCthermalfast}
\end{figure}
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 \SI{120}{\kilo\hertz}, so the optical field intensity is modulated by the $\sin^2$ term of Eq.~\eqref{EOM}. 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.~\ref{PhCthermalfast}\textbf{a}.

The incident heat flux modulates the temperature of the membrane with an amplitude of about ~\SI{0.3}{\kelvin} 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.~\ref{PhCthermalfast}\textbf{b}, 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$~\si{\mega\pascal}). 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$~\si{\kilo\hertz}, which corroborates with the fact that only observe combs based on the first few eigenmodes.

\subsubsection{Self-oscillation}\label{Selfoscillation}
\begin{figure}
\includegraphics[width = 0.5\textwidth]{FigureS9.pdf}
\caption{Simulated optical power necessary to create a limit-cycle of amplitude $R$, using the absorption expected for a bare Si$_3$N$_4$ 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, \SI{3.6}{\milli\watt}.}
\label{fighopfbifurcation}
\end{figure}
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$~\si{\milli\watt} drive, we follow the model of \cite{Aubin2004}. There, the authors derive the minimum incident laser power to bring a resonator into the self-oscillation regime (Hopf bifurcation),
\begin{equation}
P_\mathrm{Hopf} = \frac{3(B^2+1)(B^2+4)}{2\pi^2 AQ\psi K},
\label{Hopfpower}
\end{equation}
where
\begin{equation}
\begin{aligned}
K &= \eta_1 \lambda^4 R^4 + \eta_2 \lambda^2 R^2 + \eta_3, \\
\eta_1 &= -4\pi^2(B^2+1)C, \\
\eta_2 &= -24\pi^2 (B^2+1)x_\mathrm{off}^2C/\lambda^2 - 24\pi^2(B^2+4)x_\mathrm{off}D/\lambda\\
&\quad+ 3(B^2+1)C, \\
\eta_3 &= -4D(B^2+4)(8\pi^2 x_\mathrm{off}^2/\lambda^2-3)x_\mathrm{off}/\lambda.
\end{aligned}
\end{equation}
Here, $A$ describes the heating due to the laser (units \si{\kelvin\per\watt}), $B$ is a dimensionless constant, $C$ represents the stiffness change due to temperature (units \si{\per\kelvin}), $D$ is the optothermal forcing term (units \si{\per\kelvin}), $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 \cite{Aubin2004} give helpful guide on how to estimate the parameters of Eq.~\eqref{Hopfpower} the simulations of Secs.~\ref{Fastthermal} and \ref{Slowthermal}. To estimate the thermal parameters $A$ and $B$, we require the steady state temperature at the point of illumination (center) of Fig.~\ref{PhCthermal}, $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.~\ref{PhCthermalfast}, $T_\mathrm{ac}$. Using
\begin{equation}
\begin{aligned}
A &= B T_\mathrm{dc} \\
B &= \frac{T_\mathrm{ac} \omega}{\sqrt{T_\mathrm{dc}^2 - T_\mathrm{ac}^2}}
\end{aligned}
\end{equation}
and the values of $T_\mathrm{dc}$ and $T_\mathrm{ac}$ reported in the previous subsections, we get $A = 222425$~\si{\kelvin\per\watt} (independent of thermal conductivity $k$), and $B = 1257$ ($k = 2$~\si{\watt\per\meter\per\kelvin}) or $B = 3971$ ($k = 20$~\si{\watt\per\meter\per\kelvin}).

The forcing parameters $C$ (parametric) and $D$ (direct) can be estimated from the simulations. By calculating the difference between eigenfrequencies of our structure at the minimum and maximum temperatures of Fig.~\ref{PhCthermalfast}, $\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$~\si{\hertz} from a simulation where we impose only the stress change due to the thermal behavior, which gives $C\simeq 0.01$~\si{\per\kelvin} 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~\cite{Metzger2008}, 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.~\ref{Fastthermal}, 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 \cite{Aubin2004} leading to Eq.~\eqref{Hopfpower} are somewhat different than for our system, since they use small displacement $x \ll \lambda$.

We evaluate Eq.~\eqref{Hopfpower} and plot the power necessary for a limit-cycle of radius $R$ in Fig.~\ref{fighopfbifurcation}. We distinguish between the absorption that we expect for an unpatterned Si$_3$N$_4$ 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, \SI{3.6}{\milli\watt} 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. 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~\cite{Aubin2004}.

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. This we show can happen due to our photonic crystal pattern.

\subsection{Overtone comb with piezo driving}\label{Combwithpiezodriving}
\begin{figure*}
\includegraphics[width = \textwidth]{FigureS10.pdf}
\caption{\textbf{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. \textbf{b:} Fundamental mechanical mode without and with driving. \textbf{c:} Displacement of the resonator as the drive power is slowly increased from $0$ to \SI{3}{\volt}. The displacement shows a less smooth signal than those obtained through thermal parametric driving in the main text. \textbf{d:} Extracted strengths of the individual comb teeth over time from the signal shown in \textbf{c}. Jumps around \SI{25}{\second} likely originate from the waveform generator internal switches. \textbf{e:} spectral maps of fundamental and first two overtones. The jumps observed in \textbf{d} cause the vertical lines.}
\label{figNoisedriving}
\end{figure*}
In this section, we demonstrate that it is possible to generate the overtone frequency comb through white-noise driving with a piezo shaker. To distinguish the mechanism that generates the comb from the thermal parametric driving, we move to a membrane that does not show any sign of the latter, for any laser position and for at least five minutes. We use the white-noise driving on an Agilent 33220A arbitrary waveform generator, and send it through a FLC A400 voltage amplifier with 20x gain. For various noise peak-peak amplitude settings on the waveform generator, we record the mechanical spectra using the laser Doppler vibrometer (Fig.~\ref{figNoisedriving}\textbf{a}). The frequency comb is absent for lower drive powers, but clearly present once the drive is sufficiently strong. Fig.~\ref{figNoisedriving}\textbf{b} shows the fundamental mode for the lowest and highest drive cases.

We record the displacement of the mechanical resonator over time as we increase the drive power in Fig.~\ref{figNoisedriving}\textbf{c}, starting from \SI{0}{\volt} at $t\leq 2.5$~\si{\second} to \SI{3}{\volt} at $t\geq 25$~\si{\second}. 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 is 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.~\ref{figNoisedriving}\textbf{d}, we see the fundamental mode increase in power first, the first \& 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, for the thermal parametric driving. The bumps between $t=20$~\si{\second} and $t=28$~\si{\second} are due to the waveform generator internal switches.

We have shown that we can generate the overtone frequency comb using different driving mechanisms. Using white noise applied to 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.~\ref{figNoisedriving}\textbf{e}, 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.

\subsection{Extension via comb interactions}\label{Seccombinteractions}
\begin{figure}
\includegraphics[width = 0.5\textwidth]{FigureS11.pdf}
\caption{\textbf{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 \SI{2083.3}{\kilo\hertz}. Readout position is identical to Fig.~\ref{Frequencycomb}\textbf{c} of the main text, where the (5,5) mode has considerable amplitude. \textbf{b}: The teeth of the secondary comb follow the \SI{9}{\hertz} shift of the fundamental comb mode as in Fig.~\ref{Slowshift}, and shift symmetrically away from the center peak.}
\label{Combinteractions}
\end{figure}
The frequency range over which a frequency comb exists is an important property for spectroscopic (frequency difference) measurements. While this range for mechanical frequency combs is much smaller than conventional (optical) frequency combs (e.g.~\cite{Delhaye2007}), our overtone frequency comb already offers a much larger span than conventional mechanical frequency combs~\cite{Mahboob2012}. Here, we show that we can extend the span further 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}$.

We observe a copy of the fundamental mode overtone comb ($\omega_0$ spacing) symmetrically around the higher-order mode of the system. In Fig.~\ref{Combinteractions}\textbf{a}, we plot the spectrum around a higher order at $\omega_\mathrm{h} =$~\SI{2083.3}{\kilo\hertz} which we designate as the (5,5) mode shown in the inset. The simulated frequency of the (5,5) mode is at \SI{2046}{\kilo\hertz}, and this is the closest mode with significant mode amplitude at the readout position (same as in Fig.~\ref{Frequencycomb}\textbf{c} 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.~\ref{Slowshift} of the main text, the fundamental mode at $\omega_0$ has a \SI{9}{\hertz} upwards frequency shift over the duration of the measurement. In Fig.~\ref{Combinteractions}\textbf{b}, 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$~\si{\hertz} 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.

\subsection{Effect of optics on comb and membrane}\label{Comblaser}

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.

\subsubsection{Comb prevalence and visible deformation}
\begin{figure}[!ht]
\includegraphics[width = 0.5\textwidth]{FigureS12.pdf}
\caption{\textbf{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. \textbf{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.}
\label{Structuralchange}
\end{figure}
We observe the overtone frequency comb behavior on multiple devices, shown in Fig.~\ref{Structuralchange}\textbf{a}. 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 \cite{deJong2021}. All of the designs feature the same internal modes that allow for strong absorption of $\lambda = 633$~\si{\nano\meter} light as in Sec.~\ref{Thermalparametricdriving}, albeit at slightly different wavelengths. Each chip contains the same set of $5$ photonic crystal designs, oriented as marked in Fig.~\ref{Structuralchange}\textbf{a}.

Of the $45$ surviving devices on the two chips, $27$ demonstrated the frequency comb behavior based on optothermal parametric driving (Sec.~\ref{Thermalparametricdriving}), and the remainder could be driven into the overtone comb regime by sufficient piezoelectric shaking (Sec.~\ref{Combwithpiezodriving}). 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.~\ref{Structuralchange}\textbf{b}, left column), and $>30$~\si{\second} later when the system has reached steady state (Fig.~\ref{Structuralchange}\textbf{b}, 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 Si$_3$N$_4$. The laser beam heats the membrane and via the thermal expansion of the Si$_3$N$_4$, the stress changes such that this edge deforms and the interference pattern changes.

\begin{figure}
\includegraphics[width = 0.5\textwidth]{FigureS13.pdf}
\caption{\textbf{Comb dependence on laser focus. a,c,e:} Microscope images of the same device with the laser focus before (\textbf{a}), at (\textbf{c}) and after (\textbf{e}) the plane of the membrane structure. \textbf{b,d,f:} Frequency comb spectra measured with the respective laser foci, showing no major change. \textbf{g:} Simulated trapping force as a function of beam waist for constant total beam power.}
\label{Focuschange}
\end{figure}

We show three devices to illustrate different changes observable in this interference pattern in Fig.~\ref{Structuralchange}\textbf{b}. 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.~\ref{Structuralchange}\textbf{b}, 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.

\subsubsection{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.~\ref{Focuschange}\textbf{a,c,e}, we show microscope images of the laser focus before, at and after the plane of the membrane. In \textbf{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.~\ref{Focuschange}\textbf{g}) confirm this behavior: The trapping force does not change for beam waists $> 3$~\si{\micro\meter}. 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.

\subsubsection{Comb dependence on optical power}
\begin{figure}
\includegraphics[width = 0.5\textwidth]{FigureS14.pdf}
\caption{\textbf{Comb dependence on optical power density. a,c,e:} Microscope images of the same device using three different objectives to focus the laser. \textbf{b,d,f:} Spectra showing the frequency comb for the different objectives. \textbf{g:} Simulated trapping force as a function of total beam power.}
\label{Objectivechange}
\end{figure}
We also study the effect of laser power on the frequency comb. In our setup, we have limited control over the laser power, but we can utilize the results of the previous subsection to our advantage: Since neither the optical trapping force nor the thermal parametric driving depend on the focus/beam waist, we can use a set of different objectives to apply different optical powers. We show microscope images obtained using the (magnification/numerical aperture) objectives (2x/0.55), (5x/0.14) and (20x/0.4) in Fig.~\ref{Objectivechange}\textbf{a,c,e} respectively. We measure the total optical power in focus after the objective using a Thorlabs PM100USB power meter.

The measured velocity power spectra for each of the three objectives are shown in Fig.~\ref{Objectivechange}\textbf{b,d,f}. The difference in the number of overtones is clearly visible, but the trend is unexpected. Based on the simulations, Fig.~\ref{Objectivechange}\textbf{g}, higher beam power results in a stronger trapping force and a stronger thermal parametric driving (assuming absorption linear with power). So the objective with the most optical power (5x/0.14, \textbf{c,d}) should result in the most and strongest overtones, but the opposite is true: it results in the least overtones.

An explanation for this behavior can be found by extending the simulations of Sec.~\ref{Overtones} to start at $x_0 = 0$ and including a driving force. For sufficiently strong $F_\mathrm{o}$, optical potential becomes strong enough to confine the resonator to a single well, and the number of comb lines is limited because the resonator cannot cross multiple optical extrema anymore. It is difficult to verify whether the membrane is truly confined by the optical potential, as it would require shifting of the optical potential extrema and simultaneous measurement of the position of the membrane. Confining macroscopic objects such as these membranes by an optical trap is an active field of research~\cite{Guccione2013,Marago2013,Ilic2019}, and this could provide an interesting opportunity for future research.

\subsection{Comb properties}\label{Combstability}
\begin{figure*}
\centering
\includegraphics[width = 0.65\textwidth]{FigureS15.pdf}
\caption{\textbf{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. \textbf{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.}
\label{Figcombstability}
\end{figure*}
In this section, we discuss several properties of the overtone frequency comb and compare them with current state-of-the-art mechanical frequency combs based on other mechanisms.

\subsubsection{Frequency spacing and bandwidth}
The frequency spacing of our comb is determined only by the mechanical eigenfrequency of the membrane that we choose, typically low-kHz to MHz. We have have demonstrated in-situ selectable mode-spacings of $118$, $257$ and \SI{359}{\kilo\hertz} in the same device. This spacing is five orders of magnitude larger than bistability-based mechanical combs~\cite{Mahboob2012,Mahboob2016}, and two orders of magnitude larger than parametric resonance combs \cite{Ganesan2018}. There is a large amount of freedom to design the membrane resonator to have the desired eigenfrequencies and thus frequency spacing, but there is an expected tradeoff between mechanical frequency and displacement (thus number of overtones achievable). For the optothermal parametric driving, there is also likely an upper frequency limit based on the time delay in optical absorption. This limits the optothermal parametric driving range to $\lesssim 500$~\si{\kilo\hertz} in our material platform. We have also demonstrated some minor tuning capabilities of the frequency spacing based on the thermal shift of the mode frequency.

The bandwidth of an overtone frequency comb is based on the combination of frequency spacing and the number of overtones. We report bandwidths of at least \SI{4.167}{\mega\hertz}, limited by the setup. This is two orders of magnitude higher than other mechanical frequency combs~\cite{Ganesan2018}, and can potentially be extended further by interactions with high-order mechanical modes(SI Sec.~\ref{Combinteractions}). Furthermore, overtone frequency combs are spectrally flat, which is unique among mechanical frequency combs, as other mechanisms result in a Lorentzian distribution~\cite{Mahboob2012,Ganesan2018} and soliton-based combs typically follow a hyperbolic secant around a central mode~\cite{Grelu2012}.

\subsubsection{Comb uniformity}
The uniformity of a frequency comb is of importance for optical frequency combs, and often requires dispersion engineering. For mechanical combs, this is less relevant dependent on the exact mechanism of generating the frequency comb; uniformity usually follows from the frequency relation between two tones (drive and mechanical eigenmode, or two drives, etc.). For the overtone mechanism we describe, the uniformity of our comb lines should in principle be perfect, i.e. we do not require dispersion engineering. We verify the relative uniformity of the frequency spacing to $2.1\times 10^{-7}$, limited by the resolution of our measurement. This does not yet improve on the best reported values for mechanical frequency combs, $10^{-9}$~\cite{Mahboob2012}. The uniformity of the frequency spacing is retained through thermal frequency shifts.

\subsubsection{Comb stability}
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.~\ref{Figcombstability}.

In Fig.~\ref{Figcombstability}\textbf{a}, we plot the power of the overtones at $\omega_0$, $5\omega_0$, $15 \omega0$ and $30\omega_0$, offset vertically for clarity. On the timescales $>1$~\si{\minute}, the amplitude of the individual overtones fluctuates around a constant mean with a 95\% confidence interval between \SI{2.0}{\decibel} and \SI{2.5}{\decibel}. On shorter timescales, a correlation between the amplitudes is visible which can be reproduced by simulations if decay is included.

In Fig.~\ref{Figcombstability}\textbf{b}, 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 \emph{exact}, to within the resolution of our measurement. On short timescales ($< 16$~\si{\second}), the frequency of the fundamental mode is stable to within the resolution of our measurement (\SI{0.78}{\hertz}). However, we can use the exact frequency correspondence and observe a slight shift of the overtone at $30\omega_0$, \SI{4.68}{\hertz} over the measurement time, or \SI{7.9e-8}{\per\second} 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~\SI{10}{\minute} 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 \SI{7.5e-10}{\per\second}. The spacing between the modes is uniform up to the measurement uncertainty of $2.1\times 10^{-7}$.

\subsection{Noise in overtone frequency combs}\label{Phasenoisederivation}
It is interesting to study fluctuations in the frequency comb as a response to (thermal) noise. We do this by starting from Eq.~\eqref{EOM} and including a Langevin force noise term $\xi(t)$,
\begin{equation}
\ddot{x} + \gamma \dot{x} + \omega_0^2 x = F_\mathrm{o} \sin^2\left(\frac{2\pi}{\lambda}(x-x_\mathrm{off})\right) + F_\mathrm{d} e^{i\omega_d t} + \xi(t),
\label{EOMnoise}
\end{equation}
Where the noise has zero mean $\langle \xi(t) \rangle = 0$ and is correlated $\langle \xi(t) \xi(t') \rangle = 2\gamma T \delta(t-t')$ such that we have assumed Markovian noise at temperature $T$, valid since resonator $Q = \frac{\gamma}{\omega_0} \gg 1$~\cite{Kac1981}. We expand the sine-squared nonlinear term in the same way as in Sec.~\ref{Overtones}, where we keep all the terms since the displacement is not small. We get
\begin{equation}
\ddot{x} + \gamma \dot{x} + \omega_0^2 x = \frac{F_\mathrm{o}}{2} \sum_{n=1}^\infty \frac{(-1)^{n-1} \left(\frac{4\pi}{\lambda} (x-x_\mathrm{off} \right)^{2n}}{(2n)!} + F_\mathrm{d} e^{i\omega_d t} + \xi(t).
\label{EOMnoiseexpanded}
\end{equation}
We can rewrite this to
\begin{equation}
\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),
\label{EOMnoisereduced}
\end{equation}
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( A e^{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
\begin{equation}
\begin{aligned}
\alpha_1 x &= \frac{\alpha_1}{2}\left( Ae^{i\omega_d t} + A^*e^{-i\omega_d t} \right)\Rightarrow \alpha_1 x \\
\alpha_2 x &= \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 \\
\alpha_3 x &= \frac{\alpha_3}{8}\left(\splitfrac{A^3 e^{3i\omega_d t} + 3|A|^2Ae^{i\omega_d t}}{+3|A|^2A^* e^{-i\omega_d t} + A^{*3}e^{-3i\omega_d t}} \right) \Rightarrow \frac{3\alpha_3}{4}|A|^2 x \\
\alpha_4 x &= \frac{\alpha_4}{16} \left( \splitfrac{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.
\end{aligned}
\end{equation}
We can thus rewrite Eq.~\eqref{EOMnoisereduced} as
\begin{multline}
\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).
\end{multline}
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
\begin{equation}
\begin{aligned}
\ddot{x} + \gamma \dot{x} + \left[1 - \sum_n \alpha_n \beta_n |A|^{n-1}\right]\omega_0^2 x &= F_\mathrm{d} e^{i\omega_d t} + \xi(t) \\
\ddot{x} + \gamma \dot{x} + \Omega_0^2 x &= F_\mathrm{d} e^{i\omega_d t} + \xi(t).
\end{aligned}
\end{equation}
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
\begin{equation}
|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}.
\end{equation}

Similarly, we can find the amplitude response to fluctuations $\xi(t)$ by turning off the coherent drive,
\begin{equation}
\ddot{x} + \gamma \dot{x} + \Omega_0^2 x = \xi(t).
\end{equation}
In Fourier space, we then get
\begin{equation}
X(\omega) = \frac{\xi(\omega)}{\Omega_0^2 - \omega^2 + i\gamma\omega}
\end{equation}
with autocorrelation
\begin{equation}
\left\langle X(\omega) X(\omega') \right\rangle = \frac{2\gamma k_\mathrm{B}T \delta(\omega + \omega')}{(\Omega_0^2 - \omega^2 + i\gamma\omega)(\Omega_0^2 - \omega'^2 + i\gamma\omega')}.
\end{equation}
Then we can find the mean square amplitude response to fluctuations as
\begin{equation}
\begin{aligned}
\overline{\delta x^2} &= \left\langle X(t)^2 \right\rangle = \int \frac{\mathrm{d}\omega}{2\pi} \frac{\mathrm{d}\omega'}{2\pi} \left\langle X(\omega) X(\omega') \right\rangle e^{i(\omega + \omega')t} \\
&= \frac{2k_\mathrm{B}T}{\Omega_0^2}.
\end{aligned}
\end{equation}
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.~\eqref{EOMnoisereduced}) and make the ansatz that our solution consists of a set of harmonics with a perturbation $\delta x$
\begin{equation}
x = \sum_{n=1}^\infty A_n e^{in\omega_0 t} + \sum_{n=1}^\infty A_n^* e^{-in\omega_0t} + \delta x = B + \delta x.
\end{equation}
The nonlinear term $\sum_n \alpha_n x^n$ from Eq.~\eqref{EOMnoisereduced} 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 + n B^{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,
\begin{equation}
\begin{aligned}
\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& \\
= \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}}&.
\end{aligned}
\end{equation}
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,
\begin{equation}
\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.
\end{equation}

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

\begin{widetext}
\begin{equation}
\sum_{n=1}^\infty A_n\left(\omega_0^2 - n^2 \omega_d + in\gamma\omega_d\right) e^{in\omega_dt} + \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)
\end{equation}
Based on the previous analysis, we can linearize this to
\begin{equation}
\sum_{n=1}^\infty A_n\left(\omega_0^2 - n^2 \omega_d + in\gamma\omega_d\right) e^{in\omega_dt} + \ddot{\delta x} + \gamma \dot{\delta x} + \Omega_0^2 \delta x + F_d e^{i\omega_d t} + \xi(t)
\end{equation}
\end{widetext}
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.~\ref{Overtones}, 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.~\ref{Slowshift}\textbf{f} 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.

\subsection{Phase-coherence}\label{Secphasecoherence}
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
\begin{equation}
\begin{aligned}
A(t) &= a_1\cos(2\pi\omega_0 t + \phi_1) + a_2\cos(2\pi\omega_0 t + \phi_2) + ... \\
 &+ a_{16}\cos(2\pi\omega_0 t + \phi_{n}).
\label{eqphasecoherence}
\end{aligned}
\end{equation}
When all phase offsets are the same, $\phi_1 = \phi_2 = ... = \phi_{n}$, the overtone comb is phase-coherent.

To fit Eq.\eqref{eqphasecoherence}, we extract the amplitudes of the first $n=16$ overtones from the measurement of the whole time signal ($\sim 36$~\si{\second}). 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$~\si{\micro\second}). We calculate two curves, one where all phase offsets $\phi_1, ..., \phi_{16}$ are the same (orange in Fig.~\ref{Figphasecoherence}\textbf{a,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.

\begin{figure}
\includegraphics[width = 0.5\textwidth]{FigureS16.pdf}
\caption{\textbf{Phase coherence of overtone comb. a,b}: Single period (blue) of the measured time signals of Fig.~\ref{Timedomaincomb}\textbf{c}, 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.}
\label{Figphasecoherence}
\end{figure}

\bibliographystyle{apsrev4-2}