Xinxin Chen

Waiting times for particles in a branching Brownian motion to reach the rightmost position by

Keywords:

It has been proved by Lalley and Sellke [13] that every particle born in a branching Brownian motion has a descendant reaching the rightmost position at some future time. The main goal of the present paper is to estimate asymptotically as s goes to infinity, the first time that every particle alive at the time s has a descendant reaching the rightmost position.

1 Introduction

The model

We consider a branching Brownian motion (BBM) on the real line R, which evolves as follows. Starting at time t = 0, one particle located at 0, called the root, moves like a standard Brownian motion until an independent exponentially distributed time with parameter 1. At this time it splits into two particles, which, relative to their birth time and position, behave like independent copies of their parent, thus moving like Brownian motions and branching at rate 1 into two copies of themselves. Let N (t) denote the set of all particles alive at time t and let N(t) := #N (t). For any v ∈ N (t) let X v (t) be the position of v at time t; and for any s < t, let X v (s) be the position of the unique ancestor of v that was alive at time s. We define R(t) := max u∈N (t)

X u (t) and L(t) := min u∈N (t)

X u (t), which stand for the rightmost and leftmost positions, respectively. The positions of the extremal particles of a BBM, R(t), have been much studied both analytically and probabilistically. Kolmogorov et al. [START_REF] Kolmogorov | Étude de l'équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique[END_REF] proved that R(t)/t converges almost surely to √ 2. Bramson [START_REF] Bramson | Maximal displacement of branching Brownian motion[END_REF] [START_REF] Bramson | Convergence of solutions of the Kolmogorov equation to travelling waves[END_REF] showed that R(t)-√ 2t+(3/2 √ 2) log t converges in law. These results hold as well for a wide class of branching random walks under mild conditions: see for example Biggins [START_REF] Biggins | The first-and last-birth problems for a multitype age-dependent branching process[END_REF], Addario-Berry and Reed [START_REF] Addario-Berry | Minima in branching random walks[END_REF], Hu and Shi [START_REF] Hu | Minimal position and critical martingale convergence in branching random walks, and directed polymers on disordered trees[END_REF], Aïdékon [START_REF] Aïdékon | Convergence in law of the minimum of a branching random walk[END_REF]. In particular, we state the following fact, which is first given by Hu and Shi [START_REF] Hu | Minimal position and critical martingale convergence in branching random walks, and directed polymers on disordered trees[END_REF] for branching random walks, and is recently proved by Roberts [START_REF] Roberts | A simple path to asymptotics for the frontier of a branching Brownian motion[END_REF]:

lim inf t→∞ R(t) - √ 2t log t = - 3 2 √ 2 almost surely; (1.1) lim sup t→∞ R(t) - √ 2t log t = - 1 2 √ 2 almost surely. (1.2)
In [START_REF] Lalley | A Conditional Limit Theorem for the Frontier of a Branching Brownian Motion[END_REF], Lalley and Sellke showed the following interesting property: every particle born in a BBM has a descendant reaching the rightmost position at some future time. Such a particle was thought of having a prominent descendant "in the lead" at this time. This property is in agreement with the branching-selection particle systems investigated in [START_REF] Brunet | Shift in the velocity of a front due to a cutoff[END_REF], [START_REF] Brunet | Microscopic models of traveling wave equations[END_REF] and [START_REF] Bérard | Brunet-Derrida Behavior of branching-selection particle systems on the line[END_REF]. These papers bring out the fact that the extremal positions of a branching system on the line cannot always be occupied by the descendants of some "elite" particles.

In the present work, we give some quantitative understanding of this behavior, and precisely speaking, about how long we have to wait so that every particle alive at time s has a descendent that has occupied the rightmost position.

The main problem

Let us make an analytic presentation for our problem. For any s > 0 and each particle u ∈ N (s), the shifted subtree generated by u is

(1.3) N u (t) := v ∈ N (t + s), u ≤ v , ∀t ≥ 0,
where u ≤ v indicates that v is a descendant of u or is u itself. Further, for any v ∈ N u (t), let

(1.4) X u v (t) := X v (t + s) -X u (s), be its shifted position. We set R u (t) := max v∈N u (t) X u v (t) and L u (t) := min v∈N u (t) X u v (t). Moreover, Let F t ; t ≥ 0 be the natural filtration of the branching Brownian motion. The branching property implies that, given F s , {R u (•); u ∈ N (s)} are independent copies of R(•). Moreover, we denote by F u ∞ the sigma-field generated by the shifted subtree started from the time s rooted at u.

For every u ∈ N (s), let

(1.5) τ u := inf{t > 0 : R(t + s) = X u (s) + R u (t)}.

The random variable τ u stands for the first time that started from time s, the particle u has a descendant reaching the rightmost position in the system. It is the object in which we are interested. We define

(1.6) Θ s := max u∈N (s) τ u ,
which represents the first time when every particle in N (s) has had a descendant occupying the rightmost position.

According to Lalley and Sellke [START_REF] Lalley | A Conditional Limit Theorem for the Frontier of a Branching Brownian Motion[END_REF], for any s > 0, P[Θ s < ∞] = 1. Since Θ s → ∞ almost surely as s → ∞, we intend to determine the rate at which Θ s increases to infinity.

The main results

To estimate Θ s = max u∈N (s) τ u , an intuitive idea consists in saying that, the further a particle is away from the rightmost one, the longer it has to wait for a descendant to be located on the rightmost position. We thus first focus on the leftmost particle. Let ℓ(s) be the leftmost particle alive at time s. By (1.5), τ ℓ(s) is defined as the shortest time needed for ℓ(s) to wait to have a descendant occupying the rightmost position.

Theorem 1.1 The following convergence holds almost surely

(1.7) lim s→∞ log τ ℓ(s) s = 4.
However, the leftmost particle is not the one who "drags the feet" of the whole population N (s). By considering the positions of all particles alive at time s, as well as their evolutions, we obtain our main result as follows.

Theorem 1.2 The following convergence holds almost surely

(1.8) lim s→∞ log Θ s s = 2 + 2 √ 2 > 4.
Remark 1.3 The proof of the theorems will reveal that the largest τ u for u ∈ N (s) is achieved by some particle located at a position around -(2 -√ 2)s which does not split until time s + 1 √ 2 s and moves towards to the left as far as possible.

The rest of this paper is organized as follows. Section 2 is devoted to discussing the behaviors of the extremal position R(•), which leads to two propositions. In Section 3, we consider the case of two independent branching Brownian motions and state another proposition. We prove Theorem 1.1 in Section 4 by means of these propositions. Finally, in Section 5, we prove Theorem 1.2.

The behavior of the rightmost position

We recall Proposition 3 in Bramson's work [START_REF] Bramson | Maximal displacement of branching Brownian motion[END_REF]. It is shown that for all 0 ≤ y ≤ t 1/2 and t ≥ 2, there exists a positive constant c which is independent of t and y, such that

(2.1) P R(t) > m(t) + y ≤ c(1 + y) 2 exp(- √ 2y),
where

m(t) := √ 2t - 3 2 √ 2 log t.
Therefore, with c 1 := c + 1, we get the following inequality, which will be applied several times in our arguments.

Fact 2.1 (Bramson [6]) For any t ≥ 2 and y ≤ √ t,

(2.2) P R(t) > m(t) + y ≤ c 1 1 + y + 2 e - √ 2y ,
with y + := max{y, 0}. Let (B s ; s ≥ 0) be a standard Brownian motion on R. We state the following lemma, which can be found in several papers (e.g. [15] [10]). It is also of frequent use.

Lemma 2.2 (many-to-one) For any measurable function F and each t > 0,

(2.3) E u∈N (t) F (X u (s), s ∈ [0, t]) = e t E F (B s , s ∈ [0, t]) ,
where, for each u ∈ N (t) and s ∈ [0, t], X u (s) denotes the position, at time s, of the ancestor of u.

Let us present the following inequality as well, which is Equation (57) in Bramson [START_REF] Bramson | Maximal displacement of branching Brownian motion[END_REF].

Fact 2.3 For any s ≥ 1 and any a > 0,

1 - s a 2 s 2π a -1 exp - a 2 2s ≤ P[B s ≥ a] ≤ s 2π a -1 exp - a 2 2s .
It immediately follows that

(2.4) P[B s ≤ -a] = P[B s ≥ a] ≤ √ s a exp - a 2 2s .
Moreover, if a = αs with some constant α > 0, we have

(2.5) P[B s ≤ -αs] = P[B s ≥ αs] = exp - α 2 2 + o s (1) s ,
where o s (1) → 0 as s goes to infinity. We define, for any y > 0,

T (y) := inf t ≥ 1; R(t) -m(t) > y .
Because of (1.2), one immediately sees that P T (y) < ∞ = 1 for any y > 0. Moreover, T (y) ↑ ∞ almost surely as y ↑ ∞.

Proposition 2.4

The following convergence holds almost surely

(2.6) lim y→∞ log T (y) y = √ 2.
Proof: First, we prove the lower bound. Let 2 ≤ y ≤ √ t, and set

Λ := E t+1 1 1 (R(s)>m(s)+y-1) ds .
Clearly, Λ = hold. Plugging them into the integration of (2.7) yields that (2.9)

t+1 1 P R(s) > m(s) + y -1 ds. Hence, Λ = y 2 1 P R(s) > m(s) + y -1 ds + t+1 y 2 P R(s) > m(s) + y -1 ds ≤ y 2 1 E u∈N (s)
y 2 1 E u∈N (s)
P R(s) > m(s) + y -1 ds ≤ t+1 y 2 c 1 y 2 e - √ 2y+ √ 2 ds ≤ c 1 ty 2 e - √ 2y+ √ 2 .
Combining (2.9) with (2.10), we have with c 3 > 0 a constant independent of (y, t).

On the other hand,

Λ ≥ E t+1 1 1 (R(s)>m(s)+y-1) ds; T (y) ≤ t = t 1 P T (y) ∈ dr E t+1 1 1 (R(s)>m(s)+y-1) ds T (y) = r .
Conditionally on the event {T (y) = r ≤ t}, the rightmost particle in N (r), denoted by ω, is located at m(r) + y. Started from the time r, ω moves according to a Brownian motion and splits into two after an exponential time. By ignoring its branches, we observe that

R(s + r) > [m(r) + y] + [ √ 2s -1] ≥ m(s + r) + y -
ds ≥ st t+1-r 0 1 Bs> √ 2s-1 ds ≥ min{1, T (- √ 2) -1 },
where ≥ st denotes stochastic dominance and T

(- √ 2) -1 := inf{t ≥ 0; B t < √ 2t -1}.
These arguments imply that (2.12) Λ ≥

t 1 P[T (y) ∈ dr]E min{1, T (- √ 2) -1 } =: c 13 P[T (y) ≤ t],
where c 13 := E min{1, T

(- √ 2) -1
} ∈ (0, ∞). Compared with (2.11), this tells us that (2.13)

P T (y) ≤ t ≤ c 5 ty 2 e - √ 2y , for 2 ≤ y ≤ √ t,
where c 5 := c 3 c 13 ∈ (0, ∞). Taking t = e √ 2y(1-δ) with δ ∈ (0, 1) yields that

∞ k=1 P T (k) ≤ e √ 2k(1-δ) < ∞.
According to the Borel-Cantelli lemma, lim inf

y→∞ log T (y) y ≥ √ 2, almost surely,
proving the lower bound in the proposition.

To prove the upper bound, we recall that

R u (t) = max{X u v (t); v ∈ N u (t)}, u ∈ N (s).
Obviously, R u (t); u ∈ N (s) are i.i.d. given F s , and are distributed as R(t).

We fix a y ∈ (0, y) and define the measurable events 1+δ) .

