Controlling Resonator Nonlinearities and Modes through Geometry Optimization

Controlling the nonlinearities of MEMS resonators is critical for their successful implementation in a wide range of sensing, signal conditioning, and filtering applications. Here, we utilize a passive technique based on geometry optimization to control the nonlinearities and the dynamical response of MEMS resonators. Also, we explored active technique i.e., tuning the axial stress of the resonator. To achieve this, we propose a new hybrid shape combining a straight and initially curved microbeam. The Galerkin method is employed to solve the beam equation and study the effect of the different design parameters on the ratios of the frequencies and the nonlinearities of the structure. We show by adequately selecting the parameters of the structure; we can realize systems with strong quadratic or cubic effective nonlinearities. Also, we investigate the resonator shape effect on symmetry breaking and study different linear coupling phenomena: crossing, veering, and mode hybridization. We demonstrate the possibility of tuning the frequencies of the different modes of vibrations to achieve commensurate ratios necessary for activating internal resonance. The proposed method is simple in principle, easy to fabricate, and offers a wide range of controllability on the sensor nonlinearities and response.


Introduction
The Resonant microstructures have demonstrated their potential in sensing [1][2][3], filtering [4,5], logic devices [6,7], energy harvesting [8,9], and signal processing and also timing applications [10,11]. Traditionally, many of these applications require operating the resonator in the linear regime. However, as the dimensions of the resonator are shrunk, their response becomes nonlinear due to the size and surface effects and quantum forces. These nonlinearities can be attributed to the midplane stretching effect, material nonlinearity, and the nonlinear nature of the actuation force [12,13]. Although at the nano and micro-scale these nonlinearities are inherently presented, one can introduce nonlinearities at the macro scale by using materials that have a nonlinear stress-strain relation or imposing certain geometrical constraints or flexible boundary conditions. Usually, nonlinear behavior is considered undesirable and might lead to the failure of the device. Yet, advances in the field of nonlinear dynamics led to the creation of a new class of devices that make use of the nonlinear phenomena in applications such as imaging [14][15][16], ultra-sensitive mass sensor [3], energy harvesting [17], filtering [18] and frequency divider [19].
Recently, researchers investigated different methods for tailoring nonlinearities to achieve targeted responses. These methods can be classified into two categories: active and passive techniques. In active techniques, the nonlinearities can be controlled over the lifetime of the device using electrostatic [20,21], electrothermal [22,23], and electromagnetic actuation [24]. Active methods are implemented to overcome fabrication imperfection and impedance mismatching, but they suffer from the need for additional circuits, high power consumption, and increases the device size and cost. On the other hand, passive originated from midplane stretching), Figure 1, to induce a geometrical nonlinearity that can be tuned to achieve the targeted nonlinear phenomenon. The presented concept is simple in principle, scalable (macro, micro, and nano scales), and can be employed to achieve a wide range of targeted responses.

Materials and Methods
We propose a new structure that combines the straight and buckled shape in the same design to take advantage of both configurations on the resonator behavior. We introduce a buckled shape dome of length DL × Lhb in the beam configuration, Figure 1. The initial shape ̂ℎ ,0 (̂) of the clamped-clamped hybrid beam, of length Lhb, is described by Equation (1).
Here ̂ denotes the position along the beam length, and ̂0 represents the dome midpoint rise.
[̂] denotes the unit-step function. The cross-sectional area of the hybrid beam is assumed to be rectangular = ℎ ℎ ℎ with a moment of inertia, = ℎ ℎ ℎ 3 /12 where bhb and hhb denote the width and the thickness, respectively.
As a case study, we consider a microbeam made from silicon and with dimensions given in Table 1. The equation of motion governing the transverse vibration of the hybrid beam is given by [23]: The resonating beam is subjected to the following clamped-clamped boundary conditions:

