I'm reading the Theorem2 in UNIFORM ESTIMATES AND BLOW-UP BEHAVIOR FOR SOLUTIONS OF $-\Delta u=V(x) e^u$ IN TWO DIMENSIONS They prove that for the solution of $$ -\Delta u= V(x)\exp u \text { in } \mathbb{R}^2 $$ if $V \in L^{p}(\mathbb{R}^2)$, $\exp u \in L^{p^{\prime}}\left(\mathbb{R}^2\right)$, with $1<p\le+\infty$, then $u \in L^{\infty}\left(\mathbb{R}^2\right)$. (Here $1/p + 1/p^{\prime}=1$)
First fix $0<\epsilon<1 $, note that for any $1 \leq p<q<r \leq \infty, L^q \subseteq L^p+L^r$, so they split $V(x)\exp u $ as $V(x)\exp u=f_1+f_2$ with $\left\|f_1\right\|_{L^1\left(\mathbb{R}^2\right)}<\epsilon$ and $f_2 \in L^{\infty}\left(\mathbb{R}^2\right)$. Let $B_r$ be the ball of radius $r$ centered at $x_0$.
They denote by $C$ various constants independent of $x_0$ (but possibly depending on $\epsilon$ ). Then let $u_i$ be the solution of $$ \left\{\begin{array}{rlll} -\Delta u_i & =f_i & \text { in } & B_1 \\ u_i & =0 & \text { on } & \partial B_1 \end{array}\right. $$
By a lemma they proved before, they proved that
Assume $\Omega \subset \mathbb{R}^2$ is a bounded domain and let $u$ be a solution of
$$ \left\{\begin{array}{ccc} -\Delta u=f(x) & \text { in } \quad \Omega, \\ u=0 & \text { on } \quad \partial \Omega, \end{array}\right. $$ Then for every $\delta \in(0,4 \pi)$ we have $$ \int_{\Omega} \exp \left[\frac{(4 \pi-\delta)] u(x) \|}{\|f\|_1}\right] \mathrm{dx} \leq \frac{4 \pi^2}{\delta}(\operatorname{diam} \Omega)^2 . $$
So applied with $\delta=4 \pi-1$, they have $$ \int_{B_1} \exp \left[\frac{1}{\epsilon}\left|u_1\right|\right] \leq \mathrm{C} $$ and in particular $\left\|u_1\right\|_{L^1\left(B_1\right)} \leq C$. We also have $\left\|u_2\right\|_{L^{\infty}\left(B_1\right)} \leq C$. Let $u_3=u-u_1-u_2$ so that $\Delta u_3=0$ on $B_1$. The mean value theorem for harmonic functions implies that \begin{equation}\tag{1} \left\|u_3^{+}\right\|_{L^{\infty}\left(B_{1 / 2}\right)} \leq \mathrm{C}\left\|u_3^{+}\right\|_{L^1\left(B_1\right)} \end{equation} (This can be simply verified by the fact that let $\Omega \subset \mathbb{R}^n$ open and let $u$ be a harmonic function in $\Omega$. If $K \subset \Omega$ is compact, then $$ \sup _{x \in K}|u(x)| \leq \frac{n}{\omega_n \operatorname{dist}(K, \partial \Omega)^n} \int_{\Omega}|u(x)| \mathrm{d} x .) $$
On the other hand they have $$ u_3^{+} \leq u^{+}+\left|u_1\right|+\left|u_2\right| $$ and since $$ p^{\prime}\int_{R^2} u^{+}\mathrm{d} x \leq \int_{R^2} \exp p^{\prime}u\mathrm{d} x \leq C $$ we see that $\left\|u_3^{+}\right\|_{L^1\left(B_1\right)} \leq \mathrm{C}$. Combining this with (1) they find that $$\left\|u_3^{+}\right\|_{L^{\infty}\left(B_{1 / 2}\right)} \leq C.$$ Finally they write \begin{equation}\tag{2} -\Delta u= V(x)\exp u= V(x)\exp u_1 \exp (u_2+u_3)=g \end{equation} (with $V \in L^p(B_1)$, $\|g\|_{L^{1+\delta}\left(B_{1 / 2}\right)} \leq C$ for some $\delta>0$ (since $e^{u_2+u_3} \in L^{\infty}\left(B_{1 / 2}\right)$, $\exp u_1 \in L^{1 / \epsilon}\left(B_1\right)$ with $\left.1 / \epsilon> p^{\prime}\right)$.
What I'm confused is the last step, they said they use once more the mean value theorem and standard elliptic estimates they deduce from (2) that
$$ \left\|u^{+}\right\|_{L^{\mathbb{\infty}}\left(B_{1 / 4}\right)} \leq C\left\|u^{+}\right\|_{L^1\left(B_{1 / 2}\right)}+C \|u\|_{L^{1+\delta}\left(B_{1 / 2}\right)} \leq C . $$ Since $C$ is independent of $x_0$ they conclude that $u^{+} \in L^{\infty}\left(\mathbb{R}^2\right)$.
I'm very confused about how they use the mean value property again and what kind of regularity they use ?
