INTEGRAL EQUATION AND DISCONTINUOUS GALERKIN METHODS
FOR THE ANALYSIS OF LIGHT-MATTER INTERACTION
By
Andrew David Baczewski

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Electrical Engineering - Doctor of Philosophy
Physics - Doctor of Philosophy
2013

ABSTRACT
INTEGRAL EQUATION AND DISCONTINUOUS GALERKIN METHODS
FOR THE ANALYSIS OF LIGHT-MATTER INTERACTION
By
Andrew David Baczewski
Light-matter interaction is among the most enduring interests in the physical sciences. As
a result, the understanding and control of this physics is of paramount importance to the
design of myriad technologies ranging from stained glass, to molecular sensing and characterization techniques, to quantum computers. The development of complex engineered
systems that exploit this physics is predicated at least partially upon in silico design and
optimization that properly capture the light-matter coupling. In this thesis, details of computational frameworks that enable this type of analysis, based upon both Integral Equation
and Discontinuous Galerkin formulations, will be explored. The primary focus will be on
the development of efficient and accurate software, with results corroborating both. The
secondary focus will be on the use of these tools in the analysis of a number of exemplary
systems.

ACKNOWLEDGMENT

First and foremost, I am grateful for the guidance of my committee, Shanker Balasubramaniam, Bhanu Mahanti, Phil Duxbury, Tim Hogan, and Leo Kempel. Shanker has been
a tremendously supportive, brilliant, and motivating adviser over these past 6 years. I am
grateful for his patience, friendly mentorship, and the many hours he invested in making
me a better computational scientist. Bhanu has shared a great deal about thinking like
and otherwise operating as a physicist. While this Dissertation is largely free of our work
together, it has been an intellectual joy to learn about electronic structure and more general
condensed matter theory under his supervision. Phil is largely responsible for introducing
me to computational physics as an undergraduate nuking nanowires, and played a big role
in helping me confidently transition into Physics. Tim has been a perpetual source of inspiration since my undergraduate years, when he taught me that scientific research should be a
means for finding and fixing the hurts that plague the world. Leo has been a very practical
mentor, who has always encouraged me to push beyond simulating spheres and PEC strips.
Thank you all for your time and support.
Other faculty who have been supportive of me as a researcher and a person include Ben
Levine, Carlo Piermarocchi, Ed Rothwell, Prem Chahal, and Virginia Ayres. My interactions
with them have been immensely useful, and all played an important role in shaping me as
a young scientist. My labmates in the CEML have played a similar role. Thank you to Lu
Chuan, Jorge Villa-Giron, Vikram Melapudi, Naveen Nair, Ozgur Tuncer, Andrew Pray, Dan
Dault, Collin Meierbachtol, Nick Miller, Steve Hughey, Jie Li, and Mohammad Rawashdeh. I
am proud to have you all as colleagues, and to have shared so many successes and frustrations
iii

over the years.
Lest I forget my cohort in Physics, thank you you to Matt DeNinno, Brent Barker,
Emi Wang, Aaron Hoffer, Brad Barquest, Josh Vredevoogd, Alex Wittig, Ravi Jagasia,
River Huang, Dat Do, and Mal-Soon Lee, for struggling with me, and sharing so much of
yourselves.
I am fortunate to have a number of friends outside of science, who have tolerated my
unusual schedule and various stresses over the years. Thank you to Devin Schnepp, Greg
Ondrus, Anthony Recca, Pat Aderhold, Michael Davidson, Steve Marangos, and Viktor
Nikolovski for keeping me grounded and providing much needed distractions. My girlfriend,
Melissa Grant, has been exceptionally patient and supportive of me during the preparation
of this document. Her love and dedication have been absolutely vital to my success.
My parents, David and Monique Baczewski, have been an uninterrupted source of love,
patience, knowledge, and support, even when I informed that I want to have a job where
‘someone pays me to think’ at the age of 7. Having essentially accomplished that goal, I am
glad that they showed me nothing but encouragement along the way. My sister Margaret
Baczewski and my grandparents, Pat Leveque and Tom Langelier, are also in some measure
responsible for getting me to this point.
This thesis is dedicated to my dear friend, Carl Coppola, whom I met on my first day as
a Freshman in 2003. Over our undergraduate years we studied, lived, and laughed together,
with both of us continuing on at Michigan State for our graduate education. Carl passed
in 2010, and the halls of the Engineering Building are poorer in the absence of his wit,
perseverance, and selflessness. I am truly privileged to have known him, and am honored to
share his memory in these pages.

iv

TABLE OF CONTENTS

LIST OF FIGURES

. . . . . . . . . . . . . . . . . . . . . . . . . . . viii

KEY TO SYMBOLS AND ABBREVIATIONS . . . . . . . . . . . . . .

xii

Chapter 1: Introduction
1.1 Physics of Interest and Technological Relevance . . . . . . . . . . . . . .
1.1.1 Periodic Structures: FSS, ‘Metamaterials’, and Photonic Crystals
1.1.2 Energy Transfer, CQED, and Quantum Information . . . . . . . .
1.2 Overview of Numerical Methods for the Maxwell Equations . . . . . . . .
1.2.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.2 Galerkin Methods . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.3 Comparison of IE and PDE-Based Methods . . . . . . . . . . . .
1.2.4 Fast Methods for IEs . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.5 Discontinuous Galerkin vs. Continuous Galerkin . . . . . . . . . .
1.3 Layout of this Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

1
3
3
4
6
7
12
16
18
20
23
23

Chapter 2: IEs for the Time-Harmonic Maxwell Problem
2.1 Overview of IEs . . . . . . . . . . . . . . . . . . . . . . .
2.2 Scattering from PEC Targets . . . . . . . . . . . . . . .
2.2.1 Statement of the Problem . . . . . . . . . . . . .
2.2.2 Numerical Details . . . . . . . . . . . . . . . . . .
2.3 Scattering from Dielectric Targets . . . . . . . . . . . . .
2.3.1 Statement of the Problem . . . . . . . . . . . . .
2.3.2 Numerical Details . . . . . . . . . . . . . . . . . .
2.4 Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

26
26
26
27
29
31
31
33
34

.
.
.
.
.
.
.

37
37
37
37
39
41
41
45

Chapter 3: Periodic ACE
3.1 Overview of Periodic ACE . . . . . . . . .
3.2 Hierarchical Fast Methods . . . . . . . . .
3.2.1 Tree Codes, FMM, and ACE . . . .
3.2.2 Specialization To Periodic Systems
3.3 Statement of the Problem . . . . . . . . .
3.3.1 Continuous Problem . . . . . . . .
3.3.2 Discrete Problem . . . . . . . . . .
v

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

3.4

3.5
3.6

3.7

3.8

Details of the ACE Algorithm . . . . . . . . . . . . .
3.4.1 Tree Construction . . . . . . . . . . . . . . . .
3.4.2 Interaction Lists . . . . . . . . . . . . . . . . .
3.4.3 Precomputation . . . . . . . . . . . . . . . . .
3.4.4 Tree Traversal . . . . . . . . . . . . . . . . . .
Cost Analysis . . . . . . . . . . . . . . . . . . . . . .
Kernel Results . . . . . . . . . . . . . . . . . . . . . .
3.6.1 Error Convergence . . . . . . . . . . . . . . .
3.6.2 Computational Cost . . . . . . . . . . . . . .
ACE for Scattering from Periodic PEC and Dielectric
3.7.1 PEC Targets . . . . . . . . . . . . . . . . . .
3.7.2 Dielectric Targets . . . . . . . . . . . . . . . .
Figures . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
Targets
. . . . .
. . . . .
. . . . .

Chapter 4: DG for the Transient Maxwell Problem
4.1 Overview of Time Domain DG . . . . . . . . . . . . . .
4.2 Time Domain Maxwell Equations . . . . . . . . . . . .
4.2.1 Statement of the Problem . . . . . . . . . . . .
4.2.2 Numerical Details . . . . . . . . . . . . . . . . .
4.2.2.1 Modal and Nodal DG Basis Functions
4.2.2.2 DG Discretization . . . . . . . . . . .
4.2.2.3 Boundary Conditions . . . . . . . . . .
4.2.2.4 Non-Conformal Discretizations . . . .
4.2.3 Validation . . . . . . . . . . . . . . . . . . . . .
4.3 Figures . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.

46
48
49
51
52
56
59
59
62
63
65
67
70

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

89
89
89
89
91
91
94
97
100
101
105

Chapter 5: DG-QED for Cavity Mediated Energy Transfer
5.1 Overview of Energy Transfer . . . . . . . . . . . . . . . . .
5.2 Cavity QED Formalism . . . . . . . . . . . . . . . . . . . .
5.2.1 Model Hamiltonian and Master Equation . . . . . .
5.2.2 Extraction of Model Parameters . . . . . . . . . . .
5.3 DG-QED Results . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 Spontaneous Emission Above a PEC Half-Space . .
5.3.2 Energy Transfer in a Semiconductor Nanocavity . .
5.4 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . .
5.5 Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

113
113
114
114
116
118
118
120
123
127

. . .
. . . .
. . . .
. . . .
. . .
. . . .

.
.
.
.
.
.

.
.
.
.
.

. .
. .
. .
. .
. .
. .

.
.
.
.
.
.

140
140
140
141
143
145

APPENDICES . . . . . . . . . . . . . . .
A.1 Bloch-Periodic Green’s Functions . . . .
A.1.1 Bloch-Periodicity Condition . . .
A.1.2 Ewald Representation . . . . . .
A.2 Periodic Translation Operator Components
A.3 Time-Domain Field of a Dipole . . . . .
vi

.
.
.
.

.
.
.
.
.
. .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

. . . .
.
. .
. . . .
. . .
. . .
. . . .

BIBLIOGRAPHY

. . . . . . . . . . . . . . . . . . . . . . . . . . .

vii

149

LIST OF FIGURES

1.1

1.2

1.3

2.4

2.5

3.6

3.7

3.8

The Lycurgus cup, a piece of 1700 year old Roman nanotechnology exhibiting
optical dichroism. c Trustees of the British Museum. For interpretation of
the references to color in this and all other figures, the reader is referred to
the electronic version of this dissertation. . . . . . . . . . . . . . . . . . . .
(Left) The domain upon which some continuous problem is defined. (Right)
A tetrahedral mesh of that domain, upon which the discrete problem will be
solved. Each of the tetrahedra into which this geometry is decomposed, is one
of the subdomains, Te . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

24

Illustration of a non-conformal mesh comprised of two hexahedra resolved at
different discretization rates. . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

An exemplary triangular mesh used in calculating scattering from a periodic
array of PEC split ring resonators. . . . . . . . . . . . . . . . . . . . . . . .

35

An exemplary tetrahedral mesh used in calculating scattering from a periodic
array of dielectric spheres. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

An exemplary periodic domain for µ = 2. The central primitive cell, as well
as its 8 closest image cells are illustrated. The support of some arbitrary ρ(r)
is represented by the blue shaded areas. . . . . . . . . . . . . . . . . . . . . .

71

The primitive unit cell (blue) and the 15 implicit image cells (red) involved in
the enlarged tree. The original tree has Nl = 2, whereas the enlarged tree has
Nl = 4. Rather than octal addresses, decimal addresses are given to provide
a simpler view of the ordering. . . . . . . . . . . . . . . . . . . . . . . . . .

72

‘Near’ (blue) and ‘far’ (light blue) interactions for source boxes (dark blue) at
different levels of the tree. The contents of the red image cells are not actually
stored, but need to be considered in constructing interaction lists. . . . . . .

73

viii

3.9

Illustration of the exemplary scenario described above. Only the primitive cell
is represented. The leaf level domains are highlighted in dark blue, whereas
the parent domains are highlighted in lighter blue. . . . . . . . . . . . . . .

74

3.10 Convergence of the relative error in reconstructing the farfield contribution to
the potential for both Helmholtz and Yukawa systems with µ = 1, µ = 2, and
µ = 3, with κ = π. 8192 sources are distributed over the unit cube, which is
also the primitive cell for all tests. . . . . . . . . . . . . . . . . . . . . . . .

75

3.11 Farfield error as a function of λ = 2π/κ for 1000 points randomly distributed
over the unit cube at a fixed ACE order, P = 6, and µ = 2. . . . . . . . . . .

76

3.12 Demonstration of sublinear scaling of the precomputation time for ACE expansions of different orders. These results are generated for a Helmholtz
problem with µ = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

3.13 Demonstration of linear scaling of the tree traversal timing for ACE expansions
of different orders. These results are generated for a Helmholtz problem with
µ = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

3.14 (Top) Plot of surface currents on Jerusalem cross. (Bottom) Power transmis√
x ˆ
sion spectrum at oblique incidence (θ = 45◦ , φ = 0◦ ) with E0 = 2/2 (ˆ − z ).
TFQMR with a diagonal preconditioner was utilized to achieve a residual
≤ 0.001%. P = 7 harmonics and Nl = 3 levels were used in the ACE algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

3.15 (Top) Plot of surface currents on an Omega split ring structure. (Bottom)
Power transmission spectrum for an Omega SRR at normal incidence (θ =
0◦ , φ = 0◦ ) for E0 = x. GMRES with a diagonal preconditioner was utilized
ˆ
to achieve a residual ≤ 0.01%. P = 7 harmonics and 3 levels were used in the
ACE algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

3.16 (Top) Surface currents on the Minkowski fractal FSS. (Bottom) Power transmission spectrum for a Minowski fractal FSS. . . . . . . . . . . . . . . . . .

81

3.17 (Left) Power transmission spectrum for a 4 layer structure at normal incidence
with E0 = x using Nl = 4 and P = 7. (Right) plot of surface currents as
ˆ
viewed from the side of a single unit cell. . . . . . . . . . . . . . . . . . . .

82

3.18 (Top) N=7,552 unknown structure excited at normal incidence with E0 =
x. (Bottom) N=20,992 unknown structure excited at normal incidence with
ˆ
E 0 = x. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ˆ

83

ix

3.19 Validation of the ACE accelerated code against an analytic solution for scattering from a homogeneous dielectric slab of width 20 nm with = ±4. . . .
3.20 Validation against scattering from an electromagnetic bandgap structure solved
using FE-BI in [116]. The discretization has N = 6, 030 unknowns. Average
time to solution (without acceleration): ∼ 813 minutes. Average time to
solution (ACE acceleration): ∼ 23 minutes. Total speedup: ∼ 37×. . . . . .
3.21 Validation against scattering from an array of polystyrene spheres ( r = 2.56)
solved analytically in [117]. The discretization has N = 7, 328 unknowns.
Average time to solution (without acceleration): ∼ 2006 minutes. Average
time solution (with acceleration): ∼ 43 minutes. Total speedup: ∼ 46×. . . .
3.22 (Top) Scale geometry with N = 10, 782 unknowns. (Bottom) Calculated
reflectivity of a single model butterfly scale. Average time to solution (without
acceleration): ∼ 1430 minutes. Average time to solution (ACE acceleration):
∼ 31 minutes. Total speedup: ∼ 46×. . . . . . . . . . . . . . . . . . . . . . .

84

85

86

87

3.23 Demonstration of capability in solving a large scattering problem (N = 147, 374)
inspired by a metamaterial design presented in [122]. An average of 188 iterations per frequency was required, with an average time per iteration of 3.18
minutes, and an average total solution time of 896 minutes. . . . . . . . . . . 88
4.24 Total field/scattered field separation for imposing fields due to a dipole soucee.
The location of the source is indicated by the arrow, the red annulus is the
total field region and the blue circle is the scattered field region. At the
boundary between the two regions, the condition in Equation 4.84 is enforced. 106
4.25 Four distinct types of fragments that can occur in a tetrahedral non-conformal
mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
4.26 Quadrature at a non-conformal interface. Triangles on one side of the interface
are outlined in black, triangles on the other are outlined in red. Blue points
indicate the vertices of fragments, and green points indicate the quadrature
points used for integration. Fragments have been further decomposed into
triangles upon which a symmetric higher order quadrature is defined [130].
Some pathologies are also illustrated, as a result of finite precision arithmetic. 107
4.27 Error in E(r, t) of a TM1,1 cavity mode in the L∞ norm after 5 cycles.
Convergence is demonstrated with respect to edge length and polynomial
order. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
x

4.28 Eigenmodes of an empty 1m cubic cavity.

. . . . . . . . . . . . . . . . . . . 109

4.29 Eigenmodes of a 1m cubic cavity for which the bottom 40% is filled with a
material with εr = 2, with free space elsewhere. . . . . . . . . . . . . . . . . 110
4.30 Convergence of the relative L2 error of the fields due to a point dipole excited
by a Gaussian modulated pulse with center frequency 300 MHz and bandwidth
300 MHz at the center of a 1m spherical domain. Boundary conditions are
imposed using pure Dirichlet, Dirichlet plus total field/scattered field, and
absorbing plus total field/scattered field boundary conditions. . . . . . . . . 111
4.31 Absolute L2 error of the fields due to a point dipole excited in the same
manner as in Figure 4.30. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
5.32 (Left) A single quantum emitter above a PEC half-space. (Right) Its equivalent image problem. Arrows indicate the direction of the associated dipole
moments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
5.33 Analytical and numerical spontaneous emission rate of a two-level system
above a PEC half-space with its dipole moment parallel and perpendicular to
the interface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
5.34 Mesh of the GaAs nanocavity in which energy transfer between two 2LSs is analyzed. All exterior boundaries are terminated with a Silver-M¨ller boundary
u
condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
5.35 Images of the electric field magnitude at three points in time. All images are
taken from a plane cut halfway through the short axis of the mesh in Figure
5.34. (Top) At the peak of the exciting current. The black line indicates
the path along which the fields are sampled to generate the Green’s function.
(Middle) Half way through the simulation, the exciting current pulse is 4
orders weaker than at peak. (Bottom) Three quarters of the way through the
simulation, the exciting current pulse is 20 orders weaker than at peak. . . . 131
5.36 Distance dependence of the Fourier spectrum of the electric field component
parallel to an exciting dipole in free space. . . . . . . . . . . . . . . . . . . . 132
5.37 Distance dependence of the Fourier spectrum of the electric field component
parallel to an exciting dipole in the center of the DBR cavity. . . . . . . . . 133
5.38 Log10 of the modulus of the inter-2LS Green’s function in free space. . . . . 134
xi

5.39 Log10 of the modulus of the inter-2LS Green’s function in the cavity. . . . . 135
5.40 The inter-2LS Green’s function in the cavity for two systems with a transition
below cutoff. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
5.41 The inter-2LS Green’s function in the cavity for two systems with a transition
resonant with the cutoff. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
5.42 The inter-2LS Green’s function in the cavity for two systems with a transition
above the cutoff. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

xii

KEY TO SYMBOLS AND ABBREVIATIONS
• ACE = Accelerated Cartesian Expansions
• CFL = Courant-Friedrichs-Levy
• CG = Continuous Galerkin
• CQED = Cavity Quantum Electrodynamics
• DG = Discontinuous Galerkin
• EFIE = Electric Field Integral Equation
• EM = Electromagnetic
• FMM = Fast Multipole Method
• FSS = Frequency Selective Surface
• IE = Integral Equation
• LHS = Left-Hand Side
• PDE = Partial Differential Equation
• PEC = Perfectly Electrically Conducting
• PMC = Perfectly Magnetically Conducting
• QED = Quantum Electrodynamics
• RHS = Right-Hand Side

xiii

Chapter 1: Introduction
The interaction of light and matter is among the most enduring interests of the physical
sciences. The Lycurgus cup, a product of 4th century CE Roman glassmakers, is exemplary
of the longevity of this subject. The cup itself is made of conventional Roman glass (i.e.,
silica, soda, and lime) doped with silver and gold nanoparticles that give rise to optical
dichroism, its most distinct feature [1]. The cup is illustrated in Figure 1.1. Subsequent
artisans have used metal nanoparticles in glass coloring as evidenced by the deep reds of the
windows of the Notre Dame de Paris cathedral [2].
These aesthetic effects are a consequence of surface plasmon resonance, a phenomena in
which the electromagnetic field couples strongly to the conduction electrons near the surface
of the nanoparticles mixed into the glass. What makes this effect remarkable is that it is
geometry-dependent, in so far as the color of the light that is most strongly scattered is a
function of the size and shape of the nanoparticles [3]. What is even more remarkable is
that these structures are typically smaller than the wavelength of the visible light exciting
them, which places them severely outside of the realm of simple physical intuition. To this
end, one of the goals of nanoscience has been to construct effective theories from which new
intuition can be derived - in this case, for the manner in which light interacts with material
systems that are at or below its wavelength.
The study of light in nanoscopic systems has generated tremendous advances both in
terms of our fundamental understanding of the interaction between light and matter, (i.e.,
QED) and an astonishing range of novel technologies that exploit effects ranging from surface
plasmon resonance [3] and the formation of photonic bandgaps [4], to the Purcell effect [5]
and resonant energy transfer [6, 7]. The applications that rely upon this physics cover a
1

diverse array of fields, ranging from sensing and characterization techniques that utilize
quantum dot emission in biology [8], to single photon sources [9] and systems that generate
and preserve entanglement [10] in quantum information science.
Broadly speaking, the focus of this Dissertation is on the development of fast and accurate computational methods for the in silico analysis of light-matter interaction in the
broad range of applications discussed above. This physics typically presents a number of
challenges in terms of numerical modeling. Given the highly dispersive materials attendant
to many of these systems, and the importance of resolving subwavelength geometric features,
these systems present discretizations that have a high electrical density (i.e., there are many
unknowns per free space wavelength). To this end, fast methods for studying scattering
from these densely discretized subwavelength structures are of great relevance and will be
explored in Chapter 2.4. In some cases, the light-matter interaction cannot be modeled
using a simple set of constitutive relations, and often quantum equations of motion must be
resolved in concert with the EM fields to reach a full understanding of the system. This is
the case in cavity-mediated energy transfer physics, which will be the focus of Chapter 4.3.
The remaining Chapters are more germane to numerical modeling in and of itself, and will
contain important implementation details for those seeking to reproduce or extend the work
in this Dissertation. Throughout, the primary focus will remain on establishing accurate and
efficient methods for the solution of PDEs and IEs that arise in the modeling of light-matter
physics. Specific applications will be left to future work, particularly the problem of using
the developed tools for optimization, the natural next step for a mature piece of scientific
code.
The remainder of this Chapter will consist of an overview of some of the relevant physics
2

and applications that may be modeled using the techniques presented herein (Section 1.1), a
brief introduction to the numerical methods that will be explored for IEs and PDEs (Section
1.2), and a brief outline of the contents of this Dissertation (Section 1.3).

1.1 Physics of Interest and Technological Relevance
In this Section, an overview of light-matter phenomenology is presented, with some indication of its associated technological relevance. This is by no means comprehensive, but
should serve as inspiration for the wide range of applications that can be further enabled by
the development of computational tools for analysis, and eventually design and optimization.

1.1.1 Periodic Structures: FSS, ‘Metamaterials’, and Photonic Crystals
Periodic structures will be of interest in Chapter 3, wherein fast methods for studying
scattering of electromagnetic waves from these systems will be explored. Periodicity has long
been used as a tool in microwave engineering, with Frequency Selective Surfaces (FSS), being
one example [11]. Here, periodic arrays of subwavelength metallic patterns are distributed
across a surface to achieve some sort of resonant scattering response, wherein fields are filtered
based upon spectral content, polarization, or both. More recently, this resonant scattering
physics has been exploited in the context of ‘metamaterials’, wherein the gross scattering
response of a periodic distribution of subwavelength resonant structures is mapped onto some
effective medium with exotic properties that are either not found in nature, e.g., negative
or near-zero index of refraction [12], or not often found in nature, e.g. spatially dispersive
constitutive parameters [13]. While the theoretical utility of this approach to engineered
media has been questioned [14], the response of these systems remains interesting regardless
3

of the framework that one uses to rationalize phenomenology. In Chapter 3, scattering from
a number of exemplary FSS and metamaterial structures will be analyzed in the context of
demonstrating the utility of the ACE algorithm. One of these exemplary structures will mix
periodicity with plasmon resonance in the context of achieving negative index behavior at
optical frequencies, and will serve as a an indicator of the high electrical densities that arise
in modeling technologically interesting light-matter systems.
Photonic crystals represent another area in which periodicity is used as a means of achieving an interesting light-matter effect. Rather than periodic metallization, photonic crystal
systems consist of a periodic variation in the refractive index of some dielectric medium.
A careful engineering of the scattering resonance of the unit cell, and the Bragg scattering
resonance of the entire periodic system can lead to the formation of a photonic bandgap
in which the existence of propagating modes of the EM field is either strongly restricted
in terms of direction, or altogether forbidden [4]. This type of physics is responsible for
what is known as structural color, wherein the visible color of an object is determined by
the formation of photonic bandgaps that strongly reflect certain colors of light, rather than
chemical mechanisms in which it is determined by selective electronic absorption of light
(i.e., pigmentation) or fluorescence [15]. Photonic bandgaps arise due to naturally occurring
periodic nanostructures in the coverings of numerous animals [16], and an example in which
this is demonstrated for a butterfly’s wing [17, 18] is presented in Chapter 3.

1.1.2 Energy Transfer, CQED, and Quantum Information
Another form of light-matter interaction that will be discussed in Chapter 5 is the interaction of spatially separated systems through their mutual coupling to the EM field [6, 7].
4

The vacuum fluctuation induced dipole-dipole interaction has been of fundamental interest
in molecular physics for many years, wherein this coupling gives rise to the hierarchy of
Casimir-Polder, London, and Dispersion forces [19]. For many years, this interaction was
studied in the context of energy transfer in polymers, dyes, and biological entities [20]. More
recently, it has been studied between quantum dots [21], and quantum dots and biological
macromolecules [8], yielding a number of putative technological applications. For example,
energy transfer processes have been proposed as a means of realizing quantum logic gates
[22], and have been successfully applied to the measurement of distances in proteins [23].
While the theory for this type of interaction is well-understood in free space or disordered
environments, it becomes considerably more interesting and challenging when it occurs in
the presence of some structured matter that alters the local density of photonic modes. This
is a fundamental problem of CQED, the quantum field theory that governs the interaction
of confined light with matter [24]. Rather than considering a pure Maxwell problem with
some phenomenological treatment of matter, as will be the case in Chapter 3, the Maxwell
Equations will have to be self-consistently coupled to some other model that ‘generates’ the
constitutive relations that link light and matter. This presents new modeling challenges, in so
far as solutions to the classical Maxwell Equations must be used to inform a material model
that couples to quantum fields. A number of theoretical frameworks exist for describing
cavity-mediated interactions [25, 26, 27], but only recently have some have been employed
in a computational setting on non-canonical systems using commercial software [28, 29].
Quantum information applications provide ample motivation for studying light-matter
interaction at this level of detail. Qubits, the fundamental unit of information in quantum
computing, can be realized in terms of a number of systems that can be controlled optically,
5

