Multistability of the Vibrating System of a Micro Resonator

Abstract

Multiple attractors and their fractal basins of attraction can lead to the loss of global stability and integrity of Micro Electro Mechanical Systems (MEMS). In this paper, multistability of a class of electrostatic bilateral capacitive micro-resonator is researched in detail. First, the dynamical model is established and made dimensionless. Second, via the perturbating method and the numerical description of basins of attraction, the multiple periodic motions under primary resonance are discussed. It is found that the variation of AC voltage can induce safe jump of the micro resonator. In addition, with the increase of the amplitude of AC voltage, hidden attractors and chaos appear. The results may have some potential value in the design of MEMS devices.

1. Introduction

Multistability, i.e., the coexistence of multiple attractors, is a common dynamical phenomenon in MEMS/NEMS [1,2]. Based on it, there are many applications such as MEMS-based memory [3] and switches [4]. In addition, considering the loss of global stability that multistability may trigger, there are some devices that should avoid the appearance of multiple attractors in their vibrating systems, such as filters [5], microvalves [6], and micro-relays [7]. As one of the fastest developing MEMS products [8], electrostatic micro-resonators should assure that the resonators undergo periodic vibration whose amplitudes vary continuously with the driven voltages. However, in practical applications of electrostatic micro-resonators [9], there are many complex dynamic behaviors such as multistability [10,11], quasi-periodic motion [12] and periodic-n motion [13], chaos [14,15], and pull-in instability.

It is of great significance to study the multistability and necessary conditions for inducing it either for avoiding this phenomenon or making use of it. Thus, multistability of vibrating systems of micro resonators has been studied experimentally and numerically during these decades [16]. Via experiments, Mohammadreza investigated the dynamic response of an electrostatic micro-actuator in the vicinity of the primary resonance and the parametric one [17]. Siewe et al. [18] studied the vibration of a double-side MEMS resonator numerically and found the variation of the driven voltage could induce the coexistence of chaos and quasi-periodic motions. Shang et al. [19] found the coexisting chaos and dynamical pull-in in the vibrating system of a single-side electrostatic micro sensor. Haghighi et al. [20] found the coexisting periodic-n motion and the chaotic motion of micromechanical resonators with electrostatic forces on both sides, and then discussed the global bifurcation of its vibrating system by approximately expressing its homoclinic orbits as the ones of a typical duffing equation. When amplifying signals of a nanomechanical duffing resonator, Almog et al. [21] found that multistability was an interesting dynamical phenomenon of nonlinear systems and could be explored for many applications. Gusso et al. [22] studied chaos of a typical micro/nanoelectromechanical beam resonator with two-sided electrodes experimentally and observed multiple attractors in a significant region of the relevant parameter space, involving periodic and chaotic attractors. By applying cell-mapping method to depict the basins of attraction for all the attractors, they also found that the basin boundaries were fractal under certain conditions of the excitations, indicating that the attractors are strongly intermingled. Liu et al. [23] applied the method of multiple scales (MMS) to analyze the multiple periodic motions induced by the local bifurcation, and used the Melnikov method to predict necessary conditions for chaos and its control. The corresponding numerical results were also presented by the basins of attraction and spectrum diagrams. Angelo et al. [24] investigated the effect of the linear and nonlinear stiffness terms and damping coefficients on dynamical behaviors of a microelectromechanical resonator and controlled the chaotic motion by forcing it into an orbit obtained analytically via the harmonic balance method. However, most study concentrated on describing or observing the phenomenon itself rather than studying its mechanism, which is still not that clear yet.

To this end, we consider a typical electrostatic driven bilateral capacitive micro-resonator and study the possible multistability and its mechanism in its vibrating system. The paper is organized as follows. In Section 2, the dynamical model is constructed and made dimensionless. In Section 3 and Section 4, two different cases for coexisting multiple periodic attractors, fractal basins of attraction, and other complex attractors of the systems are discussed both theoretically and numerically. In Section 5, the conclusions are presented.

2. Dynamical Model

