\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{amssymb}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2016 (2016), No. 33, pp. 1--14.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2016/33\hfil Subelliptic $p$-Laplace type systems]
{Morrey estimates for subelliptic $p$-Laplace type systems with VMO coefficients
 in Carnot groups}

\author[H. Yu, S. Zheng \hfil EJDE-2016/33\hfilneg]
{Haiyan Yu, Shenzhou Zheng}

\address{Haiyan Yu \newline
Department of Mathematics, 
Beijing Jiaotong University, 
Beijing 100044, China. \newline
College of Mathematics, 
Inner Mongolia University for Nationalities, 
Tongliao, 028043, China}
\email{12118381@bjtu.edu.cn}

\address{Shenzhou Zheng (corresponding author)\newline
Department of Mathematics, Beijing Jiaotong University,
Beijing 100044, China}
\email{shzhzheng@bjtu.edu.cn}

\thanks{Submitted August 16, 2015. Published January 21, 2016.}
\subjclass[2010]{35H20, 35B65, 35D30}
\keywords{Subelliptic p-Laplace; VMO coefficients; controllable growth;
\hfill\break\indent  Morrey regularity; Carnot group}

\begin{abstract}
 In this article, we study  estimates in Morrey spaces to the horizontal
 gradient of weak solutions for a class of quasilinear sub-elliptic systems
 of  $p$-Laplace type with VMO  coefficients under the controllable growth
 over Carnot group if $p$ is not too far from 2. We also show a local
 H\"older continuity  with an optimal exponent to the solutions.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{definition}[theorem]{Definition}
\allowdisplaybreaks

\section{Introduction}

Let $\mathcal{G}$ be a Carnot group of step $r\ge 2$, that is, a simply
connected Lie group with Lie algebra $\mathfrak{g}$ admits a decomposition
$\mathfrak{g} =\oplus_{i=1}^r V_j$ such that $[V_1, V_j ]= V_{j+1}$ for
 $1\le j\le r-1$  and $[V_1, V_{r}] = 0$. The homogeneous dimension of
$\mathcal{G}$ is defined as $Q =\sum_{i=1}^r im_i$, where $m_i = \dim V_i$
is the topological dimension with $m_1=m$.
For a family of vector fields $X=(X_{1},\,X_{2},\dots,\,X_{m})$ satisfying
the H\"ormander's finite rank condition rank Lie$[X_{1},\,X_{2},\dots,\,X_{m}]=r$,
 we assume that each component
$b_{ij}(x)$ ($i=1,\,2,\dots,m$; $j=1,2,\dots,n$) of vector-field
$X_{i}=\sum_{j=1}^nb_{ij}(x)\partial_j$ is a smooth function defined
in Carnot group $\mathcal{G}$ for $i=1,2,\dots,m$. Therefore,
$Xu=(X_{1}u,X_{2}u,\dots,X_{m}u)$ may be called the horizontal gradient
of $u$, which is understood as $X_i u=\langle X_i, \nabla u
\rangle=\Sigma_{j=1}^n b_{ij}(x)\partial_j u$ for $i=1,2,\dots,m$ if
 $u\in C^1(\mathcal{G})$, see \cite{DoM,Fol,Lu,CiGL,Zh3}.
To describe our assumptions and main results better, we first recall
some relevant notations and basic facts.

\begin{definition} \label{def1.1} \rm
  An absolutely continuous path $\gamma:[0,T]\to \mathcal{G} $ is said
$X$-subunit if
$$
\dot{\gamma}(t)=\sum_{i=1}^{m}c_{i}(t)X_{i}(\gamma(t))
$$
with
$\sum_{i=1}^{m}c_{i}(t)\leq 1$, for almost every $t\in [0,T]$.
\end{definition}

With $X$-subunit in hand, we can define  the Carnot-Caratheodory's metric
(the C-C distance) $d_X(x,y)$ as follows, see \cite{JoX,XuZ}.
$$
d_X(x,y)=\inf \{T>0: \exists
\gamma:[0,T]\to \mathcal{G}\ X \text{-subunit with }\
\gamma(0)=x,\,\gamma(T)=y\}.
$$

Note that these vector-fields $(X_1,\dots,X_m)$  are free up to the order
$r$, then there exists a positive constant $C>0$ satisfying the following
 relation between the C-C distance and the Euclidean metric, see \cite{XuZ,DiZ};
\begin{equation*}\label{metric-ineq}
C^{-1}|x-y|\le d_X(x,y)\le C|x-y|^{1/r}.
\end{equation*}

In this context, all balls centered at $x$ of radius $R$ with respect to
$d_X(x,y)$ are called metric balls and denoted still by $B_{R}(x)$.
The distance function $d_X(\cdot,\cdot)$ satisfies the local doubling property,
that is, for  $B_{2R}(x)\Subset  \mathcal{G} $ there exists a positive
 constant $R_0$ depending on vector fields $X$ and $ \mathcal{G}$ such that
for all $0<R \le R_0$ there holds
\begin{equation}\label{doubling}
|B_{2R}(x)|\le C_d|B_{R}(x)|,
\end{equation}
where the number $Q =\log_2 C_d$ is called the local homogeneous dimension of
$\mathcal{G}$ with respect
to the vector fields $X_1, X_2,\dots , X_m$. In fact,  $Q$ will play
a role of the dimension in the local analysis involving what we are
considering problems.


Let $\Omega$ be a bounded open subset in  Carnot group $\mathcal{G}$.
Let us recall the following horizontal Sobolev space with respect to
given a family of vector fields $X=(X_{1},X_{2},\dots,X_{m})$.
For any $1<p<\infty$ and $k\in \mathbb{N}$, we define
\begin{equation*} %\label{Sobolev}
HW^{k,p}(\Omega):=\{ u\in L^p(\Omega): X_{i_1}\dots X_{i_k} u\in L^p(\Omega)
\text{ for all } \{i_1,\dots,i_k\}\subset\{1,\dots,m\}\}
\end{equation*}
with the norm
$\|u\|_{HW^{k,p}(\Omega)}=\|u\|_{L^p(\Omega)}+\|X^ku\|_{L^p(\Omega)}$.
Furthermore, the closure of $C^{\infty}_0 (\Omega)$ in $HW^{k,p}(\Omega)$ is
denoted by $HW^{k,p}_0(\Omega)$.

In this article, we  consider the estimates in Morrey spaces to the horizontal
gradient of  weak solutions in $HW^{1,p}_0(\Omega)$ for the following
degenerate subelliptic systems of $p$-Laplace type.
\begin{equation}\label{eq1}
-\sum_{i=1}^{m}X_{i}\left(\langle A(x)X_{i}u,X_{i}u
 \rangle^{\frac{p-2}{2}}A(x)X_{i}u\right)=B(x,u,Xu),  \quad\text{a. e. }x\in \Omega,
\end{equation}
where $A(x)\in VMO\bigcap L^{\infty}(\Omega)$ and  $B(x,u,Xu)$ satisfies a
controllable growth. In order to more precisely impose structural assumptions
on $A(x), B(x,u,Xu)$ and state our main results, we need recall two useful
notations (see\cite{BrB,Gi,YuZ}).

\begin{definition}[BMO functions] \label{BMO} \rm
Let $\Omega(x,r)=\Omega\cap B_{r}(x)$. For any $0<s<+\infty$, we say
$u\in L^{1}_{\rm loc}(\Omega)$ belongs to $BMO(\Omega)$ if
\begin{equation*}
  M_{u}(s):=\sup_{x\in \Omega, 0<r<s}
\frac{1}{|\Omega(x,r)|}\int_{\Omega(x,r)}|u(y)-\bar{u}_{\Omega(x,r)}|dy<+\infty,
\end{equation*}
where
\begin{equation*}
\bar{u}_{\Omega(x,r)}=\hbox{--}\hskip-9pt\int_{\Omega(x,r)}u(y)dy
=\frac{1}{|\Omega(x,r)|}\int_{\Omega(x,r)}u(y)dy.
\end{equation*}
\end{definition}

\begin{definition}[VMO functions] \label{VMO} \rm
Let $M_{u}(s)$ be defined as above. We say $u\in BMO(\Omega)$ belongs to
$VMO(\Omega)$ if
\begin{equation*}
 \lim_{s\to 0}  M_{u}(s)=0,
\end{equation*}
where $M_{u}(s)$ is called  the $VMO$ modulus of $u$.
\end{definition}

\begin{definition}[Morrey space]  \rm
 Let $p\geq 1,\lambda > 0 $. For $u\in L^p_{\rm loc}(\Omega)$, if