ranging from excitons [30], to spins [31], to the polarization states of single photons [32].
The realization of practical technologies that actually control these systems requires the
engineering of light at the single photon level. Reliable single photon [33] and single plasmon
sources [34] have been enabled by the careful understanding of surface plasmon physics,
wherein the high energy density achieved in plasmonic modes ensures a strong light-matter
coupling. Somewhat more challenging is the generation of entangled pairs of photons. This
is typically achieved at low efficiency by way of parametric down conversion, though recent
advances have shown that the careful application of the Purcell effect can be used to engineer
the coupling between biexcitons and excitons in two-cavity systems as to yield much more
efficient generation of entangled photons [10]. In excitonic cubits, cavity-plasmon mediated
interactions have been proposed as a means of generating entanglement between spatially
separated quantum dots [35, 29].
While this is necessarily an incomplete overview of the many avenues through which
light-matter interaction promises to enable quantum information technology, it should provide sufficient evidence that this is a field that will benefit greatly from the development of
the computational methods being proposed in this Dissertation. Chapter 5 presents preliminary work to this end, with a focus on energy transfer applications, but with proposals for
ongoing/future work that are more strongly focused on quantum computation applications.

1.2 Overview of Numerical Methods for the Maxwell Equations
The accurate discretization and numerical solution of the Maxwell Equations has been
the focus of significant intellectual effort for many years. The following Sections serve to give
an overview of some of the results of this effort germane to the contents of this Dissertation,
6

as well as to establish a common notation and framework for the diverse methods that will
be discussed. Basic notations are established in Section 1.2.1. Given the breadth of methods
utilized in this work, generality is maintained with respect to establishing a notation that
is suitable for both IE- and PDE-based formulations. Similarly, Section 1.2.2 gives a brief
overview of Galerkin’s method for discretization at a level of generality that is suitable
for both formulations. In Section 1.2.3, the relative advantages and disadvantages of both
IEs and PDEs will be discussed, justifying their application to the different applications
presented in this thesis. Sections 1.2.4 and 1.2.5 will focus on details specific to IEs and
PDEs, respectively. In the former, fast methods for the solution of IEs, which take the
form of dense linear systems upon discretization, will be discussed. In the latter, differences
between Discontinuous and Continuous Galerkin discretizations will be explored, with a
focus on elaborating the advantages of DG.

1.2.1 Notation
Before discussing the details of various numerical methods, it is useful to develop the
notations that will be used throughout. Regardless of whether IE or PDE formulations are
used, much of the useful notations remain the same. To this end, the discretization and
solution of some IE or PDE of interest can be decomposed into the following steps:
1. Discretizing the domain of the problem (i.e., meshing).
2. Constructing basis functions used to represent the problem solution and data.
3. Establishing a norm for measuring the ‘correctness’ of possible solutions.
4. Defining a residual that captures the physics of the problem, and projecting the residual
7

onto a set of test functions to form a system of equations.
5. Resolving the system of equations numerically.
Notations attendant to each step of this generic process will now be developed, while details
specific to various formulations of the Maxwell problem will be discussed later.
In meshing, there will be some idealized domain Ω, with boundary ∂Ω, upon which
the continuous problem is defined. This domain will be approximated in terms of some
tessellation of non-overlapping subdomains:
Ne
Ω ≈ Ωh =

e=1

Te

(1.1)

The subdomains, Te , will also be referred to as elements. The boundary of Te is indicated
by Fe . An exemplary tetrahedral mesh is illustrated in Figure 1.2
Upon discretizing the domain, appropriate basis functions can be defined relative to
geometric primitives in Ωh . In the IE formulations that will be utilized for PEC scatterers,
linear vector basis functions will be associated with interior edges of triangular meshes. For
IE formulations for dielectric scatterers, linear vector basis functions will be associated with
all faces of tetrahedral meshes. In the PDE formulations used in this thesis, higher order
basis functions will be assigned to all tetrahedra. Both modal and nodal basis sets will be
employed, i.e., occasionally it will be convenient to work with orthogonal polynomials, and
occasionally it will be convenient to work with interpolating polynomials [36]. In all cases,
a solution to the problem will be specified as:
Nq Np (q)
f (r) ≈ fh (r) =

q=1 p=1
8

q q
cp lp (r)

(1.2)

Here, the subscript h indicates that fh (r) is an approximation of f (r) relative to Ωh . The
summation on q is over the geometric primitives (i.e., edges, faces, volumes) to which the
basis functions are tied, and the summation on p is over the Np (q) basis functions on the
q
qth primitive. The pth basis function on the qth primitive is lp (r), which is weighted by an
q
q
unknown coefficient, cp . Prior to resolving cp , however, a means for establishing how ‘good’
such an approximation is, is necessary.

A global inner product, and norm must be defined to establish the quality of an approximation. For functions defined on Ωh , f (r) and g(r), the global inner product is defined
as:
f, g Ω =
h

drf (r)g(r)

(1.3)

Ωh
For complex-valued functions, the complex conjugate of f (r) is typically used on the RHS of
Equation 1.3. For vector-valued functions, the vector inner product of the two arguments is
taken, in addition to integration. Such an inner product is useful for defining a global norm,
which is defined as:
||f || 2
= f, f Ω
L (Ω)
h

(1.4)

Similarly, it may be useful to work with an L∞ (Ω) norm:

||f ||L∞ (Ω ) = sup |f (r)|
h
r∈Ωh
9

(1.5)

Using these norms, error metrics for approximation can be constructed:

L2

= ||f − fh || 2
L (Ωh )

L∞ = ||f − fh ||L∞ (Ωh )

(1.6a)
(1.6b)

Similarly, the convergence of solutions will be validated in a discrete L2 norm.

The solution to our unspecified problem can be formulated in terms of defining a residual,
a functional that that measures the degree to which a solution fails to satisfy the IE or PDE
of interest, and forcing it to vanish in some sense. The Galerkin projection of the residual
will be forced to vanish, a process that will be described in more detail in Section 1.2.2.
Nq
This process will result in a system of N =
Np (q) equations in N unknowns, rendered
q=1
as a matrix equation for numerical representation. The entries of the attendant matrices
will require the evaluation of a number of local inner products that will now be defined.
The aforementioned concerns associated with complex- and vector-valued functions are still
assumed. For IEs, the local inner product will be taken with respect to the support of the
testing function:
f, g supp(f ) =

drf (r)g(r)

(1.7)

supp(f )
In the DG framework that will be used for PDEs, a local inner product will be defined on
each element:
f, g Te =

drf (r)g(r)

(1.8)

Te

In so far as the support of f and g will typically be Te , this isn’t necessarily distinct from
Equation 1.7; however, the notation will be more convenient in context. Additionally a local
10

inner product on the trace (boundary) of Te , Fe , will be essential:
f, g Fe =

drf (r)g(r)

(1.9)

Fe

In the DG framework, it is essential to consider inner products with multi-valued functions
on some exemplary Fe . Values are defined on both the positive and negative traces, f +
and f − , where convention dictates that the positive trace is defined on the outside of Te ,
and the negative trace is defined on the inside of Te . To quantify the behavior of multivalued functions at boundaries, the average and jump operators are defined. The average is
indicated:
{{f }} =

f− + f+
2

(1.10)

Here, f + and f − refer to the value of f on the positive and negative traces of a particular
element. The jump operator is defined as:

[[f ]] = f + − f −

(1.11)

Upon evaluating the necessary inner products, the discrete problem will be rendered in
terms of a finite linear system of equations, which takes following well-known form:

Ax = b

(1.12)

Here, A is an N × N matrix, the entries of which are generated using the inner products
defined above, x is an N vector of unknown coefficients, and b is a an N vector the entries
of which are inner products with the problem data. The mapping of both IE and PDE
11

formulations of the EM field equations onto problems of this form will be the focus of the
numerical aspects of this Dissertation.

1.2.2 Galerkin Methods
Among the most prevalent means of discretizing continuous linear problems are Galerkin
methods. These exist among the broader class of weighted residual methods, which includes
collocation and least squares methods, among others [37]. The essence of these methods
is the construction of an approximate solution in terms of basis functions that span some
space in which the regularity/physics of the problem is well-represented. A residual with
respect to the exact solution is then defined, and forced to vanish in some method-dependent
manner. In collocation methods, the residual is forced to vanish at a discrete set of points,
whereas in Galerkin methods the residual is forced to be orthogonal to the span of the basis
functions used in representing the approximate solution. It is this constraint equation on the
residual that is used to generate a linear system of equations that can be solved using matrix
methods. To give an overview of the Galerkin method that is sufficiently general to cover
both IE and PDE-based problem, an abstract problem is presented first. It is relevant to note
that the sense in which Galerkin method is being used here, is not as rigorous as would be
employed by many applied mathematicians. For example, the Method of Moments, that will
be applied to the discretization of IEs is not generally a Galerkin method. It can, however,
be realized as one, and there is little practical consequence to this distinction germane to
the heuristic discussion presented in this Section.
Consider the resolution of some unknown function, ψ, supported by a d−dimensional
open subset of Rd , Ω. The (d − 1)-dimensional boundary of Ω, denoted ∂Ω is sufficiently
12

smooth that it admits the definition of trace operators that meaningfully restrict ψ to ∂Ω.
ψ verifies the following equations:

Dψ(r) = f (r),

∀r ∈ Ω

(1.13a)

Bψ(r) = g(r),

∀r ∈ ∂Ω

(1.13b)

Here, D is a linear operator acting on some space of functions supported by Ω, and B is
some trace operator that meaningfully restricts ψ to ∂Ω. The data, f and g, as well as
the operators, D and B, are specified in sufficient detail to ensure that Equation 1.13 is
well-posed, i.e., it has a unique solution that exists and is not pathologically sensitive to the
data [38].
The first step in solving this problem with a Galerkin method is the selection of a function
space that approximates the solution of Equation 1.13, W(Ω) = {wi | supp(wi ) ⊆ Ω, i ∈
1, . . . , N }, here N is a finite number. The exact details of W(Ω) will be contingent upon the
specifics of the problem being solved. To this end, a very general list of details concerning
the choice of W(Ω) follows:
• Degree of regularity consistent with the solution of the continuous problem, i.e., if D
or B require integration/differentiation, elements of W(Ω) should be commensurately
integrable/differentiable.
• Systematic improvability, i.e., W(Ω) can be enlarged in such a way that approximation
error gets better, or at the very least does not get worse.
• Existence of convenient norms for estimating approximation error.
• Easy integration/differentiation of elements of W(Ω) (numerical or analytical)
13

Having chosen W(Ω), the approximate solution of Equation 1.13, ψh , is expressed as a linear
combination of functions in W(Ω), as illustrated previously in Equation 1.2:
N
ψh (r) =

ci wi (r)

(1.14)

i=1
Assuming that the form of the wi (r) are known, ψh (r) is practically specified by the set of
expansion coefficients, {ci |i ∈ 1, . . . , N }.
The next step of the Galerkin method is the resolution of these unknowns via the construction of a residual - a functional that quantifies the extent to which ψh fails to satisfy
Equation 1.13. Moving the data to the left hand side of Equation 1.13, a simple pointwise residual can be constructed by considering the extent to which the right hand side is
non-zero.

Dψh (r) − f (r) = Rh (r),
Bψh (r) − g(r) = Rh,∂ (r),

∀r ∈ Ω
∀r ∈ ∂Ω

(1.15a)
(1.15b)

It is apparent that ψh (r) for which Rh and Rh,∂ vanish almost everywhere over the closure
of Ω is the solution to Equation 1.13. Knowing this, however, does little to actually solve
for ψh (r), or equivalently, the set {ci |i ∈ 1, . . . , N }. To this end, Rh and Rh,∂ must be
driven to zero in some manner such that {ci |i ∈ 1, . . . , N } can be practically resolved. In so
far as N basis functions are being used to construct ψh (r), naively it should suffice to find
N unique vanishing projections of the residuals to fully determine {ci |i ∈ 1, . . . , N }. While
this approach will not, in general, result in residuals that vanish almost everywhere over
the closure of Ω, i.e., the true solution, approximate solutions computed in this manner can
14

typically be made to converge to the actual solution in some meaningful way. The rigorous
understanding of this convergence is a problem of functional analysis beyond the scope of
the heuristic outline provided in this Dissertation.
The defining feature of Galerkin methods, over other weighted residual methods, is the
specific choice of projections for which the residual vanishes. In general, the residual is forced
to vanish when projected onto the problem’s dual space, W ∗ (Ω). For a self-adjoint problem,
this is W(Ω), though it is not in general. In either case, the discrete problem is rendered:
vi , Dψh Ω − vi , f Ω = vi , Rh Ω = 0,
h
h
h

∀vi ∈ W ∗

vi , Bψh ∂Ω − vi , g ∂Ω = vi , Rh,∂ ∂Ω = 0,
h
h
h

∀vi ∈ W ∗

(1.16a)
(1.16b)

By the linearity of D and B, this can be cast in terms of a linear system of N equations in
N unknowns, of the form given in Equation 1.12. Here, the inner products in Equation 1.16
constitute elements of some matrix A some vector b, with the unknown vector x consisting
of the basis function coefficients. More specifically, for the IE problems presented in this
Dissertation, the elements of the matrix A will be the tested IE operators, with boundary
conditions implicitly enforced by the Green’s function, and the elements of the vector b will
be the tested data, i.e., the incident field. The nature of the problem will be somewhat
different for the DG problems, as they will be first rendered in a semi-discrete form. Here,
a Galerkin residual will be formed, but discretization will proceed first in space, to render
the Maxwell PDEs as a system of ODEs. Applying an explicit time marching scheme to the
system of ODEs will yield an Equation of this form, but A will be a simple block-diagonal
matrix, and b will be some complicated function of the data and boundary conditions. More
precise details concerning both will be given in Chapters 1.4 and 3.8.
15

1.2.3 Comparison of IE and PDE-Based Methods
As both IE and PDE-based methods are discussed in this Dissertation, it is worth addressing the relative advantages and disadvantages of both, to make apparent the need for
this diversity of methodology.
IE formulations of the Maxwell problem typically rely on equivalence theorems that
replace materials with equivalent currents radiating into an environment with a response that
is exactly captured by an appropriate Green’s function. Boundary conditions are implicitly
satisfied by the solution fields generated by this Green’s function. This is especially important
when considering scattering problems, in which asymptotic boundary condition must be
satisfied by scattered fields as the propagate to infinity. However, this formulation will
give rise to a dense system of equations when discretized, due to the non-local character
of the associated integral operators. For a system with N discrete degrees of freedom, IE
discretizations will give rise to matrices with O(N 2 ) non-zero entries, and will require O(N 3 )
or O(N 2 ) operations to solve directly or iteratively, respectively.
PDE formulations of the Maxwell problem typically deal with the electromagnetic fields
directly, and require the explicit treatment of fields throughout the computational domain.
While boundary conditions can be satisfied exactly on finite domains, a well-known deficiency
of PDE formulations is that radiation problems cannot be treated exactly, and require augmentation. Some of the methods for realizing proper asymptotic boundary conditions for
PDE problems include absorbing boundary conditions [39], perfectly matched layers [40], infinite elements [41], hybridization with a boundary integral [42, 43], or pole condition/Hardy
space methods [44, 45, 46]. The benefit to using PDE-based formulations of the Maxwell
problem primarily lie in the sparse nature of the linear algebra problems that they spawn.
16

For a system with N discrete degrees of freedom, PDE discretizations will give rise to matrices with O(N ) non-zero entries that can be solved rapidly using either direct or iterative
methods. Here, the cost of solution is more complex, and depends upon factors including
the dimensionality of the problem, the ordering of unknowns, and the presence or absence
of other sparse structure that may be exploited.
The relative computational cost of PDEs and IEs is a bit misleading as presented, as the
number of degrees of freedom will change dramatically between the two methods. One of
the primary benefits of IE methods is that only the scattering surfaces/volumes need be discretized. In PDE formulations, the surrounding vacuum must also be discretized, increasing
the number of unknowns. Further, IEs can also be formulated in terms of equivalent surface problems for volumetric scatterers, using the PMCHWT [47] or M¨ller [48] formalisms,
u
increasing the intrinsic disparity in the number of unknowns even further.
Time domain simulations provide another set of criteria for comparison. In PDE-based
methods, numerical dispersion becomes a significant issue, wherein propagating waves no
longer obey the proper physical dispersion relation due to the discretization process [49].
While this is not an issue in IE-based methods, there are significant issues related to stability
that are much more readily addressed in PDE-based methods. It is worth noting, however,
that there is a growing body of research on the stability of time domain IEs [50, 51].
For the purposes of this Dissertation, IE methods will be applied to the analysis of
scattering problems in Chapter 3. In particular, scattering from doubly periodic planar
structures will be of interest, wherein the use of the appropriate periodic Green’s function
yields solution currents that satisfy both Bloch-periodic conditions in plane, and the proper
asymptotic boundary conditions out of plane. However, the PDE-based DG framework will
17

be utilized for time domain simulations of cavity structures in Chapters 4 and 5. Here, the
easy realization of a fast stable solver for closed domain problems is the primary motivating
factor.

1.2.4 Fast Methods for IEs
In Section 1.2.3, the high computational cost of using IE methods was briefly addressed.
Since the 1980s, a topic of considerable research interest has been the development of fast
methods to amortize this cost. While fast direct solvers are a topic of intense current research,
the majority of the work that has been done over the past 25 years has been focused on
fast iterative solvers. To this end, the primary cost being reduced is that of matrix-vector
multiplication between the discretized IE operator and a sequence of trial solutions that
converge to the correct solution, if the problem is well-posed. Methods for accelerating the
iterative solution of IEs in this manner can be broadly split into two classes: those that effect
this through the use of FFTs, and others that exploit some manifest rank deficiency in the
IE operator.
The flagship FFT-based method is the Adaptive Integral Method due to Bleszynski [52].
Since their development in the 60s [53], it has been established that FFTs can be used to
accelerate the discrete convolution of two sequences of length N from O(N 2 ) to O(N log(N )).
The basis for AIM is the notion that the application of the discrete integral operator to a
trial solution in an iterative solver can be reduced to such a discrete convolution. However,
FFT methods can only be applied if the data lies on a uniform grid. To this end, an essential
step of the AIM algorithm is to anterpolate currents from a non-uniform mesh to a uniform
grid, and the subsequent interpolation of fields from the uniform grid back onto the non18

uniform mesh after the FFT is applied. Overall, AIM yields fast O(N log(N )) iterative
solutions. The AIM algorithm has seen success in a number of problems in computational
electromagnetics, ranging from integration into time domain codes [54], fast FE-BI codes
[55], and highly scalable parallel codes [56].

While a number of methods that exploit manifest rank deficiency exist (ACA, H-matrices,
etc), one of the primary concerns of this Dissertation will be on the ACE algorithm, which
is a variation on the FMM and its descendants [57, 58]. These methods are predicated upon
the construction of an addition theorem that permits the separation of source and testing
domains in the representation of the IE kernel. Rank deficiency arises in considering interactions between aggregate source and testing domains that are well separated in some
metric, i.e., while N sources may be tested by N observers, only O(N log(N )) or O(N )
interactions need be evaluated, rather than O(N 2 ). Intuitively, one can think of the multipole expansion of the potential due to some charge distribution. Receding away from the
distribution, progressively fewer multipoles will give appreciable contributions to the total
potential. Consequently, at a sufficiently large distance, it is more economical to represent
the potential in terms of the multipole moments of the distribution about some centroid,
rather than the individual charges comprising the distribution. The specific details of the
ACE algorithm [59], a hierarchical fast method similar to FMM, will be discussed in Chapter 3. In particular, its extension to periodic systems with long-range interactions, namely
low-frequency Helmholtz, weakly-screened Yukawa, and Coulomb potentials.
19

1.2.5 Discontinuous Galerkin vs. Continuous Galerkin
The solver that has been developed as a component of Chapters 4 and 5 is based upon
a DG framework. To this end, it is interesting to consider the relative advantages and
disadvantages of DG relative to its Continuous Galerkin (CG) counterparts.
Discretization schemes based upon a CG framework are guided by the principle that the
underlying function space should be ‘conforming’. This means, that given the function space
in which the solution to the continuous problem exists, the discretization space should be
chosen as one of its finite dimensional subspaces. Preferably, this is associated with some
framework for selecting sequentially ‘better’ subspaces, that realize a solution that is in some
sense closer to the continuous one, and this gives rise to the notion of h, p, and hp adaptivity.
This choice should also be mindful of regularity requirements, and will typically have some
means of being tied to a mesh. To this end, nodal CG discretizations rely upon the definition
of a set of interpolating polynomials defined in terms of nodes tied to mesh elements in
some regular way [60]. Within this framework, it is possible to generate a basis set that
is polynomial complete to a particular order that continuously interpolates over the entire
domain. A useful consequence of this is that the approximate solution will be single-valued
at all points, by construction. Other choices of CG basis functions include Nedelec/Whitney
elements [61, 62], that are conformal to the appropriate space of vector-valued functions in
which Maxwell fields exist.
While the basis functions in a CG framework meet many of the associated criteria, (i.e.,
regularity, improvability) there are a number of drawbacks. As CG basis functions are defined
in terms of the global geometry, the underlying mesh is required to maintain some degree
of regularity. To this end, these methods are restricted in their ability be applied to meshes
20

constructed from many different parts, with geometric primitives that do not necessarily
overlap. These meshes are deemed ‘non-conformal’, and one is illustrated in Figure 1.3
This may arise in a number of scenarios - perhaps the physics of the problem leads
to a solution with different length scales in different regions, as might be the case in an
EM problem in which two domains with very different refractive indices exist. Otherwise,
it may be the case that many different engineers designed various parts of the mesh of a
larger structure, and it is technically challenging to force them to be conformal. Otherwise,
another negative consequence of the global nature of the basis function definition arises in
time-dependent problems, wherein strong restrictions are placed on the nature of the time
integration scheme that can be utilized. Typically, one is restricted to implicit methods for
CG discretizations, which require the solution of a sparse linear algebra problem at each
time step. For higher-order methods, conditioning of the associated system worsens, and
this restriction becomes increasingly burdensome. The DG framework seeks to remedy some
of these problems.
In a DG framework the associated function spaces still typically consist of polynomial
spaces, however they are completely local, and the global functions space is described as
‘broken’. By local, it is meant that each basis function is defined with respect to exactly
one geometric primitive (typically a simplex), remaining agnostic to the global structure
of the mesh. The implication of the associated function space being broken is that it no
longer guarantees an approximate solution that is single-valued throughout the domain by
construction. It is in fact the case that the global solution will be discontinuous across
elements. To this end, these functions spaces are not conforming in the same sense as the
ones used in CG methods. To guarantee the construction of solutions with minimal numerical
21

discontinuity the concept of a ‘numerical flux’ is introduced to DG. Here, the multi-valued
nature of the solution at inter-element boundaries is introduced into the associated matrices
in such a way as to either solve some locally defined Riemann problem [63], or to penalize
the breaking of boundary conditions in a variational sense, i.e., interior penalty methods
[64].

Here, the restrictions that were previously outlined are largely lifted. Non-conformal
meshes [65], as well as non-conformal approximation order [36] are possible, thanks to the new
found global agnosticism of the basis. Similarly, in time-dependent problems, the restriction
on time integration methods is lifted as the relevant matrices are rendered block diagonal, and
explicit methods can be applied to circumvent issues of sparse linear algebra and conditioning
[66, 67]. To deal with the introduction of a CFL-like condition, DG also supports the use
of local time stepping algorithms, in which different regions of a mesh can be evolved at
different rates, to prevent multiscale discretizations from being dominated by their smallest
scales [68]. Other advantages include the ability to construct essentially arbitrary basis
sets can be constructed, particularly from non-polynomial functions that better describe the
problem physics [69]. DG also becomes a natural means of realizing domain decomposition,
and this has been used to realize adaptive local basis functions [70]. The benefits of DG
do not come without some disadvantages, though, and the element-wise definition of basis
functions means that there are more unknowns than in a CG framework [71]. However,
this has been recently combatted by mapping a CG method into the DG framework for
advection-diffusion problems [72].
22

1.3 Layout of this Dissertation
The contents of the remainder of this Dissertation are as follows:
• Chapter 2: A brief outline of the surface and volumetric EFIEs that will be used to
study scattering from periodic structures.
• Chapter 3: The use of the ACE algorithm for accelerating the computation of pairwise
potentials on periodic domains, and the solution of periodic EFIEs.
• Chapter 4: An overview of the DG framework as applied to the time domain Maxwell
Equations.
• Chapter 5: A formalism for studying field-mediated interactions between quantum
emitters using a time domain Maxwell solver, and proposals for future work.
Many of the results in Chapter 3 have been published elsewhere in [73, 74].

1.4 Figures

23

Figure 1.1: The Lycurgus cup, a piece of 1700 year old Roman nanotechnology exhibiting
optical dichroism. c Trustees of the British Museum. For interpretation of the references
to color in this and all other figures, the reader is referred to the electronic version of this
dissertation.

Figure 1.2: (Left) The domain upon which some continuous problem is defined. (Right) A
tetrahedral mesh of that domain, upon which the discrete problem will be solved. Each of
the tetrahedra into which this geometry is decomposed, is one of the subdomains, Te .
24

Figure 1.3: Illustration of a non-conformal mesh comprised of two hexahedra resolved at
different discretization rates.

25

Chapter 2: IEs for the Time-Harmonic Maxwell
Problem
2.1 Overview of IEs
In IE-based formulations of the electromagnetic field equations, the unknown quantity is
some equivalent source distribution, in terms of which potential integrals and, equivalently,
fields can be constructed. From these integral expressions of the fields in terms of the equivalent source distribution, equations can be constructed and resolved in terms of unknowns
distributed over some subset of the computational domain. Of concern in this Chapter will
be two forms of the Electric Field Integral Equation (EFIE) - one in terms of unknowns defined on a surface, for perfectly electrically conducting (PEC) targets, and another in terms
of unknowns defined over a volume, for dielectric targets.

2.2 Scattering from PEC Targets
The surface EFIE is typically used in the analysis of scattering from open PEC surfaces,
i.e., metallic screens, plates, etc [75]. This stands in contrast to its counterpart, the Magnetic
Field Integral Equation (MFIE), which can only be applied to closed surfaces [76]. When
possible, the two are often taken as a linear combination, the Combined Field Integral
Equation (CFIE) [77]. As the operators forming the two terms have disjoint null spaces,
this formulation avoids becoming ill-conditioned in the vicinity of the null space of either
taken in isolation. Still, the form of the EFIE (Fredholm first type) leads to its becoming
ill-conditioned at low frequencies - though numerous remedies exist [78, 79, 80] involving
26

