Simplified model and analysis of a pair of coupled thermo-optical MEMS oscillators

TL;DR

This paper develops a simplified third-order model for coupled thermo-optical MEMS oscillators, analyzing their dynamics and bifurcation structure with analytical and numerical methods.

Contribution

It introduces a new simplified model capturing key features of thermo-optical MEMS oscillators and analyzes the dynamics of coupled oscillators using analytical techniques.

Findings

Analytical results agree with numerical simulations.

Bifurcation diagram reveals the system's dynamic transitions.

Model effectively captures key oscillator behaviors.

Abstract

Motivated by the dynamics of micro-scale oscillators with thermo-optical feedback, a simplified third order model, capturing the key features of these oscillators is developed, where each oscillator consists of a displacement variable coupled to a temperature variable. Further, the dynamics of a pair of such oscillators coupled via a linear spring is analyzed. The analytical procedures used are the variational equation method and the two-variable expansion method. It is shown that the analytical results are in agreement with the results of numerical integration. The bifurcation structure of the system is revealed through a bifurcation diagram.

Tables

Table 1. Table 1 : Stability of the system to different values of α 𝛼 \alpha .

Value of $\alpha$	Stability
$>0.9$	only IP
$0.88$	OP becomes stable
$0.057$	IP loses stability in
	pitchfork to modified IP
$0.0468$	modified IP loses stability
	in Hopf to quasiperiodic
$0.0436$	two LCs coalesce in homoclinic
$0.042$	stable LC collides with separatrix
	and vanishes
	only OP

Equations

$$\ddot{z}+\frac{1}{Q}\dot{z}+\left(1+CT\right)z+\beta z^3=DT$$

$$\dot{T}=H{P}_L\left[\alpha +\gamma {\sin }^2\left(2\pi (z-z_0)\right)\right]-BT$$

$$\ddot{z}+z=T$$

$$\dot{T}+T=z^2-z{z}_0$$

$$\ddot{z}+z=T$$

$$\dot{T}+T=\epsilon z^2-{\epsilon }^{2}pz$$

$${z_{(0)}}^{\prime \prime }+z_{(0)}=T_{(0)}$$

$${T_{(0)}}^{\prime }+T_{(0)}=0$$

$${z_{(1)}}^{\prime \prime }+z_{(1)}=-2k_1{z_{(0)}}^{\prime \prime }+T_{(1)}$$

$${T_{(1)}}^{\prime }+T_1=z_{(0)}^2-p{z}_{(0)}-k_1{T_{(0)}}^{\prime }$$

$${z_{(2)}}^{\prime \prime }+z_{(2)}-2k_1{z_{(1)}}^{\prime \prime }+\left(-2k_2-k_1^2\right){z_{(0)}}^{\prime \prime }+T_{(2)}$$

$${T_{(2)}}^{\prime \prime }+T_{(2)}=(2z_{(2)}-p)z_{(1)}-k_1{T_{(2)}}^{\prime }-k_2{T_{(0)}}^{\prime }$$

$$z=A\cos \omega t+A^2\left(\frac{1}{2}-\frac{1}{15}\sin 2\omega t-\frac{1}{30}\cos 2\omega t\right)$$

$$T=A^2\left(\frac{1}{2}+\frac{1}{5}\sin 2\omega t+\frac{1}{10}\cos 2\omega t\right)$$

$$A=\sqrt{10p}/3$$

$$\omega =1-p/27$$

$${\ddot{z}}_1+z_1=T_1+\alpha \left(z_2-z_1\right)$$

$${\dot{T}}_1+T_1=z_1^2-z_{1}p$$

$${\ddot{z}}_2+z_2=T_2+\alpha \left(z_1-z_2\right)$$

$${\dot{T}}_2+T_2=z_2^2-z_{2}p$$

$$x=z_1+z_2$$

$$y=z_1-z_2$$

$$u=T_1+T_2$$

$$v=T_1-T_2$$

$$\ddot{x}+x=u$$

$$\dot{u}+u=\frac{x^2+y^2}{2}-px$$

