Computer Physics Communications 181 (2010) 687–702

Contents lists available at ScienceDirect

Computer Physics Communications
www.elsevier.com/locate/cpc

Meep: A ﬂexible free-software package for electromagnetic simulations by the
FDTD method ✩Ardavan F. Oskooi a,c,∗ , David Roundy b , Mihai Ibanescu a,c,d , Peter Bermel c , J.D. Joannopoulos a,c,d ,Steven G. Johnson a,c,e,∗∗aCenter for Materials Science & Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, United States
Department of Physics, Oregon State University, Corvallis, OR 97331, United StatescResearch Laboratory of Electronics, Massachusetts Institute of Technology, Cambridge, MA 02139, United StatesdDepartment of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, United Statese
Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139, United Statesba r t i c l e

i n f oa b s t r a c t

Article history:
Received 2 October 2009
Received in revised form 13 November 2009
Accepted 17 November 2009
Available online 20 November 2009
Keywords:
Computational electromagnetism
FDTD
Maxwell solverThis paper describes Meep, a popular free implementation of the ﬁnite-difference time-domain (FDTD)
method for simulating electromagnetism. In particular, we focus on aspects of implementing a full-featured FDTD package that go beyond standard textbook descriptions of the algorithm, or ways in which
Meep differs from typical FDTD implementations. These include pervasive interpolation and accurate
modeling of subpixel features, advanced signal processing, support for nonlinear materials via Padé approximants, and ﬂexible scripting capabilities.
Program summary
Program title: Meep
Catalogue identiﬁer: AEFU_v1_0
Program summary URL:: http://cpc.cs.qub.ac.uk/summaries/AEFU_v1_0.html
Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland
Licensing provisions: GNU GPL
No. of lines in distributed program, including test data, etc.: 151 821
No. of bytes in distributed program, including test data, etc.: 1 925 774
Distribution format: tar.gz
Programming language: C++
Computer: Any computer with a Unix-like system and a C++ compiler; optionally exploits additional free
software packages: GNU Guile [1], libctl interface library [2], HDF5 [3], MPI message-passing interface [4],and Harminv ﬁlter-diagonalization [5]. Developed on 2.8 GHz Intel Core 2 Duo.
Operating system: Any Unix-like system; developed under Debian GNU/Linux 5.0.2.
RAM: Problem dependent (roughly 100 bytes per pixel/voxel)
Classiﬁcation: 10
External routines: Optionally exploits additional free software packages: GNU Guile [1], libctl interface
library [2], HDF5 [3], MPI message-passing interface [4], and Harminv ﬁlter-diagonalization [5] (which
requires LAPACK and BLAS linear-algebra software [6]).
Nature of problem: Classical electrodynamics
Solution method: Finite-difference time-domain (FDTD) method
Running time: Problem dependent (typically about 10 ns per pixel per timestep)
References:[1] GNU Guile, http://www.gnu.org/software/guile[2] Libctl, http://ab-initio.mit.edu/libctl[3] M. Folk, R.E. McGrath, N. Yeager, HDF: An update and future directions, in: Proc. 1999 Geoscience andRemote Sensing Symposium (IGARSS), Hamburg, Germany, vol. 1, IEEE Press, 1999, pp. 273–275.

✩
This paper and its associated computer program are available via the Computer Physics Communications homepage on ScienceDirect (http://www.sciencedirect.com/
science/journal/00104655).
Corresponding author at: Center for Materials Science & Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, United States.
Principal corresponding author.
E-mail addresses: ardavan@mit.edu (A.F. Oskooi), roundyd@physics.oregonstate.edu (D. Roundy), michel@alum.mit.edu (M. Ibanescu), bermel@mit.edu (P. Bermel),joannop@mit.edu (J.D. Joannopoulos), stevenj@math.mit.edu (S.G. Johnson).

***0010-4655/$ – see front matterdoi:10.1016/j.cpc.2009.11.008© 2009 Elsevier B.V. All rights reserved.

688

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702[4] T.M. Forum, MPI: A Message Passing Interface, in: Supercomputing 93, Portland, OR, 1993, pp. 878–
883.
[5] Harminv, http://ab-initio.mit.edu/harminv.
[6] LAPACK, http://www.netlib.org/lapack/lug.

© 2009 Elsevier B.V. All rights reserved.

1. Introduction
One of the most common computational tools in classicalelectromagnetism is the ﬁnite-difference time-domain (FDTD) algorithm, which divides space and time into a regular grid andsimulates the time evolution of Maxwell’s equations [1–5]. This
paper describes our free, open-source implementation of the
FDTD algorithm: Meep (an acronym for MIT Electromagnetic Equation Propagation), available online at http://ab-initio.mit.edu/meep.
Meep is full-featured, including, for example: arbitrary anisotropic,
nonlinear, and dispersive electric and magnetic media; a variety of boundary conditions including symmetries and perfectly
matched layers (PML); distributed-memory parallelism; Cartesian(1d/2d/3d) and cylindrical coordinates; and ﬂexible output andﬁeld computations. It also includes some unusual features, suchas advanced signal processing to analyze resonant modes, accurate
subpixel averaging, a frequency-domain solver that exploits thetime-domain code, complete scriptability, and integrated optimization facilities. Here, rather than review the well-known FDTD algorithm itself (which is thoroughly covered elsewhere), we focus on
the particular design decisions that went into the development of
Meep whose motivation may not be apparent from textbook FDTD
descriptions, the tension between abstraction and performance in
FDTD implementations, and the unique or unusual features of oursoftware.
Why implement yet another FDTD program? Literally dozens
of commercial FDTD software packages are available for purchase,but the needs of research often demand the ﬂexibility provided
by access to the source code (and relaxed licensing constraints tospeed porting to new clusters and supercomputers). Our interactions with other photonics researchers suggest that many groups
end up developing their own FDTD code to serve their needs (ourown groups have used at least three distinct in-house FDTD implementations over the past 15 years), a duplication of effort that
seems wasteful. Most of these are not released to the public, andthe handful of other free-software FDTD programs that could be
downloaded when Meep was ﬁrst released in 2006 were not nearlyfull-featured enough for our purposes. Since then, Meep has been
cited in over 100 journal publications and has been downloaded
over 10,000 times, reaﬃrming the demand for such a package.
FDTD algorithms are, of course, only one of many numerical tools that have been developed in computational electromagnetism, and may perhaps seem primitive in light of other sophisticated techniques, such as ﬁnite-element methods (FEMs) withhigh-order accuracy and/or adaptive unstructured meshes [6–8],
or even radically different approaches such as boundary-element
methods (BEMs) that discretize only interfaces between homogeneous materials rather than volumes [9–12]. Each tool, of course,
has its strengths and weaknesses, and we do not believe that any
single one is a panacea. The nonuniform unstructured grids of
FEMs, for example, have compelling advantages for metallic structures where micrometer wavelengths may be paired with nanometer skin depths. On the other hand, this ﬂexibility comes at a priceof substantial software complexity, which may not be worthwhile
for dielectric devices at infrared wavelengths (such as in integrated
optics or ﬁbers) where the refractive index (and hence the typicalresolution required) varies by less than a factor of four between
materials, while small features such as surface roughness can be

accurately handled by perturbative techniques [13]. BEMs, basedon integral-equation formulations of electromagnetism, are especially powerful for scattering problems involving small objects ina large volume, since the volume need not be discretized and no
artiﬁcial “absorbing boundaries” are needed. On the other hand,
BEMs have a number of limitations: they may still require artiﬁcial
absorbers for interfaces extending to inﬁnity (such as input/output
waveguides) [14]; any change to the Green’s function (such as introduction of anisotropic materials, imposition of periodic or symmetry boundary conditions, or a switch from three to two dimensions) requires re-implementation of large portions of the software(e.g. singular panel integrations and fast solvers) rather than purely
local changes as in FDTD or FEM; continuously varying (as opposed to piecewise-constant) materials are ineﬃcient; and solution
in the time domain (rather than frequency domain, which is inadequate for nonlinear or active systems in which frequency is notconserved) with BEM requires an expensive solver that is nonlocal
in time as well as in space [11]. And then, of course, there are specialized tools that solve only a particular type of electromagnetic
problem, such as our own MPB software that only computes eigenmodes (e.g. waveguide modes) [15], which are powerful and robust
within their domain but are not a substitute for a general-purpose
Maxwell simulation. FDTD has the advantages of simplicity, generality, and robustness: it is straightforward to implement the fulltime-dependent Maxwell equations for nearly arbitrary materials (including nonlinear, anisotropic, dispersive, and time-varying
materials) and a wide variety of boundary conditions, one canquickly experiment with new physics coupled to Maxwell’s equations (such as populations of excited atoms for lasing [16–20]), andthe algorithm is easily parallelized to run on clusters or supercomputers. This simplicity is especially attractive to researchers whose
primary concern is investigating new interactions of physical processes, and for whom programmer time and the training of newstudents is far more expensive than computer time.
The starting point for any FDTD solver is the time-derivative
parts of Maxwell’s equations, which in their simplest form can be
written:∂B= −∇ × E − J B ,∂t∂D= +∇ × H − J,
∂t

(1)
(2)

where (respectively) E and H are the macroscopic electric andmagnetic ﬁelds, D and B are the electric displacement and magnetic induction ﬁelds [21], J is the electric-charge current density, and J B is a ﬁctitious magnetic-charge current density (sometimes convenient in calculations, e.g. for magnetic-dipole sources).
In time-domain calculations, one typically solves the initial-value
problem where the ﬁelds and currents are zero for t < 0, andthen nonzero values evolve in response to some currents J(x, t )and/or J B (x, t ). (In contrast, a frequency-domain solver assumes atime dependence of e −i ωt for all currents and ﬁelds, and solves theresulting linear equations for the steady-state response or eigenmodes [22, Appendix D].) We prefer to use dimensionless units
ε0 = μ0 = c = 1. From our perspective, this choice emphasizes
both the scale invariance of Maxwell’s equations [22, Chapter 2]
and also the fact that the most meaningful quantities to calculate are almost always dimensionless ratios (such as scattered

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702

689

Fig. 1. The computational cell is divided into chunks (left) that have a one-pixel overlap (gray regions). Each chunk (right) represents a portion of the Yee grid, partitioned
into owned points (chunk interior) and not-owned points (gray regions around the chunk edges) that are determined from other chunks and/or via boundary conditions.
Every point in the interior of the computational cell is owned by exactly one chunk, the chunk responsible for timestepping that point.

power over incident power, or wavelength over some characteristic lengthscale). The user can pick any desired unit of distancea (either an SI unit such as a = 1 μm or some typical lengthscale of a given problem), and all distances are given in unitsof a, all times in units of a/c, and all frequencies in units of c /a.
In a linear dispersionless medium, the constituent relations are
D = ε E and B = μH, where ε and μ are the relative permittivity
and permeability (possibly tensors); the case of nonlinear and/or
dispersive media (including conductivities) is discussed further in
Section 4.
The remaining paper is organized as follows. In Section 2, wediscuss the discretization and coordinate system; in addition to
the standard Yee discretization [1], this raises the question of how
exactly the grid is described and divided into “chunks” for parallelization, PML, and other purposes. Section 3 describes a central principle of Meep’s design, pervasive interpolation providing (asmuch as possible) the illusion of continuity in the speciﬁcation ofsources, materials, outputs, and so on. This led to the development of several techniques unique to Meep, such as a scheme for
subpixel material averaging designed to eliminate the ﬁrst-ordererror usually associated with averaging techniques or stairstepping of interfaces. In Section 4, we describe and motivate ourtechniques for implementing nonlinear and dispersive materials,
including a slightly unusual method to implement nonlinear materials using a Padé approximant that eliminates the need to solve
cubic equations for every pixel. Section 5 describes how typicalcomputations are performed in Meep, such as memory-eﬃcienttransmission spectra or sophisticated analysis of resonant modes
via harmonic inversion. This section also describes how we haveadapted the time-domain code, almost without modiﬁcation, tosolve frequency-domain problems with much faster convergenceto the steady-state response than merely time-stepping. The userinterface of Meep is discussed in Section 6, explaining the considerations that led us to a scripting interface (rather than a GUI or
CAD interface). Section 7 describes some of the tradeoffs between
performance and generality in this type of code and the speciﬁccompromises chosen in Meep. Finally, we make some concluding
remarks in Section 8.
2. Grids and boundary conditions
The starting point for the FDTD algorithm is the discretizationof space and time into a grid. In particular, Meep uses the standard
Yee grid discretization which staggers the electric and magnetic

