Finding first foliation tangencies in the Lorenz system

Classical studies of chaos in the well-known Lorenz system are based on reduction to the one-dimensional Lorenz map, which captures the full behavior of the dynamics of the chaotic Lorenz attractor. This reduction requires that the stable and unstable foliations on a particular Poincaré section are transverse locally near the chaotic Lorenz attractor. We study when this so-called foliation condition fails for the first time and the classic Lorenz attractor becomes a quasi-attractor. This transition is characterized by the creation of tangencies between the stable and unstable foliations and the appearance of hooked horseshoes in the Poincaré return map. We consider how the three-dimensional phase space is organized by the global invariant manifolds of saddle equilibria and saddle periodic orbits — before and after the loss of the foliation condition. We compute these global objects as families of orbit segments, which are found by setting up a suitable two-point boundary value problem (BVP). We then formulate a multi-segment BVP to find the first tangency between the stable foliation and the intersection curves in the Poincaré section of the two-dimensional unstable manifold of a periodic orbit. It is a distinct advantage of our BVP set-up that we are able to detect and readily continue the locus of first foliation tangency in any plane of two parameters as part of the overall bifurcation diagram. Our computations show that the region of existence of the classic Lorenz attractor is bounded in each parameter plane. It forms a slanted (unbounded) cone in the three-parameter space with a curve of terminal-point or T-point bifurcations on the locus of first foliation tangency; we identify the tip of this cone as a codimension-three T-point-Hopf bifurcation point, where the curve of T-point bifurcations meets a surface of Hopf bifurcation. Moreover, we are able to find other first foliation tangencies for larger values of the parameters that are associated with additional T-point bifurcations: each tangency adds an extra twist to the central region of the quasi-attractor.

