Matching geometric and expansion characteristics of wild chaotic attractors

Wild chaotic attractors exhibit chaotic dynamics with a robustness property that cannot be destroyed with small perturbations. We consider a discrete-time system with the smallest possible dimension, namely, defined by a non-invertible map on the complex plane. For this map, wild chaos has been proven to exist in a small parameter region. Recently, it was conjectured to exist in a much larger region of parameter space, past a so-called backward critical tangency, at which a sequence of pre-images of a critical point converges to a saddle fixed point. Geometrically, a backward critical tangency leads to an abundance of homoclinic and heteroclinic tangencies between invariant manifolds of different dimensions, generating precisely what are believed to be the necessary ingredients for wild chaos. In this paper, we present corroborating evidence for this conjecture by computing Lyapunov exponents associated with the attractor. When the sum of the two (largest) Lyapunov exponents is positive, the dynamics is wild chaotic for this non-invertible map. We find that the zero-sum locus matches the locus of backward critical tangency, confirming its role as a boundary of existence of wild chaos.


Introduction
The definition of chaos that is commonly used in dynamical systems theory is that a system is chaotic if there exists a compact invariant set that satisfies the following three properties [5]: (i) trajectories in the invariant set exhibit sensitive dependence on initial conditions; (ii) the invariant set contains a non-periodic trajectory (the system is topologically transitive); (iii) periodic orbits are dense in the invariant set.
The invariant set need not be an attractor, but for the purpose of this paper, we are only interested in the case of a chaotic attractor. Attracting chaotic dynamics can also be defined in terms of the Lyapunov exponents associated with the (attracting) invariant set [25,32]. Lyapunov exponents are the expansion or contraction rates of the linearised system along an arbitrary nonperiodic trajectory on the attractor. In many cases, the attractor will be chaotic in the sense as defined above, if it has a positive Lyapunov exponent. Computing Lyapunov exponents is a practical way of demonstrating the presence of a chaotic attractor in experiments, or from numerical simulation. In dynamical systems theory, chaotic attractors are studied with a combination of geometrical and numerical tools. It is well a e-mail: h.m.osinga@auckland.ac.nz (corresponding author) known that typical trajectories on a chaotic attractor exhibit an expanding direction, which is measured by a positive Lyapunov exponent [25,32]. Geometrically, one is interested in the so-called backbone or skeleton of the dynamical system, which comprises fixed points, periodic orbits, and other compact invariant manifolds, along with their stable and unstable manifolds defined by families of trajectories that converge to these compact invariant manifolds in forward or backward time, respectively. A chaotic attractor contains infinitely many periodic orbits that are all of saddle type; their stable and unstable manifolds intersect in homoclinic and heteroclinic trajectories. The precise topological nature of the attractor is controlled via tangencies between such manifolds, which give rise to an accumulation of tangencies in one-parameter intervals [31,36]. The number of Lyapunov exponents for a dynamical system is the same as the dimension of its phase space. Hence, if chaos exists already when a single Lyapunov exponent is positive, what can one say about the other Lyapunov exponents? Or perhaps phrased differently, are there different kinds of chaos that can be classified by properties of the Lyapunov exponents?
The idea of a chaotic attractor with two or more positive Lyapunov exponents goes back to Rössler [33] who calls this type of chaotic behaviour hyperchaos. In this paper, we focus on the weaker property that only the sum, and not both, of the largest two Lyapunov exponents is positive. This property is associated with a different form of higher-dimensional chaos that is called wild chaos in the more theoretical litera- ture [26]; wild chaos is defined by the robust existence of homoclinic and heteroclinic tangencies, that is, stable and unstable manifolds of a hyperbolic set remain tangent to each other over an entire parameter interval. Hence, wild chaos persists under C 1 -perturbation (in particular, there are no so-called periodic windows) [3,37]. As yet, it is an open question which geometric ingredients, or in other words, which structure of the backbone of the dynamical system are necessary and sufficient for it to exhibit wild chaos. In previous work by one of the authors, the geometric ingredients necessary for generating wild chaos were investigated for a smooth non-invertible map on the complex plane [17], which we call the Lorenz-like double-cover map. This example was chosen because the Lorenz-like doublecover map has the lowest dimension possible and wild chaos was proven to exist in a small neighbourhood of a specific point in parameter space [2]. The theory of wild chaos tends to consider only saddle-type hyperbolic sets, rather than attractors with this property. More recently, wild chaotic attractors have been a focus of study, which are called pseudo-hyperbolic attractors in [14,15].
In this paper, we focus specifically on the parameter region for which the Lorenz-like double-cover map has an attractor and we compute the two Lyapunov exponents along representative trajectories over a range of two parameters. From [15], we know that wild attractors not only have a positive Lyapunov exponent, but the sum of the two largest Lyapunov exponents is also positive; in general, this is not a sufficient condition for wild chaos, but due to the non-invertibility of the Lorenz-like double-cover map, a positive Lyapunov sum is both necessary and sufficient for our study.
The Lorenz-like double-cover map from [2] is defined on the complex plane as the non-invertible map with real parameters a > 0 and λ ∈ (0, 1). In [2], wild chaos was proven to exist for system (1) provided (a, λ) ≈ (1, 1). However, wild chaos is conjectured to exist for a much larger region in the (a, λ)-plane [17,29]. Note that f is not defined for z = 0. This point is called the critical locus. The conjecture made in [17] is that wild chaos arises from a specific global bifurcation, called the (first) backward critical tangency, defined as the moment when the critical locus is contained in the unstable set of a saddle periodic orbit, namely, a fixed point of f . We wish to test this hypothesis more globally by comparison between the zero-sum locus of the two Lyapunov exponents associated with the chaotic attractors of system (1) and the conjectured locus of onset of wild chaos given by the curve of backward critical tangency in the (a, λ)-plane.
In the next section, we provide basic definitions relevant for a non-invertible system like system (1). In Sect. 3, we discuss a well-known route to chaos, which is representative of the parameter-dependent transforma-tions of the dynamics for system (1). We explain how we compute Lyapunov exponents in Sect. 4 and compare the sign of the maximal Lyapunov exponent in the (a, λ)-plane with the existence of a classical chaotic attractor. In Sect. 5, we consider the sum of the two Lyapunov exponents and overlay the curve of first backward critical tangency in the (a, λ)-plane. Our findings suggest that the zero-sum locus of the Lyapunov exponents agrees well with the locus of first backward critical tangency. We end with conclusions in Sect. 6 and discuss future directions of research.

