Dynamical systems analysis of the Maasch-Saltzman model for glacial cycles

TL;DR

This paper analyzes the internal dynamics of the Maasch-Saltzman model for glacial cycles, revealing how invariant manifolds and bifurcations organize the model's ability to produce glacial oscillations.

Contribution

It provides a dynamical systems analysis of the Maasch-Saltzman model, identifying invariant manifolds and bifurcation structures that explain glacial cycle behavior.

Findings

Long-term dynamics occur on invariant manifolds.

Bogdanov-Takens bifurcations govern the reduced dynamics.

Parameter regions for glacial cycles are organized by bifurcation curves.

Abstract

This article is concerned with the internal dynamics of a conceptual model proposed by Maasch and Saltzman [J. Geophys. Res., 95, D2 (1990) 1955-1963] to explain central features of the glacial cycles observed in the climate record of the Pleistocene Epoch. It is shown that, in most parameter regimes, the long-term system dynamics occur on certain intrinsic two-dimensional invariant manifolds in the three-dimensional state space. These invariant manifolds are slow manifolds when the characteristic time scales for the total global ice mass and the volume of North Atlantic Deep Water are well- separated, and they are center manifolds when the characteristic time scales for the total global ice mass and the volume of North Atlantic Deep Water are comparable. In both cases, the reduced dynamics on these manifolds are governed by Bogdanov-Takens singularities, and the bifurcation curves associated to these singularities organize the parameter regions in which the model exhibits glacial cycles. This work was submitted March 30, 2017.

Equations

$$\begin{array}{rl}\dot{x}& =-x-y,\\ \dot{y}& =ry-pz+s{z}^2-y{z}^2,\\ \dot{z}& =-qx-qz.\end{array}$$

$$\begin{array}{rl}\dot{x}& =-x-y,\\ \dot{y}& =ry-pz-y{z}^2,\\ \dot{z}& =-qx-qz,\end{array}$$

$$\begin{array}{rl}\dot{x}& =-x-y,\\ \dot{y}& =ry-pz-y{z}^2,\\ \epsilon \dot{z}& =-x-z,\end{array}$$

$${\mathcal{M}}_{\epsilon }=\{(x,y,z):z=h_{\epsilon }(x,y)\}.$$

$$\epsilon \frac{d}{dt}{h}_{\epsilon }(x,y)=-x-h_{\epsilon }(x,y),$$

$$h_{\epsilon }(x,y)=h_0(x,y)+\epsilon h_1(x,y)+{\epsilon }^2{h}_2(x,y)+\cdots .$$

$$\begin{array}{rl}{h}_0(x,y)& =-x,\\ h_1(x,y)& =-(x+y),\\ h_2(x,y)& =-(x+y)+(ry+px-x^{2}y),\\ h_3(x,y)& =-(1-2p+4xy)(x+y)+(1-r+x^2)(ry+px-x^{2}y).\end{array}$$

$$\begin{array}{rl}\dot{x}& =-x-y,\\ \dot{y}& =ry-p{h}_{\epsilon }(x,y)-(h_{\epsilon }(x,y){)}^{2}y,\end{array}$$

$$\begin{array}{rl}{x}^{\prime }& =-\epsilon (x+y),\\ y^{\prime }& =\epsilon (ry-pz-y{z}^2),\\ z^{\prime }& =-x-z,\end{array}$$

$$e0=\{r=1-\epsilon p+{\epsilon }^{3}p(p-1)+\mathcal{O}({\epsilon }^4),p>1-\epsilon +{\epsilon }^2-{\epsilon }^3+\mathcal{O}({\epsilon }^4)\}.$$

$$\text{mbox}e1-e2=\{p=1+\epsilon (1-2r)+{\epsilon }^2(1-2r)+{\epsilon }^3(3-8r+4r^2)+\mathcal{O}({\epsilon }^4),r>1-\epsilon +{\epsilon }^2-{\epsilon }^3+\mathcal{O}({\epsilon }^4)\}.$$

$$\overline{x}=\lim \underset{t\to \infty }{\sup }x(t),$$

$$\begin{array}{rl}\dot{x}& =-x-y,\\ \dot{y}& =ry-pz+s{z}^2-y{z}^2,\\ \epsilon \dot{z}& =-x-z,\end{array}$$

$$h_{\epsilon }(x,y)=h_0(x,y)+\epsilon h_1(x,y)+{\epsilon }^2{h}_2(x,y)+{\epsilon }^3{h}_3(x,y)+\mathcal{O}({\epsilon }^4),$$

$$h_2$$

$$h_3$$

$$\begin{array}{rl}\dot{x}& =-x-y,\\ \dot{y}& =ry-p{h}_{\epsilon }(x,y)+s{\left(h_{\epsilon }(x,y)\right)}^2-y{\left(h_{\epsilon }(x,y)\right)}^2.\end{array}$$

$$\begin{array}{rl}\dot{x}& =-x-y,\\ \dot{y}& =ry+px+s{x}^2-y{x}^2.\end{array}$$

$$x_1^{\ast }=\frac{1}{2}[-s+\sqrt{s^2+4(r-p)}],x_2^{\ast }=\frac{1}{2}[-s-\sqrt{s^2+4(r-p)}].$$

$$Q_0=(1,1)\text{ and }{Q}_1=(1+\frac{1}{2}{s}^2,1+\frac{1}{4}{s}^2),$$

$$\begin{array}{rl}\dot{x}& =y,\\ \dot{y}& =(r-p)x+(r-1)y-(s+y)x^2-x^3.\end{array}$$

$$x(t)=\eta u(\tilde{t}),y(t)={\eta }^{2}v(\tilde{t}),\tilde{t}=\eta t,\mu =\frac{r-p}{{\eta }^2},\lambda =\frac{r-1}{{\eta }^2},\delta =\frac{s}{\eta }.$$

$$\begin{array}{rl}\dot{u}& =v,\\ \dot{v}& =\mu u-\delta u^2-u^3+\eta (\lambda -u^2)v,\end{array}$$

$$\begin{array}{rl}\dot{u}& =v,\\ \dot{v}& =u-\delta u^2-u^3,\end{array}$$

$$u_1^{\ast }(\delta )=\frac{1}{2}\left(-\delta +\sqrt{{\delta }^2+4}\right),u_2^{\ast }(\delta )=\frac{1}{2}\left(-\delta -\sqrt{{\delta }^2+4}\right).$$

$$\begin{array}{rl}{\Gamma }_0^{\pm }:\left(u_0^{\pm }(t),v_0^{\pm }(t)\right)& =\left(\frac{\pm 3\alpha }{\cosh (t)\pm \alpha \delta },\frac{\mp 3\alpha \sinh (t)}{{\left(\cosh (t)\pm \alpha \delta \right)}^2}\right),\end{array}$$

$$\begin{array}{rl}D(\lambda )& =\eta {\int }_{-\infty }^{\infty }\left(\lambda -(u_0^{\pm }(t){)}^2\right)(v_0^{\pm }(t){)}^{2}dt+\mathcal{O}({\eta }^2)=\eta \left(\lambda I_0^{\pm }(\delta )-I_2^{\pm }(\delta )\right)+\mathcal{O}({\eta }^2),\end{array}$$

$$\begin{array}{rl}{I}_0^{\pm }(\delta )& ={\int }_{-\infty }^{\infty }(v_0^{\pm }(t){)}^{2}dt=\frac{2}{3}\left(3{\varphi }^2-1\right)+4\psi {\varphi }^2{\tan }^{-1}\left(\psi \mp \varphi \right),\\ I_2^{\pm }(\delta )& ={\int }_{-\infty }^{\infty }(u_0^{\pm }(t)v_0^{\pm }(t){)}^{2}dt=\frac{16}{15}+\frac{23}{3}{\psi }^2+7{\psi }^4+2\psi {\varphi }^2\left(3{\varphi }^2+4{\psi }^2\right){\tan }^{-1}\left(\psi \mp \varphi \right),\end{array}$$

$$\lambda ={\lambda }_0^{\pm }(\delta )=\frac{I_2^{\pm }(\delta )}{I_0^{\pm }(\delta )}+\mathcal{O}(\eta ),$$

$$\lambda ={\lambda }_0^{\cup }(\delta )=\frac{I_2^{+}+I_2^{-}}{I_0^{+}+I_0^{-}}+\mathcal{O}(\eta ).$$

$$(\lambda -\mu )(r-1)=\lambda (p-1).$$

Full text

1 Introduction

