Pressure-induced nonlinear resonance frequency changes for extracting Young’s modulus of nanodrums

The resonance frequency of ultra-thin layered nanomaterials changes nonlinearly with the tension induced by the pressure from the surrounding gas. Although the dynamics of pressurized nanomaterial membranes have been extensively explored, recent experimental observations show significant deviations from analytical predictions. Here, we present a multi-mode continuum model that captures the nonlinear pressure-frequency response of pre-tensioned membranes undergoing large deflections. We validate the model using experiments conducted on polysilicon nanodrums excited opto-thermally and subjected to pressure changes in the surrounding medium. We demonstrate that considering the effect of pressure on the nanodrum tension is not sufficient for determining the resonance frequencies. In fact, it is essential to also account for the change in the membrane’s shape in the pressurized configuration, the mid-plane stretching, and the contributions of higher modes to the mode shapes. Finally, we show how the presented high-frequency mechanical characterization method can serve as a fast and contactless method for determining Young’s modulus of ultra-thin membranes.

Recently, it was also shown that the tension-induced shift in the resonance frequency of 2D membranes as a function of applied pressure could be used to characterize elastic properties [21]. In contrast to AFM which requires large indentations to enter the nonlinear regime for estimating Young's modulus [22], it was demonstrated that the geometrically nonlinear regime may be easily reached using this approach by sweeping the pressure across the membrane. However, it was found that the fit of the pressure-frequency response, results in estimations of Young's modulus that are an order of magnitude lower than the well-accepted values in the literature [21,23]. Therefore, in order to use this method to characterize suspended ultra-thin materials that undergo large deflections, a more comprehensive model is required to describe the underlying physics.
Unlike models commonly used for characterizing pressurized nanodrums [24,25], here, we develop a model that accounts for a change in the static equilibrium shape of the vibrating membrane, mid-plane stretching, and in-plane and out-of-plane mode coupling. By comparing the model to the finite element method (FEM) simulations of nanodrums subjected to large deflections, we validate the analytic formulation. Finally, we acquire the nonlinear pressure-frequency response of polysilicon nanodrums experimentally and demonstrate that the suggested multi-mode model can be utilized to characterize (extract Young's modulus) from the high-frequency dynamic response of pressurized ultra-thin membranes.

Governing equations
We consider a thin circular drum with a diameter of 2R and thickness h R clamped at the outer edge, as shown in Fig. 1a. The drum is assumed to be made of a homogeneous and isotropic material of density ρ, Young's modulus E, and Poisson's ratio ν. The drum is also assumed to be subjected to an axisymmetric tension n 0 , and the pressure difference alongside it is Fig. 1 a Schematic of the membrane and its cross-sections. b Side-view of the undeformed membrane, and c deformed configuration under transverse loading due to pressure difference alongside the membrane. g is the cavity depth, and P ext and P cav are pressure above the membrane and inside the cavity, respectively p (see Fig. 1b, c). We only account for axisymmetric vibrations and suppose that the aspect ratio is very small (i.e., h/R < 0.001 [26]) such that the effect of bending rigidity can be discarded, and the motion can be modeled using membrane theory [27]. The equations of motion for a membrane subjected to large transverse displacement w, moderate rotations, are the dynamic equivalents of the von Kármán equations. In this context, the governing equations of a pre-stressed membrane with negligible bending rigidity are [27]: where: Here, u = [u; v] is the axial displacements; moreover, ∇ 2 w, ∇w and ∇u denote the Laplacian of scalar field w, the vector gradient of scalar field w and the tensor gradient of the vector field u, respectively. divu and div N are the scalar and vector divergences of vector field u and tensor field N. Finally, ∇w ⊗ ∇w corresponds to the tensor product between vectors ∇w and ∇w. Overdot indicates derivation with respect to time. Considering only axisymmetric vibrations the following set of equations can be obtained in cylindrical coordinates and in terms of radial (u) and transverse (w) displacements: where ,r and ,rr denote the first and second derivatives with respect to r .
To write Eqs. (3a, b) in dimensionless form, we introduce the following set of dimensionless variables: This would yield the following non-dimensional form of Eq. (3): (5b) In order to solve Eqs. (5a, b), we first expand the transverse displacement w and radial displacement u as follows: where q k and η p denote the unknown modal coordinates, and N w and N u are the number of generalized coordinates that will be retained in the analysis. Moreover, k and p are mode shapes associated with the transverse and radial displacements, respectively. These modes are chosen such that they satisfy the following eigenvalue problems that originate from the linear counterpart of Eqs. (5): with ω k and γ p being the transverse and in-plane dimensionless eigenfrequencies, respectively. We note that Eqs. (7) are Bessel equations, and their solution can be expressed as: where S is the domain of integration. By applying Galerkin technique, namely inserting Eqs. (8), and Eqs. (6) in Eqs. (5), multiplying Eq. (5a) by k and Eq. (5b) by p , and then integrating over the entire domain, the nonlinear partial differential Eqs. (5a) and (5b) reduce to: k r dr is the projection of p on the k-th mode of vibration, and the nonlinear modal coefficients a k pi , b p i j , and c k i jl are: In order to find η p in terms of q k , we rewrite Eq. (10b) as: Inserting Eq. (12) in Eq. (10a) [28] leads to: Eq. (13) is a reduced-order model, which accounts for mid-plane stretching and higher-mode interactions and can be used for investigating nonlinear pressurefrequency response of nanodrums.

Solution procedure
Generally, the force vector applied on a membrane can be split into a constant, and a time-varying part. The constant force leads to a static deflection, and the time-varying force vibrates the structure around the new configuration. Considering small-amplitude oscillations around large pressure-induced deflections, it is a decent approximation to separate the contribution into static (q s k ) and dynamic (q d k ) components as follows [12]: This separation allows us to separate Eq. (13) to a static and dynamic part as well. We can obtain the static deflection due to Q k , and obtain the dynamic equations around the statically deflected configuration which we will then use for obtaining the natural frequencies of the nanodrum in the deformed state and as a function of the applied pressure, as follows: where In order to solve Eq. (16a) for q s k , we use the software package MANLAB. The software uses the Asymptotic Numerical Method (ANM) [29], a continuation algorithm based on high-order Taylor series expansions, to solve nonlinear systems of equations. The ANM continuation method allows for solving nonlinear algebraic equations written in the form: where F denotes a set of nonlinear algebraic equations, is a vector of unknown coefficients, and λ is the continuation parameter. The method also requires the system of equations to be written with quadratic nonlinearities [29]. To perform continuation, we thus recast the static part of Eq. (16a) to: By solving Eq. (19) for q s k , one can find the shift in resonance frequencies due to pressure. We note that the fundamental frequency of the pressurized nanodrum u can be obtained by considering the pre-factors of q d k in Eq. (16b) which is denoted with the matrix A, defined as: where δ uv is the kronecker delta. Diagonalizing the matrix A leads to the modified frequencies and mode shapes of the bulged nanodrum with an initial bulged shape: where P are the eigenvectors associated with the frequencies u for the pressurized membrane.

Convergence study
As mentioned in Sect. 2, in order to better capture small-amplitude oscillations of a pressurized membrane, one must (i) linearize the dynamic response around the new static configuration, (ii) account for mid-plane stretching effect, and (iii) include contributions from out-of-plane modes. Thus, to emphasize the significance of these assumptions and to facilitate a more accurate comparison of these parameters in a pressurized configuration versus a flat configuration, in Fig. 2, we simulate the effect of pressure on the static and dynamic in-plane and out-of-plane displacements of a nanodrum with Figure 2a shows that the original flat configuration changes significantly with increasing pressure. This result suggests that when estimating frequency shifts due to pressure, the dynamic governing equation must be obtained around the new bulged shape (q s k = 0) rather than the flat configuration (q s k = 0). The application of pressure causes a nonlinear in-plane displacement field at high pressures, as depicted in Fig. 2b. In addition, Fig. 2c illustrates the slight differences in the first mode shape that emerge due to statically deformed shape as a function of external pressure obtained from the multimode model (w d 0 denotes the dynamic transverse deflection at the center of the membrane). One must bear in mind that the membrane is flat at p = 0 mbar configurations. As we will discuss next, neglecting these effects could yield wrong estimation of the resonance frequencies of pressurized ultra-thin membranes.
We also note that the model is obtained in modal coordinates. Nonetheless, the governing partial differential equations (Eq. (3)) were formulated with in-plane (u) and out-of-plane (w) displacement fields. Moreover, Fig. 2 depicts the results for u, w s , and w d , rather than q s , and q d . Thus, to retrieve u from the modal equations, one must use Eqs. (4), (6b), and (12) and obtain u(r, t) as follows: However, for the out-of-plane deflection, one should note that separation of static and dynamic modal coordinates results in the separation of static and dynamic transverse displacements. Using Eqs. (15), (6a) and (4), the transverse displacement can be then obtained as follows: where: In order to highlight the influence of the effects depicted in Fig. 2 on the estimation of the resonance frequencies, we benchmark in Fig. 3a-c the fundamental frequency derived from our model against FEM results. The simulations are performed for the same drum specifications as Fig. 2. In Fig. 3a, we simulate the frequency shift where only one single out-of-plane mode is retained in the analysis (N u = 0, N w = 1). Frequency shifts are found, by using Eq. (16b) for the undeformed (flat) as well as deformed configurations. The figure confirms our earlier prediction in Fig. 2a that using the equations around the deformed configuration substantially increases accuracy. When compared to FEM results, the model still has a maximum error of 9% at 1 bar.
To illustrate the effect of in-plane displacements on the predicted frequencies, we increase N u from 0 to 4 in Fig. 3b while retaining only one out-of-plane generalized coordinate (N w = 1) in the analysis. It can be seen that including more in-plane modes in the model results in a more accurate result (the error relative to FEM simulations decreases to 4% at 1 bar). This confirms the important role of mid-plane stretching when tracing the fundamental resonance frequency of nanodrums as a function of pressure.
Although the response from this model is already close to the FEM simulations, Fig. 3c shows that the accuracy can be even further improved by including more transverse degrees-of-freedom in the basis functions by increasing N w in Eqs. (16). This finally leads to a negligible difference between simulations based on the present multi-mode model and FEM results, which highlights the slight contribution of higher order out-of-plane modes to the estimated resonance frequency [30]. Notably, each of these three criteria contributes to model convergence, and neglecting them would reduce the accuracy the model. It should be also noted that, considering more out-of-plane displacements while ignoring in-plane displacements would not produce accurate results, as in-plane displacements are always required to account for the mid-plane stretching effect.

Single-transverse-mode approximation
Although increasing the number of transverse modes improves model accuracy, it also increases the complexity of the governing equations, necessitating the employment of numerical techniques to solve them. However, the model can be analytically solved using just one transverse mode. Figure 3 demonstrates that the single-transverse-mode approximation is precise enough for tracing the pressure induced resonance frequency shifts up to 1 bar. Considering only N w = 1 and N u = 4, and converting equations back to dimensional form using Eq. (4) one can find the following equations as a simplified version of Eq. (16): where (ν) is plotted in Fig. 4. Moreover,q s 1 = hq s 1 denotes dimensional static deflection at the center of the nanodrum andq d 1 = hq d 1 is the dimensional dynamic deflection at the center of the membrane. For this simplified scenario, the nonlinear pressure-deflection and pressure-frequency relationships are derived analytically as follows: where A, B, and C are defined for the current study in We finally compare the result of our reduced-order model (Eq. (26)) with analytic models available in the literature, which are often used for the dynamic In Table (1) we benchmark this new formulation against the analytic solutions available in the literature for pressurized nanodrums [24,25]. Significant discrepancies between our model and those of Refs. [24,25] can be attributed to the simplifying assumptions lifted in our study, including linearization about flat configuration, absence of mid-plane stretching, and out-of-plane modal interactions. As demonstrated in Fig. 3, the first two assumptions have a greater impact on the model's accuracy than the third one.

Experimental validation
To verify the applicability of the proposed formulation for the dynamic characterization of ultra-thin drums, we measure 33 nm thick polysilicon membranes measured by Atomic Force Microscopy (AFM), which we transfer over SiO2 substrate cavities, that are created using reactive ion etching with a depth of 350 nm. The 10 µm-diameter drums are placed in a vacuum chamber with variable pressure ranging from 50  (26) to 1000 mbar. We employ a modulated blue laser diode (λ = 405 nm) to push the polysilicon nanodrums into resonance (Fig. 6a). The suspended nanodrum modulates the intensity of the reflected red laser (λ = 633 nm), which is collected at a photodiode and analyzed by a Vector Network Analyser (VNA). Our setup includes a PID controller that monitors chamber pressure using a vacuum pump and gas supply (nitrogen). Figure 6b depicts the frequency spectra of a polysilicon drum at different pressures P = P ext − P cav , where P cav and P ext are the pressure outside and inside the cavity, respectively (see Fig. 1). As a result of increased pressure, we observe an apparent tension-induced increase in the nanodrum's resonance frequency. To obtain the pressure-frequency response, we sweep the pressure from 50 to 1000 mbar and fit the resonance peaks with the damped linear harmonic oscillator model to obtain resonance frequencies at each pressure values (See Fig. 6b). In Fig. 7a we show typical pressurefrequency measurement data which we fit using our model (E = 155 GPa, n 0 = 1.6 N/m). In contrast to Fig. 3, we see a minimum in the pressure-frequency response at P cav = 248 mbar, which corresponds to a flat configuration. In this configuration, the drum's  fundamental frequency is determined solely by pretension. By increasing pressure, the nanodrum deforms statically, and the resonance frequency varies nonlinearly as f ∝ 6 √ E p 1/3 . Knowing this relation, such pressure-frequency measurements can also be used to determine the nanodrum's Young's modulus.
To demonstrate the versatility of the method for characterizing Young's modulus of ultra-thin membranes, we repeated the same measurements on 13 nanodrums that were adequately sealed (or the leak rate was low enough relative to the experiment's data collection speed) and did not exhibit hysteresis in the pressure-frequency response when the pressure was swept up and down. From these measurements, we determined n 0 = 1.1 ± 0.5 N/m for the nanodrums. Figure 7b depicts the retrieved Young's moduli histogram. We see that the obtained average value E = 148 ± 7 GPa from our technique is comparable to the literature-reported uni-axial stretching test findings for thin polysilicon beams [31,32], further validating the accuracy of our model. We note that in our experiments, owing to the low aspect ratio of the fabricated devices, the bending rigidity is very small and the linear stiffness mediated by bending rigidity is proportional to k bending = Eh 3 /(12R 2 (1 − ν 2 )) = 0.02 N/m which compared to the linear stiffness due to the pre-tension of the nanodrum (k stretching = n 0 = 1.1 N/m) in the flat configuration is negligible. Therefore, consistent with our modelling approach, the motion of the nanodrum is dominated by its tension.
At high pressures, the damping also rises, which may make it difficult to observe resonance peaks. Yet, our experiments demonstrate that even at ambient pressure the oscillation is under-damped and the resonance peak is apparent. In addition, as depicted in Fig. 7a, the membranes reach the nonlinear regime at p ≈ 400 mbar, which is sufficient for characterization purposes.

