Anne-Laure Basdevant
email: .anne-laure.basdevant@ens.fr

Arvind Singh
email: arvind.singh@ens.fr

On the speed of a cookie random walk

Keywords: Law of large numbers, cookie or multi-excited random walk, branching process with migration

We consider the model of the one-dimensional cookie random walk when the initial cookie distribution is spatially uniform and the number of cookies per site is finite. We give a criterion to decide whether the limiting speed of the walk is non-zero. In particular, we show that a positive speed may be obtained for just 3 cookies per site. We also prove a result on the continuity of the speed with respect to the initial cookie distribution.

Introduction

We consider the model of the multi-excited random walk, also called cookie random walk, introduced by Zerner in [START_REF] Zerner | Multi-excited random walks on integers[END_REF] as a generalization of the model of the excited random walk described by Benjamini and Wilson in [START_REF] Benjamini | Excited random walk[END_REF] (see also Davis [3] for a continuous time analogue). The aim of this paper is to study under which conditions the speed of a cookie random walk is strictly positive. In dimension d ≥ 2, this problem was solved by Kozma [START_REF] Kozma | Excited random walk in three dimensions has positive speed[END_REF][START_REF] Kozma | Excited random walk in two dimensions has linear speed[END_REF] who proved that the speed is always non-zero. In the one-dimensional case, the speed can either be zero or strictly positive. We give here a necessary and sufficient condition to determine if the walk's speed is strictly positive when the initial cookie environment is deterministic, spatially uniform and with a finite number of cookies per site. Let us start with an informal definition of such a process:

Let us put M ≥ 1 "cookies" at each site of Z and let us pick p 1 , p 2 , . . . , p M ∈ [ 1 2 , 1). We say that p i represents the "strength" of the i th cookie at any given site. Then, a cookie random walk X = (X n ) n≥0 is simply a nearest neighbour random walk, eating the cookies it finds along its path by behaving in the following way:

• If X n = x and there is no remaining cookie at site x, then X jumps at time n + 1 to x + 1 or x -1 with equal probability 1 2 . • If X n = x and there remain the cookies with strengths p j , p j+1 , . . . , p M at this site, then X eats the cookie with attached strength p j (which therefore disappears from this site) and then jumps at time n + 1 to x + 1 with probability p j and to x -1 with probability 1 -p j .

This model is a particular case of self-interacting random walk: the position of X at time n + 1 depends not only of its position at time n but also on the number of previous visits to its present site. Therefore, X is not a Markov process.

Let us now give a formal description of the general model. We define the set of cookie environments by Ω = [ 1 2 , 1] N * ×Z . Thus, a cookie environment is of the form ω = (ω(i, x)) i≥1,x∈Z where ω(i, x) represents the strength of the i th cookie at site x. Given x ∈ Z and ω ∈ Ω, a cookie random walk starting from x in the cookie environment ω is a process (X n ) n≥0 on some probability space (Ω, F , P ω,x ) such that:

   P ω,x {X 0 = z} = 1, P ω,x {|X n+1 -X n | = 1} = 1, P ω,x {X n+1 = X n + 1 -X 1 , . . . , X n } = ω(j, X n ) where j = ♯{0 ≤ i ≤ n , X i = X n }.
In this paper, we restrict our attention to the set of environments Ω u M ⊂ Ω which are spatially uniform with at most M ≥ 1 cookies per site:

ω ∈ Ω u M ⇐⇒    for all x ∈ Z and all i ≥ 1 ω(i, x) = ω(i, 0), for all i > M ω(i, 0) = 1 2 , for all i ≥ 1 ω(i, 0) < 1.
The last condition ω(i, 0) < 1 is introduced only to exclude some possible degenerated cases but can be relaxed (see Remark 2.4). A cookie environment ω ∈ Ω u M may be represented by (M, p) where p = (p 1 , . . . , p M ) = (ω(1, 0), . . . , ω(M, 0)).

In this case, we shall say that the associated cookie random walk is an (M, p)-cookie random walk and we will use the notation P (M,p) instead of P ω .

The question of the recurrence or transience of a cookie random walk was solved by Zerner in [START_REF] Zerner | Multi-excited random walks on integers[END_REF] for general cookie environments (even in the case where the initial cookie environment may itself be random). In particular, he proved that, if X is an (M, p) cookie random walk, there is a phase transition according to the value of

α = α(M, p) def = M i=1 (2p i -1) -1.
(1)

• If α ≤ 0 then the walk is recurrent i.e. lim sup X n = -lim inf X n = +∞ a.s.

• If α > 0 then X is transient toward +∞ i.e lim X n = +∞ a.s.
In particular, for M = 1, the cookie random walk is always recurrent for any choice of p. However, as soon as M ≥ 2, the cookie random walk can either be transient or recurrent depending on p. Zerner [START_REF] Zerner | Multi-excited random walks on integers[END_REF] also proved that the speed of a (M, p)-cookie random walk X is always well defined (but may or not be zero). Precisely,

• there exists a constant v(M, p) ≥ 0 such that

X n n -→ n→∞ v(M, p) P (M,p) -almost surely.
• The speed is monotonic in p: if p = (p 1 , . . . , p M ) and q = (q 1 , . . . , q M ) are two cookie environments such that p i ≤ q i for all i, then v(M, p) ≤ v(M, q).

• The speed of a (2, p)-cookie random walk is always 0.

The question of whether one can construct a (M, p)-cookie random walk with strictly positive speed was affirmatively answered by Mountford, Pimentel and Valle [START_REF] Mountford | On the speed of the one-dimensional excited random walk in the transient regime[END_REF] who considered the case where all the cookies have the same strength p ∈ [ 1 2 , 1) i.e. the cookie vector p has the form [p] M def = (p, . . . , p). They showed that:

• For any p ∈ ( 1 2 , 1), there exists an M 0 such that for all M > M 0 the speed of the (M, [p] M )-cookie random walk is strictly positive.

• If M(2p -1) < 2, then the speed of the (M, [p] M )-cookie random walk is zero.

They also conjectured that when M(2p -1) > 2, the speed should be non-zero. The aim of this paper is to prove that such is indeed the case.

Theorem 1.1. Let X denote a (M, p)-cookie random walk, then

lim n→∞ X n n = v(M, p) > 0 ⇐⇒ α(M, p) > 1
where α(M, p) is given by [START_REF] Antal | The excited random walk in one dimension[END_REF].

In particular, we see that a non-zero speed may be achieved for as few as 3 cookies per site. Comparing this result with the transience/recurrence criteria, we have a second order phase transition at the critical value α = 1. In fact, it shall be proved in a forthcoming paper that, for 0 < α < 1, the rate of transience of X n is of order n α+1 2 . One would certainly like an explicit calculation of the limiting velocity in term of the cookie environment (M, p) but this seems a challenging problem (one can still look at the end of Section 3 where we give an implicit formula for the speed). However, one can prove that the speed is continuous in p and has a positive right derivative at all its critical points: Theorem 1.2.

• For each M, the speed v(M, p) is a continuous function of p in Ω u M . • For any environment (M, pc ) with α(M, pc ) = 1, there exists a constant C > 0 (depending on (M, pc )) such that

lim p→pc p∈Ω u M α(p)>1 v(M, p) α(M, p) -1 = C.
In particular, for M ≥ 3, the (unique) critical value for an (M, [p] M )-cookie random walk is p c = 1 M + 1 2 and the function v(p) is continuous, non-decreasing, zero for p ≤ p c , and admits a finite strictly positive right derivative at p c .

The remainder of this paper is organized as follow. In the next section, we construct a Markov process associated with the hitting time of the cookie random walk. The method is similar to that used by Kesten, Kozlov and Spitzer [START_REF] Kesten | A limit law for random walk in a random environment[END_REF] for the determination of the rates of transience of a random walk in a one-dimensional random environment. It turns out that, in our setting, the resulting process is a branching process with random migration. The study of this process and of its stationary distribution is done in Section 3. This enables us to complete the proof of Theorem 1.1. Finally, the last section is dedicated to the proof of Theorem 1.2.