The dynamics of glacial cycles during the Pleistocene Epoch—the period from approximately 2.6 million years before present (2.6 Myr BP) until approximately 11.7 thousand years before present (11.7 Kyr BP)—are of great current interest in the geosciences community, see [42], [10, §11] and [33, §12.3]. The geological record shows cycles of advancing and retreating continental glaciers, mostly at high latitudes and high altitudes, and especially in the Northern Hemisphere. The typical temperature pattern inferred from proxy data resembles that of a sawtooth wave, where a slow glaciation is followed by a rapid deglaciation. In the early Pleistocene (until approximately 1.2 Myr BP), the period of a glacial cycle averaged 40 Kyr; after a transition period of approximately 400 Kyr, the glacial cycles had a noticeably greater amplitude and their period averaged 100 Kyr. Although the periods appear to correlate to the cycles of the orbital forcing (Milankovitch theory [35]), the evidence is subject to debate [26, § 11.8], and there is currently no widely-accepted explanation for the mid-Pleistocene transition, when the period of the cycles changed from 40 Kyr to 100 Kyr. Several models have been proposed to explain the various observations; see, for example, [1, 2, 7, 19, 24, 22, 23, 32, 36, 37, 38, 39, 40, 41, 43, 44, 45, 47]. We refer the reader to [8] for an overview of these various modeling efforts and to [42] for a general introduction to paleoclimate modeling. The present investigation focuses on the internal dynamics of the conceptual model developed by Maasch and Saltzman [32, 43].

1.1 The Maasch–Saltzman Model

The Maasch and Saltzman (MS) model is based on physical arguments and emphasizes the role of atmospheric $\mathrm{CO}_{2}$ in the development and evolution of the glacial cycles. In nondimensional form, it consists of the following three ordinary differential equations:

[TABLE]

The state variables $x$, $y$, and $z$ represent the anomalies (deviations from long-term averages) of the total global ice mass, the atmospheric $\mathrm{CO}_{2}$ concentration, and the volume of the North Atlantic Deep Water (NADW), respectively. The latter is a measure of the strength of the North Atlantic overturning circulation and thus of the strength of the oceanic $\mathrm{CO}_{2}$ pump. The parameters $p$, $q$, $r$, and $s$ are combinations of various physical parameters. They are all positive, with $q>1$ for physical reasons. The derivation of the model is given in [43, §2].

In [32], Maasch and Saltzman showed computationally that the model (1.1) exhibits oscillatory behavior with dominating periods of 40 Kyr in response to insolation forcing with such periods, and limit cycles with 100 Kyr periods if $p=1,q=1.2,r=0.8,$ and $s=0.8$ in the absence of forcing. They also showed that a transition from a 40 Kyr cycle to a 100 Kyr cycle can be achieved by slowly varying the parameters $p$ and $r$ across a certain threshold.

Figure 1.1 shows a representative 100 Kyr limit cycle. Each cycle is clearly asymmetric: a rapid deglaciation is followed by a slow glaciation. This asymmetry arises in (1.1) for $s>0$. Also, the three variables are properly correlated: as the concentration of the atmospheric $\mathrm{CO}_{2}$ (a greenhouse gas) increases, the climate gets warmer, and the total ice mass decreases (deglaciation); as the volume of NADW increases, the strength of the North Atlantic overturning circulation increases, more atmospheric $\mathrm{CO}_{2}$ is absorbed by the ocean and, consequently, the atmospheric $\mathrm{CO}_{2}$ concentration decreases.

In this article, we present a dynamical systems analysis of the internal dynamics of the Maasch–Saltzman (MS) model (1.1). We identify the Bogdanov-Takens (BT) points [3, 4, 20, 30, 46] that act as organizing centers in the parameter space for all of the equilibria, limit cycles, homoclinic orbits, and their bifurcations. In addition to being of intrinsic interest, our analysis of the internal dynamics of (1.1) will be instrumental for investigations of the effects of time-dependent forcing, especially of orbital (Milankovitch) forcing, and of the effects of slowly-varying parameters $p$ and $r$.

Two observations are useful for the analysis of the MS model. First, in the special case $s=0$, the system (1.1) reduces to

[TABLE]

which possesses a reflection symmetry; if $(x,y,z)$ is a solution then so is $(-x,-y,-z)$. Hence, it will be useful to study the symmetric MS model (1.2) and to use the results to understand how the dynamics change for $s>0$, as the physical symmetry is broken.

Second, the parameter $q$ is essentially the ratio of the characteristic time scales for the total global ice mass ($x$) and the volume of NADW ($z$). It turns out to be useful to consider first the asymptotic case where $q$ is finite but large, denoted by $q\gg 1$. In this case, the MS model (either in its original form or in its symmetric form) becomes a slow–fast system, where the fast variable $z$ is slaved to the slow variables $x$ and $y$. The insights gained from this analysis will then be a useful guide for understanding the dynamics for all finite $q>1$.

1.2 Summary of the Results

The first results are for the slow–fast regime of the symmetric MS model (1.2). In this regime, the system is (2+1)-dimensional, with two slow variables and one fast variable. We refer to this as the slow–fast symmetric model. We show that there is a family of two-dimensional slow invariant manifolds to which all solutions quickly relax, and we study the dynamics on the slow manifolds. The central feature is a $\mathbb{Z}_{2}$-symmetric BT bifurcation point, from which all bifurcation curves emanate. The curves of Hopf bifurcations, homoclinic bifurcations, and saddle-node bifurcations of limit cycles determine the regions in parameter space where the stable limit cycles exist. In addition, since all solutions relax quickly to the slow manifolds, one can determine the basins of attraction of the various limit cycles. These first results build naturally on the recent analysis of the symmetric MS model (1.2) in the limit $q=\infty$ [15].

The second results concern the effects of asymmetry ($s>0$). In the regime of finite but large values of $q$ ($q\gg 1$), the system (1.1) with $s>0$ is also a slow–fast system. We refer to it as the slow–fast asymmetric system. There is again a family of exponentially attracting two-dimensional, invariant slow manifolds, but the symmetry-breaking makes the dynamics on the slow manifolds more complex. With $s>0$, the limit cycles observed in (1.1) are asymmetric, exhibiting a relatively rapid deglaciation and a relatively slow glaciation, as shown in Figure 1.1.

With these results in hand, we are then in a position to analyze and visualize the dynamics of the full, asymmetric ($s>0$) MS model (1.1) for all $q>1$. We show that for all $q>1$ the system possesses a family of two-dimensional center manifolds toward which solutions relax. Moreover, the solutions of the full system may be accurately approximated by those of the reduced systems on the center manifolds for all $q>1$, and the manifold is at least $C^{1}$-smooth for all $q$ greater than a critical value $q_{c}(p,r,s)$. On the center manifolds, the system has a pair of BT singularities, and the bifurcation curves emanating from them organize the system dynamics, including the boundaries of the regions where stable limit cycles exist.

1.3 Outline of the Article

The article is organized as follows. In Section 2, we present the analysis of the slow–fast symmetric system. In Section 3, we present the analysis of the slow–fast asymmetric system. Section 4 presents the analysis of the full three-dimensional MS model. We conclude with a discussion in Section 5. Appendix A provides details of the unfolding and Melnikov analysis for the persistence of homoclinic orbits in the slow–fast asymmetric system. Appendix B presents essential information about the center manifolds in (1.1).

2 Slow–Fast Dynamics of the Symmetric Model

In this section, we analyze the symmetric MS model (1.2) for large $q$ ($q\gg 1$). The system is readily formulated as a slow–fast system,

[TABLE]

where $\varepsilon=1/q$ is the small parameter, which measures the separation of time scales. We refer to (2.1) as the slow–fast symmetric model. Here, $x$ and $y$ are slow, and $z$ is fast. We assume that $p$ and $r$ are $\mathcal{O}(1)$ with respect to $\varepsilon$.

2.1 Slow Manifolds

With $\varepsilon=0$, the system (2.1) has a critical manifold $\mathcal{M}_{0}=\{(x,y,z):z=-x\}$, which is invariant under the flow. Since $(\partial/\partial z)(-x-z)=-1$ for all $(x,y)$, $\mathcal{M}_{0}$ is normally attracting. By Fenichel’s Geometric Singular Perturbation Theory [18, 25, 27], normal hyperbolicity of $\mathcal{M}_{0}$ implies that, for any sufficiently small and positive $\varepsilon$, there exists a family of persistent normally attracting invariant slow manifolds,

[TABLE]

The functions $h_{\varepsilon}$ are $C^{k}$ for any $k>0$. They are solutions of the invariance equation

[TABLE]

which satisfy $\displaystyle\lim_{\varepsilon\to 0}h_{\varepsilon}(x,y)=-x$. The functions $h_{\varepsilon}$ are identical to all orders in powers of $\varepsilon$ and differ only by terms of $\mathcal{O}(e^{-c/\varepsilon})$ as $\varepsilon\to 0$, for some $c>0$. The expansion may be represented by

[TABLE]

The terms in this expansion are found by substituting (2.4) into (2.3) and equating coefficients of like powers of $\varepsilon$. All terms in $h_{\varepsilon}$ are of odd degree in the variables, due to the $\mathbb{Z}_{2}$ symmetry of (2.1). The first few terms are

[TABLE]

On the slow manifolds $\mathcal{M}_{\varepsilon}$, the system (2.1) reduces to the planar system,

[TABLE]

where $h_{\varepsilon}$ is given by (2.4) and (2.5). The system (2.6) is the object of study in this section.

For completeness, we observe that the fast dynamics, along which solutions relax to $\mathcal{M}_{\varepsilon}$, may be analyzed by introducing the fast time $\tau=t/\varepsilon$ and rewriting (2.1) as a fast system

[TABLE]

