The effects of nonlinear damping on degenerate parametric amplification

This paper considers the dynamic response of a single degree of freedom system with nonlinear stiffness and nonlinear damping that is subjected to both resonant direct excitation and resonant parametric excitation, with a general phase between the two. This generalizes and expands on previous studies of nonlinear effects on parametric amplification, notably by including the effects of nonlinear damping, which is commonly observed in a large variety of systems, including micro- and nano-scale resonators. Using the method of averaging, a thorough parameter study is carried out that describes the effects of the amplitudes and relative phase of the two forms of excitation. It is found that the frequency response can, in addition to the isolae and dual peaks that are known to occur in these systems, exhibit loops when nonlinear damping is present. The transitions among various topological forms of the frequency response are determined, and bifurcation analyses in parameter spaces of interest are carried out. In general, these results provide a complete picture of the system response and allow one to select drive conditions of interest that avoid bistability while providing maximum amplitude gain, maximum phase sensitivity, or a flat resonant peak, in systems with nonlinear damping.


Introduction
Parametric amplification (PA) is the use of a resonant parametric excitation to enhance the response of a resonantly driven oscillator. This approach allows one to alter the effective damping of the system, even to the limit of zero damping at the point of parametric instability, thus bringing benefits of spectral narrowing and higher frequency selectivity to resonant systems [1,2]. Specifically, the amplification, deamplification, and thermal noise squeezing have been analyzed for a Josephson parametric amplifier (JPA) [3,4]. In a classic study * lid2016@my.fit.edu † sshaw@fit.edu of a mechanical device, PA and thermomechanical noise squeezing were observed in a vibrating microcantilever and were analyzed using a linear model [5]. These studies were based on a linear model, which demonstrated the main effects. Studies on the impacts of stiffness nonlinearity on PA and similar systems have shown that the response can be quite rich [6,7,8]. The Duffing nonlinearity is especially of interest, as it is oftentimes exhibited in systems with large vibration amplitude, resulting in both opportunities for improved performance and challenges due to the complexity in dynamic responses [9,10].
PA is used in a wide variety of applications, especially in the realm of nano-and micro-electro-mechanicalsystems (N/MEMS) [11]. Due to their small size, N/MEMS devices have the advantages of high sensitivity, high frequency range, low power consumption, low noise, and excellent integration with electronics [12]. Specifically, PA is used for nano-scale applications including piezoelectrically pumped parametric amplifiers [13,14], carbon nanotubes [15], and graphene-based resonators [16]. Furthermore, it is also used in micro-scale applications including multi-analyte mass sensors [17], mass sensing arrays [18], Coriolis mass flow sensors [19], phase-modulated microscopy [20], parametric symmetry breaking transducers [21], parametron (a resonatorbased logic device) [22], and parametrically pumped thermomechanical-noise-driven resonators [2]. In addition, the effect of PA can also arise from mode coupling, such as 1 : 1 and 1 : 2 internal resonances. As an example of 1 : 1 resonance, the combined effect of Coriolis force and nonlinear resonance coupling is observed in MEMS vibratory gyroscopes, such as vibrating ring gyroscopes (VRG) [23] and disk resonator gyroscopes (DRG) [24], from which PA arises [25,26]. Many different kinds of mechanical structures have been considered to achieve PA, including torsional oscillators [27], metallized siliconnitride diaphragms [28], silicon disk oscillators [29], single-crystal silicon in-plane arch microbeams [30], stoichiometric silicon nitride (SiN) membranes [31], and double-clamped cantilever beams with a concentrated mass at the center [32]. Several analyses have been carried out in experiments, including parametric noise squeezing for a microcantilever [33], effects of geometric nonlinearity for a integrated piezoelectric actuation and sensing system [34], closed-loop stability in the presence of nonlinearity [35], and distortion of the actuation waveform by the displacement-dependent electrostatic nonlinearity [36].
Moreover, PA can also be employed in macroscopic devices. This can be achieved by using base excited cantilever beams [37], wherein the effects of cubic nonlinearity [38], cubic parametric stiffness [39], and parametric bistability [40] have been demonstrated. Other macroscopic systems include sheet metal plates [8], horizontal wind turbine blades in steady rotation enduring cyclic transverse loading [41], dual-frequency parametric amplifiers (DFPA) with macroscopic modular mass and linear voice-coil actuator [42], thin stretched strings carrying an alternating electric current in a non-uniform magnetic field [43], and doubly clamped strings [44].
Physics applications also constitute an important field for parametric amplifiers. Josephson parametric amplifiers are of particular interest [3,4], which have oftentimes been used for superconducting quantum interference devices (SQUID) [45,46]. For example, terahertz Josephson plasma waves amplified in a cuprate superconductor [47] and quasiparticles flowing through a superconductor-insulator-superconductor junction [48] are analyzed. PA has also been used to study the reheating of an inflationary Universe [49].
Parametric suppression, also known as deamplification, attenuation, or splitting in the spectrum, which can be achieved by modulating the relative phase between the direct and parametric excitations, has also gained interest in both theory [9,37] and experiments [1,27,29,31,38,50]. It can be utilized to enhance the phase response of an oscillator by increasing the steepness of the phase slope near resonance [51].
Recent theoretical analysis for PA include the effects of quadratic and cubic nonlinearities [7], dual frequency parametric amplifiers with quadratic and cubic nonlinearities [52], regular and chaotic vibrations with time delay [53], and frequency comb responses [54].
In this present work, nonlinear damping, also known as nonlinear dissipation or nonlinear friction, is taken into account, in addition to stiffness nonlinearity. This is of practical interest since nonlinear damping is frequently observed in a large variety of structures. For instance, nonlinear damping has commonly been observed in NEMS resonators based on carbon nanotubes [55], graphene [55,56,57], and diamond [58]. Likewise, it has also been observed in micro-structures including noncontacting atomic force microscope (AFM) microbeams [59] and MEMS clamped-clamped beams [60,61]. In addition, nonlinear damping has been observed in macroscopic mechanical systems, such as large-amplitude ship rolling motions [62], concrete structures [63], stainless steel rectangular plates, stainless steel circular cylindrical panels, and zirconium alloy hollow rods [64]. Nonlinear damping of a given mode can also result from mode interactions such as induced two-phonon processes [65] and internal resonances [57,66,67]. In addition to the experimental observations, many theoretical works have also been completed, covering the topics of the relaxation of nonlinear oscillators interacting with a medium [68], estimation using Melnikov theory [69], estimation using analytic wavelet transform [70], dynamic response to harmonic drive [60], and characterization using the ringdown response [61].
In this paper, a degenerate parametric amplifier with nonlinearities in both stiffness and damping is considered. One quantity of particular interest is the parametric gain, defined as the ratio of the peak amplitudes (near resonance) with and without the parametric pump [5], and expressed as wherer peak is the steady-state peak amplitude. When G > 1, the oscillator is amplified, and when G < 1, the oscillator is suppressed, by the parametric pump. The relative phase between the direct drive and the parametric pump, denoted by ψ, plays a pivotal role in the nature of the system response. It affects the parametric gain of the system, an important quantity that will be defined in Sect. 2.2. It also has a significant effect on the structure of the steady-state response curves. Two special values of the relative phase, −π/4 and +π/4, will be considered in detail, as these values provide the system with maximum and minimum parametric gains, respectively. This paper is organized as follows. In Sect. 2, we formulate the problem, preview the general effects of nonlinear damping on PA, and provide an analysis of how the PA gain is affected by nonlinear damping. In Sect. 3, we analyze the case where the relative phase provides maximum parametric gain (ψ = −π/4), for which a bifurcation analysis is presented, and provide a discussion of both steady-state and transient responses. In Sect. 4, we analyze the case where the relative phase provides minimum parametric gain (ψ = +π/4). In addition to the bifurcation analysis, a special condition corresponding to the infinite phase slope at resonance is also discussed. In Sect. 5, we consider the system with an arbitrary relative phase. Finally, some conclusions are drawn in Sect. 6.

