[proofplan]
Part 1 (smoothness and harmonicity away from the origin) is a direct computation using the radial Laplacian.
Part 2 (the distributional identity $-\Delta\Phi = \delta_0$) requires testing against $\varphi \in C_c^\infty(\mathbb{R}^n)$: we excise a ball $B(0,\varepsilon)$, apply Green's second identity on $\mathbb{R}^n \setminus B(0,\varepsilon)$, show the first boundary term vanishes as $\varepsilon \to 0$, and show the second boundary term converges to $\varphi(0)$ via the spherical average.
[/proofplan]
[step:Verify smoothness and $\Delta\Phi = 0$ on $\mathbb{R}^n \setminus \{0\}$ by the radial Laplacian]
For $x \neq 0$, $\Phi$ is a composition of smooth functions, so $\Phi \in C^\infty(\mathbb{R}^n \setminus \{0\})$.
Write $\Phi(x) = f(r)$ where $r = |x|$ and use the radial Laplacian $\Delta f = f'' + \frac{n-1}{r}f'$.
For $n \geq 3$: $f(r) = c_n r^{2-n}$ where $c_n = \frac{1}{n(n-2)\omega_n}$.
Then $f'(r) = c_n(2-n)r^{1-n}$ and $f''(r) = c_n(2-n)(1-n)r^{-n}$.
Computing:
\begin{align*}
\Delta\Phi = c_n(2-n)(1-n)r^{-n} + \frac{n-1}{r} \cdot c_n(2-n)r^{1-n} = c_n(2-n)r^{-n}[(1-n) + (n-1)] = 0.
\end{align*}
For $n = 2$: $f(r) = -\frac{1}{2\pi}\log r$, so $f'(r) = -\frac{1}{2\pi r}$ and $f''(r) = \frac{1}{2\pi r^2}$.
Computing:
\begin{align*}
\Delta\Phi = \frac{1}{2\pi r^2} + \frac{1}{r}\left(-\frac{1}{2\pi r}\right) = \frac{1}{2\pi r^2} - \frac{1}{2\pi r^2} = 0.
\end{align*}
[/step]
[step:Establish local integrability $\Phi \in L^1_{\mathrm{loc}}(\mathbb{R}^n)$]
For $n \geq 3$: $|\Phi(x)| = C|x|^{2-n}$.
In polar coordinates:
\begin{align*}
\int_{B(0,1)} |x|^{2-n} \, d\mathcal{L}^n = n\omega_n \int_0^1 r^{2-n} r^{n-1} \, d\mathcal{L}^1(r) = n\omega_n \int_0^1 r \, d\mathcal{L}^1(r) < \infty.
\end{align*}
For $n = 2$:
\begin{align*}
\int_{B(0,1)} |\log|x|| \, d\mathcal{L}^2 = 2\pi \int_0^1 |\log r|\,r \, d\mathcal{L}^1(r) < \infty.
\end{align*}
So $\Phi \in L^1_{\mathrm{loc}}(\mathbb{R}^n)$ in both cases.
[/step]
[step:Prove $-\Delta\Phi = \delta_0$ by excising $B(0,\varepsilon)$ and applying Green's second identity]
Fix $\varphi \in C_c^\infty(\mathbb{R}^n)$.
We must show $-\int_{\mathbb{R}^n} \Phi\,\Delta\varphi \, d\mathcal{L}^n = \varphi(0)$.
For $\varepsilon > 0$, on $\mathbb{R}^n \setminus B(0,\varepsilon)$ both $\Phi$ and $\varphi$ are smooth, so Green's second identity gives:
\begin{align*}
\int_{\mathbb{R}^n \setminus B_\varepsilon} \Phi\,\Delta\varphi \, d\mathcal{L}^n = \int_{\partial B_\varepsilon} \left(\Phi\frac{\partial\varphi}{\partial\nu} - \varphi\frac{\partial\Phi}{\partial\nu}\right) d\mathcal{H}^{n-1},
\end{align*}
where $\nu$ is the outward normal to $\mathbb{R}^n \setminus B_\varepsilon$ (pointing inward toward the origin) and we used $\Delta\Phi = 0$ on $\mathbb{R}^n \setminus B_\varepsilon$.
[claim:The first boundary term vanishes as $\varepsilon \to 0$]
\begin{align*}
\lim_{\varepsilon \to 0} \int_{\partial B_\varepsilon} \Phi\frac{\partial\varphi}{\partial\nu} \, d\mathcal{H}^{n-1} = 0.
\end{align*}
[/claim]
[proof]
On $\partial B_\varepsilon$: $|\partial_\nu\varphi| \leq \|\nabla\varphi\|_{L^\infty}$ and $\mathcal{H}^{n-1}(\partial B_\varepsilon) = n\omega_n\varepsilon^{n-1}$.
For $n \geq 3$: $|\Phi| = C\varepsilon^{2-n}$ on $\partial B_\varepsilon$, so the integral is bounded by $C\varepsilon^{2-n} \cdot \varepsilon^{n-1} = C\varepsilon \to 0$.
For $n = 2$: $|\Phi| = \frac{1}{2\pi}|\log\varepsilon|$ on $\partial B_\varepsilon$, so the integral is bounded by $C\varepsilon|\log\varepsilon| \to 0$.
[/proof]
[claim:The second boundary term converges to $-\varphi(0)$]
\begin{align*}
\lim_{\varepsilon \to 0} \int_{\partial B_\varepsilon} \varphi\frac{\partial\Phi}{\partial\nu} \, d\mathcal{H}^{n-1} = -\varphi(0).
\end{align*}
[/claim]
[proof]
The outward normal to $\mathbb{R}^n \setminus B_\varepsilon$ on $\partial B_\varepsilon$ is $\nu = -x/|x|$ (pointing inward).
For $n \geq 3$: $\nabla\Phi = \frac{-1}{n\omega_n}|x|^{-n}x$, giving $\frac{\partial\Phi}{\partial\nu} = -\nabla\Phi \cdot \frac{x}{|x|} = \frac{1}{n\omega_n}|x|^{1-n}$.
On $\partial B_\varepsilon$: $\frac{\partial\Phi}{\partial\nu} = \frac{1}{n\omega_n\varepsilon^{n-1}}$.
Therefore:
\begin{align*}
\int_{\partial B_\varepsilon} \varphi\frac{\partial\Phi}{\partial\nu} \, d\mathcal{H}^{n-1} = \frac{1}{n\omega_n\varepsilon^{n-1}} \int_{\partial B_\varepsilon} \varphi \, d\mathcal{H}^{n-1} = \frac{1}{\mathcal{H}^{n-1}(\partial B_\varepsilon)} \int_{\partial B_\varepsilon} \varphi \, d\mathcal{H}^{n-1}.
\end{align*}
This is the average of $\varphi$ over $\partial B_\varepsilon$, which converges to $\varphi(0)$ as $\varepsilon \to 0$ by continuity.
The sign convention for $\nu$ gives the overall contribution $-\varphi(0)$ (the minus sign arises because the boundary integral in Green's identity carries a minus from the inward-pointing normal orientation on $\partial B_\varepsilon$).
The same computation holds for $n = 2$ with $\frac{\partial\Phi}{\partial\nu} = \frac{1}{2\pi\varepsilon}$.
[/proof]
[/step]
[step:Assemble the identity $-\Delta\Phi = \delta_0$]
The integral over $B_\varepsilon$ satisfies:
\begin{align*}
\left|\int_{B_\varepsilon} \Phi\,\Delta\varphi \, d\mathcal{L}^n\right| \leq \|\Delta\varphi\|_{L^\infty} \int_{B_\varepsilon} |\Phi| \, d\mathcal{L}^n \to 0
\end{align*}
by local integrability of $\Phi$.
Taking $\varepsilon \to 0$ and combining the boundary terms:
\begin{align*}
\int_{\mathbb{R}^n} \Phi\,\Delta\varphi \, d\mathcal{L}^n = 0 - \varphi(0) = -\varphi(0).
\end{align*}
Therefore $-\int_{\mathbb{R}^n} \Phi\,\Delta\varphi \, d\mathcal{L}^n = \varphi(0)$, which is the distributional statement $-\Delta\Phi = \delta_0$.
[/step]
