Symmetry-breaking induced frequency combs in graphene resonators

Nonlinearities are inherent to the dynamics of two-dimensional materials. Phenomena like intermodal coupling already arise at amplitudes of only a few nanometers, and a range of unexplored effects still awaits to be harnessed. Here, we demonstrate a route for generating mechanical frequency combs in graphene resonators undergoing symmetry-breaking forces. We use electrostatic force to break the membrane's out-of-plane symmetry and tune its resonance frequency towards a two-to-one internal resonance, thus achieving strong coupling between two of its mechanical modes. When increasing the drive level, we observe splitting of the fundamental resonance peak, followed by the emergence of a frequency comb regime. We attribute the observed physics to a non-symmetric restoring potential, and show that the frequency comb regime is mediated by a Neimark bifurcation of the periodic solution. These results demonstrate that mechanical frequency combs and chaotic dynamics in 2D material resonators can emerge near internal resonances due to symmetry-breaking.

Amongst different nonlinear phenomena that emerge in 2D material membranes, mode coupling is particularly interesting as it allows for the transfer of energy between vibrational states of single [15] or coupled 2D resonators [17]. Mode coupling is also closely linked to nonlinear dissipation [9,18], and can be tuned utilizing internal resonance (IR); a condition at which two or more resonance frequencies become commensurate. The application of IR in mechanical resonators spans from frequency division [19] and time-keeping [20,21] to enhancing the sensitivity of scanning probe techniques [22].
Here, we present a mechanism for generating frequency combs by symmetry-breaking, that exploits internal resonances of a few-nm-thick graphene resonator. We make use of the extreme flexibility of graphene to controllably break its out-of-plane symmetry by bending it using electrostatic force, and achieve two-to-one (2:1) IR between its distinct mechanical modes. Unlike recent demonstrations of frequency comb generation in graphene that require strong coupling of the suspended membrane with a high quality factor SiN x substrate [23], here we show that by careful tuning of the intermodal coupling between two modes of vibration in a single resonator, frequency combs can be generated. As a result of this 2:1 modal interaction, we observe splitting of the resonance peak at a critical gate voltage and drive level, leading to equally spaced spectral lines near the fundamental resonance. By using an analytical model that accounts for the broken symmetry and comprises quadratic coupling terms, we account for the characteristic dependence of the frequency comb region on the membrane tension and deflection amplitude, and confirm that symmetry-broken mechanics lies at the root of the observations.

