Self-sensing, tunable monolayer MoS2 nanoelectromechanical resonators

Excellent mechanical properties and the presence of piezoresistivity make single layers of transition metal dichalcogenides (TMDCs) viable candidates for integration in nanoelectromechanical systems (NEMS). We report on the realization of electromechanical resonators based on single-layer MoS2 with both piezoresistive and capacitive transduction schemes. Operating in the ultimate limit of membrane thickness, the resonant frequency of MoS2 resonators is primarily defined by the built-in mechanical tension and is in the very high frequency range. Using electrostatic interaction with a gate electrode, we tune the resonant frequency, allowing for the extraction of resonator parameters such as mass density and built-in strain. Furthermore, we study the origins of nonlinear dynamic response at high driving force. The results shed light on the potential of TMDC-based NEMS for the investigation of nanoscale mechanical effects at the limits of vertical downscaling and applications such as resonators for RF-communications, force and mass sensors.

The arrow shows the direction of increasing built-in strain. We note that increasing the built-in strain results in higher , while the sensitivity to the applied gate voltage decreases, which is expected as the electrostatically induced strain becomes a smaller contribution to the total strain in comparison to the built-in strain. Additionally, the curvature of the tuning curve changes from convex shape for small built-in strain to a concave shape as the built-in strain increases. (b) Resonant frequency plotted as a function of gate voltage for increasing mass density (the direction of the arrow). The higher the mass density, the lower , would be and the tunability decreases. In contrast to the effect of built-in strain, the curvature of the tuning curve does not change for higher mass density. The higher the built-in strain, the lower the sensitivity to the applied gate voltage. Also, the curvature of the tuning curve changes with increasing the built-in strain. (b) Resonator with a higher mass density have a lower resonant frequency, and, the heavier the resonator gets, the less sensitive it becomes to the applied gate voltage. The curvature of the tuning curve is independent of the mass density. (c) The mode shapes and resonant frequencies of the first three flexural modes for the same geometry as R2. (d) and (e) show the FEM simulation of the piezoresistive and capacitive transduction of mechanical vibrations. In order to compare with the calculation in Supplementary Note 4, we modeled resonators with the geometry, built-in strain and mass density as extracted for device R1. The results show that the piezoresistive transduction mechanism is expected to be a factor of 4 stronger than the capacitive transduction mechanism. Figure 4. The effect of nonlinear restoring force on the frequency response of doubly clamped MoS2 resonators, and the dependence of resonant width on the nonlinear damping mechanism (a) Schematic illustration of the frequency response of a Duffing resonator for 0. For the small drive the peak has a symmetric shape with a maximum at (frequency of linear harmonic oscillator). As the drive increases, the peak shape becomes asymmetric. The schematic illustration of the response corresponding to the highest drive (red), illustrates the appearance of bistability in the frequency response. (b) The frequency response of the mixing current for device R2 from the main text. The local gate voltage 1 V. The dashed red line shows the Lorentzian fit to the resonant peak. The extracted values of resonant width and quality factor are 1.6 MHz and 75, respectively. (c) The resonant width and the quality factor as a function of local gate voltage, measured for asymmetric peaks of device R2 ( Figure 4a from the main text). The red markers are the experimental results and the black line is a fit to the model described by Eichler et al. 12 Supplementary Figure 5. Peak current as a function of input RF power. The experimentally measured (red markers) peak current increases with the input RF power, and is proportional to the square root of the input RF power (black line). The mixing current mapped as a function of the driving frequency and the local gate voltage. The dashed lines are guides to the eye. The resonant peak from the monolayer MoS2 membrane is defined by the applied tension, and thus, is tunable by changing the applied electrostatic strain. In contrast, the resonant behavior of the thick metal contacts is in the plate limit, and therefore, the resonant frequency is independent of the applied tension (local gate voltage). (c) Further increasing the local gate voltage increases the resonant frequency of the MoS2 vibrations to surpass the resonant frequency of the suspended contacts. For gate voltages up to 12.5 V the resonant frequency of MoS2 membrane remains tunable and the resonant frequency of the metal contacts remains independent of the applied electrostatic strain. (d) The frequency response of the mixing current for 10.5 V. The resonant peak from the MoS2 membrane is now at a higher frequency compared to the resonant frequency of the suspended contacts. Additionally, as a result of high drive amplitude due to the high , the peak becomes asymmetric which confirms that the MoS2 membrane has entered the nonlinear regime.