Model
We consider a single degree of freedom system consisting of a weakly nonlinear, weakly damped oscillator with eigenfrequency ω 0 that is excited by both a near resonant direct drive at frequency ω ≈ ω 0 and by a parametric pump at frequency 2ω, and a relative phase between the drive and the pump. Specifically, in addition to the usual Duffing nonlinearity, we assume that the oscillator is also subjected to nonlinear damping.
The equation of motion for this degenerate nonlinear parametric amplifier is given bÿ where Γ 1 and Γ 2 represent the linear and nonlinear damping coefficients and are both assumed to be positive, ω 2 0 and γ denote the linear and nonlinear stiffness coefficients, λ indicates the amplitude of the parametric pump, f specifies the direct drive amplitude, ω dictates the drive frequency, and ψ describes the relative phase of the two drives. The effects of damping, nonlinearities, and drives are assumed to be small, in the sense that when both the amplitude and the eigenfrequency are normalized to O(1), all other coefficients are O(ε), where 0 < ε ≪ 1 is a small scaling parameter.

Averaged system
Since only primary and principal parametric resonances are of interest, it is assumed that ω is close to ω 0 . Thus, a nondimensional frequency detuning can be introduced as which is used to illustrate how the drive frequency deviates from the natural frequency, for example, during a frequency sweep. Under the stated assumptions, the amplitude and phase of x (t) will be slowly varying functions of time.
The method of averaging is employed to obtain the timeinvariant equations that govern these quantities. To obtain equations suitable for averaging, we first apply the van der Pol transformation where r (t) and φ (t) are the slowly-varying polar coordinates representing the amplitude and the phase of x (t). We next average time over one period, 2π/ω, and implement the detuning definition above, which converts Eq.
(2) into to the following approximate averaged equationṡ Under the stated assumptions, all terms on the right hand side are small so that r (t) and φ (t) vary slowly in time, consistent with the physical assumptions. All the subsequent analyses will be based on the autonomous dynamical system governed by Eqs. (6)- (7). The steady-state condition is obtained by solvingṙ = 0 andφ = 0 simultaneously, which provides solutions for the steady-state amplituder and phaseφ, representing a fixed point of the averaged dynamical system. The stability of a fixed point can be determined by the local Jacobian matrix. To facilitate further analytical calculations, the steady-state phaseφ is first eliminated from the steady-state condition, leaving a single equation for the steady-state amplituder. With only the leadingorder terms kept, this equation is given by Since Eq. (8) is quintic inr 2 when ignoring the trivial solutions, which have no physical meaning, it suggests that there can be a maximum of five fixed points for a given frequency. From Eq. (8), it can be seen that nonlinear stiffness and detuning appear only in the collective term (4σω 2 0 −3γr 2 ), which implies that in the amplitudefrequency space, the nonlinear stiffness has the effect of bending the response curves horizontally, yet has no effect on the amplitude or the topological structure of the curves, other than those effects associated with the bending of the curves. This term can be rewritten as the backbone curve equation which quantifies the amount that the curves bend horizontally as a function of the amplitude. This is, of course, the same backbone curve as that for the usual Duffing equation. The amplitude on the backbone curve can then be written asr b = ω 0 σ/ (3γ). For sufficiently prominent effects of nonlinear stiffness, the usual Duffing-type bistability can be exhibited (plus more, as will be seen subsequently). For the special cases of maximum parametric gain (ψ = −π/4) and minimum parametric gain (ψ = +π/4), the detuning parameter can be solved explicitly from Eq.
(8) as a function ofr. These expressions, given here, are convenient for plotting the steady-state response in the amplitude-frequency space, where the subscripts of the two plus-minus signs in each expression indicate their mutual independence, implying that there can be up to four different drive frequencies resulting in the same amplitude. Furthermore, the structure of this equation indicates that the frequency response can have up to four saddle-node bifurcations. The role of the backbone curve, given by Eq. (9), is evident in the conditions for the response curves given by Eqs. (10)- (11).
The points at which the response curves intersect the backbone curve can be expressed in closed form. To calculate these amplitudes, let Eq. (9) apply to Eq. (10) and Eq. (11), respectively. Each can be simplified into a factored form, given by Note that there are three factors and, in fact, at most three such points, in contrast to the usual Duffing equation with direct drive which has only a single such point. Some of these factors do not have physical meaning. These two equations will be discussed in detail later in Sect. 3.1 and Sect. 4.1, respectively.

