Ali Taheri

published or not. The documents may come

• - - - • • • • - - - - - • • • • - - - - - • • • • • • • • 1. 2. • ⊆ ⊆ ⊆ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ∀ ∀ ∀ ∀ ∀ ∏ ∏ ∏ • ∑ ∑ • • • • ∀ ∀ ∀ • • • • • • • • ∀ ∀ ∀ ∀ • • β β β β β β β β
Tab. 49. The ranking of different conditions that emerge by new changes occurring in a fbs-chain;

μ μ μ μ μ μ μ μ μ μ μ μ μ μ μ μ μ μ μ μ μ μ μ μ μ μ μ λ λ λ λ λ λ λ λ λ λ λ λ θ θ τ τ τ τ β ∑ μ β ∑ λ μ β ∑ θ τ β β β μ λ μ θ τ μ μ  μ λ λ λ λ  μ μ  λ λ λ λ  μ • • • ∑ ∑ ∑ ∑ • • • • • • 𝑢𝐸𝑃 𝑖 = 1 -( |𝐴𝑉| |𝑈𝐼𝑉| ) 𝑢𝐸𝑃 𝑖 = ( |𝐴𝑉| |𝑈𝐼𝑉 -𝐼𝑉| ) 𝑒𝐸𝑃 𝑖 = 1 -𝑢𝐸𝑃 𝑖 ∑ ∑ N(t) = κ 2 { 𝑡 -𝑡 𝑚 |𝑡 -𝑡 𝑚 | . [1 -exp ( 𝑙𝑛(0.04) ∆𝑡 . |𝑡 -𝑡 𝑚 |)] + 1} N(t) = κ (1 + 𝑒 -αt-β ) κ N(t) = 𝑒 αt+β κ κ α β κ κ κ κ κ α α α κ κ β 1. 2. 3. 4. 5. 6. α • • • κ κ κ • • • • μ TUD a = ∑ (TPM i . μ g i n 1 ) ∑ μ g i n 1 μ ∑ η η η ∑ ∑ ∑ • • • • • • • • • • • • 1.     2

Il y a 37 questions dans ce questionnaire.


Vincent Bruneau

Adrien Doradoux

Pierre Fabrie

Convergence of a Vector

come

Introduction

Numerical simulation of unsteady incompressible flows has always been the subject of important investigations. The difficulty arises from the coupling between the velocity and the pressure at each time step due to the incompressibility constraint. Fractional step algorithms are widely used to deal with this problem. Among them, pressure projection methods introduced by Chorin and Temam ( [6], [START_REF] Temam | Sur l'approximation de la solution des équations de navier-stokes par la méthode des pas fractionnaires (ii)[END_REF]) reduce the saddle point problem to two distinct elliptic problems on the velocity and the pressure. Recently, Vector Penalty Projection schemes have been introduced by Angot et al. ( [START_REF] Angot | Vector penalty projection methods for the solution of unsteady incompressible flows[END_REF], [START_REF] Angot | A fast vector penalty projection for incompressible nonhomogeneous or multiphase navier-stokes problems[END_REF]) and avoid many drawbacks of other projection methods. In particular, the pressure scalar does not need to be computed, which do not impose the resolution of a Poisson type equation that introduce boundary conditions on the pressure. In [START_REF] Angot | Vector penalty-projection methods for incompressible fluid flows with open boundary conditions, in Algoritmy[END_REF] the authors obtained a second order convergence rate for pressure and velocity in space and time for a 1 second order backward temporal scheme. The convergence towards the Navier-Stokes problem when the penalty parameter (on the divergence) tends to 0 has been studied in [START_REF] Angot | Analysis for the fast vector penalty-projection solver of incompressible multiphase navier-stokes/brinkman problems[END_REF]. In these previous works the domain is fixed and no space penalization is included.

Let Ω be a simply connected bounded domain of IR d (d = 2 or 3) and T > 0.

With Dirichlet boundary condition and a source term f given, the penalized Navier Stokes problem reads in presence of a body moving at the velocity v s :

∂v ∂t + (v.∇)v -div(2µD(v)) + 1 η χ ω(t) (v -v s ) + ∇p = f on Ω div(v) = 0 on Ω v = 0 on ∂Ω v(0, x) = v 0 (x) on Ω (1)
Where χ ω(t) is the characteristic function of the solid domain ω(t), v s is the velocity of the moving body and D is the strain rate tensor. Hypothesis (H1): We suppose that v s is the restriction to t<T ω(t) of a function ψ such that:

             div(ψ) = 0 on Ω ψ ∈ L ∞ (]0, T [; L ∞ ) ∩ L 2 (]0, T [; H 2 ) ∇ψ ∈ L ∞ (]0, T [; L ∞ ) ∂ψ ∂t ∈ L 2 (]0, T [; L 2 ) (2)
The existence of such function is ensured as the regularity of v s is sufficient and the moving body does not meet ∂Ω.

(H2): We assume that f belongs to the space L 2 (]0, T [; L 2 ). The Vector Penalty Projection Scheme is a fractional step method:

• A predicted velocity ṽn+1 is first computed (v n → ṽn+1 ): considering the pressure gradient at the previous time step t n . At the end of this step, the velocity does not respect the free divergence condition.

• The velocity is then corrected such that div(v n+1 ) is approximately 0 at the end of the time step.

• We finally actualize the pressure gradient ∇p n+1 .

For all n ∈ N such that nδt ≤ T the numerical scheme reads:

ṽn+1 -v n δt + B(v n , ṽn+1 ) -div(2µD(ṽ n+1 )) + 1 η χ ω(t n+1 ) (ṽ n+1 -v s ) + ∇p n = f n+1 (3)
ε δt vn+1 -∇ (div(v n+1 ) + div(ṽ n+1 )) = 0 (4)

∇(p n+1 -p n ) + 1 ε ∇(div(v n+1 )) = 0 (5)
Where v n+1 = ṽn+1 + vn+1 . It is completed by the following initial and boundary conditions on ∂Ω:

ṽn+1 = 0 on ∂Ω, vn+1 .ν = 0 on ∂Ω (6) ṽ0 = v 0 in Ω, v0 = 0 in Ω (7)
where ν is the outward unit normal vector on ∂Ω.

The space Ω is a bounded simply connected domain of IR d . Furthermore, p 0 ∈ L 2 0 the L 2 -space of null average functions as well as div(v n+1 ) since v n+1 .ν = 0. Therefore we can show recursively that p n+1 has a null average on Ω, for all n ∈ N, solving (5) in the space of null average functions. We finally obtain:

ε(p n+1 -p n ) + div(v n+1 ) = 0 ( 8
)
Proposition 0.1 (Existence of the iterates).

We suppose that v 0 ∈ H 1 and p 0 ∈ L 2 0 . Then , for all n ∈ N such that nδt ≤ T , we have: 3) is a linear Stokes problem, therefore we easily show that ṽn+1 exists and lies in H 1 0 . We show now the existence of vn+1 . For a simply connected domain, the norm [START_REF] Foias | Remarques sur les équations de navier-stokes stationnaires et les phénomènes successifs de bifurcation[END_REF]). By energy estimates, we obtain that for any n ∈ N there exists a constant C > 0 such that:

(v n , p n ) ∈ H 1 × L 2 0 Proof. Suppose that v n ∈ H 1 and show v n+1 ∈ H 1 . Equation (
|| • || 2 H 1 is equivalent to the norm ||div(•)|| 2 L 2 + ||curl(•)|| 2 L 2 (see
ε δt ||v n+1 || 2 L 2 + 1 2 ||div(v n+1 )|| 2 L 2 ≤ C
Moreover, vn+1 is a gradient, therefore curl(v n+1 ) = 0. We deduce that vn+1 lies in H 1 . Using (5), we deduce that ∇p n+1 lies in H -1 . As p n+1 has a null average on Ω, we conclude using the Poincaré Wirtinger inequality that p n+1 ∈ L 2 0 .

Main result

In this paper, to each sequence (v k ) k of functions defined on Ω we will associate a sequence of functions (v δt ) δt which are the step functions in time v δt defined by:

v δt (t) = v k if t ∈ [t k , t k+1 [. (9)
For the step functions contructed from the sequences of Proposition 0.1 we have the following convergence result.

Theorem 1.1 (Convergence when ε and δt tend to 0).

Let Ω ∈ IR d be a simply connected bounded domain and

d = 2 or 3. Let µ > 0 and v 0 ∈ H where H = v ∈ L 2 , div(v) = 0, v.ν = 0 on ∂Ω .
We assume that there exists λ > 0 such that ε = λδt and we suppose that hypothesis (H1) and (H2) are verified.

Then, up to a subsequence, (v δt , p δt ) defined by (3)-( 7) and (9) converges towards (v, p), a weak solution of the penalized Navier Stokes problem (1) when ε = λδt tends to 0. Furthermore, v and p satisfy:

v ∈ L ∞ (]0, T [; L 2 ) ∩ L 2 (]0, T [; H 1 ) p ∈ W -1,∞ (]0, T [; L 2 0 ) Moreover, this solution is unique for d = 2.
The proof of this theorem is given in section 4. We complete this last result by giving a new proof of existence of weak solutions to the Navier Stokes equation in presence of a moving body. To do so, we focus on the behavior of the solution when η tends to 0 and denote v η the weak limit obtained in Theorem 1.1. Theorem 1.2 (Convergence when η tends to 0). When η tends to 0, the sequence (v η ) η weakly converges towards a limit v which satisfies

v |∂ω(t) = v s (t)
Furthermore, v |Ω\ω(t) is a weak solution of the Navier Stokes equations on Ω \ ω(t). Moreover, there exists h ∈ W ′ such that:

1 η χ ω(t) (v -v s ) -⇀ h in W ′ ( 10
)
Where W is defined in section 2.1.

The proof of this result is given in Section 5. Before proving these theorems, we first prove the stability of the numerical scheme. In order to perform energy estimates, we need to build a lifting of the velocity. In the following to the function v (resp. v n , ṽn , vn , v δt ) we will associate the lifting functions w (resp. w n , wn , ŵn , w δt ) defined by substraction of the function ψ:

w = v -ψ, w n = v n -ψ, wn = ṽn -ψ, ŵn = vn , w δt = v δt -ψ. (11)
The system (3)-( 7) becomes:

wn+1 -w n δt + B(w n , wn+1 ) + B(ψ n , wn+1 ) -div(2µD( wn+1 )) + 1 η χ ω(t n+1 ) wn+1 + ∇p n = F n+1 -B(w n , ψ n+1 ) (12) ŵn+1 δt - 1 ε ∇(div( wn+1 + ŵn+1 )) = 0 (13) ∇(p n+1 -p n ) + 1 ε ∇(div(w n+1 )) = 0 ( 14
)
w 0 (x) = v 0 -ψ 0 (15)
Where

F n+1 = f n+1 - ψ n+1 -ψ n δt + div(2µD(ψ n+1 )) -B(ψ n , ψ n+1 ). Proposition 1.1 (Stability). For µ > 0, v 0 ∈ L 2 , p 0 ∈ L 2 0 , ∇p 0 ∈ L 2 and f ∈ L 2 (]0, T [; L 2 )
given, if we assume that ψ verifies the hypothesis (H1). Then, for any T > 0 there exists

C > 0 such that i) w δt is bounded in L ∞ (]0, T [; L 2 ). ii) wδt (• + δt) -w δt (•) is bounded in L ∞ (]0, T -δt[; L 2 ) and || wδt (• + δt) -w δt (•)|| L ∞ (]0,T -δt[;L 2 ) ≤ C √ δt iii) ∇ wδt is bounded in L 2 (]0, T [; L 2 ).
The proof of this proposition can be found in section 3.

2 Mathematical recalls

Notations

Let us define the main notations of this paper.

• For p > 0, L p = L p (Ω).

• For p > 0,

L p 0 = v ∈ L p ; Ω v dx = 0 .
• For p ∈ IR, H p = H p (Ω), the classical Sobolev space.

• For p > 0,

H p 0 = v ∈ H p ; v |∂Ω = 0 . • H = v ∈ L 2 ; div(v) = 0 on Ω; γ ν (v) = 0 . • V = v ∈ H 1 0 ; div(v) = 0 on Ω . • W = v ∈ L 2 (]0, T [; V); ∂v ∂t ∈ L 2 (]0, T [; V ′ ); v(T ) = 0 .
Where γ ν is the trace operator.

We recall the definition of the Nikolskii spaces.

Definition 2.1. Let E be a Banach space and f ∈ L 1 (]0, T [; E). For 1 ≤ q < ∞ and 0 < σ < 1, we define the Nikolskii space N σ q (]0, T [; E) by

N σ q = f ∈ L q (]0, T [; E), sup 0<h<T ||τ h f -f || L q (]0,T -h[;E) h σ < ∞
and for f ∈ N σ q (]0, T [; E) we define the associated norm

||f || N σ q (]0,T [;E) = ||f || q L q (]0,T [;E) + sup 0<h<T 1 h σ ||τ h f -f || L q (]0,T -h[;E) q 1 q
Where τ h is the translation operator defined by

τ h f (•) = f (• + h)
We denote ( • , • ) L 2 the usual scalar product on L 2 and < • , • > E ′ ,E the duality bracket. In all the paper the constants are always denoted C.

Mathematical properties

To deal with the nonlinear convective term, we use the bilinear form B introduced by Temam (see [START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier Stokes[END_REF] and [START_REF] Temam | Navier-Stokes Equations. Theory and Numerical Analysis[END_REF]).

For u ∈ H 1 and v ∈ H 1 0 , B(u, v) = (u.∇)v + 1 2 div(u)v (16)
Taking the scalar product of B(u, v) by w ∈ H 1 0 and integrating by part the second term, we obtain the associated trilinear form b:

b(u, v, w) = 1 2 Ω (u.∇)v.w dx - 1 2 Ω (u.∇)w.v dx (17)
The trilinear form b satisfies the antisymmetry property b(u, v, w) = -b(u, w, v) and b(u, v, v) = 0. We now recall the discrete Gronwall Lemma (see [START_REF] Holte | Discrete gronwall lemma and applications[END_REF], [START_REF] Gil | Difference equations in normed spaces: stability abd oscillations[END_REF], [START_REF] Quarteroni | Numerical Models for Differential Problems[END_REF]).

Lemma 2.1 (Discrete Gronwall lemma [START_REF] Holte | Discrete gronwall lemma and applications[END_REF]).

Let (y n ), (f n ) and (g n ) three non-negative sequences such that:

y n ≤ f n + n-1 k=0 g k y k for n ≥ 0
Then,

y n ≤ f n + n-1 k=0 f k g k exp   n-1 j=k+1 g j   for n ≥ 0
In order to prove the convergence of the velocity we will need the following analysis result ( [5], p99). Let X and Y two Banach spaces such that X is embedded in a continuous and dense way into Y , and let T > 0 and p, q satisfy 1 ≤ p, q ≤ +∞. We denote:

E p,q = u ∈ L p (]0, T [, X), du dt ∈ L q (]0, T [, Y ) Proposition 2.1.
Any element u of E p,q (defined almost everywhere) possesses a continuous representation on [0, T ] with values in Y , and the embedding of

E p,q into C 0 ([0, T ], Y ) is continuous. Moreover, for all t 1 , t 2 ∈ [0, T ] we have u(t 2 ) -u(t 1 ) = t2 t1 du dt dt
where it is understood that we have identified u and its continuous representation.

Finally, let us formulate an important compactness theorem, which will be useful to prove the strong convergence of the velocity in section 4.

Theorem 2.1 (Simon [START_REF] Simon | Compact sets in the space l p (0, t; b)[END_REF]). Let B 0 , B 1 and B 2 three Banach spaces with B 0 ⊂ B 1 ⊂ B 2 . We suppose that the embedding of B 0 in B 1 is compact and the embedding of B 1 in B 2 is continuous. Then, for all 1 ≤ q ≤ +∞ and 0 < σ < 1, the embedding

L q (]0, T [; B 0 ) ∩ N σ q (]0, T [; B 2 ) ֒→ L q (]0, T [; B 1 )
is compact.

Stability analysis

In this section we prove the stability of the numerical scheme on the problem ( 12)- [START_REF] Temam | Sur l'approximation de la solution des équations de navier-stokes par la méthode des pas fractionnaires (ii)[END_REF]. To do so, energy estimates are performed on the prediction, the correction and pressure equations respectively given by ( 12), ( 13) and ( 14). Then we will obtain an upper bound on div(w) in L 2 (]0, T [; L 2 ). By definition of w δt (see [START_REF] Holte | Discrete gronwall lemma and applications[END_REF] and ( 11)) Proposition 1.1 will be a direct consequence of the following result

Proposition 3.1 (Stability). For µ > 0, v 0 ∈ L 2 , p 0 ∈ L 2 0 , ∇p 0 ∈ L 2 and f ∈ L 2 (]0, T [; L 2 )
given, if we assume that ψ verifies the hypothesis (H1). Then, for any T > 0 there exists

C > 0 such that i) the sequence (w k ) k is bounded in L 2 (Ω).
ii) the sequence

