Edoardo Bocchi

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

Christophe Prange
email: christophe.prange@cyu.fr

ANISOTROPY AND STRATIFICATION EFFECTS IN THE DYNAMICS OF FAST ROTATING COMPRESSIBLE FLUIDS

Keywords: November 4, 2020. 2010 Mathematics Subject Classification. 35Q86, 35Q30, 76D50, 76N10, 35B40, 76M45 Compressible Navier-Stokes equations, rotating fluids, anisotropy, stratification, Ekman boundary layers, relative entropy, maximal regularity

The primary goal of this paper is to develop robust methods to handle two ubiquitous features appearing in the modeling of geophysical flows: (i) the anisotropy of the viscous stress tensor, (ii) stratification effects. We focus on the barotropic Navier-Stokes equations with Coriolis and gravitational forces. Two results are the main contributions of the paper. Firstly, we establish a local well-posedness result for finite-energy solutions, via a maximal regularity approach. This method allows us to circumvent the use of the effective viscous flux, which plays a key role in the weak solutions theories of Lions-Feireisl and Hoff, but seems to be restricted to isotropic viscous stress tensors. Moreover, our approach is sturdy enough to take into account non constant reference density states; this is crucial when dealing with stratification effects. Secondly, we study the structure of the solutions to the previous model in the regime when the Rossby, Mach and Froude numbers are of the same order of magnitude. We prove an error estimate on the relative entropy between actual solutions and their approximation by a large-scale quasi-geostrophic flow supplemented with Ekman boundary layers. Our analysis holds for a large class of barotropic pressure laws.

Introduction

This paper is devoted to the study of a class of barotropic Navier-Stokes systems. Our focus is on developing robust methods to handle two ubiquitous effects appearing in the modeling of geophysical fluids: (i) the strong anisotropy of the viscous stress tensor, (ii) stratification effects. As for the first aspect, we consider systems with anisotropic viscosity tensors (1.1) ∆ µ,ε := µ∆ h + ε∂ 2 3 = µ(∂ 2 1 + ∂ 2 2 ) + ε∂ 2 3 , where the parameters µ and ε are dimensionless numbers such that ε µ. Concerning the second point, we consider highly rotating fluids with a strong Coriolis force 1 ε e 3 ×ρu and a strong gravitational potential 1 ε 2 ρ∇G = -1 ε 2 ρe 3 , where u = u(x, t) ∈ R 3 denotes the velocity field of the fluid, ρ = ρ(x, t) ∈ R represents the density of the fluid and e 3 = (0, 0, 1) T is the unit vertical vector. We handle barotropic fluids, so the pressure of the fluid is assumed to be a function of the density only, i.e. P = P (ρ); see (1.10) for precise assumptions on the pressure law. We consider the simplest possible geometrical set-ups and scalings enabling to study non trivial phenomena, such as boundary layers and vertical stratification.

Our goal in this paper is twofold. The first part of the paper is concerned with the existence and uniqueness of strong solutions for the barotropic Navier-Stokes equations with potential force ∇G:

(1.2)

∂ t ρ + ∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇P (ρ) = ∆ µ,ε u + λ∇(∇ • u) + ρ∇G .
The term ∇G is responsible for stratification effects, thus the equilibrium density ρ of the previous system becomes non-constant. As a consequence, there are two main difficulties to handle in the well-posedness study for system (1.2): on the one hand, the strong anisotropy of the viscous tensor, see (1.1); on the other hand, the fact that ρ introduces variable coefficients in the equations (see more details below).

Our main result in this direction is Theorem 1 on page 8. We set system (1.2) in the simple space-time domain R 3 × (0, T ), with T > 0; domains with boundaries, such as (R 2 × (0, 1)) × (0, T ), may be handled by the same method, at the price of more technical difficulties. Here we work with fixed values of the parameters (µ, λ, ε), and do not keep track of how the estimates depend on them: the existence of solutions to (1.2) is a challenge in itself. Our well-posedness result still holds if one incorporates a Coriolis force ρe 3 × u in the left hand side of (1.2).

The second part of the paper is devoted to the asymptotic analysis of the barotropic Navier-Stokes equations in presence of fast rotation and gravitational stratification. We consider the following system:

(1.3)    ∂ t ρ + ∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ρ ε e 3 × u + ∇P (ρ) ε 2 = ∆ µ,ε u + λ∇(∇ • u) + ρ ε 2
∇G. set on the strip (R 2 × (0, 1)) × (0, T ), with the no-slip boundary conditions (1.4) u = 0 at x 3 = 0, 1 .

Here G is the gravitational potential, i.e. G(x 3 ) = -x 3 . We describe the structure of weak solutions in the limit ε → 0 for well-prepared data, analyze the Ekman boundary layers and their effect on the limit quasi-geostrophic flow, and prove quantitative bounds based on relative entropy estimates. Our main result in this direction is Theorem 2 on page 34. Our results hold for a large class of monotone pressure laws.

Here our focus is on the asymptotic behavior for a family of weak solutions, under the assumption that such global-in-time weak solutions exist. The two parts are connected. Indeed, the quantitative stability estimates obtained in the second part lay the ground for a large-time well-posedness result for system (1.3) in the limit when ε → 0. Inspired by previous results for incompressible flows [START_REF] Gallagher | The Tridimensional Navier-Stokes Equations with Almost Bidimensional Data: Stability; Uniqueness; and Life Span[END_REF][START_REF] Charve | Global Well-Posedness and Asymptotics for a Geophysical Fluid System[END_REF][START_REF] Ponce | Global stability of large solutions to the 3D Navier-Stokes equations[END_REF][START_REF] Scrobogna | Highly rotating fluids with vertical stratification for periodic data and vanishing vertical viscosity[END_REF], we believe that strong solutions may be constructed for large data close to the two-dimensional limit quasi-geostrophic flow, by using a variation of the well-posedness result of the first part. Though, dealing with the compressible system (1.3) requires much more work, which is left for a subsequent paper.

1.1. Modelling of geophysical flows. We consider mid-latitudes and high-latitudes motions of fast rotating compressible fluids, typically the Earth's atmosphere or oceans. System (1.3) is a particular case of the general non-dimensional system (1.5)

           ∂ t ρ + ∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) - 1 Re h ∆ h u - 1 Re 3 ∂ 2 3 u -
As usual, the Reynolds number measures the ratio of inertial forces to viscous forces. System (1.5) modeling large-scales geophysical flows has, in particular, a horizontal Reynolds number Re h = U L ν h and a vertical Reynolds number Re 3 = U L ν 3 that may be of different orders of magnitude; see below the comments about the anisotropy in the viscous stress tensor. Here U represents the typical speed of the flow, L the typical length, ν h is the horizontal viscosity and ν 3 the vertical one. The Mach number is defined as Ma = U c , where the constant c is the propagation speed of acoustic waves. For strong jet streams near the tropopause U = 50 m • s -2 , which corresponds to Ma = 0.15, see [START_REF] Klein | Multiple spatial scales in engineering and atmospheric low Mach number flows[END_REF]. The Rossby number, defined as Ro = U 2ΩL , measures the effect of the Earth's rotation; here Ω = 7.3 • 10 -5 s -1 is the module of the Earth's angular velocity. In the case of the Gulf Stream, the length L = 100 km and U = 1 m • s -1 are smaller than the typical oceanic scales, the Rossby number is about Ro = 0.07. Notice that here we consider the f-plane approximation of the Coriolis force.

Stratification. The gravitational force deriving from the geopotential G = -gx 3 tends to lower regions of fluid with higher density and raise regions of fluid with lower density. In the equilibrium configuration, the density profile decreases with respect to the vertical direction. The Froude number is then defined as Fr = U √ gH , where g = 9.81 m • s -2 is the acceleration of gravity. It measures the ratio of inertial forces of a fluid element to its weight. The centrifugal force also derives from a potential. It is often neglected in models for the atmosphere [START_REF] Cushman-Roisin | Introduction to Geophysical Fluid Dynamics[END_REF][START_REF] Feireisl | Multiple scales and singular limits for compressible rotating fluids with general initial data[END_REF], but in certain regimes it can have a dramatic effect. The mathematical analysis of the centrifugal force poses different challenges that we do not dwell upon in the present work.

Anisotropy. In general the horizontal and vertical viscosities are not equal, in particular when dealing with large-scale motions of geophysical flows. For instance, in the ocean the horizontal turbulent viscosity ν h ranges from 10 3 to 10 8 cm 2 • s -1 , while the vertical viscosity ν 3 is much smaller and ranges from 1 to 10 3 cm 2 s -1 . A justification of this fact can be seen in: (i) the anisotropy between the horizontal and the vertical scales of the flows, and (ii) the stabilizing effect of the Coriolis force, which makes the large-scale motion be almost two-dimensional (see also Section 3). Hence, a frequent (and crude) modeling assumption is to suppose that the diffusion in the vertical direction is much weaker than the horizontal viscosity, which is enhanced by turbulent phenomena. For further insights about the physics of anisotropic diffusion, we refer to [START_REF] Temam | Some mathematical problems in geophysical fluid dynamics[END_REF], in particular to equation (2.122), to Chapter 4 of [START_REF] Pedlosky | Geophysical fluid dynamics[END_REF], to [START_REF] Bresch | Rotating fluids in a cylinder[END_REF] and [START_REF] Chemin | Ekman boundary layers in rotating fluids[END_REF].

Scaling. In (1.3), the dimensionless number ε denotes the Rossby number, which measures the strength of the rotation. We consider the scaling where the Rossby, the Mach and the Froude numbers are of the same order of magnitude, i.e. Ro = Ma = Fr = ε. This is the richest scaling, since the effects due to the rotation are in balance with the compressible and gravitational effects. Notice however that other scalings are considered in the physical literature [START_REF] Klein | Scale-Dependent Models for Atmospheric Flows[END_REF], for application to meteorology, as well as in the mathematical literature, see e.g. [START_REF] Feireisl | Multi-scale analysis of compressible viscous and rotating fluids[END_REF], [START_REF] Feireisl | Multiple scales and singular limits for compressible rotating fluids with general initial data[END_REF].

Ekman layers. Ekman boundary layers are regions near horizontal boundaries with no-slip boundary condition where viscous effects balance the Coriolis force. The thickness of these boundary layers (see Part I of [START_REF] Chemin | Mathematical geophysics. An introduction to rotating fluids and the Navier-Stokes equations[END_REF]) is

δ E = ν 3 Ω 1 2 ,
which does not depend on the velocity. Remark that the faster the rotation, the smaller is the layer affected by viscosity. See also Chapter 8 of [START_REF] Cushman-Roisin | Introduction to Geophysical Fluid Dynamics[END_REF].

Pressure laws. We consider barotropic flows, for which the pressure is a function of the density only. A typical example is that of Boyle's laws P (ρ) = aρ γ , with a > 0 and γ ≥ 1. For the precise definition of the pressure law, we refer to (1.10).

1.2. Mathematical challenges related to anisotropy and stratification. We outline here some aspects of the study of compressible viscous fluids. We focus on two points in particular: (i) well-posedness and the difficulties related to the anisotropy in the viscosity, (ii) asymptotic analysis in the presence of stratification.

Well-posedness. For fluids modelled by the incompressible Navier-Stokes equations, the fact that the viscous stress tensor is isotropic or anisotropic does not affect the well-posedness theory. One can prove the existence of weak (Leray-Hopf), mild or strong solutions regardless of the structure of the viscous stress tensor.

For compressible fluids modelled by the compressible Navier-Stokes system though, anisotropy represents a major hurdle. Interestingly, the obstacle has similar roots for several well-posedness theories of weak solutions.

In the isotropic case, µ = ε in (1.2), one has a pointwise relation between the pressure and the divergence of u. Indeed, applying the divergence to the momentum equation and using the algebraic relation ∇

• (∇∇ • u) = ∆∇ • u = ∇ • ∆u, one gets (1.6) ∆ (µ + λ)∇ • u -P (ρ) + P (1) = ∇ • ρ∂ t u + ρu • ∇u .
This suggests to define the following quantity

F := (µ + λ)∇ • u -P (ρ) + P (1),
which is dubbed the effective viscous flux. According to (1.6), one obtains the relation

(1.7) ∆F = ∇ • ρ∂ t u + ρu • ∇u .
The property (1.7) is key to the existence of global finite-energy weak solutions in the Lions [START_REF] Lions | Compressible models[END_REF] and Feireisl [START_REF] Feireisl | Singular limits in thermodynamics of viscous fluids[END_REF][START_REF] Novotny | Introduction to the mathematical theory of compressible flow[END_REF] theory on the one hand, and to the existence of global weak solutions with bounded density in the Hoff theory [START_REF] Hoff | Global solutions of the Navier-Stokes equations for multidimensional compressible flow with discontinuous initial data[END_REF] (see also [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF]).

In the anisotropic case µ = ε, the above analysis breaks down, and the relation between ∇ • u and P (ρ) becomes nonlocal. Indeed, (1.7) becomes

(∆ µ,ε + λ∆)∇ • u -∆P = ∇ • ρ∂ t u + ρu • ∇u ,
so that the definition of a modified effective flux should read

F ani := ∆ -1 (∆ µ,ε + λ∆)∇ • u -P (ρ) + P (1).
The nonlocal operator ∆ -1 (∆ µ,ε + λ∆) changes the picture dramatically. In view of the existence of Hoff-type solutions with bounded density [START_REF] Hoff | Global solutions of the Navier-Stokes equations for multidimensional compressible flow with discontinuous initial data[END_REF][START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF], one of the major flaws of F ani is the lack of boundedness of ∆ -1 (∆ µ,ε + λ∆) on L ∞ . Similarly, for the existence of global finite-energy weak solutions, the nonlocality is a major obstacle to getting the compactness of a sequence of approximate solutions; see [START_REF] Bresch | Rotating fluids in a cylinder[END_REF][START_REF] Bresch | Global existence of weak solutions for compressible Navier-Stokes equations: thermodynamically unstable pressure and anisotropic viscous stress tensor[END_REF]. In the breakthrough work of Bresch and Jabin [START_REF] Bresch | Global existence of weak solutions for compressible Navier-Stokes equations: thermodynamically unstable pressure and anisotropic viscous stress tensor[END_REF], a totally new compactness criterion was proved that enables to prove the existence of global weak solutions to the compressible system (1.2) in the anisotropic case. There remains a restriction that |µ -ε| < λ -µ 3 which is compatible with the modeling of large-scale geophysical flows. A more important limitation of the result in [START_REF] Bresch | Global existence of weak solutions for compressible Navier-Stokes equations: thermodynamically unstable pressure and anisotropic viscous stress tensor[END_REF] is on the pressure law:

γ ≥ 2 + √ 10
2 , where γ is defined in (1.10). On that subject Bresch and Burtea proved recently in [START_REF] Bresch | Global existence of weak solutions for the anisotropic compressible Stokes system[END_REF] the global existence of weak solutions for the quasi-stationary compressible Stokes equations and to the stationary compressible Navier-Stokes system [START_REF] Bresch | Weak Solutions for the Stationary Anisotropic and Nonlocal Compressible Navier-Stokes System[END_REF] with an anisotropic viscous tensor. Their approach is based on the control of defect measures.

The anisotropy prompts us to consider the framework of strong solutions, in particular those with minimal regularity assumptions as in [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF]. In that perspective, a further challenge for the well-posedness theory is the presence of non constant reference density states, due to the gravitational term. Indeed, the fact that the densities are perturbations of a constant state plays a major role in the analysis of [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF]. This question seems to remain broadly unexplored in general. In our work, we are able to incorporate non constant reference density states in the approach of [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF].

Asymptotic analysis. We review here some of the literature concerned with the quantitative analysis of viscous barotropic fluids in high rotation (for an overview of this topic for incompressible fluids, we refer to [START_REF] Chemin | Mathematical geophysics. An introduction to rotating fluids and the Navier-Stokes equations[END_REF]). For systems of the type (1.3), results on the combined low Mach and low Rossby limits were obtained in several directions: well-prepared data (data close to the kernel of the penalization operator, see (3.57) for what it means in our context), ill-prepared data, slip or no-slip boundary data, with or without stratification (centrifugal force or gravitational force), different scalings regimes. We do not attempt to be exhaustive here, but select some works which are relevant to our study here.

In the well-prepared case, paper [START_REF] Bresch | Rotating fluids in a cylinder[END_REF] studies the limit for the same scaling as in (1.3), with no-slip boundary conditions and G = 0. They study the Ekman layers and prove stability estimates in the limit ε → 0 in the case of pressure laws with γ = 2. As far as we know, this seems to be the only work concerned with the study of Ekman layers for compressible fluids. Extending it to more general pressure laws, and taking into account gravitational stratification effect is an obvious motivation for our present work. The analysis with a gravitational potential was considered in [START_REF] Feireisl | Rotating compressible fluids under strong stratification[END_REF], for fluids slipping on the boundary, using the relative entropy method.

In the ill-prepared case, the fact that the initial data is away from the kernel of the penalization operator is responsible for the propagation of high frequency acoustic-Poincaré waves. Different scaling regimes of Ma and Ro were analyzed in the works [START_REF] Feireisl | A singular limit for compressible rotating fluids[END_REF][START_REF] Feireisl | Multi-scale analysis of compressible viscous and rotating fluids[END_REF][START_REF] Feireisl | Multiple scales and singular limits for compressible rotating fluids with general initial data[END_REF]. All these works are concerned with fluids satisfying the slip boundary condition, hence no boundary layers are needed in the asymptotic expansions. Let us point out that the work [START_REF] Feireisl | Multi-scale analysis of compressible viscous and rotating fluids[END_REF] manages to handle the centrifugal force.

1.3. Novelty of our results. We comment here on the main theorems of the paper.

Well-posedness of system (1.2): Theorem 1 on page 8. We prove the short-time existence of finite-energy solutions to system (1.2). We introduce a simple and sturdy method based on a priori bounds obtained via maximal regularity estimates, following the approach of [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF][START_REF] Shibata | Global existence of classical solutions and optimal decay rate for compressible flows via the theory of semigroups[END_REF]. However, our techniques are robust enough to deal with both effects mentioned above: (i) the anisotropy in the viscous stress tensor, (ii) the presence of a non-constant reference density state. Thus, our result represents a generalization of [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF][START_REF] Shibata | Global existence of classical solutions and optimal decay rate for compressible flows via the theory of semigroups[END_REF] in both directions. However, as already explained in Subsection 1.2, the anisotropy makes it impossible to use Hoff's effective viscous flux as in [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF], so we need to look for higher regularity estimates for the density, in order to get compactness for passing to the limit in the pressure term. For this reason, we need to require some regularity on the initial data: roughly, the initial density ρ in is sufficiently well-localised around the reference profile ρ, namely ρ in -ρ ∈ H 2 , while the initial velocity u in is taken in H 3/2+ . So, the solutions that we construct are strong solutions with finite energy; roughly, they are halfway between the Hoff solutions [START_REF] Hoff | Global solutions of the Navier-Stokes equations for multidimensional compressible flow with discontinuous initial data[END_REF][START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF] with bounded density and the strong solutions of Matsumura-Nishida [START_REF] Matsumura | The initial value problem for the equations of motion of viscous and heat-conductinve gases[END_REF] (which require ρ in -1 ∈ H 3 ).

We emphasize that the initial density ρ in can be taken close to an arbitrary reference density profile ρ = ρ(x 3 ) ∈ W 3,∞ (R). Thus, setting G = G(x 3 ) = H ρ(x 3 ) , with H defined in (1.11) below, we infer that ρ is a static state of system (1.2), namely ρ satisfies the logistic equation (1.8) ∇P (ρ) = ρ∇G.

However, the smallness condition rests only on the quantity ρ in -ρ L ∞ , and not on higher-order derivatives. Finally, we point out that the assumption ρ = ρ(x 3 ) is made only for modelling purposes (we have in mind the case when G is the gravitational potential), but our method works also in the more general situation ρ = ρ(x), x ∈ R 3 . So, our theorem opens the way to achieving the well-posedness of systems with stratification effects, such as (1.3), in the strip R 2 × (0, 1) with the physical gravitational potential G = -x 3 . Asymptotic analysis of system (1.3): Theorem 2 on page 34. We build an asymptotic expansion for the solutions of the compressible system (1.3) when ε → 0 and prove quantitative estimates on the errors. We consider well-prepared initial data, i.e. close to the kernel of the penalization operator. The scaling

