\mathbf{e}gin{document}

\numberwithin{equation}{section}

\title[On a viscoelastic Stokes problem for salt rocks]{On the existence, uniqueness and regularity of solutions of a viscoelastic Stokes problem modelling salt rocks}

\author{R. A. Cipolatti}
\author{I.-S. Liu}
\author{L. A. Palermo}
\author{M. A. Rincon}
\author{R. M. S. Rosa}

\address[Rolci A. Cipolatti]{Instituto de Matem\'atica, Universidade Federal do Rio de Janeiro, Brazil}
\address[I-Shih Liu]{Instituto de Matem\'atica, Universidade Federal do Rio de Janeiro, Brazil}
\address[Luiz A. Palermo]{CENPES/Petrobr\'as, Rio de Janeiro, Brasil}
\address[Mauro A. Rincon]{Instituto de Matem\'atica, Universidade Federal do Rio de Janeiro, Brazil}
\address[Ricardo M. S. Rosa]{Instituto de Matem\'atica, Universidade Federal do Rio de Janeiro, Brazil}

\email[R. A. Cipolatti]{cipolatti@im.ufrj.br}
\email[I.-S. Liu]{liu@im.ufrj.br}
\email[L. A. Palermo]{luizpalermo@petrobras.com.br}
\email[M. A. Rincon]{rincon@dcc.ufrj.br}
\email[R. M. S. Rosa]{rrosa@im.ufrj.br}

\date{January 25, 2017}

\thanks{All the authors acknowledge the financial support of CENPES/PETROBR\'AS. The last author was also partly supported by CNPq, Bras\'{\i}lia, Brazil.}

\subjclass[2000]{35Q30, 76D06, 35B40, 37A60}
\keywords{generalized Stokes problem, elliptic regularization, viscoelastic fluids, salt modelling}

\mathbf{e}gin{abstract}
A Stokes-type problem for a viscoelastic model of salt rocks is considered, and existence, uniqueness and regularity are investigated in the scale of $L^2$-based Sobolev spaces. The system is transformed into a generalized Stokes problem, and the proper conditions on the parameters of the model that guarantee that the system is uniformly elliptic are given. Under those conditions, existence, uniqueness and low-order regularity are obtained under classical regularity conditions on the data, while higher-order regularity is proved under less stringent conditions than classical ones. Explicit estimates for the solution in terms of the data are given accordingly.
\end{abstract}

\maketitle

\section{Introduction}
We are interested in the viscoelastic Stokes-type problem 
\mathbf{e}gin{equation}
  \label{viscoelasticStokesorg}
   \mathbf{e}gin{cases}
      - \mu_1 \Delta \mathbf{v} - \mu_2 \mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v})\mathsf{B} + \mathsf{B}\mathsf{D}(\mathbf{v})) \\ 
      \qquad \qquad \qquad -\mu_ 3\mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v})\mathsf{B}^{-1} + \mathsf{B}^{-1} \mathsf{D}(\mathbf{v})) + \mathbf{n}abla p = \mathbf{b}f, & \text{in } \Omega,\\
      \mathbf{n}abla \cdot \mathbf{v} = 0, & \text{in } \Omega,\\
      \mathbf{v} = 0, & \text{on } \partial \Omega,
    \end{cases}
\end{equation}
on a three-dimensional bounded, smooth domain $\Omega\subset {\mathbb{R}}^3$, where $\mathbf{x}=(x_1,x_2,x_3)\in {\mathbb{R}}^3$ is the spatial variable and the unknowns are the scalar field $p=p(\mathbf{x})$ and the three-dimensional vector field $\mathbf{v}=\mathbf{v}(\mathbf{x})$, with coordinates $\mathbf{v}=(v_1,v_2,v_3)$. In this system of equations, $\mathsf{D}=\mathsf{D}(\mathbf{v})$ is the symmetric gradient tensor given by
\mathbf{e}gin{equation}
  \label{straintensor}
  \mathsf{D}(\mathbf{v}) = \frac{1}{2}(\mathbf{n}abla \mathbf{v} + (\mathbf{n}abla \mathbf{v})^\text{\rm tr}),
\end{equation}
where the superscript ``${}\text{\rm tr}$'' denotes the transpose of a matrix; $\mathsf{B}=\mathsf{B}(\mathbf{x})$ is a given symmetric, positively defined tensor with coordinates $\mathsf{B}=(b_{ij})_{ij}$; and the parameters $\mu_1, \mu_2, \mu_3$ are scalar fields.

Our aim is to obtain conditions on the parameters $\mu_1, \mu_2$, and $\mu_3$, and on the tensor $\mathsf{B}$ that guarantee that the problem \eqref{viscoelasticStokesorg} is uniformly elliptic and that we have existence, uniqueness and regularity properties of the solutions. We also aim to give explicit estimates for the regularity of the solutions depending on the regularity of $\mathsf{B}$, $\mathbf{b}f$, and the parameters $\mu_1, \mu_2, \mu_3$. 

Problem \eqref{viscoelasticStokesorg} does not seem to have been treated directly before. We rewrite \eqref{viscoelasticStokesorg} into a form that fits previously considered general system and verify the necessary conditions that guarantee the existence, uniqueness and the desired regularity up to $H^2$ in $\mathbf{v}$ and $H^1$ in $p$. For higher-order regularity, the known general results require too much regularity from $\mathsf{B}$ and we present a proof that yields the desired regularity upon less stringent conditions. Then we work out the estimates in detail. These are particularly important for the application of this result to the time-dependent problem mentioned below and addressed in a follow-up work.

The parameters $\mu_1, \mu_2$, and $\mu_3$ need not be constant, they are assumed to be general scalar fields. In applications, this dependence is typically on invariants of the model, such as the trace  of $\mathsf{B}$, or of powers of $\mathsf{B}$, and they may also depend on the temperature distribution of the medium. The actual dependence is important when analysing the corresponding evolutionary problem coupled with  equations for $\mathsf{B}$ and for the temperature. Here, however, we consider $\mu_1, \mu_2$, and $\mu_3$ as given general scalar fields. 

Problem \eqref{viscoelasticStokesorg} arises as the dissipative part for the linear momentum equation of a homogeneous and incompressible model of salt rocks, for which the constitutive law for the stress tensor $\mathsf{\sigma}$ is of the form
 \mathbf{e}gin{equation}
  \label{sigmaconstitutive}
   \mathsf{\sigma} = -p\mathsf{I} + 2\mu_1 \mathsf{D} + s_1 \mathsf{B} + s_{-1} \mathsf{B}^{-1} + \mu_2 (\mathsf{D}\mathsf{B} + \mathsf{B}\mathsf{D}) + \mu_3(\mathsf{D}\mathsf{B}^{-1} + \mathsf{B}^{-1}\mathsf{D}).
\end{equation}
In \eqref{sigmaconstitutive}, the tensor $\mathsf{I}$ is the identity in three-dimension; $p$ is the kinematic pressure; $\mathbf{v}$ is the velocity field; $\mathbf{n}abla \mathbf{v}$ is the velocity gradient; $\mathsf{D}$ is the strain-rate tensor defined in \eqref{straintensor}, which is the symmetric part of the velocity gradient;  $\mathsf{B}=\mathsf{F}\mathsf{F}^\text{\rm tr}$ is the left Cauchy-Green tensor; $\mathsf{F}$ is the deformation gradient tensor; $s_1$ and $s_{-1}$ are Mooney-Rivlin elastic parameters; the $\mu_1$ term is a Newtonian viscous stress; and $\mu_2$ and $\mu_3$ are Mooney-Rivlin-type parameters for the combined viscoelastic terms in the constitutive law (see \cite{truesdellnoll2004, haupt2002, Liu2002, temammiranville2005} for the fundamentals of continuum mechanics). 

The material model given by the stress tensor form \eqref{sigmaconstitutive} was introduced in 
\cite{LiuCipolattiRinconPalermo2014} and was referred to as \emph{Mooney-Rivlin type isotropic viscoelastic solid}. Notice that the purely elastic part of the constitutive equation is the same as that for Mooney-Rivlin isotropic elastic solids, while the viscous part is linear in the strain-rate tensor $\mathsf{D}$ and depends on the Cauchy-Green stress tensor $\mathsf{B}$ in a form similar to the purely elastic part.  This constitutive equation was successfully used to model the salt migration and the associated formation of salt diapirs (or salt domes) in a two-layer sediment-salt body problem, in which a denser elastic sediment material lies on top of a lighter layer of viscoelastic salt rock. The influence of the temperature in the formation of these salt domes was then investigated in \cite{teixeiraliurinconcipolattipalermo2014}. For other models and more information on this important geophysical problem, see the literature referenced in these two articles.

The associated evolutionary system of equations are written in the form
\mathbf{e}gin{equation}
  \label{LiueulerB}
    \mathbf{e}gin{cases}
      \mathsf{B}_t + (\mathbf{v}\cdot \mathbf{n}abla) \mathsf{B} = (\mathbf{n}abla \mathbf{v}) \mathsf{B} + ((\mathbf{n}abla \mathbf{v})\mathsf{B})^\text{\rm tr}, \\
      \rho(\mathbf{v}_t + (\mathbf{v}\cdot \mathbf{n}abla) \mathbf{v}) + \mathbf{n}abla p =  \nabla \cdot \mathsf{T}(\mathsf{B},\mathsf{D}(\mathbf{v})) - \rho g \mathbf{e}_3, \\
      \theta_t + (\mathbf{v}\cdot\mathbf{n}abla) \theta = \mathbf{n}abla \cdot (\mathsf{Q}(\mathsf{B}) \mathbf{n}abla \theta),  \\
      \nabla \cdot \mathbf{v} = 0, 
    \end{cases}
\end{equation}
on the domain $\Omega$ and time interval $[0,T]$, with $T>0$, along with appropriate boundary conditions on $\partial \Omega$ and initial conditions at $t=0$. In the system \eqref{LiueulerB}, $g$ is the acceleration of gravity, $\mathbf{e}_3=(0,0,1)$ is the vertical direction, the viscoelastic stress tensor $\mathsf{T}(\mathsf{B},\mathsf{D})$ is given according to $\mathsf{T} = \mathsf{\sigma} + p\mathsf{I}$, where $\mathsf{\sigma}$ is as is \eqref{sigmaconstitutive}, the thermal conductivity tensor is given by $\mathsf{Q}(\mathsf{B}) = \kappa_1 \mathsf{I} + \kappa_2 \mathsf{B} + \kappa_3 \mathsf{B}^{-1}$, and the viscoelastic parameters $s_1$, $s_{-1}$, $\mu_1$, $\mu_2$ and $\mu_3$ and the thermal conductivity coefficients $\kappa_1, \kappa_2$, and $\kappa_3$ may depend on invariants of the left Cauchy-Green tensor $\mathsf{B}$ and on the temperature field $\theta$.

Notice that in the formulation \eqref{LiueulerB} the left Cauchy-Green stress tensor $\mathsf{B}$ becomes an unknown. This is the trick to write an Eulerian formulation for the model, instead of the classical Lagrangian formulation. 

The formulation \eqref{LiueulerB} is inspired by a similar formulation given in \cite{linliuzhang2005, linzhang2008} for a particular type of Oldroyd-B viscoelastic model (see also \cite{guillopesaut1990, lionsmasmoudi2000, cheminmasmoudi2001, liuwalkington2001, sideristhomases2005} for Eulerian formulations of more general types of Oldroyd-B models). In the model considered in \cite{linliuzhang2005, linzhang2008} the stress tensor is of the form $\sigma = -p\mathsf{I} + 2\mu\mathsf{D} + s\mathsf{F}\mathsf{F}^\text{\rm tr}$, so that the linear part is simply the classical Stokes problem, and the unknown is actually the deformation tensor $\mathsf{F}$. Our formulation, besides including a more complicated constitutive law, uses instead $\mathsf{B}$ as an unknown, which seems to have some computational advantages. The Lagrangian position $\mathbf{x}(t,\mathbf{X})$ at time $t$ of a particle which initially was at $\mathbf{x}(0,\mathbf{X}) =\mathbf{X}=(X_1,X_2,X_3)$ can be recovered from the vector field $\mathbf{v}$, and then it can be showed that $\mathsf{B}$ is precisely $\mathsf{F}\mathsf{F}^\text{\rm tr}$, for $\mathsf{F}=\mathsf{F}(t,\mathbf{x})$ given according to $\mathsf{F}(t,\mathbf{x}(t,\mathbf{X})) = D_\mathbf{X} \mathbf{x}(t,\mathbf{X})=(\partial x_i/\partial_{X_j}(t,\mathbf{X}))_{ij}$ (see \cite[Lemma 1.1]{linliuzhang2005} for the corresponding result in the context of the model in \cite{linliuzhang2005}, in which the unkown is $\mathsf{F}$ itself). 

The details for the initial-value problem \eqref{LiueulerB} will be given in a forthcoming paper. The current work is an important step towards that, since proving existence, uniqueness, and regularity properties of the viscoelastic Stokes problem \eqref{viscoelasticStokesorg} is critical for proving local well-posedness of the evolutionary problem \eqref{LiueulerB}. In particular, we will exploit the explicit estimates obtained here to derived suitable a~priori estimates for the solutions of the evolutionary system.

The model of salt dynamics that we are interested in corresponds to particular combinations of the parameters $\mu_1, \mu_2$, and $\mu_3$, which may be different for different types of salts (see \cite{LiuCipolattiRinconPalermo2014} and \cite{teixeiraliurinconcipolattipalermo2014}). Nevertheless, we analyse the mathematical properties of the viscoelastic Stokes problem \eqref{viscoelasticStokesorg} for more general choices of the parameters, for possibly other types of materials.

In the Eulerian formulation \eqref{LiueulerB}, since the left Cauchy-Green tensor $\mathsf{B}$ is taken to be an unknown variable of the system, it is also of fundamental importance for us to derive the explicit dependence of the Sobolev estimates for the solution $\mathbf{v}$ of problem \eqref{viscoelasticStokesorg} in terms not only on $\mathbf{b}f$ but also on $\mathsf{B}$. This is part of the difficulty of the problem.

Regarding the notation, the operator $\mathbf{n}abla = (\partial_{x_i})_i$ acts to yield the vector field $\mathbf{n}abla p = (\partial_{x_i}p)_i$, the tensor $(\mathbf{n}abla \mathbf{v})_{ij}=(\partial_{x_j}v_i)_{ij}$, the scalar field $\mathbf{n}abla \cdot \mathbf{v} = \sum_i \partial_{x_i} v_i$, and the vector field $\mathbf{n}abla\cdot \mathsf{A} = (\sum_j\partial_{x_j}a_{ij})_i$, when $p$ is a scalar field, $\mathbf{v}=(v_i)_i$ is a vector field, and $\mathsf{A} = (a_{ij})_{ij}$ is a tensor.  The Laplace operator can be written as $\Delta = \mathbf{n}abla \cdot \mathbf{n}abla$. For vectors $\mathbf{u}, \mathbf{v}$ and matrices $\mathsf{A}, \mathsf{B}$, we write $\mathbf{u}\cdot\mathbf{v} = \sum_i u_iv_i$, $\mathbf{u}\cdot \mathsf{A} = (\sum_j u_j a_{ij})_{i}$, and $\mathsf{A} : \mathsf{B} = \sum_{ij} a_{ij}b_{ij}$, as well as $\mathsf{A}\mathsf{B} = (a_{ik}b_{kj})_{ij}$ and $\mathbf{u}\otimes \mathbf{v} = (u_iv_j)_{ij}$.

Now, concerning the system \eqref{viscoelasticStokesorg}, we address it in the following way. First, notice that, due to the divergence-free condition, we have
\[ \mathbf{n}abla \cdot (\mathbf{n}abla \mathbf{v})^\text{\rm tr} = \left(\sum_i \partial_{x_i} \partial_{x_j} v_i \right)_j = \mathbf{n}abla (\mathbf{n}abla \cdot \mathbf{v}) = 0,
\]
so that
\[ \Delta \mathbf{v} = \mathbf{n}abla \cdot \mathbf{n}abla \mathbf{v} = \mathbf{n}abla \cdot ((\mathbf{n}abla \mathbf{v}) + (\mathbf{n}abla \mathbf{v})^\text{\rm tr})  = 2\mathbf{n}abla \cdot \mathsf{D}(\mathbf{v}) = \mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v}) \mathsf{I} + \mathsf{I} \mathsf{D}(\mathbf{v})).
\]
Hence, we write
\mathbf{e}gin{multline*}  \mu_1 \Delta \mathbf{v} + \mu_2 \mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v})\mathsf{B} + \mathsf{B}\mathsf{D}(\mathbf{v})) +\mu_ 3\mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v})\mathsf{B}^{-1} + \mathsf{B}^{-1} \mathsf{D}(\mathbf{v})) \\
  = \mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v})(\mu_1\mathsf{I} + \mu_2 \mathsf{B} + \mu_3 \mathsf{B}^{-1}) + (\mu_1\mathsf{I} + \mu_2 \mathsf{B} + \mu_3 \mathsf{B}^{-1})\mathsf{D}(\mathbf{v})).
\end{multline*}
In order to simplify the notation and also the analysis of the problem, we introduce the operator
\mathbf{e}gin{equation}
  \label{defAcal}
  {\mathcal A}(\mathsf{B}) = \mu_1 \mathsf{I} + \mu_2 \mathsf{B} + \mu_3 \mathsf{B}^{-1},
\end{equation} 
and we write
\[  \mu_1 \Delta \mathbf{v} + \mu_2 \mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v})\mathsf{B} + \mathsf{B}\mathsf{D}(\mathbf{v})) +\mu_ 3\mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v})\mathsf{B}^{-1} + \mathsf{B}^{-1} \mathsf{D}(\mathbf{v}))
  = \mathsf{D}(\mathbf{v}){\mathcal A}(\mathsf{B}) + {\mathcal A}(\mathsf{B})\mathsf{D}(\mathbf{v}).
\]

This allows us to rewrite the viscoelastic Stokes problem \eqref{viscoelasticStokesorg} in the form
\mathbf{e}gin{equation}
  \label{viscoelasticStokes}
   \mathbf{e}gin{cases}
      - \mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v}){\mathcal A}(\mathsf{B})+ {\mathcal A}(\mathsf{B})\mathsf{D}(\mathbf{v})) + \mathbf{n}abla p = \mathbf{b}f, & \text{in } \Omega,\\
      \mathbf{n}abla \cdot \mathbf{v} = 0, & \text{in } \Omega,\\
      \mathbf{v} = 0, & \text{on } \partial \Omega.
    \end{cases}
\end{equation}
We then separate the analysis of problem \eqref{viscoelasticStokes} into two parts. First we consider the generalized Stokes problem
\mathbf{e}gin{equation}
  \label{generalizedStokes}
   \mathbf{e}gin{cases}
      - \mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v})) + \mathbf{n}abla p = \mathbf{b}f, & \text{in } \Omega,\\
      \mathbf{n}abla \cdot \mathbf{v} = 0, & \text{in } \Omega,\\
      \mathbf{v} = 0, & \text{on } \partial \Omega.
    \end{cases}
\end{equation}
for a given tensor $\mathsf{A}$, and investigate the existence, uniqueness and regularity of solutions of problem \eqref{generalizedStokes} in the case that $\mathsf{A}$ is an essentially bounded and uniformly positive definite symmetric tensor. 

The existence and uniqueness of solutions of the weak form of \eqref{generalizedStokes} together with the regularity $\mathbf{v}\in H_0^1(\Omega)$, $p\in L^2(\Omega)$ follows from a classical application of Lax-Milgram Theorem, once the ellipticity is guaranteed (see Theorem \ref{weaksolutionweakgeneralizedStokeswithpressure}). The regularity $\mathbf{v}\in H^2(\Omega)$, $p\in H^1(\Omega)$ follows by applying a result of Mach\'a \cite[Theorem 4.2 with $k=0$]{macha2011} for a slightly more general system, once all the conditions are verified, and along with that we obtain more explicit estimates in terms of all the data of the problem (see Theorem \ref{h2regularity}). For higher-order regularity, the case $k\geq 1$ in \cite[Theorem 4.2]{macha2011} is not optimal, and neither is any other result that we are aware of (the main reason is that the conditions on $\mathsf{A}$ are imposed independently of the space dimension and the Sobolev embeddings for particular dimensions are not explored; see Remark \ref{dontneedsomuchregularityonA}). We, therefore, work out the details of the proof for higher-order regularity and obtain the desired result under less stringent regularity assumptions on $\mathsf{A}$, along with the associated explicit estimates (see Theorem \ref{hkplustworegularity}).

Once the generalized Stokes problem \eqref{generalizedStokes} has been understood, we study the transformation \eqref{defAcal}, i.e. $\mathsf{B} \longmapsto {\mathcal A}(\mathsf{B})$, and investigate the conditions on $\mathsf{B}$ and on $\mu_1, \mu_2, \mu_3$ that yield the proper conditions on ${\mathcal A}(\mathsf{B})$ to ensure the regularity and the uniform ellipticity of \eqref{generalizedStokes} with $\mathsf{A} = {\mathcal A}(\mathsf{B})$. More precisely, we give conditions on the parameters $\mu_1$, $\mu_2$, $\mu_3$ and $\mathsf{B}$ that guarantee that ${\mathcal A}(\mathsf{B})$ is uniformly positive definite (Propositions \ref{propessinfglambdacondition} and \ref{propmuscenarios}) and sufficiently regular (Propositions \ref{proplinftyboundonacalb}, \ref{propwoneinftyboundonacalb}, \ref{propw23boundonacalb}, and \ref{prophmboundonacalb}). This yields a condition for \eqref{viscoelasticStokesorg} to be uniformly elliptic and for the existence, uniqueness and regularity results to hold for the viscoelastic Stokes problem as well (Theorems \ref{weaksolutionweakviscoelasticdStokeswithpressure}, \ref{h2regularityforviscoelasticStokes}, and \ref{hkplustworegularityforviscoelasticStokes}).

\section{Function spaces}

In this work, we consider a bounded domain $\Omega\subset {\mathbb{R}}^3$, i.e. an open and connected bounded set in ${\mathbb{R}}^3$. We denote by $L^p(\Omega)$ the classical Lebesgue spaces with $1\leq p \leq \infty$, with norm $\|\cdot \|_{L^p(\Omega)}$, and by $W^{k,p}(\Omega)$ the classical Sobolev spaces with $k\in {\mathbb{N}}$ and $1\leq p \leq \infty$, of functions with derivatives up to order $k$ in $L^p(\Omega)$, with norm $\|\cdot\|_{W^{k,p}(\Omega)}$. In the case $p=2$, we denote $H^k(\Omega) = W^{k,2}(\Omega)$, with norm $\|\cdot\|_{H^{k}(\Omega)}$. 

We consider the three-dimensional vector-valued versions $L^p(\Omega)^3$, $W^{k,p}(\Omega)^3$, and $H^k(\Omega)^3$ of those spaces, although for notational simplicity we keep the notation for the norms and the inner products the same. For instance, the $L^p$ norm, $1\leq p < \infty$, of a vector field is such that
\[ \|\mathbf{v}\|_{L^p(\Omega)} = \left(\sum_{i=1}^3 \|v_i\|_{L^p(\Omega)}^p\right)^{1/p},
\]
while the $W^{k,p}$ norm, $k\in {\mathbb{N}}$, $1\leq p \leq \infty$ is such that
\[ \|\mathbf{v}\|_{W^{k,p}(\Omega)} = \left( \sum_{i=1}^3 \|v_i\|_{W^{k,p}(\Omega)}^p \right)^{1/p}.
\]
The vector valued version of the space of test functions and distributions and their duality product is also considered. Similarly, we consider the tensor-valued version of the spaces $L^p(\Omega)$, $W^{k,p}(\Omega)$, and $H^k(\Omega)$, but keep the notation for the norms and inner products the same as those for scalars and vector fields.

We identify the dual of $L^2(\Omega)$ with itself and denote the dual space of $H_0^1(\Omega)$ by $H^{-1}(\Omega)$. For the sake of notational simplicity, we write the duality product as an integral even when the functional does not belong to $L^2(\Omega)$, i.e.
\[ \dual{f,\varphi}_{H^{-1}(\Omega),H_0^1(\Omega)} = \int_\Omega f \varphi \;{\text{\rm d}}\mathbf{x}.
\]

For higher-order derivatives, we consider a multi-index $\mathbf{g}amma = (\gamma_1, \gamma_2, \gamma_3)\in {\mathbb{N}}_0^3$, where each $\gamma_i\in {\mathbb{N}}_0$ is a non-negative integer. The order of the multi-index is the number $|\mathbf{g}amma| = \gamma_1 + \gamma_2 + \gamma_3$. The associated derivative $D^\mathbf{g}amma$ of a scalar function $\varphi$ is 
\[ D^\mathbf{g}amma \varphi = \partial_{x_1}^{\gamma_1}\partial_{x_2}^{\gamma_2}\partial_{x_3}^{\gamma_3}\varphi = \frac{\partial^{|\mathbf{g}amma|}\varphi}{\partial_{x_1}^{\gamma_1}\partial_{x_2}^{\gamma_2}\partial_{x_3}^{\gamma_3}}. 
\]
Similarly for the $D^\mathbf{g}amma$ derivative of vector fields and tensors. Notice that we can write $\mathbf{g}amma = \gamma_1\mathbf{e}_1 + \gamma_2\mathbf{e}_2 + \gamma_3\mathbf{e}_3$, where $\mathbf{e}_1=(1,0,0)$, $\mathbf{e}_2=(0,1,0)$, and $\mathbf{e}_3=(0,0,1)$ form the canonical basis of ${\mathbb{R}}^3$ and can also be regarded as the multi-indices of the first-order derivatives. Combining all the derivatives of a certain order $k$ yields the tensor 
\[  D^k\varphi = (D^\mathbf{g}amma \varphi)_{|\mathbf{g}amma| = k}, 
\]
where the index $\mathbf{g}amma$ runs through all the ${\mathbb{N}}_0^3$ indices with order $|\mathbf{g}amma|=k$.

For dimensional reasons, we define the equivalent norms
\mathbf{e}gin{equation}
  \label{dimensionalhigherordernorm}
  \| \varphi\|_{W^{k,p}_\lambda(\Omega)} = \sum_{j=0}^k \lambda^{(k-j)/2}\|D^j\varphi\|_{L^p(\Omega)},
