Three-to-One Internal Resonance in MEMS Arch Resonators

We present an investigation of the nonlinear dynamics of a microelectromechanical system (MEMS) arch subjected to a combination of AC and DC loadings in the presence of three-to-one internal resonance. The axial force resulting from the residual stress or temperature variation is considered in the governing equation of motion. The method of multiple scales is used to solve the governing equation. A four first-order ordinary differential equation describing the modulation of the amplitudes and phase angles is obtained. The equilibrium solution and its stability of the modulation equations are determined. Moreover, we also obtain the reduced-order model (ROM) of the MEMS arch employing the Galerkin scheme. The dynamic response is presented in the form of time traces, Fourier spectrum, phase-plane portrait, and Poincare sections. The results show that when there is an internal resonance, the energy transfer occurs between the first and third modes. In addition, the response of the MEMS arch presents abundant dynamic behaviors, such as Hopf bifurcation and quasiperiodic motions.


Introduction
Microelectromechanical system (MEMS) devices are widely used in various engineering applications due to their low energy consumption, high sensitivities, and small size. Examples of these include resonators [1], bio-sensors [2], actuators [3], airbag accelerometers [4], and so on. In the design of these MEMS devices, various structures are used, such as fixed-fixed microbeams [5], cantilever microbeams [6], H-shaped microbeams [7], double microbeams [8], microplates [9] and so on. Due to their attractive features, arch structures have attracted a growing amount of attention from the microelectromechanical system community. For example, it has been recently reported that arch structures are used as mechanical memories [10] and logic devices [11]. However, because of the nonlinear nature aroused by its own initial curvature, mid-plane stretching, and electric load, MEMS arches can present abundant nonlinear dynamic behaviors.
There has been much theoretical research carried out on the nonlinear dynamics of MEMS arches. For example, Ouakad et al. [12] investigated the primary resonance of the first natural frequency of a MEMS arch using the method of multiple scales. They found that the arch showed softening behavior. Younis et al. [13] proposed a reduced-order model of a MEMS arch based on the Galerkin method and then examined its primary resonance and subharmonic resonance. Krylov et al. [14], using same technique, investigated the pull-in and snap-through instability of MEMS arches. Farokhi et al. [15] constructed a high-dimensional reduced-order Galerkin model and then studied the effect of various system parameters on pull-in instabilities. Laura et al. [16,17]  published. For example, Ahmad et al. [18] experimentally examined the primary resonance of the first natural frequency of an MEMS arch. Nonlinear phenomena such as jumps, hysteresis, and softening behaviors were observed. Laura et al. [19] detected super harmonic resonances in frequency sweeps experiments. Younis et al. [20] investigated the response of a MEMS arch under a large AC harmonic load experimentally. They found that the dynamic snap-through and dynamic pull-in can be induced for some specific parameters. Zhang et al. [21] examined the snap-through instability and pull-in instability of a MEMS arch with different configurations. Hajjaj et al. [22] examined the effect of initial shapes, arc, and cosine wave on the dynamic behaviors of MEMS arch resonators. However, it should be noted that the above literature only covers single-mode dynamics. Little attention had been devoted to the dynamics of multi-modes.
Recently, modal interaction has been highly emphasized by researchers. Abdallah et al. [23] experimentally studied an arch resonator in the presence of internal resonance. The results showed that the first and third modes coupled with each other. Hajjaj et al. [24] experimentally investigated an electrothermally tuned and electrostatically actuated arch resonator. The results showed that one-to-one internal resonance between the first and third modes can be activated for some specific parameters. Samanta et al. [25] experimentally detected three different internal resonances in a few-layer MoS 2 resonator. Vyas et al. [26] designed a T-beam nano-resonator using its nonlinear modal interaction due to 1:2 internal resonance. Damian et al. [27] and Antonio et al. [28] demonstrated that for a MEMS arch clamped at both ends the internal resonance can reduce the frequency-amplitude interdependence and stabilized frequency. In addition, analytical research concerning nonlinear modal interaction has also been performed. Younis and Nayfeh [29] explored the possibility of modal interactions between the first and second modes for a MEMS resonator via the method of multiple scales. They found that the first two modes are uncoupled and that energy transfer cannot occur between the involved modes. Pourkiaee et al. [30] examined the nonlinear modal interaction of a doubly clamped piezoelectric MEMS arch. They found that the interaction between the first and third mode can be activated due to the three-to-one internal resonance.
It can be noted that analytical studies on the modal interaction in MEMS arches are scarce. Therefore, this investigation aims to thoroughly study the dynamics of a MEMS arch in the presence of three-to-one internal resonance. To accomplish this, the multiple scale method is used to solve the governing equation of motion. Frequency response curves and force response curves are obtained. In addition, a reduced-order model, using the Galerkin scheme, is derived. The dynamics of the MEMS arch is investigated with respect to the time history, phase-plane portrait, Poincare sections, and Fourier spectrum. Figure 1 shows a MEMS arch resonator with electrostatic actuation which is composed of a deformable top electrode and a stationary bottom electrode. The deformable electrode is actuated by a distributed electrostatic force provided by the stationary electrode. The deformable electrode is modeled as a fixed-fixed arch subjected to a combined AC and DC voltage load. The initial shape of the arch is described by the function w 0 (x) = 0.5b 0 (1 − cos2πx), where b 0 is the initial rise. Recently, modal interaction has been highly emphasized by researchers. Abdallah et al. [23] experimentally studied an arch resonator in the presence of internal resonance. The results showed that the first and third modes coupled with each other. Hajjaj et al. [24] experimentally investigated an electrothermally tuned and electrostatically actuated arch resonator. The results showed that oneto-one internal resonance between the first and third modes can be activated for some specific parameters. Samanta et al. [25] experimentally detected three different internal resonances in a fewlayer MoS2 resonator. Vyas et al. [26] designed a T-beam nano-resonator using its nonlinear modal interaction due to 1:2 internal resonance. Damian et al. [27] and Antonio et al. [28] demonstrated that for a MEMS arch clamped at both ends the internal resonance can reduce the frequency-amplitude interdependence and stabilized frequency. In addition, analytical research concerning nonlinear modal interaction has also been performed. Younis and Nayfeh [29] explored the possibility of modal interactions between the first and second modes for a MEMS resonator via the method of multiple scales. They found that the first two modes are uncoupled and that energy transfer cannot occur between the involved modes. Pourkiaee et al. [30] examined the nonlinear modal interaction of a doubly clamped piezoelectric MEMS arch. They found that the interaction between the first and third mode can be activated due to the three-to-one internal resonance.

