\begin{document}

 \title{Finite-Energy Solutions \  for  Compressible Two-Fluid Stokes System } 

\abstract
We prove existence of global in time weak solutions to a compressible two-fluid Stokes system  with a single velocity field and algebraic closure for the pressure law. The constitutive relation involves  densities of both fluids through  an implicit function. The system appears to be outside the class of problems that can be treated using the classical Lions-Feireisl approach. Adapting the novel compactness tool developed by the first author and P.--E. {\sc Jabin} in the mono-fluid compressible Navier-Stokes setting, we first prove the weak sequential stability of solutions. Next, we construct weak solutions via unconventional approximation using  the Lagrangian formulation, truncations and stability  result of trajectories for rough velocity fields.

\nablaoindent {\bf Keywords.} Compressible Stokes system, two-fluid model,  global weak solutions,  Lagrangian coordinates, stability of trajectories.

\section{Introduction}
Multi-component fluid models arise in various applications including studies of water wave impact on  coastal structures (violent aerated flows) \cite{DiDuGh}, petroleum industry \cite{Evje08}--\cite{Evje16}, cancer cell migration in compressible media \cite{Ev}, or turbulent mixing in nuclear industry, reactive flows, propulsion and sprays \cite{De}, to name only few. Classical derivation of multi-fluid models begins writing the equilibrium equations for each component of the flow at the microscopic level  with small scale interfaces. 
The second step is to  perform a volume averaging under suitable  closure assumptions. Averaged models  bypass the local geometrical complexity of the interphase at the cost of including new variables --
the volumetric rate of presence of each fluid/phase, see \cite{De}. This formal derivation can be found in the monographs of {\sc M. Ishii} and {\sc T. Hibiki} \cite{IsHi},  and of {\sc D. Drew} and {\sc S. L. Passman} \cite{DrPa}. 
Mathematically rigorous derivation of several models from mono-fluid systems  may be found in \cite{BrHi}, \cite{BrHi2}, \cite{BrHu}, \cite{AmZl}, \cite{Se}.  The  reader is also referred to the recent chapter \cite{BrDeGhGrHi} for discussion on modelling and mathematical studies of multi-fluid systems in the compressible setting. 

 

    In the present paper, we analyze a bi-fluid compressible system  in the semi-stationary Stokes regime. 
    We assume a common velocity field and pressure for both fluids (algebraic pressure closure). Our system of equations reads:
\begin{subequations}\label{spec}
\begin{align}
&\partialtb{\alpha^\partialm\varrho^\partialm}+\partial_{xx}iv(\alpha^{\partialm}\varrho^\partialm\vc{u})=0,\label{spec_1}\\
&-\mu\Delta\vc{u} - (\lambda+\mu) \nablaabla \partial_{xx}iv \vc{u} +\nabla p=\vc{0},\\
&\alpha^++\alpha^-=1,\\
& p = p^+ = p^- ,
\varepsilonnd{align}
\varepsilonnd{subequations}
with constant shear and bulk viscosities $\mu$ and $\lambda$ such that  $\lambda+2\mu>0$ {and $\mu>0$}. The unknowns of the system \varepsilonqref{spec} are  the volumetric rates of presence of fluid $+$ {and} $-$,  $\alpha^+,\alpha^-$, respectively,  with  
\varepsilonq{\label{alpha}
0\le \alpha^\partialm \le 1,
}
the two mass densities ${\varrho^+}$, ${\varrho^-}$, and the common velocity field $\vc{u}$. By $p^+$, $p^-$ we denote the internal barotropic pressures for each fluid with the explicit form:
\varepsilonq{\label{CP}
p^+=a^+{\varrho^+}^{\int_{\mathbb{R}}amma_+}, q_{l}uad p^-=q_{l}uad a^-{\varrho^-}^{\int_{\mathbb{R}}amma-}, 
}
where $a^\partialm>0$, $\int_{\mathbb{R}}amma^\partialm>1$ are given constants. 
The purpose of this paper is to prove the existence of solutions "\`a la {\sc Leray}" (finite energy) to this system on the $d$-dimensional torus $\mathbb{T}^d$, $d=2,3$,  under the following constraint 
\varepsilonq{\label{av}
\intO{ \vc{u}(t,x)}=0 q_{l}uad \mbox{ for } \  t\in (0,T),}
and with the initial conditions
\varepsilonq{\label{Ic_1}
\alpha^+ \varrho^+ {\cal E}rt_{t=0} = R_0, q_{l}uad \alpha^- \varrho^- {\cal E}rt_{t=0} = Q_0, q_{l}quad 
R_0 \int_{\mathbb{R}}e 0, q_{l}uad Q_0 \int_{\mathbb{R}}e 0.\\
}
Moreover, we ask for the following compatibility condition for the initial data:
$$\alpha^+ {\cal E}rt_{t=0} = \alpha^+_0,q_{l}quad \alpha^- {\cal E}rt_{t=0} = \alpha^-_0$$
\varepsilonq{\label{Ic_2}
\alpha^+_0 + \alpha^-_0 = 1, q_{l}quad \alpha^\partialm_0 \int_{\mathbb{R}}e 0, q_{l}quad
p^+(\varrho^+_0) = p^-(\varrho^-_0).
}
\begin{rmk}\label{Rmk_1}
If $\varrho_0^+$ is nonzero, then $\alpha_0^+ = R_0/\varrho_0^+$. Because of the algebraic pressure law closure which provides $\varrho_0^-= (a_+/a_-)^{1/\int_{\mathbb{R}}amma_-} (\varrho_0^+)^{{\frac{\int_{\mathbb{R}}amma^+}{\int_{\mathbb{R}}amma^-}}} \nablaot = 0$ we can also define  $\alpha_0^- = Q_0/\varrho_0^-$.
If, on the other hand, $\varrho_0^+ = \varrho_0^- = 0$, we choose $\alpha_0^\partialm  = 1/2$ for each phase.  
\varepsilonnd{rmk}
 
 Our main result for system \varepsilonqref{spec} reads:
 \begin{thm}. \label{Mainbifluid}
Let  $\int_{\mathbb{R}}amma^\partialm>1$, {$\int_{\mathbb{R}}amma^+\nablaeq \int_{\mathbb{R}}amma^-$,} $a^\partialm>0$, $\lambda + 2\mu>0$, {$\mu>0$}, and let the initial data  \varepsilonqref{Ic_1} satisfy \varepsilonqref{Ic_2}. 
Assume that 
$$\intO{ (R_0^{\int_{\mathbb{R}}amma^+} + Q_0^{\int_{\mathbb{R}}amma^-}) }< \infty, q_{l}quad
   0< \intO{ R_0 } <  \infty, q_{l}quad
   0 < \intO{ Q_0 } < \infty.$$ 
Then there exists a global weak solution $(\alpha^\partialm, \varrho^\partialm,\vc{u})$ of system \varepsilonqref{spec}--\varepsilonqref{av} 
satisfying
$$\alpha^\partialm\varrho^\partialm \in  L^\infty(0,T;L^{\int_{\mathbb{R}}amma^\partialm}(\mathbb{T}^d)) \cap L^{2\int_{\mathbb{R}}amma^\partialm}((0,T)\times \mathbb{T}^d) \cap {\cal C}([0,T]; L^{\int_{\mathbb{R}}amma^\partialm}(\mathbb{T}^d)),
$$
$$\varrho^\partialm \in  L^\infty(0,T;L^{\int_{\mathbb{R}}amma^\partialm}(\mathbb{T}^d)) \cap L^{2\int_{\mathbb{R}}amma^\partialm}((0,T)\times \mathbb{T}^d) \cap {\cal C}([0,T]; L^{\int_{\mathbb{R}}amma^\partialm}(\mathbb{T}^d)),
$$
$$ \vc{u} \in L^2(0,T; H^1(\mathbb{T}^d)),$$
where Equations \varepsilonqref{spec} are satisfied in ${\cal D}'((0,T)\times \mathbb{T}^d)$ and the initial conditions \varepsilonqref{Ic_1}
are satisfied in ${\cal D}'(\mathbb{T}^d)$. Moreover, the equations \varepsilonqref{spec_1} are satisfied in the renormalized sense.
\varepsilonnd{thm}

   Note  that if $a^+=a^-=1$ and $\int_{\mathbb{R}}amma^+ = \int_{\mathbb{R}}amma^-$, system \varepsilonqref{S} is reduced to the 
 semi-stationary Stokes version of  \varepsilonqref{entropytransport2}. In this case our technique provides stronger results that the classical technique of Lions. Indeed comparing \cite[Theorem 2]{MaMiMuNoPoZa} (see also \cite[Theorem 3.1]{Michalek}) with the proof of Theorem \varepsilonqref{Thm:main} we see that we are able to prove compactness of sequences approximating $R,Q$ and $Z$, while the results from \cite{MaMiMuNoPoZa} provide strong convergence of sequence approximating $Z$, but only a weak convergence of $\varrho_n\to\varrho$.
 To get Theorem \ref{Mainbifluid}, we first prove global existence of weak solutions for a reformulation of the system \varepsilonqref{spec}.  Introducing the notation
$$R =\varrho^+\alpha^+,q_{l}quad Q =\varrho^-\alpha^-, q_{l}quad Z= \varrho^+,$$
we check that the pressure $p$ is expressed in terms of $R,Q$. In fact we have
\begin{equation}\label{pZ}
 p=P(R,Q)= a^+ {Z}^{\int_{\mathbb{R}}amma_+},
\varepsilonnd{equation}
for $Z=Z(R,Q)$ such that 
\varepsilonq{\label{TZ}
\lr{\frac{a^-}{a^+}}^{1/\int_{\mathbb{R}}amma^-} Q = \lr{1-\frac{R}{Z} } Z^\int_{\mathbb{R}}amma,q_{l}uad  \mbox{with}q_{l}uad \int_{\mathbb{R}}amma=\frac{\int_{\mathbb{R}}amma_+}{\int_{\mathbb{R}}amma_-},
}
and
\varepsilonq{\label{RleqZ}R\leq Z.}
The system \varepsilonqref{spec}-\varepsilonqref{av}  can be therefore transformed to the following form
\begin{subequations}\label{S}
\begin{align}
  &\partialt R+ \partial_{xx}iv(R\vc{u} )=0,\label{SR}\\ 
  &\partialt Q + \partial_{xx}iv(Q\vc{u} )=0 ,\label{ST}\\%& \mbox{ in } & \mathbb{T}^d \times (0,T), \\
&- (\lambda + 2 \mu)  \partial_{xx}iv \vc{u} + a^+\lr{Z(R,Q)^{\int_{\mathbb{R}}amma^+}  - \{ Z(R,Q)^{\int_{\mathbb{R}}amma^+} \} }=0, \label{Div}\\
&{\rm rot} \, \vc{u} = 0, q_{l}uad \intO{ \vc{u}(t,x) }= 0, \label{rotu}
\varepsilonnd{align}
\varepsilonnd{subequations} 
where $\{f \}= \bigl( \int_{\mathbb{T}^d} f(x) \, dx\bigr)/|\mathbb{T}^d|$ and $Z$ is related to $R$ and $Q$ through the non-explicit formula \varepsilonqref{TZ}.
The initial conditions for the previous system reads
\varepsilonq{\label{IcNew}
R {\cal E}rt_{t=0} = R_0, q_{l}uad Q {\cal E}rt_{t=0} = Q_0,q_{l}quad  R_0 \int_{\mathbb{R}}e 0, q_{l}uad Q_0 \int_{\mathbb{R}}e 0,
}
with the following compatibility condition on  $ Z{\cal E}rt_{t=0} = Z_0$
\varepsilonq{\label{RelationInit}
\lr{\frac{a^-}{a^+}}^{1/\int_{\mathbb{R}}amma^-} Q_0 = \lr{1-\frac{R_0}{Z_0} } Z_0^\int_{\mathbb{R}}amma, 
q_{l}quad  \hookrightarrowbox{ with } q_{l}quad R_0 \le Z_0.
}
For system \varepsilonqref{S} we prove the following:
\begin{thm}\label{Thm:main}
Let $\int_{\mathbb{R}}amma^\partialm>1$, {$\int_{\mathbb{R}}amma^+\nablaeq \int_{\mathbb{R}}amma^-$}, $a^\partialm>0$, $\lambda + 2\mu>0$, {$\mu>0$}, and let the initial data be given by \varepsilonqref{IcNew} with $Z_0$
defined through \varepsilonqref{RelationInit}. Further, assume that
\varepsilonq{\label{eq_ini}
\intO{ (R_0^{\int_{\mathbb{R}}amma^+} + Q_0^{\int_{\mathbb{R}}amma^-}) }< \infty,  q_{l}quad
   0<  \intO {R_0} < \infty, q_{l}quad  0< \intO {Q_0} < \infty.
}
Then there exists $(R,Q,Z,\vc{u})$ -- a global in time weak solution to the system \varepsilonqref{TZ}--\varepsilonqref{S}  for $(t,x)\in (0,T)\times \mathbb{T}^d$,
 such that
$$R  \in L^\infty(0,T;L^{\int_{\mathbb{R}}amma^+}(\mathbb{T}^d)) \cap L^{2\int_{\mathbb{R}}amma^+}((0,T)\times \mathbb{T}^d) \cap {\cal C}([0,T]; L^{\int_{\mathbb{R}}amma^+}(\mathbb{T}^d)),
$$
$$
   Q \in L^\infty(0,T;L^{\int_{\mathbb{R}}amma^-}(\mathbb{T}^d)) \cap L^{2\int_{\mathbb{R}}amma^-}((0,T)\times \mathbb{T}^d) \cap {\cal C}([0,T]; L^{\int_{\mathbb{R}}amma^-}(\mathbb{T}^d))