Calderon identities [81]. As the PEC targets of interest for the periodic systems discussed
in Chapter 2.4 are primarily FSS, which consist of conducting strips and plates, only the
EFIE is explored in what follows. When computational results are presented, it is worth
noting that only simple measures (matrix-based preconditioning) have been taken to address
problems with conditioning.

2.2.1 Statement of the Problem
Consider a domain, Ω ⊂ R3 , in which a PEC surface, ∂ΩC , equipped with the surface
normal vector, n, is immersed. The material composition of the background in which ∂ΩC
ˆ
is immersed is homogeneous and linear, and without loss of generality, it is considered to be
vacuum characterized by the constitutive parameters ε0 and µ0 . Incident on this domain is
a plane wave of the form:
Einc (r) = E0 e−iκ0 ·r

(2.17)

Here, κ0 is the associated wavevector with κ = |κ0 |, and E0 ∈ C3 describes the incident field strength and polarization. The interrogating field gives rise to a scattered field,
Escat (r), such that the total field can be expressed as a linear superposition, E(r) =
Einc (r) + Escat (r). Using the surface equivalence principle [82, 76], this scattered field
can be expressed in terms of a convolution of an equivalent source distribution JS (r) supported by ∂ΩC , radiating into an homogeneous domain:
d2 r G(r, r ) · JS (r )

Escat (r) = −iκη
∂ΩC

27

(2.18)

Here, η =

µ0 /ε0 and G(r, r ) is a dyadic Green’s function, defined in terms of the scalar

Green’s function G(r, r ) as follows [82]:

G(r, r ) =

I+

κ2

G(r, r )

(2.19)

Here, I is the unit dyad. At this point, it is worth emphasizing that the Green’s function
in Equation 2.19 gives the response of the system at r due to a point excitation at r , in
the absence of the PEC surface, ∂ΩC . The Green’s function is responsible for satisfying
the Sommerfeld radiation boundary condition in free space, though often more complicated
boundary conditions may be desired. Of particular interest in Chapter 3 will be Green’s
functions that generate Bloch-periodic fields over the unit cell boundaries of some array,
while also satisfying radiation conditions receding away to infinity [83]. Details of G(r, r ) for
such domains are given in Appendix 5.5. What remains of this section will give consideration
to the statement of the surface EFIE with an arbitrary G(r, r ).

A Fredholm integral equation of the first type can be posed by enforcing boundary conditions on the tangential electric field over the support of the equivalent currents. This
condition can be stated physically as the vanishing of the tangential component of the electric field over ∂ΩC , or mathematically as:

n × Escat (r) = −ˆ × Einc (r),
ˆ
n

∀r ∈ ∂ΩC

n
d2 r G(r, r ) · JS (r ) = −ˆ × Einc (r). ∀r ∈ ∂ΩC

− iκηˆ ×
n

(2.20)
(2.21)

∂ΩC
TS {JS } [ r ] = −ˆ × Einc (r). ∀r ∈ ∂ΩC
n
28

(2.22)

This equality is the surface EFIE. In Section 2.2.2, the details necessary for its discretization
are presented, and then later in Section 3.7.1 it will be utilized to study scattering from
periodic arrays of PEC structures.

2.2.2 Numerical Details
Having stated the continuous problem, the remaining task is the specification of a set
of basis functions for the discrete problem. To this end, an approximate JS,h (r) will be
defined in terms of functions supported by a simplicial tesselation of ∂ΩC , i.e., a triangular
mesh of the PEC surface. An exemplary triangular mesh for a periodic scattering problem
is illustrated in Figure 2.4.
Linear vector basis functions due to Rao, Wilton, and Glisson [84] (RWG basis functions)
will be defined on each interior edge of the mesh. Some of the useful features of these functions
are as follows:
• Basis functions are associated with all interior edges of the mesh and compactly supported by the two triangles sharing such an edge.
• Continuity of currents is satisfied across all edges.
• The charge density, i.e., divergence, on the interior of the support is piecewise constant,
and globally neutral.
• The charge density is free of line charges on all boundaries of the support.
Per the list of desirable features presented in Section 1.2.2, these functions have the
appropriate degree of regularity, are easily integrated/differentiated, and have the added
benefit of being physically intuitive. If a more accurate representation of the current is
29

needed, higher-order vector basis functions [85] may employed, though they will not be
explored in this Dissertation.

Having defined an appropriate basis set, the approximate equivalent surface current is
given as:
Nq
JS (r) ≈ JS,h (r) =

cq uq (r)

(2.23)

q=1

Here, the summation is taken over all Nq interior edges of the mesh upon which the qth
RWG basis function is defined. The discrete form of the EFIE can be posed by i) computing
the matrix associated the Galerkin projection of the integral operator onto the RWG basis
and ii) computing the column vector associated with the projection of the right hand side
onto the RWG basis. The resultant discrete form is given as:

ur , TS {JS,h } supp(u ) = − ur , n × Einc supp(u )
ˆ
r
r

(2.24)

Here the equation is verified for all ur , and the unknown quantity is the set of Nq complex
coefficients cq . The non-local nature of the integral operator makes this a dense system of
2
equations that will, in general, have have Nq non-zero entries. Using matrix factorization
3
methods [86], the cost of solution will scale as O(Nq ), whereas using iterative methods [87]
2
the cost of solution will scale as O(nit Nq ), where nit is the number of iterations required
to achieve a fixed error. The focus of Section 3.7 will be on reducing the cost of iterative
solution to O(nit Nq ) using the ACE algorithm.
30

2.3 Scattering from Dielectric Targets
The volumetric formulation of the EFIE is typically applied to scattering from inhomogeneous or thin dielectric structures [88]. Frequently, it is advantageous to instead use
the Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT) [47] or M¨ller [48] formulations
u
when studying scattering from dielectrics, as they only require the representation of interfaces between regions in which the medium is piecewise homogeneous. However, they require
the doubling of unknowns at these interfaces, as both interior and exterior problems must be
posed. Further, as the M¨ller and PMCHWT formulations require the use of the in-medium
u
wavenumber in the Green’s function, whereas the volumetric EFIE does not, fast methods
reliant upon the representation of the Green’s function as some sort of truncated series may
require the retention of more terms to achieve the requisite accuracy in lossy or high index
media. As ACE is one such method and scattering from such media will be demonstrated
in Chapter 2.4, namely noble metals at optical frequencies, the volumetric EFIE will be
considered in what follows. Unlike the surface EFIE, this form of the EFIE takes the form
of a perturbation on the identity (Fredholm second type), and is typically free of issues
with conditioning in practice [89, 90, 91]. However, because it is convenient, matrix-based
preconditioning measures will still be applied when generating results.

2.3.1 Statement of the Problem
Consider a domain, Ω ⊂ R3 , in which a set of dielectric scatterers is embedded, characterized by the spatially varying dielectric function ε(r, ω). A homogeneous background of
vacuum is assumed to exist, such that ΩD = supp(ε(r, ω) − ε0 ), is compact. The domain
31

is excited by an incident planewave of the same form as Equation 2.17. The interrogating
field gives rise to a scattered field, Escat (r), such that the total field can be expressed as a
linear superposition, E(r) = Einc (r) + Escat (r). Using the volumetric equivalence principle
[76], this scattered field can be expressed in terms of convolutions with an equivalent source
distribution JV (r) supported by ΩD , radiating into an homogeneous domain:
d3 r G(r, r ) · JV (r )

Escat (r) = −iκη

(2.25)

ΩD
Here, G(r, r ) is defined in Equation 2.19 and JV (r) is related to the electric displacement,
D(r) = ε(r, ω)E(r), as follows:

JV (r) = iωCD (r)D(r),

where CD (r) =

ε(r, ω) − ε0
ε(r, ω)

(2.26)

From this substitution, it should be evident that the support of JV (r) is ΩD . Another
convenient feature of this substitution is that any discontinuities in the normal component
of JV (r) across material interfaces is due to CD (r) rather than D(r), which is normally
continuous in the absence of impressed charges.
A Fredholm integral equation of the second type can be posed by enforcing the identity
of the total field.

E(r) − Escat (r) = Einc (r),
D(r)/ε(r, ω) − iκη

∀r ∈ ΩD

(2.27)

d3 r G(r, r ) · iωCD (r )D(r ) = Einc (r),

∀r ∈ ΩD

(2.28)

ΩD
D(r)/ε(r, ω) + TV {D} [ r ] = Einc (r),
32

∀r ∈ ΩD

(2.29)

This equality is the volumetric EFIE. In Section 2.3.2, the details necessary for its numerical
solution are presented, and then later in Section 3.7.2 it will be utilized to study scattering
from periodic arrays of dielectric structures.

2.3.2 Numerical Details
The basis function space that will be utilized for the discrete problem is very similar
to the one for the surface EFIE in Section 2.2.2. While the approximate Dh (r) will be
defined in terms of functions supported by a simplicial tesselation of ΩD , this domain is of
a volumetric extent, so the associated simplices are tetrahedra. An exemplary tetrahedral
mesh is illustrated in Figure 2.5.
Linear vector basis functions will be used, due to Schaubert, Wilton, and Glisson [88]
otherwise known as SWG basis functions. Some of the useful features of these functions are
as follows:
• Basis functions are associated with all interior/exterior faces of the mesh and compactly
supported by the one or two adjacent tetrahedra.
• Continuity of the electric displacement is satisfied across all faces.
• The charge density, i.e., divergence, on the interior of the support is piecewise constant,
and globally neutral.
Similar to RWG basis functions, SWG basis functions have the appropriate regularity, are
easy to manipulate, and intuitive. As with the surface EFIE, higher-order vector basis
functions [85] are also available, but will not be necessary for the results presented in Chapter
2.4.
33

Using the SWG basis, the approximate electric displacement is rendered as:
Nq
D(r) ≈ Dh (r) =

cq uq (r)

(2.30)

q=1

Here the summation is taken over all Nq faces, both interior and exterior, upon which the qth
SWG basis functions is defined. A discrete form of the volumetric EFIE is derived similar
to the discrete form of the surface EFIE, leading to:

ur , Dh /ε supp(u ) + ur , TV {Dh } supp(u ) = ur , Einc supp(u )
r
r
r

(2.31)

This equation is verified for all ur , and the unknown quantities are the Nq complex coefficients. As with the surface EFIE, this is a dense system of equations that will generally
3
2
require O(Nq ) operations for a direct solution or O(nit Nq ) for an iterative solution. Section
3.7 will include results that demonstrate the O(nit Nq ) solution of this form of the EFIE.

2.4 Figures

34

Figure 2.4: An exemplary triangular mesh used in calculating scattering from a periodic
array of PEC split ring resonators.

35

Figure 2.5: An exemplary tetrahedral mesh used in calculating scattering from a periodic
array of dielectric spheres.

36

Chapter 3: Periodic ACE
3.1 Overview of Periodic ACE
Chapter 2 presented an overview of the IE formulations that will be used to analyzing
scattering from periodic structures. The focus of this Chapter will be fast methods for
solving them - namely, the ACE algorithm. At first, a generic approach will be taken
in which extensions of the ACE algorithm will be presented for a number of long-range
interactions in periodic systems - weakly screened Yukawa and low-frequency Helmholtz.
The attendant theoretical and implementation-oriented details will be described, including
modifications made to the octree structure - a hierarchical decmomposition of the domain
common to hierarchical fast methods. Results will be presented that validate claims of
linear scaling and convergence to arbitrary accuracy. After validating the performance of
the kernel evaluation, the ACE algorithm will be applied to the analysis of scattering from
PEC and dielectric structures, with details affirming agreement with established codes and
measurement, as well as significant acceleration over an unaccelerated reference code.

3.2 Hierarchical Fast Methods
3.2.1 Tree Codes, FMM, and ACE
In Section 1.2.4, the basic philosophy of fast methods for IEs was discussed in terms of
reducing the cost of effecting a discrete convolution between a dense N × N kernel matrix
and a vector of length N . Classically, the cost of this operation is O(N 2 ), however, fast
methods promise to reduce this cost to O(N logα (N )) for α ∈ (0, 1), a geometry-dependent
37

constant. While this was discussed in the context of the iterative solution of IEs, more
generally one might consider the evaluation of the energy of interaction between N bodies
due to a pairwise potential. Here, the unknown to be resolved is simply the convolution
itself, rather than its operand. This was in fact the motivation that drove the development
of the first tree codes, due to Barnes and Hut [92], and the FMM of Greengard and Rokhlin
[93] - the rapid evaluation of the gravitational/Coulombic potential due to a distribution of
masses/charges.
In spite of their frequent conflation, tree codes and FMMs represent two distinct approaches to the efficient evaluation of pairwise potentials that are best explained in [94].
The Barnes-Hut tree code [92] was based upon the observation that the gravitational potential at a point due to masses far away can be accurately represented in terms of a truncated
multipole expansion of the mass distribution. More generally, in tree codes the interaction
between a source-observer pair is computed using one of three options: (i) directly, (ii) at
each observation point using the multipole expansion due to a cluster of sources, or (iii) using
local expansions at a cluster of observers. FMMs take this set of principles further, introducing aggregation and disaggregation operators, that permit the computation of potentials
in an optimal manner [94].
The ACE algorithm is included under the umbrella of FMM-type methods, as it includes
both aggregation and disaggregation operators. Classical FMMs rely on constructing a representation of the Green’s function in terms of special functions; for instance, the Coulombic
Green’s function is represented in terms of spherical harmonics. The ACE algorithm, on the
other hand, is based upon a representation of the Green’s function as a generalized Taylor
expansion that is expressed in terms of totally symmetric tensors which reduce the overall
38

cost relative to other Cartesian methods.
In 2007, ACE was introduced for the evaluation of potentials of the form r−ν [59]. While
it is based on Taylor series expansions, it was shown that it is possible to develop exact aggregation and disaggregation operators. Due to this feature, no further error is introduced
in translating the origins of multipole or local expansions, in contrast to FMMs. More interestingly, using the well-known equivalence between traceless Cartesian tensors and Legendre
polynomials, it was shown that it is possible to develop relationships, both in terms of cost
and operations, between ACE and classical FMM. As ACE is not wedded to addition theorems for special functions, it is possible to apply these to the rapid evaluation of many
different potentials. To date, this has also been done for Lienard-Wiechert potentials [95],
and diffusion, Klein-Gordon, and lossy wave potentials [96]. Likewise, ACE has also been implemented together with FMM for the wideband analysis of electromagnetic phenomena [97],
with analytical error bounds that have been demonstrated via numerical experimentation.
In what follows, the extension of ACE to potentials in periodic systems will be described.

3.2.2 Specialization To Periodic Systems
While FMMs, including ACE, have been studied extensively in the context of nonperiodic problems [98, 99, 100], their adaptation to periodic problems is encountered less
frequently in the available literature. In particular, virtually all of the previous work on
periodic FMMs is for Coulomb systems, with only a few examples of papers on periodic
Helmholtz [101, 102] or Yukawa systems [103]. This can be partially attributed to difficulties
commonly encountered in periodizing long-range potentials. This can be best seen in terms
of the development of periodic Coulomb FMMs. In their seminal paper on the FMM [93],
39

Greengard and Rokhlin implement periodic boundary conditions, as well as Dirichlet and
Neumann, for the two-dimensional Coulomb potential. In this context, the influence of an infinite array of image cells must be computed within some original primitive unit cell. To this
end, periodicity is enabled by constructing a multipole expansion for the entire primitive unit
cell, acknowledging that each image cell in the periodic lattice will have the same form, and
evaluating a lattice sum to construct a periodized local expansion at the original primitive
unit cell. For periodic boundary conditions, Schmidt and Lee later extended this approach
to three-dimensions incorporating rapidly convergent Ewald summations [104, 105], whereas
subsequent work by Berman and Greengard saw the development of a renormalization-based
identity for evaluating such summations in two- and three-dimensions [106]. One particularly interesting extension due to Lambert, et al [107] utilizes the hierarchical structure of
the FMM to accelerate the aggregation of multipole expansions of periodic image cells recursively, yielding very efficient summations over the periodic lattice. The evaluation of
lattice sums becomes more challenging in the case of Yukawa or Helmholtz potentials. In
[103], for example, integral representations of the lattice sums are utilized and evaluated by
quadrature.

The majority of the previously referenced FMMs incorporate periodic boundary conditions only at the top of the tree, corresponding to a multipole expansion of the entire
unit cell. The local expansion due to the influence of the rest of the lattice, excepting the
nearfield of the unit cell, is evaluated using lattice sums. At all levels below the top, free
space translation operators are used, taking into account not only the effect of boxes lying
inside of the unit cell, but nearfield boxes outside of unit cell, as well. In this Dissertation,
a slightly different approach is followed, i.e., an addition theorem is developed for the full
40

periodic Green’s function for multipole-to-local translations. This is very similar in spirit to
what has been done using interpolatory methods [108, 109, 110], and in fast time domain
methods [111] in electromagnetics. By using a peroidized translation operator, all translations will be restricted to the interior of a single unit cell, reducing the total number of
operations required per tree traversal relative to these other methods, at the expense of requiring more complex translation operators and the evaluation of periodic Green’s functions
in the nearfield, both of which are relegated to preprocessing. Given differences between
test architectures and implementation, it is difficult to draw direct comparisons between
timings for our method and others available in the literature. However, extensive results will
be provided that demonstrate an exceptional acceleration relative to direct evaluation, and
breakeven points that are clearly very competitive with extant methods.

3.3 Statement of the Problem

3.3.1 Continuous Problem
Consider a domain, Ω ⊂ R3 , supporting a source distribution, ρ(r) : Ω → R, that gives
rise to an unknown potential, ψ(r), governed by the following Helmholtz equation:

2 + κ2 ψ(r) = −ρ(r)

(3.32)

Here, κ ∈ C, such that Equation 3.32 can be continued to a Yukawa problem for purely
imaginary κ, a lossy Helmholtz problem for general complex κ, or a Coulomb problem for
κ = 0. Assume that ρ(r) is Bloch-periodic with respect to a µ-dimensional lattice, Lµ ,
41

defined as:

µ
Lµ := {t(nµ ) =

i=1

ni ai | ni ∈ Z}

(3.33)

Here, nµ is used as a short-hand notation for the µ-tuple (n1 , . . . , nµ ), and ai are the
primitive vectors associated with some Bravais lattice. One can then define the associated
reciprocal lattice, L∗ :
µ
L∗ = {k(nµ ) =
µ

µ
i=1

ni bi | ni ∈ Z}

(3.34)

Here, bi are the primitive reciprocal lattice vectors, which together with the primitive lattice
vectors satisfy ai · bj = 2πδij .

Using these definitions, the Bloch-periodicity of ρ(r) can be expressed as follows:

ρ(r − t(nµ )) = e−iκ0 ·t(nµ ) ρ(r) : t(nµ ) ∈ Lµ

(3.35)

Here, κ0 ∈ Rµ is the Bloch wavevector describing the phase relationship between values
of ρ(r) when the argument is translated by a vector in Lµ . This condition is derived in
Section A.1.1 of the Appendix. It is evident that the specific case of κ0 = 0 corresponds to
conventional periodicity, wherein ρ(r) is left invariant under a translation of the argument
by a vector in Lµ . Given this constraint on ρ(r), it can be characterized over a reduced
volume namely the primitive cell, Ωµ ⊂ Ω, the minimal set which reflects the translational
symmetry of the lattice. One can completely reconstruct Ω from a union of primitive cells
shifted by all lattice vectors in Lµ :
Ω = {Ωµ + t(nµ ) | t(nµ ) ∈ Lµ }
42

(3.36)

Here, the set Ωµ + t(nµ ) = Ωn is referred to as the nµ th image cell, where Ω(0,...,0) = Ωµ
µ
is the central primitive cell. The measure of space occupied by a single primitive cell in the
subspace of the lattice is denoted as Aµ (i.e., A1 = length of the primitive cell, etc.). Figure
3.6 provides a pictorial representation of Ω for µ = 2.

By the translational invariance of the coefficients in Equation 3.32, that ρ(r) is Blochperiodic is sufficient to guarantee the Bloch-periodicity of ψ(r), as well.

ψ(r − t(nµ )) = e−iκ0 ·t(nµ ) ψ(r) : t(nµ ) ∈ Lµ

(3.37)

While this will suffice to constrain all boundary conditions for µ = 3, for µ = 1 and µ = 2
there are 3 − µ unspecified boundaries. These boundaries are chosen to behave as if unbounded space lies beyond them, as this is frequently the physically relevant choice for
numerous applications. More specifically, the situation in which the potential either decays
as an observer recedes away from the lattice for complex κ, or satisfies a radiation boundary
condition for real κ is under consideration.

Having fully specified the problem, i.e., PDE and boundary condition, potential solutions
for ψ(r) may be constructed.

dr Gµ (|r − r |)ρ(r )

ψ(r) =

(3.38)

Ωµ
Here, Gµ (|r − r |) is the appropriate Bloch-periodic Green’s function. Given its rapid and
absolute convergence in and away from Lµ , the Ewald representation [112] of the Bloch43

periodic Green’s function is utilized.

Gµ (|r−r |) =

t(nµ )∈Lµ

e−iκ0 ·t(nµ ) Er (|r−r − t(nµ )|)+

k(nµ )∈L∗
µ

Ek (r−r , k(nµ )+κ0 )
(3.39)

Following the usual convention, the first summation will be referred to as the ‘real sum’ and
the second as the ‘reciprocal sum’. The functional form of the terms in the real sum are
independent of µ:

Er (|r − r − t(nµ )|) =

e±iκ|r−r −t(nµ )|
κ
erfc η|r − r − t(nµ )| ± i
2η
± 8π|r − r − t(nµ )|

(3.40)

Here, ‘erfc’ is the complementary error function, and η ∈ R+ is referred to as the ‘splitting
parameter’, the value of which controls the relative rate of convergence of the two sums
in Equation 3.39. The significance of this parameter is addressed in Section A.1.2 of the
Appendix.
The form of the terms in the reciprocal sum depend on the value of µ:

Ek (r − r , k(nµ ) + κ0 ) =

∞
ei(k(n1 )+κ0 )·rl
(−1)ν
(|rt |η)2ν E1+ν
4πA1
ν!
ν=0

α2 (n1 )
4η 2

(µ = 1)
(3.41a)

=

ei(k(n2 )+κ0 )·rl
4A2 α(n2 )

±

e±α(n2 )|rt | erfc

α(n2 )
± η|rt |
2η

(µ = 2)
(3.41b)

2
2
ei(k(n3 )+κ0 )·rl −α (n3 )/4η
=
A3 α2 (n3 )

(µ = 3)
(3.41c)

44

Here, En is the exponential integral of nth order, and rl and rt are projections of r − r
parallel and transverse to the span of Lµ , respectively. The form of α(nµ ) is independent
of µ:
|k(nµ ) + κ0 |2 − κ2

α(nµ ) =

(3.42)

Given the form of Gµ (|r − r |), and some arbitrary ρ(r) on Ωµ , ψ(r) can be evaluated for
all r in Ω. An exemplary discretization of this calculation is considered next, to which the
ACE algorithm can be applied to achieve the efficient evaluation of a discrete form of ψ(r).

3.3.2 Discrete Problem

Without loss of generality, consider the case where ρ(r) can be expressed in terms of N
point sources distributed throughout Ωµ :
N
ρ(r) =
s=1

ws δ(r − rs )

(3.43)

Here, rs ∈ Ωµ is the location of the sth source, and ws is its associated weight. Equation
3.38 can be used to write ψ(r) as:
N
ψ(r) =
s=1

ws Gµ (|r − rs |)
45

(3.44)

The calculation of the mutual interaction between the sources is considered, i.e., Equation
3.38 is evaluated at all source points. This can be recast as a matrix equation as follows:

ψ(ro ) =
s=o

Gµ (|ro − rs |) → Ψo = (1 − δos )Gµ (|ro − rs |)Qs
Ψo = Gos Qs

(3.45a)

(3.45b)

Here, Gos is an N ×N matrix in which the self-terms are explicitly excluded and Qs and Ψo
are N × 1 column vectors. From this form of the potential, it is evident that the evaluation
of Ψo will require O(N 2 ) operations and O(N 2 ) storage using direct methods. It is the goal
of Sections 3.4, 3.5 and 3.6 to demonstrate that the ACE algorithm can be used to reduce
the total cost of this computation to (N ) in both time and storage. In Section 3.7 it will
be demonstrated that this same computational framework can be adapted to the iterative
solution of integral equations that arise in electromagnetic scattering problems. It is worth
remarking here that the latter Section will be concerned with the resolution of Qs , given
some known Ψo .

3.4 Details of the ACE Algorithm

ACE and other FMM-type methods achieve an O(N ) cost in timing and storage by
approximating Equation 3.45 as follows:

Ψo = Gos Qs ≈ Gnear Qs + T ACE (Qs )
os
46

(3.46)

Here, Gnear is a sparse matrix, carrying only the entries of Gos that arises due to sourceos
observer pairs which are in some metric, ‘near’, and T ACE is some composition of operations
which approximately effects the remaining ‘far’ source-observer pairs. T ACE will be defined
explicitly in Section 3.4.4. The ACE algorithm describes the metric that demarcates ‘near’
and ‘far’ interactions based upon a hierarchical decomposition of Ωµ , and provides rules for
evaluating both terms in Equation 3.46 in a manner requiring O(N ) operations and storage.

Such a hierarchical decomposition is achieved by mapping all source and observer points
onto a regular octree data structure, henceforth referred to as ‘the tree’. This structure
provides a natural metric by which ‘near’ and ‘far’ interactions can be separated, as well as
a means by which the application of T ACE can be mapped onto the traversal of the tree.
With this in mind, the ACE algorithm proceeds along the following steps:

1. Construct the tree based upon discretization of ρ(r) and ψ(r) in Ωµ .
2. Fill nearfield/farfield interaction lists based upon the structure of the tree.

3. Precompute nearfield matrix elements and ACE translation operators.

4. Compute Ψo via sparse nearfield matrix multiplication and tree traversal.
Steps 1-3 constitute pre-processing, which need be performed only once for a fixed source
distribution. Step 4 is the only recurring step, required for the repeated evaluation of a
potential in which the entries of Q vary, as may be required for the iterative resolution of a
linear system of equations. In Sections 3.4.1-3.4.4 each of these steps is described in detail.

47

