# Weak solutions for the Stokes system for compressible non-Newtonian fluids
with unbounded divergence

Milan Pokorný111Charles University, Faculty of Mathematics and Physics,
Sokolovská 83, 186 75 Praha 8, Czech Republic, e-mail address:
<EMAIL_ADDRESS>& Maja Szlenk222University of Warsaw, Faculty of
Mathematics, Informatics and Mechanics, Banacha 2, 02-097 Warsaw, Poland,
e-mail address<EMAIL_ADDRESS>
###### Abstract

We investigate the existence of weak solutions to a certain system of partial
differential equations, modelling the behaviour of a compressible non-
Newtonian fluid for small Reynolds number. We construct the weak solutions
despite the lack of the $L^{\infty}$ estimate on the divergence of the
velocity field. The result was obtained by combining the regularity theory for
singular operators with a certain logarithmic integral inequality for $BMO$
functions, which allowed us to adjust the method from [13] to more relaxed
conditions on the velocity.

Keywords: compressible Stokes system, weak solutions, non-Newtonian fluids,
power law fluids

##  1 Introduction

Our aim is to investigate the existence of weak solutions to equations,
modelling a special case of compressible, non-Newtonian fluid. In the most
general setting, the motion of such a fluid without the presence of the
external forces is described by the system of partial differential equations

$\displaystyle\varrho_{t}+\mathrm{div}\,(\varrho u)=0,$ (1.1)
$\displaystyle(\varrho u)_{t}+\mathrm{div}\,(\varrho u\otimes
u)-\mathrm{div}\,\mathbb{S}=0,$

where $\varrho$ is the density, $u$ is a velocity vector and $\mathbb{S}$ is
the stress tensor; we assume that it is given by

$\mathbb{S}(\mathbb{D}u,\varrho)=\mu\mathbb{D}u+\big{(}\lambda\mathrm{div}\,u-p(\varrho)\big{)}\mathbb{I},$

where $\mu>0$ and $\lambda$ are the viscosity coefficients,
$\mathbb{D}=\frac{1}{2}(\nabla+\nabla^{T})$ is the symmetric gradient and
$p(\varrho)$ is the pressure. In the case of constant viscosity (i.e., the
resulting system is called the compressible Navier–Stokes equations)
$d\lambda+\mu\geq 0$, where $d$ is the space dimension.

We will focus on the case where the Reynolds number
$\mathrm{Re}\sim\frac{\varrho|u|}{\mu}$ is small. As in this situation the
advective forces are small compared to the viscous ones, we can approximate
system (1.1) by the compressible Stokes-like system

$\displaystyle\varrho_{t}+\mathrm{div}\,(\varrho u)=0,$ (1.2)
$\displaystyle-\mathrm{div}\,\mathbb{S}=0.$

Our aim is to obtain weak solutions to a special case of system (1.2). We
assume that the shear viscosity $\mu$ is in the form

$\mu=\mu_{0}(|\mathbb{D}u|)+2\mu_{1},\quad\mu_{1}>0\;\;\text{constant}$ (1.3)

and the bulk viscosity $\lambda=\lambda(|\mathrm{div}\,u|)$, where

$0\leq\mu_{0}(z),\lambda(z)\leq\frac{C}{z},\quad z>0.$ (1.4)

Furthermore, the functions $\mathbb{R}^{d\times d}\ni B\mapsto\mu_{0}(|B|)B$
and $\mathbb{R}\ni s\mapsto\lambda(|s|)s$ are assumed to be monotone.

For the pressure we assume the barotropic case with
$p(\varrho)=\varrho^{\gamma}$ for $\gamma\geq 1$. (Indeed, we could replace
this precise form just by asymptotic growth conditions, similarly as in [13],
i.e. $p(0)=0$, $p^{\prime}(z)>0$ for $z>0$ and
$\lim_{z\to\infty}\frac{p(z)}{z^{\gamma}}\in(0,\infty)$ for some $\gamma\geq
1$, but we skip such unnecessary complications). For simplicity we consider
the space-periodic boundary conditions, hence
$u\colon[0,T]\times\mathbb{T}^{d}\to\mathbb{R}^{d}$ and
$\varrho\colon[0,T]\times\mathbb{T}^{d}\to\mathbb{R}$, where $\mathbb{T}^{d}$
is the $d$-dimensional torus. In conclusion, the analysed system of equations
yield

$\displaystyle\varrho_{t}+\mathrm{div}\,(\varrho u)$ $\displaystyle=0,$ (1.5)
$\displaystyle-\mathrm{div}\,(\mu_{0}(|\mathbb{D}u|)\mathbb{D}u)-\mu_{1}\Delta
u-\nabla((\mu_{1}+\lambda(\mathrm{div}\,u))\mathrm{div}\,u)+\nabla\varrho^{\gamma}$
$\displaystyle=0,$

with the initial condition

$\varrho_{|_{t=0}}=\varrho_{0}\in L^{\infty}(\mathbb{T}^{d})$ (1.6)

and

$\int_{\mathbb{T}^{d}}u(t,x)\,\mathrm{d}x=0\quad\forall_{t>0}.$

Our system describes a fluid belonging to a class of power-law fluids. They
are characterized by the behavior of the shear viscosity, which satisfies the
relation

$\mu\sim|\mathbb{D}u|^{r-2}$ (1.7)

for some exponent $r\geq 1$. Typically, it is assumed that
$\mu=\mu_{0}|\mathbb{D}u|^{r-2}$ or $\mu=\mu_{0}(a+|\mathbb{D}u|)^{r-2}$,
$a>0$, to ensure that the viscosity is strictly positive and does not have
singularities. For $r=2$, the fluid becomes Newtonian, whereas it is shear-
thinning for $r<2$ and shear-thickening for $r>2$. The power-law fluids are
used in many fields, for example glaciology [19, 14] and to analyze the
dynamics in the Earth’s Mantle [24] or blood flow [7, 22]. For more
information we refer the reader, e.g., to [4]. Our situation corresponds
specifically to a Herschel-Bulkley fluid, where the shear viscosity is in the
form

