Energy density with complex permittivity? Question
What is the correct form of the energy density when we have a complex permittivity (such as in a dielectric)? 
Additional information
A complex permittivity means that we have: 
$$\vec D_r +i\vec D_i =(\epsilon_r+i\epsilon_i)  (\vec E_r+i\vec E_i) $$
Giving: 
$$\vec D_r=\epsilon_r \vec E_r-\epsilon_i \vec E_i$$
To me this would indicate an energy density of:
$$u=\frac{1}{2}\vec D_r \cdot \vec E_r=\frac{1}{2}(\epsilon_r \vec E_r\cdot \vec E_r-\epsilon_i \vec E_i\cdot \vec E_r) $$
However, this form is inconsistent with the energy lost in a lossy dielectric whist the form:
$$u=\frac{1}{2} \epsilon_r \vec E_r \cdot \vec E_r$$
Is consistent, showing that it is most likely the last one which is the correct form of the energy density but this makes no sense. So what is the correct form of energy density with complex permittivity and why?
 A: Jackson's E&M book covers this in Chapter 6, Section 8 (3rd Edition, or the blue covered book).  The following is a summary of that section.
Answer (not much explanation)
In a lossy medium, where $\Im \left[ \varepsilon \right] \neq 0$ and/or $\Im \left[ \mu \right] \neq 0$, we can still use linear approximations but there are some modification.  In Fourier space, we can define:
$$
\mathbf{D}\left( \mathbf{x}, \omega \right) = \varepsilon \left( \omega \right) \ \mathbf{E}\left( \mathbf{x}, \omega \right) \\
\mathbf{B}\left( \mathbf{x}, \omega \right) = \mu \left( \omega \right) \ \mathbf{H}\left( \mathbf{x}, \omega \right)
$$
Then Poynting's theorem goes to (after a bunch of algebra):
$$
\frac{ \partial w_{eff} }{ \partial t } + \nabla \cdot \mathbf{S} = -\mathbf{j} \cdot \mathbf{E} - 2 \omega_{o} \Im \left[ \varepsilon \left( \omega_{o} \right) \right] \langle \mathbf{E}\left( \mathbf{x}, t \right) \cdot \mathbf{E}\left( \mathbf{x}, t \right) \rangle_{\omega_{o}}  \\
- 2 \omega_{o} \Im \left[ \mu \left( \omega_{o} \right) \right] \langle \mathbf{H}\left( \mathbf{x}, t \right) \cdot \mathbf{H}\left( \mathbf{x}, t \right) \rangle_{\omega_{o}}
$$
where $\langle \rangle_{\omega_{o}}$ is an ensemble average over the narrow frequency $\omega_{o}$ that both $\mathbf{E}$ and $\mathbf{H}$ were assumed to oscillate in order to get the above expression.  That expression also assumes that both $\mathbf{E}$ and $\mathbf{H}$ are slowly varying with respect to both $1/\omega_{o}$ and the frequency range over which $\varepsilon \left( \omega \right)$ is assumed to vary.  Note that the effective energy density, $w_{eff}$, is given by:
$$
w_{eff} = \Re \left[ \frac{ d \left( \omega \varepsilon \right) }{ d\omega } \left( \omega_{o} \right) \right] \langle \mathbf{E}\left( \mathbf{x}, t \right) \cdot \mathbf{E}\left( \mathbf{x}, t \right) \rangle_{\omega_{o}} + \Re \left[ \frac{ d \left( \omega \mu \right) }{ d\omega } \left( \omega_{o} \right) \right] \langle \mathbf{H}\left( \mathbf{x}, t \right) \cdot \mathbf{H}\left( \mathbf{x}, t \right) \rangle_{\omega_{o}}
$$
In the above version of Poynting's theorem, the terms are physically described as follows:


*

*$\partial_{t} w_{eff} \neq 0$ $\rightarrow$ rate of change of electromagnetic energy density

*$\nabla \cdot \mathbf{S} \neq 0$ $\rightarrow$ rate of energy density flowing into/out of region of interest (more generally the surface used to define the region)

*$-\left( \mathbf{j} \cdot \mathbf{E} \right)$ $\rightarrow$ explicit Ohmic losses

*(2nd and 3rd terms on right-hand side) $\rightarrow$ absorptive dissipation

*$\Im \left[ \varepsilon \right] \neq 0$ and/or $\Im \left[ \mu \right] \neq 0$ $\rightarrow$ conversion of electromagnetic energy into average particle kinetic energy (i.e., heat) or different forms of radiation


Notes on Poynting's Theorem
The full complex version describes the following:


*

*$\Re \left[ \right]$ $\rightarrow$ conservation of energy for the ensemble average of the quantities