Experimental Characterization of Frequency Comb
Experiments are performed on a 15 nm thick exfoliated graphene flake, transferred over a circular cavity of 8 µm diameter and 220 nm depth forming a drum resonator. The motion of graphene is read-out in a Fabry-Pérot interferometer where a red helium-neon laser (λ = 633 nm) is used to probe the motion [24,25], (see Figure 1-a, c). The drum is driven opto-thermally using a power modulated blue laser (λ = 485 nm), and to control the static deflection of the drum, a local gate electrode is placed at the bottom of the cavity, see Figure 1-b. Moreover, to reduce damping by the surrounding gas, the sample is measured in a vacuum chamber with pressure ≤ 10 −4 mbar.
By sweeping the modulation frequency f of the blue laser using a Vector Network Analyzer (VNA), we observe multiple directly-driven resonances, appearing as pronounced peaks in the resonator's spectral response (Figure 1-d). Among them, the primary and secondary axisymmetric modes of the drum can be readily identified at f 0,1 = 7.0 MHz and f 0,2 = 15.8 MHz, with f 0,2 /f 0,1 = 2.25, close to the theoretically predicted ratio of 2.29 for a membrane [26]. We note that the resonance frequencies depend strongly on the membrane tension, which we can tune via the electrostatic force generated by the electrostatic gate electrode.
By sweeping the gate voltage V g , we control the tension in the membrane and alter the out of plane offset (see Supplementary Information 1). The electrostatic force pulls the drum out of its initial flat configuration, and breaks its out-of-plane symmetry [27]. This broken-symmetry can have significant influence on the dynamics of the resonator, especially in the nonlinear regime, where the resonant response deviates from the common Duffing model, because it introduces quadratic terms in the nonlinear stiffness [28].
We note that increasing V g causes the resonance frequencies of the drum to shift at different rates (see Supplementary Figure 1). At a certain critical voltage V IR = 7 V we observe (Figure 1-e) splitting of the fundamental resonance peak at f IR =22.73 MHz, which we attribute to the occurrence of a 2:1 internal resonance with a higher mode, since it occurs when the frequency of a higher mode at 44 MHz is exactly twice that of the fundamental mode, see Figure 2-a. Besides splitting, the height of both resonance peaks also diminishes close to V IR , providing evidence for the presence of 2:1 IR and energy redistribution between the interacting modes.
By driving the drum at elevated blue laser powers and performing upward frequency sweeps, we observe in Figure 2-b a butterfly-shaped response, consisting of two Duffing-like asymmetric resonances, one of which bending to lower and the other to higher frequency, indicating that one of the split peaks experiences a spring softening, and the other a spring hardening nonlinearity. Interestingly, at the maximum drive level (10 dBm), the strong coupling between the resonant modes yields the emergence of a third peak in the middle of the split region at frequency f IR = 22.73 MHz (see Figure 2-c). This unexpected additional peak at f IR was not observed in nonlinear resonators undergoing a similar IR [29,30].
In order to investigate this unconventional response in depth, we drive the graphene drum to the critical voltage V IR required to observe the split peak at f IR , and used a Zurich UHFLI to analyze the fast oscillations of the drum at high drive powers. By simultaneously tracing the response spectrum while sweeping the driving frequency around f IR we noticed that for driving frequencies outside the region where the middle peak was spotted, the motion is harmonic. However, close to f IR , the spectral response suddenly changes and a frequency comb is observed consisting of multiple equally spaced peaks near f IR (see

Nonlinear Model of Frequency Comb
To explain the nonlinear physics associated with the observed dynamics and frequency comb near IR in a system with broken-symmetry, we present an analytical model to derive the system's Lagrangian and obtain the governing equations of motion (See Supplementary Information 3). We approximate the coupled motion by only considering the drum's first two axisymmetric modes of vibration with frequencies f 0,1 and f 0,2 . For an ideal circular membrane, the ratio of these first two axisymmetric modes can be tuned to approach f 0,1 /f 0,2 ≈ 2 by changing the tension distribution. These variations in tension distribution might originate from variations in the electrostatic force if the distance to the gate electrode is non-uniform due to membrane deflection, wrinkling or buckling. Moreover, to account for the broken-symmetry mechanics, we model the drum with a static deflection from its undeformed state, that has the shape of its fundamental mode shape [31], with an amplitude W 0 . This leads to the presence of both quadratic and cubic coupling terms in the equations of motion. However, we note that not all the terms in a 2:1 IR are resonant [9], and retain only the relevant terms to obtain the following set of simplified equations near the IR (See Supporting Information 4): Here, x and q are the generalized coordinates, k x and k q are the intrinsic mode stiffness and T x and T q represent added stiffness due to electrostatic tuning of the tension. τ x and τ q are the linear damping coefficients of the generalized coordinates. Moreover, γ is the Duffing coefficient, and α is the coupling strength that can be determined analytically in terms of the offset shape and modes of vibration (See Supplementary Information 3). Finally, F is the forcing amplitude and Ω = 2πf d is the excitation frequency. All the terms in Equation 1 and Equation 2 are mass normalized.
In order to investigate the resonant interaction numerically, we time-integrate the equations of motion. We start by recording the time response of the system at Ω far from resonance and sweep Ω through the 2:1 IR condition. Simulations are performed first at a low driving level that is associated with the linear harmonic oscillator response and then F is increased until the specific characteristics of the nonlinear interaction such as mode splitting appear. We perform our simulations using nonlinear parameters γ = 5.78 × 10 30 (Hz/m) 2 , α = 1.97 × 10 24 (Hz 2 /m). These values correspond to the analytical model of a 15 nm thick drum with a diameter of 8 µm, Young's modulus of E = 1 TPa, and initial axisymmetric offset amplitude of 90 nm. Figure 3-a shows the modelled variation of the resonance frequency as a function of the applied tension (T x ). By changing the tension T x , the fundamental resonance frequency f 0,1 is tuned and a peak splitting, similar to that in Figure 1-e, is observed near the internal resonance frequency f IR .
The splitting phenomenon becomes more apparent at elevated drive powers (see Figure 3-b), similar to the experimental observations in Figure 2-c. This leads to the emergence of a similar butterflyshaped responsivity x/F , as the nonlinear coupling becomes stronger at higher drive levels, where energy leaks to the interacting mode. Interestingly, we also note the presence of the third middle peak in our simulation. In Figure 3-c, it can be seen that this peak indeed appears within the split region, at zero detuning from IR condition, confirming that 2:1 IR that follows from the equations of motion (1) and (2), can be held accountable for our experimental observations.
In Figure 3-c it can be also noted that when driving near f IR the second generalized coordinate q shows an enhanced amplitude, with a response that resembles that of x. It is important to note that in the experiments, the middle peak observed at f IR is only due to the fundamental amplitude x, since our measurements are performed in a homodyne detection scheme.
To better understand the mechanism that lies at the centre of our observation, we investigated the stability of the solution branches using a numerical continuation software package (AUTO). We found that the middle peak appears in a region that is confined between two Neimark bifurcations (red dashed lines in Figure 3-c and Supplementary Information 4). Similar to the Hopf bifurcation, at which a fixed point becomes a limit cycle, at a Neimark bifurcation a periodic orbit becomes a quasi-periodic orbit [32]. Quasi-periodic motion is characterized by a closed invariant curve in the Poincaré map of the phase space that is known to result in amplitude-modulated motion, and thus the emergence of frequency combs in the spectral response [31].
To investigate the spectral characteristics of the quasi-periodic oscillations, we swept the excitation frequency Ω in the spectral neighborhood of the region confined by the two Neimark bifurcations, and analyzed the time response of the nonlinear equations, similar to [33]. Figure 3-d shows the frequency content of the simulated time signal inside and outside this region. It can be observed that the frequencies around f IR are discretely separated from each other, creating a frequency comb that was nonexistent before reaching the onset of Neimark bifurcation, resembling the frequency comb in Figure 2-d. We also show that the time-dependent motion becomes amplitude modulated when entering the Neimark bifurcation region(see Figure Figure 3-e), while having constant amplitude outside of that region. Interestingly, numerical simulations also show signatures of chaotic states upon amplification of the drive level, suggesting that 2:1 IR and broken-symmetry mechanics can represent the onset of a transition from quasi-periodic to chaotic oscillations in 2D material resonators(see Figure Figure 4  By simulating the equations of motion at IR while sweeping the parameters, it is also possible to show that the Neimark bifurcations and thus frequency comb generation is sensitive to mechanical parameters of the system. At 2:1 IR, where the Neimark bifurcation is activated, any change in mechanical properties of the drum will be reflected in the frequency spectrum, as a change in the comb intensity, spacing, and population. Figures 4-b and 4-c reveal the sensitivity of these combs to the drum offset and tension, which were obtained by sweeping the initial offset (broken-symmetry) amplitude and T x . Combs only appear if there is sufficient nonlinear coupling induced by the broken symmetry; increases in the offset influences comb spacing and population. Furthermore, near IR, the frequency comb can be used as a sensitive probe for changes in the parameters of the two interacting modes. Any shift in the resonance frequency of the coupled modes results in changes in comb spacings, making it possible to simultaneously probe changes in both frequencies by solely measuring the response of the fundamental mode at the Neimark bifurcation. External parameters like drive power and drive frequency are also observed to influence frequency comb region, and serve as controls for tuning comb intensity, spacing, and population (See Supplementary Information 5). showing the evolution of phase spaces and sensitivity of frequency comb generation in a graphene drum with broken-symmetry and 2:1 IR to (b) offset amplitude (c) tension variation at the drive levelF = 0.0025. a) Bifurcation diagram of the graphene drum at 2:1 IR, showcasing a quasi-periodic route to chaos. b) Offset amplitude W 0 has been swept while the FFT of the time signal is being extracted in each step. As the offset increases, so does the boundaries of Neimark bifurcation and comb population. c) Added stiffness due to the tension change, T x , has been swept while the FFT of the time signal is being extracted in each step, as the added tension moves the resonance frequencies with respect to the 2:1 IR condition.
In summary, we demonstrate a route for generating frequency combs in the nonlinear response of graphene drums that utilizes broken symmetry and 2:1 internal resonance. Unlike other methods that use multiple wave mixing [34,35], resonant nonlinear friction [36], or SNIC bifurcation [37], to generate mechanical frequency combs, the presented method makes use of an electrostatic gate to controllably tune frequency combs that are mediated by broken-symmetry. When the drum is brought close to the broken-symmetry induced 2:1 IR, we observe strong splitting of the fundamental resonance peak, exhibiting both softening and hardening nonlinearity. Between the split peaks, we observe resonant interactions when driving at relatively high powers, that are generated by Neimark bifurcations of the periodic motion. This regime hosts quasi-periodic oscillations that are held accountable for the observed frequency combs. The experimentally observed phenomena were explained using a continuum mechanics model of a deflected drum with 2:1 IR between its first two axisymmetric modes.
Emerging from the inherent geometric nonlinearities, mechanical frequency combs are closely linked to the mechanical properties of 2D materials, including tension, Young's modulus and brokensymmetry, and thus can be utilized for probing these properties and tracing their variations with frequency and drive levels [23]. The frequency comb generation mechanism described here also provides a platform for demonstrating quasi-periodic route to chaos [32] with nanomechanical resonators and paves the way towards controllable use of IR for sensing applications with 2D materials.

