\section{Introduction} \label{sec:intro}

\bigbreak

We are interested in the following {\em inhomogeneous incompressible} Navier-Stokes system:
\begin{equation} \label{sys}
\left\{ \begin{array}{lr}
\varrho_t+v\cdot \nabla\varrho=0 & {\rm in } \ \:\mathbb{R}_+\times\Omega,  \\[1ex]
\varrho v_t+\varrho v \cdot\nabla v-\Delta v +\nabla P=0 & {\rm in }\  \:\mathbb{R}_+\times \Omega,\\[1ex]
\div v=0 & {\rm in }\ \: \mathbb{R}_+\times \Omega.
\end{array}
\right.
\end{equation}
This system describes the motion of incompressible fluids with variable density and originates from simplified models in geophysics.
The unknowns in the above system are the velocity $v$, the density $\varrho$ and the pressure $P,$ depending on the time variable $t\geq0$
and on the space variable $x\in\Omega$ where the fluid domain $\Omega$
 is  a smooth bounded  subset  of $\mathbb{R}^d$ with  the physical  dimensions $d=2,3.$

The system is supplemented with the initial data
\begin{equation} \label{ic}
    \varrho|_{t=0}=\varrho_0 \quad\hbox{and}\quad  v|_{t=0} =v_0.
\end{equation}
At the boundary, we prescribe the no-slip condition
\begin{equation} \label{bc}
v|_{\partial\Omega} = 0.
\end{equation}

The existence of weak solutions to \eqref{sys} is nowadays well understood and the state of the art on this issue is rather similar to that of the classical incompressible Navier-Stokes system (i.e. with constant density). The analysis goes back to the work of A. Kazhikhov \cite{Kaz} who showed global existence of weak solutions for initial density bounded away from zero. This constraint was removed by J. Simon in \cite{Simon}. Later, P.-L. Lions \cite{PLL} showed that the density is a renormalized solution to the continuity equation, this allowed in particular to treat the case of density-dependent viscosity in \cite{Des}.
Still in the framework of weak solutions,
F. Fanelli and I. Gallagher investigated recently in \cite{FG} the fast rotation limit of \eqref{sys} supplemented with
a Coriolis force.

Producing  `strong  solutions' (by strong, we mean solutions
having  the uniqueness property) requires more constraints on the data: enough regularity and no vacuum, typically.
Roughly speaking, according to the classical literature,
 for smooth enough data and provided
the density does not vanish,
we have  global existence of strong solutions
\emph{even for  large data}  in dimension two,
and, like for the constant density
case, for small enough initial velocity in dimension three.
For such results in the bounded domain case, one can
refer to the pioneering work by O.A. Ladyzhenskaya and
V.~Solonnikov in \cite{LS} (further extended to less regular data
by the  first author in \cite{D2}).

A number of works have been dedicated to solving \eqref{sys} in  $\Omega=\mathbb{R}^d$ in so-called
`critical regularity frameworks'. The underlying idea
(that originates from Fujita and Kato's work \cite{FK}
for the constant density case) is that `optimal'
functional spaces for well-posedness of \eqref{sys}
have to share its scaling invariance, namely, for all
$\ell>0,$
\begin{equation}\label{eq:critical}
(\varrho,v,P)(t,x) \leadsto (\varrho,\ell v,\ell^2P)(\ell^2t,\ell x)\quad\hbox{and}\quad (\varrho_0,v_0)(x) \leadsto (\varrho_0,\ell v_0)(\ell x).
\end{equation}
Observing that the couple of homogeneous Besov space
$\dot B^{\frac d2}_{2,1}(\mathbb{R}^d)\times
\bigl(\dot B^{\frac d2-1}_{2,1}(\mathbb{R}^d)\bigr)^d$
indeed possesses this invariance, the first author proved in
 \cite{D0} the well-posedness of \eqref{sys} supplemented
 with initial velocity $v_0$  in $\dot B^{\frac d2-1}_{2,1}(\mathbb{R}^d)$
and initial density  $\varrho_0$  close
to some positive constant in $\dot B^{\frac d2}_{2,1}(\mathbb{R}^d).$
Note that, owing to the embedding $\dot B^{\frac d2}_{2,1}(\mathbb{R}^d)\hookrightarrow {\mathcal C}_b(\mathbb{R}^d),$
this forces the density to be continuous.
Subsequent  improvements have been brought to this approach
(see e.g. \cite{AP}) but,
still, the density  has to be `almost' continuous.
In particular, one cannot consider initial densities
that have  a jump across an interface, even a smooth one.
\medbreak
Toward considering less regular densities, a first
breakthrough has been brought by the first two authors
in \cite{DM0} and \cite{DM5}:
 taking advantage of Lagrangian coordinates (that will be
presented below), they
established well-posedness results
for densities that are possibly discontinuous
along interfaces, provided the jump is small enough.

Then, in \cite{PZZ}, by a totally different approach, M. Paicu,
P. Zhang
and Z. Zhang succeeded in proving the global existence in $\mathbb{R}^2$ for $v_0 \in H^s, \; s>0$ and in $\mathbb{R}^3$ for $v_0 \in H^1$ with $\|v_0\|_2\|\nabla v_0\|_2$ sufficiently small,
provided the initial density satisfies
$$0<c_0\leq\varrho_0\leq C_0<\infty.$$
In dimension $3,$ this work was extended in \cite{CZZ} to initial
velocities that are only in $H^s$ for some $s>1/2.$
Still, the density  has to be  bounded  away from zero, and the solution in not time continuous with values in $H^s.$
Very recently, P. Zhang in \cite{Zhang20} achieved the critical regularity  $\dot B^{1/2}_{2,1}$ for the initial velocity,
but did not address the uniqueness issue in this setting.

For more results where the initial density is allowed
to be discontinuous but, still, strictly positive,
the reader may  refer among others to \cite{GGJ,HPZ,FZ} and to a recent result \cite{CJ}, where an inflow boundary condition is considered. Let us also mention that global well-posedness in the half-space $\mathbb{R}^d_+$ with initial density only bounded but close to a positive constant was shown in \cite{DZ1}.
\medbreak
All the above results require the strict positivity of the initial density.
To our knowledge, the existence of unique solutions in presence of
vacuum has been  first proved in \cite{CK} for rather high
regularity of the initial density and velocity,
(namely $\varrho_0 \in L^{3/2} \cap H^2$ and $u_0\in H^2$)
and provided the following compatibility condition is satisfied:
\begin{equation} \label{comp}
-\Delta v_0+\nabla P_0=\sqrt{\varrho_0}g \quad \textrm{for some} \; g \in L_2 \; {\rm and} \; P_0 \in H^1.
\end{equation}
Global existence of unique solutions in a 3D bounded domain or in  $\mathbb{R}^3$ under the same compatibility condition and smallness of $\|u_0\|_{\dot H^{1/2}}$ was shown in \cite{CHW}.

Condition \eqref{comp} was removed in \cite{Li}, where local well-posedness in a bounded domain is shown, but still for sufficiently smooth initial density. Global existence in the whole space $\mathbb{R}^3$, again under sufficient regularity of initial density, was proved recently in \cite{HLL}.
  \medbreak
An important place in the theory of \eqref{sys} is taken by
the so-called {\em density patch problem}: assuming that
\begin{equation} \label{ic0}
    \varrho_0=\alpha_1 \chi_{A_0} + \alpha_2 \chi_{\Omega \setminus A_0}
    \end{equation}
 for some nonnegative constants $\alpha_1, \alpha_2$ and a measurable set $A_0,$ can we say that $\varrho(t)$ has the same structure
 for all time,
 with persistence of the regularity of the interface~?
 This problem seems to have been first raised
 by P.-L. Lions in \cite{PLL} in the specific case
 where   $\varrho_0=\chi_{A_0}$ with $A_0 \in \mathbb{R}^2,$ and $\sqrt{\varrho_0}u_0 \in L_2.$ The original question was whether for all time $\varrho(t)=\chi_{A(t)}$
for some domain  $A(t)$ with the same regularity  as $A_0.$

It turns out that a positive answer is obtained for the $C^1$ regularity as a consequence of the works of the first two authors in
\cite{DM0,DM5}  if  $\alpha_1, \alpha_2>0$  are close to each other.
Much more complete results have been obtained in the two-dimensional case in \cite{LZ1} (case $\alpha_1-\alpha_2$ small), and then in \cite{LZ2} for any positive constants. There, the
authors actually establish the persistence of high `striated' Sobolev regularity for the density.
Similar results have been proved in the 3D case in \cite{LL}.
The propagation of striated regularity has been adapted to the case
where the viscosity depends on the density in \cite{PZ2}.
By a different approach, persistence of H\"older continuity
of the interface if $\alpha_1,\alpha_2$ are close to each other was shown in \cite{DZ2}.
\smallbreak
Requiring that the initial density is away from zero
precludes to consider the original Lions' problem, namely
the case when  $\alpha_2=0$ in \eqref{ic0}.
Recently in \cite{DM1}, the first and second authors proved the
well-posedness of \eqref{sys} for only bounded initial density
\begin{equation} \label{rho0}
0 \leq \varrho_0 \leq \varrho^*
\end{equation}
and initial velocity satisfying
\begin{equation} \label{v0}
v_0 \in H^1_0(\Omega), \quad \div v_0 =0.
\end{equation}
In the two-dimensional case, the solutions are global
 without any additional condition while,  in the 3D case,
  $v_0$  has to satisfy some smallness condition
  (as the results of \cite{DM1} are of particular importance for our analysis, they will be recalled  precisely below).
As a by-product, the authors obtained a positive answer to the Lions question in the case $\varrho_0=\chi_{A_0}$:
persistence of H\"older regularity $C^{1,\alpha}$
holds true for any $0<\alpha<1$ in 2D and $0<\alpha<\frac{1}{2}$ in 3D.
\smallbreak
However, the question concerning the stability of the solutions
was left open in \cite{DM1}. In fact, if the density is bounded away from zero
then  the stability can be proved in the same way as uniqueness, but this is no longer the case if the initial density  is allowed to vanish
(this has to do with  the parabolic character of the momentum equation, which degenerates if the density vanishes). This typically happens if we consider
the following model configuration:
 the original density is  $\varrho^{(org)}_0$
and the perturbation  is  $\varrho^{(per)}_0,$ that is
\begin{equation}\label{ic2}
    \varrho_0^{(org)}=\chi_{A(0)}\quad\hbox{and}\quad\varrho_0^{(per)}=\chi_{A(0)} + \phi.
\end{equation}

\begin{figure}[h!] \label{model}
\centering
\includegraphics[width=0.55\textwidth]{rys1.pdf}
\caption{Support of $\varrho_0^{(org)} = {\bf A}$ and support of $\varrho_0^{(per)}={\bf B}$.}
\end{figure}
The problem happens whenever
\begin{equation*}
    ({\rm supp}\; \varrho_0^{(org)} \setminus {\rm supp}\; \varrho_0^{(per)} )\cup
    ({\rm supp}\; \varrho_0^{(per)} \setminus {\rm supp}\; \varrho_0^{(org)})
    \ \hbox{ is not empty}.\end{equation*}
Indeed, unless $\phi=0,$ the perturbation is large in $L_\infty(\Omega)$
but small in $L_p(\Omega)$ for $p<\infty,$ and we need a functional framework for stability that captures this situation.
\medbreak
The goal of the present paper
is to show the stability of solutions to the inhomogeneous Navier-Stokes \eqref{sys} with respect to initial condition of type \eqref{ic}, in a regularity setting that includes the density patches problem \eqref{ic0}
\emph{even if one of the parameters $\alpha_1,\,\alpha_2$ vanishes}.
In particular, we allow the density to be a characteristic function of some set. In the most pathological case the supports of $\varrho_0^{(org)}$
and $\varrho_0^{(per)}$ can be even disjoint.

\medbreak
Before stating our main stability result, we have to recall
the state-of-the-art concerning the global well-posedness
for  \eqref{sys} supplemented with general data satisfying \eqref{rho0} and \eqref{v0}.

In dimension $d=2,$   \cite[Thm 2.1]{DM1} states the following result:
\begin{thm} \label{thm:d=2}
Let $\Omega$ be a smooth bounded domain of $\mathbb{R}^2,$ or a two-dimensional torus. Let $\varrho_0 \in L_\infty(\Omega)$ satisfy \eqref{rho0} and let $v_0$ satisfy \eqref{v0}. Then, System \eqref{sys} admits a unique solution $(\varrho,v,P)$ such that
\begin{equation}\label{reg:1}
\begin{aligned}
&\varrho\in L_\infty(\mathbb{R}_+;L_\infty),\quad  v \in L_\infty(\mathbb{R}_+;H^1), \\ &\sqrt{\varrho}v_t, \nabla^2 v,\nabla P \in L_2(\mathbb{R}_+;L_2), \quad
\nabla v\in L_{1,loc}(\mathbb{R}_+;L_\infty),\\
&\sqrt\varrho v\in {\mathcal C}(\mathbb{R}_+;L_2)\quad\hbox{and}\quad \varrho\in{\mathcal C}(\mathbb{R}_+;L_p) \quad \hbox{for all}\ 1\leq p <\infty.
\end{aligned}
\end{equation}
For arbitrarily large but finite time $T>0$ these solutions satisfy
in addition for all $1 \leq r <2, \; 1 \leq m<\infty, s<{1}/{2}$ and $1 \leq p <\infty,$
\begin{equation} \label{reg:2d}
\begin{aligned}
&\nabla(\sqrt{t}P), \nabla^2(\sqrt{t}v) \in L_{\infty}(0,T;L_r(\Omega)) \cap L_2(0,T;L_m(\Omega)), \; v \in H^s(0,T;L_p(\Omega)),\\
& \sqrt{t \varrho} \,v_t \in L_\infty(0,T;L_2) \mbox{ \ and \ } \nabla v_t \in L_2(0,T;L_2).
\end{aligned}
\end{equation}
\end{thm}
\noindent
In the three dimensional case, we know from  \cite[Thm 2.2]{DM1} the following result:
\begin{thm} \label{thm:d=3}
Let $\Omega$ be a smooth bounded domain of $\mathbb{R}^3,$ or a three-dimensional torus. Let $\varrho_0 \in L_\infty(\Omega)$ satisfy \eqref{rho0} and $v_0$ satisfy \eqref{v0}.
There exists $c>0$ such that if, in addition,
\begin{equation} \label{small1}
(\varrho^*)^{3/2}\|\sqrt{\varrho_0}v_0\|_2\|\nabla v_0\|_2 \leq c\mu^2,
\end{equation}
then System \eqref{sys} admits a unique solution $(\varrho,v,P)$ satisfying \eqref{reg:1} and, for any finite $T>0$,
$s<{1}/{2}$ and $1 \leq p <\infty,$
\begin{equation} \label{reg:3d}
\begin{aligned}
&\nabla(\sqrt{t}P), \nabla^2(\sqrt{t}v) \in L_{\infty}(0,T;L_2(\Omega)) \cap L_2(0,T;L_6(\Omega)), \; v \in H^s(0,T;L_p(\Omega)), \\
& \sqrt{t \varrho} \,v_t \in L_\infty(0,T;L_2) \mbox{ \ and \ } \nabla v_t \in L_2(0,T;L_2).
\end{aligned}
\end{equation}
\end{thm}
Although the above solutions are  unique, the question of their stability remains open so far.
Here we aim at supplementing the above statements
with a stability result.
In order to obtain the most accurate information,
it is natural  to use Lagrangian coordinates
since, in this setting,  the density is time independent (it only depends on the position of the particles initially).
Therefore,  the problem is reduced to the control of the difference of the velocities
which, somehow, satisfies a parabolic equation.
\medbreak
Let us shortly recall how to define  Lagrangian coordinates
in our setting. First, we introduce the flow
$X:(t,y)\mapsto X(t,y)$ of $v,$
that is  the unique solution to the following ODE:
\begin{equation}\label{L-E}
\left\{    \begin{array}{lr}
   \displaystyle \frac{dX}{dt}=v(t,X(t,y))     & \mbox{in \ }  \mathbb{R}_+\times \Omega, \\[10pt]
   \displaystyle X(0,y)=y     &  \mbox{in \ } \Omega.
    \end{array}\right.
\end{equation}
Integrating \eqref{L-E} yields the following relation between the Eulerian $``x"$ and Lagrangian $``y"$ coordinates:
\begin{equation} \label{Lag1}
    x= X(t,y)=y+\int_0^t v(t',X(t',y))\,{\rm d}t'.
\end{equation}
By the standard theory of ODEs, the above change of coordinates
is well defined
whenever $v \in L_{1,(loc)}(\mathbb{R}_+;C^{1,0})$. In the coordinates system $(t,y),$  the unknown functions are named as follows:
\begin{equation*}
    u(t,y)=v(t,X(t,y)), \;\;\; \eta(t,y)=\varrho(t,X(t,y)), \;\;\; Q(t,y)=P(t,X(t,y)).
\end{equation*}
Let us denote
\begin{equation} \label{Au}
A_u(t)=\left(\frac{dX}{dy}\right)^{-1}=\left(\hbox{\rm Id}+\int_0^t \nabla_y u(t',y)\,{\rm d}t'\right)^{-1}=\left[{\rm cof}\left(\hbox{\rm Id}+\int_0^t \nabla_y u(t',y)\,{\rm d}t' \right)\right]^T,
\end{equation}
where $\rm cof(\cdot)$ denotes the cofactor matrix and, in the second equality, we used the fact that ${\rm det}A_u=1$ (see e.g. \cite{DM1}).
For a function $f(t,x),$ denote $\tilde f(t,y)=f(t,X(t,y))$.
Then, owing to the chain rule,
\begin{equation} \label{trans1}
\nabla_x f(t,x)=A_u^T\nabla_y \tilde f(t,y)=:\nabla_u \tilde f(t,y), \quad
\partial_t f(t,x)+v(t,x)\cdot\nabla_xf(t,x)= \partial_t \tilde f(t,y).
\end{equation}
In order to transform the divergence operator, observe that
$\, \hbox{\rm div}_y\,   A_u = 0$ since $\det A_u=1$ .
Therefore, for any vector field $z(t,y)$, one may write
\begin{equation} \label{trans:div1}
\, \hbox{\rm div}_y\,  (A_u z)=A^T : \nabla_y z + z \cdot \, \hbox{\rm div}_y\,  (A_u)=A^T : \nabla_y z.
\end{equation}
Hence, if we denote $\tilde w(t,y)=w(t,X(t,y))$
 for any vector field $w(t,x),$ we discover that
\begin{equation} \label{div:id}
\, \hbox{\rm div}_x\,   w(t,x)=A_u : \nabla_y \tilde w(t,y)=\, \hbox{\rm div}_y\,  (A_u \tilde w)=:\, \hbox{\rm div}_u\,   \tilde w.
\end{equation}
Taking all above into account we see that in coordinates $(t,y)$
system \eqref{sys} reads
\begin{equation} \label{sys:lag}
\left\{\begin{array}{llr}
    \eta_t=0     & \mbox{in}& \mathbb{R}_+\times\Omega,  \\[1ex]
    \eta u_t -\, \hbox{\rm div}_u\,  \nabla_u u + \nabla_u Q =0 &\mbox{in}& \mathbb{R}_+\times\Omega,\\   [1ex]
    \, \hbox{\rm div}_u\,   u=0 &\mbox{in}& \mathbb{R}_+\times\Omega,\\[1ex]
    u|_{t=0}=v_0, \;\; \eta|_{t=0}=\varrho_0 &\mbox{in}& \Omega.
    \end{array}\right.
\end{equation}
The main achievement of this paper  is the
following stability result in the Lagrangian coordinates setting.
\begin{thm} \label{thm:lag}
Let $\Omega$ be a smooth bounded domain of $\mathbb{R}^d$ with $d=2,3$.
Let $(\varrho^1,v^1)$ and $(\varrho^2,v^2)$ be two solutions of \eqref{sys}
with initial data $(\varrho_0^1,v_0^1)$ and $(\varrho_0^2,v_0^2)$, respectively,  with non identically zero bounded
 $\varrho_0^1$ and $\varrho_0^2,$ and $v_0^1,$ $v_0^2$ in $H^1_0(\Omega),$
 given either by Theorem \ref{thm:d=2} or by Theorem \ref{thm:d=3} (depending on the dimension).
Denote by $u^1$ and $u^2$ the
corresponding velocities in Lagrangian coordinates.
Finally, set $\delta\! u:=u^2-u^1,$ $\delta\! v:=v^2-v^1$ and
$\delta\! \vr:=\varrho^2-\varrho^1.$

Then, there exists a positive constant $\beta$ depending
only on the shape of $\Omega,$ and such that
\begin{equation} \label{decay:final}
\sup\limits_{t\in\mathbb{R}_+} \Bigl(e^{\frac{\beta\mu t} {\varrho^*\delta^2}}[\|{\rm \min}\Big\{\sqrt{\varrho^1_0},\sqrt{\varrho^2_0}\Big\}\delta\! u(t)\|_2\Bigr)
+ \|e^{\frac{\beta\mu t} {\varrho^*\delta^2}}\nabla \delta\! u\|_{L_2(\mathbb{R}_+,L_2)}\leq C \left( \|\sqrt{\varrho^1_0}\,\delta\! v_0\|_2 + \|\delta\! \vr_0\|_2^{1/2} \right),
\end{equation}
where $\delta$ stands for the diameter of $\Omega,$
and $\varrho^*:=\|\varrho_0\|_{\infty}.$
\end{thm}

Coming back from Lagrangian to Eulerian coordinates, we obtain

\begin{cor} \label{cor:eul}
Under the assumptions of Theorem \ref{thm:lag}, in Eulerian coordinates, we have
 $$  \sup_{t\in\mathbb{R}_+} \|\delta\! \vr(t)\|_{W^{-1}_p}
   + \|e^{\frac{\beta\mu t} {\varrho^*\delta^2}}\nabla \delta\! v\|_{L_2(\mathbb{R}_+,L_2)}
   \leq C_0 \left(
    \|\sqrt{\varrho^1_0}\,\delta\! v_0\|_2 + \|\delta\! \vr_0\|_2^{1/2}\right)$$
for $1<p<\infty$ if $d=2$ and $1<p\leq 6$ if $d=3$.
\end{cor}
In our low regularity setting, the key difficulty for proving
stability of solutions to \eqref{sys} comes from the partially
hyperbolic nature of the system. In fact,  if writing
the system satisfied by the difference
$(\delta\! \vr,\delta\! v,dP)$ of two solutions $(\varrho^1,v^1,P^1)$
and $(\varrho^2,v^2,P^2)$ of \eqref{sys}, then  the mass equation gives:
$$(\partial_t+v^2\cdot\nabla)\delta\! v=\delta\! v\cdot\nabla\varrho^2.$$
In our framework where the density is only in $L_\infty(\Omega),$
this forces us to perform estimates for $\delta\! \vr$ in
a space with regularity index equal to $-1.$
Following the duality approach initiated by D. Hoff in \cite{Hoff} and recently
renewed in \cite{DM2} and \cite{MP} (for
the related compressible Navier-Stokes system), we shall
actually prove stability estimates for
the density in  $W^{-1}_p,$  and in $L_2$ for the velocity.

A key point in order to get rid of any smallness condition, is to first establish a sufficiently strong time-decay of the solutions that have been constructed in Theorems \ref{thm:d=2} and \ref{thm:d=3}.  This will be achieved
in Section \ref{s:decay}, where by means of rather classical energy arguments, we will even get \emph{exponential} decay,
  owing to the boundedness of the fluid domain.
Before that, in Section \ref{s:stability}, we shall  compare two global solutions satisfying {\em a priori} those decay properties, and estimate their difference in Lagrangian coordinates
in terms of the difference of the data, getting the result of Theorem \ref{thm:lag}. Finally, we rewrite those estimates in the Eulerian coordinates to obtain Corollary \ref{cor:eul}.

Proving the exponential decay estimates of  Section \ref{s:decay}
 can be achieved by means of a remarkably simple  energy method
 that is performed directly on the nonlinear problem.
 This is in sharp contrast with the proof of decay estimates
 for compressible Navier-Stokes and related models
 which requires a refined analysis of the linearized
 system combined with a perturbation argument
 (see among others \cite{DHP,ShiShi,ES,PSZ} and the references therein).

When  comparing  Theorem \ref{thm:lag} with results of \cite{DM1}, a remark is in order concerning the domain. Although Theorems \ref{thm:d=2} and \ref{thm:d=3} hold both for a bounded domain with Dirichlet condition and a torus, we here restricted our analysis to the case of a bounded domain
to  avoid further technical complications. In fact, in  the case of Dirichlet boundary conditions,
we have  the  basic Poincar\'e inequality at our disposal, which is fundamental  to close the estimates
globally in time. In the torus case,  the corresponding
 Poincar\'e inequality has an additional term (namely the total momentum of the solution, see Lemma A.1 in \cite{DM1})
 which, although probably harmless, entails serious complications in the
 proof of  decay estimates. Therefore, we leave the torus case for future research.
\smallbreak

{\bf Notation.} We use standard notations $L_p$ and $W^k_p$ for Lebesque and Sobolev spaces (the dependency
with respect to the fluid domain $\Omega$ is omitted).
For the corresponding norms, we use the short notation:
$$
\|\cdot\|_p := \|\cdot\|_{L_p}, \qquad \|\cdot\|_{k,p} := \|\cdot\|_{W^k_p}.  $$
In some computations, we will agree
that  $f_p(t)$ denotes a generic function of time which is in $L_p(\mathbb{R}_+)$, and  that  $f_{p,q}(t)$ stands for a function which is in $L_p(\mathbb{R}_+) \cap L_q(\mathbb{R}_+)$. The precise form of these functions may vary from line to line, but the property of
integrability is preserved.

\section{Stability under given decay properties}\label{s:stability}

Throughout this section, we are given two solutions
$(\varrho^1,v^1,P^1)$ and $(\varrho^2,v^2,P^2)$ pertaining to data  $(\varrho_0^1,v_0^1)$ and $(\varrho_0^2,v_0^2),$ and satisfying
the following properties for some $\beta_0>0$:

\begin{subequations}\label{eq:regularity}
\begin{align}
&\sqrt{\varrho}\, e^{\beta_0 t} v_t^i \in L_2(\mathbb{R}_+,L_2), \label{eq:regularity:8}\\[3pt]
& e^{\beta_0 t} \nabla v^i\in L_1(\mathbb{R}_+;L_\infty)\cap L_4(\mathbb{R}_+;L_3)\cap L_2(\mathbb{R}_+;L_6),\label{eq:regularity:2}\\[3pt]
&\sqrt t\, e^{\beta_0 t} v^i\in L_1\cap L_\infty(\mathbb{R}_+;L_\infty),\label{eq:regularity:3}\\[3pt]
&\sqrt t\, e^{\beta_0 t}(\nabla^2v^i,\nabla  P^i)\in L_2(\mathbb{R}_+;L_6),\label{eq:regularity:4}\\[3pt]
&\sqrt t\,e^{\beta_0 t} v_t^i \in L_{4/3}(\mathbb{R}_+;L_6), \label{eq:regularity:5}\\[3pt]
& e^{\beta_0 t} v^i  \in L_1\cap L_\infty(\mathbb{R}_+;L_6), \label{eq:regularity:6}\\[3pt]
&\sqrt t\, e^{\beta_0 t} \nabla v^i\in L_2(\mathbb{R}_+;L_\infty),\label{eq:regularity:1}\\[3pt]
& e^{\beta_0 t} \nabla^2  v^i \in L_1(\mathbb{R}_+;L_r) \; \textrm{with} \; r>d. \label{eq:regularity:7}
\end{align}
\end{subequations}
\noindent
We denote by $(\eta^i,u^i,Q^i)$ the corresponding solutions in Lagrangian coordinates (hence $\eta^i=\varrho_0^i$).

\begin{lem} \label{l:est_lag}
Let $(\varrho,v,P)$ solve \eqref{sys} and satisfy Conditions \eqref{eq:regularity:8} to \eqref{eq:regularity:7}.
 Then, the corresponding Lagrangian solution  $(\eta,u,Q)$ also satisfies  \eqref{eq:regularity:8} to \eqref{eq:regularity:7}.
\end{lem}
\begin{proof}
All  properties except for the ones involving a time derivative
and the last one just follow from
the corresponding ones  for $v,$ and from the fact that  the matrix $A_u$ is bounded.

For proving  \eqref{eq:regularity:8}, we start from the identity
$$\sqrt\eta e^{\beta_0t} u_t=\sqrt{\rho\circ X_t}\,e^{\beta_0t}(v_t+v\cdot \nabla v)\circ X_t.
$$
As $X_t$ is measure preserving, the term  with $v_t$ may be bounded by means of
\eqref{eq:regularity:8}. For bounding the other term, it suffices to observe that
$$\|e^{\beta_0t}\sqrt{\rho\circ X_t}\,v\cdot \nabla v)\circ X_t\|_{L_2(\mathbb{R}_+\times\Omega)}
\leq \rho^*\|e^{\beta_0t/2} v\|_{L_\infty(\mathbb{R}_+;L_6)}\|e^{\beta_0t/2}\nabla v\|_{L_2(\mathbb{R}_+;L_3)}.$$
The first term of the right-hand side may be bounded thanks to \eqref{eq:regularity:6}
while the second one, according to \eqref{eq:regularity:2} and the boundedness of $\Omega.$
As regards  \eqref{eq:regularity:5},  one can use again
$u_t=(v_t+v\cdot \nabla v)\circ X_t,$ and  properties    \eqref{eq:regularity:5} to  \eqref{eq:regularity:1} for $v.$

In order to prove \eqref{eq:regularity:7}, we differentiate the identity $\nabla_y u(t,y)={}^T\!(A_u)^{-1}\nabla_x v(t,X(t,y))$ with respect to $y$. By \eqref{Au} we obtain $$
\|\nabla^2_y u\|_{L_1(0,T;L_r)}\leq C\|\nabla_x^2 v\|_{L_1(0,T;L_r)} +
C\|\nabla_y A\|_{L_\infty(0,T;L_r)}\|\nabla_y u\|_{L_1(0,T;L_\infty)},
$$
which implies \eqref{eq:regularity:7} for $u$ due to  \eqref{eq:regularity:7} for $v$  and to \eqref{eq:regularity:2} for $u$.
\end{proof}

In the rest of this section, we aim at estimating
$$
\delta\! u:= u^2-u^1 \mbox{ \ \  and  \ \ } \delta\! Q:= Q^2-Q^1
$$
in terms of the difference of the data.  Obviously, denoting $\Delta_u:=\, \hbox{\rm div}_u\,  \nabla_u$ and $
\delta\! v_0=v^2_0-v^1_0$, the couple $(\delta\! u,\delta\! Q)$ satisfies
\begin{equation}\label{NSL-uniq}
 \begin{aligned}
  &{\varrho_0^1}\delta\! u_t - \Delta_{u^1}\delta\! u + \nabla_{u^1}\delta\! Q=
(\Delta_{u^2} -\Delta_{u^1}) u^2 - (\nabla_{u^2}-\nabla_{u^1}) Q^2-\delta\! \vr_0 \,u^2_t,\\
&\div_{\!u^1}\delta\! u=(\div_{\!u^1}-\div_{\!u^2})u^2,\\
&\delta\! u|_{t=0}=\delta\! v_0.
\end{aligned}
\end{equation}
Note that
\begin{equation*}
    \delta\! \vr_0:=\varrho_0^2 - \varrho_0^1=\eta^2-\eta^1.
\end{equation*}
By $(\ref{sys:lag})_1,$ functions $\eta^i$ are constant in time, so the perturbation of the density stays the same in time. It is one of the main reasons why we chose the Lagrange coordinates approach to deal with the
stability issue of system \eqref{sys}.

\subsection{The case of a nice control of vacuum}

Compared to the proof of uniqueness that has been performed in \cite{DM1}, the troublemaker is
 the term $\delta\! \vr_0 \,u^2_t$    in the equation for $\delta\! u$ since  Theorems \ref{thm:d=2}
 or \ref{thm:d=3} only provide us with an information on $\sqrt{\varrho^2}\, v_t^2$ (hence on $\sqrt{\varrho_0^2}\,u_t^2$)
while we do not  necessarily have
 \begin{equation}\label{x3}
 {\rm supp}\; \delta\! \vr_0 \subset {\rm supp}\; \varrho_0^2.
\end{equation}
In this part, we  assume that the initial densities satisfy
\begin{equation} \label{def:xnorm}
\|\delta\! \vr_0\|_X \,:= \, \|{\delta\! \vr_0}/{\sqrt{\varrho_0^2}}\|_4 < \infty
\end{equation}
and  derive a differential inequality which is crucial for proving Theorem \ref{thm:lag}.
The general case, when \eqref{x3} is not valid,  will be investigated in the  next subsection.
\medbreak
The idea is to decompose $\delta\! u$ into
\begin{equation}\label{l5}
 \delta\! u = w+z,
\end{equation}
where $w$ stands for a suitable solution to the divergence equation
\begin{equation}\label{l6}
 \div_{\!u^1} w = (\div_{\!u^1}-\div_{\!u^2})u^2={}^T\!\delta\! A: \nabla u^2=\div (\delta\! A \, u^2),
\end{equation}
with $\delta\! A:=A^1-A^2$ and $A^i:=A_{u^i}.$
\smallbreak
Since at $t=0$ by \eqref{Au} we have $\delta\! A=0$, we put $w=0$ at $t=0$.
Note that $A^1$ and $A^2$ need not to be close to $\hbox{\rm Id}$, but  are invertible and uniformly bounded in time.

\begin{lem}\label{Lem:9}
There exists a solution to \eqref{l6} such that
\begin{equation} \label{est:w}
\begin{aligned}
& \|e^{\beta_0 t}w(t)\|_{2} \leq  f_{1,\infty}(t) \biggl(\int_0^t \|\nabla \delta\! u(\tau)\|_2^2 \,{\rm d}\tau\biggr)^{1/2}, \\
& \|e^{\beta_0 t}\nabla w(t)\|_{2} \leq f_2(t)  \biggl(\int_0^t \| \nabla \delta\! u(\tau)\|_2^2\, d\tau\biggr)^{1/2}, \\
& \|e^{\beta_0 t}w_t(t)\|_{3/2} \leq f_{4/3}(t) \biggl(\int_0^t \| \nabla \delta\! u(\tau)\|_2^2\, d\tau\biggr)^{1/2}+
f_{4}(t)\|\nabla \delta\! u(t)\|_{2},
\end{aligned}
\end{equation}
where the notation $f_p(t)$ and $f_{p,q}(t)$ is explained at the end of Section \ref{sec:intro}.
\end{lem}
\begin{proof}
The vector $w$ will be sought in the form
$    w=(A^1)^{-1} A^1 w= (A^1)^{-1} \bar w $,
where  $\bar w$ is given as a solution to
\begin{equation}\label{l7}
 \div \bar w \big(= \div (A^1 w) = \div_{u^1} w\big) = (\div_{\!u^1}-\div_{\!u^2})u^2= {}^T\!\delta\! A: \nabla u^2 =\div (\delta\! A \, u^2).
\end{equation}

\medbreak
As a first step in the proof of our claim, let us establish the following bounds:
\begin{equation}\label{eq:www}\begin{aligned}
\|e^{\beta_0 t}\bar w\|_{L_4(0,T;L_2)}&\leq C\|e^{\beta_0 t}\delta\! A \,u^2\|_{L_4(0,T;L_2)}, \\
\|e^{\beta_0 t}\nabla \bar w\|_{L_2(0,T;L_2)}&\leq C\|e^{\beta_0 t}\;{}^T\!\delta\! A:\nabla u^2\|_{L_2(0,T;L_2)}\\
\hbox{\rm and }\
\|e^{\beta_0 t}\bar w_t\|_{L_{4/3}(0,T;L_{3/2})}&\leq C\|e^{\beta_0 t}(\delta\! A \,u^2)_t\|_{L_{4/3}(0,T;L_{3/2})}.
\end{aligned}\end{equation}
The existence of a vector field $\bar w$ satisfying \eqref{l7}-\eqref{eq:www} is assured by the following Lemma:
\begin{lem} \label{lem:Bog}
Let $A$ be a matrix valued function with $\det A\equiv1.$
Consider the following divergence equation in a bounded domain with smooth boundary:
\begin{equation*}
    \div b =f \mbox{ \ \ in \ }  \mathbb{R}_+\times\Omega,  \qquad b=0 \mbox{ \ \ at \ }  \mathbb{R}_+\times\partial\Omega,
\end{equation*}
where $f=A^T:\nabla d=\div(Ad)$ and the average of $f$ is equal zero.

Then, there exists a  constant $C$  such that for any $\beta\geq0,$
there exists a solution $b$  to the above equation, such that for all
$t\geq0,$ we have
\begin{equation*}
\begin{aligned}
&\|e^{\beta t} b(t)\|_2 \leq C\|e^{\beta t} A(t) d(t)\|_2, \qquad \|e^{\beta t}\nabla b(t)\|_2 \leq C \|e^{\beta t}A(t)^T \nabla d(t)\|_2,\\[5pt]
&\|e^{\beta t}b_t(t)\|_{3/2} \leq C\|e^{\beta t}(Ad)_t(t)\|_{3/2}.
\end{aligned}
\end{equation*}
\end{lem}
Lemma \ref{lem:Bog} has been proved without exponential weight by the first two authors in \cite{DM3} (see also
\cite[Lemma A.3]{DM1}). In their proof the function $b$ is given by an explicit formula, so a time weight can be treated as a multiplicative parameter (as the norms in Lemma \ref{lem:Bog} only involve
the space variable).
\medbreak
Now, let us bound the right-hand sides of \eqref{eq:www}.
In order to emphasize that we do not need any smallness of $\int_0^t\nabla_yu\,\,{\rm d}\tau,$
let us derive an explicit formula for $\delta\! A.$

In the two dimensional case, starting from \eqref{Au},  we immediately obtain
\begin{equation} \label{dA:2d}
\delta\! A(t) = \left[\begin{array}{cc}
\int_0^t \delta\! u_{2,y_2}\,{\rm d}t' & -\int_0^t \delta\! u_{1,y_2}\,{\rm d}t'\\[1ex]
-\int_0^t \delta\! u_{2,y_1}\,{\rm d}t' & \int_0^t \delta\! u_{1,y_1}\,{\rm d}t'
\end{array}
\right]\cdotp
\end{equation}
In three dimensions we can also use \eqref{Au}, let us compute precisely one entry of $\delta\! A$. We have
$$
(A_u^i)_{11}(t)=\biggl(1+\int_0^t u^i_{2,y_2}\,{\rm d}t'\biggr)\biggl(1+\int_0^t u^i_{3,y_3}\,{\rm d}t'\biggr)-
\int_0^t u^i_{3,y_2}\,{\rm d}t'\:\int_0^t u^i_{2,y_3}\,{\rm d}t'
$$
for $i=1,2$, therefore
$$\displaylines{\quad
(\delta\! A(t))_{11}= \int_0^t \delta\! u_{2,y_2}\,{\rm d}t' + \int_0^t \delta\! u_{3,y_3}\,{\rm d}t'
+\int_0^t \delta\! u_{2,y_2}\,{\rm d}t'\int_0^t u^2_{3,y_3}\,{\rm d}t'+\int_0^t \delta\! u_{3,y_3}\,{\rm d}t'\int_0^t u^1_{2,y_2}\,{\rm d}t'
\hfill\cr\hfill
-\int_0^t \delta\! u_{3,y_2}\,{\rm d}t'\int_0^t u^2_{2,y_3}\,{\rm d}t'
-\int_0^t \delta\! u_{2,y_3}\,{\rm d}t'\int_0^t u^1_{3,y_2}\,{\rm d}t'.\quad}$$
Other entries will have a similar structure, we can write them as
\begin{equation} \label{dA:3d}
(\delta\! A(t))_{ij} = \sum_{1 \leq k,l \leq 3} a^{ij}_{kl}\int_0^t\delta\! u_{k,y_l}\,{\rm d}t'
+ \sum_{1 \leq k,l,m,n \leq 3, \; s\in\{1,2\} }b^{ij}_{k,l,m,n,s}\int_0^t \delta\! u_{k,y_l}\,{\rm d}t'\:
\int_0^t u^s_{m,y_n}\,{\rm d}t',
\end{equation}
where $a^{ij}_{kl},b^{ij}_{k,l,m,n,s} \in \{0,1\}$.
Now if $u^1,u^2$ satisfy \eqref{eq:regularity:2} then, by H\"older inequality, we obtain
\begin{equation}\label{l128}
  \|t^{-1/2} \delta\! A(t)\|_2 \leq C\Big\|t^{-1/2}\int_0^t \nabla \delta\! u\,{\rm d}\tau\Bigr\|_{2}
  \leq C\biggl(\int_0^t\|\nabla\delta\! u\|_2^2\,{\rm d}\tau\biggr)^{1/2}.
\end{equation}
Therefore, by \eqref{eq:regularity:1}, we have
\begin{equation}\label{eq:new}
\|e^{\beta_0 t}\;{}^T\!\delta\! A:\nabla u^2(t)\|_2\leq
   \|t^{-1/2} \delta\! A(t)\|_2\|e^{\beta_0 t}t^{1/2}\nabla u^2\|_{\infty} \leq  C
   f_2(t) \left(\int_0^t\|\nabla\delta\! u\|_2^2\,{\rm d}\tau\right)^{1/2},
 \end{equation}
 which implies \eqref{est:w} for $\|\nabla \bar w\|_{L_2}$.
Similarly, by \eqref{eq:regularity:3} and \eqref{l128}
\begin{align*}
\|e^{\beta_0 t}\delta\! A\, u^2\|_2
\leq \|t^{-1/2}\delta\! A\|_2\|e^{\beta_0 t}t^{1/2}u^2\|_{\infty}
\leq f_{1,\infty}(t) \biggl(\int_0^t\|\nabla\delta\! u(\tau)\|_2^2\,{\rm d}\tau\biggr)^{1/2},
\end{align*}
whence, applying \eqref{eq:www} gives
\begin{equation}\label{l9} \|e^{\beta_0 t} \bar w(t)\|_2
  \leq f_1(t) \biggl(\int_0^t\|\nabla\delta\! u(\tau)\|_2^2\,{\rm d}\tau\biggr)^{1/2}.
\end{equation}
In order to bound $\bar w_t,$ it suffices to derive an appropriate estimate
in $L_{4/3}(0,T;L_{3/2})$ for
$$(\delta\! A\,u^2)_t=\delta\! A\,u^2_t+(\delta\! A)_t\,u^2.$$
For the first term, thanks to \eqref{eq:regularity:5} and \eqref{l128} we have
$$
\begin{aligned}
\|e^{\beta_0 t}\delta\! A\,u^2_t\|_{3/2}&\leq \|t^{-1/2}\delta\! A\|_2\|e^{\beta_0 t}t^{1/2}u^2_t\|_6
\leq f_{4/3}(t)\biggl(\int_0^t\|\nabla\delta\! u(\tau)\|_2^2\,{\rm d}\tau\biggr)^{1/2}.
\end{aligned}
$$
The other term can be bounded as follows:
$$\|e^{\beta_0 t}(\delta\! A)_t\,u^2\|_{3/2}\leq
\|(\delta\! A)_t\|_2\|e^{\beta_0 t}u^2\|_6.$$
Differentiating \eqref{dA:3d} with respect to $t$ and using \eqref{eq:regularity:1} for $u^1$ and $u^2,$
we see that
$$
\|\delta\! A_t(t)\|_{2}\leq C\biggl(\|\nabla\delta\! u(t)\|_{2}+\Big\|t^{-1/2}\int_0^t \nabla \delta\! u(\tau)\,{\rm d}\tau\Big\|_2
\bigl(\|t^{1/2}\nabla u^1(t)\|_\infty+\|t^{1/2}\nabla u^2(t)\|_\infty\bigr)\biggr).
$$
In the 2D case, owing to \eqref{dA:2d}, one can skip the second term on the right-hand side
of the above inequality. Thus,
owing to \eqref{eq:regularity:6}, for both $d=2,3$ we obtain
$$
\|\delta\! A_t\,u^2(t)\|_{3/2}\leq f_{4/3}(t) \biggl(\int_0^t\|\nabla\delta\! u(\tau)\|_2^2\,{\rm d}\tau\biggr)^{1/2}+
f_4(t)\|\nabla \delta\! u(t)\|_{2},
$$
which, by \eqref{eq:www} implies \eqref{est:w} for $\|\bar w_t\|_{3/2}$.
Altogether, this gives the thesis of Lemma \ref{Lem:9}
for $\bar w$, but not yet for $w$.
\smallbreak

In order to get \eqref{est:w} for $w$ from the estimates of $\bar w,$ we
first observe that
\begin{equation} \label{A1}
\sup_{t\leq T}\|(A_u)^{-1}(t)\|_\infty \leq C \|\nabla u\|_{L_1(0,T;L_\infty)}.
\end{equation}
Therefore
$$    \|e^{\beta_0 t}w\|_{L_4(0,T;L_2)} \leq C\|A^1\|_{\infty}\|e^{\beta_0 t}\bar w \|_{L_4(0,T;L_2)} \leq C \|e^{\beta_0 t}\bar w \|_{L_4(0,T;L_2)}. $$
In order to estimate $\|\nabla w\|_{L^2(0,T;L^2)}$ we proceed as follows:
\begin{align*}
\|e^{\beta_0 t}\bar w \nabla((A^1)^{-1})\|_{L_2(0,T;L_2)} & \leq
\|\nabla((A^1)^{-1})\|_{L_\infty(0,T;L_r)}\|e^{\beta_0 t}\bar w\|_{L_2(0,T;L_{r*})}\\
& \leq C \|\nabla^2 u^1\|_{L_1(0,T;L_r)}\|e^{\beta_0 t}\bar w\|_{L_2(0,T;L_{r*})}
\end{align*}
with $\frac12 =\frac{1}{r} + \frac{1}{r^*}$, where in the last passage we used \eqref{Au}. Therefore by \eqref{eq:regularity:2} and \eqref{eq:regularity:7}
$$
\begin{aligned}
    \|e^{\beta_0 t}\nabla w\|_{L_2(0,T;L_2)} &\leq \|e^{\beta_0 t}(A^1)^{-1} \nabla\bar w\|_{L_2(0,T;L_2)} + \|e^{\beta_0 t}\bar w \nabla((A^1)^{-1})\|_{L_2(0,T;L_2)}\\
    & \leq C [\|e^{\beta_0 t}\nabla \bar w\|_{L_2(0,T;L_2)} + \|e^{\beta_0 t} \bar w\|_{L_2(0,T;L_{r^*})} ].
\end{aligned}$$
\noindent
Since, in  \eqref{eq:regularity:7}, one can take $r>d$,  one can always ensure that  $r^*<6.$
\smallbreak
Hence, Finally we have
$$
\|(A^1)^{-1})_t \bar w\|_{3/2} \leq \|(A^{-1})_t\|_6 \|\bar w\|_2
\leq \|\nabla u^1\|_6 \|\bar w\|_2,
$$
which together with \eqref{A1} implies
$$\begin{aligned}
\|e^{\beta_0 t}w_t\|_{L_{4/3}(0,T;L_{3/2})}
&\leq C\|(A^1)^{-1}\|_{\infty} \|e^{\beta_0 t}\bar w_t\|_{L_{4/3}(0,T;L_{3/2})} +
    C\|\nabla u^1\|_{L_2(0,T;L_6)} \|e^{\beta_0 t}\bar w\|_{L_4(0,T;L_2)}\\
&\leq C \bigl(\|e^{\beta_0 t}\bar w_t\|_{L_{4/3}(0,T;L_{3/2})} + \|e^{\beta_0 t}\bar w\|_{L_4(0,T;L_2)}\bigr),
\end{aligned}$$
which completes the proof of Lemma \ref{Lem:9}.\end{proof}

\smallbreak
Next, let us restate the equations \eqref{NSL-uniq} as the following system for $(z,\delta\! Q)$:
\begin{equation}\label{l21}
 \left\{\begin{aligned}
&{\varrho^1_0} z_t -\Delta_{u^1} z + \nabla_{u^1}\delta\! Q =
 (\Delta_{u^2}-\Delta_{u^1})u^2+(\nabla_{u^1}-\nabla_{u^2})Q^2
   - {\varrho^1_0} w_t + \Delta_{u^1} w - \delta\! \vr_0 \, u^2_t, \\[1ex]
&\div_{\!u^1} z=0, \qquad z|_{t=0}=\delta\! v_0.
 \end{aligned}\right.
\end{equation}
Observe that for a vector field $z$ and functions $f,g$ defined in Lagrangian coordinates we have,
according to \eqref{trans1} and integration by parts,
\begin{equation} \label{parts:div}
\begin{aligned}
-\int_{\Omega}f \, \hbox{\rm div}_u\,   z \,{\rm d}y & = -\int_{\Omega}f\div (A_u z)\,{\rm d}y
=\int_{\Omega}A_u z \cdot \nabla_y f\,{\rm d}y\\[5pt]
& = \int_{\Omega} z \cdot  A_u^T\nabla_yf \,{\rm d}y = \int_{\Omega} z \cdot \nabla_u f\,{\rm d}y,
\end{aligned}
\end{equation}
which implies
\begin{equation} \label{parts:lap}
-\int_{\Omega} f \Delta_u g \,{\rm d}y = -\int_{\Omega} f \, \hbox{\rm div}_u\,   (\nabla_u g) \,{\rm d}y
= \int_{\Omega} \nabla_u f \cdot \nabla_u g \,{\rm d}y.
\end{equation}
These identities allow to test \eqref{l21} by $z$ while
\eqref{parts:div} implies the following crucial property
thanks to which one does not have to care about the difference of the pressures:
\begin{equation}\label{l22}
 \int_{\Omega} (\nabla_{u^1} \delta\! Q)\cdot z \,{\rm d}y = -\int_{\Omega} \div_{\!u^1} z \, \delta\! Q\,{\rm d}y =0.
\end{equation}
Therefore, using also \eqref{parts:lap}, we obtain
\begin{equation}\label{l23}
 \frac 12 \frac{d}{dt} \int_{\Omega} {\varrho_0^1} |z|^2 \,{\rm d}x + \int_{\Omega} |\nabla_{u^1} z|^2\,{\rm d}x = \sum_{k=1}^5 I_k,
\end{equation}
where
$$
\begin{array}{lll}
 &\displaystyle I_1:= \int_{\Omega} \bigl((\Delta_{u^2} -\Delta_{u^1}) u^2\bigr) \cdot z \,{\rm d}y,  \quad
&  \displaystyle I_2:= \int_{\Omega}  ((\nabla_{u^1}-\nabla_{u^2})Q^2)\cdot z \,{\rm d}y,\\[3ex]
 &\displaystyle  I_3:= - \int_{\Omega} \, {\varrho^1_0} \, w_t \cdot z \,{\rm d}y, \quad & \displaystyle  I_4: = \int_{\Omega} (\Delta_{u^1} w)\cdot z\,{\rm d}y,\\[3ex]
 & \displaystyle I_5:= -\int_{\Omega} \, \delta\! \vr_0\,u_t^2 \cdot \,z \,{\rm d}y. &
\end{array}$$
In order to bound $I_1,$ we combine   \eqref{eq:regularity:1} and \eqref{l128} to write for all $\varepsilon>0,$
\begin{eqnarray}
| I_1| &\!\!\!=\!\!\!& \left| \int_{\Omega} \div\bigl(((\delta\! A)^T\!A_2+ (A_1^T)\,\delta\! A)\nabla u^2\bigr)\cdot z\,{\rm d}y \right|
\nonumber\\
 &\!\!\!\leq\!\!\!&\int_{\Omega} \bigl|(\delta\! A^T)\,A_2 + (A_1^T)\, \delta\! A\bigr|\:  |\nabla u^2 |\: |\nabla z| \,{\rm d}x
 \nonumber\\
 &\!\!\!\leq\!\!\!& C\|t^{-1/2} \delta\! A\|_2 \|t^{1/2} \nabla u^2\|_\infty \|\nabla z\|_2\nonumber\\\label{est:i1}
  &\!\!\!\leq\!\!\!&\varepsilon \|\nabla z\|^2_2 + C\varepsilon^{-1}\|t^{1/2} \nabla u^2\|_\infty^2  \biggl(\int_0^t \|\nabla \delta\! u\|_2^2\,{\rm d}\tau\biggr).
\end{eqnarray}
Next, according to \eqref{eq:regularity:4},\eqref{l128} and to Sobolev embedding
$H^1_0\hookrightarrow L_4,$
\begin{equation*
 |I_2(t)| \leq  \bigg|\int_{\Omega} \delta\! A \, \nabla Q^2 \cdot z \,{\rm d}x \bigg|\leq C\|t^{-1/2} \delta\! A\|_2 \|t^{1/2} \nabla Q^2\|_4 \|z\|_4,
\end{equation*}
Therefore, for all $\varepsilon>0,$
\begin{equation} \label{est:i2}
  I_2(t)

\leq \varepsilon\|\nabla z\|^2_2 + C\varepsilon^{-1}\|t^{1/2} \nabla Q^2\|_4^2  \biggl(\int_0^t \|\nabla \delta\! u\|_2^2 \,{\rm d}\tau\biggr)\cdotp
 \end{equation}
 Note that  from  H\"older inequality and the Sobolev embedding $H^1(\Omega)\hookrightarrow L_6(\Omega)$, we have
$$
 \|{\varrho_0}^{1/4}z\|_{{3}}\leq \|\sqrt{\varrho_0}\, z\|_{2}^{1/2}\|z\|_{6}^{1/2}
 \leq C \|\sqrt{\varrho_0}\, z\|_{2}^{1/2}\|z\|_{H^1}^{1/2}.
 $$
 Therefore using H\"older inequality, one can write that
 \begin{equation} \label{I3:5}
 \begin{aligned}
  I_3(t)& \leq \|w_t\|_{3/2}| \|\varrho_0^1 z\|_{3}
  \\&\leq C \|w_t\|_{3/2}| \|(\varrho_0^1)^{1/4} z\|_{3}  \\
  &\leq \|w_t\|_{3/2}\|\sqrt{\varrho_0^1} z\|_2^{1/2}\| z\|^{1/2}_{H^1}.
  \end{aligned}
 \end{equation}

 %

 Next, integrating by parts,  we get for all $\varepsilon>0,$
\begin{equation} \label{est:i4}
 I_4 \leq  \int_{\Omega} |\nabla_{u^1} w|\, | \nabla_{u^1} z| \,{\rm d}x
  \leq    \varepsilon \|\nabla_{u^1} z\|_2^2+ C\varepsilon^{-1}\|\nabla_{u^1} w\|_2^2.

\end{equation}

Finally, in order to estimate $I_5,$ we apply
H\"older  and Young inequality, as well as the embedding
$H^1_0\hookrightarrow L_4$ to get
\begin{align} \label{est:i5}
    I_5 &\leq \left| \int_{\Omega} \frac{\delta\! \vr_0}{\sqrt{\varrho^2_0}} \, \sqrt{\varrho_0^2} u^2_t\, z \,{\rm d}y \right|\nonumber \\
    &\leq
    \biggl\|\frac{\delta\! \vr_0}{\sqrt{\varrho^2_0}}\biggr\|_{4} \|\sqrt{\varrho^2_0} u_t^2\|_{2}\|z\|_4\nonumber\\ &\leq
    C\varepsilon^{-1}\|\sqrt{\varrho^2_0} u_t^2\|_2^2 \|\delta\! \vr_0\|^2_X + \varepsilon  \|\nabla z\|_2^2,
\end{align}
where $\|\delta\! \vr_0\|_X$ is defined in \eqref{def:xnorm}.
Let us point out that our `nice control of vacuum' hypothesis \eqref{def:xnorm} comes
into play only for handling~$I_5.$
\medbreak
In the end, plugging \eqref{est:i1}, \eqref{est:i2}, \eqref{I3:5}, \eqref{est:i4} and \eqref{est:i5} in  \eqref{l23},  we obtain
$$
\begin{aligned}
\frac12\frac{d}{dt} \|\sqrt{{\varrho_0^1}} z(t)\|_2^2 + \|\nabla_{u^1} z\|_2^2  \leq \; & \varepsilon \|\nabla z\|_2^2
+\left(\|t^{1/2}\nabla u^2\|_\infty^2+\|t^{1/2}\nabla Q^2\|_4^2\right) \int_0^t \|\nabla \delta\! u\|^2_2 \,{\rm d}t\\
&+ \|\sqrt{\varrho_0^2}u_t^2\|_2^2 \|\delta\! \vr_0\|_X^2
+C\|\nabla_{u_1}w\|_2^2
+\|w_t\|_{3/2}\|\sqrt{\varrho_0^1}z\|_2^{1/2}\|z\|_{H^1}^{1/2}.
\end{aligned}
$$
Since $\nabla u_1={}^T\,A_{u^1}\nabla z,$ we have
$$\|\nabla_{u^1}z\|_2\geq\|(A_{u^1})^{-1}\|_\infty^{-1}\|\nabla z\|_2.$$
Hence, taking $\varepsilon$ small enough, the above inequality implies
for some $c_0>0,$
\begin{multline}\label{l33}
\frac{d}{dt} \|\sqrt{{\varrho_0^1}} z(t)\|_2^2 +
c_0\|\nabla z\|_2^2
\leq  \left(\|t^{1/2}\nabla u^2\|_\infty^2+\|t^{1/2}\nabla Q^2\|_4^2\right) \int_0^t \|\nabla \delta\! u\|^2_2 \,{\rm d}t\\
+ \|\sqrt{\varrho_0^2}u_t^2\|_2^2 \|\delta\! \vr_0\|_X^2
+C\|\nabla_{u_1}w\|_2^2
+\|w_t\|_{3/2}\|\sqrt{\varrho_0^1}z\|_2^{1/2}\|z\|_{H^1}^{1/2}.
\end{multline}
Now, using the fact that $\|\nabla \delta\! u\|_2^2 \leq 2(\|\nabla z\|_2^2 + \|\nabla w\|_2^2),$ multiplying \eqref{l33} with
$e^{2\beta t}$ (with $\beta\leq\beta_0$) and adding up the second inequality of \eqref{est:w},
we arrive at
\begin{multline}\label{l33b}
\frac{d}{dt} \|e^{\beta t}\sqrt{{\varrho_0^1}} z(t)\|_2^2 + (\|e^{\beta t}\nabla z\|_2^2+\|e^{\beta t}\nabla w\|_2^2)  \leq 2\beta e^{2\beta t}\|\sqrt{{\varrho_0^1}} z(t)\|_2^2  \\[5pt]
\left(\|e^{\beta t}t^{1/2}\nabla u^2\|_\infty^2+\|e^{\beta t}t^{1/2}\nabla Q^2\|_4^2+f_1(t)\right) \int_0^T (\|\nabla z\|_2^2+\|\nabla w\|_2^2) \,{\rm d}t
+ \|e^{\beta t}\sqrt{\varrho_0^2}u_t^2\|_2^2 \|\delta\! \vr_0\|_X^2 \\[5pt]
+\|e^{\beta t}w_t\|_{3/2}\|e^{\beta t}\sqrt{\varrho_0^1}z\|_2^{1/2}\|e^{\beta t}z\|_{H^1}^{1/2}.
\end{multline}
For small $\beta$ one can absorb the first term on the right-hand side due to Poincar\'e inequality. By \eqref{eq:regularity}, provided $\beta \leq \beta_0$ we have
\begin{equation} \label{l33c1}
\|e^{\beta t}t^{1/2}\nabla u^2\|_\infty^2+\|e^{\beta t}t^{1/2}\nabla Q^2\|_4^2+\|e^{\beta t}\sqrt{\varrho_0^2}u_t^2\|_2^2 \leq f_1(t),
\end{equation}
while, by Lemma \ref{Lem:9},
$$\displaylines{
\|e^{\beta t}w_t\|_{3/2}\|e^{\beta t}\sqrt{\varrho_0^1}z\|_2^{1/2}\|e^{\beta t}z\|_{1,2}^{1/2} \hfill\cr\hfill\leq
\biggl(f_{4/3}(t) \biggl(\int_0^t \|\nabla \delta\! u\|_2^2\,{\rm d}\tau\biggr)^{1/2} +
 f_4(t)\|\nabla \delta\! u\|_2 \biggr)\|e^{\beta t}\sqrt{\varrho_0^1}z\|_2^{1/2}\|e^{\beta t}z\|_{1,2}^{1/2}.}$$

Taking advantage of Poincar\'e inequality and Lemma \ref{Lem:9}, we bound the right-hand side of the above inequality as follows:
$$\displaylines{
f_{4/3}(t) \biggl(\int_0^t\! \|\nabla \delta\! u\|_2^2 \,{\rm d}\tau\biggr)^{1/2} \|e^{\beta t}\sqrt{\varrho_0^1} z\|_2^{1/2}\|e^{\beta t}\nabla z\|^{1/2}_2 \hfill\cr\hfill\leq
 f_{4/3}^{4/3}(t) \int_0^t\! \|\nabla \delta\! u\|_2^2 \,{\rm d}\tau+ f_{4/3}^{2/3}
 \|e^{\beta t}\sqrt{\varrho_0^1} z\|_2\|e^{\beta t}\nabla z\|_2,}$$
 whence
\begin{multline} \label{l33c}
f_{4/3}(t) \biggl(\int_0^t \|\nabla \delta\! u\|_2^2 \,{\rm d}\tau\biggr)^{1/2} \|e^{\beta t}\sqrt{\varrho_0^1} z\|_2^{1/2}\|e^{\beta t}\nabla z\|^{1/2}_2 \\\leq f_1(t)\int_0^t \|\nabla \delta\! u\|_2^2\,{\rm d}\tau+ \varepsilon  \|e^{\beta t}\nabla z\|_2^2 + f_1(t) \|e^{\beta t}\sqrt{\varrho_0^1} z\|_2^2,
\end{multline}
and
\begin{equation} \label{l33d}
\begin{aligned}
    f_4(t)\|\nabla \delta\! u\|_2  \|e^{\beta t}\sqrt{\varrho_0^1} z\|_2^{1/2}\|e^{\beta t}\nabla z\|^{1/2}_2&\leq
    \varepsilon \|\nabla \delta\! u\|_2^2 +
    f_4^2(t) \|e^{\beta t}\sqrt{\varrho_0^1} z\|_2\|e^{\beta t}\nabla z\|_2\\ &\leq
    f_1(t)\|e^{\beta t}\sqrt{\varrho_0^1} z\|_2^2+
    \varepsilon (\|e^{\beta t}\nabla z\|_2^2 + \|\nabla \delta\! u\|_2^2).
\end{aligned}
\end{equation}
The $\varepsilon$ terms from \eqref{l33c} and \eqref{l33d} can again be absorbed by the left-hand side of \eqref{l33b}.
Then, plugging \eqref{l33c1}-\eqref{l33d} in \eqref{l33b} we obtain
\begin{multline} \label{l33e}
\frac{d}{dt} (e^{2\beta t}\|\sqrt{{\varrho_0}} z(t)\|_2^2) +  e^{2\beta t}(\|\nabla z\|_2^2
+\|\nabla w\|_2^2) \\ \leq
f_1(t) \left(e^{2\beta t}\|\sqrt{{\varrho_0}} z(t)\|_2^2) +  \int_0^t e^{2\beta s}(\|\nabla z\|_2^2
+\|\nabla w\|_2^2)\,{\rm d}s \right) + f_1(t) \|\delta\! \vr_0\|_X^2.
\end{multline}
Denoting
\begin{equation*}
    G(t)=e^{2\beta t}\|\sqrt{{\varrho_0}} z(t)\|_2^2 + \int_0^t e^{2\beta s}(\|\nabla z\|_2^2
+\|\nabla w\|_2^2) \,{\rm d}s,
\end{equation*}
we rewrite \eqref{l33e} as
\begin{equation} \label{eq:G}
    \frac{d}{dt}G(t) \leq f_1(t)G(t) + f_1(t)\|\delta\! \vr_0\|_X.