Problem Formulation
It can be noted that analytical studies on the modal interaction in MEMS arches are scarce. Therefore, this investigation aims to thoroughly study the dynamics of a MEMS arch in the presence of three-to-one internal resonance. To accomplish this, the multiple scale method is used to solve the governing equation of motion. Frequency response curves and force response curves are obtained. In addition, a reduced-order model, using the Galerkin scheme, is derived. The dynamics of the MEMS arch is investigated with respect to the time history, phase-plane portrait, Poincare sections, and Fourier spectrum. Figure 1 shows a MEMS arch resonator with electrostatic actuation which is composed of a deformable top electrode and a stationary bottom electrode. The deformable electrode is actuated by a distributed electrostatic force provided by the stationary electrode. The deformable electrode is modeled as a fixed-fixed arch subjected to a combined AC and DC voltage load. The initial shape of the arch is described by the function w0(x) = 0.5b0(1 − cos2πx), where b0 is the initial rise. In the present study, the equation of motion was obtained by assuming that: (1) the arch is shallow, i.e., dw0/dx << 1, and, hence, the parallel-plate assumption is valid; (2) the Euler-Bernoulli  In the present study, the equation of motion was obtained by assuming that: (1) the arch is shallow, i.e., dw 0 /dx << 1, and, hence, the parallel-plate assumption is valid; (2) the Euler-Bernoulli beam theory may be used, neglecting the effect of shear and rotary inertia; (3) the simplest viscous damping model may be adopted to model the dissipative mechanisms of the resonator.

