Global manifold structure of a continuous-time heterodimensional cycle
Andy Hammerlindl Bernd Krauskopf Gemma Mason Hinke M. Osinga




Global manifold structure of a 
 continuous-time heterodimensional cycle




Andy Hammerlindl^1, Bernd Krauskopf^2, Gemma Mason^2 
 and Hinke M. Osinga^2






^1 School of Mathematical Sciences, Monash University, Melbourne VIC 3800, Australia 
 
^2 Department of Mathematics, The University of Auckland, Auckland 1142, New Zealand



















June 2019







A heterodimensional cycle consists of a pair of heteroclinic connections between two saddle periodic orbits with unstable manifolds of different dimensions. Recent theoretical work on chaotic dynamics beyond the uniformly hyperbolic setting has shown that heterodimensional cycles may occur robustly in diffeomorphisms of dimension at least three. We study a concrete example of a heterodimensional cycle in the continuous-time setting, specifically in a four-dimensional vector field model of intracellular calcium dynamics. By employing advanced numerical techniques, Zhang, Krauskopf and Kirk [Discr. Contin. Dynam. Syst. A 32(8) 2825–2851 (2012)] found that a heterodimensional cycle exists in this model. 

We investigate the geometric structure of the associated stable and unstable manifolds in the neighbourhood of this heterodimensional cycle, consisting of a single connecting orbit of codimension one and an entire cylinder of structurally stable connecting orbits between two saddle periodic orbits. We employ a boundary-value problem set-up to compute their stable and unstable manifolds, which we visualize in different projections of phase space and as intersection sets with a suitable three-dimensional Poincaré section. We show that, locally near the intersection set of the heterodimensional cycle, the manifolds interact as described by the theory for three-dimensional diffeomorphisms. On the other hand, their global structure is more intricate, which is due to the fact that it is not possible to find a Poincaré section that is transverse to the flow everywhere. Our results show that the abstract concept of a heterodimensional cycle arises and can be studied in continuous-time models from applications.














§ INTRODUCTION



One of the best-studied forms of chaotic dynamics is uniform hyperbolicity, which arose from the work of Anosov and Smale in the 1960s <cit.>. In a uniformly hyperbolic system, every orbit tends toward one of a finite number of basic sets in phase space. Each point in a basic set has well-defined directions of expansion and contraction, and tangent to these directions are stable and unstable manifolds through the point. This form of dynamics is robust, that is, any small perturbation of a uniformly hyperbolic system is still uniformly hyperbolic. Hence, it is likely that such chaotic systems may be observed in physical processes, and evidence of such behaviour has been found in a number of settings; for example, see <cit.>. Shortly after its discovery, Smale and others conjectured that uniform hyperbolicity was the “default” form of chaotic dynamics. That is, they conjectured that any chaotic dynamical system could, after a small perturbation, be turned into a uniformly hyperbolic system. Consequently, most chaotic behaviour identified in a mathematical model of a physical process would be uniformly hyperbolic. Unfortunately, this hypothesis turned out to be false and other forms of robust chaotic dynamics have since been identified; see <cit.>. Since these more complicated forms of non-hyperbolic chaotic dynamics also persist under perturbation, they could conceivably be observed in physical processes as well.





One mechanism that creates non-hyperbolic chaotic dynamics is a heterodimensional cycle <cit.>. Consider a discrete-time dynamical system, a diffeomorphism, with two hyperbolic fixed (or periodic) points γ_1 and γ_2. These points form a cycle if the stable manifold W^s(γ_1) of γ_1 intersects the unstable manifold W^u(γ_2) of γ_2 and the stable manifold W^s(γ_2) of γ_2 intersects the unstable manifold W^u(γ_1) of γ_1. The cycle is heterodimensional if γ_1 has a different index from γ_2, that is, their stable manifolds have different dimensions. This requires a diffeomorphism on a phase space of dimension at least three. fig:hdcycle shows a sketch of a heterodimensional cycle for this case of lowest dimension, where W^s(γ_1) and W^u(γ_2) are two dimensional and W^u(γ_1) and W^s(γ_2) are one dimensional. The intersection B between the surfaces W^s(γ_1) and W^u(γ_2) is transverse and, therefore, persists and is structurally stable under perturbation of the dynamics. The one-dimensional manifolds W^u(γ_1) and W^s(γ_2), on the other hand, intersect non-transversely in the heteroclinic orbit  ( a_k )_k ∈. Indeed,  this intersection is of codimension one and can be destroyed by a perturbation of the system. Nevertheless, it has been proven that, near a given heterodimensional cycle as shown in fig:hdcycle, the existence of a heterodimensional cycle is a C^1-open property <cit.>. That is, the existence of any specific heterodimensional cycle is not itself a C^1-open property, but for any diffeomorphism f with a heterodimensional cycle and any sufficiently small ε, there is a C^1-open set U in the space of diffeomorphisms ε-close to f where every element of U has a heterodimensional cycle. This phenomenon is closely associated with objects called blenders, which are invariant sets of diffeomorphisms with the property, in the present context of dimension three, that their one-dimensional stable or unstable manifold cannot be avoided by one-dimensional curves from (an open neighbourhood of) a certain direction <cit.>. Blenders are in some sense a generalization of Smale's horseshoe construction to higher dimensions, and their robust defining property is a form of partially hyperbolic dynamics that has been used to prove the existence of heterodimensional cycles. 

Another mechanism to create non-hyperbolic chaotic dynamics is a homoclinic tangency between a stable and an unstable manifold of a single fixed or periodic point of a diffeomorphism <cit.>. Such a non-transverse intersection may also be a robust phenomenon when the phase space of the diffeomorphism is at least three <cit.>. Homoclinic tangencies and heterodimensional cycles are closely related; see <cit.>. A conjecture of Palis states that every dynamical system may be perturbed (in the C^r-topology of diffeomorphisms) to produce a system that is either uniformly hyperbolic or has a homoclinic tangency or a heterodimensional cycle <cit.>. Limited versions of this conjecture have been verified in a number of settings <cit.>.

The progress made in the abstract study of possible forms of chaotic dynamics has primarily focussed on diffeomorphisms, that is, on discrete-time dynamical systems. One major reason is that the phase-space dimension needed for discrete-time systems to create a particular type of chaos is always smaller than that needed for continuous-time systems. In many applications, on the other hand, mathematical models have continuous time and are given in the form of a system of ordinary differential equations (ODEs), that is, their dynamics is determined by a flow; see, for example, <cit.> as entry points to the extensive literature on models from applications. Hence, it is an important question how different types of dynamics that have been found in diffeomorphisms manifest themselves in flows.

