In Fig. 2(a), the blue laser is driven at 2 f, while detecting the photodiode signal at f. When increasing the blue laser driving voltage V ac,in a remarkable effect is observed. One-by-one, the parametric resonances of graphene appear, up to 7 different modes. Each mode reaches resonance at a different threshold driving amplitude V ac,in , due to differences in quality factor and the frequency dependence of the parametric driving parameter δ40. The experiment is repeated on a different drum in Fig. 2(b). Interestingly, in this case overlap between parametric resonances is observed at high driving levels. When overlap occurs, a direct transition between the high-amplitude solution of two adjacent parametric resonances is observed, e.g. at V ac,in = 382.7 mV (RMS) between the second and third resonance. Moreover, in some cases transitions between the high-amplitude and low-amplitude solutions are observed, e.g. at V ac,in = 489.6 mV (RMS) between the same 2 modes. This random process is attributed to a strong dependence of the basin of attractions of the parametric high-amplitude and low-amplitude solutions on the initial conditions41. 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 (V ac,in = 576.2 and 707.1 mV (RMS)).

Figure 2 Multi-mode response of a parametrically driven graphene resonators. (a) Waterfall plot of the multimode response at different driving amplitudes. Each mode appears at different driving levels due to variations in quality factor and effective driving force between them. The scale bar indicates the root mean square value (RMS) of V ac,out and the labels on the right indicate the RMS driving amplitude V ac,in . (b) Waterfall plot for a different drum, showing more mechanical modes and modal interactions. (c) Forward and backward frequency sweep at the highest parametric driving amplitude for the drum in (b) revealing 14 distinct mechanical modes in parametric resonance. Full size image

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. 2(c), one can observe the period doubling bifurcations of the resonator. As many as 14 parametric resonances are observed in this system, whereas previously only 7 modes could be excited in cryogenic environments14.

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. 3(a,c)). Sweeping the frequency forward (Fig. 3(a)) and backward (Fig. 3(b)) results in a hysteresis, that grows as the driving level is increased. This is typical for the geometric nonlinearity of the Duffing-type, 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 V ac,out is detected at half of the driving frequency V ac,in . Similar to the directly driven case, a hysteresis occurs between the forward (Fig. 3(b)) and backward (Fig. 3(d)) sweeps in frequency. Below an RMS drive amplitude of 0.11 mV, δ < δ t and no response is observed. For δ > δ i the parametric resonance obtains two stable phases of resonance separated by 180 degrees13, an additional measurement that measures this behavior is shown in the Supplementary Information S2.

Figure 3 Frequency response of the fundamental mode to direct and parametric drive, for forward and backward frequency sweeps. (a) Direct drive with the frequency swept forwards. (b) Parametric drive with the frequency swept forwards. Below a driving threshold near V ac,in ≈ 0.11 mV (RMS) no mechanical response is observed. (c) Direct drive with the frequency swept backwards. (d) Parametric drive with the frequency swept backwards. Full size image

Figure 4(a–d) 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 degree-of-freedom model is derived that describes the motion of the resonator (see Supplementary Information S5) and this is fitted to the response curves in Fig. 4(a–d) (see Supplementary Information S6). The model is a combination of the Duffing, van der Pol and Matthieu-Hill equation also used in other works32,34,42,43:

$$\ddot{x}+\mu \dot{x}+

u {x}^{2}\dot{x}+(\beta +\delta \,\cos \,\omega t)x+\gamma {x}^{3}=F\,\cos \,\omega t,$$ (1)

where x is the displacement (which is approximately proportional to V ac,out ), μ is the damping coefficient, ν the nonlinear damping coefficient, β the linear stiffness coefficient, γ the nonlinear stiffness coefficient, δ cos ωt the parametric driving and F cos ωt the direct driving term. By setting γ = 0 and ν = 0 one can fit the direct response at low drive level (Fig. 4(a)) which is used to obtain values for μ, β and F. Then, γ is used to fit the large amplitude direct response in Fig. 4(a,c). Finally it is found that the large amplitude parametric resonances in Fig. 4(b,d) can only be fitted using a non-zero value of the nonlinear damping term ν. Numerical values for the fit parameters are provided in the Supplementary Information S1.

Figure 4 Comparison of experimental mechanical responses to theory. (a) Directly driven response at 7.1 and 250.9 mV RMS driving voltage and the fit obtained from Eq. 1. (b) Parametric response and fit at 250.9 mV RMS driving voltage and the fit from Eq. 1. (c) Directly driven response at 446.2 mV RMS driving level, the fit from Eq. 1 shows a disagreement with the backward sweep, highlighted by black arrows. (d) Parametric response at 446.2 mV (RMS). Black arrows highlight the disagreement between Eq. 1 and experiment. (e) Parametric resonance instability map for the fundamental mode of drum 2, compared to the prediction from Eq. 1. (f) Parametric resonance instability map for the fundamental mode of drum 1 (Fig. 2). Full size image