Nonlinear dynamic characterization of two-dimensional materials

Owing to their atomic-scale thickness, the resonances of two-dimensional (2D) material membranes show signatures of nonlinearities at forces of only a few picoNewtons. Although the linear dynamics of membranes is well understood, the exact relation between the nonlinear response and the resonator’s material properties has remained elusive. Here we show a method for determining the Young’s modulus of suspended 2D material membranes from their nonlinear dynamic response. To demonstrate the method, we perform measurements on graphene and MoS2 nanodrums electrostatically driven into the nonlinear regime at multiple driving forces. We show that a set of frequency response curves can be fitted using only the cubic spring constant as a fit parameter, which we then relate to the Young’s modulus of the material using membrane theory. The presented method is fast, contactless, and provides a platform for high-frequency characterization of the mechanical properties of 2D materials.


Supplementary Note 1: Equations of motion
In this section the mathematical derivation of equations (1) and (3) of the main text is presented. The strain energy of the circular membrane can be obtained as [1] where E is the Young's modulus, ν is the Poisson's ratio, h is the thickness and R is the radius of the membrane. Moreover, rr , θθ , and γ rθ are the normal and shear strains that are determined as where u, v and w are the radial, tangential and transverse displacements respectively.
For a membrane with fixed edges u and w shall vanish at r = R. Moreover, u should be zero at r = 0 for continuity and symmetry. Assuming only axisymmetric vibrations (v = 0 and ∂w/∂θ = ∂v/∂θ = ∂u/∂θ = 0) and fixed edges, the solution is approximated as [2] w = x(t)J 0 α 0 r R , (5a) Here it should be noted that for axisymmetric vibrations the shear strain γ rθ would become zero. In Supplementary Equations (5a,b), x(t) is the generalized coordinate associated with the transverse motion of the fundamental axisymmetric mode and q k (t) are the generalized coordinates associated with the radial motion. Moreover, J 0 is the Bessel function of order zero, and α 0 = 2.40483. In addition,N is the number of terms in the expansion of radial displacement, and u 0 is the initial displacement due to pre-tension n 0 that is obtained from the initial stress σ 0 = n 0 /h as follows: where the overdot indicates differentiation with respect to time t.
In the presence of transverse harmonic distributed pressure (p = F el cos(ωt)/R 2 π) of constant direction, the virtual work done is where ω is the excitation frequency and F el gives the force amplitude, positive in transverse direction. Higher-order terms in w are neglected in Supplementary Equation (8) [3].
The Lagrange equations of motion are and q = [x(t), q k (t)], k = 1, ...,N , is the vector including all the generalized coordinates.
Since radial inertia has been neglected, Supplementary Equation (9) leads to a system of nonlinear equations comprising of a single differential equation associated with the generalized coordinate x(t) andN algebraic equations in terms of q k (t). By solving theN algebraic equations it is possible to determine q k (t) in terms of x(t). This will reduce theN + 1 set of nonlinear equations to a single Duffing oscillator as follows: where it is found that and c is the damping coefficient that has been added to the equation of motion to introduce can be expressed in the form: where C 3 is dimensionless constant which is a function of the Poisson's ratio. The solutions for C 3 as a function of ν are plotted in Supplementary Figure 1 for values of the Poisson's ratio between 0 and 0.35.
The relation between C 3 and ν is best described with the inverse of a second-order polynomial, namely: This functional dependence is similar to the one used for AFM nanoindentation measurements, often referred to as q(ν) [4].
Next, the following dimensionless parameters are introduced: By using Supplementary Equations (10) and (14) the following dimensionless equation of motion can be obtained: where Supplementary Equation (15) is valid for studying nonlinear vibrations of membranes subjected to external harmonic excitation in the frequency neighbourhood of the fundamental mode if the fundamental mode of vibration is not involved in an internal resonance with other modes. If such condition retains, then other modes accidentally excited will decay with time to zero due to the presence of damping [5]. In this work, it is assumed that this condition is preserved and therefore the response of the membrane is described by a single Duffing oscillator for performing nonlinear parameter estimation.

Supplementary Note 2: Nonlinear identification
In order to obtain the coefficients of the Duffing oscillator (Supplementary Equation (10)) from the measured nonlinear response curves, here we utilize the harmonic balance method.
This method is a suitable and accurate mathematical technique that entails the solution of nonlinear equations to be approximated by a truncated Fourier series. In case of the dimensionless Duffing oscillator (i.e. Supplementary Equation (15)), a first order truncation has been shown to provide accurate results [5]. Hence, Substituting Supplementary Equation (17) into Supplementary Equation (15) yields: where A = x 2 1 + x 2 2 is the amplitude of motion. Moreover, x 1 = A sin φ and x 2 = A cos φ, φ being the phase difference between the excitation and the response. From Supplementary Equations (18a) and (18b) the following analytic frequency-amplitude relation could be found: The idea of harmonic balance based parameter estimation is to follow a reverse path [6].
In other words, the identification is conducted by assuming that the vibration amplitude A and frequency ratio r are already known for every frequency step from experiments.
Therefore, in order to obtain unknown parameters, the following system should be solved for every jth frequency step, r (j) : System (20) can be compactly written asĀ h · Y =B h . This system is over constrained since it contains 2 × m equations. In order to solve (20) and to estimate system parameters, least squares technique is used and the norm of the error should be minimized. Accordingly, here the pseudo-inverse of matrix A h is calculated and the solution is obtained as follows: A problem in utilizing the least squares technique is that the identified peak amplitudes in the frequency-response curves do not correspond to the ones obtained from the experiments.
In order to resolve this issue, a correction on the quality factor is made by making use of the following expression (see [6] for the derivation details): in which A max is the experimentally measured peak amplitude for each frequency-amplitude curve. This will yield the nonlinear identification procedure to a single fit parameter algorithm for the estimation of η 3 .

