Kneadings, Symbolic Dynamics and Painting Lorenz Chaos. a Tutorial
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 bi-parametric structures and detect their organizing centers – codimension-two T-points and separating saddles in the kneading-based 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, T-points, Lorenz attractor, chaos, homoclinic and heteroclinic orbits
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, Morse-Smale, 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 low-order 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 bi-parametric 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 time-discrete [Lorenz, 2008] and time-continuous 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 saddle-focus [Shilnikov & Shilnikov, 2007].
Application of the Lyapunov exponents technique fails, in general, to reveal fine structures embedded in the bi-parametric scans of Lorenz-like 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 Lorenz-like 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 bi-parametric parameter plane of the system in question are organized around the so-called T[erminal]-points corresponding to codimension-two or -higher heteroclinic connections between two or more saddle equilibria. For -symmetric systems with the Lorenz attractor, the degeneracy is reduced to a cod-2 bifurcation of a closed heteroclinic connection involving two saddle-foci and a saddle at the origin. T-points have been located in various models of diverse origins including electronic oscillators [Bykov, 1998; Fernández-Sá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 T-point in the parameter plane in a symmetric system [Bykov, 1993; Glendinning & Sparrow, 1986]; various T-points configurations for other heteroclinic connections were examined in detail in [Bykov, 1993]. Here, the heteroclinic connection is formed by a pair of symmetric saddle-foci and a saddle whose one-dimensional stable (incoming) and unstable (outgoing), respectively, separatrices merge only at the codimension-two T-point in the bi-parametric space of the Lorenz model. It follows from the theoretical analysis that the unfolding of this bifurcation includes three main curves ending at the T-point: a spiraling bifurcation curve, (1-hom ) corresponding to two simultaneous (due to the symmetry) homoclinic loops of the saddle, a curve, (1-hom ), corresponding to two simultaneous homoclinic loops of the saddle-foci, and another codimension-one bifurcation curve, (1-het ), corresponding to a heteroclinic correction between both saddle-foci. In addition, the unfolding includes infinitely many subsidiary T-points in between 1-hom and 1-het .
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 Lorenz-like systems that employs the kneading invariants [Milnor & Thurston, 1988]. We will then show how the toolkit can be used for detecting various structures in bi-parametric 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 Shimizu-Morioka model. Finally, In Section 6 the proposed technique will be tested on a 6D model of the optically pumped, far infrared red three-level 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 de-facto proof of deterministic chaos. The butterfly-shaped 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:
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 Uni-parametric cut through the Lorenz equation
An exploration of primary bifurcations in the Lorenz equation begins with a single-parameter 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 uni-parametric 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 pitch-fork 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 pre-turbulence 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 cylindrically-shaped 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 Andronov-Hopf 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 bi-parametric 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 pitch-fork 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 non-leading or strongly stable direction, in . After the pitch-fork 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 codimension-one, 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 so-called 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 figure-8 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 so-called 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 codimension-2 bifurcations, termed resonant saddle, orbit- and inclination-flip, correspondingly [Shilnikov et al., 1998,2001], as it occurs in the Shimizu-Morioka 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 cylinder-shaped 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 saddle-foci of the (1,2)-type through a subcritical Andronov-Hopf bifurcation [Roshchin, 1978]. The topological (1,2)-type means that each saddle-focus has 2D unstable and 1D stable manifolds; the latter is formed by two incoming separatrices. Some local properties of the saddle-foci 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 saddle-foci meet the Shilnikov condition [Shilnikov, 1965, 1967; Shilnikov & Shilnikov, 2007]. Therefore, as soon as the saddle-focus possesses a homoclinic loop, such a bifurcation causes the abundance of periodic orbits nearby. Those periodic orbits constantly undergo saddle-node and period doubling bifurcations as the parameters are varied. Moreover, since the divergence of the vector field of the Lorenz equation is always negative, saddle-node bifurcations give rise to stable periodic orbits near the homoclinic saddle-focus bifurcation. Under fulfillment of some global conditions, a single Shilnikov saddle-focus bifurcation can lead to the formation of a spiral or screw-like attractor. However, a strange attractor due to the Shilnikov saddle-focus 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 quasi-attractor thus referring to that because besides stable periodic orbits with weak basins it may have structurally unstable or non-transverse 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 so-called wild chaotic attractors [Shilnikov, 1994, 1997; Turaev & Shilnikov, 1998; Shilnikov, 2002].
The Lorenz attractor is non-hyperbolic 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 saddle-node 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] saddle-foci, , 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 finite-time transient on the chaotic attractor.
After the stable foci have lost the stability through the subcritical Andronov-Hopf 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 saddle-foci 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 one-dimensional 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 flip-flop pattern of switching between the butterfly wings centered around the saddle-foci. 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 flip-flops around the equilibrium states and in the 3D phase space of the Lorenz equation (Fig. 5). Or, alternatively, can be reduced to the time-evolution of the -coordinate of the separatrix, as shown in panel B of Fig. 5. The sign-alternation 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 Lorenz-like 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 Lorenz-like 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 cross-section 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:
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 flip-flopping 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 Lorenz-like 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:
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:
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 cusp-shaped 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:
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 off-lying 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 cross-sections, 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 saddle-foci 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 bi-parametric 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 bi-parametric scan of the Lorenz equation within a specific range in the -plane. The resolution of scans is set by using mesh grids of equally-distanced points. Next by accurately integrating the separatrix using Taylor series software TIDES111freeware 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 bi-parametric mapping: with some appropriately chosen , the value that determines the depth of the scan. The mapping is then colorized in Matlab by using various built-in non-linear 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 bi-parametric scans.
Figure 8 represents the kneading-based 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 Andronov-Hopf bifurcation and a pitch-fork bifurcation describing continuous transitioning between stable symmetric and asymmetric periodic orbits. This holds true for a period-doubling 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 saddle-node bifurcation, homoclinic bifurcations, and some other.
A borderline between two solid-color 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 bi-parametric screening based on the finite-time Lyapunov exponent approach (see below), and what we developed it for, is the detection of bifurcations within the Lorenz strange attractor. The corresponding yellow-reddish 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 T-points. Panel B of this figure depicts the kneading mapping near the left-bottom 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 yellow-reddish 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 pseudo-hyperbolic 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 T-point. 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 saddle-focus , followed by three revolutions around the saddle-focus , and then returns to the saddle at the origin. While moving along the spiraling curve toward the T-points, the separatrix makes progressively more turns around , or more precisely around the 1D incoming separatrix of this saddle-focus. With each incremental turn around , the separatrix comes closer to while the bifurcation curve becomes one scroll closer to the T-point simultaneously. Due to this feature the T-point is called a Terminal point. The T-point corresponds to the following symbolic sequence . In virtue of the symmetry of the Lorenz equation, the T-point actually corresponds to a closed heteroclinic connection involving all three saddle-equilibria, see Figs. 1 and 5. The merger of the right (left) separatrix of the saddle with the incoming separatrix of the saddle-focus , increases the codimension (degeneracy) of this heteroclinic bifurcation to two; note that intersections of the 2D unstable manifolds of the saddle-foci, 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 saddle-focus, as well as to a heteroclinic connection between both saddle-foci (see the sketch of the bifurcation unfolding of the T-point in Fig. 1). The corresponding bifurcation curves of these homoclinic and heteroclinic bifurcations originating rightward from the T-point bound a sector containing subsidiary T-points [Bykov, 1993; Glendinning & Sparrow, 1984]. Each new T-point produces other self-similar structures scaled like fractals. Panel C shows two such identified T-points: primary and secondary . The primary codimension-two T-point in the Lorenz equation was originally discovered by Yudovich [Pertovksaya & Yudovich, 1980].
As soon as the saddle-foci and their bifurcations become involved in the dynamics of the Lorenz equation near the primary T-point, the Lorenz attractor loses the purity of the genuine chaotic attractor that used to have neither stable periodic orbits nor non-transverse homoclinic trajectories. It transforms into a quasi-chaotic attractor with weakly stable orbits of small basins and nontransverse homoclinic orbits. The idea of non-transversality 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 T-point, thereby separating the existence region of the Lorenz attractor from bifurcations of the saddle-foci, and consequently from all subsidiary T-points 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 T-points. Note that basins of spiraling T-points are separated by corresponding saddles. One can spot self-similar smaller-scale spiral structures within large-scale ones and so forth. The richness of such fractal structures in the parameter plane results from the synergy of the Lorenz-like dynamics amplified by chaos induced by the Shilnikov saddle-foci.
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 bi-parametric 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 (yellow-reddish) is characterized by a small positive Lyapunov exponent, variations of which are not significant enough to reveal fine structures, as spiraling T-points. 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 large-scale spiral hubs around T-points, 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 Shimizu-Morioka model
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 stable-foci or saddle-foci of the (1,2)-topological type.
This model was originally introduced to examine a pitch-fork bifurcation of the stable figure-8 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 codimension-three 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 T-points 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 codimension-two: resonant saddle with the zero saddle value or the saddle index , and the orbit-flip bifurcation where the separatrix value in Eq. (3) vanishes. Recall that the sign of the separatrix value determines whether the homoclinic loop, here double-pulsed, is oriented or flipped. These codimension-two points globally organize the structure of the compact ()-parameter region of the Shimizu-Morioka 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 Shimizu-Morioka model and depicts primary bifurcation curves and the corresponding phase portraits of the separatrix behaviors. Among those are the Andronov-Hopf bifurcation curve, , above which the equilibrium states, are stable, and saddle-foci 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 figure-8 at larger values of . The transition between the branches of the curve is no bifurcation like the inclination-switch because the saddle value, is negative here. The pathway at demonstrates the evolution of the simple Morse-Smale 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 Andronov-Hopf curve form the homoclinic butterfly with . The saddle value becomes positive to left from the codimension-two 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 Shimizu-Morioka 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 saddle-node of merging stable (through the supercritical Andronov-Hopf bifurcation) and saddle (through the homoclinic bifurcation on ) periodic orbits. The curve labeled by corresponds to a pair of the double-pulsed homoclinic loops. Continued further away from the codimension-two point, , the curve frames the chaotic region in the parameter plane. The point, , on this curve is a codimension-two orbit-flip 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 Shimizu-Morioka model, that separates wild chaos with homoclinic tangencies (below). This like passes through the primary T-point in the -parameter plane. Figure 13 depicts a few more bifurcation curves originating from the point, : two four-pulses homoclinic curve terminating at subsidiary T-points 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 orbit-flip bifurcation and so forth.
Indeed the skeleton of the bifurcation set of the Shimizu-Morioka is more complex. The detailed bifurcation diagram is shown in the top panel of Fig. 13. It reveals several T-points, the pitch-fork bifurcation curve, , among other bifurcation curves for various homoclinic and heteroclinic connections. The detailed description of the bifurcation structure of the Shimizu-Morioka 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 home-made software also based on the symbolic kneading toolbox.
The bottom panel of Fig. 13 is a de-facto proof of the new kneading invariant mapping technique. The panel represents the bi-parametric scan in color of the dynamics of the Shimizu-Morioka 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 high-end workstation without any parallelization efforts. The color scan reveals a plethora of large-scale T-points, as well as nearby smaller ones (Fig. 14) invisible in the given parameter range, as well as the saddles separating spiral structures.
The solid-color zones in the mapping correspond to simple Morse-Smale 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 figure-8 periodic orbits (dark blue region). The borderlines between the simple and complex dynamics in the Shimizu-Morioka model are clearly demarcated. On the top it is the curve, , (see the top panel of Fig. 13). The transition from the stable 8-shaped periodic orbits to the Lorenz attractor (through the boundary, ) is similar though more complicated as it involves a pitch-fork bifurcation and bifurcations of double-pulsed 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 smaller-scale spiral structures within scrolls of the primary T-vortices, as predicted by the theory [Bykov, 1993]. Figure 14 shows a magnification of the scan of the Shimizu-Morioka model near the primary T-point that depicts several other small-scale T-points.
Finally, we compare the new kneading scanning apparatus with the customary bi-parametric sweeping (shown in Fig. 15) of the Shimizu-Morioka 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 red-colored region that correspond to stability windows in chaos-land. In such windows the Lorenz attractor has an emergent lacuna containing, initially, a single symmetric saddle periodic orbit. The orbit then undergoes a pitch-fork 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 T-points 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 Lorenz-like 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 three-level molecular laser [Moloney et al., 1989; Forysiak et al., 1991] given by
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 pitch-fork 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 Lorenz-like 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 Lorenz-like systems, including the Lorenz attractor per se (with lacunae as well), similar Andronov-Hopf, pitchfork, various homoclinic and heteroclinic bifurcations including codimension-two ones T-points, 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 Shimizu-Miorioka 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 super-critical Andronov-Hopf bifurcation of the central equilibrium state that gives rise to a stable figure-8 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 codimension-two Khorozov-Takens bifurcation of an equilibrium state with two zero Lyapunov exponents. The bifurcation is an extension of the Bogdadov-Takens bifurcation on a symmetric central manifold. The unfolding of this bifurcation includes two more curves: standing for the Andronov-Hopf bifurcation for the secondary equilibrium states, ; and standing for a homoclinic connection made of two symmetric separatrix loops bi-asymptotic 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 figure-8 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 figure-8 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 Andronov-Hopf bifurcation gives rise to the stable symmetric figure-8 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 Shimizu-Morioka 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 codimension-two 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 codimension-two point, and becomes positive further down on . As the curve is continued, the homoclinic butterfly undergoes another codimension-two orbit-flip bifurcation so that the separatrices loops of the saddle, become non-orientated. As the result, further down the curve, each loop gains an extra turn around the incoming separatrix of the opposite saddle-focus, i.e. becomes a double-pulsed 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 codimension-two T-point 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 double-pulsed separatrix loops ends up at the codimension-two point of the resonant saddle with the zero saddle-value (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 T-point, there is another bifurcation curve corresponding to the same kneading that connects the codimension-two orbit-switch 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 T-points. Observe also a thin chaotic layer bounded away from the curve by a curve of double-pulsed homoclinics with the kneading connecting the codimension-two points: the resonant saddle and the orbit-flip 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 figure-8 periodic orbit in the phase space of the given laser model is shown in panel A of Fig.16.
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 bi-parametric plane. The core of the approach is the evaluation of the kneading invariants for regularly or chaotically varying flip-flop 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 ready-for-use 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 Lorenz-like systems such as a plethora of underlying spiral structures around T-points, 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 Shimizu-Morioka model shows that a fine scan based on the finite-time 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 saddle-foci.
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 in-class presentation to reveal the fascinating complexity of low-order models in the cross-disciplinary field of nonlinear dynamics. The bi-parametric mapping technique can be easily adopted for parallel multicore GPU platforms allowing for ultra-fast simulations of models in questions. Additional implementation of high-precision computations of long transients shall thoroughly reveal multi-layer complexity of self-similar fractal patterns near T-point 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 saddle-foci, that require two and more kneading invariants for the comprehensive symbolic description.
This work is supported by the Spanish Research project MTM2009-10767 (to R.B.), and by NSF grant DMS-1009591, RFFI Grant No. 08-01-00083, 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.
- 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.  Afraimovich, V., Bykov, V. V. & Shilnikov, L. P.  “The origin and structure of the Lorenz attractor,” Sov. Phys. Dokl. 22, 253–255.
- Afraimovich et al.  Afraimovich, V., Bykov, V. V. & Shilnikov, L. P.  “On structurally unstable attracting limit sets of Lorenz attractor type,” Trans. Moscow Math. Soc. 44, 153–216.
- Afraimovich & Shilnikov  Afraimovich, V. S. & Shilnikov, L. P.  “Strange attractors and quasiattractors,” Nonlinear dynamics and turbulence, Interaction Mech. Math. Ser. (Pitman, Boston, MA), pp. 1–34.
- Barrio et al.  Barrio, R., Blesa, F. & Serrano, S.  ‘‘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 bi-parameter space of dissipative systems with Shilnikov saddle-foci,” 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  Barrio, R. & Serrano, S.  “A three-parametric study of the Lorenz model,” Phys. D 229, 43–51.
- Barrio & Serrano  Barrio, R. & Serrano, S.  “Bounds for the chaotic region in the Lorenz model,” Phys. D 238, 1615–1624.
- Bykov  Bykov, V.  “In bifurcations leasing to chaos in chua’s circuit.” Inter. J. Bif. Chaos 8, 685–699.
- Bykov  Bykov, V. V.  “The bifurcations of separatrix contours and chaos,” Phys. D 62, 290–299.
- Bykov & Shilnikov  Bykov, V. V. & Shilnikov, A.  “Boundaries of the domain of existence of a Lorenz attractor,” Methods in qualitative theory and bifurcation theory (in Russian) , 151–159.
- Bykov & Shilnikov  Bykov, V. V. & Shilnikov, A. L.  “On the boundaries of the domain of existence of the Lorenz attractor,” Selecta Math. Soviet. 11, 375–382.
- Fernández-Sánchez et al.  Fernández-Sánchez, F., Freire, E. & Rodríguez-Luis, A. J.  “T-points in a -symmetric electronic oscillator. I. Analysis,” Nonlinear Dynam. 28, 53–69.
- Forysiak et al.  Forysiak, W., Moloney, J. & Harrison, R.  “Bifurcations of an optically pumped three-level laser model,” Physica D 53, 162–186.
- Gallas  Gallas, J. A. C.  “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.  Gaspard, P., Kapral, R. & Nicolis, G.  “Bifurcation phenomena near homoclinic systems: A two-parameter analysis,” J. Stat. Phys. 35, 697–727.
- Glendinning & Hall  Glendinning, P. & Hall, T.  “Zeros of the kneading invariant and topological entropy for Lorenz maps,” Nonlinearity 9, 999–1014.
- Glendinning & Sparrow  Glendinning, P. & Sparrow, C.  “Local and global behavior near homoclinic orbits,” J. Stat. Phys. 35, 645–696.
- Glendinning & Sparrow  Glendinning, P. & Sparrow, C.  “T-points: a codimension two heteroclinic bifurcation,” J. Stat. Phys. 43, 479–488.
- Glendinning & Sparrow  Glendinning, P. & Sparrow, C.  “Prime and renormalisable kneading invariants and the dynamics of expanding lorenz maps.” Physica D: Nonlinear Phenomena 1–4, 22–50.
- Gonchenko et al.  Gonchenko, S., Ovsyannikov, I., Simo, C. & Turaev, D.  “Three-dimensional Hénon-like maps and wild Lorenz-like attractors,” Inter. J. Bif. Chaos 15(11), 3493–3508.
- Gonchenko et al.  Gonchenko, S. V., Shil’nikov, L. P. & Turaev, D. V.  “On models with non-rough poincarï¿½ homoclinic curves,” Physica D 1-4, 1–14.
- Gonchenko et al.  Gonchenko, S. V., Shil’nikov, L. P. & Turaev, D. V.  ‘‘Dynamical phenomena in systems with structurally unstable Poincare homoclinic orbits.” Chaos 6, 15–31.
- Guckenheimer & Williams  Guckenheimer, J. & Williams, R. F.  “Structural stability of Lorenz attractors,” Inst. Hautes Études Sci. Publ. Math. 50, 59–72.
- Haken  Haken, H.  “Analogy between higher instabilities in fluids and lasers,” Phys. Letters A 53, 77–78.
- Hubner et al.  Hubner, U., Weiss, C., Abraham, N. & Tang, D.  “Lorenz-like chaos in nh3-fir lasers,” Infrared Phys. Technol. 36(1), 489–512.
- Kaplan & Yorke  Kaplan, J. L. & Yorke, J. A.  “Preturbulence: a regime observed in a fluid flow model of Lorenz,” Comm. Math. Phys. 67, 93–108.
- Li & Malkin  Li, M.-C. & Malkin, M.  “Smooth symmetric and Lorenz models for unimodal maps,” Inter. J. Bif. Chaos 13, 3353–3371.
- Lorenz  Lorenz, E.  “Deterministic nonperiodic flow,” J. Atmospheric Sci. 20, 130–141.
- Lorenz  Lorenz, E. N.  “Compound windows of the Hénon-map,” Phys. D 237, 1689–1704.
- Malkin  Malkin, M.  “Rotation intervals and dynamics of lorenz type mappings.” Selecta Math. Sovietica 10, 265–275.
- Milnor & Thurston  Milnor, J. & Thurston, W.  “On iterated maps of the interval,” Lecture Notes in Math. 1342, 465–563.
- Moloney et al.  Moloney, J., Forysiak, W., Uppal, T. & Harrison, R.  “Regular and chaotic dynamics of optically pumped molecular lasers.” Phys Rev A 39, 1277–1285.
- Pertovksaya & Yudovich  Pertovksaya, N. & Yudovich, V.  ‘‘Homolinic loops on the Zaltzman-Lorenz system,” Methods of Qualitative Theory of Differential Equations, Gorky University. , 73–83.
- Rand  Rand, D.  “The topological classification of Lorenz attractors,” Mathematical Proceedings of the Cambridge Philosophical Society 83(03), 451–460.
- Robinson  Robinson, C.  “Homoclinic bifurcation to a transitive attractor of Lorenz type.” Nonlinearity 2, 495–518.
- Roshchin  Roshchin, N.  “Dangerous stability boundaries in Lorenz’s model,” Prikl. Mat. Mekhanika 42, 950.
- Rychlic  Rychlic, M.  “Lorenz attractor through Shil’nikov type bifurcation I.” Erof. Theory and Dyn. Systems 10, 793–821.
- Shilnikov  Shilnikov, A.  “Bifurcations and chaos in the Marioka-Shimizu model. Part I,” Methods in qualitative theory and bifurcation theory (in Russian) , 180–193.
- Shilnikov  Shilnikov, A.  “Bifurcations and chaos in the Marioka-Shimizu model. Part II,” Methods in qualitative theory and bifurcation theory (in Russian) , 130–138.
- Shilnikov  Shilnikov, A.  “Bifurcation and chaos in the Morioka-Shimizu system,” Selecta Math. Soviet. 10(2), 105–117.
- Shilnikov  Shilnikov, A.  “On bifurcations of the Lorenz attractor in the Shimizu-Morioka model,” Physica D 62(1-4), 338–346.
- Shilnikov & Shilnikov  Shilnikov, A. & Shilnikov, L.  “On the nonsymmetric Lorenz model,” Inter. J. Bif. Chaos 1, 773–776.
- Shilnikov et al.  Shilnikov, A. L., Shilnikov, L. P. & Turaev, D. V.  “Normal forms and Lorenz attractors,” Inter. J. Bif. Chaos 3, 1123–1139.
- Shilnikov  Shilnikov, L.  ‘‘The existence of a denumerable set of periodic motions in four-dimensional space in an extended neighborhood of a saddle-focus.” Soviet Math. Dokl. 8(1), 54–58.
- Shilnikov  Shilnikov, L.  “On the birth of a periodic motion from a trajectory bi-asymptotic to an equilibrium state pf the saddle type.” Soviet Math. Sbornik. 35(3), 240–264.
- Shilnikov  Shilnikov, L.  “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  Shilnikov, L.  “The theory of bifurcations and quasiattractors,” Uspeh. Math. Nauk 36(4), 240–242.
- Shilnikov  Shilnikov, L.  “Chua’s circuit: Rigorous results and future problems,” Inter. J. Bif. Chaos 4(3), 489–518.
- Shilnikov  Shilnikov, L.  “Mathematical problems of nonlinear dynamics: A tutorial. Visions of nonlinear mechanics in the 21st century,” Journal of the Franklin Institute. 334(5-6), 793–864.
- Shilnikov  Shilnikov, L.  “Bifurcations and strange attractors,” Proc. Int. Congress of Mathematicians, Beijing (China) (Invited Lectures). 3, 349–372.
- Shilnikov & Shilnikov  Shilnikov, L. & Shilnikov, A.  “Shilnikov bifurcation,” Scholarpedia 2(8), 1891.
- Shilnikov  Shilnikov, L. P.  “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  Shimizu, T. & Morioka, N.  “On the bifurcation of a symmetric limit cycle to an asymmetric one in a simple model,” Physics Letters A 76, 201 – 204.
- Sparrow  Sparrow, C.  The Lorenz equations: bifurcations, chaos, and strange attractors, Applied Mathematical Sciences, Vol. 41 (Springer-Verlag, New York).
- Tresser & Williams  Tresser, C. & Williams, R.  “Splitting words and Lorenz braids,” Physica D: Nonlinear Phenomena 62(1–4), 15–21.
- Turaev & Shilnikov  Turaev, D. & Shilnikov, L.  “An example of a wild strange attractor,” Sbornik. Math. 189(2), 291–314.
- Vladimirov & Volkov  Vladimirov, A. & Volkov, D.  “Low-intensity chaotic operations of a laser with a saturable absorber,” Optics Communications 100(1-4), 351–360.
- Wieczorek & Krauskopf  Wieczorek, S. & Krauskopf, B.  “Bifurcations of -homoclinic orbits in optically injected lasers,” Nonlinearity 18, 1095–1120.