$\mu=\left\\{\begin{aligned} \mu_{0},&\quad|\mathbb{D}u|<\delta,\\\
\frac{\tau_{0}}{|\mathbb{D}u|}+k|\mathbb{D}u|^{n-1},&\quad|\mathbb{D}u|\geq\delta\end{aligned}\right.$

for some $n\geq 1$ and the parameters $\mu_{0},\tau_{0},k$ are chosen in such
way that $\mu$ remains continuous. Fluids of this type were thoroughly
analysed in the incompressible case, and have many industrial applications,
see e.g. [3, 8, 11].

The mathematical theory concerning weak solutions to systems describing
incompressible non-Newtonian fluids have been thoroughly developed in the
past. There is a large number of papers dealing with several aspects of these
problems. As it turns out, the existence and regularity of solutions to
incompressible Navier–Stokes equations with the power-law relation (1.7) for
viscosity depends on the value of $r$. For $r>\frac{2d}{d+2}$ the existence of
weak solutions for the problem

$\displaystyle\mathrm{div}\,u$ $\displaystyle=0,$ (1.8) $\displaystyle
u_{t}+\mathrm{div}\,(u\otimes u)-\mathrm{div}\,\mathbb{S}$ $\displaystyle=0$

with Dirichlet boundary conditions was shown for the first time in [10]; its
uniqueness is known for $r\geq\frac{3d+2}{d+2}$, see [5]. As a matter of fact,
the problem for $r<\frac{2d}{d+2}$ is ill-posed, see [6]. However, existence
of more general, dissipative solutions can be shown also in this case, see
[1].

Contrary to the incompressible case, the literature on the compressible non-
Newtonian fluids is very limited. In [17, 18] Mamontov proved the existence of
weak solutions to the system with linear pressure term and in the framework of
Orlicz spaces with exponential growth, see also [2] for further properties of
these solutions. The results for more general form of the pressure were
obtained in [13], where the authors considered the system (1.1) with $\mu$ of
the form (1.7) and a special form of $\lambda$, which provided the
$L^{\infty}$ bound on $\mathrm{div}\,u$. Using the classical Lions & Feireisl
method [12, 16], they proved the existence of weak solutions for the same
range of $r$ as the uniqueness and regularity theory is developed for
incompressible fluids ($r\geq\frac{11}{5}$ in three dimensions). The
additional bound on the divergence was crucial to obtain the strong
convergence of the density in the final limit passage.

###  1.1 Main result and structure of the paper

Since the definition of the weak solution for our system (1.2) is
straightforward, we may directly formulate our main result.

######  Theorem 1.1.

Let $\varrho_{0}\in L^{\infty}(\mathbb{T}^{d})$, $\gamma\geq 1$ and let (1.4)
hold. Then for any $T>0$ there exists a weak solution to the system (1.5),
satisfying

$\|\nabla
u\|_{L^{2}((0,T)\times\mathbb{T}^{d})}+\|\mathrm{div}\,u\|_{L^{\infty}(0,T;L^{p})}+\|\varrho\|_{L^{\infty}(0,T;L^{p})}\leq
C(p,T),$

for any $1\leq p<\infty$, where $C$ approaches $\infty$ if $p$ or $T$ does so.

The proof of the theorem uses the technique from [13]. However, due to the
absence of the convective term we are able to obtain the result for relaxed
assumptions on $\mathrm{div}\,u$. In particular, we do not have from the very
beginning that $\mathrm{div}\,u$ is bounded. Instead of it, we obtain the
$BMO$ regularity in space of the term

$(\mu_{1}+\lambda)\mathrm{div}\,u-\varrho^{\gamma}.$

This allows us to replicate the main step in the limit passage by using
certain integral inequality from [23].

###  1.2 Preliminaries

In the paper, we use standard notation for the Lebesgue and Sobolev space as
well as the corresponding norms. For the notational simplicity, we omit the
subscript while integrating over the torus, namely

$\int\,\mathrm{d}x\ :=\ \int_{\mathbb{T}^{d}}\,\mathrm{d}x.$

By $\\{\cdot\\}_{Q}$ we denote the mean integral over $Q$, while in the case
of the whole torus we again omit the subscript. As the results do not depend
on the values of $\mu_{1}$, for simplicity we set $\mu_{1}=1$.

A few other, slightly nonstandard spaces are defined below.

#### Functions of bounded mean oscillation.

A function $f\in L^{1}(\Omega)$ belongs to space of bounded mean oscillation
$BMO(\Omega)$ iff

$\|f\|_{BMO}=\sup_{Q\subset\Omega}\frac{1}{|Q|}\int_{Q}|f-\\{f\\}_{Q}|\,\mathrm{d}x<\infty,$

where the supremum is taken over all cubes in $\Omega$.

Note that $\|\cdot\|_{BMO}$ is not a norm, as $\|f\|_{BMO}=0$ for $f$
constant. However, we can endow the space $BMO(\Omega)$ with the norm

$\|\cdot\|_{L^{1}}+\|\cdot\|_{BMO}$

and then it becomes the Banach space. The important result concerning $BMO$
functions in terms of our work is the following inequality:

$\displaystyle\left|\int_{\mathbb{T}^{d}}fg\,\mathrm{d}x\right|\leq$
$\displaystyle C\|f\|_{BMO}\|g\|_{L^{1}}$ (1.9)
$\displaystyle\times\big{(}|\ln\|g\|_{L^{1}}|+\ln(e+\|g\|_{L^{q}})+(1+|\ln\|g\|_{L^{1}}|)\|g\|_{L^{q}}^{\frac{q-2}{2}}\big{)}$

for $g\in L^{q}$, $q>2$. The inequality was first obtained for $g\in
L^{\infty}$ in [20] and then adjusted in [23] for a wider range of functions.
In [9] the similar inequality was obtained for $f$ in the Orlicz space
$L^{\exp}$.

The rest of the article is devoted to prove Theorem 1.1. First, in Section 2
we derive the a priori estimates, in particular the crucial $BMO$ estimate for
the quantity $(2+\lambda(\mathrm{div}\,u))\mathrm{div}\,u-\varrho^{\gamma}$.
Then, in Section 3 we prove the existence of solutions to an approximated
system with the regularized continuity and momentum equations. In Section 4 we
finish the proof by passing to the limit in the weak formulation of the
approximate system and in consequence we obtain the weak solution of the
original one.

##  2 A priori estimates

######  Lemma 2.1.

Under the assumptions of Theorem 1.5, if the solution to (1.5) is sufficiently
smooth, it satisfies

$\|\nabla
u\|_{L^{2}((0,T)\times\mathbb{T}^{d})}+\|\varrho\|_{L^{\infty}(0,T;L^{\gamma})}\leq
C$

and

$\|(2+\lambda(\mathrm{div}\,u))\mathrm{div}\,u-\varrho^{\gamma}\|_{L^{\infty}(0,T;BMO)}\leq
C$

for $C$ depending only on $\|\varrho_{0}\|_{\infty}$. Furthermore,

$\|\mathrm{div}\,u\|_{L^{\infty}(0,T;L^{p})}+\|\varrho\|_{L^{\infty}(0,T;L^{p})}\leq
C(p,T)$

for any $p<\infty$, where $C(p)\to\infty$ as $p\to\infty$ or $T\to\infty$.

###### Proof.

Multiplying the second equation of (1.5) by $u$ and integrating over the
torus, we obtain (if $\gamma=1$, the last integral is replaced by
$\int\varrho\ln\varrho\,\mathrm{d}x$)

$\int\mu_{0}(|\mathbb{D}u|)|\mathbb{D}u|^{2}\,\mathrm{d}x+\int|\nabla
u|^{2}\,\mathrm{d}x+\int(\mathrm{div}\,u)^{2}\,\mathrm{d}x+\int\lambda(\mathrm{div}\,u)(\mathrm{div}\,u)^{2}\,\mathrm{d}x+\frac{1}{\gamma-1}\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int\varrho^{\gamma}\,\mathrm{d}x=0.$

Integrating the above equality from $0$ to $T$, we obtain the first desired
estimate.

To obtain the $L^{p}$ estimate of the density, we use as test function in
(1.2)2 the function
$\psi=-\Delta^{-1}\nabla\left(\varrho^{\theta}-\\{\varrho^{\theta}\\}\right)$
for some $\theta>1$. We have

$\mathrm{div}\,\psi=\varrho^{\theta}-\\{\varrho^{\theta}\\}$

and

$\|\nabla\psi\|_{L^{r}((0,T)\times\mathbb{T}^{d})}\leq
C(r)\|\varrho^{\gamma}\|_{L^{r}((0,T)\times\mathbb{T}^{d})}\;\text{for
any}\;1<r<\infty.$

Note that $C(r)\to\infty$ if $r\to 1^{+}$ or $r\to\infty$. Moreover,

$-\int\Delta u\cdot\psi\,\mathrm{d}x=\int
u\cdot\nabla\left(\varrho^{\theta}-\\{\varrho^{\theta}\\}\right)\,\mathrm{d}x=-\int\varrho^{\theta}\mathrm{div}\,u\,\mathrm{d}x.$

Then

$\displaystyle\int\varrho^{\gamma+\theta}\,\mathrm{d}x-2\int\varrho^{\theta}\mathrm{div}\,u\,\mathrm{d}x$
$\displaystyle=-\int\mu_{0}(|\mathbb{D}u|)\mathbb{D}u:\nabla\psi\,\mathrm{d}x-\int\varrho^{\theta}\lambda(\mathrm{div}\,u)\mathrm{div}\,u\,\mathrm{d}x$
$\displaystyle+\frac{1}{|\mathbb{T}^{d}|}\int\varrho^{\gamma}\,\mathrm{d}x\int\varrho^{\theta}\,\mathrm{d}x.$

As

$\int\varrho^{\theta}\mathrm{div}\,u\,\mathrm{d}x\,\mathrm{d}t=\frac{-1}{\theta-1}\frac{{\rm}d}{{\rm
d}t}\int\varrho^{\theta}\,\mathrm{d}x,$

we obtain

$\frac{2}{\theta-1}\frac{{\rm}d}{{\rm
d}t}\int\varrho^{\theta}\,\mathrm{d}x+\int\varrho^{\gamma+\theta}\,\mathrm{d}x\leq
C\left(\int\varrho^{\theta}\,\mathrm{d}x+\left(\int\varrho^{(1+\delta)\theta}\,\mathrm{d}x\right)^{\frac{1}{1+\delta}}\right);$

whence, for a suitably chosen $\delta$

$\|\varrho\|_{L^{\infty}(0,T;L^{p}}\leq C(p,T,\varrho_{0}).$

The last estimate comes from the Calderón–Zygmund estimates. By taking the
divergence of the momentum equation, we get

$\Delta((2+\lambda(\mathrm{div}\,u))\mathrm{div}\,u-\varrho^{\gamma})=-\mathrm{div}\,\mathrm{div}\,\left(\mu_{0}(|\mathbb{D}u|)\mathbb{D}u\right).$

Therefore in consequence

$\displaystyle(2+\lambda(\mathrm{div}\,u(t,x)))\mathrm{div}\,u(t,x)-\varrho^{\gamma}(t,x)$
$\displaystyle=-\int\mathrm{div}\,\mathrm{div}\,\left(\mu_{0}(|\mathbb{D}u|)\mathbb{D}u\right)\bar{K}(x-y)\,\mathrm{d}y$
$\displaystyle=-\int\mathrm{div}\,\left(\mu_{0}(|\mathbb{D}u|)\mathbb{D}u\right)\cdot\nabla\bar{K}(x-y)\,\mathrm{d}y$
$\displaystyle=-\int\mu_{0}(|\mathbb{D}u|)\mathbb{D}u:\nabla^{2}\bar{K}(x-y)\,\mathrm{d}y,$

where $\Delta\bar{K}(x)=\delta_{x}$ in $\mathbb{T}^{d}$. Note that we can
derive $\bar{K}$ by periodizing the fundamental solution to the Laplace
equation on the whole space. If $K(x)=\frac{C_{d}}{|x|^{d-2}}$, then we simply
put

$\bar{K}(x)=K(x)+\sum_{k\in\mathbb{Z}^{d}\setminus\\{0\\}}(K(x+k)-K(k)).$

As $\nabla^{2}K$ is a singular kernel, so is $\nabla^{2}\bar{K}$. Therefore
from the Calderón–Zygmund theorem we conclude that

$\|(2+\lambda(\mathrm{div}\,u))\mathrm{div}\,u-\varrho^{\gamma}\|_{L^{\infty}(0,T;BMO)}\leq\left\|\mu_{0}(|\mathbb{D}u|)\mathbb{D}u\right\|_{L^{\infty}((0,T)\times\mathbb{T}^{d})}\leq
C$

for some constant $C$ independent of $T$. Since $\varrho^{\gamma}$ is bounded
in $L^{\infty}(0,T;L^{p})$ and $\lambda(\mathrm{div}\,u))\mathrm{div}\,u$ in
$L^{\infty}((0,T)\times\mathbb{T}^{d})$, we finish the proof. ∎

##  3 Existence of approximate solutions

In this section we construct a unique solution to the system

$\displaystyle\varrho_{t}+\mathrm{div}\,(\varrho u)+\delta\varrho^{\beta}$
$\displaystyle=\delta\Delta\varrho,$ (3.1)
$\displaystyle-\mathrm{div}\,(\mu_{0}(\mathbb{D}u)\mathbb{D}u)-\Delta
u-\nabla(1+\lambda(\mathrm{div}\,u))\mathrm{div}\,u+\nabla\varrho^{\gamma}$
$\displaystyle=-\varepsilon\Delta^{2m}u,\quad\int u(t,x)\,\mathrm{d}x=0$

for a sufficiently small $\delta,\varepsilon>0$, sufficiently large
$m\in\mathbb{N}$ and $\beta\geq\max\\{\gamma+1,4\\}$ being an even integer,
with the initial condition

$\varrho_{|_{t=0}}=\varrho_{0,\delta}\in
C^{\infty}(\mathbb{T}^{d}),\quad\varrho_{0,\delta}>0,\quad\varrho_{0,\delta}\to\varrho_{0}\;\text{in
any}\;L^{p},\;p<\infty.$

To prove the existence of solutions, we will employ the following version of
the Schauder fixed point theorem:

######  Theorem 3.1.

Let $X$ be a Banach space and $\Phi\colon X\to X$ be continuous and compact.
If the set

$\\{x\in X:\quad x=s\Phi(x)\quad\text{for some}\quad s\in[0,1]\\}$

is bounded, then $\Phi$ has a fixed point.

Let us define the map $\Phi\colon C([0,T];L^{2\gamma})\to
C([0,T];L^{2\gamma})$ in the following way:

  1. 1.

For $\tilde{\varrho}\in C([0,T];L^{2\gamma})$, let $u$ be the unique solution
to the equation

$-\mathrm{div}\,(\mu_{0}(|\mathbb{D}u|)\mathbb{D}u)-\Delta
u-\nabla(1+\lambda(\mathrm{div}\,u))\mathrm{div}\,u+\nabla(\tilde{\varrho})^{\gamma}=-\varepsilon\Delta^{2m}u,\quad\int
u(t,x)\,\mathrm{d}x=0.$ (3.2)

  2. 2.

Then, let $\varrho$ be the solution to

$\varrho_{t}+\mathrm{div}\,(\varrho
u)+\delta\varrho^{\beta}=\delta\Delta\varrho,\quad\varrho_{|_{t=0}}=\varrho_{0,\delta}.$
(3.3)

We set $\Phi(\tilde{\varrho}):=\varrho$. It is easy to see that the fixed
point $\varrho$ and the corresponding $u$ solve our problem (3.1).

First, let us show that the operator $\Phi$ is well-defined.

######  Proposition 3.1.

If $\tilde{\varrho}^{\gamma}\in L^{\infty}(0,T;L^{2})$, then there exists a
unique solution $u$ to equation (3.2), satisfying

$\|\nabla
u\|_{L^{\infty}(0,T;L^{2})}+\sqrt{\varepsilon}\|\Delta^{m}u\|_{L^{\infty}(0,T;L^{2})}\leq
C\|\tilde{\varrho}^{\gamma}\|_{L^{\infty}(0,T;L^{2})}.$

In particular, if $m$ is large enough, then

$\|u\|_{L^{\infty}(0,T;W^{1,\infty})}\leq\frac{C}{\sqrt{\varepsilon}}\|\tilde{\varrho}^{\gamma}\|_{L^{\infty}(0,T;L^{2})}.$

###### Proof.

By multiplying the equation by $u$ and integrating over the torus, we get

$\displaystyle\int\mu_{0}(|\mathbb{D}u|)|\mathbb{D}u|^{2}+|\nabla
u|^{2}+(1+\lambda(\mathrm{div}\,u))(\mathrm{div}\,u)^{2}+$
$\displaystyle\varepsilon|\Delta^{m}u|^{2}\;\,\mathrm{d}x=\int\tilde{\varrho}^{\gamma}\mathrm{div}\,u\;\,\mathrm{d}x$
$\displaystyle\leq\eta\int(\mathrm{div}\,u)^{2}\,\mathrm{d}x+\frac{C}{\eta}\|\tilde{\varrho}^{\gamma}\|_{L^{\infty}(0,T;L^{2})}^{2},$

hence picking $\eta$ small enough and taking supremum over time, we get the
desired estimate.

For existence, we consider the functional $I$ defined in
$\overline{H^{2m}}(\mathbb{T}^{d})=\\{v\in H^{2m}(\mathbb{T}^{d}):\int
v\,\mathrm{d}x=0\\}$ given by

$I[v]:=\int\Big{(}F(\nabla v)+\frac{1}{2}|\nabla
v|^{2}+\Lambda(\mathrm{div}\,v)+\frac{\varepsilon}{2}|\Delta^{m}v|^{2}-\tilde{\varrho}^{\gamma}(t,\cdot)\mathrm{div}\,v\Big{)}\,\mathrm{d}x,$

where $F$ satisfies

$\frac{\partial}{\partial b_{i,j}}F(B)=\mu_{0}(|B|)b_{i,j},$

for $B=(b_{i,j})_{i,j}\in\mathbb{R}^{d\times d}$ and $\Lambda$ is such that
$\Lambda^{\prime}(s)=s+\lambda(s)s$. In particular, the assumptions on
$\mu_{0}$ and $\lambda$ imply that $F$ and $\Lambda$ are convex and bounded
from below.

From the definitions of $F$ and $\Lambda$ it follows that any minimizer of $I$
corresponds to a weak solution to (3.2). By the convexity of $F$ and
$\Lambda$, the functional $I$ is convex. Moreover, for certain $C$ and
$\eta<C$,

$I[v]\geq\varepsilon\|\Delta^{m}v\|_{L^{2}}^{2}+C\|\nabla
v\|_{L^{2}}^{2}-\eta\|\nabla
v\|_{L^{2}}^{2}-\frac{C}{\eta}\|\tilde{\varrho}^{\gamma}(t,\cdot)\|_{L^{2}}^{2}\geq
C\|v\|_{H^{2m}}^{2}-C\|\tilde{\varrho}^{\gamma}\|_{L^{\infty}(0,T;L^{2})}^{2},$

and thus $I$ is coercive. Therefore $I$ has at a.e. time level a unique
minimizer $v(t,\cdot)\in\overline{H^{2m}}(\mathbb{T}^{d})$ and in consequence
there exists a unique $u\in L^{\infty}(0,T;H^{2m})$ with zero mean value over
the torus solving (3.2). ∎

Now we use the following classical result for the heat equation:

######  Proposition 3.2.

Let $h\in L^{2}(0,T;L^{q})$ for $1<q<\infty$. Then the solution to

$\partial_{t}\varrho-\varepsilon\Delta\varrho=h,\quad\varrho_{|_{t=0}}=\varrho_{0}$

satisfies the estimate

$\varepsilon^{1/2}\|\varrho\|_{L^{\infty}(0,T;W^{1,q})}+\|\partial_{t}\varrho\|_{L^{2}(0,T;L^{q})}+\varepsilon\|\varrho\|_{L^{2}(0,T;W^{2,q})}\leq
C\left(\varepsilon^{1/2}\|\varrho_{0}\|_{W^{1,q}}+\|h\|_{L^{2}(0,T;L^{q})}\right).$
(3.4)

Moreover, if $h=\mathrm{div}\,w$, $w\in L^{2}(0,T;L^{q})$, then

$\varepsilon^{1/2}\|\varrho\|_{L^{\infty}(0,T;L^{q})}+\varepsilon\|\nabla\varrho\|_{L^{2}(0,T;L^{q})}\leq
C\left(\varepsilon^{1/2}\|\varrho_{0}\|_{L^{q}}+\|w\|_{L^{2}(0,T;L^{q})}\right).$
(3.5)

From the previous proposition, we can also conclude

######  Proposition 3.3.

If $u\in L^{\infty}(0,T;W^{1,\infty})$, then for equation (3.3) there exists a
unique nonnegative solution $\varrho\in L^{\infty}(0,T;W^{1,r})$ with
$\partial_{t}\varrho\in L^{2}(0,T;W^{-1,r})$ for any $r<\infty$.

###### Proof.

We construct the solution $\varrho\in L^{\infty}(0,T;W^{1,2})$ by the Galerkin
approximation. The nonnegativity of solutions is obtained by testing by
negative part of $\varrho$ and is a conclusion of the fact that
$\varrho^{\beta}\geq 0$. Then we improve the regularity by bootstrapping
argument. We skip the details of these steps since they are based on classical
arguments. Next, testing equation (3.3) by $p\varrho^{p-1}$, we have

$\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int\varrho^{p}\,\mathrm{d}x\leq(p-1)\|\mathrm{div}\,u\|_{L^{\infty}}\int\varrho^{p}\,\mathrm{d}x$

and therefore

$\|\varrho\|_{L^{\infty}(0,T;L^{p})}\leq\|\varrho_{0,\varepsilon}\|_{L^{p}}e^{\frac{p-1}{p}\|u\|_{L^{1}(0,T;W^{1,\infty})}}.$

Taking $p\to\infty$ we have $\varrho\in
L^{\infty}((0,T)\times\mathbb{T}^{d})$. In consequence $\varrho u\in
L^{2}(0,T;L^{r})$ for any $r<\infty$ and we can use (3.5) to obtain

$\nabla\varrho\in L^{2}(0,T;L^{r}).$

Employing the fact that $u\in L^{\infty}(0,T;W^{1,\infty})$, we have
$\mathrm{div}\,(\varrho u)\in L^{2}(0,T;L^{r})$ and by (3.4) with
$h=-\mathrm{div}\,u-\delta\varrho^{\beta}$, $\varrho\in
L^{\infty}(0,T;W^{1,r})$ for any $r<\infty$, whereas the estimate for
$\partial_{t}\varrho$ comes directly from the equation (3.3). The uniqueness
follows directly from the estimate

$\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int(\varrho_{1}-\varrho_{2})^{2}\,\mathrm{d}x\leq
C\int(\varrho_{1}-\varrho_{2})\mathrm{div}\,u\,\mathrm{d}x,$

where $\varrho_{1}$ and $\varrho_{2}$ are two possibly distinct solutions to
problem (3.3). ∎

We will now show the properties of $\Phi$, which will allow us to apply
directly the Schauder fixed point theorem (Theorem 3.1).

######  Proposition 3.4.

The operator $\Phi$ is continuous and compact from $C([0,T];L^{2\gamma})$ to
itself. Moreover, the set

$\\{\varrho\in C([0,T];L^{2\gamma}):\varrho=s\Phi(\varrho)\quad\text{for
some}\quad s\in[0,1]\\}$

is bounded.

###### Proof.

Let $\tilde{\varrho}_{1},\tilde{\varrho}_{2}\in C([0,T];L^{2\gamma})$ and
$u_{1},u_{2}$ be the corresponding solutions to (3.2). As before, denote
$\Phi(\tilde{\varrho}_{i})=\varrho_{i}$, $i=1,2$.

Compactness. From the previous propositions we know that $\varrho\in
L^{\infty}(0,T;W^{1,2\gamma})$ and $\partial_{t}\varrho\in
L^{2}(0,T;W^{-1,2\gamma})$ and the bounds are uniform for bounded sets of
$\tilde{\varrho}$ in the given spaces. Therefore, the compactness of $\Phi$ in
$C([0,T];L^{2\gamma})$ follows straight from a variant of the Aubin–Lions
lemma from [21].

Continuity. We will estimate $u_{1}-u_{2}$ in terms of
$\tilde{\varrho}_{1}-\tilde{\varrho}_{2}$. We have

$-\mathrm{div}\,\Big{(}\mu_{0}(|\mathbb{D}u_{1}|)\mathbb{D}u_{1}-\mu_{0}(|\mathbb{D}u_{2}|)\mathbb{D}u_{2}\Big{)}-\nabla\Big{(}(1+\lambda(\mathrm{div}\,u_{1}))\mathrm{div}\,u_{1}-(1+\lambda(\mathrm{div}\,u_{2}))\mathrm{div}\,u_{2}\Big{)}\\\
-\Delta(u_{1}-u_{2})+\varepsilon\Delta^{2m}(u_{1}-u_{2})=-\nabla(\tilde{\varrho}_{1}^{\gamma}-\tilde{\varrho}_{2}^{\gamma}).$

Multiplying the above equality by $u_{1}-u_{2}$ and integrating over
$\mathbb{T}^{d}$, we get

$\mathcal{A}(u_{1},u_{2})+\int\varepsilon|\Delta^{m}(u_{1}-u_{2})|^{2}+|\nabla(u_{1}-u_{2})|^{2}\;\,\mathrm{d}x\\\
=\int(\tilde{\varrho}_{1}^{\gamma}-\tilde{\varrho}_{2}^{\gamma})(\mathrm{div}\,u_{1}-\mathrm{div}\,u_{2})\,\mathrm{d}x,$

where

$\mathcal{A}(u_{1},u_{2})=\int\left(\mu_{0}(|\mathbb{D}u_{1}|)\mathbb{D}u_{1}-\mu_{0}(|\mathbb{D}u_{2}|)\mathbb{D}u_{2}\right):(\mathbb{D}u_{1}-\mathbb{D}u_{2})\\\
+\Big{(}(1+\lambda(\mathrm{div}\,u_{1}))\mathrm{div}\,u_{1}-(1+\lambda(\mathrm{div}\,u_{2}))\mathrm{div}\,u_{2}\Big{)}(\mathrm{div}\,u_{1}-\mathrm{div}\,u_{2})\,\mathrm{d}x\geq
0$

from the monotonicity of the functions $B\mapsto\frac{B}{a+|B|}$ and
$s\mapsto\lambda(s)s$. In consequence, we have

$\displaystyle\|\nabla(u_{1}-u_{2})\|_{L^{2}(\mathbb{T}^{d})}^{2}+$
$\displaystyle\varepsilon\|\Delta^{m}(u_{1}-u_{2})\|_{L^{2}(\mathbb{T}^{d})}^{2}\leq\|\tilde{\varrho}_{1}^{\gamma}-\tilde{\varrho}_{2}^{\gamma}\|_{L^{2}(\mathbb{T}^{d})}\|\mathrm{div}\,(u_{1}-u_{2})\|_{L^{2}(\mathbb{T}^{d})}$
$\displaystyle\leq
C\left(\|\tilde{\varrho}_{1}\|_{L^{\infty}(\mathbb{T}^{d})}^{\gamma-1}+\|\tilde{\varrho}_{2}\|_{L^{\infty}(\mathbb{T}^{d})}^{\gamma-1}\right)\|\tilde{\varrho}_{1}-\tilde{\varrho}_{2}\|_{L^{2}(\mathbb{T}^{d})}\|\nabla(u_{1}-u_{2})\|_{L^{2}(\mathbb{T}^{d})}$
$\displaystyle\leq
C(\eta)\|\tilde{\varrho}_{1}-\tilde{\varrho}_{2}\|_{L^{2}(\mathbb{T}^{d})}^{2}+\eta\|\nabla(u_{1}-u_{2})\|_{L^{2}(\mathbb{T}^{d})}^{2}.$

Hence, choosing $\eta$ small enough,

$\displaystyle\|u_{1}-u_{2}\|_{L^{\infty}(0,T;W^{1,\infty})}\leq
C\|\Delta^{m}(u_{1}-u_{2})\|_{L^{\infty}(0,T;L^{2})}$ $\displaystyle\leq
C\|\tilde{\varrho}_{1}-\tilde{\varrho}_{2}\|_{L^{\infty}(0,T;L^{2})}$
$\displaystyle\leq
C\|\tilde{\varrho}_{1}-\tilde{\varrho}_{2}\|_{L^{\infty}(0,T;L^{2\gamma})}.$

Moreover, $\varrho_{1}-\varrho_{2}$ satisfy

$\partial_{t}(\varrho_{1}-\varrho_{2})+\delta(\varrho_{1}^{\beta}-\varrho_{2}^{\beta})-\delta\Delta(\varrho_{1}-\varrho_{2})=-\mathrm{div}\,(\varrho_{1}u_{1}-\varrho_{2}u_{2})$
(3.6)

with

$(\varrho_{1}-\varrho_{2})_{|_{t=0}}=0.$

Let us now estimate $\|\varrho_{1}-\varrho_{2}\|_{L^{\infty}(0,T;L^{p})}$.
First, we write $\mathrm{div}\,(\varrho_{1}u_{1}-\varrho_{2}u_{2})$ as

$\mathrm{div}\,(\varrho_{1}u_{1}-\varrho_{2}u_{2})=u_{1}\nabla(\varrho_{1}-\varrho_{2})+\nabla\varrho_{2}(u_{1}-u_{2})+\mathrm{div}\,u_{1}(\varrho_{1}-\varrho_{2})+\varrho_{2}(\mathrm{div}\,u_{1}-\mathrm{div}\,u_{2}).$

Then, multiplying (3.6) by
$p|\varrho_{1}-\varrho_{2}|^{p-2}(\varrho_{1}-\varrho_{2})$ and integrating,
we obtain

$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int|\varrho_{1}-\varrho_{2}|^{p}\,\mathrm{d}x+$
$\displaystyle\delta
p\int\left(|\varrho_{1}-\varrho_{2}|^{p-2}(\varrho_{1}^{\beta}-\varrho_{2}^{\beta})(\varrho_{1}-\varrho_{2})+(p-1)|\varrho_{1}-\varrho_{2}|^{p-2}|\nabla(\varrho_{1}-\varrho_{2})|^{2}\right)\,\mathrm{d}x$
$\displaystyle=$
$\displaystyle-(p-1)\int|\varrho_{1}-\varrho_{2}|^{p}\mathrm{div}\,u_{1}\,\mathrm{d}x$
$\displaystyle-\int(\nabla\varrho_{2}(u_{1}-u_{2})+\varrho_{2}(\mathrm{div}\,u_{1}-\mathrm{div}\,u_{2}))|\varrho_{1}-\varrho_{2}|^{p-2}(\varrho_{1}-\varrho_{2})\,\mathrm{d}x.$

In consequence, as
$(\varrho_{1}^{\beta}-\varrho_{2}^{\beta})(\varrho_{1}-\varrho_{2})\geq 0$,

$\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int|\varrho_{1}-\varrho_{2}|^{p}\,\mathrm{d}x\leq(p-1)\|u_{1}\|_{W^{1,\infty}}\int|\varrho_{1}-\varrho_{2}|^{p}\,\mathrm{d}x\\\
+\|\varrho_{2}\|_{W^{1,p}(\mathbb{T}^{d})}\|u_{1}-u_{2}\|_{W^{1,\infty}(\mathbb{T}^{d})}\|\varrho_{1}-\varrho_{2}\|_{L^{p}(\mathbb{T}^{d})}^{p-1}.$

Therefore from Gronwall’s lemma

$\|\varrho_{1}-\varrho_{2}\|_{L^{\infty}(0,T;L^{p})}\leq
C\|u_{1}-u_{2}\|_{L^{\infty}(0,T;W^{1,\infty})},$

where $C$ depends on $T$, $\|u_{1}\|_{L^{1}(0,T;W^{1,\infty})}$ and
$\|\varrho_{2}\|_{L^{2}(0,T;W^{1,p})}$. In particular,

$\|\varrho_{1}-\varrho_{2}\|_{L^{\infty}(0,T;L^{2\gamma})}\leq
C\|u_{1}-u_{2}\|_{L^{\infty}(0,T;W^{1,\infty})}\leq
C\|\tilde{\varrho}_{1}-\tilde{\varrho}_{2}\|_{L^{\infty}(0,T;L^{2\gamma})}.$

Estimates for the fixed points. To complete the proof of the Proposition, we
need to check if the points satisfying $\varrho=s\Phi(\varrho)$ are bounded in
$L^{\infty}(0,T;L^{2\gamma})$ for any $s\in[0,1]$. Throughout the proof we
will denote by $C$ various constants independent on $s$. We have for $s>0$ (if
$s=0$, the proof is trivial)

$\frac{1}{s}\partial_{t}\varrho+\frac{1}{s}\mathrm{div}\,(\varrho
u)+\frac{\delta}{s^{\beta}}\varrho^{\beta}=\frac{1}{s}\delta\Delta\varrho$

and

$-\mathrm{div}\,\left(\mu_{0}(|\mathbb{D}u|)\mathbb{D}u\right)-\Delta
u-\nabla(1+\lambda(\mathrm{div}\,u))\mathrm{div}\,u+\varepsilon\Delta^{2m}u+\nabla\varrho^{\gamma}=0.$

Multiplying the momentum equation by $u$ and integrating, we obtain
analogously as for the a priori estimates

$\|u\|_{L^{2}(0,T;W^{1,2})}^{2}+\varepsilon\|\Delta^{m}u\|_{L^{2}((0,T)\times\mathbb{T}^{d})}^{2}+\|\varrho\|_{L^{\infty}(0,T;L^{\gamma})}^{\gamma}\\\
+\frac{\delta}{s^{\beta-1}}\frac{\gamma}{\gamma-1}\int_{0}^{T}\int\varrho^{\beta+\gamma-1}\,\mathrm{d}x\,\mathrm{d}t+\delta\gamma\int_{0}^{T}\int|\nabla\varrho|^{2}\varrho^{\gamma-2}\,\mathrm{d}x\,\mathrm{d}t\leq\int\varrho_{0,\varepsilon}^{\gamma}\,\mathrm{d}x\leq
C.$

Repeating the estimate from Proposition 3.3, we also get

$\|\varrho\|_{L^{\infty}(0,T;L^{2\gamma})}\leq\|\varrho_{0,\varepsilon}\|_{L^{\infty}(\mathbb{T}^{d})}e^{\|u\|_{L^{1}(0,T;W^{1,\infty})}}\leq
C.$

∎

In consequence, the assumptions of Theorem 3.1 are satisfied and there exists
at least one solution to (3.1) on $[0,T]\times\mathbb{T}^{d}$ for arbitrary
$T>0$.

##  4 Compactness

We will now prove that we can pass to the limit with $\delta,\varepsilon\to 0$
to obtain the solutions to system (1.5). First, we will pass to the limit with
$\varepsilon\to 0$ and then we improve the estimates on $\varrho$ uniform in
$\delta$ and perform the second limit passage. Below, by $\overline{f}$ we
will denote the weak limit of $f$.

###  4.1 Limit passage with $\varepsilon\to 0$

Let $(\varrho_{\delta,\varepsilon},u_{\delta,\varepsilon})$ be a solution to
(3.1). We have the following estimates uniform in $\varepsilon$ (here we use
the lower bounds on $\beta)$:

$\displaystyle\|u_{\delta,\varepsilon}\|_{L^{2}(0,T;W^{1,2})}^{2}$
$\displaystyle+\|\varrho_{\delta,\varepsilon}\|_{L^{\infty}(0,T;L^{\gamma})}^{\gamma}+\delta\|\nabla\varrho_{\delta,\varepsilon}^{\gamma/2}\|_{L^{2}((0,T)\times\mathbb{T}^{d})}^{2}$
$\displaystyle+\delta\|\varrho_{\delta,\varrho}\|_{L^{\gamma+\beta-1}((0,T)\times\mathbb{T}^{d})}^{\gamma+\beta-1}+\delta\|\nabla\varrho_{\delta,\varepsilon}\|_{L^{2}((0,T)\times\mathbb{T}^{d})}^{2}\leq
C.$

In particular, at least up to a subsequence,

$u_{\delta,\varepsilon}\rightharpoonup u_{\delta}\quad\text{in}\quad
L^{2}(0,T;W^{1,2}).$

Moreover, as $\beta+\gamma-1\geq 2\gamma$, we know that

$\|\varrho_{\delta,\varepsilon}^{\gamma}\|_{L^{2}((0,T)\times\mathbb{T}^{d})}\leq
C(\delta)$

and

$\|\varrho_{\delta,\varepsilon}u_{\delta,\varepsilon}\|_{L^{p}((0,T)\times\mathbb{T}^{d})}\leq
C(\delta)$

for some suitable $p<2$. In consequence,

$\|\nabla\varrho_{\delta,\varepsilon}\|_{L^{2}((0,T)\times\mathbb{T}^{d})},\;\|\varrho_{\delta,\varepsilon}\|_{L^{p}((0,T)\times\mathbb{T}^{d})},\;\|\partial_{t}\varrho_{\delta,\varepsilon}\|_{L^{p}(0,T;W^{-1,p})}\leq
C(\delta).$

Therefore from the Aubin–Lions lemma
$\varrho_{\delta,\varepsilon}\to\varrho_{\delta}$ in
$L^{p}((0,T)\times\mathbb{T}^{d})$, (at least up to a subsequence). Then we
also have $\varrho_{\delta,\varepsilon}^{\gamma}\to\varrho_{\delta}^{\gamma}$
and $\varrho_{\delta,\varepsilon}^{\beta}\to\varrho_{\delta}^{\beta}$ in
suitable $L^{q}$ spaces. In consequence, we are able to pass to the limit in
the continuity equation. For the momentum equation, first note that the
regularizing term satisfies

$\varepsilon^{1/2}\|\Delta^{m}u_{\delta,\varepsilon}\|_{L^{2}((0,T)\times\mathbb{T}^{d})}\leq
C$

and thus in the weak formulation

$\varepsilon\int_{0}^{T}\int\Delta^{m}u_{\delta,\varepsilon}\cdot\Delta^{m}\phi\,\mathrm{d}x\,\mathrm{d}t\to
0$

for $\phi\in C_{0}^{\infty}((0,T)\times\mathbb{T}^{d})$. Therefore in the weak
formulation we obtain

$\int_{0}^{T}\int\overline{\mu_{0}(|\mathbb{D}u_{\delta}|)\mathbb{D}u_{\delta}}:\mathbb{D}\varphi+\nabla
u_{\delta}:\nabla\varphi+\mathrm{div}\,u_{\delta}\mathrm{div}\,\varphi+\overline{\lambda(\mathrm{div}\,u_{\delta})\mathrm{div}\,u_{\delta}}\mathrm{div}\,\varphi-\varrho_{\delta}^{\gamma}\mathrm{div}\,\varphi\,\mathrm{d}x\,\mathrm{d}t=0.$

Testing by $u_{\delta}$, we get

$\int_{0}^{T}\int\overline{\mu_{0}(|\mathbb{D}u_{\delta}|)\mathbb{D}u_{\delta}}:\mathbb{D}u_{\delta}+|\nabla
u_{\delta}|^{2}+(\mathrm{div}\,u_{\delta})^{2}+\overline{\lambda(\mathrm{div}\,u_{\delta})\mathrm{div}\,u_{\delta}}\mathrm{div}\,u_{\delta}-\varrho_{\delta}^{\gamma}\mathrm{div}\,u_{\delta}\,\mathrm{d}x\,\mathrm{d}t=0.$

On the other hand,

$\int_{0}^{T}\int\overline{\mu_{0}(|\mathbb{D}u_{\delta}|)|\mathbb{D}u_{\delta}|^{2}}+\overline{\lambda(\mathrm{div}\,u_{\delta})(\mathrm{div}\,u_{\delta})^{2}}\,\mathrm{d}x\,\mathrm{d}t+\limsup_{\varepsilon\to
0}\int_{0}^{T}\int|\nabla
u_{\delta,\varepsilon}|^{2}+(\mathrm{div}\,u_{\delta,\varepsilon})^{2}\,\mathrm{d}x\,\mathrm{d}t\\\
-\int_{0}^{T}\int\varrho_{\delta}^{\gamma}\mathrm{div}\,u\,\mathrm{d}x\,\mathrm{d}t\leq
0.$

Therefore using the monotonicity of $\mu_{0}(|\cdot|)\cdot$ and
$\lambda(|\cdot|)\cdot$ and weak lower semicontinuity of the norm, we obtain
the convergence $\nabla u_{\delta,\varepsilon}\to\nabla u_{\delta}$ in
$L^{2}((0,T)\times\mathbb{T}^{d})$, which allows us to pass to the limit in
the remaining nonlinear terms.

###  4.2 Limit passage with $\delta\to 0$

Let $(\varrho_{\delta},u_{\delta})$ be the function obtained in the previous
section, solving

$\displaystyle\varrho_{t}+\mathrm{div}\,(\varrho u)+\delta\varrho^{\beta}$
$\displaystyle=\delta\Delta\varrho,$ (4.1)
$\displaystyle-\mathrm{div}\,((\mu_{0}(|\mathbb{D}u|)+1)\mathbb{D}u)-\nabla(\lambda(\mathrm{div}\,u)\mathrm{div}\,u)+\nabla\varrho^{\gamma}$
$\displaystyle=0.$

Note that repeating the calculations from Section 2, we get the estimate

$\|\mathrm{div}\,u_{\delta}+\lambda(\mathrm{div}\,u_{\delta})\mathrm{div}\,u_{\delta}-\varrho_{\delta}^{\gamma}\|_{L^{\infty}(0,T;BMO)}\leq
C.$

Moreover, using the uniform estimates on $\|u_{\delta}\|_{L^{2}(0,T;W^{1,2})}$
and $\|\varrho_{\delta}\|_{L^{\infty}(0,T;L^{\gamma})}$, we will improve the
integrability of $\varrho_{\delta}$ uniformly in $\delta$.

Let $p>1$. Define the function $P_{k}(\varrho)$ as

$P_{k}(\varrho)=\varrho\int_{0}^{\varrho}\frac{T_{k}(z)^{p}}{z^{2}}\,\mathrm{d}z,$

where $T_{k}\in C^{\infty}([0,\infty))$ is the truncation operator, namely
$T_{k}(z)=z$ for $z<k$, $T_{k}(z)=k+1$ for $z>2k$, $T_{k}^{\prime}(z)\geq 0$
as well as $T_{k}(z)\nearrow z$ as $k\to\infty$. It is easy to see that
$P_{k}(\varrho)\nearrow\varrho^{p}$ as well. Using the renormalized continuity
equation, we get

$\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int
P_{k}(\varrho_{\delta})\,\mathrm{d}x+\delta\int\left(\varrho_{\delta}^{\beta}P_{k}^{\prime}(\varrho_{\delta})+\frac{pT_{k}(\varrho_{\delta})^{p-1}T_{k}^{\prime}(\varrho_{\delta})}{\varrho_{\delta}}|\nabla\varrho_{\delta}|^{2}\right)\,\mathrm{d}x=-\int
T_{k}(\varrho_{\delta})^{p}\mathrm{div}\,u_{\delta}\,\mathrm{d}x.$ (4.2)

Now, let us test the momentum equation by the function

$\psi=\Delta^{-1}\nabla\left(T_{k}(\varrho_{\delta})^{p}-\\{T_{k}(\varrho_{\delta})^{p}\\}\right).$

We have

$\displaystyle\int_{0}^{T}\int$
$\displaystyle\varrho_{\delta}^{\gamma}T_{k}(\varrho_{\delta})^{p}\,\mathrm{d}x$
$\displaystyle\leq
C\big{(}\|\mu_{0}(|\mathbb{D}u_{\delta}|)\mathbb{D}u_{\delta}\|_{L^{\infty}((0,T)\times\mathbb{T}^{d})}+\|\lambda(\mathrm{div}\,u_{\delta})\mathrm{div}\,u_{\delta}\|_{L^{\infty}((0,T)\times\mathbb{T}^{d})}\big{)}\|\nabla\psi\|_{L^{1+\frac{\gamma}{p}}((0,T)\times\mathbb{T}^{d})}$
$\displaystyle+2\int_{0}^{T}\int
T_{k}(\varrho_{\delta})^{p}\mathrm{div}\,u_{\delta}\,\mathrm{d}x+\|\varrho_{\delta}\|_{L^{\infty}(0,T;L^{\gamma})}\int_{0}^{T}\int
T_{k}(\varrho_{\delta})^{p}\,\mathrm{d}x.$

By Cauchy inequality,

$\displaystyle\int_{0}^{T}\int\varrho_{\delta}^{\gamma}T_{k}(\varrho_{\delta})^{p}\,\mathrm{d}x\leq$
$\displaystyle
C\|T_{k}(\varrho_{\delta})\|_{L^{p+\gamma}((0,T)\times\mathbb{T}^{d})}^{p}+2\int_{0}^{T}\int
T_{k}(\varrho_{\delta})^{p}\mathrm{div}\,u_{\delta}\,\mathrm{d}x$
$\displaystyle\leq\eta\|T_{k}(\varrho_{\delta})\|_{L^{p+\gamma}((0,T)\times\mathbb{T}^{d})}^{p+\gamma}+C(\eta)+2\int_{0}^{T}\int
T_{k}(\varrho_{\delta})^{p}\mathrm{div}\,u_{\delta}\,\mathrm{d}x.$

As $T_{k}(\varrho)^{p+\gamma}\leq\varrho^{\gamma}T_{k}(\varrho)^{p}$, for
sufficiently small $\eta$ we get

$\int_{0}^{T}\int\varrho_{\delta}^{\gamma}T_{k}(\varrho_{\delta})^{p}\,\mathrm{d}x\,\mathrm{d}t-2\int_{0}^{T}\int
T_{k}(\varrho_{\delta})^{p}\mathrm{div}\,u_{\delta}\,\mathrm{d}x\leq C.$

Therefore using (4.2), we get

$\displaystyle\int_{0}^{T}\int\varrho_{\delta}^{\gamma}T_{k}(\varrho_{\delta})^{p}\,\mathrm{d}x\,\mathrm{d}t+\sup_{t\in(0,T]}\int
P_{k}(\varrho_{\delta}(t,\cdot))\,\mathrm{d}x$
$\displaystyle+\delta\int\left(\varrho_{\delta}^{\beta}P_{k}^{\prime}(\varrho_{\delta})+\frac{pT_{k}(\varrho_{\delta})^{p-1}T_{k}^{\prime}(\varrho_{\delta})}{\varrho_{\delta}}|\nabla\varrho_{\delta}|^{2}\right)\,\mathrm{d}x\leq
C(T,p).$

We pass to the limit with $k\to\infty$ using monotone convergence theorem and
in consequence

$\|\varrho_{\delta}\|_{L^{\infty}(0,T;L^{p})}\leq C(T,p)\quad\text{for
any}\quad p<\infty.$

Having that estimate, we are ready to pass to the limit with $\delta\to 0$.

From the estimates uniform in $\delta$, we know that in particular

$\displaystyle u_{\delta}\rightharpoonup u$ $\displaystyle\text{in}\quad
L^{2}(0,T;W^{1,2}),$ $\displaystyle\varrho_{\delta}\rightharpoonup^{*}\varrho$
$\displaystyle\text{in}\quad L^{\infty}(0,T;L^{\gamma}),$
$\displaystyle\varrho_{\delta}^{\gamma}\rightharpoonup^{*}\overline{\varrho^{\gamma}}$
$\displaystyle\text{in}\quad L^{\infty}(0,T;L^{p})$

and

$\mu_{0}(\mathbb{D}u_{\delta})\mathbb{D}u_{\delta},\;\lambda(\mathrm{div}\,u_{\delta})\mathrm{div}\,u_{\delta}\rightharpoonup^{*}\overline{\mu_{0}(|\mathbb{D}u|)\mathbb{D}u},\;\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}\quad\text{in}\quad
L^{\infty}((0,T)\times\mathbb{T}^{d}).$