with three positive dimensionless parameters: the Rayleigh number ρ, the Prandtl number σ, and the coupling strength β. System (1) is invariant under rotation by π about the z-axis, so that its dynamics are symmetric under the transformation (2) (x, y, z) → (−x, −y, z).
The origin 0 = (0, 0, 0) is always an equilibrium of system (1); it is stable for 0 < ρ < 1 and becomes a saddle in a pitchfork bifurcation at ρ = 1. For ρ > 1, the origin is hyperbolic and has two negative real eigenvalues and a third positive eigenvalue. Therefore, by the Stable Manifold Theorem [51], there exists a two-dimensional stable manifold W s (0) that consists of all trajectories that tend to 0 in forward time and a one-dimensional unstable manifold W u (0) that consists of all trajectories that tend to 0 in backward time. Under the symmetry (2), the stable manifold W s (0) is invariant, and the two branches of W u (0) map to one another. The stable manifold W s (0) plays an important role in organizing the dynamics in phase space and is referred to as the Lorenz manifold [44]. System (1) also has a pair of secondary equilibria, which exist for ρ > 1 and are each other's image under the symmetry (2). The equilibria p ± arise from the pitchfork bifurcation at ρ = 1 and are initially stable; they become saddles in a Hopf bifurcation at Hence, as saddles, their stable and unstable manifolds, denoted W s (p ± ) and W u (p ± ), have dimensions one and two, respectively. Famously, Lorenz found sensitive dependence on initial conditions for ρ = 28, σ = 10, and β = 8 3 [47,61]. At these now classical parameter values, the well-known Lorenz attractor is the only attractor.
The route to chaos and the birth of the chaotic attractor, denoted by L, has been the subject of much study, particularly for 0 < ρ ≤ 28 with σ = 10 and β = 8 3 fixed; for example, see [11,24,25,26,41,54]. In 1977, Afraimovich, Bykov, and Shilnikov studied the analogous route to chaos for β and ρ fixed at their classic values and 0 < σ < 10 [1]. The occurrence in the system of trajectories with arbitrarily long transients was investigated in 1979 by Kaplan and Yorke, who dubbed this behavior preturbulence and found it to be a precursor of the more complicated chaotic dynamics discovered by Lorenz [41]. The concept of turbulence, meaning irregular and complicated chaotic behavior, was coined by Ruelle and Takens in the 1970s in their study of viscous fluids [55]. They indicated that this behavior could be the result of the presence of a strange attractor and went on to investigate the relationship between turbulence and the Lorenz attractor in [54].
The structure of the strange attractor in the Lorenz system was further investigated by Guckenheimer and Williams [33,65]. In their work and independently in [1], one finds the first geometric description of the dynamics on the Lorenz attractor and the introduction of the geometric Lorenz model. The analysis is done via reduction to a one-dimensional map with equivalent behavior. This reduction is based on the two-dimensional (invertible) Poincaré return map that is defined on the planar section Σ ρ := {z = ρ − 1} (chosen to contain the equilibria p ± ) in the region where the flow is down. The behavior of contraction and expansion on Σ ρ is organized by so-called stable and unstable foliations, each consisting of curve segments, or leaves, that are mapped in forward or backward time to other leaves in the foliation; in particular, L is contained in the unstable foliation. The one-dimensional map, called the Lorenz map [34], is obtained by projection along leaves of the stable foliation. The Lorenz map completely represents the dynamics of L provided that the leaves of the stable foliation are locally transverse to L, they map injectively under the flow, and the dynamics on the stable leaves is a contraction toward L; this is called the foliation condition.
The foliation condition is satisfied for the geometric Lorenz model [1,2,33,35,65], and the dynamics on L are accurately described by the one-dimensional Lorenz map representing the dynamics from leaf to leaf in the stable foliation. Tucker [61] proved in 1999 that the Lorenz system (1) for the classical parameter values satisfies the conditions of the geometric Lorenz model; hence, the foliation condition holds, and the chaotic dynamics are described by the one-dimensional Lorenz map. It is generally believed that this is also true for ρ ≤ 30 with σ = 10 and β = 8 3 . However, for ρ past 30, there is a first tangency between the stable and unstable foliations, and the foliation condition fails [17,37,38,60]. As a result, for larger ρ-values, the one-dimensional map no longer captures the full dynamics on L, and much less is known about the dynamics of the Lorenz system.
In this paper, we build on the recent approach taken in [23,25], where the computation of stable and unstable manifolds, via the continuation of suitable two-point boundary value problems (BVPs), revealed the geometric mechanisms in the full three-dimensional phase space that are behind the transition to chaotic dynamics in the Lorenz system (1). We are interested here in the precise moment when the foliation condition fails. The only specific numerical estimate available for the loss of foliation condition was computed by Sparrow [60] as lying in the interval ρ ∈ [30.1, 30.2] for fixed σ = 10 and β = 8 3 . He found this interval by using numerical integration to determine how vectors in Σ ρ align with the attractor on return to the section. Furthermore, Bykov and Shilnikov [16] computed, for fixed β = 8 3 , a locus in the (ρ, σ)-plane along which L changes to a nonorientable or quasi-attractor ; see also the translation [17]. A quasi-attractor is a complex limit set that contains a dense set of stable periodic orbits with narrow basins of attraction [2,58]. The quasi-attractor contains infinitely many windows of stability related to the creation and destruction of stable periodic orbits; hence, the dynamics in this region are far more complicated. This topological change corresponds to the loss of the foliation condition, and it was detected by Bykov and Shilnikov as a change of orientability of homoclinic orbits. We refer to subsection 5.1 for details on their computational approach and Figure 9(b) for a reproduction of their sketch [17, Figure 3]. Although no precise values are given there, the value for σ = 10 appears to be ρ ≈ 31.0.
Our computational approach to investigating the loss of the foliation condition is to identify and calculate the onset of a tangency between the stable and unstable foliations of the Poincaré section Σ ρ , which we call a first foliation tangency and denote by F 1 . We approximate leaves of the stable and unstable foliations in Σ ρ directly as intersection curves with Σ ρ of two-dimensional manifolds of equilibria and periodic orbits. We detect the parameter values at which the first foliation tangency occurs by solving BVPs with the software package Auto [20,21]; for the general theory, see [22,44]. The BVP setup has the particular benefit that it is straightforward to continue specific families of solutions and, hence, the locus of the loss of the foliation condition in any pair of the system parameters.
More specifically, our results are as follows. Our estimate for the value of the loss of the foliation condition when σ = 10 and β = 8 3 is ρ ≈ 31.01, which is correct to two decimal places. We also compute the locus F 1 in the (ρ, σ)-plane for fixed β = 8 3 and find good agreement with the sketch by Bykov and Shilnikov. We find that F 1 is a smooth curve in this plane with a codimension-two terminal-point, or T-point, also known as a Bykov cycle [18,31,52]. At this T-point, the one-dimensional unstable manifold W u (0) is contained in the one-dimensional stable manifolds W s (p ± ), creating a heteroclinic connection from 0 to p ± ; moreover, the two-dimensional manifolds W u (p ± ) intersect W s (0) transversally, forming a heteroclinic connection from p ± to 0. There are many other such T-points in the Lorenz system, corresponding to increasingly complicated windings of W u (0) around p ± before connecting to p ± ; see [18] for details. The bifurcation structure near a T-point has been studied, for example, in [40,43,59,66] but not in the context of the loss of the foliation condition. The bifurcation curves of the main homoclinic and heteroclinic bifurcations of p ± both terminate at a T-point. These bifurcation curves lie extremely close to each other in the (ρ, σ)-plane for fixed β = 8 3 . As was also shown in [18], each T-point in the Lorenz system is associated with a phenomenon called an α-flip bifurcation, where the α-limits of the respective branches of W s (p ± ) switch sides. The first, or principal, T-point was discovered by Petrovskaya and Yudovich [52] in 1980 and independently by Alsfen and Frøyland [4] in 1985. The first 25 T-points in the Lorenz system are shown in [18] in relation to the homoclinic and heteroclinic bifurcations of p ± and the α-flip bifurcation of W s (p ± ).
The chaotic attractor L is created in a so-called EtoP bifurcation, at which W u (0) connects to a (symmetric) pair of periodic orbits; we refer to section 2 and [1,41,55] for details. We find that the bifurcation curve EtoP and the locus F 1 have two intersection points, meaning that the region of existence of the Lorenz attractor L is bounded in the (ρ, σ)-plane for fixed β = 8 3 . This discovery already goes beyond the work of Bykov and Shilnikov, who found only one intersection point. Moreover, we determine that the region of existence of L is also bounded in the (ρ, β)-and (σ, β)-planes, with σ = 10 and ρ = 28, respectively. By computing the curves EtoP and F 1 in the (ρ, σ)-plane for different values of β, we show that the region of existence of L has the shape of a slanted cone in the three-dimensional (ρ, σ, β)-space, which is bounded by surfaces EtoP and F 1 . There exists a curve T 1 of principal T-points on the surface F 1 that ends at the bottom or tip of the cone in a degenerate T-point called a T-point-Hopf (TH) bifurcation; at this point, the equilibria p ± involved in the T-point bifurcation undergo a Hopf bifurcation. The codimension-three bifurcation point TH has been studied in [5,28] in general terms as an organizing center for homoclinic bifurcations of the origin in the Lorenz system; however, its role for the loss of the foliation condition has not been recognized.
After the first foliation tangency, along the locus F 1 , the stable and unstable foliations are no longer transverse, and the foliation condition is lost. We find that, as ρ and σ are increased, additional first foliation tangencies occur, and we adapt our BVP setup to detect them. Specifically, we detect a further two first foliation tangencies and compute their loci F 2 and F 3 in the (ρ, σ)-plane for fixed β = 8 3 . We observe that the loci F 2 and F 3 also each have T-points T 2 and T 3 on them, respectively; see Figure 12. In [13], Barrio, Shilnikov, and Shilnikov computed a color code for the kneading sequence based on the number of revolutions of W u (0) around p ± . Their image of the (ρ, σ)-plane with β = 8 3 shows regions of stability and the structure of homoclinic orbits around T-points. Our computation of homoclinic bifurcation curves and boundaries of the regions of stability in the (ρ, σ)-plane are in good agreement with their findings, and we discuss the loci of the first foliation tangency in this context. Finally, we show how each additional first foliation tangency manifests itself as additional folds in the unstable foliation.
The layout of this paper is as follows. In the next section, we briefly review the bifurcations involved in the route to chaos and the birth of L, and we discuss known properties of the Lorenz map for parameters outside the region where the foliation condition is satisfied. We illustrate transverse and nontransverse foliations in Σ ρ in section 3 and present more specific evidence for the existence of a foliation tangency in subsection 3.2. In section 4, we formulate and implement as a BVP the computation of the first foliation tangency associated with the loss of the foliation condition and determine its ρ-value for σ = 10 and β = 8 3 fixed. In section 5, we compute bifurcation diagrams near the principal T-point in the planes defined by each possible parameter pair, and in subsection 5.2, we present the locus F 1 of the first foliation tangency in the full three-parameter space. Subsection 5.3 presents two additional loci of first foliation tangencies in the (ρ, σ)-plane with β = 8 3 and discusses their manifestation in terms of W u (Γ rl ). We finish with a discussion in section 6. The list the parameters used for the BVP computations and further computational details can be found in Appendix A.
2. Review of the bifurcations en route to and beyond the creation of L. We briefly recall the main bifurcations that occur in the Lorenz system (1) in the transition to chaos and the birth of L; see also [23,25] and references therein. Here, we assume that σ = 10 and β = 8 3 are fixed at their classical values and that only ρ is varied. As mentioned in the introduction, the origin 0 is stable for 0 < ρ < 1; at ρ = 1, a pitchfork bifurcation occurs that gives rise to two stable equilibria p ± that are rendered saddles in a Hopf bifurcation at ρ = ρ H ≈ 24.7368. A homoclinic bifurcation of 0 occurs at ρ ≈ 13.9265 when W u (0) is contained in W s (0) and consists of a loop h r and its symmetric counterpart h l around p + and p − , respectively; we will use h r to denote this homoclinic bifurcation throughout. The homoclinic bifurcation h r generates a chaotic saddle S containing infinitely many saddle periodic orbits; this bifurcation is also called the homoclinic explosion. In particular, the homoclinic explosion creates the symmetric pair of periodic orbits Γ + and Γ − , which loop once around p + and p − , respectively. The periodic orbits Γ ± have twodimensional stable manifolds W s (Γ ± ) and two-dimensional unstable manifolds W u (Γ ± ). A heteroclinic connection from 0 to Γ ± forms when W u (0) is contained in W s (Γ ± ) at ρ ≈ 24.0579; this is the EtoP connection referred to in the introduction. The parameter interval between h r and EtoP is the preturbulent regime identified by Kaplan and York [41,42]. In this regime, the equilibria p ± remain the only attractors, but trajectories with arbitrarily long chaotic transients can be found. The complicated dynamics can be attributed to the Cantorlike structure of the stable manifold of S [25]. The EtoP bifurcation creates L [1,41,55], which initially coexists with the attractors p ± until their stability changes at the Hopf bifurcation. Except for Γ ± , the saddle periodic orbits created in the homoclinic explosion persist through the Hopf bifurcation and are dense in L [23,47,63].
The classic Lorenz attractor L is an attracting set in the Lorenz system that consists of an infinite number of two-dimensional surfaces that join along the one-dimensional unstable manifold of the origin W u (0) [1,34,47]. A sketch of L was already presented in Lorenz's original paper [47]. The point 0 and the unstable manifold W u (0) are contained in L [1,65], and W u (0) forms a natural boundary of L. Furthermore, L does not enter a certain neighborhood of p ± at the classic parameter values. Points on the one-dimensional stable manifolds W s (p ± ) are the only points in phase space that do not tend to L. Since saddle periodic orbits are dense in L [47,63], their unstable manifolds accumulate on L due to the stretching nature of the dynamics. In fact, L contains the sheets of all the two-dimensional unstable manifolds of the periodic orbits, which lie very close together in phase space due to the strong contraction in the system, and connect up along W u (0) [33]. Therefore, the twodimensional unstable manifold of any saddle periodic orbit gives a good representation of L.
The geometric Lorenz model and the associated Lorenz map provide a way of accurately describing the dynamics on L provided that the technical conditions of the reduction hold. In particular, each point of each intersection set of L with Σ ρ needs to intersect exactly one leaf (locally) of the stable foliation [34]. Numerical investigations by Bykov and Shilnikov [17], Hénon and Pomeau [37,38], and Sparrow [60] show that the two-dimensional return map appears to develop hooked horseshoes that persist for large values of ρ. Then the above oneto-one correspondence in Σ ρ is violated, the foliation condition is no longer satisfied, and, therefore, L is no longer completely represented by the geometric Lorenz model or the Lorenz map. The one-dimensional Lorenz map develops additional maxima and minima that enter through the boundaries of the interval. Luzzatto and Viana studied Lorenz-like families of one-dimensional maps with such extra critical points in an attempt to describe the change in dynamics before and after the loss of the foliation condition [48,49]; however, they do not use them to find the moment of the loss of the foliation condition. Hao, Liu, and Zheng created a series of one-dimensional maps of the Lorenz system, called first-return maps, for 28 < ρ < 200 [36]. Their first-return maps develop additional maxima and minima as ρ is increased, corresponding to tangency points between the stable and unstable foliations in Σ ρ . The first of these tangencies appears to correspond to the loss of the foliation condition. Although Hao, Liu, and Zheng show some curves of the stable and unstable foliations of Σ ρ , they do not use these to identify the loss of the foliation condition.
Dullin et al. show diagrams of the (ρ, σ)-plane with β = 8 3 , for the range (ρ, σ) ∈ [0, 2000] × [0, 2000], that are colored according to the asymptotic behavior of trajectories in phase space [26]. These diagrams clearly show a layered periodic pattern due to the alternating stability of symmetric and nonsymmetric periodic orbits. Barrio and Serrano use maximum Lyapunov exponents and look for chaotic orbits at points in the planes given by each pair of parameters ρ, σ and β [11,12]. Their computations show the regions in the parameter planes where chaos exists. The scan methods of Dullin et al. [26] and Barrio and Serrano [11,12] do not give any indication of a topological change in L.
3. Transverse and nontransverse foliations in Σ ρ . The stable and unstable foliations in Σ ρ are not known analytically and can only be approximated numerically. We compute them as intersection curves of stable and unstable manifolds of periodic orbits and equilibria for fixed σ = 10 and β = 8 3 . Our computations are based on continuation of a two-point boundary value problem, and we refer to [18,24,25,44] for more information on the setup.
3.1. Transverse foliations for ρ = 28. As described in section 2, the chaotic attractor L consists of infinitely many surfaces each of which intersects Σ ρ in a one-dimensional curve that is part of one of the leaves of the unstable foliation [32]. Furthermore, the two-dimensional unstable manifold of any periodic orbit lies dense in L. Therefore, its intersection curves with Σ ρ form a dense subset of the intersection curves of L with Σ ρ . We denote each periodic orbit in L by Γ s 1 ,...,sn , where s 1 , . . . , s n is the symbol sequence of the orbit such that s i = r indicates a loop around p + and s i = l indicates a loop around p − . Throughout this paper, we work mainly with the topologically simplest periodic orbit in L, the symmetric periodic orbit Γ rl that loops once around p + and then once around p − ; however, we also consider Γ rll to confirm our calculations. The periodic orbit Γ rll loops once around p + and then twice around p − ; hence, it is not symmetric and coexists with its symmetric counterpart Γ lrr . Viswanath [63] uses a modified Newton method to find periodic orbits in L up to a high period and provides a list of initial conditions in Σ ρ for each periodic orbit to high precision (14 significant figures). We take the appropriate initial condition for Γ rl (or Γ rll ) from this list as sufficiently accurate for the computation and continuation with Auto [20,21]. More recently Barrio, Dena, and Tucker [14] created a comprehensive database of periodic orbits in the Lorenz system, with initial conditions verified by interval-arithmetic techniques up to an accuracy of 1000 significant figures.
We approximate part of W u (Γ rl ) for ρ = 28 as a one-parameter family of orbit segments that start close to Γ rl in its linear unstable direction and are computed up to a maximum chosen integration time τ max . We then compute a sufficient number of the intersection curves directly as end points of a one-parameter family of solutions of a suitable BVP. The intersection curves W u (Γ rl ) are part of the unstable foliation.
The unstable manifold W u (0) forms the outer boundary to W u (Γ rl ). We approximate W u (0) as two orbit segments that are computed up to a maximum arclength L max from their starting points near 0 on the linear unstable direction v u . We compute the intersection points Values of τ max , L max and other accuracy parameters used to create each figure are given in Appendix A. Figure 1 shows L for ρ = 28 as represented by W u (Γ rl ), the intersection curves W u (Γ rl ), and as a sketch of the one-dimensional Lorenz map. Panel (a) shows the computed part of the two-dimensional surface W u (Γ rl ) with Γ rl and W u (0). The classic section Σ ρ through p ± is also shown. Panel (b) shows curves in W u (Γ rl ); in this and subsequent similar panels, we rotate Σ ρ by π 2 so the horizontal axis is the diagonal where x = y. The four intersection points in γ rl := Γ rl ∩ Σ ρ lie on W u (Γ rl ) and are shown in panel (b). Also shown in panel (b) are the first eight intersection points of W u (0), which correspond to the end points of the main segments of W u (Γ rl ). Furthermore, panel (b) shows the tangency locus C that contains p ± and separates Σ ρ into regions where the flow is upward ( ) and downward (⊗) [46]. A sketch of the one-dimensional Lorenz map is shown in Figure 1(c) with the line of discontinuity at Figure 1. Representation of the classic Lorenz attractor L for ρ = 28, σ = 10, and β = 8 3 . Panel (a) shows W u (Γ rl ) with the periodic orbit Γ rl (yellow), the first part of W u (0) (brown), and the Poincaré section Σρ (green). The intersection curves W u (Γ rl ) (red) in the plane Σρ are shown in panel (b), together with the four intersection points of γ rl (yellow) and the first eight of W u (0) (brown). Also shown are p ± and the tangency locus C (gray) that divides Σρ into regions where the flow is up and where it is down ⊗; note that Σρ has been rotated so the diagonal (where x = y) is the horizontal axis. Panel (c) shows a sketch of the one-dimensional Lorenz map with the vertical line of discontinuity at 0 (blue). 0 corresponding to points that lie on W s (0) and, hence, do not return to the plane Σ ρ . The boundaries ±a = ± β(ρ − 1) of the Lorenz map are the x-components of p ± , respectively. The empty square regions in the corners of panel (c) correspond to the local regions around p ± that are not part of L; see also [34].
On the level of Figure 1(b), it appears that W u (Γ rl ) consists of only four disjoint curves. In fact, the intersection of the Lorenz attractor with the Poincaré section has a fractal structure [47], but the distance between the curves is less than 10 −5 [64] because of the strong contraction toward the Lorenz attractor L. We actually computed and plot 60 intersection curves, which cannot be distinguished at the scale of panel (b). Further points in W u (0), beyond the first eight that are shown, correspond to end points of the shorter curves of W u (Γ rl ). The intersection set W u (Γ rl ) does not cross C when ρ = 28, and so the curves of W u (Γ rl ) in the outer regions, where the flow is up, are the images under the flow of the central curves of W u (Γ rl ), in which the flow points down. Historically, the Lorenz map is constructed from the central region of Σ ρ , where the flow is down, in between the two curves that form the tangency locus C. Figure 2 compares W u (Γ rl ) in panel (a) with two alternative representations of L, generated by the unstable manifolds W u (Γ rll ) of Γ rll in panel (b) and W u (p + ) of the equilibrium p + in panel (c); all these panels show the projection onto the (x, z)-plane along the y-direction. Note that W u (p + ) only accumulates on L because the local manifold W u loc (p + ) near p + is not part of L. Panel (d) shows all three manifolds together with the plane Σ ρ in the three-dimensional (x, y, z)-space. Also shown in each panel are the first parts of the two symmetrically related branches of W u (0), which can be seen to bound W u (Γ rl ), W u (Γ rll ) and W u (p + ).
Each of the invariant manifolds shown in Figure 2 was computed as a separate family of orbit segments that start in the linear unstable direction of the periodic orbit or equilibrium. Due to the strong contraction in the Lorenz system, it is impossible to distinguish these unstable manifolds at the level of the surfaces shown in panel (d); however, they can be distinguished by the associated families of orbit segments. The close proximity of the computed surfaces in phase space constitutes a confirmation of the accuracy of our computations.
The theory also tells us that the stable foliation contains the intersection curves of the two-dimensional stable manifolds of all equilibria and periodic orbits [36]. To find curves of the stable foliation, we consider the intersection curves W s (0) := W s (0) ∩ Σ ρ of the two-dimensional Lorenz manifold W s (0) with Σ ρ . The Lorenz manifold lies dense in phase space due to the sensitivity on initial conditions [61], and therefore its intersection curves with Σ ρ form a dense subset of the curves of the stable foliation [36]. To compute these intersection curves, we start with the same approach as in [18,24], where W s (0) is computed inside a suitably large sphere S R with radius R. Here, we choose R = 70 to ensure that the entire attractor for ρ = 28 and, hence, W u (Γ rl ) is fully contained in S R . The sphere is centered at (0, 0, ρ − 1), and therefore Σ ρ bisects S R at the equator into two equal halves. We first compute an orbit segment with one end point on an ellipse in the stable eigenspace of the origin and the other end point on S R . We then continue the orbit segment around the ellipse and detect each instance that the end point on S R intersects Σ ρ at the equator. For each of these detected solutions, we then constrain the end point to lie on Σ ρ and continue around the ellipse to generate a family of orbit segments whose end points in Σ ρ trace out curves in W s (0). Due to the finite-time nature of our computations, we can compute only a finite number of curves in W s (0). Figure 2. Representation of L for ρ = 28, σ = 10, and β = 8 3 by three different unstable manifolds (red), namely, of Γ rl (yellow) (a), Γ rll (orange) (b), and p + (black) (c). Panel (d) shows all three unstable manifolds together with Σρ (green). Also shown in all panels is the first piece of W u (0) (brown). Figure 3. Section Σρ and the stable manifold W s (0) for ρ = 28, σ = 10, and β = 8 3 . Panel (a) shows W s (0) (blue) computed inside the sphere SR for R = 70 with Σρ (green). The surface W s (0) is rendered transparent above Σρ and solid below, and the computed part of W s (0) ∩ SR has been highlighted in black. Panel (b) shows curves of W s (0) (blue) and W u (Γ rl ) (red) and points of γ rl (yellow) and W u (0) (brown) within the circle with radius R = 70 around (0, 0, ρ − 1) in Σρ. Panel (c) is an enlargement of W s (0), W u (Γ rl ), γ rl , and W u (0) in the rotated view of Σρ; compare with Figure 1(b). Figure 3 shows the computed part of the Lorenz manifold W s (0) and its intersection curves with Σ ρ . Panel (a) shows the initial piece of W s (0) within S R for R = 70, together with Σ ρ in the three-dimensional phase space. Here, W s (0) is rendered as a solid blue surface beneath Σ ρ and a transparent blue surface above it. The equator of S R is the black circle on Σ ρ . The intersection curves W s (0) in Σ ρ lie inside the equator and are shown as blue curves in each panel. Panel (b) also shows W u (Γ rl ), and panel (c) is an enlargement near the curves W u (Γ rl ) in the rotated view. Panels (b) and (c) also show the sets of points p ± , γ rl , and W u (0) and the tangency locus C; compare with Figure 1(b). Figure 3(c) shows that, locally near the attractor, the computed intersection curves in W s (0) and W u (Γ rl ) are transverse, meaning that each curve in W s (0) intersects each curve in W u (Γ rl ) exactly once near L. Figure 3 suggests that the foliation condition is satisfied, which is indeed the case for ρ = 28 [17,35,60].
3.2. Nontransverse foliations for ρ = 60. Next, we consider the intersection curves W u (Γ rl ) and W s (0) in Σ ρ for ρ = 60. This value of ρ was chosen well past the estimated ρ-value at which the foliation condition is lost. Figure 4 shows W u (Γ rl ) and the corresponding intersection curves W u (Γ rl ) and W s (0). Panel (a) is a view of W u (Γ rl ) with Σ ρ in the threedimensional phase space; compare with Figure 1(a). Panel (b) shows curves in W u (Γ rl ) ⊂ Σ ρ with the points γ rl that lie on W u (Γ rl ) and the points of W u (0) that lie at the end of the visible curve segments. In panel (b), W u (Γ rl ) appears to consist of two disjoint curve segments, but again 60 curves have actually been plotted. The points p ± and locus C are also shown in panel (b). The sheets of the two-dimensional manifold W u (Γ rl ) have been pulled up and folded back near p ± in Figure 4(a); compare with Figure 1(a). The boundary of W u (Γ rl ) is still formed by W u (0), but it is no longer the curve that lies closest to p ± .
The structure of the intersection curves W u (Γ rl ) in Figure 4(b) is also different from those in Figure 3(c). For ρ = 60, some curves in W u (Γ rl ) cross the tangency locus C. Intersections of W u (Γ rl ) and C are brought about by a tangency of W u (0) with Σ ρ that, by definition, occurs on C. The tangency occurs at ρ C ≈ 30.4318, which we found by continuation of the first-return point of W u (0) to Σ ρ in ρ until a fold was detected. The tangency occurs near but does not correspond to the loss of the foliation condition [24]. The second and more important difference between Figures 3(c) and 4(b) is that some of the curves of W s (0) in the central region of Σ ρ in Figure 4(b) now intersect curves of W u (Γ rl ) twice. Moreover, there are points of (approximate) tangency between curves of W u (Γ rl ) and curves of W s (0). Note that W s (0) does not undergo a dramatic change between Figures 3(c) and 4(b).
Indeed, Figure 4(b) illustrates that W u (Γ rl ) and W s (0) are not locally transverse for ρ = 60. Rather, there are structurally stable, quadratic tangencies arbitrarily close to L, which means that the foliation condition does not hold. More specifically, each curve of W u (Γ rl ) is diffeomorphic to another curve of W u (Γ rl ) under the flow (1). Hence, each curve of W u (Γ rl ) has a point of tangency with W s (0); most tangency points are impossible to distinguish by eye due to the strong contraction of the system. Indeed, the sheets of the twodimensional unstable manifolds of periodic orbits in L connect along W u (0) [34], and so all intersection curves of these unstable manifolds with Σ ρ have a point of tangency with W s (0).
4. The first foliation tangency. We wish to find the ρ-value that corresponds to the onset of quadratic tangencies between W u (Γ rl ) and W s (0); we call this onset a first foliation tangency bifurcation. From the estimates of Sparrow [60] and Bykov and Shilnikov [17], the foliation condition fails between ρ = 30 and ρ = 32. Figure 5 shows the stable and unstable foliations for ρ = 30 and ρ = 32, just before and just after the loss of the foliation condition, respectively. Panel (a1) shows computed curves in W u (Γ rl ) and W s (0) in Σ ρ for ρ = 30, and panel (b1) shows them for ρ = 32. Since p + is the image of p − under the symmetry (2), it suffices to consider a neighborhood of, say, p + . Panels (a2) and (b2) are enlargements that show W u (Γ rl ) near p + for ρ = 30, and ρ = 32, respectively, and just two curves of W s (0) in   the bottom-left corner of each panel. Points of W u (0) are shown at the end of the main curve segments as before; the point q u ∈ W u (0) is the second intersection point of W u (0) with Σ ρ . The tangency locus C is shown in each panel as before. Figure 5 shows that W u (Γ rl ) comes very close to p ± for values of ρ near the loss of the foliation condition. Unfortunately, the computed intersections curves in W s (0) do not come sufficiently close to p + to allow for the comparison with W u (Γ rl ), also not for much higher integration times. Since W s (0) winds around W s (p + ) without intersecting it, we approximate the direction of W s (0) locally near p + by the stable eigenspace of p ± projected orthogonally onto Σ ρ ; that is, we define The direction vector V can readily be computed and is shown as the light-blue curve through p ± and q u in Figure 5. Observe that the curves in W s (0) closest to p ± appear to be parallel to V for all practical purposes. A qualitative change in W u (Γ rl ) with respect to V occurs near the equilibria p ± for ρ between ρ = 30 and ρ = 32; compare panels (a2) and (b2) of Figure 5. Figure 5(a2) shows that the direction V at q u has only one intersection with each of the curves in W u (Γ rl ), which indicates that the stable and unstable foliations are transverse to each other and the foliation condition holds for ρ = 30. Figure 5(b2) shows that the direction V at q u has additional intersections with curves in W u (Γ rl ) for ρ = 32. The curves in W u (Γ rl ) have formed visible hooks in Figure 5(b2). This means that the foliation condition no longer holds.

