Dynamics near a periodically forced
robust heteroclinic cycle
July 19, 2019
Abstract.
There are few examples of nonautonomous vector fields exhibiting complex dynamics that may be proven analytically. We analyse a family of periodic perturbations of a weakly attracting robust heteroclinic network defined on the twosphere. We derive the first return map near the heteroclinic cycle for small amplitude of the perturbing term, and we reduce the analysis of the nonautonomous system to that of a twodimensional map on a cylinder.
When the perturbation strength is small there is an attracting invariant closed curve not contractible on the cylinder. Near the centre of frequency locking there are parameter values for which this curve coexists with an attracting periodic solution. Increasing the perturbation strength there are periodic solutions that bifurcate into a closed contractible invariant curve and into a region where the dynamics is conjugate to a full shift on two symbols. These last two dynamical features appear at a discretetime BogdanovTakens bifurcation.
Keywords: periodic forcing, heteroclinic cycle, global attractor, bifurcations, stability
2010 — AMS Subject Classifications Primary: 37C60; Secondary:34C37, 34D20, 37C27, 39A23, 34C28
1. Introduction
Archetypal examples of robust heroclinic cycles have been studied by Guckenheimer and Holmes [13], and May and Leonard [19], using a system of LotkaVolterra equations. The authors found saddleequilibria on the axes and attracting heteroclinic cycles and networks.
A heteroclinic cycle in an autonomous dynamical system consists of a connected union of saddletype invariant sets and heteroclinic trajectories connecting them. A heteroclinic network is a connected union of heteroclinic cycles. In equivariant systems the existence of invariant subspaces may force the existence of connecting trajectories between flowinvariant sets; heteroclinic cycles become robust in the sense that the connections persist under small symmetrypreserving perturbations. In generic dynamical systems without symmetry or other constraints, such configurations are structurally unstable.
In classical mechanics, dissipative nonautonomous systems received only limited attention, in part because it was believed that, in these systems, all trajectories tend toward Lyapunov stable sets (fixed points or periodic solutions). Evidence that second order equations with a periodic forcing term can have interesting behavior first appeared in the study of van der Pol’s equation, which describes an oscillator with nonlinear damping. Results of [8] pointed out an attracting set more complicated than a fixed point or an invariant curve. Levinson obtained detailed information for a simplified model [18].
Examples from the dissipative category include the equations of Lorenz, Duffing equation and Lorentz gases acted on by external forces [10]. The articles [9, 32] deal with heteroclinic tangles in timeperiodic perturbations in the dissipative context and show, for a set of parameters with positive Lebesgue measure, the existence of an attracting torus, of infinitely many horseshoes and of strange attractors with SRB measures.
While some progress has been made, both numerically and analytically, the number of differential equations with periodic forcing whose flows exhibit attracting heteroclinic networks, for which a rigorous global description of the dynamics is available, has remained small. To date there has been very little systematic investigation of the effects of perturbations that are timeperiodic, despite being natural for the modelling of many biological effects, see Rabinovich et al [24].
2. The object of study
For , the focus of this paper is on the following set of the ordinary differential equations with a periodic forcing:
(2.1) 
where
The amplitude of the perturbing term is governed by . We have chosen the perturbing term
It appears only in the first coordinate for two reasons: first, it simplifies the computations. Secondly, it allows comparison with previous work by other authors [3, 6, 12, 24, 31]. We denote the vector field associated to (2.1) by .
Remark 1.
The perturbation term may be replaced by where is any periodic and continuously differentiable function. In some places we use the property .
2.1. The unperturbed system ()
The dynamics associated to this equation has been systematically studied in [6]. For the sake of completeness, we recall its main properties. The vector field has two symmetries of order 2:
forming a symmetry group isomorphic to . The symmetry remains after the perturbation governed by . The unit sphere is flowinvariant and globally attracting. The points and are equilibria. From the symmetries , it follows that the planes and are flowinvariant, and hence they meet in two flowinvariant circles connecting the equilibria – see Figure 1. Since and , then these two equilibria are saddles, and there are heteroclinic trajectories going from each equilibrium to the other one. More precisely, the expanding and contracting eigenvalues and of the derivative of the vector field at are:
with . The origin is a repellor.
Considering the system restricted to , equivariance forces the invariant manifolds in of and to be in a very special position: they coincide. In , the invariant manifolds of and are onedimensional and contained in the invariant circles , , giving rise to a heteroclinic network . In the restriction to each of the invariant planes , the equilibria , have a saddlesink connection, so this network is persistent under perturbations that preserve the symmetry, and in this sense it is robust.
For all nonequilibrium points , we have , whereas for , we have , as in Figure 1.
The heteroclinic network is asymptotically stable by the Krupa and Melbourne criterion [16, 17]. Note that:
The constant measures the strength of attraction of the cycle in the absence of perturbations. There are no periodic solutions near because . Typical trajectories near the heteroclinic network spend increasing amounts of time near each saddle point on each occasion they return close to it. In some places we assume that we are in the weakly attracting case , we make the assumption explicitly when it is used. The case , corresponds to a resonant bifurcation of the robust heteroclinic cycle – this case has been explored by Postlethwaite and Dawes in [21, 22].
2.2. A pullback attractor for small
Kloeden and Rasmussen [14] have results connecting attractors for autonomous systems and their perturbations, that may be applied here. We have that is a global attractor of the autonomous flow (), the vector field is uniformly Lipschitz and the periodic perturbation term is limited. Then Section 11 of [14] allows us to conclude that the nonautonomous system (2.1) generates a process which has a pullback attractor such that
where is the euclidean distance on . Moreover, for a given , the sets have the same Haussdorf dimension for all . This suggests that solutions of the perturbed system (2.1) should make excursions around the ghost of . In this article we explore the resulting dynamics.
Dawes and T.L. Tsai [12, 30, 31], presented preliminary results on the perturbation of the example studied in [13]. They identified distinct dynamical regimes depending on whether is or not close to 1. Here we deal with the case in general. Our results provide insight into the dynamics of a nonautonomous periodic forcing of an autonomous equation with a compact attractor.
Our purpose in writing this paper is not only to point out the range of phenomena that can occur when simple nonlinear equations are periodically forced, but to bring to the foreground the techniques that have allowed us to reach these conclusions in a relatively straightforward manner. These techniques are clearly not limited to the systems considered here. It is our hope that they will find applications in other dynamical systems, particularly those that arise naturally from mechanics or physics.
2.3. Main results and structure of the paper
We now describe briefly the main results and the contents of this paper. We provide informal statements of the results, that will be made rigorous in the course of the exposition. Since there is a large number of constants and parameters used in the article, we have included a list of notation as an appendix.
Expressions for the first return map to a section transverse to the connection are obtained in Section 3 for the case and in Section 4 for the general case.
We linearise the autonomous equations around the equilibria in § 3.1 to construct a first return map. In Section 4, we begin by presenting a systematic calculation of the first return map for the robust heteroclinic cycle subjected to the nonnegative periodic forcing function . By including the timedependent terms through all steps in the calculation, we obtain a first return map, that can be quantitatively compared to the dynamics of the original differential equations associated to – see Remarks 3, 4 and 6. These comparisons show that the new return map associated to captures the dynamics well. The Poincaré map for (2.1), described in Theorem A, yields a description of the dynamics in terms of a twodimensional map for the coordinate and the return time at which trajectories reach the crosssection:
Theorem A.
In the weakly attracting case , the first return map to a given crosssection, of the flow defined by (2.1), may be approximated by a map of the form:
The values of the constants for Theorem A are given in the Appendix. Theorem A is proved in Section 4. Although this result is valid for the weakly attracting case, we may study the dynamics of for all . The expression of coincides with that obtained by Tsai and Dawes [31] for a different system, here we clarify and extend their results. Note that, because is periodic in the variable , the natural phase space for is the cylinder
(2.2) 
The consequences of the interaction between and in the original system are not trivial. In the weakly attracting case, if is small and , then the map has both stable and unstable fixed points (corresponding to stable/unstable periodic solutions for the flow), say and . For the nonautonomous equation (2.1), solutions might fluctuate by small amounts near the periodic orbit of an amplitude that would be anticipated from considering the timeaveraged value of the forcing function. Based on the Annulus Principle [4], we obtain:
Theorem B.
There is an open set of parameters for which the maximal attractor for is an invariant closed curve not contractible on the cylinder.
Section 6 is concerned with bifurcation and stability of periodic solutions of (2.1). For this, we first introduce an auxiliary parameter and look for fixed points of . A first step is the following:
Proposition C.
Consider the problem of finding fixed points of with , , . The curves and and with separate the first quadrant of the plane into five regions (see Figure 2) corresponding to the following behaviour:

