1 Introduction

We consider a dynamical system defined in discrete time by a diffeomorphism that maps a given phase space smoothly back to itself in forward and backward time. The diffeomorphism may be known explicitly, as is the case for the well-known Hénon map [26], or be given as the Poincaré return map of a smooth flow. In the latter case, the map must be computed numerically, which is generally done by numerical integration of a point in a chosen Poincaré section (of codimension one) until the trajectory returns to the section. An important role for organising the overall behaviour, especially chaotic dynamics, is played by the stable and unstable invariant manifolds of saddle fixed and periodic points; see Sect. 2 for the definition. For diffeomorphisms of the plane, these are curves and their intersections famously give rise to homoclinic tangle and associated hyperbolic sets, which are saddle Cantor sets with Smale’s horseshoe dynamics [43]. This type of dynamics is an example of uniform hyperbolicity, which is one of the most studied forms of chaotic dynamics introduced by Anosov [2] and Smale [43]. Its characterising feature is that the nonwandering set has a continuous decomposition of its tangent bundle into stable and unstable subspaces and, thereby, admits stable and unstable manifolds.

Uniformly hyperbolic dynamics is structurally stable under \(C^1\)-perturbations, and it was conjectured that this is the default or generic form of \(C^1\)-robust chaotic dynamics [12, 13]. This conjecture is believed to be true for two-dimensional diffeomorphisms [1], because tangencies between one-dimensional stable and unstable invariant manifolds in the plane cannot occur robustly in the \(C^{1}\)-topology [37]. However, in diffeomorphisms of dimensions at least three and, equivalently, flows of dimensions at least four, homoclinic tangencies or heterodimensional cycles—which are specific examples of non-uniformly hyperbolic dynamics—may occur robustly and can, hence, be found in open intervals in parameter space  [12, 40]; see also [4, 22,23,24, 27] for examples. This type of dynamics is also referred to as wild chaos [37].

A heterodimensional cycle is formed by the intersections of the stable and unstable manifolds of two hyperbolic sets—fixed points in the simplest case—that have different unstable dimensions [15, 19]. Each individual heterodimensional cycle is of codimension at least one, which means that it is broken by any generic \(C^1\)-perturbation. Bonatti and Díaz showed that, nevertheless, there exists some heterodimensional cycle for any diffeomorphism in a \(C^1\)-neighbourhood of the initial diffeomorphism [9]. To this end, they introduced the concept of a blender. Blenders play a similar role as the thick horseshoes introduced by Newhouse [39, 40]. A blender is a transitive hyperbolic set \(\Lambda \) for which, say, the k-dimensional stable manifold \(W^s(\Lambda )\) acts as an object of dimension larger than k. In the lowest-dimensional case of a three-dimensional diffeomorphism, which we focus on here, this means that the one-dimensional invariant manifold \(W^s(\Lambda )\) cannot be avoided by one-dimensional curves from (an open neighbourhood of) a certain direction, as is illustrated in [8]. This property is \(C^1\)-robust and implies the existence of robust heterodimensional cycles when, say, the one-dimensional unstable manifold of a hyperbolic set of different unstable dimension intersects \(W^s(\Lambda )\); see also [10, 12, 15].

The blenders defined in [9] are referred to as blender-horseshoes, and they are, essentially, generalisations of Smale’s horseshoe construction to higher dimensions; in particular, they are conjugate to the complete shift on two symbols; see also [10, 11, 17, 18]. Subsequently, other types of blenders were also considered, including symbolic blenders [6], which have central bundles that are not necessarily one dimensional; double-blenders in dimension at least four [5, 7], which consist of two nested blenders in such a way that both their stable and unstable manifolds behave as if they were of higher dimension; and super-blenders, which are particular constructions capable of yielding invariant manifolds of any desired perceived dimension [3, 38, 41].

We consider the simplest case of a blender-horseshoe [12, 16] in \(\mathbb {R}^3\), which is known to exist in the family of Hénon-like diffeomorphisms introduced in [28]; see also [29, 30]. The family is given by an extension to \(\mathbb {R}^3\) of the planar Hénon map [26], which we refer to as h and discuss in detail in Appendix A, via the addition of an affine shear in the third variable. We write this family of diffeomorphisms as

$$\begin{aligned} \mathcal {H}(x,y,z) = (y,\, \alpha - y^2 - \beta x,\, \xi z +y). \end{aligned}$$
(1)

Throughout, we fix \(\alpha =4.2\) and \(\beta =0.3\), in which case h is orientation preserving and has a full Smale horseshoe with dynamics conjugate to a full shift on two symbols on the associated hyperbolic set \(\Lambda _h\). Moreover, we consider the case of expansion in the z-direction, as given by the eigenvalue \(\xi > 1\). Further relevant properties of the family, as needed here, are reviewed in Appendix A.

According to the theory [16], the associated hyperbolic set \(\Lambda \) of \(\mathcal {H}\) is expected to be a blender when \(\xi > 1\) is sufficiently close to 1. For large expansion rates \(\xi \gg 1\), on the other hand, \(\Lambda \) will not be a blender but rather have a Cantor-like structure. These general statements raise the following practical and natural questions:

Q1.:

Over which range of the expansion rate \(\xi \) is the hyperbolic set \(\Lambda \) of the family \(\mathcal {H}\) a blender?

Q2.:

How exactly is the blender lost as \(\xi \) is increased from 1?

The defining property of the blender in this three-dimensional setting is that the one-dimensional stable manifold \(W^s(\Lambda )\) behaves as a surface when viewed from an open set of directions. More precisely, the closure of the projection of \(W^s(\Lambda )\) along a fixed, appropriate direction has non-empty interior; moreover, this property persists in a \(C^1\)-open neighbourhood of such a projection. It is not possible to compute \(W^s(\Lambda )\) in its entirety to check this property, because \(W^s(\Lambda )\) comprises the one-dimensional stable manifolds of all the uncountably many points in \(\Lambda \). However, since periodic points are dense in \(\Lambda \), one may equivalently check that the one-dimensional stable manifold \(W^s(p)\) of any periodic point p has a dense projection. We refer to this characterising property of a blender as the carpet property [28, 29], because the curve \(W^s(p)\) essentially forms a surface near the hyperbolic set \(\Lambda \) by weaving back and forth in phase space.

The carpet property can be investigated numerically by computing a sufficiently long piece of a selected stable manifold. Since such a manifold makes increasingly long excursions away from the bounded set \(\Lambda \), further and further towards infinity and back, such a computation is only feasible after compactification of phase space. This approach was introduced and taken in [28, 29] to produce the first images of blenders, as rendered from the computed stable manifolds \(W^s(p^{\pm })\) of the fixed points \(p^-\) and \(p^+\) of \(\mathcal {H}\) after compactification of \(\mathbb {R}^3\) to a cylinder \(\mathcal {C}\). These computations were performed for different values of \(\xi \) with an adaptation of the manifold growth algorithm from [33].

Fig. 1
figure 1

Phase portrait of \(\mathcal {H}\) for \(\alpha =4.2\), \(\beta =0.3\), and \(\xi =1.2\) in the cylinder \(\mathcal {C}\), shown in projection onto the \((\bar{x},\bar{y})\)-plane in a and in \((\bar{x},\bar{y},\bar{z})\)-space in b. Shown are the fixed points \(p^\pm \) (golden dots), the period-two orbit \(p_2=\{p^1_2, p^2_2\}\) (red dots), and the period-three orbit \(p_3=\{p^1_3, p^2_3, p^3_3\}\) (green dots), as well as the one-dimensional stable manifold \(W^s(p^-)\) (blue curve) and the unstable manifold \(W^u(p^-)\) (red surface), which intersect in the hyperbolic set \(\Lambda \approx W^s(p^-)\cap W^u(p^-)\) (black points). Also shown are the boundary circles of \(\mathcal {C}\) with the sinks \(q^{\pm }\) and the source \(s^+\) (black dots) at infinity

Full size image

To set the stage, Fig. 1 shows a phase portrait with the blender of \(\mathcal {H}\) for \(\xi =1.2\) in the cylinder

