Quasineutral limit of the Euler-Poisson system for ions in a domain with boundaries II

In this paper, we study the quasineutral limit of the isothermal Euler-Poisson equation for ions, in a domain with boundary. This is a follow-up to our previous work \cite{GVHKR}, devoted to no-penetration as well as subsonic outflow boundary conditions. We focus here on the case of supersonic outflow velocities. The structure of the boundary layers and the stabilization mechanism are different.


Introduction
This work is about the quasineutral limit of the isothermal Euler-Poisson equation, describing the dynamics of ions in a plasma. We focus on the role of the boundary of the plasma domain, and the associated boundary layer. It complements our previous work on the topic [5] (see also [4], [10]). We refer to the introduction of [5] for a substantial bibliography on the quasineutral limit and related issues.
In the Euler-Poisson model, ions are described by their density n 0 and their velocity field u ∈ R 3 , while electrons are assumed to follow a Maxwell-Boltzmann law, i.e., their density is given by e −φ , where φ denotes the electric potential. We consider that the plasma is contained in the domain R 3 + := {x = (y, z) ∈ R 2 × R + } (yet, general domains can actually be considered, see [5,Section 5.1]). The isothermal Euler-Poisson equation under study then reads, for (t, x) ∈ R + × R 3 + , ∂ t n + div(nu) = 0, ∂ t u + u · ∇u + T i ∇ ln(n) = ∇φ, where ε > 0 is loosely speaking the ratio between the Debye length of the plasma and the typical length of observation. In all practical situations, it is a small parameter. The parameter T i > 0 is the temperature of the ions. We consider in this work the case of supersonic outflow velocities. It means that we consider initial velocity fields u(0) = (u 1 (0), u 2 (0), u 3 (0)) satisfying We are interested in the behaviour of the solutions of (1.1) as ε goes to 0, that is, in the so-called quasineutral limit. We refer once again to the introduction of [5] for an extensive discussion. The expected formal limit, obtained by taking directly ε = 0, is the isothermal Euler system: (1.2) ∂ t n + div(nu) = 0, ∂ t u + u · ∇u + (T i + 1) ∇ ln(n) = 0, together with the neutrality relation n = e −φ . With regards to this last relation and the Dirichlet condition (1.1d) on φ, one would like to impose the boundary condition n| x3=0 = e −φ b . However, this condition is not a priori compatible with the hyperbolic system (1.2), therefore we expect the formation of a boundary layer in the vicinity of the boundary {x 3 = 0}.
In our previous work [5], we have considered subsonic outflow velocities. In this setting, a boundary condition on the normal velocity has to be added: u · n = u.
In [5], we have notably focused on the non-penetration boundary condition (that is, u = 0), and we have also considered the subsonic case − √ T i < u < 0.
The supersonic condition (1.3) is usually called Bohm condition (or criterion), while the boundary layer is usually referred to as the sheath in the physics literature. It plays an important role in the study of confined plasmas: see for instance the book of Lieberman and Lichtenberg [7,Chapter 6] and the review paper of Riemann [9].  Furthermore, the rate of convergence is O( √ ε).
Note that the supersonic condition forbids y → u 0 3 (t, y, 0) to vanish at infinity: more precisely, it behaves asymptotically in y like w ref , which necessarily satisfies w ref < − √ T i + 1. The smallness condition (1.4) is used in particular to ensure the stability of a boundary layer. The study of this boundary layer, which we build beforehand in Theorem 2, is crucial in the proof of Theorem 1, though it does not appear explicitly. We refer to Theorem 3 and Corollary 1 for more precise statements.
To prove Theorem 1 (and its refined version), we shall use a classical two-step argument as in our previous work [5]. The first step is a consistency step where we first build approximate solutions for (1.1) at a sufficiently high order, see Theorem 2. In the second step we combine linear estimates and a continuation argument to deduce the stability of those approximate solutions. The main differences compared to the analysis of [5] are as follows.
(1) Contrary to [5], the existence of the sequence of approximate solutions is not unconditional, and relies on the first part of the smallness assumption (1.4). Furthermore, at the main order, compared to the subsonic case treated in [5], there is a boundary layer for the third (i.e., the normal) component of the velocity.
(2) Because of this additional singular term, the linear estimates of [5] are not relevant. Instead, we shall consider weighted estimates, inspired by the classical work of Goodman [6]. The idea is to use the stabilizing effect of convection. In order to obtain relevant energy estimates, we shall also borrow some ideas from a recent paper of Nishibata, Ohnawa and Suzuki [8] on a related problem. Namely, this problem is the stability of a special solution of the unscaled Euler-Poisson system, that is, when ε = 1. Roughly, this special solution corresponds to the main order part of the boundary layer constructed in Theorem 2, but without any y-or ε-dependence. On that topic, one can also refer to the previous work of Suzuki [11], as well as the works of Ambroso, Méhats and Raviart [2], and Ambroso [1]).
Our analysis will combine three types of L 2 weighted estimates. In [5], only the "physical" unweighted energy estimate was needed. This was due to the fact that the physical energy, which is well adapted to symmetrize the singular term coming from the Poisson equation in the quasineutral limit, was compatible with the structure of the boundary layers. Here this is not the case and we have to use the stabilizing effect of convection via the weights and other energy functionals in order to control the new singular terms that arise.
T i can be handled. In this regime, no boundary condition can be imposed when ε > 0, but one boundary condition can (and actually must) be imposed when ε = 0, that is, for system (1.2). We can notably endow (1.2) with the boundary condition n| x3=0 = e −φ b . This additional boundary condition allows to satisfy the Dirichlet condition φ| x3=0 = φ b in the limit ε → 0, hence no boundary layer appears. Convergence of solutions of (1.1) to the solution of (1.2) is in this context straightforward to prove: one can build an accurate approximate solution of (1.1), whose first term is the solution of (1.2). As explained before, this approximate solution will not have any boundary layer part. Then, one can perform an energy estimate between the approximate and exact solutions of (1.1) (with the same initial data) to show convergence. In this way, we get Furthermore, the rate of convergence is O(ε).
Note that we get convergence in L ∞ due to the absence of boundary layers. In the setting of Theorem 1, boundary layers have to be included in order to describe the asymptotic behavior in L ∞ (see the refined version of the convergence in Corollary 1).