Background and definitions
Let us first explain some properties of dynamical systems described by a non-invertible map that are different from the more standard invertible case. Recall that system (1) is not well defined on all of C, because f is not defined for z = 0. The image of f covers the complex plane twice, except for a (closed) disk with radius 1−λ centred at 1, which is not in its range. Consequently, points inside this disk have no pre-images, while points outside it have two; hence, trajectories are not, or not uniquely defined in backward time. We remark that f (z) and f (z) are complex conjugate, so the system has the symmetry z →z, and the real axis is invariant. Following the theory for non-invertible maps [22,27], we refer to 0 as the critical point, also denoted J 0 , and the boundary of the closed disk is the critical set, denoted J 1 , that separates the points in the complex plane that have two pre-images of f from those that have none; note that J 1 can be thought of as the muli-valued image of the critical point 0. The sequences of pre-images (backward iterates) of 0 and images (forward iterates) of J 1 form the backward critical set and the forward critical set, respectively.
Observe that f is continuous and its derivative is also continuous. Hence, locally and away from z = 0, we can define a unique inverse for f . This means that locally we can apply the theory for real invertible systems. To this end, we identify C with R 2 and write f in terms of coordinates (x, y) ∈ R 2 with z = x + i y. For a sufficiently small, the map f has a saddle fixed point p on the positive real axis; hence, f (p) = p and the Jacobian matrix of f evaluated at p has one eigenvalue with modulus less than 1 and the other with modulus larger than 1. Since, locally near p, say, in a neighbourhood U p ⊂ R 2 of p, we can identify a unique local inverse of f that maps p to itself, the Stable Manifold Theorem [30] guarantees the existence of local stable and unstable manifolds of p, denoted W s loc (p) and W u loc (p), respectively. They are defined as

