## Abstract

In this paper, the presence of *U*(Universal)-sequence (a sequence of periodic windows that appear beyond the period doubling (PD) route to chaos) in electrostatic microelectromechanical systems (MEMS) is reported. The MEM system is first brought to a nonlinear steady state by the application of a large dc bias close to the dynamic pull-in voltage of the device. An ac voltage (the bifurcation parameter) is next applied to the system and increased gradually. A sequence of PD bifurcations leading to chaos is observed for resonant and superharmonic excitations (frequency of the ac voltage). On further increase in the ac voltage (beyond where chaos sets in), *U*-sequence is observed in the system. Under superharmonic excitation, the sequence is found to be a modified form of the *U*-sequence referred to as the ‘*UM*-sequence’ in this paper. The appearance of a periodic window with *K* oscillations per period or *K*-cycles in the normal *U*-sequence is replaced by a corresponding periodic window with *KM*-cycles in the *UM*-sequence. *M* stands for the *M*th superharmonic frequency of excitation. The formation of the periodic windows from a chaotic state in the *UM*-sequence takes place through intermittent chaos as the ac voltage is gradually increased. On the other hand, the periodic states/cycles formed through intermittent chaos transform back into a chaotic state through the period doubling route. A sequence of period doubling bifurcations of the *UM*-sequence cycles result in the formation of -cycles in electrostatic MEMS. *n* corresponds to the *n*th period doubling bifurcation in the sequence. A simplified mass–spring–damper (MSD) model for MEMS is used to understand the physical mechanism that gives rise to these nonlinear dynamic properties in MEMS. The nonlinear nature of the electrostatic force acting on the MEM device is found to be responsible for the reported observations.

## 1. Introduction

The nonlinear coupling between the displacement of the microstructure and the electrostatic force acting on it is an inherent source of nonlinearity in electrostatic microelectromechanical systems (MEMS) (El-Bassiouny & Eissa 2003). This inherent nonlinearity has been shown to be responsible for various interesting dynamic properties in electrostatic MEMS like jump phenomena (El-Bassiouny & Eissa 2003), spring hardening/softening (Tilmans & Legtenberg 1994), hysteresis in resonant strain gauges (Gui *et al*. 1998), multiple resonant peaks (De & Aluru 2004) and parametric resonance in microelectromechanical probes (Turner *et al*. 1998). The existence of the period doubling (PD) route to chaos at resonant excitation (Wang *et al*. 1998; Liu *et al*. 2003) and at superharmonic excitation (De & Aluru 2005) has been observed experimentally/theoretically due to the nonlinear nature of the electrostatic force. These nonlinear dynamic properties and chaos in electrostatic MEMS can be exploited for various applications like MEMS filters with shiftable resonant frequencies (Adams *et al*. 1996), chaotic microfluidic mixers (Lee *et al*. 2001), and secure (chaotic) communications (Wang *et al*. 1998). As a result, it is important to thoroughly understand and analyse the nonlinear/chaotic regions in MEMS for better and efficient applications. It has been shown by the authors in De & Aluru (2005) (and more details in De & Aluru (2006)) using numerical simulations, that under normal atmospheric conditions, as the dc bias applied to a MEM device is increased in the absence of any ac voltage, a characteristic change in the phase plot of the system termed as the ‘dc symmetry breaking’ takes place close to the dynamic pull-in voltage of the system. Once a steady state was reached due to a dc bias exhibiting dc symmetry breaking, an ac voltage was applied at the *M*th superharmonic frequency of excitation. For an ac voltage around ‘ac symmetry breaking’ (De & Aluru 2005, 2006; another symmetry breaking in the phase plot of the system), *M* oscillations per period or *M*-cycles were observed in the MEM device. On further increasing the amplitude of the ac voltage, PD bifurcations (Thompson & Stewart 1986) took place resulting in the formation of -cycles and ultimately chaos. *n* stands for the *n*th PD bifurcation in the sequence and *M*=1 implies resonant excitation.

