onyx2.tex

计算光子晶体能带的程序

\[
\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)
\]
where,
\[
\hat{\sigma}^{ij}(\mathbf{r})=\sigma(\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}}
\hspace{1cm}
\hat{\sigma}_{m}^{ij}(\mathbf{r})=\sigma_m(\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,
\[
\Delta^{+}_{\tau}\hat{\mathbf{E}}(\mathbf{r},t)=
-\frac{\delta t \sigma}{\varepsilon_{0}\varepsilon(\mathbf{r})}
\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)=
-\frac{\delta t \sigma_{m}}{\mu_{0}\mu(\mathbf{r})}
\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)
\]
Writing 
\[
\hat{\sigma}=\frac{\delta t \sigma}{\varepsilon_{0}\varepsilon(\mathbf{r})}
\] 
and 
\[
\hat{\sigma}_{m}=\frac{\delta t \sigma_{m}}{\mu_{0}\mu(\mathbf{r})}
\]
Then,
\begin{eqnarray*}
\hat{E}_{1}(\mathbf{r},t+\delta t)=\left[1-\hat{\sigma}\right]
\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)=\left[1-\hat{\sigma}\right]
\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)=\left[1-\hat{\sigma}\right]
\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)&=&
\left[1+\hat{\sigma}_{m}\right]^{-1}\left[
\hat{H}'_{1}(\mathbf{r},t)\right.\\
&-&\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\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\}\\
&-&\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\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\}\\
&-&\left.\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\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)&=&
\left[1+\hat{\sigma}_{m}\right]^{-1}\left[
\hat{H}'_{2}(\mathbf{r},t)\right.\\
&-&\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\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\}\\
&-&\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\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\}\\
&-&\left.\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\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)&=&
\left[1+\hat{\sigma}_{m}\right]^{-1}\left[
\hat{H}'_{3}(\mathbf{r},t)\right.\\
&-&\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\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\}\\
&-&\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\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\}\\
&-&\left.\left(\frac{\delta t\ c_{0}}{Q_{0}}\right)^{2}
\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*}
These are the update equations for the fields that we require.

\subsection{Stability Criterion}
These equations give a stable updating procedure for the fields if
the time step is kept sufficiently small. The criterion is easy to find.
Starting from the dispersion on the discrete mesh,
\[
\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\}
\]
The condition that the maximum value of the right hand side must
correspond to a real frequency gives,
\[
(\delta t)^2<\left(\frac{c_0^2}{Q_1^2}+\frac{c_0^2}{Q_2^2}
+\frac{c_0^2}{Q_3^2}\right)^{-1}
\]

\subsection{The real-space, discrete frequency transfer matrix}\label{sec:tmm}
The transfer matrix~\cite{J+A,JBP} is an operator which given the electric
 and magnetic
fields in one layer of a discrete lattice allows us to find the fields in the
adjacent layer. It is normally defined on a lattice which is discrete in space
but continuous in the frequency domain, that is to say no approximations are
made to the frequency $\omega$. However, for our purposes it is useful to
re-derive the transfer matrix for the finite difference time-domain case where
the frequency $\omega$ is also discrete. This is simply 
achieved by replacing the frequency $\omega$ in Maxwell's equations with the
relevant approximations and then deriving the transfer matrix as before.
Note that we will only be concerned with the transfer matrix in a uniform
dielectric with dielectric constant $\varepsilon_{\mathrm{r}}$.