$$\begin{aligned} \mathcal {C} = \mathcal {D} \times \mathcal {I}= \left\{ (\bar{x}, \bar{y}, \bar{z}) \in \mathbb {R}^3 \bigm | \Vert (\bar{x}, \bar{y})\Vert \le 1 \text { and } \Vert \bar{z}\Vert \le 1 \right\} , \end{aligned}$$
(2)

which is the relevant compactification of \(\mathbb {R}^3\); see (9) in Appendix A for the coordinate transformation. Here, the carpet property is illustrated with a reasonably long, computed piece of the one-dimensional stable manifold \(W^s(p^-)\) of the fixed point \(p^-\) only. Also shown is a first piece of the two-dimensional unstable manifold \(W^u(p^-)\), which (due to the skew product nature of the map \(\mathcal {H}\)) can be rendered from the data of the corresponding one-dimensional unstable manifold of the Hénon map h as a ruled and flat surface that extends from the bottom to the top of \(\mathcal {C}\), represented by its two boundary circles. The projection in Fig. 1a onto the \((\bar{x},\bar{y})\)-plane of compactified coordinates shows that \(W^s(p^-)\) and \(W^u(p^-)\) intersect transversely everywhere and well inside the bounding circle, confirming that the Hénon map h, indeed, has a full horseshoe. The three-dimensional view in panel (b) in compactified \((\bar{x},\bar{y},\bar{z})\)-space illustrates how this intersection gives rise to the hyperbolic set \(\Lambda \) (black points), which also contains fixed points \(p^\pm \) (golden dots), the period-two orbit \(p_2=\{p^1_2, p^2_2\}\) (red dots), and the period-three orbit \(p_3=\{p^1_3, p^2_3, p^3_3\}\) (green dots) that are shown in Fig. 1. Note that the curve \(W^s(p^-)\) makes repeated excursions to very near a source \(s^+\) at infinity, where it turns sharply to complete what we call a weaving through the central region with the hyperbolic set \(\Lambda \). Similarly, \(W^u(p^-)\) repeatedly approaches the line between two sinks \(q^{\pm }\) at infinity.

The curve \(W^s(p^-)\) shown in Fig. 1b appears to form a surface in certain areas; on the other hand, there are places where one clearly sees strands of the curve with obvious gaps between them. To determine whether a one-dimensional manifold \(W^s(p)\) actually has the carpet property or not requires checking whether such gaps in projection close or remain open as the manifold is computed up to increasingly larger arclengths. In [28, 29], this is done by considering the projection onto an interval of (unordered) intersection points of \(W^s(p)\) with a suitably chosen plane \(\Sigma \). The carpet property is confirmed when the largest gap or gaps between these points of the interval show convergence to zero as the arclength of \(W^s(p)\) is increased. On the other hand, it is not satisfied when the largest gaps reach a constant, non-zero size. By checking this criterion numerically for the manifolds \(W^s(p^\pm )\) of the two fixed points of \(\mathcal {H}\), an estimate was given in [28] of the \(\xi \)-range over which the hyperbolic set \(\Lambda \) is a blender (for the different parameter choices \(\alpha = 9.5\) and \(\beta = -0.1\)). This range turned out to be considerably larger than expected, well beyond the value of \(\xi = 1.185\) for which an analytical proof was provided. A different approach is taken in [30] by constructing a parameter-dependent box for the family \(\mathcal {H}\), successive backward (when \(\xi > 1\)) images of which intersect this initial box in a nested sequence of sub-boxes containing \(\Lambda \) and local segments of \(W^s(\Lambda )\). Tracing the relevant edges of sub-boxes as a function of \(\xi \) provides complementary insight into where the hyperbolic set \(\Lambda \) is a blender. This geometric illustration of the three-dimensional horseshoe nature of \(\mathcal {H}\) also motivated the characterisation in [14] of blenders and heterodimensional cycles via covering relations and cone conditions, which can be verified rigorously with established computer-assisted techniques; in this way, a proof of existence of a blender (for \(\alpha = 9.5\) and \(\beta = -0.1\)) was given also for \(\xi \in [1.01, 1.125]\).

As this discussion shows, the existence of a blender for \(\xi \) not too far from 1 can successfully be established by checking the convergence of gaps numerically, as well as with rigorous methods of proof—be they by hand or computer assisted. However, for increasing values of \(\xi > 1\), near where the carpet property is lost, these approaches run into difficulties. In particular, it is not clear from previous work in [28, 29] how exactly gaps open or close in the transition between blender and non-blender as \(\xi \) is increased or decreased, respectively. In short, previous approaches are not able to answer Q1 fully and are not really suitable for addressing Q2. This realisation motivates the work presented here, which builds on the approach in [28, 29] of considering the intersection set of a given stable manifold \(W^s(p)\) with a section \(\Sigma \). The new aspect is that we present and use a much refined algorithm for computing \(W^s(p)\) that allows us to find accurately a large number of points in \(W^s(p) \cap \Sigma \) that, moreover, are ordered reliably for any \(\xi \) in a parameter interval of interest. This ordering is induced by the fact that we compute \(W^s(p)\) as two branches, which are each parametrised by arclength. Our algorithm is based on ideas from [33], specifically, the efficient representation of a one-dimensional manifold as an arclength-parametrised curve with ordered mesh points that are distributed as determined by curvature. This allows for the control of the error (in the Hausdorff metric) of a computed piece of manifold of some finite arclength, by specifying corresponding accuracy settings. In contrast to this earlier work, which involves adding point after point until a pre-specified arclength is reached, we grow the manifold by taking iterates of an initial (and already reasonably large) fundamental domain. In light of the expansion along the invariant manifold (in forward or backward time), this requires considerable mesh refinement, and we use essentially the same accuracy criteria as in [33] to achieve this by interpolation in the pre-image (rather than the image). This modification of the growth algorithm significantly improves the speed of the computation, and this allows us to compute a given manifold up to a large arclength for many values of \(\xi \). We present our algorithm in Sect. 2 in general terms for an unstable manifold \(W^u(p)\) of a fixed point of some diffeomorphism f on a phase space of dimension at least two; see also Appendix B for information on and reference to the source code and a demo in Matlab. The role of the different accuracy parameters is discussed in Sect. 3 for the illustrative example of the stable manifold \(W^s(p^-)\) as computed with the fourth inverse iterate \(\mathcal {H}^{-4}\).

Our ability to control the accuracy of these manifold computations implies that we can reliably find any pre-specified intersection point, denoted \(w(l) = w(l;\xi )\) with \(\Sigma \) along either of the two branches of \(W^s(p)\) for any given value of \(\xi \); here, we define the intersection index \(l \in \mathbb {Z} \setminus \{0\}\) to represent the |l|th intersection point, where l is taken to be positive along the positive branch and negative along the negative branch of \(W^s(p)\) (as defined appropriately in context). The intersection points are computed accurately via local cubic spline interpolation of the six nearest mesh points representing \(W^s(p)\) locally, three on each side of \(\Sigma \). In this way, we define \(\xi \)-parametrised curves

$$\begin{aligned} \varvec{w}(l) := \left\{ w(l;\xi ) \bigm | \xi \in [\xi _a, \xi _b] \right\} \subset \Sigma \times [\xi _a, \xi _b], \end{aligned}$$

of a given intersection index \(l \in \mathcal {L}\) for a suitable index set \(\mathcal {L} \subset \mathbb {Z} \setminus \{0\}\) and \(\xi \)-range \([\xi _a, \xi _b]\) of interest. In practice, we compute \(W^s(p)\) for a suitably fine mesh of \(\xi \)-values in \([\xi _a, \xi _b]\) and then render the curve \(\varvec{w}(l)\) from the respective |l|th intersection points with \(\Sigma \) on the appropriate branch of \(W^s(p)\).

Fig. 2
figure 2

