Giant Tunable Mechanical Nonlinearity in Graphene-Silicon Nitride Hybrid Resonator

TL;DR

This paper demonstrates a hybrid graphene-silicon nitride resonator system where large, tunable mechanical nonlinearity is achieved, enabling novel phononic frequency combs and advancing control over mechanical resonator nonlinearities.

Contribution

It introduces a hybrid platform combining SiNx and graphene resonators to induce and control giant mechanical nonlinearity, including the observation of a new phononic frequency comb.

Findings

Induced large nonlinear response on SiNx resonator via coupling with graphene.

Achieved tunable nonlinearity controlled by gate voltage.

Observed a novel phononic frequency comb.

Abstract

High quality factor mechanical resonators have shown great promise in developing classical or quantum technologies. Simultaneously, progress has been made in developing controlled mechanical nonlinearity. Here we combine these two directions of progress in a single platform consisting of coupled Silicon Nitride (SiNx) and graphene mechanical resonators. We show that nonlinear response can be induced on a large area SiNx resonator mode and can be efficiently controlled by coupling it to a gate-tunable, freely suspended graphene mode. The induced nonlinear response of the hybrid modes, as measured on the SiNx resonator surface is giant, with one of the highest measured Duffing constants. We observe a novel phononic frequency comb which we use as an alternate validation of the measured values, along with numerical simulations which are in overall agreement with measurements.

Tables

Table 1. Table 1: Parameters for numerical simulation

Parameter	unit	Fig. 2c, main text	Fig. 3c, main text
$m_g$	kg	$10\times 0.625\times {10}^{-16}$	$10\times 0.625\times {10}^{-16}$
$m_s$	kg	$2.38\times {10}^{-11}$	$2.38\times {10}^{-11}$
${\omega }_g$	$s^{-1}$	$2\pi \times 2.864\times {10}^6$	$2\pi \times 3.005\times {10}^6$
${\omega }_s$	$s^{-1}$	$2\pi \times 2.866\times {10}^6$	$2\pi \times 3.007\times {10}^6$
${\gamma }_{bare}^g$	$s^{-1}$	$2\pi \times 25050$	$2\pi \times 8045.8$
${\gamma }_{bare}^s$	$s^{-1}$	$2\pi \times 744$	$2\pi \times 496.4$
$\alpha$	${\mathrm{kgs}}^{-2}$	$4{\pi }^2\times 2.328\times {10}^{-3}$	$4{\pi }^2\times 3.4\times {10}^{-3}$
${\beta }_{\mathrm{bare}}^g$	$N/m^3$	$5.8\times {10}^{13}$	$1.07\times {10}^{12}$
${\eta }_{\mathrm{bare}}^g$	$\mathrm{Ns}/m^3$	$9.8\times {10}^6$	$7.5\times {10}^4$

Equations

$$\begin{gathered} {\ddot{x}}_g=-\frac{{\omega }_g}{Q_g}{\dot{x}}_g-\frac{{\eta }_{\mathrm{bare}}^g}{m_g}{x}_g^2{\dot{x}}_g-\frac{{\beta }_{\mathrm{bare}}^g}{m_g}{x}_g^3 \\ -[{\omega }_g^2+{ϵ}_p\cos ({\omega }_{p}t)]x_g-\frac{\alpha }{m_g}{x}_s, \end{gathered}$$

$${\ddot{x}}_s=-\frac{{\omega }_s}{Q_s}{\dot{x}}_s-{\omega }_s^2{x}_s-\frac{\alpha }{m_s}{x}_g.$$

$${\ddot{x}}_g+{\gamma }_{bare}^g{\dot{x}}_g+{{\omega }_g}^2{x}_g-\frac{\alpha }{m_g}{x}_s=\frac{F_g^{th}}{m_g}$$

$${\ddot{x}}_s+{\gamma }_{bare}^s{\dot{x}}_s+{{\omega }_s}^2{x}_s-\frac{\alpha }{m_s}{x}_g=\frac{F_s^{th}}{m_s},$$

$$S_{x,g}^{1/2}=\Bigg{[}\kappa\bigg{(}\frac{\frac{S_{F,g}^{th}}{m_{g}^{2}}\{(\omega_{s}^{2}-\omega^{2})^{2}+(\gamma^{s}_{bare})^{2}\omega^{2}\}+\{\frac{S_{F,s}^{th}\alpha^{2}}{m_{s}^{2}m_{g}^{2}}\}}{[(\omega_{g}^{2}-\omega^{2})(\omega_{s}^{2}-\omega^{2})-\gamma^{g}_{bare}\gamma^{s}_{bare}\omega^{2}-\frac{\alpha^{2}}{m_{g}m_{s}}]^{2}+[(\omega_{g}^{2}-\omega^{2})\gamma_{s}\omega+(\omega_{s}^{2}-\omega^{2})\gamma^{g}_{bare}\omega]^{2}}\bigg{)}+S_{noise}\Bigg{]}^{1/2},$$

$$F=slope\frac{{\omega }_s^2{m}_s}{Q_s}{V}_{ac}.$$

$$\zeta F=(A{x}_s^2+B{x}_s^4+C{x}_s^6{)}^{1/2}$$

$${\beta }_{hyb}^{\mathrm{s}}=1.54\frac{m_s{\omega }_s^2}{Q_s{x_{hyb}^{s,cr}}^2}$$

$$x_{bare}^{s,cr}=x_{hyb}^{s,cr}\sqrt{\frac{{\beta }_{hyb}^s}{{\beta }_{bare}^s}}$$

$${\ddot{x}}_g=-{\gamma }_{\mathrm{bare}}^g{\dot{x}}_g-\frac{{\eta }_{\mathrm{bare}}^g}{m_g}{x_g}^2{\dot{x}}_g-\frac{{\beta }_{\mathrm{bare}}^g}{m_g}{x_g}^3-[{\omega }_g^2+{ϵ}_p\cos ({\omega }_{p}t)]x_g+\frac{\alpha }{m_g}{x}_s$$

$${\ddot{x}}_s=-{\gamma }_{\mathrm{bare}}^s{\dot{x}}_s-{\omega }_s^2{x}_s+\frac{\alpha }{m_s}{x}_g$$

$$x_g=\frac{1}{2\sqrt{\alpha /m_s}}\left(x_2+x_1\right)=x_2^g+x_1^g,$$

$$x_s=\frac{1}{2\sqrt{\alpha /m_g}}\left(x_2-x_1\right)=x_2^s-x_1^s$$

$${\ddot{x}}_g+{\omega }_0^2{x}_g+\frac{{\beta }_{\mathrm{bare}}^g}{m_g}{x}_g^3-{\alpha }_g{x}_s=0,$$

$${\ddot{x}}_s+{\omega }_0^2{x}_s-{\alpha }_s{x}_g=0.$$

$$x_g^{(0)}=A_g\cos \left[\left({\omega }_0+{\beta }_{\mathrm{bare}}^g\frac{3A_g^2}{8{\omega }_0{m}_g}\right)t\right]-\frac{{\beta }_{\mathrm{bare}}^g{A}_g^3}{32{\omega }_0^2{m}_g}\left(\cos {\omega }_{0}t-\cos 3{\omega }_{0}t\right).$$

$${\ddot{x}}_s+{\omega }_0^2{x}_s-\frac{{\beta }_{\mathrm{hyb}}^s}{m_s}{x}_s^3=B\cos {\omega }_{0}t,$$

$$B\equiv {\alpha }_s{A}_g-\frac{{\alpha }_s{\beta }_{\mathrm{bare}}^g{A}_g^3}{8{\omega }_0^2{m}_g},$$

