Graphene multimode parametric oscillators
Abstract
In the field of nanomechanics, parametric excitations are of interest since they can greatly enhance sensing capabilities and eliminate crosstalk. However, parametric excitations often rely on externally tuned springs, which limits their application to high quality factor resonators and usually does not allow excitation of multiple higher modes into parametric resonance. Here we demonstrate parametric amplification and resonance of suspended singlelayer graphene membranes by an efficient optothermal drive that modulates the intrinsic spring constant. With a large amplitude of the optical drive, a record number of 14 mechanical modes can be brought into parametric resonance by modulating a single parameter: the pretension. In contrast to conventional mechanical resonators, it is shown that graphene membranes demonstrate an interesting combination of both strong nonlinear stiffness and nonlinear damping.
\ul
The history of parametric oscillations dates back to the 19th century and the observation of surface waves in the famous singing wineglass experiment of Michael Faraday Faraday (1831). The advent of micro and nano engineering brought to life new ideas for exploiting parametric excitation for enhancing force and mass sensitivity Turner et al. (1998); Rugar and Grütter (1991); Karabalin et al. (2009, 2011); Zhang and Turner (2005, 2004); Zhang et al. (2002), effective quality factor Mahboob and Yamaguchi (2008a), and signal to noise ratio Rugar and Grütter (1991) of tiny resonators. To date, many sensors, including gyroscopes OropezaRamos and Turner (2005); Hu et al. (2011); Harish et al. (2008), mass sensors Zhang and Turner (2005, 2004); Zhang et al. (2002) and even mechanical memories Mahboob and Yamaguchi (2008b); Mahboob et al. (2014); Roukes (2004); Freeman and Hiebert (2008) employ parametric excitation for improved performance.
Parametric micromechanical oscillators normally rely on stiffness modulation. The mechanical stiffness of microbeams is however fully determined by their fixed Young’s modulus and geometry; parametric excitation can therefore usually only be achieved by modulating the stiffness of an externally applied spring, for example a voltage controlled electrostatic spring Turner et al. (1998); Rugar and Grütter (1991); Carr et al. (2000); Rhoads et al. (2010); Mathew et al. (2016); Prasad et al. (2017). However, this method has a small modulation amplitude, limiting its application to resonators with a high quality factor. Besides, the external spring force usually cannot excite all the higher mechanical modes of a single mechanical element.
Here, we overcome these limitations by using tension modulation for parametric driving of suspended graphene membranes. Since graphene membranes have negligible bending rigidity, their stiffness is dominated by their pretension that can be efficiently modulated by heat Ye et al. (2017); Yang et al. (2017); Davidovikj et al. (2017a). The tension modulation is given by , where is the thermal expansion coefficient, the 2D Young’s modulus and the temperature modulation. Using approximate values from literature Lee et al. (2008); Yoon et al. (2011), one finds that NmK, which means that a temperature modulation of 1 K already results in a tension modulation of the order of the intrinsic pretension (estimated to be between 0.003 N/m and 0.03 N/m Dolleman et al. ()) of the graphene membranes studied here. One can define the relative shift of the resonance frequency per unit of temperature as a figure of merit for the efficiency of the optothermal parametric drive: to 1 K. This estimated value for graphene is 5005000 times larger than in other optically pumped oscillators Zalalutdinov et al. (2001), thus illustrating that the parametric driving scheme for graphene membranes is possibly the most efficient method for reaching parametric oscillation in mechanical systems. Another advantage of tension modulation in membranes is that all mechanical degrees of freedom can be parametrically excited, since the effective stiffness of every mechanical mode of the system is proportional to .
We demonstrate the optothermal tension modulation technique experimentally on single layer graphene membranes and show that as many as 14 mechanical modes can be parametrically excited. A detailed analysis of both the directly and parametrically driven response of the fundamental mode reveals that nonlinear damping is essential to describe the mechanical motion. We experimentally demonstrate parametric signal amplification in these resonators, raising the amplitude of motion and increasing the effective quality factor of resonance. Finally, by analyzing the parametric resonances, we obtain information on the dissipation mechanism and loss tangent in graphene membranes, a subject that has recently received much attention Güttinger et al. (2017); Eichler et al. (2011); Singh et al. (2016); Croy et al. (2012).
Optothermal parametric driving of graphene
Experiments are performed on singlelayer CVD graphene drum resonators with a diameter of 5 m and a cavity depth of 300 nm. The drums have venting channels to the environment to prevent the trapping of gas in the cavity (Fig. 1a,b, see Methods section for details on the fabrication). To achieve parametric drive, we use the experimental setup shown in Fig. 1c. The light from a blue diode laser is focused on the membrane and its intensity is modulated by an input voltage . This periodically heats up the membrane and creates a parametric drive due to the thermal strain. Parametric resonance occurs if the parametric driving term exceeds a threshold , determined by the resonance frequency and quality factor of resonance Rugar and Grütter (1991). At the same time, imperfections such as initial outofplane deformations, wrinkles and ripples in the membrane geometry enable this laser to directly drive the resonator by thermal expansion force, because thermal expansion will enhance these deformations and thus actuate the membrane. A more detailed discussion on this mechanism can be found in the Supplementary Information S3.
A red heliumneon laser is used to read out the motion by the optical interference between the graphene membrane and the fixed substrate Bunch et al. (2007); CastellanosGomez et al. (2013). The ratio of the AC voltage amplitude generated by the photodetector and the AC driving voltage of the blue laser / is determined by a vector network analyzer (VNA). The VNA has the ability to perform frequency conversion measurements, hence both homodyne and heterodyne detection can be performed in this setup such that both direct and parametric resonances can be analyzed.
Graphene multimode parametric oscillators
In Fig. 2a, the blue laser is driven at 2, while detecting the photodiode signal at . When increasing the blue laser driving voltage a remarkable effect is observed. Onebyone, the parametric resonances of graphene appear, up to 7 different modes. Each mode reaches resonance at a different threshold driving amplitude , due to differences in quality factor and the frequency dependence of the parametric driving parameter Dolleman et al. (). The experiment is repeated on a different drum in Fig. 2b. Interestingly, in this case overlap between parametric resonances is observed at high driving levels. When overlap occurs, a direct transition between the highamplitude solution of two adjacent parametric resonances is observed, e.g. at mV (RMS) between the second and third resonance. Interestingly, in some cases also transitions between the highamplitude and lowamplitude solutions are observed, e.g. at mV (RMS) between the same 2 modes. This pseudorandom process is attributed to a strong dependence of the basin of attractions of the parametric highamplitude and lowamplitude solutions on the initial conditions Sanchez and Nayfeh (1990). Hence, the amplitude can fall into two stable solutions: either the high amplitude solution of the third mode or the zero amplitude solution of the third mode which is also observed at higher driving amplitudes ( and mV (RMS)).
Due to the overlap of parametric resonances in this drum, some resonances are skipped and not all resonances are found by sweeping from low to high frequency. Instead, when sweeping the frequency backward as shown in Fig. 2c, as many as 14 parametric resonances are observed in this system. To our knowledge, this is the largest number of parametrically excited modes in a single mechanical element, as previously only 7 modes could be excited in cryogenic environments Mahboob et al. (2014).
Mechanical nonlinearities in graphene resonators
For a more detailed analysis of the physics, we focus on the frequency response of the fundamental mode to both direct and parametric drives. Figure 3 shows direct and parametric resonance of the fundamental mode as function of driving level, on a different drum than Fig. 2. The VNA is configured to detect the directly driven frequency response (Fig. 3a,c). Sweeping the frequency forward (Fig. 3a) and backward (Fig. 3b) results in a hysteresis, that grows as the driving level is increased. This is typical for the geometric nonlinearity of the Duffingtype resonator, where the stiffness becomes larger at high amplitudes. In order to detect the parametric resonance, the VNA was configured in a heterodyne scheme at which is detected at half of the driving frequency . Similar to the directly driven case, a hysteresis occurs between the forward (Fig. 3b) and backward (Fig. 3d) sweeps in frequency. Below an RMS drive amplitude of 0.11 mV, no response is observed. To show that the parametric resonance shows two stable phases of resonance separated by 180 degrees Mahboob and Yamaguchi (2008b), we performed an additional measurement which is shown in the Supplementary Information S2.
Figure 4ad shows both directly and parametrically driven responses at different driving levels. In order to eludicate the effect of nonlinearities on the observed mechanical responses, a single degreeoffreedom model is derived that describes the motion of the resonator (see Supplementary Information S4) and this is fitted to the response curves in Fig. 4ad (see Supplementary Information S5). The model is a combination of the Duffing, van der Pol and MatthieuHill equation also used in other works Eichler et al. (2011); Houri et al. (2017); Lifshitz and Cross (2008); Aubin et al. (2004):
(1) 
where is the displacement (which is approximately proportional to ), is the damping coefficient, the nonlinear damping coefficient, the linear stiffness coefficient, the nonlinear stiffness coefficient, the parametric driving and the direct driving term. By setting and one can fit the response at low drive level (Fig. 4a) and obtain an initial value for , and . Initially fitting at high driving levels was attempted by setting , however it is found that such a model cannot account for the observed parametric response: nonlinear damping is indispensable to describe the maximum amplitude. Next, , and are used as fitting parameters to describe the nonlinear response (Fig. 4ad). Numerical values for the fit parameters are provided in the Supplementary information S1.
Figure 4 compares the fitted model and the experimental data for the directly and parametrically driven fundamental resonance. This shows excellent agreement at lower driving levels (Fig. 4a). We note that the fitting parameters , , and are the same in the direct and parametric response within the error of the fitting procedure (see Supplementary information S1) and that both and are very nearly proportional to driving voltage . It is however observed that the region of instability (Figs. 4e, f) is narrower in our experiments than what is expected from eq. 1.
Parametric signal amplification in graphene
Now we investigate the effects of parametric drive at low driving levels () by examining parametric amplification of the directly driven resonance. To measure parametric amplification, it is required to simultaneously drive the system at and (where is near the resonance frequency ). This is realized by splitting the driving circuit connected to the diode laser into two parts (Fig. 1c). One path provides a small direct drive that excites the primary resonance of the membrane in the linear regime. The second path contains a frequency doubler, amplifier and phase shifter to enable parametric driving with controllable phase and gain with respect to the direct drive. A harmonic oscillator model is fitted to the response to extract the amplitude and the effective quality factor. The relation between amplitude gain , parametric drive amplitude and phase shift of the direct drive is given by Mahboob and Yamaguchi (2008a); Rugar and Grütter (1991):
(2) 
First, the amplification effect as function of parametric pumping amplitude in Fig. 5a was examined by keeping the phase fixed at 45 degrees. Increasing the amplitude of parametric drive increases the amplitude at resonance by a factor of 34 (Fig. 5b) and the effective quality factor of resonance by almost a factor of 3 (Fig. 5c). Figure 5d shows that shifting the phase of the parametric drive significantly changes the amplitude of harmonic resonance. Figure 5ef shows that the gain and effective Qfactor depend strongly on the phase of the parametric drive with respect to the direct drive. Fits of the data in Fig. 5b, e show that the drive and phasedependence of the parametric amplification is in accordance with theory.
Discussion
We have shown optothermal tension modulation as a mechanism to achieve large parametric excitations in graphene. The large parametric excitation enables low modes to be brought into parametric resonance. This allows parametric resonance study at room temperature (where graphene’s factor is a factor 100 lower than in the cryogenic regime Chen et al. (2013); Zande et al. (2010)), and also allows higher order modes (with lower Q) to be brought into parametric resonance. These beneficial conditions enable us to realize the graphene multimode parametric oscillators as shown in Fig. 2. Due to their small mass and low spring constant, graphene membranes are very force sensitive. Here we demonstrate parametric amplification of graphene, which can further push graphene’s limits of force sensitivity, since it enhances both the gain and effective quality factor of the resonators.
The asymmetry around observed in the region of instability (Fig 4) is a surprising result: such an asymmetry should not arise for the equation of motion (Eq. 1) used in the analysis. Something similar is observed in the directly driven response, where the lower saddle node bifurcation in the downward frequency sweep is always found at a lower frequency than simulated (Fig. 4c) at high driving levels. Possibly, this indicates that the forcing terms are nonlinear Rhoads et al. (2006). However, we find that both forcing terms and extracted from the fits are linear with the applied modulation amplitude and the forward frequency sweeps are welldescribed by this model (see Supplementary Information S1). The observed deviations (e.g. in Fig. 4) can therefore not be explained by forcing nonlinearities.
The asymmetry and apparent decrease in resonance linewidth (Fig. 4) thus suggest that a more unconventional dissipation model should be considered, including further terms to describe the amplitudedependence of the dissipation. Similar deviations from conventional dissipation models have been previously found in multilayered graphene resonators Singh et al. (2016), where it was concluded that the vanderPol term does not describe the nonlinear damping. Here we conclude that the vanderPol term is generally in agreement with the experiments, since it describes the saddle node bifurcation of the parametric resonances well, however additional dissipation terms might be needed to account for the asymmetry and narrowing of the parametric stability region (Fig. 4e,f).
The fit to the nonlinear response of the membrane allows us to extract a number for the Duffing () and vanderPol terms () in our resonators. As shown in the Supplementary information S6, the mechanical loss tangent of graphene at the resonance frequency can be determined from the ratio of these terms, . From the values of the fits we obtain for drum 2 and for drum 3. The values of these loss tangents are in the same range as found by Jinkins et al.Jinkins et al. (2015). The obtained values for the loss tangent are relatively high for a crystalline material as graphene, therefore the observed nonlinear damping is likely not due to the intrinsic material properties but to other effects, such as sidewall adhesion Bunch and Dunn (2012) or unzipping of wrinkles RuizVargas et al. (2011).
Multimode parametric oscillators are interesting for applications where accurate frequency tracking of multiple modes is necessary. All of these modes can be parametrically amplified up to relatively high amplitudes. This is difficult to achieve with conventional oscillators that require a feedback loop and special filters or actuation schemes to make sure only the desired resonance is brought into oscillation. A possible application is inertial imaging Hanay et al. (2015), where accurate tracking of multiple resonances allows one to determine the mass, location and shape of a particle on top of a resonator. In addition, these parametric oscillators can be used to build a binary information and computation system Mahboob et al. (2014), where information is stored in the phase of the resonator. Multimode resonators have the potential of enabling parallel processing and data storage. The high resonance frequencies and relatively low Q of the graphene membranes can increase computation speed. Moreover, since the driving frequency is double the readout frequency, parametric driving is less sensitive to crosstalk that is often present in directly driven resonators Turner et al. (1998) which can improve sensor performance.
Conclusions
In conclusion, we report on multimode parametric resonance and amplification in single layer graphene resonators by an optothermal tension modulation technique. It is demonstrated that the tensiondominated restoring force results in parametric excitation of multiple resonance modes in the system when the system is optothermally driven. The parametrically and directly driven resonances are compared to a single degreeoffreedom model based on the Duffing, van der Pol and Matthieu equations, with good agreement at low driving levels. This allows simultaneous determination of nonlinear stiffness and damping coefficients and results in a highfrequency determination of graphene’s mechanical loss tangent. It was demonstrated that weak parametric drives can be used to amplify the motion and enhance the effective quality factor of resonance, that can potentially enhance the force sensitivity of graphene resonators. Graphene resonators are thus an interesting platform to study parametric excitations and their utilization for sensors with improved performance.
Acknowledgements.
We gratefully acknowledge Applied Nanolayers B.V. for the growth and transfer of the single layer graphene used in this study. We thank Y.M. Blanter, B.J. Gallacher, W.J. Venstra and S.W. Shaw for useful discussions and D. Davidovikj for help with scanning electron microscopy. This work is part of the research programme Integrated Graphene Pressure Sensors (IGPS) with project number 13307 which is financed by the Netherlands Organisation for Scientific Research (NWO). The research leading to these results also received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 649953 Graphene Flagship.Methods
Fabrication
Graphene resonators are fabricated by etching dumbbellshaped cavities in a thermally grown, 285 nm SiO layer on a silicon wafer. The etching did not fully stop at the silicon layer, resulting in cavities that are 300 nm thick. Circular membranes are formed by transfer of single layer chemical vapor deposited (CVD) graphene (Fig. 1a). During the transfer process one side of the dumbbell is broken while the other side remains intact, creating a circular resonator on one side with a venting channel to the environment (Fig. 1b). This prevents gas from being trapped in the cavity when the pressure in the surroundings changes. In the main section of this work four identical drums with a diameter of 5 micrometer are used; results obtained on drum 1 are shown in Figs. 2a, 3 and 4f, drum 2 in Figs. 2b,c, drum 3 in Figs. 4ae and drum 4 in Fig. 5. A fifth drum was used in the Supplementary Information S2. More details on the fabrication and transfer process of the drum resonators can be found in ref. Dolleman et al. ().
Experiments
All measurements are performed at room temperature in a high vacuum environment with a pressure less than mbar to minimize the effects of gas damping. The blue diode laser (Thorlabs LP405SP10) has a wavelength of 405 nm and is biased with a 32 mA current, resulting in 0.76 mW of incident power measured before the objective. The red laser illuminates the sample with 1.2 mW of incident power (measured before the objective lens). The vector network analyzer is of type Rohde & Schwarz ZNB4 with the frequency conversion option (k4) installed.
Author contributions
R.J.D., S.H., F.A., A.C. and P.G.S. analyzed and interpreted the measurements. R.J.D. and S.H. performed the experiments. F.A., A.C. and P.G.S. derived the equation of motion and fits in Fig. 4. H.S.J.v.d.Z. and P.G.S. supervised the project. The manuscript was jointly written by all authors with a main contribution from R.J.D.
Supplementary Information: Graphene multimode parametric oscillators
In section S1, we show the complete dataset obtained while analyzing the nonlinearities of the graphene membrane. Section S2 shows an additional experiment that demonstrates the parametric oscillator has two stable phases. In section S3, an additional discussion is added that proposes why both direct and parametric excitations are observed in the same setup. Section S4 derives the equations of motion that has been used to perform the fitting and section S5 describes the numerical simulations used for the fitting procedure. Finally in section S6 the expression for the mechanical loss tangent of graphene is derived.
S1: Complete datasets for the analysis of mechanical nonlinearities
Figure 4 in the main text shows the fitting of the nonlinear mechanical response of the resonator (drum 3). In this section the remainder of this analysis is presented and the complete dataset from the fundamental mode of drum 1 is shown (Figs. 2a, 3, 4f in the main text).
Direct drive  Parametric drive  

