Measurements

Experiments are performed on a 10 nm thick multilayer graphene nanodrum with a diameter of 5 µm, which is transferred over a cavity etched in a layer of SiO₂ with a depth of 285 nm. A blue laser is used to thermomechanically actuate the membrane, where a red laser is being used to detect the motion, using interferometry (see “Methods” for details). A schematic of the setup is shown in Fig. 1a.

Fig. 1

Nonlinear dynamic response of a graphene nanodrum near 2:1 internal resonance.

a Fabry–Pérot interferometry with thermomechanical actuation and microscope image of the graphene. Experiments are performed in vacuum at 10⁻³ mbar. Red laser is used to detect the motion of the graphene drum and the blue laser is used to optothermally actuate it. BE beam expander, QWP quarter wave plate, PBS polarized beam splitter, PD photodiode, DM dichroic mirror, VNA vector network analyzer, $$ {\mathrm{V}}_{a{c}_{in}}$$ analyzer port, $$ {\mathrm{V}}_{{\mathrm{ac}}_{\mathrm{out}}}$$ excitation port. In the device schematic, Si and SiO₂ layers are represented by orange and blue colors, respectively. b Direct frequency response curve of the device (motion amplitude vs. drive frequency), showing multiple resonances (drive level = −12.6 dBm). The mode shapes are simulated by COMSOL. Resonance peaks are associated with $$ f_{m,n}^{\left(k\right)}$$ where m represents the number of nodal diameters, n nodal circles, and k = 1, 2 stand for the first and second asymmetric degenerate modes. Dashed line shows 2f_0,1, which is the drive frequency where the parametric resonance of mode f_0,1 is activated. c Parametric resonance curves (calibrated motion amplitude vs. drive frequency), driven at twice the detection frequency. As the parametric resonance curves approach the 2:1 internal resonance (IR), f_SNB first locks to 2:1 IR frequency and consecutively saddle-node bifurcation surges to a higher frequency and amplitude. A_SNB and f_SNB stand for the amplitude and frequency of saddle-node bifurcation. d Variation of the nonlinear damping τ_nl1 as a function of drive F₁. Dashed lines represent different regimes of nonlinear damping. White region represents a constant nonlinear damping, purple region an increase in nonlinear damping in the vicinity of 2:1 IR and orange region an increase in nonlinear damping due to IR with a higher mode.

By sweeping the drive frequency, we obtain the frequency response of the nanodrum in which multiple directly driven resonance modes can be identified (Fig. 1b). We find the fundamental axisymmetric mode of vibration at f_0,1=20.1 MHz and several other modes, of which the two modes, at $$ f_{2,1}^{\left(1\right)}$$ = 47.4 MHz and $$ f_{2,1}^{\left(2\right)}$$ = 50.0 MHz, are of particular interest. This is because, to study the effect of IR on nonlinear damping, we aim to achieve a 2:1 IR by parametrically driving the fundamental mode, such that it coincides with one of the higher frequency modes. The frequency ratios $$ f_{2,1}^{\left(1\right)}/f_{0,1}\approx 2.3$$ and $$ f_{2,1}^{\left(2\right)}/f_{0,1}\approx 2.4$$ are close to the factor 2, however additional frequency tuning is needed to reach the 2:1 IR condition.

The parametric resonance can be clearly observed by modulating the tension of the nanodrum at frequency ω_F with the blue laser while using a frequency converter in the vector network analyzer (VNA) to measure the amplitude at ω_F/2 as shown in Fig. 1c. By increasing the parametric drive, we observe a Duffing-type geometric nonlinearity over a large frequency range, such that the parametrically driven fundamental resonance can be tuned across successive 2:1 IR conditions with modes $$ f_{2,1}^{\left(1\right)}$$ and $$ f_{2,1}^{\left(2\right)}$$, respectively.

In Fig. 1c, we observe that the parametric resonance curves follow a common response until they reach the saddle-node bifurcation frequency f_SNB above which the parametric resonance curve reaches its peak amplitude A_SNB and drops down to low amplitude. We note that the value of A_SNB can be used to determine the degree of nonlinear damping³³. Therefore, to extract the nonlinear damping coefficient τ_nl1 of mode f_0,1 from the curves in Fig. 1c, we use the following single-mode model to describe the system dynamics