Supplementary Note 1: Device fabrication
We have fabricated suspended channel MoS2 devices featuring a local bottom gate. The presence of local gate minimizes the overlap area between the gate and the source/drain electrodes, and reduces the source-gate and the drain-gate capacitances. This is important since at higher frequencies the impedance of the parasitic capacitances between the gate and the source/drain electrodes are no longer negligible and could introduce filtering effects on the output electrical signal.
The fabrication starts by depositing local gate electrodes on an intrinsic Si substrate covered by 270 nm of oxide. Subsequently, 250 nm of Al2O3 is deposited on top of the gate electrodes, using atomic layer deposition (ALD), which acts as the gate dielectric. In the next step, single crystalline domains of CVD grown monolayer MoS2 are transferred on top of the local gates. The material synthesis is described elsewhere. 1 Using CVD materials allows controlling the reproducibility of device properties and increases the throughput of the fabrication process.
The MoS2 flakes were then shaped into nanoribbons using a polymethyl methacrylate (PMMA) mask patterned in an electron-beam lithography (EBL) step followed by a reactive ion etching process in oxygen plasma and mask removal in acetone. Metal electrodes (Cr/Au) were deposited by standard bilayer technique (methyl methacrylate (MMA)/PMMA) in electron-beam lithography (EBL), followed by thermal evaporation and lift-off in acetone. The metal contacts where then thermally annealed in Ar at 200 ºC for 2 hours.
Wet etching using orthophosphoric acid undercuts the Al2O3 to yield suspended MoS2 channel. The suspended devices were successively transferred to water and ultrapure ethanol, followed by critical point drying (CPD) to prevent the suspended structure from collapsing due to capillary forces and surface tension. Prior to electromechanical measurements, samples are imaged using atomic force microscope (AFM) in tapping mode to verify their suspension. (Figure 1b from the main text).

Supplementary Note 2: DC transport characterization of the devices
Following the realization of suspended structure, the resonators are bonded on a chip carrier and loaded in high vacuum where all the electrical/electromechanical characterization is carried out. The devices are then annealed at 170℃ for more than 24 hours to remove water and other adsorbates from both surfaces of the suspended 2D channel.
After annealing and prior to RF electromechanical measurements, we characterize the electrical transport behavior of resonant channel transistors in the DC regime using Agilent E5270B parameter analyzer. Supplementary Figure 1 shows the DC transfer characteristics for devices R1and R2 from the main text.
The transfer characteristics show a relatively low conductivity which is expected since the devices are operating in the subthreshold regime. We avoid higher values of gate voltage to prevent the collapse of suspended structure before electromechanical characterization.

Supplementary Note 3: The RF electromechanical characterization setup
During the DC characterization of the suspended channel transistors, the chip carriers are already mounted on a holder which is connected to SMA connectors via a printed circuit board. The RF signals are transmitted using high frequency coaxial cables and all the connectors on the cables and vacuum feedthrough are of the SMA type with a frequency range up to 3 GHz and a voltage standing wave ratio (VSWR) of less than 1.3. We used as a microwave signal generator Keysight N5183B and the mixed-down current was detected using the current input of a Stanford Research Systems SR830 DSP lock-in amplifier. The DC gate voltage was applied using a Keithley 2450 source measure unit. The matching load on the signal input consists of a 50 Ω resistor and 100 pF capacitor connected in series. The 100 nF capacitor connected to the signal output provides an RF ground and eliminates fluctuations in the frequency response of the mixing current due to the floating output. Figure 1c in the main text shows the circuit diagram for the RF measurement setup.

Detection technique
The frequency response is detected indirectly, using a mixing technique which facilitates the detection by transferring the output signal to lower frequencies without the loss of the information stored in its amplitude. To implement the signal mixing we choose frequency modulation (FM) over amplitude modulation (AM) or 2-source techniques, because FM only detects the terms with electromechanical origin and suppresses the purely electrical terms. 2 Another advantage of FM or AM or 2-source technique is that the absence of background signal results in enhanced signal-to-background ratio and signal-to-noise ratio. 2 We apply to the source electrode, an FM modulated signal in the form of where is the amplitude of the applied voltage, is the carrier frequency in the RF range, t is time, 2π ⁄ is the frequency deviation (typically 50 kHz to 100 kHz), and 2π ⁄ (typically 1 kHz) is the reference frequency at which the demodulated signal is detected using the lock-in technique.