An associated branching process with migration

In the remainder of this paper, X = (X n ) n≥0 will denote a (M, p)-cookie random walk. Since the speed of a recurrent cookie random walk is zero, we will also assume that we are in the transient regime i.e.

α(M, p) = M i=1 (2p i -1) -1 > 0.
(

) 2
For the sake of brevity, we simply write P x for P (M,p),x and P instead of P 0 (the process starting from 0). Let T n stand for the hitting time of level n ≥ 0 by X:

T n = inf(k ≥ 0 , X k = n). (3)
For 0 ≤ k ≤ n, let U n i denote the number of jumps of the cookie random walk from site i to site i -1 before reaching level n

U n i = ♯{0 ≤ k < T n , X k = i and X k+1 = i -1}.
Let also K n stand for the total time spent by X in the negative half-line up to time

T n K n = ♯{0 ≤ k ≤ T n , X k < 0}.
A simple combinatorial argument readily yields

T n = K n -U n 0 + n + 2 n k=0 U n k .
Notice that, as n tends to infinity, the random variable K n increases almost surely toward K ∞ , the total time spent by the cookie random walk in the negative half line. Similarly,

U n 0 increases toward U ∞ 0 the total number of jumps from 0 to -1. Since X is transient, K ∞ + U ∞
0 is almost-surely finite and therefore

T n ∼ n→∞ n + 2 n k=0 U n k . (4)
Let us now prove that for each n, the reverse process (U n n , U n n-1 , . . . , U n 1 , U n 0 ) has the same law as the n first steps of some branching process Z with random migration. We first need to introduce some notations. Let (B i ) i≥1 denote a sequence of independent Bernoulli random variable under P with distribution:

P{B i = 1} = 1 -P{B i = 0} = p i if i ≤ M, 1 2 if i > M. (5)
For j ∈ N, define

k j = min(k ≥ 1, ♯{1 ≤ i ≤ k, B i = 1} = j + 1)
and

A j = ♯{1 ≤ i ≤ k j , B i = 0} = k j -j -1.
We have the following easy lemma:

Lemma 2.1.

• For any i, j ≥ 0, we have P{A j = i} > 0.

• For all j ≥ M, we have

A j law = A M -1 + ξ 1 + . . . + ξ j-M +1 (6)
where

(ξ i ) i≥0 is a sequence of i.i.d. geometrical random variable with parameter 1 2 independent of A M -1 .
Proof. The first part of the lemma is a direct consequence of the assumption that p is such that p k < 1 for all k. To prove the second part, we simply notice that k M -1 ≥ M so that for j ≥ M, the random variable A j -A M -1 has the same law as the random variable

min(k ≥ 1, ♯{1 ≤ i ≤ k, B i = 1} = j + 1 -M) -j -1 + M (7)
where ( B i ) i≥0 is a sequence of i.i.d. random variables independent of A M -1 and with common Bernoulli distribution

P{ B i = 0} = P{ B i = 1} = 1 2 .
It is clear that [START_REF] Kozma | Excited random walk in two dimensions has linear speed[END_REF] has the same law as ξ 1 + . . . + ξ j-M +1 .

By possibly extending the probability space, we now construct a process Z = (Z n , n ≥ 0) and a family of probability (P z ) z≥0 such that, under P z , the process Z is a Markov chain starting from z, with transition probability:

P z {Z 0 = z} = 1, P z {Z n+1 = k | Z n = j} = P{A j = k}.
Since the family of probabilities (P z ) depends on the law of the cookie environment (M, p), we should rigourously write P (M,p),z instead of P z . However, when there is no possible confusion we will keep using the abbreviated notation. Furthermore, we will simply write P for P 0 and E will stand for the expectation with respect to P.

Let us now notice that, in view of the previous lemma, Z n under P z may be interpreted as the number of particles alive at time n of a branching process with random migration starting from z, that is a branching process which allows immigration and emigration (see Vatutin and Zubkov [START_REF] Vatutin | Branching processes[END_REF] for a survey on these processes). Indeed:

• If Z n = j ≥ M -1,
then according to Lemma 2.1, Z n+1 has the same law as j-M +1 k=1 ξ k + A M -1 , i.e. M -1 particles emigrate and the remaining particles reproduce according to a geometrical law with parameter 1 2 and there is also an immigration of A M -1 new particles.

• If Z n = j ∈ {0, . . . , M -2} then Z n+1 has the same law as A j i.e. all the j particles emigrate and A j new particles immigrate.

We can now state the main result of this section:

Proposition 2.2. For each n ∈ N, (U n n , U n n-1 , . . . , U n 0 ) under P has the same law as

(Z 0 , Z 1 , . . . , Z n ) under P.
Proof. The argument is similar to the one given by Kesten et al. in [START_REF] Kesten | A limit law for random walk in a random environment[END_REF]. Recall that U n i represents the numbers of jumps of the cookie random walk X from i to i -1 before reaching n. Then, conditionally on (U n n , U n n-1 , . . . , U n i+1 ), the number of jumps U n i from i to i -1 depends only on the number of jumps from i + 1 to i, that is, depends only of U n i+1 . This shows that (U n n , U n n-1 , . . . , U n 0 ) is indeed a Markov process. By definition, Z 0 = 0 P-a.s. and U n n = 0 P-a.s. It remains to compute

P{U n i = k | U n i+1 = j}.
Note that the number of jumps from i to i -1 before reaching level n is equal to the number of jumps from i to i -1 before reaching i + 1 for the first time plus the sum of the number of jumps from i to i -1 between two consecutive jumps from i + 1 to i which occur before reaching level n. Thus, conditionally on {U n i+1 = j}, the random variable U n i has the same law as the number of failures (i.e. B k = 0) in the Bernoulli sequence (B 1 , B 2 , B 3 , . . .) defined by ( 5) before having exactly j + 1 successes. This is precisely the definition of A j and therefore

P{U n i = k | U n i+1 = j} = P j {Z 1 = k}. Since U n
0 is the number of jumps from 0 to -1 of the cookie random walk X before reaching level n and since we assumed that the cookie random walk X is transient, U n 0 increases almost surely toward the total number U ∞ 0 of jumps of X from 0 to -1. In view of the previous proposition, this implies that under P, Z n converges in law toward a random variable which we denote by Z ∞ .

Let us also note that Z is a irreducible Markov chain (this is a consequence of part 1 of Lemma 2.1). Since Z converges in law toward a limiting distribution, this shows that Z is in fact a positive recurrent Markov chain. In particular, Z n converges in law toward Z ∞ independently of its starting point (i.e. the law of Z ∞ is the same under any P x ) and the law of Z ∞ is also the unique invariant probability for Z.

Corollary 2.3. Recall that v(M, p) denotes the limiting speed of the cookie random walk X. We have

v(M, p) = 1 1 + 2E[Z ∞ ]
(with the convention 0 = 1 +∞ ).

In particular, the speed of an (M, p)-cookie random walk is non zero i.i.f. the limiting random variable Z ∞ of its associated process Z has a finite expectation.

Proof. Since X is transient, we have the well known equivalence valid for v ∈ [0, ∞] :

X n n -→ n→∞ v P-a.s. ⇐⇒ T n n -→ n→∞ 1 v P-a.s. (8)
On the one hand, this equivalence and (4) yield

1 n n k=0 U n k -→ n→∞ 1 2v(M, p) - 1 2 P-a.s. (9)
On the other hand, making use of an ergodic theorem for the positive recurrent Markov chains Z with stationary limiting distribution Z ∞ , we find that

1 n n i=1 Z k → n→∞ E[Z ∞ ] P-a.s. ( 10
) (this result is valid even if E[Z ∞ ] = ∞). Proposition 2.
2 implies that the limits in ( 9) and ( 10) are the same. This completes the proof of the corollary.

Remark 2.4. We assumed in the definition of an (M, p) cookie environment that

