onyx2.tex

计算光子晶体能带的程序

in the frequency domain as well as spatially.

\subsection{Defining a plane wave basis set} \label{sec:pw}
A vital element of any transmission and reflection calculation, if we 
wish to be able to calculate the elements of the transmission matrix rather
than just the transmitted intensity at some point, is the ability to project
the transmitted / reflected wavefield onto a complete basis set of plane waves.
Fortunately, we can use the transfer matrix to define just such a basis set.

A plane wave solution to Maxwell's equations has the following form,
\[
\mathbf{E}(\mathbf{r},t)=\mathbf{E}_{0}\exp{[i(\mathbf{k\cdot r}-\omega t)]}
\hspace{2cm}
\mathbf{H}(\mathbf{r},t)=\mathbf{H}_{0}\exp{[i(\mathbf{k\cdot r}-\omega t)]}
\]
So if we substitute into Maxwell's equations on the lattice,
\[
\mathbf{k}^{-}\times \mathbf{H}_{0}=-\varepsilon_{0}\omega^{+}\mathbf{E}_{0}
\hspace{2cm}
\mathbf{k}^{+}\times \mathbf{E}_{0}=\mu_{0}\omega^{-}\mathbf{H}_{0}
\]
The second equation allows us to generate the components of $\mathbf{H}$
given $\mathbf{E}$. 
For a given $\mathbf{k}$ there are two polarisation
states which we must consider. For the s-polarisation,
\[
\mathbf{E}_{s}\propto \mathbf{k}^{-} \times \mathbf{n}
\]
and for the p-polarisation,
\[
\mathbf{E}_{p}\propto \mathbf{k}^{-} \times \mathbf{k}^{+} \times \mathbf{n}
\]
where $\mathbf{n}$ is an arbitrary unit vector. Once we have 
$\mathbf{E}$, finding $\mathbf{H}$ is a straight-forward application of our
discrete Maxwell's equations. At a given frequency $\omega$ we loop over
all $\mathbf{k}_{\parallel}$ to generate the complete set.