Ro = Ma = Fr = ε which is considered in (1.
3) is the richest scaling, in the sense that the Coriolis force, the pressure and the gravitational force balance each other. In addition, we analyze the effect of boundary layers on the limiting two-dimensional quasi-geostrophic equation. This limit equation, see (3.36) below, represents the large-scale dynamics of the bulk flow. It is a two-dimensional incompressible Navier-Stokes equation written in terms of the stream function Q, where u = ∇ ⊥ Q is the limit velocity. Frictional effects dissipate energy in the Ekman boundary layer flow, so a damping term appears in (3.36). Such effects were already pointed out for incompressible [START_REF] Masmoudi | Ekman layers of rotating fluids: the case of general initial data[END_REF][START_REF] Chemin | Ekman boundary layers in rotating fluids[END_REF][START_REF] Chemin | Mathematical geophysics. An introduction to rotating fluids and the Navier-Stokes equations[END_REF][START_REF] Gérard-Varet | Highly rotating fluids in rough domains[END_REF] or compressible [START_REF] Bresch | Rotating fluids in a cylinder[END_REF] fluids in high rotation.

As far as we know, that paper [START_REF] Bresch | Rotating fluids in a cylinder[END_REF] is the only work dealing with the asymptotic analysis of compressible fluids in high rotation in the presence of Ekman layers. Our result extends the state of the art in two main directions: (i) we take into account stratification effects due to gravitation, hence we handle non constant reference density states; (ii) we consider general pressure laws P (ρ) ∼ ρ γ , with γ ≥ 3/2 (see (1.10) below). Concerning (i), let us stress that our work seems the first one able to tackle the combined effect of the gravitational force and boundary layers. As for (ii), work [START_REF] Bresch | Rotating fluids in a cylinder[END_REF] handles the case γ = 2. Incidentally, the threshold 3/2 for the number γ is the same as the one for the Lions-Feireisl theory of weak solutions theory.

Our approach is based on relative entropy estimates for system (1.3), see [START_REF] Germain | Weak-strong uniqueness for the isentropic compressible Navier-Stokes system[END_REF][START_REF] Feireisl | Relative entropies, suitable weak solutions, and weak-strong uniqueness for the compressible Navier-Stokes system[END_REF][START_REF] Feireisl | Sun: Suitable weak solutions to the Navier-Stokes equations of compressible viscous fluids[END_REF][START_REF] Feireisl | Rotating compressible fluids under strong stratification[END_REF]. The progress achieved in this paper is made possible thanks to the introduction of a simple tool, which seems to be new in this context: we rely on anisotropic Sobolev embeddings, see Lemma 3.11 below. This enables to compensate for the lack of coercivity for ∂ 3 u h , due to the strong anisotropy in the Lamé operator (1.9). Doing so, we are able to extend the range of values for the parameter γ, and also to improve the quantitative error bounds.

1.4. Outline of the paper. The paper consists of two parts. The first one, treated in Subsection 2, is devoted to the proof of the well-posedness result for system (1.2). The main result is Theorem 1 on page 8. The proof relies on maximal regularity estimates for a parabolic equation related to an anisotropic and variable coefficients Lamé operator; see Proposition 2.2. The second part, which is the matter of Subsection 3, is concerned with the proof of the quantitative estimates for (1.3) in the limit ε → 0. We first build an expansion based on a formal multi-scale analysis. Second, we derive a relative entropy inequality. Finally, we carry out the quantitative estimates, using in particular the anisotropic Sobolev embeddings of Lemma 3.11.

1.5. Main notations and definitions. Since our interest is in real fluid flows, the whole paper is written in space dimension d = 3. Notice though that some results, in particular those of Subsection 2 such as the maximal regularity statement, can easily be extended to higher-dimensional systems. The domain Ω denotes an open set, usually Ω = R 3 or R 2 × (0, 1) in this paper.

When appropriate, we use Einstein's convention on repeated indices for summation. Given a two-dimensional vector v = (v 1 , v 2 ), we define v ⊥ := (-v 2 , v 1 ). Given a vector field v ∈ R 3 , we will often use the notations ∇•v and ∇×v to denote respectively div v and curl v. For a vector x ∈ R 3 , we often use the notation x = (x h , x 3 ) ∈ R 3 to denote the horizontal component x h ∈ R 2 and the vertical component x 3 ∈ R. According to this decomposition, we define the horizontal differential operators ∇ h , ∆ h and ∇ h • as usually; we also set ∇ ⊥ h := (-∂ 2 , ∂ 1 ). These operators act just on the x h variables. Notice that the third component of the vector

∇ × v is ∂ 1 v 2 -∂ 2 v 1 :
this quantity will be denoted by ∇ ⊥ h • v h . We introduced ∆ µ,ε as in (1.1). Similarly, we define the modified gradient operator

∇ µ,ε := √ µ∂ 1 , √ µ∂ 2 , √ ε∂ 3 . The anisotropic Lamé operator L is defined by (1.9) Lu = -∆ µ,ε u -λ∇∇ • u .
Throughout this paper, given a Banach space X and a sequence (a ε ) ε of elements of X, the notation (a ε ) ε ⊆ X is to be understood as the fact that the sequence (a ε ) ε is uniformly bounded in X. We will often denote, for any p ∈ [1, ∞] and any Banach space X, L p T (X) := L p (0, T ); X(Ω) . When T = +∞, we will simply write L p (X). For the definition and basic properties of Besov spaces, we refer to Chapter 2 of [START_REF] Bahouri | Fourier Analysis and Nonlinear Partial Differential Equations[END_REF], or to Subsections 2.2 and 2.3 of [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF].

Pressure law. We consider barotropic flows, for which the pressure P is supposed to be a smooth function of the density only. We assume (see e.g. [START_REF] Feireisl | Sun: Suitable weak solutions to the Navier-Stokes equations of compressible viscous fluids[END_REF], [START_REF] Feireisl | A singular limit for compressible rotating fluids[END_REF], [START_REF] Feireisl | Relative entropies, suitable weak solutions, and weak-strong uniqueness for the compressible Navier-Stokes system[END_REF], [START_REF] Feireisl | Multiple scales and singular limits for compressible rotating fluids with general initial data[END_REF]) that P ∈ C [0, ∞) ∩ C 3 (0, ∞) enjoys (1.10) P (0) = 0 , P (ρ) > 0 ∀ ρ > 0 , lim ρ→∞ P (ρ) ρ γ-1 = a > 0 , for some γ ≥ 1. Given P , we define the internal energy function H by the formula

(1.11) H(ρ) := ρ ˆρ 1 P (z) z 2 dz for all ρ ∈ (0, ∞) .
Notice that the relation ρ H (ρ) = P (ρ) holds for all ρ > 0.

2. A well-posedness result in the presence of a strongly anisotropic viscous stress tensor and stratification

In this section we show a well-posedness result for the barotropic Navier-Stokes system (1.2). As explained in Subection 1.2, there are two main difficulties. The first one is to handle the anisotropy of the viscous stress tensor. It prevents one from using classical compactness techniques to prove the existence of weak solutions having finite energy (see e.g. [START_REF] Lions | Compressible models[END_REF][START_REF] Novotny | Introduction to the mathematical theory of compressible flow[END_REF][START_REF] Feireisl | On the existence of globally defined weak solutions to the Navier-Stokes equations of compressible isentropic fluids[END_REF] and the comments in [START_REF] Bresch | Rotating fluids in a cylinder[END_REF][START_REF] Bresch | Global existence of weak solutions for compressible Navier-Stokes equations: thermodynamically unstable pressure and anisotropic viscous stress tensor[END_REF]). It is also a major obstacle to the use of the effective flux (1.6) as in [START_REF] Hoff | Global solutions of the Navier-Stokes equations for multidimensional compressible flow with discontinuous initial data[END_REF][START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF]. The second one is the non-constant reference density state ρ, due to the potential force ∇G.

Our approach is reminiscent from the works [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF][START_REF] Shibata | Global existence of classical solutions and optimal decay rate for compressible flows via the theory of semigroups[END_REF]. It is based on maximal regularity estimates for the velocity field. This approach enables us to fully exploit the parabolic gain of regularity due to the momentum equation, and to use it in the mass equation in order to transport higher-order Sobolev norms of the density function. Moreover, it allows us to consider non constant density reference states, which is crucial in view of studying stratification effects, see system (1.3) and Section 3.

It is not clear whether or not the whole method of the paper [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF] works in presence of an anisotropic Lamé operator, since it deeply uses Hoff's effective flux (1.6) and algebraic cancellations appearing in its equation, see (1.7). Thus, compared to [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF], we will work with solutions in the energy space, so with less integrability, but which are more regular. We are also able to deal with parabolic operators with variable coefficients when applying the maximal regularity results.

Our main result in this direction is the following statement. Our aim is to give a streamlined method for well-posedness, appropriate for proving the existence and uniqueness of finite-energy solutions in the presence of: (i) an anisotropic viscous stress tensor given by (1.1), (ii) non constant reference density states. A small compromise consists on the one hand on the fact that we do not strive for optimality in the assumptions of the theorem, and on the other hand on the fact that we work in the space domain Ω = R 3 . Notice however that the case of more general domains can be treated with some technical adaptations.

Theorem 1. Let γ ≥ 1. Let ρ ∈ W 3,∞ (R).
Assume that ρ is uniformly bounded from below, i.e. ρ ≥ κ > 0. We define the potential by G = H (ρ), where H is defined by (1.11). Then for any η > 0 which verifies

(2.1) η ≤ min(1/(8C 0 ), κ/8),
where C 0 is the constant given by Proposition 2.2, the property below holds.

Consider system (1.2), supplemented with the initial datum ρ, u |t=0 = ρ in , u in .

For any ρ in , u in , with ρ in -ρ ∈ H 2 and u in ∈ B 3/2 2,4/3 and such that

ρ in -ρ L ∞ ≤ η , there exist a time T * (γ, µ, ε, λ, ρ W 3,∞ , κ, u in B 3/2 2,4/3 , ρ in -ρ H 2 ) > 0 and a unique solution ρ, u to (1.2) on [0, T * ] × Ω, such that: (1) ρ -ρ ∈ C T * (H 2 ), with ρ -ρ L ∞ T * (L ∞ ) ≤ 4 η; (2) u ∈ L ∞ T (L 2 ) ∩ L 2 T * (L ∞ ), with in addition ∇u ∈ L 4 T * (L 2 ) ∩ L 2 T * (L ∞ ), ∇ 2 u ∈ L 4 T * (L 2 ), ∇ 3 u ∈ L 4/3 T * (L 2 ) and ∂ t u ∈ L 4/3
T * (H 1 ); (3) ρ, u satisfies the classical energy inequality (see estimate (3.42) below). Remark 2.1. Notice that we have some freedom on the regularity of the initial datum. In this sense, we take B 3/2 2,4/3 regularity for u 0 for simplicity of presentation, but this condition can be somehow weakened. Notice also that H s → B 3/2 2,4/3 for any s > 3/2. The rest of this section is devoted to the proof of Theorem 1. The proof is based only on elementary methods, namely energy estimates and maximal regularity. We will limit ourselves to showing a priori estimates for smooth solutions to system (1.2). The proof of existence by an approximation scheme is rather standard, see e.g. [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF] and references therein.

2.1. Maximal regularity for an anisotropic Lamé operator. Here we prove the following maximal regularity result for a parabolic equation with an anisotropic Lamé operator and vertical stratification due to the coefficient ρ. Pay attention to the fact that, in the definition of W 2,1 p 2 ,r 2 below, the indices for space and time are in reverse order, in order to stick to the classical definition (see e.g. [START_REF] Seregin | New version of the Ladyzhenskaya-Prodi-Serrin condition[END_REF]).

Proposition 2.2. Let ρ ∈ W 1,∞ (R). Assume that ρ is uniformly bounded from below, i.e. ρ ≥ κ > 0. Let (p j , r j ) j=0,1,2 satisfy 1 < p 2 , r 2 < +∞, r 2 < r 0 < +∞, r 2 < r 1 < +∞, p 0 ≥ p 2 , p 1 ≥ p 2 , together with the relations 2 r 2 + 3 p 2 = 1 + 2 r 1 + 3 p 1 and 2 r 2 + 3 p 2 = 2 + 2 r 0 + 3 p 0 .
Let h in be in Ḃs 2 p 2 ,r 2 with s 2 := 2 -2/r 2 , and let f be in L r 2 loc R + ; L p 2 (R d ) . Then, there exists a constant C 0 = C 0 (µ, ε, λ, r 0 , p 0 , r 1 , p 1 , r 2 , p 2 , ρ W 1,∞ , κ) > 0 such that, for all T > 0 and all h ∈ L r 2 T (W 1,p 2 ), with h solution to the Lamé system

(2.2) ρ(x 3 )∂ t h + Lh = f h |t=0 = h in
where L is defined by (1.9), one has the properties h ∈ L r 0 [0, T ]; L p 0 and ∇h ∈ L r 1 [0, T ]; L p 1 , together with the estimate

h L ∞ T ( Ḃs 2 p 2 ,r 2 ) + h L r 0 T (L p 0 ) + ∇h L r 1 T (L p 1 ) + ∂ t h, ∇ 2 h L r 2 T (L p 2 ) ≤ C 0 h in Ḃs 2 p 2 ,r 2
+ f L r 2 T (L p 2 ) + h L r 2 T (L p 2 ) + ∇h L r 2 T (L p 2 ) . (2.3)
To prove this result, we follow the idea of [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF]. We apply the Leray projector P to the system (2.2) and rely on the maximal regularity for a divergence-form parabolic equation. Because of the stratification, additional commutators involving P and the reference density ρ have to be analyzed. We do not strive for making the dependence of C 0 explicit in ε, µ, λ, ρ W 1,∞ and κ, because our aim is to obtain a well-posedness theorem for a fixed set of parameters, see Theorem 1.

Remark 2.3. Notice that the two last terms in the right hand side of estimate (2.3) cannot be get rid of. Their appearance is due to the commutator terms involving the Leray projector. It is not possible to swallow them in the left hand side of (2.3). In spite of this, estimate (2.3) can be used as such to prove the well-posedness result stated in Theorem 1. Indeed, the terms h L r 2 T (L p 2 ) and ∇h L r 2 T (L p 2 ) are of lower order and can be handled by interpolating with the finite energy. Such an analysis will also be done for the source term f , which contains in particular the nonlinear term ρu • ∇u, see below in Section 2.4.

Proof of Proposition 2.2. We proceed in several steps.

Step 1: reduction to the heat equation. We first rewrite system (2.2) in divergence form: we have

∂ t h -∇ µ,ε • 1 ρ(x 3 ) ∇ µ,ε h -λ∇ 1 ρ(x 3 ) ∇ • h = F , (2.4)
where, for simplicity of notation, we have defined

F := 1 ρ(x 3 ) f + ε ρ (x 3 ) ρ(x 3 ) 2 ∂ 3 h + λ ρ (x 3 ) ρ(x 3 ) 2 ∇ • h e 3 . (2.5)
As in [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF], we apply the Leray projector P on the equation. The difference with [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF] is that we now have commutator terms appearing. Recall that (2.6)

P = Id -Q , Q := -∇(-∆) -1 ∇•
, where the formulas have to be interpreted in the sense of Fourier multipliers. Applying P to (2.4), we get

∂ t Ph -∇ µ,ε • 1 ρ(x 3 ) ∇ µ,ε Ph = PF + C , (2.7)
where the commutator term C is defined by

C := -∇ µ,ε • 1 ρ(x 3 ) , P ∇ µ,ε h . (2.8)
Here above, the symbol [A, B] := AB -BA denotes the commutator between two operators A and B. We now compute the equation for Qh = h -Ph. Notice that, using (2.6), we have

∇∇ • h = ∇∇ • Qh = ∆Qh Therefore, rewriting (2.4) in the form ∂ t h -∇ µ,ε • 1 ρ(x 3 ) ∇ µ,ε h -λ 1 ρ(x 3 ) ∇∇ • h = 1 ρ(x 3 ) f + ε ρ (x 3 ) ρ(x 3 ) 2 ∂ 3 h
and taking the difference of this equation with (2.7), we immediately find

∂ t Qh -∇ µ+λ,ε+λ • 1 ρ(x 3 ) ∇ µ+λ,ε+λ Qh = Q 1 ρ(x 3 ) f -P ε ρ (x 3 ) ρ(x 3 ) 2 ∂ 3 h + λ ρ (x 3 ) ρ(x 3 ) 2 ∇ • h e 3 -C + λ ρ (x 3 ) ρ(x 3 ) 2 ∂ 3 Qh.
(2.9)

Both equations (2.7) and (2.9) are heat-type equations with a variable coefficient. The next step will show that all the quantities appearing in their right-hand sides are of lower order in h.

Step 2: computation of the commutators. Our goal is now to compute the commutator term C, defined by (2.8). Owing to (2.6) and switching the position of the constant factors, we have

C := -∇ µ,ε • 1 ρ(x 3 ) , P ∇ µ,ε h = ∇ • 1 ρ(x 3 ) , Q ∇ µ 2 ,ε 2 h .
Recall the definition of Q in (2.6). We start by applying the divergence operator to the various quantities: after observing that

∇ • ∇ µ 2 ,ε 2 = ∆ µ,ε , we find 1 ρ , Q ∇ µ 2 ,ε 2 h = - 1 ρ ∇(-∆) -1 ∆ µ,ε h + ∇(-∆) -1 ∇ • 1 ρ ∇ µ 2 ,ε 2 h .
Therefore, we find

C = ∇ • 1 ρ(x 3 ) , Q ∇ µ 2 ,ε 2 h = -∇ • 1 ρ(x 3 ) ∇(-∆) -1 ∆ µ,ε h + ∇ • ∇(-∆) -1 ∇ • 1 ρ(x 3 ) ∇ µ 2 ,ε 2 h = ρ (x 3 ) ρ(x 3 ) 2 ∂ 3 (-∆) -1 ∆ µ,ε h + 1 ρ(x 3 ) ∆ µ,ε h -∇ • 1 ρ(x 3 ) ∇ µ 2 ,ε 2 h , which in the end leads us to (2.10) C = ρ (x 3 ) ρ(x 3 ) 2 ∂ 3 (-∆) -1 ∆ µ,ε h + ε ∂ 3 h .
Since the operator (-∆) -1 ∆ µ,ε is a singular integral operator, whose symbol is homogeneous of degree 0, the previous computations show that C is indeed a lower order term, as claimed.

Step 3: estimates via maximal regularity. We apply maximal regularity estimates of [START_REF] Hieber | Prüss: Heat kernels and maximal L p -L q estimates for parabolic evolution equations[END_REF] for the divergence-form parabolic operator

(2.11) ∂ t ( • ) -∇ µ,ε • ρ -1 ∇ µ,ε ( • )
to equation (2.7). We obtain

∂ t Ph, ∇ 2 Ph L r 2 T (L p 2 ) ≤ C Ph in Ḃs 2 p 2 ,r 2
+ PF + C L r 2 T (L p 2 )
. Using the continuity of P and (-∆) -1 ∆ µ,ε on L p 2 (since 1 < p 2 < +∞) and keeping in mind definition (2.5) of F and (2.10), it is easy to bound

PF + C L r 2 T (L p 2 ) ≤ C f L r 2 T (L p 2 ) + ∇h L r 2 T (L p 2 )
. So, in the end we find

∂ t Ph, ∇ 2 Ph L r 2 T (L p 2 ) ≤ C 0 Ph in Ḃs 2 p 2 ,r 2 + f L r 2 T (L p 2 ) + ∇h L r 2 T (L p 2 ) . (2.12)
Next, thanks to the gaussian bounds on the fundamental solution of the parabolic equation (2.11), see [START_REF] Nash | Continuity of solutions of parabolic and elliptic equations[END_REF] and [START_REF] Aronson | Bounds for the fundamental solution of a parabolic equation[END_REF], arguing as in Lemma 2.4 of [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF] we infer that (2.13)

Ph L r 0 T (L p 0 ) ≤ C 0 Ph in Ḃs 2 p 2 ,r 2 + f L r 2 T (L p 2 ) + ∇h L r 2 T (L p 2
) . Thus, in order to complete the proof, it remains us to bound the L r 1 T (L p 1 ) norm of ∇Ph. This is the goal of the next step.

Step 4: a functional inequality. We claim that, for any w ∈ W 2,1 p 2 ,r 2 R 3 × (0, T ) , one has the estimate

(2.14) ∇w L r 1 T (L p 1 ) ≤ C ∇ 2 w L r 2 T L p 2 (R 3 ×(0,T )) + ∂ t w L r 2 T L p 2 (R 3 ×(0,T )) .
In fact, estimate (2.14) is a functional inequality, which does not use the parabolic equation. We are going to deduce it from the corresponding inequality in a bounded domain, see in particular [START_REF] Solonnikov | Estimates of solutions of the Stokes equations in S. L. Sobolev spaces with a mixed norm[END_REF]Proposition 2.1]. In order to prove (2.14), we take