Significance of nonlinear damping
We include nonlinear damping in the model for three important reasons: (i) it is commonly observed in experiments across many fields [55,56,57,58,59,60,61,62,63,64,65,66,67], (ii) it arises from fundamental microscopic considerations in micro/nano resonators [57,68], and (iii) it allows for the closure of the nontrivial response branches in parametric resonance by saddle-node bifurcations that can occur even near resonance (in fact, it can limit the response even in the absence of the Duffing nonlinearity). In order to make the third point, consider an oscillator being excited only parametrically, in which case the parametric instability threshold can be observed from Eq. (8) by consideringr = 0, which re-duces to This is the condition for the well-known Arnold tongue, which is a pitchfork bifurcation condition in the equations and corresponds to a period doubling in the original system. This stability condition is independent of nonlinearities since it relates to the linear stability of the trivial response. Note that the parametric instability threshold at zero detuning has the lowest value, denoted here as λ AT,0 = 4Γ 1 /ω 0 . If the oscillator with only linear damping is driven parametrically above λ AT,0 , the frequency response, as predicted by first order perturbation methods, has four branches that do not merge by saddle-node bifurcations at any frequency, even far from resonance where the perturbation analysis is not valid.
To see the effect of nonlinear damping, consider removing the nonlinear stiffness terms 3γr 2 from Eq. (8), as they do not display any effect pertinent to the range of possible amplitudes, but only to the frequencies at which specific amplitudes occur. In this case, the presence of nonlinear damping, Γ 2 > 0, always results in a finite am- (b) η versus nonlinear damping using a nondimensionalized Γ 2 for three levels of the parametric pump below the instability threshold λ AT,0 , where the solid lines represent ψ = −π/4 (maximum gain case) and the dashed lines represent ψ = +π/4 (minimum gain case) plitude, even at resonance, for Γ 2 is the coefficient of the highest order term inr. Another important effect of nonlinear damping is that it competes with nonlinear stiffness in terms of the system exhibiting Duffing bistability. To clarify this, we consider a simple Duffing oscillator with nonlinear damping, for which the oscillator will not undergo a cusp bifurcation (i.e., the condition for the onset of bistability [71]) if the nonlinear damping is greater than the following value [60,72] For a Duffing oscillator with only direct drive, nonlinear damping with Γ 2 > Γ * 2 eliminates any possibility of the system exhibiting Duffing bistability. For the system with both direct and parametric drives, however, the bistability may still be observed under conditions to be demonstrated in the following. However, as will be shown, the system can undergo a cusp bifurcation which terminates the Duffing bistability for sufficiently large direct drive.