\end{equation}
where $\lambda$ is a parameter with the physical dimensional of $(\text{length})^{-2}$ (such as the constant in the Poincar\'e inequality mentioned in \eqref{poincareineq} below). Notice that $\|\varphi\|_{W^{k,p}_\lambda(\Omega)}$ has the same physical dimension (unit) as $\|D^k\varphi\|_{L^p(\Omega)}$. Similarly for the vector-valued and tensor-valued versions of these spaces. In the case $p=2$, we write $W^{k,2}_\lambda(\Omega)=H^k_\lambda(\Omega)$. For negative powers, $H^{-k}(\Omega)$, for $k\in {\mathbb{N}}$, is the dual of $H_0^k(\Omega)$, and $H^{-k}_{\lambda}(\Omega)$ is the dual norm associated with the norm $H^k_{\lambda}(\Omega)$ in the subspace $H_0^k(\Omega)$, i.e. 
\[ \| \xi \|_{H^{-k}_{\lambda}(\Omega)} = \sup_{\varphi\in H_0^k(\Omega), \;\|\varphi\|_{H^k_\lambda(\Omega)=1}}\dual{\xi,\varphi}, \quad \forall \xi\in H^{-k}_{\lambda}(\Omega),
\]
where $\dual{\cdot,\cdot}$ is the corresponding duality product.

A useful inequality is the following interpolation inequality, balancing different orders of derivatives:
\mathbf{e}gin{equation}
  \label{interpolationineq}
  \|\varphi\|_{H^m_\lambda(\Omega)} \leq c(\Omega)\|\varphi\|_{H^j_\lambda(\Omega)}^\theta\|\varphi\|_{H^k_\lambda(\Omega)}^{(1-\theta)},
\end{equation}
for any $\varphi\in H^k_\lambda(\Omega)$ and for any $j\leq m\leq k$, where $c(\Omega)$ is a non-dimensional constant that only depends on the shape of $\Omega$, and $\theta\in ]0,1[$ is given by
\mathbf{e}gin{equation}
  \label{interpolationtheta}
  m = \theta j + (1-\theta)k.
\end{equation}
The same inequality holds with $\varphi$ replaced by a vector field or a tensor.

Concerning the constant that appears in \eqref{interpolationineq} and in many other inequalities involving functional norms given in this work, we call a ``universal constant'' any term that is independent of any data of the problem, and by a ``shape constant'' any term that may depend only on the shape of $\Omega$, i.e. it is independent of translations, rotations and dilations of the spatial domain. The universal constants will be denoted by $c$, while shape constants will be denoted by $c(\Omega)$. From a physical point of view, they are non-dimensional constants. All the estimates in this work have been properly dimensionalized, in the sense that the constants are non-dimensional, so that all the dimensional quantities appear explicitly and with the proper dimensionalization (e.g. the exact power).

As usual, the pressure term in \eqref{viscoelasticStokes} or in \eqref{generalizedStokes} satisfies an equation of the form $\mathbf{n}abla p = \mathbf{g}$ and is, therefore, obtained up to a constant. In this regard, it is useful to consider the quotient space $L^2(\Omega)/{\mathbb{R}}$, which is the space $L^2(\Omega)$ quotient equivalence up to translation by a constant. More precisely, two functions $\varphi$ and $\psi$ in $L^2(\Omega)$ are equivalent up to translations by a constant if $\varphi-\psi$ is constant in $\Omega$. In particular, $\varphi$ is in the same equivalent class as its zero mean representative $\varphi - \int_\Omega \varphi$. Such a quotient space is also a Banach space. Given a function $\varphi$ in $L^2(\Omega)$, this function belongs to an equivalence class in $L^2(\Omega)/{\mathbb{R}}$ which with a certain abuse of notation we still denote by $\varphi$, for simplicity, and is such that its norm in the equivalence class is given by
\[ \|\varphi\|_{L^2(\Omega)/{\mathbb{R}}} = \inf_{c\in {\mathbb{R}}} \|\varphi + c\|_{L^2(\Omega)} = \left\|\varphi - \int_\Omega \varphi\right\|_{L^2(\Omega)}.
\]

We denote by ${\mathcal C}^0(\mathbf{a}r\Omega)$ the space of continuous functions defined on $\mathbf{a}r\Omega$ and by ${\mathcal C}^k(\mathbf{a}r\Omega)$, the subspace of ${\mathcal C}^0(\mathbf{a}r\Omega)$ with derivatives up to order $k$ defined and continuous on $\Omega$ and with well defined extensions up to the boundary. We define the space ${\mathcal C}^{k,\alpha}(\mathbf{a}r\Omega)$, $k\in {\mathbb{Z}}$, $k\geq 0$, $0<\alpha\leq 1$, as the subspace of ${\mathcal C}^k(\mathbf{a}r\Omega)$ such that the derivatives of order $k$ are H\"older-continuous function with H\"older exponente $\alpha$, or a Lipschitz function, in the case $\alpha=1$.

We follow \cite[page 77]{adn2} (see also \cite{adamsfournier2003} for related and more general notions) and say that the domain $\Omega$ has a boundary of class ${\mathcal C}^k$, $k\in {\mathbb{Z}}$, $k\geq 0$, when its boundary $\partial \Omega$ can be covered by a finite number of open sets $U_m$, $m=1,\ldots, M$, $M\in {\mathbb{N}}$, such that each portion $\overline{U_m}\cap\overline{\Omega}$ of the domain is homeomorphic to a hemisphere, the portion $\overline{U_m}\cap\partial \Omega$ is the image of the flat face of the hemisphere, and the homeomorphism and its inverse are of class ${\mathcal C}^k$, with all derivatives up to order $k$ bounded by a constant independent of $m$. Similarly, we say that $\Omega$ has a boundary of class ${\mathcal C}^{k,\alpha}$, $k\in {\mathbb{Z}}$, $k\geq 0$, $0<\alpha\leq 1$, when $\Omega$ has a boundary of class ${\mathcal C}^k$ and the $k$-th derivatives of the associated homeomorphism and of its inverse are H\"older continuous with exponent $\alpha$ and with modulus of H\"older continuity bounded independently of $m$. In the case ${\mathcal C}^{0,1}$, we say that $\Omega$ has a Lipschitz boundary.

When $\Omega$ has a boundary of class ${\mathcal C}^1$, the space $W^{1,\infty}(\Omega)$ can be identified with the space of Lipschitz functions ${\mathcal C}^{0,1}(\mathbf{a}r\Omega)$, and the spaces $W^{1,p}(\Omega)$, $p>3$, can be identified with the space of H\"older functions ${\mathcal C}^{0,1 - 3/p}(\mathbf{a}r\Omega)$ (see \cite[Theorem 5.6.5, pg 283]{evans2010}).\footnote{If $\varphi\in W^{1,\infty}(\Omega)$, let $\varphi_\varepsilon$ be a mollification of $\varphi$, which converges to $\varphi$ in any $W^{1,p}_{\text{\rm loc}}(\Omega)$, for $1\leq p < \infty$. Thus, a subsequence of $\varphi_\varepsilon$ converges to $\varphi$ and is such that $\mathbf{n}abla\varphi_\varepsilon$ converges to $\mathbf{n}abla\varphi$, almost everywhere in $\Omega$. Since $|\delta_i^h\varphi_\varepsilon|\leq \sup_\Omega|\mathbf{n}abla\varphi_\varepsilon|$, we find at the limit that $\|\delta_i^h\varphi\|_{L^\infty(\Omega_0)} \leq \|\partial_{x_i}\varphi\|_{L^\infty(\Omega)}$, for every subdomain $\Omega_0$ compactly included in $\Omega$ and for every $h\neq 0$ sufficiently small such that $|h|\leq \textrm{dist}(\Omega_0,\partial\Omega)$.}

For the sake of the boundary conditions, we also consider the Hilbert space $H_0^1(\Omega)$, which is the completion, under the norm $H^1(\Omega)$, of the space of test functions ${\mathcal D}(\Omega)={\mathcal C}_0^\infty(\Omega)$ of compactly supported, infinitely-many-times continuously differentiable functions on $\Omega$.  The dual ${\mathcal D}(\Omega)'$ is the space of distributions, and the duality product between them is denoted $\dual{\cdot,\cdot}_{{\mathcal D}(\Omega)', {\mathcal D}(\Omega)'}$.

Since $\Omega$ is assumed to be bounded, we have the Poincar\'e inequality
\mathbf{e}gin{equation}
  \label{poincareineq}
  \int_\Omega \varphi^2(\mathbf{x}) \;{\text{\rm d}} \mathbf{x} \leq \frac{1}{\lambda_1} \int_\Omega |\mathbf{n}abla \varphi(\mathbf{x})|^2 \;{\text{\rm d}}\mathbf{x}, \quad \forall \varphi \in H_0^1(\Omega).
\end{equation}
The notation comes from the fact that $\lambda_1$ ends up being the first eigenvalue of the Laplace operator. The Poincar\'e inequality extends, with the same constant, to the vector-valued functions in $H_0^1(\Omega)^3$.

In view of the Poincar\'e inequality \eqref{poincareineq}, the norm in $H_0^1(\Omega)$ is taken to be 
\mathbf{e}gin{equation}
  \label{norminH01}
  \|\varphi\|_{H_0^1(\Omega)} = \|\mathbf{n}abla \varphi\|_{L^2(\Omega)}, 
\end{equation}
and similarly in the vector-valued case $H_0^1(\Omega)^3$. For higher-order derivatives, we have, by induction and using integration by parts and Cauchy-Schwarz inequality, that
\mathbf{e}gin{equation}
  \label{higherorderpoincare}
  \|D^j\varphi\|_{L^2(\Omega)}^2 \leq \frac{1}{\lambda_1^{k-j}} \|D^k\varphi\|_{L^2(\Omega)}^2, \quad \forall \varphi\in H_0^1(\Omega)\cap H^k(\Omega),
\end{equation}
for any integers $0\leq j \leq k$. In particular, the $H_{\lambda_1}^k$ norm in the subspace $H_0^1(\Omega)\cap H^k(\Omega)$ is equivalent to the $D^k$-norm, i.e.
\[ \|D^k\varphi\|_{L^2(\Omega)} \leq \|\varphi\|_{H_{\lambda_1}^k(\Omega)} \leq (k+1)\|D^k\varphi\|_{L^2(\Omega)}, \qquad \forall \varphi \in H_0^1(\Omega)\cap H^k(\Omega).
\]

In our case of a three-dimensional space, we have the embeddings $H^2(\Omega)\subset L^\infty(\Omega)$ (see \cite[Theorem 3.9]{agmon1965}) and $H^1(\Omega)\subset L^6(\Omega)$. In the particular case of vanishing boundary conditions, and using the Poincar\'e inequality \eqref{poincareineq}, we write, with the proper dimensionalization,
\mathbf{e}gin{equation}
  \label{agmoninequality} 
  \|D^j\varphi\|_{L^\infty} \leq \frac{c_j(\Omega)}{\lambda_1^{1/4}}\|D^{j+2}\varphi\|_{L^2(\Omega)}, \quad \forall \varphi\in H_0^1(\Omega) \cap H^{j+2}(\Omega)
\end{equation}
and
\mathbf{e}gin{equation}
  \label{l6inequality} 
  \|D^j\varphi\|_{L^6} \leq c_j(\Omega)\|D^{j+1}\varphi\|_{L^2(\Omega)}, \quad \forall \varphi\in H_0^1(\Omega)\cap H^{j+1}(\Omega),
\end{equation}
for any integer $j\geq 0$, where $c_j(\Omega)$ is a shape constant, which may be different for different $j$. Interpolating between \eqref{higherorderpoincare} and \eqref{l6inequality} we also find the inequality
\mathbf{e}gin{equation}
  \label{l3inequality} 
  \|D^j\varphi\|_{L^3} \leq \frac{c_j(\Omega)}{\lambda_1^{1/4}}\|D^{j+1}\varphi\|_{L^2(\Omega)}, \quad \forall \varphi\in H_0^1(\Omega) \cap H^{j+1}(\Omega),
\end{equation}

In the case that $\varphi$ does not vanish on the boundary, a version of \eqref{higherorderpoincare} holds with the full norm,
\mathbf{e}gin{equation}
  \label{higherorderscaledembedding}
  \|\varphi\|_{H^j_\lambda(\Omega)}^2 \leq \frac{c_{j,k}(\Omega)}{\lambda^{k-j}}\|\varphi\|_{H^k_\lambda(\Omega)}^2,
\end{equation}
for any $\varphi\in H^k_\lambda(\Omega)$ and any $j\leq k$, where $c_{j,k}(\Omega)$ is a shape constant.

For divergence-free vector fields, we also consider the space
\[ {\mathcal V}(\Omega) =  \{ \mathbf{v} \in {\mathcal C}_0^\infty(\Omega); \; \mathbf{n}abla \cdot \mathbf{v} = 0 \}.
\]
The completion of the space ${\mathcal V}(\Omega)$ under the norm of $L^2(\Omega)$ is denoted by $H(\Omega)$ and the completion of the space ${\mathcal V}(\Omega)$ under the norm of $H_0^1(\Omega)^3$ is denoted by $V(\Omega)$. The norms and inner products in $H(\Omega)$ and in $V(\Omega)$ are those inherited from $L^2(\Omega)^3$ and $H_0^1(\Omega)^3$, respectively. The Poincar\'e inequality also holds in $V(\Omega)$.

If we denote
\mathbf{e}gin{equation}
   \mathsf{D}(\mathbf{v}) = (d_{ij})_{ij} = \frac{1}{2}\left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i}\right)_{ij}, \quad d_{ij}= d_{ji},
\end{equation} 
then, 
\[ \sum_{ij} d_{ij}^2 = \frac{1}{4} \sum_{ij} \left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i}\right)^2 = \frac{1}{4} \sum_{ij} \left( \left(\frac{\partial v_i}{\partial x_j}\right)^2 + 2 \frac{\partial v_i}{\partial x_j}\frac{\partial v_j}{\partial x_i} + \left(\frac{\partial v_j}{\partial x_i}\right)^2 \right).
\]
Thus, 
\mathbf{e}gin{multline*}
   \|\mathsf{D}(\mathbf{v})\|_{L^2(\Omega)}^2 = \sum_{ij} \int_\Omega d_{ij}^2\;{\text{\rm d}} \mathbf{x} = \frac{1}{4} \sum_{ij} \int_\Omega\left( \left(\frac{\partial v_i}{\partial x_j}\right)^2 + 2 \frac{\partial v_i}{\partial x_j}\frac{\partial v_j}{\partial x_i} + \left(\frac{\partial v_j}{\partial x_i}\right)^2 \right) \;{\text{\rm d}}\mathbf{x} \\
   =  \frac{1}{2}\sum_{ij} \int_\Omega \left(\frac{\partial v_i}{\partial x_j}\right)^2 \;{\text{\rm d}}\mathbf{x} +  \frac{1}{2}\sum_{ij} \int_\Omega \frac{\partial v_i}{\partial x_j}\frac{\partial v_j}{\partial x_i} \;{\text{\rm d}} \mathbf{x} =  \frac{1}{2}\|\mathbf{n}abla\mathbf{v}\|_{L^2(\Omega)}^2 + \frac{1}{2}\sum_{ij} \int_\Omega \frac{\partial v_i}{\partial x_j}\frac{\partial v_j}{\partial x_i} \;{\text{\rm d}} \mathbf{x}.
  \end{multline*} 

If $\mathbf{v}\in {\mathcal D}(\Omega)^3$, then by two integration by parts and using that $\mathbf{v}$ vanishes on the boundary we find that
\[ \sum_{ij} \int_\Omega \frac{\partial v_i}{\partial x_j}\frac{\partial v_j}{\partial x_i} \;{\text{\rm d}} \mathbf{x} = \sum_{ij} \int_\Omega \frac{\partial v_i}{\partial x_i}\frac{\partial v_j}{\partial x_j} \;{\text{\rm d}} \mathbf{x} = \int_\Omega \left(\mathbf{n}abla \cdot \mathbf{v}\right)^2 \;{\text{\rm d}} \mathbf{x} = \|\mathbf{n}abla \cdot \mathbf{v}\|_{L^2(\Omega)}^2.
\]
Thus, by the density of ${\mathcal D}(\Omega)^3$ in $H_0^1(\Omega)^3$, we find that (c.f. Korn's inequality \cite[Theorem 6.3-3]{ciarlet1988})
\mathbf{e}gin{equation}
  \label{Dandgradvinhzero1}
  \|\mathsf{D}(\mathbf{v})\|_{L^2(\Omega)}^2 = \frac{1}{2}\|\mathbf{n}abla\mathbf{v}\|_{L^2(\Omega)}^2 + \frac{1}{2}\|\mathbf{n}abla \cdot \mathbf{v}\|_{L^2(\Omega)}^2, \qquad \forall \mathbf{v}\in H_0^1(\Omega)^3.
\end{equation}
In particular,
\mathbf{e}gin{equation}
  \label{Dandgradvcoerciveinhzero1}
  \frac{1}{2}\|\mathbf{n}abla\mathbf{v}\|_{L^2(\Omega)}^2 \leq \|\mathsf{D}(\mathbf{v})\|_{L^2(\Omega)}^2 \leq \|\mathbf{n}abla\mathbf{v}\|_{L^2(\Omega)}^2, \qquad \forall \mathbf{v}\in H_0^1(\Omega)^3.
\end{equation}

If, moreover, $\mathbf{v}\in V(\Omega)$, then $\mathbf{n}abla \cdot \mathbf{v} = 0$ and we are left with
\mathbf{e}gin{equation}
  \label{Dandgradv}
  \|\mathsf{D}(\mathbf{v})\|_{L^2(\Omega)}^2 = \frac{1}{2}\|\mathbf{n}abla\mathbf{v}\|_{L^2(\Omega)}^2, \qquad \forall \mathbf{v}\in V(\Omega).
\end{equation}

For the proof of higher-order regularity, we use the finite-difference operators
\mathbf{e}gin{equation}
  \label{finitedifference}
  \delta_i^h\mathbf{w}(\mathbf{x}) = \frac{\mathbf{w}(\mathbf{x} + h\mathbf{e}_i) - \mathbf{w}(\mathbf{x})}{h},
\end{equation}
for $i=1,2,3$ and $h>0$. A useful estimate related to these operators is the following:
  \mathbf{e}gin{align} 
    \label{finitedifferenceintestimate}
    \|\delta_i^h \varphi\|_{L^2(\Omega_0)} & \leq \|\partial_{x_i} \varphi\|_{L^2(\Omega)}, \quad \forall \varphi\in H^1(\Omega),
  \end{align}
for $h\neq 0$ sufficiently small such that $|h| < \textrm{dist}(\Omega_0, \partial\Omega)$, for any subdomain $\Omega_0$ compactly included in $\Omega$ (see \cite[Lemma 15.1]{friedman1969}). By extending the function $\varphi$ to zero outside $\Omega$, the extended function $\tilde\varphi$ satisfies
  \mathbf{e}gin{align} 
    \label{finitedifferenceintestimate2}
    \|\delta_i^h \tilde\varphi\|_{L^2(\Omega)} & \leq \|\partial_{x_i} \varphi\|_{L^2(\Omega)}, \quad \forall \varphi\in H^1(\Omega),
  \end{align}
for any $h\neq 0$. Of course, these estimates extend to vector fields and tensors. 

In the case $\varphi\in W^{1,\infty}(\Omega)$ (so that $\varphi$ is Lipschitz), then it is not difficult to see that $\delta_i^h\varphi$ is uniformly bounded in the interior, i.e.
\mathbf{e}gin{equation}
  \label{finitedifferenceintestimateinfty}
  \|\delta_i^h\varphi\|_{L^\infty(\Omega_0)} \leq \|\partial_{x_i}\varphi\|_{L^\infty(\Omega)},
\end{equation}
for any $h\neq 0$ sufficiently small such that $|h| < \textrm{dist}(\Omega_0, \partial\Omega)$, for any subdomain $\Omega_0$ compactly included in $\Omega$. For the estimates near the boundary, we localize and rectify the domain, so that we may assume that a portion of the boundary is locally flat and given by $\{x_3=0\}$. More precisely, assume there is a ``rectified domain'' $\tilde\Omega$ and a ball $B$ such that $B\cap \tilde\Omega = \{\mathbf{x}\in B; x_3<0\}$ and extend $\varphi$ to a function $\tilde\varphi$ that vanishes on $B\cap ({\mathbb{R}}^3\setminus \tilde\Omega)$. Then, we bound any ``tangential'' finite difference by the corresponding ``tangential derivative'', i.e.
\mathbf{e}gin{equation}
  \label{finitedifferencebdryestimateinfty}
  \|\delta_i^h\varphi\|_{L^\infty(B_0)} \leq \|\partial_{x_i}\varphi\|_{L^\infty(B\cap\Omega)},
\end{equation}
for every $i=1,2$ and every $h\neq 0$ such that $|h| < \textrm{dist}(B_0, \partial B)$, where $B_0$ is a ball concentric to $B$ and with strictly smaller radius.

\section{Weak formulations}

If $\mathbf{v}$ and $p$ are smooth solutions of \eqref{generalizedStokes} and $\mathbf{w}$ is another smooth vector field vanishing at the boundary, and assuming that $\mathsf{A}$ and $\mathbf{b}f$ are sufficiently smooth, then by multiplying scalar-wise the equation by $\mathbf{w}$ and integrating that by parts, using the boundary conditions, we arrive at
\[ \int_\Omega (\mathsf{D}(\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v})) : \mathbf{n}abla \mathbf{w} \;{\text{\rm d}} \mathbf{x} - \int_\Omega p \mathbf{n}abla \cdot \mathbf{w} \;{\text{\rm d}} \mathbf{x} = \int_\Omega \mathbf{b}f \cdot \mathbf{w} \;{\text{\rm d}} \mathbf{x}.
\]
The above equation still makes sense for less regular vector fields and tensors, which leads us to the following weak formulation for the generalized Stokes problem \eqref{generalizedStokes}.

\mathbf{e}gin{prob}[Weak formulation for the generalized Stokes problem with pressure]
  Given a tensor $\mathsf{A}\in L^\infty(\Omega)^{3\times 3}$ and a vector field $\mathbf{b}f\in H^{-1}(\Omega)^3$, find a vector field $\mathbf{v}\in V(\Omega)$ and a scalar field $p\in L^2(\Omega)$ such that
\mathbf{e}gin{equation}
    \label{weakgeneralizedStokeswithpressure}
    \int_\Omega (\mathsf{D}(\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v})) : \mathbf{n}abla \mathbf{w} \;{\text{\rm d}} \mathbf{x} - \int_\Omega p \mathbf{n}abla \cdot \mathbf{w} \;{\text{\rm d}} \mathbf{x} = \int_\Omega \mathbf{b}f \cdot \mathbf{w} \;{\text{\rm d}} \mathbf{x}, \qquad \forall \mathbf{w}\in H_0^1(\Omega)^3.
  \end{equation}  
\end{prob}

If we further restrict the test functions $\mathbf{w}$ to the space $V(\Omega)$ of divergence-free vector fields, then the pressure term disappears and we arrive at the following weak formulation.

\mathbf{e}gin{prob}[Weak formulation for the generalized Stokes problem]
  Given a tensor $\mathsf{A}\in L^\infty(\Omega)^{3\times 3}$ and a vector field $\mathbf{b}f\in H^{-1}(\Omega)^3$, find a vector field $\mathbf{v}\in V(\Omega)$ such that
  \mathbf{e}gin{equation}
    \label{weakgeneralizedStokes}
    \int_\Omega (\mathsf{D}(\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v})) : \mathbf{n}abla \mathbf{w} \;{\text{\rm d}} \mathbf{x} = \int_\Omega \mathbf{b}f \cdot \mathbf{w} \;{\text{\rm d}} \mathbf{x}, \qquad \forall \mathbf{w}\in V(\Omega).
  \end{equation}
\end{prob}

In the case $\mathsf{A}$ is replaced by ${\mathcal A}(\mathsf{B})$, we have the following weak formulation for the viscoelastic Stokes problem \eqref{viscoelasticStokes} with the pressure term (see Proposition \ref{proplinftyboundonacalb} to justify the hypotheses on $\mathsf{B}$).

\mathbf{e}gin{prob}[Weak formulation for the viscoelastic Stokes problem with pressure]
  Given a tensor $\mathsf{B}\in L^\infty(\Omega)^{3\times 3}$ with $(\det \mathsf{B})^{-1} \in L^\infty(\Omega)$, coefficients $\mu_1, \mu_2, \mu_3$ in $L^\infty(\Omega)$, and a vector field $\mathbf{b}f\in H^{-1}(\Omega)^3$, find a vector field $\mathbf{v}\in V(\Omega)$ and a scalar field $p\in L^2(\Omega)$ such that
  \mathbf{e}gin{equation}
    \label{weakviscoelasticStokeswithpressure}
    \int_\Omega (\mathsf{D}(\mathbf{v}){\mathcal A}(\mathsf{B}) + {\mathcal A}(\mathsf{B})\mathsf{D}(\mathbf{v})) : \mathbf{n}abla \mathbf{w} \;{\text{\rm d}} \mathbf{x} - \int_\Omega p \mathbf{n}abla \cdot \mathbf{w} \;{\text{\rm d}} \mathbf{x}  = \int_\Omega \mathbf{b}f \cdot \mathbf{w} \;{\text{\rm d}} \mathbf{x}, \qquad \forall \mathbf{w}\in H_0^1(\Omega)^3.
  \end{equation}
\end{prob}

Similarly, we also have the following weak formulation in $V(\Omega)$.

\mathbf{e}gin{prob}[Weak formulation for the viscoelastic Stokes problem]
  Given a tensor $\mathsf{B}\in L^\infty(\Omega)^{3\times 3}$ with $(\det \mathsf{B})^{-1} \in L^\infty(\Omega)$ and given a vector field $\mathbf{b}f\in H^{-1}(\Omega)^3$, find a vector field $\mathbf{v}\in V(\Omega)$ such that
  \mathbf{e}gin{equation}
    \label{weakviscoelasticStokes}
    \int_\Omega (\mathsf{D}(\mathbf{v}){\mathcal A}(\mathsf{B}) + {\mathcal A}(\mathsf{B})\mathsf{D}(\mathbf{v})) : \mathbf{n}abla \mathbf{w} \;{\text{\rm d}} \mathbf{x} = \int_\Omega \mathbf{b}f \cdot \mathbf{w} \;{\text{\rm d}} \mathbf{x}, \qquad \forall \mathbf{w}\in V(\Omega).
  \end{equation}
\end{prob}

Of course, if $\mathsf{B} = \mathsf{F}\mathsf{F}^\text{\rm tr}$, where $\mathsf{F}$ is the deformation tensor, then $\det \mathsf{B}=1$ thanks to the incompressibility condition. Nevertheless, the previous problems are stated for a general tensor $\mathsf{B}$.

In order to characterize the problems above as elliptic problems and to obtain suitable existence, uniqueness and continuous dependence on the data, an important condition is that of uniform ellipticity. In that regard, we make the following definition, suitable in our context.
\mathbf{e}gin{defs}
  \label{Stokesuniformlyelliptic}
  We say that, for a given tensor $\mathsf{A}$, the generalized Stokes problem \eqref{weakgeneralizedStokes} is \textbf{uniformly elliptic} in $V(\Omega)$, or, more generally, in $H_0^1(\Omega)^3$, if there exist $\delta' \geq \delta >0$ such that
  \mathbf{e}gin{multline}
    \label{probgenStokesunifellip}
    \delta \int_\Omega |\mathbf{n}abla \mathbf{v}|^2\;{\text{\rm d}}\mathbf{x} \leq \int_\Omega (\mathsf{D}(\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v})) : \mathbf{n}abla \mathbf{v} \;{\text{\rm d}} \mathbf{x} \\ \leq \delta' \int_\Omega |\mathbf{n}abla \mathbf{v}|^2\;{\text{\rm d}}\mathbf{x}, \quad \forall \mathbf{v} \in V(\Omega), \text{ or, respectively, } \forall \mathbf{v} \in H_0^1(\Omega)^3.
  \end{multline}  
\end{defs}

As we will see in the Section \ref{secgenstokes}, the key condition to show the uniform ellipticity of the generalized Stokes problem is that the tensor $\mathsf{A}$ be essentially bounded and uniformly positive definite (see Proposition \ref{probgenStokesisunifelliptic}). With that purpose in mind, we make the following definition.
\mathbf{e}gin{defs}
  \label{defuniformpositivedef}
  We say that a symmetric tensor $\mathsf{A}$ on $\Omega$ is \textbf{uniformly positive definite} if there exists 
\mathbf{e}gin{equation}
  \label{alpha}
  \alpha >0
\end{equation}
such that, for almost every $\mathbf{x}\in \Omega$,
\mathbf{e}gin{equation}
  \label{unifellip}
  \mathsf{A}(\mathbf{x})\mathbf{x}i \cdot \mathbf{x}i \geq \alpha |\mathbf{x}i|^2, \qquad \forall \mathbf{x}i \in {\mathbb{R}}^3.
\end{equation}
\end{defs}

Notice that, if \eqref{unifellip} holds, then, by taking $\mathbf{x}i$ to be one of the vectors of the canonical basis, we find that $\alpha$ is bounded by every element of the diagonal of the tensor $\mathsf{A}$. In particular,
\mathbf{e}gin{equation}
  \label{alphaboundedLinfty}
  \alpha \leq \|\mathsf{A}\|_{L^\infty}.
\end{equation}

Similarly, concerning the viscoelastic Stokes problem, we have the following definition of uniform ellipticity:
\mathbf{e}gin{defs}
  \label{viscoelasticuniformlyelliptic}
  We say that, for a given tensor $\mathsf{B}$, the viscoelastic Stokes problem \eqref{weakviscoelasticStokes} is \textbf{uniformly elliptic} in $V(\Omega)$, or, more generally, in $H_0^1(\Omega)^3$, if there exist $\delta'\geq \delta>0$ such that
  \mathbf{e}gin{multline}
    \label{probviscoelasticStokesunifellip}
    \delta \int_\Omega |\mathbf{n}abla \mathbf{v}|^2\;{\text{\rm d}}\mathbf{x} \leq \int_\Omega (\mathsf{D}(\mathbf{v}){\mathcal A}(\mathsf{B}) + {\mathcal A}(\mathsf{B})\mathsf{D}(\mathbf{v})) : \mathbf{n}abla \mathbf{v} \;{\text{\rm d}} \mathbf{x} \\
    \leq \delta' \int_\Omega |\mathbf{n}abla \mathbf{v}|^2\;{\text{\rm d}}\mathbf{x}, \quad \forall \mathbf{v} \in V(\Omega), \text{ or, respectively, } \forall \mathbf{v} \in H_0^1(\Omega)^3.
  \end{multline}  
\end{defs}

The condition on $\mathsf{B}$ that ensures the uniform ellipticity of the viscoelastic Stokes problem is investigated in Section \ref{conditionsellipticity}. For the moment, we only make the following definition, for the sake of notational simplification.

\mathbf{e}gin{defs}
  \label{defproper}
  A positive definite symmetric tensor $\mathsf{B}=\mathsf{B}(\mathbf{x})$ is called ${\mathcal A}$-\textbf{positive} when $\mathsf{A}={\mathcal A}(\mathsf{B})$ is a uniformly positive definite symmetric tensor in the sense of Definition \ref{defuniformpositivedef}.
\end{defs}

Hence, it will follow from Propositions \ref{probgenStokesisunifelliptic} and \ref{proplinftyboundonacalb} that every ${\mathcal A}$-positive tensor $\mathsf{B}$ in $L^\infty(\Omega)^{3\times 3}$ yields a uniformly elliptic viscoelastic Stokes problem, and the question then is to characterize the tensors $\mathsf{B}$ which are ${\mathcal A}$-positive. As we just mentioned, this is done in Section \ref{conditionsellipticity}.

\mathbf{e}gin{rmk}[Entropy condition]
The ellipticity conditions in Definitions \ref{viscoelasticuniformlyelliptic} and \ref{defproper} are closely related to the entropy principle of the second law of thermodynamics (see \cite{Liu2008}), which concerns the viscous part of the stress tensor and, in our model, reads
\[ (\mathsf{D}{\mathcal A}(\mathsf{B}) + {\mathcal A}(\mathsf{B})\mathsf{D}) \cdot \mathsf{D} \geq 0.
\]
The only difference is that the condition of uniform ellipticity requires a strict inequality, uniformly on the material.
\end{rmk}

Once we find a weak solution, we may recover the pressure using the following classical result (see \cite[Proposition 1.1.1]{temam}).

\mathbf{e}gin{prop}
  \label{proprecoverpressure}
  Let $\Omega\subset {\mathbb{R}}^d$ be open and $\mathbf{g}\in {\mathcal D}'(\Omega)^d$, $d\in {\mathbb{N}}$. A necessary and sufficient condition for having $\mathbf{g} = \mathbf{n}abla p$, in the distribution sense, for some $p\in {\mathcal D}'(\Omega)$ is that 
  \[ \dual{\mathbf{g},\mathbf{v}}_{{\mathcal D}'(\Omega)^d,{\mathcal D}(\Omega)^d} = 0, \quad \forall \mathbf{v}\in {\mathcal V}(\Omega).
  \]
\end{prop}

The regularity of $p$ can be inferred from the regularity of $\mathbf{g}=\mathbf{n}abla p$ as follows (see \cite[Proposition 1.1.2]{temam} and \cite{necas}).

\mathbf{e}gin{prop}
  \label{regpressurel2}
  Let $\Omega\subset{\mathbb{R}}^d$ be an open bounded Lipschitz domain, with $d\in {\mathbb{N}}$. Suppose $p\in {\mathcal D}(\Omega)'$ is a distribution with $\mathbf{n}abla p \in H^{-1}(\Omega)^3$. Then, $p\in L^2(\Omega)$ and
  \mathbf{e}gin{equation}
    \label{pressureestl2}
    \|p\|_{L^2(\Omega)/{\mathbb{R}}} \leq c(\Omega) \|\mathbf{n}abla p\|_{H^{-1}(\Omega)},
  \end{equation}
  where $c(\Omega)$ is a shape constant.
\end{prop}

Once $\mathbf{n}abla p$ is more regular, higher regularity for $p$ itself follows directly. In fact, if $\mathbf{n}abla p$ belongs to $L^2(\Omega)^3$, then clearly $p\in H^1(\Omega)$, and so on for higher regularities of $\mathbf{n}abla p$.

\section{The generalized Stokes problem}
\label{secgenstokes}

In this section, we study problem \eqref{generalizedStokes} for a given tensor $\mathsf{A}$ in $L^\infty(\Omega)^{3\times 3}$, which is assumed to be symmetric and uniformly positive definite in the sense of Definition \ref{defuniformpositivedef}. As we will see later on (Section \ref{subsecunifelliptic}), this ensures that the generalized Stokes problem \eqref{weakgeneralizedStokes} is uniformly elliptic in the sense of Definition \ref{Stokesuniformlyelliptic}. We then use this to obtain the existence and uniqueness of the solutions of this problem (Section \ref{subsecexistuniqgenStokes}) and higher-order regularity (Sections \ref{subsecreggenStokes} and \ref{subsechighreggenStokes}).

\subsection{Uniform ellipticity of the generalized Stokes problem}
\label{subsecunifelliptic}

We start with the following useful lemma.
\mathbf{e}gin{lem}
  \label{lempositivityofL}
  Let $\mathsf{A}\in L^\infty(\Omega)^{3\times 3}$ be a symmetric tensor on $\Omega$ which is uniformly positive definite and let $\alpha>0$ be the corresponding constant such that \eqref{unifellip} holds. For $\mathsf{M}\in L^\infty(\Omega)^{3\times 3}$, consider the operator
  \mathbf{e}gin{equation}
    \label{defopL}
    \mathsf{L}(\mathsf{M}) = \mathsf{S}(\mathsf{M})\mathsf{A} + \mathsf{A}\mathsf{S}(\mathsf{M}),
  \end{equation} 
  where $\mathsf{S}$ is the symmetrization operator
  \mathbf{e}gin{equation}
    \label{defopS}
    \mathsf{S}(\mathsf{M}) = \frac{1}{2}(\mathsf{M} + \mathsf{M}^\text{\rm tr}).
  \end{equation}
  Then, $\mathsf{L}$ satisfies
  \mathbf{e}gin{equation}
    \label{symposopL}
    \mathsf{L}(\mathsf{M}): \mathsf{M} \geq 2\alpha \mathsf{S}(\mathsf{M}):\mathsf{S}(\mathsf{M}), \quad \forall \mathsf{M}\in L^\infty(\Omega)^{3\times 3},
  \end{equation}
  everywhere on $\Omega$. 
\end{lem}

\mathbf{e}gin{proof}
  Since $\mathsf{A}$ and $\mathsf{S}(\mathsf{M})$ are symmetric, so is $\mathsf{L}(\mathsf{M})$. Thus, omitting $\mathbf{x}\in \Omega$ for the sake of simplicity,
  \[ \mathsf{L}(\mathsf{M}) : \mathsf{M} = \frac{1}{2}\mathsf{L}(\mathsf{M}):\mathsf{M} + \frac{1}{2}\mathsf{L}(\mathsf{M})^\text{\rm tr}: \mathsf{M} = \frac{1}{2}\mathsf{L}(\mathsf{M}):\mathsf{M} + \frac{1}{2}\mathsf{L}(\mathsf{M}): \mathsf{M}^\text{\rm tr} = \mathsf{L}(\mathsf{M}): \mathsf{S}(\mathsf{M}).
  \]
  Write $\mathsf{S}(\mathsf{M}) = (s_{ij})_{ij}$ and $\mathsf{A}=(a_{ij})_{ij}$. Then,
\[ \mathsf{L}(\mathsf{M}) = \mathsf{S}(\mathsf{M})\mathsf{A}+\mathsf{A}\mathsf{S}(\mathsf{M}) = \left( \sum_l (s_{il}a_{lj} + a_{il}s_{lj})\right)_{ij}.
\]
Using again that $\mathsf{A}$ and $\mathsf{S}(\mathsf{M})$ are symmetric,
\mathbf{e}gin{multline*}
  \mathsf{L}(\mathsf{M}): \mathsf{M} = \mathsf{L}(\mathsf{M}):\mathsf{S}(\mathsf{M}) = \sum_{ijl} (s_{il}a_{lj} + a_{il}s_{lj}) s_{ij} \\
  = \sum_{ijl} (a_{lj}s_{il}s_{ij} + a_{li} s_{lj} s_{ij}) = \sum_{ijl} (a_{lj}s_{il}s_{ij} + a_{li} s_{jl} s_{ji}).
\end{multline*}
Switching $i$ with $j$ in the second term in the right hand side above and using the uniform strict positivity of $\mathsf{A}$ we find that
\[ \mathsf{L}(\mathsf{M}): \mathsf{M} = 2 \sum_{ijl} a_{lj}s_{il}s_{ij} \geq 2\alpha \sum_{ij} s_{ij}^2 = 2\alpha \mathsf{S}(\mathsf{M}): \mathsf{S}(\mathsf{M}),
\]
which completes the proof.
\end{proof}

\mathbf{e}gin{rmk}
  \label{rmkopLnotpositive}
  Lemma \ref{lempositivityofL} says that $\mathsf{L}$ is uniformly strictly positive on the space of essentially bounded symmetric tensors, which we may denote by $L_{\textrm{sym}}^\infty(\Omega)^{3\times 3}$. Indeed, from \eqref{symposopL} and using that $\mathsf{S}(\mathsf{M}) = \mathsf{M}$ for a symmetric tensor, we have
  \[ \mathsf{L}(\mathsf{M}): \mathsf{M}  \geq 2\alpha \mathsf{M}: \mathsf{M}, \quad \forall \mathsf{M}\in L_{\textrm{sym}}^\infty(\Omega)^{3\times 3}.
  \]
  However, $\mathsf{L}$ is not strictly positive on all $L^\infty(\Omega)^{3\times 3}$. Indeed, for any $\mathsf{M}$ anti-symmetric, we have $\mathsf{S}(\mathsf{M}) = 0$, so that $\mathsf{L}(\mathsf{M}) = 0$. In this work, we are actually interested in the case that $\mathsf{M} = \mathbf{n}abla \mathbf{v}$, which can be written as $\mathsf{M} = (\mathbf{n}abla \otimes \mathbf{v})^\text{\rm tr}$. If we extrapolate on this idea and look at tensors of the form $\mathsf{M}=\mathbf{x}i\otimes\boldsymbol{\eta}=(\xi_i\eta_j)_{ij},$ for vectors $\mathbf{x}i$ and $\boldsymbol{\eta}$, then we notice that
  \[ \mathsf{S}(\mathsf{M}) = \mathsf{S}(\mathbf{x}i\otimes\boldsymbol{\eta}) = \frac{1}{2}(\mathbf{x}i\otimes\boldsymbol{\eta} + \boldsymbol{\eta}\otimes\mathbf{x}i),
  \]
  and
  \[ \mathsf{S}(\mathsf{M}):\mathsf{S}(\mathsf{M}) = \frac{1}{4}(\mathbf{x}i\otimes\boldsymbol{\eta} + \boldsymbol{\eta}\otimes\mathbf{x}i): (\mathbf{x}i\otimes\boldsymbol{\eta} + \boldsymbol{\eta}\otimes\mathbf{x}i)
 = \frac{1}{2}\left(|\mathbf{x}i|^2|\boldsymbol{\eta}|^2 + (\mathbf{x}i\cdot\boldsymbol{\eta})^2 \right) = \frac{1}{2}\mathsf{M}:\mathsf{M} + \frac{1}{2}(\text{\rm tr}\mathsf{M})^2.
  \]
  This expression is directly related to the identity \eqref{Dandgradvinhzero1}. In particular,
  \[  \mathsf{S}(\mathsf{M}):\mathsf{S}(\mathsf{M})  \geq \frac{1}{2}\mathsf{M}:\mathsf{M}.
  \]
  Thus,
  \[ \mathsf{L}(\mathsf{M}): \mathsf{M} \geq \alpha \mathsf{M}:\mathsf{M},
  \]
  for any $\mathsf{M}=\mathbf{x}i\otimes\boldsymbol{\eta}$, so that $\mathsf{L}$ is uniformly strictly positive on the space of tensors of this form.
\end{rmk}

The motivation for considering the operators $\mathsf{L}$ and $\mathsf{S}$ above follows from the fact that, by taking $\mathsf{M}=\mathbf{n}abla \mathbf{v}$, we see that
\mathbf{e}gin{equation}
  \label{DvasSnablav} 
  \mathsf{D}(\mathbf{v}) = \mathsf{S}(\mathbf{n}abla\mathbf{v}),
\end{equation}
and
\mathbf{e}gin{equation}
  \label{linearpartasL}
  \mathsf{D}(\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v}) = \mathsf{L}(\mathbf{n}abla \mathbf{v}),