We begin, as usual, from Maxwell's equations,
\[
\mathbf{k}\times\mathbf{E}=+\mu_{0}\omega\mathbf{H}
\hspace{1cm};\hspace{1cm}
\mathbf{k}\times\mathbf{H}=-\varepsilon_{0}\varepsilon_{\mathrm{r}}
\omega\mathbf{E}
\]
then we make approximations for $\mathbf{k}$ and $\omega$.
\[
\mathbf{k}^{+}\times\mathbf{E}=+\mu_{0}\omega^{-}\mathbf{H}
\hspace{1cm};\hspace{1cm}
\mathbf{k}^{-}\times\mathbf{H}=-\varepsilon_{0}\varepsilon_{\mathrm{r}}
\omega^{+}\mathbf{E}
\]
where, as before, for the electric field terms,
\[
k_j\mapsto k^{+}_j=\frac{\exp{[ik_j Q_j]}-1}{iQ_j}
\hspace{1cm}
\omega\mapsto \omega^{+}=\frac{\exp{[-i\omega\delta t]}-1}{-i\delta t}
\]
And for the magnetic field terms,
\[
k_j\mapsto k^{-}_j=\frac{1-\exp{[-ik_j Q_j]}}{iQ_j}
\hspace{1cm}
\omega\mapsto \omega^{-}=\frac{1-\exp{[i\omega\delta t]}}{-i\delta t}
\]
Next we rescale the magnetic field,
\[
\mathbf{H}'=\frac{\delta t}{\varepsilon_{0}Q_{0}}\mathbf{H}
\]
so that,
\[
\mathbf{k}^{+}\times\mathbf{E}=+\frac{Q_{0}}{c_{0}^{2}\delta t}
\omega^{-}\mathbf{H}'
\hspace{1cm};\hspace{1cm}
\mathbf{k}^{-}\times\mathbf{H}'=-\frac{\delta t \varepsilon_{\mathrm{r}}}
{Q_{0}}\omega^{+}\mathbf{E}
\]
If we now write out component by component,
\[
E_{z}=-\frac{Q_{0}}{\delta t\  \varepsilon_{\mathrm{r}}\ \omega^{+}}
\left[-k_{y}^{-}H_{x}'+k_{x}^{-}H_{y}'
\right]
\hspace{1cm};\hspace{1cm}
H_{z}'=\frac{c_{0}^{2}\delta t}{Q_{0} \omega^{-}}\left[
-k_{y}^{+}E_{x}+k_{x}^{+}E_{y}\right]
\]
\[
k_{z}^{+}E_{x}=k_{x}^{+}E_{z} + \frac{Q_{0}}{c_{0}^{2}\delta t}\omega^{-}H_{y}'
\hspace{1cm};\hspace{1cm}
k_{z}^{-}H_{x}'=k_{x}^{-}H_{z}' - \frac{\delta t\ \varepsilon_{\mathrm{r}}}
{Q_{0}}\omega^{+}E_{y}
\]
\[
k_{z}^{+}E_{y}=k_{y}^{+}E_{z} - \frac{Q_{0}}{c_{0}^{2}\delta t}\omega^{-}H_{x}'
\hspace{1cm};\hspace{1cm}
k_{z}^{-}H_{y}'=k_{y}^{-}H_{z}' + \frac{\delta t\ \varepsilon_{\mathrm{r}}}
{Q_{0}}\omega^{+}E_{x}
\]
We eliminate the z-components to leave us with a set of four equations,
\[
k_{z}^{+}E_{x}=k_{x}^{+}\left[-\frac{Q_{0}}{\delta t \ 
\varepsilon_{\mathrm{r}}\ \omega^{+}}
\left\{-k_{y}^{-}H_{x}'+k_{x}^{-}H_{y}'\right\}
\right]+\frac{Q_{0}}{c_{0}^{2}\delta t}\omega^{-}H_{y}'
\]
\[
k_{z}^{+}E_{y}=k_{y}^{+}\left[-\frac{Q_{0}}{\delta t \ 
\varepsilon_{\mathrm{r}} \ \omega^{+}}
\left\{-k_{y}^{-}H_{x}'+k_{x}^{-}H_{y}'\right\}
\right]-\frac{Q_{0}}{c_{0}^{2}\delta t}\omega^{-}H_{x}'
\]
\[
k_{z}^{-}H_{x}'=k_{x}^{-}\left[\frac{c_{0}^{2}\delta t}{Q_{0} \omega^{-}}
\left\{-k_{y}^{+}E_{x}+k_{x}^{+}E_{y}\right\}
\right]-\frac{\delta t\ \varepsilon_{\mathrm{r}}}{Q_{0}}\omega^{+}E_{y}
\]
\[
k_{z}^{-}H_{y}'=k_{y}^{-}\left[\frac{c_{0}^{2}\delta t}{Q_{0} \omega^{-}}
\left\{-k_{y}^{+}E_{x}+k_{x}^{+}E_{y}\right\}
\right]+\frac{\delta t\ \varepsilon_{\mathrm{r}}}{Q_{0}}\omega^{+}E_{x}
\]
Expanding the brackets and simplifying gives,
\[
E_{x}(\mathbf{k_{\parallel}},z+c,\omega)=E_{x}(\mathbf{k_{\parallel}},z,\omega)
+\left[\frac{iQ_{0}^{2}}{\delta t\ \varepsilon_{\mathrm{r}}\ 
 \omega^{+}}k_{x}^{+}k_{y}^{-}\right]
H_{x}'(\mathbf{k_{\parallel}},z,\omega)
-\left[\frac{iQ_{0}^{2}}{\delta t\ \varepsilon_{\mathrm{r}}\ 
\omega^{+}}k_{x}^{+}k_{x}^{-}
      -\frac{iQ_{0}^{2}\omega^{-}}{c_{0}^{2}\delta t}\right]
H_{y}'(\mathbf{k_{\parallel}},z,\omega)
\]
\[
E_{y}(\mathbf{k_{\parallel}},z+c,\omega)=E_{y}(\mathbf{k_{\parallel}},z,\omega)
+\left[\frac{iQ_{0}^{2}}{\delta t\ \varepsilon_{\mathrm{r}}\ 
\omega^{+}}k_{y}^{+}k_{y}^{-}
      -\frac{iQ_{0}^{2}\omega^{-}}{c_{0}^{2}\delta t}\right]
H_{x}'(\mathbf{k_{\parallel}},z,\omega)
-\left[\frac{iQ_{0}^{2}}{\delta t\ \varepsilon_{\mathrm{r}}\ 
\omega^{+}}k_{y}^{+}k_{x}^{-}\right]
H_{y}'(\mathbf{k_{\parallel}},z,\omega)
\]
\[
H_{x}'(\mathbf{k_{\parallel}},z+c,\omega)=
H_{x}'(\mathbf{k_{\parallel}},z,\omega)
-\left[\frac{i c_{0}^{2}\delta t}{\omega^{-}}k_{x}^{-}k_{y}^{+}
\right]E_{x}(\mathbf{k_{\parallel}},z+c,\omega)
+\left[
\frac{i c_{0}^{2}\delta t}{\omega^{-}}
k_{x}^{+}k_{x}^{-} -i\delta t\ \varepsilon_{\mathrm{r}}\ \omega^{+}\right]
E_{y}(\mathbf{k_{\parallel}},z+c,\omega)
\]
\[
H_{y}'(\mathbf{k_{\parallel}},z+c,\omega)=
H_{y}'(\mathbf{k_{\parallel}},z,\omega)
+\left[i\delta t\ \varepsilon_{\mathrm{r}}\ \omega^{+}
-\frac{i c_{0}^{2}\delta t}{\omega^{-}}k_{y}^{+}k_{y}^{-}\right]
E_{x}(\mathbf{k_{\parallel}},z+c,\omega)
+\left[\frac{i c_{0}^{2}\delta t}{\omega^{-}} k_{x}^{+}k_{y}^{-}\right]
E_{y}(\mathbf{k_{\parallel}},z+c,\omega)
\]
These equations define the transfer matrix for a system which is discrete