where the prime denotes differentiation with respect to $\tau$. Systems (2.1) and (2.7) are equivalent for all $\varepsilon\not=0$. In the limit as $\varepsilon\to 0$, the fast system reduces to a single equation for $z$, with $x$ and $y$ constant in time. Let $(x_{0},y_{0},z_{0})$ denote an arbitrary initial condition. In the fast time, $\tau$, the solution with initial condition $(x_{0},y_{0},z_{0})$ relaxes to the point $(x_{0},y_{0},-x_{0})\in\mathcal{M}_{0}$. Then, for $\varepsilon$ small and positive, the $z$ component again relaxes quickly on the $\tau$ time scale, now to $\mathcal{M}_{\varepsilon}$, and the $x$ and $y$ components will only change slowly guided by the dynamics on $\mathcal{M}_{\varepsilon}$, see Section 2.3.

2.2 Global Bifurcations

The trivial state, $P_{0}=(0,0)$, is an equilibrium of the system (2.6) on $\mathcal{M}_{\varepsilon}$ for all values of the parameters $(p,r)$. In the regime $r>p$, there are also equilibria at $P_{1}=(x_{1}^{*},-x_{1}^{*})$ and $P_{2}=(x_{2}^{*},-x_{2}^{*})$ on $\mathcal{M}_{\varepsilon}$, where $x_{1}^{*}=\sqrt{r-p}$ and $x_{2}^{*}=-\sqrt{r-p}$. The equilibria $P_{1}$ and $P_{2}$ emerge in a symmetric pitchfork bifurcation from $P_{0}$ along the diagonal $r=p$.

At the point $(p,r)=(\frac{1}{1+\varepsilon},\frac{1}{1+\varepsilon})$ in the parameter space, the equilibrium $P_{0}$ of (2.6) undergoes a $\mathbb{Z}_{2}$-symmetric BT bifurcation, since the Jacobian has a zero eigenvalue of geometric multiplicity two. This point is referred to as an organizing center [21], and we denote it by $Q$. All bifurcation curves emanate from this point, including a pair of Hopf bifurcation curves, a homoclinic bifurcation curve, and a curve of saddle-node bifurcations of limit cycles; see [30, §8.4] and [31]. The organizing center and these bifurcation curves are shown in Figure 2.1. The representations of the global bifurcation curves emanating from the organizing centers, shown here, as well as those shown throughout this article, are obtained from numerical continuation using the software package AUTO [11, 12, 13].

The equilibrium $P_{0}$ undergoes a supercritical Hopf bifurcation along the (red) curve, which emanates to the right from the organizing center. The curve is given by

[TABLE]

This curve is quadratic in $p$ and concave up, and the tangent line has slope $-\varepsilon$ at the organizing center to leading order. Similarly, $P_{1}$ and $P_{2}$ undergo subcritical Hopf bifurcations along the (red) curve emanating to the left from $Q$. The curve is given by

[TABLE]

This curve is quadratic in $r$ and concave up, and the tangent line has slope $-1/(2\varepsilon)$ at the organizing center to leading order.

The homoclinic bifurcation curve is shown in blue. Its existence is established using standard BT theory, which entails an unfolding procedure followed by a Melnikov analysis. This analysis shows that, to leading order, at the organizing center, the tangent line to the homoclinic bifurcation curve is $r-(1-\varepsilon)=-4(p-(1-\varepsilon))$, see also Section 3.3.

The curve of saddle-node bifurcations of limit cycles is shown in black. The existence of this second global bifurcation curve is a consequence of the $\mathbb{Z}_{2}$ symmetry and also follows from BT theory [30, 31]. The analysis also shows that, to leading order, this curve is tangent to the line $r-(1-\varepsilon)\approx-3.03(p-(1-\varepsilon))$ at the organizing center.

Remark 2.1**.**

In the limit $\varepsilon\to 0$ ($q\to\infty$), the slow manifolds $\mathcal{M}_{\varepsilon}$ approach the critical manifold $\mathcal{M}_{0}$, and (2.6) reduces to the two-dimensional symmetric model studied in [15, § 4].

Figure 2.2 shows the dynamics of (2.1) in terms of the function $(p,r)\mapsto\overline{x}(p,r)$ for $(p,r)\in(0,1.5)\times(0,2)$ as a color map. Here, $\overline{x}$ is the function defined by

[TABLE]

and the initial conditions were chosen randomly for each $(p,r)$. Below the diagonal and e0, the color map is light green since all solutions approach $P_{0}$. As one crosses e0 with increasing $r$, the color changes from light green to orange and then to pink, indicating the presence of periodic orbits which are created in the supercritical Hopf bifurcation. Between the Hopf bifurcation curve e1-e2 and the homoclinic bifurcation curve, there is a similar shift to pink as $r$ increases. One also sees some green patches, indicating that some of the randomly chosen initial conditions lie in the basin of attraction of the stable equilibrium $P_{2}$. Next, in the region between the homoclinic bifurcation curve and the curve of saddle-node bifurcations of limit cycles, the color map has largely the same color, because solutions with initial conditions that lie inside the large unstable limit cycle approach one of the stable equilibria (green or orange), and those with initial conditions outside the unstable limit cycle approach the large stable limit cycle (pink). Finally, below the curve of saddle-node bifurcations of limit cycles, the color map consists entirely of green and orange, indicating that all of the solutions are attracted either to $P_{1}$ or to $P_{2}$, as expected since there are no stable limit cycles.

2.3 Basins of Attraction

An important feature of slow–fast systems like (2.1) is that Geometric Singular Perturbation Theory provides insight into the dynamics not only of solutions on the slow manifold $\mathcal{M}_{\varepsilon}$, but also of solutions in a neighborhood of $\mathcal{M}_{\varepsilon}$ and thus provides a way to explore the basins of attraction of different solutions—such as equilibria and limit cycles—on $\mathcal{M}_{\varepsilon}$.

In a neighborhood of $\mathcal{M}_{\varepsilon}$, any solution $X(t)=(x(t),y(t),z(t))$ of (2.1) decomposes into fast and slow components. The fast component is directed along the fast fiber $\mathcal{F}_{\varepsilon}(b(t))$, which stands above the base point $b(t)$ of the fiber on $\mathcal{M}_{\varepsilon}$. It decays exponentially fast toward the manifold $\mathcal{M}_{\varepsilon}$, $\|X(t)-b(t)\|\leq e^{-ct/\varepsilon}\|X(0)-b(0)\|$ for some positive constant $c$, which is of $\mathcal{O}(1)$ as $\varepsilon\to 0$. The family of fast fibers is invariant in the sense that $\mathcal{F}_{\varepsilon}(b(t))=\phi_{t}\left(\mathcal{F}_{\varepsilon}(b(0))\right)$ for all $t$, where $\phi_{t}$ is the time-$t$ flow map of (2.1). The slow component of $X(t)$ lies in the tangent plane to the slow manifold $\mathcal{M}_{\varepsilon}$ at the base point $b(t)$.

We now demonstrate —in two cases— how one can use this decomposition to investigate the basins of attraction of solutions on $\mathcal{M}_{\varepsilon}$. The first case is shown in Figure 2.3. Here, $(p,r)$ is in the region between e0 and the diagonal where $P_{0}$ is an unstable spiral and there is a unique stable limit cycle $\gamma_{\varepsilon}$ of (2.1) around $P_{0}$.

Let $\mathcal{B}$ denote the basin of attraction of $\gamma_{\varepsilon}$. The set $\mathcal{B}|_{\mathcal{M}_{\varepsilon}}$ of all initial conditions on $\mathcal{M}_{\varepsilon}$ that lie in $\mathcal{B}$ is completely determined by the analysis of the slow flow (2.6) on $\mathcal{M}_{\varepsilon}$. For any initial condition $X(0)=(x(0),y(0),z(0))$ that lies near $M_{\varepsilon}$ but not on it, there is a unique fast stable fiber that contains $X(0)$. Let $b(0)$ be the base point of this fiber on $\mathcal{M}_{\varepsilon}$, and denote the fiber by $\mathcal{F}_{\varepsilon}(b(0))$. Since any solution $X(t)$ near $\mathcal{M}_{\varepsilon}$ decomposes into a fast component and a slow component, we know immediately that $X(0)$ lies in $\mathcal{B}$ whenever $b(0)$ lies in $\mathcal{B}|_{\mathcal{M}_{\varepsilon}}$. Moreover, since the fast stable fibers completely foliate the neighborhood of $\mathcal{M}_{\varepsilon}$ and each initial condition $X(0)$ near $\mathcal{M}_{\varepsilon}$ lies on a unique fiber, one may apply the above analysis to each initial condition $X(0)$ near $\mathcal{M}_{\varepsilon}$. Therefore, the basin of attraction $\mathcal{B}$ of $\gamma_{\varepsilon}$ contains all solutions near $\mathcal{M}_{\varepsilon}$, except those that start exactly on the stable fiber $\mathcal{F}_{\varepsilon}(P_{0})$ with base point $P_{0}$. ($P_{0}$ is not in $\mathcal{B}|_{\mathcal{M}_{\varepsilon}}$.)