p i = 1 for all 1 ≤ i ≤ M.
This hypothesis is intended only to ensure that Z starting from 0 is not almost surely bounded (for instance, if p 1 = 1 then 0 is a absorbing state for Z). More generally, one may check from the definition of the random variables A j that Z starting from 0 is almost surely unbounded i.i.f.

♯{1 ≤ j ≤ i , p j = 1} ≤ i 2 for all 1 ≤ i ≤ M. (11)
When this condition fails, Z starting from 0 is almost surely bounded by M -1, thus E[Z ∞ ] < ∞ and the speed of the associated cookie random walk is strictly positive. Otherwise, when [START_REF] Zerner | Recurrence and transience of excited random walks on Z d and strips[END_REF] is fulfilled, Z ultimately hits any level x ∈ N with probability 1 and the proof of Theorem 1.1 remains valid.

3 Study of Z ∞ .

We proved in the previous section that the strict positivity of the speed of the cookie random walk X is equivalent to the existence of a finite first moment for the limiting distribution of its associated Markov chain Z. We shall now show that, for any cookie environment (M, p) (with α(M, p) > 0), we have

E[Z ∞ ] def = E (M,p) [Z ∞ ] < ∞ ⇐⇒ α(M, p) > 1.
This will complete the proof of Theorem 1.1. We start by proving that Z ∞ cannot have moments of any order.

Proposition 3.1. We have

E Z M -1 ∞ = +∞.
Proof. Let us introduce the first return time to 0 for Z:

σ = inf(n ≥ 1 , Z n = 0).
Since Z is a positive recurrent Markov chain, we have 1 ≤ E 0 [σ] < ∞ and the invariant probability measure is given for any y ∈ N by

P{Z ∞ = y} = E 0 σ-1 k=0 1l Z k =y E 0 [σ] .
A monotone convergence argument yields

E 0 σ-1 k=0 Z M -1 k = E 0 [σ]E[Z M -1 ∞ ] (12)
(where both side of this equality may be infinite). We can find n 0 ∈ N * such that P 0 {Z n 0 = M, n 0 < σ} > 0 (in fact, since we assume that p i < 1 for all i, we can choose n 0 = 1). Therefore, making use of the Markov property of Z, we find that

E 0 σ-1 k=0 Z M -1 k ≥ P 0 {Z n 0 = M, n 0 < σ}E M σ-1 k=0 Z M -1 k = P 0 {Z n 0 = M, n 0 < σ} ∞ k=0 E M Z M -1 k∧σ . (13)
In view of ( 12) and ( 13), we just need to prove that

∞ k=0 E M Z M -1 k∧σ = ∞. (14)
We now use a coupling argument. Let us define a new Markov chain Z such that, under P z , the process evolves in the following way

• Z 0 = z, • if Z n = k ∈ {0, 1, . . . , M -1} then Z n+1 = 0, • if Z n = k > M -1 then Z n+1 has the same law as k-(M -1)
i=1 ξ i where (ξ i ) i≥1 is a sequence of i.i.d. geometrical random variables with parameter 1 2 .

Thus, Z is a branching process with emigration: at each time n, there are min( Z n , M -1) particles which emigrate the system and the remaining particles reproduce according to a geometrical law of parameter 1 2 . Recall that Z is a branching process with migration, where at most M -1 particles emigrate at each unit of time, and has the same offspring reproduction law as Z. Therefore, for any z ≥ 0, under P z , the process Z is stochastically dominated by Z. Since 0 is an absorbing state for Z, this implies that, for all n ≥ 0 and all z ≥ 0,

E z [ Z M -1 n ] ≤ E z [Z M -1 n∧σ ]. (15)
Our process Z belongs to the class of processes studied by Kaverin [START_REF] Kaverin | Refinement of limit theorems for critical branching processes with emigration[END_REF]. Moreover, all the hypotheses of Theorem 1 of [START_REF] Kaverin | Refinement of limit theorems for critical branching processes with emigration[END_REF] are clearly fulfilled (in the notation of [START_REF] Kaverin | Refinement of limit theorems for critical branching processes with emigration[END_REF], we have here λ = θ = M -1 and B = 1). Therefore, for any z ≥ M, there exists a constant c > 0 (depending on z) such that

E z [ Z M -1 n ] ∼ n→∞ c n . ( 16
)
The combination of ( 15) and ( 16) yield ( 14).

Remark 3.2. In view of the last proposition and Corollary 2.3, we recover the fact that for M = 2, the speed of the cookie random walk is always zero.

In order to study more precisely the distribution of Z ∞ , we will need the following lemma Lemma 3.3. We have

E [A M -1 ] = 2 M i=1 (1 -p i ).
Proof. Recall that (B i ) i≥1 denotes a sequence of independent Bernoulli random variables with distribution given by [START_REF] Kesten | A limit law for random walk in a random environment[END_REF]. Recall also that

k M -1 = min(k ≥ 1, ♯{1 ≤ i ≤ k, B i = 1} = M), A M -1 = k M -1 -M.
In particular, for any

j ∈ N * P{A M -1 = j} = P{k M -1 = M + j} = P ♯{1 ≤ i ≤ M + j -1, B i = 1} = M -1 P{B M +j = 1}.
Hence,

P{A M -1 = j} = 1 2 M l=1 P ♯{1 ≤ i ≤ M, B i = 1} = M -l P ♯{M +1 ≤ i ≤ M +j -1, B i = 1} = l-1 = M ∧j l=1 P ♯{1 ≤ i ≤ M, B i = 1} = M -l C l-1 j-1 1 2 j . Let L = ♯{1 ≤ i ≤ M, B i = 0}, therefore P{A M -1 = j} = M ∧j l=1 P{L = l}C l-1 j-1 1 2 j for j ∈ N * .
Making use of the relation jC l-1 j-1 = lC l j , we get

E[A M -1 ] = ∞ j=1 M ∧j l=1 P{L = l} 1 2 j lC l j = M l=0 lP{L = l} ∞ j=l 1 2 j C l j = 2 M l=0 lP{L = l} = 2E[L].
We now compute E[L] by induction on the number of cookies.

P{L = l} = 1≤i 1 <i 2 <...<i l ≤M l j=1 (1 -p i j ) j / ∈{i 1 ,...,i l } p j for 0 ≤ l ≤ M.
Decomposing the last sum according to whether i 1 = 1 or i 1 = 1, we obtain

E[L] = M l=1 l 2≤i 2 <...<i l ≤M (1 -p 1 ) l j=2 (1 -p i j ) j / ∈{1,i 2 ,...,i l } p j + M -1 l=0 l 2≤i 1 <...<i l ≤M l j=1 (1 -p i j )p 1 j / ∈{1,i 1 ,...,i l } p j = M -1 l=0 (l -(1 -p 1 )) 2≤i 1 <...<i l ≤M l j=1 (1 -p i j ) j / ∈{1,i 1 ,...,i l } p j = M -1 l=0 (l -(1 -p 1 ))P{ L = l} = E[ L] + 1 -p 1 ,
where L = ♯{2 ≤ i ≤ M, B i = 0}. Finally, we conclude by induction that

E[L] = M i=1 (1 -p i ).
We now study the law of the limiting distribution Z ∞ of the Markov chain Z. This is done via the study of its probability generating function (p.g.f.)

G(s) = E s Z∞ for s ∈ [0, 1].
Lemma 3.4. The p.g.f. G of Z ∞ is the unique p.g.f. solution of the following equation

1 -G 1 2 -s = a(s)(1 -G(s)) + b(s) for all s ∈ [0, 1], (17) with a
(s) = 1 (2 -s) M -1 E [s A M -1 ]
, and

b(s) = 1 - 1 (2 -s) M -1 E [s A M -1 ] + M -2 k=0 G (k) (0) E s A k (2 -s) M -1 E [s A M -1 ] - 1 (2 -s) k .
Proof. The law of Z ∞ is a stationary distribution for the Markov chain Z, therefore