3.4.1 Tree Construction
Tree construction is based upon the specification of the desired number of levels, Nl , chosen to optimize the cost and/or error. The primitive cell, Ωµ is recursively subdivided into
boxes of equal volume Nl − 1 times. A single level of the tree is defined by a set of boxes of
equal volume, with 1 being the level with the smallest (‘leaf’) boxes. At a given level, a box
subordinate to a larger box at the level above is deemed the child to the larger box, which is
deemed its parent. Every box is assigned an address in octal, based upon a canonical octree
decomposition, allowing us to readily acquire the ‘genealogy’ of a particular box given its
address alone. This method of addressing boxes is called Morton ordering or Z-space filling
curves [113]. Each box at the leaf level carries a list of the point sources/observers lying
inside its boundaries. Given the hierarchical structure of the tree, it is trivial to determine
the point sources/observers subordinate to boxes at lower levels by recursively aggregating
children until the leaf level is reached.

As will be explicitly discussed in the next section, the construction of the interaction lists
requires the consideration of not just boxes that lie inside of Ωµ , but their nearest images as
well. This is done to properly catalog near and far interactions, as is elucidated in the next
section. The images of the cell are included by adding two fictitious levels together with the
tree hierarchy that represents the nearest images of Ωµ . Figure 3.7 illustrates this modification to the overall structure. Note, these image boxes will not contain any explicit sources
or observers, and are not involved in tree traversal. Instead, they only serve as placeholders used in addressing boxes, incurring a negligible computational overhead. Consequently,
their contribution to the height of the tree is not considered when it is referenced. As the
48

addressing scheme of a regular octree follows Morton ordering, boxes inside the primitive cell
will fall within a contiguous address space [113]. This provides a simple means by which real
boxes can be distinguished from image boxes. As will be discussed in Section 3.4.2, knowledge of the nearest images of the primitive cell are essential in constructing the necessary
interaction lists, which this extension of the conventional tree structure provides.

3.4.2 Interaction Lists
Among the most crucial steps in tree-based methods is choosing an appropriate rule for
separating ‘near’ and ‘far’ interactions. In non-periodic domains, this demarcation is a relatively straightforward task. Two boxes are considered to be ‘far’ from one another if (i)
they are separated by at least one box length and (ii) their parents are not ‘far’ from one
another. The first portion of this criterion is naturally tied to the distance between boxes,
and ensures that spatial variations in the relative potential of the two domains are limited.
The latter portion, on the other hand, ensures that ‘far’ interactions are computed in O(N )
operations.

One of the primary challenges in adapting ACE to periodic domains is constructing a rule
that meets both of these needs: minimal spatial variation in the relative potential as well as
O(N ) scaling. This is somewhat non-trivial, in that the periodic boundary conditions map
the problem onto a domain with a toroidal topology. Consequently, the distance between
two boxes is no longer unique, as each box will possess images that effectively contribute to
the relative potential through the periodic Green’s function. In determining whether or not
49

two boxes in the primitive cell are ‘far’ from one another, a rule that gives consideration to
all image boxes must be constructed. Such a rule is as follows:

1. The original boxes and all of their images are separated by at least one box length.
2. Among the parents of both the original and image boxes, at least one pair is not ‘far’.
Figure 3.8 illustrates the ‘near’ and ‘far’ boxes for exemplary boxes at different levels of
a tree with Nl ≥ 4. The contents of boxes outlined in blue are stored explicitly in the tree,
whereas boxes outlined in red are part of the fictitious extra levels discussed in the previous
Section. For both exemplary boxes, it is worth noting the manner in which the nearfield
wraps around the unit cell. The source box residing in the upper left corner of the primitive
cell will participate in nearfield interactions with observer boxes residing at each of the other
3 corners, in spite of their apparent distance. Poorer error convergence has been observed
when using the usual, non-periodic rules for parsing interactions, so this is an important
subtlety to keep in mind when adapting extant codes to periodic problems.

Mechanically, the construction of the interaction lists is straightforward. For a given
level, the range of addresses which lie inside the primitive cell is known, so an iteration over
boxes can be explicitly restricted to only these boxes. This iteration begins at the (Nl − 1)th
level, just below the ‘root’ of the real tree, at which all boxes are ‘near’ each other. At the
next level down, the children of these boxes are iterated over, and the ‘near’ versus ‘far’ rule
outlined above is applied. As the nearest image boxes are addressed using Morton ordering,
both the relative location of an image box given its corresponding real box and its entire
‘genealogy’ can be located. The descent of the tree continues, iterating over only real boxes
50

at each level, until all boxes at the leaf level are visited, at which point the interaction lists
are complete. These lists are next employed in precomputation, in which they are used
to construct a list of the necessary nearfield matrix elements and unique ACE translation
operators that need be computed.

3.4.3 Precomputation

Precomputation can be broken into two stages, (i) nearfield matrix elements and (ii) ACE
translation operators. In constructing nearfield matrix elements, all unique source-observer
pairs are considered, and elements of Gnear are computed directly. Precomputation of
ACE translation operators is slightly more complex. All farfield interactions are sorted
by the distance separating source and observer parent domains in each of the Cartesian
coordinates. Given the regularity that our decomposition of the domain imposes (i.e. all
boxes at a given level have the same dimensions), and that the ACE translation operators
only depend upon the relative separation of two boxes, and not their absolute position, there
is significant degeneracy among the ACE translation operators. Consequently, in sorting by
domain separation, a subset of unique translation operators can be identified. These can
be computed once, stored, and recycled. This is particularly important in periodic ACE, as
the calculation of a single unique translation operator will possess (P + 3)(P + 2)(P + 1)/6
unique components, each of which is an Ewald-like infinite sum (see Equation 3.51).
51

3.4.4 Tree Traversal
The evaluation of T ACE (Qs ), in Equation 3.46 proceeds along the five steps presented
below. In general, these steps are common to any FMM, though the mathematical details
differ from algorithm to algorithm.
1. Charge-to-Multipole (C2M): Construct Multipole expansions from point sources
in each leaf box.
2. Multipole-to-Multipole (M2M): Aggregate Multipole expansions at higher levels
of the tree.
3. Multipole-to-Local (M2L): Translate Multipole expansions about source domains
to Local expansions about observer domains.
4. Local-to-Local (L2L): Disaggregate Local expansions at lower levels of the tree.
5. Local-to-Observer (L2O): Compute potential at observers from Local expansions
in each leaf box.
In what follows, these steps are mathematically codified in terms of Theorems, proofs of
which are available in [59]. As these Theorems are reliant upon the language of Cartesian
tensors, a brief overview of the necessary notation is provided in exposition.

A Cartesian tensor of rank n is denoted by A(n) . In general, such a tensor consists of 3n
components from C, indexed by the set {αi | i ∈ {1, . . . , n}, αi ∈ {1, 2, 3}}, where an indi(n)
(n)
vidual component is given as Aα ...αn . A totally symmetric tensor is one in which Aα ...αn
1
1
is independent of any permutation on the indices, and can be represented by (n + 1)(n + 2)/2
52

independent components. Consequently, one can index components of the totally symmetric
tensor, A(n) , as A [n1 , n2 , n3 ], where ni is the number of times that index i occurs and
n1 + n2 + n3 = n. An n-fold contraction between two tensors is denoted by ·n·, such that
C(m−n) = A(m) · n · B(n) . Finally, an n-fold product of a rank 1 tensor, r, with itself is
denoted by r(n) . More details concerning Cartesian tensors and notation can be found in
References [59] and [114].

To provide context for the ACE Theorems, an exemplary scenario is constructed in which
consideration is given to the computation of potentials observed in a domain Ωobs ⊂ Ωµ ,
due to sources in a domain Ωsrc ⊂ Ωµ . The two domains are well-separated in the sense
that they are in each others’ farfield. The linearity of the ACE framework implies that
potentials in Ωobs due to other farfield source clusters can be found using the same procedure.
Superordinate to the two domains of interest are their respective parent domains, Ωsrc ⊂
p
p
p
Ωsrc and Ωobs ⊂ Ω . The centroids of these domains are located at rsrc , rsrc , robs , and
obs
p
r . This scenario is illustrated in Figure 3.9.
obs
Using these notations the steps of the tree traversal process, as outlined in the beginning
of this Section, are effected using the following five Theorems.

Theorem 3.4.1 (Charge-to-Multipole Expansion (C2M)). The potential, ψ(ro ), at
any point ro ∈ Ωobs , due to S sources at points rs ∈ Ωsrc with strength qs (s ∈ {1, 2, . . . , S})
53

can be expressed in terms of a Cartesian multipole expansion.
S
ψ(ro ) =
M(n) =

s=1
S

ws Gµ (|ro − rs |) =

∞
n=0

M(n) · n ·

(n) G (r − r ),
µ o
src

qs
(−1)n (rs − rsrc )(n)
n!
s=1

(3.47a)

(3.47b)

Theorem 3.4.2 (Multipole-to-Multipole Expansion (M2M)). A multipole expansion
of S sources about rsrc , O(n) , can be expressed in terms of a multipole expansion about the
p
point rsrc .

M(n) =

O(n) =

n
m! p
n qs (r − rp )(n) =
(r
− rsrc )(n−m) O(m) (3.48a)
(−1)
src
s
n!
n! src
m=0 P (m,n)
s=1
S

S

qs
(−1)n (rs − rsrc )(n)
n!
s=1

(3.48b)

Here, it is important to note that Equation 3.48a is mathematically equivalent to Equation 3.47b. In other words, no additional error is incurred in shifting the origin of a multipole
expansion, using Equation 3.48a instead of Equation 3.47b. It is this detail of the ACE algorithm which leads to an error which is completely independent of the height of the tree, a
feature which is demonstrated to machine precision in [59, 97, 96].

Theorem 3.4.3 (Multipole-to-Local Translation (M2L)). For a multipole expansion
p
p
about rsrc , M(n) , a local expansion L(n) that produces the same potential in Ω
is given
obs
54

by:
∞

p
(ro − r )(n) · n · L(n)
obs
n=0
∞
1 (m−n)
p
p
(n) =
M
· (m − n) · (m) Gµ (|r
− rsrc |)
L
obs
m=n n!
ψ(ro ) =

Here, the set of tensors,

(3.49a)

(3.49b)

(m) G (|rp − rp |), where m ∈ N, is required. This set is
µ obs
src

otherwise known as the translation operator. The elements of these tensors correspond to
p
p
coefficients of a Taylor series expansion of the Green’s function about |R| = |r
− rsrc |.
obs
For a rank p tensor component of the translation operator, its individual components are
written as:

(p) G (|R|) p , p , p = ∂ px ∂ py ∂ pz G (|R|) : p + p + p = p
µ
x y z
x
y
z
x y z µ

(3.50)

As the Ewald representation is utilized for Gµ (|R|), the necessary partial derivatives are
applied to each term of Equation 3.39, yielding:

(p) G (|R|) p , p , p =
µ
x y z
t(nµ )
...+
k(nµ )

p py p
eiκ0 ·t(nµ ) ∂xx ∂y ∂z z Er (|R + t(nµ )|) . . .
p py p
∂xx ∂y ∂z z Ek (R, k(nµ ) + κ0 )

(3.51)

Relatively convenient closed-form expressions for each term of these summations are provided in Section A.2 of the Appendix to facilitate implementation.

Theorem 3.4.4 (Local-to-Local Expansion (L2L)). A Local expansion O(n) centered
55

p
about r
can be expressed in terms of a local expansion about the point robs , using:
obs
L(n) =

∞

m
p
O(m) · (m − n) · (robs − r )(m−n)
obs
m=n m − n

(3.52)

Finally, the potential at a point rα in Ωo can be computed from the local expansion
centered about ro .

Theorem 3.4.5 (Local-to-Observer (L2O)). The potential, ψ(ro ) can be expressed in
terms of L(n) , a local expansion centered about robs , using:

ψ(ro ) =

∞
n=0

L(n) · n · (ro − robs )(n)

(3.53)

3.5 Cost Analysis
The computational cost of the ACE algorithm has been previously analyzed in detail
for non-periodic problems [59]. It is worth noting that this analysis is predicated upon the
consideration of a randomly uniform distribution, and deviations from non-ideal scaling, i.e.,
O(N log(N )), will be evident for highly skewed distributions. In what follows, a brief review
of the dominant costs is provided and minor differences which arise in adapting ACE to
periodic systems are highlighted.
Consider a primitive cell in which N co-located source/observer points are randomly distributed. These points are mapped onto an Nl level tree, where the number of boxes at level
l is Bl and Bl−1 = 8Bl . The average number of unknowns per leaf box is denoted s, i.e.
Nl
N/s = B1 , and the total number of boxes at all levels
B ∼ N ∼ B1 . The total cost
s
l=1 l
56

of evaluating Equation 3.46, truncating all expansions at P th order, can be broken down
into nearfield (NF), C2M, M2M, M2L, L2L, and L2O costs, each of which is summarized
below:

1. NF: s2 operations per ‘nearfield’ neighbor × 27 ‘nearfield’ neighbors per leaf box ×
B1 leaf boxes: CN F = 27N s
3
2. C2M: s operations per component × ∼ P unique multipole components × B1 leaf
6
boxes: CC2M = N P 3
6
P6
3. M2M: ∼ 720 operations per M2M translation × ∼ B1 M2M translations: CM 2M =
N P6
s 720
P6
4. M2L: ∼ 720 operations per ‘farfield’ neighbor × ∼ 56 ‘farfield’ neighbors per box ×
P6
∼ B1 boxes: CM 2L = 56 N 720
s
P6
5. L2L: Same as M2M: CL2L = N 720
s
6. L2O: Same as C2M: CL2O = N P 3
6
This leaves us with the total cost:

Ctot = N

27s +

P3
P6
+ 58
3
720s

∼ O(N )

(3.54)

3
The multiplicative factor behind N is minimized by a density of s ≈ P . As discussed
18
in Section 3.4.1, rather than specifying a minimum box size to achieve this ideal density
of unknowns prior to constructing the tree, the number of levels, Nl , is specified. Knowing
57

that B1 = 8Nl −1 , this ideal density of unknowns is then realized for Nl = 1 + log8 ( 18N ) .
P3

Excepting an improvement made since the publication of [59], this is essentially the same
cost structure. In the original publication, the number of ‘farfield’ neighbors per box is given
as 189, here it has been reduced to 56. This is achieved by taking advantage of the exact
up/down tree traversal operations in ACE. In situations in which a box is in the ‘farfield’ of
all of the children of a given parent, the M2L translation will occur at the level of the parent.
For a given box, there will be 27 − 8 = 19 such parent level interactions, and 64 − 27 = 37
child level interactions, leading to 56 such interactions in aggregate. This change in the way
‘farfield’ interactions are treated is independent of whether ACE is being used for periodic
or non-periodic problems.

The only significant difference in cost between periodic and non-periodic ACE is in the
prefactors in the C2M and M2L costs. The 27 ‘nearfield’ neighbors per leaf box and 56
‘farfield’ neighbors per box are upper bounds, and in the case of non-periodic ACE, only
realized for boxes that have no faces touching the boundary of the computational domain.
For boxes on the boundary, however, both the number of near and far interactions will be
reduced, the net effect of which is a practically negligible change in the optimal value for s
and the overall prefactor. In periodic ACE, however, boundary boxes essentially participate
in near and far interactions as if they were on the interior of the computational domain,
given our revised rule for filling interaction lists. In this sense, periodic ACE is actually
closer to the idealized cost given above.

58

3.6 Kernel Results
In Section 3.6.1 and 3.5 results are presented validating the accuracy and linear scaling
of the ACE algorithm. In both Sections results were obtained using a kernel version of the
code that only evaluates Equation 3.46.

3.6.1 Error Convergence
In demonstrating error convergence, only the contributions to the total potential arising
due to farfield interactions are considered. ΨACE = T ACE (Qs ) is computed using the
o
ACE algorithm whereas Ψdirect = Gos Qs − Gnear Qs is computed directly. This differos
o
ence is not explicitly computed, rather nearfield pairs are simply ignored in evaluating the
potential. This is done to allay any concerns about numerical cancellation in the reference
data. The relative error in the L2 -norm is reported in all numerical experiments:
N
εf ar =

o=1

|ΨACE − Ψdirect |2
o
o
N
o=1

(3.55)

|Ψdirect |2
o

As farfield contributions will tend to be slightly smaller than those due to nearfield interactions, this is in some sense a ‘worst-case’ metric for error, and one can typically expect an
order of magnitude improvement in the error for the entire potential.
In Figure 3.10, convergence in the number of ACE harmonics is illustrated for µ = 1,
µ = 2, and µ = 3. This result was achieved using a numerical experiment with the following
conditions: 8192 point sources were distributed randomly over the unit cube, which also
59

served as the primitive unit cell, and the farfield contribution to the total potential arising
due to their mutual interaction is computed both with ACE and without ACE. The order
at which the ACE expansions were truncated was varied from P = 1 to P = 10 and the
relative error in only the farfield contribution to the potential was recorded.
The convergence is monotonic in P , and achieves an acceptable level of accuracy even
for small values of P , e.g., the error is < 1% for P ≥ 3. For P ≥ 6, the rate of convergence
decreases. This can be attributed to error in reconstructing higher order terms in the Ewald
sum. As |t(nµ )| and |k(nµ )| increase in the real and reciprocal sums, the individual terms
will exhibit a slower convergence in P . This is most intuitively accessible in terms of the
reciprocal sum, in the individual terms become increasingly oscillatory as |k(nµ )| increases.
However, as these terms will make an increasingly smaller contribution to the overall summation, the rate of convergence of the sum in aggregate will be dominated by terms with
small |t(nµ )| and |k(nµ )|. However, as the reconstruction error approaches the level at
which these contributions enter the total sum, their poorer convergence will give rise to the
observed non-uniformity in P -convergence. Otherwise, the rate of convergence is essentially
independent of the type (Helmholtz or Yukawa) as well as µ, and further results will be
restricted to µ = 2 for clarity of presentation.
In demonstrating the accuracy and applicability of ACE to a wide range of practical
problems, it is interesting to analyze convergence at a fixed value of P as κ varies. Recall
that κ is inversely proportional to the screening length for the Yukawa potential, and the
wavelength for the Helmholtz potential, motivating the introduction of a unitless length
scale, λ = 2π/(κ(Aµ )1/µ )). This scale roughly defines how rapidly the potential will vary
over a unit cell, i.e., the number of wavelengths in a unit cell for the Helmholtz potential.
60

That the ACE algorithm maintains a high degree of accuracy at a fixed value of P for a
broad range of relevant λ values for both Helmholtz and Yukawa potentials is made evident
in Figure 3.11. Here 1000 unknowns are distributed over a cube of unit volume and the
farfield error is evaluated at a fixed value of P = 6, for both the Helmholtz and Yukawa
potentials with µ = 2 and λ ranging from 0.5 to 1024. The error decreases uniformly for
λ > 0.5 from ∼ 10−1 − 10−2 at λ = 0.5 to ∼ 10−7 at λ = 1024. This is essentially
due to the fact that as λ increases, slowly varying terms become increasingly dominant in
the Green’s function. To examine the degradation in error as λ decreases, it is useful to
consider previous results and analysis for the non-periodic Helmholtz kernel, as presented in
[97]. Here, the authors allude to a decrease in the efficiency of the ACE algorithm at a fixed
accuracy, i.e. as λ becomes small relative to a box length, P must be increased to achieve the
same relative error. This type of behavior is similarly evident in the periodic Helmholtz and
Yukawa potentials. Numerical experimentation indicates that the periodic and non-periodic
ACE expansions degrade at approximately the same rate as wavelength is decreased at fixed
P . While increasing P somewhat alleviates this, as λ decreases, P will become prohibitively
large at some point in maintaining a desired level of accuracy, limiting the viability of ACE.

To this end, it is worthwhile to examine this apparent limitation in the context of the
physical relevance of small values of λ. Consider λ = 0.5 (κ = 4π), the value for which
error is highest in Figure 3.11. For the Yukawa potential, λ = 0.5 corresponds to a strongly
screened system, in as far as second nearest neighbor cells will introduce a correction to the
potential less than 0.001% of the correction due to the nearest neighbor cells. In such a
strongly screened system, the interactions are better handled by short-range acceleration,
61

that may be as simple as maintaining interaction lists for a truncated potential. Physically
relevant values of λ will instead be on the order of or greater than a primitive cell length
1/µ
(Aµ ), i.e. a weakly screened system for which the ACE algorithm is demonstrably accurate. For the Helmholtz potential, rather than strong screening, small values of λ correspond
to increasingly oscillatory behavior. Fortunately, most physically relevant conditions for the
periodic Helmholtz potential involve values of λ which correspond to systems in which the
period of oscillation is on the order of, or greater than a primitive cell length. Typically
frequency selective structures, metamaterials, and photonic crystals, to name but a few applications, consist of subwavelength unit cells (i.e. λ ≥ 1 on our scale), to which the ACE
algorithm is also very well-suited.

3.6.2 Computational Cost
In Figures 3.12 and 3.13 scaling of the precomputation and tree traversal times with the
number of point sources, N , are presented. Timings are presented for 5 different values of
P to demonstrate that linear scaling can be achieved at arbitrary precision. In generating
these results, N sources were distributed randomly over a unit cube on a µ = 2 lattice with
|a1 | = |a2 | = 1 and κ = 2π. As the cost has a very weak µ dependence, only µ = 2 is
presented as it will be of interest for applications, and the scaling is essentially the same
for µ = 1 and µ = 3. The number of levels at each value of N is increased, from 3 to 6,
maintaining an average density of unknowns per leaf box of 80. The precomputation times
in Figure 3.12 include the time required for tree construction as well as the calculation of
all unique translation operators that will be later applied during tree traversal. This is a
one time cost that is presented to demonstrate two things: i.) precomputation scales weakly
62

with the number of source/number of levels and ii.) the scaling is independent of P .
The tree traversal timings in Figure 3.13 include the time required for all 5 steps of
the tree traversal process. These results are intended to demonstrate that linear scaling is
achieved independent of P . A linear regression indicates an N 1.01 scaling for all values of
P , except for P = 1, in which the scaling is N 1.09 .

3.7 ACE for Scattering from Periodic PEC and Dielectric Targets
Throughout the remainder of this Chapter, the integration of the periodized ACE algorithm into both surface and volumetric EFIE solvers will be presented. Up until now, all
of the results presented are for simple kernel tests in which scalar convolutions of the form
given in Equation 3.38 are accelerated. Two fundamental concerns must be addressed in
adapting what has been presented up to this point to the EFIE:
• The potential integrals in the EFIE are vector quantities.
• In Equation 3.38, the unknown quantity is the result of the convolution. In an integral
equation, the unknown quantity is the operand upon which the convolution is acting.
Each of these concerns will be addressed in sequence.
In the context of an EFIE solver, the role of the potential is played by a vector quantity,
namely:
dr G(r, r ) · J(r )

Escat (r) = −iκη

(3.56)

As is evident in Equations 2.22 and 2.29, the Galerkin projection of convolutions of this
form are essential components of both EFIEs. As the source (J) and observed (Escat ) fields
are vector quantities, three separate trees are necessary, to accommodate each of the vector
63

components. Otherwise, as the Green’s function is a dyadic quantity, one might think a
‘dyadic M2L translation operator’ necessary. In practice, however, this is dealt with in one
of two ways. As the component of G coefficient to I is diagonal, it can be disregarded as
being potentially problematic. The only modification that is necessary is due to the term
that is coefficient to

, i.e., the off-diagonal contribution. To effect this portion of the

Green’s function, either the derivative terms can be transferred to the source and testing
basis functions, yielding a four tree scheme, or the derivatives can be applied to the Local
Expansion polynomials. The former will be more computationally intensive, whereas the
latter will reduce the accuracy of the ACE expansion, by effectively removing two polynomial
orders from the Local expansion.
Turning to the issue of the unknown quantity being the operand of the convolution, rather
than the result, it suffices to again consider the scalar convolution in Equation 3.45. Iterative
solvers based upon Krylov subspace methods [87] can be used to solve for Q. Such a solver
generates a minimal sequence of vectors, {Qi |i ∈ 1, . . . , nit }, in the range of G from which an
n
approximate solution to Equation 3.45 can be constructed, such that ||Gos Qs it − Ψo || ≤ ε,

where ε is some designated tolerance for error. The most costly part of the iterative solution
process is the matrix-vector product, through which test solutions are generated and tested,
leading to an O(nit N 2 ) cost. To this end, the ACE algorithm is used to accelerate the
necessary O(N 2 ) matrix-vector products to O(N ), yielding an O(nit N ) overall cost.
In the following two subsections, results for an ACE accelerated EFIE solver are presented. In so far as the problem statements given in Sections 2.2.1 and 2.3.1 were agnostic
to the choice of Green’s function, it is important to note that all of the results are presented
for Bloch-periodic Helmholtz problems on a square lattice with µ = 2. Section 2.2 uses a
64

surface formulation of the EFIE for the resolution of scattering from PEC targets, as stated
in 2.2. Section 2.3 uses a volumetric formulation for the resolution of scattering dielectric
targets, as stated in 2.3. In both Sections, the focus is on first demonstrating accuracy in
reproducing either data from the literature or a reference code, and then on demonstrating
that a considerable speedup is achieved using the ACE algorithm.

3.7.1 PEC Targets
As validation of the ACE accelerated surface EFIE code for PEC targets, power transmission spectra of a number of interesting structures are presented. In all cases, the transmission
coefficient is given for the (0,0) Floquet mode on a dB scale. When making comparisons
between the unaccelerated reference code and the ACE accelerated code, the same iterative
methods, preconditioners, and convergence criteria are utilized and all infinite sums are evaluated to a residual relative error ≤ 0.1%. With the exception of the Minkowski fractal, which
is compared to a structure with dimensions found in [111], all geometries are specified in
normalized units, wherein each cell is taken to be a square of area unity, and all frequencies
are scaled such that the unit frequency corresponds to unity wavelength.
In Figure 3.14, the accuracy of the ACE accelerated solver is established by comparing
it to an unaccelerated reference code converged to a low residual error (≤ 0.001%). The
geometry under consideration is a Jerusalem cross with N = 85 unknowns, which fits inside of
a 0.87×0.87 square in our normalized units. Across all frequencies, the ACE accelerated code
and the reference code are in excellent agreement, and the absolute error in the transmission
coefficient is strictly less than 1% (typically < 0.1%). In general, even for such a low residual,
the ACE algorithm and the unaccelerated reference code require the same number of TFQMR
65

