A hybrid averaging and harmonic balance method for weakly nonlinear asymmetric resonators

Models for nonlinear vibrations commonly employ polynomial terms that arise from series expansions about an equilibrium point. The analysis of symmetric systems with cubic stiffness terms is very common, and the inclusion of asymmetric quadratic terms is known to modify the effective cubic nonlinearity in weakly nonlinear systems. When using low (second, in this case)-order perturbation methods, the net effect in these cases is found to be a monotonic dependence of the free vibration frequency on the amplitude squared, with a single term that depends on the coefficients of the quadratic and cubic terms. However, in many applications, such a monotonic dependence is not observed, necessitating the use of techniques for strongly nonlinear systems, or the inclusion of higher-order terms and perturbation methods in weakly nonlinear formulations. In either case, the analysis involves very tedious and/or numerical approaches for determining the system response. In the present work, we propose a method that is a hybrid of the methods of averaging and harmonic balance, which provides, with relatively straightforward calculations, good approximations for the free and forced vibration response of weakly nonlinear asymmetric systems. For free vibration, it captures the correct amplitude–frequency dependence, including cases of non-monoticity. The method can also be used to determine the steady-state response of damped, harmonically driven vibrations, including information about stability. The method is described, and general results are obtained for an asymmetric system with up to quintic nonlinear terms. The results are applied to a numerical example and validated using simulations. This approach will be useful for analyzing a variety of system models with polynomial nonlinearities.


