## File: meep.tex

package info (click to toggle)
meep 0.10-2.1
 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331 % 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 \mathbf D = \mathbf \epsilon_{\infty} \mathbf E + \mathbf P 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 \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} where $\gamma$, $\omega_0$ and $\Delta\epsilon$ are material parameters. The energy lost due to the absorption by this resonance is simply \Delta U = \mathbf P \frac{d\mathbf E}{dt} 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$. \mathbf D = \left(\epsilon + \xi \left|\mathbf E\right|^2 \right)\mathbf E \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 \sigma_\phi(r) = \frac1r \int_0^r \sigma_r(r)dr 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. \frac{d^2\mathbf{P}}{dt^2} + \gamma \frac{d\mathbf{P}}{dt} + \omega^2 \mathbf{P} = \Delta\epsilon\ \omega^2 \mathbf{E} To this, we need add one more term to maxwell's equation for $\mathbf E$: c \nabla \times \mathbf{H} = \epsilon_\infty \frac{d\mathbf{E}}{dt} + \frac{d\mathbf{P}}{dt} 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}