iterations.
Next, a result is provided that gives an indication of low breakeven point that can be
achieved with an efficient implementation. In Figure 3.15, the ACE code is compared to the
unaccelerated code for an Omega split ring structure with N = 356 unknowns. This structure
fits inside of a 0.95 × 0.8 rectangle, with the longer dimension parallel to the split, and a
gap of 0.11 units between layers. In this test, a slightly more permissive residual is achieved
(≤ 0.01%), and the observables remain in good agreement. This particular geometry was
chosen to show the conditions under which the total time to solution for both the ACE and
unaccelerated codes are commensurate, in this case an average of 30 seconds per solution.
The next example (Figure 3.16) demonstrates a considerable acceleration in the time to
solution for a larger problem - a Minkowski fractal FSS with dimensions taken from [111] with
N = 2, 148 unknowns. Here, GMRES is used with a residual error ≤ 1%. For frequencies
above 1 GHz, a diagonal preconditioner is utilized, whereas an ILU preconditioner is utilized
below. Using an unaccelerated reference code, the average time to solution is 4, 440 seconds
(74 minutes) per frequency. For the ACE accelerated code, the average time to solution is
108 seconds (1.8 minutes) per frequency. From this results, it is evident that a speedup of
40x can be achieved for a structure that is relatively small by the standards of free space
fast methods, while retaining accuracy.
In Figure 3.17, accuracy is validated for a layered structure with N = 3, 206 unknowns.
This geometry is inspired by a Minkowski fractal structure, and consists of 2 pairs of complementary patterns, the largest of which is bounded by a square of dimension 0.95 × 0.95,
for a total of 4 layers separated by 0.24 units in the z direction. Here, TFQMR with an ILU
ˆ
preconditioner is utilized to achieve a residual ≤ 0.1%. The average time to solution for the
66

ACE code is 314 seconds per frequency (∼ 5.2 minutes), relative to 6, 530 seconds (∼ 109
minutes) per frequency for the unaccelerated code.
In Figure 3.18, surface currents for a number of structures with exceptionally high electrical densities (N=7,552 and N=20,992) are shown. Both structures have transverse dimensions that fit inside of a square of dimension 0.95 × 0.95 and were generated via fractalization
of a simple primitive structure. The first structure is based upon a triply iterated cross aperture, whereas the second structure is based upon a triply iterated Sierpinski carpet. While
the first structure consists of a single layer, the second consists of 2 complementary pairs
for a total of 4 layers separated by 0.16 units in the z direction. For the first structure, the
ˆ
average time to solution is 1, 220 seconds (∼ 20 minutes) using ILU preconditioned TFQMR
with a tolerance ≤ 0.1% for P = 7 harmonics and Nl = 5 levels. For the second structure,
it is 11, 200 seconds (∼ 3.1 hours) using diagonal preconditioned GMRES with a tolerance
≤ 1% for P = 7 harmonics and Nl = 6 levels. Extrapolation based upon the unaccelerated
code yields average solution times of ∼ 11.8 and ∼ 91.3 hours, respectively, without ACE.

3.7.2 Dielectric Targets
As with the surface EFIE solver in the previous section, validation is based upon the
reproduction of power transmission/reflection spectra for the volumetric EFIE solver for
dielectric targets. All results in this section were generated by a diagonal-preconditioned
TFQMR solver with a tolerance of ≤ 0.1%. All periodic Green’s function evaluations, in
both the Ewald and Floquet representations, are converged to an accuracy of ≤ 0.001%. All
ACE translation operators were evaluated to a similar accuracy. Unless otherwise indicated,
a tree with Nl = 5 and P = 7 ACE harmonics are utilized, to guarantee an error of ∼ 0.01%
67

in reconstructing potentials.
The first result is given in Figure 3.19. In it, the accuracy of the ACE accelerated code is
established in reproducing power transmission through an infinite dielectric slab of width 20
nm as the angle of incidence and polarization is varied at a fixed wavelength, λ = 400 nm.
A similar result is presented in [115], wherein a layered-medium formulation is compared
against an analytic calculation. Here, a periodic volumetric formulation is utilized, wherein
the unit cell has an edge length of 80 nm. As is evident, our ACE accelerated code reliably
reproduces the analytic solution at both TE and TM polarization, at arbitrarily shallow
incidence, for both positive and negative relative permittivities.
The results in Figure 3.20 and Figure 3.21 are intended to illustrate i) agreement of the
ACE accelerated VIE code with established results from the literature and ii) acceleration
relative to an unaccelerated reference code. The structure in Figure 3.20 is an electromagnetic
bandgap structure with applications in microwave engineering. The reference data is taken
from [116], wherein a hybrid finite element-boundary integral formulation is utilized. The
ACE accelerated code is in good agreement with the established result, and a 37× speedup
relative to the unaccelerated code is achieved. The result in Figure 3.21 is taken from [117]
wherein an analytic calculation is performed using Mie theory. Good agreement is again
found between ACE and the established result and a 46× speedup relative to our reference
code is achieved. This structure has been used elsewhere in the literature as validation of a
periodic FMM code, wherein the authors report 4−18 minutes per solution on an 8 processor
platform [118]. As a surface integral equation formulation is utilized, as well as parallelism,
it is difficult to draw a direct comparison between our results. Taking these differences into
consideration, it is safe to at least claim that the ACE accelerated code is competitive with
68

the state of the art.

The next result is intended to demonstrate the utility of the ACE algorithm in studying
models for novel experimental applications. It has been understood for many years that a
number of the optical effects common in the wings of butterflies, such as their deep coloring
or iridescence, arise due to periodic nanostructures that occur naturally in the scales covering
these wings, on top of any chemical coloring (i.e. pigment) that may exist [17, 119, 120].
The blue coloring of the Morpho butterfly is partially due to a photonic bandgap (PBG) in
the blue portion of the visible spectrum. This PBG is supported by a periodic ribbing in the
material covering the butterfly’s wings on the scale of visible wavelengths [17]. Previously,
Huang, et al.[18] have performed an extensive analysis of the optical properties of not just
naturally occurring wings, but coated and synthetic wings made from Al2 O3 . A simple model
that replicates the blue structural coloring of these types of structures has been constructed,
albeit one that is both thinner and oriented differently than the one in [18]. The reflectance
of this structure are given in Figure 3.22. Here, a single unit cell with |a1 | = |a2 | = 320
nm is outlined in red. The height of the structure out of plane is 280 nm and the diameter
of the larger cylinders is 130 nm, with a center-center spacing of 160 nm between nearest
neighbors. The smallest cylinders of diameter 20 nm with axes in the plane of periodicity are
oriented along the polarization vector of the incident field, with a center-center spacing of
70 nm out of the plane of periodicity. Dielectric constants for Al2 O3 were interpolated from
[121] and the resultant mesh has N = 10, 782 unknowns. It is worth noting that the unit
cell is artificially enlarged to ensure that it is cubic. This both simplifies tree construction
and increases the number of degrees of freedom for demonstrative purpose, and is not, in
principal, necessary. In this result, the unaccelerated time to solution was extrapolated based
69

upon the nearfield matrix fill time for ACE alone, neglecting the time for iterative solution,
so the speedup factor of 46× is actually representative of a lower bound on the acceleration.
The final result, in Figure 3.23, is a demonstration of capability in solving a very densely
discretized structure with a highly dispersive dielectric response. This structure is inspired
by [122], wherein a yellow-light fishnet metamaterial made of layered Al2 O3 and silver on
an ITO substrate is simulated and fabricated. Here, a structure with the same periodicity
and nanostrip widths but without tapering and an unmodified dielectric function for Ag is
modeled. The well-established Johnson and Christy permittivity for Ag [123] is utilized here,
while the authors of [122] have phenomenologically increased the collision rate to account for
surface roughness, grain boundaries, and size effects. Thusly, a direct comparison of results
is not possible. The resultant mesh has N = 147, 374 unknowns, and was solved using ACE
with 6 levels and P = 5. TFQMR was applied with diagonal preconditioning to arrive at an
error of < 3% with an average of 188 iterations per frequency at 20 frequencies. The average
time per iteration was 3.18 minutes, and the average total time to solution was 896 minutes.
Extrapolating the nearfield matrix fill time indicates that the matrix fill process alone would
require 256, 500 minutes using the unaccelerated reference code with the same parameters,
yielding a lower bound on the speedup of ∼ 286×.

3.8 Figures

70

Ω(−1,1)

Ω(−1,0)

Ω(0,1)

a2

Ω2

Ω(1,1)

Ω(1,0)

a1
Ω(−1,−1)

Ω(0,−1)

Ω(1,−1)

Figure 3.6: An exemplary periodic domain for µ = 2. The central primitive cell, as well as
its 8 closest image cells are illustrated. The support of some arbitrary ρ(r) is represented by
the blue shaded areas.

71

43

44

47

48

59

60

63

64

41

42

45

46

57

58

61

62

35

36

39

40

51

52

55

56

33

34

37

38

49

50

53

54

11

12

15

16

27

28

31

32

9

10

13

14

25

26

29

30

3

4

7

8

19

20

23

24

1

2

5

6

17

18

21

22

Figure 3.7: The primitive unit cell (blue) and the 15 implicit image cells (red) involved in the
enlarged tree. The original tree has Nl = 2, whereas the enlarged tree has Nl = 4. Rather
than octal addresses, decimal addresses are given to provide a simpler view of the ordering.

72

Figure 3.8: ‘Near’ (blue) and ‘far’ (light blue) interactions for source boxes (dark blue) at
different levels of the tree. The contents of the red image cells are not actually stored, but
need to be considered in constructing interaction lists.

73

p
rsrc
rsrc

robs
p
robs
Figure 3.9: Illustration of the exemplary scenario described above. Only the primitive cell is
represented. The leaf level domains are highlighted in dark blue, whereas the parent domains
are highlighted in lighter blue.

74

Error Convergence, Random Distribution Inside Unit Cube, κ=π
-1

Relative Farfield Error, L2 Norm

10

Helmholtz, µ=1
Yukawa, µ=1
Helmholtz, µ=2
Yukawa, µ=2
Helmholtz, µ=3
Yukawa, µ=3

10-2

10-3

10-4

10-5

10-6

1

2

3

4
5
6
7
Number of ACE Harmonics (P)

8

9

10

Figure 3.10: Convergence of the relative error in reconstructing the farfield contribution to
the potential for both Helmholtz and Yukawa systems with µ = 1, µ = 2, and µ = 3, with
κ = π. 8192 sources are distributed over the unit cube, which is also the primitive cell for
all tests.

75

-1

Helmholtz
Yukawa

Log10 Farfield Error

-2
-3
-4
-5
-6
-7
-8
10-1

100

101

102

103

104

λ
Figure 3.11: Farfield error as a function of λ = 2π/κ for 1000 points randomly distributed
over the unit cube at a fixed ACE order, P = 6, and µ = 2.

76

2

Precomputation Time (seconds)

10

101

P=1
P=3
P=5
P=7
P=9

100

-1

10

103

104
105
106
Number of Point Sources

107

Figure 3.12: Demonstration of sublinear scaling of the precomputation time for ACE expansions of different orders. These results are generated for a Helmholtz problem with µ = 2.

77

4

Tree Traversal Time (seconds)

10

3

10

2

10

P=1
P=3
P=5
P=7
P=9

1

10

100
10-1
-2

10

-3

10

-4

10

103

104
105
106
Number of Point Sources

107

Figure 3.13: Demonstration of linear scaling of the tree traversal timing for ACE expansions
of different orders. These results are generated for a Helmholtz problem with µ = 2.

78