Problem Formulation
Therefore, the equation of motion governing the transverse deflection w(x,t) of the arch with length L, width b, and thickness h is expressed as [12,31] and the boundary condition can be expressed as where E is Young's modulus, G is the shear modulus, c is the viscous damping coefficient, ρ is the material density, d is the gap between electrodes, A is the cross-sectional area, I is the moment of inertia, l is the material length scale parameter, P is the axial load due to residual stress, and ε is the dielectric constant of the air. The governing equation can be transformed to the dimensionless forms by introducing the following parameters.
The governing equation and the boundary condition can be cast into the dimensionless form ..
where the overdot indicates the derivative with respect to t, the prime indicates the derivative with respect to x, and

Method of Analysis
In this section we used the method of multiple scales to solve the response of the arch under a combined AC and DC voltage load. To this end, we introduced the time scales T 0 = t, T 1 = εt, and T 2 = ε 2 t, where ε is a small dimensionless bookkeeping parameter. Moreover, in order for nonlinearity to balance the effects of the damping c and the excitation V AC , we rescaled them as ε 2 c and ε 3 V AC . Hence, we sought a solution for the governing Equation (3) in the form  (3) and (4) and equating coefficients of like powers of ε, we obtained where D n ≡ ∂/∂T n . The boundary conditions for all orders are.
One can note that Equation (7) is the static deflection problem of the arch under the DC voltage. As a case study, we considered a MEMS arch resonator made of silicon. Its geometric and material parameters were GPa, and ρ = 2320 kg/m 3 [32]. Here, the Galerkin method was employed to solve Equation (7) [12]. Figure 2b shows variation of the midpoint static deflection of the arch with the DC voltage (V DC ). To validate the results using the Galerkin method, the finite element method (FEM) was also employed to simulate the static deformation of the MEMS arch. It can be noted from Figure 2b that the results via the two methods are in agreement. Here, the results of the finite element method were obtained by using the software COMSOL, as shown in Figure 2a.
where Dn ≡ ∂/∂Tn. The boundary conditions for all orders are.
One can note that Equation (7) is the static deflection problem of the arch under the DC voltage. As a case study, we considered a MEMS arch resonator made of silicon. Its geometric and material parameters were L = 1000 μm, h = 2.4 μm, b = 30 μm, d = 10.1 μm, b0 = 3.6 μm, P = 0 N, E = 167.8 GPa, and ρ = 2320 kg/m 3 [32]. Here, the Galerkin method was employed to solve Equation (7) [12]. Figure  2b shows variation of the midpoint static deflection of the arch with the DC voltage (VDC). To validate the results using the Galerkin method, the finite element method (FEM) was also employed to simulate the static deformation of the MEMS arch. It can be noted from Figure 2b that the results via the two methods are in agreement. Here, the results of the finite element method were obtained by using the software COMSOL, as shown in Figure 2a.
where A m (T 2 ) denotes the complex function to be determined, ω m is nth natural frequency of the system, φ m (x) is the mode function, and cc indicates the complex conjugate of the preceding terms. Much literature has been published on the linear eigenvalue problem of the arch. For more details of the solution methodology, refer to works [12,29]. Here, Figure 3a shows variation of the first three natural frequencies with the DC voltage and Figure 3b shows variations of three times the first natural frequency and the third natural frequency with the DC voltage. In these figures, the results obtained by FEM are given. It is found that the third natural frequency is approximately equal to three times the first natural frequency, being near V DC = 27.95. Hence, there may be a three-to-one internal resonance for a specific range of parameters. The internal resonance is perfectly tuned when V DC ≈ 27.95. In the present investigation, we explored the dynamics of the system in the presence of three-to-one internal resonance between the first and third modes. Here, we assumed that neither of these two modes was involved in an internal resonance with other modes. Due to the presence of damping and internal resonance, only the first and third modes contribute to the long-time dynamic response. As a result, Equation (8) can be expressed as where Am(T2) denotes the complex function to be determined, ωm is nth natural frequency of the system, ϕm(x) is the mode function, and cc indicates the complex conjugate of the preceding terms. Much literature has been published on the linear eigenvalue problem of the arch. For more details of the solution methodology, refer to works [12,29]. Here, Figure 3a shows variation of the first three natural frequencies with the DC voltage and Figure 3b shows variations of three times the first natural frequency and the third natural frequency with the DC voltage. In these figures, the results obtained by FEM are given. It is found that the third natural frequency is approximately equal to three times the first natural frequency, being near VDC = 27.95. Hence, there may be a three-to-one internal resonance for a specific range of parameters. The internal resonance is perfectly tuned when VDC ≈ 27.95. In the present investigation, we explored the dynamics of the system in the presence of three-to-one internal resonance between the first and third modes. Here, we assumed that neither of these two modes was involved in an internal resonance with other modes. Due to the presence of damping and internal resonance, only the first and third modes contribute to the long-time dynamic response. As a result, Equation (8) can be expressed as Next, we offer a detailed investigation on the dynamic of the MEMS arch resonator in the presence of the internal resonance between the first and third modes.

