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

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 [30]): 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.
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. 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 [29]. The limiting system takes the name of anelastic approximation. The physical importance of the anelastic approximation is discussed, for example in [22] in the context of various atmospheric flows, and in [1] in the context of astrophysics models.
1.1. 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: 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 (1), 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 .
The pressure function p is assumed of the following form: (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; they will arise naturally in the computations of Sects. 3, 4 and 5.
Before going on, let us briefly comment on the cold pressure component p c ( ) = − 1 vacuum regions ≈ 0. We mention that, for similar reasons, a cold pressure term already appeared in works [4] (in the context of heat conducting fluids), [28] (for mixture of fluids with chemical reactions) and [19] (for some lubrication models in one space dimension), and was recently employed in the study of [6] on the compressible Navier-Stokes equations.
From a more physical perspective (see for instance the discussion in [4]), the term p c ( ) is associated with the zero Kelvin isothermal curve for heat conducting fluids; also, singular pressures naturally appears when e.g. Van der Waals type laws are considered. As a matter of fact, for densities and temperatures close to 0 the properties of the medium drastically change, damaging the validity of the equations of motion: the presence of p c ( ) may be seen as a way of preserving stability of the model. Observe that different terms could be added to the system instead of p c ( ), in order to improve the stability of the model in vacuum regions: for instance, one could use some drag forces (like in [3] by Bresch and Desjardins for a 2-D shallow water model), or one could impose some additional integrability conditions on the initial velocity, as done in [27] by Mellet and Vasseur. We refer to papers [6] and [32] for a complete discussion on this subject, as well as for recent developments. However, we have to remark that, because of technical reasons, those approaches seem not to work in our context: while it is not clear to us that the Mellet-Vasseur estimates can be obtained uniformly in the small parameter ε ∈ ]0, 1], the turbulent drag term from [3] seems to create difficulties in passing to the limit (see more details in Sect. 1.3 below). For this reason, we resort here to the cold component p c ( ) in the pressure law.
When ε → 0 + in equations (1), 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 Smoothness of G(x) and strict monotonicity of p imply that there exists a smooth function b ∈ C 3 ( ) satisfying (4). Monotonicity of p implies convexity of the pressure potential H (defined in (9) below) which provides existence of constants b * , b * ∈ R such that ∀ x ∈ , 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 [7]. On the other hand, if G(x) = −gx 3 is the gravitational potential, it is easy to see that b = b(x 3 ) verifies (5).
Since ≈ b for ε → 0 + , linking the pressure and external force terms according to (4) and assuming that we can identify the limits of the non-linear terms appearing in (1), we formally check that the limiting system is an anelastic approximation with variable viscosity coefficient, namely We refer to system (6) 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 [23]. 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. 1.2. 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 [26], 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 [15]. Finally, we refer to the book by Feireisl and Novotný [16] 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 [10] by Ebin and [20,21] 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 [14], 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 [14] 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 [26], 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 [9], 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 [13], 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 [7] 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 [3] by Bresch and Desjardins and in [18] 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 [11,12] for a generalisation of these results to the 3-D setting and to the case of general ill-prepared initial data.
1.3. The content of the paper. In the context depicted above, our work can be seen as a generalisation of the result from [7], 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 [28,33]. The case without this assumption has been completed much more recently in [32] 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 [3].
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 (6). It is well known that passing to the limit in the weak formulation of equations (1), 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 [7]. Our approach is inspired by the previous works [25] and [26] 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 (1). 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 Sect. 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 Sect. 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 Sect. 5 we rigorously perform the convergence in the weak formulation of equations (1), and conclude the derivation of the anelastic approximation (6). 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. Initial data. Problem (1)-(3) is supplemented by general ill-prepared initial data. Namely, for any small parameter ε ∈ ]0, 1] fixed, we pick initial data 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 assume 1 that (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 Energy functionals. Next, we need to introduce various energy functionals. The internal energy function (sometimes called pressure potential) is defined by the ODE which implies in particular that Notice that H is defined up to the sum of an affine function. Here, we fix the classical choice to be the classical energy and the BD entropy functions. We also set E , , 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.
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: (3) The equations of system (1) are satisfied in the sense of distributions: more precisely, we have for any test function ξ ∈ D [0, T [ × , and for any test function ψ ∈ D [0, T [ × ; R 3 ; (4) For almost every t ∈ [0, T [ , the following energy inequalities hold true: 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 [33] and [28], 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 [28][29][30][31][32][33] in order to include the force.
Main result. Before stating the main result of this paper, we need some additional tools and notation. Following [15][16][17][18][19][20][21][22][23][24][25] (see also [24]), 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 where is the unique solution to the Neumann problem , which is defined as the space of functions f : −→ R 3 which are L 2summable with respect to the measure 1 b dx. Similarly to the case of the classical Leray-Helmholtz projector P = P 1 and its L 2orthognal 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 < +∞.
We can now state the main result of this paper, which is contained in the following theorem.
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 .
for any T > 0 and any test function

Remark 2.3. Note that the initial condition equals
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 We also point out that, in the computations below, the Lebesgue exponent p ∈ [1, +∞] is allowed to vary from an inequality to another. Each time we use it, we put in evidence the right range of values of p for which the corresponding inequality holds true. Thus, for instance, the exponent p appearing below is completely uncorrelated to the one from Theorem 2.2.

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 where the function E , u b has been defined in (10) and we recall that we have set 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 (17), we now derive first uniform bounds for the family ε , u ε ε . In fact, owing to our assumptions on the initial data, and in particular to (15), the right-hand side of (17) is uniformly bounded: specifically, we have Then, it is easy to deduce the following uniform bounds: Let us now focus on the density functions. To begin with, following the approach of [16], 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 where the constants b * and b * have been defined in (5) Then, given a function h, we can write Here above, 1 A denotes the characteristic function of a set A ⊂ . For later use, it is convenient to divide the residual set ε as the regions where ε , respectively, stays bounded and may become unbounded. After this preparation, let us come back to (17) and derive uniform bounds for the family ε ε . In the essential set, we can perform Taylor's expansion of the function H ; we thus get which implies that 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. [17] for details) the following bounds on the residual set: The previous estimate immediately implies that where we have denoted by L(A) the Lebesgue measure of a set A ⊂ . At this point, let us define the quantity From the uniform bound (20), we may deduce that On the other hand, using (21) we can compute that for any p ≤ γ , we have Thus, we finally infer that 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. [3,5]). 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 [33] and [28] 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.
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 [7], 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 [7] it follows that for sufficiently smooth solutions of the continuity equation in (1), the following equality holds true Step 2. In this step, one multiplies the momentum equation by ν∇ ln b , we get Let us now rewrite each term from the above expression. First, note that the transport term gives For the diffusion term, after noticing that ∇ 2 is always symmetric, we have Finally, for the pressure and force term, using (4), we obtain Step 3. Now, we sum up the equalities from the previous steps along with the energy estimate. In our case, after setting and so, we finally get 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 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 ( ) for some z between and b, and the fact that is bounded in ess to write 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 (23). Case 2: integral over ε res . In the residual set, it is convenient to split J 3 into two pieces, as follows: The bound of the first term in the last equality is easy. We start by writing Then, we observe that uniformly in ε > 0, owing to (21). Indeed, one has p ( ) Thus, using (21) quantitatively and arguing as in Case 1, we deduce that 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 (28). We start by writing First note that Thus, on the one hand we get On the other hand, we also have Now, in view of (21), for 2 − γ ≥ 0 we have that If, on the other hand 2 − γ < 0, then we have In the latter case, in order to have q ≥ 2, we need the condition In any case, and under the assumption that the previous requirement is satisfied in the case γ > 2, we can thus bound In the end, using (21), (22) and the definition of m, after an application of Young's inequality we arrive at Finally, putting all those inequalities together, we see that Hence, for δ > 0 small enough, we can absorbe the last term of the previous inequality into the left-hand side of (27). The proof of the lemma is thus completed.
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.

3.3.
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 (17). Those additional bounds will be fundamental in order to prove convergence in our setting.
Next, combining (25) with (17), we gather that with uniform inclusion (of course, uniform with respect to ε ∈ ]0, 1]). At this point, an easy computation shows that (21). Putting this information together with (30), by Sobolev embeddings we finally deduce that 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 The term inside the parentheses can be written as Thus, we get the following uniform embeddings: 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 Therefore, the uniform bound (33) implies that Let us derive a couple of properties from (34). We notice that we can write On the one hand, from the previous relation we immediately deduce that 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: This inequality, together with (34) above, implies that On the other hand, on the set ε res,U B we have that ε ≥ 2 b * . Hence, on that set we can write Using (34) again, we gather that 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 We conclude this part by showing a third consequence of the BD entropy estimate (25), 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 Since on the set ε ess ∪ ε res,B we have 0 ≤ ε ≤ 2b * , we easily get which implies that, for all T > 0 fixed, one has 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 Case 2. For 1 − γ /2 ≥ 0 instead, we can write Since, for any T > 0 fixed, ε res ε is uniformly bounded in L ∞ T (L γ ), we have that In turn, combining this information with (38) implies that 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 and, when γ ≤ 2, we bound φ ε ε uniformly in L ∞ T (L γ ), thanks to (23) and (24). When γ > 2, instead, we use (23) again, together with the fact that, by definition of φ ε , we have where we have used also the bounds provided by (21). We thus conclude that, when γ > 2, we have φ ε ε ⊂ L ∞ T (L 2 ) for any T > 0 fixed. In the end, the proposition is proved.
We also notice that, as a consequence of Proposition 3.2 and Sobolev embedding, we get that

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 Sect. 4.1 we establish first convergence properties for the sequence of solutions ε , u ε ε towards some targe profile , u . Then, in Sect. 4.2 we derive some constraints that the limit point , u has to satisfy. Finally, in Sect. 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 Sects. 3.1, 3.2 and 3.3, we can derive first convergence properties for the family of solutions ε , u ε ε of our primitive system (1). 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 (6).
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 (41) Since is of finite measure, we can interpolate those convergence properties with (31) to deduce that The uniform bounds of Sect. 3 allow us to find more quantitative convergence properties. Indeed, from (23) 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 where the symbol * stands for the weak- * convergence in the respective functional space. On the other hand, owing to (24), we know that 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

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 (18), (29). Here, we decide to work with the momentum However, the first step is to get some uniform bounds on the velocity fields u ε . This is the goal of the next proposition. ∇u ε ε ⊂ L 2 loc R + ; L p 1 ( ) , p 1 := 2κ κ + 1 .
In particular, we also have that Remark 4.2. Note that, due to our assumption κ > 3, one has p 1 > 3/2 and p 2 > 3 in the above proposition.
Finally, the last uniform bound for the family of u ε 's follows from the previous property and Sobolev embeddings.

Proposition 4.3. Let ε , u ε ε be a family of global in time finite energy weak solutions to system
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 with the uniform embedding properties (iii) After setting V ε := b u ε , we can also write 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 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 ∇ 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 At this point, we use that Hence, on the one hand, by using (35), we see that on the other hand, since in view of (36) we gather that also 1 ε

The proof of item (iii) is similar. This time, we write
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 , for all T > 0 fixed. For this, owing to (18), the last part of Proposition 4.1 and Remark 4.2, it is enough to show that As a matter of fact, since (20) and (21) we immediately get that 1 ε ess + 1 ε is uniformly bounded in L ∞ T (L 2 ). Finally, as done above, we have which, in view of (21) again, implies that In the end, we have shown (44), which in turn implies the sought bound for the vector fields W ε .
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 ε * To conclude this part, we define the target velocity field U as 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. Then U has to satisfy the anelastic constraint

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.
we have the equality In addition, one has that 1 Proof. The claimed identity follows from a simple algebraic computation: 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 The argument is similar to the proof of Lemma 4.1 in [13]; 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 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 Finally, we notice that for which one can deduce the sought bounds by using the controls in (21).
In the end, the lemma is completely proved.
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 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 . Thus, passing to the limit in the previous equality yields the constraint 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 Sect. 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 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 Using the definition of φ ε , we can recast the continuity equation as 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 (8). 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.
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 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: where we have defined By the uniform bounds (18), 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 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 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 [2]) imply that 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 [26], thanks to the uniform bounds provided by the BD entropy estimates.
In addition, as an immediate consequence of (49), we get that 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 Sect. 5 below, it is important to introduce a fine decomposition of those terms. (ii) For any M ∈ N and any ε ∈ ]0, 1], we can write where we have the uniform estimates 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 with the uniform bounds with the uniform estimates Finally, let us prove item (ii). We start by considering the case 1 < γ < 2. In this case, we decompose Thus, if we set ϕ ε,M := S M φ ε ess , estimates (23) 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 (24). 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.

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 (6).
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 Sect. 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 Also, it is useful to introduce the following notation: we denote by R ε,M any remainder term, that is any term such that 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 Typically, we will have 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 [26] 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. 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 where all terms are well-defined, and the last term on the right can be bounded as follows: . At this point, taking advantage of Bernstein's inequalities of Lemma A.1 in the Appendix, we can bound 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 As a result, the previous computations show that, in the sense of (58), one has 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 Indeed, first note that, due to (40), we can write, for γ < 2, the following estimate: 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 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 This last equality finally ends the proof of the lemma.

Compensated compactness.
Owing to Lemma 5.1, it is enough to compute the limit of the approximate convective term 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
where the symbol × denotes the usual 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 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 Owing to (52) and item (i) of Lemma 4.7, it is clear that the second term in the righthand 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 Thanks to (51) and item (ii) of Lemma 4.7, we can estimate 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 As a result of the previous computations, we infer that

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
Proof. We start the proof by recalling that, owing to Proposition 4.3, we have V ε ε ⊂ L 2 Observe that, by dual Sobolev embeddings (see e.g. Theorem 0.5 of [16]), 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 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 [31]), we gather that for any compact set K ⊂ . This implies that, for any M ∈ N and any compact K ⊂ fixed, the sequence Next, we write By what we have just said, we have that On the other hand, using the embedding V ε ε ⊂ L 2 T (W 1, p 3 ) again, we have whereas item (iv) of Lemma 4.7 and Bernstein's inequalities (see Lemma A.1 in the Appendix) imply for a constant C M which blows up when M → +∞, but which is uniform in ε > 0. Observe that Thus, we deduce that The proof of the lemma is now completed.

5.2.3.
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 where ( ε ; b) has been defined in Lemma 4.6. We will use a fundamental cancellation (appearing after regularisation), which is already present in [26]. For this, we need the following preparatory lemma.
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 (21), and the bounds provided by (64). We find that 1 ε Hence, owing to (40) again, when γ > 2 we can set θ := γ −2 2(γ −1) and estimate 1 A simple computation shows that 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: 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 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: Finally, when γ ≥ 2, we can simply use Proposition 3.2 and relation (68) of the Appendix to find that Id − S M φ ε L 2 Owing to item (i) of Lemma 4.7 and estimate (22), we easily see that the integral over the residual set is small. Therefore, after integrating also in time, in the end we get The fundamental point, here, is that the first term on the right-hand side exactly cancels out with the term coming from (66).

5.2.5.
Limit of the convective term: conclusion Putting together Lemma 5.1 and the computations of Sect. 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 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: 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 and performing computations in (61) backwards, we deduce that, for any test function ψ such that div (b ψ) = 0, we have 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 which immediately implies that 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 5.3. Deriving the asymptotic system: final computations. In Sects. 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 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 Theorem2.2 is now proven.
Data sharing Data sharing not applicable to this article as no datasets were generated or analysed during the current study.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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 [2] for details. However, everything can be adapted (see e.g. reference [8]) 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 Next, we introduce the so called Littlewood-Paley decomposition, based on a nonhomogeneous 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 ∀ r ∈ R\{0} , j∈Z ϕ 2 − j r = 1 .
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 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 ): Finally, we introduce the following low frequency cut-off operators: for any j ∈ Z, we define S j u := u 0 + m≤ j−1 m u .
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 [2] 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 where we have denoted |α| := j α j . By use of Littlewood-Paley decomposition, we can now define the class of Besov spaces. It is well known that, for all s ∈ R, the space B s 2,2 coincides with H s , with equivalent norms: 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) .
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 We conclude this appendix by recalling Lemma 2.73 of [2].