\end{equation}
so that the generalized Stokes problem \eqref{generalizedStokes} can be written as
\mathbf{e}gin{equation}
  \label{generalizedStokeswithL}
   \mathbf{e}gin{cases}
      - \mathbf{n}abla \cdot (\mathsf{L}(\mathbf{n}abla \mathbf{v})) + \mathbf{n}abla p = \mathbf{b}f, & \text{in } \Omega,\\
      \mathbf{n}abla \cdot \mathbf{v} = 0, & \text{in } \Omega,\\
      \mathbf{v} = 0, & \text{on } \partial \Omega.
    \end{cases}
\end{equation}
and similarly for the elliptic problem \eqref{weakgeneralizedStokes}. With that in mind, we obtain, using Lemma \ref{lempositivityofL}, the uniform ellipticity of the Problem \eqref{weakgeneralizedStokes}.

\mathbf{e}gin{prop}
  \label{probgenStokesisunifelliptic}
  If $\mathsf{A}\in L^\infty(\Omega)^{3\times 3}$ is a symmetric tensor on $\Omega$ which is uniformly positive definite then the associated problem \eqref{weakgeneralizedStokes} is uniformly elliptic in $H_0^1(\Omega)^3$. More precisely, if \eqref{unifellip} holds for $\mathsf{A}$ with $\alpha>0$, then \eqref{probgenStokesunifellip} holds with $\delta=\alpha$ in either $V(\Omega)$ and $H_0^1(\Omega)^3$ and with $\delta' = \sqrt{2}\|\mathsf{A}\|_{L^\infty(\Omega)}$, in $V(\Omega)$, and $\delta' = 2\|\mathsf{A}\|_{L^\infty(\Omega)}$, in $H_0^1(\Omega)^3$.
\end{prop}

\mathbf{e}gin{proof}
   Taking \eqref{linearpartasL} and \eqref{DvasSnablav} into consideration and applying Lemma \ref{lempositivityofL} we find from \eqref{symposopL} that
  \[ (\mathsf{D}(\mathbf{v}) \mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v})) : \mathbf{n}abla \mathbf{v} = \mathsf{L}(\mathbf{n}abla\mathbf{v}): \mathbf{n}abla\mathbf{v} \geq 2\alpha \mathsf{S}(\mathbf{n}abla\mathbf{v}) : \mathsf{S}(\mathbf{n}abla \mathbf{v}) = 2\alpha \mathsf{D}(\mathbf{v}): \mathsf{D}(\mathbf{v}).
  \]
  Integrating over $\Omega$ and using the lower bound in \eqref{Dandgradvcoerciveinhzero1} we find that
  \[ \int_\Omega (\mathsf{D}(\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v})) : \mathbf{n}abla \mathbf{v} \;{\text{\rm d}} \mathbf{x} \geq  2\alpha \|\mathsf{D}(\mathbf{v})\|_{L^2(\Omega)}^2 \geq \alpha \|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)}^2,
  \]
for every $\mathbf{v}$ in $H_0^1(\Omega)^3$. This proves the lower bound in \eqref{probgenStokesunifellip} with $\delta = \alpha$, in $H_0^1(\Omega)^3$, and, in particular, in $V(\Omega)$.

For the upper bound in \eqref{probgenStokesunifellip}, we use H\"older's inequality and the assumption that $\mathsf{A}\in L^\infty(\Omega)^{3\times 3}$ to find that
\[ \int_\Omega (\mathsf{D}(\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v})) : \mathbf{n}abla \mathbf{v} \;{\text{\rm d}} \mathbf{x} 
  \leq 2\|\mathsf{A}\|_{L^\infty(\Omega)}\|\mathsf{D}(\mathbf{v})\|_{L^2(\Omega)}\|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)},
\]
for all $\mathbf{v}\in H_0^1(\Omega)^3$. Then, using the upper bound in \eqref{Dandgradvcoerciveinhzero1} for $\mathsf{D}(\mathbf{v})$ yields the upper bound in \eqref{probgenStokesunifellip}, in $H_0^1(\Omega)^3$, with $\delta'=2\|\mathsf{A}\|_{L^\infty(\Omega)}$, while the upper bound in \eqref{Dandgradv} for $\mathsf{D}(\mathbf{v})$ yields the upper bound in \eqref{probgenStokesunifellip}, in $V(\Omega)$, with $\delta'=\sqrt{2}\|\mathsf{A}\|_{L^\infty(\Omega)}$.
\end{proof}

\subsection{Existence and uniqueness of solutions of the weak form of the generalized Stokes problem}
\label{subsecexistuniqgenStokes}

Thanks to the ellipticity property proved in Proposition \ref{probgenStokesisunifelliptic}, an application of the Lax-Milgram theorem solves the question of existence and uniqueness of the weak formulation of the generalized Stokes problem.

\mathbf{e}gin{prop}
  \label{weaksolutionweakgeneralizedStokes}
  Let $\Omega\subset {\mathbb{R}}^3$ be a bounded domain, let $\mathsf{A}\in L^\infty(\Omega)^{3\times 3}$ be a uniformly positive definite symmetric tensor on $\Omega$ and let $\mathbf{b}f\in H^{-1}(\Omega)^3$. Then, there exists a unique $\mathbf{v}\in V(\Omega)$ that solves the weak formulation \eqref{weakgeneralizedStokes} of the generalized Stokes problem. Moreover, if $\alpha>0$ is such that \eqref{unifellip} holds, then
  \mathbf{e}gin{equation}
    \label{weaksolutionweakgeneralizedStokesestimate}
    \|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)} \leq \frac{1}{\alpha} \|\mathbf{b}f\|_{H^{-1}(\Omega)}.
  \end{equation}
\end{prop}

\mathbf{e}gin{proof}
  Consider the bilinear form $a: V(\Omega)\times V(\Omega) \rightarrow {\mathbb{R}}$ given by
  \[ a(\mathbf{v},\mathbf{w}) = \int_\Omega (\mathsf{D}(\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v})) : \mathbf{n}abla \mathbf{v} \;{\text{\rm d}} \mathbf{x}, \quad \forall \mathbf{v}, \mathbf{w} \in V(\Omega),
  \]
  and the linear form $\ell:V(\Omega)\rightarrow {\mathbb{R}}$ given by
  \[ \ell(\mathbf{w}) = \int_\Omega \mathbf{b}f\cdot \mathbf{w} \;{\text{\rm d}}\mathbf{x}, \quad \forall \mathbf{w}\in V(\Omega).
  \]
  Since $\mathsf{A}\in L^\infty(\Omega)^{3\times 3}$, it follows that $a(\cdot,\cdot)$ is a well-defined and continuous bilinear form on $V(\Omega)$. Moreover, thanks to Proposition \ref{probgenStokesisunifelliptic}, the estimate \eqref{probgenStokesunifellip} holds in $V(\Omega)$ with $\delta=\alpha$, and $a(\cdot,\cdot)$ is a coercive bilinear form on $V(\Omega)$. On the other hand, since $\mathbf{b}f\in H^{-1}(\Omega)^3\subset V'(\Omega)$, then $\ell(\cdot)$ is a continuous linear form on $V(\Omega)$. Thus, it follows from the Lax-Milgram Theorem (see e.g. \cite{laxmilgram}) that there exists one, and only one, $\mathbf{v}\in V(\Omega)$, such that $a(\mathbf{v},\mathbf{w}) = \ell(\mathbf{w})$, for all $\mathbf{w}\in V(\Omega)$, which proves the existence and uniqueness part of the statement. For the estimate, we use the ellipticity condition \eqref{probgenStokesunifellip} to write, with $\mathbf{w}=\mathbf{v}$,
  \[ \alpha \|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)}^2 \leq a(\mathbf{v},\mathbf{v}) = \ell(\mathbf{v}) = \int_\Omega \mathbf{b}f\cdot \mathbf{v} \;{\text{\rm d}}\mathbf{x} \leq \|\mathbf{b}f\|_{H^{-1}(\Omega)}\|\mathbf{v}\|_{H_0^1(\Omega)} = \|\mathbf{b}f\|_{H^{-1}(\Omega)}\|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)},
  \]
  so that
  \[ \alpha \|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)} \leq \|\mathbf{b}f\|_{H^{-1}(\Omega)},
  \]
  which yields \eqref{weaksolutionweakgeneralizedStokesestimate}. This completes the proof.
\end{proof}

\mathbf{e}gin{rmk}
  \label{rmkweaksolutionweakgeneralizedStokesdualforces}
  In the case $\mathbf{b}f\in L^2(\Omega)^3$ in Proposition \ref{weaksolutionweakgeneralizedStokes}, the estimate \eqref{weaksolutionweakgeneralizedStokesestimate} implies  $\alpha \|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)} \leq \lambda_1^{-1/2}\|\mathbf{b}f\|_{L^2(\Omega)}$, where $\lambda_1$ is the constant in the Poincar\'e inequality \eqref{poincareineq}.
\end{rmk}

Using Proposition \ref{weaksolutionweakgeneralizedStokes}, along with Propositions \ref{proprecoverpressure} and \ref{regpressurel2}, we recover the pressure and obtain the following existence and uniqueness result.
\mathbf{e}gin{thm}[Existence and uniqueness]
  \label{weaksolutionweakgeneralizedStokeswithpressure}
  Let $\Omega\subset {\mathbb{R}}^3$ be a bounded Lipschitz domain, $\mathsf{A}\in L^\infty(\Omega)^{3\times 3}$ be a uniformly positive definite symmetric tensor on $\Omega$, and $\mathbf{b}f\in H^{-1}(\Omega)^3$. Then, there exist a unique vector field $\mathbf{v}\in V(\Omega)$ and a scalar field $p\in L^2(\Omega)$ that solve the weak formulation \eqref{weakgeneralizedStokeswithpressure} of the generalized Stokes problem. The field $p$ is unique up to a constant, in $L^2(\Omega)$. Moreover, if $\alpha>0$ is such that \eqref{unifellip} holds, then
  \mathbf{e}gin{equation}
    \label{weaksolutionweakgeneralizedStokesestimaterepeat}
    \|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)} \leq \frac{1}{\alpha} \|\mathbf{b}f\|_{H^{-1}(\Omega)},
  \end{equation}
  and  
  \mathbf{e}gin{equation}
    \label{weaksolutionweakgeneralizedStokesestimatewithpressure}
    \|p\|_{L^2(\Omega)/{\mathbb{R}}} \leq \frac{c(\Omega)}{\alpha} \|\mathsf{A}\|_{L^\infty(\Omega)}\|\mathbf{b}f\|_{H^{-1}(\Omega)},
  \end{equation}
  where $c(\Omega)$ is a shape constant.
\end{thm}

\mathbf{e}gin{proof}
  From Proposition \ref{weaksolutionweakgeneralizedStokes} we have a unique solution $\mathbf{v}\in V(\Omega)$ of the generalized Stokes problem \eqref{weakgeneralizedStokes}, which can be written as
  \[ \dual{-\mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v})) - \mathbf{b}f, \mathbf{w} }_{{\mathcal D}'(\Omega)^3, {\mathcal D}(\Omega)^3} = 0, \qquad \forall \mathbf{w}\in {\mathcal V}(\Omega).
  \]
  Then, from Propositions \ref{proprecoverpressure} and \ref{regpressurel2}, there is a distribution $p\in L^2(\Omega)/{\mathbb{R}}$, which is clearly unique up to a constant, such that
  \[ \mathbf{n}abla p = -\mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v})) - \mathbf{b}f,
  \]
  with
  \[ \|p\|_{L^2(\Omega)/{\mathbb{R}}} \leq c(\Omega) \|\mathbf{n}abla(\mathsf{D}(\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v})) + \mathbf{b}f\|_{H^{-1}(\Omega)}.
  \]
  Thus,
  \[ \|p\|_{L^2(\Omega)/{\mathbb{R}}} \leq c(\Omega) \left(\|\mathsf{D}(\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v})\|_{L^2(\Omega)} + \|\mathbf{b}f\|_{H^{-1}(\Omega)}\right).
  \]
  In three dimensions, we estimate
  \[ \|\mathsf{D}(\mathbf{v})\mathsf{A}\|_{L^2(\Omega)} \leq \sqrt{3} \|\mathsf{A}\|_{L^\infty(\Omega)}\|\mathsf{D}(\mathbf{v})\|_{L^2(\Omega)},
  \]
  and similarly for $\mathsf{A}\mathsf{D}(\mathbf{v})$. Then, using \eqref{Dandgradv}, we find
  \[ \|p\|_{L^2(\Omega)/{\mathbb{R}}}  \leq c(\Omega)\left(\sqrt{3} \|\mathsf{A}\|_{L^\infty(\Omega)}\|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)} + \|\mathbf{b}f\|_{H^{-1}(\Omega)}\right). 
  \]
  Using the previous estimate \eqref{weaksolutionweakgeneralizedStokesestimaterepeat} for $\mathbf{n}abla \mathbf{v}$ we obtain
  \[ \|p\|_{L^2(\Omega)/{\mathbb{R}}} \leq c(\Omega)\left(\frac{\sqrt{3} \|\mathsf{A}\|_{L^\infty(\Omega)}}{\alpha} + 1\right)\|\mathbf{b}f\|_{H^{-1}(\Omega)}.
  \]  
Using \eqref{alphaboundedLinfty} we finally arrive at \eqref{weaksolutionweakgeneralizedStokesestimatewithpressure}, with another shape constant $c(\Omega)$.
\end{proof}

\subsection{Regularity of the solutions of the generalized Stokes problem}
\label{subsecreggenStokes}

As usual, if we have more regularity on the data, then we obtain more regularity on the solution. 

The first regularity result that we consider follows from the paper of Mach\'a \cite[Theorem 4.2]{macha2011}, about a general class of systems that includes ours, or from the result of Giaquinta and Modica \cite[Part II, Theorem 1.2 in page 198 and Remark 1.5 in page 199]{giaquintamodica1982}, about a slightly more general system, and recalling that, with smooth boundaries, we have $W^{1,\infty}={\mathcal C}^1$. See also the related result of Hyu and Star\'a \cite[Theorem 4.1]{huystara2006}, which unfortunately requires that the operator $\mathsf{L}$, defined in \eqref{defopL}, be strictly positive for any $\mathsf{M}$ and not only for symmetric tensors $\mathsf{M}$ (see Remark \ref{rmkopLnotpositive}), although this is clearly unnecessary as discussed in \cite[Part II, Remark 1.5]{giaquintamodica1982}. The proof of those results is based on a now classical idea of estimating the finite-difference quotients of the solution, considered for instance in \cite{adn2}, and used here in the proof of higher-regularity in Theorem \ref{hkplustworegularity} (see also \cite[Theorem 8.12, page 186]{gilbargtrudinger2001} for a result with similar hypotheses in a classical elliptic problem).

\mathbf{e}gin{thm}[Regularity]
  \label{h2regularity}
  Let $\Omega$ be a bounded domain in ${\mathbb{R}}^3$ with a boundary of class ${\mathcal C}^2$ and suppose $\mathsf{A}\in W^{1,\infty}(\Omega)^{3\times 3}$ is a uniformly positive definite symmetric tensor on $\Omega$ and $\mathbf{b}f \in L^2(\Omega)^3$. Then, the solution $(\mathbf{v}, p)$ of the weak formulation \eqref{weakgeneralizedStokeswithpressure} of the generalized Stokes problem is such that $\mathbf{v} \in V(\Omega)\cap H^2(\Omega)$ and $p\in H^1(\Omega)$, with 
  \mathbf{e}gin{equation}
    \label{h2regularityforh2v}
    \|D^2 \mathbf{v}\|_{L^2(\Omega)} \leq \frac{c(\Omega)}{\alpha}\left( \|\mathbf{b}f\|_{L^2(\Omega)} + \frac{1}{\alpha}\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}\|\mathbf{b}f\|_{H^{-1}(\Omega)}\right),
  \end{equation}
  and
  \mathbf{e}gin{equation}
    \label{h2regularityforh1p}
    \|\mathbf{n}abla p\|_{L^2(\Omega)} \leq \frac{c(\Omega)}{\alpha}\|\mathsf{A}\|_{L^\infty(\Omega)}\left( \|\mathbf{b}f\|_{L^2(\Omega)} + \frac{1}{\alpha}\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}\|\mathbf{b}f\|_{H^{-1}(\Omega)}\right),
  \end{equation}
  where $\lambda_1$ is the Poincar\'e constant and $c(\Omega)$ is a shape constant.  Moreover, $(\mathbf{v},p)$ solves the generalized Stokes problem \eqref{generalizedStokes} almost everywhere in $\Omega$, with the first two equations in \eqref{generalizedStokes} holding in $L^2(\Omega)^3$ and $H^1(\Omega)^3$, respectively, and with the last equation holding for the trace of $\mathbf{v}$ on $\partial\Omega$ in the space $H^{3/2}(\partial \Omega)$.
\end{thm}

\mathbf{e}gin{proof} 
  We apply the result \cite[Theorem 4.2 with $k=0$]{macha2011} to obtain the desired regularity. In \cite{macha2011}, the system considered is of the form \eqref{generalizedStokeswithL} (even more general, allowing the divergence to be any given function in $H^1(\Omega)$). It is assumed that $\mathsf{L}$ is weakly coercive in the sense defined in \cite{macha2011}, which is precisely the lower bound in \eqref{probgenStokesunifellip}, which, as proved in Proposition \ref{probgenStokesisunifelliptic}, holds under our hypothesis that $\mathsf{A}$ is a uniformly positive definite symmetric tensor. It is also assumed that the boundary of $\Omega$ is of class ${\mathcal C}^2$ and that $\mathbf{b}f \in L^2(\Omega)^3$. Finally, it is assumed that $\mathsf{L}\in W^{1,\infty}(\Omega)^{3\times 3}$, which follows immediately from the assumption that $\mathsf{A}\in W^{1,\infty}(\Omega)^{3\times 3}$ and the fact, which is clear from \eqref{defopL}, that $\mathsf{L}$ has the same regularity as $\mathsf{A}$. Hence, we find that $\mathbf{v}\in H^2(\Omega)^3$ and $p\in H^1(\Omega)$. In particular, we deduce, from the weak formulation, that the strong form \eqref{generalizedStokes} holds almost everywhere, with the regularities as we described in the statement of this theorem.
  
  Now that we know that the solution is more regular, we address the Sobolev estimates for the solution in terms of the data. 
  
  
  Since we now know that $\mathbf{v}$ is in $H^2(\Omega)^3$, we have that $\mathbf{n}abla \mathbf{v}$ is in $H^1(\Omega)$ and hence it has sufficient regularity to be a test function in the weak formulation, but it lacks the condition of vanishing on the boundary. Nevertheless, since $\mathbf{v}$ vanishes on the boundary, the tangential derivatives of $\mathbf{v}$ along the boundary also vanish, i.e. $\mathbf{n}abla\mathbf{v}(\mathbf{x})\cdot \mathbf{x}i = 0$, for every $\mathbf{x}i$ tangent to $\Omega$ at $\mathbf{x}\in \partial\Omega$. This will be sufficient for our purposes.  
  
  The idea is to look at the equations for the partial derivatives of $\mathbf{v}$. By taking the distributional derivative of \eqref{generalizedStokes} with respect to $x_i$, $i=1,2,3$, we find the equations
\mathbf{e}gin{equation}
  \label{generalizedStokesforderivative}
   \mathbf{e}gin{cases}
      - \mathbf{n}abla \cdot (\mathsf{D}(\partial_{x_i}\mathbf{v})\mathsf{A} + \mathsf{A}\mathsf{D}(\partial_{x_i}\mathbf{v})) -\mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v})(\partial_{x_i}\mathsf{A}) + (\partial_{x_i}\mathsf{A})\mathsf{D}(\mathbf{v})) \\
      \qquad \qquad \qquad + \mathbf{n}abla (\partial_{x_i}p) =  \partial_{x_i}\mathbf{b}f, & \text{in } \Omega,\\
      \mathbf{n}abla \cdot (\partial_{x_i}\mathbf{v}) = 0, & \text{in } \Omega,
    \end{cases}
\end{equation}
valid in the distribution sense on $\Omega$. Since we already know that $\mathbf{v} \in V(\Omega)\cap H^2(\Omega)$ and $p\in H^1(\Omega)$, then the first equation in \eqref{generalizedStokesforderivative} actually holds in $H^{-1}(\Omega)^3$, while the second one holds in $L^2(\Omega)$. The system \eqref{generalizedStokesforderivative} can be written as
\mathbf{e}gin{equation}
  \label{generalizedStokesforderivativealt}
   \mathbf{e}gin{cases}
      - \mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v}_i)) - \mathbf{b}_i + \mathbf{n}abla p_i = \mathbf{b}f_i, & \text{in } \Omega,\\
      \mathbf{n}abla \cdot \mathbf{v}_i= 0, & \text{in } \Omega,
    \end{cases}
\end{equation}
where
\[ \mathbf{v}_i = \partial_{x_i}\mathbf{v}, \qquad p_i = \partial_{x_i}p, \qquad \mathbf{b}f_i = \partial_{x_i}\mathbf{b}f, \qquad \mathbf{b}_i = \mathbf{B}(\partial_{x_i}\mathsf{A},\mathbf{v}),
\]
with
\mathbf{e}gin{equation}
  \label{defbB}
  \mathbf{B}(\mathsf{M},\mathbf{w}) =  \mathbf{n}abla \cdot (\mathsf{D}(\mathbf{w})\mathsf{M} + \mathsf{M}\mathsf{D}(\mathbf{w})),
\end{equation}
for any tensor $\mathsf{M}$ and any vector field $\mathbf{w}$.

Equation \eqref{generalizedStokesforderivativealt} has the same form of \eqref{generalizedStokes}, if viewing the additional lower order term $\mathbf{b}_i$ as an extra forcing term. Thus, if we had the same no-slip boundary conditions for the derivatives $\mathbf{v}_i = \partial_{x_i}\mathbf{v}$, Theorem \ref{weaksolutionweakgeneralizedStokeswithpressure} would apply to \eqref{generalizedStokesforderivativealt} with the right hand side as $\mathbf{b}f_i + \mathbf{b}_i$, and estimates \eqref{weaksolutionweakgeneralizedStokesestimaterepeat} and \eqref{weaksolutionweakgeneralizedStokesestimatewithpressure} would yield estimates for $\|\mathbf{n}abla\mathbf{v}_i\|_{L^2(\Omega)}$ and $\|p_i\|_{L^2(\Omega)/{\mathbb{R}}}$ in terms of $\alpha$, $\|\mathbf{b}f_i + \mathbf{b}_i\|_{H^{-1}(\Omega)}$ and $\|\mathsf{A}\|_{L^\infty(\Omega)}$. These estimates, together  with an estimate for $\mathbf{b}_i = \mathbf{B}(\partial_{x_i}\mathsf{A},\mathbf{v})$ using the bound for $\partial_{x_i}\mathsf{A}$ in $L^\infty$ and the estimate \eqref{weaksolutionweakgeneralizedStokesestimaterepeat} for $\mathbf{n}abla \mathbf{v}$, would yield estimates similar to \eqref{h2regularityforh2v} and \eqref{h2regularityforh1p}, with the differences that $c$ above is independent of $\Omega$ and only the higher order term $\mathbf{n}abla \mathsf{A}$ of $\mathsf{A}$ would appear. 

However, of course, $\mathbf{v}_i = \partial_{x_i}\mathbf{v}$ does not vanish on the boundary, so the above argument does not apply. Nevertheless, the essence of the argument can be recovered by localizing the equations to the interior and to the boundary of the domain and using that the derivatives of $\mathbf{v}$ tangential to the boundary vanish. At the end, this will give us \eqref{h2regularityforh2v} and \eqref{h2regularityforh1p} with the lower order terms of $\mathsf{A}$ and the multiplicative shape constant $c(\Omega)$.
  
More precisely, we use a partition of unit based on a finite covering of the boundary by open sets $U_m$, $m=1, \ldots, M$, for some $M\in{\mathbb{N}}$, and a complementing open set $U_0$ compactly included in $\Omega$, so that $\Omega\subset\bigcup_{m=0}^M U_m.$ We write $\mathbf{v}_i=\sum_m \eta_m\mathbf{v}_i$, where each $\eta_m$ is supported on $U_m$, with $0\leq \eta_m\leq 1$, and $\sum_m\eta_m = 1$ on $\Omega$. We then estimate each portion $\eta_m\mathbf{v}_i$. 
  
The idea of the proof is that, for the component $m=0$, $\eta_0$ is supported on the interior of the domain, so that $\eta_0\mathbf{v}$ vanishes on a neighborhood of the boundary, and so does $\eta_0^2\partial_i\mathbf{v}$, for every $i=1,2,3$. Moreover, since we already know that $\mathbf{v}\in H^2(\Omega)^3$, then $\eta_0^2\partial_i\mathbf{v}$ belongs to $H_0^1(\Omega)^3$. Thus, the estimate for $\eta_0\partial_{x_i}\mathbf{v}$ is obtained directly from the equation \eqref{generalizedStokesforderivative}, by testing it with $\eta_0^2\partial_i\mathbf{v}$. 

For the boundary components $m=1,\ldots, M$, we first rectify the domain to a hemisphere by a ${\mathcal C}^2$ change of variables, transforming the system into a similar system with coefficients with the same regularity as $\mathbf{A}$, and then use that, in the rectified system, all the derivatives in the directions tangential to the boundary are well defined and do vanish on the boundary, which allows us to use the ellipticity of the equation and the weak formulation of \eqref{generalizedStokesforderivative} tested with a localized version of each tangential derivative to get uniform bounds for the tangential derivatives. Then, we use the differential equation itself to express the derivative in the direction normal to the boundary in terms of the derivatives in the directions tangential to the boundary and other lower order derivatives (just like we write $v_{zz} = - f - v_{xx}-v_{yy}$ over the half-plane domain $z<0$, in the case of the elliptic equation $-\Delta v = f$), hence obtaining an estimate of the same order for the derivative normal to the boundary as well. Let us go into that with more details. And we work with \eqref{generalizedStokesforderivativealt} for notational simplicity. 

For the interior regularity, we test the equation \eqref{generalizedStokesforderivativealt}, valid in $H^{-1}(\Omega)^3$, with the function $\mathbf{w} = \eta_0^2 \mathbf{v}_i$, which belongs to $H_0^1(\Omega)^3$, so that
  \mathbf{e}gin{multline}
    \label{weakgeneralizedStokeswithpressureinterior} 
    \int_\Omega (\mathsf{D}(\mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v}_i)) \cdot \mathbf{n}abla (\eta_0^2\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x}  + \int_{\Omega}  (\mathsf{D}(\mathbf{v})(\partial_{x_i}\mathsf{A}) + (\partial_{x_i}\mathsf{A})\mathsf{D}(\mathbf{v}))\cdot \mathbf{n}abla (\eta_0^2\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x} \\
    - \int_\Omega  p_i \mathbf{n}abla \cdot (\eta_0^2\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x} = \int_\Omega \mathbf{b}f_i \cdot (\eta_0^2\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x}.
  \end{multline}