Moreover,
$\|\mathrm{div}\,u+\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}-\overline{\varrho^{\gamma}}\|_{L^{\infty}(0,T;BMO)}\leq
C$. Note that from the continuity equation it also follows that
$\varrho_{\delta}\to\varrho$ in $C([0,T];W^{-1,r})$ for a suitable $r$, and in
consequence $\overline{\varrho u}=\varrho u$. Having the above estimates and
testing the continuity equation by $\varrho_{\delta}$, we also obtain

$\delta^{1/2}\|\nabla\varrho_{\delta}\|_{L^{2}((0,T)\times\mathbb{T}^{d})}\leq
C.$

Then for $\phi\in C_{0}^{\infty}((0,T)\times\Omega)$, together with the
estimate on $\|\varrho_{\delta}\|_{L^{\infty}(0,T;L^{p})}$,

$\delta\int_{0}^{T}\int\left(\varrho_{\delta}^{\beta}\varphi+\nabla\varrho_{\delta}\cdot\nabla\phi\right)\,\mathrm{d}x\,\mathrm{d}t\to
0\quad\text{with}\quad\delta\to 0.$

In consequence, $(\varrho,u)$ satisfies the continuity equation (1.5) in the
renormalized sense.

Next, we will pass to the limit in the momentum equation and apply an argument
from [13]. Passing to the limit in the weak formulation, we get for any
$\phi\in C_{0}^{\infty}((0,T)\times\mathbb{T}^{d})$ and $t\leq T$