Discussion
Employing the pressure-frequency response as a material characterization method has several advantages over commonly used techniques. Compared to AFM that requires considerable deflections to estimate Young's modulus, i.e., the force-deflection curve should reach a cubic regime [22] (F ∝ δ 3 ; F: tip force, δ: membrane's center deflection), pressurization using the surrounding gas requires a lower force for determining Young's modulus and distributes the load over the membrane, minimizing the chance of membrane failure. The method also reaches the nonlinear regime suitable for estimating Young's modulus with substantially less deflections and stresses, thus making it a more practical method for characterizing brittle nanomaterials such as polysilicon (see supplementary information S1). Furthermore, unlike nonlinear dynamic characterization, which uses large driving signals to push the membrane into a nonlinear Duffing regime and requires proper calibration of the vibration amplitude [11], this method is independent of the vibration amplitude and only measures the resonance frequency as a function of the applied pressure.
It is worth noting that the low aspect ratio of ultrathin membranes is an inherent feature of nanomechanical systems, which accounts for their display of nonlinear dynamics even at small amplitudes. In cases where a higher aspect ratio is present, plate models can be utilized instead of the membrane model proposed in this study. Yet, it is important to note that the bending rigidity essentially re-scales the linear stiffness of these systems, and in high-pressure regimes, the nonlinear stretching factors caused by pressure continue to dominate the Young's modulus characterization.
Several points should be kept in mind when employing the proposed method for characterization. First, in an optical detection scheme, the cavity depth has to be optimized so that the photodiode voltage is still linearly related to the motion at high amplitudes [33]. Moreover, gas leakage through the clamped boundaries and poor adhesion forces between the nanodrum and substrate should be avoided to ensure negligible leakage rates and to prevent significant variations of the internal cavity pressure during the measurement [21,34,35]. We note that the present model does not account for morphological imperfections such as wrinkles [36] and boundary slippage [37], that are proven to affect the stiffness of 2D material membranes.
Furthermore, when P cav > 0, the cavity pressure change with the membrane's static deflection, and the squeeze film effect [38] may influence the dynamics of pressurized ultra-thin drums. Although both effects contribute significantly to the tension induced in the membrane, their cumulative effect on the drum's frequency shift near 1 bar and in the high-pressure regime is negligible for our polysilicon nanodrums, keeping the extracted value of Young's modulus unchanged. We also note that neglecting the squeeze-film effect and the bending rigidity of a relatively thick membrane may result in an overestimation of pre-tension; therefore, accounting for these two effects can enhance the accuracy of the anticipated pre-tension value which was estimated to be high in our analysis due to neglecting these effects. As discussed in the supplementary information S2, by accounting for these two effects, the pre-tension value drops from 1.6 to 1.1 N/m.
Other than material characterization, the model developed in this work can be used for estimating the deflection and resonance frequency shifts of pressurized membranes. These include 2D material resonant pressure sensors [4], microphones [39], and biosensors [6], where correct estimation of deflection due to surrounding fluid is important for probing the stiffness of the membrane. Furthermore, the analytical models provided here are also useful for accurate estimation of the mass-density of 2D membranes [40], which is often achieved by fitting the resonance frequency curves of the membrane as a function of pressure with tension and mass-density as the fit parameters.

Conclusion
In summary, attempts to use frequency shifts of pressurized nanodrums for characterization purposes were inaccurate previously due to lack of appropriate mathematical models. To resolve this, a multi-mode continuum model for predicting the dynamics of ultra-thin pre-tensioned drums subjected to high pressures is presented here. By maintaining in-plane and out-of-plane membrane deformations, we investigate the pressuredependent resonance frequency and mode shapes and illustrate the crucial importance of mid-plane stretching in the high-pressure regime. We highlight the discrepancies by comparing the model's accuracy to previously reported analytical models and FEM simulations. By fitting the resonance frequency measurements of a series of polysilicon drums as a function of the applied pressure as determined by laser interferometry, we demonstrate the model's validity and applicability for determining Young's modulus of ultra-thin nanomaterial membranes using on-resonance high-frequency characterization. Since pressurization results in a distributed stress field compared to indentation, this method could be utilized as an effective toolset for determining the Young's modulus of brittle nanomaterials.