The computed piece of \(W^s(p^-)\) (blue curve) that generates the first \(2^{11}\) intersection points \(\{w^-(l)\}_{1 \le l \le 2^{11}}\) with \(\Sigma \) (yellow) as a function of \(\xi \). a and b show how \(W^s(p^-)\) intersects \(\Sigma \) for \(\xi =1.2\) and \(\xi =2.4\), respectively, with two views of \((\bar{x},\bar{y},\bar{z})\)-space, the plane \(\Sigma \) with the points \(\{w^-(l)\}\), and their projection onto the \(\bar{z}\)-coordinate. c shows how these intersection points form the set of curves \(\{\varvec{w}^-(l)\}_{1 \le l \le 2^{11}}\) as a function of \(\xi =[1.02,3]\). The (curves of) intersection points are coloured from light to dark blue according to their intersection index l

Full size image

Figure 2 illustrates this approach for the manifold \(W^s(p^-)\) from Fig. 1 and the section \(\Sigma \) defined as the half-plane

$$\begin{aligned} \Sigma := \left\{ (\bar{x}, \bar{y}, \bar{z}) \in \mathcal {C} \bigm | \bar{x} = \bar{y} \text { and } \bar{x} < 0 \right\} , \end{aligned}$$

which we use throughout; note that \(p^- \in \Sigma \). Panels (a) and (b) of Fig. 2 illustrate how \(W^s(p^-)\) intersects \(\Sigma \) for \(\xi =1.2\) and for \(\xi =2.4\), respectively. The fixed point \(p^-\) has positive eigenvalues, and this implies that one branch of \(W^s(p^-)\) immediately goes to infinity and does not intersect \(\Sigma \) at all, while the other branch of \(W^s(p^-)\) weaves back and forth through \(\Sigma \) infinitely many times. We show this branch as the curve computed up to the arclength needed to generate \(2^{11}\) ordered intersection points \(\{w^-(l)\}:= \{w^-(l)\}_{1 \le l \le 2^{11}} \subset W^s(p^-) \cap \Sigma \). The two panels (a) and (b) show \(\Sigma \) in the cylinder \(\mathcal {C}\) with the computed piece of \(W^s(p^-)\) on the left and in a cut-away view in the middle column. The resulting intersection set \(\{w^-(l)\}\) in \(\Sigma \) is shown on the right, and its projection onto the \(\bar{z}\)-axis on the very right; here, the points in \(\{w^-(l)\}\) are coloured with a gradient from light to dark blue as the intersection index \(1 \le l \le 2^{11}\) increases.

Row (a) of Fig. 2 constitutes numerical evidence that \(W^{s}(p^-)\) has the carpet property, meaning that \(\Lambda \) is a blender for \(\xi =1.2\). In particular, the projection of \(\{w^-(l)\}\) onto the z-axis strongly indicates that the limit is an interval with no gaps. Notice further that the shades of blue appear to be quite mixed along this interval. In row (b), on the other hand, \(W^s(p^-)\) does not seem to have the carpet property. In particular, there are now clear gaps in the projection of \(\{w^-(l)\}\) onto the z-axis. This is indicative of a Cantor set construction and strongly suggests that \(\Lambda \) is not a blender for \(\xi =2.4\).

The observations are very similar to, and in agreement with, those in [28, 29]. However, in this previous work, the intersection set \(W^s(p^-) \cap \Sigma \) was not considered as an ordered set, while we are now able to compute them as the set \(\{w^-(l)\}\), with the same order given by \(1 \le l \le 2^{11}\) for any \(\xi \). In other words, the two intersection points for \(\xi =1.2\) and for \(\xi =2.4\) with \(\Sigma \) that have the same shade of blue lie on, and are connected by, one and the same curve \(\varvec{w}^-(l)\) for the respective intersection index l. As Fig. 2c shows, we are able to compute and plot \(2^{11}\) separate curves \(\{\varvec{w}^-(l)\}_{1 \le l \le 2^{11}}\) over the range \(\xi =[1.02,3]\); this was done for a uniform mesh in \(\xi \) with stepsize 0.01. The curves in \(\{\varvec{w}^-(l)\}\) are shown in the relevant projection onto the z-axis, where each curve is coloured according to its intersection index l, in the same way as in panels (a) and (b); the vertical lines indicate the values of \(\xi = 1.2\) and \(\xi = 2.4\) from rows (a) and (b). Also shown in projection onto the \(\bar{z}\)-coordinate are the curves traced out by the fixed points \(p^-\) and \(p^+\), which form the top and the bottom boundaries of the set of curves \(\{\varvec{w}^{-}(l)\}\). As \(\xi \) approaches 1, the finite set of curves \(\{\varvec{w}^-(l)\}_{1 \le l \le 2^{11}}\) appears to approach \(p^-\) and, hence, ‘covers’ an increasingly smaller interval. However, this phenomenon is caused by the weaker and weaker expansion rate in the \(\bar{z}\)-direction, which means that far more intersection points of \(W^s(p^-)\) with \(\Sigma \) would need to be computed to compensate.

The representation in Fig. 2c demonstrates that we are, indeed, able to ‘track’ each ordered intersection point in \(W^s(p^-) \cap \Sigma \) individually as a curve in dependence on the expansion rate \(\xi \). This new capability allows us to determine how their ordering in projection onto the \(\bar{z}\)-coordinate changes, and how they move to close the gaps for decreasing \(\xi \), so that the set of curves \(\{\varvec{w}^{-}(l)\}\) appears to fill a solid region in the \((\xi , \bar{z})\)-plane for about \(1< \xi < 1.6\). As we will demonstrate in Sect. 4, we are able to choose accuracy parameters for our refined method in such a way that, over the same \(\xi \)-range, curves of ordered intersection points with \(\Sigma \) can also be computed reliably for the stable manifolds of the fixed point \(p^+\), the period-two orbit \(p_2\), and the period-three orbit \(p_3\)—in spite of the fact that this requires the second, fourth, and third iterate of the inverse of \(\mathcal {H}\), respectively. Indeed, taking higher iterates makes the computation of an invariant manifold much more challenging, because it means working with larger contraction or expansion rates. Note that the eigenvalues of \(p^+\) and \(p_2\) are negative, so that one needs to consider the double cover, that is, the second iterate of \(\mathcal {H}^{-1}\) and \(\mathcal {H}^{-2}\), respectively.

We conclude in Sect. 5 by plotting the computed manifolds \(W^s(p^-)\), \(W^s(p^+)\), \(W^s(p_2)\), and \(W^s(p_2)\) for \(\xi = 1.2\) all together in the cylinder \(\mathcal {C}\), which provides an unprecedented impression of the properties of a blender. Moreover, our computations of ordered intersection sets of the one-dimensional manifolds of various different fixed and periodic points suggest their rather dynamic organisation in dependence on \(\xi \). We believe this opens the way to the study of how blenders bifurcate, and we point out some ongoing and future research directions in this regard.

2 Formulation of the algorithm

In this section, we consider an arbitrary diffeomorphism

$$\begin{aligned} f : \mathbb {R}^n \rightarrow \mathbb {R}^n \end{aligned}$$

that has a hyperbolic fixed point \(p = f(p) \in \mathbb {R}^n\) with a single real and positive unstable eigenvalue \(\lambda ^u > 1\). According to the Stable Manifold Theorem [42], the unstable manifold

$$\begin{aligned} W^u(p) = \left\{ x\in \mathbb {R}^2 \bigm | f^{-k}(x)\rightarrow p \text { as } k \rightarrow \infty \right\} \end{aligned}$$

is then a one-dimensional smooth curve that is tangent at p to the unstable eigenvector \(\textbf{v}^u\). Similarly, the stable manifold \(W^s(p)\) of p is an \((n-1)\)-dimensional manifold, but we do not consider it here.

As is common in the field, we formulate our algorithm for the computation of \(W^u(p)\) for notational convenience. Since f is a diffeomorphism, any one-dimensional stable manifold can simply be computed as the unstable manifold of the inverse \(f^{-1}\) of f. We also assume for convenience that \(\lambda ^u\) is a positive eigenvalue, which means that f is orientation preserving on \(W^u(p)\) and, hence, consists of two branches (on either side of p) that are each invariant under f. This assumption is again not a restriction: in the orientation-reversing case when \(\lambda ^u < -1\), one considers p as a fixed point under the second iterate \(f^{2}\) with eigenvalue \((\lambda ^u)^2 > 1\) instead. Similarly, formulating the method for a fixed point is not a restriction either, because any periodic point is a fixed point under the appropriate iterate of the map.

