Numerically Stable Light Absorption Density Calculation I am reading through the documentation provided by Lumerical concerning light absorption per unit volume : https://kb.lumerical.com/layout_analysis_pabs_simple.html
Lumerical says that: 

It can be shown that the above formula is equivalent to
$P a b s = 0.5 \operatorname { real } \left( i \omega \vec { E } \cdot
\vec { D } ^ { * } \right)$

Where can I find a good derivation of this formula ? 
Unfortunately, Lumerical does not provide any reference. 
 A: I made a wikipedia article 10 years ago with a compendium of formulas related to light absorption:
https://en.wikipedia.org/wiki/Mathematical_descriptions_of_opacity#Complex_electric_permittivity
There are textbook references throughout. You may also need definitional formulas like $D=\epsilon E$ and $D = \epsilon_0 E + P$, and things you learn from experience like replacing $d/dt$ by $-i\omega$ for a parameter proportional to $e^{-i\omega t}$, and of course there's Maxwell's equations.
A: Let us consider the conservation of energy in a small volume element in space. The work done per unit volume by an electromagnetic field is $ \vec{\text{J}}\cdot \vec{\text{E}} $ which is the energy dissipation in a unit volume. This is associated with net decrease in the enrgy density and power flow out of the volume.
\begin{equation}
\vec{J} \cdot \vec{E} = \vec{E}\cdot({\nabla}\times \vec{H})-\vec{E}\cdot\frac{\partial \vec{D}}{\partial t}
\end{equation}
\begin{equation}
\nabla\cdot(\vec{E}\times\vec{H})=\vec{H}\cdot(\nabla\times\vec{E})-\vec{E}\cdot(\nabla\times\vec{H})
\end{equation}
Implies,
\begin{equation}
\vec{J} \cdot \vec{E}=-\nabla\cdot(\vec{E}\times\vec{H})-\vec{H}\cdot\frac{\partial\vec{B}}{\partial t}-\vec{E}\cdot\frac{\partial\vec{D}}{\partial t}
\end{equation}
If the material medium is assumed to be linear in its electromagnetic properties, that is, $ \epsilon$ and $\mu$ are independent of the field strength, then,
\begin{equation}
\frac{\partial U}{\partial t}+\nabla\cdot\vec{S}=-\vec{J}\cdot\vec{E}
\end{equation}
where,
\begin{equation}
U = \frac{1}{2}(\vec{E}\cdot\vec{D}+\vec{B}\cdot\vec{H})
\end{equation}
is the energy density of the EM field with units $ \text{Jm}^{-3} $ and
\begin{equation}
\vec{S} = \vec{E}\times\vec{H}
\end{equation}
is the energy flux called the Poynting vector with units $\text{Js}^{-1}\text{m}^{-2}$ or $ \text{Wm}^{-2}$
If there is no current flow, then, 
\begin{equation}
\frac{\partial U}{\partial t}+\nabla\cdot\vec{S}=0
\end{equation}
Since $\vec{B}=\mu_0(\vec{H}+\vec{M})$ and since $\vec{D}=\epsilon_0\vec{E}+\vec{P},$
\begin{equation}
\vec{J} \cdot \vec{E}=-\nabla\cdot(\vec{E}\times\vec{H})-\mu_0\vec{H}\cdot\frac{\partial\vec{H}}{\partial t}-\mu_0\vec{H}\cdot\frac{\partial\vec{M}}{\partial t}-\epsilon_0\vec{E}\cdot\frac{\partial\vec{E}}{\partial t}-\vec{E}\cdot\frac{\partial\vec{P}}{\partial t}
\end{equation}
Now, the dipolar dissipation is mainly focused on the last term, $\vec{E}\cdot\frac{\partial \vec{P}}{\partial t}$.
For a complex quantity, 
\begin{equation}
A=|A|e^{i(\omega t+\alpha)}
\end{equation}
Let the real part be
\begin{equation}
a(t) = \mathfrak{R}(A)=\mathfrak{R}\Big(|A|e^{i(\omega t+\alpha)}\Big)=|A|\cos(\omega t+\alpha)
\end{equation}
\begin{equation}
\dot{a}(t)=-|A|\omega\sin(\omega t+\alpha)
\end{equation}
Using the time averaging of sinusoidal products,
\begin{equation}
\langle a(t)b(t)\rangle=\int|A|\cos(\omega t+\alpha)|B|\cos(\omega t+\beta)=\frac{1}{2}|A||B|\cos(\alpha-\beta)=\frac{1}{2}\mathfrak{R}(A^*B)
\end{equation}
This gives the time-averaged Poynting vector as
\begin{equation}
\vec{S}=\frac{1}{2}\mathfrak{R}\big(\vec{E}\times\vec{H}^*\big)
\end{equation}
Thus, the time averaged dipolar dissipation is
\begin{equation}
P_D = \frac{1}{2}\Bigg(\vec{E}^*\cdot\frac{\partial\vec{P}}{\partial t}\Bigg)
\end{equation}
\begin{equation}
P_D = \mathfrak{R}\Bigg(\frac{1}{2}\Bigg(\frac{\vec{D}^*}{\epsilon}\cdot\frac{\epsilon_0\chi\partial\vec{E}}{\partial t}\Bigg)\Bigg)
\end{equation}
\begin{equation}
P_D = \mathfrak{R}\Bigg(\frac{1}{2}\Bigg(\frac{\vec{D}^*}{\epsilon_0(\chi+1)}\cdot\frac{\epsilon_0\chi\partial\vec{E}}{\partial t}\Bigg)\Bigg)
\end{equation}
For highly susceptible materials,
\begin{equation}
\frac{\chi}{\chi+1}\approx 1
\end{equation}
Thus,
\begin{equation}
P_D = \mathfrak{R}\Bigg(\frac{1}{2}\Bigg(\vec{D}^*\cdot\frac{\partial\vec{E}}{\partial t}\Bigg)\Bigg)
\end{equation}
Since the partial derivative with respect to time can be written in the phasor or the frequency domain as $ s = \pm i\omega $ depending on the phase i.e. $e^{\pm i\omega t}$,
\begin{equation}
P_D = \mathfrak{R}\Bigg(\frac{1}{2}\left|(\vec{D}^*\cdot i\omega\vec{E})\right|\Bigg)
\end{equation}
\begin{equation}
P_D = \mathfrak{R}\Bigg(\frac{1}{2}\left|(i\omega \vec{E}\cdot\vec{D}^*)\right|\Bigg)
\end{equation}
