Multistable internal resonance in electroelastic crystals with nonlinearly coupled modes

Nonlinear modal interactions have recently become the focus of intense research in micro- and nanoscale resonators for their use to improve oscillator performance and probe the frontiers of fundamental physics. However, our understanding of modal coupling is largely restricted to clamped-clamped beams, and lacking in systems with both geometric and material nonlinearities. Here we report multistable energy transfer between internally resonant modes of an electroelastic crystal plate and use a mixed analytical-numerical approach to provide new insight into these complex interactions. Our results reveal a rich bifurcation structure marked by nested regions of multistability. Even the simple case of two coupled modes generates a host of topologically distinct dynamics over the parameter space, ranging from the usual Duffing bistability to complex multistable behaviour and quasiperiodic motion.


Equations of motion for finite deformations
The nonlinear theory of continuous media has been developed at great length by a number of authors [1][2][3] . Here, we collect some of the main results and sketch a derivation of the equations of motion for an electroelastic crystal. We represent the material body as a Euclidean manifold embedded in E 3 . In order to describe finite deformations we must distinguish between the reference (undeformed) and current (deformed) configurations. We refer these to Cartesian systems of rectangular coordinates with bases {E K } and {e i }, where upper and lowercase indices, respectively, serve to distinguish material and spatial coordinates. The motion of a material point X is given by x i = x i (X K , t) or inversely X K = X K (x k , t). The deformation gradient is then defined by where a comma denotes partial differentiation. The determinant of the deformation gradient is given by J = det x i,K . The Green and Lagrangian strain tensors are then defined, respectively, as where u is the displacement vector that extends from a material point P in the undeformed body to its spatial position in the deformed body, δ ij is the Kronecker delta, and the summation convention over repeated indices is assumed hereafter. The displacement is related to the reference and current positions by u K = δ iK x i − X K , where δ iK is the shifter that relates components of vectors in the material and spatial coordinates. When the two frames are identical the shifter reduces to the Kronecker delta.
We now collect some key results from a series of papers by Tiersten (see refs 31-33 in the main article) and enforce a notational consistency useful in what follows. The papers (1) derive the equations of electroelasticity based on a macroscopic model of an electronic charge continuum coupled to a lattice continuum and (2) proceed to expand these equations to cubic order in the displacement and electric potential gradients. This lays the groundwork necessary to treat weakly nonlinear dynamics of an electroelastic crystal.
Applying the conservation of mass and linear momentum, along with Maxwell's equations under the quasistatic electric field approximation, we get where ρ is the mass density, ρ 0 is the mass density in the reference configuration, v j = dx j /dt is the velocity, τ ij is the Cauchy stress tensor, and P i , E i , D i , and φ are the the electric polarization, field, displacement, and potential, respectively. The first equation above can be recast in terms of the electrostatic Maxwell stress tensor Now, by defining a thermodynamic state function χ = E − E i P i /ρ, where E is the stored energy per unit mass, the conservation of energy can be applied to show that Furthermore, by applying invariance principles it can be shown that χ = χ(E KL , W L ), where To put these equations in a useful form we need to convert them to the reference are symmetric mechanical and Maxwell electrostatic stress tensors, respectively, we arrive at where K Lj = F Lj + M Lj is the total first Piola-Kirchhoff stress tensor, whose mechanical and electrical contributions are defined as Piola transforms of τ S ij and T ES ij .
and D L is the Piola transform of electric displacement.
Finally, to render the equations tractable we approximate χ using a polynomial expansion about the reference configuration: where c ABCD , e ABC , χ AB , c ABCDEF , k ABCDE , b ABCD , χ ABC , and c ABCDEF GH are the secondorder elastic, piezoelectric, second-order electric permeability, third-order elastic, first odd electroelastic, electrostrictive, third-order electric permeability, and fourth order elastic constants, respectively.