Introduction
Recently, the vibrations of micro-and nano-electrome chanical resonators have provided a new realm of applications for nonlinear vibrations [1][2][3][4]. Many significant advances have been made in theory [5][6][7], experiments [8][9][10], and applications [11][12][13], including exploring new phenomena uncovered in previously unattainable parameter regimes [14][15][16]. The dynamics of these systems range from simple Duffing responses [17][18][19] to more exotic responses, sometimes involving two or more interacting vibration modes [20][21][22]. Even for single-mode nonlinear vibrations, a gap remains in the modeling and analysis of systems that exhibit higher-order nonlinear effects, particularly those resulting in non-monotonic amplitude-frequency characteristics [23][24][25]. Much of the present work of such systems relies on idealized models with reflection-symmetric stiffness involving both hardening and softening at different orders, the simplest version of which includes cubic and quintic nonlinearities with coefficients of different signs, which can be experimentally measured [26,27]. While these models adequately describe the frequency response of such systems, they are typically phenomenological and not fully descriptive of the underlying physics, which generically has asymmetric effects. That is, for asymmetric physical systems, the effective cubic and quintic coefficients of the symmetric model are not directly linked with the physics of the device. This gap leads to an inadequate connection between device physics and the dynamical models. Closing this gap allows one to better understand the sources of the nonlinear effects and predict how they depend on physical parameters.
The main issue in dealing with asymmetric systems that exhibit non-monotonic amplitude-frequency dependence is that the situation requires models that are cumbersome to analyze. In fact, to consistently capture the frequency response of such a system with polynomial nonlinearities, one must include all terms up to fifth order in stiffness and then carry out a higherorder analysis using one of the conventional techniques, say the method of multiple scales, which will involve four time scales-an extremely tedious task; see, for example, a recent study of an asymmetrical structural model [28]. Other similar methods that are equally taxing at higher orders are the methods of averaging and normal forms, each of which require several successive coordinate transformations to obtain consistent results for such systems [29]. Global approaches, for example, using action-angle coordinates, are valid for strong nonlinearities and can be employed to achieve the desired results [7,30,31], but these are less accessible to many researchers and often require numerical methods to execute, except in the case of simple potentials that yield results in terms of elliptic functions. Even in the relatively simple case with only quadratic and cubic stiffness terms, capturing non-monoticity requires going beyond second-order perturbation, but as shown below, a fully consistent result from a model with expanded stiffness terms requires all terms up to fifth order.
In this short paper, we introduce an approach that combines averaging and harmonic balance that allows one to determine the most essential features of the dynamic response of weakly nonlinear systems with non-monotonic amplitude frequency relationships. The method is a blend of a first-order perturbation calculation with a high-order harmonic balance (HB) calculation, carried out by amplitude expansions. The HB aspect of the method captures information about higher harmonic overtones and the constant offset from equilibrium. The results differ from what is obtained using the more involved calculations of standard techniques, but are accurate and attained with much less effort. We use the results to examine the conditions under which a system exhibits monotonic or non-monotonic frequency responses for asymmetric systems. These results are the missing link between physical models of asymmetric systems and the various types of frequency responses observed in many experiments. The price paid for these savings in effort is a lack of information about the transient response, as subsequently described.
To motivate our analysis, we present in Fig. 1 four generic forms of free vibration amplitude-frequency relationships, commonly referred to as backbone curves, that are frequently observed for a vibration mode. These are pure softening, pure hardening, mixed softening to hardening, and mixed hardening to softening. Often, one employs a simple effective cubic stiffness nonlinearity to describe the first two cases. For example, a well-known analysis for the monotonic backbone behavior for systems with quadratic and cubic nonlinearities can be found in classic physics and engineering books [32,33]. The results show that the system has an effective Duffing coefficient in which the quadratic term has a softening effect and the effect of the cubic term depends on its sign. Using secondorder perturbation methods, this combination results in a monotonic amplitude-frequency dependence that can be described by an effective cubic term, that is, a simple Duffing resonator. It is important to note that the Duffing coefficient depends on the level of asymmetry in the system and can vary as one varies the level of asymmetric bias. The description of a non-monotonic amplitude-frequency dependence (i.e., the latter two backbone curves in Fig. 1) with polynomial nonlinearities can be described by a simple model with reflection Fig. 1 Four generic forms of amplitude-frequency backbone curves: pure softening (red) γ < 0, σ ≤ 0, pure hardening (orange) γ > 0, σ ≥ 0, mixed softening to hardening (blue) γ < 0, σ > 0, and mixed hardening to softening γ > 0, σ < 0 (purple) symmetry, where the signs of γ and σ are opposite. In fact, this model can describe all four cases depicted in Fig. 1, as indicated in the caption. With this model, the backbone curves can be captured using a first-order perturbation approach, cf. [26], although a fully consistent result requires higher-order calculations, as demonstrated below. This description is adequate from a phenomenological point of view, but the symmetry of the nonlinear terms is inconsistent with the physics of asymmetric systems, which are far more generic than symmetric systems. The non-monotonic amplitude-frequency relationships can also be captured using higher-order perturbation methods, which can exhibit this feature even for systems with only quadratic and cubic nonlinearities; see [28]. In fact, in a consistent formulation, all the coefficients of lower-order nonlinear terms (quadratic, cubic, and quartic) will affect the effective quintic coefficient, analogous to the fact that the quadratic coefficient affects the effective cubic coefficient. A key feature of the present approach is that it gives accurate results using only a first-order perturbation method, coupled with higher-order harmonic balance.
It should be noted that perturbation methods such as the method of multiple scales and averaging are analytical formulations that provide convenient solutions for simple prototypical systems. However, for most systems, such methods must be complemented by numerics, for example, to solve for fixed points, compute eigenvalues, etc. This combination of analy-sis and numerics is powerful and certainly worthwhile in many cases. In other situations, brute force numerical methods are more convenient, and direct numerical simulation of the differential equations are often used to benchmark results of analytical approximations. However, parameter studies by direct simulation are very tedious. The HB method offers a combination that is generally more numerically intensive than perturbation methods, but is more convenient than simulations for parameter studies of periodic responses. When using HB, a series of nonlinear equations is generally solved for the steady state harmonic amplitudes, and stability is computed using numerical evaluation of eigenvalues of a monodromy matrix, or some equivalent numerical calculation. The literature on HB is vast. Here we point to a recent book [34] that covers topics relevant for the present work, and two open-source HB software packages, MANLAB [35] and HarmonicBalance.jl [36]. These packages are general and quite powerful.
The contribution of the present work is twofold: (i) the introduction of a new hybrid approach for analyzing the free and forced vibration for asymmetric systems, and (ii) general results for systems with up through quintic nonlinearities, which are used to categorize the possible frequency responses and how they depend on the nonlinear stiffness coefficients. The paper begins by introducing the method and provides its detailed execution for quintic systems. This is followed by a discussion of the results and a demonstration using an example of the non-monotonic response obtained from a resonator with quadratic and cubic nonlinearities.