G(s) = E E Z∞ s Z 1 = ∞ k=0 P{Z ∞ = k}E k s Z 1 = M -2 k=0 P{Z ∞ = k}E k s Z 1 + ∞ k=M -1 P{Z ∞ = k}E k s Z 1 .
By definition of Z, for 0 ≤ k ≤ M -2, Z 1 under P k has the same law as A k under P. For k ≥ M -1, Z 1 under P k has the same law as A M -1 + ξ 1 + . . . + ξ k-M +1 where (ξ i ) i≥1 is a sequence of i.i.d. random variables independent of A M -1 and with geometric distribution with parameter 1 2 . Thus,

G(s) = M -2 k=0 P{Z ∞ = k}E s A k + ∞ k=M -1 P{Z ∞ = k}E s A M -1 +ξ 1 +...+ξ k+1-M = M -2 k=0 P{Z ∞ = k}E s A k + E s A M -1 E [s ξ ] M -1 ∞ k=M -1 P{Z ∞ = k}E s ξ k = M -2 k=0 P{Z ∞ = k} E s A k -E s A M -1 E s ξ k+1-M + E s A M -1 E [s ξ ] M -1 G E s ξ . Since E s ξ = 1 2-s , and P{Z ∞ = k} = G k (0), we get G(s) = M -2 k=0 G k (0) E s A k -E s A M -1 (2 -s) M -1-k +E s A M -1 (2-s) M -1 G 1 2 -s ,
from which we deduce that G solves (17). Furthermore, the uniqueness of the solution of this equation amongst the class of probability generating function is a direct consequence of the uniqueness of the stationary law for the irreducible Markov chain Z.

Given two functions f and g, we use the classical notation f (x) = O(g(x)) in the neighbourhood of zero if |f (x)| ≤ C|g(x)| for some constant C and all |x| small enough.

Lemma 3.5. The functions a and b of Lemma 3.4 are analytic on (0, 2). In particular, they admit a Taylor expansion of any order near point 1 and, as x goes to 0:

a(1 -x) = 1 -αx + O(x 2 ), b(1 -x) = O(x).
Proof. Recall the definitions of the random variables A k given in Section 2. Since a geometric random variable with parameter 1 2 admits exponential moments of order strictly smaller than 2, it follows that the p.g.f. s → E[s A k ] are strictly positive and analytic on (0, 2). From the explicit form of the functions a and b given in the previous lemma, we conclude that these two functions are indeed analytic on (0, 2). A Taylor expansion of a near 1 gives

a(1 -x) = 1 -(M -1 -E[A M -1 ]) x + O(x 2 ) = 1 -αx + O(x 2 ), (18)
where we used Lemma 3.3 for the last equality. Since G is a p.g.f. we have G(1) = 1 which, in view of (17), yields b(1) = 0 and therefore b(1

-x) = O(x).
The following proposition relies on a careful study of equation ( 17) and is the key to the proof of Theorem 1.1.

Proposition 3.6. Recall that α = M i=1 (2p i -1) -1 > 0.
The p.g.f. G of Z ∞ is such that, as x goes to 0:

• if 0 < α < 1, then 1 -G(1 -x) ∼ c 1 x α , for some constant c 1 > 0.
In particular E[Z ∞ ] = +∞.

• if α = 1, then 1 -G(1 -x) ∼ c 2 x| ln x|, for some constant c 2 > 0.
In particular E[Z ∞ ] = +∞.

• if α > 1, then 1 -G(1 -x) = c 3 x + O(x 2∧α
) for some constant c 3 > 0.

In particular E[Z ∞ ] < +∞.

Proof. Since G is a p.g.f, it is completely monotonic and we just need to prove the proposition along the sequence x = 1 n with n ∈ N * . Making use of Lemma 3.4 with s = 1 -1 n , we get, for all n ≥ 1

1 -G 1 - 1 n + 1 = a 1 - 1 n 1 -G 1 - 1 n + b 1 - 1 n .
Let us define the sequence (u n ) n≥1 by

u 1 = 1 -G(0) = 1 -P(Z ∞ = 0) > 0, u n = 1-G(1-1/n) n-1 i=1 a(1-1/i) for n ≥ 2. ( 19
)
Hence, (u n ) is a sequence of positive numbers and satisfies the equation

u n+1 = u n + b(1 -1/n) n i=1 a(1 -1/i) , hence u n = u 1 + n-1 j=1 b(1 -1/j) j i=1 a(1 -1/i)
.

This equality may be rewritten

1 -G 1 - 1 n = n-1 i=1 a 1 - 1 i 1 -G(0) + n-1 j=1 b(1 -1/j) j i=1 a(1 -1/i) . ( 20
)
Using Lemma 3.5, we easily obtain

n i=1 a 1 - 1 i = c 4 n α 1 + O 1 n , with c 4 > 0. ( 21
)
Lemma 3.5 also states that, when b is not identically 0 then there exists a unique k ∈ {1, 2, . . .} such that

b(1 -x) = D k x k + O(x k+1 ), with D k = 0.
If b is identically 0, we use the convention k = +∞. In particular, when k is finite, using (21) we deduce that

b(1 -1/n) n i=1 a(1 -1/i) = D k c -1 4 n α-k + O(n α-k-1 ). (22)
Let us now suppose that k = 1. Combining (20), ( 21) and ( 22) we find that 1 -G(1 -1 n ) converges towards D 1 α = 0 as n tends to infinity but this cannot happen because G is continuous at 1 -with G(1) = 1. Thus, we have shown that in fact k ≥ 2.

We now consider the three cases α > 1, α = 1, α < 1 separately.

α > 1

We have three sub-cases: either α > k -1, or α < k -1, or α = k -1 with k ≥ 3.

• α > k -1: Making use of ( 22), we have

n-1 j=1 b(1 -1/j) j i=1 a(1 -1/i) = D k c -1 4 α -k + 1 n α-k+1 + O(1 ∨ n α-k ).
By (20) and (21), we deduce that

1 -G 1 - 1 n = D k (α -k + 1)n k-1 + O 1 n k∧α .
If k was strictly larger that 2, we would have

lim n→∞ n(1 -G(1 -1/n)) = 0
and therefore G ′ (1) = E[Z ∞ ] = 0 which cannot be true because Z is a positive random variable which is not equal to zero almost surely. Thus k must be equal to 2 and

1 -G 1 - 1 n = D 2 (α -1)n + O 1 n 2∧α . (23)
• α < k -1: We prove that this case never happens. Indeed, in view of (22) we find that, for any ε

∈ (0, k -1 -α) b(1 -1/n) n i=1 a(1 -1/i) = O 1 n 1+ε (24)
(this result also trivially holds when k = ∞), thus

∞ j=1 b(1 -1/j) j i=1 a(1 -1/i) < ∞.
Combining this with (20) and ( 21) we see that

1 -G 1 - 1 n = O 1 n α .
Since α > 1, just as in the previous case, this implies that E[Z ∞ ] = 0 which is absurd.

• α = k -1 and k ≥ 3: Again, we prove that this case is empty. Using (22), we now get

b(1 -1/n) n i=1 a(1 -1/i) ∼ D k c -1 4 n .
And, by (20) and (21), we conclude that

1 -G 1 - 1 n ∼ D k ln n n k-1 .
Since k ≥ 3, we obtain E[Z ∞ ] = 0 which is unacceptable.

Thus, we have completed the proof of the proposition when α > 1 and we proved by the way that k must be equal to 2.

α = 1

We first prove, just as in the previous cases, that k = 2. Let us suppose that k ≥ 3. In view of Lemma 3.5, for any l ≥ 3, we can write the Taylor expansion of b of order l near 1 in the form b

(1 -x) = D 3 x 3 + . . . + D l x l + O(x l+1 ) (25)
where D i ∈ R for i ∈ {3, 4, . . . , l}. Similarly,

a(1 -x) = 1 -x + a 2 x 2 + . . . + a l x l + O(x l+1 ),
from which we deduce that, as n goes to infinity