$$\ddot{y}+\left(1+2\alpha \right)y=v$$

$$\dot{v}+v=xy-py$$

$$\delta \stackrel{...}{y}+\delta \ddot{y}+\left(1+2\alpha \right)\delta \dot{y}+\left(1+2\alpha +p-F(t)\right)\delta y=0$$

$$F(t)=2A\cos \omega t+2A^2\left(\frac{1}{2}-\frac{1}{15}\sin 2\omega t-\frac{1}{30}\cos 2\omega t\right).$$

$${\omega }^3{q}^{\prime \prime \prime }+{\omega }^2{q}^{\prime \prime }+\nu \omega q^{\prime }+\left(\nu +\mu \cos \tau +{\mu }^{2}p\right)q=0.$$

$$q_0^{\prime \prime \prime }+q_0^{\prime \prime }+q_0^{\prime }+q_0=0$$

$$q_0=m\cos \tau +n\sin \tau$$

$$q_1^{\prime \prime \prime }+q_1^{\prime \prime }+q_1^{\prime }+q_1=-\frac{m}{2}\left(1+\cos 2\tau \right)-\frac{n}{2}\sin 2\tau .$$

$$q_1=\frac{(m-2n)\cos 2\tau }{30}+\frac{(2m+n)\sin 2\tau }{30}-\frac{m}{2}$$

$$q_2^{\prime \prime \prime }+q_2^{\prime \prime }+q_2^{\prime }+q_2=-q_1\cos \tau -(p/9)q_0^{\prime \prime \prime }+(2/27)p{q}_0^{\prime \prime }$$

$$-({\nu }_2-(p/27))q_0^{\prime }-(p+{\nu }_2)q_0$$

$$\frac{2np}{27}+\frac{29mp}{27}+{\nu }_{2}n-\frac{n}{30}+{\nu }_{2}m-\frac{29m}{60}=0$$

$$\frac{29np}{27}-\frac{2mp}{27}+{\nu }_{2}n+\frac{n}{60}-{\nu }_{2}m+\frac{m}{60}=0$$

$$\frac{67600p^2+133920{\nu }_{2}p-29520p+116640{\nu }_2^2-31104{\nu }_2-405}{58320}=0$$

Full text

††journal: Nonlinear Dynamics

1 Introduction

The dynamics of resonant MEMS, or micro-electro mechanical systems, are of interest due to both the large number of current and potential applications Younis and for the new questions that arise in their design and analysis. Of interest here are MEMS systems consisting of a resonator suspended above a substrate and illuminated with a focused laser. The suspended resonator and substrate form a Fabry-Pérot interferometer group1 with the net result that the amount of laser heating of the resonator is modulated by the resonator/substrate gap. The resulting feedback of heating, temperature and thermal stress can result in limit cycle oscillations group2 . Prior work has shown that these oscillators can be entrained Pandey1 to inertial and optical signals Pandey2 .

Accounting for laser heating and thermo-mechanical coupling the oscillator can be described by the following equation of motion group4

[TABLE]

where $z$ is the displacement of the oscillator from its equilibrium position and $T$ is the temperature of the oscillator. In Eq. (1a), $z$ is the out-of-plane displacement of the oscillator, normalized by the wavelength of the laser light, $Q$ is the quality factor, $C$ is the thermal coefficient for linear stiffness, $\beta$ is the cubic stiffness, $D$ is the static position per unit change in temperature, $E$ is the inertial drive amplitude, and $\omega$ is the inertial drive frequency. In Eq. (1b), $T$ is the average temperature of the oscillator, $B$ and $H$ are thermal constants, $P_{L}$ is the continuous-wave laser power, $\alpha$ is the minimum absorption, $\gamma$ is the contrast in absorption, and $z_{0}$ is the equilibrium position of the oscillator with respect to the absorption curve. In Eq. (1a), the linear stiffness and the thermal driving force depend on the temperature and in Eq. (1b), the laser absorption is modulated by the displacement of the oscillator group3 .