ϕ ∈ C ∞ c (B(0, 1)) such that ϕ = 1 on B(0, 1 2 ). For k ∈ N, we consider w k := ϕ(•/k)w. Notice that w k is supported in B(0, k).
Then, we consider the rescaled function w λ k for λ ≥ k, according to the parabolic scaling:

w λ k = w k (λ•, λ 2 •) = ϕ( λ k •)w(λ•, λ 2 •). Notice that w λ k is supported in B(0, 1)
. By the mixed norm parabolic Sobolev embedding in the bounded domain B(0, 1), we have

∇w λ k L r 1 T L p 1 (B(0,1)×(0,T /λ 2 )) ≤ C w λ k W 2,1 p 2 ,r 2 (B(0,1)×(0,T /λ 2 )) ,
where C is a constant independent of T . Then, by rescaling, we obtain

λ λ 3 p 1 + 2 r 1 ∇w k L r 1 T L p 1 (B(0,k)×(0,T )) ≤ C 1 λ 3 p 2 + 2 r 2 w k L r 2 T L p 2 (B(0,k)×(0,T )) + λ λ 3 p 2 + 2 r 2 ∇w k L r 2 T L p 2 (B(0,k)×(0,T )) + λ 2 λ 3 p 2 + 2 r 2 ∇ 2 w k L r 2 T L p 2 (B(0,k)×(0,T )) + λ 2 λ 3 p 2 + 2 r 2 ∂ t w k L r 2 T L p 2 (B(0,k)×(0,T ))
, and letting λ → ∞ we deduce

∇w k L r 1 T L p 1 (R 3 ×(0,T )) (2.15) ≤ C ∇ 2 w k L r 2 T L p 2 (R 3 ×(0,T )) + ∂ t w k L r 2 T L p 2 (R 3 ×(0,T )) .
From estimate (2.15) and the fact that w ∈ W 2,1 p 2 ,r 2 , we get that ∇w k is a Cauchy sequence in L r 1 T (L p 1 ). Hence we can pass to the limit in (2.15): thus, we obtain (2.14), as claimed.

Step 5: end of the proof. Owing to the fact that 1 < p 2 < +∞ and to (2.12), we have that Ph belongs to W 2,1 p 2 ,r 2 (R 3 × (0, T )). Therefore, we can apply inequality (2.14) to w = Ph. We find

(2.16) ∇Ph L r 1 T (L p 1 ) ≤ C 0 Ph in Ḃs 2 p 2 ,r 2 + f L r 2 T (L p 2 ) + ∇h L r 2 T (L p 2 )
, for a possibly new constant C 0 > 0. In addition, since Qh solves equation (2.9), which is analogous to (2.7), estimates similar to (2.12), (2.13) and (2.16) hold true for Qh and its derivatives. Then, writing h = Ph + Qh completes the proof of the proposition.

Basic energy estimates.

Let us take a smooth solution (ρ, u) to system (1.2), such that ρ -ρ and u decay sufficiently fast at space infinity. We want to find a priori estimates in suitable norms. For this, we are going to work with the variables r(t) := ρ(t) -ρ and u.

In the same way, we set r in = ρ in -ρ. Recall that (2.17)

r in L ∞ ≤ η .
We start by performing classical energy estimates, which provides us with a bound for the low frequencies of the velocity field. Namely, by multiplying the momentum equation in (1.2) by u and integrating by parts, we get, in a standard way, the control

(2.18) √ ρ u L ∞ (L 2 ) + ∇u L 2 (L 2 ) ≤ C energy .
See Section 3.2 for similar bounds. Of course, the constant C energy depends also on the (fixed) values of the coefficients (µ, ε, λ). From this control and

ˆΩ ρ|u| 2 = ˆΩ ρ|u| 2 -ˆΩ r|u| 2 ,
we infer that, for any T > 0,

(2.19) u L ∞ T (L 2 ) ≤ κ -1 C energy + r L ∞ T (L ∞ ) u L ∞ T (L 2 )
, where we recall that κ is a lower bound for ρ. Hence, if T > 0 is such that

(2.20) r L ∞ T (L ∞ ) ≤ 4 η < κ 2 ,
where η > 0 is the size of the initial datum in the L ∞ norm, recall (2.17), we deduce

(2.21) u L ∞ T (L 2 ) ≤ 2C energy κ -1 .
The next goal is to exhibit a control on the density variation function r. We will work in higher Sobolev norms, namely in H 2 . However, we are going to bound its L ∞ norm independently (i.e. without using Sobolev embeddings), in order to get, in view of (2.20), a smallness condition only on r in L ∞ , and not on the higher order norm of r in .

2.3. Estimates for the density function. In this subsection, we find transport estimates for the density variation function r. First of all, from the mass equation in (1.2) we find that r fulfills

(2.22) ∂ t r + u • ∇r + r ∇ • u = -u 3 ρ -ρ ∇ • u .
A basic L p estimate for this equation gives, for any t > 0,

r(t) L p ≤ r in L p + 1 - 1 p ˆt 0 r(τ ) L p ∇ • u(τ ) L ∞ dτ + ˆt 0 ρ L ∞ u 3 (τ ) L p + ∇ • u(τ ) L p dτ .
On the one hand, in the case p = 2, we get

r L ∞ T (L 2 ) ≤ r in L 2 + 1 2 ˆT 0 r L 2 ∇ • u L ∞ dτ + ˆT 0 ρ L ∞ u 3 (τ ) L 2 + ∇ • u(τ ) L 2 dτ (2.23)
for any time T > 0. In turn, Grönwall's lemma gives the bound

(2.24) r L ∞ T (L 2 ) ≤ e ´T 0 ∇•u(τ ) L ∞ dτ × r in L 2 + ˆT 0 ρ L ∞ u 3 (τ ) L 2 + ∇ • u(τ ) L 2 dτ .
On the other hand, by letting p → +∞, from Grönwall's lemma again we deduce

(2.25) r L ∞ T (L ∞ ) ≤ e ´T 0 ∇•u(τ ) L ∞ dτ × r in L ∞ + ˆT 0 ρ L ∞ u 3 (τ ) L ∞ + ∇ • u(τ ) L ∞ dτ .
We now define the time T > 0 as

(2.26) T := sup t > 0 ˆt 0 ρ L ∞ u 3 (τ ) L ∞ dτ + ˆt 0 (1 + ρ L ∞ ) ∇ • u(τ ) L ∞ dτ ≤ min η, log 2 ,
where η > 0 is as in (2.17). From the previous bound we gather then

(2.27) r L ∞ T (L ∞ ) ≤ 2 r in L ∞ + 2 η ≤ 4 η .
We now differentiate equation (2.22) with respect to x j , for any j ∈ {1, 2, 3}. We get that ∂ j r verifies the following equation:

(2.28) ∂ t ∂ j r + u • ∇∂ j r + ∂ j r ∇ • u = -∂ j u • ∇r -r ∇ • ∂ j u -∂ j u 3 ρ -u 3 ∂ j ρ -∂ j ρ∇ • u -ρ∇ • ∂ j u .
The same L 2 energy estimate for the continuity equation as above yields

∂ j r(t) L 2 ≤ ∂ j r in L 2 + 1 2 ˆt 0 ∂ j r L 2 ∇ • u L ∞ dτ + ˆt 0 r L ∞ ∇ • ∂ j u L 2 + ∇r L 2 ∂ j u L ∞ + ∂ j ρ L ∞ u 3 L 2 + ρ L ∞ ∂ j u 3 L 2 + ∂ j ρ L ∞ ∇ • u L 2 + ρ L ∞ ∇ • ∂ j u L 2 dτ ,
for any t ≥ 0. Hence, for all t ≥ 0, one has

∇r(t) L 2 ≤ ∇r in L 2 (2.29) + C ˆt 0 r L ∞ ∇ 2 u L 2 + ∇r L 2 ∇u L ∞ + ρ W 2,∞ u H 2 dτ .
We do not apply Grönwall's lemma directly on this inequality. It is better to bound first the second-order derivatives of r. For this, let us differentiate equation (2.28) with respect to x k , for any k ∈ {1, 2, 3}. We deduce the following equation for the quantity r jk := ∂ 2 kj r:

∂ t r jk + u • ∇r jk + r jk ∇ • u = -∂ 2 jk u • ∇r + ∂ j u • ∇∂ k r + ∂ k r • ∇∂ j u -r∇ • ∂ 2 kj u -∂ k ∂ j u 3 ρ + u 3 ∂ j ρ + ∂ j ρ∇ • u + ρ∇ • ∂ j u .
Performing another energy estimate, we infer, for any t ≥ 0, the inequality

r jk (t) L 2 ≤ r in,jk L 2 + 1 2 ˆt 0 r jk L 2 ∇ • u L ∞ dτ + ˆt 0 ∇r L 4 ∇ 2 u L 4 + ∇u L ∞ ∇ 2 r L 2 + r L ∞ ∇ 3 u L 2 + ρ W 3,∞ u H 3 dτ ,
where we have defined r in,jk := ∂ 2 jk r in . After using the interpolation inequality

(2.30) φ L 4 ≤ φ 1/4 L 2 φ 3/4 L 6
together with the Sobolev embedding Ḣ1 (R 3 ) → L 6 (R 3 ), we obtain

∇r L 4 ∇ 2 u L 4 ≤ ∇r 1/4 L 2 ∇ 2 r 3/4 L 2 ∇ 2 u 1/4 L 2 ∇ 3 u 3/4 L 2 ≤ ∇r L 2 + ∇ 2 r L 2 ∇ 2 u L 2 + ∇ 3 u L 2 .
In view of this inequality, the previous bound on r jk yields, on the time interval [0, T ],

∇ 2 r L ∞ T (L 2 ) ≤ ∇ 2 r in L 2 + C ˆt 0 r L ∞ ∇ 3 u L 2 + ρ W 3,∞ u H 3 + ∇r L 2 + ∇ 2 r L 2 ∇u L ∞ + ∇ 2 u L 2 + ∇ 3 u L 2 dτ .
(2.31)

Let us now introduce the notation

R(t) := r L ∞ t (H 2 ) + r L ∞ t (L ∞ ) .
Summing up (2.24), (2.27), (2.29) and (2.31) we deduce the following bound:

(2.32) R(T ) ≤ C r in H 2 + ˆT 0 ρ W 3,∞ + R(τ ) ( ∇u L ∞ + u H 3 ) dτ .
To end this part, we observe that, for almost all t ∈ [0, T ] and x ∈ R 3 , |r(x, t)| ≤ 4η. So, by assumptions (1.10) on the pressure, there exists a constant C > 0, depending on the function P and ρ L ∞ , such that, for any parameter s ∈ [0, 1], on [0, T ] × R 3 one has (2.33) |P (ρ + rs)| + P (ρ + rs) + P (ρ + rs) + P (ρ + rs) ≤ C .

2.4.

Maximal regularity estimates for the velocity field. It remains to find a bound on u. For this, we mimic an approach used in [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF], based on the maximal regularity estimates of Subsection 2.1. To begin with, we recast the equation for u as (2.34)

ρ∂ t u + Lu = -r ∂ t u -f ,
where the anisotropic Lamé operator L is defined by (1.9), and where we have set

(2.35) f := (ρ + r) u • ∇u + ∇P (ρ) -ρ∇G.
In the computations below, when convenient, we will resort to the notation ρ + r = ρ, and use the bounds (2.18) and (2.27) for √ ρu and √ ρ respectively.

In view of (2.32) above, we are interested in H 3 bounds for u: for this, we will apply Proposition 2.2 to both u and ∇u. To this end, we fix the following values of the parameters:

p 2 , r 2 = (2, 4/3) , (p 0 , r 0 ) = (+∞, 2) , (p 1 , r 1 ) = (2, 4) . (2.36)
With these choices, all the hypotheses in Proposition 2.2 are satisfied. Indeed, thanks to the energy estimates (2.18) and (2.21) and maximal regularity (2.12), one can easily verify that, a priori, both u and ∇u belong to the space W 2,1 2,4/3 . We choose p 0 = +∞ in order to have, thanks to (3.53) above, a control on u in any L p , p ∈ [2, +∞]. Notice also that we have some freedom for the values of r 2 and p 1 , which then determine r 0 and r 1 . Here, we take the simple choice r 2 = 4/3 and p 1 = 2. This implies s 2 = 1/2. For t ≥ 0, let us introduce the quantity

U(t) := u L ∞ t (L 2 )∩L 2 t (L ∞ ) + ∇u L 4 t (L 2 )∩L 2 t (L ∞ ) + ∇ 2 u L 4 t (L 2 ) + ∇ 3 u L 4/3 t (L 2 ) + ∂ t u L 4/3 t (H 1 ) . On [0, T ],
where T > 0 is the time defined in (2.26), we have

(2.37) U(T ) ≤ C 0 u in Ḃ1/2 2,4/3 + ∇u in Ḃ1/2 2,4/3 + r∂ t u + f L 4/3 T (H 1 ) + u L 4/3 T (L 2 ) + ∇u L 4/3 T (L 2 ) + ∇ 2 u L 4/3 T (L 2 )
, where f is defined in (2.35). Notice that

u in Ḃ1/2 2,4/3 + ∇u in Ḃ1/2 2,4/3 ≤ u in Ḃ1/2 2,4/3 ∩ Ḃ3/2 2,4/3 ≤ u in B 3/2 2,4/3
, where the last inequality holds in view of the fact that B s p,r = Ḃs p,r ∩ L p for any s > 0, see also Chapter 2 of [START_REF] Bahouri | Fourier Analysis and Nonlinear Partial Differential Equations[END_REF].

Next, we bound the term containing the time derivative. Notice that, for any j ∈ {1, 2, 3}, we have

∂ j (r ∂ t u) = r ∂ t ∂ j u + ∂ j r ∂ t u .
Therefore, we have

r ∂ t u L 4/3 T (L 2 ) + r ∂ t ∇u L 4/3 T (L 2 ) ≤ 4 η ∂ t u L 4/3 T (H 1 )
, so this term can be absorbed in the left-hand side of (2.37), if η > 0 is fixed so that condition (2.1) is fulfilled. As for the remaining term ∂ j r ∂ t u, we notice first of all that from (2.27), the lower bound for ρ, the momentum equation in (1.2) and (2.18), we can bound

∂ t u L 2 T (L 2 ) ≤ C ρ∂ t u L 2 T (L 2 ) ≤ C Lu + ρu • ∇u + ∇P -ρ∇G L 2 T (L 2 )
. At this point, we obviously have

Lu L 2 T (L 2 ) ≤ T (1/2)-(1/4) Lu L 4 T (L 2 ) ≤ C T 1/4 U(T ). In addition, one has ρu • ∇u L 2 T (L 2 ) ≤ C T (1/2)-(1/2) √ ρ u L ∞ T (L 2 ) ∇u L 2 T (L ∞ ) ≤ C U(T ) . (2.38)
Let us now turn to the pressure and potential force terms. In view of (1.8), we have

(2.39) ∇P (ρ) -ρ∇G = ∇ P (ρ) -P (ρ) -(ρ -ρ)∇G = ∇ P (ρ + sr)r -r∇G ,
for some s ∈ ]0, 1[ . On the one hand, since G = H (ρ), G is bounded; so, for a constant C > 0 depending on ρ W 1,∞ and on the function P , we have

r∇G L 2 T (L 2 ) ≤ C T 1/2 r L ∞ T (L 2
) . On the other hand, by writing (2.40)

∇ P (ρ + sr)r = P (ρ + sr)∇r + P (ρ + sr)r∇ρ + sP (ρ + sr)r∇r and making use of (2.33) and (2.27), direct computations show that

∇ P (ρ + sr)r L 2 T (L 2 ) ≤ C T 1/2 r L ∞ T (H 1 )
. Here, the constant C > 0 depends on the constant appearing in (2.33), on η and ρ W 1,∞ . Assuming without loss of generality that T ≤ 1, we have thus proved that (2.41)

∂ t u L 2 T (L 2 ) ≤ C R(T ) + U(T ) .
With this control at hand, let us return to the bound of ∂ j r ∂ t u in L 4/3 T (L 2 ), for j ∈ {1, 2, 3}. By resorting again to the interpolation inequality (2.30), in view of the Sobolev embedding Ḣ1 → L 6 and of Young's inequality we get

∂ j r ∂ t u L 4/3 T (L 2 ) ≤ ∇r 1/4 L ∞ T (L 2 ) ∇ 2 r 3/4 L ∞ T (L 2 ) ∂ t u 1/4 L 4/3 T (L 2 ) ∂ t ∇u 3/4 L 4/3 T (L 2 ) ≤ C R(T ) T (3/4)-(1/2) ∂ t u L 2 T (L 2 ) 1/4 ∂ t ∇u 3/4 L 4/3 T (L 2 ) ≤ C T 1/16 R(T ) R(T ) + U(T ) 1/4 U(T ) 3/4 ≤ C T 1/16 R(T ) + U(T ) 2 .
Inserting all the previous bounds in (2.37), under condition (2.1) we deduce that

U(T ) ≤C u in B 3/2 2,4/3 + T 1/16 R(T ) + U(T ) 2 + f L 4/3 T (H 1 ) + u L 4/3 T (H 2 ) ≤C u in B 3/2 2,4/3 + T 1/16 1 + R(T ) + U(T ) 2 + f L 4/3 T (H 1 ) .
(2.42)

Our next goal is to bound f in the L 4/3
T (H 1 ) norm. First of all, let us focus on the L 4/3 T (L 2 ) norm. Repeating the computations which led to (2.41), we easily find that

f L 4/3 T (L 2 ) ≤ C T 3/4 R(T ) + T 1/4 U(T ) ≤ C T 1/4 R(T ) + U(T )
, where again, without loss of generality, we have assumed that T ≤ 1. Next, for j ∈ {1, 2, 3}, from the definition (2.35) of f , we get

∂ j f = ∂ j ρ u • ∇u + ρ ∂ j u • ∇u + ρ u • ∇∂ j u + ∂ j ∇P (ρ) -ρ∇G . (2.43)
We first bound the convective terms in (2.43). By (2.27), interpolation between Lebesgue spaces and (2.30), we have

ρ u • ∇∂ j u L 4/3 T (L 2 ) ≤ C T (3/4)-(1/2) u L 4 T (L 4 ) ∇ 2 u L 4 T (L 4 ) ≤ C T 1/4 u 1/2 L ∞ T (L 2 ) u 1/2 L 2 T (L ∞ ) ∇ 2 u 1/4 L 4 T (L 2 ) ∇ 3 u 3/4 L 4 T (L 2 ) ≤ C T 1/4 U(T ) 2 .
Arguing in a similar way, we obtain

ρ ∂ j u • ∇u L 4/3 T (L 2 ) ≤ C T (3/4)-(1/2) ∇u 2 L 4 T (L 4 ) ≤ C T 1/4 ∇u 1/4 L 4 T (L 2 ) ∇ 2 u 3/4 L 4 T (L 2 ) 2 ≤ C T 1/4 U(T ) 2 .
Finally, writing ρ = ρ + r, we split the term

∂ j ρ u • ∇u = ∂ j ρ u • ∇u + ∂ j r u • ∇u.
We estimate on the one hand, similarly to (2.38),

∂ j ρ u • ∇u L 4/3 T (L 2 ) ≤ CT 1 4 U(T ),
and on the other hand

∂ j r u • ∇u L 4/3 T (L 2 ) ≤ C T (3/4)-(1/4)-(5/12) ∇r L ∞ T (L 4 ) ∇u L 4 T (L 6 ) u L 12/5 T (L 12 ) ≤ C T 1/12 ∇r 1/4 L ∞ T (L 2 ) ∇ 2 r 3/4 L ∞ T (L 2 ) × ∇ 2 u L 4 T (L 2 ) u 1/6 L ∞ T (L 2 ) u 5/6 L 2 T (L ∞ ) ≤ C T 1/12 R(T ) U(T ) 2 .
We now turn to the last term in (2.43). Thanks to (2.39) and (2.40), we have

∂ j ∇P (ρ) -ρ∇G = P (∂ j ρ + s∂ j r)∇r + P ∂ j ∇r + P (∂ j ρ + s∂ j r)r∇ρ + P ∂ j r∇ρ + P r∂ j ∇ρ + sP (∂ j ρ + s∂ j r)r∇r + sP ∂ j r∇r + sP r∂ j ∇r -∂ j r∇G -r∂ j ∇G
where all the functions P , P and P are computed at the point ρ + sr. Making repeated use of (2.33), (2.27) and (2.30),

∂ j ∇P (ρ) -ρ∇G L 4/3 T (L 2 ) ≤ CT 3/4 R(T ) + R(T ) 2 ,
where

C = C(P, ρ W 3,∞ ) > 0.
Putting all those bounds together, assuming again T ≤ 1, we deduce that

∇f L 4/3 T (L 2 ) ≤ C T 1/12 1 + R(T ) + U(T ) 3 ,
which in turn implies that

(2.44) f L 4/3 T (H 1 ) ≤ C T 1/12 1 + R(T ) + U(T ) 3 .
In the end, inserting (2.44) into (2.42), we find

(2.45) U(T ) ≤ C u in B 3/2 2,4/3 + T 1/16 1 + R(T ) + U(T ) 3 .
On the other hand, the integral term in (2.32) can be bounded by

ρ W 3,∞ + R(T ) T 1/2 ∇u L 2 T (L ∞ ) + T 3/4 u L 4 T (H 2 ) + T 1/4 ∇ 3 u L 4/3 T (L 2 ) , which implies that (2.46) R(T ) ≤ C r in H 2 + T 1/4 1 + R(T ) + U(T ) 2 .
Define now, for all t ≥ 0, the quantity

N (t) := R(t) + U(t) .
Summing up estimates (2.45) and (2.46), we infer that

N (T ) ≤ C u in B 3/2 2,4/3 + r in H 2 + T 1/16 1 + N (T ) 3 .
From this inequality, it is a standard matter to deduce the existence of a time

T * (γ, µ, ε, λ, ρ W 3,∞ , κ, u in B 3/2 2,4/3 , r in H 2 ) > 0, with T * ≤ min{1, T }, such that N (T * ) ≤ 2 C u in B 3/2 2,4/3 + r in H 2 .
The a priori estimates are hence proved in the interval [0, T * ].

2.5. Uniqueness. The uniqueness of solutions, claimed in Theorem 1, is a consequence of the following statement.

Proposition 2.4. Let γ ≥ 1. Let ρ ∈ W 3,∞ (R).
Assume that ρ is uniformly bounded from below, i.e. ρ ≥ κ > 0. We define the potential by G = H (ρ), where H is defined by (1.11). Let ρ in , u in be such that ρ in -ρ ∈ H 2 and u in ∈ B 3/2 2,4/3 , with ρ in -ρ L ∞ ≤ η, for some η > 0 satisfying (2.1). Assume that ρ 1 , u 1 and ρ 2 , u 2 are two solutions to system (1.2), related to the same initial datum ρ in , u in and belonging to the space 1)-( 2) of Theorem 1 hold true , for some T > 0.