The hybrid averaging method
We consider a model for a mode of vibration with up to fifth-order stiffness and damping nonlinearities, where ω 0 is the mode eigenfrequency, is the linear decay rate (damping), g(x) = 1+ 4 n=1 g n x n includes nonlinear damping effects, f (x) = 5 n=2 α n x n represents nonlinear stiffness, h(x) = 1 + h 1 x + h 2 x 2 allows for displacement-dependent excitation, and F and ω are the drive amplitude and frequency. We note that Eq. (2) is quite general and is a truncated version of the model that one gets from the microscopic theory of a resonator that interacts with a medium (see the The numerically computed spectral content (black curve) and analytical approximations (red, orange, yellow, green, blue, and purple) of the spectral components of a damped, driven resonator with quadratic and cubic nonlinearities. The system parameters are ω 0 = 10, α 2 = 9.64, Appendix of Ref. [37] and the references therein). This model can also arise from projection of a distributed parameter model onto a single mode. Here we use Eq.
(2) only as a prototypical example, and our method is not restricted to this model. To be specific, the functions g(x), f (x), and h(x) in Eq. (2) are intimately related via the fluctuation-dissipation theorem. Nevertheless, our method can be applied to other systems, where the functions g(x), f (x), and h(x) have different forms that can also depend onẋ.
We focus on energy levels that are associated with weakly nonlinear dynamics. The response of the system in this regime is dominated by the fundamental harmonic and can be expressed as x = a cos(ω 0 t + φ)+HOT. The HOT are the higher harmonic overtones and a constant offset (herein referred to as DC, from "direct current" in electronics) generated from the nonlinearities, and their energy is significantly lower than the energy of the fundamental harmonic (see Fig. 2).
To make the previous statement quantitative, we write the response as where the a k 's are polynomial functions of a that start from quadratic order, and are given by a 2k = 4 n=2k,n =0 c kn ( a) n for even harmonics, and by a 2k+1 = 5 n=2k+1 c kn ( a) n for odd harmonics. The bookkeeping parameter has been introduced to keep track of small quantities.
To capture the effects of linear damping, forcing, and nonlinearities near the primary resonance (ω ≈ ω 0 ) with a single slowly varying timescale, we make the following assumptions: (i) The amplitude and phase of the fundamental harmonic are slowly varying functions of time, specifically, such that a = a( 4 t) and φ = φ( 4 t). (ii) The damping, detuning, and forcing amplitude are small and can be scaled as = 4˜ , ω − ω 0 = 4 ω, and F = 5F . In the spirit of averaging, we assume the forṁ so that Eqs. (3) and (4) represent a time-varying change of coordinates from (x,ẋ) to (a, φ), one that requires the standard constraint equation, a cos(ωt + φ) − aφ sin(ωt + φ) = 0. (5) Note that there is no need for additional constraints on the zero (DC) and higher harmonics as their time derivatives are of higher order, i.e.,ȧ k | k =1 ∼ O( 6 ). From Eq. (4), we find thaẗ To obtain the dependence of the DC and overtone amplitudes on the primary harmonic amplitude a, we substitute Eqs. (3)-(6) into Eq. (2) and equate the coefficients of cos[k(ωt + φ)] for k = 0, 2, 3, 4, 5 to zero (this is HB), and for each harmonic, we equate the coefficient of the same power of to zero (amplitude expansions). Solving the resulting linear algebraic equations and recombining the expansion terms, we determine Note that, due to the scaling used, these expressions for the amplitudes of the DC term and the overtones are independent of the damping and forcing terms and are therefore solely relevant to capturing the DC offset and frequency shift from nonlinear stiffness effects. Also note that higher-order terms in are neglected in these expansions and are not indicated here. The above expressions are substituted into Eq. (3) and terms up to order 5 are retained. That result is used in Eq. (2), along withẋ from Eq. (4) andẍ from Eq. (6) to obtain a lengthy equation that involves a and φ, their time derivatives, and the coefficients from the model. That result, together with the constraint equation Eq. (5), provides the usual equations for the method of averaging. This pair of equations is solved forȧ(t) anḋ φ(t) and the results are averaged over the period 2π/ω, while assuming that a and φ do not change appreciably over that period. Setting = 1 in the result from this process yields the following equations that govern the slow dynamics of (a(t), φ(t)): where γ eff = α 3 − 10 9 represent the effective cubic and quintic coefficients that are typically used in models. A key point is that these expressions inform how they depend on the nonlinear coefficients of the original model. Note that Eqs. (8)- (9) are precisely what one obtains using firstorder averaging on the symmetric system, Eq. (1), with γ eff = γ and σ eff = σ [26]. For the case with quadratic and cubic terms, the expression for γ eff is well known and can be found in textbooks, e.g., [32,38]. Under the given assumptions, Eqs. (8)- (9) indicate that a and φ are slowly varying in time. However, an inspection of Eq. (9) reveals that there is an intermediate timescale The backbone curves exhibit mixed softening-to-hardening behavior, accompanied by a zero dispersion (ZD) point, denoted by the symbol , at which dω res /da res = 0. The locus of the ZD points is indicated by the dashed black curve. The position of the ZD point on the backbone curve changes considerably as α 2 is changed from α 2 = 9.488 (red backbone curve) to α 2 = 9.64 (purple backbone curve). All other parameters of the conservative system are held fixed (ω 0 = 10, α 3 = 1). The dashed translucent curves are calculated from the frequency ω exact of the exact analytical solution (see Ref. [41]), which is in excellent agreement with the frequency ω res (solid dark curves) of our approximate solution.
Right panel: Response curves for ω 0 = 10, α 2 = 9.64, α 3 = 1, F = 10 −2 , = F/(2ω 0 a res ), where a res is the resonance amplitude on the backbone curve (dashed black curve) that varies from a res = 0.2 (red response curve) to a res = 2 (purple response curve), which effectively varies the damping coefficient. Only when a res ≥ a zd = 1.069 does the response curve become nonmonotonic. The stable/unstable branches of the response curves are denoted by solid/dashed curves. The symbols and indicate the steady-state fundamental harmonic amplitudes reached during sweep-up (increasing ω) and sweep-down (decreasing ω), respectively, obtained by performing numerical time integration of Eq. (2) and Fourier analysis at the steady-state. Note that several of the and overlap proportional to a 2 , i.e., with the bookkeeping parameter φ = φ( 2 t, 4 t), which is neglected in the present analysis. If an accurate transient response is needed, the dynamics associated with this intermediate timescale can be captured with more sophisticated methods, such as multiple time scales or higher-order averaging [38]. Nevertheless, the present results provide valid information about the backbone curve and the steady-state response, including stability. Consequently, we can use this significantly simpler analysis that, strictly speaking, is restricted to a = a( 2 t → ∞, 4 t) and φ = φ( 2 t → ∞, 4 t) and take advantage of the fact that, while a(t) and φ(t) are still evolving in their slowest timescale 4 t, they reach their steady-state values in the intermediate timescale 2 t. As a result, the DC and overtones adiabatically follow the primary harmonic on the slowest timescale. In other words, the steadystate effects of the intermediate timescale are implicitly accounted for, and 2 t is only needed for the transient dynamics [see Appendix A:]. Relevant to this point is that the averaged rotating frame phase planes generated from Eqs. (8)- (9) do not accurately describe the dynamics except at, and very near, the fixed points. This is the price that is paid for the simplified calculations.
To determine the steady-state response of the system, we proceed in the standard way by setting (ȧ,φ) = (0, 0) in Eqs. (8)- (9). Eliminating the phase in the usual way, we seek solutions of the steady-state amplitude a ss from the following algebraic equation for which the attendant steady-state phase is given by In this way the method mimics the usual method of first-order averaging.
The resonance amplitude a res , is the largest amplitude on the response curve for a fixed value of F (see Fig. 3). This amplitude can be readily obtained from Eq. (11) and (13) and is given by a res = F/(2ω 0 ), which is known to be independent of the nonlinearity for models with only linear damping [39]. The corresponding resonance phase is φ res = −π/2 [40]. Thus, we find from Eq. (12) that the resonance frequency ω res satisfies the following backbone curve relation Furthermore, Eq. (13) also serves as an approximation for the amplitude-dependent frequency of free oscillation of the conservative system, i.e., when = 0, F = 0. The excellent agreement between this approximation for ω res and the exact frequency of free oscillations ω exact for a system with only quadratic and cubic nonlinearities is demonstrated in the left panel of Fig. 3, which demonstrates how the present method captures non-monotonic behavior even in this case. 1 If one were to use simple second-order perturbation methods, this system would predict a monotonic backbone curve, whereas the present approach captures the amplitude dependent softening of the quadratic term and the hardening of the cubic term. It should be noted that one must include quartic and quintic coefficients from a physical model to be strictly consistent to the same order, but their inclusion does not alter the approach nor the general conclusions.
We note that a non-monotonic response can only occur if the backbone curve has a zero-dispersion (ZD) point at which the derivative dω res /da res vanishes at 1 Note that the exact solution of the free oscillation given in Eq. (11) of Ref. [41] is obtained from an elliptic integral of the conservative system, and is expressed in terms of the total energy of the system E tot rather than in terms of the amplitude of oscillation a. However, the relation between these two quantities is simply E tot = ω 2 0 x 2 max /2 + α 2 x 3 max /3 + α 3 x 4 max /4, where x max = a + 5 k=0,k =1 a k , and therefore, we can easily compute ω exact in terms of a. a non-zero amplitude. From Eq. (13), we find that the amplitude and the frequency at the ZD point are given by a zd = √ (−3γ eff )/(5σ eff ) and ω zd = ω 0 − 9γ 2 eff /(80ω 0 σ eff ). For resonators with only quadratic and cubic nonlinearities, these expressions simplify to As shown in Fig. 3, the existence of a ZD point in the backbone curve is not the only requirement for a non-monotonic frequency response curve in the damped-driven case. The additional requirement is that the drive amplitude must be sufficiently large such that a res ≥ a zd . In particular, for resonators with only quadratic and cubic nonlinearities, this requirement is satisfied when To assess the local stability of the different branches of the response curve, we superimpose a perturbation δu = (δa, δφ) T on the steady-state response u ss = (a ss , φ ss ) T , carry out the corresponding substitutions in Eqs. (8)- (9), linearize in the perturbation terms, and obtain a pair of linear equations δu = J · δu, where J is the Jacobian matrix evaluated at u ss , ⎦ whose trace and determinant are given by respectively. For a non-conservative system, = 0, so that tr(J) = 0 and thus no Hopf bifurcations can occur [29]. However, static bifurcations in which det(J) = 0 can occur. From a simple analysis of the amplitude response curve given by Eq. (11), it can be shown that only saddle-node bifurcations, which correspond to an infinite slope in the response curve (da ss /dω → ∞), are possible [39]. Moreover, from Eqs. (11) and (15), we find the following simple relation must hold at a bifurcation point (ω sn , a sn ), From Eq. (16), we can calculate the critical drive level F cr at which saddle node bifurcations emerge. To calculate F cr , we set d F 2 /da 2 sn = 0, and find the critical amplitude a sn,cr from the resulting cubic equation in a 2 sn . Then, a sn,cr is substituted into Eq. (15) with the condition det(J) = 0 to eliminate the frequency dependency of a sn,cr , and finally, a sn,cr is substituted into Eq. (16) to yield F cr . For systems with σ eff = 0, that is, with only effective cubic nonlinearities, the above process significantly simplifies, and we obtain the well-known result F cr = [256ω 2 0 3 /(9 √ 3|γ eff |)] 1/2 [32].