Materials and Methods
We propose a new structure that combines the straight and buckled shape in the same design to take advantage of both configurations on the resonator behavior. We introduce a buckled shape dome of length D L × L hb in the beam configuration, Figure 1. The initial shapeŵ hb,0 (x) of the clamped-clamped hybrid beam, of length L hb , is described by Equation (1).ŵ Herex denotes the position along the beam length, andb 0 represents the dome midpoint rise. U[x] denotes the unit-step function. The cross-sectional area of the hybrid beam is assumed to be rectangular A = b hb h hb with a moment of inertia, I = b hb h 3 hb /12 where b hb and h hb denote the width and the thickness, respectively.
As a case study, we consider a microbeam made from silicon and with dimensions given in Table 1. The equation of motion governing the transverse vibration of the hybrid beam is given by [23]: The resonating beam is subjected to the following clamped-clamped boundary conditions:ŵ hb 0,t =ŵ hb L hb ,t = 0; and ∂ŵ hb ∂x (0,t) Equations (2) and (3) are normalized by using the following nondimensional variables: where r = I A denotes the gyration radius of the cross-section, and T s = ρb hb h hb L 4 hb EI is the timescale. Substituting Equation (4) into Equations (2) and (3), we obtain the nondimensionalized equation of motion and associated boundary conditions: The nondimensional parameters appearing in Equation (5) are defined as follow: To study the variation of the natural frequencies of the hybrid beam due to the alteration of the axial load, N non , we start by linearizing the equation of motion describing the small dynamic behavior of the arch beam around the new static configuration governed by Equation (9)  with the associated boundary conditions w hb,s (0) = w hb,s (1) = 0 and dw hb,s dx x=0 = dw hb,s dx x=1 = 0 (9) Next, after dropping the damping terms from the equation of motion, Equation (5), the total deflection is assumed to be equal to the sum of the static configuration, w hb,s (x), and a small dynamic deflection of the hybrid beam, w hb,d (x,t), around w hb,s (x). The outcome equation becomes with the associated boundary conditions.
The eigenvalue problem is solved using the Galerkin discretization. Toward this, we let where q i (t) (i = 1, 2..., n) denotes the nondimensional modal coordinates and ϕ i (x) (i = 1, 2..., n) denotes the mode shape of the unforced and undamped clamped-clamped standard arch beam [16].
Next, we substitute Equation (12) into Equation (10), multiply the outcome by the mode shape ϕ j , and integrate over the beam domain, which yields the following equation: ..
Using 5 modes, the static deflection and then the Jacobian of the five obtained equations are computed for each N non . Then, we calculate the resonators' natural frequencies, at constant N non , by taking these eigenvalues' square root. A convergence analysis was conducted in previous works to determine the number of exact modes needed to simulate the curved resonator behavior [23,49].

Results
To explore the potential of the proposed hybrid design in controlling the nonlinearities and the microbeam frequencies, we solve the beam equation for a different set of dome length, initial rise, and axial stress. The geometric and mechanical parameters of the hybrid beam used in the simulations are given in Table 1. To calculate the static deflection and eigenvalue, we used five modes of the arch configuration needed to reach convergence.
Firstly, to analyze the effect of the different geometry parameters on the nonlinearity's type and magnitude, we compute the quadratic and cubic nonlinearity coefficients of the proposed structure (the derivation is presented in the Appendix A). As shown in Figure 2a-c, the quadratic nonlinearity is more sensitive to the dome length and changes rapidly for a fixed initial rise and thickness. Also, increasing the initial rise expands the quadratic nonlinearity range with minimal effect on the cubic nonlinearity. Figure 2d,e show that increasing the thickness from 1 µm to 3 µm, for an initial rise equals to 2 µm, reduces the cubic nonlinearity by order of magnitude. As revealed in Figure 2f, designing a microstructure with strong quadratic, cubic, or even zero effective nonlinearity can be easily achieved by appropriately selecting the dome length, beam thickness, and initial rise. One should mention that the dynamic response of beams with strong quadratic nonlinearity will experience a softening effect while it will show a hardening effect for beams with strong cubic nonlinearities. where qi(t) (i = 1,2..n) denotes the nondimensional modal coordinates and φi(x) (i = 1,2..n) denotes the mode shape of the unforced and undamped clamped-clamped standard arch beam [16].
Next, we substitute Equation (12) into Equation (10), multiply the outcome by the mode shape φj, and integrate over the beam domain, which yields the following equation: Using 5 modes, the static deflection and then the Jacobian of the five obtained equations are computed for each Nnon. Then, we calculate the resonators' natural frequencies, at constant Nnon, by taking these eigenvalues' square root. A convergence analysis was conducted in previous works to determine the number of exact modes needed to simulate the curved resonator behavior [23,49].

