Francesco Fanelli
email: fanelli@math.univ-lyon1.fr

Ewelina Zatorska
email: e.zatorska@imperial.ac.uk

Low Mach number limit for degenerate Navier-Stokes equations in presence of strong stratification

Keywords: compressible Navier-Stokes equations, density-dependent viscosity, low Mach and Froude numbers, strong stratification, cold pressure. 2020 Mathematics Subject Classification: 35Q35 (primary), 35B40, 76M45, 35B25 (secondary)

In this paper, we investigate the low Mach and low Froude numbers limit for the compressible Navier-Stokes equations with degenerate, density-dependent, viscosity coefficient, in the strong stratification regime. We consider the case of a general pressure law with singular component close to vacuum, and general ill-prepared initial data. We perform our study in the three-dimensional periodic domain. We rigorously justify the convergence to the generalised anelastic approximation, which is used extensively to model atmospheric flows.

Introduction

Flows in the atmosphere are typically characterised by two main features (see [START_REF] Pedlosky | Geophysical fluid dynamics[END_REF]): first of all, they are weakly compressible, moreover they undergo the combined effect of a strong stratification (due to the action of gravity) and of a strong Coriolis force (due to the rotation of the Earth, which is very fast if compared to the space-time scales of the flows).

Neglecting the effects of the Earth rotation, the importance of the other two factors, i.e. weak compressibility and strong stratification, may be assessed by introducing two physical a-dimensional parameters, the Mach number and the Froude number, respectively. The smaller these parameters are, the more predominant weak compressibility and strong stratification become. Thus, as usual in Physics, it is natural to look at the regime where both parameters vanish, to find reduced models for atmospheric flows. They are simpler to deal with than the corresponding primitive system, both from the analytical and numerical point of view.

When the Mach number and the Froude number go to zero with the same speed, the flow becomes incompressible and stratified at the same rate. Formally, this asymptotic regime was considered already by Ogura and Phillips in [START_REF] Ogura | Scale analysis for deep and shallow convection in the atmosphere[END_REF]. The limiting system takes the name of anelastic approximation. The physical importance of the anelastic approximation is discussed, for example in [START_REF] Klein | Scale-dependent asymptotic models for atmospheric flows[END_REF] in the context of various atmospheric flows, and in [START_REF] Almgren | Low Mach number modeling of type Ia supernovae[END_REF] in the context of astrophysics models.

The primitive system and the limit dynamics

In this paper we will give a rigorous derivation of what we call the generalised anelastic approximation, namely an anelastic approximation with variable viscosity. The starting system (referred to as the primitive system) is the barotropic Navier-Stokes equations, with bulk viscosity coefficient equal to 0 and the shear viscosity coefficient proportional to the density of the fluid. In particular, the system strongly degenerates close to vacuum. This choice of the viscosity coefficients is physically relevant, as viscosity is, in general, hardly expected to be uniform for flows on large scales. Their specific form allows one to exploit a certain mathematical structure of the system, called the BD-entropy (see more details in the discussion below).

Assuming that both the Mach and Froude numbers are equal to a small parameter ε > 0, the system of equations reads as follows:

∂ t + div( u) = 0, ∂ t ( u) + div( u ⊗ u) + 1 ε 2 ∇p( ) -ν div( Du) = 1 ε 2 ∇G. (1)
The unknowns are the mass density = (t, x) ≥ 0, and the velocity vector field u = u(t, x) ∈ R 3 .

The function p = p( ) denotes the internal pressure, the constant ν > 0 is the viscosity coefficient, and D = 1 2 ∇ + ∇ t is the symmetric part of the gradient. Finally, G = G(x) is a smooth function (say G ∈ C 3 (Ω)) describing a scalar external force acting on the flow. G typically encodes the action of gravity, in which case G = -gx 3 , where g is the gravitational acceleration constant.

Due to the present state-of-the-art of the mathematical theory for system [START_REF] Almgren | Low Mach number modeling of type Ia supernovae[END_REF], we assume that the fluid occupies the periodic box in R 3 , i.e. we consider the equations (1) on the space domain

Ω := T 3 . (2)
The pressure function p is assumed of the following form:

p( ) = p e ( ) + p c ( ) = 1 γ γ - 1 κ -κ , γ > 1, κ ≥ γ -2, κ > 3. (3)
The first part is the standard barotropic pressure, while the second part is the so-called "cold pressure", because it is most significant in the region of temperatures close to zero. The constants 1/γ and 1/κ are just normalisation factors; their presence is useful in some computations. The restriction on the adiabatic exponent γ > 1 is somehow classical in the theory of compressible Navier-Stokes equations. The conditions on the exponent κ, instead, are of technical nature.

When ε → 0 + in equations [START_REF] Almgren | Low Mach number modeling of type Ia supernovae[END_REF], we observe a competition between the large size of the pressure term (low Mach number effect), which tends to drive the flow to incompressibility, and the large size of the forcing term (low Froude number effect), which increases the stratification of the flow. Due to the choice of scaling, those two terms are in balance in the limit process. Therefore, it is easy to see that, when ε → 0 + , will tend, at least formally, to the profile b = b(x) satisfying

∇p(b) = b ∇G . (4)
Smoothness of G(x) and strict monotonicity of p imply that there exists a smooth function b ∈ C 3 (Ω) satisfying [START_REF] Bresch | On some compressible fluid models: Korteweg, lubrication, and shallow water systems[END_REF]. Monotonicity of p implies convexity of the pressure potential H (defined in [START_REF] Fanelli | Highly rotating viscous compressible fluids in presence of capillarity effects[END_REF] below) which provides existence of constants b * , b * ∈ R such that

∀ x ∈ Ω , 0 < b * ≤ b(x) ≤ b * . (5)
Note that, for p( ) = 2 2 , one gets G = b up to the choice of an irrelevant additive constant, which is the case considered in [START_REF] Bresch | An example of low Mach (Froude) number effects for compressible flows with nonconstant density (height) limit[END_REF]. On the other hand, if G(x) = -gx 3 is the gravitational potential, it is easy to see that b = b(x 3 ) verifies [START_REF] Bresch | An example of low Mach (Froude) number effects for compressible flows with nonconstant density (height) limit[END_REF].

Since ≈ b for ε → 0 + , assuming that we can identify the limits of the non-linear terms appearing in [START_REF] Almgren | Low Mach number modeling of type Ia supernovae[END_REF], it is not difficult to check that (formally) the limiting system is an anelastic approximation with variable viscosity coefficient, namely

div(b U) = 0 ∂ t U + (U • ∇)U + ∇π -b -1 ν div b DU = 0 . (6)
We refer to system [START_REF] Danchin | Zero Mach number limit for compressible flows with periodic boundary conditions[END_REF] as the generalised anelastic approximation. In the above system, π = π(t, x) denotes an unknown scalar function, and the term ∇π plays the role of a Lagrangian multiplier associated with the anelastic constraint div (b U) = 0. The limiting system can be also regarded as the viscous counterpart of the so-called lake equation, whose study was initiated in [START_REF] Levermore | Global well-posedness for models of shallow water in a basin with a varying bottom[END_REF].

The goal of this paper is to rigorously justify the above formal derivation in the framework of global in time finite energy weak solutions to the primitive system (1)-(3) for general ill-prepared initial data.

An overview of related results

Due to the physical importance of the anelastic approximation, its rigorous derivation has been the object of intensive studies in the past years.

In [START_REF] Masmoudi | Rigorous derivation of the anelastic approximation[END_REF], Masmoudi proved the rigorous derivation of the anelastic approximation, starting from the classical barotropic Navier-Stokes system. He considered a bounded domain in R 3 , supplemented with Dirichlet boundary conditions, and the limit was performed for ill-prepared initial data via a compensated compactness argument. Soon after that, Feireisl, Málek, Novotný and Straškraba proved an analogous result on a periodic box, and for pressure laws which are small variations of the ideal gas law, see [START_REF] Feireisl | Anelastic approximation as a singular limit of the compressible Navier-Stokes system[END_REF]. Finally, we refer to the book by Feireisl and Novotný [START_REF] Feireisl | Singular limits in thermodynamics of viscous fluids[END_REF] for a complete account of the mathematical literature on the low Mach number limit, in the presence of both low and high stratification effects. They presented the theory for the full Navier-Stokes-Fourier system and in the framework of global in time finite energy weak solutions. The literature related to the incompressible limit of compressible fluid equations is of course much more extensive (see e.g. the pioneering works [START_REF] Ebin | The motion of slightly compressible fluids viewed as a motion with strong constraining force[END_REF] by Ebin and [18]- [START_REF] Klainerman | Compressible and incompressible fluids[END_REF] by Klainerman and Majda) and recalling all relevant results goes far beyond the scope of this introduction. We thus limit ourselves to quote a couple of recent works.

In [START_REF] Feireisl | On singular limits arising in the scale analysis of stratified fluid flows[END_REF], a variant of the anelastic approximation was derived, starting from a version of the Navier-Stokes-Fourier system with neglected thermal diffusion: the potential temperature is assumed to be just transported by the velocity field. The limit system that is identified in [START_REF] Feireisl | On singular limits arising in the scale analysis of stratified fluid flows[END_REF] reads as a coupling of the anelastic approximation system (6) with a transport equation for the limiting temperature. The convergence is proven in the infinite slab R 2 × ]0, 1[ through a spectral analysis of the singular perturbation operator and an application of the celebrated RAGE theorem from scattering theory. The advantage of that technique, in comparison to the one used in [START_REF] Masmoudi | Rigorous derivation of the anelastic approximation[END_REF], is that it allows to get the compactness of the sequence of velocity fields.

On the side of the incompressible limit (with no stratification effects, though), another interesting result is [START_REF] Danchin | Compressible Navier-Stokes system: large solutions and incompressible limit[END_REF], where the the authors deal with weakly compressible viscous fluids in a critical regularity functional framework. In that paper, weak compressibility is obtained by taking a large bulk viscosity coefficient limit, instead of the classical low Mach number limit. More recently, in [START_REF] Fanelli | Incompressible and fast rotation limit for barotropic Navier-Stokes equations at large Mach numbers[END_REF], a similar idea was implemented for fast rotating fluids.

For the degenerate Navier-Stokes system (1), as considered in our paper, the relevant results are much more sparse. The first one to mention is the incompressible limit in a strong stratification regime considered in [START_REF] Bresch | An example of low Mach (Froude) number effects for compressible flows with nonconstant density (height) limit[END_REF] by Bresch, Gisclon and Lin. Their system includes two artificial drag terms in the momentum equation, in order to improve the available information for the velocity field close to vacuum. The convergence to the anelastic approximation is proven using the relative energy method, for a special choice of pressure law p( ) = 2 /2 and for well-prepared initial data.

A similar method was used in [START_REF] Bresch | Existence of global weak solution for a 2D viscous shallow water equation and convergence to the quasi-geostrophic model[END_REF] by Bresch and Desjardins and in [START_REF] Jüngel | An asymptotic limit of a Navier-Stokes system with capillary effects[END_REF] by Jüngel, Lin and Wu, for the low Mach number and low Rossby number limit in a two-dimensional geometry. The external force in these papers is replaced by the Coriolis force (whence the low Rossby number regime) and a capillarity term. The resulting limiting system is a quasi-geostrophic type equation for the stream function of the limit velocity field. We also refer to [START_REF] Fanelli | Highly rotating viscous compressible fluids in presence of capillarity effects[END_REF]- [START_REF] Fanelli | A singular limit problem for rotating capillary fluids with variable rotation axis[END_REF] for a generalisation of these results to the 3-D setting and to the case of general ill-prepared initial data.

The content of the paper

In the context depicted above, our work can be seen as a generalisation of the result from [START_REF] Bresch | An example of low Mach (Froude) number effects for compressible flows with nonconstant density (height) limit[END_REF], to the case of ill-prepared initial data and of more general pressure laws (and hence, more general external forces G). We work in the framework of global in time finite energy weak solutions to the primitive system. Their existence, in presence of a cold part of the pressure (3), has been established in [START_REF] Zatorska | On the flow of chemically reacting gaseous mixture[END_REF], [START_REF] Mucha | Approximate solutions to a model of two-component reactive flow Discrete Contin[END_REF]. The case without this assumption has been completed much more recently in [START_REF] Vasseur | Existence of global weak solutions for 3D degenerate compressible Navier-Stokes equations[END_REF] by Vasseur and Yu. In all these results, the finite energy condition plays, of course, a major role. However, the degenerate Navier-Stokes system (1) possesses also a second energy inequality, usually named BD entropy inequality after Bresch and Desjardins, who investigated this second energy conservation law in [START_REF] Bresch | Existence of global weak solution for a 2D viscous shallow water equation and convergence to the quasi-geostrophic model[END_REF].

The BD entropy estimate provides a control on the gradient of a certain function of the density, whose exact form depends on the form of the viscosity coefficient. For system (1), this function is ∇ √ . The BD entropy also allows to control the skew-symmetric part A = 1 2 ∇ -∇ t of the gradient of the velocity. This, when combined with the classical energy, provides the corresponding bound for the full gradient of the velocity.

The classical energy inequality, the BD entropy inequality, and all the bounds that follow, are essential also in the present paper. As a matter of fact, for any value of parameter ε ∈ ]0, 1], we consider a finite energy weak solution ε , u ε to system (1), which satisfy both those energy inequalities. However, proving that the BD entropy estimate is satisfied uniformly with respect to ε requires some effort, especially when a general pressure law is considered: this is one of the first problems solved in our paper.

Having all these estimates satisfied uniformly for a sequence of finite energy weak solutions ε , u ε ε , the rest of the proof of the derivation of the generalised anelastic approximation (6) boils down to showing that the weak limit (b, U) is indeed a solution to [START_REF] Danchin | Zero Mach number limit for compressible flows with periodic boundary conditions[END_REF]. It is well known that passing to the limit in the weak formulation of equations [START_REF] Almgren | Low Mach number modeling of type Ia supernovae[END_REF], especially in its nonlinear parts, is problematic. This is because the singular terms, along with the ill-prepared initial data, are responsible for fast time oscillations of the solutions (the so-called acoustic waves), which may prevent, in the end, the convergence of the nonlinear terms to the expected limit. Showing that this does not happen is the core of the whole proof.