wk+1 -w k √ δt k is bounded in L 2 (Ω). iii) n-1 k=0 δt||∇ wk+1 || 2 L 2 is uniformly bounded with respect to δt.
Proof of Proposition 3.1.

We obtain the result using several energy estimates as in [START_REF] Angot | Analysis for the fast vector penalty-projection solver of incompressible multiphase navier-stokes/brinkman problems[END_REF] for homogeneous Navier-Stokes flows. In our estimates, an additional term appears due to the penalization term on the moving body. Taking wn+1 as a test function in [START_REF] Quarteroni | Numerical Models for Differential Problems[END_REF], we obtain:

1 δt wn+1 -w n , wn+1 L 2 + 2µ D( wn+1 ), D( wn+1 ) L 2 + 1 η Ω χ ω(t n+1 ) | wn+1 | 2 dx + b(w n , wn+1 , wn+1 ) + b(ψ n , wn+1 , wn+1 ) + ∇p n , wn+1 L 2 = F n+1 , wn+1 L 2 -b(w n , ψ n+1 , wn+1 )
The diffusion term is integrated by parts, the Korn inequality ( [START_REF] Korn | Die eigenschwingungen eines elastichen korpers mit ruhender oberflache[END_REF], [START_REF] Korn | Ubereinige ungleichungen, welche in der theorie der elastichen und elektrischen schwingungen eine rolle spielen[END_REF]) is then used to obtain the lower bound:

||∇ wn+1 || 2 L 2 ≤ 2||D( wn+1 )|| 2 L 2 .
The convective terms b(w n , wn+1 , wn+1 ) and b(ψ n , wn+1 , wn+1 ) vanish by antisymmetry of the trilinear form b. By definition of b given in (17) and using standard norm estimates, there exists a constant C > 0 such that for any µ > 0:

|b(w n , ψ n+1 , wn+1 )| ≤ 1 2µ ||w n || 2 L 2 C||∇ψ n+1 || 2 L ∞ + ||ψ n+1 || 2 L ∞ + µ 4 ||∇ wn+1 || 2 L 2 ≤ C||w n || 2 L 2 + µ 4 ||∇ wn+1 || 2 L 2
We finally use the following equality:

(a -b, a) = 1 2 ||a|| 2 -||b|| 2 + ||a -b|| 2 . ( 18
)
The following estimate is obtained:

1 2δt || wn+1 || 2 L 2 -||w n || 2 L 2 + || wn+1 -w n || 2 L 2 + µ 2 ||∇ wn+1 || 2 L 2 + 1 η Ω χ ω(t n+1 ) | wn+1 | 2 dx + (∇p n , wn+1 ) L 2 ≤ C||F n+1 || 2 L 2 + C||w n || 2 L 2 . (19)
From ( 13) and ( 14), we have:

w n+1 -wn+1 δt + ∇(p n+1 -p n ) = 0. ( 20
)
By taking w n+1 as a test function in (20) and using again (18) we obtain:

1 2δt ||w n+1 || 2 L 2 -|| wn+1 || 2 L 2 + ||w n+1 -wn+1 || 2 L 2 + ∇p n+1 -∇p n , w n+1 L 2 = 0. ( 21
)
Choosing p n+1 as a test function in (8), we have:

ε 2 ||p n+1 || 2 L 2 -||p n || 2 L 2 + ||p n+1 -p n || 2 L 2 -∇p n+1 , w n+1 L 2 = 0. ( 22
)
At last, taking ∇p n+1 as a test function in (20) we obtain:

δt 2 ||∇p n+1 || 2 L 2 -||∇p n || 2 L 2 + ||∇p n+1 -∇p n || 2 L 2 + ∇p n+1 , w n+1 -wn+1 L 2 = 0 (23)
Finally, these four estimates ( 19), ( 21)-( 23) are summed up. The sum of the scalar products reduces to ∇p n+1 -∇p n , w n+1 -wn+1 , which is bounded using Young inequality:

| ∇(p n+1 -p n ), w n+1 -wn+1 L 2 | ≤ δt 2 ||∇p n+1 -∇p n || 2 L 2 + 1 2δt ||w n+1 -wn+1 || 2 L 2
Therefore,

1 2δt ||w n+1 || 2 L 2 -||w n || 2 L 2 + || wn+1 -w n || 2 L 2 + µ 2 ||∇ wn+1 || 2 L 2 + ε 2 ||p n+1 || 2 L 2 -||p n || 2 L 2 + ||p n+1 -p n || 2 L 2 + δt 2 ||∇p n+1 || 2 L 2 -||∇p n || 2 L 2 + 1 η Ω χ ω(t n+1 ) | wn+1 | 2 dx ≤ C||F n+1 || 2 L 2 + C||w n || 2 L 2
This last equation is multiplied by 2δt and written in k instead of n. Finally, the equations is summed from

k = 0 to n -1 with n ≤ N = E T δt
where E denotes the floor function, and we deduce:

||w n || 2 L 2 + δtε||p n || 2 L 2 + δt 2 ||∇p n || 2 L 2 + n-1 k=0 || wk+1 -w k || 2 L 2 + µ n-1 k=0 δt||∇ wk+1 || 2 L 2 + ε n-1 k=0 δt||p k+1 -p k || 2 L 2 + 2 η n-1 k=0 δt Ω χ ω(t k+1 ) | wk+1 | 2 dx ≤ ||w 0 || 2 L 2 + εδt||p 0 || 2 L 2 + δt 2 ||∇p 0 || 2 L 2 + 2C n-1 k=0 δt||F k+1 || 2 L 2 + 2C n-1 k=0 δt||w k || 2 L 2 (24)
It implies:

||w n || 2 L 2 ≤ f n + n-1 k=0 g k ||w k || 2 L 2 with,        f n = ||w 0 || 2 L 2 + εδt||p 0 || 2 L 2 + δt 2 ||∇p 0 || 2 L 2 + 2C n-1 k=0 δt||F k+1 || 2 L 2 g k = 2Cδt
The discrete Gronwall lemma 2.1 thus gives the following upper bound on ||w n ||:

||w n || 2 L 2 ≤ f n (1 + 2CT exp(2CT ))
Going back to (24) we deduce there exists C > 0 such that:

δtε||p n || 2 L 2 + δt 2 ||∇p n || 2 L 2 + n-1 k=0 || wk+1 -w k || 2 L 2 + µ n-1 k=0 δt||∇ wk+1 || 2 L 2 + ε n-1 k=0 δt||p k+1 -p k || 2 L 2 + 2 η n-1 k=0 δt Ω χ ω(t k+1 ) | wk+1 | 2 dx ≤ C (25)
which concludes the proof of Proposition 3.1 and then of Proposition 1.1 as, by assumptions, the quantities f n are uniformly bounded with respect to n < N and to δt.

Lemma 3.1. Under the hypothesis of Proposition 1.1, we have i) The divergence of w δt lies in L 2 (]0, T [; L 2 ) and there exists C > 0 such that for any ε > 0,

||div w δt || L 2 (]0,T [;L 2 ) ≤ C √ ε ( 26
)
As ψ is divergence free, the same inequality holds for v δt which implies the strong convergence of div(v δt ) towards 0 when ε tends to 0.

ii) ŵδt is bounded in L 2 (]0, T [; H -1 ) and || ŵδt || 2 L 2 (]0,T [;H -1 ) ≤ Cδt δt ε
Proof of i). From the pressure equation ( 8) we have ε(p n+1p n ) = -div(w n+1 ). Then we have

N -1 k=0 δt||div(w k+1 )|| 2 L 2 = ε N -1 k=0 εδt||p k+1 -p k || 2 L 2 .
and we deduce (26) exploiting the stability result (25).

Proof of ii). The second point is proven using the correction equation [START_REF] Simon | Compact sets in the space l p (0, t; b)[END_REF]. Taking the H -1 -norm we obtain:

|| ŵk+1 || H -1 ≤ δt ε ||div(w k+1 )|| L 2 . ( 27
)
Therefore, summing the square of this inequality from k = 0 to N -1 and using the bound of the velocity's divergence (26), we finally obtain:

N -1 k=0 δt|| ŵk+1 || 2 H -1 ≤ Cδt 2 ε .
Lemma 3.2. Under the hypothesis of Proposition 1.1, we have i) The velocity translation satisfies:

N -1 k=0 ||w k+1 -w k || 2 H -1 ≤ 2C δt ε ε δt + 1 (28) ii) ∇w δt is bounded in L 2 (]0, T [; L 2 ) and N -1 k=0 δt||∇w k+1 || 2 L 2 ≤ C 1 + δt ε
Proof of i). The stability result (25) gives a bound on the difference between the predicted velocity at the current time step and the velocity at the previous time step. Using the embedding of L 2 in H -1 we deduce:

N -1 k=0 || wk+1 -w k || 2 H -1 ≤ C.
Then, combining (27) with Lemma 3.1 the following inequality holds:

N -1 k=0 ||w k+1 -w k || 2 H -1 ≤ 2 N -1 k=0 ||w k+1 -wk+1 || 2 H -1 + || wk+1 -w k || 2 H -1 ≤ 2C δt ε 1 + ε δt
Proof of ii). To prove this point, we take ŵk+1 as a test function in the correction step ( 13)

|| ŵk+1 || 2 L 2 + δt ε ||div( ŵk+1 )|| 2 L 2 = - δt ε div( wk+1 ), div( ŵk+1 ) L 2 ≤ δt 2ε ||div( ŵk+1 )|| 2 L 2 + δt 2ε ||div( wk+1 )|| 2 L 2 .
We thus obtain an estimate on the corrected velocity ŵk+1 and its divergence:

|| ŵk+1 || 2 L 2 + δt 2ε ||div( ŵk+1 )|| 2 L 2 ≤ δt 2ε ||div( wk+1 )|| 2 L 2 . ( 29
)
Using that the norm

|| • || H 1 is equivalent to the norm || • || 2 L 2 + ||div(•)|| 2 L 2 + ||curl(•)|| 2 L 2 1
2 and curl( ŵk+1 ) = 0, we obtain:

|| ŵk+1 || 2 H 1 ≤ C || ŵk+1 || 2 L 2 + ||div( ŵk+1 )|| 2 L 2 ≤ δt 2ε + 1 ||∇ wk+1 || 2 L 2
The previous inequality is summed up from k = 0 to N -1. The predicted velocity gradient is bounded using the stability result (25). Then, we can find an upper bound on the total velocity gradient:

N -1 k=0 ||∇w k+1 || 2 L 2 δt ≤ 2 N -1 k=0 δt||∇ wk+1 || 2 L 2 + 2 N -1 k=0 δt||∇ ŵk+1 || 2 L 2 ≤ 2C + 2 1 + δt 2ε N -1 k=0 δt||∇ wk+1 || 2 L 2 ≤ 2 2 + δt 2ε C
4 Convergence analysis when ε and δt tend to 0

A stability result has been obtained in the previous section. The main purpose of this section is to establish Theorem 1.1 which can be write as the following convergence theorem, when ε and δt tend to 0 with ε = λδt.

Theorem 4.1 (Convergence when ε and δt tend to 0). We assume that the hypothesis of Theorem 1.1 are satisfied. Then, up to a subsequence, (v n , p n ) n solution of (3)-( 7) converges towards (v, p) weak solution of the penalized Navier Stokes problem (1) when ε and δt tend to 0 with ε = λδt. Furthermore, v and p satisfy:

v ∈ L ∞ (]0, T [; L 2 ) ∩ L 2 (]0, T [; H 1 ), p ∈ W -1,∞ (]0, T [; L 2 0
). Moreover, this solution is unique in two dimensional space.

Weak convergence of the velocity

We first establish the following result: Lemma 4.1. Under the hypothesis of Proposition 1.1, there exists v ∈ L 2 (]0, T [; L 2 ) (respectively ṽ ∈ L 2 (]0, T [; L 2 )) such that, up to a subsequence, v δt (resp. ṽδt ) weakly converges towards v (resp. ṽ) when ε and δt tend to 0

with ε = λδt : i) (v δt ) δt -⇀ v weakly in L 2 (]0, T [; H 1 ).
ii) (ṽ δt ) δt -⇀ ṽ weakly in L 2 (]0, T [; H 1 ) Moreover, at the limit ṽ = v Proof. This result directly comes from the stability study. Indeed, the regularity of ψ (2), Proposition 1.1 and lemma 3.2 ensures that (v δt ) δt is uniformly bounded in L 2 (]0, T [; H 1 ) provided that ε = λδt. Therefore, we can extract a subsequence still denoted v δt that weakly converges towards a function v in L 2 (]0, T [; H 1 ). In the same way, the sequence (ṽ δt ) δt is bounded in L 2 (]0, T [; H 1 ). Then we can extract a subsequence that weakly converges towards ṽ in L 2 (]0, T [; H 1 ) when δt and ε tend to 0 with ε = λδt.

Let us show that ṽ = v. From lemma 3.1 when ε = λδt ŵδt := w δt -wδt -→ δt→0 0 strongly in L 2 (]0, T [; H -1 ). Moreover, Lemma 3.2 and Proposition 1.1 ensures that w δt -wδt is bounded in L 2 (]0, T [; H 1 ). Therefore,

w δt -wδt -⇀ 0 weakly in L 2 (]0, T [; H 1 ) (30)
By unicity of the limit, we conclude that w = w and ṽ = v.

Strong convergence of the velocity

In order to use Simon's results (Theorem 2.1), let us estimate w δt in an appropriate Nikolskii space. The following lemma ensures that w δt belongs to the Nikolskii space N 1 2

2 (]0, T [; H -1 ) Lemma 4.2. Let C M be a positive constant and h > 0. Let u defined on a time interval [0, T ) with values in H -1 . We denote u k the value of u at the time t k and u δt the step function defined as in [START_REF] Holte | Discrete gronwall lemma and applications[END_REF]. Moreover we assume that the following uniform upper bounds hold:

N -1 k=0 ||u k+1 -u k || 2 H -1 ≤ C M sup k≤N ||u k || 2 H -1 ≤ C M and, ||(v δt .∇)ṽ δt || L 1 (]0,T [;L 3 2 ) ≤ ||v δt || L 2 (]0,T [;L 6 ) ||∇ṽ δt || L 2 (]0,T [;L 2 )
Thanks to Proposition 3.1 and Sobolev embeddings in three dimensions, the r.h.s. in the above estimates is uniformly bounded. Then using interpolation theorems, we deduce that (v δt .∇)ṽ δt is uniformly bounded in 5 ). Therefore, in this space, the sequence ((v δt .∇)ṽ δt ) δt weakly converges towards a function g that remains to determine. Combining the continuity of the application (w, v)

L 4 3 (]0, T [; L 6
→ (w.∇)v from L 2 (]0, T [; L 2 )×L 2 (]0, T [; H 1 ) in L 1 (]0, T [; L 1 )
with the convergences of Lemma 4.1 and Lemma 4.3, we deduce the convergence of the inertia term in the space

L 1 (]0, T [; L 1 ): (v δt .∇)ṽ δt -⇀ (v.∇)v in L 1 (]0, T [; L 1 )
The weak convergence in the smaller space 5 ) yields g = (v.∇)v and:

L 4 3 (]0, T [; L 6
(v δt .∇)ṽ δt -⇀ (v.∇)v weakly in L 4 3 (]0, T [; L 6 5 ) ( 35
)
We now study the convergence of div(v δt )ṽ δt . First, we know that:

||ṽ δt div(v δt )|| L 2 (]0,T [;L 1 ) ≤ ||∇v δt || L 2 (]0,T [;L 2 ) ||ṽ δt || L ∞ (]0,T [;L 2 ) . (36)
From the pressure equation, the following inequality holds:

|| ŵn+1 || 2 L 2 ≤ 2δt 2 ||∇p n || 2 L 2 + ||∇p n+1 || 2 L 2 (37)
Then, the stability result (25) ensures for all n ≤ N that δt 2 ||∇p n || 2 ≤ C and we obtain a bound of ŵδt in L ∞ (]0, T [; L 2 ). Since w δt also lies in L ∞ (]0, T [; L 2 ), we deduce that wδt is bounded in L ∞ (]0, T [; L 2 ). Furthermore, from Hölder's inequality, a second inequality holds:

||ṽ δt div(v δt )|| L 1 (]0,T [;L 3 2 ) ≤ ||∇v δt || L 2 (]0,T [;L 2 ) ||ṽ δt || L 2 (]0,T [;L 6 )
Finally, from interpolation theorems we obtain a bound of ṽδt div(v δt ) in 5 ) and consequently the weak convergence of div(v δt )ṽ δt towards a function g that remains to determine. The operator (u, v) → 1 2 div(u)v is continuous from H 1 × L 2 to L 1 so using the previous convergences, div(v δt )ṽ δt converges towards div(v)v in the space L 1 (]0, T [; L 1 ). As the product is bounded in the smaller space L

L 4 3 (]0, T [; L 6

The two dimensional case

For a two dimensional space, we can demonstrate the convergence in a higher regularity space. Indeed, as 2 ). Yet, 3 ). Using the same arguments as above we deduce 3 ). It concludes the proof of Lemma 4.4 in the two dimensional case.

