Computing connecting orbits to infinity
associated with a homoclinic flip bifurcation

Homoclinic flip bifurcations are bifurcations of codimension two that occur in families of continuous-time dynamical systems, given by ODEs or vector fields, whose phase space is of dimension at least three. This type of bifurcation of a homoclinic orbit to a real hyperbolic saddle—a special trajectory that converges both in forward and backward time to the saddle equilibrium—occurs when a stable or unstable manifold transitions, when followed along the homoclinic orbit, from being orientable to being non-orientable, or vice versa. While such a change of orientability may occur in higher-dimensional phase spaces, the characterization of homoclinic flip bifurcations and their unfoldings have been studied in detail mostly for the lowest-dimensional case of a three-dimensional systems, both from a theoretical [10, 18, 19, 20, 22, 31] and a numerical point of view [1, 8, 15, 16, 21].

In three-dimensions, which is the case we also consider here, the orientability of a homoclinic orbit is determined by the orientablility of the two-dimensional (un)stable manifold. The saddle equilibrium is assumed to be a hyperbolic saddle, meaning that it has one or two stable and two or one unstable eigenvalues, respectively. In the case of one stable eigenvalue, which we encounter in the example vector field below, its stable manifold is one dimensional, that is, a curve consisting of two trajectories that converge to the saddle in forward time; its unstable manifold is two dimensional, that is, a surface formed by all trajectories that converge to the saddle in backward time. Generically, this surface, when followed locally along the homoclinic orbit in backward time to the equilibrium, closes up along the one-dimensional strong unstable manifold, which is tangent to the strongest unstable eigendirection of the saddle, to form either a cylinder in the orientable case, or a Möbius strip in the non-orientable case. The orientability of the homoclinic orbit, that is, of the two-dimensional unstable manifold (in this case), can change in three different ways:

• 1.

  the two unstable eigenvalues become complex conjugate and the equilibrium turns into a saddle-focus;

• 2.

  orbit flip: the one-dimensional stable manifold returns (in backward time) to the equilibrium tangent to the strong unstable eigendirection (instead of the weakest unstable eigendirection);

• 3.

  inclination flip: the two-dimensional unstable manifold when followed along the homoclinic orbit is tangent to the plane spanned by the stable and weak unstable eigendirections (instead of the plane spanned by the stable and strong unstable eigendirections).
  

The first case, when the saddle equilibrium has a double leading eigenvalue, is known as a Belyakov bifurcation [6, 7] and a numerical study of its simplest unfolding was performed in [8, 25]; see also [2]. Both the orbit flip and inclination flip bifurcations have similar unfoldings, which are split into three different generic cases, referred to as $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$, depending on the eigenvalues of the saddle equilibrium; see [18, 19, 20, 21] for the actual eigenvalue conditions. In a two-parameter unfolding, case $\mathbf{C}$, which is the most complicated one, gives rise to infinitely many curves of secondary bifurcations, including saddle-node, period-doubling, and $n$-homoclinic bifurcations (which involve so-called $n$-homoclinic orbits that make $n-1$ close passes of the equilibrium before returning to it). Moreover, there are two different unfoldings with quite different arrangements of the associated secondary bifurcations, called inward twist $\mathbf{C}_{\rm in}$ and outward twist $\mathbf{C}_{\rm out}$; which of the two unfoldings occurs is determined by global geometric properties of the two-dimensional manifold [10, 18].

We previously conducted numerical studies of the unfoldings of the different homoclinic flip bifurcations, with a particular focus on clarifying changes of two-dimensional global invariant manifolds [1, 15, 16]. To this end, we studied a model vector field developed by Sandstede [32, 33]: a system of three ordinary differential equations with eight parameters, which contains all three cases $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$ of both orbit flip and inclination flip bifurcations for suitable choices of the parameters; the underlying homoclinic orbit is always to the saddle located at the origin. However, all unfoldings of case $\mathbf{C}$ in Sandstede’s model are outward twisted for the inclination flip [33]. Furthermore, we considered the case of the orbit flip in Sandstede’s model and did not find a parameter regime where the unfolding of case $\mathbf{C}$ is inward twisted. In fact, no explicit example of a vector field with the inward-twisted case $\mathbf{C}_{\rm in}$ of a flip bifurcation was known.

This has changed very recently, when Algaba, Domínguez-Moreno, Merino, and Rodríguez-Luis [3] found an example of a three-dimensional quadratic system with an inward-twisted homoclinic flip bifurcation. More precisely, they presented the system