Impact of nonlinear damping on parametric gain
To examine how the parametric gain, defined in Eq. (1), is affected by nonlinear damping, it is of interest to compare this gain with the gain of the linear system (Γ 2 = 0). For λ < λ AT,0 , where the parametric pump level is below the instability threshold, η is introduced to denote the ratio of the gains with and without nonlinear damping at the peak, defined by where G linear is the linear gain given in [5]. For ψ = ±π/4, this ratio can be obtained from Eqs. (12)- (13) (which are subsequently elucidated by Eq. (21) and Eq. (27), respectively), and written in the form of the depressed cubic equation whose real root is given by where From Eqs. (18)- (20), it can be seen that the ratio η is characterized by a single variable ξ, and this dependence is shown in Fig. 1a. As ξ increases, η decreases monotonically, indicated by the linear-log plot.
In a similar manner, as shown in Fig. 1b, the parametric gain attenuates as the nonlinear damping (nondimensionalized in the figure) increases. The impact of nonlinear damping on the ψ = −π/4 case, shown in solid lines, compared to the ψ = +π/4 case, shown in dashed lines, is significantly more pronounced, due to its larger response amplitude, for which nonlinear damping will be more prominent. For this same reason, for other values of the relative phase, the corresponding curve will lie between the solid and dashed lines. Moreover, it is also suggested by Eqs. (18)-(20) that all the curves in Fig.  1b will coincide with the curve in Fig. 1a under horizontal stretching. As a consequence, the steepness of the solid line near Γ 2 = 0 strongly depends on the system parameters, as suggested by Eq. (20). For instance, if the parametric pump level is close to the threshold, that is, 4Γ 1 − λω 0 ≈ 0, then η declines rapidly, indicating that nonlinear damping has a very strong impact on the parametric gain when operating near the threshold. For the majority of applications related to PA, this value of the relative phase is of primary interest since it generates the highest effective quality factor (as measured by the width of the resonance peak in linear systems with a parametric pump) [1]. This makes the oscillator optimally energy-efficient and frequency-selective when compared to other values of relative phases. One interesting phenomenon that can be observed in the frequency domain is that, when nonlinear damping is present, for sufficiently large levels of the parametric pump, an isolated response branch, called an isola, appears in the frequency response curve, as demonstrated below.