Note that \(W^u(p)\) consists of infinitely many orbits that converge to p under iteration of \(f^{-1}\). The smooth curve \(W^u(p)\) may have a very complicated geometry in \(\mathbb {R}^n\), and the task is to compute and represent it in a suitable way. The computation of one-dimensional (un)stable manifolds of maps is already well established, and there are a number of algorithms with various degrees of sophistication; for example, see [33] as an entry point to the literature. All these approaches have in common that one starts from a local approximation of \(W^u(p)\) near p and then ‘grows’ the manifold by considering suitable iterates under f of this initial approximation. The methods differ in whether the linear approximation generated by \(\textbf{v}^u\) is considered, or a higher-order Taylor expansion; whether the data points on \(W^u(p)\) are ordered dynamically, as a family of orbits up to a certain iterate, or geometrically, as sequential points on a curve; and how the accuracy is controlled when obtaining additional points from local interpolation.

For our purposes, we require that each computed branch of \(W^u(p)\) is computed as a curve, represented by an ordered and suitably fine mesh, that is parametrised by arclength from p. This is achieved by the growth method introduced in [33], which grows \(W^u(p)\) point by point subject to accuracy conditions that adapt the distance between mesh points according to the curvature of the manifold. This arclength-based approach has subsequently also been used to compute one-dimensional stable manifolds of maps that do not have a (unique) inverse [20], and it is now implemented in [32] as part of the continuation software package MatContM [36]. This growth method has been extended further (via the definition of suitable two-point boundary value problems) for the computation of one-dimensional (un)stable manifolds of Poincaré maps [21], as well as one-dimensional isochrons (which are invariant manifolds of a fixed-time return map) [25].

The growth method from [33] is not really fast enough for the purpose of producing extremely long pieces of one-dimensional manifolds as is required to check the carpet property of a blender. This was already noted in [28, 29], where an initial and already quite long part of the manifold was computed and then iterated under the map or its inverse. However, this ad-hoc approach led to issues with the ordering of mesh points, but that was not a big deal in this previous work, because the required intersection points were ordered according to the lengths of their successive differences (with respect to some projection). Here, on the other hand, we require very long and consistently ordered meshes to detect the |l|th intersection point reliably along a curve, even for very large |l| and for many points in the specified \(\xi \)-interval. To achieve this, we implement the iteration of a fundamental domain with the curvature-dependent mesh adaptation from [31, 33] to ensure that the resulting ordered mesh still represents a branch of \(W^u(p)\) as an arclength-parametrised curve in a well-resolved way.

To start, we choose a fundamental domain of \(W^u(p)\), that is, a half-open segment \(U_0 \subset W^u(p)\) that intersects every orbit of f on \(W^u(p)\) exactly once; the end points of \(U_0\) satisfy that the point closest to p (when measured with respect to arclength along the manifold) is contained in \(U_0\), and f maps this end point to the other end point of \(U_0\). This means that the respective branch U of \(W^u(p)\) can be represented as the ordered union

$$\begin{aligned} U = \bigcup _{k=-\infty }^{\infty } U_k = \bigcup _{k=-\infty }^{\infty } f^k(U_0). \end{aligned}$$

Note that \(f^k(U_0)\) is also a fundamental domain for any \(k \in \mathbb {Z}\).

In practice, we select \(U_0\) after computing a reasonably long initial piece of the branch \(U \subset W^u(p)\) with the growth method from [33], which starts from the linear approximation given by the eigenvector \(\textbf{v}^u\), very close to the fixed point p. The initial fundamental domain \(U_0\) should be a reasonably long segment that lies already sufficiently far from p, so as to avoid the need for a lot of iterations of \(U_0\) to grow \(W^u(p)\). We then extend or globalise the computed part of the branch U up to a desired arclength by computing a suitable number K of forward images \(U_k = f^k(U_0)\) for \( 1 \le k \le K\).

The key issue for this approach is that f is expanding on \(W^u(p)\), which means that the image of any initially chosen mesh on \(U_0\) under some iterate \(f^k\) is generally not an accurate representation of \(f^k(U_0)\). To guarantee the mesh quality at each step k, we combine the accuracy conditions defined in [31] (developed for planar systems) and [33] and require these to be satisfied for all mesh points on \(f^k(U_0)\). This is achieved by refining the mesh for \(f^{k-1}(U_0)\) until the image of this mesh satisfies the different accuracy conditions as a representation of \(f^k(U_0)\). Crucially, we control the overall interpolation error, and thus, the global computation error, by constructing the refined mesh for \(f^k(U_0)\) from f-images of suitable points on \(f^{k-1}(U_0)\), which ensures the quality of the overall mesh of the computed branch.

2.1 Mesh generation and accuracy conditions

A first part of \(W^u(p)\) is computed with the MatContM implementation [36] of the growth algorithm [33]. The initial fundamental domain \(U_0\) is then chosen to be represented by the ordered list \(M_0= \{q_0,\ldots , q_{m_0}\}\) of mesh points, where \(m_0 \ge 1\) is sufficiently large and \(f(q_{0}) = q_{m_0}\). As mentioned, \(M_0\) should be such that the corresponding segment \(U_0\) has a decent arclength \(\ell _0 > 0\); here, and throughout, this arclength is given (to good approximation) by the sum of the distances between consecutive mesh points in \(M_0\), that is, \(\ell _0 \approx \sum _{j=1}^{m_0} \Vert q_j - q_{j-1}\Vert \).

Suppose now that we computed a piecewise-linear representation of \(U_{k-1} = f^{k-1}(U_0)\) represented by the ordered list \(M_{k-1} = \{q_{m_{k-1}},\ldots , q_{m_k}\}\) that satisfies, at each mesh point, the accuracy conditions discussed below. At step k, we form the mesh \(M^{(0)}_k = f(M_{k-1}):= \{f(q_{m_{k-1}}),\ldots , f(q_{m_k})\}\) and then add and/or remove mesh points in stages producing meshes \(M^{(i)}_k\), until the specified accuracy conditions defined below are met. The new mesh \(M_k\) representing \(U_k\) is then the mesh obtained in the last such stage. Once the desired arclength has been reached, the piece of U in question is represented by the overall mesh \(M = \{q_0, \ldots q_{m_0}, \ldots q_{m_1}, \ldots , q_{m_K}\}\), which is the concatenation of the meshes representing \(U_k\) for \(0 \le k \le K\).

Fig. 3
figure 3

Illustration of the curvature conditions at a mesh point \(f(q_{j})\) with the neighbouring points \(f(q_{j-1})\) and \(f(q_{j+1})\) at distances \(\Delta _{j-1}\) and \(\Delta _{j}\), respectively, that form the angle \(\alpha _j\). Also shown is the image \(f(\hat{q}_{j-1})\) of an interpolated point that decreases \(\Delta _{j-1}\) and/or \(\Delta _{j-1}\alpha _j\)

Full size image

For each stage i, we work with a refined mesh \(\widehat{M}^{(i)}_{k-1}\), such that \(\widehat{M}^{(0)}_{k-1} = M_{k-1}\) and every mesh point in \(M^{(i)}_k\) is of the form \(f(q_{j})\) for some \(q_{j} \in \widehat{M}^{(i)}_{k-1}\). For each such \(f(q_{j}) \in M^{(i)}_k\), we consider the distances \(\Delta _{j-1}\) and \(\Delta _j\) to the previous point \(f(q_{j-1})\) and next point \(f(q_{j+1})\) in \(M^{(i)}_k\), and the angle \(\alpha _j\) between consecutive points \(f(q_{j-1})\), \(f(q_{j})\) and \(f(q_{j+1})\). This setup is illustrated in Fig. 3, and it allows us to check the following criteria:

$$\begin{aligned} \Delta _j< & {} \Delta _{\max },\end{aligned}$$
(3)
$$\begin{aligned} \alpha _j< & {} \alpha _{\max },\end{aligned}$$
(4)
$$\begin{aligned} \Delta _j\alpha _j< & {} (\Delta \alpha )_{\max },\end{aligned}$$
(5)
$$\begin{aligned} \Delta _{j-1}\alpha _{j}< & {} (\Delta \alpha )_{\max },\end{aligned}$$
(6)
$$\begin{aligned} \Delta _{\min }< & {} \Delta _j, \end{aligned}$$
(7)