The second case is shown in Figure 2.4. Here, $(p,r)$ is in the region between the Hopf bifurcation curve e1-e2 and the homoclinic bifurcation curve; $P_{0}$ is a saddle, $P_{1}$ and $P_{2}$ are stable spirals, each surrounded by an unstable limit cycle. The limit cycle around $P_{1}$ is labeled $\gamma_{1,\varepsilon}$; the limit cycle $\gamma_{2,\varepsilon}$ around $P_{2}$ is not shown. Also, there is a large stable limit cycle, labeled $\gamma_{3,\varepsilon}$, which surrounds $P_{0},P_{1},$ and $P_{2}$, and the unstable limit cycles.

There are three primary basins of attraction, one for each of the stable equilibria $P_{1}$ and $P_{2}$, and one for the large stable limit cycle $\gamma_{3,\varepsilon}$, denoted by $\mathcal{B}_{1}$, $\mathcal{B}_{2}$, and $\mathcal{B}_{3}$, respectively. The sets of initial conditions in these primary basins are determined as follows.

First, we observe that all of the initial conditions on the slow manifold $\mathcal{M}_{\varepsilon}$ that lie inside the unstable limit cycle $\gamma_{1,\varepsilon}$ are in $\mathcal{B}_{1}$. Denote this set of initial conditions by $\mathcal{B}_{1}|_{\mathcal{M}_{\varepsilon}}$. Then, any initial condition $X(0)$ not on $\mathcal{M}_{\varepsilon}$ is in $\mathcal{B}_{1}$ if it lies on the fast stable fiber $\mathcal{F}_{\varepsilon}(b(0))$ of an initial condition $b(0)$ in $\mathcal{B}_{1}|_{\mathcal{M}_{\varepsilon}}$. This completely determines the basin $\mathcal{B}_{1}$, see Figure 2.4.

Similarly, the basin of attraction $\mathcal{B}_{2}$ of the stable equilibrium $P_{2}$ consists of (i) the set of all initial conditions on $\mathcal{M}_{\varepsilon}$ that lie inside the unstable limit cycle $\gamma_{2,\varepsilon}$, a set which we denote by $\mathcal{B}_{2}|_{\mathcal{M}_{\varepsilon}}$, and (ii) the set of all initial conditions $X(0)$ (not on $\mathcal{M}_{\varepsilon}$) that lie on fast stable fibers $\mathcal{F}_{\varepsilon}(b(0))$ with base points $b(0)$ inside $\mathcal{B}_{2}|_{\mathcal{M}_{\varepsilon}}$.

Finally, the basin of attraction $\mathcal{B}_{3}$ of the large stable limit cycle $\gamma_{3,\varepsilon}$ consists of (i) the set of all initial conditions on $\mathcal{M}_{\varepsilon}$ that are exterior to $\gamma_{1,\varepsilon}$ and $\gamma_{2,\varepsilon}$ and that do not lie on the stable manifold of $P_{0}$, a set which we denote by $\mathcal{B}_{3}|_{\mathcal{M}_{\varepsilon}}$, and (ii) the set of all initial conditions $X(0)$ not on $\mathcal{M}_{\varepsilon}$ which lie on fast stable fibers whose base points are in $\mathcal{B}_{3}|_{\mathcal{M}_{\varepsilon}}$.

The second case is more complicated than the first due to the presence of the unstable saddle at $P_{0}$. However, the two cases are representative, and the methods are similar for finding the basin of attraction of any other limit cycle on $\mathcal{M}_{\varepsilon}$.

3 The Slow–Fast Asymmetric System

In this section, we study the effect of asymmetry ($s>0$) in (1.1) in the limit of large $q$,

[TABLE]

where $\varepsilon=1/q$, as before. We refer to (3.1) as the slow–fast asymmetric model. Just as was the case for the symmetric model (2.1), this slow–fast system has a family of slow invariant manifolds, on which the long–term system dynamics occur and on which the limit cycles lie. By studying the flow on these slow manifolds, we establish the existence of two nondegenerate BT points. For each $s>0$, the two nondegenerate BT points are organizing centers in the $(p,r)$ plane. The geometry of these organizing centers and bifurcation curves (and hence also the geometry of the regions where the limit cycles lie) may be understood as a symmetry-breaking of the lone $\mathbb{Z}_{2}$-symmetric BT point studied in Section 2.2.

3.1 Slow Manifolds

With $\varepsilon=0$, the slow–fast asymmetric model has the same normally hyperbolic critical manifold $\mathcal{M}_{0}=\{(x,y,z):z=-x\}$ as the slow–fast symmetric model (2.1). This critical manifold persists for all $0<\varepsilon\ll 1$. In particular, there is a family of invariant slow manifolds $\mathcal{M}_{\varepsilon}$, which are given to all orders by the graph of a function $h_{\varepsilon}(x,y)$,

[TABLE]

and this function depends on $s$. The first two coefficients, $h_{0}$ and $h_{1}$, in the expansion of $h_{\varepsilon}$ are the same as in (2.5). Then, one finds

[TABLE]

On the slow manifolds $\mathcal{M}_{\varepsilon}$, the system (3.1) reduces to the planar system

[TABLE]

To simplify the presentation, we focus on the dynamics of (3.3) on the critical manifold. That is, we represent $h_{\varepsilon}$ by the leading order term $h_{0}=-x$, which gives

[TABLE]

The dynamics of (3.3) for $0<\varepsilon\ll 1$ are regular perturbations of those of (3.4).

3.2 Organizing Centers

The origin $P_{0}=(0,0)$ is an equilibrium of (3.4) for all $p$, $r$, and $s$. If $r>p-\tfrac{1}{4}s^{2}$, there are two additional equilibria, namely $P_{1}=(x_{1}^{*},-x_{1}^{*})$ and $P_{2}=(x_{2}^{*},-x_{2}^{*})$, where

[TABLE]

Figure 3.1 shows the results of a linear stability analysis of the equilibria of (3.4).

We refer to the line $r=p-\tfrac{1}{4}s^{2}$ as the shifted diagonal (marked d2 and sd in Figure 3.1). Note that $x_{2}^{*}<x_{1}^{*}<0$ if $p-\tfrac{1}{4}s^{2}<r<p$, and $x_{2}^{*}<0<x_{1}^{*}$ if $r>p$.

At the points

[TABLE]

the system (3.4) has a zero eigenvalue of geometric multiplicity two (at $P_{0}$ and $P_{1}$, respectively). Thus, $Q_{0}$ and $Q_{1}$ are organizing centers.

Figure 3.2 shows the organizing centers and the branches of global bifurcations of (3.4). There are two branches of homoclinic bifurcations emanating from $Q_{0}$ (solid blue), one to the left and one to the right, and a single branch emanating from $Q_{1}$ to the right (dashed blue). The companion Figure 3.3 shows the types of homoclinic orbits along these branches. In Figure 3.3, the equilibria $P_{0},P_{1},$ and $P_{2}$ are marked by black, red, and green dots, respectively. These branches, and the homoclinic orbits along them, are as follows:

• •

The branch emanating from $Q_{0}$ to the left consists of right-homoclinic orbits to $P_{0}$, which enclose $P_{1}$. A representative orbit is shown in Figure 3.3(a).

• •

The branch emanating from $Q_{0}$ to the right consists of three segments:

• –

A segment from $Q_{0}$ to the shifted diagonal, most clearly seen in Figure 3.2(b); this branch consists of right-homoclinic orbits to $P_{1}$, which enclose $P_{0}$. A representative orbit is shown in Figure 3.3(b).

• –

A segment from the shifted diagonal to the diagonal, most clearly seen in Figure 3.2(a) and (c); this branch consists of large-amplitude double-loop homoclinic orbits to $P_{1}$, which enclose $P_{0}$ and $P_{2}$. A representative orbit is shown in Figure 3.3(c).

• –

A segment beyond the diagonal, most clearly seen in Figure 3.2(a); this branch consists of large-amplitude double-loop homoclinic orbits to $P_{0}$, which enclose $P_{1}$ and $P_{2}$. A representative orbit is shown in Figure 3.3(d).

• •

The branch of homoclinic bifurcations emanating from $Q_{1}$ (dashed blue) is made up of two segments:

• –

A segment from $Q_{1}$ to the diagonal, most clearly seen in Figure 3.2(c); this branch consists of left-homoclinics to $P_{1}$, which enclose $P_{2}$. A representative orbit is shown in Figure 3.3(e).

• –

A segment beyond the diagonal, most clearly seen in Figure 3.2(a); this branch consists of left-homoclinics to $P_{0}$, which enclose $P_{2}$. A representative orbit is shown in Figure 3.3(f).

The existence of all six types of homoclinics will be formally proven in Section 3.3. Figure 3.2 also shows two branches of saddle-node bifurcations of limit cycles (black curves). We note that these curves are no longer attached to the organizing centers due to the asymmetry, i.e., since the BT points are non-degenerate. The local and global bifurcation curves serve as boundaries between the regions of different dynamical behaviour.

Figure 3.4 shows the sequence of phase portraits that are encountered for a fixed value of $p>1+\tfrac{1}{2}s^{2}$ as $(p,r)$ moves through the different regions for increasing values of $r$. There is one phase portrait for each of the six small black diamonds in Figure 3.2.