Fig. 1
Backbone of system (1) for a = 0.6 and λ = 0.8. Shown are the backward critical set (green dots), containing 0 and its pre-images, and the forward critical set (green closed curves), containing the circle J1 and its forward images, together with the fixed point p (bright-green cross), its stable set W s (p) (blue curves) and unstable set W u (p) (red curve), and the two attracting fixed points q andq (light-blue triangles) and are tangent at p to the corresponding stable and unstable eigenvectors of p; here, the local inverse is assumed for f (−n) in the definition of W u loc (p). These local manifolds can be extended to global stable and unstable sets, denoted W s (p) and W u (p), respectively. The stable set W s (p) contains all points in the complex plane that converge to p under forward iteration, which means that W s (p) comprises W s loc (p) and all possible pre-images of W s loc (p) [11,27]. Similarly, W u (p) contains all points in the complex plane for which there exists a sequence of pre-images that converge to p; this means that W u (p) may have self-intersections at points for which both pre-images lie on sequences of pre-images converging to p [12,22]. For this reason, W s (p) and W u (p) are typically not manifolds, which is why we refer to them as stable and unstable sets. Figure 1 shows the backbone of system (1) generated by the map f with a = 0.6 and λ = 0.8. The origin 0 and the other points in the backward critical set are shown as green dots, and J 1 with the other curves in the forward critical set are the green closed curves. The stable set W s (p) (blue curves) of the fixed point p is composed of disjoint manifolds that are preimages of the positive real axis (which contains p and W s loc ). The closure of W s (p) forms a regular tree of degree four with the points in the backward critical set as vertices. Observe that W s (p) includes points inside the disk bounded by J 1 , because such points do have forward iterates and some converge to p; similarly, W s (p) contains points in the interior of all (green) closed curves in the forward critical set. The unstable set W u (p) of p consist of two complex-conjugate branches and is a manifold for this choice of parameters; indeed, each branch is associated with a uniquely defined inverse of f . The unstable set W u (p) cannot intersect J 1 (or images of J 1 ), because such points do not have (a sequence of) pre-images that could converge to p. Rather than intersect, the forward critical set accumulates on W u (p). Each branch of W u (p) spirals in and accumulates onto an attracting fixed point q orq (light-blue triangles); the complex-conjugate pair q andq is given by q = exp (iπ/3) = 1 2 + 1 2 √ 3 i and, hence, does not depend on a or λ. The fixed points q and q are stable (attractors) when a and λ are small enough; they become repelling in a Neȋmark-Sacker bifurcation [23] when a = 1 2λ , that is, when the eigenvalues of the associated Jacobian matrix lie on the unit circle in the complex plane.
Variation of a or λ may induce an interaction between the stable and unstable sets of system (1). A homoclinic tangency occurs when W s (p) and W u (p) have a nontransverse intersection, which forms a trajectory that converges to p in both forward and backward time. As soon as the intersection is transverse, we speak of W s (p) and W u (p) forming a homoclinic tangle [31]. Intersections between stable and unstable sets of two different fixed points or periodic points also occur and this is called a heteroclinic tangency or tangle. (A periodic point is a fixed point of the k-th iterate of f for some k ∈ N with k > 1, and its stable and unstable sets are defined similarly by using the k-th iterate of f .) The presence of a bounded homoclinic orbit forces the stable and unstable sets to accumulate on themselves, leading to the existence of an accumulating sequence of parameter values of further homoclinic tangencies between the same stable and unstable sets [31].
The geometric study of wild chaos in [17] focused on the interactions of W s (p) and W u (p) with each other and with the forward and backward critical sets. Indeed, the proof of existence of wild chaotic dynamics for (a, λ) sufficiently close to (1, 1) assumes that the chaotic invariant set contains p and W u (p) [2]. A classical chaotic attractor in this situation is given by the closure of W u (p). Hence, a bifurcation that changes the chaotic attractor to a wild chaotic attractor is expected to involve W u (p).   Figure 2 shows phase portraits of this transition for λ = 0.8, with panels (a-d) showing phase portraits for a = 0.63, a = 0.65, a = 0.67, and a = 0.7, respectively. As in Fig. 1, we plot the fixed point p (bright-green cross) with its stable and unstable sets W s (p) (blue curves) and W u (p) (red curve), respectively, together with the forward critical set (green closed curves), the backward critical set (green dots), and the complex-conjugate fixed points q andq (darkred squares). Indeed, a > 1 2λ = 0.625 for all four cases, so that q andq are unstable (sources). Furthermore, we also plot the different attractors that exist for the specified parameter values. Figure 2a illustrates that the Neȋmark-Sacker bifurcation initially gives rise to two complex-conjugate attracting invariant circles (light-blue curves) that surround q orq. The two complexconjugate branches of W u (p) (red curves) each accumulate on the respective invariant circles; the dynamics on the two attractors is quasi-periodic and there is no chaos. Here, the basins of attraction of the two invariant circles are separated by the closure of W s (p). From the theory, we know that such quasi-periodic invariant circles are not likely to persist over a large range of parameters [34]. Indeed, as a increases to a = 0.65, we appear to be close to their break-down: surrounding each invariant circle is a (complex-conjugate) pair of period-8 orbits, one of saddle type (green crosses) and the other attracting (light-blue triangles). Figure 2b shows that the invariant circle comes close to the period-8 saddle; there are now a total of four coexisting attractors and each branch of W u (p) appears to accumulate onto both an invariant circle and a period-8 attractor, suggesting the existence of a heteroclinic connection between W u (p) and the stable set of the saddle period-8 orbit (not shown).
When a = 0.67, shown in Fig. 2c, both the invariant circles and the period-8 attractors have disappeared; the sequence of bifurcations that transforms or destroys the non-chaotic attractors occurs in an extremely small parameter interval and we consider the precise details beyond the scope of this paper. For a = 0.67, each branch of W u (p) accumulates on a chaotic attractor, which we represent by 10 4 points of a single orbit (orange), obtained from iterating the (arbitrary) initial condition 1 + i and discarding the first 1000 iterates. Note that the complex-conjugate pair of chaotic attractors lies close, in places, to the stable set W s (p), which forms their shared basin boundary. As a increases further, a homoclinic tangency between W s (p) and W u (p) occurs where the two attractors merge into one large chaotic attractor that is itself complex conjugate; this transition is a form of interior crisis [16]. An example is shown in Fig. 2d for a = 0.7. We remark here that the homoclinic tangle between W s (p) and W u (p) necessarily implies that the unstable set W u (p) has selfintersections and is no longer a manifold; the homoclinic orbits arising from the intersection of W s (p) and W u (p) include points on the real axis that will, by definition, lie on both branches of W u (p). The route to chaos described above for fixed λ = 0.8 is representative for the dynamics exhibited by system (1). (However, the periods of the periodic orbits near the invariant circle will vary.) We remark that the perhaps more common route of period-doubling to chaos [5] also occurs in this system, but only leads to a chaotic attractor if the system is embedded in a larger three-parameter class of non-invertible maps. For example, consider the map f such that a parameter c ∈ C defines the location of the centre of the circle J 1 ; the map (1) is then given by c = 1. In this generalised setting, if a = 2 and λ = 1, the map is equal to the complex quadratic map and has associated period-doubling sequences organised by the Mandelbrot set; see [18] for a more detailed study of the dependence on λ ∈ [0, 1] and c ∈ C for this more generally defined map f with a = 2 fixed.