In this paper, the chaotic region formed after the sequence of PD bifurcations in electrostatic MEMS is investigated for both resonant and superharmonic excitations for the presence of periodic windows. It is often desirable to establish non-chaotic behaviour in the vicinity of parameter values (the ac voltage in this case) that give rise to chaos (Barreto *et al*. 1997), making it important to analyse the periodic windows that exist in the chaotic region of a dynamical system. On gradually increasing the ac voltage in the chaotic region at resonant excitation, the *U*(Universal)-sequence (Thompson & Stewart 1986) is observed and is reported for the first time in this paper for electrostatic MEMS. Periodic windows suddenly appear and then disappear within the chaotic region following a definite sequence (popularly known as the *U*-sequence). A modified form of the *U*-sequence is found to exist in the chaotic region of the MEM system under superharmonic excitations and is also reported in this paper for the first time. The appearance of a periodic window with *K*-cycles in the normal *U*-sequence is replaced by the appearance of a periodic window with *KM*-cycles in the modified sequence termed as the ‘*UM*-sequence’, for an exciting frequency of (*M*th superharmonic frequency). is the resonant frequency of the device. The appearance of the periodic windows (i.e. the transition from chaos to the periodic states) in the *U*- and *UM*-sequences is found to take place through intermittent chaos (Thompson & Stewart 1986). The density of the chaotic regions in the time-series plot for the MEM device gradually decreases and the density of the periodic regions gradually increases as the ac voltage is increased from that at a chaotic state. Finally, the formation of the periodic window is marked by the total disappearance of the chaotic regions in the time-series plot. On the other hand, the disappearance of the periodic windows (i.e. the transition from the periodic states to chaos) takes place through a sequence of PD bifurcations. A *KM*-cycle oscillation gives rise to -cycles due to PD bifurcations and ultimately goes to chaos . In this paper, these observations are made in two different MEM devices namely, a MEMS fixed–fixed beam and a MEMS torsion micromirror, using numerical simulations.

The rest of the paper is outlined as follows: §2 presents the theory of MEMS dynamics, where physical models for the different energy domains present in MEMS, the computational procedure used for solving the physical models and a simpler mass–spring–damper (MSD) model for MEMS are discussed. Section 3 presents numerical simulation results (using the physical models) and §4 presents a discussion (using the physical/MSD models) on the new nonlinear dynamic properties of electrostatic MEMS reported in this paper. Conclusions are presented in §5.

## 2. Theory of MEMS dynamics

Electrostatic MEMS dynamics is governed by coupled electrical, mechanical and fluidic (electromechanical–fluidic) energy domains. Figure 1 shows a typical MEM device—a deformable fixed–fixed beam over a fixed ground electrode. This device and the other MEM device considered in this paper, a micromirror, are characterized by a thin film of air in between the electrodes. We will explain the coupling among the three energy domains by considering the example shown in figure 1, but the discussion is applicable to the micromirror and other MEM devices as well. Considering figure 1, at *t*=0, a potential difference *V* is applied between the two conductors (the beam and the ground electrode) which induces electrostatic charges on the surface of the conductors. The distribution of electrostatic charges on the surface of the conductors depends on the relative position of the two conductors. These electrostatic charges give rise to an electrostatic pressure, which acts normal to the surface of the conductor as shown in figure 1*a*. The ambient fluid pressure also acts on the beam as shown in figure 1*a* at *t*=0. When the beam deforms due to the electrostatic pressure, the charge redistributes on the surface of the conductors and consequently, the electrostatic pressure field changes. At the same time, the displacement of the fluid/air in between the electrodes due to the deformation of the beam gives rise to squeeze film damping. The electrostatic pressure and the squeeze film damping pressure acting on the lower side of the beam and the ambient fluid pressure acting on the top surface of the beam together cause the beam to deform to a state where they are balanced by the internal stiffness and the inertial forces at that time instant (see figure 1*b*). The stiffness force depends on the displacement of the beam at the given time instant and on the material properties of the beam. The inertial force depends on the acceleration of the beam at that time instant and on the density of the beam. The fluid pressure due to squeeze film damping depends on the velocity and position of the beam with respect to the ground electrode at the time instant. In summary, at each time instant, the inertial, stiffness, electrostatic and the fluidic damping forces need to be computed and the position of the beam is obtained by considering the balance of all the forces. The physical models to compute the various forces, the computational procedure for solving the physical models and a MSD model for MEMS are described next. The MSD model is used for the analysis of the nonlinear dynamic simulation results (obtained using the physical models).

### (a) Physical models

The physical models for coupled electromechanical–fluidic analysis of MEMS dynamics are presented in this section. The mechanical deformation of the MEM structure can be computed by a two-dimensional large deformation elastic analysis of the microstructure. The transient governing equations for an elastic body using a Lagrangian description are given by (Chandrasekharaiah & Debnath 1994)(2.1)(2.2)(2.3)(2.4)(2.5)Equation (2.1) is the governing equation where *ρ* is the material density in the undeformed (initial) configuration, is the deformation gradient, and is the second Piola–Kirchhoff stress. , and are the displacement, velocity and acceleration vectors, respectively. Equations (2.2) and (2.3) are the displacement and the surface traction boundary conditions, respectively. is the prescribed displacement and is the unit outward normal vector in the initial configuration. is the surface traction on the structure due to the electrostatic and fluid pressures and is the first Piola–Kirchhoff stress tensor. Equations (2.4) and (2.5) are the initial conditions for displacement and velocity, respectively. and are the initial displacement and velocity, respectively. A Newmark scheme with an implicit trapezoidal rule (see Bathe 1995 for details) is used to solve the nonlinear dynamical system posed in equations (2.1)–(2.5). Numerical discretization is done by using the finite cloud method (FCM) (see Aluru & Li 2001; Jin *et al*. 2004, 2005 for details on FCM and the accuracy/convergence of the method).