and — there are no fixed points;

and — there are fixed points for each in a closed interval;

and — there are fixed points for each in the union of two disjoint closed intervals;

and — there are fixed points for each ;

and — there are fixed points for each in the complement of a limited open interval.
Moreover, when fixed points exist for in an interval, then there are two fixed points for each in the interior of the interval and only one for in the boundary.
Besides the finding of fixed points of , in Section 6 we also discuss their stability and bifurcations between the different regions of the parameter space .
Corollary D.
For , , , fixed points of undergo saddlenode bifurcations on the surfaces in the three dimensional parameter space given by
In this context, we find in § 6.3 an organising centre for the dynamics of as follows:
Theorem E.
There are curves in the three dimensional parameter space where fixed points of undergo a discretetime BogdanovTakens bifurcation, a single curve for , two curves for .
In particular, we conclude that there are contractible closed invariant curves arising at Hopf bifurcations and there exist small regions in parameter space where has chaotic and nonhyperbolic dynamics. This is a consequence of the discretetime BogdanovTakens bifurcation studied by Broer et al [7] and Yagasaki [33, § 3]. The stability of bifurcating solutions is studied in § 6.4.
In § 6.5 these results are used to find periodic solutions of the original differential equations whose period is an integer multiple of the period of the forcing terms — frequency locked solutions. This agrees well with numerics presented in [30]. We show their existence for different values of , as well as the existence of invariant tori and chaotic regions. We also show that there is no gain in looking for different multiples , of the period, because we obtain essentially the same solution for all . This is summarised as follows:
Theorem F.
If , then there are two frequency locked solutions of (2.1) with period , , for the following values of , according to the regions in Proposition C:

