The Nonlinear Dynamics of a MEMS Resonator with a Triangular Tuning Comb

The nonlinear dynamic response of a MEMS resonator with a triangular tuning comb is studied. The motion equation with dis-smooth tuning electrostatic force is derived according to Newton’s second law. The analytical solution of the periodic response is obtained using the harmonic balance method and section integral method. The singularity theory is then applied to investigate the bifurcation of the periodic response of the untuned system. The transition sets on the DC-AC voltage plane dividing the planes into several persistent regions are obtained. The bifurcation diagrams’ topological structures and jump phenomena corresponding to different parameter regions are analyzed. We explore the effects of tuning voltage on the response. This demonstrates that the amplitude–frequency curves present more hardening characteristics with increased tuning voltage. Many twists, bifurcation points, and unstable solutions appear, leading to complicated jump phenomena. Two bifurcation points exist on the response curves: the smooth and dis-smooth bifurcation points, with the latter occurring on the switching plane of non-uniform fingers.


Introduction
MEMS comb resonators are widely used in signal processing, gyroscopes, oscillators, etc., due to their advantages in terms of signal-to-noise ratio, quality factor, and sensitivity [1].Most MEMS comb resonators typically feature fingers of equal cross-section and length.Recently, researchers have explored alternatives, including combs with non-uniform finger lengths which can produce a displacement-dependent electrostatic force.This has gained significant attention due to its potential applications in modifying the stiffness and resonance frequency of resonators.
Lee et al. [2] first presented a triangular comb array of linearly varied finger lengths for frequency tuning, and reduced the resonant frequency by changing the tuning voltage.Dai et al. [3] investigated the fabrication of a micromechanical tunable resonator.The tuning unit has a movable comb with constant finger length and a fixed comb with linearly varied finger length.Kao et al. [4] presented the design and fabrication of a micromechanical tunable in-plane resonator with non-uniform finger lengths.The effects of the geometric shape of the tuning part and the tuning voltage-to-stiffness ratio on the tuned frequency were discussed.Shmulevich et al. [5,6] believed that the comb drives with tapered fingers produce an electrostatic force proportional to the applied voltage and the motion.They used the device as an electrostatic anti-spring for frequency tuning of the first instability window of a MEMS Meissner parametric resonator.Kavitha et al. [7] presented a design of a MEMS accelerometer with varied finger lengths for sensing low-frequency low-amplitude accelerations.Taherian et al. [1,8] proposed a MEMS comb tunable resonator featuring a triangular comb with non-uniform finger lengths to expand the tuning resonant frequency range while minimizing the required tuning voltage.These studies have been conducted on the design and fabrication of a resonator with combs with variable finger lengths and the variability of the tuned frequency.However, few studies have focused on the equally important aspect of the resonator's response.
The response of the MEMS resonator can be greatly affected by nonlinearity factors.The nonlinearity of the suspension and electrostatic force may result in specific response characteristics, such as hardening/softening/mixture behaviors and pull-in.These characteristics could eventually impact the resonator's performance [9].Therefore, much work has been done on the nonlinear dynamics of the MEMS comb resonator.Nguyen et al. [10] reported a monolithic comb drive oscillator with high-Q and modeled its operation and performance.Jeong et al. [11] showed a theoretical method to obtain the maximum linear displacement of the actuating comb using a nonlinear dynamic model.Zhang et al. [12] studied a MEMS resonator's nonlinear responses and dynamics under two-frequency parametric and external excitations.The effect of varying the DC bias, the squeeze film damping, cubic stiffness, and AC excitation amplitude on the frequency response curves, resonant frequencies, and nonlinear dynamic characteristics, were given.Elshurafa et al. [13] discussed the softening and hardening phenomena of a resonator.They derived the electrostatic force considering both the transverse and longitudinal capacitances of the combs and gave the linear and nonlinear spring constants of the support beams by solving a nonlinear boundary value problem.Tusset et al. [14] studied the chaotic behaviors of the MEMS resonator in a fractional order by historical time and phase portraits.Khirallah [15] added a new electrode that generates an electrostatic force on the truss for the parametric excitation, parametric amplification, and linear and nonlinear tuning of a folded-beam comb drive oscillator.Zhong et al. [16] considered the effects of the inclination of the comb fingers caused by fabrication error.Han et al. [17] investigated the primary resonance of a folded-MEMS comb drive resonator and presented a zoning diagram depicting the extreme vibration amplitude versus DC voltage regarding the hardening and softening bending of response branches and the ideal estimation of dynamic pull-in instability.Zhang et al. [9] investigated the bifurcation characteristics of a resonator and gave the transition sets of different bifurcation behaviors on the parameter plane.The jump phenomena and the influences of the parameters on the amplitude-frequency response were discussed.By disconnecting the rotor and keeping it electrostatically floating, Kassie et al. [18] turned the common resonator into a parametric resonator and gave theoretical and experimental explanations for the considerably large response.Ramanan et al. [19] derived a nonlinear dynamic model with a damping constant obtained from a Monte Carlo simulation and described the vibration responses of micro-resonators operating in the nonlinear region.
In these investigations, researchers typically considered the cubic nonlinear restoring force of the suspension structure and the fractional electrostatic force.However, in resonators with non-uniform finger lengths, there is an additional dis-smooth component of the electrostatic force.Previous studies on the nonlinear dynamics of MEMS resonators have not fully considered this dis-smooth electrostatic force despite its presence in variablelength comb resonators.In refs.[2][3][4][5][6][7][8], the dis-smooth electrostatic force is simplified to a linearly varying force with displacement, which can lead to solution errors.Therefore, it is essential to study the nonlinear dynamics of variable-length comb resonators while taking into account the influence of dis-smooth electrostatic forces.
This study examines the nonlinear dynamics of a MEMS resonator that contains a triangular tuning comb.The innovation lies in considering the dis-smooth electrostatic force of the comb with non-uniform finger lengths and the analytical solution of periodic response based on the harmonic balance method and section integration.We also investigate the impact of parameters on response bifurcation.The results provide a theoretical foundation for designing and utilizing MEMS resonators with non-uniform finger lengths.The paper's structure is as follows: In Section 2, the dynamic equations of the MEMS resonator are established based on Newton's second law, and the electrostatic force of the comb with non-uniform finger lengths is derived.In Section 3, the analytical solution of the periodic response is obtained using the harmonic balance method and section integral method.The stability of the periodic solution is analyzed.In Section 4, the effects of the DC and AC voltages on the bifurcation behaviors of the untuned system are discussed using the singularity theory.The dynamic response of the tuning system is analyzed in Section 5. Finally, the conclusion is summarized.

