1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331

% 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 021111307, 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://abinitio.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{Frequencydependent epsilon}
In the case of a frequencydependent $\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 gainthis 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 locationsthis 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}$ conductivitiesyes, 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}