X T := ρ, u ∈ L ∞ T (L ∞ ) × L ∞ T (L 2 ) Properties (
Then

ρ 1 = ρ 2 and u 1 = u 2 almost everywhere in [0, T ] × Ω.
We show first a simple lemma.

Lemma 2.5. Let ρ in , u in be as in Proposition 2.4 above, and let ρ, u a solution to system (1.2), related to the initial datum ρ in , u in and belonging to the space X T , for some T > 0.

Then ρ ∈ C [0, T ] × Ω and u ∈ C [0, T ]; H 1 (Ω) . In addition, the following estimate holds true, for a "universal" constant C > 0:

u 2 L ∞ T (H 1 ) ≤ C u in 2 B 3/2 2,4/3 + ∂ t u L 2 T (L 2 ) u L 2 T (H 2 )
.

Proof of Lemma 2.5. By definition of the space X T , we know that r := ρ -ρ belongs to C T (H 2 ), which, by Sobolev embeddings, is continuously embedded in C T (C ∩ L ∞ ). Thus ρ = ρ + r is continuous with respect both x and t.

Next, let us consider u. We have shown in (2.41) above that ∂ t u ∈ L 2 T (L 2 ) and that u ∈ L 2 T (H 2 ). From those properties, one easily derives (see e.g. Section 5.9 of [START_REF] Evans | Partial differential equations[END_REF]) that u ∈ C T (H 1 ). The quantitative estimate is a simple consequence of the bound given in [START_REF] Evans | Partial differential equations[END_REF] (see Theorem 3 page 305 therein), combined with the embedding B 3/2 2,4/3 → H 1 . The lemma is hence proved.

We can now prove uniqueness.

Proof of Proposition 2.4. For j = 1, 2, let us define r j := ρ j -ρ. We also set

δr := r 1 -r 2 = ρ 1 -ρ 2 and δu := u 1 -u 2 .
Let us deal with δr first. Its equation reads as follows:

∂ t δr + u 1 • ∇δr + δr ∇ • u 1 = -ρ∇ • δu + δu 3 ρ + δu • ∇r 2 + r 2 ∇ • δu .
An L 2 estimate for this equation (see the beginning of Section 2.3 above) yields, for all t ∈ [0, T ], the bound

δr(t) L 2 ≤ δr in L 2 + 1 2 ˆt 0 δr L 2 ∇ • u 1 L ∞ dτ + ˆt 0 r 2 L ∞ ∇ • δu L 2 + δu L ∞ ∇r 2 L 2 dτ + ˆt 0 ρ L ∞ δu 3 L 2 + ρ L ∞ ∇ • δu L 2 dτ.
We also look for an L 2 bound on ∇δr. To this end, we differentiate the equation for δr with respect to the space variables and we perform an energy estimate. We get

∇δr(t) L 2 ≤ ∇δr in L 2 + 1 2 ˆt 0 ∇δr L 2 ∇ • u 1 L ∞ dτ + ˆt 0 ∇u 1 L ∞ ∇δr L 2 + ∇∇ • u 1 L 4 δr L 4 + ∇∇ • δu L 2 r 2 L ∞ + ∇δu L 4 ∇r 2 L 4 + δu L ∞ ∇ 2 r 2 L 2 + ρ W 2,∞ δu H 2 dτ .
(2.47)

Using the embedding H 1 (R 3 ) → L 4 (R 3 ) (see (2.30) above), and summing up the previous estimate with (2.47), we finally get

δr(t) H 1 ≤ δr in H 1 + C ˆt 0 δr H 1 ∇u 1 L ∞ + ∇ 2 u 1 L 4 dτ (2.48) + C ˆt 0 r 2 L ∞ + ∇r 2 L 4 ∇δu L 2 + ∇ 2 δu L 2 + δu L ∞ ∇ 2 r 2 L 2 + ρ W 2,∞ δu H 2 dτ .
Observe that, by definition of the space X T , we have ∇u 1 ∈ L 2 T (L ∞ ). In addition, from the control (2.30) applied to ∇ 2 u 1 L 4 and Young's inequality, we infer that

∇ 2 u 1 ∈ L 4/3
T (L 4 ). By the same token we get ∇r 2 ∈ L ∞ T (L 4 ), while we already know verifies the following bound, for a suitable exponent α > 0:

that r 2 ∈ L ∞ T (L ∞ ) ∩ L ∞ T (H 2 ).
θ(t) ≤ C t α θ(t) + ˆt 0 ( δu H 2 + δu L ∞ ) dτ , (2.49)
for all t ∈ 0, min{T, 1} . Notice that here we used the fact that δr in = 0.

We now turn to velocity estimates. Taking the difference of the equations for u 1 and u 2 yields the following equation for δu:

ρ 1 ∂ t δu + ρ 1 u 1 • ∇δu + Lδu (2.50) = -δr ∂ t u 2 -∇ P (ρ 1 ) -P (ρ 2 ) + δr∇G -ρ 1 u 1 -ρ 2 u 2 ) • ∇u 2 .
We observe that

ρ 1 u 1 -ρ 2 u 2 ) • ∇u 2 = ρ 1 δu • ∇u 2 + δr u 2 • ∇u 2 ∇ P (ρ 1 ) -P (ρ 2 ) = P (ρ 1 ) ∇δr + P (ζ) δr ∇(ρ + r 2 ) ,
for some ζ = ζ(x, t) in between the values of ρ 1 (x, t) and ρ 2 (x, t). A basic energy estimate for the equation then gives, for any t ∈ [0, T ], the bound 1 2

d dt ρ 1 δu 2 L 2 + ˆΩ |∇ µ,ε δu| 2 + λ ˆΩ |∇ • δu| 2 ≤ δu L 2 δr L 4 ∂ t u 2 L 4 + ρ 1 L ∞ δu L 2 ∇u 2 L ∞ + δr L 4 u 2 L 4 ∇u 2 L ∞ + C ∇δr L 2 + C δr L 2 + C δr L 4 ∇r 2 L 4
, where we have also used the L ∞ T (L ∞ ) boundedness of ρ 1 and ρ 2 in order to control the terms involving derivatives of the pressure function. Let us forget about the viscosity terms for a while. A simple argument allows one to deduce the following control, for any t ∈ [0, T ]:

δu L ∞ t (L 2 ) ≤ C ˆt 0 δu L 2 ∇u 2 L ∞ + θ(τ ) 1 + ∂ t u 2 H 1 + u 2 L 4 ∇u 2 L ∞ + ∇r 2 L 4
dτ ,

where we have used the fact that δu in = 0. At this point, we observe that u 2 and ∇u 2 both belong to L 2 T (L ∞ ). Moreover, by (2.30) and the fact that

u 2 ∈ L ∞ T (L 2 ) and ∇u 2 ∈ L 4 T (L 2 ), one gathers that u 2 ∈ L 16/3 T (L 4 ). Finally, ∂ t u 2 ∈ L 4/3
T (H 1 ). In the end, similarly to what is done in (2.49), we deduce the existence of a positive exponent, that we keep calling α without loss of generality, such that

δu L ∞ t (L 2 ) ≤ C t α θ(t) + δu L ∞ t (L 2 )
(2.51) for all t ∈ 0, min{T, 1} . As a last step, we rewrite equation (2.50) in the following form:

(2.52)

ρ∂ t δu + Lδu = -r 1 ∂ t δu + δr ∂ t u 2 + δf ,
where we have defined

δf := ρ 1 u 1 • ∇δu + ρ 1 δu • ∇u 2 + δr u 2 • ∇u 2 + P ρ 1 ∇δr + P (ζ) δr ∇(ρ + r 2 ) -δr∇G .
Applying Proposition 2.2 to equation (2.52), with the choice (2.36) of the parameters, and using the smallness condition (2.1), we get, for any t ∈ [0, T ], the inequality

δu L 2 t (L ∞ ) + ∇δu L 4 t (L 2 ) + ∂ t δu, ∇ 2 δu L 4/3 t (L 2 ) ≤ C δu in Ḃ1/2 2,4/3 + δr ∂ t u 2 L 4/3 t (L 2 ) + δf L 4/3 t (L 2 ) + δu L 4/3 T (L 2 ) + ∇δu L 4/3 T (L 2 ) .
(2.53)

Recall that ∂ t u 2 ∈ L 2 T (L 2 )
δr ∂ t u 2 L 4/3 t (L 2 ) ≤ t 3/4-11/16 δr L ∞ t (L 4 ) ∂ t u 2 L 16/11 t (L 4 ) ≤ t 1/16 δr L ∞ t (H 1 ) ∂ t u 2 1/4 L 2 T (L 2 ) ∂ t ∇u 2 3/4 L 4/3 T (L 2 ) ≤ C t 1/16 θ(t) .
Next, we are going to bound δf in L 4/3 t (L 2 ). This was already done for the energy estimate (2.51) above, but here we have to take special care of the integrability in time. For the pressure terms and the potential force term we have

P (ζ) δr ∇(ρ + r 2 ) L 4/3 t (L 2 ) ≤ C t 3/4 δr L ∞ t (L 2 ) ∇ρ L ∞ T (L ∞ ) + δr L ∞ t (L 4 ) ∇r 2 L ∞ T (L 4 ) ≤ C t 3/4 θ(t) P ρ 1 ∇δr L 4/3 t (L 2 ) ≤ C t 3/4 ∇δr L ∞ t (L 2 ) ≤ C t 3/4 θ(t) δr∇G L 4/3
t (L 2 ) ≤ Ct 3/4 θ(t) . Thus, ut remains us to bound the convective terms. First of all, recall that we showed above that u 2 L 4 ∇u 2 L ∞ belongs to L 16/11 T

. Therefore,

δr u 2 • ∇u 2 L 4/3 t (L 2 ) ≤ t (3/4)-(11/16) δr L ∞ t (L 4 ) u 2 • ∇u 2 L 16/11 T (L 4 ) ≤ Ct 1/16 θ(t) .
In addition, we have

ρ 1 δu • ∇u 2 L 4/3 t (L 2 ) ≤ t (3/4)-(1/2) ρ 1 L ∞ T (L ∞ ) δu L ∞ t (L 2 ) ∇u 2 L 2 T (L ∞ ) ≤ C t 1/4 δu L ∞ t (L 2 )
. Finally, for the last term we use the following Gagliardo-Nirenberg type inequality (see e.g. Lemma II.3.3 of [START_REF] Galdi | An introduction to the mathematical theory of the Navier-Stokes equations. Steady-state problems[END_REF]):

u L ∞ ≤ u 1/4 L 2 ∇ 2 u 3/4 L 2 . This, combined with Young's inequality, implies u 2 ∈ L 4
T (L ∞ ) (recall the definition of the space X T ). Using this bound, we can estimate

ρ 1 u 1 • ∇δu L 4/3 t (L 2 ) ≤ t (3/4)-(1/2) ρ 1 L ∞ T (L ∞ ) u 1 L 4 T (L ∞ ) ∇δu L 4 t (L 2 ) ≤ C t 1/4 ∇δu L 4 t (L 2 )
. Putting all the previous estimates into (2.53), we have shown that there exists a positive exponent α > 0 for which, for all t ∈ 0, min{1, T } , one has

δu L 2 t (L ∞ ) + ∇δu L 4 t (L 2 ) + ∂ t δu, ∇ 2 δu L 4/3 t (L 2 )
(2.54)

≤ C t α θ(t) + δu L ∞ t (L 2 ) + ∇δu L 4 t (L 2 ) .
Let us now introduce the quantity

D(t) := δu L ∞ t (L 2 ) + δu L 2 t (L ∞ ) + ∇δu L 4 t (L 2 ) + ∂ t δu, ∇ 2 δu L 4/3 t (L 2
) . Summing up inequalities (2.49), (2.51) and (2.54), we finally deduce that, for all t ∈ 0, min{1, T } , we have

θ(t) + D(t) ≤ C t α θ(t) + D(t) ,
for a suitable exponent α > 0. Therefore, if t is now small enough, we can absorbe the right hand side into the left-hand side, deducing that both θ(t) and D(t) have to be 0. In particular, we deduce that ρ 1 ≡ ρ 2 and u 1 ≡ u 2 almost everywhere on [0, t] × Ω. In this way, we also infer that

∇δu L 2 t (L ∞ ) + ∇ 3 δu , ∂ t ∇δu L 4/3 t (L 2 ) = 0 , from which we deduce that δu L ∞ t (H 1 ) = 0
, where we have used also (2.41) and the estimate in Lemma 2.5.

To complete the argument, let us define the set

I := t ∈ [0, T ] ρ 1 (τ ) -ρ 2 (τ ) H 1 + u 1 (τ ) -u 2 (τ ) H 1 = 0, ∀ τ ∈ [0, t] .
Of course, I = ∅, since 0 ∈ I. In addition, the previous argument, combined with Lemma 2.5, shows that I is open. On the other hand, by continuity in time of the norms appearing in the definition of the space X T (see again Lemma 2.5 above), we infer that I is also closed. Then, by connectedness, we must have I = [0, T ]. This completes the proof of the proposition.

Quantitative asymptotic analysis with stratification effects and anisotropic diffusion

The goal of this section is to analyze the structure of the solutions to the highly rotating compressible system (1.3) with vertical stratification. The main result of this part is contained in Theorem 2, see Subsection 3.3: there we derive an asymptotic expansion and quantify the error in terms of the parameter ε. This stability result relies on relative entropy estimates.

3.1. Formal asymptotic expansion. In this section we perform formal computations in order to have a grasp on the structure of the solutions to sytem (1.3). Because of the no-slip boundary conditions (1.4), boundary layers appear in the limit ε → 0 both in the vicinity of the top boundary R 2 × {1} and of the bottom boundary R 2 × {0}.

We will specify later on the precise hypotheses on the initial conditions

ρ |t=0 = ρ in and u |t=0 = u in at time t = 0.
For the purpose of the formal analysis, let us say in a loose way that we impose the following far field conditions, for |x| → ∞:

ρ in (x) -→ ρ(x 3 ) and u in (x) -→ 0 ,
where ρ is a strictly positive function satisfying the logistic equation

(3.1) ∇P (ρ) = ρ ∇G .
In addition, the initial densities are assumed to be far away from vacuum. Moreover, we focus on well-prepared initial data, in the sense specified in Subsection 3.3 below.

3.1.1. Construction of the ansatz. We start with expanding the solution (u ε , ρ ε ) to (1.3) as

u ε = u 0 (x h , x 3 , t) + u bl 0,b (x h , x 3 ε , t) + u bl 0,t (x h , 1-x 3 ε , t) + ε u 1 (x h , x 3 , t) + u bl 1,b (x h , x 3 ε , t) + u bl 1,t (x h , 1-x 3 ε , t) + O(ε 2 ) ρ ε = ρ 0 (x h , x 3 , t) + ρ bl 0,b (x h , x 3 ε , t) + ρ bl 0,t (x h , 1-x 3 ε , t) + ε ρ 1 (x h , x 3 , t) + ρ bl 1,b (x h , x 3 ε , t) + ρ bl 1,t (x h , 1-x 3 ε , t) + ε 2 ρ 2 (x h , x 3 , t) + ρ bl 2,b (x h , x 3 ε , t) + ρ bl 2,t (x h , 1-x 3 ε , t) + O(ε 3 ) . (3.2)
The superscript bl stands for "boundary layer", while the subscripts b and t stand for "bottom" and "top" respectively. For simplicity of the presentation, in the next computations we are going to consider only the boundary layer near the bottom, since the terms related to the top boundary layer are dealt with in the exact same way. Therefore, from now on we omit the subscript b for the boundary layer terms. However, when needed, we will explicitly write t or b subscripts to avoid confusion. Below we denote by ζ = x 3 ε the fast vertical variable in the boundary layer. The boundary layer profiles are supposed to decay to 0 at exponential rate when ζ → ∞, since their effect is almost negligible in the interior of the domain: we will use this fact repeatedly in the following computations. We remark that, at this level, (3.2) is just a formal ansatz. As is usual, we will first formally derive the equations for the profiles: this is the purpose of the present section. After that, we will prove quantitative estimates for the difference between the solution and the profiles we have constructed, using the relative entropy method: this will be done in Section 3.

Identification of the profiles. In order to identify the profiles, we plug the ansatz (3.2) into (1.3) and identify the terms of the same order of magnitude in ε. We immediately notice that the highest order term is a term of order ε -3 , which appears in the third component of the momentum equation:

P (ρ 0 + ρ bl 0 ) ∂ ζ ρ bl 0 = 0 .
We assume that ρ 0 + ρ bl 0 stays bounded away from zero. This hypothesis is fully justified here below. In view of the hypothesis (1.10) on the pressure and the fact that ρ bl 0 has to vanish for ζ → ∞, we immediately deduce that ρ bl 0 ≡ 0. Thanks to that property, and ignoring the terms of order O(ε 2 ), which have been neglected in (3.2) in the expansion of the velocity fields, we find the following cascade of equations: from the conservation of mass equation, we get

ρ 0 ∂ ζ u bl 0,3 = 0 (mass-ε -1 ) ∂ t ρ 0 + ∇ h • ρ 0 (u 0,h + u bl 0,h ) + ∂ 3 (ρ 0 u 0,3 ) + ∂ 3 ρ 0 u bl 0,3 (mass-ε 0 ) + ρ 1 ∂ ζ u bl 0,3 + ∂ ζ (ρ bl 1 u bl 0,3 ) + ∂ ζ ρ bl 1 u 0,3 + ρ 0 ∂ ζ u bl 1
,3 = 0 , and from the momentum equation we get ∇P (ρ 0 ) + 0

