text:
# Bifurcation structure of two Coupled Periodically driven double-well Duffing Oscillators

## I Introduction

In recent years, a large number of theoretical calculations, numerical simulations and experiments has been carried out on systems of coupled anharmonic oscillators which provide fundamental models of the dynamical problems in several disciplines. Such models describe a wide range of processes, helping in general to understand the routes to chaos that take place in biological, chemical, physical, electromechanical and electronical systems. Here we refer only to recent papers close to the subject of our own\[1-18\]. Coupling may arise between oscillators of the same types or oscillators of different types. Amongst these systems, the most intensively investigated examples are the Duffing oscillators\[1-8\] (strictly dissipative systems which converge to quiescence when not driven), the Van der Pol oscillators\[9-11\] (self-excited systems which final behavior when not driven is a limit cycle) and the coupled Van der Pol - Duffing oscillators \[12-13\].
In this paper we present the results of an investigation of two identical coupled double-well Duffing oscillators subjected to a periodically driven force. Our study was stimulated by the earlier work of Kozlowski et al.  who have studied the model of the same type, but considered only single-well potentials. Using bifurcation diagrams and phase diagrams these authors have demonstrated that the global pattern of bifurcation curves in parameter space consists of repeated subpatterns similar to the superstructure found in a one-dimensional driven Duffing oscillator. They have also proven the existence of a Hopf Bifurcation which does not appear in a model of one-dimensional Duffing oscillator. Other interesting results have been reported in the context of coupled Duffing Oscillators. In ref.  Kunick et al showed that two coupled oscillators can behave regularly with nonzero coupling and chaotically with zero coupling. In a model of two coupled Duffing oscillators driven with incommensurate frequencies and coupled additively, Stagliano et al  observed the period-doubling of an attracting $`T^2`$ torus and its destruction in parameter space. Migration control in two Duffing oscillators by open-plus-closed-loop control method and adaptive control algorithm were studied by Paul et al . These authors, in another system of two coupled Duffing oscillators, have also investigated a basin of attraction with several coexisting chaotic attractors and synchronization of chaos . Yagasaki , after having described an Auto driver Hommap for numerical analysis of homoclinic in maps of periodically forced Duffing systems, modified a version of the homoclinic Melnikov method for orbits homoclinic to two types of periodic orbits. These theories have been successfully applied to a weakly coupled Duffing oscillators. Moreover Yin et al  investigated the effect of phase difference in the driving forces of two coupled identical Duffing oscillators. They observed desynchronization and lag synchronization of chaotic attractors. Bifurcation behaviors, showing a kind of Hopf bifurcation, from synchronous chaos of a chain of Duffing oscillators have been explored by Ma Wen-Qi et al.  using the generalized winding number in tangent space.
The equations of motion for the two identical coupled double-well Duffing oscillators that we are interested in are the following:

$`{\displaystyle \frac{d^2x}{dt^2}}`$ $`=`$ $`\alpha {\displaystyle \frac{dx}{dt}}+xx^3+k(yx)+f\mathrm{cos}(\mathrm{\Omega }t)`$ (1)
$`{\displaystyle \frac{d^2y}{dt^2}}`$ $`=`$ $`\alpha {\displaystyle \frac{dy}{dt}}+yy^3k(yx)`$ (2)
where $`\alpha `$ is the damping parameter, $`k`$ the coupling parameter, $`f`$ and $`\mathrm{\Omega }`$ the amplitude and the frequency of the driving force, respectively. In system (1), we represent the oscillators by the state variables $`(x,dx/dt)`$ as oscillator $`A`$ or subsystem $`A`$ subjected to the periodic force $`f(t)=f\mathrm{cos}(\mathrm{\Omega }t)`$, which is coupled to the other one with the state variables $`(y,dy/dt)`$ as oscillator $`B`$ or subsystem $`B`$. Here $`f`$ and $`\mathrm{\Omega }`$ are the control parameters of our system. As evident, for $`k=0`$, equation (1) describes two uncoupled Duffing Oscillators with their motion governed by $`\alpha `$ and $`f`$. The coupling considered can be interpreted as a perturbation of each oscillator through a signal proportional to the difference of their amplitudes. A natural question to ask is how chaotic attractors arise as these parameters are varied. In the case of one dimensional double-well Duffing oscillator, quasiperiodic and period-doubling routes to chaos have been found \[19-27\].
The rest of the paper is structured as follows. In section II, we establish a linear stability analysis of the system projected on the Poincaré map. Section III is devoted to the survey of bifurcation diagrams and the chaotic behavior is characterized in section IV. We end with the conclusions of the work in section V.

## II Stability analysis

For stability analysis and even for purposes of numerical simulations, it is convenient to transform the second order differential equation (1) into an autonomous system of first order differential equations of the following form:

$`{\displaystyle \frac{dx}{dt}}`$ $`=`$ $`v`$ (3)
$`{\displaystyle \frac{dv}{dt}}`$ $`=`$ $`\alpha v+xx^3+k(yx)+f\mathrm{cos}(\theta )`$ (4)
$`{\displaystyle \frac{dy}{dt}}`$ $`=`$ $`w`$ (5)
$`{\displaystyle \frac{dw}{dt}}`$ $`=`$ $`\alpha w+yy^3k(yx)`$ (6)
$`{\displaystyle \frac{d\theta }{dt}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }}{2\pi }}`$ (7)
or equivalently

$`{\displaystyle \frac{dV}{dt}}=F(V,\mu )`$ (8)
where $`T=2\pi /\mathrm{\Omega }`$ is the period of the driving force, $`\theta `$ the cyclic variable, $`V(x,v,y,w,\theta )`$ an autonomous vector field and $`\mu (\alpha ,k,f,\mathrm{\Omega })`$ an element of the parameter space. This system generates a flow $`\varphi =\{\varphi ^T\}`$ on the phase space $`𝐑^4\times S^1`$and there exists a global map:

$`P:\mathrm{\Sigma }_c\mathrm{\Sigma }_c`$ (9)
$`V_p(x,v,y,w)P(V_p)`$ $`=`$ $`\{\varphi ^T\}|_{\mathrm{\Sigma }_c(x,v,y,w,\theta _0)}`$ (10)
with $`\theta _0`$ being a constant determining the location of the Poincaré cross-section and $`(x,v,y,w)`$ the coordinates of the attractors in the Poincaré cross-section $`_c`$ defined by: $`_c=\{(x,v,y,w,\theta )𝐑^4\times S^1|\theta =\theta _0\}`$. System (3) is symmetric since the transformation $`S:(x,v,y,w,\theta )(x,v,y,w,\theta +\pi )`$ leaves (3) invariant. This system can support symmetric orbits and also asymmetric ones which are not invariant with the transformation S. We use a perturbation analysis to proceed analysis of solution stability in the poincaré map $`\mathrm{\Sigma }_c`$. This method consists of differentiation and translation. Thus if we add a small perturbation $`\delta V_p`$($`\delta x`$, $`\delta v`$, $`\delta y`$, $`\delta w`$) to the steady state $`V_{p0}`$($`x_0,v_0,y_0,w_0`$) we obtain equations for perturbed motion

$$x=x_0+\delta x,y=y_0+\delta y,v=v_0+\delta v,w=w_0+\delta w$$
(11)
Substituting equations (11) into equations (3) and developing into powers of the above small perturbations yields the following variational differential equation:

$$\frac{d\delta V_p}{dt}=DF(V_{p0})\delta V_p,$$
(12)
with

$`DF(V_{p0})=\left(\begin{array}{cccc}0& 1& 0& 0\\ 1k3x_0^2& \alpha & k& 0\\ 0& 0& 0& 1\\ k& 0& 1k3y_0^2& \alpha \end{array}\right)`$

where $`DF(V_{p0})`$ is the Jacobian $`4\times 4`$ matrix describing the vector field along the solution $`\delta V_p(t)`$ and $`V_{p0}`$ being an equilibrium point or the steady state. Next for the very small amplitudes, we neglect in the matrix $`DF(V_{p0})`$ the powers higher than one, that is, $`x_0^2`$ and $`y_0^2`$. The solution after one period $`T`$ of the oscillations in the linearised Poincaré map is simply expressed as:

$$\delta V_p(T)=\delta V_{p0}\mathrm{exp}(DF(V_{p0})T),\delta V_p(0)=\delta V_{p0}$$
(13)
where $`DF(V_{p0})`$ is the monodromy matrix of a periodic orbit connecting arbitrary infinitesimal variations in the initial conditions $`\delta V_{p0}`$ with corresponding changes $`\delta V_p(T)`$, after one period T. The stability of the periodic motion is therefore determined according to the real parts of the roots of the characteristic equation $`det(DF(V_{p0})I\lambda )=0`$ written in the following form:

$$\lambda ^4+A_3\lambda ^3+A_2\lambda ^2+A_1\lambda +A_0=0$$
(14)
where $`(A_i)`$,$`i=0,1,2,3`$ are the coefficients depending on the two parameters $`k`$, $`\alpha `$ and $`\lambda `$ =$`(\lambda _i)`$ the eigenvalues of $`DF(V_{p0})`$. The roots of equation (14) are obtained using the Bairstow-Newton-Raphson algorithm  and for different values of $`k[5,5]`$ and $`\alpha [0.2,0.2]`$. The results as the spectrum of the eigenvalues are depicted in a complex plane $`𝐂`$ (see fig.1). This figure enables us to get an idea of different types of bifurcation likely to appear in the system. Since $`DF(V_{p0})`$ is a real matrix, complex eigenvalues occur in a complex conjugate pairs responsible of the symmetry observed along the real axis. Thus if $`\lambda _i`$ is real, it is clear from equation (13) that the eigenvalues is the rate of contraction ($`\lambda _i<0`$) or expansion ($`\lambda _i>0`$) near the steady state. Next if $`\lambda _i`$ is complex, the real part of $`\lambda _i`$ gives the rate of contraction $`(Re(\lambda _i)<0)`$ or expansion $`(Re(\lambda _i)>0)`$ of the spiral while the imaginary part $`Im(\lambda _i)`$ is of the frequency rotation. This can be well understood with the following expression of the eigenvalues of the linearized Poincaré map:

$$\sigma _i=exp(Re(\lambda _i))[cos(Im(\lambda _i)T)+jsin(Im(\lambda _i)T)]$$
(15)
Globally, if $`Re(\lambda _i)<0`$ for all $`\lambda _i`$, then all sufficiently small perturbations tend toward zero as t goes to infinity and the steady state (nodes (n), saddle nodes (sn), spiral (sp)) is stable. If $`Re(\lambda _i)>0`$ for all $`\lambda _i`$, then any small perturbation grows with the time and the steady state $`(n,sn,sp)`$ is unstable (A stable or unstable equilibrium state with no complex eigenvalues is often called node). On the other hand if there exists $`i`$ and $`l`$ such that $`Re(\lambda _i)<0`$ and $`Re(\lambda _l)>0`$, therefore the equilibrium state is non-stable (A non-stable equilibrium state is often called saddle; an equilibrium point whose eigenvalues all have a non zero real part are called hyperbolic ). Having regards to the precedent with again a look on the spectrum of fig.1, it follows that our system of the given parameter k and $`\alpha `$ can undergo many types of bifurcations namely: saddle-node $`(Re(\lambda _i)=1)`$, period-doubling $`(\lambda _i=1)`$, Hopf bifurcations ($`\lambda _i=\gamma \pm j\beta ,(j^2=1`$) with $`\gamma <0`$ and the motion has the form of exponential increasing beats with period $`2\pi /\beta `$ ) and of course symmetry-breaking bifurcation which is often a prerequisite for the first period-doubling bifurcation . Except Hopf bifurcation, such bifurcations have been successfully found in one-dimensional double-well Duffing oscillators subjected to a periodically driven force  using the eigenvalues called Floquet multipliers of the linearized Poincaré map. In a coupled identical single-well Duffing oscillators, Kozlowski et al  showed by linear analysis a similar eigenvalues scenario and they conjectured that the resulting alternating bifurcation sequences have to be expected for all systems of coupled strictly dissipative oscillators that are driven periodically. The above analysis allows us to know how the steady state becomes locally unstable and to be aware of the type of bifurcation expected in the system. Here, no exploding amplitude is possible since the shape of the potential offers globally bounded solutions $`(V(x,y)=(1/2)x^2(1/2)y^2)+(1/4)x^4+(1/4)y^4)+\mathrm{}as|x|,|y|\mathrm{})`$.

## III Local Bifurcations