1.3.
Overview and classification. -Summarizing the results obtained in [5] and in the present paper, we obtain the following classification of outflow boundary conditions for the study of the quasineutral limit of the Euler-Poisson equation.
and Potential and Velocity (under smallness (non-characteristic) assumption (1.4)) The rest of the paper is now entirely devoted to the proof of Theorem 1 (more precisely of the refined Theorem 3).

Derivation of the boundary layers
We construct in this section accurate approximate solutions of the Euler-Poisson system, of boundary layer type. They are expansions in powers of ε, of the form: with K an arbitrarily large integer. These approximations split into two parts: -a regular part, with coefficients (n i , u i , φ i ) depending on (t, x). This regular part models the macroscopic behaviour of the solutions.
-a singular part, with coefficients (N i , U i , Φ i ) depending on the regular variables t, y, but also on a rescaled variable z = x 3 /ε ∈ R + . It models a boundary layer correction, of size ε near the boundary. Accordingly, we shall impose In order for the whole approximation φ a to satisfy the proper Dirichlet condition, we shall further impose Note that, far away from the boundary, the term (n 0 , u 0 , φ 0 ) will drive the dynamics of these approximate solutions, and will satisfy the quasineutral limit system: together with the relation n 0 = e −φ 0 . Our construction is summarized in the Theorem 2. -There exists δ 0 > 0 such that: for any K ∈ N * , for any (n 0 0 , u 0 0 ) such that inf n 0 0 > 0 satisfying the regularity assumption sup u 0 0,3 < 0, (sup u 0 0,3 )) 2 > T i + 1, and the smallness condition one can find T > 0 and an expansion of type (2.1) with the following properties.
) and decays together with its derivatives uniformly exponentially in z.
(iv) Let us consider a solution (n ε , u ε , φ ε ) to (1.1) and define: Then (n, u, φ) satisfies the system of equations: ∂ t u + (u a + u) · ∇u + u · ∇u a + T i ∇n n a + n − ∇n a n a n n a + n = ∇φ + ε K R u , where R n , R u , R φ are remainders satisfying: The whole section is devoted to the proof of this theorem. As already said in the introduction, the main difference with the boundary layer problem analyzed in [5] is that the third component of U 0 does not vanish: the velocity field has a boundary layer part of amplitude O(1). This modifies all the boundary layer equations, compared to [5]. In particular, well-posedness is not unconditional anymore, which explains our condition (2.7). In fact, the construction of the first boundary layer term only requires a lower bound (first inequality in (2.7)), whereas the construction of the next terms requires an additional upper bound. Moreover, the constant δ 0 in this condition can be made explicit: see Section 2.2. Let us finally point out that the regularity assumption on the initial data (n 0 0 , u 0 0 ) can be easily lowered to H s , with s large enough. (a) Away from the boundary, we find as expected that (n 0 , u 0 ) satisfies (2.4), plus the neutrality relation φ 0 = − ln(n 0 ). The local existence and uniqueness of (n 0 , u 0 ), starting from our supersonic data, will be stated rigorously in Section 2.2.
(b) In the boundary layer, with a Taylor expansion of the regular part of (2.1), the third line of (1.1) yields whereas the first line yields Taking into account (2.2), we obtain We then consider the vertical component of the momentum equation in (1.1). This gives We can integrate in z, and use (2.2) again to obtain Inserting relation (2.11) in the last identity, we end up with This last equation can be written (2.13) Φ 0 (t, y, z) = F t, y, Γn 0 (t, y) + N 0 (t, y, z) Γn 0 (t, y) where (2.14) The derivative with respect to N is One has ∂ N F (t, y, N ) < 0 over (0, N F (t, y)), N F (t, y) := Γu 0 3 (t, y) 2 /T i . Hence, for all t, y, N → F (t, y, N ) has a smooth inverse, and its inverse is decreasing from (Φ F (t, y), +∞) to (0, N F (t, y)), where Φ F (t, y) := F (t, y, N F (t, y) < 0. We make a slight abuse of notation, and denote F −1 (t, y, Φ) such an inverse.
Back to (2.10), we obtain This is an autonomous ODE in z, (t, y) playing the role of parameters. Of course, the r.h.s. is only defined as long as F −1 is. The ODE is completed by boundary conditions:  This parameter δ measures how far the quasineutral solution is from satisfying the Dirichlet condition. It will determine the size of the boundary layer.
(c) From there, one can derive (n 1 , u 1 , φ 1 ), and (N 1 , U 1 , Φ 1 ). More generally, knowing (n k , u k , φ k ), and (N k , U k , Φ k ), k i − 1, one can derive equations on (n i , u i , φ i ), and then on (N i , U i , Φ i ). We shall stick to i = 1 for the sake of brevity, and refer to [5] for a full derivation in a close context. Clearly, The well-posedness of this system over (0, T ) (where T is the time up to which the supersonic boundary condition on u 0 is satisfied) will be reminded in the next paragraph.
As regards the boundary layer terms, we look as before at the Poisson equation of (1.1), which leads to . Then, the mass equation yields (remember that U 0 Then, we write down ε 0 terms in the vertical component of the momentum equation: One can as before integrate with respect to z, to get Using (2.11) and (2.20) to transform the first term at the l.h.s., we end up with is known from the previous steps of the construction. This expression allows to express Γn 1 + N 1 in terms of Φ 1 , and to go back to (2.19) to have a closed equation on Φ 1 .
A tedious but straightforward calculation shows then that y, Φ 0 (t, y, z))Φ 1 (t, y, z) + F 6 (t, y, z) with F 6 (t, y, z) := n 0 (t, y, 0) ∂ N F t, y, n 0 (t, y, 0) + N 0 (t, y, z) n 0 (t, y, 0) This equation is completed with the boundary condition The well-posedness of this equation will be discussed in the next section.
As soon as Φ 1 is determined, the expression of N 1 follows, and from (2.20), we get the expression of U 1 3 . The horizontal part U 1 y is obtained as for U 0 y by considering the horizontal components of the momentum equation in (1.1). We leave all details to the reader. As mentioned before, the next order terms in the expansion satisfy similar problems, and we omit their construction for the sake of brevity.

2.2.
Well-posedness. -We discuss here the well-posedness of the reduced models derived formally in the previous section.
(a) The supersonic condition (2.6) expresses that all characteristics of the Euler system (2.4) are outgoing, so that no boundary condition is necessary for solvability. More precisely, starting from the initial data satisfying (2.6)-(2.5), there exists a small time T > 0 and a smooth solution [3] ). Moreover, one can assume that the supersonic condition is satisfied globally over [0, T ]:

Under this last condition, one can solve the equations on
We remind that these equations are of the type (2.18) in the special case i = 1). These are linear hyperbolic systems with unknowns (n i , u i ), and with outgoing characteristics over [0, T ] by (2.23). Starting (for instance) from zero initial data More precisely, one can check inductively that all source terms f i n,u,φ belong to , and the same holds for φ i .
(b) We now have to address rigorously the construction of boundary layer systems. Note that the N i 's and U i 's are given in an extrinsic manner from the Φ i 's, so that we only need to focus on the latter. We begin by discussing the well-posedness of (2.15)-(2.16)-(2.17). We already noticed that t, y are just parameters, which makes it an ODE problem. This ODE problem has been addressed in [2], see also [11]. Omitting the dependence with respect to t, y, it reduces to the search of a trajectory z → Φ(z) of a second-order ODE in order for F −1 and X to be defined: see the discussion below (2.14). This problem can be solved using the hamiltonian structure of the equation. Intro- The variations of V can be determined using that V (Φ) = X (Φ) has the same sign as Y (N ) = N − e −F (N ) , which in turn has the same sign as ln N + F (N ). We refer to [11] for more details. As a result, there exists δ 0 > 0 such that: for all Φ 0 ∈ [−δ 0 , +∞), there is a unique (monotonic) solution of (2.24) connecting Φ 0 to 0. Moreover, as V (0) > 0, (0, 0) is a hyperbolic point for the 2 × 2 system associated to (2.24), which yields the exponential decay of Φ by the stable manifold theorem. More precisely, the rate of decay is given by Back to the original problem (with dependence in t, y), this analysis provides a solution of (2.15)-(2.16)-(2.17) under the first inequality in (2.7). The regularity of Φ 0 with respect to (t, y) follows from standard arguments: using (for instance) the implicit function theorem, one can show that Moreover, one can check that Φ 0 (n ref , ·) = 0. The smoothness in (t, y) and H ∞ integrability in y follow. For the sake of brevity, we leave the details to the reader.
Remark 2. -Note that the decay rate (2.25) is independent of the amplitude parameter defined in Remark 1. We thus get for Φ 0 an estimate under the form: there exists C > 0, γ 0 > 0 such that for every δ ∈ (0, δ 0 ), we have the estimate As regards the next order boundary layer terms, they all (formally) satisfy equations of the type where F i depends on the lower order profiles. Their wellposedness (in the space of smooth functions, exponentially decaying in z) is shown inductively: freezing t, y, it reduces to show the well-posedness of a linear ODE of the type Note that, up to consider Φ(z) = Φ(z)−ψχ(z), with χ ∈ C ∞ c (R + ) satisfying χ(0) = 0, we can always assume that ψ = 0.
As X (Φ) < 0 for large Φ, the well-posedness of these systems is not clear under a mere lower bound on Φ 0 . But as X (0) > 0, up to take a smaller δ 0 and impose the additional upper bound in (2.7), we ensure that X (φ 0 ) α > 0 uniformly in z. The existence of a variational solution (say in H 1 0 (R + )) follows then from Lax-Milgram's lemma. Smoothness in z is a consequence of standard elliptic regularity results, whereas the exponential decrease follows again from the stable manifold theorem.
The proof of Theorem 2 is thus complete.