Characterization of the first foliation tangency.
To detect numerically the first foliation tangency, we select the curve from W u (Γ rl ) that ends at the point q u and corresponds to orbit segments on W u (Γ rl ) with the fewest returns to Σ ρ ; we denote this curve by W u F (Γ rl ). We define the direction vector Z as the unit vector tangent to W u F (Γ rl ) at q u . The direction vectors V , defined above in (4), and Z both depend on ρ; indeed, the angle between Z and V has a regular sign change between ρ = 30 and ρ = 32. We define the first foliation tangency as the moment at which the angle is zero.
Our approach of detecting a zero angle between V and Z is based on the continuation of multiple orbit segments and motivated by Lin's method [45]. The multi-segment BVP setup for the first foliation tangency is well defined and has user-specified accuracy parameters to allow for the convenient continuation of its solutions in any of the parameters. This BVP can be set up readily for invariant objects other than Γ rl , and we consider also the two-dimensional unstable manifolds W u (p + ) of the equilibrium p + and W u (Γ rll ) of the nonsymmetric periodic orbit Γ rll ; a detailed comparison of the different approximations can be found in Appendix B, Table 2. We also checked that V is indeed a good approximation to the stable foliation in the vicinity of p + and refer to Table 3 in Appendix B for details.
We calculate an approximation of Z numerically as a secant in the following way. First, we calculate an orbit segment u 0 (t) that lies on W u (0) with one end point in the unstable eigendirection of 0 and the other end point at q u in Σ ρ . We then consider a circle C κ ⊂ Σ ρ of radius κ centered at q u and compute a second orbit segment, u Γ (t) on W u (Γ rl ), with one end point on the unstable eigendirection of Γ rl and the other end point on C κ ; we refer to the latter point on W u F (Γ rl ) as q κ . The vector direction Z is then approximated by the normalized unit vector The direction Z κ is a good approximation of Z provided that κ = |q u − q κ | is sufficiently small. In our calculations, we choose Z κ with κ = 10 −4 . Full implementation details are given in subsection 4.2. Figure 6. Setup of the multi-segment BVP at ρ = 32 with the orbit segments u 0 (t) ⊂ W u (0) and u Γ (t) ⊂ W u (Γ rl ) that determine the direction Zκ (orange); also shown are the periodic orbit Γ rl (yellow) and the direction vector V (light blue) at p ± . Figure 6 shows the setup of the multi-segment BVP used to compute Z κ for ρ = 32; here, κ was chosen very large to illustrate the different components. The orbit segment u 0 (t) on W u (0) (dark red) starts near 0, and its other end point is q u in Σ ρ , which is not labeled but lies quite close to p + . The orbit segment u Γ (t) on W u (Γ rl ) (red) starts near Γ rl , and its other end point is q κ in Σ ρ , at distance κ from q u . The secant Z κ is shown in orange between q u and q κ , and the direction V in Σ ρ is shown at p ± .
We continue the multi-segment BVP that defines Z κ (for sufficiently small, fixed κ) in ρ; the effect is that q κ moves along the circle C κ . We monitor the test function where V ⊥ ∈ Σ ρ is perpendicular to V ; initially, σ = 10 and β = 8 3 are fixed. The function φ κ (ρ) is a regular test function in ρ, meaning that it is smooth and, generically, its roots are regular and isolated. During the continuation of the BVP in the direction of decreasing ρ, a root ρ F κ of the function φ κ will be detected. Since Z κ → Z, the root ρ F κ → ρ F as κ → 0; for convenience of notation, we write ρ F instead of ρ F κ from now on. Figure 7 illustrates that the detection of a regular zero of φ κ corresponds to a tangency between W u F (Γ rl ) and V at q u . The curve W u F (Γ rl ) and the direction V are shown for ρ = 30 in column (a), for ρ = ρ F ≈ 31.01 in column (b), and for ρ = 32 in column (c). Panels (a1)-(c1) show the situation locally near p + , and panels (a2)-(c2) show substantial enlargements of the y (a1) circle C κ with κ = 10 −4 centered at q u . The point q κ ∈ W u F (Γ rl ) ∩ C κ defines the secant Z κ . Figure 7(b1) corresponds to the moment of first tangency between W u (Γ rl ) and V , which can be confirmed in the enlargement in panel (b2) that shows the alignment of Z κ and V , with W u F (Γ rl ) tangent to V in very good approximation. 4.2. Implementation as a boundary value problem. We now give the specific boundary conditions for the computation of the multi-segment BVP to find the point of tangency between the stable and unstable foliations in Σ ρ ; for the general theory, we refer to [27,45]. Without loss of generality and for ease of notation, we describe our implementation for the unstable manifold W u (Γ) of a periodic orbit Γ.
In all computations, we use the time-rescaled Lorenz system where f is given by (1). A solution u(t) of (7) is an orbit segment on the time interval t ∈ [0, 1] with end points u(0) and u(1). This orbit segment is also a solution of (1) but with respect to the unscaled time t τ , where τ is the total integration time for (1). In our formulation, τ is positive when computing unstable manifolds and negative when computing stable manifolds. For fixed ρ = 32, σ = 10, and β = 8 3 , we compute an orbit segment u 0 (t) ⊂ W u (0) satisfying (7) by continuation in τ subject to the boundary condition where v u is the normalized unstable eigenvector of 0 and δ 0 = 10 −7 . We monitor the zcomponent u 0 z (1) of u 0 (1). Whenever u 0 z (1) = ρ − 1, the end point u 0 (1) lies in Σ ρ . We choose the end point that corresponds to q u ∈ Σ ρ and impose the boundary condition (9) u 0 z (1) − (ρ − 1) = 0, which constrains u 0 z (1) to lie in Σ ρ . We compute a second orbit segment u Γ (t) on W u (Γ) in the following way, where we assume that we have available Γ and its unstable Floquet bundle (as solutions of an additional BVP [44,45]). Since trajectories on W u (Γ) spiral away from Γ, it suffices to consider only one normalized vector v u (γ) of the unstable Floquet bundle. For convenience, we choose γ as the intersection point of Γ with Σ ρ . The point w u 0 is defined as where δ Γ is a fixed distance from Γ along vector v u (γ). We define the point w u 1 as the first return of the trajectory through w u 0 to the local planar section spanned by v u (γ) and a vector v Σ ∈ Σ ρ that is perpendicular to v u (γ). The line segment between w u 0 and w u 1 can be viewed as an approximate fundamental domain for W u (Γ), and we consider the boundary condition (11) u where ζ ∈ [0, 1). The orbit segment u Γ (t) is computed as satisfying (7) and (11), and we again monitor its z-component u Γ z (1) to detect the returns to Σ ρ . We choose the point u Γ (1) = q κ that lies on W u F (Γ) near p + and corresponds to the orbit segment with the lowest number of returns to Σ ρ , namely, 16 returns. We then impose the boundary condition (12) u Γ z (1) − (ρ − 1) = 0. The distance between the two points u 0 (1) and u Γ (1) in Σ ρ is initially large, as in Figure 6. We continue the BVP for u Γ satisfying (7) and boundary conditions (11) and (12) with ζ as a free continuation parameter and monitor We then impose the additional boundary condition (13) |q u − q κ | − κ = 0 and continue the BVP for u 0 and u Γ with (8), (9), and (11)-(13), where the system parameter ρ is decreasing. We detect ρ F as the moment when (6) is zero, which means that Z κ is aligned with V .
The radius κ of C κ determines the accuracy of Z κ with respect to the actual tangent Z of W u F (Γ) at the end point q u . The difference between Z κ and Z is O(κ 2 ) since W u F (Γ) is smooth. On the other hand, there is a trade-off, and κ should not be chosen too small: as κ decreases, the point q κ approaches the point q u , which implies that u Γ (t) passes very close to 0. In the limit of κ → 0, the orbit segment u Γ (t) connects to 0, and the integration time τ goes to ∞. Numerically, this means that it is difficult to calculate u Γ (t) accurately for very small κ. We found that κ = 10 −4 is the most suitable value for the calculations we present here, and we have confidence in our value of ρ F ≈ 31.01 to two decimal places. This statement is corroborated by performing the above computation of ρ F for each of the manifolds W u (Γ rl ), W u (Γ rll ), and W u (p + ) and for several choices of κ. We found that the ρ F -values for each κ only change in the fourth decimal place; see Appendix B and the data in Table 2 for details.
5. The locus of the first foliation tangency. The advantage of our method is that, once a zero of (6) has been detected, it can be continued in any pair of the system parameters. To this end, we impose the additional boundary condition (14) V ⊥ · Z κ = 0 and continue the BVP for u 0 and u Γ with (8), (9), and (11)-(14) with either ρ and σ, ρ and β, or σ and β as continuation parameters. In this way, we find the locus F 1 of the loss of the foliation condition in any two-parameter plane. When starting the computation of F 1 in the (ρ, σ)-plane from the detected value of ρ F for σ = 10, we find that ρ decreases and σ increases. As we continue F 1 in the direction of decreasing ρ, the point q u moves closer to p + in Σ ρ , and the integration time τ increases dramatically. In the limit when q u coincides with p + , the computed branch of the one-dimensional unstable manifold W u (0) coincides with the respective branch of the one-dimensional stable manifold W s (p + ); by the symmetry (2), the other branch of W u (0) connects to p − . This configuration is known as a T-point [31,52], and the locus F 1 is associated with the first, or principal, T-point, denoted T 1 , that is located at (ρ T , σ T ) ≈ (30.8680, 10.1673). There exists an infinite sequence T i of T-points, each of which can be found accurately and systematically with a BVP setup; we refer to [18] for details.
We cannot continue the curve F 1 of first foliation tangency through the point T 1 , but we can find the other side of F 1 past T 1 by setting up the BVP anew. This is illustrated in Figure 8. Panel (a) shows that the locus F 1 has two branches that meet at the point T 1 . From the point T 1 emerge two curves, labeled α 1 and h 1 . Along α 1 , there is a sudden switch in the direction of escape to infinity for the one-dimensional stable manifolds W s (p ± ). This phenomenon, referred to as an α-flip [18], leads to an additional half turn for one branch of W s (p ± ). There are other curves α i for any integer i that end at the respective principal T-points T i ; see [18] for details. The curve h 1 represents the main homoclinic and heteroclinic connections of p ± , which lie extremely close together so that they cannot be distinguished in the (ρ, σ)-plane. Also shown is a circle C η with center T 1 and radius η = 2.8729, which is parameterized as (15) ρ = ρ T + η cos(θ), σ = σ T + η sin(θ).  Figure 8. Computing the locus F1 past the T-point T1. Panel (a) shows the curves F1 (brown), h1 (green), and α1 (teal) in the (ρ, σ)-plane with the circle Cη for η = 2.8729 centered at T1. Panel (b) shows intersection curves of W u (Γ rl ) in Σρ at θ = 0.3142, indicated by the black dot on Cη, together with the direction V at p + ; the range shown is the square with all sides at distance 5 from p + .
The curve F 1 intersects C η at θ F ≈ 2.1231 and θ F ≈ −0.7680, while α 1 intersects C η at θ α ≈ 3.28720 and h 1 at θ h ≈ 0.14097. The black dot on the circle C η labeled (b) corresponds to θ = 0.3142; the situation on Σ ρ for this parameter point is illustrated in Figure 8(b), which shows intersection curves in W u (Γ rl ) near p + , the point q u , the direction V at p + , and the locus C. We observe that the phase portraits in Figure 8(b) and Figure 5(b2), for (ρ, σ)-pairs on either side of h 1 , are topologically different. Nevertheless, the curve W u F (Γ rl ) has the end point q u in Σ ρ also in Figure 8(b). This allows us to set up a BVP as in subsection 4.2 but now with the additional condition that the values of ρ and σ lie on a circle C η . That is, we continue the BVP for u 0 and u Γ with conditions (8), (9), (11)- (13), and (15), where θ is a decreasing (or increasing) continuation parameter, until (6) is again satisfied. We then compute the new branch of F 1 as before with ρ and σ as continuation parameters and without requiring (15). Indeed, the top branch of F 1 in Figure 8(a) was computed in this way. Note that the overall curve F 1 is smooth at T 1 . (ρ, σ)-plane. Figure 9 is a direct comparison of our computed bifurcation diagram in panel (a) with the corresponding sketch by Bykov and Shilnikov [16,17] in panel (b). The main elements of comparison are the curve EtoP of heteroclinic connections from the origin to Γ rl and the curve F 1 of first foliation tangency, with the T-point on it; they are denoted EtoP, F 1 , and T 1 , respectively, in panel (a) and l a , l k , and Q, respectively, in panel (b). The agreement between these objects is very good. Notice, in particular, that the intersection  The agreement of the two panels of Figure 9 over the parameter range of panel (b) confirms simultaneously the validity of our continuation approach for finding F 1 as well as Bykov and Shilnikov's method for finding the equivalent curve l k . Indeed, their approach is quite different and not based on continuation. Rather, to compute the locus l k , Bykov and Shilnikov consider W u (0) at the boundary of the chaotic attractor as a separatrix. Specifically, they consider the separatrix value A, which describes the orientation of the Poincaré return map: for A < 0, a hooked horseshoe exists, and for A > 0, it does not. Therefore, the locus A = 0 are the points at which there is a topological change in the chaotic attractor that corresponds to the loss of the foliation condition. To approximate the value of A, Bykov and Shilnikov determine in [17] a change in the orientation of homoclinic orbits of 0. The manifold W s (0), when followed around a homoclinic loop of W u (0), may form either an orientable surface, homeomorphic to a cylinder, or a nonorientable one, homeomorphic to a Möbius band. The homoclinic orbit switches from being orientable to being nonorientable when A = 0, which occurs at a codimension-two bifurcation point called an inclination flip [40]. For the Lorenz system, there are many different homoclinic loops to 0, and the associated inclination flip points appear to lie on a curve; it is not clear, however, that l k is indeed a curve. Bykov and Shilnikov computed l k by following a trajectory close to a homoclinic loop of W u (0) and checking on which side of W s (0) this trajectory returns to a cross section close to 0 and transverse to W s (0); the cross section used by Bykov and Shilnikov had to be chosen carefully, especially for parameter values near T-points [57]. Further computational details are not given in [17] of how and for which homoclinic orbits the inclination flip points were identified to determine and sketch the curve l k in the (ρ, σ)-plane. Figure 9(a) shows the (ρ, σ)-plane over a larger range; moreover, apart from the curves EtoP and F 1 , we also computed and show the curve h r of the first homoclinic bifurcation, the curve H of Hopf bifurcation, and the curves α 1 and h 1 associated with the T-point T 1 . In this way, we are able to distinguish regions of different dynamics of the Lorenz system. In the white region to the left of the curve h r the equilibria p ± are stable and the only attractors. In the gray-shaded region between h r and EtoP, one finds the preturbulent regime [41], where p ± are still the only attractors but arbitrarily long chaotic transients can be found. To the right of EtoP, one finds a chaotic attractor that coexists with the attractors p ± to the left of the curve H. As was already mentioned, the Lorenz attractor L exists in the purple-shaded region bounded by the curves EtoP and F 1 . As Figure 9(a) shows, we find a second intersection point between these two curves at K 2 ≈ (48.02, 4.02) that was not known to Bykov and Shilnikov. Hence, we conclude that the region of existence of L is, in fact, bounded in the (ρ, σ)-plane. In the blue-shaded regions, a so-called quasi-attractor exists [2,58]. In this region of the parameter space, the quasi-attractor undergoes many more bifurcations, including perioddoubling cascades, which are not shown in the figure; see [11,26,60] for more details.