Steady-state response and local bifurcation analysis
The equation for the steady-state amplitude is given by Eq. (10). Fig. 2 shows the transition of the steady-state response curves in the amplitude-frequency space as the parametric pump level is varied for systems without and with nonlinear stiffness. As the pump level increases, an isola emerges under the main response curve, consisting of two branches bounded by a pair of saddle node bifurcations, as shown in Fig. 2b for a system with linear stiffness. Since the nonlinearity in stiffness bends the curves horizontally, when it is sufficiently large, the Duffing nonlinearity can lead to bistability in the main response branch, whose frequency range is bounded by another pair of saddle-node bifurcations, as shown in Fig.  2c, d. When the pump level leads to an isola and the nonlinear stiffness leads to bistability, four saddle-node bifurcations occur, as depicted in Fig. 2d. There have been several experiments where the observation of the isola is reported [21,22,40,44]. From Fig. 2, it is known that all the three possible amplitude local extrema occur on the backbone curve, given by Eq. (12). It can be seen that, in this equation, the first factor always has exactly one positive solution, which corresponds to the peak amplitude on the main response branch given by the depressed cubic equation The second factor in Eq. (12) may have up to two positive solutions, depending on the system parameters, which corresponds to the highest and the lowest amplitudes on the isola When there is no nonlinear stiffness (γ = 0), as shown in Fig. 3a, b, the dynamic response at zero detuning is symmetric about φ = π/4 + nπ. The peak amplitude has the steady-state phase ofφ = 5π/4, while both amplitude extrema of the isola have the steady-state phase ofφ = π/4.
The condition for which the isola appears/disappears, here referred to as the isola onset condition, can be obtained by simultaneously solving Eq. (22), ∂/∂r of Eq. Amplitudes on the backbone curve versus the direct drive level. The parametric pump level of the blue and orange curves are above λ AT,0 , while the green curve is below λ AT,0 . For both cases, amplification can be observed (22), and eliminatingr, which yields When f > f isola , there is no isola in frequency response, as shown in Fig. 2a, c. When f < f isola , there may exist an isola in frequency response, as shown in Fig. 2b, d. Note that Eq. (23) suggests that the pump level satisfies λ > 4Γ 1 /ω 0 . Therefore, the conditions for an isola to exist are λ > λ AT,0 and f < f isola , independent of the Duffing nonlinearity. It is important to point out that these conditions provide useful information regarding the selection of the drive levels needed to achieve maximum gain without encountering an isola. Fig. 3 demonstrates how these aforementioned amplitudes change by varying the two drive levels. The saddle-node bifurcations seen from the orange curves signify the appearance or the disappearance of the isola in the frequency response curve, which demonstrate the isola condition described by Eq. (23). Additionally, by comparing the blue and green curves (λ > 0) to the gray curve (λ = 0), the effect of PA can be seen. Moreover, the green curve in Fig. 3b has a steeper slope than the gray curve, suggesting that it can be used to significantly improve the sensitivity in force detection and signal enhancement.
Local bifurcation conditions are also of interest, as they provide information regarding regions of multistability in various parameter spaces. Bifurcations of both codimension one and codimension two occur in parameter planes of interest. Obtaining the saddle-node bifurcation conditions in the frequency domain is similar to that of the isola condition. They can be calculated numerically by simultaneously solving Eq. (10), ∂/∂r of Eq. (10), and eliminatingr. The cusp bifurcation condition, on the other hand, can be obtained by utilizing a third equation, ∂/∂r 2 of Eq. (10), and further eliminating σ. Obtaining these conditions allows for a detailed analysis of multistability in various parameter spaces. Fig. 4 provides examples of bifurcation diagrams in some of the parameter planes of interest. In Fig. 4a-c, the numbers 1,2,3 indicate the number of stable equilibria inside the respective parameter regions. The curves show the saddle-node bifurcation conditions, with orange color being associated with the isola and blue color being associated with the Duffing bistability. In these figures, γ > 0 is assumed. For the case of γ < 0, the figures are simply reflected about σ = 0. It can be seen from these figures that increasing the direct drive amplitude leads to the contraction and disappearance of the isola, while increasing the parametric pump level leads to the appearance and expansion of the isola and the Duffing bistability. Between Fig. 4a, b, however, there is a fundamental distinction. The former case corresponds to a system with a nonlinear damping coefficient that is less than the critical value given in Eq. (15), that is, Γ 2 < Γ * 2 , while in the latter case, Γ 2 > Γ * 2 . For Γ 2 < Γ * 2 , increasing the direct drive amplitude leads to the appearance and expansion of the Duffing bistability, while for Γ 2 > Γ * 2 , increasing the direct drive amplitude leads to the contraction and disappearance of the Duffing bistability.
Furthermore, it is also of interest to explore the parameter space pertinent to the parametric pump level and the direct drive amplitude, that is, (f, λ), as shown in Fig. 4d. The orange, dashed curve represents the isola for the blue curves at f = 0 in Fig. 4d, and for the blue curves at λ = 0 in Fig. 4d. From Eq. (25), it is clear that its denominator verifies the expression for Γ * 2 given by Eq. (15). In Fig. 4d, for Γ 2 > Γ * 2 , instead of eventually reaching the abscissa as f is increased, the curve of the cusp bifurcation condition increases monotonically, which is consistent with the phenomenon shown in Fig. 4b.

Phase portraits and a global bifurcation
The transient dynamics of the averaged system are governed by Eqs. (6)-(7), with ψ = −π/4. From Fig. 2d, it is known that the system can have either one, two, or three stable steady states. Samples of the three generic cases of the possible phase portraits are shown in Fig. 5. The transient process and the phase portraits for the system with other relative phases are very similar to these. Note that the system is globally bounded, due to positive damping, which is apparent from Eq. (6). Fig. 5 shows that for generic cases, there are only three distinct possibilities. The first case is to have a single fixed point, as shown in Fig. 5a. This equilibrium is globally, asymptotically stable. Fig. 5b evolves from Fig. 5a by adding a pair of fixed points, one stable and one unstable, via a saddle-node bifurcation; this is the appearance of the isola in the frequency response. In a similar manner, the case of Fig. 5c emerges from a second saddle-node bifurcation, resulting in three stable equilibria and two saddle points; this is the first saddlenode on the Duffing response branch. The stable (orange) and unstable (blue) manifolds of the saddle points are shown in the figure. The stable manifolds serve as separatrices, which form boundaries among the basins of attraction of the stable fixed points. The unstable manifolds of saddle points, on the other hand, are asymptotic to the stable fixed points. As the detuning is continually increased, the isola fixed points will remerge in another saddle-node bifurcation and then two of the remaining fixed points will undergo the typical Duffing saddle-node bifurcation. This full transition, from a single stable response through the four saddle-node bifurcations to another single stable response, may require a global bifurcation, as considered next.
For the situation in which there are two saddle points, it is worth noting that a global bifurcation will occur, as shown in Fig. 6. This global bifurcation is topologically required in order for the phase portraits to transition as they do as parameters are varied. The bifurcation is a standard planar saddle connection, and the transition across this bifurcation alters the basins of attraction of the stable steady states.