The aim of this section is to seek numerically the routes which lead to chaotic solutions of the system when the control parameter $`\mathrm{\Omega }`$ evolves, for different values of the driving force $`f`$. Bifurcation diagrams are very good for numerical as well as for experimental studies when there is a tunable parameter. When a control parameter is varied and bifurcation takes place a qualitative change of the system happens. For the numerical computations of these diagrams the control parameter $`\mathrm{\Omega }`$ is increased from an initial value $`\mathrm{\Omega }_i`$ to a final value $`\mathrm{\Omega }_f`$ and then decreased from $`\mathrm{\Omega }_f`$ to $`\mathrm{\Omega }_i`$ in a very small step ($`\mathrm{\Delta }\mathrm{\Omega }=10^5`$). The last computed cyclic point for a given value of $`\mathrm{\Omega }`$ is always used as a new initial value for the next value of $`\mathrm{\Omega }`$. Starting with initial conditions $`(x,v)=(1,0),(y,w)=(1,0)`$ at $`\mathrm{\Omega }_i`$, system (3) is integrated, using the standard Runge Kutta Algorithm , for $`100`$ periods of the driving force until the transient has died out; the trajectories expected are attractors and the local calculation error is sufficiently small. Then to find out whether the trajectory is periodic (quasiperiodic) or chaotic, the system is integrated for the next $`180`$ periods in order to catch the maximum of coexisting attractors. This procedure allows us to make sure that coexisting attractors are realised due to different initial conditions which are handed throughout from one parameter to another. We chose $`k=5`$ and $`\alpha =0.1`$ for all examples in the subsequent section. The bifurcation diagrams obtained here show the projection of the attractors in the Poincaré section onto $`x`$ or $`v`$ versus the control parameter $`\mathrm{\Omega }`$ and for different values of $`f`$. In a one-dimensional double-well Duffing oscillator \[19-21,23-27\], it is shown that chaotic solutions arise through quasiperiodic, period-doubling cascade with a sequence of symmetry-breaking and saddle-nodes. Hopf bifurcations have not been found, but are abundant in our model.
We computed a great number of diagrams, and some of them showing the main features of the system are presented in figures 2-4. It is important to note that no hysteresis has been detected while checking all these diagrams. We can therefore avoid any confusion over the periods of attractors when reading this type of diagram; for instance one single period-2 yields to 2 points for one frequency value whereas two coexisting attractors each of period-1 also yield to two points for one frequency. These diagrams show a great number of coexisting attractors (chaotic domains) intermingled with imbricated windows made up of periodic solutions of different periodicity, period-doubling of both types, sudden Chaos and Hopf bifurcation. The so called chaotic domains (or chaotic sea) will be characterized in the following section.
At lowest frequencies ($`0<\mathrm{\Omega }0.4`$) quasiperiodic as well as chaotic solutions occur together with a number of symmetry-breaking (sb), period doubling (pd) and saddle-nodes (sn). In fig.2(a) where $`f=15`$, it appears in the neighborhood of $`\mathrm{\Omega }=0.385`$, 6 pairs of sb followed by 3 reversed pd till the chaotic domain. When the driving force is increased $`f=20`$, overlapping pd cascade and sn appear as can be seen in fig.2(b). Similar features as in the latter case are also observed for small values of $`f`$ ($`f<15`$). In other ranges of the driving frequency bigger than $`\mathrm{\Omega }=0.4`$ the aformentioned features still exist and additional ones such as resonances and Hopf bifurcations appear in the system; the results are displayed in figures 3 and 4. In figure 3 it can be seen throughout overlapping pd cascade of both types, sb, resonances (R) and even sudden chaos (SC). When $`f=4.5`$ in fig.3(f) at $`\mathrm{\Omega }0.785`$, period-8 attractor undergoes reversed period-doubling cascade and period-2 is created as from $`\mathrm{\Omega }0.805`$. Besides, for $`f7`$ appear resonances (R) in fig.3(a,b,d,e), fig.6(a) and sudden chaos in fig.3(c). However in fig.4 when $`f`$ becomes much larger and $`\mathrm{\Omega }0.5`$, Hopf Bifurcation (H) is abundant and often appear more than once at a given frequency. This notation $`nH`$ means that there are $`n`$, ($`n3`$), integer Hopf bifurcation at a given frequency. One can therefore observe -H- in (a), -3H- in (b), -2H-..-2H- in (c,d), -2H-..-3H-..-2H- in (e,f) and -3H-3H- in (g,h). In fig.4(h) for instance, at $`\mathrm{\Omega }=0.5925`$, a period-6 attractor undergoes reversed period-doubling to become a period-3; next after 3 pairs of symmetry-breaking near $`\mathrm{\Omega }0.6`$, appears again period-3 attractor as from $`\mathrm{\Omega }0.601`$. And then at $`\mathrm{\Omega }=0.603`$ this attractor of period-3 undergoes 3 Hopf bifurcations. Windows of periodic solutions separated by quasiperiodic ones are clearly visible. The Poincaré cross-section of quasiperiodic attractors of this kind (-3H-) consists of three invariant tori (destroyed or not). Thus at any $`\mathrm{\Omega }`$ where $`nH`$ exist, the projection of the attractor corresponds to $`n`$ Tori on the Poincaré cross-section. The torus attractor arises from Hopf bifurcation which is a bifurcation from a fixed point to an invariant curve. At the transition to quasiperiodic motion, this fixed point (projection of the limit cycle in to the cross-section) loses its stability and gives birth to an invariant circle, which is a cross-section of the torus in the flow. The three pairs of sb appearing at the neigborhood of $`\mathrm{\Omega }0.6`$ (fig.4(g,h)) seem to be apparently a window of period-6 attractor. In fact it is just a period-3 attractor as can be illustrated further in fig.8(b).This orbit, born asymmetrically, coexists with its inversion symmetric orbit of the same period-3 and are caught simultaneously.
Next we pursue our investigation for large values of $`\mathrm{\Omega }`$ ($`\mathrm{\Omega }3`$) and for different values of $`f`$. It turns out that the features observed so far have disappeared and the system becomes regular providing attractors with a small number of periods as can be seen in figure 5. The phase diagrams which often show bifurcations, curves, surfaces and which is out of the scope of this paper could also help to match the features obtained in this analysis.

## IV Chaotic behavior