Lyapunov exponents as a measure to detect wild chaos
Lyapunov exponents represent the average contraction or expansion rates on an attractor that admits an invariant measure [25,28]. An m-dimensional discretetime system will associate m Lyapunov exponents to each trajectory. Hence, an attractor of system (1) has two Lyapunov exponents, Λ 1 ≥ Λ 2 . For a fixed-point attractor, we have Λ i := ln (σ i ) for i = 1, 2, where σ i are the singular values of the Jacobian matrix associated with the fixed point. (The singular values of a matrix A are the square-roots of the eigenvalues of the positive semi-definite matrix product given by A and its transpose.) Similarly, for an attracting periodic orbit of period k, the Lyapunov exponents Λ 1 and Λ 2 are derived from the singular values, say, μ 1 and μ 2 of the Jacobian matrix of the kth iterate. The Jacobian matrix of the kth iterate is, in fact, a chain of k Jacobian matrices for the map evaluated at successive iterates of the initial point. Hence, on average in each of the k iterations, the contraction or expansion rates will be Λ i := 1 k ln (μ i ) for i = 1, 2. Considering now an arbitrary non-periodic trajectory on the attractor, one obtains the Lyapunov exponents from the timevarying linear system given by the linearisation along this trajectory [19]; effectively, we have an infinite chain of Jacobian matrices, and the average contraction or expansion rates must be calculated as a limit where the number of iterates goes to infinity.
The calculation of Lyapunov exponents can numerically be rather delicate, particularly when both expansion and contraction are involved; see [13] for an overview of different approaches. We follow the description in [4] of the method from [10,38], which is essentially a (time-varying) power method with orthonormalisation to avoid alignment of vectors along the direction associated with the largest Lyapunov exponent; see also [7][8][9]. More precisely, assuming z 0 ∈ C is an initial condition on the attractor, we generate the linearised system where z n = f (n) (z 0 ) is the nth iterate of z 0 and Hence, in theory, both the maximal Lyapunov exponent and the sum of the two Lyapunov exponents is obtained as a limit: In practice, we use the values of 1 N Λ (N ) 1 and 1 N D (N ) for large, finite N to estimate these limits. When computing Lyapunov exponents, it is implicitly assumed that the trajectory generated from the initial condition z 0 ∈ C lies on the attractor.
Recall that we must also take into account that system (1) exhibits multi-stability, and each attractor will be associated with different Lyapunov exponents. For the computations that follow, we generate a uniform grid of 25 points (x, y) ∈ (0, 10) × (0, 10) and iterate each pair of points ±x + i y a total of 10 6 times to obtain 50 initial conditions z 0 ∈ C that are assumed to lie on an attractor. We then use the above algorithm, starting from randomly selected matrices M 0 ∈ SO (2), and determine the Lyapunov exponents as the computed values at iteration N = 10 4 . We checked the accuracy of our computations against known Lyapunov exponents for the case when the attractor is a fixed or periodic point. In such cases, we find that our approximations for Λ 1 and Λ 1 + Λ 2 are correct to two or three decimal places. The approximation error is larger for a chaotic attractor, especially when the attractor covers a large portion of phase space and may not be sufficiently represented by a trajectory of N = 10 4 iterates. Nevertheless, based on the significant variation of the 50 approximations of Λ 1 and Λ 2 at parameter points for which we know that a single chaotic attractor exists, we estimate that we can determine the zero-locus of either Λ 1 or Λ 1 + Λ 2 in the (a, λ)-plane up to an error of 0.05 in the parameters. Figure 3 shows approximations of the maximal Lyapunov exponent Λ 1 computed on a square grid of 100×100 parameter pairs (a, λ); here, both a and λ take the values 0.5 + k 0.005 with k = 0, 1, . . . , 99. Panel (a) shows the largest and panel (b) the smallest values found for Λ 1 at each parameter point, where we used shades of magenta and yellow for negative and positive values, respectively, as indicated by the colour bar. Overlayed are the curve NS (dark green) of Neȋmark-Sacker bifurcation (given by a = 1 2λ ) and the curve HC (dark blue) of (first) homoclinic tangency between W s (p) and W u (p) (computed with the software Mat-ContM [24]). The parameter points on the line λ = 0.8 for the phase portraits in Figs. 1 and 2 are indicated by crosses (red). Figure 3 confirms our expectation that Λ 1 changes sign in between the curves NS, to the left of which the complex-conjugate fixed points q andq are the only attractors, and HC, which is a limit of accumulating homoclinic tangencies in the (a, λ)-plane that guarantees the onset of classical chaos [31,36].