$$
$$Z  \in L^\infty(0,T;L^{\int_{\mathbb{R}}amma^+}(\mathbb{T}^d)) \cap L^{2\int_{\mathbb{R}}amma^+}((0,T)\times \mathbb{T}^d) \cap {\cal C}([0,T]; L^{\int_{\mathbb{R}}amma^+}(\mathbb{T}^d)),
$$
$$\vc{u} \in L^2(0,T;H^1(\mathbb{T}^d)),
$$
where Equations \varepsilonqref{pZ}-\varepsilonqref{rotu} are satisfied in ${{\cal D}'((0,T)\times \mathbb{T}^d})$, and the initial conditions \varepsilonqref{IcNew}  with 
the constraint \varepsilonqref{RelationInit} are satisfied in  ${\cal D}'(\mathbb{T}^d)$. Moreover, equations \varepsilonqref{SR} and \varepsilonqref{ST} are satisfied in the renormalized sense.
\varepsilonnd{thm} 
 
Our paper provides the first proof of existence of global-in-time weak solutions for a bi-fluid system with constant viscosities under algebraic pressure law closure in
physical dimensions 2 and 3. The only other available results treat the density dependent viscosities with two velocity fields, see \cite{BrDeGhGr}, and \cite{BrHuLi} (see also \cite{Evje16}--\cite{Evje16a} for specific
linear pressure laws). The main difficulty in analysis of the system 
\varepsilonqref{spec} is due to, roughly speaking, complex form of the pressure. Indeed, by $\alpha^++\alpha^-=1$
 and \varepsilonqref{CP}, the pressure $p$ depends on  $\alpha^+\varrho^+$ and $\alpha^-\varrho^-$ ($R$ and $Q$, respectively) in  nonlinear implicit, see the relation \varepsilonqref{TZ}. It causes that the nowadays standard approach developed by {\sc P.-L. Lions} \cite{Lions2} and {\sc E. Feireisl} \cite{EF2001}, see also \cite{MaMiMuNoPoZa},  seems to be not applicable in all generality. 
 Therefore, we adapt a brand new technique  from \cite{BJ,BJ_short, BrJa} that could be used
 to cover more complicated case in a future work  for instance with viscosities depending
on the volumic rates $\alpha^\partialm$.   
  By reformulating system \varepsilonqref{spec} in terms of the pressure argument $Z$ and one of the conserved  quantities, for example $R=\alpha^+\varrho^+$, we show that the new technique provides compactness for sequences approximating both unknowns. 

Let us briefly discuss various contributions related to global weak solutions to the study of compressible fluid equations with intricate pressure law.

\nablaoindent {\bf Mono-fluid systems.} Note that putting $\alpha^+ =1$ in \varepsilonqref{spec}, we get the usual semi-stationary compressible Stokes system
\varepsilonq{
&\partialt \varrho +\partial_{xx}iv(\varrho\vc{u})=0,\\
&-\mu\Delta\vc{u} - (\lambda+\mu) \nablaabla \partial_{xx}iv \vc{u} +\nabla p(\varrho)=\vc{0},
\label{Quasi-Stokes}}
that has been studied for instance in \cite{Lions2} with a monotone pressure law $p(\varrho) = a \varrho^\int_{\mathbb{R}}amma$ with $\int_{\mathbb{R}}amma>1$, and
more recently in \cite{BJ_short} with non-monotone, locally Lipschitz, pressure law $p(\varrho)$, such that $p(0)=0$
and 
$$C^{-1} \varrho^\int_{\mathbb{R}}amma - C \le p(\varrho) \le C\varrho^\int_{\mathbb{R}}amma + C, q_{l}quad
   |p'(s)| \le \overline p s^{\int_{\mathbb{R}}amma-1},
$$
for some constants $C>0$, $\overline p >0$ and $\int_{\mathbb{R}}amma>1$.

\nablaoindent {\bf Compressible systems with two continuity equations.}  Two-components compressible systems have been studied in the density dependent viscous case in \cite{BrDeGhGr}, \cite{BrHuLi} for instance. 
Concerning the constant viscosity case, 
the existence of weak solutions to the two--phase model 
{\varepsilonq{
&\partialt {\varrho^+} + \partial_{xx}iv({\varrho^+}\vc{u})=0\\
&\partialt{\varrho^-} + \partial_{xx}iv({\varrho^-} \vc{u}) = 0\\
& \partialt ([{\varrho^+}+{\varrho^-}]\vc{u})+ \partial_{xx}iv ([{\varrho^+}+{\varrho^-}] \vc{u}\otimes \vc{u}) + \nablaabla p({\varrho^+},{\varrho^-}) - \mu \partial_{xx}elta \vc{u} - (\lambda + \mu) \nablaabla\partial_{xx}iv \vc{u} = \vc{0}
\label{two-phaseAV}}}
was recently proven by A. {\sc Vasseur}  {\it et al.} in \cite{VaWeYu}, for the pressure law equal to {$p({\varrho^+},{\varrho^-}) = \lr{{\varrho^+}}^\int_{\mathbb{R}}amma + \lr{{\varrho^-}}^\alpha$} with $\int_{\mathbb{R}}amma>9/5$ and $\alpha \int_{\mathbb{R}}e 1$. For the existence of strong solutions close 
to equilibrium we refer, for example, to \cite{Hao12}. For the extensive analysis of two-component models in the one-dimensional setting we refer to papers of {\sc S. Evje}  {\it et al.} \cite{Evje08, Evje15} in the constant viscosity case
and with specific linear pressure law.

In another recent paper of D. {\sc Maltese}  {\it et al.} \cite{MaMiMuNoPoZa}, see also \cite{Michalek, FeKlNoZa}, the authors considered the system used in the geophysical flow modelling
\varepsilonq{
&\partialt \varrho + \partial_{xx}iv(\varrho \vc{u})=0\\
&\partialt s  + \vc{u}\overlineot\nabla s  = 0\\
& \partialt (\varrho \vc{u}) + \partial_{xx}iv (\varrho  \vc{u}\otimes \vc{u})+ \nablaabla p(\varrho,s) - \mu \partial_{xx}elta \vc{u} - (\lambda + \mu) \nablaabla\partial_{xx}iv \vc{u} = \vc{0}
\label{entropytransport},}
where $s$ denotes the entropic variable, with the pressure law given by $p(\varrho,s) = \varrho^\int_{\mathbb{R}}amma {\cal T}(s)$ with $\int_{\mathbb{R}}amma>1$, $s>0$,
and ${\cal T}(\overlineot)$ a given smooth and strictly monotone function. They proved the existence of weak solution to the following reformulation of \varepsilonqref{entropytransport}  
\varepsilonq{\label{entropytransport2}
&\partialt \varrho + \partial_{xx}iv(\varrho \vc{u})=0\\
&\partialt Z  + \partial_{xx}iv(Z \vc{u}) = 0\\
& \partialt (\varrho \vc{u}) + \partial_{xx}iv (\varrho  \vc{u}\otimes \vc{u})+ \nablaabla Z^\int_{\mathbb{R}}amma - \mu \partial_{xx}elta \vc{u} - (\lambda + \mu) \nablaabla\partial_{xx}iv \vc{u} = \vc{0},}
where $Z$ denotes $\varrho[{\cal T}(s)]^{1/\int_{\mathbb{R}}amma}$, after which they proved the equivalence between solutions to systems \varepsilonqref{entropytransport} and \varepsilonqref{entropytransport2} for  $\int_{\mathbb{R}}amma\int_{\mathbb{R}}eq \frac95$. For analysis and numerical simulations of system \varepsilonqref{entropytransport2} with the so-called congestion constraint we refer to \cite{DeMiZa, DeMiNaZa}.
 
 

Note that the two systems mentioned before \varepsilonqref{two-phaseAV} and \varepsilonqref{entropytransport} include pressure laws which are monotone with respect to variables satisfying continuity equations. This allows the authors to adapt the tools already developed by {\sc P.--L. Lions} \cite{Lions1, Lions2} and {\sc E.~Feireisl} \cite{EF2001} for mono-fluid systems.  It is not the case for the two-fluid Stokes system \varepsilonqref{S}. For this system, although $\partialartial_RZ,\ \partialartial_Q Z\int_{\mathbb{R}}eq 0$, we do not have the crucial property
$$\Ov{Z(R,Q)}\int_{\mathbb{R}}eq Z(R,Q),$$
for $\Ov{Z(R,Q)}$ denoting the weak limit of  $Z_n(R_n,Q_n)$ when the approximation parameter $n$ goes to $+\infty$.
Looking at the reformulation \varepsilonqref{S}, the pressure argument $Z$ is given implicitly in terms of $R$ and $Q$. As a consequence, $Z$ does not satisfy a simple continuity equation, {but} we have
\varepsilonqh{
\partialt Z+\partial_{xx}iv(Z\vc{u})+\frac{(1-\int_{\mathbb{R}}amma) (Z-R)Z}{ \int_{\mathbb{R}}amma(Z-R)+R}\partial_{xx}iv\vc{u}=0.}
The additional friction term $\frac{(1-\int_{\mathbb{R}}amma) (Z-R)Z}{ \int_{\mathbb{R}}amma(Z-R)+R}\partial_{xx}iv\vc{u}$ causes that 
its seems that we have no compensation of compactness between $\partial_{xx}iv\vc{u}$ and the pressure $p$. Continuity equation with production term in the Navier-Stokes type of system has been recently investigated by {\sc N. Vauchelet} and {\sc E. Zatorska} \cite{VauchZator}.
In this case, the so-called effective viscous flux equality does not imply a strong convergence of the sequence approximating $Z$. 

In this work we  show that the recent development proposed by {\sc D. Bresch} and {\sc P.--E. Jabin} \cite{BJ} may be adapted to treat the bi-fluid system \varepsilonqref{S}. Our  work therefore provides a  generalization  of this result to the pressure law that depends on two transported quantities, as in semi-stationary compressible Stokes system. 

In \cite{BJ}, the authors explain how to handle the non-monotone truncated pressure in the heat-conducting Navier-Stokes system. In this system the pressure  depends on two variables: the density and the temperature. The density satisfies the continuity equation, while the temperature is given by the heat equation, and hence some properties providing  compactness in space  of the second unknown are available.  Our result {in this paper} covers the pressure laws depending on two quantities without knowing any {\it a-priori} compactness in space for any of them. 

 
 

In our proof we rescale the unknowns  and the  viscosity  coefficient $\lambda + 2 \mu$ so that  $a^+=a^-=1$. Since we keep $\int_{\mathbb{R}}amma^+\nablaeq \int_{\mathbb{R}}amma^-$ this assumption does not lead  to loss of generality. For the sake of brevity, we will always consider $\int_{\mathbb{R}}amma^+\le \int_{\mathbb{R}}amma^-$, equivalently $\int_{\mathbb{R}}amma\leq 1$. However,  due to the symmetry of the problem, the result will remain in force also if $\int_{\mathbb{R}}amma^+>\int_{\mathbb{R}}amma^-$. 

The paper is divided into two parts:

\nablaoindent {\bf Part I.} In Section \ref{Sect:2}, we first prove energy estimates and extra integrability properties on the solutions of the system \varepsilonqref{S}. We also study the nonlinear relation
between $Q$, $R$ and $Z$ and present the equation satisfied by $Z$. Then, in Section \ref{Sect:3},  we prove the weak sequential stability of solutions to \varepsilonqref{S}. This means that the hypothetical sequence of sufficiently smooth solutions
$\{R_n,Q_n,Z_n,\vc{u}_n\}_{n=1}^\infty$ satisfying  the energy and extra integrability estimates uniformly w.r.t. $n$, has a limit when $n\to\infty$,
that is a weak solution to \varepsilonqref{S}.

   
\nablaoindent {\bf Part II.} In Section \ref{Sect:4}, we construct the approximate solutions and show that they converge to solutions of system \varepsilonqref{S}.  The starting point is the Lagrangian reformulation of \varepsilonqref{S} with truncation of the pressure. Using  {\sc G. Crippa} and {\sc C. De Lellis}' stability of the flow  result (see \cite{CrDe} and \cite{CoCrSp}) we show that approximate solutions constructed in
the Lagrangian coordinates define suitable approximate solutions of the system in the Eulerian coordinates. These solutions satisfy the uniform bounds requested  in the first part of the paper. It is worth to emphasize that our construction does not introduce any parabolic regularization of the continuity equation commonly used in compressible setting. In a sense it is similar to construction of regular solutions \cite{Mucha2003, MuchaZaj2002, MuchaZaj2006}.

\section{Preliminary observations}\label{Sect:2}
Here we provide basic a-priori estimates for the sequence of solutions $\{R_n, Q_n, Z_n, \vc{u}_n\}_{n=1}^{\infty}$, uniformly with respect to $n$. We assume that for any $n\int_{\mathbb{R}}eq 1$ $(R_n, Q_n, Z_n, \vc{u}_n)$ is a
smooth solution to \varepsilonqref{S}, defined on $(0,T)\times\mathbb{T}^d$. We drop the index $n$, when no confusion can arise, and we recall that we assume $\int_{\mathbb{R}}amma\leq 1$.
Moreover, the results from this section do not depend on the value of the viscosity coefficient if only $\lambda+2\mu>0$. Therefore, without loss of generality we take
{$$\lambda+2\mu=1.$$}

\begin{lemma}\label{Lemma_posit}
 Let $R,Q,\vc{u}$ be sufficiently smooth solutions to (\ref{S}), then 
 \begin{equation}\label{RQ0}
 0\leq R,Q q_{l}uad \mbox{and}q_{l}uad R,Q\in L^\infty(0,T; L^1(\mathbb{T}^d)).
\varepsilonnd{equation}
Moreover, assuming $R$, $Q \int_{\mathbb{R}}e 0$ given,  there exists a unique $Z$ solving \varepsilonqref{TZ} and \varepsilonqref{RleqZ}.
\varepsilonnd{lemma}
\partialf Integrating  equation \varepsilonqref{SR} over $\mathbb{T}^d$ we deduce that
\begin{equation*}
\partial_{xx}t \intO{R}=0,
\varepsilonnd{equation*}
i.e.  $\intO{R_{0}(x)}=M_R$, implies $\intO{R(t,x)}=M_R$ for any $t\in[0,T]$.\\
Moreover, since $R$ is smooth and $R_0\int_{\mathbb{R}}eq0$, we have the following estimate
\begin{equation}\label{ch2:ron0}
R(\tau,x)\int_{\mathbb{R}}eq \inf_{x\in\mathbb{T}^d}R_0(x)\varepsilonxp\left(-\int_{0}^\tau\|\partial_{xx}iv\vc{u}\|_{L^\infty(R)}{\rm d}t\right),
\varepsilonnd{equation}
in particular $R\int_{\mathbb{R}}eq0$. Repeating the same procedure for $Q$ we obtain \varepsilonqref{RQ0}.

For fixed nonnegative $R$ and $Q$, we find
a candidate $Z\int_{\mathbb{R}}e R$ satisfying the equation \varepsilonqref{TZ}, i.e.
$$f_{R,Q}(Z) = Z^\int_{\mathbb{R}}amma - R Z^{\int_{\mathbb{R}}amma-1} - Q = 0.$$
Note that $f_{R,Q}(R) = - Q \le 0$ and  $\partialartial_Z f_{R,Q}(Z)= \int_{\mathbb{R}}amma Z^{\int_{\mathbb{R}}amma-2}(Z-R) + R Z^{\int_{\mathbb{R}}amma-2} $. 
Therefore in the range $Z\int_{\mathbb{R}}eq R$ there exists a unique $Z$ solving \varepsilonqref{TZ}.  $\Box$

Our next goal is to derive estimates for $R,\,Q,\, Z$, and $\vc{u}$ following from boundedness of the energy associated with system \varepsilonqref{S}.
\begin{lemma} \label{Lemma_est}
 Let $R,Q,\vc{u}$ be sufficiently smooth solutions to (\ref{S}), then the following estimates are valid
 \begin{equation}\label{est_ener}
  \sup_{t<T} \intO{ (Z^{\int_{\mathbb{R}}amma^+}+R^{\int_{\mathbb{R}}amma^+} + Q^{\int_{\mathbb{R}}amma^-}) } +  \intTO{|\nabla\vc{u}|^2 }\leq C,
 \varepsilonnd{equation}
and 
\begin{equation}\label{est_bog}
 \intTO{ (Z^{2\int_{\mathbb{R}}amma^+}+R^{2\int_{\mathbb{R}}amma^+} + Q^{2\int_{\mathbb{R}}amma^-})} \leq C(1+T),
\varepsilonnd{equation}
for any $T>0$.
\varepsilonnd{lemma}

\partialf 
First note that from $\intO{Z_0^{\int_{\mathbb{R}}amma^+}} <  \infty$ and \varepsilonqref{RelationInit} it follows that
\varepsilonq{\label{R0Q0}
\intO{R_0^{\int_{\mathbb{R}}amma^+}} <  \infty,q_{l}uad \intO{Q_0^{\int_{\mathbb{R}}amma^-}} <  \infty.}
On the opposite, assuming that \varepsilonqref{R0Q0} holds then using the H\"older inequality and \varepsilonqref{RelationInit}  it follows that $\intO{Z_0^{\int_{\mathbb{R}}amma^+}} <  \infty$.
  We next define $\alpha$ as a solution to
  {\varepsilonqh{
  \alpha=\left\{
  \begin{array}{ll}
  \frac{R}{Z}&  \hookrightarrowbox{ if }  \  Z\nablaot = 0\\
  \frac{1}{2} & \hookrightarrowbox{otherwise}.
  \varepsilonnd{array}
  \right.
  } }
 By \varepsilonqref{RleqZ} we deduce that $0\leq \alpha\leq 1$. 
From \varepsilonqref{TZ} we  have that 
{\varepsilonq{\label{equiv}
Z^{\int_{\mathbb{R}}amma^+}=\lr{\frac{Q}{1-\alpha}}^{\int_{\mathbb{R}}amma^-} \ \hookrightarrowbox{if}\ \alpha\nablaeq1q_{l}uad \text{ or equivalently}q_{l}uad  Z^{\int_{\mathbb{R}}amma^{+}}=\lr{\frac{R}{\alpha}}^{\int_{\mathbb{R}}amma^+}\ \hookrightarrowbox{if}\ \alpha\nablaeq 0.} }
Therefore, the gradient of the pressure can be written as
\varepsilonq{
\nabla P&=\nabla Z^{\int_{\mathbb{R}}amma^+}=\alpha\nabla \lr{\frac{R}{\alpha}}^{\int_{\mathbb{R}}amma^+}+(1-\alpha)\nabla\lr{\frac{Q}{1-\alpha}}^{\int_{\mathbb{R}}amma^-}
\\
&=\frac{\int_{\mathbb{R}}amma^+}{\int_{\mathbb{R}}amma^+-1}R \nabla \lr{\frac{R}{\alpha}}^{\int_{\mathbb{R}}amma^+-1}+\frac{\int_{\mathbb{R}}amma^-}{\int_{\mathbb{R}}amma^--1}Q\nabla \lr{\frac{Q}{1-\alpha}}^{\int_{\mathbb{R}}amma^--1}.
}
{Multiplying  the last term on the l.h.s. of \varepsilonqref{Div} by $-\partial_{xx}iv\vc{u}$, and integrating by parts, we obtain}
\varepsilonqh{
&-\intO{(P-\{P\})\partial_{xx}iv\vc{u}}=\intO{\nabla P\overlineot\vc{u}}\\
&=\intO{\lr{\frac{\int_{\mathbb{R}}amma^+}{\int_{\mathbb{R}}amma^+-1}R \nabla \lr{\frac{R}{\alpha}}^{\int_{\mathbb{R}}amma^+-1}+\frac{\int_{\mathbb{R}}amma^-}{\int_{\mathbb{R}}amma^--1}Q\nabla \lr{\frac{Q}{1-\alpha}}^{\int_{\mathbb{R}}amma^--1}}\overlineot\vc{u}}\\
&=-\frac{\int_{\mathbb{R}}amma^+}{\int_{\mathbb{R}}amma^+-1}\intO{\partial_{xx}iv(R\vc{u})\lr{\frac{R}{\alpha}}^{\int_{\mathbb{R}}amma^+-1}}
-\frac{\int_{\mathbb{R}}amma^-}{\int_{\mathbb{R}}amma^--1}\intO{\partial_{xx}iv(Q\vc{u})\lr{\frac{Q}{1-\alpha}}^{\int_{\mathbb{R}}amma^--1}}\\
&=\frac{\int_{\mathbb{R}}amma^+}{\int_{\mathbb{R}}amma^+-1}\intO{\partialt R\lr{\frac{R}{\alpha}}^{\int_{\mathbb{R}}amma^+-1}}+\frac{\int_{\mathbb{R}}amma^-}{\int_{\mathbb{R}}amma^--1}\intO{\partialt Q\lr{\frac{Q}{1-\alpha}}^{\int_{\mathbb{R}}amma^--1}}\\
&=\frac{1}{\int_{\mathbb{R}}amma^+-1}\partial_{xx}t\intO{ \lr{\frac{R}{\alpha}}^{\int_{\mathbb{R}}amma^+}\alpha}+\intO{ \lr{\frac{R}{\alpha}}^{\int_{\mathbb{R}}amma^+}\partialt \alpha}\\
&q_{l}quad+\frac{1}{\int_{\mathbb{R}}amma^--1}\partial_{xx}t\intO{ \lr{\frac{Q}{1-\alpha}}^{\int_{\mathbb{R}}amma^-}(1-\alpha)}-\intO{ \lr{\frac{Q}{1-\alpha}}^{\int_{\mathbb{R}}amma^-}\partialt \alpha}\\
&= \partial_{xx}t\intOB{\frac{1}{\int_{\mathbb{R}}amma^+-1} \lr{\frac{R}{\alpha}}^{\int_{\mathbb{R}}amma^+}\alpha+\frac{1}{\int_{\mathbb{R}}amma^--1}\lr{\frac{Q}{1-\alpha}}^{\int_{\mathbb{R}}amma^-}(1-\alpha)}.
}
Multiplying the l.h.s. of the momentum equation \varepsilonqref{Div} by $-\partial_{xx}iv\vc{u}$ we therefore get
\varepsilonq{\label{ener_est}
&\sup_{t\in (0,T)}\intOB{\frac{1}{\int_{\mathbb{R}}amma^+-1}\lr{\frac{R}{\alpha}}^{\int_{\mathbb{R}}amma^+}\alpha+ \frac{1}{\int_{\mathbb{R}}amma^--1}\lr{\frac{Q}{1-\alpha}}^{\int_{\mathbb{R}}amma^-}(1-\alpha)}\\
&q_{l}quad+ \intTO{ |\partial_{xx}iv\vc{u}|^2}\\
& \leq \intOB{\frac{1}{\int_{\mathbb{R}}amma^+-1}\lr{\frac{R_0}{\alpha_0}}^{\int_{\mathbb{R}}amma^+}\alpha_0+\frac{1}{\int_{\mathbb{R}}amma^--1} \lr{\frac{Q_0}{1-\alpha_0}}^{\int_{\mathbb{R}}amma^-}(1-\alpha_0)}.
}
On account of \varepsilonqref{equiv} and \varepsilonqref{R0Q0} the r.h.s. is bounded. Using \varepsilonqref{ener_est}, and \varepsilonqref{equiv} again, we obtain
\varepsilonqh{
\sup_{t\in (0,T)}\intO{Z^{\int_{\mathbb{R}}amma^+}}= {\sup_{t\in (0,T)}}\intOB{Z^{\int_{\mathbb{R}}amma^+}\alpha+ Z^{\int_{\mathbb{R}}amma^+}(1-\alpha)}\\=
\sup_{t\in (0,T)}\intOB{Z^{\int_{\mathbb{R}}amma^+}\alpha+ \lr{\frac{Q}{1-\alpha}}^{\int_{\mathbb{R}}amma^-}(1-\alpha)}\leq C,
}
which provides the uniform estimate for $Z^{\int_{\mathbb{R}}amma^+}$, as stated in \varepsilonqref{est_ener}.
The uniform estimate for $R^{\int_{\mathbb{R}}amma^+}$ can be obtained using \varepsilonqref{RleqZ}, and then \varepsilonqref{equiv} provides the uniform estimate for $Q^{\int_{\mathbb{R}}amma^-}$. 

In order to estimate the full gradient of $\vc{u}$ we notice that
\varepsilonqh{\intTO{|\nabla\vc{u}|^2}=\intTO{|\partial_{xx}iv\vc{u}|^2}+\intTO{|{\rm rot}\,\vc{u}|^2},
}
the first term is bounded due to \varepsilonqref{ener_est}, and the second one is equal to $0$ on account of \varepsilonqref{rotu}.

In order to prove \varepsilonqref{est_bog} we multiply momentum equation \varepsilonqref{Div} by $P$ and integrating over time and space we get that
\varepsilonq{\label{P2}
\intTO{P^2}&=\intTO{\partial_{xx}iv\vc{u}\, P}+\intTO{\{P\} P}\\
&=\intTO{\partial_{xx}iv\vc{u}\, P}+\intT{\{P\}^2}.
}
The last term is bounded due to \varepsilonqref{est_ener}, and we use the Cauchy inequality to estimate
\varepsilonq{\label{divuP}
\left|\intTO{\partial_{xx}iv\vc{u} P}\right|&\leq \frac12\intTO{P^2}+\frac12\intTO{|\partial_{xx}iv\vc{u}|^2}\\
&\leq C+ \frac12\intTO{P^2}.
}
So, combining \varepsilonqref{divuP} with \varepsilonqref{P2} we obtain the uniform estimate for $Z^{2\int_{\mathbb{R}}amma^+}$ as in \varepsilonqref{est_bog}. The rest of bounds from \varepsilonqref{est_bog} follows, as previously, from the relations between $Z$, $R$, and $Q$, see  \varepsilonqref{RleqZ}, and \varepsilonqref{equiv}. $\Box$

We can now use the above estimates in order to deduce that $R$ and $Q$ satisfy equations \varepsilonqref{SR} and \varepsilonqref{ST} in the renormalized sense.
\begin{lemma}\label{Lemma_ren}
Assume $X \in L^q((0,T)\times\mathbb{T}^d)$ with $q\int_{\mathbb{R}}eq 2$, and $\vc{u} \in L^2(0,T;W^{1,2}(\mathbb{T}^d))$. Let $(X,\vc{u})$ solve 
$$\partialt X+\partial_{xx}iv(X\vc{u})=f~\text{in}~{\cal D}'((0,T)\times \mathbb{T}^d),$$
where $f\in L^p((0,T)\times\mathbb{T}^d)$ for some $p>1$, $p'\lr{\frac{q}{2}-1}\leq 1$.
Then $(X,\vc{u})$ is also a renormalized solution i.e. it solves 
\begin{equation}\label{renorent}
\partialartial_t b(X) + \partial_{xx}iv(b(X)\vc{u}) + \big(b'(X)X-b(X)\big) \partial_{xx}iv\vc{u} = f {b'(X)}~\text{in}~{\cal D}'((0,T)\times \mathbb{T}^d),
\varepsilonnd{equation}
where
\begin{equation}\label{regb}
b \in C^1([0,\infty)),\ | b'(s)|\leq Cs^\lambda,\ for\ s>1,\ where \ \lambda\leq \frac{q}{2}-1.
\varepsilonnd{equation}
\varepsilonnd{lemma}
The proof of this lemma is a consequence of the DiPerna-Lions theory \cite{DL} of renormalized solutions to the transport equation.

We now derive the equation satisfied by $Z$. If $R,Q,\vc{u}$ are smooth, the evolution equation for $Z=Z(R,Q)$ can be deduced from the continuity equations for $R$ and for $Q$, and the formula \varepsilonqref{TZ}. However, we will use the the equivalency between the equations \varepsilonqref{SR}, \varepsilonqref{ST} and the evolution equation for $Z$ at the level of weak solutions, i.e. solutions with regularity specified in Lemma \ref{Lemma_est}. For such solutions we have the following result
\begin{lemma}\label{Lemma_equiv}
Let $\vc{u}\in L^2(0,T; W^{1,2}(\mathbb{T}^d))$, $R\in L^{2\int_{\mathbb{R}}amma^+}((0,T)\times \mathbb{T}^d)\cap L^\infty(0,T; L^{\int_{\mathbb{R}}amma^+}(\mathbb{T}^d))$, 
$Q\in L^{2\int_{\mathbb{R}}amma^-}((0,T)\times \mathbb{T}^d)\cap L^\infty(0,T; L^{\int_{\mathbb{R}}amma^-}(\mathbb{T}^d))$,
and let $(R, Q,\vc{u})$ solve \varepsilonqref{SR} and \varepsilonqref{ST} in the sense of distributions, and {\varepsilonqref{Div} and \varepsilonqref{rotu} a.e. in $(0,T)\times\Omega$}. Then $Z$ defined by \varepsilonqref{TZ} belongs to $L^{2\int_{\mathbb{R}}amma^+}((0,T)\times \mathbb{T}^d)\cap L^\infty(0,T; L^{\int_{\mathbb{R}}amma^+}(\mathbb{T}^d))$ and it satisfies 
\varepsilonq{\label{SZ}
\partialt Z+\partial_{xx}iv(Z\vc{u})+\frac{(1-\int_{\mathbb{R}}amma) (Z-R)Z}{ \int_{\mathbb{R}}amma(Z-R)+R}\partial_{xx}iv\vc{u}=0,}
in the sense of distributions.
{Conversly}, let $(R,Z,\vc{u})$ solve \varepsilonqref{SR} and \varepsilonqref{SZ} in the sense of distributions. Then $Q$ defined by \varepsilonqref{TZ} satisfies \varepsilonqref{ST} in the sense of distributions. 
\varepsilonnd{lemma}
\partialf The fact that $Z\in L^{2\int_{\mathbb{R}}amma^+}((0,T)\times \mathbb{T}^d)\cap L^\infty(0,T; L^{\int_{\mathbb{R}}amma^+}(\mathbb{T}^d))$ follows from Lemma \ref{Lemma_est}. Testing the equations \varepsilonqref{SR} and \varepsilonqref{ST} by ${\rm x}i_\varepsilonta(x-\overlineot)$, where ${\rm x}i_\varepsilonta$ is a standard periodized mollifier, we obtain
\varepsilonq{\label{RQeta}
\partialt R_\varepsilonta+\partial_{xx}iv(R_\varepsilonta\vc{u})=r^1_\varepsilonta,\\
\partialt Q_\varepsilonta+\partial_{xx}iv(Q_\varepsilonta\vc{u})=r^2_\varepsilonta,
}
satisfied a.e. in $(0,T)\times\mathbb{T}^d$, where $a_\varepsilonta$ denotes $a\ast{\rm x}i_\varepsilonta$. From the Friedrichs commutator lemma we know that 
\varepsilonq{\label{r1r2}
r^1_\varepsilonta\to 0 q_{l}uad\mbox{in} \ L^{p_1}((0,T)\times \mathbb{T}^d),\\
r^2_\varepsilonta\to 0 q_{l}uad\mbox{in} \ L^{p_2}((0,T)\times \mathbb{T}^d)}
for $\frac{1}{p_1}=\frac{1}{2}+\frac{1}{2\int_{\mathbb{R}}amma^+}$, $\frac{1}{p_2}=\frac{1}{2}+\frac{1}{2\int_{\mathbb{R}}amma^-}$.
We now define $Z_\varepsilonta$ via
$$Q_\varepsilonta=\lr{1-\frac{R_\varepsilonta}{Z_\varepsilonta}}Z_\varepsilonta^\int_{\mathbb{R}}amma,$$
and as previously we find that $R_\varepsilonta\leq Z_\varepsilonta$. Let us now apply $\partialartial_{Q_\varepsilonta}$, $\partialartial_{R_\varepsilonta}$ to both sides of the above formula,  we obtain respectively:
\[1=\int_{\mathbb{R}}amma Z_\varepsilonta^{\int_{\mathbb{R}}amma-1}\partialartial_{Q_\varepsilonta} Z_\varepsilonta-R_\varepsilonta(\int_{\mathbb{R}}amma-1)Z_\varepsilonta^{\int_{\mathbb{R}}amma-2}\partialartial_{Q_\varepsilonta} Z_\varepsilonta,\]
\[0=\int_{\mathbb{R}}amma Z_\varepsilonta^{\int_{\mathbb{R}}amma-1}\partialartial_{R_\varepsilonta} Z_\varepsilonta-Z_\varepsilonta^{\int_{\mathbb{R}}amma-1}-R_\varepsilonta(\int_{\mathbb{R}}amma-2)Z_\varepsilonta^{\int_{\mathbb{R}}amma-1}\partialartial_{R_\varepsilonta} Z_\varepsilonta,
\]
therefore,
\varepsilonq{\label{pTR}
\partialartial_{Q_\varepsilonta} Z_\varepsilonta=\frac{1}{\int_{\mathbb{R}}amma Z_\varepsilonta^{\int_{\mathbb{R}}amma-1}-R_\varepsilonta(\int_{\mathbb{R}}amma-1)Z_\varepsilonta^{\int_{\mathbb{R}}amma-2}},q_{l}uad
\partialartial_{R_\varepsilonta} Z_\varepsilonta= \frac{Z_\varepsilonta^{\int_{\mathbb{R}}amma-1}}{\int_{\mathbb{R}}amma Z_\varepsilonta^{\int_{\mathbb{R}}amma-1}-R_\varepsilonta(\int_{\mathbb{R}}amma-1)Z_\varepsilonta^{\int_{\mathbb{R}}amma-2}}.
}
Using the assumption $\int_{\mathbb{R}}amma\leq1$  and inequality $R_\varepsilonta\leq Z_\varepsilonta$, we check that we can estimate the partial derivatives of $Z_\varepsilonta(R_\varepsilonta,Q_\varepsilonta)$ by the integrable function, more precisely
\varepsilonqh{
|\partialartial_{Q_\varepsilonta} Z_\varepsilonta|\leq\frac{Z_\varepsilonta^{1-\int_{\mathbb{R}}amma}}{\int_{\mathbb{R}}amma},q_{l}uad |\partialartial_{R_\varepsilonta} Z_\varepsilonta|\leq \frac{1}{\int_{\mathbb{R}}amma}.
}
It means that $\partialartial_{Q_\varepsilonta} Z_\varepsilonta$ and $\partialartial_{R_\varepsilonta} Z_\varepsilonta$ are suitable test functions for  \varepsilonqref{RQeta}.
We now check that
\varepsilonqh{
\partialt Z_\varepsilonta&=\partialartial_{Q_\varepsilonta}Z_\varepsilonta\partialt Q_\varepsilonta+\partialartial_{R_\varepsilonta} Z_\varepsilonta\partialt R_\varepsilonta\\
&= -\partialartial_{Q_\varepsilonta}Z_\varepsilonta\partial_{xx}iv(Q_\varepsilonta\vc{u})-\partialartial_{R_\varepsilonta} Z_\varepsilonta\partial_{xx}iv(R_\varepsilonta\vc{u})+r^1_\varepsilonta\partialartial_{R_\varepsilonta} Z_\varepsilonta+r^2_\varepsilonta\partialartial_{Q_\varepsilonta} Z_\varepsilonta\\
&= -[\partialartial_{Q_\varepsilonta}Z_\varepsilonta\nabla Q_\varepsilonta+\partialartial_{R_\varepsilonta} Z_\varepsilonta\nabla R_\varepsilonta]\overlineot\vc{u}-[Q_\varepsilonta\partialartial_{Q_\varepsilonta} Z_\varepsilonta +R_\varepsilonta\partialartial_{R_\varepsilonta} Z_\varepsilonta ]\partial_{xx}iv\vc{u}\\
&q_{l}uad+r^1_\varepsilonta\partialartial_{R_\varepsilonta} Z_\varepsilonta+r^2_\varepsilonta\partialartial_{Q_\varepsilonta} Z_\varepsilonta\\
&=-\nabla Z_\varepsilonta\overlineot\vc{u}-[Q_\varepsilonta\partialartial_{Q_\varepsilonta} Z_\varepsilonta +R_\varepsilonta\partialartial_{R_\varepsilonta} Z_\varepsilonta ]\partial_{xx}iv\vc{u}+r^1_\varepsilonta\partialartial_{R_\varepsilonta} Z_\varepsilonta+r^2_\varepsilonta\partialartial_{Q_\varepsilonta} Z_\varepsilonta,
}
therefore
\varepsilonq{\label{ZA}
\partialt Z_\varepsilonta+\partial_{xx}iv(Z_\varepsilonta\vc{u})+[Q_\varepsilonta\partialartial_{Q_\varepsilonta} Z_\varepsilonta +R_\varepsilonta\partialartial_{R_\varepsilonta} Z_\varepsilonta -Z_\varepsilonta]\partial_{xx}iv\vc{u}=r^1_\varepsilonta\partialartial_{R_\varepsilonta} Z_\varepsilonta+r^2_\varepsilonta\partialartial_{Q_\varepsilonta} Z_\varepsilonta.
}
Substituting \varepsilonqref{pTR} we compute
\varepsilonqh{
 &Q_\varepsilonta\partialartial_{Q_\varepsilonta} Z_\varepsilonta +R_\varepsilonta\partialartial_{R_\varepsilonta} Z_\varepsilonta-Z_\varepsilonta
 = \frac{Q_\varepsilonta+R_\varepsilonta Z_\varepsilonta^{\int_{\mathbb{R}}amma -1}-Z_\varepsilonta(\int_{\mathbb{R}}amma Z_\varepsilonta^{\int_{\mathbb{R}}amma-1} - R_\varepsilonta(\int_{\mathbb{R}}amma -1) Z_\varepsilonta^{\int_{\mathbb{R}}amma-1})}{\int_{\mathbb{R}}amma Z_\varepsilonta^{\int_{\mathbb{R}}amma-1} - R_\varepsilonta(\int_{\mathbb{R}}amma -1) Z_\varepsilonta^{\int_{\mathbb{R}}amma-1}} \\
&= \frac{Z_\varepsilonta^\int_{\mathbb{R}}amma - R_\varepsilonta Z_\varepsilonta^{\int_{\mathbb{R}}amma-1} + R_\varepsilonta Z_\varepsilonta^{\int_{\mathbb{R}}amma-1} - \int_{\mathbb{R}}amma Z_\varepsilonta^\int_{\mathbb{R}}amma + R_\varepsilonta(\int_{\mathbb{R}}amma -1) Z_\varepsilonta^{\int_{\mathbb{R}}amma -1}}{\int_{\mathbb{R}}amma Z_\varepsilonta^{\int_{\mathbb{R}}amma-1} - R_\varepsilonta(\int_{\mathbb{R}}amma -1) Z_\varepsilonta^{\int_{\mathbb{R}}amma-1}} =
 \frac{ (1-\int_{\mathbb{R}}amma)(Z_\varepsilonta-R_\varepsilonta)Z_\varepsilonta}{\int_{\mathbb{R}}amma(Z_\varepsilonta-R_\varepsilonta)+R_\varepsilonta}.
}
Using this, the equation for $Z_\varepsilonta$ can be written as
\varepsilonq{\label{Zs}
\partialt Z_\varepsilonta+\partial_{xx}iv(Z_\varepsilonta\vc{u})+\frac{(1-\int_{\mathbb{R}}amma) (Z_\varepsilonta-R_\varepsilonta)Z_\varepsilonta}{ \int_{\mathbb{R}}amma(Z_\varepsilonta-R_\varepsilonta)+R_\varepsilonta}\partial_{xx}iv\vc{u}=r^1_\varepsilonta\partialartial_{R_\varepsilonta} Z_\varepsilonta+r^2_\varepsilonta\partialartial_{Q_\varepsilonta} Z_\varepsilonta.
}
Note that since $\int_{\mathbb{R}}amma <1$, we easily show that $\frac{(1-\int_{\mathbb{R}}amma) (Z_\varepsilonta-R_\varepsilonta)Z_\varepsilonta}{ \int_{\mathbb{R}}amma(Z_\varepsilonta-R_\varepsilonta)+R_\varepsilonta}\int_{\mathbb{R}}eq 0$, moreover $\frac{(1-\int_{\mathbb{R}}amma) (Z_\varepsilonta-R_\varepsilonta)Z_\varepsilonta}{ \int_{\mathbb{R}}amma(Z_\varepsilonta-R_\varepsilonta)+R_\varepsilonta}\leq \frac{1-\int_{\mathbb{R}}amma}{ \int_{\mathbb{R}}amma} (Z_\varepsilonta-R_\varepsilonta)\in L^{2\int_{\mathbb{R}}amma^+}((0,T)\times\mathbb{T}^d)$. This allows us to let $\varepsilonta\to 0$ in the l.h.s. of \varepsilonqref{Zs}. The r.h.s. of \varepsilonqref{Zs} vanishes provided that $p_1\int_{\mathbb{R}}eq 1$ and $\frac{1}{p_2}+\frac{1-\int_{\mathbb{R}}amma}{2\int_{\mathbb{R}}amma^+}\leq 1$, which is fulfilled provided that $\int_{\mathbb{R}}amma^+\int_{\mathbb{R}}eq1$.  

In order to recover  the equation for $Q$ from the equations for $R$ and $Z$, one derives the equation on $Q^{1/\int_{\mathbb{R}}amma}$ first. It is easy to observe that \varepsilonqref{TZ} yields
$$
\left|\frac{\partialartial Q^{1/\int_{\mathbb{R}}amma}}{\partialartial Z} \right|\leq 1 q_{l}uad\mbox{ and } q_{l}uad\left|\frac{\partialartial Q^{1/\int_{\mathbb{R}}amma}}{\partialartial R} \right|\leq 1/\int_{\mathbb{R}}amma.
$$
Therefore, the rigorous procedure involving mollifying  and passage to the limit $\varepsilonta\to 0$, described above, can be repeated for $Q$. In this manner we obtain a renormalized version of the equation for  $Q$ 
$$
\partialt Q^{1/\int_{\mathbb{R}}amma} + \partial_{xx}iv (Q^{1/\int_{\mathbb{R}}amma} \vc{u}) + Q^{1/\int_{\mathbb{R}}amma}\partial_{xx}iv \vc{u}=0.
$$
This finishes the proof. $\Box$

As a consequence of this Lemma and Lemma \ref{Lemma_ren} we have:
\begin{cor}
The couple $(Z,\vc{u})$ is a renormalized solution to \varepsilonqref{SZ}.
\varepsilonnd{cor}
\partialf Indeed, for the assumptions of Lemma \ref{Lemma_ren} to be fulfilled we note that
$$\left|\frac{(1-\int_{\mathbb{R}}amma) (Z-R)Z}{ \int_{\mathbb{R}}amma(Z-R)+R}\partial_{xx}iv\vc{u}\right|\leq C Z|\partial_{xx}iv\vc{u}|,$$
and the r.h.s. is bounded in $L^1((0,T)\times\mathbb{T}^d)$ on account of \varepsilonqref{est_ener} and \varepsilonqref{est_bog}. $\Box$

\section{Sequential stability of solutions}\label{Sect:3}
The purpose of this section is to pass to the limit $n\to\infty$ in the sequence  $\{R_n, Q_n, Z_n, \vc{u}_n\}_{n=1}^{\infty}$ and  to verify that the limit $(R,Q,Z,\vc{u})$ satisfies the system \varepsilonqref{S} in the weak sense. We  prove the following theorem
\begin{thm}\label{Thm_main2}
Let $T>0$. Assume that for any $n$ the quadruple 
 $(R_n, Q_n, Z_n,\vc{u}_n)$ satisfies \varepsilonqref{TZ}--\varepsilonqref{S} with the initial conditions 
 \varepsilonqh{
 R_n {\cal E}rt_{t=0} = R_{0,n}, q_{l}uad Q_n {\cal E}rt_{t=0} = Q_{0,n}, q_{l}uad R_{0,n} \int_{\mathbb{R}}e 0, q_{l}uad Q_{0,n} \int_{\mathbb{R}}e 0,
 }
 with $Z_n {\cal E}rt_{t=0} = Z_{0,n}$ satisfying \varepsilonqref{RelationInit},  and s.t.
 \varepsilonqh{
 R_{n,0}\to R_0q_{l}uad strongly \ in\ L^1(\mathbb{T}^d),\\
 Z_{n,0}\to Z_0q_{l}uad strongly \ in\ L^1(\mathbb{T}^d).
 }
Let the estimates from the Lemmas \ref{Lemma_posit} and \ref{Lemma_est} hold uniformly with respect to $n$. Then up to the subsequence
\varepsilonqh{
&R_n\to Rq_{l}uad strongly\ in \ L^{2\int_{\mathbb{R}}amma^+-\varepsilonp}((0,T)\times\mathbb{T}^d),\\
&Q_n\to Qq_{l}uad strongly \ in \  L^{2\int_{\mathbb{R}}amma^--\varepsilonp}((0,T)\times\mathbb{T}^d),\\
&Z_n\to Zq_{l}uad strongly\ in \ L^{2\int_{\mathbb{R}}amma^+-\varepsilonp}((0,T)\times\mathbb{T}^d),\\
&\vc{u}_n\to \vc{u}q_{l}uad weakly\ in\  L^2(0,T; H^1(\mathbb{T}^d)),
}
for any $\varepsilonp>0$. 
Moreover, $(R,Q,Z,\vc{u})$ satisfies \varepsilonqref{TZ}--\varepsilonqref{S} in the sense of distributions. 
\varepsilonnd{thm}
Passage to the limit $n\to\infty$ in the two first equations of system \varepsilonqref{S} requires at least weak convergence of the sequences $R_n, Q_n$ and $\vc{u}_n$. This can be deduced directly from the a-priori estimates from Lemmas  \ref{Lemma_posit} and \ref{Lemma_est} using nowadays classical techniques (see, for example {\sc P.-L. Lions} \cite{Lions2} or {\sc E. Feireisl} \cite{EF2001}), and we skip this part.
The  core of the proof of Theorem \ref{Thm_main2} is to pass to the limit in the nonlinear term of the  momentum equation \varepsilonqref{Div}. Indeed, identification of the limit $\lim_{n\to\infty} p(R_n, Q_n)= p(R,Q)$ requires some sort of strong convergence of sequences $R_n$, and $Q_n$. Instead of proving the strong convergence of these sequences directly, we use  the equivalence between the system \varepsilonqref{S} and its reformulation in terms of $(R,Z,\vc{u})$, as stated in the Lemma \ref{Lemma_equiv}.
The proof of Theorem \varepsilonqref{Thm_main2} can be therefore reduced to the proof of compactness of the sequence $Z_n$ and justification that the limit quantities $R,\ Q,\ Z$ satisfy the relation \varepsilonqref{TZ}. We follow the strategy proposed by {\sc D. Bresch \& P.-E. Jabin} \cite{BJ}  (see also \cite{BJ_short}) in the context of compressible Navier-Stokes equations with the non-monotone pressure law.  As a byproduct of this approach, we obtain a compactness result for the sequence $R_n$, and using \varepsilonqref{TZ} the strong convergence of the sequence $Q_n$, concluding the proof of Theorem \ref{Thm_main2}. The rest of this section will be therefore devoted to the proof of the following result:
\begin{prop}\label{prop:rho}
Let $T>0$. Assume that $\{(R_n, Z_n,\vc{u}_n)\}_{n=1}^{\infty}$ satisfies 
\begin{subequations}\label{S_Z}
\begin{align}
 &\partialt R+ \partial_{xx}iv(\vc{u} R)=0,\label{SR_Z}\\ 
&\partialt Z+\partial_{xx}iv(Z\vc{u})+\frac{(1-\int_{\mathbb{R}}amma) (Z-R)Z}{ \int_{\mathbb{R}}amma(Z-R)+R}\partial_{xx}iv\vc{u}=0 ,\label{ST_Z}\\%& \mbox{ in } & \mathbb{T}^d \times (0,T), \\
 &{\partial_{xx}iv\vc{u} =Z^{\int_{\mathbb{R}}amma^+}- \{Z^{\int_{\mathbb{R}}amma^+}\}\label{Su_Z}},\\ 
&{\rm rot}\ \vc{u}=0,q_{l}uad \intO{ \vc{u}(t,x) }= 0,\label{Sur_Z}
\varepsilonnd{align}
\varepsilonnd{subequations}
 with the initial conditions \varepsilonqref{IcNew} satisfying \varepsilonqref{RelationInit}, and let the estimates from the Lemmas \ref{Lemma_posit} and \ref{Lemma_est} hold uniformly with respect to $n$. 
Then the sequences $\{R_n\}_{n=1}^\infty$ $\{Z_n\}_{n=1}^\infty$ are compact in $L^1((0,T)\times\mathbb{T}^d)$.
\varepsilonnd{prop}

\subsection{Preliminaries}\label{ssec:prelim}
In order to prove the strong convergence of $\{R_n,Z_n\}_{n=1}^{\infty}$ necessary to pass to the limit in the momentum equation, we will use the compactness criterion introduced in the context of Navier-Stokes equations in \cite{BJ}. 
First let us introduce the necessary notation.
We define the positive, bounded and symmetric function $\{K_h\}_{h>0}$ such that  
$$ K_h(x)=\frac{1}{(|x|+h)^{a}}$$ 
with $|x|=\sqrt {\sum_{i=1}^d x_i^2}$ for $|x|\leq 1/2$ with some $a>d$ and $K_h$ positive, independent of $h$ for $|x|\int_{\mathbb{R}}e 2/3$, $K_h$ positive constant outside $B(0,3/4)$ so that $K_h\in{\cal C}^\infty(\mathbb{T}^d\backslash B(0,3/4))$ and it is a periodic function.
Further, we denote
\varepsilonq{\label{notationK}
\overline{K_h}(x) = \frac{K_h(x)}{\|K_h\|_{L^1(\mathbb{T}^d)}},
q_{l}quad
{\mathcal K}_{h_0}(x) = \int_{h_0}^1 \overline{K_h}(x) \frac{dh}{h}.
}
We  also use the following properties of the kernel $K_h$: 
\begin{equation}\label{ineqK}
K_h(x)= K_h(-x), q_{l}quad |x||\nablaabla K_h(x)| \leq C K_h(x),
\varepsilonnd{equation}
for some constant $C>0$ independent of $h$ and
\varepsilonq{\label{Klogh}
\|{\cal K}_{h_0}\|_{L^1(\mathbb{T}^d)} \sim |\log h_0|.}

We recall the following compactness criterion, for the proof see \cite{Belgacem}, Lemma 3.1.
\begin{lemma}\label{lem:compact}
Let $\{X_n\}_{n=1}^\infty$ be a sequence of functions uniformly bounded in $L^p((0,T)\times\mathbb{T}^d)$
with $1\leq p<+\infty$. Assume that ${\cal K}_h$ is a sequence of positive, bounded functions s.t.
\begin{enumerate}
\item[i)] $\forall \varepsilonta>0$, \ $\sup_h \partialisplaystyle \intO{{\cal K}_h(x)\vc{1}_{{\{x:\> |x|\int_{\mathbb{R}}eq\varepsilonta\}}}}<\infty$,
\item[ii)] $\|{\cal K}_h\|_{L^1(\mathbb{T}^d)}\to+\infty$\ as\ $h\to 0$.
\varepsilonnd{enumerate}
If $\{\partialt X_n\}_{n=1}^\infty$ is uniformly bounded in $L^r([0,T],W^{-1,r}(\mathbb{T}^d))$ with 
$r\int_{\mathbb{R}}eq 1$ and  
$$
\underset{n}{\lim\sup} \left( \frac{1}{\|{\cal K}_h\|_{L^1}}  \int_0^T\!\!\!\iintO{
{\cal K}_h(x-y)|X_n(t,x)-X_n(t,y)|^p} \, {\rm d} t \right) \to 0, q_{l}uad \mbox{ as } h\to 0,
$$
then, $\{X_n\}_{n=1}^\infty$ is compact in $L^p([0,T]\times\mathbb{T}^d)$.
Conversely, if $\{X_n\}_{n=1}^\infty$ is compact in $L^p([0,T]\times\mathbb{T}^d)$, 
then the above $\lim\sup$ converges to $0$ as $h$ goes to $0$.
\varepsilonnd{lemma}
\subsection{Propagation of oscillations }
Having  transport equations for $R$ and $Z$ together with necessary a-priori bounds, our next goal is to derive the equations for perturbations of both of these quantities. 
Perturbations are described by the evolution of $|R(t,x)-R(t,y)|$ and $|Z(t,x)-Z(t,y)|$, respectively, for any couple of points $x,y\in\mathbb{T}^d$. To obtain them we first subtract the equations for $R(t,x)$ and $R(t,y)$  
\varepsilonqh{
&\partialt (R_x-R_y) + \partial_{xx}iv_x (\vc{u}_x\lr{R_x-R_y}) +
\partial_{xx}iv_y (\vc{u}_y\lr{R_x-R_y})  \\
&=\frac 12 (\partial_{xx}iv_x \vc{u}_x + \partial_{xx}iv_y \vc{u}_y) \lr{R_x-R_y}  - \frac 12 (\partial_{xx}iv_x \vc{u}_x-\partial_{xx}iv_y \vc{u}_y)(R_x+R_y),
}
where we denoted $R_x=R(t,x), \ R_y=R(t,y)$.
Multiplying this equation by the sign of their difference $s_R=sign(R_x-R_y)$ we get
\varepsilonq{\label{eqR}
&\partialt |R_x-R_y| + \partial_{xx}iv_x (\vc{u}_x|R_x-R_y|) +
\partial_{xx}iv_y (\vc{u}_y|R_x-R_y|)  \\
&=\frac 12 (\partial_{xx}iv_x \vc{u}_x + \partial_{xx}iv_y \vc{u}_y) |R_x-R_y|  - \frac 12 (\partial_{xx}iv_x \vc{u}_x-\partial_{xx}iv_y \vc{u}_y)(R_x+R_y)s_R.
}
By the similar token, we obtain the equation for $|Z_x-Z_y|$, namely
\varepsilonq{\label{eqZ}
&\partialt|Z_x-Z_y|+
\partial_{xx}iv_{x}(\vc{u}_x|Z_x-Z_y|)+ \partial_{xx}iv_{y}(\vc{u}_y|Z_x-Z_y|)\\
&=\frac 12 (\partial_{xx}iv_x \vc{u}_x + \partial_{xx}iv_y \vc{u}_y) |Z_x-Z_y| - \frac 12 (\partial_{xx}iv_x \vc{u}_x-\partial_{xx}iv_y \vc{u}_y)(Z_x+Z_y)s_Z\\
&q_{l}uad-\left[
\frac{(1-\int_{\mathbb{R}}amma) (Z_x-R_x)Z_x}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x} \partial_{xx}iv_x \vc{u}_x-
\frac{(1-\int_{\mathbb{R}}amma) (Z_y-R_y)Z_y}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}\partial_{xx}iv_y\vc{u}_y
\right] s_Z.
}
We  now multiply \varepsilonqref{eqR} and \varepsilonqref{eqZ} by $w_x+w_y$, where $w_x=w(t,x)$ denotes the solution to the transport equation
\varepsilonq{\label{def_w}
\left\{\begin{array}{l}
\partialt w+\vc{u}\overlineot\nabla w+ \theta {\cal D} w=0,\\
w(0,x)=1,
\varepsilonnd{array}\right.
}
where $\theta$ is a constant parameter that will be chosen later on and ${\cal D}$ depending on $\vc{u}$ and $Z$.
Our next step is to write equation for 
\varepsilonqh{
S(t):=\iintO{ K_h(x-y)O_{x-y} (w_x+w_y)},
}
where $O_{x-y}=|R_x-R_y|+ |Z_x-Z_y|$, we have
\varepsilonq{\label{Rt0}
\partial_{xx}t S(t)&= \iintO{ \nabla K_h(x-y)(\vc{u}_x-\vc{u}_y)O_{x-y} (w_x+w_y)}\\
&q_{l}uad+\frac 12 \iintO{K_h(x-y)(\partial_{xx}iv_x \vc{u}_x + \partial_{xx}iv_y \vc{u}_y)O_{x-y} (w_x+w_y)}\\
&q_{l}uad - \frac 12 \iintO{K_h(x-y)(\partial_{xx}iv_x \vc{u}_x-\partial_{xx}iv_y \vc{u}_y)[(R_x+R_y)s_R+(Z_x+Z_y)s_Z](w_x+w_y)}\\
&q_{l}uad-\iintO{K_h(x-y)\left[
\frac{(1-\int_{\mathbb{R}}amma) (Z_x-R_x)Z_x}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x} \partial_{xx}iv_x \vc{u}_x-
\frac{(1-\int_{\mathbb{R}}amma) (Z_y-R_y)Z_y}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}\partial_{xx}iv_y\vc{u}_y
\right] s_Z(w_x+w_y)}\\
&q_{l}uad+ \iintO{K_h(x-y)O_{x-y} \lr{\partialt w_x+\vc{u}_x\overlineot\nabla w_x}}\\
&q_{l}uad+ \iintO{K_h(x-y)O_{x-y} \lr{\partialt w_y+\vc{u}_y\overlineot\nabla w_y}}
}
Using the symmetry of $K_h(x-y)$, $O_{x-y}$, and the symmetry of the second, third and fourth integrals on the r.h.s. we obtain
 \varepsilonq{\label{Rt1}
&\partial_{xx}t S(t)= \iintO{ \nabla K_h(x-y)(\vc{u}_x-\vc{u}_y)O_{x-y} (w_x+w_y)}\\
&q_{l}uad+ \iintO{K_h(x-y)(\partial_{xx}iv_x \vc{u}_x + \partial_{xx}iv_y \vc{u}_y)O_{x-y} w_x}\\
&q_{l}uad - \iintO{K_h(x-y)(\partial_{xx}iv_x \vc{u}_x-\partial_{xx}iv_y \vc{u}_y)[(R_x+R_y)s_R+(Z_x+Z_y)s_Z]w_x}\\
&q_{l}uad-2\iintO{K_h(x-y)\left[
\frac{(1-\int_{\mathbb{R}}amma) (Z_x-R_x)Z_x}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x} \partial_{xx}iv_x \vc{u}_x-
\frac{(1-\int_{\mathbb{R}}amma) (Z_y-R_y)Z_y}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}\partial_{xx}iv_y\vc{u}_y
\right] s_Zw_x}\\
&q_{l}uad+2 \iintO{K_h(x-y)O_{x-y} \lr{\partialt w_x+\vc{u}_x\overlineot\nabla w_x}}.
}
Finally, writing $(\partial_{xx}iv_x \vc{u}_x + \partial_{xx}iv_y \vc{u}_y)=-(\partial_{xx}iv_x \vc{u}_x - \partial_{xx}iv_y \vc{u}_y)+2\partial_{xx}iv_x \vc{u}_x$ and combining the second and the third term on the r.h.s. we obtain
 \varepsilonq{\label{Rt2}
&\partial_{xx}t S(t)= \iintO{ \nabla K_h(x-y)(\vc{u}_x-\vc{u}_y)O_{x-y} (w_x+w_y)}\\
&q_{l}uad - \iintO{K_h(x-y)(\partial_{xx}iv_x \vc{u}_x-\partial_{xx}iv_y \vc{u}_y)[R_xs_R+Z_x s_Z]w_x}\\
&q_{l}uad-2\iintO{K_h(x-y)\left[
\frac{(1-\int_{\mathbb{R}}amma) (Z_x-R_x)Z_x}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x} \partial_{xx}iv_x \vc{u}_x-
\frac{(1-\int_{\mathbb{R}}amma) (Z_y-R_y)Z_y}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}\partial_{xx}iv_y\vc{u}_y
\right] s_Zw_x}\\
&q_{l}uad+2 \iintO{K_h(x-y)O_{x-y} \lr{\partialt w_x+\vc{u}_x\overlineot\nabla w_x+\partial_{xx}iv_x\vc{u}_x w_x}}\\
&=I_1+I_2+I_3+I_4.
}
We now estimate each term in \varepsilonqref{Rt2}, separately.