Σ 1 := L(a y ) ≥ -2a y ; N(a y ) ≥ exp 1 2 a y , Σ := Σ 1 ∩ T (y) > e √ 2y (
Then, (2.14)

P T (y) > e √ 2y(1+δ) ≤ P Σ c 1 + P Σ 1 ∩ T (y) > e √ 2y (1+δ)
≤ P N(a y ) ≤ exp 1 2 a y + P L(a y ) ≤ -2a y + P Σ .

We .

At this stage, it is convenient to recall the proof of Proposition 15 of Roberts [START_REF] Roberts | A simple path to asymptotics for the frontier of a branching Brownian motion[END_REF], saying that there exists a constant c ′ > 0 such that for y large enough,

P ∃r ∈ e √ 2y(1+ 1 2 δ) , 2e √ 2y(1+ 1 2 δ) : R(r) ≥ m(r) + 1 √ 2 log r > c ′ > 0.
Thus, P Σ ≤ (1c ′ ) e ay/2 ≤ exp(-c ′ e δ 1 y/2 ).

It remains to estimate P N(a y ) ≤ exp( 1 2 a y ) and P L(a y ) ≤ -2a y . On the one hand, the branching mechanism tells us that for any s ≥ 0, N(s) follows the geometric distribution with parameter e -s (for example, see Page 324 of [START_REF] Mckean | Application of Brownian motion to the equation of Kolomogrov-Petrovskii-Piskunov[END_REF]). It thus yields that P N(a y ) ≤ exp( 1 2 a y ) ≤ e -δ 1 y/2 . On the other hand, as shown in Proposition 1 of Lalley and Sellke [START_REF] Lalley | Travelling waves in inhomogeneous branching Brownian motions II[END_REF], for any µ ≥ √ 2 and s > 0,

(2.15)

P L(s) ≤ -µs = P R(s) ≥ µs ≤ µ -1 (2πs) -1/2 exp -s µ 2 2 -1 .
Consequently, (2.14) becomes that

P T (y) > e √ 2y (1+δ)
≤ e -δ 1 y/2 + e -δ 1 y + exp(-c ′ e δ 1 y/2 ) (2.16)

≤ c 6 e -δ 1 y/2 .
By the Borel-Cantelli lemma again, we conclude that almost surely lim sup

y→∞ log T (y) y ≤ √ 2,
which completes the proof of the proposition. For α > 0 and β > 0, set p(z, α, β) := P ∃r ≤ e αz : R(r) ≤ m(r)βz .

Proposition 2.5 There exists a positive constant C 1 , independent of (α, β, z), such that for any z ≥ z(α, β),

(2.17) p(z, α, β) ≤ C 1 exp - βz 6 √ 2 .
Proof:

It follows from (1.1) that as z → ∞, p(z, α, β) = P ∃r ≤ e αz : R(r) ≤ m(r) -βz -→ 0.
Hence, there exists z 0 (α, β) large enough, such that for all z ≥ z 0 (α, β),

(2.18) P ∃r ≤ e αz : R(r) ≤ m(r) -βz/2 ≤ 1/2.
For any b z < e αz , we have

(2.19) p(z, α, β) ≤ P ∃u ∈ N (b z ), s.t. min s≤bz X u (s) ≤ √ 2b z -βz/2 + P ∃r ≤ e αz : R(r) ≤ m(r) -βz ∩ L(b z ) ≥ √ 2b z -βz/2 .
On the one hand, by the many-to-one lemma,

P ∃u ∈ N (b z ), s.t. min s≤bz X u (s) ≤ √ 2b z -βz/2 ≤ E u∈N (bz) 1 (min s≤bz Xu(s)≤ √ 2bz-βz/2) = e bz P min s≤bz B s ≤ √ 2b z -βz/2 .
On the other hand, by simple observations,

P ∃r ≤ e αz : R(r) ≤ m(r) -βz ∩ L(b z ) ≥ √ 2b z -βz/2 ≤ P u∈N (bz) ∃t ≤ e αz , s.t. R u (t) < m(t) -βz/2 = E u∈N (bz) P ∃t ≤ e αz , s.t. R(t) < m(t) -βz/2 ,
where the last equality follows from the branching property. Going back to (2.19), one has

p(z, α, β) ≤ e bz P min s≤bz B s ≤ √ 2b z -βz/2 +E u∈N (bz) P ∃t ≤ e αz , s.t. R(t) < m(t)-βz/2 . Let b z = β 6 √ 2 z. Then, by (2.18), for all z ≥ z(α, β) := max{z 0 (α, β), 1 β }, p(z, α, β) ≤ e bz P min s≤bz B s ≤ -βz/3 + E 1 2 N (bz ) ≤ c 7 e -3bz + e -bz ≤ C 1 exp - βz 6 √ 2 ,
with C 1 := c 7 + 1, which completes the proof of the proposition.

Corollary 2.6 For any δ ∈ (0, 1), there exists some s(δ) ≥ 1, such that for all s ≥ s(δ),

(2.20) P R(s) ≤ √ 2(1 -δ)s ≤ C 1 exp - δs 12 √ 2 .
Proof: Since we always have m(s)

-δs/2 ≥ √ 2(1 -δ)s when s is sufficiently large, P R(s) ≤ √ 2(1 -δ)s ≤ P ∃r ≤ e s : R(r) ≤ m(r) -δs/2 .
which by Proposition 2.5 is bounded by

C 1 exp -δs 12 √
2 for all s large enough.

The case of two independent branching Brownian motions

We consider two independent branching Brownian motions, denoted by X A (•) and X B (•).

Suppose that P[X A (0) = 0] = P[X B (0) = z] = 1 with z > 0, where X A (0) and X B (0) represent the position of the roots, respectively. We write R A (•) (R B (•), respectively) for the position of rightmost particle of the BBM X A (•) (X B (•), respectively). We define, for any y > 0,

T A (y) := inf{t ≥ 1; R A (t) > m(t) + y}; T B (y) := inf{t ≥ 1; R B (t) > m(t) + y}.
Let T A>B be the first time when the rightmost point of X A exceeds that of X B , i.e.,

T A>B = T A>B (z) := inf{t ≥ 0; R A (t) > R B (t)}.
We immediately observe that the distribution of T A>B (z) merely depends on the parameter z. Actually, we can take another pair of independent standard BBM's (both rooted at the origin), namely, X I (•) and X II (•). Their rightmost positions are denoted by R I (•) and R II (•), respectively. For any positive z, let

T (z) := inf{t ≥ 0 : R I (t) -R II (t) > z}.
Then T A>B (z) is distributed as T (z). Besides, z → T (z) is increasing.

Proposition 3.1 The following convergence holds almost surely

(3.1) lim z→∞ log T (z) z = √ 2.
Proof: For any δ ∈ (0, 1),

P T (z) ≤ e √ 2z(1-δ) = P T A>B (z) ≤ e √ 2z(1-δ) ≤ p 1 + p 2 ,
where

p 1 := P ∃t ≤ e √ 2z(1-δ) , s.t. R B (t) < m(t) + z -δz/2 , p 2 := P T A>B ≤ e √ 2z(1-δ) ∩ R B (t) ≥ m(t) + (1 -δ/2)z, ∀t ≤ e √ 2z (1-δ)
.

Clearly,

p 1 = P ∃t ≤ e √ 2z(1-δ) , s.t. R(t) < m(t) -δz/2 = p(z, √ 2(1 -δ), δ/2)
. By Proposition 2.5, for all z ≥ z(δ),

p 1 ≤ C 1 exp - δz 12 √ 2 .
At the same time, we notice that

(3.2) T A>B ≤ e √ 2z(1-δ) ∩ R B (t) ≥ m(t) + (1 -δ/2)z, ∀t ≤ e √ 2z(1-δ) ⊂ ∃t ≤ e √ 2z(1-δ) : R A (t) ≥ R B (t) ≥ m(t) + (1 -δ/2)z ⊂ T A (1 -δ/2)z ≤ e √ 2z(1-δ) .
This yields that

p 2 ≤ P T (1 -δ/2)z ≤ e √ 2z(1-δ) ≤ c 5 z 2 e -δz/ √ 2 ,
because of the inequality (2.13).

As a result,

(3.3) P T (z) ≤ e √ 2z(1-δ) ≤ C 1 exp - δz 12 √ 2 + c 5 z 2 e -δz/ √ 2 ≤ c 8 exp - δz 12 √ 2 ,
for some constant c 8 > 0 and all z large enough. Thus, by the Borel-Cantelli lemma, lim inf

z→∞ log T (z) z ≥ √ 2 almost surely.
To prove the upper bound, we observe that

(3.4) P T (z) > e √ 2z(1+δ) = P T A>B (z) > e √ 2z(1+δ) ≤ q 1 + q 2 ,
where

q 1 := P T A z(1 + δ/2) > e √ 2z(1+δ) ∪ T A z(1 + δ/2) < e √ 2z
,

q 2 := P e √ 2z ≤ T A z(1 + δ/2) ≤ e √ 2z(1+δ) < T A>B (z) .
Notice that T A (y) is distributed as T (y) for any y > 0. According to the inequalities (2.13) and (2.16), there exists δ 2 := δ 2 (δ) > 0 such that q 1 ≤ e -δ 2 z for z large enough. It remains to estimate q 2 :

q 2 ≤ e √ 2z (1+δ)
e

√ 2z P T A z(1 + δ/2) ∈ dr P T A>B > r T A z(1 + δ/2) = r ≤ e √ 2z (1+δ)
e

√ 2z P T A z(1 + δ/2) ∈ dr P R B (r) > m(r) + z(1 + δ/2) .
By the inequality (2.2) again, this tells that

q 2 ≤ e √ 2z(1+δ) e √ 2z P T A z(1 + δ/2) ∈ dr c 2 (z + 1) 2 e - √ 2δz/2 ≤ c 2 (z + 1) 2 e - √ 2δz/2 .
Thus, recalling (3.4), we obtain that for all z large enough,

P T (z) > e √ 2z(1+δ) ≤ e -δ 2 z + c 2 (z + 1) 2 e - √ 2δz/2 .
It follows that almost surely lim sup z→∞

log T (z) z ≤ √ 2. Proposition 3.1 is proved.
4 Proof of Theorem 1.1

For any k ∈ N + and δ ∈ (0, 1/20), we define

N δ (k) := {u ∈ N (k) : X u (k) ≤ - √ 2(1 -δ/2)k}.
In order to study the asymptotic behavior of τ ℓ(s) for s ∈ R + , we first look for a lower bound for min u∈N δ (k) τ u and an upper bound for max u∈N δ (k) τ u .

Recall the definitions (1.3) and (1.4) of the shifted subtrees. For any particle u ∈ N δ (k), we use X u (•) to represent the branching Brownian motion generated by u started from the time k. Meanwhile, we use X r (•) to represent the branching Brownian motion generated by the rightmost point at time k. Accordingly, the random variable T u>r is defined to be the first time when u has a descendant exceeding all descendants of the rightmost particle at time k.

Considering that T u>r ≤ τ u for each u ∈ N δ (k), one sees that

P u∈N δ (k) {τ u ≤ e 4k(1-10δ) } ≤ p ′ 1 + p ′ 2 ,
where

p ′ 1 := P R(k) ≤ √ 2(1 -δ/2)k , p ′ 2 := P 1 R(k)≥ √ 2(1-δ/2)k u∈N δ (k) 1 T u>r ≤e 4k(1-10δ)
.

Given F k , the BBM's X u and X r are independent. Then,

p ′ 2 ≤ E u∈N δ (k) 1 R(k)≥ √ 2(1-δ/2)k P T u>r ≤ e 4k(1-10δ) F k = E u∈N δ (k) 1 R(k)≥ √ 2(1-δ/2)k P T (R(k) -X u (k)) ≤ e 4k(1-10δ) F k .
By the monotonicity of T (•), this gives that

p ′ 2 ≤ E u∈N δ (k) P T (2 √ 2k(1 -δ/2)) ≤ e 4k(1-10δ) = E u∈N δ (k) 1 P T (2 √ 2k(1 -δ/2)) ≤ e 4k(1-10δ) .
Using the inequality (3.3), for all k sufficiently large,

p ′ 2 ≤ E u∈N δ (k) 1 c 8 exp - 3δk 2 (1 -δ/2) .
Then by the many-to-one lemma and by (2.4), we obtain that (4.1)

p ′ 2 ≤ e k P[B k ≤ - √ 2k(1 -δ/2)]c 8 exp -3δk 2 (1 -δ/2) ≤ e -c 9 δk ,
where c 9 is a positive constant independent of (δ, k).

In view of Corollary 2.6, for large k, one has

(4.2) p ′ 1 ≤ C 1 exp - δk 24 √ 2 .
Combining (4.2) with (4.1) yields that for k large enough,

P u∈N δ (k) {τ u ≤ e 4k(1-10δ) } ≤ C 1 exp - δk 24 √ 2 + e -c 9 δk .
By the Borel-Cantelli lemma, almost surely,

(4.3) lim inf k→∞ log min u∈N δ (k) τ u k ≥ 4(1 -10δ),
which gives the lower bound for min u∈N δ (k) τ u .

To obtain an upper bound for max u∈N δ (k) τ u , let us estimate P ∪ u∈N δ (k) τ u ≥ e 4k(1+10δ) . We consider the subtree generated by any particle u ∈ N δ (k). Recall that the shifted positions of its descendants are denoted by

X u v (•) := X v (• + k) -X u (k) for any v ∈ N (• + k) satisfying u < v,
and that R u (•) := max X u v (•). We set T u (y) := inf{t ≥ 1; R u (t)m(t) > y} for any y > 0, which is obviously distributed as T (y). Let y = 2 √ 2k(1 + δ/2), then (4.4)

P u∈N δ (k) τ u ≥ e 4k(1+10δ) ≤ q ′ 1 + q ′ 2 + q ′ 3 ,
where

q ′ 1 := P u∈N δ (k) T u (y) ≥ e 4k(1+10δ) ∪ T u (y) ≤ e k , q ′ 2 := P L(k) ≤ - √ 2k , q ′ 3 := P u∈N δ (k) {e k < T u (y) < e 4k(1+10δ) ≤ τ u }; L(k) > - √ 2k .
First, we observe that

q ′ 1 ≤ E u∈N δ (k) 1 P T (y) ≥ e 4k(1+10δ) ∪ T (y) ≤ e k .
Using the many-to-one lemma for the first term on the right-hand side,

q ′ 1 ≤ e k P B k ≤ - √ 2k(1 -δ/2) P T (y) ≥ e 4k(1+10δ) ∪ T (y) ≤ e k .
According to the inequalities (2.4) (2.13) and (2.16), there exists

δ 4 := δ 4 (δ) > 0 such that q ′ 1 ≤ e -δ 4 k for k large enough. Meanwhile, by (2.2), q ′ 2 ≤ 2c 2 (log k + 1) 2 k -3/2 . It remains to bound q ′ 3 . Since T u (y) is independent of F k , it follows that q ′ 3 ≤ E   u∈N δ (k) e 4k(1+10δ) e k P T u (y) ∈ dr P τ u > r; L(k) ≥ - √ 2k T u (y) = r, F k   .
Given {T u (y) = r} and F k , the event

{τ u > r}∩{L(k) ≥ - √ 2k} implies that ∪ w∈N (k)\{u} {R ω (r)+ X w (k) > R u (r)+X u (k) ≥ m(r)+y - √ 2k}, whose probability is less than w∈N (k)\{u} c 1 1+ (y - √ 2k -X w (k)) 2 + e - √ 2y+2k+ √ 2Xw(k) (see (2.2)
). This yields that k) .