\begin{equation}\label{Morrey}
\|u\|_{L^{p,\lambda}_{X}{(\Omega)}}:=\sup_{x_0\in\Omega,\, 0<r\leq d_0}
 \Big(\frac{r^{\lambda}}{|\Omega(x_{0},r)|}\int_{\Omega(x_{0},r)}|u|^p\,dx
\Big)^{1/p}<+\infty,
\end{equation}
then $u\in L^{p,\lambda}_{X}(\Omega)$, where
$d_0=\min\{ \operatorname{diam}(\Omega),R_D\}$, and the norm of $u$ is
$\|u\|_{L^{p,\lambda}_{X}{(\Omega)}}$.
\end{definition}


\begin{definition}[Campanato space] \label{campanato} \rm
Let $p\ge 1,\lambda>-p$. If $u\in L^p_{\rm loc}(\Omega)$ satisfies
\begin{equation}\label{Campanato}
|u|_{\mathcal {L}^{p,\lambda}_X(\Omega)}:=\sup_{x_0\in\Omega,\, 0<r\leq d_0}
\Big(\frac{r^{\lambda}}{|\Omega(x_0,r)|}
\int_{\Omega(x_0,r)}|u-u_{x_0,r}|^p\,dx\Big)^{1/p}<+\infty,
\end{equation}
then $u\in \mathcal{L}^{p,\lambda}_X{(\Omega)}$,  and has norm
$\|{u}\|_{{\mathcal L}^{p,\lambda}_X{(\Omega)}}
:=\|{u}\|_{L^p{(\Omega)}}+|u|_{\mathcal{L}^{p,\lambda}_X(\Omega)}$.
\end{definition}

On the basis of the above notation, we are now in a position to impose some
structure assumptions on $A(x)$ and $B(x,u,Xu)$ as follows.
\begin{itemize}
\item[(H1)]  (Uniform ellipticity)
For $A(x)=(a_{ij}^{\alpha\beta}(x))$, there exist $L$ and $\nu$,
$L\geq\nu>0$, such that for a.e. $x\in \Omega$ and for any
$\xi\in\mathbb{R}^{nN}$ we have
$$
\nu|\xi|^2\le a_{ij}^{\alpha\beta}(x)\xi_i^\alpha\xi_j^\beta\leq L|\xi|^2;
$$

\item[(H2)] $a_{ij}^{\alpha\beta}(x)\in L^{\infty}(\Omega)\bigcap VMO_X$;

\item[(H3)] (Controllable growth) The inhomogeneity $B(x,u,Xu)$ satisfies
$$
|B(x,u,Xu)|\leq \mu(|Xu|^{p(1-\frac{1}{\gamma})}+g(x)),
$$
where
\[
\gamma=\begin{cases}
\frac{pQ}{Q-p}, & 1<p<Q\\
\text{any } \gamma\ge p, &p\ge Q,
\end{cases}
\]
 $g(x)\in L^{q,\mu}(\Omega,\mathbb{R}^{nN})$ with $q>\frac {\gamma}{\gamma-1}$,
and $Q$ is the homogenous dimension.
\end{itemize}

We say that $u\in HW^{1,p}(\Omega,\mathbb{R}^{nN})$ is a weak solution
of \eqref{eq1}, if
\begin{equation}\label{eq2}
\int_{\Omega}\big\langle\langle A(x)Xu,Xu\rangle^{\frac{p-2}{2}}A(x)Xu,
X\varphi\big\rangle\,dx =\int_{\Omega}B(x,u,Xu)\varphi\,dx,
\end{equation}
for all $\varphi\in HW^{1,p}_{0}(\Omega)$.

Recently several studies on subelliptic PDEs arising from non-commuting vector
fields have been well developed based on the H\"ormander's fundamental work \cite{Ho};
see \cite{CaG,CaDG,DiZ,DoM, Fol, NaSW, RoS, XuZ, MZZ, WaL, Zh3, ZhF15}.
Many important results about the fundamental solution to subelliptic operators
and the Harmonic analysis theory on stratified nilpotent Lie groups have
been obtained by Folland \cite{Fol}, Rothschild-Stein \cite{RoS} and
 Nagel-Stein-Wainger \cite{NaSW}.
These results laid a solid foundation for further investigation of subelliptic
Partial Differiential Equations theory. Up to the 90s, the function theory and
harmonic analysis tools on Carnot groups, such as the Sobolev embedding
inequality of $X$-gradient and the isoperimetric inequality, become increasingly
mature,  cf. \cite{CaDG,Fol,GaN,NaSW,RoS}. In fact, such subelliptic problems
have received continuous attention due to their significant applications in
geometry and physics \cite{BeF, RoS}. In the case of Euclidean spaces
(i.e. $m = n$, $X_i =\frac{\partial}{\partial x_i}$), it was an important
observation by Uhlenbeck \cite{Uh} that there exists the interior
$C^{1,\alpha}$-regularity by using the classical De Giorgi-Moser-Nash
iteration to the homogeneous $p$-harmonic systems as a prototype.
However, it is not true for subelliptic systems of $p$-Laplacian with $p>1$.
Actually, Domokos-Manfredi in \cite{DoM,DoM1} and Domokos in \cite{DoM2,DoM3}
have derived $\Gamma^{1,\alpha}$ regularity for $p$-harmonic systems only
while $p$ is in a neighborhood of $2$ in Heisenberg group and in Carnot group,
respectively. Very recently, Zheng-Feng \cite{Zh3, ZhF15} also got the estimates
and an application of  the Green functions for  subelliptic A-harmonic operators,
and $\Gamma^{1,\alpha}$ regularity for weak solutions to subelliptic $p$-harmonic
systems under the subcritical growth with $p$ close to $2$, respectively.
Notice that Fazio-Fanciullo \cite{DiF} and Dong-Niu \cite{DoN} recently
established the estimates of the gradient in Morrey spaces to nonlinear
subelliptic systems for $p=2$.  Therefore, this is a natural thought what
happens if one consider a regularity of the gradient in Morrey spaces to
subelliptic A-harmonic systems. In this article, we are devoted to local
Morrey regularity of the horizontal gradient to a class of subelliptic A-harmonic
systems with VMO coefficients under the controllable growth.
More precisely, we have the following result.