Supplementary Note 3: Estimation of the electrostatic spring softening
In this section the effect of nonlinearities in the electrostatic force on the nonlinear spring constant is estimated.
The electrostatic force acting on the membrane is given by where U el = 1 2 C g V 2 is the electrostatic energy, V = V dc + V ac cos(ωt) is the applied voltage and C g is the gate capacitance. Assuming x << g, where g is the gap between the membrane and the backgate, the gate capacitance can be approximated using a parallel plate capacitor model: where R is the radius of the membrane and ε 0 is the vacuum permittivity. The resulting electrostatic force is given by If we expand this expression around x = 0, using x/g << 1, we get The first term ( 1 g 2 ) is the dominant electrostatic actuation term and the second term is what is usually described as a spring softening term ( 2x g 3 ). This term influences only the linear spring constant of the resonator. The term including x 3 will have a softening effect on the cubic spring constant: k 3,soft = 2 g 5 ε 0 R 2 πV 2 dc (for V dc >> V ac and g >> x). The resulting cubic spring constant k 3,tot will be given by The ratio of the two contributions (using the expression for k 3 from Supplementary Equation (12)) is For a Young's modulus E = 594 GPa, radius of R = 2.5 µm, thickness of h = 5 nm, gap size g = 385 nm and V dc ≤ 3 V, this ratio becomes which means that the electrostatic softening will have a negligible effect on the extracted Young's modulus (resulting in an error of < 0.1%). It should be noted that the cavity depth has a significant influence on the effect of electrostatic softening of the cubic spring constant.
To get reasonable error margins (< 5%), the ratio of Supplementary Equation (29) should be kept above 20.

Supplementary Note 4: The effect of the Young's modulus on the strength of the nonlinear dynamic response
In Supplementary Figure 2

Supplementary Note 5: Determination of the driving force
In this section the electrostatic force is determined by two methods: (i) based on the geometry and the voltage applied on the drum and (ii) based on the measured amplitudes at resonance.
The electrostatic force acting on the membrane at a frequency ω, assuming a parallel plate capacitor model, can be calculated as: where g is the gap size and ξ = 0.432 is a correction factor that accounts for the fundamental mode shape of a fixed circular membrane (see Section 1 and main text).
In a homodyne detection scheme, only the motion resulting from a force at frequency ω can be detected. Therefore, the effective force on the drum will be: It is worth mentioning that due to the fact that the output of the VNA has 50Ω impedance and the drum is a capacitor (high impedance), the ac voltage felt by the drum is double the voltage that is output by the VNA (V ac = 2V ac,VNA ).
Knowing the applied driving force is crucial for fitting the nonlinear response curves and finding the Young's modulus of the membranes, and it therefore needs to be determined Finally, due to the presence of the AuPd electrode underneath the drum, the electric field will be screened by the edges of the electrode around the perimeter of the drum. This will lower the effective capacitance of the drum and therefore also the effective electrostatic force felt by it. For a flat drum, using a gap size of 385 nm and a 100 nm-thick AuPd electrode, such fringe field effects can account for up to 15 % of overestimation of the force.
To correct for these effects, we calculate the expected peak amplitude at resonance (x calc,RMS | ω=ω 0 ) using the force calculated using Supplementary Equation (31). In the linear vibration regime, this function is expected to be linear with the force, hence with the applied ac voltage: The slope of the x RMS vs. V ac,RMS curve at low ac voltages (linear regime) will be a constant β, defined as: Knowing β, we can calculate the force for any driving ac voltage: The actual force felt by the drum (F RMS ) will be modified depending on the uncertainties in determining V dc (V dc ) and g (g ), possibly resulting in a different slope of the x RMS vs.
V ac,RMS curve: β . We can use this value to cross-check and correct for the actual force (F RMS ): To visualize the procedure, in Supplementary Figure 3, we present data from a 4-µm in diameter, 8-nm thick FLG drum. The calibrated amplitude at resonance (or x RMS,max in the nonlinear regime) is plotted as blue dots (similar data can be seen in [11]). The calculated expected amplitude based on Supplementary Equation (33) is represented by the red line.
To correct for the difference in slopes in the linear part (left hand side), we fit the slope of the measured amplitude to extract β . The resulting error in the force estimation is roughly 10 %, which, assuming that the error is due to trapped charges, translates to a 0.1 V offset over the applied V dc = 1 V. This offset can easily stem from residual charges on the graphene flake [7].
As long as the force scales linearly with the applied driving voltage (and there are no significant nonlinearities in the actuation), this method can be easily extended to calibrate the force of other actuation methods, like optothermal actuation, since the force can be extracted from the calibrated amplitude data. This is provided that the Q-factor is known and the displacement x can be properly calibrated.