We choose to study a class of bilateral micro resonator whose simplified diagram is shown in Figure 1. The driven forces on the resonator are electrostatic ones between the moving electrode and the fixed electrode [25]. The driven voltage in Figure 1 is the combination of alternate current (AC) and direct current (DC) actuation. In the figure, x is the vertical displacement of the moving electrode at moment t, d the initial gap width between the moving electrode and each fixed one, V_b the DC bias voltage, V_ACsinΩt the AC voltage where V_AC is the amplitude and Ω the frequency. Suppose that the amplitude of the AC voltage V_AC is much lower than the bias DC voltage V_b, i.e., $V_{AC}\ll V_b$. According to the Second Law of Newton, the vibrating system of the moving electrode can be expressed as a nonlinear system as follows:

$$m\frac{d^{2}x}{d{t}^2}+c\frac{dx}{dt}+k_{1}x+k_2{x}^3=\frac{C_0}{2{\left(d-x\right)}^2}{\left(V_b+V_{AC}sin\Omega t\right)}^2-\frac{C_0{V}_b^2}{2{\left(d+x\right)}^2}$$

(1)

where m represents the effective lumped mass of the moving electrode, k₁ its linear mechanical stiffness, k₂ its cubic nonlinear stiffness, c the damping coefficient, C₀ the initial capacitance of the parallel-plate structure.

Introducing the following dimensionless variables

$${\omega }_0=\sqrt{\frac{k_1}{m}},\omega =\frac{\Omega }{{\omega }_0},\mu =\frac{c}{m{\omega }_0},\alpha =\frac{k_2{d}^2}{m{\omega }_0^2},\beta =\frac{C_0{V}_b^2}{2k_1{d}^3},\gamma =\frac{V_{AC}}{V_b},T={\omega }_{0}t,u=\frac{x}{d},\dot{u}=\frac{du}{dT}$$

(2)

and substituting Equation (2) into Equation (1), one can obtain that

$$\ddot{u}+\mu \dot{u}+u+\alpha u^3=\frac{\beta }{{\left(1-u\right)}^2}{\left(1+\gamma sin\omega T\right)}^2-\frac{\beta }{{\left(1+u\right)}^2}$$

(3)

which is a dimensionless system. Since in the original system (1), the viscous damping coefficient of air $c$ is very tiny, and $V_{AC}\ll V_b$, the parameters $\mu$ and $\gamma$ in (3) will be both small and can be considered as perturbed parameters. Thus, considering $\mu =0$ and $\gamma =0$ in Equation (3), one has the unperturbed system that can be expressed as below:

$$\dot{u}=v,\hspace{1em}\dot{v}=-u-\alpha u^3+\frac{\beta }{{\left(1-u\right)}^2}-\frac{\beta }{{\left(1+u\right)}^2}.$$

(4)

Letting the right side of Equation (4) be zero, one can determine equilibria of the dimensionless system (3). Equation (5) is a Hamilton system with the Hamiltonian

$$H\left(u,v\right)=\frac{1}{2}{v}^2+\frac{1}{2}{u}^2+\frac{\alpha }{4}{u}^4-\frac{\beta }{1-u}-\frac{\beta }{1+u}+2\beta$$

(5)

and the function of potential energy (P.E.)

$$V\left(u\right)=\frac{1}{2}{u}^2+\frac{\alpha }{4}{u}^4-\frac{\beta }{1-u}-\frac{\beta }{1+u}+2\beta .$$

(6)

Concerning Equation (4), the number of the equilibria, and the shapes and positions of the possible potential wells of the unperturbed system (4) depend on the parameters $\alpha$ and $\beta$. The same as in [20], the values of the parameters in the system (1) are given by:

$$m=5\times {12}^{-12} \mathrm{k}\mathrm{g},c=5\times {12}^{-8} \mathrm{k}\mathrm{g}/s, k_1=5\hspace{0.17em} \mathsf{\mu }\mathrm{N}/\mathsf{\mu }\mathrm{m}, k_2=15\hspace{0.17em} \mathsf{\mu }\mathrm{N}/{\mathsf{\mu }\mathrm{m}}^3, d=2\mathsf{\mu }\mathrm{m}, C_0=1.875\times {10}^{-18} \mathrm{m}\mathrm{F}.$$

(7)

Accordingly, in system (4), $\alpha =12$.

