Introduction

High-Q nanomechanical resonators play a central role in sensing and enable ultrasmall mass, acceleration, and force detection. Nonetheless, due to their nanoscale size and exceptional isolation from the surrounding environment, even minute forces, as small as a few piconewtons, can induce large-amplitude oscillations in them and result in a plethora of nonlinear phenomena that include bi-stability1,2,3,4, parametric resonance5,6,7, self-oscillations8,9,10, and mode-coupling11,12,13,14. Many of these nonlinear phenomena can provide new information that is absent in the linear regime of operation. For instance, nonlinear resonances can be used to characterize nanomaterial properties15,16, enhance frequency stability12,17, or generate mechanical frequency combs18,19. Some of these nonlinear phenomena have been engineered by leveraging the interplay between geometric and electrostatic nonlinearities6,20. However, the introduction of external competing nonlinear forces may give rise to a series of unwanted side-effects, including noise21,22 or back-action effects23 that can further complicate the nonlinear dynamic behavior, device fabrication, and operation. Therefore, methodologies that can tailor the dynamic characteristics of nanomechanical devices solely through geometric design in the fabrication stage are highly desirable.

Although numerous studies have already demonstrated the design optimization of resonance frequencies and Q-factor of nanomechanical resonators24,25,26,27,28,29, the influence of geometric design on nonlinear dynamics has rarely been investigated11,30,31. Here, we show that soft-clamping techniques that are utilized to realize high-Q nanomechanical resonators can also be engineered to tune nonlinear dynamics. By manipulating the support boundary in high-stress Si3N4 string resonators, we can tune the stress field and induce strong in-plane to out-of-plane coupling to simultaneously increase the Q-factor and the onset of nonlinearity over three times that of a doubly clamped string. Furthermore, by changing the support angle, we show that it is possible to engineer compressive forces in the soft-clamped resonators and achieve buckled configurations in a controllable manner. These buckled states allow us to maximize geometric nonlinearity and change the response from hardening to softening. To understand the conditions required for strain engineering and buckling, we develop reduced-order models from finite element (FE) simulations, which highlight the role of the support angle in tuning the nonlinear dynamic response. Our results thus provide experimental evidence of controllable nonlinear dynamic engineering of nanomechanical resonators solely by geometric design and pave the way for integrating arrays of highly tunable nonlinear nanodevices on a single chip32,33,34,35,36.

Results

Nonlinear dynamic characterization

Figure 1a shows a scanning electron microscope (SEM) image of a nanomechanical resonator we fabricated with high-stress Si3N4 (see “Methods” for more details). The central string resonator of all studied devices had the same length (L = 200 μm) and width (w = 2 μm). However, the geometric parameters of the support beams, namely the support length Ls, the width ws, and the angle θ (see Fig. 1b), are varied to investigate their effect on the nonlinear response. To characterize the nonlinear dynamics of the string resonators, we fix a chip with nanodevices to a piezo actuator that provides harmonic drive in the out-of-plane direction. We use a Zurich Instruments HF2LI lock-in amplifier connected to an MSA400 Polytec Laser Doppler Vibrometer that is focused at the center of the string resonator to carry out frequency sweeps (see Fig. 1c). We perform all the measurements at room temperature and in vacuum chamber with pressure below 2 × 10−6 mbar to improve the signal-to-noise ratio and minimize air damping.

Fig. 1: Nonlinear dynamic characterization of string nanoresonators.
figure 1

a Scanning electron microscope (SEM) image, colored in blue, of a string resonator with Ls = 110 μm, θ = 0.2 rad. b Illustration of design parameters. c Schematic of the measurement set-up comprising a Micro System Analyzer (MSA) Laser Doppler Vibrometer (LDV) for reading out the motion and a piezo-actuator for generating the excitation force. d Duffing nonlinear response of the string resonator with ws = 1 μm, Ls = 130 μm, and θ = 0 as a function of the drive level. e Sensitivity of the Duffing response to the support length Ls for string resonators with ws = 1 μm and θ = 0. The measurements are conducted on four string resonators with different Ls marked with different colors. The fitted backbone curves are shown in red, and the fitted β values for Ls = (30, 70, 110, 150) μm are β = (8.28, 2.12, 0.96, 0.62) ×  1022 m−2 s−2, respectively. The mass normalized excitation levels Fexc in (c) and (d) are indicated by the color scale.