In Figure 3.4(a), the parameters are chosen in the region bounded by e0 and the shifted diagonal, where $P_{0}$ is an unstable spiral surrounded by a stable limit cycle (black curve). In frame (b), $P_{1}$ (red dot) and $P_{2}$ (green dot) lie inside the large stable limit cycle, for $(p,r)$ values in the region between the shifted diagonal and e2. In frame (c), we see the unstable limit cycle (red) surrounding $P_{2}$ that lies in the region bounded by e2 and the dashed curve of homoclinics. Next, frame (d) shows a representative phase portrait obtained after the homoclinic bifurcation curve (dashed blue curve in Figure 3.2) is crossed. One sees that all solutions (not on the stable manifold of $P_{0}$) are forward asymptotic to $P_{2}$ (green dot) or to the large stable limit cycle (black curve). In frame (e), there are large-amplitude unstable (red) and stable (black) limit cycles, in the narrow region between the homoclinic bifurcation curve (blue) and the curve of saddle-node bifuractions of limit cycles (black). These disappear as one crosses the curve of saddle-node bifurcations of limit cycles (upper black curve in Figure 3.2). In frame (f), only the three equilibria remain.

Remark 3.1**.**

The bifurcation structure of (3.4) collapses to that of the symmetric case (2.6) as $s\to 0$ (see Figure 3.5). More specifically, the shifted diagonal line of saddle-node bifurcations of $P_{1}$ and $P_{2}$ collapses onto the diagonal line of transcritical bifurcations of $P_{0}$ and $P_{1}$, thus creating the diagonal line of pitchfork bifurcations of system (2.6). The organizing centers $Q_{0}$ and $Q_{1}$ merge and become the $\mathbb{Z}_{2}$-symmetric BT point at $Q=(1,1)$ as $s\to 0$. Concomitantly, the curves, e1 and e2, of Hopf bifurcations merge to the single curve labelled e1-e2 in Figure 2.2. Additionally, the curves of homoclinic bifurcations merge into a single curve (which is tangent to the line $r-1=-4(p-1)$ at $Q$) as $s\to 0$. Similarly, the curves of saddle-node bifurcations of limit cycles collapse to a single curve (which is tangent to the line $r-1\approx-3.03(p-1)$ at $Q$) as $s\to 0$.

3.3 Bogdanov-Takens Unfolding Analysis

In this section, we present the unfolding analysis of the non-degenerate BT points $Q_{0}$ and $Q_{1}$ in system (3.4). Specifically, we study the homoclinic orbits of an appropriate Hamiltonian system and use Melnikov theory [20, 34] to determine the parameter sets for which these homoclinics persist under small perturbations. In this manner, we formally prove the persistence of the six types of homoclinic orbits, which lie along the blue branches that emanate from $Q_{0}$ and $Q_{1}$ in Figure 3.3. Figure 3.6 summarizes the results of this subsection.

3.3.1 Rescaling and Partition of the Parameter Plane

First, we make the change of variables $(x,-(x+y))\mapsto(x,y)$, so that (3.4) becomes

[TABLE]

The equilibria are $P_{0}=(0,0)$, $P_{1}=(x_{1}^{*},0)$, and $P_{2}=(x_{2}^{*},0)$, where $x_{1}^{*}$ and $x_{2}^{*}$ are again given by (3.5). We rescale the dependent and independent variables and the parameters by

[TABLE]

With these rescaled variables and parameters, system (3.7) is equivalent with

[TABLE]

where the overdot denotes $\frac{d}{d{\tilde{t}}}$, and we drop the tildes.

There are four distinct regions in the $(p,r)$ plane, depending on $\mu$ and $\lambda$ (Figure 3.6):

1. I.

The set $\{\mu>0,\,\lambda>0\}$ corresponds to the region of the $(p,r)$ plane above the main diagonal and above the line $r=1$ (red shaded region). 2. II.

The set $\{\mu>0,\,\,\lambda<0\}$ is the region above the main diagonal and below the line $r=1$ (blue shaded region). 3. III.

The set $\{-\tfrac{1}{4}\delta^{2}<\mu<0,\,\,\lambda>0\}$ is the region enclosed by the two diagonals, and lies above the line $r=1$ (purple shaded region). 4. IV.

The set $\{-\tfrac{1}{4}\delta^{2}<\mu<0,\,\,\lambda<0\}$ corresponds to the region enclosed by the two diagonals and lies below the line $r=1$ (green shaded region).

The organizing centers are $Q_{0}=(0,0)$ and $Q_{1}=(-\tfrac{1}{4}\delta^{2},\tfrac{1}{4}\delta^{2})$ in the $(\mu,\lambda)$ plane, corresponding to $P_{0}=(0,0)$ and $P_{1}=(-\tfrac{1}{2}\delta,0)$, respectively.

3.3.2 Persistence of the Homoclinics to $P_{0}$

To study the homoclinics to $P_{0}$, we consider the region of $\mu>0$. Without loss of generality, we set $\mu=1$ in (3.9) and begin with the $\eta=0$ limit of (3.9)

[TABLE]

which is Hamiltonian. The Hamiltonian function is $H(u,v)=\tfrac{1}{2}v^{2}-\tfrac{1}{2}u^{2}+\tfrac{1}{3}\delta u^{3}+\tfrac{1}{4}u^{4}$. System (3.10) has equilibria at $(0,0)$, $(u_{1}^{*},0)$, and $(u_{2}^{*},0)$, where

[TABLE]

The saddle $P_{0}$ is connected to itself by a pair of asymmetric homoclinic orbits, $\Gamma^{-}_{0}$ and $\Gamma^{+}_{0}$, which surround $(u_{2}^{*},0)$ and $(u_{1}^{*},0)$, respectively. They are given explicitly by

[TABLE]

where $\alpha=\sqrt{\frac{2}{9+2\delta^{2}}}$. The phase portrait of (3.10) is shown in Figure 3.7, and we note that the solutions $\Gamma_{0}^{\pm}$ limit exactly, as $\delta\to 0$, onto the homoclinic orbits of the unperturbed reduced system on $\mathcal{M}_{0}$ in the slow–fast symmetric system.

Now, for each of the unperturbed homoclinic orbits, $\Gamma^{+}_{0}$ and $\Gamma^{-}_{0}$, we develop the distance function, $D(\lambda)$, as an asymptotic expansion in $\eta$:

[TABLE]

where we define the integrals $I_{0}^{\pm}(\delta)$ and $I_{2}^{\pm}(\delta)$ as

[TABLE]

where $\psi=\tfrac{1}{3}\delta\sqrt{2}$ and $\phi=\tfrac{1}{3}\sqrt{2\delta^{2}+9}$ (and $\phi^{2}-\psi^{2}=1$). We note that $I_{0}^{\pm}(\delta)$ is the area enclosed by $\Gamma_{0}^{\pm}$ in the $(u,v)$ plane, since $\int v^{2}\,dt=\int v\,du$. The bifurcation equation $D(\lambda)=0$ gives the parameter set for which the homoclinics persist. Thus, the homoclinics, $\Gamma_{0}^{\pm}$, to $P_{0}$ persist for

[TABLE]

which is shown in Figure 3.8. This demonstrates the existence of the curves of homoclinic bifurcations to $P_{0}$. More specifically, $\lambda_{0}^{+}$ corresponds to the solid branch of right-homoclinics in region I that emanates to the left of $Q_{0}$, and these are the right-homoclinic orbits shown in Figure 3.3(a). Also, $\lambda_{0}^{-}$ corresponds to the dashed branch of left-homoclinics in region I that emanates to the right of $Q_{1}$ and lies above the main diagonal, and these are the left-homoclinic orbits shown in Figure 3.3(f).

To establish the existence and persistence of the large-amplitude double-loop homoclinics to $P_{0}$ (see Figure 3.3(d)), we make the following numerical and analytical observations. Numerically, the large-amplitude homoclinics to $P_{0}$, as illustrated in Figure 3.6 for $s=0.8$ for instance, converge in the symmetric (i.e., $s\to 0$) limit to the pair of homoclinics of the slow–fast symmetric model (see Figure 3.9). As such, the appropriate homoclinic orbit along which to measure the distance function is the concatenation, $\Gamma_{0}^{-}\cup\Gamma_{0}^{+}$, of the left and right homoclinics to $P_{0}$. Analytically, it is known from multi-pulse Melnikov theory that the Melnikov functions for double-loop homoclinic orbits are given to leading order by the sum of the individual single-loop Melnikov functions, with the contribution from the passage near the saddle being of higher order, see for example [5]. Thus, the bifurcation equation $D(\lambda)=0$ for the persistence of the large-amplitude homoclinics to $P_{0}$ gives

[TABLE]

The function $\lambda_{0}^{\cup}(\delta)$ is shown in Figure 3.8 (green curve).

The definition of $\mu$ and $\lambda$ (see (3.8)) generates a linear relation between $p$ and $r$

[TABLE]

This formula yields the slope of the tangent lines to the curves of homoclinic bifurcations at the organizing centers. Here, $\mu=1$, and $\lambda$ is positive along both branches of homoclinics. This analysis agrees with the numerical continuation results seen previously.

