Resonant Tori, Transport Barriers, and Chaos in a Vector Field with a Neimark-Sacker Bifurcation

TL;DR

This paper numerically investigates a three-dimensional dissipative vector field undergoing a Neimark-Sacker bifurcation, analyzing the invariant torus, manifold embeddings, and their roles in transport barriers and global bifurcations.

Contribution

It provides a detailed numerical analysis of the invariant torus and manifold structures in a system with a Neimark-Sacker bifurcation, highlighting their influence on global dynamics.

Findings

Identification of the invariant torus and its parameter-dependent behavior.

Analysis of stable and unstable manifold embeddings as transport barriers.

Insights into the role of manifolds in global bifurcations.

Abstract

We make a detailed numerical study of a three dimensional dissipative vector field derived from the normal form for a cusp-Hopf bifurcation. The vector field exhibits a Neimark-Sacker bifurcation giving rise to an attracting invariant torus. Our main goals are to (A) follow the torus via parameter continuation from its appearance to its disappearance, studying its dynamics between these events, and to (B) study the embeddings of the stable/unstable manifolds of the hyperbolic equilibrium solutions over this parameter range, focusing on their role as transport barriers and their participation in global bifurcations. Taken together the results highlight the main features of the global dynamics of the system.

Figures

Equations

$$f(x,y,z)=\left(\begin{array}{c}(z-\beta )x-\delta y\\ \delta x+(z-\beta )y\\ \tau +\alpha z-\frac{z^3}{3}-(x^2+y^2)(1+\epsilon z)+\zeta z{x}^3\end{array}\right),$$

$$z^{\prime }=\tau +\alpha z-\frac{z^3}{3}=\text{mathchar}58g(z).$$

$$\text{real}({\lambda }_1^s)\le \cdots \le \text{real}({\lambda }_{k_s}^s)<0<\text{real}({\lambda }_1^u)\le \cdots \le \text{real}({\lambda }_{k_u}^u).$$

$${\lambda }_1^s{\theta }_1\frac{\partial }{\partial {\theta }_1}P({\theta }_1,\dots ,{\theta }_{k_s})+\dots +{\lambda }_{k_s}^s{\theta }_{k_s}\frac{\partial }{\partial {\theta }_{k_s}}P({\theta }_1,\dots ,{\theta }_{k_s})=f(P({\theta }_1,\dots ,{\theta }_{k_s})),$$

$$P(0,\dots ,0)=\hat{p},\text{mbox}and\frac{\partial }{\partial {\theta }_j}P={\xi }_j^s,1\le j\le k_s.$$

$$\Lambda_{s}=\left(\begin{array}[]{ccc}\lambda_{1}&\ldots&0\\ \vdots&\ddots&\vdots\\ 0&\ldots&\lambda_{k_{s}}\end{array}\right),$$

$$DP(\theta ){\Lambda }_s\theta =f(P(\theta )).$$

$$L(\theta ,t)=e^{{\Lambda }_{s}t}\theta .$$

$$\varphi \left(P(\theta ),t\right)=P\left(e^{{\Lambda }_{s}t}\theta \right),$$

$$P({\theta }_1,{\theta }_2)=\sum _{m=0}^{\infty }\sum _{n=0}^{\infty }\left(\begin{array}{c}{p}_{mn}^1\\ p_{mn}^2\\ p_{mn}^3\end{array}\right){\theta }_1^m{\theta }_2^n,$$

$${\lambda }_1{\theta }_1\frac{\partial }{\partial {\theta }_1}P({\theta }_1,{\theta }_2)+{\lambda }_2{\theta }_2\frac{\partial }{\partial {\theta }_2}P({\theta }_1,{\theta }_2)=f(P({\theta }_1,{\theta }_2)),$$

$${\lambda }_1{\theta }_1\frac{\partial }{\partial {\theta }_1}P({\theta }_1,{\theta }_2)+{\lambda }_2{\theta }_2\frac{\partial }{\partial {\theta }_2}P({\theta }_1,{\theta }_2)=\sum _{m=0}^{\infty }\sum _{n=0}^{\infty }(m{\lambda }_1+n{\lambda }_2)p_{mn}{\theta }_1^m{\theta }_2^n,$$

$$P_j({\theta }_1,{\theta }_2)=\sum _{m=0}^{\infty }\sum _{n=0}^{\infty }{p}_{mn}^j{\theta }_1^m{\theta }_2^n,$$

$$(p^3\ast p^1{)}_{mn}=\sum _{j=0}^m\sum _{k=0}^n{p}_{(m-j)(n-k)}^3{p}_{jk}^1,$$

$$(p^3\ast p^1\ast p^1\ast p^1{)}_{mn}=\sum _{i_1=0}^m\sum _{i_2=0}^{i_1}\sum _{i_3=0}^{i_2}\sum _{k_1=0}^n\sum _{k_2=0}^{k_1}\sum _{k_3=0}^{k_2}{p}_{(m-i_1)(n-k_1)}^3{p}_{(i_1-i_2)(k_1-k_2)}^1{p}_{(i_2-i_3)(k_2-k_3)}^1{p}_{i_3{k}_3}^1.$$

$$(m\lambda_{1}+n\lambda_{1})\left[\begin{array}[]{c}p_{mn}^{1}\\ p_{mn}^{2}\\ p_{mn}^{3}\end{array}\right]=$$

$${\tiny\left[\begin{array}[]{c}(p^{3}*p^{1})_{mn}-\beta p_{mn}^{1}-\delta p_{mn}^{2}\\ (p^{3}*p^{2})_{mn}-\beta p_{mn}^{2}+\delta p_{mn}^{1}\\ \alpha p_{mn}^{3}-\displaystyle\frac{1}{3}(p^{3}*p^{3}*p^{3})_{mn}-(p^{1}*p^{1})_{mn}-(p^{2}*p^{2})_{mn}-\varepsilon(p^{1}*p^{1}*p^{3})_{mn}-\varepsilon(p^{2}*p^{2}*p^{3})_{mn}+\zeta(p^{1}*p^{3}*p^{3}*p^{3})_{mn}\end{array}\right],}$$

$$(p^3\ast p^1{)}_{mn}=p_{00}^3{p}_{mn}^1+p_{00}^1{p}_{mn}^3+(p^3\hat{\ast }{p}^1{)}_{mn},$$

$$(p^3\hat{\ast }{p}^1{)}_{mn}=\sum _{j=0}^m\sum _{k=0}^n{\hat{\delta }}_{jk}^{mn}{p}_{(m-j)(n-k)}^3{p}_{jk}^1,$$