Supplementary Information
Symmetry-breaking induced frequency combs in graphene resonators February 8, 2022 1 Evolution of the overall frequency response with gate voltage In Supplementary Figure 1-a (bottom panel) we show the evolution of the resonance frequencies of our graphene drum as a function of the applied gate voltage, and in top panel we show the frequency spectrum at V g =1.9 V. We observe that change of the gate voltage changes the resonance frequencies at different rates, providing the possibility to obtain IR conditions. ∼ 3 V is a striking feature of this mapping where there is an abrupt jump in the resonance frequencies. The observation coincides with a rapid change in the observed color of the drum (Supplementary Figure 1-b), which points towards a rapid adjustment in the equilibrium position of the oscillator, i.e. snap-through instability (buckling). This is expected, if the membrane cavity is sealed very well in the atmospheric pressure, where it bulges upwards inside the vacuum chamber due to the pressure difference. This means that center of the membrane is above the surface level of the substrate confirming the presence of initial static deflection. The electrostatic attraction between the membrane and the bottom of the cavity, above the threshold voltage, slowly lowers the center of the membrane, until a critical voltage where an instant jump in the resonance frequencies occurs. Supplementary

Additional experimental result
Here, we provide additional experimental result of the same device with different tension and pressure levels, where the 2:1 IR condition was met ( Supplementary Figure 2-a). Similar to the main text, we see a nonlinear splitting of the peak as the drive level is increased. Furthermore, in the center of the splitting, the third peak is observed due to Neimark bifurcaitons. The bifurcation gives birth to quasiperiodic motion, which is measured, observed as an amplitude modulated response (see Supplementary  Figure 2

Equations of motion
To obtain the governing equations of the circular drum, we use the Rayleigh-Ritz approach. The elastic strain energy of the circular drum is obtained as in which h is the thickness of the drum, R is the radius, E is the Young's modulus, and ν is the Poisson's ratio. Moreover, rr , θθ , and γ rθ are the normal and shear strains. During axisymmetric oscillations γ rθ = 0, and the normal strains are obtained in terms of the transverse deflection (w) and radial displacement (u) of the drum as follows In Eq.(2), w 0 is the initial transverse displacement of the drum associated with zero initial stress. Assuming modal interactions between the first two axisymmetric eigenmodes and fixed boundary conditions (u = w = 0), the solution is approximated as where u 0 is the initial radial displacement due to pre-tension n 0 , and x and q are generalized coordinates associated with the first and the second axisymmetric mode of the drum, respectively. Moreover, J 0 is the Bessel function of order zero and β 1 and β 2 are its first two roots. In addition, y j are the radial generalized coordinates and n is the number of these coordinates retained in the approximation. Moreover, we assume that the initial offset w 0 has the same form of the first axisymmetric mode with known amplitude W 0 , thus: The kinetic energy of the drum is then obtained as in which the overdot represents differentiation with respect to time t, and ρ is the mass density. Next, Lagrange equations of motion are obtained [1] leading to the following set of coupled equations 12 xq + α in which α (k) lm and γ k lm are the quadratic and cubic nonlinear terms. It is worth noting that in equations (8) and (9) not all terms are resonant. To recover the resonant terms under 2:1 IR condition (ω 2 2ω 1 ), we assume harmonic motion of the form x ≈ cos(ω 1 t) and q ≈ cos(2ω 1 t) as a first approximation. Inserting these relations in equations 8 and 9 shows that the terms x 3 ≈ 3 4 cos(ω 1 t) + 1 4 cos(ω 1 t) and xq ≈ 1 2 (cos(ω 1 t) + cos(3ω 1 t)) in the first equation of motion are trivially resonant. The same holds for the term x 2 ≈= 1 2 (1 + cos(2ω 1 t)) which can be viewed as a resonant term in equation 9. In a similar fashion it can be shown that the cubic coupling terms xq 2 and q 2 x are dispersive terms. Therefore, the governing equations of motion can reduce tö Eqs. 10 are the normal form of a coupled oscillator undergoing 2:1 IR, and assuming that the second mode does not surpass its Duffing nonlinearity. It is interesting to note that frequencies ω 1 and ω 2 , the coupling term α in which α and γ are evaluated assuming ν = 0.16. We note that the Duffing nonlinearity γ depends on the Young's modulus and Poisson's ratio of the drum [1] while the quadratic coupling α in addition depends on the initial deformation w 0 . In other words, in the absence of W 0 no mechanical coupling exists between the first two eigenmodes of the drum.

Additional simulations
Following on simulations in the main text (Figure 3), here we provide additional simulations. Supplementary Figure 3 shows the frequency response at the internal resonance condition, obtained by the numerical continuation software (AUTO). The Neimark bifurcations are obtained at the points where the motion becomes quasi-periodic in the time integration simulations, also shown in the main text. Furthermore, we investigate the evolution of the phase space as the drive frequency passes the 2:1 IR condition (Supplementary Figure 4). It is possible to observe that, the periodic motion of the resonator turns to quasi-periodic oscillation (which can be also seen from the Poincaré maps) when the bifurcation is triggered (at f d /f IR = 0.9951). The phase space afterwards shows an ergodic behavior till f d /f IR = 1.0043 where the system regains its stable periodic motion through the second Neimark bifurcation. Increasing the drive levels further increases the complexity of the Poincaré maps (Supplementary Figure 5)  5 Sensitivity of the frequency combs to the external drive at 2:1 IR The sensitivity of the comb spacing to the drive frequency can be also studied using our model. In Supplementary Figure 6a it is possible to trace the evolution of the frequency combs as the drive level is increased. Higher drive levels enlarge the Neimark bifurcation region and enrich the spectral response by increasing the comb population, until a certain drive level where the response of the system becomes chaotic. In Supplementary Figure 6b and c we show the variation of comb spacing (∆f ) with respect to the drive frequency in the region mediated by Neimark bifurcations, in experiment (from the measurements in Figure 2-d from the main text) and simulations (from the results in Figure 3-d from the main text). It can be seen that the comb spacing is sensitive to the drive frequency.