Different equilibria and potential energy diagrams of the unperturbed system under different values of the parameter $\beta$ can be seen in Figure 2. It shows that there are three P.E. poles when β = 0.211, five P.E. poles when $\beta$ increases to 0.338, and only one P.E. pole when β increases to 0.6. Under different values of $\beta$, the potential wells and unperturbed orbits are shown in Figure 3. When β = 0.211, there are three equilibria (two non-trivial equilibria are saddles and the origin is a center) as well as one well surrounded by heteroclinic orbits (see Figure 3a). As β increases to 0.338, there will be five equilibria among which two non-trivial equilibria S₁ (−0.196339,0) and S₂ (0.196339,0) are centers of the two wells surrounded by homoclinic orbits; the other three equilibria are unstable. When β = 0.6, no wells or non-trivial equilibria of the unperturbed system (4) exist. The P.E. poles in Figure 2 correspond to the fixed points shown in Figure 3. Therefore, according to Equations (2) and (7), when the structural parameters are fixed, the number of centers will depend on the value of DC bias voltage $V_b$: when the DC bias voltage is very low, there will be a center of the system (4) as well as a stable point attractor of the system (3) without AC voltage. Under a higher DC bias voltage, there may be two centers of the system (4). As is well known, periodic vibration can often be attributed to the perturbation of the centers. Since the number and the location of the centers in Figure 3a,b are totally different, the mechanism for the possible multiple periodic attractors of the vibrating system of the micro-resonator can be different as well. Therefore, in Section 3 and Section 4, we discuss the different mechanism of multi-stability for these two different cases, i.e., the only center (the origin) and the two non-trivial centers, respectively.

3. Multiple Periodic Attractors in the Neighborhood of the Origin

Considering the case where the DC bias voltage is low, and the periodic vibration of the microstructure is induced by the perturbation of the only center (see Figure 3a, where V_b = 3 V), one may use the Method of Multiple Scales (MMS) to analyze the periodic solutions in the neighborhood of the origin. Expanding the fractional terms of the dimensionless system (3) as Taylor series in the neighborhood of u = 0, and neglecting the higher-order-than-three terms of u, one has:

$$\ddot{u}+\mu \dot{u}+u+\alpha u^3=2\beta \gamma \sin \omega T+4\beta u+4u\beta \gamma \sin \omega T+6u^2\beta \gamma \sin \omega T+8\beta u^3+8u^3\beta \gamma \sin \omega T.$$

(8)

As mentioned in Section 2, the values of the parameters $\mu$ and $\gamma$ in the above system are small; one can introduce a small parameter $\epsilon$ satisfying $0<\epsilon \ll 1$, and can re-scale the two parameters in the system (8) as:

$$\mu ={\epsilon }^2\tilde{\mu }, \gamma ={\epsilon }^2\tilde{\gamma }.$$

(9)

Then Equation (8) becomes

$$\ddot{u}+{\tilde{\omega }}^{2}u=-{\epsilon }^2\tilde{\mu }\dot{u}+2{\epsilon }^2\beta \tilde{\gamma }\sin \omega T+4u{\epsilon }^2\beta \tilde{\gamma }\sin \omega T+6u^2{\epsilon }^2\beta \tilde{\gamma }\sin \omega T-P_1{u}^3+8u^3{\epsilon }^2\beta \tilde{\gamma }\sin \omega T.$$

(10)

where

$${\tilde{\omega }}^2=1-4\beta , P_1=\alpha -8\beta .$$

(11)

To apply MMS, one may rescale some terms in the system (10) that

$$\omega =\tilde{\omega }+\epsilon \sigma ,\hspace{0.17em}u=\epsilon u_1+{\epsilon }^2{u}_2+\cdots ,\hspace{0.17em}\sigma =O(1).$$

(12)

and

$$T_i={\epsilon }^{i}T,D_i=\frac{\partial }{\partial T_i},\frac{d}{dT}={\displaystyle \sum _{i=0}^n{\epsilon }^i{D}_i}\hspace{1em}(i=0,1,2,\cdots )$$

(13)

Comparing the coefficients of ${\epsilon }^1,{\epsilon }^2$, and ${\epsilon }^3$ in the system (10), respectively, one obtains that

$${\epsilon }^1:\hspace{1em}\hspace{0.17em}{D}_0{}^2{u}_1+{\omega }^2{u}_1=0,$$