Linear Stability
In this section, we establish L 2 type estimates for linear systems of the following form: where r = (r n , r u , r φ ) is a given source term, and where h ∈ {h 0 , h 1 }, where cf. (4.6). We add to the system the boundary condition To prove our main linear stability estimate, we shall use Goodman type weights in order to use the stabilizing effect of convection in the problem. From the construction of the approximate solution in the previous section, in particular, from Remarks 1 and 2, we know that the main boundary layer part of the approximate solution satisfies an estimate under the form where δ will be assumed small whereas γ 0 is fixed. We recall that, as in the previous section, for any function f = f (t, x), we denote by Γf the function (t, y, z) → f (t, y, 0). We shall also assume that the trace on the boundary of the tangential velocity of the limit system is small, i.e., where µ, γ > 0 are some fixed parameters, to be specified later. Observe that η satisfies the following properties: The main idea is that the parameter µ > 0, µ γ 0 will be chosen small enough. We also assume that the parameter δ that measures the strength of the boundary layers is sufficiently small such that δ/µ 2 is also small. This yields in particular For example, we can choose µ = δ 1/4 . The assumptions on µ further imply This inequality will be crucial in many estimates of the paper. We also note for later purposes that We shall eventually assume that the forcing terms (r n , r u , r φ ) can be split into a small singular part and a regular part: Note that since the derivative of η is of order 1/ε, the first term in the r.h.s. of (3.9) is singular in ε. Concretely, the linearized system (3.1) will be obtained by taking ε-derivatives of (2.8). The remainders will contain commutators, and satisfy (3.9).
The crucial weighted L 2 estimate is given in the following proposition.
-Let (n a , u a , φ a ) be the approximate solution constructed in Theorem 2 and consider some smooth (n, u, φ) on [0, T ] such that for some M > 0 and such that the Bohm condition is verified on the boundary Moreover, let us assume that the source term (r n , r u , r φ ) of (3.1) verifies the assumption (3.9). Then, there exists δ 0 > 0, µ 0 > 0 and C(C a , M ) (C a only depends on the approximate solution) such that for every ε ∈ (0, ε 0 ], for every δ ∈ (0, δ 0 ] and for every µ ∈ (0, µ 0 ] with δ/µ 2 sufficiently small, we have for the solution (ṅ,u,φ) of (3.1) the estimate . This whole section is dedicated to the proof of Proposition 2.
Proof of Proposition 2. -We first gather some useful bounds satisfied by the approximate solution (n a , u a , φ a ). In the proof, we shall denote by C a and C(C a , M ) numbers which depend only on the estimates of the approximate solution and on the number M defined in (3.10) and that may change from line to line. The important thing is that they are uniformly bounded for ε ∈ (0, 1], δ ∈ (0, 1) and T ∈ (0, T 0 ] where T 0 is the interval of time on which the approximate solution is defined. We first have, by construction (see Theorem 2): with C a independent of ε. In the other hand, we have for all x 3 > 0: where we recall that 0 < δ 1 can be considered as a small parameter. We shall combine many energy estimates for the proof of Proposition 2. As already pointed out in the introduction, our computations share features with those led in [8, Section 2.1]. The starting point of all the energy estimates will be the following lemma: -Under the assumptions of Proposition 2, we have the estimate where the quadratic form Q 0 is associated to the symmetric matrix and (3.14) The main difficulty in proving Proposition 2 will be to handle the term I , that involves the potentialφ. We note that in the left hand side of the above estimate, the quadratic form Q 0 is positive thanks to the Bohm condition (3.11).
Proof of Lemma 1. -First, multiplying the velocity equation by (n a + n)u η, and performing standard manipulations, we obtain: where We recall that I was defined in (3.14). The last two terms at the r.h.s. of (3.16) can be easily estimated through (3.12), (3.10), (3.6) and (3.9). We find . Integrating by parts, we have the identity To write the boundary term, we have used that η(0) = 1, see (3.5). We can then use the evolution equation onṅ to express divu in terms ofṅ. We find For the last term J 2 , we use again (3.12), (3.10), (3.6) and (3.9) to get where we have used (3.7) to go from the first to the second line. We insert (3.18)- (3.19) in (3.17) to obtain (3.20) 1 (2) Treatment of I 2 .
-We write Relying once again on (3.3)-(3.12), and (3.7), we infer that: Conclusion. -Gathering (3.16), (3.20) and (3.21), we obtain To conclude, it suffices to observe that the only significant contribution in the terms involving η comes from the traces of the coefficients. Indeed, we have that to obtain Since εη and x 3 η are uniformly bounded (see (3.6)), we end up with By the Bohm condition, the quadratic form Q 0 is positive definite: this allows to absorb the last term at the r.h.s. for µ and δ small enough. The estimate of Lemma 1 follows.
A. The estimate for the weighted physical energy.
-We now use Lemma 1 in order to prove Proposition 2. The first estimate that we shall establish is the following . Note that the above estimate is not enough by itself because the singular term √ η (ε ∇ṅ, ε divu) 2 L 2 T L 2 (R 3 + ) in the right hand side is not controlled by the left hand side.
The first step in the proof of (3.25) consists in manipulating the integral I (see Lemma 1): it will give us some control on the potentialφ. More precisely, we first integrate by parts to write: There is no boundary term sinceφ satisfies an homogeneous Dirichlet condition (see (3.2)). The first term at the r.h.s can be controlled thanks to (3.3)-(3.7) and (3.10). We get ηφ (n a + n) divu. As above, we use the evolution equation onṅ to express (n a + n) divu in terms ofṅ: We integrate by parts the second integral, and obtain Note that we have used implicitly the bound to handle the last term coming from the integration by parts. (3.28) One has by using (3.12) and (3.10) the straightforward estimate and by using also (3.9) As regards the last term, we claim Lemma 2. -The following inequality holds.
We postpone the proof of the lemma to the end of the section. Combining all the previous estimates, we deduce: (2) Treatment of J 2 . -We use the Poisson equation to expressṅ in terms ofφ. We get One has (3.32) Again, the bound on the first term at the right hand side is a consequence of (3.7). Thanks to (3.12) and (3.9) one has also We point out the simplification of the boundary term, due to the homogeneous boundary condition onφ. Once again, (3.7) leads to For L 2 , by integration by parts, we can write: where we have used that ∇ 1,2φ = 0 on the boundary. The first term has a "bad sign" but will be compensated by another contribution (coming from J 3 ). The second term satisfies: Finally, the third (boundary) term of L 2 has also a bad sign, but will be compensated by the other boundary term of (3.33).
For L 3 , we rely on the fact that Together with our assumption (3.4) on the trace (|(Γu 0 ) 1,2 | δ), this implies Note that the last term in the right hand side of the above estimate is bounded by C(M ) ε∇φ 2 L 2 thanks to (3.10). Putting together the estimates on L 1 , L 2 and L 3 , we find Eventually, we derive an inequality on J 2 : (3) Treatment of J 3 . -First, we use the Poisson equation to replaceṅ: Let us start with the estimate of K 1 . With the help of an integration by parts, we get: The first term has a "good sign", and will be used to absorb the first term at the r.h.s. of (3.35). The second and third terms can be bounded using properties of the approximate solution and of the weight η (notably (3.8)). We end up with K 2 is also a good term, that will absorb the third term at the r.h.s. of (3.35). For K 3 , we use as usual (3.9) and write Finally, these last bounds yield Conclusion. -We can now collect the estimates (3.30), (3.35) and (3.37), and insert them into (3.27), followed by (3.26). We get By combining this last estimate with the estimate of Lemma 1 and using the substitution (3.23) in the terms involving η , we obtain the estimate and Q A is the positive quadratic form determined by the following symmetric matrix: Since µ can be chosen as small as we want, it suffices to prove that the leading above matrix is positive. According to Sylvester's criterion, we only have to check that the leading principal minors, denoted by (∆ i ) 1 i 5 are positive. We compute: which is also positive as a straightforward consequence of Bohm condition (3.11).
Finally the determinant is equal to: which is positive as soon as |(u a + u) 3 | > √ T i + 1 as ensured by the Bohm condition (3.11). Indeed, we have that |(u, n, φ)| = O(ε) thanks to (3.10), and that since n 0 = e φ 0 that n a = e −φa + O(δ + ε), therefore for ε and δ sufficiently small. The quadratic form Q 0 is also positive as already observed. To derive (3.25) from (3.38), it remains to note that, by positivity of Q A , we can write the Young inequality L 2 , and absorb the term µ √ η φ 2 L 2 in the left hand side by choosing µ small enough. In the same spirit, the term µ √ η ε∇φ 2 L 2 can be absorbed by the quadratic terms in ε∇φ at the left hand side of (3.38) for µ small enough.
To conclude this section, we still have to provide the proof of Lemma 2.
Proof of Lemma 2. -To estimate the term ε 2 η φ ∂ t ∂ x3φ in (3.28), we shall rely on an elliptic estimate on the Poisson equation. We first have the straightforward bound: Differentiating the Poisson with respect to time and taking the product with η ∂ tφ , we obtain: We have chosen the parameters of the weight η so that δ/µ 2 is small enough. Hence, We also have by using (3.10) and the Young inequality that and by using (3.9) that We thus end up with the elliptic estimate: . Now using the transport equation satisfied byṅ, we can estimate √ η ε∂ tṅ L 2 . This yields: For the last term, we have used that εη µ with the choice of the parameters. We have thus proven that B. The second energy estimate. -The previous computation suggests to look for a control on ε ∇ṅ and ε divu, and this motivates the next energy estimate. We shall prove that . Let us stress that this estimate alone does not allow to conclude: the right hand side involves the full ε∇u, while the left hand side controls only ε divu. Later, we shall establish a third estimate, in order to handle ε∇u and close the argument.
To prove (3.41), we shall again rely on Lemma 1, but I will be handled in a different way. By integration by parts: with The bilinear singular term I 2 will be "compensated" in the end (thanks to the Bohm condition). Therefore we first focus on I 1 .
(a) Study of J 1 . -To evaluate J 1 , the idea is to take the divergence in the equation satisfied byu, in order to express ∆φ in terms ofu andṅ. This reads: ε∂ i (u a + u) · ∇u i + εT i div ∇ṅ n a + n = ε∆φ + ε div r u .
e φa ε div ∇ṅ n a + n (ε divu), The terms K 2 and K 4 can be bounded through standard arguments. With (3.7) and (3.9) in mind, we obtain Note the presence of ε∇u at the r.h.s., due to K 2 . We now evaluate the contribution of the convection term K 1 . Standard manipulations yield The second term can be bounded by a crude L 2 estimate, the third one can be bounded thanks to (3.3)-(3.7): we get The second term is non-positive, and will be one of the crucial terms to control all singular terms (thanks to the Bohm condition). Note finally that the boundary contribution is non-positive: this will also help later. Now we turn to the integral that involves the pressure term, for which we perform an integration by parts e φa (n a + n)(ε divu) · ε∇ṅ, We use the equation satisfied byṅ in order to rewrite the ∇[...] part in L 1 . To this end, we take the ε∇ operator on the transport equation for n to get: where e 1 = (1, 0, 0) , e 2 = (0, 1, 0) , e 3 = (0, 0, 1) . As before the convection term on (ε∇ṅ) is very useful. With the usual manipulations, we obtain: By now standard manipulations, we also have From the estimates on the L i 's, one can deduce an inequality on K 3 . Combining this inequality with (3.44) and (3.45), we end up with and (b) Study of J 2 . -Let us now work on J 2 , which has to be treated very carefully. Indeed, a rough L 2 estimate only shows that this is a singular term in 1/ε, which does not seem small. Instead, we use the transport equation satisfied byṅ to express (n a + n) divu in terms of other quantities. With similar manipulations as before, we obtain: e φa (u a + u) 3 |ṅ| 2 + 1 2 x3=0 Observe here that the last two terms are non-positive. Together with (3.46) and (3.47), the bound (3.43) leads to and (2) Treatment of I 2 . -By using (3.23), we have the straightforward bound: Note that the last term is actually not singular, thanks to (3.10). We now claim thatφ satisfies the following bound: The term 1 2 R 3 + η |(ua+u)3| na+n |ṅ| 2 has the bad sign but will be compensated by other terms later.
To prove estimate (3.49), we write the Poisson equation satisfied byφ: We multiply by mφ, with the weight m := η |(Γu 0 ) 3 |. By a standard energy estimate: We use the Young inequality for the first term, resulting in Next, we observe that thanks to (3.9), we have By the choice of the weight function, we have Finally, since n 0 = e −φ0 , relations (3.22), (3.23) give We thus get that Applying (3.24) to the first term at the r.h.s., we obtain the desired inequality (3.49). It follows that (3) Conclusion of estimate B. -We can then inject (3.48) and (3.51) into (3.42). To ease reading, we set We get e φa (ε divu) ε∂ 3ṅ Here B.T. denotes the boundary contributions, G stands for a good term, that is, and R 3 is for . Let us stress that some singular terms vanish at leading order: indeed, the relation e φ0 = 1/n 0 implies By combining the previous estimate and Lemma 1, we conclude that we have e φa |ṅ| 2 2 with Q B 1 , Q B 2 are quadratic forms determined by the two following symmetric matrices: We recall that µ > 0 is a parameter that can be taken arbitrarily small. Observe that M B 1 is positive for µ sufficiently small, as soon as the Bohm condition (3.11) is satisfied. We also get that M B 2 is positive thanks to the Bohm condition. To obtain (3.41), we still have to get rid of the term µ √ η φ 2 L 2 at the r.h.s. of (3.53). But as the quadratic form Q B 1 controls thanks to (3.50), up to add some good term and some µ √ η ṅ 2 L 2 term at the r.h.s. By taking µ small enough, this singular term is controlled by the l.h.s., which gives (3.41).
This estimate by itself allows us almost to conclude; the only problematic term comes from the singular contribution R 3 + η |ε ∇u| 2 in the right hand side which involves a full derivative ofu, and that will be treated below.
C. The third energy estimate. -The previous computations suggest that we also need an estimate of a full ε-derivative ofu. Note that we can not work with the curl ofu (together with the previous estimate on the divergence), because we would not be able to control the trace ofu on the boundary {x 3 = 0}.
We shall prove that (3.54) (ṅ,u, ε∇ṅ, ε∇u) 2 . We shall use again the estimate of Lemma 1, but through a different approach of the integral I . The idea will be to combine the estimate of Lemma 1 with an estimate on the ε-derivatives of the system to cancel the singular term through the Poisson equation.
Let us first first apply ε-derivatives on the equation on the velocity fieldu: for i = 1, 2, 3, this yields n a + n We then take the scalar product with η ε ∂ iu and make the sum in i. We find The term S above will refer from now and on to any admissible singular term, that is, a quantity which is singular but bounded as The first term in S will be absorbed by positive terms in the left hand side for µ small enough, while the second term will be absorbed by the left hand sides of Estimates A and B. The term G above will refer from now and on to a good term, that is, satisfying (3.52).
(1) Treatment of I 1 . -We use the identity: We shall now use the transport equation satisfied by ε ∂ iṅ , which reads: Therefore, we have, relying on the standard manipulations: We end up with (2) Treatment of I 2 .
-After some integration by parts, we can rewrite I 2 as Then we have where we set Similarly to B stands for an admissible boundary term which is bounded by boundary terms in the left hand sides of Estimates A and B, that is to say Putting together (3.60), (3.62) and (3.56), we deduce (3) Treatment of I . -We recall that the I term defined in Lemma 1 can be written as with J := − R 3 + ηφ (n a + n) divu. We shall study I (and J) in a different way than that followed for the first two energy estimates. The idea is to combine J with I 2 through the term L (see (3.61)-(3.21)). Therefore, we decompose J as follows: Since by construction of the approximation solution, we have e −φ 0 = n 0 , J 2 can be bounded as follows: Considering J 1 , by the Poisson equation, we have: One can rely on the transport equation satisfied byṅ to replace divu, which yields: which we shall recast, after the usual manipulations, as We obtain Consequently, by combining (3.63) and (3.64), we obtain Here B.T. denotes the important boundary contributions: (3) Conclusion of estimate C. -By combining, (3.65) and Lemma 1, we obtain where Q C is the quadratic form of the symmetric matrix: which is positive thanks to the Bohm condition (3.11).
To conclude, we can handle the two first terms in the right hand side of the above estimates by using the Young inequality. Indeed, we write and for µ sufficiently small, the term √ η ε∇u 3 2 L 2 can be absorbed in the left hand side since Q C is positive while the other term can be incorporated in S . We proceed in the same way for the boundary term by writing and we absorb the first term in the left hand side by the positivity of Q C while the other term can be incorporated in B.
To obtain (3.54), it suffices to observe that the term µ √ η ε∇u 2 L 2 can be absorbed in the left hand side for µ sufficiently small and to integrate in time. left hand side of (3.41). This yields (3.66) (ṅ,u,φ, ε∇ṅ, ε divu, ε∇φ) 2 To handle the singular term µ √ η ε∇u 2 L 2 T L 2 (R 3 + ) in the right hand side of (3.66), we consider (3.66) + √ µ(3.54). To ease reading, we set , and To conclude, we observe that by taking µ sufficiently small, the singular terms can be absorbed in the right hand side. We end the proof by integrating in time, by applying a Gronwall estimate and by using (3.10).