RMS Drive (mV) 



250.9  0.0045  70  225  8  0.0045  76  215  1.06  
281.5  0.0045  72  220  9.2  0.0045  76  220  1.22  
354.4  0.0045  74  230  11.7  0.0045  79  225  1.6  
397.6  0.0045  76  230  14.2  0.0046  80  225  1.8  
446.2  0.0045  76  230  15.5  0.0046  80  225  1.93 
.1 Dataset of drum 1
Direct drive  Parametric drive  

RMS Drive (mV) 



149.4  0.0030  36  250  1.42  0.0030  36  250  0.74  
167.7  0.0030  37  245  1.6  0.0030  36  220  1.01  
188.1  0.0030  37  245  2.0  0.0030  34  225  1.18  
211.1  0.0030  37  245  2.5  0.0030  34  225  1.31  
236.9  0.0030  37  250  2.8  0.0030  33  225  1.46  
265.8  0.0030  36  250  3.3  0.0030  34  225  1.81  
298.2  0.0030  35  250  3.9  0.0030  35  225  2.05  
334.6  0.0030  35  250  4.5  0.0030  35  225  2.25 
S2: Additional experiment showing two stable phases
The frequency conversion option on the vector network analyzer loses information on the phase at which the resonator is oscillating. To show that the parametrically excited resonance has two stable phases separated by 180 degrees, the experiment was repeated by using the frequency doubler in the circuit used for the parametric amplification experiment (Fig. 1(c) in the main text). Using this, the VNA does not require to perform a frequency conversion and phase information is preserved. This results in the mechanical responses shown in Fig. S4.
S3: Additional discussion: mechanism for direct and parametric driving
Optothermal driving leads to two mechanisms that can excite the resonance in the graphene resonators. Parametric drive (Fig. S5a) occurs due to the modulation of pretension in the membrane via laser heating and thermal expansion, since the stiffness term for the outofplane deflection field of the membrane determined by the pretension. Parametric driving will only activate the parametric resonance if the modulation of the blue laser is near twice the mechanical resonance frequency.
As demonstrated in the main text (Figs. 3, 4), the experiments also show a direct driving component. This can be explained Aubin et al. (2004) by assuming a small initial membrane displacement from equilibrium (Fig. S5b). In graphene resonators rippling, wall adhesion or outofplane crumples could lie at the root of such an initial displacement.
In order to analyze the data, we will derive the equations of motion (Eq. S10) using a Lagrangian approach by including this initial deflection field. In this manner, the equations are reduced to a singledegreeoffreedom model that can be used to fit to the data, significantly simplifying the analysis. The derivation of this single degree of freedom model is shown below in section S4.
S4: Equations of motion
A Lagrangian approach is used for obtaining equations of motion of an optothermally excited monolayer graphene membrane. In this respect, the potential energy of the thermally actuated circular membrane is obtained as Amabili (2008):
(S1) 
where is the thickness, is the radius, is the thermal expansion coefficient, and is the temperature change in the membrane. Moreover, , , , are the Kirchhoff stresses that can be obtained as follows:
(S2)  
in which , , and are the Green strains and are derived as:
(S3)  
where , and are the radial, tangential and transverse displacements, respectively. Moreover, is the deviation of the membrane from flat configuration, is the Young’s modulus and is the Poisson’s ratio.
The temperature difference can be obtained by solving the following heat conduction equation:
(S4) 
in which is the power absorbed by the membrane, is the thermal time constant Dolleman et al. (), is the thermal capacitance, and represents the time variable.
For a membrane with fixed edges and shall vanish at . Moreover, should be zero at for continuity and symmetry. Furthermore, assuming only axisymmetric vibrations ( and ), the solution can be approximated as Davidovikj et al. (2017b):
(S5) 
(S6) 
Here it should be noted that for axisymmetric vibrations the shear strain would become zero. In equation (S5), is the generalized coordinate associated with the fundamental mode of vibration. Furthermore, in equation (S6), ’s are the generalized coordinates associated with the radial motion. Moreover, is the zeroth order Bessel function of the first kind and . In addition, is the number of necessary terms in the expansion of radial displacement and is the initial displacement due to pretension that is obtained from the initial stress as follows :
(S7) 
The kinetic energy of the membrane neglecting inplane inertia, is given by:
(S8) 
The Lagrange equations of motion are given by:
(S9) 
and q=[,( )], is the vector containing all the generalized coordinates. Equation (S9) leads to a system of nonlinear equations comprising of a single differential equation associated with the generalized coordinate and algebraic equations in terms of () . By solving the algebraic equations it is possible to determine in terms of Davidovikj et al. (2017b). This will reduce the +1 set of nonlinear equations to the following DuffingMatthieuHill equation:
(S10) 
where represents derivative with respect to time and is the mass. and are the linear viscous damping coefficient and nonlinear material damping coefficient, respectively Zaitsev et al. (2012); Lifshitz and Cross (2008). They are added to the equation of motion explicitly to introduce dissipation. represents the linear stiffness term dominated by the pretension and is the amplitude of parametric drive resulting from temperature variation . Moreover, represents the quadratic nonlinear stiffness coefficient due to imperfection and denotes the cubic nonlinear stiffness coefficient arising from geometric nonlinearity. Finally, is the amplitude of direct drive term due to the presence of imperfection , and is the excitation frequency. Indeed for a flat membrane, .
S5: Numerical simulations
In order to perform the numerical simulations, equation (S10) is normalized with respect to the mass of the membrane and the fundamental frequency () as follows:
(S11) 
Introducing an effective stiffness nonlinearity , whose value is given by =) Nayfeh and Mook (1995), equation (S11) is reduced to:
(S12) 
where the normalized coefficients are given in table 3.
Normalized parameter  