The nonlinearities and presence of a large number of parameters make the equations intractable and thus we develop a model that has the essential features of the system but is amenable to analysis. The premise of our study is that understanding the dynamics of the simplified model will help in explaining phenomena associated with the more complex model, and the corresponding experiments. The dynamics of a single oscillator are first discussed followed by an analysis of the dynamics of a coupled pair of oscillators.

To the best of one’s knowledge this is one of the first studies where coupled third order oscillators have been considered. The archetypical system for synchronization is the phase only oscillator, where each oscillator has the structure $\dot{\theta}=\omega+f(\theta,t)$. Here $\theta$ denotes the phase of each oscillator while $f(\theta,t)$ represents various forms of coupling and/or nonlinear term. This kind of study was pioneered by Kuramoto Kuramoto and those results were later applied to biological situations by Mirollo and Strogatz Strogatz . The case of two coupled van der Pol oscillators was first studied by Storti and Rand Storti . They obtained the regions of in-phase and out-of-phase locking as well as drifting.

In recent years, interest has arisen in a multitude of other kinds of oscillators. Valente et. al. Valente have considered a system consisting of two masses coupled via a spring, which can rigidly impact a fixed stop. Hybrid periodic orbits are found as well as trajectories where there are infinitely many impacts with the stop occurring in a finite time. Fradkov and Andrevsky Fradkov have studied synchronization of two simple pendula coupled by a linear spring. Sliwa and Grygiel Sliwa have considered a pair of coupled Kerr oscillators. These have the equation of the structure $dz/dt=-j\omega z-j\epsilon a*a^{2}+Fe^{-j\Omega t}-\gamma z$ where $z$ is a complex variable, $j$ denotes the imaginary unit and all other symbols denote constants. This oscillator arises in nonlinear optics. Switching of the system between different semi-stable attractors as well as chaotic beats have been observed. Suchorsky and Rand Suchorsky have considered van der Pol oscillators coupled by fractional derivative. The authors have considered regions of locking and drifting and have demonstrated the reduction to known results in the limits where the fractional derivative is replaced by an integer. Finally, Chavez et. al. Chavez have considered a pair of forced Duffing oscillators coupled via soft impacts. A bifurcation diagram has been presented by these authors.

The previous paragraphs give some indication of the breadth and variety of oscillators and coupling mechanisms which have been considered in literature. However, all of these either consider phase-only or second order models. The introduction of a higher order system, in this work, presents new and interesting phenomena which will be discussed. The complexity of the system requires the use of a multitude of techniques. We will use a variational equation method, a regular perturbation theory on a Mathieu-like equation as well as a two-variable expansion on the system equations to generate a set of slow flow equations.

2 Simplified model of one oscillator

A key observation from the analysis of Eq. (1), is the presence of a limit cycle for a large enough value of $P_{L}$ group1 . Hence to emulate the full equations the reduced model must also support a limit cycle. The damping in Eq. (1a) is inessential as there is already a damping term in Eq. (1b), and that alone provides all the damping which we require for off-cycle motions to die out. The cubic nonlinearity term and temperature coefficient in the linear stiffness term in Eq. (1a) are not essential for the formation of a limit cycle. In Eq. (1b), $\sin^{2}(z-z_{0})$ can be Taylor expanded as

${{\left[\left(z-{{z}_{0}}\right)-{{\left(z-{{z}_{0}}\right)}^{3}}/6+...\right]}^{2}}$. Neglecting terms of order $O((z-z_{0})^{4})$ and higher, this term can be further simplified to $z^{2}-zz_{0}$ without losing any essential features. Thus Eq. (1) can be reduced to the following:

[TABLE]

A Lindstedt Poincaré type analysis Rand of this system is performed, assuming the amplitude of oscillation to be small, of order $\varepsilon$. Further, the assumption is made that the parameter $z_{0}$ is of the second order of smallness, and it is scaled as $z_{0}=\varepsilon^{2}p$ to get

[TABLE]