;

and ;

and ;

all ;
where the for have the values of Proposition C.
There are values of and values and such that, for each :
Finally we analyse in Section 7 the dynamics near special solutions with frequency locked periods, called the centres of frequency locking.
Theorem G.
Near the centre of the frequency locking, the dynamics of is approximated by the discretisation of the equation for a damped pendulum with constant torque
An explicit expression for and in terms of the parameters in is given in Section 7. The effects of timeperiodic forcing on the damped pendulum was extensively studied by Andronov et al [1] and we use their results to obtain information about our problem.
As a corollary we conclude that there is an open region in the space of parameters where solutions of (2.1) are attracted either to a stable nontrivial periodic solution or to an attracting torus. These two flowinvariant sets may coexist and attract open sets of trajectories. We say that the system exhibits bistability. The physical meaning is explained in [11].
Dynamics similar to that which we have described is expected to occur generically near periodically forced robust weakly attracting heteroclinic cycles. This is why the result of our computations may be applied to other similar cases. The understanding of the effects of different classes of perturbation is far from being done. We finish the article with a short discussion in Section 8 of the consequences of our findings, both for the map and for the equation (2.1).
3. Timeindependent first return map
We will define four crosssections transverse to all trajectories in a neighbourhood of . Repeated intersections to a fixed crosssection define a return map from the section to itself; studying the dynamics of this map enables us to understand the dynamics of trajectories near .
We construct the return map as the composition of two types of map: local maps within neighbourhoods of the saddletype equilibrium points where the dynamics can be well approximated by the flow of the linearised equations, and transition maps from one neighbourhood to another (also called global maps).
Near the equilibrium , the crosssections are denoted and .
3.1. Linearisation
3.2. The crosssections
Consider cubic neighbourhoods and in of and , respectively, given in the coordinates of (3.3) by (see Figure 3):
for small. The boundary of consists of six connected components.

Two squares, the top and the bottom of the cube parametrised by , where the flow enters in the radial direction.

The set of points that go inside in positive time, that has two components, one of which is

The set of points that go outside in positive time, with two components, one of which is
Similarly, the boundary of consists of six connected components.

Two squares, the top and the bottom of the cube parametrised by , where the flow enters in the radial direction.

The set of points that go inside in positive time, that has two components, one of which is

