\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{amssymb}
\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2019 (2019), No. 109, pp. 1--29.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2019 Texas State University.}
\vspace{8mm}}
\begin{document}
\title[\hfilneg EJDE-2019/109\hfil Variable Lorentz estimate]
{Variable Lorentz estimate for generalized Stokes systems in non-smooth domains}
\author[S. Liang, S. Zheng, Z. Feng \hfil EJDE-2019/109\hfilneg]
{Shuang Liang, Shenzhou Zheng, Zhaosheng Feng}
\address{Shuang Liang \newline
Department of Mathematics,
Beijing Jiaotong University,
Beijing 100044, China}
\email{shuangliang@bjtu.edu.cn}
\address{Shenzhou Zheng (corresponding author) \newline
Department of Mathematics,
Beijing Jiaotong University,
Beijing 100044, China}
\email{shzhzheng@bjtu.edu.cn}
\address{Zhaosheng Feng \newline
Department of Mathematics,
University of Texas Rio Grande Valley,
Edinburg, TX 78539, USA}
\email{zhaosheng.feng@utrgv.edu}
\thanks{Submitted March 2, 2018. Published September 26, 2019.}
\subjclass[2010]{35D30, 35J47, 76D07}
\keywords{Generalized Stokes systems; Lorentz estimates with variable power;
\hfill\break\indent small BMO; Reifenberg flatness; large-M-inequality principle}
\begin{abstract}
We prove a global Calder\'on-Zygmund type estimate in the framework of
Lorentz spaces for the variable power of the gradient of weak solution
pair $(u,P)$ to the generalized steady Stokes system over a bounded
non-smooth domain. It is assumed that the leading coefficients satisfy the
small BMO condition, the boundary of domain belongs to Reifenberg flatness,
and the variable exponent $p(x)$ is log-H\"older continuous.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks
\section{Introduction}
Let $\Omega\subset \mathbb{R}^n (n\ge 2)$ be a given bounded domain with a
rough boundary specified later.
The aim of this article is to study a global Lorentz estimate for the variable
power of the gradient of weak solution to the generalized steady Stokes problem
\begin{equation}\label{Stokes-systems}
\begin{gathered}
\operatorname{div}(\mathbf{A}(x)\nabla u)-\nabla P
= \operatorname{div} \mathbf{F}, \quad \text{in }\Omega,\\
\operatorname{div} u=0, \quad \text{in }\Omega,\\
u=0, \quad \text{on }\partial\Omega.
\end{gathered}
\end{equation}
Throughout this article, as usual we assume that the fourth-order tensor
$\mathbf{A}(x)=\big(A^{\alpha\beta}_{ij}\big)^n_{i,j,\alpha,\beta=1}:
\Omega\to \mathbb{R}^{n^2\times n^2}$ satisfies a uniform boundedness
and ellipticity for constants $0 < \nu \leq \Lambda < +\infty $:
\begin{equation} \label{uniform-bound-ellipticity}
\nu|\xi|^2 \leq \mathbf{A}(x)\xi\cdot\xi \leq \Lambda |\xi|^2,
\end{equation}
where $x\in \Omega$ a.e., $\xi \in \mathbb{R}^{n^2} $, and
$\mathbf{F}(x)=(F^{\alpha}_{i})^n_{i,\alpha=1}$.
The unknown velocity of vectorial-value functions is denoted by
$u=(u^1,u^2,\dots,u^n): \Omega\to \mathbb{R}^n$
and $P:\Omega\to \mathbb{R}$ is the pressure.
Let us recall some recent progresses of the Calder\'on-Zygmund theory concerning
partial differential equations with discontinuous coefficients.
The interior and global $W^{2,p}$ estimates for nondivergence linear elliptic
equations with the VMO discontinuous coefficients
were presented by Chiarenza-Frasca-Longo \cite{ChFL91,ChFL93}.
Since then, there has been continuous attention on the Calder\'{o}n-Zygmund
theory for various elliptic and parabolic problems with discontinuous coefficients.
Apart from an earlier technique by using singular integral
operators and its commutators, there are three kinds of important arguments to
deal with the Calder\'{o}n-Zygmund theory concerning elliptic and parabolic
problems with the VMO or small BMO discontinuous coefficients.
The first one is the so-called geometrical approach originally traced from
Byun-Wang's work \cite{ByW04}, which is used to attain the global $L^p$
estimate by way of weak compactness based on the boundedness of the
Hardy-Littlewood maximal operators and the modified Vitali covering theory
of distributional functions regarding the gradients of solutions.
Indeed, this can also be regarded as a development from Caffarelli-Peral's
work \cite{CaP98} to obtain local $W^{1,p}_{\rm loc}$-estimates for solutions
of a class of elliptic problems of $p$-Laplace type.
Secondly, Dong-Kim-Krylov (for examples, see \cite{DoK11,Kr07}) presented a
unified approach of studying $L^p$ solvability for elliptic and parabolic
problems on the basis of the Fefferman-Stein theorem on sharp functions and the
Hardy-Littlewood maximal function theorem for the spatial derivatives of solutions.
The third technique is called the large-M-inequality principle originated
from Acerbi-Mingione's work \cite{AcM05,AcM07}, which is directly applied to
argue on certain Calder\'on-Zygmund-type coverings instead of the maximal
function operator and other harmonic
techniques such as the good-$\lambda$-inequality. Here, we would like
to mention that recently Byun et al have obtained
numerous global Calder\'{o}n-Zygmund-type results to various nonlinear
elliptic and parabolic problems over non-smooth domains by combining the
large-M-inequality principle with a geometrical approach \cite{ByOW14, ByW04, ByS17}.
As we have seen, Byun-Ok-Wang \cite{ByOW14} attained a global Calder\'on-Zygmund
estimate with the variable exponent of gradients
of solution to the zero Dirichlet problem for linear elliptic systems in the
divergence form with partial BMO coefficients
and log-H\"older continuity $p(x)$, which implies that
\[
\mathbf{F}\in L^{p(x)}(\Omega,\mathbb{R}^n) \Rightarrow
Du\in L^{p(x)}(\Omega,\mathbb{R}^n).
\]
Later, Tian-Zheng \cite{TiZh17} further extended it to a global
Calder\'on-Zygmund-type estimate for variable power of the gradient of
solution in Lorentz spaces for the same problem with partial BMO coefficients.
We know that solvability and optimal regularity of Stokes system under minimal
regular datum have been a classical and important problem
in the theory of partial differential equations and fluid dynamics.
In the past decades, we have seen a great deal of literature concerning
the interior regularity of Stokes system \cite{Bre13,DiK13,Die14}
and global regularity of the generalized Stokes system in the Lipschitz domain
(cf. \cite{BuBS16, DaJSK11, DiLR11, DoK16, GuS15, Mac11}), the domain with
the Reifenberg flatness boundary \cite{ByS17,HuLW14}.
It is quite necessary to mention some recent advances concerning the generalized
steady Stokes problems \eqref{Stokes-systems}
with discontinuous coefficients.
Dan\v{e}\v{c}ek-John-Star\'{a} \cite{DaJSK11} investigated the Morrey regularity
for the gradient of weak solution $(u, P)$ to a generalized Stokes system with
symmetric elliptic coefficients whose entries
satisfy the boundedness and VMO discontinuity.
Gu-Shen \cite{GuS15} considered a uniform $W^{1,p}$-regularity in the
homogenization theory of the generalized Stokes system with setting of rapidly
oscillating periodic coefficients, and obtained the global $W^{1,p}$-estimates
with $1
p$ and $q\in(0,\infty]$. Later, Zhang-Zhou \cite{ZhZ14} extended
the result of \cite{MeP12} to the quasilinear elliptic $p(x)$-Laplacian
equations by using a geometrical argument. Adimurthi-Phuc \cite{AdP15} showed
global Lorentz and Lorentz-Morrey estimates below the natural exponent
of quasilinear equations. Zhang-Zheng \cite{ZhZ16, ZhZ17} studied Lorentz
estimates of fully nonlinear parabolic and
elliptic equations with small BMO nonlinearities, and obtained weighted Lorentz
estimates of the Hessian of strong solution for nondivergence linear elliptic
equations with partial BMO coefficients.
Motivated by the progresses mentioned above, in this work we focus on a global
Calder\'on-Zygmund type estimate for the variable power
of the gradient of weak solution in the framework of Lorentz spaces to
the Dirichlet problem of the generalized steady Stokes
systems \eqref{Stokes-systems}
in the non-smooth domain. Here, we allow the coefficient tensor
$ \mathbf{A}(x)$ to be discontinuous, but it suffices to impose
a small BMO regular condition,
the variable exponent $p(x)$ satisfies log-H\"older continuity, and the
boundary of domain belongs to Reifenberg flatness.
This study is also inspired by elegant results presented in
\cite{AcM05,AcM07,Bar13,Bar14, ByS17}. That is, we apply the mixed argument
of large-M-inequality principle and the geometric approach to prove global
Lorentz estimates for the variable power of the gradient of weak solution
to the Dirichlet problem of \eqref{Stokes-systems} over a bounded Reifenberg
flatness domain.
Indeed, the key ingredient is based on making use of Calder\'on-Zygmund type
covering, approximate estimate and iteration arguments to obtain an estimate
of the measure of the super-level
set for the variable power of the gradient of weak solution.
The rest of this artivcle is organized as follows.
In section 2, we recall the definition of weak solution to problem
\eqref{Stokes-systems} and state our main result.
In Section 3 we present some technical lemmas.
In Section 4 we prove of our main result.
\section{Preliminaries and statement of main results}
In this section, we present some related definition and notations and state
our main result on the Dirichlet problem of
the generalized steady Stokes systems \eqref{Stokes-systems}.
Denote by $c(n,\nu,\Lambda,\dots)$ a universal constant
depending only on prescribed quantities and possibly varying from line to
line in the following context.
Let us recall the Lorentz space $L(t,q)(U)$ with an open subset
$U\subset \mathbb{R}^n$ for any
parameters $1\leq t<\infty$ and $0\mu\}|\Big)^{q/t}\frac{d\mu}{\mu}<\infty;
\]
while the Lorentz space $L(t,\infty)$ for $1\leq t<\infty$ and $q=\infty$ is
defined by the Marcinkiewicz space $\mathcal{M}^{t}(U)$ as usual, which is
the space of measurable functions $g$ with
\[
\|g\|_{L(t,\infty)}
=\|g\|_{\mathcal{M}^{t}(U)}:=\sup_{\mu>0}
\Big(\mu^{t}|\{\xi \in U: |g(\xi)|>\mu\}|\Big)^{1/t}<\infty.
\]
The local variance of such spaces is defined in the usual way. We remark that
if $t=q$, then the Lorentz space $L{(t,t)}(U)$ is
nothing but a classical Lebesgue space. Indeed, by Fubini's theorem it gives
\[
\|g\|^{t}_{L^{t}(U)}=t\int^{\infty}_{0}\mu^{t}
|\{\xi \in U: |g(\xi)|>\mu \} |\frac{d\mu}{\mu}
=\|g\|^{t}_{L(t,t)(U)},
\]
which implies $L^{t}(U)=L(t,t)(U)$, see \cite[4,5,26,31,\cite{ ZhZ16}]{ZhZ16}.
\begin{remark}\label{Lorentz-q-nor} \rm
Because of the lack of sub-additivity, the quantity $\|\cdot\|_{L(t,q)(U)}$
is just a quasi-norm. Nevertheless, the mapping $g\mapsto\|g\|_{L(t,q)(U)}$
is still weak lower semi-continuous, for details see \cite[Remark 3]{Min10}.
\end{remark}
\begin{definition}\label{def-weak-solution} \rm
Let $\mathbf{F} \in L^2(\Omega)^{n^2}$.
If $u \in W_0^{1,2}(\Omega)^n$, $\operatorname{div} u=0$ and satisfies
\begin{equation}\label{weak-ellip-equa-a}
\int_\Omega \langle \mathbf{A}(x)\nabla u, \nabla v\rangle\,dx
= \int_\Omega \langle \mathbf{F}, \nabla v\rangle\,dx
\end{equation}
for all $v \in W_0^{1,2}(\Omega)^n$ and $\operatorname{div} v=0$,
then $u$ is called a weak solution of the zero
Dirichlet problem of the generalized Stokes system \eqref{Stokes-systems}.
If $u$ is such a weak solution and $ P \in L^2(\Omega)$ satisfies
\begin{equation}\label{weak-ellip-equa-b}
\int_\Omega \big(\langle \mathbf{A}(x)\nabla u, \nabla v\rangle
-\langle P, \operatorname{div} v \rangle \big)dx
= \int_\Omega \langle \mathbf{F}, \nabla v\rangle\,dx
\end{equation}
for all $v \in W_0^{1,2}(\Omega)^n$, then $(u,P)$ is said to be a weak
solution pair of the generalized Stokes system \eqref{Stokes-systems}
and $P$ is called an associated pressure of $u$.
\end{definition}
The traditional assumption on the variable exponent $p(\cdot)$ is
log-H\"older continuity, which ensures that
the Hardy-Littlewood maximal operator is bounded within the framework of
generalized Lebesgue space. Briefly,
we recall that $p(x)$ is log-H\"older continuous, denoted by $p(x)\in LH(\Omega)$,
if there exist constants $c_{0}$ and $\delta>0$ such that for all
$x,y \in \Omega$ with $|x-y|< \delta$, it holds
\[
|p(x)-p(y)| \le \frac{c_{0}}{-\log(|x-y|)}.
\]
In this context, we assume that $p(x): \Omega \to \mathbb{R}$ is a
log-H\"older continuous function, and there exist positive constants
$\gamma_{1}$ and $\gamma_{2}$ such that
\begin{gather}\label{bound-exponent}
2 <\gamma_{1} \le p(x) \le \gamma_{2} < \infty, \quad \forall x\in \Omega, \\
\label{exponent-modulus-continu}
|p(x)-p(y)| \le \omega(|x-y|), \quad \forall x, y \in \Omega,
\end{gather}
where $\omega: [0,\infty)\to [0,\infty)$ is a modulus of continuity of $p(x)$.
Without loss of generality, we suppose that $\omega$ is a non-decreasing
continuous function with
$$
\omega(0)=0,\quad \operatorname{lim\,sup}_{r \to 0} \omega(r)
\log\big( \frac{1}{r} \big) <\infty.
$$
With the above assumptions, it is clear that $p(x)\in LH(\Omega)$ and there
exists a positive constant $A$ such that
\begin{equation}\label{log-equal}
\omega(r) \log(\frac{1}{r}) \le A \;\Leftrightarrow\;
r^{-\omega(r)} \le e^{A} \text{ for any } r\in (0,1).
\end{equation}
On the other hand, the generalization of the classical steady Stokes
system consists of general second order elliptic equations
in the divergence form instead of the standard Stokes equation.
It is rather necessary to impose certain regular assumptions on the leading
coefficients $\mathbf{A}(x)$ and the geometric structure on the boundary
of domain. To this end, we let
\[
B_r(y)=\{x\in \mathbb{R}^n:|x-y|0$. Denote that
\begin{gather*}
\Omega_r(y)=B_r(y)\cap\Omega, \quad
\partial_\omega\Omega_r(y)=B_r(y)\cap\partial\Omega,\\
B_r^+=B_r(0)\cap\{x_n>0\},\quad
T_r=B_r(0)\cap\{x_n=0\}.
\end{gather*}
For any bounded domain $U\subset \mathbb{R}^n$, we denote
\[
f_U= {-\hspace{-0.38cm}\int}_Uf(x)dx = \frac{1}{|U|}\int_Uf(x)dx.
\]
In what follows, a key assumption is that the coefficient tensor
$\mathbf{A}(x)$ is allowed to be sufficiently small BMO discontinuous,
and the boundary $\partial\Omega$ of domain is Reifenberg flat.
\begin{definition}\label{small BMO} \rm
We say that the pair $(\mathbf{A}, \Omega)$ is $(\delta,R_0)$-vanishing,
if for any $x \in \Omega$ and for each number $r\in (0,R_0]$ with
\[
\operatorname{dist}(x, \partial \Omega)
= \min_{z\in \partial \Omega} \operatorname{dist}(x, z) > \sqrt{2} r,
\]
there exists a coordinate system depending on $x$ and $r$ such that in the
new coordinate system $x$ is the origin and satisfies
\begin{equation}\label{small-BMO-condition-a}
{-\hspace{-0.38cm}\int}_{B_r(x)}\big|\mathbf{A}(y)-\mathbf{A}_{B_r(x)}\big|^2dy
\leq \delta^2.
\end{equation}
While, for any $x \in \Omega$ and for each number $r\in (0,R_0]$ with
\[
\operatorname{dist}(x, \partial \Omega)
= \min_{z \in \partial \Omega} \operatorname{dist}(x, z)
=\operatorname{dist}(x, z_0)\le \sqrt{2} r
\]
for some $z_0\in \partial\Omega$, there exists a coordinate system depending
on $x$ and $r$ such that
in the new coordinate system $z_{0}$ is the origin and satisfies
\begin{gather}\label{Ref-bou}
B_{3r}(z_0) \cap \{x_{1} \ge 3 \delta r \}
\subset B_{3r}(z_0)\cap \Omega \subset B_{3r}(z_0)
\cap \{x_{1}\ge - 3\delta r \}, \\
{-\hspace{-0.38cm}\int}_{B_{3r}(z_{0})}
|\mathbf{A}(y)-\mathbf{A}_{B_{3r}(z_{0})}|^2 \,dx
\le \delta^2,
\end{gather}
where $A(x)$ is a zero-extension from $B_{3r}(z_{0})\cap \Omega$ to $B_{3r}$,
and the parameters $\delta>0$ and $R_0$ will be specified later.
\end{definition}
We assume that $\delta$ is a small positive constant, saying $0<\delta<1/8$.
Notice that $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain.
Obviously, it is an $A$-type domain,
which implies the following measure density condition \cite{ByW04}:
\begin{equation}\label{thick-a}
\sup_{00$.
Suppose that $(u,P)$ satisfying $u \in W_0^{1,2}(\Omega)^n$,
$\operatorname{div} u=0$ and is a weak solution pair of
the generalized Stokes system
\eqref{Stokes-systems}--\eqref{uniform-bound-ellipticity}. If
\[
|F|^{p(x)} \in L(t,q)(\Omega),\quad \text{for $t\geq 1$ and $q \in (0,+\infty]$},
\]
then there exists a small constant
$\delta_0=\delta_0(n,\gamma_1,\gamma_2,\nu,\Lambda)>0$ such that for every
$\delta \in (0,\delta_0]$, we have $ (|Du|+|P|)^{p(x)} \in L(t,q)(\Omega)$
with the estimate
\begin{equation}\label{Lorentz-estimate-a}
\big\| (|Du|+|P|)^{p(x)}\big\|_{L(t,q)(\Omega)}
\leq c\Big(\||F|^{p(x)}\|_{L(t,q)(\Omega)}+1\Big)^{\gamma_2/\gamma_1},
\end{equation}
where the constant $c$ depends only on
$n, \gamma_1, \gamma_2, \nu, \Lambda, t, q, \delta_0, R_0, \Omega,
\omega(\cdot)$ and $|\Omega|$ (except in the case $q=\infty$).
\end{theorem}
\section{Technical tools}
In this section, we present some useful technical lemmas, which will play
an essential role in proving our main result. We start
with recalling the existence and energy estimate of weak solution pair to
the generalized Stokes system \eqref{Stokes-systems},
see \cite[Lemma 2.9]{ByS17}.
\begin{lemma}\label{existence-energy-estimate}
Let $\Omega \subset \mathbb{R}^n $ for $ n\geq 2$ be an open bounded
$(\delta,R_0)$-Reifenberg flat domain with
sufficiently small $\delta>0$, and $\mathbf{F} \in L^2(\Omega)^{n^2}$.
Then there exists a unique solution
pair $(u,P)\in W_0^{1,2}(\Omega)^n\times L^2(\Omega)$ to the generalized
Stokes system \eqref{Stokes-systems}
with $\operatorname{div} u=0$ and $\int_\Omega P\,dx =0 $ such that the
following standard estimate holds
\begin{equation}\label{existence-energy-estimate-a}
\|\nabla u\|_{L^2(\Omega)^{n^2}}+\|P\|_{L^2(\Omega)}
\leq c \|\mathbf{F}\|_{L^2(\Omega)^{n^2}},
\end{equation}
where $c=c(n,\nu,\Lambda,\Omega)$.
In addition, if $u\in W_0^{1,q}(\Omega)^n$ with $\operatorname{div}u=0$
and $\mathbf{F} \in L^q(\Omega)^{n^2}$ for $q \in [2,+\infty)$, then
\begin{equation}\label{existence-energy-estimate-b}
\|P\|_{L^q(\Omega)}
\leq c\big( \|\mathbf{F}\|_{L^q(\Omega)^{n^2}}+\|\nabla u\|_{L^q(\Omega)^{n^2}}\big),
\end{equation}
where $c=c(n,\nu,\Lambda,\Omega)$.
\end{lemma}
The following lemma is regarding the principle of local reverse H\"older's
inequality, which can be obtained form \cite[Prop. 1.2 Chapter 5]{Gia83}.
\begin{lemma}\label{reverse-Holder-b-1}
Let $0p$ such that there is $g\in L_{\rm loc}^{p'}(\Omega,\mathbb{R}^N)$
with the estimate
\[
\Big({-\hspace{-0.38cm}\int}_{B_R}|g(x)|^{p'}dx\Big)^{\frac{1}{p'}}
\leq c \Big({-\hspace{-0.38cm}\int}_{B_{2R}}|g(x)|^{p}dx\Big)^{1/p}
+c\Big({-\hspace{-0.38cm}\int}_{B_{2R}}|h(x)|^{p}dx\Big)^{1/p}.
\]
\end{lemma}
Based on Lemmas \ref{existence-energy-estimate} and \ref{reverse-Holder-b-1},
we now prove a higher integrability for the gradients of weak solutions
and the pressure $P$ to the
generalized Stokes system \eqref{Stokes-systems} in the admissible
set $ W_0^{1,2}(\Omega)^n\times L^2(\Omega)$ with $\operatorname{div} u=0$.
\begin{lemma}\label{reverse-Holder-a}
Let $(u,P)\in W_0^{1,2}(B_{2r})^n\times L^2(B_{2r})$ be a weak solution pair
to the generalized Stokes system \eqref{Stokes-systems}
with $\operatorname{div} u=0$ under the usual assumption
\eqref{uniform-bound-ellipticity}. If $|\mathbf{F}|^{p(x)} \in L^t(B_{2r})$
with $p(x)>\gamma_1>2 $ and $t> 1$ with $B_{2r} \Subset \Omega$,
then there exist positive constants $c=c(n,\gamma_1,\gamma_2,\nu,\Lambda)$ and
small $\sigma_0>0$ with
\[
\sigma_0< \frac{t\gamma_1}{2}-1,
\]
such that for any $0<\sigma \leq\sigma_0$, it holds
\begin{equation}\label{reverse-Holder-ineq-a}
\begin{aligned}
&\Big({-\hspace{-0.38cm}\int}_{B_r}|\nabla u|^{2(1+\sigma)}dx
\Big)^{\frac{1}{1+\sigma}}
+\Big({-\hspace{-0.38cm}\int}_{B_r}|P|^{2(1+\sigma)}dx\Big)^{\frac{1}{1+\sigma}} \\
&\leq c{-\hspace{-0.38cm}\int}_{B_{2r}}|\nabla u|^2dx
+c{-\hspace{-0.38cm}\int}_{B_{2r}}|P|^2dx
+c\Big({-\hspace{-0.38cm}\int}_{B_{2r}}|\mathbf{F}|^{2(1+\sigma)}dx
\Big)^{\frac{1}{1+\sigma}}.
\end{aligned}
\end{equation}
\end{lemma}
\begin{proof}
Let $\eta \in C_0^\infty(B_{2r})$ be a cut-off function such that
$0\leq \eta \leq 1$, $\eta\equiv 1$ on $B_r$ and
$|\nabla\eta| \leq 2/r$. Taking $v=\eta^2(u-(u)_{2r})$
into \eqref{weak-ellip-equa-a} as a test function, we obtain
\begin{equation}\label{weak-ellip-equa-a-1}
\int_{B_{2r}} \langle \mathbf{A}(x)\nabla u,\nabla\big(\eta^2(u-(u)_{2r})\big)
\rangle\,dx
=\int_{B_{2r}} \langle \mathbf{F},\nabla\big(\eta^2(u-(u)_{2r})\big)\rangle\,dx.
\end{equation}
In view of the ellipticity and boundedness \eqref{uniform-bound-ellipticity},
for $0<\varepsilon_1<1$, using Young's inequality we have
\begin{align*}
&\nu\int_{B_{2r}}\eta^2|\nabla u|^2dx \\
&\leq \int_{B_{2r}} \langle \mathbf{A}(x)\nabla u,\nabla\big(\eta^2(u-(u)_{2r})\big)
\rangle\,dx
+2\Lambda \int_{B_{2r}}|\eta \nabla u|\cdot|\nabla\eta(u-(u)_{2r})|dx \\
&\leq \int_{B_{2r}} \langle\mathbf{A}(x)\nabla u,\nabla\big(\eta^2(u-(u)_{2r})\big)
\rangle\,dx
+ c\varepsilon_1\int_{B_{2r}}\eta^2|\nabla u|^2dx \\
&\quad + c(\varepsilon_1)\int_{B_{2r}}|\nabla\eta |^2|u-(u)_{2r}|^2dx.
\end{align*}
It follows from the Sobolev-Poinc\'{a}re inequality with $2_*=\frac{2n}{n+2}$ that
\begin{equation}\label{left-hand-side}
\begin{aligned}
&\nu\int_{B_{2r}}\eta^2|\nabla u|^2dx\\
&\leq \int_{B_{2r}} \langle \mathbf{A}(x)\nabla u,\nabla\big(\eta^2(u-(u)_{2r})\big)
\rangle\,dx
+c\varepsilon_1\int_{B_{2r}}\eta^2|\nabla u|^2dx \\
&\quad + c(\varepsilon_1)r^{n-\frac{2n}{2_*}}
\Big(\int_{B_{2r}}|\nabla u|^{2_*}dx\Big)^{2/2_*}.
\end{aligned}
\end{equation}
Similarly, we can find the estimate to the right-hand side of formula
\eqref{weak-ellip-equa-a-1}:
\begin{equation}\label{right-hand-side}
\begin{aligned}
&\int_{B_{2r}} \langle \mathbf{F},\nabla\big(\eta^2(u-(u)_{2r})\big)\rangle\,dx \\
&\leq \varepsilon_2\int_{B_{2r}}\eta^2|\nabla u|^2dx \\
&\quad + c(\varepsilon_2,\varepsilon_3)\int_{B_{2r}}|\mathbf{F}|^2dx
+ c(\varepsilon_3)r^{n-\frac{2n}{2_*}}
\Big(\int_{B_{2r}}|\nabla u|^{2_*}dx\Big)^{2/2_*}
\end{aligned}
\end{equation}
for $0<\varepsilon_2$ and $\varepsilon_3<1$.
Using \eqref{weak-ellip-equa-a-1}-\eqref{right-hand-side} and choosing
$c\varepsilon_1+\varepsilon_2=\nu/2$ yields
\[
{-\hspace{-0.38cm}\int}_{B_{r}}|\nabla u|^2dx
\leq c\Big({-\hspace{-0.38cm}\int}_{B_{2r}}|\nabla u|^{2_*}dx\Big)^{2/2_*}
+c{-\hspace{-0.38cm}\int}_{B_{2r}}|\mathbf{F}|^2dx.
\]
From Lemma \ref{reverse-Holder-b-1} by taking
$g(x)=\nabla u,\ h(x)=\mathbf{F}$, $p=2$, $\theta=0$ and $ s=2_*<2$,
it follows that there exists a $\sigma$ satisfying $2<2(1+\sigma)< t\gamma_1$
such that
\begin{equation}\label{Du-a}
\begin{aligned}
\Big({-\hspace{-0.38cm}\int}_{B_{r}}|\nabla u|^{2(1+\sigma)}dx
\Big)^{\frac{1}{1+\sigma}}
&\leq c{-\hspace{-0.38cm}\int}_{B_{2r}}|\nabla u|^2dx
+c{-\hspace{-0.38cm}\int}_{B_{2r}}|\mathbf{F}|^2dx \\
&\leq c{-\hspace{-0.38cm}\int}_{B_{2r}}|\nabla u|^2dx
+c\Big({-\hspace{-0.38cm}\int}_{B_{2r}}|\mathbf{F}|^{2(1+\sigma)}dx
\Big)^{\frac{1}{1+\sigma}}.
\end{aligned}
\end{equation}
By \eqref{Du-a} and \eqref{existence-energy-estimate-b} with
$q=2((1+\sigma))> 2$, we obtain
\begin{equation}\label{p-a}
\begin{aligned}
&\Big({-\hspace{-0.38cm}\int}_{B_{r}}|P|^{2(1+\sigma)}dx\Big)^{\frac{1}{2(1+\sigma)}}\\
&\leq c\Big({-\hspace{-0.38cm}\int}_{B_{2r}}|\nabla u|^{2(1+\sigma)}dx
\Big)^{\frac{1}{2(1+\sigma)}}
+c\Big({-\hspace{-0.38cm}\int}_{B_{2r}}|\mathbf{F}|^{2(1+\sigma)}dx
\Big)^{\frac{1}{2(1+\sigma)}} \\
&\leq c\Big({-\hspace{-0.38cm}\int}_{B_{2r}}|\nabla u|^2dx\Big)^{\frac{1}{2}}
+c\Big({-\hspace{-0.38cm}\int}_{B_{2r}}|\mathbf{F}|^{2(1+\sigma)}dx
\Big)^{\frac{1}{2(1+\sigma)}}.
\end{aligned}
\end{equation}
Combining \eqref{Du-a} and \eqref{p-a}, we arrive at the desired estimate
\eqref{reverse-Holder-ineq-a}.
\end{proof}
The following higher integrability on the boundary version is a self-improving
result due to the Reifenberg flatness condition of domain being an $A$-type
domain as shown in \eqref{thick-a}.
\begin{lemma}\label{reverse-Holder-b}
Let $(u,P)\in W_0^{1,2}(\Omega_{2r})^n\times L^2(\Omega_{2r})$ be a weak
solution pair to the generalized Stokes system \eqref{Stokes-systems}
with $\operatorname{div} u=0$ under condition \eqref{uniform-bound-ellipticity}.
Suppose that $|\mathbf{F}|^{p(x)} \in L^t(\Omega_{2r})$ with $p(x)>\gamma_1>2 $
and $t> 1$, and the boundary of $\Omega$ satisfies local
$(\delta,R_0)$-Reifenberg flatness: there exist $R_0>0$
and $\delta_0>0$ such that for $0-4\delta r\}.
\]
Then there exists a small positive constant $\sigma_0>0$ with
$\sigma_0< \frac{t\gamma_1}{2}-1$
such that for any $0<\sigma \leq\sigma_0$, we have
$\nabla u\in L^{2(1+\sigma)}(\Omega_r)$ with the estimate
\begin{equation}\label{reverse-Holder-ineq-b}
\begin{aligned}
&\Big({-\hspace{-0.38cm}\int}_{\Omega_r}|\nabla u|^{2(1+\sigma)}dx
\Big)^{\frac{1}{1+\sigma}}
+\Big({-\hspace{-0.38cm}\int}_{\Omega_r}|P|^{2(1+\sigma)}dx\Big)^{\frac{1}{1+\sigma}}\\
&\leq c{-\hspace{-0.38cm}\int}_{\Omega_{2r}}|\nabla u|^2dx
+c{-\hspace{-0.38cm}\int}_{\Omega_{2r}}|P|^2dx
+c\Big({-\hspace{-0.38cm}\int}_{\Omega_{2r}}|\mathbf{F}|^{2(1+\sigma)}dx
\Big)^{\frac{1}{1+\sigma}},
\end{aligned}
\end{equation}
where $c=c(n,\gamma_1,\gamma_2,\nu,\Lambda,\Omega)>0$.
\end{lemma}
\begin{proof}
Without loss of generality, for any fixed boundary point $y \in \partial\Omega$
we set $\Omega_{2r}=\Omega_{2r}(y)$.
By a zero extension of $u$ in $B_{2r}(y)\setminus \Omega_{2r}$, we can take
$v=\eta^2(u-(u)_{2r})$
as a test function in the neighbourhood of boundary point. Using the arguments
analogous to the proof of Lemma \ref{reverse-Holder-a}
due to the measure density property of $\Omega$ (cf. \eqref{thick-b}),
we can obtain the estimate \eqref{reverse-Holder-ineq-b} immediately.
\end{proof}
Here, let us recall a basic property that the generalized Stokes problem
\eqref{Stokes-systems} is invariant under
scaling transformation and normalization, see \cite[Lemma 2.6]{ByS17}.
\begin{lemma}\label{scaling-normalization}
For fixed $K>1$ and $0<\rho<1$, let
\[
\tilde{\mathbf{A}}(x)=\mathbf{A}(\rho x),\quad
\tilde u(x)=\frac{u(\rho x)}{K\rho}, \quad \tilde P(x)=\frac{P(\rho x)}{K}, \quad
\tilde{\mathbf{F}}(x)=\frac{\mathbf{F}(\rho x)}{K}
\]
and let $\tilde\Omega=\{\frac{x}{\rho}: x\in \Omega\}$.
Then the following three statements are true.
\begin{itemize}
\item[(i)] If $(u,P)$ is a weak solution pair to system \eqref{Stokes-systems},
then $(\tilde u,\tilde P)$ is a weak solution pair to
\begin{gather*}
\operatorname{div}(\tilde{\mathbf{A}}(x)\nabla\tilde u)-\nabla\tilde P
= \operatorname{div} \tilde{\mathbf{F}}, \quad \text{in } \tilde\Omega \\
\operatorname{div} \tilde u=0, \quad \text{in }\tilde\Omega,\\
\tilde u=0, \quad \text{on } \partial\tilde\Omega.
\end{gather*}
\item[(ii)] If $\mathbf{A}$ satisfies condition \eqref{uniform-bound-ellipticity},
then so dose $\tilde{\mathbf{A}}$ with the same constants $\nu$ and $\Lambda$.
\item[(iii)] If $(\mathbf{A},\Omega)$ is $(\delta,R_0)$-vanishing, then
$(\tilde{\mathbf{A}},\tilde\Omega)$ is
$\big(\delta,\frac{R_0}{\rho}\big)$-vanishing.
\end{itemize}
\end{lemma}
Let us consider the comparison estimates at the interior point and boundary point.
For simplicity, we fix $x_i, z_i \in \Omega$ and let $r_i< \frac{R}{250}$ for
each $i$. Set
\begin{equation}\label{R-1}
R \leq \min \big\{\frac{R_0}{2}, \frac{R_0}{c^*}, 1\big\},
\end{equation}
where $c^*=c^*(n,\gamma_1,\gamma_2,\nu,\Lambda,\omega(\cdot),|\Omega|))
\geq |\Omega|+1$.
Let
\begin{gather*}
B_i^0=B_{r_i}(x_i), \quad B_i^j=B_{5jr_i}(x_i), \quad j=1,\dots,6;\\
\Omega_i^0=\Omega_{r_i}(z_i), \quad \Omega_i^j=\Omega_{25jr_i}(0),
\quad j=1,\dots ,6.
\end{gather*}
For a boundary point $ y_i \in B_{30r_i}(x_i) \cap \partial\Omega$,
there exists a new coordinate system in
$z_i =( z^1_i,z^2_i,\dots,z^n_i)$-variables such that
\begin{equation}\label{coordinate-system}
\begin{gathered}
\ z_i=x_i; \quad y_i +150\delta r_i(0,\dots,0,1) \quad \text{is the origin}, \\
B_{150r_i}^+(0) \subset \Omega_{150r_i}(0) \subset B_{150r_i}(0)
\cap \{z^n_i>-300\delta r_i\}.
\end{gathered}
\end{equation}
Choose $0<\delta <1/300$ such that
$\Omega_{r_i}(z_i) \subset \Omega_{40r_i}(0)$
in the new $z$-coordinate system. Then we have
\begin{gather}\label{B-omega}
B_i^{j+}=B_{25jr}^{+}(0) \subset \Omega_i^j \subset B_{25jr}(0)
\cap \{y_n>-300\delta r\}, \\
\label{omega-a}
\Omega_i^j \subset \Omega_{150r_i}(0) \subset \Omega_{190r_i}(z_i) .
\end{gather}
By a scaling transformation to \eqref{B-omega}, without loss of generality,
we set
\[
B_6^+ \subset \Omega_6 \subset B_6 \cap \{x_n>-12\delta\}.
\]
We now consider a series of localizing problems in the neighbourhood of
boundary point as follows:
\begin{equation}\label{comparison-equ-1}
\begin{gathered}
\operatorname{div}(\mathbf{A}(x)\nabla u)-\nabla P= \operatorname{div} \mathbf{F},
\quad \text{in } \Omega_6,\\
\operatorname{div} u=0, \quad \text{in } \Omega_6,\\
u=0, \quad \text{on } \partial_\omega\Omega_6;
\end{gathered}
\end{equation}
\begin{equation}\label{comparison-equ-2}
\begin{gathered}
\operatorname{div}(\mathbf{A}(x)\nabla v)-\nabla P_v= 0, \quad
\text{in } \Omega_5,\\
\operatorname{div} v=0, \quad \text{in } \Omega_5,\\
v=u, \quad \text{on }\partial\Omega_5;
\end{gathered}
\end{equation}
\begin{equation}\label{comparison-equ-3}
\begin{gathered}
\operatorname{div}(\mathbf{A}_{B_4^+}\nabla w)-\nabla P_w= 0, \quad
\text{in } \Omega_4,\\
\operatorname{div} w=0, \quad \text{in } \Omega_4,\\
w=v, \quad \text{on } \partial\Omega_4;
\end{gathered}
\end{equation}
and
\begin{equation}\label{comparison-equ-4}
\begin{gathered}
\operatorname{div}(\mathbf{A}_{B_4^+}\nabla h)-\nabla P_h= 0, \quad
\text{in $B_4^+$},\\
\operatorname{div} h=0, \quad \text{in } B_4^+,\\
h=w, \quad \text{on } T_4.
\end{gathered}
\end{equation}
Let us recall an approximating estimate in accordance with Byun-So's work,
see \cite[Lemma 3.6]{ByS17}.
\begin{lemma}\label{limiting-problem}
For any $0< \epsilon <1$, if $(w,P_w)$ is a weak solution pair of
problem \eqref{comparison-equ-3} and satisfies
\[
{-\hspace{-0.38cm}\int}_{\Omega_4}\big(|\nabla w|^2+|P_w|^2\big)dx \leq c.
\]
Then, there exists a weak solution pair $(h,P_h) $ of
problem \eqref{comparison-equ-4} with
\begin{gather*}
{-\hspace{-0.38cm}\int}_{B_4^+}\left(|\nabla h|^2+|P_h|^2\right)dx \leq c,\\
{-\hspace{-0.38cm}\int}_{B_4^+}|w-h|^2dx \leq \epsilon^2.
\end{gather*}
\end{lemma}
A local Lipschitz regularity in the neighbourhood of any boundary point
for the Dirichlet problem \eqref{comparison-equ-4}
is well-established, see \cite[Lemma 3.5]{ByS17}.
\begin{lemma}\label{boundary-bound-estimate-a}
Let $(h,P_h)$ be the weak solution pair of problem \eqref{comparison-equ-4}.
Then we have
\begin{gather*}
\|\nabla \bar h\|_{L^\infty(\Omega_3)^{n^2}}
=\|\nabla h\|_{L^\infty(B_{3}^{+})^{n^2}}
\leq c\|\nabla h\|_{L^2(B_{4}^{+})^{n^2}}, \\
\|P_{\bar h}\|_{L^\infty(\Omega_3)}
=\|P_{ h}\|_{L^\infty(B_{3}^{+})}
\leq c\Big(\|\nabla h\|_{L^2(B_{4}^{+})^{n^2}}+\|P_h\|_{L^2(B_4^+)}\Big),
\end{gather*}
where $\bar h$ is the zero extension of $h$ from $B_3^{+}$ to $\Omega_3$,
and $P_{\bar h}$ is an associated pressure of $\bar h$ by extending it
from $B_3^{+}$ to $\Omega_3$.
\end{lemma}
\begin{lemma}\label{comparison-1}
Let $(u,P)$ be a weak solution pair of problem \eqref{comparison-equ-1}.
If for any $0< \epsilon <1$ there exists a constant
$\delta=\delta(\epsilon,\gamma_1,\gamma_2)$
such that
\begin{equation}\label{CZ-operator-estimate-a}
{-\hspace{-0.38cm}\int}_{\Omega_5}\big(|\nabla u|^2+|P|^2\big) dx \leq 1, \quad
{-\hspace{-0.38cm}\int}_{\Omega_5}|{\mathbf{F}}|^2dx
\leq\delta^{\gamma_1/\gamma_2}, \quad
{-\hspace{-0.38cm}\int}_{\Omega_6}|\mathbf{A}-\mathbf{A}_{\Omega_6}| ^2\,dx
\leq \delta^2,
\end{equation}
then there exists a weak solution pair $(h,P_h)$ of
problem \eqref{comparison-equ-4} such that
\begin{gather*}
{-\hspace{-0.38cm}\int}_{B_4^+}\big(|\nabla h|^2+|P_h|^2\big)dx \leq c, \\
{-\hspace{-0.38cm}\int}_{\Omega_3}
\big(|\nabla u-\nabla\bar h|^2+|P-P_{\bar h}|^2\big)dx \leq \epsilon^2,
\end{gather*}
where $\bar h$ and $P_{\bar h}$ are the same as given in
Lemma \ref{boundary-bound-estimate-a}.
\end{lemma}
\begin{proof}
Let $(v,P_v)$ and $(w,P_w)$ be weak solution pairs to
problems \eqref{comparison-equ-2} and \eqref{comparison-equ-3}, respectively.
It follows Lemma \ref{existence-energy-estimate} and
\eqref{CZ-operator-estimate-a} that
\begin{equation}\label{u-v-a}
{-\hspace{-0.38cm}\int}_{\Omega_5}\big(|\nabla u- \nabla v|^2+|P-P_v|^2\big)dx
\leq c{-\hspace{-0.38cm}\int}_{\Omega_5}|{\mathbf{F}}|^2dx
\leq c\delta^{\gamma_1/\gamma_2}.
\end{equation}
Combining \eqref{comparison-equ-3} and \eqref{comparison-equ-2} leads to
\begin{gather*}
\operatorname{div}(\mathbf{A}_{B_4^+}\nabla (v- w))-\nabla (P_v-P_w)
= -\operatorname{div}(\mathbf{A}(x)-\mathbf{A}_{B_4^+})\nabla v), \quad
\text{in } \Omega_4,\\
\operatorname{div} v-w=0, \quad \text{in } \Omega_4,\\
v-w=0, \quad \text{on } \partial\Omega_4.
\end{gather*}
In view of higher integrability on the boundary
and normalization conditions \eqref{CZ-operator-estimate-a}, similar to the
proof of \cite[Lemma 3.7]{ByS17},
using \eqref{existence-energy-estimate-a} we obtain
\begin{equation}\label{v-w-a}
{-\hspace{-0.38cm}\int}_{\Omega_4}
\Big(|\nabla v- \nabla w|^2+|P_v-P_w|^2\Big)dx
\leq c\big(\delta^2+\delta^3\big)^{\frac{r_1-2}{r_1}} ,
\end{equation}
where $r_1>2$ is the same as given in Lemma \ref{reverse-Holder-b}.
Combining \eqref{u-v-a} and \eqref{v-w-a} leads to
\begin{equation}\label{u-w-a}
{-\hspace{-0.38cm}\int}_{\Omega_4}\Big(|\nabla u- \nabla w|^2+|P-P_w|^2\Big)dx
\leq c\big(\delta^{\gamma_1/\gamma_2}+\delta^{2-\frac{4}{r_1}}
+\delta^{3-\frac{6}{r_1}}\big),
\end{equation}
from which and \eqref{CZ-operator-estimate-a} we obtain
\[
{-\hspace{-0.38cm}\int}_{\Omega_4}\Big(|\nabla w|^2+|P_w|^2\Big)dx \leq c.
\]
It follows Lemma \ref{limiting-problem} that there exists a weak solution
pair $(h,P_h)$ of problem \eqref{comparison-equ-4} with
\[
{-\hspace{-0.38cm}\int}_{B_4^+}\Big(|\nabla h|^2+|P_h|^2\Big)dx \leq c, \quad
{-\hspace{-0.38cm}\int}_{B_4^+}|w-h|^2dx \leq \epsilon_*^2,
\]
where $\epsilon_*^2$ is to be determined later.
Notice that $(\bar h,P_{\bar h})$ is a weak solution pair of
\begin{gather*}
\operatorname{div}(\mathbf{A}_{B_4^+}\nabla \bar h)-\nabla P_{\bar h}
= -\frac{\partial}{\partial x_n}
\Big(\bar a_{nn}^{\alpha\beta}\frac{\partial h^\alpha}{\partial x_n}(x',0)
\chi_{\mathbb{R}^n_-}(x)\Big),
\quad \text{in } \Omega_4,\\
\operatorname{div} \bar h=0, \quad \text{in } \Omega_4,\\
\bar h=0, \quad \text{on } \partial_\omega\Omega_4,
\end{gather*}
where $\mathbf{A}_{B_4^+}=\bar a_{ij}^{\alpha\beta}$,
$x'=(x_1,x_2,\dots,x_{n-1})$ and $\chi$ is a standard characteristic function.
Processing in a similar manner to the proof of \cite[Lemma 3.7]{ByS17},
we deduce that
\begin{equation}\label{w-h-a-3}
{-\hspace{-0.38cm}\int}_{\Omega_2}|\nabla w-\nabla \bar h|^2dx
\leq \epsilon_*^2+c\delta+c\delta^{\frac{2}{n}},
\end{equation}
and
\begin{equation}\label{p-w-h-a}
\begin{aligned}
&{-\hspace{-0.38cm}\int}_{\Omega_2}|P_w-P_{\bar h}|^2\,dx\\
& \leq c{-\hspace{-0.38cm}\int}_{\Omega_2}|\nabla w-\nabla \bar h|^2\,dx
+c{-\hspace{-0.38cm}\int}_{\Omega_2}\big|\bar a_{nn}^{\alpha\beta}
\frac{\partial h^\alpha}{\partial x_n}(x',0)\chi_{\mathbb{R}^n_-}(x)\big|^2\,dx \\
& \leq c\big(\epsilon_*^2+\delta+\delta^{\frac{2}{n}}\big).
\end{aligned}
\end{equation}
Taking $\epsilon_*, \delta >0$ small enough and making use of
\eqref{u-w-a}-\eqref{p-w-h-a} such that
\[
c\Big(\delta^\frac{\gamma_1}{\gamma_2}+\delta^{2-\frac{4}{r_1}}
+\delta^{3-\frac{6}{r_1}}+\epsilon_*^2+\delta
+\delta^{\frac{2}{n}}\Big)<\epsilon^2,
\]
consequently, by \eqref{w-h-a-3} and \eqref{p-w-h-a} we arrive at the desired result.
\end{proof}
For the interior case, one can process in an analogous but simple way as
for the boundary case. Similar to Lemmas \ref{boundary-bound-estimate-a}
and \ref{comparison-1}, we replace $\Omega_j$ by $B_{j-1}$ with
$B_6\Subset \Omega$. The first one is a local Lipschitz regularity
of $\nabla w$ to the Dirichlet problem for a limiting system of local
constant coefficients, and the second one is a comparison between
limiting system and the local version of system \eqref{Stokes-systems}
under the normalization. For details, see
\cite[Lemmas 3.1 and 3.2]{ByS17}.
\begin{lemma}\label{interior-bound-estimate-b}
Suppose that $(w,P_w)$ is a weak solution pair to
\begin{equation}\label{interior-comparison-equ-3}
\begin{gathered}
\operatorname{div}(\mathbf{A}_{B_3}\nabla w)-\nabla P_w= 0, \quad \text{in $B_3$},\\
\operatorname{div} w=0, \quad \text{in $B_3$}\\
w=v, \quad \text{on } \partial B_3,
\end{gathered}
\end{equation}
and $u$ is a local weak solution of system \eqref{Stokes-systems} with
\[
{-\hspace{-0.38cm}\int}_{B_4}\big(|\nabla u|^2+|P|^2\big)dx\leq 1 \quad
\text{and } {-\hspace{-0.38cm}\int}_{B_4}|{\mathbf{F}}|^2dx
\leq \delta^{\gamma_1/\gamma_2}.
\]
Then
\[
\|\nabla w\|_{L^\infty(B_2)^{n^2}}+ \|P_w\|_{L^\infty(B_2)} \leq c_2,
\]
where $c_2>1$.
\end{lemma}
\begin{lemma}\label{interior-bound-estimate-a}
Let $u$ be a weak solution of \eqref{Stokes-systems}.
If for any $0< \epsilon <1$, there exists a constant
$\delta=\delta(\epsilon,\gamma_1,\gamma_2)$ with
\[
{-\hspace{-0.38cm}\int}_{B_4}|{\mathbf{F}}|^2dx
\leq \delta^{\gamma_1/\gamma_2} \quad \text{and} \quad
{-\hspace{-0.38cm}\int}_{B_4}|\mathbf{A}(x)-\mathbf{A}_{B_4}| ^2\,dx
\leq \delta^2,
\]
then there exists a weak solution pair $(w,P_w)$ to \eqref{interior-comparison-equ-3}
such that
\[
{-\hspace{-0.38cm}\int}_{B_3}\left(|\nabla u-\nabla w|^2+|P-P_{w}|^2\right)dx
\leq \epsilon^2.
\]
\end{lemma}
\begin{proof}
Let $(v,P_v)$ be a weak solution pair to \eqref{comparison-equ-2} in $B_4$.
From \eqref{comparison-equ-2} and \eqref{comparison-equ-1} it follows that
\begin{gather*}
\operatorname{div}(\mathbf{A}(x)\nabla (u-v))-\nabla (P-P_v)
= \operatorname{div} \mathbf{F}, \quad \text{in } B_4,\\
\operatorname{div} (u-v)=0, \quad \text{in } B_4 \\
u-v=0, \quad \text{on } \partial B_4.
\end{gather*}
By Lemma \ref{existence-energy-estimate} and inequality
\eqref{uniform-bound-ellipticity}, and H\"older's inequality it follows that
\begin{equation}\label{interior-u-v}
{-\hspace{-0.38cm}\int}_{B_4}\big(|\nabla (u-v)|^2+|P-P_v|^2\big)dx
\leq {-\hspace{-0.38cm}\int}_{B_4}|\mathbf{F}|^2\,dx
\leq c\delta^{\gamma_1/\gamma_2}.
\end{equation}
Similarly, let $(w,P_w)$ be a weak solution pair to \eqref{interior-comparison-equ-3}.
Regarding \eqref{comparison-equ-2} in $B_4$, we have
\begin{gather*}
\operatorname{div}(\mathbf{A}_{B_3}\nabla (v-w))-\nabla (P_v-P_w)
=- \operatorname{div} \big((\mathbf{A}(x)-\mathbf{A}_{B_3})\nabla v\big), \quad
\text{in } B_3,\\
\operatorname{div} (v-w)=0, \quad \text{in } B_3,\\
v-w=0, \quad \text{on } \partial B_3.
\end{gather*}
Processing in a way analogous to the proof of \cite[Lemma 3.3]{ByS17}, we obtain
\begin{equation}\label{interior-v-w}
{-\hspace{-0.38cm}\int}_{B_3}\left(|\nabla (v-w)|^2+|P_v-P_w|^2\right)dx
\leq c\delta^{2-\frac{4}{r_2}}.
\end{equation}
Combining \eqref{interior-u-v} and \eqref{interior-v-w} gives
\[
{-\hspace{-0.38cm}\int}_{B_3}\left(|\nabla u-\nabla w|^2+|P-P_w|^2\right)dx
\leq c \big( \delta^{\gamma_1/\gamma_2}+\delta^{2-\frac{4}{r_2}} \big).
\]
Taking $\delta>0$ small enough such that
$c (\delta^{\gamma_1/\gamma_2}+\delta^{2-\frac{4}{r_2}})
=\epsilon^2$, hence we obtain the desired result.
\end{proof}
Let us recall the embedding relation with respect to the Lorentz spaces,
see \cite[Proposition 3.9]{MeP12}.
\begin{proposition}\label{preliminary-results}
Let $U$ be a bounded measurable subset of $\mathbb{R}^n$. Then the following
three statements are true.
\begin{itemize}
\item[(i)] If $00$ it holds
\[
\int_0^\infty \lambda^r \Big(\int_\lambda^\infty f(\mu) d\mu \Big)^\alpha
\frac{d\lambda}{\lambda}
\leq \big(\frac{\alpha}{r}\big)^\alpha\int_0^\infty \lambda^r
\Big(\lambda f(\lambda) \Big)^\alpha \frac{d\lambda}{\lambda}.
\]
\end{lemma}
\begin{lemma}\label{reverse-Holder-c}
Let $h: [0,+\infty)\to [0,+\infty)$ be a non-increasing measurable function.
Suppose that $\alpha_1\leq \alpha_2$ and $r>0$. If $\alpha_2< \infty$, then
\[
\Big(\int_\lambda^\infty \big(\mu^r h(\mu) \big)^{\alpha_2}
\frac{d\mu}{\mu}\Big)^{1/\alpha_2}
\leq \varepsilon \lambda^r h(\lambda)
+\frac{c}{\varepsilon^{\frac{\alpha_2}{\alpha_1-1}}}
\Big(\int_\lambda^\infty \big(\mu^r h(\mu)\big)^{\alpha_1} \frac{d\mu}{\mu}
\Big)^{1/\alpha_1}
\]
for $\varepsilon \in (0,1]$ and $\lambda\geq 0$. If $\alpha_2=\infty$, then
\[
\sup_{\mu> \lambda} \left(\mu^r h(\mu) \right)
\leq c\lambda^r h(\lambda)
+c\Big(\int_\lambda^\infty \left(\mu^r h(\mu)\right)^{\alpha_1}
\frac{d\mu}{\mu}\Big)^{1/\alpha_1},
\]
where $c$ depends on $\alpha_1,\, \alpha_2$ and $r$.
\end{lemma}
The following lemma is regarding an iteration argument,
see \cite {Gia83} or \cite[Lemma 4.1]{Min11}.
\begin{lemma}\label{iteration-argument}
Let $\varphi:[r_1,2r_1]\to [0,\infty)$ be a function such that
\[
\varphi(\rho_1)\leq \frac{1}{2}\varphi(\rho_2)+B_0(\rho_2-\rho_1)^{-\beta}+L
\]
for $r_1<\rho_1<\rho_2< 2r_1$, where $B_0,\,L\geq 0$ and $\beta>0$. Then we have
\[
\varphi(r_1)\leq c(\beta)B_0r_1^{-\beta}+cL.
\]
\end{lemma}
\section{Proof of Theorem \ref{main-result}}
With aid of the lemmas presented in the preceding section,
now we are in a position to prove Theorem \ref{main-result} by means of the
large-M-inequality principle \cite{AcM07} and a geometric argument \cite{ByOW14}.
To present our discussion in a straightforward and lucid manner, we
separate our proof into six steps.
\begin{proof}
We only treat the boundary case. For the interior case, one can process in
a similar but much simpler way.
For the boundary case, a proper translation and rotation of the original
coordinates does not change the corresponding features.
Without loss of generality, we may assume
that $R_0\leq 1$. For any fixed $x_0\in \Omega$, we set
\begin{gather*}
p^-=\inf_{\Omega_{2R}(x_0)} p(x), \quad
p^+=\sup_{\Omega_{2R}(x_0)} p(x), \\
p_i^-=\inf_{\Omega_i^5(x_0)} p(x), \quad
p_i^+=\sup_{\Omega_i^5(x_0)} p(x).
\end{gather*}
\smallskip
\noindent\textbf{Step 1.}
In this step, we present a modified Vitali's covering. Let $u$ be the weak
solution of system \eqref{Stokes-systems}.
For $\Omega_{R}=\Omega_{R}(x_0)$, we define
\begin{equation}\label{lambda-0}
\lambda_0:={-\hspace{-0.38cm}\int}_{\Omega_{2R}}
(|\nabla u|+|P|)^{2p(x)/p^-}\,dx
+ \frac{1}{\delta}\Big({-\hspace{-0.38cm}\int}_{\Omega_{2R}}
\Big(|\mathbf{F}|^{2p(x)/p^-}+1\Big)^\eta \,dx\Big)^{1/\eta},
\end{equation}
where $\delta >0 $ and $ \eta>1$ will be specified later.
We now introduce the super-level set
\[
E\left(\lambda,\Omega_{R}\right)
:=\big\{x\in \Omega_{R}: (|\nabla u|+|P|)^{2p(x)/p^-}>\lambda
\big\}
\]
for some $\lambda>M\lambda_0 \geq 1 $ with $M=(\frac{8000}{7})^n$.
For $ x_i\in E\left(\lambda,\Omega_{R}\right)$ and radii $0 \lambda$.
Thus, in view of absolute continuity of the integral, we can pick a maximal
radius $r_i=r_{x_i}$ such that
\begin{equation}\label{equ-lambda}
\begin{aligned}
&CZ(\Omega_{r_i}( x_i))\\
&= {-\hspace{-0.38cm}\int}_{\Omega_{r_i}( x_i)}
(|\nabla u|+|P|)^{2p(x)/p^-}\,dx
+ \frac{1}{\delta}\Big({-\hspace{-0.38cm}\int}_{\Omega_{r_i}(x_i)}
|\mathbf{F}|^{\frac{2p(x)}{p^-}\eta}\,dx\Big)^{1/\eta}
=\lambda
\end{aligned}
\end{equation}
for any point $x_i\in E\left(\lambda,\Omega_{R}\right)$.
Moreover, for $r \in (r_i, R]$ one has
\begin{equation}\label{CZ-operator-estimate-c}
CZ(\Omega_{r}( x_i))< \lambda.
\end{equation}
From \eqref{equ-lambda}, we obtain the following alternatives:
\begin{equation} \label{alternatives}
\frac{\lambda}{2} \leq {-\hspace{-0.38cm}\int}_{\Omega_{r_i}( x_i)}
(|\nabla u|+|P|)^{2p(x)/p^-}\,dx\quad \text{or} \quad
\left(\frac{\delta\lambda}{2}\right)^{\eta}
\leq {-\hspace{-0.38cm}\int}_{\Omega_{r_i}( x_i)}|\mathbf{F}|^{\frac{2p(x)}{p^-}\eta}
\,dx.
\end{equation}
Suppose that the first case of \eqref{alternatives} is valid,
and split the integral as
\begin{align*}
&{-\hspace{-0.38cm}\int}_{\Omega_{r_i}( x_i)}
(|\nabla u|+|P|)^{2p(x)/p^-}\,dx \\
& \leq \frac{|\Omega_{r_i}( x_i)\backslash E(\frac{\lambda}{4},\Omega_{2R})|}
{|\Omega_{r_i}( x_i)|}{-\hspace{-0.38cm}\int}_{\Omega_{r_i}( x_i)
\backslash E(\frac{\lambda}{4},\Omega_{2R})}
(|\nabla u|+|P|)^{2p(x)/p^-}\,dx\\
&\quad + \frac{1}{|\Omega_{r_i}( x_i)|}
\int_{\Omega_{r_i}( x_i)\cap E(\frac{\lambda}{4},\Omega_{2R})}
(|\nabla u|+|P|)^{2p(x)/p^-}\,dx\\
&\leq \frac{\lambda}{4}+c\frac{|\Omega_{r_i}( x_i)
\cap E(\frac{\lambda}{4},\Omega_{2R})|}{|\Omega_{r_i}( x_i)
\cap E(\frac{\lambda}{4},\Omega_{2R})|^{\frac{1}{1+\sigma_1}}}
\frac{1}{|\Omega_{r_i}( x_i)|} \\
&\quad\times \Big(\int_{\Omega_{r_i}( x_i)}(|\nabla u|+|P|)^{\frac{2p(x)}{p^-}
(1+\sigma_1)}\,dx\Big)^{\frac{1}{1+\sigma_1}}\\
&\leq \frac{\lambda}{4}+c\Big(\frac{|\Omega_{r_i}( x_i)
\cap E(\frac{\lambda}{4},\Omega_{2R})|}{|\Omega_{r_i}( x_i)|}
\Big)^{1-\frac{1}{1+\sigma_1}}
\Big({-\hspace{-0.38cm}\int}_{\Omega_{r_i}( x_i)}
(|\nabla u|+|P|)^{\frac{2p(x)}{p^-}(1+\sigma_1)}\,dx
\Big)^{\frac{1}{1+\sigma_1}}.
\end{align*}
Take
\[
0<\sigma_1 \leq \frac{\gamma_1(1+\sigma)}{\gamma_1+\omega(2R)}-1,
\]
where $\sigma$ is the same as in Lemma \ref{reverse-Holder-b}. Then
\[
\frac{p(x)}{p^-}(1+\sigma_1)
=\Big(1+\frac{p(x)-p^-}{p^-}\Big)(1+\sigma_1)
\leq \Big(1+\frac{\omega(2R)}{p^-}\Big)(1+\sigma_1)
\leq (1+\sigma).
\]
Let $\eta=1+\sigma_1$. In view of \eqref{CZ-operator-estimate-c}, it follows by
H\"older's inequality (given in Lemma \ref{reverse-Holder-b})
that
\begin{align*}
& \Big({-\hspace{-0.38cm}\int}_{\Omega_{r_i}( x_i)}
(|\nabla u|+|P|)^{\frac{2p(x)}{p^-}(1+\sigma_1)}\,dx
\Big)^{\frac{1}{1+\sigma_1}}\\
& \leq c \Big[{-\hspace{-0.38cm}\int}_{\Omega_{2r_i}( x_i)}|\nabla u|^{2p(x)/p^-}
\,dx+{-\hspace{-0.38cm}\int}_{\Omega_{2r_i}( x_i)}|P|^{2p(x)/p^-}\,dx \\
&\quad + \Big({-\hspace{-0.38cm}\int}_{\Omega_{2r_i}( x_i)}|\mathbf{F}|^{\frac{2p(x)}{p^-}
(1+\sigma_1)}\,dx \Big)^{\frac{1}{1+\sigma_1}}+1\Big]\\
&\leq c \lambda.
\end{align*}
Hence, using \eqref{alternatives} in combination with
\[
\frac{\lambda}{4}
\leq c\Big(\frac{|\Omega_{r_i}( x_i)\cap E(\frac{\lambda}{4},\Omega_{2R})|}
{|\Omega_{r_i}( x_i)|}\Big)^{1-\frac{1}{1+\sigma_1}}\lambda
\]
we have
\begin{equation}\label{omega-estimate-a}
|\Omega_{r_i}( x_i)|\leq c \big|\Omega_{r_i}( x_i)\cap
E\big(\frac{\lambda}{4},\Omega_{2R}\big) \big|,
\end{equation}
where the constant $c$ depends on $n,\gamma_2,\gamma_2,\nu,\Lambda$ and $t$.
For the second estimate in \eqref{alternatives}, it follows from Fubini's theorem
that
\begin{align*}
\Big(\frac{\lambda \delta}{2}\Big)^{\eta}
& \leq {-\hspace{-0.38cm}\int}_{\Omega_{r_i}( x_i)}
|\mathbf{F}|^{\frac{2p(x)}{p^-}\eta}dx\\
&= \frac{\eta}{|\Omega_{r_i}( x_i)|}\int_0^{\infty}\mu^{\eta}
\big|\big\{x\in \Omega_{r_i}( x_i): |\mathbf{F}|^{2p(x)/p^-}>\mu
\big\}\big|\frac{d\mu}{\mu}\\
&= \frac{\eta}{|\Omega_{r_i}( x_i)|}\int_0^{\zeta\lambda}\mu^{\eta}
\big|\big\{x\in \Omega_{r_i}( x_i): |\mathbf{F}|^{2p(x)/p^-}>\mu
\big\}\big|\frac{d\mu}{\mu}\\
&\quad +\frac{\eta}{|\Omega_{r_i}( x_i)|}\int_{\zeta\lambda}^{\infty}\mu^{\eta}
\big|\big\{x\in \Omega_{r_i}( x_i): |\mathbf{F}|^{2p(x)/p^-}>\mu
\big\}\big|\frac{d\mu}{\mu}\\
&\leq (\zeta\lambda)^{\eta}+ \frac{\eta}{|\Omega_{r_i}( x_i)|}
\int_{\zeta\lambda}^{\infty}\mu^{\eta}
\big|\big\{x\in \Omega_{r_i}( x_i): |\mathbf{F}|^{2p(x)/p^-}>\mu
\big\}\big|\frac{d\mu}{\mu}.
\end{align*}
Let $\delta=4 \zeta$. We derive that
\[
(\zeta \lambda )^{\eta}
\leq \frac{\eta}{|\Omega_{r_i}( x_i)|}\int_{\zeta\lambda}^{\infty}\mu^{\eta}
\big|\big\{x\in \Omega_{r_i}( x_i): |\mathbf{F}|^{2p(x)/p^-}>\mu\big\}\big|
\frac{d\mu}{\mu}
\]
and
\begin{equation} \label{omega-estimate-b}
|\Omega_{r_i}( x_i) | \leq \frac{\eta}{(\zeta \lambda)^{\eta}}
\int_{\zeta\lambda}^{\infty}\mu^{\eta}
\big|\big\{x\in \Omega_{r_i}( x_i): |\mathbf{F}|^{2p(x)/p^-}>\mu\big\}\big|
\frac{d\mu}{\mu}.
\end{equation}
Combining \eqref{omega-estimate-a} and \eqref{omega-estimate-b} we have
\begin{equation}\label{omega-estimate-c}
\begin{aligned}
|\Omega_{r_i}( x_i)|
&\leq c \big|\Omega_{r_i}( x_i)\cap E\big(\frac{\lambda}{4},\Omega_{2R} \big) \big|\\
&\quad + \frac{c\eta}{(\zeta \lambda)^{\eta}}\int_{\zeta\lambda}^{\infty}\mu^{\eta}
\big|\big\{x\in \Omega_{r_i}( x_i): |\mathbf{F}|^{2p(x)/p^-}>\mu\big\}\big|
\frac{d\mu}{\mu}.
\end{aligned}
\end{equation}
Note that $\Omega$ is $(\delta,R_0)$-Reifenberg flat.
It follows from \eqref{omega-a} that
\begin{align*}
&{-\hspace{-0.38cm}\int}_{\Omega_i^j}
(|\nabla u|+|P|)^{2p(z)/p^-}dz
+\frac{1}{\delta}\Big({-\hspace{-0.38cm}\int}_{\Omega_i^j}\mathbf{F}|^{\frac{2p(z)}
{p^-}\eta}dz\Big)^{1/\eta} \\
&\leq \frac{|\Omega_{250r_i}(z_i)|}{|\Omega_{25r_i}(0)|}
\Big[{-\hspace{-0.38cm}\int}_{\Omega_{250r_i}(z_i)}
(|\nabla u|+|P|)^{2p(z)/p^-}dz
+ \Big({-\hspace{-0.38cm}\int}_{\Omega_{250r_i}(z_i)}
|\mathbf{F}|^{\frac{2p(z)}{p^-}\eta} dz\Big)^{1/\eta}\Big]\\
& \leq \frac{|B_{250r_i}(z_i)|}{|B^+_{25r_i}(0)|}
\Big[{-\hspace{-0.38cm}\int}_{\Omega_{250r_i}(z_i)}
(|\nabla u|+|P|)^{2p(z)/p^-}dz
+\frac{1}{\delta}\Big( {-\hspace{-0.38cm}\int}_{\Omega_{250r_i}(z_i)}
|\mathbf{F}|^{\frac{2p(z)}{p^-}\eta} dz\Big)^{1/\eta}\Big]\\
& \leq 2 \cdot 10^n
\Big[ {-\hspace{-0.38cm}\int}_{\Omega_{250r_i}(z_i)}
(|\nabla u|+|P|)^{2p(z)/p^-}dz+ \frac{1}{\delta}
\Big({-\hspace{-0.38cm}\int}_{\Omega_{250r_i}(z_i)}
|\mathbf{F}|^{\frac{2p(z)}{p^-}\eta} dz\Big)^{1/\eta}\Big].
\end{align*}
Employing \eqref{CZ-operator-estimate-c} and using change of variables
of \eqref{coordinate-system}, we deduce that
\begin{equation}\label{CZ-operator-estimate-boundary}
{-\hspace{-0.38cm}\int}_{\Omega_i^j} (|\nabla u|+|P|)^{2p(z)/p^-}dz
+\frac{1}{\delta}\Big({-\hspace{-0.38cm}\int}_{\Omega_i^j}
|\mathbf{F}|^{\frac{2p(z)}{p^-}\eta} dz\Big)^{1/\eta} \leq 2 \cdot 10^n \lambda.
\end{equation}
\smallskip
\noindent\textbf{Step 2.}
We consider various comparison estimates. Since $(\mathbf{A},\Omega)$ is
$(\delta,R_0)$-vanishing for some $R_0>0$, we have
\begin{equation}\label{Reifenberg-flat-a}
{-\hspace{-0.38cm}\int}_{\Omega_{150r_i}(0)}
|\mathbf{A}-\mathbf{A}_{\Omega_{150r_i}(0)}|^2dx \leq \delta^2.
\end{equation}
Taking \eqref{CZ-operator-estimate-boundary} into account, we have
\begin{equation}\label{CZ-operator-estimate-d}
\begin{gathered}
{-\hspace{-0.38cm}\int}_{\Omega_i^5}
(|\nabla u|+|P|)^{2p(x)/p^-}\,dx \leq c\lambda, \\
\Big( {-\hspace{-0.38cm}\int}_{\Omega_i^5}|\mathbf{F}|^{\frac{2p(z)}{p^-}\eta}
\,dx\Big)^{1/\eta} \leq c\lambda\delta.
\end{gathered}
\end{equation}
Let us first show that
\begin{equation}\label{CZ-operator-estimate-e}
\begin{gathered}
{-\hspace{-0.38cm}\int}_{\Omega_i^5}\big(|\nabla u|^2+|P|^2\big) \,dx
\leq c_3\lambda^{p^-/p_i^+}, \\
{-\hspace{-0.38cm}\int}_{\Omega_i^5}|\mathbf{F}|^2\,dx
\leq c_3\lambda^{p^-/p_i^+}\delta^{\gamma_1/\gamma_2},
\end{gathered}
\end{equation}
for a constant $c_3\geq 1$. We claim that
\begin{equation}\label{u-p-2-a}
\Big({-\hspace{-0.38cm}\int}_{\Omega_i^5}\big(|\nabla u|^2+|P|^2\big) dx
\Big)^{p_i^+-p_i^-} \leq c,
\end{equation}
where $c\geq 1$ is a universal constant.
Notice that $|\mathbf{F}|^{p(x)}\in L(t,q)(\Omega)$ implies
$\int_\Omega |\mathbf{F}|^{p(x)}\,dx \leq c$.
From \eqref{existence-energy-estimate-a}, we have
\[
\int_{\Omega}|\mathbf{F}|^2dx \leq \int_{\Omega}\left(|\mathbf{F}|^{p(x)}+1\right)dx
\leq c+|\Omega|,
\]
and
\begin{equation}\label{u-p-2-b}
\int_{\Omega}\big(|\nabla u|^2+|P|^2\big) dx \leq c(1+|\Omega|).
\end{equation}
Given $p_i^+-p_i^- \leq \omega(250r_i)$, it yields
\begin{align*}
\Big({-\hspace{-0.38cm}\int}_{\Omega_i^5}
\big(|\nabla u|^2+|P|^2\big) dx\Big)^{p_i^+-p_i^-}
&= \Big(\frac{1}{|\Omega_i^5|}\Big)^{p_i^+-p_i^-}
\Big(\int_{\Omega_i^5}\big(|\nabla u|^2+|P|^2\big) dx\Big)^{p_i^+-p_i^-}\\
& \leq c \Big(\frac{1}{|B_{250r_i}|}\Big)^{p_i^+-p_i^-}
\Big(\int_{\Omega_i^5}\big(|\nabla u|^2+|P|^2\big) dx\Big)^{p_i^+-p_i^-}\\
& \leq c\Big(\frac{1}{250r_i}\Big)^{n\omega(250r_i)}
\Big(\int_{\Omega_i^5}\big(|\nabla u|^2+|P|^2\big) dx\Big)^{p_i^+-p_i^-}\\
& \leq c\Big(\int_{\Omega_i^5}\Big(|\nabla u|^2+|P|^2\Big)dx\Big)^{p_i^+-p_i^-}.
\end{align*}
On the other hand, using \eqref{R-1}, \eqref{u-p-2-b} as well as
$\frac{1}{250r_i}\geq \frac{1}{R}\geq\frac{c^*}{R_0}\geq |\Omega|+1$, we find
\begin{align*}
\Big(\int_{\Omega_i^5}\Big(|\nabla u|^2+|P|^2\Big)dx\Big)^{p_i^+-p_i^-}
& \leq \Big(\int_{\Omega}\Big(|\nabla u|^2+|P|^2\Big)dx\Big)^{p_i^+-p_i^-} \\
& \leq c(|\Omega|+1)^{p_i^+-p_i^-}\\
& \leq c\big(\frac{1}{250r_i}\big)^{\omega(250r_i)}
\leq c.
\end{align*}
So, we arrive at \eqref{u-p-2-a} due to the log-H\"older continuity of $p(x)$.
Recalling $\gamma_1\leq p_i^-$ and \eqref{u-p-2-a} with $\lambda>1$, we obtain
\begin{align*}
{-\hspace{-0.38cm}\int}_{\Omega_i^5}\Big(|\nabla u|^2+|P|^2\Big)dx
& = \Big({-\hspace{-0.38cm}\int}_{\Omega_i^5}
\big(|\nabla u|^2+|P|^2\Big)dx\Big)^{\frac{p_i^+-p_i^-}{p_i^+}}
\Big({-\hspace{-0.38cm}\int}_{\Omega_i^5}\big(|\nabla u|^2+|P|^2\big)dx
\Big)^{\frac{p_i^-}{p_i^+}} \\
& \leq c^{1/\gamma_1}\Big({-\hspace{-0.38cm}\int}_{\Omega_i^5}
\Big(|\nabla u|^2+|P|^2\Big)dx\Big)^{p^-/p_i^+}\\
& \leq c\Big({-\hspace{-0.38cm}\int}_{\Omega_i^5}
\big(|\nabla u|+|P|\big)^{\frac{2p_i^-}{p^-}}dx\Big)^{p^-/p_i^+}\\
& \leq c\Big({-\hspace{-0.38cm}\int}_{\Omega_i^5}
\big(|\nabla u|+|P|\big)^{2p(x)/p^-}dx+1\Big)^{p^-/p_i^+} \\
&\leq c\lambda^{p^-/p_i^+}.
\end{align*}
Similarly, recalling $\delta\lambda_0\geq 1$ and $\lambda \geq M\lambda_0$ we find
\begin{align*}
{-\hspace{-0.38cm}\int}_{\Omega_i^5}|\mathbf{F}|^2dx
&\leq c\Big({-\hspace{-0.38cm}\int}_{\Omega_i^5}|\mathbf{F}|^{2p(x)/p^-}dx
+1\Big)^{p^-/p_i^+}\\
&\leq c\left(\delta\lambda+1\right)^{p^-/p_i^+}\\
&\leq c\left(\delta\lambda+\delta\lambda_0\right)^{p^-/p_i^+}\\
&\leq c\lambda^{p^-/p_i^+}\delta^{\gamma_1/\gamma_2}.
\end{align*}
We now define
\begin{gather*}
\tilde{\mathbf{A}}_i(x)=\mathbf{A}(25r_i x),\quad
\tilde u_i(x)=\frac{u(25r_i x)}{25r_i\sqrt{c_3\lambda^{p^-/p_i^+}}},\\
\tilde {P}_i(x)=\frac{P(25r_i x)}{\sqrt{c_3\lambda^{p^-/p_i^+}}}, \quad
\tilde{\mathbf{F}}_i(x)=\frac{\mathbf{F}(25r_i x)}{\sqrt{c_3\lambda^{p^-/p_i^+}}}.
\end{gather*}
By Lemma \ref{scaling-normalization}, we obtain that $(\tilde u,\tilde P)$
is a pair of weak solution of
\begin{gather*}
\operatorname{div}({\tilde{\mathbf{A}}}_i(x)\nabla{\tilde u}_i)
-\nabla{\tilde P}_i= \operatorname{div} \tilde{\mathbf{F}}_i,
\quad \text{in } \Omega_5, \\
\operatorname{div} {\tilde u}_i=0, \quad \text{in } \Omega_5,\\
{\tilde u}_i=0, \quad \text{on } \partial\Omega_5.
\end{gather*}
Moreover, using \eqref{Reifenberg-flat-a} and \eqref{CZ-operator-estimate-e}
leads to
\begin{gather*}
\ {-\hspace{-0.38cm}\int}_{\Omega_{6}}|{\tilde{\mathbf{A}}}_i
-({\tilde{\mathbf{A}}}_i)_{\Omega_{6}}|^2dx \leq \delta^2, \\
{-\hspace{-0.38cm}\int}_{\Omega_5}\left(|\nabla {\tilde u}_i|^2
+|{\tilde P}_i|^2\right)dx \leq 1,\\
{-\hspace{-0.38cm}\int}_{\Omega_5}|\tilde{\mathbf{F}}_i|^2dx
\leq \delta^{\gamma_1/\gamma_2}.
\end{gather*}
In accordance with Lemmas \ref{boundary-bound-estimate-a} and \ref{comparison-1},
we obtain
\begin{gather*}
{-\hspace{-0.38cm}\int}_{\Omega_2}\left(|\nabla {\tilde u}_i
-\nabla {\tilde h}_i|^2+|{\tilde P}_i-{\tilde P_{{\tilde h}_i}}|^2\right)dx
\leq \epsilon^2, \\
\|\nabla{\tilde h}_i\|_{L^\infty(\Omega_1)} \leq c_1, \quad
\|\tilde P_{{\tilde h}_i} \|_{L^\infty(\Omega_1)} \leq c_1.
\end{gather*}
Scaling back with
\[
{\tilde h}_i(x)= \frac{{\bar h}_i(250r_ix)}{25r_i\sqrt{c_3\lambda^{p^-/p_i^+}}},
\]
where ${\bar h}_i$ is the weak solution of \eqref{comparison-equ-4},
replacing $B_4^+$ and $T_4$ by $B_i^{4+}$ and $T_i^{4}$ respectively,
and extending $B_i^{4+}$ to $\Omega_i^{4}$, we obtain
\begin{gather}\label{boundary-comparision-a}
{-\hspace{-0.38cm}\int}_{\Omega_i^2}
\left(|\nabla u-\nabla {\bar h}_i|^2+|P- P_{{\bar h}_i}|^2\right)dx
\leq c_3\lambda^{p^-/p_i^+}\epsilon^2, \\
\label{boundary-comparision-b}
\|\nabla{\bar h}_i\|_{L^\infty(\Omega_i^1)} \leq c_1\lambda^{p^-/p_i^+}, \quad
\| P_{{\bar h}_i}\|_{L^\infty(\Omega_i^1)} \leq c_1\lambda^{p^-/p_i^+}
\end{gather}
for $c_1 >1$.
Similar to \eqref{boundary-comparision-a} and \eqref{boundary-comparision-b},
one can deduce the interior estimates
\begin{gather}\label{interior-comparision-a}
{-\hspace{-0.38cm}\int}_{B_i^2}\left(|\nabla u-\nabla w_i|^2
+|P- P_{ w_i}|^2\right)dx \leq c_4\lambda^{p^-/p_i^+}\epsilon^2, \\
\label{interior-comparision-b}
\|\nabla w_i\|_{L^\infty(B_i^1)} \leq c_2\lambda^{p^-/p_i^+}, \quad
\| P_{ w_i}\|_{L^\infty(B_i^1)} \leq c_2\lambda^{p^-/p_i^+}
\end{gather}
for constants $c_2,c_4 >1$, where $w_i$ is the weak solution of
\eqref{interior-comparison-equ-3} when replacing $B_3$ by $B_i^{3}$.
\smallskip
\noindent\textbf{Step 3.}
We want to make an estimate of the super-level
$E\left(A\lambda, \Omega_{R}\right)$.
For any fixed point $x \in \Omega$, we select a universal constant $R$
satisfying $ R \leq \min \{\frac{R_0}{2},\frac{R_0}{|\Omega|+1},1\}$,
and there exists a constant $\delta=\delta(\epsilon)>0$ as given in
Lemmas \ref{comparison-1} and \ref{interior-bound-estimate-a}. Let
\begin{equation}\label{C0-K-Def}
c_0=\max\big\{c_1,c_2\big\}.
\end{equation}
For any $x\in E(A\lambda,\Omega_{R})$, we consider the collection
$\mathcal{B}_\lambda$ of all subsets of $\Omega_{r}(x)$.
By the Vitali covering argument, we extract a countable sub-collection
$\{\Omega_i^0\} \in \mathcal{B}_\lambda$ such that five times larger ball
$\Omega_{5r_i}(x_i)$ covers almost all $E(A\lambda,\Omega_{R})$ and
the balls $\{\Omega_i^0 \}_{i=1}^{\infty}$ are pointwise disjoints with
$\Omega_i^0=\Omega_{r_i}(x_i)$ for $i\in \mathbb{N}$. This leads to
\begin{gather*}
\Omega_i^0 \cap \Omega_j^0=\emptyset, \quad \text{whenever } i\neq j, \\
E(A\lambda,\Omega_{R})\subset \cup_{i\in \mathbb{N}} \Omega_{5r_i}(x_i)
\cup \mathcal{N}_\lambda,
\end{gather*}
with $|\mathcal{N}_\lambda|=0$.
Let $A=(8c_0)^{\gamma_2/\gamma_1}$. We separate the above resulting estimation
into the two cases of interior and boundary, and deduce that
\begin{equation}\label{E-estimate-a}
\begin{aligned}
|E(A\lambda,\Omega_{R})|
&= |E((8c_0)^{\gamma_2/\gamma_1}\lambda,\Omega_{R}) | \\
& = \big| \big\{x\in \Omega_R: (|\nabla u|+|P|)^{2p(x)/p^-}
\geq (8c_0)^{\gamma_2/\gamma_1}\lambda \big\} \big| \\
& \leq \sum_{i \geq 1} \big| \big\{x\in \Omega_{5r_i}:
(|\nabla u|+|P|)^2 \geq 8c_0\lambda^{p^-/p(x)} \big\} \big| \\
& \leq \sum_{i \geq 1} \big| \big\{x\in \Omega_{5r_i}: |\nabla u|^2+|P|^2
\geq 4c_0\lambda^{p^-/p(x)} \big\} \big| \\
& = \sum_{\rm interior case} \big| \big\{x\in \Omega_{5r_i}:
|\nabla u|^2+|P|^2 \geq 4c_0\lambda^{p^-/p(x)} \big\} | \\
&\quad + \sum_{\rm boundary case} | \big\{x\in \Omega_{5r_i}:
|\nabla u|^2+|P|^2 \geq 4c_0\lambda^{p^-/p(x)} \big\} \big|.
\end{aligned}
\end{equation}
For the interior setting, from \eqref{interior-comparision-a}-\eqref{C0-K-Def}
we see that
\begin{equation}\label{E-interior-estimate}
\begin{aligned}
&\big|\big\{x\in \Omega_{5r_i}: |\nabla u|^2+|P|^2
\geq 4c_0\lambda^{p^-/p(x)} \big\}\big| \\
& = \big|\big\{x\in B_i^1: |\nabla u|^2+|P|^2 \geq 4c_2\lambda^{p^-/p(x)} \big\}\big|
\\
& \leq \big|\big\{x\in B_i^1: |\nabla u-\nabla w_i|^2
+|P- P_{ w_i}|^2\geq c_2\lambda^{p^-/p_i^+} \big\}\big| \\
&\quad + \big|\big\{x\in B_i^1: |\nabla w_i|^2+|P_{ w_i}|^2
\geq c_2\lambda^{p^-/p_i^+} \big\}\big| \\
& \leq \frac{1}{ c_2\lambda^{p^-/p_i^+}}\int_{B_i^1}
\big(|\nabla u-\nabla w_i|^2+|P- P_{ w_i} |^2\big)\,dx \\
& \leq c\epsilon^2 |B_i^1 |
\leq c\epsilon^2 |B_i^0|,
\end{aligned}
\end{equation}
where we applied the weak $(1,1)$-type estimate
\[
\big|\big\{x \in E: f(x)> \lambda\big\}\big| \leq \frac{1}{\lambda} \int_E f(x)dx.
\]
Similarly, for the boundary setting one can derive that
\begin{equation}\label{E-boundary-estimate}
\begin{aligned}
&\big|\big\{x\in \Omega_{5r_i}: |\nabla u(x)|^2+|P(x)|^2
\geq 4c_0\lambda^{p^-/p(x)} \big\}\big| \\
& = \big|\big\{z\in \Omega_{5r_i}: |\nabla u(z)|^2+|P(z)|^2
\geq 4c_0\lambda^{\frac{p^-}{p(z)}} \big\}\big| \\
& \leq \big|\big\{z\in \Omega_i^1: |\nabla u(z)|^2+|P(z)|^2
\geq 4c_1\lambda^{\frac{p^-}{p(z)}} \big\}\big| \\
& \leq \big|\big\{z\in \Omega_i^1: |\nabla u-\nabla {\bar h}_i |^2
+|P- P_{{\bar h}_i} |^2 \geq c_1\lambda^{p^-/p_i^+} \big\}\big| \\
&\quad + \big|\big\{z\in \Omega_i^1: |\nabla {\bar h}_i |^2
+ |P_{{\bar h}_i} |^2 \geq c_1\lambda^{p^-/p_i^+} \big\}\big| \\
& \leq \frac{1}{c_1\lambda^{p^-/p_i^+}}
\int_{\Omega_i^1}\big( |\nabla u-\nabla {\bar h}_i |^2
+|P- P_{{\bar h}_i}|^2\big) dz \\
& \leq c\epsilon^2 |\Omega_i^6| \\
& \leq c\epsilon^2 \big(\frac{16}{7}\big)^n|\Omega_i^0|.
\end{aligned}
\end{equation}
This is so because
\begin{gather*}
c\epsilon^2 |\Omega_i^6| \leq c\epsilon^2|B_{150r_i}|=c\epsilon^2|B_{r_i}|,\\
c\epsilon^2 |B_{r_i} | \leq c\epsilon^2 \frac{| B_i^0 |}{|\Omega_i^0|}
|\Omega_i^0 |
\leq c\epsilon^2 \big(\frac{1}{1-\delta}\big)^n |\Omega_i^0|
\leq c\epsilon^2 \big(\frac{16}{7}\big)^n |\Omega_i^0 |.
\end{gather*}
By \eqref{E-estimate-a}-\eqref{E-boundary-estimate}, we have
\begin{equation}\label{vitali-cover-0}
|E(A\lambda,\Omega_R)|
\leq c\epsilon^2 \sum_{i \geq 1} |\Omega_i^0|.
\end{equation}
Using the Vitali covering argument and \eqref{omega-estimate-c}, we obtain
\begin{equation}\label{E-estimate-b}
\begin{aligned}
&|E(A\lambda,\Omega_R)| \\
& \leq c\epsilon^2 \sum_{i \geq 1} \big|\Omega_{r_i}( x_i)\cap
E (\frac{\lambda}{4},\Omega_{2R} ) \big| \\
&\quad + c\epsilon^2\frac{\eta}{(\zeta \lambda)^{\eta}}
\sum_{i \geq 1} \int_{\zeta\lambda}^{\infty}\mu^{\eta}
\big|\big\{x\in \Omega_{r_i}( x_i): |\mathbf{F}|^{\frac{2p(x)}{P^-}}>\mu
\big\}\big|\frac{d\mu}{\mu} \\
& \leq c\epsilon^2 \big|E( \frac{\lambda}{4},\Omega_{2R}) \big|
+c\epsilon^2\frac{\eta}{(\zeta \lambda)^{\eta}}
\int_{\zeta\lambda}^{\infty}\mu^{\eta}
\big|\big\{x\in \Omega_{2R}: |\mathbf{F}|^{\frac{2p(x)}{P^-}}>\mu\big\}\big|
\frac{d\mu}{\mu}.
\end{aligned}
\end{equation}
\smallskip
\noindent\textbf{Step 4.}
We prove that $\|(|\nabla u|+|P|)^{p(x)}\|_{L(t,q)(\Omega_{2R})}<\infty$
in the case $0A\lambda\big\}\big|
\Big)^{q/t}\frac{d\lambda}{A\lambda} \\
&\leq c \epsilon^{\frac{2q}{t}}(I_1+I_2),
\end{aligned}
\end{equation}
where $c$ depends on $n,\gamma_1,\gamma_2,\nu,\Lambda,q,t$ and
$\omega(\cdot)$, and
\begin{equation}\label{E-integral2}
\begin{gathered}
I_1= \frac{tp^-}{2} \int_{0}^{\infty}
\Big(\lambda^{tp^-/2}\big|\big\{x\in \Omega_{2R}:
(|\nabla u|+|P|)^{2p(x)/p^-} >\frac{\lambda}{4}\big\}\big|\Big)^{q/t}
\frac{d\lambda}{\lambda} \\
I_2= \frac{tp^-}{2} \int_{0}^{\infty} \lambda^{q(\frac{p^-}{2}
-\frac{\eta}{t})}\Big(\int_{\zeta\lambda}^{\infty} \mu^{\eta}
\big|\big\{x\in \Omega_{2R}:
|\mathbf{F}|^{2p(x)/p^-}>\mu\big\}\big|\frac{d\mu}{\mu}\Big)^{q/t}
\frac{d\lambda}{\lambda}
\end{gathered}
\end{equation}
Thanks to \eqref{result-3}, we have
\begin{equation}\label{result-4}
\begin{aligned}
&\|(|\nabla u|+|P|)^{p(x)}\|^q_{L(t,q)(\Omega_{R})} \\
&= \|(|\nabla u|+|P|)^{2p(x)/p^-}\|^{\frac{p^-}{2}q}_{L(\frac{tp^-}{2},
\frac{qp^-}{2})(\Omega_{R})} \\
&= \frac{tp^-}{2} \int_0^\infty \Big(\mu^{tp^-/2}\big|x \in \Omega_{R}:
(|\nabla u|+|P|)^{2p(x)/p^-}>\mu\big|\Big)^{q/t} \frac{d\mu}{\mu}.
\end{aligned}
\end{equation}
Making a change of variables yields
\[
I_1 =c(q)\big\|(|\nabla u|+|P|)^{p(x)}\big\|^q_{L(t,q)(\Omega_{2R})}.
\]
For estimating $I_2$, we consider two cases.
\smallskip
\noindent\textbf{Case 1.}
If $q \geq t$, note that \eqref{Hardy's inequality-condition} is satisfied
because of $|\mathbf{F}|^{2p(x)/p^-}\in L^{\eta}(\Omega_{2R})$.
By making the change of variables $\bar \lambda= \zeta\lambda$ and
$\zeta=\frac {\delta}4$, in view of $\alpha=q/t \geq 1$ and
$ r=q\big(\frac{p^-}{2}-\frac{\eta}{t}\big)>0$,
and using
\[
f(\mu) =\mu^{\eta-1}\big|\big\{x\in \Omega_{2R}:
|\mathbf{F}|^{\frac{2p(x)}{P^-}}>\mu\big\}\big|,
\]
it follows from Lemma \ref{Hardy's inequality} that
\begin{align*}
I_2
&= c\frac{tp^-}{2} \int_{0}^{\infty}\bar{\lambda}^{q(\frac{p^-}{2}
-\frac{\eta}{t})}\Big(\int_{\bar\lambda}^{\infty}
\mu^{\eta}
\big|\big\{x\in \Omega_{2R}: |\mathbf{F}|^{2p(x)/p^-}>\mu \big\}\big|
\frac{d\mu}{\mu}\Big)^{q/t}\frac{d\bar\lambda}{\bar\lambda}\\
&\leq c\frac{tp^-}{2} \int_{0}^{\infty} \bar\lambda^{\frac{qp^-}{2}}
\big|\big\{x\in \Omega_{2R}: |\mathbf{F}|^{2p(x)/p^-}>\bar\lambda
\big\}\big|^{q/t}\frac{d\bar\lambda}{\bar\lambda}\\
& = c\big\||\mathbf{F}|^{p(x)}\big\|_{L(t,q)(\Omega_{2R})}^q,
\end{align*}
where $c=c(\gamma_1,\gamma_2,q,t)$.
\smallskip
\noindent\textbf{Case 2.}
If $0\mu
\big\}\big|^{q/t},\quad
r=\frac{\eta q}{t},\ \alpha_1=1<\frac{t}{q}=\alpha_2, \quad
\varepsilon=1,
$$
we have
\begin{align*}
& \Big(\int_{\lambda}^{\infty}\mu^{\eta}
\big|\big\{x\in \Omega_{2R}: |\mathbf{F}|^{2p(x)/p^-}>\mu\big\}\big|
\frac{d\mu}{\mu}\Big)^{q/t}\\
&\leq \lambda^{\frac{\eta q}{t}}
\big|\big\{x\in \Omega_{2R}: |\mathbf{F}|^{2p(x)/p^-}>\lambda
\big\}\big|^{q/t} \\
&\quad +c \int_{\lambda}^{\infty}\mu^{\frac{\eta q}{t}}
\big|\big\{x\in \Omega_{2R}: |\mathbf{F}|^{2p(x)/p^-}>\mu \big\}\big|^{q/t}
\frac{d\mu}{\mu}.
\end{align*}
Hence, after changing the variable $\zeta \lambda \to \lambda$, it follows
by Fubini's theorem that
\begin{align*}
I_2
&\leq c\frac{tp^-}{2} \int_{0}^{\infty}\lambda^{q(\frac{p^-}{2}-\frac{\eta}{t})}
(\lambda)^{\frac{\eta q}{t}}\big|\big\{x\in \Omega_{2R}:
|\mathbf{F}|^{2p(x)/p^-}>\lambda\big\}\big|^{q/t}\frac{d\lambda}{\lambda}\\
&\quad + c\frac{tp^-}{2} \int_{0}^{\infty}\lambda^{q(\frac{p^-}{2}
-\frac{\eta}{t})}\int_{\lambda}^{\infty}\mu^{{\frac{\eta q}{t}}-1}
\big|\big\{x\in \Omega_{2R}: |\mathbf{F}|^{2p(x)/p^-}>\mu\big\}\big|^{q/t} d\mu
\frac{d\lambda}{\lambda}\\
&\leq c\||\mathbf{F}|^{p(x)}\|_{L(t,q)(\Omega_{2R})}^q \\
&\quad +c\frac{tp^-}{2}\int_{0}^{\infty}\lambda^{q(\frac{p^-}{2}-\frac{\eta}{t})}
\Big(\int_{\lambda}^{\infty}\mu^{{\frac{\eta q}{t}}-1}
\big|\big\{x\in \Omega_{2R}: |\mathbf{F}|^{2p(x)/p^-}>\mu\big\}\big|^{q/t} d\mu
\Big)\frac{d\lambda}{\lambda}\\
& \leq c\||\mathbf{F}|^{p(x)}\|_{L(t,q)(\Omega_{2R})}^q,
\end{align*}
where $c=c(\gamma_1,\gamma_2,q,t)$.
Substituting the estimates of $I_1$ and $I_2$ into \eqref{E-integral},
for $t\geq 1$ we derive that
% \label{gradient-u-estimate-a}
\begin{align*}
&\|(|\nabla u|+|P|)^{p(x)}\|_{L(t,q)(\Omega_{R})} \\
&\leq c\Big[\frac{tp^-}{2}\int_{M\lambda_0}^{\infty}
\Big((A\lambda)^{tp^-/2}\big|\big\{x\in \Omega_R:
(|\nabla u|+|P|)^{2p(x)/p^-}>A\lambda\big\}\big|\Big)^{q/t}
\frac{d(A\lambda)}{A\lambda}\Big]^{1/q} \\
&\quad + c\Big[\frac{tp^-}{2}\int_{0}^{M\lambda_0}
\Big((A\lambda)^{tp^-/2} \big|\big\{x\in \Omega_R:
(|\nabla u|+|P|)^{2p(x)/p^-}>A\lambda \big\}\big|\Big)^{q/t}
\frac{d(A\lambda)}{A\lambda}\Big]^{1/q} \\
&\leq c\Big[\frac{tp^-}{2}\int_{0}^{M\lambda_0}
\Big((A\lambda)^{tp^-/2}\big|\big\{x\in \Omega_R:
(|\nabla u|+|P|)^{2p(x)/p^-}>A\lambda \big\}\big|\Big)^{q/t}
\frac{d(A\lambda)}{A\lambda}\Big]^{1/q} \\
&\quad + \bar c \epsilon^{2/t}
\Big(\|(|\nabla u|+|P|)^{p(x)} \|_{L(t,q)(\Omega_{2R})}+
\| |\mathbf{F}|^{p(x)}\|_{L(t,q)(\Omega_{2R})}\Big) \\
&\leq c\lambda_0^{p^-/2}|\Omega_{2R}|^{1/t}
+ \bar c \epsilon^{2/t}\Big(\big\|(|\nabla u|+|P|)^{p(x)}
\big\|_{L(t,q)(\Omega_{2R})}
+ \big\||\mathbf{F}|^{p(x)}\big\|_{L(t,q)(\Omega_{2R})}\Big),
\end{align*}
where $\bar c= \bar c(n,\gamma_1,\gamma_2,\nu,\Lambda,q,t,\omega(\cdot))$.
It suffices to choose $\epsilon>0$ small enough such that
$\bar c \epsilon^{2/t} \leq \frac{1}{2}$. Once the selection
of $\epsilon$ is fixed, we can find the corresponding constant
$\delta=\delta(\epsilon,\gamma_1,\gamma_2)$ such that
\begin{equation} \label{gradient-u-estimate-b}
\begin{aligned}
& \big\|(|\nabla u|+|P|)^{p(x)}\big\|_{L(t,q)(\Omega_{R})} \\
&\leq c\lambda_0^{p^-/2}|\Omega_{2R}|^{1/t}
+ \frac{1}{2} \big\|(|\nabla u|+|P|)^{p(x)}\big\|_{L(t,q)(\Omega_{2R})}
+c\big\||\mathbf{F}|^{p(x)}\big\|_{L(t,q)(\Omega_{2R})}.
\end{aligned}
\end{equation}
By a standard iteration argument, we can attain an estimate as
\eqref{Lorentz-estimate-a} in the case of $t>1$ and $0\lambda\}$ in line
with \eqref{E-estimate-b}, we have
\[
E_k(A\lambda,\Omega_R) \leq c \epsilon^2
\big|E_k\big(\frac{\lambda}{4},\Omega_{2R} \big)\big|
+c\epsilon^2\frac{\eta}{(\zeta \lambda)^{\eta}}
\int_{\zeta\lambda}^{\infty}\mu^{\eta}\big|\big\{x\in \Omega_{2R}:
|\mathbf{F}|^{2p(x)/p^-}>\mu \big\}\big|\frac{d\mu}{\mu}
\]
for $k\in \mathbb{N}\cap [M\lambda_0,\infty)$.
Indeed, for $k\leq A\lambda$ we have $E_k(A\lambda,\Omega_R)=\emptyset$, which
implies that the above estimate holds trivially. For $k>A\lambda$, it is
also valid because
$$
E_k(A\lambda,\Omega_R)=E(A\lambda,\Omega_R)
=\big\{x\in \Omega_R,\,(|\nabla u|+|P|)^{p(x)}>A\lambda \big\}
$$
and
$$
E_k\big(\frac{\lambda}{4},\Omega_{2R}\big)
=E\big(\frac{\lambda}{4},\Omega_{2R}\big).
$$
Proceeding in the same manner as the above, one can see that
\eqref{gradient-u-estimate-b} holds with
$\left|(|\nabla u|+|P|)^{p(x)}\right|_k$ in place of $(|\nabla u|+|P|)^{p(x)}$.
Let
\begin{gather*}
B_0=0,\quad
L=c\lambda_0^{p^-/2}|\Omega_{2R}|^{1/t}
+c\big\||\mathbf{F}|^{p(x)}\big\|_{L(t,q)(\Omega_{2R})}, \\
\varphi(\rho)=\big\|\, |(|\nabla u|+|P|)^{p(x)}|_k\big\|_{L(t,q)(\Omega_\rho)}.
\end{gather*}
because
\[
\big\|\,|(|\nabla u|+|P|)^{p(x)}|_k\big\|_{L(t,q)(\Omega_R)}<\infty,
\]
using the iteration argument we have
\[
\big\|\,|(|\nabla u|+|P|)^{p(x)}|_k\big\|_{L(t,q)(\Omega_R)}
\leq c\lambda_0^{p^-/2}|\Omega_R|^{1/t}+
c\big\|\,|\mathbf{F}|^{p(x)}\big\|_{L(t,q)(\Omega_{2R})}.
\]
In what follows, we use a standard finite covering argument to achieve the
global estimate.
Note that $\Omega$ is a bounded domain in $\mathbb{R}^n$.
There exist $N \in\mathbb{N}$ and $x_l \in \Omega$ for $l=1,2,\dots,N $ such that
\[
\overline{\Omega} \subset \cup_{l=1}^N B_R(x_l).
\]
Then we have
\begin{align*}
\big\|\big|(|\nabla u|+|P|)^{p(x)}\big|_k\big\|_{L(t,q)(\Omega)}
& \leq \sum_{l=1}^N\big\|\,|(|\nabla u|+|P|)^{p(x)}|_k\big\|_{L(t,q)(\Omega_R)}\\
& \leq c\sum_{l=1}^N\Big(\lambda_0^{p^-/2}|\Omega_R|^{1/t}
+ \big\||\mathbf{F}|^{p(x)}\big\|_{L(t,q)(\Omega_{2R})}\Big)\\
& \leq c N \Big(\lambda_0^{p^-/2}|\Omega_R|^{1/t}
+\big\||\mathbf{F}|^{p(x)}\big\|_{L(t,q)(\Omega_{2R})}\Big).
\end{align*}
Recalling the definition of $\lambda_0$, we obtain
\begin{equation}\label{gradient-u-estimate-c}
\begin{aligned}
& \big\|\,|(|\nabla u|+|P|)^{p(x)}|_k\big\|_{L(t,q)(\Omega)} \\
&\leq cN |\Omega_R|^{1/t}\Big({-\hspace{-0.38cm}\int}_{\Omega_{2R}}
(|\nabla u|+|P|)^{2p(x)/p^-}dx \\
&\quad +\Big({-\hspace{-0.38cm}\int}_{\Omega_{2R}}
\Big(|\mathbf{F}|^{2p(x)/p^-}+1\Big)^\eta\,dx\Big)^{1/\eta}\Big)^{p^-/2}+
cN \big\|\mathbf{F}|^{p(x)}\big\|_{L(t,q)(\Omega_{2R})} \\
& \leq cN |\Omega_R|^{1/t}\Big({-\hspace{-0.38cm}\int}_{\Omega_{2R}}
(|\nabla u|+|P|)^{2p(x)/p^-}dx \\
&\quad +\Big({-\hspace{-0.38cm}\int}_{\Omega_{2R}}
\Big(|\mathbf{F}|^{2p(x)/p^-}+1\Big)^\eta\,dx\Big)^{1/\eta}\Big)^{p^-/2}+
cN \big\|\mathbf{F}|^{p(x)}\big\|_{L(t,q)(\Omega)}.
\end{aligned}
\end{equation}
Note that
\[
\frac{2p^+}{p^-}
=2\Big(1+\frac{p^+-p^-}{p^-}\Big)
\leq 2\Big(1+\frac{\omega(2R)}{\gamma_1}\Big) \leq 2(1+\sigma),
\]
where $\sigma$ is the same as given in Lemma \ref{reverse-Holder-b}.
Then, it follows from the reverse H\"{o}der's inequality of
Lemma \ref{reverse-Holder-b} that
\begin{equation}\label{gradient-u-estimate-d}
\begin{aligned}
{-\hspace{-0.38cm}\int}_{\Omega_{2R}}(|\nabla u|+|P|)^{2p(x)/p^-}dx
&\leq {-\hspace{-0.38cm}\int}_{\Omega_{2R}}(|\nabla u|+|P|)^{\frac{2p^+}{p^-}}dx+1 \\
& \leq c \Big({-\hspace{-0.38cm}\int}_{\Omega_{4R}}\big(|\nabla u|^2+|P|^2\big) dx
+1 \Big)^{p^+/p^-} \\
&\quad +{-\hspace{-0.38cm}\int}_{\Omega_{4R}}|\mathbf{F}|^{\frac{2p^+}{p^-}}dx.
\end{aligned}
\end{equation}
Using \eqref{existence-energy-estimate-a} and H\"{o}der's inequality we have
\begin{equation}\label{gradient-u-estimate-e}
\begin{aligned}
&\Big({-\hspace{-0.38cm}\int}_{\Omega_{4R}}\big(|\nabla u|^2+|P|^2\big) dx
\Big)^{p^+/p^-} \\
&\leq \Big(\frac{1}{|\Omega_{4R}|}\Big)^{p^+/p^-}
\Big(\int_{\Omega}\big(|\nabla u|^2+|P|^2\big) dx \Big)^{p^+/p^-} \\
& \leq c \Big(\frac{1}{|\Omega_{4R}|}\Big)^{p^+/p^-}
\Big(\int_{\Omega}|\mathbf{F}|^2dx \Big)^{p^+/p^-} \\
& \leq c \Big(\frac{1}{|\Omega_{4R}|}\Big)^{p^+/p^-}|\Omega|^{1-\frac{p^-}{p^+}}
\int_{\Omega}|\mathbf{F}|^{\frac{2p^+}{p^-}}dx \\
& \leq c \Big(\frac{1}{|\Omega_{4R}|}\Big)^{p^+/p^-}\int_{\Omega}
\Big(|\mathbf{F}|^{\frac{2p(x)}{p^-}\frac{p^+}{p^-}}+1 \Big)dx.
\end{aligned}
\end{equation}
From \eqref{gradient-u-estimate-d}-\eqref{gradient-u-estimate-c}, we deduce that
\begin{equation}\label{gradient-u-estimate-c2}
\begin{aligned}
& \|\left|(|\nabla u|+|P|)^{p(x)}\right|_k\|_{L(t,q)(\Omega)} \\
&\leq cN |\Omega_R|^{1/t}\Big\{ \Big(\frac{1}{|\Omega_{4R}|}\Big)^{p^+/p^-}
\int_{\Omega}\Big(|\mathbf{F}|^{\frac{2p(x)}{p^-}\frac{p^+}{p^-}}+1\Big)dx \\
&\quad + \Big[{-\hspace{-0.38cm}\int}_{\Omega_{2R}}\big(|\mathbf{F}|^{2p(x)/p^-}
+1\big)^\eta\,dx\Big]^{1/\eta}\Big\}^{p^-/2}
+ cN \big\||\mathbf{F}|^{p(x)}\big\|_{L(t,q)(\Omega)} \\
&\leq cN |\Omega_R|^{1/t}\Big\{ \Big(\frac{1}{|\Omega_{4R}|}\Big)^{p^+/p^-}
\int_{\Omega}\big(|\mathbf{F}|^{\frac{2p(x)}{p^-}\frac{p^+}{p^-}}+1\big)dx\\
&\quad + \big(\frac{1}{|\Omega_{2R}|} \big)^{1/\eta} \Big[\int_{\Omega}
\big(|\mathbf{F}|^{2p(x)/p^-}+1\big)^\eta\,dx \Big]^{1/\eta}
\Big\}^{p^-/2}
+ cN \big\|\mathbf{F}|^{p(x)}\big\|_{L(t,q)(\Omega)}.
\end{aligned}
\end{equation}
Using a standard Hardy's inequality in the Marcinkiewicz spaces
\cite[Lemma 2.3]{Min07}) and
the reverse H\"older's inequality of Lemma \ref{reverse-Holder-c}, we obtain
\begin{equation} \label{F-estimate-a}
\begin{aligned}
\int_{\Omega}|\mathbf{F}|^{\frac{2p(x)}{p^-}\frac{p^+}{p^-}}\,dx
&\leq \frac{t(p^-)^2}{t(p^-)^2-2p^+}|\Omega|^{1-\frac{2p^+}{t(p^-)^2}}
\| \,|\mathbf{F}|^{2p(x)/p^-}\|_{{\mathcal{M}}^{tp^-/2}(\Omega)}^{p^+/p^-} \\
&= \frac{t(p^-)^2}{t(p^-)^2-2p^+}|\Omega|^{1-\frac{2p^+}{t(p^-)^2}} \\
&\quad\times \Big\{\sup_{h>0} \Big[h^{tp^-/2} \big|\big\{x\in \Omega:
|\mathbf{F}|^{2p(x)/p^-}>h \big\}
\big|\Big]^{\frac{2}{tp^-}}\Big\}^{p^+/p^-} \\
& \leq c|\Omega|^{1-\frac{2p^+}{t(p^-)^2}}
\| \,|\mathbf{F}|^{2p(x)/p^-}\|^{p^+/p^-}_{L(\frac{tp^-}{2},\frac{qp^-}{2})(\Omega)} \\
& \leq c|\Omega|^{1-\frac{2p^+}{t(p^-)^2}} \||\mathbf{F}|^{p(x)}
\|^{\frac{2p^+}{(p^-)^2}}_{L(t,q)(\Omega)}.
\end{aligned}
\end{equation}
Similarly, we can derive that
\begin{align*}
\Big(\int_{\Omega}|\mathbf{F}|^{\frac{2p(x)}{p^-}\eta}\,dx\Big)^{1/\eta}
&\leq c(\gamma_1,\gamma_2,q,t)|\Omega|^{\frac{1}{\eta}-\frac{2}{tp^-}}
\| \,|\mathbf{F}|^{2p(x)/p^-}\|_{L(\frac{tp^-}{2},\frac{qp^-}{2})(\Omega)}\\
&\leq c(\gamma_1,\gamma_2,q,t)|\Omega|^{\frac{1}{\eta}-\frac{2}{tp^-}}
\big\|\,|\mathbf{F}|^{p(x)}\big\|^{\frac{2}{p^-}}_{L(t,q)(\Omega)}.
\end{align*}
In the case $q<\infty$, it follows from \eqref{gradient-u-estimate-c} that
\begin{align*}
\big\| \big|\big(|\nabla u|+|P|\big)^{p(x)}\big|_k\big\|_{L(t,q)(\Omega)}
&\leq cN \Big\{\Big(\frac{|\Omega|}{|\Omega_{2R}|}\Big)^{\frac{p^+}{2}
-\frac{1}{t}}\Big(\||\mathbf{F}|^{p(x)}
\|^{p^+/p^-}_{L(t,q)(\Omega)}+1\Big) \\
&\quad + \Big(\frac{|\Omega|}{|\Omega_{2R}|}\Big)^{\frac{p^-}{2\eta}-\frac{1}{t}}
\big[\||\mathbf{F}|^{p(x)}\|_{L(t,q)(\Omega)}+1\big]\Big\}\\
& \leq cN \Big\{\Big(\frac{|\Omega|}{|B_{2R}|}\frac{|B_{2R}|}{|\Omega_{2R}|}
\Big)^{\frac{p^+}{2}-\frac{1}{t}}
\Big(\| |\mathbf{F}|^{p(x)}\|^{p^+/p^-}_{L(t,q)(\Omega)}+1\Big)\\
&\quad + \Big(\frac{|\Omega|}{|B_{2R}|}\frac{|B_{2R}|}{|\Omega_{2R}|}
\Big)^{\frac{p^-}{2\eta}-\frac{1}{t}}
\big[\| \,|\mathbf{F}|^{p(x)}\|_{L(t,q)(\Omega)}+1\big]\Big\}\\
& \leq cN \Big\{\big[\frac{|\Omega|}{|B_{2R}|}\big(\frac{2}{1-\delta}
\big)^n\big]^{\frac{p^+}{2}
-\frac{1}{t}}\Big(\| |\mathbf{F}|^{p(x)}\|^{p^+/p^-}_{L(t,q)(\Omega)}+1\Big)\\
&\quad + \big[\frac{|\Omega|}{|B_{2R}|}\big(\frac{2}{1-\delta}\big)^n
\big]^{\frac{p^-}{2\eta}-\frac{1}{t}}
\Big(\| |\mathbf{F}|^{p(x)}\|_{L(t,q)(\Omega)}+1\Big)\Big\}\\
& \leq cN \big(\| |\mathbf{F}|^{p(x)}\|_{L(t,q)(\Omega)}+1\big)^{p^+/p^-}\\
& \leq cN \big(\| |\mathbf{F}|^{p(x)}\|_{L(t,q)(\Omega)}+1
\big)^{\gamma_2/\gamma_1}.
\end{align*}
Taking $k\to \infty$, by the lower semi-continuity of Lorentz quasi-norm we have
\[
\big\| \,|(|\nabla u|+|P|)^{p(x)}|_k\big\|_{L(t,q)(\Omega)}
\leq cN \Big(\big\||\mathbf{F}|^{p(x)}\big\|_{L(t,q)(\Omega)}
+1\Big)^{^{\gamma_2/\gamma_1}},
\]
where $c$ depends only on $n,\gamma_1,\gamma_2,\nu,\,\Lambda,t,q,\omega(\cdot),R_0$
and $|\Omega|$.
\smallskip
\noindent\textbf{Step 6.}
Finally, for the case $q=\infty$, we obtain back to the second inequality
in \eqref{alternatives} and split it into two parts with a small $\iota>0$
to be determined later:
\[
\big(\frac{\lambda}{2}\big)^{\eta}
\leq \frac{1}{\delta^\eta}{-\hspace{-0.38cm}\int}_{\Omega_i^0}
|\mathbf{F}|^{\frac{2p(x)}{p^-}\eta}dx
\leq \frac{(\iota\lambda)^{\eta}}{\delta^\eta}+\frac{1}{\delta^\eta
|\Omega_i^0|}\int_{\{x\in \Omega_i^0: |\mathbf{F}|^{2p(x)/p^-}>\iota\lambda
\}}|\mathbf{F}|^{\frac{2p(x)}{p^-}\eta}dx.
\]
Set
\[
G(\iota\lambda,\Omega_i^0)=\{x\in \Omega_i^0: |\mathbf{F}|^{2p(x)/p^-}
>\iota\lambda \}, \quad
G(\mu,\Omega_i^0)=\{x\in \Omega_i^0: |\mathbf{F}|^{2p(x)/p^-}>\mu \}.
\]
Similar to \eqref{F-estimate-a}, by using H\"older's inequality we obtain
\begin{align*}
&\big(\frac{\lambda}{2}\big)^{\eta}- \big(\frac{\iota\lambda}{\delta}\big)^{\eta}\\
& \leq \frac{1}{\delta^\eta |\Omega_i^0|}
\int_{\{x\in \Omega_i^0: |\mathbf{F}|^{2p(x)/p^-}>\iota\lambda \}}
|\mathbf{F}|^{\frac{2p(x)}{p^-}\eta}dx\\
& \leq \frac{t}{(t-\eta)\delta^\eta }
\frac{|G(\iota\lambda,\Omega_i^0)|^{1-\frac{\eta}{t}}}{ |\Omega_i^0|}
\sup_{\mu>0} \mu^{\eta} \big|\big\{x \in G(\iota\lambda,\Omega_i^0 ):
|\mathbf{F}|^{\frac{2p(x)}{p^-}\eta}\geq \mu \big\}\big|^{\eta/t}\\
& \leq \frac{t |G(\iota\lambda,\Omega_i^0) |^{1-\frac{\eta}{t}}}
{(t-\eta)\delta^\eta |\Omega_i^0|}
\Big[(\iota\lambda)^{\eta} | G(\iota\lambda,\Omega_i^0)
|^{\eta/t}+\sup_{\mu>\iota\lambda}
\mu^{\eta} |G(\mu,\Omega_i^0) |^{\eta/t}\Big]\\
& = \frac{t}{(t-\eta)\delta^\eta}\Big[(\iota\lambda)^{\eta}+
\frac{|G(\iota\lambda,\Omega_i^0)|^{1-\frac{\eta}{t}}}{ |\Omega_i^0|}
\sup_{\mu>\iota\lambda} \mu^{\eta}
|G(\mu,\Omega_i^0)|^{\eta/t}\Big].
\end{align*}
Choose $\iota>0$ appropriately small to satisfy
\[
(\frac{\lambda}{2})^{\eta}-(\frac{\iota \lambda}{\delta})^\eta
- \frac{t}{t-\eta}(\frac{\iota \lambda}{\delta})^\eta
= (\frac{\lambda}{2})^{\eta} - (\frac{\iota \lambda}{\delta})^\eta
(1+\frac{t}{t-\eta})
\geq (\frac{\lambda}{4})^{\eta}.
\]
Then there exists a positive constant $c(t)$ depending only on $t$ such that
$\iota\le c(t)\delta$.
Thus, we obtain
\begin{equation}\label{second-infty}
\begin{aligned}
|\Omega_i^0|
&\leq \frac{c t}{t-\eta}
\frac{|G(\iota\lambda^i,\Omega_i^0)|^{1-\frac{\eta}{t}}}
{(\iota\lambda)^{\eta}} \Big(\sup_ {\mu>\iota\lambda}
\mu^t |G(\mu,\Omega_i^0) |\Big)^{\eta/t} \\
& \leq \frac{c t(\iota\lambda)^{-t}}{t-\eta}
\Big((\iota\lambda)^{t} |G(\iota\lambda,\Omega_i^0) |\Big) ^{1-\frac{\eta}{t}}
\Big(\sup_ {\mu>\iota\lambda} \mu^{t} |G(\mu,\Omega_i^0) |\Big)^{\eta/t} \\
& \leq \frac{ct(\iota\lambda)^{-t}}{t-\eta}
\sup_{\mu>\iota\lambda} \mu^{t} |G(\mu,\Omega_i^0)|.
\end{aligned}
\end{equation}
Substituting \eqref{alternatives} into \eqref{vitali-cover-0}
(namely, plugging \eqref{omega-estimate-a} and \eqref{second-infty}
into \eqref{vitali-cover-0}) yields
\begin{equation}\label{E-estimate-c}
|E(A\lambda,\Omega_R)|\leq c\epsilon
\big|E \big(\frac{\lambda}{4},\Omega_{2R} \big) \big|
+c\epsilon (\iota\lambda)^{-t}\sup_{\mu>\iota\lambda} \mu^{t}|G(\mu,\Omega_{2R})|.
\end{equation}
Multiplying both sides of \eqref{E-estimate-c} by $(A\lambda)^{tp^-/2}$
and taking the supremum with respect to $\lambda$ over $(M\lambda_0,\infty)$,
we deduce that
\begin{align*}
&\sup_{\lambda >M\lambda_0} (A\lambda)^{tp^-/2}
\big|\big\{x\in \Omega_R:(|\nabla u|+|P|)^{2p(x)/p^-}>A\lambda\big\}\big|\\
& \leq c \epsilon^2 A^{tp^-/2}
\Big(\sup_{\lambda>M\lambda_0} \lambda^{tp^-/2}
\big|\big\{x\in \Omega_{2R}:(|\nabla u|+|P|)^{2p(x)/p^-}>
\frac{\lambda}{4}\big\}\big| \\
&\quad +\sup_ {\lambda >M\iota\lambda_0} \lambda^{\frac{tp^-}{2}-t}
\Big(\sup_{\mu>\lambda} \mu^{t}|G(\mu,\Omega_{2R})|\Big)\Big)\\
& \leq c \epsilon^2 \Big(\sup_{\lambda>M\lambda_0}
\lambda^{tp^-/2}\big|\big\{x\in \Omega_{2R}:
(|\nabla u|+|P|)^{2p(x)/p^-}>\frac{\lambda}{4}\big\}\big| \\
&\quad +\sup_ {\lambda >M\iota\lambda_0}
\Big(\sup_{\mu>\lambda} \mu^{tp^-/2}|G(\mu,\Omega_{2R})|\Big)\Big).
\end{align*}
Note that
$$
\sup_{\lambda>M\iota\lambda_0}\sup_{\mu>\lambda}\mu^{tp^-/2}
|G(\mu,\Omega_{2R})|
\leq \big\| |\mathbf{F}|^{p(x)}\big\|_{\mathcal{M}^{t}(\Omega_{2R})}^t.
$$
Let $\epsilon>0$ be so small that $c \epsilon^{2/t}\leq \frac{1}{2 }$.
Then
\begin{align*}
& \big\|(|\nabla u|+|P|)^{p(x)}\big\|_{\mathcal{M}^{t}(\Omega_{R})}\\
& \leq c \epsilon^{2/t}\Big(\big\|(|\nabla u|+|P|)^{p(x)}
\big\|_{\mathcal{M}^{t}(\Omega_{2R})}
+ c(\gamma_1,\gamma_2,q,t)\big\||\mathbf{F}|^{p(x)}
\big\|_{\mathcal{M}^{t}(\Omega_{2R})}\Big) \\
&\quad +c |\Omega_{2R}|^{1/t}M\lambda_0^{p^-/2}\\
& \leq \frac{1}{2}\big\|(|\nabla u|+|P|)^{p(x)}
\big\|_{\mathcal{M}^{t}(\Omega_{2R})}
+ c \big\||\mathbf{F}|^{p(x)}\big\|_{\mathcal{M}^{t}(\Omega_{2R})}\\
&\quad +c |\Omega_{2R}|^{1/t} \Big[{-\hspace{-0.38cm}\int}_{\Omega_{2R}}
(|\nabla u|+|P|)^{2p(x)/p^-}dx \\
&\quad +\Big({-\hspace{-0.38cm}\int}_{\Omega_{2R}}
\Big(|\mathbf{F}|^{2p(x)/p^-}+1\Big)^\eta\,dx\Big)^{1/\eta}\Big]^{p^-/2}.
\end{align*}
The rest of the proof is closely similar to the argument of Step 5.
Consequently, we arrive at the desired result for the case $q=\infty$.
\end{proof}
\subsection*{Acknowledgements}
This work is supported by NSFC grant 11371050 and NSFC-ERC grant 11611530539,
and partially supported by NSF grant DMS-1204497.
\begin{thebibliography}{00}
\bibitem{AcM05} Acerbi, E.; Mingione, G.;
\emph{Gradient estimates for the $p(x)$-Laplacean system},
J. Reine Angew. Math., \textbf{584} (2005), 117--148.
\bibitem{AcM07} Acerbi, E.; Mingione, G.;
\emph{Gradient estimates for a class of parabolic systems},
Duke Math. J., \textbf{136} (2) (2007), 285--320.
\bibitem{AdP15} Adimurthi, K.; Phuc, N. C.;
\emph{Global Lorentz and Lorentz-Morrey estimates below the natural exponent
for quasilinear equations}, Calc. Var. Partial Differential Equations,
\textbf{54} (2015), 3107--3139.
\bibitem{AdMZ17} Adjal, H.; Moussaoui, M.; Zine, A.;
\emph{ Exitence and regularity of solutions for Stokes systems with non-smooth
boundary data in a polyhedron}, Electron J. Differ. Equ., \textbf{2017 } (147) (2017), 1--10.
\bibitem{Bar13} Baroni, P.;
\emph{Lorentz estimates for degenerate and singular evolutionary systems},
J. Differential Equations, \textbf{255} (2013), 2927--2951.
\bibitem{Bar14} Baroni, P.;
\emph{Lorentz estimates for obstacle parabolic problems},
Nonlinear Anal., \textbf{96} (2014), 167--188.
\bibitem{Bre13} Breit, D.;
\emph{Smoothness properties of solutions to the nonlinear Stokes problem
with nonautonomous potentials}, Comment. Math. Univ. Carolin.,
\textbf{54} (4) (2013), 493--508.
\bibitem{BrF11} Breit, D.; Fuchs, M.;
\emph{The nonlinear Stokes problem with general potentials having superquadratic
growth}, J. Math. Fluid Mech., \textbf{13} (3) (2011), 371--385.
\bibitem{BuBS16} Bul\'{i}\v{c}ek, M.; Burczak, J.; Schwaarzacher, S.;
\emph{A unified theory for some non-Newtonian fluids under singular forcing},
SIAM J. Math. Pura. Anal., \textbf{48} (6) (2016), 4241--4264.
\bibitem{ByOW14} Byun, S. S.; Ok, J.; Wang, L. H.;
\emph{$W^{1,p(x)}$-Regularity for elliptic equations with measurable
coefficients in nonsmooth domains}, Commun. Math. Phys., \textbf{329} (2014),
937--958.
\bibitem{ByW04} Byun, S. S.; Wang, L. H.;
\emph{Elliptic equations with BMO coefficients in Reifenberg domains},
Comm. Pure Appl. Math., \textbf{57} (10) (2004), 1283--1310.
\bibitem{ByS17} Byun, S. S.; So, H.;
\emph{Weighted estimates for generalized steady Stokes systems in nonsmooth domains},
J. Math. Phys., \textbf{58} (2) (2017), 023101--19.
\bibitem{CaP98} Caffarelli, L. A.; Peral, I.;
\emph{On $W^{1,p}$ estimates for elliptic equations in divergence form},
Comm. Pure Appl. Math., \textbf{51} (1998), 1--21.
\bibitem{ChFL91} Chiarenza, F.; Frasca, M.; Longo, P.;
\emph{Interior $W^{2,p}$ estimates for nondivergence elliptic equations with
discontinuous coeffcients}, Ricerche Mat., \textbf{40} (1991), 149-168.
\bibitem{ChFL93} Chiarenza, F.; Frasca, M.; Longo, P.;
\emph{$W^{2,p}$ solvability of the Dirichlet problem for nonlinear elliptic
equations with VMO coeffcients}, Trans. Amer. Math. Soc., \textbf{336} (2) (1993),
841--853.
\bibitem{DaJSK11} Dan\v{e}\v{c}ek, J.; John, O.; Star\'{a}, J.;
\emph{Morrey space regularity for weak solutions of Stokes systems with VMO
coefficients}, Annali di Matematica, \textbf{190} (2011), 681--701.
\bibitem{DiK13} Diening, L.; Kaplick\'y, P.;
\emph{$L^q$ theory for a generalized Stokes system}, Manuscripta Math.,
\textbf{141} (1-2) (2013), 333--361.
\bibitem{Die14} Diening, L.;
\emph{Campanato estimates for the generalized Stokes System},
Anal. Mat. Pura. Appl., \textbf{193} (6) (2014), 1779--1794.
\bibitem{DiLR11} Diening, L.; Lengeler, D.; R$\mathring{u}$\v{z}i\v{c}ka, M.;
\emph{The stokes and poisson problem in variable exponent spaces}, Complex Var.
Elliptic Equ., \textbf{56} (7-9) (2011), 789--811.
\bibitem{DoK11} Dong, H. J.; Kim, D.;
\emph{On the $L_{p}$-solvability of higher order parabolic and elliptic systems
with BMO coefficients}, Arch. Ration. Mech. Anal., \textbf{199} (2011), 889--941.
\bibitem{DoK16} Dong, H. J.; Kim, D.;
\emph{$L_{q}$-estimates for stationary Stokes system with coefficients measurable
in one direction}, arXiv:1604.02690, 25 pages.
\bibitem{GuS15} Gu, S.; Shen, Z. W.;
\emph{Homogenization of Stokes systems and uniform regularity estimates},
SIAM J. Math. Anal., \textbf{47} (5) (2015), 4025--4057.
\bibitem{Gia83} Giaquinta, M.;
\emph{Multiple integral in the calculus of variations and nonlinear elliptic systems},
Princeton Univ. Press, Princeton, 1983.
\bibitem{HuLW14} Hu, F. J.; Li, D. S.; Wang, L. H.; \emph{A new proof of $L^p$ estimates of Stokes equations}, J. Math. Anal. Appl., \textbf{420} (2014), 1251--1264.
\bibitem{Kr07} Krylov, N. V.;
\emph{Parabolic and elliptic equations with VMO coefficients},
Comm. Partial Differential Equations, \textbf{32} (1-3) (2007), 453--475.
\bibitem{Mac11} M\'{a}cha, V.;
\emph{On a generalized Stokes problem}, Cent. Eur. J. Math., \textbf{9} (2011),
874--887.
\bibitem{MeP12} Mengesha, T.; Phuc, N. C.;
\emph{Global estimates for quasilinear elliptic equations on Reifenberg flat domains},
Arch. Rational Mech. Anal., \textbf{203} (2012), 189--216.
\bibitem{Min07} Mingione, G.;
\emph{The Calder\'on-Zygmund theory for elliptic problems with measure data},
Ann. Sc. Norm. Super. Pisa. Cl. Sci., \textbf{6} (5) (2007), 195--261.
\bibitem{Min10} Mingione, G.;
\emph{Gradient estimates below the duality exponent},
Math. Ann., \textbf{346} (2010), 571--627.
\bibitem{Min11} Mingione, G.;
\emph{Gradient potential estimates}, J. Eur. Math. Soc., \textbf{13} (2011),
459--486.
\bibitem{TiZh17} Tian, H.; Zheng, S. Z.;
\emph{Lorentz estimates for the gradient of weak solutions to elliptic obstacle
problems with partially BMO coefficients}, Bound. Value Probl.,
\textbf{128} (2017), 1--27.
\bibitem{ZhZ16} Zhang, J. J.; Zheng, S. Z.;
\emph{ Lorentz estimates for asymptotically regular elliptic equations in
quasiconvex domains}, Electron. J. Differ. Equ., \textbf{2016 } (142) (2016), 1--13.
\bibitem{ZhZ16} Zhang, J. J.; Zheng, S. Z.;
\emph{Lorentz estimates for fully nonlinear parabolic and elliptic equations},
Nonlinear Anal., \textbf{148} (2017), 106--125.
\bibitem{ZhZ17} Zhang, J. J.; Zheng, S. Z.;
\emph{Weighted Lorentz estimates for nondivergence linear elliptic equations
with partially BMO coefficients}, Commun. Pure Appl. Anal.,
\textbf{16} (3) (2017), 899--914.
\bibitem{ZhZ14} Zhang, C.; Zhou, S. L.;
\emph{Global weighted estimates for quasilinear elliptic equations with
non-standard growth}, J. Funct. Anal., \textbf{267} (2014), 605--642.
\end{thebibliography}
\end{document}