$\int_{0}^{t}\int\overline{\mu_{0}(|\mathbb{D}u|)\mathbb{D}u}:\mathbb{D}\phi+\nabla
u\cdot\nabla\phi+\mathrm{div}\,u\
\mathrm{div}\,\phi+\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}\
\mathrm{div}\,\phi\;\,\mathrm{d}x\,\mathrm{d}s=\int_{0}^{t}\int\overline{\varrho^{\gamma}}\mathrm{div}\,\phi\,\mathrm{d}x\,\mathrm{d}s.$
(4.3)

The regularity of $u$ allows us to put $\phi=u$ in (4.3) and then

$\int_{0}^{t}\int\overline{\mu_{0}(|\mathbb{D}u|)\mathbb{D}u}:\mathbb{D}u+|\nabla
u|^{2}+(\mathrm{div}\,u)^{2}+\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}\
\mathrm{div}\,u\;\,\mathrm{d}x\,\mathrm{d}s=\int_{0}^{t}\int\overline{\varrho^{\gamma}}\mathrm{div}\,u\,\mathrm{d}x\,\mathrm{d}s.$
(4.4)

On the other hand, the solutions to approximate equation (4.1) satisfy

$\int_{0}^{t}\int\mu_{0}(|\mathbb{D}u_{\delta}|)|\mathbb{D}u_{\delta}|^{2}\,\mathrm{d}x\,\mathrm{d}s+|\nabla
u_{\delta}|^{2}+(\mathrm{div}\,u_{\delta})^{2}+\lambda(\mathrm{div}\,u_{\delta})(\mathrm{div}\,u_{\delta})^{2}\;\,\mathrm{d}x\,\mathrm{d}s\\\
+\frac{1}{\gamma-1}\int\varrho_{\delta}^{\gamma}(t,\cdot)\,\mathrm{d}x\leq\frac{1}{\gamma-1}\int\varrho_{0,\delta}^{\gamma}\,\mathrm{d}x.$
(4.5)