To better support the results obtained above, we intend in this section to characterize chaotic behavior. For this purpose some indicators are used namely the Lyapunov exponents and the power spectrum density. Also, Poincaré cross-sections and phase portraits are shown to illustrate some attractors present in the system. The Lyapunov exponent is one of the most important tools for understanding chaotic behavior, obtained by examining “a very sensitive dependence of initial conditions”. In particular it is generally well established that a rigorous measure of chaos may be given in terms of Lyapunov spectrum of the dynamical system. A positive Lyapunov exponent is characteric of chaos while zero and negative values of the exponent signify a marginally stable or quasiperiodic orbit and periodic orbit, respectively. Solving numerically the variational equation (12) together with system (3) we calculate the maximum Lyapunov exponent $`\lambda _{max}`$ in the Poincare cross-section with $`\lambda _{max}=lim_{\tau +\mathrm{}}\frac{1}{\tau }\mathrm{ln}(L(\tau ))`$, where $`L(\tau )=\left(\delta x^2+\delta v^2+\delta y^2+\delta w^2\right)^{1/2}`$. The spectrum is displayed in fig.6 together with a bifurcation diagram. Windows of chaotic domains and periodic solutions are clearly justified. Power Spectrum density using Fast Fourier Transformation (FFT) has also been computed with some values from fig.6 for a periodic solution ($`\mathrm{\Omega }=0.495`$, $`f=7`$) and a chaotic solution ($`\mathrm{\Omega }=0.47`$,$`f=7`$) (see fig.7). Fig.8 shows the phase portraits and Poincaré cross-section. In fig.8(a,b) plotted are the phase portraits, related to fig.4(g,h), showing a period-3 attractor for $`f=30`$, $`\mathrm{\Omega }=0.6025`$(a) and a coexistence of two asymmetric attractors, each of period-3 for $`f=30`$, $`\mathrm{\Omega }=0.6`$ (b). The dots ($``$) represent the corresponding points in the Poincaré cross-sections. In fig.8(c,d) two chaotic attractors are plotted in the Poincaré cross-section for $`\mathrm{\Omega }=0.5`$, $`f=2`$ and for $`\mathrm{\Omega }=0.375`$, $`f=15`$, respectively. Furthermore, we observe in fig.8(e,f) three $`T^2`$ tori doubling at $`\mathrm{\Omega }=0.605`$, $`f=30`$ and at $`\mathrm{\Omega }=1.14`$, $`f=15`$, two $`T^2`$ tori doubling (destroyed) corresponding for -3H- and -2H-, respectively (values from fig.4(c,h)). In the torus route to chaos, the original torus appears to split into two circles at the torus doubling bifurcation point. This route is reminiscent to the period-doubling route to chaos, although there are finite number of torus doublings before the onset of chaos motion. Also, the torus doubling route to chaos is a higher-dimensional phenomenon, requiring at least a four-dimensional flow, or a three-dimensional map. It is not observed in one-dimensional maps, unlike the period-doubling route to chaos .

## V Conclusion

To sum up we have investigated coupled double-well Duffing oscillators subjected to a periodic driving force. Our main conclusion is that the shape of the potential in the coupled Duffing Oscillators system has a profound influence on the routes to chaos taken by the system. By a linear analysis we have examined the stability of the steady state of the system leading to different types of bifurcation likely to appear for $`k[5,5]`$ and $`\alpha [0.2,0.2]`$. With the help of bifurcation diagrams, it has been shown that the bifurcation structure, which is complicated, depends strongly on the values of the control parameter $`\mathrm{\Omega }`$. For a simple electronic realization of our model, it is obvious that the oscillators may be very difficult to control for small values of $`\mathrm{\Omega }`$ due to the large windows for chaos. This system becomes regular for large values of $`\mathrm{\Omega }`$ ($`\mathrm{\Omega }3`$). Apart from the routes already encountered in a one-dimensional double-well Duffing oscillator such as quasiperiodic and period-doubling cascade, sudden chaos and mostly Hopf bifurcations have been newly found in this work. Besides, resonances (R), imbricated sb, sn, pd intermingled by chaotic domains are observed, rending the structure highly chaotic. Two and three $`T^2`$ tori doubling, have been also observed. To characterize chaotic behavior of this system, Lyapunov exponents and the power spectrum density using FFT have been employed. Poincaré cross-sections and phase portraits have been also used for better visualisation of periodic and chaotic attractors. An analytical approach using approximative technics could be an interesting topic for future work, and would probably shed further light on the rich behaviour of this model.

## VI ACKNOWLEDGMENTS

It is a pleasure to thank Prof. Dr. Jan M. Rost for illuminating discussions and Dr. Nic Shannon for a careful reading of the manuscript. Financial support by Alexander von Humboldt (AvH) Foundation/Bonn-Germany, under the grant of Research fellowship No.IV.4-KAM 1068533 STP is gratefully acknowledged.

## VII Figure caption

* Figure 1: Eigenvalues spectrum for the stability of the system with $`\alpha [0.2,0.2]`$ and $`k[5,5]`$. Note that when $`\alpha <0`$, the eigenvalues are in the half-plane left and in the half-plane right otherwise.
* Figure 2: Bifurcation diagrams showing the first coordinate $`x`$ for the driving amplitude $`f=15`$ $`(a)`$ and the second coordinate $`v`$ for $`f=20`$ $`(b)`$ of the Poincaré section versus the driving frequency $`\mathrm{\Omega }`$. Windows of periodic solutions and chaotic domains are visible. Near $`\mathrm{\Omega }=0.385`$ (a), $`6`$ pairs of $`sb`$ occur followed by the reversed $`pd`$ till the chaotic domain. Imbricated $`pd`$, $`sb`$, $`sn`$ are mentioned. Note that similar behavior as (b) is observed when $`\mathrm{\Omega }0.4`$ and for small values of $`f`$ ($`f<15`$).
* Figure 3: Bifurcation diagrams for the driving force $`f=4.5`$ (f), $`f=15`$ (a,b,e), $`f=20`$ (d), $`f=25`$ (c) with $`\mathrm{\Omega }>0.4`$. Throughout imbricated $`sb`$, $`sn`$, $`pd`$ of both types and chaotic domains are observed. Besides, resonances (R) appear at $`\mathrm{\Omega }0.475`$ (a,b), $`\mathrm{\Omega }0.445`$ (d), $`\mathrm{\Omega }0.555`$ (e)and sudden chaos at $`\mathrm{\Omega }0.408`$ (c).
* Figure 4: Bifurcation diagrams for $`f=30`$ (a,g,h), $`f=25`$ (b), $`f=15`$ (c,d), $`f=20`$ (e,f) showing many Hopf bifurcations (H). One can observe -H-(a), -3H-(b), -2H-..-2H (c,d),-2H-..-3H-..-2H-(e,f), -3H-3H- (g,h) (see also text)
* Figure 5: Bifurcation diagram for $`f=5`$ showing attractors of small number of periods. Note that the system becomes regular and this behavior is observed for $`\mathrm{\Omega }3`$ and for any value of $`f`$.
* Figure 6: Bifurcation diagram for $`f=15`$ (a) and the corresponding Lyapunov spectrum (b) in the Poincaré map. windows of chaotic domains are clearly justified with positive values of $`\lambda _{max}`$.
* Figure 7: Power spectra densities (PSD). The values of $`f=7`$ and $`\mathrm{\Omega }`$ are from the bifurcation diagrams shown in fig.6 with $`\mathrm{\Omega }=0.495`$ (a) periodic motion and $`\mathrm{\Omega }=0.47`$ (b) chaotic motion.
* Figure 8 : Phase portraits and projection of the attractors in the Poincaré section onto the two first coordinates $`x,v`$ of the Poincaré cross-section. At $`\mathrm{\Omega }=0.6025`$ (a), one period-3 orbit and at $`\mathrm{\Omega }=0.6`$ (b) coexistence of two asymmetric orbits of the same period-3. The values of $`\mathrm{\Omega }`$ are from fig.4 (g,h) and $`f=30`$. The dots ($``$) represent the corresponding points in the Poincaré cross-section. Two chaotic attractors shown in (c,d) for $`\mathrm{\Omega }=0.5`$, $`f=2`$ and $`\mathrm{\Omega }=0.375`$, $`f=15`$, respectively. Three $`T^2`$ tori doubling for $`\mathrm{\Omega }=0.605`$, $`f=30`$ (e) and two $`T^2`$ tori doubling (destroyed) for $`\mathrm{\Omega }=1.14`$, $`f=15`$ (f). The values of $`\mathrm{\Omega }`$ are from fig.4(c,h).