Frequency deviation
is proportional to the amplitude of the detected signal and could be adjusted to get a reasonable signal to noise ratio. On the other hand, high values of result in the contributions from higher-order terms to be detected in the measured signal 2 which leads to nonlinearity and broadening of the peak shape that does not originate from the physics governing the mechanical behavior of the resonator, but is due to the nonlinearities in the actuation and detection mechanisms. Additionally, it is crucial to choose values of and that result in a transient regime that disappears before being detected. This allows detection of the steady state response which is the regime of interest in the behaviour of NEMS resonators. Both concerns could be addressed by choosing ω and that are much smaller than the resonant width Δω. 2 This is especially important for studying the nonlinear behavior of our resonators in the absence of contributions from additional sources of nonlinearity that come from the experimental detection technique. In practice, it is possible to find the appropriate range for by monitoring the evolution of the peak shape with increasing . 2

Supplementary Note 4: Comparison between piezoresistive transduction and electrostatic transduction
When using FM detection, the source-drain current is composed of several high-frequency components ( , 2 ,…), a DC component and a mixing component at , with the latter being the only electromechanical term in and the only term that will be detected by the lock-in amplifier.
The mixed-down current at frequency is proportional to | / |, with the amplitude of the mechanical motion of the resonator. Gouttenoire  where is the conductance of the device and translates the mechanical vibrations into an electrical signal and, therefore, it is defined by the transduction mechanism.
The traditional capacitive transduction relies on the capacitive coupling between the suspended membrane and the local gate electrode. The vibrations modulate the gate capacitance (Cg), which in turn modulates the charge carrier density in the semiconducting MoS2 channel and results in the modulation of the conductance. Because of the oscillating conductance, the mechanical vibrations of the resonator is transduced into electrical signals. Therefore, the capacitive transduction mechanism can be formulated as: 2 Considering that: 3 where is the membrane area and is the charge carrier mobility, and since is constant for the vibrating membrane: 4 We estimated the field-effect mobility from the DC transfer characteristics of the device R1 (Supplementary Figure 1a) to be ~ 0.42 cm 2 V −1 s −1 . This relatively low value of mobility is expected since the device is operating in the sub-threshold regime for the range of used to measure the transfer characteristics. We avoided higher values of for DC characterization in order to minimize the risk of membrane collapsing due to the electrostatic force from the gate.
Using the extracted and ~ 17 pF m (derived assuming a parallel plate capacitance), we can estimate for the measured peak presented in Figure 2a in the main text, ~ 8.8 10 -4 S m -1 .
Alternatively, thanks to the presence of strong piezoresistivity in monolayer MoS2 3 , an additional transduction mechanism is available in our resonators. Because of the geometry of a doubly clamped thin membrane, the mechanical vibrations induce oscillating strain due to the periodic elongation of the membrane. The oscillating strain in turn modulates conductivity of the MoS2 membrane, meaning that, the piezoresistive transduction mechanism can be extracted from: where ~ 3 nS is the initial conductance of the MoS2 channel when it is not vibrating and is estimated from the DC transport characterization. We approximate the shape of the first mode with two straight lines resulting in ~ . Using the average value of for monolayer MoS2 3 and equal to the vibration amplitude at the onset of nonlinearity (Supplementary Note 10), we derive ~ 3 10 -3 S m -1 .
We find that for R1, the piezoresistive transduction is more than 3 times stronger than the capacitive transduction: The piezoresistive transduction mechanism is not available in semimetallic graphene where the piezoresistive gauge factor is two orders of magnitude smaller compared to MoS2 due to the absence of a bandgap. 4