q ′ 3 ≤ E   u∈N δ (k) e 4k(1+10δ) e k P T u (y) ∈ dr w∈N (k)\{u} c 2 (y + 1) 2 e - √ 2y+2k+ √ 2Xw(k)   ≤ E   u∈N δ (k) w∈N (k)\{u} c 2 (y + 1) 2 e - √ 2y+2k+ √ 2Xω(k)   = c 2 (y + 1) 2 e - √ 2y+2k E u∈N δ (k) w∈N (k)\{u} e √ 2Xw (
By integrating with respect to the last time at which the most recent common ancestor of u and ω was alive (see e.g. [START_REF] Harris | The many-to-few lemma and multiple spines[END_REF] for more details),

E u∈N δ (k) ω =u e √ 2Xω(k) is equal to 2 k 0 e 2k-s ds R P B s ∈ dx P B k ≤ - √ 2(1 -δ/2)k B s = x E e √ 2B k B s = x = 2 k 0 e 2k-s ds R P B s ∈ dx P B k ≤ - √ 2(1 -δ/2)k B s = x e √ 2x e k-s ,
where the second equivalence follows from the Markov property of Brownian Motion. We rearrange the integration as follows:

E u∈N δ (k) w∈N (k)\{u} e √ 2Xw(k) = 2 k 0 e 3k-2s E e √ 2Bs ; B k ≤ - √ 2(1 -δ/2)k ds = 2 k 0 e 3k-2s E e √ 2B k e - √ 2(B k -Bs) ; B k ≤ - √ 2(1 -δ/2)k ds ≤ 2 k 0 e 3k-2s e -2(1-δ/2)k E e - √ 2(B k -Bs) ds,
which is bounded by e (2+δ)k by simple computation. Thus, q ′ 3 ≤ c 10 k 2 e -δk for some constant c 10 > 0.

Going back to (4.4),

P u∈N δ (k) τ u ≥ e 4k(1+10δ) ≤ e -δ 4 k + 2c 2 (log k + 1) 2 k -3/2 + c 10 k 2 e -δk ,
for all k sufficiently large.

Therefore, by the Borel-Cantelli lemma, lim sup

N∋k→∞ log max u∈N δ (k) τ u k ≤ 4(1 + 10δ) almost surely.
We now turn to study {τ ℓ(s) ; s ≥ 0}. On the one hand, for any δ > 0, we claim that almost surely for s large enough, the leftmost particle ℓ(s) at time s must have at least one descendant belonging to N δ (⌊s⌋ + 1). In fact, let us write

Υ k := {∃u ∈ N (k + 1) : u / ∈ N δ (k + 1); ∃s ∈ [k, k + 1], X u (s) ≤ - √ 2s + δ ′ s} with δ ′ := ( √ 2 -1)δ/2
. By the many-to-one lemma, we get that for k ≥ 100/δ,

P[Υ k ] ≤ 4 δk √ 2π e 1+k-δ 2 k 2 8 ,
which is summable over k. It follows that (4.5) P[Υ k infinitely often] = 0.

In view of (1.2), when s is large enough, L(s) always lie below -√ 2s + δ ′ s almost surely. Combining with (4.5), we obtain that almost surely for k sufficiently large, max

s∈[k,k+1] τ ℓ(s) ≤ max u∈N δ (k+1) τ u + 1.
On the other hand, using similar arguments, one can say that almost surely for s sufficiently large, the leftmost located particle ℓ(s) at time s must come from one particle in N δ (⌊s⌋). This gives that almost surely for k sufficiently large,

The lower bound of Theorem 1.2

This subsection is devoted to checking that: almost surely, lim inf

k→∞ log Θ k k ≥ 2 + 2 √ 2.
For 0 < a < √ 2, we define

Z a (k) := u ∈ N (k); X u (k) ≤ -ak and Z a (k) := #Z a (k).
For 0 < ε < (1 -a 2 2 )/2 and 0 < δ < 1, we denote

E k := Z a (k) ≥ exp[k(1 - a 2 2 -ε)] , D k := Θ k ≤ exp[(2 + 2 √ 2 -δ)k] . Let us estimate P[D k ∩ E k ].
For any s > 0 and β > 0, we write Γ = Γ(s, β) := N(s) = 1, L(s) ≤ -βs . Similarly, let Γ u := N u (s) = 1, L u (s) ≤ -βs for every u ∈ N (k). Then, (5.1)

P D k ∩ E k ≤ P u∈Za(k) Γ c u ∩ E k + P u∈Za(k) Γ u ∩ D k .
By the branching structure, we obtain that (5.2)

P u∈Za(k) Γ c u ∩ E k ≤ P 1 -P[Γ] Za(k) ; E k ≤ e -P[Γ] exp[k(1-a 2 2 -ε)] .
Clearly, P[Γ] = e -s P[B s ≤ -βs]. By (2.5), one sees that, for ε > 0 small and s large enough, (5.3)

P u∈Za(k) Γ c u ∩ E k ≤ exp -exp -s(1 + β 2 2 + ε) + k(1 - a 2 2 -ε) ,
which is bounded by e -e kε if we choose s =

1-a 2 2 -2ε 1+ β 2 2 +ε
k with k sufficiently large.

It remains to bound Ω :

= P u∈Za(k) Γ u ∩ D k for s = 1-a 2 2 -2ε 1+ β 2 2 +ε
k. Recalling the definition of Θ k , one sees that for any ρ ∈ (0, 2),

(5.4) Ω ≤ P u∈Za(k)

τ u < e ρk + P u∈Za(k) Γ u ∩ u∈Za(k) e ρk ≤ τ u ≤ e (2+2 √ 2-δ)k =: Ω a + Ω b .
We choose now ρ = 1 -2ε and z = ( In view of the inequalities (2.4) and (2.13), one immediately has (5.6)

√ 2 -a 2 2 √ 2 -ε √
P u∈Za(k) {T u (z) < e ρk } ≤ e k √ k ak e -a 2 k/2 c 5 e ρk z 2 e - √ 2z ≤ e -ηk ,
for some η := η(ε) > 0 small enough.

For the second term of the right-hand side in (5.5), we observe that for any u ∈ Z a (k), {τ u < e ρk ≤ T u (z)} implies that at time τ u < e ρk , the rightmost position R(k + τ u ) is exactly equal to R u (τ u ) + X u (k), which is less than m(τ u ) + zak. Hence, the event ∪ u∈Za(k) {τ u < e ρk ≤ T u (z)} ensures that there exists some time r < e ρk such that the rightmost position R(k + r) is less than m(r) + zak. This gives that (5.7)

P u∈Za(k) {τ u < e ρk ≤ T u (z)} ≤ P ∃r ≤ e ρk s.t. R(k + r) ≤ m(r) + z -ak .
Notice that with our choice of ρ and z, Proposition 2.5 can be applied to show that for all k sufficiently large, (5.8)

P u∈Za(k) {τ u < e ρk ≤ T u (z)} ≤ e -ηk .
Combined with (5.6), the inequality (5.5) becomes Ω a ≤ 2e -ηk .

As shown in (5.4), it remains to study Ω b . For the particles u ∈ Z a (k) such that N u (s) = 1, we focus on the subtree rooted at u but started from time k + s. Define

R u (t) := max X v (k + s + t) -X u (k + s); v ∈ N (k + s + t), u < v , ∀t ≥ 0; and T u (y) := inf t ≥ 1; R u (t) ≥ m(t) + y , ∀y > 1.
Since ( R u (t), t ≥ 0) is distributed as (R(t), t ≥ 0), T u (y) has the same law as T (y). Let us take

√ 2x = k(2 + 2 √ 2 -δ/2). Comparing T u (x) with e k(2+2 √ 2-δ) yields that (5.9) Ω b ≤ P ∃ ω ∈ Z a (k), s.t. N ω (s) = 1, L ω (s) ≤ -βs, T ω (x) ≤ e k(2+2 √ 2-δ) +P ∃u ∈ Z a (k) s.t. N u (s) = 1, L u (s) ≤ -βs, e (1-2ε)k ≤ τ u ≤ e k(2+2 √ 2-δ) < T u (x) =: Ω b1 +Ω b2 .
By first conditioning on F k+s and then on F k , one has

Ω b1 ≤ E Z a (k) P Γ P T (x) ≤ e k(2+2 √ 2-δ) (5.10) = e k P B k ≤ -ak e -s P B s ≤ -βs P T (x) ≤ e k(2+2 √ 2-δ) .
For 0 < ε < min{δ/8, (1 -a 2 2 )/2}, by (2.13) and (2.4),

Ω b1 ≤ e 3εk c 5 x 2 e - √ 2x e k(2+2 √ 2-δ) ≤ c 11 k 2 e -εk . (5.11) On the other hand, the event {∃u ∈ Z a (k) s.t. N u (s) = 1, L u (s) ≤ -βs, e (1-2ε)k ≤ τ u ≤ e k(2+2 √ 2-δ) < T u (x)} implies that there exists some time r ∈ [e (1-2ε)k -s, e k(2+2 √ 2-δ) -s]
such that the rightmost position R(k + s + r) is less than -akβs + m(r) + x. Thus,

Ω b2 = P ∃u ∈ Z a (k) s.t. N u (s) = 1, L u (s) ≤ -βs, e (1-2ε)k ≤ τ u ≤ e k(2+2 √ 2-δ) < T u (x) ≤ P ∃r ∈ [e (1-2ε)k , e k(2+2 √ 2-δ) + k], s.t. R(r) ≤ m(r -k -s) + x -ak -βs . By taking a = β = 2 - √ 2, we obtain s = 1-a 2 2 -2ε 1+ β 2 2 +ε k = ( 1 √ 2 -ε 1 )k for some sufficiently small ε 1 = ε 1 (ε) > 0. Let δ ≥ 8 √ 2ε 1 , then m(r -k -s) + x-ak -βs ≤ m(r) -ε 1 k. By Proposition 2.5, for k large enough, (5.12) Ω b2 ≤ C 1 e -ε 1 k/6 √ 2 .
Since Ω b ≤ Ω b1 + Ω b2 , it follows from (5.11) and (5.12) that

Ω b ≤ c 11 k 2 e -εk + C 1 e -ε 1 k/6 √ 2 .
Combined with the fact that Ω a ≤ 2e -ηk , (5.4) implies that

Ω 1 ≤ Ω a + Ω b ≤ 2e -ηk + c 11 k 2 e -εk + C 1 e -ε 1 k/6 √ 2 .
According to the inequality (5.1), for 0

< ε < δ/8, η(ε) > 0, 0 < ε 1 (ε) ≤ δ/8 √ 2 with δ
sufficiently small, and for all k sufficiently large,

P[D k ∩ E k ] ≤ exp -e k(1-a 2 2 -ε)-s(1+ β 2 2 +ε) + Ω 1 ≤ e -e εk + 2e -ηk + c 11 k 2 e -εk + C 1 e -ε 1 k/6 √ 2 .
Consequently, [START_REF] Biggins | Martingale convergence in the branching random walk[END_REF] showed that almost surely (5.13) lim

k P[D k ∩ E k ] < ∞. By Borel-Cantelli lemma, P[D k ∩ E k i.o.] = 0. Recall that E k := Z a (k) ≥ exp k(1 -a 2 2 - ε) . Biggins
k→∞ log Z a (k) k = 1 - a 2 2 .
Therefore, for any δ > 0 small, lim inf k→∞

log Θ k k ≥ 2 + 2 √ 2 -δ almost surely.
This implies the lower bound of Theorem 1.2.

The upper bound of Theorem 1.2

It remains to prove the upper bound, namely, almost surely, (5.14) lim sup

k→∞ log Θ k k ≤ 2 + 2 √ 2.
Before bringing out the proof of (5.14), let us state some preliminary results first.

For M ∈ N + , define

σ M := inf{s > 0; N(s) = M + 1}.
Clearly, σ M is a stopping time with respect to {F s ; s ≥ 0}. Since N(s) follows the geometric distribution, one sees that for any s ≥ 0,

P[σ M ≤ s] = P[N(s) ≥ 1 + M] = (1 -e -s ) M .
Moreover, σ M has a density function, denoted by f M , as follows:

(5.15)

f M (s) := 1 (s≥0) Me -s (1 -e -s ) M -1 ≤ 1 (s≥0) Me -s .
Recall that L(s) = inf u∈N (s) X u (s). Let L(σ M ) denote the leftmost position at time σ M .

Notice that at time σ M , there are M +1 particles which occupy at most M different positions. This tells us that for any s, µ > 0, (5.16)

P L(s) ≤ -µs σ M = s ≤ MP B s ≤ -µs ≤ M µs e -µ 2 s/2 ,
where the last inequality holds because of (2.4). Let ε ∈ (0, 1/2). For r > 1/ε and 0 < s < r, we set λ := λ(s, r) > 0 such that

s(1 + λ 2 2 ) = r. Let Φ(r, λ) := {σ M > r -1} ∪ {εr ≤ σ M ≤ r -1, L(σ M ) ≤ -λ(σ M , r)σ M }, (5.17) Ψ(r, λ) := εr ≤ σ M ≤ r -1, L(σ M ) ≥ -λ(σ M , r)σ M . (5.18)
We have the following lemma, which gives some results of the random vector (σ M , L(σ M )).

Lemma 5.1 (i)

There exists a constant c 12 > 0 such that

(5.19) P Φ r, λ ≤ c 12 M 2 re -r .
(ii) There exists a constant c 13 > 0 such that

(5.20) E e 2σ M - √ 2L(σ M ) ; Ψ(r, λ) ≤ c 13 M 2 r 2 e √ 2r .
Proof of Lemma 5.1. (i) Observe that

P Φ r, λ ≤ P σ M > r -1 + P εr ≤ σ M ≤ r -1, L(σ M ) ≤ -λ(σ M , r)σ M = ∞ r-1 f M (s)ds + r-1 εr P L(s) ≤ -λ(s, r)s σ M = s f M (s)ds ≤ ∞ r-1 Me -s ds + r-1 εr M λ(s, r)s e -λ(s,r) 2 s/2 Me -s ds,
where the last inequality follows from (5.15) and (5.16). A few lines of simple computation yield (5.19).

(ii) Let us prove the inequality (5.20). By Fubini's theorem, we rewrite the expectation

E e 2σ M - √ 2L(σ M ) ; Ψ(r, λ) as follows: r-1 εr e 2s f M (s)E e - √ 2L(s) ; L(s) ≥ -λ(s, r)s σ M = s ds (5.21) = r-1 εr e 2s f M (s) +∞ -λ(s,r)s √ 2e - √ 2x P -λ(s, r)s ≤ L(s) ≤ x σ M = s dxds.
For εr ≤ s ≤ r, one sees that λ(s, r) = 2(rs)/s ≥ 2/(r -1) > 0. We choose

0 < λ 0 = min{1 - √ 2/2, 2/(r -1)} so that +∞ -λ(s,r)s √ 2e - √ 2x P -λ(s, r)s ≤ L(s) ≤ x σ M = s dx ≤ +∞ -λ 0 s √ 2e - √ 2x dx + -λ 0 s -λ(s,r)s √ 2e - √ 2x P -λ(s, r)s ≤ L(s) ≤ x σ M = s dx ≤ e ( √ 2-1)s + -λ 0 s -λ(s,r)s √ 2e - √ 2x P -λ(s, r)s ≤ L(s) ≤ x σ M = s dx. (5.22)
The last term on the right-hand side of (5.22), by a change of variable x = -µs, becomes

λ(s,r) λ 0 √ 2se √ 2µs P -λ(s, r)s ≤ L(s) ≤ -µs σ M = s dµ ≤ λ(s,r) λ 0 √ 2se √ 2µs P L(s) ≤ -µs σ M = s dµ ≤ λ(s,r) λ 0 √ 2se √ 2µs M µs e -µ 2 s/2 dµ,
where the last inequality comes from (5.16). Going back to (5.21),

E e 2σ M - √ 2L(σ M ) ; Ψ(r, λ) ≤ r-1 εr e 2s f M (s) λ(s,r) λ 0 √ 2s M µ e ( √ 2µ-µ 2
2 )s dµ + e ( √ 2-1)s ds.

By (5.15), this is bounded by

r-1 εr λ(s,r) λ 0 √ 2s M 2 µ exp (1 + √ 2µ - µ 2 2 )s dµds + M r-1 εr e √ 2s ds. Notice that if s ≤ r/2, then λ(s, r) = 2( r s -1) ≥ √ 2. It follows that (5.23) max 0<µ≤λ(s,r) (1 + √ 2µ - µ 2 2 )s = 2s ≤ r < √ 2r. Otherwise, r/2 < s < r implies λ(s, r) < √ 2. Hence, max 0<µ≤λ(s,r) (1+ √ 2µ-µ 2
2 )s is achieved when µ = λ(s, r), which equals

(5.24) 1 + √ 2λ(s, r) - λ(s, r) 2 2 s = 1 + √ 2λ(s, r) -λ(s,r) 2 2 1 + λ(s,r) 2 2 r.
It is bounded by

√ 2r since max z≥0 1+ √ 2z-z 2 /2 1+z 2 /2 = √ 2.
Combing the two cases, we obtain that max ǫr≤s≤r-1;0<µ≤λ(s,r)

e (1+ √ 2µ-µ 2 2 )s ≤ e √ 2r .
This implies that

E e 2σ M - √ 2L(σ M ) ; Ψ(r, λ) ≤ r-1 εr λ(s,r) λ 0 √ 2s M 2 µ e √ 2r dµds + Mre √ 2r ≤ M 2 e √ 2r r-1 εr √ 2s λ(s, r) λ 0 ds + Mre √ 2r .
As

λ 0 = min{1- √ 2/2, 2/(r -1)}, we deduce that E e 2σ M - √ 2L(σ M ) ; Ψ(r, λ) ≤ c 13 M 2 r 2 e √ 2r ,
which completes the proof of (ii) in Lemma 5.1.

Let us turn to prove the upper bound of Θ k .

Proof of (5.14). For any u ∈ N (k) and t > 0, we denote A u (t) := {τ u > t}. Then

{Θ k > t} = ∪ u∈N (k) A u (t).
For any θ ∈ Q ∩ (0, 1), let

a j := √ 2 -jθ, b j := √ 2 -(j -1)θ, for j = 1, . . . , K := K(θ) = ⌊ √ 2 θ ⌋, so that 0 < a j < √ 2 for all j ≤ K. Let I k (a, b) := {u ∈ N (k); ak ≤ X u (k) ≤ bk} for -∞ < a < b < ∞. Given the event Ξ := - √ 2k ≤ L(k) ≤ R(k) ≤ √ 2k , we can write N (k) = I k (-θ, √ 2) ∪ 1≤j≤K I k (-b j , -a j ) , so that {Θ k > t} = u∈N (k) A u (t) = u∈I k (-θ, √ 2) A u (t) ∪ 1≤j≤K u∈I k (-b j ,-a j )
A u (t) .

As a consequence, (5.25)

P Θ k > t ∩ Ξ ≤ P u∈I k (-θ, √ 2) A u (t) ∩ Ξ + K j=1 P u∈I k (-b j ,-a j )
A u (t) ∩ Ξ .

We first estimate

P u∈I k (-θ, √ 2) A u (t) ∩ Ξ . For any particle u ∈ N (k), let σ u M := inf{s > 0; N u (s) = 1 + M}. Recall that L u (t) = min{X v (t + k) -X u (k); v ∈ N (k + s), u < v} for any t > 0. By the branching property, conditioned on F k , {σ u M , L u (σ u M )} u∈N (k) are i.i.d. copies of (σ M , L(σ M )). Similarly, We define Φ u (r, λ) := {σ u M > r -1} ∪ {εr ≤ σ u M ≤ r -1, L u (σ u M ) ≤ -λ(r, σ u M )σ u M } and Ψ u (r, λ) := εr ≤ σ u M ≤ r -1, L u (σ u M ) ≥ -λ(σ u M , r)σ u M .
One immediately observes that (5.26)

P u∈I k (-θ, √ 2) A u (t) ∩ Ξ ≤ P u∈I k (-θ, √ 2)
Φ u (r, λ)

+ P u∈I k (-θ, √ 2) A u (t) ∩ Φ u (r, λ) c ; Ξ .
Conditioning on F k yields that (5.27)

P u∈I k (-θ, √ 2) Φ u (r, λ) ≤ E u∈I k (-θ, √ 2) 1 Φ u (r,λ) = E u∈I k (-θ, √ 2
)
1 P Φ(r, λ) .

Clearly,

(5.28)

E u∈I k (-θ, √ 2) 1 ≤ E[N(k)] = e k .
It then follows from (5.19) that (5.29)

P ∪ u∈I k (-θ, √ 2) Φ u (r, λ) ≤ c 12 M 2 re -r+k .
We now choose r = k(1 + ε) and set Λ 0 := P ∪ u∈I k (-θ, √

2) A u (t) ∩ Φ u (r, λ) c ; Ξ . Then for all k large enough, (5.26) becomes (5.30)

P u∈I k (-θ, √ 2)
A u (t) ∩ Ξ ≤ c 14 M 2 ke -εk + Λ 0 .

It remains to estimate Λ 0 . Since Φ u (r, λ)

c ⊂ {σ u M < εr} ∪ Ψ u (r, λ), we write (5.31) Λ 0 = P u∈I k (-θ, √ 2) A u (t) ∩ Φ u (r, λ) c ; Ξ ≤ Λ 1 + Λ 2 ,
where

Λ 1 := P u∈I k (-θ, √ 2) A u (t) ∩ {σ u M < εr}; Ξ Λ 2 := P u∈I k (-θ, √ 2)
A u (t) ∩ Ψ u (r, λ); Ξ . T v (y) + εr.

Recall that A u (t) = {τ u > t}. By comparing S u (y) with t, we obtain that for any t 1 ∈ (0, t),

Λ 1 ≤ P u∈I k (-θ, √ 2)
t 1 ≤ S u (y) ≤ t < τ u ; Ξ (5.33)

+ P u∈I k (-θ, √ 2) S u (y) > t; σ u M < εr + P u∈I k (-θ, √ 2)
S u (y) < t 1

=: Λ 1a + Λ 1b + Λ 1c .
For δ ∈ (0, 1) and θ ∈ Q ∩ (0, 1), we take y = ( 1+2δ) .

√ 2 + 2)k(1 + θ), t 1 = e k , t = e √ 2y (
As {σ u M < εr} implies {N u (εr) > M}, P[ S u (y) > t; σ u M < εr|F k+εr ] is less than P[T (y) > tεr] M . Conditionally on F k+εr , then by (5.28), we get , since r = (1 + ε)k with ǫ ∈ (0, 1/2). By (2.16), Λ 1b ≤ e k × e -M δ 1 y/3 for k large enough. We take M = 6 δ 1 to ensure that Λ 1b ≤ e -k . On the other hand, we observe that

Set Ξ 1 := R(k + εr) ≤ 2(k + εr) ∩ Ξ. Then, (5.36) Λ 1a ≤ P R(k + εr) > 2(k + εr) + E u∈I k (-θ, √ 2)
1 (τu>t≥ S u (y)≥t 1 ) ; Ξ 1 .

We define Λ 1rest := E u∈I k (-θ, √

2) 1 (τu>t≥ Su (y)≥t 1 ) ; Ξ 1 for convenience. On the one hand, P[R(k + εr) > 2(k + εr)] ≤ e -k because of (2.15). On the other hand, we have

Λ 1rest ≤ E u∈I k (-θ, √ 2) t t 1 1 ( S u (y)∈dr ′ ) P τ u > r ′ ; Ξ 1 F k , F u ∞ .
Since { S u (y) = r ′ } ⊂ {R u (r ′ ) ≥ L u (εr) + m(r ′εr) + y}, the event {τ u > r ′ } conditioned on { S u (y) = r ′ } implies ∪ w∈N (k)\{u} X w (k)+R w (r ′ ) > X u (k)+L u (εr)+m(r ′ -εr)+y . Further,

this set is contained in ∪ w∈N (k)\{u} R w (r ′ ) > m(r ′ ) + X u (k) + L u (εr) + y - √ 2εr -X w (k) . As Ξ 1 guarantees that X u (k) + L u (εr) + y - √ 2εr -X w (k) + ≤ C 2 k, the inequality (2.2)
can be applied to show that εr) ,

E u∈I k (-θ, √ 2) t t 1 1 ( S u (y)∈dr ′ ) P τ u > r ′ ; Ξ 1 F k , F u ∞ ≤ E u∈I k (-θ, √ 2) t t 1 1 ( Su (y)∈dr ′ ) w∈N (k)\{u} c 2 (1 + C 2 k) 2 e - √ 2 Xu(k)+L u (εr)+y- √ 2εr-Xw(k) ≤ c 15 k 2 e - √ 2y+2εr E u∈I k (-θ, √ 2) w∈N (k)\{u} e √ 2(Xw(k)-Xu(k)) E e - √ 2L (
where the last inequality follows from the fact that for u ∈ N (k), L u (εr) are independent of F k and are independent copies of L(εr).

Whereas by the estimation of

E u∈N δ (k) w∈N (k)\{u} e √ 2Xw(k) in Section 4, E u∈I k (-θ, √ 2) w∈N (k)\{u} e √ 2(Xw(k)-Xu(k)) = 2 k 0 e 2k-s ds R P B s ∈ dx E e - √ 2B k ; -θk ≤ B k ≤ √ 2k B s = x E e √ 2B k B s = x = 2 k 0 e 3k-2s E e √ 2Bs- √ 2B k ; -θk ≤ B k ≤ √ 2k ds. Because E e √ 2Bs- √ 2B k ; -θk ≤ B k ≤ √ 2k ≤ e √ 2θk E e √ 2Bs = e √ 2θk+s
, we obtain that (5.37)

E u∈I k (-θ, √ 2) w∈N (k)\{u} e √ 2(Xw(k)-Xu(k)) ≤ 2e 3k+ √ 2θk .
Besides, E e - √ 2L(εr) ≤ E v∈N (εr) e - √ 2Xv(εr) = e 2εr . As a result,

Λ 1rest ≤ c 15 k 2 e - √ 2y+2εr × 2e 3k+ √ 2θk+2εr . Recall that 0 < θ < 1, y = ( √ 2 + 2)k(1 + and that r = k(1 + ε) with k large enough so that r > 1/ε. Hence, Λ 1rest ≤ c 15 k 2 e (1-2 √ 2)k for ε ∈ (0, θ 3

It remains to estimate Λ

2 = P u∈I k (-θ, √ 2) A u (t) ∩ Ψ u (r, λ); Ξ where Ψ u (r, λ) = εr ≤ σ u M ≤ r -1, L u (σ u M ) ≥ -λ(σ u M , r)σ u M . For any particle u ∈ N (k) satisfying Ψ u (r, λ), define (5.39) S u (y) := min v∈N u (σ u M )
T v (y) + σ u M , for any y > 0.

Comparing S u (y) with t yields that

(5.40) Λ 2 ≤ P u∈I k (-θ, √ 2) t 1 ≤ S u (y) ≤ t < τ u ∩ Ψ u (r, λ); Ξ + P u∈I k (-θ, √ 2) S u (y) > t ∩ Ψ u (r, λ) + P u∈I k (-θ, √ 2) S u (y) < t 1 ∩ Ψ u (r, λ) =: Λ 2a + Λ 2b + Λ 2c .
According to the definition of S u (y), one sees that P[ S u (y) > t; Ψ u (r, λ)] ≤ P[T (y) > tr] M and that P[ S u (y) < t 1+2δ) . For any -∞ < a < b < ∞, by (2.16),

1 ; Ψ u (r, λ)] ≤ 2MP[T (y) < t 1 ]. Recall that r = (1 + ε)k, M = 6/δ 1 , y = ( √ 2 + 2)k(1 + θ), t 1 = e k and t = e √ 2y (
(5.41)

E u∈I k (a,b) 1 { S u (y)>t}∩Ψ u (r,λ) ≤ e k × P[T (y) > e √ 2y(1+δ) M ≤ e -k .
Meanwhile, by (2.13),

(5.42)

E u∈I k (a,b) 1 { S u (y)<t 1 }∩Ψ u (r,λ) ≤ e k × 2MP[T (y) < t 1 ] ≤ 2c 5 Me -2k .
Hence, taking a = -θ and b

= √ 2 implies that Λ 2b + Λ 2c ≤ e -k + 2c 5 Me -2k . Let Ξ 2 := {max 0≤r 0 ≤r R(k + r 0 ) ≤ 6k} ∩ Ξ. We get (5.43) Λ 2a ≤ P max 0≤r 0 ≤r R(k + r 0 ) > 6k + E u∈I k (-θ, √ 2) 1 {t 1 ≤ S u (y)≤t<τu}∩Ψ u (r,λ) ; Ξ 2 .
By the many-to-one lemma, for k large enough, (5.44) P max

0≤r 0 ≤r R(k + r 0 ) > 6k ≤ e -k .
For the second term on the right-hand side of (5.43), we need to recount the arguments to estimate Λ 1rest . Let (5.45)

Λ 2rest := E u∈I k (-θ, √ 2) 1 {t 1 ≤ S u (y)≤t<τu}∩Ψ u (r,λ) ; Ξ 2 .
It immediately follows that

(5.46) Λ 2rest ≤ E   u∈I k (-θ, √ 2) t t 1 1 ( S u (y)∈dr ′ ) 1 Ψ u (r,λ) × P {τ u > r ′ } ∩ Ξ 2 F k , F u ∞   .
Comparing τ u with S u (y) tells that

Λ 2rest ≤ E u∈I k (-θ, √ 2) t t 1 1 S u (y)∈dr ′ 1 Ψ u (r,λ) × E 1 Ξ 2 w∈N (k)\{u} 1 R w (r ′ )>m(r ′ )- √ 2σ u M +y+Xu(k)-Xw(k)+L u (σ u M ) F k , F u ∞ . On the event Ξ 2 ∩ Ψ u (r, λ), we have 1 + - √ 2σ u M + y + X u (k) -X w (k) + L u (σ u M ) + ≤ C 3 k.
Applying the inequality (2.2) for R w (r ′ ) yields that In view of (5.41), (5.42) and (5.44),

Λ 2rest ≤ E u∈I k (-θ, √ 2
(5.53) Λ ′ 2 ≤ e -k + 2c 5 Me -2k + e -k + E u∈I k (-b j ,-a j )

1 {t 1 ≤ S u (y)≤t<τu}∩Ψ u (r,λ) ; Ξ 2 .

We define Λ ′ 2rest := E u∈I k (-b j ,-a j ) 1 {t 1 ≤ S u (y)≤t<τu}∩Ψ u (r,λ) ; Ξ 2 . Thus applying the analogous arguments to the estimation of Λ 2rest gives that This implies the upper bound in Theorem 1.2.

2

 2 It suffices to show that almost surely lim N∋k→∞ log Θ k k = 2 + 2 √ 2, as the sequence {Θ s ; s > 0} is monotone.

2

  )k. Then comparing T u (z) and e ρk for every u ∈ Z a (k) tells us that (5.5) Ω a ≤ P u∈Za(k) {T u (z) < e ρk } + P u∈Za(k) {τ u < e ρk ≤ T u (z)} . It follows from the branching property that the first term of the right-hand side is bounded by E[Z a (k)]P[T (z) < e ρk ], which is e k P[B k ≤ -ak]P[T (z) < e ρk ] by the many-to-one lemma.

For

  any particle u ∈ I k (-θ, √ 2) such that {σ u M < εr}, for any y ≥ 1, we define (5.32) S u (y) := min v∈N u (εr)

( 5 . 2 ) 1 P

 521 34) Λ 1b ≤ E u∈I k (-θ, √ T (y) > tεr M ≤ e k P T (y) > e √ 2y(1+δ) M

1 T

 1 v (y)<e k ≤ e k+εr P T (y) < e k . By (2.13), Λ 1c ≤ e k+εr c 5 y 2 e - √ 2y e k . Thus Λ 1b + Λ 1c ≤ 2e -k for sufficiently large k.

2 √ 2 ) 2 √ 2 )

 2222 ). Going back to (5.36), we get Λ 1a ≤ e -k + c 15 k 2 e (1-k for all k sufficiently large. Consequently, (5.33) becomes (5.38) Λ 1 ≤ e -k + c 15 k 2 e (1-k + 2e -k ≤ c 16 e -εk .

) 1 Ψ 2 (

 12 u (r,λ) × Xw(k)-Xu(k)) E e 2σ M - √ 2L(σ M ) ; Ψ(r, λ) ,by the fact that {σ u M , L u (σ u M )} are i.i.d. and independent of F k . Recall that r = (1 + ε)k with ε ∈ (0, θ/3). It then follows from (5.37) and (5.20) that(5.47) Λ 2rest ≤ C 2 3 k 2 e -(2+2 √ 2)k(1+θ) × 2e 3k+ √ 2θk × c 13 M 2 r 2 e √ 2r ≤ c 17 k 4 M 2 e -( √ 2-1)k .

( 5 .√ 2 (√ 2 ( 2 k 0 e√ 2 (

 522202 54) Λ ′ 2rest ≤ c 21 k 2 e - √ 2y E u∈I k (-b j ,-a j ) w∈N (k)\{u} e Xw(k)-Xu(k)) × E e 2σ M - √ 2L(σ M ) ; Ψ u (r, λ) .Once again, by means of integrating with respect to the last time at which the most recent common ancestor of u and v was alive,E u∈I k (-b j ,-a j ) w∈N (k)\{u} e Xw(k)-Xu(k)) equals 2k-s ds R P B s ∈ dx E e - √ 2B k ; -b j k ≤ B k ≤ -a j k B s = x E e √ 2B k B s = x 2B k ; -b j k ≤ B k ≤ -a j k ds k ∈ dx E e √ 2Bs B k = x .Let (b s (x); 0 ≤ s ≤ k) denote a Brownian bridge from 0 to x of length k. Then E e √ 2Bs B k = x equals E e √ 2bs(x) , which turns out to be exp(s(ks + √ 2x)/k). Note that a j > 0 and that b j = a j + θ. This gives that E u∈I k (-b j ,-a j ) w∈N (k)\{u} e Xw(k)-Xu(k))

2 + 2 θP 2 ≤ C 1

 2221 c 20 e -εk + c 23 M 2 e -θk .Summing over j ∈ {1, • • • , K = ⌊ √ ⌋} ∪ u∈I k (-b j ,-a j ) A u (t); Ξ ≤ C(θ)M 2 e -εθk/2 ,where C(θ) is a positive constant associated with θ (but independent of k, δ and M) and k is large enough.Going back to (5.25), we combine (5.49) and (5.57) to say thatP τ (k) > t ∩ -√ 2k ≤ L(k) ≤ R(k) ≤ √ 2k ≤ c 19 M 2 e -εk/2 + C(θ)M 2 e -εθk/(θ) δ 2 e -εθk/2 , where θ ∈ Q ∩ (0, 1), ε ∈ (0, θ 3 ), δ > 0 and t = exp[k(2 + 2 √ 2)(1 + θ)(1 + 2δ)] and k is sufficiently large.According to the Borel-Cantelli Lemma, we conclude that for any θ ∈ Q ∩ (0, 1) and any θ)(1 + 2δ), almost surely.

  choose a y = δ

	4+2 √ 2y(1+ 1 2 y =: δ 1 y from now on to evaluate P[Σ]. Since Σ 1 ∈ F ay , for y large 2 δ) ≤ e √ 2y(1+δ) -a y , we have P Σ F ay ≤ 1 Σ 1 enough so that 2e √ u∈N (ay) P R u (r) ≤ m(a y + r) + y -X u (a y ), ∀r ≤ e √ 2y(1+δ) -a y F ay
	≤ P R(r) ≤ m(r) + y + 2a y +	√	2a y , ∀r ∈ e √	2y(1+ 1 2 δ) , 2e √	2 δ) 2y(1+ 1	e ay /2
	≤ P R(r) ≤ m(r) +	1 √ 2	log r, ∀r ∈ e √	2y(1+ 1 2 δ) , 2e √	2 δ) 2y(1+ 1	e ay /2

  ′ 2rest ≤ c 22 k 9/2 M 2 e -2θk, and thus for all k sufficiently large,(5.56) Λ ′ 2 ≤ c 23 M 2 e -θk .Consequently, by (5.52),P ∪ u∈I k (-b j ,-a j ) A u (t); Ξ ≤ c 12

	by	√	ke 2k+(1-	a 2 j 2 +	√	2a j )k+	√	2θk . Going back to (5.54) and applying (5.20),
			Λ ′ 2rest ≤ c 21 k 2 e -(2+2 √ ≤ c 22 k 9/2 M 2 e -2θk exp k ( 2)k(1+θ) √ ke 2k+(1-√ 2 + 1)(1 -a 2 a 2 j 2 + √ 2a j )k+ √ 2θk c 13 M 2 r 2 e √ j /2) + √ 2a j -2 2r √	2 ,
	as r = k(1 -(5.55)	a 2 j 2 )(1 + ε) with ε ∈ (0, θ 3 ). Observe that ( √ 2 + 1)(1 -a 2 j /2) + √ 2a j -2 √ 2 = -	√	2 + 1 2	a j -(2 -	√	2)	2	≤ 0.
	We get Λ
												√	1 2πk	exp -	x 2 2k	-	√	2x +	√	2	s k	x dxds
			≤ 2	0	k	e 3k-2s+s(k-s)/k θk √ 2πk	exp -	a 2 j 2	k +	√	2b j k -	√	2a j s ds,
	which is bounded by is less than 1. One hence sees that E √ k exp 3k -a 2 j 2 k +	√	2b j k since	k 0 exp -2s + s(k -s)/k -√ 2(Xw(k)-Xu(k)) is bounded √ 2a j s ds

u∈I k (-b j ,-a j ) w∈N (k)\{u} e

Acknowledgements

I would like to thank Zhan Shi for advice and help.

Combined with (5.38), Λ 0 ≤ c 16 e -εk + c 18 M 2 e -εk . Going back to (5.30), we conclude that (5.49)

To complete the proof, we still need to evaluate P ∪ u∈I k (-b j ,-a j ) A u (t) ∩ Ξ . Recall that for any particle u ∈ N (k), σ u M = inf{s > 0; N u (s) = 1 + M} and

for any r > 1/ε and λ(s, r) = 2( r s -1) with 0 < s < r. Clearly,

On the one hand, (5.50)

1 P Φ(r, λ) .

29

We now take r = k(1 -

θ ⌋, with θ ∈ Q ∩ (0, 1). Then note that each a j is strictly positive. Thus, by the many-to-one lemma and (5.19), (5.50) becomes that

where the last inequality follows from (2.4). P u∈I k (-b j ,-a j ) Φ u (r, λ) is hence bounded by

On the other hand, recalling that Ψ u (r, λ)

where

A u (t) ∩ Ψ u (r, λ); Ξ .

Furthermore, by an argument similar to the one used in estimating Λ 1 , we have Λ ′ 1 ≤ c 20 e -εk . Thus, (5.52)

{ S u (y) > t} ∩ Ψ u (r, λ) + P u∈I k (-b j ,-a j )

{ S u (y) < t 1 } ∩ Ψ u (r, λ)

+ P max 0≤r 0 ≤r R(k + r 0 ) > 6k + P u∈I k (-b j ,-a j ) t 1 ≤ S u (y) ≤ t < τ u ∩ Ψ u (r, λ); Ξ 2 .


Pierre-Olivier Malaterre
email: pierre-olivier.malaterre@cemagref.fr

David Dorchies

Jean-Pierre Baume

New method for coupling 1D unsteady flow hydraulic models using RST equations

Keywords: Coupling, Numerical modeling, Hydraulic, Preissmann, Saint-Venant, 1D, RST, double-sweep, PALM

This paper is devoted to the presentation and evaluation of different alternatives to couple two 1D hydraulic models based on the Saint-Venant's equations. This presentation is first giving a short description of the coupling platform we use for these tests, namely the PALM coupling platform developed by Cerfacs. Then, we present 2 classical methods for this coupling: 1-way and 2-way. The 2-way method has itself different options depending on the variables exchanged at the coupling interface. Since we are using the implicit Preissmann scheme to solve numerically the Saint-Venant's equations, we can imagine another option using the so-called RST equations. This method is not very intrusive into the source code since it does not require any additional calculation, and can improve a lot the results obtained by the 2 coupled models compared to the global one. This method can be implemented in an implicit way, being equivalent to the global model, or explicitly with some degradation of the results.

Introduction

The coupling of different types of models is gaining more and more attention recently. This is due, in particular, to the needs of more global model encompassing different disciplines (eg.: multi-physics) and different approaches (eg.: multi-scale, nesting, 1D-2D-3D) (see [?], [?], [?], [?], [?], [?], [?]). Also, the possibility to assemble different modeling units inside a friendly modeling software platform is an attractive solution compared to developing more and more global complex models. This coupling problem can be seen from at least 2 very different perspectives: coding or software aspects (how to put the models together, link them, run simulations and get results) and physical or algorithmic aspects (what variables must be exchanged at the coupling interfaces and how).

Software aspects of coupling

On the coding or software aspects of coupling there have been lots of recent developments with different alternatives and available packages. We can, for example, quote the OpenMi European initiative defining a standard communication protocol between different softwares implementing this interface. This approach implies the development of specific routines (OpenMi API) into the model code(s). The different models can then exchange data at their interfaces without using other third party platforms. But this approach has also several limitations, for example not allowing to couple several instances of the same program. Also, since the two models are communicating directly between them, without being included into a supervisory platform, the coupled simulations cannot benefit from high level toolboxes such as optimization, sensitivity analysis or data assimilation ones.

We can also quote other approaches such as the PALM coupling platform developed by Cerfacs [?]. The original model codes have to be "Palmerized", which means slight local modifications to declare the variables that can be exchanged (input or/and output) with the other models. Then, a user friendly interface named PrePALM (Fig. If the flow conditions are supercritical (Froude number F > 1) the calculation can be done for each reach separately, in cascade from upstream to downstream. Indeed, two upstream boundary conditions are used for the calculation of either reach and no downstream boundary condition is required. In this specific case, the calculation in cascade of reach 1 and then reach 2 using a 1-way coupling (providing 2 variables at the coupling interface) would give the same results as the global calculation of the entire 2-reach river. The same remark is valid for the coupling of hydrologic and hydraulic models when there is no feedback from the hydraulic model towards the hydrologic one (see [?] [?]).

But in most cases, irrigation canals or rivers are under subcritical flow conditions, at least in most of their spatial sections, and during most of the time. In this case, we know that the hydraulic calculation of a reach needs one upstream and one downstream boundary condition. This is due to the hyperbolic nature of the equations and with two opposite directions for the characteristic lines. It is then impossible to make the calculation in cascade (in either order) of one reach after the other without simplification hypothesis and therefore degradation of the original global solution.

1-way coupling

We can first check what happens if we handle the coupling interface only 1-way. In this case, we propagate the flow at the coupling interface from upstream to downstream. For the downstream boundary condition of the first upstream reach we keep a local rating curve Q(Z). We observe in this case, that the approximate solution obtained with this 1-way coupling has errors compared to the global true solution, either in terms of discharges Q and water elevations Z (Fig. ??). This error does exist even if we reduce the calculation time step ∆t since this is a physical modification of the system, and not only a numerical one. This error is plotted in Fig. ?? at the coupling interface, but also propagates further upstream and downstream.

2-way coupling

In this case we exchange two variables at the coupling interface, so as to provide one required boundary condition for each reach. The information used for one reach is coming from the other one. Eventhough we know the number and locations of the required boundary conditions to get a well-posed problem, there exist many different options for these variables. The basic option presented here (Fig. ??) consists in providing the flow Q from reach 1 to 2 and the water elevation Z from reach 2 to 1. The physics is respected, but the exchange of coupling variables being made with one time step delay, we can observe some oscillatory behavior.

The above results can be improved by reducing the calculation time steps ∆t 1 and ∆t 2 of the two coupled models (Fig. ??). But in this case we loose the advantages of using an implicit scheme such as the Preissmann scheme.

In order to overcome the problems illustrated in this section, we have again several options. The first one is to make iterations of the coupled models, until the variables exchanged at the interface converge towards the same values. This needs a supervision platform, such has PALM allowing the management of this type of iterations. The advantage of this method is to be non-intrusive, since this does not necessitate any further modification of the several codes involved. But there are several drawbacks. As for any iterative methods, the questions of convergence should be addressed. The required CPU time can be large, eventhough the number of iterations can be minimized for example using methods based on the Riemann invariant variables (cf. [?]).

New solution for coupling

In order to present a new method for coupling 1D models, we give some more details on the equations and numerical scheme we are using. The dynamic behavior of water in an open-channel is well described by the so-called Saint-Venant's equations:

∂Q ∂x + ∂S ∂t = 0 ∂Q ∂t + ∂ Q 2 S ∂t + g.S ∂Z ∂t + g.S.J = 0 ( 1
)
where Q is the discharge (m 3 /s), S the wetted cross-area (m 2 ), Z the water elevation (m), J the friction slope, x the longitudinal abscissa (m) and t the time (s). The friction slope J is usually obtained from the Manning-Strickler formula:

J = n 2 Q 2 S 2 R 4 3
, where n is the Manning's coefficient (0.02) and R is the hydraulic radius (m) (R = S P , where P is the wetted perimeter). The two hyperbolic, first-order, non-linear, partial-derivative equations Eq. (??) are discretized and linearized in time (∆t time step) and space (∆x space step) through the implicit Preissmann finite difference scheme ([?]). The obtained set of numerical equations are solved using the double-sweep method.

The two Saint-Venant's equations after the discretization between two sections j and j + 1 give the following relations ([?]):

∆Z j = e∆Z j+1 + d∆Q j+1 + f (2)
∆Q j = b∆Z j+1 + a∆Q j+1 + c (3)
The upstream boundary condition is linearized:

r 1 ∆Q 1 + s 1 ∆Z 1 = t 1 (4)
where r 1 , s 1 , t 1 are three known coefficients corresponding to the first section of the reach (subscript 1). This equation is named impedance relation for the section 1. The first sweep consists in triangularizing the matrix of 2n equations (n sections within the reach, so 2(n -1) Saint-Venant equations and 2 boundary conditions). That means to compute the impedance relation in each section by recurrence:

If r j , s j , t j are known values at section j, the impedance relation for section j can be written as:

r j ∆Q j + s j ∆Z j = t j (5)
Replacing equation (??) by r j (??) -(??) + s j (??) gives:

r j+1 ∆Q j+1 + s j+1 ∆Z j+1 = t j+1 (6)
with

   r j+1 = r j a + s j d s j+1 = r j b + s j e t j+1 = t j -r j c + s j f (7)
The impedance relation in the last section n is then:

r n ∆Q n + s n ∆Z n = t n (8)
The downstream condition is linearized, the same way as for the upstream condition:

r ′ n ∆Q n + s ′ n ∆Z n = t ′ n (9)
where r ′ n , s ′ n , t ′ n are known values, depending on the type of downstream boundary condition.

The solution of equations ?? and ?? gives ∆Z n and ∆Q n . The second sweep allows calculating ∆Z j and ∆Q j for the (n -1) remaining sections j going from downstream to upstream.

Let's suppose we want to couple 2 unsteady flow 1D models at some interface section k. The idea of the RST coupling method we propose in this paper is to exchange the r k , s k , t k variables of the impedance relation for section k. This can be done explicitly with a one time step delay. But this can also be done keeping the implicit feature of the double sweep method. In this case the procedure is the following:

• run the first sweep on reach 1

• provide the 3 r k , s k , t k scalar variables of the impedance relation though the coupling interface, from reach 1 to reach 2

• run the first sweep on reach 2

• run the second sweep on reach 2

• provide the 3 r k , s k , t k scalar variables of the impedance relation though the coupling interface from reach 2 to reach 1

• run the second sweep on reach 1

This method does not require any additional calculation, since the impedance equations are already computed. It just needs exchange of 3 scalar variables through the coupling interface, upward and downward. It is equivalent to the global calculation, but allowing modularity of different instances of 1D hydraulic model to be coupled.

Conclusion

A new method for coupling different instances of 1D hydraulic models has been presented, based on RST equations. This is equivalent to the global calculation in simple cases (2 reaches in cascade). But at the same time, the modularity offered by handling such coupling method allows many interesting applications, using a coupling platform such as PALM. For example, coupling 1D models with other types of models (handling supercritical flows, open channel 1D and sewer 1D models, etc). We can also imagine extensions of this method for the coupling of 1D-2D models since the RST method is providing a boundary condition at the interface.

  ??) allows to construct the global model by linking the different units together. PALM is the platform we use for the tests presented in this paper. In Fig. ?? we have coupled (2-way coupling) two instances of the 1D hydraulic model SIC (or more precisely its unsteady flow module Sirene) developed by Cemagref ([?]).

Figure 1 :

 1 Figure 1: Coupling of models in the Palm platform

Figure 2 :

 2 Figure 2: Discharge and Elevation at the coupling interface, 1-way coupling, ∆t = 600s

Figure 3 :

 3 Figure 3: Discharge and Elevation at the coupling interface, 2-way coupling, ∆t = 600s

Figure 4 :

 4 Figure 4: Discharge at the coupling interface, 2-way coupling, ∆t = 120s

Acknowledgements

This work is supported by the Montpellier Agglomération, Agence de l'Eau RMC and Région de Montpellier within the context of the "Lez multi-purpose management project". This project involves the contributions of BRGM, UMR G-eau, UMR TETIS, UMR HydroScience, Cerfacs and Biotope research institutes.


F Benmakrouha

C Hespel
email: hespel@insa-rennes.fr

E Monnier
email: monnier@insa-rennes.fr

D Quichaud
email: quichaud@insa-rennes.fr

Experiments and design of an inference fuzzy system

Keywords: fuzzy system, diabetes, datum plane covering

), (0). Our method has the advantage of occuring before initialization step and therefore avoiding unfired rules which make difficult to produce an accurate output. .

INTRODUCTION

The lack of data is important and is pointed out by M Lutaud-Brunet in (0). Sugeno and Yasukawa underline in (0) how difficult it is to build a fuzzy model when data are scarce and membership functions don't sweep over all the universe of discourse. This may improve the generalisation ability of fuzzy systems (0). Several methods have been proposed to obtain automatic fuzzy rules from sparse training data. In (0), the authors first construct fuzzy rules from collect data. Then, they use kernel regressions for generate training data. Another technique used when classical inference methods produce sparse fuzzy rules is a diffusion procedure based on interpolation to initialize incomplete rules (0), (0),(0). The objective of this paper is to measure the impact of datum plane covering on the outcome of a fuzzy inference system. Most of optimization methods make the assumption that datum plane is sufficiently covered. If this assumption no longer holds, we will see that these methods cannot work, since it implies that, before optimization, the fuzzy system gives acceptable results. In (0), we analysed the relationship between learning set Ω and labels of fuzzy inference system. In this paper, we take into account a data density repartition, by the measure of number of available data on intervals of each input variable domain.

All these tecniques take place during optimisation step when our method is used before initialization step. This in turn allows one to isolate unfired rules and to proceed, if necessary, to a partial remodelling of the FIS before training and optimization.

The Tagaki-Sugeno Model

The model under consideration in this section is a Tagaki-Sugeno model, which corresponds to discretized linear models of order 1, combined with nonlinear functions.

y(t) = r ∑ i=1 h i (z(t))(a i 1 .y(t -idecal) + b i 1 .u(t -idecal))
y(t) is the output, r is the number of linear models, z(t) a vector which depends lineary or not on the state, h i (z(t)) >= 0 , i = 1, • • • , r nonlinear functions verifying the convex sum property. idecal is the time lag between input and its effect, which is specially interesting in our application.

We have measures about every five minutes and we admit that the effect of insulin (considered in our application)is fast and noticeable ten (idecal=2) minutes later, up to half an hour(idecal = 6). The determination of unknown parameters a i 1 and a i 2 is done by the algorithm of recursive least square.

2 Application to the insulin/glycaemia behavior of diabetics

The available data

The correlated data "'insulin infusion delivery/glycaemia"' has been provided by the team of Pr. Pinget, CHU of Strasbourg. They concern the same person and the same insulin.The insulin infusion has been done by an intra-peritoneal route and the glycaemia has been checked by a subcutaneous sensor. Measures of glycaemia have been made every five minutes during 7 days, which corresponds to 1700 measures.

A bolus is a dose of insulin infused manually, in addition to the basic dose, since postprandial glycemia cannot be regulated satisfactorily. The insulin file contains crude data about basic insulin doses as well as boluses. So, a pretraitement of the insulin file has been necessary to produce a file of insulin delivery for the same person every five minutes.

Experiments and Validation of the model

The learning set is composed of the first measures(280 points) that corresponds to insulin infusion and blood glucose concentration of a patient during a day. We take 7 (r = 7)linear models, considering that each model is valid about three and half hours. The mean square error(MSE)is calculated on the totality of the measures(1700 points).

We make experiments by changing the parameter idecal of the model, time lag between an input and its effect.

The test of our modeling method shows that we can predict the glycaemia over a long period (7 days), by considering glycaemia and insulin delivery 15minute (resp 30-minute) before with an error of about 6%(resp 16%), which is a good result compared with current results. However, we see that results obtained are not so good in the last case, when we consider slow effect insulin (with 2-hours delay). In this case, our model has to be refined, by increasing its order.

Datum Plane Covering

We propose a measure used to pre-validate a fuzzy model. We suppose that there exists a learning set Ω = {(x j , d j )}, where x j is an input vector and d j , the corresponding output. We also assume that the desired function f is defined in

V = [a 1 , b 1 ] × [a 2 , b 2 ] × ... × [a p , b p ]
Usually, to validate a fuzzy inference system, the mean square error (MSE) is calculated on a test set. If the MSE exceeds a threshold, then training is done, using a gradient method. This consists in modifying C j at each presentation of examples from the error (y(x j )d j ). Unfortunately, in case of model invalidation, we cannot determine never learned rules that cause the gap between the model and the real system. Moreover, if there is an insufficient covering of datum plane, training and finer splitting of input space are inefficient and useless. With the criterion proposed below, we estimate the datum plane coverage and we are able to isolate inactivated rules. Then, partial remodeling of the fuzzy inference system is possible. The study is investigating the relationship between a quantitative variable X, number of available data for each input, and a qualitative variable Y , labels of membership functions.

When designing a fuzzy system, we attribute to each input I r modalities (or labels) noted y 1 , • • • y l , • • • y r . We note X I the variable for the input I of average xI and variance σ 2 X I . We note Ω I the corresponding learning set . Each label y l of I defines a subset Ω I l of Ω I : we obtain a partition of Ω I in m classes. We note n I l = card(Ω I l ) and n I = card(Ω I ). We have n I = Σ m l=1 n I l . Then, if we consider the restriction de X I to Ω I l ) (l = 1, • • • , m), we may define the average (noted xI l )and the variance (noted σ 2 X I l of X I on this subset:

xI l = ∑ ω∈Ω I l X(ω) n I l σ 2 X I l = ∑ ω∈Ω I l (X(ω) -xI l ) 2 n I l
We have an index of connection between the datum plane coverage (for an input I) and the learning set defined by :

s I = σ 2 E σ 2 XI where σ 2 X I = σ 2 E + σ 2 R and σ 2 E = ∑ r l=1 n I l * ( xI l -xI ) 2 n I and σ 2 R = ∑ r l=1 n I l * σ 2 X I l n I and xI = ∑ r l=1 (n I l * xI l ) n I
This index of connection consists in detecting relationships between number of data of the learning set Ω I and r I labels. This index is low if the features of these labels are not so different (Test 1). When this index is high, it points to that there is a bad repartition of membership functions (Test 3). This gives an information about the repartition of data of learning set Ω I between membership functions.

Experiments

We have made 3 tests for the first input (with two triangular membership functions) using our application. We give in these array the features of these functions. We obtain for the first test a low index(0.035), the second a medium index (0.28) and for the last test a high index (0.84). In the third test, first membership function is useless and the corresponding rules are inactived. So, we can suppress them without affecting results. We have made a 4th test where three (out of four) membership functions and the associated rules were unnecessary.

graphic representation

We represent sets of data by Box & Whiskers Plots to underline the relation between number of data and labels of membership functions. The first figure (resp figure 2 and figure 3) corresponds to Test1 (resp Test 2 and Test 3).

Conclusion

We have proposed a measure for detecting useless rules and thus pre-validating a fuzzy inference system. When the model is not pre-validated, we have not to carry out next steps, particularly optimization step. We have shown that this criterion gives useful information about datum plane coverage.

  Figure1

  Figure 2

Table 1 :

 1 First Table.

	r idecal MSE
	7	2	0.04
	7	3	0.06
	l
	7	6	0.16
	7	24	1.02

Table 2 :

 2 First membership function.

	Test Center Left corner Right corner
	1	1.4	1.6	1.0
	2	1.0	0.6	1.0
	3	0.2	0.2	0.2
	Test Center Left corner Right corner
	1	2.6	1.0	0.9
	2	2.6	1.0	0.9
	3	3.0	2.8	0.5

Table 3 :

 3 Second membership function.


Irinel-Constantin Morȃrescu
email: constantin.morarescu@univ-lorraine.fr

Wim Michiels
email: wim.michiels@cs.kuleuven.be

Marc Jungers
email: marc.jungers@univ-lorraine.fr

Keywords: Coupled oscillators, distributed delays, synchronization, delayed couplings

We derive conditions for the synchronization of coupled nonlinear oscillators affected by distributed delays in the interconnections. The distributed delays are characterized by a gamma-distribution kernel with a gap. The approach is based on the stability analysis of synchronized equilibria in the (delays,gain) parameter space and the characterization of the structure of the emanating solutions at the bifurcations. The results are applied to networks of coupled Lorenz system. In particular it is shown that, independently of the network topology, for sufficiently large coupling gains the distribution of the delay has a stabilizing effect on the stability of the synchronized equilibrium.

Synchronization of coupled nonlinear oscillators with shifted gamma-distributed delays I. INTRODUCTION

Systems of coupled nonlinear oscillators have been used to describe a variety of phenomena in different fields of science ranging from electric network modeling to brain activity behavior. The emergence of synchronization in networks of oscillators with couplings given by various types of topologies has received much attention [START_REF] Kuramoto | Chemical Oscillations, Waves and Turbulence[END_REF], [START_REF] Pecora | Master stability functions for synchronized coupled systems[END_REF], [START_REF] Strogatz | Sync: the emerging science of spontaneous order[END_REF]. The influence of discrete fixed delays on the consensus and synchronization of coupled subsystems is considered in [START_REF] Oguchi | Synchronization in networks of chaotic systems with time-delay coupling[END_REF], [START_REF] Michiels | Synchronization of delay-coupled nonlinear oscillators: An approach based on the stability analysis of synchronized equilibria[END_REF], [START_REF] Morȃrescu | Consensus with constrained convergence rate and time-delays[END_REF]. However, as pointed out in several studies [START_REF] Macdonald | Biological delay systems: linear stability theory[END_REF], [START_REF] Insperger | Semi-discretization of time-dely systems[END_REF], [START_REF] Michiels | Consensus problems with distributed delays, with application to traffic flow models[END_REF], the behavior of interconnected systems is often more accurately reflected by considering distributed delays affecting their interactions.

This paper focuses on the synchronization problem in networks of identical oscillators with couplings affected by distributed delays. Throughout the study we consider a set of p identical oscillators coupled via a fixed, directed network topology. Each oscillator has a behavior described by the following nonlinear dynamics

ẋi (t) = f (x i (t)) + Bu i (t), y i (t) = Cx i (t), i = 1, . . . , p (1)
where x i ∈ R nx , i = 1, . . . , p, B, C ∈ R nx×m , and f : R n → R n is twice continuously differentiable.

Assumption 1: We assume that for u i = 0 the system has at least one unstable equilibrium of focus type, which we denote by x * in what follows. The network topology is described by the directed graph

G(V, E, G) (2)
characterized by the node set V = {1, . . . , p}, a set of edges E where (k, l) ∈ E if and only if α k,l = 0, and a weighted adjacency matrix G with zero diagonal entries and non-diagonal entries equal to α k,l . Obviously, each vertex corresponds to an oscillator and the edges describe the coupling strength between them. Precisely, we couple the systems (1) by means of the control law

u i (t) = k p l=1 α i,l ∞ 0 g(θ)(y l (t -θ) -y i (t))dθ , i = 1, . . . , p, (3)
where k > 0 represents the gain parameter and g(•) stands for the gamma-distributed delay with a gap kernel given by

g(ξ) =    0, ξ < τ, (ξ-τ ) n-1 e - ξ-τ T T n (n-1)! , ξ ≥ τ, (4)
where n ∈ N, T > 0 and τ ≥ 0. Note that g(ξ) ≥ 0 for all ξ ≥ 0 and ∞ 0 g(ξ)dξ = 1. A motivation for the consideration of this type of delay distribution can be found in [START_REF] Morarescu | Stability crossing curves of shifted gamma-distributed delay systems[END_REF].

The coupled system (1)-( 3) has an equilibrium

(x * , . . . , x * ) ∈ R n×p , (5)
which we denote the synchronized equilibrium. In what follows we perform a stability analysis of the equilibrium [START_REF] Michiels | Synchronization of delay-coupled nonlinear oscillators: An approach based on the stability analysis of synchronized equilibria[END_REF] in the parameter space (k, T, τ ), characterize the emanating solutions in the bifurcations and apply these results to networks of coupled Lorenz systems.

II. PRELIMINARIES

The following assumptions are supposed to hold. Assumption 2: The graph G is strongly connected. Assumption 3: The adjacency matrix G is stochastic i.e.

p l=1 α i,l = 1, i = 1, . . . , p.
The first assumption is natural in the context of synchronization while the second assumption will be motivated in the sequel.

We recall that Perron-Frobenius theorem for stochastic matrices assures that 1) G has a simple eigenvalue equal to 1 and [1, . . . , 1] is the corresponding right-eigenvector.

2) The left eigenvector γ = [γ 1 , . . . , γ p ] associated with the eigenvalue 1 may be chosen such that