Full size image

To probe the geometric nonlinearity of our devices, we conduct frequency sweeps at different drive levels and measure the vibrations of the central string. Fig. 1d shows the frequency response at various drive levels for a device with θ = 0. We note the presence of a hardening-type nonlinearity at large amplitudes that arises from the elongation of the string during vibrations. To quantify the observed nonlinearity, we use the Duffing equation:

$$\ddot{q}+\mu \dot{q}+\alpha q+\beta {q}^{3}={F}_{{{{{{{{\rm{exc}}}}}}}}}\sin (2\pi ft).$$
(1)

where q is the displacement of the center of the string, \({F}_{{{{{{{{\rm{exc}}}}}}}}}\sin (2\pi ft)\) is the mass normalized harmonic drive force, \(\alpha ={(2\pi {f}_{0})}^{2}\), μ = 2πf0/Q are the mass-normalized linear stiffness damping coefficient, respectively. Furthermore, β is the Duffing constant that we extract by fitting the backbone of the experimental frequency response curves using the expression: \({f}_{\max }^{2}={f}_{{{{{{{{\rm{0}}}}}}}}}^{2}+\frac{3}{16{\pi }^{2}}\beta {A}_{\max }^{2}\), where \({f}_{\max }\) is the drive frequency corresponding to the maximum amplitude \({A}_{\max }\)37,38 (see Supplementary Note 1). In Fig. 1e, we quantify the change in the Duffing constant β when varying the support length Ls with θ = 0. We observe a factor 13 reduction in the extracted β when increasing the support length Ls from 30 μm to 150 μm. This is because long supports offer less rigidity to the in-plane motion, thus allowing the central string to relax when vibrating at large amplitudes and consequently reducing the overall geometric nonlinearity.

Theoretical model showing the influence of soft-clamping

The reduction in the Duffing constant observed in Fig. 1e highlights the role of the in-plane stiffness on the geometric nonlinearity. Hence, to better understand the influence of support stiffness, we perform additional experiments on a large number of string resonators with different support length Ls and width ws and extract their Duffing constant β. Moreover, we develop a simplified model consisting of a string connected to in-plane springs at both ends, denoted as kin (see Fig. 2a), to quantitatively capture the influence of Ls and ws on kin and thus on β when θ = 0. We particularly model the boundary springs as doubly clamped beams with pre-tension σ0 = 1.06 GPa and use their central deflection to analytically estimate kin as follows29:

$${k}_{{{{{{{{\rm{in}}}}}}}}} ={\left({k}_{{{{{{{{\rm{Ei}}}}}}}}}^{-1}+{k}_{{{{{{{{\rm{si}}}}}}}}}^{-1}\right)}^{-1}+{k}_{\sigma {{{{{{{\rm{i}}}}}}}}}\\ ={\left[{\left(\frac{16Eh{w}_{{{{{{{{\rm{s}}}}}}}}}^{3}}{{L}_{{{{{{{{\rm{s}}}}}}}}}^{3}}\right)}^{-1}+{\left(\frac{2Eh{w}_{{{{{{{{\rm{s}}}}}}}}}}{(1+\nu )\eta {L}_{{{{{{{{\rm{s}}}}}}}}}}\right)}^{-1}\right]}^{-1}+\frac{4{\sigma }_{0}h{w}_{{{{{{{{\rm{s}}}}}}}}}}{{L}_{{{{{{{{\rm{s}}}}}}}}}},$$
(2)

where kEi is the bending stiffness, ksi is the shear stiffness, kσi is the additional contribution from the pre-tension. Moreover, η is a geometric factor that, for a rectangular cross-section, is η = 1.239. In contrast to the central part, which is assumed to be a vibrating string, the support beam is modeled as a moderately thick beam (ws/h ≥ 2.94), for estimating kin. To validate our analytical estimation of the in-plane stiffness kin, we also obtained it numerically using FE simulations. The analytical (lines) and the FE (circles) results in Fig. 2b closely match one another and confirm our earlier prediction that wider and shorter support beams offer more rigidity against the deflection, particularly in the in-plane direction. Next, to capture the effect of in-plane springs on the in-plane to out-of-place coupling and thus nonlinear dynamics, we obtain the Lagrangian \({{{{{{{\mathcal{L}}}}}}}}=T-{U}_{{{{{{{{\rm{s}}}}}}}}}-{U}_{{{{{{{{\rm{k}}}}}}}}}\), where T is the kinetic energy of the string, Us and Uk are the potential energy of the string and the two springs kin, respectively, and use Lagrange equations to obtain the updated α and β as follows (see Supplementary Note 2 for details):

