THE QUEST FOR ACTIVE MEDIA MODELS: A SELF-CONSISTENT FRAMEWORK FOR

SIMULATING WAVE PROPAGATION IN NONLINEAR SYSTEMS

By

Connor Adrian Glosser

A DISSERTATION

Michigan State University

in partial fulfillment of the requirements

Submitted to

for the degree of

Physics – Doctor of Philosophy

Electrical Engineering – Dual Degree

2018

THE QUEST FOR ACTIVE MEDIA MODELS: A SELF-CONSISTENT FRAMEWORK FOR

SIMULATING WAVE PROPAGATION IN NONLINEAR SYSTEMS

ABSTRACT

By

Connor Adrian Glosser

This work presents new approaches to simulations of active media at the level of individual
particles. Active systems contain internal, nonlinear, processes beyond those of simple
scattering systems; thus these new models afford high degrees of fidelity in exploring the
underlying physics without recourse to continuum or spatially-averaged approximations.
First, I examine the dynamics of microspheres set into motion by ambient acoustic
radiation in a fluid described by potential flow in the long-wavelength limit. Variations
in the local surface pressure caused by scattering from each microsphere set each micro-
sphere into motion following Newton’s second law. By expanding this pressure in terms
of spherical harmonics—natural eigenfunctions of the unretarded radiation kernel—I
recover an analytic description of the force on individual microspheres due to an incident
waveform. High-order numerical integrations then relate the surface potential on one
microsphere to the surface pressure on the others, thereby coupling the microspheres’
trajectories. These simulations predict a dominant translational effect along the direction
of propagation of the incident waveform, though they also reveal significant dipolar
interactions between microspheres that produce secondary expansions and contractions of
the collective microsphere system.

Extending my approach from acoustic to electromagnetic systems, I apply it to a
collection of quantum dots: “artificial” two-level atoms with a size-dependent energy
structure. The optical Maxwell-Bloch equations give the evolution of quantum dots under
the influence of electromagnetic fields; this evolution then produces secondary radiation
that couples a collection of quantum dots together. In my computational model, I cast
my secondary electromagnetic fields in terms of a point-to-point integral operator that

accurately recovers both near- and far-field effects. These fields, then, drive a set of implicitly
coupled Bloch equations (solved with an exponentially-fitted predictor/corrector scheme)
to give the dynamics of the system as a whole. In ensembles of up to 10 000 quantum dots,
my model predicts synchronized multiplets of particles that exchange energy, quantum
dots that dynamically couple to screen the effect of incident external radiation, localization
of the polarization due to randomness and interactions, as well as wavelength-scale regions
of enhanced and suppressed polarization.

The remainder of the work uses the same physical quantum dot system while moving
towards efficient computer-aided device design. I detail an improved propagation algorithm
to reduce the time and space complexity of the simulation dramatically, thereby facilitating
rapid analysis of promising device structures. The algorithm makes use of physical and
numerical approximations to effect large-scale calculations in reasonable CPU time. A
rotating-frame approximation removes high-frequency components in the evolution of
the system while simultaneously preserving accurate interference phenomena in space,
thereby affording far larger simulation timesteps. Additionally, projecting the source
current distribution onto a regular spatial grid makes use of a low-rank approximation
to the field propagator to communicate radiation information between distant groups of
particles via fast Fourier transforms in a manner reminiscent of fast multipole methods.

Copyright by
CONNOR ADRIAN GLOSSER
2018

I dedicate this thesis to my loving mom and dad, Cindy and Ron, and their eternal,

unconditional support.

v

ACKNOWLEDGMENTS

As I reflect on closing this chapter of my life, I can’t help but beam at the memories of
everyone who has helped shape me into the person I am today. Thanks in no small part to
my studies, I’ve gotten to meet many truly incredible people who guided me when I felt
lost, encouraged my grandiose schemes, or simply shared a portion of their story with me.
Some of you remain with me and some of you have left only marks, but I am a far better
person for having spent time with all of you and I offer my dearest thanks here.

First and foremost, I need to thank my incredible parents, Cindy Roach and Ron Glosser.
I truly don’t have enough words to express my gratitude for the love and support you’ve
given me in everything I’ve decided to pursue, academic or otherwise. Whereas school
taught me facts, you gave me what I needed to think and examine the world around me, no
doubt setting me on a lifelong journey of scientific and personal discovery on which this
Ph.D is but a stop.

Next, I’d like to thank my two wonderful advisers-turned-friends, Carlo Piermarocchi
and Shanker Balasubramaniam. I’m so happy my interests put me on a course to work with
both of you daily and I don’t think I could have asked for better mentors as I concurrently
learned to be both a physicist and an engineer. Together you’ve taught me a great deal
about doing science on a computer and a great deal more about determination and
professionalism. It saddens me that I have to leave MSU now that I’ve just started to come
into my own as a researcher, but it makes me much happier to think that we’ll find ways to
collaborate for a long time to come.

I’d also like to thank the members of my committee—Phil Duxbury, Stuart Tessmer,
John Albrecht, and John Luginsland—for your encouragement and insights into the work
I present here. Phil, I need to thank you in particular for giving my high-school self the
chance to experience what being a scientist is like. When we started working together
in 2008 (if you could call a high-schooler’s attempts at research “work”) I had no idea I

vi

was taking my first steps on a road that would lead me across a couple of continents and
through a Ph.D. I realize now I’ve learned so much about the world and the people that live
in it for having walked that road with you, and I sincerely hope I can pay your kindness
forward to many students of my own someday.

Xander van Kooten, Liz Fujita, Kim & Erin McGarrity, Scott O’Connor, and Bill Martinez,
thank you for being such noble and selfless companions. You all have shown me exceptional
warmth and kindness when I’ve needed it most, and every day you inspire me to learn, do,
or try something new and to become the best person I can be. While we may sometimes
live a great distance apart, I find tremendous joy in knowing our shared sense of adventure
will reunite us with striking regularity. You will always have a home with me.

Jenni & Giuseppe Zerilli, Dan & Marilyn Olds, Shannon Nicley, Chelsea Johnson, and
Gift Pimcharoen, I cashed in all my undergrad points for this Ph.D, but somehow I only
ended up with more undergrad points. Thank you for adopting me into Game Night
and showing me the real ropes of graduate school. Our evenings spent betraying each
other and hurling (false) accusations of betraying each other remain some of my favorite
memories of my time at MSU and I hope we can find time for numerous reunions soon.
Jos Thijssen, Chris Verzijl, Laurens Janssen, Félix Carlier, Joost de Gussem, Max
Huisman, Witek Nawara, Daan Kuitenbrouwer, and all of the other ICCP-ers, thank you for
sharing so much of your culture with this loudmouth American. Everything I ever needed
to know about grad school I learned in ICCP alongside you all, and the class simply would
not have existed without the many international exchanges you helped facilitate.

Finally, thanks to some of the friends I’ve had along the way: Charlie & JoAnn Steele,
Juan Manfredi, Tony Szedlak, Chris Sullivan, Robin de Kivit, Håkon Jensen, Julia Rohrer,
Louis Garcia, Mike Bennett, Chris Morse, and Steve Hughey. You all continue to brighten
my days with discussions of everything from the latest philosophies on doing science to
the most entertaining (and possibly least effective) game strategies.

vii

You all have routinely helped me find my way into truly magnificent experiences greater
than any I could have imagined, and for that I am profoundly grateful. I consider myself
beyond fortunate to count you amongst my friends and family and I sincerely hope our
bonds last long, long into the future. Thank you all for everything you have given me.

(cid:44)

viii

TABLE OF CONTENTS

LIST OF TABLES .

. .

LIST OF FIGURES . .

.

.

.

.

.

.

.

.

.

.

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

xi

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii

CHAPTER 1

.
.

.
.

INTRO .
.

1.1 Active media .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
.
1.1.1 Acoustophoresis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
.
1.1.2 Light in disordered media . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Numerical simulation of wave propagation . . . . . . . . . . . . . . . . . . .

.
.

.
.

1
1
2
3
4

.

.

.

. .

Introduction .

2.1
2.2 Continuum problem statement
2.3 Discretization of the Integral Equations
2.4 Analytic results .

CHAPTER 2 ACOUSTOPHORESIS IN RIGID MICROSPHERES . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . .

6
6
7
9
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
Single microsphere solution . . . . . . . . . . . . . . . . . . . . . . . . 12
. 14
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.5.1
. 16
2.5.2 Many-particle simulations . . . . . . . . . . . . . . . . . . . . . . . . . 18
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.4.1
2.4.2 Low-order interactions . . . . . . . . . . . . . . . . . . . . . . . . . .

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

2.5 Numerical Results .

Single microspheres

2.6 Conclusions .

. .

.

.

.

.

.

.

.

.

.

.
.

.
.

Introduction .

CHAPTER 3 RADIATION IN COUPLED QUANTUM DOT SYSTEMS . . . . . 23
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.1
. .
3.2 Problem Statement
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.3 Computational Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.4 Numerical Results .
Stability & adjacency effects . . . . . . . . . . . . . . . . . . . . . . . . 31
. . . . . . . . . . . . . . . . . . . . . . . . . 35
Inhomogeneous broadening . . . . . . . . . . . . . . . . . . . . . . . . 37
3.5 Conclusions & future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.4.1
3.4.2 Polarization enhancement
3.4.3

.

.

CHAPTER 4 AN FFT-ACCELERATED SOLUTION METHODOLOGY FOR

.

.

.
.

.
.
.

.
.
.

.
Introduction .

RADIATIVE MAXWELL-BLOCH SYSTEMS . . . . . . . . . . . . . 40
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
.
.
4.1 Abstract .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
. .
4.2
4.3 Problem statement
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.4 Computational approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.4.1 Expansion matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. 47
4.4.2 Complexity . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. 49
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
.
4.5.1 Convergence . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. 50

4.5 Results .

. .

. .

.
.
.

.

.

.

ix

4.5.2 Timing .
. .
4.5.3 Maxwell-Bloch .

. .

.
.

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

. 52
. 53

5.1 Conclusion .
.
5.2 Discussion . .

CHAPTER 5 CONCLUSIONS AND FUTURE PERSPECTIVE . . . . . . . . . . . 55
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
. 57
. 58
. 59

5.2.1 Advanced dissipators . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2 Micromagnetics .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Technical details . . . . . . . . . . . . . . . . . . . . . . . .

5.2.2.1

.
.
. .

.
.

.
.

.
.

. .

APPENDICES .
APPENDIX A
APPENDIX B
APPENDIX C
APPENDIX D

. .

. .

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

. .
. 61
QuEST Manual
. . . . . . . . . . . . . . . . . . . . . . . . . . . 62
The predictor/corrector integration scheme . . . . . . . . . . . 67
Chebyshev polynomials . . . . . . . . . . . . . . . . . . . . . . 72
Conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

BIBLIOGRAPHY .

. .

.

.

.

.

.

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

. 82

x

LIST OF TABLES

Table 2.1: Typical simulation parameters.

. . . . . . . . . . . . . . . . . . . . . . . . . 15

Table 3.1: Simulation parameters (unless otherwise stated); e and a0 denote the
elementary charge and Bohr radius. The decoherence times here,
while shorter than those typical of optical resonance experiments,
afford a shorter computational time but preserve dynamical emission
phenomena.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

. .

. .

Table 4.1: Power-law regressions and their coefficeints of determination (r2) for
the data in fig. 4.8. Here, t denotes the walltime of the simulation (in
seconds) and n the number of spatial basis functions.

. . . . . . . . . . . 53

xi

LIST OF FIGURES

Figure 1.1: Comparison of differential-equation and integral-equation methodolo-
gies. As both approaches have several advantages and disadvantages,
the choice of which technique to use largely depends on the problem
under consideration.

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

4

Figure 2.1: Coordinate notation.

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

. 10

Figure 2.2: Perpendicular configuration.

. . . . . . . . . . . . . . . . . . . . . . . . . 14

Figure 2.3: Translation of a single microsphere interacting with an incident pulse
( f0 (cid:3) 0.5 MHz, σ (cid:3) 7 µs). Microspheres interacting with the pulse
translate a finite distance along k due to the Gaussian envelope in eq. (2.9). 16

Figure 2.4: Confinement of non-interacting spheres to planes; identical counter-
propagating pulses ( f0 (cid:3) 20 MHz, σ (cid:3) 23.8 µs) initially displaced
along ˆz tend to align objects in ∇P (cid:3) 0 planes at λ/2 intervals. Here,
we have sampled the field and trajectories every 30 timesteps and
smoothed the resulting data with a 16-sample windowed average.

. . . 17

Figure 2.5: Calculated isopotential contours near a lone microsphere. Red and
blue colorations represent regions of positive and negative potential.
The motion of each microsphere through the background medium
serves primarily to produce a dipolar field of velocity potential with
vs serving as the sphere’s dipole moment.

. . . . . . . . . . . . . . . . . 18

Figure 2.6: Scaling behavior of two microspheres arranged perpendicularly to
an incident pulse for various radii and initial separations. The ( ,
) symbols on each axis denote data associated with that axis. The
follow a regression of ∆|v|d (cid:3) 0.250 754d
follow
∆|vmax|r (cid:3) 3.133 28 × 10−5a2.99814
. These trends strongly indicate dom-
inant dipolar interactions between microspheres.

−3.00077
12

0

, and the

. . . . . . . . . . . . . 19

Figure 2.7:

Isosurfaces of velocity potential (arb. units) calculated by evaluating
the ˆS and ˆD terms in eq. (2.8) for a N (cid:3) 16 particle simulation. Red,
blue, and yellow surfaces denote regions of positive, negative, and
zero potential respectively, with holes appearing due to intersections
with the bounding box. Rendered with VisIt[21].

. . . . . . . . . . . . . 20

xii

Figure 2.8: Fractional change in the volume of 20 randomly-initialized microsphere
clouds subject to the same incident pulse, smoothed with a 128-sample
moving average. Positive and negative values denote expansion and
contraction. σ (cid:3) 1.5 cm.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

Figure 3.1: Bloch sphere representation of a two-level quantum system. Equa-

tion (3.1) describes the evolution of(cid:12)(cid:12)ψ(cid:11) on the surface of this sphere in

response to an external electric field.

. . . . . . . . . . . . . . . . . . . . 25

Figure 3.2: Nonorthogonal and C0-continuous temporal basis function T(t/∆t)

constructed from intervals of third-order Lagrange polynomials.

. . . . 28

Figure 3.3: Flowchart of one iteration of the coupled quantum/electromagnetic

solution procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

Figure 3.4: Comparison of ρ00(t) using both fixed (requiring 59959 timesteps) and
rotating (requiring 160 timesteps) reference frames for two interacting
quantum dots. Both frames produce similar trajectories, however the
inset spy reveals the fixed frame contains a minute oscilliatory term
that the rotating frame does not.

. . . . . . . . . . . . . . . . . . . . . . . 32

Figure 3.5: Population dynamics of adjacent quantum dots (top, shown with the
trajectory of a quantum dot driven exclusively by the external laser
in black) and a 1024-particle neighborhood (dots (1) and (2) maintain
their separation and orientation between simulations). The interaction
between particles gives a greatly diminished response to the external
pulse through a dynamical detuning of the two-dot system. The
majority of the neighboring particles in the 1024-dot system follow
trajectories nearly identical to that prescribed by the laser (omitted
for clarity); many-dot effects, however, produce significant oscillatory
modes in the evolution of selected quantum dots (shown in gray). Note
that ρ(1)
00 appear to have two coherent modes in the presence of
multiple dots: a high-frequency oscillation between the pair, as well
as a low-frequency oscillation of the group about a decaying envelope.

Figure 3.6: Spatial distribution of(cid:12)(cid:12) ˜ρ01(cid:12)(cid:12) for 1024 quantum dots as an indicator of

00 and ρ(2)

33

polarization. Recorded at t (cid:3) 2 ps relative to the peak of a 1 ps-wide
π-pulse, the color and size of each sphere indicates the location of
each quantum dot and its polarization. Following a π-pulse, a single
quantum dot would have no remnant polarization; here, due to the
near-field interactions between particles, clusters of quantum dots
remain in highly-polarized states depending on their separation.

. . . . 34

xiii

Figure 3.7:

Inverse participation ratio (IPR) for the system shown in fig. 3.6. The
dashed line indicates the sample time for fig. 3.6. . . . . . . . . . . . . . . 34

Figure 3.8: Coloration of(cid:12)(cid:12) ˜ρ01(cid:12)(cid:12) as an indicator of(cid:12)(cid:12)˜P(cid:12)(cid:12) at t1 (cid:3) −0.05 ps, t2 (cid:3) 0 ps, t3 (cid:3)

0.05 ps, and t4 (cid:3) 0.10 ps relative to the peak of a 1 ps-wide pulse. 10 000
quantum dots randomly distributed throughout a 0.2 µm (radius) ×
4 µm cylinder oriented along kL demonstrate near-field the effects of
fig. 3.5 as distinct, outlying bright/dark quantum dots. Additionally,
the size of the system allows for wavelength-scale phenomena that
appear here as five standing regions of enhanced polarization. (Note
that we model quantum dots as point objects; the size of the spheres
here has no physical interpretation.)

Figure 3.9: Coloration of(cid:12)(cid:12) ˜ρ01(cid:12)(cid:12) for 10 000 quantum dots arranged in a finite planar
Figure 3.10: z-distribution of polarization(cid:12)(cid:12) ˜ρ01(cid:12)(cid:12) for the geometry in fig. 3.8. In each

geometry (in units of λ). The slab displays a prominent polarization
pattern 1.25 ps after the peak of a 1 ps pulse.

. . . . . . . . . . . . . . . . . . . . . 36

. . . . . . . . . . . . . . . . 37

Figure 4.1:

Figure 4.2:

simulation, 10 000 quantum dots had random Gaussian noise (with
width parameters of 0.0 meV, 0.2 meV, 0.4 meV and 0.6 meV) added to
a resonant ω0 (cid:3) 1500 meV. The induced long-range patterns in the
polarization remain for mild detunings (cid:46) 0.5 meV but have completely
washed out at 0.6 meV. Additionally, each simulation displayed the
characteristic near-field coupling effects of fig. 3.5 (not visible here).
the interaction matrix for g(r, t) (cid:3) |r|−1 between points randomly
distributed throughout a unit cube appears to have very little structure
(left). while we cannot fully order the points in three dimensions,
we can permute the matrix so as to give points within the same
neighborhood adjacent indices. points ordered thusly produce a
diagonally-dominant interaction matrix with a hierarchical Toeplitz
substructure (right). such structures indicate portions of the matrix
have very accurate low-rank approximations that offer the possibility
of significant compression.

. . . . . . . . . . . . . . . . . . . . . . . . . . 43

. . 38

Illustration of the grid structure and related terminology. All of the
sources within a box (shown as the central shaded square) map to the
same set of expansion points (shown as open circles) indexed relative
to rbox.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

.

.

.

.

.

.

xiv

Figure 4.3:

Figure 4.4:

Illustration of the nearfield criterion for a third order expansion (γ (cid:3) 2).
The dashed line indicates the complete nearfield of the box associated
with r0—i.e. all boxes that have an expansion point within ∆s of the
expansion around r0. Consequently, all of the s(cid:96)(r) within the central
dark blue square have a pairwise interaction with the s(cid:96)(cid:48)(r) inside the
dashed box.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

.

.

.

.

.

Illustration of nearfield corrections between close boxes. Expansions
A and B overlap, but only box B lies in the nearfield of box C. The
grid-based propagation strategy only remains accurate for distant
source/observer pairs. To avoid incurring undue error, we remove
the interaction “through the grid” between the BC pair (red line) and
replace it with a more accurate “direct” interaction (dashed blue line).
The AC pair requires no such treatment as they have well-separated
expansion regions.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

Figure 4.5: Spatial expansion pattern for s(cid:96)(r) (cid:3) ˆd(cid:96)δ(r − r(cid:96)) through order five.
The numbers on the grid indicate the points added in the equivalent
order expansion such that an expansion of order m requires (m + 1)3
grid points. Moreover, increasing the expansion order incorporates
new points in such a way as to keep s(cid:96)(r) as close to the center of the
expansion region as possible.

. . . . . . . . . . . . . . . . . . . . . . . . . 48

Figure 4.6: Polynomial interpolation of g(x − xsrc) near xobs. Here, the green
curve represents the actual g(x − xsrc) and the dashed black line its
approximation. Evaluating the mth-order approximation requires
samples of the signal at m + 1 grid points surrounding xobs.

. . . . . . . 51

Figure 4.7: (cid:96)2 error calculation of g(r − r(cid:48)) (cid:3) e ik|r−r(cid:48)|/|r − r(cid:48)| for expansion orders
zero through four. For each of the points above, a source and obser-
vation box separated by ∆r (cid:3) 10λˆx + 10λ ˆy + 10λˆz each contain 64
randomly placed points, thus the points within each box all map to
identical nodes in the auxiliary grid. The moment-matching expansion
scheme dictates that the overall error for an expansion of order m
scales as O(∆sm+1).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Figure 4.8: Time to evaluate retarded potential vs. number of spatial basis functions
for several system sizes. Simulated on an Intel(R) Core(TM) i5-4690K
CPU @ 3.50 GHz.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

.

.

xv