Closing remarks
In this brief exposition, we introduced a method that combines higher-order harmonic balance and firstorder averaging. This technique provides a relatively simple approach to dealing with systems with higherorder polynomial nonlinearities when compared to conventional methods such as the method of multiple scales. The main assumption of the method is that the constant offset (DC) and overtone amplitudes and phases adiabatically track the amplitude and phase of the primary harmonic, whose dynamics are described by evolution equations that capture only the slowest timescale in the system. The results provide equations that govern the very slow dynamics of the primary amplitude and phase, along with predictions for the DC and overtone amplitudes in terms of the primary amplitude. This allows one to determine the backbone curve for free vibrations and the steady-state response for the damped, driven system, along with its stabil-ity characteristics. The assumptions and simplifications are restrictive only in terms of the transient dynamics, a minor inconvenience when considering the benefits of the approach. It should also be noted that one obtains distinct expressions for the higher-order frequency coefficient σ eff using the present method, higher-order MMS, and a series expansion of the exact integral for free vibrations (this can be most easily confirmed for the case with only quadratic and cubic nonlinearities); however, all three approximations give accurate results when compared to the exact integral.
The results illuminate why experimentally observed frequency responses show generic forms that can be described by equivalent symmetric systems with only cubic and quintic terms, and they provide the explicit connection between these forms and the coefficients from asymmetric physical system models. A typical example of such an asymmetric system is one that is nominally symmetric, but with a DC bias that breaks the symmetry [7]. Such bias can arise, for example, from gravity or, in the case of MEMS, a bias voltage. In such cases, the coefficients of the nonlinear stiffness terms in the dynamic model will vary with the level of bias and, in turn, the coefficients in the harmonic terms given in Eq. (7) and in the averaged equations Eq.(8)-(9) will also vary, providing an explicit connection between the physics and the predicted response. In this way, the present results will be of use for parameter identification in asymmetric systems [42]. Also, while not considered here, with appropriate scaling this method can be applied to parametric, subharmonic, superharmonic, and internal resonances of multimode systems. These topics will be considered in future studies.
Funding This work is supported by BSF under Grant No. 2018041. SWS is also supported by NSF grant CMMI-1662619. OS is also supported by the Pearlstone Center of Aeronautical Engineering Studies at Ben-Gurion University of the Negev. S.R. acknowledges the financial support of the Kreitman school of advanced graduate studies at Ben-Gurion University of the Negev under the Negev-Tzin Scholarship. Just prior to submission, we learned of a concurrent formulation of a very similar technique by Mark Dykman that will be used in a forthcoming paper; discussions with him helped to clarify and expand our understanding of the approach.
Data Availability Data will be made available by reasonable requests.

Conflicts of interest
The authors declare that they have no conflict of interest.

Appendix A: Transient dynamics
Here, we demonstrate the main shortcoming of the proposed approach, namely its inability to correctly capture the full transient dynamics. To this end, we numerically integrated both Eq. (2) and Eqs. (8)-(9) with identical initial conditions and compared their transient dynamics. To extract the amplitude a(t) and phase φ(t) of the fundamental harmonic from Eq. (2), we heterodyned the numerically obtained signal x(t) with the in-phase and quadrature components at the drive frequency ω, and passed the mixed signals through a bandpass filter to remove the overtones and DC components [26]. We then use the resulting slowly-varying (with respect to ω −1 ) quadratures of the fundamental harmonics, u(t) and v(t), to construct the amplitude and phase. The top panel of Fig. 4