Physical Model of the Resonator and its Mathematical Description
We consider a tunable MEMS resonator, as shown in Figure 1 [1].A nested folded beam suspension supports the proof mass of the resonator.The resonator has a driving part, a sensing part, and a tuning part.The driving and sensing parts consist of movable, fixed combs with constant finger lengths.The tuning comb contains a movable comb of constant finger length and a fixed comb of varied finger length.The DC voltage V D is applied to the movable structure, the AC voltage V S = V A cosωt is applied to the fixed driving comb, and the proof mass is actuated under the action of the electrostatic force generated by the voltages.The sensing part generates variation in capacitance during the operation of the resonator.When applying a voltage V T to the fixed tuning comb, the electrostatic force can tune the resonant frequency and the response.According to Newton's second law, the motion equation of the resonator can be given as follows: where m is the proof mass, k 1 and k 3 are the linear and cubic nonlinear stiffness coefficients, and F d , F s , and F t denote the driving, sensing, and tuning electrostatic force, respectively.Figure 2 shows the local enlargement of (I) in Figure 1 and the deformation of the suspension structure.Assuming that the length of the supporting beams is equal, when the proof mass has a displacement x, the deformation of each supporting beam is According to ref. [13], the restoring force of the structure in Figure 2 can be obtained as follows: where E is the Young's modulus of the resonator, I = hb 3 /12 is the section modulus of the supporting beams, h is the structural thickness, and b and L are the width and length of the beams.Then, the restoring force of the suspension of the resonator is To obtain the expression of the electrostatic forces, the capacitance of the combs should be analyzed.The comb capacitance consists of three parts (see Figure 3a): the capacitance between the fixed finger and the movable finger C f , the capacitance at the tip of the fixed finger C t f , and the capacitance at the tip of the movable finger C tm .For the combs with uniform finger length [13] where w is the finger width, d is the spacing between the fingers, and x 0 and l are the static separation and overlap between the fixed and moving combs. is the dielectric constant.When the proof mass has a displacement x (shown in Figure 3b), the separation and the overlaps between the combs are changed, and the driving and sensing capacitance are where N is the number of fingers on a single side of the driving and sensing combs; thus, the electrostatic forces in the driving and sensing directions can be expressed as The movable fingers of the tuning comb are uniform (see Figure 4).The tip capacitance of the ith movable finger is The length of the fixed fingers is non-uniform, and the separation of each finger is different.Thus, the tip capacitance of the ith fixed finger is where x 0(i) = L − l i , L is the static separation between fixed and movable electrodes, and l i denotes the length of the ith finger of the fixed comb.The capacitance between the fixed and movable fingers depends on the overlap.There is a capacitance when the movable and fixed fingers overlap (see finger i in Figure 4).The capacitance should be zero if the fingers do not overlap (see finger i − 1 in Figure 4).Thus, the capacitance between ith movable and fixed fingers can be given as where Therefore, the tuning capacitance as the proof mass has a displacement x is where Then, the tuning electrostatic force is obtained as follows: Note that l i − x 0 + x = 0 is a critical state, and when x crosses the critical state, F t will change abruptly.Finite element analysis is performed to verify the tuning electrostatic force.Figure 5 presents the electric potential distribution of the tuning electrode.Figure 6 shows the influence of the displacement on the tuning capacitance, where C t denotes the tuning capacitance and x is the displacement.It can be seen that the theoretical result is approximate to the FEA result.The difference may stem from the insufficient consideration of the edge field in the theoretical result.Figure 7 gives the effect of the displacement on the tuning electrostatic force, where the electrostatic force F t is calculated using the differential of the tuning capacitance so that the FEA result appears to be fluctuating, but the sudden jumps in the force are clearly visible.
By setting x = Xx 0 , τ = ωt, the motion equation can be written into the nondimensional form as where where In our study, the fixed tuning comb is set to be triangular.For the given data in Table 1, f (δ) can be rewritten as where j is an integer.Figure 8 shows the f (δ) curve versus the non-dimensional displacement X.Notice that f (δ) is dis-smooth, and the switching plane of f (δ) is j/21.