Using the perturbation theory, we expand $z=z_{(0)}+\varepsilon z_{(1)}+\varepsilon^{2}z_{(2)}$, $T=T_{(0)}+\varepsilon T_{(1)}+\varepsilon^{2}T_{(2)}$ and use the time dilation, $\tau=\omega t$, where $\omega=1+k_{1}\varepsilon+k_{2}\varepsilon^{2}$. At the lowest order we get,

[TABLE]

where prime denotes differentiation relative to the stretched time, $\tau$. Since Eq. (4b) has exponentially decaying solutions, they are not of interest in driving a limit cycle so we take $T_{(0)}=0$. With this substitution, Eq. (4a) has the standard oscillatory solutions. Since the system is autonomous, the phase is arbitrary and the oscillation can be treated as $z_{(0)}=A\cos\tau$. Balancing the coefficients of $\varepsilon$ we get,

[TABLE]

In Eq. (5b), since $T_{(0)}=0=T^{\prime}_{(0)}$, the last term drops out. Removal of resonance terms from Eq. (5a) gives $k_{1}=0$ but yields no information about $A$, which carries along as a parameter, giving $z_{(1)}$ and $T_{(1)}$ in terms of $A$. Balancing the coefficients of $\varepsilon^{2}$ in the perturbation analysis we get,

[TABLE]

Removing secular terms from Eq. (6a) yields the following solution:

[TABLE]

where

[TABLE]

Since the perturbation theory is developed in terms of the amplitude, one is not required to explicitly report $\varepsilon$ in these expressions. It is noted that $A$ and $\omega$ are the amplitude and frequency of the fundamental motion $z_{(0)}$, respectively. Thus Eq. (2) exhibits limit cycle oscillations for all positive values of $p$. The amplitude increases and frequency decreases with increasing $p$.

3 Coupled pair of oscillators and variational equation method

We now turn to the analysis of two coupled MEMS oscillators, here assumed to be identical. The case where detuning is present will be considered in a later study. Linear coupling is assumed, which could be achieved either via a spring or via electrostatic coupling. The equations of the coupled oscillators are

[TABLE]

where $\alpha$ is the coupling stiffness. The model consists of two identical coupled oscillators, and hence in-phase (IP) and out of phase (OP) modes are expected. Here , the IP mode is defined as $z_{1}(t)=z_{2}(t)$ and $T_{1}(t)=T_{2}(t)$ whereas the OP mode is defined as $z_{1}(t)=-z_{2}(t)$ and $T_{1}(t)=-T_{2}(t)$. Indeed, linearization about the origin yields all six eigenvectors to have the structure corresponding to IP and OP modes. Two of the modes are associated with negative real eigenvalues for all $p$ and $\alpha$. The four remaining eigenvalues are two complex conjugate pairs, which are purely imaginary for $p=0$ and have positive real parts when $p>0$. One pair of these eigenvalues attaches to a pair of eigenvectors which has the IP mode shape while the other pair attaches to an eigenvector pair which has the OP mode shape. It is noted that in the full nonlinear system, the change in sign of the eigenvalues as $p$ is increased through zero corresponds to the Hopf bifurcation which gives rise to limit cycle oscillations.

Eq. (9) is transformed in the following way :

[TABLE]

Substituting Eq. (10) into Eq. (9), one gets,

[TABLE]

The presence of IP mode is clearly demonstrated by these equations - in this mode $y=v=0$, and $x$ and $y$ satisfy Eq. (2). An exact OP mode however does not exist; setting $x=u=0$ leads to an apparent contradiction. The problem of persistence of linearized modes in nonlinear systems has been considered by Hennig Hennig .

3.1 Stability of the IP mode

The stability of the IP mode is now analyzed by constructing the variational equation, as in Stoker Stoker . Small perturbations $\delta x$, $\delta y$, $\delta u$ and $\delta v$ are introduced on top of the steady state solutions of $x$,$y$,$u$ and $v$, where the steady state solutions are characterized by $x=2z$ and $u=2T$, where $z$ and $T$ are given by Eq. (7), and $y$ and $v$ are zero. This variational equation for the $\delta y$-$\delta v$ system can be combined into a single third order equation as

[TABLE]