(14)

$${\epsilon }^2:\hspace{1em}{D}_0{}^2{u}_2+{\omega }^2{u}_2=-2D_1{D}_0{u}_1+2\omega \sigma u_1+2\beta \tilde{\gamma }\sin \omega T,$$

(15)

and

$${\epsilon }^3:\hspace{1em}{D}_0{}^2{u}_3+{\omega }^2{u}_3=-2D_1{D}_0{u}_2-\tilde{\mu }{D}_0{u}_1-D_1{}^2{u}_1+2u_2\omega \sigma -{\sigma }^2{u}_1-2D_2{D}_0{u}_1-P_1{u}_1{}^3+4\beta \tilde{\gamma }{u}_1\sin \omega T.$$

(16)

To solve Equation (14), one can assume that

$$u_1=A_1(T_1,T_2)e^{i\omega T_0}+{\overline{A}}_1(T_1,T_2)e^{-i\omega T_0},$$

(17)

where

$$A_1=\frac{a(T_1,T_2)}{2}{e}^{i\theta (T_1,T_2)}.$$

(18)

Substituting Equations (17) and (18) into Equation (15), and eliminating the secular terms of Equation (15), one will have:

$$D_1{A}_1=-\frac{\beta \tilde{\gamma }}{2\omega }-i\sigma A_1.$$

(19)

Solving Equation (15), one may assume:

$$u_2=A_2(T_2)e^{i\omega T_0}+{\overline{A}}_2(T_2)e^{-i\omega T_0}.$$

(20)

Substituting Equation (20) into Equation (16), and eliminating secular terms of Equation (16), one will obtain:

$$D_2{A}_1=-\frac{\tilde{\mu }}{2}{A}_1+\frac{\beta \tilde{\gamma }}{2\omega }-\frac{\sigma \beta \tilde{\gamma }}{4{\omega }^2}+\frac{3i{P}_1{A}_1^2{\overline{A}}_1}{2\omega }.$$

(21)

Since

$${\dot{A}}_1\approx D_0{A}_1+\epsilon D_1{A}_1+{\epsilon }^2{D}_2{A}_1,$$

(22)

Substituting Equations (19) and (21) into Equation (22), and expressing it by the original dimensionless parameters of Equation (3), one has:

$$\begin{array}{l}\epsilon \dot{a}=-\frac{\mu }{2}(\epsilon a)-P_2\cos \theta ,\\ (\epsilon a)\dot{\theta }=-(\omega -\tilde{\omega })(\epsilon a)+\frac{3P_1{(\epsilon a)}^3}{8\omega }+P_2\sin \theta .\end{array}$$

(23)

where

$$P_2=\frac{(3\omega -\tilde{\omega })\beta \gamma }{2{\omega }^2}.$$

(24)

According to Equation (18), it is obvious that the amplitude of the periodic solution $a$ is the function of the time scale $T_1$ where $T_1$ is a one-order term of $\epsilon$; thus, one can assume the amplitude of the solution u of Equation (10) $\tilde{a}$ as:

$$\epsilon a=\tilde{a}.$$

(25)

Letting $\dot{a}=0,$ and $\dot{\theta }=0,$ one can obtain:

$$-\frac{\mu }{2}\tilde{a}=P_2\cos \theta ,\hspace{1em}(\omega -\tilde{\omega })\tilde{a}-\frac{3P_1{\tilde{a}}^3}{8\omega }=P_2\sin \theta .$$

(26)

Eliminating the triangulation function of Equation (26), one can get:

$$\frac{{\mu }^2}{4}{a}^2+{(\omega -\tilde{\omega }-\frac{3P_1{a}^2}{8\omega })}^2{a}^2=\frac{{(3\omega -\tilde{\omega })}^2}{4{\omega }^4}{\beta }^2{\gamma }^2.$$

(27)

According to Equations (17), (18) and (25), the periodic solution can be expressed as:

$$u\approx \tilde{a}\cos (\omega T+\theta ).$$

(28)

To determine the stability of the periodic solutions, one can get the corresponding characteristic equation of the periodic solution based on Equation (23). It shows that the periodic solution will lose its stability when its amplitude $\tilde{a}$ satisfies:

$$(8(\omega -\tilde{\omega })-\frac{9{\tilde{a}}^2{P}_1}{\omega })(8(\omega -\tilde{\omega })-\frac{3{\tilde{a}}^2{P}_1}{\omega })\ge 16{\mu }^2.$$

(29)

Based on Equations (26)–(29), the variation of the amplitude of the periodic solutions of the system (3) and their stability with AC voltage is shown in Figure 4 where the frequency and the amplitude of AC voltage are considered as the control parameters in Figure 4a and Figure 4b, respectively. In Figure 4a, where V_AC = 0.01 V, when ω is lower than 0.45, there is only one periodic attractor in the system (3) whose amplitude changes continuously with the increase in ω. Comparatively, when ω ranges from 0.46 to 0.69, the global dynamical behaviors of the system (3) will change to bistable periodic attractors, which can be attributed to Hopf bifurcation. As ω continues to increase from 0.7, the periodic attractor with the higher amplitude will disappear, only the periodic attractor with the lower amplitude will exist, and its amplitude will decrease continuously with the increase of ω. Similarly, the change in global dynamical behaviors in Figure 4b also shows that in a certain range of V_AC, there will be two periodic attractors coexisting, which can be due to Hopf bifurcation of the system (3). The accuracy of the theoretical prediction in Figure 4 are verified by the numerical results.

Figure 4 demonstrates that the parameters ω and VAC can induce the coexistence of bistability, meaning that under fixed values of parameters of system (3), different initial conditions may lead to different periodic attractors. Accordingly, it is necessary to classify the basins of attraction for the two different periodic attractors. Here, the 4th order Runge-Kutta approach and the cell-mapping method are applied to depict the basins of attraction of the system (3). The time step is taken as 1/10² of the period of excitation. To investigate the long-term dynamical behaviors, it is supposed that an initial condition will be safe if the vibration in this initial condition keeps satisfying |u(T)| < 1 within 10⁵ excited circles; otherwise, the micro resonator will undergo pull in [20]. The union of all initial conditions leading to the same periodic motion will be the basin of attraction for that attractor which will surely be marked in the same color in the initial plane. The basins of attraction of system (3) are drawn in sufficiently large ranges for the initial position and velocity of the proof mass defined as |u(0)| < 1 and |$\dot{\mathrm{u}}\left(0\right)$| < 1.5 by generating an 200 × 100 array of initial points. The change of attractors and the area and nature of their basins of attraction with frequency ω is shown in Figure 5 where the amplitude of AC voltage V_AC is fixed as 0.01 V.

According to Figure 5, with the increase in parameter ω, the number of attractors and the boundary of basins of attraction will both change. When ω = 0.45 in system (3), there will only be one periodic attractor whose basin of attraction is comparatively bigger with a smooth boundary (see Figure 5a1,a2). However, with a small increase of ω, i.e., ω = 0.45, the global dynamics are totally different, as shown in Figure 5b1,b2 where two periodic attractors coexist, whose basins of attraction mix each other and are both fractal. It means that the dynamical behavior of system (3) is highly sensitive to initial conditions. In other words, system (3) may undergo a safe jump. A similar phenomenon can be seen in Figure 5c1,c2 under a higher ω. As ω increases, the basin of attraction of the periodic attractor with the higher amplitude becomes small (see the red regions in Figure 5b2,c2,d2). Specifically, in Figure 5d2, the regions of attraction are almost blue, and there is very little area of basin of attraction for that periodic attractor. As ω = 0.70 (see Figure 5e1,e2), the periodic attractor with the higher amplitude disappears, and there is only the other attractor whose basin of attraction is almost the same as that in Figure 4b, showing that when the frequency ω increases enough, the periodic attractor with the lower amplitude replaces the initial one.

4. Multistability in the Neighborhood of Non-Trivial Equilibria

In this section, the case that the DC voltage is higher, and the periodic vibration of the microstructure is induced by the perturbation of the two non-trivial centers (see Figure 3b) is considered; thus, we set β = 0.338, i.e., $V_b=3.8\mathrm{V}$. In addition, we consider the effect of AC voltage on the global dynamics of the system (3). To begin with, setting

$$\epsilon \widehat{u}=u\mp u_c$$

