A route to computational chaos revisited: noninvertibility and the breakup of an invariant circle Christos E. Frouzakis 3 0 Combustion Research Laboratory, Paul Scherrer Institute 0 2 CH-5232, Villigen, Switzerland n a Ioannis G. Kevrekidis ∗ J 5 Department of Chemical Engineering 2 and Program in Applied and Computational Mathematics, ] Princeton University, Princeton, NJ 08544 S D Bruce B. Peckham . h t Department of Mathematics and Statistics a m University of Minnesota Duluth, Duluth, MN 55812 [ 1 February 1, 2008 v 1 0 Abstract 3 1 In a one-parameter study of a noninvertible family of maps of the plane aris- 0 ing in the context of a numerical integration scheme, Lorenz studied a sequence of 3 0 transitions from an attracting fixed point to “computational chaos.” As part of the h/ transition sequence, he proposed the following as a possible scenario for the breakup t ofaninvariantcircle: theinvariantcircledevelopsregionsofincreasingly sharpercur- a m vature until at a critical parameter value it develops cusps; beyond this parameter : value, the invariant circle fails to persist, and the system exhibits chaotic behavior v on an invariant set with loops [16]. We investigate this problem in more detail and i X show that the invariant circle is actually destroyed in a global bifurcation before it r has a chance to develop cusps. Instead, the global unstable manifolds of saddle-type a periodicpoints arethe objects which develop cuspsand subsequently “loops” or “an- tennae.” The one-parameter study is better understood when embedded in the full two-parameter space and viewed in the context of the two-parameter Arnold horn structure. Certain elements of the interplay of noninvertibility with this structure, theassociated invariant circles, periodicpoints andglobal bifurcations are examined. Keywords: Noninvertible maps, bifurcation, chaos, integration, invariant circles. ∗ Corresponding author, Department of Chemical Engineering, Princeton University, Princeton, NJ 08544,Phone: (609) 258 2818, Fax: (609) 258 0211, e-mail: [email protected] 1 1 Introduction In an insightful 1989 paper entitled “Computational Chaos: a prelude to computational instability” E. N. Lorenz [16] reported on a one-parameter computational study of the dynamics of the noninvertible map x = (1+aτ)x −τx y L : R2 → R2 = n+1 n n n (1) (a,τ) y = (1−τ)y +τx2, ( n+1 n n which arises from a simple forward Euler integration scheme (τ being the time step of the integration) of the two coupled nonlinear ODEs dx/dt = ax−xy dy/dt = −y +x2. (2) These ODEs are obtained by starting with the familiar Lorenz system [15]: x˙ = −σ(x − y),y˙ = −xz + ρx − y,z˙ = xy − βz, letting σ → ∞, and rescaling the variables and the remaining parameters. Lorenz fixed the value of a at 0.36, and varied the time step τ. Al- though the corresponding differential equations (eq. (2)) exhibit an attracting equilibrium point, computer simulations indicated that the discrete approximation (eq. (1)) progressed from exhibiting an attracting fixed point, to an attracting invariant circle (IC), to a chaotic attractor (termed computational chaos) as the time step τ was gradually increased. Se- quences of bifurcations similar to those described by Lorenz have been observed for other noninvertible maps of the plane [24, 31], suggesting that this might be a universal “non- invertible route to chaos.” Some of the bifurcations along the route — local changes of stability, homoclinic and heteroclinic tangencies, crises caused by an attractor interacting with its basin boundary — are also observed in families of invertible maps. New bifur- cations, however, unique to noninvertible families, as well as “invertible bifurcations with noninvertible complications,” are also observed along the route. In this paper we perform a more detailed numerical study of Lorenz’s family, but focus primarilyonanarrowrangeofparameterswhichincludesthebreakupoftheinvariantcircle and the identification of τ (the greatest lower bound on the values of the parameter CD τ for which the corresponding maps have “chaotic dynamics”), and τ , the greatest chaos lower bound on the values of the parameter τ for which the corresponding maps exhibit a “chaotic attractor.” Lorenz was interested in τ (labelled τ in [16]) as an indication chaos b of how poorly the Euler map L approximated the original differential equation. He (0.36,τ) observed that, in the transition from smooth IC to chaotic attractor, the process began with the circle developing features with increasingly high curvature. There appeared to be a critical τ value at which the “IC developed cusps.” After this value, the “chaotic attractor” suggested by computer simulations appeared to have loops and a Cantor-like structure. It certainly was no longer a topological circle. He claimed (correctly) that if this scenario did in fact occur, then this critical τ value had to be an upper bound on τ . chaos Our numerical investigations suggest that a smooth IC does not persist all the way to a “cusp” parameter. Instead, the IC is destroyed in a heteroclinic tangency initiating 2 the crossing of a branch of the stable manifold and a branch of the unstable manifold of a periodic (here period-37) saddle point. (As for invertible maps, chaotic dynamics are apparently present during this crossing, so this tangency is an upper bound on τ .) It CD is the unstable manifold that subsequently develops cusps at a critical parameter value we labelτ , andloopsbeyond that value. Asecond manifoldcrossing results intheapparent cusp appearance of the chaotic attractor, with loops inherited from the unstable manifold. We conclude thatτ isnotasingle isolatedbifurcationvalueseparating smoothICattractors cusp from chaotic attractors. Rather, τ is strictly above τ and strictly below τ . All three parameter values cusp CD chaos are part of the transition from a smooth IC attractor to a chaotic attractor. Because a cusp on an IC with an irrational rotation number would necessarily force a dense set of cusp points on the IC, it would make the existence of the IC itself unlikely. We therefore expect the existence of an entire transition interval of parameters, as opposed to a single bifurcation parameter, to be the generic scenario for the transition between smooth IC and chaotic attractor in the presence of noninvertibility. Mechanisms for the breakup of ICs are relatively well established for invertible maps of the plane (see for example [4, 6, 26]). In the invertible scenario which appears to us to most closely match the noninvertible scenario of this paper, an IC is born in a Hopf (also called Neimark-Sacker) bifurcation, grows in size, coexists with a periodic orbit (after a saddle-node birth of the periodic orbit off the IC), is destroyed in a first global manifold crossing (also referred to in the literature as a crisis [25] where the IC collides with its basin boundary), and is reconstituted after a second global manifold crossing. The periodic orbit which persists through the destruction and reconstitution of the IC typically switches fromoutside (inside) the ICbeforeit isdestroyed toinside (outside) theIC after it isreconstituted. Chaoticinvariant sets necessarily exist only during themanifoldcrossings. In contrast, the breakup of the IC in the noninvertible scenario is part of a transition to a chaotic attractor. In particular, the “chaotic attractor with loops,” which appears in this noninvertible case, is not a feature which appears in invertible maps. In fact, cusps and loops on iterates of smoothly embedded curves are not possible for smooth invertible maps, but are common features on (invariant) curves of noninvertible maps [11]. Further, bifurcations such as the manifold crossings (crises) can have additional complications due to the presence of noninvertibility. For example, the unstable “manifolds” involved in the global crossings may have self intersections and cusps. Or stable and unstable “manifolds” may cross transversely at one homoclinic point, but fail to preserve transversality at other points along the homoclinic orbit. In any case, the transition mechanism for the maps in eq. (1) is different from the invertible transitions, and the noninvertible nature of the map plays an important role: the IC — for parameter values for which it exists — interacts with (actually crosses) the locus on which the Jacobian of the linearized map becomes singular. This locus (and its images and preimages) is crucial in organizing the dynamics of noninvertible maps; it is termed “critical curve” in the pertinent literature [14, 18, 24, 1] and constitutes the generalization (in two dimensions) of the critical point in unimodal maps of the interval [8,7]. Inparticular, asexplained in[11], theangleofintersection between asaddle unstable 3 manifold and the critical curve underpins the transition of the local image of the manifold from being smooth and injective to nonsmooth and injective (with a cusp) to smooth but not injective (with loops). As for invertible Hopf bifurcations (even in the noninvertible setting the dynamics near a Hopf bifurcation are necessarily locally invertible), an understanding of Lorenz’s a = 0.36 one-parameter family is possible only by viewing it in the two-parameter Arnold horncontext [3,6,4,27,28]. We dothis withL , using thesecond parameter (a)already (a,τ) provided in eq. (1). This leads to an examination of the internal Arnold resonance horn structure, to be compared to and contrasted with the invertible case. The interiors of the resonance horns have features which vary significantly from those in the horns of invertible maps (see preliminary results in [13, 12]). Although our understanding of this internal structure is far from complete, examination of the partial picture still provides insight into the bifurcations observed in the Lorenz one-parameter cut. Developing computational tools to further investigate the internal structure is part of our ongoing research. In this paper the observations of the transition in [16] are briefly summarized in Sec. 2, and then revisited in detail in Sec. 3. In Sec. 4, we place our one-parameter cut in con- text by examining the “resonance horns” in the larger two-parameter space. We discuss related noninvertible issues, including computational challenges, in Sec. 5, and state final conclusions in Sec. 6. 2 Observations Figure 1 briefly summarizes the relevant observations in [16]. For a fixed at 0.36, the differential equation (eq. (2)) has an attracting equilibrium point at (x,y) = (0.6,0.36). Its basin of attraction is the right half plane. There is a symmetric attracting equilibrium for the left half plane at (−0.6,0.36). We will deal only with the right half plane attractors in this paper. The corresponding maps L all have a fixed point at (x,y) = (0.6,0.36). (0.36,τ) In a one-parameter cut with respect to τ, the fixed point is stable for small τ, but loses its stability via Hopf bifurcation at τ = 1.38889. The resulting attractor is initially a smooth IC (fig. 1a, τ = 1.55) but then develops progressively sharper features (fig. 1b, τ = 1.775, where the circle is still smooth even though it appears to have cusps), then self-intersects, suggesting chaos (fig. 1c, τ = 1.785), and eventually exhibits clear signs of chaotic behavior (fig. 1d, τ = 1.91). Lorenz supports his claim of chaos with the computation of a positive Lyapunov exponent forτ = 1.91 [16]. The last two figures are unusual in that the attractor apparently shows self-intersections, a clear sign of noninvertibility. All four figures were created by following a single orbit, after dropping an initial transient, from an initial condition near the fixed point. We use the term attractor loosely to denote the object generated by the computer via such a simulation. More formal definitions appear in the next section. Lorenz proposed the following scenario as the simplest possible transition of the at- tractor from a smooth IC to a chaotic attractor. As a parameter is increased (τ in this example) and the IC crosses successively farther over the critical curve (called “J ” below), 0 4 Figure 1: Phase portraits of the map L for a = 0.36 and (a) τ = 1.55, smooth IC (b) τ = 1.775, regions of sharp curvature appear on the IC (c) τ = 1.785, self-intersections appear (d) τ = 1.91, “fully chaotic” attractor causing it to develop successively sharper features. Simultaneously, its rotation number changes (increasing with τ in this case). By restricting to τ values for which the IC is quasiperiodic (and avoiding periodic lockings, defined below in Sec. 3.1), one obtains a sequence of τ values, each corresponding to a map with a smooth quasiperiodic IC, and limiting to τ , the lower limit on chaotic behavior. At τ , the attractor, which may chaos chaos or may not still be a topological circle, would develop cusps. Beyond τ , the attractor chaos would be chaotic, or at least have parameter values accumulating to τ from above for chaos which the corresponding maps would exhibit chaos. We will now examine the implications of such a scenario. 3 Observations Revisited By performing a more detailed numerical study, we show below that, while Lorenz’s sce- nario is the correct basic description of the transition, the IC actually breaks apart before the “cusp parameter” value is reached. The transition from smooth IC to chaotic attractor thus occurs over a range of parameter values, rather than at a single critical parameter value. We emphasize that Lorenz did not claim that his scenario did happen, only that there could be no simpler scenario. Our study shows that the proposed “simple” scenario doesnothappeninthisparticularfamily, norshoulditbeexpected inthegeneraltransition from smooth invariant circle to chaotic attractor in the presence of noninvertibility. 5 Before beginning the description of our numerical investigation, we recall some termi- nology and preliminary results. 3.1 Noninvertible preliminaries Let f : ℜ2 → ℜ2 be a smooth (C∞) map. Let x be a period-q point for f. As for invertible maps, we define the stable manifold of x to be Ws(x) = {y ∈ ℜ2 : fkq(y) → x as k → ∞}. Since the inverse map is not necessarily uniquely defined, we define the unstable manifold of x as Wu(x) = {y ∈ ℜ2 : there exists a biinfinite orbit {y } with y = f(y ),y = j i+1 i 0 y and y−kq → x as k → ∞}. The use of the term “manifold” is an abuse of terminology for noninvertible maps since, while the local stable and unstable manifolds are guaranteed to be true manifolds, the global manifolds are not [20, 30, 11]. (Some authors address this problem by using “set” instead of “manifold” [11, 17, 19, 21, 32].) By iteration of the local smooth unstable manifold, the global unstable manifold can be computed as for invertible maps. If the manifold (or any curve) ever crosses the critical set (J , defined 0 below) tangent to the zero eigenvector of the map at the crossing point, its image can have a cusp [11, 33]. More degenerate cases are also possible, but not considered here. The global unstable manifold can also have self intersections. Thus, it is globally a smoothly immersed submanifold, with smoothness violated only at forward images of these special crossingpoints. Thisisdiscussedfurtherbelow. Theglobalstablemanifold,duetomultiple preimages, can be disconnected, fail to be smooth, and even increase in dimension [33, 11]. The multiple preimages lead to interesting basins of attraction (see Section 5 and [2, 23]). A set S ⊂ ℜ2 is invariant if f(S) = S. When S is invariant and a topological circle we call it an invariant circle (IC). An invariant circle on which there exists a periodic orbit is called a frequency locked circle or circle in resonance. (The terminology is borrowed from return maps of periodically forced systems.) Typically we see a single attracting and a single repelling periodic orbit on a frequency locked circle, but multiple orbits are certainly possible. A periodic orbit is also called a periodic locking or a locked solution, although we usually reserve the use of periodic locking to refer to a periodic orbit on an invariant circle. We have already indicated above that we are using the term attractor informally to indicate whatever is suggested via computer simulations. More formally, we define a com- pact set K to be an attractor block if f(K) ⊂ K◦ (the interior of K). Then we define Λ, the attracting set associated with attractor block K by Λ = ∞ fn(K) [4, 30]. In the n=1 cases we consider, including for τ in the range depicted in fig. 1, we take the attractor T block K to be an annulus containing the “attractor” in each picture: inner radius a small circle around the repelling fixed point O, and outer radius a (topological) circle which is big enough to contain the attractor, but still in the right half plane. The attracting set will always be taken to be associated with this annulus. Any invariant circles or periodic points in the attractor block will also be contained in the attracting set Λ. So will the unstable manifolds of the periodic points. The attractors displayed by computer simulation are generally (an approximation of) some subset of Λ. If Λ exhibits at least some approximate recurrence property, such as a dense orbit (for example, a quasiperiodic IC), or a numerical orbit which appears dense (for example, a 6 Figure 2: (a) Invariant circle, Γ, one of its preimages, Γ−1, and the critical curves J and 2 0 J (τ = 1.775). The image of a “curve” crossing J (at A and B) is tangent to J (at 1 0 1 L(A) and L(B), respectively). (b) “weak chaotic ring” possessing loops (τ = 1.785); (c) enlargement showing the loop on the attractor at L(C). locked circle with a highperiodlocking), then theattractor might bea goodapproximation of all of Λ. In other cases, such as a locked circle with a low period locking, simulation would onlydisplay thesinks. Computationoftheunstable manifolds ofthesaddles (crucial to our numerical investigation) often gives a better picture of Λ than does a simulation. Note that a locked circle is still an invariant circle, even though the whole circle is not the attractor. The circle is smooth by the unstable manifold theorem [20] through the saddles, but as with invertible maps, its smoothness where two branches of the unstable manifold meet at the sinks depends on the eigenvalues [4]. Unlike invertible maps, as indicated above, the smoothness of an IC can also be violated by the development of cusps. We now define the critical curve J0 (termed “curve C” in [16] and often LC−1 for Ligne Critique in the recent noninvertible literature [11, 24]) as the locus in phase space where the linearization of the map becomes singular. We call its image J (“curve D” and LC, 1 respectively). The two sets J and J , as well as additional images and preimages, are 0 1 known to be a key to understanding the dynamics of a noninvertible map. Consider, for example, Fig. 2a showing an enlargement of the IC, Γ for τ = 1.775, with three important additions: J , J , and Γ−1, one additional first-rank preimage of the invariant circle Γ. 0 1 2 The map L (abbreviated as L) has either one or three inverses depending on the (0.36,1.775) phase-plane point in consideration (the term “Z −Z map” has been proposed for such 1 3 maps, [23]). The invariant circle Γ has three first-rank preimages: itself, the curve Γ−1 2 shown in the figure, and a third preimage Γ−1, further away in phase space (not shown). 3 The geometry of the map can be visualized by first folding the left side of the phase plane along J and onto the right side. Γ and Γ−1 should now coincide. Next rotate roughly 0 2 ◦ 90 so that J maps to J . The image of Γ lands exactly on itself (before the folding), 0 1 sending A to L(A) and B to L(B) in the process. Note that L maps points that are in 7 the intersection of the two regions bounded by Γ and Γ−1 to the region bounded by the 2 pieces of Γ and J which connect at L(A) and L(B). That is, some points inside Γ are 1 mapped outside Γ. Similarly, points inside Γ−1, but not Γ are mapped from outside Γ to 2 inside Γ. Contrast this with the property that interiors of invariant circles for (orientation preserving) invertible maps are invariant. Note also that Γ is tangent to J at L(A) and L(B). The “folding” along J , typ- 1 0 ically creates a quadratic tangency along J at images of curves crossing J , unless the 1 0 curve crossing J does so at an angle tangent to the zero eigenvector at that point. This 0 phenomenon has been extensively discussed in [11]. Our explanation here is based on the understanding reached in that paper. As a parameter (in this case τ) is varied, the map L (again abbreviated as L), the entire attracting set Λ, the curves J and J as well (0.36,τ) 0 1 as the intersection points of Λ with J (e.g. A) and with J (e.g. L(A)) vary. In our case, 0 1 between τ values for figs. 2a and b, the invariant curve develops a “fjord’ which pushes across J , creating two new intersection points of Λ with J , and two new tangencies with 0 0 J at the two image points. We call the lower intersection point of the fjord C (see fig. 2b). 1 The points C and L(C) vary with τ as well. When the tangent of the “curve” in Λ at the intersection point C becomes coincident with the null vector of the Jacobian of the map at C, the image of Λ (Λ itself) acquires a cusp touching J at L(C) [11]. As the parameter 1 value τ is further varied, the cusp becomes again a quadratic tangency, but Λ acquires a loop (figs. 2b and c). Iteration of this “loop” gives rise to an infinity of such loops (called “antennae” in [16]) on Λ, suggesting chaotic behavior. Further, Lorenz gives a heuristic justification of a simple criterion for an attractor (not formally defined, but assumed by us to be a compact invariant set with some sort of recur- rence) to be a chaotic attractor (an attractor exhibiting sensitivity to initial conditions): the existence of two distinct points on the attracator which map to the same image point (also on the attractor) [16]. In this context, the tracing out of an object with a single orbit suggests recurrence, and existence of a loop implies the existence of two points which map to the same image point. Thus, if this criterion is correct, attractors with loops would necessarily be chaotic attractors. 3.2 Numerical investigation We set out to examine the transition from smooth IC to chaotic attractor in finer detail in figures 3 to 11. Before we embark on a detailed description, we should state that the real picture is yet more complex, involving similar phenomena on an even finer scale. For example, the same phenomena which cause the destruction of the “large” IC could also destroy the “small” IC’s which we refer to below. Nevertheless, this sequence of figures includes the essential features of the transition: coexistence of a periodic orbit with a nearby IC, destruction of the IC through a global bifurcation, appearance of loops on an unstable manifold, and the reappearance of an attractor, this time chaotic with loops. It might be useful to look ahead to the two-parameter bifurcation diagram in fig. 12 to put the one-parameter cut described in the current section in context. We start from the smooth IC at τ = 1.775 (fig. 2a). At τ ≈ 1.776243 a period-37 8 Figure 3: Coexisting attractors for τ = 1.7765: the attracting period-37 foci (•), and the IC. The period-37 saddles (×) are also shown, along with both branches of the unstable manifold S . One branch approaches the focus F , while the other branch approaches the 5 5 IC. saddle and node pair appears in a saddle-node (turning point) bifurcation. The bifurcation occurs off of the IC, and results in the coexistence of two attractors: the stable invariant circle, Γ, and an attracting period-37 sink. Although the sinks have real eigenvalues at the saddle-node bifurcation, the eigenvalues very quickly become complex and they have been labelled in fig. 3b for τ = 1.7765 by F for “focus.” The separatrix between these i two coexisting attractors consists of the stable manifolds of the period-37 saddles. Global unstable manifolds have been systematically constructed as (37th) forward iterates of local unstable manifolds as for invertible maps. In this case, we have plotted only the two branches emanating from saddle S . The right branch of the global unstable manifold 5 directly approaches F . This is clear in the blowup of fig. 3b. The left branch very quickly 5 approaches the IC Γ. Successive 37th iterates travel slowly in a clockwise direction around the IC, eventually showing a good approximation to the whole IC. The approach to the IC is somewhat clearer in fig. 3b, where the branch appears to “join” the IC just to the left of S . The unstable manifold then continues almost “on” the IC: past F , S , and 5 28 28 F , evolving around the IC, reentering from the right, just below F , before completing 14 5 the (approximation to the) IC just to the left of S . 5 Notice also that 18 of the period-37 saddle-node pairs are “inside” (those marked by numbers inside Γ) and the rest “outside” the invariant circle; this is another feature of noninvertible maps [11]). We discuss this again in Sec. 5. As τ increases further, the consequences of global bifurcations are observed; these bi- 9 Figure 4: During the “first” manifold crossing (τ = 1.77688). a) the period-37 saddle and sink solutions have been plotted with one branch of the global unstable manifold (the branch involved in the crossing) which emanates from each of the 37 saddles. b) Both branches of the unstable manifold emanating from the saddle labelled S have been 19 plotted (but not the branches emanating from S , or any other saddle). 5 furcations involve interactions between the stable and unstable manifolds of the saddle period-37 points. The first such global bifurcation is evident in the sequence of figures 3 through 5: the left branch of the unstable manifold of S misses S to the left in fig. 3 and 19 5 approaches the IC (the curve coming into the inset from the right and below F is prac- 5 tically coincident with this unstable manifold branch); it makes the “saddle connection” with S in fig. 4; it misses S to the right in fig. 5 and approaches F . 5 5 5 The global crossing actually “starts” with a heteroclinic tangency (between a branch of a stable manifold of a period-37 saddle and and branch of an unstable manifold of another saddle on the same orbit) which destroys the IC at τ ≈ 1.776878. This tangency is also called a crisis [25], since as τ increases toward the tangency parameter value, the IC (approached by the unstable manifold) and its basin boundary (the stable manifold) approach each other, colliding at the tangency parameter value. The attracting set during the crossing appears — as for invertible maps — to be a complicated topological object which at least includes the closure of the unstable manifolds of the period-37 saddles. Numerically, most orbits appear to approach the period-37 sink (fig. 4). Thirty-seventh iteratesofeachofthe37branchesoftheunstablemanifoldplottedinfig.4aevolveclockwise toward the “next” saddle, suggesting that we are at least close to a heteroclinic crossing of stable and unstable manifolds. Although the stable manifolds of the saddles cannot be realistically computed because oftheexponential explosion (“arborescence”) ofthe number of period-37 preimages, the crossing of the stable manifold is evident in the inset (fig. 4b). 10