where $F(t)$ is the limit cycle motion exhibited by $x(t)$ whose stability is being sought. Since $x=z_{1}+z_{2}$ which equals $2z_{1}$ for the IP mode, and since the limit cycle for a single oscillator was found in Eq. (7), one may write

[TABLE]

A perturbative approach will now be used to determine the stability of the IP mode. For this approach it is sufficient to retain only the first term in the right hand side of Eq. (13). From this point onwards, $\delta y$ is written as $q$ for notational convenience. To carry out a perturbation procedure the variables are rescaled as follows :

[TABLE]

where $\tau=\omega t$, the prime denotes differentiation with respect to the variable $\tau$, $\mu$ is the amplitude of the parametric excitation and $\nu=1+2\alpha$. We use the expansion, $\nu=1+\mu\nu_{1}+\mu^{2}\nu_{2}$, let $\omega=1-\mu^{2}p/27$ and then expand $q=q_{0}+\mu q_{1}+\mu^{2}q_{2}$. Before proceeding further, an outline of the philosophy behind the perturbation is given. To obtain the transition curves of Eq. (12), one seeks the motion that is periodic on the curve itself, with a period equal to that of $F(t)$. The choice of period is motivated by the fact that in the absence of the perturbations, Eq. (14) possesses oscillatory solutions with natural frequency 1; since the frequency of the parametric excitation is also perturbatively close to 1, we are close to a 1:1 resonance of the system. Substituting the perturbation expansions into Eq. (14), at zero order one has

[TABLE]

which has the fundamental solutions $\cos\tau$, $\sin\tau$ and $\operatorname{e}^{-\tau}$. The exponentially decaying solution is not of interest so we let

[TABLE]

where $m$ and $n$ are arbitrary. Equating terms of order $\mu$ in the expanded form of Eq. (14), and substituting Eq. (16) into the result, we get

[TABLE]

The resonant terms can be removed from this equation if one sets $\nu_{1}=0$. Using the method of undetermined coefficients, one finds

[TABLE]

At order $\mu^{2}$ we get

[TABLE]

where $q_{0}$ and $q_{1}$ have been determined above. This contains resonance terms on the right hand side; removal of which requires

[TABLE]

Equations (20) are a pair of simultaneous linear homogeneous equations; a non-trivial solution exists if and only if the determinant of the matrix vanishes. This condition leads to

[TABLE]

For $p=0.1$ the value of $\nu_{2}$ obtained from Eq. (21) is 0.2455. Now, $\nu$ is $1+\mu^{2}\nu_{2}$ which also equals $1+2\alpha$. Since $\mu$ is twice the amplitude of the single oscillator limit cycle, determined earlier to be 0.333 using Eq. (8a), the transitional value of $\alpha$, above which the IP mode is stable and below unstable, is found to be 0.0545. The theoretical predictions are compared with the numerical simulations. At $\alpha=0.1$, with initial conditions $\{z_{1}(0),z_{2}(0),\dot{z}_{1}(0),\dot{z}_{2}(0),T_{1}(0),T_{2}(0)\}=\{0.1,0.09,0,0,0,0\}$, close to the IP mode, the numerical solution shows an IP mode as shown in Fig. 1. For the same initial conditions, as $\alpha$ is decreased to 0.055, the straight line of the IP mode branches out into an ellipse as show in Fig 2. The critical value of $\alpha$ is very close to the value predicted by the variational method. For initial conditions, $\{z_{1}(0),z_{2}(0),\dot{z}_{1}(0),\dot{z}_{2}(0),T_{1}(0),T_{2}(0)\}=\{0.1,-0.09,0,0,0,0\}$, close to the OP mode, the numerical solution resembles an OP motion as shown in Fig. 3. Although the IP mode is an exact straight line, the OP mode is not.