(30)

where $u_c$ is the abscissa of the right center (see S₁ in Figure 3), rescaling he two parameters $\mu$ and $\gamma$ in the system (3) by Equation (9), expanding the fractional terms of the dimensionless system (3) as a Taylor series in the neighborhood of the non-trivial equilibria and ignoring the higher-order-than-cubic terms of $\widehat{u}$, the system (3) becomes

$$\ddot{\widehat{u}}=-{\epsilon }^2\tilde{\mu }\dot{\widehat{u}}-{\widehat{\omega }}^2\widehat{u}+\epsilon Q_1{\widehat{u}}^2+{\epsilon }^2{Q}_2{\widehat{u}}^3+\frac{2\epsilon \beta \tilde{\gamma }\sin \omega T}{{(1\mp u_c)}^2}+\frac{4{\epsilon }^2\beta \tilde{\gamma }\sin \omega T}{{(1\mp u_c)}^3}\widehat{u}+\frac{6{\epsilon }^3\beta \tilde{\gamma }\sin \omega T}{{(1\mp u_c)}^4}{\widehat{u}}^2,$$

(31)

where

$${\widehat{\omega }}^2=1+3\alpha u_c{}^2-\frac{4\beta (1+3u_c{}^2)}{{(1-u_c{}^2)}^3},\hspace{0.17em}{Q}_1=\pm 3u_c(-\alpha +\frac{8\beta (1+u_c{}^2)}{{(1-u_c{}^2)}^4}),\hspace{0.17em}{Q}_2=-\alpha +\frac{8\beta (1+10u_c{}^2+5u_c{}^4)}{{(1-u_c{}^2)}^5}.$$

(32)

To apply the Method of Multiple Scale in Equation (32), one can assume in this equation that:

$$\widehat{\omega }=\omega +\epsilon \widehat{\sigma },\hspace{0.17em}\widehat{u}={\widehat{u}}_0+\epsilon {\widehat{u}}_1+{\epsilon }^2{\widehat{u}}_2+\cdots ,\hspace{0.17em}\widehat{\sigma }=O(1).$$

(33)

Comparing the coefficients of ${\epsilon }^1,{\epsilon }^2$ and ${\epsilon }^3$, one has:

$${\epsilon }^0:\hspace{1em}{D}_0{}^2{\widehat{u}}_0+{\omega }^2{\widehat{u}}_0=0,$$

(34)

$${\epsilon }^1:\hspace{1em}{D}_0{}^2{\widehat{u}}_1+{\omega }^2{\widehat{u}}_1=-2D_1{D}_0{\widehat{u}}_0+2\omega \widehat{\sigma }{\widehat{u}}_0+\frac{2\beta \tilde{\gamma }\sin \omega T}{{(1\mp u_c)}^2}+Q_1{\widehat{u}}_0{}^2,$$

(35)

and

$${\epsilon }^2:\hspace{1em}{D}_0{}^2{\widehat{u}}_2+{\omega }^2{\widehat{u}}_2=-2D_1{D}_0{\widehat{u}}_1-\tilde{\mu }{D}_0{\widehat{u}}_0-D_1{}^2{\widehat{u}}_0+2{\widehat{u}}_1\omega \widehat{\sigma }-{\widehat{\sigma }}^2{\widehat{u}}_0-2D_2{D}_0{\widehat{u}}_0+2Q_1{\widehat{u}}_0{\widehat{u}}_1+Q_2{\widehat{u}}_0{}^3+\frac{4\beta \tilde{\gamma }\sin \omega T}{{(1\mp u_c)}^3}{\widehat{u}}_0.$$

(36)

One can set the solution of Equation (34) as:

$${\widehat{u}}_0=B_1(T_1)e^{i\omega T_0}+{\overline{B}}_1(T_1)e^{-i\omega T_0}.$$

(37)

Substituting Equation (37) into Equation (35), and eliminating the secular terms of Equation (35), one can obtain that:

$$\begin{array}{l}{D}_1{B}_1=\frac{-\beta \tilde{\gamma }}{2\omega {(1\mp u_c)}^2}-i\sigma B_1,\\ {\widehat{u}}_1=-\frac{B_1^{2}P}{3{\omega }^2}{e}^{i2\omega T_0}-\frac{{\overline{B}}_1^{2}P}{3{\omega }^2}{e}^{-i2\omega T_0}+\frac{2B_1{\overline{B}}_{1}P}{{\omega }^2}.\end{array}$$