Remark 3.2**.**

The unfolding and Melnikov analysis for the homoclinic orbits that emanate from $Q_{1}$ is similar. The main difference is that the hyperbolic saddle of the unperturbed Hamiltonian is located at the equilibrium $P_{1}$ instead of $P_{0}$. As such, the first step in the analysis is a translation to place $P_{1}$ at the origin. The details are presented in Appendix A.

Remark 3.3**.**

System (3.7) with $h_{\varepsilon}$ given to leading order by $h_{0}=-x$ is a special case of general four-parameter planar vector fields studied in [9, 28]. These four-parameter vector fields arise as part of the unfolding of vector fields, $\dot{x}=y,\,\,\dot{y}=-x^{3}-x^{2}y,$ which have a $\mathbb{Z}_{2}$-symmetric BT singularity, see [6, 29]. In [9], the four-parameter vector fields are of the form $\dot{x}=y,\,\,\dot{y}=-(x^{3}+ax^{2}+bx+c)+(d-x^{2})y,$ where $a,b,c,$ and $d$ are real numbers, and in [28] the vector field is the same except for one difference, namely the quadratic term is $xy$. We refer to [9] for an extensive analysis, based on elliptic integrals, of the bifurcation curves in this system, including the Hopf bifurcations, homoclinic bifurcations, saddle-node bifurcations of limit cycles, the effects of the $\mathbb{Z}_{2}$ symmetry-breaking, and many other bifurcations. Also, we observe that with higher–order terms in $h_{\varepsilon}$, the system (3.7) involves higher–order polynomials than those studied in [9, 28]. We refer to [14] for the analysis of a related co-dimension three singularity.

4 The Maasch–Saltzman Model

In this section, we bring the results from the previous sections together to study the full asymmetric MS model (1.1) for all $q>1$ and $s>0$,

[TABLE]

The system (4.1) has two organizing centers

[TABLE]

We show that, for values of $(p,r)$ near $Q_{0}$ and all $q>1$, the system (4.1) has a family of two-dimensional center manifolds to $P_{0}$, and there is a critical value $q_{c}(p,r,s)$ such that the manifolds are at least $C^{1}$-smooth for all $q>q_{c}(p,r,s)$. There is a similar result for values of $(p,r)$ near $Q_{1}$, and these center manifolds near $Q_{1}$ coincide with those near $Q_{0}$.

4.1 Equilibria and Bifurcations

The equilibria of (4.1) are the same as those in the previous section. The trivial state $P_{0}=(0,0,0)$ is again a solution for all parameter values. If $s^{2}+4(r-p)>0$, then there are two additional equilibria, $P_{1}=(x_{1}^{*},-x_{1}^{*},-x_{1}^{*})$ and $P_{2}=(x_{2}^{*},-x_{2}^{*},-x_{2}^{*})$, where

[TABLE]

recall (3.5). Also recall that $x_{2}^{*}<x_{1}^{*}<0$ if $p-\tfrac{1}{4}s^{2}<r<p$, and $x_{2}^{*}<0<x_{1}^{*}$ if $r>p$.

Let $P=(x^{*},-x^{*},-x^{*})$ be any of the equilibria, with $x^{*}=0$, $x_{1}^{*}$, or $x_{2}^{*}$. The characteristic equation of the Jacobian at $P$ is $\lambda^{3}+b\lambda^{2}+c\lambda+d=0$, with

[TABLE]

By the Routh–Hurwitz conditions, $P$ is (linearly) stable if $b>0$, $c>0$, $d>0$, and $e=bc-d>0$. We analyze these conditions for fixed $q$ and $s$, considering $b$, $c$, $d$, and $e$ as functions of $p$ and $r$. Figure 4.1 illustrates the results for $q=1.2$ and $s=0.8$ (the values chosen in [32]), using the same color scheme as in Figure 3.2.

Let $\Omega$ be the domain in the positive quadrant of the $(p,r)$ plane where all four inequalities are satisfied. The Routh–Hurwitz conditions cease to be satisfied when at least one of the inequalities becomes an equality. Hence, $\partial\Omega$ consists of segments where at least one of the coefficients $b$, $c$, $d$, or $e$ vanishes.

Consider the points of $\partial\Omega$ where at least one of the coefficients $b$, $c$, $d$, or $e$ vanishes. We claim that at any such point, with possibly finitely many exceptions, $b>0$, $c>0$, and either $d=0$ and $e>0$ or $d>0$ and $e=0$. To prove the claim, we assume that at such a point exactly one of the coefficients $b$, $c$, $d$, or $e$ vanishes, while the remaining three coefficients are all strictly positive. This assumption is true generically—that is, with possibly finitely many exceptions. If $b=0$ or $c=0$, then $e=-d$, so $d$ and $e$ cannot both be positive. Therefore, it must be the case that $b>0$, $c>0$, and either $d=0$ or $e=0$, as claimed.

If $d=0$, the characteristic polynomial reduces to $\lambda(\lambda^{2}+b\lambda+c)=0$, which yields a simple root at the origin and two roots with negative real parts. If $e=0$, it reduces to $(\lambda+b)(\lambda^{2}+c)=0$, which yields a negative real root $-b$ and a conjugate pair of purely imaginary roots $\pm i\sqrt{c}$, indicating that the equilibrium loses stability due to a Hopf bifurcation. The quantity $\sqrt{c}$ is the natural frequency.

4.1.1 Stability of $P_{0}$

If $x^{*}=0$, the coefficients of the characteristic polynomial are $b_{0}=1+q-r$, $c_{0}=q-(1+q)r$, and $d_{0}=q(p-r)$. The condition $d_{0}>0$ is satisfied if $p>r$. The condition $e_{0}=b_{0}c_{0}-d_{0}>0$ is satisfied if

[TABLE]

Since $p>r$ is a necessary condition for stability, the expression in the right member must be positive—that is, either $0<r<\frac{q}{1+q}$ which corresponds to the case where both factors in the right member are positive, or $r>1+q$ which corresponds to the case where both factors in the right member are negative. The latter inequality violates the Routh–Hurwitz condition $b_{0}>0$; hence, stability can occur only if $0<r<\frac{q}{1+q}$.

The zero-level set of $e_{0}$ constitutes a curve, labeled e0 in Figure 4.1, of supercritical Hopf bifurcations of $P_{0}$. The equilibrium $P_{0}$ loses stability, and a stable limit cycle is created as $(p,r)$ crosses e0 going upward. The curve e0 is a parabola in the $(p,r)$ plane, which is open to the right. The lower branch of the parabola is given by

[TABLE]

where $f_{0}(p)$ is monotonically decreasing. (We note that the condition $0<r<\frac{q}{1+q}$ is satisfied if and only if $\frac{q}{1+q}<p<1+q$.) Therefore, the Routh-Hurwitz criteria are satisfied in the region enclosed by the $p$ axis, the diagonal, and the curve e0.

In addition, the Jacobian of (4.1) at $P_{0}$ has a zero eigenvalue of geometric multiplicity two at $Q_{0}$, and a third eigenvalue $\lambda_{3}=-(1+q-\frac{q}{1+q})$. Hence, $Q_{0}$ is an organizing center.

4.1.2 Stability of $P_{1}$

Recall that $P_{1}$ exists if and only if $r\geq p-\tfrac{1}{4}s^{2}$. With $x^{*}=x_{1}^{*}$, we have

[TABLE]

The zero-level set of $d_{1}$ is the diagonal $r=p$. On the diagonal, $P_{0}$ and $P_{1}$ exchange stability in a transcritical bifurcation. The zero-level set, e1, of $e_{1}=b_{1}c_{1}-d_{1}$ emerges from $Q_{0}$ (see Figure 4.1) and corresponds to a curve of supercritical Hopf bifurcations of $P_{1}$. Therefore, the Routh-Hurwitz criteria are satisfied in the region enclosed by the $r$ axis, the diagonal, and the curve e1.

In addition, the Jacobian of (4.1) at $P_{1}$ has a zero eigenvalue of geometric multiplicity two at $Q_{1}$, and a third eigenvalue $\lambda_{3}=-(1+q-\frac{q}{1+q})$. Hence, $Q_{1}$ is an organizing center.

4.1.3 Stability of $P_{2}$

Recall that $P_{2}$, like $P_{1}$, exists if and only if $r\geq p-\tfrac{1}{4}s^{2}$. With $x^{*}=x_{2}^{*}$, we have

[TABLE]

The zero-level set of $d_{2}$ is the shifted diagonal $r=p-\tfrac{1}{4}s^{2}$. On the shifted diagonal, $P_{1}$ and $P_{2}$ are created in a saddle-node bifurcation. The zero-level set, e2, of $e_{2}=b_{2}c_{2}-d_{2}$ is the continuation of e1 (after a gap between $Q_{0}$ and $Q_{1}$), and corresponds to a curve of subcritical Hopf bifurcations of $P_{2}$. At $Q_{1}$, the zero-level set of $e_{2}$ is tangent to the shifted diagonal. Therefore, the Routh-Hurwitz criteria are satisfied in the region enclosed by the axes, the shifted diagonal, and e2.