Scaled time derivative  
Nondimensional excitation frequency  
Scaled linear damping coefficient  
Scaled nonlinear damping coefficient  
Scaled linear stiffness coefficient  
Scaled parametric excitation amplitude  
Scaled nonlinear quadratic stiffness coefficient  
Scaled nonlinear cubic stiffness coefficient  
Scaled effective nonlinear stiffness coefficient  
Scaled direct excitation amplitude 
Here it should be noted that, mass of the single layer graphene membrane is unknown. Without the exact mass value, optical transduction factors present between the voltage signal measured by the VNA during the experiment and the actual motion of the membrane in physical units cannot be calibrated. Thus, the normalized coefficients shown in table 3 include a linear transduction factor ’’ for the oscillation amplitude (), for the parametric drive amplitude () and for the direct drive amplitude (). Where and are voltage signals measured in the experiment.
Finally, the equation (S12) is simulated using a pseudo arc length continuation and collocation technique Doedel et al. (1998) to detect bifurcations and obtain periodic solutions. The simulations are performed as follows:

The bifurcation analysis is carried out with the coefficient as the first continuation parameter and is incremented to the desired value in order to match the experimental direct response.

Once the desired value of is obtained, the parametric drive amplitude is used as the second continuation parameter and a value is chosen to replicate the experimental parametric response.