$${\hat{\delta }}_{jk}^{mn}=\{\begin{array}{ll}0& \text{mbox}ifj=k=0\\ 0& \text{mbox}ifj=m\text{mbox}andk=n\\ 1& \text{mbox}otherwise\end{array}.$$

$$g(x,z)=xz,$$

$$(g\circ P)_{mn}=\nabla g(p_{00}^{1},p_{00}^{3})\left[\begin{array}[]{c}p_{mn}^{1}\\ p_{mn}^{3}\end{array}\right]+(p^{1}\hat{*}p^{3})_{mn}.$$

$$(m\lambda_{1}+n\lambda_{2})p_{mn}^{1}=\nabla g(p_{00}^{1},p_{00}^{3})\left[\begin{array}[]{c}p_{mn}^{1}\\ p_{mn}^{3}\end{array}\right]+(p^{1}\hat{*}p^{3})_{mn}-\beta p_{mn}^{1}-\delta p_{mn}^{2}.$$

$$\nabla g(p_{00}^{1},p_{00}^{3})\left[\begin{array}[]{c}p_{mn}^{1}\\ p_{mn}^{3}\end{array}\right]-\beta p_{mn}^{1}-\delta p_{mn}^{2}-(m\lambda_{1}+n\lambda_{2})p_{mn}^{1}=-(p^{1}\hat{*}p^{3})_{mn},$$

$$\nabla g(p_{00}^{1},p_{00}^{3})\left[\begin{array}[]{c}p_{mn}^{1}\\ p_{mn}^{3}\end{array}\right]-\beta p_{mn}^{1}-\delta p_{mn}^{2}=\nabla f_{1}(\hat{p})\left[\begin{array}[]{c}p_{mn}^{1}\\ p_{mn}^{3}\end{array}\right].$$

$$Df(\hat{p})=\left[\begin{array}[]{c}\nabla f_{1}(\hat{p})\\ \nabla f_{2}(\hat{p})\\ \nabla f_{3}(\hat{p})\end{array}\right],$$

$$(f\circ P{)}_{mn}=Df(\hat{p})p_{mn}+$$

$${\tiny\left[\begin{array}[]{c}(p^{3}\hat{*}p^{1})_{mn}\\ (p^{3}\hat{*}p^{2})_{mn}\\ -\displaystyle\frac{1}{3}(p^{3}\hat{*}p^{3}\hat{*}p^{3})_{mn}-(p^{1}\hat{*}p^{1})_{mn}-(p^{2}\hat{*}p^{2})_{mn}-\varepsilon(p^{1}\hat{*}p^{1}\hat{*}p^{3})_{mn}-\varepsilon(p^{2}\hat{*}p^{2}\hat{*}p^{3})_{m}n+\zeta(p^{1}\hat{*}p^{3}\hat{*}p^{3}\hat{*}p^{3})_{mn}\end{array}\right].}$$

$$(m{\lambda }_1+n{\lambda }_2)p_{mn}=Df(\hat{p})p_{mn}+$$

$${\tiny\left[\begin{array}[]{c}(p^{3}\hat{*}p^{1})_{mn}\\ (p^{3}\hat{*}p^{2})_{mn}\\ -\displaystyle\frac{1}{3}(p^{3}\hat{*}p^{3}\hat{*}p^{3})_{mn}-(p^{1}\hat{*}p^{1})_{mn}-(p^{2}\hat{*}p^{2})_{mn}-\varepsilon(p^{1}\hat{*}p^{1}\hat{*}p^{3})_{mn}-\varepsilon(p^{2}\hat{*}p^{2}\hat{*}p^{3})_{m}n+\zeta(p^{1}\hat{*}p^{3}\hat{*}p^{3}\hat{*}p^{3})_{mn}\end{array}\right],}$$

Full text

Keywords: Invariant tori, Neimark-Sacker bifurcation, parameterization method

AMS Subject Classifications: 34C45, 37G35, 37M05, 37C55

1 Introduction

Interactions between equilibrium and oscillating states provide a basic mechanism for generating complicated dynamics in nonlinear systems. Such interactions are the focus of the present investigation, where we study the global dynamics of a one parameter family of three dimensional vector fields whose main features are stable and saddle type equilibrium solutions and a periodic orbit with a complex conjugate pair of Floquet exponents. The frequency of the periodic orbit together with the frequency of the complex exponent constitute two competing natural modes of oscillation. Tension between these internal frequencies gives rise to a number of interesting dynamical phenomena. In particular the system admits a Neimark-Sacker bifurcation, where the real part of the complex conjugate Floquet exponents crosses the imaginary axis as the parameter is changed [1, 2]. The loss of stability of the periodic orbit triggers the appearance of a smooth attracting invariant torus supporting quasiperiodic motions. Global bifurcations of the torus lead to resonant motions and eventually to the appearance of a chaotic attractor.

The local theory describing the appearance, evolution, and disappearance of invariant tori in dissipative multi-frequency systems is well developed and we refer the reader to the works of [3, 4, 5, 6, 7, 8, 9, 10] on dissipative dynamics, the related work of [11, 12, 13, 14] on area and volume preserving systems, and to the numerical studies of [5, 15, 16, 17, 18, 19, 20]. Global questions about the dynamics of dissipative systems with attracting tori lead to difficult analytical and computational problems. While many important theoretical questions have by now been settled – see for example [21, 22, 23, 24] and the references therein – there remains much to be learned from careful qualitative studies of important special cases.

While many of the canonical examples of dynamical systems theory come from specific physical or engineering applications, another source of compelling problems is to study the normal form of an interesting bifurcation. Such systems caricature the universal features of an entire class of problems, and this is precisely the setting of the present paper. We study, from the numerical point of view, a model derived from the normal form unfolding the cusp-Hopf bifurcation. This system, which is described in detail in Section 1.1, was first introduced in [20] and is referred to hereafter as the Langford system. As already mentioned in the opening paragraph, a main feature is that model undergoes a supercritical Neimark-Sacker bifurcation resulting in the appearance of a smooth attracting invariant torus. We provide detailed computations of the torus, monitoring it as its dynamics change from quasi-periodic to resonant – and as it changes from a $C^{k}$ to a $C^{0}$ invariant manifold – before finally breaking up in a global bifurcation resulting in the appearance of a chaotic attractor.

In addition to undertaking a detailed description of the attracting invariant torus, the present work aims also to describe the dynamics nearby. We are especially interested in any dramatic changes in the organization of the phase space as the bifurcation parameter is varied. Such changes may be triggered by either local or global bifurcations. More precisely we have the following distinction.

Definition 1.1**.**

We say that a bifurcation is *local *if it occurs due to a change in linear stability of an invariant object.

Definition 1.2**.**

We say that a bifurcation is global if it is triggered by the formation of tangencies between invariant manifolds.

In the present work we mainly observe local bifurcations of equilibrium and periodic solutions – and global bifurcations where the invariant manifolds do not intersect at all prior to, and intersect transversally after the global bifurcation.

The discussion just presented makes it clear that the goals of the present work require careful examination of the embeddings of some hyperbolic invariant objects like stable/unstable manifolds of equilibrium and periodic orbits. Much of the analysis is simplified by considering an appropriate surface of section, as this reduces the invariant torus and the stable/unstable manifolds of periodic orbits to one dimensional curves. Embeddings of stable/unstable manifolds attached to equilibrium solutions on the other hand are often difficult to characterize in a fixed section, and studying their structure is more delicate. We employ the parameterization method of [25, 26, 27] to compute high order representations of the two dimensional local stable/unstable manifolds in the full three dimensional phase space. The parameterization method is a functional analytic framework for studying invariant manifolds and in particular provides a natural notion of a-posteriori error analysis. The local representations obtained using the parameterization method are extended using standard adaptive numerical integration schemes.

The detailed numerical calculations performed in the main body of the paper provide insights into the dynamics of the system which are summarized in Section 5, and which give a coarse qualitative description of the global dynamics as a function of the bifurcation parameter. Since the Langford system is derived from a normal form, it is reasonable to expect qualitatively similar dynamics in an appropriately restricted region for any system undergoing the sequence of bifurcations unfolded by this vector field. Moreover, the approach of using the parameterization method in conjunction with geometric analysis in Poincaré sections could be applied to the study of a wide variety of dynamical systems.

The remainder of the paper is organized as follows. In the next two subsections we first describe the three dimensional model under consideration, and then discuss briefly some related literature. In Section 2 we review the main ideas of the parameterization method for an equilibrium solution and apply them to the Langford system. In Section 3 we study the Neimark-Sacker bifurcation and the resulting attracting invariant torus in an appropriate Poincaré section. We provide numerical evidence for a global bifurcation from a quasi-periodic torus to a resonant one, and for a second global bifurcation which destroys the torus and appears to create a chaotic attractor. In Section 4 we study the invariant manifolds of the equilibrium solutions before and after the Neimark-Sacker bifurcation, with an emphasis on the omega limit sets of two dimensional unstable manifolds and on the role of the manifolds as separatrices. We also study their role in further global bifurcations. We conclude the paper in Section 5 with a summary of our observations about the global dynamics of the system and a few further conclusions and observations.

1.1 The Langford system

We study the dynamical system generated by the 3D vector field

[TABLE]

where $\varepsilon=0.25$, $\tau=0.6$, $\delta=3.5$, $\beta=0.7$, $\zeta=0.1$, and with $\alpha>0$ treated as a bifurcation parameter. The system was derived by Langford in [20] by truncating to second order the normal form unfolding a simultaneous Hopf/cusp bifurcation. A third order term is then added to the vector field, breaking the axial symmetry of the second order truncation. This symmetry breaking is important for describing the dynamics following a generic bifurcation. Since the Hopf bifurcation creates a periodic orbit, and the cusp bifurcation creates three nearby equilibrium solutions, interesting interactions between these states are to be expected.

We begin with some elementary observations which inform the numerical study to follow. Note that the $z$-axis is an invariant sub-system as $x=y=0$ implies that $x^{\prime}=y^{\prime}=0$. The dynamics on the $z$-axis are governed by the scalar differential equation

[TABLE]

The function $g(z)$ is illustrated in Figure 1, and since $\tau>0$, $g$ has one, two, or three zeros depending on the parameter $\alpha$. Moreover, equilibria of $f$ occur at $(0,0,z_{*})$ where $z_{*}$ is a zero of $g$. Observe that for large positive $z$, $z^{\prime}<0$. While for large negative $z$, $z^{\prime}>0$. That is, the field tends to diminish the $z$ value of a phase point whose $z$ value happens to be large.

For all $\alpha\in\mathbb{R}$ Equation (1) has at exactly one equilibrium solution with $x=y=0$ and $z>0$, which we denote by $p_{0}\in\mathbb{R}^{3}$. This equilibrium has one stable eigenvalue, whose eigenvector coincides with the $z$-axis. The remaining eigenvalues are complex conjugate unstable. At $\alpha\approx 0.9321697517861$ there is a saddle node bifurcation giving rise to a new pair of equilibrium points $p_{1},p_{2}\in\mathbb{R}^{3}$. These equilibria persist for all larger values of $\alpha$. One of the equilibrium points appearing out of the saddle node bifurcation is fully stable, with three eigenvalues having negative real parts, and we denote it by $p_{2}\in\mathbb{R}^{3}$. The other new equilibrium, which we denote by $p_{1}\in\mathbb{R}^{3}$, is a saddle-focus with a complex conjugate pair of stable eigenvalues and one real unstable eigenvalue. The unstable eigenvector again coincides with the $z$-axis. Indeed, since the $z$-axis is invariant, the stable manifold of $p_{0}$ and the unstable manifold of $p_{1}$ coincide, and are contained in the $z$-axis. This intersection is not transverse, and is rather forced by a rotational symmetry of the problem.

Now consider the plane $z=\beta$, and note that when the field is projected onto this plane the nonlinear terms vanish from the first two components giving a pure rotation. The plane is however not invariant, as $z^{\prime}$ does not vanish there. Nevertheless there is a periodic orbit $\gamma$ near the $z=\beta$ plane. This periodic orbit, and the invariant $z$ axis organize the dynamics of the system. The vector field along with the periodic orbit and the dynamics on the $z$-axis are illustrated in the left frame of Figure 2.

As we will see below, the periodic orbit $\gamma$ has a pair of complex conjugate Floquet exponents, hence solutions of the differential equation tend to circulate around $\gamma$. The orbit may be either attracting or repelling depending on the value of $\alpha$. This circulation about the periodic orbit is a dominant feature of the dynamics.

Further insight into the dynamics is obtained by numerically integrating some trajectories (phase space sampling), as was done in the work of Langford [20]. We provide, for the sake of completeness, the results of a few such simulations. The results illustrated in Figure 3 make clear the typical behavior of the system, and suggest the existence of a “torus-like” attractor. Simulations were run for roughly one hundred time units. The periodic orbit $\gamma$ runs through the center of the torus but is, as we will see, repelling for these parameter values. The saddle focus points $p_{0}$ and $p_{1}$ are at the top and bottom of the torus.

1.2 Some remarks on the literature

Roughly speaking, the dynamics described above suggests the system as a toy model for dissipative vortex dynamics, or for a rotating viscus fluid. There is a rich literature on the dynamics of vortex bubbles, and the interested reader might consult the works of [14, 28, 29, 30, 31] for a more thorough discussion of the literature. We remark that the torus bifurcations seen in the Langford system are similar to those seen in the piecewise linear electronic circuit of [32], the commodity distribution model of [33], and the mechanical oscillators of [34, 35] to name only a few. The appearance and destruction of invariant tori, as well as resonance phenomena and routes to chaos are discussed much more generally in [5, 36] and the references found therein.

One further remark is in order. The system given by Equation (1) has been called the Aizawa system by some researchers, and is the subject of some other recent work on visualization. For example researchers interested in computer animation [37], three dimensional printing [38], and even in graphical arts [39] have made interesting studies and use this name for the equations. This nomenclature seems to be a misnomer, as the equations do not appear in the works of Yoji Aizawa, and a more appropriate name for Equation (1) would seem to be the Langford system, due to the fact that – as already mentioned above – the system was proposed in [20].

2 Review of the parameterization method

The parameterization method is a general functional analytic framework for studying invariant manifolds of discrete and continuous time dynamical systems, first developed in [25, 26, 27] in the context of stable/unstable manifolds attached to fixed points of nonlinear mappings on Banach spaces, and later extended in [40, 41, 42] for studying whiskered tori. There is a thriving literature devoted to computational applications of the parameterization method, and the interested reader may want to consult [10, 18, 19, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54], though the list is far from being exhaustive. A much more complete discussion is found in the book [55].

This section provides a practical overview of the parameterization method with a strong emphasis on numerical aspects utilized in the sequel. The discussion focuses on analytic vector fields, and requisite formal series calculations are carried out for the specific example of Equation (1). Since this material is not completely standard outside a certain circle of practitioners, it is included primarily so that the present work stands alone for a broad readership. The reader either already familiar with or uninterested in these developments is encouraged to skip ahead to Section 3, referring back to this section only as needed.

The Langford system admits equilibrium solutions with complex eigenvalues, so that it is best to present the entire theory for complex vector fields. Later we explain how to recover parameterizations of real invariant manifolds associated with complex conjugate eigenvalues of real vector fields. So, let $f\mathrel{\mathop{\mathchar 58\relax}}\mathbb{C}^{k}\to\mathbb{C}^{k}$ be an analytic vector field and $\hat{p}\in\mathbb{C}^{k}$ have $f(\hat{p})=0$ so that $x(t)=\hat{p}$ is an equilibrium solution of the differential equation $x^{\prime}=f(x)$. Assume for the sake of simplicity that $Df(\hat{p})$ is diagonalizable over $\mathbb{C}$ having $k_{s}$ stable (and $k_{u}$ unstable) eigenvalues of multiplicity one. We do not necessarily assume that $k_{s}+k_{u}=k$, that is we do not rule out the possibility of some center directions at $\hat{p}$ (though this situation will not occur in the present work).

Label the stable eigenvalues as $\lambda_{1}^{s},\cdots,\lambda_{k_{s}}^{s}$ and the unstable ones as $\lambda_{1}^{u},\cdots,\lambda_{k_{u}}^{u}$ and order them according to the convention that

[TABLE]

Since $Df(\hat{p})$ is diagonalizable there are linearly independent eigenvectors $\xi_{1}^{u},\cdots,\xi_{k_{u}}^{u}\in\mathbb{C}^{k}$ and $\xi_{1}^{s},\cdots,\xi_{k_{s}}^{s}\in\mathbb{C}^{k}$ associated with the unstable and stable eigenvalues respectively.

Remark 2.1**.**

The assumption that $Df(\hat{p})$ is diagonalizable is made only for the sake of convenience. See [25] for a much more general theoretical setup. See also [47] for a complete description of the functional analytic set up and examples of the numerical implementation when there are repeated eigenvalues. Nevertheless, the assumption holds in the examples considered throughout the present work.

2.1 Invariance equation

Given the setup introduced in the previous section we are interested in computing an accurate representation of the $k_{s}$ dimensional local stable manifold attached to $\hat{p}$. The parameterization method seeks a smooth surjective map $P$ satisfying the first order system of partial differential equations

[TABLE]

for $\theta=(\theta_{1},\ldots,\theta_{k_{s}})\in\mathbb{C}^{k_{s}}$, and subject to the first order constraints

[TABLE]

Equation (2) is referred to as the invariance equation for $P$. A map $P$ solving Equation (2) subject to the first order constraints of Equation (3) is a parameterization of the local stable manifold, as we explain below. Making the obvious adjustments for the unstable eigenvalues/eigenvectors leads to a parameterization method for the $k_{u}$ dimensional unstable manifold.

To explain the meaning of Equation (2), let

[TABLE]

so that the invariance equation becomes

[TABLE]

In the language of differential geometry, this equation says that the push forward by $P$ of the linear vector field $\Lambda_{s}\theta$ is equal to the vector field $f$ restricted to the image of $P$. Where the vector fields are equal they generate the same dynamics. But the dynamics generated by $\Lambda_{s}\theta$ are completely understood: all orbits converge exponentially to the origin. It follows that all orbits on the image of $P$ converge to $\hat{p}$. Since the image of $P$ is a smooth $k_{s}$ dimensional disk, it is a local stable manifold for $\hat{p}$. The situation is illustrated in Figure 4.

The observation is made more precise as follows. Denote by $\phi\colon\mathbb{C}^{k}\times\mathbb{C}\to\mathbb{C}^{k}$ the flow generated by $f$. The flow generated by $\Lambda_{s}\theta$ is given explicitly by

[TABLE]

One checks that $P$ satisfies Equation (2) if and only if

[TABLE]

for all $t\geq 0$. This flow conjugacy is illustrated in Figure 5. Elementary proofs of these claims are found in any of the references [48, 55, 56]. Moreover, replacing the stable by the unstable eigenvalues and eigenvectors in the discussion above and reversing time, the entire discussion carries through for the unstable manifold.

Remark 2.2** (Real analytic vector fields and manifolds).**

If $f\colon\mathbb{R}^{k}\to\mathbb{R}^{k}$ is a real analytic vector field with a real equilibrium $\hat{p}\in\mathbb{R}^{k}$ then the discussion above applies to an analytic extension of the vector field in a neighborhood of $\hat{p}$. In this case any complex eigenvalues of $Df(\hat{p})$ appear in complex conjugate pair, and the associated eigenvectors can be taken complex conjugate. We look for a solution $P$ of Equation (2) taking real values on complex conjugate variables. This condition imposes a symmetry on the Taylor coefficients of the parameterization $P$, as illustrated explicitly in the examples below.

2.2 Formal series solution of Equation (2) for the Langford system

In this section we further restrict to the case of interest in the present work, where $\hat{p}\in\mathbb{C}^{3}$ and $\lambda_{1},\lambda_{2}\in\mathbb{C}$ are a pair of stable (or unstable) complex conjugate eigenvalues and $\lambda_{3}$ has the opposite stability. Let $\xi_{1},\xi_{2}\in\mathbb{C}^{3}$ be an associated pair of linearly independent complex conjugate eigenvectors. Since the field is analytic, we look for an analytic parameterization

[TABLE]

satisfying Equation (2), which in this case is reduced to

[TABLE]

where $f\colon\mathbb{C}^{3}\to\mathbb{C}^{3}$ is the Langford vector field given in Equation (1). Here $p_{mn}^{j}\in\mathbb{C}$ for all $j=1,2,3$. Imposing the linear constraints of Equation (3) gives that $p_{00}=\hat{p}$, $p_{10}=\xi_{1}$ and $p_{01}=\xi_{2}$.

Now we would like to expand Equation (2) in terms of the power series. The left hand side of Equation (2) is

[TABLE]

on the level of power series. To expand the right hand side we begin by writing

[TABLE]

for $j=1,2,3$ to denote the component power series. The field contains the nonlinear terms $zx$, $zy$, $z^{3}$, $x^{2}z$, $y^{2}z$, and $zx^{3}$ (see again Equation (1)). Computing the power series for $f\circ P$ requires expanding these monomials of components of $P$, which is accomplished using Cauchy products. For example the coefficients of $P_{3}\cdot P_{1}$ are

[TABLE]

while the coefficients of $P_{3}\cdot P_{1}^{3}$ are

[TABLE]

Other products are similar.

Substituting these power series expansions into the invariance equation (2) and matching like powers of $\theta_{1}$ and $\theta_{2}$ leads to

[TABLE]

[TABLE]

for $m+n\geq 2$. To isolate terms of order $(m,n)$ consider that

[TABLE]

where

[TABLE]

and

[TABLE]

The point here is that $(p^{3}\hat{*}q^{1})_{mn}$ is precisely the sum left when terms containing $p_{mn}$ are extracted from the Cauchy product.

This expression is directly related to the derivative of $f$. To see this, let

[TABLE]

and note that Equation (6) becomes

[TABLE]

Using this notation the first component of Equation (5) is

[TABLE]

Isolating terms of order $(m,n)$ on the left and lower order terms on the right gives

[TABLE]

which is linear in $p_{mn}^{1}$. Comparing the right hand side in the equation above with the vector field $f$, and recalling that $\hat{p}=p_{00}$, we see that

[TABLE]

Combining the equation above with a nearly identical computation for the second component, and a somewhat lengthier computation for the third component, and noting that

[TABLE]

we obtain the expansion

[TABLE]

[TABLE]

Substituting this expansion into Equation (5) gives

[TABLE]

[TABLE]

and by isolating terms of order $(m,n)$ on the left we obtain the linear homological equations

[TABLE]

for $p_{mn}$, where

[TABLE]

with

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

We make the following observations:

• •

While our derivation of Equation (7) is particular to the Langford system of Equation (1), we remark that the form of the homological equations is always the same. Only the right hand side depends on the particular nonlinearity of the given system.

• •

The matrix acting on $p_{mn}$ is the characteristic matrix for the differential at $\hat{p}$. Then the equation is uniquely solvable at order $(m,n)$ if $m\lambda_{1}+n\lambda_{2}$ is not an eigenvalue.

• •

Since $\lambda_{3}$ has the opposite stability of $\lambda_{1},\lambda_{2}$, we obtain the non-resonance condition

[TABLE]

If the non-resonance conditions are satisfied for all $m,n\in\mathbb{N}$ with $m+n\geq 2$, then the formal series solution of Equation (2) is formally well defined to all orders.

• •

If $\lambda_{2}=\overline{\lambda_{1}}$, that is if we consider the complex conjugate case, then there is no possibility of a resonance and we can compute the power series coefficients of the parameterization to any desired finite order.

• •

When $\lambda_{1},\lambda_{2}$ are complex conjugates, the coefficients of $P$ have the symmetry $\overline{p_{nm}}=p_{mn}$ for all $m+n\geq 2$. This is seen by taking complex conjugates of both sides of the homological equation, and using the fact that $Df(\hat{p})$ is a real matrix.

Since $\hat{p}$ is real, choosing complex conjugate eigenvectors $\xi_{2}=\overline{\xi_{1}}$ enforces the symmetry to all orders. The power series solution $P$ has complex coefficients, but we obtain the real image of $P$ by taking complex conjugate variables. That is, we define for the real parameters $s_{1},s_{2}$ the function

[TABLE]

which parameterizes the real stable/unstable manifold.

2.3 Numerical considerations

The homological equations derived in the previous section allow us to recursively compute the power series coefficients of the stable/unstable manifold parameterization $P$ to any desired order $m+n=N$. The coefficients are uniquely determined up to the choice of the scaling of the eigenvectors. In practical applications we have to decide how to answer the following questions:

• •

To what order $N$ will we compute the approximate parameterization?

• •

What scale to choose for the eigenvectors?

• •

On what domain do we to restrict the polynomial $P^{N}$?

In practice we proceed as follows. First we choose a convenient value for $N$, based on how long we want to let the computations run. Then, we always restrict $P$ to the unit disk for the sake of numerical stability. Finally, we choose the eigenvector scaling so that the last coefficients, the coefficients of order $N$, are smaller than some prescribed tolerance. A good empirical rule of thumb is that the truncation error is roughly the same magnitude as the $N$-th order coefficients.

In practice we can prescribe the size of the $N$-th order terms as soon as we know the exponential decay rate of the coefficients. In the next section we describe the relationship between the scale of the eigenvectors and the exponential decay rate.

2.3.1 Rescaling the eigenvectors

In Section 2.2 we saw that the power series coefficients of the parameterization are uniquely determined up to the choice of the eigenvector. Since the eigenvectors are unique up to the choice of length, we have that the length determines uniquely the coefficients. In fact the effect of rescaling the eigenvectors is made completely explicit as follows. The material in this section is discussed in greater detail in [46].

Suppose that

[TABLE]

is the formal solution of Equation (2), with

[TABLE]

where $\|\xi_{1}\|=\|\xi_{2}\|=1$. Assuming that $P$ is bounded and analytic on the complex poly-disk with radii $R_{1},R_{2}>0$, there is a $C>0$ so that

[TABLE]

by the Cauchy estimates.

Now choose non-zero $s_{1},s_{2}\in\mathbb{R}$ and define the rescaled eigenvectors

[TABLE]

The new parameterization associated with the rescaled eigenvectors is given by

[TABLE]

where

[TABLE]

See [46] for a proof of this identity, and also the discussion in [25, 27]. The coefficients of the rescaled parameterization have the new exponential decay rate given by

[TABLE]

These observations lead to a practical algorithm. First compute the parameterization $P$ with an arbitrary choice of eigenvector scaling (for example scaled to have length one). Then solve the homological equations to some order $N_{0}$ using this scaling, and compute $C$, $R_{1}$ and $R_{2}$ using an exponential best fit. Suppose that $\varepsilon_{0}>0$ is the desired tolerance, that is the desired size of the order $N\geq N_{0}$ coefficients. Then we choose $s_{1}$ and $s_{2}$ so that

[TABLE]

Finally we recompute the coefficients $q_{mn}$ for $2\leq m+n\leq N$. The rescaled coefficients could be computed from the old coefficients using the formula of Equation (8). In practice however better results are obtained by recomputing the coefficients $q_{mn}$ from scratch via the homological equations.

We remark that in the case of complex conjugate eigenvalues we want the eigenvectors to be complex conjugates. Assuming that $\xi_{2}=\overline{\xi}_{1}$ we take $s_{1}=s_{2}\in\mathbb{R}$ so that $\eta_{2}=\overline{\eta}_{1}$. Also note that by choosing our domain to be the unit poly-disk, we have that $R_{1}=R_{2}=1$, further simplifying the analysis.

2.3.2 A-posteriori error

Once we have chosen the polynomial order $N$ and the scaling of the eigenvectors, that is once we have uniquely specified our parameterization to order $N$, we would like a convenient measure of the truncation error. As mentioned above, a good heuristic indicator is that the error is roughly the size of the highest order coefficients (assuming we take the unit disk as the domain of our approximate parameterization). In this section we discuss a more quantitative indicator.

We remark that there exist methods of a-posteriori error analysis for the parameterization method, which – when taken to their logical conclusion – lead to mathematically rigorous computer assisted error bounds on the truncation errors. The interested reader will find fuller discussion and more references to the literature in [27, 46, 47, 55, 56, 57] and discussion of related techniques in [58, 59, 60].

The analysis in the present work is qualitative and we don’t require the full power of mathematically rigorous error bounds. Instead we employ an error indicator inspired by the fact that the parameterization satisfies the flow invariance property given in Equation (4). We choose $T{\not=}0$, and a partition of the interval $[0,2\pi]$ into $K$ angles, $\alpha_{j}=2\pi j/(K+1)$, for $0\leq j\leq K$. Since we are interested in the case of complex conjugate eigenvalues $\lambda,\overline{\lambda}\in\mathbb{C}$, we define complex conjugate parameters

[TABLE]

and the linear mapping

[TABLE]

which maps complex conjugate inputs to complex conjugate outputs. The a-posteriori indicator is

[TABLE]

Here $T>0$ if the complex conjugate eigenvalues $\lambda,\overline{\lambda}$ are stable and $T<0$ if they are unstable. In practice the flow map $\phi(x,t)$ will be evaluated using a numerical integration scheme, and the accuracy of the indicator is limited by the accuracy of the integrator.

2.3.3 A numerical example

As an example of the performance of the method, consider the parameterization of the two dimensional unstable manifold of the Langford system (Equation (1)) at the equilibirum $p_{0}$, computed to order $N=20$. Figure 6 illustrates the effect of the choice of the eigenvector scaling on the decay rate of the Taylor coefficient. We remark that the magnitude of the last Taylor coefficient computed is a good heuristic indicator of the size of the truncation error. For example if we choose eigenvectors scaled to length one, we obtain the decay rate illustrated in the left frame of Figure 6, and we see that the norm of the largest coefficient of order twenty is about $10^{-6}$. On the other hand if we rescale to eigenvector to have length $1/2$ then the coefficients decay as in the right frame of Figure 6, and the largest norm of any coefficient of order twenty is now about $10^{-12}$.

To visualize the parameterized local manifold we evaluate the polynomial approximation on the unit disk. First we take a Delaunay triangulation of the unit disk as illustrated in the left frame of Figure 7. This triangulation of the unit disk is pushed forward to the phase space $\mathbb{R}^{3}$ by the polynomial parameterization, resulting in a triangulation of the two dimensional local unstable manifold as illustrated in the top left frame of Figure 8.

To “grow” a larger representation of the unstable manifold we choose a fundamental domain, for example by taking $\tau=0.25$ and considering the annulus in parameter space formed by the boundary of the unit disk and by the circle of radius $R=\left|e^{\lambda_{u}\tau}\right|\approx 0.733$. We mesh this annulus using $100$ angular subdivisions and $40$ radial subdivisions, as illustrated in the right frame of Figure 7. We lift this fundamental domain to the phase space and repeatedly apply the time $\tau=0.25$ map via numerical integration of the vertices of the triangulation. We refine the mesh whenever any side of a triangle in phase space gets too large. In the present work we measure “too large” just by looking at the resulting picture.

The top right, and bottom frames of Figure 8 illustrate the results of iterating a triangulation of a fundamental domain for the local unstable manifold at $p_{0}$, and we see that the “bubble” grows in a quite regular way. However, by the time we take $60$ iterates the embedding of the initial annulus is becoming quite complicated.

2.4 Numerical approximation methods for stable/unstable manifolds

The literature devoted to numerical approximation of stable/unstable manifolds is substantial, and we take a moment to reframe the techniques just discussed in this light. A classic general reference is the work of [61]. The essential remark is that computational methods for studying stable/unstable manifolds decompose naturally into two independent tasks:

• •

Step 1: Calculate an approximation of the local invariant manifold.

• •

Step 2: Advect the local approximation, “growing” the representation of the manifold.

A natural approach to Step 1 is to approximate the local manifold to first order, reducing the problem to linear algebra. That is, by computing the eigenvalues/eigenvectors of the differential at the equilibrium we can approximate the local stable/unstable manifolds by the stable/unstable eigenspaces. Step 2 is in general much more difficult, due to the fact that nonlinearities cause the manifold to grow in a highly nonuniform way. For this reason, much work focuses on the development of powerful methods for Step 2. We refer the interested reader to the works of [62, 63, 64, 65, 66, 66, 67, 68, 69], and also to the survey paper of [70] for much fuller discussion of the topic. The woks just cited develop sophisticated adaptive subdivision schemes to control the accuracy and complexity of the advection problem, growing the stable/unstable manifolds in a uniform way.

Another way to fight the nonuniformity encountered at Step 2 is to employ a higher order approximation scheme, and hence to compute a larger portion of the stable/unstable manifold at Step 1. The idea is that a manifold approximation holding in a large neighborhood of the equilibrium reduces the dramatic expansion which results from integrating a very small polygonal manifold patch until it describes a large portion of the manifold.

The parameterization method as discussed in the present section accomplishes this. Indeed, deriving the homological equations for the system facilitates the implementation of programs which compute the Taylor coefficients of the local parameterization to any desired order. We refer back to the calculations discussed in Section 2.3.3 where we saw the parameterized local manifold grow quite uniformly after the initial high order computation. The error from computing the local invariant manifold from Step 1 can be estimated even in a large neighborhood of the equilibrium using the a-posteriori indicator.

Of course, even when the parameterization method is used in Step 1, we have to employ advection schemes to see a larger portion of the manifold. The parameterization method simply provides a high order approach to Step 1: it does not eliminate the need for Step 2. In fact any of the Step 2 schemes mentioned above could be used in conjunction with a high order approximation computed at Step 1 using the parameterization method. Far from being competitors, the various techniques complement one another. See [71, 72, 73] for examples of calculations which combine the parameterization method in Step 1 with adaptive advection schemes in Step 2.

3 The invariant torus

The first of our two main goals is to study the appearance of the smooth attracting invariant torus, the major changes in its dynamics as the bifurcation parameter increases – including interestingly enough the loss of differentiability – and finally the disappearance of the torus in a global bifurcation resulting in the appearance of a new chaotic attractor. The discussion takes place in a Poincaré section, where periodic orbits are reduced to collections of points, their stable/unstable manifolds are reduced to curves, and invariant tori are reduced to invariant circles.

3.1 Neimark-Sacker bifurcation in the return map

We begin by studying the dynamics near the periodic orbit $\gamma$ as the bifurcation parameter $\alpha$ varies. To this end we fix as a surface of section the half plane $\Sigma$ given by $x=0,y>0$ (with $z$ free) and consider the first return map $R\colon\Sigma\to\Sigma$, which is well defined in a (possibly quite large) neighborhood of the periodic orbit $\gamma$. In the discussion that follows all fixed points and $k$-cycles of $R$ are computed using standard Newton schemes, and derivatives of the Poincaré map are computed by integrating the variational equations of the flow.

We first observe that for $0<\alpha\leq 0.65$ the first return map has an attracting fixed point $p_{*}\in\Sigma$ corresponding to the attracting periodic orbit $\gamma$ discussed in the introduction. At $\alpha_{1}\approx 0.697144898322973$ the fixed point looses stability, triggering a super-critical Neimark-Sacker bifurcation (see [74] for precise definitions). This results in the appearance of a smooth attracting invariant circle $\Gamma$ near $p_{*}$ in $\Sigma$, which is of course an invariant Torus $\mathcal{T}$ for the flow. The bifurcation value is computed using a Newton scheme for an appropriate augmented system where the parameter $\alpha$ is treated as one of the unknowns, hence the bifurcation parameter is known to roughly machine precision. Such techniques are discussed at length in the classic works of [75, 76, 77].

The dynamics in the section just at and just after the bifurcation are illustrated in Figure 9. For $\alpha>\alpha_{1}$ the fixed point $p_{*}$ is repelling and the invariant circle is attracting. The general theory of the Neimark-Sacker bifurcation [74] dictates that for small enough $\epsilon>0$, the invariant circle at $\alpha=\alpha_{1}+\epsilon$ is smoothly conjugate to an irrational rotation. The four frames of Figure 9 illustrate the initially attracting fixed point (top left frame), the Neimark-Sacker bifurcation (top right frame), and the attracting invariant circle surrounding the now repelling fixed point where the size of the circle gets larger as $\alpha$ increases (bottom two frames).

3.2 Resonant tori

When $\alpha>0$ is large enough the dynamics on the invariant circle $\Gamma$ change in a fundamental way, as we now discuss. We say that $q$ is a period $k$ point for $R\colon\mathbb{R}^{2}\to\mathbb{R}^{2}$ if $q,R(q),R^{2}(q),\ldots,R^{k-1}(q)$ is a collection of $k$ distinct points having $R^{k}(q)=q$. We say that the point set

[TABLE]

is a $k$-cycle for $R$. Notions like stability and stable/unstable manifolds of $k$-cycles are defined in the obvious way after observing that $q$ and each of its iterates are fixed points of the composition map $R^{k}$. See for example [50] for precise definitions and references to the literature. The following notion is critical in the discussion to follow.

Definition 3.1**.**

Let $R\colon\mathbb{R}^{2}\to\mathbb{R}^{2}$ be a smooth map of the plane and $\Gamma\subset\mathbb{R}^{2}$ be a topological circle invariant under $R$. We say that $\Gamma$ is a simple resonant invariant circle if there is an attracting $k$-cycle $\mathcal{Q}_{1}$ and a saddle $k$-cycle $\mathcal{Q}_{2}$ so that

[TABLE]

The situation is that the one dimensional unstable manifold of the saddle cycle is completely absorbed into the basin of attraction of the stable cycle, in such a way that a circle is formed. In this case the dynamics on the invariant circle are conjugate to a gradient system. Observe that the unstable manifold of the saddle cycle is smooth (analytic if the map is), even if – as we will see below – the regularity of the invariant circle is another matter completly. The situation is illustrated in Figure 10 for the simple case of a one-cycle.

Remark 3.2**.**

If the decomposition of $\Gamma$ requires multiple stable and saddle cycles of different periods we say that we have a compound resonant invariant circle. However we do not encounter this situation in the present study, and for this reason we usually drop the term “simple” and say simply that we have a resonant invariant circle.

Remark 3.3**.**

Suppose that $R$ is a Poincaré map for a 3-dimensional smooth vector field $f$, and that $R$ has a simple resonant invariant circle $\Gamma$. Then the flow generated by $f$ has a resonant invariant torus given by

[TABLE]

We further remark that the global regularity of a resonant invariant circle (or torus) is determined by the linearization of $R$ at $q_{1}$ or any of its iterates. So for example if $DR(q_{1})$ has real distinct stable eigenvalues then the resulting invariant circle is finitely differentiable, with regularity determined by the ratio of these eigenvalues. If on the other hand $DR(q_{1})$ has complex conjugate eigenvalues then the torus in phase space is only $C^{0}$, as the unstable manifold of $\mathcal{Q}_{2}$ is forced to approach $\mathcal{Q}_{1}$ in a spiraling fashion and the resulting curve cannot be differentiable or even Lipschitz at $q_{1}$ or any of its iterates.

3.3 Resonant tori in the Langford system

The formation of a resonant invariant torus in the Langford system of Equation (1) involves a global bifurcation which can be observed in the Poincaré section, as we now describe. We begin with the observation that at $\alpha=0.82$ there is an attracting $3$-cycle, which we denote by $\mathcal{Q}_{1}$, and which lies outside the invariant circle $\Gamma$. The basin of attraction of the $3$-cycle is fairly small, as there is a saddle type $3$-cycle nearby, which we denote by $\mathcal{Q}_{2}$.

For parameter values near $\alpha=0.82$, the unstable manifold of $\mathcal{Q}_{2}$ has the following behavior: half of $W^{u}(\mathcal{Q}_{2})$ accumulates on the attracting invariant circle $\Gamma$ while the other half accumulates to the attracting $3$-cycle $\mathcal{Q}_{1}$. Things remain much the same for nearby parameter values, for example at $\alpha=0.8224$, with the caveat that the saddle $3$-cycle $\mathcal{Q}_{2}$ has moved closer to $\Gamma$. The situation is illustrated in Figure 11. The stable manifold of $\mathcal{Q}_{2}$ appears to form a separatrix between the basins of attraction of $\Gamma$ and $\mathcal{Q}_{1}$. See for example the left frame of Figure 11.

For some parameter value $0.8224<\alpha_{2}<0.825$ there appears to be a global bifurcation where $\mathcal{Q}_{2}$ collides with the invariant circle $\Gamma$. At this point $W^{u}(\mathcal{Q}_{2})$, rather than accumulating on the invariant circle $\Gamma$, has become $\Gamma$. Both halves of the unstable manifold accumulate at $\mathcal{Q}_{1}$, which is now inside the invariant circle as well. See Figure 12 for an illustration of the phase space configuration just before and just after the global bifurcation. The situation remains for parameter values $\alpha>\alpha_{2}$ as illustrated in Figure 13.

Let $q_{1}$ and $q_{2}$ be points on the $3$ cycles $\mathcal{Q}_{1}$ and $\mathcal{Q}_{2}$ respectively. By numerical calculation we find that the eigenvalues of $DR^{3}(q_{1})$ are complex conjugate stable. Then the resonant invariant circle appearing in this bifurcation is only $C^{0}$. This is a dramatic change, as for $\alpha<\alpha_{2}$ the simulations indicate that the torus is smooth (at least finitely differentiable). The global bifurcation just described gives a vivid natural example of a low regularity invariant manifold for a smooth (in fact analytic) vector field.

3.4 Transient chaotic motions

Increasing the bifurcation parameter $\alpha$ past the global bifurcation at $\alpha_{2}$ shows that the embedding of the attracting $C^{0}$ resonant invariant circle $\Gamma$ appears to get even “wilder”. See the three frames of Figure 14. The blue curve illustrates the unstable manifold of $\mathcal{Q}_{2}$ and in the left two frames is contained in the invariant circle $\Gamma$, indicating that the circle is loosing regularity.

A more quantitative discussion about dynamical complexity in the system begins by observing that just before the global bifurcation at $\alpha_{2}$, as $\mathcal{Q}_{2}$ is approaching $\Gamma$, there is the appearance of chaotic dynamics in the system. To see this we observe that in the top left corner of the left hand frame of Figure 12, the unstable manifold of $\mathcal{Q}_{2}$ (vertical blue curve) is to the right side of the stable manifold (vertical red curve). In the top left corner of the right hand frame in the same figure we see that the situation is reversed: the unstable manifold now being to the left of the stable one (again vertical blue and red respectively). Since the curves move continuously this suggest that there should be a range of parameters where they intersect.

Figure 15 shows that this is indeed the case. At $\alpha=0.8225$ we can see that the resonant torus in phase space has not formed yet, as the unstable manifold (blue curve) does not accumulate at the stable $3$-cycle (green point). Here, on close inspection we see that $W^{s}(\mathcal{Q}_{2})$ and $W^{u}(\mathcal{Q}_{2})$ do intersect, in fact transversally. Then there is a Smale horse shoe and hence chaotic dynamics in a neighborhood of $\Gamma$ [78]. Note however that the invariant circle $\Gamma$ is still attracting and that the dynamics on $\Gamma$ are simple before and after the global bifurcation at $\alpha_{2}$. This suggests that the chaotic motions are only transient, in the sense that the horse shoe is not in the attractor.

3.5 Destruction of the torus and appearance of a chaotic attractor

By further increasing the bifurcation parameter we eventually observe the destruction of the invariant torus, as we now describe. We begin by recalling a classical result concerning the disappearance of an invariant circle. We refer to the discussion in [36] for the details of the proof and generalizations to higher dimensions. See also [79] and [5, 15, 16]. The set up is as follows.

Suppose that a one parameter family of smooth discrete time dynamical systems $R_{\alpha}\colon\mathbb{R}^{2}\to\mathbb{R}^{2}$ has at $\alpha=\alpha_{0}$ a resonant $C^{0}$ invariant circle $\Gamma\subset\mathbb{R}^{2}$ formed by the closure of the unstable manifold of a saddle cycle $\mathcal{Q}_{2}$ accumulating to a stable cycle $\mathcal{Q}_{1}$ as in Figure 10. We then have a resonance region in the sense of [36] and, while it may be obvious it is nevertheless important to note that, the resonant torus $\Gamma$ is robust under small perturbations including small changes in $\alpha$. This is because the saddle cycle, the attracting cycle, and any local unstable manifold of the saddle are structurally stable objects. Since a local unstable manifold is all that is required to reach the basin of attraction of $\mathcal{Q}_{1}$, we have that the resonant torus is robust. It follows that there is a one parameter family of attracting invariant circles $\Gamma(\alpha)$ for $\alpha$ near $\alpha_{0}$. Indeed the tori vary continuously in $\alpha$, again see [36].

Now, suppose that at some parameter $\hat{\alpha}>\alpha_{0}$ the invariant torus no longer exists. Then, by the least upper bound property of $\mathbb{R}$ there is an $\alpha^{*}$ so that $\Gamma$ is a continuous attracting invariant circle on the interval $[\alpha_{0},\alpha^{*})$, but that for $\alpha>\alpha^{*}$ this fails to be true. The theorem on torus breakdown [36] asserts the following three possible mechanisms for the destruction of $\Gamma$: (i) loss of cycle stability– i.e. a local bifurcation at $\mathcal{Q}_{1}$, (ii) occurrence of a tangency bifurcation of the stable and saddle cycles on $\Gamma$, or (iii) formation of a homoclinic tangency between $W^{u}(\mathcal{Q}_{2})$ and $W^{s}(\mathcal{Q}_{2})$. Mechanism (iii) is illustrated schematically in Figure 17.

Of course the theorem gives only a trichotomy. It does not say which alternative actually occurs in a given example, and it is with this in mind that we investigate the fate of the invariant torus in the Langford system. The situation is illustrated in the two frames of Figure 17, where we observe the formation of a homoclinic tangle for $W^{u,s}(\mathcal{Q}_{2})$. Since these manifolds do not intersect in the frame on the left, and do intersect in the frame on the right, they must develop a tangency at some point $0.92<\alpha_{3}<0.93$. While the closure of the unstable manifold (blue) on the left is still an attracting invariant circle, in the right frame $W^{u}(\mathcal{Q}_{2})$ is no longer contained in the attractor suggesting that the torus was destroyed in the homoclinic tangency – that is that we have in alternative (iii).

Further numerical evidence for this claim is given in the three frames of Figure 19. Here the bifurcation parameter is increased slightly further to $\alpha=0.95$ so that the picture opens up a little. The frame on the right is obtained by iterating a large number of initial conditions until they numerically converge to the attractor, represented by the black curve. Note that the stable and saddle $3$-cycles (green and red collections of dots) appear to have moved off the attractor as they do not touch the black curve. This indicates multi-stability in the system as the green orbit is itself attracting. Moreover the left frame shows the numerically computed unstable manifold $W^{u}(\mathcal{Q}_{2})$, and it is clear by comparing the left with the right frame that – while $W^{u}(\mathcal{Q}_{2})$ is accumulating on the chaotic attractor union $\mathcal{Q}_{1}$, the unstable manifold of $\mathcal{Q}_{2}$ is no longer contained in the attractor. This becomes even more clear when we look again at the right frame of the figure and see no black line from the red dots to the green: hence the attractor is no longer a resonant torus.

Zooming in suggests in fact that the attractor is now a quite complicated shape, as illustrated in the middle frame of Figure 19, and also in Figure 18. Indeed Figure 18 reveals the fractal structure of the unstable manifold of $\mathcal{Q}_{2}$. The actual structure of the set is even more complicated than any picture can reveal, as results from topological dynamics imply that once there is a transverse homoclinic for $\mathcal{Q}_{2}$ the closure of the unstable manifold, which contains the attractor, is an indecomposable continuum [80, 81]. We refer also the the work of [82] on the persistence of normally hyperbolic invariant manifolds in the absence of uniform rates and to the much more recent and constructive work of [83].

3.6 Visualizing the torus in phase space

Studying the dynamics of a three dimensional system in a two dimensional Poincaré section helps us to locate bifurcations of the invariant tori by reducing them to invariant circles. Nevertheless, it is still desirable to visualize dynamical structures of the original system in the full phase space, and to this end we provide several images which show side by side the invariant objects found in the Poincaré section and the corresponding invariant objects for the full Langford system.

See for example Figure 20. The left frame illustrates the attracting invariant circle in the section for the parameter value $\alpha=0.85$. The three red and green dots represent the saddle type and stable $3$-cycles $\mathcal{Q}_{2}$ and $\mathcal{Q}_{1}$ respectively. The magenta dot in the center of the frame represents the repelling fixed point, while the blue curve is the unstable manifold of the saddle. The blue curve is clearly absorbed into the basin of attraction of the stable $3$-cycle, forming a resonant invariant circle. The $3$-cycles $\mathcal{Q}_{1},\mathcal{Q}_{2}$ give rise to periodic solutions in $\mathbb{R}^{3}$, which we denote by $\gamma_{1}$ and $\gamma_{2}$ respectively.

The right frame of the same figure illustrates the embedding in phase space of the same objects. Here the red curve represents the repelling periodic orbit $\gamma$, green curve the attracting periodic orbit $\gamma_{1}$, and the blue curve is the saddle periodic orbit $\gamma_{2}$. The unstable manifold of $\gamma_{2}$ accumulates at $\gamma_{1}$ forming the resonant torus. Half of the torus is cut-away so that the skeleton given by the periodic orbits stands our clearly.

Figures 21 and 22 depict the same information at $\alpha=0.9$ and $\alpha=0.929$. Taken together the three images provide much more insight into the structure of the invariant dynamics than can be gained by studying simulations of individual orbits like those illustrated in Figure 3.

4 Global dynamics

We now come to the second part of the present work, and study the dynamics not on but near the attracting invariant torus. The remaining important features of the surrounding phase space are the equilibrium solutions on the $z-$ axis, and their invariant manifolds. For this part of the study we abandon the Poincaré section and consider features of the three dimensional phase space.

4.1 Accumulation of $W^{u}(p_{0})$ on a component of the global attractor

One of the most important features in the phase space of the Langford system is the equilibrium point $p_{0}\in\mathbb{R}^{3}$, which for all $\alpha\geq 0$ is located on the positive $z$ axis and which has two dimensional unstable manifold and one dimensional stable manifold. For all $\alpha\geq 0$ the stable manifold of $p_{0}$ is a subset of the $z$-axis. The two dimensional unstable manifold on the other hand is much more interesting. All the calculations in this section are preformed using the parameterization method/continuation scheme as discussed in the Section 2.

Recall that for $0\leq\alpha\leq\alpha_{4}\approx 0.9321697517861$, the point $p_{0}$ is the only equilibrium of the Langford system. The manifold $W^{u}(p_{0})$ is illustrated in Figure 23 for six such values of $\alpha$. We see that for $\alpha<\alpha_{1}$ – the value of the Neimark- Sacker bifurcation – $W^{u}(p_{0})$ accumulates on the attracting periodic orbit $\gamma$ as seen in frame (a) of Figure 23. The periodic orbit $\gamma$ appears to be the global attractor in this parameter range.

After the Neimark- Sacker bifurcation at $\alpha_{1}\approx 0.697144898322973$ the periodic orbit $\gamma$ is repelling and $W^{u}(p_{0})$ appears to accumulate on the smooth attracting invariant torus $\mathcal{T}$ which was discussed at length in Section 3. This is seen in frames (b), (c) and (d) of Figure 23. Frames (d) and (e) illustrate the situation after the appearance of the attracting period three cycle in the Poincaré section (see Figure 11), and there is an attracting periodic orbit in phase space which we denote by $\gamma_{1}$. The system is bistable, with the attractor being the union of the invariant torus $\mathcal{T}$ and the periodic orbit $\gamma_{1}$. The manifold $W^{u}(p_{0})$ accumulates on the union of these two objects – a disjoint set – and this is what introduces the rough folds in the embedding seen in Frames (d) and (e).

Frame (f) illustrates the embedding of $W^{u}(p_{0})$ for $\alpha>\alpha_{2}$ but just before the global bifurcation at $\alpha=\alpha_{4}$ which destroys the torus. Here the torus is resonant and only $C^{0}$, a fact which is clearly reflected in the embedding of $W^{u}(p_{0})$.

For $\alpha>\alpha_{4}$ we are past the saddle node bifurcation, and $p_{0}$ is no longer the unique equilibrium. This has dramatic consequences for the global dynamics, and we illustrate the unstable manifold for two such parameter values in Figure 24. Somewhere between $\alpha=0.95$ (frame (a) of the figure) and $\alpha=1.1022$ (frame (b)) something dramatic happens. The change however is not easily understood by looking only at $W^{u}(p_{0})$, and we must consider the embedding of new invariant objects which appear only after the occurrence of the saddle node bifurcation.

4.2 $W^{s}(p_{1})$ as a separatrix

At $\alpha_{4}\approx 0.9321697517861$ the system undergoes a saddle node bifurcation, resulting in the appearance of two new equilibrium solutions denoted $p_{1}$ and $p_{2}$. For all $\alpha>\alpha_{4}$ the point $p_{2}$ is a stable equilibrium, making $p_{2}$ a new component of the global attractor. The equilibrium solution at $p_{1}$ on the other hand is a saddle, with one dimensional unstable manifold on the $z$-axis and a two dimensional stable manifold associated with a pair of complex conjugate eigenvalues.

Our simulations suggest that for some range of $\alpha>\alpha_{4}$, the two dimensional invariant manifold $W^{s}(p_{1})$ is a separatrix for the basin of attraction of $p_{2}$ and the attractor onto which $W^{u}(p_{0})$ accumulates. In this sense, $W^{s}(p_{1})$ forms a kind of “bubble”, where inside the bubble we have an attractor comprised of either the resonant torus $\mathcal{T}$ or the chaotic set which appears after the destruction of $\mathcal{T}$. The inside of the bubble is the basin of attraction of this attractor, and the outside of the bubble is the basin of the stable equilibrium $p_{2}$. The situation is illustrated in Figure 25.

4.3 Heteroclinic intersections and the loss of bistability

Studying $W^{u}(p_{0})$ and $W^{s}(p_{1})$ reveals yet another global bifurcation which dramatically alters the phase space dynamics of the system. It appears that for some $\alpha_{5}\approx 1.05$ these manifolds develop a tangency, and that after this tangency there are transverse heteroclinic connections from $p_{0}$ to $p_{1}$. The situation is illustrated in Figures 26 and 27, where we see the transverse intersections of the manifolds and the resulting heteroclinic connections respectively.

Once intersections appear between $W^{u}(p_{0})$ and $W^{s}(p_{1})$, the latter ceases to function as a separatrix, and orbits can pass from inside the bubble to outside. The equilibrium $p_{2}$ remains attracting and its basin appears now to extend into the inside of the bubble. Indeed our numerical simulations suggest that for $\alpha>1$, that is after the formation of intersections between the unstable/stable manifolds of $p_{0}$ and $p_{1}$, the equilibrium $p_{2}$ becomes the global attractor for the system. That is, all orbits which start inside the bubble eventually accumulate there. This finally explains the “openness” of Figure 24 (b) remarked upon earlier: the occurrence of the heteroclinic tangency between $W^{u}(p_{0})$ and $W^{s}(p_{1})$ appears to destroy the attractor which had previously dominated the dynamics inside the bubble.

5 Conclusions and Discussion

We will summarize the results of the present work by sketching the main features of the global dynamics of the Langford system (1) as suggested by our analysis. First recall the main local and global bifurcations undergone by the system.

• •

At $\alpha_{1}\approx 0.697144898322973$ the periodic orbit $\gamma$ undergoes a Neimark-Sacker bifurcation. This is a local bifurcation of $\gamma$ which results in the appearance of the invariant torus $\mathcal{T}$.

• •

At $\alpha_{2}\approx 0.823$ the invariant torus $\mathcal{T}$ develops a resonance. After this $\mathcal{T}$ is the union of two periodic orbits $\gamma_{2}$ (saddle stability), and $\gamma_{1}$ (attracting), and the unstable manifold of $\gamma_{2}$. The resonance is triggered by the collision of a saddle periodic orbit with the invariant torus. Since the torus is an attractor before and after the collision, this bifurcation involves no change in the stability of $\mathcal{T}$ and is hence a global bifurcation.

• •

At $\alpha_{3}\approx 0.925$ there is a global bifurcation triggered by the formation of a tangency between $W^{s}(\gamma_{1})$ and $W^{u}(\gamma_{1})$.

• •

At $\alpha_{4}\approx 0.9321697517861$ there is a saddle node bifurcation resulting in the appearance of the equilibrium points $p_{1}$ (saddle-focus stability) and $p_{2}$ (attracting). This is a local fold bifurcation for the equilibrium points.

• •

At $\alpha_{5}\approx 1$ a global bifurcation is triggered by the development of a tangency between $W^{u}(p_{0})$ and $W^{s}(p_{1})$.

Between the bifurcation values listed above, we conjecture based on our numerical studies that the system has the following properties.

Conjecture 5.1** (Sketch of the global dynamics).**

The flow generated by the Langford vector field, given in Equation (1), has that;

For $0<\alpha<\alpha_{1}$ the periodic orbit $\gamma$ is the global attractor. 2. 2.

For $\alpha_{1}<\alpha<\alpha_{2}$ the global attractor is either $\mathcal{T}$ or $\mathcal{T}\cup\gamma_{1}$. 3. 3.

For $\alpha_{2}<\alpha<\alpha_{3}$ the global attractor is $\mathcal{T}$. 4. 4.

For $\alpha_{3}<\alpha<\alpha_{4}$ the global attractor is $\mathcal{T}\cup p_{2}$. 5. 5.

For $\alpha_{4}<\alpha<\alpha_{5}$ there is multi-stability. The global attractor is comprised of at least the components $\tilde{\mathcal{T}}$ - the chaotic attractor appearing after the break-up of the invariant torus, the attracting periodic orbit $\gamma_{2}$, and the attracting equilibrium solution $p_{2}$. 6. 6.

For $\alpha>\alpha_{5}$ the equilibrium solution $p_{2}$ is the global attractor. 7. 7.

For $0<\alpha<\alpha_{3}$ the unstable manifold $W^{u}(p_{0})$ accumulates on the global attractor, which is either $\gamma$ (until $\alpha=\alpha_{1}$) or $\mathcal{T}$. 8. 8.

For $\alpha_{3}<\alpha<\alpha_{5}$ the stable manifold $W^{s}(p_{1})$ is a separatrix. The basin of attraction of $p_{2}$ is outside the bubble formed by $W^{s}(p_{1})$.

It is essential to stress that the eight points above are still just conjectures, however well informed. The numerical work carried out in the present work is not sufficient to rule out other components of the global attractor, for example other attracting periodic orbits near the resonant torus or the chaotic attractor. This point is elaborated on below.

It is also worth remarking that softer sorts of conclusions are encapsulated in the paper’s many figures. The deliberate calculations and three dimensional renderings of invariant manifolds throughout our work provide more delicate insights into the dynamics of the system than are obtained by straight forward simulations of ensembles of initial conditions. As a final illustration of this point we give in Figure 28 a side by side comparison of the results obtained using the methodology of the present work with the results obtained by direct integration, for the parameter value $\alpha=0.9321$. Simulation results cannot illuminate the full attractor, as numerical integrations will never reveal the unstable periodic orbit $\gamma_{2}$ which lies inside the invariant torus $\Gamma$.

More generally, since the Langford system is derived by truncating the normal form of a cusp-Hopf singularity, we expect qualitatively similar behavior in systems exhibiting this bifurcation. An interesting topic for future research would be to repeat the numerical analysis performed in the present work for other modifications of the normal form. For example one could perturb the system in such a way that the $z$-axis is no longer invariant. Or, one could modify the system so that the Neimark-Sacker bfurcation is subcritical rather than supercritical.

As remarked already in [20] (and in the introduction of the present work) complex dynamics are often generated by interactions between equilibrium and periodic solutions. The fact that the Langford system is close to a simultaneous cusp-Hopf bifurcation is precisely what provides the multiple equilibrium solutions in close proximity to a limit cycle. This is the basic mechanism organizing many of the interesting dynamical phenomena discussed in the present work.

The normal form unfolding a pitchfork-Hopf bifurcation was also studied by Langford in [84], and it would be a nice project to apply the methods of the present work to this system, or to systems derived from the fold-Hopf bifurcation. Other important normal forms are discussed for example in the woks of [85, 86, 87].

Another interesting topic of future work would be to prove – possibly with computer assistance – as much of Conjecture 5.1 as possible. For example, the following Theorem is found in the Author’s work with Maciej Capinski [88].

Theorem 5.2** (Existence of a $C^{0}$ invariant torus).**

For $\alpha=0.85$ Equation (1) has a $C^{0}$ resonant invariant torus, which is not even globally Lipschitz much less $C^{1}$. More precisely, the torus contains exactly two periodic orbits – one attracting and the other a saddle. The Floquet exponents of the attracting periodic orbit are complex conjugates. The saddle periodic orbit has one stable and one unstable Floquet exponent. The one dimensional unstable manifold of the saddle periodic orbit is completely captured in the basin of attraction of the stable periodic orbit, so that torus is the union of the stable periodic orbit, the saddle periodic orbit, and the unstable manifold of the saddle.

The proof of this theorem is based on the techniques developed in [56, 59, 89] for validating bounds on local manifold parameterizations and computer assisted proofs for heteroclinic connections, the methods developed in [90, 91, 92, 93] for rigorous integration of vector fields and computer assisted proof in Poincaré sections, and the methods of [94, 95] for obtaining validated error bounds on stable/unstable manifolds in Poincaré sections. This one theorem provides a glimpse of what could be accomplished in this and similar systems by constructing computer assisted arguments.

For example, the techniques developed in [91, 92, 96, 97] could be used to prove the existence of the global bifurcations observed above. Using the techniques developed in [98, 99], it should be possible to study in a mathematically rigorous way the global attractor of the Langford system over a large parameter range, and verify and/or clarify many of the claims of Conjecture 5.1. For example parts (1), (2), (3), (4), and (5) appear to us susceptible to this kind of analysis. Combining the techniques of the references just cited with the mathematically rigorous methods for computing stable/unstable invariant manifold atlases developed in [49] could provide means of verifying parts (6) and (7) of the conjecture. We also remark that the recent work of [83] could be applied to give computer assisted proofs for the chaotic attractor after breakdown.

Another project could be to combine the parameterization method for hyperbolic invariant tori developed in [10, 19] with the methods of computer assisted proof developed in [100] to prove the existence of the invariant tori studied in the present work for $\alpha<\alpha_{3}$. That is, to study the tori before the onset of resonance. As mentioned already in [88], the methods of computer assisted proof developed there appear to struggle in the parameter range $\alpha_{1}<\alpha<\alpha_{2}$ due to the apparent lack of uniform contraction rates near the torus. If arguments like the ones suggested in this paragraph could succeed for the Langford system, they then could also be applied to systems coming from other important normal forms.

6 Acknowledgments

The authors would like to thank Jordi-Lluís Figueras, Maciej Capiński, and Vincent Naudot for many helpful suggestions and invaluable insights. We owe also a special thanks to Takahito Mitsui for bringing the paper of Langford [20] to our attention after reading a much rougher early version of this manuscript. The published version of the manuscript benefitted greatly from the suggestions of two anonymous referees. The second author was partially supported by NSF grant DMS-1813501.