Along the Hopf bifurcation curve e2, there is a Bautin bifurcation point, marked by the black triangle in Figure 4.1. The Hopf bifurcations are subcritical below the Bautin point and supercritical above it. Also, a (black) branch of saddle-node bifurcations of limit cycles emerges from the Bautin point. For each fixed $s\geq 0$, the Bautin point is close to $Q_{1}$ for $q$ close to one, and then it moves up along the curve e2 as $q$ increases.

Remark 4.1**.**

In the symmetric limit ($s\to 0$), the curves e1 and e2 collapse to the same parabola in the $(p,r)$ plane, given by

[TABLE]

That is, the curves e1 and e2 are different $s$-unfoldings of the parabola (4.8).

4.1.4 Global Bifurcations

The curves of global bifurcations that emanate from the organizing centers $Q_{0}$ and $Q_{1}$ are shown in Figure 4.1. The blue curve is the curve of homoclinic bifurcations of (4.1). The system (4.1) possesses the same six types of single-loop and double-loop (large-amplitude) homoclinic orbits to the saddles $P_{0}$ and $P_{1}$ as the system (3.4) on the slow manifold $\mathcal{M}_{\varepsilon}$. They lie in regions I and III (recall Figure 3.6), where the lower boundaries of these regions are now given by the Hopf bifurcation curves e1 and e0, respectively.

The black curves in Figure 4.1 indicate curves of saddle-node bifurcations of limit cycles. Along these black curves, a pair of limit cycles of (4.1) of opposite stability merge and annihilate each other. The existence of these global bifurcation curves follows directly from the existence of center manifolds (see the next subsection) and the presence of BT points of the reduced equations on the center manifolds.

Remark 4.2**.**

Figure 4.2 shows how the bifurcation structure of (4.1) collapses as $s\to 0$. The shifted diagonal of saddle-node bifurcations of $P_{1}$ and $P_{2}$ converges to the diagonal of transcritical bifurcations, resulting in a diagonal line of pitchfork bifurcations in the symmetric limit. The organizing center $Q_{1}$ collapses onto the organizing center $Q_{0}$, leaving the $\mathbb{Z}_{2}$-symmetric BT point, $Q$, seen in the symmetric MS model (1.2). Consequently, the curves e1 and e2 of Hopf bifurcations of $P_{1}$ and $P_{2}$ merge to the curve (4.8). Similarly, the curves of homoclinic bifurcations merge to a single curve of homoclinics, and the curves of saddle-node bifurcations of limit cycles coalesce and become a single curve.

4.2 Center Manifolds

In this section, we establish the existence of the center manifolds associated to the equilibria $P_{0}$ and $P_{1}$ in the system (4.1) for parameter values near the organizing centers $Q_{0}$ and $Q_{1}$. We follow the general approach described in [20, §3.2]. The derivation of the center manifolds of $P_{0}$ for parameters near $Q_{0}$ is given in detail. Since $P_{1}$ lies on the center manifolds to $P_{0}$, the center manifolds attached to $P_{1}$ are in fact the same as those attached to $P_{0}$ to all orders.

As a first step to establish the center manifolds of $P_{0}$, we transfer $Q_{0}$ to the origin in parameter space by introducing $(\tilde{p},\tilde{r})=(p-\frac{q}{1+q},r-\frac{q}{1+q})$. Thus, (4.1) transforms to

[TABLE]

where

[TABLE]

The matrix $A$ has eigenvalues $\lambda_{1}=0$, $\lambda_{2}=0$, and $\lambda_{3}=-\left(1+q-\frac{q}{1+q}\right)$. Next, we change coordinates to reduce $A$ to its Jordan normal form,

[TABLE]

The columns of $F$ are the (generalized) eigenvectors of $A$,

[TABLE]

If $(x,y,z)$ satisfies (4.9), then $(u,v,w)^{T}=F^{-1}(x,y,z)^{T}$ satisfies

[TABLE]

where we note that the variables $u$ and $v$ here are distinct from those used in the unfolding analysis of the previous section, we have introduced the abbreviation $\ell=1+q+q^{2}$, and $n(u,v,w,\tilde{p},\tilde{r},s)$ is

[TABLE]

By standard center manifold theory [6, 20, 30], there exists a family of two-dimensional center manifolds, $W^{c}(0)$, which for any $k>0$ are given by

[TABLE]

The center manifolds are not unique, but they are $C^{k}$ equivalent. The function $h$ satisfies the invariance equation,

[TABLE]

with $h=0$, $\frac{\partial h}{\partial u}=0$, and $\frac{\partial h}{\partial v}=0$ at $(0,0,\tilde{p},\tilde{r})$. Also, $h$ may be represented by a series of the form

[TABLE]

where $h_{0}$ is a constant and $h_{i}$ ($i=1,2,\ldots$) is a homogeneous polynomial function of degree $i$ of the variables $u$, $v$, $\tilde{p},\tilde{r}$, and the coefficients in these polynomials depend on $q$ and $s$. The first two terms vanish identically since the center manifolds are tangent to the center subspace at the origin. The expressions for $h_{2}$ and $h_{3}$ are given in Appendix B.

Remark 4.3**.**

It is useful to compare the expression (4.17) for $h(u,v,\tilde{p},\tilde{r})$ with the expression (3.2) for $h_{\varepsilon}(x,y,p,r)$ for the slow–fast system. Using the transformation $(x,y,z)^{T}=F(u,v,w)^{T}$ and the slow manifold expansion (3.2) with terms up to and including $\mathcal{O}(\varepsilon^{3})$, one finds that

[TABLE]

where the quantity $u+(1-\varepsilon)v+h_{\varepsilon}$ is the representation of the slow manifold in the $(u,v,w)$ coordinates. Thus, the slow manifold and the center manifold are $\mathcal{O}(\varepsilon^{4})$ close as $\varepsilon\to 0$ ($q\to\infty$), which is an important consistency check. This shows that the center manifolds naturally continue to the slow manifolds from the regime $q\gg 1$.

On the center manifold, the dynamics are governed by

[TABLE]

where $n$ is given by (4.14). The center manifold $W^{c}(0)$ and the reduced equations (4.18) on it are illustrated in Figure 4.3, with $s=0$ and $(p,r)$ chosen such that $P_{0}$ is an unstable focus enclosed by a stable limit cycle. Figure 4.3 also shows how a solution through a nearby initial condition rapidly approaches $W^{c}(0)$ and then winds towards the stable limit cycle.

The system (4.18) has a BT at $Q_{0}$ corresponding to the equilibrium $P_{0}$, and the unfolding of the BT point is similar to that in Section 3.3 (details not shown). Also, as shown in Figures 4.1 and 4.2, the bifurcation curves that emanate from the BT point on $W^{c}(0)$ are of the same type as those that emanate from $Q_{0}$ on the slow manifold.

The smoothness of the center manifolds is determined by the Lyapunov type numbers that measure the growth rates of solutions in the directions tangential to the center manifolds and the growth rates normal to the manifolds [16]. For each $(p,r,s)$, there is a critical value of $q$, which we label $q_{c}(p,r,s)$, given by the largest value of $q>1$ for which both Lyapunov type numbers are less than one. The center manifolds are at least $C^{1}$-smooth for all $q>q_{c}(p,r,s)$.

Numerically, the Lyapunov type numbers can be approximated from the eigenvalues at the equilibria of (4.1) and the Floquet multipliers of the limit cycles of (4.1). Let $\sigma_{ij}(p,r,s)$, $j=1,2$, denote the eigenvalues at the equilibrium $P_{i}$, $i=0,1,2$, corresponding to the two eigenvectors in the tangent plane to the manifolds. For $s=0$ and for each $(p,r)$ in an array of values above the diagonal, we calculated $M(p,r)=\max_{i=0,1,2;j=1,2}|\operatorname{Re}(\sigma_{i,j}(p,r,0))|$. This quantity is a function of $q$, and we determined the smallest value of $q$ for which $M(p,r)<|\lambda_{3}|$, where we recall that $\lambda_{3}=-(1+q-(q/1+q))$ is the transverse eigenvalue at $Q_{0}$. This smallest value marks the value at which the tangential growth rate at the equilibria on $W^{c}(0)$ first equals the normal growth rate. This enables a simple numerical approximation of $q_{c}$, since for all $q$ greater than this value, the tangential growth rate is less than the normal growth rate, as approximated by $\lambda_{3}$, and hence shows that the reduced equations (4.18) are an approximation to (4.1). Numerically, for $(p,r)\in\left\{0<p\leq 2,0<r\leq\tfrac{3}{2}\right\}$, we find that $q_{c}(p,r,0)$ is less than or equal to one. For $p$ fixed at $1$, we find that $q_{c}(p,r,0)$ increases monotonically with $r$ (with $q_{c}(1,2,0)\approx 1.4$). Similarly, $q_{c}$ increases monotonically with $r$ for other fixed values of $p$. We also computed the Floquet multipliers of the limit cycles on $W^{c}(0)$ for sample points $(p,r)$ in the same array, and verified that the maximal tangential growth rates associated to the limit cycles are less than $|\lambda_{3}|$ for $q>q_{c}(p,r)$.

