A Generic Rotating-Frame-Based Approach to Chaos Generation in Nonlinear M/NEMS Resonators

Chaotic dynamics have received interest owing to their extraordinary ability to generate complex behaviors, such as synchronization patterns, even in simple and fixed arrangements of coupled nodes. Countless applications have been discussed, spanning control systems, telecommunications and neuroscience Hilborn et al. (2000); Ott (2002); Buscarino et al. (2014); Boccaletti et al. (2018). Recently, the field has witnessed a resurgence of interest due to the possibility of building large-scale hardware reservoirs from coupled nonlinear oscillators. To meet the requirements for practical application, high-integration and low-power implementation are necessary Minati (2018); Tanaka et al. (2019).
 Micro- and nano-electromechanical systems (M/NEMS) provide experimental platforms for investigating and generating such dynamics, as they are easily amenable to very large-scale integration, low-power operation, and they inherently exhibit rich nonlinear behavior Lifshitz and Cross (2008). Indeed, chaos generation is well-reported in the M/NEMS scientific literature Güttinger et al. (2017); Karabalin et al. (2009); Kenig et al. (2011); Bienstman et al. (1998); Ashhab et al. (1999); Park et al. (2008); De and Aluru (2006); Wang et al. (1998); DeMartini et al. (2007); Towfighian et al. (2011). However, in the existing studies, chaos generation was obtained in a manner that is system/device specific and neither provides a minimalist approach nor a general one to generating chaos in such devices. What is meant here by general is a procedure that is device-independent, i.e., an approach that would reliably generate chaos without detailed knowledge of the material, actuation, dimensions, or design of M/NEMS devices. The term minimalist indicates the smallest number of phase-space dimensions ($n$) necessary for the onset of chaos (i.e., $n=3$).
 Chaos generation reported in the literature for M/NEMS devices resorts to non-general properties, e.g. non-smooth nonlinearity Bienstman et al. (1998), a high number of phase-space dimensions ($n>3$) obtained via inter-device or inter-modal coupling Karabalin et al. (2009); Güttinger et al. (2017); Kenig et al. (2011), or extreme nonlinearity by operating near the electrostatic pull-in Park et al. (2008); De and Aluru (2006). However, a common approach is based on the creation of a static double-well potential either through electrostatic forces or by using buckled structures Wang et al. (1998); DeMartini et al. (2007); Towfighian et al. (2011); Emam and Nayfeh (2004). None of these approaches hinges around a common denominator property: they are, therefore, not transferable to nonlinear M/NEMS resonators in general. In fact, they commonly require a large actuation voltage, negating one of the main appeals of M/NEMS devices.
 An interesting case, is that of the static double-well potential, which is described by a Duffing equation (a system possessing a cubic nonlinearity)where the linear component is negative and the cubic component is positive. Period-doubling bifurcations and chaos in such systems have been studied Sharma et al. (2012) and been subject to experimental investigations Moon (1980); Holmes (1979); Moon and Holmes (1979). While such systems can be reproduced in M/NEMS devices Wang et al. (1998); DeMartini et al. (2007); Towfighian et al. (2011); Emam and Nayfeh (2004), nearly all M/NEMS devices inherently exhibit a different type of cubic nonlinearity Kaajakari et al. (2004), equally captured by a variant of the Duffing equation, whereby the linear component is positive and the cubic component can be either positive, or negative. Such M/NEMS resonators can exhibit dynamic bistability Cleland (2013).
 In contrast to the static double-well systems, the approach presented here relies on the dynamically-generated double-well pseudo-potential that is created when a generic nonlinear Duffing resonator is driven into the bistable regime Cleland (2013). Since bistability is accessed when the resonator’s vibration amplitude is on the order of some scaling parameter (e.g. thickness Kaajakari et al. (2004)), such pseudo-potentials can be generated and manipulated by changing the drive conditions without requiring device or setup redesign.
 The main purpose of this work is to demonstrate chaos generation in a perturbed “dynamic double-well” in a manner that parallels chaos generation in driven “static double-well” systems, albeit with significantly lower drive amplitudes. Thus, two drive tones are applied to the system, whereby the first creates the “dynamic double-well” and the second perturbs it.
 Since the displacement of the M/NEMS resonator studied herein is moderate, i.e. on the order of the structural thickness, a perturbation-based approach to analyse the dynamics is justified Cleland (2013); Greywall et al. (1994); Kozinsky et al. (2007). To represent the underlying dynamics we employ the rotating frame approximation (RFA), whereby a “slow flow” (a time-varying envelope) dynamics is overlaid on top of an otherwise purely sinusoidal response, and the timescales ($\mathcal{\tau}$) associated with this “slow flow” are on the order of the resonator line-width ($\mathcal{\gamma}$), i.e., $\mathnormal{\tau\approx}~{}\mathcal{O}\mathnormal{(1/\gamma)}$ Greywall et al. (1994); Cleland (2013).
 When a Duffing resonator is driven near resonance, its steady-state response as seen in the RFA corresponds to a fixed-point in the phase space; in case the resonator is driven in the bistable regime, the response shows two distinct stable fixed points and a saddle-node point Kozinsky et al. (2007). This latter configuration implies a homoclinic connection (i.e., a trajectory that starts and ends in the saddle-node and orbits one of the stable fixed points) may exist in the RFA phase-space of a Duffing resonator. Thus, just as a “static double-well” potential provides a homoclinic connection in the rest frame, so does the dynamic Duffing bistability provide a homoclinic connection in the rotating-frame, and it is the perturbation of such homoclinics that is responsible for the generation of chaos Mel’nikov (1963); Holmes (1979); Sharma et al. (2012).
 This argument can be demonstrated by considering the Duffing equation for a two tone-driven M/NEMS resonator, given as Cleland (2013)