n i=1 a 1 - 1 i = a ′ 1 n + a ′ 2 n 2 + . . . + a ′ l n l + O 1 n l+1 with a ′ 1 > 0. (26)
From ( 25) and (26) we also deduce that

b(1 -1/n) n i=1 a(1 -1/i) = d ′ 2 n 2 + . . . + d ′ l-1 n l-1 + O 1 n l . Thus, n-1 j=1 b(1 -1/j) j i=1 a(1 -1/i) = g 0 + g 1 n + g 2 n 2 + . . . + g l-2 n l-2 + O 1 n l-1 . (27)
Therefore, in view of (20), ( 26) and ( 27), we get

1 -G 1 - 1 n = n-1 i=1 a 1 - 1 i 1 -G(0) + n-1 j=1 b(1 -1/j) j i=1 a(1 -1/i) = λ 1 n + λ 2 n 2 + . . . + λ l-1 n l-1 + O 1 n l .
Comparing with the Taylor expansion of the p.g.f. G, we conclude that E(Z l-1 ∞ ) < ∞ for all l which contradicts Proposition 3.1. Thus, k = 2 and (22) yields

b(1 -1/n) n i=1 a(1 -1/i) ∼ D 2 c -1 4 n with D 2 = 0.
And, by (20) and (21), we conclude that

1 -G 1 - 1 n ∼ D 2 ln n n ,
and therefore

E[Z ∞ ] = +∞. ( 28
) α < 1 Since k ≥ 2, the relation (22) yields ∞ j=1 b(1 -1/j) j i=1 a(1 -1/i) < ∞
(of course, this is trivially true when k = ∞). Thus, the sequence (u n ) defined by (19) converges to a constant c 5 ≥ 0. Suppose first that c 5 = 0. In this case, k cannot be infinite (because when k = ∞, the sequence (u n ) is constant and then c 5 = u 1 > 0). From (22) we deduce that

u n = - ∞ j=n b(1 -1/j) j i=1 a(1 -1/i) ∼ D k (k -α -1)c 4 n k-α-1 ,
therefore, with the help of (21) we get that

1 -G 1 - 1 n = u n n-1 i=1 a 1 - 1 i ∼ D k (k -α -1)n k-1 .
Since k ≥ 2, this implies that n(1 -G(1 -1/n)) converges to a finite constant and so E[Z ∞ ] < ∞. We already notice that this implies a strict positive speed for the cookie random walk in the associated cookie environment (M, p). But (by possibly extending the value of M) we can always construct a cookie environment (M, q) such that p ≤ q and α(q) = 1. In view of (28), the associated cookie random walk has a zero speed and this contradicts a monotonicity result of Zerner (c.f. Theorem 17 of [START_REF] Zerner | Multi-excited random walks on integers[END_REF]). Therefore c 5 cannot be 0 and by ( 19) and (21), we get that

1 -G 1 - 1 n = u n n-1 i=1 a 1 - 1 i ∼ c 5 c 4 n α .
As we already noticed, Theorem 1.1 is now a direct consequence of the last proposition and Corollary 2.3. We also proved that, when the speed is strictly positive, its value is given by the formula

v = α -1 α -1 + 2D 2
where 2D 2 = b ′′ (1) > 0.

Remark 3.7. In the transient case and when the limiting speed is zero, Proposition 3.6 gives with the help of a classical Abelian/Tauberian Theorem the asymptotic of the distribution tail of Z ∞ i.e. the distribution tail of the total number of jumps from 0 to -1:

P {Z ∞ > n} ∼ n→∞ c 6 n α if 0 < α < 1, c 7 ln n n if α = 1. ( 29
)
The functional equation given in Lemma 3.4 for the p.g.f. of Z ∞ also gives a similar equation for the total number of returns R to the origin for the cookie random walk. Indeed, recall that U n 0 (resp. U n 1 ) stands for the respective total number of jumps from 0 to -1 (resp. from 1 to 0) before reaching level n. Thus, the total number of returns to the origin before reaching level n is U n 0 + U n 1 which, under P has the same distribution as Z n + Z n-1 under P. Therefore, we can express the p.g.f. H of the random variable R in term of G:

H(s) = E s Z∞ E Z∞ s Z∞ = 1 a(s) G s 2 -s + M -2 k=0 G (k) (0)s k E s A k - 1 a(s)(2 -s) k .
In particular, Proposition 3.6 also holds for H and the tail distribution of the total number of returns to the origin when α ≤ 1 has the same form as in (29).

Acknowledgments. The authors would like to thank Yueyun Hu for all his precious advices.

Remark 3.8. In the particular case M = 2 (there are at most 2 cookies per site), the only unknown in the definition of the function b is G(0). Since we know that b ′ (1) = 0 ( c.f. the beginning of the proof of Proposition 3.6) we can therefore explicitly calculate G(0), that is the probability that the cookie random walk never jumps from 0 to 1 which is also the probability that the cookie random walk never hits -1. According to the previous remark, we can also calculate the probability that the cookie random walk never returns to 0. Hence, we recover Theorem 18 of [START_REF] Zerner | Multi-excited random walks on integers[END_REF] in the case of a deterministic cookie environment.

Continuity of the speed and differentiability at the critical point

The aim of this section is to prove Theorem 1.2. Recall that

where b ′′ (1) stands for the second derivative at point 1 of the function b defined in Lemma 3.4:

Furthermore, we also proved in Proposition 3.6 that, when α(M, p) = 1, then b ′′ (1) is strictly positive. Hence, in order to prove Theorem 1.2, we just need to show that b ′′ (1) = b ′′ (M,p) (1) is a continuous function of p in Ω u M . It is also clear from the definition of the random variables A k that the functions p → E (M,p) s A k (i) (1) (i.e. the i th derivative at point 1) are continuous in p in Ω u M for all k ≥ 0 and all i ≥ 0 (it is a rational function in p 1 , . . . , p M ). Therefore, it simply remains to prove that, for any k ≥ 0, the function

The following lemma is based on the monotonicity of the hitting times of a cookie random walk with respect to the environment. Lemma 4.1. Let (M, p) be a cookie environment such that α(M, p) > 0. Then there exist ε > 0 and f :

Proof. Let us fix (M, p) with α(M, p) > 0. For ε > 0, define the vector pε = (p ε 1 , . . . , p ε M ) by p ε i = max( 1 2 , p i -ε). We can choose ε > 0 such that α(M, pε ) > 0. Then, for all q ∈ B(p, ε), we have pε ≤ q (30) (where ≤ denotes the canonical partial order on R M ). Let now pick q ∈ B(p, ε), j ∈ N and n ∈ N. Recall that U ∞ 0 denotes the total number of jump of the cookie random walk from 0 to -1 and P (M,q) {Z ∞ = j} = P (M,q) {U ∞ 0 = j} = P (M,q) {X jumps j times from 0 to -1}, and P (M,q) {Z n = j} = P (M,q) {U n 0 = j} = P (M,q) {X jumps j times from 0 to -1 before reaching n}.

Hence

where A is the event "X visits -1 at least once after reaching level n". Recall the notation ω = ω(i, x) i≥1,x∈Z for a general cookie environment given in the introduction. Let now ω X,n denote the (random) cookie-environment obtained when the cookie random walk X hits level n for the first time and shifted by n, i.e. for all x ∈ Z and i ≥ 1, if the initial cookie environment is ω, then

With this notation we have

