\section{Introduction and main results}
\label{sec:intro}
The main goal of the present paper is to obtain the gradient estimates of solutions to following quasilinear divergence form elliptic equations:
\begin{align}\label{eq:diveq}
\begin{cases}
\mathrm{div}(A(x,\nabla u)) &= \ \mathrm{div}(f) \quad \text{in} \ \ \Omega, \\
\hspace{1.2cm} u &= \ g \quad \qquad \text{on} \ \ \partial \Omega,
\end{cases}
\end{align}
in the weighted Lorentz spaces. Moreover, the problem is considered under the assumption that domain $\Omega$ is Reifenberg flat domain of $\mathbb{R}^n$ ($n \ge 2$), and the small BMO (Bounded Mean Oscillation) condition on $A$. Here, the nonlinearity $A$ is a Carath\'edory vector valued function only satisfying the growth and monotonicity conditions: there holds
\begin{align*}
\left| A(x,\xi) \right| &\le \Lambda_1 |\xi|^{p-1},\\
\langle A(x,\xi)-A(x,\eta), \xi - \eta \rangle &\ge \Lambda_2 \left( |\xi|^2 + |\eta|^2 \right)^{\frac{p-2}{2}}|\xi - \eta|^2,
\end{align*}
for every $(\xi,\eta) \in \mathbb{R}^n \times \mathbb{R}^n \setminus \{(0,0)\}$ and a.e. $x \in \mathbb{R}^n$, $\Lambda_1$ and $\Lambda_2$ are positive constants. This operator and its properties are emphasized in Section~\ref{sec:pre}. Additionally, boundary data $g \in W^{1,p}(\Omega;\mathbb{R})$ for some $p>1$ and $f \in L^{p'}$, $p'$ is the exponent conjugate to $p$.

The Reifenberg flat domain was first introduced in 1960
by E. Reifenberg \cite{ER60} when he studied on the solution of the Plateau problem. It appears and plays an important role in the theory of minimal surface and free boundary problems. In the Calder\'on-Zygmund theory for various elliptic and parabolic partial differential equations, when the roughness of the domain is discussed or some boundary regularity results more than H\"older is expected, the Reifenberg domains are then studied and considered.  For instance, the regularity of the elliptic kernel (for example, Poisson kernel) and harmonic measure are earlier shown in many fine papers by C. Kenig and T. Toro in~\cite{KT1997,Toro1997,KT1999,KT2003}. Geometrically, Reifenberg domains include domains with rough fractal boundaries, that can be confined by two hyperplanes inside and outside of the domains. The Reifenberg domains are \emph{flat} in the sense that their boundaries can be well approximated by planes/hyperplanes and they are geometrically invariant with respect to the chosen scale. This class of domains include all $C^1$-domains, Lipschitz domains with small Lipschitz constants, and domains with fractal boundaries. There are also many published papers working on Reifenberg flat domains, their properties and applications, that can be found in~\cite{LM3,LM2,LM1,DT1999,MT2010} and related references therein. 

Regularity of solutions to elliptic problems with homogeneous Dirichlet problem (zero Dirichlet boundary data) is an interesting topic. The regularity theory in various function spaces has attracted attention of mathematicians for many years. Under various assumptions of domain $\Omega$, the vector field $A$ and functional data $f,g$, there have been increasing interests in studying on gradient estimates of solutions. The reader can find a plenty of materials related to this topic by E. DiBenedetto in~\cite{DiBenedetto1983, Lieberman1984, Lieberman1986, Lieberman1988,DiBenedetto1993}, by T.~Iwaniec in~\cite{Iwaniec}, S.~S. Byun et al. in~\cite{SSB2, BW1_1, BW2, SSB4} and further generalization to this type of homogeneous equation are the subjects of~\cite{CM2014, Phuc2015, CoMi2016, BCDKS, BW1, KZ} and their related references. Later, more general extensions of regularity to the non-homogeneous quasilinear elliptic equations of the form $\mathrm{div}(A(x,\nabla u)) = \mathrm{div}(f)$ in Lorentz-type spaces, as well as in Morrey-type spaces were discussed and addressed in many papers, such as~\cite{Truyen2018,Tuoc2018,MPTN2019}. 


Our aim in this paper is that to present the estimates and point-wise estimates on gradient of solutions to~\eqref{eq:diveq} in weighted Lorentz spaces, where $\Omega$ is Reifenberg flat domain, the vector field $A$ satisfies the growth, monotonicity conditions with small BMO. More precisely, our approach relies on the gradient estimates and point-wise estimates for solutions in terms of fractional maximal functions and Riesz potentials. To our knowledge, the technique using maximal functions was first presented by F. Duzaar et al. in their famous papers~\cite{Duzamin2,55DuzaMing}, and later there have been a lot of works developed following this introduced method. The main tool of our work is that we rely on the methods introduced in~\cite{55QH4} and those of our previous papers in~\cite{MPT2018,MPT2019,MPTN2019} relate to this subject. The results obtained here are comparable with such an approach in~\cite{Truyen2018}, but in our present paper, there is no differentiation assumed for the nonlinearity $A(x,\xi)$ with respect to $\xi$. It also emphasizes that in our regularity results, we obtain the weighted Lorentz gradient norm estimates in term of fractional maximal functions ${\mathbf{M}_\alpha}$, for $0 \le \alpha<n$, stronger and more general than those considered in previous works. 

Let us now give precise statements of our main results, via some main theorems presented as below.The first theorem is also called the ``good-$\lambda$'' theorem for weighted context. For the unweighted case, the proof of this theorem was well-proved in our previous work~\cite{MPTN2019}.

\begin{taggedtheorem}{A}\label{theo:goodlambda}
Let $p>1$, $\alpha \in [0,n)$, $\omega \in \mathcal{A}_\infty$, $f \in L^{\frac{p}{p-1}}(\Omega; \mathbb{R}^n)$, $g \in W^{1,p}(\Omega;\mathbb{R})$ and $u$ be a weak solution to \eqref{eq:diveq}. For any $\varepsilon>0$, $\lambda>0$ and $R_0>0$, there exist some constants $\delta = \delta(n,\varepsilon,[\omega]_{\mathcal{A}_\infty})$, $s_0 = s_0(n,p,\Lambda_1,\Lambda_2)$, ${\sigma} = \sigma(n,p,\alpha)$ and $\kappa = \kappa(n,p,\alpha,\varepsilon, [\omega]_{\mathcal{A}_\infty}, \mathrm{diam}(\Omega)/R_0)$ such that if $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain satisfying $[A]_{s_0}^{R_0} \le \delta$, then
\begin{multline}\label{eq:lambdabound}
\omega\left(\{\mathbf{M}\mathbf{M}_{\alpha}(|\nabla u|^{p}))>{\sigma}\lambda, \mathbf{M}_{\alpha}(|f|^{\frac{p}{p-1}}+|\nabla g|^p) \le \kappa\lambda \}\cap \Omega \right) \\ \leq C \varepsilon \omega\left(\{\mathbf{M} \mathbf{M}_{\alpha}(|\nabla u|^{p})> \lambda\}\cap \Omega \right).
\end{multline}
Here, we note that the constant $C$ depends only on $n,p,\alpha, \varepsilon,\mathrm{diam}(\Omega)/R_0,[\omega]_{\mathcal{A}_\infty}$.
\end{taggedtheorem}

Throughout the paper, the denotation $diam(\Omega)$ is the diameter of a set $\Omega$ defined as:
\begin{align*}
diam(\Omega) = \sup\{d(x,y) \  : \ x,y \in \Omega\},
\end{align*}
and $\mathcal{A}_\infty$ the Muckenhoupt weights will be described in Section \ref{sec:pre} later. Moreover, here and hereforth, for simplicity, the set $\{x \in \Omega: |g(x)| > \Lambda\}$ is denoted by $\{|g|>\Lambda\}$. The following theorem \ref{theo:regularity} establishes the estimates on gradients of solutions in term of fractional maximal functions.

\begin{taggedtheorem}{B}\label{theo:regularity}
Let $p>1$,  $f \in L^{\frac{p}{p-1}}(\Omega; \mathbb{R}^n)$, $g \in W^{1,p}(\Omega;\mathbb{R})$ and $u$ be a weak solution to \eqref{eq:diveq}. For any $\alpha \in [0,n)$, $\omega \in \mathcal{A}_\infty$, $0<q<\infty$, $0<s\le \infty$, there exist some constants $\delta = \delta(n,\varepsilon,[\omega]_{\mathcal{A}_\infty})$, $s_0 = s_0(n,p,\Lambda_1,\Lambda_2)$ such that if $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain satisfying $[A]_{s_0}^{R_0} \le \delta$ for some $R_0>0$, then
\begin{align}\label{eq:reg}
\|\mathbf{M}\mathbf{M}_{\alpha}(|\nabla u|^p)\|_{L^{q,s}_{\omega}(\Omega)} \le C \|\mathbf{M}_{\alpha}(|f|^{\frac{p}{p-1}} + |\nabla g|^p)\|_{L^{q,s}_{\omega}(\Omega)}.
\end{align}
The constant $C$ here depends only on $n,p,\alpha, \varepsilon,\mathrm{diam}(\Omega)/R_0,[\omega]_{\mathcal{A}_\infty}$.
\end{taggedtheorem}

In recent papers~\cite{KMin2013,KMin2014,Mi4}, T. Kuusi and G. Mingione firstly proved the point-wise gradient estimates for solutions to elliptic equations with measure data, in terms of the so-called Riesz potential. Our motivations here arise from their research papers, and that the problem of deriving point-wise estimates needs the appearance of linear/nonlinear potentials. These give rise to action, we also prove in this paper the point-wise gradient estimates of solutions in term of Riesz potential $\mathbf{I}_\beta$ stated via Theorem~\ref{theo:I-alpha} as follows.

\begin{taggedtheorem}{C}\label{theo:I-alpha}
Let $p>1$,  $f \in L^{\frac{p}{p-1}}(\Omega; \mathbb{R}^n)$, $g \in W^{1,p}(\Omega;\mathbb{R})$ and $u$ be a weak solution to \eqref{eq:diveq} in a $(\delta,R_0)$-Reifenberg flat domain $\Omega$ for sufficiently small $\delta>0$, with $[A]_{s_0}^{R_0} \le \delta$ for some $R_0, s_0>0$. Then for any $\alpha \in [0,n)$, $\beta \in (0,n)$, $0<q<\infty$, the following point-wise estimate
\begin{align}\label{eq:I-alpha}
\mathbf{I}_{\beta}\left(|\mathbf{M}_{\alpha}(|\nabla u|^p)|^q \chi_{\Omega}\right)(x) \le C \mathbf{I}_{\beta}\left(|\mathbf{M}_{\alpha}(|f|^{\frac{p}{p-1}} + |\nabla g|^p)|^q \chi_{\Omega} \right)(x)
\end{align}
holds for almost everywhere $x \in \mathbb{R}^n$.
\end{taggedtheorem}

Furthermore, the problem context also includes many applications. We believe that the method introduced in our work can also be extended to address the study on existence of solutions for some more general problems of this type. Particularly, in this manuscript, as an application, we are concerned with the existence of solutions to equations of the type:
\begin{align}\label{eq:existI}
\begin{cases} -\mathrm{div}(A(x,\nabla u)) & = \ \mathbf{I}_{\beta}(|\nabla u|^p)^q + \mathrm{div}(f), \ \mbox{ in } \ \Omega, \\ \hspace{1.5cm} u & = \  g, \ \mbox{ on } \partial \Omega,\end{cases}
\end{align}
and our proof rests on the well-known Riesz potentials and the Riesz capacity condition. Let us now state the existence theorem that we are currently interested in. 
\begin{taggedtheorem}{D}\label{theo:existence_Ialpha}
Let $\alpha, \, \beta \in (0,n)$, $p>1$, $\frac{p-1}{p}<q<\frac{n}{n-\beta}$, $f \in L^{\frac{p}{p-1}}(\Omega; \mathbb{R}^n)$ and $g \in W^{1,p}(\Omega;\mathbb{R})$. There exist some positive constants $\delta$, $s_0$, $\varepsilon$ such that if $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain satisfying $[A]_{s_0}^{R_0} \le \delta$ for some $R_0>0$ and the functional data $|\mathcal{F}|^p: = |f|^{\frac{p}{p-1}} + |\nabla g|^p$ satisfies the following inequality
\begin{align}\label{eq:ReiszCapa_cond}
\mu(K) \le \varepsilon \ \mathrm{Cap}_{\mathbf{I}_{\beta + \frac{1}{q}},\frac{pq}{pq-p+1}}(K),
\end{align}
for any compact set $K \subset \mathbb{R}^n$ with $d\mu  = |\mathcal{F}(x)|^p dx$, then the equation~\eqref{eq:existI} admits at least a solution $u \in W^{1,p}(\Omega)$ and there holds
\begin{align}\label{eq:theo-D}
\mathbf{I}_{\alpha}(|\nabla u|^p)(x) \le \Lambda \mathbf{I}_{\alpha} (|\mathcal{F}|^p)(x), \quad  \mbox{ for a.e. } x \in \mathbb{R}^n,
\end{align}
for a constant $\Lambda>0$.
\end{taggedtheorem}

 The mentioned condition \eqref{eq:ReiszCapa_cond} is well-known in \cite{V2} as the so-called Riesz capacity condition, where the $(\alpha,p)$-capacity ${\mathrm{Cap}}_{\mathbf{I}_\alpha,p}(K)$ corresponds to the Sobolev spaces $W^{\alpha,p}(\mathbb{R}^n)$ of a compact set $K$ is defined by
\begin{align*}
{\mathrm{Cap}}_{\mathbf{I}_\alpha,p}(K) = \inf \left\{\int_{\mathbb{R}^n}{|\phi(x)|^pdx}: \ \phi \in L_+^p(\mathbb{R}^n), \, \mathbf{I}_\alpha[\phi] \ge \chi_K \right\}.
\end{align*}
It refers to the most popular book authored by D. R. Adams and L. I. Hedberg in~\cite{AH}, see also~\cite{V2,Phuc2008,55QH4,BNV2018} for further information.

It is of interest to find a necessary and sufficient condition to obtain the convergence result. Unfortunately, our sufficient condition~\eqref{eq:ReiszCapa_cond} in Theorem~\ref{theo:existence_Ialpha} does not adapt to the necessary condition~\eqref{eq:ReiszCapa_cond-b} which is considered in the next theorem. However, two these conditions will be the same in a special case when $\frac{1}{q} + \frac{\beta}{p} = 1$. 