The next task is to properly normalise this basis set. This is done by
ensuring the current carried by each mode is normalised to unity. The
z-component of the current can be shown to be,
\[
S_{z}=\Exp{i k_{z}c} E_{x} H_{y}^{\prime *}
     -\Exp{i k_{z}c} E_{y} H_{x}^{\prime *}
    +\Exp{iw \delta t} H_{y}^{\prime}(\Exp{i k_{z}c} E_{x})^{*}
    -\Exp{iw \delta t} H_{x}^{\prime}(\Exp{i k_{z}c} E_{y})^{*}
\]
It is very easy to show that the plane waves which we have just defined are
the right eigenvectors of the transfer matrix for a homogeneous dielectric. 
By definition the effect 
of the transfer matrix on an arbitrary set of fields is,
\[
\hat{T} \mathcal{F}_{r}(\mathbf{r}) = \mathcal{F}_{r}(\mathbf{r+c})
\]
where,
\[
\mathcal{F}_{r} = \left(
\begin{array}{c}
E_{x}\\E_{y}\\H^{\prime}_{x}\\H^{\prime}_{y}
\end{array}
\right)
\]
But for plane waves,
\[
\exp{[ik_{z}c]}\mathcal{F}_{r}(\mathbf{r}) = \mathcal{F}_{r}(\mathbf{r+c})
\]
hence,
\[
\hat{T} \mathcal{F}_{r}(\mathbf{r}) 
= \exp{[ik_{z}c]}\mathcal{F}_{r}(\mathbf{r})
\]
Unfortunately for us, the right eigenvectors of the transfer matrix do not
in general form an orthonormal set. So
if we are to project out the components of particular plane waves from an
arbitrary set of fields, we need a set of vectors which are the reciprocal set
to the right eigenvectors. We could generate the set by inverting the matrix
formed by the right eigenvectors but fortunately there is a simpler way of
doing this. The transfer matrix also has a distinct set of left eigenvectors
and these are the reciprocal set which we need.
\[
\mathcal{F}_{l}(\mathbf{r})\hat{T} 
= \exp{[ik_{z}c]}\mathcal{F}_{l}(\mathbf{r})
\]
and
\[
\mathcal{F}_{l}^{i} \cdot \mathcal{F}_{r}^{j}=\delta^{ij}
\]
By examining the elements of the transfer matrix (which we derived in
section~\ref{sec:tmm}) we can show how to generate the left vectors
from the elements of the right ones.
\[
\mathcal{F}_{l}^{s+}=\left(-E_{y}^{p+},+E_{x}^{p+},+\alpha H_{y}^{\prime p+},
-\alpha H_{x}^{\prime p+}
\right)
\]
\[
\mathcal{F}_{l}^{p+}=\left(-E_{y}^{s+},+E_{x}^{s+},+\alpha H_{y}^{\prime s+},
-\alpha H_{x}^{\prime s+}
\right)
\]
\[
\mathcal{F}_{l}^{s-}=\left(-E_{y}^{p-},+E_{x}^{p-},+\alpha H_{y}^{\prime p-},
-\alpha H_{x}^{\prime p-}
\right)
\]
\[
\mathcal{F}_{l}^{p-}=\left(-E_{y}^{s-},+E_{x}^{s-},+\alpha H_{y}^{\prime s-},
-\alpha H_{x}^{\prime s-}
\right)
\]
where,
\[
\alpha=\left(\frac{Q_{0}}{c_{0}\delta t}\right)^{2}
             \frac{\omega^{-}}{\varepsilon_{\mathrm{r}}\ \omega^{+}}
\]
Note that $s+$ denotes an s-polarised wave coming from $+\infty$, $p-$ denotes 
a p-polarised wave coming from $-\infty$ \textit{etc}. This left eigenvector
set can now be used to project out the individual modes from an arbitrary set
of fields. If we write the arbitrary fields as a sum over modes,
\[
F = \sum_{i} \alpha_{i}\mathcal{F}^{i}_{r}
\]
then,
\[
\alpha_{j} = \mathcal{F}^{j}_{l} \cdot F
\]

\subsection{Defining an absorbing boundary condition} \label{sec:abs}
In this section we will derive formulae for an absorbing 
boundary region which
we can use to pad ends of our real-space lattice and absorb any spurious 
reflections. This absorbing layer is essential if we are to successfully 
perform a transmission/reflection calculation where even small
reflections from the
ends of the system could cause serious errors. 
Figure~\ref{fig:PML} gives a schematic to illustrate the set-up for
performing a transmission/reflection calculation. On one side of the system
some incident field is set up, probably a Gaussian pulse but certainly 
something which contains all of the frequencies which we are interested in.
The fields are then time evolved and at every time step the fields are stored
on the planes indicated with the dashed lines in the figure. At the end of the
integration time the stored fields are Fourier transformed into the frequency
domain and projected onto a complete basis set of plane waves to give the
amplitude of each mode present in the scattered fields.

