Idea of the proof: We interpret $\Phi$ as a distribution by pairing it with test functions and compute its Laplacian via duality. We split the domain near the singularity and use Green's identity to reduce the problem to a boundary integral over spheres of radius $\varepsilon$, then pass to the limit using dominated convergence.

Let $\varphi \in C_c^\infty(\mathbb{R}^n)$. Define the distribution $\Delta \Phi$ by duality: \begin{align*} \langle \Delta \Phi, \varphi \rangle := \langle \Phi, \Delta \varphi \rangle = \int_{\mathbb{R}^n} \Phi(x) \Delta \varphi(x) \, dx. \end{align*} We fix $R > 0$ such that $\operatorname{supp}(\varphi) \subseteq B(0, R)$. For $0 < \varepsilon < R$, we write: \begin{align*} \int_{\mathbb{R}^n} \Phi(x) \Delta \varphi(x) \, dx = \int_{B(0, \varepsilon)} \Phi(x) \Delta \varphi(x) \, dx + \int_{B(0, R) \setminus B(0, \varepsilon)} \Phi(x) \Delta \varphi(x) \, dx. \end{align*}

**Step 1: Limit of the integral over $B(0, \varepsilon)$ is zero**

Let $f_\varepsilon := \mathbf{1}_{B(0, \varepsilon)} \cdot \Phi \cdot \Delta \varphi$. Since $\Phi \in L^1_{\mathrm{loc}}(\mathbb{R}^n)$ and $\Delta \varphi$ is bounded, $f_\varepsilon$ is integrable. Also, for almost every $x \in B(0, R)$, \begin{align*} \lim_{\varepsilon \to 0} f_\varepsilon(x) = 0. \end{align*} We have: \begin{align*} |f_\varepsilon(x)| \leq \|\Delta \varphi\|_{L^\infty(B(0,R))} \cdot |\Phi(x)| \cdot \mathbf{1}_{B(0, \varepsilon)}(x), \end{align*} and the right-hand side is integrable on $B(0,R)$ by local integrability of $\Phi$. Hence, by the **dominated convergence theorem**: \begin{align*} \lim_{\varepsilon \to 0} \int_{B(0, \varepsilon)} \Phi(x) \Delta \varphi(x) \, dx = 0. \end{align*}

**Step 2: Use of Green’s identity on the annulus**

Since $\Phi \in C^\infty(\mathbb{R}^n \setminus \{0\})$ and $\Delta \Phi = 0$ on $\mathbb{R}^n \setminus \{0\}$, we may use Green's identity on the open set $B(0,R) \setminus \overline{B(0,\varepsilon)}$: \begin{align*} \int_{B(0, R) \setminus B(0, \varepsilon)} \Phi \Delta \varphi \, dx = \int_{\partial B(0, \varepsilon)} \left( \Phi \frac{\partial \varphi}{\partial \nu} - \varphi \frac{\partial \Phi}{\partial \nu} \right) \, d\mathcal{H}^{n-1}. \end{align*} Here, $\nu$ denotes the outward unit normal and $\mathcal{H}^{n-1}$ is $(n-1)$-dimensional Hausdorff measure on the sphere $\partial B(0, \varepsilon)$. Since $\varphi$ has compact support in $B(0,R)$, the outer boundary term at $\partial B(0,R)$ vanishes.

Let us evaluate the two terms separately.

**Step 3: Evaluation of $\int_{\partial B(0,\varepsilon)} \varphi \frac{\partial \Phi}{\partial \nu} \, d\mathcal{H}^{n-1}$**

Let $x = \varepsilon \omega$ for $\omega \in \partial B(0,1)$. Then $|x| = \varepsilon$, and since $\nu(x) = \omega$, we have: - For $n = 2$: \begin{align*} \frac{\partial \Phi}{\partial \nu}(x) = \nabla \Phi(x) \cdot \nu(x) = -\frac{1}{2\pi} \cdot \frac{1}{\varepsilon}. \end{align*} - For $n \geq 3$: \begin{align*} \frac{\partial \Phi}{\partial \nu}(x) = \nabla \Phi(x) \cdot \nu(x) = -\frac{1}{\alpha(n)} \cdot \frac{1}{\varepsilon^{n-1}}. \end{align*}

Next, write: \begin{align*} \int_{\partial B(0, \varepsilon)} \varphi(x) \frac{\partial \Phi}{\partial \nu}(x) \, d\mathcal{H}^{n-1}(x) &= \frac{\partial \Phi}{\partial \nu}(\varepsilon \omega) \cdot \int_{\partial B(0, \varepsilon)} \varphi(x) \, d\mathcal{H}^{n-1}(x). \end{align*} We expand $\varphi$ in Taylor series: \begin{align*} \varphi(\varepsilon \omega) = \varphi(0) + \varepsilon D\varphi(0) \cdot \omega + \varepsilon^2 R_\varepsilon(\omega), \end{align*} where $R_\varepsilon(\omega)$ is bounded uniformly in $\omega$ as $\varepsilon \to 0$. Since the integral of $\omega$ over the unit sphere vanishes, we obtain: \begin{align*} \int_{\partial B(0, \varepsilon)} \varphi(x) \, d\mathcal{H}^{n-1}(x) &= \varphi(0) \cdot \mathcal{H}^{n-1}(\partial B(0, \varepsilon)) + o(\varepsilon^{n-1}). \end{align*}

