Frequency-comb response of a parametrically-driven Duffing oscillator to a small added ac excitation

The study of the effect of small signals on non-autonomous nonlinear dynamical systems near bifurcation points was pioneered by K. Wiesenfeld and B. McNamara in the mid 80’s Wiesenfeld and McNamara (1985, 1986). They showed that several different dynamical systems are very sensitive to coherent perturbations near the onset of codimention-one bifurcations, such as period doubling, saddle node, transcritical, Hopf, and pitchfork (symmetry-breaking) bifurcations. They developed a general linear response theory, based on perturbation and Floquet theories, explaining the effects of small coherent signals perturbing limit cycles of nonlinear systems near the onset of bifurcation points. One of the systems to which they applied their theory was the ac-driven Duffing oscillator. It was found by them that nonlinear dynamical systems could be used as narrow-band phase sensitive amplifiers.

Parametric amplification has been studied in electronic systems since at least from late 50’s and early 60’s by P. K. Tien Tien (1958), R. Landauer Landauer (1960), and Louisell Louisell (1960). It has been used for its desirable characteristics of high gain and low noise Batista and R. S. N. Moreira (2011). Parametrically-driven Duffing oscillators have been used to model many different physical systems such as the nonlinear dynamics of buckled beams Nayfeh and Sanchez (1989); Abou-Rayan et al. (1993) or driven Rayleigh-Bénard convection Lücke and Schank (1985). Homoclinic bifurcations were found in the parametrically-driven Duffing oscillator Parthasarathy (1992). Further bifurcations and chaos were found in Jiang et al. (2017). None of these papers investigated the type of spectral response we study here.

With the advent and development of micro-electromechanical systems (MEMS) technology in the 90’s new mechanical resonators were developed, such as the doubly clamped beam resonators that could reach very high quality factors. The dynamics of the fundamental mode of these resonators is well approximated by the Duffing equation. Furthermore, these micromechanical devices might exhibit, if properly tuned, a bistable response that can be quantitatively modeled by the bistability obtained in Duffing oscillators Aldridge and Cleland (2005). More recently, amplifiers that operate near the threshold of bifurcations were shown to present very high-gain amplification R. Almog, S. Zaitsev, O. Shtempluck, and Buks (2006, 2007); Almog et al. (2007). A degenerate parametrically excited Duffing amplifier was proposed by Rhoads et al. Rhoads and Shaw (2010).

One might be interested in a very sensitive high gain amplifier with a spread spectrum, so as to sample selectively a broad band of frequencies. This is particularly interesting for applications requiring spectrum sensing, and one way to achieve this is through the use of frequency combs Jamali and Babakhani (2019); Yang et al. (2019). The response of frequency combs to a narrowband small signal is given by a series of narrow peaks equally spaced in frequency. Cao et al. Cao et al. (2014) proposed a phononic frequency-comb generation based on an ac driven Fermi-Pasta-Ulam-Tsingou chain. In 2017, Ganesan et al Ganesan et al. (2017, 2018) created the first mechanical resonators that had spectral lines similar to optical frequency-combs. Their experimental apparatus is based on two symmetrical cantilevers mechanically coupled to one another at their bases. They used two coupled nonlinear normal modes to model the dynamics in which the first mode is resonantly excited by an added external drive, while the second mode is parametrically excited by the first mode. These experiments were followed by Czaplewski et al Czaplewski et al. (2018) in 2018, who used another two-mode coupling nonlinear model to describe experimental data from a mechanical resonator with both flexural and torsional vibrations. In Ref. Houri et al. (2019), the authors propose a theoretical model based on an approximation to a nonlocal Euler-Bernoulli model with many normal modes to describe frequency combs. More recently, Singh et al. Singh et al. (2020) observed a mechanical frequency comb in a graphene-silicon nitride hybrid resonator. They modeled the observed phenomena using two coupled normal modes, one nonlinear (graphene) and the other linear (SiN).