Supplementary Note 5: Finite element modeling of atomically thin MoS2 resonators
Finite element models are built using COMSOL 5.3 to investigate the dynamic mechanical behavior of MoS2 membranes and the effects of built-in strain, mass density and electrostatically induced strain (via an external gate voltage). Additionally, we simulated the electrical response of the MoS2 NEMS resonators to compare the piezoresistively transduced response with the capacitively transduced one.
A multiphysics simulation is implemented that would allow the study of the interplay between electrostatic actuation, static and dynamic mechanical response, as well as, the electrical currents resulting from the transduction of dynamic mechanical response.
Devices R1 and R2 are modeled as doubly clamped membranes suspended over an airgap. The membrane geometry is extracted from the AFM topography image in non-contact mode. The AFM images also give the depth of the airgap making it possible to implement the equivalent gate capacitance corresponding to the thicknesses of air dielectric and the remaining Al2O3 dielectric.
The material parameters used for MoS2 are Poisson's ratio = 0.27 and Young's modulus = 270 GPa. 5 It is important to notice that even though there is resist residue on top of the MoS2 membrane, the Young's modulus of MoS2 is at least two orders of magnitude larger than that of polymer resists and, therefore, it dominantly defines the equivalent stiffness for the combination of MoS2 and resist adsorbates. The relative dielectric constants of 1 and 9 are used for air and for Al2O3, respectively. 6 MoS2 is modeled as an isotropic linear elastic material since the applied strain levels are below the range where the nonlinear elastic behavior becomes considerable. 5 Dielectrics are modeled as linear elastic dielectrics with the airgap defined as a nonsolid material that can deform in the direction normal to the membrane. The membrane is meshed using maximum element size of 20 nm in the in-plane direction and 0.5 nm in the vertical direction. The dielectric material is meshed using maximum element size of 20 nm in the in-plane direction and 50 nm in the vertical direction.
The mechanical boundary conditions are set in a way that the membrane is clamped at the two edges. The lower boundary of the dielectric material is connected to a voltage terminal which applies a DC gate voltage. A small-signal perturbative voltage is applied between the clamped edges of the membrane to simulate the AC source-drain voltage.
In order to study the effect of built-in strain, mass density and electrostatically induced strain, we use an electromechanics multiphysics (solid mechanics and electrostatics) package and solve the model with a stationary analysis followed by an eigenfrequency analysis.
From the stationary study, we determine the static vertical deflection and the static strain distribution induced by the electrostatic DC force. The static strain obtained from FEM simulations is in agreement with the analytically calculated strain when approximating the shape of the first mode with two straight lines: Next, we simulate the piezoresistive and capacitive transduction mechanisms by adding electric currents to the above mentioned multiphysics. The model is, first, solved using a prestressed analysis to obtain the stationary response to the DC gate force. Once the stationary response is obtained a frequency domain analysis solves the dynamic behaviour for a defined frequency range.
To simulate the piezoresistive transduction, the conductivity of the membrane at each point is defined by the strain distribution according to: 9 here is the conductivity of the membrane when it is not vibrating, is the piezoresistive gauge factor for monolayer MoS2, 3 is the mechanical strain induced by the vibrations of the membrane and is the parameter that couples the mechanical and electrical parts of the multiphysics package.
On the other hand, the capacitive transduction mechanism is modelled by defining the conductivity of the membrane as: 10 here represents the capacitive coupling to the gate terminal which is modulated with the mechanical vibrations of the membrane and is the parameter that couples the mechanical and electrical parts of the multiphysics package.
By implementing the models corresponding to each transduction mechanism, we simulate the vibration induced modulation of conductivity. In the case of piezoresistive transduction, the conductivity changes periodically as a result of vibration-induced strain, while in capacitive transduction the vibrations change the gate capacitance which in turn changes the charge carrier density in the channel.
The solution to the frequency domain analysis determines the current density. The total current is derived by integrating the current density over the cross section of the device. In order to compare the results with the calculations in Supplementary Note 4, the simulated frequency response of device R1 is shown in Supplementary Figure 3d and Supplementary Figure 3e for the piezoresistive transduction and capacitive transduction, respectively. The simulations predict that the piezoresistive transduction is ~4 times stronger than the capacitive transduction, which is in line with the results obtained by analytical estimations in Supplementary Note 4.

Supplementary Note 6: Strain-modulated electromechanical response of the atomically thin MoS2 resonators
The voltage applied across the gate capacitance results in the application of a force on the suspended resonator which is written as:

11
where is the spatial derivative of the gate capacitance, is the static voltage applied to the local gate electrode and is the AC voltage applied between the source and drain electrodes of the device. Linearization of the force for small results in where is the electrostatic DC force applied on the membrane and is described as: and is the driving force that actuates the mechanical vibrations and is described as: While the electrostatic force results in a static deflection of the membrane from its rest position, the driving force leads to a time-varying bending profile.
The dynamic behavior of resonators can be modeled using a common approximation to the equation of motion for a doubly clamped resonator, given by the Euler-Bernoulli Equation 7 . A perturbative treatment of the Euler-Bernoulli equation results in a continuum mechanics model for a doubly clamped resonator with the resonant frequency can be described as 7,8 :

14
Where is the resonant frequency of the fundamental mode, , and are the length, width and thickness of the resonator, is the Young's modulus, is the built-in tension in the resonator, and is the mass density. Supplementary Equation 14 shows that the resonant frequency of a doubly clamped resonator is a function of both its geometry and its intrinsic material properties. Two limits of the dynamic mechanical behavior happen when either one of the arguments under the square root become dominant over the other. For a thick resonator the first argument is dominant, the resonant frequency could be described as /12 and the resonator behaves in the so called "plate limit". The other extreme is the "membrane limit" where the resonator is very thin and Supplementary Equation 14 is reduced to / , being the effective sheet density. Therefore, for resonators in the membrane limit, the resonant frequency is predominantly set by the tension rather than the bending rigidity (∝ of the material.
In this work, in order to normalize for the different ribbon geometries, we use built-in strain in place of built-in tension. The built-in strain is derived from the built-in tension using assuming = 270 GPa. 5 Monolayer MoS2 resonators present the ultimate limit of membrane thinness and are well in the membrane limit. Therefore, the presence of the built-in strain on the as-fabricated resonators accounts for the high resonant frequencies observed, even in the absence of any electrostatic DC force. In addition to the built-in strain, the gate voltage induced results in additional electrostatic strain in the resonator. As a result of the additional strain, the resonant frequency shifts towards higher values. The shift in the resonant frequency could be described by the function Δ: Using the built-in strain and the mass density as the fitting parameters, we fit this model to the experimentally measured resonant frequency as a function of local gate voltage for 7 devices. The deduced values of mass density and built-in strain are shown in the Figure 3c in the main text. The limit to the tunability of the resonant frequency is set by the pull-in voltage at which is large enough for the suspended resonator to snap in and collapse.