Primary Resonance of the First Mode
To describe the nearness of the ω3 to 3ω1 and the excitation frequency Ω to ω1, we let where σ1 and σ2 are detuning parameters. Substituting Equation (13)   Next, we offer a detailed investigation on the dynamic of the MEMS arch resonator in the presence of the internal resonance between the first and third modes.

Primary Resonance of the First Mode
To describe the nearness of the ω 3 to 3ω 1 and the excitation frequency Ω to ω 1 , we let where σ 1 and σ 2 are detuning parameters. Substituting Equation (13) into Equation (9) yields where A particular solution of Equation (9) can be expressed as The ψ ij are solutions of the following boundary-value problem: This has the boundary conditions ψ ij = 0 and ψ ij = 0 at x = 0, 1 where the operator M is defined as the following equation: Substituting Equations (13) and (16) into Equation (10) and considering Equation (14), we obtained where the overdot indicates the differential with respect to T 2 , NST denotes terms that do not produce secular effects, and the functions χ's and F(x) are defined in Appendix A.
Due to the fact that the homogeneous part of Equation (23) has a nontrivial solution, Equation (23) has a nontrivial solution only if the solvability conditions are satisfied. Because this problem is self-adjoint, it can be demanded that the right hand of Equation (23) be orthogonal to φ m (x)exp(−iω m T 0 ) and φ n (x)exp(−iω n T 0 ). Performing these manipulations, we were able to obtain the modulation equations where f 1 , Λ i , and S ij (i,j = 1,3) are given in Appendix. f 1 is the projection of the excitation force onto the first mode. S ij are nonlinear coefficients due to electric and geometric sources and the Λ i are the nonlinear interaction coefficients between the first and third modes. If the nonlinear interaction coefficients Λ i are identically zero, the involving modes are orthogonal in the nonlinear sense and there is no potential for energy transfer between them. Next, we introduced the Cartesian transformation into Equation (24), separated the obtained results into those real and imaginary, and yielded where s 1 = σ 2 and s 2 = 3σ 2 − σ 1 . The equilibrium solution of the system was obtained by setting p 1 = p 3 = q 1 = q 3 = 0. Then, the resulting algebraic equation was numerically solved and the eigenvalue of the corresponding Jacobian matrix found; these were used for examining the stability of the equilibrium solution.
where a 1 and a 3 are the amplitudes of the first and second modes, respectively.