ﬁelds in time and in space, with each ﬁeld component sampled at
different spatial locations offset by half a pixel, allowing the timeand space derivatives to be formulated as center-difference approximations [23]. This much is common to nearly every FDTD implementation and is described in detail elsewhere [1]. In order toparallelize Meep, eﬃciently support simulations with symmetries,
and to eﬃciently store auxiliary ﬁelds only in certain regions (forPML absorbing layers), Meep further divides the grid into chunks
that are joined together into an arbitrary topology via boundary
conditions. (In the future, different chunks may have different resolutions to implement a nonuniform grid [24–27].) Furthermore,
we distinguish two coordinate systems: one consisting of integer
coordinates on the Yee grid, and one of continuous coordinates in“physical” space that are interpolated as necessary onto the grid(see Section 3). This section describes those concepts as they are
implemented in Meep, as they form a foundation for the remaining sections and the overall design of the Meep software.
2.1. Coordinates and grids
The two spatial coordinate systems in Meep are described by
the vec, a continuous vector in Rd (in d dimensions), and theivec, an integer-valued vector in Zd describing locations on the
Yee grid. If n is an ivec, the corresponding vec is given by
0.5xn, where x is the spatial resolution (the same along x, y,
and z)—that is, the integer coordinates in an ivec correspondto half -pixels, as shown in the right panel of Fig. 1. This is torepresent locations on the spatial Yee grid, which offsets different ﬁeld components in space by half a pixel as shown (in 2d)in the right panel of Fig. 1. In 3d, the E x and D x components
are sampled at ivecs (2 + 1, 2m, 2n), E y and D y are sampled
at ivecs (2, 2m + 1, 2n), and so on; H x and B x are sampled
at ivecs (2, 2m + 1, 2n + 1), H y and B y are sampled at ivecs(2 + 1, 2m, 2n + 1), and so on. In addition to these grids for thedifferent ﬁeld components, we also occasionally refer to the centered grid, at odd ivecs (2 + 1, 2m + 1, 2n + 1) corresponding to
the “center” of each pixel. (The origin of the coordinate systems is
an arbitrary ivec that can be set by the user, but is typically thecenter of the computational volume.) The philosophy of Meep, as
described in Section 3, is that as much as possible the user should
be concerned only with continuous physical coordinates (vecs),
and the interpolation/discretization onto ivecs occurs internallyas transparently as possible.

690

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702

Fig. 2. Meep can exploit mirror and rotational symmetries, such as the 180-degree (C 2 ) rotational symmetry of the S-shaped structure in this schematic example. Although
Meep maintains the illusion that the entire structure is stored and simulated, internally only half of the structure is stored (as shown at right), and the other half is inferredby rotation. The rotation gives a boundary condition for the not-owned grid points along the dashed line.

2.2. Grid chunks and owned points
An FDTD simulation must occur within a ﬁnite volume of space,
the computational cell, terminated with some boundary conditions
and possibly by absorbing PML regions as described below. This
(rectilinear) computational cell, however, is further subdivided into
convex rectilinear chunks. On a parallel computer, for example, different chunks may be stored at different processors. In order tosimplify the calculations for each chunk, we employ the common
technique of padding each chunk with extra “boundary” pixels that
store the boundary values [28] (shown as gray regions in Fig. 1)—
this means that the chunks are overlapping in the interior of thecomputational cell, where the overlaps require communication tosynchronize the values.
More precisely, the grid points in each chunk are partitioned
into owned and not-owned points. The not-owned points are determined by communication with other chunks and/or by boundary
conditions. The owned points are time-stepped within the chunk,
independently of the other chunks (and possibly in parallel), andevery grid point inside the computational cell is owned by exactly one
chunk.
The question then arises: how do we decide which points
within the chunk are owned? In order for a grid point to beowned, the chunk must contain all the information necessary fortimestepping that point (once the not-owned points have beencommunicated). For example, for a D y point (2, 2m + 1, 2n) tobe owned, the H z points at (2 ± 1, 2m + 1, 2n) must both be in
the chunk in order to compute ∇ × H for timestepping D at thatpoint. This means that the D y points along the left (minimum-x)
edge of the chunk (as shown in the right panel of Fig. 1) cannot
be owned: there is no H z point to the left of it. An additional dependency is imposed by the case of anisotropic media: if there isan εxy coupling E x to D y , then updating E x at (2 + 1, 2m, 2n) requires the four D y values at (2 + 1 ± 1, 2m ± 1, 2n) (these are thesurrounding D y values, as seen in the right panel of Fig. 1). Thismeans that the E x (and D x ) points along the right (maximum-x)
edge of the chunk (as shown in the right panel of Fig. 1) cannot
be owned either: there is no D y point to the right of it. Similarlyfor ∇ × D and anisotropic μ.
All of these considerations result in the shaded-gray region ofFig. 1 (right) being not-owned. That is, if the chunk intersects k + 1
pixels along a given direction starting at an ivec coordinate of 0(e.g. k = 5 in Fig. 1), the endpoint ivec coordinates 0 and 2k + 1
are not-owned and the interior coordinates from 1 to 2k (inclusive)
are owned.
2.3. Boundary conditions and symmetries
All of the not-owned points in a chunk must be determinedby boundary conditions of some sort. The simplest boundary conditions are when the not-owned points are owned by some otherchunk, in which case the values are simply copied from that chunk(possibly requiring communication on a multiprocessor system)
each time they are updated. In order to minimize communications
overhead, all communications between two chunks are batched
into a single message (by copying the relevant not-owned points
to/from a contiguous buffer) rather than sending one message perpoint to be copied.
At the edges of the computational cell, some user-selectedboundary condition must be imposed. For example, one canuse perfect electric or magnetic conductors where the relevant
electric/magnetic-ﬁeld components are set to zero at the boundaries. One can also use Bloch-periodic boundary conditions, wherethe ﬁelds on one side of the computational cell are copied fromthe other side of the computational cell, optionally multiplied by
a complex phase factor e iki Λi where ki is the propagation constantin the ith direction, and Λi is the length of the computational cellin the same direction. Meep does not implement any absorbing
boundary conditions—absorbing boundaries are, instead, handled
by an artiﬁcial material, perfectly matched layers (PML), placed
adjacent to the boundaries [1].
Bloch-periodic boundary conditions are useful in periodic systems [22], but this is only one example of a useful symmetry that
may be exploited via boundary conditions. One may also have mirror and rotational symmetries. For example, if the materials andthe ﬁeld sources have a mirror symmetry, one can cut the computational costs in two by storing chunks only in half the computational cell and applying mirror boundary conditions to obtainthe not-owned pixels adjacent to the mirror plane. As a more unusual example, consider an S-shaped structure as in Fig. 2, which
has no mirror symmetry but is symmetric under 180-degree rotation, called C 2 symmetry [29]. Meep can exploit this case as well(assuming the current sources have the same symmetry), storing
only half of the computational cell as in Fig. 2 and inferring the
not-owned values along the dashed line by a 180-degree rotation.
(In the simple case where the stored region is a single chunk, thismeans that the not-owned points are determined by owned points
in the same chunk, requiring copies, possibly with sign ﬂips.) Depending on the sources, of course, the ﬁelds can be even or odd
under mirror ﬂips or C 2 rotations [22], so the user can specify anadditional sign ﬂip for the transformation of the vector ﬁelds (andpseudovector H and B ﬁelds, which incur an additional sign ﬂip
under mirror reﬂections [21,22]). Meep also supports four-fold rotation symmetry (C 4 ), where the ﬁeld can be multiplied by factorsof 1, i, −1, or −i under each 90-degree rotation [29]. (Other rotations, such as three-fold or six-fold, are not supported because
they do not preserve the Cartesian Yee grid.) In 2d, the xy-plane
is itself a mirror plane (unless in the presence of anisotropic materials) and the symmetry decouples TE modes (with ﬁelds E x , E y ,and H z ) from TM modes (H x , H y , and E z ) [22]; in this case Meep
only allocates those ﬁelds for which the corresponding sources arepresent.
A central principle of Meep is that symmetry optimizations betransparent to the user once the desired symmetries are speci-

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702

691