$$\alpha =\frac{{k}_{1}}{{m}_{{{{{{{{\rm{eff}}}}}}}}}}=\frac{{\pi }^{2}(1-\nu ){\sigma }_{0}}{\rho {L}^{2}}{\left(1+\frac{2EA}{{k}_{{{{{{{{\rm{in}}}}}}}}}L}\right)}^{-1},$$
(3)
$$\beta =\frac{{k}_{3}}{{m}_{{{{{{{{\rm{eff}}}}}}}}}}=\frac{{\pi }^{4}E}{4\rho {L}^{4}}{\left(1+\frac{2EA}{{k}_{{{{{{{{\rm{in}}}}}}}}}L}\right)}^{-1},$$
(4)

where k1 and k3 are the linear stiffness and the Duffing constant before mass normalization, respectively, meff = ρAL/2 is the effective mass, and A = hw is the area of the string’s cross-section. It is worth noting that \({(1+2EA/{k}_{{{{{{{{\rm{in}}}}}}}}}L)}^{-1}\) serves as a tuning factor introduced by the finite in-plane stiffness kin, which changes both α and β of a doubly clamped string with stress at the same rate. The effect of kin on the nonlinear stiffness can intuitively be understood by realizing that the geometric nonlinear stiffness of a string resonator is due to the increase of average string length and the resulting tension proportional to q2. If the clamping at the ends of the string is weakened by reducing the value of springs kin, then the tension increase for the same length increase will be less, such that the nonlinear stiffness β will reduce, as shown by Eq. (4).

Fig. 2: Influence of in-plane stiffness on the geometric nonlinearity of string resonators.
figure 2

a Simplified model of a string vibrating in the out-of-plane direction, denoted as z, with finite in-plane stiffness kin. b Analytical (solid lines) and finite element (FE) based (circles and inset) results of the in-plane stiffness kin for a support beam with varying Ls and ws. c Analytical (solid lines) and measured (diamonds) β of a string resonator with varying Ls and ws corresponding to (b). The dashed line shows the analytical estimation of β for a doubly clamped string resonator.

Full size image

In Fig. 2c, we show the extracted Duffing constant β from experiments via fitting the backbone curves of the frequency responses and compare those to the analytical model predictions from Eq. (4). From Fig. 2b and Fig. 2c, it is apparent that the variation in β matches the model quite well. We shall note that for the tested device with the slenderest support (Ls = 150 μm, ws = 1 μm) and the central string stress of 5.08 MPa, the analytically derived f0 = 94 kHz and β = 7.21 × 1021 m−2 s−2 compare well with the measured counterparts f0 = 103 kHz and β = 6.20 × 1021 m−2 s−2. Therefore, the assumptions of having high-stress strings with sinusoidal eigenmode are valid for obtaining the analytical expressions in the case of θ = 0. By comparing the dashed line to the experimental values in Fig. 2c, it can be seen that the geometric nonlinearity β can be reduced by up to two orders of magnitude using the presented support design. This substantial reduction in geometric nonlinearity highlights the role of support design in tailoring nonlinear dynamic behaviors in nanomechanical resonators with soft-clamping. It is also noteworthy that the β calculated by using the simplified model converges to the value of a doubly clamped string (1.33 × 1024 m−2 s−2) when the support length Ls tends to zero (kin approaches infinity in Eq. (4)), thus further confirming our model. However, in the case of wider and shorter supports, we notice that the simplified model deviates from the measurements. We attribute this to the fact that the dimensions of the supports become comparable to the underetch distance of the Si3N4 (~5 μm) such that the assumptions used for the derivation of Eq. (2) might not be valid anymore.

Engineering the nonlinearity from hardening to buckling-induced softening