5.2.
The locus in the full parameter space. An advantage of our BVP continuation approach is that, once we have identified points on either of the two branches of the curve  1 in the (ρ, σ)-plane, we can follow them as curves in any combination of two parameters. Figure 10 shows the locus F 1 of first foliation tangency as part of the bifurcation diagram in the (ρ, β)-plane for fixed σ = 10 in panel (a) and in the (σ, β)-plane for fixed ρ = 28 in panel (b); the labels and shading are as in Figure 9(a). Note that the principal T-point T 1 also appears in each of these parameter planes, at T 1 ≈ (30.4744, 2.6232) in panel (a) and at T 1 ≈ (8.9466, 2.3490) in panel (b), and again divides the curve F 1 into two branches. Moreover, Figure 11. The surfaces EtoP (blue) and F1 (brown) up to their intersection curves K1 and K2 form a cone in (ρ, σ, β)-space, shown here for (ρ, σ, β) ∈ [10, 120] × [0, 50] × [1,7]; also shown is the surface H (red) inside the cone. The curve T1 lies on F1, and it meets K1 and K2 at the T-point-Hopf bifurcation point TH at (ρ, σ, β) ≈ (16.4907, 3.9297, 1.0324). The (ρ, σ)-plane for β = 8 3 from Figure 9(a) is shown in gray with the corresponding intersections of F1, EtoP, and H highlighted as curves.
we find that the curves EtoP and F 1 have two intersection points, at K 1 ≈ (18.86, 1.64) and K 2 ≈ (155.98, 9.03) in the (ρ, β)-plane of panel (a) and at K 1 ≈ (3.23, 1.53) and K 2 ≈ (18.10, 2.70) in the (σ, β)-plane of panel (b). Hence, we conclude that the region of existence of L between the curve F 1 and EtoP is a bounded region in all three parameter planes in Figures 9(a) and 10. Notice further that the general features of the three two-parameter bifurcation diagrams are very similar. Figure 11 shows that, in the full (ρ, σ, β)-space, the region of existence of the Lorenz attractor L is a slanted cone bounded by the surfaces EtoP and F 1 that intersect transversally in curves K 1 and K 2 . Also shown are the surface H of Hopf bifurcation inside the cone and the curve T 1 of T-points on the surface F 1 . The slanted conical shape explains why a bounded region of existence of L can be found in each cross section defined by a parameter pair, provided that the respective third parameter is large enough. Moreover, we conjecture that in these three cross sections, the relative positions of the curves EtoP, F 1 , and H are always topologically as shown in Figures 9(a), 10(a), and 10(b). We illustrate this observation in Figure 11 by showing the corresponding segments of the bifurcation curves in the (ρ, σ)-plane for β = 8 3 . In fact, the surfaces in Figure 11 were obtained by computing the curve segments of EtoP, F 1 , and H in the (ρ, σ)-plane for 60 different values of β ∈ [1,7] and then rendering the corresponding data as individual surfaces; intersection data were used to obtain the curves K 1 and K 2 and the intersection between H and F 1 , while the curve T 1 was computed directly by continuation.
The tip of the slanted cone in Figure 11 is formed by a point TH of T-point-Hopf bifurcation at (ρ, σ, β) ≈ (16.4907, 3.9297, 1.0324). As the name suggests, at the codimension-three bifurcation point TH, the equilibria p ± involved in the T-point heteroclinic connection undergo a Hopf bifurcation and, hence, are nonhyperbolic. In (ρ, σ, β)-space, this corresponds to the curve T 1 ending on the surface H. Past the point TH, the (pair of) connecting orbits from the origin 0 to p ± are no longer of codimension two but rather structurally stable because p ± are now attractors. The T-point-Hopf bifurcation point in the Lorenz system has been found and studied by Algaba et al. in [5]. They show that a curve of codimension-two T-point bifurcation emerges from this codimension-three point and that codimension-one connections from 0 to the bifurcating saddle periodic orbits exist nearby. The latter correspond to the surface EtoP in Figure 11, which indeed meets the surface H exactly at the point TH.
The new observation in Figure 11, which goes well beyond the results in [5], is that the T-point-Hopf bifurcation point TH also involves the locus F 1 of first foliation tangency and, hence, the emergence of a region of existence of the Lorenz attractor L. In fact, we find that the two intersection curves K 1 and K 2 of the surfaces EtoP and F 1 are two branches of a smooth curve K that is transverse to T 1 at TH; the intersection curve of F 1 and H is also smooth and tangent to K at TH. Our results show that the codimension-three point TH emerges as the organizing center that gives rise to the Lorenz attractor L and the cone that bounds its region of existence. It is a natural conjecture from Figure 11 that the region of existence of L is unbounded in (ρ, σ, β)-space. This is corroborated by the work of Barrio and Serrano in [8,11,12], who use numerical integration and the computation of Lyapunov exponents to identify regions of chaotic dynamics via the computation of Lyapunov exponents. More specifically, these authors consider where chaotic dynamics exist in the (σ, β)-plane of the Lorenz system for various values of ρ; they find that the corresponding regions grow as ρ increases, yielding an overall conical shape of the region of existence of chaotic dynamics in the (ρ, σ, β)-space. Barrio and Serrano do not distinguish the region of existence of L within the overall parameter set of chaotic dynamics. Yet for low values of the parameters ρ, σ, and β, the cone of existence of L in Figure 11 appears to agree well with the parameter region of chaotic dynamics found in [12]. (ρ, σ)-plane. The first foliation tangency created by crossing F 1 corresponds to the loss of the foliation condition and provides a boundary to the region of existence of L. However, we find that there are further first foliation tangencies F i associated with further T-points T i . Figure 12 shows a larger region of the (ρ, σ)plane with the additional loci F 2 and F 3 that pass through the T-points T 2 and T 3 , respectively. As discussed earlier, each T-point T i comes with associated curves h i of homo/heteroclinic connections and α i of α-flips. We also show the first homoclinic bifurcation h r , the heteroclinic