Using the monotonicity of $\mu_{0}(|\mathbb{D}u|)\mathbb{D}u$ and
$\lambda(\mathrm{div}\,u)\mathrm{div}\,u$, we know that

$\overline{\mu_{0}(|\mathbb{D}u|)|\mathbb{D}u|^{2}}\geq\overline{\mu_{0}(|\mathbb{D}u|)\mathbb{D}u}:\mathbb{D}u$
(4.6)

and

$\overline{\lambda(\mathrm{div}\,u)(\mathrm{div}\,u)^{2}}\geq\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}\
\mathrm{div}\,u.$ (4.7)

Therefore, taking $\liminf_{\delta\to 0}$ in the energy inequality (4.5), we
obtain

$\int_{0}^{t}\int\overline{\mu_{0}(|\mathbb{D}u|)\mathbb{D}u}:\mathbb{D}u+|\nabla
u|^{2}+(\mathrm{div}\,u)^{2}+\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}\
\mathrm{div}\,u\;\,\mathrm{d}x\,\mathrm{d}s\\\
+\frac{1}{\gamma-1}\int\overline{\varrho^{\gamma}}(t,\cdot)\,\mathrm{d}x\leq\frac{1}{\gamma-1}\int\varrho_{0}^{\gamma}\,\mathrm{d}x.$