$${\beta }_{\mathrm{hyb}}^s\equiv \frac{{\alpha }_s{\beta }_{\mathrm{bare}}^g{A}_g^3}{8{\omega }_0^2{m}_g{A}_s^3}.$$

$${\beta }_{\mathrm{hyb}}^s\equiv \frac{\alpha {\beta }_{\mathrm{bare}}^g{m}_s^3{Q}_g^3}{8{\omega }_0^2{m}_g^4{Q}_s^3}$$

$${\beta }_{\mathrm{hyb}}^s\propto \frac{\alpha {\beta }_{\mathrm{bare}}^g{m}_s^3}{m_g^4}.$$

$$x_k=(A_1{e}^{-i\Delta t}+A_1^{\prime }{e}^{i\Delta t}+A_2{e}^{-i3\Delta t}+A_2^{\prime }{e}^{i3\Delta t})e^{i{\omega }_{0}t}$$

$$\begin{gathered} {\beta }_{hyb}^k⟨x_k^2⟩x_k+{\eta }_{hyb}^k⟨x^2⟩\dot{x_k}= \\ {\beta }_{hyb}^k(A_1{e}^{-i\Delta t}+A_1^{\prime }{e}^{i\Delta t}+A_2{e}^{-i3\Delta t}+A_2^{\prime }{e}^{i3\Delta t}{)}^3{e}^{i{\omega }_{0}t} \\ +i{\eta }_{hyb}^k(A_1{e}^{-i\Delta t}+A_1^{\prime }{e}^{i\Delta t}+A_2{e}^{-i3\Delta t}+A_2^{\prime }{e}^{i3\Delta t}{)}^2 \\ \{({\omega }_0-\Delta )A_1{e}^{-i\Delta t}+({\omega }_0+\Delta )A_1^{\prime }{e}^{i\Delta t}\}{e}^{i{\omega }_{0}t} \end{gathered}$$

$$(\beta^{k}_{hyb})^{2}+(\omega_{0}-\Delta)^{2}(\eta^{k}_{hyb})^{2}=36\bigg{(}\frac{A_{2}}{A_{1}^{3}}\bigg{)}^{2}\Delta^{2}\omega_{0}^{2}m_{g}^{2}$$

$$(\beta^{k}_{hyb})^{2}+(\omega_{0}+\Delta)^{2}(\eta^{k}_{hyb})^{2}=36\bigg{(}\frac{A^{\prime}_{2}}{{A^{\prime}}_{1}^{3}}\bigg{)}^{2}\Delta^{2}\omega_{0}^{2}m_{g}^{2}$$

$$9(\beta^{k}_{hyb})^{2}+(3\omega_{0}-5\Delta)^{2}(\eta^{k}_{hyb})^{2}=100\bigg{\{}\frac{A_{3}^{2}}{(A_{1}^{2}A_{2}+A_{1}^{{}^{\prime}}A_{2}^{2})^{2}}\bigg{\}}\Delta^{2}\omega_{0}^{2}m_{g}^{2}$$

$$9(\beta^{k}_{hyb})^{2}+(3\omega_{0}+5\Delta)^{2}(\eta^{k}_{hyb})^{2}=100\bigg{\{}\frac{{A^{\prime}}_{3}^{2}}{({A^{\prime}}_{1}^{2}{A^{\prime}}_{2}+A_{1}{A^{\prime}}_{2}^{2})^{2}}\bigg{\}}\Delta^{2}\omega_{0}^{2}m_{g}^{2}$$

$$\beta^{k}_{hyb2}=\frac{\omega_{0}\Delta m_{g}}{\sqrt{12\omega_{0}\Delta}}\bigg{[}100(\omega_{0}^{2}-2\omega_{0}\Delta)\frac{A_{3}^{2}}{(A_{1}^{2}A_{2}+{A^{\prime}}_{1}A_{2}^{2})^{2}}\\ -36(9\omega_{0}^{2}-30\omega_{0}\Delta)\bigg{(}\frac{A_{2}}{A_{1}^{3}}\bigg{)}^{2}\bigg{]}^{1/2}$$

$$\eta^{k}_{hyb2}=\frac{\omega_{0}\Delta m_{g}}{\sqrt{12\omega_{0}\Delta}}\bigg{[}-100\frac{A_{3}^{2}}{(A_{1}^{2}A_{2}+{A^{\prime}}_{1}A_{2}^{2})^{2}}\\ +324\bigg{(}\frac{A_{2}}{A_{1}^{3}}\bigg{)}^{2}\bigg{]}^{1/2}$$

$$\beta^{k}_{hyb1}=\frac{\omega_{0}\Delta m_{g}}{\sqrt{12\omega_{0}\Delta}}\bigg{[}-100(\omega_{0}^{2}+2\omega_{0}\Delta)\frac{A_{3}^{2}}{(A_{1}^{2}A_{2}+{A^{\prime}}_{1}A_{2}^{2})^{2}}\\ +36(9\omega_{0}^{2}+30\omega_{0}\Delta)\bigg{(}\frac{A_{2}}{A_{1}^{3}}\bigg{)}^{2}\bigg{]}^{1/2}$$

$$\eta^{k}_{hyb1}=\frac{\omega_{0}\Delta m_{g}}{\sqrt{12\omega_{0}\Delta}}\bigg{[}100\frac{A_{3}^{2}}{(A_{1}^{2}A_{2}+{A^{\prime}}_{1}A_{2}^{2})^{2}}\\ -324\bigg{(}\frac{A_{2}}{A_{1}^{3}}\bigg{)}^{2}\bigg{]}^{1/2}$$

Full text

††preprint: XXX

Introduction

For more than a century, mechanical resonators Braginsky92 have played a central role in measuring forces Cavendish-1798 ; Abbott-16 ; Weber-16 and testing fundamental physical principles Clerk10 ; Wollman-15 ; Ockeloen-Korppi-18 ; Marinkovic-18 . With the advent of micro and nanoscale mechanical resonators, and in particular, after experimental observation of their quantum mechanical behavior OConnell-10 ; Taufel-11 ; Schliesser-11 ; Aspelmeyer-14 , there has been a renewed interest in usage of such resonator modes in classical Mahboob-11 and quantum technologies Rips-16 ; Mika-07 . Significant progress has been made in two broad directions over the last decade. On one hand, there has been progress in developing high quality factor ($Q$) mechanical resonators modes at high resonant frequency ($f_{0}$) towards quantum devices at room temperature  kippenberg-07 ; chakram-14 . Silicon Nitride (SiNx) has emerged as a dominant material of choice Fink-16 for such resonators, demonstrating mechanical $Q$’s in excess $10^{7}$ at MHz frequencies Reinhardt-16 ; Norte-16 ; Patil-15 . On the other hand, progress has been made in developing tunable mechanical nonlinearity, towards conditional phase shifts of mechanical modes Mahboob-11 . Freely suspended graphene resonator, with low mass and high Young’s modulus Lee-08 ; Eichler-11 ; Davidovikj-17 ; Storch-18 resulting in exceptional nonlinear elastic properties Eichler-11 ; Davidovikj-17 along with gate tunable resonant frequency, emerged as an efficient choice; mixing Karabalin-09 ; Mahboob-12 ; Cao-14 ; Mahboob-16 ; Ganesan-17 ; Seitner-17 ; Ganesan-18 ; Czaplewski-18 and side-band cooling of its strongly coupled modes have been observed Alba-16 ; Mathew-16 . A platform that can integrate these two directions of progress, combining high quality factors of SiNx resonators with gate tunable response of graphene resonators can be a logical next step Fink-16 .

