What is the energy density of a linear dielectric and why is it confusing?

[dEdt]

I've examined and re-examined the proof that the energy density of a linear dielectric is $u=\frac{1}{2}\mathbf{D}\cdot\mathbf{E}$, but I still don't understand it. I'll rewrite it here and try to explain what I'm having trouble grasping.

We first imagine that the free charge (as a function of position) changes by an amount $\delta \rho_{\mathrm{free}}(\mathbf{x})$. The change in the potential energy $\delta U$ is then
$$\int \delta\rho_{\mathrm{free}}(\mathbf{x})V(\mathbf{x})\mathrm{d}^3 \mathbf{x}.$$
This can be rewritten as
$$\int\nabla\cdot(\delta \mathbf{D}(\mathbf{x}))V(\mathbf{x})\mathrm{d}^3 \mathbf{x}=\int \nabla\cdot(\delta \mathbf{D} \ V)\mathrm{d}^3 \mathbf{x} + \int \delta \mathbf{D}\cdot\mathbf{E}\mathrm{d}^3 \mathbf{x}.$$
The first term vanished because we're integrating over all of space. Also, if the dielectric is linear we have that $\delta \mathbf{D}\cdot\mathbf{E} =\frac{1}{2}\delta(\mathbf{D}\cdot\mathbf{E})$. Hence
$$\delta U = \frac{1}{2} \int \delta(\mathbf{D}\cdot\mathbf{E})\mathrm{d}^3 \mathbf{x} \rightarrow U=\frac{1}{2}\int \mathbf{D}\cdot\mathbf{E}\mathrm{d}^3 \mathbf{x},$$
proving the desired result.

My problem is with step 1. I'm imagining that time is frozen, and that we come along and sprinkle a little bit of free charge in the dielectric, like dust. Then potential energy changes by the amount $\int \frac{1}{2} \delta (\mathbf{D}\cdot\mathbf{E})\mathrm{d}^3 \mathbf{x}$. So far so good. But if we start time up again, the whole system is going to change as the polarization and bound charges adjust to the introduction of this new free charge.

"But that's not a problem," you say. "During this readjustment, energy is still conserved, so the change in $\delta U$ is still given by that formula."

Fine, but here's the issue: after this readjustment, the change in $\mathbf{D}\cdot\mathbf{E}$ is not the $\delta(\mathbf{D}\cdot\mathbf{E})$ that appears in the formula above. Why? Because during the readjustment, $\mathbf{D}\cdot\mathbf{E}$ changes as well. That means that as we incrementally add free charge, we cannot conclude that the incremental change in energy is $\int \frac{1}{2} \Delta (\mathbf{D}\cdot\mathbf{E})\mathrm{d}^3 \mathbf{x}$ (where $\Delta (\mathbf{D}\cdot\mathbf{E})$ is the actual change in D dot E), because we haven't taken into account that $\mathbf{D}\cdot\mathbf{E}$ changes from readjustment.

Also, if we want to calculate the change in energy when free charge is added, why don't we take into account the change in energy from the fact that the potential changes ie why don't we include a term like $\delta U = \int \rho \delta V\mathrm{d}^3 \mathbf{x}$?

 

[dEdt]

Anyone?