\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2016 (2016), No. 323, pp. 1--24.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2016/323\hfil Elliptic system with non-standard growth]
{Partial continuity for a class of elliptic systems with non-standard growth}

\author[J. Ok \hfil EJDE-2016/323\hfilneg]
{Jihoon Ok}

\address{Jihoon Ok \newline
School of Mathematics,
Korea Institute for Advanced Study, Seoul 02455, Korea}
\email{jihoonok@kias.re.kr}

\thanks{Submitted October 24, 2016. Published December 20, 2016.}
\subjclass[2010]{35B65, 35J60}
\keywords{Partial continuity; elliptic system; non-standard growth;
\hfill\break\indent  variable exponent}

\begin{abstract}
 We study partial H\"older continuity of weak solutions to elliptic systems
 with variable non-standard growth, which are related to the function
 $\Phi(x,t):=t^{p(x)}\log(e+t)$.  We prove that weak solutions are H\"older
 continuous for any H\"older exponent,  except Lebesgue measure zero sets,
 if systems satisfy certain continuity assumptions.
 In particular, the variable exponent functions $p(\cdot)$ are assumed to satisfy
 so-called vanishing log-H\"older continuity.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\allowdisplaybreaks

\section{Introduction} \label{sec1}

The aim of this article is to  study  partial H\"older continuity of solutions to
elliptic systems with non-standard growth conditions, and, in particular, 
with a certain type of Orlicz growth conditions. Problems with non-standard 
growth conditions have initially studied in the fundamental papers of
 Marcellini \cite{M2, M3, M4, M5} and have subsequently been the object 
of intensive investigation. The new methods here will be explained for 
a special, yet significant case, incorporating both the features of 
purely Orlicz structures and the ones of variable exponent functionals. 
The same methods can be considered as a substantial starting point to treat 
more general structures including both Orlicz structures and the 
non-autonomous case, that in the case of non-standard growth conditions 
poses additional difficulties.

Specifically, let $\Omega\subset\mathbb{R}^n$ be a bounded open set in $\mathbb{R}^n$ 
with $n\geq 2$, and $p(\cdot):\Omega\to\mathbb{R}$ be a variable exponent function 
satisfying
\begin{equation}\label{pxgamma}
2<\gamma_1\leq p(\cdot)\leq \gamma_2<\infty.
\end{equation}
Moreover, we assume that $p(\cdot)$ satisfies the so-called vanishing 
log-H\"older continuity:
\begin{equation}\label{vlogholder}
\lim_{r\to0}\omega(r)\log\big(\frac{1}{r}\big)=0,
\end{equation}
where $\omega(\cdot):[0,\infty)\to[0,\infty)$ is the modulus of continuity of 
$p(\cdot)$, that is, it is a nondecreasing and concave function satisfying 
that $\omega(0)=0$ and $|p(x)-p(y)|\leq \omega(|x-y|)$ for every $x,y\in\Omega$. 
 With this $p(\cdot)$, we consider the elliptic system
\begin{equation}\label{maineq}
\operatorname{div} \mathbf{a}(x,u,Du)=0 \quad \text{in }\Omega.
\end{equation}
Here,  $\mathbf{a}:\Omega\times \mathbb{R}^N\times \mathbb{R}^{Nn}\to\mathbb{R}^{Nn}$, $N\geq1$, satisfies 
growth and ellipticity conditions:
\begin{equation}\label{acondition}
\begin{gathered}
|\mathbf{a}(x,\zeta,\xi)|+|\partial\mathbf{a}(x,\zeta,\xi)|(1+|\xi|)
\leq  \Lambda \Phi_{p(x)-1}(1+ |\xi|),\\
\partial\mathbf{a}(x,\zeta,\xi)\eta\cdot \eta  \geq  \nu \Phi_{p(x)-2}(1+|\xi|) |\eta|^2,
\end{gathered}
\end{equation}
and
\begin{equation} \label{xconti}
\begin{aligned}
& |\mathbf{a}(x_1,\zeta,\xi)-\mathbf{a}(x_2,\zeta,\xi)|&&\\
& \leq \Lambda \omega(|x_1-x_2|)[\Phi_{p(x_1)-1}(1+|\xi|)
+ \Phi_{p(x_2)-1}(1+|\xi|)]\log(1+|\xi|),
\end{aligned}
\end{equation}
for every $x,x_1,x_2\in\Omega$, $\zeta \in\mathbb{R}^N $ and $\xi,\eta\in\mathbb{R}^{nN}$
and for some $0<\nu\leq \Lambda$, where $\Phi_p:[0,\infty)\to[0,\infty)$, $p>0$,
is denoted by
\begin{equation}\label{phip}
\Phi_p(t):=t^p\log(e+t),
\end{equation}
and $\partial\mathbf{a}(x,\zeta,\xi):=D_\xi\mathbf{a}(x,\zeta,\xi)$.
We note from the ellipticity condition, that is the second inequality in
\eqref{acondition}, that one can see the  monotonicity
\begin{equation}\label{monotonicity}
\begin{aligned}
&(\mathbf{a}(x,\zeta,\xi_1)-\mathbf{a}(x,\zeta,\xi_2)):(\xi_1-\xi_2) \\
&\geq \tilde \nu  \Phi_{p(x)-2}(1+|\xi_1|+|\xi_2|) |\xi_1-\xi_2|^2\\
&\geq\ \frac{\tilde \nu}{2} \{\Phi_{p(x)-2} (1+|\xi_1|) |\xi_1-\xi_2|^2
+ \Phi_{p(x)} ( |\xi_1-\xi_2|)\}.
\end{aligned}
\end{equation}

For a function $p(\cdot)$, we set
\begin{equation}\label{phi}
\Phi(x,t):=\Phi_{p(x)}(t)=t^{p(x)}\log(e+t),
\end{equation}
and define the Musielak-Orlicz space $L^{\Phi}(\Omega,\mathbb{R}^N)$ by the set 
of all measurable vector valued functions $f:\Omega\to\mathbb{R}^N$ satisfying
$$
\int_\Omega\Phi(x,|f|)\, dx=\int_\Omega|f|^{p(x)}\log(e+|f|)\, dx<\infty.
$$
and $W^{1,\Phi}(\Omega,\mathbb{R}^N)$ by the set of all $f\in W^{1,1}(\Omega,\mathbb{R}^N)$
 with $|f|,|Df|\in L^{\Phi}(\Omega,\mathbb{R}^1)$.

Our purpose in this article is to investigate partial H\"older continuity of 
weak solutions to  \eqref{maineq}. 
We say $u\in W^{1,\Phi}(\Omega,\mathbb{R}^N)$ 
is a weak solution to \eqref{maineq} if it satisfies
$$
\int_{\Omega}\mathbf{a}(x,u,Du):D\varphi\, dx=0\quad \text{for every }
\varphi\in W^{1,\Phi}_0(\Omega,\mathbb{R}^N),
$$
where $W^{1,\Phi}_0(\Omega,\mathbb{R}^N)$ is the closure of $C_0^\infty(\Omega,\mathbb{R}^N)$ 
in $W^{1,\Phi}(\Omega,\mathbb{R}^N)$. To obtain the desired regularity, we impose 
continuity assumptions on the nonlinearity $\mathbf{a}$ as follows. 
There exists a nondecreasing and concave function $\mu:[0,\infty)\to[0,\infty)$ 
with $\mu(0)=0$ such that
\begin{gather}\label{zetaconti}
|\mathbf{a}(x,\zeta_1,\xi)-\mathbf{a}(x,\zeta_2,\xi)  |
\leq \Lambda\, \mu(|\zeta_1-\zeta_2|)\, \Phi_{p(x)-1}(1+|\xi|), \\
\label{xiconti}
|\partial\mathbf{a}(x,\zeta,\xi_1)-\partial\mathbf{a}(x,\zeta,\xi_2)|
\leq \Lambda \mu \Big(\frac{|\xi_1-\xi_2|}{1+|\xi_1|+|\xi_2|}\Big) 
\Phi_{p(x)-2}(1+|\xi|).
\end{gather}
Note that the prototype of $\mathbf{a}$ is
\begin{equation}\label{proto}
\mathbf{a}(x,\zeta,\xi)\equiv \mathbf{a}(x,\xi)=D_{\xi}[\Phi(x,1+|\xi|)],
\end{equation}
for which one can see that the nonlinearity $\mathbf{a}$ in \eqref{proto} 
satisfies \eqref{xiconti} with $\mu(r)=r^\alpha$ for any $\alpha\in(0,1)$.

The following theorem is the main result in this article; 
while the notation used here will be introduced in the next section.

\begin{theorem}\label{mainthm}
Let $u\in W^{1,\Phi}(\Omega,\mathbb{R}^N)$ be a weak solution to \eqref{maineq}. 
Then there exists $\Omega_u\subset \Omega$ such that $|\Omega\setminus\Omega_u|=0$ 
and
$u\in C^\alpha(\Omega_u,\mathbb{R}^N)$ for every $\alpha\in(0,1)$. Moreover,
\begin{align*}
\Omega\setminus\Omega_u
&\subset \Big\{x_0\in\Omega: \liminf_{r\downarrow0} 
 -\hspace{-0.38cm}\int_{B_r(x_0)}|Du-(Du)_{x_0,r}|\,dx>0\Big\}\\
&\cup \Big\{x_0\in\Omega: \limsup_{r\downarrow0}
 -\hspace{-0.38cm}\int_{B_r(x_0)}\Phi(x, |Du|)\,dx=\infty\Big\}.
\end{align*}
\end{theorem}

Partial H\"older continuity for elliptic systems is one of classical regularity 
issues, which was first presented by Campanato \cite{Cam1,Cam2}, after 
several papers where partial regularity has been proved for the gradient 
(rather than for the solution itself) under stronger assumptions. 
Note that there has been reported various systems, for which there exists a 
weak solution having a singularity, see for example \cite{Min1}, 
to which we also refer for results concerning regularity and vectorial problems.
 For systems with standard $p$-growth, Foss and Mingione \cite{FM1} proved it
 when $p\geq2$ and $\mathbf{a}$ is continuous for the space variable $x$, 
see also \cite{BFM1} for the parabolic counterpart. On the other hand, 
Habermann \cite{Ha1} considered systems with $p(x)$-growth and also $p$-growth
 with $1<p<2$. In addition, B\"ogelein, Duzaar, Harbermann and Scheven \cite{BDHS1} 
considered $p$-growth systems with VMO coefficients. Finally, it is worth 
mentioning the papers \cite{KM1,KM2}, which treat elliptic system with measure data.

As mentioned at the  beginning the problems considered here are significant 
since they do combine features of problems with non-standard growth typical 
of the Orlicz setting (see for instance \cite{Cia1, Cia2}) with those coming 
from variable exponent spaces (see \cite{DHHR1} and references) and more in 
general of functionals with non-standard growth conditions of non-autonomous type. 
The main point of the energies considered here is that the type of growth with 
respect to the gradient variable strongly depends on the variable $x$; this type 
of phenomenon poses new problems, leads to new techniques and conditions for 
regularity. This is a recent direction that only now starts being exploited 
in detail \cite{BCM1,BCM2,ColM1,ColM2,ColM3,EMM1, EMM2, HHK1,HHT1,Ok1,Ok2,Ok3,Ok4}.
The main common feature of all these examples is the analysis of functionals
\begin{equation}\label{integra}
v \mapsto \int_{\Omega}f(x, Dv), dx
\end{equation}
with $\Psi$-growth, i.e.,
$$
f(x, Dv)\approx \Psi(x, |Dv|)
$$
and for every fixed $x$, $z \mapsto \Psi(x, |z|)$ is an Orlicz function.
 In many cases, this maps has non-standard growth conditions of polynomial type, 
as for instance considered in \cite{M2, M3, Le1, Le2}.

In this paper, we consider systems with $\Phi$-growth, with $\Phi$ being 
the function introduced in \eqref{phi}; anyway, as already mentioned at 
the beginning, some of the methods considered can be extended to more general 
choices of $\Phi$. These systems are related to functionals as \eqref{integra} 
via the Euler-Lagrange equation. For the single equation case, i.e. $N=1$, 
and $\mathbf{a}(x,\zeta,\xi)\equiv \mathbf{a}(x,\xi)$, in view of \cite{GP1,Ok2}, one can 
deduce that $u\in C^{0,\alpha}$ for all $\alpha\in(0,1)$ if $p(\cdot)$ 
satisfies \eqref{vlogholder} without any further continuity assumption 
like \eqref{xiconti}. In addition, For systems with $p(x)$-growth, 
partial H\"older continuity was proved by Harbermann \cite{Ha1}, 
in which the assumptions on $p(\cdot)$ and $\mathbf{a}$ are the exactly same to ones 
in this paper (of course, those are modified in the setting of $p(x)$-growth).

However, we point out that compared with \cite{FM1,BDHS1,Ha1}, our case has a 
different behavior. Note that the function $\Phi$ is vary with respect to the 
position $x$. Moreover, even on fixed positions, we have 
$\Phi(x_0,t)=\Phi_{p(x_0)}(t)=t^{p(x_0)}\log(e+t)$, which does not  have 
standard $p$-growth structure. In fact, $\Phi_p$, $1<p<\infty$, is the typical 
example of  the Orlicz function. We refer to \cite{Li1,DE1} and related 
references for problems with Orlicz growth. This makes obtaining partial 
H\"older continuity difficult. We point out that although the paper \cite{Ha1} 
considers systems with variable growth, estimates in there are obtained in 
space with fixed exponent, i.e., classical Lebesgue space.