Weakly nonlinear resonance
We can now expand the equations to cubic order in the displacement gradient and to linear order in the electric potential gradient. Terms involving the electric potential are kept linear because of the relatively small piezoelectric coupling of quartz. First, note that in this case the total Piola-Kirchhoff stress tensor K Lj reduces to F Lj , since M Lj is quadratic in the electric field. We now let the two Cartesian coordinate systems that describe the reference and current configurations coincide. This means that the shifter δ jM is simply the Kronecker delta. Carrying out the expansion to cubic order in u A,B and to linear order in φ ,A gives where KL = 0 δ KL + χ KL is the electric permittivity tensor. To simplify these equations we can define effective material constants in terms of fundamental material constants: We have developed Mathematica code to expand out a given component of K LM . It is a general code that accommodates tensors of arbitrary rank and automatically symmetrizes them and converts their components to condensed notation. For essentially thickness-shear vibrations in thin plates the boundary conditions on the minor surfaces can be neglected, and in many applications of Y-cut quartz K 21 and K 22 are the relevant stress components 6 .
In this case only the u 1,2 and u 2,2 displacement gradients are retained: Now, for our specific application of thickness-shear resonance u 1 dominates and we consider only K 21 , which reduces to Invoking equation (1) and accounting for boundary conditions at the faces of the plate and V and ω are the amplitude and angular frequency of the voltage applied across the crystal thickness. Hats have been added to denote dimensional quantities and will be removed after nondimensionalization. A linear damping term with coefficientr 66 has been added to account for the small losses of quartz 7 . If we introduce nondimensional quantities defined by then the nondimensional equations are where we have dropped the subscript on u 1 and It is worth noting that

Perturbation analysis
To obtain equations governing the complex modal amplitudes we employ the method of multiple scales. We seek the solution of equation (6) in the form where 1 is a nondimensional bookkeeping parameter and T k = k t, such that In order that the damping, nonlinearities, and resonances balance each other, we scale r 66 , Γ, and F as r 66 , Γ, and F . Substituting equations (7) and (8) into equations (6) and equating coefficients of like powers of , we obtain Order 0 : Order 1 : where the prime indicates differentiation with respect to y. Equations (9) are the linear homogeneous boundary-value problem associated with equations (6). Its general solution can be expressed as where c.c. means complex conjugate of the preceding terms and the spatial eigenmodes are Substituting equations (11) into equations (10) gives where NST stands for terms that do not produce secular terms. Since the homogeneous part of equations (12)-(15) has a nontrivial solution, the nonhomogeneous system has a solution only if solvability conditions are met 8,9 . Imposing the solvability conditions gives equations governing the complex modal amplitudes A n and A m :  (14) and (15) are first shifted onto equations (12) and (13). In this case the solvability conditions require that the right-hand sides of equations (12) and (13)  Note that ω n γ nm = ω m γ mn , which is necessary if equations (16) are to be derivable from a Lagrangian.
We convert equations (16) to a set of real-valued autonomous ODEs via the transforma- and set = 1, which gives p n (t) = γ nn q n p 2 n + q 2 n + γ nm q n p 2 m + q 2 m + γ n0 p n (q m p n − 2p m q n ) − q m q 2 n − µ n p n − ν 1 q n q n (t) = f − γ nn p n p 2 n + q 2 n − γ nm p n p 2 m + q 2 m + γ n0 p m q 2 n − p n (p m p n + 2q m q n ) − µ n q n + ν 1 p n p m (t) = γ mm q m p 2 m + q 2 m + γ mn q m p 2 n + q 2 n + γ m0 q n 3p 2 n − q 2 n − µ m p m − ν 2 q m q m (t) = −γ mm p m p 2 m + q 2 m − p m γ mn p 2 n + q 2 n + γ m0 p n 3q 2 n − p 2 n − µ m q m + ν 2 p m (17) where ν 1 = σ 1 and ν 2 = 3σ 1 − σ 2 . The modal amplitudes are then given by a l = p 2 l + q 2 l Finally, before carrying out numerical continuation we redimensionalize these equations as follows:

Numerical continuation
All numerical continuation results were obtained using the AUTO software. Plots of the primary and coupled mode amplitudes (a n and a m ) were obtained by simulating equations (17), after redimensionalization via equations (18). The amplitudes were extracted as a l = p 2 l + q 2 l , l = n, m. In plots containing both equilibria and periodic solutions (e.g. Fig. S1 below), the L 2 norm is used. From the AUTO-07p manual this is calculated as where the U k (t) are the p n , p m , q n , and q m in equation (17). Note that the integration limits of [0,1] occur since AUTO transforms the independent variable to range from 0 to 1 for periodic solutions.
As explained in the main article, the finite dimensions and electrode geometry of experimental crystals cause the thickness-shear resonances to couple with lateral eigenmodes.
This does not appear to cause qualitative changes in the bifurcation structure for AT-cut quartz, as shown in analyses using approximate Mindlin plate equations that account for lateral resonances [10][11][12][13] . However, it does cause multistable behavior to occur at lower driving voltages than predicted theoretically. Since the theoretical results also predict multistable behavior at progressively lower voltages as σ 2 decreases (see Fig. 4 in the main article), we used an effective value for σ 2 of 1.8 KHz to compute our theoretical frequency responses, instead of the experimentally determined value of σ 2 = 2.25 KHz. As a result, the theoretical drive voltages more closely match those at which multistability occurs in experiment.
Accounting for the effect of Joule heating would also slightly modify the effective spring constants of each mode, and in turn the effective value of σ 2 . Even without these extra considerations the agreement is still quite close, though in future work we intend to incorporate the effects of Joule heating and coupling with lateral eigenmodes.

Homoclinic orbits and quasiperiodic motion
One advantage of numerical continuation methods is the ability to continue periodic solutions emanating from Hopf bifurcations of equilibria. Figure S1 demonstrates this for the isolated loop of a representative frequency response. The periodic branches undergo a series of period doubling bifurcations and rapidly increase in amplitude. This accounts for the |S 11 | dips in the experimental response (see Fig. 5 in the main article). As the amplitude of the periodic branches increases, their period likewise increases, effectively tending to infinity.
Such behavior indicates the presence of homoclinic orbits. Furthermore, the continuation detects a torus bifurcation, and both suggest the existence of complex dynamics in the region between the two Hopf points. While these results anticipate the more complex transient behaviour witnessed in Fig. 5, it should be emphasized that the specific nature of the transient dynamics cannot be directly read from Fig. S1. Proving the existence of periodic orbits only tells us where more complex transient behavior might occur and not exactly how it will manifest. Specifically, we know such behaviour should occur in the region bounded by the two Hopf points in Fig. S1. Indeed, the experimental results corroborate these predictions as noted in the main article. Lastly, the main purpose of the model is to predict the steady state solutions and whether multistability is possible. The further prediction of experimentally observed transient behaviour is an added bonus, but future work is needed to unravel the specifics.

Experimental Setup
Signal Generator Spectrum Analyzer DUT (Quartz Resonator) Figure S2: Diagram of experimental setup. The spectrum analyzer uses built-in software to control the signal generator and collect the frequency response of the device. An RF circulator routes the signal from the source to the quartz resonator and the reflected signal is directed to the spectrum analyzer.
A schematic of the test setup is given in Fig. S2. To measure the frequency response of the quartz resonator the signal generator is controlled by the spectrum analyzer. A single sinusoidal signal is swept around 30 MHz at various drive voltages to produce the plots in Figs. 1a, 2 and 5 of the main article. For large enough driving amplitudes two modes are excited: the directly forced thickness-shear mode near 30 MHz and the internally resonant mode near 90 MHz (i.e. a 1:3 internal resonance between the third and ninth harmonics).