Analysis for minimum parametric gain (ψ = +π/4)
The minimum parametric gain occurs for the relative phase of ψ = +π/4. In certain contexts, this is known as parametric suppression because the parametric gain can be less than unity. This special value of the relative phase is also of interest in certain applications, for example, to enhance the phase sensitivity near resonance, due to the strong impact that the parametric pump has on the steepness of the phase slope at resonance [51].
Just as increasing the parametric pump level can lead to the emergence of an isola in the frequency response, as discussed in Sect. 3.1, here it results in the appearance of a loop, via the initial formation of a dimple in between a pair of dual peaks. At the transition point between a dimple and a loop, the frequency response curve exhibits a cusp, and if the stiffness nonlinearity is zero, the phase slope will be infinite at resonance, a feature of interest.

Steady-state response and local bifurcation analysis
The equation for the steady-state amplitude is given by Eq. (11). Fig. 7 shows the transition of the frequency response curves as the parametric pump level is varied for a system with zero nonlinear stiffness. As the pump level is initially increased, the peak transforms into a dimple which is between a pair of peaks with equal amplitudes. When the pump level further increases, the dimple becomes more acute and transforms into a cusp. (This is not to be confused with a cusp bifurcation.) Beyond this, the cusp transforms into a loop. It should be pointed out that the two steady states at the crossing point in the loop have have distinct phases. When introducing the nonlinear stiffness to the system, these response curves may exhibit Duffing bistabilities, yet neither the amplitudes nor the general structure change, which is similar to what has been shown in Fig. 2.
When there are dual maxima, as shown in Fig. 7b-d, a condition for the amplitude of the dual peaks can be obtained by setting the inner radicand of Eq. (11) equal to zero, given by It is apparent that all other possible amplitude extrema are on the backbone curve. It can be seen that the first factor of Eq. (13) always has exactly one positive solution, which corresponds to the amplitude of the peak or the dimple on the response curve, as The third factor of Eq. (13) can have up to two equal positive solutions, depending on the system parameters, which corresponds to the amplitude of the cusp point shown in Fig. 7c or the amplitude of the crossing point shown in Fig. 7d, given bȳ It can be seen that the amplitude of the crossing point is independent of the direct drive amplitude. When the stiffness nonlinearity is zero, as shown in Fig. 7, the dynamic response at zero detuning is symmetric about φ = 3π/4 + nπ. The amplitude extremum on the backbone curve has the steady-state phase of The condition for which the response curve transitions from a single peak to dual peaks and a dimple is calculated by setting both radicands of Eq. (13) equal to zero, yielding This condition is interesting since it provides a relatively flat frequency response at resonance, that is, a minimal sensitivity to the amplitude to variations in the frequency.
The cusp condition can be obtained in a manner similar to that used for determining the isola condition. Specifically, it is obtained by simultaneously solving Eq. (13), ∂/∂r of Eq. (13), and eliminatingr, which yields For f ≥ f flat , the frequency response has only a single peak, as shown in Fig 7a. For f cusp < f < f flat , it has a dimple and dual peaks, as shown in Fig 7b. For f = f cusp , the dimple has become a cusp, as shown in Fig  7c. Finally, for f < f cusp , it has is a loop, as shown in Fig 7d. Similar as before, in order for the response curve to form a cusp or a loop, it is implied that λ > λ AT,0 . Fig. 8 demonstrates how these aforementioned amplitudes change by varying the drive levels. The intersection between the orange curves and the purple curves demonstrate a supercritical pitchfork bifurcation at zero detuning. In Fig. 8 (a), as the parametric pump level is increased, the response curve transitioning from a single peak to a dimple and dual peaks, and finally a loop can be clearly seen. In Fig. 8 (b), it can be seen thatr cross is independent of the direct drive amplitude. Additionally, by comparing to the gray curve (λ = 0), parametric suppression is also shown from the blue and orange solid curves.
The process of obtaining the local bifurcation conditions is also very similar to that described in Sect. 3.1. Fig. 9 provides examples of bifurcation diagrams in some two parameter spaces of interest. In Fig. 9a-c, the numbers 1,2,3 indicate the number of stable equilibria in the respective parameter regions. The curves show the saddle-node bifurcation conditions, with orange being associated with the dimple or loop, and blue being associated with the Duffing bistability. It can seen that increasing the direct drive amplitude leads to the contraction and disappearance of the dimple/loop, while increasing the parametric pump level leads to the appearance and expansion of the dimple/loop and the Duffing bistability. The distinction between Fig. 9a, b) is also very similar to that of Fig. 4a, b, where the former case corresponds to Γ 2 < Γ * 2 , while the latter case corresponds to Γ 2 > Γ * 2 , where Γ * 2 is given in Eq. (15). For Γ 2 < Γ * 2 , increasing the direct drive amplitude leads to the appearance and expansion of the Duffing bistability, while for Γ 2 > Γ * 2 , increasing the direct drive amplitude leads to the contraction and disappearance of the Duffing The bifurcations in the parameter space involving the parametric pump and the direct drive are shown in Fig.  9d. The orange curve represents the cusp bifurcation condition associated with the dimple or loop, the blue curve represents the SN bifurcation condition associated with the Duffing bistability, the magenta dashed curve indicates the appearance of the dimple, which is a condition described by Eq. (30), and the green dashed curve indicates the appearance of the loop, which is a condition described by Eq. (31). The intercepts of the blue curve are the same as those provided in Sect. 3.1, which are given by Eqs. (24)-(25).