The two-dimensional governing equation for the electrostatic analysis can be written in a boundary integral form in the Lagrangian frame (see Li & Aluru 2002*a* for details)(2.6)(2.7)(2.8)where *ϵ* is the dielectric constant of the medium, *ϕ* is the electrostatic potential and *σ* is the electrostatic surface charge density. *P* and *Q* are the source and field points in the initial configuration corresponding to the source and field points *p* and *q* in the deformed configuration, and *G* is the Green's function. In two dimensions, , where is the distance between the source point and the field point . is the total charge of the system and *C* is an unknown variable which can be used to compute the potential at infinity. is the tangential unit vector at field point *Q* and is the Green deformation tensor. Equations (2.6)–(2.8) are solved using the boundary cloud method (BCM) (see Li & Aluru 2002*b*, 2003 for details on BCM and the accuracy/convergence of the method) to obtain the distribution of surface charge density *σ* on the conductors. The electrostatic pressure normal to the surface of the conductors is given by(2.9)

Reynold's squeeze film equation (RSFE) is used to compute the fluid/air damping pressure acting on the MEM structure. RSFE is applicable for structures where a small gap between two plates opens and closes with respect to time (Starr 1990). This assumption holds for structures where the seismic mass moves perpendicular to a fixed wall, for plates with tilt around horizontal axes and for clamped beams where the flexible part moves against a fixed wall (Hung & Senturia 1999; Schrag & Wachutka 2002) as shown in figure 2. Under the assumption of a small Reynolds number and a sufficiently large ratio of structure length (*l*) to fluid film (air in this case) thickness (*h*), the Navier–Stokes equations can be approximated by the much simpler Reynolds equation (Gross *et al*. 1980). The isothermal RSFE for a compressible slip flow is given by (Burgdorfer 1959; Hung & Senturia 1999)(2.10)where *h* is the gap between the movable structure and the ground electrode of the MEM device (same as the fluid film thickness), is the fluid pressure acting on the structure, and *η* is the viscosity of the surrounding fluid. is the Knudsen number, where *λ* is the mean free path of the surrounding fluid. The mean free path *λ* is related to the ambient temperature and pressure by the relation (Vincenti & Kruger 1967)(2.11)where is the Boltzmann's constant, *T* is the absolute temperature, *p* is the ambient pressure and *d* is the collision diameter of the fluid molecules (*d*=3.66 Å for air (Vincenti & Kruger 1967)). The fluid domain where the RSFE is solved is the projection of the MEM structure on the *X*–*Z* plane (Hung & Senturia 1999) as shown in figure 2. As the moving structure deforms due to the application of an external electric field between it and the ground plane, the projected fluid domain also changes. Hence, a Lagrangian form of equation (2.10) is developed(2.12)where *X*, *Z* are the coordinates of the fluid domain corresponding to the undeformed state of the movable structure and *u* is the deformation of the movable structure in the *X*-direction. As the two-dimensional mechanical equations are solved in the *X*–*Y* domain, mechanical deformation and its variation in the *Z*-direction are assumed to be zero. The fluid pressure obtained from equation (2.12) is integrated along the *Z*-direction to compute an effective fluid pressure , which is applied as a boundary condition in the two-dimensional mechanical analysis in the *X*–*Y* domain (Hung & Senturia 1999).

The effective fluid pressure from the fluidic analysis and the electrostatic pressure obtained from the electrostatic analysis are used to compute in equation (2.3) using(2.13)where . A full-Lagrangian Newton scheme (De & Aluru 2004) is used to obtain a self-consistent solution of the coupled electromechanical–fluidic analysis at each time step. The full-Lagrangian Newton scheme uses a Lagrangian description of all the energy domains and a Newton method for coupling the different energy domains (see De & Aluru 2004 for details). The physical models presented above for the coupled electromechanical–fluidic analysis of MEMS have been validated in De & Aluru (2006) considering different MEM devices.

### (b) Mass–spring–damper model