\nablaoindent{\bf Estimate of $I_1$.} Recall that $K_h$ satisfies \varepsilonqref{ineqK}, 
we also know that
\varepsilonq{\label{est_Du}
|\vc{u}_x-\vc{u}_y|\leq C|x-y|(D_{|x-y|}\vc{u}_x+D_{|x-y|}\vc{u}_y),q_{l}uad\mbox{where}q_{l}uad D_h \vc{u}_x=\frac{1}{h}\int_{|z|\leq h}\frac{|\nabla \vc{u}_{x+z}|}{|z|^{d-1}}\,{\rm d}z.}
Recall that $D_h \vc{u}  \leq M |\nablaabla \vc{u}|$, where $M$  denotes the maximal operator:
$$
M f(x) = \sup_{r\leq 1} \frac{1}{|B(0,r)|} \int_{B(0,r)} f(x+z)\,dz.
$$
For the proof of this fact we refer the reader to \cite[Lemma 3.1 and Eq. 3.3]{Jabin2010}.
Combining estimate \varepsilonqref{est_Du} with \varepsilonqref{ineqK}, we have
\varepsilonqh{
I_1&=\iintO{ \nabla K_h(x-y)(\vc{u}_x-\vc{u}_y)O_{x-y} (w_x+w_y)}\\
&\leq C\iintO{K_h(x-y)(D_{|x-y|}\vc{u}_x+D_{|x-y|}\vc{u}_y)O_{x-y} w_x}.
}
Next, writing
$$
D_{|x-y|}\vc{u}_x+D_{|x-y|}\vc{u}_y =  D_{|x-y|}\vc{u}_y-D_{|x-y|}\vc{u}_x + 2D_{|x-y|}\vc{u}_x,
$$ 
and estimating $D_{|x-y|}\vc{u}(x)$ by the Maximal operator $M|\nabla\vc{u}|(x)$ we get
\varepsilonq{\label{est_I1}
I_1\leq &C\iintO{K_h(x-y)(D_{|x-y|}\vc{u}_y-D_{|x-y|}\vc{u}_x)O_{x-y} w_x}\\
&+C\iintO{K_h(x-y)M|\nabla\vc{u}_x|O_{x-y} w_x}.
}

