Computational Photonics: Frequency and Time Domain Methods

1
Computational PhotonicsFrequency and Time
Domain Methods

• Steven G. Johnson
• MIT Applied Mathematics

2
Nano-photonic media (l-scale)
strange waveguides
microcavities
B. Norris, UMN
Assefa Kolodziejski, MIT
3d structures
Mangan, Corning
synthetic materials
optical phenomena
hollow-core fibers
3
Photonic Crystals
periodic electromagnetic media
can have a band gap optical insulators
4
Electronic and Photonic Crystals
atoms in diamond structure
Periodic Medium
Bloch waves Band Diagram
electron energy
wavevector
5
Electronic Photonic Modeling
Electronic
Photonic
strongly interacting entanglement, Coulomb
tricky approximations
non-interacting (or weakly), simple
approximations (finite resolution) any
desired accuracy
lengthscale dependent (from Plancks h)
scale-invariant e.g. size ?10 ? ? ?10
(except materials may change)
6
Computational Photonics Problems
Time-domain simulation start with current
J(x,t) run numerical experiment to simulate
E(x, t), H(x, t)
Frequency-domain linear response start with
harmonic current J(x, t) ei?t J(x) solve
for steady-state harmonic fields E(x), H(x)
involves solving linear equation Axb
Frequency-domain eigensolver solve for
source-free harmonic eigenfields E(x), H(x)
ei?t involves solving eigenequation Ax?2x
7
Numerical Methods Basis Choices
finite difference
finite elements
in irregular elements, approximate unknowns by
low-degree polynomial
discretize unknowns on regular grid
boundary-element methods
discretize only the boundaries between homogeneous
media solve integral equation via Greens
functions
spectral methods
complete basis of smooth functions (e.g. Fourier
series)


..
8
Numerical Methods Basis Choices
finite difference
finite elements
in irregular elements, approximate unknowns by
low-degree polynomial
discretize unknowns on regular grid
boundary-element methods
spectral methods
discretize only the boundaries between homogeneous
media solve integral equation via Greens
functions
complete basis of smooth functions (e.g. Fourier
series)


..
Much easier to analyze, implement, generalize,
parallelize, optimize,
Potentially much more efficient, especially
for high resolution
9
Computational Photonics Problems
Numerical Methods Basis Choices
finite difference
Time-domain simulation start with
current J(x,t) run numerical experiment
to simulate E(x, t), H(x, t)
spectral methods

Frequency-domain linear response start with
harmonic current J(x, t) ei?t J(x) solve
for steady-state harmonic fields E(x), H(x)
involves solving linear equation Axb