We present one parametrically-driven Duffing oscillator model and make a nonlinear analysis of it based on the averaging method. When there is no external added ac excitation, this nonlinear system can present a bistable region (in which the quiescent solution and one limit cycle are stable as seen in the original non-autonomous system). From the perspective of the averaged equations one has a tristable region, whose onset corresponds to a dual saddle-node bifurcation. There is also a threshold for parametric instability, which corresponds to a pitchfork bifurcation in the averaged equations, either supercritical or subcritical. With the application of an external ac drive, the nonlinear system may present a frequency-comb-like behavior as the parametric pump amplitude is increased past the dual saddle-node bifurcation. It occurs just after a sharp increase in the spectral component corresponding to parametric instability. This corresponds to a symmetry-breaking transition, in which the number of spectral peaks is doubled with peaks equally spaced in the frequency spectrum. All these responses are captured by the first-order averaging method.

In this article we show that a weaker form of mechanical frequency combs (MFCs) gradually appears when the parametric pump amplitude is increased. In addition to that, a stronger form of MFCs occurs after a symmetry-breaking bifurcation. In this form, there are twice as many peaks in the Duffing oscillator response to the added ac excitation as in the weaker form of the frequency comb. Furthermore, unlike the models proposed in Refs. Ganesan et al. (2017, 2018); Czaplewski et al. (2018), we show that only one parametrically-driven nonlinear mode is needed to present the frequency-comb-like behavior in the spectral response of the resonator. Different from the analysis developed by Bryant and Wiesenfeld Bryant and Wiesenfeld (1986), there is no typical period-doubling bifurcation here.

The one-degree of freedom model we use to describe the dynamics of a parametrically-driven nonlinear resonator is given by

$$\ddot{x}(t)=-\gamma\dot{x}(t)-x(t)-\alpha x^{3}(t)+F_{p}\cos(2\omega t)x(t)$$

(1)

in dimensionless units. Here $\gamma$ is the dissipation rate, $\alpha$ is the nonlinear coefficient, $F_{p}$ is the parametric pump amplitude, and $2\omega$ is the parametric pump angular frequency. Assuming $\gamma,\alpha,$ and $F_{p}$ are $O(\epsilon)$, with $0<\epsilon<<1$, we can apply the averaging method Guckenheimer and Holmes (1983) to obtain a slow autonomous dynamics. This is accomplished via the transform

	$\displaystyle x(t)$	$\displaystyle=u(t)\cos(\omega t)-v(t)\sin(\omega t),$		(2)
	$\displaystyle\dot{x}(t)$	$\displaystyle=-\omega\left[u(t)\sin(\omega t)+v(t)\cos(\omega t)\right].$

After applying this change of variables and neglecting fast oscillating terms, via Poincaré weakly non-linear transformation, we obtain

	$\displaystyle\dot{u}$	$\displaystyle=-\frac{1}{2\omega}\bigg{\{}\gamma\omega u+\left[\Omega+\frac{F_{p}}{2}+\frac{3\alpha}{4}(u^{2}+v^{2})\right]v\bigg{\}},$		(3)
	$\displaystyle\dot{v}$	$\displaystyle=-\frac{1}{2\omega}\left\{\left[-\Omega+\frac{F_{p}}{2}-\frac{3\alpha}{4}(u^{2}+v^{2})\right]u+\gamma\omega v\right\},$

where $\Omega=1-\omega^{2}=O(\epsilon)$. The fixed points are obtained from the solution of

		$\displaystyle\gamma\omega u+\left[\Omega+\frac{F_{p}}{2}+\frac{3\alpha}{4}r^{2}\right]v=0,$		(4)
		$\displaystyle\left[-\Omega+\frac{F_{p}}{2}-\frac{3\alpha}{4}r^{2}\right]u+\gamma\omega v=0,$

where $r^{2}=u^{2}+v^{2}$. The characteristic equation based on Eq. (4) can be written as

$$\frac{F_{p}^{2}}{4}=(\gamma\omega)^{2}+\left(\Omega+\frac{3\alpha}{4}r^{2}\right)^{2}\\$$

(5)

The steady-state squared amplitude is given by