γ i > 0, 1 ≤ i ≤ p and p i=1 γ i = 1.
3) All eigenvalues of G have modulus smaller than or equal to 1. In what follows we denote the eigenvalues of G as λ i (G), i = 1, . . . , p, where we take the following convention.

Convention 1: λ 1 (G) = 1. Following the coordinate transformation in [START_REF] Michiels | Synchronization of delay-coupled nonlinear oscillators: An approach based on the stability analysis of synchronized equilibria[END_REF] we can bring the system (1) and (3) in the form

ẋ1 (t) = f (x 1 (t)) + BkC p l=1 α 1,l ∞ 0 g(θ)(x 1 (t -θ) -x 1 (t))dθ + BkC p l=1 α 1,l ∞ 0 g(θ)e l (t -θ)dθ,    ė2 (t)
. . .

ėp (t)    =    f (x 1 (t) + e 2 (t)) -f (x 1 (t)) . . . f (x 1 (t) + e p (t)) -f (x 1 (t))    - k                   p l=1 α 2,l . . . p l=1 α p,l          ⊗ BC             e 2 (t) . . . e p (t)    + k G ⊗ BC    ∞ 0 g(θ)e 2 (t -θ)dθ . . . ∞ 0 g(θ)e p (t -θ)dθ    -          k          p l=1 α 1,l - p l=1 α 2,l . . . p l=1 α 1,l - p l=1 α p,l          ⊗ BC          • ∞ 0 g(θ)(x 1 (t -θ) -x 1 (t))dθ, (6)
where e i = x i -x 1 , i = 2, . . . , p. As in [START_REF] Michiels | Synchronization of delay-coupled nonlinear oscillators: An approach based on the stability analysis of synchronized equilibria[END_REF], Assumption 3 guarantees the existence of synchronized solutions (e 2 = . . . = e p = 0 implies ė2 = . . . ėp = 0).