Maximal Lyapunov exponent and classical chaos
For a more quantitative discussion of the computed values for Λ 1 , we consider again the specific transition to chaos along the line λ = 0.8. As expected, the parameter point for a = 0.6, associated with Fig. 1a, lies in the region where Λ 1 < 0. Indeed, the two attractors q andq both have the same pair of stable complexconjugate eigenvalues, with associated Lyapunov exponents Λ 1 = Λ 2 ≈ log (0.9798) ≈ −0.0204. For each of the 50 initial conditions z 0 , our numerical approximations give Λ 1 = −0.02 when rounded to two decimal places. When a = 0.63, (almost) all initial conditions converge to one of the two invariant circles, which means that Λ 1 = 0. Our computations consistently give 0 < Λ 1 ≤ 0.0005. As a increases to a = 0.65, two different complex-conjugate pairs of attractors exists, which means that Λ 1 depends on the chosen initial condition: for points ±x + i y inside an invariant circle, or just outside of it, iterates end up at initial conditions z 0 (almost) on the invariant cirlce and Λ 1 = 0, and the same happens for any preimage of such points; for (almost) all other points, the iterates end up at initial conditions z 0 (almost) on one of the complexconjugate period-8 orbits and Λ 1 ≈ −0.0342 is determined by their (equal) eigenvalues. When approximating Λ 1 by simulation, we get Λ 1 ∈ [−0.0342, 0.0004]. Hence, the pale magenta colour in Fig. 3a at this parameter point is determined by the invariant circles, while it has a darker shade in panel (b), corresponding to the period-8 attractors. When a = 0.67, the approximations give Λ 1 ∈ [0.1272, 0.1399], even though the two complex-conjugate chaotic attractors have identical maximal Lyapunov exponents. However, all values are positive, which confirms that classical chaos appears before the first homoclinic tangency between the stable and unstable sets W s (p) and W u (p) of the real saddle p. Finally, when a = 0.7, our computed values for Λ 1 range from about 0.2323 to about 0.2442, suggesting that this parameter point lies firmly in the chaotic regime.
For the five parameter points discussed above, we found that Λ 1 changes sign as a increases, but the sum

