meep.tex

来自mit的fdtd开放性源代码meep

%  Copyright (C) 2003 David Roundy
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2, or (at your option)
%  any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software Foundation,
%  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  
\documentclass[floats]{book}
\usepackage{color}
\usepackage{graphicx}
\usepackage{amsmath, amsthm, amssymb}
\usepackage{hyperref}
\usepackage{verbatim}

\begin{document}

% Definition of title page:
\title{
    Meep\\
{\Large \it MIT Electromagnetic Equation Propogation}\\
{\Large (this is the old and somewhat obsolete {\LaTeX} manual)}\\
{\Large ---See http://ab-initio.mit.edu/meep/doc for the latest documentation---}
}
\author{
    David Roundy, Mihai Ibanescu, Peter Bermel, \\
    Steven G. Johnson, and Ardavan Farjadpour
}

\maketitle 

\tableofcontents

\chapter{Discretizing Maxwell's equations}

Maxwell's equations in the absence of sources are:
\begin{align}
\frac{d\mathbf H}{dt} &= -c \mathbf \nabla \times \mathbf E\\
\frac{d\mathbf D}{dt} &= c \mathbf \nabla \times \mathbf H
\end{align}
If the material is a simple isotropic dielectric, we can simply write
$\mathbf D = \epsilon \mathbf E$ and get on with our lives.  Alas, all too
often this is not the case!  We need to be able to deal with anisotropic
dielectrics in which $\mathbf \epsilon$ is a tensor quantity, nonlinear
materials in which $\mathbf \epsilon$ is a function of $\mathbf E$ itself,
and polaritonic and polaronic materials in which $\mathbf \epsilon$ is a
function of frequency.

\subsection{Frequency-dependent epsilon}

In the case of a frequency-dependent $\mathbf \epsilon$, we write
\begin{equation}
\mathbf D = \mathbf \epsilon_{\infty} \mathbf E + \mathbf P
\end{equation}
where $\mathbf P$ is the polarization as a function of time associated with
the frequency dependence of $\mathbf \epsilon$.  Actually, in general there
will be a set of polarizations, and we'll need a summation here.  For
simplicity we'll only describe the case of a single polarization in this
section.  The time dependence of a single polarization is given by
\begin{equation}
\frac{d^2\mathbf{P}}{dt^2} + \gamma \frac{d\mathbf{P}}{dt}
+ \omega_0^2 \mathbf{P} = \Delta\epsilon\ \omega_0^2 \mathbf{E}
\end{equation}
where $\gamma$, $\omega_0$ and $\Delta\epsilon$ are material parameters.
The energy lost due to the absorption by this resonance is simply
\begin{equation}
\Delta U = \mathbf P \frac{d\mathbf E}{dt}
\end{equation}
In fact, if one sets $\Delta\epsilon$ to be negative, we can model gain
effectively in this way, and in this case keeping track of the energy
allows us to model a situation in which there is a depleteable population
inversion which is causing the gain---this is the situation of gain with
saturation.

\subsection{Nonlinear dielectrics}

In nonlinear dielectrics $\mathbf D$ is typically given by a cubic function
of $\mathbf E$.
\begin{equation}
  \mathbf D = \left(\epsilon + \xi \left|\mathbf E\right|^2 \right)\mathbf E
\end{equation}

\subsection{Anisotropic dielectrics}
In anisotropic dielectrics the dielectric constant is a tensor quantity
rather than a scalar quantity.  In this case we write (FIXME: how to do a
tensor in latex?)
\begin{align}
  \mathbf D = \bar{\bar \epsilon} \mathbf E\\
  \mathbf E = \bar{\bar \epsilon}^{-1} \mathbf D\\
\end{align}

\subsection{Putting it all together}

Putting it all together, we get a simplified time stepping of something
like
\begin{align}
\frac{d\mathbf H}{dt} &= -c \mathbf \nabla \times \mathbf E\\
\frac{d\mathbf E}{dt} &= \left( \bar{\bar{\epsilon}}_\infty +
                                \mathbf E \cdot \bar{\bar \xi} \cdot \mathbf E
                         \right)^{-1}
  \left(
  c \mathbf \nabla \times \mathbf H
  - \sum_i \frac{d\mathbf P_i}{dt}
  \right)
\end{align}
\begin{align}
\frac{d^2\mathbf{P}_i}{dt^2} + \gamma_i \frac{d\mathbf{P}_i}{dt}
+ {\omega_0}_i^2 \mathbf{P}_i = \Delta\epsilon_i\ {\omega_0}_i^2 \mathbf{E}
\end{align}

\section{The Yee lattice}

In discretizing Maxwell's equations, we need to put $\mathbf E$ and
$\mathbf H$ on a grid.  Because we only need to calculate the curl of these
quantities, we only need to know them at limited locations--this gives us
the accuracy of a fine grid while only requiring as much data as a grid
twice as coarse.  This trick is called the Yee lattice.
Figure~\ref{yee_fig} shows the Yee lattice in cylindrical coordinates (with
$\hat z$ being to the right).  The gray squares indicate the locations at
which $\epsilon$ is stored.

\begin{figure}
\caption{Yee lattice in cylindrical coordinates.
\label{yee_fig}}
\centering
\includegraphics[width=7.8cm,clip=true]{Yee_bulk}
\end{figure}

The Yee lattice has the property that all the derivatives needed for
$\mathbf \nabla \times \mathbf H$ are known at the Yee lattice points of
$\mathbf E$.  For example, if you look at the $H_\phi$ location.
$\frac{dH_\phi}{dt}$ depends on $\frac{dE_r}{dz}$ and $\frac{dE_z}{dr}$.
This is great, because $E_r$ is known to the right and left of $H_\phi$,
and $E_z$ is known above and below $H_\phi$.

The same principle that the Yee lattice does with space, we also do with
time.  $\mathbf E$ and $\mathbf H$ are known at different times, so that
the time derivative of $\mathbf E$ is known at a $\mathbf H$ time and vice
versa.

\section{Maxwell's equations in cylindrical coordinates}

Here are Maxwell's equations in cylindrical coordinates.  We take the
fields to be of the form:
\begin{equation*}
\mathbf{E}(r,\phi,z) = \mathbf{E}_m(r,z)e^{i m \phi} 
\end{equation*}

Without further ado:
\begin{align}
\frac1c\frac{dH_r}{dt} &= \frac{dE_\phi}{dz} - \frac{im}r E_z\\
\frac1c\frac{dH_\phi}{dt} &= \frac{dE_z}{dr} - \frac{dE_r}{dz}\\
\frac1c\frac{dH_z}{dt} &= \frac{im}r E_r - \frac1r\frac{d(rE_\phi)}{dr}
\end{align}
\begin{align}
\frac\epsilon c\frac{dE_r}{dt} &= \frac{im}r H_z - \frac{dH_\phi}{dz} \\
\frac\epsilon c\frac{dE_\phi}{dt} &= \frac{dH_r}{dz} - \frac{dH_z}{dr} \\
\frac\epsilon c\frac{dE_z}{dt} &= \frac1r\frac{d(rH_\phi)}{dr} - \frac{im}r H_r
\end{align}


\chapter{PML}

PML (Perfectly Matched Layers) is used to provide absorbing boundary
conditions in either the $z$ or $r$ direction.  PML consists of a material
in which some of the field components are split into two fields, each of
which has a conductivity associated with it, which is responsible for the
absorption of the PML.

PML is a sort of material that contains a set of conductivities $\sigma_r$,
$\sigma_\phi$ and $\sigma_z$.  These conductivities are both $\mathbf{E}$
and $\mathbf{H}$ conductivities---yes, we have magnetic monopoles moving
around in our PML.  $\ddot\smile$ Each $\sigma$ causes absorption of
radiation in the direction it is named after.  Thus $\sigma_\phi$ is small,
and almost unnecesary, and is only needed because of the curvature of the
radial surface.  The value of $\sigma_\phi$ at a given radius is equal to
\begin{equation}
\sigma_\phi(r) = \frac1r \int_0^r \sigma_r(r)dr
\end{equation}

If we had a IDTD (Infinitesimal Difference Time Domain) code, PML would be
perfectly absorbing, regardless of the variation of $\sigma$ with position.
However, since meep is a lowly FDTD code, we have to make sure that
$\sigma$ varies only slowly from one grid point to the next.  We do this by
making $\sigma_z$ (for example) vary as $z^2$, with a maximum value of
$\sigma_{max}$ right in front of the boundary.  At the edge of the PML
region is a metalic boundary condition.  The optimal value of
$\sigma_{max}$ is determined by a tradeoff between reflection off the
metallic boundary, caused by too little a $\sigma_{max}$, and reflection
off the sigma itself, caused by too large a $\sigma_{max}$, which makes for
a large variation of $\sigma$ from one grid point to the next.

Here are the field equations for a PML material:
\begin{align}
\frac{dH_{r\phi}}{dt} &= - c \frac{im}r E_z             - \sigma_\phi H_{r\phi} &
\frac{dH_{rz}}{dt} &= c \frac{dE_\phi}{dz}              - \sigma_z H_{rz}\\
\frac{dH_{\phi z}}{dt} &= - c \frac{dE_r}{dz}           - \sigma_z H_{\phi z} &
\frac{dH_{\phi r}}{dt} &= c \frac{dE_z}{dr}             - \sigma_r H_{\phi r} \\
\frac{dH_{zr}}{dt} &= - c \frac1r\frac{d(rE_\phi)}{dr}  - \sigma_r H_{zr}  &
\frac{dH_{z\phi}}{dt} &= c \frac{im}r E_r               - \sigma_\phi H_{z\phi} \\
\epsilon\frac{dE_{r\phi}}{dt} &=   c \frac{im}r H_z             - \sigma_\phi E_{r\phi} &
\epsilon\frac{dE_{rz}}{dt} &= -c\frac{dH_\phi}{dz}              - \sigma_z E_{rz}\\
\epsilon\frac{dE_{\phi z}}{dt} &=   c \frac{dH_r}{dz}           - \sigma_z E_{\phi z} &
\epsilon\frac{dE_{\phi r}}{dt} &=-c \frac{dH_z}{dr}             - \sigma_r E_{\phi r} \\
\epsilon\frac{dE_{zr}}{dt} &=   c \frac1r\frac{d(rH_\phi)}{dr}  - \sigma_r E_{zr}  &
\epsilon\frac{dE_{z\phi}}{dt} &=-c \frac{im}r H_r               - \sigma_\phi E_{z\phi} 
\end{align}

\chapter{Of polaritons and plasmons}\label{polaritons}

Most real materials, at least in some frequency range, have polarizations
that are not actually instantaneously proportional to the local electric
field.  We model these polaritonic and plasmonic effects by introducing one
or more additional polarization fields, to be propogated along with the
electric and magnetic field.  The polarization field, $\mathbf{P}$, is a
vector field which exists on the electric field Yee lattice points.

The polarization field obeys a second order differential equation, which
means that we need to keep track of the polarization at two time steps, in
order to integrate it.

\begin{equation}
\frac{d^2\mathbf{P}}{dt^2} + \gamma \frac{d\mathbf{P}}{dt}
+ \omega^2 \mathbf{P} = \Delta\epsilon\ \omega^2 \mathbf{E}
\end{equation}

To this, we need add one more term to maxwell's equation for $\mathbf E$:

\begin{equation}
c \nabla \times \mathbf{H} = \epsilon_\infty \frac{d\mathbf{E}}{dt}
 + \frac{d\mathbf{P}}{dt}
\end{equation}

So far, the polarization is beautifully simple.  However, we would love to
be able to put polaritonic materials into our PML regions, and
unfortunately in the PML region the electric field has been split into two
components, so we need to figure out which of the two components gets the
contribution from $\frac{d\mathbf{P}}{dt}$.  The obvious solution to this
(well, maybe not exactly obvious, but it is the solution) is to split the
polarization field also into two components in the PML region, just as we
split the electric and magnetic fields.

The electric field propogation equations in PML then become:
\begin{align}
\epsilon\frac{dE_{r\phi}}{dt} &=   c \frac{im}r H_z
             - \sigma_\phi E_{r\phi} - \frac{dP_{r\phi}}{dt} \\
\epsilon\frac{dE_{\phi z}}{dt} &=   c \frac{dH_r}{dz}
             - \sigma_z E_{\phi z} - \frac{dP_{\phi z}}{dt} \\
\epsilon\frac{dE_{zr}}{dt} &=   c \frac1r\frac{d(rH_\phi)}{dr}
             - \sigma_r E_{zr} - \frac{dP_{zr}}{dt}
\end{align}
\begin{align}
\epsilon\frac{dE_{rz}}{dt} &= -c\frac{dH_\phi}{dz}
             - \sigma_z E_{rz} - \frac{dP_{rz}}{dt}\\
\epsilon\frac{dE_{\phi r}}{dt} &=-c \frac{dH_z}{dr}
             - \sigma_r E_{\phi r} - \frac{dP_{\phi r}}{dt} \\
\epsilon\frac{dE_{z\phi}}{dt} &=-c \frac{im}r H_r \label{polariton_pml}
             - \sigma_\phi E_{z\phi} - \frac{dP_{z\phi}}{dt}
\end{align}

\chapter{Hints for writing finite difference time domain code}

(Or \emph{Things I forgot many times, so I wrote down so maybe I won't make
  the same mistake again.})

There is just one rule to remember when writing time domain code, and that
is (as Lefteris has repeatedly told me) ``Always know \emph{when} each
equation is evaluated.''  The trick, of course, lies in knowing how to
apply this rule, and remembering to actually apply it (and I think the
latter is perhaps harder than the former).

As an example, I'll convert a PML polariton equation into a finite
difference equation taken from equation~\ref{polariton_pml} of
chapter~\ref{polaritons}.
\begin{equation*}
\epsilon\frac{dE_{z\phi}}{dt} = -c \frac{im}r H_r
             - \sigma_\phi E_{z\phi} - \frac{dP_{z\phi}}{dt}
\end{equation*}
If we consider the $\mathbf{E}$ timesteps to be at times $n$, $n+1$ etc.,
then this equation needs to be evaluated at time $n+\frac12$.  This is no
problem for most of the terms, but it means that the $\sigma_\phi
E_{z\phi}$ term needs to be an average of its values at time $n$ and
$n+1$.  In short (taking $\Delta t$ to be unity)...
\begin{equation*}
\epsilon (E_{z\phi}^{n+1} - E_{z\phi}^n) = -c \frac{im}r H_r^{n+\frac12}
  - \sigma_\phi ( E_{z\phi}^{n+1} + E_{z\phi}^n) - (dP_{z\phi}^{n+1} - dP_{z\phi}^n)
\end{equation*}
Simplifying a tad gives
\begin{equation*}
E_{z\phi}^{n+1} - E_{z\phi}^n = \frac1{\epsilon + \frac12\sigma_\phi}
    \left(-c \frac{im}r H_r^{n+\frac12}
  - \sigma_\phi E_{z\phi}^n - (dP_{z\phi}^{n+1} - dP_{z\phi}^n)\right)
\end{equation*}
Basically, that is all there is to it.  You now have the equation to
determine $E_{z\phi}^{n+1}$ from $E_{z\phi}^n$, $\frac{im}r
H_r^{n+\frac12}$, $dP_{z\phi}^{n+1}$ and $dP_{z\phi}^n$.

\chapter{Tutorial}

\input{simple.tex}

\input{complicated.tex}

\input{simplebands.tex}

\input{omniguide.tex}

\input{polaritonbands.tex}

\input{energy_cons.tex}

\input{energy_cons_1d.tex}

\input{epsilon_polariton_1d.tex}

\input{lossgain_epsilon.tex}

\input{nonlinear.tex}

\appendix

\input{gpl.tex}

\end{document}