Now, \begin{align*} \mathcal{H}^{n-1}(\partial B(0, \varepsilon)) = n \alpha(n) \varepsilon^{n - 1}. \end{align*}

Hence: - For $n = 2$: \begin{align*} \int_{\partial B(0,\varepsilon)} \varphi(x) \frac{\partial \Phi}{\partial \nu}(x) \, d\mathcal{H}^{1}(x) &= -\frac{1}{2\pi \varepsilon} \cdot (2\pi \varepsilon \varphi(0) + o(\varepsilon)) = -\varphi(0) + o(1). \end{align*} - For $n \geq 3$: \begin{align*} \int_{\partial B(0,\varepsilon)} \varphi(x) \frac{\partial \Phi}{\partial \nu}(x) \, d\mathcal{H}^{n-1}(x) &= -\frac{1}{\alpha(n) \varepsilon^{n-1}} \cdot (n \alpha(n) \varepsilon^{n - 1} \varphi(0) + o(\varepsilon^{n-1})) = -n \varphi(0) + o(1). \end{align*}

But recall that the coefficient in the definition of $\Phi$ for $n \geq 3$ is normalized so that this integral equals $-\varphi(0)$ (by the computation of the normal derivative). Hence, in all dimensions $n \geq 2$, \begin{align*} \int_{\partial B(0, \varepsilon)} \varphi(x) \frac{\partial \Phi}{\partial \nu}(x) \, d\mathcal{H}^{n-1}(x) = -\varphi(0) + o(1). \end{align*}

**Step 4: Estimate of $\int_{\partial B(0, \varepsilon)} \Phi \frac{\partial \varphi}{\partial \nu} \, d\mathcal{H}^{n-1}$**

We bound this term: \begin{align*} \left| \int_{\partial B(0, \varepsilon)} \Phi(x) \frac{\partial \varphi}{\partial \nu}(x) \, d\mathcal{H}^{n-1}(x) \right| &\leq \sup_{x \in B(0, R)} \left| \nabla \varphi(x) \right| \cdot \int_{\partial B(0, \varepsilon)} |\Phi(x)| \, d\mathcal{H}^{n-1}(x). \end{align*}

For $n = 2$, $\Phi(x) = -\frac{1}{2\pi} \log \varepsilon$, so: \begin{align*} \int_{\partial B(0, \varepsilon)} |\Phi(x)| \, d\mathcal{H}^{1}(x) = \frac{1}{2\pi} |\log \varepsilon| \cdot 2\pi \varepsilon = \varepsilon |\log \varepsilon| \to 0 \quad \text{as } \varepsilon \to 0. \end{align*}

For $n \geq 3$, $\Phi(x) = \frac{1}{n(n - 2)\alpha(n)} \cdot \varepsilon^{-(n-2)}$, and: \begin{align*} \int_{\partial B(0, \varepsilon)} \Phi(x) \, d\mathcal{H}^{n-1}(x) &= \Phi(\varepsilon) \cdot \mathcal{H}^{n-1}(\partial B(0, \varepsilon)) = \frac{1}{n(n-2)\alpha(n)} \cdot \varepsilon^{-(n-2)} \cdot n \alpha(n) \varepsilon^{n-1} = \frac{\varepsilon}{n-2} \to 0. \end{align*}

Thus in both cases, \begin{align*} \int_{\partial B(0,\varepsilon)} \Phi(x) \frac{\partial \varphi}{\partial \nu}(x) \, d\mathcal{H}^{n-1}(x) \to 0. \end{align*}

**Step 5: Conclusion**

Combining all steps: \begin{align*} \langle \Delta \Phi, \varphi \rangle &= \lim_{\varepsilon \to 0} \left( \int_{\partial B(0, \varepsilon)} \Phi \frac{\partial \varphi}{\partial \nu} \, d\mathcal{H}^{n-1} - \int_{\partial B(0, \varepsilon)} \varphi \frac{\partial \Phi}{\partial \nu} \, d\mathcal{H}^{n-1} + \int_{B(0, \varepsilon)} \Phi \Delta \varphi \, dx \right) \\ &= -\varphi(0). \end{align*}

Therefore, \begin{align*} \Delta \Phi = - \delta_0 \quad \text{in } \mathcal{D}'(\mathbb{R}^n). \end{align*}