# Global-in-time well-posedness of the compressible Navier-Stokes equations with striated density

Xian Liao, Sagbo Marcel Zodji

<figure><img src="image_1.png" /></figure>

## To cite this version:

Xian Liao, Sagbo Marcel Zodji. Global-in-time well-posedness of the compressible Navier-Stokes equations with striated density. 2024. hal-04586708

HAL Id: hal-04586708

https://hal.science/hal-04586708

Preprint submitted on 24 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Global-in-time well-posedness of the compressible Navier-Stokes equations with striated density

Xian Liao, Sagbo Marcel Zodji

<figure><img src="image_1.png" /></figure>

## To cite this version:

Xian Liao, Sagbo Marcel Zodji. Global-in-time well-posedness of the compressible Navier-Stokes equations with striated density. Revista Matemática Iberoamericana, In press, 10.4171/RMI/1561. hal-04586708v2

HAL Id: hal-04586708

https://hal.science/hal-04586708v2

Submitted on 17 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
GLOBAL-IN-TIME WELL-POSEDNESS OF THE COMPRESSIBLE
NAVIER-STOKES EQUATIONS WITH STRIATED DENSITY

XIAN LIAO AND SAGBO MARCEL ZODJI

ABSTRACT. We first show local-in-time well-posedness of the compressible Navier-Stokes equations, assuming striated regularity while no other smoothness or smallness conditions on the initial density. With these local-in-time solutions served as blocks, for less regular initial data where the vacuum is permitted, the global-in-time well-posedness follows from the energy estimates and the propagated striated regularity of the density function, if the bulk viscosity coefficient is large enough in the two dimensional case. The global-in-time well-posedness holds also true in the three dimensional case, provided with large bulk viscosity coefficient together with small initial energy. This solves the density-patch problem in the exterior domain for the compressible model with $W^{2,p}$-Interfaces. Finally, the singular incompressible limit toward the inhomogenous incompressible model when the bulk viscosity coefficient tends to infinity is obtained.

# 1. INTRODUCTION

In this paper, we establish the existence and uniqueness of global-in-time weak solutions of compressible viscous flows, and at the same time we investigate the dynamics of density-interfaces in dimension $d \in \{2, 3\}$. More precisely, we consider the following compressible Navier-Stokes equations describing the motion of compressible viscous fluids:

$$ (1.1) \qquad \begin{cases} \partial_t \rho + \operatorname{div}(\rho u) = 0, \\ \partial_t (\rho u) + \operatorname{div}(\rho u \otimes u) + \nabla P(\rho) = \mu \Delta u + (\mu + \lambda) \nabla \operatorname{div} u. \end{cases} $$

Here $\mu > 0$ represents the dynamic viscosity, and $\lambda > 0$ stands for the kinetic viscosity. In the present paper, $\mu$ is some fixed positive constant while the constant $\lambda$ may become very large. For notational simplicity, we introduce the so-called bulk viscosity coefficient

$$ \nu = 2\mu + \lambda, $$

which tends to infinity when $\lambda \to \infty$. We always assume that our fluids are (strictly) viscous:

$$ \nu \geq \underline{\nu} > 0, $$

where $\underline{\nu}$ is a fixed positive constant.

In the above, $t > 0$, $x \in \mathbb{R}^d$, $d = 2, 3$ denote the time and space variables respectively. The notations $\rho = \rho(t, x) \geq 0$ and $u = u(t, x) \in \mathbb{R}^d$ represent, respectively, the density and velocity of the compressible fluid, which serve as the unknowns in the problem. Meanwhile, $P = P(\rho)$ is a given smooth (in this paper we assume $P \in \mathcal{C}^2(\mathbb{R}; \mathbb{R})$) and increasing function of the density (that is, $P'(\rho) > 0$). The system (1.1) is supplemented with initial data

$$ (1.2) \qquad (\rho, \rho u)|_{t=0} = (\rho_0, \rho_0 u_0), $$

which satisfy

$$ (1.3) \qquad \rho_0 \geq 0, \quad \rho_0 \in L^\infty(\mathbb{R}^d; [0, \infty)), \quad \rho_0 - \tilde{\rho} \in L^2(\mathbb{R}^d; \mathbb{R}), \quad u_0 \in H^1(\mathbb{R}^d; \mathbb{R}^d), $$

where $\tilde{\rho} > 0$ is some given positive equilibrium state of the density.

<footnote>Date: May 21, 2024.</footnote>

<footnote n="2020">2020 Mathematics Subject Classification. 35A01, 35A02, 35Q30, 35R35, 76N10.</footnote>

<footnote>Key words and phrases. Compressible Navier-Stokes Equations, Tangential Regularity, Incompressible limit, Vacuum, Density Patch Problem, Free Boundary Problem in Exterior Domain.</footnote>
# GLOBAL-IN-TIME WELL-POSEDNESS OF THE COMPRESSIBLE NAVIER-STOKES EQUATIONS WITH STRIATED DENSITY

XIAN LIAO AND SAGBO MARCEL ZODJI

**ABSTRACT.** We first show local-in-time well-posedness of the compressible Navier-Stokes equations, assuming striated regularity while no other smoothness or smallness conditions on the initial density. With these local-in-time solutions served as blocks, for less regular initial data where the vacuum is permitted, the global-in-time well-posedness follows from the energy estimates and the propagated striated regularity of the density function, if the bulk viscosity coefficient is large enough in the two dimensional case. The global-in-time well-posedness holds also true in the three dimensional case, provided with large bulk viscosity coefficient together with small initial energy. This solves the density-patch problem in the exterior domain for the compressible model with $W^{2,p}$-Interfaces. Finally, the singular incompressible limit toward the inhomogenous incompressible model when the bulk viscosity coefficient tends to infinity is obtained.

## 1. INTRODUCTION

In this paper, we establish the existence and uniqueness of global-in-time weak solutions of compressible viscous flows, and at the same time we investigate the dynamics of density interface in dimension $d \in \{2, 3\}$. More precisely, we consider the following compressible Navier-Stokes equations describing the motion of compressible viscous fluids:

$$ (1.1) \qquad \begin{cases} \partial_t \rho + \operatorname{div}(\rho u) = 0, \\ \partial_t (\rho u) + \operatorname{div}(\rho u \otimes u) + \nabla P(\rho) = \mu \Delta u + (\mu + \lambda) \nabla \operatorname{div} u. \end{cases} $$

Here $\mu > 0$ represents the dynamic viscosity, and $\lambda > 0$ stands for the kinetic viscosity. In the present paper, $\mu$ is some fixed positive constant while the constant $\lambda$ may become very large. For notational simplicity, we introduce the so-called bulk viscosity coefficient

$$ \nu = 2\mu + \lambda, $$

which tends to infinity when $\lambda \to \infty$. We always assume that our fluids are (strictly) viscous:

$$ \nu \geq \underline{\nu} > 0, $$

where $\underline{\nu}$ is a fixed positive constant.

In the above, $t > 0$, $x \in \mathbb{R}^d$, $d = 2, 3$ denote the time and space variables respectively. The notations $\rho = \rho(t, x) \geq 0$ and $u = u(t, x) \in \mathbb{R}^d$ represent, respectively, the density and velocity of the compressible fluid, which serve as the unknowns in the problem. Meanwhile, $P = P(\rho)$ is a given smooth (in this paper we assume $P \in \mathcal{C}^2(\mathbb{R}, \mathbb{R})$).

The system (1.1) is supplemented with initial data

$$ (1.2) \qquad (\rho, \rho u)|_{t=0} = (\rho_0, \rho_0 u_0), $$

which satisfy

$$ (1.3) \qquad \rho_0 \geq 0, \quad \rho_0 \in L^\infty(\mathbb{R}^d; [0, \infty)), \quad \rho_0 - \tilde{\rho} \in L^2(\mathbb{R}^d; \mathbb{R}), \quad u_0 \in H^1(\mathbb{R}^d; \mathbb{R}^d), $$

where $\tilde{\rho} > 0$ is some given positive equilibrium state of the density.

<footnote>Date: August 17, 2025.</footnote>

<footnote n="2020">2020 Mathematics Subject Classification. 35A01, 35A02, 35Q30, 35R35, 76N10.</footnote>

<footnote>Key words and phrases. Compressible Navier-Stokes Equations, Tangential Regularity, Incompressible limit, Vacuum, Density Patch Problem, Free Boundary Problem in Exterior Domain.</footnote>
1.1. **Striated regularity.** We assume further striated regularity with respect to a given nondegenerate family of vector fields for the initial density $\rho_0$ in this paragraph.

We first introduce some notations, based on [10]. For some $p \in (d, \infty)$, $L^{\infty,p}(\mathbb{R}^d; \mathbb{R}^d)$ denotes the vector space of bounded vector fields with gradients in $L^p(\mathbb{R}^d; \mathbb{R}^{d \times d})$. From now on we denote the Lebesgue spaces $L^p(\mathbb{R}^d; \mathbb{R}^n)$ resp. Sobolev spaces $H^s(\mathbb{R}^d; \mathbb{R}^n)$ with $p \in [1, \infty]$, $s \in \mathbb{R}$ and $n \in \mathbb{N}^*$, simply by $L^p(\mathbb{R}^d)$ or $L^q$, resp. $H^s(\mathbb{R}^d)$ or $H^q$ with an abuse of notations. We have defined

$$L^{\infty,p}(\mathbb{R}^d) = \{ Y \in L^{\infty}(\mathbb{R}^d) \mid \|Y\|_{L^{\infty,p}(\mathbb{R}^d)} := \|Y\|_{L^{\infty}(\mathbb{R}^d)} + \|\nabla Y\|_{L^p(\mathbb{R}^d)} < \infty \}.$$

We define the norm $\| \cdot \|_{L^{\infty,p}}$ for a family of vector fields $\mathcal{Y} = (Y_1, Y_2, \dots, Y_m) \subset L^{\infty,p}(\mathbb{R}^d)$, $m \in \mathbb{N}$ as

$$\|\mathcal{Y}\|_{L^{\infty,p}(\mathbb{R}^d)} := \sup_{1 \le v \le m} \|Y_v\|_{L^{\infty,p}(\mathbb{R}^d)}.$$

**Definition 1.1** (Nondegeneracy). Let $\mathcal{Y} = (Y_1, Y_2, \dots, Y_m) \subset L^{\infty,p}(\mathbb{R}^d)$ be a family of $m$ vector fields with $m \ge d-1$ and $p \in (d, \infty)$. We say that $\mathcal{Y}$ is nondegenerate if it satisfies the following property

$$I(\mathcal{Y}) := \inf_{x \in \mathbb{R}^d} \sup_{\Upsilon \in \Upsilon_{d-1}^m} \left| \bigwedge_{i=1}^{d-1} Y_{\Upsilon(i)} \right|^{\frac{1}{d-1}} > 0.$$

Above $\Upsilon \in \Upsilon_{d-1}^m$ means that $\Upsilon = (v_1, v_2, \dots, v_{d-1})$ with each $v_i \in \{1, \dots, m\}$ and $v_i < v_j$ for $i < j$, $Y_{\Upsilon} := (Y_{v_1}, Y_{v_2}, \dots, Y_{v_{d-1}})$, while the symbol $\bigwedge_{i=1}^{d-1} Y_{\Upsilon(i)}$ stands for the unique element of $\mathbb{R}^d$ such that

$$(\bigwedge_{i=1}^{d-1} Y_{\Upsilon(i)}) \cdot Z = \det(Y_{v_1}, Y_{v_2}, \dots, Y_{v_{d-1}}, Z), \quad \forall Z \in \mathbb{R}^d.$$

**Definition 1.2** (Striated regularity with respect to a nondegenerate family of vector fields). Let $Y \in L^{\infty,p}(\mathbb{R}^d)$, $p \in (d, \infty)$ be a (single) vector field, and $\mathcal{Y} = (Y_1, Y_2, \dots, Y_m) \subset L^{\infty,p}(\mathbb{R}^d)$ be a nondegenerate family of vector fields with $m \ge d-1$.

(a) A function $g \in L^{\infty}(\mathbb{R}^d)$ is said to be of class $L^p(\mathbb{R}^d)$ along $Y$, if

$$g \in L_Y^p(\mathbb{R}^d) := \{g \in L^{\infty}(\mathbb{R}^d) \mid \mathrm{div}(gY) \in L^p(\mathbb{R}^d)\}.$$

We define the derivative of the function $g$ along $Y$ as follows

$$\partial_Y g := \mathrm{div}(gY) - g \mathrm{div} Y,$$

and hence we can equivalently define

$$L_Y^p(\mathbb{R}^d) = \{g \in L^{\infty}(\mathbb{R}^d) \mid \partial_Y g \in L^p(\mathbb{R}^d)\}.$$

(b) A function $g \in L^{\infty}(\mathbb{R}^d)$ is said to be of class $L^p(\mathbb{R}^d)$ along the family $\mathcal{Y}$, if

$$g \in L_Y^p(\mathbb{R}^d) := \bigcap_{1 \le v \le m} L_{Y_v}^p(\mathbb{R}^d),$$

and we equip the space $L_Y^p(\mathbb{R}^d)$ with the following norm

$$\|g\|_{L_Y^p(\mathbb{R}^d)} := \frac{1}{I(\mathcal{Y})} \sup_{1 \le v \le m} \left[ \|g\|_{L^{\infty}(\mathbb{R}^d)} \|Y_v\|_{L^{\infty,p}(\mathbb{R}^d)} + \|\mathrm{div}(gY_v)\|_{L^p(\mathbb{R}^d)} \right],$$

which is equivalent to the norm with $\mathrm{div}(gY_v)$ above replaced by $\partial_{Y_v} g$.

We now continue with the assumption of the initial density $\rho_0$ given in (1.2)-(1.3) associated with the compressible Navier-Stokes equations (1.1). We assume further that

$$(1.4) \qquad \rho_0 \in L_{X_0}^p(\mathbb{R}^d),$$

where $X_0 = (X_{0,1}, \dots, X_{0,m}) \subset L^{\infty,p}(\mathbb{R}^d)$ is a given nondegenerate family of vector fields for some $m \ge d-1$ and $p \in (d, \infty)$.
1.1. **Striated regularity.** We assume further striated regularity with respect to a given nondegenerate family of vector fields for the initial density $\rho_0$ in this paragraph.

We first introduce some notations, based on [10]. For some $p \in (d, \infty)$, $L^{\infty,p}(\mathbb{R}^d; \mathbb{R}^d)$ denotes the vector space of bounded vector fields with gradients in $L^p(\mathbb{R}^d; \mathbb{R}^{d \times d})$. From now on we denote the Lebesgue spaces $L^p(\mathbb{R}^d; \mathbb{R}^n)$ resp. Sobolev spaces $H^s(\mathbb{R}^d; \mathbb{R}^n)$ with $p \in [1, \infty]$, $s \in \mathbb{R}$ and $n \in \mathbb{N}^*$, simply by $L^p(\mathbb{R}^d)$ or $L^q$, resp. $H^s(\mathbb{R}^d)$ or $H^q$ with an abuse of notations. We have defined

$$L^{\infty,p}(\mathbb{R}^d) = \left\{ Y \in L^{\infty}(\mathbb{R}^d) \mid \|Y\|_{L^{\infty,p}(\mathbb{R}^d)} := \|Y\|_{L^{\infty}(\mathbb{R}^d)} + \|\nabla Y\|_{L^p(\mathbb{R}^d)} < \infty \right\}.$$

We define the norm $\| \cdot \|_{L^{\infty,p}}$ for a family of vector fields $\mathcal{Y} = (Y_1, Y_2, \dots, Y_m) \subset L^{\infty,p}(\mathbb{R}^d)$, $m \in \mathbb{N}$ as

$$\|\mathcal{Y}\|_{L^{\infty,p}(\mathbb{R}^d)} := \sup_{1 \le v \le m} \|Y_v\|_{L^{\infty,p}(\mathbb{R}^d)}.$$

**Definition 1.1** (Nondegeneracy). Let $\mathcal{Y} = (Y_1, Y_2, \dots, Y_m) \subset L^{\infty,p}(\mathbb{R}^d)$ be a family of $m$ vector fields with $m \ge d-1$ and $p \in (d, \infty)$. We say that $\mathcal{Y}$ is nondegenerate if it satisfies the following property

$$I(\mathcal{Y}) := \inf_{x \in \mathbb{R}^d} \sup_{\Upsilon \in \Upsilon_{d-1}^m} \left| \bigwedge^{d-1} Y_{\Upsilon}(x) \right|^{\frac{1}{d-1}} > 0.$$

Above $\Upsilon \in \Upsilon_{d-1}^m$ means that $\Upsilon = (v_1, v_2, \dots, v_{d-1})$ with each $v_i \in \{1, \dots, m\}$ and $v_i < v_j$ for $i < j$, $Y_{\Upsilon} := (Y_{v_1}, Y_{v_2}, \dots, Y_{v_{d-1}})$, while the symbol $\bigwedge^{d-1} Y_{\Upsilon}$ stands for the unique element of $\mathbb{R}^d$ such that

$$(\bigwedge^{d-1} Y_{\Upsilon}) \cdot Z = \det(Y_{v_1}, Y_{v_2}, \dots, Y_{v_{d-1}}, Z), \quad \forall Z \in \mathbb{R}^d.$$

**Definition 1.2** (Striated regularity with respect to a nondegenerate family of vector fields). Let $Y \in L^{\infty,p}(\mathbb{R}^d)$, $p \in (d, \infty)$ be a (single) vector field, and $\mathcal{Y} = (Y_1, Y_2, \dots, Y_m) \subset L^{\infty,p}(\mathbb{R}^d)$ be a nondegenerate family of vector fields with $m \ge d-1$.

(a) A function $g \in L^{\infty}(\mathbb{R}^d)$ is said to be of class $L^p(\mathbb{R}^d)$ along $Y$, if

$$g \in L_Y^p(\mathbb{R}^d) := \{h \in L^{\infty}(\mathbb{R}^d) \mid \operatorname{div}(hY) \in L^p(\mathbb{R}^d)\}.$$

We define the derivative of the function $g$ along $Y$ as follows

$$\partial_Y g := \operatorname{div}(gY) - g \operatorname{div} Y,$$

and hence we can equivalently define

$$L_Y^p(\mathbb{R}^d) = \{h \in L^{\infty}(\mathbb{R}^d) \mid \partial_Y h \in L^p(\mathbb{R}^d)\}.$$

(b) A function $g \in L^{\infty}(\mathbb{R}^d)$ is said to be of class $L^p(\mathbb{R}^d)$ along the family $\mathcal{Y}$, if

$$g \in L_Y^p(\mathbb{R}^d) := \bigcap_{1 \le v \le m} L_{Y_v}^p(\mathbb{R}^d),$$

and we equip the space $L_Y^p(\mathbb{R}^d)$ with the following norm

$$\|g\|_{L_Y^p(\mathbb{R}^d)} := \frac{1}{I(\mathcal{Y})} \sup_{1 \le v \le m} \left[ \|g\|_{L^{\infty}(\mathbb{R}^d)} \|Y_v\|_{L^{\infty,p}(\mathbb{R}^d)} + \|\operatorname{div}(gY_v)\|_{L^p(\mathbb{R}^d)} \right],$$

which is equivalent to the norm with $\operatorname{div}(gY_v)$ above replaced by $\partial_{Y_v} g$.

We now continue with the assumption of the initial density $\rho_0$ given in (1.2)-(1.3) associated with the compressible Navier-Stokes equations (1.1). We assume further that

$$(1.4) \qquad \rho_0 \in L_{X_0}^p(\mathbb{R}^d),$$

where $X_0 = (X_{0,1}, \dots, X_{0,m}) \subset L^{\infty,p}(\mathbb{R}^d)$ is a given nondegenerate family of vector fields for some $m \ge d-1$ and $p \in (d, \infty)$.
**Remark 1.3** (Initial density of density-patch type). It's interesting to notice that the initial density of the form

$$ (1.5) \qquad \rho_0 = \alpha \mathbb{1}_{D_0} + \tilde{\rho} \mathbb{1}_{D_0^c}, \quad \alpha \ge 0, $$

satisfies the assumptions for $\rho_0$ in (1.3)-(1.4), if $D_0$ is a $W^{2,p}(\mathbb{R}^d)$ (with $p > d$) bounded, simply connected domain in $\mathbb{R}^d$. Indeed, (1.4) holds for a nondegenerate (divergence-free) family of vector fields $\mathcal{X}_0 = (X_{0,1}, \dots, X_{0,m}) \subset \mathbb{L}^{\infty,p}(\mathbb{R}^d)$ which is tangent to $\partial D_0^1$, and this means that the initial density given by (1.5) persists tangential regularity with respect to the boundary $\partial D_0$.

## 1.2. Statement of the main results. The purpose of this paper is threefold.

1. We establish the local-in-time well-posedness of the system (1.1) for positive density function with striated regularity, under the compatibility condition. We thus remove the smallness condition required on the density in Danchin, Fanelli, Paicu's paper [10].

2. These local-in-time solutions become global-in-time unique solutions of the Cauchy problem (1.1)-(1.2)-(1.3)-(1.4), if

* $d=2$, and the bulk viscosity coefficient is large enough $\nu \ge \nu_0$ with $\nu_0$ depending on the norms of the initial data given in (1.2)-(1.3). This result is inspired by the work by Danchin and Mucha [16].

* $d=3$, the initial energy is small and the bulk viscosity coefficient $\nu \ge \nu_0$ is large enough. Here although $\|\rho_0 - \tilde{\rho}\|_{L^2(\mathbb{R}^d)}$ is assumed to be small, $\rho_0$ may have large variation in $L^\infty(\mathbb{R}^d)$. This result supplements the local-in-time well-posedness work [10] with global-in-time well-posedness result and the work by Shibata and Zhang [46] with less regular initial data.

3. Additionally, by letting the bulk viscosity tend to infinity $\nu \to \infty$, we establish a singular limit toward the incompressible inhomogeneous model on the whole space, in the spirit of Danchin and Mucha's work [16] where the considered domain has finite measure.

### 1.2.1. Local-in-time well-posedness and continuation criterion.
We begin by providing the statement of the local-in-time result, which technically further assumes the strict positivity of the initial density function and the compatibility condition on the initial data.

**Theorem 1.4** (Local-in-time well-posedness and continuation criterion). We consider the Cauchy problem of the compressible Navier-Stokes equations (1.1) supplemented with the initial data (1.2) verifying (1.3) and (1.4). We further assume the strict positivity of the initial density and the compatibility condition as follows

$$ (1.6) \qquad 0 < \underline{\rho} \le \rho_0(x) \quad \text{and} \quad \mu \Delta u_0 + (\mu + \lambda) \nabla \text{div} \, u_0 - \nabla P(\rho_0) \in L^2(\mathbb{R}^d). $$

Then, there exists a time $T > 0$ and a unique solution $(\rho, u)$ to the Cauchy problem (1.1)-(1.2), satisfying the following properties:

1) Energy bounds: $u \in \mathcal{C}([0, T], H^1(\mathbb{R}^d))$, $\dot{u} \in \mathcal{C}([0, T], L^2(\mathbb{R}^d))$, $\sqrt{\sigma} \nabla \dot{u}, \sigma \ddot{u} \in L^\infty((0, T), L^2(\mathbb{R}^d))$, and $\nabla \dot{u}, \sqrt{\sigma} \ddot{u}, \sigma \nabla \ddot{u} \in L^2((0, T) \times \mathbb{R}^d)$.

<footnote n="1">¹Indeed, for $d=2$ the existence of such a nondegenerate family of tangential vector fields is obvious since we can take $X_{0,1} = \begin{pmatrix} \partial_{x_2} f \\ -\partial_{x_1} f \end{pmatrix} =: \nabla^\perp f$ to be the tangent vector field close to $\partial D_0$ with $f \in W^{2,p}(\mathbb{R}^2)$ and $f|_{\partial D_0} = 0$ and $\nabla f|_{\partial D_0} \ne 0$, while $X_{0,2} = \nabla^\perp(\chi x_1)$ to be a non-zero vector field with $\chi$ a smooth cutoff function away from the boundary, see e.g. [43, (1.10)] (with $m=3$). The existence result for $d=3$ with $m=5$ follows from the similar idea, see e.g. [24, Proposition 3.2], where $X_{0,1} = \begin{pmatrix} 0 \\ -\partial_{x_3} f \\ \partial_{x_2} f \end{pmatrix}$, $X_{0,2} = \begin{pmatrix} \partial_{x_3} f \\ 0 \\ -\partial_{x_1} f \end{pmatrix}$, $X_{0,3} = \begin{pmatrix} -\partial_{x_2} f \\ \partial_{x_1} f \\ 0 \end{pmatrix}$ are generated by the function $f \in W^{2,p}(\mathbb{R}^3)$ with $f|_{\partial D_0} = 0$ and $\nabla f|_{\partial D_0} \ne 0$, while $X_{0,4} = \begin{pmatrix} \partial_{x_3}(\chi x_3) \\ 0 \\ -\partial_{x_1}(\chi x_3) \end{pmatrix}$, $X_{0,5} = \begin{pmatrix} -\partial_{x_2}(\chi x_1) \\ \partial_{x_1}(\chi x_1) \\ 0 \end{pmatrix}$ forms a nondegenerate family away from $\partial D_0$ with $\chi$ a smooth cutoff function away from the boundary $\partial D_0$.</footnote>
**Remark 1.3** (Initial density of density-patch type). It's interesting to notice that the initial density of the form

$$ (1.5) \qquad \rho_0 = \alpha 1_{D_0} + \tilde{\rho} 1_{D_0^c}, \quad \alpha \ge 0, $$

satisfies the assumptions for $\rho_0$ in (1.3)-(1.4), if $D_0$ is a $W^{2,p}(\mathbb{R}^d)$ (with $p > d$) bounded, simply connected domain in $\mathbb{R}^d$. Indeed, (1.4) holds for a nondegenerate (divergence-free) family of vector fields $\mathcal{X}_0 = (X_{0,1}, \dots, X_{0,m}) \subset \mathbb{L}^{\infty,p}(\mathbb{R}^d)$ which is tangent to $\partial D_0$<fnref n="1" />, and this means that the initial density given by (1.5) possesses tangential regularity with respect to the boundary $\partial D_0$.

## 1.2. Statement of the main results. The purpose of this paper is threefold.

1. We establish the local-in-time well-posedness of the system (1.1) for positive density function with striated regularity, under some compatibility condition. We thus remove the smallness condition required on the density fluctuation in Danchin, Fanelli, Paicu's paper [10].

2. These local-in-time solutions become global-in-time unique solutions of the Cauchy problem (1.1)-(1.10)-(1.2)-(1.3)-(1.4), if

* $d = 2$, and the bulk viscosity coefficient is large enough $\nu \ge \nu_0$ with $\nu_0$ depending on the norms of the initial data given in (1.2)-(1.3). This result is inspired by the work by Danchin and Mucha [16].

* $d = 3$, the initial energy is small and the bulk viscosity coefficient $\nu \ge \nu_0$ is large enough. Here although $\|\rho_0 - \tilde{\rho}\|_{L^2(\mathbb{R}^d)}$ is assumed to be small, $\rho_0$ may have large variation in $L^\infty(\mathbb{R}^d)$. This result supplements the local-in-time well-posedness work [10] with global-in-time well-posedness result and the work by Shibata and Zhang [46] with less regular initial data.

3. Additionally, by letting the bulk viscosity tend to infinity $\nu \to \infty$, we establish a singular limit toward the incompressible inhomogeneous model on the whole space, in the spirit of Danchin and Mucha's work [16] where the considered domain has finite measure.

### 1.2.1. Local-in-time well-posedness and continuation criterion.
We begin by providing the statement of the local-in-time result, which technically further assumes the strict positivity of the initial density function and the compatibility condition on the initial data.

**Theorem 1.4** (Local-in-time well-posedness and continuation criterion). We consider the Cauchy problem of the compressible Navier-Stokes equations (1.1) supplemented with the initial data (1.2) verifying (1.3) and (1.4). We further assume the strict positivity of the initial density and the compatibility condition as follows

$$ (1.6) \qquad 0 < \underline{\rho} \le \rho_0(x) \quad \text{and} \quad \mu \Delta u_0 + (\mu + \lambda) \nabla \text{div} \, u_0 - \nabla P(\rho_0) \in L^2(\mathbb{R}^d). $$

Then, there exists a time $T > 0$ and a unique solution $(\rho, u)$ to the Cauchy problem (1.1)-(1.2), satisfying the following properties:

1) Energy bounds: $u \in \mathcal{C}([0, T], H^1(\mathbb{R}^d))$, $\dot{u} \in \mathcal{C}([0, T], L^2(\mathbb{R}^d))$, $\sqrt{\sigma} \nabla \dot{u}, \sigma \ddot{u} \in L^\infty((0, T), L^2(\mathbb{R}^d))$, and $\nabla \dot{u}, \sqrt{\sigma} \ddot{u}, \sigma \nabla \ddot{u} \in L^2((0, T) \times \mathbb{R}^d)$.

<footnote n="1">¹Indeed, for $d = 2$ the existence of such a nondegenerate family of tangential vector fields is obvious since we can take $X_{0,1} = \begin{pmatrix} \partial_{x_2} f \\ -\partial_{x_1} f \end{pmatrix} =: \nabla^\perp f$ to be the tangent vector field close to $\partial D_0$ with $f \in W^{2,p}(\mathbb{R}^2)$ and $f|_{\partial D_0} = 0$ and $\nabla f|_{\partial D_0} \ne 0$, while $X_{0,2} = \nabla^\perp(\chi x_1)$ to be a non-zero vector field with $\chi$ a smooth cutoff function away from the boundary, see e.g. [43, (1.10)] (with $m = 3$). The existence result for $d = 3$ with $m = 5$ follows from the similar idea, see e.g. [24, Proposition 3.2], where $X_{0,1} = \begin{pmatrix} 0 \\ -\partial_{x_3} f \\ \partial_{x_2} f \end{pmatrix}$, $X_{0,2} = \begin{pmatrix} \partial_{x_3} f \\ 0 \\ -\partial_{x_1} f \end{pmatrix}$, $X_{0,3} = \begin{pmatrix} -\partial_{x_2} f \\ \partial_{x_1} f \\ 0 \end{pmatrix}$ are generated by the function $f \in W^{2,p}(\mathbb{R}^3)$ with $f|_{\partial D_0} = 0$ and $\nabla f|_{\partial D_0} \ne 0$, while $X_{0,4} = \begin{pmatrix} \partial_{x_3}(\chi x_3) \\ 0 \\ -\partial_{x_1}(\chi x_3) \end{pmatrix}$, $X_{0,5} = \begin{pmatrix} -\partial_{x_2}(\chi x_1) \\ \partial_{x_1}(\chi x_1) \\ 0 \end{pmatrix}$ forms a nondegenerate family away from $\partial D_0$ with $\chi$ a smooth cutoff function away from the boundary $\partial D_0$.</footnote>
Here and in what follows we use the notations

$$
(1.7) \qquad \sigma = \sigma(t) := \min\{1, t\}, \quad \dot{v} := (\partial_t + u \cdot \nabla)v, \quad \ddot{v} := (\partial_t + u \cdot \nabla)\dot{v}.
$$