(38)

Substituting the equation above into Equation (36) and eliminating its secular terms, one can have:

$$D_2{B}_1=-\frac{\tilde{\mu }{B}_1}{2}-\frac{\widehat{\sigma }\beta \tilde{\gamma }}{4{\omega }^2{(1\mp u_c)}^2}-i(\frac{5Q_1^2}{3{\omega }^3}+\frac{3Q_2}{2\omega })B_1^2{\overline{B}}_1.$$

(39)

Now setting

$$B_1=\frac{1}{2}\epsilon b(T_1,T_2)e^{i\phi (T_1,T_2)},$$

(40)

considering

$${\dot{B}}_1\approx D_0{B}_1+\epsilon D_1{B}_1+{\epsilon }^2{D}_2{B}_1,$$

(41)

and substituting Equations (38) and (39) into Equation (41), and expressing Equation (41) by the original dimensionless parameters of Equation (3), one has:

$$\begin{array}{l}\dot{b}=\frac{-(3\omega -\widehat{\omega })\beta \gamma \cos \phi }{2{\omega }^2{(1\mp u_c)}^2}-\frac{\mu b}{2},\\ b\dot{\phi }=\frac{(3\omega -\widehat{\omega })\beta \gamma \sin \phi }{2{\omega }^2{(1\mp u_c)}^2}-(\omega -\widehat{\omega })b-\frac{5b^3{Q}_1^2}{12{\omega }^3}-\frac{3b^3{Q}_2}{8\omega }.\end{array}$$

(42)

The periodic solution of the system (3) satisfies $\dot{b}=0,$ and $\dot{\phi }=0,$ i.e.,

$$-\frac{\mu b}{2}=\frac{(3\omega -\widehat{\omega })\beta \gamma \cos \phi }{2{\omega }^2{(1\mp u_c)}^2},\hspace{1em}(\omega -\widehat{\omega })b+(\frac{5Q_1^2}{12{\omega }^3}+\frac{3Q_2}{8\omega })b^3=\frac{(3\omega -\widehat{\omega })\beta \gamma \sin \phi }{2{\omega }^2{(1\mp u_c)}^2}.$$

(43)

The periodic solution can be expressed analytically as:

$$u=\pm u_c+\frac{2b^2{Q}_1}{3{\omega }^2}+b\cos (\omega T+\phi )-\frac{b^2{Q}_1}{3{\omega }^2}{\cos }^2(\omega T+\phi ).$$

(44)

According to the characteristic solutions of Equation (42), it shows that the theoretical periodic solution expressed by Equation (44) will become unstable if:

$$(\omega -\widehat{\omega }-(\frac{5Q_1^2}{12{\omega }^3}+\frac{9Q_2}{8\omega })b^2)(\omega -\widehat{\omega }+(\frac{5Q_1^2}{12{\omega }^3}+\frac{3Q_2}{8\omega })b^2)\ge \frac{{\mu }^2}{4}.$$

(45)

Based on Equations (43)–(45), the evolution of the periodic solutions of system (3) with the amplitude of AC voltage when $\omega$ = 0.6 is shown in Figure 6. Obviously, when V_AC increases from 0, the two non-trivial equilibria lose their stability; instead, there are two periodic attractors coexisting. The amplitudes of the two periodic attractors increases with the amplitude of AC voltage. The coexistence of multiple periodic attractors can be attributed to the disturbance of the bistable non-trivial equilibria of the system (3) when V_AC = 0 V.

In Figure 6, when V_AC varies from 0 to 0.055 V, the numerical simulation is in great agreement with the theoretical solution. However, when V_AC exceeds 0.056 V, the theoretical prediction of the periodic attractor in the neighborhood of the right non-trivial equilibria is not that accurate, which may be due to the limitation of the Method of Multiple Scale. It will then be essential for us to apply numerical simulation to investigate the evolution of the attractors with the change in AC voltage. The basic settings for the simulation, such as the time step and initial plane, are the same as that in Section 3. The change of the attractors and the area and nature of their basins of attraction with V_AC are shown in Figure 7, where $\omega$ = 0.6. The evolution of global dynamics of system (3) with the increase in V_AC can be separated into the following five stages.