\end{equation}
We have $G(0)=\|\sqrt{\varrho_0}\delta\! v_0\|_2^2$, therefore \eqref{eq:G} implies
\begin{equation} \label{decay:G}
    G(t) \leq \|\sqrt{\varrho_0} \delta\! v_0\|_2^2 e^{\int_0^t f_1(\tau) \,{\rm d}\tau} +
    \|\delta\! \vr_0\|_X \int_0^t e^{\int_s^t f_1(\tau) d\tau} f(s)\,{\rm d}s.
\end{equation}
Notice that the first inequality of \eqref{est:w} implies that
\begin{equation} \label{decay:w}
\sup_{t\in \mathbb{R}_+} e^{2\beta t}\|w(t)\|_2^2 \leq f_1(t)\int_0^t \|\nabla \delta\! u(\tau)\|_2^2\,{\rm d}\tau \leq f_1(t)G(t).
\end{equation}
Combining \eqref{decay:G} and \eqref{decay:w}, we obtain the following result:
\begin{lem}
Assume that \eqref{x3} holds. Then, there exist a positive
constant $\beta<\beta_0$ and $C_0$ depending only on the data
such that
\begin{equation} \label{decay:deltau}
    \sup_{t\in [0,\infty)}
    e^{2\beta t} \|\sqrt{\varrho_0} \delta\! u \|_2^2 +\int_0^\infty
    e^{2\beta t} \|\nabla \delta\! u\|_2^2 \,{\rm d}t \leq C_0 \bigl(\|\sqrt{\varrho_0} \delta\! v_0\|_2^2 + C_2 \|\delta\! \vr_0\|_X^2\bigr)\cdotp