Figure 4.9: Comparison of a Maxwell-Bloch quantum dot simulation with direct
(red) and FFT-accelerated (blue) propagation strategies. Here, we show
the off-diagonal matrix elements of 25 random quantum dots (out of
1024) as a function of time. The quantum dots experience an incident
π-pulse and their interactions induce significant secondary oscillations
(see the second-to-last row). The FFT-accelerated algorithm developed
here readily captures these effects, making it eminently suitable for
large-scale quantum dot simulations.

. . . . . . . . . . . . . . . . . . . . 54

Figure A.1: Hierarchy of the major pieces of QuEST. While a true class diagram
would depict the inheritance relationships of the objects that comprise
the QuEST library, These blocks roughly indicate the interdependence
of the program pieces that form QuEST, akin to a class hierarchy
(though avid software developers should not take this literally).

If ∂tx•(t) (cid:3) f(cid:8)(cid:0)(cid:9)x•(t − |r••|) + x•(t − |r••|))

Figure B.1: Illustration of the “retardation problem” in solving coupled delay
differential equations.
then determining x•(t − |r••|) and x•(t − |r••|) requires some sort of
interpolation scheme. Evaluating ∂tx•(t) in hopes of obtaining x•(t)
therefore requires knowledge of x•(t). The converse holds as well;–
evaluating ∂tx•(t) requires knowledge of x•(t) which poses problems
in evaluating the derivatives simultaneously. The predictor step of the
integrator overcomes this issue by providing a stable extrapolation of
past quantities to provide a guess value of x•(t) and x•(t).

. . . . . . . . 69

. . . . . 62

Figure B.2: Selection of the λn parameters used in determining the predictor/cor-

rector coefficients.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

Figure C.1: Chebyshev polynomials T0(x) through T6(x). The roots of Tm+1(x)
give the sample points for an interpolation scheme of order m on the
interval [−1, 1].

Figure C.2: Interpolation of a pathological function, f(x) (cid:3)(cid:0)1 + 16x2(cid:1)−1/2, using

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

.

.

Chebyshev and equally-spaced Lagrange interpolation schemes. Be-
cause the Chebyshev scheme clusters sample points in the tails of
the interpolation window, it does not suffer from high-order ringing
effects (Runge’s phenomenon).

. . . . . . . . . . . . . . . . . . . . . . . . 74

xvi

CHAPTER 1

INTRO

“Begin at the beginning,” the King said, very gravely, “and go on till you come to the
end: then stop.”

—The King of Hearts, Alice in Wonderland

As our ability to resolve physical systems grows, so, too, do our demands of the
models we build. Unfortunately, we rapidly approach an age where approximate analytic
descriptions of physical systems simply cannot recover important features in systems of
interest. Substituting numerical—i.e. in silico—techniques into higher resolution models
alleviates this problem somewhat (particularly as computer hardware evolves “below”
the model), though the computational overhead of naïve algorithms limits their utility
beyond nearly analytic calculations, particularly in multi-scale or multiphysics applications.
Computational scientists, then, have a tricky job: they must have expertise in their own
domain so as to begin formulating testable hypotheses about the (multi-)physical systems,
as well as expertise in “computation” to understand the constraints inherent in exploring
such hypotheses numerically.

1.1 Active media

Active media—dynamical systems that couple together through radiative processes—
encompass a large class of problems across length and time scales and exhibit behavior
fundamentally different from their passive counterparts. These behaviors offer new avenues
towards novel device designs, though engineering devices to exploit them presents a
considerable challenge. In particular, the underlying dynamical nature of active systems
makes them inherently nonlinear, and accounting for these dynamics alongside the radiative
coupling requires a multiphysical approach. As such, this thesis primarily concerns itself
with the development of computational techniques for use in modeling active systems both

1

efficiently and accurately with a particular focus on acoustophoresis (acoustically-induced
motion) and quantum optics.

1.1.1 Acoustophoresis

The placement of biological material or microscopic objects plays a large role in a number
of technological applications related to biosensing [78]. For example, droplet-based
microfluidics entails the careful construction, manipulation, and sensing of µL (or smaller)
droplets of fluid for myriad lab-on-a-chip applications. Similarly, DNA microarrays employ
large numbers of micron-sized spots—each filled with a particular DNA sequence—to
perform an experiment on thousands of genes simultaneously. Manipulating individual
droplets and fabricating a microarray both require careful placement of apparatus; though,
given the small length scale of these problems, conventional manipulation techniques
(such as pipetting) become problematic. Instead, researchers require alternative means of
precisely controlling the motion of such objects.

Because of their low energy and minimal invasiveness, so-called acoustic tweezers offer
an ideal mechanism for precise microscopic positioning in both biological and nonbiological
applications. By varying the pressure in the surrounding environment through the precise
application of ultrasonic pulses, acoustic tweezers effect acoustophoresis for a variety of
purposes including object placement, fractionation [59], and flow shaping. Present models
of such systems often fail to account for multiple scattering events or do so only by way of
time- or volume-averaged techniques. This poses a significant problem in extremely small
systems where such effects necessarily become important. To overcome this, chapter 2
develops the machinery to simulate the motion of an ensemble of hard spheres within a
fullwave scattering framework.

2

1.1.2 Light in disordered media

Disordered nanostructures give rise to striking optical phenomena. Cyphochilus spp. beetles,
for instance, owe their intense white coloration to an aperiodic arrangement of scales that
scatters light across wavelengths [80]. Whereas man-made objects of similar whiteness and
intensity require structures as thick as 100 µm, the scales of cyphochilus spp. can span as
few as 4 µm by exploiting an interior structure of disordered filaments. Additional effects
such as Anderson localization [5]—wherein light remains spatially confined for extended
periods of time—and Lévy flights [10]—in which light experiences superdiffusive, i.e.
superlinear (in time), behavior—all follow from similarly disordered structures. Because
of their dynamical absorption and emission effects, these phenomena all fall under the
broad category of active light/matter interactions.

Quantum dots—colloquially called “tuneable atoms” because of the ease with which
engineers can manipulate their spectra—stand as a premier element in modern optoelec-
tronic devices and have enjoyed widespread success in areas such as light harvesting and
quantum information storage/retrieval. Because of this, a number of analytic methods
based on ensemble-averaged single particle and higher order Green’s functions have arisen
to describe light propagation in such systems. These methods essentially cast the active
quantum system as a homogeneous continuum source in the governing electromagnetic
equations (Maxwell’s equations for semiclassical systems) [6], though such a treatment
cannot account for effects arising from long-range order, partial order, or anisotropies.
Models of these effects, then, require a detailed computational framework that captures
scattering/radiation effects between discrete, disordered elements. In chapter 3 we develop
just such a framework; one that eschews a phenomenological treatment of matter, prefer-
ring instead to “generate” the constitutive relations that govern the sources in Maxwell’s
equations.

3

Differential
equation
methods

Properties:

± Local operators
+ Straightforward mathematics
+ Incorporates material properties
- “Volumetric” discretization
- Artificial radiation boundary condition

Properties:

Propagation
techniques

Integral
equation
methods

± Global operators
-
- Exclusively models free/homogeneous

Intricate mathematics

space

+ “Surface” discretization
+ Automatic radiation boundary condition

Figure 1.1: Comparison of differential-equation and integral-equation methodologies. As
both approaches have several advantages and disadvantages, the choice of which technique
to use largely depends on the problem under consideration.

1.2 Numerical simulation of wave propagation

Wave phenomena influence every aspect of our lives across every scale of the universe,
thus it should come as little surprise that techniques to explore and analyze these phe-
nomena lie at the heart of many scientific and engineering disciplines. Broadly, these
techniques fall into differential-equation-based and integral-equation-based categories
with each having several benefits and drawbacks. Differential-equation-based models have
a large, sparse character that follows directly from a discretization of the wave equation
throughout a volume of interest. Such a volumetric discretization explicitly describes
fields within the entirety of the simulation domain and works well to capture cavity
effects or effects involving variations in material parameters. Open systems where waves
propagate “to infinity” require significant care, however, and employ corrections such
as perfectly-matched layers [11] or absorbing boundary conditions [57] to approximate
correct asymptotic behavior.

Integral-equation-based methods, on the other hand, remedy many of the problems

4

associated with differential methods. By employing a Green’s function to directly deter-
mine the response of interacting sources, these methods automatically capture radiation
boundary conditions and eliminate the numerical issues associated with intermediate
meshes entirely. These features, however, do not come without cost. Integral-equation
methodologies have a great deal more mathematical sophistication than their differential-
equation analogues and their discretizations inherently produce dense interaction matrices.
The O(n2) solution cost of these matrices scales poorly for systems with many degrees
of freedom, thus a large body of work has emerged to reduce this cost to O(np) (where
p < 2) or even O(n log n). These so-called “fast methods” exploit the underlying structure
of the interaction matrix to develop compressed representations of the interaction between
distant sources with fully-controllable numerical error. Because of its accuracy and speed,
the development of such a method stands as the principal contribution of chapter 4.

5

CHAPTER 2

ACOUSTOPHORESIS IN RIGID MICROSPHERES

“Where are you going?”
“Which way should I go?”
“That depends on where you are going.”
“I don’t know.”
“Then it doesn’t matter which way you go.”

—The Cheshire Cat and Alice, Alice in Wonderland

2.1 Introduction

Computational approaches that employ an integral equation formalism to examine
acoustic scattering from particles typically assume a static environment in which scatterers
remain stationary. At present, a large body of work details such scattering problems [81,
24, 83]. While these stationary integral equation methods offer a large degree of accuracy
in capturing the underlying physics, many problems of interest require a fully dynamical
treatment. For instance, in biomedical physics, gas-filled microspheres exposed to ultrasonic
beams have demonstrated effectiveness as a contrast imaging agent [13] and as a drug
delivery method [3, 42], and Ding et al. have demonstrated their manipulation using
acoustic tweezers in microfluidic channels [23]. Moreover, composite materials consisting
of colloidal in-fluid suspensions have peculiar sound propagation properties that can
deviate from the ones of homogeneous liquids [82]. In each of these applications, the
unconstrained motion of scatterers requires a self-consistent description of their dynamics
in conjunction with a description of the acoustic field propagation.

Here, we demonstrate the applicability of coupling particle kinetics to a time-domain
integral equation (TDIE) scattering framework to model rigid-sphere motion induced by a

6

time-dependent acoustic potential. Specifically, we consider the case of an acoustic pulse
acting on microspheres that move in a fluid. Effective Langevin time-averaged radiation
pressure forces [50, 17], which consider the case of a steady radiation flux incident upon
a body kept in static equilibrium, do not provide an appropriate model in this case as
they cannot accommodate inter-particle scattering effects. While many theoretical and
computational descriptions of higher-order acoustic interactions exist [40, 25, 26, 44, 8],
few actually make use of computed fields to predict particle trajectories. As we consider
only short-duration pulses, we refrain from time-averaging in favor of using a time-
domain scattering formulation to explicitly calculate particle trajectories resulting from
a prescribed pulse. By adopting a weakly-compressible potential formulation of the
fluid media, our scalar wave problem inherits a number of similarities and solution
techniques from scattering problems in electromagnetic theory, a topic previous works
discuss extensively [77, 40, 54]. Moreover, our time-domain formulation readily allows the
study of transient phenomena (such as acoustic tweezing), a convenience not shared with
more common frequency domain approaches.

We structure the remainder of this chapter as follows: we first provide a formal
mathematical description of the problem—including details on both the kinetic and field
methods—followed by data obtained from various pulse and microsphere configurations,
demonstrating both attractive and repulsive regimes suitable for subtle control of spherical
systems in a homogeneous fluid. Finally, we offer concluding remarks on the effectiveness
of the simulation as well as our thoughts on possible future extensions.

2.2 Continuum problem statement

Consider a collection of N rigid, non-intersecting spherical scatterers (microspheres),
each having radius ak, position rk, and enclosing volume Vk ⊂ R3. The microspheres
move in a homogeneous exterior fluid occupying VE, where we denote the boundary of
each microsphere as Ωk (cid:3) ∂Vk and thus may ascribe to each an outward-pointing normal

7

(cid:0)θ, φ(cid:1), where θ and φ represent colatitude and azimuthal angles with respect to the

ˆnk
local origin (microsphere center). We wish to investigate the reaction of the system to
an incident acoustic pulse, thus the fluid carries a prescribed (band-limited) waveform
through the microsphere system in which it interacts with each of the ∂Ωk according to the
“sound-hard” regime presented in [54]. The incident acoustic pulse, in combination with
the acoustic field scattered from each microsphere and the hydrodynamic field induced
by the relative velocity of each microsphere, acts as a perturbation to the initially at-rest
uniform ideal fluid [58, 53]. We consider here the linear regime, in which the perturbation
induced by the acoustic and aerodynamic contribution remain sufficiently small so that the
velocity field v(r, t) satisfies the condition |v(r, t)| (cid:28) cs, where cs represents the speed of
sound in the fluid. In this limit, the velocity potential, defined by v(r, t) (cid:3) ∇ϕ(r, t), satisfies
the scalar wave equation:

ϕ(r, t) (cid:3) 0,

(2.1)

(cid:19)

(cid:18)∇2 − 1

∂2
∂t2

c2
s

and we may express the pressure perturbation at any point in the exterior medium as

p(r, t) (cid:3) −ρ0

∂ϕ(r, t)

∂t

,

(2.2)

where ρ0 denotes the equilibrium density of the fluid. Rigidity of the Ωk necessarily
prescribes boundary conditions on the normal velocity components at each interface,
namely,

· ˆnk.

(2.3)

(cid:12)(cid:12)(cid:12)(cid:12)r∈Ωk

∂ϕ(r, t)
∂ ˆnk

drk
dt

(cid:3)

where rk represents the center-of-mass coordinate for the kth microsphere.

Using these relations, we apply the Kirchoff-Helmholtz theorem to define the following

(cid:18)

N−1


i(cid:3)0

system of integral equations,
(cid:120)
ϕ(r, t) (cid:3) ϕinc(r, t)+
(cid:48) , (2.4)
where A(cid:48) ∈ Ωk(t(cid:48)) and gr(r, t; r(cid:48), t(cid:48)) denotes the Green’s function for a retarded potential,
(2.5)

δ(t − t(cid:48) − |r − r(cid:48)|/cs)

(cid:48)) ∂gr(r, t; r(cid:48), t(cid:48))

(cid:48)) ∂ϕ(r(cid:48), t(cid:48))
∂ ˆnk