P (ρ 0 )∂ ζ ρ bl 1 = 0 λ∂ 2 ζ u bl 0,3 + ρ 0 ∇G (mom-ε -2 ) ρ 0 (u 0,3 + u bl 0,3 ) • ∂ ζ u bl 0 + e 3 × ρ 0 (u 0 + u bl 0 )+ (mom-ε -1 ) ∇ h (P (ρ 0 )(ρ 1 + ρ bl 1 )) ∂ 3 (P (ρ 0 )ρ 1 ) + ∂ 3 (P (ρ 0 ))ρ bl 1 + P (ρ 0 ) ρ bl 1 ∂ ζ ρ bl 1 + P (ρ 0 ) ∂ ζ ρ bl 2 = ∂ 2 ζ u bl 0,h 0 + λ ∇ h ∂ ζ u bl 0,3 ∂ ζ ∇ h • u bl 0,h + ∂ 2 ζ u bl 1,3 + (ρ 1 + ρ bl 1 )∇G.
We will examine the equation at order O(ε 0 ) coming from the momentum equation later. Let us first infer some properties for the profiles.

The terms in the interior. Recall that the boundary layer profiles are expected to go to zero when ζ → ∞. Therefore, it follows from (mom-ε -2 ) that

(3.3) ∇P (ρ 0 ) = ρ 0 ∇G ,
which yields, by using (1.10), the properties

(3.4) H (ρ 0 ) = G + c(t) and ∇ h ρ 0 = 0 .
Hence ρ 0 is independent of x h , namely ρ 0 = ρ 0 (x 3 , t), and satisfies the ODE (3.5) P (ρ 0 )∂ 3 ρ 0 = -ρ 0 .

Since P ∈ C 1 ((0, ∞)) and non-zero, we can use Cauchy-Lipschitz theorem to get that ρ 0 (t) ∈ C 1 ((0, 1)), and hence bounded. Moreover, from (mom-ε -1 ) we infer

ρ 0 u ⊥ 0,h 0 + P (ρ 0 )∇ h ρ 1 ∂ 3 (P (ρ 0 )ρ 1 ) = ρ 1 ∇G .
This equation is called geostrophic balance; it implies the Taylor-Proudman theorem (see Section 1). In particular, its third component reads

∂ 3 (P (ρ 0 )ρ 1 ) = -ρ 1 .
Using the previous relation together with (3.3), we get

(3.6) ∂ 3 P (ρ 0 ) ρ 0 ρ 1 = 0 , hence the quantity (3.7) Q := P (ρ 0 ) ρ 0 ρ 1 verifies Q = Q(t, x h ) ,
i.e. Q is independent of the vertical variable. From the horizontal component, instead, we get (recall that ∇ h ρ 0 = 0)

(3.8) u 0,h = ∇ ⊥ h P (ρ 0 ) ρ 0 ρ 1 = ∇ ⊥ h Q .
In particular, we deduce that u 0,h = u 0,h (x h , t), which justifies the introduction of boundary layer terms in order to enforce the no-slip boundary conditions on x 3 = 0, 1.

In addition, applying the horizontal divergence we obtain

∇ h • u 0,h = ∇ h • ∇ ⊥ h P (ρ 0 ) ρ 0 ρ 1 = 0 ,
so that u 0,h is a 2-D horizontal divergence-free vector field.

We now exploit (mass-ε 0 ): considering it in the interior of the domain (i.e., neglecting the boundary terms) and using the inequalities just proved, after an integration in the vertical variable we infer that (3.9)

ˆ1 0 ∂ t ρ 0 dx 3 = - ˆ1 0 ∂ 3 (ρ 0 u 0,3 ) dx 3 = 0 .
By taking the time derivative of (3.4) and using (1.10) we have

∂ t ρ 0 = ∂ t c H (ρ 0 )
We integrate in the vertical variable and, from (3.9), we get ∂ t c = 0, hence ∂ t ρ 0 = 0. This implies that ρ 0 has to be independent also of time, and hence it is equal to a positive function ρ(x 3 ), solution of (3.3), or equivalently (3.1). Thanks to this fact, we have now that ∂ 3 (ρ 0 u 0,3 ) = 0. Using the no-slip boundary condition and the positivity of ρ 0 = ρ, we find u 0,3 ≡ 0. From now on ρ denotes ρ 0 . Let us now consider the equations outside the boundary layers, at order O(ε) in the mass equation, (mass-ε 1 )

∂ t ρ 1 + ∇ • (ρ u 1 ) + ∇ h • (ρ 1 u 0,h ) = 0 ,
and at order O(ε 0 ) in the momentum equation,

ρ ∂ t u 0 + ∇ • (ρ u 0 ⊗ u 0 ) + e 3 × (ρ 1 u 0 + ρ u 1 ) (mom-ε 0 ) + ∇ P (ρ) 2 ρ 2 1 + P (ρ) ρ 2 = µ∆ h u 0 + λ∇(∇ • u 0 ) + ρ 2 ∇G.
Recall that u 0 = (u 0,h (t, x h ), 0). Taking the curl of the horizontal component in (mom-ε 0 ), we obtain an equation for the horizontal vorticity

ω 0 = ∇ ⊥ h • u 0,h : ρ ∂ t ω 0 + ρ u 0,h • ∇ h ω 0 + ∇ h • (ρ 1 u 0,h ) + ∇ h • (ρ u 1,h ) -µ∆ h ω 0 = 0 .
Notice that, by (3.8), we get

ω 0 = ω 0 (t, x h ) = ∆ h Q ,
where Q is defined in (3.7); from the previous relation it follows that

(3.10) ρ ∂ t ∆ h Q + ρ ∇ ⊥ h Q • ∇ h ∆ h Q + ∇ h • (ρ u 1,h ) -µ∆ 2 h Q = 0
, where we have used the cancellation

(3.11) ∇ h • (ρ 1 ∇ ⊥ h ρ 1 ) = 1 2 ∇ h • ∇ ⊥ h (ρ 2 1 ) = 0
in order to get rid of the term ∇ h • (ρ 1 u 0,h ). In order to compute the term ∇ h • (ρ u 1,h ) in (3.10), we use equation (mass-ε 1 ) and the cancellation (3.11) again: we find (3.12)

∇ h • (ρ u 1,h ) = -∂ t ρ 1 -∇ h • (ρ 1 u 0,h ) -∂ 3 (ρ u 1,3 ) = -∂ t ρ 1 -∂ 3 (ρ u 1,3 ) .
After integrating in x 3 both (3.10) and (3.12) and summing up the resulting expressions, we eventually obtain

(3.13) ∂ t ρ ∆ h Q - ρ P (ρ) Q + ρ ∇ ⊥ h Q • ∇ h ∆ h Q -µ∆ 2 h Q = ρ(1) u 1,3 (x h , 1, t) -ρ(0) u 1,3 (x h , 0, t) ,
where f = ´1 0 f (x 3 ) dx 3 denotes the vertical mean of f . Boundary layer terms. We now consider the boundary layer terms. These terms are crucial to compute the right hand side of (3.13): indeed u bl j,3,b (x h , 0, t) = -u j,3 (x h , 0, t) and u bl j,3,t (x h , 0, t) = -u j,3 (x h , 1, t) (3.14) for j = 0, 1, in order to enforce the no-slip boundary condition on the bottom and top boundaries. First of all, (mass-ε -1 ) yields u bl 0,3 = u bl 0,3 (x h , t), and hence (3.15) u bl 0,3 ≡ 0 . Using (3.15), we obtain from (mom-ε -2 ) that

P (ρ) ∂ ζ ρ bl 1 = λ ∂ 2 ζ u bl 0,3 = 0 .
Hence, thanks to (1.10), ρ bl 1 = ρ bl 1 (x h , t) is constant in the boundary layer and goes to zero when ζ → ∞, therefore ρ bl 1 ≡ 0. Taking into account this last equality and reading the horizontal component of (mom-ε -1 ), one has [START_REF] Danchin | Paicu: A well-posedness result for viscous compressible fluids with only bounded density[END_REF]), x h is a parameter. We use Taylor formula at first order

