Explanation of the identities $\rho=\rho(\delta)$ ($\delta$ function) and $\rho=q\,\delta(\bar{r})$

Question

Poisson's equation is $$\boldsymbol{\nabla}^{2}\varphi=-4\pi k_{e}\,\rho, \tag{*}$$ that in the case of a point charge $q$, already with spherical symmetry, has as solution

\begin{equation} \varphi(r)=k_{e} \frac{q}{r}, \tag{**} \end{equation}

Replacing the (**) in the (*) we get:

\begin{equation} k_{e}\, q \boldsymbol{\nabla}_{r}^{2}\left(\frac{1}{r}\right)=-4\pi \, k_{e}\, \rho \end{equation}

1. Why is the charge density $\rho$ also considered as a delta-function ($\delta$) over the whole classical space-time $\mathbb{R}^4$?
2. Why exist this identity

\begin{equation} \rho=q\,\delta(\bar{r})\quad ? \end{equation}

Answer 1

As you noticed the Laplacian is ill-defined at the origin: such equations are usually solved by lifting them to distribution theory and using Green's functions.

Let $$ -\nabla^2 \phi(\mathbf{x}) = -4\pi \rho(\mathbf{x})\tag{1} $$

one can show$^1$ that the solution $\phi(\mathbf{x})$ can always be written as$^2$ $$ \phi(\mathbf{x}) = \frac{1}{4\pi} \int_V d^3 x' \rho(\mathbf{x})G(\mathbf{x}, \mathbf{x}') + \frac{1}{4\pi}\int_{\partial V}d\sigma\Big[G(\mathbf{x}, \mathbf{x}')\frac{\partial \phi}{\partial n} - \phi(\mathbf{x}') \frac{\partial}{\partial n}G(\mathbf{x}, \mathbf{x}')\Big] \tag{2} $$ where the Green's function $G(\mathbf{x}, \mathbf{x}')$ solves the associated Green's equation $$ -\nabla^2 G(\mathbf{x}, \mathbf{x}') = -4\pi \delta(\mathbf{x}-\mathbf{x}')\tag{3} $$ Using appropriate boundary conditions (in the simple cases one can demand the functions to vanish at the boundaries and the first order derivative to be proportional to the surface element$^3$) equation $(2)$ can be solved by plugging the solution of $(3)$, which you already have recognised to be the potential of a single point charge.

Details on all of this can be found in the standard textbook of J. D. Jackson on Classical Electrodynamics.

---

$^1$ In order to show why this holds multiply $(1)$ by $G(\mathbf{x}, \mathbf{x}')$ and $(3)$ by $\phi(\mathbf{x})$ and integrate against $\delta(\mathbf{x}-\mathbf{x}')$.

$^2$ There might be some $2\pi$ being forgotten left and right, somewhere.

$^3$ I might be wrong about this.