Transcript Computational Photonics: Frequency and Time Domain Methods

Computational Photonics:
Frequency and Time Domain Methods
Steven G. Johnson
MIT Applied Mathematics
Nano-photonic media (l-scale)
strange waveguides
& microcavities
[B. Norris, UMN]
[Assefa & Kolodziejski,
MIT]
3d
structures
[Mangan,
Corning]
synthetic materials
hollow-core fibers
optical phenomena
Photonic Crystals
periodic electromagnetic media
1887
1987
1-D
2-D
periodic in
one direction
periodic in
two directions
3-D
periodic in
three directions
can have a band gap: optical “insulators”
Electronic and Photonic Crystals
dielectric spheres, diamond lattice
photon frequency
electron energy
Bloch waves:
Band Diagram
Periodic
Medium
atoms in diamond structure
wavevector
wavevector
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 Planck’s h)
• scale-invariant
—e.g. size 10  l 10
(except materials may change)
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) = e–it J(x)
— solve for steady-state harmonic fields E(x), H(x)
— involves solving linear equation Ax=b
• Frequency-domain eigensolver
— solve for source-free harmonic eigenfields
E(x), H(x) ~ e–it
— involves solving eigenequation Ax=2x
Numerical Methods: Basis Choices
finite difference
finite elements
in irregular “elements,”
approximate unknowns
by low-degree polynomial
discretize
unknowns
on regular grid
df f (x  x)  f (x  x)

 O(x 2 )
dx
x
boundary-element methods
spectral methods
complete basis of
smooth functions
(e.g. Fourier series)
+
+
+ …..
discretize only the
boundaries between
homogeneous media
…solve
integral equation
via Green’s functions
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
complete basis of
smooth functions
(e.g. Fourier series)
+
+
+ …..
Much easier to analyze, implement,
generalize, parallelize, optimize, …
discretize only the
boundaries between
homogeneous media
…solve
integral equation
via Green’s functions
Potentially much more efficient,
especially for high resolution
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)
• Frequency-domain linear response
— start with harmonic current J(x, t) = e–it J(x)

— solve for steady-state harmonic fields E(x), H(x)
— involves solving linear equation Ax=b
spectral methods
df f (x  x)  f (x  x)

 O(x 2 )
dx
x
+
+
+ …..
finite elements
• Frequency-domain eigensolver
— solve for source-free harmonic eigenfields
E(x), H(x) ~ e–it
— involves solving eigenequation Ax=2x
in irregular “elements,”
approximate unknowns
by low-degree polynomial
boundary-element methods
discretize only the
boundaries between
homogeneous media
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
H
1
   E
t

E 1
J
  H 
t 

dielectric function (x) = n2(x)
The Yee Discretization (1966)
a cubic “voxel”: ∆x  ∆y  ∆z
(i, j+1)
(i+1, j+1, k+1)
Hz
Ey
(i, j, k+1)
Ez
Hx
Hy
(i, j)
(i+1, j+1, k)
Hz
Ex
(i+1, j)
Ey
(i, j, k)
Ex
(i+1, j, k)
Staggered grid in space:
— every field component is stored on a different grid
The Yee Discretization (1966)
H
1
   E 
t

H z
t
1
1
i , j 
2
2
1 Ey Ex 
  


  x
y 
(i, j+1)
Ey
(i, j)
Hz
Ex
(i+1, j)


1
1
1
1
E (i  1, j  )  Ey (i, j  ) Ex (i  , j  1)  Ex (i  , j) 
1  y
2
2 
2
2
  

 
x
y



+ O(∆x2) + O(∆y2)
all derivatives become center differences…
The Yee Discretization (1966)
all derivatives become center differences…
including derivatives in time
H
t

t nt
1

 E
t nt
1
1
H(n  )  H(n  )
2
2

