On Mon, 7 Aug 2006, Ben Cowan wrote:
There's a stickier issue as well - a PML overlapping a periodic structure
can be unstable: Suppose the mode propagating along the grating has Bloch
wavenumber k, so the amplitude goes as exp(-ikz), assuming positive omega,
if the grating is periodic in the z direction. Now a PML on the +z side
essentially multiplies k by a complex number with a negative imaginary
component, so that the mode attenuates as it propagates in the +z direction.
But in a periodic structure, there can be modes with anomalous dispersion:
the group velocity is in the opposite direction to k. As such a mode
propagates in +z, it has negative k. Then, as it enters the PML, it gets a
positive imaginary component, so it amplifies. Then it reflects off the
boundary behind the PML, and Im(k) is negative, so it amplifies again as it
propagates in -z. Some of the power can reflect off the interface between
the PML and the normal region, so the process feeds back. Modes can build
up from numerical noise this way, leading to late-time instability in the
simulation.
I don't think this analysis is correct. When you add an imaginary part to
the dielectric constant, as PML effectively does (in a
frequency-dependent, anisotropic fashion), the result is either gain or
loss, depending upon the sign of the imaginary part, but independent of
the direction of propagation (or the sign of k). The basic error in your
analysis is that you are using the phase velocity rather than the group
velocity to determine the sign of the imaginary part of k due to the
imaginary change in epsilon. (Equivalently, you are ignoring the z
dependence from the Bloch envelope, in addition to the phase factor.)
There are a couple of ways to see this explicitly from perturbation
theory, supposing that you add a small imaginary part to the dielectric
constant. You can do perturbation theory in the frequency, getting a
small imaginary part in the frequency at a fixed real k. Then, using the
analtic properties of the dispersion relation, you can switch to get the
imaginary part of k at a fixed real frequency, and it is easy to show that
this involves simply dividing by the group velocity (independent of the
sign of k). Or, you can do perturbation theory directly in k, and you
find the same expression because the eigenproblem in k changes the
normalization by a factor of the power over the energy, which again gives
the group velocity.
In fact, this must be the case since the sign of k in a photonic crystal
is completely arbitrary and unphysical---k is ambiguous because you can
always add a reciprocal lattice vector and get the same eigenstate [*].
It may well be that there are long-time numerical instabilities---there
often are in FDTD, in discontinuous structures where the usual Von-Neumann
stability analysis breaks down---but if so this is not the source of them.
Steven
[*] On the other hand, it is possible to get k opposite to the group
velocity in a *uniform* waveguide structure, in which case the sign of k
is not arbitrary. See e.g. Ibanescu et al, PRL 92, 063903 (2004).
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