Notice that $\mathbf{n}abla(\eta_0^2\mathbf{v}_i) = \eta_0 \mathbf{n}abla(\eta_0\mathbf{v}_i) + \eta_0(\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i$, and 
$D(\eta_0\mathbf{v}_i) = \eta_0D(\mathbf{v}_i) + \mathsf{S}((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i)$, so that the first term in \eqref{weakgeneralizedStokeswithpressureinterior} can be written as
\mathbf{e}gin{align*}
  \int_\Omega (\mathsf{D}(\mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v}_i)) & \cdot \mathbf{n}abla (\eta_0^2\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x} \\
    & = \int_\Omega \eta_0(\mathsf{D}(\mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v}_i)) \cdot \mathbf{n}abla (\eta_0\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x}  \\
  & \qquad + \int_\Omega \eta_0(\mathsf{D}(\mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v}_i)) \cdot ((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i) \;{\text{\rm d}} \mathbf{x} \\
  & =  \int_\Omega (\mathsf{D}(\eta_0\mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{D}(\eta_0\mathbf{v}_i)) \cdot \mathbf{n}abla (\eta_0\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x} \\
  & \qquad - \int_\Omega (\mathsf{S}((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{S}((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i)) \cdot \mathbf{n}abla (\eta_0\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x}  \\
  & \qquad + \int_\Omega (\mathsf{D}(\eta_0\mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{D}(\eta_0\mathbf{v}_i)) \cdot ((\mathbf{n}abla\eta_0)\otimes (\eta_0\mathbf{v}_i)) \;{\text{\rm d}} \mathbf{x} \\
  & \qquad - \int_\Omega (\mathsf{S}((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{S}((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i)) \cdot ((\mathbf{n}abla\eta_0)\otimes (\eta_0\mathbf{v}_i)) \;{\text{\rm d}} \mathbf{x}.
\end{align*}
Since $p_i = \partial_{x_i}p$ and $\mathbf{v}_i$ is divergence-free, we rewrite the pressure term as
\[ - \int_\Omega  \partial_{x_i} p \mathbf{n}abla \cdot (\eta_0^2\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x} = 2\int_\Omega  p \partial_{x_i} ((\mathbf{n}abla \eta_0) \cdot (\eta_0\mathbf{v}_i)) \;{\text{\rm d}} \mathbf{x}.
\]
Using that $\mathbf{b}f_i=\partial_{x_i}\mathbf{b}f$, we rewrite the forcing term as
\[ \int_\Omega \mathbf{b}f_i \cdot (\eta_0^2\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x} = - \int_\Omega \mathbf{b}f \cdot \partial_{x_i}(\eta_0^2\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x}.
\]
Hence, we find that
  \mathbf{e}gin{equation}
    \label{weakgeneralizedStokeswithpressureinteriorb} 
    \mathbf{e}gin{aligned}
    \int_\Omega & (\mathsf{D}(\eta_0\mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{D}(\eta_0\mathbf{v}_i)) \cdot \mathbf{n}abla (\eta_0\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x} \\
    & = - \int_{U_0}  \mathbf{b}_\mathbf{g}amma \cdot (\delta_i^{-h}(\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma)) \;{\text{\rm d}} \mathbf{x} \\
    & \qquad + \int_\Omega (\mathsf{S}((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{S}((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i)) \cdot \mathbf{n}abla (\eta_0\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x} \\
    & \qquad - \int_\Omega (\mathsf{D}(\eta_0\mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{D}(\eta_0\mathbf{v}_i)) \cdot ((\mathbf{n}abla\eta_0)\otimes (\eta_0\mathbf{v}_i)) \;{\text{\rm d}} \mathbf{x} \\
  & \qquad + \int_\Omega (\mathsf{S}((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{S}((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i)) \cdot ((\mathbf{n}abla\eta_0)\otimes (\eta_0\mathbf{v}_i)) \;{\text{\rm d}} \mathbf{x} \\
  & \qquad - \int_{\Omega}  (\mathsf{D}(\mathbf{v})(\partial_{x_i}\mathsf{A}) + (\partial_{x_i}\mathsf{A})\mathsf{D}(\mathbf{v}))\cdot \mathbf{n}abla (\eta_0^2\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x} \\
  & \qquad - 2\int_\Omega  p \partial_{x_i} ((\mathbf{n}abla \eta_0) \cdot (\eta_0\mathbf{v}_i)) \;{\text{\rm d}} \mathbf{x} \\
   & \qquad  - \int_\Omega \mathbf{b}f \cdot \partial_{x_i}(\eta_0^2\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x}.
  \end{aligned}
  \end{equation}
  
  Now we estimate each term in \eqref{weakgeneralizedStokeswithpressureinteriorb}. Thanks to the uniform ellipticity of \eqref{probgenStokesunifellip} with $\delta=\alpha$, proved in Proposition \ref{probgenStokesisunifelliptic}, the term in the left hand side of \eqref{weakgeneralizedStokeswithpressureinteriorb} is bounded below as follows
  \mathbf{e}gin{equation}
    \int_\Omega (\mathsf{D}(\eta_0\mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{D}(\eta_0\mathbf{v}_i)) \cdot \mathbf{n}abla (\eta_0\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x} \geq \alpha \|\mathbf{n}abla(\eta_0\mathbf{v}_i)\|_{L^2(\Omega)}^2.
  \end{equation}
  
  The first term in the right hand side of \eqref{weakgeneralizedStokeswithpressureinteriorb} is bounded using H\"older's inequality and the $L^\infty$ bounds on $A$ and on $\mathbf{n}abla \eta_0$:
  \mathbf{e}gin{multline*}
    - \int_\Omega (\mathsf{S}((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{S}((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i)) \cdot \mathbf{n}abla (\eta_0\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x} \\
    \leq 6\|\mathsf{A}\|_{L^\infty(\Omega)}\|\mathbf{n}abla \eta_0\|_{L^\infty(\Omega)} \|\mathbf{n}abla\mathbf{v}\|_{L^2(\Omega)}\|\mathbf{n}abla (\eta_0\mathbf{v}_i)\|_{L^2(\Omega)}.
  \end{multline*}
  The second term is estimated using again H\"older's inequality and the $L^\infty$ bounds on $A$ and on $\mathbf{n}abla \eta_0$, together with the fact that $0\leq \eta_0 \leq 1$:
  \mathbf{e}gin{multline*}
     \int_\Omega (\mathsf{D}(\eta_0\mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{D}(\eta_0\mathbf{v}_i)) \cdot ((\mathbf{n}abla\eta_0)\otimes (\eta_0\mathbf{v}_i)) \;{\text{\rm d}} \mathbf{x}  \\
    \leq 6\|\mathsf{A}\|_{L^\infty(\Omega)}\|\mathbf{n}abla \eta_0\|_{L^\infty(\Omega)} \|\mathbf{n}abla\mathbf{v}\|_{L^2(\Omega)}\|\mathbf{n}abla (\eta_0\mathbf{v}_i)\|_{L^2(\Omega)}.
  \end{multline*}

  The estimate on the third term also depends on the $L^\infty$ bounds on $A$ and on $\mathbf{n}abla \eta_0$, using again the fact that $0\leq \eta_0 \leq 1$:
  \mathbf{e}gin{multline*}
    - \int_\Omega (\mathsf{S}((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i)\mathsf{A} + \mathsf{A}\mathsf{S}((\mathbf{n}abla\eta_0)\otimes \mathbf{v}_i)) \cdot ((\mathbf{n}abla\eta_0)\otimes (\eta_0\mathbf{v}_i)) \;{\text{\rm d}} \mathbf{x} \leq 6\|\mathsf{A}\|_{L^\infty(\Omega)}\|\mathbf{n}abla \eta_0\|_{L^\infty(\Omega)}^2 \|\mathbf{n}abla\mathbf{v}\|_{L^2(\Omega)}^2.
  \end{multline*}
  
  The fourth term in the right hand side depends on the $L^\infty$ bound on $\partial_{x_i}\mathsf{A}$:
  \mathbf{e}gin{multline*} 
    \int_{\Omega}  (\mathsf{D}(\mathbf{v})(\partial_{x_i}\mathsf{A}) + (\partial_{x_i}\mathsf{A})\mathsf{D}(\mathbf{v}))\cdot \mathbf{n}abla (\eta_0^2\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x}  \\
     \leq 6\|\partial_{x_i}\mathsf{A}\|_{L^\infty(\Omega)}\|\mathbf{n}abla\mathbf{v}\|_{L^2(\Omega)}(\|\mathbf{n}abla\eta_0\|_{L^\infty(\Omega)}\|\mathbf{n}abla\mathbf{v}\|_{L^2(\Omega)} + \|\mathbf{n}abla(\eta_0\mathbf{v}_i)\|_{L^2(\Omega)}).
  \end{multline*}  
    
  The pressure term is estimated using H\"older's inequality and $L^\infty$ bounds on $\mathbf{n}abla \eta_0$ and on $D^2\eta_0$, along with the fact that $0\leq \eta_0 \leq 1$:
  \mathbf{e}gin{multline*} 
  2\int_\Omega  p \partial_{x_i} ((\mathbf{n}abla \eta_0) \cdot (\eta_0\mathbf{v}_i)) \;{\text{\rm d}} \mathbf{x} \leq 2\|p\|_{L^2(\Omega)} \|\mathbf{n}abla((\mathbf{n}abla \eta_0) \cdot (\eta_0\mathbf{v}_i))\|_{L^2(\Omega)} \\
  \leq 2\|p\|_{L^2(\Omega)} \left(\|\mathbf{n}abla (\mathbf{n}abla\eta_0)\|_{L^\infty(\Omega)} \|\eta_0\mathbf{v}_i\|_{L^2(\Omega)} + \|\mathbf{n}abla\eta_0\|_{L^\infty(\Omega)} \|\mathbf{n}abla (\eta_0\mathbf{v}_i)\|_{L^2(\Omega)}\right) \\
  \leq 2\|p\|_{L^2(\Omega)} \left(\|D^2\eta_0\|_{L^\infty(\Omega)} \|\mathbf{n}abla\mathbf{v}\|_{L^2(\Omega)} + \|\mathbf{n}abla\eta_0\|_{L^\infty(\Omega)} \|\mathbf{n}abla (\eta_0\mathbf{v}_i)\|_{L^2(\Omega)}\right).
  \end{multline*}
 
 The estimate on the forcing term in the right hand side uses H\"older's inequality and the $L^\infty$ bound on $\mathbf{n}abla\eta_0$, using again that $0\leq \eta_0\leq 1$:
  \mathbf{e}gin{multline*}
    - \int_\Omega \mathbf{b}f \cdot \partial_{x_i}(\eta_0^2\mathbf{v}_i) \;{\text{\rm d}} \mathbf{x} \leq \|\mathbf{b}f\|_{L^2(\Omega)} \|\mathbf{n}abla(\eta_0^2\mathbf{v}_i)\|_{L^2(\Omega)} \\
    \leq \|\mathbf{b}f\|_{L^2(\Omega)}(\|\mathbf{n}abla\eta_0\|_{L^\infty(\Omega)}\|\eta_0\mathbf{v}_i\|_{L^2(\Omega)} + \|\eta_0\mathbf{n}abla(\eta_0\mathbf{v}_i)\|_{L^2(\Omega)}) \\
    \leq \|\mathbf{b}f\|_{L^2(\Omega)}(\|\mathbf{n}abla\eta_0\|_{L^\infty(\Omega)}\|\mathbf{n}abla\mathbf{v}\|_{L^2({\Omega})} + \|\mathbf{n}abla(\eta_0\mathbf{v}_i)\|_{L^2(\Omega)}).    
  \end{multline*}

  Now, using Young's inequality we put all the estimates together and arrive at 
\mathbf{e}gin{multline}
    \label{interiorestimateh2}
    \frac{\alpha}{2} \|\mathbf{n}abla(\eta_0\mathbf{v}_i)\|_{L^2(\Omega)}^2 \leq 6\|\mathsf{A}\|_{L^\infty(\Omega)}\|\mathbf{n}abla \eta_0\|_{L^\infty(\Omega)}^2 \|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)}^2 + \frac{2}{\alpha} 12^2\|\mathsf{A}\|_{L^\infty(\Omega)}^2\|\mathbf{n}abla \eta_0\|_{L^\infty(\Omega)}^2 \|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)}^2 \\
    + 6\|\partial_{x_i}\mathsf{A}\|_{L^\infty(\Omega)} \|\mathbf{n}abla \eta_0\|_{L^\infty(\Omega)}\|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)}^2 + \frac{2}{\alpha}6^2\|\partial_{x_i}\mathsf{A}\|_{L^\infty(\Omega)}^2 \|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)}^2 \\
    + 2\|p\|_{L^2(\Omega)} \|\mathbf{n}abla(\mathbf{n}abla\eta_0)\|_{L^\infty(\Omega)}\|\mathbf{n}abla\mathbf{v}\|_{L^2(\Omega)} + \frac{2}{\alpha}2^2\|p\|_{L^2(\Omega)}^2 \|\mathbf{n}abla\eta_0\|_{L^\infty(\Omega)}^2 \\
    + \|\mathbf{b}f\|_{L^2(\Omega)}\|\mathbf{n}abla\eta_0\|_{L^\infty(\Omega)}\|\mathbf{n}abla\mathbf{v}\|_{L^2(\Omega)} + \frac{2}{\alpha}  \|\mathbf{b}f\|_{L^2(\Omega)}^2.
\end{multline}
This can be simplified to the following form:
\mathbf{e}gin{multline}
  \label{interiorestimateh2b}
  \frac{\alpha}{2} \|\mathbf{n}abla(\eta_0\mathbf{v}_i)\|_{L^2(\Omega)}^2 \leq \left(\alpha \lambda_1 c_1(\eta_0) + \frac{\lambda_1c_2(\eta_0)}{\alpha}\|\mathsf{A}\|_{L^\infty(\Omega)}^2 + \frac{c}{\alpha}\|\partial_{x_i}\mathsf{A}\|_{L^\infty(\Omega)}^2\right)\|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)}^2 \\
  + \frac{c}{\alpha}\|\mathbf{b}f\|_{L^2(\Omega)}^2 + \frac{\lambda_1 c_1(\eta_0)}{\alpha}\|p\|_{L^2(\Omega)}^2,
\end{multline}
where $c$ is a universal constant and $c_1(\eta_0) = c\lambda_1^{-1}(\|\mathbf{n}abla \eta_0\|_{L^\infty(\Omega)}^2 + \|\mathbf{n}abla(\mathbf{n}abla\eta_0)\|_{L^\infty(\Omega)})$, $c_2(\eta_0) = c\lambda_1^{-1}\|\mathbf{n}abla \eta_0\|_{L^\infty(\Omega)}^2$. Notice that $c_1(\eta_0)$ and $c_2(\eta_0)$ are non-dimensional and can, in fact, be chosen independently of translations, rotations and dilations of $\Omega$, so they are shape constants.

Using \eqref{alphaboundedLinfty}, the inequality \eqref{interiorestimateh2b} yields
\mathbf{e}gin{multline}
  \label{interiorestimateh2c}
  \|\mathbf{n}abla(\eta_0\mathbf{v}_i)\|_{L^2(\Omega)}^2 \\ \leq \frac{c(\Omega)}{\alpha^2}\left( \|\mathbf{b}f\|_{L^2(\Omega)}^2 + \lambda_1\|p\|_{L^2(\Omega)}^2 + \left(\lambda_1\|\mathsf{A}\|_{L^\infty(\Omega)}^2 + \|\mathbf{n}abla\mathsf{A}\|_{L^\infty(\Omega)}^2\right)\|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)}^2\right),
\end{multline}
for $i=1,2,3$, where $c(\Omega)$ is a shape constant.

  As for the estimates near the boundary, they are obtained, as mentioned above, by first rectifying each portion of the covering of the boundary. For each $m=1,\ldots, M$, this amounts to a change of variables $\mathbf{v}arphi: B \rightarrow {\mathbb{R}}^3$ that is a diffeomorphism of class ${\mathcal C}^2$ from a ball $B\subset{\mathbb{R}}^3$ onto the open set $U_m$. The transformation $\mathbf{x} = \mathbf{v}arphi(\mathbf{y})$ changes the coefficients of $\mathsf{A}$ and of the forcing term $\mathbf{b}f$ to terms with the same regularity, and bounded by the $W^{1,\infty}(\Omega)^{3\times 3}$ norm of $\mathsf{A}$ and the $L^2(\Omega)^3$ norm of $\mathbf{b}f$, respectively (see \cite[Theorem 13.49]{AB}). Thus,  the weak formulation changes to 
\mathbf{e}gin{multline}
    \label{weakgeneralizedStokeswithpressureinhemisphere}
    \int_{B^-} (\mathsf{D}(\mathbf{v}(\mathbf{v}arphi(\mathbf{y})))\mathsf{A}(\mathbf{v}arphi(\mathbf{y})) + \mathsf{A}(\mathbf{v}arphi(\mathbf{y}))\mathsf{D}(\mathbf{v}(\mathbf{v}arphi(\mathbf{y})))) : \mathbf{n}abla \mathbf{w}(\mathbf{v}arphi(\mathbf{y})) \det D\mathbf{v}arphi(\mathbf{y})\;{\text{\rm d}} \mathbf{y} \\ 
    + \int_{B^-} p(\mathbf{v}arphi(\mathbf{y})) (\mathbf{n}abla \cdot \mathbf{w})(\mathbf{v}arphi(\mathbf{y})) \det D\mathbf{v}arphi(\mathbf{y})\;{\text{\rm d}} \mathbf{y} = \int_{B^-} \mathbf{b}f(\mathbf{v}arphi(\mathbf{y})) \cdot \mathbf{w}(\mathbf{v}arphi(\mathbf{y})) \det D\mathbf{v}arphi(\mathbf{y})\;{\text{\rm d}} \mathbf{y},
\end{multline}  
for all $\mathbf{w}\in H_0^1(\Omega)^3$, where $D\mathbf{v}arphi$ is the differential of $\mathbf{v}arphi$ and $B^-= B \cap \{y_3<0\} = \mathbf{v}arphi^{-1}(\Omega\cap U_m)$ is an open hemisphere (the ``southern'' hemisphere) with flat boundary $B\cap \{y_3=0\}$. 

We introduce the rectified variables 
\mathbf{e}gin{equation}
  \label{localizedrectifiedvariables}
  \tilde\mathbf{v}(\mathbf{y}) = \mathbf{v}(\mathbf{v}arphi(\mathbf{y})), \qquad \tilde p(\mathbf{y}) = p(\mathbf{v}arphi(\mathbf{y})),
\end{equation}
and use that
\[ \mathbf{n}abla \tilde\mathbf{v}(\mathbf{y}) = \mathbf{n}abla\mathbf{v}(\mathbf{v}arphi(\mathbf{y}))D\mathbf{v}arphi(\mathbf{y}),
\]
and
\[ \mathsf{D}(\mathbf{v}(\mathbf{v}arphi(\mathbf{y})) = \mathsf{S}(\mathbf{n}abla\tilde\mathbf{v}(\mathbf{y}) D\mathbf{v}arphi(\mathbf{y})^{-1}),
\]
where $\mathsf{S}$ is the symmetrization operator defined in \eqref{defopS}. Similarly for $\tilde p$, $\tilde \mathbf{w}$, with
\[ (\mathbf{n}abla \cdot \mathbf{w})(\mathbf{v}arphi(\mathbf{y})) = D\tilde\mathbf{w}(\mathbf{y}) : D\mathbf{v}arphi(\mathbf{y})^{-tr},
\]
where $D\mathbf{v}arphi(\mathbf{y})^{-tr} = (D\mathbf{v}arphi(\mathbf{y})^{-1})^\text{\rm tr}$. We also define $\tilde \mathsf{A}(\mathbf{y}) = \mathsf{A}(\mathbf{v}arphi(\mathbf{y}))$ and $\tilde\mathbf{b}f(\mathbf{y}) = \mathbf{b}f(\mathbf{v}arphi(\mathbf{y}))$. Then, the weak formulation is rewritten as
\mathbf{e}gin{multline}
    \label{weakgeneralizedStokeswithpressureinhemispherenewvar}
    \int_{B^-} (\mathsf{S}(\mathbf{n}abla\tilde\mathbf{v}(\mathbf{y}) D\mathbf{v}arphi(\mathbf{y})^{-1})\tilde \mathsf{A}(\mathbf{y}) + \tilde \mathsf{A}(\mathbf{y})\mathsf{S}(\mathbf{n}abla\tilde\mathbf{v}(\mathbf{y}) D\mathbf{v}arphi(\mathbf{y})^{-1}) : (\mathbf{n}abla \tilde\mathbf{w}(\mathbf{y}) D\mathbf{v}arphi(\mathbf{y})^{-1}) \det D\mathbf{v}arphi(\mathbf{y})\;{\text{\rm d}} \mathbf{y} \\ 
    + \int_{B^-} \tilde p(\mathbf{y}) (D\tilde\mathbf{w}(\mathbf{y}) : D\mathbf{v}arphi(\mathbf{y})^{-tr}) \det D\mathbf{v}arphi(\mathbf{y}))\;{\text{\rm d}} \mathbf{y} = \int_{B^-}\tilde\mathbf{b}f(\mathbf{y}) \cdot \tilde\mathbf{w}(\mathbf{y}) \det D\mathbf{v}arphi(\mathbf{y})\;{\text{\rm d}} \mathbf{y},
\end{multline}  
for any ``rectified'' test function $\tilde\mathbf{w}\in H_0^1(B^-)^3$.  The divergence-free condition on $\mathbf{v}$ is transformed to the condition
\[ D\tilde\mathbf{v}(\mathbf{y}) : D\mathbf{v}arphi(\mathbf{y})^{-tr} = 0
\]
on $\tilde\mathbf{v}$.

The ellipticity of the higher order term with respect to the vector field $\mathbf{n}abla\tilde\mathbf{v} D\mathbf{v}arphi^{-1}$ in the problem \eqref{weakgeneralizedStokeswithpressureinhemispherenewvar} follows immediately from the estimate \eqref{symposopL} in Lemma \ref{lempositivityofL} and the fact that $\mathbf{v}arphi$ is a diffeomorphism. Moreover, since $\mathbf{v}arphi$ is of class ${\mathcal C}^2$, all the coefficients have the appropriate regularity to yield an estimate for the tangential derivatives $\partial_{y_1}\tilde\mathbf{v}$ and $\partial_{y_2}\tilde\mathbf{v}$ of the same form as the interior estimates \eqref{interiorestimateh2c}. Here, we exploit the fact that, since $\tilde\mathbf{v}$ vanishes on the boundary and the boundary is flat, all the derivatives on the directions tangential to the boundary are well defined and vanish on the boundary as well, so that they can be used as test functions. In particular, this implies that the tangential derivatives also vanish on the boundary, which is expected. The estimates rely on the inequalities \eqref{finitedifferenceintestimate2} and \eqref{finitedifferencebdryestimateinfty} and are of the same form as \eqref{interiorestimateh2c}.

Finally, using \eqref{weakgeneralizedStokeswithpressureinhemispherenewvar}, we isolate the distributional derivatives $\partial_{y_3}\tilde\mathbf{v}$, in the direction normal to the boundary, in terms of the derivatives $\partial_{y_1}\tilde\mathbf{v}$ and $\partial_{y_2}\tilde\mathbf{v}$ and also lower order terms, to obtain an estimate of the form \eqref{interiorestimateh2c} for the derivative in the normal direction too. 

With all the local estimates ready, we use the partition of unity to write
\[D^2 \mathbf{v} = \mathbf{n}abla(\mathbf{n}abla\mathbf{v}) = \mathbf{n}abla\left(\sum_{m=0}^M \eta_m \mathbf{n}abla\mathbf{v}\right) = \sum_{m=0}^M \mathbf{n}abla(\eta_m \mathbf{n}abla\mathbf{v}),
\]
so that, putting all the estimates together, we arrive at the following global estimate for $D^2\mathbf{v}$:
\mathbf{e}gin{equation}
  \label{estimateh2}
  \|D^2 \mathbf{v}\|_{L^2(\Omega)}^2 \leq \frac{c(\Omega)}{\alpha^2}\left( \|\mathbf{b}f\|_{L^2(\Omega)}^2 + \lambda_1\|p\|_{L^2(\Omega)}^2 +\left(\lambda_1\|\mathsf{A}\|_{L^\infty(\Omega)}^2 + \|\mathbf{n}abla\mathsf{A}\|_{L^\infty(\Omega)}^2\right)\|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)}^2\right),
\end{equation}
where $c(\Omega)$ is a shape constant.

Using the estimates \eqref{weaksolutionweakgeneralizedStokesestimaterepeat} and  \eqref{weaksolutionweakgeneralizedStokesestimatewithpressure} for the derivatives of the velocity field and the pressure, we obtain precisely  \eqref{h2regularityforh2v}, for another shape constant $c(\Omega)$.

An estimate for $\mathbf{n}abla p$ is obtained directly from the equation \eqref{generalizedStokes}: 
\mathbf{e}gin{equation}
  \label{basicpestimatefromStokes}
  \|\mathbf{n}abla p \|_{L^2(\Omega)}  \leq \|\mathbf{b}f\|_{L^2(\Omega)} + \|\mathsf{A}\|_{L^\infty(\Omega)}\|D^2\mathbf{v}\|_{L^2(\Omega)} + \|\mathbf{n}abla\mathsf{A}\|_{L^\infty(\Omega)}\|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)}.
\end{equation}  
Using the estimates \eqref{weaksolutionweakgeneralizedStokesestimaterepeat} and \eqref{h2regularityforh2v}, we obtain
\mathbf{e}gin{align*}
  \|\mathbf{n}abla p \|_{L^2(\Omega)} & \leq \|\mathbf{b}f\|_{L^2(\Omega)}  + \frac{1}{\alpha}\|\mathbf{n}abla\mathsf{A}\|_{L^\infty(\Omega)}\|\mathbf{b}f\|_{H^{-1}(\Omega)}  \\
    & \qquad + \frac{c(\Omega)}{\alpha}\|\mathsf{A}\|_{L^\infty(\Omega)}\left( \|\mathbf{b}f\|_{L^2(\Omega)} + \frac{1}{\alpha}\left(\lambda_1^{1/2}\|\mathsf{A}\|_{L^\infty(\Omega)} + \|\mathbf{n}abla\mathsf{A}\|_{L^\infty(\Omega)}\right)\|\mathbf{b}f\|_{H^{-1}(\Omega)}\right),
\end{align*}  
Then, using \eqref{alphaboundedLinfty} and replacing $c(\Omega)$ by $c(\Omega)+1$, we find the estimate \eqref{h2regularityforh1p}, which completes the proof.
\end{proof}

\mathbf{e}gin{rmk}
  The proof of regularity in Theorem \ref{h2regularity} may be done directly, without resorting to \cite[Theorem 4.2 with $k=0$]{macha2011}, by working with the finite differences defined in \eqref{finitedifference}. Instead of taking the differential of the equation \eqref{generalizedStokes} to arrive at \eqref{generalizedStokesforderivative} and test the localized versions of this equation against $\eta_0^2\mathbf{v}_i$, where $\mathbf{v}_i = \partial_{x_i}\mathbf{v}$, we test the localized versions of \eqref{generalizedStokes} against $\delta_i^{-h}(\eta_0^2\delta_i^h\mathbf{v})$ and use a change of variables that works like integration by parts. This leads to estimates like \eqref{interiorestimateh2c}, but with $\delta_i^h\mathbf{v}$ at the left hand side in the place of $\mathbf{v}_i$. Then, passing to the limit as $h\rightarrow 0$, we obtain the regularity and the corresponding $H^1$-bound for $\partial_{x_i}\mathbf{v}$, for every $i=1,2,3$. This process is done in more details for the higher-order derivatives obtained in the following section, for which we did not find any regularity result with the desired smoothness conditions on the coefficients of the system.
\end{rmk}

\mathbf{e}gin{rmk}
In the derivation of the localized and rectified weak formulation \eqref{weakgeneralizedStokeswithpressureinhemispherenewvar}, since $\mathbf{v}$ is to be interpreted as a velocity field, say corresponding to the speed of a trajectory $\mathbf{x}(t)=\mathbf{v}arphi(\mathbf{y}(t))$, it would be more natural do defined the velocity $\tilde\mathbf{v}(\mathbf{y})$ following the idea that 
\[ \mathbf{v}(t) = \frac{{\text{\rm d}} \mathbf{x}(t)}{{\text{\rm d}} t} = D\mathbf{v}arphi(\mathbf{y}(t))\frac{{\text{\rm d}} \mathbf{y}(t)}{{\text{\rm d}} t} = D\mathbf{v}arphi(\mathbf{y}(t))\tilde\mathbf{v}(t).
\]
This would lead to the rectified velocity field $\tilde\mathbf{v}(\mathbf{y}) = (D\mathbf{v}arphi(\mathbf{y}))^{-1}\mathbf{v}(\mathbf{v}arphi(\mathbf{y}))$, 
and the localized and rectified weak formulation 
\mathbf{e}gin{multline}
    \label{weakgeneralizedStokeswithpressureinhemispherenewvarold}
    \int_{B^-} (D\mathbf{v}arphi(\mathbf{y})\mathsf{D}_{\mathbf{y}}(\tilde\mathbf{v}(\mathbf{y}))D\mathbf{v}arphi(\mathbf{y})^{-1}\mathsf{A}(\mathbf{v}arphi(\mathbf{y})) + \mathsf{A}(\mathbf{v}arphi(\mathbf{y}))D\mathbf{v}arphi(\mathbf{y})\mathsf{D}_{\mathbf{y}}(\tilde\mathbf{v}(\mathbf{y}))D\mathbf{v}arphi(\mathbf{y})^{-1}) \\
      : D\mathbf{v}arphi(\mathbf{y})\mathbf{n}abla_{\mathbf{y}} \tilde\mathbf{w}(\mathbf{y})D\mathbf{v}arphi(\mathbf{y})^{-1} \det D\mathbf{v}arphi(\mathbf{y})\;{\text{\rm d}} \mathbf{y} \\ 
    + \int_{B^-} \tilde p(\mathbf{y}) \text{\rm tr}(D\mathbf{v}arphi(\mathbf{y})\mathbf{n}abla_{\mathbf{y}} \tilde\mathbf{w}(\mathbf{y})D\mathbf{v}arphi(\mathbf{v})^{-1}) \det D\mathbf{v}arphi(\mathbf{y})\;{\text{\rm d}} \mathbf{y} = \int_{B^-} \tilde\mathbf{b}f(\mathbf{y}) \cdot \tilde\mathbf{w}(\mathbf{y}) \det D\mathbf{v}arphi(\mathbf{y})\;{\text{\rm d}} \mathbf{y},
\end{multline}
This works perfectly fine, except that we would need one more degree of regularity from the boundary. For the velocity field $\mathbf{v}$ to be in $H^2$, we would need $\mathbf{v}arphi$ of class ${\mathcal C}^3$. In order to avoid this unnecessary extra regularity, we perform the change of variables according \eqref{localizedrectifiedvariables}, as done is most classical works (see e.g. \cite{adn2, giaquintamodica1982, huystara2006, macha2011}).
\end{rmk}

\subsection{Higher-order regularity of the solutions of the generalized Stokes problem}
\label{subsechighreggenStokes}

For higher-order regularity, the result of \cite[Theorem 4.2 with $k>0$]{macha2011} is not optimal, requiring too much regularity from $\mathsf{A}$ than we can afford for the well-posedness of the associated evolutionary problem \eqref{LiueulerB}. (The result \cite[Theorem 4.2 with $k>0$]{macha2011} requires $\mathsf{A}\in W^{k+1,\infty}(\Omega)$, besides $\mathbf{b}f \in H^k(\Omega)^3$, to have $\mathbf{v}\in H^{k+2}(\Omega)^3$ and $p\in H^{k+1}(\Omega)$). Similarly in the works \cite{giaquintamodica1982, huystara2006}.) 

The fundamental work of Agmon, Douglis and Niremberg \cite{adn2} also requires too much regularity from $\mathsf{A}$ (see the condition just before Theorem 10.5, in page 78, of \cite{adn2}, in which it is required that $\mathsf{A}\in {\mathcal C}^{l-s_i}(\mathbf{a}r\Omega)$, with $l$ equals to our $k$, and where $s_i$ can be chosen, in our case, to be $s_1=s_2=s_3=0$ and $s_4=-1$, so that we need at least some coefficients of $\mathsf{A}$ to be ${\mathcal C}^{k+1}$, and others to be ${\mathcal C}^k$, to obtain the desired higher-order regularity for $\mathbf{v}$ and $p$). 

It is worth mentioning the work of Solonnikov \cite{solonnikov2001}, which treats a general Stokes-type problem that includes ours, but only addresses H\"older-continuity, and the work of Ghidaglia \cite{ghidaglia1984}, which has a more direct approach to the classical Stokes problem, but does not apply to our system. 

For general and similar elliptic systems, but not of Stokes type, the works of Friedman \cite{friedman1969}, Lions and Magenes \cite{lionsmagenes1972} and Gilbarg and Trudinger \cite{gilbargtrudinger2001} have results that also assume more regularity from the coefficients of the linear term. 

Most of the works omit the proof of higher-order regularity, limiting themselves to say that the proof follows as in the $H^2$ regularity. The work \cite{friedman1969}, however, does give some details of the proof, and the reason for the regularity conditions becomes clear. The main reason for the high-order regularity assumption on coefficients of the linear term is that the theorems do not take into full consideration the space dimension of the system and the associated Sobolev embeddings. In our case, we exploit the fact that the space dimension is three in order to obtain a more forgiven condition on $\mathsf{A}$. See Remark \ref{dontneedsomuchregularityonA} for more discussion on this issue. With that in mind, we present, in Theorem \ref{hkplustworegularity} below, a result assuming less regularity on the coefficients than the previously mentioned works, and which will be perfect for our needs. 

\mathbf{e}gin{thm}[Higher-order regularity]
  \label{hkplustworegularity}
    Let $\Omega$ be a bounded domain in ${\mathbb{R}}^3$ with a boundary of class ${\mathcal C}^{k+2}$ and suppose $\mathsf{A}\in  W^{1,\infty}(\Omega)^{3\times 3} \cap W^{2,3}(\Omega)^{3\times 3} \cap H^{k+1}(\Omega)^{3\times 3}$ is a uniformly positive definite symmetric tensor on $\Omega$ and $\mathbf{b}f \in H^k(\Omega)^3$, where $k\in {\mathbb{N}}$. Then, the solution $(\mathbf{v}, p)$ of the weak formulation \eqref{weakgeneralizedStokeswithpressure} of the generalized Stokes problem is such that $\mathbf{v} \in V(\Omega)\cap H^{k+2}(\Omega)$ and $p\in H^{k+1}(\Omega)$. Moreover, $(\mathbf{v},p)$ solves the generalized Stokes problem \eqref{generalizedStokes} almost everywhere in $\Omega$, with the first two equations in \eqref{generalizedStokes} holding in $H^{k-1}(\Omega)^3$ and $H^k(\Omega)$, respectively, and with the last equation holding for the trace of $\mathbf{v}$ on $\partial\Omega$, in the space $H^{k+1/2}(\partial \Omega)$. Finally, for $k=1$, we have the estimates
   \mathbf{e}gin{multline} 
  \label{h3egularityforv}
  \|D^3 \mathbf{v}\|_{L^2(\Omega)} \leq \frac{c_3(\Omega)}{\alpha} \left(\|\mathbf{b}f\|_{H^1_{\lambda_1}(\Omega)}  \right. \\
    \left.  + \frac{1}{\alpha}\left(\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)} + \|D^2\mathsf{A}\|_{L^3(\Omega)}\right)\left( \|\mathbf{b}f\|_{L^2(\Omega)} + \frac{1}{\alpha}\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}\|\mathbf{b}f\|_{H^{-1}(\Omega)}\right)\right),
   \end{multline}    
    and
    \mathbf{e}gin{multline}
      \label{h2egularityforp}
      \|D^2p\|_{L^2(\Omega)} \leq \frac{c_3(\Omega)}{\alpha}\|\mathsf{A}\|_{L^\infty(\Omega)} \left(\|\mathbf{b}f\|_{H^1_{\lambda_1}(\Omega)}  \right. \\
     \left. + \frac{1}{\alpha}\left(\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)} + \|D^2\mathsf{A}\|_{L^3(\Omega)}\right)\left( \|\mathbf{b}f\|_{L^2(\Omega)} + \frac{1}{\alpha}\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}\|\mathbf{b}f\|_{H^{-1}(\Omega)}\right) \right),
    \end{multline}     
 while, for any integer $k\geq 2$,
 \mathbf{e}gin{equation} 
  \label{hkplus2regularityforv}
  \|D^{k+2} \mathbf{v}\|_{L^2(\Omega)} \leq \frac{c_{k+2}(\Omega)}{\alpha}R_k(\alpha, \lambda_1, \mathsf{A}, \mathbf{b}f),
\end{equation}    
    and
    \mathbf{e}gin{equation}
      \label{hkplus1regularityforp}
      \|D^{k+1}p\|_{L^2(\Omega)} \leq \frac{c_{k+2}(\Omega)}{\alpha}\|\mathsf{A}\|_{L^\infty(\Omega)}R_k(\alpha, \lambda_1, \mathsf{A}, \mathbf{b}f),
    \end{equation}
    where 
    \mathbf{e}gin{equation}
      \label{Rkdef}
      \mathbf{e}gin{aligned}
      R_k(\alpha, & \lambda_1, \mathsf{A}, \mathbf{b}f) = \|\mathbf{b}f\|_{H^k_{\lambda_1}(\Omega)} + \sum_{j=2}^{k-1} \left( \sum_{i=1}^{k-j} \frac{1}{\alpha^{(k-j)/i}}\frac{1}{\lambda_1^{(k-j)/4i}}\|\mathsf{A}\|_{H^{i+2}_{\lambda_1}(\Omega)}^{(k-j)/i}\right)  \|\mathbf{b}f\|_{H^j_{\lambda_1}(\Omega)}  \\
   &  + \left( \sum_{i=1}^{k-1} \frac{1}{\alpha^{(k-1)/i}}\frac{1}{\lambda_1^{(k-1)/4i}}\|\mathsf{A}\|_{H^{i+2}_{\lambda_1}(\Omega)}^{(k-1)/i}\right)\left( \|\mathbf{b}f\|_{H^1_{\lambda_1}(\Omega)} \right. \\
  & \quad \quad + \frac{1}{\alpha}\left( \|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)} + \|D^2\mathsf{A}\|_{L^3(\Omega)}\right)\left(\|\mathbf{b}f\|_{L^2(\Omega)} + \frac{1}{\alpha}\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)} \|\mathbf{b}f\|_{H^{-1}(\Omega)}\right),
    \end{aligned}
    \end{equation}
    and the $c_{k+2}(\Omega)$, for $k\geq 1$, are non-dimensional shape constants.