An intriguing feature is that the variational analysis on the $x$,$u$ equations, considering perturbations $\delta x$ and $\delta u$ off their steady state values, leads to the same equation as Eq. (12) but with $\alpha=0$. This is supposed to be unstable. However, we recall that a variational analysis gives Lyapunov stability whereas the limit cycles possesses only orbital stabilityStoker . In the same way, a variational analysis of the limit cycle of van der Pol oscillator also yields a ‘spurious’ instability. The source of the instability is phase shear. As the trajectories approach the limit cycle, their period change, reaching a limiting value at the cycle itself. If the motion of a point exactly on the limit cycle is compared with that of a point, $\varepsilon$ distance away from it, then it will be seen that the the separation between the two points initially increases from $\varepsilon$ to a finite value due to the unequal periods of their trajectories.

3.2 Stability of the OP mode

It is now desired to determine the stability of the OP-like mode. It was mentioned earlier that the strict OP mode $z_{1}=-z_{2}$ and $T_{1}=-T_{2}$ does not exist for Eqs. (9). In order to investigate the stability of the motion similar to the OP mode (cf. Fig. 3), one first needs to obtain an expression for the motion, which one will do using LINDSTEDT’s method. The variable scaling is the same as that of a single oscillator Eq. (3), and the solutions $z_{1}=A\cos(\omega t)$ and $z_{2}=-A\cos(\omega t)$ are imposed at the $\varepsilon^{0}$ level. Computer algebra is used for the calculation, owing to intractability of the mathematical manipulations. As in Section I, the starting step is to introduce the parameter $\varepsilon$, in a manner consistent with the scalings for a single oscillator Eq. (3). Time is rescaled as $\tau=(1+k_{1}\varepsilon+k_{2}\varepsilon^{2})t$ and the variables are expanded as $z_{1}=z_{10}+\varepsilon z_{11}+\varepsilon^{2}z_{12}$, and so on. The zero order equations are

[TABLE]

It is found that the $z$’s are oscillatory and the $T$’s are damped. Neglecting the exponentially decaying terms, the frequency of oscillation of Eqs. (22a) and (22c) is obtained as $(1+2\alpha)^{1/2}$. Now, we impose externally that $z_{10}=A\cos\omega\tau$ and $z_{20}=-A\cos\omega\tau$. The order $\varepsilon$ equations are

[TABLE]

Plugging the zero order solutions into these equations leads to a coupled inhomogeneous system of ODEs for the 1-level variables. To uncouple this system the transformation to sum and difference coordinates is performed with $z_{11}=(u+v)/2$ and $z_{21}=(u-v)/2$. Removal of the resonance term gives $k_{1}=0$. Solving the system in the new variables and inverting gives $z_{11}$ and $z_{21}$ as trigonometric functions of $t$, with the amplitude $A$ still as a parameter.

Determination of amplitude is achieved from the next level. The equations at this level are

[TABLE]

Substitution of the already determined 0- and 1- level quantities into Eq. (24), followed by yet another sum and difference transformation, leads to a situation where the resonance terms can be removed. Computer algebraic manipulations yield the following :

[TABLE]

where

[TABLE]

The trajectory predicted by Eqs. (25) and(26) is plotted in in Fig. 4 as a phase portrait, which has the desired figure-of-eight shape but is wider than the actual, as in Fig. 3. This solution will now be substituted into Eqs. (9) to construct the linear variational equation for perturbations $\delta z_{1}$, $\delta z_{2}$, $\delta T_{1}$ and $\delta T_{2}$ added onto the “steady state” solutions $z_{1}^{\ast}(t)$, $z_{2}^{\ast}(t)$, $T_{1}^{\ast}(t)$ and $T_{2}^{\ast}(t)$. The linear variational equation turns out to be

[TABLE]

Unfortunately, these equations do not uncouple upon introducing the substitution Eq. (10). This is because a strict OP mode is a not a solution of Eqs. (9). Hence, numerical simulation is resorted to for the solution of the above system. It yields that the perturbations remain bounded for $\alpha<0.88$ while they grow exponentially for $\alpha>0.88$. Numerical integration of the full system Eqs. (9) yields a transition from stable to unstable at approximately $\alpha=0.82$ which is in good agreement with the above prediction.

Fig. 5 summarizes the foregoing sections by showing the stability of the IP and OP modes as a function of $\alpha$.