Infinite phase slope condition
In this section, we consider the case where the system parameters satisfy the condition given by Eq. (31) and without nonlinear stiffness, which is shown in Fig. 7c. For the present analysis, we assume that the stiffness nonlinearity is zero. The conditions at the cusp point are given in Sect. 4.1 asr cusp = (λω 0 − 4Γ 1 ) /Γ 2 and φ cusp = −π/4. This is a case of special interest because it provides infinite phase slope versus the frequency detuning and therefore infinite sensitivity to fluctuations in phase. This is shown in Fig. 10a, which can be used to maximize phase sensitivity at resonance. This is also the supercritical pitchfork bifurcation condition, as demonstrated in Fig. 8.
In order to analyze the dynamics at the cusp point, it is convenient to use an alternative formulation for the averaged equations, in particular, a Cartesian form, in contrast to the polar form in Eqs. (4)- (5). This formulation, which leads to simpler calculations in this case, is given by x (t) = −ω 0 a (t) sin (ω 0 t) + ω 0 b (t) cos (ω 0 t).
The phase plane is shown in Fig. 10b. With the standard center manifold approach, v = h (u) is computed to the leading-order term where the coefficients for all odd-order terms are identically zero due to symmetry. The dynamics on the center manifold is then governed bẏ where all the even-order terms vanish due to symmetry. Because the condition given by Eq. (31) must satisfy λ > λ AT,0 , the sign ofu near u = 0 is negative, proving that the dynamics of the system on the slow manifold is stable near this equilibrium. Therefore, when operating at this special condition with infinite phase slope, the system is weakly, that is, nonlinearly, dynamically stable.

