[proofplan]
We construct the Green's function for the ball using Kelvin inversion (the method of images): for each $x \in B(0,r)$, the reflected point $\tilde{x} = r^2 x / |x|^2$ lies outside $\overline{B(0,r)}$, and a corrector built from $\Phi(\cdot - \tilde{x})$ makes the Green's function vanish on $\partial B(0,r)$.
The Poisson kernel is obtained as the normal derivative of the Green's function on the boundary.
We verify that the resulting integral formula defines a harmonic function with the correct boundary values by checking the three properties of an approximate identity: positivity, unit mass, and concentration.
[/proofplan]
[step:Construct the Green's function via Kelvin inversion]
For $x \in B(0,r) \setminus \{0\}$, define the Kelvin reflection $\tilde{x} := \frac{r^2}{|x|^2}x$, which maps $B(0,r) \setminus \{0\}$ to $\mathbb{R}^n \setminus \overline{B(0,r)}$.
The Green's function is:
\begin{align*}
G(x, y) := \Phi(y - x) - \Phi\!\left(\frac{|x|}{r}(y - \tilde{x})\right).
\end{align*}
The corrector $\phi^x(y) := \Phi(|x|(y - \tilde{x})/r)$ is harmonic in $y$ on $B(0,r)$ since $\tilde{x} \notin \overline{B(0,r)}$.
For $|y| = r$, the identity $|y - \tilde{x}| = \frac{r}{|x|}|y - x|$ holds:
\begin{align*}
|y - \tilde{x}|^2 = |y|^2 - 2y \cdot \tilde{x} + |\tilde{x}|^2 = r^2 - \frac{2r^2}{|x|^2}y \cdot x + \frac{r^4}{|x|^2} = \frac{r^2}{|x|^2}(|x|^2 - 2y \cdot x + r^2) = \frac{r^2}{|x|^2}|y-x|^2.
\end{align*}
Therefore $\Phi(y - x) = \Phi(|x|(y - \tilde{x})/r)$ when $|y| = r$, giving $G(x, y) = 0$ on $\partial B(0,r)$.
[/step]
[step:Derive the Poisson kernel as the normal derivative of $G$]
The Poisson kernel is $K(x, y) := -\frac{\partial G}{\partial\nu_y}(x, y)$ for $y \in \partial B(0,r)$, where $\nu = y/r$ is the outward normal.
Computing $\frac{\partial}{\partial\nu_y}\Phi(y - x)$ on $\partial B(0,r)$ and subtracting the corrector contribution (using the identity $|y - \tilde{x}| = \frac{r}{|x|}|y - x|$ from the previous step), the simplification yields:
\begin{align*}
K(x, y) = \frac{r^2 - |x|^2}{n\omega_n r\,|y - x|^n} \quad \text{for } x \in B(0,r),\, y \in \partial B(0,r).
\end{align*}
[/step]
[step:Verify the Poisson kernel properties and the boundary behaviour]
For $y \in \partial B(0,r)$ and $x \in B(0,r)$: $K(x, y) > 0$ since $|x| < r$ and $|y - x| > 0$.
The unit mass property $\int_{\partial B(0,r)} K(x, y) \, d\mathcal{H}^{n-1}(y) = 1$ follows by applying the representation formula to $u \equiv 1$, which is harmonic with boundary data $g \equiv 1$.
As $x \to y_0 \in \partial B(0,r)$, the numerator $r^2 - |x|^2 \to 0$ while $|y - x|^n$ remains bounded below for $y$ away from $y_0$.
This concentrates $K(x, y)$ at $y = y_0$, so $K$ behaves as an approximate identity.
Define $u(x) := \int_{\partial B(0,r)} K(x, y)\,g(y) \, d\mathcal{H}^{n-1}(y)$.
The function $u$ is harmonic in $B(0,r)$: differentiation under the integral sign is justified by the smoothness of $K$ in $x$ for $x$ bounded away from $\partial B$, and $\Delta_x K = 0$ because $G$ is harmonic in $x$ for $x \neq y$.
For the boundary limit: as $x \to x_0 \in \partial B(0,r)$, split $\int_{\partial B} = \int_{|y - x_0| < \delta} + \int_{|y - x_0| \geq \delta}$.
On the far piece, $K(x, y) \to 0$ uniformly.
On the near piece, $|g(y) - g(x_0)| < \varepsilon$ by continuity.
Using $\int K = 1$ and the standard approximate identity argument gives $u(x) \to g(x_0)$.
[/step]