Figure 3.14: (Top) Plot of surface currents on Jerusalem cross.√
(Bottom) Power transmission
◦ , φ = 0◦ ) with E = 2/2 (ˆ − z ). TFQMR with a
spectrum at oblique incidence (θ = 45
x ˆ
0
diagonal preconditioner was utilized to achieve a residual ≤ 0.001%. P = 7 harmonics and
Nl = 3 levels were used in the ACE algorithm.
79

Figure 3.15: (Top) Plot of surface currents on an Omega split ring structure. (Bottom)
Power transmission spectrum for an Omega SRR at normal incidence (θ = 0◦ , φ = 0◦ ) for
E0 = x. GMRES with a diagonal preconditioner was utilized to achieve a residual ≤ 0.01%.
ˆ
P = 7 harmonics and 3 levels were used in the ACE algorithm.
80

Figure 3.16: (Top) Surface currents on the Minkowski fractal FSS. (Bottom) Power transmission spectrum for a Minowski fractal FSS.
81

Figure 3.17: (Left) Power transmission spectrum for a 4 layer structure at normal incidence
with E0 = x using Nl = 4 and P = 7. (Right) plot of surface currents as viewed from the
ˆ
side of a single unit cell.
82

Figure 3.18: (Top) N=7,552 unknown structure excited at normal incidence with E0 = x.
ˆ
(Bottom) N=20,992 unknown structure excited at normal incidence with E0 = x.
ˆ

83

Figure 3.19: Validation of the ACE accelerated code against an analytic solution for scattering from a homogeneous dielectric slab of width 20 nm with = ±4.

84

Figure 3.20: Validation against scattering from an electromagnetic bandgap structure solved
using FE-BI in [116]. The discretization has N = 6, 030 unknowns. Average time to solution
(without acceleration): ∼ 813 minutes. Average time to solution (ACE acceleration): ∼ 23
minutes. Total speedup: ∼ 37×.

85

Figure 3.21: Validation against scattering from an array of polystyrene spheres ( r = 2.56)
solved analytically in [117]. The discretization has N = 7, 328 unknowns. Average time to
solution (without acceleration): ∼ 2006 minutes. Average time solution (with acceleration):
∼ 43 minutes. Total speedup: ∼ 46×.
.

86

Figure 3.22: (Top) Scale geometry with N = 10, 782 unknowns. (Bottom) Calculated reflectivity of a single model butterfly scale. Average time to solution (without acceleration):
∼ 1430 minutes. Average time to solution (ACE acceleration): ∼ 31 minutes. Total speedup:
∼ 46×.

87

Figure 3.23: Demonstration of capability in solving a large scattering problem (N = 147, 374)
inspired by a metamaterial design presented in [122]. An average of 188 iterations per
frequency was required, with an average time per iteration of 3.18 minutes, and an average
total solution time of 896 minutes.
88

Chapter 4: DG for the Transient Maxwell Problem

4.1 Overview of Time Domain DG

The primary focus of this Chapter is the solution of the time domain Maxwell Equations.
This is in contrast to the previous Chapters on IEs in a number of ways. Whereas the
unknowns of interest were previously equivalent currents that give rise to fields that satisfy
the Maxwell Equations, here the unknowns of interest are the fields themselves. Further,
while time harmonic excitations were previously considered, now transient excitations are
considered that permit the extraction of observables over a band of frequencies from a single
simulation.

4.2 Time Domain Maxwell Equations

4.2.1 Statement of the Problem

Consider a domain, Ω ⊂ R3 , in which some distribution of non-dispersive matter characterized by isotropic constitutive parameters ε(r) and µ(r) is immersed. The EM fields will
89

satisfy the Maxwell Equations:

∂ E(r, t)
= × H(r, t) + J(r, t)
∂t
∂ H(r, t)
µ(r)
= − × E(r, t)
∂t
ε(r)

(4.57a)
(4.57b)

· D(r, t) = ρ(r, t)

(4.57c)

· B(r, t) = 0

(4.57d)

Additionally, the source charge and current satisfy the continuity equation:

∂ρ(r, t)
+
∂t

· J(r, t) = 0

(4.58)

It is assumed that the sources are quiescent and autonomous, as their time evolution does
not need to be solved for, but is instead impressed at t = 0. Of primary concern will be the
resolution of E(r, t) and H(r, t) from Equations 4.57b and 4.57a, assuming that the initial
conditions, E(r, 0) and H(r, 0) are given and satisfy Equations 4.57d and 4.57c.

Considerable economy can be realized by introducing the conservation form of Equations
4.57b and 4.57a. To this end, E(r, t) and H(r, t) can be united into a 6-vector q(r, t) that
obeys the following:
Q(r)

∂q(r, t)
=−
∂t

· F(q(r, t)) + S(r, t)

(4.59)

The constitutive parameters are contained in the 6 × 6 diagonal matrix, Q(r):




0 
 ε(r)I
Q(r) = 

0
µ(r)I
90

(4.60)

The conserved quantities (i.e., the fields) form a 6-vector, q(r, t):




 E(r, t) 
q(r, t) = 

H(r, t)

(4.61)

These conserved quantities are tied to a tensor current:




u
 −ˆi × H(r, t) 
F(q(r, t)) = 

ui × E(r, t)
ˆ

(4.62)

where ui is the ith Cartesian unit vector, as well as a 6-vector source term:
ˆ


S(r, t) = 


0
J(r, t)




(4.63)

The remainder of this Chapter will focus on the discretization of this form of the Maxwell
Equations within the DG framework (Section 4.2.2) and results that validate its proper
implementation (Section 4.2.3) germane to the applications being explored in Chapter 5.

4.2.2 Numerical Details
4.2.2.1 Modal and Nodal DG Basis Functions
The basis functions utilized in discretizing Equation 4.59 are a set of globally discontinuous higher-order polynomials. To this end, both Eh (r, t) and Hh (r, t) are defined in terms
of a collection of local expansions restricted to each element of a simplicial tessellation of
Ne
Ω ≈ Ωh =
Te . These expansions may be of either a ‘nodal’ or ‘modal’ character. While
e=1
91

these expansions are equivalent by construction, each will have attendant advantages that
will be utilized in the discretization process.

The nodal expansion is of the form:
Ne Np (e)
E(r, t) ≈ Eh (r, t) =

H(r, t) ≈ Hh (r, t) =

e=1 p=1
Ne Np (e)

e
Eh (re , t)lp (r)
p

(4.64)

e
Hh (re , t)lp (r)
p

(4.65)

e=1 p=1

Here, the sum on e is taken over each of the Ne tetrahedra, upon which a set of Np (e)
e
interpolatory polynomials, lp (r), are defined. These polynomials interpolate on a set of
Np (e) points defined by the global coordinates re ⊂ Te . Each of the Ne local expansions
p
is complete to polynomial order Π(e), with the following relationship between order and
cardinality:
Np (e) =

(Π(e) + 1)(Π(e) + 2)(Π(e) + 3)
6

(4.66)

In the nodal framework, the unknown will be the 3-vector-valued functions of time Eh (re , t)
p
and Hh (re , t), i.e., the solution fields evaluated at the interpolation nodes. The specific
p
interpolating nodes that will be used in practice minimize the Lebesgue constant associated
with the interpolation, and are derived in [36] - the Lebesgue constant being a prefactor that
describes the distance a given interpolant is from the best possible polynomial approximation.
92

The modal expansion is of the form:
Ne Np (e)
Eh (r, t) =
e=1 p=1
Ne Np (e)
Hh (r, t) =
e=1 p=1

e
ce (t)ψp (r)
E,p

(4.67)

e
ce (t)ψp (r)
H,p

(4.68)

Again, the sum on e is taken over all tetrahedra, upon which a set of Np (e) modal polynoe
mials, ψp (r), are defined. These polynomials are element-wise orthogonal, such that:

e e
ψp , ψq Te ∼ δp,q

(4.69)

Each of the local expansions is complete to polynomial order Π(e) and obeys the ordercardinality relationship in Equation 4.66. The unknowns, ce (t) and ce (t) are 3-vectorE,p
H,p
valued functions of time that are the local expansion coefficients for Eh (r, t) and Hh (r, t)
restricted to tetrahedron e. The specific modal polynomials that are utilized are a minimal
set of products of Jacobi polynomials derived in [36].
It is useful to consider why these two equivalent expansions are necessary. In what follows, the modal basis will be used whenever elementary calculus is required, e.g., integration
or differentiation. It is possible to perform these operations on the modal polynomials as
they have a simple closed form. For simplicial elements, there is in fact no closed form
expression for the nodal interpolating polynomials, so these operations cannot be applied
directly. However, as will be discussed in Section 4.2.2.2, the off-diagonal matrix elements
will consist of surface integrals over the faces of the Ne tetrahedra. In the modal basis,
these integrals will generally give rise to dense off-diagonal blocks. The nodal polynomials,
93

however, have the nice property that these blocks will be sparse due to the fact that interpolating polynomials associated with nodes that do not lie on a given face, are identically
zero on that face. To this end, the general philosophy will be to compute the unknowns
in the nodal basis, but to evaluate matrix elements in the modal basis and then transform
them into the nodal basis. The form of the requisite elements will be given below in Section
4.2.2.2, whereas explicit details concerning their evaluation are contained in [36].

4.2.2.2 DG Discretization
The discretization process in DG is somewhat more involved than it was with the IEs
of Chapter 2. As before, however, it begins with the projection of Equation 4.59 onto an
e
arbitrary nodal testing function, lp (r):
e ∂q
e
lp , Q h Te = − lp ,
∂t

e
· F(qh ) Te + lp , S Te

(4.70)

The key ‘trick’ in the DG discretization scheme is to apply integration by parts to the
divergence term twice. The first integration by parts yields:
e ∂q
e
lp , Q h Te = − lp , n · F∗ (qh ) Fe +
∂t

e
e
lp , F(qh ) Te + lp , S Te

(4.71)

Here, F∗ (qh ) is called the ‘numerical flux’. It takes the place of the usual boundary term
arising from integration by parts, which is ill-defined, given the element-wise definition of
basis functions. In other words, because the basis functions do not explicitly constrain
continuity across faces, the fields at the boundary are inherently multi-valued, making the
definition of a boundary term in the integration by parts ambiguous. A second integration
94

by parts yields the following:
e ∂q
e
lp , Q h Te = − lp ,
∂t

e ˆ
e
· F(qh ) Te + lp , n · F(qh ) − F∗ (qh ) Fe + lp , S Te

(4.72)

This returns the problem to a strong form, and provides insight into how F∗ (qh ) should be
defined. As the boundary term in Equation 4.72 should vanish for solutions to the Maxwell
equations, the numerical flux should be designed in such a way as to numerically enforce
inter-element boundary conditions, as they are not explicitly built into the function space.
To this end it is useful for F∗ (qh ) to take the form of an upwind flux, first derived for Finite
Volume Methods [124, 67]:


∗ (q ) = 

n · F(qh ) − F

h



n×(Z + Hh −n× Eh )

−

2{{Z}}

+ E +n× H ) 
n×(Y

h
h
2{{Y }}

(4.73)

It is important to note that this choice of the flux is not unique. In fact, in so far as this
flux introduces a term proportional to Eh in Ampere’s law, and Hh in Faraday’s law, this
particular form will render the final discrete system dissipative. Consequently, one might
consider removing these terms to yield what is known as a ‘central flux’. While this flux can
be used with symplectic integrators to guarantee discrete energy conservation [125], it does
not damp spurious modes that arise in a nodal DG framework [126].
e
If Equation 4.72 is made to hold for all lp , it describes a system of 6N ODEs in 6N
Ne
unknowns, where N =
Np (e). Here, the unknowns are Eh (re , t) and Hh (re , t), the field
p
p
i=1
values at all interpolation nodes as a function of time. To this end, Equation 4.72 is called
the ‘semi-discrete form’ of the problem, in so far as it remains un-discretized in time. For
95

convenience, it is re-written in terms of the field components Eh and Hh :
e ∂E
e
e
lp , ε h Te = − lp , × Hh Te + lp , FE {Eh , Hh } Fe
∂t
e
e
e
e ∂ Hh
= lp , × Eh Te + lp , FH {Eh , Hh } Fe + lp , J Te
lp , µ
∂t Te

(4.74a)
(4.74b)

Here FE and FH are the relevant components of the upwind flux in Equation 4.73. The
features of the matrices that this form generates are as follows:

1. The matrix elements associated with the elemental inner products taken alone form a
block-diagonal matrix. Each element generates an Np (e) × Np (e) block.
2. The matrix elements associated with the trace inner products taken alone contribute to
both diagonal and off-diagonal blocks. Each interior face produces gives diagonal and
off-diagonal contributions, exterior faces produce only diagonal contributions. Boundary conditions are enforced through these terms, details of which are given in Section
4.2.2.3.

Specific details of the matrix elements themselves can be found in [36].
The final step is the time discretization of Equation 4.74. As both matrices on the
LHS are block-diagonal, it is convenient to integrate this system of ODEs using an explicit
method. In particular, a low-storage fourth order Runge-Kutta [127] is applied. Attendant
to this choice is the limitation that the time step used in this integration, ∆t, is limited by a
CFL-like condition [67] that gives the maximum stable time step. In practice, the time step
is chosen as:
∆t = C

h
c max(Np (e)2 )
96

(4.75)

Here, C ≤ 1 is a scale factor that can be used to decrease the time step beyond the minimum
required for stability, c is the speed of light in the units of the problem, the factor of Np (e)2
accounts for the reduction in the distance between interpolation nodes at higher order, and
h is a mesh-dependent quantity given as:

2
h = min
Ωh 3

Volume of Te
minFe (Area of face)

(4.76)

This will hold for straight-edged tetrahedra, for which the factor of 2/3 comes from the
volumetric and surface Jacobians associated with transforming onto reference elements. For
all results, C = 1 to achieve the largest stable time step.

4.2.2.3 Boundary Conditions

The final piece of the DG framework worth remarking on is the enforcement of boundary
conditions. Simply put, these are imposed by restricting the form of

Eh

and

Hh

at

interior or exterior boundaries. In particular, in so far as E− and H− are set by the problem
h
h
solution, E+ and H+ must be specified. For PEC boundaries, these are given as:
h
h
E+ = −E−
h
h

(4.77a)

H+ = H−
h
h

(4.77b)

Heuristically, one might view this as enforcing the vanishing of the tangential Eh field by
constraining the field on the positive trace as to cancel it. Similarly, for PMC boundaries
97

these are given as:

E+ = E−
h
h

(4.78a)

H+ = −H−
h
h

(4.78b)

This can be heuristically justified in an analogous way. To realize a simple Silver-M¨ller
u
boundary condition, the following conditions are enforced [128, 129]:

E+ = −E−
h
h

(4.79a)

H+ = −H−
h
h

(4.79b)

This form is perhaps not as physically transparent as the PEC and PMC boundary conditions. However, one can see that the upwind flux in Equation 4.73 vanishes when the usual
Silver-M¨ller relationship between tangential Eh and Hh is satisfied. It is also frequently
u
useful to enforce a Dirichlet boundary condition when the fields at a boundary are known
analytically. In this case, the positive traces are given as:

E+ = −E− + 2Eenf
h
h

(4.80)

H+ = −H− + 2Henf
h
h

(4.81)

Here, Eenf and Henf are the enforced fields, and this condition can be physically reconciled
as forcing the positive traces of both fields to cancel the negative trace, and adding in the
imposed fields.

The final type of boundary condition to be discussed arises in the total field/scattered
98

field formalism [49]. This formalism is typically applied in the context of imposing the
incident field in scattering problems. Here, it will be useful in representing the fields due to
a point-dipole excitation, J(r, t), of the form:

J(r, t) = ds js (t)δ(r − rs )

(4.82)

Here, ds is the orientation of the dipole, rs is its location, and js (t) is the temporal signature
of the current driving it. The nearer one gets to the point excitation, the more poorly
approximated the fields will be due to the smooth polynomial basis functions being employed.
The total field/scattered field formalism can be used to remove the singular portion of the
field from numerical consideration, and has been applied to this end for the treatment of line
sources in DG [128]. To achieve this, the fields are decomposed into incident and scattered
components:

Eh = Eh,inc + Eh,scat

(4.83a)

Hh = Hh,inc + Hh,scat

(4.83b)

A boundary is demarcated around the point source, outside of which the total fields Eh and
Hh are solved, and inside which only the scattered components Escat and Hscat are solved.
This scenario is illustrated in Figure 4.24. Eh,inc and Hh,inc are the fields generated by
the source satisfying the Maxwell Equations in free space. Explicit expressions for these
fields are given in Section A.3 of the Appendix. At the boundary between the the total and
scattered field regions, a boundary condition must be imposed that ensures that the two
regions remain consistent. Rather than restricting the form of

99

Eh

and

Hh

through

the positive traces, however, this boundary condition is enforced by augmenting the existing
value by the incident field at the boundary. In other words, the jumps are instead given as:

Eh
Hh

= E+ − E− ± Eh,inc
h
h

(4.84a)

= H+ − H− ± Hh,inc
h
h

(4.84b)

Here, the + sign is utilized when the positive trace is within the scattered field region, and
the − sign is utilized when the positive trace is in the total field region.

4.2.2.4 Non-Conformal Discretizations
In Section 1.2.5, the ability of the DG framework to accommodate non-conformal meshes
was briefly discussed. To this end, when working with such meshes there are no changes
to the discretization, and Equation 4.74 continues to yield the appropriate matrix elements.
However, the evaluation of the matrix elements arising due to the trace terms will change,
and requires non-trivial modifications as far as implementation is concerned. To see this, one
might consider the many ways in which two meshes might meet at a non-conformal interface.
Rather than being aligned, as would be the case in a conformal mesh, the triangular faces
at the non-conformal boundary can overlap in four different ways, illustrated in Figure 4.25.
To this end, new data structures must be introduced to deal with a number of assumptions
that can no longer be made:

1. An element may now have arbitrarily many neighbors, rather than one per face (i.e.,
a maximum of four). Similarly, each face is no longer associated with at most two
abutting tetrahedra.
100

2. The interpolation nodes on two overlapping faces will no longer be exactly degenerate.
This degeneracy is used to simplify the evaluation of matrix elements a great deal in
conformal meshes, as they can be computed on a reference element, stored, and scaled
when they are applied in the Runge-Kutta integration.
In modifying an existing DG implementation to support these assumptions, this new data
type will consist of a list of the separate polygonal fragments that tessellate the non-conformal
interface, each node of which includes the following information:
• Indices of the faces and elements that comprise the fragment.
• Area of the fragment.
• Interpolation operators that map from the associated DG basis nodes on both traces
onto a set of degenerate quadrature nodes, derived from [130].
• Matrix elements associated with the overlap of basis functions on the fragment.
First, however, one must generate a description of the fragments themselves. This is most
easily accomplished using polygon clipping algorithms [131]. Such an algorithm has been
implemented, and an exemplary non-conformal interface processed by it is illustrated in
Figure 4.26.

4.2.3 Validation
In this Section, results validating a Fortran implementation of the scheme discussed above
are presented. These results are intended to demonstrate the correctness of the code for later
use. In sequence, results will be demonstrated that illustrate the following:
101

• Higher order accuracy of the DG method in terms of convergence of analytically known
cavity fields in the L∞ norm.
• Solution of an inhomogeneous Maxwell Equation in which a point dipole excites a PEC
cavity, yielding its spectrum of modes.
• Implementation of piecewise homogeneous materials, as indicated by the modal specturm of a partially filled cavity.
• Convergent representation of the fields due to a dipole in free space using Dirichlet,
total field/scattered field, and absorbing boundary conditions.
The first result validates the proper higher order convergence of the DG implementation,
and is inspired by an example in [36]. Here initial conditions for a TM1,1 cavity mode
are imposed as the boundary data on a PEC cavity, and the resultant E(r, t) and H(r, t)
are integrated out for many cycles. Results are presented in Figure 4.27 that demonstrate
convergence in both polynomial order and the edge length, h. The error in the L∞ norm is
estimated from a sampling of the field at 10,000 randomly chosen points on the interior of
the cavity at the peak of the fifth cycle of the cavity mode. It is evident that the error can
be reduced by either increasing the polynomial order, Np , or decreasing h.
The next sequence of results is the reconstruction of a set of cavity modes from the fields
generated by a dipolar current source excited by a broadband pulse, inspired by results
in [132]. A cubic 1 m PEC cavity is meshed, and 5 dipolar current sources are distributed
throughout, with random positions and orientations. Each source is excited by the derivative
of a Gaussian-modulated pulse with center frequency 540 MHz and an effective bandwidth of
660 MHz. The fields are evolved for 1000 ns, and recorded at 20 random positions throughout
102

the cavity. The discrete Fourier spectrum is computed for each electric field component at
each position. The resultant cavity modes are identified by an averaging of the spectra over
all observers and electric field components, and the results are illustrated in Figure 4.28. It
is evident that the location of the first 21 modes are in good agreement with analytic results.
This validates the ability of this implementation to reproduce the spectral response of a PEC
structure.
The next result in this sequence validates the implementation of piecewise homogeneous
dielectric materials. The eigenmodes of a cubic 1m cavity filled partially (40%) with a
material for which εr = 2 were computed analytically and numerically. Results are illustrated
in Figure 4.29. Again, the cavity is excited by 5 dipole point excitations, with randomly
chosen positions and orientations, each driven by the derivative of a Gaussian-modulated
pulse with center frequency 270 MHz and an effective bandwidth 400 MHz. The fields are
evolved for 800 ns, and the discrete Fourier spectrum is computed from 20 samples at random
positions throughout the cavity. As before, the resultant cavity modes are identified by an
averaging of the spectra over all observers and electric field components. The location of the
first 20 modes are in good agreement with analytic results. In the case of nearly degenerate
modes, the numerically computed peaks are not always split. To see this splitting, a longer
simulation could be run to achieve a better spectral resolution.
Next, a result is presented illustrating the correct implementation of Dirichlet, total
field/scattered field, and Silver-M¨ller boundary conditions. A sphere of radius 1m is meshed
u
with a desired edge length of λ/5 and three separate tests are run. In the first, a 0.5m sphere
is removed from the center of the domain, and Dirichlet boundary conditions are applied on
both the interior and exterior faces. Electric and magnetic fields due to a dipole source at
103

the center of the sphere are imposed on both surfaces, for a current with the time signature
of a Gaussian pulse with center frequency 300 MHz and bandwidth 300 MHz. The pulse
is allowed to propagate through the computational domain, and the relative L2 error is
computed after it has passed for polynomial orders 1 through 5. In the next test, the same
source is utilized, but the Dirichlet boundary condition is only applied at the exterior face,
where it should exactly capture the outgoing character of the fields propagating into an
unbounded medium. At the interior face, the total field/scattered field boundary condition
is applied, and the 0.5m sphere that was removed is again inserted. This region is treated as
a scattered field region, and the exterior spherical shell is treated as a total field region. L2
error is computed in both regions under the same conditions. The final test uses the same
geometry, but rather than applying a Dirichlet boundary condition at the exterior face, a
Silver-M¨ller boundary condition is imposed. Again, the L2 error is computed throughout
u
the domain. The results of this sequence of tests are given in Figure 4.30

Here, there are two points worth noticing. First, the rate of convergence for the total
field/scattered field plus Dirichlet test is greater than that of the pure Dirichlet test. This is
due to the fact that the computational domains are different. In the pure Dirichlet test, there
is no interior to the spherical shell, for which the scattered field region is defined in the other
test. Second, the convergence for the total field/scattered field plus absorbing boundary
boundary condition saturates, whereas neither of the other tests do. This is due to the
presence of spurious scattering from the boundary of the computational domain. More can
be learned about the nature of this error by considering the time signature of the absolute
error integrated over only the spatial extent of the domain. This is illustrated in Figure 4.31.
104

While the total field/scattered field plus absorbing boundary condition error converges
uniformly in polynomial order up until t ≈ 10 ns, beyond this time, all orders have the
same time signature. This error is not due to a deficiency in the choice of representation,
as it eventually overtakes the error in the 1st order total field/scattered field plus Dirichlet
data, but is instead due to spurious reflections at the approximate absorbing boundary
condition. Normalizing the error in terms of the amplitude of the exciting field indicates
that slightly less than 1% of the power incident on this boundary is spuriously reflected,
which is consistent with the level at which the relative error was found to saturate in Figure
4.30. This is a well-known issue associated with utilizing a Silver-M¨ller condition of this
u
order. To reduce the size of these reflections, the termination of the domain could be placed
further from the source, or an alternative transparent boundary condition could be applied,
such as a perfectly matched layer, or a pole condition-based truncation. The implementation
and exploration of these types of boundary conditions are left to future work, and the level
of reflection associated with this simple boundary condition must simply be left as a caveat
to the user.
Having established the validity of numerous components of the DG implementation, all of
the necessary components of the code are suitably vetted for the application to be presented
in Chapter 5.

4.3 Figures

105

Figure 4.24: Total field/scattered field separation for imposing fields due to a dipole soucee.
The location of the source is indicated by the arrow, the red annulus is the total field region
and the blue circle is the scattered field region. At the boundary between the two regions,
the condition in Equation 4.84 is enforced.

3 vertices

4 vertices

5 vertices

6 vertices

Figure 4.25: Four distinct types of fragments that can occur in a tetrahedral non-conformal
mesh.

106

Figure 4.26: Quadrature at a non-conformal interface. Triangles on one side of the interface
are outlined in black, triangles on the other are outlined in red. Blue points indicate the
vertices of fragments, and green points indicate the quadrature points used for integration.
Fragments have been further decomposed into triangles upon which a symmetric higher
order quadrature is defined [130]. Some pathologies are also illustrated, as a result of finite
precision arithmetic.

107

Error Convergence for a TM1,1 Cavity Mode
0

Absolute RMS Error, Frobenius Norm

10

10-2
10-4
-6

10

h=0.500
h=0.333
h=0.250
h=0.200
h=0.160
h=0.140
h=0.125
h=0.110
h=0.100

-8

10

10-10
10-12
-14

10

1

2

3

4

5
6
7
Polynomial Order

8

9

10

Figure 4.27: Error in E(r, t) of a TM1,1 cavity mode in the L∞ norm after 5 cycles. Convergence is demonstrated with respect to edge length and polynomial order.

108

Eigenmodes of an Empty Cavity
2

Average Modal Amplitude (a.u.)

10

Numerical Modes
Analytic Modes

101
0

10

10-1
10-2
10-3
10-4
-5

10

100

200

300

400
500
Frequency (MHz)

600

Figure 4.28: Eigenmodes of an empty 1m cubic cavity.

109

700

800

Eigenmodes of a Partially Filled Cavity

Average Modal Amplitude (a.u.)

102

Numerical Modes
Analytic Modes

101
0

10

10-1
10-2
10-3
10-4
-5

10

100

150

200

250
300
Frequency (MHz)

350

400

Figure 4.29: Eigenmodes of a 1m cubic cavity for which the bottom 40% is filled with a
material with εr = 2, with free space elsewhere.

110

Convergence of Fields due to a Point Dipole
100

Dirichlet BC
TF/SF + Dirichlet BC
TF/SF + ABC

-1

Relative L2 Error

10

-2

10

10-3
10-4
-5

10

-6

10

1

2

3
Polynomial Order

4

5

Figure 4.30: Convergence of the relative L2 error of the fields due to a point dipole excited
by a Gaussian modulated pulse with center frequency 300 MHz and bandwidth 300 MHz at
the center of a 1m spherical domain. Boundary conditions are imposed using pure Dirichlet, Dirichlet plus total field/scattered field, and absorbing plus total field/scattered field
boundary conditions.

111

Spurious Reflections Associated With Absorbing Boundary Condition
0

10

10-1

Absolute L2 Error

10-2
-3

10

-4

10

10-5
-6

10

-7

10

TF/SF + Dirichlet (Order 1)
TF/SF + ABC (Order 1)
TF/SF + ABC (Order 2)
TF/SF + ABC (Order 3)
TF/SF + ABC (Order 4)

-8

10

10-9
-10

10

0

5

10

15
20
Time (ns)

25

30

35

Figure 4.31: Absolute L2 error of the fields due to a point dipole excited in the same manner
as in Figure 4.30.

112

Chapter 5: DG-QED for Cavity Mediated Energy Transfer

5.1 Overview of Energy Transfer

In this Chapter, the time domain Maxwell solver discussed in Chapter 4 is applied to the
analysis of systems in which a number of quantum mechanical subsystems interact through
their mutual coupling to the EM field. This type of modeling can be applied to the spectrum
of applications identified in Section 1.1.2, including energy transfer and quantum information.
Here, results will be limited to proofs of concept. In Section 5.2, details of the modeling
formalism will be identified. The quantum mechanical subsystems will be described in terms
of a set of two-level quantum systems (2LSs) that are coupled to the EM field through a
multipolar Hamiltonian truncated at the dipole level. The quantum master equation that
describes the evolution of the reduced density matrix associated with only the 2LSs will be
given, and coefficients will be identified in terms of quantities that can be computed using
a classical Maxwell solver. In Section 5.2, an algorithm for extracting the relevant model
parameters from a sequence of classical time domain Maxwell solutions is presented. Section
5.3 will include results that illustrate two proofs of concept, namely spontaneous emission
in the vicinity of a PEC plane, and the coupling of two quantum dots in a semiconductor
microcavity. Finally, Section 5.4 will outline future work.
113

5.2 Cavity QED Formalism

5.2.1 Model Hamiltonian and Master Equation

Consider an open domain, Ω ⊂ R3 , with boundary ∂Ω, some subset of which may be PEC,
∂ΩC , and a distribution of classical dielectric matter characterized by the dielctric function
ε(r, ω) = ε0 ε (r, ω) + iε (r, ω)

with analytic properties consistent with the Kramers-

Kronig relations [133]. For simplicity, it is assumed that the matter is non-magnetic and
characterized by permeability µ0 , noting that the extension to magnetic media will require
nothing that is not analogously encountered in accommodating a dispersive dielectric. Some
distribution of Nd two-level systems (2LSs) is contained within Ω, such that the mth 2LS
is characterized by the position of its centroid, rm , transition frequency ωm , and dipole
transition moment dm . The Hilbert space of each 2LS consists of the usual complex 2vectors of unit norm, which are operated upon by the usual local SU(2) basis, i.e., the
x y z
Pauli matrices plus the identity {Im , σm , σm , σm }. Raising/lowering operators are built
y
±
1 x
from these matrices, σm = 2 (σm ± iσm ), and used to define dipole transition operators
−
+
dm = dm σm + d∗ σm . The 2LSs will be coupled to the electromagnetic field within a
m
multipolar formalism, truncated at the electric dipole level [19, 134].

H=

dr
Ω

∞

†
dω ω f (r, ω)f (r, ω) +

0

Nd ∞
+ −
ωm σm σm −
dω dm · E(rm , ω) + h.c.
m=1
m=1 0
(5.85)
Nd

Here, f (r, ω) and its adjoint are bosonic fields that obey the usual equal-frequency, equalspace commutation relations that generate the EM field, with the quantum fields, E(rm , ω)
114

and B(rm , ω), obeying:

E(rm , ω) = i

B(rm , ω) =

i
ω

ω2
π c2

ε (r , ω)G(rm , r , ω) · f (r , ω)

dr

(5.86)

Ω

× E(rm , ω)

(5.87)

Here, G(r, r , ω) is the classical dyadic Green’s function for the electric field that satisfies
the following inhomogeneous vector Helmholtz equation:

×

× G(r, r , ω) −

ω2
(ε (r, ω) + iε (r, ω))G(r, r , ω) = δ(r − r )
c2

(5.88)

subject to the constraint that the generated fields obey classical boundary conditions on ∂Ω,
and throughout materials described by ε(r, ω).

Using this Hamiltonian, the Heisenberg equations of motion can be derived in the usual
way for some operator, O:
dO
i
=
H, O
dt

(5.89)

− +
Here, O may be any of the relevant dynamical variables, e.g., σm , σm , and f (r, ω). As there
is a continuum of EM field variables, it is favorable to integrate them out, and to instead
study the reduced density matrix, ρs (t) = Tr(ρ)EM . The subscript indicates that the trace
is evaluated over the EM field variables, and that only the 2LS sector of the Hilbert space
remains in the reduced density matrix. A master equation that governs the evolution of this
quantity has been derived in the limit that the relative detuning between 2LSs is small over
115

the scale of frequency variations in the system Green’s function [134]. It is given as:

∂ ρs
i
=−
∂t
2

Nd

z
ωm σm , ρs + i

m=1

...

Nd
m=n

+ −
gm,n σm σn , ρs − . . .

1
−
+ −
+
+ −
Γm,n σm σn ρs − 2σn ρs σm + ρs σm σn
2 m,n

(5.90)

Here, ωm = ωm − gm,m is the Lamb-shifted transition frequency, and the quantities gm,n ,
and Γm,n can be derived from the classical dyadic Green’s function governed by Equation
5.88:
2
ωn
dm · Re{G(rm , rn , ωn )} · dn
ε 0 c2
2
2ωn
Γm,n =
dm · Im{G(rm , rn , ωn )} · dn
ε 0 c2
gm,n =

(5.91)
(5.92)

For canonical systems in which G(r, r , ω) has a closed-form expression, these parameters
can be computed, and the master equation can be integrated numerically using something
simple, such as Runge-Kutta. In general G(r, r , ω) has no closed-form expression, and it is
necessary to consider methods for numerically generating relevant samples of G(r, r , ω). In
what follows, the manner in which the time domain Maxwell solver from Chapter 4 can be
used to this end will be discussed.

5.2.2 Extraction of Model Parameters
In general, a complete description of G(r, r , ω) would require the solution of Equation
5.88 for all source-observer pairs in Ω, at all frequencies. Fortunately, Equation 5.90 only
2
requires samples for a finite number of source-observer pairs, namely Nd , at frequencies in
116

the set of 2LS transition frequencies, {ωm |m ∈ 1 . . . Nd }. These samples can be extracted
from field data generated from Nd separate solutions of the classical time domain Maxwell
Equations. The mth simulation will consist of the following:

1. A point dipole excitation at rm , oriented along the direction defined by dm is excited
by the derivative of a Gaussian modulated pulse. To this end, the source term is
rendered:

(t−td )2
−
d 
2σ 2  δ(r − rm )
Jm (r, t) = dm sin(ωc (t − td ))e

dt


Jm (r, t) = dm j(t)

(5.93)

(5.94)

Here, td is some shift in time chosen to ensure that the system is approximately
quiescent at t = 0, and ωc and σ define the center frequency and bandwidth of the
frequency interval over which G(r, r , ω) is desired.

2. The total classical electric field will be recorded at each of the remaining Nd − 1 2LSs,
and the scattered electric field will be recorded at the mth 2LS. These data will be
Fourier transformed to yield E(rn , ω) and Escat (rm , ω).

3. Entries of the dyadic Green’s function can then be extracted from:

E(rn , ω) = −iωµ0 G(rn , rm , ω) · Jm (rm , ω)

(5.95)

Escat (rm , ω) = −iωµ0 Gscat (rm , rm , ω) · Jm (rm , ω)

(5.96)

117

Here the scattered field dyadic Green’s function, Gscat (r, r , ω) is defined as:

G(r, r , ω) = G0 (r, r , ω) + Gscat (r, r , ω)

(5.97)

Where G0 (r, r , ω) is the free space dyadic Green’s function. This quantity is so-defined
as to regularize the singular behavior of G(r, r , ω) as r → r .

4. From these samples, gn,m and Γn,m can be computed. The self-terms, gm,m and
Γm,m , must be augmented by the free space contribution due to G0 (r, r , ω) in Equation 5.97. These contributions are simply the Lamb shift and emission rate of an
isolated emitter in free space.

5.3 DG-QED Results

5.3.1 Spontaneous Emission Above a PEC Half-Space

Analytical results due to Wylie and Sipe [135, 136] for a single quantum emitter above
a PEC half-space are used to validate the presented method of calculation. Here, Equation
5.90 reduces to a single parameter, Γ, where the subscript has been removed. Γ sets the rate
of spontaneous emission, which will changing depending upon whether or not the emitter’s
dipole moment is parallel or perpendicular to the plane of the half-space. The respective
rates, Γ and Γ⊥ , are most easily expressed normalized against the free space spontaneous
118

ω 3 |d|2
emission rate, Γ0 =
:
3 ε 0 c3
1
1
3 cos(η)
+
−
Γ /Γ0 = 1 −
2
η η3
η2
sin(η) cos(η)
Γ⊥ /Γ0 = 1 + 3
−
η3
η2

sin(η)

(5.98)
(5.99)

Here, η = 2zω/c is a unitless parameter describing the distance of the emitter from its own
image in the half-space, rationalized by the emission wavelength.
The equivalent image problem is simulated, using the DG solver from Chapter 4. Here,
a dipolar current source is placed at the center of a spherical domain and excited with
the derivative of a Gaussian pulse. Free space termination is simulated using an absorbing
boundary condition, and a 2nd order complete basis is utilized. The original problem and its
equivalent image problem are illustrated in Figure 5.32. The Green’s function is extracted
using the procedure outlined in Section 5.2.2, from which the normalized emission rates are
computed. The numerical results are compared against the analytic results in Equation
5.98 in Figure 5.33. The agreement between the computed rate and the analytic rate is as
expected. The L2 error in the fields is on the order of 1% throughout the simulation domain.
Slightly larger deviations exist in the computed rates, and are most evident in the tail of the
parallel-oriented decay rate. This is attributed to additional error incurred in interpolating
the discrete Fourier spectrum.
The calculation was attempted for η < 1, corresponding to the emitter being less than
0.08λ from the PEC surface. While this type of nearfield calculation is possible, it becomes
very difficult to resolve the imaginary part of the Green’s function. This can be rationalized
in terms of the singular behavior of the half-space dyadic Green’s function. The real part
119

diverges as ∼ (z/λ)−1 , whereas the imaginary part goes to a constant. Consequently, while
the fields may be resolved well in the L2 norm, this may be dominated by the representation
of the large real part, rather than the small imaginary part, requiring a significant increase
in the order of the DG basis. As this will also require a decrease in the element size to
represent the scattered field region, the CFL-like condition will become even more restrictive, and the accurate calculation of the extreme nearfield of the emitter rapidly becomes a
computationally expensive endeavor.

5.3.2 Energy Transfer in a Semiconductor Nanocavity
In the next result, the coupling parameters associated with two quantum dots inside of a
semiconductor nanocavity are analyzed. To this end, each quantum dot is modelled as a 2LS
with a simple excitonic transition yielding the gap between the ground and excited states.
Here, a cavity is formed from by a 260 nm thick layer of GaAs (n = 3.5) sandwiched between
two AlAs/GaAs Distributed Bragg Reflectors (DBR) (n = 3.0, n = 3.5), each consisting of 6
alternating quarter-wave slabs. Previously, others have analytically characterized long-range
energy transfer processes between quantum dots in free space [137, 138], as well as in such
a cavity with more DBR layers [139], and at shorter ranges in a planar PEC microcavity
[140]. A mesh of the computational domain, with dimensions, is illustrated in Figure 5.34.
The procedure for extracting inter-2LS couplings as outlined in Section 5.2.2 is implemented to study the coupling of a quantum dot at the center of this structure, with another
placed a distance d away from it along the axis of the cavity. The moments of both dots
are aligned, and oriented such that they point along the short edge of the computational
domain, and negligible fields are radiated in this direction. The coupling parameters are
120

investigated over a spectrum of energies centered at 1.36 eV, the cutoff of a guided mode of
the system, with a FWHM of 0.68 eV. Fields are evolved for a period of time sufficient for
the entire pulse to pass through the domain, and for the cavity fields to return to a quiescent
state. Snapshots of the electric field intensity are illustrated in Figure 5.35.
In this image, the first snapshot is taken as the current source amplitude is peaked. The
second and third are taken after the source has died down to a negligible level and give some
indication of the guided mode structure being excited. By the end of the simulation, the fields
have died down to the extent that they are unobservable on the scale of the visualization,
i.e., they are at least 3 orders of magnitude below the peak.
The electric field component parallel to the exciting dipole is sampled along the central
axis of the GaAs nanocavity, and Fourier transformed to give an indication as to how the
cavity dresses propagation. The interval over which the field is sampled is indicated in the
top part of Figure 5.35. Figure 5.37 illustrates the Fourier transformed field component over
the spectral range for which the energy of the exciting current is above 10% of its peak, from
a distance of 40 nm from the exciting source, to the edge of the computational domain 1040
nm away. For reference, Figure 5.36 illustrates the same quantity over the same scale, but
in free space.
The Gaussian nature of the excitation is evident in the free space data in Figure 5.36.
That the cavity strongly modifies this field distribution is then evident in Figure 5.37. While
it is still evident that the fields are generated by a source with a Gaussian spectrum, it is
clear that the cavity dramatically modifies this even at short range. One of the more relevant
features is that the spectral content at and below 1.36 eV appears to be reduced relative to
the content above this energy. This is a cutoff effect, in so far as 1.36 eV corresponds to an
121

in-GaAs wavelength of 260 nm, the height of the cavity, for which TE and TM guided modes
exist [139]. Spectral content above this energy/below this wavelength will be able to match
the dispersion relation of these modes and propagate, while content below this energy/above
this wavelength will evanesce. In so far as this mode has a finite Q, this cutoff is not strict,
and energy is coupled into leaky modes rather than evanescing as they would in a lossless
cavity. These leaky modes exist independent of wavelength over the interval of interest. As
the DBRs being modeled only consist of 3 layers above and below the cavity, the Q is not
as high as is often experimentally realized, and some fields are still evident below the guided
mode cutoff. Even so, the strong field between 1.36 and 1.9 eV from 250 nm to 540 nm can
still be attributed to the dipole coupling to some low- to moderate-Q TE guided mode.

Next, the inter-2LS light-matter couplings are extracted using the procedure outlined in
Section 5.2.2. In Figures 5.39 and 5.38, the modulus of the inter-2LS Green’s functions are
illustrated for the DBR cavity and in free space for reference. In free space, a strong shortrange interaction is evident, corresponding to the tail of the usual F¨rster energy transfer,
o
scaling as d−3 . Further away, this scaling transitions to a slower d−1 , as the radiated fields
begin to dominate. The situation is very much different in the DBR cavity. Here, a cutoff
effect is evident in the short-range portion of the interaction, in which the strength of the
inter-2LS coupling is diminished over a band of energies near the cavity cutoff. Further
away, the interaction is enhanced by at least an order of magnitude relative to free space,
both above and below cutoff. Below cutoff, the interaction strength decays more slowly than
in free space. Above cutoff, the interaction is oscillatory, with a sharp suppression of the
interaction corresponding to nodes of the cavity field evident in Figure 5.37, on top of a weak
decay. For a higher Q structure, it is expected that the relative enhancement/suppression of
122

the interaction above cutoff would be more significant, with a maximum on resonance with
the guided mode, and the coupling to the leaky modes below cutoff would yield a slowly
decaying but weaker interaction. In particular, the results in [139] indicate that in a higher
Q DBR system with the system dipoles aligned in the same way, the resonant interaction is
mediated by a TE mode that gives rise to a characteristic d−1/2 decay in the interaction
strength. A similar long-range interaction is evident in the results of the calculation in this
Section - while the relative size of the enhancement/suppression above and below cutoff is
different, the distance-dependence of the cavity-mediated interaction scales similarly.
To consider the interaction in more detail, the real and imaginary parts of the Green’s
function are plotted for 2LSs with transitions in three regimes - below cutoff, at cutoff, and
above cutoff in Figures 5.40, 5.41, and 5.42. Again, it is evident that the oscillations in the
net interaction strength become stronger at or above cutoff, corresponding to suppression at
the nodes of the guided mode observed in Figure 5.37. By resolving the real and imaginary
parts of the Green’s function individually, the relative strength of the coherent and incoherent
terms in Equation 5.90 interactions can be determined. This type of knowledge is especially
useful for quantum information applications in which the preservation of coherence is vital.

5.4 Future Work
Ultimately, the analysis of higher Q structures is desirable, both to see more dramatic
effects in terms cavity-mediated interactions, and to explore systems of relevance to quantum
information applications. To increase the Q of the structure studied in Section 5.3.2, a
number of steps might be taken. The simplest would be to add more layers to the DBRs
forming the cavity. This would increase the computational effort required for the solution of
123

the problem, but is otherwise tractable given enough CPU hours. Another, would be to use
a higher contrast in index between DBR layers. The GaAs/AlAs example presented has a
relatively low dielectric contrast simply because it is experimentally possible to realize this
type of structure. While computational methods have no such restrictions, this does raise
questions of relevance. However, even with these strategies, some other issues arise in the
accurate analysis of high Q systems using time domain methods. One is evident in Figure
4.28. Here, the modal structure of a PEC cavity was resolved using a technique similar
to the one employed in this Chapter. Theoretically, the modes being interrogated have an
infinite Q, in so far as there are no losses in the system being modeled. Instead, a spurious
numerical linewidth is evident, that makes it appear as if these modes have a finite lifetime
due to the finite duration of the time signature from which they are extracted. Further,
the dissipative nature of the upwind flux introduces numerical losses that contribute to this
broadening, especially as longer simulations are required to resolve higher Q resonances.
Ultimately, this is not a significant deficiency for energy transfer applications, and further
explorations of systems similar to the one in Section 5.3.2 will be carried out. So long as the
light-matter coupling is weak to moderate, this approach is eminently practical. However,
as the transition into strongly coupled CQED is approached, as would be the case with
quantum information applications, the high fidelity representation of cavity Q will become
increasingly important.

To this end, a reasonable solution is to model the Q factor directly, by solving a Maxwell
eigenproblem on an open domain. The discretization of the Maxwell eigenproblem is possible
within the DG framework, so the underlying basis functions and code base can be reused.
Immediately, one may be concerned that the eigenanalysis may be polluted by spurious
124

modes 1 , as is the case in nodal CG frameworks [141], however, the use of interior penalty
(IP) factors [142] has been shown to alleviate this in DG. To this end, Buffa and Perugia
[143] have outlined a set of conditions necessary for declaring that a DG discretization is free
of spurious modes, and have subsequently shown that numerous IP and local DG schemes
meet these criteria. Even in the case of non-conformal meshing, spurious solutions have been
empirically demonstrated to be of little concern [65].
Even so, this analysis has been carried out for problems on closed domains. A structure
with a finite Q is necessarily defined on an open domain. While the Silver-M¨ller boundary
u
condition used earlier in this Dissertation can be used to effect this, it has some disadvantages
as it must be placed sufficiently far from a given structure that the fields incident on it must
be transverse, lest spurious reflections arise. A PML [40] can be used to bring the effective
boundary closer, and will generally result in smaller spurious reflections. However, i.) PML
optimization will become an issue, ii.) the PML may itself give rise to spurious modes, and
iii.) error may be less controllable than would be desired for high Q systems. Even with a
PML, if the open domain has some structure, as was the case with the DBR result presented
in Section 5.3.2, it is lost in the finite truncation of the computational domain, i.e., the PML
does not preserve layering of media, and neither does the Silver-M¨ller boundary condition.
u
In practice, this can be mitigated by constructing a mesh that is sufficiently large that the
absence of infinite stratification is not noticeable for a point source excitation at its center,
and this is exactly what was done in Section 5.3.2. Ideally, however, a method that can
achieve the same result with a smaller mesh while preserving this structure is desirable.
To this end, future work will see the development of an IPDG framework for a pole
1 Fundamentally, this would also be an issue for the transient analysis, as was explored in
the last section of [126].
125

