SYNCHRONISATION AND EXTREME VALUE THEORY FOR COUPLED MAP LATTICES Abstract. We show that the probability of appearance of synchronisation in chaotic coupled map lattices is related to the distribution of the maximum of a certain observable evaluated along almost all orbit. We show that such distribution belongs to the family of extreme value laws, whose parameters, namely the extremal index, allow us to get a detailed description of the probability of synchronisation. Theoretical results are
arXiv:1708.00191v1 [math.DS] 1 Aug 2017
supported by robust numerical computations that allow to go beyond the theoretical framework provided and are potentially applicable to physically relevant systems.
D. Faranda 1, H. Ghoudi2, P. Guiraud3, S. Vaienti4 1. Introduction We present in this paper a new application of Extreme Value Theory (EVT) to Coupled Map Lattices (CML) on a finite lattice. Actually a first result in this direction was given in the paper [10], although not explicitly related to EVT, where the authors considered two coupled interval maps and applied their theory of open systems with holes to investigate the first entrance of the two components into a small strip along the diagonal, which is equivalent to the synchronisation of the two-components lattice up to a certain accuracy. We will show that synchronisation processes could be interpreted and quantified by computing the asymptotic distribution of the maximum of a suitable random process. The rigorous generalisation to lattices with more than two components requires strong assumptions on the system; in more general situations we will present arguments to sustain the existence of a limit distribution for the maxima related to synchronisation, and we will discuss a formula approximating the extremal index for lattices with arbitrary number of components. We therefore estimate the behavior of such an index when the number of components is large. We will then generalize the theory to CML which are randomly perturbed with additive noise and show, with numerical evidence, that the Date: August 2, 2017. 1 LSCE-IPSL, CEA Saclay l’Orme des Merisiers, CNRS UMR 8212 CEA-CNRS-UVSQ, Universit´e Paris-Saclay, 91191 Gif-sur-Yvette, France . E-mail:
[email protected]. 2 Laboratory of Dynamical Systems and Combinatorics, Sfax University Tunisia and Aix Marseille Universit´e, CNRS, CPT, UMR 7332, 13288 Marseille, France and Universit´e de Toulon, CNRS, CPT, UMR 7332, 83957 La Garde. Email:
[email protected] . 3 CIMFAV, Facultad de Ingeniera, Universidad de Valparaso, Valpara`ıso, Chile.
E-mail:
[email protected]. 4 Aix Marseille Universit´e, CNRS, CPT, UMR 7332, Marseille, France and and Universit´e de Toulon, CNRS, CPT, UMR 7332, 83957 La Garde. E-mail:
[email protected]. 1
extremal index will be 1 for any dimension of the lattice. Synchronisation is usually intended to last for a while once it started and this is what usually happens for some kinds of chains of synchronised oscillators. This is not the case of course for our chaotic CML, since almost every orbit is recurrent by Poincar´e theorem. What we actually investigate is therefore the probability of a first synchronisation and how long we should wait to get it with a prescribed accuracy. We will see that EVT provides all these kind of quantitative information. Although we could not get a global synchronisation persisting in time, we could ask about the distribution of the number of successive synchronisation events when the systems evolves up to a certain time. We will see that after a suitable rescaling, the distribution of that number will follow a compound Poisson statistics: it is worth mentioning that for two uncoupled expanding maps of the circle, this result dates back to a paper by Coelho and Collet, [5]. A few words on the structure of the paper. In Section 2 we will present a powerful and general approach based on perturbation of the transfer operator, and which has the advantage of being applicable to a large class of observables arising in the study of EVT. In Section 3, after Proposition 3.1, we will give a very short insight into EVT, especially when it is applied to recurrence in dynamical systems and we will define the extremal index. Our analytic results and estimates are supported by numerical computations in the last part of the paper; they confirm the existence of an extreme value distribution for a different kind of synchronisation, which we called local in section 7, and they validate the expected compound Poisson statistics for the distribution of the number of successive visits. It is an interesting numerical discovery the fact that the extremal index for local synchronisation seems not depend on the size of the lattice. The kind of CML which we studied is particularly simple; it represents a global coupling where each map interacts with the mean value of all lattice components, see [15] where this mean-field model was introduced and investigated 5. Our aim in this work is to provide a first approach to coupled map lattices by using extreme value theory and to show how to get a certain number of rigorous results. In forthcoming papers we will study more general CML with non-local form of coupling including the important case 5Kaneko
was one of the first, together with Kapral and Kuznetsov, to introduce and investigate CML.
The first contribution which looked at these maps in the framework and with the tools of ergodic theory and statistical properties of dynamical systems, was the work by Bunimovich and Sinai [3]. Since then the progress to understand CML have been enormous with contributions of several people. We defer to the book [4] for a wide panorama on the different approaches to CML and for exhaustive references. 2
of diffusive or Laplacian interaction; a few other possible developments will be presented at the end of the paper (see section 7.2). We hope that our approach could be helpful to understand and quantify those phenomena where bursts of synchronisation, like in neuronal spikes or in business cycles of financial markets, happen, disappear, happen again, apparently in a disordered manner, but very often following the extreme distributions arising in chaotic systems. 2. The map and the operators As we said in the Introduction we will consider a very simple example of CML: Tˆ : I n → I n , where I n := [0, 1]n and6 n
γX xi = Tˆ(x)i = (1 − γ)T (xi ) + T (xj ) n j=1 0
(2.1)
where x = (x1 , . . . , xn ) ∈ I n , γ ∈ (0, 1) and T is a piecewise expanding map of the unit interval onto itself with a finite number of branches, say q, and which we take of class C 2 on the interiors of the domains of injectivity A1 , . . . , Aq , and extended by continuity to the boundaries. This system represents a global coupling where each component is coupled to all the others; it was first introduced by Kaneko in [15] and successively studied, among others, by P. Ashwin [2] (and references therein). Let us denote with Uk , k = 1, . . . , q n , the domains of local injectivity of Tˆ; by the previous assumption on T there will be open sets Vk ⊃ Uk such that Tˆk := Tˆ|Vk is a C 2 diffeomorphism (on the image). We will require that sup sup ||DTˆi−1 (x)|| < sn < λ < 1 i
x∈Tˆi Vi
where λ = supi supx∈Ti Ai |DTi−1 | and || · || stands for the euclidean norm. We will write dist for the distance with respect to this norm. An important tool for our further considerations is the transfer, or Perron-Fr¨obenius (PF) operator; we design with Pˆ the Perron-Fr¨obenius operator of the map Tˆ: it is simply defined by the duality integral relation Z Z ˆ P (f )gdLeb = f g ◦ TˆdLeb, where Leb denotes the Lebesgue measure on I n , f ∈ L1 , g ∈ L∞ .7 The spectral properties of the PF operator become interesting when it acts on suitable Banach spaces. Let us will not index the map Tˆ with n, hoping it will be clear from the context. the following we will use the same symbol Leb for any n; moreover L1 , Lp and L∞ will be taken R with respect to Leb. Finally integral with respect to Lebesgue measure will be denoted with dLeb(x) R or dx. 6We
7In
3
therefore suppose that there exists a Banach space B with norm ||·||B , which is compactly injected in L1 and the following properties hold8: • P1 For any f ∈ B there exists η < 1 and C > 0 such that ||Pˆ f ||B ≤ η||f ||B + C||f ||1 . This is called the Lasota-Yorke inequality; it implies that Pˆ has an isolated eigenvalue equal to 1 and which is also the spectral radius of Pˆ (spectral gap property). We will often call η the contraction factor in the Lasota-Yorke inequality. Sometimes one fails to get such an inequality for the map Tˆ; in this case it will be enough to recover it for a given iterate of Pˆ . • P2 We suppose that 1 is a simple eigenvalue with no other eigenvalues on the unit circle; this implies that Pˆ preserves a mixing measure µ ˆ which is the unique absolutely continuous invariant measure with respect to Lebesgue. We moreover ˆ ∈ L∞ . assume that the density h Properties P1, P2 could be achieved by taking the map T with its transfer operator having a spectral gap and preserving a mixing absolutely continuous invariant measure. Then by perturbation theory [11], it will be enough to check that, for any f ∈ B: ||(Pˆ − P0 )f ||1 ≤ pγ ||f ||B , where P0 is the PF operator of the uncoupled system (γ = 0), and pγ is a monotone upper semi-continuous function going to zero when γ converges to 0. If this happen P2 will follow. The aforementioned perturbation theory will be used in a more stringent way and in the version presented in [10], to control a different perturbation of the operator Pˆ , which arises naturally in the context of the EVT as we will see in the next section. In order to define this perturbation, we consider two connected and disjoints domains D1 , D2 ∈ I n , which are the closures of their interiors and with piecewise C ∞ and codimension 1 boundaries; set D = D1 ∪ D2 . Moreover D will depend upon the index l, and for that we rename it as Dl , in such a way that Leb(Dl ) → 1, l → ∞. The perturbed operator is defined as, for any h ∈ B: Pel (h) := Pˆ (h1Dl ) • P3 The operator Pel verifies a Lasota-Yorke inequality, uniform in l, on the space B, namely the factors η and C (see above) will be the same for all sufficiently large l. This is the first step of the comparison between the operators Pˆ and Pel with the objective to get a perturbative formula for the spectral radius of Pel close to 1 for large values of l; we will see that it will give us the extremal index in the limiting distribution of the Gumbel’s law. Two more properties are required: 8
We will call a Banach space with this property adapted (to L1 ). 4
• P4 For any h ∈ B, the quantity Z rl :=
|
sup
(Pˆ h − Pel h)dLeb|
h,||h||B ≤1
must go to zero when l → ∞. ˆ of the unperturbed invariant measure µ • P5 For the density h ˆ, we need: ˆ B ≤ C 0µ rl ||(Pˆ − Pel )h|| ˆ(Dlc ),
(2.2)
where C 0 is a constant independent of l. We will moreover assume that the density ˆ is strictly positive, namely its infimum is larger that h ˆ (inf ) > 0 on a set of full h measure. Under the assumptions (P1)-(P5), we have that the spectral radius ρl of Pel verifies 1 − ρl = µ ˆ(Dlc )θ(1 + o(1)), in the limit l → ∞
(2.3)
∞ X 1 − ρl lim qk = θ := 1 − l→∞ µ ˆ(Dlc ) k=0
(2.4)
where
being R qk := lim qk,l := lim l→∞
l→∞
ˆ (Pˆ − Pel )Pelk (Pˆ − Pel )(h)dLeb µ ˆ(Dlc )
(2.5)
and provided the limits exist. We stress that ρl is the largest eigenvalue of Pel , that there are no other eigenvalues on the circle of radius ρl , and that there are functions gˆl ∈ B and measures µ ˆl for which the operators Pel satisfy, for h ∈ B : Z e Pl h = ρl gˆl hdˆ µl + Ql h. (2.6) Moreover
R
gˆl dˆ µ = 1,
R
hdˆ µl →
R
hdˆ µ, for l → ∞ and finally Ql is a linear operator with
spectral radius strictly less than ρl and verifying: ||Qnl ||B ≤ ςln , for a suitable 0 < ςl < 1, see again [10] for the derivation of these formulae. It is a remarkable fact that this approach automatically provides the good scaling exponents for the asymptotic distribution of the maxima, see eq. (4.21) below and therefore it gives a new proof of the existence of that distribution. The quantity θ is called the extremal index (EI) and it will play an important role in the following; we will see in particular that it gives a correction to the pure exponential law for the distribution of the maxima. In that respect it coincides with the extremal index as it is defined in EVT, see [7], [17]. Our next task will therefore be to look for a Banach space which verifies the preceding five properties and successively prove the existence of the limits (2.4) and (2.5) for our particular systems. One natural candidate would be the space of functions of bounded variation on Rn restricted to the L1 functions supported on Ien := interior(I n ), but in this case it seems difficult to obtain (P5). This space was used in [10] in dimension 5
2, but it seems not generalizable in higher dimensions 9. We therefore turn our attention to another functional space, the quasi-H¨older space, and we start by defining, for all functions h ∈ L1 (I n ), and 0 < α ≤ 1 the seminorm Z 1 |h|α := sup a osc(h, Bε (x))dLeb 0 um ) → τ, m → ∞. We say again that the sequence Mn has an Extreme Value Law, if there exists a nondegenerate distribution function H : R → [0, 1], with H(0) = 0 such that µ ˆ(Mm ≤ um ) → 1 − H(τ ), m → ∞. By using the expression of ψ we can rewrite (3.11) as (n) mµ ˆ(Sm )→τ
(4.17)
(n) Sm := {x ∈ I n ; max |xi − xj | ≤ νm }; where νm := e−um
(4.18)
i6=j
and consequently (3.12) can be restated as (n) µ ˆ(x ∈ I n ; Tˆk (x) ∈ / Sm , k = 0, . . . , m − 1) → 1 − H(τ )
12
(4.19)
The limit (4.19) could also be interpreted as the probability that the n components have synchronized for the first time after m iterations with accuracy ac of order e−um . We cannot use the Pure Probabilistic Approach to prove the existence of the limit (4.19). The reason is that our new observable becomes infinite on a line (the diagonal), and for the moment rigorous results are avalaible when the set of points where the observable is maximised is at most countable, see [6] for a discussion of these problems. The Spectral Approach will bypass that issue by using the Banach space B given by quasi-H¨oder functions, since for such a space we can check properties (P1)-(P5). Nevertheless there will be a problem remaining, namely prove the existence of the limits (2.4,), (2.5); we will return on that in the next section. We now show, as promised, how the spectral approach allows us to get the asymptotic distribution functions of the extreme value theory. Let us begin by rewriting the maximum given in (4.19) as (we put ˆ the density of the absolutely continuous invariant measure µ again h ˆ.) Z µ ˆ(Mn ≤ um ) =
ˆ ˆ h(x)1 (n) (x)1 (n) c (T x) . . . 1 (n) c (T (Sm )c (Sm ) (Sm )
ˆm−1
Z x)dLeb =
ˆ Pemm (h)dLeb (4.20)
where, from now on, Pem (·) := Pˆ (1(Sm (n) c ·) ) (n)
Notice that (Sm )c plays the role of the set Dl in section 2. By invoking the spectral representation (2.6) we have with obvious interpretation of the symbols Z Z Z m ˆ m ˆ ˆ e Pm (h)dLeb = ρm hdˆ µm + Qm m hdLeb, ˆ hdLeb = 1, as m → ∞, and the spectral radius of Qm is strictly less than ρm . We now need a good frame to bound ρm , the largest eigenvalue of Pem , for
where
R
ˆ µm → hdˆ
R
increasing m and it is given by (2.3), namely, with obvious adaptation of notations (n) 1 − ρm = µ ˆ(Sm )θ∆ (1 + o(1)), in the limit m → ∞;
then Z
(n)
ˆ Pemm (h)dLeb = e−(θ∆ mˆµ(Sm
(n)
)+mo(ˆ µ(Sm ))
Z
ˆ µm + O(ρ−m ||Qm ||B ) hdˆ m m
(4.21)
which converges to e−τ θ∆ under the assumptions on µ ˆm , the spectral radius of Qm and the condition (4.17).
The exponent θ∆ denotes now the EI along the diagonal set
n
∆ := {x ∈ R ; x1 = x2 , · · · = xn } and its existence will follow if we prove the limits (2.4) and (2.5). From now on we will simply write θn for the EI along the diagonal set for lattices with n components. 13
We now return to (4.19) since we now know that 1 − H(τ ) = e−θn τ . If we suppose 1 (n) τ n−1 n−1 13 and therefore the probability of the first that µ ˆ(Sm ) = O(νm ) , then e−um ∼ m 1 τ n−1 synchronisation after m iterations and with accuracy ac ∼ m , is e−θn τ .14 If the components of the vector Tˆk (x) are seen as the positions of different particles on a lattice or on a network, we have a quantitative estimate of the probability of synchronisation of the lattice after a prescribed time and with a given accuracy. Example 4.1.
• Ex 1. Suppose we use the data in Section 7, with an EI θ3 ∼
10 2 1 − ( 27 ) ∼ 0.86, having chosen λ = 1/3 and γ = 0.1, and take 3 particles each
living on the unit interval. If we want to synchronize them with a probability larger than 1/2 and an accuracy ac = 0, 01 before m iterations, then we have to iterate the lattice around m = 8 100 times. • Ex. 2 Analogously, if we want to observe with a probability larger than 1/2 the synchronization of 100 particles each living on the unit interval with an accuracy ac = 0, 01 and before m iterations of the CML, then m has to be larger than 100100 . 5. Computation of the extremal index The extremal index is given by formula (2.4); Keller showed in [9] that it coincides with that given in Proposition 3.1 for the process (3.9) and the proof is exactly the computation we performed at the end of the previous section. As we said in the Introduction, the rigorous computation of the EI for two coupled maps was given in [10]. Their map was slightly different from ours in the sense that for the i-th component the averaged term γ Pn j=1 T (xj ) does not contain the contribution of T (xi ). They first observed that in (2.5), n all the qk but q0 , are zero due to the fact that the diagonal is invariant and q0 reads: q0 = lim
(2) (2) µ ˆ(Sm ∩ Tˆ−1 Sm ) (2)
m→∞
(5.22)
µ ˆ(Sm )
This quantity was explicitly computed giving the formula [10]: Z ˆ 1 h(x, x) 1 θ2 = 1 − dx, R ˆ x)dx 1 − 2δ h(x, |DT (x)| ˆ has bounded variation and for almost every x ∈ I the value h(x, ˆ x) where the density h ˆ − u, x + u) and h(x ˆ + u; x − u) as u → 0. is the average of the limits of h(x 13
Actually this is a very crude approximation.
In fact what is possible to prove easily is an
upper bound on the Lebesgue measure of the domain {x ∈ I n , |xi − xj | < νm ; i 6= j} which is simply (2νm )n−1 . We sketch the argument for n = 3; in this case the measure we are looking for R R R is dx1 dx2 1{|x1 −x2 |≤νm } (x) dx3 1{|x1 −x3 |≤νm } (x)1{|x2 −x3 |≤νm } (x). The last integral will contribute with 2νm and so the second one. 14We defer to the discussion after Proposition 3.1 for the validity of this argument and its approximations. 14
We get a similar result and still for n = 2, with a modification due to the fact that our map is different, see formula (5.30) in the remark below; instead the density along the diagonal is defined again as a bounded variation function with the same averaging procedure. It seems difficult to extend such a result in higher dimensions without much stronger assumptions; before doing that we will explore how the EI θn behaves for large n in a quite general setting with the objective to show that for large n such an index approaches 1 and therefore the Gumbel’s law will emerge as the extreme value distribution. ˆ n , while we continue to use the symbol µ We will index with n the invariant densities h ˆ ˆ n . We for the invariant measure, despite the fact µ ˆ depend on n too, via the density h need first an additional property, namely: • P0: The one-dimensional map T is continuous. Notice that with this assumption, the map is not necessarily onto on its domains of injectivity and it remains not differentiable on the boundaries of such domains. We need continuity to show that all the qk but q0 are zero. We first notice that the quantities qk,l introduced in (2.5), read: R ˆ (Pˆ − Pel )Pelk (Pˆ − Pel )(h)dLeb qk,l := =µ ˆDlc {x ∈ Dlc ; tDic (x) = k + 1} µ ˆ(Dlc )
(5.23)
where µ ˆDlc is the conditional measure to Dlc , and tDic (x) denotes the first return time of the point x ∈ Dlc to Dlc (we will come back on this equality in the next section). We sketch the argument for k = 1, the others being similar. By replacing l with m in (5.23) we show that for m large enough the quantity Z (n) (n) c (n) ˆ dLeb = µ (Pˆ − Pem )Pem (Pˆ − Pem )(h)f ˆ(Sm ∩ Tˆ−1 (Sm ) ∩ Tˆ−2 Sm )
(5.24)
is zero. By continuity of the maps Tˆ and Tˆ2 , there will be an open set O ⊃ ∆ such that whenever x ∈ O and Tˆ−1 is the inverse branch of Tˆ such that Tˆ−1 Tˆ2 (x) = Tˆ(x) ∈ O, then β
β
Tˆβ−1 ∆ ⊂ ∆ (remember the diagonal is Tˆ-invariant). We now take m large enough in such a (n) (n) (n) way that Sm ⊂ O and choose a point x ∈ Sm and such that Tˆ2 (x) ∈ Sm . We now show (n) that Tˆ(x) must be in Sm too. Take a point y ∈ ∆ and such that dist(y, Tˆ2 (x)) < νm . The preimages of those two points under Tˆβ−1 will be at a distance less that νm since the inverse branches strictly contract, and the point Tˆ−1 (y) will be on the diagonal: so β
(n) Tˆ(x) ∈ Sm .
Proposition 5.1. Let us suppose our CML verifies properties (P1)-(P5) on a Banach ˆ n ∈ L∞ and h ˆ (inf) ˆ n > 0. space B with λ = inf |DT |−1 < 1 − γ, the density h := inf I n h n Then lim sup m→∞
(n) (n) µ ˆ(Sm ∩ Tˆ−1 Sm ) (n)
µ ˆ(Sm ) 15
ˆ n ||∞ λn−1 ||h ≤ ˆ (inf) (1 − γ)n−1 h n
Remark 5.2. The upper bound makes sense of course when the right hand side of the above inequality is less or equal to 1. Moreover the EI θn will converge to 1, under the additional P0 assumption, when n → ∞ if the ratio with υ > (λ/(1 − γ))−1 .
ˆ n ||∞ ||h (inf) ˆ hn
will not grow faster than υ n−1
Proof. We start by writing µ ˆ
(n) Sm
(n) ∩ Tˆ−1 Sm = Z
Z dx1
I n−1
I
1Sm(n)
Z In
ˆ n (x)1 (n) (x)1 (n) (Tˆx) = dxh Sm Sm
ˆ n (x1 , . . . , xn )1 (n) (x)· dx2 . . . dxn h Sm
n n γX γX (1 − γ)T (x1 ) + T (xi ), . . . , (1 − γ)T (xn ) + T (xi ) . n i=1 n i=1
We have now to reduce the domain of integration of I 3 x1 in two steps: the first, 0 changing I into Im , consists in removing intervals of length 2νm on the left and on the
right on each boundary point of the Al , l = 1, · · · , q. Clearly the difference between the 0 will converge to zero when m → ∞ since the integrand functions integrals over I and Im
are bounded (remember the density is in L∞ ), but this argument should be made more precise later on together with the reason of that reduction. For the moment we simply 0 0 write I(I/Im ) for the integral over I/Im . By introducing the operator Pl acting on the
variable xl , l ≥ 2, we could continue as: Z Z (n) ˆ −1 (n) 0 µ ˆ Sm ∩ T Sm = I(I/Im )+ dx1 0 Im
I n−1
h i ˆ n (x1 , . . . , xn )1 (n) (x) · dx2 . . . dxn P2 ◦· · ·◦Pn h Sm
γ γ 1Sm(n) (1 − γ)T (x1 ) + T (x1 ) + x2 + · · · + xn , . . . , (1 − γ)xn + T (x1 ) + x2 + · · · + xn . n n If we now introduce the sets νm νm (n) Sm,γ , j = 2, . . . , n; | xi −xj |≤ i 6= j 6= 1}; (T x1 ) = {(x2 , x3 , · · · , xn ) ∈ I n ; | T (x1 )−xj |≤ 1−γ 1−γ
and (n) Sm (x1 ) = {(x2 , x3 , · · · , xn ) ∈ I n ; |x1 −xj | ≤ νm ; j = 2, · · · , n; | xi −xj |≤ νm
we have
(n) (n) µ ˆ(Sm ∩ Tˆ−1 Sm ) (n)
i 6= j 6= 1},
≤
µ ˆ(Sm ) R
0 Im
dx1
R
h
i 0 ˆ dx2 . . . dxn P2 ◦ · · · ◦ Pn hn (x1 , . . . , xn )1Sm(n) (x) + I(I/Im ) (n) Sm,γ (T x1 ) R R ˆ 1 , · · · , xn ) dx1 Sm dx2 · dxn h(x (n) I 000 (x1 ) m
000 We reduced the domain of integration in the integral in the denominator from I to Im : 0 this kind of reduction will also affect Im and it will be explained in the forthcoming
footnote. Let us now consider for simplicity the structure of the operators when n = 3: ˆ 3 (x1 , x2 , x3 )1 (3) (x1 , x2 , x3 ) = P2 ◦ P 3 h Sm 16
ˆ 3 (x1 , T −1 x2 , T −1 x3 )1 (3) (x1 , T −1 x2 , T −1 x3 ) XX h j j k k S m
j
| DT (Tj−1 x2 ) || DT (Tk−1 x3 ) |
k
1T Aj (x2 )1T Ak (x3 ),
(5.25)
where {Ak } denotes the intervals of monotonicity of the map T . The preceding constraints and the assumption γ < 1 − λ imply that: | Tj−1 x2 − x1 |< νm , | Tk−1 x3 − x1 |< νm ; since the original partition is finite, if we take first m large enough and having removed the intervals of length 2νm around the boundary point of the domain of monotonicity of T , 0 for any x1 ∈ Im there will be only one preimage which can contribute in each sum. By
generalizing to n components we could therefore bound the term (5.25) by λn−1 ||h||∞ . Moreover a simple geometrical inspection shows that the Lebesgue measures of the sets (n)
(n)
Sm,γ (T x1 ) and Sm (x1 ) are independent of the point x1 and also the ratio of the two measures is independent of m and gives15 (n) (n) µ ˆ Sm ∩ Tˆ−1 Sm
(n)
µ ˆ(Sm )
(n)
Leb(Sm,γ ) (n)
Leb(Sm )
=
1 . (1−γ)n−1
We therefore get
(n) ˆ n ||∞ + I(I/I 00 ) Leb(Sm,γ )λn−1 ||h m ≤ . (n) ˆ (inf) Leb(Sm )hn
(5.26)
(n) 00 00 ˆ ∞ Leb(Sm,γ ); ) can be immediately bounded by ||h|| )Leb(I/Im We now notice that I(I/Im (n)
this allows us to factorize the term Leb(Sm,γ ) in the denominator and divide it by (n)
Leb(Sm ). By taking the lim sup we finally get our result.
We can now strengthen the previous result by adding one further assumption; it seems in fact hard to use trace formulas for quasi-H¨older or bounded variation densities in higher dimensions. ˆ is continuous on I n . • P6 The density h This condition is for instance satisfied in the uncoupled case for smooth and locally onto maps T of the unit circle. Proposition 5.3. Let us suppose that our CML verifies properties (P1)-(P6) on a Banach space B with λ = inf |DT |−1 < 1 − γ, then lim
(n) (n) µ ˆ(Sm ∩ Tˆ−1 Sm ) (n)
m→∞ 15This
µ ˆ(Sm )
R hˆ n (x,··· ,x) dx 1 I |DT (x)|n−1 = . R ˆ (x, · · · , x)dx (1 − γ)n−1 h I n
point needs a clarification; take for simplicity n = 2. There is in fact dependence of the two
sets on x1 since they intersect I 2 3 (x2 , x3 ) and as a consequence their measure will depend on the location of x1 . It will be therefore enough to evaluate the external integrals in x1 on a even smaller 00 0 000 domain Im ⊂ Im and on Im in the denominator, in such a way they will not contain a (disconnected)
neighborhood U of 0 and 1 and its preimages T −1 U. As a consequence, we can keep the full amount (2)
(2)
(2)
of the area of the two sets Sm,γ (T x1 ) and Sm (x1 ), which from now on we simply write as Leb(Sm,γ ) (2)
00 000 and Leb(Sm ). Clearly the difference between the integrals over I and Im , Im will converge again to zero (n)
when m → ∞. About the other issue: write Sm,γ (T x1 ) as the integral of obvious characteristic functions in the variables x2 , · · · , xn . Then make the change of variables: x0k = xk (1 − γ) + γT (x1 ); in this way we (n)
get the measure of Sm (x1 ) multiplied by (1 − γ)1−n . 17
Proof. We will write the proof for n = 3, the generalisation being immediate, and this will allows us to use the simple formulas in the previous demonstration. By the same arguments in the latter and by denoting with Tx−1 the inverse branch of T such that 1 Tx−1 (T (x1 )) = x1 , we have 1 (3) (3) Sm ∩ Tˆ−1 Sm
µ ˆ
Z
Z
=
dx1
(3)
00 Im
Sm,γ (T x1 )
ˆ 3 (x1 , T −1 x2 , T −1 x3 ) h x1 x1 00 dx2 dx3 +I(I/Im ) (5.27) −1 x )| |DT (Tx−1 x )||DT (T 2 x1 3 1
(3)
(3)
and we have a lower bound for µ ˆ Sm ∩ Tˆ−1 Sm
00 00 ) ) term. We call I(Im without the I(I/Im
the first integral on the right hand side. ˆ 3 is continuous on I 3 and then uniformly continuous, having fixed εe > 0, it will be Since h x2 −x1 | ≤ νm ), x2 −x1 | ≤ νm ; |Tx−1 enough to choose νm small enough (remember that |Tx−1 1 1 −1 −1 ˆ 3 (x1 , x1 , x1 ) + O(e ˆ 3 (x1 , T x2 , T x3 ) = h ε). to have h x1
x1
For the derivative we can use the fact that our map is C 2 on the interior of the Al , l = 1, · · · , q and extendable with continuity on the boundaries to get by the mean value theorem DT (Tx−1 x2 ) = DT (x1 )+D2 T (ˆ x2 )|Tx−1 x2 −x1 |; DT (Tx−1 x3 ) = DT (x1 )+D2 T (ˆ x3 )|Tx−1 x3 −x1 | 1 1 1 1 where xˆ2 belongs to the interval with endpoints Tx−1 x2 and x1 , and xˆ3 belongs to the 1 interval with endpoints Tx−1 x3 and x1 and these two intervals are in the domains where 1 00 T is locally injective. By replacing in I(Im ) we have: Z ˆ 3 (x1 , x1 , x1 ) Z dx2 dx3 h 00 I(Im ) = dx1 2 T (ˆ 2 D x ) 2 (3) −1 |DT (x1 )| 00 |Tx x2 − x1 |][1 + Im Sm,γ (T x1 ) [1 + DT (x1 )
Z
1 dx1 |DT (x1 )|2 00 Im
Z
D2 T (ˆ x3 ) |Tx−1 x3 1 DT (x1 )
1
O(e ε)
(3) Sm,γ (T x1 )
[1 +
D2 T (ˆ x2 ) |Tx−1 x2 1 DT (x1 )
− x1 |][1 +
D2 T (ˆ x3 ) |Tx−1 x3 1 DT (x1 )
− x1 |]
− x1 |]
dx2 dx3
We now rewrite the first piece as I1,m :=
(3) Leb(Sm,γ )
Z dx1 00 Im
Z
ˆ 3 (x1 , x1 , x1 ) h 1 (3) 2 |DT (x1 )| Leb(Sm,γ ) dx2 dx3
(3)
Sm,γ (T x1 )
[1 +
D2 T (ˆ x2 ) |Tx−1 x2 1 DT (x1 )
− x1 |][1 +
D2 T (ˆ x3 ) |Tx−1 x3 1 DT (x1 )
− x1 |]
(5.28)
where we have suppressed the dependence on T x1 in the Lebesgue measure of the external (3)
00 Sm,γ , which are independent of T x1 when x1 ∈ Im , and the second piece as Z 1 1 (3) I2,m := Leb(Sm,γ ) dx1 (3) |DT (x1 )|2 Leb(Sm,γ 00 ) Im
Z (3) Sm,γ (T x1 )
O(e ε) dx2 dx3 [1 +
D2 T (ˆ x2 ) |Tx−1 x2 1 DT (x1 )
− x1 |][1 + 18
D2 T (ˆ x3 ) |Tx−1 x3 1 DT (x1 )
− x1 |]
+
Using same arguments we have: Z (3) (3) µ ˆ Sm ) = Leb(Sm )
ˆ 3 (x1 , x1 , x1 ) dx1 h
000 Im
(3) Leb(Sm )
Z dx1 000 Im
Z
1 (3) Leb(Sm )
(3) Sm (x1 )
Z
1 (3)
Leb(Sm )
dx2 dx3 +
(3)
Sm (x1 )
000 000 O(e ε)dx2 dx3 +I(I/Im ) = I3,m +I4,m +I(I/Im ) (5.29)
(3)
000 ) term. Hence we get and with a lower bound for µ ˆ Sm ) without the I(I/Im (n) (n) 00 I1,m + I2,m ) µ ˆ(Sm ∩ Tˆ−1 Sm ) I1,m + I2,m + I(I/Im ≤ ≤ (n) 000 I3,m + I4,m + I(I/Im ) I3,m + I4,m µ ˆ(Sm ) (3) 00 00 ˆ ∞ Leb(Sm,γ ) ) ≤ ||h|| )Leb(I/Im As in the proof of Proposition 5.1, we have that I(I/Im (3) ˆ ∞ Leb(Sm )Leb(I/I 000 ). We can then factorize in the numerator and and I(I/I 000 ) ≤ ||h|| m
m
(3)
(3)
in the denominator the Lebesgue measures of the sets Sm,γ and Sm and remember that (3) Leb(Sm,γ ) (3) Leb(Sm )
=
part of I1,m
1 . (1−γ)2
After this factorization and when m goes to infinity, the remaining R ˆ 3 (x1 ,x1 ,x1 ) converges to I dx1 h|DT by the dominated convergence theorem and (x1 )|2
the fact that |Tx−1 xj − x1 | ≤ νm , j = 2, 3, while the remaining part of I2,m converges to 1 a O(e ε) term. Still after the previous factorization, the remaining part of I3,m goes to R ˆ 3 (x1 , x1 , x1 ), while the remaining part of I4,m goes to a O(e dx1 h ε) term. The result I then follows sending εe to zero.
Corollary 5.4. As a consequence of the Proposition, the extremal index θn for maps satisfying P0 too, is given by R hˆ n (x,··· ,x) dx 1 I |DT (x)|n−1 θn = 1 − . R ˆ n (x, · · · , x)dx (1 − γ)n−1 h I
(5.30)
and it will converge to 1 when n → ∞. Very likely, the previous result holds under less stringent assumptions. 5.1. Random perturbations. As for localisation, we expect that the extremal index be one when we keep n fixed and we put additional noise. In the paper [1] we extended the spectral approach to randomly perturbed dynamical systems, mostly with additive noise. Even if we assume properties (P1)-(P5) on some Banach space B, there will be a new difficulty related to the computation of the quantities qk in (2.5) in the random setting. Such a computation as it was done in Proposition 5.3 in [1] strongly relies on the fact that the observable becomes infinite in a single point, the center of a ball: we do not know how to adapt it in the neighborhood of the diagonal ∆. We will present nevertheless numerical evidences in Section 7 that in presence of noise the EI is 1. 19
6. Distribution of the number of successive visits We anticipated in the Introduction that once the synchronisation is turned on for the first time, it cannot last since almost every orbit is recurrent. However the orbit Tˆn (x0 ) will visit for almost every point x0 infinitely often the neighborhood of the diagonal. We could therefore reasonably expect that the exponential law e−τ given by the EVT describes the time between successive events in a Poisson process. To formalize this, let (n)
us take a neighborhood Sς
of the diagonal ∆ with accuracy ac = ς and introduce the following quantity (remember that the map Tˆ and the measure µ ˆ depend on n too): "
# t (n) µ(S ˆ ς )
X
Nς(n) (t) =
1Sς(n) (Tˆl (x))
l=1
and consider the following distribution N (n, ς, t, k) := µ ˆ(Nς(n) (t) = k) If the target set was a ball of radius ς around a generic point z or a dynamical cylinder set converging to this point, one can prove under reasonably mixing assumptions like those invoked in this paper, that in the limit of vanishing radius or infinite length for the cylinder, then N (n, ς, t, k) converges to the Poisson distribution
tk e−t , k!
see for instance
[12],[13]. Instead if we take the target point z periodic of minimal period q, one get the so-called compound Poisson distribution, see [14] and [8], which in our situation reads, for k ≥ 1 : −t(1−p)
N (n, ς, t, k) = e
k X
p
k−j
j j+1 t (1
(1 − p)
j=0
where p=
− p)j k − 1 j! j−1
1 | det(Dz Tˆq )|
(6.31)
(6.32)
Remark 6.1. We do not dispose for the moment of analogous formulae when a ball is replaced by a strip along our diagonal set ∆. At our knowledge the only known result is in dimension 2 for the uncoupled systems given by the direct product of two piecewise expanding and smooth maps of the circle, see [5], and it is consistent with our results. Nevertheless a few preliminary considerations16 seem to indicate that the compound distribution (6.31) still holds with p in (6.32) replaced by 1−θn in (5.30), and more generally with the EI given by formulae (2.4), (2.5), with the quantities qk,l given by the right hand side of (5.23) when transfer operator is not available. In particular one should recover a pure Poisson distribution when the size n of the lattice tends to infinity. Example 6.2. (Ex. (4.1) revisited) 16At
this regard see also the discussion in the last part of Section 7. 20
• Suppose we consider as in the example (4.1), Ex. 2, n = 100 particles living in the unit interval and take the accuracy ς = 0.01. With that value of n and taking the coupling γ sufficiently small, we could consider that the previous number of (n)
visits Nς (t) follows a Poisson distribution. Since the probability of entering the neighborhood of the diagonal is of order 100−100 , the probability to observe exactly 5 synchronization events during m iterations of the lattice is maximal for m = 5 100100 and is of the order of 18%. • If instead we consider Ex. 1 with 3 particles and the same accuracy, the probability to observe 5 synchronisation is maximal after 50 000 iterations and it is again of 18%. Comment 6.3. In the case of large n the extremely high number of iterations needed to get synchronisation or a given number of successive synchronisation could surprise. One reason is surely due to the fact we considered lattices which are globally coupled and we looked at global synchronisation. It would be interesting, and it will be the objects of future investigations, to explore CML where only the nearest-neighbors of a given site contribute to the coupling term (diffusive coupling), and also synchronisation of the closest neighbors. About the latter we will give a few preliminary numerical results in the next section. 7. Extensions and Numerical computations The goal of numerical computations will be to show that in the situations considered above we have effectively convergence toward an extreme value law and moreover the extremal index satisfies the behavior we predicted theoretically. We will be mostly interested in synchronisation, since for localisation we have plenty of analytic results. But there is one aspect where the comparison with localisation is particularly useful. In order to explain that, we have first to introduce a new observable to depict a different kind of synchronisation. 7.1. Local synchronisation. Up to now synchronisation was defined by asking that all the components of the evolutionary state become close to each other with a given accuracy ac . We could ask instead that each component synchronise only with the close neighbors. This is done by introducing the following observable Θ(x) := − log{max |xi − xj |, i 6= j; j = i ± 1}
(7.33)
(of course on the boundary points j will take only one value). We could generalize to more than one neighbors j = i ± 2, ±3, etc., but we limit ourselves here to the case ±1. It is not immediately obvious to have a geometrical description of the set the orbit of a point will visit for the first time (and therefore to give analytical results in terms of the EI), 21
although the ”physical” interpretation will be the same, namely we get the probability that the lattice will have for the first time and after a given number of iteration m, all the components synchronised with the close neighbors and with a given accuracy ac . We call that local synchronisation, to distinguish from the global synchronisation described in the preceding section. It seems intuitive from a physical point of view, that for m large enough and for a given accuracy ac , the probability to get local synchronisation for the first time (from now on we write it as P1 (·) for the different cases), is larger than that to get global synchronisation, P1 (glob. sync.) ≤ P1 (loc. sync.), and this will be confirmed by the numerical simulation as we will see in a moment. On the other hand as soon as the global synchronisation occurs, all the components of the lattice will be aligned in a narrow strip around all of them, and this is close to localisation. Therefore we will expect that the probability to get localisation is larger than the probability of global synchronisation. This is also confirmed by an easy application of the theory. Suppose we fix m and the accuracy ac ; we have also fixed n. By using the scaling formulae in the preceding sections and supposing a pure exponential law for the asymptotic distribution of the maximum, we have 1
n
τ n • For localisation: ac ∼ ( m ) , which gives P1 (local.) ∼ e−τ ∼ e−mac . 1
τ n−1 • For global synchronisation: ac ∼ ( m ) , which gives P1 (glob. sync.) ∼ e−τ ∼ n−1
e−mac . We see that P1 (glob. sync.) ≤ P1 (local.).
7.2. Blocks of synchronisation. The observable (7.33) could be modified further by introducing a new one which we are going to define. Let us first construct N blocks of L successive integer indices: Bq := {iq , · · · , iq + L} and take these blocks disjoint and possibly scattered along the lattice. Then we define: Υ(x) := − log{max |xi − xj |, i 6= j; (i, j) ∈ Bq ; q = 1, · · · , N }. The distribution of the maximum of this observable will give us the probability that the particles in the N blocks will synchronise for the first time with a given accuracy. On the other hand we do not require any synchronisation of the particles outside those blocks. If such a limiting distribution would exist, it could be consistent with the appearance of chimeras in chains of coupled particles, namely patterns of synchronized sets which emerge as a consequence of the self-organization of the entire network, see e.g. [18]. If our claim would be confirmed, such a self-organization would be another statistical property of chaotic systems with several degrees of freedom. 22
7.3. Simulations. Let us now analyse the results of numerical procedure. The experiment performed is the following: we consider the one-dimensional map T in (2.1) as T (x) = 3x mod 1. Once we have constructed the CML Tˆ we will perturb it with additive noise: Tˆω (x)i = Tˆ(x)i + εωi
mod 1
where ε is here the noise intensity and ω with components ωi is a random variable drown from a uniform distribution between -0.5 and +0.5. The stationary measure for such a map will be L1 close to that for γ = 0 which is the direct product of the uniform Lebesgue measures on the unit circle for each component and this independently of the value of ε. Let us notice that we are considering now a one-dimensional map on the circle; this is not a restriction to our previous considerations and moreover it allows us to define correctly the additive noise. Numerically we produce trajectories of 104 iterations for γ < 2/3 and 0.02 increments. The range 3 < n < 53 is analysed. We consider the two observable ψ, see (4.16) and Θ, see (7.33), correspondingly respectively to global and local synchronisation cases and in the following we will refer to them as the global and local cases. We analyze also the role of small noise ε = 10−4 and moderate noise ε = 10−2 . We first assess the convergence of the maxima of ψ and Θ to the Gumbel law by analyzing the tail index ξ, see section 3. Here we chose to consider the complementary approach to the block-maxima selection, i.e. the peak over threshold. The two approaches are equivalent in chaotic systems as shown in [17]. The maxima of the observable are defined as the exceedances over the 0.98 quantile of ψ and Θ distributions. If a good convergence towards the Gumbel law is reached, than ξ ' 0. The values of ξ as a function of γ and n are reported in Fig. 3. A maximum likelihood estimator has been used for computation. The left panels show the global case ψ while the local case Θ is reported on the right. From top to bottom we switch on the noise. In general, the convergence towards the Gumbel law is satisfactory although some differences exist between global and local cases. For the global case the convergence is slower as the global synchronization event is more rare then the local one. Moreover, the quality of the fits is lower when n and γ are larger. The addition of noise helps the convergence to the Gumbel law as for the systems analyzed in [17]. We now study whether global and local synchronisation imply for the extremal index θ. For the analysis presented in this paper, we adopt the estimator by S¨ uveges (see the book [17] for explanation and to retrieve the codes for the computation). For fixed quantile q, S¨ uveges’ estimator reads: 23
PNc i
(1 − q)Si + N + Nc −
θ=
r P 2
Nc i (1
− q)Si + N + Nc
PNc
(1 − q)Si
i
2
− 8N c
PNc i
(1 − q)Si ,
where N is the number of recurrences above the chosen quantile, Nc is the number of observations which form a cluster of at least two consecutive recurrences, and Si the length of each cluster i. From the numerical point of view, this estimator is the expected value of the compound distribution N (n, ς, t, k) with Si being the empirical equivalent of (n)
the quantity Nς (x). We begin by checking the theoretical results predicted in Remark 5.4: for the 3x mod 1 map, θn can now be estimated by taking the trace of the density on the diagonal 1 1 reasonably of order 1 in (5.30), so that in dimension 2: θ2 ∼ 1 − (1−γ) and in dimension 3
3: θ3 ∼ 1 −
1 1 . (1−γ)2 9
The comparison between the theoretical curves and the numerical computations are shown in Fig. 1. For each case n = 2, 3 and γ < 2/3 we produce 10 simulations of the map consisting of 104 iterations and we estimate the extremal index as a function of γ. The numerical estimates indeed match the theoretical curves (bold magenta lines). We now check the asymptotic formula for large n and still with the same assumption λ n−1 ) , with λ = 1/3. For each 3 < n < 53 on the trace of the density, namely θn ∼ 1 − ( 1−γ
and γ < 2/3 we perform one simulation of the deterministic 3x mod 1 map and compare the obtained extremal index θn with the previous asymptotic formula. Results are shown in Figure 2. There is indeed very good agreement between asymptotic and numerical results. The largest divergence is obtained for γ ' 2/3 which correspond to the limit value for the map. We then perform a numerical analysis of the extremal index in the cases not covered by the theory, namely for the observable Θ. The results are presented in Fig. 4. The top-left panel is repeated for convenience and show the global case results. The latter show that global and local cases are substantially different. For the global, the synchronization depends on both n and γ: in particular, it is easier to synchronize systems with n small because the probability of find all the particles in the same state decrease quickly with n. On the other hand, in the local case the extremal index θn is substantially independent on n. In fact, whether n is small or large, the particle sees only the nearest neighbors for synchronization so that is insensitive to the size of the lattice. The only dependence left is in γ: in particular, for all n, we see the dependence is compatible with the case n = 2 1 1 . The addition of the noise destroys clusters as of the global coupling case: θ2 ∼ 1 − (1−γ) 3
24
observed in [1]. Qualitatively, the structure of the extremal index is quite robust with respect to small perturbations. To fully destroy the clusters, large intensity of the noise are needed. The results for ε = 10−4 also demonstrate that our results are stochastic stable because one recover the deterministic structure of the extremal index for low noise values. Although the numerical estimates of the extremal index are done by computing the expected values of the Compound Poisson distribution (S¨ uveges’ estimator), we can also check that the waiting times qk,ς 17 defined in (5.23), between consecutive entrances in the neighborhood of the diagonal with accuracy ς, provides the same information. Actually this is what we get for recurrence in balls as we discussed above, see [14] and [8], and which should just be proved in our situation as well. Therefore we give some examples of time series of ψ and Θ in Figs. 5 and 6 respectively. The noise increases from top to bottom. The histograms of the waiting times in cluster are renormalised to sum-up to 1 (empirical probability density function EPDF) and are in y-log scale. No clustering corresponds to an exponential law (sequence of linearly decreasing boxes in log scale), whereas the clustering case is characterized by an higher EPDF for lower waiting times. As one can see from the deterministic cases, the higher the EPD for short waiting times, the lower θ. Effectively the fraction of waiting times equals to 1 which excess the standard exponential law is exactly the extremal index θ. We stress again that although we cannot demonstrate this identity theoretically, the numerical evidence suggests that one can use directly qk,ς as defined in (5.23), for the estimation of the extremal index θ.
Acknowledgments SV and PG were supported by the MATH AM-Sud Project Physeco. SV was supported by the Leverhulme Trust for support thorough the Network Grant IN-2014-021 and by the project APEX Syst`emes dynamiques: Probabilit´es et Approximation Diophantienne PAD funded by the R´eegion PACA (France). DF was partially supported by the ERC grant A2C2 (No. 338965). PG thanks FONDECYT project 1171427. SV warmly thanks J. M. Freitas, P. Giulietti, N. Haydn and B. Saussol for illuminating discussions. References [1] H. Aytach, J.M. Freitas, S. Vaienti, Laws of rare events for deterministic and random dynamical systems, Trans. Amer Math. Soc., 36 , 8229-8278, 2015. [2] P. Ashwin, Riddled Basis and Coupled Dynamical System, in [4] [3] L.A. Bunimovich, Ya.G. Sinai, Space-time chaos in coupled map lattices, Nonlinearity, 1, 491-519, 1988. 17We
now index this quantity with the size ς → 0 of the neighborhood of the diagonal. 25
[4] J-R. Chazottes, B. Fernandez Editors, Dynamics of Coupled Map Lattices and of Related Spatially Extended Systems, Volume 671 of the series Lecture Notes in Physics pp 265-284, 2005. [5] Z. Coelho, P. Collet,Asymptotic limit law for the close approach of two trajectories in expanding maps of the circle, Probability Theory and Related Fields, 99, 237-250, 1994. [6] D. Azevedo, A.C.M. Freitas, J.M. Freitas, F.B. Rodrigues, Extreme Value Laws for dynamical systems with countable extremal sets, J. Stat. Phys., 167, no. 5, 1244-1261, 2017. [7] A. C.M,Freitas, J. M. Freitas and M. Todd, The extremal index, hitting time statistis and periodicity, Adv. Math., 231(5): 2626-2665, 2012. [8] A. C.M,Freitas, J. M. Freitas and M. Todd , The compound poisson limit ruling periodic extreme behaviour of non-uniformly hyperbolic dynamics, Comm. Math. phys., pages 1-45, 2013. [9] G. Keller, Rare events, exponential hitting times and extremal indices via spectral perturbation, Dyn. Syst. 27, no. 1, 1127 2012. [10] G. Keller, C. Liverani, Rare events, escape rates and quasistationarity: some exact formulae, J. Stat. Phys. 135 (2009), no. 3, 519534. [11] G. Keller, C. Liverani, Stability of the spectrum for transfer operators, Ann. Scuola Norm. Sup. Pisa Cl. Sci.(4), 28:141-152, 1999. [12] M. Hirata, B. Saussol, S. Vaienti,
Statistics of return times: A general framework and new
applications, Communication in Mathematical Physics, 206, 33-55, 1999. [13] N.Haydn, S. Vaienti, The limiting distribution and error terms for return time of dynamical systems, Discrete and Continuous Dynamical Systems A, 10, 584-616, 2004. [14] N. Haydn, S. Vaienti, The compound Poisson distribution and return times in dynamical systems, Probability Theory and Related Fields, 144, 517-542, 2009. [15] K. Kaneko,On the Strength of Attractors in a High-dimensional System: Milnor Attractor Network, Robust Global Attraction, and Noise-induced Selection, Physica D, 124, 322-344, 1998 [16] M. R. Leadbette G. Lindgren and H Rootzn, Extremes and related properties of random sequences and processes, Springer Series in Statistics. Springer-Verlag, New York, 1983. [17] V. Lucarini, D. Faranda, A. M. Freitas, J. M. Freitas, M. Holland, T. Kuna, M. Nicol, M. Todd, S. Vaienti, Extremes and Recurrence in Dynamical Systems, Wiley, New York, 2016, ISBN: 9781-118-63219-2 [18] Y. Sul Cho, T. Nishikawa, A. E. Motter, Stable Chimeras and Independently Synchronizable Clusters, to appear in Phys. Rev. Lett., https://arxiv.org/pdf/1707.06657.pdf. [19] B. Saussol, Absolutely continuous invariant measures for multidimensional expanding maps, Israel J. Math. 116: 223248, 2000,
26
Figure 1. Extremal index θ of the Generalized Pareto distribution as a function of the coupling parameter γ. Thin lines indicate estimates for 10 different realization of the maps 3x-mod1. Bold magenta lines indicate the expected theoretical values. Left: n = 2, Right: n = 3.
27
Figure 2. Extremal index θ as a function of the number of variables n and the coupling parameter γ. Top: global case ψ. Bottom: theoretical asymptotic formula.
28
Figure 3. Shape parameter ξ of the Generalized Pareto distribution as a function of the number of variables n and the coupling parameter γ. Left: global case ψ. Right: local case Θ. From top to bottom: Deterministic, additive noise with intensity ε = 10−4 , additive noise with intensity ε = 10−2 .
29
Figure 4. Extremal index θ as a function of the number of variables n and the coupling parameter γ. Left: global case ψ. Right: local case Θ. From top to bottom: Deterministic, additive noise with intensity ε = 10−4 , additive noise with intensity ε = 10−2 .
30
Figure 5. Example of global case. Right: ψ time series (red) and exceedances (black). Left: empirical probability distribution (EPDF) of waiting time in the clusters. From top to bottom: Deterministic, additive noise with intensity ε = 10−4 , additive noise with intensity ε = 10−2 .
31
Figure 6. Example of local case. Right: Θ time series (red) and exceedances (black). Left: empirical probability distribution (EPDF) of waiting time in the clusters. From top to bottom: Deterministic, additive noise with intensity ε = 10−4 , additive noise with intensity ε = 10−2 .
32