Comparing the last equation with (4.4), we get

$\frac{1}{\gamma-1}\int\overline{\varrho^{\gamma}}(t,\cdot)\,\mathrm{d}x-\frac{1}{\gamma-1}\int\varrho_{0}^{\gamma}\,\mathrm{d}x\leq-\int_{0}^{t}\int\overline{\varrho^{\gamma}}\mathrm{div}\,u\,\mathrm{d}x\,\mathrm{d}s.$

We would like to estimate
$\displaystyle\int\overline{\varrho^{\gamma}}(t,\cdot)-\varrho^{\gamma}(t,\cdot)\,\mathrm{d}x$.
As we already know that $\varrho$ satisfies the continuity equation in the
renormalized sense, we have

$\frac{1}{\gamma-1}\int\varrho^{\gamma}(t,\cdot)\,\mathrm{d}x-\frac{1}{\gamma-1}\int\varrho_{0}^{\gamma}\,\mathrm{d}x=-\int_{0}^{t}\int\varrho^{\gamma}\mathrm{div}\,u\;\,\mathrm{d}x\,\mathrm{d}s.$
(4.8)

Therefore

$\frac{1}{\gamma-1}\int\big{(}\overline{\varrho^{\gamma}}(t,\cdot)-\varrho^{\gamma}(t,\cdot)\big{)}\,\mathrm{d}x\leq-\int_{0}^{t}\int\big{(}\overline{\varrho^{\gamma}}-\varrho^{\gamma}\big{)}\mathrm{div}\,u\,\mathrm{d}x\,\mathrm{d}s.$