format: text

text:
# Bifurcations without Parameters:Some ODE and PDE Examples

## 1 . Applied motivation

In this article we sketch and illustrate some elements of the nonlinear dynamics near equilibrium manifolds. Denoting the equilibrium manifold by $`x=0`$, in local coordinates $`(x,y)\mathrm{𝖨𝖱}^n\times \mathrm{𝖨𝖱}^k`$, we consider systems

$$\begin{array}{ccc}\hfill \dot{x}& =& f(x,y)\hfill \\ \hfill \dot{y}& =& g(x,y)\hfill \end{array}$$
(1.1)
and assume

$$\begin{array}{ccccc}\hfill f(0,y)& =& g(0,y)\hfill & =& 0,\hfill \end{array}$$
(1.2)
for all $`y`$. For simplicity we will only address the cases $`k=1`$ of lines of equilibria, and $`k=2`$ of equilibrium planes. Sufficient smoothness of $`f,g`$ is assumed. The occurrence of equilibrium manifolds is infinitely degenerate, of course, in the space of all vector fields $`(f,g)`$ – quite like many mathematical structures are: equivariance under symmetry groups, conservation laws, integrability, symplectic structures, and many others. The special case

$$\begin{array}{ccc}\hfill g& & 0\hfill \end{array}$$
(1.3)
in fact amounts to standard bifurcation theory, in the presence of a trivial solution $`x=0`$; see for example . Note that condition (1.3), which turns the $`k`$-dimensional variable $`y`$ into a preserved constant parameter, is infinitely degenerate even in our present setting (1.2) of equilibrium manifolds. Due to the analogies of our results and methods with bifurcation theory, we call our emerging theory bifurcation without parameters. This terminology emphasizes the intricate dynamics which arises when normal hyperbolicity of the equilibrium manifold fails; see the sections below.

To motivate assumption (1.2), we present several examples. First, consider an “octahedral” graph $`\mathrm{\Gamma }`$ of $`2(m+1)`$ vertices $`\{\pm 1,\mathrm{},\pm (m+1)\}`$. The graph $`\mathrm{\Gamma }`$ results from the complete graph by eliminating the “diagonal” edges, which join the antipodal vertices $`\pm j`$, for $`j=1,\mathrm{},m+1`$. For $`m=1`$ we obtain the square, for $`m=2`$ the octahedron, and so on. Consider the system

$$\begin{array}{ccc}\hfill \dot{u}_j& =& f_j(u_j,\underset{k\pm j}{}u_k)\hfill \end{array}$$
(1.4)
of oscillators $`u_j\mathrm{𝖨𝖱}^n^{}`$ on $`\mathrm{\Gamma }`$, additively coupled along the edges by $`f_j`$. We assume an antipodal oddness symmetry of the individual oscillator dynamics

$$\begin{array}{ccc}\hfill f_j(u_j,0)& =& f_j(u_j,0).\hfill \end{array}$$
(1.5)
As a consequence, the antipode space

$$\begin{array}{ccc}\hfill \mathrm{\Sigma }& :=& \{u=(u_j)_{j\mathrm{\Gamma }};u_j=u_j\}\hfill \end{array}$$
(1.6)
is invariant under the flow (1.4). Moreover, the flow on $`\mathrm{\Sigma }`$ completely decouples into a direct product flow of the $`m+1`$ diagonally antipodal, decoupled pairs

$$\begin{array}{ccc}\hfill \dot{u}_{\pm j}& =& f_{\pm j}(u_{\pm j},0).\hfill \end{array}$$
(1.7)
For the square case $`m=1`$, this decoupling phenomenon was first observed in . For more examples see also .

An $`m`$-plane of equilibria arises from periodic solutions of the decoupled system (1.7). Assume (1.7) possesses time periodic orbits $`u_j(t+\phi _j)`$ of equal period $`T_j=2\pi `$, for $`j=1,\mathrm{},m+1`$. Choose arbitrary phases $`\phi _jS^1`$ and let $`u_j(t):=u_j(t),\phi _j=\phi _j`$. Then

$$\begin{array}{ccccc}\hfill u^\phi (t)& :=& (u_j(t+\phi _j))_{j\mathrm{\Gamma }}\hfill & & \mathrm{\Sigma }\hfill \end{array}$$
(1.8)
is a $`2\pi `$-periodic solution of (1.4), (1.7), for arbitrary phases $`\phi T^{m+1}`$. Eliminating one phase angle $`\phi _{m+1}`$ by passing to an associated Poincaré map, an $`m`$-dimensional manifold of fixed points arises, parametrized by the remaining $`m`$ phase angles. Assuming, in addition to the diagonal oddness symmetry (1.5), equivariance of (1.4) with respect to an $`S^1`$-action, the Poincaré map can in fact be obtained as a time-$`2\pi `$ map of an autonomous flow within the Poincaré section. In suitable notation, $`y=\{\phi _1,\mathrm{},\phi _m\}`$, the fixed point manifold then becomes an $`m`$-dimensional manifold of equilibria, as presented in (1.1), (1.2) above. For more detailed discussions of this example in the context of bifurcations without parameters see .