where $x$ is the displacement, and $\gamma$, $\omega_{0}$, $\alpha$ are respectively the damping, natural frequency, and Duffing nonlinearity of the resonator. ${F_{1}}$, ${F_{2}}$, $\omega_{1}$, and $\omega_{2}$ are the amplitudes and frequencies of two-externally applied driving forces, and $\eta$ is the transduction coefficient.We introduce dimensionless constants as ${\bar{t}=t\times\omega_{0}}$, ${\bar{\gamma}=\gamma/\omega_{0}}$, ${\bar{\alpha}=\alpha/\omega_{0}^{2}}$, ${\bar{{F_{1}}}={\eta F_{1}}/\omega_{0}^{2}}$, and ${\bar{{F_{2}}}={\eta F_{2}}/\omega_{0}^{2}}$. Hereon, all equations are written using this form, however, the bars are dropped for convenience.
 The application of the RFA, in which the modal displacement takes the form ${x(t)=A(t)\cos(\omega_{1}t+\phi(t))}$, where ${A(t)}$ and ${\phi(t)}$ are slowly varying amplitude and phase envelopes, gives the following rotating-frame system

where ${X=A\cos(\phi)}$ and ${Y=A\sin(\phi)}$ are the rotating-frame quadratures, ${A^{2}=X^{2}+Y^{2}}$, and ${\delta=(\omega_{1}-\omega_{0}})/\omega_{0}$ (details are provided in the supplementary materials).
 Since Eqs. (1)-(2) are generically applicable to Duffing-type resonators, the results below can, in principle, be implemented in various physical realizations of nonlinear resonators, such as optical Hansson and Wabnitz (2014) and superconducting resonators Manucharyan et al. (2007), without loss of generality.
 Initially, consider the conventional case with only one applied force, i.e., ${{F_{2}}=0}$, whereby the system in Eq. (2) is reduced to the first two equations only. The fixed points of the system, obtained by setting the time-derivative in Eq. (2) to zero, exhibit a bistable response in a region of the dimensionless parameter space shown in Fig. 1(a) for a lossless and a low-loss ($\gamma=10^{-3}$) driven Duffing resonators. To visualize the phase space and associated homoclinic orbit, we select a constant-force cut through the parameter space, Fig. 1(b), and then a constant-detuning cut where bisatbility exists, Fig. 1(c). It is convenient to plot the RFA Hamiltonian ($\mathcal{H}_{\mathrm{RFA}}$) Dykman et al. (2019); Huber et al. (2019), shown in Fig. 1(c), along with the fixed points. For the case $\gamma=0$, trajectories on the $\mathcal{H}_{\mathrm{RFA}}$ surface follow closed orbits around the fixed points, the so-called libration orbits Houri et al. (2020a). A homoclinic orbit is then the limit case in which the libration orbit intersects the saddle-node point, as shown in Fig. 1(c).