\end{thm}

\mathbf{e}gin{proof}
  The proof is obtained in a way similar to the proof of the estimates obtained in Theorem \ref{h2regularity} for the regularity of $(\mathbf{v},p)$ in $H^2(\Omega)\times H^1(\Omega)$. The difference being that, since we do not have yet the desired regularity, we are not allowed to take the full higher-order derivative in the equation \eqref{generalizedStokes}. Instead we use the classical technique of \cite{adn2} of estimating the finite differences, which are then showed to converge to the desired higher-order derivative, obtaining the desired bounds along the way.
  
  We argue by induction that, for a given $k\in{\mathbb{N}}_0$, the regularity has been obtained up to $(\mathbf{v},p)\in (V(\Omega)\cap H^{k+1}(\Omega)^3) \times H^k(\Omega)$, and we intend to prove that it holds in the Sobolev spaces of order $k+2$ and $k+1$, respectively. The case $k=0$ has been proved in Theorem \ref{h2regularity}, so we continue here with the induction process for $k\in {\mathbb{N}}$. We need to show that $\mathbf{v}\in H^{k+2}$ and $p\in H^{k+1}$.
  
  By taking the derivative of \eqref{generalizedStokes} up to order $k$, we find the equations
\mathbf{e}gin{equation}
  \label{generalizedStokesforderivativegamma}
   \mathbf{e}gin{cases}
      - \mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v}_\mathbf{g}amma)\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v}_\mathbf{g}amma)) - \mathbf{b}_\mathbf{g}amma + \mathbf{n}abla p_\mathbf{g}amma =  \mathbf{b}f_\mathbf{g}amma, & \text{in } \Omega,\\
      \mathbf{n}abla \cdot \mathbf{v}_\mathbf{g}amma = 0, & \text{in } \Omega,
    \end{cases}
\end{equation}
for any multi-index $\mathbf{g}amma\in {\mathbb{N}}_0^3$ with $|\mathbf{g}amma|=k$, where 
\[ \mathbf{v}_\mathbf{g}amma = D^\mathbf{g}amma\mathbf{v}, \qquad p_\mathbf{g}amma = D^\mathbf{g}amma p, \qquad \mathbf{b}f_\mathbf{g}amma = D^\mathbf{g}amma \mathbf{b}f,
\]
and, by recursion,
\mathbf{e}gin{equation}
  \label{bgammadef}
  \mathbf{b}_\mathbf{g}amma = \partial_{x_j}\mathbf{b}_{\mathbf{g}amma'} + \mathbf{B}(\partial_{x_j}\mathsf{A},\mathbf{v}_{\mathbf{g}amma'}) = \partial_{x_j}\mathbf{b}_{\mathbf{g}amma'} + \mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v}_{\mathbf{g}amma'})(\partial_{x_j}\mathsf{A}) + (\partial_{x_j}\mathsf{A})\mathsf{D}(\mathbf{v}_{\mathbf{g}amma'})),
\end{equation}
where $\mathbf{g}amma'\in {\mathbb{N}}_0^3$ is another multi-index, one-degree lower than $\mathbf{g}amma$, given by $\mathbf{g}amma = \mathbf{g}amma' + \mathbf{e}_j$, with the convention that $\mathbf{b}_{\mathbf{g}amma'} = 0$, if $|\mathbf{g}amma'|=0$. Since we are assuming that $(\mathbf{v},p)\in V(\Omega)\cap H^{k+1}(\Omega)^3 \times H^k(\Omega)$, the first equation in \eqref{generalizedStokesforderivativegamma} holds in $H^{-1}(\Omega)^3$ and the divergence-free equation in \eqref{generalizedStokesforderivativegamma} holds in $L^2(\Omega)$.
  
If we were allowed to take one more derivative of the equation \eqref{generalizedStokesforderivativegamma}, we would see that
\mathbf{e}gin{equation}
  \label{generalizedStokesforhighorderderivative}
   \mathbf{e}gin{cases}
      - \mathbf{n}abla \cdot (\mathsf{D}(\partial_{x_i}\mathbf{v}_\mathbf{g}amma)\mathsf{A} + \mathsf{A}\mathsf{D}(\partial_{x_i}\mathbf{v}_\mathbf{g}amma)) - \mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v}_\mathbf{g}amma)(\partial_{x_i}\mathsf{A}) + (\partial_{x_i}\mathsf{A})\mathsf{D}(\mathbf{v}_\mathbf{g}amma)) - \partial_{x_i} \mathbf{b}_\mathbf{g}amma \\ 
      \qquad \qquad \qquad + \mathbf{n}abla (\partial_{x_i}p_\mathbf{g}amma)  = \partial_{x_i}\mathbf{b}f_\mathbf{g}amma, & \text{in } \Omega, \\
      \mathbf{n}abla \cdot (\partial_{x_i}\mathbf{v}_\mathbf{g}amma) = 0, & \text{in } \Omega,
    \end{cases}
\end{equation}  
for any $i=1,2,3$. The equation \eqref{generalizedStokesforhighorderderivative} has the same form of the equation \eqref{generalizedStokesforderivative} used for the $H^2$ estimate of $\mathbf{v}$, if viewing the additional term $\partial_{x_i}\mathbf{b}_\mathbf{g}amma$ as an extra forcing term. Naturally, the $H^{k+2}$ estimate would follow a similar argument, by localization of the equation on the interior and on patches of the boundary and rectification of the boundary. The estimates in the interior for the second-order derivatives of $D^\mathbf{g}amma\mathbf{v}$ would be essentially the same as those for the second-order derivatives of $\mathbf{v}$, by testing the localized equation with $\eta_0\partial_{x_i}\mathbf{v}_\mathbf{g}amma$. The estimates on the boundary would change slightly, because $\mathbf{v}_\mathbf{g}amma$ does not vanish on the boundary. Nevertheless, any combination of tangential derivatives vanishes on the boundary and that would be sufficient. 

However, we do not have yet sufficient regularity to derive \eqref{generalizedStokesforhighorderderivative}, so we must also prove this regularity.  The idea is to work with \eqref{generalizedStokesforderivativegamma} and use the finite-difference operators defined in \eqref{finitedifference} to obtain $H^1$-estimates for $\eta_m\delta_i^h\mathbf{v}_\mathbf{g}amma$, uniformly with respect to $h$, and then, at the limit as $h\rightarrow 0$, arrive at the regularity and an estimate for the partial derivative $\eta_m\partial_{x_i}\mathbf{v}_\mathbf{g}amma$ of $\mathbf{v}_\mathbf{g}amma$ and for $\eta_m\partial_{x_i}p_\mathbf{g}amma$. This is done for all the derivatives in the interior and for all the combinations of tangential derivatives near the boundary. The derivatives near the boundary that contain normal directions are obtained directly from the equation, by induction on the degree of normal derivatives (see below). Then, by combining all the partitions, we obtain a global estimate for $\mathbf{n}abla\mathbf{v}_\mathbf{g}amma$ and $\mathbf{n}abla p_\mathbf{g}amma$, proving, along the way, that those derivatives have indeed the desired regularity.

More precisely, we consider a covering of $\Omega$ by open sets $\{U_m\}_{m=0}^M$, with $U_0$ in the interior of $\Omega$ and the remaining open sets covering the boundary, as done in the proof of Theorem \ref{h2regularity}, along with a partition of unity by functions $\{\eta_m\}_{m=0}^M$, subordinated to this covering. 

For the interior regularity, we test the weak formulation \eqref{generalizedStokesforderivativegamma} with the function $\mathbf{w} = \delta_i^{-h}(\eta_0^2 \delta_i^h\mathbf{v}_\mathbf{g}amma)$, which does vanish on the boundary, and in fact it vanishes outside $U_0$, so that, after integrating by parts the leading order term and the pressure term, we find
  \mathbf{e}gin{multline}
    \label{weakgeneralizedStokeswithpressureinteriorgammadifference} 
    \int_{U_0} (\mathsf{D}(\mathbf{v}_\mathbf{g}amma)\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v}_\mathbf{g}amma)) : \mathbf{n}abla (\delta_i^{-h}(\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma)) \;{\text{\rm d}} \mathbf{x} - \int_{U_0}  \mathbf{b}_\mathbf{g}amma \cdot (\delta_i^{-h}(\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma)) \;{\text{\rm d}} \mathbf{x}\\
    - \int_{U_0}  p_\mathbf{g}amma \mathbf{n}abla \cdot (\delta_i^{-h}(\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma)) \;{\text{\rm d}} \mathbf{x} = \int_{U_0} \mathbf{b}f_\mathbf{g}amma \cdot (\delta_i^{-h}(\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma)) \;{\text{\rm d}} \mathbf{x}.
  \end{multline}
By changing variables we rewrite the first term in the form
\[ \int_{U_0} (\mathsf{D}(\mathbf{v}_\mathbf{g}amma)\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v}_\mathbf{g}amma)) : \mathbf{n}abla (\delta_i^{-h}(\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma)) \;{\text{\rm d}} \mathbf{x} = \int_{U_0} \delta_i^h(\mathsf{D}(\mathbf{v}_\mathbf{g}amma)\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v}_\mathbf{g}amma)) : \mathbf{n}abla (\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma) \;{\text{\rm d}} \mathbf{x}.
\]

The product rule works for difference operators just like it does for derivatives:
  \mathbf{e}gin{equation}
    \label{finitedifferenceonlinearpart}
    \delta_i^h(\mathsf{D}(\mathbf{v}_\mathbf{g}amma)\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v}_\mathbf{g}amma))  = \mathsf{D}(\delta_i^h\mathbf{v}_\mathbf{g}amma)\mathsf{A} +\mathsf{A}\mathsf{D}(\delta_i^h\mathbf{v}_\mathbf{g}amma) + \mathsf{D}(\mathbf{v}_\mathbf{g}amma)(\delta_i^h\mathsf{A}) + (\delta_i^h\mathsf{A})\mathsf{D}(\mathbf{v}_\mathbf{g}amma).
  \end{equation}
Moreover, we can write
\mathbf{e}gin{equation}  
    \label{breakingetasquare} 
    \mathbf{n}abla(\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma) = \eta_0(\mathbf{n}abla\eta_0)\otimes (\delta_i^h\mathbf{v}_\mathbf{g}amma) + \eta_0 \mathbf{n}abla(\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma),
\end{equation}
and
\mathbf{e}gin{equation}
    \label{breakingDetadelta}
    \mathsf{D}(\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma) = \frac{1}{2}(\mathbf{n}abla (\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma) + (\mathbf{n}abla (\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma))^\text{\rm tr}) = \mathsf{S}((\mathbf{n}abla\eta_0) \otimes (\delta_i^h\mathbf{v}_\mathbf{g}amma)) + \eta_0\mathsf{D}(\delta_i^h\mathbf{v}_\mathbf{g}amma),
  \end{equation}
  where $\mathsf{S}(\cdot)$ is the symmetrization operator defined in \eqref{defopS}. Thus, we rewrite the leading order term as
  \mathbf{e}gin{equation}
    \label{lineartermdeltaihvgamma}
 \mathbf{e}gin{aligned}
     \int_{U_0} & \delta_i^h(\mathsf{D}(\mathbf{v}_\mathbf{g}amma)\mathsf{A} + \mathsf{A}\mathsf{D}(\mathbf{v}_\mathbf{g}amma)) : \mathbf{n}abla (\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma) \;{\text{\rm d}} \mathbf{x} \\
       & = \int_{U_0} (\mathsf{D}(\delta_i^h\mathbf{v}_\mathbf{g}amma)\mathsf{A} +\mathsf{A}\mathsf{D}(\delta_i^h\mathbf{v}_\mathbf{g}amma)) : \mathbf{n}abla (\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma) \;{\text{\rm d}} \mathbf{x} \\
       & \qquad + \int_{U_0} (\mathsf{D}(\mathbf{v}_\mathbf{g}amma)(\delta_i^h\mathsf{A}) + (\delta_i^h\mathsf{A})\mathsf{D}(\mathbf{v}_\mathbf{g}amma)) : \mathbf{n}abla (\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma) \;{\text{\rm d}} \mathbf{x} \\
  & = \int_{U_0} \eta_0(\mathsf{D}(\delta_i^h\mathbf{v}_\mathbf{g}amma)\mathsf{A} +\mathsf{A}\mathsf{D}(\delta_i^h\mathbf{v}_\mathbf{g}amma)) : \left((\mathbf{n}abla\eta_0)\otimes (\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma) + \mathbf{n}abla (\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma) \right)\;{\text{\rm d}} \mathbf{x} \\
   & \qquad    + \int_{U_0} (\mathsf{D}(\mathbf{v}_\mathbf{g}amma)(\delta_i^h\mathsf{A}) + (\delta_i^h\mathsf{A})\mathsf{D}(\mathbf{v}_\mathbf{g}amma)) : \mathbf{n}abla (\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma) \;{\text{\rm d}} \mathbf{x} \\
  & = \int_{U_0} (\mathsf{D}(\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma)\mathsf{A} +\mathsf{A}\mathsf{D}(\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma)) : \left((\mathbf{n}abla\eta_0)\otimes (\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma) + \mathbf{n}abla (\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma) \right)\;{\text{\rm d}} \mathbf{x} \\
   & \qquad   - \int_{U_0} (\mathsf{S}((\mathbf{n}abla\eta_0) \otimes (\delta_i^h\mathbf{v}_\mathbf{g}amma))\mathsf{A} +\mathsf{A}\mathsf{S}((\mathbf{n}abla\eta_0) \otimes (\delta_i^h\mathbf{v}_\mathbf{g}amma))) : \left((\mathbf{n}abla\eta_0)\otimes (\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma) + \mathbf{n}abla (\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma) \right)\;{\text{\rm d}} \mathbf{x} \\
   & \qquad  + \int_{U_0} (\mathsf{D}(\mathbf{v}_\mathbf{g}amma)(\delta_i^h\mathsf{A}) + (\delta_i^h\mathsf{A})\mathsf{D}(\mathbf{v}_\mathbf{g}amma)) : \mathbf{n}abla (\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma) \;{\text{\rm d}} \mathbf{x}.
\end{aligned}
\end{equation}  

Since $\mathbf{v}_\mathbf{g}amma$ is divergence free, we rewrite the pressure term in \eqref{weakgeneralizedStokeswithpressureinteriorgammadifference} as
\[ -\int_{U_0}  p_\mathbf{g}amma \mathbf{n}abla \cdot (\delta_i^{-h}(\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma)) \;{\text{\rm d}} \mathbf{x} = - \int_{U_0}  p_\mathbf{g}amma (\delta_i^{-h}(2\eta_0(\mathbf{n}abla\eta_0)\cdot(\delta_i^h\mathbf{v}_\mathbf{g}amma))) \;{\text{\rm d}} \mathbf{x}.
\]
Thus, \eqref{weakgeneralizedStokeswithpressureinteriorgammadifference} becomes
\mathbf{e}gin{equation}
\label{weakgeneralizedStokeswithpressureinteriorgammadifferenceb}
\mathbf{e}gin{aligned}
  \int_{U_0} & (\mathsf{D}(\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma)\mathsf{A} +\mathsf{A}\mathsf{D}(\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma)) : \left( \mathbf{n}abla (\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma) \right)\;{\text{\rm d}} \mathbf{x} \\
  & = \int_{U_0}  \mathbf{b}_\mathbf{g}amma \cdot (\delta_i^{-h}(\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma)) \;{\text{\rm d}} \mathbf{x} \\
  & \qquad + \int_{U_0} (\mathsf{S}((\mathbf{n}abla\eta_0) \otimes (\delta_i^h\mathbf{v}_\mathbf{g}amma))\mathsf{A} +\mathsf{A}\mathsf{S}((\mathbf{n}abla\eta_0) \otimes (\delta_i^h\mathbf{v}_\mathbf{g}amma))) : \left(\mathbf{n}abla (\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma) \right)\;{\text{\rm d}} \mathbf{x} \\ 
   & \qquad   - \int_{U_0} (\mathsf{D}(\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma)\mathsf{A} +\mathsf{A}\mathsf{D}(\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma)) : \left((\mathbf{n}abla\eta_0)\otimes (\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma) \right)\;{\text{\rm d}} \mathbf{x}  \\
   & \qquad   + \int_{U_0} (\mathsf{S}((\mathbf{n}abla\eta_0) \otimes (\delta_i^h\mathbf{v}_\mathbf{g}amma))\mathsf{A} +\mathsf{A}\mathsf{S}((\mathbf{n}abla\eta_0) \otimes (\delta_i^h\mathbf{v}_\mathbf{g}amma))) : \left((\mathbf{n}abla\eta_0)\otimes (\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma) \right)\;{\text{\rm d}} \mathbf{x} \\
   & \qquad  - \int_{U_0} (\mathsf{D}(\mathbf{v}_\mathbf{g}amma)(\delta_i^h\mathsf{A}) + (\delta_i^h\mathsf{A})\mathsf{D}(\mathbf{v}_\mathbf{g}amma)) : \mathbf{n}abla (\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma) \;{\text{\rm d}} \mathbf{x} \\
   & \qquad + 2 \int_{U_0}  p_\mathbf{g}amma (\delta_i^{-h}((\mathbf{n}abla\eta_0)\cdot(\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma))) \;{\text{\rm d}} \mathbf{x} \\
   & \qquad + \int_{U_0} \mathbf{b}f_\mathbf{g}amma \cdot (\delta_i^{-h}(\eta_0^2\delta_i^h\mathbf{v}_\mathbf{g}amma)) \;{\text{\rm d}} \mathbf{x}.
\end{aligned}
\end{equation}
Notice how similar \eqref{weakgeneralizedStokeswithpressureinteriorgammadifferenceb} is to \eqref{weakgeneralizedStokeswithpressureinteriorb}, with $\mathbf{v}$, $\mathbf{b}f$ and $p$ replaced by $\mathbf{v}_\mathbf{g}amma$, $\mathbf{b}f_\mathbf{g}amma$ and $p_\mathbf{g}amma$, the derivative $\partial_{x_i}$ replaced by the finite difference $\delta_i^h$, and the addition of the term with $\mathbf{b}_\mathbf{g}amma$. Using \eqref{finitedifferenceintestimate} and \eqref{finitedifferenceintestimateinfty}, we estimate \eqref{weakgeneralizedStokeswithpressureinteriorgammadifferenceb} just like \eqref{weakgeneralizedStokeswithpressureinteriorb}, except that we need to handle the extra term $\mathbf{b}_\mathbf{g}amma$. 

The idea is to consider $\mathbf{b}_\mathbf{g}amma$ as an extra forcing term. Once $\mathbf{b}_\mathbf{g}amma$ belongs to $L^2(\Omega)^3$, it can be treated just like $\mathbf{b}f_\mathbf{g}amma$. The fact that $\mathbf{b}_\mathbf{g}amma\in L^2(\Omega)^3$ will be discussed further on. 

Then, estimating the terms in \eqref{weakgeneralizedStokeswithpressureinteriorgammadifferenceb} like we did for \eqref{weakgeneralizedStokeswithpressureinteriorb}, we arrive at an interior estimate analogous to \eqref{interiorestimateh2c},
\mathbf{e}gin{multline}
  \label{interiorestimatehkplus2a}
  \|\mathbf{n}abla(\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma)\|_{L^2(\Omega)}^2 \leq \frac{c(\Omega)}{\alpha^2}\left( \|\mathbf{b}f_\mathbf{g}amma\|_{L^2(\Omega)}^2 + \|\mathbf{b}_\mathbf{g}amma\|_{L^2(\Omega)}^2 + \lambda_1\|p_\mathbf{g}amma\|_{L^2(\Omega)}^2 \right. \\
  \left. + \left(\lambda_1\|\mathsf{A}\|_{L^\infty(\Omega)}^2 + \|\mathbf{n}abla\mathsf{A}\|_{L^\infty(\Omega)}^2\right)\|\mathbf{n}abla \mathbf{v}_\mathbf{g}amma\|_{L^2(\Omega)}^2\right),
\end{multline}
for $i=1,2,3$, where $c(\Omega)$ is a shape constant. With $\mathbf{b}f_\mathbf{g}amma=D^\mathbf{g}amma\mathbf{b}f$ and $\mathbf{b}_\mathbf{g}amma$ in $L^2(\Omega)^3$, $\mathsf{A}\in W^{1,\infty}(\Omega)^{3\times 3}$, and $\mathbf{v}_\mathbf{g}amma = D^\mathbf{g}amma\mathbf{v} \in H^1(\Omega)^3$, the right hand side above is bounded and independent of $h$. Thus, $\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma$ has a subsequence that converges weakly in $H^1(\Omega)^3$. Since, on the other hand, $\mathbf{v}_\mathbf{g}amma$ belongs to $H^1(\Omega)$ and, hence, $\delta_i^h\mathbf{v}$ converges weakly to $\partial_{x_i}\mathbf{v}$ in $L^2_{\text{\rm loc}}(\Omega)^3$, then $\eta_0\delta_i^h\mathbf{v}_\mathbf{g}amma$ must converge in $H^1(\Omega)^3$ precisely to $\eta_0\partial_{x_i}\mathbf{v}_\mathbf{g}amma$, with $\eta_0\partial_{x_i}\mathbf{v}_\mathbf{g}amma\in H^1(\Omega)^3$. Since this applies to any $i=1,2,3$ and $\eta_0$ is of class ${\mathcal C}^2$, we find that $\eta_0\mathbf{v}_\mathbf{g}amma\in H^2(\Omega)^3$ and, passing to the limit as $h\rightarrow 0$ in \eqref{interiorestimatehkplus2a}, we find that
\mathbf{e}gin{multline}
  \label{interiorestimatehkplus2b}
  \|D^2(\eta_0\mathbf{v}_\mathbf{g}amma)\|_{L^2(\Omega)}^2 \leq \frac{c(\Omega)}{\alpha^2}\left( \|\mathbf{b}f_\mathbf{g}amma\|_{L^2(\Omega)}^2 + \|\mathbf{b}_\mathbf{g}amma\|_{L^2(\Omega)}^2 + \lambda_1\|p_\mathbf{g}amma\|_{L^2(\Omega)}^2 \right. \\
  \left. + \left(\lambda_1\|\mathsf{A}\|_{L^\infty(\Omega)}^2 + \|\mathbf{n}abla\mathsf{A}\|_{L^\infty(\Omega)}^2\right)\|\mathbf{n}abla \mathbf{v}_\mathbf{g}amma\|_{L^2(\Omega)}^2\right).
\end{multline}
  
  The estimates on the boundary are a bit more involved than in the case of low-order regularity. The difficulty is that, here, $\mathbf{v}_\mathbf{g}amma$ may not vanish on the boundary. Recall that in Theorem \ref{h2regularity} we have, instead of $\mathbf{v}_\mathbf{g}amma$, the vector field $\mathbf{v}$, which does vanish on the boundary. In that case, we could take any tangential derivative of $\mathbf{v}$ and still have it vanishing on the boundary.   Here, this procedure needs to be slightly modified, as follows.
  
  In general, the vector field $\mathbf{v}_\mathbf{g}amma$ does not vanish on the boundary because $\mathbf{g}amma$ may contain directions which are not tangential to $\partial \Omega$. Fortunately, it suffices to start by allowing $\mathbf{g}amma$ to contain only tangential derivatives. More precisely, if a multi-index $\mathbf{g}amma$ only contains tangential directions (say $\mathbf{g}amma\in \{1,2\}^k$ in the case the boundary has been rectified to $x_3=0$), then $\tilde\mathbf{v}_\mathbf{g}amma$ vanishes on the boundary. Hence, if $i=1$ or $2$, the finite difference $\delta_i^h\tilde\mathbf{v}_\mathbf{g}amma$ vanishes on the boundary, as well. Thus, by testing the localized and rectified equation against $\delta_i^{-h}(\eta_m\delta_i^h\tilde\mathbf{v}_\mathbf{g}amma)$, we obtain an $H^1$-estimate for $\eta_m\delta_i^h\tilde\mathbf{v}_\mathbf{g}amma$. This means that we obtain $L^2$ estimates for all the $(k+2)$-order derivatives with at most one derivative in the normal direction (the one coming from the $H^1$ norm). This leaves out all the $(k+2)$-order derivatives which have two or more derivatives in the normal direction. 
  
  Luckily, they can all be estimated directly from the localized and rectified version of the equation \eqref{generalizedStokesforderivativegamma}. This is done step by step. With $\tilde\mathbf{v}_\mathbf{g}amma$ having only tangential derivatives, we now obtain, directly and explicitly from the localized and rectified version of \eqref{generalizedStokesforderivativegamma}, an estimate for the second-order derivative of $\tilde\mathbf{v}_\mathbf{g}amma$ in the normal direction in terms of all the previously estimated second-order derivatives of $\mathbf{v}_\mathbf{g}amma$ with at most one normal direction. This gives an estimate for all the $(k+2)$-order derivatives with at most two derivatives in the normal direction. Then we allow $\tilde\mathbf{v}_\mathbf{g}amma$ to have at most one derivative in the normal direction and look again at \eqref{generalizedStokesforderivativegamma} to obtain estimates for all the $(k+2)$-order derivatives with at most three derivatives in the normal direction. We continue this way, until we find estimates for any $(k+2)$-order derivative. Coming back to the original variables, we obtain an estimate similar to \eqref{interiorestimatehkplus2a} and \eqref{estimateh2}, namely
\mathbf{e}gin{multline}
  \label{estimateh2vgamma}
  \|D^2 (\eta_m\mathbf{v}_\mathbf{g}amma)\|_{L^2(\Omega)}^2 \leq \frac{c(\Omega)}{\alpha^2}\sum_{|\hat\mathbf{g}amma|=k}\left( \|\mathbf{b}f_{\hat\mathbf{g}amma} + \mathbf{b}_{\hat\mathbf{g}amma}\|_{L^2(\Omega)}^2 + \lambda_1\|p_{\hat\mathbf{g}amma}\|_{L^2(\Omega)}^2 \right. \\
  \left. +\left(\lambda_1\|\mathsf{A}\|_{L^\infty(\Omega)}^2 + \|\mathbf{n}abla\mathsf{A}\|_{L^\infty(\Omega)}^2\right)\|D \mathbf{v}_{\hat\mathbf{g}amma}\|_{L^2(\Omega)}^2\right),
\end{multline}
where $c(\Omega)$ is a shape constant. Notice that, in the process described above, we need to consider all the multi-indices $\hat\mathbf{g}amma$ of order $k$, i.e. all the lower-order derivatives, to obtain the desired bound. This is reflected in \eqref{estimateh2vgamma} in the summation over all the multi-indices with $|\hat\mathbf{g}amma|=k$.

Combining all the localized estimates in the partition of unity we obtain a global bound for $\mathbf{v}_\mathbf{g}amma = \sum_m \eta_m \mathbf{v}_\mathbf{g}amma = \sum_m \eta_m D^\mathbf{g}amma \mathbf{v}$, for any $|\mathbf{g}amma|=k$, which can be written as
\mathbf{e}gin{multline}
  \label{estimatehkplus2}
  \|D^{k+2} \mathbf{v}\|_{L^2(\Omega)}^2 \leq \frac{c(\Omega)}{\alpha^2}\left( \|D^k\mathbf{b}f\|_{L^2(\Omega)}^2 + \sum_{|\mathbf{g}amma|=k}\|\mathbf{b}_\mathbf{g}amma\|_{L^2(\Omega)}^2 + \lambda_1\|D^k p\|_{L^2(\Omega)}^2 \right. \\
    \left. +\left(\lambda_1\|\mathsf{A}\|_{L^\infty(\Omega)}^2 + \|\mathbf{n}abla\mathsf{A}\|_{L^\infty(\Omega)}^2\right)\|D^{k+1} \mathbf{v}\|_{L^2(\Omega)}^2\right),
\end{multline}
where $c(\Omega)$ is another shape constant.    
    