Further first foliation tangencies in the
σ 0 40 Figure 12. Bifurcation diagram in the (ρ, σ)-plane with the first three first foliation tangency curves Fn for n = 1, . . . , 3 (brown); they contain the T-points (purple) T1 ≈ (30.8680, 10.1673), T2 ≈ (85.0292, 11.8279), and T3 ≈ (164.1376, 12.9661), respectively, with the associated curves of homo/heteroclinic connections hi (green) and α-flips αi (teal). Also shown are the curves hr, EtoP, and H, with the same colors and shadings as in Figure 9(a), and further homoclinic bifurcation curves h rl , h rlr , h rlrr , and h rlrll (black). The curve PD1 bounds the dark-gray region, where the periodic orbit Γ rl is an attractor. The limit of a cascade of period doublings in the light-gray region is indicated by the curve PD5 of fifth period-doubling bifurcation. The black diamonds labeled (a)-(c) correspond to columns of Figure 13. connection EtoP, and the Hopf bifurcations H from Figure 9(a) with the same shading of the different regions. The periodic orbit Γ rl is an attractor in the dark-gray shaded region, bounded by the curve PD 1 of a period-doubling bifurcation, which is part of a cascade of period-doubling bifurcations; we plot the fifth period-doubling bifurcation curve PD 5 , which indicates the limit of this cascade. Finally, we show four further curves of homoclinic bifurcations, labeled h rl , h rlr , h rlrr , and h rlrll , indicating the respective symbol sequence of the associated homoclinic orbit. Notice that h rlrr and h rlrll spiral into T 2 and T 3 , respectively, while h rl and h rlr pass in between T-points.
A successive crossing of the curves F i as ρ and/or σ increase corresponds to the onset of additional tangencies between the intersection curves W u (Γ rl ) and W s (0) with Σ ρ . How this works geometrically is illustrated in Figure 13 with images of the respective invariant manifolds at the parameter points (black diamonds) labeled (a)-(c) in Figure 12. Specifically, the top panels of Figure 13 show Γ rl with its two-dimensional unstable manifold W u (Γ rl ), the one-dimensional unstable manifold W u (0), and the plane Σ ρ ; the bottom panels show the corresponding intersection curves W u (Γ rl ) and points W u (0) with Σ ρ , together with W s (0). Observe from panels (a2)-(c2) that none of the curves in W u (Γ rl ) crosses the tangency locus C. Figure 13(a) is for (ρ, σ) = (50, 15), a parameter pair in between F 1 and F 2 , and panel (a2) clearly shows that there are tangencies between W u (Γ rl ) and W s (0) in the central region Figure 13. The unstable manifold W u (Γ rl ) (red) and the two foliations in Σρ after F1 for (ρ, σ) = (50,15) in column (a), after F2 for (ρ, σ) = (100, 20) in column (b), and after F3 for (ρ, σ) = (170, 255) in column (c); see the correspondingly labeled points (black diamonds) in Figure 12. Panels (a1)-(c1) show Γ rl (yellow) with W u (Γ rl ) (red), W u (0) (brown), and Σρ (green), and panels (a2)-(c2) show W u (Γ rl ) (red), W s (0) (blue), and the tangency locus C (gray) in Σρ; compare with Figures 3 and 4. (b1) where the flow is down; compare with Figure 4(b). These tangencies are due to the U-shape of curves in W u (Γ rl ); there is a single point of tangency with a curve of the stable foliation on each of the curves in W u (Γ rl ). When F 2 is crossed, the central part of W u (Γ rl ) appears to scroll around the z-axis; see Figure 13(b1). Panel (b2) shows that the curves in W u (Γ rl ) now make a full turn around the origin of Σ ρ ; hence, each curve in W u (Γ rl ) now has two points of tangency with the stable foliation. When F 3 is crossed, one finds an additional half scroll of W u (Γ rl ) around the z-axis in Figure 13(c1), so that curves in W u (Γ rl ) now make one and a half turns around the origin of Σ ρ and each curve has three points of tangency with the stable foliation.
We remark that Dellnitz et al. [19] computed the closure of the unstable manifold W u (p + ) of p + , which is a global attractor, for different values of β ∈ [0.4, 8 3 ] and fixed ρ = 28 and σ = 10 to illustrate their box covering method. As β is decreased, their computed attractor shows an increasing amount of spiraling behavior around the z-axis. With their focus on computation and visualization, the authors of [19] did not further investigate this phenomenon. Their global attractor is reminiscent of W u (Γ rl ), as shown in Figure 13(c1), since decreasing β for fixed ρ = 28 and σ = 10 is analogous to increasing ρ for fixed β = 8 3 and σ = 10; compare Figures 9 and 10.
We detected and continued the first foliation tangencies along the loci F 2 and F 3 with our multi-segment BVP setup from subsection 4.2 as follows. The computation is started at fixed θ = 0.3142 on the circle C η with radius η, shown in Figure 8(a). We continue the BVP with conditions (8), (9), (11)- (13), and (15) and increase η as a continuation parameter until we detect further zeros in (6). Each such zero corresponds to the moment of onset of another foliation tangency. Once it is detected, the new first foliation tangency is continued in ρ and σ as the solution of the BVP with conditions (8), (9), and (11)- (14). In this way, we compute the loci F 2 and F 3 above and up to the T-points T 2 and T 3 , respectively. To detect and compute the other sides of F 2 and F 3 , we continue the BVP with conditions (8), (9), (11)- (13), and (15) and increase η until η = 60 and η = 140, respectively. We then fix η at each of these values and continue the same BVP but now with θ as a decreasing continuation parameter. In this way, a zero of (6) is detected in each case, and continuation in ρ and σ of the BVP with conditions (9) and (11)- (14) yields the other sides of F 2 and F 3 , below T 2 and T 3 . Notice from Figure 12 that the loci F 1 , F 2 , and F 3 end very near the curve PD 5 , that is, near the limit of the period-doubling cascade starting from the periodic orbit Γ rl that is used in our computational setup.
It is interesting to compare the bifurcation diagram in Figure 12 with the figure by Barrio, Shilnikov and Shilnikov from [13, Figure 8(a)] that illustrates the type of attractor of the Lorenz system. Specifically, these authors considered a grid of parameter points in the (ρ, σ)plane and determined numerically "solid-color regions associated with constant values of the kneading invariant [that] correspond to simple dynamics dominated by stable equilibria or stable periodic orbits" [13]. The curves h r and EtoP are readily identified. More importantly, one can observe changes of color between neighboring regions that appear to lie along curves associated with T-point bifurcations T i . These color changes correspond to the addition of a half twist in the template of the attractor, and they align very well with the loci F 2 and F 3 in Figure 12.
6. Conclusions and outlook. We characterized and found the moment of loss of the foliation condition, when the geometric Lorenz map no longer faithfully describes the dynamics of the classic Lorenz attractor L. More specifically, we identified and computed the onset of tangencies between the stable and unstable foliations in the classic Poincaré section Σ ρ through the secondary equilibria p ± . To this end, we computed the two-dimensional unstable manifold W u (Γ rl ) of the symmetric periodic orbit Γ rl and its intersection curves W u (Γ rl ) with Σ ρ , which lie in the unstable foliation. We also computed curves of the stable foliation as the intersection curves W s (0) of the stable manifold W s (0) of the origin 0 with Σ ρ . To find the moment of first foliation tangency, we formulated and implemented a multi-segment boundary value problem to define a regular test function that is zero when the tangent vector at the end point of a curve in W u (Γ rl ) is collinear with the direction of the stable foliation. This allowed us to determine the value of ρ = ρ F ≈ 31.01 as the moment when the foliation condition fails for σ = 10 and β = 8 3 . This value was confirmed as accurate to two decimal places by independent computations with unstable manifolds W u (Γ rll ) of Γ rll and W u (p + ) of p + . Our approximation of ρ F lies outside the interval suggested by Sparrow [60] (from simulations of how vectors return to Σ ρ ) but matches very well the approximation found in the work of Bykov and Shilnikov [17].
An advantage of our approach is that the locus F 1 of the first foliation tangency can be continued readily as a curve in two parameters. There is very good agreement between the curve F 1 in the (ρ, σ)-plane with β = 8 3 fixed and the corresponding curve l k found by Bykov and Shilnikov [17], who identified numerically when flip bifurcations of various homoclinic orbits occur; we view this agreement as a mutual verification of two complementary approaches. Our implementation of simply the condition that the foliations in Σ ρ become tangent does not require making assumptions about the underlying dynamics, such as the existence of homoclinic orbits. We were able to continue the locus F 1 further in the (ρ, σ)plane and show that it intersects the curve of EtoP connections twice; hence, the region of existence of the Lorenz attractor L is, in fact, bounded. We find that it is bounded for all sufficiently large values of β as well as in the (ρ, β)-plane and the (σ, β)-plane for sufficiently large values of σ and ρ, respectively. By computing segments of the respective bifurcation curves in the (ρ, σ)-plane for 60 values of β ∈ [1,7], we show that the region of existence of L has the shape of a slanted cone in the three-dimensional (ρ, σ, β)-parameter space; this is consistent with the findings of Barrio and Serrano [8,11,12] regarding parameter regions with chaotic dynamics.
In each of the different two-parameter planes, we find that the curve of first foliation tangency has two branches on either side of a point of principal T-point bifurcation, where one finds a codimension-two heteroclinic cycle between the origin 0 and the equilibria p ± . Hence, in (ρ, σ, β)-space, there exists a curve of principal T-point bifurcation that lies on the surface of the loss of the foliation condition, and we find that the tip (ρ, σ, β) ≈ (16.4907, 3.9297, 1.0324) of this cone is a codimension-three T-point-Hopf bifurcation point. The bifurcation structure near a T-point-Hopf bifurcation in the Lorenz system has recently been studied by Algaba et al. [5] but is not yet complete. Our results are consistent with their analysis and, furthermore, show that the locus of the loss of the foliation condition and the Lorenz attractor L are effectively created by this codimension-three point. Determining the exact connection between the T-point-Hopf bifurcation and the loss of the foliation condition is an interesting subject for further research.
There exists a sequence of further T-point bifurcations in the Lorenz system that can be found systematically by identifying and continuing the phenomenon we refer to as an α-flip [18]. It is a natural conjecture that a first foliation tangency F i is associated with any of the infinitely many T-points T i . Indeed, we computed the curves F 2 and F 3 of first foliation tangencies in the (ρ, σ)-plane that are associated with the T-points T 2 and T 3 , respectively. Each such tangency between the stable and unstable foliations leads to changes of the observed attractor. Our bifurcation diagram with loci of first foliation tangencies in the (ρ, σ)-plane is in very good agreement with the numerical images by Barrio,Shilnikov,and Shilnikov [13] that illustrate the symbolic type of the attractor. The boundary of the quasiattractor was also found by Shilnikov et al. in the related Shimizu-Morioka system, again with the method from [17] of finding and connecting inclination flip points of homoclinic orbits in a two-parameter plane [56,58]. A sequence of additional T-points with segments of the associated "hook formation curves" are shown in [58, Figure 13]. It would be interesting to compute these loci by continuation as curves of first foliation tangencies and place them in the context of the wider bifurcation diagram of the Shimizu-Morioka system.
A related interesting question is to study further how other global bifurcations of the Lorenz system relate to the overall bifurcation diagram in the (ρ, σ)-plane presented in this paper. In particular, many (structurally stable) heteroclinic connections from p + to 0 were identified in [23] as part of an infinite sequence. All of these emerge from the first homoclinic bifurcation h r ; when continued in ρ, each of these connecting orbits, except for the basic one, has a fold for some high value of ρ and, when continued past the fold for decreasing ρ, ends at a particular secondary homoclinic bifurcation. The many secondary homoclinic bifurcations found in [23] give rise to curves in the (ρ, σ)-plane in the range 0 ≤ ρ ≤ 200 that accumulate on the curve EtoP. In the process, these curves cross the loci F i of first foliation tangency and therefore, the corresponding homoclinic bifurcations are expected to undergo inclination flip bifurcations [17]. Inclination flip bifurcations have been linked to emerging stability windows of stable equilibria and periodic orbits in the Shimizu-Morioka system [67]. These numerous inclination flip bifurcations, of which there exist several types [40], and their role for the organization of nearby chaotic dynamics and stability windows can be studied in the spirit of recent investigations [3,30]. This is an interesting direction for future work.
The BVP approach presented here allows one to find and investigate systematically (first) tangencies between invariant manifolds that lead to topological changes of chaotic attractors, including their creation and destruction. We expect that this will also be relevant for the study of other systems, especially of other three-dimensional dissipative systems. Examples are the Rössler system and the Rosenzweig-MacArthur models considered by Barrio,Blesa,and Serrano [7,9,10]. In particular, they [10, Figure 1] identify a boundary curve in a twoparameter plane that "determines a change in the topological structure of chaotic attractors [of the Rössler system and a tritrophic food chain] from spiral . . . to screw shaped" [10]. This curve passes through a central or focal point of a bifurcation structure that is strikingly similar to that near a T-point of the Lorenz system. The Rössler system [53] has different attractors for different parameter values, which are referred to in [29] as the simple single-scroll Rössler attractor and the more complicated funnel. The appearance of additional extrema in the first-return map corresponds to a topological change in the structure of the attractor via the creation of additional funnel structures [7,10]; a transition from scroll to funnel attractor has also been observed in a Chua circuit experiment in this way [50]. It would be interesting to find and visualize the stable and unstable manifolds in these systems, to determine how their interactions organize the observed topological changes, and to compute loci of first foliation tangencies directly by continuation.
Finally, we mention that it will be an interesting challenge to consider with similar methods bifurcations of Lorenz-like attractors in higher dimensions. Higher-dimensional geometric Lorenz attractors have been constructed, for example, by Bamón, Kiwi, and Rivera-Letelier [6] and Shilnikov and Turaev [62]. The system constructed in [6] is a noninvertible map of the plane that has been shown to contain a so-called wild Lorenz-like chaotic attractor, which is a higher-dimensional version of a Lorenz attractor that contains a hyperbolic set with robust homoclinic tangencies. Bifurcations that generate this type of dynamics in this map were studied recently in [39]; in particular, tangency bifurcations of stable and unstable manifolds of different hyperbolic sets were identified as important ingredients. In the vector-field context, the role of a T-point might be played by a heterodimensional cycle formed by intersections of the stable and unstable manifolds of saddle periodic orbits with different unstable dimensions. This type of structure has been shown to imply the robustness of homoclinic tangencies in higher dimensions [15]. The minimal example of a heteroclinic cycle between two saddle periodic orbits in a four-dimensional vector field has been identified recently in a model of intracellular calcium dynamics [68]. The investigation of the bifurcation structure of invariant manifolds in this vector-field model is ongoing. Table 1 Accuracy parameters used in the computations for the given figures. foliation, and the direction V , which represents the tangent at q u to the corresponding leaf in the stable foliation. Leaves of the unstable foliation can be approximated by using the unstable manifolds of periodic orbits or of p ± . We computed ρ F with three different unstable manifolds, namely, W u (Γ rl ), W u (Γ rll ), and W u (p + ), and increasingly smaller values of κ. The results are shown in Table 2. The computations were performed with other accuracy parameters as listed in Table 1. We approximated ρ F for κ = 10 −k , where k = 1, . . . , 5, except for W u (Γ rll ), for which we only got κ down to 10 −4 ; as mentioned in subsection 4.2, smaller values of κ are impractical due to the associated very close passage near the origin of the orbit segment u Γ . Table 2 suggests that the relative error appears to decrease quadratically with κ for each of the three manifolds. This is in line with the theoretical expectation because we are approximating a tangent of a generically quadratic curve with a secant through points at distance κ. We checked that an increase in the number of mesh points along each orbit segment affects the values of ρ F in Table 2 only in the fourth decimal place. Overall, Table 2 provides the evidence that estimates for ρ F are accurate to two decimal places when using κ = 10 −4 .
As discussed in section 4, we use the line V , given as the direction of W s (p + ) projected onto Σ ρ , as a convenient local approximation of the tangent to the stable foliation W s (0) at the point q u near p + . To check the influence of the choice of V on the accuracy of the computed value of ρ F , we proceed in the following way. For fixed ρ = 31.0, we compute the curves in W s (0) up to total integration time τ max = −5 and consider the curve nearest p + in Table 2 Estimates of ρ F computed with each of the unstable manifolds W u (Γ rl ), W u (Γ rll ), and W u (p + ) and various choices of κ. κ ρ F using W u (Γ rl ) ρ F using W u (Γ rll ) ρ F using W u (p + ) Since the stable direction changes linearly along W u F (Γ rl ), to first approximation, we interpolate linearly between V F and V to find the improved vector V int = −0.85751589 0.34066195 at q u ≈ (8.8294691, 8.75783413). Instead of V , we then use V int , as well as the arguably less accurate vector V F in our computation of ρ F , as described in subsection 4.1. The first row in Table 3 summarizes the estimates of ρ F for each of these vectors. The next two rows are the angle difference from V in both radians and degrees. We find that both V int and V F lead to the same approximation of ρ F to two decimal places in spite of an angle difference at ρ = 31 of up to about 2.26 degrees. This suggests that our BVP setup is quite insensitive to the exact choice of direction vector for the stable foliation at q u .