$$r^{2}=-\frac{4}{3\alpha}\left[\Omega\pm\sqrt{\frac{F_{p}^{2}}{4}-\gamma^{2}\omega^{2}}\,\right].$$

(6)

In Fig. 1, we present the bifurcation diagram in the $\omega\times F_{p}$ parameter space of the averaged equations of motion, Eq. (3), of the parametrically-driven Duffing oscillator. When $\omega<1$ ($\Omega>0$) and $|F_{p}|<2\sqrt{(\gamma\omega)^{2}+\Omega^{2}}$ the only possible solution is $r=0$. When $\omega>1$ ($\Omega<0$) and $|F_{p}|<2\gamma\omega$, we also find the fixed-point solution with $r=0$. These conditions characterize region I of the bifurcation diagram. In region II, we have $|F_{p}|>2\sqrt{(\gamma\omega)^{2}+\Omega^{2}}$. In this region, the steady-state amplitude $r$ admits two solutions: $r=0$ and

$$r_{+}=\sqrt{\frac{4}{3\alpha}\left(-\Omega+\sqrt{\frac{F_{p}^{2}}{4}-\gamma^{2}\omega^{2}}\,\right)}.$$

(7)

With $u=r\cos\theta$ and $v=r\sin\theta$, we obtain the steady-state angle $\theta$ from

$$\tan\theta_{+}=-\dfrac{\gamma\omega}{\Omega+F_{p}/2+\frac{3\alpha r^{2}_{+}}{4}}.$$

In region III, we have $\Omega<0$ and $2\gamma\omega<|F_{p}|<2\sqrt{(\gamma\omega)^{2}+\Omega^{2}}$. In this region, there are three different solutions for the steady-state amplitude $r$: $0$, and

$$r_{\pm}=\sqrt{-\frac{4}{3\alpha}\left(\Omega\pm\sqrt{\frac{F_{p}^{2}}{4}-\gamma^{2}\omega^{2}}\,\right)}.$$

(8)

The corresponding values of the steady-state angle $\theta$ can be obtained from

$$\tan\theta_{\pm}=-\dfrac{\gamma\omega}{\Omega+F_{p}/2+\frac{3\alpha r^{2}_{\pm}}{4}}.$$

We determine the stability of the fixed points based on the eigenvalues of the Jacobian matrix. At a fixed point $(\bar{u},\bar{v})$, it is given by

	$\displaystyle Df$	$\displaystyle=\left(\begin{array}[]{cc}\frac{\partial f_{1}}{\partial u}&\frac{\partial f_{1}}{\partial v}\\ \frac{\partial f_{2}}{\partial u}&\frac{\partial f_{2}}{\partial v}\\ \end{array}\right)_{(\bar{u},\bar{v})}$		(9)
		$\displaystyle=-\frac{1}{2\omega}\left(\begin{array}[]{cc}\gamma\omega+\frac{3\alpha}{2}\bar{u}\bar{v}&\Omega+\frac{F_{p}}{2}+\frac{3\alpha}{4}(\bar{u}^{2}+3\bar{v}^{2})\\ -\Omega+\frac{F_{p}}{2}-\frac{3\alpha}{4}(3\bar{u}^{2}+\bar{v}^{2})&\gamma\omega-\frac{3\alpha}{2}\bar{u}\bar{v}\\ \end{array}\right).$

For the $(0,0)$ fixed point, we obtain

$$Df(0,0)=-\frac{1}{2\omega}\left(\begin{array}[]{cc}\gamma\omega&\Omega+\frac{F_{p}}{2}\\ -\Omega+\frac{F_{p}}{2}&\gamma\omega\\ \end{array}\right).$$

(10)

The eigenvalues are given by

$$\lambda_{\pm}=-\frac{\gamma}{2}\pm\frac{1}{2\omega}\sqrt{\frac{F_{p}^{2}}{4}-\Omega^{2}}\,$$

(11)