2) **Striated regularity:** For all $0 < t < T$, we have $\rho(t) \in \mathbb{L}_{\mathcal{X}(t)}^p(\mathbb{R}^d)$, where $\mathcal{X}(t) = (X_v(t))_{1\le v\le m} \subset \mathbb{L}^{\infty,p}(\mathbb{R}^d)$ is the nondegenerate family of vector fields transported by the fluid flow, in the sense that each vector field $X_v(t)$, $1 \le v \le m$ solves uniquely the following Cauchy problem

$$
(1.8) \qquad \begin{cases} \partial_t X_v + u \cdot \nabla X_v = \partial_{X_v} u, \\ X_v|_{t=0} = X_{0,v}. \end{cases}
$$

Here the directional derivative was defined in **Definition 1.2**: $\partial_{X_v} u^j = \operatorname{div}(u^j X_v) - u^j \operatorname{div} X_v$, $1 \le j \le d$.

The velocity field is Lipschitz continuous when integrated in time and persists also the striated regularity for positive times:

• For $d=2$ or for $d=3$ and $3<p\le 6$, $\nabla u \in L^2((0,T), \mathbb{L}_{\mathcal{X}}^p(\mathbb{R}^d))$;

• For $d=3$ and $6<p<\infty$, $\sigma^{\frac{3}{4}-\frac{1}{r}-\frac{3}{2p}}\nabla u \in L^r((0,T), \mathbb{L}_{\mathcal{X}}^p(\mathbb{R}^3))$ and $\sigma^{\frac{3}{4}-\frac{1}{r}}\nabla \dot{u} \in L^r((0,T), L^3(\mathbb{R}^3))$, for any $2 \le r \le \infty$.

3) Continuation criterion: If $(\rho, u)$ is the solution defined up to a maximal time $T^* > 0$ and

$$
(1.9) \qquad \limsup_{t \to T^*} \left\{ \|\mathcal{X}(t)\|_{L^{\infty,p}(\mathbb{R}^d)} + \frac{1}{I(\mathcal{X}(t))} + \left\| \frac{1}{\rho(t)} \right\|_{L^{\infty}(\mathbb{R}^d)} + \|\rho(t)\|_{L^{\infty}(\mathbb{R}^d)} \right. \\
\qquad \left. + \|\partial_{\mathcal{X}(t)}\rho(t)\|_{L^p(\mathbb{R}^d)} + \|\nabla u(t)\|_{L^2(\mathbb{R}^d)} + \|\dot{u}(t)\|_{L^2(\mathbb{R}^d)} \right\} < \infty,
$$

then $T^* = \infty$.

The solution of **Theorem 1.4** is constructed in the spirit of a recent contribution of the second author [50], which deals with the more involved case of density-dependent viscosity coefficient. Thus, we only present a sketch of the proof of **Theorem 1.4** in **Appendix B**.

**Remark 1.5.** (a) This result supplements the contribution by Danchin, Fanelli and Paicu [10] by removing the smallness condition on the density deviation. Unlike the maximum regularity argument used in [10], which requires a critical regularity for one part of the initial velocity, our method relies on the change into Lagrangian coordinates along with energy estimation methods.

(b) The compatibility condition in (1.6), $\mu\Delta u_0 + (\mu + \lambda)\nabla \operatorname{div} u_0 - \nabla P(\rho_0) \in L^2(\mathbb{R}^d)$ expresses the continuity of the stress tensor, and does not require (explicitly) smoothness of the density. The parabolic effect of the momentum equations ensures that this condition holds true at positive times even for less regular initial data, see [27].

(c) The velocity field persists further regularity property which are stated in **Corollary B.2** below, thanks to the decomposition of the velocity gradient (**B.13**).

1.2.2. Definitions of energy functionals. The global-in-time well-posedness result will follow from the above local-in-time well-posedness, continuation criterion and a series of energy estimates. We define in this subsection the relevant energy functionals.

Recall the positive density equilibrium $\tilde{\rho}$ in (1.3). We define first

$$
(1.10) \qquad \text{The pressure equilibrium } \tilde{P} := P(\tilde{\rho}),
$$

$$
(1.11) \qquad \text{The } \rho\text{-dependent functions } H_l(\rho) = \rho \int_{\tilde{\rho}}^{\rho} s^{-2} |P(s) - \tilde{P}|^{l-1} (P(s) - \tilde{P}) ds, \quad l \in [1, \infty),
$$

$$
(1.12) \quad \text{The pressure deviation } G(t, x) := (P(\rho))(t, x) - \tilde{P},
$$

$$
(1.13) \quad \text{The effective flux } F(t, x) = \nu(\operatorname{div} u)(t, x) - G(t, x).
$$

We define the associated energy function of the compressible Navier-Stokes equations (1.1)

$$
(1.14) \quad E(t) = \int_{\mathbb{R}^d} \left[ \rho \frac{|u|^2}{2} + H_1(\rho) \right] (t, x) dx + \int_0^t \left[ \mu \|\nabla u(t')\|_{L^2(\mathbb{R}^d)}^2 + (\mu + \lambda) \|\operatorname{div} u(t')\|_{L^2(\mathbb{R}^d)}^2 \right] dt',
$$
Here and in what follows we use the notations

$$
(1.7) \qquad \sigma = \sigma(t) := \min\{1, t\}, \quad \dot{v} := (\partial_t + u \cdot \nabla)v, \quad \ddot{v} := (\partial_t + u \cdot \nabla)\dot{v}.
$$

2) Striated regularity: For all $0 < t < T$, we have $\rho(t) \in \mathbb{L}_{\mathcal{X}(t)}^p(\mathbb{R}^d)$, where $\mathcal{X}(t) = (X_v(t))_{1\le v\le m} \subset \mathbb{L}^{\infty,p}(\mathbb{R}^d)$ is the nondegenerate family of vector fields transported by the fluid flow, in the sense that each vector field $X_v(t)$, $1 \le v \le m$ solves uniquely the following Cauchy problem

$$
(1.8) \qquad \begin{cases} \partial_t X_v + u \cdot \nabla X_v = \partial_{X_v} u, \\ X_v|_{t=0} = X_{0,v}. \end{cases}
$$

Here the directional derivative was defined in *Definition 1.2*: $\partial_{X_v} u^j = \text{div}(u^j X_v) - u^j \text{div} X_v$, $1 \le j \le d$.

The velocity field is Lipschitz continuous when integrated in time: $\nabla u \in L^1((0,T), L^\infty(\mathbb{R}^d))$, and enjoys further the striated regularity for positive times:

• For $d=2$ or for $d=3$ and $3<p\le 6$, $\nabla u \in L^2((0,T), \mathbb{L}_{\mathcal{X}}^p(\mathbb{R}^d));$

• For $d=3$ and $6<p<\infty$, $\sigma^{\frac{3}{4}-\frac{1}{r}-\frac{3}{2p}}\nabla u \in L^r((0,T), \mathbb{L}_{\mathcal{X}}^p(\mathbb{R}^3))$ and $\sigma^{\frac{3}{4}-\frac{1}{r}}\nabla \dot{u} \in L^r((0,T), L^3(\mathbb{R}^3))$, for any $2 \le r \le \infty$.

3) Continuation criterion: If $(\rho, u)$ is the solution defined up to a maximal time $T^* > 0$ with $T^* < \infty$, then

$$
(1.9) \qquad \limsup_{t \to T^*} \left\{ \begin{aligned}
& \|\mathcal{X}(t)\|_{L^{\infty,p}(\mathbb{R}^d)} + \frac{1}{I(\mathcal{X}(t))} + \left\| \frac{1}{\rho(t)} \right\|_{L^{\infty}(\mathbb{R}^d)} + \|\rho(t)\|_{L^{\infty}(\mathbb{R}^d)} \\
& + \|\partial_{\mathcal{X}(t)}\rho(t)\|_{L^p(\mathbb{R}^d)} + \|\nabla u(t)\|_{L^2(\mathbb{R}^d)} + \|\dot{u}(t)\|_{L^2(\mathbb{R}^d)}
\end{aligned} \right\} = \infty.
$$

The solution of **Theorem 1.4** is constructed in the spirit of a recent contribution of the second author [49], which deals with the more involved case of density-dependent viscosity coefficient. Thus, we only present a sketch of the proof of **Theorem 1.4** in Appendix B.

While we do not pursue optimal local-in-time well-posedness in this paper, we instead employ the approximation argument to the local theory established above to prove our main result concerning global-in-time well-posedness in **Theorem 1.6** below. We just mention here that the strict positivity requirement in assumption (1.6) can be relaxed through standard approximation techniques, with the estimates and results being corrected correspondingly (e.g. with density weights as in **Proposition B.1**).

**Remark 1.5.** (a) This result supplements the contribution by Danchin, Fanelli and Paicu [10] by removing the smallness condition on the density deviation. Unlike the maximum regularity argument used in [10], which requires a critical regularity for one part of the initial velocity, our method relies on the change into Lagrangian coordinates along with energy estimation methods.

(b) The compatibility condition in (1.6), $\mu\Delta u_0 + (\mu + \lambda)\nabla \text{div} u_0 - \nabla P(\rho_0) \in L^2(\mathbb{R}^d)$ expresses the continuity of the stress tensor, and does not require (explicitly) smoothness of the density. The parabolic effect of the momentum equations ensures that this condition holds true at positive times even for less regular initial data, see [27].

(c) The velocity field possesses indeed further regularity properties which are stated in **Corollary B.2** below, thanks to the decomposition of the velocity gradient (**B.13**).

1.2.2. *Definitions of energy functionals.* The global-in-time well-posedness result will follow from the above local-in-time well-posedness, continuation criterion and a series of energy estimates. We define in this subsection the relevant energy functionals. In the following of this paper we assume the general non-decreasing pressure law $P = P(\rho)$ as follows:

$$
(1.10) \qquad P \in \mathcal{C}^2(\mathbb{R}, \mathbb{R}) \text{ and } P'(\rho) \ge 0 \text{ for } \rho \ge 0.
$$
which consists of

• The kinetic energy $\frac{1}{2}\|\sqrt{\rho}u(t)\|_{L^2(\mathbb{R}^d)}^2$;

• The potential energy $\int_{\mathbb{R}^d} (H_1(\rho))(t, x)dx$ with $H_1(\rho)$ defined in (1.11): $H_1(\rho) = \rho \int_{\tilde{\rho}}^{\rho} \frac{P(s) - \tilde{P}}{s^2}ds$;

• The energy dissipation $\mu\|\nabla u\|_{L^2((0,t)\times\mathbb{R}^d)}^2 + (\mu + \lambda)\|\text{div } u\|_{L^2((0,t)\times\mathbb{R}^d)}^2$.

The energy $E(t)$ is conserved for regular enough solutions of (1.1).

Recall the notations in (1.7), and we introduce energy functionals of higher order

$$ (1.15) \qquad \begin{aligned} A_1(t) &= \frac{\mu}{2} \|\nabla u(t)\|_{L^2(\mathbb{R}^d)}^2 + \frac{\mu + \lambda}{2} \|\text{div } u(t)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \|\sqrt{\rho} \dot{u}(t')\|_{L^2(\mathbb{R}^d)}^2 dt', \\ A_2(t) &= \sigma(t) \|\sqrt{\rho} \dot{u}(t)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \sigma(t') \left[ \mu \|\nabla \dot{u}(t')\|_{L^2(\mathbb{R}^d)}^2 + \frac{\mu + \lambda}{\nu^2} \|\dot{F}(t')\|_{L^2(\mathbb{R}^d)}^2 \right] dt'. \end{aligned} $$

The hierarchy of energy functionals $E(t)$, $A_1(t)$, $A_2(t)$ encode $L^2(\mathbb{R}^d)$-norm, $\dot{H}^1(\mathbb{R}^d)$-norm for $u(t)$ and (time-weighted) $L^2(\mathbb{R}^d)$-norm for the material derivative $\dot{u}(t)$, respectively. Although trivially $|\text{div } u| \le d|\nabla u|$, we will make efforts to get the (large) viscosity coefficient $\lambda$ before $\text{div } u$ in the definition of $E$, $A_1$, such that intuitively $\text{div } u \to 0$ as $\lambda \to \infty$ if $E, A_1$ is bounded uniformly in time. The review of their history can be found in Section 1.3 below.

Recall the initial data (1.2) and we denote $G_0(x) = G(0, x) = (P(\rho_0))(x) - \tilde{P}$. For notational simplicity we denote

$$ (1.16) \qquad \text{The first initial energy } E_0 := E(0) = \int_{\mathbb{R}^d} \left[ \rho_0 \frac{|u_0|^2}{2} + H_1(\rho_0) \right] (x) dx, $$

$$ (1.17) \qquad \text{The total initial energy } E_0^\nu := E_0 + \mu \|\nabla u_0\|_{L^2(\mathbb{R}^d)}^2 + \nu \|\text{div } u_0\|_{L^2(\mathbb{R}^d)}^2 + \frac{1}{\nu} \|G_0\|_{L^2(\mathbb{R}^d)}^2, $$

$$ (1.18) \qquad \text{The upper bound of the initial density } \rho_0^* := \sup_{x \in \mathbb{R}^d} \rho_0(x). $$

We observe that for initial data given in (1.2)-(1.3),

$$ E_0 \le C(\rho_0^*) \|(\rho_0 - \tilde{\rho}, u_0)\|_{L^2(\mathbb{R}^d)} < +\infty, \quad E_0^\nu \le C(\mu, \underline{\nu}, \rho_0^*)(E_0 + \|\nabla u_0\|_{L^2(\mathbb{R}^d)}^2) + \nu \|\text{div } u_0\|_{L^2(\mathbb{R}^d)}^2 < +\infty. $$

We aim to bound $A_1(t)$, $A_2(t)$ globally in time in terms of $E_0^\nu, \rho_0^*$, if $\nu \ge \nu_0$ is large enough (and if the initial energy is small enough for $d=3$). The following quantity $A_3(t)$ captures the striated regularity of the density function along the family of vector fields $\mathcal{X}(t) = (X_\nu(t))_{1\le\nu\le m}$ transported by the flow as in (1.8)

$$ (1.19) \qquad \mathcal{A}_3(t) = \|\mathcal{X}(t)\|_{L^\infty, p(\mathbb{R}^d)} + \sup_{1\le\nu\le m} \|(\partial_{X_\nu}\rho)(t)\|_{L^p(\mathbb{R}^d)}. $$

It is straightforward to see that $\mathcal{A}_3(t)$ grows exponentially in $\|\nabla u\|_{L^1(\mathbb{R}^d)}$. We aim to show that the striated regularity encoded in $\|\log \mathcal{A}_3\|_{L^1((0,t))}$ together with the energy functionals $A_1(t)$, $A_2(t)$ controls $\|\nabla u\|_{L^1(\mathbb{R}^d)}$. Gronwall's inequality hence implies the exponential-in-time control of $\|\nabla u\|_{L^1(\mathbb{R}^d)}$.

1.2.3. Global-in-time well-posedness. We now state our global-in-time result for less regular initial data on which the assumption (1.6) is not assumed.

**Theorem 1.6.** Assume the Cauchy problem (1.1)-(1.2)-(1.3)-(1.4) and the following conditions

$$ (1.20) \qquad \text{either } d=2 \text{ and } \nu \ge \nu_0, $$

$$ (1.21) \qquad \text{or } d=3, p \in (3,6), E_0^\nu E_0 \le c \text{ and } \nu \ge \nu_0, $$

where c is a fixed constant depending only on $\mu, \underline{\nu}$ while $\nu_0$ is a constant depending additionally on the initial norms: $E_0, \|\nabla u_0\|_{L^2(\mathbb{R}^d)}, \rho_0^*$. Then the Cauchy problem has a unique global-in-time solution $(\rho, u)$ verifying
Recall the positive density equilibrium $\tilde{\rho}$ in (1.3). We define first

$$
(1.11) \quad \text{The pressure equilibrium } \tilde{P} := P(\tilde{\rho}),
$$

$$
(1.12) \quad \text{The } \rho\text{-dependent functions } H_l(\rho) = \rho \int_{\tilde{\rho}}^{\rho} s^{-2} |P(s) - \tilde{P}|^{l-1} (P(s) - \tilde{P}) ds, \quad l \in [1, \infty),
$$

$$
(1.13) \quad \text{The pressure deviation } G(t, x) := (P(\rho))(t, x) - \tilde{P},
$$

$$
(1.14) \quad \text{The effective flux } F(t, x) = \nu(\operatorname{div} u)(t, x) - G(t, x).
$$

We remark that due to the monotonicity property of the pressure law $P = P(\rho)$ in (1.10), the function $H_l(\rho)$ is always nonnegative for nonnegative $\rho$.

We define the associated energy function of the compressible Navier-Stokes equations (1.1)-(1.10)

$$
(1.15) \quad E(t) = \int_{\mathbb{R}^d} \left[ \rho \frac{|u|^2}{2} + H_1(\rho) \right] (t, x) dx + \int_0^t \left[ \mu \|\nabla u(t')\|_{L^2(\mathbb{R}^d)}^2 + (\mu + \lambda) \|\operatorname{div} u(t')\|_{L^2(\mathbb{R}^d)}^2 \right] dt',
$$

which consists of

• The kinetic energy $\frac{1}{2} \|\sqrt{\rho}u(t)\|_{L^2(\mathbb{R}^d)}^2$;

• The potential energy $\int_{\mathbb{R}^d} (H_1(\rho))(t, x)dx$ with $H_1(\rho)$ defined in (1.12): $H_1(\rho) = \rho \int_{\tilde{\rho}}^{\rho} \frac{P(s) - \tilde{P}}{s^2} ds$;

• The energy dissipation $\mu \|\nabla u\|_{L^2((0,t)\times\mathbb{R}^d)}^2 + (\mu + \lambda) \|\operatorname{div} u\|_{L^2((0,t)\times\mathbb{R}^d)}^2$.

The energy $E(t)$ is conserved for regular enough solutions of (1.1).

Recall the notations in (1.7), and we introduce energy functionals of higher order

$$
\begin{equation} \tag{1.16}
\begin{aligned}
\mathcal{A}_1(t) &= \frac{\mu}{2} \|\nabla u(t)\|_{L^2(\mathbb{R}^d)}^2 + \frac{\mu + \lambda}{2} \|\operatorname{div} u(t)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \|\sqrt{\rho} \dot{u}(t')\|_{L^2(\mathbb{R}^d)}^2 dt', \\
\mathcal{A}_2(t) &= \sigma(t) \|\sqrt{\rho} \dot{u}(t)\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \sigma(t') \left[ \mu \|\nabla \dot{u}(t')\|_{L^2(\mathbb{R}^d)}^2 + \frac{\mu + \lambda}{\nu^2} \|\dot{F}(t')\|_{L^2(\mathbb{R}^d)}^2 \right] dt'.
\end{aligned}
\end{equation}
$$

The hierarchy of energy functionals $E(t)$, $\mathcal{A}_1(t)$, $\mathcal{A}_2(t)$ encode $L^2(\mathbb{R}^d)$-norm, $\dot{H}^1(\mathbb{R}^d)$-norm for $u(t)$ and (time-weighted) $L^2(\mathbb{R}^d)$-norm for the material derivative $\dot{u}(t)$, respectively. Although trivially $|\operatorname{div} u| \le d|\nabla u|$, we will make efforts to get the (large) viscosity coefficient $\lambda$ before $\operatorname{div} u$ in the definition of $E$, $\mathcal{A}_1$, such that intuitively $\operatorname{div} u \to 0$ as $\lambda \to \infty$ if $E$, $\mathcal{A}_1$ are bounded uniformly in time. A review of the history of these energy functionals can be found in Section 1.3 below.

Recall the initial data (1.2) and we denote $G_0(x) = G(0, x) = (P(\rho_0))(x) - \tilde{P}$. For notational simplicity we denote

$$
(1.17) \quad \text{The first initial energy } E_0 := E(0) = \int_{\mathbb{R}^d} \left[ \rho_0 \frac{|u_0|^2}{2} + H_1(\rho_0) \right] (x) \,dx,
$$

$$
(1.18) \quad \text{The total initial energy } E_0^\nu := E_0 + \mu \|\nabla u_0\|_{L^2(\mathbb{R}^d)}^2 + \nu \|\operatorname{div} u_0\|_{L^2(\mathbb{R}^d)}^2 + \frac{1}{\nu} \|G_0\|_{L^2(\mathbb{R}^d)}^2,
$$

$$
(1.19) \quad \text{The upper bound of the initial density } \rho_0^* := \sup_{x \in \mathbb{R}^d} \rho_0(x).
$$

We observe that for initial data given in (1.2)-(1.3),

$$
E_0 \le C(\rho_0^*) \|(\rho_0 - \tilde{\rho}, u_0)\|_{L^2(\mathbb{R}^d)} < +\infty, \quad E_0^\nu \le C(\mu, \nu, \rho_0^*)(E_0 + \|\nabla u_0\|_{L^2(\mathbb{R}^d)}^2 + \nu \|\operatorname{div} u_0\|_{L^2(\mathbb{R}^d)}^2) < +\infty.
$$

We aim to bound $\mathcal{A}_1(t)$, $\mathcal{A}_2(t)$ globally in time in terms of $E_0^\nu, \rho_0^*$, if $\nu \ge \nu_0$ is large enough (and if the initial energy is small enough for $d=3$). The following quantity $\mathcal{A}_3(t)$ captures the striated regularity of the density function along the family of vector fields $\mathcal{X}(t) = (X_\nu(t))_{1\le\nu\le m}$ transported by the flow as in (1.8)

$$
(1.20) \qquad \mathcal{A}_3(t) = \|\mathcal{X}(t)\|_{\ell^{\infty,p}(\mathbb{R}^d)} + \sup_{1 \le \nu \le m} \|(\partial_{X_\nu} \rho)(t)\|_{L^p(\mathbb{R}^d)}.
$$
1) Energy bounds: For all $t \ge 0$, we have

$$ (1.22) \qquad \begin{cases} E(t) + A_1(t) + A_2(t) \le C_0^\nu, \\ \|\rho(t) - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)}^2 \le \|\rho_0 - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)}^2 + C_0^\nu, \end{cases} \quad \text{where the constant } C_0^\nu \text{ depends on } \mu, \underline{\nu}, \rho_0^*, \text{ and (superlinearly) on } E_0^\nu. $$

2) Striated regularity: For all $t \ge 0$, $\rho(t) \in L_{\mathcal{X}(t)}^p(\mathbb{R}^d)$, where $\mathcal{X}(t) = (X_\nu(t))_{1\le v \le m} \subset L^{\infty,p}(\mathbb{R}^d)$ is a nondegenerate family of vector fields defined to solve the Cauchy problem (1.8).

Moreover $\nabla u \in L_{\text{loc}}^1([0, \infty), L^\infty(\mathbb{R}^d))$ with the following estimates:

$$ (1.23) \qquad \begin{cases} A_3(t) \le A_3(0) \exp \left( C_0 \int_0^t \left[ 1 + \sqrt{t} + \|\nabla u(t')\|_{L^\infty(\mathbb{R}^d)} dt' \right] \right), \int_0^t \|\nabla u(t')\|_{L^\infty(\mathbb{R}^d)} dt' \le C_0 \left(1 + \frac{A_3(0)}{I(\mathcal{X}_0)}\right) \exp(C_0 t), \\ \int_0^t \|\partial_{\mathcal{X}(t')} \nabla u(t')\|_{L^p(\mathbb{R}^d)} dt' \le C_0 \left(1 + t + t A_3(t)\right) A_3(t), \end{cases} \quad \text{where } C_0 \text{ depends on } \mu, \underline{\nu}, m, d, p, \rho_0^*, E_0^\nu. $$

**Remark 1.7** (Bounds for $\text{div } u$). We have assumed some uniform bounds (with respect to $\nu$) for $\text{div } u_0$ implicitly: The conditions in (1.20) and (1.21) imply that

$$ \nu \|\text{div } u_0\|_{L^2(\mathbb{R}^d)}^2 \le \begin{cases} E_0^\nu < \infty, & \text{if } d=2, \\ E_0^\nu \min \left\{ 1, \frac{c}{E_0} \right\} < \infty, & \text{if } d=3. \end{cases} $$

This boundedness is propagated over time:

$$ \nu \|\text{div } u(t)\|_{L^2(\mathbb{R}^d)}^2 \le C_0^\nu. $$

**Theorem 1.6** and **Remark 1.3** imply immediately

**Corollary 1.8** (Density patch problem in the exterior domain). The Cauchy problem (1.1)-(1.2) with initial density of density-patch type (1.5) and $u_0 \in H^1(\mathbb{R}^d)$, under the assumption (1.20) or (1.21), has a unique global-in-time solution $(\rho, u)$ with $\rho(t)$ persisting tangential regularity with respect to the boundary $\partial D_t$ which is transported by the flow of $u$ and keeps its $W^{2,p}(\mathbb{R}^d)$-regularity.

**Remark 1.9.** If $\alpha > 0$, deriving a uniformly positive lower bound for the density is straightforward (see [17, Section 2.1 (Step 6)]). This results in an exponential-in-time decay of the jump in the density $\rho(t)$ across $\partial D_t$, as observed in [29, 30].

Intuitively, thanks to the uniform bound in (1.22): $A_1(t) \le C_0^\nu$, letting $\nu \to \infty$ yields a couple $(\varrho, v)$ that satisfies the incompressible inhomogeneous model:

$$ (1.24) \qquad \left\{ \begin{aligned} \partial_t \varrho + \text{div}(\varrho v) &= 0, \\ \partial_t (\varrho v) + \text{div}(\varrho v \otimes v) + \nabla \Pi - \mu \Delta v &= 0, \\ \text{div } v &= 0. \end{aligned} \right. $$

**Corollary 1.10** (Incompressible limit). Let $(\rho_0, u_0)$ be the initial data given in (1.2) verifying (1.3), (1.4) and $\text{div } u_0 = 0$. Let $(\rho^{(\nu)}, u^{(\nu)})$ be the corresponding unique solution constructed in **Theorem 1.6**, under the assumption (1.20) or (1.21).

Then $(\rho^{(\nu)}, u^{(\nu)})_\nu$ converges weakly-$^*$ to $(\varrho, v)$ in $L^\infty((0, \infty) \times \mathbb{R}^d) \times L^\infty((0, \infty), H^1(\mathbb{R}^d))$ as $\nu$ goes to infinity, and $(\varrho, v)$ solves (uniquely) the inhomogeneous, incompressible model (1.24) with initial data $(\rho_0, u_0)$ in the distribution sense. Moreover, we have

$$ (1.25) \qquad \begin{cases} \text{div } u^{(\nu)} = \mathcal{O}(\nu^{-1/2}) & \text{in } L^2 \cap L^\infty((0, \infty), L^2(\mathbb{R}^d)), \\ \partial_t (\rho^{(\nu)} u^{(\nu)}) + \text{div}(\rho^{(\nu)} u^{(\nu)} \otimes u^{(\nu)}) - \nabla F^{(\nu)} - \mu \Delta u^{(\nu)} = \mathcal{O}(\nu^{-1/2}) & \text{in } L^\infty((0, \infty), \dot{H}^{-1}(\mathbb{R}^d)), \end{cases} \quad \text{where } F^{(\nu)} = \nu \text{div } u^{(\nu)} - G^{(\nu)} \text{ with } G^{(\nu)} = P(\rho^{(\nu)}) - \tilde{P}. $$
It is straightforward to see that $\mathcal{A}_3(t)$ grows exponentially in $\|\nabla u\|_{L_t^1 L^\infty}$. We aim to show that the striated regularity encoded in $\|\log \mathcal{A}_3\|_{L^1((0,t))}$ together with the energy functionals $\mathcal{A}_1(t)$, $\mathcal{A}_2(t)$ controls $\|\nabla u\|_{L_t^1 L^\infty}$. Grönwall's inequality hence implies the exponential-in-time control of $\|\nabla u\|_{L_t^1 L^\infty}$.

### 1.2.3. Global-in-time well-posedness.
We now state our global-in-time result for less regular initial data on which the assumption (1.6) is not assumed.

**Theorem 1.6.** Assume the Cauchy problem (1.1)-(1.10) with initial data (1.2)-(1.3)-(1.4) and the following conditions

$$ (1.21) \qquad \text{either } d = 2 \text{ and } \nu \ge \nu_0, $$

$$ (1.22) \qquad \text{or } d = 3, p \in (3,6), E_0^\nu E_0 \le c \text{ and } \nu \ge \nu_0, $$

where c is a fixed constant depending only on $\mu, \nu$ while $\nu_0$ is a constant depending additionally on the initial norms: $E_0, \|\nabla u_0\|_{L^2(\mathbb{R}^d)}$, $\rho_0^*$. Then the Cauchy problem has a unique global-in-time solution $(\rho, u)$ verifying

1) Energy bounds: For all $t \ge 0$, we have

$$ (1.23) \qquad \begin{cases} E(t) + \mathcal{A}_1(t) + \mathcal{A}_2(t) \le C_0^\nu, \\ \|\rho(t) - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)}^2 \le \|\rho_0 - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)}^2 + C_0^\nu, \end{cases} $$

where the constant $C_0^\nu$ depends on $\mu, \nu, \rho_0^*$, and (superlinearly) on $E_0^\nu$.

2) Striated regularity: For all $t \ge 0$, $\rho(t) \in \mathbb{L}^p_{\mathcal{X}(t)}(\mathbb{R}^d)$, where $\mathcal{X}(t) = (X_\nu(t))_{1 \le \nu \le m} \subset \mathbb{L}^{\infty,p}(\mathbb{R}^d)$ is a nondegenerate family of vector fields defined to solve the Cauchy problem (1.8).

Moreover $\nabla u \in L_{\text{loc}}^1([0, \infty), L^\infty(\mathbb{R}^d))$ with the following estimates:

$$ (1.24) \qquad \begin{cases} \mathcal{A}_3(t) \le \mathcal{A}_3(0) \exp \left( C_0 \int_0^t \left[ 1 + \sqrt{t} + \|\nabla u(t')\|_{L^\infty(\mathbb{R}^d)} dt' \right] \right), \\ \int_0^t \|\nabla u(t')\|_{L^\infty(\mathbb{R}^d)} dt' \le C_0 \left(1 + \frac{\mathcal{A}_3(0)}{I(\mathcal{X}_0)}\right) \exp(C_0 t), \\ \int_0^t \|\partial_{\mathcal{X}(t')} \nabla u(t')\|_{L^p(\mathbb{R}^d)} dt' \le C_0 \left(1 + t + t \mathcal{A}_3(t)\right) \mathcal{A}_3(t), \end{cases} $$

where $C_0$ depends on $\mu, \nu, m, d, p, \rho_0^*, E_0^\nu$.

**Remark 1.7** (Bounds for div $u$). We have assumed some uniform bounds (with respect to $\nu$) for div $u_0$ implicitly: The conditions in (1.21) and (1.22) imply that

$$ \nu \|\text{div } u_0\|_{L^2(\mathbb{R}^d)}^2 \le \begin{cases} E_0^\nu < \infty, & \text{if } d=2, \\ E_0^\nu \min \left\{ 1, \frac{c}{E_0} \right\} < \infty, & \text{if } d=3. \end{cases} $$

This boundedness is propagated over time:

$$ \nu \|\text{div } u(t)\|_{L^2(\mathbb{R}^d)}^2 \le C_0^\nu. $$

**Theorem 1.6** and **Remark 1.3** imply immediately

**Corollary 1.8** (Density patch problem in the exterior domain). The Cauchy problem (1.1)-(1.10) with initial data (1.2) of density-patch type (1.5) and $u_0 \in H^1(\mathbb{R}^d)$, under the assumption (1.21) or (1.22), has a unique global-in-time solution $(\rho, u)$ with $\rho(t)$ enjoying tangential regularity with respect to the boundary $\partial D_t$ which is transported by the flow of $u$ and keeps its $W^{2,p}(\mathbb{R}^d)$-boundary regularity.

**Remark 1.9.** If $\alpha > 0$, deriving a uniformly positive lower bound for the density is straightforward (see [17, Section 2.1 (Step 6)]). This results in an exponential-in-time decay of the jump in the density $\rho(t)$ across $\partial D_t$, as observed in [29, 30].
The proofs of Theorem 1.6 and Corollary 1.10 are presented in Section 2.3, based on the a priori estimates in Section 2.1 and their proofs in Section 2.2.

**Remark 1.11.** This result in Corollary 1.10 is a partial continuation of the work by Danchin and Mucha [13, 14, 16], and Danchin and Wang [17], and stands, as far as we know, as the first one dealing with discontinuous initial data in the whole space. We notice that, except for the work [13] dealing with the whole space case and initial data in the critical Besov space, the other studies rely heavily on the assumption that the domain has finite measure. The extension to the whole space, especially for $d = 2$, is not obvious, and it requires some refined computations, e.g. the compensated result by Coifman, Lions, Meyer, Semmes in [8].

**1.3. Review of known results.** Classical solutions for the Navier-Stokes equations (1.1) with regular initial data are known to exist, since the work by Nash [42], Itaya [33, 34], Solonnikov [47], Tani [48], just to cite a few examples. These solutions are defined up to a positive time which depends on the (norms of) initial data. The first result addressing the global-in-time well-posedness of classical solutions is provided by Matsumura and Nishida [41] for small initial data in $L^1(\mathbb{R}^3) \cap H^3(\mathbb{R}^3)$. Nowadays, global-in-time classical solutions are known to exist under smallness assumption on the initial data in critical Besov space [4, 7, 26].

**Weak solutions and estimates for $E(t)$ and $G$.** Similar to the solutions constructed by Leray [35] for the incompressible Navier-Stokes equations, there are well-established results that investigate the existence of global-in-time weak solutions for the compressible Navier-Stokes equations (1.1), with finite initial energy. The first result was obtained by P.-L. Lions [40], followed by Feireisl, Novotný, Petzeltová [23], for pressure laws of the form $P(\rho) = a\rho^\gamma$, $a > 0$, with some limitations on $\gamma$. These weak solutions verify the following classical energy inequality:

$$ (1.26) \qquad E(t) \le E(0) = E_0, $$

where the functional $E$ has been given in (1.15).

The introduction of the (generalized) specific energy $H_l(\rho)$, $l \in [1, \infty)$, in (1.11) helps (technically) to estimate the pressure deviation $G$. As observed in e.g. [3], the so-defined $H_l(\rho)$ is non-negative: $H_l(\rho) \ge 0$, since the pressure $P(\rho)$ is an increasing function of the density.

For the classical case $l = 1$, $H_1(\rho)$ appears in the definition of $E(t)$, which is integrable in space uniformly in time due to (1.26): $\int_{\mathbb{R}^d} (H_1(\rho))(t, x) dx \le E(t) = E_0$. Consequently, under the a priori assumption

$$ \rho(t, x) \le \rho^*, $$

we have the estimates for $G$ uniformly in time by the energy $E_0$ below

$$ (1.27) \qquad \sup_{[0,t]} \|G\|_{L^q(\mathbb{R}^d)}^q \le C^* \sup_{t' \in [0,t]} \int_{\mathbb{R}^d} H_1(\rho(t',x)) dx \le C^* E_0, \quad \text{with } q \in [2, \infty), $$

where the constant $C^*$ depends only on $\rho^*$ and $q$. In the above, the first inequality follows from the definition of $H_1(\rho)$ in (1.11).

General $H_l(\rho)$, $l \ge 1$, as a function of $\rho$, satisfies the following ordinary differential equation

$$ \rho H_l'(\rho) - H_l(\rho) = |P(\rho) - \tilde{P}|^{l-1} (P(\rho) - \tilde{P}), $$

and hence, by virtue of the mass equation (1.1)$_1$, the function $(H_l(\rho))(t, x)$ satisfies the following time evolutionary equation

$$ \partial_t H_l(\rho) + \operatorname{div}(H_l(\rho)u) + |P(\rho) - \tilde{P}|^{l-1} (P(\rho) - \tilde{P}) \operatorname{div} u = 0, $$

which is in the same spirit of the renormalized continuity equation appearing in e.g. [40]. By view of the definitions (1.12) and (1.13), it is equivalent to

$$ \partial_t H_l(\rho) + \operatorname{div}(H_l(\rho)u) + \frac{1}{\nu} |G|^{l+1} = -\frac{1}{\nu} |G|^{l-1} GF. $$
Intuitively, thanks to the uniform bound in (1.23): $A_1(t) \le C_0^\nu$, letting $\nu \to \infty$ yields a couple $(\varrho, v)$ that satisfies the incompressible inhomogeneous model:

$$ (1.25) \qquad \left\{ \begin{array}{l} \partial_t \varrho + \operatorname{div}(\varrho v) = 0, \\ \partial_t (\varrho v) + \operatorname{div}(\varrho v \otimes v) + \nabla \Pi - \mu \Delta v = 0, \\ \operatorname{div} v = 0. \end{array} \right. $$

**Corollary 1.10** (Incompressible limit). Let $(\rho_0, u_0)$ be the initial data given in (1.2) verifying (1.3), (1.4) and $\operatorname{div} u_0 = 0$. Let $(\rho^{(\nu)}, u^{(\nu)})$ be the corresponding unique solution constructed in Theorem 1.6, under the assumption (1.21) or (1.22).

Then $(\rho^{(\nu)}, u^{(\nu)})_\nu$ converges weakly-* to $(\varrho, v)$ in $L^\infty((0, \infty) \times \mathbb{R}^d) \times L^\infty((0, \infty), H^1(\mathbb{R}^d))$ as $\nu$ goes to infinity, and $(\varrho, v)$ solves (uniquely) the inhomogeneous, incompressible model (1.25) with initial data $(\rho_0, u_0)$ in the distribution sense. Moreover, we have

$$ (1.26) \quad \left\{ \begin{array}{ll} \operatorname{div} u^{(\nu)} = \mathcal{O}(\nu^{-1/2}) & \text{in } L^2 \cap L^\infty((0, \infty), L^2(\mathbb{R}^d)), \\ \partial_t(\rho^{(\nu)}u^{(\nu)}) + \operatorname{div}(\rho^{(\nu)}u^{(\nu)} \otimes u^{(\nu)}) - \nabla F^{(\nu)} - \mu\Delta u^{(\nu)} = \mathcal{O}(\nu^{-1/2}) & \text{in } L^\infty((0, \infty), \dot{H}^{-1}(\mathbb{R}^d)), \end{array} \right. $$
where $F^{(\nu)} = \nu \operatorname{div} u^{(\nu)} - G^{(\nu)}$ with $G^{(\nu)} = P(\rho^{(\nu)}) - \tilde{P}$.

The proofs of Theorem 1.6 and Corollary 1.10 are presented in Section 2.3, based on the a priori estimates of Section 2.1 and their proofs in Section 2.2.

**Remark 1.11.** This result in Corollary 1.10 is a partial continuation of the work by Danchin and Mucha [13, 14, 16], and Danchin and Wang [17], and stands, as far as we know, as the first one dealing with discontinuous initial data in the whole space. We notice that, except for the work [13] dealing with the whole space case and initial data in the critical Besov space, the other studies rely heavily on the assumption that the domain has finite measure. The extension to the whole space, especially for $d=2$, is not obvious, and it requires some refined computations, e.g. the compensated result by Coifman, Lions, Meyer and Semmes in [8].

**1.3. Review of known results.** Classical solutions for the Navier-Stokes equations (1.1) with regular initial data are known to exist, since the works by Nash [42], Itaya [33, 34], Solonnikov [47], Tani [48], just to cite a few examples. These solutions are defined up to a positive time which depends on the (norms of) initial data. The first result addressing the global-in-time well-posedness of classical solutions is provided by Matsumura and Nishida [41] for small initial data in $L^1(\mathbb{R}^3) \cap H^3(\mathbb{R}^3)$. Nowadays, global-in-time classical solutions are known to exist under smallness assumption on the initial data in critical Besov space [4, 7, 26].

**Weak solutions and estimates for $E(t)$ and $G$.** Similar to the solutions constructed by Leray [35] for the incompressible Navier-Stokes equations, there are well-established results that investigate the existence of global-in-time weak solutions for the compressible Navier-Stokes equations (1.1), with finite initial energy. The first result was obtained by P.-L. Lions [40], followed by Feireisl, Novotný, Petzeltová [23], for pressure laws of the form $P(\rho) = a\rho^\gamma$, $a > 0$, with some limitations on $\gamma$. These weak solutions verify the following classical energy inequality:

$$ (1.27) \qquad E(t) \le E(0) = E_0, $$

where the functional $E$ has been given in (1.16).

The introduction of the (generalized) specific energy $H_l(\rho)$, $l \in [1, \infty)$, in (1.12) helps (technically) to estimate the pressure deviation $G$. As observed in e.g. [3], the so-defined $H_l(\rho)$ is non-negative: $H_l(\rho) \ge 0$, since the pressure $P(\rho)$ is an non-decreasing function of the density (1.10).

For the classical case $l=1$, $H_1(\rho)$ appears in the definition of $E(t)$, which is integrable in space uniformly in time due to (1.27): $\int_{\mathbb{R}^d} (H_1(\rho))(t, x) dx \le E(t) = E_0$. Consequently, under the a priori assumption

$$ \rho(t, x) \le \rho^*, $$
Consequently, integrating the above in space yields, after Hölder's inequality, the following

$$ (1.28) \qquad \frac{d}{dt} \|H_l(\rho)(t, x)\|_{L^1(\mathbb{R}^d)} + \frac{1}{\nu} \|G\|_{L^{l+1}(\mathbb{R}^d)}^{l+1} \le \frac{1}{\nu} \|G\|_{L^{l+1}(\mathbb{R}^d)}^l \|F\|_{L^{l+1}(\mathbb{R}^d)}^{l+1}, $$

and hence by Young's inequality and integration in time, $G$ can be controlled by $F$ in the following way

$$ (1.29) \qquad \frac{1}{\nu} \|G\|_{L^{l+1}((0,t)\times\mathbb{R}^d)}^{l+1} \le 2 \|H_l(\rho_0)(x)\|_{L^1(\mathbb{R}^d)} + \frac{C}{\nu} \|F\|_{L^{l+1}((0,t)\times\mathbb{R}^d)}^{l+1}, \quad \forall t \in (0, \infty). $$

**Density patch problem.** In the last three decades, there has been growing interests in exploring the properties of weak solutions to models arising from fluid mechanics that enable tracking down discontinuities of some quantities such as density or vorticity. We refer to the density patch problem for incompressible models stated in [39]: Consider the incompressible model (1.24) in two dimension with initial density as the characteristic function of some regular domain $D_0 \in \mathbb{R}^2$: $\rho_0 = \mathbb{1}_{D_0}$. The density-patch problem asks whether or not, at positive times, the density is still some characteristic function $\mathbb{1}_{D(t)}$ with the domain $D(t) \subset \mathbb{R}^2$ preserving the initial regularity of $D_0$. This problem is almost solved for incompressible models, even for density-dependent viscosity (under some smallness assumption) or higher Sobolev regularity of $D_0$, see [11, 12, 15, 18, 19, 20, 25, 36, 37, 38].

However, for a similar problem in the context of compressible fluids, there are not so many results. On one hand, the global classical solutions constructed by Matsumura and Nishida, or in critical Besov space, are too strong in a way that they do not allow for discontinuous solutions. On the other hand, the weak solutions constructed by P.-L. Lions or Feireisl, Novotný, Petzeltová only require that the initial energy is finite, allowing for discontinuous density. However, the velocity is relatively weak, with $\nabla u \in L^2((0, \infty) \times \mathbb{R}^d)$, and this is insufficient to track down discontinuities in the density. A natural idea is to construct weak solutions in a class that allows for tracking down the discontinuity of the interface. The first result addressing this issue is, as far as we know, [29] by Hoff, where the author considered an initial density with Hölder regularity on both sides of a suitable curve, allowing for jumps across this curve. The initial curve is transported by the flow of the velocity into a curve that maintains its initial regularity. The density also remains Hölder continuous on both sides of the transported curve, and moreover, its jump through the latter decays exponentially over time. This result pertains only in the case of linear pressure law and small bulk viscosity. Recently, these restrictions were removed in [49], even in the more challenging case of density-dependent viscosity. Theorem 1.6 is thus added to this list, in the constant viscosity setting, with domains having Sobolev regularity, and the density can be large in $L^\infty(\mathbb{R}^d)$, unlike the cited results.

**1.3.1. Hoff's strategy.** We review briefly some key concepts in Hoff's works [27, 28, 29, 30].

**Energy functionals.** In [27], Hoff introduced the following energy functionals which can be compared with our definitions in (1.15)

$$ \left\{ \begin{aligned} A_1^H(t) &= \sup_{[0,t]} \sigma \|\nabla u\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \sigma(t') \|\sqrt{\rho} \dot{u}(t')\|_{L^2(\mathbb{R}^d)}^2 dt', \\ A_2^H(t) &= \sup_{[0,t]} \sigma^d \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \sigma^d(t') \|\nabla \dot{u}(t')\|_{L^2(\mathbb{R}^d)}^2 dt', \\ \mathcal{B}(t) &= \sup_{[0,t]} \|\rho - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)}^2, \end{aligned} \right. $$

where the time weight $\sigma$ and the material derivative $\dot{u}$ are defined as in (1.7). He provides bounds for these functionals by requiring that the initial velocity is small in $L^2(\mathbb{R}^d)$ but can be large in $L^{2^d}(\mathbb{R}^d)$. Additionally, he requires that the initial density is bounded away from zero and bounded from above, along with some technical assumptions.

**Effective flux and vorticity.** Hoff's computations, mainly while propagating the lower and upper bounds of the density, rely strongly on the *effective viscous flux* $F$ given in (1.13): $F = \nu \operatorname{div} u - G$. It plays a crucial role by connecting the momentum equations and the mass equation, which was discovered by Hoff and Smoller in [31]. It presents its power also in the study of the propagation of oscillations in [45], and in the constructions of weak solutions in [23, 28, 40].
we bound $G(t, x) := P(\rho)(t, x) - \tilde{P}$ uniformly in time by the energy $E_0$ as below

$$ (1.28) \quad \sup_{[0,t]} \|G\|_{L^q(\mathbb{R}^d)}^q \le C^* \sup_{t' \in [0,t]} \int_{\mathbb{R}^d} H_1(\rho(t', x)) dx \le C^* E_0, \quad \text{with } q \in [2, \infty), $$

where the constant $C^*$ depends only on $\rho^*$ and $q$. In the above, the first inequality follows from the definition of $H_1(\rho)$ in (1.12).

General $H_1(\rho)$, $l \ge 1$, as a function of $\rho$, satisfies the following ordinary differential equation

$$ \rho H_l'(\rho) - H_l(\rho) = |P(\rho) - \tilde{P}|^{l-1}(P(\rho) - \tilde{P}), $$

and hence, by virtue of the mass equation (1.1)$_1$, the function $(H_l(\rho))(t, x)$ satisfies the following time evolutionary equation

$$ \partial_t H_l(\rho) + \operatorname{div}(H_l(\rho)u) + |P(\rho) - \tilde{P}|^{l-1}(P(\rho) - \tilde{P}) \operatorname{div} u = 0, $$

which is in the same spirit of the renormalized continuity equation appearing in e.g. [40]. By view of the definitions (1.13) and (1.14), it is equivalent to

$$ \partial_t H_l(\rho) + \operatorname{div}(H_l(\rho)u) + \frac{1}{\nu}|G|^{l+1} = -\frac{1}{\nu}|G|^{l-1}GF. $$

Consequently, integrating the above in space yields, after Hölder's inequality, the following

$$ (1.29) \quad \frac{d}{dt} \|H_l(\rho)(t, x)\|_{L^1(\mathbb{R}^d)} + \frac{1}{\nu} \|G\|_{L^{l+1}(\mathbb{R}^d)}^{l+1} \le \frac{1}{\nu} \|G\|_{L^{l+1}(\mathbb{R}^d)}^l \|F\|_{L^{l+1}(\mathbb{R}^d)}, $$

and hence by Young's inequality and integration in time, $G$ can be controlled by $F$ in the following way

$$ (1.30) \quad \frac{1}{\nu} \|G\|_{L^{l+1}((0,t) \times \mathbb{R}^d)}^{l+1} \le 2 \|H_l(\rho_0)(x)\|_{L^1(\mathbb{R}^d)} + \frac{C}{\nu} \|F\|_{L^{l+1}((0,t) \times \mathbb{R}^d)}^{l+1}, \quad \forall t \in (0, \infty). $$

**Density patch problem.** In the last three decades, there has been growing interests in exploring the properties of weak solutions to models arising from fluid mechanics that enable tracking down discontinuities of some quantities such as density or vorticity. We refer to the density patch problem for incompressible models stated in [39]: Consider the incompressible model (1.25) in two dimension with initial density as the characteristic function of some regular domain $D_0 \in \mathbb{R}^2$: $\rho_0 = \mathbb{1}_{D_0}$. The density-patch problem asks whether or not, at positive times, the density is still some characteristic function $\mathbb{1}_{D_0(t)}$ with the domain $D(t) \subset \mathbb{R}^2$ preserving the initial regularity of $D_0$. This problem is almost solved for incompressible models, even for density-dependent viscosity (under some smallness assumption) or higher Sobolev regularity of $D_0$, see [11, 12, 15, 18, 19, 20, 25, 36, 37, 38].

However, for a similar problem in the context of compressible fluids, there are not so many results. On one hand, the global classical solutions constructed by Matsumura and Nishida, or in critical Besov space, are too strong in a way that they do not allow for discontinuous solutions. On the other hand, the weak solutions constructed by P.-L. Lions or Feireisl, Novotný, Petzeltová only require that the initial energy is finite, allowing for discontinuous density. However, the regularity of the velocity is relatively weak, with $\nabla u \in L^2((0, \infty) \times \mathbb{R}^d)$, and this is insufficient to track down discontinuities in the density. A natural idea is to construct weak solutions in a class that allows for tracking down the discontinuity of the interface. The first result addressing this issue is, as far as we know, [29] by Hoff, where the author considered an initial density with Hölder regularity on both sides of a suitable curve, allowing for jumps across this curve. The initial curve is transported by the flow of the velocity into a curve that maintains its initial regularity. The density also remains Hölder continuous on both sides of the transported curve, and moreover, its jump through the latter decays exponentially over time. This result pertains only to the case of linear pressure law and small bulk viscosity. Recently, these restrictions were removed in [10, 50], where the later reference treated even more challenging case of density-dependent viscosity. Theorem 1.6 is thus added to this list, in the constant viscosity setting, with domains having Sobolev regularity, and the density can be large in $L^\infty(\mathbb{R}^d)$, unlike the cited results.
In fact, recall the momentum equations (1.1)$_2$, which can be written by virtue of the mass conservation law (1.1)$_1$ as

$$ \rho \dot{u} - \mu \Delta u - (\mu + \lambda) \nabla \operatorname{div} u + \nabla (P(\rho) - \tilde{P}) = 0. $$

We apply the divergence operator to it to obtain the Poisson equation for $F$ as follows

$$ (1.30) \qquad \Delta F = \operatorname{div}(\rho \dot{u}). $$

Similarly we can apply the curl operator to the momentum equations, to derive the Poisson equation for the vorticity, curl $u$, as follows

$$ (1.31) \qquad \mu \Delta \operatorname{curl} u = \operatorname{curl}(\rho \dot{u}). $$

Consequently, the regularity of the material derivative of the velocity $\dot{u}$, as provided by functionals $A_1^H$ and $A_2^H$, allows the effective flux $F$ and the vorticity curl $u$ to be regular at positive time, even for rough density. This means that there is some cancellation between the divergence of the velocity and the pressure at positive times. In particular, the fact that $F \in L^{8/3}((1, \infty), L^\infty(\mathbb{R}^d))$ allows him to propagate the $L^\infty(\mathbb{R}^d)$ estimate for the density.

Thanks to this observation, under smallness condition on the initial data, Hoff proved existence of global weak solutions for the system (1.1) with a linear pressure law in a first paper [27]. He later considered pressure laws of the form $P(\rho) = a\rho^\gamma$, with $\gamma > 1$ in a second paper [28], in which, again, the effective flux played a crucial role in proving compactness for the density.

**Velocity gradient expression involving Riesz operators.** In order to study the dynamics of discontinuous surfaces, Hoff in [29] used the following decomposition of the velocity gradient:

$$ (1.32) \qquad \mu \nabla u = -(-\Delta)^{-1} \nabla(\rho \dot{u}) + \frac{\mu + \lambda}{\nu} RRF + \frac{\mu}{\nu} RRG =: \mu \nabla \tilde{u} + \mu \nabla u_G, $$

which is nothing but a rewriting of the above momentum equations, where $R_j = (\frac{1}{t}\partial_j)(-\Delta)^{-\frac{1}{2}}, 1 \le j \le d$, are the Riesz operators.

By assuming more regularity on the velocity $u_0 \in H^\beta(\mathbb{R}^2)$, he reduces the singularity of time weighs in the definitions of functionals $A_1^H$ and $A_2^H$. Namely, in dimension two, the time weights $\sigma$ and $\sigma^2$ are replaced, respectively, by $\sigma^{1-\beta}$ and $\sigma^{2-\beta}$. Thus, $\nabla \tilde{u}$ and the effective flux $F$ belong to $L_{\text{loc}}^1([0, \infty), C^\alpha(\mathbb{R}^2))$ for all $0 < \alpha < \beta$. With the help of the regularity of $F$, Hoff propagated the piecewise Hölder regularity of the density, resulting piecewise Hölder continuity of $\nabla u_G$ on both sides of a time-dependent curve. This time-dependent curve is the transport of an initial suitable curve with some geometric assumptions, and only provided with bounded velocity gradient can the structure of the density and of the curve be propagated.

However, since Riesz operators fail to be continuous on $L^\infty(\mathbb{R}^d)$, additional regularity must be assumed on the density to obtain $RRG \in L^\infty(\mathbb{R}^d)$. In [30], Hoff and Santos observed that in the configuration of the previous works (see [27, 28]), the rough part of the velocity gradient $\nabla u_G$ belongs to $L^\infty((0, \infty), BMO(\mathbb{R}^d))$. In this case, the initial interface $\gamma_0 \in C^\alpha$, $\alpha > 0$, is transported to an interface $\gamma_t \in C^{\alpha_t}$ at time $t > 0$, with $\alpha_t$ decaying exponentially to zero in time.

Hence to propagate interface regularity (more than continuity) requires a Lipschitz velocity. For the incompressible model with constant viscosity, this regularity is directly obtained from energy computations and interpolations. In contrast, for the compressible case with discontinuous density, the problem is more delicate, and the issue is to find an appropriate functional space framework such that even-order Riesz operators maps into $L^\infty(\mathbb{R}^d)$.

### 1.3.2. The strategy by use of tangential regularity.
Apart from the tools used in [29, 49] to handle the rough part of the velocity gradient, there exists another framework that allows for the same. It is referred as tangential/striated regularity space, which goes back to Chemin's study (see e.g. [5, 6]) of the vortex patch problem for the ideal incompressible model. See also [2] for another interesting geometric proof for the persistence of regularity in the vortex patch problem. Chemin's idea has been further developed to higher dimensional cases in [9, 24], to the inhomogeneous case in [22], as well as the density patch problem of the inhomogeneous incompressible Navier-Stokes model in e.g. [36, 37, 38, 43]. However,
1.3.1. *Hoff's strategy.* We review briefly some key concepts in Hoff's works [27, 28, 29, 30].

**Energy functionals.** In [27], Hoff introduced the following energy functionals which can be compared with our definitions in (1.16)

$$
\left\{
\begin{aligned}
A_1^H(t) &= \sup_{[0,t]} \sigma \|\nabla u\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \sigma(t') \|\sqrt{\rho} \dot{u}(t')\|_{L^2(\mathbb{R}^d)}^2 dt', \\
A_2^H(t) &= \sup_{[0,t]} \sigma^d \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^d)}^2 + \int_0^t \sigma^d(t') \|\nabla \dot{u}(t')\|_{L^2(\mathbb{R}^d)}^2 dt', \\
\mathcal{B}(t) &= \sup_{[0,t]} \|\rho - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)}^2,
\end{aligned}
\right.
$$

where the time weight $\sigma$ and the material derivative $\dot{u}$ are defined as in (1.7). He provides bounds for these functionals by requiring that the initial velocity is small in $L^2(\mathbb{R}^d)$ but can be large in $L^{2^d}(\mathbb{R}^d)$. Additionally, he requires that the initial density is bounded away from zero and bounded from above, along with some technical assumptions.

**Effective flux and vorticity.** Hoff's computations, mainly while propagating the lower and upper bounds of the density, rely strongly on the *effective viscous flux* $F$ given in (1.14): $F = \nu \operatorname{div} u - G$. It plays a crucial role by connecting the momentum equations and the mass equation, which was discovered by Hoff and Smoller in [31]. It presents its power also in the study of the propagation of oscillations in [45], and in the constructions of weak solutions in [23, 28, 40].

In fact, recall the momentum equations (1.1)$_2$, which can be written by virtue of the mass conservation law (1.1)$_1$ as

$$ \rho \dot{u} - \mu \Delta u - (\mu + \lambda) \nabla \operatorname{div} u + \nabla (P(\rho) - \tilde{P}) = 0. $$

We apply the divergence operator to it to obtain the Poisson equation for $F$ as follows

$$ (1.31) \qquad \Delta F = \operatorname{div}(\rho \dot{u}). $$

Similarly we can apply the curl operator to the momentum equations, to derive the Poisson equation for the vorticity, curl $u$, as follows

$$ (1.32) \qquad \mu \Delta \operatorname{curl} u = \operatorname{curl}(\rho \dot{u}). $$

Consequently, the regularity of the material derivative of the velocity $\dot{u}$, as provided by functionals $A_1^H$ and $A_2^H$, allows the effective flux $F$ and the vorticity curl $u$ to be regular at positive time, even for rough density. This means that there is some cancellation between the divergence of the velocity and the pressure at positive times. In particular, the fact that $F \in L^{8/3}((1, \infty), L^\infty(\mathbb{R}^d))$ allows Hoff to propagate the $L^\infty(\mathbb{R}^d)$ estimate for the density.

Thanks to this observation, under a smallness condition on the initial data, Hoff proved existence of global weak solutions for the system (1.1) with a linear pressure law in a first paper [27]. He later considered pressure laws of the form $P(\rho) = a\rho^\gamma$, with $\gamma > 1$ in a second paper [28], in which, again, the effective flux played a crucial role in proving compactness for the density.

**Velocity gradient expression involving Riesz operators.** In order to study the dynamics of discontinuous surfaces, Hoff in [29] used the following decomposition of the velocity gradient:

$$ (1.33) \qquad \mu \nabla u = -(-\Delta)^{-1} \nabla (\rho \dot{u}) + \frac{\mu + \lambda}{\nu} \mathcal{R} \mathcal{R} F + \frac{\mu}{\nu} \mathcal{R} \mathcal{R} G := \mu \nabla \tilde{u} + \mu \nabla u_G, $$

which is nothing but a rewriting of the above momentum equations, where $\mathcal{R}_j = (\frac{1}{i}\partial_j)(-\Delta)^{-\frac{1}{2}}, 1 \le j \le d$, are the Riesz operators.

By assuming more regularity on the velocity $u_0 \in H^\beta(\mathbb{R}^2)$, he reduces the singularity of time weights in the definitions of functionals $A_1^H$ and $A_2^H$. Namely, in dimension two, the time weights $\sigma$ and $\sigma^2$ are replaced, respectively, by $\sigma^{1-\beta}$ and $\sigma^{2-\beta}$. Thus, $\nabla \tilde{u}$ and the effective flux $F$ belong to $L_{\text{loc}}^1([0, \infty), \mathcal{C}^\alpha(\mathbb{R}^2))$ for all $0 < \alpha < \beta$. With the help of the regularity of $F$, Hoff propagated the piecewise Hölder regularity of the density, resulting piecewise Hölder continuity of $\nabla u_G$ on both sides of a time-dependent curve. This time-dependent curve is the transport of an initial suitable curve with some geometric assumptions,
there are very few results in this direction for the compressible case. To the best of our knowledge, the only one work in the literature is [10] by Danchin, Fanelli, and Paicu. They establish the local-in-time well-posedness of the compressible equations (1.1) with a striated initial density, and we now delve into a brief discussion of their methods. From the momentum equations (1.1)$_2$, they express the velocity as:

$$ (1.33) \qquad u = w - \nabla(I_d - \Delta)^{-1}G, $$

where $w$ solves a parabolic equation with source term in some suitable space $L^r((0,T), L^p(\mathbb{R}^d))$. They employ maximal regularity tools to establish Lipschitz bounds for $w$. Meanwhile, Lipschitz bound for the second term of the velocity's expression (1.33), associated with the pressure, is obtained through tangential regularity estimates. The maximal regularity argument requires smallness assumption on the density in $L^\infty(\mathbb{R}^d)$, and the global-in-time result is still missing. Toward this, we establish local-in-time well-posedness of the system (1.1) without imposing any smallness condition on the initial data (see Theorem 1.4), and global-in-time well-posedness (see Theorem 1.6) without any smallness assumption of the initial density fluctuation in $L^\infty(\mathbb{R}^d)$, and the vacuum is allowed. This is accomplished through a coupling mechanism that involves the effective flux. By achieving this objective, we propagate the Sobolev regularity of interfaces over time.

**Incompressible limit.** We aim also to establish an incompressible limit in the spirit of the work of Danchin and Mucha [16]. Let us look briefly at this question. The work by Matsumura and Nishida [41] paved the way for attempts to relax the assumptions on the initial data. Despite reducing the regularity assumption to critical Besov space or even Lebesgue space, the condition of smallness is frequently encountered in the literature. In their work [13], Danchin and Mucha introduced a new framework that enables them to bypass the smallness condition on the initial data, namely, replacing the smallness in the initial data by large enough bulk viscosity coefficient. In particular, as the bulk viscosity $\nu \to \infty$, the solution converges to a limit that satisfies the incompressible model. This has been down for initial data in critical Besov space. For less regular initial data, they work on the torus in [14, 16], where they rely technically on the finite-measure of the domain, particularly on the logarithmic interpolation inequality, which proves to be crucial in handling vacuum states in [16]. We also refer to the work by Danchin and Wang [17] where exponential decay rate of the solutions of the compressible model on torus has been investigated. However, the exponential decay does not generally hold in the whole space. For instance, the work by Hu and Wu [32] provides lower bound for the norms of solutions in certain cases. We obtain similar results to those in [16] in the presence of vacuum on the whole space (see Corollary 1.10), where we apply technically, in dimension two. The extension to the whole space, particularly in two dimension, the compensated result by Coifman, Lions, Meyer, Semmes [8].