To gain a deeper understanding of the full potential of support design on nonlinear dynamics, we also look into the influence of the support angle θ on the nonlinear frequency response curves. By changing θ from positive to negative, we are able to tune the tilting direction of the backbone curves around the resonance, from the common hardening nonlinearity to softening (see Fig. 3a). To understand the physical mechanism behind this observation; we use Keyence Digital Microscope VHX-6000 to focus at the middle of central string and the unreleased Si3N4 layer to measure the difference H of their focal heights, as shown in Fig. 3b. We note a maximum deviation of H =  22.17 μm for the device with Ls =  150 μm and θ = −0.1 rad, which suggests the presence of broken-symmetry in nanomechanical resonators with θ < 0. We attribute this to a change of built-in stress in the Si3N4 resonator from tension to compression, which, upon surpassing the buckling bifurcation point, breaks the out-of-plane symmetry and yields a buckled configuration (see SEM image in Fig. 3b). To verify these observations, we simulate the buckled resonator response by nonlinear reduced-order modeling of full FE models40 and numerical continuation41 (see Supplementary Note 3), which had been successfully applied to model the nonlinear dynamics of graphene drums40. It is worth mentioning that for a buckled string, the maximum amplitude does not always occur at the center of its first out-of-plane symmetric mode. Accordingly, we use FE simulations to obtain the amplitude ratio between the center of the mode and where it has the maximum amplitude, thus scaling the measured amplitudes at the center. The simulated results are shown as solid curves in Fig. 3a and demonstrate that the buckled configurations can account for the experimentally observed softening response (see Supplementary Note 3)42.

Fig. 3: Influence of support angle θ on the geometric nonlinearity.
figure 3

a Finite element (FE) based (bold lines) and experimentally measured (triangles) frequency response curves of resonators with ws = 1 μm and Ls = 150 μm, showcasing the shifting between hardening and softening induced by the support angle θ. The solid parts of bold lines represent the stable branch of the simulated response, while the dotted parts are unstable. The legends show the mass-normalized drive level. The fitted β values for θ = (−0.3, −0.5)rad are β = (−0.56, −0.27) × 1022 m−2 s−2. For θ = (0.3, 0.5)rad, the Duffing constant β = (7.00, 9.01) ×  1022 m−2 s−2. b Finite element (FE) based results (dots) and measurements (diamonds) of the buckling induced static displacement height H of the string at its center with ws = 1 μm for different support angles θ and lengths Ls. The inset shows the SEM image, colored in blue, of an array of buckled string resonators with ws = 1 μm, θ = −0.2 rad, and different Ls from 150 μm to 50 μm.

Full size image

To further investigate the role of the support angle θ on the tunability of the dynamical properties, we perform additional measurements on string resonators with values of θ ranging from -0.5 rad to 0.5 rad while keeping Ls = 150 μm and ws = 1 μm constant. In Fig. 4, we show the variation of the resonance frequency f0, Q-factor, and Duffing constant β of the first symmetric out-of-plane mode as a function of θ. We note that supports with positive θ significantly increase the values of the dynamical parameters shown in Fig. 4. This is attributed to higher kin that results in higher tension and thus translates into higher values of f0, Q-factor, and the Duffing constant β (see Eq. (4)). The most intriguing observation, however, lies in the region where the devices transition near θ = 0 from a flat configuration to the buckled state. Here, we notice a sudden increase in the resonance frequency and the maximum ∣β∣ of our devices, which is due to the large offset from the flat state (See Fig. 4a, c). By reducing θ towards -0.5 rad, however, the post-buckling offset is found to decrease again, and subsequently, both f0 and ∣β∣ decrease monotonically. We noticed that near θ = 0, the Q-factor of our devices drops to a similar level to that of stress-free string resonators43, whose dissipation dilution disappears with the relaxation of high tension (see Fig. 4b)27,29. However, we note that the FE-simulated Q-factor of devices close to the onset of buckling is lower than the intrinsic Q-factor Q0 = 9864. This could be attributed to the anti-spring behavior of the buckled strings that can result in Utotal/Ubending < 1 and thus Q < Q029,44. Included in Fig. 4c, we compare experimental data to numerical results from reduced-order modeling of FE simulations that offer a higher θ resolution than experiments. The simulations (dots) are in good agreement with our experiments, except around θ = −0.1 rad, where we faced numerical instability. We shall note that the simulations were conducted after a careful convergence study of the corresponding multi-mode reduced-order model, which also accounted for the broken symmetry induced by buckling. Furthermore, all the simulated results in Fig. 4 were obtained using the same material and geometric properties, and θ is the only parameter that is changed. The agreement between experiments and simulations is evident in the reliability of the modeling approach.

Fig. 4: Tuning of the resonance frequency, Q-factor, and the effective geometric nonlinearity by varying support angles.
figure 4