As we gradually increase the drive level, f_SNB increases until it reaches the vicinity of the IR, where we observe an increase in τ_nl1 (Fig. 1d). Whereas f_SNB increases with parametric drive F₁, Fig. 1c shows that its rate of increase $$ \frac{\mathrm{d}{f}_{\mathrm{SNB}}}{\mathrm{d}{F}_1}$$ slows down close to $$ f_{2,1}^{\left(1\right)}$$, locking the saddle-node bifurcation frequency when f_SNB ≈ 45 MHz. At the same time, τ_nl1 increases significantly at the associated parametric drive levels, providing the possibility to tune nonlinear damping up to twofolds by controlling F₁, as seen in Fig. 1d.

Figure 1c also shows that above a certain critical parametric drive level F_1,crit, the frequency locking barrier at f_SNB ≈ 45 MHz is broken and f_SNB suddenly jumps to a higher frequency (≈5 MHz higher), and a corresponding larger A_SNB. We label this increase in the rate $$ \frac{\mathrm{d}{f}_{\mathrm{SNB}}}{\mathrm{d}{F}_1}$$ by “surge” in Fig. 1c, where an abrupt increase in the amplitude–frequency response is observed to occur above a critical drive level F_1,crit. Interestingly, even above F_1,crit a further increase in τ_nl1 is observed with increasing drive amplitude, indicating that a similar frequency locking occurs when the parametric resonance peak reaches the second IR at $$ f_{\mathrm{SNB}}\approx f_{2,1}^{\left(2\right)}$$. Similar nonlinear phenomena are also showcased in a second nanodrum undergoing parametric–direct modal interaction, confirming the reproducibility of the observed physics (see Supplementary Note II).

Theoretical model

Although the single-mode model in Eq. (1) can capture the response of the parametric resonance, it can only do so by introducing a nonphysical drive level dependent nonlinear damping coefficient τ_nl1(F₁) (Fig. 1d). Therefore, to study the physical origin of our observation, we extend the model by introducing a second mode whose motion is described by generalized coordinate x₂. Moreover, to describe the coupling between the interacting modes at the 2:1 IR, we use the single term coupling potential $$ U_{\mathrm{cp}}=\alpha x_1^2{x}_2$$ (see Supplementary Note III). The coupled equations of motion in the presence of this potential become

The two-mode model describes a parametrically driven mode with generalized coordinate x₁ coupled to x₂ that has eigenfrequency $$ {\omega }_2=2\pi f_{2,1}^{\left(1\right)}$$, damping ratio τ₂, and is directly driven by a harmonic force with magnitude F₂.

To understand the dynamics of the system observed experimentally and described by the model in Eq. (2), it is convenient to switch to the rotating frame of reference by transforming x₁ and x₂ to complex amplitude form (see Supplementary Note IV). This transformation reveals a system of equations that predicts the response of the resonator as the drive parameters (F₁, F₂, and ω_F) are varied. Solving the coupled system at steady state yields the following algebraic equation for the amplitude a₁ of the first mode

To understand the physics associated with the frequency locking and amplitude–frequency surge, we use the experimentally extracted fit parameters from the two-mode model and numerically generate parametric resonance curves using Eq. (3) for a large range of drive amplitudes (see Fig. 2a). We see that for small drive levels, an upward frequency sweep will follow the parametric resonance curve and then will lock and jump down at the first saddle-node bifurcation (SNB1) frequency, which lies close to $$ f_{\mathrm{SNB}}\approx f_{2,1}^{\left(1\right)}$$. At higher parametric drive levels, the parametric resonance has a stable path to traverse the IR toward a group of stable states at higher frequencies.

Fig. 2

Parametric–direct internal resonance.

a Color map of the analytical model response curves obtained by using the fitted parameters from experiments. Colors correspond to frequency response (motion amplitude vs. drive frequency) solutions with a certain parametric drive level. Black lines show samples from these solutions where solid lines are stable and dashed lines are unstable solutions. White dashed line is where parametric resonance meets with interacting mode and undergoes internal resonance. b The underlying route of the amplitude–frequency surge is revealed by tracing the evolution of saddle-node bifurcations (green and red squares represent theoretical SNB1 and SNB2, whereas experimental SNB1 is represented by crosses) of the parametric resonance curves.

A more extensive investigation of this phenomenon can be carried out by performing bifurcation analysis of the steady-state solutions (see Supplementary Note IV). The bifurcation analysis reveals two saddle-node bifurcations near the singular region of the IR, one at the end of the first path (SNB1) and another at the beginning of the second path (SNB2) (Fig. 2b). As the drive amplitude increases, the bifurcation pair starts to move toward each other until they annihilate one another to form a stable solution at the connecting point, which we labeled as “surge.” It is also possible to observe that the rate at which saddle-node pairs approach each other dramatically drops near the IR condition, demonstrating the “locking” which we also observed in the experiments.