\end{equation}
\end{lem}

\subsection{The general case}

The estimate \eqref{decay:deltau} has been obtained
under assumption  \eqref{def:xnorm}. It may happen however that the denominator of \eqref{def:xnorm}
is zero on a subset
with positive measure, while the numerator is not. To overcome this obstacle, instead of comparing solutions
emanating from $(\varrho_0^1,v_0^1)$ and $(\varrho_0^2,v_0^2)$ directly we compare each of them to an appropriate intermediate solution satisfying \eqref{def:xnorm}.

\smallbreak
To proceed, let us denote by $(\rho^{3/2},v^{3/2},P^{3/2})$ the  velocity solution to \eqref{sys} emanating from $\left( \frac12 (\varrho_0^1+\varrho_0^2),v_0^2 \right)$ and by $(\eta^{3/2},u^{3/2},Q^{3/2})$ the corresponding solution in Lagrangian coordinates. Then, we look at the following  differences between solutions (recall that the density component of the solution in Lagrangian coordinates is constant in time):
\begin{equation*}
    \begin{aligned}
    &(\delta\! \vr^I,\delta\! u^I)=(\frac12 (\varrho_0^1+\varrho_0^2)-\varrho_0^1,\, u^{3/2}-u^1 )\\
    \quad\hbox{and}\quad&(\delta\! \vr^{II},\delta\! u^{II})=(\frac12 (\varrho_0^1+\varrho_0^2)- \varrho_0^2,\, u^{3/2}-u^2).
    \end{aligned}
