onyx.tex

order发计算光子晶体的能带结构

\documentclass{article}
\pagestyle{myheadings}
\markright{Order N program Documentation}
\oddsidemargin 5mm
\textwidth 15cm
%\textheight 24cm
\begin{document}

\title{Order N program Documentation}
\author{A. J. Ward \\ Imperial College London}

\begin{titlepage}
\maketitle
\thispagestyle{empty}
\begin{abstract}
The following document is intended 
to accompany the Order N method (Finite Difference Time
Domain) electromagnetism program that we have been developing. 
The program is written in Fortran 90.

At the moment the code can perform two basic types
of calculation. The first is  a band structure calculation,
 which starts off with an
arbitrary set of initial fields corresponding to some wavevector $\mathbf{k}$,
integrates in time storing the fields at some random sampling points and then
Fourier transforms into the frequency domain. The peaks in the frequency 
spectrum correspond to the frequencies the the allowed modes.
The second type of calculation finds the time dependent Green's function
by setting the initial fields to be a delta function in space for one of the
field components, zero for all the others. The fields are integrated in time
and stored at each time step to give $g(t,\mathbf{r},\mathbf{r'})$. Fourier
transforming gives $G(\omega,\mathbf{r},\mathbf{r'})$.

Program code is presented in \texttt{typewriter} font. Fortran 90 keywords
 are in \textbf{bold}.
\end{abstract}
\end{titlepage}

\tableofcontents
\newpage

\section{Background Theory}
\subsection{Deriving the discrete Maxwell's equations}

We wish to derive discrete versions of Maxwell's equations which
can be used as the basis for a finite difference time domain calculation
of a form similar to that of Yee~\cite{Yee}.
Starting from Maxwell's equations,
\[
\nabla\times\mathbf{H}=\varepsilon_{0}\varepsilon_{r}
\frac{\partial\mathbf{E}}{\partial t}
\hspace{1cm}
\nabla\times\mathbf{E}=-\mu_{0}\mu_{r}
\frac{\partial\mathbf{H}}{\partial t}
\]

Next we want to use the results for writing Maxwell's equations in
generalised co-ordinates that were developed in~\cite{AJW2} and
explained in detail in chapter four of my thesis~\cite[pages
66-93]{mythesis}.
The basic
results from that work are very simple and mean that in the generalised
co-ordinates can be written,
\[
\nabla_{q}\times\hat{\mathbf{H}}=\varepsilon_{0}\hat{\varepsilon}(\mathbf{r})
\frac{\partial\hat{\mathbf{E}}}{\partial t}
\hspace{1cm}
\nabla_{q}\times\hat{\mathbf{E}}=-\mu_{0}\hat{\mu}(\mathbf{r})
\frac{\partial\hat{\mathbf{H}}}{\partial t}
\]
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}}
\]
The effective permittivity $\hat{\varepsilon}$ and permeability $\hat{\mu}$
can be shown to be tensors of the following form,
\[
\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}}
\]
\[
\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}}
\]
Where the vectors $\mathbf{u}_i$ are unit vectors pointing along the
axes of the generalised co-ordinate system.
The point of this result is that all the details of the co-ordinate
system have been `hidden' inside the effective $\varepsilon$ and $\mu$.
All we need do is ensure that we use $\hat{\varepsilon}^{ij}$ and
 $\hat{\mu}^{ij}$ at each point.

Now we want to write down the discrete versions of these time dependent
equations. To this end we replace the temporal and spacial derivatives
with finite differences. We choose to use forward differences for 
derivatives of the E-field and backwards differences for the H-fields. 
This choice leads to conditionally stable equations.
So\ldots
\[
\frac{\partial\hat{E}}{\partial q_{i}}\mapsto 
\Delta^{+}_{i}\hat{E}=
\hat{E}(\mathbf{r+Q_{i}},t)-\hat{E}(\mathbf{r},t)
\]
\[
\frac{\partial\hat{H}}{\partial q_{i}}\mapsto 
\Delta^{-}_{i}\hat{H}=
\hat{H}(\mathbf{r},t)-\hat{H}(\mathbf{r-Q_{i}},t)
\]
and
\[
\frac{\partial\hat{E}}{\partial t}\mapsto 
\frac{\Delta^{+}_{\tau}\hat{E}}{\delta t}=
\frac{\hat{E}(\mathbf{r},t+\delta t)-\hat{E}(\mathbf{r},t)}{\delta t}
\]
\[
\frac{\partial\hat{H}}{\partial q_{i}}\mapsto 
\frac{\Delta^{-}_{\tau}\hat{H}}{\delta t}=
\frac{\hat{H}(\mathbf{r},t)-\hat{H}(\mathbf{r},t-\delta t)}{\delta t}
\]
Where $\mathbf{Q}_i=Q_i\mathbf{u}_i$.
We can use the discrete spacial derivatives to define a discrete version
of the curl,
\[
\left[\nabla_{q}^{+}\times\mathcal{F}\right]_{3}=
\Delta_{1}^{+}\mathcal{F}_{2}-\Delta_{2}^{+}\mathcal{F}_{1}
\]
etc. Substituting into Maxwell's equations gives,
\[
\nabla_{q}^{+}\times\hat{\mathbf{E}}=-\frac{\mu_{0}}{\delta t}
\hat{\mu}(\mathbf{r})\Delta^{+}_{\tau}\hat{\mathbf{H}}
\]
\[
\nabla_{q}^{-}\times\hat{\mathbf{H}}=+\frac{\varepsilon_{0}}{\delta t}
\hat{\varepsilon}(\mathbf{r})\Delta^{-}_{\tau}\hat{\mathbf{E}}
\]
Here we should note that on the discrete mesh the $Q_i$ factors simply
correspond to the spacings between the mesh points in each direction.

This choice of finite differences is equivalent to making the following
approximations to $\mathbf{k}$ and $\omega$. For the electric field
terms,
\[
k_j\mapsto k^E_j=\frac{\exp{[ik_j Q_j]}-1}{iQ_j}
\hspace{1cm}
\omega\mapsto \omega^E=\frac{\exp{[-i\omega\delta t]}-1}{-i\delta t}
\]
And for the magnetic field terms,
\[
k_j\mapsto k^H_j=\frac{1-\exp{[-ik_j Q_j]}}{iQ_j}
\hspace{1cm}
\omega\mapsto \omega^H=\frac{1-\exp{[i\omega\delta t]}}{-i\delta t}
\]
For the free-space dispersion relation, $\omega/k=c_0$ these 
approximations give,
\[
\frac{4}{\delta t^2}\sin^2\left(\frac{\omega \delta t}{2}\right)=
4c_0^2\left\{
 \frac{1}{Q_1^2}\sin^2\left(\frac{Q_1\ k_x}{2}\right)
+\frac{1}{Q_2^2}\sin^2\left(\frac{Q_2\ k_y}{2}\right)
+\frac{1}{Q_3^2}\sin^2\left(\frac{Q_3\ k_z}{2}\right)\right\}
\]

For reasons of numerical stability it is desirable to make the electric
and magnetic fields of roughly equal magnitude. This can be done by 
introducing a rescaled magnetic field.
\[
\hat{\mathbf{H}}'=\frac{\delta t}{\varepsilon_{0}Q_{0}}\hat{\mathbf{H}}
\]
Where we introduce the constant $Q_0$ to have the dimensions of length
and be roughly equal in length to the $Q_i$'s.
If we also introduce a rescaling to $\hat{\varepsilon}$ and $\hat{\mu}$
then Maxwell's equations become,
\[
\Delta^{+}_{\tau}\hat{\mathbf{E}}(\mathbf{r},t)=
\left[\hat{\varepsilon}(\mathbf{r})\right]^{-1}
\nabla^{-}_{q}\times \hat{\mathbf{H}}'(\mathbf{r},t)
\]
\[
\Delta^{-}_{\tau}\hat{\mathbf{H}}'(\mathbf{r},t)=
-\left\{\frac{\delta t\ c_{0}}{Q_{0}}\right\}^{2}
\left[\hat{\mu}(\mathbf{r})\right]^{-1}
\nabla^{+}_{q}\times \hat{\mathbf{E}}(\mathbf{r},t)
\]
where $\hat{\varepsilon}$ and $\hat{\mu}$ are now,
\[
\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}}
\]
So if we know the electric and magnetic fields at time $t$ we can
calculate the fields at $t+\delta t$
\begin{eqnarray*}
\hat{E}_{1}(\mathbf{r},t+\delta t)=\hat{E}_{1}(\mathbf{r},t)
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{11}\left\{
 \hat{H}'_{3}(\mathbf{r},t)-\hat{H}'_{3}(\mathbf{r-b},t)
-\hat{H}'_{2}(\mathbf{r},t)+\hat{H}'_{2}(\mathbf{r-c},t)\right\}\\
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{12}\left\{
 \hat{H}'_{1}(\mathbf{r},t)-\hat{H}'_{1}(\mathbf{r-c},t)
-\hat{H}'_{3}(\mathbf{r},t)+\hat{H}'_{3}(\mathbf{r-a},t)\right\}\\
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{13}\left\{
 \hat{H}'_{2}(\mathbf{r},t)-\hat{H}'_{2}(\mathbf{r-a},t)
-\hat{H}'_{1}(\mathbf{r},t)+\hat{H}'_{1}(\mathbf{r-b},t)\right\}\\
\end{eqnarray*}
\begin{eqnarray*}
\hat{E}_{2}(\mathbf{r},t+\delta t)=\hat{E}_{2}(\mathbf{r},t)
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{21}\left\{
 \hat{H}'_{3}(\mathbf{r},t)-\hat{H}'_{3}(\mathbf{r-b},t)
-\hat{H}'_{2}(\mathbf{r},t)+\hat{H}'_{2}(\mathbf{r-c},t)\right\}\\
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{22}\left\{
 \hat{H}'_{1}(\mathbf{r},t)-\hat{H}'_{1}(\mathbf{r-c},t)
-\hat{H}'_{3}(\mathbf{r},t)+\hat{H}'_{3}(\mathbf{r-a},t)\right\}\\
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{23}\left\{
 \hat{H}'_{2}(\mathbf{r},t)-\hat{H}'_{2}(\mathbf{r-a},t)
-\hat{H}'_{1}(\mathbf{r},t)+\hat{H}'_{1}(\mathbf{r-b},t)\right\}\\
\end{eqnarray*}
\begin{eqnarray*}
\hat{E}_{3}(\mathbf{r},t+\delta t)=\hat{E}_{3}(\mathbf{r},t)
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{31}\left\{
 \hat{H}'_{3}(\mathbf{r},t)-\hat{H}'_{3}(\mathbf{r-b},t)
-\hat{H}'_{2}(\mathbf{r},t)+\hat{H}'_{2}(\mathbf{r-c},t)\right\}\\
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{32}\left\{
 \hat{H}'_{1}(\mathbf{r},t)-\hat{H}'_{1}(\mathbf{r-c},t)
-\hat{H}'_{3}(\mathbf{r},t)+\hat{H}'_{3}(\mathbf{r-a},t)\right\}\\
&+&\left[\hat{\varepsilon}^{-1}(\mathbf{r})\right]^{33}\left\{
 \hat{H}'_{2}(\mathbf{r},t)-\hat{H}'_{2}(\mathbf{r-a},t)
-\hat{H}'_{1}(\mathbf{r},t)+\hat{H}'_{1}(\mathbf{r-b},t)\right\}\\
\end{eqnarray*}
\begin{eqnarray*}
\hat{H}'_{1}(\mathbf{r},t+\delta t)&=&\hat{H}'_{1}(\mathbf{r},t)
-\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}\times\\
\left[\right.&+&\left.\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{11}\left\{
 \hat{E}_{3}(\mathbf{r+b},t)-\hat{E}_{3}(\mathbf{r},t)
-\hat{E}_{2}(\mathbf{r+c},t)+\hat{E}_{2}(\mathbf{r},t)\right\}\right.\\
&+&\left.
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{12}\left\{
 \hat{E}_{1}(\mathbf{r+c},t)-\hat{E}_{1}(\mathbf{r},t)
-\hat{E}_{3}(\mathbf{r+a},t)+\hat{E}_{3}(\mathbf{r},t)\right\}\right.\\
&+&\left.
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{13}\left\{
 \hat{E}_{2}(\mathbf{r+a},t)-\hat{E}_{2}(\mathbf{r},t)
-\hat{E}_{1}(\mathbf{r+b},t)+\hat{E}_{1}(\mathbf{r},t)\right\}\right]
\end{eqnarray*}
\begin{eqnarray*}
\hat{H}'_{2}(\mathbf{r},t+\delta t)&=&\hat{H}'_{2}(\mathbf{r},t)
-\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}\times\\
\left[\right.&+&\left.\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{21}\left\{
 \hat{E}_{3}(\mathbf{r+b},t)-\hat{E}_{3}(\mathbf{r},t)
-\hat{E}_{2}(\mathbf{r+c},t)+\hat{E}_{2}(\mathbf{r},t)\right\}\right.\\
&+&\left.
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{22}\left\{
 \hat{E}_{1}(\mathbf{r+c},t)-\hat{E}_{1}(\mathbf{r},t)
-\hat{E}_{3}(\mathbf{r+a},t)+\hat{E}_{3}(\mathbf{r},t)\right\}\right.\\
&+&\left.
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{23}\left\{
 \hat{E}_{2}(\mathbf{r+a},t)-\hat{E}_{2}(\mathbf{r},t)
-\hat{E}_{1}(\mathbf{r+b},t)+\hat{E}_{1}(\mathbf{r},t)\right\}\right]
\end{eqnarray*}
\begin{eqnarray*}
\hat{H}'_{3}(\mathbf{r},t+\delta t)&=&\hat{H}'_{3}(\mathbf{r},t)
-\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}\times\\
\left[\right.&+&\left.\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{31}\left\{
 \hat{E}_{3}(\mathbf{r+b},t)-\hat{E}_{3}(\mathbf{r},t)
-\hat{E}_{2}(\mathbf{r+c},t)+\hat{E}_{2}(\mathbf{r},t)\right\}\right.\\
&+&\left.
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{32}\left\{
 \hat{E}_{1}(\mathbf{r+c},t)-\hat{E}_{1}(\mathbf{r},t)
-\hat{E}_{3}(\mathbf{r+a},t)+\hat{E}_{3}(\mathbf{r},t)\right\}\right.\\
&+&\left.
\left[\hat{\mu}^{-1}(\mathbf{r})\right]^{33}\left\{
 \hat{E}_{2}(\mathbf{r+a},t)-\hat{E}_{2}(\mathbf{r},t)
-\hat{E}_{1}(\mathbf{r+b},t)+\hat{E}_{1}(\mathbf{r},t)\right\}\right]
\end{eqnarray*}
where we have written $\mathbf{a}\equiv\mathbf{Q}_1$, 
$\mathbf{b}\equiv\mathbf{Q}_2$ and $\mathbf{c}\equiv\mathbf{Q}_3$.

\subsection{Introducing the current terms}

All that remains now is to include the current terms.
This is done by means of a finite electrical conductivity. For historical
reasons a ficticious magnetic
current in also introduced by including a magnetic conductivity. 
\[
\nabla\times\mathbf{H}=\sigma\mathbf{E}+\varepsilon_{0}\varepsilon_{r}
\frac{\partial\mathbf{E}}{\partial t}
\hspace{1cm}
\nabla\times\mathbf{E}=-\left[\sigma_{m}\mathbf{H}+\mu_{0}\mu_{r}
\frac{\partial\mathbf{H}}{\partial t}\right]
\]
These conductivities are not
really necessary and are really only included for historical reasons.
In the generalised co-ordinates,
\[
\frac{\nabla^{-}_{q}\times \hat{\mathbf{H}}(\mathbf{r},t)}{Q_0}
=\varepsilon_0\hat{\varepsilon}(\mathbf{r})\frac{\Delta^{+}_{\tau}
\hat{\mathbf{E}}(\mathbf{r},t)}{\delta t}
+\hat{\sigma}\hat{\mathbf{E}}(\mathbf{r},t)
\hspace{1cm}
\frac{\nabla^{+}_{q}\times \hat{\mathbf{E}}(\mathbf{r},t)}{Q_0}=
-\mu_0\hat{\mu}(\mathbf{r})
\frac{\Delta^{-}_{\tau}\hat{\mathbf{H}}(\mathbf{r},t)}{\delta t}
-\hat{\sigma}_{m}\hat{\mathbf{H}}(\mathbf{r},t)
\]