III. STABILITY ANALYSIS OF SYNCHRONIZED

EQUILIBRIA

Linearizing the system (1) and (3) around the synchronized equilibrium (x * , . . . , x * ), we obtain    v1 (t) . . .

vp (t)    = I ⊗ (A -BkC)    v 1 (t) . . . v p (t)    + kG ⊗ BC    ∞ 0 g(θ)v 1 (t -θ)dθ . . . ∞ 0 g(θ)v p (t -θ)dθ    , (7)
where A = ∂f ∂x (x * ).

A. Characteristic equation

The characteristic function of ( 7) is given by

f (s, k, T, τ ) = det F (s, k, T, τ ), (8)
where the characteristic matrix F is defined by

F (s, k, T, τ ) = I ⊗ (sI -A + kBC) -G ⊗ kBC e -sτ (1 + sT ) n .
(9) Using the factorization G = P ΛP -1 , where Λ ∈ C p×p is triangular and P ∈ C p×p , the characteristic function may be decomposed as

f (s, k, T, τ ) = p i=1 f i (s, k, T, τ ), (10)
where

f i (s, k, T, τ ) = det F i (s, k, T, τ ), F i (s, k, T, τ ) = sI -A + kBC -kBCλ i (G) e -sτ (1 + sT ) n , i = 1, . . . , p. (11)
We note that each function f i (s, k, T, τ ) may be rewritten