Firstly, when V_AC = 0 V, there are two point attractors coexisting whose basins of attraction are fractal and trigger each other (see Figure 7a1,a2). According to Figure 7a2, in a small neighborhood of each point attractor, the attractor of system (3) is locally stable. Otherwise, a small disturbance of initial conditions will lead to a different point attractor, meaning that it is easy to induce a safe jump.

Secondly, when V_AC increases from 0 to 0.01 V (see Figure 7b1–d2), the number of the periodic attractors increases with V_AC. At V_AC = 0.005 V, the two point attractors become two periodic attractors; apart from these two periodic attractors predicted theoretically, a new periodic attractor appears suddenly, marked by the yellow curve in Figure 7b1, and its basin of attraction is discrete (see the yellow regions in Figure 7b2). It shows that the new periodic attractor is a hidden attractor [26]. When V_AC increases to 0.006 V, another hidden attractor appears, which is almost symmetric to the former one (see the blue curve of Figure 7c1 and the blue regions of Figure 7c2). When V_AC = 0.01 V, there are five periodic attractors coexisting, as shown in Figure 7d1. A new periodic attractor appears (see the green curve in Figure 7d1), whose amplitude is much bigger than the other ones.

Thirdly, as V_AC increases from 0.01 V to 0.116 V, the number of attractors will decrease. Comparing Figure 7e1 with Figure 7d1, it is obvious that when V_AC increases to 0.02 V, the yellow periodic attractor disappears whose basin of attraction is eroded by that of the green attractor; thus, the basin of attraction of the green attractor can be much bigger in Figure 7e2 than in Figure 7d2. When V_AC continues to increase, the other three periodic attractors, i.e., the blue attractor, the red one, and the black one, disappear successively (see Figure 7f1,h1,j1) whose basins of attraction are aggressed by the basin of attraction of the green attractor, as shown in Figure 7e1–j2. Till V_AC becomes 0.116 V, there will be a single periodic attractor left whose basin of attraction is not fractal but with a smooth boundary (see Figure 7j1,j2).

Besides, when V_AC increases to 0.128 V, there will be a new complex attractor coexisting with the former green periodic attractor. It is a period-3 attractor (see the purple curve in Figure 7k1) whose basin of attraction is fractal and eroded to the basin of attraction of the periodic attractor (see Figure 7k2). It follows that a small change of initial conditions possibly shifts the dynamical behavior of the system (3) from a periodic motion to a period-3 motion, which is another type of safe jump.

Finally, as V_AC continues to increase, another type of complex dynamical behavior is induced. According to the phase map, Poincare map, and frequency spectrum in Figure 8a–c, there is only a chaotic attractor when V_AC = 0.28 V, and the boundary of its basin of attraction is not fractal (see Figure 8d).

5. Conclusions

In this paper, a typical electrostatic bilateral micro-resonator is considered. The theory of local bifurcation and numerical approaches are applied to analyze the global dynamics of the vibrating system of the micro resonator. The main conclusions are presented as follows:

(1)

DC bias voltage has some effect on the dynamics of the micro resonator. Without AC voltage, when the DC bias voltage is low, there will be only one stable point attractor in its vibrating system; when the DC bias voltage increases, there may be bistable point attractors.

(2)

In the case of a low bias DC voltage, multiple periodic attractors and the corresponding safe jump occur due to Hopf bifurcation when varying the frequency or amplitude of AC voltage in certain ranges.

(3)

Under a higher bias DC voltage that can induce bistable point attractors, when increasing the value of the amplitude of AC voltage, there will be multiple periodic attractors attributed to the loss of stability of the two non-trivial point attractors; apart from this, there will be some other complex dynamical behaviors of the micro-resonator vibrating system, such as safe jump, hidden attractors, period-n attractor, and chaos.

Our results provide some theoretical reference in avoiding complex dynamics of micro resonators, thus having some potential values in the design of micro sensors. The hidden attractors are depicted numerically, but their mechanism is still not that clear, which will be discussed in our future study.