\nablaoindent{\bf Estimate of $I_2$.}  {Due to \varepsilonqref{Su_Z}, we can write that}
\varepsilonqh{
 -\lr{ \partial_{xx}iv_x\vc{u}_x-\partial_{xx}iv_y\vc{u}_y}[R_xs_R+Z_x s_Z]
 =- \lr{ Z_x^{\int_{\mathbb{R}}amma^+}-Z_y^{\int_{\mathbb{R}}amma^+}}R_xs_R- \lr{ Z_x^{\int_{\mathbb{R}}amma^+}-Z_y^{\int_{\mathbb{R}}amma^+}}Z_x s_Z.
 }
 Note that the last term is always nonpositive, hence the contribution to $I_1$ coming from this term can be moved to the l.h.s. of \varepsilonqref{Rt2}. Concerning the first term, it can only be estimated from above. Using $R_x=\alpha_x Z_x$ we get
 \varepsilonq{\label{est_I2}
I_2= &- \iintO{K_h(x-y)(\partial_{xx}iv_x \vc{u}_x-\partial_{xx}iv_y \vc{u}_y)[R_xs_R+Z_x s_Z]w_x}\\
\leq&- \iintO{K_h(x-y)|Z_x^{\int_{\mathbb{R}}amma^+}-Z_y^{\int_{\mathbb{R}}amma^+}|Z_xw_x}\\
&+
\iintO{K_h(x-y)|Z_x^{\int_{\mathbb{R}}amma^+}-Z_y^{\int_{\mathbb{R}}amma^+}|\alpha_xZ_xw_x}\\
=&\iintO{K_h(x-y)|Z_x^{\int_{\mathbb{R}}amma^+}-Z_y^{\int_{\mathbb{R}}amma^+}|(\alpha_x-1)Z_xw_x}.
 }
 Since $\alpha_x\leq1$, the integral $I_2$ is nonpositive and can be moved to the left hand side of \varepsilonqref{Rt2}.

\nablaoindent{\bf {Estimate of $I_3$}.}
 The estimate of this term is most lengthy and requires splitting the term in a big bracket several times. We first note that
\varepsilonq{\label{re1}
&-\left[
\frac{(1-\int_{\mathbb{R}}amma) (Z_x-R_x)Z_x}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x} \partial_{xx}iv_x \vc{u}_x-
\frac{(1-\int_{\mathbb{R}}amma) (Z_y-R_y)Z_y}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}\partial_{xx}iv_y\vc{u}_y
\right] s_Z\\
&=(\int_{\mathbb{R}}amma-1)\left[
\frac{ (Z_x-R_x)Z_x^{\int_{\mathbb{R}}amma^++1}}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}-
\frac{(Z_y-R_y)Z_y^{\int_{\mathbb{R}}amma^++1}}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}
\right] s_Z\\
&-(\int_{\mathbb{R}}amma-1)\left[ \frac{ (Z_x-R_x)Z_x}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}-
\frac{(Z_y-R_y)Z_y}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y} 
\right] s_Z \{Z^{\int_{\mathbb{R}}amma^+}\}\\
&=(\int_{\mathbb{R}}amma-1)\frac{(Z_y-R_y)}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}\lr{Z_x^{\int_{\mathbb{R}}amma^++1}-Z_y^{\int_{\mathbb{R}}amma^++1}} s_Z\\
&+(\int_{\mathbb{R}}amma-1)
\left[
\frac{ (Z_x-R_x)}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}-
\frac{(Z_y-R_y)}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}
\right] Z_x^{\int_{\mathbb{R}}amma^++1} s_Z\\
&-(\int_{\mathbb{R}}amma-1)\frac{(Z_y-R_y)}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}\lr{Z_x-Z_y} s_Z \{Z^{\int_{\mathbb{R}}amma^+}\} \\
&-(\int_{\mathbb{R}}amma-1)
\left[
\frac{ (Z_x-R_x)}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}-
\frac{(Z_y-R_y)}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}
\right] Z_x\{Z^{\int_{\mathbb{R}}amma^+} \}s_Z,
}
{where we used \varepsilonqref{Su_Z} to substitute} for $\partial_{xx}iv\vc{u}$, and the observation that $\{Z_x^{\int_{\mathbb{R}}amma^+}\}=\{Z_y^{\int_{\mathbb{R}}amma^+}\}$.

Observe that
\varepsilonq{\label{sign_frac}
\frac{(Z_y-R_y)}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}=\frac{(1-\alpha_y)}{ \int_{\mathbb{R}}amma(1-\alpha_y)+\alpha_y}\ \in  \lr{0,1/{\int_{\mathbb{R}}amma}},
}
 therefore, the first term on the r.h.s.  of \varepsilonqref{re1} is non-positive for $\int_{\mathbb{R}}amma<1$ and can be eventually considered on the l.h.s. of \varepsilonqref{Rt2}. For the second term on the r.h.s.  of \varepsilonqref{re1} we continue to write
 \varepsilonqh{
 &(\int_{\mathbb{R}}amma-1)
\left[
\frac{ (Z_x-R_x)}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}-
\frac{(Z_y-R_y)}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}
\right] Z_x^{\int_{\mathbb{R}}amma^++1} s_Z\\
&=(\int_{\mathbb{R}}amma-1)\frac{R_y Z_x-R_xZ_y}{\lr{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}\lr{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}}Z_x^{\int_{\mathbb{R}}amma^++1} s_Z\\
&=(\int_{\mathbb{R}}amma-1)\frac{R_y (Z_x-Z_y)}{\lr{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}\lr{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}}Z_x^{\int_{\mathbb{R}}amma^++1} s_Z\\
&+(\int_{\mathbb{R}}amma-1)\frac{(R_y -R_x)Z_y}{\lr{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}\lr{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}}Z_x^{\int_{\mathbb{R}}amma^++1} s_Z,
}
 again, the first term has a good sign, while the second one can be transformed to
 \varepsilonqh{
 &(\int_{\mathbb{R}}amma-1)\frac{(R_y -R_x)Z_y}{\lr{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}\lr{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}}Z_x^{\int_{\mathbb{R}}amma^++1} s_Z\\
 &=(\int_{\mathbb{R}}amma-1)\frac{(R_y -R_x)}{\lr{ \int_{\mathbb{R}}amma(1-\alpha_x)+\alpha_x}\lr{ \int_{\mathbb{R}}amma(1-\alpha_y)+\alpha_y}}Z_x^{\int_{\mathbb{R}}amma^+} s_Z\\
& \leq\frac{1-\int_{\mathbb{R}}amma}{\int_{\mathbb{R}}amma^2}|R_x-R_y| Z_x^{\int_{\mathbb{R}}amma^+}.
 }
The third and fourth terms on the r.h.s. of \varepsilonqref{re1} can be treated similarly and estimated by 
$\frac{1-\int_{\mathbb{R}}amma}{\int_{\mathbb{R}}amma}|Z_x-Z_y|\{Z^{\int_{\mathbb{R}}amma^+} \}$, and  by $\frac{1-\int_{\mathbb{R}}amma}{\int_{\mathbb{R}}amma^2}O_{x-y}\{Z^{\int_{\mathbb{R}}amma^+} \}$, respectively.

 Putting all the terms together, we conclude that
 \varepsilonq{\label{est_I3}
&I_3=-2\iintO{K_h(x-y)\left[
\frac{(1-\int_{\mathbb{R}}amma) (Z_x-R_x)Z_x}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x} \partial_{xx}iv_x \vc{u}_x-
\frac{(1-\int_{\mathbb{R}}amma) (Z_y-R_y)Z_y}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}\partial_{xx}iv_y\vc{u}_y
\right] s_Zw_x}\\
&\leq -2(1-\int_{\mathbb{R}}amma)\iintO{K_h(x-y)\frac{(Z_y-R_y)}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}\left|Z_x^{\int_{\mathbb{R}}amma^++1}-Z_y^{\int_{\mathbb{R}}amma^++1}\right|w_x}\\
&q_{l}uad-2(1-\int_{\mathbb{R}}amma)\iintO{K_h(x-y)\frac{R_y |Z_x-Z_y|}{\lr{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}\lr{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}}Z_x^{\int_{\mathbb{R}}amma^++1}w_x}\\
&q_{l}uad+C\iintO{K_h(x-y)|R_x-R_y| Z_x^{\int_{\mathbb{R}}amma^+}w_x}+C\{Z^{\int_{\mathbb{R}}amma^+} \}\iintO{K_h(x-y)O_{x-y}w_x},
 }
 with some constant $C$ depending on $\int_{\mathbb{R}}amma$.

\nablaoindent{\bf Estimate of $I_4$.}
  For the last term in \varepsilonqref{Rt2} we simply use the definition of the measure $w$ from \varepsilonqref{def_w}, we therefore get
  \varepsilonq{\label{est_I4}
  I_4=&2 \iintO{K_h(x-y)O_{x-y} \lr{\partialt w_x+\vc{u}_x\overlineot\nabla w_x+\partial_{xx}iv_x\vc{u}_x w_x}}\\
  =&2 \iintO{K_h(x-y)O_{x-y} (\partial_{xx}iv_x\vc{u}_x-\theta {\cal D}_x) w_x}.
  }
  As a conclusion, we obtain from \varepsilonqref{Rt2}, using estimates \varepsilonqref{est_I1}, \varepsilonqref{est_I2}, \varepsilonqref{est_I3}, and \varepsilonqref{est_I4} that
  \varepsilonq{\label{mainS}
  \partial_{xx}t S(t)&+\iintO{K_h(x-y)|Z_x^{\int_{\mathbb{R}}amma^+}-Z_y^{\int_{\mathbb{R}}amma^+}|(1-\alpha_x)Z_xw_x}\\
  &+2(1-\int_{\mathbb{R}}amma)\iintO{K_h(x-y)\frac{(Z_y-R_y)}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}\left|Z_x^{\int_{\mathbb{R}}amma^++1}-Z_y^{\int_{\mathbb{R}}amma^++1}\right|w_x}\\
  &+2(1-\int_{\mathbb{R}}amma)\iintO{K_h(x-y)\frac{R_y |Z_x-Z_y|}{\lr{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}\lr{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}}Z_x^{\int_{\mathbb{R}}amma^++1}w_x}\\
  \leq
  &C\iintO{K_h(x-y)(D_{|x-y|}\vc{u}_y-D_{|x-y|}\vc{u}_x)O_{x-y} w_x}\\
&+C\iintO{K_h(x-y)M|\nabla\vc{u}_x|O_{x-y} w_x}\\
&+C\iintO{K_h(x-y)|R_x-R_y| Z_x^{\int_{\mathbb{R}}amma^+}w_x}+C\{Z^{\int_{\mathbb{R}}amma^+} \}\iintO{K_h(x-y)O_{x-y}w_x}\\
&+2 \iintO{K_h(x-y)O_{x-y} (\partial_{xx}iv_x\vc{u}_x -\theta {\cal D}_x)w_x}.
  }
Note that the only assumptions are: $\int_{\mathbb{R}}amma<1$, i.e. $\int_{\mathbb{R}}amma^+<\int_{\mathbb{R}}amma^-$, and $\alpha\in[0,1]$.

\subsection{Compactness criterion with weights}
At this point of the proof it is convenient to chose the ${\cal D}$ from the definition of the weight $w_x$ \varepsilonqref{def_w}. Taking for example
\varepsilonq{\label{def_D}
{\cal D}=M|\nabla\vc{u}|+|\partial_{xx}iv\vc{u}|+Z^{\int_{\mathbb{R}}amma^+}+ \{Z^{\int_{\mathbb{R}}amma^+}\},}
and $\theta$ sufficiently large we obtain from \varepsilonqref{mainS} 
\varepsilonqh{\partial_{xx}t S(t)\leq C\iintO{K_h(x-y)(D_{|x-y|}\vc{u}_y-D_{|x-y|}\vc{u}_x)O_{x-y} w_x}.}
Note that, thanks to uniform estimates from Lemmas \ref{Lemma_posit} and \ref{Lemma_est},
we have that ${\cal D}$ defined in \varepsilonqref{def_D} is uniformly bounded in $L^2((0,T)\times\mathbb{T}^d)$ (note that $\|M|\nabla\vc{u}|\|_{L^2((0,T)\times\mathbb{T}^d)}\leq C\|\nabla\vc{u}\|_{L^2((0,T)\times\mathbb{T}^d)}$). This allows us to deduce the following properties of the weight $w$. 

\begin{prop}[Proposition 7.2 in \cite{BJ}]
\label{prop:w}
Assume that ${\cal D}$ defined in \varepsilonqref{def_D} is given and that it is  bounded in $L^2(0,T\times\mathbb{T}^d)$.
Then, there exists a weight $w$ solving  \varepsilonqref{def_w}. Moreover, we have
\begin{itemize}
\item[(i)] For any $(t,x)\in(0,T)\times T^d$, $0\leq w(t,x) \leq 1$.
\item[(ii)] If we assume moreover that the pair $(X,\vc{u})$ is a solution to the continuity equation:
\varepsilonqh{\partialt X+\partial_{xx}iv(X\vc{u})=0,}
 and $X$ is bounded in $L^2((0,T)\times\mathbb{T}^d)$, there exists $C\int_{\mathbb{R}}eq 0$, such that
\begin{equation}\label{boundlogw}
\intO{ X |\log w| } \leq C \theta.
\varepsilonnd{equation}
\varepsilonnd{itemize}
\varepsilonnd{prop}
Let us now introduce
\varepsilonq{\label{Sfini}
S_{h_0}(t)=\int_{h_0}^1 \frac{S(t)}{\|K_h\|_{L^1}}\,\frac{dh}{h}=\int_{h_0}^1\int_{\mathbb{T}^{2d}} \overline{K_h}(x-y)O_{x-y} (w_x+w_y)\, {\rm d} {x}\, {\rm d} {y}\,\frac{{\rm d}h}{h}.}
Recalling the notation \varepsilonqref{notationK}, and changing the variables $z=x-z$ we get from \varepsilonqref{Sfini} that
\varepsilonqh{\partial_{xx}t S_{h_0}(t)\leq &C \int_{h_0}^1 \iintO{\overline{K_h}(x-y) \left|D_{|x-y|}\vc{u}(y)-D_{|x-y|}\vc{u}(x)\right|O_{x-y} w_x}\,\frac{{\rm d}h}{h}\\
\leq &C \int_{h_0}^1 \int_{\mathbb{T}^{d}}\overline{K_h}(z) \|D_{|z|}\vc{u}(\overlineot)-D_{|z|}\vc{u}(\overlineot+z)\|_{L^2} \,{\rm d}z\,\frac{{\rm d}h}{h}
,}
where to get the last inequality we have used the Cauchy-Schwarz inequality, $L^2$ integrability of $R$ and $Z$ together with boundedness of the weight $w$. Integrating with respect to time we obtain
\varepsilonq{\label{Rt0}
S_{h_0}(t)-S_{h_0}(0)\leq  &C\int_0^t\int_{h_0}^1 \int_{\mathbb{T}^{d}}\overline{K_h}(z) \|D_{|z|}\vc{u}(\overlineot)-D_{|z|}\vc{u}(\overlineot+z)\|_{L^2} \,{\rm d}z\,\frac{{\rm d}h}{h}\,{\rm d}s.
}
We now use the following lemma from \cite{BJ}.
\begin{lemma}[Lemma 6.3 in \cite{BJ}]
{There exists $C>0$ such that} for any $\vc{u}\in H^1(\mathbb{T}^d)$,
\begin{equation}\label{ineq:lem63}
\int_{h_0}^1 \int_{\mathbb{T}^{d}} \overline{K_h}(z) \|D_{|z|}\vc{u} - D_{|z|}\vc{u}(\overlineot+z)\|_{L^2(\mathbb{T}^d)} \,{\rm d}z \frac{{\rm d}h}{h}
\leq C |\log h_0|^{1/2} \|\vc{u}\|_{H^1(\mathbb{T}^d)}.
\varepsilonnd{equation}
\varepsilonnd{lemma}
With this at hand \varepsilonqref{Rt0} gives
\varepsilonqh{
S_{h_0}(t)-S_{h_0}(0)\leq  C|\log(h_0)|^{1/2} \int_0^t \|\vc{u}(s)\|_{H^1(\mathbb{T}^d)}\, {\rm d}s.
}
From this, boundedness of $\vc{u}$ in $L^2(0,T; H^1(\mathbb{T}^d))$ we get that
\varepsilonq{\label{Sg}
S_{h_0}(t)\leq C|\log(h_0)|^{1/2}+ S_{h_0}(0).}
Using the reversed Lemma \ref{lem:compact} and strong convergence of initial data we get that  
\varepsilonqh{
&\limsup_n\frac{S_{h_0}(0)}{|\log h_0|}\\
&=\limsup_n\lr{\frac{\iintO{ {\mathcal K}_{h_0}(x-y)\Big[ \left|Z_{0,n}(x)-Z_{0,n}(y)\right|+ |R_{0,n}(x)-R_{0,n}(y)|\Big] }}{|\log h_0|}}\to 0,
}
as $h_0\to 0$. Therefore, \varepsilonqref{Sg} yields that
\varepsilonq{\label{Sw}
\sup_{t\in[0,T]}\frac{S_{h_0}}{|\log h_0|}\to 0,q_{l}uad \mbox{as}q_{l}uad h_0\to0.
}