We now use the fact that
$\mathrm{div}\,u+\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}-\overline{\varrho^{\gamma}}\in
L^{\infty}(0,T;BMO)$ and the logarithmic inequality (1.9). As
$\overline{\varrho^{\gamma}}\geq\varrho^{\gamma}$ and
$\varrho^{\gamma},\overline{\varrho^{\gamma}}\in L^{\infty}(0,T;L^{p})$ for
any $p<\infty$, we have

$\displaystyle-\int_{0}^{t}\int\big{(}\overline{\varrho^{\gamma}}-\varrho^{\gamma}\big{)}\mathrm{div}\,u\;\,\mathrm{d}x\,\mathrm{d}s=$
$\displaystyle-\int_{0}^{t}\int\big{(}\overline{\varrho^{\gamma}}-\varrho^{\gamma}\big{)}(\mathrm{div}\,u+\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}-\overline{\varrho^{\gamma}})\,\mathrm{d}x\,\mathrm{d}s$
$\displaystyle-\int_{0}^{t}\int\big{(}\overline{\varrho^{\gamma}}-\varrho^{\gamma}\big{)}\overline{\varrho^{\gamma}}\,\mathrm{d}x\,\mathrm{d}s+\int_{0}^{t}\int\big{(}\overline{\varrho^{\gamma}}-\varrho^{\gamma}\big{)}\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}\,\mathrm{d}x\,\mathrm{d}t$
$\displaystyle\leq$
$\displaystyle-\int_{0}^{t}\int\big{(}\overline{\varrho^{\gamma}}-\varrho^{\gamma}\big{)}(\mathrm{div}\,u+\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}-\overline{\varrho^{\gamma}})\,\mathrm{d}x\,\mathrm{d}s$
$\displaystyle+\int_{0}^{t}\int\big{(}\overline{\varrho^{\gamma}}-\varrho^{\gamma}\big{)}\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}\,\mathrm{d}x\,\mathrm{d}s$
$\displaystyle\leq$ $\displaystyle
C\int_{0}^{t}\int\big{(}\overline{\varrho^{\gamma}}-\varrho^{\gamma}\big{)}\,\mathrm{d}x\left(\left|\ln\left(\int\big{(}\overline{\varrho^{\gamma}}-\varrho^{\gamma}\big{)}\,\mathrm{d}x\right)\right|+1\right)\,\mathrm{d}s$
$\displaystyle+\|\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}\|_{L^{\infty}((0,T)\times\mathbb{T}^{d})}\int_{0}^{t}\int\Big{(}\overline{\varrho^{\gamma}}-\varrho^{\gamma}\Big{)}\,\mathrm{d}x\,\mathrm{d}s,$