Analysis for arbitrary relative phase
The relative phase between the direct drive and the parametric pump (ψ) is of vital importance. Not only does it affect the structure of the frequency response curves and bifurcation diagrams, it also determines the parametric gain of the system. The parametric gain, as described  Figure 12: Ratio of the amplitudes of local extrema with various parametric pump levels to those values of extrema without the pump versus the relative phase ψ, with Γ 1 = 0.005, ω 0 = 1, and f = 0.01. The highest ratio in each color is considered to represent the parametric gain of the system. The stability of the curves is not indicated. Note that the green curves in both diagrams have the parametric pump level of λ AT,0 , the instability threshold. (a) Γ 2 = 0.005. (b) Γ 2 = 0 in Eq. (1), reflects how much the parametric pump can amplify or suppress the peak response amplitude. Fig. 11 shows frequency response curves for sample values of the parametric pump level and relative phase. Generally, for ψ = π/4 + nπ/2, isolae, loops, and dimples can be off-centered and, as parameters are varied, the loops can pinch off to form isolae. While the structure of the response curves in the amplitude-frequency space can be diverse, there are still limitations on the possibilities. Specifically, there can only be up to one isola and one dimple or loop in the response curve. In addition, as mentioned in Sect. 2.1, the maximum number of fixed points is five, in which case three are stable equilibria and two are saddle points. Fig. 11a, b illustrate how an the isola becomes tangent to the main response branch and forms a loop. Fig. 11c, d illustrate how a dimple transforms into a loop via a cusp. Fig.  11e, f illustrate how a dimple is formed via an inflection point. Considered together, these transitions describe the complete evolutionary possibilities for the frequency response curves. Fig. 12a provides an overview of the topology of the frequency responses that occur as one varies the relative phase ψ and how they change from one form to another. This figure shows the amplitudes of the local extrema that exist over the entire resonant frequency range for various parametric pump levels. In the figure, these amplitudes are normalized by the corresponding amplitudes that would occur without the pump.
In order to describe the possible transitions in the frequency response curves, we start our discussion focusing on the blue curve, corresponding to λ = 0.04, a portion of which is shown in the inset of Fig. 12a, and relate these to structures in the frequency response. This case shows all possible transitions and sets the stage for more special cases. First, for ψ near −π/4 there is always a single extrema for this case and for pump levels below the isola onset condition given by Eq. (23), indicating a simple frequency response. Increasing the phase from this point, one encounters point A as shown in Fig. 12a. This point, a fold in the language of catastrophe theory [74], corresponds to the appearance of the isola, initially formed as a point (at A) and then results in two new extrema associated with the isola, which exists between points A and C. When point B is encountered, the frequency response develops an inflection point that is separate from the isola, as depicted in the inset of Fig. 11e (for different parameter values), which then develops into a dimple. Such a dimple, without the isola, is shown in Fig. 11f. This dimple then forms into a cusp, as shown in the inset of Fig. 11c, which then transitions into a loop. (Note that this cusp transition, at the condition given by Eq. (31), does not alter the number of extrema and thereby is not indicated in Fig. 11.) Such a loop, without the isola, is shown in Fig. 11d. Therefore, between points B and C there exist five extrema. For this pump level, three of these are very close to one another, as shown in the inset of Fig. 11b. Here two extrema are associated with the isola (including the rightmost and uppermost of the three shown in the inset) and two are from the loop (which is very small in the inset). At point C, the two rightmost and uppermost extrema shown in the inset of Fig. 11b merge, the outcome of which is a large loop that replaces the isola. This process is reversed as ψ is increased, due to the symmetry of the diagram about ψ = +π/4. Note that for λ = 0.04, the transitions near points B and C occur very rapidly as ψ varies.
For larger pump levels, for example, λ = 0.05 shown in Fig. 12a, point A no longer exists and the main isola or associated loop exists for all values of ψ. As the pump level increases from 0.04 to 0.05, point A moves towards ψ = −π/4, eventually reaching it at which point it merges with its symmetric counterpart at the condition given by Eq. (23), and then disappears. In this case, there are three extrema for all values of ψ, corresponding to either an isola or a dimple/loop. Points B and C still exist, resulting a dimple/loop that exists around ψ = +π/4, replacing the isola over a range of phases.
For smaller pump levels, for example, λ = 0.03 shown in Fig. 12a, point A also does not exist, but in this case it corresponds to the complete absence of an isola. Consequently, point C does not exist, either. In this case, there is a single extremum for ψ = −π/4, and as ψ is increased an inflection point is encountered, resulting in a dimple/loop that exists around ψ = +π/4, corresponding to three extrema. Here there is a swallowtail structure [74] in the diagram. As the pump level is further decreased, the swallowtail and the attendant range of ψ over which the dimple/loop exists continues to shrink until it eventually disappears at the inflection point condition given by Eq. (30), in which case the swallowtail no longer exists and there is a smooth response curve with a single maximum for all values of ψ.
An additional feature of Fig. 12a is that it demonstrates that ψ = −π/4 and ψ = +π/4 always correspond to the maximum and the minimum parametric gains, respectively. This can been seen since such amplitude ratio (or the largest amplitude ratio if multiple ratios exist) represents the parametric gain for a given value of the pump level and relative phase, and the absolute maximum always occurs at ψ = −π/4 while the minimum of the largest branch occurs where the branches cross at ψ = +π/4.
In contrast, Fig. 11b shows the parametric gain in the absence of nonlinear damping. For lower levels of the pump, the response shares features with those of the lower pump values shown in Fig. 12a, including the swallowtail. However, as the parametric pump level reaches and exceeds the instability threshold (shown in the green curve), the resonance peak no longer exists across the entire domain of the relative phase.

Conclusion
This paper describes in detail the frequency response of systems with nonlinear damping and stiffness subjected to both direct and parametric near-resonant driving. This generalizes and expands on previous work in the area to include the effects of nonlinear damping, which is relevant to PA, for example, in micro/nano-scale resonators.
The results provide a thorough description of the possible types of frequency responses that can be encountered and how these depend on system and drive parameters. Of particular interest are how isolae, which were known to occur in these systems [9] and have been experimentally observed [21,22,40,44], are formed as the drive parameters, such as the relative phase ψ, are varied. New features for these systems are also described herein, including loops and dimples that are closely related to the isolae, and amplitude response curves with degenerate flat resonance peaks. The analysis provides a complete description of how these features are related via transitions involving cusps, inflection points, and tangencies. These results are useful, for example, by allowing one to select drive conditions that provide maximum amplitude gain without bistability, maximum phase sensitivity at the desired vibration amplitude, or flat resonance peaks with potential application for reducing amplitude noise, for given device parameters.
It is expected that the features described in this paper can be experimentally observed in any lightly damped resonator with the described characteristics, which are quite generic. It is also of interest to consider the closedloop, self-oscillating version of this system, with noise included in the model. Of particular interest is the behavior of noise near the transition points, for example, the flat resonance peak, where the effects of certain noises on the system response might be amplified or attenuated.