Supplementary Note 6: Calibration of the amplitude
In this section the procedure of calibration of the motion amplitude using the resonator's Brownian noise is explained. In addition, the effect of nonlinearities in the readout method on the calibration is estimated.
To convert the measured signal (in Volts) to the motion amplitude of the resonator (in metres), we use thermal calibration (we follow closely the procedure described in [12]). This means that we can relate the time-averaged undriven (Brownian) motion of the drum to its thermal energy. Using the equipartition theorem, the following relation between the drum's amplitude and its thermal noise power spectral density (PSD) can be stated: where x 2 n (t) is the mean-square amplitude of motion of the n-th mode of the drum, f is the frequency, and S zz (f ) is the measured one-sided spectral density of the drum's motion.
The thermal noise PSD of the drum (S zz (f )) is given by [12]: where k B is the Boltzmann constant, T is the temperature (in our case, we take T = 293 K), and f n , m eff,n and Q n are the resonance frequency, effective mass (m eff = 0.2695 m at the centre of the drum [12]) and quality factor of the n-th resonance mode respectively.
We start by experimentally measuring the PSD using a spectrum analyser, S VV (f ). This represents the total PSD of the system. This includes noise from various sources in the detection system, among which the dark current noise of the photodiode (S PD VV ), which can be easily measured and eliminated. Assuming all other sources of noise to be white, we can consider them to be a flat additive contribution(S w VV ) to the total PSD. S w VV also determines the noise floor of the measurement. The remaining signal can be entirely attributed to the thermal noise PSD of the drum (S zz (f )) multiplied by a transduction factor α (V 2 /m 2 ), which depends on the mechanical-optical-electrical transduction of the system. It is worth mentioning that α is highly susceptible to changes in the measurement laser intensity, po-sition and focus of the drum, the cavity depth and optical/electrical losses in the system.
To ensure that the value of α stays constant throughout the measurements we take all measurements in one go, without changing the sample position.
We can now write a general relation between the aforementioned quantities: In α. This fit is shown by the black curve in Supplementary Figure 4 (b). From this fit, we can determine 3 fit parameters: the resonance frequency f 0 , the quality factor Q and, most importantly, the transduction factor α (in this case α = 2.56 · 10 12 V 2 /m 2 ).
Using this transduction factor, we can now determine the amplitude (in nanometres) of the driven motion as well. The magnitude of the driven motion at resonance is obtained from a two-port network analyser measurement and is given by This amplitude can be converted to metres (x RMS (f )) considering a constant transduction factor α by: Substituting γ = 4π g λ , C 0 = A+B cos 4π g λ , we get: Using homodyne detection with a very narrow bandwidth filter (in our case < 1kHz), we are only sensitive to the cos ωt component of Supplementary Equation (48) (labeled f ), hence the actual intensity detected by the VNA at the resonance frequency ω 0 is: The linear coefficient of Supplementary Equation (49) corresponds to the transduction coefficient extracted using the calibration procedure √ α. At small amplitudes, the cubic component plays little role. The error that the cubic term introduces in the amplitude calibration at large amplitudes is therefore given by: or: Hence, using a laser with λ = 632.8 nm, for amplitudes of up to x max = 10 nm, the error that we make in the estimation of the amplitude is |ε| < 0.5 %. We note that this error is independent of the the gap size g.
Another important parameter that may influence the calibration is the actual temperature of the membrane. Even though the measurements are conducted at room temperature, the measurement laser locally heats up the membrane. For this reason, the laser power throughout the measurement is kept relatively low (0.42 mW). To estimate the temperature increase due to the laser heating, we perform measurements at varying laser power. From the measured resonance frequency we can determine the strain of the membrane (ε). This measurement is shown in Supplementary Figure 6.
As expected, the strain grows linearly with the laser power, due to a finite temperature increase ∆T . This can be modelled with a simple equation: where ε 0 is the mechanical pre-strain, α T is the thermal expansion coefficient and ∆T is the temperature increase. The second term in the equation is the thermally induced strain.
Since α T for graphene is negative [14], the total strain of the drum increases with increasing laser power.
By fitting the measured data with a linear fit, we can extrapolate the curve and estimate ε 0 (see Supplementary Figure 6). Using α T ≈ −6 · 10 −6 K −1 , we find the temperature increase at 0.42 mW laser power to be less than 0.8 K. This temperature difference, following Supplementary Equations (39) and (42), results in about 0.1 % error in the amplitude calibration.