\end{equation*}
Clearly,  we have
\begin{equation} \label{delta}
\delta\! \vr^I=-\delta\! \vr^{II}=\frac12 (\varrho_0^2-\varrho_0^1)\quad\hbox{and}\quad \delta\! u=\delta\! u^I-\delta\! u^{II},
\end{equation}
and  the condition \eqref{x3} is fulfilled in both cases, i.e.
\begin{equation*}
    {\rm supp}\; \delta\! \vr_0^{I,II} \subset {\rm supp}\; \frac12 (\varrho_0^1 + \varrho_0^2).
\end{equation*}
Let us define
$$
\|\delta\! \vr_0^I\|_X = \Big\| \frac{\delta\! \vr_0^I}{\sqrt{\varrho_0^1+\varrho_0^2}}\Big\|_4=\Big\|\frac{\delta\! \vr_0^{II}}{\sqrt{\varrho_0^1+\varrho_0^2}}\Big\|_4\cdotp
$$
Note that, obviously
\begin{equation*}
    |\frac{\varrho_0^1-\varrho_0^2}{\sqrt{\varrho_0^1+\varrho_0^2}}|^2 =
    |\varrho_0^1-\varrho_0^2| |\frac{\varrho_0^1-\varrho_0^2}{\varrho_0^1+\varrho_0^2}|
    \leq |\varrho_0^1-\varrho_0^2|,