Although numerical simulations (coupled electromechanical–fluidic analysis) using the physical models for the different energy domains give accurate results (when compared with experimental data), it is easier to understand the physics behind the nonlinear dynamic properties of electrostatic MEMS through a simpler MSD model. A MSD model for electrostatic MEMS dynamics is given by (Liu *et al*. 2003)(2.14)where *m* is the mass of the movable structure, *c* is the squeeze film damping constant, *k* is the stiffness constant of the structure and *x* is the displacement of the movable structure. The analytical expression for the electrostatic force depends on the gap (*g*) in the undeformed state between the microstructure and the ground electrode, the area (*A*) of the microstructure surface facing the ground electrode and the applied voltage (*V*). Nonlinearities in the mechanical stiffness and fluid damping forces (present in the physical models) are neglected in the MSD model. The only source of nonlinearity in the MSD model is from the nonlinear expression of the electrostatic force in terms of the displacement *x*. As a result, the MSD model is primarily used to explain the effect and the significance of the electrostatic force in MEM dynamics in §4.

## 3. Nonlinear dynamic properties of microelectromechanical systems

Numerical simulations using the coupled electromechanical–fluidic physical models are done to present new nonlinear dynamic properties of electrostatic MEMS, namely, the *U*- and *UM*-sequences, intermittent chaos and the formation of -cycles. Two different MEM devices are considered to present the new properties. The first device is a MEMS torsion mirror (Pan *et al*. 1998) shown in figure 3. The mirror plate is 1500 μm long, 1400 μm wide and 3 μm thick. Four 150 μm long electrodes (of negligible thickness) are attached to the mirror plate and the ground plane. The distance from the centre of the electrodes on the mirror to the centre of the mirror is 505 μm. The gap between the mirror and the ground plane is 42 μm. *I*=5.51×10^{−15} kg m^{2} is the moment of inertia of the mirror and *K*=4.49×10^{−7} Nm rad^{−1} is the torsional stiffness of the mirror. Numerical simulations of the MEMS torsion mirror are done by applying a potential (combination of a dc bias and an ac voltage of the form ) to one of the electrodes on the mirror and all the other three electrodes (one on the mirror and two on the ground plane) are grounded as shown in figure 3. The time series, phase and Poincarè plots of the MEMS torsion mirror presented in the following sections are based on the angular rotation *θ*, angular velocity and angular acceleration of the mirror plate (which is assumed to be rigid (Pan *et al*. 1998)) obtained from the numerical simulations at each time step. The second device considered is a silicon fixed–fixed beam 80 μm long, 1 μm thick and 10 μm wide placed 1 μm over a ground plane as shown in figure 4. The Young's modulus of 169 GPa, a mass density of 2231 kg m^{−3} and a Poisson's ratio of 0.3 are employed in the simulations for the silicon beam. The time series, phase and Poincarè plots of the fixed–fixed beam presented in the following sections are based on the peak displacement, velocity and acceleration at the centre of the beam in the negative *Y*-direction (see figure 4) obtained from the numerical simulations at each time step.

Both the devices (the torsion mirror and the fixed–fixed beam) are considered to be in air at 1 atm and 293 K. The viscosity is *η*=1.82×10^{−5} kg ms^{−1} for air and the mean free path *λ* is computed from equation (2.11) for air. In all the examples, numerical discretization errors are minimized by using fine point distributions and small time steps (the adequacy of the point distribution and the time step is checked by the convergence of the solution). The time step used in all the simulations presented here is 1/50*f* (50 simulation points per time period), where *f* is the frequency of the applied voltage or the resonant frequency of the device if the applied voltage is a dc bias. The ac voltage is considered to be the bifurcation parameter for both the devices and is increased gradually in the numerical simulations by passing the end conditions of the present simulation as the initial conditions to the next simulation at a slightly higher ac voltage.

### (a) U-sequence

Numerical simulations of the MEMS torsion mirror show that the mirror exhibits dc symmetry breaking from 133 V to dynamic pull-in at 142 V dc (De & Aluru 2005). An ac voltage of the form is applied to the MEMS torsion mirror at the *M*th superharmonic frequency of excitation, i.e. at , after a steady state is reached due to the application of a dc bias of 141 V (in the dc symmetry breaking region). is the resonant frequency of the device for the applied dc voltage. As the ac voltage (the bifurcation parameter of the system) is increased, ac symmetry breaking (De & Aluru 2005) takes place. Figures 5 and 6 show the phase plots for the mirror for two different ac voltages under the 141 V dc bias at *M*=1 (resonant excitation). 0.5 V ac does not exhibit any ac symmetry breaking (phase plot is symmetric with respect to the two axes), whereas a large 3 V ac (ac pull-in takes place at 3.85 V ac) exhibits ac symmetry breaking (phase plot is not symmetric with respect to the two axes). The resonant frequency is *f*_{0}=425 Hz at 141 V dc for the mirror. On further increasing the ac voltage, a sequence of PD bifurcations takes place in the system giving rise to 2^{n}-cycles and ultimately chaos. Figures 7 and 8 show the phase plot and the Poincarè map, respectively, for the chaotic state reached at 3.7999 V ac after the sequence of PD bifurcations in the MEMS mirror.