Analytical Solution of the Periodic Response and its Stability
Periodic response is essential to the performance of resonators.In this section, the harmonic balance method is used to calculate the periodic response of the system.Considering that tuning may cause the equilibrium position of the response to deviate from 0, the solution is assumed as By substituting Equation (20) into Equation ( 16), and extracting the coefficients of cos τ, sin τ, and the constant term through Fourier expansion, the equations about the amplitudes A 0 , A 1 and phase angular θ can be given as [20] Except for F T0 , F T1 , F T2 , other terms can be obtained through the integral of residue theorem [21].The detailed calculation process and results are shown in Appendix A.
Figure 9 shows the schematic diagram of the integral of f (δ).Considering f (δ) is non-smooth, its integral must be carried out in sections equal to each section's sum.Since F T ≥ 0, A 0 ≥ 0, the maximum value for displacement X is A 0 + A 1 , and the minimum value is A 0 − A 1 .According to the expression of f (δ), we set where indicates rounding down to an integer and denotes rounding up.Obviously, when X < n 2 21 or X ≥ n 1 21 , f (δ) is constant, thus where For n 2 21 < X ≤ n 1 21 , the switching plane divides the response into n 1 − n 2 intervals.The response curves within each interval are similar, so there is where Similarly, F T1 and F T2 can be obtained.The integral results are given in Appendix B. Figure 10 gives the comparisons between the analytical solutions and the numerical solutions, where τ is the non-dimensional time, X is the non-dimensional displacement, HB denotes the solution using the harmonic balance method in this paper, and RK denotes the numerical solutions obtained by the fourth-order Runge-Kutta method.It can be observed that the analytical solutions are in good agreement with the numerical solutions.Figure 10.The comparisons between the analytical and numerical solutions, where τ is the nondimensional time, X is the non-dimensional displacement, HB is the analytical solution using the harmonic balance method, and RK denotes the numerical solutions obtained using the Runge-Kutta method.(a) The Floquet theory can determine the stability of periodic solutions.The motion equation is rewritten into the state equation as where ) is the smooth part of the electrostatic force, A is the derivative matrix of the smooth part, and B is the vector of the non-smooth part, where As shown in Figure 11, the phase trajectory of the periodic solution moving counterclockwise is divided into n 1 − n 2 + 3 sections by switching planes Σ and point P, where Σ j denotes the switching plane at X = j/21 and point P corresponds to τ = 0 and 2π.We define the section of phase trajectory from the plane Σ j−1 to Σ j as C j,+ , and the section from the plane Σ j to Σ j−1 as C j,− .The plane Σ j through which the phase trajectory crosses from left to right is Σ j,+ , and vice versa is Σ j,− .τ j,+ and τ j−1,+ correspond to the intersection of phase trajectory cross Σ j,+ and Σ j−1,+ .It is noted that, corresponding to the section of the phase trajectory between two planes, f (δ) is constant and the system is smooth.The system is dis-smooth only on the switching plane f (δ) changes.The Floquet multiplier matrix of periodic response can be determined via the continuous product of basic solution matrices of all smooth sections and switching matrices on the plane as where M and S are the basic solution matrix and switching matrix corresponding to the corresponding section C and plane Σ, respectively.M S and M E are the starting and ending basic solution matrices between P and the nearest switching plane.The basic solution matrix M is obtained using the following method [22].