**Outline of the paper.** The rest of the paper is structured as follows. In next section Section 2, we give the proofs of Theorem 1.6 and Corollary 1.10, provided with the validity of Theorem 1.4, whose proof is postponed in Appendix B. A useful density-weighted interpolation inequality is established in Appendix A.

## 2. PROOF OF THE MAIN RESULTS

This section is devoted to the proofs of Theorem 1.6 and Corollary 1.10, which goes from a priori estimates for solutions of the Navier-Stokes equations (1.1) to the proof of the compactness of approximate solutions. It is divided into three parts. In the first one, Section 2.1, we summarise all key ideas with brief explanations and give the a priori estimates in a series of lemmas. Technical details and the proofs of these lemmas are presented in the second part Section 2.2. As we will see in the final part of the proof in Section 2.3, the existence of a local-in-time solution (without any smallness condition in the density) is by no means obvious, and it is the purpose of Appendix B. The regularity of the (local-in-time) solution is sufficient in order to use $u$ and $\dot{u}$ as test functions in the subsequent computations to get energy estimates.
and only provided with bounded velocity gradient can the structure of the density and of the curve be propagated.

However, since Riesz operators fail to be continuous on $L^\infty(\mathbb{R}^d)$, additional regularity must be assumed on the density to obtain $\mathcal{R}\mathcal{R}G \in L^\infty(\mathbb{R}^d)$. In [30], Hoff and Santos observed that in the configuration of the previous works (see [27, 28]), the rough part of the velocity gradient $\nabla u_G$ belongs to $L^\infty((0, \infty), \mathcal{B}MO(\mathbb{R}^d))$. In this case, the initial interface $\gamma_0 \in \mathcal{C}^\alpha$, $\alpha > 0$, is transported to an interface $\gamma_t \in \mathcal{C}^{\alpha_t}$ at time $t > 0$, with $\alpha_t$ decaying exponentially to zero in time.

Hence to propagate interface regularity (more than continuity) requires a Lipschitz velocity. For the incompressible model with constant viscosity, this regularity is directly obtained from energy computations and interpolations. In contrast, for the compressible case with discontinuous density, the problem is more delicate, and the issue is to find an appropriate functional space, which can be mapped by even-order Riesz operators into $L^\infty(\mathbb{R}^d)$.

### 1.3.2. The strategy by use of tangential regularity.
Apart from the tools used in [29, 50] to handle the rough part of the velocity gradient, there exists another framework that allows for the same. It is referred as tangential/striated regularity space, which goes back to Chemin's study (see e.g. [5, 6]) of the vortex patch problem for the ideal incompressible model. See also [2] for another interesting geometric proof for the persistence of regularity in the vortex patch problem. Chemin's idea has been further developed to higher dimensional cases in [9, 24], to the inhomogeneous case in [22], as well as to the density-patch problem for the inhomogeneous incompressible Navier-Stokes model in e.g. [36, 37, 38, 43]. However, there are very few results in this direction for the compressible case. To the best of our knowledge, the only work in the literature is [10] by Danchin, Fanelli, and Paicu. They establish the local-in-time well-posedness of the compressible equations (1.1) with a striated initial density, and we now delve into a brief discussion of their methods. From the momentum equations (1.1)$_2$, they express the velocity as:

$$ (1.34) \qquad u = w - \nabla(I_d - \Delta)^{-1}G, $$

where $w$ solves a parabolic equation with source term in some suitable space $L^r((0, T), L^p(\mathbb{R}^d))$. They employ maximal regularity tools to establish Lipschitz bounds for $w$. Meanwhile, Lipschitz bound for the second term of the velocity's expression (1.34), associated with the pressure, is obtained through tangential regularity estimates. The maximal regularity argument requires smallness assumption on the density in $L^\infty(\mathbb{R}^d)$, and the global-in-time result is still missing. Toward this, we establish local-in-time well-posedness of the system (1.1) without imposing any smallness condition on the initial data (see Theorem 1.4), and global-in-time well-posedness (see Theorem 1.6) without any smallness assumption of the initial density fluctuation in $L^\infty(\mathbb{R}^d)$, and the vacuum is allowed. This is accomplished through a coupling mechanism that involves the effective flux. By achieving this objective, we propagate the Sobolev regularity of interfaces over time.

### Incompressible limit.
We aim also to establish an incompressible limit in the spirit of the work of Danchin and Mucha [16]. Let us briefly look at this question. The work by Matsumura and Nishida [41] paved the way for attempts to relax the assumptions on the initial data. Despite reducing the regularity assumption to critical Besov space or even Lebesgue space, the condition of smallness is frequently encountered in the literature. In their work [13], Danchin and Mucha introduced a new framework that enables them to bypass the smallness condition on the initial data, namely, replacing the smallness in the initial data by large enough bulk viscosity coefficient. In particular, as the bulk viscosity $\nu \to \infty$, the solution converges to a limit that satisfies the incompressible model. This has been done for initial data in critical Besov space. For less regular initial data, they work on the torus in [14, 16], where they rely technically on the finite-measure of the domain, particularly on the logarithmic interpolation inequality, which proves to be crucial in handling vacuum states in [16]. We also refer to the work by Danchin and Wang [17] where exponential decay rate of the solutions of the compressible model on torus has been investigated. However, the exponential decay does not generally hold in the whole space. For instance, the work by Hu and Wu [32] provides lower bound for the norms of solutions in certain cases. We obtain similar results as those in [16] in the presence of vacuum on the whole space (see Corollary 1.10). Specifically in the two-dimensional case, we do some algebraic computations and succeed in applying
## 2.1. Proof ideas and statements of lemmas.
In this section, we give the main ideas of the proof of Theorem 1.6. We state the energy estimates for the solutions of the compressible Navier-Stokes equations (1.1) with initial data (1.2) satisfying (1.3). The tangential regularity (1.4) is assumed when we show the boundedness of the Lipschitz-norm of the velocity vector field as a second step. Recall the definitions of the energy functionals

$$E(t), \mathcal{A}_1(t), \mathcal{A}_2(t), \mathcal{A}_3(t)$$

together with the notations $\tilde{P}, H_l(\rho), G, F$ and $E_0, E_0^\nu, \rho_0^*$, given in Section 1.2.2.

In the literature (see e.g. [40] where finally only the energy inequality (1.26) was established for weak solutions) the following a priori energy equality for $E(t)$ was shown for strong solutions

$$(2.1) \qquad E(t) = E(0) =: E_0.$$

More precisely, it follows from taking the scalar product of the momentum equation (1.1)$_2$ with the velocity $u$ and then integrating in time and space. This energy balance (2.1) is going to be used freely in the proof, and we aim to show the estimates for $\mathcal{A}_1, \mathcal{A}_2, \mathcal{A}_3$.

In the following we state step by step

* In Section 2.1.1: Energy estimates for $\mathcal{A}_1, \mathcal{A}_2$ together with the boundedness of the density deviation $\|\rho - \tilde{\rho}\|_{L_{t,x}^\infty}$.

Under the assumption that the density is a priori upper bounded

$$(2.2) \qquad 0 \le \rho(t, x) \le \rho^*,$$

for some $\rho^* > 0$, we show first (local-in-time) a priori energy estimates for $\mathcal{A}_1, \mathcal{A}_2$ (see Lemma 2.1 and Lemma 2.2) and then a boundedness of the density in terms of $\mathcal{A}_1, \mathcal{A}_2$ (see Lemma 2.3) for solutions of the Cauchy problem (1.1)-(1.2)-(1.3). Under the assumption (1.20) or (1.21), that is, in the case of either large bulk viscosity coefficient for $d = 2$ or with small initial energy and large bulk viscosity coefficient for $d = 3$, a bootstrap argument implies the global-in-time a priori energy estimates for $\mathcal{A}_1, \mathcal{A}_2$, and density bound estimate (see Lemma 2.4).

* In Section 2.1.2: The striated regularity estimate for $\mathcal{A}_3$ together with the boundedness of the velocity gradient $\|\nabla u\|_{L_t^1 L_x^\infty}$.

With the estimates in Section 2.1.1 at hand, we turn to the striated regularity for the density function $\mathcal{A}_3(t)$ for solutions of the Cauchy problem (1.1)-(1.2)-(1.3)-(1.4), which finally implies the Lipschitz-continuity of the velocity field (see Lemma 2.6), thanks to the $L^\infty$-estimates for the double Riesz-operators provided with extra striated regularity (see Proposition 2.5).

### 2.1.1. A priori estimates for $\mathcal{A}_1, \mathcal{A}_2$ and $\|\rho - \tilde{\rho}\|_{L_{t,x}^\infty}$.
In order to derive higher-order energy estimates for the velocity $u$ and its material derivative

$$\dot{u} := (\partial_t + u \cdot \nabla)u,$$

we use first $\dot{u}$ as a test function in the weak formulation of the momentum equation (1.1)$_2$ to establish bounds for $\mathcal{A}_1$. The functional $\mathcal{A}_2$ emerges when, first rewriting the momentum equation (1.1)$_2$ with the effective flux $F$, and then applying the operator $\partial_t \cdot + \text{div}(\cdot u)$ to the resulting equation before testing it with $\dot{u}$.

In two dimension, the following estimates are valid for these functionals $\mathcal{A}_1, \mathcal{A}_2$.

**Lemma 2.1** (Preliminary energy estimates for $d = 2$). Assume that $d = 2$ and (2.2). Then the following a priori bounds hold true for the functionals $\mathcal{A}_1$ and $\mathcal{A}_2$:

$$(2.3) \qquad \mathcal{A}_1(t) \le C^* \left( E_0^\nu + \frac{1}{\nu^{3/2}} \mathcal{A}_1(t)(E_0 + \mathcal{A}_1(t)) \right) \exp(C^*E_0),$$

$$(2.4) \qquad \mathcal{A}_2(t) \le C^* \left( (E_0 + \frac{1}{\nu^4} E_0^2) + (1 + E_0 + \mathcal{A}_1(t)) \mathcal{A}_1(t) \right),$$

where the constant $C^*$ depends on $\mu, \nu$ and (increasingly) on $\rho^*$.
Coifman, Lions, Meyer, and Semmes' compensated integrability result in [8] to achieve uniformly in $\nu$ estimates for the energy functionals.

**Outline of the paper.** The rest of the paper is structured as follows. In next section Section 2, we give the proofs of Theorem 1.6 and Corollary 1.10, by use of the results in Theorem 1.4, whose proof is postponed in Appendix B. A useful density-weighted interpolation inequality is established in Appendix A.

## 2. PROOF OF THE MAIN RESULTS

This section is devoted to the proofs of Theorem 1.6 and Corollary 1.10, which go from a priori estimates for solutions of the Navier-Stokes equations (1.1) to the proof of the compactness of approximate solutions. It is divided into three parts. In the first one, Section 2.1, we summarise all key ideas with brief explanations and give the a priori estimates in a series of lemmas. Technical details and the proofs of these lemmas are presented in the second part Section 2.2. As we will see in the final part of the proof in Section 2.3, the existence of a local-in-time solution (without any smallness condition in the density) is by no means obvious, and it is the purpose of Appendix B. The regularity of the (local-in-time) solution is sufficient in order to use $u$ and $\dot{u}$ as test functions in the subsequent computations to get energy estimates.

### 2.1. Proof ideas and statements of lemmas.
In this section, we give the main ideas of the proof of Theorem 1.6. We state the energy estimates for the solutions of the compressible Navier-Stokes equations (1.1)-(1.10) with initial data (1.2) satisfying (1.3). The tangential regularity (1.4) is assumed when we show the boundedness of the Lipschitz-norm of the velocity vector field as a second step. Recall the definitions of the energy functionals

$$E(t), \mathcal{A}_1(t), \mathcal{A}_2(t), \mathcal{A}_3(t)$$

together with the notations $\tilde{P}, H_l(\rho), G, F$ and $E_0, E_0^\nu, \rho_0^*$, given in Section 1.2.2.

In the literature (see e.g. [40] where finally only the energy inequality (1.27) was established for weak solutions) the following a priori energy equality for $E(t)$ was shown for strong solutions

$$(2.1) \qquad E(t) = E(0) =: E_0.$$

More precisely, it follows from taking the scalar product of the momentum equation (1.1)$_2$ with the velocity $u$ and then integrating in time and space. This energy balance (2.1) is going to be used freely in the proof, and we aim to show the estimates for $\mathcal{A}_1, \mathcal{A}_2, \mathcal{A}_3$.

In the following we state step by step

* In Section 2.1.1: Energy estimates for $\mathcal{A}_1, \mathcal{A}_2$ together with the boundedness of the density deviation $\|\rho - \tilde{\rho}\|_{L_{t,x}^\infty}$.

Under the assumption that the density is a priori from above bounded:

$$(2.2) \qquad 0 \le \rho(t, x) \le \rho^*,$$

for some $\rho^* > 0$, we show first (local-in-time) a priori energy estimates for $\mathcal{A}_1, \mathcal{A}_2$ (see Lemma 2.1 and Lemma 2.2) and then a boundedness of the density in terms of $\mathcal{A}_1, \mathcal{A}_2$ (see Lemma 2.3) for solutions of the Cauchy problem (1.1)-(1.2)-(1.3). Under the assumption (1.21) or (1.22), that is, in the case of either large bulk viscosity coefficient for $d=2$ or with small initial energy and large bulk viscosity coefficient for $d=3$, a bootstrap argument implies the global-in-time a priori energy estimates for $\mathcal{A}_1, \mathcal{A}_2$, and density bound estimate (see Lemma 2.4).

* In Section 2.1.2: The striated regularity estimate for $\mathcal{A}_3$ together with the boundedness of the velocity gradient $\|\nabla u\|_{L_t^1 L_x^\infty}$.

With the estimates in Section 2.1.1 at hand, we turn to the striated regularity for the density function $\mathcal{A}_3(t)$ for solutions of the Cauchy problem (1.1)-(1.2)-(1.3)-(1.4), which finally implies the Lipschitz-continuity of the velocity field (see Lemma 2.6), thanks to the $L^\infty$-estimates for the double Riesz-operators provided with extra striated regularity (see Proposition 2.5).
The proof of **Lemma 2.1** is established through refined computations, and the compensated result by Coifman *et al.* [8] turns out to be crucial for achieving a uniform bound with respect to $\lambda$. We refer to **Section 2.2.2** below for the detailed proof.

For $d=3$, the following estimates hold true for functionals $A_1$ and $A_2$.

**Lemma 2.2** (Preliminary energy estimates for $d=3$). Assume that $d=3$ and (2.2). Then, the following estimates hold true for the functionals $A_1$ and $A_2$:

$$ (2.5) \qquad A_1(t) \le C^*(E_0^\nu + \frac{1}{\nu^2} E_0^{\frac{1}{3}}) + C E_0 A_1(t)^2, $$

$$ (2.6) \qquad A_2(t) \le C^* \left( \frac{1}{\nu} E_0^{\frac{1}{3}} + E_0 + E_0^2 + (1 + A_1(t)^2) A_1(t) \right). $$

Here $C$ depends on $\mu, \nu$, and $C^*$ depends on $\mu, \nu$ and (increasingly) on $\rho^*$.

The proof is given in **Section 2.2.3**. Let us point out that the computations in [14, 16, 17] depend heavily on the fact that the domain has finite measure. **Lemma 2.1** and **Lemma 2.2** are the first to provide bounds for the solution $(\rho, u)$ uniformly with respect to $\lambda$ (large) in the whole space, with only bounded density.

Based on the above estimates, it turns out that the functionals $A_1$ and $A_2$ are under control as long as the density is upper-bounded. Therefore, the next step is devoted to estimating the upper bound of the density, whose proof is given in **Section 2.2.4**.

**Lemma 2.3** (Density upper bound in terms of energies). Assume (2.2). Then the following bounds hold true for the density:

$$ (2.7) \qquad \|\rho - \tilde{\rho}\|_{L^{\infty}([0,t] \times \mathbb{R}^d)} \le \|\rho_0 - \tilde{\rho}\|_{L^{\infty}(\mathbb{R}^d)} \\ \qquad \qquad + \frac{C^*}{\nu^{\frac{1}{3}}} \times \begin{cases} (1 + E_0^{\frac{1}{18}})(E_0^{\frac{1}{6}} + \nu^{\frac{1}{6}} A_1(t)^{\frac{1}{6}})(A_1(t)^{\frac{1}{3}} + A_2(t)^{\frac{1}{3}}), & d=2, \\ (A_1(t)^{\frac{1}{2}} + A_2(t)^{\frac{1}{2}}), & d=3. \end{cases} $$

Finally, we notice that for $d=2$ there is a small factor $\frac{1}{\nu}$ (or its positive powers) before $A_1(t)$, $A_2(t)$ in the estimates (2.3) and (2.7), while $A_2(t)$ can be bounded by $A_1(t)$ and $\rho^*$ by (2.4). We can close the estimates in **Lemma 2.1**, **Lemma 2.2** and **Lemma 2.3** by a bootstrap argument as in e.g. [3, 16], which is not repeated here.

**Lemma 2.4** (Global-in-time estimates under the assumption (1.20) or (1.21)). There exist $c$ depending only on $\mu, \nu$ and $\nu_0 \ge \nu$ depending on $\mu, \nu, E_0, \|\nabla u_0\|_{L^2(\mathbb{R}^d)}$, $\rho_0^*$ such that

(1) If $d=2$ and $\nu \ge \nu_0$, then

$$ A_1(t) + A_2(t) \le C_0^\nu \quad \text{and} \quad \|\rho - \tilde{\rho}\|_{L^\infty([0,t] \times \mathbb{R}^2)} \le \|\rho_0 - \tilde{\rho}\|_{L^\infty(\mathbb{R}^2)} + (C_0^\nu)^{\frac{1}{2}}. $$

(2) If $d=3$, $E_0^\nu E_0 \le c$ and $\nu \ge \nu_0$, then

$$ A_1(t) + A_2(t) \le C_0^\nu \quad \text{and} \quad \|\rho - \tilde{\rho}\|_{L^\infty([0,t] \times \mathbb{R}^3)} \le \|\rho_0 - \tilde{\rho}\|_{L^\infty(\mathbb{R}^3)} + (C_0^\nu)^{\frac{1}{2}}. $$

Above, $C_0^\nu$ depends on $\mu, \nu, \rho_0^*$ and (superlinearly) on $E_0^\nu$.

2.1.2. A priori estimates for $A_3$ and $\|\nabla u\|_{L^1_t L^\infty_x}$. Now we have **Lemma 2.4**, which gives the (uniform) bounds of functionals $A_1, A_2$ and $\rho$. We use the notation $C_0$ below to denote some time-independent constant depending on the initial data as follows:

$$ (2.8) \qquad C_0 = C_0(\mu, \nu, m, d, p, \rho_0^*, E_0^\nu), $$

where $m, p$ appears in the initial condition (1.4). $C_0$ may vary from lines to lines and bounds in particular $A_1, A_2, \rho$. The next step is dedicated to translating these bounds into the tangential regularity estimates for the density, together with the Lipschitz norm of the velocity. As the tangential regularity $A_3$ can be transported by Lipschitz continuous flow, we sketch the idea to show Lipschitz continuity of $u$ as follows.
2.1.1. A priori estimates for $\mathcal{A}_1, \mathcal{A}_2$ and $\|\rho - \tilde{\rho}\|_{L_{t,x}^\infty}$. In order to derive higher-order energy estimates for the velocity $u$ and its material derivative

$$ \dot{u} := (\partial_t + u \cdot \nabla)u, $$

we use first $\dot{u}$ as a test function in the weak formulation of the momentum equation $(1.1)_2$ to establish bounds for $\mathcal{A}_1$. The functional $\mathcal{A}_2$ emerges when, first rewriting the momentum equation $(1.1)_2$ with the effective flux $F$, and then applying the operator $\partial_t \cdot + \text{div}(\cdot u)$ to the resulting equation before testing it with $\dot{u}$.

In dimension two, the following estimates are valid for these functionals $\mathcal{A}_1, \mathcal{A}_2$.

**Lemma 2.1** (Preliminary energy estimates for $d=2$). Assume that $d=2$ and (2.2). Then the following a priori bounds hold true for the functionals $\mathcal{A}_1$ and $\mathcal{A}_2$:

$$ (2.3) \qquad \mathcal{A}_1(t) \le C^* \left( E_0^\nu + \frac{1}{\nu^{3/2}} \mathcal{A}_1(t)(E_0 + \mathcal{A}_1(t)) \right) \exp(C^*E_0), $$

$$ (2.4) \qquad \mathcal{A}_2(t) \le C^* \left( (E_0 + \frac{1}{\nu^4} E_0^2) + (1 + E_0 + \mathcal{A}_1(t))\mathcal{A}_1(t) \right), $$

where the constant $C^*$ depends on $\mu, \nu$ and (increasingly) on $\rho^*$.

The proof of Lemma 2.1 is established through refined computations, and the compensated result by Coifman et al. [8] turns out to be crucial for achieving a uniform bound with respect to $\nu$. We refer to Section 2.2.2 below for the detailed proof.

For $d=3$, the following estimates hold true for functionals $\mathcal{A}_1$ and $\mathcal{A}_2$.

**Lemma 2.2** (Preliminary energy estimates for $d=3$). Assume that $d=3$ and (2.2). Then, the following estimates hold true for the functionals $\mathcal{A}_1$ and $\mathcal{A}_2$:

$$ (2.5) \qquad \mathcal{A}_1(t) \le C^* \left( E_0^\nu + \frac{1}{\nu^3} E_0^{\frac{1}{3}} \right) + C E_0 \mathcal{A}_1(t)^2, $$

$$ (2.6) \qquad \mathcal{A}_2(t) \le C^* \left( \frac{1}{\nu} E_0^{\frac{1}{3}} + E_0 + E_0^2 + (1 + \mathcal{A}_1(t)^2) \mathcal{A}_1(t) \right). $$

Here $C$ depends on $\mu, \nu$, and $C^*$ depends on $\mu, \nu$ and (increasingly) on $\rho^*$.

The proof is given in Section 2.2.3. Let us point out that the computations in [14, 16, 17] depend heavily on the fact that the domain has finite measure. Lemma 2.1 and Lemma 2.2 are the first to provide high regularity bounds for the solution $(\rho, u)$ uniformly with respect to $\lambda$ (large) in the whole space, with only bounded density.

Based on the above estimates, it turns out that the functionals $\mathcal{A}_1$ and $\mathcal{A}_2$ are under control (for large $\nu$) as long as the density is upper-bounded. Therefore, the next step is devoted to estimating the upper bound of the density, whose proof is given in Section 2.2.4.

**Lemma 2.3** (Density upper bound in terms of energies). Assume (2.2). Then the following bounds hold true for the density:

$$ (2.7) \qquad \|\rho - \tilde{\rho}\|_{L^\infty([0,t] \times \mathbb{R}^d)} \le \|\rho_0 - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)} \\ \qquad \qquad + \frac{C^*}{\nu^{\frac{1}{3}}} \times \begin{cases} (1 + E_0^{\frac{1}{18}})(E_0^{\frac{1}{6}} + \nu^{\frac{1}{6}} \mathcal{A}_1(t)^{\frac{1}{6}})(\mathcal{A}_1(t)^{\frac{1}{3}} + \mathcal{A}_2(t)^{\frac{1}{3}}), & d=2, \\ (\mathcal{A}_1(t)^{\frac{1}{2}} + \mathcal{A}_2(t)^{\frac{1}{2}}), & d=3. \end{cases} $$

Finally, we notice that for $d=2$ there is a small factor $\frac{1}{\nu}$ (or its positive powers) before $\mathcal{A}_1(t), \mathcal{A}_2(t)$ in the estimates (2.3) and (2.7), while $\mathcal{A}_2(t)$ can be bounded by $\mathcal{A}_1(t)$ and $\rho^*$ by (2.4). We can close the estimates in Lemma 2.1, Lemma 2.2 and Lemma 2.3 by a bootstrap argument as in e.g. [3, 16], which is not repeated here.

**Lemma 2.4** (Global-in-time estimates under the assumption (1.21) or (1.22)). There exist $c$ depending only on $\mu, \nu$ and $\nu_0 \ge \nu$ depending on $\mu, \nu, E_0, \|\nabla u_0\|_{L^2(\mathbb{R}^d)}$, $\rho_0^*$ such that
We first recall the following decomposition of the velocity gradient

$$ (2.9) \qquad \begin{aligned} \nabla u &= \nabla \tilde{u} + \nabla u_G \\ &:= \left( -\frac{1}{\nu} \mathcal{R}\mathcal{R}(-\Delta)^{-1} \operatorname{div}(\rho \dot{u}) - \frac{1}{\mu} \mathcal{R}\mathcal{R}(-\Delta)^{-1} \cdot \operatorname{curl}(\rho \dot{u}) \right) + \left( \frac{1}{\nu} \mathcal{R}\mathcal{R}G \right). \end{aligned} $$

where $\mathcal{R}_j = (\frac{1}{i}\partial_j)(-\Delta)^{-\frac{1}{2}}$, with $1 \le j \le d$, is the Riesz transform and $G = P(\rho) - \tilde{P}$. Indeed, we notice that the following expression

$$ (2.10) \qquad \Delta u^j = \partial_j \operatorname{div} u + \partial_k \operatorname{curl}_{jk} u, \quad \operatorname{curl}_{jk} u = \partial_k u^j - \partial_j u^k, \quad j, k = 1, \dots, d, $$

and from (1.13), (1.30) and (1.31) we have

$$ (2.11) \qquad \operatorname{div} u = \frac{1}{\nu} (F + G), \quad F = -(-\Delta)^{-1} \operatorname{div}(\rho \dot{u}) \quad \text{and} \quad \mu \operatorname{curl} u = -(-\Delta)^{-1} \operatorname{curl}(\rho \dot{u}). $$

Hence the velocity gradient can be expressed as in (2.9):

$$ \begin{aligned} \nabla u &= -\nabla(-\Delta)^{-1}\Delta u \\ &= -\nabla(-\Delta)^{-1}\nabla \operatorname{div} u - \nabla(-\Delta)^{-1} \operatorname{div}(\operatorname{curl} u) \\ &= -\frac{1}{\nu}\nabla(-\Delta)^{-1}\nabla(F+G) - \frac{1}{\mu}\nabla(-\Delta)^{-1} \operatorname{div}(\mu \operatorname{curl} u) \\ &= \left( -\frac{1}{\nu}\mathcal{R}\mathcal{R}(-\Delta)^{-1} \operatorname{div}(\rho \dot{u}) - \frac{1}{\mu}\mathcal{R}\mathcal{R}(-\Delta)^{-1} \cdot \operatorname{curl}(\rho \dot{u}) \right) + \left( \frac{1}{\nu}\mathcal{R}\mathcal{R}G \right). \end{aligned} $$

Thanks to the regularity of $\dot{u}$ provided by functionals $\mathcal{A}_1$ and $\mathcal{A}_2$, $\nabla \tilde{u} \in L^1((0, t), L^\infty(\mathbb{R}^d))$. Motivated by the pioneering work of Chemin [5, 6] and Danchin, Fanelli, and Paicu [10] which show $\mathcal{R}\mathcal{R}G \in L^\infty(\mathbb{R}^d)$ provided with extra tangential regularity on $G$, the $L^\infty$-bound for $\nabla u_G$ in our case relies on the following logarithmic inequality, which is simply the application of the Sobolev embedding $L^\infty(\mathbb{R}^d) \cap \dot{W}^{1,p}(\mathbb{R}^d) \subset \mathcal{C}^{1-\frac{d}{p}}(\mathbb{R}^d)$ to Theorem 7.40 of [1].

**Proposition 2.5** ([1], $L^\infty$-bound for double Riesz transforms provided with tangential regularity). Let $\mathcal{X} = (X_\nu)_{1 \le \nu \le m} \subset L^{\infty,p}(\mathbb{R}^d)$, with $d < p < \infty$, be a non-degenerate family of $m \in \mathbb{N}^*$ vector fields as in **Section 1.1**. Let $1 \le q < \infty$.

There exists a constant $C = C(m, d, p, q) > 0$ such that for all $G \in L^q(\mathbb{R}^d) \cap L^\infty(\mathbb{R}^d) \cap L^p_{\mathcal{X}}(\mathbb{R}^d)$, the following estimate holds true:

$$ (2.12) \qquad \|\mathcal{R}\mathcal{R}G\|_{L^\infty(\mathbb{R}^d)} \le C\|G\|_{L^q(\mathbb{R}^d)} + C\|G\|_{L^\infty(\mathbb{R}^d)} \left( 1 + \log \left( e + \frac{\|G\|_{L^p_{\mathcal{X}}(\mathbb{R}^d)}}{\|G\|_{L^\infty(\mathbb{R}^d)}} \right) \right). $$

With the aid of the above logarithmic estimate, we can propagate tangential regularity of density and achieve Lipschitz regularity of the velocity at the same time.

**Lemma 2.6** (Tangential regularity for the density and Lipschitz continuity for the velocity). Assume the initial condition (1.4) that $\rho_0 \in L_{\mathcal{X}_0}^p(\mathbb{R}^d)$ where $\mathcal{X}_0 = (X_{0,v})_{1 \le v \le m} \subset L^{\infty,p}(\mathbb{R}^d)$ is a non-degenerate family of $m \in \mathbb{N}^*$ vectors fields, with $m \ge d-1$, with $2 < p < \infty$ if $d=2$ or $3 < p < 6$ if $d=3$.

Then, the family of vector fields $\mathcal{X}(t) = (X_\nu(t))_{1 \le \nu \le m}$, defined as solution of the Cauchy problem (1.8), is non-degenerate and $\mathcal{X}(t) \subset L^{\infty,p}(\mathbb{R}^d)$. Moreover, $\rho(t) \in L_{\mathcal{X}(t)}^p(\mathbb{R}^d)$ and the following bounds hold true:

$$ (2.13) \qquad \left\{ \begin{array}{l} A_3(t) \le A_3(0) \exp \left( C_0 \int_0^t \left[ 1 + \sqrt{t} + \|\nabla u(t')\|_{L^\infty(\mathbb{R}^d)} dt' \right] \right), \\ \int_0^t \|\nabla u(t')\|_{L^\infty(\mathbb{R}^d)} dt' \le C_0 \left( 1 + \frac{A_3(0)}{I(\mathcal{X}_0)} \right) \exp(C_0 t). \end{array} \right. $$

The proof of the above Lemma is the object of **Section 2.2.5**.
(1) If $d=2$ and $\nu \ge \nu_0$, then

$$A_1(t) + A_2(t) \le C_0^\nu \quad \text{and} \quad \|\rho - \tilde{\rho}\|_{L^\infty([0,t] \times \mathbb{R}^2)} \le \|\rho_0 - \tilde{\rho}\|_{L^\infty(\mathbb{R}^2)} + (C_0^\nu)^{\frac{1}{2}}.$$

(2) If $d=3$, $E_0^\nu E_0 \le c$ and $\nu \ge \nu_0$, then

$$A_1(t) + A_2(t) \le C_0^\nu \quad \text{and} \quad \|\rho - \tilde{\rho}\|_{L^\infty([0,t] \times \mathbb{R}^3)} \le \|\rho_0 - \tilde{\rho}\|_{L^\infty(\mathbb{R}^3)} + (C_0^\nu)^{\frac{1}{2}}.$$

Above, $C_0^\nu$ depends on $\mu, \nu, \rho_0^*$ and (superlinearly) on $E_0^\nu$.

2.1.2. A priori estimates for $A_3$ and $\|\nabla u\|_{L^1_t L^\infty_x}$. Now we have Lemma 2.4, which gives the (uniform) bounds of functionals $A_1, A_2$ and $\rho$. We use the notation $C_0$ below to denote some time-independent constant depending on the initial data as follows:

$$(2.8) \qquad C_0 = C_0(\mu, \nu, m, d, p, \rho_0^*, E_0^\nu),$$

where $m, p$ appears in the initial condition (1.4). $C_0$ may vary from lines to lines and controls in particular $A_1, A_2, \rho$. The next step is dedicated to translating these bounds into the tangential regularity estimates for the density, together with the Lipschitz norm of the velocity. As the tangential regularity $A_3$ can be transported by Lipschitz continuous flow, we sketch the idea to show Lipschitz continuity of $u$ as follows.

We first recall the following decomposition of the velocity gradient

$$(2.9) \qquad \begin{aligned} \nabla u &= \nabla \tilde{u} + \nabla u_G \\ &:= \left( -\frac{1}{\nu} \mathcal{R}\mathcal{R}(-\Delta)^{-1} \operatorname{div}(\rho \dot{u}) - \frac{1}{\mu} \mathcal{R}\mathcal{R}(-\Delta)^{-1} \cdot \operatorname{curl}(\rho \dot{u}) \right) + \left( \frac{1}{\nu} \mathcal{R}\mathcal{R}G \right). \end{aligned}$$

where $\mathcal{R}_j = (\frac{1}{i}\partial_j)(-\Delta)^{-\frac{1}{2}}$, with $1 \le j \le d$, is the Riesz transform and $G = P(\rho) - \tilde{P}$. Indeed, we notice that the following expression

$$(2.10) \qquad \Delta u^j = \partial_j \operatorname{div} u + \partial_k \operatorname{curl}_{jk} u, \quad \text{with} \quad \operatorname{curl}_{jk} u = \partial_k u^j - \partial_j u^k, \quad j, k = 1, \dots, d,$$

and from (1.14), (1.31) and (1.32) we have

$$(2.11) \qquad \operatorname{div} u = \frac{1}{\nu} (F + G), \quad F = -(-\Delta)^{-1} \operatorname{div}(\rho \dot{u}) \quad \text{and} \quad \mu \operatorname{curl} u = -(-\Delta)^{-1} \operatorname{curl}(\rho \dot{u}).$$

Hence the velocity gradient can be expressed as in (2.9):

$$\begin{aligned} \nabla u &= -\nabla(-\Delta)^{-1}\Delta u \\ &= -\nabla(-\Delta)^{-1}\nabla \operatorname{div} u - \nabla(-\Delta)^{-1} \operatorname{div}(\operatorname{curl} u) \\ &= -\frac{1}{\nu}\nabla(-\Delta)^{-1}\nabla(F+G) - \frac{1}{\mu}\nabla(-\Delta)^{-1} \operatorname{div}(\mu \operatorname{curl} u) \\ &= \left( -\frac{1}{\nu}\mathcal{R}\mathcal{R}(-\Delta)^{-1} \operatorname{div}(\rho \dot{u}) - \frac{1}{\mu}\mathcal{R}\mathcal{R}(-\Delta)^{-1} \cdot \operatorname{curl}(\rho \dot{u}) \right) + \left( \frac{1}{\nu}\mathcal{R}\mathcal{R}G \right). \end{aligned}$$

Thanks to the regularity of $\dot{u}$ provided by functionals $A_1$ and $A_2$, $\nabla \tilde{u} \in L^1((0, t), L^\infty(\mathbb{R}^d))$. Motivated by the pioneering work of Chemin [5, 6] and Danchin, Fanelli, and Paicu [10] which show $\mathcal{R}\mathcal{R}G \in L^\infty(\mathbb{R}^d)$ provided with extra tangential regularity on $G$, the $L^\infty$-bound for $\nabla u_G$ in our case relies on the following logarithmic inequality, which is simply the application of the Sobolev embedding $L^\infty(\mathbb{R}^d) \cap \dot{W}^{1,p}(\mathbb{R}^d) \subset \mathcal{C}^{1-\frac{d}{p}}(\mathbb{R}^d)$ to Theorem 7.40 of [1].

**Proposition 2.5** ([1], $L^\infty$-bound for double Riesz transforms provided with tangential regularity). Let $\mathcal{X} = (X_\nu)_{1 \le \nu \le m} \subset L^\infty, p(\mathbb{R}^d)$, with $d < p < \infty$, be a non-degenerate family of $m \in \mathbb{N}^*$ vector fields as in Section 1.1. Let $1 \le q < \infty$.

There exists a constant $C = C(m, d, p, q) > 0$ such that for all $G \in L^q(\mathbb{R}^d) \cap L^\infty(\mathbb{R}^d) \cap L^p_\mathcal{X}(\mathbb{R}^d)$, the following estimate holds true:

$$(2.12) \qquad \|\mathcal{R}\mathcal{R}G\|_{L^\infty(\mathbb{R}^d)} \le C\|G\|_{L^q(\mathbb{R}^d)} + C\|G\|_{L^\infty(\mathbb{R}^d)} \left( 1 + \log \left( e + \frac{\|G\|_{L^p_\mathcal{X}(\mathbb{R}^d)}}{\|G\|_{L^\infty(\mathbb{R}^d)}} \right) \right).$$
**Remark 2.7** (Improved time regularity). We have the following improved time regularity, which is required for the uniqueness result, see e.g. [10, Equation (4.31)]: For some $t_0 > 0$,

$$ \int_0^{t_0} \sigma(t')^s \|\nabla u(t')\|_{L^\infty(\mathbb{R}^d)}^2 dt' < \infty, $$

where $s = 4/9$ if $d = 2$ and $s = 1/2$ if $d = 3$. Indeed, we apply Hölder's inequality with respect to the time variable to (2.66) in the proof in Section 2.2.5 to obtain (noticing (2.67) and (2.69))

$$ \int_0^t \|\nabla u_G(t')\|_{L^\infty(\mathbb{R}^d)}^2 dt' \le C_0 t \left(1 + \frac{A_3(0)}{I(\mathcal{X}_0)} + t + \int_0^t \|\nabla u(t')\|_{L^\infty(\mathbb{R}^d)}\right)^2 $$

and similarly as in the proof of (2.68) and (2.70), we have

$$ \int_0^t \sigma^{\frac{4}{9}} \|\nabla \tilde{u}\|_{L^\infty(\mathbb{R}^2)}^2 \le C_0(1 + t^{\frac{1}{3}}), \quad \int_0^t \sqrt{\sigma} \|\nabla \tilde{u}\|_{L^\infty(\mathbb{R}^3)}^2 \le C_0. $$

To complete the proof of **Theorem 1.6**, we need to construct an approximate sequence $(\rho^\delta, u^\delta)_\delta$ globally defined in time that converges to $(\rho, u)$, the unique solution of (1.1). Once this is done, we will have obtained a sequence $(\rho^{(\nu)}, u^{(\nu)})$ of solutions to (1.1), and the last step will be to justify that this sequence converges to some $(\varrho, v)$ that solves the inhomogeneous incompressible model. This is the purpose of **Section 2.3**, and, as we will see, the local solutions constructed in [10] cannot serve as building blocks. Thus, we will need to establish local well-posedness for the system (1.1) in **Appendix B**.

**2.2. Proofs.** In this subsection we give the detailed proofs of **Lemma 2.1**, **Lemma 2.2**, **Lemma 2.3** and **Lemma 2.6**. Before that we recall some basic facts, which will be used freely in the proofs below.

**2.2.1. Basic facts.** Under the assumption (2.2), we have the $L^\infty((0, t), L^q(\mathbb{R}^d))$ estimate for the vibration of the pressure term $G(t, x) = (P(\rho))(t, x) - \tilde{P}$ given in (1.27)

$$ (2.14) \qquad \|G\|_{L^\infty([0,t], L^q(\mathbb{R}^d))} \le C^*(E_0)^{\frac{1}{q}}, \quad \text{with } q \in [2, \infty], $$

where $C^*$ depends on $q, \rho^*$. Here the case $q = \infty$ follows straightforwardly from the definition. Recall also the estimate (1.29) for $G$ by $F$: For any $l \ge 1$,

$$ (2.15) \qquad \frac{1}{\nu} \|G\|_{L^{l+1}((0,t)\times\mathbb{R}^d)}^{l+1} \le C^*(l)E_0 + \frac{C}{\nu} \|F\|_{L^{l+1}((0,t)\times\mathbb{R}^d)}^{l+1}, \quad \forall t \in (0, \infty), $$

where we estimated $\|H_l(\rho_0)\|_{L^1(\mathbb{R}^d)}$ by $C^*(l)E_0$.

Recall also the relations (2.10):

$$ \Delta u^j = \partial_j \operatorname{div} u + \partial_k \operatorname{curl}_{jk} u, \quad j, k = 1, \dots, d, $$

and (2.11) between $\operatorname{div} u, F, G, \rho \dot{u}$, $\operatorname{curl} u$:

$$ (2.16) \qquad \operatorname{div} u = \frac{1}{\nu} (F + G), \quad F = -(-\Delta)^{-1} \operatorname{div}(\rho \dot{u}) \quad \text{and} \quad \mu \operatorname{curl} u = -(-\Delta)^{-1} \operatorname{curl}(\rho \dot{u}). $$

By use of the $L^q(\mathbb{R}^d)$, $q \in (1, \infty)$, $d \ge 2$-boundedness of Riesz operators, the following estimates follow immediately

$$ (2.17) \qquad \|\nabla u\|_{L^q(\mathbb{R}^d)} \le C(q, d) (\|\operatorname{div} u\|_{L^q(\mathbb{R}^d)} + \|\operatorname{curl} u\|_{L^q(\mathbb{R}^d)}), $$

$$ (2.18) \qquad \|\nabla F\|_{L^q(\mathbb{R}^d)} + \mu \|\nabla \operatorname{curl} u\|_{L^q(\mathbb{R}^d)} \le C(q, d) \|\rho \dot{u}\|_{L^q(\mathbb{R}^d)}. $$

We now recall the compensated result by Coifman, Lions, Meyer, Semmes [8] in dimension two.

**Proposition 2.8** (Coifman-Lions-Meyer-Semmes' Estimate for $d = 2$). Let $v, w \in \dot{H}^1(\mathbb{R}^2; \mathbb{R})$ be two functions and let us define

$$ g = \det \begin{pmatrix} \partial_1 v & \partial_2 v \\ \partial_1 w & \partial_2 w \end{pmatrix}. $$
With the aid of the above logarithmic estimate, we can propagate tangential regularity of density and achieve Lipschitz regularity of the velocity at the same time.

**Lemma 2.6** (Tangential regularity for the density and Lipschitz continuity for the velocity). Assume the initial condition (1.4) that $\rho_0 \in \mathbb{L}_{\mathcal{X}_0}^p(\mathbb{R}^d)$ where $\mathcal{X}_0 = (X_{0,v})_{1\le v\le m} \subset \mathbb{L}^{\infty,p}(\mathbb{R}^d)$ is a non-degenerate family of $m \in \mathbb{N}^*$ vectors fields, with $m \ge d-1$, with $2<p<\infty$ if $d=2$ or $3<p<6$ if $d=3$.

Then, the family of vector fields $\mathcal{X}(t) = (X_v(t))_{1\le v\le m}$, defined as solution of the Cauchy problem (1.8), is non-degenerate and $\mathcal{X}(t) \subset \mathbb{L}^{\infty,p}(\mathbb{R}^d)$. Moreover, $\rho(t) \in \mathbb{L}_{\mathcal{X}(t)}^p(\mathbb{R}^d)$ and the following bounds hold true:

$$ (2.13) \qquad \begin{cases} A_3(t) \le A_3(0) \exp(C_0 \int_0^t [1 + \sqrt{t} + \|\nabla u(t')\|_{L^\infty(\mathbb{R}^d)} dt']) \\ \int_0^t \|\nabla u(t')\|_{L^\infty(\mathbb{R}^d)} dt' \le C_0(1 + \frac{A_3(0)}{I(\mathcal{X}_0)}) \exp(C_0 t). \end{cases} $$

The proof of the above Lemma is the object of Section 2.2.5.

**Remark 2.7** (Improved time regularity). We have the following improved time regularity, which is required for the uniqueness result, see e.g. [10, Equation (4.31)]: For some $t_0 > 0$,

$$ \int_0^{t_0} \sigma(t')^s \|\nabla u(t')\|_{L^\infty(\mathbb{R}^d)}^2 dt' < \infty, $$

where $s = 4/9$ if $d=2$ and $s = 1/2$ if $d=3$. Indeed, we apply Hölder's inequality with respect to the time variable to (2.66) in the proof in Section 2.2.5 to obtain (noticing (2.67) and (2.69))

$$ \int_0^t \|\nabla u_G(t')\|_{L^\infty(\mathbb{R}^d)}^2 dt' \le C_0 t \left(1 + \frac{A_3(0)}{I(\mathcal{X}_0)} + t + \int_0^t \|\nabla u(t')\|_{L^\infty(\mathbb{R}^d)}^2 \right) $$

and similarly as in the proof of (2.68) and (2.70), we have

$$ \int_0^t \sigma^{\frac{4}{9}} \|\nabla \tilde{u}\|_{L^\infty(\mathbb{R}^2)}^2 \le C_0(1 + t^{\frac{1}{3}}), \quad \int_0^t \sqrt{\sigma} \|\nabla \tilde{u}\|_{L^\infty(\mathbb{R}^3)}^2 \le C_0. $$

To complete the proof of Theorem 1.6, we need to construct an approximate sequence $(\rho^\delta, u^\delta)_\delta$ globally defined in time that converges to $(\rho, u)$, the unique solution of (1.1). Once this is done, we will have obtained a sequence $(\rho^{(\nu)}, u^{(\nu)})$ of solutions to (1.1), and the last step will be to justify that this sequence converges to some $(\rho, v)$ that solves the inhomogeneous incompressible model. This is the purpose of Section 2.3, and, as we will see, the local solutions constructed in [10] cannot serve as building blocks. Thus, we will need to establish local well-posedness for the system (1.1) in Appendix B.

**2.2. Proofs.** In this subsection we give the detailed proofs of Lemma 2.1, Lemma 2.2, Lemma 2.3 and Lemma 2.6. Before that we recall some basic facts, which will be used freely in the proofs below.

**2.2.1. Basic facts.** Under the assumption (2.2), we have the $L^\infty((0,t), L^q(\mathbb{R}^d))$ estimate for the pressure fluctuation term $G(t,x) = (P(\rho))(t,x) - \tilde{P}$ given in (1.28)

$$ (2.14) \qquad \|G\|_{L^\infty([0,t], L^q(\mathbb{R}^d))} \le C^*(E_0)^{\frac{1}{q}}, \quad \text{with } q \in [2, \infty], $$

where $C^*$ depends on $q, \rho^*$. Here the case $q = \infty$ follows straightforwardly from the definition. Recall also the estimate (1.30) for $G$ by $F$: For any $l \ge 1$,

$$ (2.15) \qquad \frac{1}{\nu} \|G\|_{L^{l+1}((0,t) \times \mathbb{R}^d)}^{l+1} \le C^*(l)E_0 + \frac{C}{\nu} \|F\|_{L^{l+1}((0,t) \times \mathbb{R}^d)}^{l+1}, \quad \forall t \in (0, \infty), $$

where we estimated $\|H_l(\rho_0)\|_{L^1(\mathbb{R}^d)}$ by $C^*(l)E_0$.

Recall also the relations (2.10):

$$ \Delta u^j = \partial_j \operatorname{div} u + \partial_k \operatorname{curl}_{jk} u, \quad j, k = 1, \dots, d, $$

and (2.11) between $\operatorname{div} u, F, G, \rho\dot{u}$, $\operatorname{curl} u$:

$$ (2.16) \qquad \operatorname{div} u = \frac{1}{\nu} (F + G), \quad F = -(-\Delta)^{-1} \operatorname{div}(\rho \dot{u}) \quad \text{and} \quad \mu \operatorname{curl} u = -(-\Delta)^{-1} \operatorname{curl}(\rho \dot{u}). $$
Then $g$ belongs to the Hardy space $\mathcal{H}^1(\mathbb{R}^2)$, whence for all $f \in BMO(\mathbb{R}^2)$ we have the estimate

$$ \left| \int_{\mathbb{R}^2} f(x)g(x)dx \right| \le \|f\|_{BMO(\mathbb{R}^2)} \|\nabla v\|_{L^2(\mathbb{R}^2)} \|\nabla w\|_{L^2(\mathbb{R}^2)}. $$

In particular, since $\dot{H}^1(\mathbb{R}^2) \hookrightarrow BMO(\mathbb{R}^2)$, for all $f \in \dot{H}^1(\mathbb{R}^2)$ we have:

$$ \left| \int_{\mathbb{R}^2} f(x)g(x)dx \right| \le \|\nabla f\|_{L^2(\mathbb{R}^2)} \|\nabla v\|_{L^2(\mathbb{R}^2)} \|\nabla w\|_{L^2(\mathbb{R}^2)}. $$

It helps in the proof, by virtue of the following equalities in dimension two

$$ (2.19) \qquad \nabla u^j \cdot \nabla u^k \partial_k u^j = \operatorname{div} u\{|\nabla u|^2 - \det(\nabla u)\}, \quad \nabla u^l \cdot \partial_l u = (\operatorname{div} u)^2 - 2 \det(\nabla u). $$

Here and in the following we use Einstein's summation convention for repeated indices, unless otherwise claimed.

### 2.2.2. Proof of Lemma 2.1 for $d=2$. This paragraph is devoted to obtaining bounds for functionals $A_1$ and $A_2$ as defined in (1.15) for $d=2$, provided with bounded density function (2.2). The constants in the following estimates may depend on the viscosity coefficient $\mu$ and the lower bound $\underline{\nu}$ for $\nu$, while not on the viscosity coefficient $\nu$ which will be chosen to be big.

**Proof of (2.3).** The functional $A_1$ arises while using $\dot{u}$ as a test functional in the weak formulation of the momentum equation (1.1)$_2$. By doing so, one obtains the following equality:

$$ \begin{aligned} A_1(t) ={}& \frac{\mu}{2} \|\nabla u_0\|_{L^2(\mathbb{R}^d)}^2 + \frac{\mu + \lambda}{2} \|\operatorname{div} u_0\|_{L^2(\mathbb{R}^d)}^2 - \mu \int_0^t \int_{\mathbb{R}^d} \nabla u^j \cdot \nabla u^k \partial_k u^j + \frac{\mu}{2} \int_0^t \int_{\mathbb{R}^d} |\nabla u|^2 \operatorname{div} u \\ & + \frac{\mu + \lambda}{2} \int_0^t \int_{\mathbb{R}^d} (\operatorname{div} u)^3 - (\mu + \lambda) \int_0^t \int_{\mathbb{R}^d} \operatorname{div} u \nabla u^l \cdot \partial_l u + \int_0^t \int_{\mathbb{R}^d} \nabla u^l \cdot \partial_l u G \\ & + \int_{\mathbb{R}^d} \operatorname{div} u(s) G(s) \Big|_{s=0}^{s=t} + \int_0^t \int_{\mathbb{R}^d} (\rho P'(\rho) - P(\rho) + \tilde{P})(\operatorname{div} u)^2. \end{aligned} \quad (2.20) $$

**Step 1: Reformulation of the energy equality.** In the following lines, we will reformulate the terms appearing in the right hand side above by use of (2.16) and (2.19).

By (2.16) and (2.19), the sum of the third and the fourth terms on the right hand side of (2.20) can be reduced as follows:

$$ (2.21) \qquad \mu \int_0^t \int_{\mathbb{R}^2} \operatorname{div} u \left[ \det(\nabla u) - \frac{1}{2} |\nabla u|^2 \right] \\ ={}& \frac{\mu}{\nu} \int_0^t \int_{\mathbb{R}^2} F \det(\nabla u) - \frac{\mu}{2\nu} \int_0^t \int_{\mathbb{R}^2} F |\nabla u|^2 + \frac{\mu}{\nu} \int_0^t \int_{\mathbb{R}^2} G \left[ \det(\nabla u) - \frac{1}{2} |\nabla u|^2 \right], $$

and similarly, the sum of the second and the third terms of the second line of (2.20) read

$$ (2.22) \qquad - \frac{\mu + \lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} F \nabla u^l \cdot \partial_l u + \frac{\mu}{\nu} \int_0^t \int_{\mathbb{R}^2} \nabla u^l \cdot \partial_l u G \\ ={}& 2 \frac{\mu + \lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} F \det(\nabla u) - \frac{\mu + \lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} F (\operatorname{div} u)^2 + \frac{\mu}{\nu} \int_0^t \int_{\mathbb{R}^2} \nabla u^l \cdot \partial_l u G. $$

Now we pack the first term of the second line of (2.20) and the middle term in the above (2.22) and use (2.16) to get

$$ (2.23) \qquad \frac{\mu + \lambda}{2} \int_0^t \int_{\mathbb{R}^2} (\operatorname{div} u)^3 - \frac{\mu + \lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} F (\operatorname{div} u)^2 = - \frac{\mu + \lambda}{2\nu^2} \int_0^t \int_{\mathbb{R}^2} (F^2 - G^2) \operatorname{div} u. $$

**Step 2: Estimates for the integrals in terms of $E_0, E_0^\nu$ and $L_{t,x}^4$-norms of $(\nabla u, G, F)$.** We are ready to estimate all the terms above.
By use of the $L^q(\mathbb{R}^d)$, $q \in (1, \infty)$, $d \ge 2$-boundedness of Riesz operators, the following estimates follow immediately

$$ (2.17) \qquad \|\nabla u\|_{L^q(\mathbb{R}^d)} \le C(q, d) (\|\operatorname{div} u\|_{L^q(\mathbb{R}^d)} + \|\operatorname{curl} u\|_{L^q(\mathbb{R}^d)}), $$

$$ (2.18) \qquad \|\nabla F\|_{L^q(\mathbb{R}^d)} + \mu \|\nabla \operatorname{curl} u\|_{L^q(\mathbb{R}^d)} \le C(q, d) \|\rho \dot{u}\|_{L^q(\mathbb{R}^d)}. $$

We now recall the compensated result by Coifman, Lions, Meyer, Semmes [8] in dimension two.

**Proposition 2.8** (Coifman-Lions-Meyer-Semmes' Estimate for $d=2$). Let $v, w \in \dot{H}^1(\mathbb{R}^2; \mathbb{R})$ be two functions and let us define

$$ g = \det \begin{pmatrix} \partial_1 v & \partial_2 v \\ \partial_1 w & \partial_2 w \end{pmatrix}. $$

Then $g$ belongs to the Hardy space $H^1(\mathbb{R}^2)$, whence for all $f \in BMO(\mathbb{R}^2)$ we have the estimate

$$ \left| \int_{\mathbb{R}^2} f(x)g(x)dx \right| \le \|f\|_{BMO(\mathbb{R}^2)} \|\nabla v\|_{L^2(\mathbb{R}^2)} \|\nabla w\|_{L^2(\mathbb{R}^2)}. $$

In particular, since $\dot{H}^1(\mathbb{R}^2) \hookrightarrow BMO(\mathbb{R}^2)$, for all $f \in \dot{H}^1(\mathbb{R}^2)$ we have:

$$ \left| \int_{\mathbb{R}^2} f(x)g(x)dx \right| \le \|\nabla f\|_{L^2(\mathbb{R}^2)} \|\nabla v\|_{L^2(\mathbb{R}^2)} \|\nabla w\|_{L^2(\mathbb{R}^2)}. $$

It helps in the proof, by virtue of the following equalities in dimension two

$$ (2.19) \qquad \nabla u^j \cdot \nabla u^k \partial_k u^j = \operatorname{div} u\{|\nabla u|^2 - \det(\nabla u)\}, \qquad \nabla u^l \cdot \partial_l u = (\operatorname{div} u)^2 - 2 \det(\nabla u). $$

Here and in the following we use Einstein's summation convention for repeated indices, unless otherwise claimed.

### 2.2.2. Proof of Lemma 2.1 for $d=2$. This paragraph is devoted to obtaining bounds for functionals $A_1$ and $A_2$ as defined in (1.16) for $d=2$, provided with bounded density function (2.2). The constants in the following estimates may depend on the viscosity coefficient $\mu$ and the lower bound $\nu$ for $\nu$, while not on the viscosity coefficient $\nu$ which will be chosen to be big.

**Proof of (2.3).** The functional $A_1$ arises while using $\dot{u}$ as a test functional in the weak formulation of the momentum equation (1.1)$_2$. By doing so, one obtains the following equality:

$$ \begin{aligned} A_1(t) = \frac{\mu}{2} \|\nabla u_0\|_{L^2(\mathbb{R}^d)}^2 &+ \frac{\mu + \lambda}{2} \|\operatorname{div} u_0\|_{L^2(\mathbb{R}^d)}^2 - \mu \int_0^t \int_{\mathbb{R}^d} \nabla u^j \cdot \nabla u^k \partial_k u^j + \frac{\mu}{2} \int_0^t \int_{\mathbb{R}^d} |\nabla u|^2 \operatorname{div} u \\ &+ \frac{\mu + \lambda}{2} \int_0^t \int_{\mathbb{R}^d} (\operatorname{div} u)^3 - (\mu + \lambda) \int_0^t \int_{\mathbb{R}^d} \operatorname{div} u \nabla u^l \cdot \partial_l u + \int_0^t \int_{\mathbb{R}^d} \nabla u^l \cdot \partial_l u G \\ (2.20) \qquad &+ \int_{\mathbb{R}^d} \operatorname{div} u(s)G(s) \Big|_{s=0}^{s=t} + \int_0^t \int_{\mathbb{R}^d} (\rho P'(\rho) - P(\rho) + \tilde{P})(\operatorname{div} u)^2. \end{aligned} $$

**Step 1: Reformulation of the energy equality.** In the following lines, we will reformulate the terms appearing in the right hand side above by use of (2.16) and (2.19).

By (2.16) and (2.19), the sum of the third and the fourth terms on the right hand side of (2.20) can be reduced as follows:

$$ \begin{aligned} & \mu \int_0^t \int_{\mathbb{R}^2} \operatorname{div} u \left[ \det(\nabla u) - \frac{1}{2} |\nabla u|^2 \right] \\ (2.21) \qquad & = \frac{\mu}{\nu} \int_0^t \int_{\mathbb{R}^2} F \det(\nabla u) - \frac{\mu}{2\nu} \int_0^t \int_{\mathbb{R}^2} F |\nabla u|^2 + \frac{\mu}{\nu} \int_0^t \int_{\mathbb{R}^2} G \left[ \det(\nabla u) - \frac{1}{2} |\nabla u|^2 \right], \end{aligned} $$
• With the help of Proposition 2.8 and (2.18), the first terms of the right hand side of (2.21) and (2.22) can be estimated, as follows:

$$ (2.24) \qquad \begin{aligned} \left| \frac{3\mu + 2\lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} F \det(\nabla u) \right| &\le C \int_0^t \|\nabla F\|_{L^2(\mathbb{R}^2)} \|\nabla u\|_{L^2(\mathbb{R}^2)}^2 \\ &= C \int_0^t \|\rho \dot{u}\|_{L^2(\mathbb{R}^2)} \|\nabla u\|_{L^2(\mathbb{R}^2)}^2 \\ &\le \eta \int_0^t \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C \rho^*}{4\eta} \int_0^t \|\nabla u\|_{L^2(\mathbb{R}^2)}^4, \quad \text{by Young's inequality,} \end{aligned} $$

for some $\eta > 0$ small enough to be determined later.

• By the energy balance (2.1) and the upper bound $\rho^*$ for the density (2.2), the last terms of (2.20) and the terms involving the pressure deviation $G$ in (2.21) and (2.22) can be bounded as follows

$$ (2.25) \qquad \left| \int_0^t \int_{\mathbb{R}^d} (\rho P'(\rho) - P(\rho) + \tilde{P})(\mathrm{div} \, u)^2 \right| + \frac{\mu}{\nu} \left| \int_0^t \int_{\mathbb{R}^2} G \left[ \det(\nabla u) - \frac{1}{2} |\nabla u|^2 \right] \right| + \frac{\mu}{\nu} \left| \int_0^t \int_{\mathbb{R}^2} \nabla u^l \partial_l u G \right| \le \frac{C^*}{\nu} E_0. $$

• Next, the middle term of (2.21) is:

$$ (2.26) \qquad \begin{aligned} \left| \frac{\mu}{2\nu} \int_0^t \int_{\mathbb{R}^2} F |\nabla u|^2 \right| &\le \frac{C}{\nu} \int_0^t \|F\|_{L^4(\mathbb{R}^2)} \|\nabla u\|_{L^4(\mathbb{R}^2)} \|\nabla u\|_{L^2(\mathbb{R}^2)} \\ &\le C \int_0^t \left[ \frac{1}{\nu^{5/2}} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^{3/2}} \|\nabla u\|_{L^4(\mathbb{R}^2)}^4 \right] + C \int_0^t \|\nabla u\|_{L^2(\mathbb{R}^2)}^2. \end{aligned} $$

• The term in (2.23) can be estimated as follows

$$ (2.27) \qquad \left| \frac{\mu + \lambda}{2\nu^2} \int_0^t \int_{\mathbb{R}^2} (F^2 - G^2) \mathrm{div} \, u \right| \le \frac{C}{\nu^3} \int_0^t \left[ \|F\|_{L^4(\mathbb{R}^2)}^4 + \|G\|_{L^4(\mathbb{R}^2)}^4 \right] + C\nu \int_0^t \|\mathrm{div} \, u\|_{L^2(\mathbb{R}^2)}^2. $$

• It only remains the first term in the last line of (2.20) which can be bounded as follows:

$$ \left| \int_{\mathbb{R}^d} \mathrm{div} \, u(s) G(s) \right|_{s=0}^{s=t} \le \eta \nu \|\mathrm{div} \, u(t)\|_{L^2(\mathbb{R}^2)}^2 + \frac{C}{4\eta\nu} \|G(t)\|_{L^2(\mathbb{R}^2)}^2 \\ + \frac{C}{\nu} \|G_0\|_{L^2(\mathbb{R}^2)}^2 + C\nu \|\mathrm{div} \, u_0\|_{L^2(\mathbb{R}^2)}^2, $$

with the second term of the right hand side controlled by the initial energy as in (2.14).

We combine all of these estimates and we choose $\eta$ small in order to obtain the following:

$$ A_1(t) \le C \left(1 + \frac{C^*}{\nu}\right) E_0^\nu + C \rho^* \int_0^t \|\nabla u\|_{L^2(\mathbb{R}^2)}^4 + \frac{C}{\nu^{3/2}} \int_0^t \left[ \|\nabla u\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^{3/2}} \|G\|_{L^4(\mathbb{R}^2)}^4 \right] $$

where $E_0^\nu$ is given in (1.17). Hence Gronwall Lemma yields:

$$ (2.28) \qquad A_1(t) \le \left[ C^* E_0^\nu + \frac{C}{\nu^{3/2}} \int_0^t \left( \|\nabla u\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^{3/2}} \|G\|_{L^4(\mathbb{R}^2)}^4 \right) \right] \exp(C^* E_0). $$

**Step 3: Final estimates.** The next step is devoted to obtaining estimate for the $L^4((0, t) \times \mathbb{R}^2)$ norm of the velocity gradient $\nabla u$, the pressure deviation $G$ and the effective flux $F$.

• $L^4$-Estimate for $G$. Recall (2.15) with $l = 3$:

$$ (2.29) \qquad \frac{1}{\nu} \|G\|_{L^4((0,t) \times \mathbb{R}^d)}^4 \le C^* E_0 + \frac{C}{\nu} \|F\|_{L^4((0,t) \times \mathbb{R}^d)}^4. $$
and similarly, the sum of the second and the third terms of the second line of (2.20) read

$$ (2.22) \qquad -\frac{\mu+\lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} F \nabla u^l \cdot \partial_l u + \frac{\mu}{\nu} \int_0^t \int_{\mathbb{R}^2} \nabla u^l \cdot \partial_l u G \\ = 2\frac{\mu+\lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} F \det(\nabla u) - \frac{\mu+\lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} F(\operatorname{div} u)^2 + \frac{\mu}{\nu} \int_0^t \int_{\mathbb{R}^2} \nabla u^l \cdot \partial_l u G. $$

Now we pack the first term of the second line of (2.20) and the middle term in the above (2.22) and use (2.16) to get

$$ (2.23) \qquad \frac{\mu+\lambda}{2} \int_0^t \int_{\mathbb{R}^2} (\operatorname{div} u)^3 - \frac{\mu+\lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} F(\operatorname{div} u)^2 = -\frac{\mu+\lambda}{2\nu^2} \int_0^t \int_{\mathbb{R}^2} (F^2 - G^2) \operatorname{div} u. $$

**Step 2: Estimates for the integrals in terms of $E_0, E_0^\nu$ and $L_{t,x}^4$-norms of $(\nabla u, G, F)$**. We are ready to estimate all the terms above.

• With the help of Proposition 2.8 and (2.18), the first terms of the right hand side of (2.21) and (2.22) can be estimated, as follows:

$$ (2.24) \qquad \left| \frac{3\mu + 2\lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} F \det(\nabla u) \right| \le C \int_0^t \|\nabla F\|_{L^2(\mathbb{R}^2)} \|\nabla u\|_{L^2(\mathbb{R}^2)}^2 \\ = C \int_0^t \|\rho \dot{u}\|_{L^2(\mathbb{R}^2)} \|\nabla u\|_{L^2(\mathbb{R}^2)}^2 \\ \le \eta \int_0^t \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C \rho^*}{4\eta} \int_0^t \|\nabla u\|_{L^2(\mathbb{R}^2)}^4, \quad \text{by Young's inequality,} \\ \text{for some } \eta > 0 \text{ small enough to be determined later.} $$

• By the energy balance (2.1) and the upper bound $\rho^*$ for the density (2.2), the last terms of (2.20) and the terms involving the pressure deviation $G$ in (2.21) and (2.22) can be bounded as follows

$$ (2.25) \qquad \left| \int_0^t \int_{\mathbb{R}^d} (\rho P'(\rho) - P(\rho) + \tilde{P})(\operatorname{div} u)^2 \right| + \frac{\mu}{\nu} \left| \int_0^t \int_{\mathbb{R}^2} G \left[ \det(\nabla u) - \frac{1}{2} |\nabla u|^2 \right] \right| \\ \qquad \qquad + \frac{\mu}{\nu} \left| \int_0^t \int_{\mathbb{R}^2} \nabla u^l \partial_l u G \right| \le \frac{C^*}{\nu} E_0. $$

• Next, the middle term of (2.21) is:

$$ (2.26) \qquad \left| \frac{\mu}{2\nu} \int_0^t \int_{\mathbb{R}^2} F |\nabla u|^2 \right| \le \frac{C}{\nu} \int_0^t \|F\|_{L^4(\mathbb{R}^2)} \|\nabla u\|_{L^4(\mathbb{R}^2)} \|\nabla u\|_{L^2(\mathbb{R}^2)} \\ \le C \int_0^t \left[ \frac{1}{\nu^{5/2}} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^{3/2}} \|\nabla u\|_{L^4(\mathbb{R}^2)}^4 \right] + C \int_0^t \|\nabla u\|_{L^2(\mathbb{R}^2)}^2. $$

• The term in (2.23) can be estimated as follows

$$ (2.27) \qquad \left| \frac{\mu + \lambda}{2\nu^2} \int_0^t \int_{\mathbb{R}^2} (F^2 - G^2) \operatorname{div} u \right| \le \frac{C}{\nu^3} \int_0^t \left[ \|F\|_{L^4(\mathbb{R}^2)}^4 + \|G\|_{L^4(\mathbb{R}^2)}^4 \right] + C\nu \int_0^t \|\operatorname{div} u\|_{L^2(\mathbb{R}^2)}^2. $$

• It only remains the first term in the last line of (2.20) which can be bounded as follows:

$$ \left| \int_{\mathbb{R}^d} \operatorname{div} u(s) G(s) \right|_{s=0}^{s=t} \le \eta \nu \|\operatorname{div} u(t)\|_{L^2(\mathbb{R}^2)}^2 + \frac{C}{4\eta\nu} \|G(t)\|_{L^2(\mathbb{R}^2)}^2 \\ + \frac{C}{\nu} \|G_0\|_{L^2(\mathbb{R}^2)}^2 + C\nu \|\operatorname{div} u_0\|_{L^2(\mathbb{R}^2)}^2, $$

with the second term of the right hand side controlled by the initial energy as in (2.14).

We combine all of these estimates and we choose $\eta$ small in order to obtain the following:

$$ A_1(t) \le C \left(1 + \frac{C^*}{\nu}\right) E_0^\nu + C \rho^* \int_0^t \|\nabla u\|_{L^2(\mathbb{R}^2)}^4 + \frac{C}{\nu^{3/2}} \int_0^t \left[ \|\nabla u\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^{3/2}} \|G\|_{L^4(\mathbb{R}^2)}^4 \right] $$
* $L^4$-Estimate for $F$. The $L^4((0,t) \times \mathbb{R}^2)$-norm of the effective flux $F$ follows from Gagliarodo-Nirenberg's inequality $\|f\|_{L^4(\mathbb{R}^2)}^2 \lesssim \|f\|_{L^2(\mathbb{R}^2)} \|\nabla f\|_{L^2(\mathbb{R}^2)}$ and (2.18) as

$$ \|F\|_{L^4(\mathbb{R}^2)} \le C \|F\|_{L^2(\mathbb{R}^2)}^{\frac{1}{2}} \|\nabla F\|_{L^2(\mathbb{R}^2)}^{\frac{1}{2}} \le C \|F\|_{L^2(\mathbb{R}^2)}^{\frac{1}{2}} \|\rho \dot{u}\|_{L^2(\mathbb{R}^2)}^{\frac{1}{2}}, $$

which can be bounded further by virtue of (2.16) and the definition of $E(t)$, $A_1(t)$ by

$$ \begin{aligned} \|F\|_{L^4((0,t) \times \mathbb{R}^2)}^4 &\le C \int_0^t (\nu \|\operatorname{div} u\|_{L^2(\mathbb{R}^2)} + \|G\|_{L^2(\mathbb{R}^2)})^2 \|\rho \dot{u}\|_{L^2(\mathbb{R}^2)}^2 \\ &\le C^*(\nu A_1(t) + E_0) A_1(t). \end{aligned} $$

* $L^4$-Estimate for $\nabla u$. Similar as above for $F$, we have $L^4$-Estimate for $\operatorname{curl} u$:

$$ \|\operatorname{curl} u\|_{L^4((0,t) \times \mathbb{R}^2)} \le C \|\operatorname{curl} u\|_{L^\infty((0,t); L^2(\mathbb{R}^2))}^{\frac{1}{2}} \|\rho \dot{u}\|_{L^2((0,t) \times \mathbb{R}^2)}^{\frac{1}{2}} \le C A_1(t)^{\frac{1}{2}}. $$

Hence, by use of (2.16)-(2.17), the following inequality holds true

$$ \begin{aligned} \|\nabla u\|_{L^4((0,t) \times \mathbb{R}^2)}^4 &\le C \left( \|\operatorname{div} u\|_{L^4((0,t) \times \mathbb{R}^2)}^4 + \|\operatorname{curl} u\|_{L^4((0,t) \times \mathbb{R}^2)}^4 \right) \\ &\le \frac{C}{\nu^4} \left( \|F\|_{L^4((0,t) \times \mathbb{R}^2)}^4 + \|G\|_{L^4((0,t) \times \mathbb{R}^2)}^4 \right) + C \|\operatorname{curl} u\|_{L^4((0,t) \times \mathbb{R}^2)}^4 \\ &\le \frac{C^*}{\nu^3} E_0 + C^* A_1(t) (E_0 + A_1(t)). \end{aligned} $$

Finally, we go back to (2.28) and we have (2.3). □

*Proof of (2.4).* We turn to providing bound for the second functional $A_2$ for $d=2$. For this purpose, by (2.16), we rewrite the momentum equation (1.1)$_2$ as follows

$$ (2.30) \qquad \rho \dot{u} = \mu \Delta u + \frac{\mu + \lambda}{\nu} \nabla F - \frac{\mu}{\nu} \nabla G. $$

We apply the operator $\partial_t \cdot + \operatorname{div}(\cdot u)$ to (2.30) and we obtain the following equation for the material derivative of the velocity:

$$ \begin{aligned} (2.31) \quad &\partial_t(\rho\dot{u}^j) + \operatorname{div}(\rho\dot{u}^j u) - \mu\Delta\dot{u}^j - \frac{\mu+\lambda}{\nu}\partial_j\dot{F} \\ &= -\mu\partial_k(\nabla u^j \cdot \partial_k u) + \mu\partial_k(\partial_k u^j \operatorname{div} u) - \mu\operatorname{div}(\partial_k u^j \partial_k u) \\ &\quad + \frac{\mu+\lambda}{\nu}\partial_j(F \operatorname{div} u) - \frac{\mu+\lambda}{\nu}\operatorname{div}(F\partial_j u) \\ &\quad + \frac{\mu}{\nu}\partial_j\left((\rho P'(\rho) - P(\rho) + \tilde{P}) \operatorname{div} u\right) + \frac{\mu}{\nu}\operatorname{div}(\partial_j u G), \quad j=1, \dots, d. \end{aligned} $$

**Step 1: Formulation of the energy equality.** To obtain the functional $A_2$, it suffices to multiply the equation above by $\sigma \dot{u}^j$, with $\sigma = \sigma(t) = \min\{1, t\}$, sum up $j$, and integrate it in time and space. The most delicate term is $-\frac{\mu+\lambda}{\nu}\partial_j \dot{F}$ on the left hand side of (2.31), which gives

$$ -\frac{\mu+\lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} \sigma \dot{u}^j \partial_j \dot{F} = \frac{\mu+\lambda}{\nu} \int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \operatorname{div} \dot{u}. $$

We first focus on this integral for a while. Applying material derivative to (2.16) gives

$$ (2.32) \qquad \operatorname{div} \dot{u} = \frac{1}{\nu} \left( \dot{F} - \rho P'(\rho) \operatorname{div} u \right) + \nabla u^k \cdot \partial_k u, $$

and hence

$$ (2.33) \qquad -\frac{\mu+\lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} \sigma \dot{u}^j \partial_j \dot{F} = \frac{\mu+\lambda}{\nu^2} \int_0^t \sigma \|\dot{F}\|_{L^2(\mathbb{R}^2)}^2 + \frac{\mu+\lambda}{\nu} \int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \nabla u^k \cdot \partial_k u \\ -\frac{\mu+\lambda}{\nu^2} \int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \rho P'(\rho) \operatorname{div} u. $$
where $E_0^\nu$ is given in (1.18). Hence Grönwall's lemma yields:

$$ (2.28) \qquad \mathcal{A}_1(t) \le \left[ C^* E_0^\nu + \frac{C}{\nu^{3/2}} \int_0^t \left( \|\nabla u\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^{3/2}} \|G\|_{L^4(\mathbb{R}^2)}^4 \right) \right] \exp(C^* E_0). $$

**Step 3: Final estimates.** The next step is devoted to obtaining estimate for the $L^4((0,t) \times \mathbb{R}^2)$ norm of the velocity gradient $\nabla u$, the pressure deviation $G$ and the effective flux $F$.

* $L^4$-Estimate for $G$. Recall (2.15) with $l=3$:

$$ (2.29) \qquad \frac{1}{\nu} \|G\|_{L^4((0,t) \times \mathbb{R}^d)}^4 \le C^* E_0 + \frac{C}{\nu} \|F\|_{L^4((0,t) \times \mathbb{R}^d)}^4. $$

* $L^4$-Estimate for $F$. The $L^4((0,t) \times \mathbb{R}^2)$-norm of the effective flux $F$ follows from Gagliarodo-Nirenberg's inequality $\|f\|_{L^4(\mathbb{R}^2)}^2 \lesssim \|f\|_{L^2(\mathbb{R}^2)} \|\nabla f\|_{L^2(\mathbb{R}^2)}$ and (2.18) as

$$ \|F\|_{L^4(\mathbb{R}^2)} \le C \|F\|_{L^2(\mathbb{R}^2)}^{\frac{1}{2}} \|\nabla F\|_{L^2(\mathbb{R}^2)}^{\frac{1}{2}} \le C \|F\|_{L^2(\mathbb{R}^2)}^{\frac{1}{2}} \|\rho \dot{u}\|_{L^2(\mathbb{R}^2)}^{\frac{1}{2}}, $$

which can be bounded further by virtue of (2.16) and the definition of $E(t)$, $\mathcal{A}_1(t)$ by

$$ \begin{aligned} \|F\|_{L^4((0,t) \times \mathbb{R}^2)}^4 &\le C \int_0^t (\nu \|\operatorname{div} u\|_{L^2(\mathbb{R}^2)} + \|G\|_{L^2(\mathbb{R}^2)})^2 \|\rho \dot{u}\|_{L^2(\mathbb{R}^2)}^2 \\ &\le C^*(\nu \mathcal{A}_1(t) + E_0) \mathcal{A}_1(t). \end{aligned} $$

* $L^4$-Estimate for $\nabla u$. Similar as above for $F$, we have $L^4$-Estimate for $\operatorname{curl} u$:

$$ \|\operatorname{curl} u\|_{L^4((0,t) \times \mathbb{R}^2)} \le C \|\operatorname{curl} u\|_{L^\infty((0,t);L^2(\mathbb{R}^2))}^{\frac{1}{2}} \|\rho \dot{u}\|_{L^2((0,t) \times \mathbb{R}^2)}^{\frac{1}{2}} \le C \mathcal{A}_1(t)^{\frac{1}{2}}. $$

Hence, by use of (2.16)-(2.17), the following inequality holds true

$$ \begin{aligned} \|\nabla u\|_{L^4((0,t) \times \mathbb{R}^2)}^4 &\le C \left( \|\operatorname{div} u\|_{L^4((0,t) \times \mathbb{R}^2)}^4 + \|\operatorname{curl} u\|_{L^4((0,t) \times \mathbb{R}^2)}^4 \right) \\ &\le \frac{C}{\nu^4} \left( \|F\|_{L^4((0,t) \times \mathbb{R}^2)}^4 + \|G\|_{L^4((0,t) \times \mathbb{R}^2)}^4 \right) + C \|\operatorname{curl} u\|_{L^4((0,t) \times \mathbb{R}^2)}^4 \\ &\le \frac{C^*}{\nu^3} E_0 + C^* \mathcal{A}_1(t) (E_0 + \mathcal{A}_1(t)). \end{aligned} $$

Finally, we go back to (2.28) and we have (2.3). $\square$

**Proof of (2.4).** We turn to providing bound for the second functional $\mathcal{A}_2$ for $d=2$. For this purpose, by (2.16), we rewrite the momentum equation (1.1)$_2$ as follows

$$ (2.30) \qquad \rho \dot{u} = \mu \Delta u + \frac{\mu + \lambda}{\nu} \nabla F - \frac{\mu}{\nu} \nabla G. $$

We apply the operator $\partial_t \cdot + \operatorname{div}(\cdot u)$ to (2.30) and we obtain the following equation for the material derivative of the velocity:

$$ \begin{aligned} (2.31) \quad &\partial_t(\rho\dot{u}^j) + \operatorname{div}(\rho\dot{u}^j u) - \mu\Delta\dot{u}^j - \frac{\mu+\lambda}{\nu}\partial_j\dot{F} \\ &= -\mu\partial_k(\nabla u^j \cdot \partial_k u) + \mu\partial_k(\partial_k u^j \operatorname{div} u) - \mu\operatorname{div}(\partial_k u^j \partial_k u) \\ &\quad + \frac{\mu+\lambda}{\nu}\partial_j(F \operatorname{div} u) - \frac{\mu+\lambda}{\nu}\operatorname{div}(F\partial_j u) \\ &\quad + \frac{\mu}{\nu}\partial_j\left((\rho P'(\rho) - P(\rho) + \tilde{P}) \operatorname{div} u\right) + \frac{\mu}{\nu}\operatorname{div}(\partial_j u G), \quad j=1, \dots, d. \end{aligned} $$

**Step 1: Formulation of the energy equality.** To obtain the functional $\mathcal{A}_2$, it suffices to multiply the equation above by $\sigma\dot{u}^j$, with $\sigma = \sigma(t) = \min\{1, t\}$, sum up $j$, and integrate it in time and space. The most delicate term is $-\frac{\mu+\lambda}{\nu}\partial_j\dot{F}$ on the left hand side of (2.31), which gives

$$ -\frac{\mu+\lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} \sigma \dot{u}^j \partial_j \dot{F} = \frac{\mu+\lambda}{\nu} \int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \operatorname{div} \dot{u}. $$
To conclude, testing (2.31) by $\sigma \dot{u}$ implies

$$ (2.34) \quad \frac{1}{2} \sigma(t) \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \mu \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{\mu + \lambda}{\nu^2} \int_0^t \sigma \|\dot{F}\|_{L^2(\mathbb{R}^2)}^2 = \frac{1}{2} \int_0^{\sigma(t)} \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \sum_{k=1}^4 I_k, $$

where

$$ I_1 = -\frac{\mu + \lambda}{\nu} \int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \nabla u^k \cdot \partial_k u + \frac{\mu + \lambda}{\nu^2} \int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \rho P'(\rho) \operatorname{div} u, $$

$$ I_2 = \mu \int_0^t \sigma \int_{\mathbb{R}^2} (\nabla u^j \cdot \partial_k u \partial_k \dot{u}^j - \partial_k u^j \operatorname{div} u \partial_k \dot{u}^j + \partial_k u^j \partial_k u \cdot \nabla \dot{u}^j), $$

$$ I_3 = \frac{\mu + \lambda}{\nu} \int_0^t \sigma \int_{\mathbb{R}^2} (-F \operatorname{div} u \operatorname{div} \dot{u} + F \partial_j u^k \partial_k \dot{u}^j), $$

$$ I_4 = -\frac{\mu}{\nu} \int_0^t \sigma \int_{\mathbb{R}^2} \left( \operatorname{div} u \operatorname{div} \dot{u} (\rho P'(\rho) - P(\rho) + \tilde{P}) + \partial_j u^k \partial_k \dot{u}^j G \right). $$

**Step 2: Estimate for $I_1$.** We focus first on the first integral in $I_1$, which can be reformulated by integration by parts (noticing $\sigma(0) = 0$) as

$$ (2.35) \quad \int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \nabla u^k \cdot \partial_k u = \sigma(t) \int_{\mathbb{R}^2} F(t) \nabla u^k \cdot \partial_k u(t) - \int_0^{\sigma(t)} \int_{\mathbb{R}^2} F \nabla u^k \cdot \partial_k u - 2 \int_0^t \sigma \int_{\mathbb{R}^2} F \partial_k u \cdot \nabla \dot{u}^k \\ + 2 \int_0^t \sigma \int_{\mathbb{R}^2} F \partial_k u \cdot \nabla u^l \partial_l u^k - \int_0^t \sigma \int_{\mathbb{R}^2} F \operatorname{div} u \nabla u^k \cdot \partial_k u. $$

The first term of the right hand side in (2.35) above can be estimated similarly as for the derivation of (2.24) and (2.27), using the equality $\nabla u^l \cdot \partial_l u = (\operatorname{div} u)^2 - 2 \operatorname{det}(\nabla u)$ in (2.19):

$$ \begin{aligned} \left| \sigma(t) \int_{\mathbb{R}^2} F(t) (\nabla u^k \cdot \partial_k u)(t) \right| &\le C \sigma(t) \|\rho \dot{u}\|_{L^2(\mathbb{R}^2)} \|\nabla u\|_{L^2(\mathbb{R}^2)}^2 + \sigma(t) \frac{1}{\nu^2} \left| \int_{\mathbb{R}^2} F(t)(F(t) + G(t))^2 \right| \\ &\le \eta \sigma(t) \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \sigma(t) \frac{C^2 \rho^*}{4\eta} \|\nabla u\|_{L^2(\mathbb{R}^2)}^4 \\ &\quad + \sigma(t) \frac{C}{\nu^2} \left( \|F(t)\|_{L^3(\mathbb{R}^2)}^3 + \|G(t)\|_{L^3(\mathbb{R}^2)}^3 \right). \end{aligned} $$

Exactly as in the derivation of (2.24) and (2.27), the second integral of the right hand side of (2.35) can be estimated as follows:

$$ (2.36) \quad \left| \int_0^{\sigma(t)} \int_{\mathbb{R}^2} F \nabla u^k \cdot \partial_k u \right| \le C E_0 + C \int_0^{\sigma(t)} \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + C \rho^* \int_0^{\sigma(t)} \|\nabla u\|_{L^2(\mathbb{R}^2)}^4 \\ + \frac{C}{\nu^3} \int_0^{\sigma(t)} \left[ \|F\|_{L^4(\mathbb{R}^2)}^4 + \|G\|_{L^4(\mathbb{R}^2)}^4 \right]. $$

In order to estimate the third term of the right hand side of (2.35), we write:

$$ (2.37) \quad \partial_j u^k \partial_k \dot{u}^j = \operatorname{div} u \operatorname{div} \dot{u} - (\partial_1 u^1 \partial_2 \dot{u}^2 - \partial_2 u^1 \partial_1 \dot{u}^2) - (\partial_2 u^2 \partial_1 \dot{u}^1 - \partial_1 u^2 \partial_2 \dot{u}^1). $$

in such a way that, after making use of the compensated result Proposition 2.8 and Young's inequality, we have:

$$ (2.38) \quad \left| \int_0^t \sigma \int_{\mathbb{R}^2} F \partial_k u \cdot \nabla \dot{u}^k \right| \le \eta \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C \rho^*}{\eta} \int_0^t \sigma \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 \|\nabla u\|_{L^2(\mathbb{R}^2)}^2 \\ + \frac{C}{\eta} \int_0^t \left[ \frac{1}{\nu^2} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2} \|G\|_{L^4(\mathbb{R}^2)}^4 \right]. $$
We first focus on this integral for a while. Applying material derivative to (2.16) gives

$$
(2.32) \qquad \operatorname{div} \dot{u} = \frac{1}{\nu} \left( \dot{F} - \rho P'(\rho) \operatorname{div} u \right) + \nabla u^k \cdot \partial_k u,
$$

and hence

$$
\begin{equation}
\begin{aligned}
-\frac{\mu + \lambda}{\nu} \int_0^t \int_{\mathbb{R}^2} \sigma \dot{u}^j \partial_j \dot{F} &= \frac{\mu + \lambda}{\nu^2} \int_0^t \sigma \|\dot{F}\|_{L^2(\mathbb{R}^2)}^2 + \frac{\mu + \lambda}{\nu} \int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \nabla u^k \cdot \partial_k u \\
&\quad - \frac{\mu + \lambda}{\nu^2} \int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \rho P'(\rho) \operatorname{div} u.
\end{aligned}
\tag{2.33}
\end{equation}
$$

To conclude, testing (2.31) by $\sigma \dot{u}$ implies

$$
(2.34) \quad \frac{1}{2} \sigma(t) \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \mu \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{\mu + \lambda}{\nu^2} \int_0^t \sigma \|\dot{F}\|_{L^2(\mathbb{R}^2)}^2 = \frac{1}{2} \int_0^{\sigma(t)} \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \sum_{k=1}^4 I_k,
$$

where

$$
I_1 = - \frac{\mu + \lambda}{\nu} \int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \nabla u^k \cdot \partial_k u + \frac{\mu + \lambda}{\nu^2} \int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \rho P'(\rho) \operatorname{div} u,
$$

$$
I_2 = \mu \int_0^t \sigma \int_{\mathbb{R}^2} \left( \nabla u^j \cdot \partial_k u \partial_k \dot{u}^j - \partial_k u^j \operatorname{div} u \partial_k \dot{u}^j + \partial_k u^j \partial_k u \cdot \nabla \dot{u}^j \right),
$$

$$
I_3 = \frac{\mu + \lambda}{\nu} \int_0^t \sigma \int_{\mathbb{R}^2} \left( -F \operatorname{div} u \operatorname{div} \dot{u} + F \partial_j u^k \partial_k \dot{u}^j \right),
$$

$$
I_4 = - \frac{\mu}{\nu} \int_0^t \sigma \int_{\mathbb{R}^2} \left( \operatorname{div} u \operatorname{div} \dot{u} (\rho P'(\rho) - P(\rho) + \tilde{P}) + \partial_j u^k \partial_k \dot{u}^j G \right).
$$

**Step 2: Estimate for $I_1$.** We focus first on the first integral in $I_1$, which can be reformulated by integration by parts (noticing $\sigma(0) = 0$) as

$$
\begin{equation}
\begin{aligned}
\int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \nabla u^k \cdot \partial_k u = \sigma(t) \int_{\mathbb{R}^2} F(t) \nabla u^k \cdot \partial_k u(t) & - \int_0^{\sigma(t)} \int_{\mathbb{R}^2} F \nabla u^k \cdot \partial_k u - 2 \int_0^t \sigma \int_{\mathbb{R}^2} F \partial_k u \cdot \nabla \dot{u}^k \\
& + 2 \int_0^t \sigma \int_{\mathbb{R}^2} F \partial_k u \cdot \nabla u^l \partial_l u^k - \int_0^t \sigma \int_{\mathbb{R}^2} F \operatorname{div} u \nabla u^k \cdot \partial_k u.
\end{aligned}
\tag{2.35}
\end{equation}
$$

The first term of the right hand side in (2.35) above can be estimated similarly as for the derivation of (2.24) and (2.27), using the equality $\nabla u^l \cdot \partial_l u = (\operatorname{div} u)^2 - 2 \operatorname{det}(\nabla u)$ in (2.19):

$$
\begin{align*}
\left| \sigma(t) \int_{\mathbb{R}^2} F(t) (\nabla u^k \cdot \partial_k u)(t) \right| &\le C \sigma(t) \|\rho \dot{u}\|_{L^2(\mathbb{R}^2)} \|\nabla u\|_{L^2(\mathbb{R}^2)}^2 + \sigma(t) \frac{1}{\nu^2} \left| \int_{\mathbb{R}^2} F(t)(F(t) + G(t))^2 \right| \\
&\le \eta \sigma(t) \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \sigma(t) \frac{C^2 \rho^*}{4\eta} \|\nabla u\|_{L^2(\mathbb{R}^2)}^4 \\
&\quad + \sigma(t) \frac{C}{\nu^2} \left( \|F(t)\|_{L^3(\mathbb{R}^2)}^3 + \|G(t)\|_{L^3(\mathbb{R}^2)}^3 \right).
\end{align*}
$$

Exactly as in the derivation of (2.24) and (2.27), the second integral of the right hand side of (2.35) can
be estimated as follows:

$$
(2.36) \quad \left| \int_0^{\sigma(t)} \int_{\mathbb{R}^2} F \nabla u^k \cdot \partial_k u \right| \le C E_0 + C \int_0^{\sigma(t)} \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + C \rho^* \int_0^{\sigma(t)} \|\nabla u\|_{L^2(\mathbb{R}^2)}^4 \\
\qquad + \frac{C}{\nu^3} \int_0^{\sigma(t)} \left[ \|F\|_{L^4(\mathbb{R}^2)}^4 + \|G\|_{L^4(\mathbb{R}^2)}^4 \right].
$$

In order to estimate the third term of the right hand side of (2.35), we write:

$$
(2.37) \qquad \partial_j u^k \partial_k \dot{u}^j = \operatorname{div} u \operatorname{div} \dot{u} - (\partial_1 u^1 \partial_2 \dot{u}^2 - \partial_2 u^1 \partial_1 \dot{u}^2) - (\partial_2 u^2 \partial_1 \dot{u}^1 - \partial_1 u^2 \partial_2 \dot{u}^1)
$$
Thanks to (2.16) and the equality $\partial_k u^l \nabla u^k \cdot \partial_l u = \operatorname{div} u \{(\operatorname{div} u)^2 - 3 \det(\nabla u)\}$ in (2.19), we have the following estimate for the last two terms of (2.35):

$$2 \left| \int_0^t \sigma \int_{\mathbb{R}^2} F \partial_k u \cdot \nabla u^l \partial_l u^k \right| + \left| \int_0^t \sigma \int_{\mathbb{R}^2} F \operatorname{div} u \nabla u^k \cdot \partial_k u \right| \\ \le C \int_0^t \sigma \left[ \|\nabla u\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2} \|G\|_{L^4(\mathbb{R}^2)}^4 \right].$$

Finally, owing to Hölder's and Young's inequalities the second integral in $I_0$ can be estimated as follows:

$$\frac{\mu + \lambda}{\nu^2} \left| \int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \rho P'(\rho) \operatorname{div} u \right| \le \eta \frac{\mu + \lambda}{\nu^2} \int_0^t \sigma \|\dot{F}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C^*}{\nu^2 \eta} E_0.$$

Gathering the above computations, we obtain the following estimate for $I_1$:

$$|I_1| \le C\left(1 + \frac{C^*}{\nu^2 \eta}\right)E_0 + \eta\sigma(t)\|\sqrt{\rho}\dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \eta\int_0^t \sigma\|\nabla\dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \eta\frac{\mu+\lambda}{\nu^2}\int_0^t \sigma\|\dot{F}\|_{L^2(\mathbb{R}^2)}^2 \\ + C\int_0^{\sigma(t)} \|\sqrt{\rho}\dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \sigma(t)\frac{C^*}{\eta}\|\nabla u\|_{L^2(\mathbb{R}^2)}^4 + \frac{C\rho^*}{\eta}\int_0^t \sigma\|\sqrt{\rho}\dot{u}\|_{L^2(\mathbb{R}^2)}^2\|\nabla u\|_{L^2(\mathbb{R}^2)}^2 \\ + C\rho^*\int_0^{\sigma(t)} \|\nabla u\|_{L^2(\mathbb{R}^2)}^4 + \sigma(t)\frac{C}{\nu^2}\left(\|F(t)\|_{L^3(\mathbb{R}^2)}^3 + \|G(t)\|_{L^3(\mathbb{R}^2)}^3\right) \\ + C\int_0^t \sigma\left[\|\nabla u\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\eta\nu^2}\|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\eta\nu^2}\|G\|_{L^4(\mathbb{R}^2)}^4\right].$$

**Step 3: Final estimates.** We now turn to the estimate of the last terms $\sum_{k=2}^4 I_k$ in (2.34). By Young's inequality it is straightforward to get

$$|I_2| \le \eta \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C}{\eta} \int_0^t \sigma \|\nabla u\|_{L^4(\mathbb{R}^2)}^4.$$

Similar as in Step 2 we have

$$|I_3| \le \eta \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C \rho^*}{\eta} \int_0^t \sigma \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 \|\nabla u\|_{L^2(\mathbb{R}^2)}^2 + \frac{C}{\nu^2 \eta} \int_0^t \sigma \left[ \|F\|_{L^4(\mathbb{R}^2)}^4 + \|G\|_{L^4(\mathbb{R}^2)}^4 \right],$$

$$|I_4| \le \eta \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C C^*}{\nu^2 \eta} E_0.$$

Summing up, we have for $\eta$ small enough that

$$\mathcal{A}_2(t) \le C\left(1 + \frac{C^*}{\nu^2}\right)E_0 + C^*(1 + E_0 + \mathcal{A}_1(t))\mathcal{A}_1(t) \\ (2.39) \quad + \sigma(t)\frac{C}{\nu^2}\left(\|F(t)\|_{L^3(\mathbb{R}^2)}^3 + \|G(t)\|_{L^3(\mathbb{R}^2)}^3\right) + C\int_0^t \sigma\left[\|\nabla u\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2}\|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2}\|G\|_{L^4(\mathbb{R}^2)}^4\right].$$

Recall Step 3 in the Proof of (2.3) above, and we get similar $L^4$-Estimates with the time weight $\sigma$. We have the following similar as (2.29)

$$\frac{1}{2\nu} \int_0^t \sigma \|G\|_{L^4(\mathbb{R}^2)}^4 \le C^* E_0 + \frac{C}{\nu} \int_0^t \sigma \|F\|_{L^4(\mathbb{R}^2)}^4,$$

which implies

$$(2.40) \quad \int_0^t \sigma \left[ \|\nabla u\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2} \|G\|_{L^4(\mathbb{R}^2)}^4 \right] \le \frac{C^*}{\nu^3} E_0 + C^*(E_0 + \mathcal{A}_1(t))\mathcal{A}_1(t).$$
in such a way that, after making use of the compensated result Proposition 2.8 and Young's inequality, we have:

$$ (2.38) \quad \left| \int_0^t \sigma \int_{\mathbb{R}^2} F \partial_k u \cdot \nabla \dot{u}^k \right| \le \eta \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C \rho^*}{\eta} \int_0^t \sigma \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 \|\nabla u\|_{L^2(\mathbb{R}^2)}^2 \\ \qquad + \frac{C}{\eta} \int_0^t \left[ \frac{1}{\nu^2} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2} \|G\|_{L^4(\mathbb{R}^2)}^4 \right]. $$

Thanks to (2.16) and the equality $\partial_k u^l \nabla u^k \cdot \partial_l u = \text{div } u \{(\text{div } u)^2 - 3 \det(\nabla u)\}$ in (2.19), we have the following estimate for the last two terms of (2.35):

$$ \begin{aligned} & 2 \left| \int_0^t \sigma \int_{\mathbb{R}^2} F \partial_k u \cdot \nabla u^l \partial_l u^k \right| + \left| \int_0^t \sigma \int_{\mathbb{R}^2} F \text{div } u \nabla u^k \cdot \partial_k u \right| \\ & \le C \int_0^t \sigma \left[ \|\nabla u\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2} \|G\|_{L^4(\mathbb{R}^2)}^4 \right]. \end{aligned} $$

Finally, owing to Hölder's and Young's inequalities the second integral in $I_1$ can be estimated as follows:

$$ \frac{\mu + \lambda}{\nu^2} \left| \int_0^t \sigma \int_{\mathbb{R}^2} \dot{F} \rho P'(\rho) \text{div } u \right| \le \eta \frac{\mu + \lambda}{\nu^2} \int_0^t \sigma \|\dot{F}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C^*}{\nu^2 \eta} E_0. $$

Gathering the above computations, we obtain the following estimate for $I_1$:

$$ \begin{aligned} |I_1| \le & C \left(1 + \frac{C^*}{\nu^2 \eta}\right) E_0 + \eta \sigma(t) \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \eta \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \eta \frac{\mu + \lambda}{\nu^2} \int_0^t \sigma \|\dot{F}\|_{L^2(\mathbb{R}^2)}^2 \\ & + C \int_0^{\sigma(t)} \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \sigma(t) \frac{C^*}{\eta} \|\nabla u\|_{L^2(\mathbb{R}^2)}^4 + \frac{C \rho^*}{\eta} \int_0^t \sigma \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 \|\nabla u\|_{L^2(\mathbb{R}^2)}^2 \\ & + C \rho^* \int_0^{\sigma(t)} \|\nabla u\|_{L^2(\mathbb{R}^2)}^4 + \sigma(t) \frac{C}{\nu^2} \left( \|F(t)\|_{L^3(\mathbb{R}^2)}^3 + \|G(t)\|_{L^3(\mathbb{R}^2)}^3 \right) \\ & + C \int_0^t \sigma \left[ \|\nabla u\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\eta \nu^2} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\eta \nu^2} \|G\|_{L^4(\mathbb{R}^2)}^4 \right]. \end{aligned} $$

**Step 3: Final estimates.** We now turn to the estimate of the last terms $\sum_{k=2}^4 I_k$ in (2.34). By Young's inequality it is straightforward to get

$$ |I_2| \le \eta \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C}{\eta} \int_0^t \sigma \|\nabla u\|_{L^4(\mathbb{R}^2)}^4. $$

Similar as in Step 2 we have

$$ \begin{aligned} |I_3| \le & \eta \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C \rho^*}{\eta} \int_0^t \sigma \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 \|\nabla u\|_{L^2(\mathbb{R}^2)}^2 + \frac{C}{\nu^2 \eta} \int_0^t \sigma \left[ \|F\|_{L^4(\mathbb{R}^2)}^4 + \|G\|_{L^4(\mathbb{R}^2)}^4 \right], \\ |I_4| \le & \eta \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C C^*}{\nu^2 \eta} E_0. \end{aligned} $$

Summing up, we have for $\eta$ small enough that

$$ \begin{aligned} \mathcal{A}_2(t) \le & C \left(1 + \frac{C^*}{\nu^2}\right) E_0 + C^*(1 + E_0 + \mathcal{A}_1(t)) \mathcal{A}_1(t) \\ (2.39) \quad + & \sigma(t) \frac{C}{\nu^2} \left( \|F(t)\|_{L^3(\mathbb{R}^2)}^3 + \|G(t)\|_{L^3(\mathbb{R}^2)}^3 \right) + C \int_0^t \sigma \left[ \|\nabla u\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2} \|G\|_{L^4(\mathbb{R}^2)}^4 \right]. \end{aligned} $$
On the other hand, we have (2.14): $\frac{1}{\nu^2} \|G(t)\|_{L^3(\mathbb{R}^2)}^3 \le \frac{C^*}{\nu^2} E_0$ and hence the following thanks to Gagliardo-Nirenberg inequality and (2.16)-(2.18)

$$ (2.41) \qquad \begin{aligned} \frac{\sigma(t)}{\nu^2} \|F(t)\|_{L^3(\mathbb{R}^2)}^3 &\le C \frac{\sigma(t)}{\nu^2} \|\nabla F(t)\|_{L^2(\mathbb{R}^2)} \|F(t)\|_{L^2(\mathbb{R}^2)}^2 \\ &\le \sigma(t) \|\rho \dot{u}\|_{L^2(\mathbb{R}^2)} \left( \frac{1}{\nu^2} \|G\|_{L^2(\mathbb{R}^2)}^2 + \|\operatorname{div} u\|_{L^2(\mathbb{R}^2)}^2 \right) \\ &\le \eta \sigma(t) \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C^*}{\eta \nu^2} \left( \frac{1}{\nu^2} E_0^2 + (A_1(t))^2 \right). \end{aligned} $$