Wild chaos and associated Lyapunov exponents
Classical chaos is associated with accumulations of homoclinic tangencies [31,36] and can occur for any two-dimensional systems given by a continuously differentiable map. On the other hand, wild chaos is associated with robust homoclinic and heteroclinic tangencies and can only be found for planar maps that are non-invertible, such as the map f given in (1). As discussed in [2], wild chaos arises from interactions of the chaotic invariant set with the critical point at 0; here, the chaotic attractor is given by the closure of the unstable set W u (p). Indeed, the existence of forward and backward critical sets are the specific features that distinguish f from an invertible map.
The conjecture in [17] states that wild chaos in system (1) arises as soon as W u (p) contains the critical point 0, which occurs at a bifurcation that the authors refer to as a backward critical tangency. Similar to the homoclinic tangency, the accumulating nature of W u (p) implies that the backward critical tangency is immediately followed by an accumulation of backward critical tangencies; the first backward critical tan-gency is the limit of this sequence. We remark that the backward critical tangency can only occur for noninvertible maps. The equivalent bifurcation for a threedimensional invertible map is more complicated, involving a heteroclinic tangency between a one-dimensional unstable manifold and a one-dimensional stable manifold of a different fixed or periodic point; we refer to [3,6] for details.
A numerical characteristic for the existence of a wild chaotic attractor is that the sum of its two largest Lyapunov exponents is positive [14,37]. In general, this is not sufficient, because the linearised dynamics transverse to the attractor must also be sufficiently (dominating) contracting; we refer to [15] for precise details. For the planar example presented in this paper, this additional condition is trivially satisfied and the requirement that Λ 1 + Λ 2 > 0 is both necessary and sufficient for the existence of wild chaos. Figure 4a shows the phase portrait of system (1) with a = 0.8 and λ again fixed at 0.8. While similar to the phase portrait in Fig. 2d (with a = 0.7), the chaotic attractor (orange dots) for this case is believed to lie in the wild chaotic regime. In both Figs. 2d and 4a, the chaotic attractor is formed by the closure of the unstable set W u (p) (red curves) of the fixed point p (bright-green cross). Note that W u (p) has 'thickened' considerably in Fig. 4(a) such that many more selfintersections are seen; moreover, the stable set W s (p) (blue curves) is more prominently filling the complex plane. The essential difference, however, is that the 'first' self-intersection of W u (p), the point with the shortest arclength distance to p along W u (p), now lies to the left of 0, on the negative real axis. Indeed, several computed points of the backward critical set (green dots) now lie inside the cardioid formed by p and the two local branches of W u (p) up to the first self-intersection. Figure 4b shows the (maximal) sum Λ 1 + Λ 2 of the Lyapunov exponents computed on the same square grid of parameter pairs (a, λ) as in Fig. 3; the magenta and yellow colour shades correspond to negative and positive values, respectively, as given by the colour bar in Fig. 3. The (red) cross corresponds to the parameter point (a, λ) = (0.8, 0.8) associated with the phase portrait in panel (a); note that it lies in the yellowshaded region. Ineed, for this point we find Λ 1 + Λ 2 ∈ [0.0647, 0.0811], which confirms the wild chaotic nature of this attractor. With an adaptation of the continuation toolbox in MatContM [24] for heteroclinic tangencies, we can numerically verify that the (first) backward critical tangency for λ = 0.8 occurs when a ≈ 0.7366, in between a = 0.7 and a = 0.8. We can subsequently compute a curve of first backward tangency in the (a, λ)-plane by continuation of a (not necessarily unique) sequence of backward iterates in the backward critical set that converges to p; we label this curve BWT and it is overlayed (brown curve) in respectively. Hence, the sum Λ 1 + Λ 2 appears to be decreasing as a increases towards the curve BWT; this counter-intuitive behaviour continues until the homoclinic tangency HC at a ≈ 0.6841, after which Λ 1 + Λ 2 appears to increase monotonically. All computed values for Λ 1 + Λ 2 lie above −0.005 as soon as a ≥ 0.735, just as a crosses the curve BWT.
In [2], a complete proof of existence of wild chaos in system (1) is given for (a, λ) close to the top-right corner (1, 1). This proof does not consider W u (p) directly, because such a global invariant set can only be found numerically. Instead, the authors impose strong conditions on the eigenvalues of p to show that wild chaos certainly exists for parameters in the (a, λ)-plane close to the point (1, 1). The numerical results in [17] give a geometric perspective and show that the first backward critical tangency bounds a region in the (a, λ)-plane that contains (1, 1) and has topologically the same backbone structure of the dynamical system. Hence, the authors conjectured that the region of wild chaos is bounded by the curve BWT of (first) backward critical tangency [17,29]. Our computations of the sum Λ 1 + Λ 2 of the two Lyapunov exponents corroborates this hypothesis; see Fig. 4b. The match is particularly good along the lower segment, although roughly for λ > 0.7, the zero-locus of Λ 1 + Λ 2 deviates away from the curve BWT in the direction of increasing a. Upon closer inspection, we find that the computed values for the sum Λ 1 + Λ 2 are almost all larger than − 0.05 when (a, λ) lies to the right of the curve BWT; the different approximations obtained for each (a, λ)pair immediately to the right of BWT always lie in the range [−0.0652, 0.0033]. For the 100 data points immediately to the right of BWT, there are 17 with a minimum of Λ 1 + Λ 2 below −0.05, but only the data point (a, λ) = (0.625, 0.995) also records a maximum sum below −0.05; the first positive Lyapunov sum is detected when (a, λ) = (0.665, 0.995). We conclude that our computations confirm the hypothesis that the curve BWT bounds the region of existence for wild chaos, to within an accuracy of O(10 −1 ). On the other hand, we note that the numerical evidence becomes slightly less convincing for λ-values very close to the (limiting) case of λ = 1.