\subsection{Removal of the weights}
We now want to remove the weights from \varepsilonqref{Sw} to prove that
\varepsilonq{\label{without_w}
\limsup_n\lr{\frac{1}{|\log h_0|}\iintO{ {\mathcal K}_{h_0}(x-y)\Big[ \left|(Z_x)_n-(Z_y)_n\right|+ |(R_x)_n-(R_y)_n|\Big] }}\to 0
}
as $h_0\to 0$, 
while for the moment we only have only convergence with the weights, see \varepsilonqref{Sfini} and \varepsilonqref{Sw}.
{We present a formal
argument leading to the estimate \varepsilonqref{boundlogw} for $X=Z_n$}. It is clear that estimate \varepsilonqref{boundlogw} holds for $X=R_n$, we will explain why one can expect the same estimate for $X=Z_n$.
Recall that from \varepsilonqref{def_w} know that the weights satisfy the equation
\varepsilonq{\label{logw}
\partialt|\log w_n|+\vc{u}_n\overlineot \nabla |\log w_n|=\theta {\cal D}_n,}
with
\varepsilonqh{
{\cal D}_n=M|\nabla\vc{u}_n|+|\partial_{xx}iv\vc{u}_n|+Z_n^{\int_{\mathbb{R}}amma^+}+ \{Z_n^{\int_{\mathbb{R}}amma^+}\},}
Multiplying \varepsilonqref{logw} by $R_n$ and using the fact that $R_n$ satisfies the continuity equation we get
\varepsilonq{\label{cwR}
\partial_{xx}t\intO{R_n|\log w_n|}= \theta\intO{R_n{\cal D}_n}.}
Mimicking this procedure for $Z_n$ satisfying equation \varepsilonqref{ST_Z} we have
\varepsilonq{\label{Zpre}
\partial_{xx}t\intO{Z_n|\log w_n|}=\theta\intO{Z_n{\cal D}_n}- \intO{\frac{(1-\int_{\mathbb{R}}amma) (Z_n-R_n)Z_n}{ \int_{\mathbb{R}}amma(Z_n-R_n)+R_n}\partial_{xx}iv\vc{u}_n|\log w_n|}.
}
The integrability of the terms $R_n{\cal D}_n$ and $Z_n{\cal D}_n$ on r.h.s. of {\varepsilonqref{cwR} and \varepsilonqref{Zpre}, respectively}, follows from the fact that $R_n{\cal D}_n\leq Z_n{\cal D}_n$ is bounded in $L^1((0,T)\times\mathbb{T}^d)$. Indeed, on account of Lemma \ref{Lemma_est} and the assumption $\int_{\mathbb{R}}amma^+\int_{\mathbb{R}}eq1$, we have $Z_n^{\int_{\mathbb{R}}amma^++1} \in L^{1}((0,T)\times \mathbb{T}^d)$. Therefore from \varepsilonqref{cwR} it follows that
\varepsilonq{\label{lwR}
\intO{R_n|\log w_n|}\leq C\theta.}
To deduce the same inequality for $Z_n$, we rewrite the term $\partial_{xx}iv\vc{u}_n$  in \varepsilonqref{Zpre} using the formula \varepsilonqref{Su_Z} 
\varepsilonqh{
\partial_{xx}t\intO{Z_n|\log w_n|}+\intO{\frac{(1-\int_{\mathbb{R}}amma) (Z_n-R_n)Z_n}{ \int_{\mathbb{R}}amma(Z_n-R_n)+R_n}Z_n^{\int_{\mathbb{R}}amma^+}|\log w_n|}\\
=\theta\intO{Z_n{\cal D}_n}- \intO{\frac{(1-\int_{\mathbb{R}}amma) (Z_n-R_n)Z_n}{ \int_{\mathbb{R}}amma(Z_n-R_n)+R_n}\{Z_n^{\int_{\mathbb{R}}amma^+}\}|\log w_n|}
.} 
 Recalling \varepsilonqref{sign_frac}, the second term on the l.h.s. is nonnegative, while the second term on the r.h.s. can be bounded, therefore we have
 \varepsilonq{\label{cwZ}
\partial_{xx}t\intO{Z_n|\log w_n|}
\leq \theta\intO{Z_n{\cal D}_n}+\frac{1-\int_{\mathbb{R}}amma}{\int_{\mathbb{R}}amma}\{Z_n^{\int_{\mathbb{R}}amma^+}\}\intO{Z_n|\log w_n|}
.} 
Applying the Gronwall Lemma we can show that
\varepsilonqh{
&\intO{Z_n|\log w_n|(t)}\\
&\leq \varepsilonxp\left\{\frac{1-\int_{\mathbb{R}}amma}{\int_{\mathbb{R}}amma}\int_0^t\{Z_n^{\int_{\mathbb{R}}amma^+}(s)\}\,{\rm d} s\right\}
\lr{\intO{Z_n|\log w_n|(0)}+ \theta\int_0^t\intO{Z_n{\cal D}_n(t)}\,{\rm d} s}.
}
So, due to \varepsilonqref{def_w}, Lemma \ref{Lemma_est} and the discussion above we get that
\varepsilonq{\label{lwZ}
\intO{Z_n|\log w_n|}\leq C\theta.}
The above reasoning may be made rigorous by following the proof of the Proposition \ref{prop:w}  presented in \cite[Proposition 7.2]{BJ} with minor changes due to presence of additional terms in the equation for $Z_n$.

Let $\varepsilonta<1$. We define $\omega_\varepsilonta = \{x: w \leq \varepsilonta\}$ and denote by 
$\omega_\varepsilonta^c$ its complementary.
We have 
\begin{align}
\iintO{{\mathcal K}_{h_0}(x-y)O_{x-y}} 
&= \int_{h_0}^1 \iintO{ \overline{K_h}(x-y)O_{x-y}}\frac{{\rm d}h}{h} = A_1 + A_2, \label{eqI12}
\varepsilonnd{align}
with
\begin{align*}
A_1 & = \int_{h_0}^1 \int_{\{x\in\omega_\varepsilonta^c\}\cup\{y\in\omega_\varepsilonta^c\}} \overline{K_h}(x-y) O_{x-y}{\rm d} {x}\,{\rm d} {y}\frac{{\rm d}h}{h} \leq \frac{1}{\varepsilonta} S_{h_0},
\varepsilonnd{align*}
due to \varepsilonqref{Sfini}, 
and
\begin{align*}
A_2 & =\int_{h_0}^1 \int_{\{x\in\omega_\varepsilonta\}\cap\{y\in\omega_\varepsilonta\}} \overline{K_h}(x-y) O_{x-y}\,{\rm d} {x}\,{\rm d} {y}\frac{{\rm d}h}{h}  \\
& \leq 2 \int_{h_0}^1 \iintO{\overline{K_h}(x-y) [R(x)+Z(x)]\frac{|\log w(x)|}{|\log \varepsilonta|}}\frac{{\rm d}h}{h},
\varepsilonnd{align*}
where we use the symmetry of $K_h$ and the fact that, by definition, 
for $\varepsilonta<1$, $|\log w(x)|\int_{\mathbb{R}}eq |\log \varepsilonta|$ for all $x\in\omega_\varepsilonta$.
Changing the variables $z=y-x$, and recalling $\|\overline{K_h}\|_{L^1}=1$,
we get
$$
A_2 \leq \frac{2 |\log h_0|}{|\log \varepsilonta|} \intO{[ R(x)+Z(x)]|\log w(x)| },
$$
where the integral on the r.h.s. is bounded due to \varepsilonqref{lwR} and \varepsilonqref{lwZ}.

Summarizing the estimates of $A_1$, $A_2$, we obtain
$$
\intO{ {\mathcal K}_{h_0}(x-y)\Big[ \left|(Z_x)_n-(Z_y)_n\right|+ |(R_x)_n-(R_y)_n|\Big] } 
\leq \frac{1}{\varepsilonta} S_{h_0} + C\frac{|\log h_0|}{|\log \varepsilonta|}.
$$
 Applying \varepsilonqref{Sg}, we obtain
$$
\iintO{ {\mathcal K}_{h_0}(x-y)\Big[ \left|(Z_x)_n-(Z_y)_n\right|+ |(R_x)_n-(R_y)_n|\Big] }
\leq C\frac{|\log h_0|^{1/2}+1}{\varepsilonta}  + C\frac{|\log h_0|}{|\log \varepsilonta|},
$$
Next, due to \varepsilonqref{Klogh}, we have
\varepsilonqh{
\frac{1}{\|{\mathcal K}_{h_0}\|_{L^1}}\iintO{ {\mathcal K}_{h_0}(x-y)\Big[ \left|(Z_x)_n-(Z_y)_n\right|+ |(R_x)_n-(R_y)_n|\Big]}
\\
\leq C\frac{|\log h_0|^{1/2}+1}{\varepsilonta|\log h_0|}  + \frac{C}{|\log \varepsilonta|},
}
so, choosing for example $\varepsilonta=|\log h_0|^{-1/4}$, and letting $h_0\to 0$, we show \varepsilonqref{without_w}, which proves the compactness criterion from Lemma \ref{lem:compact}.
\subsection{Concluding remarks}
From the previous section it follows that the sequences $R_n$ and $Z_n$ converge strongly in $L^1((0,T)\times\mathbb{T}^d)$ to $R,\ Z$, respectively. Using the uniform bounds from Lemma \ref{Lemma_est}, we therefore deduce the strong convergence of both of these sequences in $L^{2\int_{\mathbb{R}}amma^+-\varepsilonp}((0<T)\times\mathbb{T}^d)$, for any $\varepsilonp>0$. Using the equivalence relation \varepsilonqref{TZ}, and the uniform estimates on $Q_n$, we deduce that $Q_n\to Qq_{l}uad strongly \ in \  L^{2\int_{\mathbb{R}}amma^--\varepsilonp}((0<T)\times\mathbb{T}^d)$. With this information at hand, it is possible to pass to the limit in all equations of system  \varepsilonqref{TZ}--\varepsilonqref{S}. The proof of Theorem \ref{Thm_main2} is therefore complete. $\Box$

\section{Existence of approximate solutions}\label{Sect:4}
Here we prove the existence of solutions to certain approximation of system \varepsilonqref{S} involving parameter $k$.  Let us denote  ${\cal T}_k(\overlineot)$ the truncation operator $ {\cal T}_k: \mathbb{R}_+ \to \mathbb{R}_+$ for $k\in \mathbb{R}_+$ such that
\varepsilonqh{
 {\cal T}_k(t)= t \mbox{ for } t \le kq_{l}uad  \mbox{ and }q_{l}uad  {\cal T}_k(t)=k \mbox{ for }
    t\int_{\mathbb{R}}eq k.
}
We consider the following approximate system for $(t,x)\in (0,T)\times \mathbb{T}^d$:
\begin{subequations}\label{S-appr}
\begin{align}
  &\partial_t R+ \partial_{xx}iv(\vc{u} R) 

  =0,\label{S-appr1}\\
  &\partial_t Q + \partial_{xx}iv(\vc{u} Q) 

  =0 ,\label{S-appr2} \\
  &Q = \lr{1-\frac{R}{Z}} Z^\int_{\mathbb{R}}amma,q_{l}uad R\leq Z,\label{S-appr3}\\
   &{ {\rm div}\, \vc{u} = ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma_+} -\{({\cal T}_k(Z))^{\int_{\mathbb{R}}amma_+}\}} \label{S-appr4},\\
  &{\rm rot} \, \vc{u} = \vc{0},\label{S-appr5}\\
  &\intO{\vc{u}(t,x)}={0}.\label{S-appr6}
\varepsilonnd{align}
\varepsilonnd{subequations}
Note that combination of \varepsilonqref{S-appr5} with \varepsilonqref{S-appr6} defines $\vc{u}$ as a potential flow, i.e. there exists $\partialhi(t,x)$ such that $\vc{u}=\nablaabla \partialhi$, equivalently $\partial_{xx}elta \partialhi ={\rm div}\, \vc{u}$. 
The main result of this section is existence of solutions to the approximate system \varepsilonqref{S-appr} supplemented by the initial conditions
 \varepsilonq{\label{ini_ep}
 R_k {\cal E}rt_{t=0} = R_{0,k}, q_{l}uad Q_k {\cal E}rt_{t=0} = Q_{0,k},
 }
 and $ Z_k {\cal E}rt_{t=0} = Z_{0,k}$ defined by \varepsilonqref{RelationInit}.