Phase trajectory
Switching plane The schematic diagram for calculating the Floquet multiplier matrix of periodic response, where Σ denotes the switching plane, C is the section of the phase trajectory, P is the starting point (τ = 0) and the ending point (τ = 2π), and M and S are the corresponding basic solution matrice of the section and switching matrice of the switching plane.
To calculate M j,+ , we divide the non-dimensional time interval [τ j−1,+ , τ j,+ ] into Ñ equal parts.The interval is C ñ is the approximate matrix of the periodic variable coefficient matrix A in the ñ th time interval, where The basic solution matrix can be given as In the same way, we can determine M S , M E , M j,− .The switching matrix S can be obtained using [23] where n T is the normal vector of the switching surface; here, we have n T = [1, 0] T .F + (Y, τ) and F − (Y, τ) denote the vector fields behind and before the switching plane, respectively.For Σ j,+ , there is And for Σ j,− , it is the opposite.

Bifurcation Analyses of the Resonator System without Tuning
In this section, the singularity theory is used to study the influence of the parameters on the bifurcation behaviors of the untuned resonator.When the tuning voltage equals the DC voltage, A 0 is very small, so we set f 0 = 0 and A 0 = 0 in Equations ( 22) and (23).By defining A 1 = u, θ = v, ω = λ, Equations ( 22) and ( 23) become The transition sets of Equation ( 38) are defined as D = B ∪ H ∪ DL, where B, H, and DL denote the bifurcation set, the hysteresis set, and the double limit point set [24].The expressions of the transition sets are given as follows.
Bifurcation set: Hysteresis set: Double limit set: Figure 12 shows the transition sets on the V D − V A plane with different quality factors where V D is the DC voltage, V A is the amplitude of the AC voltage, Q is the quality factor of the resonator, and BS, HS, and DLS denote the bifurcation set, hysteresis set, and double limit point set.The transition sets dividing the planes into several persistent regions become complicated as the quality factor increases.The number of persistent regions also increases with the quality factor.Figure 13 provides the bifurcation diagrams for different parameter regions.To clarify the topological structure of the bifurcation curves, we have identified four key points (a through d) that mark the turning points of the curves.
The bifurcation diagrams have a single solution branch in all regions except E, F, and G. Region A does not have a turning point, and the bifurcation curve has a single solution for any ω value.Regions B and C have two turning points, with three solutions between the key points.Four key points appear in regions D, H, and I.In region D, point d is to the right of b, and there is a single solution region between points b and d.In region H, point d is to the left of b, resulting in four solutions between b and d.In region I, point d is to the left of a. Regions E, F, and G have bifurcation diagrams with two solution branches.Only one turning point exists in region E. In region F, point d is to point b's right, while region G is to the left of b.Notably, the curves display hardening in region B, softening in C and E, and a combination of both in regions D and F through I.In Figure 14, we have provided the frequency-amplitude curves for various regions, where f = ω/2π is the driving frequency, A is the non-dimensional amplitude, SS is the stable analytical solution, US is the unstable analytical solution, BP is the bifurcation point, NSFI denotes the numerical solution as the driving frequency increases, and NSFD denotes the numerical solution as the driving frequency decreases.Table 2 contains the corresponding parameter values.Different types of jump phenomena are observed in these regions.Specifically, no jump occurs in region A. In regions B and C, the response jumps once as the frequency increases and decreases.In regions F to I, jump phenomena occur twice with the increase of driving frequency and once with the decrease.In region D, the jump will occur twice when the frequency increases and decreases.However, the response will not jump down in regions E and F when the frequency is reduced.The numerical results presented in Figure 14 align well with the analytical results.