Fig. 3. A key principle of Meep is that continuously varying inputs yield continuously varying outputs. Here, an eigenfrequency of a photonic crystal varies continuously withthe eccentricity of a dielectric rod, accomplished by subpixel smoothing of the material parameters, whereas the nonsmoothed result is “stairstepped”. Speciﬁcally, the plotshows a TE eigenfrequency of 2d square lattice (period a) of dielectric ellipses (ε = 12) in air versus one semi-axis diameter of the ellipse (in gradations of 0.005a) for nosmoothing (red squares, resolution of 20 pixels/a), subpixel smoothing (blue circles, resolution of 20 pixels/a) and “exact” results (black line, no smoothing at resolution of
200 pixels/a). (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)ﬁed. Meep maintains the illusion that the entire computational
cell is computed—for example, the ﬁelds in the entire computational cell can still be queried or exported to a ﬁle, ﬂux planes andsimilar computations can still extend anywhere within the computational cell, and so on. The ﬁelds in the nonstored regions are
simply computed behind the scenes (without ever allocating memory for them) by transforming the stored chunks as needed. A key
enabling factor for maintaining this illusion eﬃciently is the loopin-chunks abstraction employed by the Meep code, described in
Section 7.
Meep also supports continuous rotational symmetry around a
given axis, where the structure is invariant under rotations and theﬁelds transform as e imφ for some m [22], but this is implemented
separately by providing the option to simulate Maxwell’s equationsin the (r , z)-plane with cylindrical coordinates, for which operators
like ∇ × change form.
3. Interpolation and the illusion of continuity
A core design philosophy of Meep is to provide the illusion ofcontinuous space and time, masking the underlying discretizationfrom the user as much as possible. There are two components tothis approach: the input and the outputs. Continuously varying inputs, such as the geometry, materials, and the source currents, leadto continuously varying outputs, as in the example of Fig. 3. Similarly, the value of any ﬁeld (or any function of the ﬁelds) can be
output at any point in space or integrated over any region. Furthermore, the effects of these inputs and the resulting outputs must
converge as quickly as possible to the exact solution as the resolution increases. In this section, we discuss how this illusion ofcontinuity is implemented for ﬁeld outputs, current inputs, and geometry/materials.
Any ﬁeld component (or any combinations such as ﬂux, energy,
and user-deﬁned functions) can be evaluated at any point in space.
In general, this requires interpolation from the Yee grid. Since theunderlying FDTD center-difference algorithm has second-order accuracy, we linearly interpolate ﬁelds as needed (which also has
second-order accuracy for smooth functions). Similarly, we provide an interface to integrate any function of the ﬁelds over any
convex rectilinear region (boxes, planes, or lines), and the integral

is computed by integrating the linear interpolation of the ﬁelds
within the integration region. This is straightforward, but there are
two subtleties due to the staggered Yee grid. First, computationof quantities like E × H that mix different ﬁeld components requires an additional interpolation: ﬁrst, the ﬁelds are interpolated
onto the centered grid (Section 2), then the integrand is computed,
and then the linear interpolation of the integrand is integrated
over the speciﬁed region. Second, the computation of quantities
like E × H mixes two ﬁelds that are stored at different times: H
is stored at times (n − 0.5)t, while E is stored at times nt [1].
Simply using these time-offset ﬁelds together is only ﬁrst-orderaccurate. If second-order accuracy is desired, Meep provides theoption to temporarily synchronize the electric and magnetic ﬁelds:
the magnetic ﬁelds are saved to a backup array, stepped by t, andthey are averaged with the backup array to obtain the magnetic
ﬁelds at nt with O (t 2 ) accuracy. (The ﬁelds are restored frombackup before resuming timestepping.) This restores second-orderaccuracy at the expense of an extra half a timestep’s computation, which is usually negligible because such ﬁeld computations
are rarely required at every timestep of a simulation—see Section 5for how Meep performs typical transmission simulations and othercalculations eﬃciently.
The conceptually reversed process is required for specifying
sources: the current density is speciﬁed at some point (for dipole
sources) or in some region (for distributed current sources) incontinuous space, and then must be restricted to a corresponding
current source on the Yee grid. Meep performs this restriction using exactly the same code (the loop-in-chunks abstraction of Section 7) and the same weights as the interpolation procedure above.
Mathematically, we are exploiting a well-known concept (originating in multigrid methods) that restriction can be deﬁned as thetranspose of interpolation [30]. This is illustrated by a 2d example in Fig. 4. Suppose that the bilinear interpolation f (blue) offour grid points (red) is f = 0.32 f 1 + 0.48 f 2 + 0.08 f 3 + 0.12 f 4 ,which can be viewed as multiplying a vector of those ﬁelds by therow-vector [0.32, 0.48, 0.08, 0.12]. Conversely, if we place a pointdipole current source J (blue) at the same point, it is restricted on
the grid (red) to values J 1 = 0.32 J , J 2 = 0.48 J , J 3 = 0.08 J , andJ 4 = 0.12 J as shown in Fig. 4, corresponding to multiplying J by

692

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702

Fig. 4. Left: a bilinear interpolation of values f 1,2,3,4 on the grid (red) to the value f at an arbitrary point. Right: the reverse process is restriction, taking a value J at anarbitrary point (e.g. a current source) and converting into values on the grid. Restriction can be viewed as the transpose of interpolation and uses the same coeﬃcients. (Forinterpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)

Fig. 5. Cerenkov radiation emitted by a point charge moving at a speed v = 1.05c /n exceeding the phase velocity of light in a homogeneous medium of index n = 1.5. Thanks
to Meep’s interpolation (or technically restriction), the smooth motion of the source current (left panel) can be expressed as continuously varying currents on the grid,whereas the nonsmooth pixelized motion (no interpolation) (right panel) reveals high-frequency numerical artifacts of the discretization (counter-propagating wavefronts
behind the moving charge).

the column vector [0.32, 0.48, 0.08, 0.12]T .1 Such a restriction has
the property of preserving the sum (integral) of the currents, andtypically leads to second-order convergence of the resulting ﬁelds
as the resolution increases (see below). An example of the utility ofthis continuous restriction process is shown in Fig. 5 via the phenomenon of Cerenkov radiation [31]: a point charge q moving at aconstant velocity v with a magnitude 1.05c /n exceeding the phasevelocity c /n in the medium emits a shockwave-like radiation pattern, and this can be directly modeled in Meep by a continuously
moving current source J = −vqδ(x − vt ) [32]. In contrast, pixelizing the motion into discrete jumps to the nearest grid point leads
to visible numerical artifacts in the radiation, as seen in the rightpanel of Fig. 5.
All of the second-order accuracy of FDTD and the above interpolations is generally spoiled to only ﬁrst-order, however, if
one directly discretizes a discontinuous material boundary [33,35].
Moreover, directly discretizing a discontinuity in ε or μ leads to“stairstepped” interfaces that can only be varied in discrete jumps
of one pixel at a time. Both of these problems are solved in
Meep by using an appropriate subpixel smoothing of ε and μ:
before discretizing, discontinuities are smoothed into continuoustransitions over a distance of one pixel x, using a carefully designed averaging procedure. Any subpixel smoothing techniquewill achieve the goal of continuously varying results as the geometry is continuously varied. In the case of Meep this is illustrated by

1

Technically, for a dipole-current source given by a delta function with amplitude

I , the corresponding current density is J = I /xd in d dimensions.

Fig. 3: in a 2d photonic crystal (square lattice of dielectric rods),
the lowest TE-polarization eigenfrequency (computed as in Section 5) varies continuously with the eccentricity of the ellipticalrods for subpixel averaging, whereas the nonaveraged discontinuous discretization produces a stairstepped discontinuous eigenfrequency. On the other hand, most subpixel smoothing techniques
will not increase the accuracy of FDTD—on the contrary, smoothing
discontinuous interfaces changes the structure, and generally introduces additional error into the simulation [33]. In order to design
an accurate smoothing technique, we exploited recent results inperturbation theory that show how a particular subpixel smoothing can be chosen to yield zero ﬁrst-order error [13,33,34,36]. The
results are shown in Figs. 6 and 7: for both computation of theeigenfrequencies (of an anisotropic photonic crystal) in Fig. 6 andthe scattering loss from a bump on a strip waveguide in Fig. 7, theerrors in Meep’s results decrease quadratically [O (x2 )], whereas
doing no averaging leads to erratic linear convergence [O (x)].
Furthermore, Fig. 6 compares to other subpixel-averaging schemes,
including the obvious strategy of simply averaging ε within eachpixel [37], and shows that they lead to ﬁrst-order convergence nobetter than no averaging at all.
The subpixel averaging is discussed in more detail elsewhere[33,34,36], so we only brieﬂy summarize it here. In order for thesmoothing to yield zero ﬁrst-order perturbation, the smoothing
scheme must be anisotropic. Even if the initial interface is between
isotropic materials, one obtains a tensor ε (or μ) which uses themean ε for ﬁelds parallel to the interface and the harmonic mean
(inverse of mean of ε −1 ) for ﬁelds perpendicular to the interface—
this was initially proposed heuristically [38] and later shown to be

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702

693

Fig. 6. Appropriate subpixel averaging can increase the accuracy of FDTD with discontinuous materials [33,34]. Here, relative error ω/ω (comparing to the “exact” ω0 froma planewave calculation [15]) for an eigenmode calculation (as in Section 5) for a cubic lattice (period a) of 3d anisotropic-ε ellipsoids (right inset) versus spatial resolution(units of pixels per vacuum wavelength λ), for a variety of subpixel smoothing techniques. Straight lines for perfect linear (black dashed) and perfect quadratic (black solid)
convergence are shown for reference. Most curves are for the ﬁrst eigenvalue band (left inset shows E x in xy cross-section of unit cell), with vacuum wavelength λ = 5.15a.
Hollow squares show Meep’s method for band 13 (middle inset), with λ = 2.52a. Meep’s method for bands 1 and 13 is shown for resolutions up to 100 pixels/a.

Fig. 7. The relative error in the scattered power from a small semicircular bump in a dielectric waveguide (ε = 12), excited by a point-dipole source in the waveguide(geometry and ﬁelds shown in inset), as a function of the computational resolution. Appropriate subpixel smoothing of the dielectric interfaces leads to roughly second-order[O (x2 )] convergence (red squares), whereas the unsmoothed structure has only ﬁrst-order convergence (blue circles). (For interpretation of the references to color in thisﬁgure legend, the reader is referred to the web version of this article.)justiﬁed via perturbation theory [13,33]. (If the initial materials areanisotropic, a more complicated formula is needed [34,36].) The
key point is that, even if the physical structure consists entirelyof isotropic materials, the discretized structure will use anisotropic
materials. Stable simulation of anisotropic media requires an FDTD
variant recently proposed in Ref. [39].
There are a few limitations to this subpixel averaging. First, thecase of perfect metals requires a different approach [40,41] thatis not yet implemented in Meep. Although Meep does not yet

implement subpixel averaging for dispersive materials, there is numerical evidence that similar accuracy improvements are obtainedin that case by the same technique [42], and we suspect that asimilar derivation can be applied (using the unconjugated form of
perturbation theory for the complex-symmetric Maxwell equations
in reciprocal media with losses [43]). Second, once the smoothingeliminates the ﬁrst-order error, the presence of sharp corners (associated with ﬁeld singularities) introduce an error intermediate
between ﬁrst- and second-order [33], which we hope to address in694

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702

 1 + ( χ (2) D̃ ) + 2( χ (3) D̃2 ) i (1) −1[ε (1) ]3εD̃ i ,
χ (2)
χ (3)
2
1 + 2( (1) 2 D̃ i ) + 3( (1) 3 D̃ )
[ε ]
[ε ]

