A rigorous proof that $\nabla \cdot E = \frac{\rho}{\epsilon_0}$

Question

Suppose that $\rho: \mathbb R^3 \to \mathbb R$ is a function that tells us the electric charge density at each point in space. According to Coulomb's law, the electric field at a point $x \in \mathbb R^3$ is $$ E(x) = \frac{1}{4 \pi \epsilon_0} \int_{\mathbb R^3} \rho(y) \frac{x - y}{\|x - y \|^3} \, dy. $$ Question: How do you prove rigorously that, for the function $E$ defined above, we have $$ \tag{1} \nabla \cdot E = \frac{\rho}{\epsilon_0} . $$

Comments: It's not easy to find a rigorous proof of this because physics textbooks tend to give non-rigorous arguments, whereas math books tend not to discuss this at all because it's about physics.

Here's a typical physics derivation of this fact: \begin{align} \nabla \cdot E(x) &= \frac{1}{4 \pi \epsilon_0} \int_{\mathbb R^3} \rho(y) \, \nabla \cdot \left(\frac{x - y}{\|x - y \|^3} \right) \, dy \\ &= \frac{1}{4 \pi \epsilon_0} \int_{\mathbb R^3} \rho(y) \, 4 \pi \,\delta(x - y) \, dy \\ &= \frac{\rho(x)}{\epsilon_0}. \end{align}

This non-rigorous argument is based on the assertion that the divergence of the function $x \mapsto \frac{x}{\|x\|^3}$ is $4 \pi \delta$, where $\delta$ is the delta function on $\mathbb R^3$. That claim is typically justified by non-rigorous arguments.

Bonus question: Where would I look to find a rigorous proof of this or similar physics equations? Is there a book on electromagnetism written for mathematicians? A similar question was asked on the physics stackexchange, but it hasn't received a good answer. Since the question is about providing a rigorous proof I think it's a better fit for math.stackexchange.

---

Edit: I think a key insight I was missing is that proving (1) is equivalent to showing that $$ \tag{2} \text{if} \quad \varphi(x) = -\frac{1}{4 \pi \epsilon_0} \int_{\mathbb R^3} \rho(y) \| x - y \|^{-1} \, dy \quad \text{then}\quad \nabla^2 \varphi = \frac{\rho}{\epsilon_0}. $$ The function $-\varphi$ is called the electric potential in electrostatics. This fact (2) is fundamental to the study of the Laplace equation and Poisson's equation, and it is of course proved rigorously in PDEs textbooks such as Evans. So, it is not that there is a gap in the mathematical literature -- rather, if you want to find a rigorous proof of (1), you just need to recognize that it's equivalent to (2).

Why does (2) imply (1)? First notice that if $y \in \mathbb R^3$ and $h(x) = -\|x - y \|^{-1}$ then $\nabla h(x) = \frac{x - y}{\|x - y \|^3}$ for all $x \in \mathbb R^3, x \neq y$. It follows that \begin{align} \tag{3} \nabla \varphi(x) &= \frac{1}{4 \pi \epsilon_0} \int_{\mathbb R^3} \rho(y) \nabla \left( -\| x - y \|^{-1} \right) \, dy \\ &= \frac{1}{4 \pi \epsilon_0} \int_{\mathbb R^3} \rho(y) \frac{x - y}{\|x - y \|^3} \, dy \\ &= E(x). \end{align} So, $$ \nabla^2 \varphi(x) = \nabla \cdot \nabla \varphi(x) = \nabla \cdot E(x). $$ Thus, if we can show that $\nabla^2 \varphi = \frac{\rho}{\epsilon_0}$, then we will have shown that $\nabla \cdot E = \frac{\rho}{\epsilon_0}$.

By the way, it is still not clear to me how to prove (3) rigorously. I think we need to use the dominated convergence theorem or one of the theorems from measure theory which imply that, under certain conditions, the limit of the integral is the integral of the limit.

Answer 1

Use of the Dirac Delta is perfectly rigorous in the context of distributions. But we don't require the machinery of generalized functions to proceed.

In the ensuing development, we assume that $\rho\in C_C^\infty$. This is not a necessary condition, but facilitates a rigorous and straightforward proof.

First, note that $\nabla_x^2 \frac{1}{\|x-y\|}=0$ for $x\ne y$ and $\nabla_x \frac{1}{\|x-y\|}=-\nabla_y \frac{1}{\|x-y\|}=-\frac{x-y}{\|x-y\|^3}$. Now, we can write

$$\begin{align} E(x)&=\frac1{4\pi\varepsilon_0}\int_{\mathbb{R}^3}\rho(y)\frac{x-y}{\|x-y\|^3}\,dy\\\\ &=\frac1{4\pi\varepsilon_0}\int_{\mathbb{R}^3}\rho(y)\nabla_y \frac{1}{\|x-y\|}\,dy\\\\ &=-\frac1{4\pi\varepsilon_0}\int_{\mathbb{R}^3}\frac{1}{\|x-y\|}\nabla_y\rho(y) \,dy\\\\ \end{align}$$

where we integrated by parts and exploited the compact support of $\rho$. Now, taking the divergence of $E(x)$, and denoting the spherical domain centered at $x$ with radius $\varepsilon$ by $B_\varepsilon(x)$ (i.e., $\|x-y\|=\varepsilon$ on $B_\varepsilon(x)$), we find that

$$\begin{align}\require{cancel} \nabla\cdot E(x)&=-\frac1{4\pi\varepsilon_0}\int_{\mathbb{R}^3}\nabla_x\left(\frac{1}{\|x-y\|}\right)\cdot \nabla_y\rho(y) \,dy\\\\ &=\lim_{\varepsilon\to 0^+}\frac1{4\pi\varepsilon_0}\int_{\mathbb{R}^3\setminus B_\varepsilon(x)}\nabla_y\left(\frac{1}{\|x-y\|}\right)\cdot \nabla_y\rho(y) \,dy\\\\ &=\lim_{\varepsilon\to 0^+}\frac1{4\pi\varepsilon_0}\int_{\mathbb{R}^3\setminus B_\varepsilon(x)}\nabla_y\cdot \left(\rho(y)\nabla_y\frac{1}{\|x-y\|}\right) \,dy\\\\ &-\lim_{\varepsilon\to 0^+}\frac1{4\pi\varepsilon_0}\int_{\mathbb{R}^3\setminus B_\varepsilon(x)}\rho(y)\cancelto{0}{\nabla_y^2\frac{1}{\|x-y\|}} \,dy\\\\ &=\lim_{\varepsilon\to 0^+}\frac1{4\pi\varepsilon_0}\oint_{\partial B_\varepsilon(x)}\rho(y) \frac{x-y}{\|x-y\|}\cdot \left(\nabla_y\frac{1}{\|x-y\|}\right) \,dy\\\\ &=\lim_{\varepsilon\to 0^+}\frac1{4\pi\varepsilon_0}\oint_{\partial B_\varepsilon(x)}\rho\left(x-\varepsilon \frac{x-y}{\|x-y\|}\right)\,d\Omega\\\\ &=\frac{\rho(x)}{\varepsilon_0} \end{align}$$

as was to be shown!