Superposition of electromagnetic waves and energy localization Sorry for the long question. I didn't know how to make it shorter.
I'm trying to understand how energy is spatially localized in an electromagnetic wave. My premises are:


*

*Electromagnetic energy is spatially localized. E and H fields store energy. An energy density can be defined as a function of E and H; and a Poynting vector can also be defined which characterizes the flow of power across an infinitesimal surface.

*Fields obey the principle of superposition. When two waves travelling in different directions (such as an incident and a reflected wave) encounter each other, they go through one another "like they were ghosts", meaning one wave doesn't "know" about the existence of the other. The total field is just the sum of the individual fields.

*Energy doesn't obey superposition. When two waves spatially coincide, the energy density is not the sum of the energies that each wave would carry by itself. There are interaction (or "crossed") terms. Those interaction terms sometimes cancel, but they exist in general.

*I consider linear media with dielectric losses. In the frequency domain, this means $\epsilon$ has a nonzero imaginary part, and then by the Kramers-Kronig relations $\epsilon$ is a function of frequency.


Consider a rectangular electromagnetic pulse of a given lenght that travels from left to right and incides on a perfect-conductor vertical wall. A total energy can be ascribed to that pulse. Consider a time instant when part of the pulse has already incided on the wall and has been reflected, and part is still travelling towards the wall.
In a region to the left of the wall there is a superposition of incident and reflected waves. In order to compute the energy density in that region, one has to consider the total fields. As said above, the energy density in general is not the sum of terms corresponding to each wave (there are interaction terms). Due to this, the energy "moves differently" (and in a more complicated manner) than the fields do.
I'd like to analyze these issues by myself. Ideally I'd like to carry out some numerical simulation in Matlab to "see" how the energy is distributed in space as time passes; and how the total energy, including that dissipated in the medium, is conserved. I did a Matlab program for that, but got stuck because $\epsilon$ is a function of frequency, and I want to analyze in time domain.
Now, my questions:


*

*Are all my premises correct? Am I missing something?

*Do you know any good reference (book, paper, web) that deals with these issues? Note that analyses of normal incidence in lossless media may not be suitable, because in that particular case energy does satisfy a kind of superposition. So that case is not general enough for the effects I'm interested in.

*Is there a way to analyze these effects (including dielectric losses) in time domain? How does one deal with the time-domain equivalent of a complex, dispersive $\epsilon$? Any reference or pointer in this direction?

 A: I imagine it is straightfoward to describe a lossy medium in the time domain, using a formulation in which the polarisation is determined not only by the field by an additional frictional term.  The problem is that you will then be faced with solving Maxwell's equations in their full form - you can't use the wave solutions, which obviously assume a single frequency.  This is of course why the frequency domain was invented - you can represent your wave-packet by a Fourier decomposition, but this introduces additional complexity and approximations.  
You are in that awkward intermediate regime - linear but dispersive.
A: With regards to your last question (analysis in the time domain), @akrasia is correct that you will basically need to go back to Maxwell's equations.  Take a look at Ampere's Law: $\nabla\times {\bf H}={\bf J}_{\rm f}+\frac{\partial {\bf D}}{\partial t}$.  In it, there are the terms which capture the medium properties: ${\bf J}_{\rm f}=\sigma {\bf E}$ and ${\bf D}={\bf E}+{\bf P}={\bf E}+\chi {\bf E}$.
The $\sigma {\bf E}$ and $\chi {\bf E}$ terms describe the linear response of the medium.  Written as such, it is easy to work in the frequency domain, with frequency-dependent quantities.  That is, ${\bf J}_{\rm f}(\omega)=\sigma (\omega) {\bf E}(\omega)$ and ${\bf P}(\omega)=\chi (\omega) {\bf E}(\omega)$.
If instead you prefer to work in the time domain, it's no problem. Just take your frequency-dependent $\chi(\omega)$ and inverse-Fourier-transform ($\mathscr{F}^{-1}$) it into the time domain: $\chi(t)=\mathscr{F}^{-1}[\chi(\omega)](t)$.  Then the linear response is a convolution ($\otimes$), e.g. ${\bf P}(t)=\chi (t)\otimes {\bf E}(t)$.  Alternatively, it may be easier to Fourier-transform your electric field pulse into frequency domain first and go from there (i.e. ${\bf P}(t)=\mathscr{F}^{-1}[\ \chi (\omega)\ \mathscr{F}[{\bf E}(t)](\omega)\, ]$).
Once you have the temporal response of the medium, ${\bf J}_{\rm f}(t,z)$ and/or ${\bf P}(t,z)$, plug them into Ampere's Law and propagate the fields in your simulation according to the resulting wave equation.