We finally combine (2.39), (2.40) and (2.41), we choose $\eta$ small to derive (2.4). $\square$

<figure><img src="image_3.png" /></figure>

### 2.2.3. Proof of Lemma 2.2.
This section is devoted to obtaining bounds for functionals $A_1$ and $A_2$ as defined in (1.15), for $d=3$. The proof is similar as in Lemma 2.1, and we adapt the estimates in three dimension, for instance, the $L^4(\mathbb{R}^2)$-norm is replaced by $L^6(\mathbb{R}^3)$-norm below. Since Proposition 2.8 does not hold in three dimension anymore, we will simply use the Sobolev embedding $\dot{H}^1(\mathbb{R}^3) \hookrightarrow L^6(\mathbb{R}^3)$ in the estimates.

**Proof of (2.5).** We recall that the first functional appears while using $\dot{u}$ as a test function in the weak formulation of (1.1)$_2$. By doing so, we obtain again (2.20):

$$ (2.42) \qquad \begin{aligned} A_1(t) ={}& \frac{\mu}{2} \|\nabla u_0\|_{L^2(\mathbb{R}^3)}^2 + \frac{\mu + \lambda}{2} \|\operatorname{div} u_0\|_{L^2(\mathbb{R}^3)}^2 - \mu \int_0^t \int_{\mathbb{R}^3} \nabla u^j \cdot \nabla u^k \partial_k u^j + \frac{\mu}{2} \int_0^t \int_{\mathbb{R}^3} |\nabla u|^2 \operatorname{div} u \\ &+ \frac{\mu + \lambda}{2} \int_0^t \int_{\mathbb{R}^3} (\operatorname{div} u)^3 - (\mu + \lambda) \int_0^t \int_{\mathbb{R}^3} \operatorname{div} u \nabla u^l \partial_l u + \int_0^t \int_{\mathbb{R}^3} \nabla u^l \partial_l u G \\ &+ \int_{\mathbb{R}^3} \operatorname{div} u(s) G(s) \Big|_{s=0}^{s=t} + \int_0^t \int_{\mathbb{R}^3} (\rho P'(\rho) - P(\rho) + \tilde{P})(\operatorname{div} u)^2. \end{aligned} $$

**Step 1: Estimates in terms of $E_0^\nu$ and $L_t^2 L_x^6$-norms of $(\nabla u, G, F)$**. With the help of Hölder's inequality the last two terms of the first line above can be straightforwardly estimated as follows:

$$ \left| - \mu \int_0^t \int_{\mathbb{R}^3} \nabla u^j \cdot \nabla u^k \partial_k u^j + \frac{\mu}{2} \int_0^t \int_{\mathbb{R}^3} |\nabla u|^2 \operatorname{div} u \right| \le C \|\nabla u\|_{L^3((0,t) \times \mathbb{R}^3)}^3. $$

Similarly, together with the relation between $\operatorname{div} u$ and $F, G$ as well as Hölder's inequality, the whole second line in (2.42) can be bounded by

$$ C \int_0^t \left( \| \operatorname{div} u \|_{L^2(\mathbb{R}^3)} \| \operatorname{div} u \|_{L^3(\mathbb{R}^3)} \|(F, G)\|_{L^6(\mathbb{R}^3)} + \| \nabla u \|_{L^2(\mathbb{R}^3)} \| \nabla u \|_{L^3(\mathbb{R}^3)} \|(F, \frac{1}{\nu} G)\|_{L^6(\mathbb{R}^3)} \right) $$

which is, by virtue of the interpolation inequality $\|f\|_{L^3(\mathbb{R}^3)} \lesssim \|f\|_{L^2(\mathbb{R}^3)}^{1/2} \|f\|_{L^6(\mathbb{R}^3)}^{1/2}$, bounded by

$$ C \int_0^t \left( \|\sqrt{\nu} \operatorname{div} u\|_{L^2(\mathbb{R}^3)}^{\frac{3}{2}} \frac{1}{\nu^{\frac{5}{4}}} \|(F, G)\|_{L^6(\mathbb{R}^3)}^{\frac{3}{2}} + \|\nabla u\|_{L^2(\mathbb{R}^3)}^{\frac{3}{2}} \|\nabla u\|_{L^6(\mathbb{R}^3)}^{\frac{1}{2}} \|(F, \frac{1}{\nu} G)\|_{L^6(\mathbb{R}^3)} \right). $$

This can be further estimated by Young's inequality by the following with some small constant $\eta > 0$:

$$ \eta \int_0^t \left\| \left( \nabla u, F, \frac{1}{\nu^{\frac{5}{6}}} G \right) \right\|_{L^6(\mathbb{R}^3)}^2 + \frac{C}{\eta} \int_0^t \left\| (\sqrt{\nu} \operatorname{div} u, \nabla u) \right\|_{L^2(\mathbb{R}^3)}^6. $$

By the same argument the whole third line in (2.42) can be bounded by

$$ \eta \nu \| \operatorname{div} u(t) \|_{L^2(\mathbb{R}^3)}^2 + \frac{C}{\eta} \| G(t) \|_{L^2(\mathbb{R}^3)}^2 + \| \operatorname{div} u_0 \|_{L^2(\mathbb{R}^3)} \| G_0 \|_{L^2(\mathbb{R}^3)} + C^* \| \operatorname{div} u \|_{L^2((0,t) \times \mathbb{R}^3)}^2 \\ \le \eta \nu \| \operatorname{div} u(t) \|_{L^2(\mathbb{R}^3)}^2 + \left( \frac{C^*}{\eta} + \frac{C^*}{\sqrt{\nu}} \right) E_0^\nu. $$
Recall Step 3 in the Proof of (2.3) above, and we get similar $L^4$-Estimates with the time weight $\sigma$. We have the following similar as (2.29)

$$ \frac{1}{2\nu} \int_{0}^{t} \sigma \|G\|_{L^4(\mathbb{R}^2)}^4 \le C^* E_0 + \frac{C}{\nu} \int_{0}^{t} \sigma \|F\|_{L^4(\mathbb{R}^2)}^4, $$

which implies

$$ (2.40) \qquad \int_{0}^{t} \sigma \left[ \|\nabla u\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2} \|F\|_{L^4(\mathbb{R}^2)}^4 + \frac{1}{\nu^2} \|G\|_{L^4(\mathbb{R}^2)}^4 \right] \le \frac{C^*}{\nu^3} E_0 + C^*(E_0 + \mathcal{A}_1(t))\mathcal{A}_1(t). $$