t
+ O(∆t2)
Explicit time-stepping:
stability requires t 
x
# dimensions
(vs. implicit time steps: invert large matrix at each step)
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!
Difficulty with a grid:
representing discontinuous materials?
“staircasing”
… how does this affect accuracy?
Field Discontinuity Degrades
Order of Accuracy
E
TE polarization (E in plane: discontinuous)
a
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
a
Sub-pixel smoothing
Can eliminate
discontinuity
by “grayscaling”
— assign some average 
to each pixel
?
= discretizing a smoothed structure
— that means we are changing geometry
— can actually add to error
Past sub-pixel smoothing methods
can make error worse!
Three previous smoothing methods
& convergence is
still only linear
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
[ Dey, 1999 ]
[ Kaneda, 1997 ]
[ Mohammadi, 2005]
A Criterion for Accurate Smoothing
1st-order errors
from
smoothing 

1
2
2 
~   E||  ( ) D 



We want the smoothing errors to be zero to 1st order
— minimizes error and 2nd-order convergent!

Use a tensor : 


(in principal axes:)
[ Meade et al., 1993 ]




1

 E||

1 
 E
Consistently the Lowest Error
a
quadratic!
quadratic accuracy!
a
[ Farjadpour et al., Opt. Lett. 2006 ]
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
A qualitatively different case: corners
still ~lowest error, but not quadratic
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
zero-perturbation
criterion
not satisfied
due to E divergence
at corner
— analytically,
error ~ ∆x1.404
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)
Transmission Spectra in FDTD
=1
a
 = 12
example: a 90° bend, 2d strip waveguide
transmitted power = energy flux here:
Transmission Spectra in FDTD
=1
a
 = 12
Gaussian-pulse
current source J
Fourier-transform the fields at each x:
E 
 E(t)e
it
dt   E(nt)e
int
n
1
P() 
2
 Re[E
*
 H  ]dx
t
Transmission Spectra in FDTD
must always do two simulations: one for normalization
P0()
electric
field Ez:
QuickTime™ and a
Y UV420 codec decompressor
are needed to see this picture.
P()
QuickTime™ and a
Y UV420 codec decompressor
are needed to see this picture.
transmission = P() / P0()
Reflection Spectra in FDTD
=1
 = 12
1
PR () 
2

Re[(E  E0 )*  (H   H 0 )]dx
for reflection, subtract incident fields
(from normalization run)
1
P() 
2
 Re[E
*
 H  ]dx
Transmission/Reflection Spectra
Grap
are
QuickTim
needed
hics e™
deco
toand
see
mp res
athis
sor
picture .
1
0.9
=1
0.8
a
 = 12
0.7
0.6
0.5
T
1–T–R
0.4
0.3
0.2
0.1
0
0.1
R
0.11
0.12
0.13
0.14
0.15
0.16
0.17
a / 2πc = a / l
0.18
0.19
0.2
Dimensionless Units
Maxwell’s equations are scale invariant
— most useful quantities are dimensionless ratios
like a / l, 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 2πc/a = a / l
– ....
Absorbing Boundaries:
Perfectly Matched Layers
“perfect” absorber: PML
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
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)
• Frequency-domain linear response
— start with harmonic current J(x, t) = e–it J(x)

— solve for steady-state harmonic fields E(x), H(x)
— involves solving linear equation Ax=b
spectral methods
df f (x  x)  f (x  x)

 O(x 2 )
dx
x
+
+
+ …..
finite elements
• Frequency-domain eigensolver
— solve for source-free harmonic eigenfields
E(x), H(x) ~ e–it
— involves solving eigenequation Ax=2x
in irregular “elements,”
approximate unknowns
by low-degree polynomial
boundary-element methods
discretize only the
boundaries between
homogeneous media
A Maxwell Eigenproblem
1

 E  
H i H
c t
c
0
1

 H  
E  J  i E
c t
c
First task:
get rid of this mess
dielectric function (x) = n2(x)
 
    H    H
 c 

1
eigen-operator
2
eigen-value
+ constraint
 H  0
eigen-state
Electronic & Photonic Eigenproblems
Electronic
Photonic
2
 2 2



1

  V   E     H    H

 c 

 2m