Here we explore a hybrid platform consisting of a large area SiNx resonator coupled to an atomically thin, freely suspended graphene that is deposited on holes etched on the SiNx. When a mechanical mode of the graphene resonator is electrostatically tuned into resonance with a SiNx mechanical mode, we observe the resulting hybrid modes develop giant nonlinear response to an external driving force, as measured on the surface of the SiNx resonator (fig. 1a). To validate the measurements, we develop an alternate, novel methodology to characterize third order (Duffing) nonlinear response and damping coefficients of these hybrid modes. By parametrically driving the coupled modes, we observe generation of novel frequency comb and use the measured amplitudes of the generated comb lines to estimate the nonlinearities and validate the results. Our measurements are in agreement with numerical simulations of a model of coupled linear and nonlinear oscillators. The model suggests that induced nonlinearity of hybrid modes, as measured on SiNx, is due to back-action force of the nonlinear graphene resonator and simple scaling estimates are in agreement with the measured giant values. These result, verified with two separate measurements, thereby combine two directions of development of electromechanical resonators in a single hybrid device – a gate tunable graphene resonator inducing giant nonlinearity in a high-$Q$ mechanical mode of a large area SiNx resonator (fig. 1a).

Nonlinear response of Graphene - Silicon Nitride Hybrid

The device consists of a 300 nm thick SiNx resonator of dimensions 320 $\times$ 320 $\mu$m2, with through holes of diameters $10,\,15,$ and $20$ $\mu$m etched on to it (fig. 1b) and monolayer (CVD) graphene is deposited on the holes Nicholl-15 . Both graphene and SiNx resonators are actuated electrostatically with a highly doped silicon back gate, separated by an insulator which results in a net separation of 10 $\mu$m between graphene/SiNx and back gate. A fiber based confocal microscope, as part of a path stabilized Michelson interferometer, is used to detect optical signals reflected from the sample (3d cartoon of fig. 1a).

When the microscope is focused on a 20 $\mu$m freely standing graphene membrane, we observe thermally driven modes, corresponding to that of a circular resonator with anisotropic tension Singh-18 ; Davidovikj-16 . The modes are tunable in excess of 1 MHz with a d.c. gate voltage (fig. 1c) and as we tune from 0 V to 250 V, we see distinct avoided level crossing that signals hybridization with modes of a SiNx resonator (fig. 1c and right panel). In a recent work we have found that a bilinear coupling model fits well with observed hybridized Brownian spectra. From the fitting, we extract the quality factors: $Q_{g}\sim 254$ for graphene and $Q_{s}\sim 3800$ for SiNx Singh-18 . The effective mass of the graphene mode, estimated from dispersion of fig. 1c is estimated at $m_{g}\sim 10m_{0}$, where $m_{0}$ is the mass of a single layer of graphene Singh-18 . We estimate the effective mass, $m_{s}$, of SiNx from its dimensions and density, and find it to be $\sim 10^{4}m_{g}$ Singh-18 .

Comparatively heavier mass of SiNx results in smaller amplitude for the Brownian power spectrum, which is below our detection sensitivity (see supporting information (SI)). However, when the device is actuated with a.c. gate voltage and the microscope is focused on SiNx surface (away from graphene), we observe dense distribution of SiNx resonator modes (fig. 1d). The mode-densities and their dispersion (fig. 1d, right panel) match well with simulated modes of a square membrane of comparable dimensions with an inbuilt tension of 80 MPa (see SI).

When SiNx modes are not hybridized with graphene, the peak amplitudes increase linearly with applied a.c. gate voltage, up to a maximum amplitude of 20 V that we can apply in our experiment. In particular, we measure three modes of SiNx at frequencies 3.071 MHz (mode 1), 3.092 MHz (mode 2), and 3.146 MHz (mode 3) as shown in fig. 1e. When the low-$Q$ fundamental mode of graphene is tuned into resonance by applying a d.c. gate voltage $V_{g}=210$ V, we observe frequency shifts of the three modes and a increased linear response due to back-action of the coupled graphene mode (see S.I.). Beyond a certain gate voltage, the response becomes nonlinear (fig. 1f). The data fits well with the model of ref. Davidovikj-17 , for the steady state amplitude ($x_{s}$) of a forced oscillator with an additional nonlinear response that is cubic (Duffing) in $x_{s}$. From fitting, we extract the effective Duffing constant, $\beta^{s}_{\rm hyb}$, for example, for the hybrid mode $2^{\prime}$ to be $\beta^{s}_{{\rm hyb2}}=8.0(\pm 0.8)\times 10^{21}$ $\rm N/m^{3}$ (see SI). This is one of the highest measured Duffing constants Huang-16 ; Gajo-20 , orders of magnitude larger than graphene’s reported values Eichler-11 ; Davidovikj-17 .

What leads to such giant Duffing constants for the hybridized SiNx modes? Aforementioned good fit of a Duffing model to the SiNx hybrid mode implies the source to be coupling graphene, a highly nonlinear Duffing oscillator. It is therefore critical to characterize the Duffing constant of the graphene resonator. For the graphene mode, we observe distinct signatures of Duffing-like hysteresis in response and asymmetric broadening due to nonlinear damping, when driven on resonance (see SI). Recent studies have characterized Duffing constant and nonlinear damping of driven graphene resonators by accurate fitting of data to theory Eichler-11 ; Davidovikj-17 . However, in our device, hybridization of graphene with multiple SiNx modes results in a complex spectra (see SI). Instead, we develop an alternative methodology to estimate the nonlinear parameters.

Mechanical frequency comb on Graphene surface

The spectra simplifies, when we drive the system parametrically at a frequency that is sum of two dominant hybrid modes (fig. 2a) Singh-18 . Parametric drive leads to gain in only two specific hybrid modes that are phase matched as opposed to the scenario of direct driving where multiple modes may interact. From the corresponding spectra, we develop a methodology to estimate Duffing constants of bare graphene as well as that of the hybrid modes as measured on graphene and on SiNx resonator surfaces.

To parametrically drive the system, we first tune the fundamental mode of graphene at a frequency, $\omega_{0}$ = 2.865 MHz, at which it strongly hybridizes with a SiNx mode. With the microscope focused on the graphene, hybridization of modes is distinctly visible in the form of a splitting into two modes at frequencies $\omega_{1}$ and $\omega_{2}$ (say). We simultaneously apply an a.c. gate voltage (parametric pump) at exactly twice the resonant frequency, $\omega_{0}$ (fig. 2a). As the amplitude of the parametric pump voltage is increased, up to a threshold voltage $V_{c}=$ 11.9 V we observe gain in both the hybridized modes (see SI and ref.  Singh-18 ). Above threshold, new frequency components on either side of the two hybridized modes develop. The number of such modes increases with increasing pump voltage, eventually spanning out into a “comb” like pattern (fig. 2b) Ganesan-17 ; Ganesan-18 . We next develop a theoretical model to account for the new generated comb lines, towards estimating nonlinearities in the system.

Theoretical Model

A model of graphene as a 1d oscillator with a quality factor $Q_{g}$, a third-order nonlinear response described by an effective Duffing constant ($\beta^{g}_{\rm bare}$) along with nonlinear damping ($\eta^{g}_{\rm bare}$)  Eichler-11 ; Lifshitz-09 ; Nayfeh-07 , coupled to a linear SiNx resonator mode explains the observations well. Specifically, we simulate the following set of equations:

[TABLE]

and

[TABLE]

Here $\alpha$ is an effective coupling constant, $x_{g,s}$ are the amplitudes of vertical displacements of graphene ($g$) and SiNx ($s$) resonators modes respectively, and $\epsilon_{p}$ denotes the magnitude of the parametric drive. $\beta^{g}_{\rm bare}$ is the Duffing constant of the bare graphene resonator mode while $\eta^{g}_{\rm bare}$ is the coefficient of nonlinear damping.

