\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2018 (2018), No. 58, pp. 1--15.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2018 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2018/58\hfil Gradient estimate for elliptic obstacle problems] {Gradient estimate in Orlicz spaces for elliptic obstacle problems with partially BMO nonlinearities} \author[S. Liang, S. Z. Zheng \hfil EJDE-2018/58\hfilneg] {Shuang Liang, Shenzhou Zheng} \address{Shuang Liang \newline Department of Mathematics, Beijing Jiaotong University, Beijing 100044, China} \email{shuangliang@bjtu.edu.cn} \address{Shenzhou Zheng (corresponding author)\newline Department of Mathematics, Beijing Jiaotong University, Beijing 100044, China} \email{shzhzheng@bjtu.edu.cn} \dedicatory{Communicated by Zhaosheng Feng} \thanks{Submitted October 2, 2017. Published March 1, 2018.} \subjclass[2010]{35D30, 35K10} \keywords{Nonlinear elliptic obstacle problems; partially BMO nonlinearities; \hfill\break\indent Reifenberg flatness; Orlicz space; the Hardy-Littlewood maximal operator} \begin{abstract} We prove a global Orlicz estimate for the gradient of weak solutions to a class of nonlinear obstacle problems with partially regular nonlinearities in nonsmooth domains. It is assumed that the nonlinearities are merely measurable in one spatial variable and have sufficiently small BMO semi-norm in the other variables, and the boundary of underlying domain is Reifenberg flat. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{definition}[theorem]{Definition} \newtheorem{assumption}[theorem]{Assumption} \allowdisplaybreaks \section{Introduction} Throughout this paper, let $\Omega$ be a bounded domain $\mathbb{R}^{n} $ for $n\ge 2$ with the non-smooth boundary specified later. Suppose that $\mathbf{F}=(f^1,f^2,\dots,f^n)$ is a given measurable vectorial-valued function, and $\mathbf{a}=\mathbf{a}(\xi,x): \mathbb{R}^n\times \Omega\to \mathbb{R}^n$ is assumed to be a Carath\'{e}odory vectorial-valued function which is measurable in $x\in \Omega$ for each $\xi\in \mathbb{R}^n$ and Lipschitz continuous in $\xi\in \mathbb{R}^n$ for each $x\in \Omega$. Let $\psi$ be an obstacle function with \begin{equation*} \psi \in W^{1,2}(\Omega), \quad \text{and} \quad \psi \leq 0 \quad \text{on } \partial\Omega; \end{equation*} and recall an admissible set $\mathcal{K}_\psi(\Omega)$ by \begin{equation*} \mathcal{K}_\psi(\Omega)=\big\{v \in W_0^{1,2}(\Omega): v\geq \psi \text{ a. e. in } \Omega \big\}. \end{equation*} We devote this present article to study a global estimate in Orlicz spaces for the gradient of weak solutions to the following variational inequalities with some weak regular assumptions on the datum in the sense of distribution. \begin{definition}\label{def-weak-solution} \rm We say that $u$ is a weak solution of variational inequalities \eqref{weak-ellip-equa}, if $u \in \mathcal{K}_\psi(\Omega)$ satisfies \begin{equation}\label{weak-ellip-equa} \int_\Omega \langle \mathbf{a}(Du,x), D(u-v)\rangle\,dx \leq \int_\Omega \langle \mathbf{F}, D(u-v)\rangle\,dx \end{equation} for all $v \in\mathcal{K}_\psi(\Omega)$. \end{definition} To ensure solvability in $L^2(\Omega)$ of \eqref{weak-ellip-equa}, it is necessary to impose some basic structural assumptions on the nolinearities with ellipticity and growth: there exist two constants $0<\nu \le \Lambda <\infty$ such that \begin{equation}\label{growth-ellipticity-condition} \begin{gathered} \langle D_{\xi}\mathbf{a}(\xi,x) \eta \cdot \eta \rangle \ge \nu |\eta|^2, \\ |\mathbf{a}(\xi,x)| + |\xi||D_{\xi}\mathbf{a}(\xi,x)| \le \Lambda |\xi| \end{gathered} \end{equation} for a.e. $x \in \Omega$ and $\xi,\eta \in \mathbb{R}^{n}$, where $D_\xi$ denotes the differentiation in $\xi \in \mathbb{R}^n$, and $\langle \cdot,\cdot\rangle$ is the standard inner product in $\mathbb{R}^n$. Consequently, the condition \eqref{growth-ellipticity-condition} readily yields that \begin{equation}\label{monotonicity-condition} \begin{gathered} \mathbf{a}(0,x) = 0,\\ \nu |\xi-\eta|^2 \leq \langle \mathbf{a}(\xi,x)-\mathbf{a}(\eta,x),\xi -\eta \rangle. \end{gathered} \end{equation} With \eqref{growth-ellipticity-condition} in hand, by the Minty-Browder argument then there exists a unique weak solution $u\in \mathcal{A}$ to the variational inequality \eqref{weak-ellip-equa} with the usual $L^2$-estimate \begin{equation}\label{energy-estimate} \int_\Omega |Du|^2 \leq c \int_\Omega (|F|^2+|D\psi|^2)\,dx \end{equation} with a positive constant $c=c(n,\nu,\Lambda)$. Let us review some recent progresses related to our research. The regularity on nonlinear elliptic problems under the weak regular datum is a classical and important topic in the aspect of partial differential equations. Indeed, the Calder\'{o}n-Zygmund theory is an extremely popular research to various elliptic and parabolic equations with some minimal regular datum, see \cite{ByK16,ByOP16,DoK10,DoK11-1,KiK07,TiZ17,ZhZ17}. In particular, the Calder\'{o}n-Zygmund theory regarding elliptic problems with partially BMO coefficients has been recently getting largely attention. As we know, Kim and Krylov in \cite{KiK07} was first to obtain the Calder\'{o}n-Zygmund theory to nondivergence linear elliptic and parabolic equations with partially VMO coefficients. Then, this study with partially regular coefficients was extended to divergence or nondivergence linear elliptic and parabolic equations/system and linear equations of higher order by Dong-Kim-Krylov, see \cite{DoK10,DoK11-1,DoK11-2} etc. Later, Byun and Palagachev \cite{ByP14} also deduced a global weighted $L^p$-theory to linear elliptic problems with small partially BMO coefficients over non-smooth domains via a rather different geometrical approach. In particular, we would like to mentioned that Byun and Kim \cite{ByK16} very recently attained the nonlinear Calder\'{o}n-Zygmund theory to elliptic equations with measurable nonlinearities in nonsmooth domains based on their usual geometrical approach. In fact, this article is also motivated by Byun and Kim's recent work. We would like to remark that this partial BMO assumption is actually a sort of minimal regular requirements on the coefficients for elliptic operators even for linear elliptic settings to ensure a satisfactory Calder\'{o}n-Zygmund theory for all $p > 1$. Indeed, this was verified by a famous counterexample due to Ural'tseva \cite{Ur67}, where he constructed an equation in $\mathbb{R}^d \, (d\ge 3)$ with the coefficients depending only on the first two coordinates so that there is no unique solvability in Sobolev spaces $W^{1,p}$ for any $p>1$. Nonlinear elliptic equations with discontinuous nonlinearities in the spatial variable are related to nonlinear problems in medium composition materials. Especially, these problems with partial regular nonlinearities are particularly attributed to the so-called laminate materials \cite{ChKV86}. Meanwhile, the relevant obstacle problems usually appear in various fields such as physics, biology, economics, computer science and so on, which is to describe practical phenomena in a situation with a kind of constraint by the so-called obstacle function. Here, we would particularly like to point out that the study of this article was inspired by a recent progress from Byun and Kim's work in \cite{ByK16}, which is first concerned with the nonlinear Calder\'{o}n-Zygmund theory involved in the measurable nonlinearities. They actually considered the following nonlinear elliptic equation with partially BMO nonlinearities in Reifenberg domains (see Definitions below): \begin{gather*} \operatorname{div}\mathbf{a}(Du,x)= \operatorname{div}\mathbf{F},\quad \text{in } \Omega, \\ u=0, \quad \text{on } \partial \Omega, \end{gather*} and obtained that \begin{equation*} \mathbf{F}\in L^{p}(\Omega,\mathbb{R}^{n}) \Rightarrow Du\in L^{p}(\Omega,\mathbb{R}^{n}), \quad 2\leq p<\infty \end{equation*} for the weak solutions $u\in W^{1,2}_0(\Omega)$. Orlicz spaces are the natural generalizations of Lebesgue spaces, and the estimates in Orlicz spaces for partial differential equations have become an extremely popular research nowadays. Areas of its applications include the study of geometric, probability, stochastic, Fourier analysis and so on, also see \cite{PKJF13,RaR02}. The regularity in Orlicz spaces is actually an extension of the classical Calder\'on-Zygmund estimates for the theory of PDEs. Just for divergence elliptic case, under some regular assumptions on the datum it implies that $f\in \mathcal{B} \Rightarrow Du\in\mathcal{B}$ for a given Orlicz space $\mathcal{B}$. For instance, Azroul et al. \cite{AzBT00} first proved that for each radial solution $u$ for Poisson equation $-\Delta u=f$, it satisfies $\frac{\partial^2u}{\partial x_i\partial x_j} \in L_{loc}^\Phi(B_r)$ if $f \in L^\Phi(B_r)$ for any Young function $\Phi$ with $\Phi(|f(x)|)\log(|x|)$ being integrable. Later, Jia et al. \cite{JiLW07,JiLW11} established Orlicz regularities for above-mentioned Poisson equation and divergence linear elliptic equations with small BMO coefficients in Reifenberg or quasicovex domains via the Hardy-Littlewood maximal functions, respectively. In particular, Byun-Ok-Palagachev\cite{ByOP16} proved the weighted Orlicz estimates for divergence linear parabolic systems while the leading coefficients are assumed to be only measurable in one spatial variable and have small BMO seminorms in the remaining variables. In addition, there were various gradient estimates in Orlicz spaces for $p$-Laplacian, quasilinear $p$-Laplacian and evolution $p$-Laplacian in Reifenberg flat domain, respectively, see \cite{YaSZ08,YaZ12}. Finally, we would also like to mention that Li-Zhang-Zheng \cite{LiZZ17} obtained a local Orlicz estimate of the Hessian strong solutions to a class of nondivergence linear elliptic equations $a_{ij}D_{ij}u=f(x)$ with partially BMO nonlinearities. Let us start with related basic notations and definitions which will be useful in this paper. \begin{definition} \label{Young-function} \rm Let $\Phi$ be a nonnegative, increasing and convex real-valued function on $[0,+\infty)$. If it satisfies \begin{equation}\label{Young-function-condition} \lim_{\rho\to 0^+} \frac{\Phi(\rho)}{\rho} =\lim_{\rho\to +\infty} \frac{\rho}{\Phi(\rho)}=0, \end{equation} where $\Phi(0)=0$, $\Phi(\infty)=\lim_{\rho\to +\infty} \Phi(\rho)$, then we say $\Phi$ is a Young function. \end{definition} \begin{definition}\label{2-condition} We say that the Young function $\Phi$ satisfies the $\Delta_2\cap\nabla_2$ condition, denoted by $\Phi\in\Delta_2\cap\nabla_2$, if \begin{itemize} \item[(i)] ($\Delta_2$ condition) there exists a positive constant $\mu_1$ such that \begin{equation}\label{2-condition-1} \Phi(2\rho)\leq \mu_1\Phi(\rho), \quad \forall\rho>0; \end{equation} \item[(ii)] ($\nabla_2$ condition) there exists a constant $\mu_2>1$ such that \begin{equation}\label{2-condition-2} \Phi(\rho)\leq \frac{1}{2\mu_2}\Phi(\mu_2\rho), \quad \forall\rho>0. \end{equation} \end{itemize} \end{definition} Indeed, the limits \eqref{Young-function-condition} along with $\Delta_2\cap\nabla_2$ mean that \begin{equation*} 0=\Phi(0)=\underset{\rho\to 0^+}{\lim}\Phi(\rho), \quad \lim_{\rho\to +\infty}\Phi(\rho)=+\infty, \end{equation*} which show that the limits are neither too slow nor too fast while $\rho\to 0^+$ and $\rho\to +\infty$, see \cite{RaR02}. We also notice that the $\Delta_2$ condition implies that there exists a constant $\mu(\lambda)>1$ such that \begin{equation*} \Phi(\lambda\rho)\leq \mu(\lambda)\Phi(\rho), \quad \forall\rho>0, \lambda>1. \end{equation*} Since $\Phi\in\Delta_2$, there exist two constants $t_1$ and $t_2$ with $10:\int_{\Omega}\Phi \Big(\frac{|f|}{\lambda}\Big)dx\leq 1\Big\}. \end{equation} As usual, the Orlicz Sobolev spaces $W^1L^{\Phi}(\Omega)$ is defined by the function spaces of all measurable functions $v\in L^{\Phi}(\Omega)$ such that its gradient vector $Dv\in L^{\Phi}(\Omega)$ with the norm \begin{equation*}\label{Orlicz-Sobolev-norm} \|v\|_{W^1L^{\Phi}(\Omega)}=\|v\|_{L^{\Phi}(\Omega)}+\|Dv\|_{L^{\Phi}(\Omega)}. \end{equation*} We can refer it to \cite{KoK91} for more details about Orlicz spaces. It is easy to observe that Orlicz spaces $L^{\Phi}$ generalize $L^p$ spaces in the sense that if we take $\Phi(x)=x^p$ with $p>1$ so that we have \begin{equation*} L^{\Phi}(\Omega)=L^{p}(\Omega), \quad W^1L^{\Phi}(\Omega)=W^{1,p}(\Omega). \end{equation*} Let us denote a type point by $x=(x^1,x')=(x^1,x^2,\dots,x^n)\in \mathbb{R}^n$, and let $B_r=\{x\in \mathbb{R}^n:|x|0\} $ with $Q_r(y)=Q_r+y, Q_r^+(y)=Q_r^++y$; and some typical boundaries by $\Omega_r(y)=Q_r(y)\cap \Omega$, $\partial_\omega\Omega_r(y)=Q_r(y)\cap\partial\Omega$, $T_r=Q_r\cap\{x^1=0\}$. We write an average of $f(x)$ in $Q_r$ for $r>0$ to be \begin{equation*} -\hspace{-0.38cm}\int_{Q_r} f(x)\,\,dx = \frac{1}{|Q_r|} \int_{Q_r} f(x)\,\,dx, \end{equation*} where $|Q_r|$ is $n$-dimensional Lebesgue measure of $Q_r$. The $(n-1)$-dimensional average of $f(x)$ in $Q'_r$ with respect to $x'$ by \begin{equation*} \bar{f}_{Q'_r}(x^{1}) = -\hspace{-0.38cm}\int_{Q'_r} f(x^{1},x')\,\,dx' = \frac{1}{|Q'_r|} \int_{Q'_r} f(x^{1},x')\,\,dx' \end{equation*} with $|Q'_r|$ as the $(n-1)$-dimensional Lebesgue measure of $Q'_r$. To impose a partially regular assumption on $\mathbf{a}(\xi,x)=\mathbf{a}(\xi,x^1,x')$ (cf. \cite[Definition 2.2]{ByK16}), we consider a function \begin{equation}\label{small BMO semi-norm} \theta(\mathbf{a},Q_r(y))=\sup_{\xi\in \mathbb{R}^n\setminus 0 } \frac{|\mathbf{a}(\xi,x^1,x')-\bar{\mathbf{a}}_{Q_r'(y')}(\xi,x^1)| }{|\xi|} \end{equation} with \begin{equation}\label{bar-a} \bar{\mathbf{a}}_{Q_r'(y')}(\xi,x^1)=-\hspace{-0.38cm}\int_{Q_r'(y')}\mathbf{a}(\xi,x^1,z')dz', \end{equation} where $\mathbf{a}(\xi,x)$ is zero extended from $\Omega\cap Q_r'$ to $Q_r'\setminus \Omega\cap Q_r'$, \begin{assumption}\label{small BMO} \rm We say that ($\mathbf{a}(\xi,x), \Omega$) is $(\delta,R)$-vanishing of codimension 1, if for every point $x_0\in \Omega$ and for any $0< r \le R$ with \begin{eqnarray*} \operatorname{dist}(x_0, \partial \Omega) = \min_{z\in \partial \Omega} \operatorname{dist}(x_0, z) > \sqrt{2} r, \end{eqnarray*} there exists a coordinate system depending only on $x_0$ and $r$, whose variables are still denoted by $x$, such that in the new coordinate system with $x_0$ as the origin and \begin{eqnarray*} -\hspace{-0.38cm}\int_{Q_r} \big| \theta(\mathbf{a}, Q_r)(x) \big|^2 \,dx \le \delta^2; \end{eqnarray*} while, for $x_0\in \Omega$ with \[ \operatorname{dist}(x_0, \partial \Omega) = \min_{z\in \partial \Omega} \operatorname{dist}(x_0, z) = \operatorname{dist}(x_0, z_0) \le \sqrt{2} r, \] where $z_0\in \partial \Omega$, one has that there exists a coordinate system depending on $x_0$ and $00$ will be specified later. \end{assumption} Now we state our main result of this paper as follows. \begin{theorem}\label{main-result} Let $u\in \mathcal{K}_\psi(\Omega)$ be a weak solution to the variational inequalities \eqref{weak-ellip-equa} with nonlinearity $\textbf{a}(\xi,x)$ satisfying the structural conditions \eqref{growth-ellipticity-condition}, and let $(\mathbf{a}(\xi,x),\Omega)$ satisfy $(\delta,R)$-vanishing of codimension 1 as Assumption \ref{small BMO}. For the Young function $\Phi\in\Delta_2\cap\nabla_2$, if $\Psi^2 \in L^\Phi(\Omega)$ with $\Psi=|F|+|D\psi|$, then there exists a small constant $\delta_0=\delta_0(n,\nu,\Lambda,|\Omega|,\Phi)>0$ such that for every $\delta \in (0,\delta_0]$, we have $|Du|^2 \in L^\Phi(\Omega)$ with the estimate \begin{equation}\label{Lorentz-estimate-a} \| |Du|^2\|_{L^\Phi}(\Omega) \leq c\left(\|\Psi^2\|_{L^\Phi(\Omega)}+1\right), \end{equation} where the positive constant $c$ is independent of $u$ and $\Psi$. \end{theorem} Here, we reach it mainly based on the Byun-Wang's geometric argument \cite{ByW12}. In particular, our key argument was actually inspired by Byun et al's recent papers \cite{ByK16,ByOP16}, and we make use of the boundedness of the Hardy-Littlewood maximal functions in Orlicz spaces, modified Vitali covering and an equivalent representation of Orlicz norm to prove the main theorem for nonlinear elliptic obstacle problem under the minimal weak regular assumptions on the nonlinearities and the boundary of domain. We would also like to point out that our consideration is twofolds: one is to extend Byun-Wang's work in \cite{ByW12} by assuming that the nonlinearities are \textit{partially BMO} instead of small BMO oscillations. The other is that our problems are involved in the obstacle constraints in more general Orlicz spaces instead of Lebesgue spaces in \cite{ByK16}. For this, however some comparison estimates for the reference problems can be cited from \cite{ByK16}. The rest of the paper is organized as follows. Section 2 is devoted to introduce some useful lemmas. In section 3, we focus on proving our main theorem. \section{Technical tools} In the section we present some useful lemmas, which will play essential roles in proving our main conclusions. We are mainly devoted to make some comparison estimates to the reference problems, and we particularly make use of Byun-Kim's important work on the interior and boundary Lipschitz regularity for limiting equations whose nonlinearities depend on the gradients of weak solutions and only one variable. Let us denote by $c(n,\nu,\Lambda,\dots)$ a universal constant depending only on prescribed quantities and possibly varying from line to line in the context. First of all, let us introduce the Hardy-Littlewood maximal function and related basic facts, see \cite{ByW12,ByK16}. \begin{definition}\label{Hardy-Littlewood-maximal-funct}\rm Let $f$ be a locally integrable function of $\mathbb{R}^n$, the Hardy-Littlewood maximal function $\mathcal{M}f$ is defined by \begin{equation*} \mathcal{M}f(x)=\underset{r>0}{\sup}-\hspace{-0.38cm}\int_{B_r(x)}|f(y)|dy. \end{equation*} If $f$ is confined in a bounded subset $U$ of $\mathbb{R}^n$, then we can define a restricted maximal function $\mathcal{M}_Uf$ in the following form \begin{equation*} \mathcal{M}_Uf=\mathcal{M}(f\chi_U), \end{equation*} where $\chi_U$ is the standard characteristic function on $U$. \end{definition} For the Hardy-Littlewood maximal function, we immediately conclude the following two useful classical properties, for details see \cite{St93}. \begin{itemize} \item[(i)] (strong $(p,p)$-type) If $f\in L^p(\mathbb{R}^n)$ for $1\alpha\}| \leq\frac{c(n)}{\alpha} \|f\|_{L^1(\mathbb{R}^n)}, \quad \forall \alpha>0. \end{equation} \end{itemize} Further, we have the following boundedness of the Hardy-Littlewood maximal function in Orlicz spaces. \begin{lemma}\label{boundedness} If $\Phi$ is a Young function satisfying the $\Delta_2\cap\nabla_2$-condition, then there exists a positive constant $c=c(n,\Phi)$ such that \begin{equation}\label{boundedness-equa-1} \int_{\mathbb{R}^n}\Phi(|f|)dx \leq \int_{\mathbb{R}^n}\Phi(\mathcal{M}f)dx \leq c\int_{\mathbb{R}^n}\Phi(|f|)dx \end{equation} for all $f\in L^\Phi(\mathbb{R}^n)$. In addition, we would like to point out that from inequality \eqref{mononicity-condition} and the Luxemburg norm \eqref{Luxemburg-norm} we have \begin{equation}\label{boundedness-equa-2} c^{-1}(\|f\|^\alpha_{L^\Phi(U)}-1) \leq \int_{U}\Phi(|f|)dx \leq c(\|f\|^\beta_{L^\Phi(U)}+1), \end{equation} where $\alpha=t_1,\beta=t_2$ satisfy \eqref{mononicity-condition} and the constant $c>1$ is independent of $f$. \end{lemma} Next, we use that the nonlinear elliptic obstacle problems under consideration is an invariant under scaling and normalization, see \cite[Lemma 3.1]{ByK16}. \begin{lemma}\label{scaling-normalization} For each $K,\rho>0$, we define \begin{equation*} \tilde{\mathbf{a}}(\xi,x)=\frac{\mathbf{a}(K\xi,\rho x)}{K},\quad \tilde u(x)=\frac{u(\rho x)}{K\rho}, \quad \tilde{\mathbf{F}}(x)=\frac{\mathbf{F}(\rho x)}{K}, \quad \tilde \psi(x)=\frac{\psi(\rho x)}{K\rho} \end{equation*} and the set $\tilde\Omega=\{\frac{x}{\rho}: x\in \Omega\}$, then we have \begin{itemize} \item[(i)] If $u \in \mathcal{K}_\psi(\Omega)$ is a weak solution of \eqref{weak-ellip-equa}, then $\tilde u\in \mathcal{K}_{\tilde\psi}(\tilde\Omega)$ is a weak solution of \[ \int_{\tilde\Omega} \langle \tilde{\mathbf{a}}(D\tilde u,x), D(\tilde u-\tilde v)\rangle\,dx \leq \int_{\tilde\Omega} \langle \tilde{\mathbf{F}}, D(\tilde u-\tilde v)\rangle\,dx, \] for all $\tilde v \in \mathcal{K}_{\tilde\psi}(\tilde\Omega)$. \item[(ii)] If the nonlinearity $\mathbf{a}(\xi,x)$ satisfies assumption \eqref{growth-ellipticity-condition}, then so dose $\tilde{\mathbf{a}}(\xi,x)$ with the same constants $\nu, \Lambda$. \item[(iii)] If the nonlinearity $(\mathbf{a}(\xi,x), \Omega)$ is $(\delta,R)$-vanishing of codimension 1 in $\Omega$, then $(\tilde{\mathbf{a}}(\xi,x), \Omega)$ is $(\delta,\frac{R}{\rho})$-vanishing of codimension 1 in $\tilde\Omega$. \end{itemize} \end{lemma} Let us now focus on some comparison estimates to a few of the related reference problems. Recalling that the domain $\Omega$ is assumed to be the $(\delta, R)$-Reifenberg flatness as a necessary minimal geometric condition in the new coordinate system. This leads to the following measure density conditions: \[ \sup_{0-16\delta \}, \\ \label{small-BMO-a} -\hspace{-0.38cm}\int_{Q_6}|\theta(a,Q_6)|^2dx\leq \delta^2. \end{gather} Now we are mainly to focus on the boundary estimates for the reference problems since the interior estimates are very similar to the boundary setting with a simpler procedure. We consider a local weak solution $u \in W^{1,2}(\Omega_6)$ of the variational inequalities \eqref{weak-ellip-equa} in $\Omega_6$ with $u=0$ on $\partial_\omega\Omega_6$. Note that it holds the measure density for Reifenberg flat domain $\Omega$, then we let $k,v\in W^{1,2}(\Omega_6), w\in W^{1,2}(\Omega_5)$ and $h\in W^{1,2}(Q_5^+)$, respectively, be the weak solutions of the following boundary value problems \begin{gather}\label{comparison-equ-obstacle} \begin{gathered} \operatorname{div}\mathbf{a}(Dk,x)= \operatorname{div}\mathbf{a}(D\psi,x), \quad \text{in } \Omega_6, \\ k=u, \quad \text{on } \partial\Omega_6; \end{gathered} \\ \label{comparison-equ-1} \begin{gathered} \operatorname{div}\mathbf{a}(Dv,x)= 0, \quad \text{in }\Omega_6, \\ v=k, \quad \text{on }\partial\Omega_6; \end{gathered} \\ \label{comparison-equ-2} \begin{gathered} \operatorname{div}\bar{\mathbf{a}}_{Q_5'}(Dw,x^1)= 0, \quad \text {in } \Omega_5, \\ w=v, \quad \text {on } \partial\Omega_5; \end{gathered}\\ \label{comparison-equ-3} \begin{gathered} \operatorname{div} \bar{\mathbf{a}}_{Q_5'}(Dh,x^1)=0, \quad \text {in } Q_5^+, \\ h=0, \quad \text{on } T_5; \end{gathered} \end{gather} where \[ \bar{\mathbf{a}}_{Q_5'}(\xi,x_1) = \begin{cases} -\hspace{-0.38cm}\int_{\Omega\cap Q_5'} \mathbf{a}(\xi,x_1,x') \,\,\,dx' & x\in \Omega\cap Q_5', \\ 0 & x\in Q_5'\backslash \Omega\cap Q_5'. \end{cases} \] Here, we would particularly like to point that it is really our starting point on the interior and boundary Lipschitz regularity for limiting problem \eqref{comparison-equ-2}, whose nonlinearities depend on the gradients of weak solutions and only one variable $x_1$, for details see Section 4 in Byun and Kim's work \cite{ByK16}. In what follows, we give some boundary comparison estimates among the above various reference problems, whose argument is vary similar to Byun et al's recent series papers \cite{ByW12,ByP14,ByK16,ByOP16}. The following comparison principle is rather necessary to ensure that each solution satisfies the admissible test functions as for the variational inequalities with an obstacle constraint. \begin{lemma}\label{comparison} Let $U\in \mathbb{R}^{n}$ be a bounded open domain. Suppose that $\psi,k \in W^{1,2}(U)$ satisfy \begin{gather*} -\operatorname{div} \mathbf{a}(D\psi,x)\leq -\operatorname{div} \mathbf{a}(Dk,x) \quad \text{in } U, \\ \psi \leq k, \quad \text {on } \partial U, \end{gather*} in the weak sense that \begin{equation}\label{comparison-estimate} \int_{U}\langle \mathbf{a}(D\psi,x)-\mathbf{a}(Dk,x),D\varphi \rangle \leq 0 \quad \text{for all } \varphi \in W_0^{1,2}(U) \text{ with } \varphi \geq 0 \end{equation} and $(\psi- k)^+ \in W_0^{1,2}(U)$. Then it holds $\psi \leq k$, a. e. in $U$. \end{lemma} \begin{proof} Taking $\varphi=(\psi- k)^+$ as a test function in \eqref{comparison-estimate} yields \begin{equation*} -\hspace{-0.38cm}\int_{\{x\in U:\psi>k\}}\langle \mathbf{a}(D\psi,x)-\mathbf{a}(Dk,x),D(\psi-k) \rangle \leq 0. \end{equation*} From the monotonic increasing \eqref{monotonicity-condition} of $\mathbf{a}(\xi,x)$, we find that \begin{equation*} \nu-\hspace{-0.38cm}\int_{\{x\in U:\psi>k\}}|D\psi-Dk|^2dx= 0, \end{equation*} which implies that $D\psi=Dk$ a.e. in $\{x\in U:\psi>k\}$ and therefore $D((\psi- k)^+)=0$ a.e. in $U$. Since $(\psi- k)^+\in W_0^{1,2}(U)$, we conclude that $\psi\leq k$ a.e. in $U$. \end{proof} Now we are in a position to show the comparison estimates among the various reference problems. \begin{lemma}\label{comparison-1} Let $u \in W^{1,2}(\Omega_6)$ be a local weak solution of \eqref{weak-ellip-equa} in $\Omega_6$ with $u=0$ on $\partial_\omega\Omega_6$, and let $v\in W^{1,2}(\Omega_6)$ be the weak solution of \eqref{comparison-equ-1}. Under the normalization conditions of \begin{equation}\label{CZ-operator-estimate-a} -\hspace{-0.38cm}\int_{\Omega_6}|Du|^2dx \leq 1, \quad -\hspace{-0.38cm}\int_{\Omega_6}\Psi^2dx \leq \delta^2, \end{equation} if for any $o<\epsilon<1$ there exists a constant $\delta=\delta(n,\epsilon,\nu,\Lambda)$ such that $(\mathbf{a}(\xi,x),\Omega)$ satisfying $(\delta,R)$-vanishing of codimension $1$ shown as \eqref{boundary-condtion-a} and \eqref{small-BMO-a}, then we derive that \begin{gather}\label{energy-estimate-v} -\hspace{-0.38cm}\int_{\Omega_6}|Dv|^2dx \leq c , \\ \label{comparison-estimate-a} -\hspace{-0.38cm}\int_{\Omega_6}|Du-Dv|^2dx \leq \epsilon. \end{gather} \end{lemma} \begin{proof} Thanks to the standard $L^2$-estimates of \eqref{comparison-equ-obstacle} and \eqref{comparison-equ-1}, it follows from \eqref{CZ-operator-estimate-a} that \begin{equation}\label{energy-estimate-v-a} -\hspace{-0.38cm}\int_{\Omega_6} |Dv|^2\,dx \leq c-\hspace{-0.38cm}\int_{\Omega_6} |Dk|^2\,dx \leq c-\hspace{-0.38cm}\int_{\Omega_6} (|Du|^2+|D\psi|^2)\,dx \leq c(1+\delta^2). \end{equation} In view of Lemma \ref{comparison} and \eqref{comparison-equ-obstacle}, we know that $k \geq \psi $, a.e. in $\Omega_6$. By extending $k$ by $u$ in $\Omega\setminus\Omega_6$, then it leads to $k\in \mathcal{K}_\psi(\Omega)$. Taking a test function $k \in W^{1,2}(\Omega_6)$ into \eqref{weak-ellip-equa}, we obtain \begin{equation}\label{weak-formula-a-1} \int_{\Omega_6} \langle \mathbf{a}(Du,x),Dk-Du\rangle\,dx \geq \int_{\Omega_6} \langle \mathbf{F},Dk-Du\rangle\,dx. \end{equation} Similarly, by taking a test function $k-u \in W_0^{1,2}(\Omega_6)$ into \eqref{comparison-equ-1} it yields \begin{equation}\label{weak-formula-a-2} \int_{\Omega_6} \langle \mathbf{a}(Dk,x),Dk-Du\rangle\,dx = \int_{\Omega_6} \langle \mathbf{a}(D\psi,x),Dk-Du\rangle\,dx. \end{equation} Then, by subtracting \eqref{weak-formula-a-1} from \eqref{weak-formula-a-2}, and by \eqref{monotonicity-condition} we obtain \begin{equation}\label{weak-formula-a} \begin{aligned} \nu \int_{\Omega_6}|Dk-Du|^2dx &\le \int_{\Omega_6} \langle \mathbf{a}(Dk,x)-\mathbf{a}(Du,x),Dk-Du\rangle\,dx \\ &\leq \int_{\Omega_6} \langle \mathbf{a}(D\psi,x)-\mathbf{F},Dk-Du\rangle\,dx. \end{aligned} \end{equation} On the other hand, by \eqref{growth-ellipticity-condition}, \eqref{CZ-operator-estimate-a} and Young inequality with $\varepsilon \in (0,1)$ it follows that \begin{equation}\label{weak-formula-esti-left-b} \begin{aligned} &\int_{\Omega_6} \langle \mathbf{a}(D\psi,x)-\mathbf{F},Dk-Du\rangle\,dx\\ &\leq \int_{\Omega_6} \left( |\mathbf{a}(D\psi,x)|+|\mathbf{F}|\right)|Dk-Du|\,dx \\ &= \int_{\Omega_6} \left( \Lambda|D\psi|+|\mathbf{F}|\right)|Dk-Du|\,dx \\ & \leq \varepsilon \int_{\Omega_6} |Dk-Du|^2\,dx +c(\varepsilon) \int_{\Omega_6} \left( \Lambda|D\psi|+|\mathbf{F}|\right)^2dx \\ &\leq \varepsilon \int_{\Omega_6} |Dk-Du|^2\,dx +c(\varepsilon,\Lambda)\delta^2. \end{aligned} \end{equation} Putting \eqref{weak-formula-a} and \eqref{weak-formula-esti-left-b} with $\varepsilon=\nu/2$ together yields \begin{equation}\label{comparison-estimate-a-1} -\hspace{-0.38cm}\int_{\Omega_6} |Dk-Du|^2\,dx \leq c\delta^2 \end{equation} with $c=c(\nu,\Lambda)$. Next, we subtract \eqref{comparison-equ-1} to \eqref{comparison-equ-obstacle}, and take a testing function $k-v \in W_0^{1,2}(\Omega_6)$ to obtain \[ \int_{\Omega_6} \langle \mathbf{a}(Dk,x)-\mathbf{a}(Dv,x),Dk-Dv\rangle\,dx = \int_{\Omega_6} \langle \mathbf{a}(D\psi,x),Dk-Dv\rangle\,dx. \] In a similar way as the above estimate \eqref{comparison-estimate-a-1}, there exists $\delta>0$ such that \begin{equation}\label{comparison-estimate-a-2} -\hspace{-0.38cm}\int_{\Omega_6}|Dk-Dv|^2dx \leq c\delta^2. \end{equation} Consequently, using \eqref{comparison-estimate-a-1} and \eqref{comparison-estimate-a-2}, it yields \begin{equation*} -\hspace{-0.38cm}\int_{\Omega_6}|Du-Dv|^2dx \leq 2-\hspace{-0.38cm}\int_{\Omega_6}|Du-Dk|^2dx +2-\hspace{-0.38cm}\int_{\Omega_6}|Dk-Dv|^2dx \leq c\delta^2. \end{equation*} Since $\epsilon>0$ is arbitrary, we choose a small $\delta>0$ such that $c\delta^2=\epsilon$, which reduces the desired estimate \eqref{comparison-estimate-a}. With the corresponding $\delta=\delta(n,\epsilon,\nu,\Lambda)$ and \eqref{energy-estimate-v-a} we obtain \eqref{energy-estimate-v}. \end{proof} Regarding the remainders of comparison estimates for the reference problems we only recall Byun and Kim's conclusions, see \cite[Lemmas 5.6 and 5.8]{ByK16}. \begin{lemma}\label{comparison-2} Assume that $u \in W^{1,2}(\Omega_6)$ is a weak solution of \eqref{weak-ellip-equa} in $\Omega_6$ with $u=0$ on $\partial_\omega\Omega_6$ under the assumptions of \eqref{boundary-condtion-a}, \eqref{small-BMO-a} and \eqref{CZ-operator-estimate-a}. If $v,w,h$ are the weak solutions of \eqref{comparison-equ-1}, \eqref{comparison-equ-2} and \eqref{comparison-equ-3}, respectively, then we have \begin{gather*} -\hspace{-0.38cm}\int_{\Omega_5}|Dw|^2dx \leq c-\hspace{-0.38cm}\int_{\Omega_6}|Dv|^2dx\leq c, \quad \|D\bar h\|_{L^\infty(\Omega_3)}\leq c_1, \\ -\hspace{-0.38cm}\int_{\Omega_5}|Dv-Dw|^2dx \leq \epsilon, \quad -\hspace{-0.38cm}\int_{\Omega3}|Dw-D\bar h|^2dx \leq \epsilon, \end{gather*} where $c_1=c_1(n,\nu,\Lambda)$, and $\bar h$ is a zero extension of $h$ from $Q_5^+$ to $Q_5$. \end{lemma} With Lemmas \ref{comparison-1} and \ref{comparison-2}, we conclude the following estimates near the boundary. \begin{lemma}\label{comparison-3} Let $u \in W^{1,2}(\Omega_6)$ of \eqref{weak-ellip-equa} in $\Omega_6$ with $u=0$ on $\partial_\omega\Omega_6$. If, for $o<\epsilon<1$ there exists a constant $\delta=\delta(n,\epsilon,\nu,\Lambda)$ with \eqref{boundary-condtion-a}, \eqref{small-BMO-a} and \eqref{CZ-operator-estimate-a}, then \begin{equation*} -\hspace{-0.38cm}\int_{\Omega_{3}}|Du-D\bar h|^2dx \leq \epsilon \quad \text{and } \quad \|D\bar h\|_{L^\infty(\Omega_3)}\leq c_1, \end{equation*} where $\bar h$ is the zero extension of $h$ from $Q_5^+$ to $Q_5$. \end{lemma} The following version of the modified Vitali covering plays an important role to the Calder\'on-Zygmund type theory over a Reifenberg flat domain, see \cite{ByW12,ByP14,ByK16}. \begin{lemma}\label{vitali} Let $C$ and $D$ be measurable sets with $C\subset D\subset \Omega$, and $\Omega$ be $(\delta, 1)$-Reifenberg flat for some small $\delta>0$. Assume that there exists a small $\varepsilon>0$ with $$ |C|<\varepsilon |B_1|, $$ further, for all $x \in \Omega$ and $r\in (0, 1]$ with $|C\cap B_r(x)| \ge |B_r(x)|$ it holds $B_r(x)\cap \Omega \subset D$. Then we have $$ |C|\le \Big(\frac {10}{1-\delta}\Big)^n |D|. $$ \end{lemma} Thanks to the above boundary comparison estimates, we conclude the following measure comparison estimate of distributions on the maximal function based on the modified Vitali covering lemma concerning the invariance property under the scaling argument and normalization. \begin{lemma}\label{boundary-maximal-function} Let $u\in W_0^{1,2}(\Omega)$ be the weak solution of \eqref{weak-ellip-equa}. If, for any $\epsilon>0$ there exists a small $\delta=\delta(n,\epsilon,\nu,\Lambda)$ such that \begin{gather*} Q_{8}^{+} \subset \Omega_{8} \subset Q\dot{}_{8} \cap \{x^1>-16\delta \},\quad -\hspace{-0.38cm}\int_{Q_8}|\theta(a,Q_8)|^2dx\leq \delta^2, \\ \{x\in \Omega_2:\mathcal{M}(|Du|^2)\leq 1\} \cap \{x\in \Omega_2:\mathcal{M}(\Psi^2)\leq \delta^2\}\neq\emptyset; \end{gather*} then there exists a positive constant $N_1=N_1(n,\nu,\Lambda)$ such that \begin{equation*} |\{x\in \Omega_2:\mathcal{M}(|Du|^2)>N_1^2\}\cap B_1| \leq \epsilon |B_1|. \end{equation*} \end{lemma} \begin{proof} This is very similar to the proof of Lemma 5.10 in \cite{ByK16} only replacing $|F|^2$ by $\Psi^2$, and we here omit its proof. \end{proof} As for the interior comparison estimates, we only state the results since it is simpler than that of the boundary case. By similar way we can also obtain the corresponding estimates as the above Lemma \ref{comparison-3} and \ref{boundary-maximal-function} without the boundary term. Without loss of generality we assume \begin{equation}\label{interior-condtion-a} Q_7 \subset B_{7\sqrt{2}} \subset \Omega. \end{equation} \begin{lemma}\label{interior-comparison-1} Let $u \in W^{1,2}(Q_6)$ be a local weak solution of \eqref{weak-ellip-equa} in $Q_6$ with the normalization of \begin{equation*} -\hspace{-0.38cm}\int_{Q_6}|Du|^2dx \leq 1, \quad -\hspace{-0.38cm}\int_{Q_6}\Psi^2dx \leq \delta^2, \end{equation*} and $w\in W^{1,2}(Q_5)$ be the weak solution of \eqref{comparison-equ-2} only replacing $\Omega_5$ by $Q_5$. If, for $0<\epsilon<1$ there exists a constant $\delta=\delta(n,\epsilon,\nu,\Lambda)>0$ such that $\mathbf{a}(\xi,x)$ satisfying $(\delta,R)$-vanishing of codimension 1 of \eqref{small-BMO-a}, then we have \begin{equation*} -\hspace{-0.38cm}\int_{Q_5}|Du-Dw|^2dx \leq \epsilon \quad \text{and} \quad \|Dw\|_{L^\infty(Q_4)}dx \leq c_2, \end{equation*} where $c_2=c_2(n,\epsilon,\nu,\Lambda)$. \end{lemma} \begin{lemma}\label{interior-maximal-function-a} Let $u\in W_0^{1,2}(\Omega)$ be a weak solution of \eqref{weak-ellip-equa}. If, for $\epsilon>0$ we find a small $\delta=\delta(n,\epsilon,\nu,\Lambda)>0$ such that $\mathbf{a}(\xi,x)$ satisfying $(\delta,R)$-vanishing of codimension 1 of \eqref{small-BMO-a} and \begin{equation*} \{x\in Q_2:\mathcal{M}(|Du|^2)\leq 1\}\cap \{x\in Q_2:\mathcal{M}(\Psi^2)\leq \delta^2\}\neq\emptyset, \end{equation*} then there exists a positive constant $N_2=N_2(n,\nu,\Lambda)$ such that \begin{equation*} |\{x\in Q_2:\mathcal{M}(|Du|^2)>N_2^2\}\cap B_1| \leq \epsilon |B_1|. \end{equation*} \end{lemma} Further, we obtain the next lemma by a scaling invariance to the above lemma, see also \cite[Lemma 5.2]{ByK16}. \begin{lemma}\label{interior-maximal-function-b} Let $u\in W_0^{1,2}(\Omega)$ be a weak solution of the variational inequalities \eqref{weak-ellip-equa}. If, for $\epsilon>0$ we find an $\delta=\delta(n,\epsilon,\nu,\Lambda)>0$ with $\mathbf{a}(\xi,x)$ satisfying $(\delta,160)$-vanishing of codimension 1 and there exists a positive constant $N_2=N_2(n,\nu,\Lambda)$ such that \begin{equation*} |\{x\in \Omega:\mathcal{M}(|Du|^2)>N_2^2\}\cap B_r(y)| \geq \epsilon |B_r(y)| \end{equation*} with $Q_{7r}(y)\subset B_{7\sqrt{2}r}(y)\subset \Omega$ for $0 1\} \cup\{x\in \Omega:\mathcal{M}(\Psi^2)> \delta^2\}. \end{equation*} \end{lemma} Now we write $ N_3=\max\{N_1,N_2,1\} $ with $N_1,N_2$ shown as in Lemma \ref{boundary-maximal-function} and \ref{interior-maximal-function-b}. By combining the interior estimate of Lemma \ref{interior-maximal-function-b} and the boundary estimate of Lemma \ref{boundary-maximal-function}, then we have the following estimates, cf. \cite[Lemmas 5.11 and 5.12]{ByK16}. \begin{lemma}\label{global-maximal-function-a} Let $u\in W_0^{1,2}(\Omega)$ be a weak solution of \eqref{weak-ellip-equa}. If for $\epsilon>0$ we can find a small $\delta=\delta(n,\epsilon,\nu,\Lambda)$ such that $(\mathbf{a},\Omega)$ satisfying $(\delta,160)$-vanishing of codimension 1 and \begin{equation*} |\{x\in \Omega:\mathcal{M}(|Du|^2)>N_3^2\}\cap B_r(y)| \geq \epsilon |B_r(y)| \end{equation*} for $y\in \Omega$ with $0 1\} \cup\{x\in \Omega:\mathcal{M}(\Psi^2)> \delta^2\}. \end{equation*} \end{lemma} By an iterating argument we conclude the following power decay estimate of the measure on the distribution concerning Hardy-Littlewood maximal functions. \begin{lemma}\label{global-maximal-function-b} Let the assumptions of Lemma \ref{global-maximal-function-a} hold and \begin{equation}\label{small} |\{x\in \Omega:\mathcal{M}(|Du|^2)>N_3^2\}| \leq \epsilon |B_r(y)|. \end{equation} Then \begin{align*} &|\{x\in \Omega:\mathcal{M}(|Du|^2)>N_3^{2k}\}| \\ &\leq \epsilon_1^k|\{x\in \Omega:\mathcal{M}(|Du|^2)> 1\}| +\sum_{i=1}^k\epsilon_1^i|\{x\in \Omega:\mathcal{M}(\Psi^2) > \delta^2N_3^{2(k-i)}\}| \end{align*} with $\epsilon_1=(\frac{10}{1-\delta})^n\epsilon$. \end{lemma} We also need the following measure equivalency to represent Orlicz spaces, which can be found in \cite{ByW12,RaR02}. \begin{lemma}\label{measure-equivalency} Let $f$ be a nonnegative measurable function in $U$, and the Young functions $\Phi\in\Delta_2\cap\nabla_2$. Then, for $\gamma>0$ and $m>1$, we have $f \in L^\Phi(U)$ if and only if \begin{equation*} S:=\sum_{k\geq 1} \Phi(m^k)|\{x\in U: f(x)>\gamma m^k\}|<\infty \end{equation*} and \begin{equation*} c^{-1}S\leq \int_{\mathbb{U}}\Phi(|f|)dx \leq c(S+|U|), \end{equation*} where the constant $c=c(\gamma,m,\Phi)>0$. \end{lemma} \section{Proof of Theorem \ref{main-result}} Note that $\Psi^2\in L^\Phi(\Omega)$ with $\Psi=|F|+|D\psi|$ for the Young function $\Phi\in \Delta_2\cap\nabla_2$. For any $q_0>1$, we use formula \eqref{mononicity-condition} with $t_1=t_2=q_0,\lambda=\Psi^2,\rho=1$ and the relation of equivalence \eqref{boundedness-equa-2}, H\"older inequality, Young inequality to get \begin{equation}\label{Psi-estimate-a} \int_{\Omega}\Psi^2\,dx \leq c\int_{\Omega}\Psi^{2q_0}dx+c|\Omega| \leq \frac{c}{\Phi(1)}\int_{\Omega}\Phi(\Psi^2)\,dx +c \leq c\big(\|\Psi^2\|^{\beta_0}_{L^\Phi(\Omega)}+1\big) \end{equation} with $c$ and $\beta_0$ depending only on $n,\Phi,|\Omega|$. According to Lemma \ref{scaling-normalization}, by a scaling argument it suffices to consider \begin{equation}\label{u-1} u_1=\frac{\delta u}{(\|\Psi^2\|_{L^\Phi(\Omega)})^{1/2}}, \quad \Psi_1=\frac{\delta \Psi}{(\|\Psi^2\|_{L^\Phi(\Omega)})^{1/2}}. \end{equation} Then, \eqref{Psi-estimate-a} yields \begin{equation}\label{Psi-estimate-b} \int_{\Omega}\Psi_1^2\,dx = \delta^2\int_{\Omega}\frac{\Psi^2}{\|\Psi^2\|_{L^\Phi(\Omega)}}\,dx \leq c\delta^2\Big(\Big\|\frac{\Psi^2}{\|\Psi^2\|_{L^\Phi(\Omega)}} \Big\|^{\beta_0}_{L^\Phi(\Omega)}+1\Big) \leq c\delta^2. \end{equation} To check condition \eqref{small}, we use the weak $(1,1)$-estimate \eqref{weak-(1,1)-estimate} on the Hardy-Littlewood maximal function, and the standard $L^2$-estimate \eqref{energy-estimate} on the variational inequalities \eqref{weak-ellip-equa} to obtain that \begin{equation}\label{small-check} \big|\big\{x\in \Omega:\mathcal{M}(|Du_1|^2)>N_3^2\big\}\big| \leq \frac {c(n)}{N^2_3}\int_{\Omega}|Du_1|^2\,dx \leq c\int_{\Omega}\Psi_1^2\,dx \leq c\delta^2 \leq \epsilon|B_1|, \end{equation} where we take $\delta>0$ sufficiently small so that the last inequality holds. Then, it follows from Lemma \ref{global-maximal-function-b} that \begin{equation}\label{Psi-estimate-b2} \begin{aligned} &\sum_{k=1}^\infty\Phi(N_3^{2k}) \big|\big\{x\in \Omega:\mathcal{M}(|Du_1|^2)>N_3^{2k}\big\}\big| \\ & \leq \sum_{k=1}^\infty\Phi(N_3^{2k})\Big(\epsilon_1^k \big|\big\{x\in \Omega:\mathcal{M}(|Du_1|^2)> 1\big\} \big| \\ &\quad +\sum_{i=1}^k\epsilon_1^i\big|\big\{x\in \Omega: \mathcal{M}(\Psi_1^2)> \delta^2N_3^{2(k-i)}\big\}\big|\Big) \\ & =\sum_{k=1}^\infty\Phi(N_3^{2k})\epsilon_1^k\big|\big\{x\in \Omega: \mathcal{M}(|Du_1|^2)> 1\big\}\big| \\ &\quad +\sum_{k=1}^\infty\Phi(N_3^{2k})\sum_{i=1}^k\epsilon_1^i \big|\big\{x\in \Omega:\mathcal{M}(\Psi_1^2)> \delta^2N_3^{2(k-i)}\big\}\big| \\ &:= I_1+I_2. \end{aligned} \end{equation} The condition $\Phi\in\Delta_2\cap\nabla_2$ implies $\Phi(N_3^2)\leq \mu\Phi(1)$ for some constant $\mu>1$ depending on $N_3^2$. Iterating this inequality, we obtain $\Phi(N_3^{2k})\leq \mu^k\Phi(1)$, then \begin{equation}\label{I-1} I_1\leq \Phi(1)|\Omega|\sum_{k=1}^\infty(\mu\epsilon_1)^k. \end{equation} Similarly, it follows from $\Phi(N_3^{2k})\leq \mu^i\Phi(N_3^{2(k-i)})$, the relation of equivalence \eqref{boundedness-equa-2}, Lemma \ref{boundedness} and Lemma \ref{measure-equivalency}, that \begin{equation}\label{I-2} \begin{aligned} I_2 &\leq \sum_{i=1}^\infty(\mu\epsilon_1)^i\sum_{k=i}^\infty\Phi(N_3^{2(k-i)}) \big|\big\{x\in \Omega:\mathcal{M}(\Psi_1^2)> \delta^2N_3^{2(k-i)}\big\}\big| \\ & = \sum_{i=1}^\infty(\mu\epsilon_1)^i\sum_{k=i}^\infty\Phi(N_3^{2(k-i)}) \big|\big\{x\in \Omega:\mathcal{M}(\frac{\Psi_1^2}{\delta^2})> N_3^{2(k-i)}\big\} \big| \\ &\leq c\sum_{i=1}^\infty(\mu\epsilon_1)^i\int_{\Omega}\Phi \Big(\mathcal{M}\Big(\frac{\Psi_1^2}{\delta^2}\Big)\Big)\,dx \\ &\leq c\sum_{i=1}^\infty(\mu\epsilon_1)^i\int_{\Omega}\Phi \Big(\frac{\Psi_1^2}{\delta^2}\Big)\,dx \\ & \leq c\sum_{i=1}^\infty(\mu\epsilon_1)^i \Big(\|\frac{\Psi_1^2}{\delta^2}\|_{L^\Phi(\Omega)}^{\beta_1}+1\Big) \\ & \leq c\sum_{i=1}^\infty(\mu\epsilon_1)^i, \end{aligned} \end{equation} where $1<\beta_1<\infty$ is a constant. Combining \eqref{Psi-estimate-b}, \eqref{I-1} and \eqref{I-2} together, then we obtain \begin{equation*} \sum_{k=1}^\infty\Phi(N_3^{2k}) \big|\big\{x\in \Omega: \mathcal{M}(|Du_1|^2)>N_3^{2k}\big\}\big| \leq c\sum_{i=1}^\infty(\mu\epsilon_1)^i\leq c, \end{equation*} where in the last inequality we take $\epsilon>0$ small enough such that $\mu\epsilon_1=\mu(\frac{10}{1-\delta})^n\epsilon\leq 1/2$. Then we find a corresponding $\delta>0$ such that Lemma \ref{global-maximal-function-a} and the estimate \eqref{small-check} hold. Finally, using Lemma \ref{boundedness} and Lemma \ref{measure-equivalency} with $\gamma=1,m=N_3^2$ again, it yields \begin{align*} \int_{\Omega}\Phi(|Du_1|^2)\,dx &\leq \int_{\Omega}\Phi(\mathcal{M}(|Du_1|^2))\,dx \\ &\leq \sum_{k=1}^\infty\Phi(N_3^{2k}) \big|\big\{x\in \Omega :\mathcal{M}(|Du_1|^2)>N_3^{2k}\big\}\big|+c \leq c. \end{align*} Furthermore, by \eqref{boundedness-equa-2}, we have \[ \||Du_1|^2\|_{L^\Phi(\Omega)} \leq \Big(c\int_{\Omega}\Phi(|Du_1|^2)+1\Big)^{1/\alpha_1}\leq c, \] where $1<\alpha_1<\infty$. By recalling the definition of $u_1$ in \eqref{u-1}, we obtain the desired estimate \eqref{Lorentz-estimate-a}. \subsection*{Acknowledgements} This research was supported by the National Natural Science Foundation of China grant No. 11371050. \begin{thebibliography}{20} \bibitem{AzBT00} Azroul, E. H.; Benkirane, A.; Tienari, M.; \emph{On the regularity of solutions to the Poisson equations in Orlicz spaces}, Bull Belg Math Soc (1)., \textbf{7} (2000), 1--12. \bibitem{ByW12} Byun, S. S.; Wang, L. H.; \emph{Nonlinear gradient estimates for elliptic equations of general type}, Calc. Var. Partial Differ. Equ., \textbf{45} (3-4) (2012), 403--419. \bibitem{ByP14} Byun, S. S.; Palagachev, D. K.; \emph{Weighted $L^p$-estimates for elliptic equations with measurable coefficients in nonsmooth domains}, Potential Anal., \textbf{41} (2014), 51--79. \bibitem{ByK16} Byun, S. S.; Kim, Y.; \emph{Elliptic equations with measurable nonlinearities in nonsmooth domains}, Adv. Math., \textbf{288} (2016), 152--200. \bibitem{ByOP16} Byun, S. S.; Ok, J.; Palagachev, D. K.; \emph{Parabolic systems with measurable coefficients in weighted Orlicz spaces}, Commun. Contemp. Math., \textbf{18} (2) (2016), 1550018, 19 pages. \bibitem{ChKV86} Chipot, M.; Kinderlehrer, D.; Vergara-Caffarelli, G.; \emph{Smoothness of linear laminates}, Arch. Ration. Mech. Anal., \textbf{96} (1986), 81--96. \bibitem{DoK10} Dong, H. J.; Kim, D.; \emph{Elliptic equations in divergence form with partially BMO coefficients}, Arch. Ration. Mech. Anal., \textbf{196} (2010), 25--70. \bibitem{DoK11-1} Dong, H. J.; Kim, D.; \emph{Parabolic and elliptic systems in divergence form with variably partially BMO coefficients}, SIAM J. Math. Anal., \textbf{43} (3) (2011), 1075--1098. \bibitem{DoK11-2} Dong, H. J.; Kim, D.; \emph{On the $L_{p}$-solvability of higher order parabolic and elliptic systems with BMO coefficients}, Arch. Ration. Mech. Anal., \textbf{199} (2011), 889--941. \bibitem{JiLW07} Jia, H. L.; Li, D. S.; Wang, L. H.; \emph{Regularity in Orlicz spaces for the Poisson equation}, Manuscripta Math., \textbf{122} (2007), 265--275. \bibitem{JiLW11} Jia, H. L.; Li, D. S.; Wang, L. H.; \emph{Global regularity for divergence form elliptic equations in Orlicz spaces on quasiconvex domains}, Nonlinear Anal., \textbf{74} (2011), 1336--1344. \bibitem{KiK07} Kim, D.; Krylov, N. V.; \emph{Elliptic differential equations with coefficients measurable with respect to one variable and VMO with respect to the others}, SIAM J. Math. Anal., \textbf{39} (2007), 489--506. \bibitem{KoK91} Kokilashvili, V.; Krbec, M.; \emph{Weighted inequalities in Lorentz and Orlicz spaces}, Singapore: World Scientific Publishing Co., 1991. \bibitem{LiZZ17} Li, H. Z.; Zhang, J. J.; Zheng, S. Z.; \emph{Orlicz estimates for nondivergence linear elliptic equations with partially BMO coefficients}, Complex Variables and Elliptic Equations, Published online: 28 Jul 2017, https://doi.org/ 10.1080/17476933.2017.1351960. \bibitem{PKJF13} Pick, L.; Kufner, A.; John, O.; Fu\v{c}\'{i}k, S.; Function Spaces, Walter de Gruyter GmbH, Berlin/Boston, 2013. \bibitem{RaR02} Rao, M. M.; Ren, Z. D.; \emph{Applications of Orlicz spaces}, New York (NY): Marcel Dekker, 2002. \bibitem{St93} Stein, E. M.; \emph{Harmonic Analysis: Real-Variable Methods, Orthogonality, and Oscillatory Integrals}, Monographs in Harmonic Analysis, III, Princeton Math. Ser., vol.43, Princeton University Press, Princeton, NJ, 1993. \bibitem{TiZ17} Tian, H.; Zheng, S. Z.; \emph{Uniformly nondegenerate elliptic equations with partially BMO coefficients in nonsmooth domains}, Nonlinear Anal., \textbf{156} (2017), 90--110. \bibitem{Ur67} Ural'tseva, N. N.; \emph{On impossibility of $W^{2,p}$ estimates for multidimensional elliptic equations with discontinuous coefficients}, Zap Nauchn Seminarov LOMI im Steklova, t. 5, ``Nauka," Leningrad (in Russian), 1967. \bibitem{YaSZ08} Yao, F. P.; Sun, Y.; Zhou, S. L.; \emph{Gradient estimates in Orlicz spaces for quasilinear elliptic equation}, Nonlinear Anal., \textbf{69} (2008), 2553--2565. \bibitem{YaZ12} Yao, F. P.; Zhou, S. L.; \emph{Global estimates in Orlicz spaces for $p$-Laplacian systems in $\mathbb{R}^N$}, J. Partial Differ. Equ., \textbf{25} (2) (2012), 1--12. \bibitem{ZhZ17} Zhang, J. J.; Zheng, S. Z.; \emph{Weighted Lorentz estimates for nondivergence linear elliptic equations with partially BMO coefficients}, Commun. Pure Appl. Anal., \textbf{16} (3) (2017), 899--914. \end{thebibliography} \end{document}