where \(\Delta _{\max }\), \(\alpha _{\max }\), \((\Delta \alpha )_{\max }\) and \(\Delta _{\min }\) are pre-specified accuracy parameters. Conditions (4), (5), and (7) are those used in [33] for growing the manifold point by point; specifically, (4) ensures an overall bound on the angle, (5) implements the mesh adaptation according to the local curvature at \(f(q_{j})\) (higher curvature of U means larger \(\alpha _j\), hence, forcing smaller \(\Delta _j\); see also [31]), and (7) sets a sufficiently small minimum distance between mesh points to prevent excessive addition of new mesh points during sharp turns of the manifold. Condition (6) is similar to condition (5) and was originally used in [31]; our adaptation of the accuracy conditions considers curvature constraints for both the next and the previous mesh points.

Rather than adding points individually, we check conditions (3)–(6) simultaneously at all mesh points of \(M^{(i)}_{k}\) by making use of vector operations in Matlab. For efficiency purposes, we first check that condition (3) is satisfied, which guarantees a suitable baseline distance of no more that \(\Delta _{\max }\) between all mesh points. Subsequently, we check with (4) that the angles \(\alpha _j\) are sufficiently small at all mesh points. Only then do we check conditions on the curvature, namely, (5) from [33] concerning the distance to the next point, as well as (6) from [31] that involves the distance to the previous point. When any of these conditions is not satisfied, a new mesh point, such as the point \(f(\hat{q}_{j-1})\) in Fig. 3, is added as follows. We determine a refinement \(\widehat{M}^{(i+1)}_{k-1}\) of the mesh \(\widehat{M}^{(i)}_{k-1}\) by adding the point \(\hat{q}_{j-1}\) that we obtain via cubic spline interpolation on \(M_{k-1}\) as the (approximate) midpoint between \(q_{j-1}\) and \(q_j\). Indeed, we add any such midpoints to \(\widehat{M}^{(i+1)}_{k-1}\) where necessary, and the next stage is then completed by setting \(M^{(i+1)}_k:= f(\widehat{M}^{(i+1)}_{k-1})\). Importantly, all points in \(M^{(i+1)}_k\) are again f-images of points obtained from cubic spline interpolation on \(M_{k-1}\), which can be done within a controlled bound on the interpolation error. The use of cubic spline interpolation, rather than linear interpolation as done in [31, 33], allows us to achieve a prespecified overall accuracy with considerably fewer mesh points, thus making it practical to compute extremely long pieces of manifolds many times.

Mesh points may be mapped too close to each other according to condition (7). In fact, this often occurs when working in a compactified coordinate space, such as the cylinder \(\mathcal {C}\) for the map \(\mathcal {H}\). Notice in Figs. 1 and 2 that the stable manifold comes ever closer to the source \(s^+ = (-1,0,1)\) on the bounding circle at infinity; hence, there is strong local contraction and large curvature of the manifold near \(s^+\), leading to a serious concentration of mesh points. A related practical issue is that the evaluation of the transformation (9) to and from compactified coordinates becomes numerically problematic for points near the boundary of \(\mathcal {C}\) (this may give NaN in Matlab), as well as for points with one of the coordinates very close to 0 (this may lead to unintended sign changes at infinity). This is why we first remove such points from the meshes \(M_{k-1}\) and \(M^{(0)}_{k} = f(M_{k-1})\), before starting the mesh refinement procedure. We then add points in stages as explained above, until conditions (4)–(6) are satisfied. Only then do we check condition (7) and remove mesh points. Since this again uses vector operations in Matlab, performing point removal only at the end considerably enhances the efficiency of our method, compared to doing so at every stage. Specifically, we simultaneously identify all points that violate (7) and then remove mesh points as follows: a local \(\Delta _{\min }\)-equidistance ‘comparing mesh’ is created by spline interpolation, and we then only keep the points that are closest to these ‘ideal’ mesh points.

The algorithm as implemented in Matlab is available with this paper; see Appendix B. While the iterative process of constructing the mesh \(M_k\) for each step k with the stages \(M_k^{(i)}\) may seem cumbersome, it can be implemented efficiently in the Matlab programming environment. To represent the mesh as it is being computed, we use a dynamic data structure in the form of an ordered list of objects/structures, which contain not only the coordinates of the mesh points but also additional information and labels. Once a mesh M has been computed up to a desired number of steps (iterations of the fundamental domain), the intersection points with a section of interest are found in a post-processing step—by detecting sign changes with respect to a normal vector and then finding the corresponding intersection point with cubic spline interpolation from nearby points. Each intersection point is then inserted at the correct position into the ordered list representing the manifold, with a flag identifying it as such. From the respective data structures computed for a suitable set of \(\xi \)-values, the curves \(\{\varvec{w}(l)\}\) can now readily be identified and plotted.

Fig. 4
figure 4

Illustration of the stages of mesh refinement in the computation of \(W^s(p^-)\) with \(f = \mathcal {H}^{-2}\). a shows the manifold \(W^s(p^-)\) (light blue) and the segments \(V_0\) and \(V_1=\mathcal {H}^{-2}(V_0)\) (dark blue). b1 and b2 show the mesh \(M_0\) with 34 points \(q_j\) (blue dots) and \(M^{(6)}_0\) with additional 57 interpolated points \(\hat{q}_j\) (red dots) of \(V_0\). Similarly, c1–c3 show, respectively, the meshes \(M_1^{(0)} = \mathcal {H}^{-2}(M_0)\), \(M_1^{(1)} = \mathcal {H}^{-2}(\widehat{M}^{(1)}_0)\) and \(M_1= M^{(6)}_1 =\mathcal {H}^{-2}(\widehat{M}^{(6)}_0)\) representing \(V_1\)

Full size image

2.2 Illustration of mesh refinement

Figure 4 illustrates a step of the algorithm for the one-dimensional stable manifold \(W^s(p^-)\) of \(\mathcal {H}\) for \(\xi = 1.2\) from Figs. 1 and 2a. For illustrative purposes, we show how an initial mesh \(M_{0}\), now representing a chosen segment \(V_{0}\) of \(W^s(p^-)\) (rather than an entire, long fundamental domain), is mapped and then refined in stages to obtain the mesh \(M_1\) of the image \(V_1\). Here, we consider the second iterate \(\mathcal {H}^{-2}\) and use the quite low accuracy setting with \(\Delta _{\max }=0.05\), \(\alpha _{\max }=0.3\), \(\Delta \alpha _{\max }=0.03\) and \(\Delta _{\min }=0.01\), which ensures that the mesh points are not too close together and, hence, visible. How the accuracy parameters should be chosen to ensure accurate results is discussed in Sect. 3.

Figure 4a shows the segments \(V_0\) and \(V_1=\mathcal {H}^{-2}(V_0)\) on \(W^s(p^-)\) in the cylinder \(\mathcal {C}\). Notice that the image \(V_1\) is much longer than \(V_{0}\) and features a close pass near the source \(s^+\). Panel (b1) shows the mesh \(M_0\) on \(V_0\) with 33 points \(q_j\) that satisfy the specified accuracy conditions, and panel (c1) its image \(M_1^{(0)} = \mathcal {H}^{-2}(M_0)\). Notice the concentration of the 34 mesh points of \(M_1^{(0)}\) near \(s^+\), leading to a very poor representation of \(V_1\) overall. At stage 1, mesh points are added in between too distant mesh points, as is illustrated in panel (c2), and the process repeats in subsequent stages until the accuracy conditions are satisfied at all mesh points. This occurs after 6 stages and results in the mesh \(M_1 = M_1^{(6)}\) in panel (c3), which is the \(\mathcal {H}^{-2}\)-image of the refined mesh \(\widehat{M}^{(6)}_0\) shown in panel (b2). Notice that the required 57 interpolated points \(\hat{q}_j \in \widehat{M}^{(6)}_0\) are very much concentrated near the two ends of \(V_0\), which illustrates the highly nonlinear nature of \(\mathcal {H}^{-2}\). Moreover, the distribution of points along \(M_1\) is, indeed, determined according to curvature as a result of conditions (5) and (6). No mesh points need to be removed in this illustrating example.