Numerically simulated spectra is in agreement with observations (see fig. 2). In particular, we find linear coupling explains the frequency comb, as opposed to nonlinear couplings in earlier works Karabalin-09 ; Mahboob-12 ; Cao-14 ; Mahboob-16 ; Ganesan-17 ; Seitner-17 ; Ganesan-18 ; Czaplewski-18 . By fitting simulated spectra to measured spectra, we get an estimate for the Duffing constant of the bare (non-hybridized) graphene mode Storch-18 to be $\beta^{g}_{\rm bare}=5.8\times 10^{13}$ $\rm N/m^{3}$ with a nonlinear damping coefficient $\eta^{g}_{\rm bare}=9.7\times 10^{6}$ $\rm Ns/m^{3}$, in close agreement with recent measurements Davidovikj-17 ; Storch-18 .

The model further suggest that the essential physical mechanism behind the comb can be understood with the normal (hybrid) modes, even in the strong driving regime (fig. 3). In particular, it is well understood that response of a parametrically driven mode becomes unstable above a threshold: beyond threshold, the instability region extends to form a tongue shaped region Hsu-63 ; Landau-82 ; Rand-12 ; Kovacic-18 . The envelope of the tongue is set by pump amplitude, nonlinear frequency, and damping. Therefore, for two hybridized modes with frequencies, $\omega_{1}$ and $\omega_{2}$, there should be two such independent instability tongues (fig. 3a). Consequently, there ought to be a region of overlap (dark region II, fig. 3a) Hansen-85 . While region I and region III correspond to self-oscillation of hybrid modes 1 and 2 respectively, in the overlap region II the system is multi-periodic. Moreover, in the instability region, large amplitude leads to strong nonlinear response. One therefore expects mixing of two accessible frequencies ($\omega_{1}$ and $\omega_{2}$) in the overlap region II. At a specific parametric drive amplitude (dotted line in fig. 3a), one thereby expects to observe these three regions.

We indeed observe these three regions when we vary the pump frequency across $\omega_{p}=2\omega_{0}$ at a fixed drive amplitude, $\rm V_{p}$ (dotted line in fig. 3a and also fig. 2a). In particular, the frequency is scanned over a range of 20 kHz around $2\omega_{0}$, keeping its amplitude fixed at $V_{p}$= 20 V (fig. 3b). We observe the single frequency self-oscillation regions I and III, on either side of region II that is characterized by the frequency comb (fig. 3b). Corresponding experimental observations match well with simulations (fig. 3c).

Estimating Graphene Nonlinearity

Fig. 3a indicates that the right and the left boundaries of region II correspond to the instability tongues of the hybridized modes, viz., mode 1 and mode 2, respectively. One can then ascribe the observed asymmetry of the envelope of fig. 2b to differing effective nonlinearities of the two modes. Furthermore, from the experimentally measured amplitudes ($A_{i}$ and $A^{\prime}_{i}$, $i$=1,2,$\cdots$) of the new comb lines generated due to cubic nonlinearity and nonlinear damping, we estimate the average nonlinear coefficients as $\beta^{g}_{\rm hyb1}$= $1.3(\pm 0.4)\times 10^{15}$ $\rm N/m^{3}$, ${\beta^{g}_{\rm hyb2}}$ = $8.6(\pm 4.6)\times 10^{14}$ $\rm N/m^{3}$, ${\eta^{g}_{\rm hyb1}}$ = $7.0(\pm 2.2)\times 10^{7}$ $\rm Ns/m^{3}$, and ${\eta^{g}_{\rm hyb2}}$ = $4.8(\pm 2.5)\times 10^{7}$ $\rm Ns/m^{3}$, for the two hybrid modes 1 and 2, as measured on the graphene surface (see SI). The estimated parameters from data are in good agreement with numerical simulations (see SI), substantiating our methodology.

There is a phase relationship of the generated modes with respect to the fundamental modes, at frequencies $\omega_{1,2}$ , and in principle, the nonlinear coefficients can also be estimated by carefully measuring relative phases of the generated modes. For the spectrum of fig. 3d, we observe pulses in time domain (fig. 3e). Repetition rate of the pulses correspond to inverse of 2$\Delta=\omega_{1}-\omega_{2}$, pulse width to inverse of the envelope of the generated comb while the carrier frequency to inverse of the carrier frequency $\omega_{0}=2.760$ MHz (fig. 3d). The frequency comb of the hybrid mode is therefore phase coherent and Fourier transform limited.

Induced frequency comb on Silicon Nitride surface

It can be noted that due to widely varying masses and quality factors of the two physical resonators, values of Duffing constants of a hybrid mode would differ when measured on graphene or on SiNx resonator surface (see SI). Interestingly, signature of the comb spectrum is also observable, when the microscope is focused on the surface of the SiNx resonator (fig. 4). However, the amplitude of oscillations is orders of magnitude smaller than that on graphene, due to significantly heavier mass of SiNx. Accordingly, signatures of measured spectrum are less pronounced. Nevertheless, we use the developed methodology to estimate Duffing constants of the hybrid modes on SiNx. In fig. 4, we focus on a SiNx mode at 2.970 MHz, while applying a parametric drive at twice its resonant frequency. When the fundamental graphene mode is off-resonant, we do not observe any parametric gain (blue region on the left of fig. 4a and black dots in the right inset). However, when the graphene mode is tuned across resonance with a gate voltage between $193\,{\rm V}$ (approx.) to $200\,{\rm V}$, we observe generation of frequency comb as well as single frequency self-oscillation regime. Furthermore, the induced nonlinearity of the hybrid mode extends over the entire SiNx surface and we observe generation of combs at distances in excess of 200 $\mu$m from the edge of the graphene drum that is only 20 $\mu$m in diameter (fig. 4b, bottom panel). Essentially, the localized mode of the graphene acts as a defect center, on the large area oscillating mode of SiNx. From amplitudes of the generated modes, we estimate the effective Duffing constant and nonlinear damping of the hybrid modes to be ${\beta^{s}_{\rm hyb1}}$= $3.4(\pm 0.1)\times 10^{23}$ $\rm N/m^{3}$, ${\beta^{s}_{\rm hyb2}}$= $6.3(\pm 5.3)\times 10^{22}$ $\rm N/m^{3}$ and ${\eta^{s}_{\rm hyb1}}$= $1.8(\pm 0.1)\times 10^{16}$ $\rm Ns/m^{3}$, ${\eta^{s}_{\rm hyb2}}$= $3.3(\pm 2.8)\times 10^{15}$ $\rm Ns/m^{3}$ respectively, as measured on SiNx surface (see SI). The estimates are in close agreement with effective Duffing constant ${\beta^{s}_{{\rm hyb}}}$ measured on SiNx in fig. 1f.

Discussion

To conclude, here we have explored nonlinear response of graphene-SiNx hybrid modes and developed a methodology to quantify corresponding nonlinear coefficients, as measured on graphene and SiNx resonator surfaces. The observations suggest that the coupled system can be described by two uncoupled, nonlinear hybrid modes. Measured Duffing constants of these hybrid modes on SiNx surface are found to be in excess of eight order of magnitude larger than that on graphene. This indicates that nonlinear response is highly efficient on SiNx surface, setting in at displacement scale that is two orders of magnitude smaller, at 30.4 pm, compared to measurement on graphene.