condition-based transparent boundary conditions that requires a simpler optimization, has
been shown to produce nice error bounds for CG methods, and accommodates unbounded
inhomogeneities. The crux of this method is the identification of the Sommerfeld radiation condition with conditions on the holomorphic extension of the Laplace transform of
the asymptotic field [44, 45]. Considerable work has been done in realizing a numerically
tractable means of generating solutions that satisfy these conditions. Perhaps the most successful approach to date is based upon the M¨bius transform of the Laplace transform of the
o
asymptotic field [46]. Further, these methods involve the use of infinite elements, which can
be constructed in such a way as to allow for the explicit modeling of infinite stratification
for structures like a DBR, and have already been applied to this end for planar waveguiding
structures for scalar wave problems [144]. These methods can also be applied in the time
domain, and implementing them provides a logical first step towards achieving a more accurate analysis of lossy cavities within the time domain framework already developed in this
Chapter, prior to shifting the theoretical focus to one based upon eigenanalysis.

In switching over to an eigenanalysis-based approach, a question naturally arises germane
to the underlying quantum mechanical model - how can the master equation coefficients be
extracted from the Green’s function if the method outlined in Section 5.2.2 determines them
from time domain data? The Green’s function can instead be constructed using a resolvent
formalism, i.e., an expansion in eigenfunctions/eigenvalues of the field. While this will suffice
for extracting parameters for the model identified in Section 5.2, alternative models may
also prove promising in eliciting more detailed physics. In particular, a modified Maxwell
eigenproblem was recently proposed in [27], in which the full light-matter eigenstates for
a system with multiple quantum subsystems can be resolved directly. This formalism also
126

includes the proper treatment of the non-local susceptibility associated with coupling to an
exciton, as would be relevant to a more first principles treatment of field-mediated coupling
between quantum dots.
This raises one final issue to be confronted in future work. In the theory in Section 5.2,
the light-matter coupling is still fundamentally phenomenological, in so far as the dipole
moments have not been related to any higher theory that derives the 2LS parameters from
something like an electronic structure model. To this end, closing the loop and performing
ab initio electronic structure calculations to inform the QED model would provide a level of
detail that has not yet been realized. Actual computation of excitonic wavefunctions may
require the retention of terms beyond the dipolar truncation of the multipolar Hamiltonian
given in Equation 5.85. Recent experiments have shown that in plasmonic systems with
strong field gradients, the light-matter coupling can become sensitive to corrections that
arise beyond the dipole approximation, justifying the application of this level of modeling
[145].

5.5 Figures

127

Figure 5.32: (Left) A single quantum emitter above a PEC half-space. (Right) Its equivalent
image problem. Arrows indicate the direction of the associated dipole moments.

128

Emission Rate for a 2-Level System Above a PEC Surface

Normalized Emission Rate

2

Parallel (Numerical)
Parallel (Analytical)
Perpendicular (Numerical)
Perpendicular (Analytical)

1.5

1

0.5

0
1

2

3

4

5

6 7 8 9 10 11 12 13 14 15 16
Normalized Distance

Figure 5.33: Analytical and numerical spontaneous emission rate of a two-level system above
a PEC half-space with its dipole moment parallel and perpendicular to the interface.

129

Figure 5.34: Mesh of the GaAs nanocavity in which energy transfer between two 2LSs is
analyzed. All exterior boundaries are terminated with a Silver-M¨ller boundary condition.
u

130

Figure 5.35: Images of the electric field magnitude at three points in time. All images are
taken from a plane cut halfway through the short axis of the mesh in Figure 5.34. (Top) At
the peak of the exciting current. The black line indicates the path along which the fields are
sampled to generate the Green’s function. (Middle) Half way through the simulation, the
exciting current pulse is 4 orders weaker than at peak. (Bottom) Three quarters of the way
through the simulation, the exciting current pulse is 20 orders weaker than at peak.
131

Field Strength (a.u.)

−5

x 10
10
9

1.1

8

1.2

7

1.3

6

1.4

5

1.5

4

1.6

3

1.7

2

1.8

Energy (eV)

1

1
0
100 200 300 400 500 600 700 800 900 1000
Distance (nm)

Figure 5.36: Distance dependence of the Fourier spectrum of the electric field component
parallel to an exciting dipole in free space.

132

Field Strength (a.u.)

−5

x 10
10
9

1.1

8

1.2

7

1.3

6

1.4

5

1.5

4

1.6

3

1.7

2

1.8

Energy (eV)

1

1
0
100 200 300 400 500 600 700 800 900 1000
Distance (nm)

Figure 5.37: Distance dependence of the Fourier spectrum of the electric field component
parallel to an exciting dipole in the center of the DBR cavity.

133

Modulus of Inter−2LS Green’s Function

−4

1
−4.5

1.1

Energy (eV)

1.2

−5

1.3
−5.5

1.4
1.5

−6

1.6
1.7

−6.5

1.8
−7
100 200 300 400 500 600 700 800 900 1000
Distance (nm)
Figure 5.38: Log10 of the modulus of the inter-2LS Green’s function in free space.

134

Modulus of Inter−2LS Green’s Function

−4

1
−4.5

1.1

Energy (eV)

1.2

−5

1.3
−5.5

1.4
1.5

−6

1.6
1.7

−6.5

1.8
−7
100 200 300 400 500 600 700 800 900 1000
Distance (nm)
Figure 5.39: Log10 of the modulus of the inter-2LS Green’s function in the cavity.

135

−5

1

Energy Transfer Coefficients

x 10

Real(G)
Imag(G)
Abs(G)

0.8
GF Couplings (a.u.)

0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
100

200

300

400 500 600 700
Distance (nm)

800

900 1000

Figure 5.40: The inter-2LS Green’s function in the cavity for two systems with a transition
below cutoff.

136

−5

1

Energy Transfer Coefficients

x 10

Real(G)
Imag(G)
Abs(G)

0.8
GF Couplings (a.u.)

0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
100

200

300

400 500 600 700
Distance (nm)

800

900 1000

Figure 5.41: The inter-2LS Green’s function in the cavity for two systems with a transition
resonant with the cutoff.

137

−5

1

Energy Transfer Coefficients

x 10

Real(G)
Imag(G)
Abs(G)

0.8
GF Couplings (a.u.)

0.6
0.4
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
100

200

300

400 500 600 700
Distance (nm)

800

900 1000

Figure 5.42: The inter-2LS Green’s function in the cavity for two systems with a transition
above the cutoff.

138

APPENDICES

139

A.1 Bloch-Periodic Green’s Functions
The intent of Sections A.1.1 and A.1.2 is the clarification of details concerning the Blochperiodic Green’s function used in Chapters 1.4 and 2.4. Derivations of the Bloch-periodicity
condition and the Ewald representation of the Green’s functions are included.

A.1.1 Bloch-Periodicity Condition
The Bloch-periodic Green’s function for the Helmholtz equation is the solution to the
following inhomogeneous PDE:

2 + κ2 G (r, r ) =
µ
t(nµ )∈Lµ

eiκ0 ·t(nµ ) δ(r − r − t(nµ ))

(100)

Physically, Gµ (r, r ) describes the scattering of a planewave excitation of unit intensity and
wavevector κ0 ∈ R3 on linear δ-localized obstructions at the points of Lµ . One is typically
concerned with elastic scattering, i.e., κ = |κ0 |, though much of the subsequent derivation
will remain agnostic to this.
Before deriving Gµ (r, r ), it is useful to apply symmetry principles to derive the Blochˆ
periodicity condition cited in Equation 3.37. Consider a linear operator, T
that transt(mµ )
lates the relative coordinates r − r by t(mµ ) ∈ Lµ . Applying this operator both sides of
Equation 100 yields:

ˆ
T
t(mµ )

2 + κ2 G (r, r ) =
µ
t(nµ )∈Lµ

eiκ0 ·t(nµ ) δ(r − r − t(nµ ) − t(nµ ))

140

(101)

Substituting t(lµ ) = t(mµ ) + t(nµ ), noting that the sum of t(nµ ) over Lµ is equivalent to
summing t(lµ ) , and restoring the dummy index l to n yields:

ˆ
T
t(mµ )

As T

t(mµ )

2 + κ2 G (r, r ) = e−iκ0 ·t(mµ )
µ
t(nµ )∈Lµ

eiκ0 ·t(nµ ) δ(r − r − t(nµ )) (102)

commutes with the Helmholtz operator, it can be determined that Gµ (r, r )

obeys:
ˆ
T
G (r, r ) = Gµ (r − r − t(mµ )) = e−iκ0 ·t(mµ ) Gµ (r − r )
t(mµ ) µ

(103)

Which is the defining characteristic of a Bloch-periodic function.

A.1.2 Ewald Representation

For brevity of the resultant expressions, assume that κ0 = 0, noting that it is trivial
to restore in the resultant expressions given the relation in Equation 103. The simplest
representation of Gµ (|r−r |) is the direct one, in which it is furnished in terms of a summation
of the associated free space Green’s functions translated over Lµ .
Gµ (|r − r |) =

t(nµ )

G(|r − r − t(nµ )|)

(104)

Here, the relevant free space Green’s function is given as:

G(|r − r |) =

e−iκ|r−r |
4π|r − r |

141

(105)

It is evident that for real κ, Equation 104 will yield a conditionally convergent sum. Consequently, an alternative representation of the Green’s function is essential for computational
purposes. The Ewald representation [112] is one such representation that converges spectrally both in and away from the µ-lattice. The starting point for this representation is the
following integral representation of the Green’s function in Equation 105:
e−iκ|r−r |
1
G(|r − r |) =
=
4π|r − r |
2π 3/2

2
−σ 2 |r−r |2 + κ 2
4σ
dσ e

∞

(106)

0

A derivation of this identity via a simple substitution is given in [146]. In so far as Equation
106 has a real-valued integrand along the real line, yet the result is ostensibly complexvalued, this is a particularly interesting identity. To this end, in order for the integral to be
well-defined, the path of integration must be carefully deformed to avoid non-analyticity as
σ → 0 on the real line. The contour appropriate for this, can be traced back to Ewald [147],
and is discussed in a bit more detail in [83].

The next step of the derivation is the splitting of the integration in Equation 106 over
[0, ∞) into separate integrals on [0, η] and [η, ∞), where η ∈ R+ . The periodic Green’s
function can then be represented as:
η
Gµ (|r−r |) =

t(nµ )∈Lµ 0

dσ G(σ, |r−r + t(nµ )|)+

∞
t(nµ )∈Lµ η

dσ G(σ, |r−r + t(nµ )|)
(107)

Here, the form of G(σ, |r − r + t(nµ )|) is evident from Equation 106. In the first term, summation and integration are exchanged, and then the Poisson summation identity is applied
142

[148].
η
t(nµ )∈Lµ 0

η
dσ G(σ, |r − r + t(nµ )|) =

dσ
k(nµ )∈L∗ 0
µ

dµ r e−ik(nµ )·r G(σ, |r |)
(108)

In this form, the integrals in both the real and reciprocal sums can be evaluated [149], leading
to the Ewald representation of the periodic Green’s function in Equation 3.39.

A.2 Periodic Translation Operator Components
In this Section, expressions are provided for the terms in the summation in Equation
3.51. For brevity of the resultant expressions, it is assumed that κ0 = 0, noting that it is
trivial to restore in the resultant expressions. First, the terms in the real sum are considered.
The integral representation of the real sum terms given in Equation 107 can be used to arrive
at an expression that avoids the differentiation of complementary error functions. To this
end, the derivatives are pulled inside of the integral and applied to the integrand prior to
integration:
(p) E (R) p , p , p =
r
x y z

1
=
2π 3/2
1
=
2π 3/2
(−1)p
=
2π 3/2

2 2 κ2
px py pz −σ |R| + 4σ 2
dσ ∂x ∂y ∂z e

∞
η
∞

2 2 κ2
p H (σR )H (σR )H (σR )e−σ |R| + 4σ 2
dσ(−σ) px
x py
y pz
z

η
∞

p

dσ

η

2 2 κ2
px ,py ,pz p+m −σ |R| + 2
4σ
Cm
σ
e

m=0
143

(109a)

(109b)

(109c)

px ,py ,pz
Here, Cm
is the coefficient of the term which is mth order in σ in the product of
Hermite polynomials in Equation 109b. To evaluate this integral in closed form, expand the
exp(κ2 /4σ 2 ) term in a Laurent series in σ.

(p) E (|R|) p , p , p =
r
x y z

(−1)p
=
2π 3/2
(−1)p
=
2π 3/2

∞

p

px ,py ,pz
Cm

dσ
m=0

η
p

∞

px ,py ,pz
Cm

m=0 ν=0

(−1)p
=
4π 3/2 |R|

p

∞

m=0 ν=0

∞

ν=0
∞
dσ

η

(κ/2σ)2ν p+m −σ 2 |R|2
σ
e
ν!

(110a)

(κ/2σ)2ν p+m −σ 2 |R|2
σ
e
ν!

(110b)

px ,py ,pz (κ2 |R|2 /4)ν Γ
Cm
ν!

p+m+1
− ν, η 2 |R|2
2
|R|p+m

(110c)

Here, Γ(n, x) is the incomplete Gamma function of nth order. It is worth noting that all
sums over incomplete Gamma functions are rapidly convergent for physically relevant values
of the arguments, and recurrence relations are utilized to rapidly compute each term.

The derivation of the translation operator components corresponding to terms in the
reciprocal sum is much simpler, and only the final expressions are presented:

(p) E (R, k(n )) p , p , p =
µ
x y z
k
144

ik(n1 )·rl ∞ (−η 2 )ν
α2 (n1 )
px e
= (ikx )
E1+ν
× ...
4πA1
ν!
4η 2
ν=0
ν
2ν−2m−py 2m−pz
ν
(2ν − 2m)!(2m)!
Ry
Rz
...
m (2ν − 2m − py)!(2m − pz)!
m=0
2ν−2m−py ≥0
2m−pz ≥0

(µ = 1)

(111a)
eik(n2 )·Rl
= (ikx )px (iky )py
(±1)pz × . . .
4A2 α(n2 )
±

α(n2 )2 2 2
p
−
−η Rz
m
 z
α(n2 )
2
pz
(−η)
4η 2
√ Hm−1
± ηRz e
+ ...
...

2η
m (α(n2 ))m−pz π
m=1
. . . (α(n2 ))pz e±α(n2 )Rz erfc

α(n2 )
± ηRz
2η

ik(n3 )·Rl −α2 (n3 )/4η 2
px (ik py )(ik pz ) e
= (ikx )
y
z
A3 α2 (n3 )

(µ = 2)

(111b)

(µ = 3)

(111c)

Here, Ri is the projection of R along the i-axis, En (x) is the nth exponential integral and
Hn (x) is the nth Hermite polynomial. For µ = 1, it is assumed that the lattice lies along
the x-axis, and for µ = 2, it is assumed that the lattice lies in xy-plane.

A.3 Time-Domain Field of a Dipole
For convenience, the time-domain electric and magnetic fields emitted by a classical dipole
are presented. To this end, a nice reference is a 2004 paper by Nevels and Jeong [150], in
which authors provide a detailed derivation of the time-domain dyadic Green’s function for
electric and magnetic currents and fields using the propagator formalism. As only electric
145

dipole currents are germane to this dissertation, the relevant results are the following two
3 × 3 blocks of the full 6 × 6 dyad:
Ge,e (|r − r |, t − t ) = Ge,e (|R|, τ ) = −

δ (cτ − |R|)
4π|R|

I−

U (cτ − |R|) − U (−|R|)
4π|R|

(112)
Gh,e (|r − r |, t − t ) = Gh,e (|R|, τ ) = η

×

δ(cτ − |R|)
4π|R|

I

(113)

Here, the second subscript indicates the type of current generating the field, and the first
subscript indicates the type of field being generated, i.e., it is of interest to consider electric
currents generating electric and magnetic fields. Using this form of the dyadic Green’s
function, the free space electric and magnetic fields generated by a current of the form given
in Equation 4.82 can be rendered as:
0
E(r, t) =
−t
0
H(r, t) =
−t

dτ Ge,e (|Rs |, τ ) · ds js (t − τ )

(114)

dτ Gh,e (|Rs |, τ ) · ds js (t − τ )

(115)

Here, Rs = r − rs and it has been assumed that js (t) was off prior to t = 0.
Evaluating the integral in Equation 114, it is clear that E(r, t) can be separated into
farfield and nearfield components, i.e, fields that decays as |Rs |−1 , and those that decay
more rapidly:
E(r, t) = Ef f (r, t) + Enf (r, t)

(116)

The farfield contributions come from the term proportional to I in Ge,e , as well as the
146

portion of the

term that differentiates the Heaviside functions twice:
j (t − |Rs |/c)
ˆ
ˆ
Ef f (r, t) = Rs (ds · Rs ) − ds s
4π|Rs |

The nearfield contribution arises due to the remaining portion of the

ˆ
ˆ
Enf (r, t) = ds − 3Rs (ds · Rs )

(117)

term:

js (t − |Rs |/c) q(t − |Rs |/c)
+
4π|Rs |2
4π|Rs |3

(118)

t
Here, q(t) is the amount of charge deposited on the dipole at time t, q(t) =
0

dtjs (t).

Evaluating the integral in Equation 115, the expression for H(r, t) is as follows:

ˆ
H(r, t) = η ds × R

js (t − |Rs |/c) js (t − |Rs |/c)
+
4π|Rs |2
4π|Rs |

Here, the nearfield/farfield separation is obvious.

147

(119)

BIBLIOGRAPHY

148

BIBLIOGRAPHY

[1] I. Freestone, N. Meeks, M. Sax, and C. Higgitt, “The lycurgus cup–a roman nanotechnology,” Gold Bulletin, vol. 40, no. 4, pp. 270–277, 2007.
[2] S. Link and M. A. El-Sayed, “Optical properties and ultrafast dynamics of metallic
nanocrystals,” Annual Review of Physical Chemistry, vol. 54, no. 1, pp. 331–366, 2003.
[3] P. Mulvaney, “Surface plasmon spectroscopy of nanosized metal particles,” Langmuir,
vol. 12, no. 3, pp. 788–800, 1996.
[4] S. John, “Strong localization of photons in certain disordered dielectric superlattices,”
Physical review letters, vol. 58, no. 23, pp. 2486–2489, 1987.
[5] E. Purcell, “Spontaneous emission probabilities at radio frequencies,” Physical Review,
vol. 69, p. 681, 1946.
[6] T. F¨rster, “Zwischenmolekulare energiewanderung und fluoreszenz,” Annalen der
o
Physik, vol. 437, no. 1-2, pp. 55–75, 1948.
[7] D. Dexter, “A theory of sensitized luminescence in solids,” The Journal of Chemical
Physics, vol. 21, p. 836, 1953.
[8] I. L. Medintz and H. Mattoussi, “Quantum dot-based resonance energy transfer and its
growing application in biology,” Phys. Chem. Chem. Phys., vol. 11, no. 1, pp. 17–45,
2008.
[9] P. Michler, A. Kiraz, C. Becher, W. Schoenfeld, P. Petroff, L. Zhang, E. Hu, and
A. Imamoglu, “A quantum dot single-photon turnstile device,” Science, vol. 290,
no. 5500, pp. 2282–2285, 2000.
[10] A. Dousse, J. Suffczy´ski, A. Beveratos, O. Krebs, A. Lemaˆ
n
ıtre, I. Sagnes, J. Bloch,
P. Voisin, and P. Senellart, “Ultrabright source of entangled photon pairs,” Nature,
vol. 466, no. 7303, pp. 217–220, 2010.
149

[11] B. A. Munk, Frequency Selective Surfaces: Theory and Design. Wiley-Interscience,
2000.
[12] J. B. Pendry, “Negative refraction makes a perfect lens,” Physical Review Letters,
vol. 85, no. 18, pp. 3966–3969, 2000.
[13] P. Belov, R. Marques, S. Maslovski, I. Nefedov, M. Silveirinha, C. Simovski, and
S. Tretyakov, “Strong spatial dispersion in wire media in the very large wavelength
limit,” Physical Review B, vol. 67, no. 11, p. 113103, 2003.
[14] B. A. Munk, Metamaterials: critique and alternatives. Wiley-Interscience, 2009.
[15] A. R. Parker, “515 million years of structural colour,” Journal of Optics A: Pure and
Applied Optics, vol. 2, no. 6, pp. R15–R28, 2000.
[16] A. R. Parker and H. E. Townley, “Biomimetics of photonic nanostructures,” Nature
Nanotechnology, vol. 2, no. 6, pp. 347–353, 2007.
[17] M. Srinivasarao, “Nano-optics in the biological world: Beetles, butterflies, birds, and
moths,” Chem. Rev., vol. 99, pp. 1935–1961, 1999.
[18] J. Huang, X. Wang, and Z. Wang, “Controlled replication of butterfly wings for achieving tunable photonic properties,” Nano Letters, vol. 6, no. 10, pp. 2325–2331, 2006.
[19] D. Craig and T. Thirunamachandran, Molecular quantum electrodynamics. Dover Publications, 1998.
[20] G. D. Scholes, “Long-range resonance energy transfer in molecular systems,” Annual
review of physical chemistry, vol. 54, no. 1, pp. 57–87, 2003.
[21] S. Crooker, J. Hollingsworth, S. Tretiak, and V. Klimov, “Spectrally resolved dynamics
of energy transfer in quantum-dot assemblies: Towards engineered energy flows in
artificial materials,” Physical Review Letters, vol. 89, no. 18, p. 186802, 2002.
[22] T. Unold, K. Mueller, C. Lienau, T. Elsaesser, and A. D. Wieck, “Optical control of
excitons in a pair of quantum dots coupled by the dipole-dipole interaction,” Physical
review letters, vol. 94, no. 13, p. 137404, 2005.
[23] A. R. Clapp, I. L. Medintz, J. M. Mauro, B. R. Fisher, M. G. Bawendi, and H. Mattoussi, “Fluorescence resonance energy transfer between quantum dot donors and dye150