Conclusions
We presented a quantitative investigation of the existence of both classical and wild chaotic dynamics for the non-invertible discrete-time system from [2] that is defined on the punctured complex plane by the Lorenzlike double-cover map (1) with parameters a > 0 and λ ∈ (0, 1). This system is well studied and a proof was provided in [2] that wild chaos exists for (a, λ) ≈ (1, 1). The geometric characteristics of wild chaos were investigated in detail in [17,29] and a conjecture was made that the first backward critical tangency is the boundary of wild chaos. We complemented the evidence for this conjecture with the computation of the maximum and sum of the Lyapunov exponents for a broad range of parameters. We found that there is a very good match between the conjectured onset of wild chaos and the moment when the sum Λ 1 + Λ 2 of the Lyapunov exponents becomes positive: the curve BWT of first backward critical tangency in the (a, λ)-plane lies very close to and always slightly to the left of the zero-locus for Λ 1 + Λ 2 , confirming that wild chaos only arises for sufficiently large values of a and λ as given by the locus BWT.
Geometrically, the map (1) is an example of a system with hetero-chaotic attractors [35]; these are characterised by the existence of (saddle) periodic points on the attractor with different stability indices (different numbers of stable and unstable eigenvalues). In future work, we plan to investigate parameter-dependent versions of the paradigm models from [35] as other possible wild chaotic maps with the conjectured bifurcation structure. This larger class of hetero-chaotic systems also includes systems with so-called hyperchaotic attractors [1,20,21], which exhibit two positive Lyapunov exponents; such systems cannot be planar and are, therefore, computationally more challenging. It would be of interest to study the transition from chaotic to hyperchaotic attractors via a wild chaotic attractor, which may well both be organised by (different) forms of backward critical tangencies; however, this is also left for future research.
In general, very few explicit examples of wild chaotic systems exist and there are no further studies of a possible backward critical tangency causing the onset of wild chaos. The authors of [15,37] present a fourdimensional Lorenz-like vector field for which the computation of Lyapunov exponents has provided insight in the existence of wild chaos. For a four-dimensional continuous-time system, it is not sufficient to check the sign of the sum of the Lyapunov exponents; the expansion characteristic becomes more elaborate and involves the computation of invariant subspaces [15]. The backward critical tangency for non-invertible maps corresponds, in vector fields, to a particular form of heteroclinic tangency between periodic orbits that have associated unstable eigenspaces with different dimensions [17]. From a computational point of view, such a heteroclinic tangency can be found in a similar way, and we still expect that the bifurcation structure is similar, although it will be challenging to identify the periodic orbits involved in this bifurcation. We plan to investigate the geometric organisation of this four-dimensional Lorenz-like system in future work to obtain further evidence of matching geometric and expansion characteristics of wild chaos.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.