The region beyond the PD route to chaos is investigated in this paper by gradually increasing the ac voltage from the chaotic state at 3.7999 V ac. *U*-sequence (a sequence of periodic windows that appear beyond the PD route to chaos (Thompson & Stewart 1986)) is observed in this region for the MEM device. The presence of *U*-sequence in electrostatic MEMS is reported for the first time in this paper. As the ac voltage is increased from 3.7999 V, stable periodic windows are observed for certain values/regions of the ac voltage. For all other values of the ac voltage, the MEM device is found to exhibit chaos. The order in which the periodic windows appear in the *U*-sequence of a dynamical system can be obtained using symbolic dynamics (Robinson 1994) and is given as (considering up to 6 cycles):(3.1)Table 1 shows the *U*-sequence and the corresponding ac voltages at which the periodic states/cycles are observed for the mirror at 141 V dc and *M*=1. Figure 9 shows the phase plot and the Poincarè map for the *K*=6-cycle oscillation (the first periodic window) formed at 3.8005 V ac in the MEMS torsion mirror. The phase plot and the Poincarè map for the subsequent *K*=5-cycle oscillation (the second periodic window) at 3.8036 V ac are shown in figure 10. Figure 11 shows the phase plot and the Poincarè map and figure 12 shows the time-series plot for the *K*=3-cycle oscillation at 3.8085 V ac. From figures 11 and 12, it can be seen that the 3-cycle is characterized by three loops per period in the phase plot or three peaks per period in the time-series plot, i.e. the 3-cycle has three oscillations per period (which is the definition of a 3-cycle (Thompson & Stewart 1986)). The time period *T* of the 3-cycle is 3*τ*, where *τ* is the time period of the applied voltage (also shown in figure 12). In short, a *K*-cycle in the *U*-sequence has *K* oscillations per period and has a time period of (indicated by *K* points in the Poincarè map). Figures 13 and 14 show the phase plots and the Poincarè maps for the *K*=6 and *K*=4-cycle oscillations following the *K*=3-cycle window at 3.8085 V ac, as shown in table 1. For ac voltages beyond the final *K*=6-cycle window at 3.83173 V ac the MEMS mirror exhibits chaos and finally dynamic pull-in takes place at 3.85 V ac. Similar observation of *U*-sequence at resonant excitation is also made in the MEM fixed–fixed beam device.

### (b) UM-sequence

A modified form of the *U*-sequence is found to exist in electrostatic MEMS under superharmonic excitations and is presented here. The appearance of a periodic window with *K*-cycles in the normal *U*-sequence is replaced by the appearance of a periodic window with *KM*-cycles in the modified sequence termed the ‘*UM*-sequence’, for an exciting frequency of (*M*th superharmonic frequency). Hence, in the *UM*-sequence, the order given in equation (3.1) is replaced by the following order at the *M*th superharmonic frequency as the exciting frequency:(3.2)An *M*-cycle formed at the *M*th superharmonic frequency of excitation (before the -cycle PD sequence) is by definition characterized by *M* loops per period in the phase plot or *M* peaks per period in the time-series plot, as shown in figures 15 and 16 for *M*=1 and *M*=3, respectively. However, the essential difference between a *K*-cycle in the *U*-sequence and an *M*-cycle at the *M*th superharmonic excitation is that the time periods of these cycles are and *τ* respectively, where *τ* is the time period of the applied voltage in each case. As a result, the Poincarè map of a *K*-cycle in the *U*-sequence has *K* points whereas the Poincarè map of an *M*-cycle at the *M*th superharmonic excitation has just one point, as the reference frequency used in the Poincarè map is the frequency of the applied voltage.