nonlinear eigenproblem
(V depends on e density ||2)
simple linear eigenproblem
(for linear materials)
—many well-known
computational techniques
Hermitian = real E & , … Periodicity = Bloch’s theorem…
A 2d Model System
dielectric “atom”
=12 (e.g. Si)
square lattice,
period a
a
a
TM
E
H
Periodic Eigenproblems
if eigen-operator is periodic, then Bloch-Floquet theorem applies:
can choose:

i k x t
H (x ,t)  e
planewave

H k (x )
periodic “envelope”
Corollary 1: k is conserved, i.e. no scattering of Bloch wave
Corollary 2: Hk given by finite unit cell,
so  are discrete n(k)
A More Familiar Eigenproblem
=1
y
 = 12
a
x

find the normal modes
of the waveguide:
H(y,t)  H k (y)e
band diagram / dispersion relation
light cone
(all non-guided modes)
i(kxt )
(propagation constant k
a.k.a. )
k
Solving the Maxwell Eigenproblem
Finite cell  discrete eigenvalues n
1
  ik    ik  H n 

Want to solve for n(k),
& plot vs. “all” k for “all” n,
constraint:
Grap
are
QuickTim
neede
hics deco
e™
d toand
see
mp res
athis
sor
picture .

1
0.9
0.8
where:
0.7
n 2
c
  ik H  0
H(x,y) ei(kx – t)
0.6
0.5
0.4
Ph otonic Ba nd Ga p
0.3
0.2
TM b ands
0.1

0
1
Limit range of k: irreducible Brillouin zone
2
Limit degrees of freedom: expand H in finite basis
3
Efficiently solve eigenproblem: iterative methods
2
Hn
Solving the Maxwell Eigenproblem: 1
1
Limit range of k: irreducible Brillouin zone
—Bloch’s theorem: solutions are periodic in k
M
first Brillouin zone
= minimum |k| “primitive cell”
G
X
2
a
ky
kx
 reduced by symmetry
irreducible Brillouin zone:
2
Limit degrees of freedom: expand H in finite basis
3
Efficiently solve eigenproblem: iterative methods
Solving the Maxwell Eigenproblem: 2a
1
Limit range of k: irreducible Brillouin zone
2
Limit degrees of freedom: expand H in finite basis (N)
N
H  H(xt )   hm bm (x t )
2
ˆ
solve: A H   H
m1
finite matrix problem:
f g   f* g
3
Ah   Bh
Am  bm Aˆ b
2
Bm  bm b
Efficiently solve eigenproblem: iterative methods
Solving the Maxwell Eigenproblem: 2b
1
Limit range of k: irreducible Brillouin zone
2
Limit degrees of freedom: expand H in finite basis
— must satisfy constraint: ( ik) H 0
Planewave (FFT) basis
H(x t )   HG e
iGxt
Finite-element basis
constraint, boundary conditions:

Nédélec elements
G
constraint:
HG  G  k  0
uniform “grid,” periodic boundaries,
simple code, O(N log N)

[ Nédélec, Numerische Math.
35, 315 (1980) ]
3
[ figure: Peyrilloux et al.,
J. Lightwave Tech.
21, 536 (2003) ]
nonuniform mesh,
more arbitrary boundaries,
complex code & mesh, O(N)
Efficiently solve eigenproblem: iterative methods
Solving the Maxwell Eigenproblem: 3a
1
Limit range of k: irreducible Brillouin zone
2
Limit degrees of freedom: expand H in finite basis
3
Efficiently solve eigenproblem: iterative methods
Ah   Bh
2
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)
Solving the Maxwell Eigenproblem: 3b
1
Limit range of k: irreducible Brillouin zone
2
Limit degrees of freedom: expand H in finite basis
3
Efficiently solve eigenproblem: iterative methods
Ah   Bh
2
Many iterative methods:
— Arnoldi, Lanczos, Davidson, Jacobi-Davidson, …,
Rayleigh-quotient minimization
Solving the Maxwell Eigenproblem: 3c
1
Limit range of k: irreducible Brillouin zone
2
Limit degrees of freedom: expand H in finite basis
3
Efficiently solve eigenproblem: iterative methods
Ah   Bh
2
Many iterative methods:
— Arnoldi, Lanczos, Davidson, Jacobi-Davidson, …,
Rayleigh-quotient minimization
for Hermitian matrices, smallest eigenvalue 0 minimizes:
“variational
theorem”
h' Ah
  min