It is remarkable that an atomically thin resonator generates a significant backaction force ($F^{\rm ba}=\alpha x_{g}$) on SiNx. Based on the $F^{\rm ba}$, a perturbative estimate yields $\beta^{s}_{\rm hyb}\propto\alpha\beta^{g}_{\rm bare}m_{s}^{3}/m_{g}^{4}$ (see SI) and indicates graphene to be a powerful candidate to induce such giant nonlinearity due to three primary factors: firstly, pristine graphene robustly couples to SiNx substrate via stable electrostatic forces resulting in a large coupling strength ($\alpha$). This also leads to better device yield. Secondly, low mass of graphene ($m_{g}$) results in a large amplitude ($x_{g}$) of oscillation, boosting the force further. Finally, exceptionally large Young’s modulus results in large nonlinear response ($\beta^{g}_{\rm bare}$) Davidovikj-17 to an applied force. Large gate tunable backaction force of graphene thereby emerges as the dominant mechanism behind observations in this work.

For our device, the tension of SiNx resonator is merely 80 MPa Fink-16 , leading to comparatively lower quality factors ($\sim 3000$ on average), along with a dense distribution of SiNx modes (fig. 1d). An immediate improvement can therefore be towards increasing inbuilt tension of the SiNx resonator, so that one can resolve mode shapes distinctly Yang-19 and observe graphene induced interaction between SiNx modes of quality factors in excess of $10^{6}$, possibly in a quantum regime. With such improvements, the hybrid device proposed here can provide a powerful platform for generating mechanical squeezed states in precision measurements and controlled interactions of mechanical modes both in classical and quantum domains at room temperature Reinhardt-16 ; Norte-16 ; Patil-15 .

Materials and Methods

Sample preparation:

Silicon Nitride membranes (thickness 300 nm) are fabricated by depositing low-stress silicon-rich silicon nitride on both sides of a silicon chip. An array of holes of 10, 15 and 20 um diameter is then patterned in the nitride using standard fabrication procedures. A metallic contact (20 nm Au) is deposited onto the top surface of the SiNx to facilitate electrical gating. Monolayer chemical vapor deposition (CVD) graphene with flake size of $\sim 90\times 110\leavevmode\nobreak\ \rm\mu m^{2}$ is then transferred onto holes in the nitride membranes. We use a high-quality atmospheric CVD growth and wet transfer. The samples are subsequently annealed in an $\rm Ar-H_{2}$ environment at $350^{\circ}C$. The graphene membranes remained clamped to the sample chip via Van der Waals interactions forming suspended circular graphene membranes.

Experimental setup:

We use a fiber based confocal microscope (see SI) with a spot size of 4 $\rm\mu m$ to optically probe our graphene-SiNx hybrid device. The microscope forms one arm of a Michelson interferometer while the reference arm is actively stabilized against drifts or fluctuations through a feedback form PI lock. A frequency and amplitude stabilized external cavity diode laser (ECDL) ($\rm\lambda$ = 780 nm) is used as an optical probe. All the measurements were conducted at probe power of $\sim$400 $\mu$W. The sample is placed inside a vacuum chamber ($10^{-2}$ $\rm mbar$) with high voltage electrical leads for gate control. The entire chamber assembly is mounted on a 3D scanning stage with active position locking. For detection, we use a balanced photo-detector with a detection bandwidth of 45 MHz. We position the sample by actively monitoring the generated 2-D confocal image, which helps in selecting the relative probe position and to lock the microscope there. The photo-current signal is analyzed with spectrum analyzer and dual-lock-in-amplifier.

Acknowledgements

We thank Srivatsan Chakram, Deb Shankar Ray, Edgar Knobloch, Siddharth Tallur, Mandar Deshmukh and Amit Agarwal for insightful discussions and comments. We also thank Om Prakash for his numerous help in construction of the experimental setup. A.S. acknowledges CSIR and K.B. acknowledges ERC grant no. 639739 and DFG TRR 227 for financial support. This work was supported under DST grant no. SERB/PHY/2015404.

I S.1. Experimental method

Experimental Setup: Fig. S.1 below illustrates the experimental setup. Fig. 1c, 2b, 3b, 3d, 4a and 4b in the main text and S.2-4, S.12 and S.15 in the supplement are taken with spectrum analyzer while Fig. 1d-f in the main text and S.5-6 and S.7 in the supplement are acquired by scanning the drive frequency from a lock-in-amplifier.

II S.2. Calibration

II.1 A. Displacement calibration from hybrid Brownian spectrum

Displacements are calibrated by fitting the Brownian spectrum to a model that is based on coupled modes of graphene and SiNx resonators, denoted by displacements $x_{g}$ and $x_{s}$, respectively. For thermally driven graphene and SiNx modes we ignore the nonlinear terms. The equations of motion are then given by:

[TABLE]

[TABLE]

[TABLE]

where $\gamma^{k}_{\rm bare}$, $\omega_{k}$ and $F^{th}_{k}$ $(k=g,s)$ represent linear damping, normal mode frequency and thermal forces acting on graphene and $\rm SiNx$ modes respectively. Coupling of graphene and SiNx modes is modeled by an effective interaction Hamiltonian, $H_{\rm int}=\alpha x_{g}x_{s}$111The bi-linear form of interaction is an approximate expression, derived from an effective Hamiltonian that can be expressed as: $H_{\rm int}=-\frac{\alpha}{2}(x_{g}-x_{s})^{2}$, for small displacements of SiNx and with renormalized resonant frequencies.. Solving the above coupled equations in Fourier space, the displacement power spectrum for the graphene resonator is:

[TABLE]

where $S_{F,k}=4k_{B}T\gamma^{k}_{bare},(k=g,s)$ is the thermal force acting on graphene ($g$) and SiNx ($s$). The calibration factor, $\kappa$ along with all other free parameters are extracted by fitting experimental data, $S_{v,g}^{1/2}$ $\rm(V/\sqrt{Hz})$ to the above equation. Using the calibration factor, the recorded spectrum is then converted into displacement spectrum.

Extracted values of the fitting parameters for the Brownian mode corresponding to Fig. S.2 are listed below.

$\kappa=1.921(\pm.074)\times 10^{10}$ $\rm V^{2}/m^{2}$, $\omega_{g}/2\pi=2.8646(\pm.0002)\times 10^{6}$ Hz, $\omega_{s}/2\pi=2.8656(\pm.0001)\times 10^{6}$ Hz, $\gamma^{g}_{bare}/2\pi$=11.237($\pm$ 0.690)$\times 10^{3}$ Hz, $\gamma^{s}_{bare}/2\pi=0.744(\pm 0.234)\times 10^{3}$ Hz, $\alpha/4\pi^{2}=1.978(\pm 0.095)\times 10^{-3}$ $\rm kgHz^{2}$.

This model can be extended for graphene interaction with multiple SiNx modes. The supplemental information of Singh, R. et.al. Singh18 can be referred for more information.

We calibrate the amplitude of graphene mode by fitting its thermal or Brownian spectrum to a model of coupled 1d oscillators as shown above. The calibrated thermal mode is then used to calibrate the displacement spectrum of the parametrically driven mode. Similarly the thermal mode of SiNx is calibrated, either by fitting it to its Brownian spectrum, or using the calibration from graphene, where we carefully maintain all other microscope parameters.

II.2 B. Mass estimation

In order to estimate the mass of graphene, we fitted the fundamental mode dispersion of the graphene drum using continuum mechanics model Chen-13 . We obtained a good fit for $m_{g}=9.982(\pm 0.008)$ $\times$ $m_{0}$ (where $m_{0}=0.625\times 10^{-16}$ kg is the mass of pristine residue-free graphene resonator) and tension is $6.91(\pm 0.01)\times 10^{-5}$ N/m (Fig. S.3).

II.3 C. Calibration with probe power