4 Two variable expansion method

The method of two variable expansionNayfehMook , also known as the method of multiple scales Nayfeh , will now be used to study the coupled system. In this process, time is split into two variables, a regular time $\xi=t$ and a slow time $\eta=\varepsilon t$. The time derivative can be written as $\mathrm{d}/\mathrm{d}t=\partial/\partial\xi+\varepsilon\partial/\partial\eta$ and the second and higher derivatives may be calculated similarly. The small parameter $\varepsilon$ is introduced into Eq. (9) in the following manner :

[TABLE]

This scaling is motivated as follows : it is reasoned that the coupling, $\alpha$ is very weak ($O(\varepsilon^{2})$) but on account of the coupling, larger values of the static position $z_{0}$ can now be permitted. The following ansatz is attempted :

[TABLE]

and additional dynamical variables are not taken for $T_{1}$ and $T_{2}$ since they are exponentially decaying at largest order. Substituting Eqs. (29) into Eqs. (28) and removing resonant terms leads to the trivial slow flow $A^{\prime}=B^{\prime}=C^{\prime}=D^{\prime}=0$ where prime denotes $\mathrm{d}/\mathrm{d}\eta$. This indicates that the method has to be carried out to one further order. Defining a super-slow time $\zeta=\varepsilon^{2}t$, $A,B,C,D$ are now taken to be functions of $\zeta$. Removal of resonance terms from the resulting equations now gives the following slow flow

[TABLE]

[TABLE]

For increased convenience one can express the above system in terms of polar coordinates i.e. $A=r_{1}\cos\theta_{1}$, $B=r_{1}\sin\theta_{1}$, $C=r_{2}\cos\theta_{2}$ and $D=r_{2}\sin\theta_{2}$. If one defines the phase difference $\varphi=\theta_{2}-\theta_{1}$ then the following system for $r_{1}$, $r_{2}$ and $\varphi$ is obtained from Eq. (30):

[TABLE]

Note that Eqs. (31) exhibit a symmetry: they are invariant under the transformation

[TABLE]

By inspection it can be seen that there are fixed points when $r_{1}=r_{2}$ and $\sin\varphi=0$ i.e. $\varphi=0$ or $\varphi=\pi$. The former corresponds to the IP mode while the latter is the OP mode. A calculation for $r$ (either 1 or 2) at the fixed point yields $\frac{\sqrt{10p}}{3}$, which is same as that obtained from the Lindstedt analysis of Section 2. To obtain the stability of the IP mode, the Jacobian is constructed for Eq. (31),

[TABLE]

For $p=0.1$, the critical value of $\alpha$ is 31/540 = 0.057, where the IP mode is stable above this critical value and unstable below it. This is in good agreement with the variational equation and the results of numerical integration. Construction of the Jacobian for the OP mode shows that it remains stable at all values of $\alpha$. Thus, the two-variable method has yielded correctly the stability transition of the IP mode. However it cannot yield the transition for the OP mode, at this level of the perturbation theory.

Recall that we have shown that there is a range of $\alpha$ values for which both the IP and OP modes are stable, Fig. 5. It is noted that these modes are limit cycles in the actual dynamical system but fixed points in the slow flow. These slow flow equilibria are separated by an unstable slow flow limit cycle which we shall refer to as a separatrix. Although it is unstable, one may nevertheless see what the separatrix looks like by choosing initial conditions to (approximately) lie on the basin boundary between the two equilibria. Moving from the three-dimensional slow flow space to the 6 dimensional space of Eqs. (9), with initial conditions $\{z_{1}(0),z_{2}(0),\dot{z}_{1}(0),\dot{z}_{2}(0),T_{1}(0),T_{2}(0)\}=\{0.1,0.00209816,$

$0,0,0,0\}$ which are very close to the threshold, the separatrix appears as a quasi periodic motion, seen in Figs. 6 and 7. Such a situation was already encountered in Storti and Rand Storti .