Of course, there is a well-known connection between continuous-time and discrete-time systems. Given a diffeomorphism one can construct its suspension, which is an abstract flow with an additional direction that represents the continuous time and whose time-one map is the diffeomorphism. Clearly, periodic points of the diffeomorphism correspond to periodic orbits of the flow. On the other hand, a periodic orbit of a vector field corresponds to fixed points of the associated Poincaré return map, which is obtained by introducing a (local) cross section transverse to the periodic orbit and, hence, the flow nearby. Indeed, the stability and bifurcation analysis of periodic orbits of vector fields is derived from the corresponding stability and bifurcations of the associated fixed points <cit.>. More precisely, the dynamics of a vector field locally in the neighbourhood of a periodic orbit can be described exactly in terms of the dynamics of a diffeomorphism locally in the neighbourhood of the corresponding fixed point of the associated Poincaré return map, and vice versa. However, when the dynamics is more complicated, and certainly in parameter regions where chaotic behaviour occurs, it may be necessary to consider a rather larger domain of definition of the Poincaré return map. Unless the vector field is periodically forced, that is, the flow is a suspension, it is not possible to find a section (codimension-one submanifold) to which the flow is transverse at any point. Hence, the Poincaré return map can only be defined locally in regions where the flow is transverse to the section and, moreover, the flow returns back to the section <cit.>. Especially when one is interested in chaotic dynamics associated with homoclinic and heteroclinic bifurcations of closed periodic orbits <cit.>, the overall picture of how their invariant manifolds intersect a codimension-one section of interest cannot be determined from the information of an equivalent diffeomorphism alone. Indeed, advanced numerical methods have an important role to play in determining the global organization of phase space, how it changes with parameters and connects with known theory <cit.>.