Results
To explore the potential of the proposed hybrid design in controlling the nonlinearities and the microbeam frequencies, we solve the beam equation for a different set of dome length, initial rise, and axial stress. The geometric and mechanical parameters of the hybrid beam used in the simulations are given in Table 1. To calculate the static deflection and eigenvalue, we used five modes of the arch configuration needed to reach convergence.
Firstly, to analyze the effect of the different geometry parameters on the nonlinearity's type and magnitude, we compute the quadratic and cubic nonlinearity coefficients of the proposed structure (the derivation is presented in the Appendix A). As shown in Figure 2a-c, the quadratic nonlinearity is more sensitive to the dome length and changes rapidly for a fixed initial rise and thickness. Also, increasing the initial rise expands the quadratic nonlinearity range with minimal effect on the cubic nonlinearity. Figure 2d,e show that increasing the thickness from 1 µm to 3 µm, for an initial rise equals to 2 µm, reduces the cubic nonlinearity by order of magnitude. As revealed in Figure 2f, designing a microstructure with strong quadratic, cubic, or even zero effective nonlinearity can be easily achieved by appropriately selecting the dome length, beam thickness, and initial rise. One should mention that the dynamic response of beams with strong quadratic nonlinearity will experience a softening effect while it will show a hardening effect for beams with strong cubic nonlinearities. Next, we investigate the influence of the dome length on the first three modes of the resonator. As shown in Figure 3a-c, the frequencies vary nonlinearly with the dome length. Also, one can note that for particular dome length values, the initial rise does not affect the frequency values; for example, at DL around 0.5, the first mode remains constant for different initial rise values. The same can be noted for the second mode at DL around 0.3 and 1 and the third mode at DL around 0.25 and 0.7. These results suggest that the geometry parameters have different effects on the effective modal stiffness, allowing a wide range of tunability.
The ratios between different vibration modes are depicted in Figure 3d-f for different initial rises of the dome. Different commensurate ratios (2, 3..., 7) can be obtained between different modes as tuning the dome length. A commensurate ratio between two modes is necessary to activate the nonlinear energy exchange among the involved modes via internal resonance. Figure 3e,f show that for a particular dome length, the ratios between the three lowest frequencies are commensurate (ω3/ω2 = ω2/ω1 = 2), leading to a potential combined internal resonance and more complex nonlinear energy exchange between involved modes. Next, we investigate the influence of the dome length on the first three modes of the resonator. As shown in Figure 3a-c, the frequencies vary nonlinearly with the dome length. Also, one can note that for particular dome length values, the initial rise does not affect the frequency values; for example, at D L around 0.5, the first mode remains constant for different initial rise values. The same can be noted for the second mode at D L around 0.3 and 1 and the third mode at D L around 0.25 and 0.7. These results suggest that the geometry parameters have different effects on the effective modal stiffness, allowing a wide range of tunability.
The ratios between different vibration modes are depicted in Figure 3d-f for different initial rises of the dome. Different commensurate ratios (2, 3..., 7) can be obtained between different modes as tuning the dome length. A commensurate ratio between two modes is necessary to activate the nonlinear energy exchange among the involved modes via internal resonance. Figure 3e,f show that for a particular dome length, the ratios between the three lowest frequencies are commensurate (ω 3 /ω 2 = ω 2 /ω 1 = 2), leading to a potential combined internal resonance and more complex nonlinear energy exchange between involved modes. Next, we investigate the influence of the dome length on the first three modes of the resonator. As shown in Figure 3a-c, the frequencies vary nonlinearly with the dome length. Also, one can note that for particular dome length values, the initial rise does not affect the frequency values; for example, at DL around 0.5, the first mode remains constant for different initial rise values. The same can be noted for the second mode at DL around 0.3 and 1 and the third mode at DL around 0.25 and 0.7. These results suggest that the geometry parameters have different effects on the effective modal stiffness, allowing a wide range of tunability.
The ratios between different vibration modes are depicted in Figure 3d-f for different initial rises of the dome. Different commensurate ratios (2, 3..., 7) can be obtained between different modes as tuning the dome length. A commensurate ratio between two modes is necessary to activate the nonlinear energy exchange among the involved modes via internal resonance. Figure 3e,f show that for a particular dome length, the ratios between the three lowest frequencies are commensurate (ω3/ω2 = ω2/ω1 = 2), leading to a potential combined internal resonance and more complex nonlinear energy exchange between involved modes. As tuning the dome length, different commensurate ratios between modes could be achieved. The above is a necessary condition for activating the nonlinear coupling between different modes via internal resonance. The DL axis starts from 0.1. The red crosses are associate with the dome length at which we have commensurate ratios.
The initial static profile is another parameter of great interest that affects the response dynamics. Here, we study the effect of the axial stress and the dome length on the static deflection by solving Equation (8). For this analysis, the case of = 2 μm is chosen. Figure  4 shows that the dome length significantly influences the hybrid resonator's static shape as tuning its axial load. For dome length equal to half beam length, the asymmetry is apparent by looking at the displacement at the midpoint, the quarter (the dome-side), and the third-quarter (the straight-side) of the beam. For dome lengths lower than half of the beam length, the static shapes combine a buckled shape on the straight side and the arch shapes on the other side. As the dome length exceeds the half beam length, the static deflection trend is similar to the normal arch. The static profile is important in determining the symmetry and asymmetry of the system (symmetry with respect to the midpoint plane). For asymmetric shapes, the classical crossing between symmetric and antisymmetric modes will be affected and changed to avoid-crossing between both modes.  (f)b 0 = 3 µm. As tuning the dome length, different commensurate ratios between modes could be achieved. The above is a necessary condition for activating the nonlinear coupling between different modes via internal resonance. The D L axis starts from 0.1. The red crosses are associate with the dome length at which we have commensurate ratios.
The initial static profile is another parameter of great interest that affects the response dynamics. Here, we study the effect of the axial stress and the dome length on the static deflection by solving Equation (8). For this analysis, the case ofb 0 = 2 µm is chosen. Figure 4 shows that the dome length significantly influences the hybrid resonator's static shape as tuning its axial load. For dome length equal to half beam length, the asymmetry is apparent by looking at the displacement at the midpoint, the quarter (the dome-side), and the third-quarter (the straight-side) of the beam. For dome lengths lower than half of the beam length, the static shapes combine a buckled shape on the straight side and the arch shapes on the other side. As the dome length exceeds the half beam length, the static deflection trend is similar to the normal arch. The static profile is important in determining the symmetry and asymmetry of the system (symmetry with respect to the midpoint plane). For asymmetric shapes, the classical crossing between symmetric and antisymmetric modes will be affected and changed to avoid-crossing between both modes.  The above is a necessary condition for activating the nonlinear coupling between different modes via internal resonance.
The DL axis starts from 0.1. The red crosses are associate with the dome length at which we have commensurate ratios.
The initial static profile is another parameter of great interest that affects the response dynamics. Here, we study the effect of the axial stress and the dome length on the static deflection by solving Equation (8). For this analysis, the case of ̂0 = 2 µm is chosen. Figure  4 shows that the dome length significantly influences the hybrid resonator's static shape as tuning its axial load. For dome length equal to half beam length, the asymmetry is apparent by looking at the displacement at the midpoint, the quarter (the dome-side), and the third-quarter (the straight-side) of the beam. For dome lengths lower than half of the beam length, the static shapes combine a buckled shape on the straight side and the arch shapes on the other side. As the dome length exceeds the half beam length, the static deflection trend is similar to the normal arch. The static profile is important in determining the symmetry and asymmetry of the system (symmetry with respect to the midpoint plane). For asymmetric shapes, the classical crossing between symmetric and antisymmetric modes will be affected and changed to avoid-crossing between both modes.   Another parameter that we explore its effect on the frequency values is the axial stress (compressive load). The axial load could be controlled by different transduction mechanisms, such as the electrothermal voltage [22,23] and sided electrostatic electrode [41]. Beams with tunable axial load show good capabilities in several potential applications, mainly thermal-conductivity-based gas and pressure sensors, force sensors, logic memories, and temperature sensors [7,50]. Figures 5-8 show the variation of the lowest three natural frequencies as tuning the axial load (compressive load) for various dome lengths and initial rises. Figure 5 shows the variation of the lowest three frequencies with axial load for the classical cases: straight beam (Figure 5a) and standard arch beam ( Figure 5b). As classically known, the frequencies of straight beam decrease at higher axial load values until reaching buckling. After buckling, the frequencies of the first and second symmetric modes increase with axial load, while the first antisymmetric mode remains constant. This leads to a crossing between the first symmetric and the first antisymmetric mode. As shown in the insets of Figure 5a, both modes do not interact with each other at the crossing zone; hence no mode hybridization is reported. The same phenomenon is reported for the arch beam case when the first symmetric and antisymmetric modes cross. One should note that the buckling instability disappears when introducing an initial curvature to the structure, Figure 5b. Also, the antisymmetric modes for the initially curved structures decrease with axial loads for the whole range of compressive axial loads contrary to the straight beam configuration.  Another parameter that we explore its effect on the frequency values is the axial stress (compressive load). The axial load could be controlled by different transduction mechanisms, such as the electrothermal voltage [22,23] and sided electrostatic electrode [41]. Beams with tunable axial load show good capabilities in several potential applications, mainly thermal-conductivity-based gas and pressure sensors, force sensors, logic memories, and temperature sensors [7,50]. Figures 5-8 show the variation of the lowest three natural frequencies as tuning the axial load (compressive load) for various dome lengths and initial rises. Figure 5 shows the variation of the lowest three frequencies with axial load for the classical cases: straight beam (Figure 5a) and standard arch beam ( Figure 5b). As classically known, the frequencies of straight beam decrease at higher axial load values until reaching buckling. After buckling, the frequencies of the first and second symmetric modes increase with axial load, while the first antisymmetric mode remains constant. This leads to a crossing between the first symmetric and the first antisymmetric mode. As shown in the insets of Figure 5a, both modes do not interact with each other at the crossing zone; hence no mode hybridization is reported. The same phenomenon is reported for the arch beam case when the first symmetric and antisymmetric modes cross. One should note that the buckling instability disappears when introducing an initial curvature to the structure, Figure 5b. Also, the antisymmetric modes for the initially curved structures decrease with axial loads for the whole range of compressive axial loads contrary to the straight beam configuration.  Another parameter that we explore its effect on the frequency values is the axial stress (compressive load). The axial load could be controlled by different transduction mechanisms, such as the electrothermal voltage [22,23] and sided electrostatic electrode [41]. Beams with tunable axial load show good capabilities in several potential applications, mainly thermal-conductivity-based gas and pressure sensors, force sensors, logic memories, and temperature sensors [7,50]. Figures 5-8 show the variation of the lowest three natural frequencies as tuning the axial load (compressive load) for various dome lengths and initial rises. Figure 5 shows the variation of the lowest three frequencies with axial load for the classical cases: straight beam ( Figure 5a) and standard arch beam ( Figure 5b). As classically known, the frequencies of straight beam decrease at higher axial load values until reaching buckling. After buckling, the frequencies of the first and second symmetric modes increase with axial load, while the first antisymmetric mode remains constant. This leads to a crossing between the first symmetric and the first antisymmetric mode. As shown in the insets of Figure 5a, both modes do not interact with each other at the crossing zone; hence no mode hybridization is reported. The same phenomenon is reported for the arch beam case when the first symmetric and antisymmetric modes cross. One should note that the buckling instability disappears when introducing an initial curvature to the structure, Figure 5b. Also, the antisymmetric modes for the initially curved structures decrease with axial loads for the whole range of compressive axial loads contrary to the straight beam configuration.   Next, we investigate the effect of the dome length and axial load on the frequency values of the first three modes for initial rise b 0 = 1 µm, as depicted in Figures 6 and A1. The results show that the first two modes are more sensitive to the dome length compared with the third mode. The third mode seems to have the same qualitative behavior as increasing the axial load without interaction with the other modes. The two lowest frequencies initially decrease with axial load, then increase after reaching a threshold axial load value. These two frequencies should cross without interacting at a certain axial load for a normal arch or buckled beam, as shown in Figure 5. However, the presence of the dome creates an asymmetry that transforms the crossing into avoided crossing between both modes. The avoided crossing frequency width varies for different dome lengths.
In contrast to the arch beam or buckled beam, the first two lowest modes hybridize around the avoid-crossing zone. The hybridization of both first symmetric and antisymmetric modes leads to the same out-of-phase mode shapes at the avoid-crossing zone. The above is similar to the case of an arch beam under asymmetric electrostatic actuation [44]. The first two lowest modes show an avoided crossing for low values of dome lengths that are smaller than 0.5 L hb . The inset of Figure 6d shows the mode shapes at different axial load values around the avoided-crossing region for dome length D L = 0.3; the two modes hybridize, and then interchange mode shapes similar to the veering phenomena. For dome lengths equal to 0.5 L hb , Figure 6e, the two modes show a similar behavior but with a higher frequency difference between them. Even if the frequency difference is high, an alteration of mode shapes is reported for this case, as seen in Figure 6e. Hence, the presence of the dome creates an asymmetry in the beam that leads to the veering between different modes as tuning their stiffness without crossing. Next, we investigate the effect of the dome length and axial load on the frequency values of the first three modes for initial rise b0 = 1 µm, as depicted in Figures 6 and A1. The results show that the first two modes are more sensitive to the dome length compared with the third mode. The third mode seems to have the same qualitative behavior as increasing the axial load without interaction with the other modes. The two lowest frequencies initially decrease with axial load, then increase after reaching a threshold axial load value. These two frequencies should cross without interacting at a certain axial load for a normal arch or buckled beam, as shown in Figure 5. However, the presence of the dome creates an asymmetry that transforms the crossing into avoided crossing between both modes. The avoided crossing frequency width varies for different dome lengths.
In contrast to the arch beam or buckled beam, the first two lowest modes hybridize around the avoid-crossing zone. The hybridization of both first symmetric and antisymmetric modes leads to the same out-of-phase mode shapes at the avoid-crossing zone. The above is similar to the case of an arch beam under asymmetric electrostatic actuation [44]. The first two lowest modes show an avoided crossing for low values of dome lengths that are smaller than 0.5 Lhb. The inset of Figure 6d shows the mode shapes at different axial load values around the avoided-crossing region for dome length DL = 0.3; the two modes hybridize, and then interchange mode shapes similar to the veering phenomena. For dome lengths equal to 0.5 Lhb, Figure 6e, the two modes show a similar behavior but with a higher frequency difference between them. Even if the frequency difference is high, an alteration of mode shapes is reported for this case, as seen in Figure 6e. Hence, the presence of the dome creates an asymmetry in the beam that leads to the veering between different modes as tuning their stiffness without crossing. For initial riseb 0 = 2 µm, the frequency tuning shows more mode coupling as varying the dome length and the axial load, as shown in Figures 7 and A2. Figure 7a shows, for the case of D L = 0.5, the ability to have veering between the second and third modes and crossing between the first and second modes. As shown in the insets of Figure 7a, the mode shapes show the hybridization of the two modes around the veering zone. The interchange between the two lowest modes is also demonstrated for lower axial load even though the difference between the two frequencies is high. Another feature that could be seen for dome length D L = 0.8, the two lowest frequencies vary with an almost constant ratio for a wide range of axial load, as shown in Figure 7b. The results show the potential of utilizing the proposed resonators as stable resonators in a harsh environment or for temperature compensation through internal resonance. The energy transfer among the coupled modes acts as a stabilizing feedback loop that suppresses the amplitude and frequency fluctuations [10].
As tuning the geometric parameters and the dome length, we showed that the type and magnitude of nonlinearity could be controlled. The ratio can also be tuned to a commensurate number. The above would be necessary to activate the nonlinear coupling between the modes via internal resonance for a wide range of axial load [46]. Operating a resonator in the internal resonance zone increases the resonator's frequency stability due to nonlinear energy transfer between the two involved modes [10,51]. The proposed system shows the potential to be operated in an internal resonance zone for a wide range of axial loads values, which mimics temperature fluctuation. For initial rise ̂0 = 2 µm, the frequency tuning shows more mode coupling as varying the dome length and the axial load, as shown in Figures 7 and A2. Figure 7a shows, for the case of DL = 0.5, the ability to have veering between the second and third modes and crossing between the first and second modes. As shown in the insets of Figure 7a, the mode shapes show the hybridization of the two modes around the veering zone. The interchange between the two lowest modes is also demonstrated for lower axial load even though the difference between the two frequencies is high. Another feature that could be seen for dome length DL = 0.8, the two lowest frequencies vary with an almost constant ratio for a wide range of axial load, as shown in Figure 7b. The results show the potential of utilizing the proposed resonators as stable resonators in a harsh environment or for temperature compensation through internal resonance. The energy transfer among the coupled modes acts as a stabilizing feedback loop that suppresses the amplitude and frequency fluctuations [10].
As tuning the geometric parameters and the dome length, we showed that the type and magnitude of nonlinearity could be controlled. The ratio can also be tuned to a commensurate number. The above would be necessary to activate the nonlinear coupling between the modes via internal resonance for a wide range of axial load [46]. Operating a resonator in the internal resonance zone increases the resonator's frequency stability due to nonlinear energy transfer between the two involved modes [10,51]. The proposed system shows the potential to be operated in an internal resonance zone for a wide range of axial loads values, which mimics temperature fluctuation. Increasing the initial rise to ̂0 = 3 µm will change the coupling behavior among modes for different dome lengths. However, the same features could be seen such as veering (avoid-crossing) and constant frequency ratios, as depicted in Figures 8 and A3. Different frequency tuning and coupling phenomena could be demonstrated by carefully choosing the initial rise and the dome length. This depends on the targeted applications. The presented approach here is scalable, so some features could be achieved for higher or lower frequency ranges by enlarging or shrinking the beam sizes. Increasing the initial rise tob 0 = 3 µm will change the coupling behavior among modes for different dome lengths. However, the same features could be seen such as veering (avoid-crossing) and constant frequency ratios, as depicted in Figures 8 and A3. Different frequency tuning and coupling phenomena could be demonstrated by carefully choosing the initial rise and the dome length. This depends on the targeted applications. The presented approach here is scalable, so some features could be achieved for higher or lower frequency ranges by enlarging or shrinking the beam sizes. Increasing the initial rise to 0 = 3 μm will change the coupling behavior among modes for different dome lengths. However, the same features could be seen such as veering (avoid-crossing) and constant frequency ratios, as depicted in Figures 8 and A3. Different frequency tuning and coupling phenomena could be demonstrated by carefully choosing the initial rise and the dome length. This depends on the targeted applications. The presented approach here is scalable, so some features could be achieved for higher or lower frequency ranges by enlarging or shrinking the beam sizes. Next, a 3D multi-physics finite-element model, using the commercial finite element software COMSOL [51], is conducted to validate the numerical model presented in this work and verify that the various assumptions made in the analytical model are valid. The hybrid beam is simulated using the parameters given in Table 1. The initial rise and dome length are equal to 2 μm and 0.8 Lhb, respectively. In order to simulate the effect of axial load, electrothermal actuation was used [23]. Electrothermal tuning of MEMS resonators Next, a 3D multi-physics finite-element model, using the commercial finite element software COMSOL [51], is conducted to validate the numerical model presented in this work and verify that the various assumptions made in the analytical model are valid. The hybrid beam is simulated using the parameters given in Table 1. The initial rise and dome length are equal to 2 µm and 0.8 L hb , respectively. In order to simulate the effect of axial load, electrothermal actuation was used [23]. Electrothermal tuning of MEMS resonators was deeply investigated in the literature in the case of the bridge [21], arch [22], V-shaped [27], and U-shaped [52] resonators, and were demonstrated for several potential applications, such as filtering, logic memories, and pressure/gas sensing. The Solid Mechanics, Electric Currents, and Heat Transfer interfaces are considered to account for the various physical domains of the problem as described in [23]. One should mention that the selection of the electrothermal tuning is based on its popularity and simplicity of this method in tuning the axial stress of microresonators. The relation between the nondimensional induced stress (N non,TH = L 2 hb EIŜ TH ) and the electrothermal voltage is plotted in Figure 9a, whereŜ TH is the thermal induced axial stress extracted from COMSOL. Figure 9b reveals the variation of the three lowest natural frequencies of the hybrid resonator with the nondimensional axial load (induced by the corresponding electrothermal voltage) showing a good agreement with the numerical data presented in Figure 7b.
Micromachines 2021, 12, x FOR PEER REVIEW 12 of 17 was deeply investigated in the literature in the case of the bridge [21], arch [22], V-shaped [27], and U-shaped [52] resonators, and were demonstrated for several potential applications, such as filtering, logic memories, and pressure/gas sensing. The Solid Mechanics, Electric Currents, and Heat Transfer interfaces are considered to account for the various physical domains of the problem as described in [23]. One should mention that the selection of the electrothermal tuning is based on its popularity and simplicity of this method in tuning the axial stress of microresonators. The relation between the nondimensional induced stress ( , = ) and the electrothermal voltage is plotted in Figure 9a, where is the thermal induced axial stress extracted from COMSOL. Figure 9b reveals the variation of the three lowest natural frequencies of the hybrid resonator with the nondimensional axial load (induced by the corresponding electrothermal voltage) showing a good agreement with the numerical data presented in Figure 7b.