After reaching the desired value, the analysis is continued with the frequency ratio as the final continuation parameter. This value is spanned around the spectral neighborhood of and in order to obtain the direct and parametric response curves.
S6: Mechanical loss tangent of graphene
In ref. Davidovikj et al. (2017b) it is shown that the Duffing term is proportional to the Young’s modulus :
(S13) 
where is a constant. In case of material damping, a complex Young’s modulus can be introduced: and the nonlinear stifness term near the resonance frequency , for becomes:
(S14) 
From this equation it can be seen that the loss tangent Lakes (2009) can be calculated by the ratio if the resonator is vibrating near its resonance frequency:
(S15) 
References
 Faraday (1831) M Faraday, ‘‘On a peculiar class of acoustical figures; and on certain forms assumed by a group of particles upon vibrating elastic surfaces,” Philisophical Transactions of the Royal Society (London) 121, 299–318 (1831).
 Turner et al. (1998) Kimberly L Turner, Scott A Miller, Peter G Hartwell, Noel C MacDonald, Steven H Strogatz, and Scott G Adams, “Five parametric resonances in a microelectromechanical system,” Nature 396, 149–152 (1998).
 Rugar and Grütter (1991) D Rugar and P Grütter, ‘‘Mechanical parametric amplification and thermomechanical noise squeezing,” Physical Review Letters 67, 699 (1991).
 Karabalin et al. (2009) RB Karabalin, XL Feng, and ML Roukes, “Parametric nanomechanical amplification at very high frequency,” Nano Letters 9, 3116–3123 (2009).
 Karabalin et al. (2011) RB Karabalin, Ron Lifshitz, MC Cross, MH Matheny, SC Masmanidis, and ML Roukes, “Signal amplification by sensitive control of bifurcation topology,” Physical Review Letters 106, 094102 (2011).
 Zhang and Turner (2005) Wenhua Zhang and Kimberly L Turner, “Application of parametric resonance amplification in a singlecrystal silicon microoscillator based mass sensor,” Sensors and Actuators A: Physical 122, 23–30 (2005).
 Zhang and Turner (2004) Wenhua Zhang and Kimberly L Turner, “A mass sensor based on parametric resonance,” in Proceedings of the Solid State Sensor, Actuator and Microsystem Workshop, Hilton Head Island, SC (2004) pp. 49–52.
 Zhang et al. (2002) Wenhua Zhang, Rajashree Baskaran, and Kimberly L Turner, “Effect of cubic nonlinearity on autoparametrically amplified resonant MEMS mass sensor,” Sensors and Actuators A: Physical 102, 139–150 (2002).
 Mahboob and Yamaguchi (2008a) I Mahboob and H Yamaguchi, “Piezoelectrically pumped parametric amplification and Q enhancement in an electromechanical oscillator,” Applied Physics Letters 92, 173109 (2008a).
 OropezaRamos and Turner (2005) Laura A OropezaRamos and Kimberly L Turner, “Parametric resonance amplification in a MEMGyroscope,” in Sensors, 2005 IEEE (IEEE, 2005) pp. 4–pp.
 Hu et al. (2011) ZX Hu, BJ Gallacher, JS Burdess, CP Fell, and K Townsend, “A parametrically amplified MEMS rate gyroscope,” Sensors and Actuators A: Physical 167, 249–260 (2011).
 Harish et al. (2008) KM Harish, BJ Gallacher, JS Burdess, and JA Neasham, “Experimental investigation of parametric and externally forced motion in resonant MEMS sensors,” Journal of Micromechanics and Microengineering 19, 015021 (2008).
 Mahboob and Yamaguchi (2008b) I Mahboob and H Yamaguchi, “Bit storage and bit flip operations in an electromechanical oscillator,” Nature Nanotechnology 3, 275–279 (2008b).
 Mahboob et al. (2014) I Mahboob, M Mounaix, K Nishiguchi, A Fujiwara, and H Yamaguchi, “A multimode electromechanical parametric resonator array,” Scientific Reports 4 (2014).
 Roukes (2004) ML Roukes, “Mechanical compution, redux?,” in Electron Devices Meeting, 2004. IEDM Technical Digest. IEEE International (IEEE, 2004) pp. 539–542.
 Freeman and Hiebert (2008) Mark Freeman and Wayne Hiebert, “NEMS: Taking another swing at computing,” Nature Nanotechnology 3, 251–252 (2008).
 Carr et al. (2000) Dustin W Carr, Stephane Evoy, Lidija Sekaric, HG Craighead, and JM Parpia, “Parametric amplification in a torsional microresonator,” Applied Physics Letters 77, 1545–1547 (2000).
 Rhoads et al. (2010) Jeffrey F Rhoads, Steven W Shaw, and Kimberly L Turner, “Nonlinear dynamics and its applications in microand nanoresonators,” Journal of Dynamic Systems, Measurement, and Control 132, 034001 (2010).
 Mathew et al. (2016) John P Mathew, Raj N Patel, Abhinandan Borah, R Vijay, and Mandar M Deshmukh, “Dynamical strong coupling and parametric amplification of mechanical modes of graphene drums,” Nature Nanotechnology (2016).
 Prasad et al. (2017) Parmeshwar Prasad, Nishta Arora, and Akshay Naik, “Parametric amplification in MoS2 drum resonator,” Nanoscale , – (2017).
 Ye et al. (2017) Fan Ye, Jaesung Lee, and Philip XL Feng, “Verywide electrothermal tuning of graphene nanoelectromechanical resonators,” in Micro Electro Mechanical Systems (MEMS), 2017 IEEE 30th International Conference on (IEEE, 2017) pp. 68–71.
 Yang et al. (2017) Rui Yang, Changyao Chen, Jaesung Lee, David A Czaplewski, and Philip XL Feng, “Localgate electrical actuation, detection, and tuning of atomiclayer MoS2 nanoelectromechanical resonators,” in Micro Electro Mechanical Systems (MEMS), 2017 IEEE 30th International Conference on (IEEE, 2017) pp. 163–166.
 Davidovikj et al. (2017a) Dejan Davidovikj, Herre SJ van der Zant, and Peter G Steeneken, “Thermal control of tension in suspended 2D materials,” Bulletin of the American Physical Society 62 (2017a), http://meetings.aps.org/link/BAPS.2017.MAR.P32.4.
 Lee et al. (2008) Changgu Lee, Xiaoding Wei, Jeffrey W Kysar, and James Hone, “Measurement of the elastic properties and intrinsic strength of monolayer graphene,” Science 321, 385–388 (2008).
 Yoon et al. (2011) Duhee Yoon, YoungWoo Son, and Hyeonsik Cheong, “Negative thermal expansion coefficient of graphene measured by Raman spectroscopy,” Nano Letters 11, 3227–3231 (2011).
 (26) Robin J Dolleman, Samer Houri, Dejan Davidovikj, Santiago J CartamilBueno, Yaroslav M Blanter, Herre SJ van der Zant, and Peter G Steeneken, “Optomechanics for thermal characterization of suspended graphene,” Physical Review B 96, 165421.
 Zalalutdinov et al. (2001) M Zalalutdinov, Anatoli Olkhovets, A Zehnder, B Ilic, D Czaplewski, Harold G Craighead, and Jeevak M Parpia, “Optically pumped parametric amplification for micromechanical oscillators,” Applied Physics Letters 78, 3142–3144 (2001).
 Güttinger et al. (2017) Johannes Güttinger, Adrien Noury, Peter Weber, Axel Martin Eriksson, Camille Lagoin, Joel Moser, Christopher Eichler, Andreas Wallraff, Andreas Isacsson, and Adrian Bachtold, “Energydependent path of dissipation in nanomechanical resonators,” Nature Nanotechnology 12, 631–636 (2017).
 Eichler et al. (2011) A Eichler, Joel Moser, J Chaste, M Zdrojek, I WilsonRae, and Adrian Bachtold, ‘‘Nonlinear damping in mechanical resonators made from carbon nanotubes and graphene,” Nature Nanotechnology 6, 339–342 (2011).
 Singh et al. (2016) Vibhor Singh, Olga Shevchuk, Ya M Blanter, and Gary A Steele, “Negative nonlinear damping of a multilayer graphene mechanical resonator,” Physical Review B 93, 245407 (2016).
 Croy et al. (2012) Alexander Croy, Daniel Midtvedt, Andreas Isacsson, and Jari M Kinaret, “Nonlinear damping in graphene resonators,” Physical Review B 86, 235435 (2012).
 Bunch et al. (2007) J Scott Bunch, Arend M Van Der Zande, Scott S Verbridge, Ian W Frank, David M Tanenbaum, Jeevak M Parpia, Harold G Craighead, and Paul L McEuen, “Electromechanical resonators from graphene sheets,” Science 315, 490–493 (2007).
 CastellanosGomez et al. (2013) Andres CastellanosGomez, Ronald van Leeuwen, Michele Buscema, Herre SJ van der Zant, Gary A Steele, and Warner J Venstra, “Singlelayer MoS2 mechanical resonators,” Advanced Materials 25, 6719–6723 (2013).
 Sanchez and Nayfeh (1990) Nestor E Sanchez and Ali H Nayfeh, “Prediction of bifurcations in a parametrically excited duffing oscillator,” International Journal of NonLinear Mechanics 25, 163–176 (1990).
 Houri et al. (2017) S Houri, S J CartamilBueno, M Poot, P G Steeneken, H S J van der Zant, and W J Venstra, “Direct and parametric synchronization of a graphene selfoscillator,” Applied Physics Letters 110, 073103 (2017).
 Lifshitz and Cross (2008) Ron Lifshitz and MC Cross, “Nonlinear dynamics of nanomechanical and micromechanical resonators,” Review of nonlinear dynamics and complexity 1, 1–52 (2008).
 Aubin et al. (2004) Keith Aubin, Maxim Zalalutdinov, Tuncay Alan, Robert B Reichenbach, Richard Rand, Alan Zehnder, Jeevak Parpia, and Harold Craighead, “Limit cycle oscillations in CW laserdriven NEMS,” Journal of microelectromechanical systems 13, 1018–1026 (2004).
 Chen et al. (2013) Changyao Chen, Sunwoo Lee, Vikram V Deshpande, GwanHyoung Lee, Michael Lekas, Kenneth Shepard, and James Hone, “Graphene mechanical oscillators with tunable frequency,” Nature Nanotechnology 8, 923–927 (2013).
 Zande et al. (2010) Arend M van der Zande, Robert A Barton, Jonathan S Alden, Carlos S RuizVargas, William S Whitney, Phi HQ Pham, Jiwoong Park, Jeevak M Parpia, Harold G Craighead, and Paul L McEuen, ‘‘Largescale arrays of singlelayer graphene resonators,” Nano Letters 10, 4869–4873 (2010).
 Rhoads et al. (2006) Jeffrey F Rhoads, Steven W Shaw, Kimberly L Turner, Jeff Moehlis, Barry E DeMartini, and Wenhua Zhang, “Generalized parametric resonance in electrostatically actuated microelectromechanical oscillators,” Journal of Sound and Vibration 296, 797–829 (2006).
 Jinkins et al. (2015) K Jinkins, J Camacho, L Farina, and Y Wu, “Examination of humidity effects on measured thickness and interfacial phenomena of exfoliated graphene on silicon dioxide via amplitude modulation atomic force microscopy,” Applied Physics Letters 107, 243107 (2015).
 Bunch and Dunn (2012) J Scott Bunch and Martin L Dunn, “Adhesion mechanics of graphene membranes,” Solid State Communications 152, 1359–1364 (2012).
 RuizVargas et al. (2011) Carlos S RuizVargas, Houlong L Zhuang, Pinshane Y Huang, Arend M van der Zande, Shivank Garg, Paul L McEuen, David A Muller, Richard G Hennig, and Jiwoong Park, “Softened elastic response and unzipping in chemical vapor deposition graphene membranes,” Nano Letters 11, 2259–2263 (2011).
 Hanay et al. (2015) M Selim Hanay, Scott I Kelber, Cathal D O’Connell, Paul Mulvaney, John E Sader, and Michael L Roukes, “Inertial imaging with nanomechanical systems,” Nature Nanotechnology 10, 339–344 (2015).
 Amabili (2008) Marco Amabili, Nonlinear Vibrations and Stability of Shells and Plates (Cambridge University Press, 2008).
 Davidovikj et al. (2017b) Dejan Davidovikj, Farbod Alijani, Santiago J CartamilBueno, Herre SJ van der Zant, Marco Amabili, and Peter G Steeneken, “Nonlinear dynamic characterization of twodimensional materials,” Nature Communications 8, 1253 (2017b).
 Zaitsev et al. (2012) Stav Zaitsev, Oleg Shtempluck, Eyal Buks, and Oded Gottlieb, “Nonlinear damping in a micromechanical oscillator,” Nonlinear Dynamics 67, 859–883 (2012).
 Nayfeh and Mook (1995) A H Nayfeh and D T Mook, Nonlinear oscillations (Wiley, 1995).
 Doedel et al. (1998) Eusebius J. Doedel, Alan R. Champneys, Thomas F. Fairgrieve, Yuri A. Kuznetsov, BjÂ¨orn Sandstede, and Xianjun Wang, “Auto 97: Continuation and bifurcation software for ordinary differential equations,” (1998).
 Lakes (2009) R. S. Lakes, Viscoelastic Materials, Vol. 1 (Cambridge University Press, New York, USA, 2009).