\begin{theorem}\label{th1}
Let $u\in HW^{1,p}(\Omega,\mathbb{R}^{N})$ is a weak solution of \eqref{eq1}
with $p$ close to 2. Suppose $A(x)$ and $B(x,u,Xu)$ satisfy {\rm (H1)--(H3)}.
Then $Xu\in L_{X}^{p,\lambda}(\Omega,\mathbb{R}^{nN})$; moreover, there exists a
constant $C> 0$ such that for any $\Omega'\Subset \Omega$ we have
\begin{equation}
\|Xu\|_{L_X^{p,\lambda}(\Omega',\mathbb{R}^{nN})}
\leq C(\|Xu\|_{L^p(\Omega,\mathbb{R}^{nN})}+\|g\|^{\frac{1}{p-1}}_{L_X^{q,\mu}
(\Omega,\mathbb{R}^{nN})}),
\end{equation}
where
\[
\lambda=\begin{cases}
\frac{p}{p-1}\frac{\mu-q}{q}, &\frac {\gamma}{\gamma-1}<q<\mu,\\
\text{any  }\lambda\in (0,Q), & q\ge \mu.
\end{cases}
\]
\end{theorem}

Here, we employ a classical disturbance argument which is compared with
subelliptic Laplacian due to the lack of regularity to subelliptic systems
of $p$-Laplace. Indeed, our approach on the singularity ($1<p<2$) and
degeneracy ($p>2$) for $A$-harmonic systems \eqref{eq1} was essentially
influenced by way of a comparison with sub-Laplacian while $p$ is close to $2$
because of the Cordes conditions. This is an important technique to
consider PDEs with wild coefficients, also see \cite{CaP98}.

With Theorem \ref{th1} in hand, as a direct consequence we can obtain an
interior H\"older continuity of weak
solutions of subelliptic systems \eqref{eq1} while $Q-n<\lambda<p$.
Let us first recall the concept of H\"older space under the Carnot-Caratheodory
metric.


\begin{definition}[H\"older space] \label{holder} \rm
Let $\Omega \Subset \mathcal{G}$ and $0<\alpha<1$. We say that
$u\in \Gamma^{0,\alpha}_X(\Omega)$ has norm
$\|u\|_{\Gamma^{0,\alpha}_X(\Omega)}$, if
\begin{equation}\label{Holder}
\|u\|_{\Gamma^{0,\alpha}_X(\Omega)}
:=\sup_{\Omega}|u|+\sup_{\Omega}\frac {|u(x)-u(y)|}{[d_X(x,y)]^{\alpha}}<\infty.
\end{equation}
\end{definition}

Now we state the interior H\"older continuity of weak solutions with a sharp
index to subelliptic systems \eqref{eq1} as follows.

\begin{theorem}\label{th2}
Let $u\in HW^{1,p}(\Omega,\mathbb{R}^{N})$ is a weak solution of \eqref{eq1} with
$p$ close to 2. Suppose $A(x)$ and $B(x,u,Xu)$ satisfy H1-H3. If $Q-n<\lambda<p$,
we have
$$
u\in \Gamma_{X,{\rm loc}}^{0,\alpha}(\Omega,\mathbb{R}^{N})
$$
 with $\alpha=1-\frac{\lambda}{p}$.
\end{theorem}


The remainder of this paper is organized as follows.
In Section 2, we present some preliminaries concerning the sub-elliptic
setting and some several technical lemmas. In Section 3, we are devoted
to proving our main results.


\section{Preliminaries}

We adopt the usual convention of denoting by $C$ a general constant,
 which may vary from line to line in the same chain of inequalities.
Now let us first recall the Sobolev embedding inequality with respect to
 the horizontal vector fields, see \cite{CaDG,CaH,GaN}.

\begin{lemma} \label{Sobolev inequality}
 Let $1\leq p < Q $ and $ 1\leq q\leq \frac{Qp}{Q-p}$, where $Q$ is the
 homogeneous dimension of $X$ in $\mathcal{G}$. Then we have:
\begin{itemize}
\item[(1)] If $u(x)\in HW^{1,p}(B_{R_0},X)$, then there exists a positive constant\\
$C=C(p,q,Q,X)$  such that for any $0<R<R_{0}$,
      \begin{equation}
      \|u-\bar{u}_{R}\|_{L^q(B_{R})}
\leq CR^{1+Q(\frac 1q-\frac 1p)}\|Xu\|_{L^p(B_{R})},
      \end{equation}
  where $\bar{u}_{R}=\frac{1}{|B_R|}\int_{B_R}u dx$.

\item[(2)] If $u\in HW_{0}^{1,q}(B_{R_0},X)$, then there exists a positive constant\\
$C=C(p,q,Q,X)$ such that for any $0<R<R_{0}$,
  \begin{equation}
       \Big(\hbox{--}\hskip-9pt\int_{B_{R}}|u|^q\Big)^{1/q}
\leq CR\Big(\hbox{--}\hskip-9pt\int_{B_{R}}|Xu|^p\Big)^{1/p}.
   \end{equation}
\end{itemize}
\end{lemma}

To obtain the higher integrability of the horizontal gradients of solutions
for \eqref{eq1} we recall the following  reverse H\"{o}lder inequality originated
from Gehring's celebrating work on quasiconformal mappings, see
\cite[ Theorem 2.3 of Chapter 5]{Gi}.

\begin{lemma}\label{Reverse}
Suppose that $h(x)$ and $f(x)$ are nonnegative measurable functions satisfying
$ h(x)\in L^t(\Omega)$ and $f(x)\in L^s(\Omega)$ with $t>s>1$.
If for  all $x_0\in \Omega$ and all
$R:0<R<R_0\leq\operatorname{dist}(x_0,\partial \Omega)$ there holds
\begin{equation}\label{Reverse1}
\hbox{--}\hskip-9pt\int_{B_{\frac R 2}(x_0)}f^s\,dx
\leq b\Big(\Big\{\hbox{--}\hskip-9pt\int_{B_R(x_{0}) }f\,dx\Big\}^s+
\hbox{--}\hskip-9pt\int_{B_R(x_{0})}h^s\,dx\Big)+\theta\hbox{--}\hskip-9pt\int_{B_R(x_{0}) }f^s\,dx,
\end{equation}
with constants  $b>1$ and $0\leq \theta <1$, then there exist positive
constants $\delta=\delta(b,Q,q,s)$ and $C=C(b,Q,q,r)$ such that
$f\in L^{t}_{\rm loc}(\Omega)$ for any $t\in [s,s+\delta)$ and
\begin{equation}\label{Reverse2}
\Big\{\hbox{--}\hskip-9pt\int_{B_{\frac R 2}(x_0)}f^t\,dx\Big\}^{1/t}
\leq C\Big\{\hbox{--}\hskip-9pt\int_{B_R(x_0)}f^s\,dx\Big\}^{1/s}
+C\Big\{\hbox{--}\hskip-9pt\int_{B_R(x_0) }h^t\,dx\Big\}^{1/t}.
\end{equation}
\end{lemma}

With the reverse H\"{o}lder inequality above in hand, We can obtain the
following higher integrability of the horizontal gradients to
systems \eqref{eq1}.

\begin{lemma}[Higher integrability] \label{higher-integrability}
 Let $u\in HW^{1,p}(\Omega)$ be any weak solution of quasilinear subelliptic
systems \eqref{eq1} with  $A(x), B(x,u,Xu)$ satisfying assumptions  {\rm(H1)} and
{\rm (H3)}. Then, there exists a higher exponent $r: p<r<p+\delta$ such that for
 $\Omega'\Subset\Omega$, we have $Xu\in HW^{1,r}(\Omega')$.
Moreover, there exists a positive constant $C=C(Q,L,\nu,p)$ such that for any
$B_R(x_0)\Subset \Omega$ with the estimate
\begin{equation}\label{higher-integr}
\begin{aligned}
&\Big (\hbox{--}\hskip-9pt\int_{B_{\frac{R}{2}}(x_0)}|Xu|^{r})dx \Big)^{1/r}\\
&\leq C\Big(\hbox{--}\hskip-9pt\int_{B_R(x_0)}|Xu|^pdx\Big)^{1/p}
+C\Big(R\Big(\hbox{--}\hskip-9pt\int_{B_R(x_0)}|g(x)|^qdx\Big)^{1/q}\Big)^{\frac{1}{p-1}}.
\end{aligned}
\end{equation}
\end{lemma}

\begin{proof}
Given any  $x_0\in \Omega$, we take $R>0$ such that
$B_{R}:=B_{R}(x_0)\Subset \Omega$.
 Let $\eta$ be a cutting-off function with $\eta\in C^{\infty}_{0}(B_R)$
such that $0\le \eta(x)\le1$, $\eta=1$ for
$x\in B_{R/2}$, $\eta=0$ for $x\in \mathbb{R}^{n}\setminus{\overline{B}_{R}}$ and
$|X\eta|\leq\frac{C}{R}$. Let us take a test function $\varphi=\eta^p(u-\bar u_R)$
in \eqref{eq1}, it follows from \eqref{eq2} that
\begin{align*}
&\int_{\Omega}\Big\langle \langle A(x)Xu,Xu\rangle^{\frac{p-2}{2}}A(x)Xu,
\eta^pXu+p\eta^{p-1}(u-\bar u_{R})X\eta\Big \rangle dx\\
&=\int_{\Omega}B(x,u,Xu)\eta^p(u-\bar u_R)dx,
\end{align*}
By considering the uniformly ellipticity (H1) and the controllable growth, it yields
\begin{equation}\label{represent}
\begin{aligned}
&\nu^{p/2} \int_{\Omega}\eta^p|X u|^pdx\\
&\le \int_{\Omega}\eta^p\langle A(x)Xu, Xu\rangle^{p/2} dx \\
&=- \int_{\Omega}\langle \langle A(x)Xu,Xu\rangle
 ^{\frac{p-2}{2}}A(x)Xu,p\eta^{p-1}(u-\bar u_{R})X\eta \rangle dx\\
&\quad +\int_{\Omega}B(x,u,Xu)\eta^p(u-\bar u_R)dx \\
&\le pL^{p/2}\int_{\Omega} |\eta X u|^{p-1}|(u-\bar u_R) X\eta|dx
 +\mu\int_{\Omega}\eta^p|u-\bar u_R||Xu|^{p(1-\frac{1}{\gamma})}dx \\
&\quad +\mu\int_{\Omega}\eta^p|u-\bar u_R||g(x)|dx \\
&:=  pL^{p/2}I_{1}+\mu I_{2}+\mu I_{3}.
\end{aligned}
\end{equation}

Next, we  estimate $I_1, I_2$ and $I_3$. For  $I_1$, using
Young inequality with $\varepsilon_{1}>0$ and Sobolev inequality we have
\begin{equation}\label{repre-1}
\begin{aligned}
I_{1}&=  \int_{\Omega}|\eta X u|^{p-1}|(u-\bar u_R) X\eta|dx \\
&\leq  \varepsilon_{1}\int_{\Omega}|\eta Xu|^pdx
 +C(\varepsilon_{1})\int_{\Omega}|(u-\bar u_R)X\eta|^pdx  \\
&\leq  \varepsilon_{1}\int_{B_{R}}| Xu|^pdx
 +\frac {C(\varepsilon_{1})}{R^p}\int_{B_{R}}|u-\bar u_R|^pdx \\
&\leq  \varepsilon_{1}\int_{B_{R}}| Xu|^pdx
 +\frac {C(\varepsilon_{1})}{R^p}
 \Big(\int_{B_{R}}|Xu|^{\frac{Qp}{Q+p}}dx\Big)^{\frac{Q+p}{Q}}.
\end{aligned}
\end{equation}
For estimating $I_2$, by Sobolev inequality and H\"older inequality, it follows that
\begin{equation}\label{repre-2}
\begin{aligned}
I_{2}&= \int_{\Omega}\eta^p|u-\bar u_R||Xu|^{p(1-\frac{1}{\gamma})}dx \\
&\leq
\Big(\int_{B_{R}}|u-u_{R}|^{\gamma}dx\Big )^{1/\gamma}
\Big (\int_{B_{R}}|Xu|^pdx\Big)^{1-\frac{1}{\gamma}}  \\
&\leq CR^{1+Q(\frac{1}{\gamma}-\frac{1}{p})}
\Big (\int_{B_{R}}|Xu|^pdx\Big )^{\frac{1}{p}-\frac{1}{\gamma}}
\Big (\int_{B_{R}}|Xu|^pdx\Big ),
\end{aligned}
\end{equation}
where $\gamma\ge p$ is defined as the assumption H3  with
$1+Q(\frac{1}{p}-\frac{1}{\gamma})\ge 0$.

Similarly,  to estimate $I_{3}$  we  use H\"older inequality, Young inequality
 with  $\varepsilon_{2}>0$ and Sobolev inequality;  it yields
\begin{align*}
I_{3}&= \int_{\Omega}\eta^p|u-\bar u_R||g(x)|dx \\
&\leq \Big(\int_{B_{R}}|u-\bar u_R|^{\gamma}dx\Big)^{1/\gamma}
 \Big(\int_{B_{R}}|g(x)|^{\frac{\gamma}{\gamma-1}}dx\Big)^{\frac{\gamma-1}{\gamma}} \\
&\leq CR^{1+Q(\frac{1}{\gamma}-\frac{1}{p})}
 \Big(\int_{B_{R}}|Xu|^pdx\Big)^{1/p}\Big(\int_{B_{R}}
 |g(x)|^{\frac{\gamma}{\gamma-1}}dx\Big)^{\frac{\gamma-1}{\gamma}} \\
&\leq  \varepsilon_{2}\int_{B_{R}}|Xu|^pdx +C(\varepsilon_{2})
 R^{\frac{p}{p-1}\left(1+Q(\frac{1}{\gamma}-\frac{1}{p})\right)}
 \Big(\int_{B_{R}}|g(x)|^{\frac{\gamma}{\gamma-1}}dx\Big)^{\frac{\gamma-1}{\gamma}
\cdot\frac{p}{p-1}}.
\end{align*}
Now let us put the estimates of $I_1,I_2,I_3$ together into \eqref{represent},
we obtain
\begin{align*}
&\int_{B_{R}}|\eta X u|^pdx \\
&\leq  \frac {C(L,p,\varepsilon_{1})}{ R^p}
\Big(\int_{B_{R}}|Xu|^{\frac{Qp}{Q+p}}dx\Big)^{\frac{Q+p}{Q}}\\
&\quad +C(\mu,\varepsilon_{2})R^{\frac{p}{p-1}
\big(1+Q(\frac{1}{\gamma}-\frac{1}{p})\big)}
\Big(\int_{B_{R}}|g(x)|^{\frac{\gamma}{\gamma-1}}dx\Big)^{\frac{\gamma-1}{\gamma}
\cdot\frac{p}{p-1}}\\
 &\quad +\Big\{\mu\varepsilon_{2}+pL^{p/2}\varepsilon_{1}
+\mu C R^{1+Q(\frac{1}{\gamma}-\frac{1}{p})}
\Big(\int_{B_{R}}|Xu|^pdx\Big )^{\frac{1}{p}-\frac{1}{\gamma}}
  \Big\}\Big (\int_{B_{R}}|Xu|^pdx\Big).
\end{align*}
Let us write
\[
\vartheta=\mu\varepsilon_{2}+pL^{p/2}\varepsilon_{1}
+\mu C R^{1+Q(\frac{1}{\gamma}-\frac{1}{p})}
\Big (\int_{B_{R}}|Xu|^pdx\Big )^{\frac{1}{p}-\frac{1}{\gamma}}.
\]
 Notice that from the absolute continuity of the Lebesgue integral, we have that
$R^{1+Q(\frac{1}{\gamma}-\frac{1}{p})}\int_{B_R}|Xu|^p\to0$
as $R\to0$. Consequently we can take small $R>0$ such that $0<\vartheta<1$, and
\begin{align*}
\int_{B_{\frac{R}{2}}}|X u|^pdx
&\leq \frac{C}{R^p}\Big(\int_{B_{R}}|Xu|^{\frac{Qp}{Q+p}}dx\Big)^{\frac{Q+p}{Q}}\\
&\quad +CR^{\frac{p}{p-1}\left(1+Q(\frac{1}{\gamma}-\frac{1}{p})\right)}
\Big(\int_{B_{R}}|g(x)|^{\frac{\gamma}{\gamma-1}}dx\Big)^{\frac{\gamma-1}{\gamma}
\cdot\frac{p}{p-1}}+\vartheta\int_{B_{R}}|Xu|^pdx,
\end{align*}
which implies
\begin{equation*}
\hbox{--}\hskip-9pt\int_{B_{\frac{R}{2}}}|X u|^pdx
\leq C\Big(\hbox{--}\hskip-9pt\int_{B_{R}}|Xu|^{\frac{Qp}{Q+p}}dx\Big)^{\frac{Q+p}{Q}}
+C\Big(\hbox{--}\hskip-9pt\int_{B_{R}}|Rg(x)|^{\gamma'}dx\Big)^{\frac{p'}{\gamma'}}
+\vartheta\hbox{--}\hskip-9pt\int_{B_{R}}|Xu|^pdx,
\end{equation*}
with $p'=\frac p{p-1}$ and $\gamma'=\frac {\gamma}{\gamma-1}$.  Therefore, we obtain
\begin{align*}
\Big(\hbox{--}\hskip-9pt\int_{B_{\frac{R}{2}}}|X u|^pdx\Big)^{1/p}
&\leq C\Big(\hbox{--}\hskip-9pt\int_{B_{R}}|Xu|^{\frac{Qp}{Q+p}}dx\Big)^{\frac{Q+p}{pQ}}
+C\Big(\hbox{--}\hskip-9pt\int_{B_{R}}(|Rg(x)|^{\gamma'/p})^p dx\Big)
^{\frac 1p{\frac{p'}{\gamma'}}}\\
&\quad +\vartheta^{1/p} \Big(\hbox{--}\hskip-9pt\int_{B_{R}}|Xu|^pdx\Big )^{1/p}.
\end{align*}
Using the reverse H\"older inequality of Lemma \ref{Reverse}, it yields
\begin{equation}
\Big(\hbox{--}\hskip-9pt\int_{B_{\frac{R}{2}}}|X u|^rdx\Big)^{1/r}
\leq C\Big(\hbox{--}\hskip-9pt\int_{B_{R}}|Xu|^pdx\Big)^{1/p}
+C\Big(\hbox{--}\hskip-9pt\int_{B_{R}}(|Rg(x)|^{\gamma'/p})^r dx\Big)^{\frac 1r{\frac{p'}{\gamma'}}},
\end{equation}
for some  $p<r\le \frac {pq(\gamma-1)}{\gamma}$ due to
$q> \frac {\gamma}{\gamma-1}$.
Note that
\begin{align*}
\Big(\hbox{--}\hskip-9pt\int_{B_{R}}(|Rg(x)|^{\gamma'/p})^r dx\Big)^{\frac  1r{\frac{p'}{\gamma'}}}
&= R^{\frac{1}{p-1}}\Big(\hbox{--}\hskip-9pt\int_{B_{R}}|g(x)|^{\frac{\gamma}{\gamma-1}\cdot\frac{r}{p}}dx\Big)^{\frac{p}{r}\cdot\frac{\gamma-1}{\gamma}\cdot\frac{1}{p-1}}\\
&\le  R^{\frac{1}{p-1}}\Big(\hbox{--}\hskip-9pt\int_{B_{R}}|g(x)|^q dx
\Big)^{\frac 1q\cdot\frac{1}{p-1}},
\end{align*}
because $r\le \frac {pq(\gamma-1)}{\gamma}$, then we obtain
\eqref{higher-integr} which completes the proof.
\end{proof}

The following elementary inequalities concerning $A(x)$ are useful to our
 main proof, see \cite{KZh}.

\begin{lemma}\label{elementary}
Suppose that $A=(A_{ij})$ is a symmetric matrix and satisfies uniform ellipticity
{\rm (H1)}. Then there exists a positive
constant $C=C(p,\nu,L)$ such that for $1<p<\infty$ we have
\begin{equation}\label{ele1}
\big\langle\langle A\xi,\xi\rangle^{\frac{p-2}{2}}A\xi
-\langle A\eta,\eta\rangle^{\frac{p-2}{2}}A\eta,\xi-\eta\big\rangle
\geq C(|\xi|^2+|\eta|^2)^{\frac{p-2}{2}} |\xi-\eta|^2.
\end{equation}
In addition, for $p\geq2$ there holds
\begin{equation}\label{ele2}
\big\langle\langle A\xi,\xi\rangle^{\frac{p-2}{2}}A\xi
-\langle A\eta,\eta\rangle^{\frac{p-2}{2}}A\eta,\xi-\eta\big\rangle
\geq \nu^{p/2}|\xi-\eta|^p;
\end{equation}
and for $1<p<2$, there exists $C=C(p,\nu,L)$ such that for every
$0<\varepsilon<1$ we have
\begin{equation}\label{ele3}
|\xi-\eta|^p\le C\varepsilon^{\frac{p-2}{p}}
\big\langle\langle A\xi,\xi\rangle^{\frac{p-2}{2}}A\xi
-\langle A\eta,\eta\rangle^{\frac{p-2}{2}}A\eta,\xi-\eta\big\rangle
+\varepsilon|\eta|^p.
\end{equation}
\end{lemma}

To use  Campanato's freezing argument, we obverse the following local
 Dirichlet problems of homogeneous elliptic systems with constant coefficients.
\begin{equation}\label{homo equation}
\begin{gathered}
 -X^{*}(\langle A_{R}Xv,Xv\rangle^{\frac{p-2}{2}}A_{R}Xv)=0, \quad x\in B_R\\
v=u, \quad x\in\partial B_R,
\end{gathered}
\end{equation}
where $A_R=\frac{1}{|B_R|}\int_{B_R}A(x)dx$ is the integral average of $A(x)$.
Now we recall \cite[Lemma 3.4]{ZhF15} while $p$ is close to 2,
and have the following perturbation estimates.

\begin{lemma}\label{homo}
Let $v\in HW^{1,p}(\Omega)$ be a weak solution to the Dirichlet problems
\eqref{homo equation} with $p$ close to 2. Then for any
$u\in HW^{1,p}(\Omega)$, there exists $C>0$ such that for any
$x_0\in \Omega$ we have
\begin{equation}\label{perturbe esti}
\int_{B_{\rho}(x_0)}|Xu|^pdx\leq C\big(\frac{\rho}{R}\big)^{Q}
\int_{B_{R}(x_0)}|Xu|^pdx+C\int_{B_{R}(x_0)}|Xu-Xv|^pdx, 
\end{equation}
\quad
for all $0<\rho<R\le R_0$.
\end{lemma}

Moreover, by a direct calculation we obtain the following conclusion,
see \cite[Lemma 6]{ZhF15}.

\begin{lemma}\label{homo1}
Let $v\in HW^{1,p}(\Omega)$ be any weak solution to the  Dirichlet problem
\eqref{homo equation} with any $x_0\in\Omega$ and $0<R\le R_0$.
Then we have
$$
\int_{B_R(x_0)}|Xv|^pdx\leq C\int_{B_R(x_0)}|Xu|^pdx.
$$
\end{lemma}

In addition, we  need the following iteration lemma from \cite{Gi} in
the proof of our main theorem.

\begin{lemma}\label{iteration}
Let $\Phi(\rho)$ be a non-negative and non-decreasing function
on $(0,R)$. Suppose that
$$
\Phi(\rho)\le A\big\{\Big(\frac
{\rho}R\Big)^{\alpha}+\epsilon\big\}\Phi(R)+BR^{\beta},\quad  \forall\,
0<\rho<R\le R_0=\operatorname{dist}(x_0,\partial\Omega),
$$
with  non-negative constants $A,B,\alpha$ and $\beta$, and $\alpha>\beta$.
Then there exist two constants $\epsilon_0=\epsilon_0(A,\alpha,\beta)$ and
$C=C(A,\alpha,\beta)$ such that for any $0<\epsilon<\epsilon_0$ we have
$$
\Phi(\rho)\le C\big\{\Big(\frac {\rho}R\Big)^{\beta}\Phi(R)+B\rho^{\beta}\big\},
$$
for any $0<\rho< R\leq R_0=\operatorname{dist}(x_0,\partial\Omega)$.
\end{lemma}

Finally, the following equivalence of spaces is useful to prove a
local H\"older continuity of the weak solutions based on the main
Theorem, see \cite{DoN,DiF}.

\begin{lemma}\label{equivalence}
If $0<\lambda<Q$, the Campanato space $\mathcal{L}^{ p,\lambda}_X(\Omega)$
is isomorphic to the Morrey space $L^{ p,\lambda}_X(\Omega)$.
If $-p<\lambda<0$, then the Campanato space $\mathcal{L}^{ p,\lambda}_X(\Omega)$
is isomorphic to the H\"{o}lder space $\Gamma^{0,\alpha}_X(\Omega)$ with
$\alpha=-\lambda/p$.
\end{lemma}

\section{Proof of the main results}

\begin{proof}[Proof of theorem \ref{th1}]
 Let $u(x)\in HW^{1,p}(\Omega)$ be any weak solution of \eqref{eq1}.
To obtain an interior estimate for the horizontal gradients to the solutions,
for any fixed point $x_0\in\Omega$ let us take a ball $B_R(x_0)\Subset \Omega$,
 and write $B_R=:B_R(x_0)$ in the context. Suppose that $v(x)$ is a weak solution
of the local Dirichlet problem \eqref{homo equation}. Computing the
difference between \eqref{eq1} and \eqref{homo equation} yields
\begin{align*}
&-X^{*}\Big (\langle A(x)Xu,Xu\rangle^{\frac{p-2}{2}}A(x)Xu
 -\langle A_R Xu,Xu\rangle^{\frac{p-2}{2}}A_R Xu\Big)\\
&-X^{*}\Big(\langle A_R Xu,Xu\rangle^{\frac{p-2}{2}}A_R Xu
 -\langle A_R Xv,Xv\rangle^{\frac{p-2}{2}}A_R Xv\Big)\\
&=-B(x,u,Xu).
\end{align*}
Let us take $\varphi=u-v$ as a test function in the weak sense to the equations
above, then we have
\begin{equation}\label{local-equ}
\begin{aligned}
&\int_{B_R}\big\langle\langle A_R Xu,Xu\rangle^{\frac{p-2}{2}}A_RXu
 -\langle A_RXv,Xv\rangle^{\frac{p-2}{2}}A_RXv,Xu-Xv\big\rangle dx  \\
&= \int_{B_R}\big\langle\langle A_RXu,Xu\rangle^{\frac{p-2}{2}}A_RXu
 -\langle A(x)Xu,Xu\rangle^{\frac{p-2}{2}}A(x)Xu,Xu-Xv\big\rangle dx  \\
&\quad +\int_{B_R}\langle B(x,u,Xu),u-v\rangle dx.
\end{aligned}
\end{equation}
Note that \eqref{perturbe esti} implies that if $p$ is close to $2$,
for any $0<\rho<R$ there holds
\begin{equation}\label{perturbe est}
\int_{B_{\rho}(x_0)}|Xu|^pdx\leq C\big(\frac{\rho}{R}\big)^{Q}
\int_{B_{R}(x_0)}|Xu|^pdx+C\int_{B_{R}(x_0)}|Xu-Xv|^pdx.
\end{equation}

Next, we focus on the estimate of $\int_{B_{R}(x_0)}|Xu-Xv|^pdx$.
To that end, we will estimate it by dividing into two cases.
\smallskip

\noindent\textbf{Case 1: $p\geq2$.}
 Let us put an elementary inequality \eqref{ele2} and  (H3) into the
equation \eqref{local-equ} above it follows that
\begin{align}
&\nu^{p/2}\int_{B_{R}}|Xu-Xv|^pdx \nonumber \\
&\leq \int_{B_{R}}\big\langle\langle A_RXu,Xu\rangle^{\frac{p-2}{2}}A_RXu
 -\langle A_RXv,Xv\rangle^{\frac{p-2}{2}}A_R Xv,Xu-Xv\big\rangle dx \nonumber \\
&=\int_{B_{R}}\big\langle\langle A_RXu,Xu\rangle^{\frac{p-2}{2}}A_RXu
 -\langle A(x)Xu,Xu\rangle^{\frac{p-2}{2}}A(x)Xu,Xu-Xv\big\rangle dx \nonumber \\
&\quad + \int_{B_R}\langle B(x,u,Xu),u-v\rangle dx \nonumber \\
&\leq C\int_{B_{R}}|A_R-A(x)||Xu|^{p-1}|Xu-Xv|dx
 +\mu\int_{B_R}|Xu|^{p(1-\frac{1}{\gamma})}\cdot|u-v|dx \nonumber \\
&\quad +\mu\int_{B_R}|g|\cdot|u-v|dx \nonumber \\
&:= J_1+J_2+J_3. \label{J123}
\end{align}
For estimating  $J_1$, the H\"{o}lder inequality and Young inequalities
with any $\varepsilon>0$ yield
\begin{align*}
&\int_{B_{R}}|A_R-A(x)||Xu|^{p-1}|Xu-Xv|dx\\
&\leq \Big(\int_{B_{R}}|A_R-A(x)|^{\frac{p}{p-1}}|Xu|^pdx\Big)^{1-\frac{1}{p}}
 \Big(\int_{B_R}|Xu-Xv|^p dx\Big)^{1/p}\\
&\leq C(\varepsilon_4)\int_{B_{R}}|A_R-A(x)|^{\frac{p}{p-1}}|Xu|^pdx+\varepsilon_4\int_{B_R}|Xu-Xv|^pdx\\
&\leq C(\varepsilon_4)|B_R|\Big(\hbox{--}\hskip-9pt\int_{B_{R}}|A_R-A(x)|^{\frac{p}{p-1}
 \cdot\frac{r}{r-p}}dx\Big)^{1-\frac{p}{r}}
 \Big(\hbox{--}\hskip-9pt\int_{B_R}|Xu|^{r}dx\Big)^{p/r}\\
&\quad +\varepsilon_4\int_{B_R}|Xu-Xv|^pdx,
\end{align*}
where an $r>p$ is the same integrable index as that in
Lemma \ref{higher-integrability}. Setting $t=\frac{pr}{(p-1)(r-p)}$,
and by a higher integrability from Lemma \ref{higher-integrability} we obtain
\begin{align*}
J_1 
&\leq C(\varepsilon_4)|B_R|\Big(\hbox{--}\hskip-9pt\int_{B_{R}}|A_R-A(x)|^tdx\Big)^{1-\frac{p}{r}}
\Big(\hbox{--}\hskip-9pt\int_{B_R}|Xu|^pdx+R^{\frac{p}{p-1}}
 \Big(\hbox{--}\hskip-9pt\int_{B_R}|g|^q dx\Big)^{\frac{1}{q}
 \cdot\frac{p}{p-1}}\Big) \\
&\quad +\varepsilon_4\int_{B_R}|Xu-Xv|^pdx \\
&\leq C(\varepsilon_4)M_{A}(R)^{1-\frac{p}{r}}\int_{B_R}|Xu|^pdx
 +CR^{Q+\frac{p}{p-1}\cdot\frac{q-Q}{q}}\|g\|_{L^q}^{\frac{p}{p-1}}\\
&\quad +\varepsilon_4\int_{B_R}|Xu-Xv|^pdx.
\end{align*}
To estimate $J_2$, we employ H\"{o}der inequality and Young inequality again,
and obtain
\begin{align*}
J_2
&\leq \mu\Big(\int_{B_R}|Xu|^p dx\Big)^{1-\frac{1}{\gamma}}
 \Big(\int_{B_R}|u-v|^{\gamma}dx\Big)^{1/\gamma}\\
&\leq \varepsilon_5\int_{B_R}|Xu|^pdx
 +C(\mu,Q,p,\varepsilon_5)\int_{B_R}|u-v|^{\gamma}dx\\
&\leq \varepsilon_5\int_{B_R}|Xu|^pdx+CR^{\gamma+Q(1-\frac{\gamma}{p})}
\Big(\int_{B_R}|Xu-Xv|^pdx\Big)^{\frac{\gamma-p}{p}}\int_{B_R}|Xu-Xv|^pdx.
\end{align*}
Observing $\delta(R):=R^{\gamma+Q(1-\frac{\gamma}{p})}
\Big(\int_{B_R}|Xu-Xv|^pdx\Big)^{\frac{\gamma-p}{p}}\to 0$ as $R\to 0$,
then there holds
\begin{equation}
J_2\leq C\delta(R)\int_{B_R}|Xu-Xv|^pdx+\varepsilon_5\int_{B_R}|Xu|^pdx.
\end{equation}
To estimate $J_3$, by using H\"{o}lder inequality, Sobolev embedding inequality and Young inequality it follows that
\begin{align*}
J_3
&\leq  \mu\Big(\int_{B_R}|g|^{\frac{\gamma}{\gamma-1}}dx\Big)
 ^{\frac{\gamma-1}{\gamma}} \Big(\int_{B_R}|u-v|^{\gamma}dx\Big)^{1/\gamma} \\
&\leq  \Big(\int_{B_R}|g|^{\frac{\gamma}{\gamma-1}}dx\Big)^{\frac{\gamma-1}{\gamma}}
  R^{1+Q(\frac{1}{\gamma}-\frac{1}{p})}\Big(\int_{B_R}|Xu-Xv|^pdx\Big)^{1/p} \\
&\leq \varepsilon_6\int_{B_R}|Xu-Xv|^pdx+C(\varepsilon_6)R^{\frac{p}{p-1}
[1+Q(\frac{1}{\gamma}-\frac{1}{p})]}
\Big(\int_{B_R}|g|^{\frac{\gamma}{\gamma-1}}dx\Big)^{\frac{\gamma-1}{\gamma}
\cdot\frac{p}{p-1}} \\
&\leq \varepsilon_6\int_{B_R}|Xu-Xv|^pdx
+CR^{Q+\frac{p}{p-1}\cdot\frac{q-Q}{q}}\|g\|_{L^q}^{\frac{p}{p-1}}.
\end{align*}

Putting   estimates of $J_1,J_2$ and $J_3$ together in \eqref{J123}, one deduces
\begin{align*}
&\nu^{p/2}\int_{B_R}|Xu-Xv|^pdx\\
&\leq C\big(\delta(R)+\varepsilon_4+\varepsilon_6\big)
 \int_{B_R}|Xu-Xv|^pdx+\Big(C(\varepsilon_4)
M_{A}(R)^{1-\frac{p}{r}}+\varepsilon_5\Big)\int_{B_R}|Xu|^pdx\\
&\quad +CR^{Q+\frac{p}{p-1}\cdot\frac{q-Q}{q}}\|g\|_{L^q}^{\frac{p}{p-1}}.
\end{align*}
Therefore, by choosing arbitrary positive constants
$\varepsilon_4,\varepsilon_6$ and $0<R<R_0$ small enough that
$C(\delta(R)+\varepsilon_4+\varepsilon_6)<\nu^{p/2}$ we obtain
\begin{equation}\label{D-u-v}
\begin{aligned}
&\int_{B_R}|Xu-Xv|^pdx \\
&\leq \left(C(\varepsilon_4)M_{A}(R)^{1-\frac{p}{r}}+\varepsilon_5\right)
 \int_{B_R}|Xu|^pdx
 +CR^{Q+\frac{p}{p-1}\cdot\frac{q-Q}{q}}\|g\|_{L^q}^{\frac{p}{p-1}} \\
&\le \left(C(\varepsilon_4)M_{A}(R)^{1-\frac{p}{r}}+\varepsilon_5\right)
 \int_{B_R}|Xu|^pdx
+CR^{Q+\frac{p}{p-1}\cdot\frac{q-\mu}{q}}\|g\|_{L^{q,\mu}_X}^{\frac{p}{p-1}}.
\end{aligned}
\end{equation}
Let us set $\varpi=C(\varepsilon_4)M_{A}(R)^{1-\frac{p}{r}}+\varepsilon_5$
in \eqref{D-u-v}, and put the estimate \eqref{D-u-v} into \eqref{perturbe est},
then for any $0<\rho<R$ we  have
\begin{equation} \label{4.3}
\int_{B_\rho}|Xu|^pdx
\leq C\big[\big(\frac{\rho}{R}\big)^Q+\varpi\big]
\int_{B_R}|Xu|^pdx+CR^{Q+\frac{p}{p-1}\cdot\frac{q-\mu}{q}}
 \|g\|_{L^{q,\mu}_X}^{\frac{p}{p-1}}.
\end{equation}
While $q\ge\mu$ such that $Q+\frac{p}{p-1}\cdot\frac{q-\mu}{q}\ge Q$,
it follows from Lemma \ref{iteration} that
\begin{equation}
\int_{B_\rho}|Xu|^pdx\le C\big(\frac{\rho}{R}\big)^{Q-\lambda}
\int_{B_R}|Xu|^pdx+C\rho^{Q-\lambda}\|g\|^{\frac{p}{p-1}}_{L_X^{q,\mu}(B_R)},
\end{equation}
for any $0<\lambda<Q$. This implies $Xu\in L_{X,{\rm loc}}^{p,\lambda}(\Omega')$
for $\Omega'\Subset\Omega$ with the estimate
$$
\|Xu\|_{L_X^{p,\lambda}(\Omega')}\leq C\Big(\|Xu\|_{L^p(\Omega)}
+\|g\|^{\frac{1}{p-1}}_{L_X^{q,\mu}(\Omega)}\Big).
$$
While $\frac {\gamma}{\gamma-1}<q<\mu$ such that
$Q+\frac{p}{p-1}\cdot\frac{q-\mu}{q}<Q$, then we deduce from
Lemma \ref{iteration} that
\begin{equation} %\label{4.3}
\int_{B_\rho}|Xu|^pdx
\leq C\big(\frac{\rho}{R}\big)^{Q-\frac{p}{p-1}\frac{\mu-q}{q}}
\int_{B_R}|Xu|^pdx+C\rho^{Q-\frac{p}{p-1}\cdot\frac{\mu -q}{q}}
\|g\|_{L^{q,\mu}_X}^{\frac{p}{p-1}},
\end{equation}
which implies $Xu\in L_{X,{\rm loc}}^{p,\frac{p}{p-1}\frac{\mu-q}{q}}(\Omega')$
with the estimate
\begin{align*} % \label{4.3}
\|Xu\|_{L_X^{p,\frac{p}{p-1}\frac{\mu-q}{q}}(\Omega')}
\leq C\Big(\|Xu\|_{L^p(\Omega)}+\|g\|^{\frac{1}{p-1}}_{L_X^{q,\mu}(\Omega)}\Big).
\end{align*}
\smallskip

\noindent\textbf{Case 2:  $1<p<2$.} Using  inequality \eqref{ele3}  yields
 \begin{align*}
&\int_{B_{\rho}}|Xu-Xv|^pdx\\
&\leq C\varepsilon^{\frac{p-2}{2}}\int_{B_{\rho}}
\big\langle\langle A_RXu,Xu\rangle^{\frac{p-2}{2}}A_RXu
 -\langle A_RXv,Xv\rangle^{\frac{p-2}{2}}A_RXv,Xu-Xv\big\rangle dx \\
&\quad +\varepsilon\int_{B_\rho}|Xu|^p dx\\
&= C(\varepsilon) \int_{B_{\rho}}
 \big\langle\langle A_RXu,Xu\rangle^{\frac{p-2}{2}}A_RXu
 -\langle A(x)Xu,Xu\rangle^{\frac{p-2}{2}}A(x)Xu,Xu-Xv\big\rangle dx\\
&\quad + B(x,u,Xu) +\varepsilon\int_{B_\rho}|Xu|^p dx\\
&\leq C(p,\nu,L)\int_{B_{R}}|A_R-A(x)||Xu|^{p-1}|Xu-Xv|dx
 +\mu\int_{B_R}|Xu|^{p(1-\frac{1}{\gamma})}\cdot|u-v|dx\\
&\quad  +\mu\int_{B_R}|g|\cdot|u-v|dx
  +\varepsilon\int_{B_R}|Xu|^p dx\\
&:= J_1+J_2+J_3+\varepsilon\int_{B_R}|Xu|^p dx.
\end{align*}
Considering the estimates of $J_1,J_2,J_3$ in Case 1, and for any
$\varepsilon>0$ we obtain
\begin{equation} %\label{4.3}
\int_{B_\rho}|Xu|^pdx
\le C\big[\big(\frac{\rho}{R}\big)^Q+\varpi'\big]
\int_{B_R}|Xu|^pdx+CR^{Q+\frac{p}{p-1}\cdot\frac{q-\mu}{q}}\|g\|_{L^{q,\mu}_X}
^{\frac{p}{p-1}},
\end{equation}
with $\varpi'=C(\varepsilon_4)M_{A}(R)^{1-\frac{p}{r}}+\varepsilon_5+\varepsilon$.
While $q\ge\mu$, it follows  form Lemma \ref{iteration} as the same as Case 1 that
$$
\int_{B_\rho}|Xu|^pdx
\le C\big(\frac{\rho}{R}\big)^{Q-\lambda}\int_{B_R}|Xu|^pdx
+C\rho^{Q-\lambda}\|g\|^{\frac{p}{p-1}}_{L^q(B_R)},
$$
for any $0<\lambda<Q$. It yields $Xu\in L_{X,{\rm loc}}^{p,\lambda}(\Omega')$
for any $\Omega'\Subset \Omega$ with the estimate
$$
\|Xu\|_{L_X^{p,\lambda}(\Omega')}
\leq C\Big(\|Xu\|_{L^p(\Omega)}+\|g\|^{\frac{1}{p-1}}_{L_X^{q,\mu}(\Omega)}\Big).
$$
While $\frac {\gamma}{\gamma-1}<q<\mu$, by Lemma \ref{iteration} we obtain
\[ %\label{4.3}
\int_{B_\rho}|Xu|^pdx
\leq C\big(\frac{\rho}{R}\big)^{Q-\frac{p}{p-1}\frac{\mu-q}{q}}
\int_{B_R}|Xu|^pdx+C\rho^{Q-\frac{p}{p-1}\cdot\frac{\mu -q}{q}}
\|g\|_{L^{q,\mu}_X}^{\frac{p}{p-1}}.
\]
which implies $Xu\in L_{X,{\rm loc}}^{p,\frac{p}{p-1}\frac{\mu-q}{q}}(\Omega')$
with the estimate
\[ %\label{4.3}
\|Xu\|_{L_X^{p,\frac{p}{p-1}\frac{\mu-q}{q}}(\Omega')}\leq C\Big(\|Xu\|_{L^p(\Omega)}+\|g\|^{\frac{1}{p-1}}_{L_X^{q,\mu}(\Omega)}\Big).
\]
This completes the proof.
\end{proof}

\begin{proof}[Proof of Theorem \ref{th2}]
For any $0<\rho<R$ with $B_R(x_0)\Subset\Omega$, by Poincare inequality
 we have
\begin{equation}\label{Poincare}
\int_{B_\rho}|u-u_\rho|^pdx
\leq C\rho^p\int_{B_\rho}|Xu|^pdx
\leq C\rho^p\rho^{Q-\lambda}\|Xu\|^p_{L_X^{p,\lambda}(B_{\rho})}.
\end{equation}
Note that by Theorem \ref{th1},
\begin{equation}\label{Thm}
\|Xu\|^p_{L_X^{p,\lambda}}(B_{\rho})
\le C\big(\frac{1}{R}\big)^{Q-\lambda}
\int_{B_R}|Xu|^pdx+C\|g\|^{\frac{p}{p-1}}_{L_X^{q,\mu}(B_R)}.
\end{equation}
Therefore, combining \eqref{Poincare} and \eqref{Thm} yields
\begin{equation}
\frac{\rho^{\lambda-p}}{|B_\rho|}\int_{B_\rho}|u-u_\rho|^pdx
\le C\big(\frac{1}{R}\big)^{Q-\lambda}\int_{B_R}|Xu|^pdx
+C\|g\|^{\frac{p}{p-1}}_{L_X^{q,\mu}(B_R)},
\end{equation}
which implies
$$
u\in \mathcal{L}_{X,{\rm loc}}^{p,\lambda-p}(\Omega).
$$
Note that $Q-n<\lambda<p$ implies $-p<Q-n-p<\lambda-p<0$, it follows from
Lemma \ref{equivalence} that
$$
u\in \Gamma_{X,{\rm loc}}^{0,\alpha}(\Omega),
$$
with $\alpha=-\frac{\lambda-p}{p}=1-\frac{\lambda}{p}$.
\end{proof}

\begin{remark} \rm
We would like to point out that Theorems \ref{th1} and  \ref{th2} are
valid only under the assumption of $p$ close to 2 when we
consider subelliptic \eqref{eq1} in Carnot group.
Indeed, we employ the perturbation inequality \eqref{perturbe esti}
in our main proof, which is attained by  a local Lipschitz boundedness
of subelliptic $p$-harmonic  only if $p$ close to 2.
However, it is not necessary to limit $p$ close to 2 if
$X=\big(\frac \partial{\partial x_1},\frac \partial{\partial x_2},\dots,
\frac \partial{\partial x_n}\big)$ is a classical usual gradient in the
Euclidian spaces $\mathbb{R}^n$.
\end{remark}

\subsection*{Acknowledgments}
This work is supported by NSF of China under 11371050, and
by the  NSF of China under grants 2013AA013702 (863) and 2013CB834205 (973).


\begin{thebibliography}{99}

\bibitem{BeF}  Bensoussan, A.; Frehse, J.;
\emph{Regularity results for nonlinear elliptic systems and applications.}
Applied Math. Sci. Vol., 151, Springer-Verlag, Berlin, 2002.

\bibitem{BrB} Bramanti, M.; Brandolini,L.;
\emph{$L^p$-estimates for nonvariational hypoelliptic operators with
 VMO coefficients.} Trans. Amer. Math. Soc., 352 (2002), 781-822.

\bibitem{CaDG} Capogna, L.; Danielli, D.; Garofalo, N.;
\emph{An embedding theorem and the Harnack inequality for nonlinear
subelliptic equations.}  Comm. Partial Differential Equations, 18(9-10)(1993),
 1765-1794.

\bibitem{CaG} Capogna, L.;  Garofalo, N.;
\emph{Regularity of minimizers of the calculus of variations in Carnot groups
via hypoellipticity of systems of H\"ormander type.}  J. Eur. Math. Soc.,
 5(1)(2003),1-40.

\bibitem{CaH} Capogna, L.; Han, Q.;
\emph{Pointwise Schauder estimates for second order linear equations in Carnot groups.}
 Proceedings for the AMS-SIAM Harmonic Analysis conference in Mt. Holyhoke, 2001.

\bibitem{CaP98} Caffarelli, L. A.; Peral, I.;
\emph{On $W^{1,p}$ estimates for elliptic equations in divergence form.}
 Comm. Pure Appl. Math., 51(1)(1998), 0001-0021.

\bibitem{CiGL}  Citti, L.; Garofalo, N.; Lanconelli, E.;
\emph{HarnackĄ¯s inequality for sum of squares of vector fields plus a potential.}
 Amer. J. Math., 115(1993), 699-734.

\bibitem{DiF} Di Fazio, G.; Fanciullo, M. S.;
\emph{Gradient estimates for elliptic systems in Carnot-Carath\'{e}odory spaces.}
 Comment. Math. Univ. Caroline., 43 (4)(2002), 605-618.

\bibitem{DoN} Dong, Y., Niu, P. C.;
 \emph{Estimates in Morrey spaces and H\"{o}lder continuity for weak solutions
to degenerate elliptic systems.}   Manuscripta Math., 138(2012), 419-437.

\bibitem{DiZ}  Di Fazio, G.; Zamboni,P.;
\emph{H\"older continuity for quasilinear subelliptic eauations in
Carnot-Carath\'eodory spaces.}   Math. Nachr., 272(2004), 3-10.

\bibitem{DoM} Domokos, A.; Manfredi, J. J.;
\emph{Subelliptic Cordes estimates.} Proc. Amer. Math. Soc.,  133(2005), 1047-1056.

\bibitem{DoM1} Domokos, A.; Manfredi, J. J.;
\emph{$C^{1,\alpha}$-regularity for $p$-harmonic functions in the Heisenberg
group for $p$ near $2$.}  Contemp. Math., 370 (2005), 17-23.

\bibitem{DoM2} Domokos, A.;
\emph{On the regularity of p-harmonic functions in the Heisenberg group.}
 Doctoral Dissertation, University of Pittsburgh, 2004.

\bibitem{DoM3}  Domokos, A.;
\emph{On the regularity of subelliptic p-harmonic functions in Carnot groups.}
Nonlinear Analysis, 69(2008), 1744-1756.

\bibitem{Fol} Folland, G.;
\emph{A fundamental solution for a subelliptic operator.}
 Bull. Amer. Math. Soc.,  79(1973), 373-376.

\bibitem{GaN}  Garofalo, N.; Nhieu, D. M.;
\emph{Isoperimetric and Sobolev inequalities for Carnot-Carath\'eodory spaces
and the existence of minmal surfaces.} Comm. Pure Appl. Math., 49(10)(1996),
1081-1144.

\bibitem{Gi} Giaquinta, M.;
\emph{Multiple integrals in the calculus of variations and nonlinear elliptic
 systems.} in 'Ann. Math. Stud.', Princeton University Press, Vol. 105, 1983.

\bibitem{Ho} H\"ormander, L.;
\emph{Hypoelliptic second order differential
equations.}   Acta Math., 119(1967), 147-171.

\bibitem{JoX} Jost, J.; Xu, C. J.;
\emph{Subelliptic harmonic maps.}  Trans. Amer. Math. Soc., 350(11)(1998), 4633-4649.

\bibitem{KZh} Kinnunen, J.; Zhou, S. L.;
\emph{A local estimate for nonlinear equations with discontinuous coefficients.}
Comm. Part. Diff. Equ., 24(11-12)(1999), 2043-2068.

\bibitem{Lu} Lu, G.;
\emph{The sharp Poincar\'e inequality for free vector fields: An endpoint result.}
Rev. Mat. Ibe., 10(1994), 453-466.

\bibitem{MZZ} Mingione, G.; Zatorska-Goldstein, A.; Zhong, X.;
\emph{Gradient regularity for elliptic equations in the Heisenberg group.}
Adv. Math., 222 (2009) 62-129.

\bibitem{NaSW} Nagel, A.; Stein, E.;  Wainger, S.;
\emph{Balls and metrics defined by
vector fields I: basic properties.} Acta Math., 155(1985), 103-147.

\bibitem{RoS} Rothschild, L.; Stein, E.;
\emph{Hypoelliptic operators and nilpotent Lie
groups.} Acta Math., 137(1977), 247-320.

\bibitem{Uh}  Uhlenbeck, K.;
\emph{Regularity for a class of nonlinear elliptic systems.}
Acta Math., 138(1977), 219-240.

\bibitem{WaL} Wang, J.; Liao, D.;
\emph{Optimal partial regularity for sub-elliptic systems with sub-quadratic
growth in Carnot groups.} Nonlinear Analysis, 75(2012), 2499-2519.

\bibitem{XuZ}  Xu, C.; Zuily, C.;
\emph{Higher interior regularity for quasilinear subelliptic systems.}
Calc. Var. Partial Differential Equations,  5(1997), 323-343.

\bibitem{YuZ}  Yu, H.; Zheng, S.;
\emph{Optimal partial regularity for quasilinear ellipyic systems with
VMO coefficients based on A-harmonic approximations.} Electronic Journal
of Differential Equation, 16(2015), 1-12.

\bibitem{Zh3}  Zheng, S.; Feng, Z.;
\emph{Green functions for a class of  nonlinear degenerate
operators with X-ellipticity.} Trans. Amer. Math. Soc.,  364(7)(2012), 3627-3655.

\bibitem{ZhF15} Zheng, S.; Feng, Z.;
\emph{Regularity of  subelliptic p-harmonic systems with subcritical growth
in Carnot group.} J. Differential Equations, 258(2015), 2471-2494.

\end{thebibliography}

\end{document}
