\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2014 (2014), No. 70, pp. 1--20.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2014 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2014/70\hfil Regularity of planar flows]
{Regularity of planar flows for shear-thickening fluids under
perfect slip boundary conditions}

\author[J. Tich\'y\hfil EJDE-2014/70\hfilneg]
{Jakub Tich\'y}  % in alphabetical order

\address{Jakub Tich\'y \newline
Charles University in Prague, Department of Mathematical Analysis,
Sokolovsk\'a 83, 186 75 Praha 8, Czech Republic. \newline
Czech Technical University in Prague, Faculty of Information
Technology, Th\'akurova 9, 160 00, Praha 6, Czech Republic}
\email{jakub.tichy@fit.cvut.cz}


\thanks{Submitted April 19, 2013. Published March 15, 2014.}
\subjclass[2000]{35B65, 35K51, 35Q35, 76D03}
\keywords{Generalized Newtonian fluid; regularity; 
\hfill\break\indent perfect slip boundary conditions}

\begin{abstract}
 For evolutionary planar flows of shear-thickening fluids in bounded
 domains we prove the existence of a solution with the H\"older continuous
 velocity gradients and pressure. The problem is equipped with perfect
 slip boundary conditions. We also show $L^q$ theory result for Stokes
 system under perfect slip boundary conditions.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}

We study flows of incompressible shear-thickening fluids, 
which in evolutionary case are governed by the following initial value problem
\begin{equation}\label{ns}
\begin{gathered}
\partial_t u -\operatorname{div} \mathcal{S}(Du) + (u\cdot\nabla)u + \nabla \pi = f,\quad 
\operatorname{div} u=0\quad \text{in }Q,\\
u(0,\cdot)=u_0\quad \text{in }\Omega,
\end{gathered}
\end{equation}
where $u$ is the velocity, $\pi$ represents the pressure, $f$ stands for the 
density of volume forces and $\mathcal{S}(Du)$ denotes the extra stress tensor. 
$Du$ is the symmetric part of the velocity gradient; i.e.,
 $Du=\frac{1}{2}[\nabla u + (\nabla u)^\top]$, $\Omega \subset \mathbb{R}^2$ 
is a bounded domain, $I=(0,T)$ denotes a finite time interval and $Q=I\times\Omega$.
We are interested in the case, when \eqref{ns} is equipped with the perfect 
slip boundary conditions
\begin{equation}\label{bcps}
u\cdot\nu=0,\quad [\mathcal{S}(Du)\nu]\cdot\tau=0\quad \text{on }I\times\partial\Omega,
\end{equation}
where $\tau$ is the tangent vector and $\nu$ is the outward normal to 
$\partial\Omega$.
The constitutive relation for $\mathcal{S}$ is given via the generalized viscosity
 $\mu$ and is of the form
$$
\mathcal{S}(Du):=\mu(|Du|)Du.
$$
The extra stress tensor $\mathcal{S}$ is assumed to possess $p$-potential
structure with $p\ge 2$. More precisely,  we can construct scalar 
potential $\Phi:[0,\infty) \mapsto [0,\infty)$ to the stress tensor $\mathcal S$;
 i.e.,
\[
\mathcal S(A)=\partial_{A}\Phi(|A|)= \Phi '(|A|)\frac{A}{|A|} \quad \forall
A\in\mathbb{R}^{2\times 2}_{\rm sym},
\]
such that $\Phi \in \mathcal{C}^{1,1}((0,\infty))\cap\mathcal{C}^{1}([0,\infty))$, 
$\Phi(0)=0$ and there exist $p\in[2,\infty)$ and $0<C_1\le C_2$ such that for 
all $A,B\in\mathbb{R}^{2\times 2}_{\rm sym}$
\begin{equation}\label{as}
C_1(1+|A|^2)^{\frac{p-2}{2}}|B|^2 \le \partial^2_A\Phi(|A|):B\otimes B 
\le C_2 (1+|A|^2)^{\frac{p-2}{2}}|B|^2.
\end{equation}

In the analysis of equations of fluid motions the question of H\"older continuity 
of velocity gradients is an important issue. For instance, in optimal control 
problems, global regularity results that guarantee boundedness of velocity 
gradients are needed in order to establish the existence of the weak solution 
for adjoint equation to the original problem and for linearised models. 
These results are closely related to the regularity of the coefficients in 
the main part of the associated differential operators and enable to derive 
corresponding optimality conditions, as is done for example in \cite{wr}. 
For optimal control of flows with shear dependent viscosities in the stationary 
case where the author is dealing with the lack of the regularity result we 
refer to \cite{arada1} and \cite{arada2}.

H\"older continuity of velocity gradients is also important when one studies 
exponential attractors. With such a regularity it is possible to show the 
differentiability of the solution operator with respect to the initial condition, 
which is the key technical step in the method of Lyapunov exponents. 
Differentiability of the solution is equivalent to the linearisation  
of the equation around particular solution which is used to study infinitesimal 
volume elements and leads to sharp dimension estimates of the global attractor. 
This is done for example in \cite{kp}.

This article closely follows \cite{kap}, where P. Kaplick\'y shows H\"older 
continuity of velocity gradients and pressure for \eqref{ns} with $p\in[2,4)$ 
under no slip boundary conditions. Based on the same structure of the proof 
and using the results from \cite{kt} we extend the result to perfect slip 
boundary conditions and $p\in[2,\infty)$. Although some steps of the proof 
in \cite{kap} can be easily modified, we have to overcome a new difficulties 
connected to the another type of boundary conditions. First of all, 
the $L^p$ theory result for the Stokes problem equipped with perfect slip 
boundary conditions has to be established. Keep at our disposal the paper 
\cite{kt}, we are able to cover the case $p\ge 4$. From the point of 
application it would be very interesting to obtain also the result for 
the case $p\in(1,2)$ for perfect slip or homogeneous Dirichlet boundary condition.

The idea of the proof goes back to \cite{ns}, where the authors show that 
every weak solution $u$ of $\partial_tu-\operatorname{div} (\mathcal S(\nabla u))=0$ 
in $Q$ has locally H\"older continuous gradient in case that 
$\Omega \subset \mathbb{R}^2$ and $p=2$. This result was extended in \cite{fs} 
to the case $p\in(1,2)$. Regularity of $\partial_t u$ is shown first and after 
moving $\partial_t u$ to the right hand side the stationary $L^q$ theory is applied.

In the case of generalized Newtonian fluids this method was modified in \cite{kms3}, 
where the authors consider the shear-thinning fluid model with periodic boundary 
conditions. In contrary to \cite{ns} the regularity of $\partial_t u$ and $\nabla u$ 
had to be obtained at once. The authors showed that velocity gradients are H\"older 
continuous for $p\in(4/3,2]$. These results were extended to electro-rheological 
fluids and non-zero initial condition in \cite{die}.

Among many works concerning regularity theory for generalized Newtonian fluids 
we would like to mention two papers dealing with the stationary case. 
In \cite{kms2} the stationary version of \eqref{ns} under homogeneous Dirichlet 
boundary conditions is considered. The same authors later in \cite{kms1} studied 
the problem equipped with non-homogeneous Dirichlet boundary conditions with two 
types of restriction on boundary data and perfect slip boundary conditions.

Let $E$ be a Banach space and $\alpha \in (0,1)$, $p,q\in[1,\infty)$, 
$s \in\mathbb{R}$. In this paper we use standard notation
for Lebesgue spaces $L^q(\Omega)$, Sobolev-Slobodecki\u\i spaces
 $W^{s,q}(\Omega)$, Bochner spaces $L^q(I,E)$ and $W^{\alpha,q}(I,E)$. 
(We do not use different notation for scalar, vector-valued or tensor-valued 
functions).

By $H^s_q(\Omega)$ we mean Bessel potential spaces and $B^s_{p,q}(\Omega)$ 
are Besov spaces.
$BUC$ stands for bounded and uniformly continuous functions.

Since the domain $\Omega$ is in our case at least $\mathcal{C}^{2,1}$, 
we can define $L^q_\sigma(\Omega)$ and $W^{1,q}_{\sigma}(\Omega)$ as follows:
\begin{gather*}
L^q_\sigma(\Omega)=\{\varphi\in L^q(\Omega),\operatorname{div} \varphi =0
\text{ in }\Omega,\,\varphi\cdot\nu=0\text{ on }\partial\Omega\},\\
W^{1,q}_{\sigma}(\Omega)=\{\varphi \in W^{1,q}(\Omega),
\operatorname{div}\varphi=0,\text{ in }\Omega,\,\varphi\cdot \nu=0 
\text{ on }\partial\Omega\}.
\end{gather*}
The duality between Banach space $E$ and its dual $E'$ is denoted 
by $\langle\cdot,\cdot\rangle$.
Set $W^{-1,p'}_\sigma(\Omega):=(W^{1,p}_\sigma(\Omega))'$.

We begin with the definition of the weak solution to the problem
 \eqref{ns} with \eqref{bcps}.

\begin{definition}\label{dws} \rm
Let $f \in L^{p'}(I,W^{-1,p'}_\sigma(\Omega))$, $p\in[2,\infty)$ and 
$u_0 \in L^2(\Omega)$. We say that the function $u: Q \mapsto \mathbb{R}^2$ 
is a weak solution to the problem \eqref{ns} with \eqref{bcps}, if 
$u \in L^\infty(I,L^2(\Omega))\cap L^p(I,W^{1,p}_\sigma(\Omega))$, 
$\partial_t u \in L^{p'}(I,W^{-1,p'}_\sigma(\Omega))$, $u(0,\cdot)=u_0$ 
in $L^2(\Omega)$ and weak formulation
\[ %\label{weak.form}
\int_I\langle\partial_t u,\varphi\rangle\,\mathrm{d}t + \int_Q\mathcal{S}(Du)\!:\!D\varphi\,\mathrm{d}x\,\mathrm{d}t
 + \int_Q(u\cdot\nabla)u\varphi\,\mathrm{d}x\,\mathrm{d}t = \int_I \langle f,\varphi\rangle\,\mathrm{d}t
\]
holds for all $\varphi \in L^p(I,W^{1,p}_\sigma(\Omega))$.
\end{definition}

If we studied also the case $p\in(1,2)$, we would have to consider only 
test functions from the space of smooth functions. It is well known that 
the weak solution exists and is unique. It could be easily proven using 
the monotone operator theory. See for example \cite[Chapter 5]{mnrr} 
for periodic boundary conditions.
Now we formulate the main results of this paper.

\begin{theorem}\label{thm1}
Let $\Omega \subset \mathbb{R}^2$ be a bounded non-circular $\mathcal{C}^{3}$ 
domain and \eqref{as} holds for some $p\in[2,\infty)$. 
Let $u_0 \in W^{2+\beta,2}(\Omega)$ for $\beta\in(0,1/4)$, 
$\operatorname{div} u_0=0$, $f\in L^{\infty}(I,L^{q_0}(\Omega))$ and 
$\partial_t f \in L^{q_0}(I,W^{-1,{q_0'}}_{\sigma}(\Omega))$ for some $q_0>2$.
Then there exists a unique solution $(u,\pi)$ of \eqref{ns} with \eqref{bcps}, 
such that for some $\alpha>0$
$$
\nabla u,\pi \in \mathcal C^{0,\alpha}(\overline Q).
$$
\end{theorem}

\begin{remark} \rm
Perfect slip boundary conditions \eqref{bcps} are, as well as homogeneous 
Dirichlet boundary conditions, the limit case of  partial slip boundary conditions which are are also often called Navier's slip boundary conditions:
$$
u\cdot\nu=0,\quad \alpha [\mathcal{S}(Du)\nu]\cdot\tau + (1-\alpha)u_\tau = 0
\quad \alpha \in[0,1] \textrm{ on }\partial \Omega.
$$
It would be very interesting to obtain the same result as in Theorem \ref{thm1} 
also for the Navier's boundary condition. In several parts of the proof of 
Theorem \ref{thm1} we apply results from \cite{kt} that are formulated only 
for perfect slip boundary conditions. We don't know how to generalize these 
results also for partial slip boundary conditions.
\end{remark}

The paper is organized as follows: Section \ref{prel} contains preliminaries 
needed later, in Section \ref{lp} we gather $L^q$ theory results for 
the classical Stokes system. Further we extend $L^q$ theory results 
to generalized Stokes system where the Laplace operator is replaced by a 
general elliptic operator in divergence form with bounded measurable coefficients.

Section \ref{qp} is devoted to the proof of the main theorem in the case 
of quadratic growth, i.e. $p=2$. In Section~\ref{section5} we introduce 
the quadratic approximation of the stress tensor $\mathcal{S}(Du)$ which 
is done by the truncation of the generalized viscosity from above; i.e.,
 $\mu^\varepsilon(|D u^{\varepsilon}|):=\min\{\mu(|Du|),1/\varepsilon\}$ for $\varepsilon \in (0,1)$. 
We prove the main result for the approximated problem and pass from the 
approximated problem to the original one at the end.

\section{Preliminary general material}\label{prel}

\subsection{Function spaces}

Let $E$ and $F$ be reflexive Banach spaces. Although it is not necessary 
to have reflexive spaces in all definitions, for convenience we assume it. 
By $\mathcal L(E,F)$ we mean the Banach space of all bounded linear operators 
from $E$ to $F$ and $\mathcal{L}(E):=\mathcal{L}(E,E)$. 
If  $E$ is a linear subspace of $F$ and the natural injection $i: x \mapsto x$ 
belongs to $\mathcal L (E,F)$, we write $E\hookrightarrow F$. In the case $E$ 
is also dense in $F$, it will be denoted by $E\overset{d}\hookrightarrow F$. 
Furthermore, $\mathcal L is(E,F)$ consists of all topological linear isomorphisms 
from $E$ onto $F$. We also write $E\doteq F$ if $E\hookrightarrow F$ and 
$F\hookrightarrow E$, i.e. $E$ equals $F$ with equivalent norms.

A Banach space $E$ is said to be of class $\mathcal{HT}$, if the Hilbert 
transform is bounded on $L^p(\mathbb{R},E)$ for some (and then for all) 
$p\in(1,\infty)$. Here the Hilbert transform $H$ of a function 
$f\in \mathcal S(\mathbb R,E)$, the Schwartz space of rapidly decreasing 
$E$-valued functions, is defined by $Hf:=\frac{1}{\pi} PV(\frac{1}{t})*f$. 
It is well known theorem that the set of Banach spaces of class $\mathcal{HT}$ 
coincides with the class of $UMD$ spaces, where the $UMD$ stands for the 
property of unconditional martingale differences. Note that all closed 
subspaces of $L^q(\Omega)$ are $UMD$ spaces provided $q\in(1,\infty)$.

\subsection{Semigroups and interpolation-extrapolation scales}

For a linear operator $A$ in $E_0$ we denote the domain of $A$ by $\mathcal D(A)$. 
$A\in\mathcal{H}(E_1,E_0)$ means that $A$ is the negative infinitesimal generator 
of a bounded analytic semigroup in $E_0$ and $E_1\doteq\mathcal D(A)$. It holds
$$
\mathcal H(E_1,E_0)=\cup_{\kappa\ge 1,\,\omega>0}\mathcal{H}(E_1,E_0,\kappa,\omega),
$$
where $A\in \mathcal{H}(E_1,E_0,\kappa,\omega)$ if 
$\omega +A \in \mathcal Lis(E_1,E_0)$ and
$$
\kappa^{-1} \le \frac{\|(\lambda+A)u\|_{E_0}}{|\lambda|\|u\|_{E_0} 
+ \|u\|_{E_1}} \le \kappa,\quad Re(\lambda) \ge \omega,\quad u\in E_1.
$$

By $\sigma(A)$ we mean the spectrum of $A$ and $\varrho(A)$ denotes the 
resolvent set. A linear operator $A$ in $E$ is said to be of positive type 
if it belongs to $\mathcal P(E):=\cup_{K>1}P_K(E)$. $A \in P_K(E)$ if it is closed, 
densely defined, $\mathbb{R}^+\subset \varrho(-A)$ and 
$(1+s)\|(s+A)^{-1}\|_{\mathcal{L}(E)}\le K$ for $s\in\mathbb{R}^+$, where $K\ge 1$.

We say that a linear operator $A$ in $E$ is of type $(E,K,\vartheta)$, 
denoted by $A\in\mathcal{P}(E,K,\vartheta)$, if it is densely defined and if
$$
\Sigma_\vartheta:=\{|\arg z|\le\vartheta\}\cup\{0\}
\subset \varrho(-A)\quad\text{and}\quad  
(1+|\lambda|)\|(\lambda+A)^{-1}\|_{\mathcal{L}(E)}\le K, 
\quad \lambda\in\Sigma_\vartheta.
$$
Put $\mathcal P(E,\vartheta):=\cup_{K>1}\mathcal{P}(E,K,\vartheta)$.


A linear operator $A$ in $E$ is said to have bounded imaginary powers, 
in symbols,
$$
A \in \mathcal{BIP}(E):=\cup_{K\ge 1,\,\theta\ge 0}\mathcal{BIP}(E,K,\theta),
$$
provided $A\in\mathcal P(E)$ and there exist $\theta\ge 0$ and 
$K\ge 1$ such that $A^{is} \in \mathcal{L}(E)$ and 
$\|A^{is}\|_{\mathcal{L}(E)}\le Ke^{\theta |s|}$ for $s\in \mathbb R$.

We introduce an interpolation-extrapolation scale which is essential in the 
proof of Theorem \ref{lptheory}. Let $p,q\in(1,\infty)$, $\theta\in(0,1)$ 
and $[\cdot,\cdot]_\theta$ denotes the complex and $(\cdot,\cdot)_{\theta,q}$ 
the real interpolation functor. Let $A\in \mathcal H(E_1,E_0)$. 
Then we denote by $[(E_\alpha,A_\alpha);\alpha\in \mathbb R]$ 
the interpolation-extrapolation scale generated by $(E,A)$ and 
$[\cdot,\cdot]_\theta$ or $(\cdot,\cdot)_{\theta,q}$, where we set 
$E_k:=\mathcal D(A^k)$ for $k\in\mathbb N$ with $k\ge 2$. Also set $E^\sharp:=E'$ 
and $A^\sharp:=A'$, where $A'$ is the dual of $A$ in $E$ in the sense of 
unbounded linear operators. Finally let $E_k^\sharp:=\mathcal D((A^\sharp)^k)$ 
for $k\in\mathbb N$. Then we define $E_{-k}$ for $k\in\mathbb N$ by 
$E_{-k}:=(E_k^\sharp)'$. We put $E_{k+\theta}:=[E_k,E_{k+1}]_{\theta}$ 
(and similarly for the real interpolation functor). If $\alpha\ge 0$ we 
denote by $A_\alpha$ the maximal restriction of $A$ to $E_\alpha$ whose domain 
equals $\{u\in E_\alpha\cap E_1;\, Au\in E_\alpha\}$. 
If $\alpha<0$ then $A_\alpha$ is the closure of $A$ in $E_\alpha$.

For the dual interpolation functor $(\cdot,\cdot)^\sharp_\theta$ 
(which is equal to $[\cdot,\cdot]_\theta$ for the complex interpolation and 
$(\cdot,\cdot)_{\theta,q'}$ for real interpolation) we abbreviate the 
interpolation-extrapolation scale generated by $(E^\sharp,A^\sharp)$ and 
$(\cdot,\cdot)_\theta^\sharp$, by 
$[(E_\alpha^\sharp,A_\alpha^\sharp);\alpha\in \mathbb R]$ and call it 
interpolation-extrapolation scale dual to 
$[(E_\alpha,A_\alpha);\alpha\in \mathbb R]$. It holds 
$(E_{-\alpha})'\doteq E_\alpha^\sharp$ and $(A_{-\alpha})'=A^\sharp_\alpha$. 
For more details see \cite[Section V.2]{amann}.

\section{$L^q$ theory for Stokes system}\label{lp}

In this section we collect facts about $L^q$ theory for the Stokes system
\begin{equation}\label{stokes}
\begin{gathered}
\partial_t u-\Delta u + \nabla \pi = f,\quad \operatorname{div} u = 0\quad 
\text{in }Q,\\
u(0,\cdot)=u_0 \quad \text{on }\Omega,
\end{gathered}
\end{equation}
equipped with the perfect slip boundary conditions
\begin{equation}\label{sps}
u\cdot\nu=0,\quad [(Du)\nu]\cdot\tau=0\quad \text{on }I\times\partial\Omega.
\end{equation}

Unlike the main theorem of this paper which is formulated for 
$\Omega \subset \mathbb R^n$, $n=2$, results of this sections are valid 
for $n \ge 2$. Let $P$ denote the projection operator from 
$L^q(\Omega)$ to $L^q_\sigma(\Omega)$ associated with the Helmholtz decomposition. 
By $Bu=0$ we mean that \eqref{sps} holds in the sense of traces. 
Using the projection $P$ we shall define the Stokes operator $\mathbb A$ 
by $\mathbb Au = -P\Delta u$ for $u\in\mathcal{D}(\mathbb A)$, where
$$ 
\mathcal D(\mathbb A) = L^q_\sigma(\Omega) \cap H^2_{q,B}(\Omega),\quad 
H^2_{q,B}(\Omega):=\{u\in H^2_q(\Omega),\,Bu=0,\text{ on }\partial\Omega\}.
$$

Applying the Helmholtz projection $P$ to \eqref{stokes} with \eqref{sps},
 we eliminate the pressure from equations and with the help of the newly 
established notation the Stokes system reduces to
\begin{equation}\label{abstracteq}
\begin{gathered}
\partial_t u + \mathbb Au = Pf,\quad \operatorname{div} u =0\quad \text{in }Q,\\
u(0,\cdot)=u_0 \quad \text{on }\Omega,\quad Bu=0 \quad\text{on }I\times\partial\Omega.
\end{gathered}
\end{equation}

At first we mention some basic properties of the Stokes operator 
$\mathbb A$. From \cite{ss} we know that 
$\mathbb A\in\mathcal H(L^q_\sigma(\Omega) \cap H^2_{q,B}(\Omega),
L^q_\sigma(\Omega))$. This also tells us that
 $\mathbb A\in\mathcal P (L^q_\sigma(\Omega),\omega)$ for $\omega\in [0,\pi/2)$ 
(see \cite[Theorem II.4.6]{en}).  Shimada later showed in \cite{shimada} 
the $L^q$-maximal regularity for $\mathbb A$. In \cite[Theorem 1]{at}
Abels and Terasawa proved the following result.

\begin{proposition}\label{abels}
Let $q\in(1,\infty)$, $n\ge 2$, $r\in(n,\infty]$ such that $q,q' \le r$. 
Let $\Omega \subset \mathbb R^n$ be a domain with $W^{2-\frac{1}{r},r}$-boundary 
and $\vartheta \in (0,\pi)$. Then there is some $R>0$ such that 
$(\lambda + \mathbb A)^{-1}$ exists and
$$
(1+|\lambda|)\|(\lambda + \mathbb A)^{-1}\|_{\mathcal L(L^q(\Omega))} \le C
$$
for all $\lambda \in \Sigma_{\vartheta}$ with $|\lambda|\ge R$. Moreover,
$$
\big\| \int_{\Gamma_R} h(-\lambda)(\lambda + \mathbb A)^{-1}\,\mathrm{d}\lambda
\big\|_{\mathcal L(L^q(\Omega))} \le C\|h\|_{L^\infty(\Sigma_{\pi-\vartheta})}
$$
for every $h \in H^\infty(\vartheta)$, where $\Gamma=\partial\Sigma_\vartheta$, 
$\Gamma_R = \Gamma \setminus \overline{B_R(0)}$ and $H^\infty(\vartheta)$ 
denotes the Banach algebra of all bounded holomorphic functions 
$h:\Sigma_{\pi-\vartheta} \to \mathbb C$. In particular, for every
$\omega \in \mathbb R$ and $\vartheta' \in(0,\vartheta]$ such that 
$\omega + \Sigma_{\vartheta '} \subset \varrho(-\mathbb A)$ the shifted Stokes 
operator $\omega + \mathbb A$ admits a bounded $H^{\infty}$-calculus with respect 
to $\vartheta'$; i.e.,
$$
h(\omega + \mathbb A):= \frac{1}{2\pi i}\int_\Gamma h(-\lambda)(\lambda 
+ \omega + \mathbb A)^{-1}\,\mathrm{d}\lambda
$$
is a bounded operator satisfying
$$
\|h(\omega + \mathbb A)\|_{\mathcal L(L^q(\Omega))} \le C\|h\|_{L^\infty(\Sigma_{\pi-\vartheta})}
$$
for all $h\in H^\infty(\vartheta')$.
\end{proposition}

Note that the class of operators with a bounded $H^{\infty}$-calculus 
is a subclass of the operators which have $\mathcal{BIP}$, therefore 
these operators admit all important properties which has operators 
with bounded imaginary powers. For another properties of a bounded 
$H^{\infty}$-calculus we refer for example to \cite[Section 2, Subsection~2.4]{dhp}.

From the result of  Shibata and  Shimada in \cite{ss} follows that 
$\omega + \Sigma_{\vartheta '} \subset \varrho(-\mathbb A)$ even 
for $\omega=0$ provided the domain $\Omega$ is bounded and non-axisymmetric
 (see Definition \ref{axi}). Thus, Proposition~\ref{abels} and
 \cite[Theorem 1.3]{ss} gives $\mathbb A \in\mathcal{BIP}$. 
The Stokes operator $\mathbb A$ has realizations $\mathbb A_\alpha$ on 
$\mathbb E_\alpha$ for some $\alpha$. Concretely, from \cite[Section~2.2]{steiger} 
we know that $\mathbb A_\alpha\in \mathcal H(\mathbb E_{\alpha+1},\mathbb E_\alpha)$ 
for $\alpha\ge-1$.  Steiger in \cite{steiger} provides the characterization 
of spaces $\mathbb E_\alpha$:

\begin{proposition}[{\cite[Corollary 2.6]{steiger}}]\label{cor26}
Set $s_\alpha:=\{-2+1/q,-1+1/q,1/q,1+1/q\}$ and $F^s_q(\Omega):=H^s_p(\Omega)$ 
for the complex interpolation functor and $F^s_q(\Omega):= B^s_{q,q}(\Omega)$ 
for the real interpolation functor. Define
\begin{equation}\label{Hspbs}
F^s_{q,B}(\Omega):= \begin{cases}
\{u\in F^s_q(\Omega),\,Bu=0\text{ on }\partial\Omega\}, &s\in(1+1/q,2], \\
\{u\in F^s_q(\Omega),\,u\cdot\nu=0\text{ on }\partial\Omega\}, &s\in[1/q,1+1/q),\\
F^s_q(\Omega), &s\in[0,1/q),\\
\big(F^{-s}_{q',B,\sigma}(\Omega)\big)', &s \in[-2,0)\setminus s_\alpha
\end{cases}
\end{equation}
and
\begin{equation}
F^s_{q,B,\sigma}(\Omega):= \begin{cases}
F^s_{q,B}(\Omega)\cap L^q_\sigma(\Omega), &s\in[0,2]\setminus s_\alpha, \\
\big(F^{-s}_{q',B,\sigma}(\Omega)\big)', &s \in[-2,0)\setminus s_\alpha.
\end{cases}
\end{equation}
Then $\mathbb E_\alpha \doteq F^{2\alpha}_{q,B,\sigma}(\Omega)$ for 
$2\alpha \in [-2,2]\setminus s_\alpha$.
\end{proposition}

This gives
\begin{equation}\label{semigrupa}
A_\alpha\in\mathcal H(F^{2\alpha+2}_{q,B,\sigma}
(\Omega),F^{2\alpha}_{q,B,\sigma}(\Omega)), \quad 
2\alpha \in [-2,2]\setminus s_\alpha.
\end{equation}

\begin{remark}[{\cite[Remark 2.3c]{steiger}}]\label{rem.stei} \rm
The Helmholtz projection $P$ enjoys following continuity properties:
\begin{equation}\label{helm.cont}
P \in \mathcal L(F^s_{q,B}(\Omega))\cap \mathcal 
L(F^s_{q,B}(\Omega),F^s_{q,B,\sigma}(\Omega)), \quad 
s\in(-1+1/q,1+1/q)\setminus s_\alpha.
\end{equation}
\end{remark}

We will use the fact, that the property of bounded imaginary powers can be 
carried over the interpolation-extrapolation scales.

\begin{proposition}[{\cite[Proposition V.1.5.5]{amann}}]\label{1.5.5}
Let $A\in \mathcal P(E)$ and let $[(E_\alpha,A_\alpha); \alpha \in(-n,\infty)]$ 
be the interpolation-extrapolation scale generated by $(E,A)$ and an exact functor. 
If $A \in \mathcal{BIP}(E,M,\sigma)$ then 
$A_\alpha \in\mathcal{BIP}(E_\alpha,M,\sigma)$.
\end{proposition}

The reiteration property will be needed.

\begin{proposition}[{\cite[Theorem V.1.5.4]{amann}}] \label{1.5.4}
Suppose that $A\in \mathcal{BIP}(E)$. Then the interpolation-extrapolation
 scale $[(E_\alpha,A_\alpha); \alpha \in[-n,\infty)]$ generated by $(E,A)$ 
and complex interpolation functor possesses the reiteration property
$$
[E_\alpha, E_\beta]_\eta\doteq E_{(1-\eta)\alpha+\eta\beta}, \quad 
-n\le\alpha\le\beta<\infty,\quad \eta \in(0,1).
$$
\end{proposition}

Let us define the maximal $L^q$-regularity for the operator $A$ 
(compare 
\cite[Section III.1, Subsection 1.5 and Section III.4, Remark 4.10.9.c]{amann})

\begin{definition} \rm
Let $A\in\mathcal{H}(E_1,E_0)$ and $q\in(1,\infty)$. We say that the pair
$\big(L^q(I,E_0), L^q(I,E_1)\cap W^{1,q}(I,E_0)\big)$ is a pair of maximal
 regularity for $A$ (or $A$ has maximal regularity), if for 
$u_0\in (E_0,E_1)_{1-1/q,q}$ and $f\in L^q(I,E_0)$ there exists a unique 
solution $u\in L^q(I,E_1)\cap W^{1,q}(I,E_0)$ of \eqref{abstracteq}, and
\begin{equation}\label{maxr}
\|\partial_t u\|_{L^q(I,E_0)} + \|u\|_{L^q(I,E_0)} + \|Au\|_{L^q(I,E_0)} 
\le C\Big(\|f\|_{L^q(I,E_0)} + \|u_0\|_{(E_0,E_1)_{1-1/q,q}}\Big).
\end{equation}
\end{definition}

Further we mention the relation between maximal regularity and negative 
infinitesimal generators of a bounded analytic semigroup.


\begin{proposition}[{\cite[Theorem III.4.10.7]{amann}}]\label{4.10.8}
Suppose that $E_0$ is a UMD space, $A\in \mathcal H (E_1,E_0)$ and there 
are constants $M>0$, $\vartheta \in(0,\pi/2)$ such that 
$\Sigma_\vartheta \subset \varrho(-A)$ and for $\lambda \in\Sigma_\vartheta$ 
and $j=0,1$ holds
\[
\|A\|_{\mathcal L(E_1,E_0)} + (1+|\lambda|)^{1-j}\|(\lambda+A)^{-1}
\|_{\mathcal L(E_0,E_j)} \le M 
\]
and suppose that there exist constants $N\ge 1$ and $\theta \in[0,\pi/2)$ 
such that $A\in\mathcal{BIP}(E_0,N,\theta)$. Then $A$ has maximal regularity 
and the estimate \eqref{maxr} holds uniformly with respect to $T$.
\end{proposition}

To specify the shape of the domain $\Omega$ we add the definition of 
axisymmetric domain in the same way as in \cite[Definition-Lemma 1]{dv}.

\begin{definition}\label{axi} \rm
Let $\Omega$ be a smooth bounded open subset of $\mathbb{R}^n$, $n\ge 2$. 
We say that $\Omega$ is axisymmetric if and only if there exists a nontrivial
 rigid motion $R$ which is tangent to $\partial\Omega$; or equivalently, 
which satisfies for all $t\in \mathbb{R}$ $e^{tR}\Omega=\Omega$. 
Here $e^{tR}$ is the isometry defined via $\frac{d}{dt}e^{tR}(x)=Re^{tR}(x)$.
\end{definition}

By rigid motions $R$ we understand affine maps $R: \Omega \to \mathbb{R}^n$
whose linear part is antisymmetric. If we consider the most common dimensions 
$n=2$ and $n=3$ we can use simpler definition. A domain in $\mathbb{R}^2$ 
is axisymmetric if it has a circular symmetry around some point. 
A domain in $\mathbb{R}^3$ is axisymmetric if it admits an axis of symmetry; i.e.,
the domain is preserved by a rotation of arbitrary angle around this axis. 
If the domain admits two nonparallel axes of symmetry, then it is spherically 
symmetric around some point.

The main result of this section is the following.


\begin{theorem}\label{lptheory}
Let $\Omega \subset \mathbb{R}^n$ be a bounded non-axisymmetric $\mathcal{C}^{2,1}$ 
domain, $q\in [2,\infty)$, $f\in L^q(I,W_\sigma^{-1,q'}(\Omega))$, 
$u_0 \in B^{1-2/q}_{q,q,B,\sigma}(\Omega)$ then there exists a constant $C>0$ 
and the unique weak solution of \eqref{abstracteq} satisfying
$$ 
\|\nabla u\|_{L^q(Q)} + \|u\|_{BUC(I,B^{1-2/q}_{q,q,B,\sigma}(\Omega))} 
\le C\Big(\|f\|_{L^q(I,W_\sigma^{-1,q'}(\Omega))} 
+ \|u_0\|_{B_{q,q,B,\sigma}^{1-2/q}(\Omega)}\Big).
$$ 
The constant $C$ is independent of $T,u,f$ and $u_0$.
\end{theorem}

\begin{proof}
 We consider the system \eqref{abstracteq} instead of \eqref{stokes} 
with \eqref{sps}.
Since for $UMD$ space $E$, $E'$ is one as well and for an interpolation 
couple of $UMD$ spaces the interpolation spaces are also $UMD$ 
(see \cite[Theorem III.4.5.2]{amann}), $\mathbb E_{-1/2}$ is a $UMD$ space. 
Proposition \ref{1.5.5} gives us $\mathbb A_{-1/2}$ has $\mathcal{BIP}$. 
Together with \eqref{semigrupa}, \cite[Corollary I.1.4.3]{amann} and 
\cite[Theorem~1.3]{ss} we can see that assumptions of Proposition \ref{4.10.8} 
are fulfilled for $\mathbb A_{-1/2}$. Therefore we obtain \eqref{maxr} 
for $\mathbb A_{-1/2}$ and $E_0=\mathbb E_{-1/2}$:
\begin{equation}\label{tojeono}
\begin{aligned}
&\|\partial_t u\|_{L^q(I,\mathbb E_{-1/2})} + \|u\|_{L^q(I,\mathbb E_{-1/2})}
 + \|\mathbb A_{-1/2}u\|_{L^q(I,\mathbb E_{-1/2})} \\
&\le C \Big(\|f\|_{L^q(I,\mathbb E_{-1/2})} 
+ \|u_0\|_{(\mathbb E_{-1/2},\mathbb E_{1/2})_{1-1/q,q}}\Big).
\end{aligned}
\end{equation}

It remains to determine the correct spaces in \eqref{tojeono}. 
For the space of initial condition $u_0$ we get by Proposition~\ref{cor26} 
for the complex interpolation functor
$$
u_0 \in (H^{-1}_{q,B,\sigma}(\Omega), H^{1}_{q,B,\sigma}(\Omega))_{1-1/q,q}.
$$
This space equals (with equivalent norms) to $B^{1-2/q}_{q,q,B,\sigma}(\Omega)$ 
since for $q\ge 2$,
\begin{equation}\label{int.ic}
\begin{split}
B^{1-2/q}_{q,q,B,\sigma}(\Omega) 
&\doteq (L^q_\sigma(\Omega), H^1_{q,B,\sigma}(\Omega))_{1-2/q,q} \\
&\doteq ([H^{-1}_{q,B,\sigma}(\Omega),H^{1}_{q,B,\sigma}(\Omega)]_{1/2},
  H^{1}_{q,B,\sigma}(\Omega))_{1-2/q,q}\\
&\doteq (H^{-1}_{q,B,\sigma}(\Omega),H^{1}_{q,B,\sigma}(\Omega))_{1-1/q,q},
\end{split}
\end{equation}
where we used Proposition~\ref{1.5.4}. The similar interpolation of the 
solenoidal functions in case of Dirichlet boundary conditions is done 
in \cite[Proof of Lemma 9.1]{pamann}. From the embedding 
\cite[Theorem V.4.10.2]{amann}
$$
L^q(I,E_1)\cap W^{1,q}(I,E_0) \hookrightarrow BUC(I,(E_0,E_1)_{1-1/q,q}),
$$
we obtain $u \in BUC(I, B^{1-2/q}_{q,q,B,\sigma}(\Omega))$. 
Due to $\|\mathbb A_{-1/2}u\|_{\mathbb E_{-1/2}}=\|u\|_{\mathbb E_{1/2}}$ 
and $\mathbb E_{1/2} \doteq W^{1,q}_\sigma(\Omega)$ we have boundedness of 
$\nabla u$ in $L^q(Q)$. It remains to find the space for $f$. 
By Proposition~\ref{cor26},
$$
f\in L^q(I,W^{-1,q'}_\sigma(\Omega)),
$$
since $H^s_{q}(\Omega)\doteq W^{s,q}(\Omega)$ for $s\in \mathbb{Z}$.
\end{proof}


Without loss of generality we may assume that there exists a symmetric 
tensor $G\in L^q(Q)$, such that the weak formulation of the right 
hand side of \eqref{stokes} can be written in the form
\begin{equation}\label{rhs}
\int_Q G:D\varphi\,\mathrm{d}x\,\mathrm{d}t = \int_I\langle f,\varphi \rangle \,\mathrm{d}t \quad 
\forall \varphi \in L^{q'}(I,W^{1,q'}_\sigma(\Omega)).
\end{equation}
To prove it, we proceed in the same way like in 
\cite[Proof of Proposition~2.1, Step~1]{kms3} where the authors are dealing with 
periodic boundary conditions. Consider the Stokes system which can be 
formulated in the weak form for a. a. $t\in I$ as follows
\begin{equation}\label{agi}
\int_\Omega Dw(t):D\varphi \,\mathrm{d}x=\langle f(t),\varphi \rangle \quad 
\forall \varphi \in W^{1,q}_\sigma(\Omega).
\end{equation}
As $f\in L^q(I,W^{-1,q}_\sigma(\Omega))$, there exists a solution 
$w(t) \in W^{1,q}_\sigma(\Omega)$ of \eqref{agi} enjoying the estimate
$$
\|w(t)\|_{W^{1,q}(\Omega)} \le C\|f\|_{W^{-1,q}_\sigma(\Omega)}
$$
with the positive constant $C$ independent of $t$. Consequently, 
$w \in L^q(I,W^{1,q}_\sigma(\Omega))$ and
$$
\|w\|_{L^q(I,W^{1,q}(\Omega))} \le C\|f\|_{L^q(I,W^{-1,q}_\sigma(\Omega))}.
$$
Defining $G=Dw$ we conclude \eqref{rhs} from \eqref{agi} by density arguments. 
Therefore, for all $f\in L^q(I,W^{-1,q}_\sigma(\Omega))$ there exists
 $G\in L^q(Q)$ such that \eqref{rhs} and  estimate
$$
\|G\|_{L^q(Q)} \le C\|f\|_{L^q(I,W^{-1,q}_\sigma(\Omega))}
$$
holds. We would like to point out that the perfect slip boundary conditions 
are hidden in the weak formulation. If $G$ is smooth enough then it holds
$$
\int_I\langle f,\varphi \rangle \,\mathrm{d}t = -\int_Q \operatorname{div}
 G \cdot \varphi \,\mathrm{d}x\,\mathrm{d}t + \int_I\int_{\partial \Omega} 
(G\nu)\tau(\varphi\cdot\tau)\,\mathrm{d}\sigma\,\mathrm{d}t \quad \forall
 \varphi \in L^q(I,W^{1,q}_\sigma(\Omega)).
$$
The Stokes system \eqref{stokes} with \eqref{sps} can be formulated in 
the weak form as follows
\begin{equation}\label{stokesweak}
\int_I \langle \partial_t u,\varphi\rangle \,\mathrm{d}t + \int_Q Du: D\varphi \,\mathrm{d}x\,\mathrm{d}t  
= \int_Q  G: D\varphi\,\mathrm{d}x\,\mathrm{d}t\quad \forall \varphi \in L^{q}(I,W^{1,q}_\sigma(\Omega)).
\end{equation}


Introducing the solution operator $S: (G,u_0) \mapsto Du$, 
we conclude first from the existence theory, that $S$ is continuous 
from $L^2(Q)\times L^2_\sigma(\Omega)$ to $L^2(Q)$ with the norm less 
or equal to 1. By Theorem \ref{lptheory} we know that $S$ is continuous 
from $L^{q_1}(Q)\times B^{1-2/{q_1}}_{q_1,q_1,B,\sigma}(\Omega)$ 
to $L^{q_1}(Q)$ with norm estimated by $C_q>1$. 
Since $S(G,u_0)=S(G,0) + S(0,u_0)$, Riesz-Thorin theorem and the real 
interpolation method implies following assertion, see for 
example \cite[Theorem~5.2.1 and Theorem~6.4.5]{bl}.

\begin{lemma}\label{riesz-thorin}
Let $\Omega$ be a bounded non-axisymmetric $\mathcal{C}^{2,1}$ domain 
and $q_1>2$. There exist constant $C>0$ and $K:=C_{q_1}^{q_1/(q_1-2)}$ 
such that for every $q\in (2,q_1)$, arbitrary $G \in L^q(I,L^{q}_\sigma(\Omega))$, 
$u_0 \in B^{1-2/q}_{q,q,B,\sigma}(\Omega)$ there exists a unique solution $u$ 
of \eqref{stokesweak} satisfying
\[
\|Du\|_{L^q(Q)} \le K^{1-\frac{2}{q}}\Big(\|G\|_{L^q(Q)} 
+ C\|u_0\|_{B^{1-2/q}_{q,q,B,\sigma}(\Omega)}\Big).
\]
\end{lemma}

For $q>2$ small enough Lemma \ref{riesz-thorin} allows us to prove the 
$L^q$ theory for a generalized Stokes system, where the Stokes operator 
is replaced by a general elliptic operator with bounded measurable coefficients. 
More precisely, let $0<\gamma_1\le\gamma_2$ and suppose that the coefficient 
matrix $\mathbb{M} \in L^{\infty}(Q)$ is symmetric in the sense 
$M^{kl}_{ij} = M^{ij}_{kl}=M^{ji}_{kl}$ for $i,j,k,l=1,2$ and fulfils for 
all $B\in\mathbb{R}^{2\times 2}$, $x\in \Omega$ and $t\in I$,
$$
\gamma_1|B|^2 \le \mathbb{M}(t,x):B\otimes B \le \gamma_2 |B|^2.
$$

We consider the  system
\begin{equation}\label{gstokes}
\begin{split}
&\int_I \langle \partial_t u,\varphi\rangle \,\mathrm{d}t 
+ \int_Q \mathbb{M}:Du\otimes D\varphi \,\mathrm{d}x\,\mathrm{d}t  \\
&= \int_Q G: D\varphi\,\mathrm{d}x\,\mathrm{d}t\quad \forall \varphi 
\in L^{q}(I,W^{1,q}_\sigma(\Omega)).
\end{split}
\end{equation}
The following lemma states the $L^q$ theory result.

\begin{lemma}\label{lptheorygs}
Let $\Omega$ be a bounded non-axisymmetric $\mathcal{C}^{2,1}$ domain and $q>2$. 
There exist constants $K,L>0$ such that if $q\in[2,2+L\frac{\gamma_1}{\gamma_2})$,
 $G\in L^q(Q)$ and $u_0 \in B^{1-2/q}_{q,q,B,\sigma}(\Omega)$ then the unique
 weak solution $u\in L^q(I,W^{1,q}_\sigma(\Omega))$ of \eqref{gstokes} satisfies
\[
\|Du\|_{L^q(Q)} + \gamma_2^{-\frac{1}{q}}\|u\|_{BUC(I,B^{1-2/q}_{q,q,B,\sigma}
(\Omega))} \le \frac{K}{\gamma_1}\Big(\|G\|_{L^q(Q)} 
+ \gamma_2^{1-\frac{1}{q}}\|u_0\|_{B^{1-2/q}_{q,q,B,\sigma}(\Omega)}\Big).
\]
\end{lemma}

\begin{proof} 
We omit the proof. It can be found in \cite[Proposition 2.1]{kms2} 
for periodic boundary conditions or in \cite[Proposition 2.1]{kap} 
for homogeneous Dirichlet boundary conditions. The only generalization 
consists of including perfect slip boundary conditions. $L^q$ theory result 
for classical Stokes system with perfect slip boundary conditions is needed, 
but it is shown in Lemma \ref{riesz-thorin}.
\end{proof}

We also use the $L^q$ theory for stationary variant of the system \eqref{gstokes}. 
For symmetric coefficient matrix $\mathbb{M} \in L^\infty(\Omega)$ fulfilling 
for all $B\in\mathbb{R}^{2\times 2}$ and $x\in \Omega$, 
$\gamma_1|B|^2 \le \mathbb{M}(x):B\otimes B \le \gamma_2 |B|^2$, 
$0<\gamma_1\le \gamma_2$ we sutdy the problem
\begin{equation}\label{stacgstokes}
\int_{\Omega}\mathbb{M}:Du\otimes D\varphi\,\mathrm{d}x = \int_{\Omega} G:D\varphi\,\mathrm{d}x \quad \forall 
\varphi \in W_\sigma^{1,q}(\Omega).
\end{equation}

\begin{lemma}\label{staclptheory}
Let $\Omega$ be a bounded non-axisymmetric $\mathcal{C}^{2,1}$ domain. 
Then there are constants $K,L>0$ such that if
 $q\in[2,2+L\frac{\gamma_1}{\gamma_2})$ and $G \in L^{q}(\Omega)$, 
then the unique weak solution of \eqref{stacgstokes} satisfies
$$
\|D u\|_{L^q(\Omega)} \le \frac{K}{\gamma_1}\|G\|_{L^q(\Omega)}.
$$
\end{lemma}

\begin{proof} 
See \cite[Lemma 2.6]{kms2} for no slip boundary conditions. For perfect slip 
boundary conditions we would proceed analogically.
\end{proof}

\section{Proof of the main results for the quadratic potential}\label{qp}

In this section we prove Theorem \ref{thm1} for $p=2$.
\smallskip

\noindent\textbf{Step 1.}
In this step we obtain a priori estimates from the existence theory. 
For $f\in W^{1,2}(I,W^{-1,2}_\sigma(\Omega))$ with $f(0) \in L^2(\Omega)$ 
and $u_0\in W^{2,2}(\Omega)\cap W^{1,2}_\sigma(\Omega)$ 
we know the existence of a unique weak solution of \eqref{ns} with 
\eqref{bcps} fulfilling
\begin{equation} \label{zexistence}
\begin{gathered}
u \in L^\infty(I,L^2(\Omega)) \cap L^2(I,W^{1,2}_\sigma(\Omega)), \\
 \partial_t u \in L^\infty(I,L^2(\Omega)) \cap L^2(I,W^{1,2}_\sigma(\Omega)),
\quad \pi\in L^2(I,L^2(\Omega)).
\end{gathered}
\end{equation}
It can be shown using Galerkin approximation.
Let $\{\omega^k\}_{k=1}^\infty$ be the orthogonal basis of $L^2_\sigma(\Omega)$
and $W^{1,2}_{\sigma}(\Omega)$ consisting of eigenvectors of the Stokes operator
with perfect slip boundary conditions. Such basis can be easily constructed
provided $\Omega$ is non-circular domain. Set
 $H^n = \operatorname{span} \{\omega_1,\dots,\omega^N\}$ and define
the continuous projection $P^N: L^2_\sigma(\Omega) \to H^N$ as follows:
$$
P^N u = \sum_{k=1}^N (u,\omega^k)\omega^k.
$$
Define $u^N(t,x)= \sum_{k=1}^N c_k^N(t)\omega^k$ where $c_k^N(t)$ solves
the Galerkin system
\begin{equation}\label{gs}
\begin{gathered}
\langle \partial_t u^N(t),\omega^k\rangle + \int_{\Omega}\mathcal{S}(Du^N):D(\omega^k)\,\mathrm{d}x
+ \int_{\Omega} (u^n\otimes u^n):\nabla\omega^k\,\mathrm{d}x = \langle f,w^k\rangle,\\
u^N(0)=u_0^N=P^N u_0, \quad 1\le k \le N.
\end{gathered}
\end{equation}
After multiplying the Galerkin system \eqref{gs} by $c^N_k(t)$,
summing up, using Gronwall's and Korn's inequalities we derive the
following a priori estimate,
$$
\sup_{t\in I}\|u^N(t)\|_2^2 + \int_I\|u^N(\tau)\|_{1,2}^2\,\mathrm{d}\tau \le C.
$$

Further we apply the time derivative to \eqref{gs}, multiply it by 
$\partial_t c^N_k(t)$ and sum up. Unlike the previous apriori estimates, 
before using Gronwall's inequality, the boundedness of $\|\partial_t u^N(0)\|_2^2$ 
needs to be shown. This can be done easily, since 
$P^N: W^{2,2}(\Omega)\cap W^{1,2}_\sigma(\Omega) \to H^N$ is bounded uniformly
 with respect to $N$ (c. f. \cite[Lemma 4.26]{mnrr}), we can use \eqref{gs}. 
Thus, after Gronwall's inequality we have
$$
\sup_{t\in I}\|\partial_t u^N(t)\|_2^2 + \int_I\|\partial_t u^N(\tau)\|_{1,2}^2\,\mathrm{d}\tau 
\le C.
$$
Passing to the limit with $N \to \infty$ (where we use the Aubin-Lions' lemma 
to obtain the strong convergence of $u^N$ in $L^2(I,L^4(\Omega))$ 
and Minty's trick to identify the limit of $\mathcal{S}(Du^N)$ with $\mathcal{S}(Du)$) 
we get the first two relations in \eqref{zexistence}.

Since $\partial_t u$, $\operatorname{div} \mathcal{S}(Du)$,
 $\operatorname{div} (u\otimes u)$ and $f$ lie in $L^2(I,W^{-1,2}_\sigma(\Omega))$, 
we can reconstruct the pressure $\pi$ at almost every time level via 
De Rham's theorem and Ne\v{c}as' theorem on negative norms and obtain 
$\pi \in L^2(\Omega)$ for almost every $t\in I$.
\smallskip

\noindent\textbf{Step 2.} We improve the regularity in space.
If we additionally assume $f \in L^\infty(I,L^2(\Omega))$ we are able to show that
\begin{equation}\label{sec4step2}
u \in L^\infty(I,W^{2,2}(\Omega)), \quad \pi\in L^\infty(I,W^{1,2}(\Omega)).
\end{equation}
From Step 1 we know that $\partial_t u$ is regular enough in order to move 
it to the right hand side of $\eqref{ns}_1$. At almost every 
time level $t\in I$ we can use the stationary theory. Boundary regularity 
in tangent direction is based on the difference quotient technique. 
In normal direction near the boundary the main tools are the operator 
curl and Ne\v{c}as' theorem on negative norms. See for example 
\cite[Section 3]{mnr} for homogeneous Dirichlet boundary conditions. 
The information about the pressure comes from the fact that the right hand 
side of $\nabla \pi = f +\operatorname{div} \mathcal{S} 
- \operatorname{div}(u\otimes u)-\partial_t u$ is in $L^2(\Omega)$ 
for a. a. $t\in I$. Adding the assumption $\int_{\Omega} \pi\,\mathrm{d}x=0$ we get by Poincar\'e 
inequality the existence of $\pi \in W^{1,2}(\Omega)$ at almost every time 
level $t\in I$ together with a bound independent of $t$.
\smallskip

\noindent\textbf{Step 3.}
We improve the regularity in time using $L^q$ theory for Stokes system.
If we moreover suppose that $f \in L^{q_1}(I,W_\sigma^{-1,q_1'}(\Omega))$ 
for some $q_1>2$ and $u_0\in W^{2+\beta,2}(\Omega)$ for $\beta\in(0,1/4)$ 
we are able to prove the existence of $q_2>2$ such that the unique weak 
solution satisfies for all $q\in(2,q_2)$

\begin{equation}\label{regtime}
\partial_t u \in L^q(I,W^{1,q}_\sigma(\Omega)) \cap BUC(I,B^{1-2/q}_{q,q,B,\sigma}
(\Omega)).
\end{equation}

Denoting $w:=\partial_t u$ and $\tau:= \partial _t \pi$ in the sense of 
distributions, we observe from \eqref{ns} that $(w,\tau)$ solves
\begin{equation}\label{timederns}
\int_I\langle \partial_t w,\varphi\rangle\,\mathrm{d}t +\int_Q \partial_{Du}^2\Phi(|Du|):
Dw\otimes D\varphi\,\mathrm{d}x\,\mathrm{d}t = \int_I\langle\partial_t (f-(u\cdot\nabla)u),
\varphi\rangle\,\mathrm{d}t,
\end{equation}
for all $\varphi\in L^q(I,W^{1,q}_\sigma(\Omega))$.
It is easy to see that $\partial_t (u\cdot\nabla u) \in L^s(I,W^{-1,s}(\Omega))$ 
for all $s\in[1,4]$.

To obtain \eqref{regtime} as a result of application of Lemma \ref{lptheorygs} 
for the system \eqref{timederns} we need to ensure that 
$\|\partial_t u(0)\|_{B^{1-2/q}_{q,q,B,\sigma}(\Omega)}$ is bounded. 
Let $\beta\in (0,1/4)$ and $\varphi \in W^{-\beta,2}(\Omega)$ with 
$\|\varphi\|_{W^{-\beta,2}(\Omega)} \le 1$ be arbitrary. 
We recall that the Helmholtz projection $P$ enjoys the continuity properties 
as mentioned in Remark \ref{rem.stei}. Thus,
\begin{equation}\label{cas.der.v0}
\begin{split}
|\langle\partial_tu(0),\varphi\rangle| 
&= |\langle\partial_tu(0),P\varphi\rangle| \\
&\leq |\langle\operatorname{div} S(Du_0) + (u_0\cdot\nabla)u_0-f(0),P\varphi\rangle| \\
&\leq C(\|u_0\|_{W^{2+\beta,2}(\Omega)} + \|u_0\|_{W^{2,2}(\Omega)}^2 
 + \|f(0)\|_{W^{\beta,2}(\Omega)})\le C.
\end{split}
\end{equation}
Since $W^{\beta,2}(\Omega) \hookrightarrow B^{1-2/q}_{q,q}(\Omega)$ if
 $q$ is close enough to $2$ we obtain 
$\|\partial_t u(0)\|_{B^{1-2/q}_{q,q,B,\sigma}(\Omega)} \le C$ for all 
$q\in(2,q_2)$ where $q_2$ is sufficiently close to 2.
\smallskip

\noindent\textbf{Step 4.}
We show that $u \in L^\infty(I,W^{2,q}(\Omega))$ due to the stationary theory.
The previous step shows us that $\partial_t u \in L^\infty(I,L^q(\Omega))$ 
for some $q>2$. Therefore we are able to move $\partial_t u$ to the right 
hand side of $\eqref{ns}_1$ and apply the result \cite[Theorem 3]{kms1} 
for $p=2$ which tells us that there exists a positive $\varepsilon$, 
such that $u\in W^{2,2+\varepsilon}(\Omega)$ and 
$\pi \in W^{1,2+\varepsilon}(\Omega)$ for \eqref{ns} with perfect slip 
boundary conditions.
\smallskip

\noindent\textbf{Step 5.}
We improve the regularity of $\pi$ in time.
There exists a $q>2$ such that for all $s\in(0,\frac{1}{2})$
$$
\pi \in W^{s,q}(I,L^q(\Omega)).
$$
We closely follow the proof of \cite[Lemma 3.4]{kap}. For a function 
$g(t)$ defined on the time interval $I$ and $(t_1,t_2) \subset I$ 
set $\delta^t g:= g(t_2)-g(t_1)$. The idea of the proof is based on 
subtracting the equation $\eqref{ns}_1$ in the time $t_2$ from the same 
equation in time $t_1$ which leads to
\begin{equation}\label{diference}
\int_{\Omega} \delta^t \pi \operatorname{div} \varphi \,\mathrm{d}x
= \int_{\Omega} [\delta^t(\partial_t u - f)\varphi - \delta^t(u\otimes u 
- \mathcal{S}(Du))D\varphi]\,\mathrm{d}x,
\end{equation}
which holds for all $\varphi \in W^{1,2}(\Omega)$ with $\varphi\cdot \nu=0$ 
on $\partial \Omega$. From \eqref{sec4step2} and \eqref{regtime} one may 
easily show the existence of $q>2$ and $s\in(0,1/2)$ such that 
$u\in W^{s,q}(I,W^{1,q}(\Omega))$ and $\partial_t u \in W^{s,q}(I,L^q(\Omega))$. 
Together with the assumptions on the right hand side $f$ we can notice that 
\eqref{diference} holds also for all $\varphi \in W^{1,q'}(\Omega)$ with 
$\varphi = 0$ at $\partial \Omega$. Consider the problem
\begin{equation}\label{bog.pi}
\begin{gathered}
\operatorname{div} \varphi^t = \delta^t \pi |\delta^t \pi|^{q-2} - \frac{1}{|\Omega|}\int_{\Omega} \delta^t \pi|\delta^t\pi|^{q-2}\,\mathrm{d}x \quad \mbox{in }\Omega,\\
\varphi^t = 0\quad \text{on }\partial\Omega.
\end{gathered}
\end{equation}
The right hand side of \eqref{bog.pi} has zero mean value over $\Omega$ 
and belongs to $L^{q'}(\Omega)$ due to \eqref{sec4step2}, therefore 
Bogovski\u{\i}'s Lemma (for the formulation and proof 
c.f. \cite[Lemma 3.3]{amrouchegirault}) guaranties the existence of $\varphi^t$ 
satisfying the estimate $\|\varphi^t\|_{1,q'} \le C\|\delta^t \pi\|_q^{q-1}$. 
Taking $\varphi^t$ as a test function in \eqref{diference} leads to
\begin{equation}\label{nejakyodhad}
\|\delta^t \pi\|_q^q \le \varepsilon \|\delta^t \pi\|_q^q 
+ C_\varepsilon(\|\delta^t\partial_t u\|_q^q + \|\delta^t f\|_{-1,q}^q 
+ \|\delta^t\nabla u\|_q^q).
\end{equation}
Dividing \eqref{nejakyodhad} by $|t_2-t_1|^{1+sq}$ and integrating twice over 
$I$ gives
$$
\|\pi\|^q_{W^{s,q}(I,L^q(\Omega))} 
= \int_I\int_I \frac{\|\delta^t \pi\|_q^q}{|t_2-t_1|^{1+sq}}\,\mathrm{d}t_1\,\mathrm{d}t_2 \le C,
$$
which completes the proof.
\smallskip

\noindent\textbf{Step 6.}
We summarize the result of this section and uses imbedding theorems 
to complete the proof.
Up to now we have shown
$$
u \in L^\infty(I,W^{2,q}(\Omega)) \cap W^{1,q}(I,L^q(\Omega)), 
\quad \pi \in L^\infty(I,W^{1,q}(\Omega))\cap W^{s,q}(I,L^q(\Omega)).
$$
As we are in two dimensions, $q>2$, $s\in(\frac{1}{q},\frac{1}{2})$, 
following imbeddings hold
\begin{gather}
L^\infty(I,W^{1,q}(\Omega)) \hookrightarrow 
 L^\infty(I,\mathcal{C}^{0,1-\frac{2}{q}}(\overline{\Omega})),\label{im1}\\
W^{1,q}(I,L^q(\Omega)) \hookrightarrow \mathcal{C}^{1-\frac{1}{q}}
 (\overline{I},L^q(\Omega)),\label{im2}\\
W^{s,q}(I,L^q(\Omega)) \hookrightarrow \mathcal{C}^{s-\frac{1}{q}}
 (\overline{I},L^q(\Omega)).\label{im3}
\end{gather}

Now we are ready to apply the following lemma.

\begin{lemma}[{\cite[Lemma 2.6]{kap}}] \label{imbd}
Let $\Omega \subset \mathbb{R}^2$ be a bounded $\mathcal{C}^2$ domain. 
Let $f\in L^\infty(I,\mathcal{C}^{0,\alpha}(\overline{\Omega}))$ and 
$f \in \mathcal{C}^{0,\beta}(\overline{I},L^s(\Omega))$ for some 
$\alpha, \beta \in (0,1)$ and $s>1$. Then 
$f\in \mathcal{C}^{0,\gamma}(\overline{Q})$ with 
$\gamma = \min\{\alpha, \frac{\alpha\beta s}{\alpha s + 2}\}$.
\end{lemma}

Using \eqref{im1} and \eqref{im2} together with Lemma \ref{imbd} we 
obtain $\nabla u \in \mathcal{C}^{0,\alpha}(\overline Q)$ for certain $\alpha>0$.
 \eqref{im1}, \eqref{im3} with Lemma \ref{imbd} gives us
$\pi \in \mathcal{C}^{0,\alpha}(\overline Q)$ for some $\alpha>0$, 
which concludes the proof of main results for $p=2$.


\section{Proof of the main results for the super-quadratic potential}\label{section5}

In this section we prove Theorem \ref{thm1} for $p>2$. 
The proof consists of several steps.
\smallskip

\noindent\textbf{Step 1.}
We introduces quadratic approximations. 
In a similar way as in \cite{kt} we are concerned with the regularized problem
\begin{equation}\label{regns}
\begin{gathered}
\partial_t  u^{\varepsilon} -\operatorname{div} \mathcal{S}^{\varepsilon}(D u^{\varepsilon}) + ( u^{\varepsilon}\cdot\nabla) u^{\varepsilon} + \nabla \pi^{\varepsilon}
= f,\quad \operatorname{div}  u^{\varepsilon}=0\quad \text{in }\,Q,\\
 u^{\varepsilon}(0,\cdot)=u_0\quad \text{in }\Omega,
\end{gathered}
\end{equation}
where we consider quadratic approximation $\mathcal S^\varepsilon$ of 
$\mathcal S$ defined for $\varepsilon \in (0,1)$ by the truncation of the 
viscosity $\mu$ from above,
\begin{equation}\label{aproximace}
\mu^\varepsilon(|D u^{\varepsilon}|):=\min\Big\{\mu(|Du|),\frac{1}{\varepsilon}\Big\},\quad 
\mathcal{S}^{\varepsilon}(D u^{\varepsilon}):= \mu^\varepsilon(|D u^{\varepsilon}|)D u^{\varepsilon}.
\end{equation}
Scalar potential $\Phi^{\varepsilon}$ to $\mathcal{S}^{\varepsilon}(D u^{\varepsilon})$ can be constructed in the following way
$$
\Phi^{\varepsilon}(s):=\int_0^s\mu^\varepsilon(t)t\,\mathrm{d}t
$$
and satisfies growth conditions \eqref{as} for $p=2$, i.e. there exists 
$C_1>0$ and $C(\varepsilon)$ such that for all 
$A,B\in\mathbb{R}^{2\times 2}_{\rm sym}$
\begin{equation}\label{reg.rust.1}
C_1|B|^2 \le \partial^2_A\Phi^{\varepsilon}(|A|):B\otimes B \le C(\varepsilon)|B|^2.
\end{equation}

The approximation \eqref{aproximace} guarantees that for a fixed 
$\varepsilon\in (0,1)$ the results of the previous section holds 
for $ u^{\varepsilon}$ and $\pi^{\varepsilon}$ solving \eqref{regns} equipped with the perfect 
slip boundary conditions.
\smallskip

\noindent\textbf{Step 2.}
We present growth conditions dependent on $\varepsilon$. 
Due to the results of the previous section we are able to use techniques 
which enable us to gain uniform estimates with respect to $\varepsilon$. 
At first we need a growth estimates of $\Phi^{\varepsilon}$ with precise dependence 
on $\varepsilon$. In other words, the constant $C(\varepsilon)$ 
in the estimate \eqref{reg.rust.1} needs to be specified. 
To this purpose we define the function $\vartheta_{\varepsilon}$ by 
$\vartheta_{\varepsilon}(s):=\min\{(1+s^2)^{\frac{1}{2}},\frac{1}{\varepsilon}\}$.
Now, there exist constants $0<C_3\le C_4$ such that for all 
$\varepsilon \in (0,1)$ and $A,B\in\mathbb{R}^{2\times 2}_{\rm sym}$
\begin{equation}\label{reg.rust.2}
C_3\vartheta_{\varepsilon}(|A|)^{p-2}|B|^2 \le \partial^2_A\Phi^{\varepsilon}(|A|):B\otimes B 
\le C_4 \vartheta_{\varepsilon}(|A|)^{p-2}|B|^2.
\end{equation}
As a corollary of \eqref{reg.rust.2}, the following estimates can be derived 
(see \cite[Lemma 2.22]{mnr} for the proof.)
\begin{gather}
C\vartheta_{\varepsilon}(|A|)^{p-2}|A|^2\le \mathcal{S}^{\varepsilon}(A):A,\label{45}\\
C|\mathcal{S}^{\varepsilon}(A)|\le \vartheta_{\varepsilon}(|A|)^{p-2}|A|.\label{46}
\end{gather}
The lower estimate in \eqref{45} can be done independent of $\varepsilon$, 
since \eqref{reg.rust.1} holds:
\begin{equation}\label{44}
C_5|A|^2 \le\mathcal{S}^{\varepsilon}(A):A.
\end{equation}
At this point we would like to emphasize that from now all 
constants in following steps are independent of $\varepsilon$.
\smallskip

\noindent\textbf{Step 3.}
We provide  $L^{\infty}(I,L^2(\Omega))\cap L^2(I,W^{1,2}(\Omega))$ estimates 
of $ u^{\varepsilon}$ and $\partial_t  u^{\varepsilon}$.
We recall estimates from the previous section which hold also for the 
approximated problem since the lower bound in \eqref{44} is independent 
on $\varepsilon$.
\begin{gather}
\| u^{\varepsilon}\|_{L^{\infty}(I,L^2(\Omega))} + \|\nabla  u^{\varepsilon}\|_{L^{2}(Q)} 
 \le C,\label{kint1} \\
\|\partial_t u^{\varepsilon}\|^2_{L^{\infty}(I,L^2(\Omega))} 
 + \|\nabla \partial_t u^{\varepsilon}\|_{L^{2}(Q)} \le C.\label{kint}
\end{gather}
The relation \eqref{kint1} is an a priori estimate obtained by taking 
solution as a test function (at the level of Galerkin approximation). 
Roughly speaking, the estimate \eqref{kint} is performed by taking 
time derivative of the equation \eqref{regns} and testing by time 
derivative of $ u^{\varepsilon}$. More precisely, it is not applied directly to 
the equation \eqref{regns}, but still to the Galerkin system. 
To estimate the time derivative of the Galerkin approximation of $ u^{\varepsilon}$ 
at the time $t=0$ we proceed in the same way like in \eqref{cas.der.v0}.

Note that \eqref{kint1} and \eqref{kint} give $ u^{\varepsilon} \in L^\infty(I,W^{1,2}(\Omega))$,
\begin{align*}
\|\nabla u^{\varepsilon}(s,\cdot)\|_2^2 - \|\nabla u^{\varepsilon}(0,\cdot)\|_2^2 
&= \int_{\Omega}\int_0^s\partial_t|\nabla u^{\varepsilon}(t,\cdot)|^2\,\mathrm{d}t\,\mathrm{d}x \\
&\le 2\|\nabla u^{\varepsilon}\|_{L^2(Q)}\|\partial_t\nabla  u^{\varepsilon}\|_{L^2(Q)} \le C.
\end{align*}

\noindent\textbf{Step 4.} We escribe the boundary $\partial \Omega$.
To discuss boundary regularity in following steps, we need a suitable 
description of the boundary $\partial\Omega$. 
Let us denote $x=(x_1,x_2)$. We suppose that $\Omega \in \mathcal{C}^3$, 
therefore there exists $c_0>0$ such that for all $a_0>0$ there exists $n_0$ 
points $P\in \partial\Omega$, $r>0$ and open smooth set 
$\Omega_0 \subset \subset \Omega$ that we have
\[
\Omega \subset \Omega_0 \cup \bigcup_P B_{r}(P)
\]
and for each point $P\in\partial \Omega$ there exists local system of coordinates 
for which $P=0$ and the boundary $\partial \Omega$ is locally described by 
$\mathcal{C}^3$ mapping $a_P$ that for $x_1 \in(-3r,3r)$ fulfils
\begin{gather*}
x \in \partial \Omega \Leftrightarrow x_2=a_P(x_1), \quad 
B_{3r}(P)\cap\Omega = \{x\in B_r(P)\,\textrm{and}\,x_2>a_P(x_1)\}=:\Omega_{3r}^P,
\\
\partial_1 a_P(0)=0, \quad |\partial_1 a_P(x_1)|\le a_0,\quad 
|\partial_1^2 a_P(x_1)|+|\partial_1^3 a_P(x_1)|\le c_0.
\end{gather*}
Point $P$ can be divided into $k$ groups such that in each group 
$\Omega_{3r}^P$ are disjoint and $k$ depends only on dimension $n$. 
Let the cut-off function $\xi_P(x) \in \mathcal{C}^{\infty}(B_{3r}(P))$ 
and reaches values
\[\label{serezavacka}
\xi_P(x)  \begin{cases}
=1  &x \in B_r(P), \\
\in (0,1) &x \in B_{2r}(P)\setminus B_r(P),\quad  \\
=0 &x\in \mathbb{R}^2\setminus B_{2r}(P).
\end{cases}
\]
Next, we assume that we work in the coordinate system corresponding to $P$. 
Particularly, $P=0$. Let us fix $P$ and drop for simplicity the index $P$. 
The tangent vector and the outer normal vector to $\partial\Omega$ are defined as
$$
\tau=\big(1,\partial_1 a(x_1)\big), \quad \nu=\big(\partial_1a(x_1),-1\big),
$$
tangent and normal derivatives as
$$
\partial_\tau = \partial_1 + \partial_1 a(x_1)\partial_2, \quad 
\partial_\nu = -\partial_2 + \partial_1 a(x_1)\partial_1.
$$
\smallskip

\noindent\textbf{Step 5.}
We show that $ u^{\varepsilon} \in L^{\infty}(I,W^{2,2}(\Omega))$ uniformly in
 $\varepsilon \in (0,1)$.
From Step 3 we obtained that $\partial_t  u^{\varepsilon} \in L^\infty(I,L^2(\Omega))$, 
therefore we can fix $t\in I$, move $\partial_t  u^{\varepsilon}$ to the right hand side 
of \eqref{regns} and at almost every time level consider the stationary problem
\begin{equation}\label{stacregns}
\begin{gathered}
-\operatorname{div} \mathcal{S}^{\varepsilon}(D u^{\varepsilon}) + ( u^{\varepsilon}\cdot\nabla) u^{\varepsilon} + \nabla \pi^{\varepsilon} 
= h,\quad \operatorname{div}  u^{\varepsilon}=0\quad \text{in }\Omega,\\
 u^{\varepsilon}\cdot\nu=0,\quad [\mathcal{S}^{\varepsilon}(D u^{\varepsilon})\nu]\cdot\tau =0\quad \text{on }\partial\Omega,
\end{gathered}
\end{equation}
where $h:=f-\partial_t  u^{\varepsilon} \in L^2(\Omega)$. Previous section provides 
$ u^{\varepsilon} \in W^{2,2}(\Omega)$, $\mathcal{S}^{\varepsilon}(D u^{\varepsilon})\in W^{1,2}(\Omega)$ and 
$\pi^{\varepsilon}\in W^{1,2}(\Omega)$. Thus we can multiply \eqref{stacregns} 
by a suitable test function which is at least in $L^2(\Omega)$ and integrate 
over $\Omega$. We focus only on the boundary regularity and work in the 
local system of coordinates.
Following \cite[Lemma 4.2, Remark 4.9]{kt} we choose as a test function 
$\varphi=(\varphi_1,\varphi_2)$,
\begin{gather*}
\varphi=(\partial_2[\Theta-\partial_\tau( u^{\varepsilon}\cdot\nu)\xi^2], 
\partial_1[-\Theta + \partial_\tau( u^{\varepsilon}\cdot\nu)\xi^2]),\\
\Theta:= \partial_\nu( u^{\varepsilon}\cdot\tau)\xi^2 -  u^{\varepsilon}\cdot(\partial_\nu\tau + \partial_\tau\nu)\xi^2.
\end{gather*}
This test function is constructed to get rid of the pressure $\pi^{\varepsilon}$ 
and to obtain optimal information from the elliptic term. 
These most difficult estimates, in which we extract from 
$-\int_{\Omega}\operatorname{div} \mathcal{S}^{\varepsilon}(D u^{\varepsilon})\cdot\varphi\,\mathrm{d}x$ boundedness of the term
 $\int_{\Omega}\mu^\varepsilon(|D u^{\varepsilon}|)|\nabla^2 u^{\varepsilon}|^2\,\mathrm{d}x$, are done in \cite[Proof of Theorem 1.7]{kt}, 
therefore we omit the calculations. It remains to estimate the convective 
term and the right hand side of \eqref{stacregns}. After long, but elementary 
calculations we are able to show that
\begin{equation}\label{con.te}
|\int_{\Omega} ( u^{\varepsilon}\cdot\nabla) u^{\varepsilon}\cdot\varphi\,\mathrm{d}x| \le C\int_{\Omega} (| u^{\varepsilon}||\nabla u^{\varepsilon}|^2 
+ | u^{\varepsilon}|^2|\nabla u^{\varepsilon}|)\,\mathrm{d}x,
\end{equation}
where we used the divergence-free constraint and the properties of the test 
function $\varphi$. Using H\"older and Young inequalities, 
$\|\cdot\|_4^2\le C\|\cdot\|_{1,2}\|\cdot\|_2$  and the information 
$ u^{\varepsilon} \in L^\infty(I,W^{1,2}(\Omega))$ we continue estimating \eqref{con.te}:
\[
C(\| u^{\varepsilon}\|_2\|\nabla u^{\varepsilon}\|_4^2 + \| u^{\varepsilon}\|_4^2\|\nabla u^{\varepsilon}\|_2) 
\le \varepsilon\|\nabla^2 u^{\varepsilon}\|_2^2 + C\|u\|_{1,2}^2 
+ C\|\nabla u^{\varepsilon}\|_2^2\| u^{\varepsilon}\|_2^2.
\]
The last estimate is easy.
\[
\big|\int_{\Omega} h\cdot\varphi\,\mathrm{d}x\big|
\le \int_{\Omega} |h|(|\nabla^2 u^{\varepsilon}| + |\nabla u^{\varepsilon}| + | u^{\varepsilon}|)\,\mathrm{d}x 
\le C\|h\|_2^2 + \varepsilon\|\nabla^2 u^{\varepsilon}\|_2^2 + C\|u\|_{1,2}^2.
\]
Since $\mu^\varepsilon (|D u^{\varepsilon}|)>1$ and $\varepsilon>0$ can be chosen arbitrarily small, 
we obtain
\begin{equation}\label{udd}
\|\nabla^2 u^{\varepsilon}\|_2^2 \le \int_{\Omega}\mu^\varepsilon(|D u^{\varepsilon}|)|\nabla^2 u^{\varepsilon}|^2\,\mathrm{d}x \le C,
\end{equation}
where $C$ does not depend on $\varepsilon$ and $t\in I$, therefore we have
\begin{equation}\label{wdvadva}
 u^{\varepsilon} \in L^{\infty}(I,W^{2,2}(\Omega)).
\end{equation}
\smallskip

\noindent\textbf{Step 6.}
We improve the information about $\partial_t  u^{\varepsilon}$. 
In the same spirit as in Step 3 from the previous section we denote 
$w:=\partial_t u$ and $\tau:= \partial _t \pi$ in the sense of distributions, 
which solves \eqref{timederns} where $\Phi$ is replaced by $\Phi^{\varepsilon}$. 
The right hand side of \eqref{timederns} is bounded uniformly with respect 
to $\varepsilon \in (0,1)$ in $L^{q_0}(I,W^{-1,{q'_0}}_\sigma(\Omega))$ 
for some $q_0>2$, since from \eqref{kint1}, \eqref{kint} and \eqref{wdvadva}
 we have $\partial_t[( u^{\varepsilon}\cdot\nabla) u^{\varepsilon}]\in L^s(I,W^{-1,s}(\Omega))$ for all 
$s\in[1,4]$.

Set $V_\varepsilon:=\sup_Q|\vartheta_{\varepsilon}(|D u^{\varepsilon}|)|$. From \eqref{reg.rust.2} we have 
for all $t\in I$, $x\in\Omega$, for all $\varepsilon \in (0,1)$ and 
$A,B\in\mathbb{R}^{2\times 2}_{\rm sym}$
$$
c|B|^2 \le \partial^2_A\Phi^{\varepsilon}(|A|):B\otimes B \le C V_\varepsilon^{p-2}(|A|)|B|^2.
$$
From Lemma \ref{lptheorygs} we have the existence of positive constants $K$ 
and $L$  such that for all $q\in (2,q_2]$, where $q_2:= 2+L/V_\varepsilon^{p-2}$ 
holds
\begin{equation}\label{est.v1}
\begin{aligned}
&\|\nabla w\|_{L^{q}(Q)} + V_\varepsilon^{\frac{2-p}{q}}
\|w\|_{BUC(I,B^{1-2/q}_{q,q,B,\sigma}(\Omega))} \\
&\le K \Big(\|f\|_{L^q(I,W^{-1,q'}(\Omega))}
 + V_{\varepsilon}^{(p-2)(1-1/q)}\|\partial_tu_0\|_{B^{1-2/q}_{q,q,B,\sigma}
(\Omega)}\Big).
\end{aligned}
\end{equation}
Without loss of generality we may assume that $q_2<q_0$. Thus, after estimating 
last norm on the right hand side of \eqref{est.v1} in the same way like in 
Step 3 in Section 3 we have
\[
\|\partial_t u^{\varepsilon}\|_{BUC(I,B^{1-2/q}_{q,q,B,\sigma}(\Omega))} 
\le C\Big(V_\varepsilon^{\frac{p-2}{q}} +V_\varepsilon^{p-2}\Big) 
\le CV_\varepsilon^{p-2}.
\]
\smallskip

\noindent\textbf{Step 7.} We improve the information about $\nabla^2  u^{\varepsilon}$. 
In this step we obtain better space regularity.
 Up to now we have $\vartheta_{\varepsilon}\in L^{\infty}(I,W^{1,2}(\Omega))$. 
We are going to show that $\vartheta_{\varepsilon}\in L^{\infty}(I,W^{1,q}(\Omega))$ for some $q>2$.

We omit estimates of $\nabla^2 u^{\varepsilon}$ in the interior of $\Omega$ and we 
focus on estimates near the boundary. We start with the tangential direction. 
Localizing the problem, we work in $\Omega_{3r}^P$, where the boundary is 
locally described by the $\mathcal{C}^3$ mapping $a_p$. 
For simplicity we drop the index $P$.

We multiply \eqref{stacregns} by $-\partial_\tau\varphi\xi$, integrate over
 $\Omega_{3r}$ and after similar steps as in \cite[Lemma 4.6]{kt} 
we derive the identity
\begin{equation}\label{identita}
\begin{split}
&\int_{\Omega_{3r}} \partial_\tau\mathcal{S}^{\varepsilon}(D u^{\varepsilon}):D\varphi\xi\,\mathrm{d}x\\
& = -\int_{\Omega_{3r}} h\cdot\partial_\tau(\varphi\xi)\,\mathrm{d}x + \int_{\Omega}( u^{\varepsilon}\cdot\nabla) u^{\varepsilon}\partial_\tau\varphi\xi\,\mathrm{d}x \\
&\quad +\int_{\Omega_{3r}} \mathcal{S}^{\varepsilon}(D u^{\varepsilon}):\Big[\partial_\tau\varphi\otimes\nabla\xi-\nabla\varphi\partial_\tau\xi 
+ (\partial_1^2a,0)\otimes \partial_2\varphi\xi 
 + \nabla\Big(\varphi\cdot\partial_\tau\nu\frac{\nu}{|\nu|^2}\xi\Big)\Big]\,\mathrm{d}x  \\
&\quad +\int_{\Omega_{3r}} \operatorname{div} \mathcal{S}^{\varepsilon}(D u^{\varepsilon})\cdot[(\varphi\cdot\partial_\tau\nu)
 \frac{\nu}{|\nu|^2}\xi-\varphi\partial_\tau\xi]\,\mathrm{d}x \\
&\quad + \int_{\Omega_{3r}} \partial_1^2 a [h_2 + (\operatorname{div} \mathcal{S}^{\varepsilon}(D u^{\varepsilon}))_2
 -( u^{\varepsilon}\cdot\nabla u^{\varepsilon})_2]\varphi_1\xi\,\mathrm{d}x  \\
&\quad  +\int_{\Omega_{3r}}[h_1 +\partial_1ah_2 + (\operatorname{div}\mathcal{S}^{\varepsilon}(D u^{\varepsilon}))_1
  + \partial_1 a(\operatorname{div}\mathcal{S}^{\varepsilon}(D u^{\varepsilon}))_2+ ( u^{\varepsilon}\cdot\nabla u^{\varepsilon})_1\\
 &\quad + \partial_1a( u^{\varepsilon}\cdot\nabla u^{\varepsilon})_2]\varphi\nabla\xi\,\mathrm{d}x
\end{split}
\end{equation}
for all $\varphi \in W^{1,q'}_{\sigma}(\Omega)$, 
$\operatorname{supp} \varphi \subset \overline{\Omega_{3r}}$. Terms on the right hand side 
of \eqref{identita} comes at first from the fact that we add and subtract 
some lower order terms in order to let the boundary term vanish while 
integrating by parts. Second, tangent derivative does not commute with the 
gradient and we use 
$\nabla\partial_\tau\varphi = \partial_\tau\nabla\varphi + (\partial_1^2a,0)\otimes\partial_2\varphi$. 
Third, we use  \eqref{stacregns} and replace 
$\partial_2\pi^{\varepsilon}$ by
$h_2 + (\operatorname{div}\mathcal{S}^{\varepsilon}(D u^{\varepsilon}))_2 
+ ( u^{\varepsilon}\cdot\nabla u^{\varepsilon})_2$ 
and similarly  for $\partial_\tau\pi^{\varepsilon}$.

We denote $w:= \partial_\tau u^{\varepsilon}\xi -(0,\partial_1^2a u^{\varepsilon}_1)\xi + z$, where $z$ 
is the solution of
\begin{gather}
\operatorname{div} z = -\partial_\tau u^{\varepsilon}\cdot\nabla\xi-\partial_1^2a u^{\varepsilon}_1\partial_2\xi\quad
\text{in } \Omega_{3r},\label{bo1} \\
z = 0 \quad \text{on } \partial\Omega_{3r}.\label{bo2}
\end{gather}
The right hand side of \eqref{bo1} was obtained from 
$\operatorname{div}\big(-\partial_\tau u^{\varepsilon}\xi +(0,\partial_1^2a u^{\varepsilon}_1)\xi\big)$ 
using the fact that $\operatorname{div}  u^{\varepsilon}=0$. The role of $z$ is to ensure 
that $\operatorname{div} w=0$. On $\partial\Omega$ it holds $w\cdot\nu=0$ since
\[
w\cdot\nu = [\partial_\tau u^{\varepsilon}\cdot\nu +\partial_1^2a u^{\varepsilon}_1]\xi + z\cdot\nu 
= \partial_\tau( u^{\varepsilon}\cdot\nu)\xi= 0.
\]
Thus, the compatibility condition holds
\begin{align*}
\int_{\partial\Omega} z\cdot\nu \,\mathrm{d}\sigma
&= \int_{\Omega}\operatorname{div} z\,\mathrm{d}x 
=\int_{\Omega} \operatorname{div}(-\partial_\tau u^{\varepsilon}\xi +(0,\partial_1^2a u^{\varepsilon}_1)\xi)\,\mathrm{d}x \\
&=-\int_{\partial\Omega}\partial_\tau( u^{\varepsilon}\cdot\nu)\xi\,\mathrm{d}\sigma = 0
\end{align*}
and $z$ solving \eqref{bo1} and \eqref{bo2} exists by Bogovski\u {\i}'s 
Lemma and enjoys the estimate $\|z\|_{1,q}\le C\|\nabla u^{\varepsilon}\|_q$ for some $C>0$.

Using the definition of $w$, from \eqref{identita} we obtain
\[
\int_{\Omega} \partial_{D u^{\varepsilon}}\mathcal{S}^{\varepsilon}(D u^{\varepsilon}): Dw\otimes D\varphi\,\mathrm{d}x 
= \langle g,\varphi\rangle \quad \forall \varphi\in W^{1,q'}_{\sigma}(\Omega),
\]
with
\begin{align*}
\langle g,\varphi\rangle &= \textrm{RHS of \eqref{identita}} 
 + \int_{\Omega}\partial_{D u^{\varepsilon}}\mathcal{S}^{\varepsilon}(D u^{\varepsilon}):[Dz + \partial_\tau u^{\varepsilon}\otimes\nabla\xi \\
&\quad + (\partial_1^2a,0)\otimes_S\partial_2 u^{\varepsilon}\xi 
- D\big((0,\partial_1^2a,0)\xi\big)]D\varphi\,\mathrm{d}x.
\end{align*}
Due to the assumption on $f$ and results from Step 4 we have 
$\|g\|_{-1,{q_2'}} \le CV_\varepsilon^{p-2}$ and after application of 
Lemma \ref{staclptheory} we obtain
\begin{equation}\label{teclp}
\|\nabla\partial_\tau u^{\varepsilon}\xi\|_{L^q(\Omega)} \le CV_\varepsilon^{p-2}.
\end{equation}

We recall that $q$ depends on $\varepsilon$ by the relation 
$q\in(2,2+L/V_\varepsilon^{p-2}]$.
In order to control whole $\nabla^2 u^{\varepsilon}$ we need an estimate of type 
\eqref{teclp} in the normal direction which is locally $x_2$. 
Since $\partial_2^2 u^{\varepsilon}_2$ can be expressed from the condition 
$\operatorname{div} u^{\varepsilon}=0$, we focus on $\partial_2^2 u^{\varepsilon}_1$. 
Following \cite[Theorem 3.19]{kms2} we can extract the desired 
estimate from the equation \eqref{stacregns} after employment of the 
operator curl. Let us shorten $\mathcal{S}^{\varepsilon}(D u^{\varepsilon})$ to $\mathcal{S}^{\varepsilon}$ and $\vartheta_{\varepsilon}(|D u^{\varepsilon}|)$ 
to $\vartheta_{\varepsilon}$. Denoting $G:=\partial_2\mathcal{S}^{\varepsilon}_{12}$ we have due to \eqref{46} 
and \eqref{reg.rust.2},
\begin{gather*}
\|\xi G\|_{-1,q} \le \|\mathcal{S}^{\varepsilon}_{12}\|_q \le \|{\vartheta_{\varepsilon}}^{p-2} D u^{\varepsilon}\|_q,\notag\\
\|\partial_1(\xi G)\|_{-1,q} \le C\|{\vartheta_{\varepsilon}}^{p-2} D u^{\varepsilon}\|_q + C'\|{\vartheta_{\varepsilon}}^{p-2}\partial_1\nabla u^{\varepsilon}\|_q.\notag
\end{gather*}
From  \eqref{stacregns} after applying curl we have
\begin{align*}
\|\partial_2(\xi G)\|_{-1,q} 
&\le C(\|\partial_1(\mathcal{S}^{\varepsilon}_{21}+\mathcal{S}^{\varepsilon}_{22}-\mathcal{S}^{\varepsilon}_{11})\|_q 
 + \|f\|_q+ \| u^{\varepsilon}\cdot\nabla u^{\varepsilon}\|_q + \|\partial_t u^{\varepsilon}\|_q)  \\
&\le C\big\{\|{\vartheta_{\varepsilon}}^{p-2} D u^{\varepsilon}\|_q + \|{\vartheta_{\varepsilon}}^{p-2}\partial_1\nabla u^{\varepsilon}\|_q 
 + V_\varepsilon^{p-2}+1\big\}:=H.
\end{align*}
Ne\v{c}as' theorem on negative norms gives
$$
\|\xi G\|_q \le C(\|\xi G\|_{-1,q} + \|\nabla(\xi G)\|_{-1,q}) \le H.
$$
From the definition of $G$ and symmetry of $Du$ we obtain
$$
\partial_{12}\mathcal{S}^{\varepsilon}_{12}\partial_2 D u^{\varepsilon}_{12} 
= \frac{G}{2} - \frac{1}{2}\partial_{11}\mathcal{S}^{\varepsilon}_{12}\partial_2D u^{\varepsilon}_{11} 
- \frac{1}{2}\partial_{22}\mathcal{S}^{\varepsilon}_{12}\partial_2D u^{\varepsilon}_{22}.
$$
Using $\partial_{12}\mathcal{S}^{\varepsilon}_{12}\ge C{\vartheta_{\varepsilon}}^{p-2}$ and the condition 
$\operatorname{div}  u^{\varepsilon}=0$ we get that
$\|\xi{\vartheta_{\varepsilon}}^{p-2}\partial_2^2 u^{\varepsilon}_1\|_{q}\le H$. Thus,
\begin{equation}\label{25}
\|\xi\vartheta_{\varepsilon}^{p-2}\nabla^2 u^{\varepsilon}\|_{q} 
\le C\|\xi G\|_q + \|\xi\vartheta_{\varepsilon}^{p-2}\nabla\partial_\tau u^{\varepsilon}\|_q 
+ \tilde{C}\sup_{x_1\in(-3r,3r)}|\partial_1a|\|\xi\vartheta_{\varepsilon}^{p-2}\nabla^2 u^{\varepsilon}\|_{q},
\end{equation}
where $\tilde{C}$ is absolute constant. Since we can choose $r$ sufficiently 
small in order to 
$\tilde{C}\max_{P\in\partial\Omega}\sup_{x_1\in(-3r,3r)}|\partial_1a|\le 1/2$, 
the last term \eqref{25} can be absorbed into the left hand side. We have
\begin{equation}\label{thetapq}
\|\xi\vartheta_{\varepsilon}^{p-2}\nabla^2 u^{\varepsilon}\|_{q_2} \le CV_\varepsilon^{p-2}V_\varepsilon^{p-2}.
\end{equation}
From \eqref{udd} the boundedness of the term 
$\int_\Omega \mu^\varepsilon(|D u^{\varepsilon}|)|\nabla^2 u^{\varepsilon}|^2\,\mathrm{d}x$ is obtained, in other words
$\|\vartheta_{\varepsilon}^{\frac{p-2}{2}}\nabla^2 u^{\varepsilon}\|_2\le C$. 
Interpolation of this result with \eqref{thetapq} gives for $q\in(2,q_2)$
\begin{equation}\label{pointerpolaci}
\|\vartheta_{\varepsilon}^{\frac{p-2}{2}}\nabla^2 u^{\varepsilon}\|_{q} \le C V_\varepsilon^{\beta 2(p-2)},
\end{equation}
where $1/q = \beta/q_2 + (1-\beta)/2$.
Since it holds 
\[
\|\vartheta_{\varepsilon}^{p/2}\|_{1,q} \le \|\vartheta_{\varepsilon}^{p/2}\|_q 
+ \|\vartheta_{\varepsilon}^{\frac{p-2}{2}}\nabla^2 u^{\varepsilon}\|_q,
\]
we want to use the following lemma for $f=\vartheta_{\varepsilon}^{p/2}$.

\begin{lemma}\label{imbd2}
Let $\Omega\subset \mathbb{R}^2$ be a bounded $\mathcal{C}^2$ domain and 
$f\in W^{1,q}(\Omega)$ for some $q>2$. 
Then $f\in \mathcal{C}(\overline{\Omega})$ and there is $C>0$ independent 
of $q$ such that
\begin{equation}\label{ooo}
\sup_{\Omega}|f| \le C\Big(\frac{q-1}{q-2}\Big)^{1-1/q}\|f\|_{1,q}.
\end{equation}
\end{lemma}

\begin{proof} It follows from the proof of \cite[Theorem 2.4.1]{ziemer}. 
The result holds also for $\Omega\subset\mathbb{R}^n$, with $q>n$ and $q-n$ 
instead of $q-2$ in the denominator of \eqref{ooo}.
 \end{proof}

Because $\frac{q-1}{q-2} \le CV^{p-2}$, we obtain
\begin{equation}\label{v}
V_\varepsilon^{p/2}\le CV^{(p-2)(1-\frac{1}{q})}V_\varepsilon^{\beta2(p-2)}.
\end{equation}
Note that $(p-2)(1-1/q) \to p/2-1$ as $q\to 2$ and the exponent containing 
the interpolation parameter $\beta$ can be made arbitrarily small, 
therefore we can rewrite \eqref{v} as $V_\varepsilon \le \hat{C}$. 
This together with \eqref{pointerpolaci} gives
$$
\sup_{t\in I} \|\nabla^2  u^{\varepsilon}\|_q \le C.
$$
\smallskip

\noindent\textbf{Step 8.}
We pass from the regularized problem to the original problem. 
In the previous step we showed $V_\varepsilon \le \hat{C}$, where 
$V_\varepsilon=\sup_Q|\vartheta_\varepsilon(|D u^{\varepsilon}|)|$.
Since $\vartheta_{\varepsilon}(s)=\min\{(1+s^2)^{1/2},\frac{1}{\varepsilon}\} 
\le \frac{1}{\varepsilon}$, it is sufficient to choose $\varepsilon$ 
in order to have $\hat{C}\le \frac{1}{\varepsilon}$. 
Thus, $ u^{\varepsilon}=u$ is the solution of the original problem \eqref{ns} 
and it holds that $\sup_Q(1+|Du|^2)^{1/2}\le C$ which leads to 
$\sup_{t\in I} \|\nabla^2 u\|_q \le C$.

Since we passed from the regularized problem to the original one, 
the regularity of pressure $\pi$ which we proved in Section \ref{qp} 
for quadratic potential holds also for the super-quadratic case.
\hfill \qed

\subsection*{Acknowledgements} 
The work was supported by  project 281211/B-MAT/MFF from  
the Grant Agency of Charles University (GAUK) and by  project
 201/09/0917 from the Czech Science Foundation (GA\v{C}R).
The author is a junior researcher in the University Centre for
 Mathematical Modelling, Applied Analysis and Computational Mathematics 
(Math MAC). The author would like to thank Petr Kaplick\'y 
for productive discussions and for his support.

\begin{thebibliography}{99}

  \bibitem{at} Abels, H., Terasawa, Y.; 
\emph{On Stokes operators with variable viscosity in bounded and unbounded domains},
  Math. Ann.  \textbf{344}  (2009),  no. 2, 381-429.

  \bibitem{amann} Amann, H.; 
\emph{Linear and Quasilinear Parabolic Problems, Volume I: Abstract Linear Theory,}
 Birkh\"auser, Basel, 1995.

  \bibitem{pamann} Amann, H.; 
\emph{On the Strong Solvability of the Navier-Stokes Equations}, 
J. Math. Fluid Mech. \textbf{2} (2000), 16--98.

  \bibitem{arada1} Arada, N.; 
\emph{Optimal control of shear-thinning fluids,} SIAM J. Control Optim. 
\textbf{50} (2012), no. 4, 2515-2542.

   \bibitem{arada2} Arada, N.; 
\emph{Optimal control of shear-thickening flows,} 
SIAM J. Control Optim. \textbf{51} (2013), no. 3, 1940-1961.

  \bibitem{amrouchegirault} Amrouche, C., Girault, V.; 
\emph{Decomposition of vector spaces and application to the Stokes problem
 in arbitrary dimension}, Czechoslovak Math. J. \textbf{44} (1994), 
no. 119, 109--140.

\bibitem{bl} Bergh, J., L\"ofstr\"om, J.;
 \emph{Interpolation spaces, an introduction.} 
Springer-Verlag, Berlin, Heidelberg, New York, 1976, x + 204 pp.

  \bibitem{dhp} Denk, R., Hieber, M., Pr\"uss, J.; 
\emph{$R$-Boundedness, Fourier Multipliers and Problems of Elliptic and Parabolic 
Type}, Mem. Amer. Math. Soc.  \textbf{166}  (2003),  no. 788, viii+114 pp.

\bibitem{dv} Desvillettes, L., Villani, C.; 
\emph{On a variant of Korn's inequality arising in statistical mechanics,} 
ESIAM: Control, Optimisation and Calculus of Variations \textbf{8} (2002), 603--619.

\bibitem{die} Diening, L.; 
\emph{Theoretical and Numerical Results for Electrorheological Fluids}, 
Freiburg: Albert-Ludwigs-Universit\"at 2002, Ph.D. Thesis.

 \bibitem{en} Engel, K.-J., Nagel, R.; 
\emph{One-Parameter Semigroups of Linear Evolution Equations}, 
With contributions by S. Brendle, M. Campiti, T. Hahn, G. Metafune, G. Nickel, 
D. Pallara, C. Perazzoli, A. Rhandi, S. Romanelli and R. Schnaubelt. 
Graduate Texts in Mathematics, 194. Springer-Verlag, New York, 2000. xxii+586 pp.

  \bibitem{fs}Frese, J., Seregin, G.; 
\emph{Full regularity for a class of degenerated parabolic systems in two 
spatial variables} Manuscripta Math. \textbf{99} (1999), 517--539.

\bibitem{kap} Kaplick\'y, P.;
 \emph{Regularity of flows of a non-Newtonian fluid subject to 
Dirichlet boundary conditions}, Journal for Analysis and its Applications, 
\textbf{24}, no. 3 (2005), 467--486.

  \bibitem{kms1} Kaplick\'y, P., M\'alek, J. and Star\'a, J.; 
\emph{On Global existence of smooth two-dimensional steady flows for a class 
of non-Newtonian fluids under various boundary conditions,} 
Applied Nonlinear Analysis, New York, Kluwer/Plenum, 1999, 213--229.

  \bibitem{kms2}Kaplick\'y, P., M\'alek, J. and Star\'a, J.; 
\emph{$C^{1,\alpha}$-solutions to a class of nonlinear fluids in two dimensions 
- stationary Dirichlet problem,} Zap. Nauchn. Sem. S.-Peterburg. Otdel. 
Mat. Inst. Steklov. (POMI) \textbf{259} (1999), 89--121.

  \bibitem{kms3}Kaplick\'y, P., M\'alek, J. and Star\'a, J.; 
\emph{Global-in-time H\"older continuity of the velocity gradients for fluids 
with shear-dependent viscosities,} NoDEA Nonlinear Differential Equations Appl. 
\textbf{9} (2002), 175--195.

  \bibitem{kp} Kaplick\'y, P., Pra\v{z}\'ak, D.; 
\emph{Differentiability of the solution operator and the dimension of 
the attractor for certain power-law fluids,} 
J. Math. Anal. Appl. \textbf{326} (2007), no. 1, 75-87.

  \bibitem{kt}Kaplick\'y, P., Tich\'y, J.; 
\emph{Boundary regularity of flows under perfect slip boundary conditions,}
 Cent. Eur. J. Math., \textbf{11(7)} (2013), 1243--1263.

  \bibitem{mnrr} M\'alek, J., Ne\v{c}as, J., Rokyta, M. and Rų\v{z}i\`eka, M.; 
\emph{Weak and Measure-valued Solutions to Evolutionary PDEs,} 
Vol. 13 of Applied Mathematics and Mathematical Computation. Chapman 
\& Hall, London, 1996.

  \bibitem{mnr} M\'alek, J., Ne\v{c}as, J. and Rų\v{z}i\`eka, M.; 
\emph{On a weak solution to a class of non-Newtonian incompressible fluids 
in bounded three-dimensional domains: the case $p \ge 2$,} 
Adv. Differential Equations \textbf{6} (2001), 257--302.

  \bibitem{ns} Ne\v{c}as, J., \v{S}ver\'ak, V.;
\emph{On regularity of solutions of nonlinear parabolic systems} Ann. Scuola Norm. 
Sup. Pisa Cl. Sci. \textbf{18} (1991)(4), 1--11.

 \bibitem{shimada} Shimada, R.; 
\emph{On the $L_p-L_q$ maximal regularity for Stokes equations with Robin 
boundary condition in a bounded domain}, Math. Methods Appl. Sci.  
\textbf{30}  (2007),  no. 3, 257-289.

  \bibitem{ss} Shibata, Y., Shimada, R.; 
\emph{On a generalized resolvent estimate for the Stokes system with Robin 
boundary condition}, J. Math. Soc. Japan \textbf{59} (2007) (2), 469--519.

  \bibitem{steiger} Steiger, O.; 
\emph{Navier-Stokes equations with first order boundary conditions}, 
J. Math. Fluid Mech. \textbf{8}  (2006),  no. 4, 456-481.

\bibitem{wr} Wachsmuth, D., Roub\'i\`eek, T.; 
\emph{Optimal control of planar flow of incompressible non-Newtonian fluids,} 
Z. Anal. Anwend. \textbf{29} (2010), no. 3, 351-376.

  \bibitem{ziemer} Ziemer, W.P.; 
\emph{Weakly differentiable functions,} Sobolev spaces and functions of 
bounded variation, vol. 120 of Graduate Texts in Mathematics. 
New York: Springer 1989.

\end{thebibliography}

\end{document}