Note that for $F_{2}=0$ the system of Eq. (2) is reduced to an autonomous two-dimensional system ($n=2$), which can not generate chaos as it lacks the necessary dimensionality. However, an additional time dependence introduced by making ${{F_{2}}\neq 0}$ increases the RFA dimensions from $n=2$ to $n=3$, thus in principle meeting the condition for chaos generation. Thus, by setting ${{F_{2}}\neq 0}$ chaos could be generated if the homoclinic orbit is sufficiently perturbed.
 The above arguments are corroborated by a combination of numerical simulations and measurements. The experiments are performed using a micro-beam GaAs piezoelectric MEMS resonator driven into the nonlinear Duffing regime, details regarding fabrication and basic properties of such resonators can be found in Ref. Yamaguchi (2017). The resonator is placed in a vacuum chamber, and its motion is measured using a laser Doppler vibrometer, whose output is fed simultaneously to a lock-in amplifier and a vector signal analyzer (for experimental details, see supplementary information). A higher harmonic mode is selected for these experiments to avoid possible inter-modal interactions Houri et al. (2020b, 2019).
 The application of a single-tone sweep produces the frequency response shown in Fig. 2(a) for the linear (${100~{}\mathrm{mV_{PP}}}$, black trace) and nonlinear Duffing regimes (${3~{}\mathrm{V_{PP}}}$, red/blue traces), which, upon fitting, give the following values for the resonator parameters ${\omega_{0}/2\pi=1.56}$ MHz, $Q~{}=~{}{1000}$, and ${\alpha=1.67\times 10^{15}\mathrm{Hz/V^{2}}}$.
As a first demonstration of the period-doubling route to chaos, a two-tone excitation is applied to the resonator with one fixed tone in the bistability region ($F_{1}~{}=~{}3~{}\mathrm{V_{PP}}$, $\omega_{1}/2\pi~{}=~{}1566.5$ kHz) and one swept tone. For large detunings between the two tones, the rotating-frame response corresponds to a low-amplitude libration oscillation having a frequency ${\Omega}$. As the tone is swept, the oscillations exhibit a quick succession of period-doubling bifurcations leading to chaos, as shown in Fig. 2(b) for both the high- and low- amplitude branches. The corresponding phase-space plots for period 1 (P1), period 2 (P2), and chaotic attractors are shown in Figs. 2(c)-(e), respectively. Note that, as will be demonstrated later, Fig. 2(b) indicates that the low-amplitude branch only generates chaos for negative detuning, i.e. $\Omega<0$, while the high-amplitude branch only generates chaos for $\Omega>0$.

Making ${{F_{2}}\neq 0}$ is clearly a necessary but not sufficient condition for chaos generation, and the question of whether and where in the 4-dimensional parameter space (${\delta}$, ${\Omega}$, ${F_{1}}$, and ${F_{2}}$) chaos exists remains to be addressed. While bounds have been set on the values of $F_{1}$ and $\delta$ such that bistability exists for $F_{2}=0$, similar bounds for $F_{2}$ and $\Omega$ have yet to be determined. This task is usually performed by applying Melnikov’s method, which sets strict conditions for the period-doubling bifurcation route to chaos to take place Moon and Holmes (1979); Sharma et al. (2012); Holmes (1979); Moon (1980, 1984); Moon and Li (1985).
 Ideally, the application of Melnikov’s method constrains chaos generation in parameter-space with an analytical bound. Unfortunately, straightforward application of Melnikov’s integral to the archetypal nonlinear resonator captured by Eq. (1) results in a relation that is not easily amenable to analytical solution; therefore, a heuristic approach analogous to the one employed in Ref. Moon (1980) is used. This approach considers that the application of a second forcing term creates the libration orbits, and the amplitude at which these libration orbits are large enough to undergo inter-well jumps is considered to be a lower bound for the onset of period-doubling bifurcation. By linearizing the libration orbit around the fixed point, an approximate closed-form bound for the period-doubling route to chaos and the associated minimum of that bound can then be expressed as (see supplementary information for detailed derivation):