(3.16) ρ (u bl 0,h ) ⊥ = ∂ 2 ζ u bl 0,h . Notice that, in (3.
ρ(x 3 ) = ρ(0) + x 3 ˆ1 0 ∂ 3 ρ(s x 3 ) ds to write (3.16) as (3.17) ρ(0) (u bl 0,h ) ⊥ + x 3 ˆ1 0 ∂ 3 ρ (s x 3 ) ds (u bl 0,h ) ⊥ = ∂ 2 ζ u bl 0,h .
Let us now consider the equation

(3.18) ρ(0) (u bl 0,h ) ⊥ = ∂ 2 ζ u bl 0,h , supplemented with the boundary condition (3.19)
u bl 0,h (x h , 0, t) = -u 0,h (x h , t) at ζ = 0, in view of (1.4) and (3.14). We remark that the system of ODEs (3.18)-(3.19) is the same (here in general ρ(0) = 1) as in the incompressible case, see e.g. Chapter 7 of [START_REF] Chemin | Mathematical geophysics. An introduction to rotating fluids and the Navier-Stokes equations[END_REF] and references therein. Its solutions are exponentially decaying and have a spiral structure. Indeed, we have the following formula:

u bl 0,h,b (x h , ζ, t) = -     e -ζ ρ(0) 2 u 0,1 (x h , t) cos ζ ρ(0) 2 + u 0,2 (x h , t) sin ζ ρ(0) 2 e -ζ ρ(0) 2 -u 0,1 (x h , t) sin ζ ρ(0) 2 + u 0,2 (x h , t) cos ζ ρ(0) 2     .
Let us move further. The vertical component in (mom-ε -1 ) is

(3.20) 0 = λ ∂ ζ ∇ h • u bl 0,h + ∂ 2 ζ u bl 1,3 + P (ρ) ∂ ζ ρ bl 2 .
Equation (mass-ε 0 ), together with the fact the ρ is strictly positive, yields

(3.21) ∇ h • u bl 0,h + ∂ ζ u bl 1,3 = 0 . Hence P (ρ) ∂ ζ ρ bl 2 =
0 and, similarly to the argument used for ρ bl 1 , we get ρ bl 2 ≡ 0. The previous equality (3.21) determines u bl 1,3 up to a constant in ζ, which we take so that u bl 1,3 converges to zero when ζ → ∞:

u bl 1,3,b (x h , ζ, t) = - e -ζ ρ(0) 2 2ρ(0) cos ζ ρ(0) 2 + sin ζ ρ(0) 2 ∇ ⊥ h • u 0,h (x h , t) .
Similar computations can be done for the top boundary layers. Indeed, denoting by η = 1-x 3 ε the fast vertical variable in the upper boundary layer, we use Taylor formula at first order (3.22) ρ(x 3 ) = ρ( 1)

-(1 -x 3 ) ˆ1 0 ∂ 3 ρ(1 -s(1 -x 3 ))ds
to define u bl 0,h,t as the solution to the equation

(3.23) ρ(1) (u bl 0,h,t ) ⊥ = ∂ 2 η u bl 0,h,t , supplemented with the boundary condition (3.24) u bl 0,h,t (x h , 0, t) = -u 0,h (x h , t)
at η = 0, recall (3.14). We have

u bl 0,h,t (x h , η, t) = -     e -η ρ(1) 2 u 0,1 (x h , t) cos η ρ(1) 2 + u 0,2 (x h , t) sin η ρ(1) 2 e -η ρ(1) 2 -u 0,1 (x h , t) sin η ρ(1) 2 + u 0,2 (x h , t) cos η ρ(1) 2     . and, from (3.21) with ∂ ζ replaced by -∂ η , u bl 1,3,t (x h , η, t) = e -η ρ(1) 2 2ρ(1) cos η ρ(1) 2 + sin η ρ(1) 2 ∇ ⊥ h • u 0,h (x h , t) .
Hence, using (3.14), one can now compute the right hand side of equation (3.13):

(3.25) ρ(1) u 1,3 (x h , 1, t) -ρ(0) u 1,3 (x h , 0, t) = -ρ(1)u bl 1,3,t (x h , 0, t) + ρ(0)u bl 1,3,b (x h , 0, t) = - √ ρ(0)+ √ ρ(1) √ 2 ω 0 = - √ ρ(0)+ √ ρ(1) √ 2 ∆ h Q .
This is the so-called Ekman pumping term, which represents the secondary (global) circulation created by the boundary layer. It appears as a damping term for the quasi-geostrophic dynamics, described by equation (3.13).

Final choices for correctors. It remains to choose the functions ρ 2 , u 1 and u bl 1,h . These terms are auxiliary terms which do not appear in the final result. We choose the interior terms in order to make the terms of order O(ε) in the mass equation and the terms of order O(ε 0 ) in the momentum equation vanish identically. Notice that (mom-ε 0 ) determines u 1,h in terms of u 0 , ρ 1 and ρ 2 , and hence, through relation (3.8), in terms of ρ 1 and ρ 2 only. Specifically,

u 1,h := 1 ρ -µ∆ h u ⊥ 0,h + ρ∂ t u ⊥ 0,h (3.26) + ρu 0,h • ∇ h u ⊥ 0,h -u 0,h ρ 1 + ∇ ⊥ h P (ρ)ρ 2 + P (ρ) 2 ρ 2 1 .
Next, the vertical component of (mom-ε 0 ) reads (3.27)

∂ 3 P (ρ) ρ 2 + P (ρ) 2 ρ 2 1 = -ρ 2
where we have used that u 0,3 ≡ 0. Since, by (1.10), P (ρ) > 0, ρ 2 can be defined as the solution of the ODE (3.28)

∂ 3 ρ 2 + ∂ 3 (P (ρ)) + 1 P (ρ) ρ 2 = - ∂ 3 P (ρ) ρ 2 1 2 P (ρ) ,
up to an arbitrary constant c(x h , t), that we take equal to zero for simplicity. We remark that this choice does not affect the choice of the other quantities since ρ 2 appears only in (mom-ε 0 ) or higher order equations. Moreover, since ρ 1 and ∇ρ 1 are bounded in time and space (Q, defined in (3.7), satisfies the quasi-geostrophic equation (3.13), which admits regular solutions, see Lemma 3.3 later), ρ 2 and ∇ρ 2 are also bounded in time and space.

Next, equation (mass-ε 1 ) determines u 1,3 up to a constant in x 3 , which we take equal to -u bl 1,3 (x h , 0, t) in order to enforce the no-slip boundary condition for the vertical component at order O(ε). Therefore, thanks to (3.12) we get

ρ(x 3 ) u 1,3 (x h , x 3 , t) = -ρ(0) u bl 1,3 (x h , 0, t) - ˆx3 0 (∂ t ρ 1 + ρ ∇ h • u 1,h ) dz . (3.29)
Differently from the case without the gravitational potential, this term does not have an affine structure as in the incompressible case (see again Chapter 7 of [START_REF] Chemin | Mathematical geophysics. An introduction to rotating fluids and the Navier-Stokes equations[END_REF]), since u 1,h does not depend only on x h . In order to enforce the no-slip boundary condition at order O(ε) also for the horizontal component, we impose

u bl 1,h (x h , 0, t) = -u 1,h (x h , 0, t) (3.
u bl 1,h,b (x h , ζ, t) := -u 1,h (x h , 0, t)e -ζ ρ(0) 2 .
Analogously, u bl 1,h,t is defined for all η ∈ [0, ∞) and

x h ∈ R 2 by u bl 1,h,t (x h , η, t) := -u 1,h (x h , 1, t)e -η ρ(1) 2 .
Remark 3.1. Contrary to the interior terms, it is not possible to make the terms of order O(ε) in the mass equation and the terms of order O(ε 0 ) in the momentum equation vanish identically. Indeed that would come down to imposing

ρ ∇ h • u bl 1,h = -∇ h ρ 1 • u bl 0,h -∂ 3 ρ u bl 1,3 λ∂ ζ (∇ h • u bl 1,h ) = -∂ 2 ζ u bl 1,3 = ∂ ζ (∇ h • u bl 0,h ) , (3.31)
which is overdetermined. This fact is due to the lack of higher-order correctors, since we truncate the expansion at order one in ε.

Notice that, due to exponential decay to zero in the interior of the domain, the boundary layer terms are small. Moreover, we can exploit their decay by relying on Hardy's inequality (see the computations in Section 3). The final stability estimate, though, will be worse than in the absence of boundary layer phenomena (as e.g. for complete slip boundary conditions). Improving this estimate would require to consider higher-order correctors in the ansatz (3.2). Notice also that, using (3.14), we have at the bottom

x 3 = 0 (3.32) u 0,h (x h , t) + u bl 0,h,b (x h , 0, t) + u bl 0,h,t (x h , 1 ε , t) = u bl 0,h,t (x h , 1 ε , t) u 1 (x h , 0, t) + u bl 1,b (x h , 0, t) + u bl 1,t (x h , 1 ε , t) = u bl 1,t (x h , 1 ε , t) ,
and at the top

x 3 = 1 (3.33) u 0,h (x h , t) + u bl 0,h,b (x h , 1 ε , t) + u bl 0,h,t (x h , 0, t) = u bl 0,h,b (x h , 1 ε , t) u 1 (x h , 1, t) + u bl 1,b (x h , 1 ε , t) + u bl 1,t (x h , 0, t) = u bl 1,b (x h , 1 ε , t).
It means that we have a (exponentially small, but still non-zero) trace of the top boundary layer on the bottom boundary and vice-versa. Hence, we will add corrector terms in the ansatz (3.34) below, in order to keep homogeneous boundary conditions. This is a technical point, but needed to apply Hardy's inequality later.

The ansatz. To put it in a nutshell, we have obtained the following ansatz for the structure of the solutions to (1.3)-(1.4):

ρ ε app (x h , x 3 , t) = ρ(x 3 ) + ε ρ 1 (x h , x 3 , t) + ε 2 ρ 2 (x h , x 3 , t) u ε app (x h , x 3 , ζ, η, t) = ∇ ⊥ h Q(x h , t) + u bl 0,h,b (x h , ζ, t) + u bl 0,h,t (x h , η, t) -u bl 0,h,1/ε (x h , x 3 , t) 0 + ε u 1,h (x h , x 3 , t) + u bl 1,h,b (x h , ζ, t) + u bl 1,h,t (x h , η, t) -u bl 1,h,1/ε (x h , x 3 , t) u 1,3 (x h , x 3 , t) + u bl 1,3,b (x h , ζ, t) + u bl 1,3,t (x h , η, t) -u bl 1,3,1/ε (x h , x 3 , t) (3.34)
with Q defined in (3.7) and

(3.35)

u bl 0,h,1/ε (x h , x 3 , t) = x 3 u bl 0,h,b (x h , 1 ε , t) + (1 -x 3 ) u bl 0,h,t (x h , 1 ε , t) u bl 1,1/ε (x h , x 3 , t) = x 3 u bl 1,b (x h , 1 ε , t) + (1 -x 3 ) u bl 1,t (x h , 1 ε , t) .
In addition, it follows from (3.13) that (3.36)

∂ t ρ P (ρ) Q -ρ ∆ h Q -ρ ∇ ⊥ h Q • ∇ h ∆ h Q + µ∆ 2 h Q - √ ρ(0)+ √ ρ(1) √ 2 ∆ h Q = 0 .
This is the quasi-geostrophic equation. Similar limit equations without damping term have been shown in e.g. [START_REF] De Anna | Global well-posedness and long-time dynamics for a higher order Quasi-Geostrophic type equation[END_REF], [START_REF] Feireisl | A singular limit for compressible rotating fluids[END_REF] and [START_REF] Feireisl | Rotating compressible fluids under strong stratification[END_REF], where the boundary layers do not appear due to the complete slip condition. Notice that in [START_REF] Feireisl | Rotating compressible fluids under strong stratification[END_REF] the parabolic term disappears, since the authors consider also the inviscid limit. We state here the well-posedness and the regularity results for the quasi-geostrophic equation (3.36), whose detailed proofs are given in [START_REF] Bocchi | Compressible-incompressible transitions in fluid mechanics: waves-structures inter-action and rotating fluids[END_REF].

Proposition 3.2. Let Q in ∈ H 1 (R 2 ).
Then, there exists a unique global weak solution Q to the quasi-geostrophic equation (3.36) with initial datum Q in , such that

Q ∈ C(R + ; H 1 (R 2 )) ∩ L ∞ (R + ; H 1 (R 2 )) and ∇ h Q ∈ L 2 (R + ; H 1 (R 2 )) .
Lemma 3.3. Let n ≥ 1 be an integer and Q in ∈ H n (R 2 ). Then, there exists a constant C n-1 > 0 such that any weak solution to (3.36) with initial datum Q in satisfies the following inequality for all t ≥ 0:

(3.37)

n-1 j=0 ∇ j h Q(t) 2 L 2 + ∇ j+1 h Q(t) 2 L 2 + n-1 j=0 ˆt 0 ∇ j+1 h Q 2 L 2 + ∇ j+2 h Q 2 L 2 ≤ C n-1 n-1 j=0 ∇ j h Q in L 2 + ∇ j+1 h Q in 2 L 2 ,
where

C 0 = C 1 = 1 and C n-1 = C n-1 ( Q in H n-1 ) for n -1 ≥ 2.
The boundary layer profiles u bl 0,h,b and u bl 0,h,t are solutions of the systems (3.18) -(3.19) and (3.23) - (3.24) respectively. We refer to the previous computations for the precise definitions of the higher-order terms. We conclude this part by remarking that, according to the previous computations, we have that (ρ ε app , u ε app ) solves the following system:

(3.38)                    ∂ t ρ ε app + ∇ • (ρ ε app u ε app ) = εR bl + ε 2 R ε ρ ε app ∂ t u ε app + ρ ε app u ε app • ∇u ε app + 1 ε e 3 × ρ ε app u ε app + 1 ε 2 ∇P (ρ ε app ) = 1 ε 2 ρ ε app ∇G + x 3 ε ˆ1 0 ∂ 3 ρ (s x 3 ) ds e 3 × u bl 0,h,b - 1 -x 3 ε ˆ1 0 ∂ 3 ρ (1 -s(1 -x 3 )) ds e 3 × u bl 0,h,t + ∆ µ,ε u ε app + λ∇(∇ • u ε app ) + S bl +
εS ε in the slab Ω with no-slip boundary conditions (1.4). The remainder terms R ε and S ε are of the form

R ε = R ε (x h , x 3 , x 3 ε , 1-x 3 ε , t) and S ε = S ε (x h , x 3 , x 3 ε , 1-x 3
ε , t) while the boundary layer terms are

R bl (x h , x 3 , x 3 ε , 1-x 3 ε , t) = ρ ∇ h • (u bl 1,h,b + u bl 1,h,t ) + ∇ h ρ 1 • (u bl 0,h,b + u bl 0,h,t ) + ∂ 3 ρ (u bl 1,3,b + u bl 1,3,t ) S bl (x h , x 3 , x 3 ε , 1-x 3 ε , t) = ρ ∂ t (u bl 0,b + u bl 0,t ) + ρ u 0,h • ∇ h (u bl 0,b + u bl 0,t ) + ρ (u bl 0,b + u bl 0,t ) • ∇ h u 0 + ρ (u bl 0,b + u bl 0,t ) • ∇ h (u bl 0,b + u bl 0,t ) + ρ (u 1,3 + u bl 1,3,b + u bl 1,3,t ) ∂ η (u bl 0,b + u bl 0,t ) -µ∆ h (u bl 0,t + u bl 0,b ) -∂ 2 η (u bl 1,b + u bl 1,t ) -λ 0 ∂ η ∇ h • (u bl 1,h,b + u bl 1,h,t ) + e 3 × ρ 1 (u bl 0,b + u bl 0,t ) + ρ(u bl 1,b + u bl 1,t ) .
The remainders ε 2 R ε and εS ε contain also the terms of order O(e -1/ε ) coming from the correctors u bl 0,h,1/ε and u bl 1,1/ε , defined in (3.35). Notice that S bl appears at order O(1), but has fast, exponential, decay inside Ω: more precisely, we have

S bl L p ≤ C ε 1 p for all p ∈ [1, ∞].
The choice of the regularity of the initial datum Q in guarantees enough regularity for the approximate solution (ρ ε app , u ε app ) in order to derive the stability estimates later in Subsection 3.3. This is stated in the following lemma, which is a straightforward consequence of Lemma 3.3 above. Lemma 3.4. The approximated density ρ ε app and velocity field u ε app can be written as

ρ ε app (x h , x 3 , t) = ρ(x 3 ) + Q(t, x h ) ρ P (ρ) (x 3 ) + Q(t, x h )l(x 3 ) u ε app (x h , x 3 , ζ, η, t) = N i=1 f i (t, x h )g i (x 3 )h i (ζ)w i (η) , for some N ≥ 1, with ρ, l, g i ∈ C 1 ([0, 1]) and h i , w i ∈ C ∞ (R + ). In addition, for Q in ∈ H 5 (R 2 ), we have Q ∈ L ∞ (R + ; H 5 (R 2 )), f i ∈ L ∞ (R + ; H k i (R 2 )) with k i ≥ 1.
3.1.2. Large-scale quasi-geostrophic equation. We recover here the equation for u 0 from (3.36). For this we need the following standard lemma, which gives the Helmholtz decomposition for two-dimensional vector fields.

Lemma 3.5. Let p ∈ (1, ∞). Let a and b be two scalar fields in L p (R 2 ).

Then there exists a unique vector field F , belonging to the homogeneous Sobolev space Ẇ 1,p (R 2 ; R 2 ), which solves the system

(3.39) ∇ ⊥ h • F = a ∇ h • F = b .
Moreover, the following formula holds:

F = -∇ ⊥ h (-∆ h ) -1 a -∇ h (-∆ h ) -1 b .
The previous result being classical, we do not give the proof here: we rather refer to [START_REF] Chemin | Fluides parfaits incompressibles[END_REF] (see Sections 1.2 and 1.3) and [START_REF] Feireisl | Singular limits in thermodynamics of viscous fluids[END_REF] (see Section 10.6) for details. We just give some explanations about the uniqueness, which will be needed below. By linearity, let us suppose that F solves (3.39) with a = b = 0. In particular ∇ × F = 0, hence (see Corollary 1.2.1 of [START_REF] Chemin | Fluides parfaits incompressibles[END_REF]) F = ∇h, for some h ∈ L p . But from ∇ • F = 0, we deduce that -∆h = 0, which admits the only solution h = 0 in L p Now, let π ∈ Ḣ1 (R 2 ) be the (unique, up to additive constants) solution to

(3.40) -∆ h π = ρ ∇ h • (u 0,h • ∇ h u 0,h ) = ρ ∇ h u 0,h : ∇ h u 0,h .
We then define F (•, t) ∈ L 2 (R 2 ; R 2 ) for almost every t ≥ 0 by the formula

F := ρ ∂ t u 0,h + ρ u 0,h • ∇ h u 0,h -µ∆ h u 0,h + √ ρ(0)+ √ ρ(1) √ 2 u 0,h + ∇ h π .
Notice that, thanks to equations (3.40) and (3.36) and the divergence-free condition ∇ h • u 0,h = 0, we have

∇ ⊥ h • F = ρ P (ρ) ∂ t Q and ∇ h • F = 0 .
Therefore, the uniqueness part of Lemma 3.5 implies that

F = ∇ ⊥ h (∆ h ) -1 ρ P (ρ) ∂ t Q = ρ P (ρ) ∂ t (∆ h ) -1 u 0,h
, where we have also used (3.8). Eventually, we find that u 0,h solves the system (3.41)

             ∂ t ρ - ρ P (ρ) (∆ h ) -1 u 0,h + ρ u 0,h • ∇ h u 0,h -µ∆ h u 0,h + √ ρ(0)+ √ ρ(1) √ 2 u 0,h + ∇ h π = 0 ∇ h • u 0,h = 0 in R 2 .
The second term appearing in the time derivative is a consequence of the combination of the effects due to density stratification and fast rotation. Notice that both (3.13) and (3.41) are averaged (in x 3 ) versions of (mom-ε 0 ).

3.2.

Weak solutions and uniform a priori bounds. We recall here some basics about finite energy weak solutions to system (1.3). We refer e.g. to [START_REF] Lions | Compressible models[END_REF], [START_REF] Novotny | Introduction to the mathematical theory of compressible flow[END_REF] and [START_REF] Feireisl | Singular limits in thermodynamics of viscous fluids[END_REF] for details.

Definition 3.6. Let ρ > 0 be the solution to the logistic equation (3.1), and let

(ρ in , u in ) verify ˆΩ 1 2 ρ in |u in | 2 + 1 ε 2 E (ρ in , ρ) dx < ∞ .
A couple (ρ, u) is a finite-energy weak solution to system (1.3) on [0, T ] × Ω, related to the initial datum (ρ in , u in ), if the following conditions are satisfied:

• ρ ≥ 0, with ρ -ρ ∈ L ∞ (0, T ); (L 2 + L γ )(Ω)
, with γ > 1 appearing in (1.10), and u ∈ L 2 (0, T ); H 1 (Ω; R 3 ) ; • the mass equation is satisfied in the weak sense: namely, for any test-function

ϕ ∈ C ∞ 0 [0, T ) × Ω , one has - ˆT 0 ˆΩ ρ ∂ t ϕ + ρ u • ∇ϕ dx dt = ˆΩ ρ in ϕ(0) dx ; • P (ρ) ∈ L 1 loc (0, T ) × Ω ,
and the momentum equation is verified in the weak sense: for any

ψ ∈ C ∞ 0 [0, T ) × Ω; R 3 , one has ˆT 0 ˆΩ -ρ u • ∂ t ψ -ρ u ⊗ u : ∇ψ + 1 ε e 3 × (ρ u) • ψ - 1 ε 2 P (ρ) ∇ • ψ + ∇ µ,ε u : ∇ µ,ε ψ + λ ∇ • u ∇ • ψ - 1 ε 2 ρ ∇G • ψ dx dt = ˆΩ ρ in u in • ψ(0) ;
• the following energy inequality holds true for almost every t ∈ (0, T ):

ˆΩ 1 2 ρ(t)|u(t)| 2 + 1 ε 2 E ρ(t), ρ dx (3.42) + ˆt 0 ˆΩ µ|∇ h u| 2 + ε|∂ 3 u| 2 + λ|∇ • u| 2 dx dτ ≤ ˆΩ 1 2 ρ in |u in | 2 + 1 ε 2 E (ρ in , ρ) dx ,
where we have defined the relative energy functional

(3.43) E ρ, ρ := H ρ -H ρ -H ρ ρ -ρ) .
The solution is said to be global if the previous conditions hold true for all T > 0.

Consider now a family of global in time finite-energy weak solutions (ρ ε , u ε ) ε to system (1.3). Recall that the existence of such a family is only assumed here, but open in general, because of anisotropy of the viscous stress tensor. We collect here some uniform bounds verified by that family. We refer e.g. to [START_REF] Feireisl | Singular limits in thermodynamics of viscous fluids[END_REF], [START_REF] Feireisl | A singular limit for compressible rotating fluids[END_REF] and [START_REF] Feireisl | Multi-scale analysis of compressible viscous and rotating fluids[END_REF] for details. These bounds will be important in the next subsection, when proving stability estimates.

By assumption, for any ε ∈ (0, 1] the energy inequality

ˆΩ ρ ε (t)|u ε (t)| 2 + 1 ε 2 E (ρ ε (t), ρ) (3.44) + ˆt 0 ˆΩ µ|∇ h u ε | 2 + ε|∂ 3 u ε | 2 + λ|∇ • u ε | 2 ≤ ˆΩ ρ ε in |u ε in | 2 + 2 ε 2 E (ρ ε in , ρ)
holds for almost every t > 0. According to inequality (4.15) of [START_REF] Feireisl | Sun: Suitable weak solutions to the Navier-Stokes equations of compressible viscous fluids[END_REF], we have the following control, which holds for any positive scalar functions ρ(x, t) and r(x, t), with 0 < r -≤ r(x, t) ≤ r + , for some real numbers r -, r + : there exist constants c 1 , c 2 > 0 such that, for almost all (x, t) ∈ Ω × R + , one has

c 1 |ρ(x, t) -r(x, t)| 2 1 {|ρ-r|(•,t)<1} + |ρ(x, t) -r(x, t)| γ 1 {|ρ-r|(•,t)≥1} ≤ E (ρ(x, t), r(x, t)) ≤ c 2 |ρ(x, t) -r(x, t)| 2 1 {|ρ-r|(•,t)<1} + |ρ(x, t) -r(x, t)| γ 1 {|ρ-r|(•,t)≥1} , (3.45)
where the notation {|ρ -r|(•, t) < 1} stands for the set of x ∈ Ω such that |ρ(x, t)r(x, t)| < 1 (and analogously for the ≥ symbol) and 1 A denotes the characteristic function of a set A ⊆ Ω. Notice that the same inequalities hold if we replace 1 by any constant M > 0, up to change the value of the constants c 1 and c 2 . Now, following Chapters 4 and 5 of [START_REF] Feireisl | Singular limits in thermodynamics of viscous fluids[END_REF], let us introduce the essential set and the residual set as follows: for almost every t > 0, we set

Ω ess (t) := x ∈ Ω |ρ ε (x, t) -ρ(x 3 )| < σ and Ω res (t) := Ω \ Ω ess (t) , (3.46)
for some σ (to be fixed later) such that 0 < σ < inf (0,1) ρ .

Accordingly, given any function h, we define its essential part and residual part as

[h] ess := h 1 Ωess and [h] res := h 1 Ωres = h -[h] ess .
Keep in mind that such a decomposition depends on ρ ε .

After this preparation, let us establish uniform bounds for (ρ ε , u ε ) ε . First of all, in view of the assumptions we will fix on the initial data (ρ ε in , u ε in ) ε in the next subsection, we can assume that the right hand side of (3.44) is uniformly bounded for ε ∈ (0, 1]. Then, using (3.45), we deduce the existence of a constant C > 0 such that, for all T > 0 fixed and all 0 < ε ≤ 1, one has

√ ρ ε u ε L ∞ T (L 2 ) ≤ C (3.47) 1 ε [ρ ε -ρ] ess L ∞ T (L 2 ) ≤ C (3.48) sup t∈[0,T ] L (Ω res (t)) + [ρ ε ] res γ L ∞ T (L γ ) ≤ C ε 2 , (3.49)
where L(A) denotes the Lebesgue measure of a set A ⊆ Ω. We refer to Section 2 of [START_REF] Feireisl | A singular limit for compressible rotating fluids[END_REF] and Section 4 of [START_REF] Feireisl | Multiple scales and singular limits for compressible rotating fluids with general initial data[END_REF] for details.

Next, let us consider the viscosity terms: recalling that µ > 0 and λ > 0, from (3.44) we immediately get

∇ h u ε L 2 T (L 2 ) + ∇ • u ε L 2 T (L 2 ) ≤ C (3.50) √ ε ∂ 3 u ε h L 2 T (L 2 ) ≤ C , (3.51)
for some universal constant C > 0 independent of ε and of the fixed time T > 0. In addition, owing to the identity

∂ 3 u ε 3 = ∇ • u ε -∇ h • u ε h ,
we also deduce that

(3.52) ∂ 3 u ε 3 L 2 T (L 2 ) ≤ C .
Finally, arguing exactly as in Section 2 of [START_REF] Feireisl | A singular limit for compressible rotating fluids[END_REF], we deduce that there exists a constant C > 0 such that, for all ε > 0 and all T > 0, one has

(3.53) u ε L 2 T (L 2 ) ≤ C .
3.3. Stability estimates. This section is devoted to estimating the error between weak solutions to (1.3) and their smooth approximation built in Subsection 3.1. We consider well-prepared initial data. Specifically, the initial density (ρ ε in ) ε and velocity fields (u ε in ) ε satisfy the following requirements:

• for all ε ∈ (0, 1], one has

(3.54) ρ ε in = ρ + ε r ε in , with (r ε in ) ε ⊆ L 2 ∩ L ∞ (Ω) ; • we have (u ε in ) ε ⊆ L 2 (Ω); • there exists Q in ∈ H 5 (R 2 ) such that, after defining (3.55) r in := ρ P (ρ) Q in and u in := (-∂ 2 Q in , ∂ 1 Q in , 0) ,
we have the strong convergence properties

r ε in -→ r in and u ε in -→ u in in L 2 (Ω) . (3.56) Remark 3.7. Condition (3.55) implies in particular that (3.57) ρ u ⊥ in,h 0 + P (ρ)∇ h r in ∂ 3 (P (ρ)r in ) = r in ∇G .
From the uniform bounds in Section 3.2, it is classical to derive that, up to extraction of a suitable subsequence, weak solutions (ρ ε , u ε ) ε converge to a limit state (ρ, u) which belongs to the kernel of the singular perturbation operator. We refer to e.g. [START_REF] Feireisl | A singular limit for compressible rotating fluids[END_REF], [START_REF] Feireisl | Multi-scale analysis of compressible viscous and rotating fluids[END_REF], [START_REF] Feireisl | Multiple scales and singular limits for compressible rotating fluids with general initial data[END_REF] for details. The goal of the present subsection is to make this convergence quantitative, to show the general structure of the solutions and to take into account the correctors due to Ekman's boundary layers. We aim at proving the following result. Recall that the relative entropy E is defined in (3.43).

Theorem 2. For γ ≥ 3/2, suppose that there exists a finite-energy weak solution (ρ ε , u ε ) ε to (1.3) with well-prepared initial data (ρ ε in , u ε in ) ε ∈ L ∞ × L 2 verifying hypotheses (3.54), (3.55) and (3.56). Let (ρ ε app , u ε app ) ε be defined as in (3.34), and define δu ε = u ε -u ε app . Then, there exist functions C 1 (t), C 2 (t) ∈ L 1 ([0, T )) for all T > 0, and constants C > 0 and ε 0 ∈ (0, 1) such that, for all ε ∈ (0, ε 0 ), the following estimate holds, for almost every t > 0:

ˆΩ ρ ε (t)|δu ε (t)| 2 dx + 1 ε 2 ˆΩ E(ρ ε (t), ρ ε app (t)) dx + ˆt 0 ˆΩ µ|∇ h δu ε | 2 + ε|∂ 3 δu ε | 2 + λ|∇ • δu ε | 2 dx ≤ Ce 2 ´t 0 C 1 (s)ds ˆΩ ρ ε in |δu ε in | 2 dx + 1 ε 2 ˆΩ E(ρ ε in , ρ ε in,app )dx + ε ˆt 0 C 2 (τ ) dτ .
Remark 3.8. The lower bound for the exponent γ comes from the control of the source term in the relative entropy inequality (3.64): in particular, in (3.83) we need γ ≥ 3/2 to apply Hölder's inequality and get the estimate (3.84).

In order to prove the previous result, we resort to the technique of the relative entropy/relative energy inequality, see e.g. [START_REF] Germain | Weak-strong uniqueness for the isentropic compressible Navier-Stokes system[END_REF], [START_REF] Feireisl | Sun: Suitable weak solutions to the Navier-Stokes equations of compressible viscous fluids[END_REF], [START_REF] Feireisl | Relative entropies, suitable weak solutions, and weak-strong uniqueness for the compressible Navier-Stokes system[END_REF], [START_REF] Feireisl | Multiple scales and singular limits for compressible rotating fluids with general initial data[END_REF] and [START_REF] Feireisl | Rotating compressible fluids under strong stratification[END_REF]. The relative entropy estimate of those works is directly applicable in our framework, but it is not immediately clear how to take advantage of the small remainders in (3.38). Instead, we directly derive the entropy inequality on the system for (δρ ε , δu ε ) and take into account from the beginning that (ρ ε app , u ε app ) is almost a solution to (1.3). On the contrary, the relative entropy inequality of e.g. [START_REF] Feireisl | Multiple scales and singular limits for compressible rotating fluids with general initial data[END_REF] holds for a much wider class of smooth functions.

3.3.1.

The relative entropy inequality. We set

δρ ε := ρ ε -ρ ε app and δu ε := u ε -u ε app .
From systems (1.3) and (3.38), it is easy to find an equation for δρ ε and δu ε : after setting δP ε := P (ρ ε ) -P (ρ ε app ), we get

∂ t δρ ε + ∇ • (u ε app δρ ε ) = -∇ • (ρ ε δu ε ) -εR bl -ε 2 R ε (3.58) ρ ε ∂ t δu ε + ρ ε u ε • ∇δu ε + 1 ε e 3 × ρ ε δu ε + 1 ε 2 ∇δP ε -∆ µ,ε δu ε -λ∇∇ • δu ε (3.59) = 1 ε 2 δρ ε ∇G -δρ ε ∂ t u ε app + (ρ ε app u ε app -ρ ε u ε ) • ∇u ε app - 1 ε e 3 × δρ ε u ε app -S bl -εS ε - x 3 ε ˆ1 0 ∂ 3 ρ(sx 3 )dse 3 × u bl 0,h,b + 1 -x 3 ε ˆ1 0 ∂ 3 ρ(1 -s(1 -x 3 ))dse 3 × u bl 0,h,t .
From the point of view of energy estimates, the main term to work on is the difference of the pressure terms. Testing it against δu ε yields

ˆΩ ∇δP ε • δu ε dx = ˆΩ ∇P (ρ ε ) • u ε dx -ˆΩ ∇P (ρ ε app ) • u ε app dx (3.60) + ˆΩ ∇ • u ε app δP ε dx -ˆΩ ∇P (ρ ε app ) • δu ε dx .
By standard computations, using the mass equation in (1.3), we get

ˆΩ ∇P (ρ ε ) • u ε dx = ˆΩ ∇ H (ρ ε ) • ρ ε u ε dx = d dt ˆΩ H(ρ ε )dx .
Similarly, from the first equation in (3.38) we gather

ˆΩ ∇P (ρ ε app ) • u ε app dx = d dt ˆΩ H(ρ ε app )dx -ε ˆΩ H (ρ ε app )(R bl + εR ε )dx .
In identity (3.60), we now add and substract the term d dt ´H (ρ ε app )δρ ε dx, in order to make the relative entropy E ρ ε (t), ρ ε app (t) appear. Then, from (3.60) and the previous computations we infer

ˆΩ ∇δP ε • δu ε dx = d dt ˆΩ E ρ ε , ρ ε app dx + ˆΩ ∇ • u ε app δP ε dx -ˆΩ ∇P (ρ ε app ) • δu ε dx + d dt ˆΩ H (ρ ε app ) δρ ε dx + ε ˆΩ H (ρ ε app ) R bl + εR ε dx .
Using again the mass equations in (1.3) and (3.38), we get

d dt ˆΩ H (ρ ε app ) δρ ε dx = ˆΩ ∂ t H (ρ ε app ) δρ ε dx + ˆΩ H (ρ ε app )∂ t δρ ε dx = ˆΩ ∂ t H (ρ ε app ) δρ ε dx + ˆΩ ∇H (ρ ε app ) • (ρ ε u ε -ρ ε app u ε app )dx -ε ˆΩ H (ρ ε app ) R bl + εR ε dx .
This relation yields

ˆΩ ∇δP ε • δu ε dx = d dt ˆΩ E ρ ε , ρ ε app dx -ˆΩ ∇P (ρ ε app ) • δu ε dx + I , (3.61)
where we have defined

I := ˆΩ ∇ • u ε app δP ε dx + ˆΩ ∂ t H (ρ ε app )δρ ε dx + ˆΩ ∇H (ρ ε app ) • (ρ ε u ε -ρ ε app u ε app )dx.
Let us work on this term for a while. We use the following Taylor expansion,

P ρ ε , ρ ε app := P (ρ ε ) -P (ρ ε app ) -P (ρ ε app ) δρ ε (3.62) = 1 2 (δρ ε ) 2 ˆ1 0 (1 -s) P ρ ε app + sδρ ε ds ,
and the fact that H (z) = P (z)/z according to (1.11), to obtain the next series of equalities:

I = ˆΩ ∇ • u ε app P (ρ ε app )δρ ε dx + ˆΩ ∇ • u ε app P ρ ε , ρ ε app dx + ˆΩ H (ρ ε app )∂ t ρ ε app δρ ε dx + ˆΩ H (ρ ε app )∇ρ ε app • (ρ ε u ε -ρ ε app u ε app )dx = ˆΩ ∇ • u ε app P (ρ ε app ) δρ ε dx -ˆΩ H (ρ ε app )ρ ε app ∇ • u ε app δρ ε dx + ˆΩ H (ρ ε app ) ∂ t ρ ε app + ∇ • (ρ ε app u ε app ) δρ ε dx + ˆΩ H (ρ ε app )∇ρ ε app • δu ε ρ ε dx + ˆΩ ∇ • u ε app P ρ ε , ρ ε app dx = ˆΩ H (ρ ε app )∇ρ ε app • δu ε ρ ε dx + ε ˆΩ H (ρ ε app )δρ ε R bl + εR ε dx + ˆΩ ∇ • u ε app P ρ ε , ρ ε app dx .
The last two terms in the above identity are small (in a sense to be made precise later). So, let us focus on the first term in the right hand side: we have

ˆΩ H (ρ ε app )∇ρ ε app • δu ε ρ ε dx -ˆΩ ∇P (ρ ε app ) • δu ε dx = ˆΩ H (ρ ε app )∇ρ ε app • δu ε δρ ε dx .
Inserting this expression into the last equality for I, from (3.61) we finally find

ˆΩ ∇δP ε • δu ε dx = d dt ˆΩ E ρ ε , ρ ε app dx + ˆΩ H (ρ ε app )∇ρ ε app • δu ε δρ ε dx + ε ˆΩ H (ρ ε app )δρ ε R bl + εR ε dx + ˆΩ ∇ • u ε app P ρ ε , ρ ε app dx . ( 3
d dt ˆΩ 1 2 ρ ε |δu ε | 2 + 1 ε 2 E ρ ε , ρ ε app dx + µ ˆΩ |∇ h δu ε | 2 dx + ε ˆΩ |∂ 3 δu ε | 2 dx + λ ˆΩ |∇ • δu ε | 2 dx ≤ 1 ε 2 ˆΩ δρ ε ∇G • δu ε dx - 1 ε 2 ˆΩ H (ρ ε app )∇ρ ε app • δu ε δρ ε dx - 1 ε ˆΩ H (ρ ε app )δρ ε R bl + εR ε dx - 1 ε 2 ˆΩ ∇ • u ε app P ρ ε , ρ ε app dx - 1 ε ˆΩ e 3 × δρ ε u ε app • δu ε dx -ˆΩ δρ ε ∂ t u ε app • δu ε dx + ˆΩ(ρ ε app u ε app -ρ ε u ε ) • ∇u ε app • δu ε dx -ˆΩ S bl • δu ε dx -ε ˆΩ S ε • δu ε dx - 1 ε ˆΩ x 3 ˆ1 0 ∂ 3 ρ (s x 3 ) ds (u bl 0,h,b ) ⊥ • δu ε h dx + 1 ε ˆΩ (1 -x 3 ) ˆ1 0 ∂ 3 ρ (1 -s(1 -x 3 )) ds (u bl 0,h,t ) ⊥ • δu ε h dx = 11 j=1 I j . (3.64)
Remark 3.9. In order to rigorously justify the relative entropy inequality (3.64), where the equality holds if the solutions are regular enough, one may either proceed as in [START_REF] Feireisl | Relative entropies, suitable weak solutions, and weak-strong uniqueness for the compressible Navier-Stokes system[END_REF], or use a regularization argument (see for instance [START_REF] Feireisl | Sun: Suitable weak solutions to the Navier-Stokes equations of compressible viscous fluids[END_REF] and [START_REF] Germain | Weak-strong uniqueness for the isentropic compressible Navier-Stokes system[END_REF]).

Our next goal is to bound each term appearing in the sum 11 j=1 I j in the right hand side of (3.64). Before doing that, let us remark that, since ρ 1 , ρ 2 ∈ L ∞ (Ω × R + ), up to restrict our attention to all ε ≤ ε 0 , with ε 0 depending on ρ 1 L ∞ t,x and ρ 2 L ∞ t,x , we can assume that -σ 2 ≤ ερ 1 + ε 2 ρ 2 ≤ σ 2 with σ > 0 as in (3.46). Consequently, we can suppose that 0 < ρ - app ≤ ρ ε app (x, t) ≤ ρ + app for all ε > 0 , with ρ - app = inf (0,1) ρ -σ and ρ + app = sup (0,1) ρ + σ. Then, in view of (3.45), we have the following control:

E ρ ε , ρ ε app (x, t) ≥ c |δρ ε (x, t)| 2 1 {|δρ ε |(•,t)<1} + |δρ ε (x, t)| γ 1 {|δρ ε |(•,t)≥1} . (3.65)
Resorting to the definitions (3.46), from (3.65) we derive the following lower bound. Lemma 3.10. There exist σ > 0 small enough (depending on inf (0,1) ρ) and a positive constant c > 0, independent of ε ∈ ]0, ε 0 ], such that, for almost all (x, t) ∈ Ω × R + , the following bound holds:

E ρ ε (x, t), ρ ε app (x, t) ≥ c [δρ ε ] 2 ess (x, t) + 1 Ωres(t) (x) . (3.66)
Proof. We divide the proof of the inequality into two steps. First we show that

E ρ ε (x, t), ρ ε app (x, t) ≥ c |δρ ε (x, t)| 2 1 {|δρ ε (x,t)|<1} implies the lower bound (3.67) E ρ ε (x, t), ρ ε app (x, t) ≥ c [δρ ε ] 2 ess (x, t).
For this, we just need to show that Ω ess (t) ⊆ {| |δρ ε (x, t)| < 1}. Let x ∈ Ω ess (t), then

- 3 2 σ ≤ -σ -ερ 1 (x, t) -ε 2 ρ 2 (x, t) < δρ ε (x, t) < -ερ 1 (x, t) -ε 2 ρ 2 (x, t) + σ ≤ 3 2 σ,
where we have used that

-σ 2 ≤ ερ 1 (x, t) + ε 2 ρ 2 (x, t) ≤ σ 2 .
By choosing σ such that σ < min 2/3 , inf (0,1) ρ , we deduce that |δρ ε (x, t)| < 1. Thus, (3.67) is proved.

Afterwards, we prove that for, x ∈ Ω res (t), one has

(3.68) E ρ ε (x, t), ρ ε app (x, t) ≥ c
, where c is a positive constant independent of ε, t and x. By the definition of Ω res (t), either ρ ε (x, t) ≤ ρ(x 3 ) -σ or ρ ε (x, t) ≥ ρ(x 3 ) + σ. Hence, since E(•, ρ ε app (x, t)) is strictly decreasing before ρ ε app (x, t) and strictly increasing after ρ ε app (x, t), we get E ρ ε (x, t), ρ ε app (x, t) ≥ E ρ(x 3 ) -σ, ρ ε app (x, t) if ρ ε (x, t) ≤ ρ(x 3 ) -σ, and E ρ ε (x, t), ρ ε app (x, t) ≥ E ρ(x 3 ) + σ, ρ ε app (x, t) if ρ ε (x, t) ≥ ρ(x 3 ) + σ. Now, by Taylor's formula, up to taking a smaller σ (which amounts to choosing a smaller ε 0 ), we have

E ρ(x 3 ) -σ, ρ ε app (x, t) ≥ H (ρ ε app (x, t)) 4 (-σ -ερ 1 (x, t) -ε 2 ρ 2 (x, t)) 2 ≥ H (ρ ε app (x, t))σ 2 16 E ρ(x 3 ) + σ, ρ ε app (x, t) ≥ H (ρ ε app (x, t)) 4 (σ -ερ 1 (x, t) -ε 2 ρ 2 (x, t)) 2 ≥ H (ρ ε app (x, t))σ 2 16 .
Then, using the uniform boundedness in time and space of ρ ε app and hypothesis (1.10), we get (3.68). The lemma is proved.

Notice that [|δρ ε |] ess is uniformly bounded. Next, we claim that there exists a constant C > 0 such that, for all T > 0 fixed, one has

(3.69) [δρ ε ] res L ∞ T (L p ) ≤ C ε 2/p ∀ p ∈ [1, γ] .
Indeed, by Hölder's inequality, the L ∞ control on ρ ε app and (3.49), we deduce ˆΩ

|[δρ ε ] res | ≤ ˆΩ [ρ ε ] res + ˆΩ ρ ε app res ≤ ˆΩ (ρ ε ) γ 1 Ωres 1/γ (L(Ω res )) 1/γ + C L(Ω res ) ≤ C ε 2 ,
which yields (3.69) for p = 1. As for the L γ norm, we write

(3.70) Ω res (t) = {0 < ρ ε (x, t) ≤ ρ(x 3 ) -σ} ∪ {ρ ε (x, t) ≥ ρ(x 3 ) + σ} .
For the first set, we just apply (3.49) again, since ρ ε is bounded therein. For the second set, we use the fact that, for a ≥ δ and b ≥ 0, with b ≤ b * , one has |a -b| γ ≤ (a + b) γ ≤ C δ,b * (a γ + 1). The case 1 < p < γ follows from interpolation.

3.3.2. Anisotropic Sobolev embedding. We introduce an anisotropic version of the standard Sobolev embedding Ḣ1 → L 6 . This estimate enables us to handle the anisotropy of the viscosity in the stability estimates below, see in particular the treatment of I 7 .

Lemma 3.11 (anisotropic Sobolev embedding). Let Ω = R 2 × (0, 1). There exists a universal constant C > 0 such that, for all κ > 0 and all u ∈ H 1 0 (Ω), one has

(3.71) u L 6 (Ω) ≤ C κ -1 2 ∇ h u L 2 (Ω) + κ ∂ 3 u L 2 (Ω) .
Proof of Lemma 3.11. Let κ > 0 and u ∈ H 1 0 (Ω). We first extend u by zero on R 3 \ Ω and still denote the extended function by u. Now u ∈ H 1 (R 3 ). We then consider the rescaled function

u κ (y h , y 3 ) = u y h κ 1 2
, κy 3 , (y h , y 3 ) ∈ R 3 .

By Sobolev's inequality [27, estimate (II.3.7)] for the whole space, there exists a universal constant C ∈ (0, ∞) such that

u κ L 6 (R 3 ) ≤ C ∇u κ L 2 (R 3 ) .
Estimate (3.71) then follows by a change of variables and the fact that u is zero outside the strip R 2 × (0, 1),

u κ L 6 (R 3 ) = u L 6 (Ω) , ∇u κ L 2 (R 3 ) = κ -1 2 ∇ h u L 2 (Ω) + κ ∂ 3 u L 2 (Ω)
. This concludes the proof.

Conclusion of the stability estimates.

Below we estimate every source term I j appearing in (3.64), for 0 < ε ≤ ε 0 (ε 0 is given by Lemma 3.10). For the terms I 1 , I 2 , I 3 , I 5 and I 6 we need to treat separately the cases γ ≥ 2 and 3/2 ≤ γ < 2, since we use different estimates, whereas the terms I 4 , I 8 , I 9 , I 10 and I 11 can be controlled in the same way for any γ ≥ 3/2. The term I 7 is more intricate; it is written as a sum of five terms: for some of them, we need to distinguish again the case γ ≥ 2 and 3/2 ≤ γ < 2. The easiest terms to handle are I 3 , I 4 , I 6 and I 9 . The terms I 1 and I 2 are combined with the Coriolis term I 5 ; For the remaining part of I 5 , we rely on Hardy's inequality, which is also useful to deal with I 7 , I 8 , I 10 and I 11 . The basic idea, borrowed from [START_REF] Bresch | Rotating fluids in a cylinder[END_REF], is that, whenever there is a boundary layer term G bl ( x 3 ε ), we gain one additional ε by using the decay of G bl in ζ. The price to pay is a ∂ 3 derivative on δu ε , which however can be swallowed by the third term in the left hand side of (3.64).

For every term, we decompose u ε app according to (3.34). The terms which require more care are those of order O(1), which involve in general u 0,h and u bl 0,h , except for I 7 where the product u ε app • ∇u ε app also involves u 1,3 and u bl 1,3 at order O(1). For the terms which are not of order O(1) the analysis can be always reduced to the case I 3 , I 4 , I 6 , I 9 , and the same estimates are used.

For every term involving a boundary layer, one has to equally consider the top and bottom boundary layers; again, for simplicity, we focus on the boundary layer at the bottom only. In the computations below, U and U bl generically denote remainder terms in the expansion for u ε app or its derivatives. The definition of these remainder terms may change from the estimate of one I i to another I j .

First, we deal with the terms for which estimates hold for any γ.

Term I 4 . We start by considering I 4 , when restricted to the essential set. Using (3.62), the assumptions on the pressure function and the fact that [|δρ ε |] ess is uniformly bounded, we can estimate

1 ε 2 ˆΩ ∇ • u ε app P (ρ ε , ρ ε app ) ess ≤ 1 ε 2 ∇ • u ε app L ∞ [δρ ε ] ess 2 L 2 ≤ C ε 1 ε 2 ˆΩ E ρ ε , ρ ε app ,
where we have used also that ∇

• u ε app = ε ∇ • u 1 + ∇ h • u bl 1,h .
Let us consider the integral over the residual set. By (3.62) again, we have [P ] res = P (ρ ε ) -P (ρ ε app ) res -P (ρ ε app )[δρ ε ] res . The second term can be easily controlled, in view of the uniform boundedness of ρ ε app and the L 1 estimate in (3.69). For the first term, we use decomposition (3.70): when ρ ε is bounded, the same argument as above applies. On the set {ρ ε ≥ ρ + σ}, instead, we use hypothesis (1.10), the uniform boundedness of ρ ε app and the controls in (3.49) to get 1

ε 2 ˆΩ |∇ • u ε app | P (ρ ε ) -P (ρ ε app ) 1 {ρ ε ≥ρ+σ} ≤ C ε ˆΩ P (ρ ε ) -P (ρ ε app ) 1 {ρ ε ≥ρ+σ} ≤ C ε [ρ ε ] res γ L γ + L (Ω res ) ≤ C ε .
Putting everything together, we finally infer that

(3.72) |I 4 | ≤ C ε + C ε 1 ε 2 ˆΩ E ρ ε , ρ ε app ,
where the last term will be handled by Grönwall's lemma.

Term I 9 . The control of I 9 is direct, as no ρ ε or δρ ε enter into play. We get

(3.73) |I 9 | ≤ ε S ε L 2 δu ε L 2 ≤ C ε K 2 (t) , where the function K 2 (t) = u ε (t) L 2 + u ε app (t) L 2 belongs to L 2 ([0, T )
) for all T > 0. Terms I 8 , I 10 , and I 11 . We deal with I 8 using Hardy's inequality. This gives

|I 8 | = ε ˆΩ x 3 ε S bl x 3 ε • δu ε x 3 (3.74) ≤ C δ ε ζ S bl 2 L 2 + δ ε ∂ 3 δu ε 2 L 2 ≤ C δ ε 2 + δε ∂ 3 δu ε 2 L 2 ,
for some small δ > 0, to be chosen later. The same holds for I 10 and I 11 : since ∂ 3 ρ is uniformly bounded, we have

|I 10 | = ε ˆΩ x 2 3 ε 2 ˆ1 0 ∂ 3 ρ(s x 3 ) ds (u bl 0,h,b ) ⊥ • δu ε h x 3 (3.75) ≤ C δ ε ζ 2 u bl 0,h,b 2 L 2 + δε ∂ 3 δu ε h 2 L 2 ≤ C δ ε 2 + δε ∂ 3 δu ε h 2 L 2 |I 11 | = ε ˆΩ (1 -x 3 ) 2 ε 2 ˆ1 0 ∂ 3 ρ(1 -s(1 -x 3 )) ds (u bl 0,h,t ) ⊥ • δu ε h 1 -x 3 (3.76) ≤ C δ ε η 2 u bl 0,h,b 2 L 2 + δε ∂ 3 δu ε h 2 L 2 ≤ C δ ε 2 + δε ∂ 3 δu ε h 2 L 2 .
In the estimates above, we have used the fact that the terms

ζ S bl 2 L 2 , ζ 2 u bl 0,h,b 2
L 2 and η 2 u bl 0,h,t 2
L 2 are O(ε 2
). We consider now the terms whose bounds must be treated differently if γ ≥ 2 or 3/2 ≤ γ < 2.

Term I 3 . First of all, observe that R bl (x 3 /ε)

2 L 2 x = O(ε). Thus, we can estimate 1 ε ˆΩ H (ρ ε app ) [δρ ε ] ess R bl + ε R ε ≤ C 1 ε R bl + ε R ε L 2 [δρ ε ] ess L 2 ≤ C ε + C ε 2 ˆΩ E ρ ε , ρ ε app ,
where we have used also (3.66). As for the residual part, in view of (3.69), we can argue in exactly the same way if γ ≥ 2. If 3/2 ≤ γ < 2, instead, we put the L ∞ norm on the remainder terms and use the L 1 bound of (3.69) to get

1 ε ˆΩ H (ρ ε app ) [δρ ε ] res R bl + ε R ε ≤ C ε .
In any case, in the end we arrive at the bound

(3.77) |I 3 | ≤ C ε + C ε 2 ˆΩ E ρ ε , ρ ε app .
Term I 6 . Once again, we use the decomposition of δρ ε into essential and residual parts. For the term involving the essential part, thanks to Young's inequality and to the controls (3.53) and (3.66), one has

ˆΩ [δρ ε ] ess ∂ t u ε app • δu ε ≤ [δρ ε ] ess L 2 ∂ t u ε app L ∞ δu ε L 2 ≤ C ε 2 K 1 (t) + 1 ε 2 E ρ ε , ρ ε app ,
where the function