Equilibrium Solutions and Their Stability
In the present section, numerical examples are shown to reveal characteristics of response. We consider again the MEMS arch resonator of Section 3. Here, we chose a DC voltage of 23.80, where a three-to-one internal resonance occurs, as shown in Section 3. The AC voltage was chosen to be 0.60. The corresponding dimensionless values were α 1 = 106.3 and α 2 = 0.023. The damping coefficient c was equal to 0.01. The first and third dimensionless natural frequencies were ω 1 = 42.9127 and ω 3 = 127.9320. These parameters were selected as fixed parameters, if no others were assigned. Figure 4 shows the frequency response curves for the first and third modes as functions of the detuning parameter σ 2 near the primary resonance of the first mode. In these figures, the solid black lines denote the stable solution, the dashed blue lines denote the unstable solution, and the thick red solid lines denote unstable foci.
It follows from Figure 4 that the curves bent to the left, indicating softening spring type behavior. As σ 2 increased from a small value of σ 2 , where two stable equilibrium solutions coexisted, the amplitude for smaller equilibrium solutions grew all the way until a saddle node bifurcation occurred at SN 1 . Increasing σ 2 beyond SN 1 , the response jumped to the other equilibrium solution. For the larger equilibrium solutions, the amplitude of the first mode decreased all the way, whereas the amplitude of the second mode decreased first and then grew until a saddle node bifurcation occurred at SN 2 . This phenomenon indicates that the energy was transferred from the first mode to the second mode. Moreover, as σ 2 increased from a small value, the response lost stability via a Hopf bifurcation at H 1 and then regained stability via a reverse Hopf bifurcation at H 2 . It is worth mentioning that the response of the system between the two Hopf bifurcations (H 1 and H 2 ) may present rich nonlinear dynamic behaviors which can be examined by a numerical method. To reveal the characteristics of the response in this region, the Galerkin truncation scheme was employed, which the dynamic solutions presented in Section 5.
In the present section, numerical examples are shown to reveal characteristics of response. We consider again the MEMS arch resonator of Section 3. Here, we chose a DC voltage of 23.80, where a three-to-one internal resonance occurs, as shown in Section 3. The AC voltage was chosen to be 0.60. The corresponding dimensionless values were α1 = 106.3 and α2 = 0.023. The damping coefficient c was equal to 0.01. The first and third dimensionless natural frequencies were ω1 = 42.9127 and ω3 = 127.9320. These parameters were selected as fixed parameters, if no others were assigned. Figure 4 shows the frequency response curves for the first and third modes as functions of the detuning parameter σ2 near the primary resonance of the first mode. In these figures, the solid black lines denote the stable solution, the dashed blue lines denote the unstable solution, and the thick red solid lines denote unstable foci.   It follows from Figure 4 that the curves bent to the left, indicating softening spring type behavior. As σ2 increased from a small value of σ2, where two stable equilibrium solutions coexisted, the amplitude for smaller equilibrium solutions grew all the way until a saddle node bifurcation occurred at SN1. Increasing σ2 beyond SN1, the response jumped to the other equilibrium solution. For the larger equilibrium solutions, the amplitude of the first mode decreased all the way, whereas the amplitude of the second mode decreased first and then grew until a saddle node bifurcation occurred at SN2. This phenomenon indicates that the energy was transferred from the first mode to the second mode. Moreover, as σ2 increased from a small value, the response lost stability via a Hopf bifurcation at H1 and then regained stability via a reverse Hopf bifurcation at H2. It is worth mentioning that the response of the system between the two Hopf bifurcations (H1 and H2) may present rich nonlinear dynamic behaviors which can be examined by a numerical method. To reveal the characteristics of the response in this region, the Galerkin truncation scheme was employed, which the dynamic solutions presented in Section 5. Figure 5 shows the influence of the internal detuning parameter σ1 (i.e., the DC voltage) on the frequency response curves. It follows from Figure 5 that with the increase of σ1, the amplitude of the indirectly third mode decreases. The interval between the two Hopf bifurcation decreases until it vanishes. The variation of the response amplitude with excitation VAC was examined when σ1 = −0.8059 and σ2 = −0.3. The force-response curves are shown in Figure 6. As the VAC increases from 0, the amplitude of the response grows until a saddle node bifurcation occurs at SN1, resulting in a jump of the response to the other stable branch. As the VAC increases further, the response loses stability via another saddle-node bifurcation at SN3. Increasing the f beyond SN3, the amplitude of the response for the first mode jumps up to another stable branch. To the left of this saddle node point, a second saddle node SN2 occurs. As VAC decrease beyond SN2, the response loses stability via a Hopf bifurcation at H1 and regained stability via a Hopf bifurcation at H2. Then, this equilibrium solution encounters a Hopf bifurcation again at H3, resulting in the response losing stability and then jumping to the other larger stability branch with increasing VAC. The variation of the response amplitude with excitation V AC was examined when σ 1 = −0.8059 and σ 2 = −0.3. The force-response curves are shown in Figure 6. As the V AC increases from 0, the amplitude of the response grows until a saddle node bifurcation occurs at SN 1 , resulting in a jump of the response to the other stable branch. As the V AC increases further, the response loses stability via another saddle-node bifurcation at SN 3 . Increasing the f beyond SN 3 , the amplitude of the response for the first mode jumps up to another stable branch. To the left of this saddle node point, a second saddle node SN 2 occurs. As V AC decrease beyond SN 2 , the response loses stability via a Hopf bifurcation at H 1 and regained stability via a Hopf bifurcation at H 2 . Then, this equilibrium solution encounters a Hopf bifurcation again at H 3 , resulting in the response losing stability and then jumping to the other larger stability branch with increasing V AC . It may be noted that the response of the system has many Hopf bifurcations and saddle node bifurcations due to internal resonance. To investigate the effect of the internal resonance on the response, Figure 7 shows force response curves for σ1 = 0.005 and σ1 = 0.6307. As can be seen from Figure 7, as σ1 increases, Hopf bifurcation diminishes. When σ1 = 0.6307, i.e., VDC = 30.80, the force response curves show behavior similar to the single-modal dynamics of the system.