h h' Bh
2
0
minimize by preconditioned
conjugate-gradient (or…)
Band Diagram
of
2ddeco
Model
System
Grap
are
QuickTim
neede
hics
e™
d to
and
see
mp res
athis
sor
picture
.
(radius 0.2a rods, =12)
a
frequency  (2πc/a) = a / l
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
Ph otonic Ba nd Ga p
0.2
TM b ands
0.1
0
irreducible Brillouin zone
M
k
G
X
G
TM
X
E
H
M
G
gap for
n > ~1.75:1
The Iteration Scheme is Important
(minimizing function of 104–108+ variables!)
h' Ah
  min
 f (h)
h h'Bh
2
0
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 ~ Newton’s method
Preconditioned conjugate-gradient: minimize (h + ad)
— d is (approximate A-1) [f + (stuff)]
The Iteration Scheme is Important
Grap
are
QuickTim
neede
hics deco
e™
d toand
see
mp re
athis
ssor
p icture .
(minimizing function of ~40,000 variables)
1000000
100000
10000
% error
1000JÑ
100
10
1
0.1
0.01
0.001
0.0001
0.00001
0.000001
1
E
E
E EE
E EEE
EEEEE
no preconditioning
JÑ
E
E
E
E
E
E
JÑ Ñ
EEEEEEEEE
J Ñ
EEEEEEEEEEEEEEEEEEEEEE
Ñ
J Ñ
EEEE
EEEEEEEEE
J ÑÑ
EEEEEEEE
J ÑÑÑ
EEEEEEE
J J ÑÑÑ
EEEEEEE
JJ JJJ ÑÑÑÑ
Ñ
Ñ
Ñ
EEEEEEEE
Ñ
Ñ
Ñ
Ñ
Ñ
JJ ÑÑ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
EEEEEEEEE
Ñ
Ñ
Ñ
Ñ
JJ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
EEEEEE
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
J
Ñ
Ñ
Ñ
EEEEEE
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
J
Ñ
EEEEE
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
EEEEE
Ñ
Ñ
Ñ
J
Ñ
Ñ
Ñ
EEEE
Ñ
Ñ
Ñ
Ñ
Ñ
EEEE
J
Ñ
Ñ
Ñ
Ñ
EEEE
Ñ
Ñ
Ñ
Ñ
EEE
preconditioned JJ
Ñ
Ñ
EEEE
EEEE
no conjugate-gradient
EEEE
conjugate-gradient JJ
J
J
10
100
# iterations
1000
The Boundary Conditions are Tricky
E|| is continuous
E is discontinuous
(D = E is continuous)
?
Any single scalar  fails:
(mean D) ≠ (any ) (mean E)
Use a tensor :
 





1

 E||

1 
 E
The -averaging is Important
Grap
are
QuickTim
needed
hics deco
e™toand
see
mp res
athis
sor
picture .
100
H
10
H
% error
J
J
1
backwards averaging
H H H
B
J
H
B
J B
J
B
B
J
B
B
B
J
tensor averaging
0.1
correct averaging
changes order
B no Baveraging
of convergence
B
B
B
from ∆x to ∆x2
J
H H H
J
H
H H
H
J
J
J J
0.01
10
resolution (pixels/period)
100
(similar effects
in other E&M
numerics & analyses)
Gap, Schmap?
Grap
are
QuickTim
neede
hics deco
e™
d toand
see
mp res
athis
sor
picture .
1
a
0.9
frequency 
0.8
0.7
0.6
0.5
0.4
0.3
Ph otonic Ba nd Ga p
0.2
TM b ands
0.1
0
G
X
M
But, what can we do with the gap?
G
Intentional “defects” are good
microcavities
3D Pho to nic C rysta l with De fe c ts
waveguides (“wires”)
Intentional “defects” in 2d
are
G
Q rap
uickTim
needed
hics decom
e™toase
ndpressor
eathis picture.
(Same computation, with supercell = many primitive cells)
Grap
are
QuickTim
hics
e™
deco
tossor
and
see
mp
res
athis
sor
picture .
are
Gra
QuickTime™
phics
neede
dec
dneeded
toompre
and
see
athis
p icture.
a
Microcavity Blues
For cavities (point defects)
frequency-domain has its drawbacks:
• Best methods compute lowest- bands,
but Nd supercells have Nd modes
below the cavity mode — expensive
• Best methods are for Hermitian operators,
but losses requires non-Hermitian
Time-Domain Eigensolvers
(finite-difference time-domain = FDTD)
Simulate Maxwell’s equations on a discrete grid,
+ absorbing boundaries (leakage loss)
• Excite with broad-spectrum dipole ( ) source

