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?