As a second example of equilibrium manifolds we consider viscous profiles $`u=u((\xi st)/\epsilon )`$ of systems of nonlinear hyperbolic conservation laws and stiff balance laws

$$\begin{array}{ccc}\hfill _tu+_\xi F(u)+\epsilon ^1G(u)& =& \epsilon _\xi ^2u.\hfill \end{array}$$
(1.9)
Viscous profiles then have to satisfy an $`\epsilon `$-independent ODE system

$$\begin{array}{ccc}\hfill \ddot{u}& =& (F^{}(u)s\mathrm{id})\dot{u}+G(u).\hfill \end{array}$$
(1.10)
Standard conservation laws, for example, require $`G0`$. The presence of $`m`$ conservation laws corresponds to nonlinearities $`G`$ with range in a manifold of codimension $`m`$ in $`u`$-space. Typically, then, $`G(u)=0`$ describes an equilibrium manifold of dimension $`m`$ of pairs $`(u,\dot{u})=(u,0)`$, in the phase space of (1.10). For an analysis of this example in the context of bifurcations without parameters see . For another example, which relates binary oscillations in central difference discretizations of hyperbolic balance laws with diagonal uncoupling of coupled oscillators, see .

We conclude our introductory excursion with a brief summary of some further examples. In , lines of equilibria have been observed for the dynamics of models of competing populations. This included a first partial analysis of failure of normal hyperbolicity.

