Beyond linear coupling in microwave optomechanics

We explore the nonlinear dynamics of a microwave optomechanical system consisting of a drumhead nano-electro-mechanical resonator (NEMS) capacitively coupled to a microwave cavity. Experiments are performed under a strong microwave Stokes pumping which triggers mechanical self-sustained oscillations. We analyze the results in the framework of an extended nonlinear optomechanical theory, and demonstrate that quadratic and cubic coupling terms in the opto-mechanical Hamiltonian have to be considered. Quantitative agreement with the measurements is obtained considering only genuine geometrical nonlinearities: no thermo-optical instabilities are observed, in contrast with laser-driven systems. Based on these results, we describe a method to quantify nonlinear properties of microwave optomechanical systems. This method is clearly a new technique available in the quantum electro-mechanics toolbox, where higher-order coupling terms are proposed as a new resource for specific quantum schemes like quantum non-demolition (QND) measurements. We also find that the motion imprints a wide comb of extremely narrow peaks in the microwave output field, which could also be exploited in specific microwave-based measurements, potentially limited only by the quantum noise of the optical {\it and} the mechanical fields for a ground-state cooled NEMS device.


I. INTRODUCTION
Combining mechanical resonators with dimensions of order a micron or less with superconducting circuit elements has led to an exciting field of research exploring the quantum properties of nanoelectromechanical systems (NEMS) [1].Mechanical components can be efficiently coupled to superconducting qubits or integrated within optomechanical resonant cavities, providing a new resource for quantum device engineering [2][3][4][5].With a mechanical mode cooled to its quantum ground state, these NEMS circuits are also a new unique tool for fundamental experiments at the frontier of quantum mechanics [6,7].
In the context of quantum electronics, non-classical mechanical states can be used as a new support for quantum information storage and processing [6,[8][9][10][11][12].By engineering the coupling between photons and phonons (i.e.bath engineering), non-reciprocal microwave quantumlimited on-chip components are being developed [13][14][15][16]; in order to couple quantum processors through optical photons, photon converters are being built around a quantum-mechanical degree of freedom [17,18].
The capabilities offered by quantum NEMS devices are thus extremely rich, but are essentially all building on the linear parametric coupling between light and motion [19].For instance, using a properly detuned pump from the resonance cavity (i.e.Anti-Stokes sideband pumping), one can actively cool down a mechanical mode to its quantum ground state, or reversely (i.e.Stokes pumping) amplify the mechanical motion [12,[20][21][22][23][24][25].From twotone schemes, one can devise back-action evading (BAE) measurements [26,27] that enable to beat the standard quantum limit by measuring one quadrature while feeding the back-action noise of the detector to the other one.
Beyond the standard optical frequency pulling proportional to mechanical position x, higher-order couplings appear to be also desirable for specific realizations: from an x 2 coupling, one can measure the energy of the mechanics and build quantum non-demolition (QND) measurements [28,29], i.e. measuring an eigenstate of the Hamiltonian while not perturbing its evolution.Successful experimental implementations of such nonlinear couplings have been realized using optics, with membranein-the-middle configurations [30] and superfluid optomechanics [31].
At very strong Stokes sideband pumping power, the mechanical mode enters self-sustained oscillations [19].A rich multi-stable attractor diagram has been described theoretically [32], with specific phase-noise and amplitude-noise properties [33][34][35].The mechanical amplitude of motion becoming very large, this coherent state dynamics is extremely sensitive to all nonlinearities present in the system; this had been discussed already in Ref. [36].Experiments in this regime have been performed in the optical domain [37][38][39]; but the strong laser pump powers always produce dominant thermo-optical nonlinearities that require specific modeling [40].
In the present Article we report on experiments per-formed at low temperatures on a microwave optomechanical setup driven in the self-oscillating regime.The dynamics imprints a comb in the microwave spectrum, from which we can measure more than 10 coherent extremely narrow (width of order a few Hz) peaks.These could potentially be used for microwave-based readouts in quantum information processing and low temperature detectors (LTDs), where combs are usually synthesized for multiplexing purposes [41].An alternative technology proposed in the literature builds on the nonlinearity of superconducting quantum interference devices (SQUIDs) [42,43].In comparison, our optomechanical combs are competitive: the peak width is of the same order (resolution of 1 part in 10 8 ), the distribution of harmonics is much more homogeneous (equal spacing), and most importantly the amplitudes are very large (while the critical current of the SQUID junctions fixes a technological limitation).Ultimately, with a NEMS cooled to its ground state the stability of the output field would be limited only by quantum noise [35].
We demonstrate that the dynamics of the selfoscillating state is imprinted by genuine geometrical nonlinearities that can be fit, and we develop the full nonlinear theory giving the tools to extract nonlinear optomechnical coupling terms up to the third order ∝ x 3 .Building on these results, self-sustained oscillations in microwave optomechanical systems become a new tool enabling the experimental definition of the full nonlinear Hamiltonian at stake.This could be employed for instance in future quantum electronics circuits with specific schemes aiming at QND measurements [28,29].

II. EXPERIMENT
We employ a standard microwave optomechanical system [3,44] consisting of a microfabricated lumped microwave cavity resonator coupled to an aluminum drumhead NEMS [6] (see Fig. 1 a, right inset).The chip is installed into a commercial dilution cryostat with base temperature 10 mK, equipped with a high electron mobility transistor (HEMT) detection circuitry.The cryogenics, thermometry and measurement techniques have been described in Ref. [25].
The chip is designed for reflection measurements.The aluminum microwave cavity resonates at ω c /2π ≈ 6.8 GHz.The cavity displays a one-directional external coupling rate of κ ext /2π ≈ 2 MHz, and a total damping rate of κ tot /2π ≈ 4 MHz.We performed the experiment using the fundamental mode of the drum NEMS device which resonates around ω m /2π ≈ 6.7 MHz and exhibits a typical damping rate of about γ m /2π ≈ 150 Hz at 50 mK.Details on the geometry and measured parameters can be found in Appendix A.
The optomechanical coupling mechanism arises from momentum transfer between light (i.e.photons) and mechanics (i.e.phonons).In the standard case of a Fabry-Perot cavity, the displacement of the movable end mirror (the mechanical degree of freedom) modulates the resonant frequency of the cavity (the optical mode).Because of the retarded nature of the radiation pressure force when the laser light is detuned from the cavity frequency, this interaction gives rise to dynamical backaction allowing either active cooling or amplification of the mechanical motion.In this respect, our experiment is analogous to optics but shifted in the microwave domain [3]; the fundamental mode of the drumhead device corresponding to the movable mirror degree of freedom modulating the capacitance C(x) of the electrical circuit [19].The Brownian motion of this mode then imprints sidebands in the optical spectrum that we measure.The motion amplitude being very small, no extra nonlinearity has to be considered and the optical damping (when cooling) and anti-damping (when amplifying) observed are linear in applied power P in [19].This is used to calibrate the linear optomechanical interaction of our setup [25].We obtain a single photon-phonon coupling strength |g 0 |/2π ≈ 10 Hz.
Blue-detuned pumping at ω c + ∆ (with ∆ > 0) gives rise to downward scattering of photons, leading to the creation of phonons in the mechanical mode, hence enhancing the Stokes sideband.This is accompanied by a narrowing of the mechanical peak due to the antidamping backaction.At very strong powers, the total mechanical damping can thus be totally canceled: this is called the parametric instability.Above this threshold, the system enters into self-sustained oscillations, the amplitude of the mechanical motion being defined self-consistently [19].In this regime the mechanical amplitude of motion is so large (reaching several nanometers) that the mechanical sidebands are not limited to a couple of peaks: a full comb appears and can be measured (see Fig. 1  b).The peaks detected are not Lorentzian anymore, and their shape is defined by phase noise in the system [33].They are extremely narrow (only a few Hz wide at GHz frequencies), essentially equally-spaced (by ω m ) and of extremely high amplitude: they can even be detected without any HEMT pre-amplification.As well, all nonlinearities in the device will impact this complex optomechanical dynamics.
At milliKelvin temperatures, microwave photon heating arising from the absorption in dielectrics does not produce any thermal expansion: there are thus no thermooptical nonlinearities in our system, in strong contrast with devices actuated by laser beams where they dominate [37][38][39][40].However, the strong pump signal required to reach the threshold of the parametric instability does give rise to heating effects.This is carefully characterized and taken into account experimentally, see Appendix B. As such, the nonlinearities that prevail in microwavebased systems are of geometrical origin.
The experiment is performed in the mechanical selfinduced oscillation regime, by measuring the output microwave signal corresponding to the Stokes peak (at frequency ω c + δ), varying the detuning δ and the power P in of the input blue-detuned pump (at frequency ω c + ω m + δ).The measured photon flux is shown in Fig. 2 (a).As a comparison, the calculation based on the basic model of Refs.[32,33] is displayed in the bottom panel.The two plots are very similar, and display strikingly a bistable region at high powers and large positive detunings.However, calculation and theory do not match perfectly, which is expected: this has to be the signature of nonlinear effects which were neglected so far.
The region of the stability diagram which seems to be the most impacted by nonlinearities is precisely the hysteretic one (Fig. 2 beyond the dashed line).Therefore, in addition to the overall topography of the measured Corresponding calculated colormap from the basic theory described in Refs.[32,33] (no nonlinearities, g1 = g2 = 0 in Section III).The region on the right of the dashed line (high powers, positive detuning) is bistable and exists only when entering from the self-oscillating state (up-sweeping frequencies).The pink cross marks the minimum power necessary for self-oscillations, while the red cross corresponds to the position of the beginning of the hysteresis.The yellow cross marks the end of this bistable region (at same power).∆P and ∆δ are discussed in the text.
signal in the (P in , δ) space, we shall measure the importance of nonlinear features by reporting the position of the bistability in powers with respect to the beginning of the self-sustained region ∆P , and its width in detuning ∆δ (see Fig. 2 b).
The question that arises is thus: which nonlinearities need to be included in a quantitative model?One would immediately think about the Duffing effect in mechanical devices [46], and correspondingly to the Kerr effect for the microwave cavity [51].Both are discussed in Section IV and within the Appendices, but are not the dominant nonlinear features essentially because small frequency shifts have only a marginal impact on the optomechanical scheme itself.We have thus to consider nonlinearities in the coupling itself, that is higher-order derivatives in the Taylor expansion of the coupling capacitance C(x).

III. THEORY
We investigate the influence of nonlinear position coupling on the dynamics of self-sustained oscillations.The wide separation of time-scales together with weak coupling and damping allows for a self-consistent approach in which the mechanical amplitude is slowly changing [32][33][34][35][36].We start by writing a modified optomechanical Hamiltonian of the form (in the rotating frame of the drive optical field): where â and b are the photon and phonon annihilation operator, respectively.g 0 ∝ dC/dx is the usual linear single photon-phonon coupling strength while we introduce g 1 ∝ d 2 C/dx 2 and g 2 ∝ d 3 C/dx 3 , respectively the quadratic and cubic coupling strengths.This expansion order is necessary and sufficient to provide quantitative fits of the data, see Section IV.C(x) is the cavity mode total capacitance while x denotes the position collective degree of freedom of the first mechanical flexural mode (see Appendix A for details).Ĥγ represents the external baths coupling Hamiltonian and Ω 2 = κ ext P in /[ (ω c + ∆)] is the normalized driving term.
In this case the equations of motion for both operators take the following form: When the amplitudes of both fields are large enough, we can neglect quantum fluctuations and use the standard semiclassical approach.We write for both the optics â → α and the mechanics b → β, leading to: This system of coupled equations can be solved by means of the ansatz for β [33]: β c being related to a static deflection x c of the drum, and Be −iφ corresponding to the (complex valued) coherent motion.In the following, we shall neglect the β c term; it is indeed responsible only for a tiny frequency shift of the mechanical resonance, which impacts only marginally the thought limit cycle state.For the same reason we did not include the mechanical (Duffing) nonlinearity in Eq. ( 1), see Appendix C for a detailed discussion on these issues.
For convenience, we now introduce a shifted detuning ∆ = ∆ + g 1 B 2 , and two renormalized coupling parameters G = 2g 0 + 3B 2 g 2 and Ḡ = 2g 0 + 6B 2 g 2 .The optical amplitude equation takes now the form: The solution can be found via the mathematical transform described in Ref. [32].We define α = αe iΘ with: leading to the simpler dynamics equation: We now use the Jacobi-Anger expansion three times (on the three terms defining Θ): where: Here, J n is the Bessel function of the first kind.Fourier transforming Eq. ( 9), we write α(t) = n∈Z αn e inωt with: and hence: where h n = i(nω − ∆ ) + κ tot /2, and we used the change of variable q = n − n in the penultimate line.Inserting Eq. ( 14) into Eq.( 5) and preserving only terms oscillating at −ω (rotating wave approximation), we obtain the dynamics equation for the mechanics: We can now recast this expression introducing the optical backaction terms, namely the optical spring term δω and the damping term γ BA : with ω = ω m + δω now explicitly defined, and: where X is written as: One can thus find all the stable states by solving selfconsistently the equation which cancels the effective damping γ m + γ BA , ensuring that: In practice, it is sufficient to solve the limit-cycle equation neglecting all kinds of mechanical shifts, assuming ω = ω m in Eq. (11).Details on the self-consistent determination of optomechanical stable states can be found in Appendix D. Following the same procedure as for α, the optical field amplitude in the cavity takes the form α = n∈Z α n e inωt with: This expression highlights the fact that the optomechanical coupling imprints a comb structure in the photon field.We can then compute the output photon flux Ṅout,n of each comb peak n as: Detuning (MHz) where we made use of the well-known input-output relation linking intra-cavity fields and output traveling fields [47].n = 0 corresponds to the pump tone at frequency ω c + ω m , n = −1 to the Stokes sideband at ω c and n = 1 to the anti-Stokes sideband at ω c + 2ω m (see black dots in Fig. 1 b).

IV. ROLE OF GEOMETRIC NONLINEARITIES
The aim is thus now to go beyond Fig. 2, and obtain quantitative agreement between theory and experiment.The theory in the previous section allows us to calculate the amplitude of the mechanical motion including geometrical nonlinearities in the couplings.However, to obtain estimates of the mechanical frequency we need to also include important contributions from other effects, especially the Duffing nonlinearity of the drum.
As soon as the system self-oscillates, the actual cavity frequency is slightly renormalized in ω c = ω c − g 1 B 2 .Besides, there is also a material-dependent shift with a logarithmic power-dependence that is attributed to Two-Level-Systems present in the dielectrics [49], which is taken into account (Appendix B).On the other hand, the cavity Kerr nonlinearity ξ c is expected to be extremely small for our device [50,51]; we give an upper bound in Tab.I, see discussions in Appendices B and C for details.The mechanical resonance is also renormalized by the optomechanical coupling, with a tiny frequency shift δω (see Section III).However, the dominant source of mechanical frequency shift is due to the Duffing effect (i.e. the mechanical nonlinearity arising from the stretching of the drum [48]), leading to ω m = ω m + δω + ξ m B 2 , with a normalized Duffing parameter ξ m in Hz per phonon.For simplicity, we will omit the prime on ω c and ω m , remembering in the following that the measured mechanical frequency shift includes all terms.
The measured output photon flux is plotted in Fig. 3 as a function of detuning δ and power P in (same data as Fig. 2 top panel, 214 mK).The amplitude of the signal is extremely large, but the most striking feature is the bistable region at high powers and positive detunings.The measured mechanical frequency shift is shown in Fig. 4; strikingly, we find that it is largest in the bistable regime.
This mechanical shift cannot be captured by the optomechanical contribution δω alone.One has to take into account the Duffing effect to quantitatively fit it (see below).However, the mechanical frequency shifts remain very small (a few kHz at most, see Fig. 4); we thus verified that they have only a marginal impact on the limit cycle dynamics (i.e. the amplitudes, B), see Appendix C.
In the hysteretic region the amplitude B becomes very large, hence the optomechanical response becomes sensitive to the nonlinear coupling coefficients, g 1 and g 2 .For symmetry reasons (see Appendix A), the sign of the g 0 parameter is irrelevant and we take it to be positive for simplicity.However then, the sign of the other coefficients is uniquely defined.
To calculate the amplitudes of the limit cycles (and hence the photon flux) using the approach in Sec.III, the only free parameters are the quadratic and cubic nonlinear coupling terms g 1 and g 2 , respectively.These two coefficients have a different impact on the calculated flux: around our best fit parameters, g 1 narrows/broadens the self-oscillating region with respect to detuning (altering the ∆δ parameter), while g 2 mostly shifts the bistable feature to higher/lower powers (∆P parameter).This is represented in Fig. 5; in Appendix E, full colormaps calculated for different nonlinear parameters are also displayed (Fig. 8).Indeed at the same time the overall shape of the theoretical maps displayed in Fig. 3  and Fig. 4 (mechanical frequency) are very sensitive to the nonlinear parameters.We can therefore reasonably well determine the values of these two terms, typically within a factor of 2 (see Appendix E).The theoretical fits are displayed as green dots in Figs. 3 and 4; as a comparison the colormap of Fig. 2 bottom panel is computed for g 1 = g 2 = 0. We performed this procedure at various cryostat temperatures.However, because of microwave absorption in the materials, the drum temperature is not homogeneous over the complete measured range of (δ, P in ).This effect is taken into account, see Appendix B. The most constrained point for the definition of the couple (g 1 , g 2 ) is the junction between the main stable region and the bistable part, defined by the red cross mark in Fig. 2. We shall thus define an effective temperature T eff characteristic of the fit at this precise point.
From the measured mechanical frequency shift (Fig.  4), we can finally fit the Duffing term ξ m .The summary of our results is given in Tab.I. Within our error bars, we can infer a unique set of parameters that fits all temperatures.This is a strong evidence that the nonlinear features g 1 , g 2 and ξ m are of geometrical origin.We give in Appendix A theoretical estimates obtained from basic arguments: a circular plate stretching nonlinearity for ξ m [48] and a corresponding plate-capacitor nonlinear expansion for g 1 , g 2 [52].The magnitudes match our findings within typically a factor of 2, apart from g 2 which prediction is the worst because of the crudeness of the plate capacitor analytic expansion.

V. CONCLUSION
We report on microwave optomechanical experiments performed in the self-sustained oscillation regime.The output spectrum of a microwave cavity resonating around 6.8 GHz coupled to a 6.7 MHz drumhead mechanical device is measured as a function of input power P in , pump frequency detuning δ and temperature.A high amplitude and narrow-peak comb structure is measured in the output spectrum, and fit to theory.
We demonstrate that the limit cycle dynamics is sensitive to nonlinearities in the optomechanical coupling.We therefore present a theory that goes beyond the standard linear optomechanical Hamiltonian, introducing quadratic and cubic terms g 1 and g 2 .Data is fit quantitatively, and we show that these g 1 and g 2 must be of geometrical origin, as opposed to the thermo-optical nonlinear features present in laser driven systems.
The work described here can thus be proposed as a new method to characterize nonlinearities in microwave nanomechanical platforms.With the development of new quantum-limited optomechanical schemes building on higher-order couplings, it represents a very useful new resource.The method is also particularly straightforward since it simply relies on the strong pumping of the mechanics via the microwave field.Besides, the generated comb itself could be used in schemes requiring microwave multiplexing.One could imagine specific designs with multiple cavities and NEMS producing much wider combs.Further work is needed to understand the extent to which the geometrical nonlinearities affect the fluctuations of the system in the limit-cycle regime and to determine whether they can give rise to squeezing of the mechanical state.
( †) Corresponding Author: eddy.collin@neel.cnrs.frThe mechanical device used in this work is a typical aluminum drumhead [6].As can be seen on the SEM picture in Fig. 1, the actual structure is rather complex; we will simply approximate it as two discs of radius R (one being fixed and the other movable) separated by a gap d.The thickness of the drum is e.These geometrical characteristics are summarized in Tab.II together with typical material parameters.
These numbers are estimated from the Kirchhoff-Love theory of plates, producing the right mechanical resonance frequency of 6.7 MHz: assuming either highstress limit (in-built stress of 60 MHz and neglecting the Young's modulus) or low-stress (0 in-built stress).Besides, from Ref. [48] we can produce a theoretical estimate for the Duffing parameter ξ m in Hz/m 2 .We obtain about 2.×10 19 Hz/m 2 for a device in the high-stress limit (a drum), and about 1.×10 19 Hz/m 2 in the low-stress case (a membrane).From the fit value quoted in Tab.I in units of Hz/phonons, we get a number in between these two numerical estimates: this validates the quantitative evaluation within ±50 %.
The linear coupling strength g 0 is defined as: with x zpf = /(2m eff ω m ) the zero-point-fluctuation, defined from the mode effective mass m eff .For simplicity we will neglect the mode-shape here and consider two planar electrodes; as such, we will take as reference for the mode mass and spring constant calculation the center of the drum (i.e.maximum of mode shape equal to 1).Numbers are given in the caption of Tab.II.
In Eq. (A2), we have dω c /dC = −ω c /(2C 0 ) with C 0 the mode effective capacitance.From standard electromagnetism we write dC/dx = + 0 πR 2 /d 2 , neglecting fringing effects which are small in the limit d/R 1 ( 0 being the vaccum permittivity) [52].By definition, we take the direction of the X-axis pointing towards the fixed electrode.Reversing the direction of the X-axis changes the sign of g 0 but also of g 2 , producing an overall (−1) n in Eq. ( 11).This has no impact on physical quantities (such as γ BA , δω and |α n | 2 ): the problem at stake is invariant under a mirror symmetry.We then obtain from Eq. (A1) a value of about 20 Hz for g 0 (choosing g 0 > 0) taking for the cavity mode C 0 ≈ 100 fF, which is consistent with the microwave design.This over-estimates g 0 (by about a factor of two) since in reality not all the drum electrode moves, the borders being clamped.
Expanding the plate capacitor expression in a Taylor series of x/d, we obtain for the cavity resonance frequency: at third order, where we identify: In our case, only the first terms in the above are relevant: the magnitude with respect to g 0 of these g n coefficients is thus fixed by (x zpf /d) n .Computing numerical estimates, we see that with the chosen value of d we under-estimate g 1 by only about 20 %, but under-estimate |g 2 | by a factor of 7 approximately.The sign of g 2 is also not captured, which shows that this crude modeling fails for high-order derivatives.the theoretical dependence of this parameter on both P in and NEMS temperature T N EM S [19], we can recalculate T N EM S for each setting, see Fig. 6.This effect being local, the absorbed power has to be proportional to the intracavity field, i.e. the photon population n cav .We can therefore extrapolate what should be the heating effects in the self-oscillating regime using the actual intracavity photon number Σ|α n | 2 .Empirical fits are shown in Fig. 6 (see lines).The curves merge when the heating effect dominates over the starting temperature; we therefore estimate that our extrapolation in the region of interest is accurate within ±20 %, see Fig. 6.Of the parameters appearing in the theory of Section III, the only temperature dependent ones are ω c , ω m and γ m .The mechanical damping is fit to measurements performed in the Brownian regime by the expression γ m /(2π) = 70.5 + 1300 T N EM S , while the mechanical resonance frequency is fit by ω m /(2π) = 6.747×10 6 +430 ln(T N EM S ).
While for this sample, the mechanical element is very sensitive to heating, the microwave cavity seems to be rather insensitive.We attribute this to the fact that the cavity is much larger than the drum, and directly coupled to the substrate instead of being suspended.However, we do measure a power dependence of the microwave resonance frequency which shifts upwards logarithmically with increasing powers.At the same time, we do not measure any change in the cavity Q factor within our resolution.These power-dependencies of superconducting microwave resonators are commonly attributed to microscopic Two Level Systems present in the devices [49].Pragmatically, we take into account this effect by adding this logarithmic frequency shift to the calculation of ω c when fitting the 3D maps: δω c = 1.8 × 10 5 ln(P in ).
Similarly to the Duffing effect of the mechanics, there is an equivalent nonlinearity in the microwave resonance called Kerr nonlinearity.This leads to an additional frequency shift ∝ n cav .This effect comes from the nonlinear behavior of the mode effective inductance L 0 when the current density J flowing in the superconductors becomes too large [50,51]: with J * = (2/3) 3/2 J C and J C the critical current density, and α l = L kin /L 0 the fraction of the total inductance of kinetic origin.For our Al film of about 100 nm, α l should be smaller than 0.1 typically.The cavity resonance frequency thus shifts as: with: and A the cross-section of the microwave cavity strip.
A crude estimate taking the bulk value for the critical current density leads to ξ c ≈ −10 −4 Rad/s, which is completely negligible.
Appendix C: Impact of static deflection, Duffing and Kerr nonlinearities The ansatz Eq. ( 6) introduces a static term β c that corresponds to a static deflection of the drum x c = x zpf 2 [β c ].It can be deduced by solving Eq. ( 5) keeping only time-independent terms.The Duffing contribution can easily be incorporated in it.This term remains always extremely small, and contributes only for a (tiny) cavity frequency shift ω Since the mechanical motion can be very large (up to about 15 nm), the nonlinear stretching effect of the membrane has to be considered; this is the so-called Duffing nonlinearity, which shifts the mechanical resonance by ξ m B 2 .This term can be taken into account recursively in the calculation of the stable states, Eq. ( 20) see Appendix D. The result is that this term has only a marginal impact on the self-oscillating states definition.However, it dominates the mechanical frequency shift over the optical spring terms.We can therefore fit ξ m on the measurement of the drum frequency, Fig. 4. The obtained value is essentially temperature-independent, and given in Tab.I.
A similar nonlinear effect exists for the cavity: this is the Kerr effect already discussed in Appendix B. At first order, this term shifts the position of the resonance by a quantity ξ c |α| 2 .The expected value for ξ c being very small, we can simply completely neglect any nonlinear effect of that sort.

Appendix D: Stable states computation
The problem is solved numerically by finding selfconsistently a stable solution B to Eq. ( 20) for any couple (δ,P in ).These stability points correspond graphically to the intersection between the function γ BA (B)/γ m + 1 and the X-axis (see Fig. 7).The output photon flux is thus calculated by means of Eq. ( 22) injecting the found value of B in Eq. (11).For simplicity, one can neglect mechanical shifts which have only a marginal effect on stable states amplitudes.The procedure is then repeated over the full range of detunings δ and input pump powers P in in order to draw the theoretical mapping of the self-oscillating state (see green points in Fig. 3).For small detunings, the curves are always monotonous.
At low powers, there is no solution since γ BA (B)/γ m + 1 > 0 (orange line in Fig. 7).Increasing the power brings eventually the curve below the X-axis, creating a single intersection γ BA (B)/γ m + 1 = 0 (green circle on the red curve).Fluctuations at small B can thus trigger the self-oscillating state as the curve smoothly goes below Y = 0.
For large positive detunings, there is a range at (large) powers where the curve displays two intersections (see magenta line in Fig. 7).For the low-B valued one (blue circle), the slope is negative which means that the state displays anti-damping: it is unstable.On the other hand, for the high-B solution the derivative is positive, which means that the state is stable (dark green circle).
However, this state is at very large amplitudes B, and has not been created by a smooth crossing of the X-axis from the whole curve, starting at the lowest B ≈ 0: this means that it can be triggered only if one comes already from high amplitude states, and not from thermal motion.This is exactly the hysteretic behavior that is seen in Fig. 3, sweeping the detuning δ upwards at constant power, and increasing the power from typically 2 nW to 30 nW.The same is true, sweeping the power downwards from the high-B state at fixed detuning.
The graphs in Fig. 7 are obtained with g 1 = +10 −7 g 0 and g 2 = −10 −13 g 0 (with g 0 > 0).The numerical calculation can be performed with the static deflection and the Duffing term taken into account, see Appendix C. The results are essentially identical.The fitting routine is explained in the next Appendix.
Appendix E: 3D Fitting procedure Measurements of the photon flux are compared to the theoretical computation in Fig. 3.The two parameters g 1 and g 2 affect the shape of the numerical (δ, P in ) colormap in different ways: we demonstrate this in Fig. 8 varying them in a dichotomic process (multiplying or dividing the optimal values by 2).By increasing g 1 the self-oscillating region is getting more narrow in the δ direction, while increasing |g 2 | up-shifts in power the starting line of the bistable region.This is discussed in the core of the paper with the parameters ∆P , ∆δ, see Fig. 5.The optimal values match the experimental findings: ∆P ≈ 16 nW (±10 %), ∆δ ≈ 3.5 MHz (±200 kHz).
We can therefore choose the red cross position in Fig. 8 (central graph, optimal g 1 and g 2 ) as a good marker for fitting these g 1 and g 2 parameters (equivalent of Fig. 2, but g 1 = g 2 = 0 value).Since the NEMS heats with applied power, this also defines the actual temperature at which the fit is essentially performed.This is summarized in Tab.I, with error bars estimated for the coupling nonlinear parameters to be about a factor of 2. Fits of the mechanical frequency shifts are discussed in Appendix C; Fig. 4 is essentially an image of the amplitude of motion squared x 2 (or equivalently B 2 ).Note that the quality of the agreement between experiment and theory in this graph also imposes strong constraints on the (g 1 , g 2 ) couple.This is also the case of the overall shape of the photon flux maps, Fig. 3.

FIG. 1 :
FIG.1:(a) Main: Power spectral density (PSD) of the Stokes peak (i.e. at the frequency of the cavity) measured in the selfsustained oscillation regime at 214 mK (blue-detuned pump power Pin of 6 nW with ∆ = +ωm).Left inset: Time domain measurement of the coherent signal (raw data units).Right inset: SEM picture of a drumhead type resonator.(b): PSD measurement of the comb produced by the strong applied power in same conditions (pump signal 6 nW, green arrow at frequency (ωc + ωm)/2π).The cavity (orange area) is displayed with an arbitrary amplitude and its linewidth κtot/2π at scale.The black points are theoretical computation of the output amplitude of each measured peak (see text).

FIG. 2 :
FIG. 2: (a): Measured output photon flux (Stokes peak) as a function of input power Pin and detuning δ at 214 mK.(b):Corresponding calculated colormap from the basic theory described in Refs.[32,33] (no nonlinearities, g1 = g2 = 0 in Section III).The region on the right of the dashed line (high powers, positive detuning) is bistable and exists only when entering from the self-oscillating state (up-sweeping frequencies).The pink cross marks the minimum power necessary for self-oscillations, while the red cross corresponds to the position of the beginning of the hysteresis.The yellow cross marks the end of this bistable region (at same power).∆P and ∆δ are discussed in the text.

6 FIG. 3 :
FIG. 3: (a): Output photon flux of the self-oscillating (Stokes) peak at frequency ωc + δ, as a function of both the power Pin and the detuning δ of the input pump signal (pump frequency ωc + ωm + δ, with -7 MHz< δ <+7 MHz) at 214 mK.The colormap is experimental data measured upsweeping both the pump detuning (from δ = −7 MHz to δ = +7 MHz) and the pump power, and green points are theoretical fits computed by solving self-consistently Eq. (20), γm + γBA = 0, see text.(b): Experimental colormap measured down-sweeping the pump detuning (from δ = +7 MHz to δ = −7 MHz) with pump power swept upwards.Green points are also theoretical computations; the hysteresis of the large power and large detuning region is clearly visible.

FIG. 4 :
FIG. 4: (a): Mechanical frequency shift of the self-oscillating (Stokes) peak as a function of both the power Pin and the detuning δ of the input pump signal (same conditions as Fig. 3, sweeping δ towards positive values).Note that in the hysteretic region, the calculated points lie slightly below the experimental ones, but obviously match the threshold position of Fig. 3. (b): Experimental colormap measured down sweeping the pump detuning.Green points are theoretical computations, see text.
FIG. 5: (a): Calculated ∆P parameter as a function of g1, g2 coefficients.(b): Calculated ∆δ parameter as a function of g1, g2.Both are essentially described by plane equations, with each nonlinear coefficient being the leading one for one of the parameters (g1 for ∆δ and g2 for ∆P , see text).Full colormaps are also presented in a matrix form in Appendix E, Fig. 8.
FIG.6: (Color online) Mechanical device temperature versus applied microwave power (expressed in terms of intra-cavity photons ncav).Dots are experimental data measured by reddetuned pumping, integrating the anti-Stokes power spectrum peak (see text).The curves are empirical expressions used for the extrapolation in the self-oscillating range (above 10 8 photons).At high enough powers, all the curves collapse, as they should.The discrepancies in the numerics in the extrapolated range is smaller than ±20 %.

FIG. 7 :
FIG.7: Theoretical stability curves giving γBA/γm + 1 as a function of (2g0/ωm) × B, calculated at 3 different positions (δ, Pin) and demonstrating the typical observed behaviors: unstable (orange line), one stable state (red line), and unstable (blue circle) plus stable states leading to hysteresis (magenta line, see text and Figs.3 and 4).The self-consistent value of B corresponds to the (light and dark) green circles.

FIG. 8 :
FIG.8: Impact of the variation of g1 and g2 on the theoretical map giving the output photon flux as a function of both the detuning δ and the input pump power Pin.The colormaps are calculated taking into account all mechanical and optical shifts, with g0 > 0. The central one is the same as in the 3D plot of Fig.3.From these graphs, one can extract the ∆P , ∆δ parameters shown in Fig.5.The red cross marks the position of the beginning of the hysteresis for the central graph (optimal g1, g2 fit parameters).

TABLE I :
Fitted parameters at different temperatures."cryo." is the cryostat measured temperature, while T eff is the characteristic fit temperature; nonlinear couplings are given in units of g0, with g0 > 0. The Kerr parameter of the cavity is estimated (see text).

TABLE II :
Typical drumhead NEMS parameters; the in-built stress is estimated to be < 60 MPa (see text).Corresponding mode effective mass m eff = 2.3×10 −14 kg and spring constant k eff = 41.N/m.
Appendix A: Drumhead characteristics