Now it is time to estimate $D^k p$ and $\mathbf{b}_\mathbf{g}amma$, with $|\mathbf{g}amma|=k$, for $k\geq 1$. In the case $k=1$,  we have $\mathbf{g}amma = \mathbf{e}_1$, $\mathbf{e}_2$, or $\mathbf{e}_3$, and, moreover, $\mathbf{b}_{\mathbf{g}amma'} = 0$ and $\mathbf{v}_{\mathbf{g}amma'} = \mathbf{v}$. Thus, we have
\[ \mathbf{b}_{\mathbf{e}_j} = \mathbf{B}(\partial_{x_j}\mathsf{A},\mathbf{v}) = \mathbf{n}abla \cdot (\mathsf{D}(\mathbf{v})(\partial_{x_j}\mathsf{A}) + (\partial_{x_j}\mathsf{A})\mathsf{D}(\mathbf{v})).
\]
We estimate each $\mathbf{b}_{\mathbf{e}_j}$ using the inequality \eqref{l6inequality} and the assumption that $\mathsf{A}\in W^{1,\infty}(\Omega)^{3\times 3}\cap W^{2,3}(\Omega)^{3\times 3}$, so that
\mathbf{e}gin{multline}
 \label{estimatebbejinl2}
  \|\mathbf{b}_{\mathbf{e}_j}\|_{L^2(\Omega)} \leq  c\left(\|\mathbf{n}abla\mathsf{A}\|_{L^\infty(\Omega)}\|D^2 \mathbf{v}\|_{L^2(\Omega)} + \|D^2\mathsf{A}\|_{L^3(\Omega)}\|\mathbf{n}abla \mathbf{v}\|_{L^6(\Omega)}\right)  \\
      \leq  c(\Omega)\left(\|\mathbf{n}abla\mathsf{A}\|_{L^\infty(\Omega)} + \|D^2\mathsf{A}\|_{L^3(\Omega)}\right)\|D^2 \mathbf{v}\|_{L^2(\Omega)},
\end{multline} 
for $j=1,2,3,$ and another shape constant $c(\Omega)$.  For $k\geq 2$, we estimate $\mathbf{b}_\mathbf{g}amma$ recursively,
  \mathbf{e}gin{align*}
    \|\mathbf{b}_\mathbf{g}amma\|_{L^2(\Omega)} & \leq \|D \mathbf{b}_{\mathbf{g}amma'}\|_{L^2(\Omega)} + c\|D^2 \mathsf{A}\|_{L^3(\Omega)}\|D^k\mathbf{v}\|_{L^6(\Omega)} + c \|\mathbf{n}abla \mathsf{A}\|_{L^\infty(\Omega)}\|D^{k+1}\mathbf{v}\|_{L^2(\Omega)}\\ 
    & \leq \|D^2\mathbf{b}_{\mathbf{g}amma''}\|_{L^2(\Omega)} + c \|D^3 \mathsf{A}\|_{L^2(\Omega)}\|D^{k-1}\mathbf{v}\|_{L^\infty(\Omega)}\\
    & \qquad \qquad\qquad + c\|D^2 \mathsf{A}\|_{L^3(\Omega)}\|D^k\mathbf{v}\|_{L^6(\Omega)} + c \|\mathbf{n}abla \mathsf{A}\|_{L^\infty(\Omega)}\|D^{k+1}\mathbf{v}\|_{L^2(\Omega)} \\
    & \leq \ldots \leq c\|D^{k+1}\mathsf{A}\|_{L^2(\Omega)}\|\mathbf{n}abla \mathbf{v}\|_{L^\infty(\Omega)} +  ... + c \|D^3 \mathsf{A}\|_{L^2(\Omega)}\|D^{k-1}\mathbf{v}\|_{L^\infty(\Omega)}\\
    & \qquad  \qquad \qquad  + c\|D^2 \mathsf{A}\|_{L^3(\Omega)}\|D^k\mathbf{v}\|_{L^6(\Omega)} + c \|\mathbf{n}abla \mathsf{A}\|_{L^\infty(\Omega)}\|D^{k+1}\mathbf{v}\|_{L^2(\Omega)},
  \end{align*}  
  where $\mathbf{g}amma''$ is a multi-index one-degree lower than $\mathbf{g}amma'$, and so on. This can be written in the form
  \mathbf{e}gin{multline*}
    \|\mathbf{b}_\mathbf{g}amma\|_{L^2(\Omega)} \leq c\sum_{j=1}^{k-1} \|D^{k+2-j}\mathsf{A}\|_{L^2(\Omega)}\|D^{j}\mathbf{v}\|_{L^\infty(\Omega)} \\
    + c\|D^2 \mathsf{A}\|_{L^3(\Omega)}\|D^k\mathbf{v}\|_{L^6(\Omega)} + c \|\mathbf{n}abla \mathsf{A}\|_{L^\infty(\Omega)}\|D^{k+1}\mathbf{v}\|_{L^2(\Omega)}.
  \end{multline*}
  
  Using  the Agmon-type inequality \eqref{agmoninequality} and the inequality \eqref{l6inequality}, we see that
  \mathbf{e}gin{multline}
    \label{estimatebgamma}
    \|\mathbf{b}_\mathbf{g}amma\|_{L^2(\Omega)} \leq  \frac{c(\Omega)}{\lambda_1^{1/4}}\sum_{j=1}^{k-1} \|D^{k+2-j}\mathsf{A}\|_{L^2(\Omega)}\|D^{j+2}\mathbf{v}\|_{L^2(\Omega)} \\
    + c(\Omega)\left(\|D^2 \mathsf{A}\|_{L^3(\Omega)}+ \|\mathbf{n}abla \mathsf{A}\|_{L^\infty(\Omega)}\right)\|D^{k+1}\mathbf{v}\|_{L^2(\Omega)}.
  \end{multline}
  This is valid for any multi-index $\mathbf{g}amma$ of order $k\geq 2$. In fact, formula \eqref{estimatebgamma} is valid even for $k=1$, since in this case the summation is empty and the estimate reduces to \eqref{estimatebbejinl2}. 
    
  Using estimate \eqref{estimatebgamma}, inequality \eqref{estimatehkplus2} yields
  \mathbf{e}gin{multline}
  \label{estimatehkplus2b}
  \|D^{k+2} \mathbf{v}\|_{L^2(\Omega)}^2 \leq \frac{c(\Omega)}{\alpha^2}\left( \|D^k\mathbf{b}f\|_{L^2(\Omega)}^2 + \frac{c(\Omega)}{\lambda_1^{1/2}}\sum_{j=1}^{k-1} \|D^{k+2-j}\mathsf{A}\|_{L^2(\Omega)}^2\|D^{j+2}\mathbf{v}\|_{L^2(\Omega)}^2 \right. \\
  \left. + \lambda_1\|D^k p\|_{L^2(\Omega)}^2 + \left(\|D^2 \mathsf{A}\|_{L^3(\Omega)}^2 + \|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}^2\right)\|D^{k+1} \mathbf{v}\|_{L^2(\Omega)}^2\right),
\end{multline}
  for $k\geq 1$.

The pressure term in \eqref{estimatehkplus2b} can be obtained by taking the derivative of \eqref{generalizedStokes} up to order $k-1$ (one degree lower than \eqref{generalizedStokesforderivativegamma}, i.e. with $\mathbf{g}amma$ replaced by $\mathbf{g}amma'$, with $|\mathbf{g}amma'|=k-1$). In doing so, we find equations for $D^k p$ valid in $L^2(\Omega)^3$, which yield an estimate analogous to \eqref{basicpestimatefromStokes}, namely
\mathbf{e}gin{multline}
  \label{estimatesdkplusoneofp}
  \|D^k p\|_{L^2(\Omega)} \leq \|D^{k-1}\mathbf{b}f\|_{L^2(\Omega)} + \sum_{|\mathbf{g}amma'|=k-1} \|\mathbf{b}_{\mathbf{g}amma'}\|_{L^2(\Omega)} + c\|\mathsf{A}\|_{L^\infty(\Omega)}\|D^{k+1}\mathbf{v}\|_{L^2(\Omega)} \\
   + c\|\mathbf{n}abla \mathsf{A}\|_{L^\infty(\Omega)}\|D^k\mathbf{v}\|_{L^2(\Omega)}.
\end{multline}
Using the Poincar\'e inequality \eqref{poincareineq} and the estimate \eqref{estimatebgamma} with $\mathbf{g}amma'$ instead of $\mathbf{g}amma$ yields
\mathbf{e}gin{multline}
  \label{estimatesdkplusoneofpb}
  \|D^k p\|_{L^2(\Omega)} \leq \|D^{k-1}\mathbf{b}f\|_{L^2(\Omega)} + \frac{c(\Omega)}{\lambda_1^{1/4}}\sum_{j=1}^{k-2} \|D^{k+1-j}\mathsf{A}\|_{L^2(\Omega)}\|D^{j+2}\mathbf{v}\|_{L^2(\Omega)} \\
    + c(\Omega)\left(\|D^2 \mathsf{A}\|_{L^3(\Omega)}+ \|\mathbf{n}abla \mathsf{A}\|_{L^\infty(\Omega)}\right)\|D^k\mathbf{v}\|_{L^2(\Omega)} + c\|\mathsf{A}\|_{L^\infty(\Omega)}\|D^{k+1}\mathbf{v}\|_{L^2(\Omega)}.
\end{multline}
Inserting \eqref{estimatesdkplusoneofpb} into \eqref{estimatehkplus2b} and using the corresponding estimate \eqref{estimatebgamma} with $\mathbf{g}amma'$ instead of $\mathbf{g}amma$ and using the Poincar\'e inequality \eqref{poincareineq} yields
\mathbf{e}gin{multline}
  \label{estimatehkplus2c}
  \|D^{k+2} \mathbf{v}\|_{L^2(\Omega)}^2 \leq \frac{c(\Omega)}{\alpha^2}\left(\|\mathbf{b}f\|_{H^k_{\lambda_1}(\Omega)}^2 + \frac{1}{\lambda_1^{1/2}}\sum_{j=1}^{k-1} \|\mathsf{A}\|_{H^{k+2-j}_{\lambda_1}(\Omega)}^2\|D^{j+2}\mathbf{v}\|_{L^2(\Omega)}^2 \right. \\
  \left. + \left(\|D^2 \mathsf{A}\|_{L^3(\Omega)}^2 + \|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}^2\right)\|D^{k+1} \mathbf{v}\|_{L^2(\Omega)}^2\right).
\end{multline}
Due to the Sobolev embeddings in our three-dimensional space, we treat the cases $k=1$ and $k\geq 2$ differently. For $k=1$, we find, from \eqref{estimatehkplus2c}, that
\mathbf{e}gin{equation}
  \label{estimateh3b}
   \|D^3 \mathbf{v}\|_{L^2(\Omega)}^2 \leq \frac{c(\Omega)}{\alpha^2}\left(\|\mathbf{b}f\|_{H^1_{\lambda_1}(\Omega)}^2 + \left(\|D^2 \mathsf{A}\|_{L^3(\Omega)}^2 + \|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}^2\right)\|D^2 \mathbf{v}\|_{L^2(\Omega)}^2\right).
\end{equation}
Estimate \eqref{estimateh3b} together with the estimate \eqref{h2regularityforh2v} for $D^2 \mathbf{v}$ yields
\mathbf{e}gin{multline}
    \label{estimateh3c}   
        \|D^3\mathbf{v}\|_{L^2(\Omega)}^2 \leq \frac{c(\Omega)}{\alpha^2}\left( \|\mathbf{b}f\|_{H^1_{\lambda_1}(\Omega)}^2  \right. \\
    \left. + \frac{1}{\alpha^2}\left(\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}^2 + \|D^2\mathsf{A}\|_{L^3(\Omega)}^2\right)\left( \|\mathbf{b}f\|_{L^2(\Omega)}^2 + \frac{1}{\alpha^2}\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}^2\|\mathbf{b}f\|_{H^{-1}(\Omega)}^2\right)\right),
      \end{multline}
which proves \eqref{hkplus2regularityforv} in the case $k=1$. Using \eqref{estimatesdkplusoneofpb} with $k=2$ together with \eqref{estimateh3b} we find that
\mathbf{e}gin{multline}
    \label{h2regularityforp}
    \|D^2 p\|_{L^2(\Omega)} \leq \frac{c(\Omega)}{\alpha}\|\mathsf{A}\|_{L^\infty(\Omega)}\left( \|\mathbf{b}f\|_{H^1_{\lambda_1}(\Omega)}  \right. \\
    \left. + \frac{1}{\alpha}\left(\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)} + \|D^2\mathsf{A}\|_{L^3(\Omega)}\right)\left( \|\mathbf{b}f\|_{L^2(\Omega)} + \frac{1}{\alpha}\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}\|\mathbf{b}f\|_{H^{-1}(\Omega)}\right)\right),
\end{multline}
which proves the estimate \eqref{hkplus1regularityforp} for the pressure in the case $k=1$.

In the case $k\geq 2$, since $H^3(\Omega)$ is included in both $W^{2,3}(\Omega)$ and $W^{1,\infty}(\Omega)$, the last term in \eqref{estimatehkplus2c} is bounded by the term with $j=k-1$ in the summation. Thus, we rewrite \eqref{estimatehkplus2c} simply as
\mathbf{e}gin{equation}
  \label{estimatehkplus2d}
  \|D^{k+2} \mathbf{v}\|_{L^2(\Omega)}^2 \leq \frac{c(\Omega)}{\alpha^2}\left(\|\mathbf{b}f\|_{H^k_{\lambda_1}(\Omega)}^2 + \frac{1}{\lambda_1^{1/2}}\sum_{j=1}^{k-1} \|\mathsf{A}\|_{H^{k+2-j}_{\lambda_1}(\Omega)}^2\|D^{j+2}\mathbf{v}\|_{L^2(\Omega)}^2 \right). 
\end{equation}
for a different constant $c(\Omega)$.  Estimate \eqref{estimatehkplus2d} is of the form
\mathbf{e}gin{equation}
  \label{vrecursion}
   v_{k+2} \leq c(\Omega)\left(d_k +  \sum_{j=1}^{k-1} a_{k+2-j}v_{j+2}\right),
\end{equation}
for $k\geq 2$, where 
\[ v_j = \|D^j\mathbf{v}\|_{L^2(\Omega)}^2, \quad d_j = \frac{1}{\alpha^2} \|\mathbf{b}f\|_{H^j_{\lambda_1}(\Omega)}^2, \quad a_j = \frac{1}{\alpha^2}\frac{1}{\lambda_1^{1/2}}\|\mathsf{A}\|_{H^j_{\lambda_1}(\Omega)}^2,
\]
for $j\geq 2$. We include the case $k=1$ in \eqref{vrecursion} by setting $d_1 = v_3$. 
Then, we show by induction that
\mathbf{e}gin{equation}
  \label{recursioninequalitycomb}
  v_{k+2} \leq c(\Omega) \left(d_k + \sum_{j=1}^{k-1} \left( \sum_{\substack{\alpha_1, \ldots, \alpha_{k-j}\in {\mathbb{N}}_0 \\ \alpha_1 + \ldots + (k-j)\alpha_{k-j} = k-j}} c_{\alpha_1,\ldots,\alpha_{k-j}}^k a_3^{\alpha_1}\cdots a_{k-j+2}^{\alpha_{k-j}}\right)d_j \right),
\end{equation}
for $k\geq 1$, where $c_{\alpha_1,\ldots,\alpha_{k-j}}^k$ are suitable numerical constants given recursively by
\mathbf{e}gin{align}
  \label{combinatoricconstants1} 
  & c_{0,\ldots, 0, 1}^k = 1, \\
  \label{combinatoricconstants2} 
  & c_{\alpha_1,\ldots,\alpha_{k-j -1},0}^k = \sum_{\ell=j+1}^{k-1} \sum_{\substack{\tilde\alpha_1, \ldots, \tilde\alpha_{\ell-j}\in {\mathbb{N}}_0 \\ \tilde\alpha_1 + \ldots + (\ell-j)\tilde\alpha_{\ell-j} = \ell-j}} c_{\tilde\alpha_1,\ldots,\tilde\alpha_{\ell-j}}^\ell,
\end{align}
for all $k\geq 1$ and $j=1,\ldots, k-1$. Moreover, in the formula \eqref{combinatoricconstants2}, the exponents $\alpha_1, \ldots, \alpha_{k-j-1}$ are related to $\ell$ and $\tilde\alpha_1, \ldots, \tilde\alpha_{\ell-j}$ by the (one-to-one) correspondences
\mathbf{e}gin{equation}
  \label{alpharelations}
  \alpha_m = \mathbf{e}gin{cases} 
      \tilde\alpha_m, & 1 \leq m \leq \ell-j, \; m\neq k-\ell, \\
      \tilde\alpha_{k-\ell} + 1, &  m = k-\ell, \text{ if } k-\ell \leq \ell - j, \\
      1, &  m = k-\ell, \text{ if } k-\ell > \ell - j, \\
      0, &  \ell-j < m \leq k-j, \; m\neq k-\ell,
    \end{cases}
\end{equation}
and
\mathbf{e}gin{equation}
  \label{alpharelations2}
  \tilde\alpha_m = \mathbf{e}gin{cases} 
      \alpha_m, & 1 \leq m \leq \ell-j, \; m\neq k-\ell, \\
      \alpha_{k-\ell} - 1, &  m = k-\ell, \text{ if } k-\ell \leq \ell - j. \\
    \end{cases}
\end{equation}
More precisely, for $k\geq 1$ and $j=1, \ldots, k-1$, the exponents $\tilde\alpha_1, \ldots, \tilde\alpha_{\ell-j}\in {\mathbb{N}}_0$ satisfy the criteria $\tilde\alpha_1 + \ldots + (\ell-j)\tilde\alpha_{\ell-j} = \ell-j$, for $\ell\in \{j+1,\ldots, k_1\}$, if, and only if, the exponents $\alpha_1,\ldots \alpha_{k-j-1}\in {\mathbb{N}}_0$ satisfy the criteria $\alpha_1 + \ldots + (k-j-1)\alpha_{k-j-1} = k-j$..

Indeed, this is trivially true for $k=1$. For $k\geq 2$, assuming it is true up to $k-1$, we find
\mathbf{e}gin{align*}
 v_{k+2} & \leq c(\Omega)\left(d_k  \right. \\
   &  \qquad  \left. + \sum_{j=1}^{k-1} a_{k+2-j}\left(d_j + \sum_{\ell=1}^{j-1} \left( \sum_{\substack{\tilde\alpha_1, \ldots, \tilde\alpha_{j-\ell}\in {\mathbb{N}}_0 \\ \tilde\alpha_1 + \ldots + (j-\ell)\tilde\alpha_{j-\ell} = j-\ell}} c_{\tilde\alpha_1,\ldots,\tilde\alpha_{j-\ell}}^ja_3^{\tilde\alpha_1}\cdots a_{j-\ell+2}^{\tilde\alpha_{j-\ell}}\right)d_\ell\right)\right) \\
   & = c(\Omega)\left(d_k +  \sum_{j=1}^{k-1} a_{k+2-j}d_j \right. \\
    & \qquad  \qquad \left. + \sum_{\ell=1}^{k-1} \sum_{j=\ell+1}^{k-1} a_{k+2-j} \left( \sum_{\substack{\tilde\alpha_1, \ldots, \tilde\alpha_{j-\ell}\in {\mathbb{N}}_0 \\ \tilde\alpha_1 + \ldots + (j-\ell)\tilde\alpha_{j-\ell} = j-\ell}} c_{\tilde\alpha_1,\ldots,\tilde\alpha_{j-\ell}}^j a_3^{\tilde\alpha_1}\cdots a_{j-\ell+2}^{\tilde\alpha_{j-\ell}}\right)d_\ell\right) \\
   & = c(\Omega)\left(d_k +  \sum_{j=1}^{k-1} a_{k+2-j}d_j \right. \\
    & \qquad  \qquad \left. + \sum_{j=1}^{k-1} \sum_{\ell=j+1}^{k-1} a_{k+2-\ell} \left( \sum_{\substack{\tilde\alpha_1, \ldots, \tilde\alpha_{\ell-j}\in {\mathbb{N}}_0 \\ \tilde\alpha_1 + \ldots + (\ell-j)\tilde\alpha_{\ell-j} = \ell-j}} c_{\tilde\alpha_1,\ldots,\tilde\alpha_{\ell-j}}^\ell a_3^{\tilde\alpha_1}\cdots a_{\ell-j+2}^{\tilde\alpha_{\ell-j}}\right)d_j\right).
\end{align*}
Putting $d_j$ in evidence we rewrite the estimate as
\mathbf{e}gin{multline}
  \label{recursioninequalitycombpre}
  v_{k+2} \leq c(\Omega)d_k \\
   + c(\Omega)\sum_{j=1}^{k-1} \left(a_{k+2-j} + \sum_{\ell=j+1}^{k-1} a_{k+2-\ell} \left( \sum_{\substack{\tilde\alpha_1, \ldots, \tilde\alpha_{\ell-j}\in {\mathbb{N}}_0 \\ \tilde\alpha_1 + \ldots + (\ell-j)\tilde\alpha_{\ell-j} = \ell-j}} c_{\tilde\alpha_1,\ldots,\tilde\alpha_{\ell-j}}^\ell a_3^{\tilde\alpha_1}\cdots a_{\ell-j+2}^{\tilde\alpha_{\ell-j}}\right)\right)d_j.
\end{multline}
We argue that
\mathbf{e}gin{multline}
  \label{combinatorics}
  a_{k+2-j} + \sum_{\ell=j+1}^{k-1} a_{k+2-\ell} \left( \sum_{\substack{\tilde\alpha_1, \ldots, \tilde\alpha_{\ell-j}\in {\mathbb{N}}_0 \\ \tilde\alpha_1 + \ldots + (\ell-j)\tilde\alpha_{\ell-j} = \ell-j}} c_{\tilde\alpha_1,\ldots,\tilde\alpha_{\ell-j}}^\ell a_3^{\tilde\alpha_1}\cdots a_{\ell-j+2}^{\tilde\alpha_{\ell-j}}\right) \\
  = \sum_{\substack{\alpha_1, \ldots, \alpha_{k-j}\in {\mathbb{N}}_0 \\ \alpha_1 + \ldots + (k-j)\alpha_{k-j} = k-j}} c_{\alpha_1,\ldots,\alpha_{k-j}}^ka_3^{\alpha_1}\cdots a_{k-j+2}^{\alpha_{k-j}},
\end{multline}
with the numerical constants satisfying \eqref{combinatoricconstants1} and \eqref{combinatoricconstants2} and the exponents satisfying \eqref{alpharelations} and \eqref{alpharelations2}.

Indeed, the first term $a_{k+2-j}$ is in the summation term in the right hand side with $\alpha_{k+j} = 1$ and the remaining exponents $\alpha_\ell =0$, $1\leq \ell < k-j$, with the associated numerical constant being 
\[ c_{0,\ldots, 0, 1}^k=1. 
\]
This is the only possible combination that contains the ``highest-order'' term $a_{k+2-j}$. 

For $\ell = j+1, \ldots, k-1$ and any combination $\tilde\alpha_1, \ldots, \tilde\alpha_{\ell-j}\in {\mathbb{N}}_0$ with $\tilde\alpha_1 + \ldots + (\ell-j)\tilde\alpha_{\ell-j} = \ell-j$, we have
\[ \tilde\alpha_1 + \ldots + (\ell-j)\tilde\alpha_{\ell-j} + (k-\ell) = \ell-j + k-\ell = k-j.
\]
Notice that $k-\ell$ may, or may not, be less than $\ell-j$. Thus, we take
\mathbf{e}gin{equation}
  \label{alphaaddone}
  \tilde\alpha_{k-\ell} = \mathbf{e}gin{cases} \alpha_{k-\ell}+1, & \text{if } k-\ell \leq \ell-j, \\ 1, & \text{ otherwise}, \end{cases}
\end{equation}
and keep the remaining exponents the same, i.e. $\alpha_m = \tilde\alpha_m$ if $m\leq \ell-j$, $m\neq k-\ell$, and $\alpha_m = 0$, if $\ell-j < m \leq k-1$, $m\neq k-\ell$. In doing so, the exponents $\alpha_1, \ldots \alpha_{k-j}$ satisfy
\[ \alpha_1 + \ldots + (k-j)\alpha_{k-j} = k-j, 
\]
with $\alpha_{k-j}=0$. Hence, this combination of exponents is in the right hand side as well, and, from \eqref{combinatorics}, we see that the associated constant $c_{\alpha_1,\ldots, \alpha_{k-j}}^k$ is given according to \eqref{combinatoricconstants2}.

On the other hand, if the exponents $\alpha_1, \ldots, \alpha_{k-j}\in {\mathbb{N}}_0$ satisfy $\alpha_1 + \ldots + (k-j)\alpha_{k-j} = k-j$, then either $\alpha_{k-j}=1$ or $0$. If $\alpha_{k-j}=1$, the remaining exponents are zero, so that this term is in the left hand side of \eqref{combinatorics} and the associated numerical constant is $c_{0,\ldots,1}=1$. 

In the case $\alpha_{k-j} = 0$, then $\alpha_{k-\ell} \neq 0$, for one or more $j + 1 \leq \ell \leq k-1$. For any such $\ell$, we see that
\[ \alpha_1 + \ldots + (k-\ell)(\alpha_{k-\ell}-1)  = k-j - (k-\ell) = \ell - j.
\]
Thus, we define $\tilde\alpha_m = \alpha_m$, for $1\leq m \leq \ell-j$, with $m\neq k-\ell$, and $\tilde\alpha_m = \alpha_m-1$, if $m=k-\ell\leq \ell-j$ (if $k-\ell> \ell-j$, there is no need to define the corresponding $\tilde\alpha_m$; see \eqref{alpharelations2}). In doing so, the exponents $\tilde\alpha_1, \ldots, \tilde\alpha_{\ell-j}$ satisfy the condition to be in the summation term in the left hand side of \eqref{combinatorics}. 

This completes the proof of \eqref{combinatorics}, with the numerical constant satisfying \eqref{combinatoricconstants1} and \eqref{combinatoricconstants2} and the exponents satisfying \eqref{alpharelations} and \eqref{alpharelations2}. Combining \eqref{recursioninequalitycombpre} and \eqref{combinatorics} we deduce \eqref{recursioninequalitycomb}.

The estimate \eqref{recursioninequalitycomb} is the sharpest we can get from \eqref{vrecursion}. But in order to have a more readable and manageable inequality, we proceed to estimate the summation terms in \eqref{recursioninequalitycomb}. We divide the combinatorial condition $\sum_{i=1}^{k-j} i\alpha_i = k-j$ by $k-j$ to see that
\[ \sum_{i=1}^{k-j} \frac{i\alpha_i}{k-j} = 1,
\]
with $0 \leq i\alpha_i/(k-j) \leq 1$, for all $i=1,\ldots, k-j$. Using Young's inequality, we find that
\mathbf{e}gin{equation}
  \label{combinatorialtermestimate}
  \prod_{i=1}^{k-j} a_{i+2}^{\alpha_i} \leq \sum_{i=1}^{k-j} \frac{i\alpha_i}{k-j} a_{i+2}^{(k-j)/i} \leq \sum_{i=1}^{k-j} a_{i+2}^{(k-j)/i}.
\end{equation}
Now, applying estimate \eqref{combinatorialtermestimate} to \eqref{recursioninequalitycomb} gives us
\mathbf{e}gin{equation}
  \label{recursioninequalityest}
  v_{k+2} \leq c(\Omega) \left(d_k + c_k \sum_{j=1}^{k-1} \left( \sum_{i=1}^{k-j} a_{i+2}^{(k-j)/i}\right)d_j \right),
\end{equation}
where
\mathbf{e}gin{equation}
  \label{defck}
  c_k = \sum_{\substack{\alpha_1, \ldots, \alpha_{k-j}\in {\mathbb{N}}_0 \\ \alpha_1 + \ldots + (k-j)\alpha_{k-j} = k-j}} c_{\alpha_1,\ldots,\alpha_{k-j}}^k.
\end{equation}
Going back to the definitions of $v_j$, $d_j$ and $a_j$, not forgetting that $d_1=v_3$, we find that
\mathbf{e}gin{multline}
  \label{recursioninequalityestfull}
  \|D^{k+2}\mathbf{v}\|_{L^2(\Omega)}^2 \leq c_k(\Omega) \left( \frac{1}{\alpha^2}\|\mathbf{b}f\|_{H^k_{\lambda_1}(\Omega)}^2 + \frac{1}{\alpha^2}\sum_{j=2}^{k-1} \left( \sum_{i=1}^{k-j} \frac{1}{\alpha^{2(k-j)/i}}\frac{1}{\lambda_1^{(k-j)/2i}}\|\mathsf{A}\|_{H^{i+2}_{\lambda_1}(\Omega)}^{2(k-j)/i}\right)  \|\mathbf{b}f\|_{H^j_{\lambda_1}(\Omega)}^2 \right. \\
  \left. + \left( \sum_{i=1}^{k-1} \frac{1}{\alpha^{2(k-1)/i}}\frac{1}{\lambda_1^{(k-1)/2i}}\|\mathsf{A}\|_{H^{i+2}_{\lambda_1}(\Omega)}^{2(k-1)/i}\right)\|D^3\mathbf{v}\|_{L^2(\Omega)}^2\right),
\end{multline}
for suitable non-dimensional shape constants $c_k(\Omega)$, for any $k\geq 2$. Using the estimate \eqref{estimateh3c} for $v_3$ we arrive at 
\mathbf{e}gin{equation} 
  \label{hkplus2regularityforvfull}
  \mathbf{e}gin{aligned}
   \|D^{k+2} & \mathbf{v}\|_{L^2(\Omega)}^2 \leq \frac{c_k(\Omega)}{\alpha^2} \left( \|\mathbf{b}f\|_{H^k_{\lambda_1}(\Omega)}^2 +  \sum_{j=2}^{k-1} \left( \sum_{i=1}^{k-j} \frac{1}{\alpha^{2(k-j)/i}}\frac{1}{\lambda_1^{(k-j)/2i}}\|\mathsf{A}\|_{H^{i+2}_{\lambda_1}(\Omega)}^{2(k-j)/i}\right)\|\mathbf{b}f\|_{H^j_{\lambda_1}(\Omega)}^2 \right. \\
   & + \left( \sum_{i=1}^{k-1} \frac{1}{\alpha^{2(k-1)/i}}\frac{1}{\lambda_1^{(k-1)/2i}}\|\mathsf{A}\|_{H^{i+2}_{\lambda_1}(\Omega)}^{2(k-1)/i}\right)\left( \|\mathbf{b}f\|_{H^1_{\lambda_1}(\Omega)}^2 \right. \\
  & \left. \quad \quad + \frac{1}{\alpha^2}\left( \|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}^2 + \|D^2\mathsf{A}\|_{L^3(\Omega)}^2\right)\left(\|\mathbf{b}f\|_{L^2(\Omega)}^2 + \frac{1}{\alpha^2}\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}^2 \|\mathbf{b}f\|_{H^{-1}(\Omega)}^2\right)\right).
  \end{aligned}
\end{equation}  
This proves \eqref{Rkdef}. The estimate \eqref{hkplus1regularityforp} for the pressure follows from \eqref{estimatesdkplusoneofpb} and \eqref{hkplus2regularityforv}.
\end{proof}