\begin{thm}\label{approx}
Let the initial conditions \varepsilonqref{ini_ep} be  such that
\varepsilonqh{0< R_{0,k},\, Q_{0,k}<\infty \ \text{a.e. in } \mathbb{T}^d,}
and let they satisfy \varepsilonqref{RelationInit} and \varepsilonqref{eq_ini}. Then there exists a global in time weak solution to \varepsilonqref{S-appr} 
such that 
\varepsilonq{\label{est_kfixed}
&R_k^{-1},\ R_k,\ Q_k^{-1},\ Q_k,\ Z_k^{-1},\ Z_k\in L^\infty((0,T)\times\mathbb{T}^d),\\
& \partialt R_k+\vc{u}_k\overlineot\nabla R_k\in L^\infty((0,T)\times\mathbb{T}^d),\\
& \partialt Q_k+\vc{u}_k\overlineot\nabla Q_k\in L^\infty((0,T)\times\mathbb{T}^d),\\
 &\partialt Z_k+\vc{u}_k\overlineot\nabla Z_k\in L^\infty((0,T)\times\mathbb{T}^d),\\
& \nablaabla_x \vc{u}\in{L^\infty(0,T;BMO(\mathbb{T}^d))},q_{l}uad \partial_{xx}iv\vc{u}\in{L^\infty((0,T)\times\mathbb{T}^d)},q_{l}uad \partialartial_t \vc{u}\in{L^\infty(0,T;L^2(\mathbb{T}^d))},
}
where Equations \varepsilonqref{S-appr1}, \varepsilonqref{S-appr2}, and \varepsilonqref{S-appr4} are satisfied in ${{\cal D}'((0,T)\times \mathbb{T}^d})$, and the initial conditions \varepsilonqref{ini_ep}   are satisfied in  ${\cal D}'(\mathbb{T}^d)$.
\varepsilonnd{thm}

The rest of this section is devoted to the proof of this Theorem. It will be divided into three main steps:
\begin{enumerate}
\item[Step 1.] Proof of existence of solutions to the Lagrangian reformulation of system \varepsilonqref{S-appr}.
\item[Step 2.] Construction of  characteristics for certain regularization of the flow.
\item[Step 3.] Passage to the limit with regularization parameter and retrieval of the Eulerian formulation.
\varepsilonnd{enumerate}

\subsection{Proof of Theorem \ref{approx} -- Step 1}

The starting point for this section is system \varepsilonqref{S-appr} written in the  Lagrangian coordinates. We omit the definition of these coordinates  at this level on purpose, we will come back to this issue in the consecutive steps of the proof of Theorem \ref{approx}.

We consider the following system for  $ (t,y)\in (0,T)\times\mathbb{T}^d$:
\begin{subequations}\label{S-L}
\begin{align}
 & \partial_t r+ r \sigma

 =0,\label{S-L1}\\
 & \partial_t q + q \sigma 

 =0, \label{S-L2}\\
  &{\sigma =  ({\cal T}_k(z))^{\int_{\mathbb{R}}amma_+}-\{({\cal T}_k(z))^{\int_{\mathbb{R}}amma_+}\}_{{\cal L}}},\label{S-L3}\\
 & q = \lr{1- \frac{r}{z}}z^\int_{\mathbb{R}}amma,q_{l}uad r\leq z, \label{S-L4}
\varepsilonnd{align}
\varepsilonnd{subequations}
where by $\{\overlineot\}_{{\cal L}}$ we denote the average on the torus 
\varepsilonq{\label{defL}
\{f\}_{{\cal L}} := \frac{1}{|\mathbb{T}^d|}\int_{\mathbb{T}^d} f(t,y) \varepsilonxp \lr{\int_0^t \sigma (s,y) \, {\rm d}s}{\rm d} {y}.}
The unknowns of the \varepsilonqref{S-L} are $r=r(t,y)$, $q=q(t,y)$, $z=z(t,y)$,  and $\sigma=\sigma(t,y)$, whereas $z=z(r,q)$ is a unique solution to  \varepsilonqref{S-L4}, see the proof of Lemma \ref{Lemma_posit}.
We supplement the system \varepsilonqref{S-L} with the following initial conditions:
\varepsilonqh{
&r(0,y)=r_0(y),q_{l}uad q(0,y)=q_0(y),q_{l}uad 0< r_0,\ q_0<\infty\ \text{a.e. in } \mathbb{T}^d,\\
&r_0\leq z_0, q_{l}uad q_0 = \lr{1- \frac{r_0}{z_0}}z_0^\int_{\mathbb{R}}amma.
}
{Note that there are} no space derivatives in system \varepsilonqref{S-L}, which is the main advantage of the Lagrangian reformulation. In fact it allows to transform the PDE system to the system of ODEs with a nonlocal term $\{({\cal T}_k(z))^{\int_{\mathbb{R}}amma_+}\}_{{\cal L}}$. In addition to that, we know a-priori that  any solution of \varepsilonqref{S-L} satisfies the estimates
\begin{equation}\label{e1}
 \|\sigma\|_{L^\infty((0,T)\times\mathbb{T}^d)} \leq k^{\int_{\mathbb{R}}amma_+},
\varepsilonnd{equation}
and 
\begin{equation}\label{e2}
\sup_{t\in(0,T)} \|r,\,q,\,r^{-1},\,q^{-1}\|_{L^\infty(\mathbb{T}^d)}(t) \leq  \|r_0,q_0,r_0^{-1},q_0^{-1}\|_{L^\infty(\mathbb{T}^d)} \varepsilonxp(T k^{\int_{\mathbb{R}}amma_+} ).
\varepsilonnd{equation}
This information suggests the choice of functional space for the fixed point theorem that will be used to prove the existence of solutions to system \varepsilonqref{S-L}. In what follows we show that the map 
\varepsilonqh{
\Phi\ :\  L^\infty((0,T)\times\mathbb{T}^d)\times L^\infty((0,T)\times\mathbb{T}^d)\to  L^\infty((0,T)\times\mathbb{T}^d)\times L^\infty((0,T)\times\mathbb{T}^d),}
\varepsilonqh{ \Phi(r,q)=(\bar r, \bar q),
}
where $\bar r, \bar q$ are the solutions to the following system in $(0,T)\times\mathbb{T}^d$
\begin{subequations}\label{Ban}
\begin{align}
 & \partial_t \bar r+ r \sigma

 =0,\label{Ban1}\\
&  \partial_t \bar q + q \sigma 

=0, \label{Ban2}\\
  &{\sigma =  ({\cal T}_k(z))^{\int_{\mathbb{R}}amma_+}-\{({\cal T}_k(z))^{\int_{\mathbb{R}}amma_+}\}_{\cal L}},\label{Ban4}\\
  &q = \lr{1-\frac{r}{z}} z^\int_{\mathbb{R}}amma,q_{l}uad r\leq z, \label{Ban3}
\varepsilonnd{align}
\varepsilonnd{subequations}
is a contraction, at least for short time $t\in(0,T)$.
 
Let us denote 
\varepsilonqh{
 \bar{\partialelta}talta \bar r = \bar r_1 - \bar r_2, q_{l}uad \bar{\partialelta}talta \bar q = \bar q_1 - \bar q_2,
 q_{l}uad \bar{\partialelta}talta r = r_1-r_2, q_{l}uad \bar{\partialelta}talta q = q_1-q_2, q_{l}uad \bar{\partialelta}talta \sigma = \sigma_1 -\sigma_2.
}
Using \varepsilonqref{Ban} we compute 
\varepsilonqh{
(\bar{\partialelta}talta \overline r, \bar{\partialelta}talta \overline q)=  \Phi(r_1,q_1)-\Phi(r_2,q_2),
}
we have
\begin{subequations}\label{Ban-dif}
\begin{align}
&  \partial_t \bar{\partialelta}talta \bar r+  \sigma_1 \bar{\partialelta}talta r + r_2 \bar{\partialelta}talta \sigma 
        =0, \label{Ban-dif1}\\
&  \partial_t \bar{\partialelta}talta \bar q + \sigma_1 \bar{\partialelta}talta  q + q_2 \bar{\partialelta}talta \sigma 
=0, \label{Ban-dif2} \\
 & \bar{\partialelta}talta \sigma =  [({\cal T}_k(z_1))^{\int_{\mathbb{R}}amma_+}-({\cal T}_k(z_2))^{\int_{\mathbb{R}}amma_+}] -[\{({\cal T}_k(z_1))^{\int_{\mathbb{R}}amma_+}\}_{\cal L} -\{({\cal T}_k(z_2))^{\int_{\mathbb{R}}amma_+}\}_{\cal L}],\label{Ban-dif3}
\varepsilonnd{align}
\varepsilonnd{subequations}
where $z_1$ and $z_2$ are functions depending respectively on $(r_1,q_1)$ and $(r_2,q_2)$ through the nonlinear
implicit relation \varepsilonqref{Ban3}. Therefore, similarly to \varepsilonqref{pTR} we can show that
\varepsilonqh{
\|\partialartial_r {\cal T}_k(z)\|_{L^\infty((0,T)\times\mathbb{T}^d)}\leq C(k),q_{l}uad \mbox{and}q_{l}uad \|\partialartial_q {\cal T}_k(z)\|_{L^\infty((0,T)\times\mathbb{T}^d)}\leq C(k),
}
this implies, in particular that
\varepsilonq{\label{dsigmab}
\sup_{t\in(0,T)}\|\bar{\partialelta}talta{\cal T}_k(z)\|_{L^\infty(\mathbb{T}^d)}(t)\leq C(k)\sup_{t\in(0,T)}\lr{\|\bar{\partialelta}talta r\|_{L^\infty(\mathbb{T}^d)}(t) + \|\bar{\partialelta}talta q\|_{L^\infty(\mathbb{T}^d)}(t)}.
}
To apply the Banach fixed point to the map $\Phi$, we need to show that  $\| \bar{\partialelta}talta \bar r,\, \bar{\partialelta}talta \bar q\|_{L^\infty((0,T)\times\mathbb{T}^d)}\leq C\| \bar{\partialelta}talta  r,\, \bar{\partialelta}talta  q\|_{L^\infty((0,T)\times\mathbb{T}^d)}$ with some constant $C<1$. To this purpose we need to estimate $\bar{\partialelta}talta \sigma$ appearing in both equations, \varepsilonqref{Ban-dif1} and \varepsilonqref{Ban-dif2}, by $\bar{\partialelta}talta r$ and $\bar{\partialelta}talta q$. Analyzing equation \varepsilonqref{Ban-dif3}  we notice that since the first term on the r.h.s. can be treated using \varepsilonqref{dsigmab},  the only  challenge is due to the nonlocal term. Using \varepsilonqref{defL} and \varepsilonqref{dsigmab} we write
\varepsilonq{\label{e6}
&\left|\mathbb{T}^d| | \{({\cal T}_k(z_1))^{\int_{\mathbb{R}}amma_+}\}_{\cal L}(t) -\{({\cal T}_k(z_2))^{\int_{\mathbb{R}}amma_+}\}_{\cal L} (t)\right|\\

&=\left| \int_{\mathbb{T}^d} ({\cal T}_k(z(r_1,q_1)))^{\int_{\mathbb{R}}amma_+} \varepsilonxp\lr{\int_0^t \sigma_1(s,y) {\rm d}s}{\rm d} {y} \right.\\
&\hookrightarrowspace{5cm} \left.- \int_{\mathbb{T}^d} \lr{{\cal T}_k(z(r_2,q_2))}^{\int_{\mathbb{R}}amma_+} \varepsilonxp\lr{\int_0^t \sigma_2(s,y) {\rm d}s} {\rm d} {y}\right| \\

&\leq \left| \int_{\mathbb{T}^d} \varepsilonxp\lr{\int_0^t \sigma_1(s,y) {\rm d}s} (({\cal T}_k(z(r_1,q_1)))^{\int_{\mathbb{R}}amma_+}-({\cal T}_k(z(r_2,q_2)))^{\int_{\mathbb{R}}amma_+}) {\rm d} {y}\right| \\

&q_{l}quad- \left|\int_{\mathbb{T}^d} ({\cal T}_k(z(r_2,q_2)))^{\int_{\mathbb{R}}amma_+} \left[ \varepsilonxp\lr{\int_0^t \sigma_2(y,s) {\rm d}s}-  \varepsilonxp\lr{\int_0^t \sigma_1(s,y) {\rm d}s}\right]{\rm d} {y}\right| \\
&\leq C(k)\lr{\|\bar{\partialelta}talta r\|_{L^\infty(\mathbb{T}^d)}(t) + \|\bar{\partialelta}talta q\|_{L^\infty(\mathbb{T}^d)}(t)}+ C(k) \int_{\mathbb{T}^d} \left|\int_0^t \bar{\partialelta}talta \sigma (s,y) {\rm d}s\right|{\rm d} {y} \\
&\leq C(k)\lr{\|\bar{\partialelta}talta r\|_{L^\infty(\mathbb{T}^d)}(t) + \|\bar{\partialelta}talta q\|_{L^\infty(\mathbb{T}^d)}(t)} \\
&q_{l}quad+ C(k)\,  t \sup_{s\in(0,t)} \left| \{({\cal T}_k(z_1))^{\int_{\mathbb{R}}amma_+}\}_{\cal L}(s) -\{({\cal T}_k(z_2))^{\int_{\mathbb{R}}amma_+}\}_{\cal L}(s)\right|,
}
where the last inequality follows from the last equation in \varepsilonqref{Ban-dif}. Taking now supremum over time $t\in (0,\tau)$ on both sides of \varepsilonqref{e6}  for $\tau C(k) \leq 1/2 |\mathbb{T}^d|$ we find
\begin{equation}\label{e7}
 \sup_{t\in(0, \tau)} | \{({\cal T}_k(z_1))^{\int_{\mathbb{R}}amma_+}\}(t) -\{({\cal T}_k(z_2))^{\int_{\mathbb{R}}amma_+}\}(t)|\leq
 C(k)\,  \sup_{t\in(0,\tau)} \lr{\|\bar{\partialelta}talta r\|_{L^\infty(\mathbb{T}^d)}(t) + \|\bar{\partialelta}talta q\|_{L^\infty(\mathbb{T}^d)}(t)}.
\varepsilonnd{equation}
With this estimate at hand we can return to (\ref{Ban-dif})  and compute  that 
\begin{equation}\label{e8}
 \sup_{0\leq t\leq \tau} \lr{\|\bar{\partialelta}talta \bar r\|_{L^\infty(\mathbb{T}^d)}(t) + \|\bar{\partialelta}talta \bar q\|_{L^\infty(\mathbb{T}^d)}(t)}\leq
 C(k)t \sup_{0\leq t\leq \tau}\lr{\|\bar{\partialelta}talta r\|_{L^\infty(\mathbb{T}^d)}(t) + \|\bar{\partialelta}talta q\|_{L^\infty(\mathbb{T}^d)}(t)}.
\varepsilonnd{equation}
So, choosing $\tau$ small such that   $\tau C(k) <1$,  map $\Phi$ is a contraction on a time interval $[0,\tau]$. Since $\tau$ depends only on the truncation parameter $k$ which is constant,  we can iterate this procedure to obtain a unique solution to  (\ref{S-L}) on the whole time interval $(0,T)$. 

\subsection{Proof of Theorem \ref{approx} -- Step 2}
\subsubsection{Explanation of the strategy}\label{Sect_strategy}
The existence of solutions proved before corresponds, in some sense, to the  Lagrangian reformulation of  system \varepsilonqref{S-appr}. We now want to define the Eulerian coordinates and show that we can recover 
the velocity vector field at the level of the Eulerian coordinates.  Using the mathematical jargon, we intend to solve the equation
\varepsilonq{\label{heur1}
\partial_{xx}iv_x \vc{u}(t,x)=\sigma(t,y),}
which means that for given $\sigma$, we will find $\vc{u}$ and $x=x(t,y)$ such that the above equality is satisfied. Our candidate $x=x(t,y)$ is a solution to an ODE defining the Lagrangian transformation
\varepsilonq{\label{heur2}
\frac{d x}{dt} = \vc{u}(t,x), q_{l}quad x{\cal E}rt_{t=0} = y.}
Combination of \varepsilonqref{heur1} and \varepsilonqref{heur2} leads to the nonlinear PDE-ODE system which we intend to solve using the following variant of the Leray-Schauder fixed point theorem.
\begin{thm}
Let $\Phi$ be a continuous, compact mapping, $X$ a Banach Space. Let for any $\zeta\in[0,1]$ the fixed point $v=\zeta\Phi(v)$, $v\in X$ be bounded. Then $\Phi$ possesses at least one fixed point in $X$.
\varepsilonnd{thm}

Let us now explain how we intend to apply the above theorem. We will consider a regularization of $\sigma$
\varepsilonqh{
 \sigma_\bar{\partialelta}talta(t,y) = \sigma(t,y) \ast \bar{\kappa}ppa_\bar{\partialelta}talta(y),
}
where  $\bar{\kappa}ppa_\bar{\partialelta}talta$ is a standard mollifier and $\sigma$ is a known function -- the solution found in the previous step of the proof.  We further define a map $\Phi(\overline \vc{u}) = \vc{u}$
\varepsilonqh{
\Phi\ :\ {\cal C}(0,T; W^{2-\alpha,p}(\mathbb{T}^d))\to \ {\cal C}(0,T; W^{2-\alpha,p}(\mathbb{T}^d)),
}
for $\alpha>0$ arbitrary small, and {$p>1$}, in the following way: 

\nablaoindent  1. For a given $\overline \vc{u}\in L^p(0,t; W^{2-\alpha,p}(\mathbb{T}^d))$ we use the Cauchy-Lipschitz theorem to find a unique $x=x(t,y)$ such that  
\varepsilonq{\label{ODE}
\frac{d x(t,y)}{dt} = \overline \vc{u}(t,x(t,y)), q_{l}quad x(t,y){\cal E}rt_{t=0} = y.}

\nablaoindent 2. We then differentiate \varepsilonqref{ODE} with respect to $y$, to check that $H(t,y) = \frac{\partialartial x}{\partialartial y}(t,y)$ satisfies the equation
$$\partialt H(t,y)=\nabla_x \bar \vc{u}(t,x(t,y)) H(t,y),q_{l}uad H(t,y){\cal E}rt_{t=0}=Id.$$
Therefore, integrating in time, we get:
\varepsilonqh{
\varepsilonxp\lr{ - \int_0^T \|\nablaabla_x \bar \vc{u} \|_{L^\infty(\mathbb{T}^d)} {\rm d} t }\leq \left\|H\right\|_{L^\infty((0,T)\times\mathbb{T}^d)} \leq \varepsilonxp\lr{  \int_0^T \|\nablaabla_x \bar \vc{u} \|_{L^\infty(\mathbb{T}^d)} {\rm d} t }.
}
Note also that the determinant of $H(t,y)$, $J(t,y)=det \,\frac{\partialartial x}{\partialartial y}(t,y)$ satisfies the equation
\varepsilonq{\label{Jf}
\partialt J(t,y)=\partial_{xx}iv_x\bar \vc{u}(t,x(t,y)) J(t,y).
}
and so
\begin{equation}\label{j_est}
 \varepsilonxp\lr{ - \int_0^T \|\partial_{xx}iv_x \bar \vc{u}\|_{L^\infty(\mathbb{T}^d)} {\rm d} t } \leq \|J\|_{L^\infty((0,T)\times\mathbb{T}^d)} \leq \varepsilonxp\lr{ \int_0^T \|\partial_{xx}iv_x \bar \vc{u}\|_{L^\infty(\mathbb{T}^d)} {\rm d} t }.
\varepsilonnd{equation}
This means that $H(t,y)$ is invertible, moreover we have
\begin{equation}\label{f11}
 \left\|\frac{\partialartial y}{\partialartial x}\right\|_{L^\infty((0,T)\times\mathbb{T}^d)} \leq \varepsilonxp\lr{ \int_0^T \|\nablaabla \bar \vc{u}\|_{L^\infty(\mathbb{T}^d)} {\rm d} t }.
\varepsilonnd{equation}

\nablaoindent 3. Because $H$ is invertible we can express $y$ as a function of $t$ and $x$. For such $y=y(t,x)$ we will look for solution
\varepsilonq{\label{phid}
\vc{u}(t,x)=\nabla\partialhi(t,x),q_{l}uad \mbox{where} q_{l}uad \Delta_x \partialhi(t,x)=\sigma_\bar{\partialelta}talta(t,y(t,x)).
}
\begin{rmk}
This approach guarantees not only that $\partial_{xx}iv_x\vc{u}(t,x)=\sigma_\bar{\partialelta}talta(t,y(t,x))$, but also, that $\vc{u}$ has a structure of a gradient flow and will satisfy \varepsilonqref{S-appr5} and \varepsilonqref{S-appr6}.
\varepsilonnd{rmk}
\begin{rmk}
The solution $\vc{u}$ constructed above depends on the parameter $\bar{\partialelta}talta$ and should be denoted $\vc{u}^\bar{\partialelta}talta$, but we omit this index at this stage of the proof.
\varepsilonnd{rmk}

\subsubsection{A priori estimates} 
We assume here that $\vc{u}=\bar\vc{u} \in {\cal C}(0,T; W^{2-\alpha,p}(\mathbb{T}^d))$ such that there exists $\vc{u}(t,x)=\nabla_x \partialhi(t,x)$ and it satisfies
\varepsilonq{\label{phidl}
\Delta_x\partialhi(t,x)= \zeta \sigma_\bar{\partialelta}talta(t,y(t,x)),q_{l}uad {\zeta}\in[0,1].}
We will show that every solution of this equation is bounded in ${\cal C}(0,T; W^{2-\alpha,p}(\mathbb{T}^d))$ uniformly with respect to $\zeta $.

\nablaoindent{\it Estimates of space derivative.} First note that the standard elliptic estimate for \varepsilonqref{phidl} gives
\varepsilonqh{\sup_{t\in(0,T)}\|\nabla_x\partialhi(t,x)\|_{W^{1,p}(\mathbb{T}^d)}\leq \zeta C(p)\|\sigma_\bar{\partialelta}talta(t,y(t,x))\|_{L^\infty((0,T)\times \mathbb{T}^d)}
\leq C(p,k)
,}
and so 
\varepsilonq{\label{uBMO}
{\sup_{t\in(0,T)}\|\nabla\vc{u}\|_{BMO(\mathbb{T}^d)}\leq C(k)}
}
for any $p<\infty$. In particular, taking $p>d$, we obtain $W^{1,p}(\mathbb{T}^d)\hookrightarrowookrightarrow \hookrightarrowookrightarrow C(\Ov{\mathbb{T}^d})$.
Differentiating equation \varepsilonqref{phidl} with respect to $x_k$ gives
\begin{equation}\label{f10}
 \partial_{xx}elta_x \frac{\partialartial\partialhi(t,x)}{\partialartial x_k} = \nablaabla_y \sigma_\bar{\partialelta}talta(t,y(x,t)) \overlineot \frac{\partialartial y}{\partialartial x_k},
\varepsilonnd{equation}
therefore using the estimate \varepsilonqref{f11} with $\vc{u}=\bar\vc{u}$ we obtain
\varepsilonqh{\sup_{t\in(0,T)} \|\nablaabla_x \vc{u}\|_{W^{1,p}(\mathbb{T}^d)}
 \leq \zeta  C(p) \|\nablaabla_y \sigma_\bar{\partialelta}talta\|_{L^\infty((0,T)\times\mathbb{T}^d)} \varepsilonxp\lr{ T \sup_{t\in(0,T)} \|\nablaabla_x \vc{u}\|_{L^\infty(\mathbb{T}^d)} } .}
We now use the following estimate for  $p >d$  
\begin{equation}\label{f12}
 \|\nablaabla f\|_{L^\infty(\mathbb{T}^d)} \leq C(p)(1+\|\nablaabla f\|_{BMO(\mathbb{T}^d)} \lr{1+\ln^+ (\|\nablaabla f\|_{W^{1,p}(\mathbb{T}^d)} + \|f\|_{L^\infty(\mathbb{T}^d)}) }^{1/2}).
\varepsilonnd{equation}
The proof of above inequality on can find in \cite{Ogawa} in Corollary 2.4 inequality (2.6). The original proof holds for the whole $\mathbb{R}^d$,  but taking an extension operator $E:\mathbb{T}^d \sim [0,1]^d \to \mathbb{R}^d$ preserving regularity in $W^{2,p}$ we get the same  result on the torus.
Together with \varepsilonqref{uBMO} and $ \zeta \leq 1$, we can write
\varepsilonqh{
 &\sup_{t\in(0,T)} \|\nablaabla_x \vc{u}\|_{W^{1,p}(\mathbb{T}^d)} \\
 &
  \leq C(p, \bar{\partialelta}talta) \varepsilonxp\lr{ T \|\nablaabla_x \vc{u}\|_{L^\infty(0,T; BMO(\mathbb{T}^d))}(1+\ln^+ \|\nablaabla_x \vc{u}\|_{L^\infty(0,T; W^{1,p}(\mathbb{T}^d))})^{1/2}} \\
& \leq C(p,\bar{\partialelta}talta) \varepsilonxp\lr{ \frac{1}{2} \ln^+ \|\nablaabla_x \vc{u}\|_{L^\infty(0,T; W^{1,p}(\mathbb{T}^d))} + C(1+T^2)} \\
& \leq C(p,\bar{\partialelta}talta) e^{C(1+T^2)} \|\nablaabla_x \vc{u}\|^{1/2}_{L^\infty(0,T; W^{1,p}(\mathbb{T}^d))}.
}
We therefore have
\begin{equation}\label{est_ud}
  \sup_{t\in(0,T)} \|\nablaabla_x\vc{u}\|_{W^{1,p}(\mathbb{T}^d)}  \leq C.
\varepsilonnd{equation}
 In this way we find the {\it a-priori} information about $\vc{u}$:
\varepsilonqh{
 \vc{u}(t,x) \in L^\infty(0,T; W^{2,p}(\mathbb{T}^d)),q_{l}quad {\rm div}\,_x \vc{u}(t,x) = \sigma \ast \bar{\kappa}ppa_\bar{\partialelta}talta (t,y(t,x)).}
The estimates uniform in $\bar{\partialelta}talta$ are
\begin{equation}\label{est_u}
 \|\nablaabla_x \vc{u}\|_{L^\infty(0,T;BMO(\mathbb{T}^d))} + \|\partial_{xx}iv_x \vc{u}\|_{L^\infty((0,T)\times\mathbb{T}^d)} \leq C.
\varepsilonnd{equation}

\nablaoindent {\it Estimate  of the time derivative.} We want to check the time-regularity of $\vc{u}(t,x)$. For fixed $\bar{\partialelta}talta$ we expect better information, however here we present only the estimates uniform with respect to $\bar{\partialelta}talta$. 
The therefore switch to the weak formulation  of (\ref{phidl}):
$$
\intO{ \partial_{xx}elta_x \partialhi (t,x) \partiali(x)} = \zeta \intO{\sigma_\bar{\partialelta}talta (t,y(x,t))\partiali(x)}= 
\zeta \int_{\mathbb{T}^d}\sigma_\bar{\partialelta}talta(t,y) \partiali(x(t,y)) J(t,y) \, {\rm d} {y},
$$
where $\partiali$ is a smooth function on $\mathbb{T}^d$. We now differentiate this identity in time 
\varepsilonq{\label{ptsigma}
\partial_{xx}t  \intO{\partial_{xx}elta_x \partialhi (t,x) \partiali(x)} =&  
\zeta \int_{\mathbb{T}^d} \partialt\left[\sigma_\bar{\partialelta}talta(t,y) \partiali(x(t,y)) J(t,y)\right]\, {\rm d} {y}\\
=& \zeta  \int_{\mathbb{T}^d}\partialt\sigma_\bar{\partialelta}talta(t,y)\partiali(x(t,y)) J(t,y)\,{\rm d} {y}
+\zeta  \int_{\mathbb{T}^d}\sigma_\bar{\partialelta}talta(t,y)\partialt\partiali(x(t,y)) J(t,y)\,{\rm d} {y}\\
&q_{l}uad+ \zeta  \int_{\mathbb{T}^d}\sigma_\bar{\partialelta}talta(t,y)\partiali(x(t,y)) \partialt J(t,y)\,{\rm d} {y}.
}
We now need to estimate all terms on the r.h.s. of \varepsilonqref{ptsigma}. For the first one, we  differentiate \varepsilonqref{S-L4}
in time, use expressions \varepsilonqref{S-L1} and \varepsilonqref{S-L2} for $\partialt r$ and $\partialt q$, and proceed as in \varepsilonqref{dsigmab} to check  that
$$
|\partialt \sigma_\bar{\partialelta}talta |=|(\partialt \sigma)\ast\bar{\kappa}ppa_\bar{\partialelta}talta| \leq |\partialt \sigma|  \leq C(k)( |r_t|+|q_t|) \in L^\infty((0,T)\times\mathbb{T}^d).
$$

Next, for the second term on the r.h.s. of \varepsilonqref{ptsigma} we use  \varepsilonqref{ODE} with $\bar\vc{u}=\vc{u}$ to write
$$
\partialt \partiali(x(t,y))= \nablaabla_x \partiali \frac{dx}{dt}= \vc{u}(t,x(t,y)) \nablaabla_x \partiali(x(t,y)),
$$
so, thanks to \varepsilonqref{est_u}, $\partialt \partiali(x(t,y))$ is bounded in $L^\infty(0,T;L^2(\mathbb{T}^d))$.  At last, the formula \varepsilonqref{Jf} with $\bar\vc{u}=\vc{u}$ together with \varepsilonqref{j_est} provides that $\partialt J(t,y)$ is bounded in $L^\infty((0,T)\times\mathbb{T}^d)$, so the third term on the r.h.s. of \varepsilonqref{ptsigma} is also bounded.

Summarizing, we  get that $\Delta_x\partialt\partialhi$ is bounded in $L^\infty(0,T;W^{-1,2}(\mathbb{T}^d))$ uniformly with respect to $\bar{\partialelta}talta$. In particular, using the Helmholtz decomposition, we deduce  that 
\varepsilonq{\label{est_ptu}
\|\partialt\vc{u}\|_{L^\infty(0,T;L^2(\mathbb{T}^d))}=\|\nabla_x\partialt\partialhi\|_{ L^\infty(0,T;L^2(\mathbb{T}^d))}\leq C,} 
with the constant $C$ that does not depend on $\bar{\partialelta}talta$.
 Obviously, combining the estimate \varepsilonqref{est_ptu} with \varepsilonqref{est_ud} and \varepsilonqref{est_u}, we verify that any fixed point satisfying \varepsilonqref{phidl} is uniformly bounded in ${\cal C}(0,T; W^{2-\alpha,p}(\mathbb{T}^d))$, independently of $\zeta $.

\subsubsection{The fixed point argument} 
We are now ready to proceed with the fixed point argument explained in Section \ref{Sect_strategy}.
We therefore take $\bar\vc{u} \in {\cal C}(0,T;W^{2-\alpha,p}(\mathbb{T}^d))$, and show that the operator $\Phi(\bar\vc{u})=\vc{u}$ defined through \varepsilonqref{phid} and \varepsilonqref{ODE} is continuous and compact in ${\cal C}(0,T;W^{2-\alpha,p}(\mathbb{T}^d))$.

Compactness is straightforward. Indeed, taking $\zeta =1$ in the system \varepsilonqref{phidl} we see that our a-priori estimates for $\vc{u}$, estimates \varepsilonqref{est_ptu} with \varepsilonqref{est_ud} and \varepsilonqref{est_u}, stay in force. Hence, on account of the Aubin-Lions lemma, the map $\Phi$ is compact. 

We now check the continuity of the map $\Phi$, by investigating the difference of two solutions $\vc{u}_1=\nabla_x\partialhi_1$ and $\vc{u}_2=\nabla_x\partialhi_2$
\varepsilonq{\label{conv_phi}
\intO{ |\Delta_x \partialhi_1(t,x) - \Delta_x \partialhi_2(t,x)|^p} &= \intO{|\sigma_\bar{\partialelta}talta(t,y_1(t,x)) - \sigma_\bar{\partialelta}talta(t, y_2(t,x))|^p}\\
 &\leq C\|\nabla_y\sigma_\bar{\partialelta}talta(t,y(t,x))\|^p_{L^\infty(\mathbb{T}^d)}\| y_1(t,x) -  y_2(t,x)\|^p_{L^\infty(\mathbb{T}^d)}.
}
Recalling \varepsilonqref{ODE}, we  note that for the regular characteristics intersecting in the point $(t,x)$ we have
\varepsilonqh{
 \|  y_1(x,t) -  y_2(x,t)\|_{L^\infty(\mathbb{T}^d)}\leq \int_0^t \|\bar\vc{u}_1 - \bar\vc{u}_2\|_{L^\infty(\mathbb{T}^d)} dt' .
}
So if $(\bar\vc{u}_1-\bar\vc{u}_2) \to 0$ in ${\cal C}(0,T;W^{2-\alpha,p}(\mathbb{T}^d))$, then of course
$(\bar\vc{u}_1-\bar\vc{u}_2) \to 0$ in ${\cal C}(0,T;L^{\infty}(\mathbb{T}^d))$, which thanks to \varepsilonqref{conv_phi} implies that 
$(\partialhi_1 - \partialhi_2) \to 0$ in ${\cal C}(0,T;W^{2,p}(\mathbb{T}^d))$. Therefore, using boundedness of $\vc{u}_1,\ \vc{u}_2$ in  $L^\infty(0,T;W^{2,p}(\mathbb{T}^d))$, and interpolation, we deduce that $(\vc{u}_1-\vc{u}_2) \to 0$ in $L^p(0,T;W^{2-\alpha,p}(\mathbb{T}^d))$. By Leray-Schauder fixed point theorem we have the existence for approximative system for $\bar{\partialelta}talta$ fixed.

\subsection{Proof of Theorem \ref{approx} -- Step 3}
 We want to let  $\bar{\partialelta}talta\to 0$ in the equation 
\begin{equation}\label{f20}
 \partial_{xx}elta_x \partialhi^\bar{\partialelta}talta(t,x^\bar{\partialelta}talta) = \sigma_\bar{\partialelta}talta(t,y(t,x^\bar{\partialelta}talta))
\varepsilonnd{equation}
remembering that  $\vc{u}^\bar{\partialelta}talta(t,x^\bar{\partialelta}talta)=\nabla_x\partialhi^\bar{\partialelta}talta(t,x^\bar{\partialelta}talta)$ and $x^\bar{\partialelta}talta $ is associated with the flow $\vc{u}^\bar{\partialelta}talta$ via \varepsilonqref{ODE} with $\bar\vc{u}=\vc{u}^\bar{\partialelta}talta$. We know that uniformly with respect to $\bar{\partialelta}talta$ we have
\varepsilonqh{
 \|\nablaabla_x \vc{u}^\bar{\partialelta}talta\|_{L^\infty(0,T;BMO(\mathbb{T}^d))} + \|\partial_{xx}iv_x \vc{u}^\bar{\partialelta}talta\|_{L^\infty((0,T)\times\mathbb{T}^d)}  + \|\partialartial_t \vc{u}^\bar{\partialelta}talta\|_{L^\infty(0,T;L^2(\mathbb{T}^d))} \leq C. 
}
We use again the weak form of (\ref{f20}), we have
\varepsilonq{\label{weak_sold}
 \intT{\!\!\!\int_{\mathbb{T}^d}\partial_{xx}elta_x \partialhi^\bar{\partialelta}talta (t,x^\bar{\partialelta}talta) {\rm x}i (t,x^\bar{\partialelta}talta) {\rm d} {x}^\bar{\partialelta}talta} = \intT{\!\!\!\int_{\mathbb{T}^d} \sigma \ast \bar{\kappa}ppa_\bar{\partialelta}talta (t,y(x^\bar{\partialelta}talta,t)) {\rm x}i(t,x^\bar{\partialelta}talta){\rm d} {x}^\bar{\partialelta}talta} \\
 = \intT{\!\!\!\int_{\mathbb{T}^d} \sigma \ast \bar{\kappa}ppa_\bar{\partialelta}talta (t,y) {\rm x}i(t,x^\bar{\partialelta}talta(t,y)) J^\bar{\partialelta}talta(t,y) {\rm d} {y}},
}
for any smooth ${\rm x}i$, where we denoted
$$J^\bar{\partialelta}talta(t,y)=\varepsilonxp\lr{ \int_0^t \partial_{xx}iv_x \vc{u}^\bar{\partialelta}talta(t',x^\bar{\partialelta}talta(t',y)) {\rm d}t'}=\varepsilonxp\lr{ \int_0^t \sigma \ast \bar{\kappa}ppa_\bar{\partialelta}talta (t',y) {\rm d}t'}.$$
We see $\partialhi^\bar{\partialelta}talta \to \partialhi$ weakly in $L^\infty(0,T;W^{2,p}(\mathbb{T}^d))$, and
$\sigma\ast \bar{\kappa}ppa_\bar{\partialelta}talta \to \sigma$ a.e. in $(0,T)\times\mathbb{T}^d$.