Note that when $|F_{p}|<2\Omega$, the eigenvalues are complex and the quiescent mode fixed-point $(0,0)$ is a stable spiral. When $2\Omega<|F_{p}|<2\sqrt{(\gamma\omega)^{2}+\Omega^{2}}$, the two eigenvalues are real and this fixed point is a stable node. Hence, in regions I and III, this fixed point is stable, whereas, in region II, it is unstable. At the transition line to parametric instability, the stable node becomes an saddle point as one goes from either region I or III into region II.

In Fig. 2 we plot the bifurcation diagrams along lines (a) and (b) of Fig. 1. In frame A, we obtain a supercritical pitchfork bifurcation as one crosses the parametric instability threshold along line (a). In frame B, we obtain a dual saddle-node bifurcation at about $F_{p}=0.2125$, and a subcritical pitchfork bifurcation at $F_{p}=0.334101$ along line (b).

If we add to Eq. (1) an ac excitation, we obtain a parametric amplifier. Here we call it a Duffing amplifier (DA). It is described by the equation

$\displaystyle\ddot{x}(t)$

$\displaystyle=$

$\displaystyle-\gamma\dot{x}(t)-x(t)-\alpha x^{3}(t)+F_{p}\cos(2\omega t)x(t)+F_{s}\cos(\omega_{s}t+\varphi_{0}),$

(12)

where $F_{s}$ is the amplitude, $\omega_{s}$ is the angular frequency, and $\varphi_{0}$ is an arbitrary phase of the external ac excitation. Assuming that $F_{s}=O(\epsilon)$, $\delta=\omega_{s}-\omega=O(\epsilon)$ and all the other coefficients are also small as in the previous section, we can apply the averaging method. After doing so, we find the slowly-varying dynamics

	$\displaystyle\dot{u}$	$\displaystyle=-\frac{1}{2\omega}\left\{\gamma\omega u+\left[\Omega+\frac{F_{p}}{2}+\frac{3\alpha}{4}(u^{2}+v^{2})\right]v\right\}+\frac{F_{s}}{2\omega}\sin\varphi(t),$		(13)
	$\displaystyle\dot{v}$	$\displaystyle=-\frac{1}{2\omega}\bigg{\{}\left[-\Omega+\frac{F_{p}}{2}-\frac{3\alpha}{4}(u^{2}+v^{2})\right]u+\gamma\omega v\bigg{\}}-\frac{F_{s}}{2\omega}\cos\varphi(t),$

where $\varphi(t)=\delta t+\varphi_{0}$.

We used the Odeint function of the Python’s scientific library package SciPy Virtanen et al. (2020) to integrate Eqs. (12) and (13). The integration time-step used in each numerical result was $h=\frac{T}{512}$, where $T=\frac{2\pi}{\omega}$. To avoid end discontinuities in the the time series we chose the total time of integration to be an integer multiple of $T_{p}=\frac{\omega T}{|\delta|}=\frac{2\pi}{|\delta|}$. For each time series in which we performed the Fourier transform (FT) (Figs. 3-5), we run the equations of motion for a time interval of $24T_{p}$, of which we discarded the first half as an equilibration time so that all transients die out. The fast Fourier transform routine used was SciPy’s fftpack. In Figs 6-10, we calculated the peak amplitudes using the following method Rajamani and Rajasekar (2017)

	$\displaystyle a_{n}$	$\displaystyle=\frac{2}{\Delta t}\int_{t_{0}}^{t_{0}+\Delta t}x(t)\cos((\omega+n\delta)t)dt,$		(14)
	$\displaystyle b_{n}$	$\displaystyle=\frac{2}{\Delta t}\int_{t_{0}}^{t_{0}+\Delta t}x(t)\sin((\omega+n\delta)t)dt,$

hence the amplitude of each peak is $r_{n}=\sqrt{a_{n}^{2}+b_{n}^{2}}$. This works because we have approximately $x(t)\approx\sum_{n=-N}^{N}\left[a_{n}\cos((\omega+n\delta)t)+b_{n}\sin((\omega+n\delta)t)\right]=\sum_{n=-N}^{N}r_{n}\cos((\omega+n\delta)t-\varphi_{n})$. Here $2N+1$ is the total number of visible peaks of the comb after the symmetry-breaking bifurcation, whereas before the symmetry-breaking bifurcation only peaks at odd values of $n$ contribute.

