Kneadings, Symbolic Dynamics and Painting Lorenz Chaos. a Tutorial
Abstract
A new computational technique based on the symbolic description utilizing kneading invariants is proposed and verified for explorations of dynamical and parametric chaos in a few exemplary systems with the Lorenz attractor. The technique allows for uncovering the stunning complexity and universality of biparametric structures and detect their organizing centers – codimensiontwo Tpoints and separating saddles in the kneadingbased scans of the iconic Lorenz equation from hydrodynamics, a normal model from mathematics, and a laser model from nonlinear optics.
Received (to be inserted by publisher)
Keywords: kneading invariant, symbolic dynamics, Tpoints, Lorenz attractor, chaos, homoclinic and heteroclinic orbits
1 Introduction
A great deal of analytical and experimental studies, including modeling simulations, have been focused on the identification of key signatures to serve as structural invariants that would allow dynamically alike nonlinear systems with chaotic dynamics from diverse origins to be united into a single class. Among these key structures are various homoclinic and heteroclinic bifurcations of low codimensions that lie at the heart of the understanding of complex behaviors because of their roles of organizing centers of dynamics in parameterized dynamical systems.
Dynamical systems theory has aimed to create purely abstract approaches that are further proceeded by creation and the development of applicable tools designed for the search and identification of such basic invariants for simple, MorseSmale, systems and ones with complex chaotic dynamics. One such a [computationally justified] approach for studying complex dynamics capitalizes on the concept of sensitivity of deterministic chaos. Sensitivity of chaotic trajectories can be quantified in terms of the divergence rate evaluated through the largest Lyapunov characteristic exponent. The approach has been proven to work exceptionally well for various systems with chaotic and simple dynamics. In several loworder dissipative systems, like the Rössler model, the computational technique based on the largest Lyapunov characteristic exponent reveals that they possess common, easily recognizable patterns involving spiral structures in biparametric planes [Afraimovich et al., 1977; Bykov, 1993; Shilnikov, 1991; Shilnikov et al., 1993; Barrio et al., 2011b; Gallas, 2010]. Such patterns turn out to be ubiquitously alike in both timediscrete [Lorenz, 2008] and timecontinuous systems [Gaspard et al., 1984; Barrio et al., 2009; Gallas, 2010], and they are easily located when the spiral patterns have regular and chaotic spiral “arms” in the systems with the Shilnikov saddlefocus [Shilnikov & Shilnikov, 2007].
Application of the Lyapunov exponents technique fails, in general, to reveal fine structures embedded in the biparametric scans of Lorenzlike systems. As such it cannot deliver the desired insights into intrinsic bifurcations because regions of chaotic dynamics appear to be uniformly. This basically means that the instability of the Lorenz attractors does not vary noticeably as control parameters of the system are varied. This holds true too when one attempts to find the presence of characteristic spiral structures that are known to exist theoretically in the Lorenzlike systems [Bykov, 1993; Glendinning & Sparrow, 1986] and therefore could only be identified using accurate bifurcation continuation approaches [Shilnikov, 1991; Shilnikov et al., 1993]. Such spirals in a biparametric parameter plane of the system in question are organized around the socalled T[erminal]points corresponding to codimensiontwo or higher heteroclinic connections between two or more saddle equilibria. For symmetric systems with the Lorenz attractor, the degeneracy is reduced to a cod2 bifurcation of a closed heteroclinic connection involving two saddlefoci and a saddle at the origin. Tpoints have been located in various models of diverse origins including electronic oscillators [Bykov, 1998; FernándezSánchez et al., 2002] and nonlinear optics [Forysiak et al., 1991; Wieczorek & Krauskopf, 2005], etc.
Figure 1 sketches an idea of the structure of the bifurcation unfolding near an ordinary Tpoint in the parameter plane in a symmetric system [Bykov, 1993; Glendinning & Sparrow, 1986]; various Tpoints configurations for other heteroclinic connections were examined in detail in [Bykov, 1993]. Here, the heteroclinic connection is formed by a pair of symmetric saddlefoci and a saddle whose onedimensional stable (incoming) and unstable (outgoing), respectively, separatrices merge only at the codimensiontwo Tpoint in the biparametric space of the Lorenz model. It follows from the theoretical analysis that the unfolding of this bifurcation includes three main curves ending at the Tpoint: a spiraling bifurcation curve, (1hom ) corresponding to two simultaneous (due to the symmetry) homoclinic loops of the saddle, a curve, (1hom ), corresponding to two simultaneous homoclinic loops of the saddlefoci, and another codimensionone bifurcation curve, (1het ), corresponding to a heteroclinic correction between both saddlefoci. In addition, the unfolding includes infinitely many subsidiary Tpoints in between 1hom and 1het .
Despite a rather overwhelming number of studies reporting the occurrence of various spiral structures, there is yet unproportionately little known about construction details and generality of underlying bifurcation scenarios giving rise to such patterns. In this paper we introduce a novel computational toolkit capitalizing on the idea of the symbolic representation for the dynamics of Lorenzlike systems that employs the kneading invariants [Milnor & Thurston, 1988]. We will then show how the toolkit can be used for detecting various structures in biparametric scans of such systems. It is our intention to enhance further the technique thus allowing for systematic studies of the stunning complexity and universality of spiral structures in models of diverse dynamics and origins.
The paper is organized as follows: in Section 2 we review the homoclinic bifurcation theory and discuss the routes to chaos and the formation stages of the strange attractor in the Lorenz model; in Section 3 we introduce basics of kneading theory; in Section 4 we demonstrate of the computational technique using kneading invariants to reveal hidden structures of biparametric chaos in the iconic Lorenz model; in Section 5 we apply the technique to uncover colorfulness of the bifurcation diagrams of the ShimizuMorioka model. Finally, In Section 6 the proposed technique will be tested on a 6D model of the optically pumped, far infrared red threelevel laser [Moloney et al., 1989; Forysiak et al., 1991] to confirm the universality of the patters produced by the deterministic chaos in the Lorenz like systems.
2 Homoclinic bifurcations in systems with the Lorenz attractor
The strange chaotic attractor in the Lorenz equation from hydrodynamics has become a defacto proof of deterministic chaos. The butterflyshaped image of the iconic Lorenz attractor, shown in Fig. 5, has become the trademark of Chaos theory and Dynamical Systems. This theory elaborates on complex trajectory behaviors in nonlinear systems from mathematics, physics, life sciences, finance, etc. Universality of the methods along with bifurcation tools has made them spread wide and deep across all other disciplines of the modern science.
The Lorenz equation [Lorenz, 1963] is a system of three differential equations:
(1) 
with three positive bifurcation parameters: being the Prandtl number quantifying the viscosity of the fluid, being a positive constant of magnitude order 1 which originates from the nonlinearity of the Boussinesq equation, and being a Reynolds number that characterizes the fluid dynamics. Notice that Eqs. (2) are symmetric, i.e. [Sparrow, 1982].
2.1 Uniparametric cut through the Lorenz equation
An exploration of primary bifurcations in the Lorenz equation begins with a singleparameter examination of the dynamics, where serves as the bifurcation parameter increasing from laminar to weekly turbulent magnitudes around 25, while the two other parameters are set to the original Saltzman values: (for the water, and for the air) and . This would give a uniparametric cut through the Lorenz equation that was originally explored in these independent studies [Afraimovich et al., 1977; Kaplan & Yorke, 1979] (see Fig. 2):

For the only equilibrium state at the origin is a global attractor in the 3D phase space of the Lorenz equation.

This equilibrium state undergoes a pitchfork bifurcation at , and for becomes a saddle so that the stability is transferred to two symmetric stable foci.