The set of points that go outside in positive time, with two components, one of which is
3.3. Local map near
The solution of the linearised system (3.3) near with initial conditions , in is:
The time of flight of a trajectory inside is the solution of
(3.4) 
Therefore, the transition map in coordinates and is:
(3.5) 
3.4. Local map near
The solution of the linearised system (3.3) near with initial conditions , in is:
The time of flight from to is the solution of :
in coordinates and is:
(3.6) 
3.5. The global maps
An approximation to the maps
is to take them as the identity. In coordinates we obtain
3.6. First return map for the unperturbed equation
The first return map to is well defined at all points . After a linear rescaling of the variables, we may assume that and obtain
(3.7) 
Either the first or the second coordinate dominates, depending on the relative size of the exponents in , i.e. depending on the sign of , see Figure 4. The transition between the boundaries of and occurs in a flowbox, hence the transition time is limited above and below. We assume the transitions far from the equilibria are instantaneous, and then the time of the first return of the point with is given by . Taking into account the transition times out of and would approximately change the value of by a constant.
4. Timedependent first return
The aim of this section is to obtain an expression for the first return map to when , that will be denoted by . When , the map should coincide with defined in (3.7). Specifically, we prove
Theorem 1.
If , and , the first return map to the crosssection , of the flow defined by (2.1), may be approximated by the map:
(4.8) 
where and .
Because of the timeperiodic perturbation, the local linearisation now includes timedependent terms that are important in the accurate calculation of the local map.
At each step, we calculate not only the point where a solution hits each crosssection but also the time the solution takes to move between crosssections. As in § 3.6, the time spent close to the connections is small compared with the time spent near the equilibria, especially when is small. This time is not taken into account in the calculations. Note that the equilibria for the vector field (associated to the equation (2.1) when ) are no longer equilibria for , but the crosssections remain transverse.
4.1. Linearization
The linearisation near and may be written as:
(4.9) 
and
(4.10) 
respectively, with .
Equation (2.1) may be written in the form for and where is any of the equations (4.9) or (4.10). In this form we have that the constant matrix has no eigenvalues with a zero real part, the perturbaton is limited and the nonlinear part is limited and uniformly Lipschitz in a compact neighbourhood of . Under these conditions, we may use Palmer’s Theorem [20, pp 754] to conclude that there is a small neighbourhood of and where the vector field is conjugate to its linear part. As before, let us label the neighbourhoods by and , respectively. The terminology for the boundary sections in and will be the same as in § 3.2.
Remark 2.
If , according to the Lagrange method of variation of parameters – see [29, pp 842], the general solution of:
is
4.2. Local map near
Let us describe the general solution of (4.9). For we get By Remark 2, the solution of the linearised system (4.9) near is:
(4.11) 
where are the initial conditions in . The time of flight is the solution of :
Therefore, the arrival time depends on and it is given by:
(4.12) 
In particular, we may write:
which is equivalent to
Replacing by of (4.12) in the first equation of (4.11), we get:
where as defined in § 2.1. Therefore, we may write:
Remark 3.
Note that when and , the first, second and third components coincide with formulas (3.4) and the second component of .
4.3. Local map near
The treatment of (4.10) is similar to § 4.2. From , it follows that . The solution of (4.10) is:
The time of flight from to is the solution of :
(4.13) 
This is difficult to solve, so we compute the Taylor expansion of as function of . It is easy to see that . Differentiating equation (4.13) with respect to , and evaluating at , it yields:
implying that:
Thus, truncating at second order in we obtain:
Since , setting , as in § 2.1, we get:
and then:
Adding up, we get:
Concerning the coordinate , we may write:
On the other hand, the following equalities hold:
Therefore, is given by:
Remark 4.
When , and , the last two components of the previous map coincide with the expression given in (3.6).
4.4. Discussion of the time dependence
We are assuming in (2.1) that , hence the term that appears inside the integrals in tends to zero as goes to – see Figure 5. For large times we may take the contribution of this integral to be time independent.
Setting
we get
On the other hand, the assumptions in (2.1) also include implying that the coefficient . Hence the exponent that appears in increases with and the integral cannot be ignored. In order to obtain estimates for the integral, let , for , be given by:
Lemma 2.
If , then:
Proof.
Integrating by parts twice, we obtain:
Hence,
which is equivalent to the expression in the statement. ∎
Lemma 3.
If and as in the expression of , we have