On the other hand, we have (2.14): $\frac{1}{\nu^2} \|G(t)\|_{L^3(\mathbb{R}^2)}^3 \le \frac{C^*}{\nu^2} E_0$ and hence the following thanks to Gagliardo-Nirenberg inequality and (2.16)-(2.18)

$$ (2.41) \qquad \begin{aligned} \frac{\sigma(t)}{\nu^2} \|F(t)\|_{L^3(\mathbb{R}^2)}^3 &\le C \frac{\sigma(t)}{\nu^2} \|\nabla F(t)\|_{L^2(\mathbb{R}^2)} \|F(t)\|_{L^2(\mathbb{R}^2)}^2 \\ &\le \sigma(t) \|\rho \dot{u}\|_{L^2(\mathbb{R}^2)} \left( \frac{1}{\nu^2} \|G\|_{L^2(\mathbb{R}^2)}^2 + \|\mathrm{div} \, u\|_{L^2(\mathbb{R}^2)}^2 \right) \\ &\le \eta \sigma(t) \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \frac{C^*}{\eta \nu^2} \left( \frac{1}{\nu^2} E_0^2 + (\mathcal{A}_1(t))^2 \right). \end{aligned} $$

We finally combine (2.39), (2.40) and (2.41), we choose $\eta$ small to derive (2.4). $\square$

### 2.2.3. Proof of Lemma 2.2.
This section is devoted to obtaining bounds for functionals $\mathcal{A}_1$ and $\mathcal{A}_2$ as defined in (1.16), for $d=3$. The proof is similar as in Lemma 2.1, and we adapt the estimates in three dimension, for instance, the $L^4(\mathbb{R}^2)$-norm is replaced by $L^6(\mathbb{R}^3)$-norm below. Since Proposition 2.8 does not hold in three dimension anymore, we will simply use the Sobolev embedding $\dot{H}^1(\mathbb{R}^3) \hookrightarrow L^6(\mathbb{R}^3)$ in the estimates.

**Proof of (2.5).** We recall that the first functional appears while using $\dot{u}$ as a test function in the weak formulation of (1.1)$_2$. By doing so, we obtain again (2.20):

$$ (2.42) \qquad \begin{aligned} \mathcal{A}_1(t) ={}& \frac{\mu}{2} \|\nabla u_0\|_{L^2(\mathbb{R}^3)}^2 + \frac{\mu + \lambda}{2} \|\mathrm{div} \, u_0\|_{L^2(\mathbb{R}^3)}^2 - \mu \int_0^t \int_{\mathbb{R}^3} \nabla u^j \cdot \nabla u^k \partial_k u^j + \frac{\mu}{2} \int_0^t \int_{\mathbb{R}^3} |\nabla u|^2 \mathrm{div} \, u \\ &+ \frac{\mu + \lambda}{2} \int_0^t \int_{\mathbb{R}^3} (\mathrm{div} \, u)^3 - (\mu + \lambda) \int_0^t \int_{\mathbb{R}^3} \mathrm{div} \, u \nabla u^l \partial_l u + \int_0^t \int_{\mathbb{R}^3} \nabla u^l \partial_l u G \\ &+ \int_{\mathbb{R}^3} \mathrm{div} \, u(s) G(s) \Big|_{s=0}^{s=t} + \int_0^t \int_{\mathbb{R}^3} (\rho P'(\rho) - P(\rho) + \tilde{P})(\mathrm{div} \, u)^2. \end{aligned} $$

**Step 1: Estimates in terms of $E_0^\nu$ and $L_t^2 L_x^6$-norms of $(\nabla u, G, F)$**. With the help of Hölder's inequality the last two terms of the first line above can be straightforwardly estimated as follows:

$$ \left| -\mu \int_0^t \int_{\mathbb{R}^3} \nabla u^j \cdot \nabla u^k \partial_k u^j + \frac{\mu}{2} \int_0^t \int_{\mathbb{R}^3} |\nabla u|^2 \mathrm{div} \, u \right| \le C \|\nabla u\|_{L^3((0,t) \times \mathbb{R}^3)}^3. $$

Similarly, together with the relation between $\mathrm{div} \, u$ and $F, G$ as well as Hölder's inequality, the whole second line in (2.42) can be bounded by

$$ C \int_{0}^{t} \left( \| \mathrm{div} \, u \|_{L^2(\mathbb{R}^3)} \| \mathrm{div} \, u \|_{L^3(\mathbb{R}^3)} \|(F, G)\|_{L^6(\mathbb{R}^3)} + \| \nabla u \|_{L^2(\mathbb{R}^3)} \| \nabla u \|_{L^3(\mathbb{R}^3)} \|(F, \frac{1}{\nu}G)\|_{L^6(\mathbb{R}^3)} \right) $$

which is, by virtue of the interpolation inequality $\|f\|_{L^3(\mathbb{R}^3)} \lesssim \|f\|_{L^2(\mathbb{R}^3)}^{1/2} \|f\|_{L^6(\mathbb{R}^3)}^{1/2}$, bounded by

$$ C \int_{0}^{t} \left( \|\sqrt{\nu} \mathrm{div} \, u\|_{L^2(\mathbb{R}^3)}^{\frac{3}{2}} \frac{1}{\nu^{\frac{5}{4}}} \|(F, G)\|_{L^6(\mathbb{R}^3)}^{\frac{3}{2}} + \|\nabla u\|_{L^2(\mathbb{R}^3)}^{\frac{3}{2}} \|\nabla u\|_{L^6(\mathbb{R}^3)}^{\frac{1}{2}} \|(F, \frac{1}{\nu}G)\|_{L^6(\mathbb{R}^3)} \right). $$
Summing up, we obtain the following by further applying the interpolation inequality and then Young's inequality to $||\nabla u||_{L^3((0,t)\times\mathbb{R}^3)}^3$:

(2.43)

$$A_1(t) \le \frac{C^*}{\eta} E_0^\nu + \eta \left( \nu \| \operatorname{div} u(t) \|_{L^2(\mathbb{R}^3)}^2 + \int_0^t \left\| \left( \nabla u, F, \frac{1}{\nu^6} G \right) \right\|_{L^6(\mathbb{R}^3)}^2 \right) + \frac{C}{\eta} \int_0^t \left\| (\sqrt{\nu} \operatorname{div} u, \nabla u) \right\|_{L^2(\mathbb{R}^3)}^6,$$

with some small parameter $\eta \le 1$ to be determined below.

**Step 2: Final Estimates.** We now turn to the estimate of the $L^2((0,t); L^6(\mathbb{R}^3))$-norm $(\nabla u, G, F)$, similar as in Step 3 in the proof of (2.3).

• $L_t^2 L_x^6$-Estimate for $G$. Recall (1.28) with $l = 5$:

$$\frac{d}{dt} \|H_5(\rho)\|_{L^1(\mathbb{R}^3)} + \frac{1}{\nu} \|G\|_{L^6(\mathbb{R}^3)}^6 \le \frac{1}{\nu} \|G\|_{L^6(\mathbb{R}^3)}^5 \|F\|_{L^6(\mathbb{R}^3)}.$$

We define the time-dependent function $h(t) := \|H_5(\rho)\|_{L^1(\mathbb{R}^3)}^{\frac{1}{3}}$. Thanks to the equivalence between $\|H_5(\rho)\|_{L^1(\mathbb{R}^3)}$ and $\|G\|_{L^6(\mathbb{R}^3)}^6$, we have $\|G\|_{L^6(\mathbb{R}^3)}^2 \sim C^* h$.

$$3 \frac{d}{dt} h + \frac{1}{\nu C^*} h \le \frac{1}{\nu} h^{\frac{1}{2}} \|F\|_{L^6(\mathbb{R}^3)},$$

which together with Young's inequality yields:

$$3 \sup_{[0,t]} h + \frac{1}{2\nu C^*} \int_0^t h \le 3h(0) + 2 \frac{C^*}{\nu} \int_0^t \|F\|_{L^6(\mathbb{R}^3)}^2.$$

Finally we find the following estimate for $G$ in terms of $E_0$ and $F$:

$$(2.44) \qquad \frac{1}{\nu} \int_0^t \|G\|_{L^6(\mathbb{R}^3)}^2 \le C^* E_0^{\frac{1}{3}} + \frac{C^*}{\nu} \int_0^t \|F\|_{L^6(\mathbb{R}^3)}^2.$$

• $L_t^2 L_x^6$-Estimate for $F$. We use the Sobolev embedding $\|g\|_{L^6(\mathbb{R}^3)} \lesssim \|\nabla g\|_{L^2(\mathbb{R}^3)}$ to bound the $L_x^6$-norm of $F$ by (2.18):

$$\int_0^t \|F\|_{L^6(\mathbb{R}^3)}^2 \le C \int_0^t \|\nabla F\|_{L^2(\mathbb{R}^3)}^2 = C \int_0^t \|\rho \dot{u}\|_{L^2(\mathbb{R}^3)}^2.$$

• $L_t^2 L_x^6$-Estimate for $\nabla u$. Similarly, by use of (2.16)-(2.17)-(2.18) the following inequality holds true

$$\int_0^t \|\nabla u\|_{L^6(\mathbb{R}^3)}^2 \le C \left( \int_0^t \|\operatorname{div} u\|_{L^6(\mathbb{R}^3)}^2 + \int_0^t \|\operatorname{curl} u\|_{L^6(\mathbb{R}^3)}^2 \right) \\ \le \frac{C}{\nu^2} \int_0^t \|(F, G)\|_{L^6(\mathbb{R}^3)}^2 + C \int_0^t \|\rho \dot{u}\|_{L^2(\mathbb{R}^3)}^2.$$

Finally (2.5) follows from (2.43) when we choose $\eta$ small enough.

□

*Proof of (2.6).* Here we derive estimate for the functional $A_2$ as defined in (1.15) for $d = 3$. We recall that it appears while rewriting the equation on the form (2.30), next applying the operator $\partial_t \cdot + \operatorname{div}(\cdot u)$ in order to obtain (2.31) which we test with the material derivative of the velocity $\dot{u}$. By doing so, we obtain (2.34):

$$(2.45) \qquad \frac{1}{2} \sigma(t) \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^3)}^2 + \mu \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^3)}^2 + \frac{\mu + \lambda}{\nu^2} \int_0^t \sigma \|\dot{F}\|_{L^2(\mathbb{R}^2)}^2 = \frac{1}{2} \int_0^{\sigma(t)} \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \sum_{k=1}^4 I_k,$$

where $I_k$, $k = 1, 2, 3, 4$ are given as in (2.34).
This can be further estimated by Young's inequality by the following with some small constant $\eta > 0$:

$$ \eta \int_0^t \left\| \left( \nabla u, F, \frac{1}{\nu^{\frac{5}{6}}} G \right) \right\|_{L^6(\mathbb{R}^3)}^2 + \frac{C}{\eta} \int_0^t \left\| (\sqrt{\nu} \operatorname{div} u, \nabla u) \right\|_{L^2(\mathbb{R}^3)}^6 $$

By the same argument the whole third line in (2.42) can be bounded by

$$ \eta \nu \| \operatorname{div} u(t) \|_{L^2(\mathbb{R}^3)}^2 + \frac{C}{\eta} \| G(t) \|_{L^2(\mathbb{R}^3)}^2 + \| \operatorname{div} u_0 \|_{L^2(\mathbb{R}^3)} \| G_0 \|_{L^2(\mathbb{R}^3)} + C^* \| \operatorname{div} u \|_{L^2((0,t) \times \mathbb{R}^3)}^2 \\ \leq \eta \nu \| \operatorname{div} u(t) \|_{L^2(\mathbb{R}^3)}^2 + \left( \frac{C^*}{\eta} + \frac{C^*}{\sqrt{\nu}} \right) E_0^\nu $$

Summing up, we obtain the following by further applying the interpolation inequality and then Young's inequality to $\|\nabla u\|_{L^3((0,t) \times \mathbb{R}^3)}^3$:

$$ A_1(t) \leq \frac{C^*}{\eta} E_0^\nu + \eta \left( \nu \| \operatorname{div} u(t) \|_{L^2(\mathbb{R}^3)}^2 + \int_0^t \left\| \left( \nabla u, F, \frac{1}{\nu^{\frac{5}{6}}} G \right) \right\|_{L^6(\mathbb{R}^3)}^2 \right) + \frac{C}{\eta} \int_0^t \left\| (\sqrt{\nu} \operatorname{div} u, \nabla u) \right\|_{L^2(\mathbb{R}^3)}^6 $$

with some small parameter $\eta \le 1$ to be determined below.

**Step 2: Final Estimates.** We now turn to the estimate of the $L^2((0,t); L^6(\mathbb{R}^3))$-norm $(\nabla u, G, F)$, similar as in Step 3 in the proof of (2.3).

* $L_t^2 L_x^6$-Estimate for $G$. Recall (1.29) with $l = 5$:

$$ \frac{d}{dt} \|H_5(\rho)\|_{L^1(\mathbb{R}^3)} + \frac{1}{\nu} \|G\|_{L^6(\mathbb{R}^3)}^6 \le \frac{1}{\nu} \|G\|_{L^6(\mathbb{R}^3)}^5 \|F\|_{L^6(\mathbb{R}^3)} $$

We define the time-dependent function $h(t) := \|H_5(\rho)\|_{L^1(\mathbb{R}^3)}^{\frac{1}{3}}$. Thanks to the equivalence between $\|H_5(\rho)\|_{L^1(\mathbb{R}^3)}$ and $\|G\|_{L^6(\mathbb{R}^3)}^6$, we have $\frac{1}{C^*}h \le \|G\|_{L^6(\mathbb{R}^3)}^2 \le C^*h$.

$$ 3 \frac{d}{dt} h + \frac{1}{\nu C^*} h \le \frac{1}{\nu} h^{\frac{1}{2}} \|F\|_{L^6(\mathbb{R}^3)} $$

which together with Young's inequality yields:

$$ 3 \sup_{[0,t]} h + \frac{1}{2\nu C^*} \int_0^t h \le 3h(0) + 2 \frac{C^*}{\nu} \int_0^t \|F\|_{L^6(\mathbb{R}^3)}^2 $$

Finally we find the following estimate for $G$ in terms of $E_0$ and $F$:

$$ (2.44) \qquad \frac{1}{\nu} \int_0^t \|G\|_{L^6(\mathbb{R}^3)}^2 \le C^* E_0^{\frac{1}{3}} + \frac{C^*}{\nu} \int_0^t \|F\|_{L^6(\mathbb{R}^3)}^2 $$

* $L_t^2 L_x^6$-Estimate for $F$. We use the Sobolev embedding $\|g\|_{L^6(\mathbb{R}^3)} \lesssim \|\nabla g\|_{L^2(\mathbb{R}^3)}$ to bound the $L_x^6$-norm of $F$ by (2.18):

$$ \int_0^t \|F\|_{L^6(\mathbb{R}^3)}^2 \le C \int_0^t \|\nabla F\|_{L^2(\mathbb{R}^3)}^2 = C \int_0^t \|\rho \dot{u}\|_{L^2(\mathbb{R}^3)}^2 $$

* $L_t^2 L_x^6$-Estimate for $\nabla u$. Similarly, by use of (2.16)-(2.17)-(2.18) the following inequality holds true

$$ \int_0^t \|\nabla u\|_{L^6(\mathbb{R}^3)}^2 \le C \left( \int_0^t \|\operatorname{div} u\|_{L^6(\mathbb{R}^3)}^2 + \int_0^t \|\operatorname{curl} u\|_{L^6(\mathbb{R}^3)}^2 \right) \\ \le \frac{C}{\nu^2} \int_0^t \|(F, G)\|_{L^6(\mathbb{R}^3)}^2 + C \int_0^t \|\rho \dot{u}\|_{L^2(\mathbb{R}^3)}^2 $$

Finally (2.5) follows from (2.43) when we choose $\eta$ small enough.

□
**Step 1 Estimates for $I_1$.** By the identity (2.35) and Hölder's inequality we achieve the estimates