− gr(r, t; r(cid:48), t

dA

(cid:48) dt

ϕ(r(cid:48), t

∂ ˆnk

gr(r, t; r(cid:48), t

(cid:48)) (cid:3)

4π|r − r(cid:48)|

.

(cid:19)

8

If the system remains localized to a region with small dimensions when compared to the
wavelength of sound, retardation effects become negligible and we may instead use the
Laplace-kernel Green’s function,

To ease notation, we define the following two integral operators,

ˆSk
ˆDk

reducing eq. (2.4) to

g(r, r(cid:48)) (cid:3)

1

4π|r − r(cid:48)| .

Þ
(cid:2)ϕ(r ∈ Ωk(t), t)(cid:3) (cid:3)
Þ
(cid:2)ϕ(r ∈ Ωk(t), t)(cid:3) (cid:3)
(cid:16) ˆDk − ˆSk
N−1


ϕ(r, t) (cid:3) ϕinc +

k(cid:3)0

g(r, r(cid:48)) ∂ˆnk
ϕ(r(cid:48), t) dA
ϕ(r(cid:48), t) ∂ˆnk g(r, r(cid:48)) dA ,

(cid:17)(cid:2)ϕ(r ∈ Ωk(t), t)(cid:3) .

(2.6)

(2.7a)

(2.7b)

(2.8)

In solving eq. (2.8), we obtain the velocity potential everywhere for a given time without
retarded scattered fields. For the incident pulse, ϕinc, we consider superpositions of wave
packets of the form

ϕinc(r, t) (cid:3) P0 cos(ω0t − k · r)e

−(cs t− ˆk·r)2/(2σ2)

(2.9)

due to their established spectral properties (i.e. known bandwidth and center frequency).
Finally, the variation in pressure (and thus ϕ) over each of the Ωk necessarily propels each
microsphere according to

Þ ∂ϕ(r, t)

(2.10)

d2rk
dt2 (cid:3) ρ0

mk

dS
where S denotes the outward-pointing normal of Ωk(t).
2.3 Discretization of the Integral Equations

∂t

To solve the integral equation scattering problem, we begin by discretizing our field in
both space and time. As we have restricted our particles to completely spherical geometries,

9

ak

ˆz

rk

Ωk

djk

rj
ˆy

a j

ˆx

Figure 2.1: Coordinate notation.

the Y(cid:96)m(θ, φ) spherical harmonics give simple eigenfunctions of the operators in eq. (2.7)
where

(cid:0)θ, φ(cid:1) (cid:3)

Y(cid:96)m

(cid:115)2(cid:96) + 1

4π

((cid:96) − m)!
((cid:96) + m)!Pm

(cid:96) (cos θ)e imφ.

As a result, they lend themselves well to an expansion of ϕ on the surface of each
microsphere with respect to the microsphere’s center,







(cid:96)(cid:62)0

|m|(cid:54)(cid:96)

ϕ(r ∈ Ωk , t) (cid:3)

(t) Y(cid:96)m

Ck
(cid:96)m

(cid:0)θ, φ(cid:1).

(2.11)

(2.12)

By considering eq. (2.2) and expressing the local velocity potential at each of the Ωk as a
linear combination of spherical harmonics, we have a complete representation of the body
force acting on each microsphere due to the orthogonality of dipole terms with the rest of
the multipoles. Hence,

body(t) (cid:3) −Þ
(cid:114)2π
3 r2(cid:16)(cid:2) (cid:219)Ck11(t) − (cid:219)Ck1−1(t)(cid:3) ˆx + i(cid:2) (cid:219)Ck11(t) + (cid:219)Ck1−1(t)(cid:3) ˆy − √

p(r, t) dS

(cid:17).

2 (cid:219)Ck10(t)ˆz

(2.13)

(cid:3) ρ0

Fk

The problem then becomes one of solving a system of linear equations that we may

compactly represent as

Z · ϕ (cid:3) F ,

10

(2.14)

Þ

with the overbar denoting a matrix quantity. We define the elements of F as projections of
the incident field onto local spherical harmonics,

F k
(cid:96)m (cid:3)

∗
(cid:96)m

Y

(θ, φ)ϕinc(r, t) dA ,

(2.15)

and detail Z jk
(cid:96)m,(cid:96)(cid:48)m(cid:48) for two cases: j (cid:3) k and j (cid:44) k. In the instances where j (cid:3) k, eq. (2.7)
propagates effects of the interaction to every point on a surface sharing a coordinate system
with the original, thus the harmonics remain orthogonal and

Z j j
(cid:96)m,(cid:96)(cid:48)m(cid:48) (cid:3)

(cid:96) + 1
2(cid:96) + 1 δ(cid:96)(cid:96)(cid:48)δmm(cid:48)




(cid:96),m

g(r, r(cid:48)) (cid:3)

1

2(cid:96) + 1

r (cid:96)
<
r (cid:96)+1

>

Y(cid:96)m

(cid:0)θ, φ(cid:1)Y

∗
(cid:96)m

(cid:0)θ(cid:48), φ(cid:48)(cid:1)

(2.16)

(2.17)

after exploiting the well-known expansion theorem for eq. (2.6),

where r< (cid:3) min(|r|, |r(cid:48)|) and r> (cid:3) max(|r|, |r(cid:48)|). A description of the off-diagonal terms
where j (cid:44) k proceeds much the same way, though the surface expansions no longer
share a local origin, complicating the projections. Translation operators for the spherical
harmonics [19, 37] allow analytic expressions for these matrix elements, though we eschew
such operators in favor of numerical integration for speed.

Thus, at every timestep of the simulation, the algorithm proceeds as follows: (i) project
the incident pulse and surface velocities onto local expansions of spherical harmonics,
(ii) propagate scattering effects through space by inverting the operators in eq. (2.8),
(iii) project these scattered fields onto local spherical harmonics to give a total representation
of ϕ on each surface, and (iv) move each microsphere according to eq. (2.10) and advance
t → t + ∆t. Only (cid:96) (cid:3) 1 terms contribute to the center-of-mass motion of rigid microspheres,
thus we use the C1m coefficients in evolving eq. (2.10).

The inversion in step (ii) above requires some care; by simply inverting the entire
ˆD − ˆS, to give a single surface pressure, eq. (2.10) reduces to a

propagation operator,
differential equation of the form

(cid:219)rk (cid:3) f(t, rk , (cid:219)rk).

11

(2.18)

This presents a number of irregularities with conventional integration schemes and will
rapidly diverge towards ±∞ due to the additional (cid:219)rk on the right if implemented naïvely.
To remedy this, we note that ˆS serves to produce only a reaction or drag term on each
ˆD and
microsphere that impedes motion. By maintaining quantities for the inversion of
ˆS separately, we remove the explicit dependence on (cid:219)rk by introducing a linear coefficient
in the form of an additional mass term—given by the (cid:219)rk-dependent contribution in the
single-layer ˆS operator—when solving eq. (2.10).

2.4 Analytic results

2.4.1 Single microsphere solution
As an example, consider a single sphere of density ρs and radius a. Taking ka (cid:28) 1, we
may approximate eq. (2.9) as ϕinc(r, t) (cid:3) v0(t)z and we wish to find the response velocity
of the sphere, u, in terms of the field velocity v (cid:3) ∇ϕinc. It follows that the expansion of
ϕinc contains only (cid:96) (cid:3) 1 terms, thus

ϕinc (cid:3) v0(t) a cos(θ)

on the surface of the sphere. Similarly, from eq. (2.3),

∂ˆnϕ (cid:3) u · ˆn

(cid:3) uz a cos(θ)

(2.19)

(2.20)

due to the symmetries present in x and y. As a result,

ϕ −Þ

(cid:3) v0a cos(θ) −Þ

ϕ(r(cid:48)) ∂ˆn(cid:48)G(r, r(cid:48)) dS

(cid:48)

auz cos(θ) G(r, r(cid:48)) dS

(cid:48) ,

(2.21)

and it becomes apparent that only (cid:96) (cid:3) 1, m (cid:3) 0 terms in eq. (2.17) remain after integrating.
Consequently, the field becomes

cos(θ)

(2.22)

(cid:18)

ϕ(r, t) (cid:3)

(cid:19)

v0(t)|r| +

a3(v0(t) − uz)

2|r|2

12

outside the microsphere and

ϕ(r ∈ Ω, t) (cid:3)

(cid:18)3
2 v0(t) − 1
2 uz

(cid:19)

a cos(θ)

(2.23)

on its surface. From this we conclude the total velocity potential in the fluid arises from a
surface-scattering term alongside a term describing the transfer of momentum from the
moving microsphere to the fluid.

Using eq. (2.10), we may then write the equation of motion for the system as

(cid:18)3

2

(cid:19)

ρsV (cid:219)uz (cid:3) ρ0V

(cid:219)v0 − 1
2

(cid:219)uz

.

(2.24)

(cid:16)ρs +

(cid:17)

where V (cid:3) 4πa3/3 gives the volume of the microsphere. The transfer of momentum from
the moving microsphere to the fluid becomes a reaction force of the fluid due to the sphere.
Landau & Lifshitz [53] initially derived this non-dissipative drag force by way of momentum
and energy conservation. Note that this drag force presents only in the case of accelerated
motion of the microsphere and we may recast its effect in the form of a virtual mass that
includes a contribution due to the mass of the displaced fluid,

(2.25)
This expression leads to a simple relation linking uz(t) and v0(t) provided that the velocity
does not remain constant and that the sphere does not move in the absence of the field:

2

ρ0
2

V (cid:219)uz (cid:3)

3ρ0V

(cid:219)v0.

uz
v0

(cid:3)

3ρ0

ρ0 + 2ρs

.

(2.26)

The idea of a virtual mass for the accelerated motion of a single sphere in an ideal
fluid readily generalizes to the case of a moving collection of mutually-interacting spheres.
Through this, we may compute the dynamics of each microsphere in the group, taking
into account the effect of the momentum exchange between the fluid and the microspheres.
This results in both drag and inter-particle forces in addition to the displacement caused
by the driving acoustic field.

13

u1

θ1

x1

k

u2

θ2

x2

L

d12

Figure 2.2: Perpendicular configuration.

2.4.2 Low-order interactions

We now consider two identical microspheres arranged perpendicularly to an incident
waveform as in fig. 2.2. Within the Born approximation, we may take eq. (2.19) as the
incident field and use it in place of the total field on the right-hand side of eq. (2.8), assuming
negligible contributions from scattering. In doing so, the field everywhere becomes

ϕ(r, t) (cid:3) v0(t)z +

By inserting this into eq. (2.10) for r1, we have

ρsus (cid:3) ρ0

In the limit of d12 → ∞, this becomes

us
v0

(cid:3)

4ρ0

ρ0 + 3ρs

.

14

Writing

Ω1

(cid:2)v0 − u2(cid:3) .

a3
3

cos(θ1)
|r − d12/2|2

m1u1 · ˆz (cid:3) 2πρ0a2Þ
Þ

ρ0

a5
3

cos θ2 (cid:3)

a
d12

(cid:2)v0 − u1(cid:3) +
cos2 θ1a3(cid:18)4

a3
3

cos(θ2)
|r + d12/2|2

(cid:19)

cos θ1

3 v0 − u1
3

d(cos θ1) +
v0 − u2
|r − d12|2 cos θ1 cos θ2 dφ1 d(cos θ1) .
(cid:114)(cid:16)1 − a
(cid:17)2
(cid:19)3
(cid:18)4
3 v0 − 1
3 us

ρ0(v0 − us)

(cid:16) a

d12

(cid:18) a

cos θ1

(cid:17)2

sin θ1

+

d12

(cid:19)

+

3

.

d12

and noting u1 (cid:3) u2 ≡ us due to symmetry in the initial configuration, we may expand
eq. (2.28) in a/d12 to give

(2.27)

(2.28)

(2.29)

(2.30)

(2.31)

Symbol Value
Quantity
Sound speed
cs
Microsphere radius
ak
Density (exterior)
ρ0
Density (interior)
ρs
Pulse amplitude
P0
Center frequency
f0
Pulse duration (st. dev.) σ

1500 m s−1
1 µm
1000 kg m−3
1 kg m−3
0.05 m2 s−1
0.5 MHz to 20 MHz
7 µs to 24 µs
Table 2.1: Typical simulation parameters.

By considering negligible scattered fields at the surface of each microsphere, we qualitatively
recover eq. (2.27) with different coefficients arising only from the Born approximation.

Moreover, the additional interaction term in eq. (2.31) scales as(cid:12)(cid:12)dij

(cid:12)(cid:12)−3; a behavior anticipated

from the dipolar nature of eq. (2.22).

2.5 Numerical Results

Here we present a series of numerically-solved systems to illustrate the utility of
the method in investigating acoustic phenomena. We perform simulations of one- and
two-particle/pulse systems to determine the principal particle-field and particle-particle
interactions, followed by simulations of larger assemblages of spheres to investigate group
phenomena and effects in systems without symmetry. Unless otherwise stated, table 2.1
gives the simulation parameters for each of the following simulations; as our interests
lie in hydrodynamic applications, we use material parameters characteristic of water
to define our external fluid medium. Similarly, we consider here the case of gas-filled
microspheres [13], and therefore set their density much smaller than that of the exterior
medium. The acoustic pulses lie in the ultrasonic regime, and the chosen frequency of
20 MHz corresponds to that of typical applications in acoustic microscopy.

15

)

m
µ
(

t
n
e
m
e
c
a
l
p
s
i
D

60

40

20

0
−20

0

5

10

15

25

20
30
Time (µs)

35

40

45

50

Figure 2.3: Translation of a single microsphere interacting with an incident pulse ( f0 (cid:3)
0.5 MHz, σ (cid:3) 7 µs). Microspheres interacting with the pulse translate a finite distance
along k due to the Gaussian envelope in eq. (2.9).

2.5.1 Single microspheres

Figure 2.3 gives the trajectory of a single microsphere initially at rest under the effects of an
incident Gaussian pulse. Under the linear and ideal fluid approximations and without
the Gaussian envelope in eq. (2.9), the microsphere merely oscillates about its origin in
accordance with eq. (2.26). In the pulsed case, however, the variation in pressure imposed
by the finite value of k modifies the system dynamics to yield a net translation of each
microsphere. Note that the regime considered here produces no net transfer of momentum
between the acoustic field and the microsphere—a consequence of the ideal fluid.

Figure 2.4 depicts smoothed results of 128 trajectories corresponding to single micro-
spheres initially spaced along ˆz and excited by identical counter-propagating pulses. By
taking the width of each pulse much greater than the radius of each microsphere, the
two pulses reproduce the effects of interfering standing waves. The confinement occurs
at ∇P (cid:3) 0 (nodal) planes where the net force on each microsphere vanishes. The half-
wavelength associated with the dominant pulse frequency gives the separation between

16

Pressure gradient, ˆz · ∇p (relative scale)

(cid:3) 37.5µm

cs2 f0

)

m
µ
(

i

e
t
a
n
d
r
o
o
c
-

ˆz

100

50

0

−50

−100

0

20

40

60

80

Time (µs)

100

120

140

Figure 2.4: Confinement of non-interacting spheres to planes; identical counter-propagating
pulses ( f0 (cid:3) 20 MHz, σ (cid:3) 23.8 µs) initially displaced along ˆz tend to align objects in ∇P (cid:3) 0
planes at λ/2 intervals. Here, we have sampled the field and trajectories every 30 timesteps
and smoothed the resulting data with a 16-sample windowed average.

neighboring planes.

Finally, fig. 2.5 shows the relative velocity potential near a single microsphere; given a
surface expansion of ϕ, we may compute the potential everywhere through application of
eq. (2.8). As predicted by eq. (2.22), this field greatly resembles that of a pointlike “velocity
dipole” with vs acting as a dipole moment.

The simulations described thus far demonstrate precise acoustic control; through careful
application of the incident field parameters, we may induce a (finite, given a finite pulse)
translation along the principal ˆk-vector with a large degree of accuracy in the overall
displacement. In addition, the application of multiple pulses serves to confine microspheres
to highly localized regions in space, offering a self-consistent model of acoustic tweezing.

17

)

m
µ
(

x

3

2

1

0
−1
−2
−3
−3

vs

−2

−1

0

1

2

3

z (µm)

Figure 2.5: Calculated isopotential contours near a lone microsphere. Red and blue
colorations represent regions of positive and negative potential. The motion of each
microsphere through the background medium serves primarily to produce a dipolar field
of velocity potential with vs serving as the sphere’s dipole moment.

2.5.2 Many-particle simulations

We now turn our attention to collections of mutually interacting microspheres. To quantify
the effects of scattering, we first decouple scattering forces from the incident pulse by
arranging two microspheres perpendicularly to the pulse’s k-vector. Figure 2.6 gives results
for such a simulation where we plot the relative change in velocity as compared with the
single-particle simulation,

∆|vmax| (cid:3) max(cid:0)(cid:12)(cid:12)vdouble(t) − vsingle(t)(cid:12)(cid:12)(cid:1).

(2.32)

In principle, describing quantities found from a complete simulation as a function of initial
separation could obfuscate scaling data considerably; forces arising from scattering could
alter the geometry of the system. In practice, however, the perpendicular configuration used
here gives scattering forces that only influence the motion along k. Consequently, ∆v ∝ z
and the microspheres’ initial separation remains a good estimator of scaling behavior. We
and the separation data exhibit strong |d12|−3
see in fig. 2.6 that the radii data scale as a3
k

18

10−0.410−0.310−0.210−0.1 100 100.1 100.2 100.3 100.4 100.5 100.6 100.7

10−2

Radius ( , µm)

)
e
v
i
t
a
l
e
r

,

(

|
x
a
m
v
∆

|

10−3

10−4

10−5

10−6

10−7

10−3

10−4

10−5

|

∆
v
m
a
x
|

(

,

r
e
l
a
t
i
v
e
)

101

Separation |d| ( , µm)

102

Figure 2.6: Scaling behavior of two microspheres arranged perpendicularly to an incident
pulse for various radii and initial separations. The ( , ) symbols on each axis denote data
associated with that axis. The follow a regression of ∆|v|d (cid:3) 0.250 754d
, and the
follow ∆|vmax|r (cid:3) 3.133 28 × 10−5a2.99814
. These trends strongly indicate dominant dipolar
0
interactions between microspheres.

−3.00077
12

scaling, again indicating a dominant dipolar interaction between microspheres as shown
by Ilinskii et al. in 2007 [44] and predicted by eq. (2.22).

Finally, we consider the dynamics of large (N (cid:3) 16) clouds of microspheres. For each
simulation, we generate a collection of microspheres initialized with zero velocity and
random positions within a 10 µm ball subject to a minimum-separation constraint to prevent
collisions. Figure 2.7 shows a snapshot of the velocity potential isosurfaces calculated
in one such simulation. Even with mutual interactions, the shape of each isosurface
remains consistent with the presence of a dipolar field oriented along the microspheres’
velocity. Again, due to the localization assumption used to justify eq. (2.6), each system
predominantly translates a finite distance in accordance with the results found for a single
microsphere in fig. 2.3. To quantify small changes in the geometry of a system, we compute
Vh, the volume of the convex hull containing each microsphere, at every timestep in the

19

Figure 2.7: Isosurfaces of velocity potential (arb. units) calculated by evaluating the ˆS and
ˆD terms in eq. (2.8) for a N (cid:3) 16 particle simulation. Red, blue, and yellow surfaces denote
regions of positive, negative, and zero potential respectively, with holes appearing due to
intersections with the bounding box. Rendered with VisIt[21].

simulation [47]. Figure 2.8 shows the fractional change in the hull volume,

∆Vh (cid:3)

Vh(t) − Vh(0)

Vh(0)

,

(2.33)

for 20 such systems after smoothing with a weighted moving average. Curves ending
above and below zero indicate larger and smaller hull volumes (system expansion and
contraction). We note from fig. 2.8 a greater tendency for random clouds to expand; the
effective dipole-dipole interaction between particles with dij ⊥ k gives purely repulsive
forces, while the interaction between particles with dij (cid:107) k gives both repulsive and
attractive effects depending on σ and the relative phase of the oscillating microsphere
velocities.

20

4

3

2

1

0
−1
−2

)
t
n
e
c
r
e
p
(

h
V
∆

·10−2

5

6

7

8

10

9
Time (µs)

11

12

13

14

15

Figure 2.8: Fractional change in the volume of 20 randomly-initialized microsphere clouds
subject to the same incident pulse, smoothed with a 128-sample moving average. Positive
and negative values denote expansion and contraction. σ (cid:3) 1.5 cm.

2.6 Conclusions

This work lends a novel, fine-grained approach to the study of acoustic response
via integral equation methods. By considering a potential representation in terms of
spherical harmonics on the surfaces of microspheres coupled to a standard molecular
dynamics scheme, we obtain a description of the microspheres’ dynamics under the effect
of ultrasound pulses without resorting to time-average approximations. Furthermore, the
confined microsphere geometries under consideration allow us to neglect small effects
arising from time-delays in scattering. We have shown that the net effect of an ultrasound
pulse on a single microsphere consists of a translation that we can tune through careful
control of pulse parameters. Additionally, systems with multiple incident waveforms tend
to confine microspheres to nodes in the pressure field governed by acoustic interference.
Finally, in the dynamics of systems with many microspheres, we have observed the effect
of weak inter-particle transient effects induced by the driving acoustic pulse. These effects
can produce both expansion and contraction of a cloud of microspheres, in addition to the

21

overall translation.

Prior work in this area [85, 75] makes use of deformable bubble boundaries about fixed
locations. Incorporation of these methodologies to our theoretical model naturally offers
possibilities for future research, as does the addition of retardation effects. Additionally, we
expect a straightforward approach to experimental confirmation of the results presented
here. Optical tracking of tracer particles[76] has demonstrated its effectiveness in similar
fluid-trajectory studies and would readily adapt to track physical analogues of our
theoretical microspheres.

22

CHAPTER 3

RADIATION IN COUPLED QUANTUM DOT SYSTEMS

“But I don’t want to go among mad people,” Alice remarked.
“Oh, you can’t help that,” said the Cat: “we’re all mad here. I’m mad. You’re mad.”
“How do you know I’m mad?” said Alice.
“You must be,” said the Cat, “or you wouldn’t have come here.”

—Alice and the Cheshire Cat, Alice in Wonderland

3.1 Introduction

Semiconductor structures containing a large number of quantum dots offer ideal
environments for exploring collective effects induced by light-matter interactions. Often,
these structures exhibit new phenomena due to geometrical randomness and nonlinearities
in the underlying system dynamics. Additionally, optical excitations (excitons) undergo
characteristic Rabi oscillations [70, 48, 43] in quantum dots analogous to those observed
in atomic systems; because quantum dots have stronger dipolar transitions than atoms,
these light-induced oscillations generate secondary fields that couple the system more
strongly than equivalent atomic species. We can therefore expect—at least in some regions
of the sample—these local secondary fields will produce modified collective behavior
in the exciton dynamics. Phenomena induced by these secondary fields have received
considerable theoretical/computational [69, 68] and experimental [7] attention as they may
provide new insight on the coherent dynamics of excitons in quantum dot systems.

In the realm of theoretical/computational investigation, researchers in atomic and solid-
state optics have developed numerous variations of the Maxwell-Bloch equations [38] to
describe features such as ringing in pulse propagation [18, 55] or emission fluctuations [41].
Early solution strategies for these equations fell to continuum models [63, 55] that recover

23

effects arising from far-field interactions or describe near-field phenomena assuming
spatial homogeneity [71]. More recently, mesh-based PDE solvers [79, 30, 9] added a
large degree of fidelity to these models, though the finite size of the mesh means they
still struggle to resolve short-range effects without unduly increasing the computational
cost. Additionally, the nature of these meshes makes them prohibitively expensive to
extend into higher dimensional geometries for optically-large systems. In this work we
develop a computational framework to discover signatures of collective effects in strongly-
driven quantum dots within a microscopic formalism. By constructing the Maxwell-Bloch
equations with an integral kernel to describe radiation, we recover near- and far-electric
fields with full fidelity across the simulation while allowing for dynamics at the level
of individual quantum dots. Our methodology—based on successful models of other
electromagnetic [67, 61, 60] and acoustic [28, 29, 36] systems—accommodates 104 particles
distributed over optically-large regions in three dimensions.

As we explicitly track the evolution of each quantum dot in the system, we will
numerically demonstrate that the collective Rabi oscillation can induce significant coupling
in sufficiently close quantum dots. This laser-induced inter-dot coupling manifests itself in
different forms: (i) The polarization generated in isolated quantum dot pairs dynamically
suppresses the Rabi rotation. We interpret this as the consequence of a time-dependent
energy shift that brings the pair temporarily out of resonance with the external driving
field.
(ii) In addition to this screening, we observe oscillations in the free-induction
decay for larger multiplets of quantum dots. (iii) Optical pulses of integer π area, for
which we expect no polarization in uncoupled systems following the pulse, produce
patterns of residual localized polarization that remain in the system. (iv) The long-range
interactions in optically-large systems produce wavelength-scale regions of enhanced and
suppressed polarization. These effects could, for instance, help identify multiplets of dots
that dynamically couple during Rabi oscillations, or help understand nonlinear pulse
propagation effects in these media.

24

ˆz (cid:3) |1(cid:105)

(cid:12)(cid:12)ψ(cid:11)

θ

ϕ

ˆx

ˆy

describes the evolution of(cid:12)(cid:12)ψ(cid:11) on the surface of this sphere in response to an external

Figure 3.1: Bloch sphere representation of a two-level quantum system. Equation (3.1)

−ˆz (cid:3) |0(cid:105)

electric field.

We structure the remainder of this chapter as follows: section 3.2 motivates the physical
model of an ensemble of two-level systems that interact through a classical electric field.
Section 3.3 presents the details of our methodology in the context of a global rotating-wave
approximation and we offer an implementation of this algorithm at [34]. Section 3.4
contains the results of our investigation where we observe polarization features not present
in noninteracting systems at both sub- and super-wavelength scales. Finally, section 3.5
contains concluding remarks where we hypothesize on the mechanisms underpinning the
observed polarization features as well as comment on our future work in this area.

3.2 Problem Statement

Consider the evolution of a set of quantum dots in response to a time-varying electric
If we concern ourselves only with electric dipole transitions in a resonant (or
field.
nearly-resonant) system, we may write the time-dependence of a given quantum dot’s
density matrix, ˆρ(t), as

(cid:3)

(3.1)
As we wish to investigate two-level systems, ˆρ(t) contains three unique unknowns (ρ00
and the real and imaginary parts of ρ01) that describe point on the Bloch sphere (fig. 3.1),
ˆH(t) represents a local Hamiltonian that governs the internal two-level structure of the

(cid:2) ˆH(t), ˆρ(cid:3) − ˆD(cid:2) ˆρ(cid:3) .

d ˆρ
dt

−i


25

quantum dot, as well as its interaction with an external electromagnetic field, and ˆD
provides dissipation terms that account for emission effects phenomenologically. Formally,

ˆH(t) ≡(cid:169)(cid:173)(cid:171)
ˆD(cid:2) ˆρ(cid:3) ≡(cid:169)(cid:173)(cid:171)

0
χ∗(t)

(cid:170)(cid:174)(cid:172)
(cid:0)ρ00 − 1(cid:1)/T1 ρ01/T2

χ(t)
ω0

ρ10/T2

ρ11/T1

(cid:170)(cid:174)(cid:172)

(3.2a)

(3.2b)

where χ(t) ≡ d · ˆE(r, t)/, d ≡ (cid:104)1|eˆr|0(cid:105), and |0(cid:105) & |1(cid:105) represent the highest valence and
lowest conduction states of the quantum dot under consideration. Finally, the T1 and T2
constants characterize average emission and relaxation times.

4π


Þ(I − ¯r⊗¯r) · ∂2

c2|r − r(cid:48)| + (I − 3¯r⊗¯r) ·(cid:18) ∂tP(r(cid:48), tR)

t P(r(cid:48), tR)

To account for the interactions between quantum dots, we turn to a semiclassical
description of the system under the assumption of coherent fields and negligible quantum
statistics effects. Such an approximation preserves the discrete two-level energy structure
of individual quantum dots though electromagnetic quantities behave like their classical
analogues. We define the total electric field at any point as E(r, t) (cid:3) EL(r, t) + F{P(r, t)}
where EL(r, t) describes an incident laser field, P(r, t) a polarization distribution arising
from the off-diagonal elements (coherences) of ˆρ, and

d3r(cid:48)

c|r − r(cid:48)|2 +

P(r(cid:48), tR)
|r − r(cid:48)|3

F{P(r, t)} ≡ −1
(3.3)
(see [53, section §72] or section D.2.1). Here, I denotes the identity dyad, ¯r ≡ (r − r(cid:48))/|r − r(cid:48)|,
⊗ represents the tensor product (i.e. (a⊗b)ij (cid:3) aib j), tR ≡ t − |r − r(cid:48)|/c, and 
 gives the
dielectric constant of the inter-dot medium. Thus, in a system composed of multiple
quantum dots, eq. (3.3) couples the evolution of each quantum dot by way of the off-
diagonal matrix elements appearing in eq. (3.2a). Note that this approach does not require
an instantaneous dipole-dipole Coulomb term between (charge-neutral) quantum dots; the
interactions between structures occur only via the electric field which propagates through
space with finite velocity. (See [22, sections Aiv and Civ] for in-depth discussions of this
point.)

(cid:19)

26

In the systems under consideration here, ω0 lies in the optical frequency band (∼
1500 meV/). As such, naïvely integrating eq. (3.1) to resolve the Rabi dynamics that occur
on the order of 1 ps becomes computationally infeasible. By introducing ˜ρ (cid:3) ˆU ˆρ ˆU† where
ˆU (cid:3) diag(1, e iωLt), we may instead write eq. (3.1) as

(cid:2) ˆU ˆH ˆU

† − i ˆV , ˜ρ(cid:3) − ˆD(cid:2) ˜ρ(cid:3) ,

d ˜ρ
dt

−i


ˆV ≡ ˆU

d ˆU†
dt

(cid:3)

(3.4)
which will contain only terms proportional to e i(ω0±ωL)t if E(t) ∼ ˜E(t) cos(ωLt). Conse-
quently, we ignore the high-frequency quantities under the assumption that such terms
will integrate to zero in solving eq. (3.4) over appreciable timescales [4]. Because the system
no longer contains any optical frequencies, one can then construct efficient numerical
strategies for solving eq. (3.4).

Due to the quantum mechanical transitions at play in producing secondary radiation,
we may assume similarly monochromatic radiated fields. As such, a similar transformation
applies to the source distribution in eq. (3.3). Writing P(r, t) (cid:3) ˜P(r, t)e iωLt and similarly
ignoring high-frequency terms, the radiated field envelope becomes

Þ(I − ¯r⊗¯r) ·

4π

˜P(r(cid:48), tR) + iωL

(cid:0)∂t
˜F{ ˜P(r, t)} ≡ −1
(I − 3¯r⊗¯r) ·

(cid:0)∂2
˜P(r(cid:48), tR)(cid:1)e−iωL|r−r(cid:48)|/c

˜P(r(cid:48), tR) + 2iωL∂t

˜P(r(cid:48), tR)(cid:1)e−iωL|r−r(cid:48)|/c
˜P(r(cid:48), tR) − ω2
c2|r − r(cid:48)|
+ (I − 3¯r⊗¯r) · ˜P(r(cid:48), tR)e−iωL|r−r(cid:48)|/c

L

t

+

c|r − r(cid:48)|2

d3r(cid:48) .
(3.5)
Critically, eq. (3.5) maintains the high-frequency phase relationship between sources
oscillating at ωL via the factors of e−iωL|r−r(cid:48)|/c that appear.

|r − r(cid:48)|3

3.3 Computational Approach

To solve eqs. (3.4) and (3.5) for each of Ns quantum dots at Nt equally-spaced timesteps,
we begin with a suitable representation of ˜P(r, t) in terms of spatial and temporal basis
functions, i.e.

Ns−1


Nt−1


(cid:96)(cid:3)0

m(cid:3)0

˜P(r, t) (cid:3)

˜A(m)
(cid:96) S(cid:96)(r)T(t − m ∆t).

(3.6)

27

)

t

∆
/

t

(
T

1
0.8
0.6
0.4
0.2
0
−0.2
−0.4

−1 −0.5

0.5

0
2
Timestep (dimensionless)

1.5

1

2.5

3

Figure 3.2: Nonorthogonal and C0-continuous temporal basis function T(t/∆t) constructed
from intervals of third-order Lagrange polynomials.

As the wavelength of any radiation in the system far exceeds the dimensions of the quantum
dots under consideration, we take S(cid:96)(r) ≡ d(cid:96)δ(r − r(cid:96)) where d(cid:96) and r(cid:96) denote the dipole
moment and position of dot (cid:96). Furthermore, we require the T(t) to have finite support
as well as causal and interpolatory properties so as to recover ˜P, ∂t
˜P at every
timestep. Accordingly, we have elected to use

˜P, and ∂2

t

where

λ j(t) (cid:3)

(1−τ)j
j!

(1+τ)p−j
(p−j)!

0

j − 1 (cid:54) τ < j
otherwise,

T(t) (cid:3)

λ j(t)



(3.7)

(3.8)

(a)k ≡ Γ(a + k)/Γ(a) denotes the Pochhammer rising factorial, and τ ≡ t/∆t. Such a T(t)
consists of shifted, backwards-looking Lagrange polynomials of order p (we require p (cid:62) 3
to recover a twice-differentiable function), forming a temporal basis set with functions
similar to the one shown in fig. 3.2. These functions reliably interpolate smooth functions
with controllable error and have a long history of use in studies of radiative systems [56, 14].

p


j(cid:3)0

28

Combining eqs. (3.5) and (3.6) and projecting the resulting field onto the spatiotemporal

basis functions produces a marching-on-in-time system of the form

k(cid:3)0
In this (block Ns × Ns) matrix equation

˜L(m)

+

(cid:3) ˜F (m).

˜Z(k) ˜A(m−k)

m

(cid:96) (cid:3)(cid:10)S(cid:96)(r), ˜EL(r, m ∆t)(cid:11)
(cid:96)(cid:96)(cid:48) (cid:3)(cid:10)S(cid:96)(r), ˜F{S(cid:96)(cid:48)(r)T(k ∆t)}(cid:11)
(cid:3)(cid:10)S(cid:96)(r), ˜E(r, m ∆t)(cid:11)

˜L(m)
˜Z(k)
˜F (m)

(cid:96)

(3.9)

(3.10a)

(3.10b)

(3.10c)

(3.12)

(3.13)

where (cid:104)·, ·(cid:105) denotes the scalar product of two functions. As we have elected to use a
uniform ∆t, the summation in eq. (3.9) represents a discrete convolution that will produce
˜F (m) given only ˜A(m(cid:48)(cid:54)m). To link ˜A(m) with the density matrix in eq. (3.4), we take

˜A(m)
(cid:96) ≡ ˜ρ(cid:96),01(t (cid:3) m ∆t)

(3.11)

as the off-diagonal matrix elements (coherences) of ˜ρ(cid:96) directly characterize the dipole
radiating at r(cid:96) under the rotating-wave approximation. Consequently, determining ˜A(m+1)
amounts to integrating eq. (3.4) from ti (cid:3) m ∆t to t f (cid:3) (m + 1) ∆t for every quantum dot
in the system. For this, we make use of the predictor/corrector scheme detailed in [33].
Defining tm ≡ m ∆t and approximating ˜ρ(t) as a weighted sum of complex exponentials,
the predictor/corrector scheme proceeds with an extrapolation predictor step,

and iterated corrector steps1,

˜ρ(cid:96)(tm+1) ← W−1

˜ρ(cid:96)(tm+1) ← W−1


w(cid:3)0

w(cid:3)−1

P(0)
w ˜ρ(cid:96)(tm−w) + P(1)

w ∂t ˜ρ(cid:96)(tm−w),

C(0)
w ˜ρ(cid:96)(tm−w) + C(1)

w ∂t ˜ρ(cid:96)(tm−w)

1The arrow in these operations denotes an operation akin to assignment; this predic-
tion/correction does not define x(tW), only the steps of an iteration that will converge on
x(tW) (or an approximation thereof).

29

Predict ˜ρ(tm+1)

˜ρ → ˜P

Calculate ˜E(tm+1)
with MoT system

U ntilconverged

m
→

1
+
m

Correct ˜ρ(tm+1)
with derivative

Use Liouville eqn.
to produce ∂t ˜ρ(tm+1)

Figure 3.3: Flowchart of one iteration of the coupled quantum/electromagnetic solution
procedure.

(C(0)
−1 ≡ 0 by construction) to advance the system by ∆t. Such an integrator has significantly
better convergence properties than Runge-Kutta integrators for equations of the type seen
in eq. (3.4) and naturally accommodates basis functions within c ∆t of each other. (See
appendix B for an in-depth discussion of this point.)

Figure 3.3 illustrates these computational steps. In short,

1. At timestep m, use eqs. (3.11) and (3.12) to predict

˜A(m+1). This prediction depends
only on the known history of the system and does not require the calculation of any
electromagnetic interactions.

2. Use eq. (3.9) to calculate ˜F (m+1). Having extrapolated ˜A(m+1) in step 1, ˜F (m+1) will

contain information from quantum dots within c ∆t of each other.

3. Produce ∂t ˜ρ(tm+1) (and thus ∂t

˜A(m+1)) by evaluating eq. (3.4) with the ˜F (m+1) found

in step 2.

4. Correct
4 until

˜A(m+1) with eq. (3.13) and ∂t
˜A(m+1) has sufficiently converged, then set m ← m + 1.

˜A(m+1) found in step 3. Repeat steps 2 through

30

Symbol Value
Quantity
Speed of light
c
Transition frequency
ω0
Transition dipole moment d
Decoherence times
Laser frequency
Laser wavevector
Pulse width
Pulse area

300 µm ps−1
1500 meV/
10 e a0 (uniform)
10 ps and 20 ps
1500 meV/
ωL/c (kL · d (cid:3) 0)
1 ps
π

T1, T2
ωL
kL
σ/ωL
–

Table 3.1: Simulation parameters (unless otherwise stated); e and a0 denote the elementary
charge and Bohr radius. The decoherence times here, while shorter than those typical of
optical resonance experiments, afford a shorter computational time but preserve dynamical
emission phenomena.

3.4 Numerical Results

Here we detail the results of investigations into coupled quantum dot behavior with the
model presented thus far. Our algorithm reliably handles tens of thousands of quantum
dots and can simulate ten picoseconds of system dynamics in two days on a single processor.
We perform simulations of systems of quantum dots randomly distributed throughout a
simulation volume experiencing a laser pulse of the form

˜EL(r, t) (cid:3) ˜E0e

−(kL·r−ωLt)2/(2σ2).

(3.14)

Table 3.1 provides the physical system parameters unless otherwise stated.

3.4.1 Stability & adjacency effects
Figure 3.5 details ρ00(t) for two-, and 1024-particle simulations shown against a solution of
eq. (3.4) for a quantum dot evolving according to ˜EL(r, t) alone. The system in fig. 3.5(a)
contains two quantum dots with a separation of 6.3 nm perpendicular to their mutual
d to ensure large contributions from the near-field term of eq. (3.5). The two quantum
dots in this system follow the same trajectory and both excite far less than either would
in response to the incident laser alone. These signatures indicate the effect arises as a
dynamical frequency shift that brings adjacent quantum dots out of resonance with the

31

)
0
0
ρ
(
n
o
i
t
a
l
u
p
o
P

1

0.8

0.6

0.4

0.2

0

Fixed frame
Rotating frame

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Time (ps)

Figure 3.4: Comparison of ρ00(t) using both fixed (requiring 59959 timesteps) and rotating
(requiring 160 timesteps) reference frames for two interacting quantum dots. Both frames
produce similar trajectories, however the inset spy reveals the fixed frame contains a minute
oscilliatory term that the rotating frame does not.

applied electric field. We note that this suppression effect occurred to varying degrees
for all near-field arrangements of two quantum dots that we investigated. In fig. 3.5(b),
the system contains the same two-dot arrangement as in (a), however we have added an
additional 1024 quantum dots randomly distributed throughout a 6.4 × 10−2 µm3 cube
centered around the original pair. Most of these additional quantum dot populations
deviate little from those produced by the laser-only pulse due to the large separation
between particles. Nevertheless, as we have filled the cube randomly, some regions of
the system contain localized clusters of quantum dots that produce the suppression effect
detailed above (for two adjacent particles) or populations with higher frequency oscillations
(in the case of clusters with three or more quantum dots). Specifically, the two “original”
quantum dots acquire an out-of-phase oscillation with respect to each other as well as a
lower frequency in-phase oscillation of the pair about a decaying envelope.

0.7

Figure 3.6 further illustrates the near-field coupling between quantum dots. Here, 1024
quantum dots randomly fill an 8 × 10−3 µm3 cube and the same 1 ps Gaussian π-pulse

32

1

0.8

0.6

0.4

0.2

0

1

0.8

0.6

0.4

0.2

0

n
o
i
t
a
l
u
p
o
P

)
y
l
n
o
r
i
a
p
(

n
o
i
t
a
l
u
p
o
P

)
d
o
o
h
r
o
b
h
g
i
e
n
(

(a)

ρ(1)
00(t)

ρ(2)
00(t)

(b)

0

ρ(1)
00 & ρ(2)
00
˜EL-only

5

10

15

20

Time (ps)

Neighborhood
30

25

Figure 3.5: Population dynamics of adjacent quantum dots (top, shown with the trajectory
of a quantum dot driven exclusively by the external laser in black) and a 1024-particle neigh-
borhood (dots (1) and (2) maintain their separation and orientation between simulations).
The interaction between particles gives a greatly diminished response to the external pulse
through a dynamical detuning of the two-dot system. The majority of the neighboring
particles in the 1024-dot system follow trajectories nearly identical to that prescribed by
the laser (omitted for clarity); many-dot effects, however, produce significant oscillatory
00 and ρ(2)
modes in the evolution of selected quantum dots (shown in gray). Note that ρ(1)
00
appear to have two coherent modes in the presence of multiple dots: a high-frequency
oscillation between the pair, as well as a low-frequency oscillation of the group about a
decaying envelope.

33

( µ m )
0

z
0.1

0.2

−0.2

−0.1

−0.2

0
x(µm)

0.5

0.4

0.3

0.2

0.1

0.2

0.1

0
−0.1
−0.2

)

m
µ
(

y

0

0.2

Figure 3.6: Spatial distribution of(cid:12)(cid:12) ˜ρ01(cid:12)(cid:12) for 1024 quantum dots as an indicator of polarization.

Recorded at t (cid:3) 2 ps relative to the peak of a 1 ps-wide π-pulse, the color and size of each
sphere indicates the location of each quantum dot and its polarization. Following a π-pulse,
a single quantum dot would have no remnant polarization; here, due to the near-field
interactions between particles, clusters of quantum dots remain in highly-polarized states
depending on their separation.

×10−2

2

1.5

R
P
I

1

0.5

0

0

5

10

15

20

25

30

35

40

Time (ps)

Figure 3.7: Inverse participation ratio (IPR) for the system shown in fig. 3.6. The dashed
line indicates the sample time for fig. 3.6.

34

illuminates the system. Without any inter-dot interactions, such a pulse would perfectly
transition every quantum dot from |0(cid:105) to |1(cid:105), leaving behind no polarization after the pulse
has passed. As the system in fig. 3.6 has experienced most of the pulse, the majority of
the quantum dots behave this way. A number of particles remain polarized well after the
pulse has passed through the system, however, due to their apparent proximity.

3.4.2 Polarization enhancement

Borrowing from standard measures of localization (in which one calculates integrals of E
over the simulation volume), we have adapted the inverse participation ratio (IPR) of the
dot polarization as


(cid:96) | ˜p(cid:96)(t)|4
(cid:16)
(cid:96) | ˜p(cid:96)(t)|2(cid:17)2

IPR(t) (cid:3)

(3.15)

to provide a quantitative description of these phenomena [66]. Figure 3.7 shows this
quantity for the system in fig. 3.6. The maximum IPR occurs some time after the pulse
has passed through the system, indicating strongly-coupled quantum dots retain their
polarization longer than their neighbors. Moreover, this dynamical localization effect
features oscillations which suggests many-dot effects contribute to the dynamics within a
narrow spectral region.

Figure 3.8 depicts the evolution of(cid:12)(cid:12)ρ01(r)(cid:12)(cid:12) as an indicator of ˜P for a cylinder containing

10 000 quantum dots. The cylinder has a radius of 0.2 µm and a length of 4 µm, and the
incident kL lies along the cylindrical axis (again perpendicular to d so as to maximize the
long-distance interaction between quantum dots). This simulation captures the suppression
effects of figs. 3.5 and 3.6 as a small number of quantum dots remain in an unexcited state
while the cylinder polarizes around them. Additionally, due to the length of the cylinder,
larger regions of enhanced polarization begin to appear as the system polarizes—an effect
we did not observe in sub-wavelength structures. We liken these nodes to standing waves
in a cavity that arise from the far-field interaction term of eq. (3.5). As the pulse varies little

35

(cid:54) 0.32 0.33

0.34

0.35

0.36

0.37

0.38

0.39 0.40 (cid:54)

e
m
T

i

Figure 3.8: Coloration of(cid:12)(cid:12) ˜ρ01(cid:12)(cid:12) as an indicator of(cid:12)(cid:12)˜P(cid:12)(cid:12) at t1 (cid:3) −0.05 ps, t2 (cid:3) 0 ps, t3 (cid:3) 0.05 ps,

and t4 (cid:3) 0.10 ps relative to the peak of a 1 ps-wide pulse. 10 000 quantum dots randomly
distributed throughout a 0.2 µm (radius) × 4 µm cylinder oriented along kL demonstrate
near-field the effects of fig. 3.5 as distinct, outlying bright/dark quantum dots. Additionally,
the size of the system allows for wavelength-scale phenomena that appear here as five
standing regions of enhanced polarization. (Note that we model quantum dots as point
objects; the size of the spheres here has no physical interpretation.)

ˆk

over the length of the cylinder, identical simulations run without interactions (i.e. Z (cid:3) 0
everywhere) produce homogeneous polarization distributions. Reducing the simulation to
a planar geometry (fig. 3.9) preserves both the short-range (dark, adjacent quantum dots)
and long-range (regions of enhanced/diminished polarization) phenomena observed in
fig. 3.8. We note these geometries produce a much weaker effect than their fully three
dimensional counterparts, though larger simulations with greater far-field contributions or
quantum dots engineered to have a larger dipole moment may produce more measurable
effects.

36

(cid:54) 0.461

0.463

0.465

0.467

0.469 (cid:54)

ˆk

1
0.8
0.6
0.4
0.2
0

Figure 3.9: Coloration of(cid:12)(cid:12) ˜ρ01(cid:12)(cid:12) for 10 000 quantum dots arranged in a finite planar geometry

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3

(in units of λ). The slab displays a prominent polarization pattern 1.25 ps after the peak of
a 1 ps pulse.

3.4.3

Inhomogeneous broadening

Thus far we have investigated only homogeneous systems (i.e. an identical ω0 for
every quantum dot in the system). To probe the effects of interaction-independent
inhomogeneities (possibly occurring due to some experimental variation in quantum dot
sizes), fig. 3.10 presents four simulations with normally-distributed transition frequencies
characterized by a width parameter.
In simulations with mild detuning, we observe
a distribution pattern characteristic of the enhancement phenomenon seen in fig. 3.8.
More variation in the detuning distribution, however, quickly serves to destroy these
phenomena, leaving only “pulse-driven” effects. At present, inhomogeneous broadening
in quantum dots typically exceeds the values used in fig. 3.10. Impurity-bound excitons
have greatly diminished inhomogeneous effects and thus could offer an avenue to observe
these phenomena [74].

3.5 Conclusions & future work

Here we developed a robust, fine-grained algorithm to solve for the dynamics of an
ensemble of quantum dots that couple in response to external light fields. By making use of

37

0 meV

0.2 meV

0.4 meV
1
2

−2

−1

0.6 meV
1
2

)
s
s
e
l
n
o
i
s
n
e
m
d
(

i

(cid:12)(cid:12)

(cid:12)(cid:12)

1
0
ρ
)
s
s
e
l

n
o
i
s
n
e
m
d
(

i

(cid:12)(cid:12)

1
0
ρ

(cid:12)(cid:12)

0.36

0.34

0.32

0.36

0.34

0.32

−2

−1

0

z (µm)

z-distribution of polarization(cid:12)(cid:12) ˜ρ01(cid:12)(cid:12) for the geometry in fig. 3.8.

Figure 3.10:
In each
simulation, 10 000 quantum dots had random Gaussian noise (with width parameters
of 0.0 meV, 0.2 meV, 0.4 meV and 0.6 meV) added to a resonant ω0 (cid:3) 1500 meV. The
induced long-range patterns in the polarization remain for mild detunings (cid:46) 0.5 meV
but have completely washed out at 0.6 meV. Additionally, each simulation displayed the
characteristic near-field coupling effects of fig. 3.5 (not visible here).

0

z (µm)

an integral equation kernel to propagate radiated fields, our model facilitates simulations
of thousands of quantum dots in three dimensions with accurate bookkeeping of both near
and far radiation fields. Our simulations predict a suppression effect between adjacent
quantum dots that screens out the incident laser pulse and we interpret this effect as a
dynamical detuning that shifts the effective ω0 of the affected quantum dots. Moreover,
we observe additional oscillatory behavior and localization effects in larger clusters of
particles. These effects could prove useful in optically identifying quantum dot “molecules”
in an extended sample by detecting residual localized polarization following integral π
pulse(s)—we expect that an experimental π-pulse calibrated to a single quantum dot with
a scanning-type polarization measurement [7] would reveal signatures of these effects
in dense samples. Finally, in larger systems of densely-packed quantum dots, we see
significant localization/enhanced polarization over length scales comparable to that of the

38

incident wavelength. These effects persist in simulations with other extended geometries
and in simulations with inhomogeneously-broadened transition frequencies, though the
effects quickly disappear with only a few meV detuning.

Semi-classical approaches can describe some superradiant effects within a continuum
formulation [38, 15, 16]. First predicted in 2005 [72] and observed in 2007 [65], superradiant
effects in quantum dot ensembles have since spurred theoretical analyses into cooperative
radiation mechanisms [73, 20]. While our semiclassical approach accounts for collective
effects due to the secondary field emission from quantum dots, we do so in the Hamiltonian

term on the right hand side of eq. (3.1) and not in the ˆD(cid:2) ˆρ(cid:3) dissipator. In future work, we

Unfortunately, the naïve O(cid:0)N2

s

plan to extend our microscopic approach to include collective dissipation effects so as to
better model superradiant phenomena. We expect that our approach—when extended to
systems containing a larger number of quantum dots—will aid in investigating the role of
many-dot interactions in systems such as nanolasers [46] that exploit these phenomena.

(cid:1) interaction calculation presented here hampers attempts

to extend these calculations to systems more than ∼ 104 spatial unknowns. Our ongoing
research includes the development of accelerated computational techniques that exploit
the structure of Z to reduce the big-O complexity of eq. (3.9). Additionally, the technique
presented here readily extends to model atomic, molecular, and semiconductor systems
with richer structure (e.g. systems with energy degeneracies or biexcitonic transitions).

39

CHAPTER 4

AN FFT-ACCELERATED SOLUTION METHODOLOGY FOR RADIATIVE

MAXWELL-BLOCH SYSTEMS

“Curiouser and curiouser!”

4.1 Abstract

—The Cheshire Cat, Alice in Wonderland

Th e Maxwell-Bloch algorithm of chapter 3 scales poorly for large numbers of quantum
dots due to an implicit O(n2) spatial complexity in the propagation calculation (eqs. (3.9)
and (3.10)). As many physical systems have the same mathematical structure, numerous
techniques exist to reduce the poor scaling of these convolutions with fully controllable
error. This chapter will explain this bottleneck, several acceleration schemes that one
might use to alleviate it, and, why we feel one particular variety of accelerator works
best. Additionally, I will detail the mathematics of such an acceleration technique, offer
a reference implementation online, and give physical results found through simulating
extended systems.

4.2 Introduction

The naïve O(n2) complexity of matrix-vector products in integral equation systems
stands as the singular bottleneck to large-scale simulations of related phenomena. Since
the 1990s, so-called “fast” solvers have played a critical role in alleviating this bottleneck;
by exploiting properties of the underlying physical system (such as smoothness and
translational invariance) these methods achieve space and time complexities of O(n log n)
or even O(n) with fully-controllable error1.

1In particular this means one can choose parameters such that the error in the calculation
falls below the floating point error of a given machine, thus fast solution techniques produce
identical results to their naïve counterparts with no discernible approximation error.

40

Kernel-type solvers [37, 27] exploit the physical similarity of sources at large distances to
develop a series representation of the Green’s function in terms of source and observation
regions independently. Specifically, these techniques lump adjacent sources together so
as to transmit their effects in aggregate across large distances (a process not dissimilar to
how a telephone network or the Internet works). These methods scale very well for surface
discretizations as they intelligently partition the simulation space into a hierarchy of boxes
and only compute interactions between regions containing sources. Unfortunately, as these
methods rely on the underlying mathematical structure of the kernel to produce a suitable
series representation, they have limited utility in systems with nonstandard (specifically
rotating-frame) kernels.

In contrast, source-type solvers [12, 84, 49] apply effective transformations to the source
distribution (leaving the actual interaction kernel unchanged) with the intention of
producing a nearly-equivalent system with exploitable low-rank properties. The Adaptive
Integral Method (AIM), which we extend here, represents a kernel matrix as a sum of
nearfield and farfield contributions, each constructed from a different set of basis functions.
For all basis functions within a prescribed distance of each other, we compute their
interactions directly, thus the nearfield matrix becomes increasingly sparse in the limit
of a large number of sources. The farfield contribution arises as the interaction between
a reduced number of auxiliary basis functions (commonly taken as δ-functions) on a
regular Cartesian grid. Each source function maps to a finite number of overlapping
auxiliary functions—thus the mapping matrices have a sparse signature and compress
the distant portions of the interaction matrix—and the translationally-invariant structure
imposed by the Cartesian grid allows for an efficient matrix diagonalization by way of a
(multidimensional) fast Fourier transform (FFT). Because AIM does not rely on the specific
mathematical form of the interaction kernel (a so-called “algebraic acceleration”), it works
equally well for rotating and fixed-frame problems.

41

4.3 Problem statement

(cid:40)

(cid:41)

We wish to calculate the field radiated from a distribution of time-dependent sources
under the action of an arbitrary (linear) differential operator, F{·}. These fields exist
throughout space and time alongside/in response to an incident field, Einc. Formally, we
may write the total field everywhere as

E(r, t) (cid:3) Einc(r, t) + F

(cid:120) g(r − r(cid:48), t − t

≡ Einc(r, t) + F(cid:8)g(r, t) ∗ P(r, t)(cid:9)

(cid:48))P(r(cid:48), t

(cid:48)) d1t

(cid:48) d3r(cid:48)

(4.1)

where g(r, t) denotes a propagation potential (most commonly the retarded potential,
g(r, t) (cid:3) δ(t − |r|/c)/|r|). To numerically evaluate eq. (4.1), we discretize P(r, t) with
separable space/time basis functions such that

P(r, t) ≈ Ns−1


Nt−1


(cid:96)(cid:3)0

m(cid:3)0

A(m)
(cid:96) s(cid:96)(r)T(t − m ∆t)

(4.2)

where ∆t denotes a fixed time interval chosen to accurately sample the dynamics of the
physical quantities involved. Furthermore, we require both s(cid:96)(r) and T(t) to have finite
support and that T(t) obey (discrete) causality (i.e. T(t) (cid:3) 0 if t (cid:60) [−∆t, Tmax]). Here, we
consider δ-functions and shifted Lagrange polynomials for the s(cid:96)(r) and T(t), respectively,
though this analysis readily extends to accommodate any similar set of functions.

Substituting eq. (4.2) into eq. (4.1) and projecting the resulting field onto the same set of

s(cid:96)(r) produces an (Ns × Ns)-block matrix equation of the form
F (m−m(cid:48)) · A(m(cid:48))

E(m)

m


(cid:3) E(m)
inc +

m(cid:48)(cid:3)0

where

E(m)
(cid:96) (cid:3) (cid:104)s(cid:96)(r), E(r, m ∆t)(cid:105)
E(m)
inc,(cid:96) (cid:3) (cid:104)s(cid:96)(r), Einc(r, m ∆t)(cid:105)
F (k)

(cid:96)(cid:96)(cid:48) (cid:3)(cid:10)s(cid:96)(r), F(cid:8)g(r, t) ∗ s(cid:96)(cid:48)(r)T(k ∆t)(cid:9)(cid:11).

42

(4.3)

(4.4a)

(4.4b)

(4.4c)

Figure 4.1: the interaction matrix for g(r, t) (cid:3) |r|−1 between points randomly distributed
throughout a unit cube appears to have very little structure (left). while we cannot
fully order the points in three dimensions, we can permute the matrix so as to give
points within the same neighborhood adjacent indices. points ordered thusly produce a
diagonally-dominant interaction matrix with a hierarchical Toeplitz substructure (right).
such structures indicate portions of the matrix have very accurate low-rank approximations
that offer the possibility of significant compression.

From this, it immediately becomes apparent that eq. (4.4c) bottlenecks the field calculation
due to the O(N2

s) complexity of the inner product.

4.4 Computational approach

To effect a sub-quadratic calculation of eq. (4.3), we approximate F (k) as a sum of
near- and far-field contributions. The near-field matrix elements follow directly from
eq. (4.4c)—sources within a prescribed distance threshold interact “directly” so as to avoid
incurring unreasonable approximation error between adjacent basis functions. Sources
beyond this threshold, however, interact via auxiliary spatial basis functions taken to reside

at the vertices of a regular Cartesian grid. These auxiliary sources recover F(cid:8)g(r, t) ∗ P(r, t)(cid:9)

at large distances and have two computational advantages: (i) they compress the interaction
matrix by representing sources within the same spatial region in terms of the same auxiliary
set (fig. 4.2) and (ii) they impose a Toeplitz structure on the resulting interaction matrix
that lends itself to efficient diagonalization through application of an FFT. One iteration of

43

Expansion region

Sources

Box

rbox

Interior grid

sy

sx

Figure 4.2: Illustration of the grid structure and related terminology. All of the sources
within a box (shown as the central shaded square) map to the same set of expansion points
(shown as open circles) indexed relative to rbox.

our algorithm, then, proceeds as follows:

1. At timestep m project each of the s(cid:96)(r)A(m)

(cid:96) onto the auxiliary sources. Aside from
discretization/sampling criteria, the operators in eq. (4.1) do not affect these projec-
tions, thus the distribution of auxiliary sources, Paux(r, t), mimics the distribution of
P(r, t) at large distances.

2. Effect the convolution in eq. (4.1) between auxiliary sources. Having imposed a
regular structure on Paux(r, t), we may efficiently diagonalize the matrix representing
this (discrete) convolution with (up to four-dimensional) blocked FFTs. Note that the
algorithm thus far has essentially evaluated the potential, g(r, t) ∗ P(r, t), at t (cid:3) m ∆t
at every point u in the grid.

3. Recover the total field under the action of F by projecting the potential on each u

44

back onto the s(cid:96)(r). These projections make use of specialized projection matrices
that depend on the particular derivatives contained inside F.

4. Subtract the fields determined by steps 1-3 for pairs of spatial basis functions within
a prescribed distance threshold and replace it with eq. (4.4c). The auxiliary grid
approximations only remain accurate at large distances, thus this step corrects large
approximation errors that occur between adjacent s(cid:96)(r).

Mathematically,

where

(cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171)

∂0
t
∂1
t
∂2
t

G(m−m(cid:48))
G(m−m(cid:48))
G(m−m(cid:48))

...

(cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)

†

Λ

R(cid:96)(cid:96)(cid:48) (cid:54) γ
otherwise,

(cid:48)) ∆t(cid:1) , g(r, t) ∗ ub(r)T(t)(cid:11)

F (m−m(cid:48)) ≈ F (m−m(cid:48))

direct + ΛF

≡ F (m−m(cid:48))

direct + F (m−m(cid:48))

FFT

− F (m−m(cid:48),τ)

(cid:96)(cid:96)(cid:48)

F (m−m(cid:48))
(cid:3)(cid:10)ua(r)δ(cid:0)t − (m − m

FFT,(cid:96)(cid:96)

F (m−m(cid:48))
direct,(cid:96)(cid:96)(cid:48) (cid:3)

0

G(m−m(cid:48))

ab

(4.5)

(4.6a)

(4.6b)

The Λ matrices in eq. (4.5) denote the (sparse) projections to and from the grid (detailed in
section 4.4.1), and ua(r) indicates an auxiliary basis function on the spatial grid indexed
by a. Finally, τmax and γ serve as adjustable input parameters to control the accuracy of
the simulation and R(cid:96)(cid:96)(cid:48) gives the minimum distance (in integral units of the grid spacing)
between the expansion regions enclosing s(cid:96)(r) and s(cid:96)(cid:48)(r) (fig. 4.3) via

grid

(cid:96)(cid:96)(cid:48) (cid:3) min{(cid:107)u − u(cid:48)(cid:107)∞ | u ∈ C(cid:96), u(cid:48) ∈ C(cid:96)(cid:48)}.

R

(4.7)

45

∆s

rb

∆s

∆s

r0

∆s

∆s + 1

ra

∆s

rc

Figure 4.3: Illustration of the nearfield criterion for a third order expansion (γ (cid:3) 2). The
dashed line indicates the complete nearfield of the box associated with r0—i.e. all boxes
that have an expansion point within ∆s of the expansion around r0. Consequently, all of
the s(cid:96)(r) within the central dark blue square have a pairwise interaction with the s(cid:96)(cid:48)(r)
inside the dashed box.

46

Expansion A

Expansion B

Expansion C

Figure 4.4: Illustration of nearfield corrections between close boxes. Expansions A and B
overlap, but only box B lies in the nearfield of box C. The grid-based propagation strategy
only remains accurate for distant source/observer pairs. To avoid incurring undue error,
we remove the interaction “through the grid” between the BC pair (red line) and replace it
with a more accurate “direct” interaction (dashed blue line). The AC pair requires no such
treatment as they have well-separated expansion regions.

4.4.1 Expansion matrices
We represent the primary s(cid:96)(r) basis functions as a weighted sum of δ-functions on the
surrounding gridpoints, thus ua(r) ∝ δ(r − ra) and

ψ(cid:96)(r) ≈


†
Λ
(cid:96)u

δ(r − u).

u∈C(cid:96)

Here, ψ(cid:96)(r) ∈ (cid:8)s(cid:96)(r) · ˆx, s(cid:96)(r) · ˆy, s(cid:96)(r) · ˆz(cid:9) and C(cid:96) denotes the collection of grid points

within the expansion region of s(cid:96)(r) (fig. 4.2). For an expansion of order2 m, this sum
contains (m + 1)3 terms corresponding to the (m + 1)3 grid points nearest to s(cid:96)(r) (fig. 4.5).
†
matrices contain few nonzero elements and we have elected to use a
Consequently, the Λ
(cid:96)u
moment-matching scheme to capture the (m + 1)3 multipole moments of s(cid:96)(r) according to

(4.8)

Þ(x − x0)mx(y − y0)my(z − z0)mz

(cid:34)

ψ(cid:96)(r) −


†
Λ
(cid:96)u

δ(r − u)

u∈c(cid:96)

d3r (cid:3) 0.

(4.9)
In this expression, 0 (cid:54) mx, my, mz (cid:54) m and r0 ≡ x0ˆx + y0 ˆy + z0ˆz denotes the origin about
†
, we then solve the least-squares
which we compute the multipoles. To determine the Λ
(cid:96)u
2In principle, one could expand to different orders in different Cartesian directions,
though this involves considerably more bookkeeping for relatively little benefit. Thus,
for convenience, we take “the expansion order” to mean the expansion order in every
direction.

(cid:35)

47

Figure 4.5: Spatial expansion pattern for s(cid:96)(r) (cid:3) ˆd(cid:96)δ(r−r(cid:96)) through order five. The numbers
on the grid indicate the points added in the equivalent order expansion such that an
expansion of order m requires (m + 1)3 grid points. Moreover, increasing the expansion
order incorporates new points in such a way as to keep s(cid:96)(r) as close to the center of the
expansion region as possible.

5

4

4

4

4

44

5

3

2

2

22

4

5

3

1

0

2

4

5

3

11

1

2

4

r(cid:96)

5

33

3

3

3

4

55

5

5

5

5

5




u∈C(cid:96)

†
WmuΛ
(cid:96)u (cid:3) Q(cid:96)m

Þ

system

where

(4.10)

(4.11a)

Wmu (cid:3) (ux − x0)mx(uy − y0)my(uz − z0)mz
Q(cid:96)m (cid:3)

ψ(cid:96)(r)(x − x0)mx(y − y0)my(z − z0)mz d3r ,

u ∈ C(cid:96), and m denotes the multi-index m (cid:3)(cid:8)mx, my, mz

(cid:9). With an infinite precision

(4.11b)

calculation, the choice of r0 (cid:3) x0ˆx + y0 ˆy + z0ˆz merely defines an origin for the polynomial
expansion system. To minimize numerical issues, we choose r0 at the center of s(cid:96)(r).

As Λ

† arises purely as a function of space with no time component, differentiating
eq. (4.9) with respect to x, y, or z can also recover spatial derivatives occurring in F. We
interpret this as as differentiating the underlying polynomial that interpolates the potential
between gridpoints at the location of s(cid:96)(r), thus the differentiation procedure removes the
high-order moments in eq. (4.11b).

48

4.4.2 Complexity

We have largely designed this acceleration algorithm to accommodate time-domain
propagators within a rotating-wave/envelope-tracking framework. Such systems notably
have decoupled length and time scales; downshifting a narrow-band signal to baseband
(in time) affords very efficient temporal discretizations due to an assumed carrier signal.
No such trick works for the spatial properties, however, due to the space/time coupling in
g(r, t) and need to preserve interference phenomena throughout space3. As such, g(r, t)
can impose a high spatial frequency on fields radiating from artificially low-frequency
sources. Consequently, the both the primary and auxiliary spatial basis functions require
a discretization on the order of this spatial frequency to accurately sample fields. The
calculation of eq. (4.5), then, remains the dominant factor in the complexity of our algorithm,
though it performs significantly better than a pairwise interaction between all of the s(cid:96)(r).
†) necessarily
Projecting sources to/from the auxiliary grid (i.e. right-multiplying by ΛF or Λ
require O(Ns) operations per timestep as each s(cid:96)(r) maps to a fixed number ua(r) auxiliary
functions, independent of Ns. Similarly, for uniformly-dense volume-filling geometries, a
volume proportional to γ3 contains a fixed number of s(cid:96)(r) independent of Ns, ergo the
calculation of F (m−m(cid:48))
direct A(m(cid:48)) also requires O(Ns) operations. Thus, the potential evaluation
between auxiliary functions—multiplying by G(m−m(cid:48)) and its derivatives—dominates the
overall complexity of the calculation. By storing a diagonal, k-space representation of
these matrices (requiring storage proportional only to the number of grid points, Ngrid, not
its square), each timestep in the algorithm requires an O(Ngrid log Ngrid) forward FFT to
equivalently transform the auxiliary source distribution, an O(Ngrid) elementwise product
to effect the “convolution”, and a backward O(Ngrid log Ngrid) to restore the potential in r-
(not k-) space. Similar Toeplitz diagonalization strategies also work in time (see [84, section
D]) though this provides almost no benefit in rotating-frame systems due to the relatively

3Because of this, rotating-wave systems share a large number of similarities with their

fully frequency-domain counterparts.

49

small number of timesteps required for radiation to traverse the system.

Additionally, we note that transferring the spatial derivatives in F onto the “receiving”
projection matrices, ΛF, results in a modest acceleration for tensorial propagators such as
that of the electric, magnetic, or combined field integral equations,

− c2∇∇ · .

F (cid:3) ∂2
t

(4.12)
Instead of projecting both s(cid:96)(r) and ∇ · s(cid:96)(r) onto the auxiliary grid and propagating
them with the relatively expensive G multiplication, we reconstruct all spatially-varying
quantities from the potential sampled on the auxiliary grid. (A similar method reconstructs
time derivatives from the history of potential samples on the auxiliary grid, though
this requires a discrete convolution with several weighted ΛF matrices.) Moreover, this
decomposition maintains the strong form of the integral relation in eq. (4.1) and makes no
assumptions of differentiability in the underlying s(cid:96)(r) discretization.

4.5 Results

4.5.1 Convergence

Consider two point particles located at xsrc and xobs. A time-independent Green’s function,
g(xobs − xsrc), describes the interaction between the two particles and we wish to construct
a polynomial approximation of g(x − xsrc) for x in the vicinity of xobs as in fig. 4.6.

To construct an interpolation polynomial over the expansion region of order M,
we define a polynomial coordinate xp in units of h such that xmin
p + M
≡ −(cid:98)M/2(cid:99). Consequently, the expansion points about xobs correspond to
where xmin
xp ∈ {−(cid:98)M/2(cid:99), −(cid:98)M/2(cid:99) + 1, −(cid:98)M/2(cid:99) + 2, . . .} with the 0th order expansion point, x0,
equivalent to xp (cid:3) 0. Such a coordinate system defines the Vandermonde linear equation

p (cid:54) xp (cid:54) xmin

p

50

xsrc

∆s

m (cid:3) 4

m (cid:3) 2

m (cid:3) 0

m (cid:3) 1

m (cid:3) 3

xobs

x4

x0

x2
x1
expansion region

x3

Figure 4.6: Polynomial interpolation of g(x − xsrc) near xobs. Here, the green curve
represents the actual g(x − xsrc) and the dashed black line its approximation. Evaluating
the mth-order approximation requires samples of the signal at m +1 grid points surrounding

xobs.


j Vijw j (cid:3) gi for the weights of an interpolating polynomial4 [62] where

(4.13a)

(4.13b)

p + i)j

(cid:16)(x0 − xsrc) + (xmin

Vij (cid:3) (xmin
gi (cid:3) g

p + i) ∆s

(cid:17)

g(xobs − xsrc) ≈ M


wi

i(cid:3)0

(cid:16) xobs − x0

(cid:17) i

∆s

and 0 (cid:54) i, j (cid:54) M. Approximating g(x − xsrc) at xobs then becomes a matter of evaluating
this polynomial at xp (cid:3) (xobs − x0)/h, i.e.

.

(4.14)

Accordingly, the polynomial approximation to g(xobs− xsrc) contains terms of order O(h−M)
and we can expect the approximation error to scale as O(∆s−(M+1)) as demonstrated in
fig. 4.7. This also motivates using the approximation to calculate interactions involving
differential operators; applying an nth-order derivative reduces the polynomial order by n,
thus the error scales like O(∆s−(M+1)+n).

4In principle this analysis works equally well in the original(xsrc, xobs) coordinate system.
The polynomial coordinate has three advantages, however: it indexes the expansion points
∈ Z and thus the Vandermonde matrix has infinite
“logically” from left to right, xmin
precision, and it makes the interpolation error in terms of h explicit.

p

51

r
o
r
r
e

2
(cid:96)

101
100
10−1
10−2
10−3
10−4
10−5
10−6
10−7

62.27(∆s)0.96
64.57(∆s)1.87
202.8(∆s)3.13
231.4(∆s)3.91
538.3(∆s)5.22
0.20λ

0.02λ

0.05λ
Grid spacing, s

0.10λ

(cid:96)2 error calculation of g(r − r(cid:48)) (cid:3) e ik|r−r(cid:48)|/|r − r(cid:48)| for expansion orders zero
Figure 4.7:
through four. For each of the points above, a source and observation box separated by
∆r (cid:3) 10λˆx + 10λ ˆy + 10λˆz each contain 64 randomly placed points, thus the points within
each box all map to identical nodes in the auxiliary grid. The moment-matching expansion
scheme dictates that the overall error for an expansion of order m scales as O(∆sm+1).

4.5.2 Timing

Having demonstrated the convergence properties of the accelerated algorithm, we now
examine its computational efficiency. As a representative simulation, we discretize using
δ-functions in space and tesselate the same arrangement of basis functions throughout
each box of a cubic grid. This ensures a predictable, volume-filling geometry that captures
the major features of the quantum dot systems we wish to investigate. We use a first-order
expansion to project to/from the grid and a fifth-order Lagrange basis in time, and the
simulation runs for 1024 timesteps, and take γ (cid:3) 1 (implying expansions separated by at
most one box require a nearfield correction). Finally, each source has a prescribed Gaussian
amplitude of width 1024 ∆t/12 and centered about 1024 ∆t/2. These parmaeters give a
relative error of ∼ 1 × 10−6 in the full fft + nearfield calculation when compared with an
analytic calculation. The FFT-accelerated simulation outpaces the direct field calculation
near Ns (cid:3) 400.

52

)
s
(

e
m

i
t
n
u
R

103

102

101

100

10−1

10−2

fft + nearfield

direct
fft only

102

103

Number of spatial basis functions

Figure 4.8: Time to evaluate retarded potential vs. number of spatial basis functions for
several system sizes. Simulated on an Intel(R) Core(TM) i5-4690K CPU @ 3.50 GHz.

Regression
Simulation
t(n) (cid:3) 2.100 86 × 10−5n2.073
direct
fft + nearfield t(n) (cid:3) 1.363 62 × 10−3n1.385
t(n) (cid:3) 2.108 96 × 10−4n1.169
fft only

r2
0.999946
0.999590
0.999577

Table 4.1: Power-law regressions and their coefficeints of determination (r2) for the data in
fig. 4.8. Here, t denotes the walltime of the simulation (in seconds) and n the number of
spatial basis functions.

4.5.3 Maxwell-Bloch

Finally, as a pragmatic test of our algorithm, we revisit a sample radiating quantum dot
system like those detailed in ??. We consider a random arrangement of 1024 quantum dots
within a 0.2 µm cube. The quantum dots have a resonant transition energy of 1500 meV
(thus ωmax (cid:3) 1500 meV/) and a dipole moment of d (cid:3) 10 e a0ˆz where e and a0 denote the
electron charge and Bohr radius, respectively. We employ a rotating-wave framework so as
to take ∆t (cid:3) 50ωmax, and we use ∆s (cid:3) c/(10ωmax) for our auxiliary grid spacing. Figure 4.9
shows the evolution of the off-diagonal matrix elements for 25 randomly sampled quantum
dots with both a fully-direct and FFT-accelerated calculation.

53

)
1
0
ρ
(
n
o
i
t
a
z
i
r
a
l
o
P

0.4
0.2
0
−0.2
−0.4
0.4
0.2
0
−0.2
−0.4
0.4
0.2
0
−0.2
−0.4
0.4
0.2
0
−0.2
−0.4
0.4
0.2
0
−0.2
−0.4

0

10

20

0

10

20 0

10

20 0

10

20 0

10

20

Time (ps)

Figure 4.9: Comparison of a Maxwell-Bloch quantum dot simulation with direct (red)
and FFT-accelerated (blue) propagation strategies. Here, we show the off-diagonal matrix
elements of 25 random quantum dots (out of 1024) as a function of time. The quantum
dots experience an incident π-pulse and their interactions induce significant secondary
oscillations (see the second-to-last row). The FFT-accelerated algorithm developed here
readily captures these effects, making it eminently suitable for large-scale quantum dot
simulations.

54

CHAPTER 5

CONCLUSIONS AND FUTURE PERSPECTIVE

“Why, sometimes I’ve believed as many as six impossible things before breakfast!”

—The White Queen, Through the Looking Glass

This thesis serves three main purposes. First, it introduces the concept of active
media—physical systems governed by ordinary differential equations that couple together
via radiative processes—and motivates the difficulty in examining these systems at a
macroscopic or continuum level. Second, it provides a discussion of numerical solution
strategies based around time-domain integral equations. Finally, it introduces and details
the QuEST software package. This chapter summarizes the main findings of each of these
parts and offers some perspective on potential future research.

5.1 Conclusion

In chapter 2 we presented a computational analysis of acoustically-driven microsphere
systems. By considering only small microspheres we cast the problem in terms of an
incompressible fluid where velocity potentials describe pressure distributions. Acoustic
radiation incident on the system scatters from rigid spherical objects, leading to pressure
variations over the objects’ surfaces; this variation necessarily sets the objects into motion,
leading to further potential/pressure perturbations within the system. By expressing the
velocity potential on the surface of each microsphere in terms of spherical harmonics, we
obtain a matrix equation relating the velocity of each microsphere to the potential on every
other microsphere. From this, we compute the net force on every body which we then
integrate numerically, noting that directly integrating the force causes numerical instabilities
due to an atypical acceleration-dependent force. To remedy this, we note that this force
serves only to impede the motion of each microsphere, thus we treat it as an additional
mass term. We find that an incident acoustic pulse predominantly serves to translate a

55

collection of microspheres a fixed distance along ˆk—the pulse’s wavevector—though the
scattering effects also produce system-scale expansions and contractions depending on the
microspheres’ initial configuration.

Chapter 3 departs from kinematic systems to develop a semiclassical description of
coupled quantum dots. Here, each of N quantum dots evolves according to the Liouville
equation in response to an ambient electric field. From this, we compute the expectation
value of a transition dipole moment. This acts as a source term in Maxwell’s equations
and gives rise to secondary radiation that couples the evolution of each quantum dot. By
employing an integral equation framework to describe this radiation, we recover fields
with a high degree of accuracy at both extremely short and extremely long length scales.
We find that pairs of exceptionally close quantum dots have a greatly diminished response
to an incident pulse; we attribute this effect to a dynamical detuning that brings the
pair out of resonance with the external driving field. Larger assemblages of quantum
dots, however, display a much richer behavior: we observe spatially-localized regions of
increased and diminished polarization as well as multiple frequency generation in the
nonlinear polarization response of particular quantum dots.

The integral equation methodology of chapter 3 has an implicit O(n2) spatial complexity
because of the pairwise interactions between quantum dots.
In practice, this limits
simulations to several tens of thousands of particles. To overcome this bottleneck, chapter 4
develops a FFT-based acceleration scheme that exploits the translational invariance of the
radiation kernel to greatly reduce the spatial complexity of the calculation. Moreover, our
acceleration scheme retains a strong-form integral formulation and makes no assumptions
about the underlying computational basis functions to effect propagation. We give
preliminary results that depict the accuracy and speed of the method and demonstrate
that we can recover fields under the action of an arbitrary linear operator to the precision
of the underlying machine.

56

5.2 Discussion

5.2.1 Advanced dissipators

Designed to provide a semblance of spontaneous emission within a semiclassical framework,
the two-level dissipator in eq. (3.2b) works to remove internal energy from each quantum
dot system and return it to the ground state. A semiclassical radiation model will never
fully describe spontaneous emission due to the lack of a quantized electromagnetic field,
but more robust operators that couple the dissipation terms of multiple quantum dots
together may more accurately model cooperative phenomena such as four-wave mixing or
superradiance. To derive the structure of one such advanced dissipator, consider first a
two-quantum dot system with particles A and B. Without rigorous derivation, a coupled
dissipator evolves the particles according to

Assuming ˆρ ≈ ˆρA⊗ ˆρB and tracing over B to remove its degrees of freedom,

where

(cid:2) ˆH 0

a + ˆH 0

b

d ˆρ
dt

(cid:3)

−i


(cid:1).

n − ˆσ+

m

ˆσ−
n

ˆρ − ˆρ ˆσ+

m

ˆσ−
n

, ˆρ(cid:3) +




m∈{A,B}
n∈{A,B}

γmn
2

m

ˆρ ˆσ+

(cid:0)2 ˆσ−
(cid:3) + ˆDAA + ˆDAB
(cid:170)(cid:174)(cid:172)

−i


(cid:3)

, ˆρA

A

d ˆρA
dt

(cid:2) ˆH 0
ˆDAA (cid:3) γAA(cid:169)(cid:173)(cid:171) −ρA11
2 (cid:169)(cid:173)(cid:171)−ρA10ρB01 − ρB10ρA01 ρA11ρB10 − ρB10ρA00
3 ω3|d|2(cid:16)
J0(rmn) + P2(cos θmn)J2(rmn)(cid:17)

ρB01ρA11 − ρA00ρB01

−ρA10/2
ρA11

ρB01ρA10 + ρA01ρB10

−ρA01/2

ˆDAB (cid:3)

γmn (cid:3)

γAB

4

(cid:170)(cid:174)(cid:172)

(5.1)

(5.2)

(5.3a)

(5.3b)

(5.3c)

Finally, J(cid:96) and P(cid:96) denote the Bessel function and Legendre polynomial of order (cid:96),
respectively. Note that the first two terms on the right side of eq. (5.2) give the evolution
ˆDAB dissipator, then, has terms proportional
equation defined in chapter 3. The additional

57

to the product of A and B’s coherences/off-diagonal terms with a long-range inverse-r
scaling due to the J0 Bessel function. These terms afford a mechanism for A and B to
influence each other’s dynamics more directly than through an intermediate radiation field
alone.

In a many-particle system, eq. (5.2) expands to include pairwise dissipators that perturb

(cid:2) ˆH 0

−i





(cid:3) +

the evolution of every particle due to the assumption of a rank-deficient ˆρ. Thus,

d ˆρm
dt

ˆDmn.

γmn

(cid:3)

m, ˆρm

(5.4)
Evaluating eq. (5.4) for each of N then scales as O(N2), making it infeasible to use in
large simulations. Fortunately, the kernel-agnostic nature of QuEST makes it relatively
straightforward to develop FFT-accelerated interactions of the form in eq. (5.3) alongside
those used to describe electromagnetic fields.

n

5.2.2 Micromagnetics

The Landau-Lifshitz-Gilbert (LLG) equation,

dM
dt

(cid:3) −γ(cid:48)M × H − λM × M × H

where

(5.5)

(5.6a)

(5.6b)

γ(cid:48)

(cid:3)

λ (cid:3)

1 + γ2η2M2
s

γ0

γ2
0 η

1 + γ2η2M2
s

γ0 (cid:3)

g|e|
2mec

details the evolution of magnetization, M, in response to an applied magnetic field, H [2].
In these equations, γ0 denotes the electron gyromagnetic ratio

,

(5.7)

η represents a material-depnedent damping parametre, and g, e, me, & c stand for the Landé
g-factor, elementary charge, electron mass, and speed of light, respectively. Simulating the
evolution of the Landau-Lifshitz-Gilbert equation at the level of individual particles stands
ready to offer novel perspectives into spintronics or other mesoscale systems.

58

5.2.2.1 Technical details

Immediately, eq. (5.5) draws a large number of parallels with the material equation of
chapter 3. Apart from the obvious first-order and vector nature of both equations, the
dual nature of electromagnetism means that M (as a source in Maxwell’s equations)
relates to H in the exact same way that P relates to E. By simply exchanging some of the
physical constants, then, eq. (3.3) also describes radiation emanating from a magnetization
distribution [64].

Solving eq. (5.5) within a semiclassical radiation framework requires a little more care
than an equivalent Bloch problem, however. Equation (3.1) describes the evolution of a
quantum wavefunction—a vector in “Bloch space”—that generates an electromagnetic
source distribution through application of an operator (see the discussion surrounding
eqs. (3.2) and (3.11) on pages 26 to 29 for an illustration of this point). As a result, the
dynamics of eq. (3.1) describe the evolution of a dipole’s amplitude but not its orientation
and we may employ separate spatial and temporal discretizations to effect a numerical
analysis. Moreover, this allow for systems with zero initial polarization (P (cid:3) 0 in an
eigenstate of the unperturbed Hamiltonian, e.g. ˆρ (cid:3) |0(cid:105)(cid:104)0| or ˆρ (cid:3) |1(cid:105)(cid:104)1|), thus the external
pulse induces all of the dynamical behavior. This makes for a far easier simulation as
sources/fields before the start of the simulation do not require any specification.

On the other hand, eq. (5.5) determines the trajectory of an electromagnetic source
vector in both (real, Cartesian, directional) space and time with no ability to decouple the
dimensions beyond

M(t) (cid:3) mx(t)ˆx + my(t)ˆy + mz(t)ˆz.

(5.8)

QuEST accommodates histories of any numeric type (scalar or tensor, real or complex)
through careful templating, and so building a propagator for eq. (5.5) does not present
much difficulty aside from some extra index tracking. Equation (5.5) also has an implicit

59

normalization condition, however, in that

d|M|
dt

(cid:3) 0.

(5.9)

This presents a larger computational challenge as it implies particles/magnetic do-
mains/“spins” always interact—even before the start of a simulation in a manner not
dissimilar to molecular dynamics simulations. The mechanism behind the predictor/cor-
rector algorithm—updating a source distribution and then reincorporating that into the
electromagnetic fields before advancing the timestep—will accommodate these interactions
without issue but meaningful simulations will require an equilbiration period of order
O(rmax/c) before collecting data.

60

APPENDICES

61

APPENDIX A

QUEST MANUAL

A.1 Design philosophy

The collection of formulas and algorithms in this thesis provide a mathematical founda-
tion upon which one can implement a simulation tool. This section gives some details on
v0.2.1 of the Quantum Electromagnetics Simulation Toolbox (QuEST) reference imple-
mentation, available at [35]. Written in C++, QuEST has minimal external dependencies,
using only the Boost Multiarray and Program Options libraries [1], the Eigen linear algebra
library [39], and FFTW [31].

As a scientific endeavor, QuEST largely emphasizes reproducibility, ease of extension,
and execution speed. Reproducibility demands modularity to allow well-tested pieces to
combine in a recordable manner reminiscent of experiments, thus the main functionality
of QuEST resides within a library of independently-constructed objects that interface

QuEST

Radiation

Pulse

ODE

Interactions

Direct

Integration

Interpolation

RHS functions

History

AIM

FFT

Figure A.1: Hierarchy of the major pieces of QuEST. While a true class diagram would
depict the inheritance relationships of the objects that comprise the QuEST library, These
blocks roughly indicate the interdependence of the program pieces that form QuEST, akin
to a class hierarchy (though avid software developers should not take this literally).

62

with each other. These objects have appropriate constructors to set the parameters of
the simulation and slot together to build an executable. For example, to isolate the ODE
solver’s behavior from the tricky numerics of the propagation code, the ODE solver only
requires the definition of a right-hand side (an RHS) (see section A.2). Coding an RHS
with a well-known solution, such as f (cid:48)(x) (cid:3) −λ f(x), can validate the numerics of the
integration step independent of the propagation calculations.

Similarly, making QuEST easy to extend necessitates encapsulating mathematical
statements as tightly as possible with these objects. As the project grows to encapsulate
new phenomena/models, only the pieces that directly represent the mathematics need to
change, not how they fit into the rest of the simulation framework. Unfortunately, these
factors often lead to seemingly-unintelligible code. What follows contains a brief discussion
of the major pieces of QuEST and why I designed them the way I did.

A.2 Descriptions

cmplx An alias for std::complex<double> with an identical bit structure to FFTW’s

fftw_complex type (hence the use of static_cast when interfacing with FFTW).

Configuration This struct effectively contains “global” simulation parameters—the speed
of light, number of particles, total time, etc. Most objects do not read these values
directly, instead preferring to have them as a parameter in their constructor to enhance
modularity. The parse_configs function uses the Boost Program Options [1] library
to provide a clean, error-checked way of reading these values from specified files.

QuantumDot Collection of parameters necessary to define a quantum dot (position, tran-
sition frequency, decay times, and dipole moment). Also provides a method to
calculate the right-hand side of eq. (3.4) given an input density matrix and electric
field. Because these right-hand sides have instance-dependent parameters (transi-
tion frequency, T1 and T2, etc.), std::bind produces the appropriate functors for
constructing Integrator::RHS instances.

63

Pulse Mathematical description of a Gaussian pulse. Provides a method to evaluate a
pulse at a given spacetime coordinate, as well as methods to read/record Pulse
objects from/to a file.

UniformLagrangeSet Tool to evaluate Lagrange polynomials and their derivatives. Upon
construction, each instance allocates a workspace table (to avoid reallocation overhead)
to contain the evaluated interpolants and their derivatives. This object uses the defini-
log f(x) (cid:3) f (cid:48)(x)/ f(x),
tion of Lagrange polynomials and logarithmic derivatives, d
dx
to efficiently fill each column (the interpolant and two derivatives) of the table
simultaneously.

TransformPair Struct to manage FFTW plans (pointers). Constructed with a forward and
backward plan and automatically de-allocates both in its destructor to avoid memory
leaks.

Interaction Abstract base-class to provide an interface for calculating an interaction
for each QuantumDot in an array. Each Interaction contains a shared pointer to a
std::vector<QuantumDot> and allocates a complex-valued Eigen row-vector (once,
to avoid reallocation overhead) upon construction. This row-vector serves as a
workspace in which to place interaction values in the evaluate method before
returning. Ultimately subclasses into PulseInteraction, DirectInteraction, and
AimInteraction.
PulseInteraction Calculates and returns d·Einc/ at the location of each QuantumDot

at a specified time.

DirectInteraction (subclass of HistoryInteraction) Calculates and returns d ·
Erad/ at the location of each QuantumDot (where Erad depends on the history of
every other QuantumDot due to retardation effects). This object stores the space/-
time interaction matrix elements between each QuantumDot directly, assuming
reciprocity and no self interactions (effectively the lower triangle of the matrix).

64

Uses Lagrange interpolation to calculate quantities “between” timesteps as well
as temporal derivatives.

AIMInteraction (subclass of HistoryInteraction) Calculates d · Erad/ between
well-separated sources (so as to avoid incurring undue approximation error).
Uses the vastly more efficient AIM (FFT) framework to reduce the storage and
accelerate the matrix-vector products. Uses Lagrange interpolation to calculate
quantities “between” timesteps as well as temporal derivatives. WARNING:
you must call AIMInteraction.evaluate with sequential arguments (starting
at zero, corresponding to sequential timesteps) due to the way this object fills
internal tables. This keeps the indexing correct as later calls to evaluate depend
on the results of previous calculations.

Grid Class to construct and manage a representation of the auxiliary AIM
grid. Provides methods to convert between grid coordinates (nx, ny, and
nz integers that index a node in the grid), a grid index (linearization of
the grid coordinates to a single integer index starting at zero), and spatial
coordinates (r (cid:3) xˆx + y ˆy + zˆz real-valued coordinates in three-dimensional
space).

ExpansionFunction Complicated black magic sorcery. Stay away.
ExpansionTable Table of dimension [number of sources][number of points each
source expands into] containing the thirteen (∂0, ∂x, . . . , ∂yz, ∂zz) expansion
weights relating each source to its immediate (i.e. not trivially zero) auxiliary
grid points.

LeastSquaresExpansionSolver Utility class to build and solve the least-squares
system in eq. (4.9). Uses a named constructor to avoid creating a persis-
tent instance that implicitly depends on a non-const Grid (only the static
get_expansions function can instantiate a LeastSquaresExpansionSolver
and thus will call its destructor when the function returns).

65

PredictorCorrector Implementation of the exponentially-fitted predictor/corrector scheme
in [33]. Implemented as a template class so as to allow solution of arbitrary ODEs
(scalar/vector, real/complex, etc.).

Weights Assistance class for solving the matrix equation that determines the predic-
tor/corrector weights. Putting the code for these matrices into its own object
automatically frees the resources used to store them (by automatically calling
the object’s destructor when it goes out of scope) after calculating the required
values.

History Wrapper around a three-dimensional Boost::multi_array to store data for
coupled first-order ODEs at fixed solution points (timesteps). Stores both a solution
and its derivative (the “RHS”) for each of Ns ODEs at Nt solution points, indexed by
[equation index][time index][derivative index]. Provides integer-equivalent typedefs
to assist with indexing. As the predictor/corrector integration scheme requires a
history to solve for any timepoint (including the first), these arrays extend the time
axis backwards to negative indices so as to always start at t (cid:3) 0.

66

APPENDIX B

THE PREDICTOR/CORRECTOR INTEGRATION SCHEME

B.1 Motivation

Given an ordinary differential equation,

(cid:3) f(t, x(t))

dx
dt

(B.1)

and some initial condition at t0, we may compactly write the unknown function at any
discrete time ti ≡ i ∆t as

Þ ti
Þ ti

ti−1

ti−1

x(ti) (cid:3) x(ti−1) +
≡ x(ti−1) +

(cid:48)

dx
dt(cid:48) dt
f(t
(cid:48), x(t

(cid:48))) dt

(cid:48)

(B.2)

Numerical integration schemes advance a solution from ti−1 to ti by approximating this
integral, often as a weighted sum of the integrand evaluated at intelligently chosen points
within the interval1. This process—particularly evaluating the integrand—becomes vastly
more complicated through the introduction of retardation effects, effectively turning
eq. (B.1) into

(cid:3) f(t, x(t), x(t − τ1), x(t − τ2), . . .).

(B.3)

dx
dt

Causality imposed by the propagation operators only requires τj > 0, thus eq. (B.3) can
relate a derivative to quantities arbitrarily far back in time, often between the discrete
ti/known solution values. Further problems arise when τj < ∆t; in such cases, the unknown
x(ti) depends on an unknown x(ti − τj) which depends on an unknown x(ti − 2τj) etc., at
least one of which lies “ahead” of the last known quantity, x(ti−1). While adjusting the
timestep to require τj (cid:62) ∆t makes the system perfectly causal, it also makes the timestep
1Many numerical integrators use much more sophisticated methods than a weighted

sum, but the idea of intelligently evaluating the integrand remains the same.

67

prohibitively small for dense systems. These constraints automatically exclude (or at
least severely limit the utility of) popular ODE algorithms and packages—the difficulty
of re-evaluating interactions prohibits variable-timestep integrators and integrators that
require midpoint evaluations (such as the popular RK4 method) cannot accommodate
small retardation factors without an exceptionally small ∆t.
Instead, we turn to the
exponentially-fitted predictor/corrector algorithm detailed in [33]. As a multistep method,
the predictor/corrector does not require any function evaluations outside of those on the
fixed timepoints, thus any causal interpolation scheme facilitates evaluating approximate
versions of all the x(t − τj) via

x(ti − τj) ≈ m


w(cid:96)x(cid:0)ti − ((cid:98)τj/∆t(cid:99) + (cid:96)) ∆t(cid:1).

(B.4)

(cid:96)(cid:3)0

(The particular interpolation scheme determines the w(cid:96).) For τj < ∆t, eq. (B.4) suggests
x(ti − τj) becomes a function of x(ti). Ordinarily, this would cause problems though the
predictor step of the algorithm ensures an approximate x(ti) exists before any interpola-
tion. The corrector step then works to correct this approximation in a convergent way,
reincorporating adjacent sources as a total solution evolves.

B.2 Integration coefficients

Let x(t) denote the solution to a given ordinary differential equation. We seek to

approximate x(t) as a linear combination of complex-valued exponentials such that

x(t) ≈ Nλ−1


(cid:96)(cid:3)0

w(cid:96)e λ(cid:96)t

(B.5)

for a given set of complex-valued λ(cid:96) arguments and w(cid:96) weights (never calculated explicitly).
As we wish to capture the behavior of physical systems—i.e. systems with exclusively
oscillating and decaying modes—we choose the λ(cid:96) to lie within the left half of the complex
plane. Moreover, assuming x(t) contains little power above some threshold ωmax, we select
the λ(cid:96) from within S (cid:3) {z ∈ C | Re{z} (cid:54) 0, |z| (cid:54) ωmax}: the left semidisk of radius ωmax
(fig. B.2a). Finally, rather than choosing λ(cid:96) from throughout S, the maximum modulus

68

5 ∆t
4 ∆t
3 ∆t
2 ∆t
1 ∆t

equations. If ∂tx•(t) (cid:3) f(cid:0)x•(t − |r••|) + x•(t − |r••|)(cid:1) then determining x•(t − |r••|) and

Figure B.1: Illustration of the “retardation problem” in solving coupled delay differential
x•(t − |r••|) requires some sort of interpolation scheme. Evaluating ∂tx•(t) in hopes
of obtaining x•(t) therefore requires knowledge of x•(t). The converse holds as well;–
evaluating ∂tx•(t) requires knowledge of x•(t) which poses problems in evaluating the
derivatives simultaneously. The predictor step of the integrator overcomes this issue by
providing a stable extrapolation of past quantities to provide a guess value of x•(t) and
x•(t).

Im λ

ωmax

Im λ

ωmax

Re λ

Re λ

(a)

(b)

Figure B.2: Selection of the λn parameters used in determining the predictor/corrector
coefficients.

69

principle ensures that arguments taken from the boundary of S in fig. B.2b will incur the
maximum approximation error in eq. (B.5). As a result, an approximation using λ(cid:96) chosen
from the boundary of S accurately recovers the behavior of all modes with arguments in S.
In particular, we choose the λ(cid:96) to have an equal spacing about the perimeter of S (fig. B.2b).
If x(t) and ∂tx(t) have known values at equidistant points t0, t1, . . . , tW−1 on the interval
[−1, 1] (thus ∆t (cid:3) 2/(W − 1)), we predict the value of x(tW) as a linear combination of past
values and derivatives,

x(tW) ← W−1


j(cid:3)0

p jx(t j) + p j+W ∂tx(t j).

Inserting eq. (B.5) into eq. (B.6) gives the matrix equation

where

Ap (cid:3) b

0 (cid:54) i < Nλ; 0 (cid:54) j < W
0 (cid:54) i < Nλ; W (cid:54) j < 2W

e λit j ,

λie λit j−W ,

Aij (cid:3)

bi (cid:3) e λitW .

(B.6)

(B.7)

(B.8a)

(B.8b)

(B.9)

Solving this equation for p, then, produces the prediction coefficients in eq. (B.6). Accord-
ingly, we use a minimum-norm least-squares procedure to determine p. Having conjured
an (approximate) value for x(tW), we may evaluate the prescribed ODE to obtain a value
for ∂tx(tW). Writing

x(tW) ←(cid:169)(cid:173)(cid:171)W−1


j(cid:3)0

c jx(t j) + cW +j ∂tx(t j)(cid:170)(cid:174)(cid:172) + c2W ∂tx(tW),

this new equation reincorporates derivative information to refine the x(tW) approximation,
requiring only a means of determining the c j. Again, inserting eq. (B.5) gives

˜Ac (cid:3) b.

70

(B.10)

Here,

˜Aij (cid:3)

Aij,

λie λitW ,

0 (cid:54) i < Nλ; 0 (cid:54) j < 2W
0 (cid:54) i < Nλ; j (cid:3) 2W

(B.11)

b remains unchanged, and the same least-squares procedure determines c.

Having derived the predictor/corrector coefficients for h (cid:3) 2/(W − 1)—i.e. a uniform
spacing between W timepoints on the interval [−1, 1]—the coefficients determined above
require only a slight modification to accommodate systems with any ∆t. Noting that the
formulation thus far remains invariant under time translation and that the substitution
τ (cid:3) αt (equivalently α (cid:3) ∆t/h) turns

∂ϕ(t)
∂t

(cid:3) λϕ(t)

(B.12a)

into

(B.12b)
we need only scale the coefficients multiplying ∂tx(t j) by α to adjust for other step sizes.

∂ϕ(τ)
∂τ (cid:3)

λ
α

ϕ(τ),

71

APPENDIX C

CHEBYSHEV POLYNOMIALS

1

0.5

0

−0.5

−1
−1.2 −1 −0.8 −0.6 −0.4 −0.2

0

0.2

0.4

0.6

0.8

1

1.2

Figure C.1: Chebyshev polynomials T0(x) through T6(x). The roots of Tm+1(x) give the
sample points for an interpolation scheme of order m on the interval [−1, 1].

The Chebyshev polynomials of the first kind, defined by

Tn(x) (cid:3) cos(n arccos(x))

(C.1)

and shown in fig. C.1, form an orthogonal (though not orthonormal) sequence well-suited
to interpolation. These polynomials have a continuous orthogonality relation over the

interval −1 (cid:54) x (cid:54) 1 with respect to w(x) (cid:3)(cid:0)1 − x2(cid:1)−1/2, giving

Such an orthogonality motivates {Tk(x) : k (cid:3) 0, 1, . . . , m} as a basis for Pm—the polynomi-
als through order m—on the interval −1 (cid:54) x (cid:54) 1. To construct a polynomial interpolation

Þ 1

−1

Ti(x)Tj(x)
√
1 − x2

dx (cid:3)

0
π/2

π

i (cid:44) j
i (cid:3) j (cid:44) 0
i (cid:3) j (cid:3) 0

(C.2)



72

(C.3)

we require samples of the function at m + 1 interpolation nodes1 {xλ : λ (cid:3) 0, 1, . . . , m}
where f(xλ) (cid:3) Pm(xλ). By choosing sample points corresponding to Chebyshev nodes,

λ (cid:3) 0, 1, 2, . . . , m,

(C.4)

of an arbitrary function f(x) in this basis such that
f(x) ≈ Pm(x)

≡ m


i(cid:3)0

ciTi(x),

(cid:18) π(λ + 1/2)

(cid:19)

m + 1

xλ (cid:3) − cos

the Chebyshev polynomials satisfy [32]

m


λ(cid:3)0

Ti(xλ)Tj(xλ) (cid:3)



0
(m + 1)/2
m + 1

i (cid:44) j
i (cid:3) j (cid:44) 0
i (cid:3) j (cid:3) 0

(C.5)

m

m


λ(cid:3)0

in addition to the continuous relationship in eq. (C.2). Inserting eq. (C.3) into eq. (C.5) then
gives

Ti(xλ) f(xλ) (cid:3)

Ti(xλ)Pm(xλ) (cid:3)

Ti(xλ)Ti(xλ)

(C.6)

ci

m


i(cid:3)0

m


λ(cid:3)0

m


λ(cid:3)0

thus

ci (cid:3)

(C.7)
Finally, having obtained a Chebyshev approximation of f(x) on the interval −1 (cid:54) x (cid:54) 1,
we may effect a Chebyshev approximation of f(x) on the interval a (cid:54) x (cid:54) b by letting

λ(cid:3)0

Ti(xλ) f(xλ)

αi (cid:3) 2 − δi0.

αi
m + 1

y ≡ x − (b + a)/2
(b − a)/2
and constructing a Chebyshev approximation in y.

(C.8)

1For convenience, Latin letters index functions while Greek letters index points.

73

(cid:0)1 + 16x2(cid:1)−1/2

Lagrange
Chebyshev

1

0.5

0

−0.5

−1.2 −1 −0.8 −0.6 −0.4 −0.2

Figure C.2: Interpolation of a pathological function, f(x) (cid:3)(cid:0)1 + 16x2(cid:1)−1/2, using Chebyshev

and equally-spaced Lagrange interpolation schemes. Because the Chebyshev scheme
clusters sample points in the tails of the interpolation window, it does not suffer from
high-order ringing effects (Runge’s phenomenon).

0.6

0.8

0.2

0.4

0

1

1.2

C.1 Notes on the Chebyshev polynomials

C.1.1 Derivatives
The Chebyshev polynomials provide two means of approximating f (cid:48)(x) given samples of
f(x) though both have numerical drawbacks. Given that we wish to repeatedly evaluate the
interpolating polynomial at a collection of static points, it becomes prudent to pre-evaluate
and cache Ti(x) in eq. (C.3) for speed. The recurrence

n(x) (cid:3) n[Tn−1(x) − xTn(x)]
(cid:48)

(C.9)

(cid:0)1 − x2(cid:1)T

analytically relates Chebyshev polynomials to their derivatives, and so we may equivalently
n(x) to repeatedly evaluate a derivative. Unfortunately, this becomes problematic
cache T(cid:48)
for x → ±1 as the 1 − x2 term can cause a catastrophic loss of precision due to division by
(near) zero. Instead,

(cid:48)
i−1 (cid:3) c

(cid:48)
i+1 + 2ici

c

(cid:48)
m+1 (cid:3) c

(cid:48)
m (cid:3) 0

c

(C.10)

74

adjusts the Chebyshev coefficients to give an approximated derivative in terms of the
original basis. While this expression does not suffer from numerical cancellation issues, it
incurs a large amount of overhead when evaluating many approximations. Moreover, this
overhead grows significantly in higher-dimensional systems with vector derivatives.

In practice, QuEST makes use of the former strategy of differentiating the “evaluating
functions” though without the recurrence of eq. (C.9). Instead, functions for the Cheby-
shev polynomials and their derivatives through a specified order, provide the necessary
evaluations. Moreover, these functions contain the requisite polynomials expressed in
Horner form to improve both speed and precision [52]. While this method does not
readily generalize to interpolations of arbitrary order (as it would require prohibitively
long enumerations of functions), it nevertheless facilitates an efficient derivative evaluation
without unduly sacrificing clarity.

C.1.2 Higher-dimensional approximations
The mechanics for constructing an approximation to f(x) extend naturally to higher-
dimensional and vector-valued systems. In such systems, the Chebyshev approximation
becomes a tensor polynomial, and the summations in eqs. (C.3) and (C.7) become mondo
tensor contractions, i.e.

f(r) ≈ cijkTi(x)Tj(y)Tk(z),
cijk (cid:3)

αii(cid:48)α j j(cid:48)αkk(cid:48)

(m + 1)3 Tii(cid:48)(xλ)Tj j(cid:48)(yµ)Tkk(cid:48)(zν)f(xλ, yµ, zν)

(C.11b)
with an implied summation over the bound indices i(cid:48), j(cid:48), k(cid:48) and λ, µ, ν. QuEST makes use
of the Tensor Algebra Compiler (taco) [51] to generate the equivalent source-expression for
eq. (C.11b) with dense tensors (listing C.1).

(C.11a)

Listing C.1: C-source expression for eq. (C.11b).

for (int32_t boxeval = 0; boxeval < num_boxes; boxeval++) {

double tbox = norm;

75

for (int32_t ialphas = 0; ialphas < size; ialphas++) {

int32_t pcoef2 = boxeval * size + ialphas;

double ti = tbox * alphas_[ialphas];

for (int32_t lambdachebS = 0; lambdachebS < size; lambdachebS++) {

int32_t pchebS2 = ialphas * size + lambdachebS;

int32_t peval2 = boxeval * size + lambdachebS;

double tlambda = ti;

double tlambda0 = chebS[pchebS2];

for (int32_t jalphas = 0; jalphas < size; jalphas++) {

int32_t pcoef3 = pcoef2 * size + jalphas;

double tj = tlambda * alphas_[jalphas];

double tj0 = tlambda0;

for (int32_t muchebS = 0; muchebS < size; muchebS++) {

int32_t pchebS20 = jalphas * size + muchebS;

int32_t peval3 = peval2 * size + muchebS;

double tmu = tj;

double tmu0 = tj0;

double tmu1 = chebS[pchebS20];

for (int32_t kalphas = 0; kalphas < size; kalphas++) {

int32_t pcoef4 = pcoef3 * size + kalphas;

double tk = tmu * alphas_[kalphas] * tmu0 * tmu1;

for (int32_t nuchebS = 0; nuchebS < size; nuchebS++) {

int32_t pchebS21 = kalphas * size + nuchebS;

int32_t peval4 = peval3 * size + nuchebS;

double tnu = tk * chebS[pchebS21];

for (int32_t dimeval = 0; dimeval < 3; dimeval++) {

int32_t peval5 = peval4 * 3 + dimeval;

76

int32_t pcoef5 = pcoef4 * 3 + dimeval;

coef[pcoef5] = coef[pcoef5] + tnu * eval[peval5];

}

}

}

}

}

}

}

}

77

APPENDIX D

CONVENTIONS

D.1 Fourier transforms

The definition of the Fourier transform leaves some room to adopt conventions. I write

Þ ∞
the fully-generalized Fourier transform f(t) ↔ F(ω) pair as
−∞ f(t)e is2ωt dt
Þ ∞
−∞ F(ω)e
−is2ωt dω

(cid:115) |s2|
(cid:115) |s2|

(2π)1+s1

(2π)1−s1

F(ω) (cid:3)

f(t) (cid:3)

(D.1a)

(D.1b)

where the s1, s2 parameters allow for different normalization and frequency conventions.
Of principal importance, eq. (D.1b) implies

(D.2)
The most common physics and engineering conventions take {s1, s2} (cid:3) {1, 1} and {s1, s2} (cid:3)
{1,−1}, though myriad other discipline-specific conventions exist.

∂t f(t) ↔ −iωs2F(ω).

D.2 Maxwell’s equations

An unfortunate ambiguity in the prescription of electromagnetic units allows for unit
systems that differ by more than simple scaling factors (mega-, kilo-, centi-, etc.)—the
placement of physical constants in equations changes as well. In a unit-independent format,
Maxwell’s equations become [45]

k2α
k1

∂E
∂t

∇ · E (cid:3) 4πk1ρ
∇ × B (cid:3) 4πk2αJ +
∇ · B (cid:3) 0
∇ × E (cid:3) −k3

∂B
∂t

78

(D.3a)

(D.3b)

(D.3c)

(D.3d)

with k1, k2, α, and k3 as system-dependent constants. One need only specify two of the
four constants to pin down the unit system, however; in particular

From eq. (D.3c), B (cid:3) ∇ × A. Then, using eq. (D.3d), it becomes apparent that

and k3 (cid:3)

1
α

.

(D.4)

k1
k2k3α (cid:3) c2

∇ ×(cid:18)

E + k3

∂A
∂t

(cid:3) 0

(cid:19)

∂A
∂t

.

so I may take

E (cid:3) −∇ϕ − k3

Thus, in electorstatic systems, E (cid:3) −∇ϕ regardless of units. By eq. (D.3a), then,

and thus

as

∇2(cid:18)

ϕ(r) (cid:3)

1
|r − r(cid:48)|

∇2ϕ (cid:3) −4πk1ρ

Þ k1ρ(r(cid:48))
|r − r(cid:48)| d3r(cid:48)
(cid:19)

(cid:3) −4πδ(r − r(cid:48)).

This defines the Laplace-kernel Green’s function as
g(r) ≡ 1
|r| .

Analogously, the Helmholtz-kernel and Wave-kernel Green’s functions become

g(r; t) ≡ δ(t − |r|/c)

|r|

g(r; k) ≡ e ik|r|
|r|

D.2.1 Vector wave equation

D.2.1.1 Operator form

From eqs. (D.3b) and (D.3d),

∇ × ∇ × E +

1
c2

∂2E
∂t2 (cid:3) −4πk2

∂J
∂t

.

79

(D.5)

(D.6)

(D.7)

Taking the divergence of both sides and removing one time derivative,

As ∇ × ∇ × A (cid:3) ∇(∇ · A) − ∇2A, the time derivative of eq. (D.7) becomes

(cid:3) −4πk2∇ · J.

1

∂t

c2∇ · ∂E
(cid:19) − ∇2 ∂E

∂t

∇(cid:18)∇ · ∂E

∂t

Immediately, the left-hand side gives the “wave operator” acting on ∂tE and thus

∂2J
∂t2

(cid:19)

.

+

1
c2

∂3E
∂t3 (cid:3) −4πk2
(cid:18) ∂2J
∂t2 − c2∇∇ · J
∂t2 − c2∇∇·(cid:19)
(cid:18) ∂2
(cid:18) ∂2
∂t(cid:48)2 − c2∇(cid:48)∇(cid:48)·(cid:19)

J(r, t)

therefore

∇2 ∂E
∂t

− 1
c2

∂3E
∂t3 (cid:3) 4πk2

E(r, t) (cid:3) g(r; t) ∗ −k2

∂
∂t

from which it follows

E(r, t) (cid:3) −k2(cid:120) δ(tR − t(cid:48))
|r − r(cid:48)|

where tR ≡ t − |r − r(cid:48)|/c and ∂tP ≡ J.

D.2.1.2 Dyadic form

P(r(cid:48), t

(cid:48)) d3r(cid:48) dt

(cid:48)

(D.8)

In addition to eq. (D.8), one may also construct an equivalent propagator without any
spatial derivatives in the resulting expression. This form of the propagator has advantages
in describing E(r, t) when derivatives become ill-defined (for example, fields arising from
point sources).

To begin, we note that for time-harmonic P(r, t) (cid:3) P(r)e−iωt eq. (D.8) becomes

(D.9)
after suppressing the time-harmonic function and taking k (cid:3) ω/c. All of the terms in this
equation commute due to the linearity of the operators, thus

E(r) (cid:3) −k2

Þ e ik|r−r(cid:48)|
|r − r(cid:48)|(cid:0)(−iω)2 − c2∇(cid:48)∇(cid:48)·(cid:1)P(r(cid:48)) d3r(cid:48)
E(r) (cid:3) −k2 (−iω)2Þ(cid:18)
≡ −k2(−iω)Þ

(cid:19) e ik|r−r(cid:48)|
|r − r(cid:48)| · P(r(cid:48)) d3r(cid:48)
I +
Gd(r − r(cid:48); ω) · P(r(cid:48); ω) d3r(cid:48)

∇∇
k2

(D.10)

80

which defines the dyadic Green’s function
Gd(r; ω) (cid:3)

where I represents the identity dyad.

(cid:18)

I +

∇∇
k2

(cid:19) e ikr

r

,

(D.11)

In spherical coordinates, Gd(r) (cid:3) Gd(r) and the spatial derivatives become particularly

easy to compute. Accordingly,

∇∇ e ikr
r

(cid:19)
(cid:3) ∇(cid:18) e ikr
r2 (ikr − 1)(cid:19)
(cid:3) ˆr∇(cid:18) e ikr
r2 (ikr − 1)ˆr
− k2(cid:19)
(cid:20)(cid:18) 3
r2 − 3ik
(cid:19)
(cid:20)(cid:18) 3
(cid:20)(I − ˆrˆr) − (I − 3ˆrˆr)(cid:18) 1

3
(ik)r

k2r2 +

− 1

ˆrˆr +

(cid:3)

r

(ik)r

Gd(r) (cid:3)

(cid:3)

ˆrˆr +

I

r

+

e ikr

ik
r

(cid:18)− 1
(cid:19)
(cid:21) e ikr
r2 (ikr − 1)(∇ˆr)
r2 +
(cid:19)
(cid:18)− 1
(cid:19)(cid:21) e ikr
k2r2 − 1
(ik)r
1
.
k2r2

+ 1

+

r

I

(D.12)

(D.13)

(cid:21) e ikr

r

Inserting this expression back into eq. (D.11) then gives

r

c

+

(cid:34)

c2
r2

(−iω)
r

· P(r(cid:48), ω) d3r(cid:48)

E(r, ω) (cid:3) −k2

Þ(cid:20)(I − ˆrˆr)(−iω)2

Having eliminated the spatial derivatives, we can perform the inverse of the Fourier
transform that gave eq. (D.9). Inserting eq. (D.13) back into eq. (D.10) and replacing k with
ω gives

(D.14)
where ˆr (cid:3) (r − r(cid:48))/|r − r(cid:48)|. Exchanging −iω for a time derivative and restoring the temporal
convolution,
E(r, t) (cid:3) −k2(cid:120)

(cid:48)) d1t
(D.15)
The derivatives acting on the δ-functions make little sense, thus we can apply them equally
to P(r, t) again due to the linearity of the operators. As a result,

+ (I − 3ˆrˆr)(cid:18)
(cid:19)(cid:21) e ikr
+ (I − 3ˆrˆr)(cid:18) c∂t δ(tR − t(cid:48))
|r − r(cid:48)| + (I − 3ˆrˆr) ·(cid:18) c(cid:219)P(r(cid:48), tR)

E(r, t) (cid:3) −k2

c2P(r(cid:48), tR)
|r − r(cid:48)|3
where dots indicate derivatives with respect to t (or, equivalently, tR).

|r − r(cid:48)|2 +

Þ(I − ˆrˆr) · (cid:220)P(r(cid:48), tR)

c2δ(tR − t(cid:48))
|r − r(cid:48)|2

δ(tR − t(cid:48))
|r − r(cid:48)|

d3r(cid:48)

(D.16)

(I − ˆrˆr) ∂2

|r − r(cid:48)|2

·P(r, t

(cid:19)(cid:35)

(cid:19)

+

t

(cid:48) d3r(cid:48) .

81

BIBLIOGRAPHY

82

BIBLIOGRAPHY

[1] BOOST C++ Libraries. http://www.boost.org.
[2] Amikam Aharoni. Introduction to the theory of ferromagnetism. Oxford University Press,

Oxford New York, 2000.
John S Allen, Donovan J May, and Katherine W Ferrara. Dynamics of therapeutic
ultrasound contrast agents. Ultrasound in medicine & biology, 28(6):805–816, 2002.

[3]

[4] L. Allen and J.H. Eberly. Optical Resonance and Two-level Atoms. Dover books on physics

and chemistry. Dover, 1975.

[5] Philip W. Anderson. The question of classical localization: A theory of white paint?

Philosophical Magazine B, 52(3):505–509, 1985.

[6] F. Arecchi and R. Bonifacio. Theory of optical maser amplifiers. IEEE Journal of

Quantum Electronics, 1(4):169–178, July 1965.

[7] Kenta Asakura, Yasuyoshi Mitsumori, Hideo Kosaka, Keiichi Edamatsu, Kouichi
Akahane, Naokatsu Yamamoto, Masahide Sasaki, and Naoki Ohtani. Excitonic rabi
oscillations in self-assembled quantum dots in the presence of a local field effect. Phys.
Rev. B, 87:241301, Jun 2013.

[8] S. Azizoglu, S. Koc, and O. Buyukdura. Time domain scattering of scalar waves by
two spheres in free-space. SIAM Journal on Applied Mathematics, 70(3):694–709, 2009.
[9] Nicolas Bachelard, Rémi Carminati, Patrick Sebbah, and Christian Vanneste. Linear
and nonlinear rabi oscillations of a two-level system resonantly coupled to an anderson-
localized mode. Phys. Rev. A, 91:043810, Apr 2015.

[10] Pierre Barthelemy, Jacopo Bertolotti, and Diederik S. Wiersma. A Lévy flight for light.

Nature, 453:495 EP –, May 2008.

[11] Jean-Pierre Berenger. A perfectly matched layer for the absorption of electromagnetic

waves. Journal of Computational Physics, 114(2):185 – 200, 1994.

[12] E. Bleszynski, M. Bleszynski, and T. Jaroszewicz. AIM: Adaptive integral method for
solving large-scale electromagnetic scattering and radiation problems. Radio Science,
31(5):1225–1251, Sept 1996.

[13] M. J. Blomley, J. C. Cooke, E. C. Unger, M. J. Monaghan, and D. O. Cosgrove.
Microbubble contrast agents: a new era in ultrasound. BMJ, 322(7296):1222–1225, May
2001.

[14] M. J. Bluck and S. P. Walker. Time-domain bie analysis of large three-dimensional
electromagnetic scattering problems. IEEE Transactions on Antennas and Propagation,
45(5):894–901, May 1997.

83

[15] R. Bonifacio, P. Schwendimann, and Fritz Haake. Quantum statistical theory of

superradiance. i. Phys. Rev. A, 4:302–313, Jul 1971.

[16] R. Bonifacio, P. Schwendimann, and Fritz Haake. Quantum statistical theory of

superradiance. ii. Phys. Rev. A, 4:854–864, Sep 1971.

[17] F. Borgnis. Acoustic radiation pressure of plane compressional waves. Rev. Mod. Phys.,

25:653–664, Jul 1953.

[18] David C. Burnham and Raymond Y. Chiao. Coherent resonance fluorescence excited

by short light pulses. Phys. Rev., 188:667–675, Dec 1969.

[19] M J Caola. Solid harmonics and their addition theorems.

Mathematical and General, 11(2):L23, 1978.

Journal of Physics A:

[20] Qing-Hu Chen, Yu-Yu Zhang, Tao Liu, and Ke-Lin Wang. Numerically exact solution

to the finite-size dicke model. Phys. Rev. A, 78:051801, Nov 2008.

[21] Hank Childs, Eric Brugger, Brad Whitlock, Jeremy Meredith, Sean Ahern, David
Pugmire, Kathleen Biagas, Mark Miller, Cyrus Harrison, Gunther H. Weber, Hari
Krishnan, Thomas Fogal, Allen Sanderson, Christoph Garth, E. Wes Bethel, David
Camp, Oliver Rübel, Marc Durant, Jean M. Favre, and Paul Navrátil. VisIt: An
End-User Tool For Visualizing and Analyzing Very Large Data. In High Performance
Visualization–Enabling Extreme-Scale Scientific Insight, pages 357–372. Oct 2012.

[22] Claude Cohen-Tannoudji, Jacques Dupont-Roc, and Gilbert Grynberg. Photons and

atoms: introduction to quantum electrodynamics. Wiley Online Library, 1989.

[23] Xiaoyun Ding, Jinjie Shi, Sz-Chin Steven Lin, Shahrzad Yazdi, Brian Kiraly, and
Tony Jun Huang. Tunable patterning of microparticles and cells using standing
surface acoustic waves. Lab Chip, 12:2491–2497, 2012.

[24] Y. Ding, A. Forestier, and T. Ha Duong. A galerkin scheme for the time domain integral
equation of acoustic scattering from a hard surface. The Journal of the Acoustical Society
of America, 86(4):1566–1572, 1989.

[25] Alexander A. Doinikov. Mathematical model for collective bubble dynamics in strong
ultrasound fields. The Journal of the Acoustical Society of America, 116(2):821–827, 2004.
[26] Alexander A. Doinikov, Richard Manasseh, and Andrew Ooi. Time delays in coupled
multibubble systems (l). The Journal of the Acoustical Society of America, 117(1):47–50,
2005.

[27] A. A. Ergin, B. Shanker, and E. Michielssen. The plane-wave time-domain algorithm

for the fast analysis of transient wave phenomena. 41(4):39–52.

[28] A. A. Ergin, Balasubramaniam Shanker, and Eric Michielssen. Analysis of transient
wave scattering from rigid bodies using a burton–miller approach. The Journal of the
Acoustical Society of America, 106(5):2396–2404, 1999.

84

[29] A. Arif Ergin, Balasubramaniam Shanker, and Eric Michielssen. Acceleration of
integral equation based transient analysis of acoustic scattering phenomena using
the plane wave time domain algorithm. The Journal of the Acoustical Society of America,
106(4):2195–2195, 1999.

[30] A Fratalocchi, C Conti, and G Ruocco. Mode competitions and dynamical frequency
pulling in Mie nanolasers: 3d ab-initio Maxwell-Bloch computations. Optics express,
16(12):8342–8349, 2008.

[31] Matteo Frigo. A fast Fourier transform compiler. SIGPLAN Not., 34(5):169–180, May

1999.

[32] Amparo Gil, Javier Segura, and Nico Temme. Numerical Methods for Special Functions.

01 2007.

[33] Andreas Glaser and Vladimir Rokhlin. A new class of highly accurate solvers for

ordinary differential equations. Journal of Scientific Computing, 38(3):368–399, 2009.

[34] Connor Glosser. cglosser/QuEST: v0.1.0, May 2017.
[35] Connor Glosser. cglosser/QuEST: v0.2.1. May 2018.
[36] Connor Glosser, Carlo Piermarocchi, Jie Li, Dan Dault, and B. Shanker. Computational
dynamics of acoustically driven microsphere systems. Phys. Rev. E, 93:013305, Jan
2016.

[37] Leslie Frederick Greengard. The Rapid Evaluation of Potential Fields in Particle Systems.

PhD thesis, New Haven, CT, USA, 1987. AAI8727216.

[38] M. Gross and S. Haroche. Superradiance: An essay on the theory of collective

spontaneous emission. Physics Reports, 93(5):301 – 396, 1982.

[39] Gaël Guennebaud, Benoît Jacob, et al. Eigen v3. http://eigen.tuxfamily.org, 2010.
[40] Nail A. Gumerov and Ramani Duraiswami. Computation of scattering from n
spheres using multipole reexpansion. The Journal of the Acoustical Society of America,
112(6):2688–2701, 2002.

[41] Fritz Haake, Harald King, Guntram Schröder, Joe Haus, and Roy Glauber. Fluctuations

in superfluorescence. Phys. Rev. A, 20:2047–2063, Nov 1979.

[42] Sophie Hernot and Alexander L. Klibanov. Microbubbles in ultrasound-triggered
drug and gene delivery. Advanced Drug Delivery Reviews, 60(10):1153 – 1166, 2008.
Ultrasound in Drug and Gene Delivery.

[43] H. Htoon, T. Takagahara, D. Kulik, O. Baklenov, A. L. Holmes, and C. K. Shih. Interplay
of rabi oscillations and quantum interference in semiconductor quantum dots. Phys.
Rev. Lett., 88:087401, Feb 2002.

85

[44] Yurii A. Ilinskii, Mark F. Hamilton, and Evgenia A. Zabolotskaya. Bubble interaction
dynamics in lagrangian and hamiltonian mechanics. The Journal of the Acoustical Society
of America, 121(2):786–795, 2007.

[45] J.D. Jackson. Classical Electrodynamics. Wiley, 2007.
[46] Frank Jahnke, Christopher Gies, Marc Aßmann, Manfred Bayer, HAM Leymann,
Alexander Foerster, Jan Wiersig, Christian Schneider, Martin Kamp, and Sven Höfling.
Giant photon bunching, superradiant pulse emission and excitation trapping in
quantum-dot nanolasers. Nature communications, 7, 2016.

[47] Eric Jones, Travis Oliphant, Pearu Peterson, et al. SciPy: Open source scientific tools

for Python, 2001–. [Online; accessed 2014-10-14].

[48] H. Kamada, H. Gotoh, J. Temmyo, T. Takagahara, and H. Ando. Exciton rabi oscillation

in a single quantum dot. Phys. Rev. Lett., 87:246401, Nov 2001.

[49] S. Kapur and D. E. Long. IES3: a fast integral equation solver for efficient 3-dimensional
extraction. In 1997 Proceedings of IEEE International Conference on Computer Aided Design
(ICCAD), pages 448–455, Nov 1997.

[50] Louis V. King. On the acoustic radiation pressure on spheres. Proceedings of the Royal
Society of London A: Mathematical, Physical and Engineering Sciences, 147(861):212–240,
1934.

[51] Fredrik Kjolstad, Shoaib Kamil, Stephen Chou, David Lugato, and Saman Amarasinghe.
The tensor algebra compiler. Proc. ACM Program. Lang., 1(OOPSLA):77:1–77:29, October
2017.

[52] Donald Knuth. The Art of Computer Programming. Addison-Wesley, Reading, Mass,

1997.

[53] L.D. Landau and E.M. Lifshitz. Fluid Mechanics. Number v. 6. Elsevier Science, 2013.
[54] Jie Li, Daniel Dault, and Balasubramaniam Shanker. A quasianalytical time domain
solution for scattering from a homogeneous sphere. The Journal of the Acoustical Society
of America, 135(4):1676–1685, 2014.

[55] J. C. MacGillivray and M. S. Feld. Theory of superradiance in an extended, optically

thick medium. Phys. Rev. A, 14:1169–1189, Sep 1976.

[56] G. Manara, A. Monorchio, and R. Reggiannini. A space-time discretization criterion for
a stable time-marching solution of the electric field integral equation. IEEE Transactions
on Antennas and Propagation, 45(3):527–532, Mar 1997.

[57] G. Mur. Absorbing boundary conditions for the finite-difference approximation of
the time-domain electromagnetic-field equations. IEEE Transactions on Electromagnetic
Compatibility, EMC-23(4):377–382, Nov 1981.

86

[58] M. K. Myers and J. S. Hausmann. Computation of acoustic scattering from a moving

rigid surface. The Journal of the Acoustical Society of America, 91(5):2594–2605, 1992.

[59] Filip Petersson, Lena Åberg, Ann-Margret Swärd-Nilsson, and Thomas Laurell.
Free flow acoustophoresis: Microfluidic-based mode of particle and cell separation.
Analytical Chemistry, 79(14):5117–5123, 2007. PMID: 17569501.

[60] A. J. Pray, Y. Beghein, N. V. Nair, K. Cools, H. Bağcı, and B. Shanker. A higher order
space-time galerkin scheme for time domain integral equations. IEEE Transactions on
Antennas and Propagation, 62(12):6183–6191, Dec 2014.

[61] 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. IEEE Transactions on Antennas and Propagation, 60(8):3772–3781,
Aug 2012.

[62] William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery.
Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge University
Press, New York, NY, USA, 3 edition, 2007.

[63] Nicholas E. Rehler and Joseph H. Eberly. Superradiance. Phys. Rev. A, 3:1735–1751,

May 1971.

[64] E.J. Rothwell and M.J. Cloud. Electromagnetics, Second Edition. Electrical Engineering

Textbook Series. Taylor & Francis, 2009.

[65] Michael Scheibner, Thomas Schmidt, Lukas Worschech, Alfred Forchel, Gerd Bacher,
Thorsten Passow, and Detlef Hommel. Superradiance of quantum dots. Nat Phys,
3(2):106–110, Feb 2007.

[66] Tal Schwartz, Guy Bartal, Shmuel Fishman, and Mordechai Segev. Transport and ander-
son localization in disordered two-dimensional photonic lattices. Nature, 446(7131):52–
55, Mar 2007.

[67] B. Shanker, A. A. Ergin, K. Aygun, and E. Michielssen. Analysis of transient electro-
magnetic scattering from closed surfaces using a combined field integral equation.
IEEE Transactions on Antennas and Propagation, 48(7):1064–1074, Jul 2000.

[68] G. Ya. Slepyan, A. Magyarov, S. A. Maksimenko, A. Hoffmann, and D. Bimberg. Rabi
oscillations in a semiconductor quantum dot: Influence of local fields. Phys. Rev. B,
70:045320, Jul 2004.

[69] G. Ya. Slepyan, S. A. Maksimenko, A. Hoffmann, and D. Bimberg. Quantum optics of

a quantum dot: Local-field effects. Phys. Rev. A, 66:063804, Dec 2002.

[70] T. H. Stievater, Xiaoqin Li, D. G. Steel, D. Gammon, D. S. Katzer, D. Park, C. Piermaroc-
chi, and L. J. Sham. Rabi oscillations of excitons in single quantum dots. Phys. Rev.
Lett., 87:133603, Sep 2001.

87

[71] C. R. Stroud, J. H. Eberly, W. L. Lama, and L. Mandel. Superradiant effects in systems

of two-level atoms. Phys. Rev. A, 5:1094–1104, Mar 1972.

[72] Vasily V. Temnov and Ulrike Woggon. Superradiance and subradiance in an inhomo-
geneously broadened ensemble of two-level systems coupled to a low-q cavity. Phys.
Rev. Lett., 95:243602, Dec 2005.

[73] Vasily V Temnov and Ulrike Woggon. Photon statistics in the cooperative spontaneous

emission. Optics express, 17(7):5774–5782, 2009.

[74] M. L. W. Thewalt. Details of the structure of bound excitons and bound multiexciton

complexes in si. Canadian Journal of Physics, 55(17):1463–1480, 1977.

[75] A. Tiwari, C. Pantano, and J. B. Freund. Growth-and-collapse dynamics of small

bubble clusters near a wall. Journal of Fluid Mechanics, 775:1–23, 7 2015.

[76] Federico Toschi and Eberhard Bodenschatz. Lagrangian properties of particles in

turbulence. Annual Review of Fluid Mechanics, 41(1):375–404, 2009.

[77] L. Tsang, K. H. Ding, S. E. Shih, and J. A. Kong. Scattering of electromagnetic waves
from dense distributions of spheroidal particles based on monte carlo simulations. J.
Opt. Soc. Am. A, 15(10):2660–2669, Oct 1998.

[78] Xander F. van Kooten, Julien Autebert, and Govind V. Kaigala. Passive removal of
immiscible spacers from segmented flows in a microfluidic probe. Applied Physics
Letters, 106(7):074102, 2015.

[79] C. Vanneste and P. Sebbah. Selective excitation of localized modes in active random

media. Phys. Rev. Lett., 87:183903, Oct 2001.

[80] Pete Vukusic, Benny Hallam, and Joe Noyes. Brilliant whiteness in ultrathin beetle

scales. Science, 315(5810):348–348, 2007.

[81] P. C. Waterman. New formulation of acoustic scattering. The Journal of the Acoustical

Society of America, 45(6):1417–1429, 1969.

[82] Ling Ye, Jing Liu, Ping Sheng, and D.A. Weitz. Sound propagation in suspensions of

solid spheres. Phys. Rev. E, 48:2805–2815, Oct 1993.

[83] Zhen Ye. Low-frequency acoustic scattering by gas-filled prolate spheroids in liquids.

The Journal of the Acoustical Society of America, 101(4):1945–1952, 1997.

[84] A. E. Yilmaz, Jian-Ming Jin, and E. Michielssen. Time domain adaptive integral
method for surface integral equations. IEEE Transactions on Antennas and Propagation,
52(10):2692–2708, Oct 2004.

[85] Zorana Zeravcic, Detlef Loshe, and Wim Van Saarloos. Collective oscillations in bubble

clouds. Journal of Fluid Mechanics, 680:114–149, 8 2011.

88