According to Theorem 1.1' in this paper we have the following estimate on classical solutions $u \in C^2(\overline{B_1^+})$ of $-\Delta u = f \text{ in } B_1^+ = B_1 \cap \{x _n \ge 0 \}$ and $u = 0 \text{ on } \partial B_1^+ \cap \{x_n = 0\}$
$$|D^2u(x) - D^2u(y)| \le C\left(r\lVert u \rVert_{L^\infty(B_{1}^+)} + \int_0^{r} \frac{\omega_f(t)}{t}\,dt + r\int_{r}^{1} \frac{\omega_f(t)}{t^2}\,dt\right) \tag{1}$$ $\forall \, x,y \in B_{1/2}^{+}$ with $r = |x-y|$ where, $\omega_f$ denotes the modulus of continuity of $f$ which we assume is Dini continuous.
The proof is supposed to be similar to proof of Theorem 1.1 which is the interior estimate. Following the lines of this proof we can get all analogous 'boundary' harmonic function estimates in $B_{1}^+$ v.i.a. applying the estimates in Theorem 1.1 to the odd extensions of the harmonic functions to the full ball $B_1$. But in the last step (going along the lines of eqn $(1.13)$ in the paper) it seems we need the estimate $$|D^2u_0(x) - D^2u_0(y)| \le C \lVert u \rVert_{L^{\infty}(B_1^+)}|x-y|, \, \forall \, x,y \in B_{1/2}^+ \tag{2}$$ where, $u_0$ satisfies $-\Delta u_0 = f(0)$ in $B_{1}^+$ and $u_0 = 0$ on $\partial B_1^+ \cap \{x_n = 0\}$ where, $C$ is supposed to be independent of $f$.
I can't seem to get the estimate $(2)$ for half-ball.
In case of interior estimate we can consider $v_0 := u_0 - \frac{f(0)}{2n}(1 - |x|^2)$ which is harmonic in $B_1$ and write \begin{align*} |D^2u_0(x) - D^2u_0(y)| = |D^2v_0(x) - D^2v_0(y)| &\le r\lVert D^3v_0\rVert_{L^{\infty}(B_{1/2})} \\ &\le r\lVert v_0\rVert_{L^{\infty}(\partial B_{1})} = r\lVert u_0\rVert_{L^{\infty}(\partial B_{1})}\end{align*}\begin{align*} |D^2u_0(x) - D^2u_0(y)| = |D^2v_0(x) - D^2v_0(y)| &\le r\lVert D^3v_0\rVert_{L^{\infty}(B_{1/2})} \\ &\le Cr\lVert v_0\rVert_{L^{\infty}(\partial B_{1})} = Cr\lVert u_0\rVert_{L^{\infty}(\partial B_{1})}\end{align*} where, $r = |x-y|$ which proves the interior analogue of $(2)$ using only gradient estimate for harmonic function $v_0$.
But a similar approach doesn't seem to be working for the boundary case (for example, considering $u_0 - \frac{f(0)}{2}x_n^2$ which is harmonic and applying odd extension to this.)
Is there a different way of approaching the estimate $(2)$?
Any help is appreciated. Thanks.