Energy, dielectrics and microscopic electric fields Let us assume all electric fields are microscopic (just to be clear I have no idea what this word means here and is one of the questions addressed below, I am just using it here because it makes the equations work in the general case). We have two electrostatic systems characterized by their electric fields. Let the total field be $\vec E$, the energy of this system will be
$$U=\frac{\epsilon_0}{2}\int E^2\, dV$$
The integrals are always taken over the whole space.
Splitting this into individual fields of the systems we get
$$U=\frac{\epsilon_0}{2}\int (\vec E_1+\vec E_2)^2\, dV$$ $$U=\frac{\epsilon_0}{2}\int E_1^2\, dV+\frac{\epsilon_0}{2}\int E_2^2\,dV+\epsilon_0\int\vec E_1\cdot\vec E_2\, dV\tag{1}$$
This can be written equivalently in terms of charge densities, where $\rho=\rho_1+\rho_2$ and $\phi=\phi_1+\phi_2$ are total charge density and potential respectively. Thus
$$U=\frac12\int\rho\phi\, dV$$ $$U=\frac12\int\rho_1\phi_1\, dV+\frac12\int\rho_2\phi_2\, dV+\frac12\int\rho_1\phi_2\, dV+\frac12\int\rho_2\phi_1\, dV\tag{2}$$
Here the last two integrals can be shown to be equal.
At this point, it seems reasonable to just assume that half of the interaction energy (the sum of the last two integrals of the equation $2$ or equivalently last integral of the equation $1$, as opposed to self-energy of the corresponding system) is in the system $1$ while another half in $2$, doing this way we can easily derive that work done on free charge in a dielectric system is $\frac12\int\vec E\cdot\vec D\, dV$, but this cannot be proved to be true in general, at least not in any way that I am aware of (My question on same).
Now with this, we can proceed toward dielectrics. We will be limiting ourselves to linear dielectric only, the fields are again microscopic. Let $\rho_f$ and $\rho_b$ be the final free charge and bound charge densities, we also have $\phi=\phi_f+\phi_b$, where $\phi$ is the final potential of the whole system, due to conditions stated above. Thus we can use equation $2$ to find the energy needed to bring the free charge into our system. But due to not being able to pinpoint in general what fraction of interaction energy will be used in bringing free charge in this system we can't use that approach. Thus we will use the property of linear dielectrics that when some fraction, $\Lambda$ of free charge is brought into our system, the potential becomes $\Lambda\phi$. This approach is the same that Schwartz followed in his book, Principles of Electrodynamics.
If we bring some charge $dq_f$, while the net potential is $\Lambda\phi$, where $\phi$ is the final potential of the configuration, then the work done is
$$dU=\Lambda\phi dq_f$$
Here $dq_f$ is $Q_fd\Lambda$, where $Q_f$ is total free charge, thus
$$dU=\Lambda\phi Q_fd\Lambda$$
Now we can write $Q_f$ as $\int\rho_f\, dV$, and taking into account that $\phi$ is a function of space we get
$$dU=\Lambda d\Lambda\int\phi\rho_f\, dV \\ ‎\\ U=\int_0^1\Lambda d\Lambda\int\phi\rho_f\, dV$$ $$ U=\frac12 \int\phi\rho_f\, dV \tag{3} $$
This can be shown to be equal to
$$U=\frac12\int\vec E\cdot\vec D\, dV\tag{4}$$
Where $$\vec D=\epsilon_0\vec  E+\vec P$$
Now here we know that this is work done on the free charge while equation $1$ gives the total work done on all the charges, thus it makes sense that if we take the difference of equation $1$ and $4$ we should get work done on charge in dielectric which should be equal to the sum of $2nd$ and $4th$ integral in equation $2$
$$\frac12\int\vec E\cdot(\epsilon_0\vec E-\vec D)\, dV=-\frac12\int\vec E\cdot \vec P\, dV$$
Using $\vec\nabla\cdot\vec P=-\rho_b$ we get
$$\frac12\int\rho_b\phi\, dV=\frac12\int\rho_b(\phi_f+\phi_b)\, dV$$
As expected.
At this point, Griffiths mentions in a footnote on page $199$ of the $4th$ edition of his book that this difference equals work done on bound charge only if $\vec E$ is taken to be a microscopic field.
Questions:

*

*Why does this work only for microscopic fields?

*In deriving equation $3$ we only assumed one system's potential to be dependent on other in a linear fashion, thus any system that follows this property has half of its interaction energy in one system and another half in the other system. Now can we generalize this statement?

 A: Das et al. "Electrostatic energy of a system of charges and dielectrics
" American Journal of Physics 63, 452 (1995); doi: 10.1119/1.18073 do a good job of examining this problem and are more clear in partitioning the energy.

*

*Microscopic vs Macroscopic - Griffiths in chapter 4 goes into some detail that "the microscopic fields inside the dielectric must be fantastically complicated" and then later shows that you don't need to know all the details of the microscopic fields within the dielectric if what you care about is the macroscopic filed - Essentially an average of the microscopic fluctuations for volumes that large compared to the atomic scale, but small compared to the size of the total object.


Notice that it all revolves around the curious fact that the average field over  any sphere (due to the charge inside) is the same as the field at the center of a  uniformly polarized sphere with the same total dipole moment. This means that no  matter how crazy the actual microscopic charge configuration, we can replace it  by a nice smooth distribution of perfect dipoles, if all we want is the macroscopic  (average) field.

This is kind of the important point according to Griffiths since depending on how you assemble your system you are either including the "stretchy springy bits" from the polarizing the electron clouds or you aren't.
The actual foot note is

$^{18}$ The “spring” itself may be electrical in nature, but it is still not included in Eq. 4.55, if E is taken to  be the macroscopic field.
(bolded added)

The Das et al. paper is more explicit and takes the approach
$$U_{Total} = U_{free} + U_{bound} + U_{spring}$$
They go on to find that
$$U_{spring} = \frac{1}{2}\int{P\cdot E dv}$$
They also point out this is different than" the well know formula for the energy of a dielectric in the presence of fixed external charges"
$$U =  -\frac{1}{2}\int{P\cdot E_0 dv}$$
"Where $E_0$ is the field in the absence of the dielectric."
and that in the case for U_{spring} the E is E=$\nabla(V_f+V_b)$
The paper seems to reach the same conclusion about the division of energy"

The energy of the assembly includes the free-charge-free-charge interaction and half of the free-charge-bound charge interaction. Bound charge-bound charge interaction and and half of the free charge - bound charge interaction is stored in (using Griffiths' nomenclature) the "springs", the "twisted and stretched molecules of the dielectric"

So I think they come to the same conclusion as you.
It is interesting that there is nothin about the mechanism of the polarization in this discussion and it might be worth thinking about high field where the P would be nonlinear if there is a nonlinear susceptibility.