We end this section with some remarks. 

 \mathbf{e}gin{rmk}
   The estimates \eqref{weaksolutionweakgeneralizedStokesestimaterepeat},  \eqref{h2regularityforh2v}, \eqref{h3egularityforv}, and \eqref{hkplus2regularityforv} for the velocity field (and similarly for the pressure estimates \eqref{weaksolutionweakgeneralizedStokesestimatewithpressure}, \eqref{h2regularityforh1p}, \eqref{hkplus2regularityforv}, and \eqref{hkplus1regularityforp}) can be put into a single expression of the form
   \mathbf{e}gin{equation} 
     \label{hkplustworegularityshort}
     \|D^{k+2}\mathbf{v}\|_{L^2(\Omega)} \leq \frac{c(\Omega)}{\alpha} \sum_{j=-1}^k \lambda_1^{(k-j)/2}a_{k,j}(\alpha,\lambda_1,\mathsf{A}) \|\mathbf{b}f\|_{H^j(\Omega)},
   \end{equation}
   for all integers $k\geq -1$, where the coefficients $a_{k,j}$ are non-dimensional quantities. In the case $k\geq 2$, these quantities are given by
   \mathbf{e}gin{equation}
   \mathbf{e}gin{aligned}
      & a_{k,k}(\alpha,\lambda_1,\mathsf{A}) = 1, \\
      & a_{k,j}(\alpha,\lambda_1,\mathsf{A}) =  \sum_{i=1}^{k-j} \frac{1}{\alpha^{(k-j)/i}}\frac{1}{\lambda_1^{(2i+1)(k-j)/4i}}\|\mathsf{A}\|_{H^{i+2}_{\lambda_1}(\Omega)}^{(k-j)/i}, \quad 1 \leq j \leq k-1, \\
      & a_{k,0}(\alpha,\lambda_1,\mathsf{A}) = \frac{1}{\alpha\lambda_1^{1/2}} \left(\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)} + \|D^2\mathsf{A}\|_{L^3(\Omega)}\right) a_{k,1}(\alpha,\lambda_1,\mathsf{A}), \\
     & \quad = \frac{1}{\alpha\lambda_1^{1/2}}\left(\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)} + \|D^2\mathsf{A}\|_{L^3(\Omega)}\right)\sum_{i=1}^{k-1} \frac{1}{\alpha^{(k-j)/i}}\frac{1}{\lambda_1^{(2i+1)(k-j)/4i}}\|\mathsf{A}\|_{H^{i+2}_{\lambda_1}(\Omega)}^{(k-j)/i}, \\
     & a_{k,-1}(\alpha,\lambda_1,\mathsf{A}) = \frac{1}{\alpha\lambda_1^{1/2}}\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}a_{k,0}(\alpha,\lambda_1,\mathsf{A}),  \\
     & \quad =  \frac{1}{\alpha^2\lambda_1}\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}\left(\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)} + \|D^2\mathsf{A}\|_{L^3(\Omega)}\right)\sum_{i=1}^{k-1} \frac{1}{\alpha^{(k-j)/i}}\frac{1}{\lambda_1^{(2i+1)(k-j)/4i}}\|\mathsf{A}\|_{H^{i+2}_{\lambda_1}(\Omega)}^{(k-j)/i},
   \end{aligned}
   \end{equation}
   while for $k=-1, 0, 1$, we have
   \mathbf{e}gin{equation}
     \mathbf{e}gin{aligned}
      & a_{-1,-1}(\alpha,\lambda_1,\mathsf{A}) = a_{0,0}(\alpha,\lambda_1,\mathsf{A}) = a_{1,1}(\alpha,\lambda_1,\mathsf{A}) = 1, \\
      & a_{0,-1}(\alpha,\lambda_1,\mathsf{A}) = \frac{1}{\alpha\lambda_1^{1/2}}\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}, \\
      & a_{1,0}(\alpha,\lambda_1,\mathsf{A}) = \frac{1}{\alpha\lambda_1^{1/2}}\left(\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)} + \|D^2\mathsf{A}\|_{L^3(\Omega)}\right), \\
      & a_{1,-1}(\alpha,\lambda_1,\mathsf{A}) = \frac{1}{\alpha^2\lambda_1}\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}\left(\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)} + \|D^2\mathsf{A}\|_{L^3(\Omega)}\right).
    \end{aligned}
  \end{equation}
  Notice, from \eqref{generalizedStokes}, that $\mathbf{b}f$ has the physical dimension (scaling) of $\alpha \lambda_1\mathbf{v}$, while $D^j$ has the physical dimension of $\lambda_1^{j/2}$, so it is interesting to rewrite \eqref{hkplustworegularityshort} as
   \mathbf{e}gin{equation} 
     \label{hkplustworegularityshort2}
     \frac{1}{\lambda_1^{(k+2)/2}}\|D^{k+2}\mathbf{v}\|_{L^2(\Omega)} \leq c(\Omega)\sum_{j=-1}^k a_{k,j}(\alpha,\lambda_1,\mathsf{A}) \frac{ \|\mathbf{b}f\|_{H^j(\Omega)}}{\alpha \lambda_1^{(j+2)/2}},
   \end{equation}    
 \end{rmk} 
 
 \mathbf{e}gin{rmk}
 We may use the interpolation inequality \eqref{interpolationineq} to slightly rewrite the summation in the right hand side of estimate \eqref{combinatorialtermestimate}, which ends up in the estimate \eqref{recursioninequalityest} and, consequently, in the final estimates \eqref{hkplus2regularityforv} and \eqref{hkplus1regularityforp} with \eqref{Rkdef}. Using \eqref{interpolationineq}, we see that (up to a multiplicative shape constant that will be accounted for at the end of the argument but which is omitted here for notational simplicity)
\[ a_{i+2} \leq a_3^{\theta_i}a_{k-j+2}^{1-\theta_i},
\]
where $\theta_i$ is given by
\mathbf{e}gin{equation}
  \label{deftheta}
  i+2 = 3\theta_i + (1-\theta_i)(k-j+2).
\end{equation}
We use Young's inequality to write
\mathbf{e}gin{equation}
  \label{aiafterinterpolationandyoung}
  a_{i+2}^{(k-j)/i} \leq a_3^{\theta_i (k-j)/i}a_{k-j+2}^{(1-\theta_i)(k-j)/i} \leq \frac{1}{p}a_3^{\theta_i (k-j)p/i} + \frac{1}{q}a_{k-j+2}^{(1-\theta_i)(k-j)q/i},
\end{equation}
for any $1\leq p, q \leq \infty$, with $p^{-1} + q^{-1} = 1$. Choosing $q$ such that 
\mathbf{e}gin{equation}
  \label{youngdefq}
  \frac{(1-\theta_i)(k-j)q}{i} = 1,
\end{equation}
we find that
\[ \frac{1}{p} = 1 - \frac{(1-\theta_i)(k-j)}{i} = \frac{i - (1-\theta_i)(k-j)}{i}.
\]
Thus, one easily checks, using the value of $\theta_i$ given in \eqref{deftheta}, that
\mathbf{e}gin{equation}
  \label{youngfindp}
  \frac{\theta_i (k-j)p}{i} = \frac{\theta_i (k-j)}{i - (1-\theta_i)(k-j)} = k-j.
\end{equation}
Combining \eqref{aiafterinterpolationandyoung}, \eqref{youngdefq} and \eqref{youngfindp} together, we see that
\[ a_{i+2}^{(k-j)/i} \leq a_3^{k-j} + a_{k-j+2}.
\]
Thus, the right hand side of \eqref{combinatorialtermestimate} is estimated as
\mathbf{e}gin{equation}
  \label{prodtoyoungineq}
  \sum_{i=1}^{k-j} a_{i+2}^{(k-j)/i} \leq c(\Omega)(k-j) \left(a_3^{k-j} + a_{k-j}\right),
\end{equation}
where now we include the suitable constant $c(\Omega)$ coming from the interpolations. Applying estimate \eqref{prodtoyoungineq} to \eqref{recursioninequalitycomb} gives us
\mathbf{e}gin{equation}
  \label{recursioninequalityestalt}
  v_{k+2} \leq c(\Omega) \left(d_k + kc_k \sum_{j=1}^{k-1} \left( a_3^{k-j} + a_{k-j}\right)d_j \right),
\end{equation}
for another suitable shape constant $c(\Omega)$, and where $c_k$ is given in \eqref{defck}. Going back to the definitions of $v_j$, $d_j$ and $a_j$, and using the estimate \eqref{estimateh3c} for $d_1=v_3$ we find that 
\mathbf{e}gin{equation} 
  \label{hkplus2regularityforvfullalt}
  \mathbf{e}gin{aligned}
   \|D^{k+2} & \mathbf{v}\|_{L^2(\Omega)}^2 \leq \frac{c_k(\Omega)}{\alpha^2} \left( \|\mathbf{b}f\|_{H^k_{\lambda_1}(\Omega)}^2 \right. \\
  & \left. + \sum_{j=2}^{k-1} \left( \frac{1}{\alpha^{2(k-j)}}\frac{1}{\lambda_1^{(k-j)/2}}\|\mathsf{A}\|_{H^3_{\lambda_1}(\Omega)}^{2(k-j)} +  \frac{1}{\alpha^2}\frac{1}{\lambda_1^{1/2}}\|\mathsf{A}\|_{H^{k-j}_{\lambda_1}(\Omega)}^2\right)  \|\mathbf{b}f\|_{H^j_{\lambda_1}(\Omega)}^2 \right. \\
   & + \left( \frac{1}{\alpha^{2(k-1)}}\frac{1}{\lambda_1^{(k-1)/2}}\|\mathsf{A}\|_{H^3_{\lambda_1}(\Omega)}^{2(k-1)} +  \frac{1}{\alpha^2}\frac{1}{\lambda_1^{1/2}}\|\mathsf{A}\|_{H^{k-1}_{\lambda_1}(\Omega)}^2\right)\left( \|\mathbf{b}f\|_{H^1_{\lambda_1}(\Omega)}^2 \right. \\
  & \left. \quad \quad + \frac{1}{\alpha^2}\left( \|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}^2 + \|D^2\mathsf{A}\|_{L^3(\Omega)}^2\right)\left(\|\mathbf{b}f\|_{L^2(\Omega)}^2 + \frac{1}{\alpha^2}\|\mathsf{A}\|_{W^{1,\infty}_{\lambda_1}(\Omega)}^2 \|\mathbf{b}f\|_{H^{-1}(\Omega)}^2\right)\right),
  \end{aligned}
\end{equation}  
for suitable non-dimensional shape constants $c_k(\Omega)$, for any $k\geq 2$. A similar estimate holds for the pressure.
 \end{rmk}
 
\mathbf{e}gin{rmk}  
   In the case $k\geq 2$, in Theorem \ref{hkplustworegularity}, solving explicitly the recurrence relation in the estimate \eqref{hkplus2regularityforv} for $D^{k+2}\mathbf{v}$ is not a trivial matter. This relation is a so-called non-homogeneous linear recurrence relation and its solution depends on finding the roots of an intricate polynomial (see \cite{GreeneKnuth1982}). Luckily, we have found the exact solution for the recurrence relation and managed to estimate it in a suitable way (if \eqref{vrecursion} were an equality, than the equality in \eqref{recursioninequalitycomb} would be an exact solution).
\end{rmk}

\mathbf{e}gin{rmk}
  \label{dontneedsomuchregularityonA}
  As discussed in the paragraph before Theorem \ref{hkplustworegularity}, all the works in the literature assume more regularity on the coefficients of the linear term of the equation than we need in our theorem. One can see the reason for that by perusing the proof of \cite[Theorem 15.1]{friedman1969} (see also \cite[Theorem 17.2]{friedman1969}). In \cite[Theorem 15.1]{friedman1969}, an equation of the form $\sum_{|\alpha|=m} a_\alpha D^\alpha u = f$ is considered, and the assumed regularity on $f$ and the coefficients $a_\alpha$ guarantee that the term $f_i = \partial_{x_i} f - A_i$ belongs to $H^{k-1}$, with $A_i$ defined in \cite[eq. (15.22)]{friedman1969} ($A_i$ contains derivatives of $a_\alpha$ and $u$, analogous to our term $\mathbf{b}_\mathbf{g}amma$, in \eqref{bgammadef}). Since $A_i$ contains terms that are products of the first-order derivatives of the solution $u$ with second-order derivatives of the coefficients of the linear term and of second-order derivatives of the solution $u$ with first-order derivatives of the coefficients of the linear term, the assumptions that $u\in H^{k+1}$ and that the coefficients are in ${\mathcal C}^{k+1}$ guarantee that $f_i\in H^{k-1}$, in any dimensions. This, in turn, is used to bootstrap and show that $\partial_{x_i}u\in H^k$, hence $u\in H^{k+2}$. In our case, however, we work specifically in the three-dimensional case, so we can exploit some useful Sobolev embeddings, such as $H^{1/2}$ into $L^3$, $H^1$ into $L^6$, and $[H^1,H^2]_{1/2}$ into $L^\infty$. In doing so, we find that $f_i\in H^{k-1}$ even with the coefficients in $W^{1,\infty}\cap H^{k+1}$ (not necessarily in $W^{k+1,\infty}$). This is precisely what is done in the proof of Theorem \ref{hkplustworegularity}. Hence, we improve the classical result for any integer $k\geq 1$, in space dimension three. Analogous improvements can also be obtained in other dimensions.
\end{rmk}

\section{The viscoelastic Stokes problem}
\label{viscoelasticsec}
In this section, we apply the results of Section \ref{secgenstokes} to the viscoelastic Stokes problem \eqref{viscoelasticStokesorg} (or \eqref{viscoelasticStokes}). We first verify the conditions on $\mathsf{B}$ that guarantee that $\mathsf{A} = {\mathcal A}(\mathsf{B})$, given according to \eqref{defAcal}, is a uniformly positive definite symmetric tensor, in the sense of Definition \ref{defuniformpositivedef}. Next we investigate under which regularity assumptions on $\mathsf{B}$ the tensor $\mathsf{A}={\mathcal A}(\mathsf{B})$ has the regularity conditions required in the results of Section \ref{secgenstokes}. Next we combine all these results to obtain the desired existence, uniqueness, and regularity properties of the viscoelastic Stokes problem.

\subsection{Conditions for ellipticity of the linear part of the viscoelastic Stokes problem}
\label{conditionsellipticity}

In applications, $\mathsf{B}$ is the left Cauchy-Green tensor $\mathsf{F}\mathsf{F}^\text{\rm tr}$, so we restrict our analysis to the case in which $\mathsf{B}$ is a positive definite symmetric tensor. We want to find conditions on $\mu_1$, $\mu_2$, $\mu_3$ and $\mathsf{B}$ for which ${\mathcal A}(\mathsf{B})$, given by \eqref{defAcal}, is also a positive define symmetric tensor, which is one of the conditions needed in Theorems \ref{weaksolutionweakgeneralizedStokeswithpressure}, \ref{h2regularity} and \ref{hkplustworegularity}. It is clear that ${\mathcal A}(\mathsf{B})$ is also symmetric, so we only need to investigate when it is positive definite.

In order to study this issue, it is natural to consider the real-valued functions
\mathbf{e}gin{equation}
  \label{defg}
  g_\mathbf{x}(\lambda) = \mu_1(\mathbf{x}) + \mu_2(\mathbf{x})\lambda + \frac{\mu_3(\mathbf{x})}{\lambda}, \quad \text{for } \lambda > 0, \; \mathbf{x}\in \Omega.
\end{equation}
Notice that if $\lambda$ is an eigenvalue of $\mathsf{B}(\mathbf{x})$, then $g_\mathbf{x}(\lambda)$ is an eigenvalue of ${\mathcal A}(\mathsf{B}(\mathbf{x}))$, with the same eigenspace. Hence, the positivity of ${\mathcal A}(\mathsf{B}(\mathbf{x}))$ is directly related to the positivity of $g_\mathbf{x}(\lambda)$, for every eigenvalue of $\mathsf{B}(\mathbf{x})$. For the sake of notational simplicity, given a symmetric matrix $\mathsf{B}$ in three dimensions, we denote by $\lambda_1(\mathsf{B}) \leq \lambda_2(\mathsf{B}) \leq \lambda_3(\mathsf{B})$ the three eigenvalues of $\mathsf{B}$, in increasing order. Then, it is straightforward to deduce that the uniform positive definiteness of the tensor ${\mathcal A}(\mathsf{B})$, in the sense of Definition \ref{defuniformpositivedef}, follows if, and only if, $g_\mathbf{x}(\lambda_i(\mathsf{B}(\mathbf{x})))$ is strictly positive uniformly on $\Omega$, for every $i=1,2,3$.