We studied dependence of thermal motion of hybrid mode of graphene fundamental mode with increasing incident probe power. The spectrum at different probe powers were recorded and fitted with equation S2. The extracted calibration factor, $\kappa$ depicts linear scaling with the probe power. Invoking the equipartition theorem, we estimated the mode temperature and have not observed significant variation. Overall, these measurements suggest a range up to 1 mW of probe power over which the graphene mode remains unperturbed. For all measurements reported in this work, we maintained the probe power at $\sim 400$ $\rm\mu W$.

III S.3. Graphene and SiNx resonator modes

III.1 A. Nonlinear modes of Graphene resonator

The sample consists of a large area SiNx resonator ($320\times 320\times 0.3$ $\rm\mu m^{3}$), with 20 $\mu$m, 15 $\mu$m and 10 $\rm\mu$m diameter circular holes etched on to it over which monolayer graphene is deposited, thereby forming suspended drums. The graphene gets clamped to SiNx at the edges via Van der Waals forces while rest of the part above the hole remains freely hanging. When the microscope is focused on graphene, we observe its thermo-mechanical spectrum from the detected signal in a electronic spectrum analyser.

The dispersion of fundamental mode of graphene as a function dc gate voltage (Fig. 1c) represents deviation of mode shape from usual Lorentzian shape, due to interaction of the graphene vibrational mode with densely packed, multiple SiNx modes. Such interactions lead to sharp dips and Fano-like asymmetry in graphene spectrum. The asymmetry gets more pronounced when the graphene mode is driven on resonance (Fig. S.5a). The corresponding phase profile shows an overall envelop corresponding to $\pi$-phase jump, as one crosses the broad graphene resonance. However, finer features in phase profile with sharp, intermediate $\pi$ phase jumps correspond to individual SiNx modes which are coupled to graphene mode. It can be noted that in general, each of these narrow SiNx modes have a unique coupling strength to the same graphene.

When driven harder, the broad graphene mode shows an asymmetric signature in spectra that is typical of a Duffing oscillator with cubic nonlinear response in displacement (Fig. S.5b). Forward and backward sweeps of drive frequency shows hysteresis in both the amplitude and phase.

Fig. S.6a shows transition from a linear response to a Duffing-like nonlinear response for the graphene mode, coupled to multiple SiNx modes. The spectrum shows an increase in the FWHM (full width at half maximum) with increasing drive voltage, pointing towards the existence of nonlinear damping. The amplitude of the modes with drive voltage show saturation after a certain drive voltage and the critical voltage for saturation is specific to a given mode (Fig. S.6b).

III.2 B. Modes of large area Silicon Nitride resonator

Silicon Nitride is a large area ($320\times 320\times 0.3$ $\rm\mu m^{3}$) resonator with through holes. Fig. S.7 shows the amplitude and phase of weakly driven SiNx modes. The modes are densely packed with quality factor in the range of 1000-4000. From COMSOL simulation, we estimate the inbuilt tension, $T_{s}$ $\sim$ 80 MPa.

III.3 C. Linear response of the hybrid Silicon Nitride mode

In Fig. 1e,f in main text, one can see a significant difference between uncoupled and hybridized SiNx mode response even in the linear response regime. This linear response behavior can be appreciated by noting (Fig. S8) that SiNx when coupled to graphene experiences two different forces: a direct capacitive force and a back action force from graphene. Even in the linear regime, the backaction force experienced by SiNx from graphene is larger than the capacitive force under direct a.c. drive. This is the reason behind the larger displacement or steeper slope in linear regime of SiNx when on resonant with the graphene. To validate our assertion, we have conducted further experiments, comparing the two differing linear regimes. In particular, we consider the case when graphene is off resonant to the SiNx mode. We record the peak amplitude of vibration for SiNx and graphene modes independently, corresponding to the a.c. voltage. The Fig. S8b and c show the response of SiNx and graphene resonator modes, with the a.c. drive voltage, respectively. From fitting, we find the slope of graphene response (in linear region) is  470 times larger than that of SiNx i.e, $x_{g}^{c}\sim 470\times x_{s}^{c}$.

When the SiNx mode is on resonant to the graphene mode, an additional coupling dependent backaction force acts on the SiNx which can be written as follow: $F_{ba}^{g\rightarrow s}=\alpha(x_{g}^{c}-x_{s}^{c})\sim\alpha_{g}^{c},x_{s}^{c}<x_{g}^{c}$. The displacement experienced by the SiNx due to this backaction force is: $x_{s}^{ba}=F_{ba}^{g\rightarrow s}\frac{Q_{s}}{m_{s}\omega_{0}^{2}}=\alpha x_{g}^{c}\frac{Q_{s}}{m_{s}\omega_{0}^{2}}$.

However, we have $x_{g}^{c}=470\times x_{s}^{c}$. Plugging this relation and typical values of other parameters ($Q_{s}=3000,\omega_{0}=2.86\leavevmode\nobreak\ \rm MHz,m_{s}=2.32\times 10^{-11}\leavevmode\nobreak\ kg,\alpha=1.19\times 10^{-3}\leavevmode\nobreak\ \rm kgHz^{2})$ in above equation, yields $\frac{x_{s}^{ba}}{x_{s}^{c}}=14$. The result suggests that the displacement of SiNx mode due to backaction force is 14 times larger than that of capacitive force. This is consistent with the observation of Fig. 1e and f, main text, where a larger linear response is observed on SiNx, when on resonant with graphene. It can be noted that such displacements are still smaller than the nonlinear threshold displacement of SiNx. However, when the graphene is driven into nonlinear regime, the response of the hybrid mode, as measured on SiNx also becomes nonlinear.

III.4 D. Nonlinear response of hybrid Silicon Nitride modes under direct driving

We observe SiNx modes, respond linearly to external drive when it is off resonant from graphene mode. However, when we couple a graphene mode to the SiNx mode, its response becomes nonlinear with applied forces and shows saturation in amplitude above a critical displacement.

To quantify nonlinear response of the SiNx modes, we follow the procedure described in D. Davidovikj et. al. David-17 . We first extract the slope from the linear region of $x_{s}$ vs $V_{ac}$ plot (Fig. 1e, main text). The rescaled force, $F$ corresponding to the $V_{ac}$ is given by,

[TABLE]

This rescaled force is plotted with $x_{s}$, the steady-state response of a Duffing oscillator, such that David-17 :

[TABLE]

where $C=\frac{9}{16}(\beta^{s}_{hyb})^{2}$. Here $\zeta$ depends on the geometry of the mode and is of the order of 1.

The value of $\beta^{s}_{hyb1}$ for hybrid mode $1^{\prime}$ extracted from fitting result (Fig. S.9) is $5.6(\pm 0.5)\times 10^{21}$ $\rm N/m^{3}$. Similarly for hybrid mode $2^{\prime}$, $\beta^{s}_{hyb2}=8.0(\pm 0.8)\times 10^{21}$ $\rm N/m^{3}$ and for hybrid mode $4^{\prime}$, $\beta^{s}_{hyb4}=7.6(\pm 0.2)\times 10^{21}$ $\rm N/m^{3}$ respectively. One can also estimate the nonlinear coefficient of the hybrid SiNx modes using critical displacement of that mode Zwickl-08 . It is given by:

[TABLE]

For hybrid mode $1^{\prime}$, $x^{s,cr}_{hyb}=38.2$ pm and $Q_{s}=1316.5$ result in $\beta^{\rm s}_{hyb1}=7.1(\pm 1.6)\times 10^{21}$ $\rm N/m^{3}$. Similarly for hybrid mode $2^{\prime}$, $\beta^{\rm s}_{hyb2}=7.6(\pm 1.5)\times 10^{21}$ $\rm N/m^{3}$ and for hybrid mode $4^{\prime}$, $\beta^{\rm s}_{hyb4}=7.1(\pm 1.2)\times 10^{21}$ $\rm N/m^{3}$. It is remarkable to note that the hybrid modes of SiNx are well described by a Duffing-oscillator model and therefore, induced SiNx nonlinearities can be effectively described by Duffing constants for hybrid modes.