However, all this relies critically on the minimisation of spurious reflections
from the ends of the computational domain if it to work effectively. To this
end we will implement a Berenger
type of perfectly matched layer (PML)~\cite{Berenger} 
as our absorbing layer and
our derivation will closely follow the work of Zhao and 
Cangellaris~\cite{PML}.
\begin{figure}
\hrule\vspace{5mm}
\begin{center}
\include{PML1}
\end{center}
\hrule
\caption{Diagram showing the typical arrangement for a 
transmission/reflection coefficient calculation. Some incident wave field is
set up and allowed to time evolve to give reflected and transmitted fields.
Absorbing regions (PML's) minimise any unphysical 
reflections from the ends of the system. The fields are stored at each time 
step on the planes shown with dashed lines and analysed to give the 
transmission and reflection coefficients as functions of frequency.}
\label{fig:PML}
\end{figure}

We begin from Maxwell's equations in an arbitrary co-ordinate system,
\[
\frac{\nabla_{q}^{-}\times\hat{\mathbf{H}}}{Q_{0}}=
-i\varepsilon_{0}\hat{\varepsilon}(\mathbf{r})
\omega^{+}\hat{\mathbf{E}}
\hspace{1cm}
\frac{\nabla_{q}^{+}\times\hat{\mathbf{E}}}{Q_{0}}=
 i\mu_{0}\hat{\mu}(\mathbf{r})
\omega^{-}\hat{\mathbf{H}}
\]
Where the renormalised fields are,
\[
\hat{E}_{i}=Q_{i}E_{i}\hspace{2cm}
\hat{H}_{i}=Q_{i}H_{i}
\]
with,
\[
Q_{i}=\sqrt{
 \left(\frac{\partial x}{\partial q_{i}}\right)^{2}
+\left(\frac{\partial y}{\partial q_{i}}\right)^{2}
+\left(\frac{\partial z}{\partial q_{i}}\right)^{2}}
\]
\[
\hat{\varepsilon}^{ij}(\mathbf{r})=\varepsilon(\mathbf{r})
g^{ij}|\mathbf{u_{1}\cdot u_{2}\times u_{3}}|\frac{Q_{1}Q_{2}Q_{3}}
{Q_{i}Q_{j}Q_{0}}
\]
\[
\hat{\mu}^{ij}(\mathbf{r})=\mu(\mathbf{r})
g^{ij}|\mathbf{u_{1}\cdot u_{2}\times u_{3}}|\frac{Q_{1}Q_{2}Q_{3}}
{Q_{i}Q_{j}Q_{0}}
\]
To make a layer which is absorbing we simply make $Q_{3}$ complex and keep 
the co-ordinate axis orthogonal. Typically, if in the free-space region,
\[
Q_1=a \hspace{5mm};\hspace{5mm} Q_2=b \hspace{5mm};\hspace{5mm} Q_3=c
\]
then in the absorbing region,
\[
Q_1=a \hspace{5mm};\hspace{5mm} Q_2=b \hspace{1cm};\hspace{1cm}
Q_3=c\ \left(1+\frac{i\omega_{z}}{\omega^{-}}\right)
\]
where $\omega_{z}$ is an arbitrary parameter which we can use to control the
strength of the absorption. Note that the absorption in $Q_{3}$ must be 
odd with respect to frequency so that it changes sign at $\omega=0$. This 
ensures that all solutions are absorbed, even non-physical negative frequency
solutions. Note also that we could equally well have chosen to use $\omega^{+}$
in the above equation just as long as we were consistent.
Substituting into $\hat{\varepsilon}$ and $\hat{\mu}$ this gives in the free
space region,
\[
\hat{\varepsilon}^{ii}=\varepsilon_{\mathrm{r}} 
\frac{Q_{1}Q_{2}Q_{3}}{Q_{i}^{2}Q_{0}}
\hspace{2cm}
\hat{\mu}^{ii}=\frac{Q_{1}Q_{2}Q_{3}}{Q_{i}^{2}Q_{0}}
\]
and in the absorber,
\[
\hat{\varepsilon}^{11}=
\varepsilon_{\mathrm{r}}\frac{b\,c}{a\,Q_{0}}
\left(1+\frac{i\omega_{z}}{\omega^{-}}\right)
\hspace{5mm};\hspace{5mm}
\hat{\varepsilon}^{22}=
\varepsilon_{\mathrm{r}}\frac{a\,c}{b\,Q_{0}}
\left(1+\frac{i\omega_{z}}{\omega^{-}}\right)
\hspace{5mm};\hspace{5mm}
\hat{\varepsilon}^{33}=
\varepsilon_{\mathrm{r}}\frac{a\,b}{c\,Q_{0}}
\left(1+\frac{i\omega_{z}}{\omega^{-}}\right)^{-1}
\]
\[
\hat{\mu}^{11}=\frac{b\,c}{a\,Q_{0}}
\left(1+\frac{i\omega_{z}}{\omega^{-}}\right)
\hspace{5mm};\hspace{5mm}
\hat{\mu}^{22}=\frac{a\,c}{b\,Q_{0}}
\left(1+\frac{i\omega_{z}}{\omega^{-}}\right)
\hspace{5mm};\hspace{5mm}
\hat{\mu}^{33}=\frac{a\,b}{c\,Q_{0}}
\left(1+\frac{i\omega_{z}}{\omega^{-}}\right)^{-1}
\]
Because the co-ordinate system remains orthogonal, the off-diagonal elements
remain zero. 
To simplify we write,
\[
\tilde{\varepsilon}^{11}=\varepsilon_{r}\,\frac{b\,c}{a\,Q_{0}}
\hspace{5mm};\hspace{5mm}
\tilde{\varepsilon}^{22}=\varepsilon_{r}\,\frac{a\,c}{b\,Q_{0}}
\hspace{5mm};\hspace{5mm}
\tilde{\varepsilon}^{33}=\varepsilon_{r}\,\frac{a\,b}{c\,Q_{0}}
\]
\[
\tilde{\mu}^{11}=\frac{b\,c}{a\,Q_{0}}
\hspace{5mm};\hspace{5mm}
\tilde{\mu}^{22}=\frac{a\,c}{b\,Q_{0}}
\hspace{5mm};\hspace{5mm}
\tilde{\mu}^{33}=\frac{a\,b}{c\,Q_{0}}
\]
When we substitute into Maxwell's equations we obtain the 
following,
\[
\left[\nabla_{q}\times\hat{\mathbf{H}}\right]_{x}=-i\varepsilon_{0}
\tilde{\varepsilon}^{11}
\left(1+\frac{i\omega_{z}}{\omega^{-}}\right)
\omega^{+}\hat{E_{x}}
\hspace{1cm};\hspace{1cm}
\left[\nabla_{q}\times\hat{\mathbf{E}}\right]_{x}=i\mu_{0}\tilde{\mu}^{11}
\left(1+\frac{i\omega_{z}}{\omega^{-}}\right)
\omega^{-}\hat{H_{x}}
\]
\[
\left[\nabla_{q}\times\hat{\mathbf{H}}\right]_{y}=-i\tilde{\varepsilon}^{22}
\left(1+\frac{i\omega_{z}}{\omega^{-}}\right)
\omega^{+}\hat{E_{y}}
\hspace{1cm};\hspace{1cm}
\left[\nabla_{q}\times\hat{\mathbf{E}}\right]_{y}=i\tilde{\mu}^{22}
\left(1+\frac{i\omega_{z}}{\omega^{-}}\right)
\omega^{-}\hat{H_{y}}
\]
\[
\left[\nabla_{q}\times\hat{\mathbf{H}}\right]_{z}=-i\tilde{\varepsilon}^{33}
\left(1+\frac{i\omega_{z}}{\omega^{-}}\right)^{-1}
\omega^{+}\hat{E_{z}}
\hspace{1cm};\hspace{1cm}
\left[\nabla_{q}\times\hat{\mathbf{E}}\right]_{z}=i\tilde{\mu}^{33}
\left(1+\frac{i\omega_{z}}{\omega^{-}}\right)^{-1}
\omega^{-}\hat{H_{z}}
\]
Clearly, what we have done is releated to pushing the frequencies off the real
axis by a small ammount determined by $\omega_{z}$.
Maxwell can now be rearranged to give,
\[
\hat{E_{x}}(\mathbf{r},t+\delta t)=\frac{1}{1+\omega_{z} \delta t}
\left\{
\hat{E_{x}}(\mathbf{r},t)+\frac{1}{\tilde{\varepsilon}^{11}}
\left[\nabla_{q}\times\hat{\mathbf{H}}'(\mathbf{r},t)
\right]_{x}