Nonlinear stability
In this section, we shall prove our main Theorem 1. As already mentioned, we shall actually prove a more precise version, see Theorem 3.
Let us write the solution (n ε , u ε , φ ε ) of (1.1) under the form where n a , u a , φ a are shorthands for n ε app , u ε app , φ ε app (defined in (2.1)). Then, we get for (n, u, φ) the system (4.1) ∂ t n + (u a + u) · ∇n + n div(u + u a ) + div(n a u) = ε K R n , ∂ t u + (u a + u) · ∇u + u · ∇u a + T i ∇n n a + n − ∇n a n a n n a + n = ∇φ + ε K R u , together with the boundary condition (4.2) φ| x3=0 = 0 and the initial condition Observe here that in (4.1), ε K R n , ε K R u and ε K+1 R φ are remainders that appear because (n a , u a , φ a ) is not an exact solution of (1.1).
The main result of this section is Let K ∈ N * , K m and (n a , u a , φ a ) an approximate solution at order K, given by Theorem 2, which is defined on [0, T 0 ]. There exists ε 0 > 0 and δ 0 > 0 such that for every δ ∈ (0, δ 0 ] and for every ε ∈ (0, ε 0 ], there is C > 0 independent of ε such that the solution of (4.1)-(4.2)-(4.3) is defined on [0, T 0 ] and satisfies the estimate We thus obtain: -Under the assumptions of Theorem 3, we have In particular, we get the L 2 and L ∞ convergences as ε → 0: Then Theorem 1 is a straightforward consequence of this corollary. From now on, our goal is to prove Theorem 3. First, for each ε > 0, one can solve the Poisson equation in (4.1), and express φ in terms of n; more precisely, by elliptic regularity, ∇φ can be seen as a semi-linear term in n and the known local existence results for the compressible Euler equation with strictly dissipative boundary conditions ( [3]) can be applied to the system (4.1). Let us assume that (n 0 , u 0 ) ∈ H m (R 3 + ) for m 3, then there exists T ε > 0 and a unique solution of (4.1) defined on [0, T ε ] such that u ∈ C ([0, T ε ), H m ) and that there exists M > 0 independent of ε such that the assumptions (3.10) and (3.11) of Proposition 2 are satisfied with Thanks to the well-posedness in H m for m 3 of the system (4.1), we can define where r is chosen such that (4.7) 5/2 < r < K and the H m ε norm is defined by The difficulty is thus to prove that the solution actually exists on an interval of time independent of ε. We shall get this result by proving uniform energy estimates combined with the previous local existence result when the initial data and the source term are sufficiently small (i.e., when the approximate solution (n a , u a ) is sufficiently accurate).
For m 3, we shall prove by using a priori estimates that under the assumptions of Theorem 3, we have for every T ∈ [0, T ε ] the control To prove this estimate, we can apply the operator Z α := (ε∂) α to (4.1) for |α| m − 1. We obtain for Z α (n, u, φ) the system (4.9)          ∂ t Z α n + (u a + u) · ∇Z α n + (n a + n) div Z α u = C n + ε K Z α R n , ∂ t Z α u + (u a + u) · ∇Z α u + T i ∇Z α n n a + n = ∇Z α φ + C u + ε K Z α R u , One has h = h 0 for |α| = 0 and h = h 1 for |α| 1. The functions C n , C u and C φ are remainders mostly due to commutators; in C n , we include the term −Z α [u · ∇n a ] − Z α [n div u a ] as well, while in C u we also include −Z α [u · ∇u a ] + T i Z α ∇na na n na+n . One can observe that this corresponds to the "abstract" system that we have studied in Proposition 2.
The assumptions (3.10), (3.11) of Proposition 2 are matched on [0, T ε ] for ε sufficiently small since they are verified by the approximate solution. The only estimate that does not follow directly from the definition of T ε is the estimate of ∂ t (n, φ) L ∞ . For ∂ t n, we immediately get by using the first line of (4.1) that for M sufficiently small. In a similar way, by taking the time derivative in the elliptic equation for φ, we have (4.10) − ε 2 ∆∂ t φ + e −(φa+φ) ∂ t φ = −∂ t n + ∂ t (e −φa )(e φ − 1) + ε K+1 ∂ t R φ .
By using the maximum principle for this elliptic equation, we thus get that To use the result of Proposition 2, we also need to check that the remainders C n , C u can be split as in the assumption (3.9).
To describe the structure of these remainders, we can start with the study of C n . We can distinguish between a linear part which corresponds to C l n = −[Z α , u a ] · ∇n − Z α (n div u a ) − [Z α , n a ] div u − Z α (u · ∇n a ) and a nonlinear part which is C nl n = −[Z α , u] · ∇n − [Z α , n] div u. By using standard tame estimates in Sobolev spaces, the nonlinear part can be seen as a regular part in the decomposition (3.9). Indeed, we easily get that for some C > 0 independent of ε, we have C nl n H 1 ε C ∇u L ∞ + ∇n L ∞ (n, u) H |α|+1 ε and hence by using the Sobolev embedding, we obtain that on [0, T ε ], C nl n H 1 ε Cε r−5/2 (n, u) H m ε . To estimate the linear part C l n , we can use the Leibniz formula and the estimates (3.12), (3.13) on the approximate solution. By using the ε weighted derivatives, the terms that lead to singular terms are the ones of the first type ε k ∂ β n a · ∇∂ γ u or the ones of the second type ε k ∂ x3 ∂ β n a · ∂ γ u with |β| + |γ| = k. Indeed, for the terms of the second type, we obtain a singular term as soon as β involves only x 3 derivatives and when they all hit the boundary layer term in n a , in this case, we can use (3.7) to estimate them by (4.11) C a µ η | ε∂ m (n, u)| + | ε∂ m (n, u)| with the notation | ε∂ m (n, u)| = |β| m |(ε∂) β (n, u)|.
We thus see the first term in (4.11) as a singular part of the source term and the second one as a regular part in the decomposition (3.9). For the terms of the first type, in order to estimate them with u in a H m ε space, we need to write them under the form ε |β|−1 ∂ β n a · (ε∇)(ε∂) γ u and hence we obtain again a singular term as soon as ∂ β is made only of x 3 derivatives and that they all hit the boundary layer terms in n a . By using (3.7), we can estimate this contribution again by C a µ η | ε∂ m (n, u)| + | ε∂ m (n, u)| We thus obtain an estimate |C l | + |ε∇C l | C a µ η | ε∂ m (n, u)| + | ε∂ m (n, u)| .
We can proceed in the same way to handle C u . For example, for the commutator C p := Z α , 1 n a + n ∇n, we can expand β,γ Z β n a ∇Z γ n + Z β n∇Z γ n , where β,γ stand for terms which are uniformly bounded in L ∞ on [0, T ε ]. By using the same arguments as before, we thus get that |C p | + |ε∇C p | C(C a , M ) µ η | ε∂ m (n, u)| + | ε∂ m (n, u)| + C(C a , M )| C p | where | C p | is bounded in L 2 by C(C a , M )ε r−5/2 n H m ε by using the usual product estimates in Sobolev spaces. For ε −1 C φ , we can use again the bounds on the approximate solution (3.12), (3.13) and standard tame estimates in Sobolev spaces. Indeed, we have that ε −1 C φ can be expanded as a sum of terms either of the form ε −1 (Z α e −φa )(e −φ − 1) which is bounded by C(C a , M ) µ η |φ| + |φ|) or of the form ε −1 Z β (e −φa ) Z γ φ e −φ with |β| + |γ| m which can also be bounded by C(C a , M ) µ η |φ| + |φ|) or ε −1 Z β (e −φa ) Z γ1 φ · · · Z γr φ e −φ with r 2, |β| + |γ 1 | + · · · + |γ r | m and β = 0. Since this term is at least quadratic in φ, we can use tame estimates to bound it in L 2 by We thus get that we can split ε −1 C φ according to (3.9): By using the above expansions of C φ , we also easily obtain that In summary, we have proven that the commutator terms can be split under the form C n , C u , ε∇C n , ε∇C u , C φ /ε, ∂ t C φ C(C a , M ) µη C s + C r with η C s Consequently, by using Proposition 2, we obtain for every T ∈ [0, T ε ] (4.12) √ µ (n, u, φ) 2 To conclude, it remains to estimate the terms involving ∂ t φ in the right hand side of the above estimate. By using the elliptic equation (4.10) and standard energy estimates, we get that For µ sufficiently small, the singular term µ √ η (n, u, φ) 2 We have thus proven (4.8).
The parameter µ can now be considered as fixed. By standard continuation arguments, we next obtain that for ε sufficiently small T ε T 0 and that on [0, T 0 ], we have This ends the proof of Theorem 3.