We summarize the above facts in the following result. 
\mathbf{e}gin{prop}
  \label{propessinfglambdacondition}
  Let $\mathsf{B}=\mathsf{B}(\mathbf{x})$ be a tensor on $\Omega$ and suppose that, for almost every $\mathbf{x}$ in $\Omega$, $\mathsf{B}(\mathbf{x})$ is symmetric and positive definite with eigenvalues $\lambda_i(\mathsf{B}(\mathbf{x}))>0$, $i=1,2,3$. Then, ${\mathcal A}(\mathsf{B})$ is a symmetric tensor and the eigenvalues of ${\mathcal A}(\mathsf{B}(\mathbf{x}))$ are given by $g_\mathbf{x}(\lambda_i(\mathsf{B}(\mathbf{x}))$, $i=1,2,3$, with the same corresponding eigenspaces of $\mathsf{B}(\mathbf{x})$. Moreover, $\mathsf{B}$ is ${\mathcal A}$-positive, in the sense of Definition \ref{defproper}, if, and only if, 
  \mathbf{e}gin{equation}
    \label{essinfglambda}
    \alpha \define \operatorname*{ess\,inf}_{\mathbf{x}\in \Omega} \left\{\min_{i=1,2,3} \left(\mu_1(\mathbf{x}) + \mu_2(\mathbf{x})\lambda_i(\mathsf{B}(\mathbf{x})) + \frac{\mu_3(\mathbf{x})}{\lambda_i(\mathsf{B}(\mathbf{x}))} \right)\right\} > 0.
  \end{equation}
  In this case, \eqref{unifellip} holds with $\mathsf{A} = {\mathcal A}(\mathsf{B})$ and $\alpha$ as given in \eqref{essinfglambda}.
\end{prop}

Typically, the parameters $\mu_1$, $\mu_2$, and $\mu_3$ may depend on some invariants of the left Cauchy-Green stress tensor $\mathsf{B}$ and/or on the temperature and other quantities, when such couplings are relevant. In some models, however, they are taken as constants. In the cases they are constant, we investigate in more details when the condition \eqref{essinfglambda} is satisfied. 

When these parameters do not change with $\mathbf{x}$, the function in \eqref{defg} is the same for every $\mathbf{x}$, i.e. $g=g_\mathbf{x}$. Since $\mathsf{B}$ is a positive definite symmetric tensor, all its eigenvalues are positive, so, given a set of real-valued parameters $\mu_1, \mu_2, \mu_3$, we look for set of positive values of $\lambda$ for which $g$ is positive, i.e.
  \mathbf{e}gin{equation}
    \label{defLambda}
     \Lambda_{\mu_1,\mu_2,\mu_3} = \left\{\lambda > 0; \;g(\lambda) = \mu_1 + \mu_2\lambda + \frac{\mu_3}{\lambda} > 0\right\}.
  \end{equation}

Since $g$ is a continuous functions, the following result is easy to obtain and its proof is omitted.
  
\mathbf{e}gin{prop}
  \label{propmuscenarios}
  Let $\mu_1, \mu_2, \mu_3\in {\mathbb{R}}$ be constants and let $\mathsf{B}$ be a positive definite symmetric tensor on a domain $\Omega\subset{\mathbb{R}}^3$. Suppose there exists a set $\Lambda_0$ which is compactly included in the set $\Lambda_{\mu_1,\mu_2,\mu_3}$ defined in \eqref{defLambda} and such that $\lambda_i(\mathsf{B}(\mathbf{x}))\in \Lambda_0$ almost everywhere on $\Omega$, for every $i=1, 2, 3$. Then, condition \eqref{essinfglambda} of Proposition \ref{propessinfglambdacondition} is fulfilled, and $\mathsf{B}$ is ${\mathcal A}$-positive in the sense of Definition \ref{defproper}.
\end{prop}

  The set  $\Lambda_{\mu_1,\mu_2,\mu_3}$ can be explicitly and easily characterized in terms of the parameters. For that purpose, it is useful to introduce the polynomial $p(\lambda) = \lambda g(\lambda) = \mu_2\lambda^2 + \mu_1\lambda +\mu_3$, which has the same sign as $g$ for $\lambda$ positive. The roots of $p$ are real when $\mu_1^2 \geq 4\mu_2\mu_3$ and are given by
  \mathbf{e}gin{equation}
    \label{proots}
    \lambda_{\pm} = \frac{-\mu_1 \pm \sqrt{\mu_1^2 - 4\mu_2\mu_3}}{2\mu_2}.
  \end{equation}
  
   One fundamental condition on the parameters comes when considering the state of the material with no displacement, in which case $\mathsf{F}$ is the identity tensor, and so is $\mathsf{B}$. Hence, all the eigenvalues of $\mathsf{B}(\mathbf{x})$ are equal to one, with $g(1)= \mu_1+\mu_2+\mu_3$. Thus, we recover the classic thermodynamic condition
  \mathbf{e}gin{equation}
    \label{sumofmuspositive}
    \mu_1 + \mu_2 + \mu_3 > 0.
  \end{equation}
  Assuming condition \eqref{sumofmuspositive}, there are a few different scenarios that may happen, as follows.
  
  \mathbf{e}gin{lem}
    \label{lemmuscenarios}
    Let $\mu_1, \mu_2, \mu_3\in {\mathbb{R}}$ be constants, and consider the function $g$ and the set $\Lambda_{\mu_1,\mu_2,\mu_3}$ given in \eqref{defLambda}. Suppose that \eqref{sumofmuspositive} holds. Then, we have the following scenarios.
  \renewcommand{(\roman{enumi})}{(\roman{enumi})}
  \mathbf{e}gin{enumerate}
    \item If $\mu_2, \mu_3> 0$ and $\mu_1> -2\sqrt{\mu_2\mu_3}$, then $\Lambda_{\mu_1,\mu_2,\mu_3} = (0,\infty)$.
    \item If $\mu_2, \mu_3> 0$ and $\mu_1\leq -2\sqrt{\mu_2\mu_3}$, then $0 <\lambda_- \leq \lambda_+$, and $\Lambda_{\mu_1,\mu_2,\mu_3} = (0,\lambda_-)\cup (\lambda_+,\infty)$.
    \item If $\mu_3>0$, $\mu_2=0$, and $\mu_1\geq 0$, then $\Lambda_{\mu_1,\mu_2,\mu_3} = (0,\infty)$.
    \item If $\mu_3>0$, $\mu_2=0$, and $\mu_1<0$, then $\lambda_0 = -\mu_3/\mu_1>0$ and $\Lambda_{\mu_1,\mu_2,\mu_3} = (0, \lambda_0)$.
    \item If $\mu_3>0$ and $\mu_2<0$, then $\Lambda_{\mu_1,\mu_2,\mu_3} = (0, \lambda_+)$.
    \item If $\mu_3=0$ and $\mu_1,\mu_2\geq 0$, then $\Lambda_{\mu_1,\mu_2,\mu_3} = (0,\infty)$.
    \item If $\mu_3=0$, $\mu_2>0$, and $\mu_1<0$, then $\Lambda_{\mu_1,\mu_2,\mu_3} = (\lambda_+,\infty)$.
    \item If $\mu_3=0$, $\mu_2<0$, and $\mu_1>0$, then $\Lambda_{\mu_1,\mu_2,\mu_3} = (0,\lambda_+)$.
    \item If $\mu_3<0$ and $\mu_2>0$, then $\Lambda_{\mu_1,\mu_2,\mu_3} = (\lambda_+,\infty)$.
    \item If $\mu_3,\mu_2<0$ and $\mu_1> 2\sqrt{\mu_2\mu_3}$, then $\Lambda_{\mu_1,\mu_2,\mu_3} = (\lambda_-,\lambda_+)$.
    \item If $\mu_3<0$, $\mu_2=0$, and $\mu_1>0$, then $\lambda_0 = -\mu_3/\mu_1>0$ and $\Lambda_{\mu_1,\mu_2,\mu_3} = (\lambda_0,\infty)$.
  \end{enumerate}
\end{lem}

We end this section with a few remarks.
\mathbf{e}gin{rmk}[Small Deformations]
  \label{rmksmalldeformations}
  Without any deformation, the tensor $\mathsf{F}$ is the identity $\mathsf{I}$, and so is the left Cauchy-Green tensor $\mathsf{B}=\mathsf{F}\mathsf{F}^\text{\rm tr} = \mathsf{I}$. In this case, $\mathsf{A}={\mathcal A}(\mathsf{I})$ is a diagonal tensor, with identical diagonal elements given by $\mu_1+\mu_2+\mu_3$. Under the assumption \eqref{sumofmuspositive}, $\mathsf{A}$ is uniformly elliptic.  Then, by continuity, it is clear that, for small perturbations, the eigenvalues of $\mathsf{B}$ are near the identity and, hence, the eigenvalues of $\mathsf{A}={\mathcal A}(\mathsf{B})$ are still strictly positive. More precisely, one has that there exists $\delta>0$ such that if $\mathsf{B}\in L^\infty(\Omega)^{3\times 3}$ is a positive definite symmetric tensor such that $\|\mathsf{B}-\mathsf{I}\|_{L^\infty(\Omega)} \leq \delta$, then $\mathsf{B}$ is still an ${\mathcal A}$-positive tensor. The value of $\delta$ can be estimated in terms of $\mu_1,\mu_2,\mu_3$ depending on the case described in Proposition \ref{propmuscenarios}. In particular, the material may, or may not, allow arbitrarily large deformations.
\end{rmk}

\mathbf{e}gin{rmk}[Salt dynamics]
In the time-dependent problem associated with the constitutive law \eqref{sigmaconstitutive}, we consider an Eulerian formulation in which $\mathsf{B}$ is a time-dependent variable of the system. In the initial value problem for this system, without any initial deformation, $\mathsf{B}$ is taken to be the identity tensor at the initial time. In this case, the associate tensor $\mathsf{A}={\mathcal A}(\mathsf{B})$ is initially a diagonal operator with diagonal element equal to $\mu_1+\mu_2+\mu_3>0$ (see Remark \ref{rmksmalldeformations}). Using energy estimates of the associated evolution equations and the regularization property of problem \eqref{generalizedStokes} proved here, we obtain a~priori estimates for $\mathsf{B}$ in $H^3(\Omega)^{3\times 3}\subset W^{1,\infty}(\Omega)$ and for $\mathbf{v}$ in $H^4(\Omega)^3$, showing that $\mathsf{B}$ remains in $H^3(\Omega)^{3\times 3}$ and stays ${\mathcal A}$-positive for a short time, while $\mathbf{v}$ remains in $H^4(\Omega)^3$. In doing so, we prove that the system is well-posed locally in time. This result will be presented elsewhere.
\end{rmk}

\subsection{Regularity of ${\mathcal A}(\mathsf{B})$}
  \label{functionalestimatesforAcalofB}

In view of the regularity conditions on $\mathsf{A}$ needed for the existence, uniqueness, and regularity results given in Section \ref{secgenstokes}, we want to obtain the corresponding regularity properties of $\mathsf{A}={\mathcal A}(\mathsf{B})$ in terms of $\mathsf{B}$. This is the aim of the current section. Along the way, some useful estimates are derived.

One step is to represent the inverse $\mathsf{B}^{-1}$ in terms of $\mathsf{B}^2$, using the representation 
\[ p_\mathsf{B}(\lambda) = \lambda^3 - \rm{I}_\mathsf{B} \lambda^2 + \rm{II}_\mathsf{B} \lambda - \rm{III}_\mathsf{B} 
\]
for the characteristic polynomial of a tridimensional matrix $\mathsf{B}(\mathbf{x})$, at each point $\mathbf{x}$. In this representation,
\mathbf{e}gin{align*} 
  \rm{I}_\mathsf{B} & = \operatorname{Tr} \mathsf{B}, \\
  \rm{II}_\mathsf{B} & =  \frac{1}{2}\left( (\operatorname{Tr} \mathsf{B})^2 - \operatorname{Tr}(\mathsf{B}^2)\right), \\
  \rm{III}_\mathsf{B} & = \frac{1}{6}\left((\operatorname{Tr} \mathsf{B})^3 - 3(\operatorname{Tr} \mathsf{B}) (\operatorname{Tr}(\mathsf{B}^2)) + 2 \operatorname{Tr}(\mathsf{B}^3)\right) = \det \mathsf{B},
\end{align*}
where $\operatorname{Tr}(\cdot)$ denotes the trace of a matrix or tensor, and $\det(\cdot)$, its determinant. We use the Cayley-Hamilton Theorem, which asserts that
\mathbf{e}gin{equation}
  \label{cayleyhamilton} 
  p_\mathsf{B}(\mathsf{B}) = \mathsf{B}^3 - \rm{I}_\mathsf{B} \mathsf{B}^2 + \rm{II}_\mathsf{B} \mathsf{B} - \rm{III}_\mathsf{B} \mathsf{I} = 0, \quad \text{ on } \Omega.
\end{equation}
Hence,
  \mathbf{e}gin{equation} 
    \label{cayleyhamiltonbinv}
    \mathsf{B}^{-1} = \frac{1}{\det \mathsf{B}}(\mathsf{B}^2 - \rm{I}_\mathsf{B} \mathsf{B} + \rm{II}_\mathsf{B} \mathsf{I}).
  \end{equation}
  
  Using this representation, we first obtain an $L^\infty$ regularity result for ${\mathcal A}(\mathsf{B})$.  
\mathbf{e}gin{prop}
  \label{proplinftyboundonacalb}
  If $\mathsf{B}\in L^\infty(\Omega)^{3\times 3}$ and $(\det B)^{-1}\in L^\infty(\Omega)$, then, $\mathsf{B}^{-1}\in L^\infty(\Omega)^{3\times 3}$. If, moreover, $\mu_1, \mu_2, \mu_3\in L^\infty(\Omega)$, then ${\mathcal A}(\mathsf{B})\in L^\infty(\Omega)^{3\times 3}$.
\end{prop}

\mathbf{e}gin{proof}
Clearly, $\mathsf{B}^2$, $\operatorname{Tr} \mathsf{B}$, and $\operatorname{Tr} \mathsf{B}^2$ are essentially bounded, with 
\[  \|\mathsf{B}^2\|_{L^\infty(\Omega)} \leq 3\|\mathsf{B}\|_{L^\infty(\Omega)}^2, \quad \|\operatorname{Tr}\mathsf{B}\|_{L^\infty(\Omega)} \leq 3\|\mathsf{B}\|_{L^\infty(\Omega)}, \quad \|\operatorname{Tr}\mathsf{B}^2\|_{L^\infty(\Omega)} \leq 9\|\mathsf{B}\|_{L^\infty(\Omega)}^2,
\]
and, hence,
\[ \|\rm{I}_{\mathsf{B}}\mathsf{B}\|_{L^\infty(\Omega)} \leq 3\|\mathsf{B}\|_{L^\infty(\Omega)}^2, \quad \|\rm{II}_{\mathsf{B}}\mathsf{I}\|_{L^\infty(\Omega)} \leq 9\|\mathsf{B}\|_{L^\infty(\Omega)}^2.
\]
Using these estimates in \eqref{cayleyhamiltonbinv} implies that $\mathsf{B}^{-1}\in L^\infty(\Omega)$, with
  \mathbf{e}gin{equation}
    \label{estimateBinvlinfty}
    \|\mathsf{B}^{-1}\|_{L^\infty(\Omega)} \leq 15\|(\det\mathsf{B})^{-1}\|_{L^\infty(\Omega)} \|\mathsf{B}\|_{L^\infty(\Omega)}^2.
  \end{equation}
  Then, using \eqref{estimateBinvlinfty} in \eqref{defAcal}, we obtain
    \mathbf{e}gin{equation}
    \label{estimateAcalBlinfty}
    \|{\mathcal A}(\mathsf{B})\|_{L^\infty(\Omega)} \leq \|\mu_1\|_{L^\infty(\Omega)} + \|\mu_2\|_{L^\infty(\Omega)}\|\mathsf{B}\|_{L^\infty(\Omega)} + 15 \|\mu_3\|_{L^\infty(\Omega)}\|(\det\mathsf{B})^{-1}\|_{L^\infty(\Omega)}\|\mathsf{B}\|_{L^\infty(\Omega)}^2,
  \end{equation}
  which proves that ${\mathcal A}(\mathsf{B})\in L^\infty(\Omega)^{3\times 3}$.
\end{proof}

For the sake of simplicity and with the applications in mind, we restrict the next estimates to the relevant case $\det\mathsf{B}=1$. We then look for $W^{1,\infty}(\Omega)$ regularity.
\mathbf{e}gin{prop}
  \label{propwoneinftyboundonacalb}
  If $\mathsf{B}\in W^{1,\infty}(\Omega)^{3\times 3}$ and $\det B = 1$ almost everywhere on $\Omega$, then, $\mathsf{B}^{-1}\in W^{1,\infty}(\Omega)^{3\times 3}$. If, moreover, $\mu_1, \mu_2, \mu_3\in W^{1,\infty}(\Omega)$, then ${\mathcal A}(\mathsf{B})\in W^{1,\infty}(\Omega)^{3\times 3}$.
\end{prop}

\mathbf{e}gin{proof}
That $\mathsf{B}^{-1}$ and ${\mathcal A}(\mathsf{B})$ belong to $L^\infty(\Omega)^{3\times 3}$ follows from Proposition \ref{proplinftyboundonacalb}. We only need to consider the derivatives of such tensors.

Since the trace is a linear operator and taking into account that $\det\mathsf{B}=1$, the derivatives $\partial_{x_i}\left(\mathsf{B}^{-1}\right)$ are easily found from \eqref{cayleyhamiltonbinv} to be
\mathbf{e}gin{equation}
  \label{partialxiofinvB}
  \partial_{x_i} \left(\mathsf{B}^{-1}\right) = 2\mathsf{B}\partial_{x_i}\mathsf{B} - \operatorname{Tr}(\partial_{x_i}\mathsf{B}) \mathsf{B} + \operatorname{Tr}(\mathsf{B})\partial_{x_i}\mathsf{B} + \left( \operatorname{Tr}(\mathsf{B})\operatorname{Tr}(\partial_{x_i}\mathsf{B}) -  \operatorname{Tr}(\mathsf{B}\partial_{x_i}\mathsf{B})\right) \mathsf{I}.
\end{equation}
for every $i=1,2,3.$ Thus,
\[ \|\partial_{x_i} \left(\mathsf{B}^{-1}\right)\|_{L^\infty(\Omega)} \leq 20\|\mathsf{B}\|_{L^\infty(\Omega)}\|\partial_{x_i}\mathsf{B}\|_{L^\infty(\Omega)},
\]
This yields 
  \mathbf{e}gin{equation}
    \label{estimateBinvwoneinfty}
    \|D(\mathsf{B}^{-1})\|_{L^\infty(\Omega)} \leq 20\|\mathsf{B}\|_{L^\infty(\Omega)} \|D\mathsf{B}\|_{L^\infty(\Omega)}.
  \end{equation}
From \eqref{defAcal} we find that
\mathbf{e}gin{multline*}
  \|\partial_{x_i}{\mathcal A}(\mathsf{B})\|_{L^\infty(\Omega)} \leq \|\partial_{x_i}\mu_1\|_{L^\infty(\Omega)} + \|\partial_{x_i}\mu_2\|_{L^\infty(\Omega)}\|\mathsf{B}\|_{L^\infty(\Omega)} + \|\mu_2\|_{L^\infty(\Omega)} \|\partial_{x_i}\mathsf{B}\|_{L^\infty(\Omega)} \\
  + \|\partial_{x_i}\mu_3\|_{L^\infty(\Omega)}\|\mathsf{B}^{-1}\|_{L^\infty(\Omega)} + \|\mu_3\|_{L^\infty(\Omega)}\|\partial_{x_i}(\mathsf{B}^{-1})\|_{L^\infty(\Omega)}.
\end{multline*}
which, together with \eqref{estimateBinvlinfty} and \eqref{estimateBinvwoneinfty}, yields   
\mathbf{e}gin{multline}
    \label{estimateAcalBwoneinfty}
    \|D{\mathcal A}(\mathsf{B})\|_{L^\infty(\Omega)} \leq \|D\mu_1\|_{L^\infty(\Omega)} + \|D\mu_2\|_{L^\infty(\Omega)}\|\mathsf{B}\|_{L^\infty(\Omega)} + \|\mu_2\|_{L^\infty(\Omega)}\|D\mathsf{B}\|_{L^\infty(\Omega)} \\
    + 9 \|D\mu_3\|_{L^\infty(\Omega)}\|\mathsf{B}\|_{L^\infty(\Omega)}^2 +20\|\mu_3\|_{L^\infty(\Omega)}\|\mathsf{B}\|_{L^\infty(\Omega)}\|D\mathsf{B}\|_{L^\infty(\Omega)}.
  \end{multline}
  This proves the desired regularity for ${\mathcal A}(\mathsf{B})$.
\end{proof}

Now we look for a $W^{2,3}(\Omega)$ estimate.
\mathbf{e}gin{prop}
  \label{propw23boundonacalb}
  If $\mathsf{B}\in W^{2,3}(\Omega)^{3\times 3}$ and $\det B = 1$ almost everywhere on $\Omega$, then, $\mathsf{B}^{-1}\in W^{2,3}(\Omega)^{3\times 3}$. If, moreover, $\mu_1, \mu_2, \mu_3\in W^{1,\infty}(\Omega)\cap W^{2,3}(\Omega)$, then ${\mathcal A}(\mathsf{B})\in W^{2,3}(\Omega)^{3\times 3}$.
\end{prop}

\mathbf{e}gin{proof}
From \eqref{cayleyhamiltonbinv}, and using that $\det \mathsf{B} =1$ almost everywhere, we see that
\[ \|\mathsf{B}^{-1}\|_{L^3(\Omega)} \leq c \|\mathsf{B}\|_{L^6(\Omega)}^2,
\] 
for some universal constant $c$. Using \eqref{partialxiofinvB}, we find that
\[ \|D\mathsf{B}^{-1}\|_{L^3(\Omega)} \leq c \|\mathsf{B}\|_{L^6(\Omega)}\|D\mathsf{B}\|_{L^6(\Omega)},
\]
for a different universal constant $c$. Taking the $\partial_{x_j}$ derivative of \eqref{partialxiofinvB} we find that
\mathbf{e}gin{multline}
  \label{partialxixjofinvB}
  \partial_{x_i x_j}^2 \left(\mathsf{B}^{-1}\right) = 2\partial_{x_j}\mathsf{B}\partial_{x_i}\mathsf{B} + 2\mathsf{B}\partial_{x_i x_j}^2\mathsf{B} \\
  - \operatorname{Tr}(\partial_{x_ix_j}^2\mathsf{B}) \mathsf{B} - \operatorname{Tr}(\partial_{x_i}\mathsf{B}) \partial_{x_j} \mathsf{B} + \operatorname{Tr}(\partial_{x_j}\mathsf{B})\partial_{x_i}\mathsf{B} + \operatorname{Tr}(\mathsf{B})\partial_{x_i x_j}^2\mathsf{B} \\
  + \left( \operatorname{Tr}(\partial_{x_j}\mathsf{B})\operatorname{Tr}(\partial_{x_i}\mathsf{B}) + \operatorname{Tr}(\mathsf{B})\operatorname{Tr}(\partial_{x_i x_j}^2\mathsf{B}) -  \operatorname{Tr}(\partial_{x_j}\mathsf{B}\partial_{x_i}\mathsf{B}) -  \operatorname{Tr}(\mathsf{B}\partial_{x_i x_j}^2\mathsf{B})\right) \mathsf{I},
\end{multline}
for every $i=1,2,3.$ Then, using H\"older's inequality, we find that
  \mathbf{e}gin{equation}
    \label{estimateBinvw23}
    \|D^2(\mathsf{B}^{-1})\|_{L^3(\Omega)} \leq c \left(\|D\mathsf{B}\|_{L^6(\Omega)}^2 +  \|\mathsf{B}\|_{L^\infty(\Omega)}\|D^2\mathsf{B}\|_{L^3(\Omega)}\right),
  \end{equation}
  for another suitable universal constant $c$. These estimates show that $\mathsf{B}^{-1}$ belongs to $W^{2,3}(\Omega)^{3\times 3}$. Now, from \eqref{defAcal} and the regularity obtained for $\mathsf{B}^{-1}$, we readily see that ${\mathcal A}(\mathsf{B})\in W^{2,3}(\Omega)^{3\times 3}$, with, in particular,
  \mathbf{e}gin{multline}
    \label{estimateAcalBw23}
    \|D^2{\mathcal A}(\mathsf{B})\|_{L^3(\Omega)} \leq c\left(\|D^2\mu_1\|_{L^3(\Omega)} 
      + \|D^2\mu_2\|_{L^3(\Omega)}\|\mathsf{B}\|_{L^\infty(\Omega)} 
      + \|D\mu_2\|_{L^\infty(\Omega)}\|D\mathsf{B}\|_{L^3(\Omega)} \right. \\ + \|\mu_2\|_{L^\infty(\Omega)}\|D\mathsf{B}\|_{L^3(\Omega)} 
+ \|D^2\mu_3\|_{L^3(\Omega)}\|\mathsf{B}^{-1}\|_{L^\infty(\Omega)} \\ \left. + \|D\mu_3\|_{L^\infty(\Omega)}\|D(\mathsf{B}^{-1})\|_{L^3(\Omega)} + \|\mu_3\|_{L^\infty(\Omega)}\|D^2(\mathsf{B}^{-1})\|_{L^3(\Omega)}\right),
  \end{multline}
  for a suitable universal constant $c$.
\end{proof}

Finally we look for $H^m(\Omega)$ estimates, for any integer $m\geq 2$.
\mathbf{e}gin{prop}
\label{prophmboundonacalb}
If $\mathsf{B}\in H^m(\Omega)$ and $\det B = 1$ almost everywhere on $\Omega$, where $m\in {\mathbb{N}}$, $m\geq 2$, then, $\mathsf{B}^{-1}\in H^m(\Omega)^{3\times 3}$. If, moreover, $\mu_1, \mu_2, \mu_3\in W^{1,\infty}(\Omega)
\cap H^m(\Omega)^{m\times m}$, then ${\mathcal A}(\mathsf{B})\in H^m(\Omega)^{3\times 3}$.
\end{prop}

\mathbf{e}gin{proof}
Since $\mathsf{B}\in H^m(\Omega)^{3\times 3} \subset L^4(\Omega)^{3\times 3}$, for $m\geq 2$, it is straightforward to see that
  \mathbf{e}gin{equation}
    \label{estimateBinvl2}
    \|\mathsf{B}^{-1}\|_{L^2(\Omega)} \leq c \|\mathsf{B}\|_{L^4(\Omega)}^2,
  \end{equation}
for a suitable universal constant $c$, so that $\mathsf{B}^{-1}$ is in $L^2(\Omega)^{3\times 3}$. Using H\"older's inequality in \eqref{partialxiofinvB} we see that
  \mathbf{e}gin{equation}
    \label{estimateBinvh1}
    \|D(\mathsf{B}^{-1})\|_{L^2(\Omega)} \leq c\|\mathsf{B}\|_{L^{12}(\Omega)}\|D\mathsf{B}\|_{L^{12/5}(\Omega)}^2,
  \end{equation}
for another universal constant $c$. Since $H^m(\Omega)\subset L^{12}(\Omega)\cap W^{1,12/5}(\Omega)$, this shows that $\mathsf{B}^{-1}\in H^1(\Omega)^{3\times 3}$. From \eqref{partialxixjofinvB}, we see that
\mathbf{e}gin{equation}
  \label{estimateBinvh2}
  \|D^2(\mathsf{B}^{-1})\|_{L^2(\Omega)} \leq c \left(\|D\mathsf{B}\|_{L^4(\Omega)}^2 +  \|\mathsf{B}\|_{L^\infty(\Omega)}\|D^2\mathsf{B}\|_{L^2(\Omega)}\right),
  \end{equation}
  for another suitable universal constant $c$. Since in three dimensions we have the embedding $H^m\subset L^\infty$, for $m\geq 2$, it follows that $\mathsf{B}^{-1}\in H^2(\Omega)^{3\times 3}$ as well. This proves the case $m=2$. For $m = 3, 4, \ldots$, we use Leibniz's rule and estimate $D^m(\mathsf{B}^{-1})$ simply by
  \[ \|D^m(\mathsf{B}^{-1})\|_{L^2(\Omega)} \leq c\sum_{i=1}^{[m/2]} \|D^i\mathsf{B}\|_{L^\infty(\Omega)} \|D^{m-i}\mathsf{B}\|_{L^2(\Omega)},
  \]
  where $[m/2]$ is the largest integer smaller than or equal to $m/2$. When $m\geq 3$ and the space is three-dimensional, the embedding $H^m(\Omega)\subset W^{i,\infty}$ holds for any integer $i\leq m/2$. Thus, $\mathsf{B}^{-1}\in H^m(\Omega)^{3\times 3}$.
  
  For the regularity of ${\mathcal A}(\mathsf{B})$, we see that
  \mathbf{e}gin{equation}
    \label{estimateAcalBinl2}
    \|{\mathcal A}(\mathsf{B})\|_{L^2(\Omega)} \leq c\left(\|\mu_1\|_{L^2(\Omega)} + \|\mu_2\|_{L^\infty(\Omega)}\|\mathsf{B}\|_{L^2(\Omega)} + \|\mu_3\|_{L^\infty(\Omega)}\|\mathsf{B}^{-1}\|_{L^2(\Omega)}\right),
  \end{equation}
  for a universal constant $c$, and
  \mathbf{e}gin{multline}
    \label{estimateAcalBinhm}
    \|D^m{\mathcal A}(\mathsf{B})\|_{L^2(\Omega)} \leq c\left(\|D^m\mu_1\|_{L^2(\Omega)} + \|\mu_2\|_{L^\infty(\Omega)}\|D^m\mathsf{B}\|_{L^2(\Omega)} + \|D\mu_2\|_{L^\infty(\Omega)}\|D^{m-1}\mathsf{B}\|_{L^2(\Omega)} \right. \\
      + \|\mu_3\|_{L^\infty(\Omega)}\|D^m(\mathsf{B}^{-1})\|_{L^2(\Omega)} + \|D\mu_3\|_{L^\infty(\Omega)}\|D^{m-1}(\mathsf{B}^{-1})\|_{L^2(\Omega)} \\
      \left. + \sum_{i=2}^m \|D^i\mu_2\|_{L^2(\Omega)}\|D^{m-i}\mathsf{B}\|_{L^\infty(\Omega)} + \sum_{i=2}^m \|D^i\mu_3\|_{L^2(\Omega)}\|D^{m-i}(\mathsf{B}^{-1})\|_{L^\infty(\Omega)} \right),
  \end{multline}
  for another suitable constant $c$. From the embedding $H^m(\Omega) \subset W^{m-i,\infty}$, for $m\geq 2$ and $2\leq i \leq m$, valid in three dimensions, it follows from \eqref{estimateAcalBinl2} and \eqref{estimateAcalBinhm} that ${\mathcal A}(\mathsf{B})\in H^m(\Omega)$.
\end{proof}

\mathbf{e}gin{rmk}
  If the regularities of ${\mathcal A}(\mathsf{B})$ in $W^{2,3}(\Omega)^{3\times 3}$ and in $H^m(\Omega)^{3\times 3}$ were the sole purposes of Propositions \ref{propw23boundonacalb} and \ref{prophmboundonacalb}, respectively, then, of course, one could ask for different regularity assumptions on $\mathsf{B}$ and on the parameters $\mu_1, \mu_2$, and $\mu_3$. But since these estimates will be used in conjunction with Theorem \ref{hkplustworegularity}, which requires that $\mathsf{A} = {\mathcal A}(\mathsf{B})$ belong to $W^{1,\infty}(\Omega)^{3\times 3}$, then we are constrained by starting with the assumptions required in Proposition \ref{propwoneinftyboundonacalb}.
\end{rmk}

\mathbf{e}gin{rmk}
  More general conditions for the regularity of ${\mathcal A}(\mathsf{B})$ can be obtained by imposing conditions not only on $\mathsf{B}$ but also on $\mathsf{B}^{-1}$ itself. However, for practical purposes, since in the evolutionary system \eqref{LiueulerB} the tensor $\mathsf{B}$ is one of the unknowns, it is useful to derive conditions on $\mathsf{B}$ alone.
\end{rmk}

\subsection{Existence and regularity results for the viscoelastic Stokes problem}
\label{existuniqregviscoelasticsec}

The existence, uniqueness, and regularity properties for the viscoelastic Stokes problem \eqref{viscoelasticStokes} are direct consequences of the corresponding results for the generalized Stokes problem \eqref{generalizedStokes}, provided the tensor $\mathsf{B}$ is ${\mathcal A}$-positive, which guarantees that the corresponding tensor $\mathsf{A}={\mathcal A}(\mathsf{B})$ is uniformly positive definite, and provided the tensor ${\mathcal A}(\mathsf{B})$ is sufficiently regular. The main condition for $\mathsf{B}$ to be ${\mathcal A}$-positive is given in Proposition \ref{propessinfglambdacondition}, in Section \ref{conditionsellipticity}, while the regularity of ${\mathcal A}(\mathsf{B})$ is investigated in Propositions \ref{proplinftyboundonacalb}, \ref{propwoneinftyboundonacalb}, \ref{propw23boundonacalb}, and \ref{prophmboundonacalb}, in Section \ref{functionalestimatesforAcalofB}. The aim of the current section is to put these things together. With that in mind, we state, without proof, the following results, corresponding to Theorems \ref{weaksolutionweakgeneralizedStokeswithpressure}, \ref{h2regularity}, and \ref{hkplustworegularity}, respectively.

\mathbf{e}gin{thm}[Existence and uniqueness]
  \label{weaksolutionweakviscoelasticdStokeswithpressure}
  Let $\Omega\subset {\mathbb{R}}^3$ be a bounded Lipschitz domain, $\mathsf{B}\in L^\infty(\Omega)^{3\times 3}$ be a uniformly positive definite symmetric tensor on $\Omega$ with $(\det \mathsf{B})^{-1}\in L^\infty(\Omega)$, $\mu_1$, $\mu_2$, and $\mu_3 \in L^\infty(\Omega)$ be scalar fields, and $\mathbf{b}f\in H^{-1}(\Omega)^3$. Assume that \eqref{essinfglambda} holds. Then, there exist a unique vector field $\mathbf{v}\in V(\Omega)$ and a scalar field $p\in L^2(\Omega)$ that solve the weak formulation \eqref{weakviscoelasticStokeswithpressure} of the viscoelastic Stokes problem. Moreover, the field $p$ is unique up to a constant, and the estimates
  \mathbf{e}gin{equation}
    \label{weaksolutionweakgeneralizedStokesestimaterepeatwithB}
    \|\mathbf{n}abla \mathbf{v}\|_{L^2(\Omega)} \leq \frac{1}{\alpha} \|\mathbf{b}f\|_{H^{-1}(\Omega)},
  \end{equation}
  and  
  \mathbf{e}gin{equation}
    \label{weaksolutionweakgeneralizedStokesestimatewithpressurewithB}
    \|p\|_{L^2(\Omega)/{\mathbb{R}}} \leq \frac{c(\Omega)}{\alpha} \|{\mathcal A}(\mathsf{B})\|_{L^\infty(\Omega)}\|\mathbf{b}f\|_{H^{-1}(\Omega)},
  \end{equation}
  hold, where $c(\Omega)$ is a shape constant.
\end{thm}

\mathbf{e}gin{thm}[Regularity]
  \label{h2regularityforviscoelasticStokes}
  Under the hypotheses of Theorem \ref{weaksolutionweakviscoelasticdStokeswithpressure}, suppose moreover that $\Omega$ has a boundary of class ${\mathcal C}^2$, $\mathsf{B}\in W^{1,\infty}(\Omega)^{3\times 3}$ with $\det \mathsf{B} = 1$, $\mu_1$, $\mu_2$, and $\mu_3\in W^{1,\infty}(\Omega)$, and $\mathbf{b}f \in L^2(\Omega)^3$. Then, the solution $(\mathbf{v}, p)$ of the weak formulation \eqref{weakviscoelasticStokeswithpressure} of the viscoelastic Stokes problem is such that $\mathbf{v} \in V(\Omega)\cap H^2(\Omega)$ and $p\in H^1(\Omega)$, with 
  \mathbf{e}gin{equation}
    \label{h2regularityforh2vwithB}
    \|D^2 \mathbf{v}\|_{L^2(\Omega)} \leq \frac{c(\Omega)}{\alpha}\left( \|\mathbf{b}f\|_{L^2(\Omega)} + \frac{1}{\alpha}\|{\mathcal A}(\mathsf{B})\|_{W^{1,\infty}_{\lambda_1}(\Omega)}\|\mathbf{b}f\|_{H^{-1}(\Omega)}\right),
  \end{equation}
  and
  \mathbf{e}gin{equation}
    \label{h2regularityforh1pwithB}
    \|\mathbf{n}abla p\|_{L^2(\Omega)} \leq \frac{c(\Omega)}{\alpha}\|{\mathcal A}(\mathsf{B})\|_{L^\infty(\Omega)}\left( \|\mathbf{b}f\|_{L^2(\Omega)} + \frac{1}{\alpha}\|{\mathcal A}(\mathsf{B})\|_{W^{1,\infty}_{\lambda_1}(\Omega)}\|\mathbf{b}f\|_{H^{-1}(\Omega)}\right),
  \end{equation}
  where $\lambda_1$ is the Poincar\'e constant, and $c(\Omega)$ is a shape constant.  Furthermore, $(\mathbf{v},p)$ solves the viscoelastic Stokes system \eqref{viscoelasticStokesorg} almost everywhere on $\Omega$, with the first two equations in \eqref{viscoelasticStokesorg} holding in $L^2(\Omega)^3$ and $H^1(\Omega)^3$, respectively, and with the last equation holding for the trace of $\mathbf{v}$ on $\partial\Omega$ in the space $H^{3/2}(\partial \Omega)$.
\end{thm}

\mathbf{e}gin{thm}[High-order regularity]
  \label{hkplustworegularityforviscoelasticStokes}
  Under the hypotheses of Theorem \ref{weaksolutionweakviscoelasticdStokeswithpressure}, suppose moreover that $\Omega$ has a boundary of class ${\mathcal C}^{k+2}$, $\mathsf{B}\in W^{1,\infty}(\Omega)^{3\times 3}\cap W^{2,3}(\Omega)^{3\times 3} \cap H^{k+1}(\Omega)^{3\times 3}$ with $\det \mathsf{B} = 1$, $\mu_1$, $\mu_2$, and $\mu_3\in W^{1,\infty}(\Omega)\cap W^{2,3}(\Omega)\cap H^{k+1}(\Omega)$, and $\mathbf{b}f \in H^k(\Omega)^3$. Then, the solution $(\mathbf{v}, p)$ of the weak formulation \eqref{weakviscoelasticStokeswithpressure} of the viscoelastic Stokes problem is such that $\mathbf{v} \in V(\Omega)\cap H^{k+2}(\Omega)$ and $p\in H^{k+1}(\Omega)$. Furthermore, $(\mathbf{v},p)$ solves the viscoelastic Stokes problem \eqref{viscoelasticStokesorg} almost everywhere on $\Omega$, with the first two equations in \eqref{viscoelasticStokesorg} holding in $H^{k-1}(\Omega)^3$ and $H^k(\Omega)$, respectively, and with the last equation holding for the trace of $\mathbf{v}$ on $\partial\Omega$, in the space $H^{k+1/2}(\partial \Omega)$. Finally, depending on the value of $k$, the estimates \eqref{h2egularityforp}, \eqref{h2egularityforp}, \eqref{hkplus2regularityforv}, and \eqref{hkplus1regularityforp} hold with $\mathsf{A}$ replaced by ${\mathcal A}(\mathsf{B})$.
\end{thm}

\mathbf{e}gin{rmk}
  The estimates needed for ${\mathcal A}(\mathsf{B})$ in Theorems \ref{weaksolutionweakviscoelasticdStokeswithpressure}, \ref{h2regularityforviscoelasticStokes}, and \ref{hkplustworegularityforviscoelasticStokes} can be taken to be those obtained in Propositions \ref{proplinftyboundonacalb}, \ref{propwoneinftyboundonacalb}, \ref{propw23boundonacalb}, and \ref{prophmboundonacalb}, but, of course, other estimates may be worked out if desired.
\end{rmk}

\mathbf{e}gin{rmk}
  In a more general form, one obtains results analogous to Theorems \ref{weaksolutionweakviscoelasticdStokeswithpressure}, \ref{h2regularityforviscoelasticStokes}, and \ref{hkplustworegularityforviscoelasticStokes} without assuming explicit regularity conditions on $\mathsf{B}$ and instead assuming directly that the tensor $\mathsf{B}$ is such that ${\mathcal A}(\mathsf{B})$ satisfies the assumptions for $\mathsf{A}$ in Theorems \ref{weaksolutionweakgeneralizedStokeswithpressure}, \ref{h2regularity}, and \ref{hkplustworegularity}, respectively. This is particularly useful when the material can be modelled with $\mu_3=0$ or in special situations or with special symmetries in which the regularity of ${\mathcal A}(\mathsf{B})$ might be improved, or even if ${\mathcal A}(\mathsf{B})$ depends on $\mathsf{B}$ on a different way.
\end{rmk}

\mathbf{e}gin{rmk}
The same kind of results (in Theorems \ref{weaksolutionweakgeneralizedStokeswithpressure}, \ref{h2regularity} and \ref{hkplustworegularity}, for the generalized Stokes problem, and in Theorems \ref{weaksolutionweakviscoelasticdStokeswithpressure}, \ref{h2regularityforviscoelasticStokes}, and \ref{hkplustworegularityforviscoelasticStokes}, for the viscoelastic Stokes problem) hold under different homogenous boundary condition, such as fully periodic conditions or combinations of no-slip and periodic boundary conditions such as in a periodic channel. Non-homogenous boundary conditions can also be treated, provided they are the trace of functions with sufficient regularity. The details are omitted. 
\end{rmk}

\section{Conclusions}

We considered in detail a particular class of viscoelastic materials, modelled by system \eqref{viscoelasticStokesorg}. Our main interest was in modelling the deformation of different types of salts, needed for simulations of the exploitation of pre-salt oil fields. 

We obtained existence, uniqueness and regularity results for system \eqref{viscoelasticStokesorg}, based on natural conditions on the parameters $\mathsf{B}$ and $\mu_1, \mu_2,\mu_3$, and on the boundary $\Omega$. In particular, we identified the condition \eqref{essinfglambda} as the main condition for the uniform ellipticity of the problem, as well as the regularity conditions given in Section \ref{viscoelasticsec}.

For the existence, uniqueness and low-order regularity (up to $H^2$), we rewrote the system into a form suitable for the application of previously known results and verified all the necessary conditions.

For higher-order regularity (in $H^k$, $k\geq 3$), we exploited the Sobolev estimates valid in three dimensions and improved the regularity needed on the coefficients of the differential operator as compared with previously available results for similar elliptic systems.

We also obtained sharper and more explicit estimates for the regularity of the solution in terms of the parameters of the problem. These estimates are important for the well-posedness of the corresponding evolutionary problem to be presented in a forthcoming paper.

In regards to the ellipticity condition \eqref{essinfglambda}, which is of practical importance for checking the validity of a numerical simulation, we identified, in the case of constant coefficients $\mu_1, \mu_2, \mu_3$,  different intervals of the positive real line in which the eigenvalues of the left Cauchy-Green deformation tensor are allowed to belong to (see Lemma \ref{lemmuscenarios}). If the deformation is such that at least one of the eigenvalues escapes those intervals, the system is no longer elliptic. This explains the breakdown of our own simulations seen in some situations.

\mathbf{e}gin{thebibliography}{25}

\bibitem{agmon1965} S. Agmon, \emph{Lectures on Elliptic Boundary Value Problems}, Van Nostrand Mathematical Studies, Van Nostrand, 1965.

\bibitem{AB} C. D. Aliprantis and K. C. Border, \emph{Infinite dimensional analysis, A hitchhiker's guide}. Third edition. Springer-Verlag, 2006.

\bibitem{adamsfournier2003} R. A. Adams and J. F. Fournier, Sobolev Spaces, Academic Press, 2003.

\bibitem{adn2} S. Agmon, A. Douglis, and L. Nirenberg, Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions II, \emph{Comm. Pure Appl. Math.} 17 (1964), 35--92.

\bibitem{cheminmasmoudi2001} J.-Y. Chemin and N. Masmoudi, About lifespan of regular solutions of equations related to viscoelastic fluids, \emph{SIAM J. Math. Anal.} 33 (2001), no. 1, 84--112

\bibitem{ciarlet1988} P. G. Ciarlet, \emph{Mathematical Elasticity: Three-dimensional Elasticity}, Vol 1 (Studies in Mathematics \& its Applications),  Elsevier Science Ltd, 1988.

\bibitem{evans2010} L. C. Evans, \emph{Partial Differential Equations}, Graduate Studies in Mathematics, vol. 19, Second Edition, American Mathematical Society, Providence, Rhode Island, 2010.

\bibitem{friedman1969} A. Friedman, \emph{Partial Differential Equations}, Holt, Rinehart and Winston, New York-Montreal-London, 1969.

\bibitem{ghidaglia1984} J.-M. Ghidaglia, Régularité des solutions de certains probl\`emes aux limites lin\'eaires li\'es aux \'equations d'Euler, \emph{Comm. Partial Differential Equations} 9 (1984), no. 13, 1265--1298. 

\bibitem{gilbargtrudinger2001} D. Gilbarg and N. S. Trudinger, \emph{Elliptic Partial Differential Equations of Second Order,} Springer-Verlag, Berlin, Heidelberg, 2001.

\bibitem{GreeneKnuth1982} D. H. Greene and D. E. Knuth, \emph{Mathematics for the analysis of algorithms.} Reprint of the third (1990) edition. Modern Birkh\"auser Classics. Birkh\"auser Boston, Inc., Boston, MA, 2008.

\bibitem{giaquintamodica1982} M. Giaquinta and G. Modica, Non linear systems of the type of the stationary Navier-Stokes system, \emph{J. Reine Angew. Math.} 330 (1982), 173--214.

\bibitem{guillopesaut1990} C. Guillop\'e and J.-C. Saut, Existence results for the flow of viscoelastic fluids with a differential constitutive law, \emph{Nonlinear Anal.} 15 (1990), no. 9, 849--869.

\bibitem{haupt2002} P. Haupt, \emph{Continuum Mechanics and Theory of Materials,} Second Edition, Springer-Verlag, 2002.

\bibitem{huystara2006} N. D. Huy and J. Star\'a, On existence and regularity of solutions to a class of generalized stationary Stokes problem, \emph{Comment. Math. Univ. Carolin.} \textbf{47}, (2006), no. 2, 241--264.

\bibitem{laxmilgram} P. D. Lax, A. N. Milgram, Parabolic equations,  Contributions to the theory of partial differential equations. \emph{Annals of Mathematics Studies,} no. 33. (1954), Princeton, N. J.: Princeton University Press. 167--190.

\bibitem{lionsmagenes1972} J. L. Lions and E. Magenes, \emph{Non-Homogeneous Boundary Value Problems and Applications}, Volume 1, Springer-Verlag, New York, Berlin, Heidelberg, 1972.

\bibitem{lionsmasmoudi2000} P. L. Lions and N. Masmoudi, Global solutions for some Oldroyd models of non-Newtonian flows, \emph{Chinese Ann. Math. Ser. B} 21 (2000), no. 2, 131--146.

\bibitem{linliuzhang2005} F.-H. Lin, C. Liu, and P. Zhang, On hydrodynamics of viscoelastic fluids, \emph{Comm. Pure Appl. Math.} 58 (2005), no. 11, 1437--1471.

\bibitem{linzhang2008} F.-H. Lin and P. Zhang,  On the initial-boundary value problem of the incompressible viscoelastic fluid system,  \emph{Comm. Pure Appl. Math.} 61 (2008), no. 4, 

\bibitem{LiuCipolattiRinconPalermo2014} I-S. Liu, R. A. Cipolatti, M. A. Rincon, L. A. Palermo, Successive linear approximation for large deformation -- instability of salt migration, \emph{Journal of Elasticity,} 114 (2014), no. 1, 19--39.

\bibitem{liuwalkington2001} C. Liu and N. J. Walkington, An Eulerian description of fluids containing visco-hyperelastic particles. \emph{Arch. Ration. Mech. Anal.} 159 (2001), no. 3, 229--252.

\bibitem{Liu2008} I-S. Liu, Entropy Flux Relation for Viscoelastic Bodies, \emph{Journal of Elasticity,} 90 (2008), 259--270.

\bibitem{Liu2002} I.-S. Liu, \emph{Continuum Mechanics.} Springer-Verlag, Berlin, 2002.

\bibitem{macha2011} V. M\'acha, On a generalized Stokes problem, \emph{Cent. Eur. J. Math.} 9 (2011), no. 4, 874--887.

\bibitem{necas} J. Ne\v{c}as, \emph{Equations aux derivees partielles,} Presses de Universit\'e de Montr\'eal, Montreal, 1965. 

\bibitem{sideristhomases2005} T. C. Sideris and B. Thomases, Global existence for three-dimensional incompressible isotropic elastodynamics via the incompressible limit. \emph{Comm. Pure Appl. Math.} 58 (2005), no. 6, 750--788.

\bibitem{solonnikov2001} V. A. Solonnikov, Initial-boundary value problem for generalized Stokes equations, Math. Bohem. 126 (2001), no. 2, 505--519.

\bibitem{teixeiraliurinconcipolattipalermo2014} M. G. Teixeira, I-S. Liu, M. A. Rincon, R. A. Cipolatti, and L. A. C. Palermo, The influence of temperature on the formation of salt domes, International Journal of Modeling and Simulation for the Petroleum Industry (Online) 8 (2014), no. 2, 35--41.

\bibitem{temam} R. Temam \emph{Navier-Stokes Equations. Theory and Numerical Analysis,} Studies in Mathematics and its Applications 2,  3rd edition, North-Holland Publishing Co., Amsterdam-New York, 1984. Reedition in 2001 in the AMS Chelsea series, AMS, Providence.
  
\bibitem{temammiranville2005} R. Temam and A. Miranville, \emph{Mathematical modeling in continuum mechanics.} Second edition. Cambridge University Press, Cambridge, 2005. 
  
\bibitem{truesdellnoll2004} C. Truesdell and W. Noll, \emph{The non-linear field theories of mechanics.}
Third edition. Edited and with a preface by Stuart S. Antman. Springer-Verlag, Berlin, 2004.  
  
\end{thebibliography}

\end{document}