\begin{taggedtheorem}{E}\label{theoE}
Let $\beta \in (0,n)$, $p>1$, $q > \max\big\{\frac{p-1}{p},\frac{\beta}{\beta+1}\big\}$ and $\mu \in \mathcal{M}^+(\Omega)$. There exist two positive constants $\delta$, $s_0$ such that  if $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain satisfying $[A]_{s_0}^{R_0} \le \delta$ for some $R_0>0$ and the following equation
\begin{align}\label{eq:existI-b}
\begin{cases} -\mathrm{div}(A(x,\nabla u)) & = \ \mathbf{I}_{\beta}(|\nabla u|^p)^q + \mu, \ \mbox{ in } \ \Omega, \\ \hspace{1.5cm} u & = \  0, \ \mbox{ on } \partial \Omega,\end{cases}
\end{align}
admits a renormalized solution $u$, then one can find a constant $C$ such that
\begin{align}\label{eq:ReiszCapa_cond-b}
\mu(K) \le C \ \mathrm{Cap}_{\mathbf{I}_{\beta + 1- \frac{\beta}{q}},\frac{pq}{pq-p+1}}(K),
\end{align}
for any compact set $K \subset \mathbb{R}^n$. 
\end{taggedtheorem}

The rest of our paper is organized as follows. We start in Section~\ref{sec:pre} by introducing/collecting some standard notations, assumptions in which our problem is formulated. Section~\ref{sec:comparison} is dedicated to the interior and boundary comparison estimates on the solutions, some preparatory lemmas in this section also present a basic idea that allows us to prove results. In Section~\ref{sec:main} we drive the so-called ``good-$\lambda$'' technique to obtain the gradient estimates for the fractional maximal operators and the point-wise gradient bounds for solutions in terms of Riesz potentials, the proofs of gradient estimate theorems are also given in this section. The last section~\ref{sec:app} is then devoted to proving Theorem~\ref{theo:existence_Ialpha} and Theorem~\ref{theoE}, an application that may interact with many mathematical or physical equations in many fields of science.

\section{Notations and Preliminaries}
\label{sec:pre}

This section consists of some necessary preliminaries in which our problem is formulated, and we also recall some well-known notations, fundamentals and results for later use herein.  

\subsection{Notations and definitions}
\label{subsec:def}