Finite element (FE) based (dots) and measured (diamonds) a resonance frequencies f0, b Q-factors, and c geometric nonlinearity β of the first symmetric out-of-plane mode; note that the error bars are smaller than the data marker size. The presented resonators share L = 200 μm, w = 2 μm, ws = 1 μm, and Ls = 150 μm but different support angles θ, which vary from − 0.5 rad to 0.5 rad. The height of string center H is shown in (c) for negative θ.

Full size image

Discussion

The large tunability of the dynamical properties provided by geometric design offers new possibilities for engineering devices that are linear over a large range. For instance, as shown in Fig. 4, for devices with θ > 0, it is possible to substantially decrease the geometric nonlinearity while maintaining a high Q-factor and thus expand the linear dynamic range of the resonator by increasing the onset of nonlinearity \({a}_{{{{{{{{\rm{1dB}}}}}}}}}\propto {Q}^{-\frac{1}{2}}{\alpha }^{\frac{1}{2}}{\beta }^{-\frac{1}{2}}\)45,46. As an example, in devices with support angle θ around 0.1 rad, we note that both a1dB and Q-factor are three times higher than in conventional doubly clamped strings (see Supplementary Note 4). The methodology also provides the possibility to minimize β, such that it is infinitely small. We expect this condition to occur around θ = −0.0063 rad according to FE simulations, where a transition from hardening to softening is observed. This condition is associated with the onset of buckling bifurcation when θ is varied, which is unstable and challenging to control, yet potentially possible to be stabilized by using external forces1. On the other hand, it is desirable to operate devices in the nonlinear regime by squeezing the dynamic range (DR) of nanomechanical resonators (defined as the ratio DR = 20lg(a1dB/ath), where ath is the thermo-mechanical noise floor). Notably, in our current designs, we observe a reduction of DR from 64dB in a double-clamped string to 51dB in devices with slender supports. We foresee that by further reducing the thickness h = 340 nm of our devices, DR can even be squeezed down to near zero (see Supplementary Note 4), enabling the study of nonlinear dynamics in the Brownian limit47. Moreover, recent studies have shown that by enhancing the ratio β/α, it is possible to realize nonlinear nanomechanical resonators that approach the quantum ground state in doubly clamped carbon nanotubes48. Our simulations suggest that the absolute ratio of β/α in buckled resonators can be increased to an order of magnitude higher than that of double-clamped ones, suggesting that buckling can be used as an effective tool for increasing nonlinearities in modes that operate close to the quantum ground state49.

In conclusion, we introduce a strategy to engineer the nonlinear dynamics of nanomechanical devices effectively by geometric design. We discuss the role of in-plane to out-of-plane coupling induced by the soft-clamping support, which we used to tune the geometric nonlinearity of high-Q string resonators down by two orders of magnitude compared to double-clamped strings. We also show that by carefully tuning the support angle, we can tune the stress in Si3N4 strings over the whole range from highly tensile to compressive, thus inducing a buckling bifurcation. We found that the post-buckled configurations show a softening nonlinearity with an absolute Duffing constant that increases with the buckling amplitude. The advantage of the presented nonlinear tuning method is that it does not require electrical or opto-thermal forces, providing a purely mechanical method for controlling nonlinear dynamics by design. This design is facilitated by the presented and validated analytic and finite element modeling techniques. Thus, the methodology allows for optimizing the dynamic range for each mechanical string resonator in the design phase, paving a robust way to realize large arrays of nanomechanical resonators with increased linearity or varying nonlinear functionalities for applications such as sensor arrays and quantum processing.

Methods

Sample fabrication

Devices are fabricated by electron beam lithography and reactive ion etching from high-stress Si3N4 layers. The layers are grown by low-pressure chemical vapor deposition (LPCVD) with a thickness of h = 340 nm on a silicon substrate. The devices are then suspended by a fluorine-based (SF6) deep reactive ion underetching step27. All nanomechanical resonators studied in this work are made of Si3N4 deposited on the same wafer, which guarantees almost identical mechanical properties, with an initial isotropic stress σ0 = 1.06 GPa, Young’s modulus E = 271 GPa, Poisson’s ratio ν = 0.23, mass density ρ = 3100 kg/m3 and the intrinsic quality factor Q0 = 9864. The resonance frequency f0 and Q-factor characterization of these devices are detailed out in our previous work29.