IV S.4. Giant, induced nonlinearity measured on SiNx surface

Here we provide few technical justifications for our usage of the term giant nonlinearity for Duffing constant as measured on SiNx surface of graphene-SiNx hybrid modes. Our justification is based on two estimates, all of which show orders of magnitude changes: (i) a comparison of nonlinear threshold of our hybrid modes to that of bare SiNx resonators without hybridization, as measured by different groups Zwickl-08 , (ii) a comparison of average thermal displacement and threshold displacement for the onset of nonlinearity for bare graphene, SiNx and hybrid SiNx modes show four orders of magnitude reduction in ratio for hybrid modes. We further discuss the importance of having induced and tunable nonlinearity of SiNx resonator modes.

IV.1 A. Comparison of nonlinear threshold displacements for bare and hybrid SiNx resonators:

The threshold displacement corresponding to the bare SiNx $(x^{s,cr}_{bare}$) is extracted using following relation:

[TABLE]

where $\beta^{s}_{hyb}$ and $\beta^{s}_{bare}(=5\times 10^{12}$ $\rm N/m^{3})$ denotes hybrid and bare (following ref. Zwickl-08 ) Duffing constant of SiNx. For mode 1, $\beta^{s}_{hyb1}=7.1\times 10^{21}$ $\rm N/m^{3}$ and $x_{hyb1}^{s}=38$ pm results in $x_{bare}^{s,cr}=1.2$ $\rm\mu m$. Similarly for mode 2, $x_{bare}^{s,cr}=1.4$ $\rm\mu m$ and for mode 3, $x_{bare}^{s,cr}=1.2$ $\rm\mu m$. One can therefore note that such estimated displacement for onset of nonlinearity is 5 order of magnitude larger than that of hybrid SiNx.

IV.2 B. Comparison of thermal displacement and displacement corresponding to nonlinear threshold:

The ratio of nonlinear threshold ($x^{g,cr}_{bare}$) and thermal displacement ($x^{g}_{th}$) for bare graphene is ($4.6\times 10^{-9}/1.04\times 10^{-12}$=) $4.4\times 10^{3}$. Using our parameters and results of  ref. Zwickl-08 , in case of bare SiNx, the ratios are $(1.4\times 10^{-6}/24.5\times 10^{-15}=)5.9\times 10^{7}$, $(1.2\times 10^{-6}/23.7\times 10^{-15}=)5.0\times 10^{7}$ and $(1.2\times 10^{-6}/25.4\times 10^{-15}=)4.7\times 10^{7}$ for mode 1, 2 and 3 respectively.

However, for the same hybrid mode measured on graphene, ratio is ($1.0\times 10^{-9}/1.04\times 10^{-12}$=) $1.0\times 10^{3}$ and ($1.2\times 10^{-9}/1.04\times 10^{-12}$=)$1.1\times 10^{3}$, same order as that of bare graphene.

In case of SiNx hybrid modes the ratio drops by four orders of magnitude to $(38.2\times 10^{-12}/24.5\times 10^{-15}=)1.6\times 10^{3}$, $(30.4\times 10^{-12}/23.7\times 10^{-15}=)1.3\times 10^{3}$ and $(31.5\times 10^{-12}/25.4\times 10^{-15}=)1.2\times 10^{3}$.

IV.3 C. Relevance of induced nonlinearity of SiNx modes:

SiNx resonators have shown significant promise of observing quantum mechanical behavior for high-Q mechanical resonators at room-temperature. However, one needs to engineer nonlinearty in such a quantum device, to make it useful. After all, fluctuations of a classical resonator in thermal state is similar in shape in phase space to that of fluctuations of a linear harmonic oscillator deep in the quantum regime, dominated by zero point motion. For the resonator to be useful for precision measurement, one requires to squeeze the fluctuations in one quadrature: this require nonlinear interactions. Similarly, for gate operations in information devices, it is necessary to have conditional switching and phase shifts, both of which require nonlinear interactions between modes.

V S.5. Theoretical model

In this section, we analyze the model and find that the nonlinear system can be described by hybrid modes to some extent, akin to that of normal modes for a corresponding linear system.

V.1 A. Coupled linear SiNx and nonlinear graphene resonator

Our model is based on coupled modes of graphene and SiNx resonators, denoted by 1-dimensional amplitudes $x_{g}$ and $x_{s}$, respectively and is described by the set of equations:

[TABLE]

where $\gamma^{k}_{\rm bare}$ and $\omega_{k}$ $(k=g,s)$ represent linear damping and frequency of graphene and SiNx modes. Nonlinearity of graphene is quantified with two parameters: nonlinear damping $\eta^{g}_{\rm bare}$ and a cubic nonlinear response, characterized by its Duffing coefficient $\beta^{g}_{\rm bare}$. The graphene mode is bilinearly coupled to a SiNx mode which is modeled by an effective interaction Hamiltonian, $H_{\rm int}=\alpha x_{g}x_{s}$, where $\alpha$ is a coupling constant. SiNx is considered to be a linear oscillator in the range of forcing that we apply in our experiments.

V.2 B. Normal modes at low-amplitudes: probe on graphene and on SiNx resonators

At low external forcing, one can ignore nonlinear terms and thereby define two normal modes $x_{1}$ and $x_{2}$. These modes extend over the entire device. However, we detect either on graphene ($x_{g}$) or on SiNx ($x_{s}$), which can be expressed as:

[TABLE]

The detected amplitudes of normal mode $x_{1}$ (or $x_{2}$) on grpahene or SiNx are scaled by the ratio of square-root of respective masses ($m_{s}/m_{g}\sim 10^{4}$). As a result, amplitude of normal mode 1 (mode 2) on graphene i.e. $x^{g}_{1}$ ($x^{g}_{2}$) is two orders of magnitude larger than the amplitude of the same mode, $x^{s}_{1}$ ($x^{s}_{2}$), detected on SiNx surface. Accordingly, we have two Duffing constants for mode 1 (mode 2): $\beta^{g}_{\rm hyb1}$ ($\beta^{g}_{\rm hyb2}$) detected on graphene and $\beta^{s}_{\rm hyb1}$ ($\beta^{s}_{\rm hyb2}$) detected on SiNx.

V.3 C. Perturbative estimation I: Difference in $\beta^{g/s}_{\rm hyb}$, measured on SiNx and graphene surfaces

The difference in scales of Duffing constants measured on SiNx and on graphene surfaces, can be understood in the following way: it can be noted that the nonlinear forcing ($F_{NL}$) of a hybrid mode is uniform all along the spatial extent of the mode. However, since the hybrid mode for our device has physically two different kinds of oscillators with varying masses and surface areas, the force can be expressed as: $F_{NL}=\beta^{g}_{\rm hyb}x^{3}_{g}=\beta^{s}_{\rm hyb}x^{3}_{s}$, as measured on graphene ($x_{g}$) or on SiNx ($x_{s}$). For a forcing $F_{0}$ and assuming a steady state amplitude of $x_{g,s}\sim F_{0}Q_{g,s}/(m_{g,s}\omega^{2}_{0})$ for graphene and SiNx, leads to an approximate ratio of the measured Duffing coefficients $\beta^{s}_{{\rm hyb}}/\beta^{g}_{{\rm hyb}}\sim(m_{s}Q_{g}/m_{g}Q_{s})^{3}\sim 10^{9}$. This is in accordance with our measured values of $\beta^{g,s}_{{\rm hyb}}$ on graphene and on SiNx and provides a simple explanation of the giant nonlinearity measured on SiNx resonator surface.

