Intermingled basins in coupled Lorenz systems Sabrina Camargo 1, Ricardo L. Viana 2, and Celia Anteneodo 1,3 1 Department of Physics, PUC-Rio, Rio de Janeiro, Brazil 2 Department of Physics, Federal University of Parana´, Curitiba, Brazil 3 National Institute of Science and Technology for Complex Systems, Rio de Janeiro, Brazil. (Dated: January 31, 2012) We consider a system of two identical linearly coupled Lorenz oscillators, presenting synchronization of chaotic motion for a specified range of the coupling strength. We verify the existence ofglobalsynchronizationandantisynchronizationattractorswithintermingledbasinsofattraction, suchthatthebasinofoneattractorisriddledwithholesbelongingtothebasinoftheotherattractor and vice versa. We investigated this phenomenon by verifying the fulfillment of the mathematical requirementsforintermingledbasins,andalsoobtainedscalinglawsthatcharacterizequantitatively 2 theriddling of both basins in thissystem. 1 0 PACSnumbers: 0.45.Xt,05.45.Df,05.45.Pq,05.45.-a 2 n a I. INTRODUCTION attractor. Thissituationcanstillbedramaticallyworsen J when the riddling phenomenon occurs. 0 The Lorenz system As a matter of fact, multistable dynamical systems 3 typically have a very complicated structure of basins ] x˙ =α(y−x), y˙ =βx−y−xz, z˙ =−γz+xy, (1) of attraction, that may be delimited by fractal bound- D aries[7]. Suppose, for instance,thata dynamicalsystem C for α = 10, β = 28, and γ = 8/3, displays a chaotic has two attractors, with the corresponding basins of at- . attractor with the familiar butterfly-like shape [1] . It traction sharing a common basin boundary in the phase n is often quoted as a paradigmatic system in nonlinear space. If a ball centered at a given initial condition and i l dynamics, since it displays many interesting dynamical witharadiusequaltotheuncertaintylevelinterceptsthe n [ properties of chaotic dissipative systems. Moreover its basin boundary, we cannot say a priori which attractor equations mimetize the dynamical behavior expected to the system will asymptote to [8]. If that boundary is 2 occur in some physically relevant systems, as convection a curve, even if fractal, the final-state sensitivity prob- v rolls in the atmosphere [1], single-mode lasers [2], and lem can be circumvented by decreasing the radius of the 1 segmented disk dynamos [3]. Coupled Lorenz systems uncertaintyball(this canbe doneinexperimentalornu- 8 5 could arise as well in the mathematical modeling of re- mericalsettings,by increasingthe precisionindetermin- 5 lated physical problems. The simplest case in the latter ing the initial condition in phase space). However, such . categoryis the coupling of two identical Lorenz systems. reductionofuncertaininitialconditionsisnotpossiblein 1 1 In identical coupled systems, even if chaotic, synchro- thelimitcaseinwhichthefractalboundaryisarea-filling, 1 nization of trajectories may occur [4]. This phenomenon i.e., the (box-counting)dimensionof the basinboundary 1 has been studied for more than two decades, motivating gets close to the dimension of the phase space itself [9]. v: a wealth of analytical, numerical, and even experimen- In that limit case, the fraction of uncertain initial con- i tal results [5]. Synchronization of chaos, besides its own ditions will likely not decrease no matter how much we X interest as a mathematical problem, finds applications decrease the uncertainty balls of each initial condition. r for instance in secure communications [6]. The chaotic The latter situation occurs for riddled basins [9]. a nature of the dynamics of one of the systems can be ex- From the mathematical point of view, riddled basins ploited to code messageswhichcould be sentto aniden- areobservedindynamicalsystemsthatexhibitaninvari- ticalsystemthroughsome formofcoupling. If the latter ant smooth hypersurface with a chaotic attractor lying systemis synchronizedwith the former,the messagecan onit,anotherasymptoticfinalstate,outoftheinvariant be securely uncoded. subspace,andnegativeLyapunovexponenttransverseto For two completely synchronized systems, either peri- the invariant subspace with positive finite-time fluctua- odicorchaotic,theirdynamicalvariablesareequalforall tions [10–12]. Under the conditions above, riddling orig- times. On the other hand, if instead of the difference, it inates from the loss of transversal stability of unstable isthe sumofsomeoftheir dynamicalvariablesthatvan- periodic orbits embedded in the chaotic attractor [13], ishes, the two systems are said to antisynchronize. Due despite the attractor being transversely stable in aver- to phase-space symmetries, coupled Lorenz systems can age. In this context, attractors must be understood in exhibit both synchronized and antisynchronized states. the weak sense of Milnor [14]. The transition associated Then,forsecurecommunicationspurposes,theexistence to the first unstable orbit on the attractor that losses ofanother,antisynchronized,stateisinprincipleasource transversal stability determines the riddling bifurcation of troubles since, depending on the initial condition, the (see for instance Ref. [15] for an overview). Depending receiver system could be tuned to the antisynchronized on the way these orbits loss stability and even on the 2 dynamics outside the invariant manifold, different bifur- II. COUPLED LORENZ SYSTEMS cationscenariosanddifferentformsofriddledbasinscan occur [13, 16–18] (to cite a few examples). Many different coupling schemes are possible for two identical Lorenz systems [24]. We have chosen, for If riddled basins exist in a multistable chaotic system, symmetry reasons, a diffusive coupling through the z- their final states are utterly unpredictable, i.e. we can- variable, as follows not say - with any degree of certainty - which attractor the system will evolve to for long times [19]. The situa- x˙ = α(y −x ), tion,inthis case,is similartothatforarandomprocess, 1 1 1 y˙ = βx −y −x z , for which there can only be determined a probability for 1 1 1 1 1 z˙ = −γz +x y +ε(z −z ), predicting the final state of the system. In fact, some 1 1 1 1 2 1 (2) phenomena formerly attributed to random variations in x˙2 = α(y2−x2), initialconditionscanbealsointerpretedasaconsequence y˙ = βx −y −x z , 2 2 2 2 2 of riddling [20]. z˙ = −γz +x y +ε(z −z ), 2 2 2 2 1 2 Thesimplestcaseofriddlingiswhenonlyoneoftheco- where we will use the same values for α, β, and γ, as in existing attractors have a riddled basin. However, when the uncoupled case, and ε is the coupling strength. there is more than one invariant subspace, then more On considering the dynamicalbehavior of the coupled thanoneattractorcanberiddled. Inthiscase,the basin system, it is convenient to perform the changes of vari- structure is called intermingled [9]. ables The aim of the present work is precisely to show the (x −x ) (y −y ) (z −z ) existenceofintermingledbasinsofattractionforthesyn- x= 2 1 , y = 2 1 , z = 2 1 , chronized and antisynchronized states of two coupled 2 2 2 (3) (x +x ) (y +y ) (z +z ) 2 1 2 1 2 1 Lorenz oscillators. In previous literature there are al- X = , Y = , Z = , 2 2 2 ready clues of such phenomenon. Kim and coworkers [21], in a work about anti-synchronization of coupled after which the coupled system (2) becomes chaotic oscillators, point to the possibility of a riddled basin of synchronization in coupled Lorenz systems, but x˙ = α(y−x), withoutgoingfurtheronthatissue. Furthermore,aone- y˙ = βx−y−(Xz+Zx), dimensional reduction of the Lorenz system (to a piece- z˙ = −(γ+2ε)z+Xy+Yx, wise approximation to the well-known Lorenz map) was X˙ = α(Y −X), (4) low-dimensional enough for an analytical treatment to Y˙ = βX −Y −(XZ+xz), be feasible and show the riddling of the synchroniza- Z˙ = −γZ+XY +xy. tion basin [22]. The verification of the transversal stability conditions through direct methods (i.e., by mak- Whenever more convenient to the analysis, we will refer ing a linear stability analysis of eachinvariantsubspace) either to the new or the old variables. is quite difficult in two coupled Lorenz systems, since FrominspectingEqs.(4)therefollowsthatthedynam- the phase space is six-dimensional. Then, we investigate ics of the coupled system is invariantwith respect to the those propertiesnumerically. We alsocharacterizequan- transformation (x,y,z)→(−x,−y,−z). Hence the con- titatively the riddledbasinsby meansofthe scalinglaws ditions x=y =z =0 define an invariant subspace M : givingtheprobabilityofmakingwrongpredictionsonthe s oneinitialconditionthatbelongsto this subspacegener- final state of the system, with respect to two quantities ates a trajectory in phase space that remains in M for of interest: (i) the phase-space distance to the invariant s any time. This three-dimensional subspace defines the subspace; and (ii) the uncertainty radius for each initial complete(orglobal)synchronizationmanifoldcharacter- condition [23]. We have verified that, for both quanti- ized by x =x , y =y , z =z . ties, the probability scales as a power-law, as required 1 2 1 2 1 2 ThedynamicsintheinvariantsubspaceM ,described for riddled basins. s by the variables (X,Y,Z), is governed by the equations Therestofthepaperisorganizedasfollows: SectionII of the uncoupled Lorenzsystem, hence there is a chaotic describesthecoupledsystemofLorenzoscillators,aswell attractor As (butterfly-like shape) lying in Ms. as the existence of both synchronized and antisynchro- Analogously, due to the symmetry (X,Y,z) → nized states. Section III presents a preliminary discus- (−X,−Y,−z),thestatesforwhichX =Y =z =0define sion of the basins of attraction of both the synchronized another invariant subspace M (anti-synchronization a and antisynchronized states. The mathematical proper- manifold), associated to the attractor A , in which a tiesrequiredforriddledbasinsandthenecessarytoolsare (x,y,Z) follows the dynamics of the uncoupled system, the object of Section IV. Section V discusses the quanti- i.e. A is a Lorenz chaotic attractor in M . a a tative characterization of riddled basins through scaling Therearealsoothersymmetriesalreadypresentinthe laws and the theoretical results supporting them. The uncoupled Lorenz system. Notice in Eqs. (2) that either last Section contains our conclusions and final remarks. (x ,y ) → (−x ,−y ) or (x ,y ) → (−x ,−y ) lead, to 1 1 1 1 2 2 2 2 3 four-dimensionalinvariantsubspaces,whilebothsymme- randomly chosen in the interval [20,24) according to a tries together lead to a two-dimensional invariant sub- uniform probability distribution. spacewithasaddlepointatthe origin. Finally,included For instance, for a coupling strength ε = 1.0 [Fig. in this two-dimensional subspace, the lines at z = z 1(a)], the section of the basin of attractor A is a se- 1 2 a and z = −z also represent invariant subsets. We did ries of thin filaments stemming from the diagonal. The 1 2 not find any other relevant attractor other than A and filaments are non-uniformly distributed and have a sug- s A , which are attractors of the dynamics in the respec- gestive self-similar appearance. In fact, successive mag- a tive subspaces M and M , and can become attractors nifications of [Fig. 1(a)] reveal similar patterns (see Ref. s a forthe wholephasespacedepending ontheir transversal [21]). Suchscenarioisalsoobservedforothervaluesofε, stability. asillustratedinFig. 1(b-d),evenifsomefeatureschange withε,suchastherelativeareaofeachbasin,orthedef- inition of the tongues anchored in the diagonal. Let us III. BASINS OF ATTRACTION note that other cuts also display a tongue structure, as depicted in Fig. 2 for ε=2.0. In dynamical systems with more than one attractor, the corresponding basins may have fractal boundaries and even more complicated structures like the Wada property[25]. Accordingly,inthecoupledLorenzsystem (2), the two coexisting attractors representing synchronized and antisynchronized states are expected to have such complex basin boundary structure. Since the phase space of the coupled system is six- dimensional, the visualizationof the basins of attraction depends on convenient phase space sections or projec- tions. Figure 1 shows a section of the basin of the anti- FIG. 2: Section at z1 = z2 = 22.0 and (a) y1 = y2 = 1.0, synchronization(synchronization)attractorAa (As), for x1,x2 random in [−1,3), (b) x1 = x2 = 1.0, y1,y2 random different values of the coupling parameter. in [−1,3) of the basins of synchronization (white pixels) and Each initial condition was integrated using a fourth- antisynchronization (black pixels) attractors of the coupled orderRunge-Kuttaschemewithfixedtimestep10−3 and Lorenz system, for ε=2.0. for a time t = 103, after which we determined to what attractorthecorrespondingorbithasasymptoted[26]. If The structure of the basins of attraction is indeed ex- anorbithasasymptotedtoanantisynchronized(synchro- pected to be altered by the coupling strength. As an nized)stateinM (M ),itsinitialconditionwaspainted example, in Fig. 3 we show that, for a given initial con- a s black (white). Hence the area painted black (white) is a dition (x1 =y1 =z1 =1.0 and x2 =y2 =z2 =0.5) inte- numerical approximation of a section of the basin of at- grateduptotime 102, thetrajectoriesinthe subspaceof tractor A (A ). We considered 105 initial conditions each oscillator are distinct for ε = 0.5, while for ε = 1.0 a s with x = x = y = y = 1.0 while z and z were trajectories tend to coincide due to synchronization. In 1 2 1 2 1 2 thelattercase,theoverlappingsegmentsreproduceacut ofthefamiliarattractorofthesingleLorenzsystem,since 24 (a) 24 (c) for synchronized orbits, the evolution proceeds towards the attractor in M which is defined by x = y = z = 0, s 23 23 and X = x = x , Y = y = y , Z = z = z follow the 1 2 1 2 1 2 z2 22 z2 22 dynamics of an uncoupled system, as described above. Differently, in the former case (ε = 0.5), the trajecto- 21 21 ries of each system depart from those of the uncoupled 2204 (b) 2204 (d) system. Moreover, the observation of synchronized or antisyn- 23 23 chronized states depends on the coupling strength. Re- z2 22 z2 22 callthattheexistenceofMs (i.e. thesynchronizedstate being a possible solution of the coupled equations) does 21 21 not mean necessarily that synchronized states, and in 20 20 particular states in its attractor A , can be observed 20 21 22 23 24 20 21 22 23 24 s z1 z1 in numerical simulations. This occurs only if there is transversal stability, in the sense that any infinitesimal FIG. 1: Section at x1 = x2 = y1 = y2 = 1.0 of the basins of displacement along directions transversal to Ms decays synchronization(whitepixels)andantisynchronization(black exponentially with time. pixels) attractors of the coupled Lorenz system, for ε = (a) Let us remark that, due to the symmetry of the 1.0, (b) 2.0, (c) 2.5 and (d) 2.8. equationswithrespecttosynchronized/antisynchronized 4 tained for the sums x +x , y +y , indicating that the 1 2 1 2 basins of the synchronized and antisynchronized states are complementary to each other. 1.00 (a) s n o diti 0.75 n o c al niti 0.50 of i n o cti 0.25 fs+fa Fra fs f a 1-(f+f) s a 0.00 0 100 200 300 400 500 600 700 800 900 1000 t 1.00 FIG. 3: Trajectories of the coupled Lorenz system for the same initial condition (x1 = y1 = z1 = 1.0 and x2 = y2 = (b) z2 = 0.5) up to t = 100 and different coupling values: (a) e=0.50 ε = 0.5; (b) 1.0. In each case, the time evolution of the 0.75 e=0.75 differences of coordinates are also shown (c)-(d) for ε = 0.5 e=0.80 and (e)-(f) for ε=1.0. +f)a ee==11..0500 1-(fs 0.50 ee==22..0500 states, comments for attractor A are also valid for A . e=2.80 s a e=3.20 Inordertovisualizetheexistenceofatransverselysta- 0.25 blesynchronizationmanifold,weconsiderthe differences x − x , y − y , and z − z , which must vanish if a 1 2 1 2 1 2 synchronized attractor is achieved. For ε & 0.7, z −z 0.00 1 2 0 500 1000 1500 2000 vanishes [Fig. 4(c)], while the other two differences may t also vanish (global synchronization) or not (local synchronization) [Fig. 4(a) and (b)]. Similar plots are ob- FIG.5: (a) Fraction ofinitial conditions(overatotal of104) yieldingtrajectories asymptotingsynchronizedfs orantisyn- chronized fa states as a function of time, for ε = 2.0. We 40 (a) alsoindicatedthefractionoftrajectoriesreachingeitherstate 20 (fs+fa)ornoneofthem(1−fs+fa). (b)Fractionofinitial x-x12-2 00 ccoonudpiltiniognsstnreontgrteha.chingthesestatesfordifferentvaluesofthe -40 60 (b) 40 We did not find any relevant attractor for the coupled y-y12- 22 000 scyosntseimderoatthioenrsthatanthAesenadndofASae.ct.BIeIs,idwees ptheerfosrymmemdetthrye -40 following numerical experiment: we considered the ini- -60 40 (c) tial conditions used to plot the sections in Fig. 1 and, 20 for each time t we computed the fraction of initial con- z-z12 0 ditions that go either to As or Aa [Fig. 5(a)]. The sum -20 of these fractions rapidly approaches 100% [Fig. 5(a)], -40 meaning that the fraction of initial conditions that do 0 1 2 3 4 5 e notasymptotetothemgoestozero(filledsquaresinFig. 5(a)), suggesting the existence of only two attractors for FIG. 4: Difference between the coordinates (a) x1−x2; (b) the coupled system. This conclusion has been observed y1 −y2; (c) z1 −z2 of two coupled Lorenz systems at time to hold for ε & 0.7 as illustrated in [Fig. 5(b)]. Oth- t=103,as a function of thecoupling strength. Onehundred erwise, neither synchronized nor antisynchronized states initialconditionswererandomlychosen(asinFig.1)foreach areapproached,asillustratedinFig.3(c)-(d)forε=0.5. valueof ε (varied in steps of 0.01). (Therefore, basin diagrams as those shown in Fig. 1 will 5 be left blank). relevant since both attractors have the same Lyapunov spectrum. Thesecondmethodweused,andwhichcanbeapplied IV. RIDDLED AND INTERMINGLED BASINS toobtainonlythelargesttransversalexponent,istocon- sider the time evolution of an infinitesimal displacement The standard requirements for the existence of a rid- along a direction transversal to the synchronized sub- dledbasinaretheexistenceof(i)asmoothinvariantsub- space Ms, which is given by [7], space(oflowerdimensionthanthephase-space)contain- ingachaoticattractor,(ii)anotherasymptoticfinalstate 1 δ(t) λ = lim ln , (5) ⊥ (not necessarily chaotic) out of the invariant subspace, t→∞ t δ(0) (iii) negativity of the Lyapunov exponents transverse to the invariantsubspace with (iv) positive finite-time fluc- where δ(t) = p(δx)2+(δy)2+(δz)2 is the norm of the tuations [9–11], which are associated to the transversal transversedisplacement, whose evolution is given by the stability properties of unstable periodic orbits (UPOs) variational equations for (x,y,z), setting x=y =z =0, embedded in the attractor. For two symmetrically inter- i.e., mingledbasins,therequirementsformutualriddlingcan be summarized as follows: δ˙x = α(δy−δx), δ˙y = βδx−δy−Xδz−Zδx, 1. There are invariantmanifolds Ss and Sa contained δ˙z = −(γ+2ε)δz+Xδy+Yδx, in the phase space H. (6) X˙ = α(Y −X), 2. The dynamics on each manifold Ss,a has a chaotic Y˙ = βX −Y −XZ, attractor. Z˙ = −γZ+XY . 3. S are transversely stable, meaning that the s,a Figure 6 shows (in gray symbols) the three largest largesttransversalLyapunov exponent λ is nega- ⊥ (infinite-time) Lyapunov exponents as a function of the tive. coupling strength ε, obtained by means of the algorithm 4. Althoughweakstabilityholdsinaverage(condition byWolfet al.,andthelargesttransversalexponentgiven 3), UPOs embedded in the chaotic attractor are by (5) is indicated by a thick black line. One of the ex- transversely destabilized. ponents is always zero, corresponding to displacements along the trajectory. The largest exponent is practically In Section III we showed that conditions 1 and 2 are alwaysequalto ∼0.9andcorrespondstothe chaoticdy- fulfilled for the coupled Lorenz systems. There exists namics on A (A ). The third exponent is the largest s a two (three-dimensional) manifolds, Ss,a = Ms,a, in the transversal exponent which we focus our attention on, six-dimensional phase space. They are invariant since both methods being in good accord in the region of in- trajectories starting in each subspace will remain there terest (as shown in Fig. 6). For the chaotic attractorsin forever. Becausethedynamicsineachsubspacecoincides bothsynchronizationandantisynchronizationmanifolds, with that of the uncoupled map, then, it will evolve towards the respective well known Lorenz attractor A s,a lying in the corresponding manifold. Moreover,foreachinvariantsubspace,thereareoutof 1.0 three transversaldirections. Condition 3 means that the 0.8 transverse Lyapunov exponents of typical orbits lying in 1 the invariant subspaces (Ma and Ms) are all negative. 0.6 The point in parameter space where they become posi- 0 0.4 tive defines the blowout bifurcation [11]. To investigate condition 3, it suffices to consider the largesttransversal l 0.2 0 20 40 60 exponent, denoted as λ =lim λ˜ (x ,t)<0, where ⊥ t→∞ ⊥ 0 x is an initial condition on the basin of attraction of 0.0 0 either A or A . a s -0.2 We computed Lyapunov exponents using two different methods. The Lyapunov spectrum was obtained fol- -0.4 0.0 1.0 2.0 3.0 4.0 5.0 lowing the algorithm of Wolf et al. [27], with a Gram- e Schmidt normalization step of 0.1. We integrated Eqs. (2) using initial conditions given by x = y = x = FIG.6: ThethreelargestLyapunovexponentsofthecoupled 1 1 2 y =1.0 and z ,z were randomly chosen in the interval Lorenz system as a function of the coupling strength. Gray 2 1 2 [20,24)fromauniformprobabilityfunction,asinFig. 1. symbols correspond to the algorithm by Wolf et al., whereas These initial conditions lead to trajectories that asymp- thethickblackcurveistheresultofthevariationalequations (6). Inset: thelargest exponent for a wider range of ε. tote to either A or A . As a matter of fact, this is not s a 6 the largesttransversalexponent vanishes, changingsign, Figure 7 shows probability distribution functions at ε ≈ 0.714 ± 0.005 and ε ≈ 3.061± 0.005, defin- (when t = 30) for different values of the coupling 1 2 ing the criticalpoints of the blowoutbifurcation. (There strength. In all the considered cases the distribution is is also another critical value for large ε, as can be seen nearly Gaussian and presents positive tails. Then, we in the inset of Fig. 6, but we will restrict our analysis computed the positive fraction of finite-time exponents, to the lower range only). The largest transversal expo- ϕ(t) = R0∞P(λ˜⊥(t))dλ˜⊥(t) > 0. The positive fraction is nent is negative for ε < ε < ε , hence, in that interval, plotted in Fig. 8 as a function of the coupling strength 1 2 condition 3 for intermingled basins is fulfilled. However, for different values of the time-t interval used to sample while the invariant subspaces M are stable in aver- the finite-time exponents. s,a age,withnegativetransversalLyapunovexponents,there may be particular unstable periodic orbits embedded in the chaoticattractorsA that arealsotransverselyun- s,a stable, with positive largest transversal Lyapunov expo- 1.00 nent [28] (condition 4). When trajectories come close to theseunstableorbits,theywillberepelledfromthevicin- 0.75 ityoftheattractor. Thiswillbe reflectedinpositiveval- ues of the finite-time largesttransversalLyapunovexpo- nent [10–12]. Then we numerically computed the finite- j 0.50 time largest transversal Lyapunov exponents λ˜ (x ,t), t=5 ⊥ 0 t=10 by means of Eq. (5) but at finite t. For large enough t, 0.25 t=20 one recovers the infinite time exponent λ⊥, which does t=50 not depend on x0, for almost all initial conditions in the 0.00 t=100 attractors A , in contrast to the finite-time ones that s,a 0.0 1.0 2.0 3.0 4.0 5.0 may depend on the initial condition. e Wequantifythecontributionsofthefinite-timelargest transversal exponent by obtaining a numerical approx- FIG.8: Positivefractionofthelargesttime-ttransversalLya- imation to the corresponding probability distribution punov exponent as a function of the coupling strength, for function P(λ˜ (x ,t)). We considered a large number of different values of t. ⊥ 0 points in M (with x=y =z =0, X =Y =1.0, and Z s randomlychosenintheinterval[20,24)),anddiscardeda Forε.ε ,finite-timeexponentssoonbecomepositive. 1 transient. These were the initial conditions usedto com- This is in agreement with the positivity of the infinite- pute the time-t largest transversal Lyapunov exponents. time exponent shown in Fig. 6 for this region, and also Alternatively, we generated a single long chaotic trajec- with the fact that trajectories do not approach the in- tory (after the transient has elapsed) and divided it into variantsubspaces,butaresoonrepelled,asalreadynoted time-tsegments,usingthenthe ergodicityofthedynam- when we tried to plot the basins in Section III, which is ics to ensure that the conditions are randomly chosen notpossible forthatparameterrange. The positivefrac- according to the natural measure of the attractor. The tiondropsrapidlyto50%,whichcorrespondstothe case results were essentially the same. for which the infinite-time exponent vanishes (consistent with symmetric P(λ˜ (t))), and then drops below 50%, ⊥ whenthe infinite-time exponentis negative. For ε≃1.4, 101 the positive fractionis minimal. Forlargervalues ofε, it e=1.0 increases, crossing the 50% level again for ε≃ε2, where e=2.0 theinfinite-timeexponentisagainzeroatthatpoint(see 100 e=2.5 Fig. 6). Afterthat,thepositivefractiontendsto1gently with ε, yielding a positive infinite-time exponent. This F smooth behavior, different from the abrupt one in the D10-1 P lowerlimitoftheregionofnegativeλ ,isconsistentwith ⊥ the observation, for ε>ε , of a basin structure reminis- 2 10-2 cent of those in Fig. 1, although the filaments from the diagonalarenotneat. Eveniftrajectoriesareultimately repelled,theycanspendlongtimeintervalsclosetoeach 10-3 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 attractor. l~^ (t=30) The fraction of positive finite-time exponents is non- null. However, for the range ε < ε < ε , that fraction 1 2 FIG. 7: Probability density functions of the time-30 largest decreases with t, as expected because the distribution transversal Lyapunov exponent for different values of ε. Ini- of finite-time exponents collapses towards a Dirac delta tial conditions were taken as in Fig. 1. The full lines are centeredatλ inthelongtimelimit. Thedecayisexpo- ⊥ Gaussian fits. nential,the faster,the closerto the minimumat ε≃1.4. 7 Hence, the absence of an abrupt decay of the positive such that ~δ(τ) = Q~δ(0), with τ the time period of the fractionindicatesanonvanishingfractionforfinitetimes. orbit and ~δ(t) = (δx,δy,δz) the column vector of trans- Then,fromthisanalysis,condition4cannotdiscardedfor verse deviations. The eigenvalue of Q, µ, with maximal any range within the interval ε1 <ε<ε2. This suggests modulusfurnishesλ⊥ =ln|µ|/τ,foraparticularperiodic thatatleastoneofUPOsshouldbetransverselyunstable orbit. Fig. 9 shows the behavior of λ as a function of ε ⊥ in that interval. for particular UPOs, up to period 5. UPOs are labeled by means of the sequence of symbols A, B denoting the turnsaroundeachunstablefixedpointC+ andC− ofthe Lorenz system. Symmetric orbits obtained by exchang- 1.0 ing A↔B or with cyclic symmetry were omitted. One observes that the lowest period orbit AB (period 0.5 2) appears to be the first in destabilizing the vicinity of ^ 0.0 ε2, hence defining a riddling bifurcation. Then, between l this point and ε riddling can occur. This interval cov- 2 -0.5 ers most of the range ε < ε < ε , except for a very 1 2 small interval in the vicinity of ε . However, note that -1.0 (a) 1 orbitsofthetypeAnB,withn=1,2,...,haveamaximal transversalLyapunovexponents that increaseswith n in 0.0 0.5 1.0 1.5 2e.0 2.5 3.0 3.5 4.0 the vicinity of ε , hence shrinking the remaining small 1 regionofstabilityaroundε≃1. Toconfirmwhetherthis 1.0 region of strong stability (with no transversely unstable (b) AB orbits) actually disappears would require the analysis of AAB AAAB higherperiodorbits,ahardtaskforthissystem,sincethe AABB 0.5 AAAAB numberofUPOsincreasesexponentiallywiththeinteger AAABB period. ^ AABAB l chaotic Near the blowout bifurcation at ε , the low-period 1 UPOs (up to period 5) destabilize for coupling strength 0.0 either weaker or stronger than the critical value, but close to it. Let us remark that, differently to the coupled Ro¨ssler system studied by Heagy et al. [28], here 0.0 0.2 0.4 0.6 e 0.8 1.0 1.2 1.4 the ordering of the exponents of the lowest period orbits inthe neighborhoodofε isinvertedwithrespectto 1 AB (c) the uncoupled case as depicted in Fig. 9. This implies AABAB bit AAAABBB that paradoxically the most stable orbits in the attrac- or AAABB tor are those responsible for the transversal destabiliza- AAAB AAAAB e e tion in this parameter region. A similar inversionoccurs 1 2 on some domains of the parameter space of a system of 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 e symmetrically coupled Ro¨ssler oscillators[15, 17]. This characteristicturnsdifficultthedeterminationoftherid- dling bifurcation (first destabilized orbit) related to the FIG. 9: (Color online) (a) Largest transversal Lyapunov blowout at ε , apparently triggered by higher period or- exponent,as afunction of ε,for theparticularunstableperi- 1 bits. odicorbits(UPOs)embeddedintheLorenzchaoticattractor (up to period 5), indicated on the figure by means of the Furthermore,our outcomes point to a different nature sequenceofsymbolsA,Bdenotingtheturnsaroundeachun- of the blowout bifurcations at ε1 and ε2. Nearby ε1, stable fixed points C+ and C− of the Lorenz system. The UPOsdestabilizeinitsvicinity. Moreover,forallthean- curve for typical chaotic trajectories in the attractor is also alyzed orbits, the multiplier µ crosses the circle |µ| = 1 shown. (b) Magnification of panel (a). (c) Stability intervals along the real positive semi-axis (associated to a pitch- for each UPO, in order of increasing stability at ε = 0 from forkbifurcation). Thisisincontrasttothescenarioatε , 2 top to bottom: stable (full segment), unstable (dotted). The where there are orbits destabilizing far from ε and with 2 vertical lines indicate ε1 and ε2. The symbols delimiting the multiplierµeither+1or−1. Inparticular,thefirstorbit stabilityintervalscorrespondtoµ=1(full)and−1(hollow). ABlosesstabilitywithµ=1. Thedifferencesareconsis- tent with the picture given by finite-time exponents, for Then we inspected the transversal stability of those instanceinconnectionwithFig. 8,whereamuchabrupt orbits along the lines of periodic orbits threshold theory behavior of the positive fraction was encountered near [28]. Once localizationinphase space andperiods oflow ε1. periodUPOsareavailableintheliteraturefortheLorenz Theintervalswhereriddlingcanoccuraredelimitedon system [29], we computed Floquet multipliers [28, 30]. one side by the blowout bifurcation and on the other by Namely, we integrated Eqs. (6), to obtain the matrix Q theriddlingbifurcation. Inourcasetherearetwoofsuch 8 z intervals and they apparently overlap,such that at least (a) one UPO has lost transversal stability in the full range x = y = 0 ε < ε < ε , although this would have to be confirmed 1 2 by the analysis of high period orbits, it is supported by the analysis of finite-time Lyapunov exponents. l/2 V. SCALING LAWS FOR RIDDLED BASINS In this Section, we will focus on the determination of thescalingpropertiesofthebasins,whichprovideamea- 0 X, Y, Z sure of their structure[10]. Let us focus on the black filaments in Fig. 1(b), which belong to the basin of the antisynchronizationattractor. They are anchored at the 100 (b) diagonalline, whichis a cut ofthe synchronizationman- ifold M , given by x = y = z = 0 and containing a s chaotic attractor As, while the antisynchronization at- 10-1 tractor A lies elsewhere. In Fig. 10(a) we portrait a a schematic picture ofthat structure. The filaments ofthe P★ basin of A are tongues anchored at points of A , and a s 10-2 the complement of the filament set belongs to the basin e=1.0 of As. If an initial condition starts within any of these ee==12..50 narrowtongues,evenifitisveryclosetoA ,theresulting e=2.4 s e=2.8 trajectory will asymptote to the other attractor. 10-310-3 10-2 10-1 100 101 The set of basin filaments for A is expected to be l a self-similar by quite general grounds. Once the riddling 100 (c) bifurcation occurs for a given periodic orbit, it also oc- cursfor everypreimageofthis orbit,yieldinga denseset of tongue-like sets anchoredat the corresponding preim- 10-1 agesonA . Thetongue-likeshapeisaconsequenceofthe s nonlineartermsinthe equationsdescribingthe transver- P★ sal dynamics. The characteristic feature of riddling is 10-2 that those tongues have widths that tend to zero as we |x|=l/2 approach A . Hence the basin of A always contains s s |z|=l/2 pieces of the basin of the other attractor, regardless the d=l/2 transversaldistance to As, so forming a fine structure of 10-310-2 10-1 100 101 basin filaments (the same applying to A ). l a Thisfine structurecanbe quantitativelycharacterized FIG.10: (a)Schematicfigureshowingthestructureofriddled by the following numerical experiment [11, 23]: let us basinsneartheinvariantsubspacethatcontainsachaoticat- consider the invariant manifold at x = y = z = 0 and, tractor. (b)Fractionoftrajectories P⋆ thatasymptotetothe depart from that manifold, for instance, by increasing z, antisynchronized state as a function of the distance l = 2|z| up to a distance l = 2|z|≡ |z −z | [as depicted by the 1 2 to thesynchronized state, for different values of the coupling red line segment in Fig. 10(a)]. Then we evaluate the strength ε. The full lines are least squares fits. (c) P⋆ for fraction Vl of points in that segment that belongs to the ε=1.5 and different orientations of the deviation l. basin of A . We obtained a numerical approximation of s thisfractionbyconsideringanumberofinitialconditions x=y =0, |z|=l/2,X =Y =1.0,andZ randomlycho- P =1−V,andisexpectedtoscalewithlasapowerlaw ⋆ l sen in the interval [20,24). If the trajectories did not P (l) ∼ lη, where η > 0 is a scaling exponent. We inte- ⋆ synchronize (within a small tolerance)up to a time such gratedseveralinitialconditionsatthe samedistancel to that transients have elapsed and stationarityholds (typ- the synchronizationsubspace and computed the fraction ically, t = 103), we consider that they asymptote to the P ofinitialconditionsthatdonotsynchronize,repeating ⋆ antisynchronization attractor A and, accordingly, they this procedure varying the distance l. The results shown a do not belong to the basin of A . If the latter is riddled inFig.10(b)confirmtheexistenceofapowerlawforthis s with tongues belonging to the basin of A , for any dis- fraction, for many values of the coupling strength. a tance l (no matter how small) there is always a nonzero The results do not vary appreciably when one departs value of V . This fraction tends to zero as l → 0. The from the synchronization manifold in other directions l fraction of length belonging to the basin of A (fraction other than z. In Fig. 10(c) we plotted, for the same a oftrajectoriesthatdonotsynchronize)canbewrittenas coupling strength (ε = 1.5), the fraction P for initial ⋆ 9 conditions with |z| = l/2, x = y = 0 (open squares), estimates η = |λ |/D are plotted in Fig. 11 together ⊥ and also for |x| =l/2, y = z =0 (filled circles) and ran- withthenumericalvalues. Numericalandestimatedval- dom values of x,y,z such that d≡px2+y2+z2 =l/2 ues are in very good agreement in the proximity of the (filledsquares),obtainingessentiallythesamescalingex- critical values, as expected [34]. For intermediate values ponent. This scaling behavior is observed within the in- (0.75<ε<2) there is a discrepancy, and the numerical terval (ε ,ε ), below ε , no trajectories synchronize as exponent remains close to one (linear behavior), as also 1 2 1 seen in the previous sections, above ε , one observes a observedforothersystemswithintermingledbasins[20]. 2 synchronizedfractionbut itdoes notchangewith l. The dependence of the numerically determined scaling expo- e nents η on the coupling strength is depicted in Fig. 11. 10-1 e==01..80 e =1.2 e 101 e=1.5 =2.0 Numerical e Theoretical )10-2 e==22..58 t 100 2l( s 10-3 10-1 h 10-4 10-2 101 102 103 t 10-3 FIG.12: Timedecayofthevarianceofthefinite-timelargest 0.5 1.0 1.5 2.0 2.5 3.0 3.5 e transversalLyapunovexponentfordifferentvaluesofthecou- plingstrengthε. Thelinesareleastsquaresfitsofthefunction f(t) = 2D/t to the numerical points, for large t, allowing to FIG. 11: Scaling exponent η for the fraction of trajectories estimate the diffusion coefficient D. that asymptote to the antisynchronized state obtained by a numericalexperiment(filledcircles). Forcomparison,thethe- Another scalinglawtypicalofriddledbasins is related oretical values, given by η = |λ⊥|/D are also plotted (open to the fraction of uncertain initial conditions, with re- squares). spect to their final-state [31]. We may regard riddled basins as an extreme case of fractal basins, for which An analytical expression for the exponent η was de- thereisfinal-statesensitivityandtheuncertaintyfraction rivedbyOttandcoworkersforasimplemodel(piecewise scales as a power-law with the uncertainty level, whose linear non-invertible map) [23]. Their theoreticalpredic- exponent gives a measure of the extreme final-state sen- tion arises from a diffusion approximation for a biased sitivity due to riddling. Consider again the points at random walk that mimics the fluctuations of finite-time x = y = 0 and |z| = l/2 drawn in the phase space por- largesttransversalLyapunovexponents λ⊥(t). They ob- trait in Fig. 13(a), as described earlier, and choose ran- tain the law P⋆ ∼ lη, with η = |λ⊥|/D where D is the domlyaninitialconditionx0 onthatregion. Nowchoose diffusion coefficient. This diffusion approximation is ex- randomlyanotherinitialconditionx′ withuniformprob- 0 pectedto be validnearthe blowoutbifurcation(λ⊥ ≃0) abilitywithin anintervaloflength2ξ andcenteredatx0 ofanattractorwithariddledbasin. Theauthorsconjec- [Fig. 13(a)]. If both points belong to different basins, turethatasimilardiffusionapproximation,henceasimi- they can be referred to as ξ-uncertain [32, 33]. larrelationinvolvingparametersλ⊥andD,mustrulethe Thefractionofξ-uncertainpoints,oruncertaintyfrac- scaling relationin a largeclass of systems. The distribution, denoted by hpi, is the probability of making a mis- tions showninFig. 7 alreadydisplaya Gaussiancharac- take when attempting to predict which basin the initial ter,whichimproveswithlargertime-tinterval,consistent condition belongs to, given a measurement uncertainty withtheprobabilitydistributionfunctionofindependent ξ. This probability scales with the uncertainty level as a random innovations, that, by the central limit theorem, power law of the form hpi∼ξφ, where φ≥0 depends on isGaussian. Additionally,weplotinFig.12thevariance bothx andl. NumericalresultsareshowninFig. 13(b). σλ2⊥(t) of the probability distribution functions for λ˜⊥(t) The sto0chastic model of Ott et al. predicts a power-law, as a function of time, for different values of the coupling with exponent given by φ = λ2/(4Dλ ), a prediction ⊥ k strength. As a matter of fact, the variance decays with that agrees with our numerical results close to the crit- time towards zero following asymptotically a power-law ical points. However, for intermediate values, while the withexponent−1,asrequiredforanormaldiffusionpro- stochastic model predicts small (though nonzero) values cesses,sovalidatingthestochasticapproachofOttetal.. (φ < 0.28), our numerical results for φ, as illustrated in Accordingly,the diffusion coefficientD canbe estimated Fig. 13(b),yieldmuchsmallervalues(byafactorgreater fromthenumericalcurves,followingσ2 ∼2D/t. The thanten). Asamatteroffact,forriddledbasins,theex- λ⊥(t) 10 z (a) bits. This is important as furnishes the sources of local transversal instability of the attractor even if stable in average. In a second place, we verified the existence of x = y = 0 two scaling laws characterizingquantitatively the degree of uncertainty related to the riddled basins. These nu- mericalresultswerecomparedtoananalyticalprediction (thestochasticmodel[23]),yieldingagoodaccordwhere l/2 expected. Beyond the characterization of the structure ξ of a riddled basin, these scaling laws allow to quantify the limitations to improve the ability in determining the 0 X, Y, Z finalstateofthesystembyincreasingtheaccuracylevel. Letusremarkthatintermingling,inparticularofsym- metric basins, has also been observed in other systems 100 (b) with either continuous (e.g., mechanical system [12] and coupled Ro¨ssler oscillators[17]) or discrete time dynamics (coupled logistic maps [18]). In the latter case, the analysisof the lowestperiod orbitwas enoughto furnish pæ 10-1 the conditions for the occurrence of riddling in certain Æ parameter region. In fact, as anticipated by the results presentedinFig. 9, deepening in thatpoint of view may e=1.0 furnish precise information on the nature of the bifurca- e=1.5 e=2.0 tions triggeringriddling,althoughthis may be adifficult 10-210-8 10-7 10-6 10-5 10-4 10-3 10-2 taskforthepresentsystem. Asotherperspectivesforfu- x tureworkonthissystem,letusalsomentiontheplausible occurrence,beyond the blowoutbifurcation, of two-state on-offintermittency[12]forwhichthereissomeevidence FIG.13: (a)Schematicfigureshowingthenumericaldetermi- [21]. Finally, it can still be worthy to explore other re- nation of the uncertainty fraction. (b) Fraction of uncertain gionsofphase space,as wellasother ranges(negativeor initial conditions as a function of the uncertainty level, for large values) of the coupling parameter. different values of the coupling strength. The solid lines are In any case, for applications, multistability is already least squares fits. a source of troubles. Still worst, the existence of inter- mingledbasinsofattractionforthesynchronizedandan- ponent φ should be rigorously zero (i.e. there would be tisynchronized chaotic states of this system jeopardizes thesolutionoftheproblemofensuringagivenfinalstate, no way to decrease the uncertainty fraction by decreas- ingtheuncertaintylevel). Indeed,ourresults[Fig. 13(b)] since the initial condition determination is always done support this scaling law, with numerically obtained ex- within a certain uncertainty level. With riddled basins, ponents close to zero. any uncertainty level, however small, lead to complete indeterminacy of the future state of the system. Hence in this case we cannot use synchronization of chaos for VI. CONCLUSIONS AND FINAL REMARKS any practical purpose, since we will always be haunted by the existence of the another, antisynchronized state, with a basin intermingled with the basin of the synchro- Riddledbasinsforthesynchronizationattractorofcou- nized state. Of course, the same difficulties concern the pled Lorenz systems have been previously suggested in predictabilityofnaturalphenomenamodeled bycoupled theliteraturebutwithoutadetailedcharacterization. In Lorenz systems. Therefore, the importance of detecting this work we offer numerical evidence that, for a spec- the regimeswhereriddlingcanoccurinadynamicalsys- ified range of the coupling parameter (ε < ε < ε ), 1 2 tem. coupled Lorenz systems exhibit symmetrically riddled basins of attraction for synchronized and antisynchronized states. Since there are only two symmetric attractors, their basins are intermingled. We firstly showed Acknowledgements: thatthemathematicalconditionsfortheexistenceofrid- dled basins are fulfilled, with the help of properties of This work was made possible with help of CNPq, finite-time largest transversal Lyapunov exponents and CAPES, FAPERJ, and Funda¸caõ Araucária (Brazilian of the largest transversal exponent for particular or- Government Agencies). [1] E. N.Lorenz, J. Atmosph.Sci. 20, 130 (1963). [2] H. Haken,Phys. Lett.A 53, 77 (1975).

Intermingled basins in coupled Lorenz systems

Similar Intermingled basins in coupled Lorenz systems