File: meep.tex

package info (click to toggle)
meep 0.10-2.1
  • links: PTS
  • area: main
  • in suites: lenny
  • size: 4,380 kB
  • ctags: 5,469
  • sloc: cpp: 50,653; sh: 8,380; haskell: 744; makefile: 367; perl: 10
file content (331 lines) | stat: -rw-r--r-- 13,775 bytes parent folder | download
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 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}