A topologically very interesting example in compact three-dimensional manifolds involves contact structures $`\eta (\xi )`$ (i.e., nonintegrable plane fields and gradient vector fields $`\dot{\xi }=V(\xi )\eta (\xi )`$. See  for an in-depth analysis. Examples include mechanical systems with nonholonomic constraints. Notably, level surfaces of regular values of the potential $`V`$ consist of tori. Under a nondegeneracy assumption, equilibria form embedded circles, that is, possibly linked and nontrivial knots.

For a detailed study of plane Kolmogorov fluid flows in the presence of a line of equilibria with a degeneracy of Takens-Bogdanov type and an additional reversibility symmetry, see .

As a caveat we repeat that lines of equilibria, which are transverse to level surfaces of preserved quantities $`\lambda `$ do not provide bifurcations, without parameters. Rather, $`\dot{y}=0`$ for $`y:=\lambda `$ exhibits this problem as belonging to standard bifurcation theory; see (1.3).

## 2 . Sample vector fields and resulting flows

In this section we collect relevant example vector fields (1.1), (1.2) with lines and planes of equilibria $`x=0`$; see  for further details. We illustrate and comment the resulting flows.

Normally hyperbolic equilibrium manifolds admit a transverse $`C^0`$-foliation with hyperbolic linear flows in the leaves. See for example ,  and the ample discussion in . As a first example, we therefore consider

$$\begin{array}{ccc}\hfill \dot{x}& =& xy,\hfill \\ \hfill \dot{y}& =& x.\hfill \end{array}$$
(2.1)
Note the loss of normal hyperbolicity at $`x=y=0`$, due to a nontrivial transverse eigenvalue zero of the linearization. Clearly $`dx/dy=y`$, and the resulting flow lines are parabolas; see Figure 1. For comparison with standard bifurcation theory, where $`y=\lambda `$, we draw the $`y`$-axis of equilibria horizontally.

As a second example, consider

$$\begin{array}{ccc}\hfill \dot{x}& =& xy,\hfill \\ \hfill \dot{y}& =& \pm x^2.\hfill \end{array}$$
(2.2)
Again, a transverse zero eigenvalue occurs – this time with an additional reflection symmetry $`yy`$. Dividing by the Euler multiplier $`x`$, the reflection becomes a time reversibility. See the left parts of Figure 2 for the resulting flows. Note the resulting integrable, harmonic oscillator case which originates from the elliptic sign $`\dot{y}=x^2`$.

As a third example, we consider $`x=(x_1,x_2)\mathrm{𝖨𝖱}^2`$, $`y\mathrm{𝖨𝖱}`$ with a line $`x=0`$ of equilibria and a purely imaginary nonzero eigenvalue $`i\omega `$ at $`x=0`$. Normal-form theory, for example as in , then generates an additional $`S^1`$-symmetry by the action of $`\mathrm{exp}(i\omega t)`$ in the $`x`$-eigenspace. This equivariance can be achieved, successively, up to Taylor expansions of any finite order. In polar coordinates $`(r,\phi )`$ for $`x`$, an example of leading order terms is given by

$$\begin{array}{ccc}\hfill \dot{r}& =& ry,\hfill \\ \hfill \dot{y}& =& \pm r^2,\hfill \\ \hfill \dot{\phi }& =& \omega .\hfill \end{array}$$
(2.3)
Since the first two equations in (2.3) coincide with (2.2), the dynamics is then obtained by simply rotating the left parts of Figure 2 around the $`y`$-axis at speed $`w`$. The right parts of Figure 2 provide three-dimensional views of the effects of higher-order terms which do not respect the $`S^1`$-symmetry of the normal forms. In the elliptic case (b), all nonstationary orbits are heteroclinic from the unstable foci, at $`y>0`$, to stable foci at $`y<0`$. The two-dimensional respective strong stable and unstable manifolds will split, generically, to permit transverse intersections.

Our fourth example addresses Takens-Bogdanov bifurcations without parameters. In suitable rescaled form it reads

$$\begin{array}{ccc}\hfill \dot{\text{}}\dot{\text{}}\dot{\text{}}y+y\dot{y}& =& \epsilon ((\lambda y)\ddot{y}+b\dot{y}^2)\hfill \end{array}$$
(2.4)
with fixed parameters $`b`$, $`\lambda `$ and $`\epsilon `$. The $`y`$-axis as equilibrium line is complemented by the two transverse directions $`x=(\ddot{y},\dot{y})`$. Note the algebraically triple zero eigenvalue, double in the transverse $`x`$-directions, for $`\lambda =y=0`$. Two examples of the resulting dynamics for small positive $`\epsilon `$ are summarized in Figure 3.

The coordinates $`\tau `$ and $`\stackrel{~}{H}`$ in Figure 3 are adapted to the completely integrable case $`\epsilon =0`$. Indeed, obvious first integrals are then given by $`\mathrm{\Theta }=\ddot{y}+\frac{1}{2}y^2`$ and $`H=\frac{1}{2}\dot{y}^2y\ddot{y}\frac{1}{3}y^3`$. Coordinates are $`\tau =\mathrm{log}\mathrm{\Theta }`$ and $`\stackrel{~}{H}=\mathrm{\Theta }^{3/2}H`$, not drawn to scale. Parameters are $`\epsilon ,\lambda >0`$ and, for the hyperbolic case, $`17/12<b<1`$. For the elliptic case we consider $`b>1`$. The equilibrium $`y`$-axis, a cusp in $`(\mathrm{\Theta },H)`$ coordinates, transforms to the top (saddles) and bottom (foci) horizontal boundaries $`\stackrel{~}{H}=\pm \frac{2}{3}\sqrt{2}`$, with $`y=0`$ relegated to $`\tau =\mathrm{}`$. Since $`\tau `$ and $`\stackrel{~}{H}`$ are constants of the flow, for $`\epsilon =0`$, they represent slow drifts on the unperturbed periodic motion, for small $`\epsilon >0`$ and $`\stackrel{~}{H}<\frac{2}{3}\sqrt{2}`$. The top value $`\stackrel{~}{H}=+\frac{2}{3}\sqrt{2}`$ also represents homoclinics to the saddles, for $`\epsilon =0`$.

Along the focus line $`\stackrel{~}{H}=\frac{2}{3}\sqrt{2}`$ we observe Hopf bifurcations without parameters, corresponding to $`y=\lambda >0`$. The value of $`b`$ distinguishes elliptic and hyperbolic cases. In addition, lines of saddle-focus heteroclinic orbits and isolated saddle-saddle heteroclinics are generated, for $`\epsilon >0`$, by breaking the homoclinic sheets of the integrable case. Note in particular the infinite swarm of saddle-saddle heteroclinics, in the hyperbolic case.

As a final, fifth example we consider a reversible Takens-Bogdanov bifurcation without parameters:

$$\begin{array}{ccc}\hfill \dot{\text{}}\dot{\text{}}\dot{\text{}}y+(13y^2)\dot{y}& =& ay\ddot{y}+b\dot{y}^2.\hfill \end{array}$$
(2.5)
Here we fix $`a,b`$ to be small. Again $`x=(\ddot{y},\dot{y})`$ denotes the directions transverse to the equilibrium $`y`$-axis. Time reversibility generates solutions $`y(t)`$ from solutions $`y(t).`$ For two examples of the resulting dynamics see Figure 4. Coordinates are the obvious first integrals, for $`\epsilon =0`$, given by $`\mathrm{\Theta }=\ddot{y}+yy^3`$ and $`H=\ddot{y}y+\frac{1}{2}\dot{y}^2+\frac{3}{4}y^4\frac{1}{2}y^2`$. Note the two Takens-Bogdanov cusps, separated by a Hopf point along the lower arc of the equilibrium “triangle”. Compare Figures 2, 3. The elliptic Hopf point (b) arises for $`a(ab)>0`$, whereas $`a(ab)<0`$ in the hyperbolic case (a). Also note the associated finite and infinite swarms of saddle-saddle heteroclinics, respectively.

## 3 . Methods

Pictures are not proofs. What has been proved, then, and how? We use ingredients involving algebra, analysis, and numerical analysis, as we outline in this section. For further details see .

In a first algebraic step, we derive normal forms for vector fields with lines or planes of equilibria, assuming spectral degeneracies of the linearization $`A`$ in transverse directions. The spectral assumptions on $`A`$ in fact coincide with those established for parameter-dependent matrix families in standard bifurcation theory. This is reflected in the naming of the five examples of section 2.

There are more or less standard procedures to derive normal forms of vector fields. By suitable polynomial diffeomorphisms, certain Taylor coefficients of the vector field are eliminated, successively, for higher and higher order. See for example  for a systematic choice of normal forms, particularly apt for introducing equivariance of the nonlinear normal-form terms under the action of $`\mathrm{exp}(A^Tt).`$ Normal forms are, however, nonunique in general and other choices are possible.

In the present cases, we adapt the normal-form procedure to preserve the locally flattened equilibrium manifolds. Although the approach in  can be modified to accommodate that requirement, it did not provide vector fields convenient for subsequent analysis of the flow. A systematic approach to this combined problem is not known, at present. All examples (2.1)–(2.5) represent truncated normal forms. For the derivation of specific normal forms, for example of (2.4), to any order, and of example (2.5), to third order, see , , respectively.

Subsequent analysis of the normal-form vector fields is based on scalings, alias blow-up constructions. This is the origin, for example, of the small scaling parameter $`\epsilon `$ in the Takens-Bogdanov example (2.4). In passing, we note a curious coincidence of two view points for (2.4), concerning the roles of the equilibrium coordinate $`y\mathrm{𝖨𝖱}`$ and the “fixed” real parameter $`\lambda \mathrm{𝖨𝖱}`$. First, we may consider $`\lambda `$ as a parameter, with a line of equilibria $`y`$ associated to each fixed $`\lambda `$. Then (2.4) describes the collision of a transverse zero eigenvalue at $`y=0`$, as in (2.1), with imaginary Hopf eigenvalues at $`y=\lambda >0`$, as in (2.3), as $`\lambda `$ decreases through zero. Alternatively, we may consider normal forms for a plane $`𝐲=(y_1,y_2)`$ of equilibria with a transverse double zero eigenvalue, at $`𝐲=0`$. It turns out that the two cases coincide, after a scaling blow-up, up to second order in $`\epsilon `$, via the correspondence $`y=y_1,\lambda =y_2`$.

The core of any successful flow analysis in bifurcation theory is an integrable vector field; see again section 2. The issues of nonintegrable perturbations, by small $`\epsilon >0`$, and of omitted higher order terms, not in normal form, both ensue. Since the underlying integrable dynamics is periodic or homoclinic, in examples (2.3)–(2.5), averaging procedures apply. Indeed, $`\epsilon >0`$ then introduces a periodically forced, slow flow on first integrals, like $`(\mathrm{\Theta },H)`$, characteristic of $`\epsilon =0`$. We therefore derive an appropriate, but autonomous Poincaré flows, on $`(\mathrm{\Theta },H)`$, such that the associated true Poincaré map can be viewed as a time discretization of first order and step size $`\epsilon `$. In the unperturbed periodic region, this amounts to averaging, while the Poincaré flow indicates Melnikov functions at homoclinic or heteroclinic boundaries. The exponential averaging results by Neishtadt , for example, then imply that the separatrix splittings in the elliptic Hopf case (b), indicated in Figure 2, are exponentially small in the radius $`r`$ of the split sphere, for analytic vector fields. See also .

Lower bounds of separatrix splittings have not been established, in our settings. This problem is related to the very demanding Lazutkin program of asymptotic expansions for exponentially small separatrix splittings. For recent progress, including the case of Takens-Bogdanov bifurcations for analytic maps, see  and the references there. In absence of rigorous lower bounds, our figures indicate only simplest possible splitting scenarios.

While the splitting near elliptic Hopf points are exponentially small, the discretization of the Poincaré flow also exhibits splittings of the unperturbed saddle homoclinic families, which are of first order in the perturbation parameter $`\epsilon `$, in example (2.4), or in the small parameters $`a,b`$, in example (2.5). Explicit expressions have been derived for the Melnikov functions associated to these homoclinic splittings, in terms of elliptic function in case (2.4), and even of elementary functions in case (2.5). Simplicity and uniqueness of zeros of the Melnikov functions, however, has only been confirmed numerically. While this does not, strictly speaking, match an analytic proof, it still at least supports the validity of the scenarios summarized in Figures 3 and 4.

## 4 . Interpretation and perspective

We indicate some consequences of the above results for the examples of coupled oscillators, viscous shock profiles, and Kolmogorov flows indicated in section 1. We conclude with a few remarks on the future perspective of bifurcations without parameters.

We first return to the example (1.4)–(1.8) of a coupled oscillator square, $`m=1.`$ Equilibria $`y`$ of the Poincaré flow then indicate decoupled antipodal periodic pairs, say with phase difference $`y+c.`$ The case of a transverse zero eigenvalue, (2.1) and Figure 1, then indicates a 50% chance of recovery of decoupling with a stable phase difference $`y<0`$, locally, even when the stability threshold $`y=0`$ has been exceeded. The hyperbolic Hopf case (2.2), Figure 2 (a), illustrates immediate oscillatory loss of decoupling stability by transverse imaginary eigenvalues. A 100% recovery of decoupling stability, in contrast, occurs at elliptic Hopf points; see (2.2), Figure 2 (b). The exponentially small Neishtadt splitting of separatrices indicates a very delicate variability in the asymptotic phase relations of this recovery, for $`t\pm \mathrm{}`$. See . The Takens-Bogdanov cases (2.4), Figures 3 (a), (b) can then be viewed as consequences of a mutual interaction, of a transverse zero eigenvalue with either Hopf case, for recovery of stable decoupling.

In the example (1.9), (1.10) of stiff balance laws, elliptic Hopf bifurcation without parameters as in (2.4), Figure 3 (b), indicates oscillatory shock profiles $`u(\tau ),\tau =(\xi st)/\epsilon .`$ Such profiles in fact contradict the Lax condition, being over-compressive, and violate standard monotonicity criteria. For small viscosities $`\epsilon >0`$, weak viscous shocks in fact turn out unstable, in any exponentially weighted norm, unless they travel at speeds $`s`$ exceeding all characteristic speeds. The oscillatory profiles can be generated, in fact, by the interaction of inherently non-oscillatory gradient flux functions $`F(u)`$ with inherently non-oscillatory gradient-like kinetics $`G(u)`$ in systems of $`dim3`$. See .

The problem of plane stationary Kolmogorov flows asks for stationary solutions of the incompressible Navier-Stokes equations in a strip domain $`(\zeta ,\eta )\mathrm{𝖨𝖱}\times [0,2\pi ]`$, under periodic boundary conditions in $`\eta `$; see . An $`\eta `$-periodic external force $`(F(\eta ),0)`$, is imposed, acting in the unbounded $`\zeta `$-direction. Kolmogorov chose $`F(\eta )=\mathrm{sin}\eta `$. The Kirchgässner reduction  captures all bounded solutions which are nearly homogeneous in $`\zeta `$, in a center manifold spirit which lets us interpret $`\zeta `$ as “time”. The resulting ordinary differential equations in $`\mathrm{𝖨𝖱}^6`$ reduce to $`\mathrm{𝖨𝖱}^3`$, by fixing the values of three first integrals. A line of $`\zeta `$-homogeneous equilibria appears, in fact, and Kolmogorov’s choice corresponds to example (2.5) with $`a=b=0`$, to leading orders. In particular note the double reversibility, then, under $`y(t)\pm y(t)`$ which is generated by

$$\begin{array}{ccc}\hfill F(\eta )& =& F(\eta +\pi ),\text{and}\hfill \\ \hfill F(\eta )& =& F(\eta ).\hfill \end{array}$$
(4.1)
As observed by Kolmogorov, an abundance of spatially periodic profiles results. The sample choice $`F(\eta )=\mathrm{sin}\eta +c\mathrm{sin}2\eta `$, in contrast, which breaks the first of the symmetries in (4.1), leads to (2.5) with $`b=0<a`$, alias an elliptic reversible Takens-Bogdanov point without parameters; see Figure 4 (b). In particular, the set of near-homogeneous bounded velocity profiles of the incompressible, stationary Navier-Stokes system is then characterized by an abundance of oscillatory heteroclinic wave fronts, which decay to different asymptotically homogeneous $`\zeta `$-profiles, for $`\zeta \pm \mathrm{}`$. The PDE stability of these heteroclinic profiles is of course wide open.

As for perspectives of our approach, we believe to have examples at hand, from sufficiently diverse origin, to justify further development of a theory of bifurcations without parameters. In fact, transverse spectra $`\{0,\pm i\omega \}`$ and $`\{\pm i\omega _1,\pm i\omega _2\}`$ still await investigation before we can claim any insight into nonhyperbolicity of even the simple case of an equilibrium plane. This assumes the absence of further structural ingredients like symplecticity, contact structures, symmetries, and the like. Certainly our example collecting activities are far from complete, at this stage.

In addition, we have not addressed the issue of perturbations, so far, which could destroy the equilibrium manifolds by small drift terms. Examples arise, for example, when slightly detuning the basic frequencies $`2\pi /T_j`$ of our uncoupled oscillators or, much more generally, in the context of multiple scale singular perturbation problems. Feedback and input from our readers will certainly be most appreciated!