K 1 = u ε 2 L 2 + u ε app 2
L 2 belongs to L 1 ([0, T )) for all T > 0. Next, let us consider the term involving the residual part: when γ ≥ 2, we can argue exactly as above, in view of (3.69). If instead 3/2 ≤ γ < 2, we start by writing

ˆΩ [δρ ε ] res ∂ t u ε app • δu ε = ˆΩ [ρ ε ] res ∂ t u ε app • δu ε - ˆΩ ρ ε app res ∂ t u ε app • δu ε .
For the second term, we use the uniform boundedness of ρ ε app and estimate (3.49) to gather, for some function K 2 ∈ L 2 ([0, T )) for all T > 0, the inequality

ˆΩ ρ ε app res ∂ t u ε app • δu ε ≤ C δu ε L 2 (L(Ω res )) 1/2 ≤ C ε K 2 (t) .
For the term involving [ρ ε ] res , we use decomposition (3.70) for the residual set. The integral over the first set can be treated exactly as just done for ρ ε app (because ρ ε is uniformly bounded therein). Concerning the integral over the second set, we have

ˆΩ ρ ε 1 {ρ ε ≥ρ+σ} ∂ t u ε app • δu ε ≤ C √ ρ ε δu ε L 2 ˆΩ ρ ε 1 {ρ ε ≥ρ+σ } 1/2 ≤ C √ ρ ε δu ε 2 L 2 + C ε 2
, since the last integral in the first line can be bounded by the integral over the residual set, for which we can use (3.49). Let us introduce the following notation: we set δ 2 -(γ) = 1 if 3/2 ≤ γ < 2, δ 2 -(γ) = 0 otherwise. In the end, from the previous computations we get

|I 6 | ≤ C ε (ε K 1 (t) + δ 2 -(γ) K 2 (t) + δ 2 -(γ) ε) (3.78) + C 1 ε 2 ˆΩ E ρ ε , ρ ε app + δ 2 -(γ) √ ρ ε δu ε 2 L 2 .
Terms I 1 , I 2 and I 5 . Terms I 1 , I 2 and I 5 have to be combined together, enabling to see a cancellation at the highest order in ε. Such a cancellation is already a key point in [START_REF] Feireisl | Multiple scales and singular limits for compressible rotating fluids with general initial data[END_REF]. After setting U := u ε app -(u 0,h + u bl 0,h , 0), we can write

I 1 + I 2 + I 5 = 1 ε 2 ˆΩ δρ ε ∇G • δu ε - 1 ε 2 ˆΩ H (ρ ε app ) ∂ 3 ρ δu ε 3 δρ ε - 1 ε ˆΩ H (ρ ε app ) ∇ρ 1 • δu ε δρ ε - 1 ε ˆΩ δρ ε P (ρ) ρ (∇ ⊥ h ρ 1 ) ⊥ • δu ε h - 1 ε ˆΩ δρ ε (u bl 0,h ) ⊥ • δu ε h -ˆΩ δρ ε e 3 × U (x h , x 3 , x 3 ε , 1-x 3 ε , t) • δu ε .
Notice that H (ρ) = P (ρ) ρ and (∇ ⊥ h ρ 1 ) ⊥ = -∇ h ρ 1 . Moreover, from (3.1) we get ρ ∇G = P (ρ) ∇ρ .

Therefore, we find

I 1 + I 2 + I 5 = - 1 ε 2 ˆΩ H (ρ ε app ) -H (ρ) ∂ 3 ρδu ε 3 δρ ε - 1 ε ˆΩ H (ρ)∂ 3 ρ 1 δu ε 3 δρ ε - 1 ε ˆΩ H (ρ ε app ) -H (ρ) ∇ρ 1 • δu ε δρ ε - 1 ε ˆΩ δρ ε (u bl 0,h ) ⊥ • δu ε h -ˆΩ δρ ε U ⊥ h (x h , x 3 , x 3 ε , 1-x 3 ε , t) • δu ε h = J 1 + J 2 + J 3 + J 4 + J 5 .
Using a Taylor expansion for h(z) = H (z) with integral remainder, we can write

J 1 + J 2 = - 1 ε ˆΩ h (ρ) ρ 1 ∂ 3 ρ + h(ρ) ∂ 3 ρ 1 δu ε 3 δρ ε -ˆΩ ρ 2 1 ˆ1 0 (1 -s)h (ρ + sερ 1 )ds ∂ 3 ρ δu ε 3 δρ ε = -ˆΩ ρ 2 1 ˆ1 0 (1 -s)h (ρ + sερ 1 )ds ∂ 3 ρ δu ε 3 δρ ε ,
where we have used (3.6) in the last equality. Since ρ 1 , ρ and ∂ 3 ρ are L ∞ t,x , in view of (1.10) the control of J 1 + J 2 becomes similar to the one exhibited for I 6 . In the same way, after noticing that ∇ h ρ 1 and ε -1 H (ρ ε app ) -H (ρ) are uniformly bounded in time and space, the control of J 3 is obtained. Then, J 1 + J 2 and J 3 verify estimate (3.78). The same can be said about J 5 , because also U h belongs to L ∞ t,x . Therefore, it remains to deal with J 4 , for which we rely on Hardy's inequality. More precisely, let us start, as usual, by dealing with the essential part: we have

1 ε ˆΩ [δρ ε ] ess (u bl 0,h ) ⊥ • δu ε h = ˆΩ [δρ ε ] ess x 3 ε (u bl 0,h ) ⊥ • δu ε h x 3 ≤ [δρ ε ] ess L 2 ζ u bl 0,h (t, x h , ζ) L ∞ t,x,ζ ∂ 3 δu ε h L 2 ≤ C ε 2 ˆΩ E ρ ε , ρ ε app + C ε 2 ζ u bl 0,h 2
L ∞ t,x,ζ ∂ 3 δu ε h 2 L 2 .
Notice that, for ε small enough, the second term can be swallowed by the third term in the left hand side of (3.64). As for the control of the residual part, suppose that γ ≥ 2 for a while: in this case, we can argue in the exact same way and obtain, in view of (3.69), that

1 ε ˆΩ [δρ ε ] res (u bl 0,h ) ⊥ • δu ε h ≤ [δρ ε ] res L 2 ζ u bl 0,h L ∞ t,x,ζ ∂ 3 δu ε h L 2 ≤ C ε 2 ˆΩ E ρ ε , ρ ε app + C ε 2 ζ u bl 0,h 2
L ∞ t,x,ζ ∂ 3 δu ε h 2 L 2 .
The case 3/2 ≤ γ < 2 is slightly more involved. The control over {0 < ρ ε ≤ ρ -σ} does not present any special difficulty, since we have uniform bounds for ρ ε (and obviously for ρ ε app ) on that set: then, we can argue as for controlling the essential part. Hence, let us focus on {ρ ε ≥ ρ + σ}. First of all, using that √ a + b ≤ √ a + √ b, we notice that (3.79)

1 ε ˆΩ [δρ ε ] res (u bl 0,h ) ⊥ • δu ε h ≤ 1 ε ˆ δρ ε 1 {ρ ε ≥ρ+σ} u bl 0,h √ ρ ε |δu ε h | + 1 ε ˆ δρ ε 1 {ρ ε ≥ρ+σ} u bl 0,h ρ ε app |δu ε h | .
For the first term in the right hand side of (3.79), we proceed in the following way:

ˆΩ δρ ε 1 {ρ ε ≥ρ+σ} u bl 0,h √ ρ ε |δu ε h | ≤ [δρ ε ] res 1/2 L γ u bl 0,h L ∞ √ ρ ε δu ε h L 2 (L(Ω res )) 1/q ,
where 1/(2γ) + 1/2 + 1/q = 1. Using (3.65) and (3.66), we deduce that 1 ε

ˆΩ δρ ε 1 {ρ ε ≥ρ+σ} u bl 0,h √ ρ ε |δu ε h | ≤ ε 1/γ+2/q-1 1 ε 2 E ρ ε , ρ ε app 1/(2γ)+1/q u bl 0,h L ∞ √ ρ ε δu ε h L 2 = 1 ε 2 E ρ ε , ρ ε app 1/2 u bl 0,h L ∞ √ ρ ε δu ε h L 2 .
After applying Young's inequality, this term can be controlled by Grönwall's lemma in the final estimate. For the last term in (3.79), we argue in the following way:

ˆ δρ ε 1 {ρ ε ≥ρ+σ} u bl 0,h ρ ε app |δu ε h | = ε ˆ δρ ε 1 {ρ ε ≥ρ+σ} x 3 ε u bl 0,h ρ ε app δu ε h x 3 ≤ ε ˆ √ ρ ε 1 {ρ ε ≥ρ+σ} x 3 ε u bl 0,h ρ ε app δu ε h x 3 + ρ ε app 1 {ρ ε ≥ρ+σ} x 3 ε u bl 0,h 1 x 3 δu ε h ≤ ε ∂ 3 δu ε h L 2 ζ u bl 0,h L ∞ × ρ ε app 1/2 L ∞ [ρ ε ] res 1/2
L γ (L(Ω res )) 1/q + ρ ε app L ∞ (L(Ω res )) 1/2 , where q is defined as above. Notice that, in view of (3.49), we have [ρ ε ] res L γ = O ε 2/γ and L(Ω res ) = O ε 2 . Therefore, we finally find

1 ε ˆ δρ ε 1 {ρ ε ≥ρ+σ} u bl 0,h ρ ε app |δu ε h | ≤ C δ ε + δ ε ∂ 3 δu ε h 2 L 2 .
In the end, we deduce the following control:

(3.80)

|I 1 + I 2 + I 5 | ≤ C ε 2 ˆΩ E ρ ε , ρ ε app + Cδ 2 -(γ) √ ρ ε δu ε h 2 L 2 + Cε(εK 1 (t) + δ 2 -(γ)K 2 (t) + δ 2 -(γ)) + (Cε 2 + δ 2 -(γ)δε ) ∂ 3 δu ε h 2 L 2
, where the last term in the right hand side can be absorbed into the left hand side of the relative entropy inequality (3.64).

Finally, let us deal with I 7 . Term I 7 . We start by considering the following decomposition:

I 7 = -ˆΩ δρ ε u ε app • ∇u ε app • δu ε - 1 ε ˆΩ ρ ε δu ε 3 ∂ ζ u bl 0,h ( x 3 ε ) • δu ε h -ˆΩ ρ ε δu ε 3 ∂ ζ u bl 1 ( x 3 ε ) • δu ε -ˆΩ ρ ε δu ε h • ∇ h (u 0,h + u bl 0,h ) • δu ε h -ε ˆΩ ρ ε δu ε • U • δu ε = J 6 + J 7 + J 8 + J 9 + J 10 ,
where εU = ε ∇u 1 + ∇ h 0 u bl 1 is the remainder term in the expansion for ∇u ε app . The first term J 6 can be handled as done with I 6 . Indeed, one has

u ε app • ∇u ε app = u 0,h + u bl 0,h • ∇ h (u 0,h + u bl 0,h ) + u 1,3 + u bl 1,3 ∂ ζ u bl 0,h ( x 3 ε ) + h.o.t. ,
where h.o.t. represents higher order terms in ε. Then u ε app • ∇u ε app is uniformly bounded in L ∞ t,x . Therefore, J 6 verifies an inequality similar to (3.78) above. The terms J 8 , J 9 and J 10 can be simply bounded as follows:

|J 8 | ≤ C ∂ ζ u bl 1 L ∞ t,x √ ρ ε δu ε 2 L 2 , |J 9 | ≤ C ∇ h (u 0,h + u bl 0,h ) L ∞ t,x √ ρ ε δu ε 2 L 2 , |J 10 | ≤ C ε U L ∞ t,x √ ρ ε δu ε 2 L 2 .
We remark that these estimates holds for γ ≥ 3/2. We now focus on the remaining term J 7 , which is the most difficult one to deal with. The difficulties come from the need to gain smallness in ε (by using Hardy's inequality as above), from the low integrability of the residual part and from the fact that this term is quadratic in δu ε . We first decompose (3.81) ρ ε = [δρ ε ] ess + [δρ ε ] res + ρ ε app .

The essential part is easy to bound: owing to the boundedness of ρ ε on that set and to an application of Hardy's inequality, we get

1 ε ˆΩ [δρ ε ] ess δu ε 3 ∂ ζ u bl 0,h ( x 3 ε ) • δu ε h ≤ C ε ζ 2 ∂ ζ u bl 0,h L ∞ t,x δu ε 3 x 3 L 2 δu ε h x 3 L 2 ≤ C ε ∂ 3 δu ε 3 L 2 ∂ 3 δu ε h L 2 ≤ C ε 3/2 ∂ 3 δu ε h 2 L 2 + C ε 1/2 ∂ 3 δu ε 3 2 L 2 .
The control of the part involving ρ ε app is similar, so let us turn to the residual part. Two different estimates are computed if γ is larger or smaller than the critical exponent 2. For γ ≥ 2, we write, for α ∈ (0, 1) to be chosen later on, We use the same technique for δu ε h with β ∈ (0, 1). Then, choosing α, β such that

δu ε 3 = (δu ε 3 ) 1-α (δu ε 3 ) α x α
α + β = 1 2 ,
we have, for all δ > 0 to be chosen later,

1 ε ˆΩ [δρ ε ] res δu ε 3 ∂ ζ u bl 0,h ( x 3 ε ) • δu ε h = ε α+β-1 ˆΩ [δρ ε ] res (δu ε 3 ) 1-α (δu ε 3 ) α x α 3 x α+β 3 ε α+β ∂ ζ u bl 0,h ( x 3 ε ) • (δu ε h ) 1-β (δu ε h ) β x β 3 ≤ ε -1/2 [δρ ε ] res L 2 ∇δu ε 3 1-α L 2 ∂ 3 δu ε 3 α L 2 ∇δu ε h 1-β L 2 ∂ 3 δu ε h β L 2 ζ α+β ∂ ζ u bl 0,h L ∞ ≤ Cε 1/2 1 ε 2 ˆΩ E(ρ ε , ρ ε app ) 1/2 ∇δu ε 3 L 2 ∇δu ε h L 2 ≤ C δ ε 2 K 1 (t) ˆΩ E(ρ ε , ρ ε app ) + δ ε ∇ h δu ε h 2 L 2 + δ ε ∂ 3 δu ε h 2 L 2 ,
where K 1 (t) = ∇u ε 3 (t) 2 L 2 + ∇u ε app,3 (t) 2 L 2 belongs to L 1 ([0, T )) for all T > 0. In the second inequality we have used the lower bound E(ρ ε (x, t), ρ ε app (x, t)) ≥ c |δρ ε (x, t)| 2 , which comes from (3.65) when γ ≥ 2. For 3/2 ≤ γ < 2, we use the same argument as in the case γ ≥ 2 for δu ε 3 . The control of δu ε h , instead, is done via the anisotropic Sobolev embedding given of Lemma 3.11. Hence, for α ∈ [0, 1] such that 2 -3 γ = α, Hölder's inequality gives, using (3.82) for δu ε 3 and (3.71) for δu ε h with κ = ε

(-1 2 + 1 γ ) + , (3.83) 1 ε ˆΩ [δρ ε ] res δu ε 3 ∂ ζ u bl 0,h ( x 3 ε ) • δu ε h = ε α-1 ˆΩ [δρ ε ] res (δu ε 3 ) 1-α (δu ε 3 ) α x α 3 x α 3 ε α ∂ ζ u bl 0,h ( x 3 ε ) • δu ε h ≤ Cε α-1 [δρ ε ] res L γ ∇δu ε 3 1-α L 2 ∂ 3 δu ε 3 α L 2 × κ -1 2 ∇ h δu ε h L 2 + κ ∂ 3 δu ε h L 2 ζ α ∂ ζ u bl 0,h L ∞ ≤ Cε 1-1 γ ∇δu ε 3 L 2 κ -1 2 ∇ h δu ε h L 2 + κ ∂ 3 δu ε h L 2 ≤ min(µ,λ) 10 ∇δu ε 3 2 L 2 + C(µ, λ)ε 2-2 γ κ -1 ∇ h δu ε h 2 L 2 + κ 2 ∂ 3 δu ε h 2 L 2
≤ min(µ,λ) 10 ∇δu ε 3 2

L 2 + C(µ, λ)ε ( 5 2 -3 γ ) -∇ h δu ε h 2 L 2 + C(µ, λ)ε 1 + ∂ 3 δu ε h 2 L 2 .
Hence we can swallow the whole right hand side on condition that γ > 6/5 (which is the case, since γ ≥ 3/2) and ε is sufficiently small. To put it in a nutshell, we obtain the following bound on J 7 :

|J 7 | ≤ C δ ε 2 K 1 (t) ˆΩ E ρ ε , ρ ε app + (δ ε + Cε 3 2 + δ 2 -(γ)C(µ, λ)ε 1 + ) ∂ 3 δu ε h 2 L 2 + δ ε + δ 2 -(γ)C(µ, λ)ε ( 5 2 -3 γ ) - ∇ h δu ε h 2 L 2
+ δ 2 -(γ) min(µ,λ) 10 + Cε

1 2 ∇δu ε 3 2 L 2 .
Therefore, we finally get the following estimate for I 7 :

(3.84)

|I 7 | ≤ Cε ε K 1 (t) + δ 2 -(γ)K 2 (t) + δ 2 -(γ)ε + C + C δ K 1 (t) ε 2 ˆΩ E ρ ε , ρ ε app + (δ 2 -(γ)C + C 1 + C 2 ε) √ ρ ε δu ε 2 L 2 + (δε + Cε 3 2 + δ 2 -(γ)C(µ, λ)ε 1 + ) ∂ 3 δu ε h 2 L 2 + δ ε + δ 2 -(γ)C(µ, λ)ε ( 5 2 -3 γ ) - ∇ h δu ε h 2 L 2 + δ 2 -(γ) min(µ,λ) 10 + Cε 1 2 ∇δu ε 3 2 L 2 .
Remark 3.12. The anisotropic Sobolev embedding in Lemma 3.11 can be used to provide better estimates only for γ small. For instance, in (3.79), using Lemma 3.11 we get a remainder term of order ε α , with 0 < α < 1 for 3/2 ≤ γ < 2 and α > 1 only for γ < 12/11, while by using the smallness of the Lebesgue measure of Ω res we get a remainder term of order ε for 3/2 ≤ γ < 2.

In the end, summing up our estimates, we get from (3.64) the following differential inequality: there exist functions C 1 (t), C 2 (t) ∈ L 1 ([0, T )), and constants C 3 > 0 and

  30) at ζ = 0. It remains to choose the boundary layer term u bl 1,h . The specifications for the boundary layer term u bl 1,h are that it is exponentially decaying to 0 for ζ → ∞ and satisfies (3.30) at the boundary ζ = 0. Hence, we define u bl 1,h,b in the following way: for all ζ ∈ [0, ∞) and x h ∈ R 2 ,

3 x α 3 and

 33 then apply Sobolev's and Hardy's inequalities: this yields (δu ε3 )

  Hence, estimate (2.48) above tells us that the quantity θ(t) := sup

	τ ∈[0,t]	δr(τ ) H 1

  by virtue of (2.41) and ∂ t ∇u 2 ∈ L T (L 2 ) by definition of X T . Hence, using Sobolev embedding and interpolation (2.30), we deduce that ∂ t u 2 ∈ L

		4/3
	16/11 T	(L 4 ). Therefore, we can bound

Acknowledgement

The three authors are partially supported by the project BORDS (ANR-16-CE40-0027-01) operated by the French National Research Agency (ANR). The second and third authors are partially supported by the project SingFlows (ANR-18-CE40-0027) operated by the French National Research Agency (ANR). The work of the second author is also partially supported by the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissement d'Avenir" (ANR-11-IDEX-0007). The first and the third authors acknowledge financial support from the IDEX of the University of Bordeaux for the BOLIDE project. The first author is partially supported by the Starting Grant project "Analysis of moving incompressible fluid interfaces"(H2020-EU.1.1.-639227) operated by the European Research Council (ERC).

ε 0 ∈ (0, 1), such that, for all ε ∈ (0, ε 0 ), all t ∈ (0, T ) and all δ > 0, one has

Let us stress that C 1 (t), C 2 (t), C 3 and ε 0 do not depend on ε. The quantities these constants depend on have been written explicitly in the computations above; in particular, C 1 (t) and C 2 (t) contains the functions K 1 (t) and K 2 (t). Choosing δ small enough and using the identity

, the last three terms in (3.85) can be swallowed in the left hand side. The estimate in Theorem 2 follows from Grönwall's lemma.