signal processing
complex n
[ Mandelshtam,
J. Chem. Phys. 107, 6756 (1997) ]
decay rate in time gives loss
Response is many
sharp peaks,
one peak per mode
Signal Processing is Tricky
signal processing
complex n
?
Gra
are
QuickTim
phics
neededec
e™
d toom
and
see
pre
athis
ssor
p icture.
a common approach: least-squares
fit
of
spectrum
Gra
are
QuickTime™
phics
neededec
d toompre
and
seeathis
ssor
p icture.
0.8
450
0.6 EEE
0.4 EE
E
0.2 E E EEEE
E
EE
E E
EEEEE EEEEE EEEE EEEE
0 EEEE
E
EEEEEE EEEEEE EEEEEEE EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
E
EEE
E
E
E EEE
-0.2 E EE
E
EE
E
E
-0.4
E
-0.6
E
-0.8
E
-1 E
0 1 2 3 4 5 6 7 8 9 10
400
Decaying signal (t)
E
fit to:
350
300
EE
250
FFT
200
150
E
A
(   0 ) 2  G 2
E
100
E
E
E
50 EE
EE
EEEEE
EEE
E
EEEEEEEEEEEEEEEEEEE
0
0 0.5 1 1.5 2 2.5 3 3.5 4
Lorentzian peak ()
Fits and Uncertainty
Gra
are
QuickTime™
phics
needede
to
compressor
and
seeathis picture.
problem: have to run long enough
todcompletely
decay
are
Gra
QuickTim
phics
needee™
de
d to
com
and
see
pre
athis
ssor
picture.
400 00
1
EE
E
E
0.8 EE EE
EE
0.6 EE
EE
0.4
EE E E
0.2
EE
E
EE
EE EE
E
EE EE
E E EE
0 EEEEE
-0.2
E EE EE
-0.4
E EE
E
E
-0.6
E
E E
E
E
E
-0.8E EE E
E
-1 EE
0 1 2
EE
EE
EE
EE
E
EE
3
EE E
EE EE
EE E
E E EE
EE EE EEE EEE EEEE EEE
E
E EE
E E E E E E E E EE EE
E E EE
EEEEEEEEEEEEE
E EE E
E E E E EE EE E
E EE
EE EEE
E
E
E
E
E EE E
EE EE EE EE E
E
E
E
EE
EE EEEE EE
E
EE
4
5
6
7
8
9
Portion of decaying signal (t)
10
350 00
actual
300 00
250 00
200 00
E
150 00
100 00
500 0
signal
portion
E E E E E
0E E E E E
0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5
Unresolved Lorentzian peak ()
There is a better way, which gets complex  to > 10 digits
Unreliability of Fitting Process
Resolving two overlapping peaks is
Grap
are
QuickTim
needed
hics e™
deco
toand
se
mpresso
eathis picture.
r
near-impossible
6-parameter nonlinear fit
(too many local minima to converge reliably)
120 0
EE
100 0
800
E
sum of two peaks
E
600
400
 = 1+0.033i