\end{equation*}
whence
\begin{equation} \label{dvr}
\|\delta\! \vr_0^I\|_X \leq \|\varrho_0^1-\varrho_0^2\|_2^{1/2}.
\end{equation}
Denoting $\delta\! Q^I=Q^{3/2}-Q^1$,  we obtain the following system for $(\delta\! u^I,\delta\! Q^I)$:
\begin{equation}\label{NSL-uniq:2}
\left\{ \begin{aligned}
  &{\varrho_0^1}\delta\! u^I_t - \Delta_{u^1}\delta\! u^I + \nabla_{u^1}\delta\! Q^I=
(\Delta_{u^{3/2}} -\Delta_{u^1}) u^{3/2} - (\nabla_{u^{3/2}}-\nabla_{u^1}) Q^{3/2}-\delta\! \vr^I \,u^{3/2}_t,\\
&\div_{\!u^1}\delta\! u^I=(\div_{\!u^1}-\div_{\!u^{3/2}})u^{3/2},\\
&\delta\! u^I|_{t=0}=\delta\! v_0.
\end{aligned}\right.
\end{equation}
Setting $\delta\! Q^{II}=Q^{3/2}-Q^2,$ we see that $(\delta\! u^{II},\delta\! Q^{II})$ satisfies an analogous system with $(\varrho^1_0,u^1,\delta\! \vr^I)$ replaced by $(\varrho^2_0,u^2,\delta\! \vr^{II})$
and with the initial condition $\delta\! u^{II}|_{t=0}=0$.
\smallbreak
Let us consider the decompositions
$$
\delta\! u^I = z^I+w^I\quad\hbox{and}\quad\delta\! u^{II} = z^{II}+w^{II},
$$
where the components are defined by \eqref{l6} and \eqref{l21} with obvious replacements of $u^1,u^2$.
Then, one can repeat the estimates from the previous subsection. In the end, defining
\begin{equation} \label{def:G12}
\begin{aligned}
    &G^I(t)=e^{2\beta t}\|\sqrt{{\varrho_0^1}} z^I(t)\|_2^2 + \int_0^t e^{2\beta s} (\|\nabla z^I(s)\|_2^2
+\|\nabla w^I(s)\|_2^2)\,{\rm d}s,\\
    &G^{II}(t)=e^{2\beta t}\|\sqrt{{\varrho_0^2}} z^{II}(t)\|_2^2 + \int_0^t e^{2\beta s}(\|\nabla z^{II}(s)\|_2^2
+\|\nabla w^{II}(s)\|_2^2)\,{\rm d}s,
\end{aligned}
\end{equation}
one obtains the following analogs of \eqref{decay:G} (recall that $\delta\! \vr_0^I=-\delta\! \vr_0^{II}$ and that the velocity component of the initial data for $\delta\! v^{II}$ is zero) for some  $f^I,f^{II} \in L^1(\mathbb{R}_+)$:
$$
\begin{aligned}
     G^I(t) &\leq \|\sqrt{\varrho_0^1} \delta\! v_0\|_2^2 e^{\int_0^tf^I(\tau) \,{\rm d}\tau} +
    \|\delta\! \vr_0^I\|_X^2 \int_0^t e^{\int_s^t f^I(\tau) d\tau} f^I(s)  \,{\rm d}s,\\
     G^{II}(t) &\leq
    \|\delta\! \vr_0^{I}\|_X^2 \int_0^t e^{\int_s^t f^{II}(\tau)\, d\tau} f^{II}(s) \,{\rm d}s,
\end{aligned}    $$
Summing  up the above inequalities and using \eqref{delta} and \eqref{dvr} we obtain
\begin{equation*}
\sup\limits_{t\in\mathbb{R}_+} \left\{ e^{\beta t}\left[\|\sqrt{\varrho^1_0}z^I(t)\|_2+\|\sqrt{\varrho^2_0}z^{II}(t)\|_2 \right] \right\}
+ \|e^{\beta t}\nabla \delta\! u\|_{L^2(\mathbb{R}_+,L^2)}\leq C \left( \|\sqrt{\varrho^1_0}\delta\! v_0\|_2 + \|\delta\! \vr_0\|_2^{1/2} \right)\cdotp
\end{equation*}
Combining this estimate with analogs of \eqref{decay:w}, where $(w,\delta\! u,G)$ has been replaced by
$(w^{I},\delta\! u^{I},G^{I})$ and $(w^{II},\delta\! u^{II},G^{II}),$ we obtain \eqref{decay:final}, which completes the proof of Theorem \ref{thm:lag}.

\subsection{Back to the Euler perspective}

In this part, we want to translate the stability result obtained in Theorem \ref{thm:lag}  in terms of Eulerian coordinates, proving Corollary \ref{cor:eul}.
As already pointed out in the introduction, in the case of only
bounded initial densities,
getting a relevant  information
on $\varrho^1(t,\cdot)-\varrho^2(t,\cdot)$ in some $L^p$ space (even for $p=1$)
is hopeless; it is more reasonable to compare  functions along two different characteristics fields.
Here we want to adopt  the language of kinetic theory,
 considering quantities along characteristics/trajectories,
defined in  terms of Wasserstein metrics like in e.g. \cite{PR,MPesz}.
\smallbreak
More precisely, denote by $X^1$ and $X^2$ the flows defined by \eqref{L-E} for, respectively, $v^1$ and $v^2,$ and take a smooth function $\phi:\Omega \to \mathbb{R}.$
Then we consider, for each $t\geq0,$ the quantity
$$I_\phi(t):= \int_\Omega \Big(\varrho^1(t,x) -  \varrho^2(t,x)\Big)\phi (x) \,{\rm d}x.$$
Then, using the change of variables (of Jacobian $1$) $x=X^1(t,y)$ and $x=X^2(t,y)$
for $\varrho^1$ and $\varrho^2,$ respectively, yields

\begin{align} \label{eq:deltarho}
I_\phi(t)&=
    \int_\Omega \Big[\varrho_0^1 (y) \phi(X^1(t,y)) -  \varrho_0^2(y) \phi (X^2(t,y)) \Big]\,{\rm d}y\nonumber\\
    &=\underbrace{\int_\Omega (\varrho_0^1(y)- \varrho_0^2(y)) \phi(X^1(t,y)) \,{\rm d}y}_{A_1(t)} +
    \underbrace{\int_\Omega  \varrho_0^2(y) (\phi(X^1(t,y)) -\phi(X^2(t,y)))\,{\rm d}y}_{A_2(t)}.
\end{align}
In order to find out the right level of regularity for $\phi,$ let us first
examine the term $A_2(t).$
For suitable functions $\underline{\eta},$ we set
$$\eta^1(t,x):=\underline{\eta} (t,Y^1(t,x)) \quad\hbox{and}\quad \eta^2(t,x):=\underline{\eta} (t,Y^2(t,x))\quad\hbox{with}\quad Y^i(t,\cdot):=(X^i(t,\cdot))^{-1}.$$
Then, using again the fact that $X^1$ and $X^2$ are measure preserving,
we have
\begin{equation} \label{eq:neg}
\begin{aligned}
    &\|\eta^1(t) - \eta^2(t)\|_{W^{-1}_p} =\sup \left\{
    \int_\Omega (\eta^1(t,x)-\eta^2(t,x))\phi(t,x)\,{\rm d}x: \phi \in W^1_{p'}(\Omega), \;
    \|\phi\|_{W^1_{p'}} \leq 1
    \right\}\\[5pt]
    &= \sup
    \left\{ \int_\Omega \underline{\eta}(y) [\phi(X^1(t,y))-\phi(X^2(t,y))]\,{\rm d}y : \phi \in W^1_{p'}(\Omega), \;
    \|\phi\|_{W^1_{p'}} \leq 1 \right\}\cdotp
\end{aligned}
\end{equation}
Following \cite{Maja}, in order to handle the term $\phi(X^1(t,y)) -\phi(X^2(t,y)),$ we consider
the  family of  intermediate measure preserving flows  $(X^s)_{1\leq s\leq2}$
between $X^1$ and $X^2$ defined as follows:
\begin{equation*}
    \frac{dX^s}{dt}(t,y)= (2-s) v^1(t,X^s(t,y)) + (s-1) v^2 (t,X^s(t,y)), \qquad
    X^s(0,y) =y.
\end{equation*}
By the chain rule and the definition of $X^s,$ we have
$$\begin{aligned
  \phi(X^2(t,y)) -\phi(X^1(t,y))&= \int_1^2 \frac{d}{ds} \phi(X^s(t,y)) \,{\rm d}s\nonumber\\  &=

  \int_1^2 \biggl(\frac d{ds} X^s(t,y)\biggr)\cdotp \nabla \phi (X^s(t,y)) \,{\rm d}s.
\end{aligned}$$
{}From the definition of $X^s,$ we discover that
$$\frac d{ds} X^s(t,y)=\int_0^t\bigl(v^2(t',X^s(t',y))-v^1(t',X^s(t',y))\bigr)\,{\rm d}t',$$
which implies that
\begin{equation}\label{eq:Xs}
\biggl\|\frac d{ds} X^s(t,\cdot)\biggr\|_p\leq \int_0^t\|\delta\! v\|_p\,{\rm d}t'.
\end{equation}
Furthermore, as said before,  $X^s(t,\cdot)$ is
measure preserving  for all $t\geq 0,$ and thus
    $\| \nabla \phi (X^s(t,\cdot)) \|_p = \|\nabla \phi\|_p.$

In the end, using the embedding $H^1_0\hookrightarrow L_p$ for all $1\leq p<\infty$ if
$d=2,$ and all $1\leq p\leq 6$ if $d=3,$ we obtain for these values of $p$ and all $t\geq0,$
\begin{equation} \label{est:A2}
    A_2(t)\leq

    C\|\varrho_0^2\|_\infty  \left( \int_0^\infty e^{2\beta t} \| \nabla \delta\! v\|_2^2 \,dt\right)^{1/2}\|\phi\|_{1,p'}.
\end{equation}
The above inequality reveals that one can take  the functions $\phi$ in the space
\begin{equation*}
    \phi \in W^{1}_{1^+}(\Omega) \mbox{ \ \ if \ \ } d=2
\quad\hbox{and}\quad
    \phi \in W^1_{6/5}(\Omega) \mbox{ \ \ if \ \ } d=3.
\end{equation*}
Bounding the term $A_1(t)$ is simpler.
Under the above assumptions on $p,$ we have $W^1_{p'}\hookrightarrow L_2.$
Hence,  using again that $X^1$ is measure preserving, we may write by Cauchy-Schwarz inequality,
\begin{equation} \label{est:A1}
    A_1(t)
    \leq  \|\delta\! \vr_0\|_2 \|\phi \|_2 \leq C \|\delta\! \vr_0\|_2 \|\phi \|_{1,p'}.
\end{equation}
Altogether, by \eqref{eq:deltarho}, \eqref{eq:neg}, \eqref{est:A2}, \eqref{est:A1} and \eqref{decay:final},  we get
for any finite $p>1$ if $d=2$ and for $p\leq6$ if $d=3,$
\begin{equation} \label{est:deltarho:final}
   \sup_{t\in\mathbb{R}_+} \|\varrho^1(t,\cdot) - \varrho^2(t,\cdot)\|_{W^{-1}_p} \leq C\biggl(
 \|\varrho_0^2\|_\infty  \left( \int_0^\infty e^{2\beta t} \| \nabla \delta\! v\|_2^2 \,dt\right)^{1/2} +\|\delta\! \vr_0\|_2\biggr)\cdotp
\end{equation}
Next, let us turn to the velocity estimate. We start from the relation
$$
\delta\! u (t,y) = v^2(t,X^2(t,y))-v^1(t,X^1(t,y)),
$$
so
$$    \nabla_y \delta\! u (t,y)=  \nabla_yX^2(t,y) \cdotp\nabla_x v^2(t,X^2(t,y))
    - \nabla_yX^1(t,y)\cdotp\nabla_x v^1(t,X^1(t,y)).$$
Hence, denoting $\delta\! X:=X^2-X^1,$ we may write
$$\begin{aligned}
    \nabla_y \delta\! u (t,y
&=\nabla_y\delta\! X(t,y)\cdotp\nabla v^2(t,X^2(t,y))+
\nabla_yX^1(t,y)\cdotp\nabla \delta\! v (t,X^2(t,y))\\&\hspace{5cm}+
\nabla_yX^1(t,y)\cdotp\left[\nabla v^1 (t,X^2(t,y))-\nabla v^1(t,X^1(t,y)) \right]\\&=:K_1+K_2+K_3.
\end{aligned}$$
First, we observe that
$$\nabla\delta\! v(t,X^2(t,y))={}^T\!A_u^1 K_2.$$
Hence, taking the $L_2(\Omega)$ norm, and using that $X^2(t,\cdot)$ is measure preserving, we get
\begin{equation}\label{eq:dv}
 \|\nabla\delta\! v\|_{2}\leq \|A_u^1\|_\infty  \|K_2\|_{2}\leq \|A_u^1\|_\infty\bigl(\|\nabla\delta\! u\|_2
 +\|K_1\|_2+\|K_3\|_2\bigr)\cdotp
\end{equation}
Next, from the definition of $X^1$ and $X^2$ in terms of the Lagrangian velocity, we have
\begin{equation*}
    \|K_1\|_2\leq Ct^{1/2} \left(\int_0^t \|\nabla\delta\! u\|_2^2\,dt'\right)^{1/2} \|\nabla v^2\|_\infty.
\end{equation*}
For bounding $K_3,$ we first notice that by the mean value theorem and the definition
of the intermediate flow $X^s$, we have
$$K_3=\nabla_yX^1(t,y)\cdotp\biggl(\int_1^2\nabla^2 v^1(t,X^s(t,y))\,{\rm d}s\biggr)\cdotp
\biggl(\frac d{ds}X^s(t,y)\biggr)\cdotp$$
Hence, since $X^s$ is measure preserving, using H\"older inequality and \eqref{eq:Xs} allows to get
$$\begin{aligned}
    \|K_3\|_2 &\leq \|\nabla_yX^1(t)\|_\infty \int_1^2\bigl\|\nabla^2 v^1(t,X^s(t))\bigr\|_3
    \Bigl\|\frac d{ds}X^s(t)\Bigr\|_6\,{\rm d}s\\
&    \leq C\|\nabla_yX^1(t)\|_\infty
\|\nabla^2 v^1(t)\|_3 \int_0^t \|\delta\! v(t')\|_6\,dt'\\
    &\leq C\|\nabla_yX^1(t)\|_\infty
    \|t^{1/2}\nabla^2v^1(t) \|_3 \left(\int_0^t \|\nabla \delta\! v\|_2^2\,dt'\right)^{1/2}.
\end{aligned}$$
Altogether, remembering \eqref{eq:dv}, we obtain for all $t\geq0,$
$$
\|\nabla\delta\! v(t)\|_2^2 \leq C \left[ \|\nabla \delta\! u(t)\|_2^2 +\|t^{1/2}\nabla v^2\|_\infty^2\int_0^t \|\nabla \delta\! u\|_2^2\,{\rm d}t'
+\|t^{1/2}\nabla^2 v^1\|_3^2\int_0^t\|\nabla \delta\! v\|_2^2\,{\rm d}t'\right]\cdotp
$$
First multiplying both sides by $e^{2\beta t}$, next integrating in time and
using Theorem \ref{thm:lag} combined with the properties \eqref{eq:regularity}, and finally applying
Gronwall lemma, we obtain for all $t\geq0,$
\begin{equation*}
    \int_0^t e^{2\beta t'} \|\nabla\delta\! v\|_2^2\,{\rm d}t'
    \leq C_0 \left( \|\sqrt{\varrho^1_0}\delta\! v_0\|_2^2 + \|\delta\! \vr_0\|_2 \right)
    \exp\left(C\int_0^t e^{2\beta t'}\|t'^{1/2}\nabla^2v^1\|_3^2\,{\rm d}t'\right)\cdotp
\end{equation*}
The term in the exponential may be bounded thanks to \eqref{eq:regularity:4}.
Hence, reverting to  \eqref{est:deltarho:final}
completes the proof of Corollary \ref{cor:eul}.

\section{Decay estimates} \label{s:decay}

The goal of this section is to show that the solutions provided by Theorems \ref{thm:d=2} and \ref{thm:d=3}
indeed satisfy Properties \eqref{eq:regularity:8} to \eqref{eq:regularity:7}. This will be an immediate
consequence of Lemmas \ref{l:decay1}, \ref{l:decay2} and \ref{l:decay3} below.
\smallbreak
Before tackling the proof, we need to recall elementary properties
of the  solutions to \eqref{sys}. The first one is the  conservation of any $L_p$ norm of the density,
a consequence of the divergence free property of the velocity field : we have
$\|\varrho(t)\|_p = \|\varrho_0\|_p$.

Of course, as $v$ in $H^1_0(\Omega),$ we have the Poincar\'e inequality:
\begin{equation} \label{Poin}
\|v\|_{2} \leq C_P \|\nabla v\|_2.
\end{equation}

\subsection{Decay of space derivatives}
The starting point is the following decay estimate which is a direct consequence of
diffusion and  Poincar\'e inequality:
\begin{lem}
Let $(\varrho,v)$ be a solution to \eqref{sys} given either by Theorem \ref{thm:d=2} or by Theorem \ref{thm:d=3}. Then
\begin{equation} \label{decay:1}
\forall t\geq0,\;
\int_{\Omega} \varrho(t)|v(t)|^2\,{\rm d}x \leq e^{-\beta_1 t}\int_{\Omega} \varrho_0|v_0|^2\,{\rm d}x,
\mbox{ \ \ where  \ \ } \beta_1 =\frac{2\mu}{\varrho^* C_P^2}\cdotp
\end{equation}
\end{lem}
\begin{proof} The proof is independent of the space dimension. We start with the classical energy estimate that is obtained testing the momentum equation by $v,$ namely
\begin{equation} \label{ene:1}
\frac12\frac{d}{d t}\int_{\Omega} \varrho|v|^2\,{\rm d}x + \mu\int_{\Omega}|\nabla v|^2\,{\rm d}x = 0.
\end{equation}
Hence, remembering  \eqref{Poin}, we get
$$\frac{C_P^2}{2\mu} \frac{d}{d t}\int_{\Omega}\varrho|v|^2\,{\rm d}x+\int_{\Omega}|v|^2\,{\rm d}x \leq 0.$$
Multiplying both sides by $\varrho^*$ and using  the obvious fact that
$\varrho^*|v|^2 \geq \varrho |v|^2,$
we obtain
$$
\frac{\varrho^*C_P^2}{2\mu} \frac{d}{d t}\int_{\Omega}\varrho|v|^2\,{\rm d}x+\int_{\Omega}\varrho|v|^2\,{\rm d}x \leq 0,
$$
from which we conclude to \eqref{decay:1}.
\end{proof}

Our next aim is to establish a  decay estimate for the gradient of the velocity.
\begin{lem} \label{l:decay1}
Let $(\varrho,v)$ be a solution to \eqref{sys} given either by Theorem \ref{thm:d=2} or by Theorem \ref{thm:d=3}. Then, there exist  positive constants  $C_0$ and $C_{0,p}$ depending only on the data
(and on $p$ for second one),
and\footnote{One can take $\beta_2$ of the
form $\beta_2=c_\Omega \mu /(\varrho^*\delta^2)$ where $\delta$ stands for the
diameter of $\Omega$ and $c$ depends only on the shape of $\Omega.$}
$\beta_2<\beta_1$  such that for all $t\geq0,$ we have
\begin{align}
&\|\nabla v(t)\|_2 \leq C_0e^{-\beta_2 t}, \label{decay:2}\\[3pt]
&\|v(t)\|_p \leq C_{0,p}e^{-\beta_2 t}, \label{decay:2c}\\[3pt]
&\int_0^{\infty} e^{2\beta_2 t}\left[\|\sqrt{\varrho(t)} v_t(t)\|_2^2 + \|\nabla^2 v(t)\|_2^2 + \|\nabla P(t)\|_2^2\right]\,{\rm d}t \leq C_0,
\label{decay:2a}
\end{align}
where, in \eqref{decay:2c}, one can take any $p\in[1,\infty)$ if $d=2,$
and any $p\in[1,6]$ if  $d=3.$
\end{lem}

\begin{proof}
Clearly, performing a suitable time, space, density and velocity rescaling
reduces the proof to the case $\varrho^*=\mu=1,$ in a domain of diameter $1.$
Hence, we shall only prove the lemma in this case for notational simplicity.

Now, testing the momentum equation by $v_t$ yields
$$\frac12\frac{d}{d t}\int_{\Omega} |\nabla v|^2\,{\rm d}x
+\int_{\Omega}\varrho |v_t|^2 \,{\rm d}x =-\int_\Omega \varrho v_t\cdot(v\cdot\nabla v)\,{\rm d}x\leq \frac12\int_{\Omega}\varrho |v_t|^2 \,{\rm d}x
+\frac12\int_{\Omega}\varrho|v\cdot\nabla v|^2\,{\rm d}x.
$$
The classical theory for the Stokes system yields for some $C_\Omega$ depending only
on the shape of $\Omega,$
$$\|\nabla^2 v\|_2^2+\|\nabla P\|_2^2\leq C_\Omega \bigl(\|\varrho v\cdot\nabla v\|_2^2+ \|\varrho v_t\|_2^2\bigr)\cdotp$$
Hence we have (remember that $\varrho^*=1$):
\begin{equation} \label{ene:2}
\frac{d}{d t}\int_{\Omega} |\nabla v|^2\,{\rm d}x+\frac12\int_{\Omega}\bigl(\varrho|v_t|^2 +
\frac1{C_\Omega}(|\nabla^2 v|^2+|\nabla P|^2)\bigr)\,{\rm d}x
\leq \frac32\int_{\Omega}\varrho|v\cdot\nabla v|^2\,{\rm d}x.
\end{equation}
Adding this inequality to \eqref{ene:1}, we get (up to a change of $C_\Omega$)
\begin{multline} \label{2.2}
\frac{d}{d t}\int_{\Omega}[|\nabla v|^2+\frac{1}{2}\varrho|v|^2]\,{\rm d}x +\int_{\Omega}[|\nabla v|^2+\frac{1}{2}\varrho|v|^2]\,{\rm d}x
+\frac12\int_{\Omega}\bigl(\varrho|v_t|^2 +\frac1{C_\Omega}(|\nabla^2 v|^2+|\nabla P|^2)\bigr)\,{\rm d}x\\
\leq \frac32\int_{\Omega}\varrho|v\cdot\nabla v|^2\,{\rm d}x + \frac{1}{2}\int_{\Omega}\varrho|v|^2\,{\rm d}x.
\end{multline}
In order to estimate the first term on the right in the two dimensional case we first write
$$\int_{\Omega}\varrho|v\cdot\nabla v|^2\,{\rm d}x
\leq \|\sqrt{\varrho}v\|_2\|v\|_\infty\|\nabla v\|_4^2.$$
Hence, using the following two interpolation inequalities:
\begin{equation} \label{Lad}
\|z\|_4 \leq C \|z\|_2^{1/2}\|\nabla z\|_2^{1/2}\ \hbox{ and }\
\|z\|_\infty\leq C\|z\|_2^{1/2}\|\nabla^2z\|_2^{1/2},
\end{equation}
and remembering that $\varrho^*=1,$ we get
$$\begin{aligned}
\frac32\int_{\Omega}\varrho|v\cdot\nabla v|^2\,{\rm d}x &\leq C \|\sqrt{\varrho}v\|_2\|v\|_2^{1/2}\|\nabla v\|_2\|\nabla^2v\|_2^{3/2}\\
&\leq \frac1{4C_\Omega}\|\nabla^2v\|_2^2 + C
\|\sqrt{\varrho}v\|_2^4\|v\|_2^{2}\|\nabla v\|_2^4.
\end{aligned}$$
The first term can be absorbed by the left-hand side of \eqref{2.2} and
one can bound $\|v\|_2$ from \eqref{Poin}. Hence,  we get
for some $\gamma>0$ depending only on the shape of the domain,
\begin{multline} \label{2.6}
\frac{d}{d t}\int_{\Omega}[|\nabla v|^2+\frac{1}{2}\varrho|v|^2]\,{\rm d}x +\int_{\Omega}[|\nabla v|^2+\frac{1}{2}\varrho|v|^2]\,{\rm d}x
+\frac12\int_{\Omega}\bigl(\varrho|v_t|^2 +\frac1{2C_\Omega}(|\nabla^2 v|^2+|\nabla P|^2)\bigr)\,{\rm d}x\\
\leq Ce^{-\gamma t} [1 + \|\nabla v\|_2^6].\end{multline}
The crucial observation is that $\|\nabla v(t)\|_2$ may  bounded uniformly in time in terms
of the data (see Propositions 3.1 and 3.3 in \cite{DM1}). Therefore, we obtain
\begin{equation} \label{2.6ba}
\frac{d}{d t}\int_{\Omega}[|\nabla v|^2+\frac{1}{2}\varrho|v|^2]\,{\rm d}x
+\frac12 \int_{\Omega}\bigl(|\nabla v|^2+\frac{1}{2}\varrho|v|^2+\varrho|v_t|^2
+\frac1{2C_\Omega}(|\nabla^2 v|^2+|\nabla P|^2)\bigr)\,{\rm d}x
\leq C_0 e^{-\gamma t}.
\end{equation}
This leads for any $\gamma_1 \in\mathbb{R}$ to
\begin{multline} \label{2.6b}
\frac{d}{d t} \left[ e^{\gamma_1 t}\int_{\Omega}\left(|\nabla v|^2+\frac{1}{2}\varrho|v|^2\right)\,{\rm d}x\right]\\
 +\Bigl(\frac12-\gamma_1\Bigr) \int_{\Omega}e^{{\gamma_1} t}\bigl(|\nabla v|^2+\frac{1}{2}\varrho|v|^2+\varrho|v_t|^2
+\frac1{2C_\Omega}(|\nabla^2 v|^2+|\nabla P|^2)\bigr)\,{\rm d}x
\leq C_0 e^{-(\gamma -\gamma_1)t}.
\end{multline}

In the three dimensional case, the only difference is the slightly more complicated treatment of the right-hand side of \eqref{2.2}. Nevertheless, by using H\"older inequality and Gagliardo-Nirenberg inequality,
we arrive (still using that $\varrho^*=1$) at:
$$\begin{aligned}
\int_{\Omega}\varrho |v\cdot\nabla v|^2 &\leq
\|\varrho^{1/4} v\|_4^2 \|\nabla v\|_4^2\\
&\leq
C\|\sqrt{\varrho}v\|_2^{1/2}\|v\|_6^{3/2}\|\nabla v\|_2^{1/2}\|\nabla^2 v\|_2^{3/2}\\
&\leq C\|\sqrt{\varrho}v\|_2^{1/2}\|\nabla v\|_2^2\|\nabla^2 v\|_2^{3/2} \\
&\leq \frac{1}{2C_\Omega}\|\nabla^2 v\|_2^2 + C\|\sqrt{\varrho}v\|_2^2\|\nabla v\|_2^8 \leq
\frac{1}{2C_\Omega}\|\nabla^2 v\|_2^2 + C_0e^{-\beta_1 t}, \end{aligned} $$
where in the last passage we have used \eqref{decay:1} and uniform boundedness of $\|\nabla v(t)\|_2$.
Again, the first term can be absorbed by the left hand side of \eqref{2.2} and we obtain \eqref{2.6b}.

Now, choosing e.g. $\gamma_1=\min(\gamma,1/4)$ and integrating \eqref{2.6b} on $[0,t]$ yields
$$\displaylines{e^{\gamma_1t}\Bigl(\|\nabla v(t)\|_2^2+\frac12\|\sqrt{\varrho(t)}\,v(t)\|_2^2\Bigr)
+\Bigl(\frac12-\gamma_1\Bigr)\int_0^te^{\gamma_1s}\Bigl(\|\nabla v\|_2^2+\frac12\|\sqrt\varrho\,v\|_2^2
\hfill\cr\hfill+\frac12\|\sqrt\varrho\,v_t\|_2^2+\frac1{2C_\Omega}(\|\nabla^2 v\|_2^2+\|\nabla P\|_2^2)\Bigr)\,{\rm d}s
\leq\frac{C_0}{\gamma-\gamma_1}
+\|\nabla v_0\|_2^2+\frac12\|\sqrt{\varrho_0}\,v_0\|_2^2,}$$
which readily gives \eqref{decay:2} and \eqref{decay:2a}.
As for \eqref{decay:2c}, it just
results from Poincar\'e inequality and Sobolev embedding.

\end{proof}

\subsection{Decay of time derivatives}
Here we  estimate  time and time-space derivatives.
\begin{lem}\label{l:decay2}
Let $(\varrho,v)$ be a solution to \eqref{sys} given either by Theorem \ref{thm:d=2} or by Theorem \ref{thm:d=3}. Then, there exists a positive constant $\beta_3<\beta_2$ (still of the form $\beta_3=c \mu/(\varrho^*\delta^2)$)
such that
\begin{equation} \label{decay:3}
\sup_{t\in \mathbb{R}_+}\int_{\Omega}te^{2\beta_3  t}\varrho|v_t|^2\,{\rm d}x +  \int_0^{\infty} e^{2\beta_3  t}\left[ \|\sqrt{t}\nabla v_t\|_2^2+
\|\sqrt{t} v_t\|_2^2
\right]\,{\rm d}t \leq C_0
\end{equation}
for some positive constant $C_0$ depending only on the data.
\end{lem}
\begin{proof}  We keep the assumption $\delta=\mu=\varrho^*=1.$
Now, differentiating the momentum equation in time and multiplying by $\sqrt{t}e^{\beta  t},$ we obtain
\begin{multline}\label{eq:ppp}
\varrho(\sqrt{t}e^{\beta  t}v_t)_t +\sqrt{t}e^{\beta  t}\varrho_tv_t
-\frac{1}{2\sqrt{t}}e^{\beta  t}\varrho v_t -
\beta  \sqrt{t}e^{\beta  t}\varrho v_t+\sqrt{t}e^{\beta  t}\varrho_tv\cdot\nabla v \\+\sqrt{t}e^{\beta  t}\varrho v_t\cdot\nabla v +\sqrt{t}e^{\beta  t}\varrho v\cdot\nabla v_t -\Delta(\sqrt{t}e^{\beta  t}v_t)+\nabla (\sqrt{t}e^{\beta  t}P_t)=0.\end{multline}
Testing \eqref{eq:ppp}  with $\sqrt{t}e^{\beta  t}v_t$ yields
\begin{equation} \label{2.7}
\frac{1}{2}\frac{d}{d t}\int_{\Omega}te^{2\beta  t}\varrho |v_t|^2\,{\rm d}x+\int_{\Omega}te^{2\beta  t}|\nabla v_t|^2 = \sum_{i=1}^5 I_i,
\end{equation}
where, using $\varrho_t=-\div(\varrho v)$ in the third term on the left-hand side of \eqref{eq:ppp}, we have
\begin{align*}
&I_1 = \frac{1}{2}\int_{\Omega}e^{2\beta  t}\varrho|v_t|^2\,{\rm d}x,
&I_2 = -\int_{\Omega}(\sqrt{t}e^{\beta  t}\varrho_t v \cdot \nabla v) \cdot(\sqrt{t}e^{\beta  t} v_t)\,{\rm d}x,\\
&I_3 = -\int_{\Omega}(\sqrt{t}e^{\beta  t}\varrho v_t \cdot \nabla v) \cdot(\sqrt{t}e^{\beta  t} v_t) \,{\rm d}x,
&I_4 = -2\int_{\Omega}(\sqrt{t}e^{\beta  t}\varrho v \cdot \nabla v_t) \cdot(\sqrt{t}e^{\beta  t} v_t) \,{\rm d}x,\\
&I_5 = \beta  \int_\Omega te^{2\beta  t}\varrho |v_t|^2. &
\end{align*}

Now, we estimate the right-hand side of \eqref{2.7}. From the continuity equation, we have
$$I_2 = \int_{\Omega} te^{2\beta  t}\div(\varrho v) (v \cdot \nabla v) \cdot v_t \,{\rm d}x =
- \int_{\Omega} t  e^{2\beta  t} \varrho v \cdot \nabla[(v\cdot\nabla v)\cdot v_t]\,{\rm d}x.$$
Therefore,
$$|I_2| \leq \int_{\Omega}t e^{2\beta  t} \varrho|v|\left( |\nabla v|^2 |v_t| + |v||\nabla^2 v||v_t| + |v||\nabla v||\nabla v_t| \right)\,{\rm d}x =: I_{21}+I_{22}+I_{23}.  $$
For $I_{21}$ in the two-dimensional case, we have
$$\begin{aligned}
|I_{21}| &\leq \int_\Omega\bigl(\sqrt{t\varrho}|v||\nabla v|^2\sqrt{t\varrho}|v_t|e^{2\beta  t}\bigr)\,{\rm d}x\\
&\leq \|v\|_{\infty}^2\|\sqrt{t\varrho} e^{\beta  t}v_t\|_2^2 + t e^{2\beta  t}\|\nabla v\|_4^4,\end{aligned}
$$
so using \eqref{Lad} and \eqref{decay:2} we get
\begin{eqnarray} \label{2.27}
|I_{21}| &\!\!\!\leq\!\!\!& \|v\|_{\infty}^2\|\sqrt{t\varrho}e^{\beta  t}v_t\|_2^2 + t\|\nabla v\|_2^2 e^{2\beta  t}\|\nabla^2 v\|_2^2\nonumber\\
&\!\!\!\leq\!\!\!& \|v\|_{\infty}^2\|\sqrt{t\varrho}e^{\beta  t}v_t\|_2^2 + Ce^{-2\beta_2 t}
\|e^{\beta  t}\nabla^2 v\|_2^2.
\end{eqnarray}

In the three dimensional case we do not have \eqref{Lad}, but we can write using H\"older inequality
and the embedding $H^1_0\hookrightarrow L_6,$
$$\begin{aligned}
|I_{21}| & \leq \sqrt{t} e^{2\beta  t} \int_{\Omega}\sqrt{t\varrho}|v_t||v||\nabla v|^2\,{\rm d}x\\ &\leq
\sqrt{t} e^{2\beta  t} \|\sqrt{t\varrho}v_t\|_4 \|v\|_6 \|\nabla v\|_{24/7}^2 \\
&\leq \sqrt{t} e^{\beta t} \|\sqrt{t\varrho}v_t\|_2^{1/4}\|\sqrt{t}e^{\beta t}v_t\|_6^{3/4}\|v\|_6\|\nabla v\|_{24/7}^2\\
&\leq \frac{1}{10}\|\sqrt{t}e^{\beta  t}\nabla v_t\|_2^2 + Ct^{4/5}e^{8\beta t/5} \|\sqrt{t\varrho}v_t\|_2^{2/5}\|\nabla v\|_{24/7}^{16/5}\|\nabla v\|_2^{8/5}.
\end{aligned}$$
Then, using the Gagliardo-Nirenberg inequality
$\|\nabla v\|_{24/7}^{16/5}\leq C \|\nabla v\|_2^{6/5}\|\nabla^2 v\|_2^2$ and \eqref{decay:2},
we discover that
\begin{equation*}
\begin{aligned}
|I_{21}| &\leq \frac{1}{10}\|\sqrt{t}e^{\beta  t}\nabla v_t\|_2^2 + t^{4/5}e^{-4\beta t/5}\|\sqrt{t\varrho} e^{\beta  t}v_t\|_2^{2/5} \|\nabla v\|_2^{14/5}\|e^{\beta  t}\nabla^2 v\|_2^2\\
&\leq \frac{1}{15}\|\sqrt{t}e^{\beta  t}\nabla v_t\|_2^2 + C_0e^{-(4\beta+14\beta_2)t/5}
\|e^{\beta  t}\nabla^2 v\|_2^2\|\sqrt{t\varrho}e^{\beta  t}v_t\|_2^{2/5}.
\end{aligned}
\end{equation*}
Therefore, there exits $c>0$ such that
\begin{equation} \label{I21}
|I_{21}|\leq \frac{1}{10}\|\sqrt{t}e^{\beta  t}\nabla v_t\|_2^2 + C_0e^{-ct}\|e^{\beta  t}\nabla^2 v\|_2^2\left(1+\|\sqrt{t\varrho}e^{\beta  t}v_t\|_2^2 \right).
\end{equation}
The remaining two parts of $I_2$ are simpler: we have
\begin{eqnarray} \label{I22}
I_{22}=\int_{\Omega}\sqrt{t\varrho} e^{\beta  t} |\nabla^2 v||v|^2\sqrt{t\varrho}e^{\beta  t}|v_t|
&\!\!\!\leq\!\!\!& \|\sqrt{t\varrho} e^{\beta  t}\nabla^2v\|_2\|v\|_6^2\|\sqrt{t\varrho} e^{\beta  t}v_t\|_6\nonumber\\
&\!\!\!\leq\!\!\!& C\|\sqrt{t} e^{\beta  t}\nabla^2v\|_2\|\nabla v\|_2^2\|\sqrt{t} e^{\beta  t}\nabla v_t\|_2\nonumber\\
&\!\!\!\leq\!\!\!&  \frac{1}{15}\|\sqrt{t}e^{\beta  t}\nabla v_t\|_2^2 + C_0e^{-4\beta_2t}
\|\sqrt{t} e^{\beta  t}\nabla^2v\|_2^2,
\end{eqnarray}
and
\begin{eqnarray} \label{I23}
I_{23} & \!\!\!\leq\!\!\!&
\|\sqrt{t\varrho} e^{\beta  t}\nabla v_t\|_2\|v\|_6^2\|\sqrt{t\varrho} e^{\beta  t}\nabla v\|_6\nonumber\\
& \!\!\!\leq\!\!\!& C\|\sqrt{t} e^{\beta  t}\nabla v_t\|_2\|\nabla v\|_2^2\|\sqrt{t} e^{\beta  t}\nabla^2 v\|_2\nonumber\\& \!\!\!\leq\!\!\!&
 \frac{1}{15}\|\sqrt{t}e^{\beta  t}\nabla v_t\|_2^2 + C_0e^{-4\beta_2t}
\|\sqrt{t} e^{\beta  t}\nabla^2v\|_2^2.
\end{eqnarray}
Hence, putting  \eqref{I21}, \eqref{I22} and \eqref{I23} together, we obtain for some $c>0,$
\begin{equation} \label{I2}
I_2 \leq  \frac{1}{5} \|\sqrt{t}e^{\beta  t}\nabla v_t\|_2^2
+ C_0e^{-ct}\|e^{\beta  t}\nabla^2 v\|_2^2\bigl(1+\|\sqrt{t\varrho}e^{\beta  t}v_t\|_2^2\bigr)\cdotp
\end{equation}
Next, we estimate $I_3$ as follows:
\begin{eqnarray} \label{I3}
I_3 &\!\!\!\leq\!\!\!& \|\nabla v\|_2 \|\sqrt{t\varrho}e^{\beta  t}v_t\|_4^2 \nonumber\\[1ex]
&\!\!\!\leq\!\!\!&
C \|\nabla v\|_2 \|\sqrt{t\varrho}e^{\beta  t}v_t\|_2^{1/2}
\|\sqrt{t}e^{\beta  t}v_t\|_6^{3/2}\nonumber\\[1ex]
&\!\!\!\leq\!\!\!& C \|\nabla v\|_2 \|\sqrt{t\varrho}e^{\beta  t}v_t\|_2^{1/2}\|\sqrt{t}e^{\beta  t}\nabla v_t\|_2^{3/2}\nonumber\\
&\!\!\!\leq\!\!\!& \frac{1}{10} \|\sqrt{t}e^{\beta  t}\nabla v_t\|_2^2 + C \|\sqrt{t\varrho}e^{\beta  t}v_t\|_2^2 \|\nabla v\|_2^4\nonumber\\
&\!\!\!\leq\!\!\!& \frac{1}{10} \|\sqrt{t}e^{\beta  t}\nabla v_t\|_2^2 + C_0e^{-4\beta_2t} \|\sqrt{t\varrho}e^{\beta  t}v_t\|_2^2
\end{eqnarray}
and
\begin{equation} \label{I4}
I_4 \leq 2\int_{\Omega}\varrho|v||\sqrt{t}e^{\beta  t}\nabla v_t||\sqrt{t}e^{\beta  t}v_t|\,{\rm d}x\leq \frac{1}{10} \|\sqrt{t}e^{\beta  t}\nabla v_t\|_2^2 + C\|v\|_{\infty}^2 \|\sqrt{t\varrho}e^{\beta  t}v_t\|_2^2.
\end{equation}
Finally, as \eqref{Poin} also applies to $v_t,$  one may write for sufficiently  small $\beta,$
\begin{equation}\label{I5} I_5 \leq \frac{1}{10} \|\sqrt{t}e^{\beta  t}\nabla v_t\|_2^2.
\end{equation}
Combining \eqref{2.7},\eqref{I2}, \eqref{I3}, \eqref{I4} and \eqref{I5}  we arrive at
\begin{multline} \label{2.30}
\frac{1}{2}\frac{d}{d t}\int_{\Omega}t e^{2\beta  t}\varrho |v_t|^2\,{\rm d}x+\frac{1}{2} \int_{\Omega}te^{2\beta  t}|\nabla v_t|^2\,{\rm d}x \\
\leq  C_0\bigl(e^{-ct}\|e^{\beta  t}\nabla^2 v\|_2^2+\|v\|_\infty^2\bigr)\|\sqrt{t\varrho} e^{\beta  t} v_t\|_2^2 \\
+C_0\bigl(\|\sqrt\varrho e^{\beta  t}v_t\|_2^2+e^{-ct}\|e^{\beta  t}\nabla^2 v\|_2^2\bigr)\cdotp
\end{multline}

By virtue of Lemma \ref{l:decay1}, the last line is integrable on $\mathbb{R}_+$ for any $\beta \leq \beta_2$, as well as the prefactor of the second line (observe that $H^2\hookrightarrow L_\infty$).
Hence, Gronwall inequality ensures that
\begin{equation} \label{2.31}
\sup_{t\in \mathbb{R}_+}\int_{\Omega} t e^{2\beta  t} \varrho|v_t|^2\,{\rm d}x + \int_{\mathbb{R}_+} e^{2\beta  t}\|\sqrt{t}\nabla v_t\|_2^2\,{\rm d}t <\infty.
\end{equation}
Now, Poincar\'e inequality implies the bound for $\|e^{\beta t}\sqrt{t}v_t\|_{L_2(\mathbb{R}_+;L_2)}$, which completes the proof.

\end{proof}

\subsection{Shift of integrability and control of  $\|\nabla v\|_{\infty}$}

Using the decay estimates we proved so far will finally enable us to establish
similar properties for higher order norms:
\begin{lem}\label{l:decay3}
Let $(\varrho,v)$ be a solution to \eqref{sys} given either by Theorem \ref{thm:d=2} or by Theorem \ref{thm:d=3}. Then, the following properties hold true:
\begin{align}
  \label{2.10e}
&\sqrt{t}e^{\beta_2 t}(\nabla^2 v,\nabla P) \in L_p(\mathbb{R}_+;L_s(\Omega)) \quad 2\leq s \leq 6\quad\hbox{and}\quad p=\frac{4s}{3s-6}\quad\hbox{if }\ d=3,\\[3pt]
 \label{2.10ee}
&\sqrt{t}e^{\beta_2 t}(\nabla^2 v,\nabla P) \in L_p(\mathbb{R}_+;L_s(\Omega)) \quad 2\leq s <\infty\quad\hbox{and}\quad p=\frac{2s}{s-2}\quad\hbox{if }\ d=2,\\[3pt]
\label{2.10f}
&e^{\beta_4  t} \nabla^2 v \in L_1(\mathbb{R}_+;L_r(\Omega)) \quad \textrm{for some} \quad r>d,
\\[3pt]
\label{2.10}
&e^{\beta_4  t}\nabla v \in L_1(\mathbb{R}_+;L_\infty(\Omega)),
\end{align}
for some $0<\beta_4=\frac{c\mu}{\varrho^*\delta^2}<\beta_2$, where $\beta_2$ is from Lemma \ref{l:decay1}.
\end{lem}
\begin{proof} We multiply \eqref{sys}$_1$ by $\sqrt{t}e^{\beta t}$, where $\beta=\beta_2$, and rewrite it as a Stokes system:
\begin{align} \label{Stokes}
-\Delta \sqrt{t}e^{\beta t}v + \nabla \sqrt{t}e^{\beta t}P = -\sqrt{t}e^{\beta t} \varrho v_t - \sqrt{t}e^{\beta t} \varrho v \cdot \nabla v, \quad
\div \sqrt{t}e^{\beta t} v=0.
\end{align}
We start with proving \eqref{2.10e}.
By the interpolation inequality
$$
\|f\|_q \leq \|f\|_2^{(6-q)/2q}\|f\|_6^{(3q-6)/2q}, \quad 2\leq q \leq 6,
$$
it is enough to prove
\begin{equation} \label{2.10e1}
\sqrt{t}e^{\beta  t}(\nabla^2 v,\nabla P) \in L_\infty(\mathbb{R}_+;L_2(\Omega)) \cap L_2(\mathbb{R}_+;L_6(\Omega)).
\end{equation}
By \eqref{decay:3} and Sobolev embedding, we have
\begin{equation} \label{2.10b}
\sqrt{\varrho t}e^{\beta  t} v_t \in L_\infty(\mathbb{R}_+;L_2(\Omega)) \cap L_2(\mathbb{R}_+;L_6(\Omega)).
\end{equation}
Therefore, the elliptic regularity of \eqref{Stokes} implies
$$\begin{aligned}
    \|\sqrt{t} e^{\beta  t}(\nabla^2 v,\nabla P)\|_{L_\infty(\mathbb{R}_+;L_2)} &\leq
    C_0\!+\! C\|\varrho v \cdot \nabla \sqrt{t} e^{\beta  t}
    v\|_{L_\infty(\mathbb{R}_+;L_2)}\\
    \leq& C_0 \!+\!C \|e^{\beta t}v\|_{L_\infty(\mathbb{R}_+;L_6)}
    \|\sqrt{t} \nabla v\|_{L^\infty(\mathbb{R}_+;L_3)}\\
    \leq& C_0\!+\! C\sqrt t \,e^{-d \beta t/6}
    \|e^{\beta  t}v\|_{L_\infty(\mathbb{R}_+;L_6)}
    \|\nabla v\|_{L_\infty(\mathbb{R}_+;L_2)}^{1-d/6}
    \|\sqrt{t}e^{\beta  t} \nabla^2 v\|_{L_\infty(\mathbb{R}_+;L_2)}^{d/6}.\end{aligned}$$
So, by Young inequality and \eqref{decay:2c}, we get
\begin{equation} \label{2.11}
     \|\sqrt{t} e^{\beta  t}(\nabla^2 v,\nabla P)\|_{L_\infty(\mathbb{R}_+;L_2)}\leq C_0.
\end{equation}
Similarly, starting from \eqref{Stokes} and thanks to Inequality \eqref{decay:3} and embedding
$H^1_\Omega\hookrightarrow L_6,$ we have
$$\begin{aligned}
    \|\sqrt{t}e^{\beta  t}(\nabla^2 v,\nabla P)\|_{L_2(\mathbb{R}_+;L_6)} &\leq C_0+C
    \|\sqrt{t}e^{\beta  t} v \cdot \nabla v\|_{L_2(\mathbb{R}_+;L_6)} \\
    &\leq C_0 + C \|e^{\beta  t/2 } \sqrt{t} v\|_{L_\infty(\mathbb{R}_+;L_\infty)}
    \|e^{\beta  t/2} \nabla v\|_{L_2(\mathbb{R}_+;L_6)}.
\end{aligned}$$
In  three-dimensional case, we have
$$
    \|\sqrt{t} e^{\beta  t/2} v\|_{L_\infty(\mathbb{R}_+;L_\infty)} \leq
    C\|\sqrt{t} e^{\beta  t/2} \nabla v\|_{L_\infty(\mathbb{R}_+;L_2)}^{1/2}\|\sqrt{t} e^{\beta  t/2} \nabla v\|_{L_\infty(\mathbb{R}_+;L_6)}^{1/2},
$$
which implies
$$
\begin{aligned}
\|\sqrt{t}e^{\beta  t} \nabla^2 v\|_{L_2(\mathbb{R}_+;L_6)} &\leq
C_0 + C\|e^{\beta t/2}\sqrt{t}\nabla v\|_{L_\infty(\mathbb{R}_+;L_2)}^{1/2}
\|e^{\beta t/2}\nabla v\|_{L_2(\mathbb{R}_+;L_6)}\|\sqrt{t}e^{\beta t/2}\nabla v\|_{L_\infty(\mathbb{R}_+;L_6)}^{1/2}\\
& \leq C_0 + C_0\|e^{\beta t}\nabla v\|_{L_\infty(\mathbb{R}_+;L_2)}^{1/2}\|\sqrt{t}e^{\beta t/2}\nabla v\|_{L_\infty(\mathbb{R}_+;L_6)}^{1/2} \\
& \leq C_0 + C_0\|e^{\beta t}\nabla v\|_{L_\infty(\mathbb{R}_+;L_2)}
+\|\sqrt{t}e^{\beta t}\nabla^2 v\|_{L_\infty(\mathbb{R}_+;L_2)},
\end{aligned}
$$
where in the first passage we used \eqref{decay:2a} and Sobolev embedding to estimate $\|e^{\beta t/2}\nabla v\|_{L_2(\mathbb{R}_+;L_6)}$.
\smallbreak
Therefore, by \eqref{decay:2} and \eqref{2.11}, choosing $\beta$ small enough, we obtain
\begin{equation*}
\|\sqrt{t}e^{\beta  t}(\nabla^2 v,\nabla P)\|_{L_2(\mathbb{R}_+;L_6)}\leq C_0.
\end{equation*}
This  completes the proof of \eqref{2.10e1}, and thus also of \eqref{2.10e}.
The easier two-dimensional case (that is \eqref{2.10ee}) is left to the reader.
\medbreak
In order to show the next part of the lemma,
we note that for $\beta_4,\delta>0$ such that $\beta_4+\delta<\beta_2$ we have %
$$
\begin{aligned}
\int_1^{+\infty} e^{\beta_4 t}\|\nabla^2 v\|_{L_s}
&\leq \left(\int_1^{+\infty} e^{-\delta p't}\,{\rm d}t\right)^{1/p'}
\left(\int_1^{+\infty} e^{(\beta_4+\delta)pt}\|\nabla^2 v(t)\|_{L_s}^p\,{\rm d}t\right)^{1/p}\\[3pt]
&\leq C \|e^{\beta_2 t}\sqrt{t} \nabla^2 v\|_{L_p(\mathbb{R}_+;L_s)}.
\end{aligned}$$
For small times we can write
$$
\int_0^1\|\nabla^2 v(t)\|_r \,{\rm d}t \leq
\left(\int_0^1 t^{-\alpha p'}\right)^{1/p'}\int_0^1 t^{\alpha p}\|\nabla^2 v(t)\|_r^p \,{\rm d}t.
$$
We can choose for instance $p=\frac{8}{5}$, which corresponds to $s=4$ in \eqref{2.10e}, then $\alpha=\frac{5}{16}$ so that $\alpha p=\frac{1}{2}$ and $\alpha p'<1$, so the first integral is again finite. This completes the proof of \eqref{2.10f}.
\smallbreak
As for  \eqref{2.10}, it results directly from \eqref{decay:2} and \eqref{2.10f} owing to a suitable Gagliardo-Nirenberg inequality that yields:
\begin{equation*}
\|e^{\beta_4 t}\nabla v\|_{L_1(\mathbb{R}_+;L_\infty)}\leq C ( \|e^{\beta_4 t}\nabla v\|_{L_1(\mathbb{R}_+;L_2)}+\|e^{\beta_4 t}\nabla^2 v\|_{L_1(\mathbb{R}_+;L_r)}.
\end{equation*}
The right-hand side is finite whenever $\beta_4<\beta_2$.
\end{proof}

\subsection* {Acknowledgments:}
The first two authors  have  been  partially supported by the ANR project INFAMIE (ANR-15-CE40-0011). The second and third author were partially supported by National Science Centre grant
No2018/29/B/ST1/00339 (Opus).