The Effects of Tuning on the Response of the Resonator
For this section, the quality factor Q is set to 20,000.We identify four groups of DC and AC voltages that correspond to amplitude-frequency curves of the untuned system with different topologies and study the effect of tuning on the response.Figure 15 shows the results of the tuning voltage on the amplitude-frequency curve of the system, where the tuning voltage is defined as V TA = |V D − V T |, and the amplitude A m = A 0 + A 1 is the maximum displacement of the response, US is the unstable analytical solution, and BP is the bifurcation point.where US is the unstable analytical solution, BP is the bifurcation point, (a) As the tuning voltage increases, the amplitude-frequency curve shifts leftwards.The amplitude-frequency curves in Figure 15a,c change from non-skewed and softening to a mixture of hardening and softening characteristics, respectively.Additionally, the tuned response curves exhibit many twists, bifurcation points, and unstable solutions, making the amplitude-frequency curve complex.Figure 16 displays the peak frequency curve as the tuning voltage varies, showing that the peak frequency drops faster with a larger tuning voltage, where f P denotes the peak frequency of the amplitude-frequency curve.
Figure 17 gives the amplitude-frequency curve and its local enlargement, with the ordinate magnified 21 times, where SS is the stable analytical solution, US is the unstable analytical solution, BP is the bifurcation point, NSFI denotes the numerical solution as driving frequency is increased, NSFD denotes the numerical solution as driving frequency is decreased, and the yellow lines denote the switching planes.The numerical results confirm the analytical solution.Two types of bifurcation points exist on the curve: smooth (P 1 in Figure 17c and P 4 in Figure 17e) and dis-smooth (P 2 in Figure 17d and P 3 in Figure 17e).The amplitude of all dis-smooth bifurcation points is an integer multiple of 1/21, indicating that they occur on the switching plane of non-uniform fingers.As the frequency changes, the system response undergoes complicated jumps.In Figure 17c,d, the response jumps down at point P 1 as the frequency increases and jumps up at point P 2 when the frequency decreases.In Figure 17e, the response does not reach the nearest solution after jumping up at point P 3 as the frequency increases but rather jumps to a solution with a larger amplitude.Considering that there is a stable solution and an unstable solution between every two switching planes, this provided the possibility to achieve the desired periodic response between the switching planes by designing the length of the fingers reasonably, and quickly switching between multiple periodic responses by jumping.However, operating the tuning resonators in the open-loop mode is not advisable, as the response can easily jump to other solutions when disturbed, resulting in an incorrect output.