The main concern is the convergence of the convective term in the momentum equation. To that purpose, we use a different technique than the one from [START_REF] Bresch | An example of low Mach (Froude) number effects for compressible flows with nonconstant density (height) limit[END_REF]. Our approach is inspired by the previous works [START_REF] Lions | Incompressible limit for a viscous compressible fluid[END_REF] and [START_REF] Masmoudi | Rigorous derivation of the anelastic approximation[END_REF] on the incompressible limit for the classical barotropic Navier-Stokes system, and is based on a compensated compactness argument. More precisely, we first regularise the primitive equations, which we recast in the form of a wave system. After that, we exploit two pieces of information coming from the wave system. First of all, we may deduce the compactness of the rotational part of the velocity fields. On the other hand, by direct but elaborated algebraic manipulations, we may infer that the interaction of the potential part of the velocity fields in the convective term gives rise to small quantities, which tend to vanish when ε → 0 + . It is worth to point that this argument is robust enough to deal with other variants of the system [START_REF] Almgren | Low Mach number modeling of type Ia supernovae[END_REF]. For instance, we could trade the cold component of the pressure function, which basically provides us with some integrability properties for ∇u, for a turbulent drag term |u|u, which would give a better integrability of the momentum V := u. In that case, most of the steps are the same, although the derivation of essential estimates becomes significantly more laborious. The only problem, and the breaking point, arises when one wants to pass to the limit in this artificial drag term. This term turns out to be even more non-linear than the convective term, because of the presence of the norm |u| of the vector u. It is not clear how to bypass this difficulty in our framework, and so, the problem remains open.

We conclude with a short outline of the paper. In the next section, we collect our main hypotheses on the initial data, we give the definition of finite energy weak solutions, and we state our main result. In Section 3, we deduce, from the energy inequality and the BD entropy inequality, a long list of uniform bounds for the family of solutions ε , u ε ε we consider. That part of the study is rather delicate, due to the degeneracy of the system close to vacuum. In Section 4, we use the previous uniform bounds to extract a weakly convergent subsequence, and to derive first basic properties on its weak limit point. At this stage we reformulate the primitive equations into the wave system, which describes the propagation of the acoustic waves. In Section 5 we rigorously perform the convergence in the weak formulation of equations [START_REF] Almgren | Low Mach number modeling of type Ia supernovae[END_REF], and conclude the derivation of the anelastic approximation [START_REF] Danchin | Zero Mach number limit for compressible flows with periodic boundary conditions[END_REF]. For the convenience of the reader, we collect some tools from Fourier analysis which we need in our study in the Appendix at the end of the paper.

Statement of the main result

In this section, we first introduce our assumptions on the initial data, then we define the notion of finite energy weak solutions to system (1)-( 3), and finally we formulate our main theorem. 3) is supplemented by general ill-prepared initial data. Namely, for any small parameter ε ∈ ]0, 1] fixed, we pick initial data