The first difficulty is the lack of normalization property, and the second 
one is that we do not expect the integral version of H\"older's continuity. 
If $G(t)=t^p$, then it is trivial that $G(st)=G(s)G(t)$ for every $s,t>0$, and
$$
\int_U|fg|\, dx
=\Big(\int_U G(|f|)\, dx\Big)^{1/p}\Big(\int_U |g|^{p'}\, dx\Big)^{1/p'}.
$$
Note that these properties have been used over all in \cite{FM1,BDHS1,Ha1}. 
However, for $G\equiv \Phi_p$, those properties do not hold anymore. 
In fact, they are not true for general Orlicz function. It seems that those 
difficulties are the reasons why the partial H\"older continuity for Orlicz 
growth has not  been proved yet. Therefore, in our best knowledge, our result 
with $p(\cdot)\equiv p$  is the first result for the partial H\"older continuity 
for systems with Orlicz growth.

The technical starting point of our proof is the approach in \cite{FM1}.
We also refer to \cite{BDHS1,Ha1}. However, we need a lot of new ideas to
overcome the difficulties mentioned above.  We also point out that in the proof 
of \cite{Ha1}, when analyzing on a ball, the author froze the variable 
exponent $p(\cdot)$ as a maximum value in a ball. This makes computations 
complicated since we consider the sequence of shrinking balls hence the frozen 
exponent is continuously changed. On the other hand, in this paper we freeze 
the variable exponent as a center point of balls. Therefore, the exponent 
is not changed with respect to shrinking balls. This makes our argument 
somehow clearer.

We organize this paper as follows. In the next section, we introduce notation
and ingredients of the proof of Theorem \ref{maineq} such as higher integrability, 
affine functions and the $\mathcal{A}$-harmonic approximation lemma. 
In Section \ref{sec3}, we prove Theorem \ref{maineq}. To do that we first 
derive a Caccioppoli type inequality, compare weak solutions with 
$\mathcal{A}$-harmonic maps, investigate a decay of a suitable excess functional 
$E$, see \eqref{E}, and finally obtain estimates of $Du$ in the Morrey space
 which implies the desired H\"older continuity.


\section{Preliminaries}\label{sec2}

\subsection{Orlicz function $\Phi_p$ and modulus of continuity}
Let $y\in\mathbb{R}^n$, $r>0$ and $U$ is a bounded domain in $\mathbb{R}^n$.
$B_r(y)$ is a standard ball in $\mathbb{R}^n$ centered at $y$ with radius $r$. 
For a locally integrable vector valued function $f$ in $\mathbb{R}^n$,  
$(f)_{y,r}$ is the integral average of $f$ in $B_r(y)$ such that
$$
(f)_{y,r}=-\hspace{-0.38cm}\int_{B_r(y)}f\, dx=\frac{1}{|B_r(y)|}\int_{B_r(y)}f\, dx.
$$
Note that if the center point is not important or well understood, we will use
the brief notation $B_r=B_r(y)$ and $(f)_{r}=(f)_{y,r}$. 
For matrix values $A=(a_{ij}),B=(b_{ij})\in \mathbb{R}^{nN}$, we define an inner product 
by $A:B:=\sum a_{ij} b_{ij}$.  $P:\mathbb{R}^n\to\mathbb{R}^N$ is always an affine function;
that is, $P(x)=Ax+b$ for some $A\in \mathbb{R}^{nN}$ and $b\in \mathbb{R}^N$. 
For simplicity, we shall write $\log^\beta t:= [\log t]^\beta$, where 
$\beta>0$ and $t\geq 1$.

Suppose $1<p<\infty$, and recall the function $\Phi_p(t)=t^p\log(e+t)$.
We define the conjugate function of $\Phi_p$, $\Phi^*_p:[0,\infty)\to[0,\infty)$, by
$$
\Phi^*_p(\tau):=\sup_{t\geq0} (\tau t-\Phi_p(t)).
$$
We then state elementary properties of $\Phi_p$ and $\Phi^*_p$. 
We refer to \cite{Ok2} for their proofs, see also \cite{DHHR1}.

\begin{proposition}\label{propbasic11}
 Let $t,\tau>0$, $1<p<\infty$, $s>1$ and $0<\theta<1$.
\begin{itemize}
\item[(1)] $\Phi_p(s t)\leq s^{p+1}\Phi_p(t)$ and
 $\Phi_p(\theta t)\leq \theta^p\Phi_p(t)$.

\item[(2)] $\Phi^*_p(s\tau)\leq s^{\frac{p}{p-1}}\Phi^*_p(\tau)$ and 
$\Phi^*_p(\theta \tau)\leq \theta^{\frac{p+1}{p}}\Phi^*_p(\tau)$.

\item[(3)] $\Phi_p(t+\tau)\leq \frac{1}{2}(\Phi_p(2t)
+\Phi_p(2\tau))\leq 2^{p}(\Phi_p(t)+\Phi_p(\tau))$.

\item[(4)] (Young's inequality) For any $\kappa\in(0,1]$, we have
\begin{gather}\label{youngs1}
t\tau\leq \Phi_p(\kappa^{\frac{1}{p}}t)+\Phi^*_p(\kappa^{-\frac{1}{p}}\tau)
\leq \kappa \Phi_p(t)+\kappa^{-\frac{1}{p-1}}\Phi^*_p(\tau), \\
\label{youngs}
t\tau \leq \Phi^*_p(\kappa^{\frac{p-1}{p}} \tau)+\Phi_p(\kappa^{-\frac{p-1}{p}}t)
\leq \kappa \Phi^*_p(\tau)+\kappa^{-\frac{p^2-1}{p}}\Phi_p(t).
\end{gather}

\item[(5)] There exists $c=c(p)>1$ such that
\begin{equation}\label{prop111}
\frac1c \Phi_p(t)\leq \Phi^*_p\big(  \Phi_p(t)t^{-1} \big) \leq c \Phi_p(t).
\end{equation}
Moreover, if $1<\gamma_1\leq p \leq \gamma_2<\infty$, then the constant $c$ 
depends only on $\gamma_1$ and $\gamma_2$, instead of $p$.
\end{itemize}
\end{proposition}

Note that (1)--(3) of the previous proposition are simple results, hence, 
we will use them without any comment except the cases that they are crucially used.

We recall the concave functions $\omega(\cdot)$ and $\mu(\cdot)$ introduced in 
Section \ref{sec1}. Then, since $\omega(0)=\mu(0)=0$, we see that for $r>0$, 
$0<\theta<1$ and $s\geq1$,
\begin{equation}\label{furtherass1}
\theta\omega(r)\leq \omega(\theta r) ,\quad
\theta \mu(r)\leq \mu(\theta r),\quad  \omega(s r)\leq s\omega(r),\quad
 \mu(s r)\leq s \mu(r).
\end{equation}
Moreover, we will ultimately consider sufficiently small value $r$ in 
$\omega(r)\log(1/r)$
and $\mu(r)$ in this paper. Hence, in view of  \eqref{vlogholder} and the 
monotonicity of $\mu$, we shall assume without loss of generality that
\begin{equation}\label{furtherass2}
\omega(r)\log\big(\frac{1}{r}\big)\leq 1\quad \text{and}\quad \mu(r)\leq 1\quad
\text{for every } r>0.
\end{equation}

We fix a weak solution $u\in W^{1,\Phi}(\Omega,\mathbb{R}^N)$ to \eqref{maineq}, 
and then define
\begin{equation}\label{M}
M:=\int_{\Omega}[\Phi(x,Du)+1]\, dx+1.
\end{equation}

For the rest of this article, we write letter $c>0$ for any constant depending
on the structure constants $n,N,\gamma_1,\gamma_2,\nu,\Lambda$.

\subsection{Basic inequalities}

First, we recall Jensen's inequality such that for any convex function 
$G:[0,\infty)\to[0,\infty)$,
$$
G\Big(-\hspace{-0.38cm}\int_{U} |f|\, dx \Big)\leq -\hspace{-0.38cm}\int_UG(|f|)\, dx,
$$
which will be used frequently in this paper.
For any vector valued function $f$, constant vector $A$ and $1<p<\infty$, we have
\begin{gather*}
-\hspace{-0.38cm}\int_{B_r(y)}|f-(f)_{y,r}|^2\,dx\leq 2-\hspace{-0.38cm}\int_{B_r(y)}|f-A|^2\,dx, \\
-\hspace{-0.38cm}\int_{B_r(y)}\Phi_p(|f-(f)_{y,r}|)\,dx\leq 2^{p+1}-\hspace{-0.38cm}\int_{B_r(y)}\Phi_p(|f-A|)\,dx.
\end{gather*}
Note that the second inequality follows from (3) of Proposition \ref{propbasic11}
and Jensen's inequality.

The next inequality related to the embedding property from $L\log^\beta L(U)$
space to $L^q(U)$ space, where $\beta>0$ and $1<q<\infty$. We refer to \cite{IV1},
 see also \cite{AM1}.
\begin{equation}\label{LlogL}
-\hspace{-0.38cm}\int_{U} |f|\log^\beta\Big(e+\frac{|f|}{-\hspace{-0.3cm}\int_U|f|\,dx}\big)\, dx
\leq c(q,\beta)\Big(-\hspace{-0.38cm}\int_{U}|f|^q\, dx \Big)^{1/q}.
\end{equation}
Here, the constant $c(q,\beta)>0$ is continuously changed with respect to $q$ 
and $\beta$.
The next lemma is the Sobolev-Poincar\'e type inequality for the function $\Phi_p$,
 which can be easily obtained by using the result in \cite[Theorem 7]{DE1}.

\begin{lemma}\label{sobolevpoincare}
Let $1<\gamma_1\leq p\leq \gamma_2< \infty$ and $f\in W^{1,1}(B_r,\mathbb{R}^N)$. 
Then there exists $\sigma_s=\sigma_s(n,N,\gamma_1,\gamma_2)>0$ such that
\begin{equation}\label{sobopoin}
-\hspace{-0.38cm}\int_{B_r}\Big[\Phi_p\Big(\frac{|f-(f)_r|}{r}\Big)\Big]^{1+\sigma_s}\, dx
\leq c\Big(-\hspace{-0.38cm}\int_{B_r}\Phi_p(|Df|)\, dx\Big)^{1+\sigma_s}.
\end{equation}
for some $c=c(n,N,\gamma_1,\gamma_2)>0$.
\end{lemma}

We finally state the  iteration lemma, see \cite[Lemma 7.3]{Gi1} and
 \cite[Lemma 2.3]{FM1}.

\begin{lemma} \label{lemtech}
Let $\phi:(0,\rho]\to \mathbb{R}$ be a positive and nondecreasing function satisfying
$$
\phi(\theta^{k+1}\rho)\leq \theta^{q}\phi(\theta^k\rho)+c_0 (\theta^k\rho)\quad
 \text{for  }k=0,1,2,\dots,
$$
where $\theta\in(0,1)$, $q\in(0,n)$ and $c_0>0$. Then there exists
$\tilde c=\tilde c(n,\theta,\gamma)>0$ such that
$$
\phi(t)\leq \tilde c\big\{(\frac{t}{\rho})^q\phi(\rho)+c_0 t^q\big\}\quad
 \text{for } t\in(0,\rho).
$$
\end{lemma}


\subsection{Higher integrability}

We introduce one of the basic regularity properties of weak solutions to \eqref{maineq},
which is the higher integrability.  For the scalar case, i.e. $N=1$, with
$$
\mathbf{a}(x,\xi)= D_\xi (\Phi(x,|\xi|)),
$$
this property was proved in \cite[Theorem 3.9]{Ok1}. 
We point out that one can also prove the same result for the system \eqref{maineq} 
in the same way. Therefore, we shall omit to prove the following result here.

\begin{theorem}\label{thmhigher}
Suppose that $p(\cdot):\Omega\to\mathbb{R}$ with \eqref{pxgamma} is vanishing 
log-H\"older continuous and satisfies the first inequality 
\eqref{furtherass2}, $\mathbf{a}:\Omega\times\mathbb{R}^N\times\mathbb{R}^{nN}\to\mathbb{R}^{nN}$ 
only satisfies the growth and ellipticity conditions in \eqref{acondition}, 
and $u\in W^{1,\Phi}(\Omega,\mathbb{R}^N)$ is a weak solution to \eqref{maineq}. 
Then, there exist small constants 
$\delta_0=\delta_0(n,N,\gamma_1,\gamma_2,\nu,\Lambda,\omega(\cdot))>0$ and 
$\sigma_h=\sigma_h(n,N,\gamma_1,\gamma_2,\nu,\Lambda)>0$ such that for any 
$B_{2r}(x_0)\subset \Omega$ with $r\leq \delta_0M^{-1}$, where $M>1$ is 
denoted in \eqref{M}, and any $\sigma\leq \sigma_h$ we have
\begin{equation}\label{higher1}
-\hspace{-0.38cm}\int_{B_r(x_0)}[\Phi(x,|Du|)]^{1+\sigma}\, dx
\leq c\Big(-\hspace{-0.38cm}\int_{B_{2r}(x_0)}\Phi(x,|Du|)\, dx+1\Big)^{1+\sigma}.
\end{equation}
for some $c=c(n,N,\gamma_1,\gamma_2,\nu,\Lambda)\geq 1$.
Moreover, for any $\delta\in(0,1]$ we also have
\begin{equation}\label{higher2}
-\hspace{-0.38cm}\int_{B_r(x_0)}[\Phi(x,|Du|)]^{1+\sigma}\, dx
\leq c(\delta)\Big\{\Big(-\hspace{-0.38cm}\int_{B_{2r}(x_0)}[\Phi(x,|Du|)]^\delta\, dx
\Big)^\frac{1+\sigma}{\delta}+1\Big\}.
\end{equation}
for some $c(\delta)=c(n,N,\gamma_1,\gamma_2,\nu,\Lambda,\delta)\geq 1$.
\end{theorem}

Note that estimates \eqref{higher2} follows from \eqref{higher1} by using 
the argument in \cite[Remark 6.12]{Gi1}.

From the above lemma we can observe the following.
 Consider $B_{4r}(x_0)\subset\Omega$ with $x_0\in\Omega$ and $r>0$ satisfying
\begin{equation}\label{rho0}
2r\leq \delta_0M^{-1}\quad \text{and} \quad
 \omega(2r)\leq \min\{\frac{\sigma_h}{8},\frac{1}{2}\},
\end{equation}
and set
$$
p_0:=p(x_0), \quad p_1:=\inf_{B_{4r}(x_0)}p(\cdot), \quad
p_2:=\sup_{B_{4r}(x_0)}p(\cdot).
$$
We first note from the fact $r\leq M^{-1}$ and the first inequality in 
\eqref{furtherass2} that
\begin{equation}\label{higher5}
\Big(-\hspace{-0.38cm}\int_{B_{2r}(x_0)}\Phi(x,|Du|)\, dx+1\Big)^{\omega(r)}
\leq c (r^{-n}M)^{\omega(r)}\leq cr^{-(n+1)\omega(r)}\leq c.
\end{equation}
Then since $p_0\leq p(x)(1+\omega(r))\leq p(x)(1+\sigma_h/8)$ for $x\in B_r(x_0)$,
applying \eqref{higher1} with $\sigma=\omega(r)$ and \eqref{higher5},  we have
\begin{align*}
-\hspace{-0.38cm}\int_{B_r(x_0)}\Phi_{p_0}(|Du|)\, dx
&\leq -\hspace{-0.38cm}\int_{B_{2r}(x_0)}[\Phi(x,|Du|)]^{1+\omega(r)}\, dx+1\\
&\leq c\Big(-\hspace{-0.38cm}\int_{B_{2r}(x_0)}\Phi(x,|Du|)\, dx+1\Big)^{1+\omega(r)}\\
&\leq c\Big(-\hspace{-0.38cm}\int_{B_{2r}(x_0)}\Phi(x,|Du|)\, dx+1\Big).
\end{align*}
Note that this implies $\Phi_{p_0}(|Du|) \in L^1_{loc}(\Omega)$.
In the same way, since 
$p_0(1+\sigma_h/2)\leq p(x)(1+\omega(r))(1+\sigma_h/2)\leq p(x)(1+\sigma_h)$ 
for $x\in B_r(x_0)$, we also have
\begin{equation}\label{higher3}
-\hspace{-0.38cm}\int_{B_r(x_0)}[\Phi_{p_0}(|Du|)]^{1+\frac{\sigma_h}{2}}\, dx
\leq c\Big\{\Big(-\hspace{-0.38cm}\int_{B_{2r}(x_0)}\Phi(x,|Du|)\, dx
\Big)^{1+\frac{\sigma_h}{2}}+1\Big\}.
\end{equation}
On the other hand, by taking $\delta=1/2\leq p_1/p_2$
 from \eqref{higher2}, H\"older's inequality and \eqref{higher5}
 with \eqref{furtherass1}
 we have
\begin{equation} \label{higher4}
\begin{aligned}
&-\hspace{-0.38cm}\int_{B_{2r}(x_0)}\Phi(x,|Du|)\, dx \\
&\leq c\Big\{\Big(-\hspace{-0.38cm}\int_{B_{4r}(x_0)}[\Phi(x,|Du|)]^{\frac{p_1}{p_2}}\, dx\Big)^{p_2/p_1}+1\Big\} \\
&\leq c\Big(-\hspace{-0.38cm}\int_{B_{4r}(x_0)}\Phi_{p_0}(|Du|)\,dx+1\Big)
\Big(-\hspace{-0.38cm}\int_{B_{4r}(x_0)}\Phi(x,|Du|)\, dx+1\Big)^{\omega(8r)}\\
&\leq c\Big(-\hspace{-0.38cm}\int_{B_{4r}(x_0)}\Phi_{p_0}(|Du|)\, dx+1\Big).
\end{aligned}
\end{equation}
Moreover, since
$$
\frac{p_2-1}{p_0-1}\Big(1+\frac{\sigma_h}{4}\Big)p_0
\leq p(x)(1+\omega(r))^2\Big(1+\frac{\sigma_h}{4}\Big)
\leq p(x)\Big(1+\frac{7}{8}\sigma_h\Big),
$$
applying \eqref{higher1} with $\sigma=\frac{p_2-1}{p_0-1}(1+\omega(r))$, or
\[
\sigma=\frac{p_2-1}{p_0-1}\big(1+\frac{\sigma_h}{4}\big)(1+\omega(r)),
\]
 \eqref{higher5} and \eqref{higher4}, we have
\begin{equation} \label{high1}
\begin{aligned}
-\hspace{-0.38cm}\int_{B_r} \Phi_{p_0}(1+|Du|)^{\frac{p_2-1}{p_0-1}} \, dx 
&\leq c\Big(-\hspace{-0.38cm}\int_{B_{2r}} \Phi(x,|Du|) \, dx+1\Big)^{(1+\omega(r))^2}\\
&\leq c\Big( -\hspace{-0.38cm}\int_{B_{4r}} \Phi_{p_0}(|Du|) \, dx+1\Big)
\end{aligned}
\end{equation}
and
\begin{equation}\label{high2}
-\hspace{-0.38cm}\int_{B_r} \Phi_{p_0}(1+|Du|)^{\frac{p_2-1}{p_0-1}\big(1+\frac{\sigma_h}{4}\big)}\, dx
\leq c\Big(-\hspace{-0.38cm}\int_{B_{4r}} \Phi_{p_0}(|Du|) \, dx+1\Big)^{1+\frac{\sigma_h}{4}}.
\end{equation}


\subsection{Affine functions} \label{affine}
For a given $u\in L^2(B_r(x_0),\mathbb{R}^N)$, we define an affine function 
$P_{x_0,r}$ by the minimizer of the functional
$$
P\mapsto -\hspace{-0.38cm}\int_{B_r(x_0)}|u-P|^2\,dx.
$$
Then we can easily see that
$P_{x_0,r}=DP_{x_0,r} (x-x_0)+(u)_{x_0,r}$, where
\[
 DP_{x_0,r}:= \frac{n+2}{r^2}-\hspace{-0.38cm}\int_{B_r(x_0)}u\otimes (x-x_0)\, dx.
\]
Then we have the following lemma.

\begin{lemma} \label{lemDP}
Let $2< p<\infty$ and $r>0$.
\begin{itemize}
\item[(1)] For any $u\in L^{\Phi_p}(B_r(x_0))$ and $\theta\in(0,1)$, we have
\begin{equation}\label{DPDP}
\Phi_p(|DP_{x_0,r}-DP_{x_0,\theta r}|)
\leq c -\hspace{-0.38cm}\int_{B_{\theta r}(x_0)}\Phi_p\Big(\frac{|u-P_{x_0, r}|}{\theta r}\Big)\,dx.
\end{equation}

\item[(2)] For any $u\in W^{1,\Phi_p}(B_ r(x_0))$, we have
\begin{equation}\label{DPDu}
\Phi_p\big(|DP_{x_0, r}-(Du)_{B_ r(x_0)}|\big)
\leq c -\hspace{-0.38cm}\int_{B_{ r}(x_0)}\Phi_p(|Du-(Du)_{B_{ r}(x_0)}|)\,dx.
\end{equation}
\end{itemize}
Here, the constants $c$ depend only on $n,N,p$.
\end{lemma}

\begin{proof} By \cite[Lemma 2]{Kr1}, we have
\begin{gather*}
|DP_{x_0, r}-DP_{x_0,\theta r}|^2 
 \leq c -\hspace{-0.38cm}\int_{B_{\theta r}(x_0)}\frac{|u-P_{x_0, r}|^2}{(\theta r)^2}\,dx, \\
|DP_{x_0, r}-(Du)_{B_ r(x_0)}|^2\leq c -\hspace{-0.38cm}\int_{B_{ r}(x_0)}|Du-(Du)_{B_ r(x_0)}|^2\,dx.
\end{gather*}
Using these and Jensen's inequality for the convex map $t\mapsto \Phi_{p}(\sqrt t)$, 
we obtain
\begin{align*}
\Phi_p(|DP_{x_0, r}-DP_{x_0,\theta r}|)
&\leq c \Phi_{p}\Big(-\hspace{-0.38cm}\int_{B_{\theta r}(x_0)}
 \frac{|u-P_{x_0, r}|^2}{(\theta r)^2}\,dx\Big)^{1/2}\\
&\leq c -\hspace{-0.38cm}\int_{B_{\theta r}(x_0)}
\Phi_p\Big(\frac{|u-P_{x_0, r}|}{\theta r}\Big)\,dx,
\end{align*}
which yileds \eqref{DPDP}. In the same way we obtain \eqref{DPDu}.
\end{proof}


\subsection{$\mathcal{A}$-harmonic approximation}
Let $\mathcal A$ be a bilinear form in $\mathbb{R}^{nN}$. Moreover $\mathcal A$ 
satisfies the so-called Legendre-Hadamard condition: there exist $0<\nu\leq\lambda$ 
such that
\begin{equation}\label{lhcondition}
\nu|\xi|^2|\eta|^2\leq \mathcal \xi\otimes\eta 
: \xi\otimes\eta\leq \Lambda|\xi|^2 |\eta|^2
\end{equation}
for every $\xi\in\mathbb{R}^n$, $\eta\in\mathbb{R}^N$. 
For this $\mathcal A$ we say that
$h\in W^{1,2t}(\Omega,\mathbb{R}^N)$ is $\mathcal A$-harmonic if
$$
\int_\Omega \mathcal A Dh:D\varphi=0
$$
for every $\varphi\in C_0^1(\Omega,\mathbb{R}^N)$. We then state the so-called 
$\mathcal A$-harmonic approximation lemma, see \cite[Lemma 3.3]{DS1}.

\begin{lemma} \label{lemharmonicapprox}
Let $\epsilon>0$ and $0<\nu\leq \Lambda$. Assume that there exists small 
$\delta=\delta(n,N,\nu,\Lambda,\epsilon)>0$ such that if the bilinear form 
$\mathcal A$ on $\mathbb{R}^{nN}$ satisfies \eqref{lhcondition}, $r>0$ and 
$w\in W^{1,2}(B_r,\mathbb{R}^N)$ satisfies 
\begin{gather*}
-\hspace{-0.38cm}\int_{B_r}|Dw|^2\, dx\leq 1, \\
\big|-\hspace{-0.38cm}\int_{B_\rho}\mathcal A Dw:D\varphi\, dx\big|
\leq \delta \|D\varphi\|_{L^\infty(B_r)}\quad \text{for every }
 \varphi\in C^1_0(B_r,\mathbb{R}^N),
\end{gather*}
then there exists $\mathcal A$-harmonic map $h$ such that
$$
-\hspace{-0.38cm}\int_{B_r}|Dw|^2\, dx\leq 1, \quad  r^{-2}-\hspace{-0.38cm}\int_{B_r}|w-h|^2\,dx\leq \epsilon.
$$
\end{lemma}


\section{Partial continuity}\label{sec3}

Let us first consider $\rho>0$  satisfying
\begin{equation}\label{rho00}
\rho\leq\rho_0:=\sup\big\{r >0: 4r\leq \delta_0M^{-1},\;
 \omega(4r)\leq  \min\{\frac{\sigma_h}{8},\frac{1}{2}\} \big\},
\end{equation}
where $M$ is denoted in \eqref{M} and $\delta_0,\sigma_h$ are determined 
in Theorem \ref{thmhigher}.  Then we see that $\rho$ satisfies \eqref{rho0} with 
$r$ replaced by $2\rho$.
We then consider $B_{4\rho}(x_0)\subset\Omega$ with $x_0\in\Omega$ and 
$\rho\leq \rho_0$ and set
$$
p_0:=p(x_0),\quad p_1:= \inf_{B_{2\rho}(x_0)} p(\cdot),\quad 
 p_2:= \sup_{B_{2\rho}(x_0)} p(\cdot),
$$
\begin{equation}\label{C}
C(x_0,\rho,P):=-\hspace{-0.38cm}\int_{B_\rho(x_0)}\big[\frac{|Du-DP|^2}{(1+|DP|)^2}
+\frac{\Phi_{p_0}(|Du-DP|)}{\Phi_{p_0}(1+|DP|)}\big]\,dx
\end{equation}
and
\begin{equation} \label{E}
 E(x_0,\rho,P):=  C(x_0,\rho,P)
+\Big[\mu\Big(-\hspace{-0.38cm}\int_{B_{\rho}} |u-P(x_0)|^2\, dx\Big)
+\omega(\rho)\log(\frac{1}{\rho}) \Big]^{\frac{1}{2p_0-1}}.
\end{equation}
Here, $P:\mathbb{R}^n\to\mathbb{R}^N$ is any affine function.
In this setting we derive a Caccioppoli type inequality.

\begin{lemma}\label{lemcaccio}
Assume that
\begin{equation}\label{lem1ass}
C(x_0,8\rho,P)\leq 1.
\end{equation}
Then 
\begin{equation} \label{caccio}
\begin{aligned}
 C(x_0,\rho,P)
&\leq c-\hspace{-0.38cm}\int_{B_{2\rho}(x_0)}
 \Big[\frac{|u-P|^2}{(2\rho)^2(1+|DP|)^2}\, dx
 +\frac{\Phi_{p_0}(|u-P|/(2\rho))}{\Phi_{p_0}(1+|DP|)}\Big]\,dx\\
&\quad +c\Big\{\mu\Big(-\hspace{-0.38cm}\int_{B_{2\rho}(x_0)} |u-P(x_0)|^2\, dx \Big)
 +\omega(2\rho)\log(\frac{1}{2\rho})\Big\}
\end{aligned}
\end{equation}
for some $c=c(n,N,\gamma_1,\gamma_2,\nu,\Lambda)>0$.
\end{lemma}

\begin{proof}
Set $\varphi=\eta^{p_0+1}(u-P)$, where $\eta \in C_0^\infty(B_{2\rho})$ 
with $0\leq \eta\leq 1$, $\eta\equiv 1$ on $B_\rho$ and 
$|D\eta|\leq c(n)/\rho$. Taking $\varphi$ as a test function in \eqref{maineq}, 
we have
$$
-\hspace{-0.38cm}\int_{B_{2\rho}}\eta^{p_0+1}\mathbf{a}(x,u,Du):D(u-P)\,dx
=-(p_0+1) -\hspace{-0.38cm}\int_{B_{2\rho}}\eta^{p_0} \mathbf{a}(x,u,Du):D\eta\otimes (u-P)\, dx.
$$
This and the trivial identity
$$
-\hspace{-0.38cm}\int_{B_{2\rho}}\mathbf{a}(x_0,P(x_0),DP):D\varphi\,dx=0
$$
imply 
\begin{equation}\label{pf101}
\begin{aligned}
 I_1&:=-\hspace{-0.38cm}\int_{B_{2\rho}}\eta^{p_0+1}(\mathbf{a}(x_0,u,Du)-\mathbf{a}(x_0,u,DP)):(Du-DP)\, dx\\
&=-\hspace{-0.38cm}\int_{B_{2\rho}}(\mathbf{a}(x_0,u,Du)-\mathbf{a}(x_0,u,DP)):D\varphi\, dx\\
&\quad-(p_0+1)-\hspace{-0.38cm}\int_{B_{2\rho}}\eta^{p_0}(\mathbf{a}(x_0,u,Du)-\mathbf{a}(x_0,u,DP)):D\eta \otimes (u-P)\, dx\\
&=-\hspace{-0.38cm}\int_{B_{2\rho}}(\mathbf{a}(x_0,u,Du)-\mathbf{a}(x,u,Du)):D\varphi\, dx\\
&\quad +-\hspace{-0.38cm}\int_{B_{2\rho}}(\mathbf{a}(x_0,P(x_0),DP)-\mathbf{a}(x_0,u,DP)):D\varphi\, dx\\
&\quad -(p_0+1)-\hspace{-0.38cm}\int_{B_{2\rho}}\eta^{p_0}(\mathbf{a}(x_0,u,Du)-\mathbf{a}(x_0,u,DP)):D\eta \otimes (u-P)\, dx\\
&=:I_2+I_3-I_4.
\end{aligned}
\end{equation}

Estimate for $I_1$. From \eqref{monotonicity} we have
\begin{equation}\label{I1}
-\hspace{-0.38cm}\int_{B_\rho}\eta^{p_0+1}\{\Phi_{p_0-2}(1+|DP|)|Du-DP|^2
+\Phi_{p_0}(|Du-DP|)\}\,dx\leq cI_1.
\end{equation}

Estimate for $I_2$. Using \eqref{xconti} and \eqref{youngs1}, for $\kappa\in(0,1)$,
we have
\begin{align*}
|I_2| 
&\leq  c -\hspace{-0.38cm}\int_{B_{2\rho}} \omega(2\rho)\log(e+1+|Du|)\,  \Phi_{p_2-1}(1+|Du|)
 \Big(|Du-DP|+\frac{|u-P|}{\rho}\Big)\, dx \\
&\leq \kappa -\hspace{-0.38cm}\int_{B_{2\rho}} 
\big[\Phi_{p_0}(|Du-DP|)+\Phi_{p_0}\big(\frac{|u-P|}{\rho}\big)\big]\, dx\\
&\quad + c(\kappa)-\hspace{-0.38cm}\int_{B_{2\rho}} \Phi_{p_0}^*
\big(\omega(2\rho)\log(e+1+|Du|) \Phi_{p_2-1}(1+|Du|)\big) \, dx.
\end{align*}
Now we estimate the last integral above:
\begin{equation}\label{tildeI2}
\tilde I_2:=-\hspace{-0.38cm}\int_{B_{2\rho}} \Phi_{p_0}^*
\left(\omega(2\rho)\log(e+1+|Du|) \Phi_{p_2-1}(1+|Du|)\right) \, dx.
\end{equation}
Using (2) and (5) of Proposition \ref{propbasic11} and that
 $1+|Du|\leq \Phi_{p_0}(1+|Du|)^{\frac{p_2-1}{p_0-1}}$, we have
\begin{align*}
\tilde I_2 
&\leq -\hspace{-0.38cm}\int_{B_{2\rho}} [\omega(2\rho)\log(e+1+|Du|)]^{i(p_0)}
 (1+|Du|)^{(p_2-p_0)\frac{p_0}{p_0-1}} \Phi_{p_0}(1+|Du|) \, dx \\
& \leq c\Big[\omega(2\rho)\log
\Big(e+\big(\Phi_{p_0}(1+|Du|)^{\frac{p_2-1}{p_0-1}}\big)_{2\rho}\Big)\Big]^{i(p_0)}
 -\hspace{-0.38cm}\int_{B_{2\rho}}[ \Phi_{p_0}(1+|Du|)]^{\frac{p_2-1}{p_0-1}}\, dx\\
&\quad + c [\omega(2\rho)]^{\frac{p_0+1}{p_0}} -\hspace{-0.38cm}\int_{B_{2\rho}} 
\log^{\frac{p_0}{p_0-1}}\Big(e+\frac{\Phi_{p_0}(1+|Du|)^{\frac{p_2-1}{p_0-1}}}
{(\Phi_{p_0}(1+|Du|)^{\frac{p_2-1}{p_0-1}})_{2\rho}}\Big) \\
&\quad\times \big[\Phi_{p_0}(1+|Du|)\big]^{\frac{p_2-1}{p_0-1}} dx,
\end{align*}
where
$$
i(p_0)=\begin{cases}
\frac{p_0}{p_0-1}& \text{if } \omega(2\rho)\log(e+1+|Du|)> 1,\\[4pt]
\frac{p_0+1}{p_0}& \text{if } \omega(2\rho)\log(e+1+|Du|)\leq 1.
\end{cases}
$$
Then, applying \eqref{furtherass2}, \eqref{LlogL} and two inequalities
 concerning the higher integrability: \eqref{high1} and \eqref{high2} with $r=2\rho$, 
we obtain
\begin{equation}\label{tildeI21}
\begin{aligned}
\tilde I_2 
&\leq  c\big[\omega(2\rho)\log(\frac{1}{\rho})\big]^{i(p_0)}
 -\hspace{-0.38cm}\int_{B_{8\rho}} \Phi_{p_0}(1+|Du|)\, dx\\
&\quad +[\omega(2\rho)]^{\frac{p_0+1}{p_0}}
\Big( -\hspace{-0.38cm}\int_{B_{2\rho}} \Phi_{p_0}(1+|Du|)^{\frac{p_2-1}{p_0-1}
\big(1+\frac{\sigma_h}{4}\big)} \, dx\Big)^{\frac{1}{1+\frac{\sigma_h}{4}}}\\
 & \leq c\big[\omega(2\rho)\log(\frac{1}{2\rho})\big]^{\frac{p_0+1}{p_0}}
 -\hspace{-0.38cm}\int_{B_{8\rho}} \Phi_{p_0}(1+|Du|)\, dx\\
& \leq  c\big[\omega(2\rho)\log(\frac{1}{2\rho})\big]^{\frac{p_0+1}{p_0}}
 \Phi_{p_0}(1+|DP|).
\end{aligned}
\end{equation}
In the last equality above, we have used the assumption \eqref{lem1ass} so that
\begin{equation}\label{341}
\begin{aligned}
 -\hspace{-0.38cm}\int_{B_{8\rho}} \Phi_{p_0}(1+|Du|)\, dx
&\leq c\Big(-\hspace{-0.38cm}\int_{B_{8\rho}} \Phi_{p_0}(|Du-DP|)\, dx
 +\Phi_{p_0}(1+|DP|)\, dx\Big) \\
&\leq c \Phi_{p_0}(1+|DP|).
\end{aligned}
\end{equation}
Therefore, 
\begin{equation}\label{I2}
\begin{aligned}
 |I_2| &\leq  \kappa -\hspace{-0.38cm}\int_{B_{2\rho}} \big[\Phi_{p_0}(|Du-DP|)
 +\Phi_{p_0}\big(\frac{|u-P|}{\rho}\big)\big]\, dx\\
&\quad +c(\kappa)\, \omega(2\rho)\log(\frac{1}{2\rho})\Phi_{p_0}(1+|DP|).
\end{aligned}
\end{equation}

Estimate for $I_3$. Applying \eqref{zetaconti},
$$
|I_3| \leq c -\hspace{-0.38cm}\int_{B_{2\rho}}\mu\left(|u-P(x_0)|^2\right)
\Phi_{p_0-1}(1+|DP|)\Big(\eta^{p_0+1}|Du-DP|+\frac{|u-P|}{\rho}\Big) \, dx.
$$
Using \eqref{youngs1},  \eqref{furtherass2} and Jensen's inequality, we see 
that for $\kappa\in(0,1)$,
\begin{equation}\label{I3}
\begin{aligned}
|I_3| 
&\leq \kappa -\hspace{-0.38cm}\int_{B_{2\rho}} 
 \Big[\eta^{p_0+1}\Phi_{p_0}(|Du-DP|)+\Phi_{p_0}
 \big(\frac{|u-P|}{\rho}\big)\Big] \, dx\\
&\quad + c(\kappa)-\hspace{-0.38cm}\int_{B_{2\rho}}\mu\left(|u-P(x_0)|^2\right)^{\frac{p_0+1}{p_0}}
 \Phi_{p_0}(1+|DP|)\, dx\\
&\leq  \kappa-\hspace{-0.38cm}\int_{B_{2\rho}}\eta^{p_0+1}\Phi_{p_0}(|Du-DP|)\, dx
 + -\hspace{-0.38cm}\int_{B_{2\rho}}\Phi_{p_0}\big(\frac{|u-P|}{\rho}\big)\, dx\\
&\quad + c(\kappa)\, \mu\Big(-\hspace{-0.38cm}\int_{B_{2\rho}} |u-P(x_0)|^2\, dx \Big)
\Phi_{p_0}(1+|DP|).
\end{aligned}
\end{equation}

Estimate for $I_4$. From the first inequality in \eqref{acondition}, 
Young's inequality and \eqref{youngs}, we have that for $\kappa\in(0,1)$,
\begin{equation}\label{I4}
\begin{aligned}
&|I_4|\\
&\leq c -\hspace{-0.38cm}\int_{B_{2\rho}}\Big(\int_0^1|\partial\mathbf{a}(x_0,u,t(Du-DP)+DP)|\,dt\Big)
 |Du-DP|\eta^{p_0}\frac{|u-P|}{\rho}\, dx\\
&\leq c -\hspace{-0.38cm}\int_{B_{2\rho}} \Phi_{p_0-2}(1+|DP|+|Du-DP|)
 |Du-DP|\eta^{p_0}\frac{|u-P|}{\rho}\, dx\\
&\leq  c \Phi_{p_0-2}(1+|DP|) -\hspace{-0.38cm}\int_{B_{2\rho}} \eta^{\frac{p_0+1}{2}} |Du-DP|\frac{|u-P|}{\rho}\, dx\\
&\quad +c -\hspace{-0.38cm}\int_{B_{2\rho}} \Phi_{p_0-1}(|Du-DP|)\eta^{p_0}\frac{|u-P|}{\rho}\, dx\\
&\leq \kappa-\hspace{-0.38cm}\int_{B_{2\rho}} \eta^{p_0+1} \left[ \Phi_{p_0-2}(1+|DP|)  |Du-DP|^2+\Phi_{p_0}(|Du-DP|)\right]\, dx\\
&\quad + c(\kappa)\Big(\Phi_{p_0-2}(1+|DP|) -\hspace{-0.38cm}\int_{B_{2\rho}} 
\frac{|u-P|^2}{\rho^2}\, dx
 + -\hspace{-0.38cm}\int_{B_{2\rho}} \Phi_{p_0}\big(\frac{|u-P|}{\rho}\big)\, dx\Big).
\end{aligned}
\end{equation}
Consequently, inserting \eqref{I1}-\eqref{I4} into \eqref{pf101} and choosing 
$\kappa$ sufficiently small, we get the estimate \eqref{caccio}.
\end{proof}


Next, we define
\begin{equation}\label{cala}
\mathcal{A}=\mathcal A(x_0,P)
:=\frac{\partial\mathbf{a}(x_0,P(x_0),DP)}{\Phi_{p(x_0)-2}(1+|DP|)},\quad
w:=\frac{u-P}{(1+|DP|)\sqrt{E(x_0,\rho,P)}}.
\end{equation}
Then one can see from \eqref{acondition}  that $\mathcal{A}$ 
satisfies \eqref{lhcondition}. The next lemma implies that if 
$E(x_0,\rho,P)$ is sufficiently small then one can apply Lemma 
\ref{lemharmonic} to $\mathcal A$ and $w$ denoted above.

\begin{lemma} \label{lemharmonic}
Suppose that \eqref{lem1ass} holds. Then there exists 
$c=c(n,N,\gamma_1,\gamma_2,\nu,\Lambda)>0$ such that
\begin{equation}\label{harmonicapprox}
\big|-\hspace{-0.38cm}\int_{B_\rho(x_0)} \mathcal ADw :D\varphi \,dx\big|
\leq c[\mu(\sqrt{E(x_0,\rho,P)})+E(x_0,\rho,P)]^{1/2}
 \sup_{B_\rho(x_0)} |D\varphi|
\end{equation}
for every $\varphi\in C_0^\infty(B_\rho(x_0))$.
\end{lemma}

\begin{proof}
We first consider $\varphi\in C_0^\infty(B_\rho(x_0))$ with 
$\sup_{B_\rho(x_0)} |D\varphi|= 1$, and set $v:= u-P$, i.e., 
$v= (1+|DP|)\sqrt{E(x_0,\rho,P)}\, w$, and $B_\rho=B_\rho(x_0)$. Then
\begin{equation}\label{pf201}
\begin{aligned}
&\Phi_{p_0-2}(1+|DP|)-\hspace{-0.38cm}\int_{B_\rho} \mathcal A Dv:D\varphi\, dx\\
&=-\hspace{-0.38cm}\int_{B_\rho}\int_0^1(\partial \mathbf{a}(x_0,P(x_0),DP)
 -\partial \mathbf{a}(x_0,P(x_0),DP+sDv))Dv:D\varphi\, dt\, dx\\
&\quad +-\hspace{-0.38cm}\int_{B_\rho}\int_0^1\partial \mathbf{a}(x_0,P(x_0),DP+sDv)Dv:D\varphi\, dt\, dx\\
&=:I_5+I_6.
\end{aligned}
\end{equation}
We first estimate $I_5$. Applying \eqref{xiconti}, we have
\begin{align*}
|I_5|
&\leq c-\hspace{-0.38cm}\int_{B_\rho}\int_0^1\mu\Big(\frac{s|Dv|}{1+|DP|}\Big)
 \Phi_{p_0-2}(1+|DP|+|DP+sDv|)|Dv|\, dt\, dx\\
& \leq c -\hspace{-0.38cm}\int_{B_\rho}\mu\Big(\frac{|Du-DP|}{1+|DP|}\Big)
 (\Phi_{p_0-2}(1+|DP|)+\Phi_{p_0-2}(|Du-DP|))|Du-DP|\, dx\\
&= c\Phi_{p_0-1}(1+|DP|)
 \times-\hspace{-0.38cm}\int_{B_\rho}\mu\Big(\frac{|Du-DP|}{1+|DP|}\Big) \\
&\quad\times \Big(\frac{|Du-DP|}{1+|DP|}+\frac{\Phi_{p_0-1}(|Du-DP|)}
{\Phi_{p_0-1}(1+|DP|)}\Big)\, dx.
\end{align*}
Set $B^+_{\rho}:= \{x\in B_\rho: |Du(x)-DP|>1+|DP| \}$ and 
$B^-_\rho:=B_\rho\setminus B^+_\rho$. Then, by \eqref{furtherass1} and 
\eqref{furtherass2}, in $B^+_\rho$,
\begin{align*}
&\mu\Big(\frac{|Du-DP|}{1+|DP|}\Big)
\Big(\frac{|Du-DP|}{1+|DP|}+\frac{\Phi_{p_0-1}(|Du-DP|)}{\Phi_{p_0-1}(1+|DP|)}
 \Big) \\
&\leq \frac{|Du-DP|^2}{(1+|DP|)^2}+\frac{\Phi_{p_0}(|Du-DP|)}{\Phi_{p_0}(1+|DP|)}.
\end{align*}
On the other hand, in $B^-_\rho$,
\begin{align*}
&\mu\Big(\frac{|Du-DP|}{1+|DP|}\Big)
 \Big(\frac{|Du-DP|}{1+|DP|}+\frac{\Phi_{p_0-1}(|Du-DP|)}{\Phi_{p_0-1}(1+|DP|)}\Big)
 \\
&\leq \mu\Big(\frac{|Du-DP|}{1+|DP|}\Big)\frac{|Du-DP|}{1+|DP|}
 +\mu\Big(\frac{|Du-DP|}{1+|DP|}\Big)
 \Big(\frac{\Phi_{p_0}(|Du-DP|)}{\Phi_{p_0}(1+|DP|)}\Big)^{\frac{p_0-1}{p_0}}.
\end{align*}
Using H\"older's inequality, \eqref{furtherass1} and Jensen's inequality for 
the concave function $\mu^{-1}$, we have from the pervious two inequalities that
\begin{align*}
&\frac{|I_5|}{\Phi_{p_0-1}(1+|DP|)} \\
&\leq c -\hspace{-0.38cm}\int_{B_\rho}\Big(\frac{|Du-DP|}{1+|DP|}+\frac{\Phi_{p_0}(|Du-DP|)}
{\Phi_{p_0}(1+|DP|)}\Big)\, dx\\
&\quad +c\Big[\mu\Big(-\hspace{-0.38cm}\int_{B_\rho}\frac{|Du-DP|}{1+|DP|}\,dx\Big)\Big]^{1/2}
 \Big(-\hspace{-0.38cm}\int_{B_\rho}\frac{|Du-DP|^2}{1+|DP|^2}\,dx\Big)^{1/2}\\
&\quad +c\Big[\mu\Big(-\hspace{-0.38cm}\int_{B_\rho}\frac{|Du-DP|}{1+|DP|}\,dx\Big)
\Big]^{1/p_0}\Big(-\hspace{-0.38cm}\int_{B_\rho}\frac{\Phi_{p_0}(|Du-DP|)}{\Phi_{p_0}(1+|DP|)}
 \,dx\Big)^{\frac{p_0-1}{p_0}}.
\end{align*}
Therefore, recalling the definitions of $C(x_0,\rho,P)$ and $E(x_0,\rho,P)$, 
and using Young's inequality, we obtain
\begin{equation}\label{I5}
\begin{aligned}
&\frac{|I_5|}{\Phi_{p_0-1}(1+|DP|)}\\
&\leq c \Big\{C(x_0,\rho,P)+\big[\mu(\sqrt{C(x_0,\rho,P)})\big]^{1/2}
[C(x_0,\rho,P)]^{1/2}\\
&\quad +\big[\mu(\sqrt{C(x_0,\rho,P)})\big]^{1/p_0}
 \left[C(x_0,\rho,P)\right]^{\frac{p_0-1}{p_0}}\Big\}\\
&\leq c\Big\{E(x_0,\rho,P)+\big[\mu(\sqrt{E(x_0,\rho,P)})\big]^{1/2}
\left[E(x_0,\rho,P)\right]^{1/2} \Big\}.
\end{aligned}
\end{equation}

Next we estimate $I_6$. Applying \eqref{xconti} and \eqref{zetaconti}, we have
\begin{align*}
I_6
&=-\hspace{-0.38cm}\int_{B_\rho} (\mathbf{a}(x_0,P(x_0),Du)-\mathbf{a}(x_0,P(x_0),DP)):D\varphi\, dx\\
&=-\hspace{-0.38cm}\int_{B_\rho} (\mathbf{a}(x_0,P(x_0),Du)-\mathbf{a}(x,P(x_0),Du)):D\varphi\, dx\\
&\quad +-\hspace{-0.38cm}\int_{B_\rho} (\mathbf{a}(x,P(x_0),Du)-\mathbf{a}(x,u,Du)):D\varphi\, dx\\
&\leq c -\hspace{-0.38cm}\int_{B_\rho} \omega(\rho) \log(e+1+|Du|) \Phi_{p_2-1}(1+|Du|)\, dx\\
&\quad +c-\hspace{-0.38cm}\int_{B_\rho} \mu(|u-P(x_0)|^2)\Phi_{p_0-1}(1+|Du|)\, dx\\
&=: c(I_{6a}+I_{6b}).
\end{align*}
For $I_{6a}$,  from Jensen's inequality for the convex function 
$\Phi_{p_0}^*$  we have 
$$
I_{6a}\leq (\Phi_{p_0}^*)^{-1} \Big(-\hspace{-0.38cm}\int_{B_\rho} 
\Phi_{p_0}^*\left(\omega(\rho)\log(e+1+|Du|) \Phi_{p_2-1}(1+|Du|)\right)\, dx\Big).
$$
Note that the integration on the right hand side above is exactly same 
to $\tilde I_2$ denoted in \eqref{tildeI2} with $2\rho$ replaced by $\rho$, 
hence it follows from \eqref{tildeI21} and (5) of Proposition \ref{propbasic11} 
that
\begin{align*}
I_{6a}
&\leq (\Phi_{p_0}^*)^{-1}\Big(\big[\omega(\rho)\log(\frac{1}{\rho})
 \big]^{\frac{p_0+1}{p_0}} \Phi_{p_0}(1+|DP|)\Big)\\
&\leq \big[\omega(\rho)\log(\frac{1}{\rho})\big]^{\frac{p_0^2-1}{p_0}} 
 \Phi_{p_0-1}(1+|DP|).
\end{align*}
In the above inequality, we   used that 
$(\Phi_p^*)^{-1}(\theta t)\leq \theta^{\frac{p-1}{p}} (\Phi_p^*)^{-1}(t)$,
 $\theta\in(0,1)$ and $t>0$, which can  be easily derived from (2) 
of Proposition \ref{propbasic11}.

As for $I_{6b}$, using \eqref{furtherass2}, H\"older's inequality, 
Jensen's inequality for the convex functions $\mu^{-1}$ and 
$\Psi(t):=\Phi_{p_0}(\Phi_{p_0-1}^{-1}(t^{\frac{2p_0-2}{2p_0-1}}))$, and 
\eqref{341}, we have
\begin{align*}
I_{6b}
&\leq  c\Big(-\hspace{-0.38cm}\int_{B_\rho} \mu(|u-P(x_0)|^2)\, dx\Big)^{\frac{1}{2p_0-1}}
 \Big(-\hspace{-0.38cm}\int_{B_\rho} \left[\Phi_{p_0-1}(1+|Du|)\right]^{\frac{2p_0-1}{2p_0-2}}
 \, dx\Big)^{\frac{2p_0-2}{2p_0-1}} \\
&\leq c \Big[\mu\Big(-\hspace{-0.38cm}\int_{B_\rho} |u-P(x_0)|^2\, dx\Big)\Big]^{\frac{1}{2p_0-1}}
\Big[\Psi^{-1}\Big(-\hspace{-0.38cm}\int_{B_\rho} \Phi_{p_0}(1+|Du|)\, dx\Big)
 \Big]^{\frac{2p_0-2}{2p_0-1}}\\
&\leq c \Big[\mu\Big(-\hspace{-0.38cm}\int_{B_\rho} |u-P(x_0)|^2\, dx\Big)\Big]^{\frac{1}{2p_0-1}}
 \left[\Psi^{-1}\left( \Phi_{p_0}(1+|DP|)\right)\right]^{\frac{2p_0-2}{2p_0-1}}\\
&\leq c \Big[\mu\Big(-\hspace{-0.38cm}\int_{B_\rho} |u-P(x_0)|^2\, dx\Big)\Big]^{\frac{1}{2p_0-1}}
 \Phi_{p_0-1}(1+|DP|).
\end{align*}
Therefore, from the previous two estimates and the definition of $E$, we have
\begin{equation}\label{I6}
  \frac{|I_6|}{\Phi_{p_0-1}(1+|DP|)}\leq c E(x_0,\rho,P).
\end{equation}
Inserting \eqref{I5} and \eqref{I6} into \eqref{pf201} and using the definition 
of \eqref{E}, we have \eqref{harmonicapprox} for every 
$\varphi\in C^\infty(B_{\rho}(x_0))$ with $\sup_{B_{\rho}}|D\varphi|= 1$. 
Therefore,  the standard normalization argument  yields \eqref{harmonicapprox} 
for every $\varphi\in C^\infty(B_{\rho}(x_0))$.
\end{proof}

We consider the affine function $P=(Du)_{x_0,\rho} (x-x_0)+(u)_{x_0,\rho}$, and set
\begin{gather}\label{CDu}
\begin{aligned}
 C(x_0,\rho)
&:= C(x_0,\rho,(Du)_{x_0,\rho} (x-x_0)+(u)_{x_0,\rho})\\
&= -\hspace{-0.38cm}\int_{B_\rho}\Big[\frac{|Du-(Du)_{x_0,\rho}|^2}{(1+|(Du)_{x_0,\rho}|)^2}
+\frac{\Phi_{p_0}(|Du-(Du)_{x_0,\rho}|)}{\Phi_{p_0}(1+|(Du)_{x_0,\rho}|)}\Big]\,dx,
\end{aligned} \\
\label{EDu}
\begin{aligned}
 \tilde E(x_0,\rho)
&:= E(x_0,\rho,(Du)_{x_0,\rho} (x-x_0)+(u)_{x_0, \rho})\\
&= C(x_0,\rho) +\Big[\mu\Big(-\hspace{-0.38cm}\int_{B_{\rho}} |u-(u)_{x_0,\rho}|^2\, dx\Big)
 +\omega(\rho)\log(\frac{1}{\rho})\Big]^{\frac{1}{2p_0-1}},
\end{aligned} \\
\label{EDubeta}
 E(x_0,\rho):=C(x_0,\rho) +\Big[\mu\Big(\rho -\hspace{-0.38cm}\int_{B_{\rho}} |Du|^2\, dx\Big) 
+\omega(\rho)\log(\frac{1}{\rho})\Big]^{\frac{1}{2p_0-1}}.
\end{gather}
Then, from Poincar\'e's inequality along with the fact $\rho<1$, we see that
\begin{equation}\label{EE}
\tilde E(x_0,\rho)\leq c E(x_0,\rho).
\end{equation}
some $c=c(n,N)\geq 1$.

\begin{lemma}\label{lemiteration}
For $\theta\in(0,1/8)$, there exists
$\epsilon_0=\epsilon_0(n,N,\gamma_1,\gamma_2,\nu,\Lambda,\theta)>0$
such that if $E(x_0,\rho)\leq\epsilon_0$ and 
$\rho\leq\epsilon_1:= \min\{\theta^n,\rho_0\}$, where $\rho_0$ is defined 
in \eqref{rho00}, then
\begin{equation}\label{CE1}
C(x_0,\theta\rho)\leq c_*\theta^2E(x_0,\rho)
\end{equation}
for some $c_*=c_*(n,N,\gamma_1,\gamma_2,\nu,\Lambda)\geq1$.
\end{lemma}

\begin{proof}
\textbf{Step 1.} We first estimate the  integrals
\begin{equation}\label{5678}
-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\frac{|u(x)-P_{2\theta\rho}|^2}{(2\theta\rho)^2}\, dx\quad
\text{and}\quad
-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\Phi_{p_0}\Big(\frac{|u-P_{2\theta\rho}|}{2\theta\rho}\Big)\, dx.
\end{equation}
where the affine function $P_{x_0,2\theta\rho}$ is defined in Section \ref{affine}.
 Recall  $\mathcal{A}$ and $w$ denoted in \eqref{cala} with 
$P=(Du)_{\rho} (x-x_0)+(u)_{\rho}$. Then
$$
w:= \frac{u-(Du)_{\rho}(x-x_0)-(u)_{\rho}}{\left(1+(Du)_{\rho}\right)
\sqrt{\tilde E(x_0,\rho)}}\quad \text{and}\quad 
-\hspace{-0.38cm}\int_{B_\rho} |Dw|^2\, dx \leq 1.
$$
Let $\epsilon\in(0,1)$ be a sufficiently small number determined later,
 for which we consider $\delta=\delta(n,N,\nu,\Lambda,\epsilon)>0$ 
determined in Lemma \ref{lemharmonicapprox}. Then by Lemma \ref{lemharmonic} 
 with \eqref{EDu}, \eqref{EDubeta} and \eqref{EE}, we have
$$
\big|-\hspace{-0.38cm}\int_{B_\rho(x_0)} \mathcal ADw :D\varphi \,dx\big|
\leq \delta \sup_{B_\rho(x_0)} |D\varphi|,
$$
by taking sufficiently small 
$\epsilon_0=\epsilon_0(n,N,\gamma_1,\gamma_2,\nu,\Lambda,\epsilon)\in(0,1)$. 
We point out that when using Lemma \ref{lemharmonic}, we need the assumption 
\eqref{lem1ass}, which is clear since 
$C(x_0,\rho)\leq E(x_0,\rho)\leq \epsilon_0\leq 1$. 
Therefore, in view of Lemma \ref{lemharmonicapprox}, there exists an
 $\mathcal A$-harmonic map $h$ such that
\begin{equation}\label{pf301}
-\hspace{-0.38cm}\int_{B_{\rho}}|Dh|^2\, dx\leq 1\quad \text{and}\quad 
-\hspace{-0.38cm}\int_{B_\rho}|w-h|^2\,dx\leq \epsilon \rho^2.
\end{equation}
We notice by a basic regularity theory for $\mathcal A$-harmonic maps, 
see for instance \cite[Chapter 10]{Gi1}, that
\begin{equation}\label{pf302}
\rho^{-2}\sup_{B_{\rho/2}}|Dh|^2+\sup_{B_{\rho/2}}|D^2h|\leq c \rho^{-2}-\hspace{-0.38cm}\int_{B_{\rho}}|Dh|^2\, dx\leq c\rho^{-2}.
\end{equation}
Moreover, the Taylor expansion of $h$ implies that for $\theta\in(0,1/4)$,
\begin{equation}\label{pf303}
\sup_{x\in B_{2\theta\rho}}|h(x)-h(x_0)-Dh(x_0)(x-x_0)|^2\leq c\theta^4\rho^2.
\end{equation}
At this point we choose
$\epsilon=\theta^{n+4}$.
Then we have from the second inequality in \eqref{pf301} and \eqref{pf303} that
$$
-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\frac{|w-h(x_0)-Dh(x_0)(x-x_0)|^2}{(2\theta\rho)^2}\, dx
\leq c\theta^2,
$$
hence, by the definitions of the affine function 
$P_{2\theta\rho}:=P_{x_0,2\theta\rho}$ denoted in Section \ref{affine} 
and $w$ and \eqref{EE}, we obtain
\begin{equation}\label{pf304}
-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\frac{|u-P_{2\theta\rho}|^2}{(2\theta\rho)^2}\, dx
\leq c\theta^2(1+|(Du)_\rho|)^2E(x_0,\rho).
\end{equation}

Next we estimate the second integral in \eqref{5678}. 
Let $t_0=t(x_0)\in(0,1)$ be a number satisfying
\begin{equation}\label{t0}
\frac{1}{p_0}=(1-t_0)+\frac{t_0}{p_0(1+\sigma_s)},
\end{equation}
where $\sigma_s$ is given in Lemma \ref{sobolevpoincare}. 
Note that since $2<\gamma_1\leq p_0\leq \gamma_2$, there exists 
$0<t_m\leq t_M<1$ depending only on $\gamma_1,\gamma_2,\sigma_s$ such that  
$t_m\leq t_0\leq t_M$. Then by H\"older's inequality, Jensen's inequality 
for the convex map $t\mapsto [(\Phi_{p_0})^{-1}(t^{p_0})]^2$, 
\eqref{pf304}, \eqref{sobopoin} and  (1) of Proposition \ref{propbasic11}, we have
\begin{align*}
&-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\Phi_{p_0}\Big(\frac{|u-P_{2\theta\rho}|}{2\theta\rho}\Big)\, 
dx \\
&\leq \Big(-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\Big[\Phi_{p_0}
 \Big(\frac{|u-P_{2\theta\rho}|}{2\theta\rho}\Big)\Big]^{1/p_0}
 \, dx\Big)^{(1-t_0)p_0} \\
&\quad\times \Big(-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\Big[\Phi_{p_0}
 \Big(\frac{|u-P_{2\theta\rho}|}{2\theta\rho}\Big)\Big]^{1+\sigma_s}\, dx
 \Big)^{\frac{t_0}{1+\sigma_s}}\\
&\leq \left[\Phi_{p_0}\left(\theta(1+|(Du)_\rho|)\sqrt{E(x_0,\rho)}\right)
\right]^{1-t_0}\Big(-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\Phi_{p_0}(|Du-DP_{2\theta\rho}|)
 \, dx\Big)^{t_0}\\
&\leq [\theta \sqrt{E(x_0,\rho)}]^{p_0(1-t_0)} 
 \left[\Phi_{p_0}(1+|(Du)_\rho|)\right]^{1-t_0}
 \Big(-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\Phi_{p_0}(|Du-DP_{2\theta\rho}|)\, dx\Big)^{t_0}.
\end{align*}
In addition,  from \eqref{DPDP}, (1) of Proposition \ref{propbasic11}, 
\eqref{sobopoin}, \eqref{DPDu} and the definition of $E$, we have
\begin{align*}
&-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\Phi_{p_0}(|Du-DP_{2\theta\rho}|)\, dx\\
&\leq c\theta^{-n}-\hspace{-0.38cm}\int_{B_{\rho}}\Phi_{p_0}(|Du-DP_{\rho}|)\, dx
 +c \Phi_{p_0}(|DP_{\rho}-DP_{2\theta\rho}|)\\
&\leq c\theta^{-n}-\hspace{-0.38cm}\int_{B_{\rho}}\Phi_{p_0}(|Du-DP_{\rho}|)\, dx
 +c\theta^{-(p+1)}-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_p
 \big(\frac{|u-P_{\rho}|}{\rho}\big)\, dx\\
&\leq c\theta^{-n-(p+1)}-\hspace{-0.38cm}\int_{B_{\rho}}\Phi_p(|Du-DP_{\rho}|)\, dx\\
&\leq c\theta^{-n-(p+1)}-\hspace{-0.38cm}\int_{B_{\rho}}\Phi_{p_0}(|Du-(Du)_{\rho}|)\, dx\\
&\leq c\theta^{-n-(p+1)}\Phi_{p_0}(1+|(Du)_{\rho}|)E(x_0,\rho).
\end{align*}
Combining the two above estimates, we obtain
\begin{align*}
&-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\Phi_{p_0}
\Big(\frac{|u-P_{2\theta\rho}|}{2\theta\rho}\Big)\, dx \\
&\leq c\theta^{p_0-(2p_0+n+1)t_0} \Phi_{p_0}
 (1+|(Du)_{\rho}|)E(x_0,\rho)^{(\frac{p_0}{2}-1)(1-t_0)+1}.
\end{align*}
Therefore, taking $\epsilon_0=\epsilon_0(n,N,\gamma_1,\gamma_2,\theta)>0$ 
sufficiently small so that
\begin{align*}
E(x_0,\rho)^{(p/2-1)(1-t_0)}
&\leq \epsilon_0^{(\frac{p_0}{2}-1)(1-t_0)}
 \leq \epsilon_0^{(\frac{\gamma_1}{2}-1)(1-t_M)}\\
&\leq  \theta^{-\gamma_1+(2\gamma_2+n+1)t_M+2}
 \leq \theta^{-p_0+(2p_0+n+1)t_0+2},
\end{align*}
we obtain
\begin{equation}\label{pf305}
-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\Phi_{p_0}
\Big(\frac{|u-P_{2\theta\rho}|}{2\theta\rho}\Big)\, dx
\leq  c\theta^2\Phi_{p_0}(1+|(Du)_{\rho}|)E(x_0,\rho).
\end{equation}

Moreover, by  assuming
$$
\sqrt{E(x_0,\rho)}\leq \sqrt{ \epsilon_0}\leq \theta^n/8,
$$
we have
\begin{equation}\label{pf306}
1+|(Du)_{\rho}|\leq 2(1+|(Du)_{\theta\rho}|)\quad \text{and}\quad
1+|(Du)_{2\theta\rho}|\leq 2(1+|(Du)_{\theta\rho}|),
\end{equation}
see \cite[page 483]{FM1}. Therefore, inserting the first inequality in 
\eqref{pf306} into \eqref{pf304} and \eqref{pf305}, we obtain
\begin{gather}\label{pf308}
-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\frac{|u-P_{2\theta\rho}|^2}{(2\theta\rho)^2}\, dx
\leq c\theta^2(1+|(Du)_{\theta\rho}|)^2E(x_0,\rho), \\
\label{pf309}
-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\Phi_{p_0}\big(\frac{|u-P_{2\theta\rho}|}{2\theta\rho}\big)
\, dx\leq c\theta^2 \Phi_{p_0}(1+|(Du)_{\theta\rho}|) E(x_0,\rho) .
\end{gather}
\smallskip

\noindent\textbf{Step 2.} Now we prove \eqref{CE1}. Suppose that
\begin{equation}\label{pf3erhoass1}
E(x_0,\rho) \leq \epsilon_0\leq \theta^{n}.
\end{equation}
Then we observe that
$$
C(x_0,8\theta\rho)\leq (8\theta)^{-n}E(x_0,\rho)\leq 1.
$$
Hence, in view of Lemma \ref{lemcaccio} with $\rho$ replaced by $\theta\rho$  
and $P=P_{2\theta\rho}$, we have
\begin{align*}
&-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_{p_0-2}(1+|DP_{2\theta\rho}|) |Du-DP_{2\theta\rho}|^2\, dx
 +-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_{p_0}(|Du-DP_{2\theta\rho}|)\,dx\\
&\leq c\Phi_{p_0-2}(1+|DP_{2\theta\rho}|) 
 -\hspace{-0.38cm}\int_{B_{2\theta\rho}}\frac{|u-P_{2\theta\rho}|^2}{(2\theta\rho)^2}\, dx
 +c-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\Phi_{p_0}\Big(\frac{|u-P_{2\theta\rho}|}{2\theta\rho}\Big)
 \,dx\\
&\quad  +c\Phi_{p_0}(1+|DP_{2\theta\rho}|)
\Big\{\mu\Big(-\hspace{-0.38cm}\int_{B_{2\theta\rho}} |u-(u)_{2\theta\rho}|^2\, dx \Big)
+\omega(2\theta\rho)\log\big(\frac{1}{2\theta\rho}\big)\Big\}.
\end{align*}
We note that Young's inequalities such as (4) and (5) of 
Proposition \ref{propbasic11} also hold for the Orlicz function 
$G(t):=t^{\frac{p}{2}}\log(e+\sqrt t)$. Using this fact and Lemma 
\ref{lemDP} with $(\rho,\theta)$ replaced by $(\theta\rho,1/2)$, we have
\begin{align*}
&-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_{p_0-2}(1+|(Du)_{\theta\rho}|) 
 |Du-(Du)_{\theta\rho}|^2\, dx \\
&\leq c-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_{p_0-2}(1+|DP_{2\theta\rho}|) 
 |Du-(Du)_{\theta\rho}|^2\, dx\\
&\quad +c-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_{p_0-2}(|(Du)_{\theta\rho}
 -DP_{\theta\rho}|)|Du-(Du)_{\theta\rho}|^2\, dx\\
&\quad +c-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_{p_0-2}(|DP_{\theta\rho}
 -DP_{2\theta\rho}|)|Du-(Du)_{\theta\rho}|^2\, dx\\
&\leq c-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_{p_0-2}(1+|DP_{2\theta\rho}|) 
 |Du-DP_{2\theta\rho}|^2\, dx+ c-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_{p_0}
 (|Du-(Du)_{\theta\rho}|)\, dx\\
&\quad +c\left(\Phi_{p_0}(|(Du)_{\theta\rho}-DP_{\theta\rho}|)
 +\Phi_{p_0}(|DP_{\theta\rho}-DP_{2\theta\rho}|)\right)\\
&\leq c-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_{p_0-2}(1+|DP_{2\theta\rho}|) 
 |Du-DP_{2\theta\rho}|^2\, dx+ c-\hspace{-0.38cm}\int_{B_{\theta\rho}}
 \Phi_{p_0}(|Du-(Du)_{\theta\rho}|)\, dx\\
&\quad +c-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\Phi_{p_0}
 \big( \frac{|u-P_{2\theta\rho}|}{2\theta\rho}\big)\,dx.
\end{align*}
Combining the two estimates above  with 
$$
-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_{p_0}(|Du-(Du)_{\theta\rho}|)\, dx
\leq c-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_{p_0}(|Du-DP_{2\theta\rho}|)\, dx,
$$
we obtain
\begin{equation}\label{caccioP}
\begin{aligned}
&-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_{p_0-2}(1+|(Du)_{\theta\rho}|) 
 |Du-(Du)_{\theta\rho}|^2\, dx
 +-\hspace{-0.38cm}\int_{B_{\theta\rho}}\Phi_{p_0}(|Du-(Du)_{\theta\rho}|)\,dx \\
&\leq c\Phi_{p_0-2}(1+|DP_{2\theta\rho}|) -\hspace{-0.38cm}\int_{B_{2\theta\rho}}
 \frac{|u-P_{2\theta\rho}|^2}{(2\theta\rho)^2}\, dx
+c-\hspace{-0.38cm}\int_{B_{2\theta\rho}}\Phi_{p_0}
 \Big(\frac{|u-P_{2\theta\rho}|}{2\theta\rho}\Big)\,dx\\
&\quad +c\Phi_{p_0}(1+|DP_{2\theta\rho}|)
\Big\{\mu\Big(-\hspace{-0.38cm}\int_{B_{2\theta\rho}} |u-(u)_{2\theta\rho}|^2\, dx \Big)
+\omega(2\theta\rho)\log\big(\frac{1}{2\theta\rho}\big)\Big\}.
\end{aligned}
\end{equation}
We further estimate the right hand side on the above inequality. 
Applying \eqref{DPDu}, \eqref{pf3erhoass1} and \eqref{pf306}, we have
\begin{align*}
\Phi_{p_0}(|DP_{2\theta\rho}|)
&\leq c\Phi_{p_0}(|DP_{2\theta\rho}-(Du)_{2\theta\rho}|)+ c\Phi_{p_0}(|(Du)_{2\theta\rho}|)\\
&\leq c\theta^{-n}-\hspace{-0.38cm}\int_{B_\rho}\Phi_{p_0}(|Du-(Du)_{\rho}|)\, dx+ c \Phi_{p_0}(|(Du)_{2\theta\rho}|)\\
&\leq c(\theta^{-n}E(x_0,\rho)+1) \Phi_{p_0}(1+|(Du)_{\rho}|)\\
&\leq c \Phi_{p_0}(1+|(Du)_{\theta\rho}|)\\
&\leq \Phi_{p_0}\left(c(1+|(Du)_{\theta\rho}|)\right),
\end{align*}
which implies 
$$
|DP_{2\theta\rho}|\leq c(1+|(Du)_{\theta\rho}|).
$$
Moreover, using Poincar\'e's inequality and the fact that $\rho\leq \theta^n$, 
we have
$$
-\hspace{-0.38cm}\int_{B_{2\theta\rho}} |u-(u)_{2\theta\rho}|^2\, dx
\leq c\theta^{-n} -\hspace{-0.38cm}\int_{B_{\rho}} |u-(u)_{\rho}|^2\, dx 
\leq c\rho -\hspace{-0.38cm}\int_{B_{\rho}} |Du|^2\, dx.
$$
Therefore, inserting the two above inequalities, \eqref{pf308} 
and \eqref{pf309} into \eqref{caccioP}, we obtain
$$
C(x_0,\theta\rho)\leq c\theta^2E(x_0,\rho)+c[E(x_0,\rho)]^{2p_0-1}
$$
Finally, assuming
$$
[E(x_0,\rho)]^{2(p_0-1)}\leq [E(x_0,\rho)]^{2(\gamma_1-1)} 
\leq \epsilon_0^{2(\gamma_1-1)} \leq \theta^2,
$$
we obtain \eqref{CE1}.
\end{proof}


Now, we are ready to prove Theorem \ref{mainthm}

\begin{proof}[Proof of Theorem \ref{mainthm}]
Fix $\alpha\in(0,1)$. Let us determine several constants such that
\begin{gather} \label{gamma}
\lambda:=n-2(1-\alpha)\in(n-2,n), \\
\label{theta}
\theta:=\min\big\{\frac{1}{8},\frac{1}{\sqrt{c_*}},\frac{1}{3^{1/(n-\lambda)}}\big\}\\
\label{epsilon2}
\epsilon_2:=\min\big\{\frac{\theta^n}{16},\frac{\epsilon_0}{2}\big\},
\end{gather}
where $c_*$ and $\epsilon_0$ are determined in Lemma \ref{lemiteration}.
We note from the continuity of $\mu(\cdot)$ and \eqref{vlogholder} that  
one can find $\delta_1=\delta_1(\mu,\omega,\epsilon_2)>0$ such that
\begin{equation}\label{delta1}
\mu(r)+\omega(r)\log\big(\frac{1}{r}\big)\leq \epsilon_2  \ \ \text{for every }r\in(0,\delta_1].
\end{equation}
Then we  denote
\begin{equation}\label{rhom}
\rho_m:=\min\{ \delta_1,\epsilon_1,\rho_0 \}.
\end{equation}
\smallskip

\noindent\textbf{Step 1.} 
In this step, we fix $B_{\rho}=B_{\rho}(x_0)\subset \Omega$ with $x_0\in\Omega$ 
and $\rho\in(0,\rho_m]$, and suppose that
\begin{equation}\label{induction1}
C(x_0,\rho)\leq \epsilon_2\quad \text{and}\quad 
M(x_0,\rho):=\rho-\hspace{-0.38cm}\int_{B_\rho}|Du|^2\,dx\leq \delta_1.
\end{equation}
Then we claim that for any $k=0,1,2,\dots$,
\begin{equation}\label{induction2}
C(x_0,\theta^k\rho)\leq \epsilon_2\quad  \text{and} \quad
 M(x_0,\theta^k\rho):=\theta^{k}\rho-\hspace{-0.38cm}\int_{B_{\theta^k\rho}}|Du|^2\,dx\leq \delta_1.
\end{equation}
We prove the claim by induction. For the sake of convenience, for $k=0,1,2,\dots$, 
we write  $\eqref{induction2}_{k,1}$ (resp. $\eqref{induction2}_{k,2}$) 
with the first (resp. second) inequality in \eqref{induction2}.

Now we suppose that the inequalities in \eqref{induction2} hold for $k$, 
and then prove that \eqref{induction2} holds for $k$ replaced by $k+1$.
 We first observe from \eqref{induction2}$_{k,1}$ and H\"older's inequality  that
\begin{equation}\label{pf401}
\begin{aligned}
 -\hspace{-0.38cm}\int_{B_{\theta^k\rho}}|Du-(Du)_{\theta^k\rho}|^2\, dx
&\leq  (1+|(Du)_{\theta^k\rho}|)^2 C(x_0,\theta^k\rho)  \\
&\leq  2 \epsilon_2\Big(1+-\hspace{-0.38cm}\int_{B_{\theta^k\rho}} |Du|^2\,dx\Big).
\end{aligned}
\end{equation}
and so, by \eqref{induction2}$_{k,2}$,
$$
\theta^k\rho -\hspace{-0.38cm}\int_{B_{\theta^k\rho}}|Du-(Du)_{\theta^k\rho}|^2\, dx  
\leq  2 \epsilon_2 \theta^k\rho+2 \epsilon_2\delta_1.
$$
This together with \eqref{theta}, \eqref{epsilon2} and \eqref{rhom} imply
\begin{align*}
M(x_0,\theta^{k+1}\rho)
&\leq 2\theta^{k+1}\rho-\hspace{-0.38cm}\int_{B_{\theta^{k+1}\rho}}|Du-(Du)_{\theta^k\rho}|^2\,dx
 +2\theta^{k+1}\rho|(Du)_{\theta^k\rho}|^2\\
&\leq 2\theta^{1-n}\theta^k\rho-\hspace{-0.38cm}\int_{B_{\theta^k\rho}}|Du-(Du)_{\theta^k\rho}|^2\,dx
 +2\theta M(x_0,\theta^k\rho)\\
&\leq 4\theta^{k+1-n} \epsilon_2\rho+4\theta^{1-n} \epsilon_2\delta_1
 +2\theta\delta_1\\
&\leq 4\theta^{-n} \epsilon_2\rho+4\theta^{-n} \epsilon_2\delta_1
 +2\theta\delta_1
\leq \delta_1,
\end{align*}
which shows \eqref{induction2}$_{k+1,2}$. It remains to prove 
\eqref{induction2}$_{k+1,1}$.  We notice from \eqref{EDubeta}, \eqref{epsilon2}, 
\eqref{delta1}, \eqref{induction2}$_{k,1}$, \eqref{induction2}$_{k,2}$ and the 
fact that $\theta^k\rho\leq \rho_m\leq \delta_1$ that
$$
E(x_0,\theta^k\rho):=\epsilon_2 
+\big[\mu(\delta_1) +\omega(\theta^k\rho)\log\big(\frac{1}{\theta^k\rho}
\big)\big]^{\frac{1}{2p_0-1}}\leq 2\epsilon_2<\epsilon_0.
$$
Therefore, applying Lemma \ref{lemiteration}, we have
$$
C(x_0,\theta^{k+1}\rho) \leq c_*\theta^2\epsilon_2\leq \epsilon_2.
$$
This shows \eqref{induction2}$_{k+1,1}$. Hence, we prove that 
\eqref{induction2} holds for every $k=0,1,2,\dots$. 
Here, we note that \eqref{pf401} also holds for every $k=0,1,2,\dots$.

Next  we observe from \eqref{epsilon2} and \eqref{pf401} that
\begin{align*}
-\hspace{-0.38cm}\int_{B_{\theta^{k+1}\rho}(x_0)}|Du|^2\,dx
&\leq  2-\hspace{-0.38cm}\int_{B_{\theta^{k+1}\rho}}|Du-(Du)_{\theta^k\rho}|^2\,dx 
 +2|(Du)_{\theta^k\rho}|^2\\
&\leq  2\theta^{-n}-\hspace{-0.38cm}\int_{B_{\theta^{k}\rho}}|Du-(Du)_{\theta^k\rho}|^2\,dx
  +2-\hspace{-0.38cm}\int_{B_{\theta^k\rho}}|Du|^2\,dx\\
&\leq  4\theta^{-n}\epsilon_2 + (4\theta^{-n}\epsilon_2+2)
 -\hspace{-0.38cm}\int_{B_{\theta^k\rho}} |Du|^2\,dx\\
&\leq  4\theta^{-n}\epsilon_2 + 3-\hspace{-0.38cm}\int_{B_{\theta^k\rho}} |Du|^2\,dx,
\end{align*}
and so, by \eqref{theta},
$$
\int_{B_{\theta^{k+1}\rho}}|Du|^2\,dx
\leq \theta^\lambda\int_{B_{\theta^k\rho}} |Du|^2\,dx+4|B_1|(\theta^k\rho)^n.
$$
Then applying Lemma \ref{lemtech} with $\phi(r)=\int_{B_r}|Du|^2\,dx$, 
 for every $r\in(0,\rho)$, we have
\begin{equation}\label{pf403}
\begin{aligned}
\int_{B_r}|Du|^2\,dx
&\leq c\big\{\big(\frac{r}{\rho}\big)^\lambda\int_{B_{\rho}} |Du|^2\,dx
 +  r^\lambda\big\}\\
&\leq \frac{c}{\rho^\lambda}\Big(\int_{\Omega} |Du|^2\,dx+ 1\Big) r^\lambda.
\end{aligned}
\end{equation}
\smallskip

\noindent\textbf{Step 2.} To complete the proof, we define
$$
\Omega_u:=\left\{x_0\in\Omega: u\in C^{\tilde \alpha}(U_{x_0},\mathbb{R}^N) 
\text{ for every $\tilde \alpha\in(0,1)$  and for some }U_{x_0}\subset\Omega\right\}.
$$
where $U_{x_0}$ is an open neighborhood of $x_0$. Suppose that 
$x_0\in\Omega$ satisfies
\begin{equation}\label{pf406}
\begin{gathered}
\liminf_{\rho\downarrow0}-\hspace{-0.38cm}\int_{B_\rho(x_0)}|Du-(Du)_{\rho}|\,dx=0,\\
M_{x_0}:=\limsup_{\rho\downarrow0} -\hspace{-0.38cm}\int_{B_\rho(x_0)}|Du|^2\,dx<\infty.
\end{gathered}
\end{equation}
For $p\in[\gamma_1,\gamma_2]$, set $t_p\in(0,1)$ such that
\begin{equation}\label{7890}
\frac{1}{p}=2t_p+\frac{(1-t_p)}{p(1+\frac{\sigma_h}{2})},
\end{equation}
where $\sigma_h$ is determined in Theorem \ref{thmhigher}.
Note that Since $p\in[\gamma_1,\gamma_2]\mapsto t_p$ is continuous there 
exists $t_m=t_m(\gamma_1,\gamma_2,\sigma_h)\in(0,1)$ such that $t_m\leq t_p$ 
for every $p\in[\gamma_1,\gamma_2]$. We define
\begin{equation}\label{sdef}
s:=\min\Big\{\Big[\frac{(\epsilon_2/2)^{\gamma_2/2}}{c_1(M_{x_0}+2)}
\Big]^{\frac{1}{\gamma_1t_m}} ,\delta_1\Big\}<1,
\end{equation}
where $c_1>0$ will be determined later. Then in view of \eqref{pf406}, 
one can find
\begin{equation}\label{rho11}
\rho<\min\{\rho_m,(2^n(M_{x_0}+1)+1)^{-1}\delta_1 \}.
\end{equation}
 with $B_{3\rho}(x_0)\subset\Omega$ such that
\begin{equation}\label{pf404}
-\hspace{-0.38cm}\int_{B_{\rho}(x_0)}|Du-(Du)_{\rho}|dx<s\quad \text{and}\quad 
-\hspace{-0.38cm}\int_{B_{2\rho}(x_0)}\Phi(x,|Du|)\,dx<M_{x_0}+1.
\end{equation}
Then by H\"older's inequality with \eqref{7890},
\begin{align*}
-\hspace{-0.38cm}\int_{B_\rho(x_0)}\Phi_{p_0}(|Du-(Du)_{\rho}|)\,dx 
&\leq \Big(-\hspace{-0.38cm}\int_{B_\rho(x_0)}[\Phi_{p_0}(|Du-(Du)_{\rho}|)]^\frac{1}{2p_0}\,dx
\Big)^{2p_0t_0}  \\
&\quad\times \Big(-\hspace{-0.38cm}\int_{B_\rho(x_0)}[\Phi_{p_0}(|Du-(Du)_{\rho}|)]^{1+\frac{\sigma_h}{2}}\,dx
 \Big)^{\frac{1-t_0}{1+\frac{\sigma_h}{2}}},
\end{align*}
where $t_0:=t_{p_0}$. From Jensen's inequality for the convex function 
$G(t):=\Phi_{p_0}^{-1}(t^{2p_0})$, we have
\begin{align*}
-\hspace{-0.38cm}\int_{B_\rho(x_0)}[\Phi_{p_0}(|Du-(Du)_{\rho}|)]^\frac{1}{2p_0}\,dx
&\leq G^{-1}\Big(-\hspace{-0.38cm}\int_{B_\rho(x_0)}|Du-(Du)_{\rho}|\,dx\Big) \\
&< G^{-1}(s) = [\Phi_{p_0}(s)]^{\frac{1}{2p_0}} \\
&\leq (2s^{p_0})^{\frac{1}{2p_0}}.
\end{align*}
On the other hand, using Jensen's inequality for the convex map 
$t\mapsto [\Phi(t)]^{1+\frac{\sigma_h}{2}}$, \eqref{higher3} and
 \eqref{pf404}, we have
\begin{align*}
-\hspace{-0.38cm}\int_{B_{\rho}(x_0)}[\Phi_{p_0}(|Du-(Du)_{\rho}|)]^{1+\frac{\sigma_h}{2}}\,dx 
&\leq c -\hspace{-0.38cm}\int_{B_{\rho}(x_0)}[\Phi_{p_0}(|Du|)]^{1+\frac{\sigma_h}{2}}\,dx \\
&\leq c \Big\{\Big(-\hspace{-0.38cm}\int_{B_{2\rho}(x_0)}\Phi(x,|Du|)\,dx\Big)^{1+\frac{\sigma_h}{2}}
 +1\Big\} \\
&< c(M_{x_0}+2)^{1+\frac{\sigma_h}{2}}.
\end{align*}
Therefore, we have
$$
-\hspace{-0.38cm}\int_{B_\rho(x_0)}\Phi_{p_0}(|Du-(Du)_{\rho}|)\,dx 
< c s^{p_0t_0}(M_{x_0}+2)^{1-t_0}
\leq   c_1(M_{x_0}+2)s^{\gamma_1 t_m},
$$
where $c_1>0$ depends only on $n,N,\gamma_1,\gamma_2,\nu,\Lambda$, and so by 
\eqref{sdef},
$$
-\hspace{-0.38cm}\int_{B_\rho(x_0)}\Phi_{p_0}(|Du-(Du)_{\rho}|)\,dx 
< \big(\frac{\epsilon_2}{2}\big)^{\gamma_2/2}.
$$
On the other hand, by \eqref{rho11} and \eqref{pf404}, we have
$$
\rho-\hspace{-0.38cm}\int_{B_\rho(x_0)}|Du|^2\,dx
\leq \rho \Big(2^n-\hspace{-0.38cm}\int_{B_{2\rho}}\Phi(x,|Du|)\,dx+1\Big)< \delta_1.
$$
In addition, by the continuity of integrals, there exists $\varrho>0$ 
such that \eqref{pf404} with $x_0$ replaced by $y$ holds for every 
$y\in B_{\varrho}(x_0)$. Without loss of generality, we can assume that 
$\varrho\leq \rho$. We note that using Jensen's inequality for the convex 
function $\Psi_{p(y)}(t):=\Phi_{p(y)}(\sqrt t)$,  for 
$y\in B_{\varrho}(x_0)$ we have
\begin{align*}
C(y,\rho)
&\leq \Psi_{p(y)}^{-1}\Big(-\hspace{-0.38cm}\int_{B_\rho(y)}\Phi_{p_0}(|Du-(Du)_{\rho}|)\,dx\Big)
+-\hspace{-0.38cm}\int_{B_\rho(y)}\Phi_{p(y)}(|Du-(Du)_{\rho}|)\,dx \\
&\leq  \Psi_{p(y)}^{-1}\big(\big(\frac{\epsilon_2}{2}\big)^{\gamma_2/2}\big)
+\big(\frac{\epsilon_2}{2}\big)^{\gamma_2/2} \\
&\leq \big(\frac{\epsilon_2}{2}\big)^{\frac{\gamma_2}{2}\frac{2}{p(y)}}
+\big(\frac{\epsilon_2}{2}\big)^{\gamma_2/2} \\
&\leq \frac{\epsilon_2}{2}+\frac{\epsilon_2}{2} \leq \epsilon_2,
\end{align*}
where $\Psi_{p}(t)=t^\frac{p}{2}\log(e+t)$ hence 
$\Psi_{p(y)}^{-1}(t)\leq t^{\frac{2}{\gamma_1}}$.
Therefore, in view of Step 1, we see that  \eqref{pf403} 
with $x_0$ replaced by $y$ hold for every $y\in B_\varrho(x_0)$ and $t\leq \rho$. 
Therefore, by Morrey-Campanato's embedding theorem, we have $u$ is
 $C^\alpha(B_{\varrho}(x_0),\mathbb{R}^N)$, that is, $x_0\in\Omega_u$.
\end{proof}

\begin{thebibliography}{00}

\bibitem{AM1} E. Acerbi and G. Mingione;
\emph{Gradient estimates for the $p(x)$-Laplacean system},
J. Reine Angew. Math. \textbf{584} (2005), 117-148.

\bibitem{BCM1} P. Baroni, M. Colombo, G. Mingione; 
\emph{Harnack inequalities for double phase functionals}, Nonlinear Anal. TMA,
 \textbf{121} (2015), 206-222.

\bibitem{BCM2} P. Baroni, M. Colombo, G. Mingione; 
\emph{Non-autonomous functionals, borderline cases and related function classes}, 
St. Petersburg Math. J. \textbf{27} (3) (2016), 347-379.

\bibitem{BFM1} V. B\"ogelein, M. Foss, G. Mingione; 
\emph{Regularity in parabolic systems with continuous coefficients},
 Math. Z. \textbf{270} (3-4) (2012), 903-938.

\bibitem{BDHS1} V.  B\"ogelein, F. Duzaar, J. Habermann, C. Scheven; 
\emph{Partial Hölder continuity for discontinuous elliptic problems with 
VMO-coefficients}, Proc. Lond. Math. Soc. (3) \textbf{103} (3) (2011), 371-404.

\bibitem{Cam1} S. Campanato; 
\emph{Partial H\"older continuity of solutions of quasilinear parabolic 
systems of second order with linear growth}, Rend. Sem. Math. Univ. Padova 
\textbf{64} (1981), 59-75.

\bibitem{Cam2} S. Campanato; 
\emph{On the nonlinear parabolic systems in divergence form H\"older 
continuity and partial H\"older continuity of the solutions}, 
Ann. Mat. Pura Appl. \textbf{137} (4) (1984), 83-122.

\bibitem {Cia1} A. Cianchi; 
\emph{Boundedness of solutions to variational
problems under general growth conditions}, Comm. Part. Diff. Equ.,
 \textbf{22} (1997), 1629-1646.

\bibitem {Cia2} A. Cianchi; 
\emph{Local boundedness of minimizers of anisotropic functionals}.
 Ann. Inst. H. Poincar\'e Anal. Non Lin\`eaire,  \textbf{17} (2000),
147-168.

\bibitem{ColM1} M. Colombo, G. Mingione; 
\emph{Regularity for double phase variational problems}, Arch. Ration. 
Mech. Anal. \textbf{215} (2) (2015), 443-496.

\bibitem{ColM2} M. Colombo, G. Mingione; 
\emph{Bounded minimizers of double phase variational integrals}, 
Arch. Ration. Mech. Anal.\textbf{218} (1) (2015), 219-273.

\bibitem{ColM3} M. Colombo, G. Mingione; 
\emph{Calder\'on-Zygmund estimates and non-uniformly elliptic operators}, 
J. Funct. Anal. \textbf{270} (4) (2016), 1416-1478.

\bibitem{DE1} L. Diening, F. Ettwein; 
\emph{Fractional estimates for non-differentiable elliptic systems with 
general growth}, Forum Math., \textbf{20} (3) (2008), 523-556.

\bibitem{DHHR1} L. Diening, P. Harjulehto, P. H\"{a}st\"{o}, M. Ru\v zi\v cka; 
\emph{Lebesgue and Sobolev Spaces with Variable Exponents}, Lecture Notes in Math., 
vol. 2017, Springer-Verlag, Berlin, 2011.

\bibitem{DS1} F. Duzaar and K. Steffen; 
\emph{Optimal interior and boundary regularity for almost minimizers to 
elliptic variational integrals}, J. Reine Angew. Math. \textbf{546} (2002), 73-138.

\bibitem{EMM1} M. Eleuteri, P. Marcellini, E. Mascolo; 
\emph{Lipschitz continuity for energy integrals with variable exponents}, 
Atti Accad. Naz. Lincei Rend. Lincei Mat. Appl. \textbf{27} (2016),  61--87.

\bibitem{EMM2} M. Eleuteri, P. Marcellini, E. Mascolo;
\emph{Lipschitz estimates for systems with ellipticity conditions at infinity}, 
Ann. Mat. Pura Appl., to appear.

\bibitem {FM1} M. Foss, G. Mingione; 
\emph{Partial continuity for elliptic problems}, Ann. Inst. H. Poincar\'e Anal. 
Non Lin\'eaire \textbf{25} (3) (2008), 471-503.

\bibitem{GP1} F. Giannetti, A. Passarelli di Napoli; 
\emph{Regularity results for a new class of functionals with non-standard 
growth conditions}, J. Differential Equations, \textbf{254} (3) (2013), 1280-1305.

\bibitem{Gi1} E. Giusti; 
\emph{Direct methods in the calculus of variations}, World Scientific, 
River Edge, NJ, 2003.

\bibitem{Ha1} J. Habermann; 
\emph{Partial regularity for nonlinear elliptic systems with continuous 
growth exponent}, Ann. Mat. Pura Appl. (4) \textbf{192} (3) (2013), 475-527.

\bibitem{HHK1} P. Harjulehto, P. H\"ast\"o, R. Kl\'en; 
\emph{Generalized Orlicz spaces and related PDE},  Nonlinear Anal. TMA,
 \textbf{143} (2016), 155-173.

\bibitem{HHT1} P. Harjulehto, P. H\"ast\"o, O. Toivanen; 
\emph{Harnack's inequality for quasiminimizers with generalized Orlicz growth 
conditions}, preprint.

\bibitem{IV1} T. Iwaniec, A Verde; 
\emph{On the operator ${L}(f)=f\log|f|$}, J. Funct. Anal., 
\textbf{169} (2) (1999), 391-420.

\bibitem{Kr1} M. Kronz; 
\emph{Partial regularity results for minimizers of quasiconvex functionals 
of higher order}, Ann. Inst. H. Poincar\'e Anal. Non Lin\'eaire \textbf{19} (2) (2002), 81-112.

\bibitem{KM1} T. Kuusi, G. Mingione; 
\emph{Nonlinear potential theory of elliptic systems}, Nonlinear Anal. TMA,
 \textbf{138} (2016), 277-299.

\bibitem{KM2} T. Kuusi, G. Mingione; 
\emph{Vectorial nonlinear potential theory}, J. Eur. Math. Soc., to appear.

\bibitem{Le1} F. Leonetti; 
\emph{Higher integrability for minimizers of integral functionals with non-standard
growth,} J. Differential Equations \textbf{112} (1994), 308--324.

\bibitem{Le2} F. Leonetti, F. Siepe; 
\emph{Integrability for solutions to some anisotropic elliptic equations}.
 Nonlinear Anal. TMA, \textbf{75} (2012), 2867-2873.

\bibitem{Li1} G. M. Lieberman; 
\emph{The natural generalization of the natural conditions of Ladyzenskaja 
and Ural\^atzeva for elliptic equations},
Comm. Partial Differential Equations, \textbf{16} (2-3) (1991), 311-361.

\bibitem{M2} P. Marcellini; 
\emph{Regularity of minimizers of integrals of the calculus of variations 
with non standard growth conditions}, Arch.~Ration.~Mech.~Anal., \textbf{105}
(1989), 267-284.

\bibitem{M3} P. Marcellini;
\emph{Regularity and existence of solutions of elliptic equations with $p,q$-growth
conditions}, J.~Differential Equations \textbf{90} (1991), 1-30.

\bibitem{M4} P. Marcellini; 
\emph{Regularity for elliptic equations with general growth conditions}, 
J.~Differential Equations \textbf{105} (1993), 296-333.

\bibitem{M5} P. Marcellini; 
\emph{Everywhere regularity for a class of elliptic systems without growth 
conditions}, Ann.~Scuola Norm.~Sup.~Pisa Cl.~Sci.~(IV), \textbf{23} (1996),
1-25.

\bibitem{Min1} G. Mingione; 
\emph{Regularity of minima: an invitation to the dark side of the calculus of 
variations}, Appl. Math., \textbf{51} (4) (2006), 355-426.

\bibitem{Ok1} J. Ok; 
\emph{Gradient estimates for elliptic equations with $L^{p(\cdot)} \log L$ growth}, 
Calc. Var. Partial Differential Equations, \textbf{55} (2) (2016), 55:26.


\bibitem{Ok2} J. Ok; 
\emph{Harnack inequality for a class of functionals with non-standard growth 
via De Giorgi's method},  Adv. Nonlinear Anal. DOI 10.1515/anona-2016-0083.

\bibitem{Ok3} J. Ok; 
\emph{Regularity results for a class of obstacle problems with non-standard growth},
 J. Math. Anal. Appl. \textbf{444} (2) (2016), 957-979.


\bibitem{Ok4} J. Ok; 
\emph{Calder\'on-Zygmund estimates for a class of obstacle problems with 
non-standard growth}, NoDEA Nonlinear Differential Equations 
Appl. \textbf{23} (4) (2016), 1-21

\end{thebibliography}

\end{document}