3 Choice of the accuracy parameters

Our goal is to compute the stable manifold and the associated \(\xi \)-parametrised intersection sets \(\{\varvec{w}(l)\}\) not only for the second fixed point \(p^-\), as in Fig. 2, but also the other invariant objects shown in Fig. 1, namely:

  • The second fixed point \(p^+\), which has negative eigenvalues;

  • The period-two orbit \(p_2=\{p^1_2,p^2_2\}\), which also has negative eigenvalues; and

  • The period-three orbit denoted \(p_3=\{p^1_3,p^2_3,p^3_3\}\), which has positive eigenvalues.

As Sect. 4 will show, this allows us to obtain a more comprehensive representation of the stable manifold \(W^s(\Lambda )\) of the invariant set \(\Lambda \) itself, thereby enhancing our understanding of when it is a blender and when not.

Table 1 Required iterates and associated local expansion rates \(1/\lambda ^s\) at the fixed points of the appropriate iterate for the computation of the manifolds \(W^s(p^-)\), \(W^s(p^+)\), \(W^s(p_2)\) and \(W^s(p_3)\); here, \(\xi = 1.2\) (with \(\alpha =4.2\) and \(\beta =0.3\))
Full size table

Our algorithm computes these manifolds as branches of orientation-preserving fixed points, obtained by considering an appropriate iterate of \(\mathcal {H}^{-1}\). Specifically, we compute the relevant branches of \(W^s(p^+)\) with \(\mathcal {H}^{-2}\), of \(W^s(p^1_2)\) and \(W^s(p^2_2)\) with \(\mathcal {H}^{-4}\), and of \(W^s(p^1_3)\), \(W^s(p^2_3)\), and \(W^s(p^3_3)\) with \(\mathcal {H}^{-3}\). Unfortunately, a higher iterate means stronger expansion along the respective branch, which makes the computations more challenging. Table 1 shows, for \(\xi =1.2\), the local expansion rates for each of these fixed points, defined as \(1/\lambda ^{s}\), where \(\lambda ^{s}\) is the stable eigenvalue of the corresponding higher iterate of \(\mathcal {H}\).

General theory states that an initial piece of an invariant manifold can be approximated as accurately (in the Hausdorff metric) as desired, provided the relevant accuracy parameters—which are \(\Delta _{\max }\), \(\alpha _{\max }\), \((\Delta \alpha )_{\max }\), and \(\Delta _{\min }\)— are chosen small enough; see [33] and also [34, 35]. However, this type of proof requires knowledge of Lipschitz constants throughout phase space, which is not available in practice. Hence, there are no a priori estimates of what the accuracy settings should be. Rather, they need to be chosen by showing that the computed (finite) piece of manifold remains the same (to within a certain precision) under further refinement of the accuracy.

As the result of extensive numerical investigations, we found that the parameter values

$$\begin{aligned} \Delta _{\max }=0.01, \quad \Delta _{\min }=10^{-6}, \quad \Delta \alpha _{\max }=0.001 \quad \text {and} \quad \alpha _{\max }=0.3, \end{aligned}$$
(8)

enable the accurate computation of all the above branches of manifolds up to the desired total number of \(2^{11}\) intersection points with the section \(\Sigma \) and over the \(\xi \)-range \(\xi =[1.02,3]\) of interest; we refer to (8) as the reference accuracy settings. Indeed, the figures in Sect. 4 were obtained with these accuracy settings and, hence, provide comprehensive numerical evidence that this choice is suitable for our purposes.

To give an impression of how such a choice can be made, we now discuss and illustrate in Figs. 5, 6 and 7 the influence of the different accuracy parameters. By way of an illustrative case study, we again consider the relevant branch of \(W^s(p^-)\) from Fig. 4, but now computed with the higher iterate \(\mathcal {H}^{-4}\). Note that the local expansion rate near \(p^-\), as a fixed point of \(\mathcal {H}^{-4}\), has the artificially high value \((18.486)^{4} \approx 1.1679 \times 10^{5}\); compare with Table 1.

Fig. 5
figure 5

Illustration of the influence of \(\Delta _{\max }\) on the accuracy of the computation of \(W^s(p^-)\) with \(\mathcal {H}^{-4}\), with the other accuracy parameters as in (8). a–c show four computations, for \(\Delta _{\max }\) as given in the legend, represented by plots of \(\bar{x}\), \(\bar{y}\), and \(\bar{z}\), which are synchronised to the intersections with \(\Sigma \) from 2040 to 2048. d1 and the enlargement d2 show that, past the 22nd intersection with \(\Sigma \), the less accurate computations (orange and green) already start to deviate in \(\bar{y}\) from the reference branch for \(\Delta _{\max } = 0.01\) (blue). The enlargement in a shows the distribution of points for \(\Delta _{\max } = 0.01\) (blue) and for \(\Delta _{\max } = 0.005\) (gold) near the last shown maximum of \(\bar{x}\)

Full size image

We first consider the effect on the accuracy of the parameter \(\Delta _{\max }\), which controls the overall density of mesh points. Figure 5 shows \(W^s(p^-)\) computed with \(\mathcal {H}^{-4}\) for the four different values of \(\Delta _{\max }\) given in the legend; all other accuracy parameters are as in (8). To allow for their comparison, the computed branches are synchronised to their respective ordered intersections \(\Sigma \) (where \(\bar{x} = \bar{y}\)). Panels (a), (b), and (c) show them in between intersection numbers 2040 to \(2048 = 2^{11}\) (the end of the computations) in terms of \(\bar{x}\), \(\bar{y}\), and \(\bar{z}\) along the respective branch. The reference branch, for \(\Delta _{\max }=0.01\) as in (8), is sufficiently accurate, which is demonstrated by the fact that the computation for the reduced value of \(\Delta _{\max }=0.005\) agrees with it to within \(O(10^{-5})\) at all (ordered) intersection points with \(\Sigma \). This is illustrated further with the enlargement near the final maximum of \(\bar{x}\): while the branch for \(\Delta _{\max }=0.005\) has a finer mesh resolution, as expected, it nevertheless lies in a small neighbourhood of the reference branch, and vice versa.

For the larger values of \(\Delta _{\max }=0.05\) and \(\Delta _{\max }=0.1\), on the other hand, the computed branches do not agree at all with the reference branch near the shown end of these computations. Figure 5d shows that the difference between these branches becomes noticeable just after the 22nd intersection with \(\Sigma \). As is to be expected, lower accuracy settings lead to a gradual deviation from the reference solution (and, hence \(W^s(p^-)\) itself), as is illustrated in panel (d1) with the plot of \(\bar{y}\); the enlargement panel (d2) shows that this process starts just past the first minimum. Note that the computed branches for \(\Delta _{\max }=0.05\) and \(\Delta _{\max }=0.1\) still appear to be very close together here, while this is clearly not the case further along \(W^s(p^-)\), as illustrated in panels (a)–(c). Overall, Fig. 5 shows that the reference accuracy settings are sufficient to compute \(W^s(p^-)\), even with the fourth iterate \(\mathcal {H}^{-4}\), up to \(2^{11}\) intersection points with \(\Sigma \).

Fig. 6
figure 6

Illustration of the influence of \(\Delta \alpha _{\max }\) on the computation of \(W^s(p^-)\) with \(\mathcal {H}^{-4}\), presented as in Fig. 5 for accuracy settings given in the legend with \(\Delta _{\min }\) and \(\alpha _{\max }\) as in (8)

Full size image
Fig. 7
figure 7

Illustration of the influence of \(\Delta _{\min }\) on the computation of \(W^s(p^-)\) with \(\mathcal {H}^{-4}\), presented as in Fig. 5 for \(\Delta _{\min }\) as given in the legend with the other accuracy parameters as in (8)

Full size image