Dynamic Solutions
When inspecting frequency responses and amplitude response curves, it can be found that many bifurcations may occur with variation of parameters. This calls for analysis of the dynamic response It may be noted that the response of the system has many Hopf bifurcations and saddle node bifurcations due to internal resonance. To investigate the effect of the internal resonance on the response, Figure 7 shows force response curves for σ 1 = 0.005 and σ 1 = 0.6307. As can be seen from Figure 7, as σ 1 increases, Hopf bifurcation diminishes. When σ 1 = 0.6307, i.e., V DC = 30.80, the force response curves show behavior similar to the single-modal dynamics of the system. It may be noted that the response of the system has many Hopf bifurcations and saddle node bifurcations due to internal resonance. To investigate the effect of the internal resonance on the response, Figure 7 shows force response curves for σ1 = 0.005 and σ1 = 0.6307. As can be seen from Figure 7, as σ1 increases, Hopf bifurcation diminishes. When σ1 = 0.6307, i.e., VDC = 30.80, the force response curves show behavior similar to the single-modal dynamics of the system.

Dynamic Solutions
When inspecting frequency responses and amplitude response curves, it can be found that many bifurcations may occur with variation of parameters. This calls for analysis of the dynamic response

Dynamic Solutions
When inspecting frequency responses and amplitude response curves, it can be found that many bifurcations may occur with variation of parameters. This calls for analysis of the dynamic response of the system. Hence, in the present section we explore the dynamic behavior caused by these bifurcations (especially the Hopf bifurcation points) using the numerical method. To this end, we employed the Galerkin scheme to discretize the governing equation into a set of nonlinear ordinary differential equations. We assumed the response of the MEMS arch in the form (31) where N is the number of the modes retained in the discrete scheme, φ i (x) is the eigen function of the straight clamped-clamped beam, and q i (t) is the generalized coordinate of the transverse displacement. Substituting Equation (31) into Equation (3), using the Taylor expansion around the w s (x), multiplying the resultant equation by φ j (x), and then integrating the result over the domain yielded the second-order nonlinear ordinary differential equation ijml q j q m q l = F i cos(Ωt) (32) where i = 1, 2, ..., N, and the overdot denotes the derivative with respective to time.
ijm , and K (1) ijml are the linear, quadratic, and cubic stiffness coefficients, respectively. F i and Ω are the amplitude and frequency of the excitation, respectively. These coefficients have been included in Appendix B. Equation (32) can be numerically integrated using the Runge-Kutta technique to obtain the dynamic solution.
For the Galerkin scheme, one key problem is its convergence. Figure 8a,b show the variation of the first two frequencies with V DC respectively, using three different values of N, namely, 3, 4, and 5. It should be noted that the first two frequencies are accurate enough when a 4-mode Galerkin truncation is employed. With regards to the parameters c = 0.01, V DC = 23.80, and V AC = 0.6, the response of the first two modes for the three different values of N are shown in Figure 8c,d, respectively. Similarly to the conclusion of natural frequency, the 4-mode Galerkin truncation is convergent. Hence, in the present paper, the dynamic solution of the system was solved by the 4-mode Galerkin scheme.
To simulate dynamic behavior, we chose the same system parameters with Figure 4 and set the excitation frequency to 42.3627, i.e., σ 2 = 0.55, in Figure 4. Figure 9a-d show the phase-plane portraits, Fourier spectrum, and projection of the response onto the phase plane of the first and third mode, respectively. It should be noted that the response contains components of the first and third modes, which indicates that the energy is transferred from the first mode to the third mode.
Decreasing the excitation frequency to 42.1127, which corresponds to a stable point at σ 2 = −0.8 in Figure 4, Figure 10a-d show the phase-plane portraits, Fourier spectrum, and projection of the response onto the phase plane of the first and third mode, respectively. Compared with Figure 9, the response amplitude of the first mode in Figure 10 is larger. Also, the contribution of the third mode to the response is smaller.
Slightly decreasing the excitation frequency to 42.1027, which corresponds the Hopf point H 2 at σ 2 = −0.81, Figure 11a-d show the response of the system in terms of time histories, phase portraits, Fourier spectrum, and Poincare sections. It should be noted that the Fourier spectrum contains some other harmonic components and the Poincare sections present a closed loop. These phenomena indicate that the response is quasiperiodic. In addition, for the other values of the detuning parameter σ 2 between H 1 and H 2 in Figure 4a,b, the response remains quasiperiodic and a similar behavior is obtained. Further, decreasing the excitation frequency to less than 40.6812, namely, σ 2 = −2.2, the To simulate dynamic behavior, we chose the same system parameters with Figure 4 and set the excitation frequency to 42.3627, i.e., σ2 = 0.55, in Figure 4. Figure 9a-d show the phase-plane portraits, Fourier spectrum, and projection of the response onto the phase plane of the first and third mode, respectively. It should be noted that the response contains components of the first and third modes, which indicates that the energy is transferred from the first mode to the third mode. Decreasing the excitation frequency to 42.1127, which corresponds to a stable point at σ2 = −0.8 in Figure 4, Figure 10a-d show the phase-plane portraits, Fourier spectrum, and projection of the response onto the phase plane of the first and third mode, respectively. Compared with Figure 9, the response amplitude of the first mode in Figure 10 is larger. Also, the contribution of the third mode to the response is smaller.
Slightly decreasing the excitation frequency to 42.1027, which corresponds the Hopf point H2 at

