(a)
\begin{equation}
u(y) \leq \frac{1}{\sigma_n r^{n-1}} \int_{S(y,r)} u(x) \, d S
\end{equation}
In $\Delta u \geqq 0$, $u$ is twice differentiable, so the spherical mean is continuous and we have:
\begin{equation}
\lim_{r \rightarrow 0^+} \frac{1}{\sigma_n r^{n-1}} \int_{S(y,r)} u(x) \, d S \rightarrow u(y) \label{2.a.1}
\end{equation}
Thus, to prove $ u(y) \leqq \frac{1}{\sigma_n r^{n-1}} \int_{S(y,r)} u(x) \, d S $, it is sufficient to show that $\frac{1}{\sigma_n r^{n-1}} \int_{S(y,r)} u(x) \, d S$ increases with $r > 0$. By the divergence theorem we have for a $C^1$ function $F$ on a compact set $\Omega$ with boundary $\partial\Omega = \Sigma$:
$$ \int_\Omega (\nabla \cdot F) \,dV = \int_\Sigma (F \cdot n) \, dS \text{, namely: } $$
$$\int_{B(y,r)} (\nabla \cdot \nabla )(x) \,dV = \int_{B(y,r)} \Delta u(x) \,dV = \int_{S(y,r)} (\nabla u \cdot n)(x) \, dS = \int_{S(y,r)} \frac{\partial u}{\partial n}(x) \, dS $$
As $\Delta u \geqq 0$, we then have:
$$ 0 \leqq \int_{B(y,r)} \Delta u(x) \,dV = \int_{S(y,r)} \frac{\partial u}{\partial n}(x) \, dS $$
Parametrizing by the normal $\nu = \frac{x-y}{|x-y|}$ on the boundary where $|x-y| = r$, we have: $\frac{\partial u}{\partial n} (x) = \partial_r u(y+ r \nu)$ on $ T : \{ \nu : | \nu | = 1 \} $, and $ dS = r^{n-1} d \nu $
Then:
$$ \int_{S(y,r)} \frac{\partial u}{\partial n}(x) \, dS = \int_{T} \partial_r( u(y+ r \nu)) r^{n-1} \, d \nu $$
By differentiating under the sign with $\partial_r$ and pulling the $r^{n-1}$ term out of the integral, we have:
$$ = r^{n-1} \partial_r \int_{T} u(y+ r \nu) \, d \nu = r^{n-1} \partial_r (r^{1-n} \int_{S(y,r)} u(x) \, d S) $$
$$ = r^{n-1} \partial_r(\sigma_n \frac{r^{1-n}}{\sigma_n } \int_{S(y,r)} u(x) \, d S ) = r^{n-1} \sigma_n \partial_r( \frac{1}{\sigma_n r^{n-1} } \int_{S(y,r)} u(x) \, d S ) $$
As $r>0$ and $\sigma_n>0$ by (2) we have:
$$ 0 \leqq \int_{S(y,r)} \frac{\partial u}{\partial n}(x) \, dS = r^{n-1} \sigma_n \partial_r( \frac{1}{\sigma_n r^{n-1} } \int_{S(y,r)} u(x) \, d S ) \implies $$
\begin{equation}
0 \leqq \partial_r( \frac{1}{\sigma_n r^{n-1} } \int_{S(y,r)} u(x) \, d S )
\end{equation}
Then by (2) as $r \rightarrow 0$ we have equality in (1), and as $r$ increases the spherical mean is non-decreasing by (3), $u(y)$ constant. Thus we have (1) for $r \geqq 0$:
$$ u(y) \leqq \frac{1}{\sigma_n r^{n-1}} \int_{S(y,r)} u(x) \, d S \phantom{O} $$
(b)
$$u(y) = \frac{\sigma_n (n r^n)}{(n \sigma_n) r^n} u(y) = \frac{\sigma_n \int_{0}^{r} t^{n-1} \, dt}{\omega_n r^n} u(y) $$
By (1) we have that:
$$ u(y) \leqq \frac{1}{\sigma_n r^{n-1}} \int_{S(y,r)} u(x) \, d S $$
So because $ \{r, \omega_n, \sigma_n\} \geqq 0 $ we have that:
$$\implies \frac{\sigma_n \int_{0}^{r} t^{n-1} \, dt}{\omega_n r^n} u(y) \leqq \frac{\sigma_n \int_{0}^{r} t^{n-1}\, dt}{\omega_n r^n} \frac{1}{\sigma_n r^{n-1}} \int_{S(y,r)} u(x) \, d S $$
$$ = \frac{1}{\omega_n r^n} \sigma_n\int_{0}^{r}t^{n-1} \,d t \frac{1}{\sigma_n r^{n-1}} \int_T u(y +r \nu) r^{n-1} \, d \nu $$
Where we parameterized as in part a. with $\nu = \frac{x-y}{|x-y|}$ on the boundary $|x-y| = r$, and $ (x) = u(y+ r \nu)$ on $ T : \{ \nu : | \nu | = 1 \} $, $ dS = r^{n-1} d \nu $. Canceling terms gives:
$$ = \frac{1}{\omega_n r^n} \int_{0}^{r}t^{n-1} \,d t \int_T u(y +r \nu) \, d \nu $$
Which allows us to change the order of integration to yield our result:
$$ = \frac{1}{\omega_n r^n} \int_{0}^{r} t^{n-1} \int_T u(y +r \nu) \, d \nu d t = \frac{1}{\omega_n r^{n}} \int_{B(y,r)} u \, d V $$
$$ \implies u(y) \leqq \frac{1}{\omega_n r^{n}} \int_{B(y,r)} u \, d V \phantom{O} $$
(c)
In part a. we used the fact that $ \Delta u \geqq 0 $, and $ \{ \sigma_n, r \} \geqq 0$ to show that the spherical mean is non-decreasing in $r$:
$$ 0 \leqq \int_{B(y,r)} \Delta u(x) \,dV = r^{n-1} \sigma_n \partial_r( \frac{1}{\sigma_n r^{n-1} } \int_{S(y,r)} u(x) \, d S )$$
$$ \implies 0 \leqq \partial_r( \frac{1}{\sigma_n r^{n-1} } \int_{S(y,r)} u(x) \, d S ) $$
Because we have equality at $r \rightarrow 0$:
$$ \lim_{r \rightarrow 0^+} \frac{1}{\sigma_n r^{n-1}} \int_{S(y,r)} u(x) \, d S \rightarrow u(y) $$
and $u(y)$ is constant with respect to $r$, we showed that for all $r > 0$
$$ u(y) \leqq \frac{1}{\sigma_n r^{n-1}} \int_{S(y,r)} u(x) \, d S $$
If instead $ \Delta u \leqq 0$, we would have the inequalities:
$$ 0 \geqq \int_{B(y,r)} \Delta u(x) \,dV = r^{n-1} \sigma_n \partial_r( \frac{1}{\sigma_n r^{n-1} } \int_{S(y,r)} u(x) \, d S )$$
$$ \implies 0 \geqq \partial_r( \frac{1}{\sigma_n r^{n-1} } \int_{S(y,r)} u(x) \, d S ) $$
And $ \frac{1}{\sigma_n r^{n-1} } \int_{S(y,r)} u(x) \, d S$ is non-increasing with $r$, which gives us by the same argument that for all $r>0$
$$ u(y) \geqq \frac{1}{\sigma_n r^{n-1}} \int_{S(y,r)} u(x) \, d S \phantom{O} \square$$
In part b. we merely integrated this result radially to extend the inequality for the interior of the sphere. Doing the same in the case of $ \Delta u \leqq 0$ yields us:
$$ u(y) \geqq \frac{1}{\omega_n r^{n}} \int_{B(y,r)} u \, d V \phantom{O} $$