where $C$ depends on
$\|\mathrm{div}\,u+\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}-\overline{\varrho^{\gamma}}\|_{L^{\infty}(0,T;BMO)}$
and $\|\overline{\varrho^{\gamma}}-\varrho^{\gamma}\|_{L^{\infty}(0,T;L^{q})}$
for some $q>2$. Thus denoting
$y(t)=\int\big{(}\overline{\varrho^{\gamma}}-\varrho^{\gamma}\big{)}\,\mathrm{d}x$,
we have the inequality

$y(t)\leq C\int_{0}^{t}y(s)(|\ln y(s)|+1)\,\mathrm{d}s\quad\text{with}\quad
y(0)=0.$

As the ordinary differential equation $z^{\prime}=Cz(|\ln z|+1)$ has a unique
solution, we can apply the comparison criterion and in consequence $y\equiv 0$
on the interval $[0,t_{*}]$ for some $t_{*}>0$. Then we can apply the same
analysis on the consecutive intervals of length $t_{*}$ and in the end
$y\equiv 0$ on the whole interval $[0,T]$. From this, as
$\overline{\varrho^{\gamma}}\geq\varrho^{\gamma}$, it follows that in fact
$\overline{\varrho^{\gamma}}=\varrho^{\gamma}$. Now taking again the limit in
(4.5) and subtracting (4.4), we get due to fact that
$\varrho^{\gamma}=\overline{\varrho^{\gamma}}$ and (4.8)

$\displaystyle 0$ $\displaystyle\geq\int_{0}^{t}\int\Big{(}\overline{|\nabla
u|^{2}}-|\nabla
u|^{2}\Big{)}\,\mathrm{d}x\,\mathrm{d}s\geq\lim_{\varepsilon\to
0}\int_{0}^{t}\int|\nabla u_{\varepsilon}-\nabla
u|^{2}\,\mathrm{d}x\,\mathrm{d}t.$

Therefore $\nabla u_{\varepsilon}\to\nabla u$ in $L^{2}((0,T)\times\Omega)$
and (for possibly another subsequence) the sequence converges also a.e. In
consequence, by virtue of the Lebesgue dominated convergence theorem,

$\overline{\mu_{0}(|\mathbb{D}u|)\mathbb{D}u}=\mu_{0}(|\mathbb{D}u|)\mathbb{D}u$

and

$\overline{\lambda(\mathrm{div}\,u)\mathrm{div}\,u}=\lambda(\mathrm{div}\,u)\mathrm{div}\,u,$

and thus $(\varrho,u)$ satisfies the weak formulation of the system (1.5).

######  Remark 1.

The assumptions on $\mu_{0}$ and $\lambda$ and the used method allows the
situation when the viscosities are singular at $0$, e.g.
$\mu_{0}=\frac{1}{|\mathbb{D}u|}$. Note, however, that in this case, while
passing to the limit in the weak formulation, the term

$\int\frac{\mathbb{D}u}{|\mathbb{D}u|}:\nabla\varphi\,\mathrm{d}x$

is well defined by the values of $\mathbb{D}u$ provided $|\mathbb{D}u|>0$. For
$|\mathbb{D}u|=0$ it is just defined as the corresponding limit, which is not
necessarily equal to zero if $|\mathbb{D}u|$ is so, cf. e.g. [15] in a similar
context.

Acknowledgement: The work of M. P. was partially supported by the Czech
Science Foundation, project No.: 22-01591S, whereas the work of M. S. was
supported by National Science Centre grant 2018/29/B/ST1/00339. The results
were obtained during the internship of M. S. in Prague.

## References

  * [1] A. Abbatiello and E. Feireisl.  On a class of generalized solutions to equations describing incompressible viscous fluids.  Ann. Mat. Pura Appl. (4), 199(3):1183–1195, 2020.
  * [2] R. Andrášik and R. Vodák.  Compressible nonlinearly viscous fluids: asymptotic analysis in a 3D curved domain.  J. Math. Fluid Mech., 21(1):Paper No. 13, 27, 2019.
  * [3] F. Belblidia, T. Haroon, and M. F. Webster.  The Dynamics of Compressible Herschel—Bulkley Fluids in Die-Swell Flows, pages 425–434.  Springer Netherlands, Dordrecht, 2009.
  * [4] J. Blechta, J. Málek, and K. R. Rajagopal.  On the classification of incompressible fluids and a mathematical analysis of the equations that govern their motion.  SIAM J. Math. Anal., 52(2):1232–1289, 2020.
  * [5] M. Bulíček, P. Kaplický, and D. Pražák.  Uniqueness and regularity of flows of non-Newtonian fluids with critical power-law growth.  Math. Models Methods Appl. Sci., 29(6):1207–1225, 2019.
  * [6] J. Burczak, S. Modena, and L. Székelyhidi.  Non uniqueness of power-law flows.  Comm. Math. Phys., 388(1):199–243, 2021.
  * [7] Y. Cho and K. Kensey.  Effects of the non-Newtonian viscosity of blood on hemodynamics of diseased arterial flows.  Advances in Bioengineering, 15:147–148, Jan. 1989.
  * [8] Y. Damianou, G. C. Georgiou, and I. Moulitsas.  Combined effects of compressibility and slip in flows of a herschel–bulkley fluid.  Journal of Non-Newtonian Fluid Mechanics, 193:89–102, 2013.  Viscoplastic Fluids: From Theory to Application.
  * [9] L. De Rosa, M. Inversi, and G. Stefani.  Weak-strong uniqueness and vanishing viscosity for incompressible Euler equations in exponential spaces.  arXiv:2204.12779, 2022.
  * [10] L. Diening, M. Růžička, and J. Wolf.  Existence of weak solutions for unsteady motions of generalized Newtonian fluids.  Ann. Sc. Norm. Super. Pisa Cl. Sci. (5), 9(1):1–46, 2010.
  * [11] H. Eberlein and M. Ruzicka.  Existence of weak solutions for unsteady motions of Herschel–Bulkley fluids.  Journal of Mathematical Fluid Mechanics - J MATH FLUID MECH, 14:1–16, 09 2011.
  * [12] E. Feireisl.  On compactness of solutions to the compressible isentropic Navier-Stokes equations when the density is not square integrable.  Comment. Math. Univ. Carolin., 42(1):83–98, 2001.
  * [13] E. Feireisl, X. Liao, and J. Málek.  Global weak solutions to a class of non-Newtonian compressible fluids.  Mathematical Methods in the Applied Sciences, 38(16), 2015.
  * [14] B. H. Kjartanson, D. H. Shields, L. Domaschuk, and C.-S. Man.  The creep of ice measured with the pressuremeter.  Canadian Geotechnical Journal, 25(2):250–261, 1988.
  * [15] M. Łasica, S. Moll, and P. B. Mucha.  Total variation denoising in $l^{1}$ anisotropy.  SIAM J. Imaging Sci., 10(4):1691–1723, 2017.
  * [16] P.-L. Lions.  Mathematical Topics in Fluid Mechanics, Vol. 2: Compressible Models.  Oxford University Press, 1998.
  * [17] A. E. Mamontov.  Existence of global solutions to multidimensional equations for Bingham fluids.  Mat. Zametki, 82:560–577, 2007.
  * [18] A. E. Mamontov.  Global solvability of the multidimensional equations of compressible non-Newtonian fluids, transport equation and the Orlicz spaces.  Sib. Elektron. Mat. Izv., 6:120–165, 2009.
  * [19] C.-S. Man and Q.-X. Sun.  On the Significance of Normal Stress Effects in the Flow of Glaciers.  Journal of Glaciology, 33(115):268–273, 1987.  Publisher: Cambridge University Press.
  * [20] P. B. Mucha and W. Rusin.  Zygmund spaces, inviscid limit and uniqueness of Euler flows.  Communications in Mathematical Physics, 280:831–841, 06 2008.
  * [21] J. Simon.  Compact sets in the space $L^{p}(0,T;B)$.  Ann. Mat. Pura Appl. (4), 146:65–96, 1987.
  * [22] H. Steffan, W. Brandstätter, G. Bachler, and R. Pucher.  Comparison of Newtonian and Non-Newtonian Blood Flow in Stenotic Vessels using Numerical Simulation.  In D. W. Liepsch, editor, Biofluid Mechanics, pages 479–485, Berlin, Heidelberg, 1990. Springer.
  * [23] M. Szlenk.  Weak solutions for the Stokes system for compressible fluids with general pressure.  Journal of Differential Equations, 312:317–346, 2022.
  * [24] D. A. Yuen and A. V. Malevsky.  Strongly Chaotic Newtonian and Non-Newtonian Mantle Convection.  In D. A. Yuen, editor, Chaotic Processes in the Geological Sciences, The IMA Volumes in Mathematics and its Applications, pages 71–88, New York, NY, 1992. Springer US.