The fixed–fixed beam MEM device described earlier exhibits dc symmetry breaking from around 69 V to pull-in at 73 V (De & Aluru 2005). An ac voltage at the *M*th superharmonic frequency of excitation is applied to the device after a steady state is reached due to the application of a dc bias of 71 V (in the dc symmetry breaking region). The resonant frequency is *f*_{0}=885 kHz at 71 V dc for the MEMS fixed–fixed beam. It has been shown in De & Aluru (2005) that for *M*=3 and 71 V dc bias, an ac voltage around 5.9 V gives rise to 3-cycle oscillations in the MEMS fixed–fixed beam. On further increasing the ac voltage, a sequence of PD bifurcations take place giving rise to -cycles and ultimately chaos at 6.01 V ac (De & Aluru 2005). In this paper, the ac voltage is gradually increased beyond 6.01 V, and the *UM*-sequence (the *U*3-sequence in this case) is observed. Table 2 shows the *U*3-sequence and the corresponding ac voltages at which the periodic states/cycles are observed at 71 V dc and *M*=3 for the fixed–fixed beam. Figure 17 shows the phase plot and the Poincarè map for the chaotic state after the sequence of PD bifurcations at 6.01 V ac. The time series and phase plot for the 15-cycle window (*K*=5, *M*=3) at 6.0137 V ac and 71 V dc are shown in figure 18. The vertical arrows above the time series are separated by one time period *T* of the response and *τ* is the time period of the applied voltage in figure 18. Within each time period *τ* of the applied voltage, a *M*=3-cycle oscillation is observed while the response repeats itself after *K*=5 voltage time periods, giving *T*=5*τ* and a 5×3=15-cycle response. Figures 19 and 20 show the time series and the phase plots for the 9-cycle window (*K*=3, *M*=3) at 6.0192 V ac and the 12-cycle window (*K*=4, *M*=3) at 6.0329 V ac, respectively, for the fixed–fixed beam at 71 V dc (refer to table 2).

In general, when a *K*-cycle in the *U*-sequence given by equation (3.1) is formed, the displacement repeats itself after every *K* voltage cycles (time period is *K* times the time period of the applied voltage). However, each voltage cycle consists of *M* loops/peaks/oscillations at the *M*th superharmonic frequency of excitation. As a result, the total number of oscillations per period of the response becomes *K*×*M* as given in equation (3.2) for an exciting frequency of (as described through figure 18). However, the time period of the *KM*-cycle remains , where *τ* is the time period of the applied voltage. The characteristic *M*=3 loops per period of the base orbit is found to be preserved in figures 18–20. The *UM*-sequence is also observed for other values of *M*. Figures 21 and 22 show the first -cycle state for an ac voltage of 6.3606 V and the first -cycle state for an ac voltage of 6.414 V in the MEMS fixed–fixed beam at 71 V dc for *M*=2. Similar observations of *UM*-sequences are also made in the MEMS torsion mirror.

### (c) Intermittent chaos

The periodic windows of the *U*-sequence are formed through intermittency (Thompson & Stewart 1986). Formation of periodic windows through intermittent chaos is observed for all the periodic windows in the *UM*-sequence. The density/proportion of the periodic states/cycles in the time-series data increases with the increase in the applied ac voltage with the complete disappearance of chaos when the cycle is formed. For the MEMS torsion mirror at 3.808 V ac, the times series is completely chaotic, whereas at 3.8085 V ac (table 1) it is completely periodic (*K*=3-cycle) as shown in figure 11 for a 141 V dc bias and *M*=1. Intermittent chaos is observed in between these two ac voltages. Figure 23 shows intermittent chaos in the time-series data at 3.8082 V ac, just before the formation of the *K*=3-cycle window. The regions in the time series (broken down into two parts, one from 100 to 450 ms and the other from 450 to 800 ms) marked by horizontal lines are periodic whereas the rest of the time series is chaotic. Figure 24 shows intermittent chaos in the time-series data at 3.8255 V ac, just before the formation of the *K*=4-cycle window in the MEMS torsion mirror at 3.8259 V ac for *M*=1 at 141 V dc (see table 1).

The periodic state with -cycles is observed at 6.0192 V ac and 71 V dc in the MEMS fixed–fixed beam as shown in table 2. At a slightly smaller voltage (e.g. 6.0188 V ac), the device is in a chaotic state. Intermittent chaos is observed during the formation of the 9-cycle from the chaotic state with the gradual increase of the applied ac voltage. Figure 25 shows intermittent chaos in the time-series data at 6.0190 V ac. There are temporal regions of periodic states (9-cycle) and regions which are chaotic in between, as shown in figure 25. Figure 26 shows similar observation of intermittent chaos just before the formation of the -cycle window for the MEMS fixed–fixed beam at 71 V dc and 6.0327 V ac. The -cycle state is formed at 6.0329 V ac (table 2).

### (d) PD bifurcations of the KM-cycles: 2^{n}KM-cycles