Using Crippa-Dellelis result  from \cite[Theorem 2.9 (stability of the flow)]{CrDe} saying
\varepsilonqh{
 \sup_{t\in(0,T)} \| x(t,y) -  x^\bar{\partialelta}talta(t,y)\|_{L^1(\mathbb{T}^d)} \leq C | \ln ( \| \vc{u} -\vc{u}^\bar{\partialelta}talta\|_{L^1((0,T) \times \mathbb{T}^d)})|^{-1}
}
we therefore get
$$
x^\bar{\partialelta}talta(t,y) \to x(t,y) \mbox{ in } L^\infty(0,T;L^1(\mathbb{T}^d)),
$$
since $\vc{u}^\bar{\partialelta}talta$ is compact in $L^1((0,T)\times \mathbb{T}^d)$.
Hence we can let $\bar{\partialelta}talta\to0$ in both sides of \varepsilonqref{weak_sold} to obtain
\varepsilonqh{
 \intTO{\partial_{xx}elta_x \partialhi (t,x) {\rm x}i (t,x)} =  \intT{\!\!\!\int_{\mathbb{T}^d} \sigma  (t,y) {{\rm x}i(t,x(y,t)) }J(t,y)\, {\rm d} {y}.}
}
Now we substitute for $\sigma$ using \varepsilonqref{Ban4}, and use Lemma 3.1. from \cite{CoCrSp} to write the weak formulation of the momentum equation in the Eulerian coordinates
{\varepsilonq{&\intTO{\partial_{xx}iv_x \vc{u} (t,x) {\rm x}i (t,x)}\label{relation} \\
&q_{l}uad= \intTO{\Big(({\cal T}_k(Z(t,x)))^{\int_{\mathbb{R}}amma_+}-\left\{({\cal T}_k(Z(t,x)))^{\int_{\mathbb{R}}amma_+}\right\}\Big){\rm x}i(t,x)},}}
where $Z(t,x)$ is defined by $Z(t,x(t,y))=z(t,y)$. Defining $R(t,x(t,y))=r(t,y)$ and $Q(t,x(t,y))=q(t,y)$ we can also pass to the limit in the weak formulation of \varepsilonqref{S-appr1} and \varepsilonqref{S-appr2}. Note that for the limit case, characteristics $x(t,y)$ are now  well defined. So, solutions given in the Lagrangian coordinates give us weak solutions to the original approximate problem (\ref{S-appr}). 

\section{Existence of solutions to system \varepsilonqref{S}}
Having proven the existence of solutions to the approximate system we now intend to pass to the limit with the approximation parameter $k\to\infty$,  recovering the system \varepsilonqref{S} and concluding the proof of Theorem \ref{Thm:main}. This section will also include rigorous justification of the compactness result from the previous section at the level of the approximate system.

From Theorem \ref{approx} it follows that there exists a sequence $\{R_k,Q_k,Z_k,\vc{u}_k\}_{k=1}^\infty$ satisfying  system \varepsilonqref{S-appr} in the sense of distributions belonging to a class \varepsilonqref{est_kfixed}.
Note, however, that this is not information uniform with respect to $k$. 
\subsection{Uniform estimates}
Although \varepsilonqref{est_kfixed} is not uniform with respect to $k$, we can still use \varepsilonqref{e2} from the previous section to deduce that
$$R_k\int_{\mathbb{R}}eq 0,q_{l}uad Q_k\int_{\mathbb{R}}eq 0,$$
uniformly with respect to $k$. Moreover, repeating the proof  of Lemma \ref{Lemma_posit}, we find unique $Z_k$ defined via \varepsilonqref{S-appr3}, therefore uniformly w.r.t. $k$ we also have that
$$R_k\leq Z_k.$$
In the rest of this section we will skip the index $k$ where no confusion can arise. 

 
 \subsubsection{Estimates of $({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+}$ uniform w.r.t.  $k$}
 
We use the renormalized equation for $Z$, \varepsilonqref{SZ}, that can be derived from \varepsilonqref{S-appr1} and \varepsilonqref{S-appr2} using \varepsilonqref{S-appr3}. Testing this equation respectively  by $(\int_{\mathbb{R}}amma_+-1) ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma_+-2} {\cal T}'_k(Z)$ and by $(\int_{\mathbb{R}}amma_+-\int_{\mathbb{R}}amma) ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma_+-\int_{\mathbb{R}}amma-1} {\cal T}'_k(Z)$, we obtain
\varepsilonqh{
&\partialartial_t ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1} + \vc{u} \overlineot \nablaabla ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1} 
     + (\int_{\mathbb{R}}amma^+ -1) ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-2}{\cal T}'_k(Z) Z\partial_{xx}iv \vc{u}\\
&  \hookrightarrowskip3cm    +  (\int_{\mathbb{R}}amma^+ -1) \frac{(1-\int_{\mathbb{R}}amma)(Z-R) ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-2}{\cal T}'_k(Z)Z  }{\int_{\mathbb{R}}amma (Z-R) + R} \partial_{xx}iv \vc{u}
= 0,
}
and
\varepsilonqh{
&\partialartial_t ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma} + \vc{u} \overlineot \nablaabla ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma} 
     + (\int_{\mathbb{R}}amma^+ -\int_{\mathbb{R}}amma) ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma-1}{\cal T}'_k(Z)Z \partial_{xx}iv \vc{u}\\
&  \hookrightarrowskip3cm    +  (\int_{\mathbb{R}}amma^+ -\int_{\mathbb{R}}amma) \frac{(1-\int_{\mathbb{R}}amma)(Z-R) ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma-1}{\cal T}'_k(Z)Z  }{\int_{\mathbb{R}}amma (Z-R) + R} \partial_{xx}iv \vc{u} 
= 0.
}
After simplification we get
\varepsilonqh{
&\partialartial_t ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1} + \vc{u} \overlineot \nablaabla ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1} 
     \\
&  \hookrightarrowskip3cm    +  (\int_{\mathbb{R}}amma^+ -1) \frac{Z^{\int_{\mathbb{R}}amma^+}  }{\int_{\mathbb{R}}amma (Z-R) + R} \partial_{xx}iv \vc{u} 1_{Z\leq k}
= 0,
}
and
\varepsilonqh{
&\partialartial_t ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma} + \vc{u} \overlineot \nablaabla ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma} \\
&  \hookrightarrowskip3cm    +  (\int_{\mathbb{R}}amma^+ -\int_{\mathbb{R}}amma) \frac{Z^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma+1}}{\int_{\mathbb{R}}amma (Z-R) + R} \partial_{xx}iv \vc{u}1_{Z\leq k}
= 0.
}
If we multiply the first equation by $R$ and the second by $Q$ and use the continuity equations \varepsilonqref{S-appr1} and \varepsilonqref{S-appr2} we get
\varepsilonq{\label{RTk}
&\partialartial_t \lr{R ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1}} + \partial_{xx}iv\lr{R ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1} \vc{u}}
     \\
&  \hookrightarrowskip3cm    +  (\int_{\mathbb{R}}amma^+ -1) \frac{RZ^{\int_{\mathbb{R}}amma^+}  }{\int_{\mathbb{R}}amma (Z-R) + R} \partial_{xx}iv \vc{u} 1_{Z\leq k}
       =0,
}
and
\varepsilonq{\label{QTk}
&\partialartial_t \lr{Q ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma}} + \partial_{xx}iv(Q \lr{{\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma} \vc{u}} \\
&  \hookrightarrowskip3cm    +  (\int_{\mathbb{R}}amma^+ -\int_{\mathbb{R}}amma) \frac{QZ^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma+1}}{\int_{\mathbb{R}}amma (Z-R) + R} \partial_{xx}iv \vc{u} 1_{Z\leq k}
       =0.
}
Multiplying equations \varepsilonqref{RTk}, \varepsilonqref{QTk} by $\frac{1}{(\int_{\mathbb{R}}amma^+-1)}$, $\frac{1}{(\int_{\mathbb{R}}amma^--1)}$, respectively and noticing that $\frac{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma}{(\int_{\mathbb{R}}amma^--1)}=\int_{\mathbb{R}}amma$ we obtain 
\varepsilonqh{
&\frac{1}{(\int_{\mathbb{R}}amma^+-1)} \partialt\lr{ R \, ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1}}  + \frac{1}{(\int_{\mathbb{R}}amma^--1)} \partialt\lr{  Q ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma} }\\
&+\frac{1}{(\int_{\mathbb{R}}amma^+-1)} \partial_{xx}iv\lr{R ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1} \vc{u}}+\frac{1}{(\int_{\mathbb{R}}amma^--1)} \partial_{xx}iv\lr{Q( {\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma} \vc{u}} \\
&   + Z^{\int_{\mathbb{R}}amma^+}\frac{R+\int_{\mathbb{R}}amma Q Z^{1-\int_{\mathbb{R}}amma}}{\int_{\mathbb{R}}amma (Z-R) + R}\partial_{xx}iv \vc{u} 1_{Z\leq k}
\leq 0,
}
which, after integrating with respect to space, and noting that $Q Z^{1-\int_{\mathbb{R}}amma} = Z- R$, gives
\varepsilonqh{
&\frac{1}{(\int_{\mathbb{R}}amma^+-1)} \partial_{xx}t\intO{ R \, ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1}}  + \frac{1}{(\int_{\mathbb{R}}amma^--1)} \partial_{xx}t \intO{  Q ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma} }\\
   &q_{l}uad+ \intO{ \partial_{xx}iv \vc{u}\,  Z^{\int_{\mathbb{R}}amma^+}1_{Z\leq k}}
   \leq 0.
}
Let us note that $Z^{\int_{\mathbb{R}}amma^+}1_{Z\leq k}=({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+}$ and $\left\{{\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+}\right\}$ is constant in space, therefore
$$\intO{ \partial_{xx}iv \vc{u}\,  Z^{\int_{\mathbb{R}}amma^+}1_{Z\leq k}}=\intO{ \partial_{xx}iv \vc{u}\,  \lr{Z^{\int_{\mathbb{R}}amma^+}1_{Z\leq k}-\left\{{\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+}\right\}}}.$$
And so, using density of smooth functions in $L^2$, and taking in  \varepsilonqref{relation} ${\rm x}i=\partial_{xx}iv\vc{u}$, we get
\varepsilonq{\label{en_approx}
&\frac{1}{(\int_{\mathbb{R}}amma^+-1)} \partial_{xx}t\intO{ R \, ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1}}  + \frac{1}{(\int_{\mathbb{R}}amma^--1)} \partial_{xx}t \intO{  Q ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma} }\\
   &q_{l}uad + \intO{\partial_{xx}iv^2 \vc{u}} \leq 0,
}
hence, 
\begin{equation}\label{en_approx-2}
\partial_{xx}t \left[ \int { \frac{1}{(\int_{\mathbb{R}}amma^+-1)} { R \, ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1}}  + \frac{1}{(\int_{\mathbb{R}}amma^--1)} {  Q ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma} }}\right]\leq 0.
\varepsilonnd{equation}
We now observe that
\varepsilonq{
&R  ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1}+Q ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma}\\
&=R  ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1}+(Z-R)Z^{\int_{\mathbb{R}}amma-1} ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma}\\\
&= Z^{\int_{\mathbb{R}}amma} ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma} + R[k^{\int_{\mathbb{R}}amma-1}-Z^{\int_{\mathbb{R}}amma-1}]1_{Z\int_{\mathbb{R}}eq k} ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma}.
}
Since $\int_{\mathbb{R}}amma-1<0$, the last term is nonnegative, therefore 
\varepsilonq{\label{ineq}
({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+} \le 
 R  ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-1}+Q ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+-\int_{\mathbb{R}}amma}.
 }
Therefore, using \varepsilonqref{ineq} and integrating \varepsilonqref{en_approx-2} over time, we obtain
\varepsilonq{\label{Tk!!!}
\|T_k(Z)\|_{L^\infty(0,T;L^{\int_{\mathbb{R}}amma^+}(\mathbb{T}^d))}\leq C,
}
uniformly w.r.t. $k$. 

\subsubsection{Uniform estimates of $R^{\int_{\mathbb{R}}amma^+}$ and on $Z^{\int_{\mathbb{R}}amma^+}$}
Recall that the equation on $R$ reads
$$\partialartial_t R + \partial_{xx}iv (R\vc{u}) 
     = 0,$$
and is satisfied in the sense of distributions on $(0,T)\times\mathbb{T}^d$. Due to \varepsilonqref{est_kfixed}, this equation is satisfied in the renormalized sense, therefore we have
$$\partialartial_t R^{\int_{\mathbb{R}}amma^+} + \partial_{xx}iv (R^{\int_{\mathbb{R}}amma^+} \vc{u}) 
       + (\int_{\mathbb{R}}amma^+-1)\, R^{\int_{\mathbb{R}}amma^+} \partial_{xx}iv\vc{u} 
     = 0.$$
{ Using \varepsilonqref{S-appr4} and integrating} with respect to space, it gives
 \varepsilonqh{
 & \frac{d}{dt} \intO{R^{\int_{\mathbb{R}}amma^+}} + (\int_{\mathbb{R}}amma^+-1) \intO{R^{\int_{\mathbb{R}}amma^+} ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+}}\\
& \hookrightarrowskip3cm =
      (\int_{\mathbb{R}}amma^+-1) \intO{R^{\int_{\mathbb{R}}amma^+}} \intO{({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+}}. 
 } 
 Using \varepsilonqref{Tk!!!}, we get  through Gronwall Lemma 
\varepsilonq{\label{R!!!}
\|R\|_{L^\infty(0,T;L^{\int_{\mathbb{R}}amma^+}(\mathbb{T}^d))}+\|R^{\int_{\mathbb{R}}amma^+} ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+}\|_{L^1((0,T)\times\mathbb{T}^d)}\leq C,
}
 uniformly w.r.t. $k$. Performing exactly the same procedure for $Z$ we obtain
 \varepsilonq{\label{Z!!!}
\|Z\|_{L^\infty(0,T;L^{\int_{\mathbb{R}}amma^+}(\mathbb{T}^d))}+\|{\cal T}_k(Z)\|_{L^{2\int_{\mathbb{R}}amma^+}((0,T)\times\mathbb{T}^d)}\leq C.
}
 These bounds are important to pass to the limit in the system (\ref{S-appr}) with respect
 to $k$. Note that the bound on $Z$ provides the same bound on
 $({\cal T}_k(Z))$. 