AUTO AUTOcite is an analytic continuation software package which we will use to plot the bifurcation diagram of the slow flow Eqs.(31). The result is shown in Fig. 8 where slow flow fixed points are displayed in the $\varphi$ vs. $\alpha$ plane. It can be seen that the IP mode is stable up to $\alpha=0.057$; thereafter it cedes stability to a pair of slow flow fixed points which are born in a pitchfork bifurcation.

Since the slow flow system Eq. (31) possesses the symmetry Eq. (32), all dynamical quantities either have this symmetry or are part of a pair of reflected twins. AUTO shows that these slow flow fixed points lose stability in a Hopf bifurcation at approximately $\alpha=0.0468$; the resulting symmetric slow flow limit cycles are stable. Numerical integration shows that these limit cycles increase in amplitude and deform in shape as $\alpha$ is decreased further, up to a further bifurcation which occurs when $\alpha$ decreases through approximately 0.0436, though this is not shown in Fig. 8. In this case there is a homoclinic bifurcation in which the two asymmetric slow flow limit cycles join to become a single slow flow limit cycle which exhibits the symmetry of Eq. (32). In Fig. 9, one sees the two limit cycles obtained by numerical integration of the slow flow with parameter value $\alpha=0.0437$ and initial conditions $r_{1}=r_{2}=1/3$ and $\varphi=+0.45$ and $\varphi=-0.45$ respectively. In Fig. 10, one can see the single symmetric limit cycle obtained at $\alpha=0.0435$ and starting from either of the two initial conditions used in obtaining the mirror symmetric limit cycles, earlier.

Another bifurcation occurs when $\alpha$ decreases through approximately 0.0415, in which the unstable separatrix limit cycle (Fig. 6) merges with the symmetric slow flow stable limit cycle which was created in the homoclinic bifurcation. For values of $\alpha$ less than approximately 0.0415, the OP mode is the only stable motion.

To visualize the bifurcation at $\alpha=0.0415$, the results of simulation are shown with IC picked on the crack between attainment of limit cycle and transition to OP mode. Fig. 11 presents the results of numerical integration of the slow flow with parameter value $\alpha=0.042$ and ICs $r_{1}=r_{2}=1/3$ and $\varphi=0.4403$ while in Fig. 12 we change the IC $\varphi$ to 0.4404. At this $\alpha$ value, there are two stable objects; the symmetric limit cycle and the fixed point corresponding to the OP mode. In these simulations it can be seen that the system remains on the separatrix for a short time before it progresses to its asymptotic state. In Fig. 13, the value of $\alpha$ is changed to 0.041. Starting at the initial conditions corresponding to Fig. 11, the system goes to the OP mode.

5 Conclusions

Thus, in this article a study of a simplified model of coupled MEMS oscillators is performed. The most interesting feature of the system is that the IP and the OP modes are both stable over a considerable range of parameter values. For sufficiently low $\alpha$ $(<0.04)$ only the OP mode is stable. At intermediate $\alpha$, both the IP and the OP modes coexist stably. At high $\alpha$ $(>0.88)$ only the IP mode is stable. This coexistent stability of the two modes is reminiscent of what was found in Reference Storti for two coupled van der Pol oscillators. Table 1 shows the bifurcations that occur in this system as the value of $\alpha$ is varied.

As regards the various methods adopted, it is noted that the two variable method yielded the most accurate results in the regions where it was effective. However, it is ineffective in predicting the stability of OP mode. It is also an interesting fact that this is one system which is driven by a third order equation, and its stability analysis led us to a third order Mathieu-like equation. It is likely that these studies can be extended to incorporate other third order parametrically excited systems and MEMS systems with detuned oscillators in the future.

6 Acknowledgements

This material is based upon work supported by the National Science Foundation under grant number CMMI-1634664. The authors wish to thank Professor John Guckenheimer for advising them on the bifurcations involved in this paper.

7 Compliance with ethical standards

There is no conflict of interest. Funding has been received from NSF as acknowledged above. The entire research presented here is the authors’ own. No part of this Article has been reproduced from other Articles or is under consideration elsewhere. Research on human and animal subjects was not necessary for this project. All authors consent to submission of this Article in its present form.