To check how closely the two-mode model captures the variation of τ_nl1 close to the IR condition, we follow a reverse path, and fit the numerically generated resonance curves of Fig. 2a using the single-mode model of Eq. (3) with τ_nl1 as the fit parameter. In this way, we track the variation of τ_nl1 in the single-mode model with the parametric drive F₁, similar to what we observed experimentally and reported in Fig. 1c. The result of this fit is shown in Fig. 3a, where a similar anomalous change of nonlinear damping is obtained for the two-mode model.

Fig. 3

Measurements and fits of the effective nonlinear damping.

a Variation of the effective nonlinear damping parameter (τ_nleff) with respect to parametric drive. The τ_nleff is obtained by fitting the numerically generated curves of Fig. 2a as the fit parameter. Dashed lines represent different regimes of nonlinear damping. White regions represent a constant nonlinear damping and purple region represents an increase in nonlinear damping in the vicinity of 2:1 IR. b Comparison of the ratio between linear damping (τ₁) and total damping (τ_tot). In the figure, blue and red dashed lines represent τ₁/τ_tot obtained from uncoupled and coupled models, whereas black crosses represent experiments.

The variation of nonlinear damping affects the total damping (sum of linear and nonlinear dissipation) of the resonator too. It is of interest to study how large this effect is. In Fig. 3b, we report the variation in the ratio of the linear damping τ₁ and the amplitude-dependent total damping τ_tot = (ω₁/Q + 0.25τ_nleff∣x₁∣²)³³ in the spectral neighborhood of $$ f_{2,1}^{\left(1\right)}$$, and observe a sudden decrease in the vicinity of IR. This abrupt change in the total damping is well captured by the two-mode model. With the increase in the drive amplitude, τ₁/τ_tot values of this model though deviate from those of the experiments due to a subsequent IR at $$ f_{2,1}^{\left(2\right)}/f_{0,1}\approx 2.4$$ that is not included in our theoretical analysis. The dependence of τ₁/τ_tot on frequency shows that near IR, the total damping of the fundamental mode increases nearly by 80%. We note that increased nonlinear damping near IR was also observed in Güttinger et al.¹¹. In that work, nonlinear damping was studied using ringdown measurements, with two modes brought close to an IR by electrostatic gating. The increased nonlinear damping was attributed to a direct–direct 3:1 IR, which as shown theoretically in Shoshani et al.²⁴ leads to a high order (quintic) nonlinear damping term. Conversely, in our work, two modes are brought into parametric–direct 2:1 IR by adjusting the parametric drive level. This results in a nonlinear damping term that already comes into play at smaller amplitudes because it is of lower (cubic) order, as discussed in Shoshani et al.²⁴. Moreover, the nonlinear damping mechanism in Güttinger et al.¹¹ is approximately described by two exponential decays with crossovers from (τ₁ + τ₂)/2 to τ₁, which implies that similar to Shoshani et al.²⁴, τ₂ > τ₁ is required to observe positive nonlinear damping. This is in contrast with the damping mechanism we describe, where the effective nonlinear damping actually increases for smaller τ₂ (see Eq. (3)).

Sample fabrication

Devices are fabricated using standard electron-beam (e-beam) lithography and dry etching techniques. A positive e-beam resist (AR-P-6200) is spin coated on a Si wafer with 285 nm of thermally grown SiO₂. The cavity patterns ranging from 2 to 10 µm in diameter are exposed using the Vistec EBPG 5000+ and developed. The exposed SiO₂ are subsequently etched away in a reactive ion etcher using a mixture of CHF₃ and Ar gas until all the SiO₂ is etched away and the Si exposed. Graphene flakes are then exfoliated from natural crystal and dry transferred on top of cavities.

Laser interferometry

The experiments are performed at room temperature in a vacuum chamber (10⁻³ mbar). A power modulated blue laser (λ = 405 nm) is used to thermomechanically actuate the nanodrum. The motion is then readout by using a red laser (λ = 633 nm) whose reflected intensity is modulated by the motion of the nanodrum in a Fabry–Pérot etalon formed by the graphene and the Si back mirror (Fig. 1a). The reflected red laser intensity from the center of the drum is detected using a photodiode, whose response is read by the same VNA that modulates the blue laser. The measured VNA signal is then converted to displacement in nanometers using a nonlinear optical calibration method³⁵ (see Supplementary Note VII).