$$ \left| \int_0^t \sigma \int_{\mathbb{R}^3} \dot{F} \nabla u^k \cdot \partial_k u \right| \lesssim \sigma(t) \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^2(\mathbb{R}^3)} \|\nabla u\|_{L^3(\mathbb{R}^3)} + \int_0^{\sigma(t)} \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^2(\mathbb{R}^3)} \|\nabla u\|_{L^3(\mathbb{R}^3)} \\ + \int_0^t \sigma \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^3(\mathbb{R}^3)} \|\nabla \dot{u}\|_{L^2(\mathbb{R}^3)} + \int_0^t \sigma \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^6(\mathbb{R}^3)}^2 \|\nabla u\|_{L^2(\mathbb{R}^3)}. $$

Similar as in the proof of (2.5) above, we use the interpolation $\|\nabla u\|_{L^3(\mathbb{R}^3)} \lesssim \|\nabla u\|_{L^2(\mathbb{R}^3)}^{1/2} \|\nabla u\|_{L^6(\mathbb{R}^3)}^{1/2}$ and Young's inequality to derive (noticing $\sigma(t) \le 1$)

$$ \left| \int_0^t \sigma \int_{\mathbb{R}^3} \dot{F} \nabla u^k \cdot \partial_k u \right| \\ \le C(\sqrt{\sigma(t)} \|F\|_{L^6(\mathbb{R}^3)}) \|\nabla u\|_{L^2(\mathbb{R}^3)}^{\frac{3}{2}} (\sqrt{\sigma(t)} \|\nabla u\|_{L^6(\mathbb{R}^3)})^{\frac{1}{2}} + \int_0^{\sigma(t)} \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^2(\mathbb{R}^3)}^{\frac{1}{2}} \|\nabla u\|_{L^6(\mathbb{R}^3)}^{\frac{3}{2}} \\ + \int_0^t \|\sqrt{\sigma} F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^6(\mathbb{R}^3)}^{\frac{1}{2}} \|\nabla u\|_{L^2(\mathbb{R}^3)}^{\frac{1}{2}} \|\sqrt{\sigma} \nabla \dot{u}\|_{L^2(\mathbb{R}^3)} + \int_0^t \sigma \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^6(\mathbb{R}^3)}^2 \|\nabla u\|_{L^2(\mathbb{R}^3)} \\ \le \eta (\sigma(t) \|(\nabla u, F)\|_{L^6(\mathbb{R}^3)}^2 + \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^3)}^2) \\ + \frac{C}{\eta} \left( \|\nabla u\|_{L^2(\mathbb{R}^3)}^6 + \int_0^t \sigma (\|F\|_{L^6(\mathbb{R}^3)}^2 \|\nabla u\|_{L^6(\mathbb{R}^3)} + \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^2(\mathbb{R}^3)}^2) \|\nabla u\|_{L^2(\mathbb{R}^3)} \right) \\ + \left( \int_0^{\sigma(t)} \|F\|_{L^6(\mathbb{R}^3)}^2 \right)^{\frac{1}{2}} \left( \int_0^{\sigma(t)} (\|\nabla u\|_{L^6(\mathbb{R}^3)}^2 + \|\nabla u\|_{L^2(\mathbb{R}^3)}^6) \right)^{\frac{1}{2}}. $$

By the Sobolev embedding $\dot{H}^1(\mathbb{R}^3) \hookrightarrow L^6(\mathbb{R}^3)$ and (2.16)-(2.17)-(2.18)-(2.44) we have

$$ \| (F, \mathrm{curl} \, u) \|_{L^6(\mathbb{R}^3)} \le C \rho^* \| \sqrt{\rho} \dot{u} \|_{L^2(\mathbb{R}^3)}, \\ \frac{1}{\nu} \int_0^t \| G \|_{L^6(\mathbb{R}^3)}^2 \le C^* E_0^{\frac{1}{3}} + \frac{C^*}{\nu} \int_0^t \| \sqrt{\rho} \dot{u} \|_{L^2(\mathbb{R}^3)}^2 \le C^* E_0^{\frac{1}{3}} + \frac{C^*}{\nu} A_1(t), \\ \| \nabla u \|_{L^6(\mathbb{R}^3)} \le C (\| \mathrm{div} \, u \|_{L^6(\mathbb{R}^3)} + \| \mathrm{curl} \, u \|_{L^6(\mathbb{R}^3)}) \le C \rho^* \| \sqrt{\rho} \dot{u} \|_{L^2(\mathbb{R}^3)} + \frac{C}{\nu} \| G \|_{L^6(\mathbb{R}^3)}, $$

which implies

$$ \sigma(t) \|(\nabla u, F)\|_{L^6(\mathbb{R}^3)}^2 \le C^* A_2(t) + \frac{C}{\nu^2} E_0^{\frac{1}{3}}, \\ \int_0^t \sigma \|F\|_{L^6(\mathbb{R}^3)} \|(F, \nabla u)\|_{L^6(\mathbb{R}^3)}^2 \|\nabla u\|_{L^2(\mathbb{R}^3)} \le C^* A_1(t)^{\frac{1}{2}} A_2(t)^{\frac{1}{2}} (A_1(t) + \frac{C^*}{\nu} E_0^{\frac{1}{3}}), \\ \left(\int_0^{\sigma(t)} \|F\|_{L^6(\mathbb{R}^3)}^2\right)^{\frac{1}{2}} \left(\int_0^{\sigma(t)} (\|\nabla u\|_{L^6(\mathbb{R}^3)}^2 + \|\nabla u\|_{L^2(\mathbb{R}^3)}^6)\right)^{\frac{1}{2}} \le C^* A_1(\sigma(t))^{\frac{1}{2}} (A_1(\sigma(t)) + \frac{1}{\nu} E_0^{\frac{1}{3}} + E_0 A_1(\sigma(t))^2)^{\frac{1}{2}}. $$

It holds by Young's inequality that

$$ \frac{\mu + \lambda}{\nu^2} \left| \int_0^t \sigma \int_{\mathbb{R}^3} \dot{F} \rho P'(\rho) \mathrm{div} \, u \right| \le \eta \frac{\mu + \lambda}{2\nu^2} \int_0^t \sigma \|\dot{F}\|_{L^2(\mathbb{R}^3)}^2 + \frac{C^*}{\nu^2 \eta} E_0. $$

Gathering all of these computations and using Young's inequality, we have:

$$ (2.46) \quad |I_1| \le C^* \left( \eta A_2(t) + \frac{1}{\nu} E_0^{\frac{1}{3}} + \frac{1}{\eta^2} A_1(t)^3 + \frac{1}{\eta} A_1(t) \frac{1}{\nu^2} E_0^{\frac{2}{3}} + A_1(\sigma(t)) + E_0 A_1(\sigma(t))^{\frac{3}{2}} \right). $$
Proof of (2.6). Here we derive estimate for the functional $A_2$ as defined in (1.16) for $d = 3$. We recall that it appears while rewriting the equation on the form (2.30), next applying the operator $\partial_t + \text{div}(\cdot u)$ in order to obtain (2.31) which we test with the material derivative of the velocity $\dot{u}$. By doing so, we obtain (2.34):

$$ (2.45) \quad \frac{1}{2} \sigma(t) \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^3)}^2 + \mu \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^3)}^2 + \frac{\mu + \lambda}{\nu^2} \int_0^t \sigma \|\dot{F}\|_{L^2(\mathbb{R}^2)}^2 = \frac{1}{2} \int_0^{\sigma(t)} \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^2 + \sum_{k=1}^4 I_k, $$

where $I_k$, $k = 1, 2, 3, 4$ are given as in (2.34).

**Step 1 Estimates for $I_1$.** By the identity (2.35) and Hölder's inequality we achieve the estimates

$$ \left| \int_0^t \sigma \int_{\mathbb{R}^3} \dot{F} \nabla u^k \cdot \partial_k u \right| \lesssim \sigma(t) \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^2(\mathbb{R}^3)} \|\nabla u\|_{L^3(\mathbb{R}^3)} + \int_0^{\sigma(t)} \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^2(\mathbb{R}^3)} \|\nabla u\|_{L^3(\mathbb{R}^3)} \\ + \int_0^t \sigma \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^3(\mathbb{R}^3)} \|\nabla \dot{u}\|_{L^2(\mathbb{R}^3)} + \int_0^t \sigma \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^6(\mathbb{R}^3)}^2 \|\nabla u\|_{L^2(\mathbb{R}^3)}. $$

Similar as in the proof of (2.5) above, we use the interpolation $\|\nabla u\|_{L^3(\mathbb{R}^3)} \lesssim \|\nabla u\|_{L^2(\mathbb{R}^3)}^{1/2} \|\nabla u\|_{L^6(\mathbb{R}^3)}^{1/2}$ and Young's inequality to derive (noticing $\sigma(t) \le 1$)

$$ \begin{aligned} \left| \int_0^t \sigma \int_{\mathbb{R}^3} \dot{F} \nabla u^k \cdot \partial_k u \right| & \le C(\sqrt{\sigma(t)} \|F\|_{L^6(\mathbb{R}^3)}) \|\nabla u\|_{L^2(\mathbb{R}^3)}^{\frac{3}{2}} (\sqrt{\sigma(t)} \|\nabla u\|_{L^6(\mathbb{R}^3)})^{\frac{1}{2}} + \int_0^{\sigma(t)} \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^6(\mathbb{R}^3)}^{\frac{1}{2}} \|\nabla u\|_{L^2(\mathbb{R}^3)}^{\frac{3}{2}} \\ & \quad + \int_0^t \|\sqrt{\sigma}F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^6(\mathbb{R}^3)}^{\frac{1}{2}} \|\nabla u\|_{L^2(\mathbb{R}^3)}^{\frac{1}{2}} \|\sqrt{\sigma}\nabla \dot{u}\|_{L^2(\mathbb{R}^3)} + \int_0^t \sigma \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^6(\mathbb{R}^3)}^2 \|\nabla u\|_{L^2(\mathbb{R}^3)} \\ & \le \eta (\sigma(t) \|(\nabla u, F)\|_{L^6(\mathbb{R}^3)}^2 + \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^3)}^2) \\ & \quad + \frac{C}{\eta} \left( \|\nabla u\|_{L^2(\mathbb{R}^3)}^6 + \int_0^t \sigma (\|F\|_{L^6(\mathbb{R}^3)}^2 \|\nabla u\|_{L^6(\mathbb{R}^3)} + \|F\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^6(\mathbb{R}^3)}^2) \|\nabla u\|_{L^2(\mathbb{R}^3)} \right) \\ & \quad + \left( \int_0^{\sigma(t)} \|F\|_{L^6(\mathbb{R}^3)}^2 \right)^{\frac{1}{2}} \left( \int_0^{\sigma(t)} (\|\nabla u\|_{L^6(\mathbb{R}^3)}^2 + \|\nabla u\|_{L^2(\mathbb{R}^3)}^6) \right)^{\frac{1}{2}}. \end{aligned} $$

By the Sobolev embedding $\dot{H}^1(\mathbb{R}^3) \hookrightarrow L^6(\mathbb{R}^3)$ and (2.16)-(2.17)-(2.18)-(2.44) we have

$$ \begin{aligned} \| (F, \operatorname{curl} u) \|_{L^6(\mathbb{R}^3)} &\le C \rho^* \| \sqrt{\rho} \dot{u} \|_{L^2(\mathbb{R}^3)}, \\ \frac{1}{\nu} \int_0^t \|G\|_{L^6(\mathbb{R}^3)}^2 &\le C^* E_0^{\frac{1}{3}} + \frac{C^*}{\nu} \int_0^t \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^3)}^2 \le C^* E_0^{\frac{1}{3}} + \frac{C^*}{\nu} A_1(t), \\ \|\nabla u\|_{L^6(\mathbb{R}^3)} &\le C (\|\operatorname{div} u\|_{L^6(\mathbb{R}^3)} + \|\operatorname{curl} u\|_{L^6(\mathbb{R}^3)}) \le C \rho^* \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^3)} + \frac{C}{\nu} \|G\|_{L^6(\mathbb{R}^3)}, \end{aligned} $$

which implies

$$ \begin{aligned} \sigma(t) \|(\nabla u, F)\|_{L^6(\mathbb{R}^3)}^2 &\le C^* A_2(t) + \frac{C}{\nu^2} E_0^{\frac{1}{3}}, \\ \int_0^t \sigma \|F\|_{L^6(\mathbb{R}^3)} \|(F, \nabla u)\|_{L^6(\mathbb{R}^3)}^2 \|\nabla u\|_{L^2(\mathbb{R}^3)} &\le C^* A_1(t)^{\frac{1}{2}} A_2(t)^{\frac{1}{2}} (A_1(t) + \frac{C^*}{\nu} E_0^{\frac{1}{3}}), \\ \left(\int_0^{\sigma(t)} \|F\|_{L^6(\mathbb{R}^3)}^2\right)^{\frac{1}{2}} \left(\int_0^{\sigma(t)} (\|\nabla u\|_{L^6(\mathbb{R}^3)}^2 + \|\nabla u\|_{L^2(\mathbb{R}^3)}^6)\right)^{\frac{1}{2}} &\le C^* A_1(\sigma(t))^{\frac{1}{2}} (A_1(\sigma(t)) + \frac{1}{\nu} E_0^{\frac{1}{3}} + E_0 A_1(\sigma(t))^2)^{\frac{1}{2}}. \end{aligned} $$

It holds by Young's inequality that

$$ \frac{\mu + \lambda}{\nu^2} \left| \int_0^t \sigma \int_{\mathbb{R}^3} \dot{F} \rho P'(\rho) \operatorname{div} u \right| \le \eta \frac{\mu + \lambda}{2\nu^2} \int_0^t \sigma \|\dot{F}\|_{L^2(\mathbb{R}^3)}^2 + \frac{C^*}{\nu^2 \eta} E_0. $$
**Step 2. Final estimates.** Recall the definitions of $I_k$, $k=2,3,4$ given in (2.34), and we can estimate them similarly as in Step 1 as follows:

$$
\begin{aligned}
|I_2 + I_3 + I_4| & \le C \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^3)} \|(\nabla u, F)\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^3(\mathbb{R}^3)} + C^* \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^3)} \|\nabla u\|_{L^2(\mathbb{R}^3)} \\
& \le \eta \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^3)}^2 + \frac{C^*}{\eta} \int_0^t \left( \sigma \|(\nabla u, F)\|_{L^6(\mathbb{R}^3)}^3 \|\nabla u\|_{L^2(\mathbb{R}^3)} + \|\nabla u\|_{L^2(\mathbb{R}^3)}^2 \right) \\
& \le \eta A_2(t) + \frac{C^*}{\eta} A_1(t)^{\frac{1}{2}} (A_2(t)^{\frac{1}{2}} A_1(t) + \frac{E_0^{\frac{2}{3}}}{\nu^2} + \frac{E_0^{\frac{1}{6}}}{\nu^3} A_1(t)) + \frac{C^*}{\eta} E_0.
\end{aligned}
$$

Gathering all of these computations, we obtain, after choosing $\eta$ small enough:

$$A_2(t) \le C^* \left( \frac{1}{\nu} E_0^{\frac{1}{3}} + E_0 + \frac{E_0^{\frac{2}{3}}}{\nu^2} A_1(t)^{\frac{1}{2}} + \left(1 + \frac{1}{\nu^2} E_0^{\frac{2}{3}}\right) A_1(t) + \left(E_0 + \frac{E_0^{\frac{1}{6}}}{\nu^3}\right) A_1(t)^{\frac{3}{2}} + A_1(t)^3 \right)$$

A further application of Young's inequality yields (2.6). □

**2.2.4.** *Proof of Lemma 2.3.* In this paragraph, we derive a priori estimate for the upper bound of the density in terms of $A_1, A_2$, under the assumption $\sup_{t,x} \rho(t,x) \le \rho^*$. The basic facts in Section 2.2.1 will be used freely.

*Proof.* With the help of the expression of $\mathrm{div} u$ in (2.16), we begin by rewriting the mass equation (1.1)$_1$ in terms of $F, G$ as follows:

$$\partial_t \rho + u \cdot \nabla \rho + \frac{\rho}{\nu} G = - \frac{\rho}{\nu} F.$$

Due to the fact that the pressure is an increasing function of the density and $G = P(\rho) - \tilde{P}$ the above equation yields:

$$(2.47) \qquad \partial_t |\rho - \tilde{\rho}| + u \cdot \nabla |\rho - \tilde{\rho}| + \frac{\rho}{\nu} |G| = - \frac{\rho}{\nu} \operatorname{sgn}(\rho - \tilde{\rho}) F.$$

It yields immediately the following $L^\infty$ estimate for the density, which we use on the short time interval $[0, \sigma(t)]$

$$(2.48) \qquad \sup_{[0, \sigma(t)]} \|\rho - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)} \le \|\rho_0 - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)} + \frac{\rho^*}{\nu} \int_0^{\sigma(t)} \|F(s)\|_{L^\infty(\mathbb{R}^d)} ds.$$

For larger time we would like to improve the $L_t^1$-norm for $\|F\|_{L_x^\infty}$ into $L_t^3$-norm, which requires less decay rate in time. From (2.47) we have

$$(2.49) \qquad \frac{1}{3} \left( \partial_t |\rho - \tilde{\rho}|^3 + u \cdot \nabla |\rho - \tilde{\rho}|^3 \right) + \frac{\rho}{\nu} |G| |\rho - \tilde{\rho}|^2 = - \frac{\rho}{\nu} \operatorname{sgn}(\rho - \tilde{\rho}) F |\rho - \tilde{\rho}|^2.$$

As well, since the pressure is an increasing function of the density such that

$$|G| |\rho - \tilde{\rho}|^2 \sim C^* |\rho - \tilde{\rho}|^3,$$

we derive the following estimate from Young's inequality on larger time interval $[\sigma(t), t]$

$$(2.50) \qquad \sup_{[\sigma(t), t]} \|\rho - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)}^3 \le \|\rho(\sigma(t)) - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)}^3 + \frac{C^*}{\nu} \int_{\sigma(t)}^t \|F(s)\|_{L^\infty(\mathbb{R}^d)}^3 ds.$$

Gathering estimates (2.48) and (2.50), we have:

$$(2.51) \qquad \sup_{[0, t]} \|\rho - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)} \le \|\rho_0 - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)} + \frac{C^*}{\nu} \int_0^{\sigma(t)} \|F(s)\|_{L^\infty(\mathbb{R}^d)} ds + \frac{C^*}{\nu^{1/3}} \left[ \int_{\sigma(t)}^t \|F(s)\|_{L^\infty(\mathbb{R}^d)}^3 ds \right]^{1/3}.$$

It only remains to estimate the norm of the effective flux $F$ in terms of the functionals $A_1$ and $A_2$ and to do so, we distinguish two cases according to the dimension.
Gathering all of these computations and using Young's inequality, we have:

$$ (2.46) \qquad |I_1| \le C^* \left( \eta A_2(t) + \frac{1}{\nu} E_0^{\frac{1}{3}} + \frac{1}{\eta^2} A_1(t)^3 + \frac{1}{\eta} A_1(t) \frac{1}{\nu^2} E_0^{\frac{2}{3}} + A_1(\sigma(t)) + E_0 A_1(\sigma(t))^{\frac{3}{2}} \right). $$

**Step 2. Final estimates.** Recall the definitions of $I_k$, $k = 2, 3, 4$ given in (2.34), and we can estimate them similarly as in Step 1 as follows:

$$ \begin{aligned} |I_2 + I_3 + I_4| &\le C \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^3)} \|(\nabla u, F)\|_{L^6(\mathbb{R}^3)} \|\nabla u\|_{L^3(\mathbb{R}^3)} + C^* \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^3)} \|\nabla u\|_{L^2(\mathbb{R}^3)} \\ &\le \eta \int_0^t \sigma \|\nabla \dot{u}\|_{L^2(\mathbb{R}^3)}^2 + \frac{C^*}{\eta} \int_0^t \left( \sigma \|(\nabla u, F)\|_{L^6(\mathbb{R}^3)}^3 \|\nabla u\|_{L^2(\mathbb{R}^3)} + \|\nabla u\|_{L^2(\mathbb{R}^3)}^2 \right) \\ &\le \eta A_2(t) + \frac{C^*}{\eta} A_1(t)^{\frac{1}{2}} (A_2(t)^{\frac{1}{2}} A_1(t) + \frac{E_0^{\frac{2}{3}}}{\nu^2} + \frac{E_0^{\frac{1}{6}}}{\nu^3} A_1(t)) + \frac{C^*}{\eta} E_0. \end{aligned} $$

Gathering all of these computations, we obtain, after choosing $\eta$ small enough:

$$ A_2(t) \le C^* \left( \frac{1}{\nu} E_0^{\frac{1}{3}} + E_0 + \frac{E_0^{\frac{2}{3}}}{\nu^2} A_1(t)^{\frac{1}{2}} + (1 + \frac{1}{\nu^2} E_0^{\frac{2}{3}}) A_1(t) + (E_0 + \frac{E_0^{\frac{1}{6}}}{\nu^3}) A_1(t)^{\frac{3}{2}} + A_1(t)^3 \right) $$

A further application of Young's inequality yields (2.6). □

**2.2.4.** *Proof of Lemma 2.3.* In this paragraph, we derive a priori estimate for the upper bound of the density in terms of $A_1, A_2$, under the assumption $\sup_{t,x} \rho(t,x) \le \rho^*$. The basic facts in Section 2.2.1 will be used freely.

*Proof.* With the help of the expression of $\text{div } u$ in (2.16), we begin by rewriting the mass equation (1.1)$_1$ in terms of $F, G$ as follows:

$$ \partial_t \rho + u \cdot \nabla \rho + \frac{\rho}{\nu} G = - \frac{\rho}{\nu} F. $$

Due to the fact that the pressure is an increasing function of the density and $G = P(\rho) - \tilde{P}$ the above equation yields:

$$ (2.47) \qquad \partial_t |\rho - \tilde{\rho}| + u \cdot \nabla |\rho - \tilde{\rho}| + \frac{\rho}{\nu} |G| = - \frac{\rho}{\nu} \operatorname{sgn}(\rho - \tilde{\rho}) F. $$

It yields immediately the following $L^\infty$ estimate for the density, which we use on the short time interval $[0, \sigma(t)]$

$$ (2.48) \qquad \sup_{[0, \sigma(t)]} \|\rho - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)} \le \|\rho_0 - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)} + \frac{\rho^*}{\nu} \int_0^{\sigma(t)} \|F(s)\|_{L^\infty(\mathbb{R}^d)} ds. $$

For larger time we would like to improve the $L_t^1$-norm for $\|F\|_{L_x^\infty}$ into $L_t^3$-norm, which requires less decay rate in time. From (2.47) we have

$$ (2.49) \qquad \frac{1}{3} \left( \partial_t |\rho - \tilde{\rho}|^3 + u \cdot \nabla |\rho - \tilde{\rho}|^3 \right) + \frac{\rho}{\nu} |G| |\rho - \tilde{\rho}|^2 = - \frac{\rho}{\nu} \operatorname{sgn}(\rho - \tilde{\rho}) F |\rho - \tilde{\rho}|^2. $$

As well, since the pressure is an increasing function of the density such that

$$ |G| |\rho - \tilde{\rho}|^2 \sim C^* |\rho - \tilde{\rho}|^3, $$

we derive the following estimate from Young's inequality on larger time interval $[\sigma(t), t]$

$$ (2.50) \qquad \sup_{[\sigma(t), t]} \|\rho - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)}^3 \le \|\rho(\sigma(t)) - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)}^3 + \frac{C^*}{\nu} \int_{\sigma(t)}^t \|F(s)\|_{L^\infty(\mathbb{R}^d)}^3 ds. $$

Gathering estimates (2.48) and (2.50), we have:

$$ (2.51) \qquad \sup_{[0,t]} \|\rho - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)} \le \|\rho_0 - \tilde{\rho}\|_{L^\infty(\mathbb{R}^d)} + \frac{C^*}{\nu} \int_0^{\sigma(t)} \|F(s)\|_{L^\infty(\mathbb{R}^d)} ds + \frac{C^*}{\nu^{1/3}} \left[ \int_{\sigma(t)}^t \|F(s)\|_{L^\infty(\mathbb{R}^d)}^3 ds \right]^{1/3}. $$
**Case d = 2.** We recall that the effective flux is given by

$$ (2.52) \qquad F = -(-\Delta)^{-1} \operatorname{div}(\rho \dot{u}) = \nu \operatorname{div} u - G, \quad G = P(\rho) - \tilde{P}. $$

Hence, interpolation inequality yields the following estimate on the small interval $[0, \sigma(t)]$

$$ \begin{aligned} \int_0^{\sigma(t)} \|F\|_{L^\infty(\mathbb{R}^2)} \le C \int_0^{\sigma(t)} \|F\|_{L^2(\mathbb{R}^2)}^{\frac{1}{2}} \|\nabla F\|_{L^4(\mathbb{R}^2)}^{\frac{2}{3}} &\le \int_0^{\sigma(t)} \|(G, \nu \operatorname{div} u)\|_{L^2(\mathbb{R}^2)}^{\frac{1}{3}} \|\rho \dot{u}\|_{L^4(\mathbb{R}^2)}^{\frac{2}{3}} \\ &\le (C^* E_0^{\frac{1}{6}} + \nu^{\frac{1}{6}} A_1(t)^{\frac{1}{6}}) \int_0^{\sigma(t)} \|\rho \dot{u}\|_{L^4(\mathbb{R}^2)}^{\frac{2}{3}}. \end{aligned} $$

Since the density contains vacuum states, we are not allowed to bound the last factor in the above inequality solely by the Gagliardo-Nirenberg inequality. On $\mathbb{T}^2$, we have at our disposal a logarithmic interpolation inequality (see [16, 17, 21]) which is not valid in the whole space. To address this issue, we prove in **Lemma A.1** an interpolation inequality that will allow us to take into account the vacuum state. Thus, from **Lemma A.1**, we have:

$$ \begin{aligned} \int_0^{\sigma(t)} \|\rho \dot{u}\|_{L^4(\mathbb{R}^2)}^{\frac{2}{3}} &\le C^* \int_0^{\sigma(t)} \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^{\frac{1}{3}} \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^{\frac{1}{3}} + C^* E_0^{\frac{1}{18}} \int_0^{\sigma(t)} \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^{\frac{2}{9}} \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^{\frac{4}{9}} \\ &\le C^* A_1(t)^{\frac{1}{6}} A_2(t)^{\frac{1}{6}} + C^* E_0^{\frac{1}{18}} A_1(t)^{\frac{1}{9}} A_2(t)^{\frac{2}{9}} \end{aligned} $$

and therefore:

$$ (2.53) \qquad \int_0^{\sigma(t)} \|F\|_{L^\infty(\mathbb{R}^2)} \le C^*(1 + E_0^{\frac{1}{18}})(E_0^{\frac{1}{6}} + \nu^{\frac{1}{6}} A_1(t)^{\frac{1}{6}})(A_1(t)^{\frac{1}{3}} + A_2(t)^{\frac{1}{3}}). $$

Similarly, interpolation inequality and **Lemma A.1** yield on the larger time interval $[\sigma(t), t]$

$$ \begin{aligned} \int_{\sigma(t)}^{t} \|F\|_{L^{\infty}(\mathbb{R}^2)}^3 &\le \int_{\sigma(t)}^{t} \|(G, \nu \operatorname{div} u)\|_{L^2(\mathbb{R}^2)} \|\rho \dot{u}\|_{L^4(\mathbb{R}^2)}^2 \\ &\le (C^* E_0^{\frac{1}{2}} + \nu^{\frac{1}{2}} A_1(t)^{\frac{1}{2}}) (C^* \int_0^t \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)} \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)} + C^* E_0^{\frac{1}{6}} \int_0^t \|\sqrt{\rho} \dot{u}\|_{L^2(\mathbb{R}^2)}^{\frac{2}{3}} \|\nabla \dot{u}\|_{L^2(\mathbb{R}^2)}^{\frac{4}{3}}) \\ &\le (C^* E_0^{\frac{1}{2}} + \nu^{\frac{1}{2}} A_1(t)^{\frac{1}{2}}) (C^* A_1(t)^{\frac{1}{2}} A_2(t)^{\frac{1}{2}} + C^* E_0^{\frac{1}{6}} A_1(t)^{\frac{1}{3}} A_2(t)^{\frac{2}{3}}). \end{aligned} $$

Hence,

$$ (2.54) \qquad \int_{\sigma(t)}^{t} \|F\|_{L^{\infty}(\mathbb{R}^2)}^3 \le C^*(1 + E_0^{\frac{1}{6}})(C^* E_0^{\frac{1}{2}} + \nu^{\frac{1}{2}} A_1(t)^{\frac{1}{2}})(A_1(t) + A_2(t)). $$

Finally (2.51), (2.53) and (2.54), lead to:

$$ (2.55) \qquad \sup_{[0,t]} \|\rho - \tilde{\rho}\|_{L^\infty(\mathbb{R}^2)} \le \|\rho_0 - \tilde{\rho}\|_{L^\infty(\mathbb{R}^2)} + \frac{C^*}{\nu^{\frac{1}{3}}} (1 + E_0^{\frac{1}{18}})(E_0^{\frac{1}{6}} + \nu^{\frac{1}{6}} A_1(t)^{\frac{1}{6}})(A_1(t)^{\frac{1}{3}} + A_2(t)^{\frac{1}{3}}). $$

**Case d = 3.** From the expression of the effective flux (2.52), we have by Gagliardo-Nirenberg inequality:

$$ \begin{aligned} \int_0^{\sigma(t)} \|F\|_{L^\infty(\mathbb{R}^3)} &\le C \int_0^{\sigma(t)} \|F\|_{L^6(\mathbb{R}^3)}^{\frac{1}{2}} \|\nabla F\|_{L^6(\mathbb{R}^3)}^{\frac{1}{2}} \le C \int_0^{\sigma(t)} \|\rho \dot{u}\|_{L^2(\mathbb{R}^3)}^{\frac{1}{2}} \|\rho \dot{u}\|_{L^6(\mathbb{R}^3)}^{\frac{1}{2}} \le C^* A_1(t)^{\frac{1}{4}} A_2(t)^{\frac{1}{4}}, \\ \int_{\sigma(t)}^{t} \|F\|_{L^\infty(\mathbb{R}^3)}^3 &\le C \int_{\sigma(t)}^{t} \|\rho \dot{u}\|_{L^2(\mathbb{R}^3)}^{\frac{3}{2}} \|\rho \dot{u}\|_{L^6(\mathbb{R}^3)}^{\frac{3}{2}} \le C^* A_1(t)^{\frac{1}{4}} A_2(t)^{\frac{5}{4}}. \end{aligned} $$

Finally, we get:

$$ \sup_{[0,t]} \|\rho - \tilde{\rho}\|_{L^\infty(\mathbb{R}^3)} \le \|\rho_0 - \tilde{\rho}\|_{L^\infty(\mathbb{R}^3)} + \frac{C^*}{\nu^{1/3}} (A_1(t)^{\frac{1}{2}} + A_2(t)^{\frac{1}{2}}), $$

and this ends the proof of **Lemma 2.3**.

□