..
finite elements
in irregular elements, approximate unknowns by
low-degree polynomial
Frequency-domain eigensolver solve for
source-free harmonic eigenfields E(x), H(x)
ei?t involves solving eigenequation Ax?2x
boundary-element methods
discretize only the boundaries between homogeneous
media
10
FDTD Finite-Difference Time-Domain methods
Divide both space and time into discrete grids
spatial resolution ?x temporal resolution
?t Very general arbitrary geometries,
materials, nonlinearities, dispersion, sources,
any photonics calculation, in principle
dielectric function e(x) n2(x)
11
The Yee Discretization (1966)
a cubic voxel ?x ? ?y ? ?z
Staggered grid in space every
field component is stored on a different grid
12
The Yee Discretization (1966)
O(?x2) O(?y2)
all derivatives become center differences
13
The Yee Discretization (1966)
all derivatives become center differences includi
ng derivatives in time
O(?t2)
Explicit time-stepping stability requires
(vs. implicit time steps invert large matrix at
each step)
14
FDTD Discretization Upshot
For stability, space and time resolutions are
proportional doubling resolution in 3d
requires at least 16 24 times the work!
But at least the error goes quadratically with
resolution right? not necessarily!
15
Difficulty with a gridrepresenting
discontinuous materials?
staircasing
how does this affect accuracy?
16
Field Discontinuity DegradesOrder of Accuracy
E
TE polarization (E in plane discontinuous)
17
Sub-pixel smoothing
Can eliminate discontinuity by grayscaling
assign some average e to each pixel
e?
discretizing a smoothed structure that means
we are changing geometry can actually add to
error
18
Past sub-pixel smoothing methods can make error
worse!
convergence is still only linear
Three previous smoothing methods
Dey, 1999 Kaneda, 1997 Mohammadi, 2005
19
A Criterion for Accurate Smoothing
1st-order errors from smoothing De
We want the smoothing errors to be zero to 1st
order minimizes error and 2nd-order convergent!
Use a tensor e
(in principal axes)
Meade et al., 1993
20
Consistently the Lowest Error
quadratic!
quadratic accuracy!
Farjadpour et al., Opt. Lett. 2006
21
in 3d too
(notice that ranking of other methods has
shuffled!)
22
A qualitatively different case corners
still lowest error, but not quadratic
zero-perturbation criterion not satisfied due to
E divergence at corner analytically, error
?x1.404
23
Yes, but what can you do with FDTD?
Some common tasks
Frequency-domain response put in harmonic
source and wait for steady-state
Transmission/reflection spectra get entire
spectrum from a single simulation (Fourier
transform of impulse response)
Eigenmodes and resonant modes get all modes
from a single simulation (some tricky signal
processing)
24
Transmission Spectra in FDTD
? 1
? 12
a
example a 90 bend, 2d strip waveguide
transmitted power energy flux here
25
Transmission Spectra in FDTD
? 1
? 12
a
Gaussian-pulse current source J
26
Transmission Spectra in FDTD
must always do two simulations one for
normalization
P0(?)
electric field Ez
P(?)
transmission P(?) / P0(?)
27
Reflection Spectra in FDTD
? 1
? 12
for reflection, subtract incident fields (from
normalization run)
28
Transmission/Reflection Spectra
? 1
? 12
a
T
1TR
R
?a / 2pc a / ?
29
Dimensionless Units
Maxwells equations are scale invariant most
useful quantities are dimensionless ratios like
a / ?, for a characteristic lengthscale a same
ratio, same ?, ? same solution regardless of
whether a 1µm or 1km Our (typical) approach
pick characteristic lengthscale a measure
distance in units of a measure time in units
of a/c measure ? in units of 2pc/a a / ?
....
30
Absorbing BoundariesPerfectly Matched Layers
Artificial absorbing material overlapping the
computation Theoretically reflectionless but
PML is no longer perfect with finite resolution,
so gradually turn on absorption over
finite-thickness PML
perfect absorber PML
31
Computational Photonics Problems
Numerical Methods Basis Choices
finite difference
Time-domain simulation start with
current J(x,t) run numerical experiment
to simulate E(x, t), H(x, t)
spectral methods

Frequency-domain linear response start with
harmonic current J(x, t) ei?t J(x) solve
for steady-state harmonic fields E(x), H(x)
involves solving linear equation Axb