Initial data. Problem (1)-(

( , u)| t=0 = 0,ε , u 0,ε (7)
satisfying the following conditions:

(i) the initial densities 0,ε ≥ 0 are assumed to be small perturbations of the static state b, defined by (4): more precisely, we assume1 that

0,ε = b + ε φ 0,ε , with φ 0,ε ε ⊂ L ∞ (Ω) and ∇ ln 0,ε b ε ⊂ L ∞ (Ω) ;
(ii) the initial velocity fields u 0,ε are such that u 0,ε ε ⊂ L ∞ (Ω).

Thus, up to extraction of a subsequence, not relabeled here, we may suppose that

φ 0,ε → φ 0 and u 0,ε → u 0 weakly- * in L ∞ (Ω) . (8)
Energy functionals. Next, we need to introduce various energy functionals. The internal energy function (sometimes called pressure potential ) is defined by the ODE

H ( ) -H( ) = p( ) ,
which implies in particular that

H ( ) = p ( ) .
Notice that H is defined up to the sum of an affine function. Here, we fix the classical choice

H( ) = 1 p(z) z 2 dz = 1 γ (γ -1) γ-1 -1 + 1 κ (κ + 1) -κ-1 -1 . (9)
We now denote

E , u b := 1 2 Ω |u| 2 dx + 1 ε 2 Ω H( ) -H(b) -H (b) ( -b) dx (10)
F , u b := Ω u + ν∇ ln b 2 dx (11)
to be the classical energy and the BD entropy functions. We also set E , u b (T ) := E (T ), u(T ) b and E , u b (0) := E 0 , u 0 b , and similarly for the function F.

Weak solutions to the primitive system. After this preparation, we are ready to give the definition of weak solutions to system (1)-(3) which are relevant for us.

Definition 2.1 Let 0 , u 0 be such that E 0 , u 0 b + F 0 , u 0 b < +∞.

We say that the couple ( , u) is a finite energy weak solution of (1)-( 3) in [0, T [ ×Ω, with the initial datum 0 , u 0 , provided the following conditions are satisfied: (1) ≥ 0 almost everywhere, with

∈ L ∞ [0, T [ ; L γ (Ω) and -1 ∈ L ∞ [0, T [ ; L κ (Ω) , ∇ √ ∈ L ∞ [0, T [ ; L 2 (Ω) and p ( ) ∇ ∈ L 2 [0, T [ ; L 2 (Ω) ; (2) √ u ∈ L ∞ [0, T [ ; L 2 (Ω) and √ ∇u ∈ L 2 [0, T [ ; L 2 (Ω) ;
(3) the equations of system (1) are satisfied in the sense of distributions: more precisely, we have

Ω 0 ξ(0) dx + T 0 Ω ∂ t ξ + u • ∇ξ dx dt = 0 (12)
for any test function ξ ∈ D [0, T [ ×Ω , and

Ω 0 u 0 • ψ(0) dx + T 0 Ω u • ∂ t ψ + u ⊗ u : ∇ψ dx dt + 1 ε 2 T 0 Ω p( ) div ψ dx dt + 1 ε 2 T 0 Ω ∇G • ψ dx dt -ν T 0 Ω Du : ∇ψ dx dt = 0 (13)
for any test function ψ ∈ D [0, T [ ×Ω; R 3 ;

(4) for almost every t ∈ [0, T [ , the following energy inequalities hold true:

E , u b (t) + ν t 0 Ω |Du| 2 dx ds ≤ E 0 , u 0 b , F , u b (t) + ν ε 2 t 0 Ω b 2 p ( ) ∇ b 2 dx ds + ν t 0 Ω |Au| 2 dx ds ≤ C 0 e C 0 (1+T ) , (14)
where the constant C 0 > 0 may depend on E 0 , u 0 b and F 0 , u 0 b but is independent of ε.

The solution ( , u) is said global in time if the previous properties hold true for any T > 0.

For any ε ∈ ]0, 1] fixed, the existence of global in time finite energy weak solutions to system (1) in the sense of previous definition was proven in [START_REF] Zatorska | On the flow of chemically reacting gaseous mixture[END_REF] and [START_REF] Mucha | Approximate solutions to a model of two-component reactive flow Discrete Contin[END_REF], in the case G = 0 (corresponding to b = const.). The argument of those papers apply in a fairly direct way also to the case considered in this paper, where G = 0 and b is non-constant: we explain in the next section how to modify the estimates of [START_REF] Zatorska | On the flow of chemically reacting gaseous mixture[END_REF]- [START_REF] Mucha | Approximate solutions to a model of two-component reactive flow Discrete Contin[END_REF] in order to include the force.

Main result. Before stating the main result of this paper, we need some additional tools and notation. Following [START_REF] Masmoudi | Rigorous derivation of the anelastic approximation[END_REF]- [START_REF] Feireisl | Anelastic approximation as a singular limit of the compressible Navier-Stokes system[END_REF] (see also [START_REF] Lions | Mathematical topics in Fluid Mechanics[END_REF]), we introduce the twisted Leray-Helmholtz projector P b , related to the smooth function b satisfying (5), as follows: for any smooth vector field v on Ω, we write b (Ω; R 3 ), which is defined as the space of functions f : Ω -→ R 3 which are L 2 -summable with respect to the measure 1 b dx. Similarly to the case of the classical Leray-Helmholtz projector P = P 1 and its L 2 -orthognal projector Q = Q 1 , it is possible to prove that both P b and Q b are bounded continuous functionals on L p (Ω; R 3 ), for any 1 < p < +∞.

v = P b [v] + b ∇Ψ ,
We can now state the main result of this paper, which is contained in the following theorem.

Theorem 2.2 Let γ > 1 and κ ≥ γ -2, κ > 3 in (3). Let 0,ε , u 0,ε ε be a family of initial data satisfying hypotheses (i)-(ii) fixed above, so in particular the condition

sup ε∈ ]0,1] E 0,ε , u 0,ε b + sup ε∈ ]0,1] F 0,ε , u 0,ε b < +∞ ( 15
)
holds. Define the couple (φ 0 , u 0 ) as in [START_REF] Ebin | The motion of slightly compressible fluids viewed as a motion with strong constraining force[END_REF]. Let ε , u ε ε be a family of global in time weak solutions to system (1)-(3), in the sense of Definition 2.1, corresponding to the previous initial data. Define the scalar quantity φ ε := ε-b ε . Then, there exists a couple of functions φ, U such that, passing to a suitable subsequence as the case may be, in the limit ε → 0 one has

ε → b strongly in L ∞ loc R + ; L p (Ω) , for any p ∈ [1, 3[ , φ ε → φ weakly in L 2 loc R + ; W 1,min{γ,2} (Ω) , u ε → U weakly- * in L ∞ R + ; L p 1 (Ω) ∩ L 2 loc R + ; W 1,p 1 (Ω) ,
where p 1 := 2κ κ + 1 .

In addition, φ = φ(b) is a function of the static profile b, while U is a solution of the target system (6) in the weak sense, related to the initial datum

U| t=0 = U 0 := 1 b P b [b u 0 ],
i.e. one has div(bU) = 0 almost everywhere in R + × Ω and

Ω b u 0 • ψ(0) dx + T 0 Ω b U • ∂ t ψ + b U ⊗ U : ∇ψ dx dt -ν T 0 Ω b ∇U : ∇ψ dx dt = 0 (16)
for any T > 0 and any test function ψ ∈ D [0, T [ ×Ω; R 3 such that div(bψ) = 0.

Remark 2.3 Note that the initial condition equals

Ω b u 0 • ζ dx = Ω P b [b u 0 ] • ζ dx = Ω b U 0 • ζ dx
for any test function ζ ∈ D Ω; R 3 such that div(bζ) = 0.

A priori estimates

In this section, we derive uniform bounds for the family of weak solutions ε , u ε ε to the original Navier-Stokes system (1). The main tools for this are the classical energy inequality and the so-called BD entropy estimate.

Here and everywhere in the text, we adopt the following notation: given a Banach space X and any p ∈ [1, +∞], we set L p T (X) := L p ([0, T ]; X); in the case T = +∞, instead, we explicitly write L p (R + ; X). When convenient, we will use also the notation L p loc (R + ; X) := T >0 L p T (X).

Bounds coming from the energy inequality

The energy inequality for ε , u ε , which is satisfied by assumption, reads as follows: for almost any time T > 0, we have

E ε , u ε b (T ) + ν T 0 Ω ε |Du ε | 2 dx dt ≤ E 0,ε , m 0,ε b , (17)
where the function E , u b has been defined in [START_REF] Fanelli | A singular limit problem for rotating capillary fluids with variable rotation axis[END_REF] and we recall that we have set

E ε , u ε b (T ) := E ε (T ), u ε (T ) b .
Notice that, due to the cold pressure, at any value of ε ∈ ]0, 1] fixed, the velocity field u ε is well-defined , thus the previous notation makes sense.

From the energy inequality [START_REF] Jüngel | An asymptotic limit of a Navier-Stokes system with capillary effects[END_REF], we now derive first uniform bounds for the family ε , u ε ε . In fact, owing to our assumptions on the initial data, and in particular to [START_REF] Feireisl | Singular limits in thermodynamics of viscous fluids[END_REF], the right-hand side of ( 17) is uniformly bounded : specifically, we have

sup ε∈ ]0,1] E 0,ε , u 0,ε b < +∞ .
Then, it is easy to deduce the following uniform bounds:

√ ε u ε ε ⊂ L ∞ R + ; L 2 (Ω) , (18)
√ ε Du ε ε ⊂ L 2 R + ; L 2 (Ω) . (19)
Let us now focus on the density functions. To begin with, following the approach of [START_REF] Feireisl | Singular limits in thermodynamics of viscous fluids[END_REF], it is convenient to decompose any function h into its essential and residual parts. Thus, for almost every time t > 0 and all ε ∈ ]0, 1], we introduce the sets

Ω ε ess (t) := x ∈ Ω b * 2 ≤ ε (t, x) ≤ 2 b * , Ω ε res (t) := Ω \ Ω ε ess (t) ,
where the constants b * and b * have been defined in (5) Then, given a function h, we can write

h = [h] ess + [h] res , where [h] ess := h 1 Ω ε ess (t) .
Here above, 1 A denotes the characteristic function of a set A ⊂ Ω.

For later use, it is convenient to divide the residual set Ω ε res (t) further: we define

Ω ε res,B (t) := x ∈ Ω ε res (t) 0 ≤ ε (t, x) < b * 2 , Ω ε res,U B (t) := x ∈ Ω ε res (t) ε (t, x) > 2 b *
as the regions where ε , respectively, stays bounded and may become unbounded.

After this preparation, let us come back to [START_REF] Jüngel | An asymptotic limit of a Navier-Stokes system with capillary effects[END_REF] and derive uniform bounds for the family ε ε . In the essential set, we can perform Taylor's expansion of the function H; we thus get

H( ) -H(b) -H (b) ( -b) ess ≥ c ε -b 2 1 Ω ε ess (t) , which implies that sup t∈R + 1 ε [ ε -b] ess L 2 (Ω) ≤ C . (20)
On the other hand, using the convexity of the function H, the fact that ε -b res ≥ b * /2 and equation ( 9), we discover (see e.g. [START_REF] Feireisl | Sun: Suitable weak solutions to the Navier-Stokes equations of compressible viscous fluids[END_REF] for details) the following bounds on the residual set:

sup t∈R + [ ε ] res γ L γ (Ω) + sup t∈R + -1 ε res κ L κ (Ω) + sup t∈R + [1] res L 1 (Ω) ≤ C ε 2 . ( 21
)
The previous estimate immediately implies that

sup t∈R + L Ω ε res (t) ≤ C ε 2 , (22)
where we have denoted by L(A) the Lebesgue measure of a set A ⊂ Ω. At this point, let us define the quantity

φ ε = 1 ε ε -b .
From the uniform bound (20), we may deduce that

φ ε ess L ∞ (R + ;L 2 ) ≤ C . ( 23
)
On the other hand, using [START_REF] Levermore | Global well-posedness for models of shallow water in a basin with a varying bottom[END_REF] we can compute that for any p ≤ γ, we have

Ω φ ε res p dx ≤ C ε p Ω ε res p dx + Ω [1] res dx ≤ C ε 2-p .
Thus, we finally infer that

∀ 1 ≤ p ≤ γ , φ ε res L ∞ (R + ;L p ) ≤ C ε (2-p)/p . ( 24
)
Of course, this estimate will be useful only in the case when p satisfies the additional restriction p ≤ 2.

The Bresch-Desjardins estimate

As it is apparent from the bounds of the previous subsection, the difficulty with system (1) is that we lose any control on the velocity fields u ε and their gradient ∇u ε close to vacuum, specifically in the region Ω ε res,B . The cold pressure term p c is of great help in order to bypass that difficulty. However, the cold pressure term alone is not strong enough to give us all the pieces of information we need to pass to the limit. On the other hand, system (1) possesses a nice underlying structure, as evidenced for the first time by Bresch and Desjardins (see e.g. [START_REF] Bresch | On some compressible fluid models: Korteweg, lubrication, and shallow water systems[END_REF], [START_REF] Bresch | Existence of global weak solution for a 2D viscous shallow water equation and convergence to the quasi-geostrophic model[END_REF]). By taking advantage of that structure, it is possible to derive, via the so-called BD entropy estimates, some uniform estimates on the gradient of the density functions ε . This is the goal of the next lemma.

A bound coming from BD entropy estimates has been required in the definition of weak solutions, see point (4) in Definition 2.1. Here we show that such a bound holds uniformly with respect to the small parameter ε ∈ ]0, 1]. Note that the result in Lemma 3.1 is stated for smooth solutions to the Navier-Stokes system (1). This is solely to justify the manipulations required to derive the inequality. Once the inequality is proven for the smooth solutions, it is possible to deduce that it is inherited also by the finite energy weak solutions considered in this paper (see [START_REF] Zatorska | On the flow of chemically reacting gaseous mixture[END_REF] and [START_REF] Mucha | Approximate solutions to a model of two-component reactive flow Discrete Contin[END_REF] for details).

In the next statement, we resort to the notation introduced in ( 11), and we recall that we denote by Au = ∇u -∇ t u /2 the skew-symmetric part of the Jacobian matrix of the vector field u. Lemma 3.1 Let ( ε , u ε ) be the smooth solution to (1)-( 3). Then we have the inequality

sup t∈ ]0,T [ F ε , u ε b (t) + 1 ε 2 Ω H( ε ) -H(b) -H (b)( ε -b) dx + ν ε 2 T 0 Ω b 2 p ( ε ) ε ∇ ε b 2 dx dt + ν T 0 Ω ε |Au ε | 2 dx dt ≤ C , (25)
where the constant C > 0 depends only on the initial data and on T . In particular, the previous bound is uniform with respect to ε ∈ ]0, 1].

Proof. The proof of this estimate follows closely [START_REF] Bresch | An example of low Mach (Froude) number effects for compressible flows with nonconstant density (height) limit[END_REF], with the only modifications associated with more general forms of the pressure and of the force. We proceed in several steps, assuming that ε , u ε are smooth enough to justify all the computations. For notational simplicity, in what follows we write ( , u) instead of ε , u ε .

Step 1. From Lemma 5.1 in [START_REF] Bresch | An example of low Mach (Froude) number effects for compressible flows with nonconstant density (height) limit[END_REF] it follows that for sufficiently smooth solutions of the continuity equation in (1), the following equality holds true

1 2 d dt Ω ∇ ln b 2 dx + Ω ∇u • ∇ ln b ∇ ln b dx + Ω ∇u • ∇ ln b∇ ln b dx + Ω u • ∇∇ ln b∇ ln b dx + Ω ∇ div u∇ ln b dx = 0.
Step 2. In this step, one multiplies the momentum equation by ν∇ ln b , we get

ν Ω (∂ t u + u • ∇u) • ∇ ln b dx + ν 2 Ω Du : ∇ 2 ln b dx + ν ε 2 Ω (∇p( ) -∇G) • ∇ ln b dx = 0.
Let us now rewrite each term from the above expression. First, note that the transport term gives

ν Ω (∂ t u + u • ∇u) • ∇ ln b dx = ν d dt Ω bu • ∇ b dx -ν Ω ∇u : ∇ t u dx + ν Ω u • ∇(∇ ln b)u dx.
For the diffusion term, after noticing that ∇ 2 is always symmetric, we have

ν 2 Ω Du : ∇ 2 ln b dx = ν 2 Ω ∇u : ∇ 2 ln b dx = ν 2 Ω b∇u : ∇ 2 b dx -ν 2 Ω ∇u • ∇ ln b ∇ ln b dx = -ν 2 Ω ∇u • ∇ ln b∇ ln b dx -ν 2 Ω ∇ div u∇ ln b dx -ν 2 Ω ∇u • ∇ ln b ∇ ln b dx.
Finally, for the pressure and force term, using (4), we obtain

ν Ω ∇p( ) -∇G ε 2 • ∇ ln b dx = ν ε 2 Ω b 2 p ( ) ∇ b 2 dx + ν ε 2 Ω p ( ) -p (b) ∇b • ∇ b dx.
Step 3. Now, we sum up the equalities from the previous steps along with the energy estimate. In our case, after setting

H( ; b) := H( ) -H (b)( -b) -H(b)
, the statement of Lemma 5.2 from [START_REF] Bresch | An example of low Mach (Froude) number effects for compressible flows with nonconstant density (height) limit[END_REF] gives

d dt Ω F , u b + 1 ε 2 H( ; b) dx + ν ε 2 Ω b 2 p ( ) ∇ b 2 dx + ν Ω |Du| 2 dx = ν Ω ∇u : ∇ t u dx -ν Ω u • (∇∇ ln b)u dx -ν 2 Ω u∇∇ ln b∇ ln b dx - ν ε 2 Ω p ( ) -p (b) ∇b • ∇ b dx . (26)
Now notice that ν

Ω |Du| 2 dx -ν Ω ∇u : ∇ t u dx = ν Ω |Au| 2 dx,
and so, we finally get

d dt Ω F , u b + 1 ε 2 H( ; b) dx + ν ε 2 Ω b 2 p ( ) ∇ b 2 dx + ν Ω |Au| 2 dx = -ν Ω u • (∇∇ ln b)u dx -ν 2 Ω u∇∇ ln b∇ ln b dx - ν ε 2 Ω p ( ) -p (b) ∇b • ∇ b dx = 3 i=1 J i . (27)
Step 4. We have to control the terms J 1 , J 2 , J 3 appearing in the right-hand side of ( 27). To begin with, we can easily estimate

T 0 J 1 dt = -ν T 0 Ω u • (∇∇ ln b)u dx dt ≤ ∇ 2 ln b L ∞ T (L ∞ ) T 0 Ω |u| 2 dx dt ≤ C . For J 2 , instead, we get T 0 J 2 dt = -ν 2 T 0 Ω u∇∇ ln b∇ ln h b dx dt ≤ ν 2 ∇ 2 ln b L ∞ T (L ∞ ) T 0 F , u b dt ,
hence J 2 can be controlled by means of a Grönwall argument. Finally, we have to deal with J 3 . This estimate is a bit more involved than the previous ones, as this term depends on the pressure. We have to distinguish some different cases.

Case 1: integral over Ω ε ess . We start by bounding the part of J 3 which is restricted to the essential set. For this, we use the Taylor expansion p ( ) = p (b) + p (z)( -b), for some z between and b, and the fact that is bounded in Ω ess to write

ν ε 2 Ω ε ess p ( ) -p (b) ∇b • ∇ b dx ν ε 2 Ω ε ess | -b| |∇b| p ( ) b ∇ b dx C δ φ ess 2 L 2 + δ ν ε 2 Ω b 2 p ( ) ∇ b 2 dx ,
where we have also applied the Young inequality in the last step. Here, δ > 0 can be taken arbitrarily small, to the price of increasing the value of C δ . Thus, taking δ small enough, after integration in time, we can absorbe the last term of the previous estimate in the left-hand side of ( 27), while the other term is bounded by a uniform constant C times T , in view of [START_REF] Lions | Incompressible limit for a viscous compressible fluid[END_REF]. Case 2: integral over Ω ε res . In the residual set, it is convenient to split J 3 into two pieces, as follows:

ν ε 2 Ω ε res p ( ) -p (b) ∇b • ∇ b dx = ν ε 2 Ω ε res p ( )∇b • ∇ b dx - ν ε 2 Ω ε res p (b)∇b • ∇ b dx (28)
The bound of the first term in the last equality is easy. We start by writing

ν ε 2 Ω ε res p ( )∇b • ∇ b dx = ν ε 2 Ω ε res p ( ) ∇ ln b • p ( ) b ∇ b dx.
Then, we observe that

p ( ) = γ + -κ =⇒ p ( ) res ∈ L ∞ R + ; L 2
uniformly in ε > 0, owing to [START_REF] Levermore | Global well-posedness for models of shallow water in a basin with a varying bottom[END_REF]. Indeed, one has p ( )

res L 2 ≤ γ/2 L γ + -1 κ/2
L κ . Thus, using (21) quantitatively and arguing as in Case 1, we deduce that

ν ε 2 Ω ε res p ( ) ∇ ln b • p ( ) b ∇ b dx C δ ε 2 p ( ) res 2 L 2 + δ ν ε 2 Ω b 2 p ( ) ∇ b 2 dx C δ + δ ν ε 2 Ω b 2 p ( ) ∇ b 2 dx .
Once again, after integration in time, the last term on the right can be absorbed in the left-hand side of ( 27), for δ > 0 small enough. It remains to deal with the last term appearing in [START_REF] Tsai | Lectures on Navier-Stokes equations[END_REF]. We start by writing

ν ε 2 Ω ε res p (b)∇b • ∇ b dx = ν ε 2 Ω ε res p ( ) p (b) b ∇b • p ( ) b ∇ b dx . First note that p ( ) = p E ( ) + p c ( ) = κ+2 γ+κ + 1 .
Thus, on the one hand we get

ν ε 2 Ω ε res,B p ( ) p (b) b ∇b • p ( ) b ∇ b dx ν ε 2 L(Ω ε res ) 1/2 p ( ) b ∇ b L 2 .
On the other hand, we also have

1 Ω ε res,U B p ( ) 1 Ω ε res,U B (2-γ)/2 ε . Now, in view of (21), for 2 -γ ≥ 0 we have that [ ] (2-γ)/2 res ∈ L ∞ (R + ; L q ) , q := 2γ 2 -γ > 2 , with [ ] (2-γ)/2 res L q = [ ] res γ/q L γ ε 2/q .
If, on the other hand 2 -γ < 0, then we have

[ ] (2-γ)/2 res ∈ L ∞ (R + ; L q ) , q := 2κ γ -2 , with [ ] (2-γ)/2 res L q = -1 res κ/q L κ ε 2/q .
In the latter case, in order to have q ≥ 2, we need the condition

κ ≥ γ -2 .
In any case, and under the assumption that the previous requirement is satisfied in the case γ > 2, we can thus bound

ν ε 2 Ω ε res,U B p ( ) p (b) b ∇b • p ( ) b ∇ b dx 1 ε 2 [ ] (2-γ)/2 res L q p ( ) b ∇ b L 2 L(Ω ε res ) 1/m
, where 1/q + 1/m = 1/2. In the end, using ( 21), [START_REF] Lions | Mathematical topics in Fluid Mechanics[END_REF] and the definition of m, after an application of Young's inequality we arrive at

ν ε 2 Ω ε res p ( ) p (b) b ∇b • p ( ) b ∇ b dx C δ + δ ν ε 2 p ( ) b ∇ b 2 L 2 .
Finally, putting all those inequalities together, we see that

T 0 J 3 dt C δ + 3 δ ν ε 2 T 0 Ω b 2 p ( ) ∇ b 2 dx dt .
Hence, for δ > 0 small enough, we can absorbe the last term of the previous inequality into the left-hand side of [START_REF] Pedlosky | Geophysical fluid dynamics[END_REF]. The proof of the lemma is thus completed. 2

Before moving on, let us observe that the constant in Lemma 3.1 depends on time. Therefore, all the bounds which we are going to derive from the BD entropy estimate will be only local in time, on any arbitrarily large but compact interval [0, T ]. This contrasts with the bounds coming from the classical energy inequality, which instead are global in time.

Consequences of the BD entropy estimate

Lemma 3.1 provides us with at least three additional pieces of information with respect to the classical energy inequality [START_REF] Jüngel | An asymptotic limit of a Navier-Stokes system with capillary effects[END_REF]. Those additional bounds will be fundamental in order to prove convergence in our setting.

First of all, estimate (25) can be used to improve [START_REF] Klainerman | Compressible and incompressible fluids[END_REF]: indeed, as ∇u = Du + Au, we have

√ ε ∇u ε ε ⊂ L 2 loc R + ; L 2 (Ω) . (29)
Next, combining [START_REF] Mucha | Approximate solutions to a model of two-component reactive flow Discrete Contin[END_REF] with [START_REF] Jüngel | An asymptotic limit of a Navier-Stokes system with capillary effects[END_REF], we gather that

√ ε ∇ ln ε b = b √ ε ∇ ε b ∈ L ∞ loc R + ; L 2 (Ω) , (30)
with uniform inclusion (of course, uniform with respect to ε ∈ ]0, 1]). At this point, an easy computation shows that [START_REF] Levermore | Global well-posedness for models of shallow water in a basin with a varying bottom[END_REF]. Putting this information together with [START_REF] Zatorska | On the flow of chemically reacting gaseous mixture[END_REF], by Sobolev embeddings we finally deduce that

b √ ε ∇ ε b = 2 ∇ √ ε - √ ε ∇ ln b = 2 ∇ √ ε - [ ε ] ess + [ ε ] res ∇ ln b . Observe that [ ε ] ess ε is uniformly bounded in L ∞ (R + × Ω), whereas [ ε ] res ε is uniformly bounded in L ∞ (R + ; L 2γ ) → L ∞ (R + ; L 2 ), in view of
(∇ √ ε ) ε ⊂ L ∞ loc R + ; L 2 (Ω) =⇒ ε ε ⊂ L ∞ loc R + ; L 3 (Ω) . (31)
The previous property provides higher integrability for ε , globally (i.e. in the whole Ω, without having to distinguish between essential and residual parts) and uniformly with respect to ε ∈ ]0, 1].

Next, we consider the quantity

1 ε 2 p ( ε ) ε b 2 ∇ ε b 2 ε ⊂ L 1 loc R + ; L 1 (Ω) . (32)
The term inside the parentheses can be written as

p ( ε ) ε b 2 ∇ ε b 2 = γ-2 ε + -κ-2 ε b 2 ∇ ε b 2 = b γ/2 ε b (γ-2)/2 ∇ ε b 2 + b -κ/2 ε b (-κ-2)/2 ∇ ε b 2 = 2 γ b γ/2 ∇ ε b γ/2 2 + - 2 κ b -κ/2 ∇ ε b -κ/2 2
.

Thus, we get the following uniform embeddings:

1 ε ∇ ε b γ/2 ε ⊂ L 2 loc R + ; L 2 (Ω) , 1 ε ∇ ε b -κ/2 ε ⊂ L 2 loc R + ; L 2 (Ω) . (33)
This is another fundamental piece of information, since it gives, roughly speaking, a uniform control on the gradient of ε /b, with a quantitative bound O(ε) in a suitable norm. However, in order to be able to fully exploit this information, some preparatory work is needed.

To begin with, we write

1 ε ∇ ε b γ/2 = 1 ε ∇ γ/2 ε -b γ/2 b γ/2 .
Therefore, the uniform bound (33) implies that

1 ε γ/2 ε -b γ/2 ε ⊂ L 2 loc R + ; L 6 (Ω) . (34)
Let us derive a couple of properties from (34). We notice that we can write

γ/2 ε -b γ/2 = γ ε -b γ γ/2 ε + b γ/2
.

On the one hand, from the previous relation we immediately deduce that

1 ε 1 Ω ε res,B ∪Ω ε ess ( γ ε -b γ ) ε ⊂ L 2 loc R + ; L 6 (Ω) .
We use a Taylor expansion of the function f (s) = s γ to get, for a suitable point z ε = z ε (t, x) between ε (t, x) and b(x), the following fact:

γ ε -b γ = γ b γ-1 ε -b + γ (γ -1) z γ-2 ε ε -b 2 ≥ γ b γ-1 ε -b .
This inequality, together with (34) above, implies that

1 Ω ε res,B ∪Ω ε ess φ ε ε ⊂ L 2 loc R + ; L 6 (Ω) . ( 35
)
On the other hand, on the set Ω ε res,U B we have that ε ≥ 2 b * . Hence, on that set we can write

γ ε -b γ γ/2 ε + b γ/2 ≥ ε 2 γ 2 γ -1 2 γ/2 ε ≥ C γ/2 ε .
Using (34) again, we gather that

1 ε 1 Ω ε res,U B γ/2 ε ε ⊂ L 2 loc R + ; L 6 (Ω) =⇒ 1 ε 1 Ω ε res,U B ε ε ⊂ L γ loc R + ; L 3γ (Ω) . ( 36
)
The main point of the last computations is that (33) does not really provide a uniform control on the gradient of the functions φ ε in L 2 T (L 2 ), whenever γ = 2. Such a control, which is true when γ = 2, would give higher integrability of the φ ε 's in L 2 T (L 6 ). Estimate (35) shows that we are not too far from getting that property, but this fact is true only when ε stays bounded. Since we do not have L ∞ bounds on ε , the density functions may grow in some part of the domain Ω (which has to be small, recall (21) above). Anyway, inequality (36) provides us with a useful control in that region.

Actually, in Ω ε res,B we can derive an even better control on the ε 's than (35). Indeed, arguing in a similar way as for getting (36), from the second piece of information in (33) we deduce

1 ε 1 Ω ε res,B -κ/2 ε ε ⊂ L 2 loc R + ; L 6 (Ω) =⇒ 1 ε 1 Ω ε res,B -1 ε ε ⊂ L κ loc R + ; L 3κ (Ω) . ( 37
)
We conclude this part by showing a third consequence of the BD entropy estimate [START_REF] Mucha | Approximate solutions to a model of two-component reactive flow Discrete Contin[END_REF], that is a control on the gradient functions ∇φ ε in suitable Lebesgue norms. Since this is a key property, we put it in the form of a proposition. Proposition 3.2 Let γ > 1 and κ > 0 such that κ ≥ γ -2. Let ε , u ε ε be a family of global in time finite energy weak solutions to system (1), in the sense of Definition 2.1, related to initial data 0,ε , u 0,ε ε such that (15) holds.

Then:

• if γ ≥ 2, one has the uniform embedding φ ε ε ⊂ L 2 loc R + ; H 1 (Ω) ; • when γ < 2, one deduces instead φ ε ε ⊂ L 2 loc R + ; W 1,γ (Ω) .
Proof. We start by observing that, as a consequence of (32), we get in particular

ε -(1-γ/2) b ∇ φ ε b ε ⊂ L 2 loc R + ; L 2 (Ω) , ε -(1+κ/2) b ∇ φ ε b ε ⊂ L 2 loc R + ; L 2 (Ω) . (38)
Since on the set Ω ε ess ∪ Ω ε res,B we have 0 ≤ ε ≤ 2b * , we easily get

1 Ω ε ess ∪ Ω ε res,B b ∇ φ ε b 1 Ω ε ess ∪ Ω ε res,B ε -(1+κ/2) b ∇ φ ε b ,
which implies that, for all T > 0 fixed, one has

1 Ω ε ess ∪ Ω ε res,B b ∇ φ ε b ε ⊂ L 2 T (L 2 ) . ( 39
)
On the subset Ω ε res,U B , we need to consider two cases. Case 1. For 1 -γ/2 < 0 we can employ the first part of (38) and an argument similar to the one used above to deduce that

∀ T > 0 , 1 Ω ε res,U B b ∇ φ ε b ε ⊂ L 2 T (L 2 ) .
Case 2. For 1 -γ/2 ≥ 0 instead, we can write

1 Ω ε res,U B b ∇ φ ε b = 1 Ω ε res,U B ε -(1-γ/2) b ∇ φ ε b ε 1-γ/2 res,U B .
Since, for any T > 0 fixed,

ε res ε is uniformly bounded in L ∞ T (L γ ), we have that ε 1-γ/2 res,U B ε is uniformly bounded in L ∞ T (L 2γ/(2-γ) ).
In turn, combining this information with (38) implies that

∀ T > 0 , 1 Ω ε res,U B b ∇ φ ε b ε ⊂ L 2 T (L γ ) .
We have thus discovered that the family of ∇ φ ε /b 's is uniformly bounded in the space L 2 T (L γ ), for any time T > 0.

In order to conclude, we simply write

∇ φ ε b = 1 b ∇φ ε - 1 b 2 φ ε ∇b ,
and, when γ ≤ 2, we bound φ ε ε uniformly in L ∞ T (L γ ), thanks to [START_REF] Lions | Incompressible limit for a viscous compressible fluid[END_REF] and [START_REF] Masmoudi | Rigorous derivation of the anelastic approximation[END_REF]. When γ > 2, instead, we use [START_REF] Lions | Incompressible limit for a viscous compressible fluid[END_REF] again, together with the fact that, by definition of φ ε , we have

Ω res,B |φ ε | 2 dx = 1 ε 2 Ω res,B | ε -b| 2 dx 1 ε 2 L Ω res 1 Ω res,U B |φ ε | 2 dx = 1 ε 2 Ω res,U B | ε -b| 2 dx 1 ε 2 Ω res,U B 2 ε dx 1 ε 2 ε res 2 L γ L Ω res (γ-2)/γ 1 ,
where we have used also the bounds provided by [START_REF] Levermore | Global well-posedness for models of shallow water in a basin with a varying bottom[END_REF]. We thus conclude that, when γ > 2, we have

φ ε ε ⊂ L ∞ T (L 2
) for any T > 0 fixed. In the end, the proposition is proved.

2

We also notice that, as a consequence of Proposition 3.2 and Sobolev embedding, we get that

if 1 < γ < 2 , φ ε ε ⊂ L 2 loc R + ; L 3γ/(3-γ) (Ω) , and if γ ≥ 2 , φ ε ε ⊂ L 2 loc R + ; L 6 (Ω) . ( 40
)
4 The singular perturbation operator and the wave system

In this section, we study in detail the singular part of the primitive equations (1). To begin with, in Subsection 4.1 we establish first convergence properties for the sequence of solutions ε , u ε ε towards some targe profile , u . Then, in Subsection 4.2 we derive some constraints that the limit point , u has to satisfy. Finally, in Subsection 4.3 we introduce the system of waves, which encodes the propagation of fast time oscillations.

Preliminary convergence properties

From the uniform bounds exhibited in Subsections 3.1, 3.2 and 3.3, we can derive first convergence properties for the family of solutions ε , u ε ε of our primitive system [START_REF] Almgren | Low Mach number modeling of type Ia supernovae[END_REF]. Of course, the convergence we are going to establish is only in weak topologies, therefore it will not be enough for deriving the limit system [START_REF] Danchin | Zero Mach number limit for compressible flows with periodic boundary conditions[END_REF].

Here below, it is convenient to work on time intervals [0, T ], for arbitrary large but fixed T > 0. Also, in the notation we imply that all the convergences are taken in the limit ε → 0 + . The density functions. We start by considering the sequence of the density functions ε . First of all, from ( 20) and ( 21), we see that we can decompose

ε = b + (1) ε + (2) ε , with (1) ε
-→ 0 in L ∞ T (L 2 ) , (2) ε
-→ 0 in L ∞ T (L γ ) . ( 41
)
Since Ω is of finite measure, we can interpolate those convergence properties with (31) to deduce that

ε -b = r ε , with r ε -→ 0 in L ∞ T (L p ) ∀ p < 3 .
The uniform bounds of Section 3 allow us to find more quantitative convergence properties. Indeed, from [START_REF] Lions | Incompressible limit for a viscous compressible fluid[END_REF] and (35) it follows that there exists a function φ ∈ L ∞ R + ; L 2 (Ω) ∩ L 2 loc R + ; L 6 (Ω) such that, up to the extraction of a suitable subsequence, one has

φ ε ess * φ in L ∞ T (L 2 ) ∩ L 2 T (L 6 ) , (42)
where the symbol * stands for the weak- * convergence in the respective functional space. On the other hand, owing to [START_REF] Masmoudi | Rigorous derivation of the anelastic approximation[END_REF], we know that

φ ε res -→ 0 in L ∞ T (L p ) , ∀ p ∈ [1, min{γ, 2}[ . (43)
Of course, if γ ≥ 2, we have the previous strong convergence only for all 1 ≤ p < 2; however, interpolating with the bounds of (40), in that case we get

φ ε res -→ 0 in L 2 T (L p ) , ∀ p ∈ [1, 6[ .

The velocity fields.

As it is apparent from equations (1), any information on the velocity fields u ε and their gradients is lost in regions close to vacuum. This is one of the main difficulties arising in the analysis of system (1). On the other hand, at least at a first sight, it is not so clear which is the right quantity to look at; for instance, keep in mind inequalities [START_REF] Klainerman | Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF], [START_REF] Vasseur | Existence of global weak solutions for 3D degenerate compressible Navier-Stokes equations[END_REF]. Here, we decide to work with the momentum

V ε := ε u ε .
However, the first step is to get some uniform bounds on the velocity fields u ε . This is the goal of the next proposition.

Proposition 4.1 Let ε , u ε ε be a family of global in time finite energy weak solutions to system (1), in the sense of Definition 2.1, related to initial data 0,ε , u 0,ε ε such that (15) holds.

Then we have the uniform estimates

u ε ε ⊂ L ∞ R + ; L p 1 (Ω) and ∇u ε ε ⊂ L 2 loc R + ; L p 1 (Ω) , p 1 := 2κ κ + 1 .
In particular, we also have that

u ε ε ⊂ L 2 loc R + ; L p 2 (Ω)
, where p 2 := 6κ κ + 3 .

Remark 4.2 Note that, due to our assumption κ > 3, one has p 1 > 3/2 and p 2 > 3 in the above proposition.

Proof. We start by writing

u ε = -1/2 ε √ ε u ε , from
which, by use of ( 18) and ( 21), we immediately deduce that u ε ε is bounded in L ∞ T (L p 1 ). For the gradient terms, the argument is analogous, since we can write

∇u ε = -1/2 ε

√

ε ∇u ε and use, this time, ( 29) and ( 21). This implies the claimed bound on the sequence ∇u ε ε .

Finally, the last uniform bound for the family of u ε 's follows from the previous property and Sobolev embeddings.

2

We are now ready to derive some important estimates for the velocity fields V ε .

Proposition 4.3 Let ε , u ε ε be a family of global in time finite energy weak solutions to system (1), in the sense of Definition 2.1, related to initial data 0,ε , m 0,ε ε such that (15) holds. For all ε ∈ ]0, 1], define V ε := ε u ε . Then, the following facts hold:

(i) the sequence V ε ε is uniformly bounded in the space L ∞ loc R + ; L 3/2 (Ω) ∩ L 2 loc R + ; W 1,p 3 (Ω) , where p 3 := 6κ/(5κ + 3);

(ii) there exist sequences V ε ε and W ε ε of vector fields such that

∀ ε ∈ ]0, 1] , V ε = V ε + ε W ε,M ,
with the uniform embedding properties

V ε ε ⊂ L ∞ loc R + ; L 2 (Ω) and W ε ε ⊂ L 2 loc R + ; L 3/2 (Ω) ;
(iii) after setting V ε := b u ε , we can also write

V ε = V ε + ε W ε ,
with the uniform embedding properties

V ε ε ⊂ L 2 loc R + ; W 1,p 1 (Ω) and W ε ε ⊂ L 2 loc R + ; L 1 (Ω) .
Proof. Let T > 0 be arbitrary, but fixed throughout the following computations. The proof of the L ∞ T (L 3/2 ) bound of item (i) is easy to get: it is enough to write V ε = √ ε √ ε u ε and use the uniform bounds given in ( 18) and (31). Next, let us focus on the bounds for the gradient. For any 1 ≤ j ≤ 3, we compute

∂ j V ε = √ ε √ ε ∂ j u ε + u ε ∂ j ε = √ ε √ ε ∂ j u ε + 2 √ ε u ε ∂ j √ ε := A ε + B ε .
Repeating the previous argument, using this time ( 29) and (31), it is easy to see that A ε ε ⊂ L 2 T (L 3/2 ). As for B ε , we start by observing that, owing to Proposition 4.1, we have u ε ε ⊂ L 2 T (L p ) for p ≤ p 2 := 6κ κ+3 . On the other hand, we have

∇ √ ε ε ⊂ L ∞ T (L 2 ) and √ ε ε ⊂ L ∞ T (L 6
). Hence, we get B ε ε ⊂ L 2 T (L p 3 ) where p 3 = 6κ 5κ+3 . Notice that p 3 > 1 if and only if κ > 3: this is precisely the place where the strongest assumption on κ appears. Item (i) is thus proven.

For showing the decomposition in item (ii), we notice that

V ε = √ b √ ε u ε + ε √ ε - √ b ε √ ε u ε .
Thus, if we define

V ε := √ b √ ε u ε , from (18) we immediately infer that V ε ε ⊂ L ∞ T (L 2
). Next, we define

W ε := √ ε - √ b ε √ ε u ε .
At this point, we use that

√ ε - √ b ≤ b -1/2 * | ε -b|.
Hence, on the one hand, by using (35), we see that

1 Ω ε res,B ∪Ω ε ess √ ε - √ b ε ε ⊂ L 2 T (L 6 ) ;
on the other hand, since

0 ≤ 1 Ω ε res,U B √ ε - √ b ≤ 1 Ω ε res,U B √ ε 1 Ω ε res,U B γ/2 ε ,
in view of (36) we gather that also

1 Ω ε res,U B √ ε - √ b ε ε is uniformly bounded in L 2 T (L 6 ). Therefore, we get that √ ε - √ b ε ε ⊂ L 2 T (L 6 ) =⇒ W ε ε ⊂ L 2 T (L 3/2 ) .
The proof of item (iii) is similar. This time, we write

ε -b = √ ε - √ b √ ε + √ b and get V ε = b u ε + √ ε - √ b √ ε u ε + √ ε - √ b √ b u ε .
If we set V ε := b u ε , Proposition 4.1 ensures us that the claimed bounds for this quantity are satisfied. Next, we claim that the sequence of

W ε := √ ε - √ b ε √ ε u ε + √ ε - √ b ε √ b u ε is uniformly bounded in L 2 T (L 1
), for all T > 0 fixed. For this, owing to [START_REF] Klainerman | Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF], the last part of Proposition 4.1 and Remark 4.2, it is enough to show that

√ ε - √ b is of order O(ε) in L ∞ T (L 2 ) . ( 44
)
As a matter of fact, since 20) and ( 21) we immediately get that

√ ε - √ b ≤ b -1/2 * | ε -b|, from (
1 Ω ε ess + 1 Ω ε res,B √ ε- √ b ε ε is uniformly bounded in L ∞ T (L 2
). Finally, as done above, we have

0 ≤ 1 Ω ε res,U B √ ε - √ b ≤ 1 Ω ε res,U B √ ε ,
which, in view of ( 21) again, implies that

1 Ω ε res,U B √ ε - √ b L 2 1 Ω ε res,U B √ ε L 2γ L(Ω ε res ) 1/2-1/2γ
ε .

In the end, we have shown (44), which in turn implies the sought bound for the vector fields W ε . 2

From the previous proposition, we immediately deduce the following corollary. The proof is rather straightforward, hence omitted.

Corollary 4.4 Under the assumptions of Proposition 4.3, there exists a vector field V, belonging to

L ∞ loc R + ; L 3/2 (Ω) ∩ L 2 loc R + ; W 1,p 3
(Ω) such that, up to the extraction of a suitable subsequence, one has V ε * V in that space.

In addition, V also belongs to

L ∞ loc R + ; L 2 (Ω) ∩ L 2 loc R + ; W 1,p 1 (Ω)
, and, up to further extractions, for all T > 0 one has that

V ε * V in L ∞ T (L 2 ) and V ε V in L 2 T (W 1,p 1
). To conclude this part, we define the target velocity field U as

U := 1 b V , ( 45
)
where V is the vector field identified in Corollary 4.4.

We notice that U is (up to a further extraction) the weak-limit point of the sequence u ε ε in the functional spaces identified in Proposition 4.1. We point out also that U belongs to the same functional spaces to which V belongs.

Constraints on the limit

In the previous subsection, we have identified the limit points b, φ and V of the families ε ε , φ ε ε and V ε ε , respectively. Our next goal is to find some properties those limit points have to satisfy. We point out that these conditions do not fully characterise the limit dynamics, which will be deduced by passing to the limit in the momentum equation. Proposition 4.5 Let ε , u ε ε be a family of global in time finite energy weak solutions to system (1), in the sense of Definition 2.1, related to initial data 0,ε , u 0,ε ε such that (15) holds. For all ε ∈ ]0, 1], define V ε := ε u ε . Let φ be the scalar function identified in (42), V the vector field identified in Corollary 4.4 and U the vector field defined in (45).

Then U has to satisfy the anelastic constraint

div b U = 0 in D R + × Ω .
On the other hand, φ = φ(b) is determined as a function of b only.

Before proving the previous proposition, let us state a simple lemma. It will be useful to understand how to deal with the singularity of the pressure and gravitational terms, and to put in evidence the right singularity in the momentum equation.

1 ε 2 ∇p( ε ) -ε ∇G = 1 ε b ∇ H (b) φ ε + 1 ε 2 ∇Π( ε ; b) .
In addition, one has that 1

ε 2 Π( ε ; b) ε is uniformly bounded in L ∞ R + ; L 1 (Ω) .
Proof. The claimed identity follows from a simple algebraic computation:

1 ε 2 ∇p( ε ) -ε ∇G = 1 ε 2 ∇ p( ε ) -p(b) -ε -b ∇G = 1 ε 2 ∇Π( ε ; b) + 1 ε 2 ∇ p (b) ε -b -ε -b ∇G = 1 ε b ∇ H (b) φ ε + 1 ε 2 ∇Π( ε ; b) ,
where, in the last step, we have used the equality ∇G = H (b) ∇b, which holds owing to (4). It remains us to show the uniform bounds for the family 1 ε 2 Π( ε ; b) ε . The argument is similar to the proof of Lemma 4.1 in [START_REF] Fanelli | Incompressible and fast rotation limit for barotropic Navier-Stokes equations at large Mach numbers[END_REF]; however, since that paper did not deal with non-constant limit density profiles b nor with the presence of a cold component p c of the pressure, we report here most of the details.

First of all, using Taylor formula at the second order for the pressure function p together with ( 23), we easily see that

1 ε 2 [Π( ε ; b)] ess ≈ φ 2 ε =⇒ 1 ε 2 [Π( ε ; b)] ess ε ⊂ L ∞ (R + ; L 1 ) .
Next, let us denote by Π E ( ε ; b) and Π c ( ε ; b) the functions defined as Π( ε ; b), but using respectively p E and p c instead of the full pressure function p. Then, a Taylor expansion again and ( 22) imply that

1 ε 2 Π E ( ε ; b) res,B + Π c ( ε ; b) res,U B ⊂ L ∞ (R + ; L 1 ) .
Finally, we notice that

Π E ( ε ; b) res,U B ε γ res and Π c ( ε ; b) res,B ε -κ res ,
for which one can deduce the sought bounds by using the controls in [START_REF] Levermore | Global well-posedness for models of shallow water in a basin with a varying bottom[END_REF].

In the end, the lemma is completely proved. 2

We can now turn our attention to the proof of Proposition 4.5.

Proof of Proposition 4.5. We start by considering the weak form of the momentum equation: given a test function ϕ ∈ D R + × Ω , with Supp ϕ ⊂ [0, T [ ×Ω for some positive time T > 0, we have

- T 0 Ω ε ∂ t ϕ dx dt - T 0 Ω ε u ε • ∇ϕ dx dt = Ω 0,ε ϕ(0) dx .
By assumption on the initial data, we know that 0,ε -b -→ 0 in L 2 ∩ L ∞ , so it is easy to pass to the limit in the right-hand side of the previous relation. Moreover, in view of (41) and Corollary 4.4, we know that ε -→ b and

ε u ε = V ε -→ V in D ([0, T [ ×Ω)
. Thus, passing to the limit in the previous equality yields the constraint

T 0 Ω V • ∇ϕ dx dt = 0 ∀ ϕ ∈ D =⇒ div V = 0 in D .
Now, using the definition of U given in (45) immediately gives the anelastic constraint.

Let us turn our attention to the momentum equation. As appears in the statement, the momentum equation does not give any relevant information on the limit points: we will discover that φ is a quantity which play no role in the limit dynamics.

Indeed, we can test the momentum equation by ε ψ, where ψ ∈ D(R + × Ω; R 3 ) is a test function with support (say) included in [0, T [ ×Ω, for some T > 0. By the uniform bounds established in Section 3, and using the identity in Lemma 4.6 above, it is possible to see that all the terms in the momentum equation tend to 0, except for the term b ∇ H (b) φ ε . Thus, passing to the limit for ε → 0 + , we find that

T 0 Ω H (b) φ div b ψ dx dt = 0 ∀ ψ ∈ D .
Since b never vanishes, this implies that ∇ H (b) φ = 0, and then H b(x) φ(t, x) = c(t) a.e. on R + × Ω, for a suitable function c(t) only depending on the time variable. We claim that c(t) ≡ c(0) is in fact constant, and does not depend on time neither. Indeed, let us denote by f the space average over

Ω of a function f = f (t, x), namely f = f (t) := 1 L(Ω) Ω f (t, x) dx .
Using the definition of φ ε , we can recast the continuity equation as

∂ t φ ε + 1 ε div V ε = 0 . ( 46
)
Taking the mean value with respect to the space variable, we discover that ∂ t φ = 0, thus the mean value of φ is preserved in time: φ (t) = φ 0 for almost all times t ≥ 0, where φ 0 is the weak limit point of the initial data φ 0,ε ε specified in [START_REF] Ebin | The motion of slightly compressible fluids viewed as a motion with strong constraining force[END_REF]. On the other hand, coming back to the equality H (b) φ = c(t) and computing the space average, we have that φ (t) = 1/H (b) c(t), which in turn implies that c(t) ≡ c 0 has to be a constant function.

The proposition is now proved. 2

As already pointed out, the previous proposition does not allow us to identify the limit dynamics yet. The main problem consists in passing to the limit in the momentum equation, which reveals a not so easy task, owing to the non-linear terms appearing in it. In order to succeed, we first need to understand the propagation of fast oscillations in time: this is the goal of the next subsection.

Acoustic equation

Lemma 4.6 and Proposition 4.5 allow us to identify the singular part of the equations of motion. Because of the ill-preparation of the initial data, this singular part is responsible for fast oscillations in time of the solution, the so called acoustic waves, which may eventually prevent the convergence of the non-linear terms. In order to study those oscillations and be able to pass to the limit in the equations, we reformulate our system (1) as a wave system.

Let us recall that we have denoted

φ ε = ε -b ε and V ε = ε u ε .
As already remarked in the proof of Proposition 4.5, in terms of those quantities the continuity equation can be rewritten as in (46). Concerning the momentum equation, instead, we can take advantage of the identity of Lemma 4.6 to combine the pressure and gravitation terms together.

In the end, we discover that system (1) can be recasted in the following form:

ε ∂ t φ ε + div V ε = 0 ε ∂ t V ε + b ∇ H (b) φ ε = ε f ε , (47)
where we have defined

f ε = ν div( ε Du ε ) -div( ε u ε ⊗ u ε ) - 1 ε 2 ∇Π( ε ; b) . ( 48
)
By the uniform bounds [START_REF] Klainerman | Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF], we get that the sequence

ε u ε ⊗ u ε ε is uniformly bounded in L ∞ T (L 1
), for all times T > 0. On the other hand, by arguing as in the proof of Proposition 4.3, it is easy to see that

ε Du ε ε is uniformly bounded in L 2
T (L 3/2 ). Finally, the term Π( ε ; b) has been dealt with in Lemma 4.6. Therefore, owing to the Sobolev embedding H s (Ω) ⊂ L ∞ (Ω) for any s > 3/2, we get that

∀ T > 0 , ∀ s > 5/2 , f ε ε ⊂ L 2 T (H -s ) . ( 49
)
For later use, it is convenient to introduce a regularised version of the wave system (47). For this, we employ the low frequency cut-off operator S M , with M ∈ N, introduced in relation (67) in the Appendix.

Since S M commutes with all partial derivatives, applying operator S M to (47) yields

ε ∂ t φ ε,M + div V ε,M = 0 ε ∂ t V ε,M + b ∇ H (b) φ ε,M = ε f ε,M + h ε,M , (50)
where we have denoted

φ ε,M = S M (φ ε ), V ε,M = S M (V ε ) and f ε,M = S M (f ε
), and we have set

h ε,M := b, S M ∇ H (b)φ ε + b ∇ H (b), S M φ ε .
Thanks to the uniform bounds for the sequence φ ε ε in L 2 loc R + ; W 1,γ (Ω) if γ < 2, in L 2 loc R + ; H 1 (Ω) when γ ≥ 2 (keep in mind Proposition 3.2 above), standard commutator estimates (see e.g. Lemma 2.97 of [START_REF] Bahouri | Fourier Analysis and Nonlinear Partial Differential Equations[END_REF]) imply that

∀ T > 0 , h ε,M L 2 T (L γ ) + curl h ε,M L 2 T (L γ ) 2 -M , (51)
where we agree that the Lebesgue exponent γ is changed into 2 whenever γ ≥ 2.

We explicitly point out that, in the above estimate (51), the multiplicative constant is uniform with respect to both M ∈ N and ε ∈ ]0, 1], but it may depend on the fixed time T > 0. Notice that, for the uniform bound on curl h ε,M , the gradient structure of the commutator terms play a key role. We remark also that the commutator term h ε,M is much better controlled than the corresponding term in [START_REF] Masmoudi | Rigorous derivation of the anelastic approximation[END_REF], thanks to the uniform bounds provided by the BD entropy estimates.

In addition, as an immediate consequence of (49), we get that

∀ T > 0 , ∀ s ≥ 0 , f ε,M L 2 T (H s ) ≤ C(s, M ) , (52)
where the constants C(s, M ) > 0 depend only on the quantities inside the parentheses. Notice that these constants blow up when M → +∞, but they have finite value for any M ∈ N fixed. We conclude this part by showing uniform bounds on the two sequences φ ε,M ε,M and V ε,M ε,M . As a matter of fact, in view of our computations in Section 5 below, it is important to introduce a fine decomposition of those terms. Lemma 4.7 Let T > 0 be arbitrarily large, but fixed. For any M ∈ N, define the quantities φ ε,M = S M (φ ε ) and V ε,M = S M (V ε ) as above.

(i) Both sequences φ ε,M ε and V ε,M ε are sequences of smooth functions in space, which are uniformly bounded (with respect to ε, but not with respect to M ) in the space L 2 T (H s ), for any s ≥ 0 arbitrarily large.

(ii) For any M ∈ N and any ε ∈ ]0, 1], we can write

φ ε,M = ϕ ε,M + ε (2-γ)/γ π ε,M , ( 53
)
where we have the uniform estimates

sup M ∈N sup ε∈ ]0,1] ϕ ε,M L ∞ T (L 2 )∩L 2 T (L 6 ) ≤ C and ∀ s ≥ 0 , sup ε∈ ]0,1] π ε,M L ∞ T (H s ) ≤ C(s, M ) ,
for suitable positive constants C and C(s, M ), which depend only on T > 0 and on the quantities in the brackets. In the case γ ≥ 2, one can simply take π ε,M ≡ 0.