Besides, X has not eaten any cookie at the sites x ≥ n before time T n . Thus, the environment ω X,n satisfies P (M,q) -almost surely ω X,n (i, x) = q i , for all x ≥ 0 and i ≥ 1 (with the convention

Hence, in view of (30), the random cookie environment ω X,n is P (M,q) -almost surely larger (for the canonical partial order) than the deterministic environment ω pε defined by

2 , for all x < 0 and i ≥ 1, ω pε (i, x) = p ε i , for all x ≥ 0 and i ≥ 1 (with the convention

Thus, Lemma 15 of [START_REF] Zerner | Multi-excited random walks on integers[END_REF] yields P ω X,n {X visits -(n + 1) at least once} ≤ P ω pε {X visits -(n + 1) at least once} P (M,q) -a.s.

In view of (31) we deduce that

where f (n) = P ω pε {X visits -(n + 1) at least once} does not depend of q. It remains to prove that f (n) tends to 0 as n goes to infinity. Let us first notice that

since these probabilities depend only on the environments on the half line [0, +∞). Recall also that the cookie random walk in the environment (M, pε ) is transient (we have chosen ε such that α(M, pε ) > 0), thus

which implies P ω pε {X n = 0 infinitely often} < 1, and a 0 -1 law (c.f. Proposition 5 of [START_REF] Zerner | Multi-excited random walks on integers[END_REF]) yields

Recall that the transition probabilities of the Markov chain Z are given by the law of the random variables A k : P (M,p) {Z n+1 = j | Z n = i} = P (M,p) {A i = j} .

It is therefore clear that for each fixed n and each k, the function p → P (M,p) {Z n = k} is continuous in p in Ω u M . In view of the previous lemma, we conclude that for each k the function p → P (M,p) {Z ∞ = k} is also continuous in p in Ω u M and this completes the proof of Theorem 1.2.


Plazolles Bastien
email: bastien.plazolles@rtech.fr

Vincent Rivola
email: vincent.rivola@rtech.fr

Baz El

Didier

MONTE-CARLO ANALYSIS OF OBJECT REENTRY IN EARTH'S ATMOSPHERE BASED ON TAGUCHI METHOD

The risk of space debris is now perceived as primordial by governments and international space agencies. Since the last decade, international space agencies have developed tools simulate the re-entry of satellites and orbital stations in order to assess casualty risk on the ground. Nevertheless, all current tools provide deterministic solutions, though models include various parameters that are not well known. Therefore, the provided results are strongly dependent on the assumptions made. One solution to obtain relevant and exploitable results would be to include uncertainties around those parameters in order to perform Monte-Carlo analysis. But such a study is very time consuming due to the large parameter space to explore (that necessitate hundreds of thousands simulations). To reduce the parameter search space, we present an application of the Taguchi Method, to model spacecraft debris re-entry in Earth's atmosphere. The Taguchi Method is a statistical analysis method that permits one to determine the parameters uncertainty that have the biggest impact on the results of the numerical simulation. We show how to use this method so as to restrain the quantity of parameters to consider for a Monte-Carlo analysis. Finally, we present the new object-oriented re-entry tool Calima developed by R.Tech. This new tool features three degree of freedom model, featuring also the perturbation of initial parameters of a numerical simulation in order to perform automatized Monte-Carlo Analysis. Calima is accelerated via Graphics Processing Unit (GPU) devices which have many cores architecture and that consume less energy than classical CPUs.

INTRODUCTION

During several decades the problem of space debris was not of primary importance for international space agencies, sending more and more satellites in orbit around the Earth, and either letting them in orbit at the end of their life, or causing their atmospheric re-entry. The first approach presents the major inconvenient that every inactive object must be monitored in order to avoid collision between those debris and active spacecraft or satellite. This strategy is finally expensive and not foolproof with for example the collision in 2009, between IRIDIUM satellite and an inactive Russian satellite, or in 2013 the collision with the Ecuadorian micro-satellite Pegasus and a piece of a Russian rocket tank, one month after launch. The main inconvenient of the second approach is that if the reentry is not well mastered, non-negligible sized objects can reach the ground and cause human casualty and damage to properties. But recently things have changed, with for example the adoption in France of the Law on Space Operations, in 2008, that imposed to French industrial the mastery of technical risks linked with space activities and among others the management of satellites end of life. Moreover, during the 6 th European conference on space debris, the European Space Agency (ESA) highlighted his will to reduce the amount of space debris orbiting around Earth with the "Clean Space" program. In this new international context, the modelling of atmospheric re-entry of spacecraft became of major interest and several space agencies and private companies have developed tools to predict the atmospheric re-entry of space debris in order to assess their casualty area. In this paper, we give an overview in Section 2 of prediction tools used by space agency and present their limitations. In Section 3 we introduce the Taguchi method, a statistical analysis method. In section 4, we present Calima a tool developed by R.Tech in collaboration with LAAS-CNRS. Section 5 deals with the use of Calima and the Taguchi method in order to determine the probability density function of the impact point of an atmospheric reentering object.

RELATED WORK

Several space agencies have developed tools to determine atmospheric re-entry of satellites or spacecraft, to assess their total casualty area, the state of the debris reaching the ground (remaining mass, kinetic Object-oriented Codes consider the satellite as an assembly of simple geometric shapes (box, sphere, cylinder or plate) contained in a parent object. Fig. 1 shows how the atmospheric re-entry of spacecraft or satellites is typically modelled in an object-oriented code. The satellite performs an atmospheric re-entry until it reaches a critical altitude (typically between 72 and 84 km), where it is decomposed into its elementary parts. Here fragmentation is prescribed; the break-up altitude is empirical. Then each component of the satellite is followed independently until it is demised or reaches the ground [START_REF] Lips | A comparison of commonly used re-entry analysis tools[END_REF]. Among the seven published tools, six are object oriented: DAS (NASA) [START_REF]Debris Assessment Software User's Guide[END_REF], ORSAT (NASA) [START_REF] Dobarco-Otero | The Object Reentry Survival Analysis Tool (ORSAT)-Version 6.0 and its Application to Spacecraft Entry[END_REF], ORSAT-J (JAXA) [START_REF] Kato | Toward the world common re-entry safety assessment procedure[END_REF], DRAMA-SESAM (ESA) [START_REF] Martin | A Debris Risk Assessment Tool Supporting Mitigation Guidelines[END_REF], DEBRISK (CNES) [START_REF] Omaly | DEBRISK, a Tool for Re-Entry Risk Analysis[END_REF] and DRAPS (CHINA) [15]. Spacecraft-oriented tools seek to model the spacecraft as close as possible to the real one. In such a tool, there is no prescribed break-up altitude, but fragmentation and ablation of the components of the spacecraft depend on the aerothermodynamics constraints during the atmospheric re-entry of the structure. In case of fragmentation, each fragment is followed independently until it demises or reaches the ground [START_REF] Lips | A comparison of commonly used re-entry analysis tools[END_REF]. To the best of our knowledge there is only one published spacecraftoriented tool: SCARAB [START_REF] Koppenwallner | SCARAB-a Multi Disciplinary Code for Destruction Analysis of Spacecraft During Reentry[END_REF] developed for ESA. The main caveat of all these tools is that they are deterministic. However atmospheric re-entry is a very chaotic system. A slight modification of any initial condition can dramatically change the result of a simulation. Moreover, the deterministic approach does not allow considering the uncertainties on simulations due to the simplification of the aero and thermodynamic models used, the uncertainties of material properties at very high temperature or the uncertainties on the initial conditions of the re-entry. It is still possible to perform Monte-Carlo analysis with existing tools but due to the amount of parameters to consider this could be very time consuming. This raises two questions: is it possible to reduce the amount of parameters to perturb and can we reduce computing time without using a supercomputer or a cluster.

TAGUCHI METHOD

The Taguchi method, introduced first by Genichi Taguchi in 1987 [START_REF] Taguchi | System of the Experimental Design[END_REF], is a statistical analysis method alternative to the traditional Monte-Carlo analysis to process probabilistic phenomenon. Unlike Monte-Carlo analysis which consists in performing a computation for each combination of input parameters, randomly selected on a given space parameter, and which therefore needs a huge amount of simulation to properly map the input space parameter [START_REF] Taguchi | Taguchi's Quality Engineering Handbook[END_REF], the Taguchi method permits one to perform only few computations (depending on the amount of input parameters: 27 computations for less than 13 input parameters, 36 computations for less than 23 parameters,…), by only taking "representatives" values of the input parameters. The counterpart is that unlike the Monte-Carlo analysis, which permits one to determine which set of input parameter drive to a given result, the Taguchi method only permits to determine which parameters are the most influential on the variation of a result and which have less influence, or none. So the Taguchi method is a good tool to reduce the amount of parameters to take into account for a Monte-Carlo analysis.

Taguchi Method principle

The Taguchi Method is based on the "Design of Experiment (DE) concept" developed by Ronald A. Fisher in the 1920s [START_REF] Taguchi | Taguchi's Quality Engineering Handbook[END_REF]. The aim of the design of experiment is to find the relationship between various input parameters of a model and how this model responds to them. The Taguchi method extends the design of experiment by considering uncertainties around input parameters. To do so, the Taguchi method generates combination of "representative" values of the uncertain input parameters thanks to Orthogonal Arrays (OA). The OA are matrices, constructed by combining Latin Square [3], and are defined by three parameters:

-The number of levels, i.e. the set of values that can be taken by elements of the matrix. -The factor number: the number of columns of the table. In the case of the Taguchi Method, each column matches an input parameter. -The observation number: the number of rows of the table. For the Taguchi Method each row corresponds to an experiment (i.e. a combination of input parameters). Taguchi has introduced many standard orthogonal arrays [START_REF] Taguchi | System of the Experimental Design[END_REF], as L 9 presented in Tab. 1 that is also called L 9_3 , with 3 levels, a factor number equal to 4 and 9 observations. Thus, L 9 permits one to define 9 experiments, each being a unique combination of the levels of the four

Experiments A B C D 1 1 1 1 1 2 1 2 2 2 3 1 3 3 3 4 2 1 2 3 5 2 2 3 1 6 2 3 1 2 7 3 1 3 2 8 3 2 1 3 9 3 3 2 1 Table 1: Standard Orthogonal Array L 9
parameters. The key point of the standard Taguchi Method is to define three "representative" values for a given parameter that will correspond to the three levels in the OA. There exist variations of this method considering combination of input parameters with three "representative" values and others with two [START_REF] Sokal | Biometry, The Principles and Practices of Statistics in Biological Research[END_REF].

Basically, Taguchi Method considers a normal uncertainty around each considered input parameter [START_REF] Mistree | Simulation Reduction using Taguchi Method[END_REF] defined via a mean value μ and a variance σ. It is commonly accepted to consider the three following values as "representative": μ, and μ ± Δ, Δ being a tolerance around the mean value and taken equal to  5 . 1 (see [START_REF] Taguchi | Taguchi's Quality Engineering Handbook[END_REF], [3], [START_REF] Phadke | Quality Engineering Using Robust Design[END_REF] and [START_REF] Suh | The Principles of Design[END_REF]). It is common practice to match levels and representative values as follows: level => µ -Δ, level 2=> µ, level 3=> µ + Δ.

Analysis of Variance

Since the Taguchi Method only samples few values of the whole input space parameters, the result analysis must include an analysis of the confidence that can be placed in the results. That is the aim of a standard technique called Analysis of Variance (ANOVA). This technique permits one to determine the variability of the results, how this variability is influenced by values of the input parameters and what level of confidence can be assigned to the previous results. ANOVA necessitates to compute many quantities and organizes them in a standard tabular format. We will not give the details in this paper all those parameters, but here is the list for informative purpose: the variance V, the sum of square SS, the pure sum of square SS', degrees of freedom (D.o.F.) f, the error e, the variance ratio F, the percent contribution P, the correction factor CF, the number of experiments n and the total degree of freedom f T (for more information interested readers should refer to chapter 5.3.4 of [START_REF] Roy | A Primer on the Taguchi Method[END_REF]). In order to determine which parameter has the biggest influence on a result and the level of confidence on these results the following parameters have to be considered: the Percent Contribution and the Variance Ratio. Percent Contribution P permits one to determine how much the variation of a given result is caused by the effects of an input parameter. Thus, the input parameter with the highest Percent Contribution for a given result, is the input parameter that has the biggest influence on this result. The percent contribution is obtained via Eq. 1

t SS SS P / ' * 100  (1)
where SS t is the total square sum. The Variance Ratio F, permits one to determine the confidence level if an input parameter really contributes to the variation of a given result. The F value computed is then compared with the values from the F-tables for a given level of significance. If the computed F value is less than the value from the F-table, then the parameter does not contribute to the variation of the result [START_REF] Roy | A Primer on the Taguchi Method[END_REF]. The variance ratio is obtained via Eq. 2.

e V V F /  (2)
where V e denotes the variance of the error.

Monte Carlo method

The Monte Carlo method consists of using a random number generator to simulate a large number of combinations of parameters within tolerances. The results of every generated combination of input parameters is computed, after which the mean and the variance of the results are computed. The Monte Carlo method also permits one to explore the space of results in order to find the combination of input parameters that provides a desired result [3].

For obtaining accurate results that are not perturbed with statistical bias due to the generation of the random numbers, Monte Carlo method necessitates to perform a large number of simulation (thousands or hundreds of thousands), all the more important that there is input parameters, and such a method could be very time consuming, especially if each simulation requires a large amount of computing time.

CALIMA

Due to their deterministic approach, and because of the uncertainties in the models they used, all the existing tools use very conservative models that can be very limiting for industrials and threaten the durability af space activities. We decided to opt for a complete different approach by developing a brand new tool using less conservative models but able to perform statistical analysis : Calima. Calima is a three Degree of Freedom (D.o.F.) objectoriented atmospheric re-entry simulation tool, developed by R.Tech, in collaboration with the C.D.A. team of LAAS-CNRS. So far, Calima is a basic tool with imposed aero-coefficients and no thermal model (no ablation). The aim of Calima is to perform statistical analysis of atmospheric re-entry debris in order to determine the casualty area. Currently, there are 11 input parameters that can be perturbed in Calima: the initial mass, altitude, latitude, longitude, azimuth, speed, angle of Once the uncertainties on the parameters, and the perturbation method are filled in Calima, the tool will automatically, ensures the coverage of the input parameters space. Also, in order to prepare a Monte Carlo analysis Calima permits one to perform an automated Taguchi analysis, which will indicate to the user which are the input parameters most infuencial. Furthermore, Calima permits one to reduce the required computing time of Monte-Carlo analysis by taking advantage of Graphics Processing Unit (GPU) to parallelize computations. Indeed, the GPUs are highly parallel, multithreaded, many-core architectures. They are better known for image processing. Nevertheless, NVIDIA introduced in 2006 CUDA (Compute Unified Device Architecture), a technology that enables users to take benefit of GPU cards to address parallel applications, multi-parametric problems and complex optimisation problems [START_REF] Ilg | Projectile Monte-Carlo Trajectory Analysis Using a Graphics Processing Units[END_REF], [START_REF] Rogers | Robust Parafoil Terminal Guidance Using Massively Parallel Processing[END_REF], | [19] [20]. As shown in Fig. 3, a parallel code on GPU (hereafter named the device), is interleaved with a serial code executed on the CPU (hereafter named the host).

The parallel threads are grouped into blocks which are organized in a grid. The grid is launched via a single CUDA program, the so-called kernel. The GPU implementation of Calima is performed using CUDA 6.0. The implementation of Calima is divided in three synchronous steps, described Fig. 2. First, the CPU part of the code loads the initial conditions, and uncertainties on the unified memory between the CPU and the GPU.

Then, each parallel thread on the GPU performs a complete re-entry simulation from the initial point until the object reaches the ground, using a unique set of initial conditions derived from the initial conditions and uncertainties provided by the user. Finally, when threads have finished, the host recovers the results from the unified memory.

COMPUTATIONAL TESTS

In this chapter we display and analyse the results relative to the application of the Taguchi method. The goal is the determination of the impact area of an object and the range of kinetic energy at the impact, using the Calima code. As the current paper focusses on the methodology the presented results are not workable. Perturbed parameters and physical models implemented in Calima will be extended in future versions.

Computational test

We are considering a sphere with a diameter of 1.0 m and we define the initial state of this sphere with the parameters provided in Tab. 2. We assume Gaussian uncertainties around each of these parameters, with the mean values and standard deviation defined in Tab. 2.

Figure 1 :

 1 Figure 1: Sketch of spacecraft atmospheric re-entry energy, etc. …) and coordinates of impact point of each debris. To date, there are seven published atmospheric re-entry tools that can be divided in two categories: Object-oriented Code and the Spacecraft Oriented code.Object-oriented Codes consider the satellite as an assembly of simple geometric shapes (box, sphere, cylinder or plate) contained in a parent object. Fig.1shows how the atmospheric re-entry of spacecraft or satellites is typically modelled in an object-oriented code. The satellite performs an atmospheric re-entry until it reaches a critical altitude (typically between 72 and 84 km), where it is decomposed into its elementary parts. Here fragmentation is prescribed; the break-up altitude is empirical. Then each component of the satellite is followed independently until it is demised or reaches the ground[START_REF] Lips | A comparison of commonly used re-entry analysis tools[END_REF]. Among the seven published tools, six are object oriented: DAS (NASA)[START_REF]Debris Assessment Software User's Guide[END_REF], ORSAT (NASA)[START_REF] Dobarco-Otero | The Object Reentry Survival Analysis Tool (ORSAT)-Version 6.0 and its Application to Spacecraft Entry[END_REF], ORSAT-J (JAXA)[START_REF] Kato | Toward the world common re-entry safety assessment procedure[END_REF], DRAMA-SESAM (ESA)[START_REF] Martin | A Debris Risk Assessment Tool Supporting Mitigation Guidelines[END_REF], DEBRISK (CNES)[START_REF] Omaly | DEBRISK, a Tool for Re-Entry Risk Analysis[END_REF] and DRAPS (CHINA)[15]. Spacecraft-oriented tools seek to model the spacecraft as close as possible to the real one. In such a tool, there is no prescribed break-up altitude, but fragmentation and ablation of the components of the spacecraft depend on the aerothermodynamics constraints during the atmospheric re-entry of the structure. In case of fragmentation, each fragment is followed independently until it demises or reaches the ground[START_REF] Lips | A comparison of commonly used re-entry analysis tools[END_REF]. To the best of our knowledge there is only one published spacecraftoriented tool: SCARAB[START_REF] Koppenwallner | SCARAB-a Multi Disciplinary Code for Destruction Analysis of Spacecraft During Reentry[END_REF] developed for ESA. The main caveat of all these tools is that they are deterministic. However atmospheric re-entry is a very chaotic system. A slight modification of any initial condition can dramatically change the result of a simulation. Moreover, the deterministic approach does not allow considering the uncertainties on simulations due to the simplification of the aero and thermodynamic models used, the uncertainties of material properties at very high temperature or the uncertainties on the initial conditions of the re-entry. It is still possible to perform Monte-Carlo analysis with existing tools but due to the amount of parameters to consider this could be very time consuming. This raises two questions: is it possible

Figure 2 :

 2 Figure 2: Steps of the GPU implementation of Calima

Table 2 :

 2 Initial conditions and uncertaintiesWe use the Orthogonal Array L 27_3 , and we assign the

	Parameter	Mean µ	Standard deviation σ
	Mass (kg)	247.224	4.12
	Altitude (km)	120.0	0.25
	Latitude (°)	0.0	0.05
	Longitude (°)	0.0	0.1
	Azimuth (°)	45.0	0.15
	Speed (m.s-1) 7272.582	48.51
	Flight Path Angle	-2.612	0.01

mass to the first column of the OA, the altitude to the second and so on. The three representatives values for each of this input parameter are defined as explained in In this manner we define 27 experiments, each being a combination of the "representative" values of the seven input parameters. These combinations are automatically performed by Calima.

Analysis of Variance

Once the 27 cases are computed by Calima, the software performs the Analysis of Variance of the results. Tab. presents the ANOVA table for the variation of the kinetic energy. We can see immediately that the parameter having the biggest impact on the variation of the impact energy is the mass, with a 92.38% contribution. We note that the variance ratio is equal to 8.43*10 4 , meaning that the level of confidence is higher than 99.5%. We can also note that the next largest contributor is the initial speed with only a 7% contribution, with a confidence level of 99.5% too. The other point to notice is that the initial latitude, longitude and F.P.A. have negligible influence on the variation of the impact energy, with contributions around 10 -4 %. While performing the same analysis for the longitude and the latitude (Tab. 4 and 5) we note that the parameters that have the most impact on their variation are the initial latitude and longitude. In the same manner we note that the initial mass, azimuth and flight path angle have almost no impact on the variation ofimpact coordinates with contributions around 10 -6 % for the initials mass and the azimuth, and around 10 -10 % for the initial flight path angle. All those results are strongly dependent on the "representative" values selected to perform the Taguchi analysis, i.e. on the estimation of the uncertainties on the initial conditions. Indeed the results presented here are only valid, for the uncertainties we assumed on the initial conditions, see Tab. 2. For example, if instead of considering the standard deviation for the input values of mass and initial speed, we consider a standard deviation of the uncertainty around the mass of 1kg, and a standard deviation for the uncertainty about the initial speed of 149 m.s -1 and the same uncertainties for the other input parameters, then the ANOVA table of the impact energy (Tab. 6) shows that now the input parameter that influence the most the variation of the impact energy is the initial speed, while the mass has almost the same contribution as the azimuth now. This highlights the importance to properly determine the uncertainties around the input parameters.

Monte-Carlo analysis

Thanks to the Taguchi Method, we have established that, for our case, we could only consider the initial altitude, latitude, longitude and speed to perform a Monte-Carlo analysis in order to determine the impact The nominal simulation performed with Calima, using these parameters, impact the ground at 0.0° latitude and 13.56° longitude, with impact energy of 245.98 kJ.

In order to avoid statistical bias in the distribution of the results, we performed 200 000 simulations, in Calima.

While such an amount of simulations would take around 300 hours of sequential computations using a classic reentry simulation tool such like Debrisk on a classic CPU (Intel Xeon E5640), Calima performed all these computations in 88 seconds on a Tesla K40 GPU accelerator, with 2880 CUDA cores at 0.745GHz and 12 GB memory, four order of magnitude time faster. It is also interesting to note that, on the same CPU, Calima computes sequentially all those computations in 37 000s. Fig. 4 and 5 present the probability distribution of the impact point in latitude and in longitude. As we can see, the small uncertainties provided in the initial conditions lead to a great dispersion of the possible impact point. Thus the 3 sigma area defined is about 10° wide in latitude for 30° wide in longitude, for an ellipse of 50 257 km.

CONCLUSION AND FUTURE WORKS

In this paper we have shown the interest of using the Taguchi Method for simulation of object re-entry in Earth's atmosphere. With few computations, this method permits one to determine which are the parameters that have the biggest impact on the variation of results of a re-entry model. One of the main advantages of this method is that the Taguchi method can be used with already existing tools.

We have also introduced Calima, the new objectoriented atmospheric re-entry simulation tool, in development by R.Tech, in collaboration with LAAS-CNRS. This new tool has the unique feature to intrinsically take into account the probabilistic nature of simulation parameters and initial condition, with automated exploration of the input space parameter, for Monte-Carlo analysis. Calima allows natively to perform a Taguchi analysis. Calima also includes the possibility to dramatically reduce the computing time by taking advantage of GPUs to massively parallelize the computation.

In future work, we will increase the amount of parameters that can be perturbed. We will also improve physical models in Calima, by including the possibility to define shapes of the object, computation of the aerocoefficient and taking into account ablation for examples. In long term evolution, we want to make a transition from object-oriented code to spacecraftoriented code, and it is also envisaged to include the possibility to compute simultaneously several objects and their interactions (wake, etc…).

ACKNOWLEDGEMENT

The authors gratefully acknowledge the support of NVIDIA Corporation with the donation of the Tesla K40 GPU used for this research work and the ANRT which allowed the funding of the PhD of M. Plazolles.