v δt ∈ L ∞ (]0, T [; L 2 ) ∩ L 2 (]0, T [; H 1 ) then by interpolation v δt lies in L 4 (]0, T [; H 1
H 1 2 is embedded into L 4 , therefore v δt ∈ L 4 (]0, T [; L 4 ). Then (v δt .∇)ṽ δt is bounded in L 4 3 (]0, T [; L 4
(v δt .∇)ṽ δt -⇀ (v.∇)v weakly in L 4 3 (]0, T [; L 4 3 ) and, ṽδt div(v δt ) -⇀ vdiv(v) weakly in L 4 3 (]0, T [; L 4

Proof of Theorem 4.1

We can now pass to the limit in the numerical scheme. Let φ ∈ V (for the notations, see Section 2.1). The numerical scheme (3)-( 7) reads in variational formulation:

d dt (v δt (t), φ) + 2µ(D(ṽ δt ), D(φ)) + (B(v δt (t -δt), ṽδt (t)), φ) + 1 η (χ ω(t) (ṽ δt (t) -v s (t)), φ) = (f (t), φ)
We multiply by a function θ ∈ C 1 (0, T ) such that θ(T ) = 0 and we integrate from 0 to T . We do not have any information on the time derivative of the velocity. Therefore, the temporal term is integrated by part so that the time derivative holds on θ:

- T 0 (v δt (t), φ) θ ′ (t) dt -(v δt (0), φ) θ(0) + 2µ T 0 (D(ṽ δt (t)), D(φ)) θ(t) dt + T 0 (B(v δt (t -δt), ṽδt (t)), φ) θ(t) dt + 1 η T 0 (χ ω(t) (ṽ δt (t) -v s (t)), φ) θ(t) dt = T 0 (f δt , φ) θ(t) dt (39)
We pass to the limit δt → 0 in this last equation with ε = λδt using Lemma 4.1 and Lemma 4.4. It gives:

- T 0 (v(t), φ) θ ′ (t) dt -(v(0), φ) θ(0) + 2µ T 0 (D(v(t)), D(φ)) θ(t) dt + T 0 ((v(t).∇)v(t), φ) θ(t) dt + 1 η T 0 (χ ω(t) (v(t) -v s (t)), φ) θ(t) dt = T 0 (f (t), φ) θ(t) dt (40)
Since ε also tends to 0 then from Lemma 3.1, we have at the limit div(v) = 0 on Ω.

Applying the above equality for θ ∈ D(0, T ), we deduce the following equality in V ′ :

- T 0 v(t) θ ′ (t) dt = T 0 div(2µD(v(t))) θ(t) dt - 1 η T 0 χ ω(t) (v(t) -v s ) θ(t) dt - T 0 (v(t).∇)v(t) θ(t) dt + T 0 f (t) θ(t) dt (41)
The operator

L : u → div(2µD(u)) (respectively B : (u, v) → B(u, v)) is continuous from V to V ′ (resp. from V × V to V ′ ).
Therefore, there exists C > 0 such that:

T 0 ||div(2µD(v))|| V ′ dt ≤ C T 0 ||v|| V dt ≤ C √ T ||v|| L 2 (]0,T [;V) . T 0 ||(v.∇)v|| V ′ dt ≤ C T 0 ||v|| 2 V dt ≤ C||v|| 2 L 2 (]0,T [;V) . (42)
As ( 41) is valid for any θ ∈ D(]0, T [) we deduce that v has a weak derivative in time which lies in L 1 (]0, T [; V ′ ) and for almost every t ∈]0, T [:

∂v ∂t -div(2µD(v)) + (v.∇)v + 1 η χ ω(t) (v -v s ) = f in V ′ . (43)
We now need to recover the initial data. Since ∂v ∂t belongs to L 1 (]0, T [; V ′ ) and v belongs to L 2 (]0, T [; V), we show using Proposition 2.1 that v is continuous with values in V ′ for the strong topology. Furthemore, by hypothesis v(0) = v 0 in the weak continuity sense with values in V ′ . Therefore, the initial condition v(0) = v 0 is verified in the strong sense because the weak limit is unique.

From De Rham theorem, we can now deduce the existence of the pressure. Let G(t) be defined by:

G(t) = -div(2µD(v)) + (v.∇)v + 1 η χ ω(t) (v -v s ) -f
Thanks to (43), for almost every t ∈]0, T [, dv dt , φ

V ′ ,V + < G(t), φ > H -1 ,H 1 0 = 0.
We integrate this last equation from 0 to t. It gives:

< v(t), φ > L 2 -< v(0), φ > L 2 + t 0 G(τ ) dτ , φ
It can be written under the form: < K(t), φ > H -1 ,H 1 0 = 0, where

K(t) = v(t) -v(0) - t 0 div(2µD(v)) dτ + t 0 (v.∇)v dτ + 1 η t 0 χ ω(τ ) (v -v s ) dτ - t 0 f dτ .
Note that K is weakly continuous in time with values in H -1 . Therefore, for all t ∈]0, T [ we deduce from De Rham theorem the existence of π(t) ∈ L 2 0 such that:

K(t) = -∇π(t)
Following the work of [START_REF] Boyer | Mathematical tools for the study of the incompressible Navier-Stokes Equations and related models[END_REF] (chapter V), we show that t → π(t) is weakly continuous in time with values in L 2 . In particular, π lies in the space L ∞ (]0, T [; L 2 0 ). Indeed, if g ∈ L 2 , there exists h ∈ H 1 0 such that div(h) = gm(g)

where m(g) denotes the mean value of g on Ω. Then, (π(t), g)

L 2 = (π(t), g -m(g)) L 2 because m(π) = 0 = (π(t), div(h)) L 2 = -(∇π(t), h) H -1 ,H 1 0 = (K(t), h) H -1 ,H 1 0
This quantity is continuous because K is weakly continuous in time with values in H -1 . We can then introduce the distribution p = ∂π ∂t which lies in the space W -1,∞ (]0, T [; L 2 0 ). Taking test functions under the form ∂φ ∂t with φ ∈ D(]0, T [×Ω), we show that the equation

∂v ∂t -div(2µD(v)) + (v.∇)v + 1 η χ ω(t) (v -v s ) + ∇p = f (44)
is satisfied in the sense of distributions.

In two dimensional space we can show the uniqueness of solutions of this equation using classical results (see [START_REF] Boyer | Mathematical tools for the study of the incompressible Navier-Stokes Equations and related models[END_REF] chapter V).

Convergence towards the Navier-Stokes Equations

The aim of this section is to study the convergence when η tends to 0. To do so, we consider the weak limit of the scheme when ε and δt tend to 0 which verifies (44) and indice the solution by η. Theorem 5.1 (Convergence when η tends to 0.). When η tends to 0, the sequence (v η ) η weakly converges towards a limit v which satisfies

v |∂ω(t) = v s (t)
Furthermore, v |Ω\ω(t) is a weak solution of the Navier Stokes equations on Ω \ ω(t). Moreover, there exists h ∈ W ′ such that:

1 η χ ω(t) (v -v s ) -⇀ h in W ′ (45)
We first prove the following lemma:

Lemma 5.1. For all t ∈]0, T [, we have

t 0 ||v η -v s || 2 L 2 (∂ω(τ )) dτ ≤ Cη 1 2
Proof. For all t ∈]0, T [,

||v η -v s || 2 L 2 (∂ω(t)) ≤ C||v η -v s || L 2 (ω(t)) ||v η -v s || H 1 (ω(t))
We integrate this last inequality from 0 to t and obtain:

t 0 ||v η -v s || 2 L 2 (∂ω(τ )) dτ ≤ C t 0 ||v η -v s || L 2 (ω(τ )) ||v η -v s || H 1 (ω(τ )) dτ ≤ C   t 0 ||v η -v s || 2 L 2 (ω(τ )) dτ   1 2   t 0 ||v η -v s || 2 H 1 (ω(τ )) dτ   1 2 ≤ C   t 0 ||v η -v s || 2 L 2 (ω(τ )) dτ   1 2 ||w η || L 2 (]0,T [;H 1 )
From Lemma 3.2, w δt is bounded in H 1 . Moreover, this bound is uniform in δt ε . Using the lower of the norm for the weak topology, we obtain that w η is also bounded in H 1 . Finally, from the energy estimates of Section 3, we obtain:

t 0 ||v η -v s || L 2 (∂ω(τ )) dτ ≤ Cη 1 2
Therefore, when η → 0, the velocity on the immersed boundary ∂ω(t) tends towards the obstacle velocity v = v s in the space L 2 (]0, T [; L 2 (∂ω(t))).

Proof of theorem 5.1. From (44) we have:

∂v η ∂t -div(2µD(v η )) + (v η .∇)v η + ∇p η + 1 η χ ω(t) (v η -v s ) = f Let φ ∈ W. We have: T 0 1 η χ ω(t) (v η -v s ), φ dt = T 0 v η , ∂φ ∂t V,V ′ dt + (v η (0), φ(0)) - T 0 2µ(D(v η ), D(φ)) + ((v η .∇)v η + ∇p η + f, φ) dt (46)
We use the continuity of the divergence and inertia operators as in (42) and that the estimates obtained through the stability study are uniform in η. The following inequality holds:

T 0 1 η (χ ω(t) (v η -v s ), φ) dt ≤ C||φ|| W
This is true for all test functions φ ∈ W, we deduce that

1 η χ ω(t) (v η -v s ) -⇀ h weakly in W ′
We can now pass to the limit in (46) using the continuity of the inertia and diffusion operators as well as the lower semi-continuity of the norm. It gives:

T 0 ∂v ∂t , φ dt + T 0 < -div(2µD(v)) + (v.∇)v + ∇p -f, φ > + < h, φ > dt = 0 ( 47
)
As h is the weak limit of 1 η χ ω(t) (v ηv s ), then for every function φ such supp(φ(t, •)) ⊂ Ω \ ω(t) for all t:

< h, φ >= 0 At the limit when η tends to 0, the velocity on the immersed boundary tends to the solid velocity (see Lemma 5.1). We find back a Dirichlet boundary condition on the obstacle boundary.

A Kolmogorov lemma

Lemma A.1. Let C M a positive constant and h > 0. Let u defined on a time interval [0, T ) with values in H -1 . We denote u k the value of u at the time t k and u δt the step function defined as in [START_REF] Holte | Discrete gronwall lemma and applications[END_REF]. Moreover we assume that the following conditions are verified:

N -1 k=0 ||u k+1 -u k || 2 H -1 ≤ C M sup k≤N ||u k || 2 H -1 ≤ C M
Then, there exists C > 0 independent of δt such that:

T -h 0 ||u δt (t + h) -u δt (t)|| H -1 dt ≤ Ch 1 2 (48)

And,  

T -h 0 ||u δt (t + h) -u δt (t)|| 2 H -1 dt   1 2 ≤ Ch 1 2 (49)
Proof. We distinguish two cases, h ≤ δt and h > δt

u δt (t + h) -u δt (t) = u k+1 -u k (50)
Consequently, on each interval [t k , t k+1 [ the function u(t + h)u(t) is non-null on an intervall of size h.

We have:

T -h 0 ||u δt (t + h) -u δt (t)|| H -1 dt ≤ N -1 k=0 h||u k+1 -u k || H -1 ≤ N -1 k=0 h 1 2 N -1 k=0 h||u k+1 -u k || 2 H -1 1 2 ≤ h 1 2 (T C M ) 1 2
(51)

And,

T -h

0 ||u δt (t + h) -u δt (t)|| 2 H -1 dt ≤ N -1 k=0 h||u k+1 -u k || 2 H -1 ≤ hC M (52)
2. h > δt:

We first consider the case where h is a multiple of δt. Let j = E h δt .

We denote by χ k the characteristic function of the interval [t k , t k+1 [. Then, on the interval [0; Th], the function u(t + h)u(t) can be expressed as:

u δt (t + h) -u δt (t) =                  j-1 k=0 u k+j χ k+j (t) -u k χ k (t) + 2j-1 k=j u k+j χ k+j (t) -u k χ k (t)
. . .

+ N -j-1 k=N -2j u k+j χ k+j (t) -u k χ k (t)                 
In the last sum many terms vanish and only the first and last terms are remaining.

u δt (t + h) -u δt =   - j-1 k=0 u k χ k (t) + N -1 k=N -j u k χ k (t)   (53)
χ k u k             (57)
If we reorganize the characteristic functions, the last equation reduces to: We finally integrate this result from 0 to Th:

u δt (t + h) -u δt (t) =                - j-1 k=0 χ k (t)u k
T -h 0 ||u δt (t + h) -u δt (t)|| H -1 dt ≤ j k=0 δt||u k || H -1 + N -1 k=N -j-1 δt||u k || H -1 ≤ 2jδt sup k≤N ||u k || H -1 ≤ 2hC (59)
Using the same arguments as above we show:

4 3 ( 6 5 4 3Thus for d = 3

 3643 ]0, T [; L ), the result below holds: ṽδt div(v δt ) -⇀ vdiv(v) weakly in L Lemma 4.4 follows from (35) and (38).

-N -j- 3 k=N -2j- 3 χ

 33 u j χ [t j ;h] (t) + u N -j-1 χ [T -2h;h1+(N -2j-1)δt] + [h1+kδt;h1+(k+1)δt] (t)u k+j+1 + u N -1 χ [h1+(N -j-2)δt;T -h]

  Therefore, if we integrate ||u δt (t+h)-u δt (t)|| H -1 between 0 and T -h we obtain considering||u δt (t + h)u δt (t)|| H -1 dt ≤ ChWe now consider the case where h is not a multiple of δt. We still denote j as E h δt and we seth 1 = (j + 1)δth.The function u δt (t + h) then reads:u j χ [t 0 ;h1] (t) + (t)u k+j+1 + u N -1 χ [h1+(N -j-2)δt;T -h] (t)Then the difference u(t + h)u(t) can be expressed as:u δt (t + h)u δt (t) =

	n-j-3
	k=0	χ [h1+kδt;h1+(k+1)δt]  u j χ [t 0 ;h1] (t)
			n-j-3
		  	+	χ [h1+kδt;h1+(k+1)δt] (t)u k+j+1
			k=0
		
		
		
		
		
		
							T -h	χ k (t) dt =
	δt for any k:						0
				j-1 k=0	δt||u k || H -1 + j-1	N -1 k=N -j	δt||u k || H -1	(54)
			≤ 2	k=0	δt sup k≤N	||u k || H -1
			≤ 2Ch
	Finally, after some manipulations:
		T -h					1 2	(55)
		0
	Moreover, we know that the characteristic functions have disjoint supports. Therefore:
	T -h 0 ||u δt (t + h) -u δt (t)|| 2 H -1 dt =	j-1 k=0	δt||u k || 2 H -1 + j-1	N -1 k=N -j	H -1 δt||u k || 2	(56)
			≤ 2	k=0	δt sup k≤N	||u k || 2 H -1
			≤ 2Ch
			20

T -h 0 ||u δt (t + h)u δt (t)|| H -1 dt = + u N -1 χ [h1+(N -j-2)δt;T -h] (t) -N -j-1 k=0

Then, there exists C > 0 independent of δt such that:

And,  

T -h

We postpone the proof of this lemma in appendix A. Proof. To prove this lemma, we first show that w δt is uniformly bounded in the Nikolskii space N 1 2

2 (]0, T [; H -1 ). Combining (28) with Proposition 3.1 i), from Lemma 4.2 we deduce the following bound on the translations:

which demonstrate that w δt belongs to the Nikolskii space N We know that w δt and ψ δt lie in L ∞ (]0, T [; L 2 ), therefore

We obtain the result for any p ∈ [2, +∞[ using interpolation properties.

Weak convergence of the inertia terms

We show the weak convergence of the inertia term B(v δt (tδt), ṽδt (δt)).

3 ) in two dimensions and (p, q) = ( 4 3 , 6 5 ) in three dimensions. We distinguish the cases of two and three dimensional spaces.


T Barois

S Perisanu

P Poncharal

P Vincent

S T Purcell

A Ayari

Quality factor enhancement of Nano electromechanical systems by capacitive driving beyond the resonance

Nano electromechanical systems are considered as ultra sensitive devices for mass and force detection. Capacitive actuation is widely used in these devices but is known to degrade the quality factor of the resonator due to DC electrostatic damping. We report the enhancement of the quality factor of SiC vibrating nanowires detected nano optomechanically and electrically by applying an AC capacitive driving at a frequency above both the resonance frequency and the electrical cut-off frequency. Self-oscillations are demonstrated for optimal conditions. We developed an analytical model of the phenomenon and showed that it can lead to an improvement of the force sensitivity.

I. INTRODUCTION

Nano electromechanical systems (NEMSs) [1], where electrons are coupled to mechanical motion, have generated a tremendous interest over the past years as they are promising candidates for testing macroscopic quantum physics [2][3][4], non-linear dynamics at the nanoscale [5,6] and for their capability in sensors applications. These systems have shown record limits of detection reaching yoctogram mass resolution [7] and a sensitivity of 12 zN/ √ Hz [8]. In these electromechanical devices, the continuous improvement of the electrostatic transduction and the elimination of different sources of damping will ultimately be limited by the electrostatic damping [9] as it is already the case for some single electron devices [10,11] or semiconducting nanowire cantilevers [12]. It is possible to reduce this damping by lowering the device resistance [12] or by applying a DC offset to compensate for residual voltages originating form work function differences and contact potentials [13], but such solutions are not always convenient. It has been recently theoretically proposed [14] that electrostatic coupling can also lead to a decrease of damping and to self-oscillations by simply applying an AC signal to a suspended electro mechanical graphene resonator and taking into account non-linear effects. The driving frequency is not related to the resonance frequency but just needs to be higher than the inverse of the charge relaxation time.

This method is appealing as (i) the operating frequency range is much less stringent than for systems coupled to a cavity [15] or a resonator [START_REF] Rocard | Dynamique générale des vibrations[END_REF][START_REF] Brown | [END_REF][18] (ii) it needs no additional elements than the already present coupling capacitor and NEMS resistance, (iii) it works for single and double clamped geometries, (iv) it gives a simple mean to tune the damping and to improve the signal to noise ratio for force sensing [19] as already demonstrated for instance for magnetometers [START_REF] Thompson | Solid-State Sensors, Actuators and Microsystems Conference[END_REF][START_REF] Kumar | Micro Electro Mechanical Systems (MEMS)[END_REF]. In this article, an experimental proof of this capacitively induced decrease of a NEMS mechanical damping is presented. We demonstrate that this effect can be large enough to induce self-oscillations in a mechanical resonator. The model of Eriksson et al. [14] is then extend as this effect raises interesting question about the influence of the electrostatic actuation itself to perform an accurate measurement of the damping in NEMS and about the noise it induces for cooling and sensing applications.

II. RESULTS

We studied 2 different SiC nanowires (samples NW1, length L = 34 µm, radius r = 75 nm; NW2, L = 86 µm, r = 90 nm) fixed at one end to W rigid tips and submitted to AC and DC capacitive forces (see Figure 1a). The measurements were performed in two different ultra high vacuum chambers one equipped with a scanning electron microscope In the electron chamber the mechanical motion of the nanowire can be imaged as in standard electron microscopy or measured by recording the secondary electron signal in spot mode [12]. Figure 1b is a SEM image of NW2 close to the counter electrode. The resonance frequencies of NW2 were measured with a lock-in amplifier for an applied AC voltage of 50 mV. The first transverse mechanical eigenmode showed two polarizations at 39.83 kHz and 40.01 kHz where the degeneracy is lifted by either the clamping condition or the nanowire asymmetry. The second and third eigenmodes were at 230 kHz and 620 kHz.

The width of the resonance peaks is about 5 Hz and outside these peaks no mechanical signal can be measured or observed with the SEM imaging. However when a higher AC signal above 18 V is applied, spontaneous mechanical vibrations where observed for driving frequencies ranging from 750 kHz to 850 kHz. The excitation of higher eigenmodes can be ruled out because the imaging of the mode doesn't show any nodes as observed in figure 1c. Non-linear high order parametric excitation of the first mode and parametric coupling between the first mode and high order mechanical mode [18] seem dubious because of the discrepancy between the theoretical expected frequency and the observed frequency range for spontaneous oscillations. For this type of sample, the resistance is of the order of 1 GΩ and the capacitance in the fF range [12,[START_REF] Barois | [END_REF]. This corresponds to an inverse of the charge relaxation time of the order of hundreds of kHz. So our observations are in agreement with the model of Eriksson et al. [14] predicting self-oscillations for a driving frequency higher than the charge relaxation time.

We performed additional measurements with NW1 in our UHV optical set up in order to get a more quantitative agreement with theory. In the optical chamber, a 35 µW and 633 nm laser beam is strongly focused on the free end of the nanowire and the transmitted light is collected using a differential measurement based on a dual quadrant photodiode [23,24] giving access to the thermo-mechanical motion of the resonator with high signal to noise ratio. So, in this experimental set up, no AC probe signal is required to measure the mechanical motion and electron beam induced self-sustained oscillations [25] will not interfere with our measurements. The first transverse mechanical eigenmode showed two polarizations at 113.434 kHz and 114.071 kHz. In the rest of the text, we will focus only on the low frequency polarization as the high frequency polarization shows the same behavior. An AC signal of 10 V is applied to the nanowire and swept back and forth between 150 kHz and 2 MHz in 25 kHz increments while recording the power spectrum density of the thermomechanical motion as plotted in figure 2a. A shift in frequency and an increase in vibration amplitude is observed when the frequency increases. The resonance peaks have been fitted with a lorentzian (see figure 1 in supplemental material). Figure 2b shows the variation of the quality factor defined as ω 0 /Γ AC where ω 0 /2π is the resonance frequency and Γ AC the resonance width. The dashed line indicates the initial quality factor Q 0 = 30 600 in the absence of AC or DC voltages. At low driving frequencies, the quality factor is strongly reduced but increases with the frequency until it reach reach a maximum above 40 000 around 650 kHz. This maximum is significantly higher than the initial quality factor. At higher frequencies, Q tends to Q 0 . At first sight, this effect seems rather counter intuitive as the electrostatic interaction usually induces an increase in damping when a DC voltage is applied and AC actuation is very often used as a reliable mean to drive NEMS. Furthermore the possibility to increase or decrease the damping is reminiscent of coupling with a blue-detuned or red-detuned cavity whereas no cavity is involved in this system. The mechanical system can be modelized as a mechanical resonator capacitively connected to its environment.

mẍ + mΓ 0 ẋ + mω 2 0 x = C ′ 2 (U tot ) 2 + F T (t) ( 1
)
where x is the position in the transverse direction of the free end of the nanowire, m is the effective mass of the mechanical resonator, the dot refers to the derivative versus time t, Γ 0 is the intrinsic mechanical damping, F T (t) is the Langevin force coming from the thermal bath, C ′ is the space derivative of the capacitance C formed between the resonator and the electrostatic environment, U tot is the total voltage difference between the point x and the ground (see Fig. 1a). This voltage is provided by a transmission line that can be roughly described by the following equation :

RC Utot + U tot = -RC ′ U tot ẋ -RC ′ Utot x + V AC cos Ωt ( 2
)
where R is the resistance of the resonator including the contact resistance, V AC is the applied AC voltage at frequency Ω/2π. The important terms are those related to the change in capacitance when the nanowire moves and changes its distance from the counter electrode.

This motion induces a motional current C ′ U tot ẋ + C ′ Utot x. These non-linear terms mix the driving voltage at angular frequency Ω with the nanowire motion at angular frequency ω.

A redshifted (Stokes) electromagnetic wave at ω -= Ωω and a blueshifted (anti-Stokes) electromagnetic wave at ω + = Ω + ω are then produced and down-mixed by the squared term of the electrostatic actuation in Eq. 1 (see Figure 3 and supplemental material). This induces a feedback mechanism that modify the damping and resonance frequency of the nanowire. At low driving frequencies, the electrical circuit should behave as if a DC voltage was applied and so an increase in damping is expected. The electrical actuation is a low pass filter, so for frequencies beyond the inverse of the charge relaxation time τ = RC (i.e.

the cut-off frequency), the effect of the applied voltage on the apex of the nanowire should vanish and the initial damping should be recovered . However, in the intermediate regime of frequency, the actuation is dephased compared to the DC actuation and the amplitude of the electromagnetic actuation although smaller can still influence the resonator. This dephasing allows a decrease of damping for a wide frequency range. The effective mechanical equation is then :

ẍ + Γ AC ẋ + ω 2 AC x = 0 (3)
with :

Γ AC = Γ 0 + g 2 1 V 2 AC ℜ Ω ℑ ω + -ℑ ω - 2ω ; ℜ Ω = 1 1 + (Ωτ ) 2 ; ℑ Ω = Ωτ 1 + (Ωτ ) 2 ; g 1 = C ′ √ 2mC (4)
where g 1 is the electrostatic coupling strength for an applied voltage of 1V. It characterizes the strength of the electrostatic coupling. The blue-shifted term ℑ ω + always leads to an increase in damping whereas the red-shifted term ℑ ω -reduces the damping as long as ω < Ω. So this mechanism can not lead to a reduction in damping for a driving lower than the resonance frequency. For ω << Ω the damping can be approximated by :

Γ AC ≈ Γ 0 + g 2 1 V 2 AC τ 1 -(τ Ω) 2 (1 + (τ Ω) 2 ) 3 (5)
For Ωτ > 1, the AC driving starts to decrease the damping. The quality factor enhancement is maximum for Ω = √ 2/τ where the damping decreases by g 2 1 V 2 AC τ /27. The data in Fig. 2b were fitted with Eq. 5. A very good agreement is observed between theory and experiment. It gives a charge relaxation time of 347 ns corresponding to a cut of frequency of 460 kHz, a maximum quality factor enhancement at 650.5 kHz and a coupling strength of 2098.5 rad/s/V corresponding to 333.98 Hz/V. Thus our system is in the weak coupling regime g 1 V AC τ ≈ 0.007 <<1 and this explains why a rather important voltage of 10 V was necessary to observe this effect. A much more spectacular effect should be expected in the strong coupling regime. As a comparison, we performed a measurement of the resonance width as a function of the DC voltage without an AC voltage. As expected [12], a square dependence in voltage was observed (see supplemental material) and fitted. From this fit and using a value of τ = 347 ns as determined previously, the coupling strength was estimated at 1622.5 rad/s/V in agreement with the g 1 value obtained from the AC experiment.

An annealing of the nanowire was performed in order to remove surface adsorbates and increase the initial quality factor up to 66 700. An AC signal of 10 V was applied to the nanowire and swept back and forth between 200 kHz and 1.5 MHz in 50 kHz increments while recording the power spectrum density of the thermo-mechanical motion as plotted in Fig. 4a. Figure 4b shows a quality factor above 150 000 between 650 kHz and 850 kHz.

The peak width could not be measured accurately in this frequency range due to the 1 Hz limitation of the spectrum analyzer bandwidth. If these points are excluded, a fit is still possible and gives similar values for g 1 and τ as before the annealing. The clear departure of the quality factor from the fit between 650 kHz and 850 kHz indicates that the nanowire self-oscillates. It can be pointed out here, that this is a true self-oscillation where the phase is free, contrary to the instability of the parametric oscillator where the phase is fixed by the AC generator.

III. DISCUSSION

An interesting point, concerning our derivation of equation 4 is that it gives a simple expression of the damping valid for any Ω (except for Ω=0 where a slight correction is needed).

In particular, it appears from Eq. 4 that when Ω is close to the resonance frequency, the damping obtained with a capacitive AC driving is overestimated by g 2 1 V 2 AC τ /(1 + (τ ω AC ) 2 ) 3 . This overestimation can be even higher when mixing technics are used. So, in the strong coupling regime and in the study of non-linear effects like the Duffing regime, the change in damping might have to be taken into account. Then, it can be expected that a change in quality factor with temperature could be simply related to for instance a change in the contact resistance. This could have implication also for the study of NEMS frequency noise by applying a capacitive driving close to the resonance and measuring the vibration noise spectrum as proposed by Zhang et al. [26].

In our particular experiment at room temperature, although the transduction was improved, this doesn't lead to an improvement in the force sensitivity as the signal and thermal vibrations are amplified by an equal amount. Nevertheless at lower temperature or in systems where thermal vibrations are not the dominant source of noise a gain in sensitivity could be expected as long as the Johnson noise from the resistance is not too high. In the supplemental material, this additional force noise is estimated to be of the same order as the thermal force noise for NW1. The effective temperature was also calculated and we showed that, contrary to the DC electrostatic damping, the AC mechanism leads to heating of the mechanical mode. This heating can be in our case as high as 4 times the room temperature.

IV. CONCLUSION

In conclusion, we have shown experimentally that applying a capacitive AC signal can improve the quality factor of a nanowire mechanical resonator and lead to self-oscillations.

The enhancement of the quality factor was maximum for a driving angular frequency at √ 2 the inverse of the charge relaxation time. We derived an analytical expression for this effect valid on a wide frequency range and calculated the influence of the AC driving on the force sensitivity of the device. This enhancement is not limited to single clamped geometry or nanowires but is more suitable for resistive sample such as semiconducting 2D sample like MoS 2 or samples contacted by a tunneling resistance. The simplicity of this mechanism opens up a new path for the simultaneous driving and ultimately synchronization [27] of NEMS array as a single AC driving can lead to self-oscillations of all the devices even when they are not perfectly identical. This technique can be seen also as a mean to measure the charge relaxation time of a NEMS in the case where cross talk make difficult a direct electrical detection.

FIG. 1 .

 1 FIG. 1. (a) Schematic of a vibrating SiC nanowire on a tungsten tip submitted to an AC voltage V AC or a DC voltage V DC . (b) SEM image of the silicon carbide nanowire NW2 close to a gold counter electrode for V AC = V DC = 0 V. (c) SEM image of NW2 in self-oscillation for V AC = 19 V, V DC = 0 V and a driving frequency Ω/2π = 810 kHz well above the resonance frequency of 40 kHz.

(FIG. 2 .

 2 FIG. 2. (a) Evolution of the displacement noise spectrum of NW1 as a function of the driving frequency Ω/2π for V AC = 10 V and V DC = 0 V. (b) Evolution of the quality factor Q as a function of the driving frequency. The solid line is a fit of the quality factor with equation 4. The dashed line is the quality factor without an AC voltage.

FIG. 3 .

 3 FIG. 3. (a) Schematic of the mechanism involved in the quality factor enhancement induced by an AC voltage. F ext is an external force like the thermal noise force. F tot is the sum of the external force and the two feedback forces. χ(ω) is the mechanical susceptibility of the nanowire. H(ω)e iϕ is the transfer function of the transmission line. U 0 and U 1 are calculated in the supplemental material. The two crossed circles represent signal mixers responsible for the non-linear effects and originating from the C' terms. (b) Evolution of the term ℜ Ω as un function of the driving frequency with the experimental parameters of NW1. The smaller the black arrow, the smaller the effect on the damping. (c) Evolution of the term ℑ Ω as un function of the driving frequency with the experimental parameters of NW1. The damping is reduced as long as the blue arrow is longer than the red arrow.

FIG. 4 .

 4 FIG. 4. (a) Evolution of the displacement noise spectrum of NW1 after annealing as a function of the driving frequency Ω/2π for V AC = 10 V and V DC = 0 V. Satellite peaks related to selfoscillation appear between 650 kHz and 850 kHz. (b) Evolution of ω 0 /Γ AC (which correspond to the quality factor Q outside the self-oscillation region) as a function of the driving frequency after annealing. The solid line is a fit of the quality factor with Eq. 4. The dashed line is the quality factor without an AC voltage. (c) Evolution of the vibration amplitude at resonance as a function of the driving frequency.

ACKNOWLEDGMENTS

The authors thank the Plateforme Nanofils et Nanotubes Lyonnaise of the University Lyon1 and are especially grateful to the NOF team for their support. This work was supported by the French National Research Agency (ANR) through its Blanc Programme (FO-CUS, ANR-13-BDS10-0012-01).


Ion Grama

AND Quansheng Liu
email: quansheng.liu@univ-ubs.fr

Eric Miqueu
email: eric.miqueu@univ-ubs.fr

Liu

ASYMPTOTIC OF THE DISTRIBUTION AND HARMONIC MOMENTS FOR A SUPERCRITICAL BRANCHING PROCESS IN A RANDOM ENVIRONMENT

Keywords: Mathematics Subject Classification. Primary 60J80, 60K37, 60J05. Secondary 60J85, 92D25 Branching processes, random environment, harmonic moments, asymptotic distribution, decay rate. Current address

Let (Z n ) be a supercritical branching process in an independent and identically distributed random environment ξ. We show the exact decay rate of the probability P(Z n = j|Z 0 = k) as n → ∞, for each j k, assuming that P(Z 1 = 0) = 0. We also determine the critical value for the existence of harmonic moments of the random variable W = lim n→∞ Zn E(Zn|ξ) under a simple moment condition.

Résumé. Soit (Z n ) un processus de branchement surcritique en environnement aléatoire ξ indépendant et identiquement distribué. Nous donnons un équivalent de la probabilité P(Z n = j|Z 0 = k) lorsque n → ∞, pour tout j k, sous la condition P(Z 1 = 0) = 0. Nous déterminons également la valeur critique pour l'existence des moments harmoniques de la variable aléatoire limite W = lim n→∞ Zn E(Zn|ξ) , sous une hypothèse simple d'existence de moments.

Introduction

A branching process in a random environment (BPRE) is a natural and important generalisation of the Galton-Watson process, where the reproduction law varies according to a random environment indexed by time. It was introduced for the first time in Smith and Wilkinson [START_REF] Smith | On branching processes in random environments[END_REF] to modelize the growth of a population submitted to an environment. For background concepts and basic results concerning a BPRE we refer to Athreya and Karlin [START_REF] Athreya | On branching processes with random environments: I: Extinction probabilities[END_REF][START_REF] Athreya | Branching processes with random environments: II: Limit theorems[END_REF]. In the critical and subcritical regime the branching process goes out and the research interest is mostly concentrated on the survival probability and conditional limit theorems, see e.g. Afanasyev, Böinghoff, Kersting, Vatutin [START_REF] Afanasyev | Limit theorems for weakly subcritical branching processes in random environment[END_REF][START_REF] Afanasyev | Conditional limit theorems for intermediately subcritical branching processes in random environment[END_REF], Vatutin [START_REF] Vatutin | A refinement of limit theorems for the critical branching processes in random environment[END_REF], Vatutin and Zheng [START_REF] Vatutin | Subcritical branching processes in a random environment without the Cramer condition[END_REF], and the references therein. In the supercritical case, a great deal of current research has been focused on large deviations, see Bansaye and Berestycki [START_REF] Bansaye | Large deviations for branching processes in random environment. Markov Process[END_REF], Bansaye and Böinghoff [START_REF] Bansaye | Upper large deviations for branching processes in random environment with heavy tails[END_REF][START_REF] Bansaye | Lower large deviations for supercritical branching processes in random environment[END_REF][START_REF] Bansaye | Small positive values for supercritical branching processes in random environment[END_REF], Böinghoff and Kersting [START_REF] Böinghoff | Upper large deviations of branching processes in a random environment -offspring distributions with geometrically bounded tails[END_REF], Huang and Liu [START_REF] Huang | Moments, moderate and large deviations for a branching process in a random environment[END_REF], Nakashima [START_REF] Nakashima | Lower deviations of branching processes in random environment with geometrical offspring distributions[END_REF]. In the particular case when the offspring distribution is geometric, precise asymptotics can be found in Böinghoff [START_REF] Böinghoff | Limit theorems for strongly and intermediately supercritical branching processes in random environment with linear fractional offspring distributions[END_REF], Kozlov [START_REF] Kozlov | On large deviations of branching processes in a random environment: geometric distribution of descendants[END_REF].

An important closely linked issue is the asymptotic behavior of the distribution of a BPRE (Z n ), i.e. the limit of P(Z n = j|Z 0 = k) as n → ∞, for fixed j 1 when the process starts with k 1 initial individuals. For the Galton-Watson process, the asymptotic is well-known and can be found in the book by Athreya [START_REF] Athreya | Branching processes[END_REF]. For the need of the lower large deviation principle of a BPRE, Bansaye and Böinghoff have shown in [START_REF] Bansaye | Small positive values for supercritical branching processes in random environment[END_REF] that, for any fixed j 1 and k 1 it holds n -1 log P(Z n = j|Z 0 = k) → -ρ as n → ∞, where ρ > 0 is a constant. This result characterizes the exponential decrease of the probability P(Z n = j|Z 0 = k) for the general supercritical case, when extinction can occur. However, it stands only on a logarithmic scale, and the constant ρ is not explicit, except when the reproduction law is fractional linear, for which ρ is explicitly computed in [START_REF] Bansaye | Small positive values for supercritical branching processes in random environment[END_REF]. Sharper asymptotic results for the fractional linear case can be found in [START_REF] Böinghoff | Limit theorems for strongly and intermediately supercritical branching processes in random environment with linear fractional offspring distributions[END_REF]. In the present paper, we improve the results of [START_REF] Bansaye | Small positive values for supercritical branching processes in random environment[END_REF] and extend those of [START_REF] Böinghoff | Limit theorems for strongly and intermediately supercritical branching processes in random environment with linear fractional offspring distributions[END_REF] by giving an equivalent of the probability P(Z n = j|Z 0 = k) as n → ∞, provided that each individual gives birth to at least one child. These results are important to understand the asymptotic law of the process, and are useful to obtain sharper asymptotic large deviation results. We also improve the result of [START_REF] Huang | Moments, moderate and large deviations for a branching process in a random environment[END_REF] about the critical value for the harmonic moment of the limit variable W = lim n→∞ Zn E(Zn|ξ) . Let us explain briefly the findings of the paper. Assume that P(Z 1 = 0) = 0. From Theorem 2.3 of the paper it follows that when Z 0 = 1, (1.1)

P (Z n = j) ∼ n→∞ γ n q j with γ = P(Z 1 = 1) > 0,
where q j ∈ (0, +∞) can be computed as the unique solution of some recurrence equations; moreover, the generating function Q(t) = ∞ j=1 q j t j has the radius of convergence equal to 1 and is characterized by the functional equation

(1.2) γQ(t) = EQ(f 0 (t)), t ∈ [0, 1),
where f 0 (t) = ∞ i=1 p i (ξ 0 )t i is the conditional generating function of Z 1 given the environment. These results extend the corresponding results for the Galton-Watson process (see [START_REF] Athreya | Branching processes[END_REF]). They also improve and complete the results in [START_REF] Bansaye | Small positive values for supercritical branching processes in random environment[END_REF] and [START_REF] Böinghoff | Limit theorems for strongly and intermediately supercritical branching processes in random environment with linear fractional offspring distributions[END_REF]: it was proved in [START_REF] Bansaye | Small positive values for supercritical branching processes in random environment[END_REF] that 1 n log P (Z n = j) → log γ, and in [START_REF] Böinghoff | Limit theorems for strongly and intermediately supercritical branching processes in random environment with linear fractional offspring distributions[END_REF] that P (Z n = 1) ∼ n→∞ γ n q 1 in the fractional linear case.

In the proofs of the above results we make use of Theorem 2.1 which shows that, with m 0 = E ξ Z 1 , we have, for any fixed a > 0,

(1.3) EW -a < ∞ if and only if E [p 1 (ξ 0 )m a 0 ] < 1, under a simple moment condition E [m p
0 ] < ∞ for some p > a, which is much weaker than the boundedness condition used in [START_REF] Huang | Moments, moderate and large deviations for a branching process in a random environment[END_REF]Theorem 1.4] (see (2.10) below).

For the proof of Theorem 2.1 our argument consists of two steps. In the first step we prove the existence of the harmonic moment of some order a > 0 using the functional relation (3.5). The key argument to approach the critical value is in the second step, which is based on the method developed in [START_REF] Liu | Asymptotic properties of supercritical age-dependent branching processes and homogeneous branching random walks[END_REF]Lemma 4.1] obtaining the decay rate of the Laplace transform φ(t) = Ee -tW as t → ∞, starting from a functional inequation of the form

(1.4) φ(t) qEφ(Y t) + Ct -a ,
where Y is a positive random variable. To prove (1.4) we use a recursive procedure for branching processes starting with k individuals and choosing k large enough. The intuition behind this consideration is that as the number of starting individuals k becomes larger, the decay rate of φ k (t) = E e -tW |Z 0 = k as t → ∞ is higher which leads to the desired functional inequation.

In the proof of Theorem 2.3, the equivalence relation (1.1) and the recursive equations for the limit values (q j ) come from simple monotonicity arguments. The difficulty is to characterize the sequence (q j ) by its generating function Q. To this end, we first calculate the radius of convergence of Q by determining the asymptotic behavior of the normalized harmonic moments EZ -r n /γ n as n → ∞ for some r > 0 large enough and by using the fact that ∞ j=1 j -r q j = lim n→∞ EZ -r n /γ n . We then show that the functional equation (1.2) has a unique solution subject to an initial condition.

The rest of the paper is organized as follows. The main results, Theorems 2.1 and 2.3, are presented in Section 2. Their proofs are given in Sections 3 and 4.

Main results

A BPRE (Z n ) can be described as follows. The random environment is represented by a sequence ξ = (ξ 0 , ξ 1 , ...) of independent and identically distributed random variables (i.i.d. r.v.'s), whose realizations determine the probability generating functions (2.1)

f n (t) = f (ξ n , t) = ∞ i=0 p i (ξ n )t i , t ∈ [0, 1], p i (ξ n ) 0, ∞ i=0 p i (ξ n ) = 1.
The branching process (Z n ) n 0 is defined by the relations

(2.2) Z 0 = 1, Z n+1 = Zn i=1 N n,i , for n 0,
where N n,i is the number of children of the i-th individual of the generation n.

Conditionally on the environment ξ, the r.v.'s N n,i (i = 1, 2, ...) are independent of each other with common probability generating function f n , and also independent of Z n .

In the sequel we denote by P ξ the quenched law, i.e. the conditional probability when the environment ξ is given, and by τ the law of the environment ξ. Then P(dx, dξ) = P ξ (dx)τ (dξ) is the total law of the process, called annealed law. The corresponding quenched and annealed expectations are denoted respectively by E ξ and E. We also denote by P k and E k the corresponding probability and expectation starting with k individuals. For n ∈ N, the probability generating function of Z n is

(2.3) G n (t) = Et Zn = E [f 0 • . . . • f n-1 (t)] = E [g n (t)] ,
where

g n (t) = f 0 • . . . • f n-1 (t)
(2.4) G k,n (t) = E k t Zn = E g k n (t)
. We also define, for n 0,

m n = m(ξ n ) = ∞ i=0 ip i (ξ n ) and Π n = E ξ Z n = m 0 ...m n-1 ,
where m n represents the average number of children of an individual of generation n when the environment ξ is given. Let (2.5)

W n = Z n Π n , n 0,
be the normalized population size. It is well known that under P ξ , as well as under P, the sequence (W n ) n 0 is a non-negative martingale with respect to the filtration

F n = σ (ξ, N j,i , 0 j n -1, i = 1, 2 . . .) ,
where by convention F 0 = σ(ξ). Then the limit W = lim n→∞ W n exists P -a.s. and EW 1. We shall assume that µ := E log m 0 ∈ (0, ∞), which implies that the BPRE is supercritical and that (2.6)

γ := P(Z 1 = 1) ∈ [0, 1).
With the extra condition E| log(1-p 0 (ξ 0 ))| < ∞ (see [START_REF] Smith | On branching processes in random environments[END_REF]), the population size tends to infinity with positive probability. We also assume in the whole paper that each individual gives birth to at least one child, i.e.

(2.7) p 0 (ξ 0 ) = 0 a.s.

Therefore, under the condition

(2.8) E Z 1 m 0 log + Z 1 < ∞,
the martingale (W n ) converges to W in L 1 (P) (see e.g. [START_REF] Tanny | A necessary and sufficient condition for a branching process in a random environment to grow like the product of its means[END_REF]) and

P(W > 0) = P(Z n → ∞) = 1.
Our first result concerns the harmonic moments of the r.v. W .

Theorem 2.1. Assume that there exists a constant p > 0 such that E [m p 0 ] < ∞. Then for any a ∈ (0, p),

E k W -a < ∞ if and only if E p k 1 (ξ 0 )m a 0 < 1.
From Theorem 2.1 we get the following corollary.

Corollary 2.2. Let a k > 0 be the solution of the equation

(2.9) E[p k 1 m a k 0 ] = 1. Assume that Em a k 0 < ∞. Then, E k W -a < ∞ for a ∈ [0, a k ), E k W -a = ∞ for a ∈ [a k , ∞).
The solution a k of the equation (2.9) is the critical value for the existence of harmonic moments of the r.v. W . Note that, when the process starts with one individual, the critical value a 1 for the harmonic moments of W has been found in Theorem 1.4 of [START_REF] Huang | Moments, moderate and large deviations for a branching process in a random environment[END_REF] under the more restrictive condition (2.10)

A 1 m 0 and ∞ i=1 i 1+δ p i (ξ 0 ) A 1+δ a.s.,
where δ > 0 and 1 < A 1 < A are some constants. Theorem 2.1 and Corollary 2.2 generalize the result of [START_REF] Huang | Moments, moderate and large deviations for a branching process in a random environment[END_REF], in the sense that we consider k initial individuals rather than just one and that the boundedness condition (2.10) is relaxed to the simple moment condition

E [m p 0 ] < ∞. The next result gives an equivalent as n → ∞ of the probability P k (Z n = j) = P (Z n = j|Z 0 = k), with k ∈ N * and j k, in the case when P(Z 1 = 1) > 0. The last condition implies that, for k 1, (2.11) γ k = P k (Z 1 = k) = E[p k 1 (ξ 0 )] > 0.
Define r k as the solution of the equation (2.12) γ k = Em -r k 0 . Theorem 2.3. Assume that P(Z 1 = 1) > 0. For any k 1 the following assertions holds. a) For any accessible state j k in the sense that P k (Z l = j) > 0 for some l 0, we have

(2.13) P k (Z n = j) ∼ n→∞ γ n k q k,j
, where q k,k = 1 and, for j > k, q k,j ∈ (0, +∞) is the solution of the recurrence relation

(2.14) γ k q k,j = j i=k p(i, j)q k,i , with q k,i = 0 for any non-accessible state i, i.e. P k (Z l = i) = 0 for all l 0. b) Assume that there exists ε > 0 such that E[m r k +ε 0 ] < ∞. Then, for any r > r k , we have ∞ j=k j -r q k,j < ∞.
In particular the radius of convergence of the power series

(2.15) Q k (t) = +∞ j=k q k,j t j is equal to 1.
c) For all t ∈ [0, 1) and k 1, we have,

(2.16) G k,n (t) γ n k ↑ Q k (t) as n → ∞,
where G k,n is the probability generating function of Z n when Z 0 = k, defined in (2.4). d) Q k (t) is the unique power series which verifies the functional equation

(2.17) γ k Q k (t) = E [Q k (f 0 (t))] , t ∈ [0, 1), with the condition Q (k) k (0) = 1.
Part a) improves the bound P (Z n j) n j γ n obtained in [START_REF] Bansaye | Large deviations for branching processes in random environment. Markov Process[END_REF] (Lemma 7) for a BPRE with P(Z 1 = 0) = 0. Furthermore, Theorem 2.3 extends the results of [START_REF] Athreya | Branching processes[END_REF] for the Galton-Watson process, with some significant differences. Indeed, when the environment is random and non-degenerate, we have, for

k 2, G k,1 (t) = Ef k 0 (t) = G k 1 (t) in general, which implies that Q k (t) = Q k (t), whereas we have the relation Q k (t) = Q k (t)
for the Galton-Watson process. Theorem 2.3 also improves the results of [START_REF] Bansaye | Small positive values for supercritical branching processes in random environment[END_REF] (Theorem 2.1), where it has been proved that for a general supercritical BPRE (2.18) lim

n→∞ 1 n log P k (Z n = j) = -ρ < 0.
Our result is sharper in the case where P (Z 1 = 0) = 0. Moreover, in the case where P (Z 1 = 0) = 0, it has been stated mistakenly in [START_REF] Bansaye | Small positive values for supercritical branching processes in random environment[END_REF] that lim n→∞

1 n log P k (Z n = j) = k log γ, whereas the correct asymptotic is lim n→∞ 1 n log P k (Z n = j) = log γ k .
Now we discuss the particular fractional linear case. The reproduction law of a BPRE is said to be fractional linear if

(2.19) p 0 (ξ 0 ) = a 0 , p k (ξ 0 ) = (1 -a 0 )(1 -b 0 ) b 0 b k 0 ,
with generating function f 0 given by

f 0 (t) = a 0 + (1 -a 0 )(1 -b 0 )t 1 -b 0 t ,
where a 0 ∈ [0, 1), b 0 ∈ (0, 1), with a 0 + b 0 1, are random variables depending on the environment ξ 0 . In this case, the mean of the offspring distribution is given by

m 0 = 1 -a 0 1 -b 0 .
The constant ρ in (2.18) was computed in [START_REF] Bansaye | Small positive values for supercritical branching processes in random environment[END_REF]: with

X = log m 0 , ρ =      -log E[e -X ] if E[Xe -X ] 0 (intermediately and strongly supercritical case), -log inf λ 0 E[e -X ] if E[Xe -X ] < 0 (weakly supercritical case).
Moreover, precise asymptotic results for the strongly and intermediately supercritical case can be found in [START_REF] Böinghoff | Limit theorems for strongly and intermediately supercritical branching processes in random environment with linear fractional offspring distributions[END_REF], where the following assertions are proved:

(1) if E[Xe -X ] > 0 (strongly supercritical case),

P(Z n = 1) ∼ ν E[e -X ] n ;
(2) if E[Xe -X ] = 0 (intermediately supercritical case),

P(Z n = 1) ∼ θ E[e -X ] n l(n)n -(1-s) ,
with θ, ν, s positive constants and l(•) a slowly varying function. In the particular case where a 0 = 0, Theorem 2.3 recovers Theorem 2.1.1 of [START_REF] Böinghoff | Limit theorems for strongly and intermediately supercritical branching processes in random environment with linear fractional offspring distributions[END_REF] with p 1 (ξ 0 ) = 1/m 0 , X = log m 0 > 0 and E Xe -X > 0. Therefore the process is strongly supercritical and

P(Z n = 1) ∼ ν E[e -X ] n = γ n .
However, since we assume P(Z 1 = 0) = 0, our result does not highlight the previous two asymptotic regimes stated in the particular case when the distribution is fractional linear. The study of the general case is a challenging problem which still remains open.

Harmonic moments of W

In this section we prove Theorem 2.1. Denote the quenched Laplace transform of W under the environment ξ by

φ k (t) = E k [φ ξ (t)] = E φ k ξ (t) = E k e -tW .
We start with a lemma which gives a lower bound for the harmonic moment of W .

Lemma 3.1. Assume that E [m p 0 ] < ∞ for some constant p > 0. For any k 1, let (3.3) α k = p 1 -log Em p 0 / log γ k , with the convention that α k = p if p 1 (ξ 0 ) = 0 a.s. (so that γ k = 0). Then, for all a ∈ (0, α k ), E k W -a < ∞. Furthermore, if P(p 1 (ξ 0 ) = 0) < 1, we have α k < α k+1 ; if additionally P (p 1 (ξ 0 ) < 1) = 1, then lim k→∞ α k = p.
Proof. We use the same approach as in [START_REF] Grama | Berry-Esseen's bound and Cramér's large deviation expansion for a supercritical branching process in a random environment[END_REF] where the case k = 1 was treated. Since W is a positive random variable, it can be easily seen that, for α > 0,

(3.4) E k W -α = 1 Γ(α) +∞ 0 φ k (t)t α-1 dt,
where Γ(α) = ∞ 0 t α-1 e -t dt is the Gamma function. Moreover, it is well-known that φ ξ (t) satisfies the functional relation

(3.5) φ ξ (t) = f 0 φ T ξ t m 0 ,
where f 0 (t) = ∞ k=1 p k (ξ 0 )t k is the generating function of Z 1 under ξ 0 , defined in (2.1). Using (3.5) and the fact that

φ k T ξ t m 0 φ 2 T ξ t m 0 for all k 2, we obtain (3.6) φ ξ (t) p 1 (ξ 0 )φ T ξ t m 0 + (1 -p 1 (ξ 0 ))φ 2 T ξ t m 0 .
Taking the k-th power in (3.6), using the binomial expansion and the fact that

φ 2k-i T ξ t m 0 φ k+1 T ξ t m 0 for all i ∈ {0, . . . , k -1}, we get φ k ξ (t) = p k 1 (ξ 0 )φ k T ξ t m 0 + k-1 i=0 C i k p 1 (ξ 0 ) i (1 -p 1 (ξ 0 )) k-i φ 2(k-i)+i T ξ t m 0 p k 1 (ξ 0 )φ k T ξ t m 0 + (1 -p k 1 (ξ 0 ))φ k+1 T ξ t m 0 = φ k T ξ t m 0 p k 1 (ξ 0 ) + (1 -p k 1 (ξ 0 ))φ T ξ t m 0 . (3.7)
By iteration, this leads to

(3.8) φ k ξ (t) φ k T n ξ t Π n n-1 j=0 p k 1 (ξ j ) + (1 -p k 1 (ξ j ))φ T n ξ t Π n .
Taking expectation and using the fact that φ T n ξ (•) 1, we have

φ k (t) E   n-1 j=0 p k 1 (ξ j ) + (1 -p k 1 (ξ j ))φ T n ξ t Π n   .
Since φ ξ (•) is non-increasing, using a truncation, we get for all A > 1,

φ k (t) E   n-1 j=0 p k 1 (ξ j ) + (1 -p k 1 (ξ j ))φ t A n   + P(Π n A n ).
As T n ξ is independent of σ(ξ 0 , ..., ξ n-1 ), and the r.v.'s p 1 (ξ i ) (i 0) are i.i.d., we obtain φ k (t)

γ k + (1 -γ k )φ t A n n + P(Π n A n ),
where [START_REF] Böinghoff | Upper large deviations of branching processes in a random environment -offspring distributions with geometrically bounded tails[END_REF]. By the dominated convergence theorem, we have lim t→∞ φ(t) = 0. Thus, for any δ ∈ (0, 1), there exists a constant K > 0 such that, for all t K, we have φ(t) δ. Consequently, for all t KA n , we have

γ k = Ep k 1 (ξ 0 ) is defined in (2.
φ t A n δ and (3.9) φ k (t) β n + P(Π n A n ),
where

(3.10) β = γ k + (1 -γ k )δ ∈ (0, 1
(t) 2β n 0 2β -1 (t/K) log β log A = C 0 t -α , with C 0 = 2β -1 K α and α = -log β log A > 0.
Thus, we can choose a constant C > 0 large enough, such that, for all t > 0, (3.12)

φ k (t) Ct -α .
Furthermore, by the definition of β, A and α, we have

α = p 1 -log Em p 0 / log (γ k + (1 -γ k )δ)
, where δ ∈ (0, 1) is an arbitrary constant and γ k = Ep k 1 (ξ 0 ). When δ → 0, we have α → α k , so that (3.12) holds for all α < α k , where α k is defined in (3.3). By (3.4) and (3.12), we conclude that EW -α < ∞ for any α < α k . Moreover, it is easily seen that if

P(p 1 (ξ 0 ) = 0) < 1, then α k < α k+1 since γ k+1 < γ k ; if additionally P(p 1 (ξ 0 ) < 1) = 1, then lim k→∞ γ k = 0 so that lim k→∞ α k = p.
The following lemma is the key technical tool to study the exact decay rate of the Laplace transform of the limit variable W .

Lemma 3.2 ([15]

, Lemma 4.1). Let φ : R + → R + be a bounded function and let Y be a positive random variable such that for some constants q ∈ (0, 1), a ∈ (0, ∞), C > 0 and t 0 0 and all t > t 0 ,

φ(t) qEφ(Y t) + Ct -a . If qE (Y -a ) < 1, then φ(t) = O(t -a ) as t → ∞.
Now we proceed to prove Theorem 2.1. We first prove the necessity. Assume that E k W -a < ∞ for some a > 0. We shall show that Ep k 1 (ξ 0 )m a 0 < 1. Note that the r.v. W admits the well-known decomposition

W = 1 m 0 Z 1 i=1 W (i),
where the r.v.'s W (i) (i 1) are i.i.d. and independant of Z 1 under P ξ , and are also independent of Z 1 and ξ 0 under P. The conditional probability law of W (i) satisfies

P ξ (W (i) ∈ •) = P T ξ (W ∈ •). Since P k (Z 1 k + 1) > 0, we have (3.13) E k W -a > E k m a 0   Z 1 i=1 W (i)   -a 1{Z 1 = k} = Ep k 1 (ξ 0 )m a 0 E k W -a ,
which implies that Ep k 1 (ξ 0 )m a 0 < 1. We now prove the sufficiency. Assume that Em p 0 < ∞ and Ep k 1 (ξ 0 )m a 0 < 1 for some a ∈ (0, p).

We first consider the case where P(p 1 (ξ 0 ) < 1) = 1. We prove that E k W -a < ∞ by showing that φ k (t) = O t -(a+ε) as t → ∞, for some ε > 0. By Lemma 3.1, there exists an integer j k large enough and a constant C > 0 such that (3.14) φ j (t) Ct -(a+ε) , with ε > 0 and a + ε < p. By (3.7), we have

(3.15) φ j-1 ξ (t) p j-1 1 (ξ 0 )φ j-1 T ξ t m 0 + φ j T ξ t m 0 .
Taking the expectation in (3.15), using (3.2), (3.14) and the independence between ξ 0 and T ξ, we obtain

φ j-1 (t) E p j-1 1 (ξ 0 )φ j-1 t m 0 + Ct -(a+ε) = γ j-1 E [φ j-1 (Y t)] + Ct -(a+ε) , (3.16)
where γ j-1 = E p j-1 1 (ξ 0 ) < 1 and Y is a positive random variable whose distribution is determined by

E [g(Y )] = 1 γ j-1 E p j-1 1 (ξ 0 )g 1 m 0 ,
for all bounded and measurable function g. By hypothesis, Ep k 1 (ξ 0 )m a 0 < 1. Then, by the dominated convergence theorem, there exists ε > 0 small enough such that Ep k 1 (ξ 0 )m a+ε 0 < 1, and since j -1 k, we have Ep j-1 a+ε) ] < 1 and using (3.16) and Lemma 3.2, we get φ j-1 (t) = O(t -(a+ε) ) as t → ∞. By induction, applying (3.15) and (3.16) to the functions φ j-2 , φ j-3 , . . . , φ k and using the same argument as in the proof for φ j-1 , we obtain

1 (ξ 0 )m a+ε 0 Ep k 1 (ξ 0 )m a+ε 0 < 1. Therefore, γ j-1 E[Y -(
(3.17) φ k (t) = O(t -(a+ε) ) as t → ∞.
Therefore, in the case where P(p 1 (ξ 0 ) < 1) = 1, we have proved that

(3.18) Ep k 1 (ξ 0 )m a 0 < 1 implies E k W -a < ∞.
Now consider the general case where P(p 1 (ξ 0 ) < 1) < 1. Denote the distribution of ξ 0 by τ 0 and define a new distribution τ0 as

(3.19) τ0 (•) = τ 0 (•|p 1 (ξ 0 ) < 1).
Consider the new branching process whose environment distribution is τ = τ0 ⊗N instead of τ = τ ⊗N 0 . The corresponding probability and expectation are denoted by P(dx, dξ) = P ξ (dx)τ (dξ) and Ẽ, respectively. Of course (W n ) is still a martingale under P. Moreover, the condition Ẽ

Z 1 m 0 log + Z 1 = E Z 1 m 0 log + Z 1 /P(p 1 (ξ 0 ) < 1) < ∞ implies that W n → W in L 1 ( P). Now we show that E k [W -a ] Ẽk [W -a ]. For 0 i n, denote A i,n = (ξ 0 , . . . , ξ n-1 ) | p 1 (ξ j 1 ) = . . . = p 1 (ξ j i ) = 1 for some 0 j 1 < . . . < j i n -1,
and p 1 (ξ h ) < 1 for all h ∈ {0, . . . n -1}\{j 1 , . . . , j i } .

Conditioning by the events A i,n (i ∈ {0, . . . , n}) and using the fact that the r.v.'s ξ 0 , . . . , ξ n-1 are i.i.d., we obtain, for all n ∈ N,

E k W -a n = n i=0 E k W -a n A i,n P (A i,n ) = n i=0 E k W -a n A i,n C i n η i (1 -η) n-i , (3.20)
with η = P(p 1 (ξ 0 ) = 1). Moreover, using (2.2), a straightforward computation leads to the decomposition

(3.21) W n = n-1 i=0 η i , with n 1 and η i = 1 Z i Z i j=1 N i,j m i .
Note that, on the event {p 1 (ξ i ) = 1} we have η i = 1. Therefore, using (3.21) and the fact that the r.v.'s ξ 0 , . . . , ξ n-1 are i.i.d., we get [START_REF] Huang | Moments, moderate and large deviations for a branching process in a random environment[END_REF] Lemma 2.1). Thus, by (3.20) and (3.22), we obtain

(3.22) E k W -a n A i,n = Ẽk W -a n-i . By the convexity of the function x → x -a , we have sup n i Ẽk [W -a n-i ] Ẽk W -a (see
(3.23) E k W -a Ẽk W -a .
Note that, conditioning by the events {p 1 (ξ 0 ) = 1} and {p 1 (ξ 0 ) < 1}, we have

Ep k 1 (ξ 0 )m a 0 = (1 -η) Ẽp k 1 (ξ 0 )m a 0 + η, with η = P(p 1 (ξ 0 ) = 1). So the condition Ep k 1 (ξ 0 )m a 0 < 1 implies that Ẽp k 1 (ξ 0 )m a 0 < 1.
Then, by (3.18) applied under the probability P, and the fact that P(p 1 (ξ 0 ) < 1) = 1, we get Ẽk [W -a ] < ∞. Therefore, by (3.23), it follows that

(3.24) Ep k 1 (ξ 0 )m a 0 < 1 implies E k W -a < ∞, which ends the proof of Theorem 2.1.

Small value probability in the non-extinction case

In this section we prove Theorem 2.3. We start with the proof of part a). For k 1 and j k, define

(4.1) a k,n (j) = P (Z n = j|Z 0 = k) γ n k , with γ k = P k (Z 1 = k).
By the Markov property, we have

P k (Z n+1 = j) P k (Z 1 = k) P k (Z n = j) . Dividing by γ n+1 k leads to (4.2) a k,n+1 (j) a k,n (j).
Therefore, by the monotone ratio theorem, we obtain

lim n→∞ ↑ a k,n (j) = q k,j ∈ R.
We shall prove that q k,j satisfies the properties claimed in the theorem. If j is such that P k (Z n = j) = 0 for any n 0, then a k,n (j) = 0 for any n 0, so that lim n→∞ a k,n (j) = 0 = q k,j . If there exists l 0 such that P k (Z l = j) > 0, then q k,j a k,l (j) = P k (Z l = j)/γ l k > 0. Now we show by induction that for all j k, we have

H(j) : sup n∈N a k,n (j) = a k (j) < ∞.
For j = k, we have a k (k) = 1. Assume that j k + 1 and that H(i) is true for all k i j -1. By the total probability formula, we obtain

P k (Z n+1 = j) γ n+1 k = 1 γ k j i=k P k (Z n+1 = j|Z n = i) P k (Z n = i) γ n k , which is equivalent to (4.3) a k,n+1 (j) = 1 γ k   j-1 i=k p(i, j)a k,n (i) + γ j a k,n (j)   with p(i, j) = P (Z 1 = j|Z 0 = i).
Using the fact that a k,n (j) a k,n+1 (j), we get by induction that

sup n∈N a k,n+1 (j)(γ k -γ j ) j-1 i=k p(i, j)a k (i) < ∞.
Thus q k,j < ∞ for all j k + 1 and k 1. Furthermore, taking the limit as n → ∞ in (4.3), leads to the following recurrent relation for q k,j :

q k,k = 1, γ k q k,j = j i=k p(i, j)q k,i (j k + 1).
This end the proof of part a) of Theorem 2.3. Now we prove part b) of Theorem 2.3. We give a proof that the radius of convergence of the power series Q k is equal to 1. The method is new even in the case of the Galton-Watson process. We start with a lemma. Lemma 4.1. Let k 1. Assume that P(p 1 (ξ 0 ) < 1) = 1 and that there exists ε > 0 such that E[m r k +ε 0 ] < ∞, where r k is the solution of the equation γ k = Em -r k 0 . Then, for any r > r k , we have

(4.4) lim n→∞ ↑ E k Z -r n γ n k < ∞.
Proof. By the Markov property,

E k Z -r n+1 E k Z -r n+1 |Z 1 = k P k (Z 1 = k) = γ k E k Z -r n ,
which proves that the sequence (E k [Z -r n ] /γ n k ) n∈N is increasing. We show that it is bounded. For n 1 and m 0, we have the following well-known branching property for Z n :

(4.5)

Z n+m = Zm i=1 Z (m) n,i ,
where, under P ξ , the random variables Z (m) n,i (i 1) are i.i.d., independent of Z m , whose conditional probability law satisfies P ξ Z (m) n,i ∈ • = P T m ξ (Z n ∈ •), with T m the shift operator defined by T m (ξ 0 , ξ 1 , . . .) = (ξ m , ξ m+1 , . . .). Intuitively, relation (4.5) shows that, conditionally on Z m = i, the annealed law of the process Z n+m is the same as that of a new process Z n starting with i individuals.

Using (4.5) with m = 1, the independence between Z 1 and Z

(1)
n,i (i 1) and the fact that

E i Z -r n E k+1 Z -r
n for all i k + 1, we have

E k Z -r n+1 = E k Z -r n+1 |Z 1 = k P k (Z 1 = k) + ∞ i=k+1 E   i h=1 Z (1) n,h -r Z 1 = i   P k (Z 1 = i) γ k E k Z -r n + E k+1 Z -r n . (4.6)
We shall use the following change of measure: for k 1 and r > 0, let P

E k+1 Z -r n = E (r) k+1 [W -r n ]c n r , with sup n∈N E (r) k+1 [W -r n ] = E (r)
k+1 [W -r ] (see [START_REF] Huang | Moments, moderate and large deviations for a branching process in a random environment[END_REF], Lemma 2.1). Moreover, we have E (r) [p k+1

1 (ξ 0 )m r 0 ] = γ k+1 /Em -r 0 < 1 for any r < r k+1 . So by Theorem 2.1 we get

E (r) k+1 [W -r ] = C(r) < ∞ and then E k+1 [Z -r
n ] C(r)c n r for any r < r k + ε < r k+1 . Coming back to (4.6) with r < r k + ε, we get by induction that (4.9)

E k Z -r n+1 γ n+1 k + C n j=0 γ n-j k c j r .
Choose r > r k such that c r < γ k . Then, we have, as n → ∞,

(4.10) E k Z -r n+1 γ n+1 k 1 + C γ k n j=0 c r γ k j → C γ k -c r . Thus the sequence (E k [Z -r n ] /γ n k ) n∈N is bounded and (4.4) holds for any r ∈ (r k , r k + ε). Using the fact that E k Z -r n+1 E k Z -r
n+1 for any r > r, the result follows for any r > r k , which ends the proof of the lemma. Remark 4.2. From the results stated above, with some additional analysis one can obtain the equivalent of the harmonic moments EZ -r n for any r > 0. However, it is delicate to have an expression of the concerned constant in the equivalence. This will be considered in a forthcoming paper. Now we show that the radius of convergence R of the power series Q k (t) = ∞ j=k q k,j t j is equal to 1. Using the fact that ∞ j=k P k (Z n = j) = 1, part a) of Theorem 2.3 and the monotone convergence theorem, we have

lim n→∞ ↑ γ -n k ∞ j=k P k (Z n = j) = ∞ j=k q k,j = +∞,
which proves that R 1. We prove that R = 1 by showing that +∞ j=k j -r q k,j < ∞ for r > 0 large enough. Using part a) of Theorem 2.3, the monotone convergence theorem and Lemma 4.1, we have, for any r > r k , (4.11)