(iii) For any M ∈ N and ε ∈ ]0, 1], one has

V ε,M = V ε.M + ε W ε,M , (54)
with the uniform bounds

sup M ∈N sup ε∈ ]0,1] V ε,M L ∞ T (L 2 ) ≤ C and ∀ s ≥ 0 , sup ε∈ ]0,1] W ε,M L 2 T (H s ) ≤ C(s, M ) .
(iv) We can also write

V ε,M = V ε.M + ε W ε,M , (55)
with the uniform estimates

sup M ∈N sup ε∈ ]0,1] V ε,M L 2 T (W 1,p 1 ) ≤ C and ∀ s ≥ 0 , sup ε∈ ]0,1] W ε,M L 2 T (H s ) ≤ C(s, M ) .
Proof. The properties claimed in item (i) are an immediate consequence of the bounds for the families φ ε ε and V ε ε , combined with the smoothing effect of the low frequency cut-off operators S M . Keep in mind ( 23), [START_REF] Masmoudi | Rigorous derivation of the anelastic approximation[END_REF] and Propositions 3.2 and 4.3. The decompositions of items (iii) and (iv), together with the corresponding uniform estimates, also follow from Proposition 4.3, simply setting V ε,M := S M (V ε ) and W ε,M = S M (W ε ), and similarly for V ε,M and W ε,M .