Conclusions
In this paper, the nonlinear response of a MEMS arch resonator with three-to-one internal resonance was investigated. Size effect and axial load were also considered. The multiple scales method was employed to solve the governing equation. The frequency and force response curves were obtained in the presence of three-to-one internal resonance. Results showed that even though the ratio between the first and second natural frequency of the MEMS arch can be three-to-one, the

Conclusions
In this paper, the nonlinear response of a MEMS arch resonator with three-to-one internal resonance was investigated. Size effect and axial load were also considered. The multiple scales method was employed to solve the governing equation. The frequency and force response curves were obtained in the presence of three-to-one internal resonance. Results showed that even though the ratio between the first and second natural frequency of the MEMS arch can be three-to-one, the two modes cannot be coupled in the nonlinear sense and hence cannot exchange energy with each other. By contrast, the first and third mode can be coupled due to three-to-one internal resonance. Moreover, many nonlinear phenomena, such as softening spring behaviors, jumping, and Hopf bifurcation, can be detected due to the existence of internal resonance. Furthermore, the effect of different parameters on internal resonance was also studied. The results show that increasing the DC voltage weakens modal interaction due to internal resonance while decreasing the DC voltage level enhances modal interactions.
In addition, to simulate the dynamic solution of the system, a reduced-order model was derived by employing the Galerkin scheme. The dynamics of the MEMS arch were investigated reharding time history, phase-plane portrait, Poincare sections, and the Fourier spectrum. Periodic, quasiperiodic, and mixed responses were observed in the motion, which demonstrates the role of modal interactions due to internal resonance. The presented results can be used in the design and optimization of novel MEMS resonators and give insight into how modal interactions can affect the stability and the resonant responses of MEMS arch resonators.