While the formation of the periodic windows takes place through intermittent chaos, their destruction (transition to chaos from the periodic windows) takes place through the PD route to chaos. The PD sequence of each periodic state in the *U*-sequence results in the formation of -cycles in a dynamical system (Thompson & Stewart 1986). For electrostatic MEM devices, due to the presence of the *UM*-sequence at the *M*th superharmonic frequency, -cycles are formed at . Figure 27 shows the first PD bifurcation of the 3-cycle (table 1) (*n*=0, *K*=3, *M*=1) into a 6-cycle (*n*=1, *K*=3, *M*=1) oscillation at 3.8107 V ac in the MEMS torsion mirror. The second PD bifurcation of the 3-cycle resulting in the formation of a 12-cycle (*n*=2, *K*=3, *M*=1) oscillation at 3.8111 V ac is shown in figure 28. Further PD bifurcations ultimately lead to chaos. Figure 29 shows the first PD bifurcation of the 9-cycle (table 2) (*n*=0, *K*=3, *M*=3) into an 18-cycle (*n*=1, *K*=3, *M*=3) oscillation at 6.0215 V ac for the MEMS fixed–fixed beam. Further PD bifurcations of the 18-cycle ultimately leads to chaos. The first PD bifurcation of the 12-cycle (table 2) (*n*=0, *K*=4, *M*=3) into a 24-cycle (*n*=1, *K*=4, *M*=3) oscillation at 6.0335 V ac for the MEMS fixed–fixed beam at 71 V dc is shown in figure 30. On further increasing the ac voltage, chaos is observed at 6.034 V ac.

The time period of a -cycle is since an *M*-cycle formed at the *M*th superharmonic frequency has the same time period as that of the exciting voltage. The time period of the response in comparison to the applied voltage is only affected by the *K*-cycle of the *U*-sequence and by the PD bifurcations. As the sampling frequency for the Poincarè maps in figures 29 and 30 is the exciting frequency, the *M*=3-cycle within each period does not appear in the Poincarè maps. As a result, the total number of cycles is obtained by multiplying the number of points in the Poincarè maps by *M*=3 leading to 18-cycles in figure 29 and 24-cycles in figure 30.

## 4. Analysis of the nonlinear dynamic properties

Nonlinearities in the coupled electromechanical–fluidic physical models can arise from all the three energy domains, namely, electrostatic, mechanical and fluidic. As a result, it would be difficult to identify the physical mechanism that gives rise to the nonlinear dynamic properties of electrostatic MEMS presented in §3 using the coupled electromechanical–fluidic physical models, although these models would predict accurate results when compared with the experimental data. In this regard, the MSD model presented in §2*b* (equation (2.14)) is very helpful as this model considers the electrostatic force as the only source of nonlinearity and neglects the nonlinearities in the mechanical and fluidic forces. The MSD model is applied to the fixed–fixed beam MEM device. The parameters in the MSD model, *m*, *c*, and *k* can be computed analytically for the fixed–fixed beam MEM device using the linearized Reynolds squeeze film damping equation (Huang *et al*. 2001) and the Euler's beam theory (Popov 1997). The parameter values are *m*=1.864×10^{−12} kg, *c*=2.832×10^{−6} Ns m^{−1} and *k*=105.625 N m^{−1}. Numerical simulations of the MSD model for the fixed–fixed beam exhibit the same dynamical properties (qualitatively) presented in §3 for the device, namely, *U*- and *UM*-sequences, intermittent chaos and the PD bifurcations of the *KM*-cycles giving rise to -cycles, but at different voltages. For example, figure 31 shows the time series and the phase plot of the first -cycle window in the fixed–fixed beam at 5.8639 V ac and 61 V dc obtained using the MSD model. The dc and the ac voltages at which the 18-cycle window is observed are significantly different from those given by the physical models (see table 2). The shift in voltages is because of the absence of mechanical and fluidic nonlinearities in the MSD model (see De & Aluru 2006 for details). Since the MSD model exhibits the same nonlinear dynamic properties as those from the physical models, it can be concluded that the nonlinear electrostatic force is the underlying physical force responsible for the dynamic properties presented in §3, as the mechanical and fluidic nonlinearities are absent in the MSD model. Besides, a harmonic balance analysis (Nayfeh & Mook 1979) of the MSD model was done in De & Aluru (2005) to show that the nonlinear nature of the electrostatic force gives rise to higher-order harmonics of the exciting frequency, resulting in the formation of *M*-cycles. The MSD model has quadratic and higher-order terms of *x* in (equation (2.14)). At small voltages, these higher-order terms can be neglected and the MSD model can be represented as a damped Mathieu equation (Turner *et al*. 1998),(4.1)In Turner *et al*. (1998), parametric resonance in MEMS has been shown at smaller voltages using equation (4.1) under low damping (small *c*). However, the nonlinear dynamic properties presented in this paper and in De & Aluru (2005, 2006) were observed at voltages where the quadratic and higher-order terms of *x* in are significant.