Finally, let us prove item (ii). We start by considering the case 1 < γ < 2. In this case, we decompose

φ ε = φ ε ess + φ ε res .
Thus, if we set ϕ ε,M := S M φ ε ess , estimates [START_REF] Lions | Incompressible limit for a viscous compressible fluid[END_REF] and (35) easily imply the uniform boundedness (both with respect to M ∈ N and ε ∈ ]0, 1]) of the sequence ϕ ε,M ε,M . Next, we define

π ε := 1 ε (2-γ)/γ φ ε res and π ε,M := S M (π ε ) ,
and we conclude with the help of [START_REF] Masmoudi | Rigorous derivation of the anelastic approximation[END_REF]. When γ ≥ 2, instead, one has that

φ ε ε is uniformly bounded in L ∞ T (L 2 ) ∩ L 2 T (H 1
), thanks to Proposition 3.2. Hence, in this case one can simply define ϕ ε,M = S M (φ ε ) and π ε,M = 0. This completes the proof of the lemma. 2

Passage to the limit

In this section, we finish the proof of our main theorem, showing the convergence (up to an extraction) of weak solutions to (1) to weak solutions to the target system [START_REF] Danchin | Zero Mach number limit for compressible flows with periodic boundary conditions[END_REF].

The main problem is to pass to the limit in the most non-linear term, namely the convective term in the momentum equation. This convergence will be proved in the following two subsections, by use of a compensated compactness argument. Finally, in Subsection 5.3 we will take care of the other terms, and complete the proof of the convergence.

