Para un sistema discreto de $N$ cargas, la energía potencial asociada a su configuración viene dada por
\begin{align}U&= \frac12\,\sum_{j=1}^N \, q_j\sum_{k\ne j}\,\frac1{4\pi\epsilon_0}\,\cdot \frac{q_k}{r_{jk}}\\ &= \frac12\,\sum_{j=1}^N \, q_j\,\varphi(\mathbf r_j)\tag 1 \end{align}
donde $\varphi$ es el potencial escalar debido a todas las cargas distintas de $q_j\,.$
La energía potencial almacenada en el sistema de carga continua puede deducirse de $(1)$ a saber
$$U= \frac12\int \rho(\mathbf r)\,\varphi(\mathbf r)\,\mathrm d^3\mathbf r \tag 2$$
Ahora, utilizando $\mathbf\nabla\cdot \mathbf E= \dfrac \rho{\epsilon_0}$ en $(2)$ obtenemos,
$$U= \frac{\epsilon_0}2\int \mathbf \nabla\cdot \mathbf E(\mathbf r)\,\varphi(\mathbf r)\,\mathrm d^3\mathbf r\tag 3$$
Utilizando ahora la identidad $ \mathbf \nabla \cdot \left( \varphi\,\mathbf A \right) = \varphi \left( \mathbf {\nabla }\cdot \mathbf A\right) +\mathbf A \cdot \mathbf {\nabla }\varphi\,^{\dagger}$ en $(3)\,,$ obtenemos
$$U= \frac{\epsilon_0}2\left[\int\, \mathbf \nabla \cdot \left( \varphi\,\mathbf E \right)\,\mathrm d^3\mathbf r- \int\, \mathbf E \cdot \mathbf {\nabla }\varphi\,\mathrm d^3\mathbf r\right]\tag 4$$
Utilizando el teorema de la divergencia, podemos reescribir $(4)$ como:
$$U= \frac{\epsilon_0}2\left[\int \varphi\,\mathbf E \,\mathrm d^2\mathbf r- \int\,\mathbf E \cdot \mathbf {\nabla }\varphi \,\mathrm d^3\mathbf r\right]\tag5$$
Ahora, utilizando $\varphi(\mathbf r\to \infty)= 0\,$ y $\mathbf E= -\mathbf \nabla \varphi$ en $(5)\,,$ obtenemos:
$$U= -\frac{\epsilon_0}2 \int\,\mathbf E\cdot (-\mathbf E)\,\mathrm d^3\mathbf r\tag 6$$
Y, de $(6)\,$ se concluye que
$$U= \frac{\epsilon_0}2 \int\,\mathbf E\cdot \mathbf E\,\mathrm d^3\mathbf r\;.$$
¿Cómo lo interpreto físicamente?
Es la energía necesaria para ensamblar el sistema de cargas en esa configuración específica .
$^\dagger$ Comprueba esto Correo electrónico: para su derivación.