labeled protein acceptors,” Journal of the American Chemical Society, vol. 126, no. 1,
pp. 301–310, 2004.
[24] S. M. Dutra, Cavity quantum electrodynamics: the strange theory of light in a box,
vol. 10. Wiley-Interscience, 2004.
[25] H. T. Dung, L. Kn¨ll, and D.-G. Welsch, “Intermolecular energy transfer in the preso
ence of dispersing and absorbing media,” Physical Review A, vol. 65, no. 4, p. 043813,
2002.
[26] D. Dzsotjan, A. S. Sørensen, and M. Fleischhauer, “Quantum emitters coupled to
surface plasmons of a nanowire: A greens function approach,” Physical Review B,
vol. 82, no. 7, p. 075427, 2010.
[27] M. Minkov and V. Savona, “Radiative coupling of quantum dots in photonic crystal
structures,” arXiv preprint arXiv:1212.4960, 2012.
[28] Y. Chen, T. R. Nielsen, N. Gregersen, P. Lodahl, and J. Mørk, “Finite-element modeling of spontaneous emission of a quantum emitter at nanoscale proximity to plasmonic
waveguides,” Physical Review B, vol. 81, no. 12, p. 125431, 2010.
[29] D. Mart´
ın-Cano, A. Gonz´lez-Tudela, L. Mart´
a
ın-Moreno, F. Garc´
ıa-Vidal, C. Tejedor,
and E. Moreno, “Dissipation-driven generation of two-qubit entanglement mediated
by plasmonic waveguides,” Physical Review B, vol. 84, no. 23, p. 235306, 2011.
[30] T. Stievater, X. Li, D. Steel, D. Gammon, D. Katzer, D. Park, C. Piermarocchi,
and L. Sham, “Rabi oscillations of excitons in single quantum dots,” Physical Review
Letters, vol. 87, no. 13, p. 133603, 2001.
[31] C. Piermarocchi, P. Chen, L. Sham, and D. Steel, “Optical rkky interaction between charged semiconductor quantum dots,” Physical Review Letters, vol. 89, no. 16,
p. 167402, 2002.
[32] L. Duan and H. Kimble, “Scalable photonic quantum computation through cavityassisted interactions,” Physical Review Letters, vol. 92, no. 12, p. 127902, 2004.
[33] R. Esteban, T. Teperik, and J. Greffet, “Optical patch antennas for single photon
emission using surface plasmon resonances,” Physical Review Letters, vol. 104, no. 2,
p. 26802, 2010.
151

[34] A. Akimov, A. Mukherjee, C. Yu, D. Chang, A. Zibrov, P. Hemmer, H. Park, and
M. Lukin, “Generation of single optical plasmons in metallic nanowires coupled to
quantum dots,” Nature, vol. 450, no. 7168, pp. 402–406, 2007.
[35] A. Gonzalez-Tudela, D. Martin-Cano, E. Moreno, L. Martin-Moreno, C. Tejedor, and
F. Garcia-Vidal, “Entanglement of two qubits mediated by one-dimensional plasmonic
waveguides,” Physical Review Letters, vol. 106, no. 2, p. 20501, 2011.
[36] J. Hesthaven and T. Warburton, Nodal discontinuous Galerkin methods: algorithms,
analysis, and applications, vol. 54. Springer, 2008.
[37] B. Finlayson and L. Scriven, “The method of weighted residuals a review,” Appl.
Mech. Rev, vol. 19, no. 9, pp. 735–748, 1966.
[38] L. Evans, Partial differential equations, vol. 19. American Mathematical Society, 2010.
[39] G. Mur, “Absorbing boundary conditions for the finite-difference approximation of the
time-domain electromagnetic-field equations,” Electromagnetic Compatibility, IEEE
Transactions on, no. 4, pp. 377–382, 1981.
[40] J. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,”
Journal of computational physics, vol. 114, no. 2, pp. 185–200, 1994.
[41] L. Demkowicz and M. Pal, “An infinite element for maxwell’s equations,” Computer
methods in applied mechanics and engineering, vol. 164, no. 1, pp. 77–94, 1998.
[42] J. Jin and J. Volakis, “Te scattering by an inhomogeneously filled aperture in a thick
conducting plane,” Antennas and Propagation, IEEE Transactions on, vol. 38, no. 8,
pp. 1280–1286, 1990.
[43] J. Jin and J. Volakis, “Tm scattering by an inhomogeneously filled aperture in a thick
conducting plane,” in Microwaves, Antennas and Propagation, IEE Proceedings H,
vol. 137, pp. 153–159, IET, 1990.
[44] T. Hohage, F. Schmidt, and L. Zschiedrich, “Solving time-harmonic scattering problems based on the pole condition i: Theory,” SIAM journal on mathematical analysis,
vol. 35, no. 1, pp. 183–210, 2003.
[45] T. Hohage, F. Schmidt, and L. Zschiedrich, “Solving time-harmonic scattering problems based on the pole condition ii: Convergence of the pml method,” SIAM journal
on mathematical analysis, vol. 35, no. 3, pp. 547–560, 2003.
152

[46] T. Hohage and L. Nannen, “Hardy space infinite elements for scattering and resonance
problems,” SIAM Journal on Numerical Analysis, vol. 47, no. 2, pp. 972–996, 2009.
[47] R. Harrington, “Boundary integral formulations for homogeneous material bodies,”
Journal of electromagnetic waves and applications, vol. 3, no. 1, pp. 1–15, 1989.
[48] C. M¨ller, Foundations of the mathematical theory of electromagnetic waves. Springer,
u
1969.
[49] A. Taflove and S. Hagness, Computational electrodynamics: The finite-difference timedomain method, vol. 160. Artech house, 2000.
[50] B. Shanker, M. Lu, J. Yuan, and E. Michielssen, “Time domain integral equation
analysis of scattering from composite bodies via exact evaluation of radiation fields,”
Antennas and Propagation, IEEE Transactions on, vol. 57, no. 5, pp. 1506–1520, 2009.
[51] A. J. Pray, N. V. Nair, and B. Shanker, “Stability properties of the time domain electric
field integral equation using a separable approximation for the convolution with the
retarded potential,” Antennas and Propagation, IEEE Transactions on, vol. 60, no. 8,
pp. 3772–3781, 2012.
[52] E. Bleszynski, M. Bleszynski, and T. Jaroszewicz, “Aim: Adaptive integral method for
solving large-scale electromagnetic scattering and radiation problems,” Radio Science,
vol. 31, no. 5, pp. 1225–1251, 1996.
[53] J. Cooley and J. Tukey, “An algorithm for the machine calculation of complex fourier
series,” Mathematics of computation, pp. 297–301, 1965.
[54] A. Yilmaz, J. Jin, and E. Michielssen, “Time domain adaptive integral method for
surface integral equations,” IEEE Transactions on Antennas and Propagation, vol. 52,
no. 10, pp. 2692–2708, 2004.
[55] T. Eibert and J. Volakis, “Adaptive integral method for hybrid fe/bi modelling of 3-d
doubly periodic structures,” IEE Proceedings - Microwaves, Antennas and Propagation, vol. 146, no. 1, pp. 17–22, 1999.
[56] F. Wei and A. Yilmaz, “A hybrid message passing/shared memory parallelization of
the adaptiveintegral method for multi-core clusters,” Parallel Computing, vol. 37, no. 9,
pp. 279–301, 2011.
153

[57] L. Greengard and V. Rokhlin, “A fast algorithm for particle simulations,” Journal of
computational physics, vol. 73, no. 2, pp. 325–348, 1987.
[58] J. Song, C. Lu, and W. Chew, “Multilevel fast multipole algorithm for electromagnetic
scattering by large complex objects,” Antennas and Propagation, IEEE Transactions
on, vol. 45, no. 10, pp. 1488–1493, 1997.
[59] B. Shanker and H. Huang, “Accelerated cartesian expansions - a fast method for computing of potentials of the form r−ν for all real ν,” Journal of Computational Physics,
vol. 226, pp. 732–753, 2007.
[60] J. Jin, “The finite element method in electromagnetics,” 2002.
[61] H. Whitney, “Geometric integration theory,” 1957.
[62] J.-C. N´d´lec, “Mixed finite elements in R3 ,” Numerische Mathematik, vol. 35, no. 3,
e e
pp. 315–341, 1980.
[63] R. J. LeVeque, Finite volume methods for hyperbolic problems, vol. 31. Cambridge
university press, 2002.
[64] I. Perugia, D. Sch¨tzau, and P. Monk, “Stabilized interior penalty methods for the
o
time-harmonic maxwell equations,” Computer methods in applied mechanics and engineering, vol. 191, no. 41, pp. 4675–4697, 2002.
[65] A. Buffa, P. Houston, and I. Perugia, “Discontinuous galerkin computation of the
maxwell eigenvalues on simplicial meshes,” Journal of computational and applied mathematics, vol. 204, no. 2, pp. 317–333, 2007.
[66] B. Cockburn and C.-W. Shu, “Tvb runge-kutta local projection discontinuous galerkin
finite element method for conservation laws ii: general framework,” Math. Comp,
vol. 52, no. 186, pp. 411–435, 1989.
[67] J. Hesthaven and T. Warburton, “Nodal high-order methods on unstructured grids:
I. time-domain solution of maxwell’s equations,” Journal of Computational Physics,
vol. 181, no. 1, pp. 186–221, 2002.
[68] E. Montseny, S. Pernet, X. Ferri´res, and G. Cohen, “Dissipative terms and local timee
stepping improvements in a spatial high order discontinuous galerkin scheme for the
time-domain maxwells equations,” Journal of Computational Physics, vol. 227, no. 14,
pp. 6795–6820, 2008.
154

[69] L. Yuan and C.-W. Shu, “Discontinuous galerkin method based on non-polynomial
approximation spaces,” Journal of Computational Physics, vol. 218, no. 1, pp. 295–
323, 2006.
[70] L. Lin, J. Lu, L. Ying, and E. Weinan, “Adaptive local basis set for kohn-sham density
functional theory in a discontinuous galerkin framework i: Total energy calculation,”
Journal of Computational Physics, 2011.
[71] T. Hughes, G. Engel, L. Mazzei, and M. G. Larson, “A comparison of discontinuous and
continuous galerkin methods based on error estimates, conservation, robustness and
efficiency,” Lecture Notes in Computational Science and Engineering, vol. 11, pp. 135–
146, 2000.
[72] T. J. Hughes, G. Scovazzi, P. B. Bochev, and A. Buffa, “A multiscale discontinuous
galerkin method with the computational structure of a continuous galerkin method,”
Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 19, pp. 2761–
2787, 2006.
[73] A. D. Baczewski, D. L. Dault, and B. Shanker, “Accelerated cartesian expansions for
the rapid solution of periodic multiscale problems,” IEEE Transactions on Antennas
and Propagation, vol. 60, pp. 4281–4290, 2012.
[74] A. D. Baczewski, N. C. Miller, and B. Shanker, “Rapid analysis of scattering from
periodic dielectric structures using accelerated cartesian expansions,” JOSA A, vol. 29,
no. 4, pp. 531–540, 2012.
[75] A. Glisson and D. Wilton, “Simple and efficient numerical methods for problems of
electromagnetic radiation and scattering from surfaces,” Antennas and Propagation,
IEEE Transactions on, vol. 28, no. 5, pp. 593–603, 1980.
[76] A. Peterson, S. Ray, and R. Mittra, Computational Methods for Electromagnetics.
Wiley-IEEE Press, 1997.
[77] “H-field, e-field, and combined-field solutions for conducting bodies of revolution,”
¨
Arch. Elek. Ubertrag., vol. 32, no. 4, pp. 157–164, 1978.
[78] H. Contopanagos, B. Dembart, M. Epton, J. Ottusch, V. Rokhlin, J. Visher, and
S. Wandzura, “Well-conditioned boundary integral equations for three-dimensional
electromagnetic scattering,” Antennas and Propagation, IEEE Transactions on,
vol. 50, no. 12, pp. 1824–1830, 2002.
155

[79] R. Adams, “Physical and analytical properties of a stabilized electric field integral
equation,” Antennas and Propagation, IEEE Transactions on, vol. 52, no. 2, pp. 362–
372, 2004.
[80] F. Andriulli, K. Cools, H. Bagci, F. Olyslager, A. Buffa, S. Christiansen, and
E. Michielssen, “A multiplicative calderon preconditioner for the electric field integral
equation,” Antennas and Propagation, IEEE Transactions on, vol. 56, no. 8, pp. 2398–
2412, 2008.
[81] G. Hsiao and R. Kleinman, “Mathematical foundations for error estimation in numerical solutions of integral equations in electromagnetics,” Antennas and Propagation,
IEEE Transactions on, vol. 45, no. 3, pp. 316–328, 1997.
[82] R. Harrington, Time-Harmonic Electromagnetic Fields. Wiley-IEEE Press, 2001.
[83] K. Jordan, G. Richter, and P. Sheng, “An efficient numerical evaluation of the green’s
function for the helmholtz operator on periodic structures,” Journal of Computational
Physics, vol. 63, no. 1, pp. 222–235, 1986.
[84] S. Rao, D. Wilton, and A. Glisson, “Electromagnetic scattering by surfaces of arbitrary
shape,” IEEE Transactions on Antennas and Propagation, vol. 30, no. 3, pp. 409–418,
1982.
[85] R. Graglia, D. Wilton, and A. Peterson, “Higher order interpolatory vector bases for
computational electromagnetics,” IEEE Transactions on Antennas and Propagation,
vol. 45, no. 3, pp. 329–342, 1997.
[86] G. Golub, Matrix Computations. The Johns Hopkins University Press, 1996.
[87] Y. Saad, Iterative Methods for Sparse Linear Systems. Society for Industrial and
Applied Mathematics, 2003.
[88] D. Schaubert, D. Wilton, and A. Glisson, “A tetrahedral modeling method for electromagnetic scattering by arbitrarily shaped inhomogeneous dielectric bodies,” IEEE
Transactions on Antennas and Propagation, vol. 32, no. 1, pp. 77–85, 1984.
[89] R. Kleinman, G. Roach, and P. Van Den Berg, “Convergent born series for large
refractive indices,” JOSA A, vol. 7, no. 5, pp. 890–897, 1990.
156

[90] J. Rahola, “On the eigenvalues of the volume integral operator of electromagnetic
scattering,” SIAM Journal on Scientific Computing, vol. 21, no. 5, pp. 1740–1754,
2000.
[91] N. Budko and A. Samokhin, “Spectrum of the volume integral operator of electromagnetic scattering,” SIAM Journal on Scientific Computing, vol. 28, no. 2, pp. 682–700,
2006.
[92] J. Barnes and P. Hut, “A hierarchical O(n log(n)) force-calculation algorithm,” Nature,
vol. 324, pp. 446–449, 1986.
[93] L. Greengard and V. Rokhlin, “A fast algorithm for particle simulations,” Journal of
Computational Physics, vol. 73, no. 2, pp. 325–348, 1987.
[94] H. Cheng, L. Greengard, and V. Rokhlin, “A fast adaptive multipole algorithm in
three dimensions,” Journal of Computational Physics, vol. 155, pp. 468–498, 1999.
[95] M. Vikram and B. Shanker, “Fast evaluation of time domain fields in sub-wavelength
source/observer distributions using accelerated cartesian expansions (ace),” Journal of
Computational Physics, vol. 227, pp. 1007–1023, 2007.
[96] M. Vikram, A. Baczewski, B. Shanker, and L. Kempel, “Accelerated cartesian expansion (ace) based framework for the rapid evaluation of diffusion, lossy wave, and
klein-gordon potentials,” Journal of Computational Physics, vol. UPDATE, p. ., 2010.
[97] M. Vikram, H. Huang, B. Shanker, and T. Van, “A novel wideband fmm for fast integral equation solution of multiscale problems in electromagnetics,” IEEE Transactions
on Antennas and Propagation, vol. 57, pp. 2094–2104, 2009.
[98] L. Greengard and J. Huang, “A new version of the fast multipole method for screened
coulomb interactions in three dimensions,” Journal of Computational Physics, vol. 180,
pp. 642–658, 2002.
[99] N. Nishimura, “Fast multipole accelerated boundary integral equation methods,” Applied Mechanics Reviews, vol. 55, pp. 299–324, 2002.
[100] M. Vikram and B. Shanker, “An incomplete review of fast multipole methods - from
static to wideband - as applied to problems in computational electromagnetics,” ACES,
vol. 24, 2009.
157

[101] Y. Otani and N. Nishimura, “A periodic fmm for maxwell’s equations in 3d and its application to problems related to photonic crystals,” Journal of Computational Physics,
vol. 227, pp. 4630–4652, 2008.
[102] V. Rokhlin and S. Wandzura, “The fast multipole method for periodic structures,” in
Proceedings of the 1994 IEEE Antennas and Propagation Society International Symposium, 1994.
[103] H. Cheng, J. Huang, and T. Leiterman, “An adaptive fast solver for the modified
helmholtz equation in two dimensions,” Journal of Computational Physics, vol. 211,
pp. 616–637, 2006.
[104] K. Schmidt and M. Lee, “Implementing the fast multipole method in three dimensions,” Journal of Statistical Physics, vol. 63, pp. 1223–1235, 1991.
[105] K. Schmidt and M. Lee, “Multipole ewald sums for the fast multipole method,” Journal
of Statistical Physics, vol. 89, pp. 411–424, 1997.
[106] C. Berman and L. Greengard, “A renormalization method for the evaluation of lattice
sums,” Journal of Mathematical Physics, vol. 35, no. 11, pp. 6036–6048, 1994.
[107] C. Lambert, T. Darden, and J. Board, “A multipole-based algorithm for efficient calculation of forces and potentials in macroscopic periodic assemblies of particles,” Journal
of Computational Physics, vol. 126, pp. 274–285, 1996.
[108] S. Li and V. Lomakin, “Fast interpolation method for field evaluation in a periodic unit
cell,” in Proceedings of the 2010 IEEE Antennas and Propagation Society International
Symposium, 2010.
[109] S. Li, D. V. Orden, and V. Lomakin, “Fast periodic interpolation method for periodic unit cell problems,” IEEE Transactions on Antennas and Propagation, vol. 58,
pp. 4005–4014, 2010.
[110] Y. Shi and C. Chan, “Multilevel green’s function interpolation method for analysis of
3-d frequency selective structures using volume/surface integral equation,” JOSA A,
vol. 27, pp. 308–318, 2010.
[111] N. Chen, M. Lu, F. Capolino, B. Shanker, and E. Michielssen, “Floquet wave–based
analysis of transient scattering from doubly periodic, discretely planar, perfectly conducting structures,” Radio Science, vol. 40, 2005.
158

[112] P. Ewald, “Die berechnung optischer und elektrostatischer gitterpotentiale,” Annalen
der Physik, vol. 369, pp. 253–287, 1921.
[113] M. Warren and J. Salmon, “A parallel hased oct-tree n-body algorithm,” in Proceedings
of the 1993 ACM/IEEE conference on Supercomputing, 1993.
[114] J. Applequist, “Cartesian polytensors,” Journal of Mathematical Physics, vol. 24,
pp. 736–742, 1983.
[115] G. Kobidze, B. Shanker, and D. Nyquist, “Efficient integral-equation-based method for
accurate analysis of scattering from periodically arranged nanostructures,” Physical
Review E, vol. 72, 2005.
[116] T. Eibert, J. Volakis, D. Wilton, and D. Jackson, “Hybrid fe/bi modeling of 3-d doubly periodic structures utilizing triangular prismatic elements and an mpie formulation
accelerated by the ewald transformation,” IEEE Transactions on Antennas and Propagation, vol. 47, no. 5, pp. 843–850, 1999.
[117] M. Inoue, K. Ohtaka, and S. Yanagawa, “Light scattering from macroscopic spherical bodies. ii. reflectivity of light and electromagnetic localized state in a periodic
monolayer of dielectric spheres,” Physical Review B, vol. 25, no. 2, pp. 689–699, 1982.
[118] Y. Otani and N. Nishimura, “A periodic fmm for maxwell’s equations in 3d and its
applications to problems related to photonic crystals,” Journ, vol. 227, pp. 4630–4652,
2008.
[119] P. Vukusic and J. Sambles, “Photonic structures in biology,” Nature, vol. 424, pp. 852–
855, 2003.
[120] A. Parker and H. Townley, “Biomimetics of photonic nanostructures,” Nature Nanotechnology, vol. 2, pp. 347–353, 2007.
[121] E. Palik, ed., Handbook of Optical Constants of Solids. Academic Press, 1998.
[122] S. Xiao, U. Chettiar, A. Kildishev, V. Drachev, and V. Shalaev, “Yellow-light negativeindex metamaterials,” Optics Letters, vol. 34, no. 22, pp. 3478–3480, 2009.
[123] P. Johnson and R. Christy, “Optical constants of the noble metals,” Physical Review
B, vol. 6, no. 12, pp. 4370–4379, 1972.
159

[124] A. H. Mohammadian, V. Shankar, and W. F. Hall, “Computation of electromagnetic
scattering and radiation using a time-domain finite-volume discretization procedure,”
Computer Physics Communications, vol. 68, no. 1, pp. 175–196, 1991.
[125] S. Piperno, “Dgtd methods using modal basis functions and symplectic local timestepping,” European Journal of Computational Mechanics/Revue Europ´enne de
e
M´canique Num´rique, vol. 15, no. 6, pp. 643–670, 2006.
e
e
[126] J. Hesthaven, T. Warburton, J. Hesthaven, and T. Warburton, “High–order nodal
discontinuous galerkin methods for the maxwell eigenvalue problem,” Philosophical
Transactions of the Royal Society of London. Series A: Mathematical, Physical and
Engineering Sciences, vol. 362, no. 1816, pp. 493–524, 2004.
[127] M. Carpenter and C. Kennedy, “Fourth-order 2n-storage runge-kutta schemes,” tech.
rep., NASA, 1994.
[128] J. Niegemann, M. K¨nig, K. Stannigel, and K. Busch, “Higher-order time-domain
o
methods for the analysis of nano-photonic systems,” Photonics and NanostructuresFundamentals and Applications, vol. 7, no. 1, pp. 2–11, 2009.
[129] K. Busch, M. K¨nig, and J. Niegemann, “Discontinuous galerkin methods in nanophoo
tonics,” Laser & Photonics Reviews, vol. 5, no. 6, pp. 773–809, 2011.
[130] S. Wandzura and H. Xiao, “Symmetric quadrature rules on a triangle,” Computers &
Mathematics with Applications, vol. 45, no. 12, pp. 1829–1840, 2003.
[131] B. R. Vatti, “A generic solution to polygon clipping,” Communications of the ACM,
vol. 35, no. 7, pp. 56–63, 1992.
[132] L. Fezoui, S. Lanteri, S. Lohrengel, and S. Piperno, “Convergence and stability of a discontinuous galerkin time-domain method for the 3d heterogeneous maxwell equations
on unstructured meshes,” ESAIM: Mathematical Modelling and Numerical Analysis,
vol. 39, no. 06, pp. 1149–1176, 2005.
[133] T. Gruner and D. Welsch, “Green-function approach to the radiation-field quantization
for homogeneous and inhomogeneous kramers-kronig dielectrics,” Physical Review A,
vol. 53, no. 3, p. 1818, 1996.
[134] H. T. Dung, L. Kn¨ll, and D.-G. Welsch, “Resonant dipole-dipole interaction in the
o
presence of dispersing and absorbing surroundings,” Physical Review A, vol. 66, no. 6,
p. 063810, 2002.
160

[135] J. Wylie and J. Sipe, “Quantum electrodynamics near an interface,” Physical Review
A, vol. 30, no. 3, p. 1185, 1984.
[136] J. Wylie and J. Sipe, “Quantum electrodynamics near an interface. ii,” Physical Review
A, vol. 32, no. 4, p. 2030, 1985.
[137] G. Parascandolo and V. Savona, “Long-range radiative interaction between semiconductor quantum dots,” Physical Review B, vol. 71, no. 4, p. 045335, 2005.
[138] G. Parascandolo and V. Savona, “Long-range radiative interaction between semiconductor quantum dots,” Superlattices and Microstructures, vol. 41, no. 5, pp. 337–340,
2007.
[139] G. Tarel, G. Parascandolo, and V. Savona, “Ultralong-range radiative excitation transfer between quantum dots in a planar microcavity,” Physica Status Solidi (B), vol. 245,
no. 6, pp. 1085–1088, 2008.
[140] K. Xu and C. Piermarocchi, “Dynamics of elastic and inelastic energy transfer between
quantum dots in a microcavity,” Physical Review B, vol. 84, no. 11, p. 115316, 2011.
[141] A. Bossavit, “Solving maxwell equations in a closed cavity, and the question ofspurious
modes’,” Magnetics, IEEE Transactions on, vol. 26, no. 2, pp. 702–705, 1990.
[142] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini, “Unified analysis of discontinuous galerkin methods for elliptic problems,” SIAM journal on numerical analysis,
vol. 39, no. 5, pp. 1749–1779, 2002.
[143] A. Buffa and I. Perugia, “Discontinuous galerkin approximation of the maxwell eigenproblem,” SIAM Journal on Numerical Analysis, vol. 44, no. 5, pp. 2198–2226, 2006.
[144] L. Nannen, T. Hohage, A. Sch¨dle, and J. Sch¨berl, “High order curl-conforming hardy
a
o
space infinite elements for exterior maxwell problems,” arXiv preprint arXiv:1103.2288,
2011.
[145] M. L. Andersen, S. Stobbe, A. S. Sørensen, and P. Lodahl, “Strongly modified plasmonmatter interaction with mesoscopic quantum emitters,” Nature Physics, vol. 7, no. 3,
pp. 215–218, 2010.
[146] V. Papanicolaou, “Ewald’s method revisited: rapidly convergent series representations
of certain green’s functions,” Journal of Computational Analysis and Applications,
vol. 1, no. 1, pp. 105–114, 1999.
161

[147] P. Ewald, “Die berechnung optischer und elektrostatischer gitterpotentiale,” Annalen
der Physik, vol. 369, pp. 253–287, 1921.
[148] J. Walker, Fourier analysis. Oxford, 1988.
[149] M. Abramowitz and I. Stegun, Handbook of mathematical functions: with formulas,
graphs, and mathematical tables, vol. 55. Dover publications, 1965.
[150] R. Nevels and J. Jeong, “The time domain green’s function and propagator for
maxwell’s equations,” IEEE Transactions on Antennas and Propagation, vol. 52,
no. 11, pp. 3012–3018, 2004.

162