What is the physical meaning of the energy density of an electrostatic field?
Actually in electrostatics energy density of E-field is not a physical observable. As you say, only when charges move will there be any work done. Since the two ways of calculating total energy end the same, you cannot distinguish whether energy is stored on the charges or in the field. Even E-field itself is more of an abstract mathematical entity, without which everything can be calculated in terms of Coulomb law.
The physical reality of E and B fields (and the energy density associated) becomes apparent only in non-static cases. For example, in electromagnetic radiation, fields can propagate in free space without being associated with charges and currents, and the radiation may do work on non-charges (for example, light pressure). Because from Maxwell equations we can derive a general formula of energy density
$$\rho = \frac{\epsilon_0}{2} |\vec E|^2 + \frac{1}{2\mu_0} |\vec B|^2$$
which coincides with the electrostatic case, we deduce that even in electrostatics energy is indeed stored in the fields.
When one has a distribution of charges $q_1,\dots,q_n$ at points $\vec{r}_1,\dots,\vec{r}_n$, the energy of the system is given by the sum of the energy of each particle due to its interaction with the others divided by two since each interaction is counted twice i.e. $$U=\sum_{i=1}^n\sum_{\substack{j=1 \\ j \neq i }}^{n}\frac{1}{4\pi\epsilon_0}\frac{q_iq_j}{\|\vec{r}_i-\vec{r}_j\|}=\frac{1}{2}\sum_{i=1}^nq_i \phi_i(\vec{r}_i)$$where $\phi_i(\vec{r}_i)$ is the potential at $\vec{r}_i$ due to all charges except $q_i$. If we go over to a continuous charge distribution with charge density $\rho(\vec{r})$, the summation is replaced by an integration over "infinitesimal chunks of charge" $dq=\rho(\vec{r}) \textrm{d}V$. Then the energy of the system is $$U=\frac{1}{2}\int\limits_V \rho(\vec{r})\phi(\vec{r})\textrm{d}V$$Now, due to Gauss's Law for electricity, $\vec{\nabla}\cdot\vec{E}=\frac{\rho}{\epsilon_0}$ we have $$U=\frac{1}{2}\int\limits_V \epsilon_0\left(\vec{\nabla}\cdot\vec{E}(\vec{r})\right)\phi(\vec{r})\textrm{d}V$$Recalling the vector identity $\vec{\nabla}\cdot\left(f\vec{F}\right)=f(\vec{\nabla}\cdot\vec{F})+\vec{F}\cdot(\vec{\nabla}f)$ we have $$U=\frac{1}{2}\epsilon_0\left[\int\limits_V \vec{\nabla}\cdot\left(\vec{E}(\vec{r})\phi(\vec{r})\right)\textrm{d}V-\int\limits_V \vec{E}(\vec{r})\cdot\vec{\nabla}\phi(\vec{r})\textrm{d}V\right]$$By replacing the first integral over the volume $V$ for one over its boundary $\partial V$ through the divergence theorem which states $\int_V\vec{\nabla}\cdot\vec{F}\textrm{d}V=\oint_{\partial V}\vec{F}\cdot\textrm{d}\vec{S}$ and by remembering the definition of potential $-\vec{\nabla}\phi=\vec{E}$ we get $$U=\frac{1}{2}\epsilon_0\left[\oint\limits_{\partial V} \vec{E}(\vec{r})\phi(\vec{r})\textrm{d}\vec{S}+\int\limits_V \vec{E}(\vec{r})\cdot\vec{E}(\vec{r})\textrm{d}V\right]$$ Now we may choose the volume of integration $V$ to be all space. Then $\partial V$ would be infinitely far away from all charges and by convention the potential $\phi$ would die out, making the first integral to dissapear. Then we would be left with the following expression for the energy of the system $$U=\frac{1}{2}\epsilon_0\int\limits_V \|\vec{E}(\vec{r})\|^2\textrm{d}V$$ From this expression follows that the energy density $\Upsilon$ at $\vec{r}$ is given by $$\Upsilon(\vec{r})=\frac{1}{2}\epsilon_0\|\vec{E}(\vec{r})\|^2$$Now we may ask "where" is this energy. Note we derived it as the energy density of the total energy of the charges in our universe. Non the less, this energy density seems to be distributed even where there may be no charges. So we ask the question, does the energy belong to the charge configuration or to the electric field? From the point of view of electrostatics, both are equivalent, and in our last equation, we are inclined to think of the electric field as being the carrier of energy. In electromagnetism, this is a necessity, since the electric field may exist and propagate quite independently of its source charges in the form of electromagnetic waves, for example light, which evidently carries energy, since most of the energy in planet earth is carried from the sun in this form. Any other question, please ask! I'll try to give an answer if I know one!
Yes, $\epsilon \vec E \cdot \vec E$ is the electrostatic part of the energy density carried by the field. The energy density of the electromagnetic field also includes the magnetic term: $$ \rho_{E,B} = \frac{\epsilon}{2} |\vec E|^2 + \frac{1}{2\mu} |\vec B|^2 $$ and this formula is valid even for arbitrary time-dependent, variable electromagnetic fields. When you mentioned the energy density $$ \frac 12 \int \rho_{\rm charge} \Phi \,\,dV, $$ one should note that one must be careful to avoid double-counting. When we assume that the energy is carried by the electromagnetic field, we should no longer add the $\rho_Q\cdot \Phi$ term separately because we could be double-counting. However, in some respects, they have to be separated and both of them have to be added.
At any rate, $\epsilon|E|^2$ is a term in the formula for the total energy, anyway. It's important to know because only the total energy, with all the terms that should be there, is conserved.
One may interpret the energy $\int dV\,\,\epsilon|E|^2/2$ as work, in the same way as for the interaction energy of the charges you mention. It's the work needed to change the electrostatic field from the situation $\vec E=0$ to the given configuration of $\vec E$. The energy may be given as an integral of the work, $$ E_{\rm energy} = \int dV \int dt\, \vec E\cdot \frac{d\vec D}{dt},\qquad \vec D \equiv \epsilon \vec E $$ Note that there is no $1/2$ in the formula above; it comes from the integration. So the larger the field is at a given point, the harder it is to increase its value there.