+∞ j=k j -r q k,j = +∞ j=k j -r lim n→∞ ↑ P k (Z n = j) γ n k = lim n→∞ ↑ E k Z -r n γ n k < ∞,
which proves part b). Now we prove part c) of Theorem 2.3. Using part a), the definition of G k,n and the monotone convergence theorem, we get (2.16). To prove the functional relation (2.17), recall that G k,1 (t) = ∞ j=k p(k, j)t j = Ef k 0 (t). By (2.14) and Fubini's theorem, we get

γ k Q k (t) = ∞ j=k ∞ i=k q k,i p(i, j)1(i j)t j = ∞ i=k q k,i ∞ j=i p(i, j)t j = ∞ i=k q k,i E f i 0 (t) = E ∞ i=k q k,i f i 0 (t) = E [Q k (f 0 (t))] .
This proves the functional relation (2.17). We now prove that the previous functional relation characterizes the function Q k . To this end it suffices to show the unicity of the solution of (2.17). Assume that there exists a power series Q(t) = ∞ j=0 qk,j t j on [0, 1) which verifies (2.17) with the initial condition q k,k = qk,k = 1. We first show by induction in n that Q(n) (0) = 0 for all n ∈ {0, . . . , k -1}. Since f 0 (0) = 0 and γ k ∈ (0, 1), by (2.17), we get γ k Q(0) = Q(0), which implies that Q(0) (0) = Q(0) = 0. By the induction hypothesis we have that Q(j) (0) = 0 for all j ∈ {0, . . . , n -1} for some n < k -1. We show that Q(n) (0) = 0.