Discussion
This study proposes a new hybrid structure combining the straight and arch beam shapes. The proposed design is simple in principle and easy to fabricate using the standard two-mask microfabrication process of silicon-on-insulator (SOI) wafers, more details can be found in [17]. A passive and active systematic methodology is demonstrated to tailor the system nonlinearities by controlling the initial shaped and geometric parameters. The results demonstrate the possibility of realizing structures with strong cubic, quadratic, or even zero effective nonlinearity. The proposed design shows different ways of tuning the frequencies by adequately selecting the dome length and initial rise values. Avoid-crossing and mode-hybridization between different modes are demonstrated for various geometric parameters. Depending on the targeted applications and performance metrics, certain geometric parameters can be selected. For example, Figure 8c shows the possibility of designing a structure with a stable second mode frequency over a wide range of axial stress. This characteristic is vital for timing and sensing devices operated in unstable temperature conditions. Also, we show the possibility of achieving commensurate ratios necessary for activating nonlinear energy transfer between the different modes. Figure 3d shows the possibility of realizing a structure with a 2:1 frequency ratio over a wide range of dome lengths, which allows large tolerance in the fabrication process. Moreover, our results show the possibility of realizing a structure with a 2:1 frequency ratio over a wide range of axial stress, inset of Figure 7b, such structure has a great potential in energy harvesting applications operated in unstable temperature conditions. One should note that more studies need to be done to characterize the nonlinear dynamics of the proposed device under different actuation forces and different axial load configurations. Driving the proposed designs using electrostatic force is expected to have rich and complex linear and nonlinear dynamics phenomena; Snap-through, pull-in bifurcations, switching between softening and hardening. Also, using partial electrodes at different positions with respect to the dome length will influence the device static and dynamic behavior. These complex and rich dynamics can be exploited to enhance the performance of the resonators or avoid it for safe implementation of the resonator.
For simplicity, we introduce S, Sq, Sc, denote the geometric nonlinear coefficient, its quadratic order, and its cubic order. Sq and Sc are defined as below: The functions χ q and χ c are defined as: The function ψ is obtained by solving the boundary value problem given by where H and Λ are defined by

Appendix B
The variation of the three lowest natural frequencies with axial load for different dome length D L and for different initial risesb 0 .
S, Sq, Sc, denote the geometric nonlinear coefficient, its quadratic order, and its cubic order. Sq and Sc are defined as below: The functions χq and χc are defined as: The function ψ is obtained by solving the boundary value problem given by

Appendix B
The variation of the three lowest natural frequencies with axial load for different dome length DL and for different initial rises .