In frame (a) of Fig. 3, we show a time series of numerical integration of Eq. (12), which corresponds to parametric amplification in the DA with the parametric pump set at $F_{p}=0.21$ and the ac excitation amplitude set at $F_{s}=0.01$. The envelopes are obtained from the averaged slowly-varying dynamics given in Eqs. (13). The envelope of the pulses is approximately symmetric in time. In frame (b) of Fig. 3, we plot the FT corresponding to the time series shown in frame (a). The semi-analytic approximation to the numerical FT spectrum was obtained from the numerical integration of the corresponding averaged system, given in Eqs. (13), with the transformation defined in Eq. (2). The first-order harmonic balance approximation is roughly accurate, since it only predicts the signal ($\omega_{s}=\omega+\delta$) and idler ($\omega-\delta$) peaks in a parametric amplifier. Here, besides the signal and idler peaks, we barely see two other sidebands due to the nonlinearity of our system.

In frame (a) of Fig. 4, we show a time series of numerical integration of Eq. (12), which corresponds to parametric amplification in the DA with the parametric pump set at $F_{p}=0.21$ and the ac excitation amplitude set at $F_{s}=0.02$. The envelopes are obtained from the averaged slowly-varying dynamics given in Eqs. (13). Note the increased time-inversion asymmetry in the envelope. In frame (b) of Fig. 4, we plot the Fourier transform (FT) corresponding to the time series shown in frame (a). The semi-analytic approximation to the numerical FT spectrum was obtained from the numerical integration of the corresponding averaged system, given in Eqs. (13), with the transformation defined in Eq. (2). The first-order harmonic balance approximation becomes inaccurate, since it only predicts the signal and idler peaks in a parametric amplifier. Here, besides the signal and idler peaks, we also have sidebands spaced from one another by $2\delta$ and symmetrically positioned around these two central peaks. One can clearly see a frequency-comb-like spectrum with ten easily seen peaks.

In frame (a) of Fig. 5, with a higher value of the parametric pump amplitude ($F_{p}=0.24$), the envelope of the DA response becomes even more asymmetric. One can see that for the same time span the number of pulses is halved as compared to the time series of Fig. 4. This could indicate that there was a period-doubling bifurcation in the dynamical system of Eq. (13), when the pump amplitude was increased from $F_{p}=0.21$, in Fig. 4(a), to $F_{p}=0.24$. In reality there is no period-doubling bifurcation there. In frame (b) of Fig. 5, in addition to the signal and idler peaks, we also have a strong peak at $\omega$, which is at half the parametric drive frequency. Further new peaks can be seen at $\omega\pm 2\delta$, $\omega\pm 4\delta$, $\omega\pm 6\delta$, …. If one increases further the parametric drive amplitude, then one gets a full transition to parametric instability in which one gets a very strong peak at $\nu=\omega$ with two small sidebands, the signal and the idler. The other spectral peaks of the frequency-comb tend to decrease gradually.

In Fig. 6, we show that a parametric instability occurs in a narrow window of amplitudes of the added ac excitation. In this region occurs the strongest forms of frequency-comb response such as occurs in Fig. 5. For $F_{s}<0.01386$ the solutions are more symmetrical such as in Fig. 3, inside the window they are least symmetrical. Whereas for $F_{s}>0.0181$, there is parametric instability suppression, less asymmetry, and the frequency-comb response is weaker, with peaks spaced by $2\delta$ (see Fig. 4).

In Figs. 7-10, we show how the seven main spectral peaks vary as a function of the parametric pump amplitude. In Figs. 7-8, we fix $\omega$ at $0.95$ and we increase $F_{p}$ along line $(a)$ of the bifurcation diagram of Fig. 1. We can see parametric instability suppression as evidenced by the delayed increase in the amplitude of the spectral peak at $\omega$. The sharp increase in this peak only occurs well past the parametric instability transition of Fig. 1, when there is no external drive ($F_{s}=0$). One sees that the higher $F_{s}$ is the more suppression in parametric instability there is.

