Intermingled basins in coupled Lorenz systems

4.3 MB·English

by  Sabrina Camargo

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 synchro-
nization 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 sta-
bility 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 synchro-
nized 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 syn-
chronization) [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 to-
wards 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 differ-
ent 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 cou-
pled 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 or-
bits 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 fil-
aments 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 distribu- tion, 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 dynam-
ics (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 antisynchro-
nized states. Since there are only two symmetric attrac-
tors, 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˜o Arauc´aria (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).