In the present paper, the considered domain $\Omega \subset \mathbb{R}^n$ is assumed to be a bounded $(\delta,R_0)$-Reifenberg flat domain, whose definition is stated as follows.
\begin{definition}[$(\delta,R_0)$-Reifenberg flat domain]\label{def:Reif} 
Let $\delta \in (0,1)$ and $R_0>0$. We say that $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain if for every $x\in \partial \Omega$ and $r \in (0,R_0]$, there exists a system coordinates $\{y_1,y_2,...,y_n\}$ such that in this coordinate $x$ coincides the original $0$ and
\begin{align*
B_r(0) \cap \{y_n > \delta r\} \subset B_r(0) \cap \Omega \subset B_r(0) \cap \{y_n > -\delta r\},
\end{align*}
where we denote the set $\{y = (y_1, y_2, ..., y_n): \ y_n > c\}$ by $\{y_n > c\}$.
\end{definition}

Throughout the study, we recall that the denotation $B_r(x)$ stands for an open ball in $\mathbb{R}^n$ with radius $r$ and centered at $x$, that is the set $B_r(x) = \{y\in \mathbb{R}^n: |y-x|<r\}$. For convenience of the reader, we use $|B|$ stands for the $n$-dimensional Lebesgue measure of a set $B \subset \mathbb{R}^n$. And in what follows, let us denote by $\displaystyle{\fint_{B_r(x)}{f(y)dy}}$ indicates the integral average (mean value) of $f$ in the variable $y$ over the ball $B_r(x)$, i.e.
\begin{align*}
\fint_{B_r(x)}{f(y)dy} = \frac{1}{|B_r(x)|}\int_{B_r(x)}{f(y)dy}.
\end{align*}

In the setting of equation $\mathrm{div}(A(x,\nabla u))= \ \mathrm{div}(f)$, the nonlinear operator $A: \mathbb{R}^n \times \mathbb{R}^n \rightarrow \mathbb{R}$ is a Carath\'eodory vector valued function (that is, $A(.,\xi)$ is measurable on $\Omega$ for every $\xi$ in $\mathbb{R}^n$, and $A(x,.)$ is continuous on $\mathbb{R}^n$ for almost every $x$ in $\Omega$) which satisfies the following growth and monotonicity conditions: for some $1<p\le n$ there exist two positive constants $\Lambda_1$ and $\Lambda_2$ such that 
\begin{align}
\label{cond:A1}
\left| A(x,\xi) \right| &\le \Lambda_1 |\xi|^{p-1},
\end{align}
and
\begin{align}
\label{cond:A2}
\langle A(x,\xi)-A(x,\eta), \xi - \eta \rangle &\ge \Lambda_2 \left( |\xi|^2 + |\eta|^2 \right)^{\frac{p-2}{2}}|\xi - \eta|^2
\end{align}
holds for almost every $x$ in $\Omega$ and every $(\xi,\eta) \in \mathbb{R}^n \times \mathbb{R}^n \setminus \{(0,0)\}$.

Moreover, in our regularity proofs, the operator $A$ is also assumed to satisfy a sufficiently small bounded mean oscillation (BMO) condition, that is described as below.

\begin{definition}[$(\delta,R_0)$-BMO condition]\label{def:BMOcond}
Let $\delta >0$ and $R_0>0$. We say that the nonlinearity $A$ satisfies a $(\delta,R_0)$-BMO condition with exponent $s>0$ if
\begin{align*
[A]_s^{R_0} = \sup_{y \in \mathbb{R}^n, \ 0<r\le R_0} \left(\fint_{B_r(y)} \left(\sup_{\xi \in \mathbb{R}^n \setminus \{0\}} \frac{|A(x,\xi) - \overline{A}_{B_r(y)}(\xi)|}{|\xi|^{p-1}}\right)^s dx\right)^{\frac{1}{s}} \le \delta,
\end{align*}
where $\overline{A}_{B_r(y)}(\xi)$ denotes the average of $A(\cdot,\xi)$ over the ball $B_r(y)$.
\end{definition}

Our work is also related to the class of Muckenhoupt's weights $\mathcal{A}_p$. This concept first appeared by B. Muckenhoupt in \cite{Muck72} and since then, numerous norm inequalities and boundedness of relevant operators have been established for the $\mathcal{A}_p$ classes in various research approaches. The Muckenhoupt classes of weighted functions are closely connected with the boundedness of Hardy-Littlewood maximal functions. Let us recall the definition of the Muckenhoupt weights and derive some of their properties for later use. 

Here and subsequently, by a weight $\omega$, we mean that $\omega$ is a non-negative measurable and locally integrable function on $\mathbb{R}^n$. For any measurable set $E\subset \mathbb{R}^{n}$ and the weight $\omega$, we denote 
$$\omega(E):=\int_{E}\omega(x)dx.$$ 

\begin{definition}[Muckenhoupt weights]\label{def:Muck}
For $1\le p \le \infty$, we say that a weight $\omega \in L^1_{\mathrm{loc}}(\mathbb{R}^n)$ belongs to the Muckenhoupt class $\mathcal{A}_p$ if there holds
\begin{align*}
[\omega]_{\mathcal{A}_p} = \sup_{B_r(x) \subset \mathbb{R}^n} \left(\fint_{B_r(x)} \omega(y) dy\right)\left(\fint_{B_r(x)}\omega(y)^{-\frac{1}{p-1}}dy\right)^{p-1}<\infty, \mbox{ when } 1<p<\infty,
\end{align*}
\begin{align*}
[\omega]_{\mathcal{A}_1} = \sup_{B_r(x) \subset \mathbb{R}^n} \left(\fint_{B_r(x)} \omega(y) dy\right) \sup_{y \in B_r(x)} \frac{1}{\omega(y)} <\infty, \mbox{ when } p=1,
\end{align*}
and there are two positive constants $C$ and $\nu$ such that
\begin{align*}
\omega(E) \le C \left(\frac{|E|}{|B|}\right)^\nu \omega(B), \mbox{ when } p=\infty,
\end{align*}
for all ball $B=B_r(x)$ in $\mathbb{R}^n$ and all measurable subset $E$ of $B$. In this case, we denote $[\omega]_{\mathcal{A}_\infty} = (C,\nu)$.
\end{definition}

\begin{remark}\label{rem:muck}
In Definition~\ref{def:Muck}, the number $[\omega]_{\mathcal{A}_p}$  is called the ${\mathcal{A}_p}$ constant of $\omega$ and it is well known that $\mathcal{A}_1 \subset \mathcal{A}_p \subset \mathcal{A}_\infty$ for all $1 \le p \le \infty$. Moreover, the Muckenhoupt class $\mathcal{A}_{\infty}$ is given by: $$\mathcal{A}_\infty = \bigcup_{p<\infty}\mathcal{A}_p.$$ 
\end{remark}

In this paper, the study will be made in the setting of weighted Lorentz spaces, defined as below. And for literature that concerning these spaces, the reader refers to~\cite{Lorentz1950,Lorentz1951,Gene98} and textbooks~\cite{55Gra,Stein} for detailed information.
 
\begin{definition}[Weighted Lorentz space]
Let $0<q<\infty$, $0<s\le \infty$ and the Muckenhoupt weight $\omega \in \mathcal{A}_{\infty}$. We define the weighted Lorentz space $L^{q,s}_{\omega}(\Omega)$ by the set of all Lebesgue measurable functions $h$ on $\Omega$ such that $\|h\|_{L^{q,s}_{\omega}(\Omega)}<+\infty$, where
\begin{align}\label{eq:w-lorentz-a}
\|h\|_{L^{q,s}_{\omega}(\Omega)} = \left[ q \int_0^\infty{ \lambda^s\omega \left( \{x \in \Omega: |h(x)|>\lambda\} \right)^{\frac{s}{q}} \frac{d\lambda}{\lambda}} \right]^{\frac{1}{s}} < +\infty,
\end{align}
if $s < \infty$ and 
\begin{align}\label{eq:w-lorentz-b}
\|h\|_{L^{q,\infty}_{\omega}(\Omega)} = \sup_{\lambda>0}{\lambda \omega\left(\{x \in \Omega:|h(x)|>\lambda\}\right)^{\frac{1}{q}}}<+\infty,
\end{align}
if $s=\infty$.
\end{definition}
In this way, when $\omega=1$, the weighted Lorentz space $L^{q,s}_{\omega}(\Omega)$ becomes the unweighted (classical) Lorentz space $L^{q,s}(\Omega)$. Moreover, in the case of weighted Lorentz spaces, when $q = s$, $L^{q,s}_{\omega}(\Omega)$ coincides the weighted Lebesgue space $L^{q}_{\omega}(\Omega)$ which is defined by the set of all measurable functions $h$ such that
\begin{align*}
\|h\|_{L^{q}_{\omega}(\Omega)} := \left(\int_{\Omega} |h(x)|^q\omega(x) dx\right)^{\frac{1}{q}} < + \infty.
\end{align*}

\begin{definition}[Riesz potential]
Let $n \ge 2$ and the Riesz potential $\mathbf{I}_\beta$ of order $\beta \in (0,n)$  of a measurable function $h \in L^1_{\mathrm{loc}}(\mathbb{R}^n;\mathbb{R}^+)$ is defined as the convolution:
\begin{align}
\label{def:Riesz}
\mathbf{I}_\beta (h)(x) \equiv (\mathbf{I}_\beta * h)(x) = \int_{\mathbb{R}^n}{\frac{h(y)}{|x-y|^{n-\beta}}dy}, \quad x \in \mathbb{R}^n.
\end{align}
\end{definition}

\begin{definition}[Wolff potential]
Let $\alpha \in (0,n)$ and $1 < \beta <\frac{n}{\alpha}$. The Wolff potential $\mathbf{W}_{\alpha,\beta}(\nu)$ of a non-negative Borel measure $\nu$ is defined as the convolution:
\begin{align*
\mathbf{W}_{\alpha,\beta} (\nu)(x) = \int_0^{\infty} \left(\frac{\nu(B_r(x))}{r^{n-\alpha\beta}}\right)^{\frac{1}{\beta -1}} \frac{dr}{r}, \quad x \in \mathbb{R}^n.
\end{align*}
\end{definition}
We write $\mathbf{W}_{\alpha,\beta} (h)$ instead of $\mathbf{W}_{\alpha,\beta} (\nu)$ if $d\nu = hdx$, where $h \in L^1_{\mathrm{loc}}(\mathbb{R}^n;\mathbb{R}^+)$. We also remark that $\mathbf{I}_{\alpha}(\nu) = \mathbf{W}_{\frac{\alpha}{2},2}(\nu)$. 


\subsection{Fractional Maximal functions}\label{sec:MMalpha}

We now recall the definition of fractional maximal function that regarding to~\cite{K1997, KS2003}. Let $0 \le \alpha \le n$, the fractional  maximal function $\mathbf{M}_\alpha$ of a locally integrable function $h \in L^1_{\mathrm{loc}}(\mathbb{R}^n;\mathbb{R})$ is defined by:
\begin{align*
\mathbf{M}_\alpha h(x) = \sup_{\rho>0}{\rho^\alpha \fint_{B_\rho(x)}{|h(y)|dy}}, \qquad x \in \mathbb{R}^n.
\end{align*}
It is worth to remark that for the case $\alpha=0$, the fractional  maximal function $\mathbf{M}_\alpha$ becomes the Hardy-Littlewood maximal function $\mathbf{M}$. The standard and classical properties of the maximal function $\mathbf{M}$ can be found in many places, see for instance~\cite{Gra97,55Gra}, and later also in~\cite{55QH10}. Here we recall some well-known properties of maximal and fractional maximal operators, that will be shown in some following lemmas.
\begin{lemma}[The principal result about maximal function, see \cite{55Gra}] \label{lem:boundM}
The maximal function $\mathbf{M}$ is bounded from $L^q(\mathbb{R}^n)$ to $L^{q,\infty}(\mathbb{R}^n)$, for $q \ge 1$, i.e., there exists a positive constant $C$ such that
\begin{align*}
\left|\{x \in \mathbb{R}^n: \ \mathbf{M}h(x)>\lambda\}\right| \le \frac{C}{\lambda^q}\int_{\mathbb{R}^n}{|h(x)|^q dx}, 
\end{align*}
for all $\lambda>0$ and  $h \in L^q(\mathbb{R}^n)$.
\end{lemma}
\begin{lemma}[see \cite{55Gra}]\label{lem:boundMlorentz}
Let $q>1$ and $0 < s \le \infty$, there exists a positive constant $C$ such that
\begin{align*
\|\mathbf{M}h\|_{L^{q,s}(\Omega)} \le C \|h\|_{L^{q,s}(\Omega)},
\end{align*}
for all $h \in L^{q,s}(\Omega)$.
\end{lemma}
Moreover, a very important property of fractional maximal function was also obtained from the boundedness property of maximal function. The detail proof of this result can be found in~\cite{MPT2019}.
\begin{lemma}\label{lem:Malpha}
Let $0\le \alpha<n$ and for any locally integrable function $h \in L^1_{\mathrm{loc}}(\mathbb{R}^n)$ there holds
\begin{equation*}
{\left|\left\{x \in \mathbb{R}^n: \ \mathbf{M}_\alpha h(x) >\lambda\right\}\right|} \leq \frac{C}{\lambda^{\frac{n}{n-\alpha}}} \left(\int_{\mathbb{R}^n}|h(x)|dx\right)^{\frac{n}{n-\alpha}},
\end{equation*}
for all $\lambda>0$.
\end{lemma}
Due to the importance of fractional maximal operators, recently in~\cite{MPT2019}, we define an additional cut-off maximal function ${\mathbf{M}}^r_{\alpha}$ of a locally integrable function $h$ that corresponding to ${\mathbf{M}}_{\alpha}$ as follows: for $0\le \alpha \le n$ and $r>0$, 
\begin{align*
{\mathbf{M}}^r_{\alpha}h(x)  = \sup_{0<\rho<r} \rho^{\alpha} \fint_{B_\rho(x)}h(y)dy, \qquad x \in \mathbb{R}^n.
\end{align*}
In the proof-of-work of the same paper~\cite{MPT2019}, we are concerned with an interesting property of the cut-off maximal function ${\mathbf{M}}^r_{\alpha}$, which we state in lemma below. This leads us to the key tools to achieve results in the sequel. We refer to~\cite{MPT2019} for the proof of this lemma.
\begin{lemma}\label{lem:MrMr}
Let $0\le \alpha <n$ and $r>0$. There exists a constant $C=C(n,\alpha)>0$ such that
\begin{align}\label{eq:MrMr}
{\mathbf{M}}^r {\mathbf{M}}^r_{\alpha} h(x) \le C {\mathbf{M}}^{2r}_{\alpha}h(x), \mbox{ for all } x \in \mathbb{R}^n,
\end{align}
for any $h \in L^1_{\mathrm{loc}}(\mathbb{R}^n)$.
\end{lemma}

\section{Comparison results}
\label{sec:comparison}

In this section, we present some local interior and boundary comparison estimates for weak solution $u$ of~\eqref{eq:diveq} that are essential to our development later. It is also remarkable that in some statements below and in what follows, we shall adopt the denotation $C$ for a suitable positive constant that is not necessary the same from line to line in each occurrence.

First of all, we can exploit the following integral estimate on gradient of solution $u$ to~\eqref{eq:diveq}, with respect to initial data $f$ and $g$.

\begin{proposition}\label{prop1}
Let $g \in W^{1,p}(\Omega), \ f \in L^{\frac{p}{p-1}}(\Omega)$ and $u$ be a weak solution of~\eqref{eq:diveq}. There exists a positive constant $C = C(n,p,\Lambda_1,\Lambda_2)$ such that
\begin{equation}\label{eq:prop1}
\int_{\Omega} |\nabla u|^p dx \le C \int_{\Omega} \left(|f|^{\frac{p}{p-1}} + |\nabla g|^p \right)dx.
\end{equation}
\end{proposition}
\begin{proof}
By using $u - g$ as a test function of equation~\eqref{eq:diveq}, we obtain
\begin{equation*}
\int_{\Omega} A(x,\nabla u) \nabla u dx = \int_{\Omega} A(x,\nabla u) \nabla g dx + \int_{\Omega} f \nabla (u - g) dx.
\end{equation*}
Taking into account both conditions of operator $A$ in~\eqref{cond:A1} and~\eqref{cond:A2}, it deduces that
\begin{equation*}
\int_{\Omega} |\nabla u|^p dx \le C \left( \int_{\Omega} |\nabla u|^{p-1} |\nabla g| dx + \int_{\Omega} |f| |\nabla u| dx + \int_{\Omega} |f| |\nabla g| dx \right).
\end{equation*}
Finally, we may easily obtain~\eqref{eq:prop1} by combining the H{\"o}lder and Young's inequalities from this estimate.
\end{proof}
\subsection{Interior estimates}
\label{subsec:int}
Let us fix a point $x_0 \in \Omega$ and $R>0$ such that $B_{2R}(x_0) \subset \Omega$. We derive herein a comparison estimate between the weak solution $u$ to \eqref{eq:diveq} and the unique solution $w$ to the homogeneous problem
\begin{equation}\label{eq:I1}
\begin{cases} \mathrm{div} A(x,\nabla w) & = \ 0, \qquad \, \mbox{ in } B_{2R}(x_0),\\ 
\hspace{1.2cm} w & = \ u - g, \ \mbox{ on } \partial B_{2R}(x_0)\end{cases}
\end{equation}
via the following lemmas.
\begin{lemma}\label{lem:I1}
Let $w$ be the unique solution to equation~\eqref{eq:I1}. Then there exists a positive constant $C=C(n,p,\Lambda_1,\Lambda_2)$ such that
\begin{align}\nonumber
\fint_{B_{2R}(x_0)} |\nabla u - \nabla w|^p dx  & \le   C \left( \fint_{B_{2R}(x_0)}  |f|^{\frac{p}{p-1}} + |\nabla g|^p dx  \right) \\ \label{eq:lem1b}
 & \qquad + C \left(\fint_{B_{2R}(x_0)}  |\nabla g|^p dx \right)^{\frac{1}{p}} \left(\fint_{B_{2R}(x_0)} |\nabla u|^{p} dx \right)^{\frac{p-1}{p}}.
\end{align}
\end{lemma}
\begin{proof}
Here, the proof is sketched in the similar way to that in~\cite[Lemma 4.3]{MPTN2019}. By choosing $u - w - \overline{g}$ as a test function of equations~\eqref{eq:diveq} and~\eqref{eq:I1}, where $\overline{g} = g$ in $\overline{B_{2R}}(x_0)$, we obtain that
\begin{align}\label{eq:L1I1}
\nonumber
\int_{B_{2R}(x_0)} \left(A(x,\nabla u) - A(x, \nabla w)\right) \nabla (u - w) dx  = \int_{B_{2R}(x_0)} \left(A(x,\nabla u) - A(x, \nabla w)\right) \nabla g dx \\ + \int_{B_{2R}(x_0)} f \nabla (u - w) dx  - \int_{B_{2R}(x_0)} f \nabla g dx.
\end{align}
Note moreover that by two conditions~\eqref{cond:A1} and~\eqref{cond:A2} of $A$, we deduce from~\eqref{eq:L1I1} that there exists a positive constant $C$ depending on $\Lambda_1, \Lambda_2$ such that
\begin{multline}\label{eq:L1I2}
\int_{B_{2R}(x_0)}  |\nabla u - \nabla w|^p dx \le C \left(\int_{B_{2R}(x_0)} |\nabla u|^{p-1} |\nabla g| dx + \int_{B_{2R}(x_0)} |\nabla w |^{p-1} |\nabla g| dx \right.  \\  \left.  + \int_{B_{2R}(x_0)} |f| |\nabla u - \nabla w| dx + \int_{B_{2R}(x_0)} |f|  |\nabla g| dx\right).
\end{multline}
Applying the following estimate
\begin{align*}
|\nabla w|^{p-1} & \le \left(|\nabla u| + |\nabla u - \nabla w|\right)^{p-1} \le 2^p \left(|\nabla u|^{p-1} + |\nabla u - \nabla w|^{p-1}\right),
\end{align*}
into~\eqref{eq:L1I2} one obtains
\begin{align}\nonumber
\int_{B_{2R}(x_0)}  |\nabla u - \nabla w|^p dx  \le  C & \left(\int_{B_{2R}(x_0)} |\nabla u|^{p-1} |\nabla g| dx + \int_{B_{2R}(x_0)} |\nabla u - \nabla w |^{p-1} |\nabla g| dx \right.  \\  \label{eq:L10}
 & \left.  + \int_{B_{2R}(x_0)} |f| |\nabla u - \nabla w| dx + \int_{B_{2R}(x_0)} |f|  |\nabla g| dx\right).
\end{align}
We use H{\"o}lder's inequality to estimate the first term on the right hand side
\begin{align}\label{eq:estI0}
\int_{B_{2R}(x_0)} |\nabla u|^{p-1} |\nabla g| dx  \le   \left(\int_{B_{2R}(x_0)}  |\nabla g|^p dx \right)^{\frac{1}{p}} \left(\int_{B_{2R}(x_0)} |\nabla u|^{p} dx \right)^{\frac{p-1}{p}}.
\end{align}
In order to estimate the other terms, we apply the following inequality~\eqref{eq:estHY} which is a consequence of H{\"o}lder's inequality and Young's inequality. More precisely, by applying H{\"o}lder and Young's inequalities, we obtain that for any $\varepsilon>0$, there exists $m = m(p,q,\varepsilon)>0$ such that
\begin{align}\label{eq:estHY}
\int_{B_{2R}(x_0)} |f_1 f_2| dx \le  \frac{\varepsilon}{2} \int_{B_{2R}(x_0)} |f_1|^p dx + m \int_{B_{2R}(x_0)} |f_2|^q dx,
\end{align}
where $q$ is the H{\"o}lder conjugate exponent corresponding to $p>1$. It is straightforward to see that these following inequalities hold:
\begin{align}\label{eq:estI1}
\int_{B_{2R}(x_0)} |\nabla u - \nabla w |^{p-1} |\nabla g| dx   \le \frac{\varepsilon}{2}  \int_{B_{2R}(x_0)} |\nabla u- \nabla w |^{p} dx + m_1 \int_{B_{2R}(x_0)} |\nabla g|^p dx,
\end{align}
\begin{align}\label{eq:estI2}
\int_{B_{2R}(x_0)} |f| |\nabla u - \nabla w| dx   \le  \frac{\varepsilon}{2} \int_{B_{2R}(x_0)} |\nabla u - \nabla w|^p dx + m_2 \int_{B_{2R}(x_0)} |f|^{\frac{p}{p-1}} dx,
\end{align}
and
\begin{align}\label{eq:estI3}
\int_{B_{2R}(x_0)} |f|  |\nabla g| dx \le  \int_{B_{2R}(x_0)} |\nabla g|^p dx + m_3 \int_{B_{2R}(x_0)} |f|^{\frac{p}{p-1}} dx,
\end{align}
where $m_1, m_2, m_3$ depend on $\varepsilon$ and $p$. For a chosen value $\varepsilon = 1$, these above estimates~\eqref{eq:L10}-\eqref{eq:estI3} conclude our interior estimate~\eqref{eq:lem1b} holds. Here, it is noticeable that the constant $C$ depends only on $n,p$ and two positive constants $\Lambda_1, \Lambda_2$ from conditions of $A$.
\end{proof}

We next recall in Lemma~\ref{lem:I1a} that the gradient of weak solution $w$ satisfies a well-known ``reverse'' H\"older integral inequality with increasing supports. It was discovered by Gehring (see~\cite{Gehring},~\cite[Theorem 6.7]{Giu}) and the applications of this inequality play an important role in development of mathematical analysis. Technique of using this inequality with small exponents to gradient estimates was proposed by G. Mingione et al. in their fine paper~\cite{Mi1} when the problem involves measure data. The reader is referred to~\cite{Phuc1, MPT2018, Mi3, 55QH4, HOk2019} and materials therein for the proof of this inequality and related results in similar research papers.

\begin{lemma}\label{lem:I1a}
Let $w$ be the unique solution to equation~\eqref{eq:I1}. Then there exist constants $\Theta_1 =\Theta_1(n,p,\Lambda_1,\Lambda_2) > p$ and $C=C(n,p,\Lambda_1,\Lambda_2)>0$ such that
\begin{equation*
\left(\fint_{B_{r}(y)} |\nabla w|^{\Theta_1} dx\right)^{\frac{1}{\Theta_1}}  \le C \left(\fint_{B_{2r}(y)} |\nabla w|^{p-1} dx\right)^{\frac{1}{p-1}},
\end{equation*}
for all $r>0$ and $y \in \Omega$ such that $B_{2r}(y) \subset B_{2R}(x_0)$.
\end{lemma}

The following lemma~\ref{lem:I2} states that in order to obtain the local bounds for solution $w$ of~\eqref{eq:I1}, the small BMO condition of $A$ is also stressed here. We refer the reader to~\cite{MP11} for the proof of this lemma.

\begin{lemma}\label{lem:I2}
Let $v$ be the unique solution to the following equation 
\begin{equation}\label{eq:I2}
\begin{cases} \mathrm{div} \overline{A}_{B_R(x_0)}(\nabla v) & = \ 0, \ \mbox{ in } B_R(x_0),\\ 
\hspace{1.2cm} v & = \ w, \ \mbox{ on } \partial B_R(x_0).\end{cases}
\end{equation}
Then there exist two positive constants $C_1=C_1(n,p,\Lambda_1,\Lambda_2)$ and $C_2=C_2(n,p,\Lambda_1,\Lambda_2)$ such that 
\begin{equation}\label{eq:lem2b}
C_1^{-1}\int_{B_R(x_0)} |\nabla v|^p dx \le \int_{B_{R}(x_0)} |\nabla w|^p dx \le C_1\int_{B_R(x_0)} |\nabla v|^p dx,
\end{equation}
and
\begin{equation}\label{eq:lem2a}
\fint_{B_R(x_0)} |\nabla w - \nabla v|^p dx \le C_2 \left[A^R_{s_1}\right]^p \fint_{B_{2R}(x_0)} |\nabla w|^p dx,
\end{equation}
where $\Theta_1$ is given in Lemma~\ref{lem:I1a} and $s_1 = \frac{p\Theta_1}{\Theta_1 - p}$. 
\end{lemma}

The comparison theorem between solution $v$ to~\eqref{eq:I2} and the weak solution $u$ to~\eqref{eq:diveq} can be stated and proved concisely, as indicated below.

\begin{theorem}\label{theo:I}
Let $v$ be the unique solution to equation~\eqref{eq:I2}. Assume moreover that $A$ satisfies a $(\delta,R)$-BMO condition with exponent $s_1$ given by Lemma~\ref{lem:I2}, for some constants $\delta\in(0,1)$. Then there exists a constant $C = C(n,p,\Lambda_1,\Lambda_2,R)$ such that 
\begin{align}\nonumber
\|\nabla v\|^p_{L^{\infty}(B_{R/2}(x_0))} & \le C \fint_{B_{2R}(x_0)} |\nabla u|^p  dx + C \fint_{B_{2R}(x_0)}  |f|^{\frac{p}{p-1}} + |\nabla g|^p dx  \\ \label{eq:theoIc}
& \qquad + C \left(\fint_{B_{2R}(x_0)} |\nabla g|^pdx \right)^{\frac{1}{p}}\left(\fint_{B_{2R}(x_0)} |\nabla u|^pdx \right)^{\frac{p-1}{p}} ,
\end{align}
and
\begin{multline}\label{eq:theoId}
\fint_{B_{R}(x_0)} |\nabla u - \nabla v|^pdx \le C \delta^p \fint_{B_{2R}(x_0)} |\nabla u|^p dx  + C   \fint_{B_{2R}(x_0)} |f|^{\frac{p}{p-1}} + |\nabla g|^p dx \\ 
+ C  \left(\fint_{B_{2R}(x_0)} |\nabla g|^pdx \right)^{\frac{1}{p}}  \left( \fint_{B_{2R}(x_0)} |\nabla u|^pdx \right)^{\frac{p-1}{p}} .
\end{multline} 
\end{theorem}
\begin{proof}
By~\eqref{eq:lem1b} in Lemma~\ref{lem:I1}, we have the following estimate
\begin{align} \nonumber
\fint_{B_R(x_0)} |\nabla w|^p dx & \le C \fint_{B_{2R}(x_0)} |\nabla u|^p  dx + C \fint_{B_{2R}(x_0)} |\nabla u - \nabla w|^p dx \\ \nonumber
& \le C \fint_{B_{2R}(x_0)} |\nabla u|^p  dx +  C \left( \fint_{B_{2R}(x_0)}  |f|^{\frac{p}{p-1}} + |\nabla g|^p dx  \right) \\ \label{eq:theoId_1}
 & \qquad + C \left(\fint_{B_{2R}(x_0)}  |\nabla g|^p dx \right)^{\frac{1}{p}} \left(\fint_{B_{2R}(x_0)} |\nabla u|^{p} dx \right)^{\frac{p-1}{p}}.
\end{align}
Moreover, by Lemma~\ref{lem:I1a} and~\eqref{eq:lem2b} in Lemma~\ref{lem:I2}, there holds
\begin{align} \label{eq:theoId_2}
\|\nabla v\|^p_{L^{\infty}(B_{R/2}(x_0))} &\le C \fint_{B_R(x_0)} |\nabla v|^p dx   \le C \fint_{B_R(x_0)} |\nabla w|^p dx. 
\end{align}
Combining this estimate together with the one inferred from~\eqref{eq:theoId_1}, we obtain~\eqref{eq:theoIc}. On the other hand, the following estimate is obtained from~\eqref{eq:lem2a} in Lemma~\ref{lem:I2},
\begin{align*}
\fint_{B_R(x_0)} |\nabla u - \nabla v|^p dx &\le C\left( \fint_{B_R(x_0)} |\nabla u - \nabla w|^p dx +  \fint_{B_R(x_0)} |\nabla w - \nabla v|^p dx \right) \\
& \le C\left( \fint_{B_R(x_0)} |\nabla u - \nabla w|^p dx + \left([A]_{s_1}^R\right)^p \fint_{B_R(x_0)} |\nabla w|^p dx \right).
\end{align*}
It concludes that the estimate~\eqref{eq:theoId} holds by using~\eqref{eq:theoId_2} and Lemma~\ref{lem:I1}, we also note here that $[A]_{s_1}^R \le \delta \in (0,1)$ by the small BMO condition.
\end{proof}

\subsection{Boundary estimates}
\label{subsec:boundary}
Now we come to the proof of the comparison estimates up to the boundary. Let us assume that $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain with $\delta \in (0,1/2]$, for $x_0 \in \partial \Omega$, $0<R<R_0/6$ and $w$ is the unique solution to the following problem
 \begin{equation}\label{eq:B1}
\begin{cases} \mathrm{div} A(x,\nabla w) & = \ 0, \ \mbox{ in } B_{6R}(x_0),\\ 
\hspace{1.2cm} w & = \ u - g, \ \mbox{ on } \partial B_{6R}(x_0).\end{cases}
\end{equation}

\begin{lemma}\label{lem:B2}
Let $w$ be the unique solution to equation~\eqref{eq:B1}. Then there exists a positive constant $C=C(n,p,\Lambda_1,\Lambda_2)$ such that
\begin{align}\nonumber
\fint_{B_{6R}(x_0)} |\nabla u - \nabla w|^p dx  & \le   C \left( \fint_{B_{6R}(x_0)}  |f|^{\frac{p}{p-1}} + |\nabla g|^p dx  \right) \\ \label{est_lem:B2}
 & \qquad + C \left(\fint_{B_{6R}(x_0)}  |\nabla g|^p dx \right)^{\frac{1}{p}} \left(\fint_{B_{6R}(x_0)} |\nabla u|^{p} dx \right)^{\frac{p-1}{p}}.
\end{align}
\end{lemma}

\begin{lemma}\label{lem:B1}
Let $w$ be the unique solution to equation~\eqref{eq:B1}. Then there exist constants $\Theta_2 = \Theta_2(n,p,\Lambda_1,\Lambda_2) > p$ and $C=C(n,p,\Lambda_1,\Lambda_2)>0$ such that
\begin{equation*
\left(\fint_{B_{\rho}(y)} |\nabla w|^{\Theta_2} dx\right)^{\frac{1}{\Theta_2}}  \le C \left(\fint_{B_{6\rho}(y)} |\nabla w|^{p-1} dx\right)^{\frac{1}{p-1}},
\end{equation*}
for all $B_{6\rho}(y) \subset B_{6R}(x_0)$.
\end{lemma}

Let us set $r = R(1-\delta)$. Since $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain with $\delta \in (0,1/2]$, there exists a coordinate system $\{y_1, y_2, ..., y_n\}$ with the origin $0 \in \Omega$ such that in this coordinate system $x_0$ and
\begin{align}\label{cond:B}
B_{r}^+(0) \subset \Omega \cap B_r(0) \subset B_r(0) \cap \{y_n > -4r\delta\},
\end{align}
where $B_{r}^+(0) = B_r(0) \cap \{y_n > 0\}$. 

\begin{lemma}\label{lem:B3}
Let $v$ be the unique solution to the following equation
\begin{equation}\label{eq:B3}
\begin{cases} \mathrm{div} \overline{A}_{B_r(0)}(\nabla v) & = \ 0, \ \mbox{ in } B_r(0),\\ 
\hspace{1.2cm} v & = \ w, \ \mbox{ on } \partial B_r(0).\end{cases}
\end{equation}
Then there exist positive constants $C_1=C_1(n,p,\Lambda_1,\Lambda_2)$ and $C_2 = C_2(n,p,\Lambda_1,\Lambda_2)$ such that 
\begin{equation}\label{est_lem:B3_a}
C_1^{-1}\int_{B_r(0)} |\nabla v|^p dx \le \int_{B_{r}(0)} |\nabla w|^p dx \le C_1\int_{B_r(0)} |\nabla v|^p dx,
\end{equation}
and
\begin{equation}\label{est_lem:B3_b}
\fint_{B_r(0)} |\nabla w - \nabla v|^p dx \le C_2 \left[A^R_{s_2}\right]^p \fint_{B_{r}(0)} |\nabla w|^p dx,
\end{equation}
where $\Theta_2$ given in Lemma~\ref{lem:B1} and $s_2 = \frac{p\Theta_2}{\Theta_2 - p}$. 
\end{lemma}

Because of the fact that $L^{\infty}$-norm of $\nabla v$ up to the boundary may not exist if $\partial \Omega$ is not regular enough, we consider $U$ as the weak solution to an another problem
\begin{equation}\label{eq:B4}
\begin{cases} \mathrm{div} \overline{A}_{B_r(0)}(\nabla U) & = \ 0, \ \mbox{ in } B_r^{+}(0),\\ 
\hspace{1.2cm} U & = \ 0, \ \mbox{ on }  B_r(0)\cap \{x_n=0\}.\end{cases}
\end{equation}
\begin{lemma}\label{lem:B5}
For any $\vartheta>0$, there exists $\delta_0 \in (0,1/2)$ such that if $v$ is a solution of~\eqref{eq:B3} under hypothesis~\eqref{cond:B} with $\delta \in (0,\delta_0)$, there exists a weak solution $U$ of~\eqref{eq:B4} satisfying
\begin{align}\label{est_lem:B5a}
\|\nabla U\|^p_{L^{\infty}(B_{r/4}(0))} &\le C \fint_{B_r(0)}|\nabla v|^p dx, \\ \label{est_lem:B5b}
\mbox{ and } \fint_{B_{r/8}(0)} |\nabla v - \nabla U|^p dx & \le \vartheta^p \fint_{B_r(0)} |\nabla v|^pdx,
\end{align}
for some constant $C=C(n,p,\Lambda_1,\Lambda_2)>0$.
\end{lemma}

\begin{theorem}\label{theo:B}
For any $\vartheta \in (0,1)$, there exists $\delta_0 \in (0,1/2)$ such that if $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain with $\delta \in (0,\delta_0)$ and $A$ satisfies a $(\delta,R_0)$-BMO condition, there exists a function $U$ such that
\begin{align}\nonumber
\|\nabla U\|^p_{L^{\infty}(B_{R/9}(x_0))} & \le C \fint_{B_{6R}(x_0)} |\nabla u|^p  dx + C \fint_{B_{6R}(x_0)}  |f|^{\frac{p}{p-1}} + |\nabla g|^p dx  \\ \label{est_theo:Ba}
& \qquad + C \left(\fint_{B_{6R}(x_0)} |\nabla g|^pdx \right)^{\frac{1}{p}}\left(\fint_{B_{6R}(x_0)} |\nabla u|^pdx \right)^{\frac{p-1}{p}} ,
\end{align}
and
\begin{multline}\label{est_theo:Bb}
\fint_{B_{R/9}(x_0)} |\nabla u - \nabla U|^pdx \le C \left(\vartheta^p + \delta^p\right) \fint_{B_{6R}(x_0)} |\nabla u|^p dx + C   \fint_{B_{6R}(x_0)} |f|^{\frac{p}{p-1}} + |\nabla g|^p dx \\ 
+ C  \left(\fint_{B_{6R}(x_0)} |\nabla g|^pdx \right)^{\frac{1}{p}}  \left( \fint_{B_{6R}(x_0)} |\nabla u|^pdx \right)^{\frac{p-1}{p}} .
\end{multline} 
\end{theorem}
\begin{proof} We first remark that we may assume $0 \in \Omega$ and there exists $\delta_0>0$ small enough such that for all $\delta \in (0,\delta_0)$ then condition~\eqref{cond:B} and the following estimate hold
\begin{align}\label{cond:B2}
B_{R/9}(x_0) \subset B_{r/8}(0) \subset B_{r/4}(0) \subset B_r(0) \subset B_{6r}(0) \subset B_{6R}(x_0),
\end{align}
where $r = R(1-\delta)$. Let $w$ and $v$ be solutions to~\eqref{eq:B1} and~\eqref{eq:B3} in Lemma~\ref{lem:B1} and Lemma~\ref{lem:B3} respectively. By Lemma~\ref{lem:B5}, for any $\vartheta>0$ there exists a small positive constant $\delta>0$ such that there exists a function $U$ satisfying~\eqref{est_lem:B5a} and~\eqref{est_lem:B5b}.

We may deduce from~\eqref{est_lem:B2} that
\begin{align} \nonumber
\fint_{B_{6R}(x_0)} |\nabla w|^p dx & \le C   \fint_{B_{6R}(x_0)} |\nabla u|^p dx + C  \fint_{B_{6R}(x_0)} |\nabla u - \nabla w|^p dx \\ \nonumber
& \le  C   \fint_{B_{6R}(x_0)} |\nabla u|^p dx + C \left( \fint_{B_{6R}(x_0)}  |f|^{\frac{p}{p-1}} + |\nabla g|^p dx  \right) \\ \label{est_lem:w}
 & \qquad + C \left(\fint_{B_{6R}(x_0)}  |\nabla g|^p dx \right)^{\frac{1}{p}} \left(\fint_{B_{6R}(x_0)} |\nabla u|^{p} dx \right)^{\frac{p-1}{p}}.
\end{align}
Combining~\eqref{est_lem:B5a} and~\eqref{est_lem:B3_a} with notice~\eqref{cond:B2}, one gets that
\begin{align*}
\|\nabla U\|^p_{L^{\infty}(B_{R/9}(x_0))}  \le C \fint_{B_{r}(0)} |\nabla v|^p dx \le C \fint_{B_{r}(0)} |\nabla w|^p dx  \le C \fint_{B_{6R}(x_0)} |\nabla w|^p dx,
\end{align*}
which follows to~\eqref{est_theo:Ba} from~\eqref{est_lem:w}. On the other hand, we obtain from~\eqref{cond:B2} that
\begin{align}\nonumber
\fint_{B_{R/9}(x_0)} |\nabla u - \nabla U|^pdx & \le C \fint_{B_{r/8}(0)} |\nabla u - \nabla U|^pdx \\ \nonumber
& \le C \fint_{B_{r/8}(0)} |\nabla u - \nabla w|^pdx + C \fint_{B_{r/8}(0)} |\nabla w - \nabla v|^pdx \\  \label{est_lem:final}
& \qquad + C \fint_{B_{r/8}(0)} |\nabla v - \nabla U|^pdx.
\end{align}
In order to obtain~\eqref{est_theo:Bb}, we need to estimate the second and third terms on the right hand side of~\eqref{est_lem:final}. The second term can be obtained by~\eqref{est_lem:B3_b} as follows
\begin{align*}
\fint_{B_{r/8}(0)} |\nabla w - \nabla v|^pdx \le C \fint_{B_{r}(0)} |\nabla w - \nabla v|^pdx \le C \delta^p \fint_{B_{6R}(x_0)} |\nabla w|^pdx.
\end{align*}
Thanks to~\eqref{est_lem:B5b} and~\eqref{est_lem:B3_a}, we may estimate the last one
\begin{align*}
\fint_{B_{r/8}(0)} |\nabla v - \nabla U|^pdx \le C\vartheta^p \fint_{B_{r}(0)} |\nabla v|^pdx \le C\vartheta^p \fint_{B_{6R}(x_0)} |\nabla w|^pdx.
\end{align*}
We can finally obtain~\eqref{est_theo:Bb} by taking into account these estimates, \eqref{est_lem:B2} and~\eqref{est_lem:w} to~\eqref{est_lem:final}.
\end{proof}

\section{Global gradient estimates}\label{sec:main}

We devote this section to proving our main results in theorems \ref{theo:goodlambda}, \ref{theo:regularity} and \ref{theo:I-alpha} that were stated in the introduction. 

To obtain the proof of theorem \ref{theo:goodlambda}, we require the following lemma \ref{lem:wReig}, the very important ingredient. The utility of this lemma normally relies on the Vitali type covering lemma (to cover a set $G$ by a countable family of pairwise disjoint closed balls) and the Lebesgue differentiation theorem (to control the size of the set on which the integral average can be large in terms of the $L^1$-norm), that are widespread in harmonic analysis. It refers to \cite{Vitali08}, the famous result in measure theory of Euclidean spaces, noticed by G. Vitali and later various literature concerning its modifications and applications \cite{CP1998,Wang2}.  The use of Vitali's covering lemma combining with maximal function techniques was first introduced by F. Duzaar and G. Mingione in \cite{Duzamin2, 55DuzaMing}. Further, several references \cite{SSB4,BW1_1,MP11} are also worth to read in solution estimates for elliptic and parabolic equations/systems.


\begin{lemma}\label{lem:wReig}
Let $\Omega$ be a $(\delta,R_0)$-Reifenberg flat domain with $\delta>0$ small enough and let $\omega$ be an $\mathcal{A}_\infty$ weight. Suppose that the sequence of balls $\{B_r(z_i)\}_{i=1}^N$ with center $z_i \in \bar{\Omega}$ and radius $r \le R_0/10$ covers $\Omega$. Let $V \subset W \subset \Omega$ be measurable sets for which there exists $0<\varepsilon<1$ such that:
\begin{itemize}
\item[(i)] $\omega(V) \le \varepsilon \omega(B_r(z_i))$ for all $i=1,2,...,N$;
\item[(ii)] $\omega(V \cap B_\rho(x)) \ge \varepsilon\omega(B_\rho(x)) \Rightarrow B_\rho(x) \cap \Omega \subset W$, for all $x \in \Omega, \rho \in (0,2r]$.
\end{itemize}
Then, there exists a constant $C = C(n,[\omega]_{\mathcal{A}_\infty})$ such that $\omega(V) \le C\varepsilon \omega(W)$.
\end{lemma}

\medskip

\begin{proof}[Proof of Theorem~\ref{theo:goodlambda}]
The proof of Theorem~\ref{theo:goodlambda} uses a technique similar to that appearing in \cite{55QH4} for the problem with measure data, but we later develop it to our framework here. Let us consider two sets 
\begin{equation*}
V = \left\{\mathbf{M}\mathbf{M}_{\alpha}(|\nabla u|^{p}))>{\sigma}\lambda, \mathbf{M}_{\alpha}(|f|^{\frac{p}{p-1}}+|\nabla g|^p) \le \kappa\lambda \right\}\cap \Omega, 
\end{equation*}
and
\begin{equation*}
W = \{\mathbf{M}\mathbf{M}_{\alpha}(|\nabla u|^{p})> \lambda\}\cap \Omega,
\end{equation*}
for any $\lambda>0$, where the constants ${\sigma}$, $\kappa$ will be specified later. Once having the Lemma \ref{lem:wReig} at hand, we outline the main steps that need to prove the sets $V$, $W$ satisfying all the assumptions, i.e., for any $\varepsilon>0$ there holds
\begin{enumerate}
\item[(i)] $\omega(V) \le C\varepsilon \omega(B_{R_0}(0))$,
\item[(ii)] for all $x \in Q = B_{2D_0}(x_0)$, $r \in (0,2R_0]$, if $\omega(V\cap B_r(x)) \ge
C \varepsilon \omega(B_r(x))$ then $B_r(x) \cap Q \subset W$, where $D_0 = \mathrm{diam}(\Omega)$.
\end{enumerate}
We outline the main steps in the proof following above conditional lemma. More precisely, we first show that $(i)$ holds.\\
 Without loss of generality, we may assume that $V \neq \emptyset$, then there exists $x_1 \in \Omega$ such that
\begin{equation}\label{eq:E1}
\mathbf{M}_{\alpha}(|f|^{\frac{p}{p-1}} + |\nabla g|^p)(x_1) \le \kappa \lambda.
\end{equation}
By Lemma~\ref{lem:boundM}, the boundedness property of maximal function $\mathbf{M}$ from $L^1(\Omega)$ to $L^{1,\infty}(\Omega)$ gives
\begin{equation}\label{eq:E4}
|V| \le \left|\left\{\mathbf{M}\mathbf{M}_{\alpha}(|\nabla u|^{p}))>{\sigma}\lambda \right\} \cap \Omega\right|  \le \frac{1}{{\sigma}\lambda} \int_{\Omega}\mathbf{M}_{\alpha}\left(|\nabla u|^p\right) dx.
\end{equation}
Applying the boundedness property of fractional maximal function $\mathbf{M}_{\alpha}$ in Lemma~\ref{lem:Malpha}, there holds
\begin{align*}
 \int_{\Omega} \mathbf{M}_\alpha(|\nabla u|^p) dx&= \int_{0}^{\infty} {\omega}\left(\left\{x\in \Omega: \ \mathbf{M}_\alpha(|\nabla u|^p)(x)>\lambda\right\}\right) d\lambda\\
 & \leq  CD_0^n\lambda_0+\int_{\lambda_0}^{\infty} {\omega}\left(\left\{x \in \Omega: \ \mathbf{M}_\alpha(|\nabla u|^p)(x)>\lambda\right\}\right) d\lambda\\
 &\leq C D_0^n\lambda_0+ C \left(\int_{\Omega}|\nabla u|^pdx\right)^{\frac{n}{n-\alpha}}\int_{\lambda_0}^{\infty}\lambda^{-\frac{n}{n-\alpha}}d\lambda\\
 &= CD_0^n\lambda_0+ C \left(\int_{\Omega}|\nabla u|^pdx\right)^{\frac{n}{n-\alpha}}\lambda_0^{-\frac{\alpha}{n-\alpha}},
\end{align*}
for any $\lambda_0>0$. In this formula, let us choose
\begin{equation*}
\lambda_0=D_0^{-n+\alpha}\int_{\Omega}|\nabla u|^pdx,
\end{equation*}
which follows that
\begin{align}\label{eq:E3}
\int_{\Omega} \mathbf{M}_{\alpha}\left(|\nabla u|^p\right) dx \le C_1 D_0^{\alpha} \int_{\Omega} |\nabla u|^p dx.
\end{align}
For the sake of readability, in some inequalities follow, the constants $C_i$ appearing might vary and must be indicated precisely. As such, this makes sense when we choose the value $\varepsilon>0$ in the statement of theorem at the end of proof depends only on a specific final constant.\\
Plugging the validity of \eqref{eq:E3} to~\eqref{eq:E4} and~\eqref{eq:prop1} from Proposition~\ref{prop1} to infer that
\begin{align}\nonumber
|V|  &\le \frac{C_2D_0^{\alpha}}{{\sigma}\lambda} \int_{\Omega} |\nabla g|^p + |f|^{\frac{p}{p-1}} dx  \\ \nonumber
 & \le  \frac{C_2D_0^{\alpha}}{{\sigma}\lambda} \int_{B_{D_0}(x_1)} |\nabla g|^p + |f|^{\frac{p}{p-1}} dx  \\ \label{eq:E6}
 & \le \frac{C_2D_0^{n}}{{\sigma}\lambda} \mathbf{M}_{\alpha}(|\nabla g|^p + |f|^{\frac{p}{p-1}})(x_1). 
\end{align}
Thanks to~\eqref{eq:E1}, it deduces from~\eqref{eq:E6} that
\begin{align}\label{eq:E7}
|V|  \le \frac{C_2D_0^{n}}{{\sigma}\lambda} \kappa \lambda \le \frac{C_3 \kappa}{{\sigma}} |B_{R_0}(0)|.
\end{align}
In view of the definition of Muckenhoupt weight $\mathcal{A}_{\infty}$, we get by~\eqref{eq:E7} that
\begin{align*}
\omega(V) \le C_4 \left(\frac{|V|}{|B_{R_0}(0)|}\right)^{\nu} \omega\left(B_{R_0}(0)\right)  \le {C_5}{({\sigma})}^{-\nu} \kappa^\nu  \omega\left(B_{R_0}(0)\right) \le C_5 \varepsilon \omega\left(B_{R_0}(0)\right),
\end{align*}
where $\kappa$ is small enough satisfying $\kappa \le {\sigma} \varepsilon^{1/\nu}$, we then immediately obtain $(i)$.\\

Let $x \in \Omega$, $r \in (0,2R_0]$ and $\lambda>0$, the remainder will be dedicated to the proof of $(ii)$, and the proof performed via a contradiction. Let us assume that $V \cap B_r(x) \neq \emptyset$ and $B_r(x) \cap \Omega \cap W^c \neq \emptyset$, i.e., there exist $x_2, x_3 \in B_r(x) \cap \Omega$ such that 
\begin{equation}\label{eq:M}
\mathbf{M}\mathbf{M}_{\alpha}(|\nabla u|^p)(x_2) \le \lambda \quad \mbox{ and } \quad \mathbf{M}_{\alpha}(|f|^{\frac{p}{p-1}}+|\nabla g|^p)(x_3) \le \kappa \lambda.
\end{equation}
We will show that 
\begin{equation}\label{eq:F1}
\omega(V \cap B_r(x)) < \varepsilon \omega(B_r(x)),
\end{equation}
which is a contradiction by Lemma~\ref{lem:wReig}. Indeed, for any $y \in B_r(x)$, it is easy to see that 
$$ B_{\rho}(y) \subset  B_{\rho+r}(x) \subset B_{\rho+2r}(x_2) \subset B_{3\rho}(x_2), \ \mbox{ for all } \ \rho \ge r,$$ which follows from \eqref{eq:M} that
\begin{align}\nonumber
\sup_{\rho\ge r}\fint_{B_{\rho}(y)}\mathbf{M}_{\alpha}(|\nabla u|^p) dx & \le 3^n  \sup_{\rho\ge r}\fint_{B_{3\rho}(x_2)}\mathbf{M}_{\alpha}(|\nabla u|^p)dx \\ \nonumber
&\le 3^n \mathbf{M}\mathbf{M}_{\alpha}(|\nabla u|^p)(x_2) \\ \label{eq:T1}
&\le 3^n \lambda.
\end{align}
Similarly, for any $y \in B_r(x)$, for all $0<\rho<r$ and  $z \in B_\rho(y)$, since $B_{\varrho}(z) \subset B_{4r}(x_2)$ for any $\varrho \ge r$, we also obtain that
\begin{align}\label{eq:T2}
\sup_{0<\rho<r} \fint_{B_\rho(y)} \left( \sup_{\varrho\ge r} \varrho^{\alpha} \fint_{B_{\varrho}(z)} |\nabla u|^p \right) dz \le 4^{n-\alpha} \mathbf{M}\mathbf{M}_{\alpha} (|\nabla u|^p)(x_2)  \le 4^{n-\alpha}\lambda.
\end{align}
Moreover, by the definitions of $\mathbf{M}$ and $\mathbf{M}_\alpha$, we can conclude from~\eqref{eq:T1} and~\eqref{eq:T2} that
\begin{align}\nonumber
\mathbf{M}\mathbf{M}_{\alpha}(|\nabla u|^p)(y) & \le \max\left\{\sup_{0<\rho<r} \fint_{B_\rho(y)} \left( \sup_{0<\varrho< r} \varrho^{\alpha} \fint_{B_{\varrho}(z)} |\nabla u|^p \right) dz; \right. \\ \nonumber
& \left. \sup_{0<\rho<r} \fint_{B_\rho(y)} \left( \sup_{\varrho\ge r} \varrho^{\alpha} \fint_{B_{\varrho}(z)} |\nabla u|^p \right) dz;  \quad \sup_{\rho\ge r}\fint_{B_{\rho}(y)}\mathbf{M}_{\alpha}(|\nabla u|^p) dx \right\} \\ \label{eq:T3}
& \le \max\left\{\sup_{0<\rho<r} \fint_{B_\rho(y)} \left( \sup_{0<\varrho< r} \varrho^{\alpha} \fint_{B_{\varrho}(z)} |\nabla u|^p \right) dz; \quad 4^n \lambda \right\},
\end{align}
for all $y \in B_r(x)$. Thanks to~\eqref{eq:MrMr} in Lemma~\ref{lem:MrMr}, we obtain that from~\eqref{eq:T3}, it provides
\begin{align}\label{eq:T4}
\mathbf{M}\mathbf{M}_{\alpha}(|\nabla u|^p)(y)  \le \max\left\{\mathbf{M}^{2r}_{\alpha} (\chi_{B_{2r}(x)}|\nabla u|^p) (y); \quad 4^n \lambda \right\}.
\end{align}
Hence if $\sigma$ is chosen satisfying ${\sigma} > 4^{n}$, then for any $\lambda>0$, by~\eqref{eq:T4} there holds
\begin{align}\label{eq:F2}
|V \cap B_r(x)| \le \left|\left\{\mathbf{M}^{2r}_{\alpha} (\chi_{B_{2r}(x)}|\nabla u|^p)  > {\sigma} \lambda\right\} \cap B_r(x) \cap \Omega\right|.
\end{align}
We remark that if $\overline{B}_{8r}(x) \subset \mathbb{R}^n\setminus \Omega$ then $V \cap B_r(x) = \emptyset$. So we need to consider two cases: $x$ is in the interior domain $B_{8r}(x) \subset \subset \Omega$ and $x$ is near the boundary $B_{8r}(x) \cap \partial \Omega \neq \emptyset$. And the proof in each case consists in matching the comparison estimates of Lemmas and Theorems in the interior domain and on the boundary.\\

Let us now consider the first case $B_{8r}(x) \subset \subset \Omega$. Let $v$ be the unique solution to~\eqref{eq:I2} with $B_{R}(x_0) = B_{4r}(x)$. Thanks to Theorem~\ref{theo:I} with the fact 
$$[A]^{4r}_{s_1} \le [A]^{R_0}_{s_0} \le \delta <1,$$ 
one has
\begin{multline}\label{est:Ia}
\|\nabla v\|^p_{L^{\infty}(B_{2r}(x))} \le C_6 \fint_{B_{8r}(x)} |\nabla u|^p dx + C_6 \fint_{B_{8r}(x)} |f|^{\frac{p}{p-1}} + |\nabla g|^p dx \\ + C_6 \left(\fint_{B_{8r}(x)} |\nabla g|^p dx \right)^{\frac{1}{p}} \left(\fint_{B_{8r}(x)} |\nabla u|^p dx \right)^{\frac{p-1}{p}}
\end{multline}
and
\begin{multline}\label{est:Ib}
\fint_{B_{4r}(x)} |\nabla u - \nabla v|^pdx \le C_6 \delta^p \fint_{B_{8r}(x)} |\nabla u|^p dx + C_6 \fint_{B_{8r}(x)} |f|^{\frac{p}{p-1}} + |\nabla g|^p dx  \\ 
+ C_6 \left(\fint_{B_{8r}(x)} |\nabla g|^p dx \right)^{\frac{1}{p}} \left(\fint_{B_{8r}(x)} |\nabla u|^p dx \right)^{\frac{p-1}{p}}.
\end{multline}
Note that $B_{8r}(x) \subset B_{9r}(x_2) \cap B_{9r}(x_3)$, we obtain from~\eqref{eq:M} and~\eqref{est:Ia} the following estimates
\begin{align} \nonumber
\mathbf{M}^{2r}_{\alpha} (\chi_{B_{2r}(x)}|\nabla v|^p)(y) & \le (2r)^{\alpha}\|\nabla v\|^p_{L^{\infty}(B_{2r}(x))} \\ \nonumber
 & \le  C_7 \left( r^{\alpha} \fint_{B_{9r}(x_2)} |\nabla u|^p dx +   r^{\alpha} \fint_{B_{9r}(x_3)} |f|^{\frac{p}{p-1}} + |\nabla g|^p dx \right)  \\ \nonumber
& \qquad \quad +  C_7   \left(r^{\alpha}\fint_{B_{9r}(x_3)} |\nabla g|^p dx \right)^{\frac{1}{p}} \left(r^{\alpha}\fint_{B_{9r}(x_2)} |\nabla u|^p dx \right)^{\frac{p-1}{p}}\\ \nonumber
& \le C_8 \left(\mathbf{M}\mathbf{M}_{\alpha}(|\nabla u|^p)(x_2) + \mathbf{M}\mathbf{M}_{\alpha}(|f|^{\frac{p}{p-1}} + |\nabla g|^p)(x_3)\right)  \\ \nonumber
& \qquad \quad + C_8  \left(\mathbf{M}\mathbf{M}_{\alpha}(|\nabla g|^p)(x_3)\right)^{\frac{1}{p}}  \left(\mathbf{M}\mathbf{M}_{\alpha}(|\nabla u|^p)(x_2)\right)^{\frac{p-1}{p}}  \\ \nonumber
& \le C_8 \left(\lambda + \kappa \lambda + \kappa^{\frac{1}{p}} \lambda\right) \\ \label{eq:M2r1}
& \le C_9 \lambda,
\end{align}
and deduce from~\eqref{est:Ib} that
\begin{align} \nonumber
(4r)^{\alpha}\fint_{B_{4r}(x)} |\nabla u - \nabla v|^p dx & \le C_{10} \delta^p r^{\alpha} \fint_{B_{9r}(x_2)} |\nabla u|^p dx \\ \nonumber 
& \ + C_{10}  \left[ r^{\alpha}  \fint_{B_{9r}(x_3)} |f|^{\frac{p}{p-1}} + |\nabla g|^p dx  \right. \\  \nonumber
& \ + \left. \left( r^{\alpha} \fint_{B_{9r}(x_3)} |\nabla g|^p dx \right)^{\frac{1}{p}} \left( r^{\alpha} \fint_{B_{9r}(x_2)} |\nabla u|^p dx \right)^{\frac{p-1}{p}} \right]\\ \label{eq:M2r2}
& \le C_{10} \left(\delta^p + \kappa + \kappa^{\frac{1}{p}}\right)\lambda. 
\end{align}
It follows easily from~\eqref{eq:M2r1} that if ${\sigma} \ge \max\{4^{n},2^pC_9\}$, then
$$\left| \left\{ \mathbf{M}^{2r}_{\alpha}\left(\chi_{B_{2r}(x)}|\nabla v|^p\right)>2^{-p}{\sigma}\lambda\right\} \cap B_r(x)\right| = 0$$
which implies from~\eqref{eq:F2} that
$$\left|V \cap B_r(x)\right| \le \left| \left\{ \mathbf{M}^{2r}_{\alpha}\left(\chi_{B_{2r}(x)}|\nabla u-\nabla v|^p\right)>2^{-p}{\sigma}\lambda\right\} \cap B_r(x)\right|.$$
Using again the bounded property of the fractional maximal function $\mathbf{M}_{\alpha}$ in Lemma~\ref{lem:Malpha}, we obtain from the above estimate and~\eqref{eq:M2r2} that
\begin{align}\nonumber
\left|V \cap B_r(x)\right| & \le \frac{C_{11}}{(2^{-p}{\sigma}\lambda)^{\frac{n}{n-\alpha}}} \left( \int_{B_{2r}(x)} |\nabla u - \nabla v|^p dx\right)^{\frac{n}{n-\alpha}} \\ \nonumber
& \le \frac{C_{11}}{(2^{-p}{\sigma}\lambda)^{\frac{n}{n-\alpha}}} (4r)^{n} \left(  (4r)^{\alpha} \fint_{B_{4r}(x)} |\nabla u - \nabla v|^p dx\right)^{\frac{n}{n-\alpha}} \\  \label{eq:M2r3}
& \le C_{12}\left(\delta^p + \kappa + \kappa^{\frac{1}{p}}\right)^{\frac{n}{n-\alpha}} |B_r(x)|.
\end{align}
By the definition of the Muckenhoupt weight $\omega \in \mathcal{A}_{\infty}$, we deduce~\eqref{eq:F1} from~\eqref{eq:M2r3}, that means
\begin{align*}
\omega(V \cap B_r(x)) & \le C \left(\frac{|V \cap B_r(x)|}{|B_r(x)|}\right)^\nu \omega(B_r(x)) \\
& \le C_{13}\left(\delta^p + \kappa + \kappa^{\frac{1}{p}}\right)^{\frac{n\nu}{n-\alpha}} \omega(B_r(x)) \\
& < \varepsilon \omega(B_r(x)),
\end{align*}
where $\kappa$ and $\delta$ are small enough such that 
$$C_{13}\left(\delta^p + \kappa + \kappa^{\frac{1}{p}}\right)^{\frac{n\nu}{n-\alpha}} < \varepsilon.$$

Let us next consider the second case when $x$ is near the boundary $B_{8r}(x) \cap \partial \Omega \neq \emptyset$. Let $x_4 \in \partial \Omega$ such that $|x_4-x| = \mbox{dist}(x,\partial\Omega)$. We remark that
$$ B_{2r}(x) \subset B_{10r}(x_4) \subset B_{540r}(x_4) \subset B_{549r}(x) \subset B_{550r}(x_2) \cap B_{550r}(x_3).$$
Applying Theorem~\ref{theo:B} with $B_{6R}(x_0) = B_{540r}(x_4)$ and $\vartheta \in (0,1)$ with notice that $[A]_{s_2}^{R_0} \le \delta$, we obtain that
\begin{multline*}
\|\nabla U\|^p_{L^{\infty}(B_{10r}(x_4))} \le C_{14} \left( \fint_{B_{540r}(x_4)} |\nabla u|^p dx + \fint_{B_{540r}(x_4)} |f|^{\frac{p}{p-1}} +|\nabla g|^p dx \right) \\ + C_{14} \left(\fint_{B_{540r}(x_4)} |\nabla g|^p dx \right)^{\frac{1}{p}}\left(\fint_{B_{540r}(x_4)} |\nabla u|^p dx \right)^{\frac{p-1}{p}},
\end{multline*}
\begin{multline*}
\fint_{B_{10r}(x_4)} |\nabla u - \nabla U|^p dx \le C_{15} \left( \vartheta^p + \delta^p \right) \fint_{B_{540r}(x_4)} |\nabla u|^p dx  \\ +   C_{15} \fint_{B_{540r}(x_4)} |f|^{\frac{p}{p-1}} + |\nabla g|^p dx \\ +  C_{15}\left(\fint_{B_{540r}(x_4)} |\nabla g|^p dx \right)^{\frac{1}{p}}\left(\fint_{B_{540r}(x_4)} |\nabla u|^p dx \right)^{\frac{p-1}{p}}. 
\end{multline*}
As in the proof of the first case, thanks to~\eqref{eq:M}, it follows from the above estimates that
\begin{align*}
(10r)^{\alpha}\|\nabla U\|^p_{L^{\infty}(B_{10r}(x_4))} & \le C_{16} \left(r^{\alpha} \fint_{B_{550r}(x_2)} |\nabla u|^p dx  + r^{\alpha}\fint_{B_{550r}(x_3)} |f|^{\frac{p}{p-1}} + |\nabla g|^p dx \right)\\
 & \qquad + C_{16} \left(r^{\alpha}\fint_{B_{550r}(x_3)} |\nabla g|^p dx \right)^{\frac{1}{p}}\left(r^{\alpha}\fint_{B_{550r}(x_2)} |\nabla u|^p dx \right)^{\frac{p-1}{p}} \\
& \le C_{17}\left(1+\kappa+\kappa^{\frac{1}{p}}\right)\lambda \\
& \le C_{18}\lambda,
\end{align*}
and
\begin{align*}
(2r)^{\alpha}\fint_{B_{2r}(x)} |\nabla u - \nabla U|^p dx 
& \le C_{19} \left( \varepsilon^p + \delta^p \right) r^{\alpha}\fint_{B_{550r}(x_2)} |\nabla u|^p dx  \\
& \qquad +   C_{19} r^{\alpha}\fint_{B_{550r}(x_3)} |f|^{\frac{p}{p-1}} + |\nabla g|^p dx \\ 
 & \qquad + C_{19} \left(r^{\alpha}\fint_{B_{550r}(x_3)} |\nabla g|^p dx \right)^{\frac{1}{p}}\left(r^{\alpha}\fint_{B_{550r}(x_2)} |\nabla u|^p dx \right)^{\frac{p-1}{p}} \\
& \le C_{20} \left(\vartheta^p + \delta^p +  \kappa +\kappa^{\frac{1}{p}} \right)\lambda.
\end{align*}
Therefore, for ${\sigma} \ge \max\{4^{n},2^pC_9, 2^pC_{18}\}$, we may conclude that
\begin{align*}
\left|V \cap B_r(x)\right| & \le \left|\left\{ \mathbf{M}^{2r}_{\alpha}\left(\chi_{B_{2r}(x)}|\nabla u-\nabla U|^p\right)>2^{-p}{\sigma}\lambda\right\}  \cap B_r(x)\right| \\
& \le \frac{C_{21}}{\left(2^{-p}{\sigma}\lambda\right)^{\frac{n}{n-\alpha}}} \left(\int_{B_{2r}(x)} |\nabla u - \nabla U|^pdx\right)^{\frac{n}{n-\alpha}} \\
& \le \frac{C_{21}}{\left(2^{-p}{\sigma}\lambda\right)^{\frac{n}{n-\alpha}}} (2r)^{n} \left((2r)^{\alpha}\fint_{B_{2r}(x)} |\nabla u - \nabla U|^pdx\right)^{\frac{n}{n-\alpha}} \\
& \le C_{22}\left(\vartheta^p + \delta^p +  \kappa +\kappa^{\frac{1}{p}} \right)^{\frac{n}{n-\alpha}}|B_r(x)|.
\end{align*}
By the definition of Muckenhoupt weight $\omega \in \mathcal{A}_{\infty}$, this follows that
\begin{align*}
\omega(V \cap B_r(x)) & \le C \left(\frac{|V \cap B_r(x)|}{|B_r(x)|}\right)^\nu \omega(B_r(x)) \\
& \le C_{23}\left(\vartheta^p + \delta^p +  \kappa +\kappa^{\frac{1}{p}} \right)^{\frac{n\nu}{n-\alpha}} \omega(B_r(x)).
\end{align*}
To complete the proof, all we need is to choose $\kappa$, $\delta$, $\vartheta$ positive and small enough such that 
$$C_{23}\left(\vartheta^p + \delta^p +  \kappa +\kappa^{\frac{1}{p}} \right)^{\frac{n\nu}{n-\alpha}} < \varepsilon.$$
\end{proof}


\medskip
\begin{proof}[Proof of Theorem~\ref{theo:regularity}]
By Theorem~\ref{theo:goodlambda}, for any $\varepsilon>0$ and $\lambda>0$, there exist some positive constants $\delta$, $s_0$, ${\sigma}$ and $\kappa$ such that if $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain satisfying $[A]_{s_0}^{R_0} \le \delta$ for some $R_0>0$ then
\begin{align}\label{eq:main-1}
\omega\left(V \right) \leq C \varepsilon \omega\left(W \right),
\end{align}
where $W = \{\mathbf{M}\mathbf{M}_{\alpha}(|\nabla u|^{p})> \lambda\}\cap \Omega$ and 
\begin{align*}
V = \left\{\mathbf{M}\mathbf{M}_{\alpha}(|\nabla u|^{p}))>{\sigma}\lambda, \ \mathbf{M}_{\alpha}(|f|^{\frac{p}{p-1}}+|\nabla g|^p) \le \kappa\lambda \right\}\cap \Omega.
\end{align*}
We deduce from~\eqref{eq:main-1} that
\begin{align}\nonumber
{\omega}\left(\{{\mathbf{M}\mathbf{M}}_\alpha(|\nabla u|^p)>\sigma\lambda\}\right) &\le C\varepsilon {\omega}\left(\{{\mathbf{M}\mathbf{M}_{\alpha}}(|\nabla u|^p)>\lambda\}\cap\Omega \right)\\ \label{eq:main-2}
& \qquad \qquad + {\omega}\left(\{{\mathbf{M}}_\alpha(|f|^{\frac{p}{p-1}}+|\nabla u|^p)>\kappa\lambda\}\cap\Omega \right).
\end{align}
By the definition of the norm given in \eqref{eq:w-lorentz-a}, one has
\begin{align*}
\|{\mathbf{M}\mathbf{M}}_\alpha(|\nabla u|^p)\|^s_{L^{q,s}_{\omega}(\Omega)} = q \int_0^\infty{\lambda^s {\omega}\left(\{{\mathbf{M}\mathbf{M}}_\alpha(|\nabla u|^p)>\lambda\} \right)^{\frac{s}{q}}\frac{d\lambda}{\lambda}}.
\end{align*}
Changing the variable $\lambda$ to $\sigma\lambda$ within the integral, we get that
\begin{align} \label{eq:main-3}
\|{\mathbf{M}\mathbf{M}}_\alpha(|\nabla u|^p)\|^s_{L^{q,s}_{\omega}(\Omega)} = \sigma^{s}q\int_0^\infty{\lambda^s{\omega}\left(\{{\mathbf{M}\mathbf{M}}_\alpha(|\nabla u|^p)>\sigma\lambda\} \right)^{\frac{s}{q}}\frac{d\lambda}{\lambda}}.
\end{align}
Thanks to~\eqref{eq:main-2}, it follows from~\eqref{eq:main-3} that
\begin{align*}
\|{\mathbf{M}\mathbf{M}}_\alpha(|\nabla u|^p)\|^s_{L^{q,s}_{\omega}(\Omega)} &\le \sigma^{s}(2C\varepsilon)^{\frac{s}{q}}q\int_0^\infty{\lambda^s{\omega}\left(\{{\mathbf{M}\mathbf{M}}_\alpha(|\nabla u|^p)>\lambda\}\cap\Omega \right)^{\frac{s}{q}}\frac{d\lambda}{\lambda}}\\
&~~~+ \sigma^{s}2^{\frac{s}{q}} q \int_0^\infty{\lambda^s{\omega}\left(\{{\mathbf{M}}_\alpha(|f|^{\frac{p}{p-1}}+|\nabla \sigma|^p)>\kappa\lambda\}\cap\Omega \right)^{\frac{s}{q}}\frac{d\lambda}{\lambda}}.
\end{align*}
Making the change of variables again in the second integral on right-hand side, we obtain that
\begin{align*}
\|{\mathbf{M}\mathbf{M}}_\alpha(|\nabla u|^p)\|^s_{L^{q,s}_{\omega}(\Omega)} &\le \sigma^{s}(2C\varepsilon)^{\frac{s}{q}}\|{\mathbf{M}\mathbf{M}}_\alpha\left(|\nabla u|^p \right)\|^s_{L^{q,s}_{\omega}(\Omega)}\\ & \qquad \qquad + \sigma^{s}2^{\frac{s}{q}} \kappa^{-s} \|{\mathbf{M}}_\alpha(|f|^{\frac{p}{p-1}}+|\nabla \sigma|^p)\|^s_{L^{q,s}_{\omega}(\Omega)}.
\end{align*}
Finally, the proof is complete by taking $\sigma^{s}(2C\varepsilon)^{\frac{s}{q}} \le \frac{1}{2}$.
\end{proof}

To complete the last studies of this section, we also prove the theorem \ref{theo:I-alpha} which relates the point-wise estimates for the Riesz potential. We first state and prove the following simple lemma.
\begin{lemma}\label{lem:Ialpha}
Let $\beta \in (0,n)$. If there exists a constant $C$ such that
\begin{align}\label{cond:w}
\int_{\mathbb{R}^n} f(x) \omega(x) dx \le C \int_{\mathbb{R}^n}g(x) \omega(x) dx, 
\end{align}
for all $\omega \in \mathcal{A}_1$, then 
\begin{align}\label{eq:wA1}
\mathbf{I}_{\beta}[f](z) \le C \mathbf{I}_{\beta}[g](z),
\end{align}
for almost everywhere $z \in \mathbb{R}^n$.
\end{lemma}
\begin{proof}
Firstly, let us consider $\omega_0(x) = |x|^{1-n}$, $ x \in \mathbb{R}^n$. It is easily seen that $\omega_0 \in \mathcal{A}_1$. In other words, there exists a constant $C_0>0$ such that $$\mathbf{M}(\omega_0)(x) \le C_0 \omega_0 (x).$$ By applying Fubini's Theorem, it concludes that for all $h \in L^1_{\mathrm{loc}}(\mathbb{R}^n;\mathbb{R}^+)$, there holds
\begin{align*}
\mathbf{M}(\mathbf{I}_{\beta}h)(x) \le C_0 \mathbf{I}_{\beta}h (x), \quad \forall x \in \mathbb{R}^n,
\end{align*}
which implies that $\mathbf{I}_{\beta}h \in \mathcal{A}_1$. For any $z \in \mathbb{R}^n$ and $\varepsilon>0$ small enough, we may choose $\omega = \mathbf{I}_{\beta}[\chi_{B_{\varepsilon}(z)}] \in \mathcal{A}_1$ as in~\eqref{cond:w} shows that
\begin{align*}
\int_{\mathbb{R}^n} f(x) \int_{\mathbb{R}^n}  \frac{\chi_{B_{\varepsilon}(z)}(y)}{|y-x|^{n-\beta}}  dy dx \le C \int_{\mathbb{R}^n} g(x)  \int_{\mathbb{R}^n}  \frac{\chi_{B_{\varepsilon}(z)}(y)}{|y-x|^{n-\beta}}  dy dx.
\end{align*} 
By Fubini's theorem, it gives
\begin{align*}
\int_{\mathbb{R}^n} \chi_{B_{\varepsilon}(z)}(y)   \int_{\mathbb{R}^n} \frac{f(x)}{|y-x|^{n-\beta}}  dx dy \le C \int_{\mathbb{R}^n} \chi_{B_{\varepsilon}(z)}(y)   \int_{\mathbb{R}^n}\frac{g(x)}{|y-x|^{n-\beta}}  dx dy.
\end{align*}
And it follows that
\begin{align}\label{est:Ifg}
\fint_{B_{\varepsilon}(z)} \mathbf{I}_{\beta}[f](y)dy \le C \fint_{B_{\varepsilon}(z)} \mathbf{I}_{\beta}[g](y)dy.
\end{align}
The assertion~\eqref{eq:wA1} of lemma completes by passing $\varepsilon$ to $0$ in \eqref{est:Ifg}.
\end{proof}
\medskip

\begin{proof}[Proof of Theorem~\ref{theo:I-alpha}]
Applying Theorem~\ref{theo:regularity}, we obtain at once that for any $\alpha \in [0,n)$, $0<q<\infty$ and $\omega \in \mathcal{A}_\infty$, there exist some positive constants $\delta$, $s_0$ such that if $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain satisfying $[A]_{s_0}^{R_0} \le \delta$ for some $R_0>0$, then
\begin{align*}
\|\mathbf{M}_{\alpha}(|\nabla u|^p)\|_{L^{q}_{\omega}(\Omega)} \le C \|\mathbf{M}_{\alpha}(|f|^{\frac{p}{p-1}} + |\nabla g|^p)\|_{L^{q}_{\omega}(\Omega)}, 
\end{align*}
which is equivalent to 
\begin{align*}
\left(\int_{\Omega} |\mathbf{M}_{\alpha}(|\nabla u|^p)|^q \omega(x) dx\right)^{\frac{1}{q}} \le C \left(\int_{\Omega} |\mathbf{M}_{\alpha}(|f|^{\frac{p}{p-1}} + |\nabla g|^p)|^q \omega(x) dx\right)^{\frac{1}{q}}. 
\end{align*}
Thanks to Lemma~\ref{lem:wReig}, it yields the following point-wise inequality:
\begin{align*}
\mathbf{I}_{\beta}\left(\chi_{\Omega}|\mathbf{M}_{\alpha}(|\nabla u|^p)|^q \right)(x) \le C \mathbf{I}_{\beta}\left( \chi_{\Omega}|\mathbf{M}_{\alpha}(|f|^{\frac{p}{p-1}} + |\nabla g|^p)|^q \right)(x), 
\end{align*}
for almost everywhere $x \in \mathbb{R}^n$. The proof is then complete.
\end{proof}

\section{Application}\label{sec:app}
In this section, we apply the point-wise estimate in Theorem~\ref{theo:I-alpha} to study the solvability of the generalized equation~\eqref{eq:existI}:
\begin{align*}
\begin{cases} -\mathrm{div}(A(x,\nabla u)) & = \ \mathbf{I}_{\beta}(|\nabla u|^p)^q + \mathrm{div}(f), \ \mbox{ in } \ \Omega, \\ \hspace{1.5cm} u & = \  g, \ \mbox{ on } \partial \Omega,\end{cases}
\end{align*}
where $\beta \in (0,n)$, $p>1$, $q>0$ and the same hypotheses of the domain $\Omega$ and the data $f$, $g$ in Theorem~\ref{theo:I-alpha}. In particular, we show that this equation admits at least one solution under an additional condition related to the Riesz capacity assumption on the data. Moreover, we also show that this type of  condition is necessary for the existence result. Our key point is based on some comparison estimates on Riesz and Wolff potentials which are firstly discussed in the next subsection.
\subsection{Comparisons on Riesz and Wolff potentials}
Let us first recall the following lemma which links a condition on the Wolff potential of a measure with a Riesz capacity assumption in the whole space.  We refer the reader to~\cite{Phuc2008} for the proof of this lemma.
\begin{proposition}\label{prop:Hung-Veron}
Let $1 < \beta_2 < \frac{n}{\beta_1}$, $s > \beta_2 -1$ and $\nu  \in \mathcal{M}^+(\mathbb{R}^n)$. Two following statements are equivalent:
\begin{itemize}
\item[(i)] The inequality
\begin{align*}
\nu (K) \le c \ \mathrm{Cap}_{\mathbf{I}_{\beta_1\beta_2},\frac{s}{s-\beta_2+1}}(K)
\end{align*}
holds for any compact set $K \subset \mathbb{R}^n$, for a constant $c$.
\item[(ii)] The inequality
\begin{align*}
\int_{\mathbb{R}^n} \left(\mathbf{W}_{\beta_1,\beta_2}(\chi_{B_r(x)}\nu )(y)\right)^s dy \le c \nu (B_r(x))
\end{align*}
holds for any ball $B_r(x) \subset \mathbb{R}^n$.
\end{itemize}
\end{proposition}

The next lemma is directly a consequence of~\cite[Lemma 2.1]{BNV2018}. The detail proof can be also found in~\cite{BNV2018}.
\begin{lemma}\label{lem:Hung-Veron}
Let $k, m \in \mathbb{R}^+$ and $s \in \mathbb{R}$. Assume that $H: \mathbb{R}^+ \to \mathbb{R}^+$ is a non-decreasing function. There exists a positive constant $C=C(k,m,s)$ such that
\begin{align*}
\int_0^{\infty} \varrho^{k} \left(\int_{\varrho}^{\infty} \frac{H(r)}{r^{s}} \frac{dr}{r}\right)^{m} \frac{d\varrho}{\varrho} \le C \int_0^{\infty} \varrho^{k} \left( \frac{H(\varrho)}{\varrho^{s}}\right)^{m} \frac{d\varrho}{\varrho}.
\end{align*}
\end{lemma}


Using Lemma~\ref{lem:Hung-Veron} and Proposition~\ref{prop:Hung-Veron}, we can perform some comparison estimates between Riesz and Wolff potentials.

\begin{lemma}\label{lem:2I} Let $\beta_1, \beta_2 \in (0,n)$. 
If $q>0$ and $\nu  \in \mathcal{M}^+(\mathbb{R}^n)$ then there exists $C_1$ depending on $\beta_1, \beta_2, n, q$ such that
\begin{align}\label{eq:lem-2Ia}
\mathbf{I}_{{\beta_2}}\left(\mathbf{I}_{{\beta_1}}(\nu )^q\right)(x) \ge C_1 \mathbf{W}_{\frac{q{\beta_1} + {\beta_2}}{q+1}, \frac{1}{q}+1}(\nu )(x), \quad \mbox{ in } \mathbb{R}^n.
\end{align}
If $0<q<\frac{n}{n-{\beta_1}}$ and $\nu  \in \mathcal{M}^+(\mathbb{R}^n)$ then there exists $C_2$ depending on $\beta_1, \beta_2, n, q$ such that
\begin{align}\label{eq:lem-2I}
\mathbf{I}_{{\beta_2}}\left(\mathbf{I}_{{\beta_1}}(\nu )^q\right)(x) \le C_2 \mathbf{W}_{\frac{q{\beta_1} + {\beta_2}}{q+1}, \frac{1}{q}+1}(\nu )(x), \quad \mbox{ in } \mathbb{R}^n.
\end{align}
\end{lemma}
\begin{proof} We remark that the Riesz potential $\mathbf{I}_{{\beta_1}}$ defined by~\eqref{def:Riesz} can be also rewritten as the following form
\begin{align*}
\mathbf{I}_{{\beta_1}}(\nu )(x) = \int_0^{\infty} \frac{\nu (B_r(x))}{r^{n-{\beta_1}}} \frac{dr}{r}.
\end{align*}
For every $x \in \mathbb{R}^n$, with notice that $B_{\varrho}(x) \subset B_{2\varrho}(y)$ for all $y \in B_{\varrho}(x)$, one gets that
\begin{align*}
\mathbf{I}_{{\beta_2}}\left(\mathbf{I}_{{\beta_1}}(\nu )^q\right)(x) & = \int_0^{\infty} \frac{1}{\varrho^{n-{\beta_2}}} \int_{B_{\varrho}(x)} \left( \int_0^{\infty} \frac{\nu (B_r(y))}{r^{n-{\beta_1}}} \frac{dr}{r} \right)^qdy \frac{d\varrho}{\varrho} \\
& \ge C \int_0^{\infty} \frac{1}{\varrho^{n-{\beta_2}}} \int_{B_{\varrho}(x)} \left( \int_0^{\infty} \frac{\nu (B_{2\varrho}(y))}{\varrho^{n-{\beta_1}}} dy \right)^q \frac{d\varrho}{\varrho} \\
& \ge C_1 \int_0^{\infty} \varrho^{{\beta_2}}  \left(  \frac{\nu (B_{\varrho}(x))}{\varrho^{n-{\beta_1}}} \right)^q \frac{d\varrho}{\varrho} \\
& = C_1 \int_0^{\infty}  \left(  \frac{\nu (B_{\varrho}(x))}{\varrho^{n-{\beta_1} -\frac{{\beta_2}}{q}}} \right)^q \frac{d\varrho}{\varrho} \\
& = C_1 \mathbf{W}_{\frac{q{\beta_1 + {\beta_2}}}{q+1},\frac{1}{q}+1}(\nu )(x),
\end{align*}
which is exactly~\eqref{eq:lem-2Ia}. On the other hand, for $0<q<\frac{n}{n-{\beta_1}}$, we recall the following estimate on Riesz potential
\begin{align*}
\|\mathbf{I}_{{\beta_1}}(\tilde{\nu })\|_{L^{\frac{n}{n-{\beta_1}},\infty}} \le C \tilde{\nu }({\mathbb{R}^n}), \quad \forall \tilde{\nu } \in \mathcal{M}^+_b(\mathbb{R}^n),
\end{align*}
which guarantees that 
\begin{align*}
\int_{B_r(x)}\mathbf{I}_{{\beta_1}}(\tilde{\nu })^q dy \le C r^n \left(\frac{\tilde{\nu }(\mathbb{R}^n)}{r^{n-{\beta_1}}}\right)^q, \quad \forall x \in \mathbb{R}^n, \ \forall r>0.
\end{align*}
Applying this inequality for $\tilde{\nu } = \chi_{B_{2r}(x)} \nu $, one has
\begin{align*}
\int_{B_r(x)}\mathbf{I}_{{\beta_1}}(\chi_{B_{2r}(x)} \nu )^q dy \le C r^n \left(\frac{\nu ({B_{2r}(x)})}{r^{n-{\beta_1}}}\right)^q, \quad \forall x \in \mathbb{R}^n, \ \forall r>0.
\end{align*}
Basing on this fact and notice that for all $r \ge \varrho>0$, since $B_r(y) \subset B_{2r}(x)$ for any $y \in B_{\varrho}(x)$, we may estimate as below
\begin{align*}
\mathbf{I}_{{\beta_2}}\left(\mathbf{I}_{{\beta_1}}(\nu )^q\right)(x) & \le C\int_0^{\infty} \frac{1}{\varrho^{n-{\beta_2}}} \int_{B_{\varrho}(x)} \left( \int_{\varrho}^{\infty} \frac{\nu (B_r(y))}{r^{n-{\beta_1}}} \frac{dr}{r} \right)^qdy \frac{d\varrho}{\varrho} \\
& \le C \int_0^{\infty} \frac{1}{\varrho^{n-{\beta_2}}} \int_{B_{\varrho}(x)} \left( \int_{\varrho}^{\infty} \frac{\nu (B_{2r}(x))}{r^{n-{\beta_1}}} \frac{dr}{r} \right)^q dy \frac{d\varrho}{\varrho} \\
& \le C \int_0^{\infty} {\varrho^{{\beta_2}}}  \left( \int_{\varrho}^{\infty} \frac{\nu (B_{2r}(x))}{r^{n-{\beta_1}}} \frac{dr}{r} \right)^q \frac{d\varrho}{\varrho}.
\end{align*}
Thanks to Lemma~\ref{lem:Hung-Veron}, we obtain from above inequality that
\begin{align*}
\mathbf{I}_{{\beta_2}}\left(\mathbf{I}_{{\beta_1}}(\nu )^q\right)(x) & \le C_2 \int_0^{\infty} {\varrho^{{\beta_2}}}  \left( \frac{\nu (B_{\varrho}(x))}{\varrho^{n-{\beta_1}}}\right)^q \frac{d\varrho}{\varrho},
\end{align*}
which leads to~\eqref{eq:lem-2I}. The proof is complete.
\end{proof}

\begin{lemma}\label{lem:cond-varphi} 
Let ${\beta_1}, \, {\beta_2}, \, \beta_3 \in (0,n)$, $0<q<\frac{n}{n-{\beta_1}}$, $qs >1$ and $\nu  \in \mathcal{M}^+(\mathbb{R}^n)$. Assume that the following inequality
\begin{align}\label{eq:cond-varphi-1}
\nu (K) \le \theta \ \mathrm{Cap}_{\mathbf{I}_{\beta_1 + \frac{\beta_2}{q}},\frac{qs}{qs-1}}(K),
\end{align}
holds for any compact set $K \subset \mathbb{R}^n$ and for a constant $\theta$. There holds
\begin{align}\label{re-lem:2I}
\mathbf{I}_{{\beta_3}}\left(\mathbf{I}_{{\beta_2}}\left(\mathbf{I}_{{\beta_1}}(\nu )^q\right)\right)^s(x) \le C \theta \ \mathbf{I}_{{\beta_3}}(\nu )(x), \quad \mbox{ in } \ \mathbb{R}^n.
\end{align}
\end{lemma}
\begin{proof}
By Proposition~\ref{prop:Hung-Veron}, the fact that the inequality~\eqref{eq:cond-varphi-1} holds for any compact set $K \subset \mathbb{R}^n$ is equivalent to
\begin{align}\label{eq:lem-2I-a}
\int_{\mathbb{R}^n} \left(\mathbf{W}_{\frac{q\beta_1 +\beta_2}{q+1},\frac{1}{q}+1}(\chi_{B_r(x)}\nu )(y)\right)^s dy \le \theta \nu (B_r(x))
\end{align}
holds for any ball $B_r(x) \subset \mathbb{R}^n$.
Applying Lemma~\ref{lem:2I}, there holds
\begin{align}\label{eq:lem-2I-b}
\mathbf{I}_{{\beta_2}}\left(\mathbf{I}_{{\beta_1}}(\nu )^q\right)(x) \le C_1\mathbf{W}_{\frac{q\beta_1 +\beta_2}{q+1},\frac{1}{q}+1}(\nu )(x), \quad \mbox{ in } \ \mathbb{R}^n.
\end{align}
It follows from~\eqref{eq:lem-2I-a} and~\eqref{eq:lem-2I-b} that
\begin{align*}
\mathbf{I}_{{\beta_3}}\left(\mathbf{I}_{{\beta_2}}\left(\mathbf{I}_{{\beta_1}}(\nu )^q\right)\right)^s(x) & = \int_0^{\infty} \frac{1}{\varrho^{n-\beta_3}}  \left(\int_{B_{\varrho}(x)} \left( \mathbf{I}_{{\beta_2}}\left(\mathbf{I}_{{\beta_1}}(\nu )^q\right)(y) \right)^s dy \right) \frac{d\varrho}{\varrho} \\
& \le  \int_0^{\infty} \frac{1}{\varrho^{n-\beta_3}}  \left(\int_{B_{\varrho}(x)} \left( C_1\mathbf{W}_{\frac{q\beta_1 +\beta_2}{q+1},\frac{1}{q}+1}(\nu )(y) \right)^s dy \right) \frac{d\varrho}{\varrho}\\
& \le C_1^s \theta\int_0^{\infty} \frac{\nu (B_{\varrho}(x)) }{\varrho^{n-\beta_3}}  \frac{d\varrho}{\varrho},
\end{align*}
which leads to~\eqref{re-lem:2I}. The proof is complete.
\end{proof}


\subsection{Existence result}
We now prove Theorem~\ref{theo:existence_Ialpha} which presents a sufficient condition for the existence of a solution to equation~\eqref{eq:existI}. We start with the following lemma which can be obtained from Theorem~\ref{theo:I-alpha}.
\begin{lemma}\label{lem:app-1}
Let $p>1$,  $f \in L^{\frac{p}{p-1}}(\Omega; \mathbb{R}^n)$, $g \in W^{1,p}(\Omega;\mathbb{R})$ and $\mathbf{I}_1(\eta ) \in  L^{\frac{p}{p-1}}(\Omega; \mathbb{R})$.  Let $u$ be a weak solution to the following equation
\begin{align}\label{eq-lem:app-1}
\begin{cases} -\mathrm{div}(A(x,\nabla u)) & = \ \eta  + \mathrm{div}(f), \ \mbox{ in } \ \Omega, \\ \hspace{1.5cm} u & = \  g, \ \mbox{ on } \partial \Omega.\end{cases}
\end{align}
There exist some positive constants $\delta$, $s_0$ such that if $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain satisfying $[A]_{s_0}^{R_0} \le \delta$ for some $R_0>0$, then for $\beta \in (0,n)$
\begin{align*
\mathbf{I}_{\beta} (|\nabla u|^p)(x) \le C^* \left(\mathbf{I}_{\beta} (\mathbf{I}_1(\eta )^{\frac{p}{p-1}}) + \mathbf{I}_{\beta}(|f|^{\frac{p}{p-1}} + |\nabla g|^p)\right)(x),
\end{align*}
for almost everywhere $x \in \mathbb{R}^n$ and for a constant $C^*>0$.
\end{lemma}
\begin{proof}
Let $B_R:=B_R(0) \supset \Omega$ and $u_0$ be the unique solution to the following equation
\begin{align*}
\begin{cases} - \Delta u_0 & = \ \eta  , \ \mbox{ in } \ B_R, \\ \hspace{0.5cm} u_0 & = \  0, \ \mbox{ on } \partial B_R.\end{cases}
\end{align*}
It is well known that the fundamental solution to this equation satisfies 
\begin{align}\label{est:fund}
|\nabla u_0(x)| \le |\nabla_x G(\eta )(x)| \le C_1\mathbf{I}_1(\eta )(x),
\end{align}
where $G$ denotes the Green kernel. Applying Theorem~\ref{theo:I-alpha} to equation~\eqref{eq-lem:app-1}, we can find some positive constants $\delta$, $s_0$ such that if $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain satisfying $[A]_{s_0}^{R_0} \le \delta$ for some $R_0>0$, then for $\beta \in (0,n)$
\begin{align*}
\mathbf{I}_{\beta}(|\nabla u|^p)  \le C_2 \left(\mathbf{I}_{\beta}(|\nabla u_0|^{\frac{p}{p-1}}) + \mathbf{I}_{\beta}(|f|^{\frac{p}{p-1}} + |\nabla g|^p) \right).
\end{align*}
It deduces from~\eqref{est:fund} to
\begin{align*}
\mathbf{I}_{\beta}(|\nabla u|^p)  \le  C^* \left(\mathbf{I}_{\beta} (\mathbf{I}_1(\eta )^{\frac{p}{p-1}}) + \mathbf{I}_{\beta}(|f|^{\frac{p}{p-1}} + |\nabla g|^p)\right),
\end{align*}
which finishes the proof.
\end{proof}

\medskip

\begin{proof}[Proof of Theorem~\ref{theo:existence_Ialpha}]
Let us first introduce a set $\mathcal{S}$ defined by
\begin{align}\label{S-def}
\mathcal{S} = \left\{v \in W^{1,p}(\Omega): \ \ \mathbf{I}_{\alpha}(|\nabla v|^p) \le \Lambda \mathbf{I}_{\alpha} (|\mathcal{F}|^p) \ \mbox{ a.e. in } \ \mathbb{R}^n\right\},
\end{align}
where the positive constant $\Lambda$ will be determined later. For every $v \in \mathcal{S}$, we define $T(v) := u$ as the unique solution to the following equation
\begin{align}\label{eq:app-2}
\begin{cases} -\mathrm{div}(A(x,\nabla u)) & = \ \mathbf{I}_{\beta}(|\nabla v|^p)^q + \mathrm{div}(f), \ \mbox{ in } \ \Omega, \\ \hspace{1.5cm} u & = \  g, \ \mbox{ on } \partial \Omega.\end{cases}
\end{align}
We next to show that one can find $\Lambda>0$ such that the mapping  $T:  \mathcal{S} \to \mathcal{S}$, $v \mapsto T(v) = u$ defined by~\eqref{eq:app-2} is well-defined.  In other words, we need to prove that $u = T(v) \in \mathcal{S}$ for all $v \in \mathcal{S}$. Indeed, thanks to Lemma~\ref{lem:app-1}, one obtains the following estimate
\begin{align}\label{est:app-1}
\mathbf{I}_{\alpha}(|\nabla u|^p) \le C^* \left(\mathbf{I}_{\alpha}(\mathbf{I}_1(\eta)^{\frac{p}{p-1}}) + \mathbf{I}_{\alpha}(|\mathcal{F}|^p)  \right), \quad \mbox{ a.e. in } \ \mathbb{R}^n,
\end{align}
where $\eta = (\mathbf{I}_{\beta}(|\nabla v|^p))^q \le \Lambda^q \mathbf{I}_{\beta}(|\mathcal{F}|^p)^q$ for $v \in \mathcal{S}$. It follows from~\eqref{est:app-1} that
\begin{align}\label{est:app-2}
\mathbf{I}_{\alpha}(|\nabla u|^p) \le C^* \left(\Lambda^{\frac{pq}{p-1}}\mathbf{I}_{\alpha}(\mathbf{I}_1(\mathbf{I}_{\beta}(|\mathcal{F}|^p)^q)^{\frac{p}{p-1}}) + \mathbf{I}_{\alpha}(|\mathcal{F}|^p)  \right), \quad \mbox{ a.e. in } \ \mathbb{R}^n.
\end{align}
Applying Lemma~\ref{lem:cond-varphi} under condition~\eqref{eq:ReiszCapa_cond} with $d\nu = |\mathcal{F}|^pdx$, $\beta_1 = \beta$, $\beta_2 = 1$, $\beta_3 = \alpha$ and $s = \frac{p}{p-1}$, there holds
\begin{align*}
\mathbf{I}_{\alpha}\left(\mathbf{I}_{1}\left(\mathbf{I}_{\beta}(|\mathcal{F}|^p)^q\right)^{\frac{p}{p-1}}\right)  \le C\varepsilon \mathbf{I}_{\alpha}(|\mathcal{F}|^p), \quad \mbox{ a.e. in } \ \mathbb{R}^n.
\end{align*}
Combining this inequality to~\eqref{est:app-2}, we can conclude that
\begin{align*}
\mathbf{I}_{\alpha}(|\nabla u|^p) \le C^* (C\varepsilon \Lambda^{\frac{pq}{p-1}} + 1) \mathbf{I}_{\alpha}(|\mathcal{F}|^p) \le \Lambda \mathbf{I}_{\alpha}(|\mathcal{F}|^p), \quad \mbox{ a.e. in } \ \mathbb{R}^n,
\end{align*}
which yields that $u \in \mathcal{S}$. We remark that in the last inequality, we may easily choose $\Lambda>0$ such that $C^* (C\varepsilon \Lambda^{\frac{pq}{p-1}} + 1) \le \Lambda$ for some $\varepsilon$ small enough. 

On the other hand, it is easy to check that the set $\mathcal{S}$ is convex, closed, the mapping $T$ is continuous and $T(S)$ is precompact under the strong topology of $W^{1,p}(\Omega)$. By the Schauder fixed point theorem,  the mapping $T$ admits at least one fixed point in $\mathcal{S}$. Finally, the estimate~\eqref{eq:theo-D} is obtained by the definition of $\mathcal{S}$ in~\eqref{S-def}. It completes the proof.
\end{proof} 
\medskip

We now give a proof of Theorem~\ref{theoE}. We emphasize that the equation~\eqref{eq:existI-b} considered in this theorem is simpler than~\eqref{eq:existI} just for simplicity of the computation.

\medskip

\begin{proof}[Proof of Theorem~\ref{theoE}]
It is well-known that if $\Omega$ is a $(\delta,R_0)$-Reifenberg flat domain satisfying $[A]_{s_0}^{R_0} \le \delta$ for some positive constants $\delta, R_0, s_0$ and $u$ is a renormalized solution to the following equation
\begin{align*}
\begin{cases} -\mathrm{div}(A(x,\nabla v)) & = \ \eta, \ \mbox{ in } \ \Omega, \\ \hspace{1.5cm} u & = \  0, \ \mbox{ on } \partial \Omega,\end{cases}
\end{align*}
then there exists a constant $C$ such that
\begin{align*}
\int_{\mathbb{R}^n} \mathbf{I}_1(\eta)^{\frac{p}{p-1}} d\omega \le C \int_{\mathbb{R}^n} |\nabla v|^p d\omega,
\end{align*}
for all $\omega \in \mathcal{A}_{\frac{p}{p-1}} \subset \mathcal{A}_1$.
Applying this fact for the solution $u$ to equation~\eqref{eq:existI-b}, we obtain that
\begin{align*}
\mathbf{I}_{\beta}(\mathbf{I}_1(\mathbf{I}_\beta(|\nabla u|^p)^q + \mu)^{\frac{p}{p-1}}) \le C \mathbf{I}_{\beta}(|\nabla u|^{p}) \quad \mbox{ a.e. in } \ \mathbb{R}^n,
\end{align*}
which follows that
$\mathbf{I}_{\beta}(\mathbf{I}_1(\nu)^{\frac{p}{p-1}})^q \le C \nu \ \mbox{ a.e. in } \ \mathbb{R}^n,$
for $\nu = \mathbf{I}_\beta(|\nabla u|^p)^q + \mu$. Applying~\eqref{eq:lem-2Ia} in Lemma~\ref{lem:2I}, one gets
\begin{align*}
\left(\mathbf{W}_{\frac{p+(p-1)\beta}{2p-1},\frac{2p-1}{p}} (\nu)\right)^q \le C \nu \quad \mbox{ a.e. in } \ \mathbb{R}^n.
\end{align*}
It implies that the following inequality holds for any ball $B_r(x) \subset \mathbb{R}^n$:
\begin{align*}
\int_{\mathbb{R}^n}\left(\mathbf{W}_{\frac{p+(p-1)\beta}{2p-1},\frac{2p-1}{p}} (\chi_{B_r(x)}\nu)(y)\right)^q dy \le C \nu(B_r(x)), 
\end{align*}
which yields that $\nu(K) \le C \mathrm{Cap}_{\mathbf{I}_{\beta + 1 - \frac{\beta}{q}},\frac{pq}{pq-p+1}}(K),$
for any compact subset $K \subset \mathbb{R}^n$. This implies to~\eqref{eq:ReiszCapa_cond-b} with notice that $\mu(K) \le \nu(K)$.
\end{proof}