Figure 6 illustrates the role of \((\Delta \alpha )_{\max }\), which controls the effect of the curvature of the manifold on the distribution of the mesh points. We again show plots of \(\bar{x}\), \(\bar{y}\), and \(\bar{z}\) at the end of the different computed branches in between intersection numbers 2040 and \(2048 = 2^{11}\). Specifically, we show the reference branch with (8) and compare it with computations for the larger value of \(\Delta _{\max }=0.05\), while varying \((\Delta \alpha )_{\max }\). This illustrates that an appropriate reduction of \((\Delta \alpha )_{\max }\) allows one to achieve an accuracy that is comparable to that of the reference branch. As we know from Fig. 5, the branch computed with \((\Delta \alpha )_{\max } = 0.001\) is not accurate, and Fig. 6 shows that this is still the case for \((\Delta \alpha )_{\max } = 0.0005\) and for \((\Delta \alpha )_{\max } = 0.0003\). For the value of \((\Delta \alpha )_{\max } = 0.0001\), on the other hand, the computed branch is again effectively indistinguishable from the reference branch. This is further illustrated with the enlargement near the last maximum of \(\bar{x}\) in panel (a). Notice that the branch computed with \(\Delta _{\max }=0.05\) and \((\Delta \alpha )_{\max } = 0.0001\) has considerably more mesh points than the reference branch. Hence, Fig. 6 also illustrates that it is computationally more efficient to strike a good balance between the values of \(\Delta _{\max }\) and \((\Delta \alpha )_{\max }\), as is the case for the reference accuracy settings (8).

In situations involving extremely sharp corners, particularly when the manifold undertakes excursions near infinity as it comes close to the source \(s^+\) on the boundary circle of the cylinder \(\mathcal {C}\), the curvature conditions (5) and (6) might result in a substantial number of mesh points being generated with extremely small distances between them. Condition (7), with the accuracy parameter \(\Delta _{\min }\), serves to avoid this, while still ensuring that sharp turns are resolved accurately enough. As a general guideline, it is important to satisfy \(\Delta _{\min }\alpha _{\max } < (\Delta \alpha )_{\max }\), implying that we must have \(\Delta _{\min } < 0.0033\) in the reference accuracy settings.

Figure 7 illustrates the impact of \(\Delta _{\min }\) on the accuracy of the computation by showing the reference branch with \(\Delta _{\min } = 10^{-6}\) in comparison with computations for larger values of \(\Delta _{\min }\). The plots of \(\bar{x}\), \(\bar{y}\), and \(\bar{z}\) in between intersection numbers 2040 and \(2048 = 2^{11}\) in panels (a)–(c) show that the branch is not computed correctly for the two largest chosen values of \(\Delta _{\min }\). For a computation with \(\Delta _{\min } = 10^{-5}\), on the other hand, there is no distinguishable difference with the reference branch. The enlarged minimum of \(\bar{x}\) in panel (a) corresponds to a close visit to \(s^+\) (while the enlarged maximum does not), and the reference branch has considerably more mesh points than the branch for \(\Delta _{\min } = 10^{-5}\); notice that the curvature condition is determined not just from \(\bar{x}\), but from the corresponding points in \(\mathcal {C}\) with components of \(\bar{y}\) and \(\bar{z}\) as well. Hence, the reference branch resolves excursions to infinity considerably better, and still with a reasonably low number of mesh points, which is why we chose \(\Delta _{\min } = 10^{-6}\) as the reference value.

4 Stable manifolds and curves of intersection points

The one-dimensional stable manifolds of different fixed and periodic points cover the stable manifold \(W^s(\Lambda )\) of a hyperbolic set \(\Lambda \) in different ways. Hence, computing more than one such curve is an efficient way to check for and illustrate the carpet property, compared to computing excessively long branches of a single manifold (which would also require even more stringent accuracy settings). In this section, we demonstrate this with the minimum required iterate of \(\mathcal {H}^{-1}\) and the reference accuracy settings (8).

Fig. 8
figure 8

The stable manifolds \(W^s(p^-)\) (blue curve) in a and \(W^s(p^+)\) (purple curve) in b in \(\mathcal {C}\) with the fixed points \(p^-\), \(p^+\), and \(s^+\), shown with a colour gradient according to arclength \(\ell \). Here, the relevant branch of \(W^{s}(p^-)\) is computed up to \(2^{11}\) and the two branches of \(W^{s}(p^+)\) up to \(2^{10}\) intersections with \(\Sigma \)

Full size image

Figure 8 directly compares the manifolds \(W^s(p^-)\) (blue) and \(W^s(p^+)\) (purple) for \(\xi =1.2\), when \(\Lambda \) is a blender. The relevant branch of \(W^s(p^-)\) (that does not immediately go to infinity) has been computed with \(\mathcal {H}^{-1}\) up to \(2^{11}\) intersections with \(\Sigma \), as in previous figures. The two branches of \(W^s(p^+)\) of the orientation-reversing fixed point \(p^+\), on the other hand, have been computed with \(\mathcal {H}^{-2}\) up to \(2^{10}\) intersections with \(\Sigma \). Hence, both manifolds have comparable total arclength of about 12,000. To illustrate how these branches visit different parts of \(\mathcal {C}\), their arclength is represented by a colour gradient from dark to light. Comparison between panels (a) and (b) of Fig. 8 illustrates how (the computed pieces of) \(W^s(p^-)\) and \(W^s(p^+)\) visit and slowly ‘fill’ different parts of the cylinder \(\mathcal {C}\) as a function of arclength from the respective fixed point; in the limit of infinite arclength, both manifolds cover \(W^s(\Lambda )\). Even though neither of these two computed manifolds is very close to this limit, together, they nicely complement each other and give a better overall impression of \(W^s(\Lambda )\).

Fig. 9
figure 9

The stable manifolds \(W^s(p^-)\) (blue) and \(W^s(p^+)\) (purple) and their intersection points with the plane \(\Sigma \). a and b show them in \(\mathcal {C}\) for \(\xi =1.2\) and \(\xi =2.4\), respectively, and c shows the corresponding sets of curves \(\{\varvec{w}^-(l)\}\) (blue) and \(\{\varvec{w}^+(l)\}\) (purple) of ordered intersection points in projection onto \(\bar{z}\) for \(\xi \in [1.02,3]\). Here, the relevant branch of \(W^s(p^-)\) is computed up to \(2^{11}\) and each branch of \(W^s(p^-)\) up to \(2^{10}\) intersection points

Full size image

In Fig. 9, we consider both \(W^s(p^-)\) and \(W^s(p^+)\) together to help identify for which values of \(\xi \) the hyperbolic set \(\Lambda \) has the carpet property. Here, the computed branches are plotted (without shading according to arclength) together with the section \(\Sigma \) and the respective intersection sets \(\{w^-(l)\}\) and \(\{w^+(l)\}\) associated with \(p^-\) and \(p^+\), respectively. Panel (a) shows the computed branches for \(\xi = 1.2\) from Fig. 8, nicely illustrating the carpet property and how \(W^s(p^-)\) and \(W^s(p^+)\) are ‘mixing’ on \(W^s(\Lambda )\). Figure 9b shows these two manifolds for \(\xi = 2.4\) in the same way, where they do not appear to have the carpet property but look locally like a Cantor set of curves from any direction. Note also that it is practically impossible to distinguish \(W^{s}(p^-)\) and \(W^s(p^+)\) for this case, because they are not covering a ‘large’ set in \((\bar{x},\bar{y},\bar{z})\)-space.

Fig. 10
figure 10

The stable manifold \(W^s(p_2)\) of the orientation-reversing period-two orbit \(p_2=\{p^1_2,p^2_2\}\), consisting of \(W^s(p^1_2)\) (dark green) and \(W^s(p^2_2)\) (light green), and their intersection points with the plane \(\Sigma \). a and b show them in \(\mathcal {C}\) for \(\xi =1.2\) and \(\xi =2.4\), respectively, and c shows the corresponding sets of curves \(\{\varvec{w}^1_2(l)\}\) (dark green) and \(\{\varvec{w}^2_2(l)\}\) (light green) in projection onto \(\bar{z}\) for \(\xi \in [1.02,3]\). Here, the respective four branches are computed up to \(2^{10}\) intersection points each