In Figs. 9-10, $\omega=1.0625$ and we vary $F_{p}$ along line $(b)$ of the bifurcation diagram of Fig. 1. One can see that the sharp increase in the spectral peak at $\omega$ of the DA response generates an even stronger frequency-comb-like behavior, with twice as many peaks, now spaced out in frequency from one another by $\delta$. By comparing the results of Figs. 9 and 10, one sees that with increasing value of $F_{s}$ there is also a suppression of parametric instability. In this case, the parametric instability transition still occurs in the bistable region, below the threshold to instability in the parametric oscillator (dashed line). Apparently, counterintuitively, the larger $F_{s}$ (the larger the external excitation) is, the more stable the smaller amplitude of the parametrically-driven Duffing oscillator response is.

Here we investigated a one-degree-of-freedom parametrically-driven Duffing oscillator with a small added ac drive that could present suppression of parametric instability or present spectral peaks similar to recent experimental results of mechanical frequency combs Ganesan et al. (2017, 2018); Czaplewski et al. (2018); Singh et al. (2020). We have seen two types of frequency comb behavior: one weaker and one stronger. The stronger frequency comb has twice as many peaks as the weaker form with peaks at half the distance from one another. We claim that the fundamental cause of the frequency-comb dynamical behavior is a symmetry breaking bifurcation that occurs near the parametric instability transition. We have also shown that the averaging method can capture this spectral response. The stronger form of the frequency comb arises after a subcritical pitchfork bifurcation of the averaged system of equations of the DA. In addition, the simple model we propose here could be used as a theoretical framework, or a toy model, for the study of the frequency comb phenomenon in mechanical oscillators.

Furthermore, we point out that our theory is similar to the one developed by K. Wiesenfeld and McNamara in the 80’s for amplification of small signals near bifurcation points. Their theory presented in Refs. Wiesenfeld and McNamara (1985, 1986) is a linear response theory based on Floquet theory, whereas here we present a nonlinear response theory based on the averaging method. Also, we construct an approximate analytical solution, whereas their model is generic. One would still have the difficult task of obtaining the Floquet eigenfunctions. In addition, the perturbing terms in their prototype example is a parametric drive, whereas in our models the perturbing terms are the added ac signals.

It is worth mentioning that Bryant and Wiesenfeld Bryant and Wiesenfeld (1986) (see their Fig. 12) obtained an effect similar to the frequency-comb spectral peaks seen here. There are several important differences between our physical systems. Their Duffing oscillator is not driven parametrically and has a very low quality factor, which effectively reduces its dimensionality. They present a suppression of a period-doubling bifurcation, whereas here we have a suppression of parametric instability due to the small added external ac signal. We saw that the suppression increases with the ac signal amplitude. In our system, the frequency-comb behavior occurs more strongly in the region III of the bifurcation diagram, where quiescent solution and limit cycle are both stable. Before the jump in amplitude of the central peak, the peaks are spaced by $2\delta$, after the jump the peaks are spaced by $\delta$. The frequency-comb-like spectrum is stronger in a narrow region after the jump of the central peak.

We investigate one simple nonlinear dynamical system to support our claim. The proposed system does not present mode coupling and, thus, has lower dimensionality than the phenomenological models used to explain the spectral signatures of the experimental mechanical frequency combs. We believe this could lead to further research to find simpler apparatuses, with less parameters and less normal modes, that could deliver a frequency-comb-like behavior. We also believe this work could spur more theoretical work to help understand better the connection between the frequency comb behavior and the various bifurcations that occur in conjunction. For example, in Ganesan et al. Ganesan et al. (2018), one can see clearly in their Fig. 2 a spectral map that we think has the footprints of a symmetry-breaking bifurcation as we have seen here. Just below half-way their spectral map, as the pump amplitude is increased, twice as many peaks appear, with each new peak in the middle of two pre-existent ones. Finally, we point out that the equations of the model proposed by Singh et al. Singh et al. (2020) for describing their mechanical frequency comb can be recast in a form similar to Eq. (13) if one integrates out the degrees of freedom of the second normal mode.