$$\left\{\begin{array}[]{rcrcl}\dot{x}_{1}&=&a&+&x_{2}\,x_{3},\\[2.84526pt] \dot{x}_{2}&=&-x_{2}&+&x_{1}^{2},\\[2.84526pt] \dot{x}_{3}&=&b&-&4x_{1}\end{array}\right.$$

(1)

and showed that it exhibits a codimension-two homoclinic orbit flip bifurcation of type $\mathbf{C}_{\rm in}$ of a saddle $p$ when $a\approx-1.20338$ and $b\approx 1.89616$. This was achieved by identifying the orbit flip homoclinic bifurcation numerically in a parameter regime where the eigenvalue condition at $p$ of case $\mathbf{C}$ is satisfied, and then computing a sufficient number of secondary bifurcation curves emanating from this codimension-two point to show that it unfolds as case $\mathbf{C}_{\rm in}$. Note that the homoclinic orbit is not to the origin but to the equilibrium $\mathbf{p}=(x_{1},x_{2},x_{3})=(b/4,\,b^{2}/16,\,-16\,a/b^{2})$, which exists provided $b\neq 0$.

Algaba et al. [3] studied the local bifurcation structure near $\mathbf{C}_{\rm in}$ in quite some detail. We are interested here in how the unfolding of $\mathbf{C}_{\rm in}$ is embedded more globally in an overall bifurcation diagram. An interesting aspect of system (1) is that it has only one finite equilibrium, the origin that is involved in the homoclinic bifurcation. By contrast, in Sandstede’s model there exists a second equilibrium, and we found that it is responsible for additional global bifurcations in the overall bifurcation diagram, including connecting orbits to the origin [1, 15, 16].

In this paper, we focus on a particular global feature of system (1), namely connecting orbits to a second equilibrium $\mathbf{q}_{\infty}$ that, intriguingly, is located at infinity. More specifically, we study the bifurcation diagram near $\mathbf{C}_{\rm in}$ in a suitable two-parameter plane and show that it features curves of $n$-homoclinic bifurcation that emanate from $\mathbf{C}_{\rm in}$. We find that these curves accumulate, as $n$ increases, on a curve of heteroclinic bifurcations involving infinity, given by the existence of a connecting orbit of codimension-one from the finite equilibrium $\mathbf{p}$ to the equilibrium $\mathbf{q}_{\infty}$ at infinity. Hence, the bifurcation diagram near the codimension-two orbit flip point $\mathbf{C}_{\rm in}$ in the quadratic system (1) features global connecting orbits to infinity.

To address the challenge of finding such connecting orbits to infinity, we adapt the numerical technique from [23], referred to as Lin’s method, for computing connecting orbits between finite objects. More precisely, we modify system (1) by translating the equilibrium $\mathbf{p}$ to the origin $\mathbf{0}$ and by introducing a third parameter that helps separate the very closely spaced bifurcations. For the transformed system, we perform a weighted directional compactification of phase space to study the behavior at infinity. The analysis at infinity involves an additional blow-up transformation to understand the behavior of solutions approaching $\mathbf{q}_{\infty}$ in backward time; these are bounded in the blow-up chart by a two-dimensional surface related to a specific periodic orbit at infinity. To set up Lin’s method, we choose a section $\Sigma$ that is well defined in the original coordinates as well as the blow-up chart near infinity. We then consider and compute two orbit segments, from this periodic orbit surrounding $\mathbf{q}_{\infty}$ to $\Sigma$ and from $\Sigma$ to $\mathbf{0}$, such that their end points in $\Sigma$ lie in the so-called Lin space. In this way, we obtain a well-defined and computable test function, which is zero exactly at the parameter values where there exists a connecting orbit to $\mathbf{q}_{\infty}$. All our computations are performed via the continuation of solutions to suitable two-point boundary value problems with the pseudo-arclength continuation package Auto [11, 12] and the homoclinic continuation toolbox HomCont [9].

This paper is organized as follows. In the next section, we introduce the transformed system with a homoclinic orbit to the origin. Furthermore, we identify the codimension-two point $\mathbf{C}_{\rm in}$ and present a bifurcation diagram in two parameters that suggests the need for the analysis of the dynamics at infinity. Section 3 presents the compactification and the blow-up analysis in different charts at infinity. We use these results in Section 4 to set up Lin’s methods by defining a suitable boundary value problem to compute the boundary of the existence of connecting orbit from $\mathbf{0}$ to the equilibrium $\mathbf{q}_{\infty}$ at infinity. Section 5 then explains how this set-up can also be used to find connecting orbits from a saddle periodic orbit to infinity. In the final Section 6 we draw conclusions and point to some directions for further research.

A homoclinic flip bifurcation of case $\mathbf{C}$ is the global bifurcation of the lowest codimension that involves a real saddle equilibrium (its eigenvalues relevant to the bifurcation are all real) and gives rise to chaotic dynamics. While its complete unfolding is not fully understood, a lot is known about the dynamics nearby [18, 19, 20, 22, 31, 16, 21]. It has been proven that there exists a nearby parameter region with Smale-horseshoe dynamics, and this means that infinitely many saddle periodic orbits are created near this codimension-two point. The precise way in which this occurs is organized by cascades of period-doubling and saddle-node bifucations as well as cascades of $n$-homoclinic bifurcations; these infinitely many different bifurcations occur arbitrarily close in parameter space to the homoclinic flip bifurcation point. The difference between the two cases $\mathbf{C}_{\rm our}$ and $\mathbf{C}_{\rm in}$ lies in the positions of these cascades relative to the primary homoclinic orbit that undergoes the flip bifurcation.

Algaba et al. [3] identified an orbit flip bifurcation of system (1), and computed and presented several bifurcation curves in the $(a,b)$-parameter plane to show that the bifurcation diagram is that of inward-twisted type $\mathbf{C}_{\rm in}$. Unfortunately, the bifurcations for system (1) occur extremely close together and it is not easy to distinguish them. Furthermore, the multi-loop periodic orbits that are created in the $n$-homoclinic bifurcations come very close to the saddle equilibrium and do not extend far in phase space. In a bid to ameliorate this, we move the unique equilibrium $\mathbf{p}$ of (1) to the origin and introduce a third parameter to obtain the system

$$\left\{\begin{array}[]{rcrcccl}\dot{x}&=&\alpha\,y&+&\gamma\,z&+&y\,z,\\ \dot{y}&=&\beta\,x&-&y&+&x^{2},\\ \dot{z}&=&-4\,x.&&&&\end{array}\right.$$

(2)

Here, the new variables are given by $(x,y,z)=(x_{1}-b/4,\,x_{2}-b^{2}/16,\,x_{3}+16\,a/b^{2})$, and we consider system (2) as a new system with three independent parameters. Note that the previous system (1) is recovered for the special choice of the new parameters $\alpha=-16a/b^{2},\,\beta=b/2\text{ and }\gamma=b^{2}/16.$ The advantage of having the specific parameter $\gamma$ is that it allows us to improve the separation of bifurcating periodic orbits from $\mathbf{0}$, the only equilibrium of system (2). For the purpose of this paper, the parameters $\alpha$ and $\beta$ are allowed to vary as the unfolding parameters of the orbit flip bifurcation, while we fix $\gamma=0.5$ throughout our investigation.

System (2) is our object of study. To find the orbit flip bifurcation in the $(\alpha,\beta)$-plane for $\gamma=0.5$, we start from the parameter values corresponding to those reported in [3] and continue the (primary) homoclinic bifurcation to $\gamma=0.5$. Next, we continue the locus of the homoclinic bifurcation as a curve in the $(\alpha,\beta)$-plane while keeping $\gamma=0.5$ fixed throughout all subsequent computations. On the curve of homoclinic bifurcations, we detect the orbit flip point, which we denote $\mathbf{C}_{\rm in}$, at $(\alpha,\beta)\approx(5.3573,2.19173)$. At this parameter point the origin $\mathbf{0}$ has eigenvalues $\lambda^{s}\approx-3.7444$, $\lambda^{u}\approx 0.2108$, and $\lambda^{uu}\approx 2.5335$. Hence, at $\mathbf{C}_{\rm in}$, and also nearby, the point $\mathbf{0}$ is a hyperbolic saddle with a one-dimensional stable manifold $W^{s}(\mathbf{0})$ and a two-dimensional unstable manifold $W^{u}(\mathbf{0})$. Moreover, the condition $\mid\!-\lambda^{uu}\!\mid<-\lambda^{s}$ on the eigenvalues for an orbit flip of case $\mathbf{C}$ is indeed satisfied at $\mathbf{C}_{\rm in}$ [21, 32].

Figure 1 shows the partial bifurcation diagram for system (2), which provides the numerical evidence that we are indeed dealing with an orbit flip of inward-twisted type $\mathbf{C}_{\rm in}$. The curve of (primary) homoclinic bifurcation in the $(\alpha,\beta)$-plane is separated by the orbit flip point $\mathbf{C}_{\rm in}$ into a branch $\mathbf{H_{o}}$ of orientable and a branch $\mathbf{H_{t}}$ of non-orientable or twisted homoclinic bifurcation. Subsequently, we found and continued other bifurcation curves emanating from $\mathbf{C}_{\rm in}$, namely, curves $\mathbf{SNP}$ of saddle-node bifurcation of periodic orbits (green), $\mathbf{PD}$ and $\mathbf{PD^{2}}$ of period-doubling bifurcation, and $\mathbf{H^{n}}$ of $n$-homoclinic bifurcation for $n=2,\dots,6$. Figure 1(a) shows the bifurcation diagram in the $(\alpha,\beta)$-plane of (2). Because the different bifurcation curves are still a bit hard to distinguish in the $(\alpha,\beta)$-plane, panel (b) shows them relative to the curve $\mathbf{H_{o/t}}$ of primary homoclinic bifurcation. More specifically, we show the $(\alpha,\hat{\beta})$-plane, where $\hat{\beta}$ represents the distance to $\mathbf{H_{o/t}}$ with respect to the $\beta$-coordinate. Hence, the curve $\mathbf{H_{o/t}}$ is now the $\alpha$-axis where $\hat{\beta}=0$. Figure 1(b) illustrates that all bifurcation curves emanate from the point $\mathbf{C}_{\rm in}$ on the side of $\mathbf{H_{o}}$; in particular, the curves $\mathbf{H^{n}}$ of $n$-homoclinic bifurcation are tangent to $\mathbf{H_{o}}$ near $\mathbf{C}_{\rm in}$, as can be seen in panel (c). Moreover, the curve $\mathbf{SNP}$, as well as the first two curves $\mathbf{PD}$ and $\mathbf{PD^{2}}$ of a cascade of period-doubling bifurcations lie on one side of $\mathbf{H_{o/t}}$, while the curves $\mathbf{H^{n}}$ lie on the other side. These are all characteristic features that distinguish the inward twist from the outward twist [16, 32, 33]. Hence, we conclude that the codimension-two point $\mathbf{C}_{\rm in}$ of (2) for $\gamma=0.5$ is of the same inward-twisted type as that of (1) found in [3].

Figure 2 illustrates the transition through the orbit flip bifurcation along the curve $\mathbf{H_{o/t}}$. On both sides of $\mathbf{C}_{\rm in}$, the one-dimensional stable manifold $W^{s}(\mathbf{0})$ returns to $\mathbf{0}$ tangent to the weak stable eigendirection to form the homoclinic orbit $\mathbf{\Gamma_{\rm HOM}}$. At the same time, the two-dimensional unstable manifold $W^{u}(\mathbf{0})$ returns back to the saddle $\mathbf{0}$ and closes up along the one-dimensional strong unstable manifold $W^{uu}(\mathbf{0})\subset W^{u}(\mathbf{0})$. The shown part of the surface consists of a family of orbit segments that start at distance of $10^{-3}$ from $\mathbf{0}$; it has been computed with the boundary-value problem set-up from [15, 16, 24]. The two typical cases of homoclinic bifurcation are that $W^{u}(\mathbf{0})$ either forms a cylinder along $\mathbf{H_{o}}$ or a Möbius strip along $\mathbf{H_{t}}$, depending on which side of $W^{uu}(\mathbf{0})$ the stable manifold $W^{s}(\mathbf{0})$ returns. This is illustrated in Fig. 2 by the different positions on the surface $W^{u}(\mathbf{0})$ of the curve $W^{uu}(\mathbf{0})$ relative to the homoclinic orbit; see especially the enlargements. The change in orientability occurs at the point $\mathbf{C}_{\rm in}$ when $W^{s}(\mathbf{0})$ returns to $\mathbf{0}$ exactly along $W^{uu}(\mathbf{0})$, which is represented in Fig. 2 by the respective branches of the two manifolds coinciding in panel $\mathbf{C}_{\rm in}$. As a result, the surface $W^{u}(\mathbf{0})$ comes back tangent to the strong direction and so is neither orientable nor non-orientable.

The top-left region of the bifurcation diagram in Fig. 1(b), to the left of $\mathbf{SNP}$ and above $\mathbf{H_{t}}$, is the only region where system (2) has no periodic orbits as a result of the flip bifurcation. Upon crossing $\mathbf{H_{t}}$, a single saddle periodic orbit $\Gamma_{t}$ is created, which is non-orientable; hence, it has negative nontrivial Floquet multipliers. When followed around the point $\mathbf{C}_{\rm in}$, the periodic orbit $\Gamma_{t}$ persists throughout the different regions in the bifurcation diagram until the curve $\mathbf{PD}$, where it merges with a repelling period-doubled orbit in a subcritical period-doubling bifurcation. This turns $\Gamma_{t}$ into an attracting periodic orbit, which exists in the region between the curves $\mathbf{PD}$ and $\mathbf{SNP}$. Since $\Gamma_{t}$ is now attracting, it can transform from a non-orientable to an orientable periodic orbit, which allows it to bifurcate at $\mathbf{SNP}$ with the orientable saddle periodic orbit $\Gamma_{o}$ that is created upon crossing $\mathbf{H_{o}}$ into the region with $\hat{\beta}>0$.

Many more periodic orbits are created and disappear again near the orbit flip point $\mathbf{C}_{\rm in}$, and we now turn our attention to an associated global feature of the bifurcation diagram: the nature of the curves $\mathbf{H^{n}}$ of $n$-homoclinic bifurcations. Observe in Fig. 1(b) that each of the curves $\mathbf{H^{2}}$ to $\mathbf{H^{6}}$ emanating from $\mathbf{C}_{\rm in}$ has a fold (a maximum) with respect to $\alpha$ and then extends towards decreasing $\alpha$ and $\hat{\beta}$, past the $\alpha$-value of the point $\mathbf{C}_{\rm in}$. Hence, all these curves also exist on the side of $\mathbf{H_{t}}$. The curve $\mathbf{PD^{2}}$ emanating from $\mathbf{C}_{\rm in}$ ends on the curve $\mathbf{H^{2}}$ at a codimension-two orbit flip bifurcation point $\mathbf{C^{2}_{\rm O}}$, quite close to the fold. We find that the bifurcation diagram in the $(\alpha,\beta)$-plane is even more complicated than was suggested in [3]. We identify codimension-two inclination flip bifurcation points $\mathbf{C^{n}_{\rm I}}$ on each of the curves $\mathbf{H^{2}}$ to $\mathbf{H^{6}}$, again very close to where they have a fold with respect to $\alpha$; see the enlargement Fig. 1(c). Also shown in all panels is the curve $\mathbf{SNP^{3}}$ of saddle-node bifurcation of periodic orbits that emanates from $\mathbf{C^{3}_{\rm I}}$. We observe that for sufficiently small values of $\alpha$ the $n$-homoclinic orbits along the curves $\mathbf{H^{3}}$ to $\mathbf{H^{6}}$ are non-orientable.

The computed curves $\mathbf{H^{2}}$ to $\mathbf{H^{6}}$ in Fig. 1 suggest that they are part of a family of curves $\mathbf{H^{n}}$ that accumulate on a well-defined limiting curve. Therefore, we now focus on the limiting behavior of the curves $\mathbf{H^{n}}$ and of the associated $n$-homoclinic orbits as the number of loops $n$ increases. The homoclinic orbits on $\mathbf{H_{t}}$ and $\mathbf{H^{2}}$ to $\mathbf{H^{6}}$ for fixed $\alpha=5.3$ are shown in Fig. 3, where they are assigned the same color as the corresponding curves in Fig. 1. Each time, from one panel to the next, the branch of $W^{s}(\mathbf{0})$ that forms the homoclinc orbit has one extra loop before closing up. Notice that, with increasing $n$, the additional loops of the homoclinic orbit extend increasingly futher along the $y$-direction. This behavior is intriguing, because it suggest that the $n$-homoclinic orbits converge with $n$ to a heteroclinic connection from $\mathbf{0}$ to an equilibrium or periodic orbit at infinity, which corresponds to the limiting case of an infinite number of larger and larger loops. This suggests that the curves $\mathbf{H^{n}}$ in the two-parameter plane accumulate onto a curve of such heteroclinic bifurcations involving infinity, which is, therefore, expected to be of codimension one.

For the purpose of finding a possible heteroclinic bifurcation involving infinity, we must identify equilibria or periodic orbits at infinity. We take advantage of the fact that system (2) is a polynomial vector field, which means that we can compactify the phase space. In general terms, the behavior at infinity is given, after a suitable compactification, by the terms of highest order. We identify and analyze different invariant objects in new coordinate charts that represent the dynamics at and near infinity. This approach makes it possible to continue equilibria or other special solutions as they interact in degenerate bifurcations at infinity [15]. The purpose here is to use charts at infinity to set up a well-posed boundary value problem with a solution that represents the heteroclinic connection to infinity.

More specifically, we follow the recent work by Matsue [28] to obtain a suitable Poincaré compactification for system (2); see also [17, 29, 30]. The underlying idea was already proposed in [13] for planar vector fields, where it is defined as a directional blow-up for so-called quasi-homogeneous vector fields. In our context, this means applying a directional compactification in the direction of positive $y$, because the $n$-homoclinic orbits extend predominantly in the $y$-direction as $n$ increases, while their $x$- and $z$-components remain relatively bounded.

Note that system (2) is not quasi-homogeneous. However, investigation of the leading terms of the right-hand side of system (2) in the limit to infinity shows that it is asymptotically quasi-homogeneous [28] to the quasi-homogeneous vector field of type $(3,4,1)$ and order $3$ given by

$$\left\{\begin{array}[]{rcl}\dot{x}&=&y\,z,\\ \dot{y}&=&x^{2},\\ \dot{z}&=&-4x.\end{array}\right.$$

(3)

The powers of the directional blow-up are then determined by the type of the quasi-homogeneous system (3), which leads to the coordinate transformation

$$(x,y,z)\mapsto(\bar{x},\bar{z},\bar{w}),\quad x=\frac{\bar{x}}{\bar{w}^{3}},\;y=\frac{1}{\bar{w}^{4}},\mbox{ and }z=\frac{\bar{z}}{\bar{w}}.$$

These coordinates define the chart with $y>0$, and $\bar{w}$ represents the distance to infinity in the $y$-direction. More precisely, let $(x_{\rm s},y_{\rm s},z_{\rm s})$ be the transformed coordinates of system (3) inside the Poincaré sphere centered at the origin, where directions of escape to infinity are represented by points on the sphere of radius one. In these coordinates, $(\bar{x},\bar{y},\bar{z})$ correspond to the projection of the positive $y_{\rm s}$-hemisphere of the two-dimensional Poincaré sphere onto the plane defined by $y_{\rm s}=1$. The resulting weighted directional compactification then becomes

$$\left\{\begin{array}[]{rcrl}\dot{\bar{x}}&=&\dfrac{1}{\bar{w}^{2}}&\left(\bar{z}+\alpha\,\bar{w}+\frac{3}{4}\bar{x}\,[\bar{w}^{2}-\beta\,\bar{w}^{3}\,\bar{x}-\bar{x}^{2}]+\gamma\,\bar{z}\,\bar{w}^{4}\right),\\[8.53581pt] \dot{\bar{z}}&=&\dfrac{1}{\bar{w}^{2}}&\left(-4\bar{x}+\frac{1}{4}\bar{z}\,[\bar{w}^{2}-\beta\,\bar{w}^{3}\,\bar{x}-\bar{x}^{2}]\right),\\[8.53581pt] \dot{\bar{w}}&=&\dfrac{1}{\bar{w}^{2}}&\left(\frac{1}{4}\bar{w}\,[\bar{w}^{2}-\beta\,\bar{w}^{3}\,\bar{x}-\bar{x}^{2}]\right).\end{array}\right.$$

It can be desingularized via a rescaling of time with the factor $\bar{w}^{2}$, yielding the desingularized vector field that contains the dynamics at infinity as

$$\left\{\begin{array}[]{rcl}\dot{\bar{x}}&=&\bar{z}+\alpha\,\bar{w}+\frac{3}{4}\bar{x}\,(\bar{w}^{2}-\beta\,\bar{w}^{3}\,\bar{x}-\bar{x}^{2})+\gamma\,\bar{z}\,\bar{w}^{4},\\[5.69054pt] \dot{\bar{z}}&=&-4\bar{x}+\frac{1}{4}\bar{z}\,(\bar{w}^{2}-\beta\,\bar{w}^{3}\,\bar{x}-\bar{x}^{2}),\\[5.69054pt] \dot{\bar{w}}&=&\frac{1}{4}\bar{w}\,(\bar{w}^{2}-\beta\,\bar{w}^{3}\,\bar{x}-\bar{x}^{2}).\end{array}\right.$$

(4)

We are now ready to analyze the dynamics at infinity and decide whether it contains equilibria or periodic orbits that could be involved in a suspected heteroclinic connection. To this end, we set $\bar{w}=0$ in system (4) and observe that the $(\bar{x},\bar{z})$-plane is indeed invariant. The resulting system

$$\left\{\begin{array}[]{rcrcl}\dot{\bar{x}}&=&\bar{z}&-&\frac{3}{4}\bar{x}^{3},\\[5.69054pt] \dot{\bar{z}}&=&-4\bar{x}&-&\frac{1}{4}\bar{x}^{2}\,\bar{z},\end{array}\right.$$

(5)

has a single equilibrium at $(\bar{x},\bar{z})=(0,0)$, which is, in fact, not hyperbolic. This equilibrium is the equilibrium $\mathbf{q}_{\infty}$ at infinity in system (2). To understand the dynamics at infinity, that is, on the $(\bar{x},\bar{z})$-plane, we convert to polar coordinates. More precisely, we consider the ellipsoidal transformation

$$(\bar{x},\bar{z})\mapsto(\bar{r},\bar{\theta}),\quad\bar{x}=\bar{r}\,\cos{\bar{\theta}},\mbox{ and }\bar{z}=2\bar{r}\,\sin{\bar{\theta}}.$$

The vector field in these ellipsoidal polar coordinates becomes

$$\left\{\begin{array}[]{rcl}\dot{\bar{r}}&=&-\frac{1}{4}\bar{r}^{3}\cos^{2}{\bar{\theta}}\left(2+\cos{2\bar{\theta}}\right),\\[5.69054pt] \dot{\bar{\theta}}&=&-2+\tfrac{1}{2}\bar{r}^{2}\cos^{3}{\bar{\theta}}\sin{\bar{\theta}}.\end{array}\right.$$

Note that $\dot{\bar{r}}<0$ for all $(\bar{r},\bar{\theta})$ with $\bar{r}>0$, and that $\dot{\bar{\theta}}<0$ and close to $-2$ as soon as $\bar{r}$ is small enough. Hence, all trajectories in the $(\bar{x},\bar{z})$-plane converge to $\mathbf{q}_{\infty}$, which lies at the origin in this planar coordinate system; moreover, locally near $\mathbf{q}_{\infty}$, trajectories will spiral clockwise towards it. This behavior is illustrated in Fig. 4, where we plot several trajectories in the $(\bar{x},\bar{z})$-plane in panel (a) and project them back onto the Poincaré sphere in panel (b); note that system (5) only describes the dynamics in the chart with $y_{\rm s}>0$, and only the corresponding half-sphere is shown.

In the full three-dimensional blown-up system (4), the point $\mathbf{q}_{\infty}$ at $(\bar{x},\bar{z},\bar{w})=(0,0,0)$ is not a hyperbolic attractor. Therefore, we perform an additional $\bar{w}$-directional blow-up by applying the transformation

$$(\bar{x},\bar{z},\bar{w})\mapsto(x_{\rm B},z_{\rm B},w_{\rm B}),\quad\bar{x}=x_{\rm B}\,w_{\rm B},\;\bar{z}=z_{\rm B}\,w_{\rm B},\mbox{ and }\bar{w}=w_{\rm B},$$

to system (4). This gives the vector field

$$\left\{\begin{array}[]{rcl}\dot{x}_{\rm B}&=&\alpha+z_{\rm B}+\gamma\,w_{\rm B}^{4}\,z_{\rm B}+\frac{1}{2}x_{\rm B}\,w_{\rm B}^{2}\,(1-\beta\,x_{\rm B}\,w_{\rm B}^{2}-x_{\rm B}^{2}),\\[5.69054pt] \dot{z}_{\rm B}&=&-4x_{\rm B},\\[5.69054pt] \dot{w}_{\rm B}&=&\frac{1}{4}w_{\rm B}^{3}\,(1-\beta x_{\rm B}\,w_{\rm B}^{2}-x_{\rm B}^{2}),\end{array}\right.$$

(6)

which further characterizes the dynamics at infinity on a local half-sphere around $\mathbf{q}_{\infty}$. Setting $w_{\rm B}=0$ in system (6), we find that the invariant $(x_{\rm B},z_{\rm B})$-plane is foliated by ellipses of the form $4x^{2}_{\rm B}+(z_{\rm B}+\alpha)^{2}=c^{2}$; see Fig. 5(a). The trajectories in the $(x_{\rm B},z_{\rm B})$-plane correspond to trajectories on the blown-up half-sphere with $\bar{w}>0$ centered at $(\bar{x},\bar{z},\bar{w})=(0,0,0)$. Figure 5(b) gives an impression of how the previously identified dynamics at infinity interacts with the blown-up half-sphere in the $(x_{\rm B},z_{\rm B},w_{\rm B})$-space.

The next step is to determine the property of system (6) for $\bar{w}>0$. First, we resort to numerical simulation and determine how initial conditions with $\bar{w}>0$ approach the $(x_{\rm B},z_{\rm B})$-plane. Figure 6 shows that there are two types of behavior. Panel (a) shows two trajectories of system (6) for $(\alpha,\beta)=(5.3,2.0)$, obtained by integration in both forward and backward time from the initial conditions $(x_{\rm B},z_{\rm B},w_{\rm B})=(1,-\alpha,0.05)$ and $(x_{\rm B},z_{\rm B},w_{\rm B})=(1.3,-\alpha,0.05)$, respectively. The former initial condition leads to a trajectory (orange) that converges in backward time in a spiraling fashion to the equilibrium $(x_{\rm B},z_{\rm B},w_{\rm B})=(0,-\alpha,0)$ of (6). The other trajectory (blue) first approaches the $(x_{\rm B},z_{\rm B})$-plane in backward time but then diverges away from it; in particular, it does not reach the equilibrium $(x_{\rm B},z_{\rm B},w_{\rm B})=(0,-\alpha,0)$. This is illustrated further in Fig. 6(b) in a local cross-section defined by $z_{\rm B}=-\alpha$. Notice that the two trajectories are very close together before they separate in backward time at about $\bar{w}=0.08$.

We conclude that there exists an invariant critical surface $S_{\rm c}$ that separates the two qualitatively different regions in phase space where trajectories converge to $(x_{\rm B},z_{\rm B},w_{\rm B})=(0,-\alpha,0)$ and where they do not. Furthermore, Fig. 6 suggests that this difference in the backward-time limit of trajectories is entirely due to the fact that ellipses near $(x_{\rm B},z_{\rm B},w_{\rm B})=(0,-\alpha,0)$ are repelling in the $\bar{w}$-direction, while beyond some distance, they are attracting in the $\bar{w}$-direction. The surface $S_{\rm c}$ is associated with the critical ellipse in the $(x_{\rm B},z_{\rm B})$-plane that is neither repelling nor attracting in the $\bar{w}$-direction, and which goes through the point $(x_{\rm B},z_{\rm B})\approx(1.1547,-\alpha)$ (magenta). Our computations indicate that the critical surface $S_{\rm c}$ is effectively a straight elliptical cylinder when $\bar{w}$ is small.

Based on these careful observations, we approximate $S_{\rm c}$ as a straight $\bar{w}$-cylinder around the $w_{B}$-axis through the point $(0,-\alpha,0)$. We denote this cylinder $C_{r^{*}}$, with a specific radius $r^{*}$, and require that

$$\int_{C_{r}}{X_{\rm B}\cdot\hat{n}_{C_{r}}\;dC_{r}}=0$$

(7)

be satisfied for $r=r^{*}$, where $\hat{n}_{C_{r}}$ is the direction normal to $C_{r}$. We are using the average zero-flux condition (7) to define $C_{r^{*}}$ because, in general, there is no cylinder that is invariant under the vector field $X_{\rm B}$ defined by (6). To find $r^{*}$ we transform system (6) to cylindrical coordinates by

$$(x_{\rm B},z_{\rm B},w_{\rm B})\mapsto(r_{\rm B},\theta_{\rm B},w_{\rm B}),\quad x_{\rm B}=r_{\rm B}\,\cos{\theta_{\rm B}}\mbox{ and }z_{\rm B}=2r_{\rm B}\,\sin{\theta_{\rm B}}-\alpha.$$

The integral can then be evaluated in a straightforward way as:

	$\displaystyle\int_{C_{r}}{X_{\rm B}\cdot\hat{n}_{C_{r}}\;dC_{r}}$	$\displaystyle=$	$\displaystyle\left.\int dw_{\rm B}\;\int_{0}^{2\pi}\dot{r}_{\rm B}\;d\theta_{\rm B}\right|_{r_{\rm B}=r}$
		$\displaystyle=$	$\displaystyle\int\frac{1}{8}\pi r\,w_{\rm B}^{2}\,(4-3r^{2})\;dw_{\rm B}=\frac{1}{24}\pi r\,w_{\rm B}^{3}\,(4-3r^{2}).$

Hence, there are two zeros of the zero-flux condition (7), namely, $r=0$ and $r=r^{*}=\frac{2}{3}\sqrt{3}$. Note that $r=0$ corresponds to the $w_{\rm B}$-axis through the equilibrium at $(0,-\alpha,0)$. We conclude that the critical cylinder $C_{r^{*}}$ with $r^{*}=\frac{2}{3}\sqrt{3}\approx 1.1547$ is the local approximation of the separating invariant surface $S_{\rm c}$. This value agrees with our numerical simulations and $C_{r^{*}}$ is a good first-order approximation of $S_{\rm c}$.

Recall that the $(x_{\rm B},z_{\rm B},w_{\rm B})$-coordinate system of (6) corresponds to a directional blow-up of the equilibrium $\mathbf{q}_{\infty}$ at infinity in the original coordinates, which corresponds to the origin in the $(\bar{x},\bar{z},\bar{w})$-coordinates of the desingularized system (4). Figure 7 illustrates in two ways the separatrix $S_{\rm c}$ (magenta surface) represented by the inverse image of the critical cylinder $C_{r^{*}}$ under the respective coordinate transformations. Panel (a) shows how $S_{\rm c}$ emanates from a corresponding periodic orbit on the blown-up half-sphere centered at the origin of the $(\bar{x},\bar{z},\bar{w})$-space. However, periodic orbits only exist on the blown-up (half-)sphere and not in the $(\bar{x},\bar{z},\bar{w})$-space itself. Deflating the blown-up sphere back to the origin, the local approximation $C_{r^{*}}$ of $S_{\rm c}$ is the cone emanating from the origin in the $(\bar{x},\bar{z},\bar{w})$-space that is shown in Fig. 7(b).

All trajectories inside the separatrix $S_{\rm c}$ converge, in backward time, to $\mathbf{q}_{\infty}$, which is the origin in the $(\bar{x},\bar{z},\bar{w})$-space. Hence, $S_{\rm c}$ acts as a kind of two-dimensional unstable manifold of the non-hyperbolic point $\mathbf{q}_{\infty}$ at infinity. For special choices of the parameters $\alpha$ and $\beta$ in system (2), the one-dimensional stable manifold $W^{s}(\mathbf{0})$ of the origin in the original $(x,y,z)$-coordinates lies in the surface $S_{\rm c}$. We refer to this well-defined phenomenon of codimension one as a heteroclinic connection between $\mathbf{0}$ and $\mathbf{q}_{\infty}$, and we denote it by $\mathbf{Het^{\infty}}$. It is our hypothesis that the curves $\mathbf{H^{n}}$ of $n$-homoclinic orbits, which have increasingly longer excursions towards infinity, accumulate in the $(\alpha,\beta)$-plane on the corresponding curve $\mathbf{Het^{\infty}}$; see Fig. 3.

Hence, the task is to find the heteroclinic connection $\mathbf{Het^{\infty}}$ and to continue it in the $(\alpha,\beta)$-plane. To this end, we employ the approach known as Lin’s method [23, 26] to set up a two-point boundary value problem (BVP) for two orbit segments such that their concatenation is the sought-after connecting orbit in $W^{s}(\mathbf{0})\cap S_{\rm c}$. The essence of Lin’s method is to choose a codimension-one plane $\Sigma$ that separates the two invariant objects involved, here $\mathbf{0}$ and $\mathbf{q}_{\infty}$, and to consider an orbit segment in $W^{s}(\mathbf{0})$ up to $\Sigma$ and an orbit segment in $S_{\rm c}$ up to $\Sigma$. For parameters that are not at the bifurcation value, these two orbit segments exhibit a gap in $\Sigma$. Lin’s theorem states that this orbit pair and, hence, the gap are uniquely determined when the difference between their end points in $\Sigma$ is constrained to lie in a fixed subspace called the Lin space [26]. The associated signed Lin gap in the Lin space is then a well-defined test function with zeros that correspond to connecting orbits; such zeros can be found via the continuation of the corresponding orbit segements as solutions of an overall BVP [23]. Once a zero is found, the associated connecting orbit can be followed in system parameters.

The challenge, here, is that one of the equilibria lies at infinity and we have an approximation for $S_{\rm c}$ in blown-up coordinates. Note that systems (2) and (6) are homeomorphic in the open sets where they coincide [28]. This allows us to define $\Sigma$ with respect to both coordinate systems. We then consider one orbit segment that is a solution of system (2) with one end point near the saddle $\mathbf{0}$ and lying in its stable eigenspace (which is the linear approximation of $W^{s}(\mathbf{0})$) and the other lying in $\Sigma$; and a second orbit segment that is a solution of system (6) with one end point near the point $(\bar{x},\bar{y},\bar{z})=(0,0,0)$ representing $\mathbf{q}_{\infty}$ lying in the linear approximation $C_{r}$ of $S_{\rm c}$ and the other lying in $\Sigma$. The respective coordinate transformations allow us to ‘glue’ the original $(x,y,z)$-coordinates of system (2) to the $(x_{\rm B},z_{\rm B},w_{\rm B})$-coordinates of the blown-up system (6), so that we can define and determine the Lin gap.

We use this adapted Lin’s method to find an initial connecting orbit in $W^{s}(\mathbf{0})\cap S_{\rm c}$, along with the relevant bifurcation value for $\beta$, where we keep $\alpha=5.3$ fixed. We define

$$\Sigma:=\{(x,y,z)\;|\;x=0\}\simeq\{(x_{\rm B},z_{\rm B},w_{\rm B})\;|\;x_{\rm B}=0\},$$

which is a suitable choice that works in both coordinate systems for $\bar{w}\neq 0$ because $x=\bar{x}/\bar{w}^{3}$ and $\bar{x}=x_{\rm B}\,w_{\rm B}$ with $w_{\rm B}=\bar{w}$.

To define the orbit segment $\mathbf{u}$ in $(x,y,z)$-coordinates that lies in $W^{s}(\mathbf{0})$ up to $\Sigma$ we define the BVP

	$\displaystyle\dot{\mathbf{u}}$	$\displaystyle=$	$\displaystyle T_{0}\,X(\mathbf{u}),$		(8)
	$\displaystyle\mathbf{u}(1)$	$\displaystyle=$	$\displaystyle\delta_{0}\,\mathbf{e}^{s},$		(9)
	$\displaystyle\mathbf{u}(0)^{*}\,\mathbf{n}$	$\displaystyle=$	$\displaystyle 0.$		(10)

Here, $X$ denotes the vector field (2) and $T_{0}$ is the total integration time between the first and last point on the orbit segment; it enters (8) in explicit form so that the orbit segment $\mathbf{u}(t)$ is defined for $t\in[0,1]$. Boundary condition (9) requires that the end point $\mathbf{u}(1)$ lies at a small distance $\delta_{0}$ from the saddle $\mathbf{0}$ along its stable eigenvector $\mathbf{e}^{s}$ (which has been normalized to have length 1). This ensures that $\mathbf{u}(1)$ lies in $W^{s}(\mathbf{0})$ to good approximation, provided $\delta_{0}$ is sufficiently small; we fix $\delta_{0}=10^{-4}$ as an appropriate value throughout. Finally, the dot product in boundary condition (10) involves the unit vector $\mathbf{n}=(1,0,0)$ normal to $\Sigma$, which ensures that the start point $\mathbf{u}(0)$ lies in $\Sigma$. We remark that the stable eigenvector $\mathbf{e}^{s}$ in  (9) needs to be continued as well when system parameters are changed; we achieve this by solving the BVP of the corresponding stable eigenvector problem [23] together with (8)–(10).

Similarly, the orbit segment $\mathbf{u}$ in $(x_{\rm B},z_{\rm B},w_{\rm B})$-coordinates in $S_{\rm c}$ up to $\Sigma$ is defined by the BVP

	$\displaystyle\dot{\mathbf{u}}_{\rm B}$	$\displaystyle=$	$\displaystyle T_{\rm B}\,X_{\rm B}(\mathbf{u}_{\rm B}),$		(11)
	$\displaystyle\mathbf{u}_{\rm B}(0)$	$\displaystyle=$	$\displaystyle(\tfrac{2}{3}\sqrt{3}\,\cos{\theta_{\rm B}},\,\tfrac{4}{3}\sqrt{3}\,\sin{\theta_{\rm B}}-\alpha,\,\delta_{\rm B}),$		(12)
	$\displaystyle\mathbf{u}_{\rm B}(1)^{*}\,\mathbf{n}$	$\displaystyle=$	$\displaystyle 0.$		(13)

In (11) the vector field (6) is denoted $X_{\rm B}$, and $T_{\rm B}$ is the total integration time. Boundary condition (12) requires that the start point $\mathbf{u}_{\rm B}(0)$ lies in the cylinder $C_{r^{*}}$, which has been parameterized by the angle $\theta_{\rm B}\in[0,\,2\pi]$ and the distance $\delta_{\rm B}$ in the $w_{\rm B}$-direction; we set $\delta_{\rm B}=0.1$ throughout. Boundary condition (13) again ensures that the end point $\mathbf{u}_{\rm B}(1)$ lies in $\Sigma$, because $\mathbf{n}=(1,0,0)$ is also the unit normal to $\Sigma$ in $(x_{\rm B},z_{\rm B},w_{\rm B})$-coordinates.

To find first orbit segments $\mathbf{u}$ and $\mathbf{u}_{\rm B}$ that satisfy (8)–(10) and (11)–(13), respectively, we fix $\beta=1.8$ and proceed as follows (recall that $\alpha=5.3$ and $\gamma=0.5$ are fixed). For $\mathbf{u}$ we require initially only (8) and (9), and start a continuation in the integration time $T_{0}$ from $T_{0}=0$; note that this constitutes solving the initial value problem from the point $\delta_{0}\,\mathbf{e}^{s}$ by continuation. During this computation we monitor the dot product and record whenever $\mathbf{u}$ satisfies condition (10), that is, $\mathbf{u}(0)$ lies in $\Sigma$. Similarly, for $\mathbf{u}_{\rm B}$ we require only (11) and (12); we start with $\theta_{\rm B}=0$ and continue in $T_{\rm B}$ from $T_{\rm B}=0$, while recording whenever (13) is satisfied and $\mathbf{u}_{\rm B}(1)$ lies in $\Sigma$. We remark that both conditions (10) and (13) are satisfied for many values of $T_{0}$ and $T_{B}$, respectively, because the trajectories that contain $\mathbf{u}$ and $\mathbf{u}_{B}$ intersect $\Sigma$ many times.

We choose orbit segments $\mathbf{u}$ and $\mathbf{u}_{\rm B}$ that have end points in $\Sigma$ which lie suitably close to each other and couple them by defining the Lin space and associated Lin gap. To this end, we define and then fix the unit vector

$$\mathbf{\Psi}:=\frac{\mathbf{u}(0)-\widetilde{\mathbf{u}}_{\rm B}(1)}{\mid\!\mid\!\mathbf{u}(0)-\widetilde{\mathbf{u}}_{\rm B}(1)\!\mid\!\mid}$$

given by the initial chosen end points of $\mathbf{u}$ and $\mathbf{u}_{\rm B}$; here $\widetilde{\mathbf{u}}_{\rm B}(1)$ is the end point of $\mathbf{u}_{\rm B}(1)$ in the original $(x,y,z)$-ccordinates of the section $\Sigma$. The vector $\mathbf{\Psi}$ is generically transverse to $S_{\rm c}\cap\Sigma$, spans the Lin space $Z$, and defines the Lin gap $\eta$ via the boundary condition

$$\widetilde{\mathbf{u}}_{\rm B}(1)=\mathbf{u}(0)+\eta\,\mathbf{\Psi}.$$

(14)

Note that the new parameter $\eta$ is the signed distance between the two end points of the orbit segments along the Lin space $Z\subset\Sigma$, which is fixed once chosen in this way.

We now consider the combined boundary value problem given by (8)–(12) and (14), which is automatically satisfied by the chosen orbit segments $\mathbf{u}$ and $\mathbf{u}_{\rm B}$ and uniquely defines the Lin gap $\eta$. When $\mathbf{u}$ and $\mathbf{u}_{\rm B}$ are continued in $\beta$, where $\theta_{\rm B}\in[0,\,2\pi]$, $T_{0}>0$, $T_{\rm B}>0$, and $\eta\in\mathbb{R}$ are free parameters (but crucially $Z\subset\Sigma$ remains fixed), the Lin gap $\eta$ is monitored. When $\beta$ changes, the orbit segment $\mathbf{u}$ as well as the $\theta$-dependent orbit segment $\mathbf{u}_{\rm B}$ vary. In light of the Lin condition (14), the angle parameter $\theta_{\rm B}$ is adjusted automatically in such a way that the end point $\mathbf{u}_{\rm B}(1)$ only varies along the direction $\mathbf{\Psi}$, either away from or towards $\mathbf{u}(0)$. When a zero of $\eta$ is detected then we have found the value of $\beta$ at which the heteroclinic connection $\mathbf{Het^{\infty}}$ occurs; the corresponding heteroclinic orbit that connects $\mathbf{q}_{\infty}$ with $\mathbf{0}$ is given as the concatenation of $\mathbf{u}$ and $\mathbf{u}_{\rm B}$.

Figure 8 illustrates the set-up with Lin’s method, shown in projection onto compactified Poincaré coordinates that represent $\mathbb{R}^{3}$ inside the unit sphere (not shown) centered at the origin $\mathbf{0}$. The plane in Fig. 8 is the common Lin section $\Sigma$ defined by $x=x_{\rm B}=0$. Notice that the chosen orbit segment $\mathbf{u}$ intersects $\Sigma$ three times, that is, we choose to work with the third intersection of the trajectory from $\mathbf{0}$. Similarly, the chosen orbit segment $\mathbf{u}_{\rm B}$ intersects $\Sigma$ many times. The orbit segment $\mathbf{u}_{\rm B}$ in Fig. 8 was chosen so that its end point $\mathbf{u}_{\rm B}(1)$ in $\Sigma$ is sufficiently close to the end point $\mathbf{u}(0)$. The Lin space $Z\subset\Sigma$, which appears curved in the compactified Poincaré coordinates of Fig. 8, remains fixed during the subsequent continuation of the BVP (8)–(12) and (14) in $\beta$. Panel (b) shows the situation when the Lin gap $\eta$ has been closed and the connecting orbit found as the concatenation of $\mathbf{u}$ and $\mathbf{u}_{\rm B}$. As is seen in Fig. 8, the orbit segment $\mathbf{u}_{\rm B}$ intersects $\Sigma$ multiple times. We remark that, from a practical perspective, it is best to choose $\mathbf{u}_{\rm B}(1)$ close to $\mathbf{u}(0)$. On the other hand, choosing any of the earlier intersection points of $\mathbf{u}_{\rm B}$ in the numerical set-up will result in the same connecting orbit, provided that a zero of the Lin gap is found.

As soon as a heteroclinic connection $\mathbf{Het^{\infty}}$ is detected as a zero of $\eta$, it can be continued with the BVP (8)–(12) and (14) in $\alpha$ and $\beta$, where $\theta_{\rm B}\in[0,\,2\pi]$, $T_{0}>0$, and $T_{\rm B}>0$ are free parameters but $\eta=0$ is now kept fixed. This continuation leads to the curve $\mathbf{Het^{\infty}}$ in the $(\alpha,\beta)$-plane that is shown in Fig. 9 together with the other curves of the bifurcation diagram from Fig. 1. As panels (a) and (b) of Fig. 9 show, the curve $\mathbf{Het^{\infty}}$ has the same general shape as the curves $\mathbf{H^{n}}$ of $n$-homoclinic bifurcation (shades of cyan) for $n=2,3,\dots,6$: it also emanates from the codimension-two flip bifurcation point $\mathbf{C}_{\rm in}$, has monotonically decreasing $\hat{\beta}$ and has a fold for a very similar value of $\hat{\beta}$. Indeed, we conclude from Fig. 1 that the curves $\mathbf{H^{n}}$ accumulate on the curve $\mathbf{Het^{\infty}}$ as $n$ tends to infinity.

Panel $\mathbf{Het^{\infty}}$ of Fig. 9 illustrates that the heteroclinic connection from $\mathbf{0}$ to $\mathbf{q}_{\infty}$ is characterized by the one-dimensional manifold $W^{s}(\mathbf{0})$ spiraling away (in backward time) from $\mathbf{0}$ towards infinity to approach $\mathbf{q}_{\infty}$ along the cone/cylinder $S_{\rm c}$. Indeed, this is the limiting case between the two generic situations that are illustrated in Fig. 10. Either $W^{s}(\mathbf{0})$ lies outside $S_{\rm c}$ and does not reach $\mathbf{q}_{\infty}$, as in panel (a), or it lies inside $S_{\rm c}$ and spirals onto $\mathbf{q}_{\infty}$, as in panel (c). The former situation occurs to the left of the curve $\mathbf{Het^{\infty}}$ in the $(\alpha,\hat{\beta})$-plane of Fig. 9, while $W^{s}(\mathbf{0})$ connects generically to $\mathbf{q}_{\infty}$ to the right of $\mathbf{Het^{\infty}}$.

The Lin’s method set-up from the previous section can be adapted to compute other types of connecting orbits to infinity. We demonstrate this here with the example of a heteroclinic connection from the orientable saddle periodic orbit $\Gamma_{o}$, which bifurcates from the curve $\mathbf{H_{o}}$ and exists for $\hat{\beta}>0$, to the point $\mathbf{q}_{\infty}$. More specifically, we compute an orbit in the intersection set $W^{s}(\mathbf{\Gamma_{o}})\cap S_{\rm c}$, which exists generically, because $W^{s}(\mathbf{\Gamma_{o}})$ and $S_{\rm c}$ are both two dimensional manifolds. As before, we concatenate two orbit segments: $\mathbf{u}$ from a common section $\Sigma$ to $\Gamma_{o}$ and $\mathbf{u}_{\rm B}$ from $\mathbf{q}_{\infty}$ to $\Sigma$, which are again found as solutions to the overall BVP (8)–(12) and (14). The difference is that the vector $\mathbf{e}^{s}$ in boundary condition (9) is now a vector in the stable Floquet bundle of $\Gamma_{o}$. The periodic orbit $\Gamma_{o}$ and its stable Floquet bundle can be computed and continued with the BVP set-up presented in [23], yielding the vector $\mathbf{e}^{s}$ (for any value of the system parameters).

A suitable initial orbit segment $\mathbf{u}$ is found by choosing and fixing $\delta_{0}$ and then, as before, continuing the initial value problem (8) and (9) in the integration time $T_{0}$ from $T_{0}=0$, while recording whenever condition (10) is satified. The initial orbit $\mathbf{u}_{\rm B}$ is found exactly as before, and the vector $\mathbf{\Psi}$, the Lin space $Z$ and the Lin gap $\eta$ are subsequently defined as in Section 4. The overall BVP (8)–(12) and (14) is then automatically satisfied and we use it to continue the two orbit segments $\mathbf{u}$ and $\mathbf{u}_{\rm B}$ to close the Lin gap $\eta$. Because the connecting orbit is generic, the continuation for this problem does not involve a system parameter, but uses the fact that the two-dimensional manifold $W^{s}(\mathbf{\Gamma_{o}})$ is a $\delta_{0}$-family of trajectories. Here, $\theta_{\rm B}\in[0,\,2\pi]$, $T_{0}>0$, $T_{\rm B}>0$, $\eta\in\mathbb{R}$, and the parameter $\delta_{0}$ are free parameters.

Figure 11 illustrates the set-up in compactified Poincaré coordinates; compare with Fig. 8. Panel (a) of Fig. 11 shows the orientable periodic orbit $\Gamma_{o}$, the equilibrium $\mathbf{q}_{\infty}$, the section $\Sigma$ and the initially chosen orbit segments $\mathbf{u}$ and $\mathbf{u}_{\rm B}$ that define the Lin space $Z\subset\Sigma$. The Lin gap $\eta$ is then closed by continuation in $\delta_{0}$, yielding the connecting orbit as the concatenation of $\mathbf{u}$ and $\mathbf{u}_{\rm B}$ as shown in Fig. 11(b); note that the system parameters $\alpha$, $\beta$ and $\gamma$ remain unchanged during this computation.

We studied a quadratic vector field, adapted from that of [3], that exhibits a homoclinic flip bifurcation of the specific inward-twisted type $\mathbf{C}_{\rm in}$. We found that the two-parameter bifurcation diagram near this special point features an accumulation of curves of secondary $n$-homoclinic bifurcations. Numerical evidence that this phenomenon involves an increasing number of loops which move closer to infinity motivated us to set up a numerical scheme based on Lin’s method to find the limiting behavior in the form of a heteroclinic connection to infinity. To this end, the orbit segment in the finite part of phase space was formulated in original coordinates, while the second orbit segment to infinity was defined in different coordinates near infinity. Both are then glued together along the Lin space in a section that is well-defined in both coordinate systems. Closing the Lin gap along the Lin space by continuation of the two coupled orbit segments yielded a first connecting orbit of codimension one between the origin and a point at infinity. A subsequent continuation gave the associated curve in the parameter plane, which was indeed found to act as the accumulation set for the curves of $n$-homoclinic orbits.

Compared to previous uses of a Lin’s method set-up to define suitable boundary value problems for finite connecting orbits, a novel element is the use of blown-up coordinate charts near infinity. Blow-up techniques for polynomial vector fields allow one to study equilibria and other invariant objects at and near infinity. When these are of saddle type in the geometric sense — meaning that they have attracting and repelling directions, but need not be hyperbolic or even semi-hyperbolic — the question arises how they interact with invariant objects in the finite part of phase space, such as equilibria and periodic orbits. Indeed, connections to infinity are a distinct possibility. As we showed, such heteroclinic phenomena involving infinity may provide important information regarding limits of finite global objects.

Our Lin’s method set-up is quite flexible and more widely applicable; this was demonstrated by computing a connecting orbit from a finite saddle periodic orbit to a point at infinity. Hence, it constitutes a new tool for the study of global properties of polynomial vector fields. The system studied here is a case in point, and its further bifurcation analysis is the subject of ongoing research. Note that this quadratic vector field is presently the only system that is known to exhibit a homoclinic flip bifurcation of the inward-twisted type of case $\mathbf{C}$; hence, it has the role of a model vector field for this specific bifurcation, much in the spirit of Sandstede’s model [32, 33] which features effectively all other types of flip bifurcations. The investigation of the outward-twisted type in the latter model shows that a flip bifurcation of case $\mathbf{C}$ gives rise to a very complicated global bifurcation structure. In light of its different local structure, we expect to find a different, yet comparably complicated overall bifurcation structure in the wider vicinity of the orbit flip of the inward-twisted type of case $\mathbf{C}$. Moreover, homoclinic flip bifurcations of all cases have been identified as organizing centers in other vector fields from the literature, specifically in mathematical models of neurons [4, 5, 27]. Their global bifurcation structure may well involve heteroclinic bifurcations with infinity. Hence, we believe that the numerical approach for the identification and continuation of connecting orbits to infinity will have a role to play in their study.

We thank Alejandro Rodríguez-Luis for alerting us to his paper [3] with Algaba, Domínguez and Merino, and for suggesting that we investigate further the inward-twisted orbit flip bifurcation of case C.