E
E
200
E
E
EE
E
E
E
E
0 EEEEEEEEEEEEEE
0.5 0.6 0.7 0.8 0.9
E
1
E  = 1.03+0.025i
E
EE
EEEEEE
EEEEEEEEEEEE
1.1 1.2 1.3 1.4 1.5
Sum of two Lorentzian peaks ()
There is a better
way, which gets
complex 
for both peaks
to > 10 digits
Quantum-inspired signal processing (NMR spectroscopy):
Filter-Diagonalization Method (FDM)
[ Mandelshtam, J. Chem. Phys. 107, 6756 (1997) ]
Given time series yn, write:
y n  y(nt)   ake
i k nt
k
…find complex amplitudes ak & frequencies k
by a simple linear-algebra problem!
Idea: pretend y(t) is autocorrelation of a quantum system:


ˆ
H i

t
say:
time-∆t evolution-operator:
iHˆ t /
ˆ
U e
yn  (0) (nt)  (0) Uˆ n (0)
Filter-Diagonalization
Method
(FDM)
[ Mandelshtam, J. Chem. Phys. 107, 6756 (1997) ]
yn  (0) (nt)  (0) Uˆ n (0)
iHˆ t /
ˆ
U e
We want to diagonalize U: eigenvalues of U are ei∆t

…expand U in basis of |(n∆t)>:
Um,n  (mt) Uˆ (nt)  (0) Uˆ mUˆ Uˆ n (0)  ymn 1
Umn given by yn’s — just diagonalize known matrix!
Filter-Diagonalization
Summary
[ Mandelshtam, J. Chem. Phys. 107, 6756 (1997) ]
Umn given by yn’s — just diagonalize known matrix!
A few omitted steps:
—Generalized eigenvalue problem (basis not orthogonal)
—Filter yn’s (Fourier transform):
small bandwidth = smaller matrix (less singular)
• resolves many peaks at once
• # peaks not known a priori
• resolve overlapping peaks
• resolution >> Fourier uncertainty
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/
Meep (FDTD) MPB (Eigensolver)
• Arbitrary(x) — including
dispersive, loss/gain,
and nonlinear [(2) and (3)]
• Arbitrary J(x,t)
• PML/periodic/metal bound.
• 1d/2d/3d/cylindrical
• power spectra • eigenmodes
Free/open-source
software (GNU)
• MPI parallelism
• exploit mirror symmetries
• Arbitrary periodic(x) —
anisotropic, magneto-optic, …
(lossless, linear materials)
• 1d/2d/3d
• band diagrams, group velocities
perturbation theory, …
• fully scriptable interface
• built-in multivariate optimization,
integration, root-finding, …
• field output (standard HDF5 format)
Unix Philosophy
combine 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, …)
Unix Philosophy
combine small, well-designed tools, via files
Input text file
MPB/Meep
GNU Guile scripting interpreter
(Scheme language)
Embed a full scripting language:
— parameter sweeps
— complex parameterized geometries
— optimization, integration, etc.
— programmable J(x, t), etc.
— … Turing complete
standard formats
(text + HDF5)
Visualization / Analysis
software
(Matlab, Mayavi [vtk],
command-line tools, …)
A Simple Example (MPB)
=1
 = 12
a
find the normal modes n(k)
of the waveguide:
i(kxt )
H(y,t)  H k (y)e
Need to specify:

• computational cell size/resolution
• geometry, i.e. (y)
• what k values
• how many modes (n = 1, 2, … ?)
y
x
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, … ?)
a
 = 12
10 (320 pixels)
=1
1 pixel
y
x
A File Format Made of Parentheses
Need to specify:
• computational cell size/resolution
• geometry, i.e. (y)
(choose units of a)
(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, … ?)
=1
a
 = 12
y
x
A File Format Made of Parentheses
Need to specify:
• computational cell size/resolution
• geometry, i.e. (y)
(units of 2π/a)
• what k values
(set! k-points
(interpolate 10 (list (vector3 0 0 0) (vector3 2 0 0))))
(built-in function)
• how many modes (n = 1, 2, … ?)
=1
a
 = 12
y
x
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
a
 = 12
y
x
Simple Example (MPB) Results
=1
 = 12
a
find the normal modes n(k)
of the waveguide:
red = even
blue = odd
y
x
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/