f i (s, k, T, τ ) = p l=0 Q l (s, k) e -sτ (1 + sT ) n l (12)
where Q l , l = 1, . . . , p, are polynomials satisfying

deg Q 0 (s) > deg Q l (s), l = 1, . . . , p. (13)

B. Eigenspaces and behavior on the onset of instability

We investigate the eigenspace of the characteristic matrix (9), corresponding to a characteristic root. For reasons of simplicity we restrict ourselves to the generic case where all the eigenvalues of the adjacency matrix G are simple. Let E i be the eigenvector of G corresponding to the eigenvalue λ i (G), i = 1, . . . , p.

If for some l ∈ {1, . . . , p}, the equation

f l (s, k, T, τ ) = 0 has a simple root at s = ŝ such that F l (ŝ, k, T, τ ) V = 0, V ∈ C n×1 , (14)
then it can be verified that

F (ŝ, k, T, τ ) (E l ⊗ V ) = 0. ( 15
)
This implies that the linearized system (7) has an exponential solution

   ν 1 (t) . . . ν p (t)    = c(E l ⊗ V )e λt = c    e l,1 V . . . e l,p V    e λt , (16)
with the constant c depending on the initial conditions. In words, in an exponential solution of ( 7) corresponding to a zero of f l (s, k, T, τ ), the relation between the state variables of an individual subsystem is determined by the null vector of F l , while the relation between the corresponding state variables of the different subsystems is solely determined by the eigenvector E l , corresponding to the l-th eigenvalue of the adjacency matrix G. This implies that all modes can be classified in at most p types, based on the relations between the behavior of the different subsystems. Note for example that bifurcations of the synchronized equilibrium due to imaginary axis zeros of f 1 are bifurcations where the synchronized behavior is preserved in the emanating solution, as