Conclusions
This paper studies the nonlinear vibration response of a MEMS resonator with a triangular tuning comb.According to Newton's second law, the dynamic equation of the resonator is derived, and it is found that the tuned electrostatic force of the triangular tuning comb is dis-smooth.The harmonic balance and section integral methods are used to obtain the analytical solution of the dis-smooth system.The stability of the periodic response was studied based on the Floquet theory.The numerical results show an excellent agreement with the analytical solutions.
The singularity theory is used to investigate the bifurcation of the periodic response of the untuned system.The transition sets on the DC-AC voltage plane corresponding to different quality factors are obtained, which divide the planes into several persistent regions.The topological structures of the bifurcation diagrams corresponding to different parameter regions are analyzed.Abundant jump phenomena of the periodic response, varying with the driving frequency, are present.
The effects of tuning voltage on the system response are studied.We demonstrated that, with the increase of the tuning voltage, the amplitude-frequency curve moves to the left and presents more hardening characteristics.Many twists, bifurcation points, and unstable solutions appear on the curve, which leads to complicated jump phenomena.Two bifurcation points exist, the smooth bifurcation points and the dis-smooth bifurcation points, where the dis-smooth bifurcation points occur on the switching plane of the nonuniform comb.

Figure 2 .
Figure 2. Local enlargement of (I) in Figure 1 and the deformation of the suspension structure when the proof mass has a displacement x.

Figure 3 .
Figure 3.The schematic diagram of the structure, capacitance, and operation of the combs with uniform finger length.(a) The structure and capacitance of the combs; (b) the operation of the driving and sensing combs.

Figure 4 .
Figure 4.The schematic diagram of the structure and capacitance of the combs with varied finger lengths.

Figure 7 .
Figure 7.The effect of the displacement on the tuning capacitance, where (b) is the local enlargement of (a), F t denotes the tuning electrostatic force, x is the displacement, and V D = 0, V T = 1 V.

Figure 9 .
Figure 9.The schematic diagram of the integral of f (δ).

Figure 12 .
Figure 12.The transition set of Equation (38) on V D − V A plane with different quality factors, where V D is the DC voltage, V A is the amplitude of the AC voltage, Q is the quality factor of the resonator, and BS, HS, and DLS denote the bifurcation set, hysteresis set, and double limit point set, (a) Q = 5000, (b) Q = 10,000, (c) Q = 20,000, (d) Q = 40,000 .

Figure 13 .
Figure 13.Bifurcation diagrams corresponding to different persist regions of Figure 12, where (A-I) correspond to regions A through I, respectively.

Figure 14 .
Figure 14.Amplitude-frequency curves corresponding to different regions of Figure8, where (A-I) correspond to regions A through I, respectively.SS is the stable analytical solution, US is the unstable analytical solution, BP is the bifurcation point, NSFI denotes the numerical solution as the driving frequency increases, and NSFD denotes the numerical solution as the driving frequency decreases.

Figure 15 .
Figure 15.The effects of the tuning voltage on the amplitude-frequency curves of the resonator, where US is the unstable analytical solution, BP is the bifurcation point, (a)V D = 1 V, V A = 10 mV, (b) V D = 2 V, V A = 15 mV, (c) V D = 5 V, V A = 8 mV, and (d) V D = 2.6 V, V A = 30 mV.

Figure 16 .
Figure 16.The curve of the peak frequency varying with the tuning voltage.

Figure 17 .
Figure17.The amplitude-frequency curves and those of the resonator under tuning.SS is the stable analytical solution, US is the unstable analytical solution, BP is the bifurcation point, NSFI denotes the numerical solution as driving frequency is increased, NSFD denotes the numerical solution as driving frequency is decreased, and the yellow lines denote the switching planes.(a)V D = 1 V, V A = 10 mV, (b) V D = 5 V, V A = 8 mV, (c,d) are the local enlargements of (a), (e,f) are the local enlargements of (b).

Table 2 .
Values of parameters corresponding to Figure 14A-I.