V.4 D. Perturbative estimation II: effective nonlinearity

To get an estimate of effective scaling of induced Duffing constant of SiNx hybrid mode, $\beta^{s}_{\rm hyb}$ to that of graphene’s mass ($m_{g}$), bare Duffing constant, $\beta^{g}_{\rm bare}$ and coupling $\alpha$, perturbatively, let us consider the following simplified equations:

[TABLE]

[TABLE]

where $\alpha_{i}=\alpha/m_{i}$, ($i=g,s$) and damping is ignored.

For uncoupled graphene mode ($\alpha_{g}=0$) and assuming $\beta^{g}_{\rm bare}x_{g}^{3}\ll 1$, standard perturbation methods yields a (zeroth-order) solution of the form:

[TABLE]

where $A_{g}$ is a constant set by initial conditions. Substituting this zeroth order expression of $x^{(0)}_{g}$ in equation S.9a, we arrive at

[TABLE]

where, we have recognized $\cos\omega_{0}t$ as $x^{(0)}_{s}/A_{s}$ ($A_{s}$ being $x^{(0)}_{s}$ when $x_{s}$ and $x_{g}$ are uncoupled) and defined:

[TABLE]

[TABLE]

Here, in the definition of $B$, we have ignored frequency correction.

Substituting the values of $A_{g,s}=F_{0}Q_{g,s}/m_{g,s}\omega_{0}^{2}$ in above expression cancels the common forcing ($F_{0}$) yielding an expression that depends only on the system parameters:

[TABLE]

where $\alpha$ is the coupling strength between two resonators, $\beta_{bare}^{g}$ is the Duffing constant of the bare grapene mode, $m_{s}(m_{g})$ is the mass of SiNx (graphene) resonator, $\omega_{0}$ is the on-resonant frequency of both the resonators, and $Q_{s}(Q_{g})$ is the quantity factor of SiNx (graphene) modes.

We further note that the sign of $\alpha$ determines whether the SiNx is effectively a soft or a hard nonlinear oscillator. From the expression of $\beta^{s}_{\rm hyb}$, one can then express an effective scaling as:

[TABLE]

VI S.6. Frequency comb I: estimating nonlinear coefficients

We first develop a numerical model that reproduce the experimental observation of the frequency comb. Simulations results thereby give us estimate of $\beta^{g}_{\rm bare}$ and $\eta^{g}_{\rm bare}$. Next, we develop a general methodology to estimate nonlinear coefficients from measured experimental spectra on a general resonator surface. We finally apply the methodology to estimate $\beta^{k}_{\rm hyb}$ and $\eta^{k}_{\rm hyb}$ ($k=g,s$), as measured on graphene or SiNx surface.

VI.1 A. Estimating nonlinear coefficient from simulated spectra

The parameters used in numerical simulation of equations S.7a,b were extracted by fitting the Brownian spectrum (Fig. S.2) of graphene (equation S.2) and are listed in Table I. By varying the free parameters i.e. $\beta^{g}_{\rm bare}$ and $\eta^{g}_{\rm bare}$, we carefully calibrate and match the spectra in the instability region, where the comb is generated. The flow diagram in Fig. S.11 describes the methodology of nonlinearity estimation.

Numerically, we observe that the asymmetry in the envelop of the generated comb increases when nonlinear coefficient ($\beta^{g}_{\rm bare}$) is increased (Fig. S.12a) while the slope of the envelop changes with non-linear damping coefficient ($\eta^{g}_{\rm bare}$) (Fig. S.12b). The overall asymmetric fan-like shape of the generated comb is therefore a result of interplay between nonlinear damping and Duffing nonlinearity.

VI.2 B. Methodology of nonlinearity estimation

Here we describe the methodology we use to estimate nonlinear coefficients from observed frequency combs on graphene and SiNx surfaces. When modes are driven at twice the resonance frequency, we observe parametric gain (Fig. S.13) in both the hybrid modes below a threshold pump voltage. Above threshold, in the self-oscillation regime, we observe mixing of modes. We attribute this mixing to nonlinearity in the system. Using amplitudes of newly generated modes, we estimate the corresponding nonlinear coefficients.

In particular, starting with amplitudes of four modes of frequency comb to be $A_{1}$, $A^{\prime}_{1}$, $A_{2}$ and $A^{\prime}_{2}$, such that the corresponding displacement (measured on graphene or SiNx surface) can be expressed as:

[TABLE]

where $k=g,s$ and $n\Delta$ $(n=1,3)$ is the separation of modes from the central frequency, $\omega_{0}$. Combining equation S.16 with equation S.7, terms corresponding to Duffing nonlinearity and nonlinear damping can be expressed as:

[TABLE]

Using rotating wave approximation and collecting the terms corresponding to modes at $\mp 3\Delta$ from central frequency, one gets:

[TABLE]

[TABLE]

Similarly, collecting terms corresponding to $\mp 5\Delta$ from equation S.7 and squaring, yields

[TABLE]

where $A_{3}$ and $A^{\prime}_{3}$ depicts amplitude of newly generated modes, emerging at $\mp 5\Delta$ from central frequency. Solving equation S.18a and S.19a for the nonlinear coefficients, we finally get:

[TABLE]

and,

[TABLE]

where $\beta^{k}_{hyb2}$ and $\eta^{k}_{hyb2}$ are nonlinear damping and Duffing nonlinear coefficient of the left ($\omega_{1}$) hybrid mode. Similarly solving equation S.18b and S.18b for the right ($\omega_{2}$) hybrid mode we get:

[TABLE]

[TABLE]

This method gives an estimate of the nonlinear coefficients of coupled hybrid graphene-SiNx mode from spectral measurements.

VI.3 C. Application I: Estimating nonlinear coefficients on graphene resonator surface

Based on the methodology discussed in appendix, we estimate values of $\beta^{g}_{hyb(2,1)}$ and $\eta^{g}_{hyb(2,1)}$ (corresponding to experimental data of Fig. 2b, main text) for every pump voltage above threshold (Fig. S.14).

Corresponding values of $\beta^{g}_{hyb(2,1)}$ and $\eta^{g}_{hyb(2,1)}$ as obtained from numerical simulation (corresponding to Fig. 2c, main text) are also plotted with pump voltage (Fig. S.15). There is an overall agreement between the experimental observations and numerical simulations.

This method gives an estimate of the nonlinear coefficients of coupled hybrid graphene-SiNx mode from spectral measurements.

VI.4 D. Application II: Estimating Duffing constant on SiNx resonator

We have already established graphene to be nonlinear resonator with many intriguing properties in parametric regime, which we expect to observe in SiNx at resonance with graphene. However due to huge mass of SiNx, the nonlinear signature is not easily detectable when probed on SiNx. The signature of parametrically driven graphene-Silicon Nitride hybrid mode is observed with small number of new generated modes (Fig. S.16). Looking at the asymmetry we can conclude about large nonlinearity.

We simulate using equations S.7a,b, where SiNx is treated as a linear resonator, i.e., $\beta_{s},\eta_{s}=0$ and observe multi-mode generation in SiNx spectrum (Fig. S.17), further validating our observation. The nonlinear coefficients estimated from simulation plots (Fig. S.17) using equations S.20-22 and S.23 turns out to be, $\beta^{s}_{hyb2}(\beta^{s}_{hyb1})$ = $4.9\times 10^{23}(4.7\times 10^{23})$ $\rm N/m^{3}$ and $\eta^{s}_{hyb2}(\eta^{s}_{hyb1})$ = $2.7\times 10^{16}(2.6\times 10^{16})$ $\rm Ns/m^{3}$, in harmony with the experimentally measured values.