Supplementary Note 7: The equation of motion and the frequency response of a driven, damped Duffing resonator
Previous studies on the mechanical properties of MoS2 have reported that the intrinsic strainstress relation becomes nonlinear with increasing deformation. 5 However, in the case of our resonators, the nonlinear behavior appears at strain levels below the nonlinear regime in intrinsic stress-strain relation of MoS2 5 and could be explained by the amplitude-dependent nonlinear effects in the equation of motion.
While for low driving force the resonators operate in the linear regime, increasing the driving force and the motional amplitude result in appearance of nonlinear response. In this regime, due to the increased , the nonlinear restoring force ( ) and nonlinear damping (η should be taken into account and the equation of motion is described by: 16 known as the van der Pol-Duffing Equation 10 Here, is the motional amplitude, is the effective spring constant, m the resonator mass, is the linear damping coefficient, is the nonlinear damping coefficient, is the Duffing force coefficient, is the amplitude of the driving force, is the driving frequency and is time. Using secular perturbation theory, Lifshitz et. al. have calculated the solution of the van der Pol-Duffing equation in the limit of week linear damping (Q smaller than 1000). 11 The obtained expressions for the steady state motional amplitude and phase are: The Duffing force coefficient could originate from two main sources; the nonlinearities arising from the presence of an external electrostatic force and the nonlinearities arising from geometric effects.
The external electrostatic drive force, applied between the resonator and the local gate electrode, results in the deviation of the resonator's equilibrium position from its rest position. The deviation from the rest position results in a cubic restoring force in the equation of motion ( ) which acts in the opposite direction of the linear elastic force ( ). Therefore, the resonator is softened and the resonant frequency reduces. 11 On the other hand, geometrical nonlinearities are expected to appear in doubly clamped atomically thin membranes where the transverse vibrations results in non-negligible elongation and stretching. In the case of geometric nonlinearities, Lifshitz et. al. have used a perturbative solution to the Euler-Bernoulli equation and derived the Duffing parameter, for the fundamental mode, which is described by for a high-tension, thin, doubly clamped resonator. 11 Using the parameters we have derived for our representative device R2, we obtain ~ 1.5 10 kg m s . Furthermore, the positive Duffing coefficient leads to the stiffening and increase in the resonance frequency which is in agreement with our observations indicating that, in the case of MoS2 resonators reported in this work, the Duffing force originates from geometric nonlinearities.
Supplementary Figure 4a shows a schematic illustration of the resulting frequency response for increasing drive amplitude. The response to the small drive resembles the Lorentzian response of a driven harmonic oscillator with a symmetric peak at , and it becomes asymmetric as the drive increases. For a strong driving force, the Duffing equation could have two stable solutions referred to as the upper and the lower stable solutions. 11 As the drive frequency is swept from lower to higher values, the experimentally measured response follows the upper stable solution until it reaches the maximum point, and then falls to the lower stable solution which emerges as the abrupt change of the motional amplitude (instead of following the dashed line in Supplementary Figure 4a).
In practice, the condition for the existence of two stable solutions (bistability) is given by 11 : √3 20 This allows extracting an upper limit for the nonlinear damping coefficient that would allow the bistability. Using the estimated value of ~ 1.5 10 kg m s , we obtain ~4 10 kg m s for device R2, which is comparable with the values that have been reported for resonators based on carbon nanotubes and graphene. 12 Bistability is also expected to give rise to hysteresis in the frequency response of Duffing resonators. However, the presence of thermal noise at room temperature, shrinks the size of the hysteresis loop and makes it undetectable in our measurements. 11 The effect of nonlinear damping ( ) appears as increased dissipation and broadening of the resonance peak. As long as the resonance bandwidth Δω is much smaller than the resonance frequency, the standard definition for the bandwidth is given by the ratio between the average energy loss rate and the average stored mechanical energy: where 〈 〉 is the time-averaged stored mechanical energy over a few oscillation periods and 〈 〉 is the average energy loss rate over the same period of time.
For a linear harmonic oscillator, the solution to the equation of motion results in .
Therefore, the bandwidth could be written as Δω which is independent of the amplitude of the mechanical vibrations, and thus, independent of the driving force. In this regime, the frequency response of the mixing current shows a symmetric lineshape, and, the peak frequency is the same for the mixing current as the stored mechanical energy. The resonators enter the nonlinear regime when the driving force is increased. In the nonlinear regime, the effective damping is a superposition of the linear and nonlinear damping expressed as . 11 It is possible to estimate the gate voltage at which the influence of nonlinear damping becomes considerable by looking at the responsivity defined as: Here is the maximum of mixing current. The resulting responsivity is plotted as a function of in Figure 4b in the main text with a peak observed around = 4 V which indicates the gate voltage at which the contribution from nonlinear damping to the effective dissipation become considerable.
At sufficiently high drive amplitudes, where the nonlinear damping becomes appreciable, the resonance width and the quality factor show strong dependence on the drive amplitude. Supplementary Figure 4c shows the resonance width as a function of the local gate voltage for R2. We used a model derived by Eichler et. al. to fit the experimental results and extracted the quality factor (right hand axis on Supplementary Figure 4c) for the asymmetric peaks, assuming a regime where 12 which holds true for 2.75 nm. The observed nonlinear damping is consistent with observations in carbon nanotubes and graphene resonators 12 and confirms that in contrast to larger resonators (from meter-scale to tens of nanometer-scale), for resonators with atomic-scale transverse dimensions, damping becomes highly dependent on the amplitude of mechanical vibrations. Further investigation of the nonlinear behavior in MoS2 resonators allows eliminating it when it is unwanted and efficiently using it when it is desirable. Additionally, the observed nonlinear behavior sheds light on the potential of MoS2 resonators for studying nonlinear mode coupling between various vibrational modes, internal resonances and energy transfer from one mode to another. 8

Supplementary Note 9: Dynamic range and mass resolution
Presence of nonlinear behavior in MoS2 NEMS resonators makes it relevant to extract the linear dynamic range of these resonators. The critical mechanical amplitude at which the nonlinear behavior sets in, can be estimated from the equation 8,13 :

23
here t is thickness, Q is the quality factor in the linear regime and T the tension experienced by the resonator at the onset of the nonlinearity (extracted from the resonant frequency for the corresponding ).
The 1 dB dynamic range (DR) is given by 13 :

24
with Δ being the measurement bandwidth and S the spectral density of the thermomechanical noise at resonant frequency given by 14 Comparing the values derived for R1 and R2 reveals that the higher values of built-in strain result in enhanced critical amplitude and improve the dynamic range. This observation is consistent with the theoretical predictions 16 and suggests that strain could be used to reduce the nonlinear behavior of NEMS resonators which is an inevitable by-product of miniaturization.

Supplementary Note 10: Effect of input RF power on the frequency response
As a result of increasing the input RF power and amplitude dependent tension, the frequency at which the maximum of the mixing current occurs ( ) deviates from the frequency of linear harmonic oscillator ( ). Lifshitz et. al. use the secular perturbation theory, to show that this deviation is described by with the maximum of the mechanical motion amplitude, α the Duffing force coefficient and the resonator mass. 11 Combining Supplementary Equations 13, 17 and 27, the dependence of on input RF power is derived as: The black line in Figure 4d in the main text is a fit to this power dependence.
Additionally, the peak value of the mixing current increases by increasing the input RF power. According to Supplementary Equation 1 the dependence of the mixing current on the input RF power is described as: ∝ 29 Supplementry Figure 5 shows the experimentally measured peak current as a function of the input RF power (red markers) and the fit to (black line).

Supplementary Note 11: Effect of the suspended contacts
Releasing the MoS2 channel by undercutting the oxide layer in a wet etching step offers the advantage that the distance between the suspended MoS2 and the substrate is uniform under all the channel area. This is due to the fact that in the lateral interface between MoS2 and the oxide, the acid diffuses rapidly and Al2O3 is etched at the same rate in all the interface area. 17 The MoS2-Al2O3 interface area includes the region of the MoS2 flake which is contacted by metal electrodes, therefore, the metal electrodes in the clamping region are also suspended.
Due to the finite rigidity of the metal contacts, vibrations of the MoS2 membrane are extended in the suspended region of the contacts. The resonance of the undercut metal electrodes result in the appearance of additional peaks in the frequency response of the mixing current. These peaks are easy to identify since they show lower quality factors and their resonant frequency is not sensitive to the gate voltage. Supplementary Figure 7a shows the frequency response for a 1.1 µm long and 92 nm wide MoS2 resonator at = 2.5 V. The prominent peak observed near 123 MHz corresponds to the mechanical resonance of the MoS2 channel and has a width of ~0.6 MHz. The broader peak observed at 129 MHz is attributed to the mechanical vibrations of the suspended portion of the electrodes and shows a width ~5 MHz. Supplementary Figure 7b shows the mixing current mapped as a function of the local gate voltage and the driving frequency. The dashed lines serve as guides to the eye. The difference between the tunability of the MoS2 resonance and the electrode's resonance with further confirms the origins of the two peaks. The MoS2 resonators operate in the membrane limit where the resonant frequency is tunable via the gate voltage induced static strain, while the thick metal electrodes behave in the plate limit where the resonant frequency is predominantly defined by the bending rigidity and is independent of the applied strain. Once it shifts close to the resonant frequency of the contacts, the resonant frequency of the MoS2 membrane could not be detected in the background. By further increasing the local gate voltage, the resonant frequency of the MoS2 membrane surpasses the resonant frequency of the suspended contacts and could be detected again (Supplementary Figure 7c). Supplementary Figure 7d shows the frequency response of the mixing current for = 10.5 V. The MoS2 resonant peak is then occurring at frequencies above the resonant frequency of the suspended contacts. Furthermore, the higher value of drive amplitude in Supplementary Figure 7d  Behaving in the plate limit due to their relatively large thickness (typically around 80 nm) and bending rigidity, the resonant frequency of the suspended gold contacts is proportional to 1/ , being the length of the suspended part of the electrode (Supplementary Note 6). The suspended part of the electrode is mainly in the clamping region meaning that ~ width of the MoS2 membrane. Due to the narrow width of the ribbons measured in our experiment, the resonance from the suspended metal contacts rarely appeared in the frequency range where we observed the electromechanical resonance from the MoS2 NEMS resonators. In the two cases where we have observed such peaks, the peaks where at 129 MHz and 114 MHz for devices of width 92 nm and 98.5 nm respectively, which is consistent with the resonant frequency being proportional to 1/ .