Full size image
Fig. 11
figure 11

The stable manifold \(W^s(p_3)\) of the orientation-preserving period-three orbit \(p_3=\{p^1_3,p^2_3, p^3_3\}\), consisting of \(W^s(p^1_3)\) (dark cyan), \(W^s(p^2_3)\) (medium cyan), and \(W^s(p^3_3)\) (light cyan), and their intersection sets with the plane \(\Sigma \). a and b show them in \((\bar{x},\bar{y},\bar{z})\)-space for \(\xi =1.2\) and \(\xi =2.4\), respectively, and c shows the corresponding sets of curves \(\{\varvec{w}^1_3(l)\}\) (dark cyan), \(\{\varvec{w}^2_3(l)\}\) (medium cyan), and \(\{\varvec{w}^3_3(l)\}\) (light cyan) in projection onto \(\bar{z}\) for \(\xi \in [1.02,3]\). Here, the relevant three branches (that do not immediately go to infinity) are computed up to \(2^{10}\) intersection points each

Full size image
Fig. 12
figure 12

Illustration of the carpet property for \(\xi =1.2\) by simultaneously plotting \(W^{s}(p^\pm )\), \(W^{s}(p_2)\), and \(W^s(p_3)\) from Figs. 9–11a; the hyperbolic set \(\Lambda \) is represented by the intersection points of these manifolds with \(W^u(p^-)\) from Fig. 1

Full size image

Figure 9c shows the projection onto \(\bar{z}\) of the two sets of \(2^{11}\) curves \(\{\varvec{w}^-(l)\}\) and \(\{\varvec{w}^+(l)\}\) over the \(\xi \)-range [1.02, 3]. Also shown are the curves traced out by \(p^-\) and \(p^+\), and vertical lines indicating the values of \(\xi = 1.2\) and \(\xi = 2.4\) from panels (a) and (b). In comparison with Fig. 2c, the plot with both sets of curves \(\{\varvec{w}^-(l)\}\) and \(\{\varvec{w}^+(l)\}\) gives a better impression of how the intersection set of \(W^s(\Lambda )\) with \(\Sigma \) changes with the expansion rate \(\xi \) of \(\mathcal {H}\); this is especially the case for weaker expansion, below \(\xi \approx 1.5\). Indeed, Fig. 9c illustrates how all visible gaps between these \(2^{11}\) curves of projected intersection points with \(\Sigma \) appear to close just to the right of \(\xi = 1.6\) as \(\xi \) is decreased from 3. For \(\xi \approx 1\), on the other hand, the expected \(\bar{z}\)-range is no longer covered by the computed sets of curves \(\{\varvec{w}^-(l)\}\) and \(\{\varvec{w}^+(l)\}\), which all lie near \(p^-\) and \(p^+\), respectively; because of the very weak expansion \(\xi >1\), the first \(2^{11}\), respectively \(2^{10}\) pairs of intersection points are not sufficient to represent \(W^s(\Lambda )\) here.

Figure 10 shows in the same way the stable manifold \(W^s(p_2)\) and its intersection sets with the section \(\Sigma \) of the orientation-reversing period-two orbit \(p_2=\{p^1_2, p^2_2\}\); here, each of the two branches of the manifolds \(W^s(p^1_2)\) and \(W^s(p^2_2)\) of the fixed points \(p^1_2\) and \( p^2_2\) under \(\mathcal {H}^{-4}\) is computed and shown up to \(2^{10}\) intersections with \(\Sigma \). Panel (a) now illustrates with \(W^s(p_2)\) that \(\Lambda \) has the carpet property for \(\xi = 1.2\). Notice that the individual manifolds \(W^s(p^1_2)\) and \(W^s(p^2_2)\) are also thoroughly mixing as they cover \(W^s(\Lambda )\). This is also the case when \(\Lambda \) is not a blender, as shown in panel (b) for \(\xi = 2.4\). The associated sets of curves \(\{\varvec{w}^1_2(l)\}\) and \(\{\varvec{w}^2_2(l)\}\) of ordered intersection points are shown for \(\xi \in [1.02,3]\) in panel (c). At this scale, these curves appear to cover \(W^s(\Lambda )\) just as well as the corresponding sets of curves \(\{\varvec{w}^\pm (l)\}\) for \(W^s(p^\pm )\) when \(\xi \) lies above about 1.2, but they also remain noticeably further from \(p^-\) and \(p^+\) for \(1< \xi < 1.2\); compare with Fig. 9c. The discrepancy is again due to the weak expansion near \(p_2\) for small \(\xi \).

Figure 11 shows the data obtained by computing the stable manifold \(W^s(p_3)\) of the period-three orbit \(p_3=\{p^1_3,p^2_3,p^3_3\}\), which is orientation preserving. Each of the respective branches of the three manifolds \(W^s(p^1_3)\), \(W^s(p^2_3)\), and \(W^s(p^3_3)\) that do not immediately go to infinity was computed up to \(2^{10}\) intersections with \(\Sigma \) as the corresponding manifold of a fixed point under \(\mathcal {H}^{-3}\). Again, these three curves are mixing well as they cover \(W^s(\Lambda )\), whether \(\Lambda \) is a blender, as in panel (a) for \(\xi =1.2\), or not, as in panel (b) for \(\xi =2.4\). Notice further that the three sets of curves \(\{\varvec{w}^1_3(l)\}\), \(\{\varvec{w}^2_3(l)\}\), and \(\{\varvec{w}^3_3(l)\}\) of ordered intersection points come considerably closer to \(p^-\) and \(p^+\) even for quite low values of \(\xi \); compare with those computed for \(W^s(p_2)\), as shown in Fig. 10c.

5 Conclusions

We presented an algorithm to compute extremely long pieces of a one-dimensional (un)stable manifold of a map as an arclength-parametrised curve. While this can be achieved, in principle, with readily available methods, such as the growth method from [33] on which our work is based, our method is much improved: a large number of ordered intersection points on the manifold with a given section can be identified efficiently and very accurately. As we showed, this enables us to find and render all of these intersection points as well-resolved curves as a system parameter is varied. Efficiency and accuracy are indeed key here, as this required manifold computations of extremely long arclength for parameter values on a sufficiently fine grid over the parameter range of interest. To achieve this, we combined fundamental domain iteration with ideas from [33] and [31] to obtain a curvature-dependent and ordered mesh representation. Our method is implemented in Matlab and more efficient because the required mesh refinement of a new fundamental domain is performed in an iterative fashion for entire vectors of mesh points.

We demonstrated the performance of our method with the example of the family of three-dimensional Hénon-like maps, which is known to feature a blender for a certain range of an important parameter, namely, the weakest linear expansion rate. We computed a large number of parametrised curves of intersection points with a plane, not only for the stable manifolds of fixed points but also for those of period-two and period-three orbits. This constitutes a considerably enhanced capability compared to earlier work in [28,29,30], where only the manifolds of the fixed points were considered, and their intersection points with a plane were not consistently ordered for different values of the expansion rate. In particular, we are able to illustrate the characterising carpet property of a blender in an unprecedented way by plotting all of the computed manifolds together in a single figure. This is done in Fig. 12, where the hyperbolic set \(\Lambda \) is approximated and represented very accurately by the 90,112 intersection points of the different stable manifolds with the first part of the unstable manifold from Fig. 1. Indeed, Fig. 12 stands as the most comprehensive illustration of a blender and its one-dimensional manifold.

We presented an advanced numerical tool for the study of complicated types of dynamics of dynamical systems defined by maps, with considerable potential to contribute to the understanding of what is known as ‘wild chaos’. In current related work, we employ this approach to identify, in a parameter-dependent way, which specific intersection points of which stable one-dimensional manifold bound certain gaps in projection. In this way, we are able to clarify in a topological fashion the process of how exactly these gaps close when the expansion rate is decreased and the blender emerges; these results will be presented elsewhere. Computing extremely long pieces of manifolds also holds considerable promise for the study of heterodimensional cycles. The construction of explicit families of maps and the subsequent study of their manifold structure near heterodimensional cycles is a challenging subject of ongoing and future research.