Using Faà di Bruno's formula, we have

(4.12) Q • f 0 (n) (t) = n j=1 Q(j) (f 0 (t))B n,j f (1) 0 (t), . . . , f (n-j+1) 0 (t) ,
where B n,j are the Bell polynomials, defined for any 1 j n by

B n,j (x 1 , x 2 , . . . , x n-j+1 ) = n! i 1 !i 2 ! • • • i n-j+1 ! x 1 1! i 1 x 2 2! i 2 • • • x n-j+1 (n -j + 1)! i n-j+1
, where the sum is taken over all sequences (i 1 , . . . , i n-j+1 ) of non-negative integers such that i 1 + • • • + i n-j+1 = j and i 1 + 2i 2 + • • • + (n -j + 1)i n-j+1 = n. In particular B n,n (x 1 ) = x n 1 . Applying (4.12) and using the fact that f 0 (0) = 0, B n,n f

(1) 0 (0) = f (1) 0 (0) n and Q(j) (0) = 0 for all j ∈ {0, . . . , n -1}, we get

(4.13) Q • f (n) (0) = Q(n) (0) f (1) 0 (0) n .
Then taking the derivative of order n of both sides of (2.17) and using (4.13), we obtain that γ k Q(n) (0) = Q(n) (0)γ n for n < k -1, which implies that Q(n) (0) = 0. Now we show that qk,j = q k,j for any j k + 1. Using Fubini's theorem, the fact that f 0 , . . . , f n-1 are i.i.d. and iterating (2.17), we get