At the unstable separatrices of the saddle return to the origin thus forming a homoclinic butterfly. This causes a “ homoclinic explosion” in the phase space of the model that becomes filled in at once with countably many saddle periodic orbits that would further compose the skeleton of the Lorenz attractor.

For the model exhibits transient chaos with two ultimate attractors: stable foci . Such transient chaos is associated with a preturbulence regime.

The value corresponds to the onset of the Lorenz attractor coexisting with stable foci and . The attraction basins of are bounded by the 2D cylindricallyshaped stable manifolds of two saddle periodic orbits that have earlier bifurcated from the homoclinic loops of the saddle at .

As increases to the saddle periodic orbits shrink the attraction basins and collapse onto the stable foci through a subcritical AndronovHopf bifurcation.

For the Lorenz equation possesses a genuinely strange chaotic attractor, known as the Lorenz attractor, containing no stable orbits.
2.2 Canonical 2D bifurcation diagram of the Lorenz equation
The pilot study of the dynamics of the Lorenz equation needs to be further enhanced by the biparametric examination of the model, including at large parameter values [Barrio & Serrano, 2007, 2009]. We will start off with the canonical bifurcation diagram, shown in the left panel of Fig. 3 (courtesy of [Shilnikov, 1980]), of the Lorenz equation that depicts the basic bifurcation curves in the parameter plane with fixed . The right panel of Fig. 3 sketches en route fragments of the formation of the Lorenz attractor on the pathway, [Afraimovich et al., 1977; Kaplan & Yorke, 1979] through the bifurcation curves. For , Eq. (2) has a single stable equilibrium state at the origin. This equilibrium state undergoes a pitchfork bifurcation at , so that for the origin becomes a saddle, , of the topological type (2,1) due to the characteristic exponents . The saddle has a 1D unstable manifold, is made of itself and a pair of 1D unstable separatrices, and (due to ) entering the saddle as , and a 2D stable manifold, , containing the leading (due to ) invariant axis; the eigenvector due to determines the nonleading or strongly stable direction, in . After the pitchfork bifurcation, the separatrices of the saddle tend to two symmetric attractors – equilibrium states, (Fig. 3(A)) that become the global attractor for all trajectories in the phase space of the Lorenz equation other than in .
A homoclinic butterfly bifurcation occurs in the Lorenz equation when both separatrices, and , of the saddle come back to the origin along the axis (Fig. 3(B)). In virtue of the symmetry of the Lorenz equation, such homoclinic loops are always formed in pairs, and therefore constitute bifurcations of codimensionone, in general. The bifurcation of the homoclinic butterfly takes place on the curve in the parameter plane. Bifurcations of the separatrices of the saddle at the origin are crucial for the Lorenz attractor.
The very first homoclinic butterfly made of two separatrices looping a single round about the equilibrium states , causes the homoclinic explosion in the phase space of the Lorenz equation. This bifurcation gives rise to a onset of countably many saddle periodic orbits that forming an unstable chaotic (saddle) set , which is not an attractor yet. For this explosion to happen, the socalled saddle value , which is the sum of the leading characteristic exponents of the saddle, must be positive; alternatively, the saddle index . Otherwise, if () the homoclinic butterfly produces a symmetric figure8 periodic orbit in the aftermath of the gluing bifurcation through which two stable periodic orbits merge after flowing into the separatrix loops. L. Shilnikov [Shilnikov, 1968] pointed out two more conditions, in addition to the primary one (1): , or , needed for the separatrix loop of the saddle in and higher dimensions to produce a single saddle periodic orbit; they are: (2) , i.e. the separatrix comes back to the saddle along the leading direction, (3) the socalled separatrix value ( in the mapping Eq. (3) below) does not vanish, its sign determines whether the separatrix loop is oriented or twisted, and hence the stable manifolds of the saddle periodic orbit are homeomorphic to a cylinder or a Möbius band. Otherwise, he predicted [Shilnikov, 1981] that the Lorenz attractor could be born right near such codimension2 bifurcations, termed resonant saddle, orbit and inclinationflip, correspondingly [Shilnikov et al., 1998,2001], as it occurs in the ShimizuMorioka and similar models [Shilnikov, 1986; Robinson, 1989; Rychlic, 1990] below.
Out of many saddle periodic orbits, which exploded the phase space of the Lorenz equation, two are the special, , ones as they demarcate the thresholds of the “interior” of the chaotic unstable set (Fig. 3(C)).
After the homoclinic butterfly bifurcation, in the region between the bifurcation curves and in the ()parameter plane, the separatrices and of the saddle switch the targets: now the right/left separatrix tends the left/right stable focus .
In order for the unstable chaotic set to become the Lorenz attractor, it must become invariant, i.e. a closed set containing all limit orbits, and hence no loose of trajectories escaping to stable foci . This occurs on the bifurcation curve, , in the parameter space (Fig. 3(D)). To the right of the curve, the basin of the Lorenz attractor is shielded away from those of the stable equilibrium states by the 2D cylindershaped stable manifolds of the two “threshold” saddle orbits, that have simultaneously emerged from both separatrix loops, at the homoclinic explosion on the curve in the parameter plane.
As one moves further to the right in the parameter plane, the saddle orbits, , keep narrowing the attraction basins of the foci , and on the bifurcation curve they collapse into the stable equilibria. The equilibrium states become saddlefoci of the (1,2)type through a subcritical AndronovHopf bifurcation [Roshchin, 1978]. The topological (1,2)type means that each saddlefocus has 2D unstable and 1D stable manifolds; the latter is formed by two incoming separatrices. Some local properties of the saddlefoci can be revealed without evaluating their characteristic exponents explicitly. Let stand for the real stable exponent of , and stand for a complex conjugate pair such that . Observe that the divergence of the vector field defined by Eqs. (2) is given by , which equals . This implies , i.e., the complex conjugate pair is closer to the imaginary axis than the real negative exponent, and hence the saddlefoci meet the Shilnikov condition [Shilnikov, 1965, 1967; Shilnikov & Shilnikov, 2007]. Therefore, as soon as the saddlefocus possesses a homoclinic loop, such a bifurcation causes the abundance of periodic orbits nearby. Those periodic orbits constantly undergo saddlenode and period doubling bifurcations as the parameters are varied. Moreover, since the divergence of the vector field of the Lorenz equation is always negative, saddlenode bifurcations give rise to stable periodic orbits near the homoclinic saddlefocus bifurcation. Under fulfillment of some global conditions, a single Shilnikov saddlefocus bifurcation can lead to the formation of a spiral or screwlike attractor. However, a strange attractor due to the Shilnikov saddlefocus in a 3D system with a negative divergence is no genuinely chaotic set in the sense that it contains stable periodic orbits within. Because of that, such a chaotic attractor is called quasiattractor thus referring to that because besides stable periodic orbits with weak basins it may have structurally unstable or nontransverse homoclinic orbits [Afraimovich & Shilnikov, 1983; Shilnikov, 1994, 1997]. Note that systems in higher dimensions can possess genuinely strange attractors with the Shilnikov loop without stable periodic orbits, the socalled wild chaotic attractors [Shilnikov, 1994, 1997; Turaev & Shilnikov, 1998; Shilnikov, 2002].
The Lorenz attractor is nonhyperbolic because it includes the singularity at the origin – the saddle equilibrium state with an 1D unstable manifold while all other saddle periodic orbits on the attractor have 2D stable and unstable manifolds. Moreover, the manifolds of those orbits in the Lorenz attractor (self) cross transversally in the 3D phase space thereby producing only structurally stable (transverse) homoclinic and heteroclinic trajectories. This condition imperative for the Lorenz attractor will not always hold for larger parameter values and lead to homoclinic tangencies of the manifolds which are followed by saddlenode bifurcations in the Newhouse regions [Gonchenko et al., 1993, 1996] of the parameter plane of the system. Thus, in order for the Lorenz attractor to be strange and chaotic with no stable orbits, it must not include the [homoclinic] saddlefoci, , as well as contain only structurally stable homoclinic orbits due to transverse intersections of the manifolds of saddle periodic orbits.
3 Kneading invariants
Chaos can be quantified by several means. One customary way is through the evaluation of the topological entropy. The greater the value of topological entropy, the more developed and unpredictable the chaotic dynamics become. Another practical approach for measuring chaos in simulations capitalizes on evaluations of the largest (positive) Lyapunov exponent of a long yet finitetime transient on the chaotic attractor.
After the stable foci have lost the stability through the subcritical AndronovHopf bifurcation, the Lorenz equation exhibits the strange attractor of the iconic butterfly shape. The wings of the butterfly are marked with two symmetric eyes containing the saddlefoci isolated from the trajectories of the Lorenz attractor. This attractor is structurally unstable [Guckenheimer & Williams, 1979; Afraimovich et al., 1983] as it undergoes bifurcations constantly as the parameters of the Lorenz equation are varied. The primary cause of structural and dynamical instability of chaos in the Lorenz equation and similar models is the singularity at the origin – the saddle with two onedimensional outgoing separatrices. Both separatrices fill in densely two spatially symmetric, [], wings of the Lorenz attractor in the 3D phase space (see Fig. 5). The Lorenz attractor undergoes a homoclinic bifurcation when the separatrices of the saddle change the flipflop pattern of switching between the butterfly wings centered around the saddlefoci. At such a change, the separatrices comes back to the saddle thereby causing additional homoclinic explosions in phase space [Afraimovich et al., 1977; Kaplan & Yorke, 1979].
The time progression of the “right” (or symmetrical “left”) separatrix of the origin can be described geometrically and categorized in terms of the number of flipflops around the equilibrium states and in the 3D phase space of the Lorenz equation (Fig. 5). Or, alternatively, can be reduced to the timeevolution of the coordinate of the separatrix, as shown in panel B of Fig. 5. The signalternation of the coordinate suggests the introduction of a based alphabet to be employed for the symbolic description of the separatrix. Namely, whenever the right separatrix turns around or , we write down or , respectively. For example, the time series shown in panel B generates the following kneading sequence starting with etc.
In what follows we will introduce and demonstrate a new computational toolkit for the analysis of chaos in the Lorenzlike models. The toolkit is inspired by the idea of kneading invariants introduced in [Milnor & Thurston, 1988]. A kneading invariant is a quantity that is intended to uniquely describe the complex dynamics of the system that admit a symbolic description using two symbols, here and . The kneading invariant is supposed to depend monotonically on the governing parameter so that any two systems can be compared and differentiated, or equivalently, ordered in terms of and , by the kneadings. Two systems with the Lorenz attractors are topologically conjugate when they share the same kneading invariant, see [Guckenheimer & Williams, 1979; Rand, 1978; Malkin, 1991; Glendinning & Sparrow, 1993; Tresser & Williams, 1993] and the references therein. Moreover, for the Lorenzlike systems, the topological entropy can be evaluated though the kneading invariants by reducing consideration to piecewise monotone mappings of a closed interval [Glendinning & Hall, 1996; Li & Malkin, 2003].
Such mappings are closely related to the geometric models of the Lorenz attractor which are 1D and 2D Poincaré return mappings defined on some crosssection transverse to trajectories of the Lorenz attractor. The basic idea behind either geometric model capitalizes on extended properties of the local Poincaré mapping near a homoclinic butterfly of the saddle [Guckenheimer & Williams, 1979; Afraimovich et al., 1983; Shilnikov et al., 1998,2001]. The mapping is assumed to meet a few other conditions related to the global properties of the flow far from the homoclinic butterfly. Such a 1D constrained mapping shown in Fig. 6(A) can be written as:
(2) 
here is the saddle index, controls the distance between the returning separatrices, , and the saddle, and is a scalar whose sign determines whether the homoclinic loops at are oriented if , or twisted when , see [Shilnikov et al., 1998,2001] for more details.
Loosely speaking, this geometric model (3) should have no critical point on both branches, and moreover possess a property of strong stretching with a rate of expansion more than [Afraimovich et al., 1983]. This would guarantee that the Lorenz attractor will be densely filled in by the forward iterates of the separatrices, with no holes – lacunas containing isolated periodic orbits inside, stable or not. Strong stretching would also imply a monotone dependence of the kneading invariants on a governing parameter.
In a symmetric system with the Lorenz attractor, the kneading invariant is assigned to quantify the symbolic description of either separatrix; in the asymmetric case one should consider two kneading invariants for both separatrices of the saddle. Thus, in respect it reflects quantitatively a qualitative change in the separatrix behavior, such as flipflopping patterns, as the parameter of the system is changed.
By construction, kneading invariants are proposed to serve as moduli of the topological equivalence that are employed to compare or contrast between any two Lorenz attractors or, equivalently, any two Lorenzlike systems. Due to the symmetry of the Lorenz mapping from Eq. (3), forward iterates of the right separatrix, , of the discontinuity point (resp. the saddle) are detected to generate a kneading sequence defined by the following rule:
(3) 
here is the th iterate of the right separatrix of the origin. The condition is interpreted as a homoclinic loop, i.e. the separatrix returns to the origin after steps.
The kneading invariant for the separatrix is defined in the form of a formal power series:
(4) 
Setting make the series (3) convergent. The smallest zero, , if any, of the graph of (3) in the interval yields the topological entropy, , of the 1D mapping (3).
Let us next draw a parallelism between the geometric model and the Lorenz equation: the kneading sequence comprised of only s of the mapping (3) corresponds to the right separatrix converging to the stable equilibrium state, (or possibly, a periodic orbit with ). The corresponding kneading invariant is maximized at . When the right separatrix converges to an limit set with , like the left stable focus, then the kneading invariant is given by because the first entry in the kneading sequence is followed by infinite s. Thus, yield the range of the kneading invariant values; for instance . Two samples of the separatrix pathways shown in Fig. 7 generating the following kneading invariants
illustrate the parallelism.
To conclude this section we point out that there is another approach for constructing 1D return mappings through the evolution of the variable of the separatrix of the Lorenz equation. The mapping generated by the turning points where on the attractor has a distinct cuspshaped form depicted in Fig. 6(B). The point corresponding to the cusp is used for the initiation of the kneading sequence. Due to its unimodal graph with two, increasing and decreasing, segments the corresponding formal power series is then defined as follows:
(5) 
i.e. .
4 Kneading scanning of the Lorenz equation
In this section we carry the concept of the kneading invariants over to numerical studies of fine structures of chaos in the Lorenz equation. For sake of simplicity in this pilot phase, we employ a rough idea for defining the kneading sequences of s and s that relies on whether the right separatrix of the saddle makes a revolution around the right equilibrium state, , or the left one, , respectively, in the ()projection of the 3D phase space. One can utilize even a simpler approach where the sign of the current kneading entry in the sequence is determined by the sign of the coordinate of the separatrix at the offlying turning points, (on the butterfly wings), given by , see the trace in Fig. 5(B). There are other alternative ways for defining kneading entries, involving crosssections, etc, that are not free of certain limitations either. In the future we plan to enhance the current kneading algorithms by utilizing the incoming 1D separatrices of the saddlefoci and finding the winding numbers around them instead.
In this computational study of the Lorenz equation and two other models below, we consider partial kneading power series truncated after the first 50 entries: . The choice of the number of entries is not motivated by numerical precision, but by simplicity, as well as by a resolution of the bitmap mappings for the biparametric scans of the models. One has also to figure the proper value of : setting it too small makes the convergence fast so that the tail of the series has a little significance and hence does not differentiate the fine dynamics of the Lorenz equation on longer time scales.
At the first stage of the routine, we perform a biparametric scan of the Lorenz equation within a specific range in the plane. The resolution of scans is set by using mesh grids of equallydistanced points. Next by accurately integrating the separatrix using Taylor series software TIDES^{1}^{1}1freeware TIDES is a versatile numerical ODE solver for integration of ODEs with an arbitrary precision, especially for chaotic systems. [Abad et al., 2011a, b; Barrio et al., 2011c] we identify and record the sequences for each point of the grid in the parameter plane. Then we define the biparametric mapping: with some appropriately chosen , the value that determines the depth of the scan. The mapping is then colorized in Matlab by using various builtin nonlinear spectra ranging between to and , respectively. In the mapping, a particular color in the spectrum is associated with a persistent value of the kneading invariant on a level curve. Such level curves densely foliate the biparametric scans.
Figure 8 represents the kneadingbased color scan of the dynamics of the Lorenz equation mapped onto a fragment of the ()parameter plane. In the scan, a window of a solid color corresponds to a constant kneading invariant. In such windows the dynamics of the Lorenz equation are dominated by simple attractors such as stable equilibria or stable periodic orbits to which the separatrix converge. A quick examination of the kneading definition (3) reveals that the kneading invariant does not vary at a supercritical AndronovHopf bifurcation and a pitchfork bifurcation describing continuous transitioning between stable symmetric and asymmetric periodic orbits. This holds true for a perioddoubling bifurcation too, because the kneading sequence, say , inherits the same block repeated twice in the code for the periodic orbit of a double period and so forth. While the kneading technique does not detect such safe bifurcation boundaries [Shilnikov et al., 1998,2001] having crossed through which the phase point does not run far from an old attractor to a new one, though it detects dangerous boundaries well, including a generic saddlenode bifurcation, homoclinic bifurcations, and some other.
A borderline between two solidcolor regions corresponds to a bifurcation through a dangerous boundary which is associated with a jump between the values of the kneading invariant. For example, the borderline in Figure 8 between the brown region with the kneading sequence and the blue region, with the kneading sequence , corresponds to the primary homoclinic butterfly of the saddle. The second borderline of the blue region corresponds to the onset of the Lorenz attractor existing on the right from it. One can see that above this border is adjoined by windows of solid colors thus indicating that the separatrices start converging to stable equilibria after chaotic transients, long or short. Indeed, crossing the curve, (or in Fig. 3), above does not imply that the Lorenz equation will have a strange attractor, but a chaotic saddle set [Bykov & Shilnikov, 1989, 1992].
What the proposed kneading technique does extraordinary well, compared to the biparametric screening based on the finitetime Lyapunov exponent approach (see below), and what we developed it for, is the detection of bifurcations within the Lorenz strange attractor. The corresponding yellowreddish regions in the parameter plane in Fig. 8 clearly demonstrate the evidence of the parametric chaos that, like in turbulence, enriched by vortices of multiple Tpoints. Panel B of this figure depicts the kneading mapping near the leftbottom corner of the bifurcation diagram in panel A. In it, the black (blue in panel A) region corresponds to the chaotic saddle dynamics transitioning into the Lorenz attractor after the mapping color shifts to the yellowreddish spectrum. Blue curves in panel B are the bifurcation curves of some homoclinic orbits with short admissible kneadings. One can see from this panel that the mapping spectrum is clearly foliated by the kneading invariant level curves of the colors gradually progressing from red to yellow. This indicates that the new born Lorenz attractor, while being structurally unstable and sensitive to parameter variations, persists initially the pseudohyperbolic property because the foliation remains uniform, and transverse to the classical pathway (Fig. 3). The homogeneous foliation starts breaking around a saddle point after which one singled bifurcation curve spirals onto the primary Tpoint. Far from this point, the curve corresponds to the formation of the homoclinic loop with the kneading ; i.e. the right separatrix makes one excursion around the saddlefocus , followed by three revolutions around the saddlefocus , and then returns to the saddle at the origin. While moving along the spiraling curve toward the Tpoints, the separatrix makes progressively more turns around , or more precisely around the 1D incoming separatrix of this saddlefocus. With each incremental turn around , the separatrix comes closer to while the bifurcation curve becomes one scroll closer to the Tpoint simultaneously. Due to this feature the Tpoint is called a Terminal point. The Tpoint corresponds to the following symbolic sequence . In virtue of the symmetry of the Lorenz equation, the Tpoint actually corresponds to a closed heteroclinic connection involving all three saddleequilibria, see Figs. 1 and 5. The merger of the right (left) separatrix of the saddle with the incoming separatrix of the saddlefocus , increases the codimension (degeneracy) of this heteroclinic bifurcation to two; note that intersections of the 2D unstable manifolds of the saddlefoci, with the 2D stable manifold of the saddle at the origin are transverse in the 3D phase space in general. Breaking the 1D heteroclinic connection gives rise to a primary homoclinic orbit to the saddlefocus, as well as to a heteroclinic connection between both saddlefoci (see the sketch of the bifurcation unfolding of the Tpoint in Fig. 1). The corresponding bifurcation curves of these homoclinic and heteroclinic bifurcations originating rightward from the Tpoint bound a sector containing subsidiary Tpoints [Bykov, 1993; Glendinning & Sparrow, 1984]. Each new Tpoint produces other selfsimilar structures scaled like fractals. Panel C shows two such identified Tpoints: primary and secondary . The primary codimensiontwo Tpoint in the Lorenz equation was originally discovered by Yudovich [Pertovksaya & Yudovich, 1980].
As soon as the saddlefoci and their bifurcations become involved in the dynamics of the Lorenz equation near the primary Tpoint, the Lorenz attractor loses the purity of the genuine chaotic attractor that used to have neither stable periodic orbits nor nontransverse homoclinic trajectories. It transforms into a quasichaotic attractor with weakly stable orbits of small basins and nontransverse homoclinic orbits. The idea of nontransversality or tangency was employed in [Bykov & Shilnikov, 1989, 1992] to numerically identify the second boundary, , in addition to the first one, ( in Fig. 3), that bounds the existence region of the Lorenz attractor in the parameter plane, see Fig. 9. Note that crosses the initial boundary, . This means that above the intersection point, crossing ( in Fig. 3) rightwards does not guarantee that the basin of the Lorenz attractor will necessarily be isolated from the basins of the stable foci, by the stable manifolds of the saddle periodic orbits, (Fig. 4). This implies that the separatrices of the saddle will demonstrate chaotic transient behavior, long or short, prior to them converging to , see [Shilnikov, 1986, 1989; Bykov & Shilnikov, 1989; Shilnikov, 1991; Bykov & Shilnikov, 1992; Shilnikov, 1993; Shilnikov et al., 1993] for full details.
The other feature of the boundary, , is that it passes through the primary Tpoint, thereby separating the existence region of the Lorenz attractor from bifurcations of the saddlefoci, and consequently from all subsidiary Tpoints existing on the right from it in the parameter plane. Here the chaotic dynamics of the Lorenz equation become even wilder and less predictable [Shilnikov, 1994, 1997; Turaev & Shilnikov, 1998; Shilnikov, 2002]. The indication to that is panel A of Fig. 8 that reveals, through the kneading scan, a parametric turbulence in the parameter plane with fractal explosions in the forms of multiple spiral structures – “tornado eyes” centered around Tpoints. Note that basins of spiraling Tpoints are separated by corresponding saddles. One can spot selfsimilar smallerscale spiral structures within largescale ones and so forth. The richness of such fractal structures in the parameter plane results from the synergy of the Lorenzlike dynamics amplified by chaos induced by the Shilnikov saddlefoci.
To conclude this section we contrast the scans of the Lorenz equation obtained using the proposed kneading technique with the sweeps based on the evaluation of the largest Lyapunov exponent, , for the separatrices of the saddle evaluated over a finite time interval [Barrio & Serrano, 2007, 2009; Barrio et al., 2011a]. Figure 10 shows a fragment of the typical biparametric sweep of the Lorenz equation: the dark region at small values of the Reynolds number, , is where is negative on the stable foci. The sweep yields a clear borderline between the regions of the simple and chaotic attractors. The chaotic region (yellowreddish) is characterized by a small positive Lyapunov exponent, variations of which are not significant enough to reveal fine structures, as spiraling Tpoints. The method can detect well stability islands in the biparametric diagram which correspond to emergent stable periodic orbits.
Two panels in Fig. 11 represent the expanded scans of the of Lorenz equation: panel (A) is the kneading invariant mapping on the grid of points and panel (B) shows based sweeping of the parameter plane. While panel A reveals a chain of largescale spiral hubs around Tpoints, panel B reveals none but stability windows (dark). We would like to point out that the stability windows can also be detected in the kneading scan in panel A showing the border of the solid color (dark yellow) island stretched horizontally at small values of the parameter .
5 The ShimizuMorioka model
Next we perform the kneadingbased biparametric scanning of another classical threedimensional system called the ShimizuMorioka model [Shimizu & Morioka, 1980; Shilnikov, 1986, 1989, 1991, 1993]:
(6) 
here, and are positive bifurcation parameters. Like the Lorenz equation, this symmetric model has 3 equilibrium states: a saddle of the (2,1)topological type at the origin, and two symmetric stablefoci or saddlefoci of the (1,2)topological type.
This model was originally introduced to examine a pitchfork bifurcation of the stable figure8 periodic orbit that gives rise to multiple cascades of period doubling bifurcations in the Lorenz equation at large values of the Reynolds number . It was proved in [Shilnikov et al., 1993] that the Eqs. (6) are a universal normal form for several codimensionthree bifurcations of equilibria and periodic orbits on central manifolds. The model turned out to be very rich dynamically: it exhibits various interesting global bifurcations [Shilnikov, 1986, 1991, 1993] including Tpoints for heteroclinic connections.
While the model inherits all basic properties of the Lorenz equation, in addition, and of special interest, are two homoclinic bifurcations of codimensiontwo: resonant saddle with the zero saddle value or the saddle index , and the orbitflip bifurcation where the separatrix value in Eq. (3) vanishes. Recall that the sign of the separatrix value determines whether the homoclinic loop, here doublepulsed, is oriented or flipped. These codimensiontwo points globally organize the structure of the compact ()parameter region of the ShimizuMorioka model, including structural transformations of the Lorenz attractor in the model, including the emergence of stability islands – lacunae inside the attractor.
Figure 13 represents a partial ()diagram of the ShimizuMorioka model and depicts primary bifurcation curves and the corresponding phase portraits of the separatrix behaviors. Among those are the AndronovHopf bifurcation curve, , above which the equilibrium states, are stable, and saddlefoci below. This bifurcation is primarily supercritical, but becomes subcritical at smaller values of the parameter . The bifurcation curve, , corresponds to the formation of the homoclinic butterfly, or figure8 at larger values of . The transition between the branches of the curve is no bifurcation like the inclinationswitch because the saddle value, is negative here. The pathway at demonstrates the evolution of the simple MorseSmale dynamics of the model from the stable equilibrium states to a stable symmetric periodic orbit. This orbit emerges through a gluing bifurcation after that two stable periodic orbits existing below the supercritical AndronovHopf curve form the homoclinic butterfly with . The saddle value becomes positive to left from the codimensiontwo point, labeled by on the bifurcation curve . The left segment of the curve is similar to the bifurcation curve, , of the Lorenz equation (Fig. 3); namely, the homoclinic butterfly with causes the homoclinic explosion in the phase space of the ShimizuMorioka model. The curve ( in Fig. 13) being an analog of the curve in the bifurcation diagram in Fig. 3 is the upper boundary of the existence region of the Lorenz attractor is the parameter plane of this model. Below the separatrices of the saddle no longer converge to stable periodic orbits but fill in the strange attractor, like in Fig. 4. The bifurcation unfolding of the homoclinic resonant saddle includes various bifurcation curves: among them in the figure we depict, in addition to , the curve, , corresponding to saddlenode of merging stable (through the supercritical AndronovHopf bifurcation) and saddle (through the homoclinic bifurcation on ) periodic orbits. The curve labeled by corresponds to a pair of the doublepulsed homoclinic loops. Continued further away from the codimensiontwo point, , the curve frames the chaotic region in the parameter plane. The point, , on this curve is a codimensiontwo orbitflip homoclinic bifurcation: to the left of it, the loops become flipped like the median of a Möbius band. The dashed line, originating from the point is alike the curve in Fig. 9 for the Lorenz equation. This curve is the second boundary of the existence region (above) of the Lorenz attractor in the ShimizuMorioka model, that separates wild chaos with homoclinic tangencies (below). This like passes through the primary Tpoint in the parameter plane. Figure 13 depicts a few more bifurcation curves originating from the point, : two fourpulses homoclinic curve terminating at subsidiary Tpoints and the curve labeled by “” corresponding to a period doubling bifurcation. An intersection of the dashed line with a homoclinic bifurcation curve (see Fig. 13) corresponds to another orbitflip bifurcation and so forth.
Indeed the skeleton of the bifurcation set of the ShimizuMorioka is more complex. The detailed bifurcation diagram is shown in the top panel of Fig. 13. It reveals several Tpoints, the pitchfork bifurcation curve, , among other bifurcation curves for various homoclinic and heteroclinic connections. The detailed description of the bifurcation structure of the ShimizuMorioka model is out of scope of this paper. The curious reader can find a wealth of information on bifurcations of the Lorenz attractor in the original papers [Shilnikov, 1986, 1989, 1993; Shilnikov et al., 1993]. We point out that those bifurcation curves were continued in the parameter plane of the model using A. Shilnikov’s homemade software also based on the symbolic kneading toolbox.
The bottom panel of Fig. 13 is a defacto proof of the new kneading invariant mapping technique. The panel represents the biparametric scan in color of the dynamics of the ShimizuMorioka model that is based on the evaluation of the first 50 kneadings of the separatrix of the saddle on the grid of points in the parameter region. Getting the mapping took about one hour on a highend workstation without any parallelization efforts. The color scan reveals a plethora of largescale Tpoints, as well as nearby smaller ones (Fig. 14) invisible in the given parameter range, as well as the saddles separating spiral structures.
The solidcolor zones in the mapping correspond to simple MorseSmale dynamics in the model. These trivial dynamics are due to either the separatrix converging to the stable focus () and emergent periodic orbit with the same kneading invariant (brown region), or to the symmetric and asymmetric stable figure8 periodic orbits (dark blue region). The borderlines between the simple and complex dynamics in the ShimizuMorioka model are clearly demarcated. On the top it is the curve, , (see the top panel of Fig. 13). The transition from the stable 8shaped periodic orbits to the Lorenz attractor (through the boundary, ) is similar though more complicated as it involves a pitchfork bifurcation and bifurcations of doublepulsed homoclinics, see [Shilnikov, 1993; Shilnikov et al., 1993] for details.
One can clearly see the evident resemblance between both diagrams found using the bifurcationaly exact numerical methods and by scanning the dynamics of the model using the proposed kneading invariant technique. The latter reveals a richer structure providing finer details. The structure can be enhanced further by examining longer tails of the kneading sequences. This allows for the detection of smallerscale spiral structures within scrolls of the primary Tvortices, as predicted by the theory [Bykov, 1993]. Figure 14 shows a magnification of the scan of the ShimizuMorioka model near the primary Tpoint that depicts several other smallscale Tpoints.
Finally, we compare the new kneading scanning apparatus with the customary biparametric sweeping (shown in Fig. 15) of the ShimizuMorioka model that is based on the evaluation of the Lyapunov exponent spectrum computed over a finite time interval [Gonchenko et al., 2005]. Likewise the case of the Lorenz model, the sweeping based on the Lyapunov exponents shows no sign of spiral or saddle structures inside the region of deterministic chaos. The regions of the solid colors are associated with the sign of the largest Lyapunov exponent, : negative values correspond to steady state attractors in the green region; corresponds to periodic attractors in the blue region; and is associated with chaotic dynamics in the model in the region. Note blue islands in the redcolored region that correspond to stability windows in chaosland. In such windows the Lorenz attractor has an emergent lacuna containing, initially, a single symmetric saddle periodic orbit. The orbit then undergoes a pitchfork bifurcation that makes it stable. The basin of the stable orbit, which is first bounded by the 2D stable manifold of two asymmetric saddle periodic orbits, increases so that the stable orbit starts to dominate over chaotic dynamics in the corresponding stability window. These bifurcations underlie the transition from simple dynamics (blue region) due to the symmetric stable periodic orbit to chaos through the curve, , as the parameter in decreased.
6 6D optically pumped laser model
The coexistence of multiple Tpoints and accompanying fractal structures in the parameter plane is a signature for systems with the Lorenz attractor. A question though remains whether the new computational technique will work for systems of dimensions higher than three. In fact, to apply the technique to a generic Lorenzlike system, only wave forms of a symmetric variable progressing in time, that consistently starts from the same initial condition near the saddle is required. Next is an example from nonlinear optics – a 6D model of the optically pumped, far infrared red threelevel molecular laser [Moloney et al., 1989; Forysiak et al., 1991] given by
(7) 
Here, and are the Rabi flopping quantities representing the electric field amplitudes at pump and emission frequencies. The parameter is a natural bifurcation parameter as it is easily varied experimentally. The second bifurcation parameter, , can be varied to some degree at the laboratory by the addition of a buffer gas. This system presents, like the Lorenz equations, a symmetry . The laser model has either a single central equilibrium state, (with ), or additionally, through a pitchfork bifurcation, a pair of symmetric equilibrium states, (with ); the stability of the equilibria depends on the parameter values.
Optically pumped, far infrared lasers are known to demonstrate a variety of nonlinear dynamic behaviors, including Lorenzlike chaos [Hubner et al., 1995]. An acclaimed example of the modeling studies of chaos in nonlinear optics is the two level laser Haken model [Haken, 1975] to which the Lorenz equation can be reduced. A validity that three level laser models would inherently persist the Lorenz dynamics were widely questioned back then. It was comprehensively demonstrated in [Forysiak et al., 1991] in 1991 that this plausible laser model possesses a plethora of dynamical and structural features of the Lorenzlike systems, including the Lorenz attractor per se (with lacunae as well), similar AndronovHopf, pitchfork, various homoclinic and heteroclinic bifurcations including codimensiontwo ones Tpoints, see Fig. 16. Similar structures were also discovered in another nonlinear optics model for a laser with a saturable absorber which can be reduced to the ShimizuMiorioka model near a steady state solution with triple zero exponents [Vladimirov & Volkov, 1993]
The laser model (7) is quite rich in bifurcations; the list also includes a supercritical AndronovHopf bifurcation of the central equilibrium state that gives rise to a stable figure8 shaped periodic orbit (in proper projections) for the parameter values to the left of the bifurcation curve, , in the bifurcation diagrams shown in panels B and C of Fig. 17. Observe from the diagram that the curve originates from the point labeled, . This point corresponds to a codimensiontwo KhorozovTakens bifurcation of an equilibrium state with two zero Lyapunov exponents. The bifurcation is an extension of the BogdadovTakens bifurcation on a symmetric central manifold. The unfolding of this bifurcation includes two more curves: standing for the AndronovHopf bifurcation for the secondary equilibrium states, ; and standing for a homoclinic connection made of two symmetric separatrix loops biasymptotic to the saddle, . The continuation of the curve, , away from reveals rather peculiar details that substantially organize the bifurcation diagram of this laser model. Near the curve, , corresponds to a homoclinic figure8 of the saddle with a negative saddle value, on the symmetric 2D central manifold in the 6D phase space of the model. Recall that the figure8 homoclinic connection stands for the case where the 1D unstable separatrices come back to the saddle, from the symmetrically opposite directions along the eigenvector corresponding to the leading stable exponent at the equilibrium state. This bifurcation gluing two stable periodic orbits emerging from through the supercritical AndronovHopf bifurcation gives rise to the stable symmetric figure8 periodic orbit existing nearby . As the curve, , is continued further away from , the stable leading direction at the saddle, , changes: it becomes the invariant axis (like the axis in the ShimizuMorioka model) so that the separatrix loops start returning tangent to each other and hence forming the homoclinic butterfly. Nevertheless, this is a gluing bifurcation, not a codimensiontwo bifurcation of the change of the leading direction (inclination switch) as the saddle value remains negative on this branch of the curve, . The saddle value vanishes, making the saddle resonant at the codimensiontwo point, and becomes positive further down on . As the curve is continued, the homoclinic butterfly undergoes another codimensiontwo orbitflip bifurcation so that the separatrices loops of the saddle, become nonorientated. As the result, further down the curve, each loop gains an extra turn around the incoming separatrix of the opposite saddlefocus, i.e. becomes a doublepulsed one with the kneading. Depending on the parameter cut there are two scenarios for termination of the curve, , in the diagram (compare the bifurcation diagrams in panels B and C of Fig. 17): first, when it terminates at the codimensiontwo Tpoint corresponding to the heteroclinic connection between all saddle equilibria, and as shown in panel B. The second scenario for is less predictable at : the branch, , of doublepulsed separatrix loops ends up at the codimensiontwo point of the resonant saddle with the zero saddlevalue (panel C). The answer to the question what makes the curve change its destination is a saddle point in the parameter diagram that the kneading scan reveals in Fig. 13. By varying the parameter cut in the 3D parameter space, this bifurcation curve is destined by the saddle to finish at either terminal point, see details in [Shilnikov et al., 1993]. In the case where it spirals onto the Tpoint, there is another bifurcation curve corresponding to the same kneading that connects the codimensiontwo orbitswitch point and the point corresponding to the resonant saddle located on the curve .
The panel A in Fig. 17 represents the kneading scans of the dynamics of the laser model which is mapped onto the parameter plane with and . The scan is done using the same 50 kneading entries. It has the regions of chaotic dynamics clearly demarcated from the solid color windows of persistent kneadings corresponding to trivial attractors such as stable equilibria and periodic orbits. The region of chaos has a vivid fractal spiral structure featuring a chain of Tpoints. Observe also a thin chaotic layer bounded away from the curve by a curve of doublepulsed homoclinics with the kneading connecting the codimensiontwo points: the resonant saddle and the orbitflip both on . One feature of these points is the occurrence of the Lorenz attractor with one or more lacunae [Afraimovich et al., 1983; Shilnikov, 1986, 1993; Shilnikov et al., 1993]. Such a strange attractor with a single lacuna containing a figure8 periodic orbit in the phase space of the given laser model is shown in panel A of Fig.16.
7 Conclusions
We have demonstrated the new proposed computational toolkit for thorough explorations of chaotic dynamics in three exemplary models with the Lorenz attractor. The algorithmically easy though powerful toolkit in based on the scanning technique that maps the dynamics of the system onto the biparametric plane. The core of the approach is the evaluation of the kneading invariants for regularly or chaotically varying flipflop patterns of a single trajectory – the separatrix of the saddle singularity in the system. In the theory, the approach allows two systems with the structurally unstable Lorenz attractors to be conjugated with the mean of a single number – the kneading invariant. By using readyforuse tools in Matlab, we could have the parameter plane of the model in question be foliated by the level curves of distinct colors corresponding to distinct values of the kneading invariants. The kneading scans revel unambiguously the key features in the Lorenzlike systems such as a plethora of underlying spiral structures around Tpoints, separating saddles in intrinsically fractal regions corresponding to complex chaotic dynamics. We point out that no other techniques, including approaches based on the Lyapunov exponents, can reveal the discovered parametric chaos with such stunning clarity and beauty. Figure 18 for the ShimizuMorioka model shows that a fine scan based on the finitetime Lyapunov exponents is able to indicate some presence of spiral and saddle structures and differentiate between the chaotic dynamics due to the Lorenz attractor and those due to additional degrees of instability brought in by saddlefoci.
The kneading based methods should be beneficial for detailed studies of other systems admitting reasonable symbolic descriptions. It bears an educational aspect as well: the kneading based scanning can be used for inclass presentation to reveal the fascinating complexity of loworder models in the crossdisciplinary field of nonlinear dynamics. The biparametric mapping technique can be easily adopted for parallel multicore GPU platforms allowing for ultrafast simulations of models in questions. Additional implementation of highprecision computations of long transients shall thoroughly reveal multilayer complexity of selfsimilar fractal patterns near Tpoint vortices. In future research we intend to enhance and refine the toolkits for exploration of other symmetric and asymmetric [Shilnikov & Shilnikov, 1991] systems of differential and difference equations, like 3D Poincaré mappings [Shilnikov et al., 1993; Gonchenko et al., 2005], including 4D models with saddlefoci, that require two and more kneading invariants for the comprehensive symbolic description.
8 Acknowledgments
This work is supported by the Spanish Research project MTM200910767 (to R.B.), and by NSF grant DMS1009591, RFFI Grant No. 080100083, GSU Brain & Behaviors program, and MESRF project 14.740.11.0919 (to A.S) as well as by the Grant 11.G34.31.0039 of the Government of the Russian Federation for state support of scientific research conducted under supervision of leading scientists in Russian educational institutions of higher professional education (to L.P.S). We thank Dima Turaev for stimulating discussions, Yuri Kuznetsov for sharing the data used in Fig. 8. We thank Aaron Kelley and Jeremy Wojcik for careful proofreading the manuscript.
References
 Abad et al. [2011a] Abad, A., Barrio, R., Blesa, F. & Rodríguez, M. [2011a] “TIDES, a Taylor Integrator for Differential EquationS,” ACM Transactions on Mathematical Software, in press.
 Abad et al. [2011b] Abad, A., Barrio, R., Blesa, F. & Rodríguez, M. [2011b] “TIDES web page: http://gme.unizar.es/software/tides”.
 Afraimovich et al. [1977] Afraimovich, V., Bykov, V. V. & Shilnikov, L. P. [1977] “The origin and structure of the Lorenz attractor,” Sov. Phys. Dokl. 22, 253–255.
 Afraimovich et al. [1983] Afraimovich, V., Bykov, V. V. & Shilnikov, L. P. [1983] “On structurally unstable attracting limit sets of Lorenz attractor type,” Trans. Moscow Math. Soc. 44, 153–216.
 Afraimovich & Shilnikov [1983] Afraimovich, V. S. & Shilnikov, L. P. [1983] “Strange attractors and quasiattractors,” Nonlinear dynamics and turbulence, Interaction Mech. Math. Ser. (Pitman, Boston, MA), pp. 1–34.
 Barrio et al. [2009] Barrio, R., Blesa, F. & Serrano, S. [2009] ‘‘Qualitative analysis of the Rössler equations: bifurcations of limit cycles and chaotic attractors,” Phys. D 238, 1087–1100.
 Barrio et al. [2011a] Barrio, R., Blesa, F. & Serrano, S. [2011a] “Behavior patterns in multiparametric dynamical systems: Lorenz model,” Inter. J. Bif. Chaos , in press.
 Barrio et al. [2011b] Barrio, R., Blesa, F., Serrano, S. & Shilnikov, A. [2011b] “Global organization of spiral structures in biparameter space of dissipative systems with Shilnikov saddlefoci,” Phys. Rev. E 84, 035201.
 Barrio et al. [2011c] Barrio, R., Rodríguez, M., Abad, A. & Blesa, F. [2011c] “Breaking the limits: the Taylor series method,” Appl. Math. Comput. 217, 7940–7954.
 Barrio & Serrano [2007] Barrio, R. & Serrano, S. [2007] “A threeparametric study of the Lorenz model,” Phys. D 229, 43–51.
 Barrio & Serrano [2009] Barrio, R. & Serrano, S. [2009] “Bounds for the chaotic region in the Lorenz model,” Phys. D 238, 1615–1624.
 Bykov [1998] Bykov, V. [1998] “In bifurcations leasing to chaos in chua’s circuit.” Inter. J. Bif. Chaos 8, 685–699.
 Bykov [1993] Bykov, V. V. [1993] “The bifurcations of separatrix contours and chaos,” Phys. D 62, 290–299.
 Bykov & Shilnikov [1989] Bykov, V. V. & Shilnikov, A. [1989] “Boundaries of the domain of existence of a Lorenz attractor,” Methods in qualitative theory and bifurcation theory (in Russian) , 151–159.
 Bykov & Shilnikov [1992] Bykov, V. V. & Shilnikov, A. L. [1992] “On the boundaries of the domain of existence of the Lorenz attractor,” Selecta Math. Soviet. 11, 375–382.
 FernándezSánchez et al. [2002] FernándezSánchez, F., Freire, E. & RodríguezLuis, A. J. [2002] “Tpoints in a symmetric electronic oscillator. I. Analysis,” Nonlinear Dynam. 28, 53–69.
 Forysiak et al. [1991] Forysiak, W., Moloney, J. & Harrison, R. [1991] “Bifurcations of an optically pumped threelevel laser model,” Physica D 53, 162–186.
 Gallas [2010] Gallas, J. A. C. [2010] “The structure of infinite periodic and chaotic hub cascades in phase diagrams of simple autonomous flows,” Inter. J. Bif. Chaos 20, 197–211.
 Gaspard et al. [1984] Gaspard, P., Kapral, R. & Nicolis, G. [1984] “Bifurcation phenomena near homoclinic systems: A twoparameter analysis,” J. Stat. Phys. 35, 697–727.
 Glendinning & Hall [1996] Glendinning, P. & Hall, T. [1996] “Zeros of the kneading invariant and topological entropy for Lorenz maps,” Nonlinearity 9, 999–1014.
 Glendinning & Sparrow [1984] Glendinning, P. & Sparrow, C. [1984] “Local and global behavior near homoclinic orbits,” J. Stat. Phys. 35, 645–696.
 Glendinning & Sparrow [1986] Glendinning, P. & Sparrow, C. [1986] “Tpoints: a codimension two heteroclinic bifurcation,” J. Stat. Phys. 43, 479–488.
 Glendinning & Sparrow [1993] Glendinning, P. & Sparrow, C. [1993] “Prime and renormalisable kneading invariants and the dynamics of expanding lorenz maps.” Physica D: Nonlinear Phenomena 1–4, 22–50.
 Gonchenko et al. [2005] Gonchenko, S., Ovsyannikov, I., Simo, C. & Turaev, D. [2005] “Threedimensional Hénonlike maps and wild Lorenzlike attractors,” Inter. J. Bif. Chaos 15(11), 3493–3508.
 Gonchenko et al. [1993] Gonchenko, S. V., Shil’nikov, L. P. & Turaev, D. V. [1993] “On models with nonrough poincarï¿½ homoclinic curves,” Physica D 14, 1–14.
 Gonchenko et al. [1996] Gonchenko, S. V., Shil’nikov, L. P. & Turaev, D. V. [1996] ‘‘Dynamical phenomena in systems with structurally unstable Poincare homoclinic orbits.” Chaos 6, 15–31.
 Guckenheimer & Williams [1979] Guckenheimer, J. & Williams, R. F. [1979] “Structural stability of Lorenz attractors,” Inst. Hautes Études Sci. Publ. Math. 50, 59–72.
 Haken [1975] Haken, H. [1975] “Analogy between higher instabilities in fluids and lasers,” Phys. Letters A 53, 77–78.
 Hubner et al. [1995] Hubner, U., Weiss, C., Abraham, N. & Tang, D. [1995] “Lorenzlike chaos in nh3fir lasers,” Infrared Phys. Technol. 36(1), 489–512.
 Kaplan & Yorke [1979] Kaplan, J. L. & Yorke, J. A. [1979] “Preturbulence: a regime observed in a fluid flow model of Lorenz,” Comm. Math. Phys. 67, 93–108.
 Li & Malkin [2003] Li, M.C. & Malkin, M. [2003] “Smooth symmetric and Lorenz models for unimodal maps,” Inter. J. Bif. Chaos 13, 3353–3371.
 Lorenz [1963] Lorenz, E. [1963] “Deterministic nonperiodic flow,” J. Atmospheric Sci. 20, 130–141.
 Lorenz [2008] Lorenz, E. N. [2008] “Compound windows of the Hénonmap,” Phys. D 237, 1689–1704.
 Malkin [1991] Malkin, M. [1991] “Rotation intervals and dynamics of lorenz type mappings.” Selecta Math. Sovietica 10, 265–275.
 Milnor & Thurston [1988] Milnor, J. & Thurston, W. [1988] “On iterated maps of the interval,” Lecture Notes in Math. 1342, 465–563.
 Moloney et al. [1989] Moloney, J., Forysiak, W., Uppal, T. & Harrison, R. [1989] “Regular and chaotic dynamics of optically pumped molecular lasers.” Phys Rev A 39, 1277–1285.
 Pertovksaya & Yudovich [1980] Pertovksaya, N. & Yudovich, V. [1980] ‘‘Homolinic loops on the ZaltzmanLorenz system,” Methods of Qualitative Theory of Differential Equations, Gorky University. , 73–83.
 Rand [1978] Rand, D. [1978] “The topological classification of Lorenz attractors,” Mathematical Proceedings of the Cambridge Philosophical Society 83(03), 451–460.
 Robinson [1989] Robinson, C. [1989] “Homoclinic bifurcation to a transitive attractor of Lorenz type.” Nonlinearity 2, 495–518.
 Roshchin [1978] Roshchin, N. [1978] “Dangerous stability boundaries in Lorenz’s model,” Prikl. Mat. Mekhanika 42, 950.
 Rychlic [1990] Rychlic, M. [1990] “Lorenz attractor through Shil’nikov type bifurcation I.” Erof. Theory and Dyn. Systems 10, 793–821.
 Shilnikov [1986] Shilnikov, A. [1986] “Bifurcations and chaos in the MariokaShimizu model. Part I,” Methods in qualitative theory and bifurcation theory (in Russian) , 180–193.
 Shilnikov [1989] Shilnikov, A. [1989] “Bifurcations and chaos in the MariokaShimizu model. Part II,” Methods in qualitative theory and bifurcation theory (in Russian) , 130–138.
 Shilnikov [1991] Shilnikov, A. [1991] “Bifurcation and chaos in the MoriokaShimizu system,” Selecta Math. Soviet. 10(2), 105–117.
 Shilnikov [1993] Shilnikov, A. [1993] “On bifurcations of the Lorenz attractor in the ShimizuMorioka model,” Physica D 62(14), 338–346.
 Shilnikov & Shilnikov [1991] Shilnikov, A. & Shilnikov, L. [1991] “On the nonsymmetric Lorenz model,” Inter. J. Bif. Chaos 1, 773–776.
 Shilnikov et al. [1993] Shilnikov, A. L., Shilnikov, L. P. & Turaev, D. V. [1993] “Normal forms and Lorenz attractors,” Inter. J. Bif. Chaos 3, 1123–1139.
 Shilnikov [1967] Shilnikov, L. [1967] ‘‘The existence of a denumerable set of periodic motions in fourdimensional space in an extended neighborhood of a saddlefocus.” Soviet Math. Dokl. 8(1), 54–58.
 Shilnikov [1968] Shilnikov, L. [1968] “On the birth of a periodic motion from a trajectory biasymptotic to an equilibrium state pf the saddle type.” Soviet Math. Sbornik. 35(3), 240–264.
 Shilnikov [1980] Shilnikov, L. [1980] “Bifurcation theory and the Lorenz model.” Appendix to Russian edition of “The Hopf Bifurcation and Its Applications.” Eds. J. Marsden and M. McCraken , 317–335.
 Shilnikov [1981] Shilnikov, L. [1981] “The theory of bifurcations and quasiattractors,” Uspeh. Math. Nauk 36(4), 240–242.
 Shilnikov [1994] Shilnikov, L. [1994] “Chua’s circuit: Rigorous results and future problems,” Inter. J. Bif. Chaos 4(3), 489–518.
 Shilnikov [1997] Shilnikov, L. [1997] “Mathematical problems of nonlinear dynamics: A tutorial. Visions of nonlinear mechanics in the 21st century,” Journal of the Franklin Institute. 334(56), 793–864.
 Shilnikov [2002] Shilnikov, L. [2002] “Bifurcations and strange attractors,” Proc. Int. Congress of Mathematicians, Beijing (China) (Invited Lectures). 3, 349–372.
 Shilnikov & Shilnikov [2007] Shilnikov, L. & Shilnikov, A. [2007] “Shilnikov bifurcation,” Scholarpedia 2(8), 1891.
 Shilnikov [1965] Shilnikov, L. P. [1965] “A case of the existence of a countable number of periodic motions,” Sov. Math. Dokl. 6, 163.
 Shilnikov et al. [1998,2001] Shilnikov, L. P., Shilnikov, A. L., Turaev, D. & Chua, L. O. [1998,2001] Methods of qualitative theory in nonlinear dynamics. Parta I and II (World Scientific Publishing Co. Inc.).
 Shimizu & Morioka [1980] Shimizu, T. & Morioka, N. [1980] “On the bifurcation of a symmetric limit cycle to an asymmetric one in a simple model,” Physics Letters A 76, 201 – 204.
 Sparrow [1982] Sparrow, C. [1982] The Lorenz equations: bifurcations, chaos, and strange attractors, Applied Mathematical Sciences, Vol. 41 (SpringerVerlag, New York).
 Tresser & Williams [1993] Tresser, C. & Williams, R. [1993] “Splitting words and Lorenz braids,” Physica D: Nonlinear Phenomena 62(1–4), 15–21.
 Turaev & Shilnikov [1998] Turaev, D. & Shilnikov, L. [1998] “An example of a wild strange attractor,” Sbornik. Math. 189(2), 291–314.
 Vladimirov & Volkov [1993] Vladimirov, A. & Volkov, D. [1993] “Lowintensity chaotic operations of a laser with a saturable absorber,” Optics Communications 100(14), 351–360.
 Wieczorek & Krauskopf [2005] Wieczorek, S. & Krauskopf, B. [2005] “Bifurcations of homoclinic orbits in optically injected lasers,” Nonlinearity 18, 1095–1120.