where $\mathrm{F_{{2min}}}$ is the minimum required ${F_{2}}$ necessary to induce period doubling, $\mathrm{\omega_{L}}$ is the libration frequency, and $A_{1}$, $A_{2}$, and $A_{3}$, are the amplitudes of the three steady-state fixed points, with $A_{2}$ being the unstable one.
 Eq. (3) sheds light on the experimental results in Fig. 2: by realising that $\mathrm{\omega_{L}}<0$ for the low-amplitude branch and $\mathrm{\omega_{L}}>0$ for the high-amplitude branch, it is easy to understand that period-doubling in Fig. 2(b) takes place mostly for $\Omega\approx\ \mathrm{\omega_{L}}$.
 A two-dimensional sweep of both $\Omega$ and ${F_{2}}$ provides further grounds for comparison between the numerical, analytical, and experimental results. Such a sweep is shown in Fig. 3 for both solution amplitude branches and the same drive conditions as in Fig. 2, i.e. $\omega_{1}/2\pi~{}=~{}1566.5$ kHz and $\mathrm{F_{1}}~{}=~{}3~{}\mathrm{V_{PP}}$. Figure 3 plots the lag of the auto-correlation maximum for the experimentally and numerically obtained time-domain signals, where an auto-correlation lag of 1 indicates P1 orbits, a lag of 2 indicates P2 orbits, and so on. Experimental and numerical data agree well in predicting the region corresponding to P2, the higher order bifurcations, and chaos. It is also interesting to note that both types of time-domain data show that for some parameter-space values the initial condition branch is unstable, and inter-branch jumps occur. Chaos is verified for a selection of experimental traces, where a correlation dimension Hegger et al. (1999) of $D_{2}=2.2\pm 0.1$ and $D_{2}=2.3\pm 0.1$ is estimated for the high- and low-amplitudes branches respectively. When a similar analysis is undertaken for simulated time-series, the corresponding largest Lyapunov exponent calculated directly from the differential equations is $\lambda_{1}=0.141\pm 0.002$ and $\lambda_{1}=0.131\pm 0.002$ (see supplementary information).
 Both results presented in Fig. 3 validate the main conclusions of the analytical model; for instance, the P2 bifurcation is mainly obtained with negative detuning for the low amplitude branch and positive detuning for the high amplitude branch, as predicted by the libration frequency analysis. Experimental and numerical results also demonstrate that the model can help constrain the necessary parameter space for chaos generation. It is equally interesting to note that the low-amplitude branch simulations have a well-formed wedge area for period doubling and chaos, and this is nearer to the analytically constrained parameter space compared to the high-amplitude branch. This can plausibly be attributed to the shape of the libration orbits around the low-amplitude branch, which more closely resembles circular ones compared to the almost banana-shaped high-amplitude branch libration orbits.

Next, we verify how closely the analytical and numerical results evolve as a function of changing drive conditions (changing $\delta$ or ${F_{1}}$). Based on the analytical model, it would be expected that, as the edge of the bistable area is approached, the distance between one of the stable fixed points and the unstable fixed point shrinks to zero, and as a consequence, the necessary ${F_{2}}$ required to achieve P2 and chaos itself is reduced to zero. This is confirmed numerically by performing 3D sweeps ($\delta$, $\Omega$, and ${F_{2}}$) and tracking the minimum necessary values of ${F_{2}}$ and $\Omega$ for the onset of P2. These are plotted against $\delta$ (for ${F_{1}}=3~{}\mathrm{V_{PP}}$), along with ${{\mathrm{F_{2min}}}}$ and $\mathrm{\Omega_{min}}$ as obtained from Eq. (3), in Fig. 4(b)-(c), respectively. Again, numerical results agree with the main features of the analytical model, in that both ${F_{2}}$ and $\Omega$ reduce to zero as the saddle node bifurcation is approached. Similarly to the results in Fig. (3), the analytical results for the low-amplitude branch adhere better to the numerical simulations compared to the high-amplitude branch ones, which show more exotic dynamics, a repeated hint that the eccentric libration orbit is more difficult to model.

It is interesting to ask whether the perturbed homoclinic-based argument presented above remains valid for the region where bistability is suppressed by the damping (note that the Melnikov approach requires bistability in the lossless version of the system). Indeed, the Melnikov method requires a homoclinic in the undamped system; however, it also assumes that the additional damping and forcing (corresponding to the terms within the brackets in Eq. 2) are small, perturbation-like terms. The authors therefore conjecture that, in the case of suppressed bistability, these terms are sufficiently large to invalidate the perturbed Hamiltonian approach.
 On a practical note, it maybe tempting to pursue chaos generation by positioning the drive near the saddle point bifurcations, as that would require a very small perturbation tone. This, however, is not an optimal experimental condition as it would be easy under the effect of noise to jump to the adjacent potential well and remain stuck there in a low-amplitude libration orbit. It is therefore more favourable for practical ends to position the main tone towards the middle of the bistable region and apply a moderate-amplitude second tone.
 In summary, this work presented an approach to chaos generation in nonlinear M/NEMS resonators that uses drive amplitudes that are typically an order of magnitude smaller than those previously reported. The method relies on applying two drive tones in order to make the rotating-frame phase-space three dimensional, and presents a model that uses linearization of perturbations within the rotating-frame to constrain the parameter space where chaos can be generated. This approach demonstrates that once bistability has been accessed via the first applied tone, chaos can in principle be generated with an arbitrarily small value of a perturbing second tone. The generality of the proposed method, and the relatively low driving forces involved underline applicability to a large range of nonlinear resonators thus potentially placing resonators, particularly M/NEMS ones, as leading candidates for the high-integration physical implementation of networks and reservoirs, as well as the experimental investigation of chaos-related phenomena.