What is the physical meaning of the energy density of an electrostatic field?

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}{2mu} |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_Qcdot 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 Ecdot frac{dvec 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.

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}^nsum_{substack{j=1 j neq i }}^{n}frac{1}{4piepsilon_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}intlimits_V rho(vec{r})phi(vec{r})textrm{d}V$$Now, due to Gauss's Law for electricity, $vec{nabla}cdotvec{E}=frac{rho}{epsilon_0}$ we have $$U=frac{1}{2}intlimits_V epsilon_0left(vec{nabla}cdotvec{E}(vec{r})right)phi(vec{r})textrm{d}V$$Recalling the vector identity $vec{nabla}cdotleft(fvec{F}right)=f(vec{nabla}cdotvec{F})+vec{F}cdot(vec{nabla}f)$ we have $$U=frac{1}{2}epsilon_0left[intlimits_V vec{nabla}cdotleft(vec{E}(vec{r})phi(vec{r})right)textrm{d}V-intlimits_V vec{E}(vec{r})cdotvec{nabla}phi(vec{r})textrm{d}Vright]$$By replacing the first integral over the volume $V$ for one over its boundary $partial V$ through the divergence theorem which states $int_Vvec{nabla}cdotvec{F}textrm{d}V=oint_{partial V}vec{F}cdottextrm{d}vec{S}$ and by remembering the definition of potential $-vec{nabla}phi=vec{E}$ we get $$U=frac{1}{2}epsilon_0left[ointlimits_{partial V} vec{E}(vec{r})phi(vec{r})textrm{d}vec{S}+intlimits_V vec{E}(vec{r})cdotvec{E}(vec{r})textrm{d}Vright]$$ 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_0intlimits_V |vec{E}(vec{r})|^2textrm{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!