In this paper, we investigate a heterodimensional cycle in a four-dimensional vector field that represents a simple model of intracellular calcium oscillations. The model is based on the so-called Atri model <cit.> and, in the form studied in <cit.>, it is given by the four differential equations


    {[                                                                                                           ċ                                                                                                           =                                                                                                          v,;                                                                                                       D  v̇                                                                                                           = s   v - ( α + k_f  c^2/c^2+ϕ_1^2  n ) 
                                      ( γ  (c_t + D   v - s   c)/s - c );                                                                                                                                                                                                                                                                                                        + k_s   c -   (J - k_p   c),;                                                                                                        ċ_̇ṫ                                                                                                           =                                                                                              (J - k_p   c),;                                                                                                        s  ṅ                                                                                                           =                                                                                      1/2( ϕ_2/ϕ_2 + c - n ) ].


for the concentration of calcium c inside the cell, its voltage v, the total calcium level c_t (including that stored in the so-called endoplasmic reticulum (ER), which is a safe internal storage of this toxic substance) and the gating variable n that describes transport of calcium through the membrane. The Atri model takes spatial variation into account, system eq:4DAtri is written in moving-frame coordinates and differentiation is with respect to the travelling-wave coordinate. Hence, periodic and solitary calcium waves are represented by periodic and homoclinic orbits of eq:4DAtri, respectively. We consider how the possible dynamics depends on the wave speed s and the parameter J, which represents a given flux of calcium entering the cell from the outside in. All other parameters are constant and their values are given in Table <ref>. 





The Atri model eq:4DAtri exhibits very complicated dynamics <cit.> and, in particular, Zhang, Krauskopf, and Kirk <cit.> realized that this system has a heterodimensional cycle between two periodic orbits with different indices. This global object was found numerically by careful computations with a new numerical method (an implementation of Lin's method) developed and presented in <cit.>. To our knowledge, system eq:4DAtri is still the only known explicit vector field arising from an application that exhibits a heterodimensional cycle between two periodic orbits, and it has the lowest required dimension. We will refer to this heterodimensional cycle also as a PtoP cycle (for Periodic orbit to Periodic orbit); note that this short-hand notation comes from the literature on computing connecting orbits of vector fields, where heteroclinic cycles between an equilibrium and a periodic orbit (such as the one that gives rise to the chaotic attractor in the Lorenz system) are also referred to as EtoP cycles <cit.>.

We present a case study of how the abstract theory of heterodimensional cycles manifests itself in a concrete vector field such as system eq:4DAtri. More specifically, we show how the respective global stable and unstable manifolds of the two periodic orbits intersect and give rise to the transverse and non-transverse parts of the heterodimensional cycle. We employ a state-of-the-art numerical approach based on boundary-value problem formulations <cit.> to compute and visualize selected families of orbit segments; these computations are conducted with the pseudo-arclength continuation package Auto <cit.>. We show three-dimensional projections of the four-dimensional phase space to illustrate one by one how the intersection sets of the different global manifolds with a suitable codimension-one section Σ arise. Taken together, this information makes clear how (the intersection set of) the heterodimensional cycle is generated, to what extent it represents the theory for diffeomorphisms, and what the differences are. More generally, the results presented here show that it is entirely feasible with state-of-the-art numerical methods to find and investigate objects as complex as a heterodimensional PtoP cycle in continuous-time models from applications. Hence, theoretical insights on new types of robust dynamics can now be investigated in terms of their relevance in specific contexts.

This paper is organized as follows. In sec:hetcycle we present the heterodimensional cycle in the four-dimensional phase space of system eq:4DAtri. sec:returnmap is devoted to the manifold structure that gives rise to the PtoP cycle. The codimension-one section Σ is introduced in sec:P, and sec:P_PtoP shows the intersection sets of the PtoP cycle with Σ. sec:ws then introduces the intersection set of the two-dimensional stable manifold and sec:wu that of the two-dimensional unstable manifold. How the two intersection sets interact in the three-dimensional section Σ is explained in sec:wswu, which also contains a comparison with the case of a diffeomorphism. We end in sec:conclusions with some conclusions and suggestions for future research.





§ THE HETERODIMENSIONAL PTOP CYCLE



In this paper, we explore the global manifold structure giving rise to a codimension-one heterodimensional PtoP cycle in system eq:4DAtri. The PtoP cycle was found in <cit.> for fixed s = s^∗ = 9.0 at J = J^∗≈ 3.02661 with an implementation of Lin's method, and its locus was then found as a curve in the (J, s)-plane. We start our investigation in the same way, where, as in <cit.>, all other parameters are as given in Table <ref>.





fig:hetcycle shows the respective part of the (J, s)-plane in panel (a) and two three-dimensional projections of the PtoP cycle at (J^∗, s^∗) in panels (b) and (c). In fig:hetcycle(a), two periodic orbits Γ_1 and Γ_2 bifurcate from the curve SL of saddle-node bifurcation of periodic orbits. These two periodic orbits are both hyperbolic, of saddle type and of different index in the region that is bounded by SL and the curves H of Hopf bifurcation, where Γ_2 bifurcates, and PD of period-doubling bifurcation, where Γ_1 bifurcates. Notice that the curve SL ends on H at a degenerate Hopf point DH of codimension two, while the curve PD is tangent to SL at a codimension-two saddle-node-period-doubling point PS. The heterodimensional cycle exists along the curve PtoP, which ends at a point (not labelled) on the curve PD. 

The point (J^∗, s^∗) is marked, and the corresponding heterodimensional PtoP cycle is shown in
panels (b) and (c) of fig:hetcycle. More specifically, these images show the PtoP cycle between the two periodic orbits Γ_1 and Γ_2, which have Floquet multipliers


  λ^u_1 ≈9.6 ×10^4,   
  λ^s_1 ≈2.8 ×10^-1,      
  λ^ss_1 ≈7.7 ×10^-3, 
 

and



  λ^uu_2 ≈7.0 ×10^3   
  λ^u_2 ≈1.3,      
  λ^s_2 ≈3.4 ×10^-1, 
 

respectively. Hence, Γ_1 has a two-dimensional unstable manifold W^u(Γ_1) and a three-dimensional stable manifold W^s(Γ_1), while Γ_2 has a three-dimensional unstable manifold W^u(Γ_2) and a two-dimensional stable manifold W^s(Γ_2). The three-dimensional manifolds W^s(Γ_1) and W^u(Γ_2) intersect transversely, as expected. The two-dimensional manifolds W^s(Γ_2) and W^u(Γ_2), on the other hand, have a non-transverse intersection at the special parameter pair (J^∗,  s^∗). We denote these non-empty intersections by


  A : = W^s(Γ_2) ∩W^u(Γ_1),      B := W^s(Γ_1) ∩W^u(Γ_2).


The intersection set A is one dimensional and consists of a single trajectory from Γ_1 to Γ_2. The intersection set B is two dimensional and forms a topological cylinder consisting of a one-parameter family of trajectories from Γ_2 to Γ_1. Note the strong contraction toward Γ_1 in backward time and the much slower approach toward Γ_2 in forward time, which is particularly visible in panel (c); this difference in speed is due to the fact that the ratio λ^s_2 / λ^ss_1 is of the order 10^2. The curve  PtoP in the (J, s)-plane in panel (a) was computed by parameter continuation of the non-transverse, codimension-one connecting orbit A, while checking that the cylinder B, that is, the entire heterodimensional cycle, also exists along PtoP.




 §.§ Numerical approach to finding the PtoP cycle



As mentioned in the introduction, the computation of the heterodimensional cycle for system eq:4DAtri requires the advanced numerical approach from <cit.>. Its major element is that it is based on Lin's method <cit.>, which has the important property that it defines a smooth test function — thereby  guaranteeing (under very mild starting conditions) that a particular homoclinic or heteroclinic orbit is found if it exists. We refer to <cit.> for the full details, but give a brief overview here of the steps involved. 

Numerical methods for approximating connecting orbits are not new; see <cit.> and references therein. Initially, the interest was in computing connections between equilibria. Depending on the dimensions of the respective eigenspaces, such connections may exist only for specific parameter values, which gives an added computational challenge. The connecting orbit is represented by an orbit segment truncated at both ends; it is assumed to start close to an equilibrium in its unstable eigenspace and end close to the same or other equilibrium in its stable eigenspace. The same approach can also be used for connections involving periodic orbits, where the begin or end point of the orbit segment is required to lie in the linear eigenbundle associated with the stable or unstable Floquet multipliers. The connecting orbit is then found as the orbit segment that solves a two-point boundary value problem (2PBVP) <cit.>; this is typically done in a continuation setting, where the specific connecting orbit is detected as a special solution from a family with a known solution. 

To achieve a well-posed and numerically robust set-up, Beyn <cit.> popularized the use of projection boundary conditions for the begin and end points of the orbit segment, which provides full error control. Various implementations exist, involving the continuation of invariant subspaces <cit.> or the (adjoint) variational equation <cit.>. The examples in these papers include computation of so-called EtoP and PtoP connections, that is, heteroclinic orbits between an equilibrium (E) and a periodic orbit (P), and homoclinic or heteroclinic orbits between one or two periodic orbits (PtoP), respectively; note that the example of a PtoP connection from <cit.> is between two periodic orbits with the same number of stable Floquet multipliers, so they are not part of a heterodimensional cycle. The continuation of the 2PBVP is started from a solution found by shooting <cit.> or via a homotopy step <cit.>; the latter is also the standard approach for the software package HomCont <cit.> that computes homoclinic and heteroclinic orbits of equilibria. 

The set-up based on Lin's method provides a natural way to find a first connecting orbit segment. The implementation for connections between equilibria was presented in <cit.> and the approach was extended to connections involving periodic orbits in <cit.>. The idea is to split the connecting orbit into two halves and require that both halves have an end point in a suitably chosen codimension-one section Π that is (locally) transverse to the flow. At the same time, it is required that the two orbit segments satisfy suitable projection boundary conditions at their other end points. Typically, the two end points in Π will not be the same, but if they are then a connecting orbit has been found. The key insight of Lin's method is that the gap between the two end points in Π, the so-called Lin gap, can be restricted to lie in a particular subspace Z ⊂Π of a dimension related to the codimension of the connecting orbit. In other words, the Lin gap can be closed in a component-wise manner with respect to a chosen basis for Z while adjusting the start or end point of one of the orbit segments within the required linear eigenspace. The computational set-up requires continuation of two orbit segments (double the system dimension), but this is outweighed by the major advantage that the two-segment 2PBVP is well posed even when the parameters are such that the particular connecting orbit does not exist; as long as Π is transverse to the flow for the family of orbit segments under consideration, a first solution to the 2PBVP can be found relatively easily and a possible special orbit segment with closed Lin gap can be searched for in the family. 

The codimension-one connecting PtoP orbit  A = W^u(Γ_1) ∩  W^s(Γ_2) from Γ_1 to Γ_2 in panels (b) and (c) of fig:hetcycle has been found in this way; here the Lin gap is of dimension one and closing it requires changing a system parameter, which was J in this case. The corresponding pair of connecting orbit segments can then be continued in system parameters with the Lin gap remaining closed, or as a single orbit segment after concatenation of the pair into a single orbit segment. The curve PtoP in fig:hetcycle(a) was computed in this way.

The cylinder B = W^u(Γ_2) ∩  W^s(Γ_1) from Γ_2 to Γ_1 in panels (b) and (c) of fig:hetcycle can be computed with effectively the same approach. One first finds a pair of orbits with projection boundary conditions in W^s(Γ_1) and W^u(Γ_2), respectively, with their other end points in Π. Again, a one-dimensional Lin gap can then be defined in Π, but to close it one does not need to change system parameters in this case. Rather, varying the distance from one of the periodic orbits closes the gap. Once a first connecting orbit has been found, the cyclinder can be swept out as a one-parameter familly while varying the distance from one of the periodic orbits (over one fundamental domain of the local return); see <cit.> for the more details.





§ GLOBAL MANIFOLD STRUCTURE OF THE PTOP CYCLE



Our goal is to investigate how the heterodimensional PtoP cycle Γ_1 ∪ A ∪Γ_2 ∪ B in the four-dimensional system eq:4DAtri arises from the intersections of the two-dimensional manifolds W^u(Γ_1) and W^s(Γ_2) and the three-dimensional manifolds of W^u(Γ_2) and W^s(Γ_1), respectively. Moreover, we wish to make the connection between this heterodimensional PtoP cycle in a flow and the abstract theory for diffeomorphisms. The challenge is to compute and visualize parts of interest of invariant manilofds that provide the insights that we seek. To this end, we again find and continue orbit segments that are defined by suitable 2PBVPs. Indeed, this approach is very versatile and accurate when it comes to computing two-dimensional (sub)manifolds. As before, the key idea is to impose a projection boundary condition near one of the periodic orbits, while there are several choices for the boundary condition at the other end point; we refer to <cit.> for more details and examples, and to <cit.> for an overview of different methods. In the present context, such computations can be started from connecting orbits, for example, by relaxing the requirement of a connection and sweeping out a relevant part of the respective invariant manifold up to a chosen section.





 §.§ The three-dimensional Poincaré section





We now consider the intersection sets of relevant invariant objects with a suitably chosen section. Recall that, in general, it is not possible to find a so-called global Poincaré section — to which the flow is  transverse everywhere and to which all points return under the flow. Unless the system is periodically forced, any three-dimensional section Σ will have a codimension-one tangency locus C ⊂Σ along which the flow, which we refer to as φ^t, is tangent to Σ.  As a result, a return map f: Σ→Σ can only be defined locally on a subset of Σ∖ C with the additional condition that the images under φ^t also exist  <cit.>. The maximal domain U of definition of a return map given by φ^t is then 


  U := {x  ∈Σ∖C  |  ∃  t > 0  φ^t(x)  ∈Σ∖C }.


Note that both U and f(U) are open subsets of Σ. Any point on the boundary of U either lies on the tangency locus C or, as it is approached from within U, the return time to Σ tends to infinity. Note that the flow direction on either side of C typically changes sign with respect to the normal to Σ; hence, one can split U into two disjoint subsets of Σ, denoted U^+ and U^-, which are separated by C. It is natural to define the return maps f^+ and f^- on these subsets U^+ and U^-, respectively; note that they are generically given by the second return to Σ, that is, by the second iterate f^2 of f. In particular, the intersection points of a periodic orbit with Σ will then be fixed points of f^+ and f^-, which both are then the usual local Poincaré maps near the periodic orbit.

Indeed, the Atri system eq:4DAtri is not periodically forced, and so any three-dimensional section will have a tangency locus. However, the system has the nice property that we can choose a hyperplane that is transverse to the heterodimensional cycle, that is, it is transverse to Γ_1 and Γ_2 as well as to A and B. The geometric idea is to choose a section that contains a joint curve of rotation of both Γ_1 and Γ_2 and, throughout, we consider here the three-dimensional section 


  Σ:= { (c,v, c_t,n) ∈^4 |  c = 0.15 }.


(Note that this is not the section Π used for finding A and B, which was chosen such that it intersects neither Γ_1 nor Γ_2 and separates these two periodic orbits.) Moreover, equation ċ = v of eq:4DAtri implies that the tangency locus of Σ is the two-dimensional plane 


  C := { (c,v, c_t,n) ∈^4 |  c = 0.15 ,  v = 0},


which divides Σ into two open sets, Σ^+ where v > 0 and  Σ^- where v < 0.





 §.§ The PtoP orbit in the section Σ







fig:Pcycle(a) shows the heterodimensional PtoP cycle with the section Σ and the corresponding intersection sets in projection onto (c, c_t, n)-space; note that Σ appears two dimensional in this projection and compare with fig:hetcycle(c). The intersection sets of the PtoP cycle are shown in fig:Pcycle(b) in the three-dimensional section Σ, that is, in (v, c_t,n)-space with c = 0.15. The periodic orbits Γ_1 and Γ_2 intersect Σ^+ and Σ^- in the points γ_1^± and γ_2^±, respectively. The cylinder B intersects Σ in two curves B^+ ⊂Σ^+ from γ_2^+ to γ_1^+ and B^- ⊂Σ^- from γ_2^- to γ_1^-. Finally, the intersection set of the connecting orbit A consists of two sequences


  ( a_k^- )_k ∈ := A ∩Σ^-      ( a_k^+ )_k ∈ := A ∩Σ^+,


with the property that f^±(a_k^±) = a_k+1^±. Notice that, owing to the strong contraction toward Γ_1 in backward time, all a_k^± for k < 0 are extremely close to γ_1^±, respectively.

fig:Pcycle shows that Σ is indeed transverse to the entire heterodimensional cycle. However, this does not imply that all trajectories that make up the manifolds W^u(Γ_1), W^s(Γ_1) W^u(Γ_2) and W^s(Γ_2) are transverse to Σ. As a consequence, and as we will show, the intersection sets of these manifolds with Σ can cross the tangency locus C and may consist of infinitely many disjoint branches. More specifically, since Γ_1 and Γ_2 are transverse to Σ, the local return maps f^± are diffeomorphisms near γ_1^± and γ_2^± and the points are hyperbolic fixed points under f^±, respectively. Hence, the Stable Manifold Theorem <cit.> guarantees the existence and smoothness of unique stable and unstable manifolds for these fixed points. However, their existence as unique smooth manifolds is guaranteed only locally. If we consider them globally, the stable and unstable manifolds of γ_1^± and γ_2^± — that is, the intersection sets of the respective manifolds of Γ_1 and Γ_2 with Σ — may cross the tangency locus C and/or consist of disjoint branches, much as invariant manifolds of endomorphisms <cit.>.

In spite of this fundamental difficulty, we now proceed with computing and visualising how the intersection sets ( a_k^±)_k ∈ of A and B^± of B arise from the interactions between the respective stable and unstable manifolds of the periodic orbits Γ_1 and Γ_2. To this end, we compute and discuss one by one (relevant parts of) two-dimensional manifolds in the four-dimensional phase space of system eq:4DAtri that lead to one-dimensional intersection curves in the three-dimensional section Σ. The final overall picture in the section Σ that we obtain in this way allows us to relate the geometry of objects of the flow to the existing theory on heterodimensional cycles of diffeomorphisms <cit.>, as visualized in fig:hdcycle.





 §.§ The intersection sets ^s, ±(Γ_2) and ^ss, ±(Γ_1)



By definition of the heterodimensional connection A, the points ( a_k^±)_k ∈ = A ∩Σ must lie on the intersection set ^s,±(Γ_2) of W^s(Γ_2) with Σ^±. We can utilize this property to find the curves of ^s,±(Γ_2) locally near points of the set ( a_k^±)_k ∈. For example, the point a^+_0 ∈Σ^+ corresponds to the part of the orbit segment representing A that starts at a^+_0 and ends very close to Γ_2 in the linear approximation of W^s(Γ_2). This part of the orbit segment can be selected and then continued by considering it as a solution of the 2PBVP with one end point in Σ^+ and the other still satisfying the projection boundary condition near Γ_2. 





fig:Ws shows that the continuation of this orbit segment starting from a^+_0 gives a single curve in Σ^+ that connects all the points a_k^+ for k > -1 to γ_2^+. Hence, this curve is (one side of) the local stable manifold of γ_2^+ under f^+, which we also refer to as the primary curve _0^s,+(Γ_2) in ^s,+(Γ_2) ⊂Σ^+. Panel (a) shows, in projection onto (c, c_t, n)-space, how _0^s,+(Γ_2) arises as the intersection of W_0^s,+(Γ_2) with Σ; shown are Γ_1, Γ_2, Σ and the side of W^s(Γ_2) that contains the connecting orbit A (which is not shown). Panel (b) shows how the primary curve _0^s,+(Γ_2) starts at γ_2^+, goes through the points a_k^+ for k > -1 and then passes γ_1^+ very closely, namely through a_-1^+; the latter point is the asterisk (not labelled) that lies practically on the diamond representing γ_1^+; this close passage is explained by the fact that Γ_1 has the very large unstable Floquet multiplier λ^u_1 ≈ 9.6 × 10^4. 

Similarly, starting from a^-_0, we find the primary curve _0^s,-(Γ_2) ⊂^s,-(Γ_2) that contains γ_2^- as well as a_k^- for k > -1 and is, hence, the local stable manifold of γ_2^- under f^-. As fig:Ws shows, both primary curves arise from the fact that the surface W_0^s,-(Γ_2) is a cylinder transverse to Σ. Both curves are entirely contained in their respective half-spaces Σ^±⊂Σ and appear to extend to infinity after passing through a^±_-1, respectively. We remark that starting the continuation from any of the points a^±_k with k ≥ -1 results in the computation of the same two primary curves _0^s,±(Γ_2); notice that we also show, in panel (b), short first pieces of the other side of the primary curves _0^s,±(Γ_2), which also appear to extend to infinity. However, neither side of these curves contain any of the points a^±_k for k ≤ -2.

The λ-lemma <cit.> implies that W^s(Γ_2) accumulates on the two-dimensional strong stable manifold W^ss(Γ_1) of Γ_1. Since the two primary curves _0^s,±(Γ_2) extend to infinity after passing through the points a_-1^±, there must be other curves in the intersection set ^s,±(Γ_2) that go through the points a_k^± for k ≤ -2; these curves must accumulate near the points γ_1^± on the two primary curves  _0^ss,±(Γ_1) through the points γ_1^±. Additional curves in ^s,±(Γ_2) can again be found by continuation of orbit segments in W^s(Γ_2) that now start at the points a_k^± for k ≤ -2; we denote these successive curves by _1^s,±(Γ_2), _2^s,±(Γ_2) and so on.





fig:Wss illustrates the accumulation of W^s(Γ_2) onto W^ss(Γ_1). Panel (a) shows, in projection onto (c, c_t, n)-space, Γ_1, Γ_2 and Σ, the two first pieces W_0^s(Γ_2) and W_1^s(Γ_2) of W^s(Γ_s), and the part of W^ss(Γ_1) that gives the intersection curve _1^s,+(Γ_2). The surface W_1^s(Γ_2) is very close to W^ss(Γ_1), which is illustrated in panel (b) by the corresponding intersection curves in Σ. Due to the strong contraction in backward time, only the curve _1^s,-(Γ_2) through a^-_-2 can be seen. Already the curve _2^s,-(Γ_2) is so close to _0^ss,-(Γ_1) that they are indistinguishable in fig:Wss(b); similarly, _1^s,+(Γ_2) through a^+_-2 and _0^ss,+(Γ_1) are indistinguishable on the scale of this figure.

We computed the curves _j^s,+(Γ_2) and _j^s,-(Γ_2) for j ≥ 1 to very long arclengths, and found that they lie on either side of the tangency locus and extend to infinity on both sides of the respective points in a_k^± for k ≤ -2. Therefore, we conclude that these infinitely many curves in ^s,+(Γ_2), which map to each other under the inverse of f^±, do not connect to form single curve with _0^s,±(Γ_2).





 §.§ The  intersection sets ^u(Γ_1) and ^uu(Γ_2)



By definition of the heterodimensional cycle A, the points in ( a_k^±)_k ∈∈ A ∩Σ must also lie on the intersection of W^u(Γ_1) with Σ. Analogously to the computations in the previous section, we now compute the intersection curves in the intersection set ^u(Γ_1) by continuation of orbit segments that start very close to Γ_1, in the linear approximation of W^u(Γ_1), and end in Σ, starting initially at one of the points in ( a_k^±)_k ∈⊂Σ. 



 

fig:Wu0 shows the primary intersection curve _0(Γ_1), which is a single closed curve with two branches that connect γ^±_1 and γ^±_2; one of these two branches contains the points a^+_0 and a^-_0 and both branches cross the tangency locus C. Panel (a) shows Γ_1 and the part of W^u(Γ_1) that was computed to find the branch of _0(Γ_1) containing a^±_0; starting from either of these points gives the same result. In contrast to fig:Ws(a), which shows the projection onto (c, c_t, n)-space, the projection in fig:Wu0(a) is now onto (c, v, n)-space. Hence, the plane C appears as a line that divides the projection of Σ into the regions Σ^+ and Σ^-; this projection illustrates that both Γ_1 and _0(Γ_1) cross C. Panel (b) shows all of _0(Γ_1) in Σ; the second branch (not containing a^±_0) was found with continuation starting near γ_1. Note that _0(Γ_1) ∩Σ^+ maps to _0(Γ_1) ∩Σ^- under the first-return map f; see also <cit.>. Similar to what we found for _0^s,±(Γ_2), the intersection curve _0(Γ_1) contains all backward images under f^± of the two intersection points a^±_0, that is, a^±_k ∈_0(Γ_1) for all k < 0; note that all of these points (marked by stars)  are again extremely close to γ^±_1.





Since _0(Γ_1) does not contain any of the points a^±_k with k > 0, this means that W^u(Γ_1) must intersect Σ in other curves. In fact, most points on _0(Γ_1) never return to Σ, but there is a small segment on _0(Γ_1) that returns as the closed intersection curve _1(Γ_1) and contains the points a^±_1; this small  segment on _0(Γ_1) contains a^+_0 and maps under the first-return map f to a similarly small segment containing a^-_0.

The curve _1(Γ_1) can be found by continuation from either a^+_0 or a^-_0, and it is shown in fig:Wu1. Panel (a) shows, in projection onto (c, v, n)-space,  Γ_1 and the part of W^u(Γ_1) that generates _1(Γ_1). Panel (b) shows _1(Γ_1) and _0(Γ_1) in Σ; also shown are γ^±_1,  γ^±_2, a^±_0 and a^±_1. Note from fig:Wu1(a) that the closed curve _1(Γ_1) has two intersections points with C, with a^+_1 and a^-_1 on either side, such that _1(Γ_1) ∩Σ^+ maps to _1(Γ_1) ∩Σ^- under the first-return map f. Geometrically, the curve _1(Γ_1) bounds a small “bump” of the surface W^u(Γ_1) that crossed Σ again. This type of intersection between a manifold and a global Poincaré section arises from a minimax quadratic tangency on C <cit.>. 

The intersection curve _1(Γ_1) only contains a^±_1 and no other points in A ∩Σ. Hence, there must exist a further return of W^u(Γ_2) to Σ, corresponding to an even smaller segment on _0(Γ_1) that contains a^+_0, maps to a segment in Σ^- that contains a^-_0, a segment in Σ^+ that contains a^+_1, one in Σ^- that contains a^-_1, and finally, produces a new intersection curve _2(Γ_1) that contains (at least) the point a^+_2. This curve _2(Γ_1) can also be found by continuation starting from the orbit segment that ends on a^+_2.





fig:Wu2 shows the result of this continuation, in one direction, for the same projections as in fig:Wu0 and <ref>. Here, panel (a) shows Γ_1 and the respective part of W^u(Γ_1) that generates the computed part of _2(Γ_1) starting at a^+_2, and panel (b) shows _2(Γ_1) with _0(Γ_1) and _1(Γ_1) in Σ. The curve _2(Γ_1) spirals and crosses C many times, which implies that it contains not only a^±_2 but also the points a^±_k with k > 2. 





To see this better, fig:Wu2zoom presents an enlargement. Its panel (a) shows only _2(Γ_1) and the local part of W^u(Γ_1) that corresponds to this spiralling intersection curve. Furthermore, we also plot one side of the stong unstable manifold W^uu(Γ_2) and its corresponding primary intersection curve _0^uu(Γ_2). The curves _2(Γ_1) and _0^uu(Γ_2) are shown in Σ in fig:Wu2zoom(b). This panel clearly shows that _0^uu(Γ_2) is a single closed curve that crosses C twice to connect the two points γ^±_2. Moreover, the curve _2(Γ_1) accumulates on _0^uu(Γ_2) in a spiralling fashion. In particular, this means that _2(Γ_1) contains all points a^±_k for k ≥ 2, which can be seen to converge to γ^±_2. When viewed locally near γ^±_2, this geometry is a result of the λ-lemma. On the other hand, the fact that all points  a^±_k for k ≥ 2 lie on _2(Γ_1) is a global property related to how W^u(Γ_1) intersects Σ.





 §.§ Interaction between ^s,±(Γ_2) and ^u(Γ_1)








Having accounted for all points in the sequence ( a_k^±)_k ∈ = A ∩Σ as part of both ^s(Γ_2) and (Γ_1), we are now able to clarify how the codimension-one heterodimensional PtoP cycle arises from the intersection W^s(Γ_2) ∩ W^u(Γ_1). fig:WsWu illustrates, in projection onto (c, c_t, n)-space, how the two surfaces W^s(Γ_2) and W^u(Γ_1) intersect in the connecting orbit A from Γ_1 to Γ_2. Also shown are the periodic orbits W^s(Γ_2) and W^u(Γ_1) and (the projection of) Σ; compare with fig:Ws(a) and note the difference in projection in fig:Wu0(a). 





fig:WsWuSelf assembles all relevant intersection sets in Σ that are associated with the heterodimensional PtoP cycle. Shown are two different views featuring the curves in ^s,+(Γ_2) and ^s,-(Γ_2) as they intersect _0(Γ_1), _1(Γ_1) and _2(Γ_1) in the heteroclinic points a^±_k in the intersection set of A with Σ. Also shown are the intersection curves W_0^ss,+(Γ_1) and W_0^ss,-(Γ_1) of the strong stable manifold near γ_1^±, and similarly, the intersection set _0^uu(Γ_2) of the strong unstable manifold near γ_2^±. Finally, fig:WsWuSelf also shows the two intersection curves B^± of the cylinder B with Σ. 

Notice that fig:WsWuSelf show a much larger part of the curve _2(Γ_1), not only the part from fig:Wu2 that spirals onto _0^uu(Γ_2) and contains the points a^+_k for k ≥ 2. The curve _2(Γ_1) keeps spiralling and has infinite arclength, even though it is confined to a bounded region. Note that _2(Γ_1) also accumulates onto the closed curve _0(Γ_1) that contains γ_1^±. Hence, ^u, ±(Γ_1) accumulates on itself, which implies the existence of homoclinic orbits to Γ_1. Notice further, that _2(Γ_1) approaches γ_1^± along the direction of B^±, that is, along the weak stable direction, as expected; see fig:Wu2(b). 





fig:WsWuAllsketch is a topological sketch of the objects shown in fig:WsWuSelf. The tangency locus C is indicated by a line that divides Σ into the two parts Σ^+ and Σ^-. Illustrated is the accumulation of the curves ^s,±(Γ_2) onto W_0^ss,±(Γ_1), and the accumulation of _2(Γ_2) onto _0^uu(Γ_2) and onto _0(Γ_1). 

fig:WsWuAllsketch is designed to make the connection between the heterodimensional PtoP cycle in the four-dimensional continuous-time system eq:4DAtri and a heterodimensional cycle between two fixed points of a diffeomorphism in ^3, as sketched in fig:hdcycle. The first return f on Σ maps the curves in Σ^+ to those in Σ^-; while the second-return maps f^± leave their respective domains Σ^± invariant. We now restrict our attention to the action of f^+ on Σ^+; of course, the respective statements hold equally for f^- on Σ^-. In fig:WsWuSelf and <ref>, the curves W_0^ss,+(Γ_1) and B^+ are tangent at γ_1^+ to the strong stable and the weak stable direction, respectively, while the curve _0(Γ_1) is tangent to the unstable direction. Similarly, at γ_2^+, the curves W_0^uu(Γ_1) and B^+ are tangent to the strong unstable and the weak unstable direction, respectively, while the curve _0^s,+(Γ_1) is tangent to the unstable direction. This agrees exactly with the local situation sketched near γ_1 and γ_2 shown in fig:hdcycle. Hence, the curve B^+ connects γ_1^+ with γ_2^+ in fig:WsWuSelf and <ref> exactly as the curve B connects γ_1 with γ_2 in fig:hdcycle. Moreover, the curves _0^s,+(Γ_2) and _0(Γ_1) in fig:WsWuSelf and <ref>  intersect in the points a^+_0 and a^+_1 and also contains the points a^+_k for k ≥ 2, which is also as sketched in fig:hdcycle for the curves W^s(γ_2) and W^u(γ_1). 

The difference is in the global nature of the interacting manifolds: although it is not shown in fig:hdcycle, this sketch suggests that W^s(γ_2) and W^u(γ_1) are single curves that form a heteroclinic tangle by intersecting in the points a_k, where W^s(γ_2) accumulates on W^ss(γ_1) and W^u(γ_1) accumulates on W^uu(γ_2). This is the logical way of completing the picture for a diffeomorphism in ^3. The diffeomorphism f^+ on Σ^+, on the other hand is not defined along the discontinuity locus C  and, moreover, certain parts of W^s(Γ_2) and W^u(Γ_1) do not intersect Σ^+; the latter results in the intersection sets _0^s,+(Γ_2) ∩Σ^+ and _0(Γ_1) ∩Σ^+ in Σ comprising of infinitely many curves. This reflects the fact that there exists no global section for system eq:4DAtri. fig:WsWuSelf shows that both W^s(Γ_2) and W^u(Γ_1) accumulate on W^ss(Γ_1) and W^uu(Γ_2), respectively, which are locally cylinders. In the section Σ, this manifests itself differently near γ_2^+ and γ_1^+ because these local cylinders are intersected by Σ differently: as two arcs in the case of W^ss(Γ_1) and as a closed curve crossing C in the case of W^uu(Γ_2). Nevertheless, locally near the heterodimensional cycle of f^+ in Σ^+, formed by the point γ_1^+, the sequence ( a_k^+ )_k ∈, the point γ_2^+, and the curve B^+, the manifold structure shown in fig:WsWuSelf and <ref> is, indeed, as sketched in fig:hdcycle. 

Finally, fig:WsWuSelf and <ref> nicely illustrate why further heterodimensional PtoP cycles between Γ_1 and Γ_2 must exist for nearby parameter values. Namely, because of its accumulation onto _0^u, ±(Γ_1), the curve of _2^u, ±(Γ_1)  comes arbitrarily close to ^s, ±(Γ_2). Hence, in any neighbourhood of the point (J^∗, s^∗) there must be infinitely many further codimension-one PtoP connections, created by different intersections between the persisting one-dimensional sets ^s, ±(Γ_2) and ^u(Γ_1). Since the cylinder B is structurally stable near (J^∗, s^∗), these will actually form further heterodimensional cycles, which involve more complicated excursions between Γ_1 and Γ_2. Moreover, the existence of homoclinic orbits to both Γ_1 and Γ_2 (see also <cit.>) implies the existence of many other periodic orbits and, hence, suggesting the existence of heterodimensional cycles involving different pairs of periodic orbits nearby in parameter space. 





§ CONCLUSIONS



We showed how stable and unstable manifolds of two saddle periodic orbits with different indices intersect to form a heterodimensional PtoP cycle of a four-dimensional vector field. The PtoP cycle of the vector field intersects the transverse three-dimensional section Σ twice on either side of a two-dimensional subspace C along which the flow is tangent, which represents the typical situation for flows that are not periodically forced. Namely, the Poincaré map can be defined only locally (as a second return to Σ) near the respective intersection sets of the PtoP cycle on either side of C. We found that the structure of the manifolds agrees locally with that of two saddle fixed points in a three-dimensional diffeomorphism as considered in the literature on heterodimensional cycles. The global nature of the intersection set, on the other hand, is more complicated and involves intersection sets crossing the tangency locus C. Geometrically, the codimension-one part of the PtoP cycle in the four-dimensional phase space arises from the fact that: (i) the two-dimensional unstable manifold of the first periodic orbit accumulates on the two-dimensional strong unstable manifold of the second periodic orbit; and (ii) the two-dimensional stable manifold of the second periodic orbit accumulates on the two-dimensional strong stable manifold of the first periodic orbit. How these manifolds intersect the section Σ near either of the two periodic orbits differs, and this explains the differences between the gobal properties of the respective intersection sets.

For these results we made extensive use of advanced numerical methods based on a two-point boundary value problem set-up to define and compute global (sub)manifolds of interest. This approach is versatile, efficient and, in particular, well suited to deal with large differences between expansion/contraction rates of the periodic orbits. It allowed us to compute and present the intersection sets of the respective manifolds in a three-dimensional section transverse to the two periodic orbits and to both heteroclinic connections that comprise the PtoP cycle. We also showed with three-dimensional projections how these intersection sets arise from the respective manifolds in the four-dimensional phase space. 



Ongoing work concerns the structure of parameter space near the heterodimensional PtoP cycle. The codimension-one PtoP connection itself occurs along a curve in the bifurcation diagram in the (J,s)-plane, while the connecting cylinder is structurally stable. Our computations suggest that there should be further curves of PtoP cycles of the same two periodic orbits nearby, which correspond to different non-transverse intersections of their two-dimensional global manifolds. Finding and continuing these and other secondary hererodimensional PtoP connections and cycles is a considerable task. After all, theory predicts that the existence of heterodimensional PtoP can occur robustly, and the question is how this manifests itself in a concrete system. In particular, robustness of heterodimensional cycles requires that there are infinitely many periodic orbits of increasingly higher periods nearby. Hence, presenting an overall consistent picture in the parameter plane also requires the computation of other bifurcation curves, especially those of saddle-node bifurcations of periodic orbits and of their period-doubling bifurcations. Indeed, each of these curves may be an end point of a curve of PtoP connections. For example, we showed that the curve of the primary codimension-one PtoP connection ends on a curve of period-doubling bifurcation. While the primary codimension-one PtoP connection is then lost, the bifurcating period-doubled periodic orbit has the correct index to form a new heterodimensional cycle with the periodic orbit that does not bifurcate. Hence, this new PtoP connection could be seen as the continuation of the primary one. Moreover, we have evidence of the existence of a period-doubling cascade, suggesting that there may be an associated cascade of heterodimensional PtoP cycles to the successive period-doubled periodic orbits. A detailed investigation of the respective bifurcation sets in the (J,s)-plane near the primary PtoP connection will be presented elsewhere. 

From a more general perspective, our study shows that abstract concepts, introduced to advance the theory of dynamical systems, can be detected and investigated in concrete models, including in ODEs (generally the models of choice) rather than in diffeomorphisms; specifically, we considered the Atri model, which is from neuroscience and describes intracellular calcium dynamics. It is an interesting direction for future research to investigate how the complicated dynamics associated with a heterodimensional PtoP cycle manifests itself in experimental measurements. In case that there is evidence of a heterodimensional PtoP cycle playing an important role in organizing the system behaviour, the conclusion must be that the ODE model for the phenomenon studied requires a dimension of at least four. 





§ ACKNOWLEDGEMENTS


We thank Katsutoshi Shinohara for helpful discussions. BK and HMO were supported by Royal Society of New Zealand Marsden Fund grant 16-UOA-286.






10

Anosov
Anosov D V 1967
Geodesic flows on closed Riemannian manifolds with negative curvature
Proc. Inst. Steklov 90 1–235

atri
Atri A, Amundsen J, Clapham D and Sneyd J 1993
A single-pool model for intracellular calcium oscillations and waves in the Xenopus laevis oocyte
Biophysical Journal 65 1727–1739 

bamon
Bamón R, Kiwi J and Rivera-Letelier J 2006 
Wild Lorenz-like attractors 
Preprint, arXiv:0508045 

beyn90
Beyn W-J 1990
The numerical computation of connecting orbits in dynamical systems
IMA J. Numer. Anal. 10(3) 379–405

beyn94
Beyn W-J 1994
On well-posed problems for connecting orbits in dynamical systems
Chaotic Numerics Contemporary mathematics 172 ed P E Kloeden
(Providence, RI: American Mathematical Society) pp 131–168

Birkhoff17
Birkhoff G D 1917
Dynamical systems with two degrees of freedom
Trans. Amer. Math. Soc.  18 199-–300.

bc2004recurrence
Bonatti C and Crovisier S 2004
Récurrence et généricité
Invent. Math. 158(1) 33–104

bcdw2016what
Bonatti C, Crovisier S Díaz L J and Wilkinson A 2016
What is… a blender?
Notices Amer. Math. Soc. 63(10) 1175–1178 

bd1996blender
Bonatti C and Díaz L J 1996 
Persistent nonhyperbolic transitive diffeomorphisms
Ann. of Math. 2nd Series 143(2) 357–396

Bonatti2008
Bonatti C and Díaz L J 2008
Robust heterodimensional cycles and C^1-generic consequences
J. Inst. Math. Jussieu 7 469–-525

bdv2005book
Bonatti C, Díaz L J and Viana M 2005 
Dynamics Beyond Uniform Hyperbolicity. A Global Geometric and Probabilistic Perspective 
Encyclopaedia Math. Sci. 102
(Berlin Heidelberg: Springer-Verlag)

homcont
Champneys A R, Kuznetsov Yu A and Sandstede B 1996
A numerical toolbox for homoclinic bifurcation analysis
Int. J. Bifurcation and Chaos 6(5) 867–887 

sturman2014
Christov I C, Lueptow R M, Ottino J M and Sturman R 2014
A study in three-dimensional chaotic dynamics: granular flow and transport in a bi-axial spherical tumbler
SIAM J. Appl. Dynam. Syst. 13(2) 901–943

cp2015essential
Crovisier S and Pujals E R 2015
Essential hyperbolicity and homoclinic bifurcations: a dichotomy phenomenon/mechanism for diffeomorphisms.
Invent. Math. 201(2) 385–517

Diaz95
Díaz L J 1995
Robust nonhyperbolic dynamics and heterodimensional cycles. 
Ergod. Theory Dyn. Syst. 15 291-–315

katsutoshi
Díaz L J, Kiriki S and Shinohara K 2014
Blenders in centre unstable Hénon-like families: with an application to heterodimensional bifurcations
Nonlinearity 27(3) 353–378

dieci04
Dieci L and Rebaza J 2004
Point-to-periodic and periodic-to-periodic connections
BIT Numerical Mathematics 44(1) 41–62; with erratum in 44(3) 617–618

AutoOrig
Doedel E J 1981
Auto, a program for the automatic bifurcation analysis of autonomous systems
Congr. Numer. 30 265–384

auto 
Doedel E J 2007
Auto-07P: Continuation and bifurcation software for ordinary differential equations,
with major contributions from Champneys A R, Fairgrieve T F, Kuznetsov Yu A, Oldeman B E, Paffenroth R C, Sandstede B, Wang X J and Zhang C;
available at <http://cmvl.cs.concordia.ca/auto/>

vanVoorn08
Doedel E J, Kooi B W, van Voorn G A K and Kuznetsov Yu A 2008
Continuation of connecting orbits in 3D-ODEs: (I) Point-to-cycle connections
Int. J. Bifurcation and Chaos 18(7) 1889–1903

vanVoorn09
Doedel E J, Kooi B W, van Voorn G A K and Kuznetsov Yu A 2009
Continuation of connecting orbits in 3D-ODEs: (II) Cycle-to-cycle connections 
Int. J. Bifurcation and Chaos 19(1) 159–169 

DKO15
Doedel E J, Krauskopf, B and Osinga, H M 2015
Global organization of phase space in the transition to chaos in the Lorenz system
Nonlinearity 28(11) R113-R139

Dullin95
Dullin H R and Wittek A 1995
Complete Poincaré sections and tangent sets
J. Phys. A 28 7157–-7180

Edelstein
Edelstein-Keshet L 1988
Mathematical Models in Biology Birkäuser Mathematics Series
(Boston: McGraw-Hill)

ErmentroutTerman10
Ermentrout G B and Terman D H 2010
Mathematical Foundations of Neuroscience Interdisciplinary Applied Mathematics 35
(New York: Springer-Verlag)

eko-ijbc2005
England J P, Krauskopf B and Osinga H M 2005
Bifurcations of stable sets in noninvertible planar maps
Int. J. Bifurcation and Chaos 15(3) 891–904

England2005
England J P, Krauskopf B and Osinga H M 2005
Computing one-dimensional global manifolds of Poincaré maps by continuation 
SIAM J. Appl. Dynam. Syst. 4(4) 1008–1041

galias1997chua
Galias Z 1997
Positive topological entropy of Chua's circuit: a computer assisted proof
Int. J. Bifurcation and Chaos 7(2) 331–349

GonchenkoMeiss2006
Gonchenko S V, Meiss J D and Ovsyannikov I I 1006 
Chaotic dynamics of three-dimensional Hénon maps that originate from a homoclinic bifurcation
Regular and Chaotic Dynamics 11(2) (2006) 191–212

GonchenkoTuraev2005
Gonchenko S V, Ovsyannikov I I, Simó C and Turaev D 2005
Three-dimensional Hénon-like map and wild Lorenz-like attractors
Int. J. Bifurcation and Chaos 15(11) 3493–3508

GonchenkoTuraev2009
Gonchenko S V, Shilnikov L P and Turaev D 2009
On global bifurcations in three-dimensional diffeomorphisms leading to wild Lorenz-like attractors
Regular and Chaotic Dynamics 14(1) 137–147

GH
Guckenheimer J and Holmes P 1983
Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields 
(New York: Springer-Verlag)

hko-wild
Hittmeyer S, Krauskopf B and Osinga H M 2013
Interacting global invariant sets in a planar map model of wild chaos
SIAM J. Appl. Dynam. Syst. 12(3) 1280–1329

hkos-blender
Hittmeyer S, Krauskopf B, Osinga H M and Shinohara K 2018
Existence of blenders in a Hénon-like family: geometric insights from invariant manifold computations
Nonlinearity 31(10) R239-R267

ho-ijbc2008
Hobbs C A and Osinga H M 2008
Bifurcations of the global stable set of a planar endomorphism near a cusp singularity
Int. J. Bifurcation and Chaos 18(8) 2207–2222

HomburgSandstede10
Homburg A J and Sandstede B 2010
Homoclinic and heteroclinic bifurcations in vector fields
Handbook of Dynamical Systems III ed B Hasselblatt, H W Broer and F Takens
(North Holland: Elsevier) pp 379–524

JacksonRadunskaya15
Jackson T and Radunskaya A 2015
Applications of Dynamical Systems in Biology and Medicine The IMA Volumes in Mathematics and its Applications 158
(New York: Springer-Verlag)

Just00
Just W and Kantz H 2000,
Some considerations on Poincaré maps for chaotic flows
J. Phys. A 33 163-–170

yorke
Kostelich E J, Kan I, Grebogi C, Ott E and Yorke J A 1997
Unstable dimension variability: a source of nonhyperbolicity in chaotic systems
Physica D 109(1–2) 81–90

redbook
Krauskopf B and Osinga H M 2007
Computing invariant manifolds via the continuation of orbit segments;
Numerical Continuation Methods for Dynamical Systems ed B Krauskopf, H M Osinga and J Galán-Vioque 
(Springer Netherlands)  pp 117–154

surveymans
Krauskopf B, Osinga H M, Doedel E J, Henderson M E, Guckenheimer J, Vladimirsky A, Dellnitz M and Junge O 2005
A survey of methods for computing (un)stable manifolds of vector fields
Int. J. Bifurcation and Chaos 15(3) 763–791

kop-siads2007
Krauskopf B, Osinga H M and Peckham P 2007
Unfolding the cusp-cusp bifurcation of planar endomorphisms 
SIAM J. Appl. Dynam. Syst. 6(2) 403–440

krauskopfriess
Krauskopf B and Rieß T 2008
A Lin's method approach to finding and continuing heteroclinic connections involving periodic orbits
Nonlinearity 21(8) 1655–1690

Kuznetsov
Kuznetsov Yu A 2004
Elements of Applied Bifurcation Theory 3rd edition
(New York: Springer-Verlag)

Lee2008
Lee C M, Collins P J, Krauskopf B and Osinga H M 2008
Tangency bifurcations of global Poincaré maps
SIAM J. Appl. Dynam. Syst. 7(3) 712–754

Lin
Lin X-B 1990
Using Melnikov’s method to solve Shilnikov’s problems
Proc. R. Soc. Edinb. A 116(3–4) 295–325

Oldeman03
Oldeman B E, Champneys A R and Krauskopf B 2003
Homoclinic branch switching: a numerical implementation of Lin's method.
Int. J. Bifurcation and Chaos 13(10) 2977–2999

pal1991homoclinic
Palis J 1991
Homoclinic bifurcations, sensitive-chaotic dynamics and strange attractors;
Dynamical systems and related topics (Nagoya, 1990) Advanced Series in Dynamical Systems 9 ed K Shiraiwa 
(Singapore: World Scientific Publishing) pp 466–473

PdM
Palis J and de Melo W 1982 
Geometric Theory of Dynamical Systems 
(New York: Springer-Verlag)

pt1993hyperbolicity
Palis J and Takens F 1993
Hyperbolicity & Sensitive Chaotic Dynamics at Homoclinic Bifurcations
(Cambridge: Cambridge University Press)

rademacher05
Rademacher J D M 2005
Homoclinic orbits near heteroclinic cycles with one equilibrium and one periodic orbit
J. Diff. Eqns. 218(2) 390–443

rademacher10
Rademacher J D M 2010
Lyapunov-Schmidt reduction for unfolding heteroclinic networks of equilibria and periodic orbits with tangencies
J. Diff. Eqns. 249(2) 305–348

rre2003explicit
Roeder R K W, Rapoport B I and Evans T E 2003
Explicit calculations of homoclinic tangles in tokamaks
Phys. Plasmas 10(9) 3796–3799

Smale
Smale S 1967
Differentiable dynamical systems
Bull. Amer. Math.  Soc. 73 747–817

turaev1998
Turaev D and Shilnikov L P 1998
An example of a wild strange attractor 
Sbornik: Mathematics 189(2) 291–314

viana2011fractal
Viana R L, Da Silva E C, Kroetz T, Caldas I L, Roberto M and Sanjuán M A F 2011
Fractal structures in nonlinear plasma physics
Philos. Trans. Royal Soc. A 369(1935) 371–395

Zhang06
Zhang W B 2006
Discrete Dynamical Systems, Bifurcations and Chaos in Economics Mathematics in Science and Engineering 204
(North Holland: Elsevier)

ZKK2012
Zhang W, Krauskopf B and Kirk V 2012 
How to find a codimension-one heteroclinic cycle between two periodic orbits
Discrete Contin. Dyn. Syst. Ser. A 32(8) 2825–2851