The concept of symmetry breaking in the phase plot of a dynamical system has been defined mathematically in Swift & Wiesenfeld (1984). Considering a periodically driven system (with a driving period of *T*) described by an *n*-dimensional system of ordinary differential equations,(4.2)the system is symmetric if *F* satisfies(4.3)This implies that a periodic solution of a dynamical system in general is symmetric if and only if the solution contains only odd harmonics (Wiesenfeld *et al*. 1984). Figure 32 shows the phase plot and the spectral density (computed using discrete Fourier transformation on the time-series data obtained from the numerical simulations) of the periodic solution of the MEMS fixed–fixed beam at 71 V dc, 1 V ac and *M*=1 obtained from the numerical simulation of the physical models. In this case, only the first-order harmonic of the ac excitation frequency is present in the solution and no ac symmetry breaking is observed. On the other hand, in figure 33, ac symmetry breaking in the phase plot, indicated by the presence of the second-order (even) harmonic of the excitation frequency in the solution at 71 V dc, 3 V ac and *M*=1 is observed. Similar observations can be made for dc symmetry breaking where the resonant frequency of the device is the driving frequency. Besides, the MSD model can be reduced to a pair of first-order ordinary differential equations as(4.4)Applying the condition for symmetry (equation (4.3)) to equation (4.4), it can be easily shown that for voltages where the higher-order terms of *x* in become significant, the MSD model is not symmetric (equation (4.3) does not hold).

The nonlinear dynamic properties of MEMS presented in this paper and in De & Aluru (2005, 2006) were observed in the MEMS fixed–fixed beam and the MEMS torsion mirror under normal atmospheric conditions (1 atm, 293 K). A large dc bias beyond dc symmetry breaking and a large ac bias beyond ac symmetry breaking were applied to observe these properties like *M*-cycles, -cycles, chaos and *U*/*UM*-sequences. However, under certain conditions like low atmospheric pressures, only one symmetry breaking may be sufficient to observe these properties. For example, figure 34 shows the formation of a *M*=3-cycle oscillation in the MEMS fixed–fixed beam using the MSD model, at a dc bias of 30 V (smaller than the dc symmetry breaking voltage of 60 V observed in the MSD model for the fixed–fixed beam). The *M*=3-cycle is formed at 20 V ac which is beyond the ac symmetry breaking voltage of around 10 V ac. Figure 35 shows the variation in the relative strength of the third harmonic (with respect to the first harmonic) of the MEMS beam at 30 V dc and 20 V ac with ambient pressure, using harmonic balance analysis on the MSD model. The relative strength of the third harmonic is found to be inversely proportional to the ambient pressure, which is also true for the other higher-order harmonics. As a result, *M*-cycles are formed more easily at smaller ambient pressures than under normal atmospheric conditions. Under normal atmospheric conditions, pull-in of the MEMS beam is found to occur before the higher-order harmonics can become stronger and form *M*-cycles at smaller dc voltages (below dc symmetry breaking).

The formation of the periodic windows (*U*-sequence) in a system through intermittent chaos and their destruction through the PD route to chaos has been well described by the theory of crisis in dynamical systems (Grebogi *et al*. 1982; Thompson & Stewart 1986). The disappearance of chaos (through intermittent chaos) leading to the formation of the periodic windows and the reappearance of chaos (through the PD route) leading to the destruction of the periodic windows is further characterized in this paper by computing the largest Lyapunov exponent *λ* of the MEMS fixed–fixed beam system using the software package TISEAN (Hegger *et al*. 1999). The time-series data from the numerical simulation of the fixed–fixed beam at 6.0190 V ac (intermittent chaos before the formation of the -cycle window (see figure 25)) is used in TISEAN and a positive value of is obtained, confirming that it is a chaotic state. The value of *λ* is found to be negative when the 9-cycle window is formed around 6.0192 V ac (table 2) and becomes positive again at 6.022 V ac after the PD bifurcation of the 9-cycle around 6.0215 V ac (figure 29). These observations are consistent with the crisis theory in dynamical systems.

## 5. Conclusion

In this paper, new nonlinear dynamic properties of electrostatic MEMS have been presented through the numerical simulation of detailed physical level models for the coupled electrical, mechanical and the fluidic energy domains. The presence of *U*-sequence in electrostatic MEMS is reported for the first time in this paper. A modified form of the *U*-sequence, the *UM*-sequence is shown to be present in electrostatic MEMS under superharmonic excitation. While the periodic windows in the *U*- and *UM*-sequences are formed through intermittent chaos, they are destroyed through the PD route to chaos. Period doubling bifurcations of the periods in the *UM*-sequence result in the formation of -cycles in electrostatic MEMS. The nonlinear nature of the electrostatic force is found to give rise to the *UM*-sequence and consequently the -cycles in MEMS under superharmonic excitations.

## Acknowledgments

This work is supported by the National Science Foundation under grant CCR-0121616 and grant ACI-0217986.

## Footnotes

- Received February 19, 2006.
- Accepted April 10, 2006.

- © 2006 The Royal Society