Q k (t) -Qk (t) = γ -n k E Q k (ḡ n (t)) -Qk (ḡ n (t)) = γ -n k ∞ j=k (q k,j -qk,j )E ḡj n (t) γ -n k ∞ j=k+1
|q k,j -qk,j | G j,n (t), (4.15) where G j,n (t) is the generating function of Z n starting with j individuals. To conclude the proof of the unicity it is enough to show that We prove (4.16) using the Lebesgue dominated convergence theorem. Note that, by (2.16), for all n ∈ N, (4.17) γ -n j G j,n (t) Q j (t). Therefore, using the fact that γ j < γ k for all j k + 1, we have

lim n→∞ γ -n k G j,n (t) = lim n→∞ γ j γ k n γ -n j G j,n (t) lim n→∞ γ j γ k n Q j (t) = 0,
and γ -n k G j,n (t) Q j (t). Now we show that ∞ j=k+1 |q k,j -qk,j | Q j (t) < ∞ for all t ∈ [0, 1). Indeed, by part b) of Theorem 2.3, we have ∞ j=k q k,j j -r < ∞ for any r > r k . In particular for a fixed r > r k , there exists a constant C > 0 such that, for all j 1, i j, it holds q j,i Ci r . Therefore,

Q j (t) C ∞ i=j i r t i Ct j ∞ i=0 (i + j) r t i C 2 r t j ∞ i=0 i r t i + 2 r t j ∞ i=0 j r t i C r j r t j .
Since Q k and Qk are power series whose radii of convergence are equal to 1, we have, for any t < 1, (4.18)

∞ j=k+1 q k,j Q j (t) C r ∞ j=k+1
q k,j j r t j < ∞, and

∞ j=k+1 qk,j Q j (t) < ∞.
Using the dominated convergence theorem, we see that

lim n→∞ γ -n k ∞ j=k+1
|q k,j -qk,j | G j,n (t) = 0.

Therefore, from (4.15) we conclude that Q k (t) = Qk (t) for all t ∈ [0, 1). This ends the proof of Theorem 2.3.

(3. 1 )

 1 φ ξ (t) = E ξ e -tW , and the annealed Laplace transform of W starting with k individuals by (3.2)

  F n -measurable random variable T , where c r = Em -r 0 . By (4.7), we obtain(4.8)

(4. 14 )

 14 E [Q k (ḡ n (t))] = γ n k Q k (t) and E Qk (ḡ n (t)) = γ n k Qk (t), where ḡn (t) = f n-1 • . . . • f 0 (t). By(4.14), for all t ∈ [0, 1) and n ∈ N, we have

  j -qk,j | γ -n k G j,n (t) = 0.

  is the conditional probability generating function of Z n when the environment ξ is given. It follows from (2.2) that the probability generating function G k,n of Z n starting with k individuals is