..
finite elements
in irregular elements, approximate unknowns by
low-degree polynomial
Frequency-domain eigensolver solve for
source-free harmonic eigenfields E(x), H(x)
ei?t involves solving eigenequation Ax?2x
boundary-element methods
discretize only the boundaries between homogeneous
media
32
A Maxwell Eigenproblem
First task get rid of this mess
0
dielectric function e(x) n2(x)
33
Electronic Photonic Eigenproblems
Electronic
Photonic
simple linear eigenproblem (for linear materials)
nonlinear eigenproblem (V depends on e density
?2)
many well-known computational techniques
Hermitian real E w, Periodicity Blochs
theorem
34
A 2d Model System
dielectric atom e12 (e.g. Si)
square lattice, period a
a
a
E
TM
H
35
Periodic Eigenproblems
if eigen-operator is periodic, then Bloch-Floquet
theorem applies
can choose
planewave
periodic envelope
Corollary 1 k is conserved, i.e. no scattering
of Bloch wave
Corollary 2 given by finite unit
cell, so w are discrete wn(k)
36
A More Familiar Eigenproblem
? 1
y
? 12
a
x
band diagram / dispersion relation
w
light cone (all non-guided modes)
find the normal modes of the waveguide
(propagation constant k a.k.a. ?)
k
37
Solving the Maxwell Eigenproblem
Finite cell ? discrete eigenvalues wn
Want to solve for wn(k), plot vs. all k for
all n,
constraint
where
H(x,y) ei(k?x wt)
Limit range of k irreducible Brillouin zone
1
Limit degrees of freedom expand H in finite basis
2
Efficiently solve eigenproblem iterative methods
3
38
Solving the Maxwell Eigenproblem 1
Limit range of k irreducible Brillouin zone
1
Blochs theorem solutions are periodic in k
M
first Brillouin zone minimum k primitive
cell
X
G
irreducible Brillouin zone reduced by symmetry
Limit degrees of freedom expand H in finite basis
2
Efficiently solve eigenproblem iterative methods
3
39
Solving the Maxwell Eigenproblem 2a
Limit range of k irreducible Brillouin zone
1
Limit degrees of freedom expand H in finite
basis (N)
2
solve
finite matrix problem
Efficiently solve eigenproblem iterative methods
3
40
Solving the Maxwell Eigenproblem 2b
Limit range of k irreducible Brillouin zone
1
Limit degrees of freedom expand H in finite basis
2
must satisfy constraint
Planewave (FFT) basis
Finite-element basis
constraint, boundary conditions
Nédélec elements
Nédélec, Numerische Math. 35, 315 (1980)
constraint
nonuniform mesh, more arbitrary
boundaries, complex code mesh, O(N)
uniform grid, periodic boundaries, simple code,
O(N log N)
figure Peyrilloux et al., J. Lightwave
Tech. 21, 536 (2003)
Efficiently solve eigenproblem iterative methods
3
41
Solving the Maxwell Eigenproblem 3a
Limit range of k irreducible Brillouin zone
1
Limit degrees of freedom expand H in finite basis
2
Efficiently solve eigenproblem iterative methods
3
Slow way compute A B, ask LAPACK for
eigenvalues requires O(N2) storage, O(N3) time
Faster way start with initial guess
eigenvector h0 iteratively improve O(Np)
storage, O(Np2) time for p eigenvectors
(p smallest eigenvalues)
42
Solving the Maxwell Eigenproblem 3b
Limit range of k irreducible Brillouin zone
1
Limit degrees of freedom expand H in finite basis
2
Efficiently solve eigenproblem iterative methods
3
Many iterative methods Arnoldi, Lanczos,
Davidson, Jacobi-Davidson, ,
Rayleigh-quotient minimization
43
Solving the Maxwell Eigenproblem 3c
Limit range of k irreducible Brillouin zone
1
Limit degrees of freedom expand H in finite basis
2
Efficiently solve eigenproblem iterative methods
3
Many iterative methods Arnoldi, Lanczos,
Davidson, Jacobi-Davidson, ,
Rayleigh-quotient minimization
for Hermitian matrices, smallest eigenvalue w0
minimizes
minimize by preconditioned conjugate-gradient
(or)
variational theorem
44
Band Diagram of 2d Model System(radius 0.2a
rods, e12)
a
frequency w (2pc/a) a / l
G
X
M
G
irreducible Brillouin zone
M
E
gap for n gt 1.751
TM
X
G
H
45
Origin of Gap in 2d Model System
orthogonal node in high e
Ez
lives in high e
G
X
M
G
Ez
E
gap for n gt 1.751
TM
H
46
The Iteration Scheme is Important
(minimizing function of 104108 variables!)
Steepest-descent minimize (h a ?f) over a
repeat
Conjugate-gradient minimize (h a d) d is
?f (stuff) conjugate to previous search dirs
Preconditioned steepest descent minimize (h a
d) d (approximate A-1) ?f Newtons
method
47
The Iteration Scheme is Important
(minimizing function of 40,000 variables)
no preconditioning
error
preconditioned conjugate-gradient
no conjugate-gradient
iterations
48
The Boundary Conditions are Tricky
e?
49
The e-averaging is Important
backwards averaging
correct averaging changes order of
convergence from ?x to ?x2
no averaging
error
tensor averaging
(similar effects in other EM numerics
analyses)
resolution (pixels/period)
50
Gap, Schmap?
a
frequency w
G
X
M
G
But, what can we do with the gap?
51
Intentional defects are good
microcavities
waveguides (wires)
52
Intentional defects in 2d
53
Microcavity Blues
For cavities (point defects) frequency-domain has
its drawbacks
Best methods compute lowest-w bands, but Nd
supercells have Nd modes below the cavity mode
expensive
Best methods are for Hermitian operators,
but losses requires non-Hermitian
54
Time-Domain Eigensolvers(finite-difference
time-domain FDTD)
Simulate Maxwells equations on a discrete
grid, absorbing boundaries (leakage loss)
Excite with broad-spectrum dipole ( ) source
Dw
Response is many sharp peaks, one peak per mode
signal processing
complex wn
Mandelshtam, J. Chem. Phys. 107, 6756 (1997)
decay rate in time gives loss
55
Signal Processing is Tricky
signal processing
complex wn
?
a common approach least-squares fit of spectrum
fit to
FFT
Decaying signal (t)
Lorentzian peak (w)
56
Fits and Uncertainty
problem have to run long enough to completely
decay
actual
signal portion
Portion of decaying signal (t)
Unresolved Lorentzian peak (w)
57
Unreliability of Fitting Process
Resolving two overlapping peaks
is near-impossible 6-parameter nonlinear fit (too
many local minima to converge reliably)
sum of two peaks
w 10.033i
w 1.030.025i
Sum of two Lorentzian peaks (w)
58
Quantum-inspired signal processing (NMR
spectroscopy)Filter-Diagonalization Method (FDM)
Mandelshtam, J. Chem. Phys. 107, 6756 (1997)
Given time series yn, write
find complex amplitudes ak frequencies wk by a
simple linear-algebra problem!
59
Filter-Diagonalization Method (FDM)
Mandelshtam, J. Chem. Phys. 107, 6756 (1997)
We want to diagonalize U eigenvalues of U are
eiw?t
60
Filter-Diagonalization Summary
Mandelshtam, J. Chem. Phys. 107, 6756 (1997)
Umn given by yns just diagonalize known matrix!
A few omitted steps Generalized eigenvalue
problem (basis not orthogonal) Filter yns
(Fourier transform) small bandwidth smaller
matrix (less singular)
resolves many peaks at once peaks not
known a priori resolve overlapping peaks
resolution gtgt Fourier uncertainty
61
Do try this at home
FDTD simulation http//ab-initio.mit.edu/meep/
Bloch-mode eigensolver http//ab-initio.mit.edu
/mpb/
Filter-diagonalization http//ab-initio.mit.edu
/harminv/
Photonic-crystal tutorials ( THIS TALK)
http//ab-initio.mit.edu/ /photons/tutorial/
62
Meep (FDTD) MPB (Eigensolver)
Arbitrary??(x) including dispersive,
loss/gain, and nonlinear ?(2) and ?(3)
Arbitrary periodic??(x) anisotropic,
magneto-optic, (lossless, linear materials)
Arbitrary?J(x,t)
1d/2d/3d
PML/periodic/metal bound.
band diagrams, group velocities perturbation
theory,
1d/2d/3d/cylindrical
power spectra eigenmodes
Free/open-source software (GNU)
fully scriptable interface
built-in multivariate optimization, integration,
root-finding,
MPI parallelism
exploit mirror symmetries
field output (standard HDF5 format)
63
Unix Philosophycombine small, well-designed
tools, via files
Input text file
MPB/Meep
standard formats (text HDF5)
Disadvantage have to learn several
programs Advantages flexibility batch
processing, shell scripting ease of development
Visualization / Analysis software (Matlab, Mayavi
vtk, command-line tools, )
64
Unix Philosophycombine small, well-designed
tools, via files
Input text file
MPB/Meep
standard formats (text HDF5)
GNU Guile scripting interpreter (Scheme language)
Visualization / Analysis software (Matlab, Mayavi
vtk, command-line tools, )
Embed a full scripting language parameter
sweeps complex parameterized geometries
optimization, integration, etc. programmable
J(x, t), etc. Turing complete
65
A Simple Example (MPB)
? 1
y
? 12
a
x
find the normal modes ?n(k) of the waveguide
Need to specify computational cell
size/resolution geometry, i.e. ?(y) what
k values how many modes (n 1, 2, ?)
66
A File Format Made of Parentheses
Need to specify computational cell
size/resolution (set! geometry-lattice (make
lattice (size no-size 10 no-size) (set!
resolution 32) geometry, i.e. ?(y) what
k values how many modes (n 1, 2, ?)
1 pixel
? 1
y
10 (320 pixels)
? 12
a
x
67
A File Format Made of Parentheses
Need to specify computational cell
size/resolution geometry, i.e. ?(y) (set!
geometry (list (make block (size infinity 1
infinity) (center 0 0 0)
(material (make dielectric (epsilon
12)))))) what k values how many modes (n
1, 2, ?)
(choose units of a)
? 1
y
? 12
a
x
68
A File Format Made of Parentheses
Need to specify computational cell
size/resolution geometry, i.e. ?(y) what
k values (set! k-points (interpolate 10 (list
(vector3 0 0 0) (vector3 2 0 0)))) how many
modes (n 1, 2, ?)
(units of 2p/a)
(built-in function)
? 1
y
? 12
a
x
69
A File Format Made of Parentheses
Need to specify computational cell
size/resolution geometry, i.e. ?(y) what
k values how many modes (n 1, 2, ?)
(set! num-bands 5)
Then run (run) or only TM polarization (run-tm
) or only TM, even modes (run-tm-yeven)
? 1
y
? 12
a
x
70
Simple Example (MPB) Results
? 1
y
? 12
a
x
find the normal modes ?n(k) of the waveguide
red even blue odd
71
Do try this at home
FDTD simulation http//ab-initio.mit.edu/meep/
Bloch-mode eigensolver http//ab-initio.mit.edu
/mpb/
Filter-diagonalization http//ab-initio.mit.edu
/harminv/
Photonic-crystal tutorials ( THIS TALK)
http//ab-initio.mit.edu/ /photons/tutorial/