Finally, renormalizing the approximate equation for $Q$ we easily deduce that
\varepsilonq{\label{Q!!!}
\|Q\|_{L^\infty(0,T;L^{\int_{\mathbb{R}}amma^-}(\mathbb{T}^d))}+\|Q^{\int_{\mathbb{R}}amma^-} ({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+}\|_{L^1((0,T)\times\mathbb{T}^d)}\leq C.
}
\subsubsection{Uniform estimates of $\vc{u}$}
With the estimate \varepsilonqref{Z!!!} at hand we can now estimate $\partial_{xx}iv\vc{u}$ directly from \varepsilonqref{S-appr4} , we have 
$$\|\partial_{xx}iv\vc{u}\|_{L^\infty(0,T;L^{1}(\mathbb{T}^d))}+\|\partial_{xx}iv\vc{u}\|_{L^{2}((0,T)\times\mathbb{T}^d)}\leq C.$$
Using the fact that ${\rm rot}\vc{u}=0$, we therefore find that
 \varepsilonq{\label{v!!!}
\|\vc{u}\|_{L^{2}(0,T; W^{1,2}(\mathbb{T}^d))}\leq C.}
 

 
\subsection{Compactness argument}
The uniform estimates from the previous section are sufficient to perform the limit passage $k\to\infty$ in all linear terms of the approximate system \varepsilonqref{S-appr}. Estimating the time derivatives of $R_k$ and $Q_k$ from equations \varepsilonqref{S-appr1}, \varepsilonqref{S-appr2}, respectively,
and using the uniform estimate \varepsilonqref{v!!!} together with a Div-Curl type argument, we justify that $R_k\vc{u}_k$ converges to $R\vc{u}$ and $Q_k\vc{u}_k$ converges to $Q\vc{u}$ in the sense of distributions.

The last problem to solve is to pass to the limit  in the nonlinear term in the momentum equation \varepsilonqref{S-appr4} -- the pressure. This requires a strong convergence of the sequence $\{Z_k\}_{k=1}^\infty$ approximating $Z$. For the moment we only know that $Z_k$ converges weakly$^*$ to $Z$ in ${L^\infty(0,T;L^{\int_{\mathbb{R}}amma^+}(\mathbb{T}^d))}$. To improve this convergence, we adapt the compactness criterion presented in Section \ref{Sect:3}. 

We first justify how to obtain the equivalent of equations \varepsilonqref{eqR} and \varepsilonqref{eqZ} on the approximate level. First, we write the equation for $R_x-R_y$ using \varepsilonqref{S-appr1} we get
\varepsilonqh{
&\partialt (R_x-R_y) + \partial_{xx}iv_x (\vc{u}_x\lr{R_x-R_y}) +
\partial_{xx}iv_y (\vc{u}_y\lr{R_x-R_y})  \\
&=\frac 12 (\partial_{xx}iv_x \vc{u}_x + \partial_{xx}iv_y \vc{u}_y) \lr{R_x-R_y} \\
&q_{l}uad - \frac 12 (\partial_{xx}iv_x \vc{u}_x-\partial_{xx}iv_y \vc{u}_y)(R_x+R_y).
}
Regularizing this equation over the space variables, i.e. testing the equation by ${\rm x}i_\varepsilonta(x-\overlineot){\rm x}i_\varepsilonta(y-\overlineot)$, where ${\rm x}i_\varepsilonta$ is a standard mollifier, and denoting $S_\varepsilonta[f]={\rm x}i\ast_y({\rm x}i\ast_x f)$ we obtain
\varepsilonq{\label{eqRdiff}
&\partialt S_\varepsilonta[R_x-R_y] + \partial_{xx}iv_x (\vc{u}_xS_\varepsilonta[R_x-R_y]) +
\partial_{xx}iv_y (\vc{u}_yS_\varepsilonta[R_x-R_y])  \\
&=r_\varepsilonta+\frac 12 S_\varepsilonta\left[(\partial_{xx}iv_x \vc{u}_x + \partial_{xx}iv_y \vc{u}_y) \lr{R_x-R_y} \right]\\
&q_{l}uad - \frac 12 S_\varepsilonta\left[(\partial_{xx}iv_x \vc{u}_x-\partial_{xx}iv_y \vc{u}_y)(R_x+R_y)\right],
}
where on account of the Friedrichs commutator lemma, $r_\varepsilonta\to 0$ in $L^p((0,T)\times\mathbb{T}^{2d})$ for any $p<\infty$.
We then multiply \varepsilonqref{eqRdiff} by
\varepsilonqh{
S_\varepsilonta[R_x-R_y] \lr{S_\varepsilonta[R_x-R_y] ^2+\bar{\partialelta}talta}^{\frac{\beta-2}{2}},
}
where $\bar{\partialelta}talta>0$. Note that on account of \varepsilonqref{est_kfixed} this function is bounded in $L^\infty((0,T)\times\mathbb{T}^d)$ uniformly w.r.t. $\varepsilonta$. 
We therefore obtain
\varepsilonq{\label{eqR_ren}
&\frac{1}{\beta}\partialt\lr{ \lr{S_\varepsilonta[R_x-R_y]}^2+\bar{\partialelta}talta}^{\frac\beta2}\\
&q_{l}uad +\frac{1}{\beta} \partial_{xx}iv_x \lr{\lr{ \lr{S_\varepsilonta[R_x-R_y]}^2+\bar{\partialelta}talta}^{\frac\beta2}\vc{u}_x} +
\frac{1}{\beta}\partial_{xx}iv_y \lr{\lr{ \lr{S_\varepsilonta[R_x-R_y]}^2+\bar{\partialelta}talta}^{\frac\beta2}\vc{u}_y}  \\
&=\Big\{r_\varepsilonta+\frac 12 S_\varepsilonta\left[(\partial_{xx}iv_x \vc{u}_x + \partial_{xx}iv_y \vc{u}_y) \lr{R_x-R_y} \right]- \frac 12 S_\varepsilonta\left[(\partial_{xx}iv_x \vc{u}_x-\partial_{xx}iv_y \vc{u}_y)(R_x+R_y)\right]\Big\}\\
       &q_{l}quad\times S_\varepsilonta[R_x-R_y] \lr{S_\varepsilonta[R_x-R_y] ^2+\bar{\partialelta}talta}^{\frac{\beta-2}{2}}.
}
This equation can now be multiplied by $K_h(x-y)\lr{w_x+w_y}$, where $K_h$ is the same as in Section \ref{ssec:prelim}, and $w_x=w(t,x)$ satisfies  the equation
\varepsilonq{\label{def_w_ep}
\left\{\begin{array}{l}
\partialt w+\vc{u}\overlineot\nabla w+ \theta {\cal D}_k w=0,\\
w(0,x)=1,
\varepsilonnd{array}\right.
}
with
\varepsilonq{\label{def_D_ep}
{\cal D}_k=M|\nabla\vc{u}|+|\partial_{xx}iv\vc{u}|+({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+}+ \{({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+}\}.
}
The choice of ${\cal D}_k$ is to accommodate the extra terms appearing in the main compactness estimate, and will be explained later,
the same applies to the choice of the constant $\theta$. The most important observation at this level is that for any $k$, 
fixed ${\cal D}_k\in L^\infty((0,T)\times\mathbb{T}^d)$ and that due to regularity of   $\vc{u}=\vc{u}_k$, see \varepsilonqref{est_kfixed}, Equation \varepsilonqref{def_w_ep} 
has a unique distributional solution, see for example Lemma 6.10 and Lemma 6.12 from \cite{BJ}.

Therefore, letting  $\bar{\partialelta}talta\to 0$, $\beta\to 1$, and $\varepsilonta\to 0$ in {$\int_{\mathbb{T}^{2d}}$\varepsilonqref{eqR_ren} ${\rm d} {x}$} and using the dominated convergence theorem, we recover \varepsilonq{\label{eqRdiff}
&\partial_{xx}t\iintO{ K_h(x-y)|R_x-R_y|\lr{w_x+w_y}}  \\
&=\iintO{ \nabla K_h(x-y)(\vc{u}_x-\vc{u}_y)|R_x-R_y| (w_x+w_y)}\\
&q_{l}uad - \iintO{K_h(x-y)(\partial_{xx}iv_x \vc{u}_x-\partial_{xx}iv_y \vc{u}_y)R_xs_Rw_x}\\
&q_{l}uad+2 \iintO{K_h(x-y)|R_x-R_y| \lr{\partialt w_x+\vc{u}_x\overlineot\nabla w_x+\partial_{xx}iv_x\vc{u}_x w_x}}.
}
To obtain the equation for $|Z_x-Z_y|$ we proceed in a similar manner. As a consequence, the analogue of \varepsilonqref{Rt2} we can be now written as
 \varepsilonq{\label{Sep}
&\partial_{xx}t S(t)= \iintO{ \nabla K_h(x-y)(\vc{u}_x-\vc{u}_y)O_{x-y}  (w_x+w_y)}\\
&= \iintO{ \nabla K_h(x-y)(\vc{u}_x-\vc{u}_y)O_{x-y} (w_x+w_y)}\\
&q_{l}uad - \iintO{K_h(x-y)(\partial_{xx}iv_x \vc{u}_x-\partial_{xx}iv_y \vc{u}_y)[R_xs_R+Z_x s_Z]w_x}\\
&q_{l}uad+2(1-\int_{\mathbb{R}}amma)\iintO{K_h(x-y)\left[\frac{ (Z_y-R_y)Z_y\partial_{xx}iv_y\vc{u}_y}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}-
\frac{ (Z_x-R_x)Z_x\partial_{xx}iv_x \vc{u}_x}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x} 
\right] s_Zw_x}\\
&q_{l}uad+2 \iintO{K_h(x-y)O_{x-y} \lr{\partialt w_x+\vc{u}_x\overlineot\nabla w_x+\partial_{xx}iv_x\vc{u}_x w_x}}\\
&=\sum_{i=1}^4 I_i,
}
where $O_{x-y}=|R_x-R_y|+ |Z_x-Z_y|$.
For abbreviation we will only comment on the changes due to the presence of new approximation terms and truncations. 
   The changes due to truncation in the momentum equation apply to terms $I_2$, and $I_3$. For the first one we use {\varepsilonqref{S-appr4}  to write}
\varepsilonq{
 &-\lr{ \partial_{xx}iv_x\vc{u}_x-\partial_{xx}iv_y\vc{u}_y}[R_xs_R+Z_x s_Z]\\
& =- \lr{ ({\cal T}_k(Z_x))^{\int_{\mathbb{R}}amma^+}-({\cal T}_k(Z_y))^{\int_{\mathbb{R}}amma^+}}R_xs_R
 - \lr{ ({\cal T}_k(Z_x))^{\int_{\mathbb{R}}amma^+}-({\cal T}_k(Z_y))^{\int_{\mathbb{R}}amma^+}}Z_x s_Z,
 }
 that gives the following contribution to the l.h.s. of \varepsilonqref{Sep} 
\varepsilonq{\iintO{K_h(x-y)\left|({\cal T}_k(Z_x))^{\int_{\mathbb{R}}amma^+}-({\cal T}_k(Z_y))^{\int_{\mathbb{R}}amma^+}\right|(-\alpha_x)Z_xw_x}.
}
For $I_3$ we write as in \varepsilonqref{re1}
\varepsilonq{\label{re1ap}
&(1-\int_{\mathbb{R}}amma)\left[\frac{ (Z_y-R_y)Z_y\partial_{xx}iv_y\vc{u}_y}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}-
\frac{ (Z_x-R_x)Z_x\partial_{xx}iv_x \vc{u}_x}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x} 
\right]s_Z\\
&=(1-\int_{\mathbb{R}}amma)\left[\frac{ (Z_y-R_y)Z_y({\cal T}_k(Z_y))^{\int_{\mathbb{R}}amma^+}}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}-
\frac{ (Z_x-R_x)Z_x({\cal T}_k(Z_x))^{\int_{\mathbb{R}}amma^+}}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x} 
\right]s_Z\\
&(1-\int_{\mathbb{R}}amma)\left[ \frac{ (Z_x-R_x)Z_x}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}-
\frac{(Z_y-R_y)Z_y}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y} 
\right] s_Z \{({\cal T}_k(Z))^{\int_{\mathbb{R}}amma^+}\}\\
&=(1-\int_{\mathbb{R}}amma)\frac{(Z_y-R_y)}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}\lr{Z_y({\cal T}_k(Z_y))^{\int_{\mathbb{R}}amma^+}-Z_x({\cal T}_k(Z_x))^{\int_{\mathbb{R}}amma^+}} s_Z\\
&+(\int_{\mathbb{R}}amma-1)
\left[
\frac{ (Z_x-R_x)}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}-
\frac{(Z_y-R_y)}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}
\right] Z_x({\cal T}_k(Z_x))^{\int_{\mathbb{R}}amma^+} s_Z\\
&-(\int_{\mathbb{R}}amma-1)\frac{(Z_y-R_y)}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}\lr{Z_x-Z_y} s_Z \{({\cal T}_k(Z_x))^{\int_{\mathbb{R}}amma^+}\} \\
&-(\int_{\mathbb{R}}amma-1)
\left[
\frac{ (Z_x-R_x)}{ \int_{\mathbb{R}}amma(Z_x-R_x)+R_x}-
\frac{(Z_y-R_y)}{ \int_{\mathbb{R}}amma(Z_y-R_y)+R_y}
\right] Z_x\{({\cal T}_k(Z_x))^{\int_{\mathbb{R}}amma^+} \}s_Z.
}
The first term has a good sign, while the rest of terms can be treated exactly as in the Section \ref{Sect:3}, using in particular the  definition of the weight \varepsilonqref{def_w_ep} and \varepsilonqref{def_D_ep} with $\theta$ sufficiently large  to absorb all the unwanted terms  by  $I_4$ in \varepsilonqref{Sep}.

The above considerations allow us to deduce 
\varepsilonqh{\label{without_w_ep}
\limsup_k\lr{\frac{1}{|\log h_0|}\iintO{ {\mathcal K}_{h_0}(x-y)\Big[ \left|(Z_x)_k-(Z_y)_k\right|+ |(R_x)_k-(R_y)_k|\Big] }}\to 0,
}
by noticing that the thesis of Proposition \ref{prop:w} stay in force also for $(X,\vc{u})$ satisfying the continuity equation with a friction term.
Indeed, since ${\cal D}_k\in L^2((0,T)\times\Omega)$ independently of $k$ (see \varepsilonqref{Z!!!} and \varepsilonqref{v!!!}),  the crucial estimate \varepsilonqref{boundlogw} 
holds for both sequences $\{R_k\}_{k=1}^\infty$, and $\{Z_k\}_{k=1}^\infty$ independently of $k$. This observation finishes the proof of Theorem \ref{Thm:main}. $\Box$
  With this at hand, using Remark \ref{Rmk_1} to define $\alpha^+$ and $\alpha^-$, we find a weak solution to our original system \varepsilonqref{spec}. 
The proof of Theorem \ref{Mainbifluid} is therefore complete. $\Box$

 

{\bf Acknowledgments.}  The author D.B. is partly supported by the 
ANR- 13-BS01- 0003-01 project DYFICOLTI, by the ANR-16-CE06-0011-02 
FRAISE and by the project  TelluS-INSMI-MI (INSU) CNRS. The  authors P.B.M. and E.Z. have been partly supported by National Science Centre grant 2014/14/M/ST1/00108 (Harmonia).
{The authors want to thank P.E. Jabin and the anonymous  referee for their valuable remarks about the paper.}

\begin{thebibliography}{99}

\bibitem{AmZl}
\sc A.A. Amosov, A.A. Zlotnik.
\rm On the error of quasi-averaging of the equations of motion of a viscous
barotoropic medium with rapidly oscillating data. 
\it Comp. Maths Phys, \rm 36 (10), 1415--1428, 1996.

\bibitem{Belgacem} 
\sc F. Ben Belgacem, P.-E. Jabin, 
\rm Compactness for nonlinear continuity equations, 
\it J. Funct. Anal. \rm 264 (1), 139--168, 2013. 

\bibitem{BrDeGhGr}
\sc D. Bresch, B. Desjardins, J.--M. Ghidaglia, E. Grenier. 
\rm On Global Weak Solu- tions to a Generic Two-Fluid Model. 
\it Arch. Rational Mech. Anal. \rm 196 (2), 599--629, 2009.

\bibitem{BrDeGhGrHi}
\sc D. Bresch, B. Desjardins, J.--M. Ghidaglia, E. Grenier, M. Hilliairet.
\rm Multifluid models including compressible fluids.
\rm Handbook of Mathematical Analysis in Mechanics of Viscous Fluids, Eds Y. Giga and A. Novotn\'y, 1--52,  2017.

\bibitem{BrHi}
\sc D. Bresch, M. Hillairet.
\rm  Note on the derivation of multicomponent flow systems. 
\it Proc. AMS, \rm 143, 3429--3443, 2015. 

\bibitem{BrHi2}
\sc D. Bresch, M. Hillairet.
\rm A compressible multifluid system with new physical relaxation terms. 
\rm To appear in \it Ann. Ecole Normale. Sup. \rm  2017.

\bibitem{BrHu}
\sc D. Bresch, X. Huang. 
\rm A multi-fluid compressible system as the limit of weak solutions of the isentropic compressible Navier--Stokes equations.
\it  Arch. Ration. Mech. Anal., \rm 201(2), 647--680, 2011.

\bibitem{BrHuLi}
\sc D. Bresch, X. Huang, J. Li.
\rm A Global Weak Solution to a One-Dimensional Non--Conservative Viscous Compressible Two-Phase System.
\it  Comm. Math. Phys. \rm 309 (3), 737--755, 2012.

\bibitem{BJ} 
\sc D. Bresch, P.-E. Jabin.
\rm Global Existence of Weak Solutions for Compresssible Navier--Stokes Equations: Thermodynamically unstable
pressure and anisotropic viscous stress tensor.
\it  Ann. of Math., 
\rm 188 (2), 577--684, 2018.

\bibitem{BJ_short} 
\sc D. Bresch, P.-E. Jabin. 
\rm Global weak solutions of PDEs for compressible media: A compactness criterion to cover new physical situations.
\it Shocks, Singularities and Oscillations in Nonlinear Optics and Fluid Mechanics, 
\rm Springer INdAM Vol. 17, 33--54, 2017.

\bibitem{BrJa}
\sc D. Bresch, P.--E. Jabin.
\rm Quantitative regularity estimates for advective equation 
       with anelastic degenerate constraint.
\it Proc. of the ICM Rio de Janeiro, \rm 2018.

\bibitem{CoCrSp}
\sc M. Colombo, G. Crippa, S. Spirito.
\rm Renormalized solutions to the continuity equation with an integrable damping term.
\it Calculs Variations and PDEs, \rm 54, 2, 1831--1845, (2015).

\bibitem{CrDe}
\sc G. Crippa, C. De Lellis.
\rm Estimates and regularity results for the DiPerna-Lions flows.
\it J. F\"ur Die Reine Und Angewandte Mathematik.
\rm 616, 15--46, 2008.

\bibitem{DeMiZa}
\sc P. Degond, P. Minakowski, E. Zatorska.
\rm Transport of congestion in two-phase compressible/incompressible flows with. 
\it  Nonlinear Anal. Real World Appl., \rm 42: 485 -- 510, 2018.
    
 \bibitem{DeMiNaZa}
\sc P. Degond, P. Minakowski, L. Navoret, E. Zatorska.   
\rm    Finite Volume approximations of the Euler system with variable congestion.
\it Computers $\&$ Fluids, \rm 169 (30), 23--39, 2018.

\bibitem{De}
 \sc L. Desvillettes. 
 \rm Some Aspects of the Modeling at Different Scales of Multiphase Flows. 
 \it Computer Methods in Applied Mechanics and Engineering, \rm 199, 1265--1267, 2010.
 
 \bibitem{DiDuGh}
 \sc F. Dias, D. Dutykh, and J.-M. Ghidaglia. 
 \rm A two-fluid model for violent aerated flows.
 \it Computers \& Fluids, \rm 39 (2), 283--293, 2010.

\bibitem{DL} 
\sc R.J. DiPerna, P.L. Lions.
\rm Ordinary differential equations, transport theory and Sobolev spaces. 
\it Invent. Math. \rm 98,  511--547, 1989.

\bibitem{DrPa} 
\sc D. Drew, S.L. Passman. 
\it Theory of multicomponent fluids. 
\rm Applied Math Sciences, 1999.

\bibitem{Ev}
\sc S. Evje.
\rm An integrative multiphase model for cancer cell migration under influence of
physical cues from the tumor microenvrionment.
\it Chemical Engineering Science, \rm 204--259, 2017.

\bibitem{Evje08}
\sc S. Evje, K.H. Karlsen.
\rm Global existence of weak solutions for a viscous two-phase model. 
\it J. Diff. Equations, \rm 245 (9), 2660--2703, 2008.

\bibitem{Evje15}
\sc S. Evje, H. Wen. 
\rm Analysis of a compressible two-fluid stokes system with constant viscosity.
 \it J. Math. Fluid Mech. \rm17 (3), 423--436, 2015.

\bibitem{Evje16a}
\sc S. Evje, H. Wen.
 \rm Stability of a compressible two-fluid hyperbolic--elliptic system arising in fluid mechanics.
 \it Nonlinear Anal.: Real World Applications, \rm 31,  610--629, 2016.
 
 \bibitem{Evje16}
\sc S. Evje, H. Wen, L. Yao.
\rm Global solutions to a one-dimensional non-conservative two-phase model.
\it Discrete and Continuous Dyn. Systems.
\rm 36 (4), 1927--1955, 2016.

\bibitem{EF2001}
\sc E.  Feireisl.
\rm On compactness of solutions to the compressible isentropic {N}avier-{S}tokes equations when the density is not square integrable.
 \it Comment. Math. Univ. Carolin. \rm 42~(1), 83--98, 2001.

\bibitem{FeKlNoZa}
\sc E. Feireisl, R. Klein, A. Novotn{{\'y}}, E. Zatorska.
\rm  On singular limits arising in the scale analysis of stratified fluid flows.
 \it Math. Models Methods Appl. Sci. \rm 26~(3), 419--443, 2016.
 
 \bibitem{Hao12}
\sc C.C Hao, H.L. Li. 
 \rm Well-posedness for a multidimensional viscous liquid-gas flow model, 
 \it SIAM J. Math. Anal., \rm 44 (3), 1304--1332, 2012.

\bibitem{IsHi}
\sc M. Ishii, T. Hibiki. 
\it Thermo-fluid dynamics of two-phase flow. 
\rm Springer, 2006.

\bibitem{Jabin2010}
 \sc P.-E. Jabin.
\rm Differential Equations with singular fields. 
\it Journal de Math\'ematiques Pures et Appliqu\'ees, 
\rm 94 (6), 597--621, 2010.

\bibitem{KoTa}
\sc H. Kozono, Y. Taniuchi.
\rm Limiting case of the Sobolev inequality in BMO, with application to the Euler equations.
\it Comm. Math. Phys., \rm 214, 191--200, 2000.
 

\bibitem{Lions1}
 \sc P.-L. Lions.
\it  Mathematical Topics in Fluid Mechanics, Vol 1: Incompressible Models.
\rm  New York: Oxford University Press, 1996.

\bibitem{Lions2} 
\sc P.-L. Lions. 
\it Mathematical Topics in Fluid Mechanics, Vol 2: Compressible Models. 
\rm New York: Oxford University Press, 1998.

\bibitem{Ogawa}
\sc T. Ogawa.
\rm Sharp Sobolev Inequality of Logarithmic Type and the Limiting Regularity Condition to the Harmonic Heat Flow.
\it  SIAM J. Math. Anal., \rm 34 (6), 1318--1330, 2003.

\bibitem{MaMiMuNoPoZa}
\sc D. Maltese, M. Michalek, P.B. Mucha, A. Novotny, M. Pokorny, E. Zatorska.
\rm Existence of weak solutions for compressible Navier-Stokes equations with entropy transport.
\it  J. Differential Equations, \rm 261 (8), 4448--4485, 2016.

\bibitem{Michalek}
\sc M. Mich\'alek.
\rm Stability result for Navier-Stokes equations with entropy transport.
\it  J. Math. Fluid Mech., \rm  17 (2), 279--285, 2015.

\bibitem{Mucha2003}
\sc P. B. Mucha.
\rm The Cauchy problem for the compressible Navier-Stokes equations in the Lp-framework. 
\it Nonlinear Anal., \rm 52 (4), 1379--1392, 2003.

\bibitem{MuchaZaj2002}
\sc P. B. Mucha, W. Zajaczkowski.
\rm On a $L_p$-estimate for the linearized compressible Navier-Stokes equations with the Dirichlet boundary conditions. 
\it J. Differential Equations, \rm 186 (2), 377--393, 2002.

\bibitem{MuchaZaj2006}
\sc P. B. Mucha, W. Zajaczkowski.
\rm Global existence of solutions of the Dirichlet problem for the compressible Navier-Stokes equations. 
\it ZAMM Z. Angew. Math. Mech., \rm 84 (6), 417--424, 2004.

\bibitem{Se}
\sc D. Serre. 
\rm Asymptotics of homogeneous oscillations in a compressible viscous fluid.
\it  Bul. Soc. Bras. Mat, \rm 32, 535--442, 2001.

\bibitem{Stein} 
\sc E. M. Stein. 
\rm Maximal functions. I. Spherical means.
\it Proc. Nat. Acad. Sci. USA \rm  73 (7), 2174--2175,  1976.

\bibitem{VauchZator}
\sc N. Vauchelet, E. Zatorska.
\rm Incompressible limit of the Navier--Stokes model with a growth term.
\it Nonlinear Analysis, \rm 163, 34--59, 2017.

\bibitem{VaWeYu}
\sc A. Vasseur, H. Wen, C. Yu.
\rm Global weak solution to the viscous two-phase model with finite energy.
\it Preprint, \rm arXiv:1704.07354.

\varepsilonnd{thebibliography}

\varepsilonnd{document}