Before going into the details, let us recall that convergence will be proved for any test function lying in the kernel of the singular perturbation operator, namely (in view of Proposition 4.5) for any test function

ψ ∈ D R + × Ω; R 3 , such that div b ψ = 0 . (56)
Also, it is useful to introduce the following notation: we denote by R ε,M any remainder term, that is any term such that lim

M →+∞ lim sup ε→0 + T 0 R ε,M , ψ D ×D dt = 0 , (57)
for some given time T > 0 and test function ψ ∈ D [0, T [ ×Ω; R 3 taken as above. Similarly, we will use the notation R ε,M for any scalar term such that lim

M →+∞ lim sup ε→0 + R ε,M = 0 . (58)
Typically, we will have

R ε,M = T 0 R ε,M , ψ D ×D dt .
In the next computations, the precise values of R ε,M and R ε,M may change from one line to another.

Approximation of the convective term

Passing to the limit in the convective term is based on a compensated compactness argument, following work [START_REF] Masmoudi | Rigorous derivation of the anelastic approximation[END_REF] of Masmoudi. This argument relies on using algebraically the structure of the wave system (47) and performing direct computations on it. Of course, for that, we need first of all to smooth out all the quantities with respect to the space variables: this is the scope of the next lemma. Notice that the approximation argument here is delicate, due to the degeneracy of the system close to vacuum. The cold pressure term p c will be of great help.

ε→0 + T 0 Ω ε u ε ⊗ u ε : ∇ψ dx dt - T 0 Ω 1 b V ε,M ⊗ V ε,M : ∇ψ dx dt = 0 .
Proof. We start by observing that, in view of Propositions 4.1 and 4.3, we have

V ε ε ⊂ L 2 T (L 3/2 ) and u ε ε ⊂ L 2 T (L p ) for all p ∈ [1, p 2 ]
, where p 2 = 6κ/(κ + 3). We also notice that, under our assumption κ > 3, we have 2/3 + 1/p 2 < 1, i.e. p 2 > 3. Thus, given any M ∈ N, we can write

I := T 0 Ω ε u ε ⊗ u ε : ∇ψ dx dt = T 0 Ω 1 b V ε ⊗ S M (b u ε ) : ∇ψ dx dt + T 0 Ω 1 b V ε ⊗ (Id -S M )(b u ε ) : ∇ψ dx dt ,
where all terms are well-defined, and the last term on the right can be bounded as follows:

T 0 Ω 1 b V ε ⊗ (Id -S M )(b u ε ) : ∇ψ dx dt V ε L 2 T (L 3/2 ) (Id -S M )(b u ε ) L 2 T (L 3 ) .
At this point, taking advantage of Bernstein's inequalities of Lemma A.1 in the Appendix, we can bound

(Id -S M )(b u ε ) L 3 ≤ j≥M ∆ j (b u ε ) L 3 j≥M 2 -j ∆ j ∇(b u ε ) L 3 (59) j≥M 2 j -1 + 3 1 p 1 -1 3 ∆ j (b u ε ) L p 1 ≤ 2 -αM ∇(b u ε ) B 0 p 1 ,∞ j≥0 2 -αj ,
where p 1 = 2κ/(κ + 1) has been defined in Proposition 4.1 and we have set α = (κ -3)/2κ > 0. Owing to the embedding L p → B 0 p,∞ for any p ∈ [1, +∞], we finally get

(Id -S M )(b u ε ) L 3 2 -αM ∇(b u ε ) L p 1 2 -αM .
As a result, the previous computations show that, in the sense of (58), one has

I = T 0 Ω 1 b V ε ⊗ S M (b u ε ) : ∇ψ dx dt + R ε,M .
Next, we use again the uniform bound

u ε ε ⊂ L 2 T (L p 2 ), together with the fact that S M (b u ε ) ε is uniformly bounded (with respect to ε, but not to M ) in the space L 2 T (L ∞ ), to get T 0 Ω 1 b V ε ⊗ S M (b u ε ) : ∇ψ dx dt (60) = T 0 Ω u ε ⊗ S M (b u ε ) : ∇ψ dx dt + ε T 0 Ω 1 b φ ε u ε ⊗ S M (b u ε ) : ∇ψ dx dt = T 0 Ω 1 b S M b u ε ⊗ S M (b u ε ) : ∇ψ dx dt + R ε,M .
Indeed, first note that, due to (40), we can write, for γ < 2, the following estimate:

T 0 Ω 1 b φ ε u ε ⊗ S M (b u ε ) : ∇ψ dx dt dt φ ε L ∞ T (L 3γ/(3-γ) ) u ε L 2 T (L p 2 ) S M (b u ε ) L 2 T (L ∞ ) 1 .
We observe that this makes sense whenever (3 -γ)/3γ + (κ + 3)/6κ ≤ 1, hence for κ ≥ 3γ/(7γ -6).

But, for γ > 1, one always has 3γ/(7γ -6) < 3, therefore the previous estimate is satisfied for all κ > 3.

Notice that, since p 2 > 3, we can repeat the exact same computations also when γ ≥ 2, up to use the right bound from (40). On the other hand, after noticing that 1/p 1 + 1/p 2 < 1 for κ > 3 and arguing in a similar way as above, we can estimate

T 0 Ω 1 b (Id -S M )(b u ε ) ⊗ S M (b u ε ) : ∇ψ dx dt (Id -S M )(b u ε ) L 2 T (L p 1 ) S M (b u ε ) L 2 T (L p 2 ) 2 -M ∇(b u ε ) 2 L p 1 2 -M .
Thus, we have proven that the last equality in (60) holds true. At this point, to conclude the argument we may use the decomposition of Lemma 4.7. Indeed, keeping in mind the definition of the vector fields V ε given in Proposition 4.3 above, we notice that

S M (b u ε ) = V ε,M = V ε,M -ε W ε,M
. Thus, using the bounds collected in item (iv) of Lemma 4.7 we easily see that

T 0 Ω 1 b S M b u ε ⊗ S M (b u ε ) : ∇ψ dx dt = T 0 Ω 1 b V ε,M ⊗ V ε,M : ∇ψ dx dt + R ε,M .
This last equality finally ends the proof of the lemma. 2

Compensated compactness

Owing to Lemma 5.1, it is enough to compute the limit of the approximate convective term

- T 0 Ω 1 b V ε,M ⊗ V ε,M : ∇ψ dx dt = T 0 Ω div 1 b V ε,M ⊗ V ε,M • ψ dx dt ,
for any test function ψ ∈ D [0, T [ ×Ω; R 3 belonging to the kernel of the singular perturbation operator, namely such that div (bψ) = 0. Here above, we have performed an integration by parts, because the vector fields V ε,M are smooth in the space variable.

Preliminary reductions

Now, we compute

div 1 b V ε,M ⊗ V ε,M = 1 b div(V ε,M ) V ε,M + V ε,M • ∇ 1 b V ε,M (61)
= 1 b div(V ε,M ) V ε,M + b 2 ∇ 1 b V ε,M 2 -curl 1 b V ε,M × V ε,M ,
where the symbol × denotes the ususal external product of vectors in R 3 and, for a 3-D vector field U, we have curl U := ∇ × U. Notice that the second term in the last line identically vanishes, whenever tested against a test function ψ as in (56). Thus, this term contributes as a remainder R ε,M to the limit, in the sense of relation (57). Therefore, resorting to the first equation in (50) for dealing with the first term, we can write

div 1 b V ε,M ⊗ V ε,M = - ε b ∂ t φ ε,M V ε,M -curl 1 b V ε,M × V ε,M + R ε,M = ε b φ ε,M ∂ t V ε,M -curl 1 b V ε,M × V ε,M + R ε,M ,
where, in the second step, we have included the total time derivative ε ∂ t φ ε,M V ε,M into the remainder R ε,M . Indeed, the time derivative can be put on the test function and the family φ ε,M V ε,M ε is uniformly bounded in e.g. L 2 T (L 2 ), owing to item (i) of Lemma 4.7. At this point, we use the second equation in (50) to deal with the term presenting the time derivative, and we get

div 1 b V ε,M ⊗ V ε,M = -φ ε,M ∇ H (b) φ ε,M + ε b φ ε,M f ε,M + 1 b φ ε,M h ε,M -curl 1 b V ε,M × V ε,M + R ε,M .
Owing to (52) and item (i) of Lemma 4.7, it is clear that the second term in the right-hand side is a remainder, in the sense of (57). Next, we claim that also the third term, i.e. φ ε,M h ε,M /b ε,M , is a remainder. For proving this claim, in the case γ ≥ 2 it is enough to employ Proposition 3.2 and recall that estimate (51) holds true with the Lebesgue exponent γ replaced by 2. When γ ∈ ]1, 2[ , instead, we notice that the piece of information coming from (40) is not good enough to cover all possible values of γ in that interval. Our approach is instead based on the use of item (ii) of Lemma 4.7, which allows us to write

1 b φ ε,M h ε,M = 1 b ϕ ε,M h ε,M + ε (2-γ)/γ 1 b π ε,M h ε,M .
Thanks to (51) and item (ii) of Lemma 4.7, we can estimate

1 b ϕ ε,M h ε,M L 1 T (L 6γ/(6-γ) ) ϕ ε,M L 2 T (L 6 ) h ε,M L 2 T (L γ ) 2 -M ,
which implies that this term satisfies (57). In addition, it is easy to see that also ε (2-γ)/γ π ε,M h ε,M /b ε verifies (57). Indeed, recall that in the case γ ≥ 2, one can simply take π ε,M ≡ 0. For 1 < γ < 2 one has only to notice that

ε (2-γ)/γ 1 b π ε,M h ε,M L 2 T (L γ ) ε (2-γ)/γ π ε,M L ∞ T (L ∞ ) h ε,M L ∞ T (L γ ) C(M ) ε (2-γ)/γ .
As a result of the previous computations, we infer that

div 1 b V ε,M ⊗ V ε,M = -φ ε,M ∇ H (b) φ ε,M -curl 1 b V ε,M × V ε,M + R ε,M . (62)

Compactness of the rotational part

The next lemma takes care of the convergence of the curl term in (62).

Lemma 5.2 Denote by V the weak-limit of V ε ε identified in Corollary 4.4, so that V M = S M V is the weak-limit of V ε,M ε when ε → 0 + . Then, for any ψ ∈ D R + × Ω; R 3 , one has lim

M →+∞ lim sup ε→0 + T 0 Ω curl 1 b V ε,M × V ε,M • ψ dx, dt - T 0 Ω S M curl 1 b V × V M • ψ dx dt = 0 .
Proof. We start the proof by recalling that, owing to Proposition 4.3, we have

V ε ε ⊂ L 2 T (W 1,p 3 ), hence curl V ε ε ⊂ L 2 T (L p 3 ).
Observe that, by dual Sobolev embeddings (see e.g. Theorem 0.5 of [START_REF] Feireisl | Singular limits in thermodynamics of viscous fluids[END_REF]), we have that, for any compact K ⊂ Ω, the space L p 3 (K) is compactly embedded into H -2 (K), for instance.

Next, consider the (not regularised) wave system (47). Dividing the momentum equation by b and then taking the curl (which simply consists in an adequate choice of the test function in the weak formulation), we deduce

∂ t curl 1 b V ε = curl f ε b .
In turn, this relation, together with (49) and the fact that b ∈ C 3 (Ω), implies that ∂ t curl 1 b V ε ε is uniformly bounded in the space L 2

T (H -s ), for any s > 7/2. Putting together these pieces of information and applying the Aubin-Lions lemma (see e.g. Lemma 3.7 of [START_REF] Tsai | Lectures on Navier-Stokes equations[END_REF]), we gather that

curl 1 b V ε ε is compact in L 2 T H -2 (K) ,
for any compact set K ⊂ Ω. This implies that, for any M ∈ N and any compact K ⊂ Ω fixed, the sequence

S M curl 1 b V ε ε is compact (with respect to ε) in L 2 T L 2 (K) . (63)
Next, we write

T 0 Ω curl 1 b V ε,M × V ε,M • ψ dx dt = T 0 Ω S M curl 1 b V ε × V ε,M • ψ dx dt + T 0 Ω curl 1 b , S M V ε × V ε,M • ψ dx dt .
By what we have just said, we have that

lim ε→0 + T 0 Ω S M curl 1 b V ε × V ε,M • ψ dx dt = T 0 Ω S M curl 1 b V × V M • ψ dx dt .
On the other hand, using the embedding

V ε ε ⊂ L 2 T (W 1,p 3 ) again, we have curl 1 b , S M V ε L 2 T (L p 3 ) 2 -M ,
whereas item (iv) of Lemma 4.7 and Bernstein's inequalities (see Lemma A.1 in the Appendix) imply

V ε,M L 2 T (L p 3 ) 2 3M ( 1 p 2 -1 p 3 ) + ε C M ,
for a constant C M which blows up when M → +∞, but which is uniform in ε > 0. Observe that

3 1 p 2 - 1 p 3 = 3 κ + 3 6κ -1 + 5κ + 3 6κ = 3 κ , with 3 κ < 1 .
Thus, we deduce that lim

M →+∞ lim sup ε→0 + T 0 Ω curl 1 b , S M V ε × V ε,M • ψ dx dt = 0 .
The proof of the lemma is now completed. 2

Handling the pressure term

Before computing the limit with respect to M → +∞ in the previous lemma, let us treat the first term in the right-hand side of (62). As a matter of fact, we need to couple it with a term coming from the pressure function, namely

- 1 ε 2 Ω Π( ε ; b) div ψ dx ,
where Π( ε ; b) has been defined in Lemma 4.6. We will use a fundamental cancellation (appearing after regularisation), which is already present in [START_REF] Masmoudi | Rigorous derivation of the anelastic approximation[END_REF]. For this, we need the following preparatory lemma.

Lemma 5.3 For any T > 0 fixed and any test function ψ as in (56), such that Supp ψ ⊂ [0, T [ ×Ω, the following relation holds true, in the sense of (58):

- 1 ε 2 T 0 Ω Π( ε ; b) div ψ dx dt = - 1 2 T 0 Ω ε ess p (b) φ 2 ε div ψ dx dt + R ε,M .
Proof. Let us start by decomposing the term on the left-hand side of the claimed equality into two integrals, one on the essential set and the other on the residual set:

1 ε 2 T 0 Ω Π( ε ; b) div ψ dx dt = 1 ε 2 T 0 Ω ε ess Π( ε ; b) div ψ dx dt + 1 ε 2 T 0 Ω ε res Π( ε ; b) div ψ dx dt .
Notice that, in order to treat the integral over Ω ε res , the bounds of Lemma 4.6 are not enough. Instead, in order to take advantage of the smallness bounds in [START_REF] Levermore | Global well-posedness for models of shallow water in a basin with a varying bottom[END_REF], we need to introduce a finer decomposition and split that term further into two pieces.

To begin with, let us consider

Π E ( ε ; b) = p E ( ε )-p E (b)-p E (b)( ε -b)
only. On the one hand, thanks the uniform bounds of (35), we have, for some z ε = z ε (t, x) belonging to the interval ] ε (t, x), b(x)[ , the estimate

1 ε 2 Ω ε res,B Π E ( ε ; b) div ψ dx 1 ε Ω ε res,B p E (z ε ) + p E (b) φ ε div ψ dx 1 ε 1 Ω ε res,B φ ε L 2 T (L 6 ) L(Ω ε res ) 5/6 ε -1+5/3 = ε 2/3 ,
which tells us that the contribution coming from this term is a remainder, in the sense of (57). On the other hand, proceeding similarly and employing the assumption p E (z) ≈ z γ-1 , we get

1 ε 2 Ω ε res,U B Π E ( ε ; b) div ψ dx 1 ε Ω ε res,U B γ-1 ε φ ε div ψ dx .
Now, owing to (36), we gather the uniform bound

1 ε 1 Ω ε res,U B γ-1 ε ε ⊂ L r T (L p ) , r := γ γ -1 and p := 3 γ γ -1 , (64)
whereas we can use (40) to control the φ ε term. After checking that

∀ 1 < γ ≤ 2 , γ -1 γ + 1 2 ≤ 1 and 1 q := 3 -γ 3 γ + γ -1 3 γ < 1 ,
we can thus bound

1 ε 2 T 0 Ω ε res,U B Π E ( ε ; b) div ψ dx dt 1 ε 1 Ω ε res,U B γ-1 ε L r T (L p ) φ ε L 2 T (L 3γ/(3-γ) ) L(Ω res ) 1-1/q ε 2(1-1/q) .
A direct computation shows that 1 -1/q > 0. Therefore, for 1 < γ ≤ 2, we have proved that the contribution coming from the integral over the residual set is indeed a remainder, in the sense of (57).

In the case γ > 2, instead, we interpolate between the L ∞ T (L γ/(γ-1) ) bound on 1 Ω ε res,U B γ-1 ε ε , which comes from [START_REF] Levermore | Global well-posedness for models of shallow water in a basin with a varying bottom[END_REF], and the bounds provided by (64). We find that

1 Ω ε res,U B γ-1 ε ε ⊂ L 2 T (L m ) , m := 3γ 2γ -3 ≥ 6 5 .
Hence, owing to (40) again, when γ > 2 we can set θ := γ-2 2(γ-1) and estimate 1) .

1 ε 2 T 0 Ω ε res,U B Π E ( ε ; b) div ψ dx dt 1 ε 1 Ω ε res,U B γ-1 ε L 2 T (L m ) φ ε L 2 T (L 6 ) L(Ω ε res ) 5/6-1/m 1 ε 1 Ω ε res,U B γ-1 ε θ L ∞ T (L γ/(γ-1) ) 1 Ω ε res,U B γ-1 ε 1-θ L r T (L p ) ε 2(5/6-1/m) ε 2/3 ε -2/m ε 2θ(γ-1)/γ ε 1-θ = ε 4/3 ε -(γ-2)/2(γ-
A simple computation shows that

∀ γ > 2 5 , 4 3 - γ -2 2(γ -1) > 0 ,
which finally implies that, also in the case γ > 2, the integral over the residual set is a remainder, in the sense of (57).

In the end, we have proved that, for any value of γ > 1, the contribution coming from the integral of Π E over the residual set is a remainder, in the sense of relation (57).

Next, let us consider the integral involving

Π c ( ε ; b) := p c ( ε ) -p c (b) -p c (b)( ε -b) over the residual set.
For this term, the roles of Ω ε res,B and Ω ε res,U B are inverted. For instance, by Taylor formula at the first order we can write, for suitable z ε (t, x) ∈ ]b, ε (t, x)[ , the following estimate:

1 ε 2 Ω ε res,U B Π c ( ε ; b) div ψ dx 1 ε 2 Ω ε res,U B p c (z ε ) + p c (b) ε div ψ dx 1 ε 2 1 Ω ε res,U B ε L γ T (L 3γ ) L(Ω ε res ) 1-1/(3γ) ε 1-2/(3γ) ,
which obviously converges to 0 when ε → 0 + . In the previous computation, we have used (36) and ( 21) to absorb the negative powers of ε.

In Ω res,B , instead, Taylor formula and the fact that z

ε (t, x) ∈ ] ε (t, x), b[ yield 1 ε 2 T 0 Ω ε res,B Π c ( ε ; b) div ψ dx dt 1 ε 2 T 0 Ω ε res,B p c (z ε ) + p c (b) (b -ε ) div ψ dx dt 1 ε 2 T 0 Ω ε res,B -κ-1 ε dx dt 1 ε 2 1 Ω ε res,B -κ/2 ε 2 L 2 T (L 6 ) -1 ε res L ∞ T (L κ ) L(Ω ε res ) 2/3-1/κ
, which, in view of (37) and ( 21), is of order O(ε 4/3 ).

To sum up, we have shown that also the contribution from Π c over the residual set is a remainder, in the sense of (57).

In light of what we have shown above, in order to complete the proof of the lemma it remains us to deal with the integral on the essential set. For this, we remark that, by Taylor formula, we can write

1 ε 2 T 0 Ω ε ess Π( ε ; b) div ψ dx dt = 1 2 T 0 Ω ε ess p (b) φ 2 ε div ψ dx dt + ε 6 T 0 Ω ε ess p (z ε ) φ 3 ε div ψ dx dt ,
for some z ε = z ε (t, x) belonging to the interval joining b(x) and ε (t, x). At this point, we notice that the second integral on the right-hand side of the previous equality can be bounded, with the help of ( 23) and (35), in the following way:

T 0 Ωess p (z ε ) φ 3 ε div ψ dx dt φ 2 ess L 1 T (L 3 ) φ ε ess L ∞ T (L 2 )
1 .

Putting everything together, we finally get the claimed relation. 2

Of course, in the equality of Lemma 5.3 there is no dependence of the remainder on the approximation parameter M . However, in order for this information to be useful, we need to introduce the regularisation in the first term of the right-hand side. This is easy, thanks to Proposition 3.2 and the fact that, on the essential set, the function φ ε possesses higher integrability in space (keep in mind property (35)).

Indeed, for any M ∈ N, we can write

φ ε 2 ess = φ ε ess φ ε,M + φ ε ess Id -S M φ ε .
Concerning the last term on the right, for any T > 0 fixed, we proceed as follows: we use the uniform bounds φ ε ess ε ⊂ L 2 T (L 6 ) and, for 6/5 ≤ γ < 2, the fact that

Id -S M φ ε L 2 T (L γ ) 2 -M Id -S M ∇φ ε L 2 T (L γ ) 2 -M . (65)
This control follows from Bernstein's inequality (see Lemma A.1 in the Appendix) and the uniform bounds of Proposition 3.2, by arguing in a similar way as done in (59). When 1 < γ < 6/5, instead, we can use an interpolation argument between Lebesgue norms (because the inequality 3γ/(3 -γ) > 6/5 holds for any γ > 6/7) to get, for a suitable θ ∈ ]0, 1[ , the following series of inequalities:

Id -S M φ ε L 2 T (L 6/5 ) Id -S M φ ε θ L 2 T (L γ ) Id -S M φ ε 1-θ L 2 T (L 3γ/(3-γ) ) 2 -θM ∇φ ε θ L 2 T (L γ ) ∇φ ε 1-θ L 2 T (L γ ) 2 -θM .
Finally, when γ ≥ 2, we can simply use Proposition 3.2 and relation (68) of the Appendix to find that

Id -S M φ ε L 2 T (L 2 ) 2 -M . Next, we further decompose φ ε ess φ ε,M = φ ε,M 2 ess + 1 Ω ε ess Id -S M φ ε φ ε,M .
At this point, we notice that all the uniform bounds satisfied by φ ε ε are also satisfied by φ ε,M ε , uniformly with respect to both M ∈ N and ε ∈ ]0, 1]. Thus, the same argument as above yields

1 Ω ε ess Id -S M φ ε φ ε,M L 1 T (L 1 ) 2 -θM ,
where θ = 1 if γ ≥ 6/5, whereas θ ∈ ]0, 1[ is as above in the case 1 < γ < 6/5.

In the end, putting everything together, we have discovered that

- 1 ε 2 T 0 Ω Π( ε ; b) div ψ dx dt = - 1 2 T 0 Ω ε ess p (b) φ 2 ε,M div ψ dx dt + R ε,M . (66)

Coupling the pressure with the convective term

We are now ready to deal with the first term appearing in the right-hand side of (62). We have to pay attention, because here the signs are important. We start by writing, for ψ as in (56),

- Ω φ ε,M ∇ H (b) φ ε,M • ψ dx = - 1 2 Ω H (b) ∇φ 2 ε,M • ψ dx - Ω φ 2 ε,M H (b) ψ • ∇b dx = - 1 2 Ω H (b) ∇φ 2 ε,M • ψ dx + Ω φ 2 ε,M H (b) b div ψ dx ,
where we have also used the fact that div(b ψ) = 0. An integration by parts shows that

- 1 2 Ω H (b) ∇φ 2 ε,M • ψ dx = 1 2 Ω φ 2 ε,M H (b) div ψ dx + 1 2 Ω φ 2 ε,M H (b) ∇b • ψ dx = 1 2 Ω φ 2 ε,M H (b) div ψ dx - 1 2 Ω φ 2 ε,M H (b) b div ψ dx .
At this point, we insert this expression into the previous one; after exploiting the definition H (z) = p (z)/z for all z > 0, we finally gather

- Ω φ ε,M ∇ H (b) φ ε,M • ψ dx = 1 2 Ω φ 2 ε,M div ψ p (b) dx .
Owing to item (i) of Lemma 4.7 and estimate [START_REF] Lions | Mathematical topics in Fluid Mechanics[END_REF], we easily see that the integral over the residual set is small. Therefore, after integrating also in time, in the end we get

- T 0 Ω φ ε,M ∇ H (b) φ ε,M • ψ dx dt = 1 2 T 0 Ω ε ess φ 2 ε,M div ψ p (b) dx dt + R ε,M .
The fundamental point, here, is that the first term on the right-hand side exactly cancels out with the term coming from (66).

Limit of the convective term: conclusion

Putting together Lemma 5.1 and the computations of Subsection 5.2, we finally discover that, for any test function ψ belonging to the kernel of the singular perturbation operator, namely such that (56) holds true, one has

lim ε→0 + T 0 Ω ε u ε ⊗ u ε : ∇ψ dx dt - 1 ε 2 T 0 Ω Π( ε ; b) div ψ dx dt dt = lim M →+∞ T 0 Ω S M curl 1 b V × V M • ψ dx dt .
At this point, remark that, since V is a weak-limit point of the sequence V ε ε , in view of Corollary 4.4 and Lemma 4.7, it enjoys the following property:

V ∈ L ∞ loc R + ; L 2 (Ω) ∩ L 2 loc R + ; W 1,p 1 (Ω) , p 1 := 2κ κ + 1 .
In particular, we also have V ∈ L 2 loc R + ; L p 2 (Ω) , where p 2 := 6κ/(κ + 3). Hence, repeating the computations used in the final part of the proof to Lemma 5.2, we get that lim

M →+∞ T 0 Ω S M curl 1 b V × V M • ψ dx dt = lim M →+∞ T 0 Ω curl 1 b V M × V M • ψ dx dt ,
and performing computations in (61) backwards, we deduce that, for any test function ψ such that div (b ψ) = 0, we have lim

M →+∞ T 0 Ω S M curl 1 b V × V M • ψ dx dt = lim M →+∞ - T 0 Ω 1 b V M ⊗ V M : ∇ψ dx dt ,
where we have also used the fact that div V = 0 (as it follows from taking the limit in the mass equation, recall Proposition 4.5 above). Now, using that V ∈ L 2 T (W 1,p 1 ) and arguing as in (59), it is easy to see that, for almost any t ∈ [0, T ], one has S M V(t) -V(t) L p 1 2 -M ∇V(t) L p 1 , which immediately implies that S M V -→ V strongly in L 2 T (L p 1 ) . In fact, this convergence holds true even in L 2 T (B 1 p 1 ,2 ) (notice that, by (69) below, we have W 1,p 1 → B 1 p 1 ,2 ), owing to Lemma A.4 and the Lebesgue dominated convergence Theorem; however, the previous weaker convergence result is enough for our scopes.

Therefore, S M V converges strongly to V in any intermediate space between L 2 T (L p 1 ) and L 2 T (L p 2 ), thus also in L 2 T (L 2 ) for instance. Thanks to this latter property, we can compute

lim M →+∞ - T 0 Ω 1 b V M ⊗ V M : ∇ψ dx dt dt = - T 0 Ω 1 b V ⊗ V : ∇ψ dx dt dt .

Deriving the asymptotic system: final computations

In Subsections 5.1 and 5.2, we have seen how passing to the limit in the convective term and the pressure term. About the latter, we recall that we have to make use of Lemma 4.6, and more precisely of the relation

1 ε 2 ∇p( ε ) -ε ∇G = 1 ε b ∇ H (b) φ ε + 1 ε 2 ∇Π( ε ; b)
, where the first term on the right disappears whenever tested again a test function satisfying (56), whereas the second term is combined with the convective term to give rise to small remainders, in the sense of relations (57) and (58).

On the other hand, the same computations performed in Proposition 4.5 show how dealing with the continuity equation for the densities ε and with the time derivative term ∂ t ( ε u ε ) ε in the momentum equation. Therefore, in order to complete the proof of Theorem 2.2, we must show convergence of the viscosity term

ν T 0 Ω ε Du ε : ∇ψ dx dt ,
where ψ is as in (56) and is such that Supp ψ ⊂ [0, T [ ×Ω. We start by observing that only the integral over Ω ε ess matters, owing to the uniform bounds

√ ε ∇u ε ε ⊂ L 2 T (L 2 ) and √ ε res ε ⊂ L ∞ T (L 2γ
), with 1/2 + 1/(2γ) < 1.

Next, on Ω ε ess we use the strong convergence ε → b in L ∞ T (L 2 ) ∩ L 2 T (L 6 ) and the weak convergence Du ε DU in L 2 T (L p 1 ). We observe that 1/6 + 1/p 1 ≤ 1. Therefore, we deduce that, for any test function ψ as above, we have A Appendix: elements of Fourier analysis

We recall here the main ideas of Littlewood-Paley theory, which we will exploit in our analysis. The classical construction is usually given in the R d setting: we refer e.g. to Chapter 2 of [START_REF] Bahouri | Fourier Analysis and Nonlinear Partial Differential Equations[END_REF] for details. However, everything can be adapted (see e.g. reference [START_REF] Danchin | Zero Mach number limit for compressible flows with periodic boundary conditions[END_REF]) to cover also the case of a d-dimensional periodic box T d a , where a ∈ R d (this means that the domain is periodic in space with, for any 1 ≤ j ≤ d, period equal to 2πa j with respect to the j-th component).

For simplicity of presentation, we focus here on the case where a j = 1 for all j, and we simply write the spacial domain as T d . We also denote by |T d | = L(T d ) the Lebesgue measure of the box T d . First of all, let us recall that, for a tempered distribution u ∈ S (T d ), we denote by Fu = u k k∈Z d its Fourier series, so that we have

u(x) = 1 |T d | 1/2 k∈Z d u k e ik•x .
Next, we introduce the so called Littlewood-Paley decomposition, based on a non-homogeneous dyadic partition of unity with respect to the Fourier variable. We fix a smooth scalar function ϕ such that 0 ≤ ϕ ≤ 1, ϕ is even and supported in the ring r ∈ R 5/6 ≤ |r| ≤ 12/5 , and such that Finally, we introduce the following low frequency cut-off operators: for any j ∈ Z, we define

∀ r ∈ R \ {0} ,
S j u := u 0 + m≤j-1 ∆ m u . (67)
We explicitly remark that, for any j ∈ Z, the operators ∆ j and S j are linear operators which are bounded on L p for any p ∈ [1, +∞], with norm independent of j and p.

At this point, we present a simplified version of the classical Bernstein inequalities, which turns out to be enough for our scopes. We refer to Chapter 2 of [START_REF] Bahouri | Fourier Analysis and Nonlinear Partial Differential Equations[END_REF] for the statement in its full generality.

Lemma A.1 There exists a constant C > 0, only depending on the space dimension d, on the size of the torus T d a and on the support of the function ϕ fixed above, such that the following properties hold true: for any j ∈ Z, for any α ∈ N d , for any couple (p, q) ∈ [1, +∞]2 such that p ≤ q, and for any smooth enough u ∈ S (T d ), we have It is well known that, for all s ∈ R, the space B s 2,2 coincides with H s , with equivalent norms:

∇ α S j u L q ≤ C |α|+1 2 j|α|+jd 1 p -1 q S j u L p
f 2 H s ∼ | u 0 | 2 + j∈Z 2 2sj ∆ j u 2 L 2 . ( 68
)
When p = 2, non-homogeneous Besov spaces are interpolation spaces between Sobolev spaces W k,p : for all p ∈ ]1, +∞[ , one has the chain of following continuous embeddings:

B 0 p,min(p,2) → L p → B 0 p,max(p,2) . (69)
As an immediate consequence of the Bernstein inequalities, one gets the following Sobolev-type embedding result.

Proposition A.3 Let 1 ≤ p 1 ≤ p 2 ≤ +∞. The, the space B s 1 p 1 ,r 1 is continuously embedded in the space B s 2 p 2 ,r 2 whenever

s 2 < s 1 -d 1 p 1 - 1 p 2 or s 2 = s 1 -d 1 p 1 - 1 p 2 and r 1 ≤ r 2 .
We conclude this appendix by recalling Lemma 2.73 of [START_REF] Bahouri | Fourier Analysis and Nonlinear Partial Differential Equations[END_REF].

  where Ψ is the unique solution to the Neumann problemdiv b ∇Ψ = div v in Ω , Ω Ψ dx = 0 .Remark that P b [v] and Q b [v] := b ∇Ψ are orthogonal in the weighted Hilbert space L 2

Lemma 4 . 6

 46 Let the assumptions of Proposition 4.5 be in force. Then, after definingΠ( ε ; b) := p( ε ) -p(b) -p (b) ε -b ,we have the equality

Lemma 5 . 1

 51 For any T > 0 fixed and any smooth test function ψ as in (56), such that Supp ψ ⊂ [0, T [ ×Ω, we have lim M →+∞ lim sup

ε

  Du ε : ∇ψ dx dt -→ T 0 Ω b DU : ∇ψ dx dt when ε → 0 + .Theorem2.2 is now proven. 2

j∈Z ϕ 2 -j r = 1 .

 21 Let us define |D| := (-∆) 1/2 as the Fourier multiplier 2 of symbol |k|, for k ∈ Z d . The dyadic blocks (∆ j ) j∈Z are then defined by∀ j ∈ Z , ∆ j u := ϕ(2 -j |D|)u = k∈Z d ϕ(2 -j |k|) u k e ik•x .Notice that, for j < 0 negative enough (in general, depending on the box T d a ), one has ∆ j ≡ 0. In addition, one has the following Littlewood-Paley decomposition in S (T d ):∀ u ∈ S (T d ) , u = u 0 + j∈Z ∆ j u in S (T d ) .

andC

  -|α|-1 2 -j|α| ∆ j u L p ≤ ∇ α ∆ j u L p ≤ C |α|+1 2 j|α| ∆ j u L p ,where we have denoted |α| := j α j .By use ofLittlewood-Paley decomposition, we can now define the class of Besov spaces. Definition A.2 Let s ∈ R and 1 ≤ p, r ≤ +∞. The non-homogeneous Besov space B s p,r = B s p,r (T d ) is the set of tempered distributions u ∈ S (T d ) for which with the standard modification in the definition of the norm in the case when r = +∞.

Lemma A. 4

 4 If 1 ≤ r < +∞, for any f ∈ B s p,r one has lim j→+∞ f -S j f B s p,r = 0 .

Here and throughout this paper, we make use of the following notation: given a normed space X and a sequence of functions fε ε all belonging to X, we write fε ε ⊂ X implicitly meaning that the sequence is also bounded in X.

Throughout we agree that f (D) stands for the pseudo-differential operator u → F -1 (f Fu).

The work of F.F. has been partially supported by the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissement d'Avenir" (ANR-11-IDEX-0007), and by the projects BORDS (ANR-16-CE40-0027-01), SingFlows (ANR-18-CE40-0027) and CRISIS (ANR-20-CE40-0020-01), all operated by the French National Research Agency (ANR).

† The work of E.Z. was partially supported by the EPSRC Early Career Fellowship no. EP/V000586/1.