E 1 = [1 1 • • • 1] T .

C. Stability crossing curves in the delay parameters space

For a fixed k = k * we describe the stability crossing curves T , which is the set of (T, τ ) values such that (8) has imaginary axis solutions.

As the parameters (T, τ ) cross the stability crossing curves, some characteristic roots cross the imaginary axis. Therefore, the number of roots on the right half complex plane are different on the two sides of the curves, from which we may describe the parameter regions of (T, τ ) in R 2 + for the system to be stable.

Another related useful concept is the frequency crossing set Ω, which is defined as the collection of all ω > 0 such that there exists a parameter pair (T, τ ) such that f (jω, k * , T, τ ) = 0. In other words, as the parameters T and τ vary, the characteristic roots may cross the imaginary axis at jω if and only if ω ∈ Ω.

Proposition 1: The crossing set Ω consists of a finite number of bounded intervals, Ω 1 , Ω 2 , ..., Ω N .

Proof: Due to decomposition [START_REF] Strogatz | Nonlinear dynamics and chaos: with applications to physics, biology, chemistry and engineering[END_REF] we may reduce the study of f (jω, k * , T, τ ) = 0 to the study of f i (jω, k * , T, τ ) = 0. Therefore we have to study the crossing set of functions given by (12). Since 1 + jωT = 0 the equation f i (jω, k * , T, τ ) = 0 can be rewritten as

Q 0 (jω, k * )(1+jωT ) pn = - p l=1 Q l (s, k * )e -jlωτ (1+jωT ) (p-l)n
Taking the modulus in the previous equation one gets

| Q 0 (jω, k * ) | | 1+jωT | pn ≤ p l=1 | Q l (s, k * ) | | 1+jωT | (p-l)n
Using | 1 + jωT |≥ 1 one can multiply each terms in the right hand side with powers of

| 1 + jωT | to obtain | Q 0 (jω, k * ) |≤ p l=1 | Q l (s, k * ) | . ( 17
)
In conclusion ω ∈ Ω implies that (17) is satisfied. The assertion follows from (13). Proposition 2: Given any ω > 0, ω ∈ Ω, all the corresponding T, τ satisfying f i (jω, k * , T, τ ) = 0 can be computed by

T = T l = 1 ω | z l (ω) | -2/n -1 1/2 (18) τ = τ l,m = 1 ω (-∠z l (ω) -n arctan(ωT ) + 2mπ), m = 0, ±1, ±2, .... ( 19
)
where

z l (ω) satisfies | z l (ω) |< 1 and is a solution of p =0 Q (jω, k * )z = 0. (20)
Proof:

Let us consider the following variable transformation z = e -jωτ (1 + jωT ) n . Then, the equation f i (jω, k * , T, τ ) = 0 rewrites as

p l=0 Q l (jω, k * )z l = 0
Solving this equation in z we get the solutions z 1 (ω), . . . , z h (ω). From definition, it is clear that | z |< 1 for any value of ω > 0. Therefore, the admissible solutions are given by those z l (ω), l = 1, . . . , h with | z l (ω) |< 1.

For such z l the corresponding pair (T, τ ) can be simply computed by ( 18) and ( 19). Note that z solving (20) can be equivalently expressed as z being an eigenvalue of the pencil

Λ(z) := (-jωI + A -k * BC) + z(k * BC).
We will not restrict ∠z l (ω) to a 2π range. Rather, we allow it to vary continuously within each interval Ω i . Thus, for each fixed m and l, (18) and ( 19) represent a continuous curve. We denote such a curve as T l,m i . Therefore, corresponding to a given interval Ω i , we have an infinite number of continuous stability crossing curves T l,m i , m = 0, ±1, ±2, .... It should be noted that, for some m, part or the entire curve may be outside of the range R 2 + and, therefore, may not be physically meaningful.

The collection of all the points in T corresponding to Ω i may be expressed as

T i = +∞ m=-∞ h l=1 T l,m i R 2 + .
Obviously, T = N i=1 T i .

D. Computation of the crossing direction

Let jω * be an imaginary axis characteristic root for parameter values (T, τ ) = (T * , τ * ), which are lying on the stability crossing curve T l,m i . Let r be such that f r (jω * , k * , T * , τ * ) = 0. Introduce the following notations:

R 0 + jI 0 = ∂fr(s,k * ,T,τ ) ∂s s=jω * ,T =T * ,τ =τ * , R 1 + jI 1 = ∂fr(s,k * ,T,τ ) ∂T s=jω * ,T =T * ,τ =τ * , R 2 + jI 2 = ∂fr(s,k * ,T,τ ) ∂τ s=jω * ,T =T * ,τ =τ * . If matrix M := R 1 R 2 I 1 I 2
is invertible then, by the implicit function theorem the equations

(f r (σ + jω, k * , T, τ )) = 0 (f r (σ + jω, k * , T, τ )) = 0
locally define a function (σ, ω) → (T, τ ) mapping (0, ω * ) to (T * , τ * ) and satisfying

∂T ∂σ (0, ω * ) ∂τ ∂σ (0, ω * ) = -M -1 R 0 I 0 , ∂T ∂ω (0, ω * ) ∂τ ∂ω (0, ω * ) = -M -1 -I 0 R 0 . (21)
We will call the direction of the curve T l,m i that corresponds to increasing ω the positive direction. Since the tangent of . Corresponding to a pair of complex conjugate solutions of (8) crossing the imaginary axis along the horizontal direction, (T, τ ) moves along the direction ∂T ∂σ , ∂τ ∂σ . So, as (T, τ ) crosses the stability crossing curve from the right hand side to the left hand side, a pair of complex conjugate solutions of ( 8) cross the imaginary axis to the right half plane, if

∂T ∂ω ∂τ ∂σ - ∂τ ∂ω ∂T ∂σ s=jω > 0, (22)
i.e. the region on the left of T l,m i gains two solutions on the right half plane. If the inequality ( 22) is reversed then the region on the left of T l,m i loses two right half plane solutions. Using (21) a direct computation shows that

∂T ∂ω ∂τ ∂σ - ∂τ ∂ω ∂T ∂σ s=jω = (R 2 0 + I 2 0 )(R 2 I 1 -R 1 I 2 ) (R 1 I 2 -I 1 R 2 ) 2 .
This brings us to the following result. Proposition 3: Let ω ∈ Ω r ⊂ Ω and (T, τ ) ∈ T i such that jω is a simple solution of the characteristic equation f r (s, k * , T, τ ) = 0. Then as (T, τ ) crosses a stability crossing curve T l,m i from the right-hand side to the left-hand side, a pair of solutions of (8) crosses the imaginary axis to the right, through the gate s

= ±jω if R 2 I 1 -R 1 I 2 > 0.
The crossing is to the left if the inequality is reversed.

IV. APPLICATION TO COUPLED LORENZ SYSTEMS

In this section the nonlinear oscillators (1) are specified as Lorenz systems: (24)

   ẋi,1 (t) = σ(x i,2 (t) -x i,1 (t)), ẋi,2 (t) = rx i,1 (t) -x i,2 (t) -x i,1 (t)x i,3 (t) + u i,1 (t), ẋi,3 (t) = -bx i,3 (t) + x i,1 (t)x i,2 (t) + u i,2 (t), y i,1 (t) = x i,2 (t), y i,2 (t) = x i,3 (t) -r, i = 1, . . . p, ( 23
) where u i = [u i,1 u i,2 ] T , y i = [y i,1 y i,2 ] T .
Note that for u i ≡ 0 each Lorenz system has three equilibria given by

(0, 0, 0), ± b(r -1), ± b(r -1), r -1 , (25)
the latter two corresponding to unstable foci. Thus, Assumption 1 holds. Furthermore, with the parameter values (24) it exhibits a chaotic attractor [START_REF] Strogatz | Nonlinear dynamics and chaos: with applications to physics, biology, chemistry and engineering[END_REF].

If we linearize the coupled system ( 23) and ( 3) around the synchronized equilibrium (x * , . . . , x * ), x * = (± b(r -1), ± b(r -1), r -1),

(26) then we obtain the linear system [START_REF] Morȃrescu | Consensus with constrained convergence rate and time-delays[END_REF], where the matrices are specified as

A =   -σ σ 0 1 1 ∓ b(r -1) ± b(r -1) ± b(r -1) -b   , B = C =   0 0 1 0 0 1   .
It is easy to show that the stability of the linearized system does not depend on which equilibrium x * in ( 26) is considered, Therefore, we will restrict ourselves to the one in the positive octant.

In what follows we analyze the stability properties of the synchronized equilibria (26) in the (k, T, τ ) parameter space. First we study the asymptotic behavior for large values of the gain parameters in §IV-A. For the standard parameters (24) this will allow us to make assertions about stability regions, stability switches and emerging behavior, which do not depend on the network topology. Next we present a numerical example in §IV-B.

A. Asymptotic behavior for large gain values

For T = 0 the behavior of the stability crossing curves for large gain values is described by the following result from [START_REF] Michiels | Synchronization of delay-coupled nonlinear oscillators: An approach based on the stability analysis of synchronized equilibria[END_REF].

Theorem 1: Consider a network of coupled Lorenz systems (23) with parameters (24) and coupling [START_REF] Macdonald | Biological delay systems: linear stability theory[END_REF]. Assume that T = 0. Then there exists a number k > 0 and a function

τ * : [ k, ∞] → R + , k → τ * (k), (27)
satisfying the following properties: 1) there is a constant k > k such that for every k > k, the synchronized equilibrium has two characteristic roots in the open right half plane for all τ ∈ [0, τ * ], while it is asymptotically stable for τ ∈ (τ * , τ * + ), with sufficiently small; 2) at τ = τ * a synchronization preserving Hopf bifurcation occurs; 3) for all k ∈ [ k, ∞] we can factor

τ * (k) = υ(k) k where lim k→∞ υ(k) = 0.586004. (28)
Furthermore, the number k and the function ( 27) are independent of the number of subsystems and of the network topology.

The following theorem shows that the qualitative behavior for large values of k is significantly different if T = 0.

Theorem 2: Consider a network of coupled Lorenz systems (23) with parameters (24) and coupling [START_REF] Macdonald | Biological delay systems: linear stability theory[END_REF]. Assume that T > 0. Then there is a constant k such that for every k > k, the synchronized equilibrium is asymptotically stable independently of τ ≥ 0.

Proof: Part 1: We first show that for τ = 0 all characteristic roots are in left half plane for sufficiently large k.

Let i ∈ {1, . . . , p}. As k → ∞ the function

f i (s, k, T, 0) k 2 = det N i (s, k), with N i (s, k) =     s + σ -σ 0 -1 k s+1 k -λ i (G) (1+st) n -1 √ b(r-1) k - √ b(r-1) k - √ b(r-1) k s+b k - λ i (G) (1+sT ) n -1     ,
uniformly converges on compact subsets of C to the function

1 -λi(G) (1+sT ) n 2 (s + σ). (29)
From Rouché's theorem it follows that, as k → ∞, 2n + 1 zeros of f i (s, k, T, 0) converge pair-wise to the 2n + 1 zeros of (29) (taking multiplicity into account). Similarly, letting

s = s k we get f i (sk, k, T, 0) k 3 = det(sI -P i (s, k)), where P i (s, k) =     -σ k σ k 0 1 k -1 k + λ i (G) (1+skT ) n -1 - √ b(r-1) k √ b(r-1) k √ b(r-1) k -b k + λ i (G) (1+skT ) n -1     .
Thus, fi(sk,k,T,0) k 3

uniformly converges on compact sets in the s-plane to s(s + 1) 2 . Rouché's theorem then implies that for sufficiently large k the function f i (s, k, T, 0) has two zeros equal to k. The function f i (s, k, T, 0) has exactly 2n + 3 zeros, induced by rank(BC) = 2. We conclude from the previous observations that all these zeros are in the open left half plane for large k if the zeros of (29) are in the open left half plane. This is the case for i ∈ {2, . . . , p}. The case i = 1 deserves special attention, since for λ 1 (G) = 1, (29) has two (rightmost) zeros equal to zero. The latter means that two zeros of f 1 (s, k, T, 0) approach the origin as k → ∞. It remains to show that these approach the imaginary axis from the left.

With ŝ = ks we define g(ŝ, k) := f 1 ŝ k , k, T, 0 . By the expansion

1 (1 + ŝT /k) n -1 = - ŝT n k + O ŝT n k 2
and, again, an application of Rouché's theorem, it follows that, as k → ∞, two zeros of g(•, k) converge to the finite eigenvalues of the pencil -A + ŝ(T n)BC. With the parameters (24) these eigenvalues are given by ŝ1,2 = 1 T n (-1.333 ± 1.926j). We conclude that the two zeros of f 1 (s, k, T, τ ) that converge to the origin as k → ∞ satisfy

s 1,2 = 1 k ŝ1,2 + O 1 k 2 ,
hence, (s 1,2 ) < 0 for large k. Part 2: We show that for large values of k eigenvalues cannot cross the imaginary axis when varying parameter τ . The theorem then follows from the stability for τ = 0.

The equation

f i (jω, k, T, τ ) = 0 is equivalent to det I -(jωI -A + kBC) -1 kBCλ i (G) e -jωτ (1 + jωT ) n = 0.
A necessary solvability conditions is given by

ρ (jωI -A + kBC) -1 kBC λ i (G) (1 + jωT ) n = 1,
with ρ(•) the spectral radius. This condition is always violated for large k. Indeed, in the complex plane the nonzero eigenvalues of the matrix

j ω k I - A k + BC -1 BCλ i (G), (30)
converge to the curve

Ω ≥ 0 → 1 jΩ + 1 λ i (G)
as k → ∞, uniformly in the parameter ω ≥ 0. Furthermore, we have:

1 1 + jΩ λ i (G) ≤ |λ i (G)| ≤ 1, ∀Ω ≥ 0 and 1 |(1 + jωT ) n | < 1, ∀ω > 0.
It follows that zeros on the imaginary axis are not possible for large values of k.

B. Numerical example

We consider a ring topology with unidirectional coupling, described by the adjacency matrix

G =      0 • • • 0 1 1 0 . . . . . . 1 0      ∈ R p×p , (31)
which has the following properties:

λ l (G) = e j 2π(l-1) p , E l =       1 e -j 2π(l-1) p . . . e -j 2π(p-1)(l-1) p      
for l = 1, . . . , p. If ( 14) is satisfied for ŝ = jω, ω > 0, then the emanating solution (16) becomes

   ν 1 (t) . . . ν p (t)    = c         V e jωt V e jωt-2π(l-1) p V e jωt-4π(l-1) p . . . V e jωt-2(p-1)π(l-1) p         . ( 32
)
It can be interpreted as a traveling wave solution, where the agents follow each other with a phase shift of 360(l -1)/p degrees. Therefore, if the characteristic root ŝ on the imaginary axis corresponds to a Hopf bifurcation of the original nonlinear system (23) and ( 3) for a critical value of some free parameter, we refer to this bifurcation as a "Hopf 360(l -1)/p" bifurcation. For p = 4 subsystems the stability regions in the (k, τ ) parameter space are displayed in Figure 1, for T = 0 (left pane) and T = 0.01 (right pane). These results are in accordance with Theorems 1-2. Figure 2 is obtained by zooming in on Figure 1 (right). Finally, Figure 3 shows the stability regions in the (T, τ ) parameter space for k = 12 and k = 25, computed by the approach of Section III-C.

V. CONCLUSION

We gained qualitative and quantitative insight in the synchronization of delay-coupled nonlinear oscillations, where the latency in the coupling is modeled by shifted gammadistributed delays, by studying the stability regions of the synchronized equilibrium in the (coupling gain-delay) parameter space and by characterizing the structure of the emanating solutions in the bifurcations. Applying the results to the particular case of coupled Lorenz systems we have shown that, independently of the network topology, for sufficiently large coupling gains the distribution of the delay has a stabilizing effect on the stability of the synchronized equilibrium. The fact that the synchronization mechanism is a high-gain mechanism, based on dominating the destabilizing terms in the synchronization error dynamics, suggests that a similar conclusion can be drawn for other types of nonlinear oscillators.

  The parameter values are given by σ = 10, r = 28, b = 8/3.

  s 1,2 , where lim k→∞ s1,2 = -1, hence they move off the infinity without leaving the open left half plane.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. Stability regions of synchronized equilibria of 4 coupled Lorenz systems in a ring configuration. The numbers refer to the number of characteristic roots in the closed right half plane. The Hopf bifurcation curves are classified based on type of emanating solution (32).