future work. Third, the ﬁelds directly on the interface are still atbest ﬁrst-order accurate even with subpixel smoothing—however,these localized errors are equivalent to currents that radiate zero
power to ﬁrst order [36,44]. The improved accuracy from smoothing is therefore obtained for ﬁelds evaluated off of the interface
as in scattered ﬂux integrated over a surface away from the interface (Fig. 7), for nonlocal properties like resonant frequencies
and eigenfrequencies (Fig. 6), and for overall integrals of ﬁelds andenergies [to which the interface contributes only O (x) of the integration domain and hence ﬁrst-order errors on the interface havea second-order effect].

where D̃ i = D i − P i . For the case of isotropic ε (1) and χ (2) = 0,
so that we have a purely Kerr (χ (3) ) material, this matches the“exact” cubic E to O ( D 7 ) error. With χ (2) = 0, the error is O ( D 4 ).
For more complicated dispersive nonlinear media or for arbitrary anisotropy in χ (2) or χ (3) , one approach that Meep may
implement in the future is to incorporate the nonlinear terms in
the auxiliary differential equations for a Lorentz medium [1].

4. Materials

4.2. Absorbing boundary layers: PML, pseudo-PML, and quasi-PML

Time-dependent methods for electromagnetism, given theirgenerality, allow for the simulation of a broad range of material systems. Certain classes of materials, particularly active andnonlinear materials which do not conserve frequency, are ideallysuited for modeling by such methods. Materials are represented
in Maxwell’s equations (1) and (2) via the relative permittivity
ε(x) and permeability μ(x) which in general depend on position,
frequency (material dispersion) and the ﬁelds themselves (nonlinearities). Meep currently supports arbitrary anisotropic material
tensors, anisotropic dispersive materials (Lorentz–Drude models
and conductivities, both magnetic and electric), and nonlinear materials (both second- and third-order nonlinearities), which taken
together permit investigations of a wide range of physical phenomena. The implementation of these materials in Meep is mostly
based on standard techniques [1], so we will focus here on twoplaces where Meep differs from the usual approach. For nonlinearities, we use a Padé approximant to avoid solving cubic equations
at each step. For PML absorbing media in cylindrical coordinates,
we only use a “quasi-PML” [46] based on a Cartesian PML, but
explain why its performance is comparable to a true PML while
requiring less computational effort.

A perfectly matched layer (PML) is an artiﬁcial absorbing
medium that is commonly used to truncate computational grids
for simulating wave equations (e.g. Maxwell’s equations), and isdesigned to have the property that interfaces between the PML
and adjacent media are reﬂectionless in the exact wave equation [1]. There are various interchangeable formulations of PML
for FDTD methods [1], which are all equivalent to a coordinate
stretching of Maxwell’s equations into complex spatial coordinates;
Meep implements a version of the uniaxial PML (UPML), expressing the PML as an effective dispersive anisotropic ε and μ [1].
Meep provides support for arbitrary user-speciﬁed PML absorption
proﬁles (which have an important inﬂuence on reﬂections due todiscretization error and other effects) for a given round-trip reﬂection (describing the strength of the PML in terms of the amplitudeof light passing through the PML, reﬂecting off the edge of thecomputational cell, and propagating back) [45]. For the case of periodic media such as photonic crystals, the medium is not analytic
and the premise of PML’s reﬂectionless property is violated; in thiscase, a “PML” material overlapped with the photonic crystal is onlya “pseudo-PML” that is reﬂectionless only in the limit of a suﬃciently thick and gradual absorber, and control over the absorption
proﬁle is important [45].
For the radial direction in cylindrical coordinates, a true PML
can be derived by coordinate-stretching, but it requires more storage and computational effort than the Cartesian UPML [52,53],
as well as increasing code complexity. Instead, we chose to implement a quasi-PML [46], which simply consists of using theCartesian UPML materials as an approximation for the true radial PML. This approximation becomes more and more accurate
as the outer radius of the computational cell increases, because
the implicit curvature of the PML region decreases with radius andapproaches the Cartesian case. Furthermore, one must recall thatevery PML has reﬂections once space is discretized [1], which canbe mitigated by gradually turning on the PML absorption over aﬁnite-thickness PML layer. The quasi-PML approximation is likewise mitigated by the same gradual absorption proﬁle, and theonly question is that of the constant factor in the reﬂection convergence: how thick does the quasi-PML need to be to achieve low
reﬂections, compared to a true PML? Fig. 8 shows that, for a typical calculation, the performance of the quasi-PML in cylindricalcoordinates (left) is comparable to that of a true PML in Cartesian
coordinates (right). Here, we plot a measure of the reﬂection fromthe PML as a function of the PML absorber length L, for a ﬁxed
round-trip reﬂection [45], using as a measure of the reﬂection the“ﬁeld convergence” factor: the difference between the E ﬁeld at agiven point for simulations with PML absorber lengths L and L + 1.
The PML conductivity σ (x) is turned on gradually as (x/ L )d for
d = 1, 2, 3, and it can be shown that this leads to reﬂections that
decrease as 1/ L 2d+2 and ﬁeld-convergence factors that decrease as
1/ L 2d+4 [45]. Precisely these decay rates are observed in Fig. 8,
with similar constant coeﬃcients. As the resolution is increased(approaching the exact wave equations), the constant coeﬃcient in

4.1. Nonlinear materials
Optical nonlinearities arise when large ﬁeld intensities inducechanges in the local ε or μ to produce a number of interesting effects: temporal and spatial soliton propagation, optical bistability,
self-focusing of optical beams, second- and third-harmonic generation, and many other effects [47,48]. Such materials are usuallydescribed by a power-series expansion of D in terms of E andvarious susceptibilities. In many common materials, or when considering phenomena in a suﬃciently narrow bandwidth (such as
the resonantly enhanced nonlinear effects [49] well-suited to FDTD
calculations), these nonlinear susceptibilities can be accurately approximated via nondispersive (instantaneous) effects [50]. Meep
supports instantaneous isotropic (or diagonal anisotropic) nonlinear effects of the form:(2 ) 2D i − P i = ε (1 ) E i + χ i(3 )E i + χi

| E |2 E i ,(3)where ε (1) represents all the linear nondispersive terms and P i is a(1)
dispersive polarization P = χdispersive (ω)E from dispersive materials such as Lorentz media [1]. (A similar equation relates B and H.)
Implementing this equation directly, however, would require one
to solve a cubic equation at each time step [1, Section 9.6], since
D is updated from ∇ × H before updating E from D.
However, Eq. (3) is merely a power-series approximation for thetrue material response, valid for suﬃciently small ﬁeld intensities,
so it is not necessary to insist that it be solved exactly. Instead, we
approximate the solution of Eq. (3) by a Padé approximant [51],which matches the “exact” cubic solution to high-order accuracy
by the rational function:Ei =[ε (1) ]2(4)

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702

695

Fig. 8. The performance of a quasi-PML in the radial direction (cylindrical coordinates, left panel) at a resolution 20 pixels/λ is nearly equivalent to that of a true PML (inCartesian coordinates, right panel). The plot shows the difference in the electric ﬁeld E z (insets) from a point source between simulations with PML thickness L and L + 1,
which is a simple proxy for the PML reﬂections [45]. The different curves are for PML conductivities that turn on as (x/ L )d for d = 1, 2, 3 in the PML, leading to different
rates of convergence of the reﬂection [45].

the Cartesian PML plot will decrease (approaching zero reﬂection),while the quasi-PML’s constant coeﬃcient will saturate at someminimum (corresponding to its ﬁnite reﬂectivity in the exact wave
equation for a ﬁxed L). This difference seems of little practical concern, however, because the reﬂection from a one-wavelength thickquasi-PML at a moderate resolution (20 pixels/λ) is already so low.

wide bandwidth including all frequencies of interest), and computeEω and Hω from the Fourier transforms of E(t ) and H(t ). There
are several different ways to compute these Fourier transforms. For
example, one could store the electric and magnetic ﬁelds throughout S over all times and at the end of the simulation perform a
discrete-time Fourier transform (DTFT) of the ﬁelds:

5. Enabling typical computations

Eω =Simulating Maxwell’s equations in the time domain enables
the investigation of problems inherently involving multiple frequencies, such as nonlinearities and active media. However, it isalso well adapted to solving frequency domain problems since it
can solve large bandwidths at once, for example analyzing resonant modes or computing transmission/reﬂection spectra. In thissection, we describe techniques Meep uses to eﬃciently compute scattering spectra and resonant modes in the time domain.
Furthermore, we describe how the time domain method can beadapted to a purely frequency domain solver while sharing almost
all of the underlying code.
5.1. Computing ﬂux spectra
A principle task of computational time-domain tools are investigations of transmission or scattering spectra from arbitrary
structures, where one wants to compute the transmitted or scattered power in a particular direction as a function of the frequency
of incident light. One can solve for the power at many frequencies in a single time-domain simulation by Fourier transforming
the response to a short pulse. Speciﬁcally, for a given surface S,
one wishes to compute the integral of the Poynting ﬂux:

P (ω) = {Eω (x)∗ × Hω (x) dA,(5)e i ωnt E(nt )t ,(6)n

for all frequencies (ω ) of interest, possibly exploiting a fast Fouriertransform (FFT) algorithm. Such an approach has the following
computational cost: for a simulation having T timesteps, F  T
frequencies to compute, N S ﬁelds in the ﬂux region and N pixels in
the entire computational cell this approach requires Θ( N + N S T )storage and Θ( N T + T log T ) time (using a FFT-based chirp-z algorithm [54]).2 The diﬃculty with this approach is that if a longsimulation (large T ) is required to obtain a high frequency resolution by the usual uncertainty relation [56], then the Θ( N S T )storage requirements for the ﬁelds E(t ) and H(t ) at each point in
S become excessive. Instead, Meep accumulates the DTFT summation of the ﬁelds at every point in S as the simulation progresses;
once the time stepping has terminated, Eq. (5) can be evaluatedusing these Fourier-transformed ﬁelds.3 The computational cost of
this approach is Θ( N + N S F ) storage [much less than Θ( N S T ) if
F  T ] and Θ( N T + N S F T ) time. Although our current approachworks well, another possible approach that we have been considering is to use Padé approximation: one stores the ﬁelds at every
timestep on S, but instead of using the DTFT one constructs a Padé
approximant to extrapolate the inﬁnite-time DTFT from a shorttime series [57]. This requires Θ( N + N S T ) storage (but T is potentially much smaller) and O ( N T + T log2 T ) time [58].

S2where Eω and Hω are the ﬁelds produced by a source at frequency ω , and  denotes the real part of the expression. The basic
idea, in time-domain, is to use a short-pulse source (covering a

Here, Θ has the usual meaning of an asymptotic tight bound [55].
It is tempting to instead accumulate the Fourier transform of the Poynting ﬂux
at each time, but this is not correct since the ﬂux is not a linear function of theﬁelds.
3

696

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702

Fig. 9. Relative error in the quality factor Q for a photonic-crystal resonant cavity (inset, period a) with Q ∼ 106 , versus simulation time in units of optical periods of theresonance. Blue circles: ﬁlter-diagonalization method. Red squares: least-squares ﬁt of energy in cavity to a decaying exponential. Filter-diagonalization requires many fewer
optical periods than the decay time Q , whereas curve ﬁtting requires a simulation long enough for the ﬁelds to decay signiﬁcantly. (For interpretation of the references tocolor in this ﬁgure legend, the reader is referred to the web version of this article.)

5.2. Analyzing resonant modes
Another major goal of time-domain simulations is analysis ofresonant phenomena, speciﬁcally by determining the resonant frequency ω0 and the quality factors Q (i.e., the number of optical
cycles 2π /ω0 for the ﬁeld to decay by e −2π ) of one or more resonant modes. One straightforward and common approach to compute ω0 and Q is by computing the DTFT of the ﬁeld at somepoint in the cavity in response to a short pulse [1]: ω0 is thenthe center of a peak in the DTFT and 1/ Q is the fractional widthof the peak at half maximum. The problem with this approach isthat the Fourier uncertainty relation (equivalently, spectral leakage from the ﬁnite time window [56]) means that resolving thepeak in this way requires a simulation much longer than Q /ω0(problematic for structures that may have very high Q , even 109
or higher [59]). Alternatively, one can perform a least squares ﬁt of
the ﬁeld time-series within the cavity to an exponentially decaying
sinusoid, but this leads to an ill-conditioned, nonconvex, nonlinear
ﬁtting problem (and is especially diﬃcult if more than one resonant mode may be present). If only a single resonant mode ispresent, one can perform a least-squares ﬁt of the energy in thecavity to a decaying exponential in order to determine Q , but along simulation is still required to accurately resolve a large Q(as shown below). A more accurate and eﬃcient approach, requiring only a short simulation even for very large Q values, isthe technique of ﬁlter diagonalization originally developed for NMR
spectroscopy, which transforms the time-series data into a small
eigenproblem that is solved for all resonant frequencies and quality factors at once (even for multiple overlapping resonances) [60].
Chapter 16 of Ref. [1] compared the DFT peak-ﬁnding method withﬁlter-diagonalization by attempting to resolve two near-degeneratemodes in a microcavity, and demonstrated the latter’s ability toaccurately resolve closely-spaced peaks with as much as a factor
of ﬁve times fewer timesteps. In our own work, we have used
ﬁlter diagonalization to compute quality factors of 108 or more using simulations only a few hundred optical cycles in length [59].
We quantify the ability of ﬁlter diagonalization to resolve a large

Q ∼ 106 in Fig. 9, comparing the relative error in Q versus simulation time for ﬁlter diagonalization and the least-squares energy-ﬁtmethod above. (The speciﬁc cavity is formed by a missing rod in a
two-dimensional photonic crystal consisting of a square lattice ofdielectric rods in air with period a, radius 0.2a, and ε = 12 [22].)
Fig. 9 demonstrates that ﬁlter diagonalization is able to identify
the quality factor using almost an order of magnitude fewer time
steps than the curve ﬁtting method. (Another possible technique toidentify resonant modes uses Padé approximants, which can alsoachieve high accuracy from a short simulation [57,61].)5.3. Frequency-domain solver
A common electromagnetic problem is to ﬁnd the ﬁelds thatare produced in a geometry in response to a source at a single frequency ω . In principle, the solution of such problems need
not involve time at all, but involve solving a linear equation
purely in the frequency domain [22, Appendix D]; this can beachieved by many methods, such as ﬁnite-element methods [6–8],
boundary-element methods [9–12], or ﬁnite-difference frequencydomain methods [62]. However, if one already has a full-featured
parallel FDTD solver, it is attractive to exploit that solver for
frequency-domain problems when they arise. The most straightforward approach is to simply run a simulation with a constantfrequency source—after a long time, when all transient effects fromthe source turn-on have disappeared, the result is the desired
frequency-domain response. The diﬃculty with this approach isthat a very long simulation may be required, especially if longlived resonant modes are present at nearby frequencies (in whichcase a time
Q /ω is required to reach steady state). Instead,
we show how the FDTD timestep can be used to directly plug afrequency-domain problem into an iterative linear solver, ﬁnding
the frequency-domain response in the equivalent of many fewer
timesteps while exploiting the FDTD code almost without modiﬁcation.
The central component of any FDTD algorithm is the time step:an operation that advances the ﬁeld by t in time. In order toextract a frequency-domain problem from this operation, we ﬁrst

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702

697

Fig. 10. Root-mean-square error in ﬁelds in response to a constant-frequency point source in vacuum (inset), for frequency-domain solver (red squares, adapted from Meep
time-stepping code) vs. time-domain method (blue circles, running until transients decay away). (For interpretation of the references to color in this ﬁgure legend, the reader
is referred to the web version of this article.)

express the timestep as an abstract linear operation: if f n represents all of the ﬁelds (electric and magnetic) at time step n, then(in a linear time-invariant structure) the time step operation canbe expressed in the form:

f n+1 = T̂ 0 f n + sn ,(7)nwhere T̂ 0 is the timestep operator with no sources and s are thesource terms (currents) from that time step. Now, suppose that
one has a time-harmonic source sn = e −i ωnt s and wish to solve
for the resulting time-harmonic (steady state) ﬁelds f n = e −i ωnt f.
Substituting these into Eq. (7), we obtain the following linear equation for the ﬁeld amplitudes f:T̂ 0 − e −i ωt f = −s.

(8)This can then be solved by an iterative method, and the key fact
is that iterative methods for Ax = b only require one to supply
a function that multiplies the linear operator A by a vector [63].
Here, A is represented by T̂ 0 − e −i ωt and hence one can simply
use a standard iterative method by calling the unmodiﬁed timestep
function from FDTD to provide the linear operator. To obtain theproper right-hand side s, one merely needs to execute a single
timestep (7), with sources, starting from zero ﬁeld f = 0. Sincein general this linear operator is not Hermitian (especially in thepresence of PML absorbing regions), we employ the BiCGSTAB-L
algorithm (a generalization of the stabilized biconjugate gradientalgorithm, where increasing the integer parameter L trades off increased storage for faster convergence) [64,65].
This technique means that all of the features implemented in
our time-domain solver (not only arbitrary materials, subpixel averaging, and other physical features, but also parallelization, visualization, and user-interface features) are immediately available
as a frequency-domain solver. To demonstrate the performance of
this frequency-domain solver over the straightforward approachof simply running a long simulation until transients have disappeared, we computed the root-mean-square error in the ﬁeld asa function of the number of time steps (or evaluations of T̂ 0
by BiCGSTAB-L) for two typical simulations. The ﬁrst simulation,
shown in Fig. 10, consists of a point source in vacuum surroundedby PML (inset). The frequency-domain solver (red squares) shows
rapid, near-exponential convergence, while the error in the timedomain method (blue circles) decreases far more gradually (in fact,
only polynomially). A much more challenging problem is to obtain the frequency-domain response of a cavity (ring resonator)
with multiple long-lived resonant modes: in the time domain,these modes require a long simulation (∼ Q ) to reach steadystate, whereas in the frequency domain the resonances correspondto poles (near-zero eigenvalues of A) that increase the condition
number and hence slow convergence [63]. Fig. 11 shows the results for a ring resonator cavity with multiple closely-spaced resonant modes, excited at one of the resonant frequencies (inset)—
although both frequency- and time-domain methods take longer toconverge than for the non-resonant case of Fig. 10, the advantageof the frequency-domain’s exponential convergence is even moreclear. The convergence is accelerated in frequency domain by using L = 10 (green diamonds) rather than L = 2 (at the expense of
more storage). In time domain, the convergence is limited by the
decay of high- Q modes at other frequencies, and the impact ofthese modes can be reduced by turning on the constant-frequency
source more gradually (magenta triangles, hyperbolic-tangent turnon of the source over 175 optical periods).
This is by no means the most sophisticated possible frequencydomain solver. For example, we currently do not use any preconditioner for the iterative scheme [63]. In two dimensions,a sparse-direct solver may be far more eﬃcient than an iterative
scheme [63]. The key point, however, is that programmer time
is much more expensive than computer time, and this techniqueallows us to obtain substantial improvements in solving frequencydomain problems with only minimal changes to an existing FDTD
program.
6. User interface and scripting
In designing the style of user interaction in Meep, we wereguided by two principles. First, in research or design one hardly
ever needs just one simulation—one almost always performs awhole series of simulations for a class of related problems, exploring the parameter dependencies of the results, optimizing some

698

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702

Fig. 11. Root-mean-square error in ﬁelds in response to a constant-frequency point source exciting one of several resonant modes of a dielectric ring resonator (inset,
ε = 11.56), for frequency-domain solver (red squares, adapted from Meep time-stepping code) vs. time-domain method (magenta triangles, running until transients decayaway). Green diamonds show frequency-domain BiCGSTAB-L solver for ﬁve times more storage, accelerating convergence. Blue circles show time-domain method for a more
gradual turn-on of source, which avoids exciting long-lived resonances at other frequencies. (For interpretation of the references to color in this ﬁgure legend, the reader isreferred to the web version of this article.)

output as a function of the input parameters, or looking at thesame geometry under a sequence of different stimuli. Second,
there is the Unix philosophy: “Write programs that do one thing
and do it well” [66]—Meep should perform electromagnetic simulations, while for additional functionality it should be combinedwith other programs and libraries via standard interfaces like ﬁlesand scripts.
Both of these principles argue against the graphical CADstyle interface common in commercial FDTD software. First, whilegraphical interfaces provide a quick and attractive route to setting up a single simulation, they are not so convenient for aseries of related simulations. One commonly encounters problems where the size/position of certain objects is determined by
the size/position of other objects, where the number of objects isitself a parameter (such as a photonic-crystal cavity surrounded
by a variable number of periods [22]), where the length of thesimulation is controlled by a complicated function of the ﬁelds,
where one output is optimized as a function of some parameter,
and many other situations that become increasingly cumbersome
to express via a set of graphical tools and dialog boxes. Second,
we don’t want to write a mediocre CAD program—if we wanted
to use a CAD program, we would use a professional-quality one,
export the design to a standard interchange format, and write aconversion program to turn this format into what Meep expects.
The most ﬂexible and self-contained interface is, instead, to allow
the user to control the simulation via an arbitrary program. Meep
allows this style of interaction at two levels: via a low-level C++ interface, and via a standard high-level scripting language (Scheme)
implemented by an external library (GNU Guile). The potential
slowness of the scripting language is irrelevant because all of theexpensive parts of the FDTD calculation are implemented in C/C++.
The high-level scripting interface to Meep is documented in detail, with several tutorials, on the Meep web page (http://ab-initio.
mit.edu/meep), so we restrict ourselves to a single short example
in order to convey the basic ﬂavor. This example, in Fig. 12, computes the (2d) ﬁelds in response to a point source located withina dielectric waveguide. We ﬁrst set the size of the computational

cell to 16 × 8 (via geometry-lattice, so-called because it determines the lattice vectors in the periodic case)—recall that theinterpretation of the unit of distance is arbitrary and up to the user
(it could be 16 μm × 8 μm, in which case the frequency units arec /μm, or 16 mm × 8 mm with frequency units of c /mm, or any
other convenient distance unit). Let us call this arbitrary unit of
distance a. Then we specify the geometry within the cell as a listof geometric objects like blocks, cylinders, etc.—in this case by asingle block deﬁning the waveguide with ε = 12—or optionally by
an arbitrary user-deﬁned function ε (x, y ) (and μ, etc.). A layer of
PML is then speciﬁed around the boundaries with thickness 1; thislayer lies inside the computational cell and overlaps the waveguide,
which is necessary in order to absorb waveguide modes when theyreach the edge of the cell. We add a point source, in this casean electric-current source J in the z-direction (sources of arbitrary
spatial proﬁle can also be speciﬁed). The time-dependence of thesource is a sharp turn-on to a continuous-wave source cos(ωt ) at
the beginning of the simulation; gradual turn-ons, Gaussian pulses,
or arbitrary user-speciﬁed functions of time can also be speciﬁed.
The frequency is 0.15 in units of c /a, corresponding to a vacuum
wavelength λ = a/0.15 (e.g. λ ≈ 6.67 μm if a = 1 μm). We set theresolution to 10 pixels per unit distance (10 pixels/a), so that theentire computational cell is 160 × 80 pixels, and then run for 200
time units (units of a/c), corresponding to 200 × 0.15 = 30 optical
periods. We output the dielectric function at the beginning, andthe E z ﬁeld at the end.
In keeping with the Unix philosophy, Meep is not a plotting
program; instead, it outputs ﬁelds and related data to the standard HDF5 format for scientiﬁc datasets [67], which can be readby many other programs and visualized in various ways. (We also
provide a way to effectively “pipe” the HDF5 output to an external program within Meep: for example, to output the HDF5 ﬁle,convert it immediately to an image with a plotting program, andthen delete the HDF5 ﬁle; this is especially useful for producing
animations consisting of hundreds of frames.)
Another important technique to maintain ﬂexibility is that ofhigher-order functions [68]: wherever it is practical, our functions

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702

699

Fig. 12. A simple Meep example showing the E z ﬁeld in a dielectric waveguide (ε = 12) from a point source at a given frequency. A plot of the resulting ﬁeld (blue/white/red =positive/zero/negative) is in the background, and in the foreground is the input ﬁle in the high-level scripting interface (the Scheme language). (For interpretation of thereferences to color in this ﬁgure legend, the reader is referred to the web version of this article.)take functions as arguments instead of (or in addition to) numbers.
Thus, for example, instead of specifying special input codes for all
possible source distributions in space and time, we simply allow
a user-deﬁned function to be used. More subtly, the arguments
output-epsilon and output-efield-z to the run-until
function in Fig. 12 are actually functions themselves: we allow
the user to pass arbitrary “step functions” to run-until that are
called after every FDTD timestep and which can perform arbitrarycomputations on the ﬁelds as desired (or halt the computation if adesired condition is reached). The output-efield-z is simply
a predeﬁned step function that outputs E z . These step-functions
can be modiﬁed by transformation functions like at-end, whichtake step functions as arguments and return a new step functionthat only calls the original step functions at speciﬁed times (atthe end of the simulation, or the beginning, or at certain intervals, for example). In this way, great ﬂexibility in the output andcomputations is achieved. One can, for example, output a givenﬁeld component only at certain time intervals after a given time,and only within a certain subvolume or slice of the computational
cell, simply by composing several of these transformations. One
can even output an arbitrary user-deﬁned function of the ﬁelds instead of predetermined components.
There is an additional subtlety when it comes to ﬁeld output,
because of the Yee lattice in which different ﬁeld components are
stored at different points; presented in this way to the user, it
would be diﬃcult to perform post-processing involving multiple
ﬁeld components, or even to compare plots of different ﬁeld components. As mentioned in Section 3 and again in Section 7.2, therefore, the ﬁeld components are automatically interpolated from theYee grid onto a ﬁxed “centered” grid in each pixel when exported
to a ﬁle.
Although at a simplistic level the input format can just be considered as a ﬁle format with a lot of parentheses, because Scheme
is a full-ﬂedged programming language one can control the simulation in essentially arbitrary ways. Not only can one write loops
and use arithmetic to deﬁne the geometry and the relationships
between the objects or perform parameter sweeps, but we also expose external libraries for multivariable optimization, integration,root-ﬁnding, and other tasks in order that they can be coupledwith simulations.

Parallelism is completely transparent to the user: exactly thesame input script is fed to the parallel version of Meep (writtenwith the MPI message-passing standard for distributed-memory
parallelism [69]) as to the serial version, and the distribution of
the data across processors and the collection of results is handled
automatically.
7. Abstraction versus performance
In an FDTD simulation, essentially just one thing has to be fast:
inner loops over all the grid points or some large fraction thereof.
Everything else is negligible in terms of computation time (but notprogrammer time!), so it can use high-level abstractions without
penalty—for example, the use of a Scheme interpreter as the userinterface has no performance consequences for a typical computation, because the inner loops are not written in Scheme.4 For these
inner loops, however, there is a distinct tension between abstraction (or simplicity) and performance, and in this section we discuss
some of the tradeoffs that result from this tension and the choicesthat have been made in Meep.
The primacy of inner loops means that some popular principles
of abstraction must be discarded. A few years ago, a colleague ofours attempted to write a new FDTD program in textbook objectoriented C++ style: every pixel in the grid would be an object, andevery type of material would be a subclass overriding the necessary timestepping and ﬁeld-access operations. Timestepping would
consist of looping over the grid, calling some “step” method ofeach object, so that objects of different materials (magnetic, dielectric, nonlinear, etc.) would dynamically apply the correspondingﬁeld-update procedures. The result of this noble experiment wasa working program but a performance failure, many times slower
than the aging Fortran software it was intended to replace: theperformance overhead of object dereferencing, virtual method dis-

4
The exception to this rule is when the user supplies a Scheme function and asksthat it be evaluated for every grid point, for example to integrate some functionof the ﬁelds. If this is done frequently during the simulation, it is slow; in thesecircumstances, however, the user can replace the Scheme function with one written
in C/C++ if needed. This is rare because most such functions that might be usedfrequently during a simulation, such as energy or ﬂux, are already supplied in C/C++
within Meep.

700

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702

patch, and function calls in the inner loop overwhelmed all otherconsiderations. In Meep, each ﬁeld’s components are stored as simple linear arrays of ﬂoating-point numbers in row-major (C) order(parallel-array data structures worthy of Fortran 66), and there are
separate inner loops for each type of material (more on this below). In a simple experiment on a 2.8 GHz Intel Core 2 CPU, merely
moving the if statements for the different material types into
these inner loops decreased Meep’s performance by a factor of two
in a typical 3d calculation and by a factor of six in 2d (where thecalculations are simpler and hence the overhead of the conditionals is more signiﬁcant). The cost of the conditionals, including thecost of mispredicted branches and subsequent pipeline stalls [70]along with the frustration of compiler unrolling and vectorization,
easily overwhelmed the small cost of computing, e.g., ∇ × H at asingle point.
7.1. Timestepping and cache tradeoffs
One of the dominant factors in performance on modern computer systems is not arithmetic, but memory: random memory
access is far slower than arithmetic, and the organization of memory into a hierarchy of caches is designed to favor locality of access [70]. That is, one should organize the computation so thatas much work as possible is done with a given datum once it isfetched (temporal locality) and so that subsequent data that areread or written are located nearby in memory (spatial locality).
The optimal strategies to exploit both kinds of locality, however,appear to lead to sacriﬁces of abstraction and code simplicity sosevere that we have chosen instead to sacriﬁce some potential performance in the name of simplicity.
As it is typically described, the FDTD algorithm has very little temporal locality: the ﬁeld at each point is advanced in time by
t, and then is not modiﬁed again until all the ﬁelds at every otherpoint in the computational cell have been advanced. In order togain temporal locality, one must employ asynchronous timestepping:essentially, points in small regions of space are advanced several
steps in time before advancing points far away, since over a shorttime interval the effects of far-away points cannot be felt. A detailed analysis of the characteristics of this problem, as well as abeautiful “cache-oblivious” algorithm that automatically exploits acache of any size for grids of any dimensionality, is described inRef. [71]. On the other hand, an important part of Meep’s usability
is the abstraction that the user can perform arbitrary computations
or output using the ﬁelds in any spatial region at any time, which
seems incompatible with the ﬁelds at different points in space being out-of-sync until a predetermined end of the computation. The
bookkeeping diﬃculty of reconciling these two viewpoints led us
to reject the asynchronous approach, despite its potential beneﬁts.
However, there may appear to be at least a small amount of
temporal locality in the synchronous FDTD algorithm: ﬁrst B is advanced from ∇ × E, then H is computed from B and μ, then D isadvanced from ∇ × H, then E is computed from D and ε . Since
most ﬁelds are used at least once after they are advanced, surely
the updates of the different ﬁelds can be merged into a single loop,
for example advancing D at a point and then immediately computing E at the same point—the D ﬁeld need not even be stored.
Furthermore, since by merging the updates one is accessing several ﬁelds around the same point at the same time, perhaps one
can gain spatial locality by interleaving the data, say by storing anarray of (E, H, ε , μ) tuples instead of separate arrays. Meep does
not do either of these things, however, for two reasons, the ﬁrstof which is more fundamental. As is well known, one cannot easily merge the B and H updates with the D and E updates at thesame point, because the discretized ∇ × operation is nonlocal (involves multiple grid points)—this is why one normally updates H
everywhere in space before updating D from ∇ × H, because incomputing ∇ × H one uses the values of H at different grid points
and all of them must be in sync. A similar reasoning, however,applies to updating E from D and H from B, once the possibilityof anisotropic materials is included—because the Yee grid stores
different ﬁeld components at different locations, any accurate handling of off-diagonal susceptibilities must also inevitably involveﬁelds at multiple points (as in Ref. [39]). To handle this, D must be
stored explicitly and the update of E from D must take place after D has been updated everywhere, in a separate loop. And sinceeach ﬁeld is updated in a separate loop, the spatial-locality motivation to merge the ﬁeld data structures rather than using parallel
arrays is largely removed.
Of course, not all simulations involve anisotropic materials—
although they appear even in many simulations with nominallyisotropic materials thanks to the subpixel averaging discussed in
Section 3—but this leads to the second practical problem withmerging the E and D (or H and B) update loops: the combinatorial explosion of the possible material cases. The update of D
from ∇ × H must handle 16 possible cases, each of which is a separate loop (see above for the cost of putting conditionals inside theloops): with or without PML (4 cases, depending upon the numberof PML conductivities and their orientation relative to the ﬁeld),with or without conductivity, and with the derivative of two H
components (3d) or only one H component (2d TE polarization).
The update of E from D involves 12 cases: with or without PML(2 cases, distinct from those in the D update), the number of offdiagonal ε −1 components (3 cases: 0, 1, or 2), and with or without
nonlinearity (2 cases). If we attempted to join these into a single loop, we would have 16 × 12 = 192 cases, a code-maintenanceheadache. (Note that the multiplicity of PML cases comes from thefact that, including the corners of the computational cell, we might
have 0 to 3 directions of PML, and the orientation of the PML directions relative to a given ﬁeld component matters greatly.)
The performance penalty of separate E and D (or H and B)
updates appears to be modest. Even if, by somehow merging theloops, one assumes that the time to compute E = ε −1 D could become zero, benchmarking the relative time spent in this operation
indicates that a typical 3d transmission calculation would be accelerated by only around 30% (and less in 2d).
7.2. The loop-in-chunks abstraction
Finally, let us brieﬂy mention a central abstraction that, while
not directly visible to end-users of Meep, is key to the eﬃciency
and maintainability of large portions of the software (ﬁeld output,
current sources, ﬂux/energy computations and other ﬁeld integrals,
and so on). The purpose of this abstraction is to mask the complexity of the partitioning of the computational cell into overlapping chunks connected by symmetries, communication, and otherboundary conditions as described in Section 2.
Consider the output of the ﬁelds at a given timestep to an HDF5
dataﬁle. Meep provides a routine get-field-pt that, given a
point in space, interpolates it onto the Yee grid and returns a desired ﬁeld component at that point. In addition to interpolation,
this routine must also transform the point onto a chunk that is
actually stored (using rotations, periodicity, etc.) and communicate
the data from another processor if necessary. If the point is ona boundary between two chunks, the interpolation process may
involve multiple chunks, multiple rotations, etc., and communications from multiple processors. Because this process involves onlya single point, it is not easily parallelizable. Now, to output theﬁelds everywhere in some region to a ﬁle, one approach is to simply call get-field-pt for every point in a grid for that regionand output the results, but this turns out to be tremendously slow
because of the repeated transformations and communications forevery single point. We nevertheless want to interpolate ﬁelds forA.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702

output rather than dumping the raw Yee grid, because it is mucheasier for post-processing if the different ﬁeld components are interpolated onto the same grid; also, to maintain transparency offeatures like symmetry one would like to be able to output thewhole computational cell (or an arbitrary subset) even if only apart of it is stored. Almost exactly the same problems arise for integrating things like ﬂux E × H or energy or user-deﬁned functionsof the ﬁelds (noting that functions combining multiple ﬁeld components require interpolation), and also for implementing volume(or line, or surface) sources which must be projected onto the gridin some arbitrary volume.
One key to solving this diﬃculty is to realize that, when theﬁeld in some volume V is needed (for output, integration, and soon), the rotations, communications, etc. for points in V are identical for all the points in the intersection of V with some chunk(or one of its rotations/translations). The second is to realize that,when interpolation is needed, there is a particular grid for whichinterpolation is easy: for owned points of the centered grid (Section 2) lying at the center of each pixel, it is always possible tointerpolate from ﬁelds on any Yee grid without any inter-chunkcommunication and by a simple equal-weight averaging of at most
2d points in d dimensions.
The combination of these two observations leads to the loopin-chunks abstraction. Given a (convex rectilinear) volume V anda given grid (either centered, or one of the Yee-ﬁeld grids),
it computes the intersection of all the chunks and their rotations/translations with V . For each intersection it invokes a callerspeciﬁed function, passing the portion of the chunk, the necessary
rotations (etc.) of the ﬁelds, and interpolation weights (if needed,
for the boundary of V ). That function then processes the speciﬁed
portion of the chunk (for example, outputting it to the corresponding portion of a ﬁle, or integrating the desired ﬁelds). All of thiscan proceed in parallel (with each processor considering only thosechunks stored locally). This is (relatively) fast because the rotations,interpolations, and so on are computed only once per chunk intersection, while the inner loop over all grid points in each chunk canbe as tight as necessary. Moreover, all of the rather complicatedand error-prone logic involved in computing V ’s intersection withthe chunks (e.g., special care is required to ensure that each conceptual grid point is processed exactly once despite chunk overlaps
and symmetries) is localized to one place in the source code; ﬁeld
output, integration, sources, and other functions of the ﬁelds are
isolated from this complexity.
8. Concluding remarks
We have reviewed in this paper a number of the unusual implementation details of Meep that distinguish our software package
from standard textbook FDTD methods. Beginning with a discussion of the fundamental structural unit of chunks that constitutethe Yee grid and enable parallelization: we provided an overview
of Meep’s core design philosophy of creating an illusion of continuous space and time for inputs and outputs; we explained andmotivated the somewhat unusual design intricacies of nonlinear
materials and PMLs; we discussed important aspects of Meep’s
computational methods for ﬂux spectra and resonant modes; wedemonstrated the formulation of a frequency-domain solver requiring only minimal modiﬁcations to the underlying time-stepping
algorithm. In addition to the inner workings of Meep’s internal
structure, we reviewed how such features are accessible to users
via an external scripting interface.
We believe that a free/open-source, full-featured FDTD package like Meep can play a vital role in enabling new research inelectromagnetic phenomena. Not only does it provide a low barrier to entry for standard FDTD simulations, but the simplicity of
the FDTD algorithm combined with access to the source code of-701

fers an easy route to investigate new physical phenomena coupledwith electromagnetism. For example, we have colleagues working
on coupling multi-level atoms to electromagnetism within Meepfor modeling lasing and saturable absorption, adapting publishedtechniques from our and other groups [16–20], but also including
new physics such as diffusion of excited gases. Other colleagues
have modiﬁed Meep for modeling gyromagnetic media in order todesign new classes of “one-way” waveguides [72]. Meep is evenbeing used to simulate the quantum phenomena of Casimir forces(from quantum vacuum ﬂuctuations, which can be computed fromclassical Green’s functions) [73,74]—in fact, this was possible without any modiﬁcations of the Meep code due to the ﬂexibility of
Meep’s scripting interface. We hope that other researchers, withthe help of the understanding of Meep’s architecture that this paper provides, will be able to adapt Meep to future phenomena that
we have not yet envisioned.
Acknowledgements
This work was supported in part by the Materials Research
Science and Engineering Center program of the National Science
Foundation under Grant Nos. DMR-9400334 and DMR-0819762, by
the Army Research Oﬃce through the Institute for Soldier Nanotechnologies under contract DAAD-19-02-D0002, and also by Dr.
Dennis Healy of DARPA MTO under award N00014-05-1-0700 administered by the Oﬃce of Naval Research. We are also grateful
to A.W. Rodriguez and A.P. McCauley for their efforts to generalize Meep for quantum-Casimir problems, S.L. Chua for his work on
dispersive and multi-level materials, and Y. Chong for early support
in Meep for gyrotropic media.
References[1] A. Taﬂove, S.C. Hagness, Computational Electrodynamics: The Finite-Difference
Time-Domain Method, 3rd ed., Artech, Norwood, MA, 2005.
[2] K.S. Kunz, R.J. Luebbers, The Finite-Difference Time-Domain Method for Electromagnetics, CRC Press, Boca Raton, 1993.
[3] D.M. Sullivan, Electromagnetic Simulation Using the FDTD Method, Wiley–IEEE
Press, New York, 2000.
[4] A. Elsherbeni, V. Demir, The Finite Difference Time Domain Method for Electromagnetics: With MATLAB Simulations, SciTech, Rayleigh, NC, 2009.
[5] W. Yu, R. Mittra, T. Su, Y. Liu, X. Yang, Parallel Finite-Difference Time-Domain
Method, Artech House, Norwood, MA, 2006.
[6] M.N. Sadiku, Numerical Techniques in Electromagnetics, 2nd ed., CRC, 2000.
[7] J. Jin, The Finite Element Method in Electromagnetics, 2nd ed., Wiley–IEEE
Press, 2002.
[8] K. Yasumoto (Ed.), Electromagnetic Theory and Applications for Photonic Crystals, CRC, 2005.
[9] S. Rao, D. Wilton, A. Glisson, Electromagnetic scattering by surfaces of arbitrary
shape, IEEE Trans. Antennas and Propagation 30 (3) (1982) 409–418.
[10] K. Umashankar, A. Taﬂove, S. Rao, Electromagnetic scattering by arbitrary
shaped three-dimensional homogeneous lossy dielectric objects, IEEE Trans.
Antennas and Propagation 34 (6) (1986) 758–766.
[11] M. Bonnet, Boundary Integral Equation Methods for Solids and Fluids, Wiley,
1999.
[12] W.C. Chew, J.-M. Jin, E. Michielssen, J. Song (Eds.), Fast and Eﬃcient Algorithms
in Computational Electromagnetics, Artech House, 2000.
[13] S.G. Johnson, M. Ibanescu, M.A. Skorobogatiy, O. Weisberg, J.D. Joannopoulos,
Y. Fink, Perturbation theory for Maxwell’s equations with shifting materialboundaries, Phys. Rev. E 65 (2002) 066611.
[14] L. Zhang, J. Lee, A. Farjadpour, J. White, S. Johnson, A novel boundary element
method with surface conductive absorbers for 3-D analysis of nanophotonics, in: Microwave Symposium Digest, 2008 IEEE MTT-S International, 2008,
pp. 523–526.
[15] S.G. Johnson, J.D. Joannopoulos, Block-iterative frequency-domain methods for
Maxwell’s equations in a planewave basis, Opt. Express 8 (3) (2001) 173–190.
[16] R.W. Ziolkowski, J.M. Arnold, D.M. Gogny, Ultrafast pulse interactions with twolevel atoms, Phys. Rev. A 52 (4) (1995) 3082–3094.
[17] A.S. Nagra, R.A. York, FDTD analysis of wave propagation in nonlinear absorbing
and gain media, IEEE Trans. Antennas and Propagation 46 (3) (1998) 334–340.
[18] S.-H. Chang, A. Taﬂove, Finite-difference time-domain model of lasing action
in a four-level two-electron atomic system, Opt. Express 12 (16) (2004) 3827–
3833.

702

A.F. Oskooi et al. / Computer Physics Communications 181 (2010) 687–702[19] Y. Huang, S.-T. Ho, Computational model of solid-state, molecular, or atomic
media for FDTD simulation based on a multi-level multi-electron system governed by Pauli exclusion and Fermi–Dirac thermalization with application tosemiconductor photonics, Opt. Express 14 (8) (2006) 3569–3587.
[20] P. Bermel, E. Lidorikis, Y. Fink, J.D. Joannopoulos, Active materials embedded
in photonic crystals and coupled to electromagnetic radiation, Phys. Rev. B 73
(2006) 165125.
[21] J.D. Jackson, Classical Electrodynamics, 3rd ed., Wiley, New York, 1998.
[22] J.D. Joannopoulos, S.G. Johnson, R.D. Meade, J.N. Winn, Photonic Crystals: Molding the Flow of Light, 2nd ed., Princeton Univ. Press, 2008.
[23] K.S. Yee, Numerical solution of initial boundary value problems involving
Maxwells Equations in isotropic media, IEEE Trans. Antennas and Propagation 14 (3) (1966) 302–307.
[24] M.J. Berger, J. Oliger, Adaptive mesh reﬁnement for hyperbolic partial differential equations, J. Comput. Phys. 53 (1984) 484–512.
[25] I.S. Kim, W.J.R. Hoefer, A local mesh reﬁnement algorithm for the time domainﬁnite difference method using Maxwell’s curl equations, IEEE Trans. Microwave
Theory Tech. 38 (6) (1990) 812–815.
[26] S.S. Zivanovic, K.S. Yee, K.K. Mei, A. Subgridding, Method for the time-domainﬁnite-difference method to solve Maxwell’s equations, IEEE Trans. Microwave
Theory Tech. 39 (3) (1991) 471–479.
[27] M. Okoniewski, E. Okoniewska, M.A. Stuchly, Three-dimensional subgridding
algorithm for FDTD, IEEE Trans. Antennas and Propagation 45 (3) (1997) 422–
429.
[28] C. Lin, L. Snyder, Principles of Parallel Programming, Addison–Wesley, 2008.
[29] T. Inui, Y. Tanabe, Y. Onodera, Group Theory and Its Applications in Physics,
Springer-Verlag, Terlos, 1996.
[30] U. Trottenberg, C.W. Oosterlee, A. Schuller, Multigrid, Academic Press, 2000.
[31] L. Landau, L. Pitaevskii, E. Lifshitz, Electrodynamics of Continuous Media, 2nded., Butterworth–Heinemann, 1984.
[32] C. Luo, M. Ibanescu, S.G. Johnson, J.D. Joannopoulos, Cerenkov radiation in photonic crystals, Science 299 (2003) 368–371.
[33] A. Farjadpour, D. Roundy, A. Rodriguez, M. Ibanescu, P. Bermel, J. Joannopoulos,
S. Johnson, G. Burr, Improving accuracy by sub-pixel smoothing in the ﬁnitedifference time domain, Opt. Lett. 31 (2006) 2972–2974.
[34] A.F. Oskooi, C. Kottke, S.G. Johnson, Accurate ﬁnite-difference time-domain simulation of anisotropic media by subpixel smoothing, Opt. Lett. 34 (18) (2009)
2778–2780.
[35] A. Ditkowski, K. Dridi, J.S. Hesthaven, Convergent Cartesian grid methods for
Maxwell’s equations in complex geometries, J. Comput. Phys. 170 (2001) 39–
80.
[36] C. Kottke, A. Farjadpour, S.G. Johnson, Perturbation theory for anisotropic dielectric interfaces, and application to subpixel smoothing of discretized numerical methods, Phys. Rev. E 77 (2008) 036611.
[37] S. Dey, R. Mittra, A conformal ﬁnite-difference time-domain technique formodeling cylindrical dielectric resonators, IEEE Trans. Microwave Theory
Tech. 47 (9) (1999) 1737–1739.
[38] R.D. Meade, A.M. Rappe, K.D. Brommer, J.D. Joannopoulos, O.L. Alerhand, Accurate theoretical analysis of photonic band-gap materials, Phys. Rev. B 48 (1993)
8434–8437;
S.G. Johnson, Phys. Rev. B 55 (1997) 15942, Erratum.
[39] G. Werner, J. Cary, A stable FDTD algorithm for non-diagonal anisotropic dielectrics, J. Comput. Phys. 226 (2007) 1085–1101.
[40] P. Mezzanotte, L. Roselli, R. Sorrentino, A. Simple, Way to model curved
metal boundaries in FDTD algorithm avoiding staircase approximation, IEEE Microwave Guided Wave Lett. 5 (8) (1995) 267–269.
[41] J. Anderson, M. Okoniewski, S.S. Stuchly, Practical 3-D contour/staircase treatment of metals in FDTD, IEEE Microwave Guided Wave Lett. 6 (3) (1996) 146–
148.
[42] A. Deinega, I. Valuev, Subpixel smoothing for conductive and dispersive media
in the ﬁnite-difference time-domain method, Opt. Lett. 32 (23) (2007) 3429–
3431.
[43] P. Leung, S. Liu, K. Young, Completeness and time-independent perturbationof the quasinormal modes of an absorptive and leaky cavity, Phys. Rev. A 49
(1994) 3982–3989.
[44] S.G. Johnson, M.L. Povinelli, M. Soljačić, A. Karalis, S. Jacobs, J.D. Joannopoulos,
Roughness losses and volume-current methods in photonic-crystal waveguides,
Appl. Phys. B 81 (2005) 283–293.
[45] A.F. Oskooi, L. Zhang, Y. Avniel, S.G. Johnson, The failure of perfectly
matched layers, and towards their redemption by adiabatic absorbers, Opt. Express 16 (15) (2008) 11376–11392.

[46] Q.H. Liu, J.Q. He, Quasi-PML for waves in cylindrical coordinates, Microwave
and Optical Tech. Lett. 19 (2) (1998) 107–111.
[47] N. Bloembergen, Nonlinear Optics, W.A. Benjamin, New York, 1965.
[48] G.P. Agrawal, Nonlinear Fiber Optics, 3rd ed., Academic Press, San Diego, 2001.
[49] A. Rodriguez, M. Soljačić, J.D. Joannopoulos, S.G. Johnson, χ (2) and χ (3) harmonic generation at a critical power in inhomogeneous doubly resonant cavities, Opt. Express 15 (12) (2007) 7303–7318.
[50] R.W. Boyd, Nonlinear Optics, Academic Press, London, UK, 1992.
[51] J. Baker, A. George, P. Graves-Morris, Padé Approximants, 2nd ed., Cambridge
University Press, 1996.
[52] F.L. Teixeira, W.C. Chew, Systematic derivation of anisotropic PML absorbing
media in cylindrical and spherical coordinates, IEEE Microwave Guided Wave
Lett. 7 (11) (1997) 371–373.
[53] J.-Q. He, Q.-H. Liu, A nonuniform cylindrical FDTD algorithm with improved
PML and quasi-PML absorbing boundary conditions, IEEE Trans. Geoscience Remote Sensing 37 (2) (1999) 1066–1072.
[54] D.H. Bailey, P.N. Swartztrauber, A. Fast, Method for the numerical evaluation ofcontinuous Fourier and Laplace transforms, SIAM J. Sci. Comput. 15 (5) (1994)
1105–1110.
[55] T.H. Cormen, C.E. Leiserson, R.L. Rivest, C. Stein, Introduction to Algorithms, 3rded., MIT Press, 2009.
[56] A.V. Oppenheim, R.W. Schafer, Discrete-Time Signal Processing, 3rd ed.,Prentice-Hall, 2009.
[57] W.-H. Guo, W.-J. Li, Y.-Z. Huang, Computation of resonant frequencies andquality factors of cavities by FDTD technique and Padé approximation, IEEE
Microwave and Wireless Comp. Lett. 11 (5) (2001) 223–225.
[58] S. Cabay, D.-K. Choi, Algebraic computations of scaled Padé fractions, SIAM J.
Comput. 15 (1) (1986) 243–270.
[59] A. Rodriguez, M. Ibanescu, J.D. Joannopoulos, S.G. Johnson, Disorder-immune
conﬁnement of light in photonic-crystal cavities, Opt. Lett. 30 (2005) 3192–
3194.
[60] V.A. Mandelshtam, H.S. Taylor, Harmonic inversion of time signals and its applications, J. Chem. Phys. 107 (17) (1997) 6756–6769.
[61] S. Dey, R. Mittra, Eﬃcient computation of resonant frequencies and quality
factors of cavities via a combination of the ﬁnite-difference time-domain technique and the Padé approximation, IEEE Microwave Guided Wave Lett. 8 (12)(1998) 415–417.
[62] A. Christ, H.L. Hartnagel, Three-dimensional ﬁnite-difference method forthe analysis of microwave-device embedding, IEEE Trans. Microwave Theory
Tech. 35 (8) (1987) 688–696.
[63] R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R.
Pozo, C. Romine, H.V. der Vorst, Templates for the Solution of Linear Systems:
Building Blocks for Iterative Methods, SIAM, Philadelphia, PA, 1994.
[64] G.L.G. Sleijpen, D.R. Fokkema, BiCGSTAB(L) for linear equations involving unsymmetric matrices with complex spectrum, Electron. Trans. Numer. Anal. 1(1993) 11–32.
[65] G.L.G. Sleijpen, H.A. van der Vorst, D.R. Fokkema, BiCGstab(L) and other hybrid
Bi-CG methods, Numer. Algorithms 7 (1994) 75–109.
[66] P.H. Salus, A. Quarter, Century of UNIX, Addison–Wesley, Reading, MA, 1994.
[67] M. Folk, R.E. McGrath, N. Yeager, HDF: An update and future directions, in:Proc. 1999 Geoscience and Remote Sensing Symposium (IGARSS), vol. 1, IEEE
Press, Hamburg, Germany, 1999, pp. 273–275.
[68] H. Abelson, G.J. Sussman, Structure and Interpretation of Computer Programs,
MIT Press, Cambridge, MA, 1985.
[69] T.M. Forum, MPI: A Message Passing Interface, in: Supercomputing ’93, Portland, OR, 1993, pp. 878–883.
[70] J.L. Hennessy, D.A. Patterson, Computer Architecture: A Quantitative Approach,
3rd ed., Elsevier, San Francisco, CA, 2003.
[71] M. Frigo, V. Strumpen, The memory behavior of cache oblivious stencil computations, J. Supercomputing 39 (2) (2007) 93–112.
[72] Z. Wang, Y.D. Chong, J.D. Joannopoulos, M. Soljačić, Reﬂection-free one-wayedge modes in a gyromagnetic photonic crystal, Phys. Rev. Lett. 100 (2008)
013905.
[73] A.W. Rodriguez, A.P. McCauley, J.D. Joannopoulos, S.G. Johnson, Casimir forcesin the time domain: Theory, Phys. Rev. A 80 (2009) 012115.
[74] A.P. McCauley, A.W. Rodriguez, J.D. Joannopoulos, S.G. Johnson, Casimir forcesin the time domain: II. Applications, arXiv.org e-Print archive, arXiv:0906.5170,
2009.