Remark 4.4**.**

The cubic approximation of $W^{c}(0)$ is valid in a neighborhood of $P_{0}$. Any extraneous equilibria generated by it lie outside the domain of validity.

5 Discussion

In this article, we presented a dynamical systems analysis of the rich behavior exhibited by (1.1), which was proposed by Maasch and Saltzman in [32] in their study of the glacial cycles observed in the climate record of the Pleistocene Epoch. We identified the regimes in which the MS model (1.1) exhibits limit cycles, and we determined the locations of the various bifurcation curves along which the limit cycles are created and disappear.

Central to understanding the limit cycles and their bifurcations is the result that the long-term system dynamics of the third-order MS model (1.1) occur on (and near) two-dimensional invariant manifolds for most values of the parameter $q$, which measures the characteristic time scales for the total global ice mass and the volume of North Atlantic Deep Water. First, by considering the regime in which $q$ is asymptotically large ($q\gg 1$), we showed in Sections 2 and 3 that the model possesses two-dimensional invariant slow manifolds. All initial conditions relax quickly to these manifolds, and the solutions approach the stable equilibria and stable limit cycles along the manifolds. These slow manifold results generalize the earlier analysis presented in [15] for the second-order system obtained by formally setting $q=\infty$ and $s=0$ in (1.1).

Second, by considering finite—but not large—values of $q$, we showed in Section 4 that the model has a family of two-dimensional invariant center manifolds for all finite values of $q>1$. These center manifolds also contain the equilibria and limit cycles, attract all nearby initial conditions, and govern how solutions approach the stable states. Moreover, the center manifolds are at least $C^{1}$-smooth for all $q$ greater than a critical value $q_{c}$. We note that the center manifolds for finite $q$ smoothly transition to the slow manifolds as $q$ becomes asymptotically large, as is consistent with the theory that slow manifolds may be viewed as special types of center manifolds [6, 17].

On both the slow and center invariant manifolds, there are Bogdanov-Takens points which act as organizing centers of the dynamics. These were first studied in the symmetric case in which the parameter $s$ was set to zero. There is a single $\mathbb{Z}_{2}$-symmetric BT point from which the three main types of bifurcation curves emanate: the Hopf bifurcation curves along which the stable and unstable limit cycles are created from the equilibria, the homoclinic bifurcation curve along which the pair of small-amplitude unstable limit cycles coalesce to form one large-amplitude unstable limit cycle, and the curve of saddle-node bifurcations of limit cycles along which the large-amplitude stable and unstable limit cycles disappear. Then, it was shown that there is a symmetry-breaking which creates two non-degenerate BT points in the physically-relevant case of $s>0$, where the limit cycles of (1.1) exhibit slow glaciation and rapid deglaciation. The bifurcation curves emanating from these two organizing centers are similar to, but more complex than, those in the symmetric case.

A summary of the local and global bifurcation curves in which the limit cycles are created and disappear is given by Figures 4.1 and 4.2. Stable limit cycles of the full MS model (1.1) are found in the region bounded by the curve e0 of supercritical Hopf bifurcations of $P_{0}$ and the curve of homoclinic bifurcations (solid blue branch emanating from $Q_{0}$ to the right). This same region is also shown in Figure 5.1 now with the isoperiod curves (alternating pink and cyan curves) for the stable limit cycles super-imposed. Maasch and Saltzman [32] focused on the parameter region near the isoperiod curve corresponding to 100 Kyr cycles.

Besides being of intrinsic interest, the results presented herein for the internal dynamics of the MS model (1.1) will be useful for studying the effects of orbital (Milankovitch) forcing and slow variation of the system parameters in (1.1). In [32], Maasch and Saltzman showed that the model can be tuned to exhibit the 40 Kyr cycles of the early Pleistocene under orbital (Milankovitch) forcing, and that slow passage of the system parameters $p$ and $r$ through Hopf bifurcations offers a mechanism to understand the mid-Pleistocene transition. We are presently pursuing these findings.

Appendix A Persistence of the Homoclinics to $P_{1}$

In this appendix, we present the Melnikov analysis for the existence and persistence of the single- and double-loop homoclinics to the saddle $P_{1}$ for the slow–fast asymmetric system. We consider (3.9) in the regions III and IV, where $-\tfrac{1}{4}\delta^{2}<\mu<0$. Without loss of generality, we set $\mu=-1$. Then (3.9) possesses three equilibria, located at

[TABLE]

In regions III and IV, the $x$-coordinates of $P_{1}$ and $P_{2}$ satisfy $x_{2}^{*}<x_{1}^{*}<0$. Thus, we identify $(u_{1}^{*},v_{1}^{*})$ with $P_{1}$ and $(u_{2}^{*},v_{2}^{*})$ with $P_{2}$. We also note that $\delta>2$ in this region (for $\mu=-1$).

First, we translate $P_{1}$ to the origin via the coordinate transformation

[TABLE]

After dropping the overlines, the translated system is

[TABLE]

where

[TABLE]

The Hamiltonian of the unperturbed version ($\eta=0$) of (A.1) is $H(u,v)=\tfrac{1}{2}v^{2}-\tfrac{1}{2}\nu u^{2}+\tfrac{1}{3}\kappa u^{3}+\tfrac{1}{4}u^{4}$, and the level set $H=0$ corresponds to a pair of homoclinic orbits, $\Gamma_{1}^{\pm}$,

[TABLE]

where $\alpha=\left(\kappa^{2}+\tfrac{9}{2}\nu\right)^{-1/2}$.

The splitting distance is obtained as an asymptotic expansion in $\eta$,

[TABLE]

The bifurcation equation $D(\lambda)=0$ gives the set of $\lambda$ for which the homoclinic orbits persist under perturbations. Thus, we have

[TABLE]

where the integrals $I_{0}^{\pm},I_{1}^{\pm},$ and $I_{2}^{\pm}$ are defined by

[TABLE]

where now $\psi=\sqrt{\frac{2\kappa^{2}}{9\nu}}$ and $\phi=\sqrt{\frac{2\kappa^{2}+9\nu}{9\nu}}$ (and $\phi^{2}-\psi^{2}=1$). Figure A.1 shows the simple zeros of $D(\lambda)$.

This demonstrates analytically the existence of the curves of homoclinic bifurcations to $P_{1}$. Here, $\lambda_{1}^{+}$ corresponds to the right-homoclinics (Figure 3.3(b)) which lie along the (flat) branch emanating from $Q_{0}$ to the right and touches the shifted diagonal, and $\lambda_{1}^{-}$ corresponds to the left-homoclinics (Figure 3.3(e)) which lie along the dashed branch that emanates from $Q_{1}$ to the right, below the diagonal (cf. Figure 3.6).

The large-amplitude homoclinics to $P_{1}$ for $\mu<0$ are obtained by measuring the splitting distance along the concatenation, $\Gamma_{1}^{-}\cup\Gamma_{1}^{+}$, of the left and right homoclinics. In this case, the bifurcation equation $D(\lambda)=0$ yields

[TABLE]

as the values of $\lambda$ for which the large-amplitude homoclinics to $P_{1}$ persist. The function $\lambda_{1}^{\cup}(\delta)$ is shown in green in Figure A.1.

The definition (3.8) of $\mu$ and $\lambda$ generates a linear relation between $p$ and $r$, namely

[TABLE]

This formula yields the slope of the tangent line to the curves of homoclinic bifurcations at $Q_{1}$. With $\mu=-1$, this analysis also agrees with the numerical continuation results.

Appendix B Center Manifold Reduction

In this appendix, we present the terms in the series (4.17) for the function $h$, whose graph represents the center manifold $W^{c}(0)$ of the Maasch–Saltzman model (4.1). Recall that $h=h_{2}+h_{3}+\cdots$, where $h_{2}$ and $h_{3}$ are homogeneous polynomials of degree 2 and 3, respectively, and the coefficients depend on $q$ and $s$,

[TABLE]

The unknown coefficients are found by substituting (B.1) in the invariance equation (4.16) and equating coefficients of like monomials. One can show by induction on the power that quadratic and cubic terms dependent only on $\tilde{p}$ and $\tilde{r}$ must be zero. To show that there are no quadratic terms of this type, we examine the invariance equation (4.16) and observe that (i) $h_{1}=0$ implies that the function $n$ in the right member of the invariance equation cannot generate such a quadratic term, and (ii) the terms in the left member cannot generate such a quadratic term either. One can then show similarly that there are no cubic terms of this type either. Hence, $b_{8},b_{9},b_{10}=0$, and $c_{17},c_{18},c_{19},c_{20}=0$. The nonzero terms are

[TABLE]

and

[TABLE]

These coefficients are used in the series representation (4.17) of $h(u,v,{\tilde{p}},{\tilde{r}})$ in Section 4.2.

Finally, we observe that the MS model (4.1) is symmetric under the reflection

[TABLE]

and all terms in the center manifold expansion respect this symmetry. However, we recall that only the regime $s>0$ is of physical relevance in order to model the asymmetry between rapid deglaciation and slow glaciation.

Acknowledgments

The authors thank Edgar Knobloch for a useful conversation.