*$\Im \left[ \right]$ $\rightarrow$ defines the "...reactive or stored energy and its alternating flow..."
Details to get Equations Above
The presence of any dispersive fields causes $\mathbf{E} \cdot \partial_{t} \mathbf{D} \neq \partial_{t} \left( \mathbf{E} \cdot \mathbf{D} \right)/2$, which is another way of saying the two share a "...temporally nonlocal connection..."  Instead we need to write this in terms of an integral over frequency as:
$$
\mathbf{E} \cdot \partial_{t} \mathbf{D} \rightarrow \frac{ 1 }{ 2 } \int \ d\omega \ \int \ d\omega' \ \mathbf{E}^{*}\left( \mathbf{x}, \omega' \right) \left[ -i \ \omega \ \varepsilon \left( \omega \right) + i \ \omega' \ \varepsilon^{*} \left( \omega' \right) \right] \cdot \mathbf{E}\left( \mathbf{x}, \omega \right) \ e^{-i \left( \omega - \omega' \right) t}
$$
where $Q^{*}$ represents the complex conjugate of $Q$.  If we use the assumption that $\mathbf{E}$ is dominated by frequency components in some narrow (compared to the frequency range over which $\varepsilon \left( \omega \right)$ varies) frequency range near $\omega_{o}$, then we can use a Taylor expansion for the $i \ \omega' \ \varepsilon^{*}\left( \omega' \right)$ factor about $\omega = \omega'$, then the terms in the brackets are given by:
$$
\left[ ... \right] \rightarrow 2 \ \omega \ \Im \left[ \varepsilon \left( \omega \right) \right] - i \ \left( \omega - \omega' \right) \frac{ d }{ d\omega } \left[ \omega \ \varepsilon^{*}\left( \omega \right) \right] + ...
$$
We can then rewrite $\mathbf{E} \cdot \partial_{t} \mathbf{D}$ as the following:
$$
\int \ d\omega \ \int \ d\omega' \ \mathbf{E}^{*}\left( \mathbf{x}, \omega' \right) \cdot \mathbf{E}\left( \mathbf{x}, \omega \right) \ \omega \ \Im \left[ \varepsilon \left( \omega \right) \right] \ e^{-i \left( \omega - \omega' \right) t} \\
\partial_{t} \frac{ 1 }{ 2 } \int \ d\omega \ \int \ d\omega' \ \mathbf{E}^{*}\left( \mathbf{x}, \omega' \right) \cdot \mathbf{E}\left( \mathbf{x}, \omega \right) \ \frac{ d }{ d\omega } \left[ \omega \ \varepsilon^{*}\left( \omega \right) \right] \ e^{-i \left( \omega - \omega' \right) t}
$$
The corresponding expression for $\mathbf{H} \cdot \partial_{t} \mathbf{B}$ is the same as above when one substitutes $\mathbf{E} \rightarrow \mathbf{H}$ and $\varepsilon \rightarrow \mu$.
A: After some thought I think I have come up with a reason why the energy density is indeed:
$$u=\frac{1}{2} \epsilon_r \vec E_r \cdot \vec E_r$$
My reasoning is that the imaginary parts of $E$ and $\epsilon$  solely come about due the propagation of the wave through a medium. The energy density of the electric field, however should not depend on how it is propagating (at a given time) and thus when finding the energy density we should ignore these imaginary parts. 
A: Actually, $w_0=\frac{1}{2} \bf{D} \cdot E$ can be right only in a linear dielectric without dissipation( permittivity is a real number).
In an arbitrary dielectric,
$$\nabla \centerdot \mathbf{D}=\rho_0 $$
$$\nabla \times \mathbf{B}=-\frac{\partial\mathbf{E}}{\partial t}$$
$$\nabla \centerdot \mathbf{B}=0 $$
$$\nabla \times \mathbf{H}=\bf{j}_0 + \frac{\partial\mathbf{D}}{\partial t}$$
The power density of EM field to the free charge and conducting current is
$$p_0=\bf{j}_0 \cdot \bf{E} =(\nabla \times \mathbf{H}-\frac{\partial\mathbf{D}}{\partial t}) \cdot \bf{E}=-\bf{E} \cdot \frac{\partial\mathbf{D}}{\partial t}-\bf{H} \cdot \frac{\partial\mathbf{B}}{\partial t} -\nabla \centerdot(\bf{E \times H})$$
If the dielectric is linear and without dissipation or dispersion,then we have
$$\bf{E} \cdot \frac{\partial\mathbf{D}}{\partial t}= \epsilon_{ij}E_i \frac{\partial\mathbf{E_j}}{\partial t}= \epsilon_{ji}E_i \frac{\partial\mathbf{E_j}}{\partial t}=\bf{D} \cdot \frac{\partial\mathbf{E}}{\partial t}$$
Similarly, we have 
$$\bf{H} \cdot \frac{\partial\mathbf{B}}{\partial t}=\bf{B} \cdot \frac{\partial\mathbf{H}}{\partial t}$$
So, we have
$$p_0=-\frac{\partial w_0}{\partial t}- \nabla \cdot \bf{S}_0$$
$$w_0=\frac{1}{2}(\bf{D} \cdot E +B \cdot H)$$
$$\bf{s}_0=E \times H$$
Then, we can assume $w_0$ as energy density of EM field. We must notice that this energy density include both true energy density of EM field and energy of polarization and magnetization. 
A: It depends on what you mean by "energy density". If you want to ascribe some energy to EM field in material medium in a general way, not requiring linear and lossless behaviour of the medium, the simplest way is to use the same Poynting theorem and formulae as in vacuum. Then density of EM energy is
$$
\frac{1}{2}\epsilon_0E^2 + \frac{1}{2\mu_0} B^2
$$
and permittivity of the medium does not enter EM energy density at all. This only gives energy ascribed to EM field and completely neglects energy of the medium.
If you want to ascribe energy to EM field in the medium in a way similar to how it is often done in textbooks for linear lossless media - that is, to include medium's energy into the expression so that there is no other energy left - then this is not easily possible in general case.
This is because in general, charge and current density, or alternatively polarizations of the medium, are not functions of the EM field. Only when they are, for example when medium is linear and lossless, one can express them via EM field vectors and introduce energy density that makes the work-energy theorem into local conservation law with energy density
$$
\frac{1}{2}\epsilon E^2 + \frac{1}{2\mu} B^2.
$$
