For a discrete system of $N$ charges, the potential energy associated with their configuration is given by
\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}
where $\varphi$ is the scalar potential due to all charges other than $q_j\,.$
Potential energy stored in the system of continuous charge can be inferred from $(1)$ viz:
$$U= \frac12\int \rho(\mathbf r)\,\varphi(\mathbf r)\,\mathrm d^3\mathbf r \tag 2$$
Now, using $\mathbf\nabla\cdot \mathbf E= \dfrac \rho{\epsilon_0}$ in $(2)$ we get,
$$U= \frac{\epsilon_0}2\int \mathbf \nabla\cdot \mathbf E(\mathbf r)\,\varphi(\mathbf r)\,\mathrm d^3\mathbf r\tag 3$$
Now using the identity $ \mathbf \nabla \cdot \left( \varphi\,\mathbf A \right) = \varphi \left( \mathbf {\nabla }\cdot \mathbf A\right) +\mathbf A \cdot \mathbf {\nabla }\varphi\,^{\dagger}$ in $(3)\,,$ we get
$$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$$
Using divergence theorem; we can rewrite $(4)$ as:
$$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$$
Now, using $\varphi(\mathbf r\to \infty)= 0\,$ and $\mathbf E= -\mathbf \nabla \varphi$ in $(5)\,,$ we get:
$$U= -\frac{\epsilon_0}2 \int\,\mathbf E\cdot (-\mathbf E)\,\mathrm d^3\mathbf r\tag 6$$
And, from $(6)\,$ it is concluded that
$$U= \frac{\epsilon_0}2 \int\,\mathbf E\cdot \mathbf E\,\mathrm d^3\mathbf r\;.$$
How do I interpret it physically?
This is the energy required to assemble the system of charges in that specific configuration.
$^\dagger$ Check this post for its derivation.