Solution (source code)

= Solution

The differential <Gauss law> is $\nabla\cdot\mathbf E=\rho/\epsilon_0$. On a large ball $B_R$, the product rule and $\nabla\phi=-\mathbf E$ give
$$
\rho\phi=\epsilon_0\nabla\cdot(\phi\mathbf E)+\epsilon_0|\mathbf E|^2.
$$
The <divergence theorem> consequently expresses the <electrostatic energy> as
$$
\frac12\int_{B_R}\rho\phi\,dV=\frac{\epsilon_0}{2}\int_{B_R}|\mathbf E|^2\,dV+\frac{\epsilon_0}{2}\int_{\partial B_R}\phi\mathbf E\cdot\mathbf n\,dS.
$$
For a localized charge distribution with the <electric potential> chosen to vanish at infinity, $\phi=O(R^{-1})$ and $\mathbf E=O(R^{-2})$, so the surface term is $O(R^{-1})$ and tends to zero. More generally, it suffices to assume that this boundary integral vanishes. Passing to all space gives
$$
\boxed{U=\frac{\epsilon_0}{2}\int_{\mathbb R^3}|\mathbf E|^2\,dV.}
$$
The physical decay and potential reference are relevant: $\mathbf E\to0$ by itself does not specify the additive constant in $\phi$. Adding a constant $C$ changes the charge-based expression by $CQ/2$ when the total charge is $Q$, while leaving the field unchanged. The zero-at-infinity potential is the usual reference implicit in this all-space <electrostatic energy> formula.