# compute-methods_cornell by liaoxiuli

VIEWS: 17 PAGES: 69

• pg 1
```									  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
optical phenomena
hollow-core fibers
Photonic Crystals
periodic electromagnetic media
1887                     1987
1-D                    2-D             3-D

periodic in            periodic in       periodic in
one direction         two directions   three directions

can have a band gap: optical ―insulators‖
Electronic and Photonic Crystals
atoms in diamond structure      dielectric spheres, diamond lattice
Medium
Periodic

photon frequency
Band Diagram
Bloch waves:

electron energy

wavevector                        wavevector
Electronic & Photonic Modeling
Electronic                  Photonic
• strongly interacting     • non-interacting (or weakly),
—entanglement, Coulomb     —simple approximations
—tricky approximations       (finite resolution)
—any desired accuracy

• lengthscale dependent    • scale-invariant
(from Planck’s h)        —e.g. size 10  l 10
(except materials may change)
Computational Photonics Problems
• Time-domain simulation
— 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,‖
discretize                                                    approximate unknowns
unknowns                                                      by low-degree polynomial
on regular grid

df   f ( x  x)  f ( x  x)
                            O(x 2 )
dx              x

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

boundary-element methods
spectral methods
discretize only the
boundaries between
complete basis of                     homogeneous media
+                     smooth functions
(e.g. Fourier series)
+                                                               …solve
integral equation
+ …..
via Green’s functions

Much easier to analyze, implement,                Potentially much more efficient,
generalize, parallelize, optimize, …               especially for high resolution
Computational Photonics Problems                          Numerical Methods: Basis Choices

finite difference
• Time-domain simulation
— run ―numerical experiment‖ to                                                                      spectral methods
simulate E(x, t), H(x, t)
df   f ( x  x)  f ( x  x)
dx

x
 O(x 2 )      +
• 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
finite elements
• Frequency-domain eigensolver                                                                 in irregular ―elements,‖
approximate unknowns
— solve for source-free harmonic eigenfields
by low-degree polynomial
E(x), H(x) ~ e–it
— involves solving eigenequation Ax=2x

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 1        J
   E            H 
t                 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)                                                                     Hz

Hx
Ez              Hy                                                (i, j)        (i+1, j)
(i+1, j+1, k)
Ex
Ey
(i, j, k)                 (i+1, j, k)
Ex

Staggered grid in space:
— every field component is stored on a different grid
The Yee Discretization (1966)
(i, j+1)
H     1
   E 
t                                         Ey
Hz
H z                   1 Ey Ex 
           
t      1
i , j 
1        x   y                  (i, j)           (i+1, j)
2      2                                               Ex
             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
1         1
H(n  )  H(n  )
H                  1                          2         2
        E            
t   t nt                  t nt
t
+ O(∆t2)
Explicit time-stepping:
x
stability requires t 
# 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?
Order of Accuracy
E
TE polarization (E in plane: discontinuous)
a

QuickTime™ and a
TIFF (Uncompressed) decompressor
are neede d 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!
& convergence is
Three previous smoothing methods
still only linear
QuickTime™ and a
TIFF (Uncompressed) decompressor
are neede d to see this picture.

[ Dey, 1999 ]
[ Kaneda, 1997 ]
A Criterion for Accurate Smoothing
1st-order errors                       1    2 
~   E||  ( ) D 
2
from
smoothing                                 

We want the smoothing errors to be zero to 1st order
— minimizes error and 2nd-order convergent!
        Use a tensor : 
                        
                           E||
(in principal axes:)
               
             1 
[ Meade et al., 1993 ]                       1       E
Consistently the Lowest Error
a

TIFF (Uncompressed) decompressor
are neede d to see this picture.

a

[ Farjadpour et al., Opt. Lett. 2006 ]
A qualitatively different case: corners
still ~lowest error, but not quadratic
QuickTime™ and a
TIFF (Uncompressed) decompressor
are neede d 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?

• 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
t
n
1
P() 
2
 Re[ E   *
 H  ]dx
Transmission Spectra in FDTD
must always do two simulations: one for normalization

electric                                                                          P0()
field Ez:

QuickTime™ and a                    QuickTime™ and a
YUV420 codec decompressor           YUV420 codec decompressor
are needed to see this picture.     are needed to see this picture.

P()
transmission = P() / P0()
Reflection Spectra in FDTD
=1

 = 12

1
PR () 
2
             0               0
Re[(E  E )*  (H   H  )]dx

for reflection, subtract incident fields
(from normalization run)

1
P() 
2
 Re[ E   *
 H  ]dx
Transmission/Reflection Spectra
QuickTime™ see this
are needed toand a picture.
Graphics decompressor

1

0.9
=1
0.8

a   0.7                      = 12
0.6

0.5      T                        1–T–R
0.4

0.3

0.2

0.1          R
0
0.1   0.11    0.12   0.13    0.14   0.15   0.16   0.17   0.18   0.19   0.2

a / 2πc = a / l
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
— run ―numerical experiment‖ to                                                                      spectral methods
simulate E(x, t), H(x, t)
df   f ( x  x)  f ( x  x)
dx

x
 O(x 2 )      +
• 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
finite elements
• Frequency-domain eigensolver                                                                 in irregular ―elements,‖
approximate unknowns
— solve for source-free harmonic eigenfields
by low-degree polynomial
E(x), H(x) ~ e–it
— involves solving eigenequation Ax=2x

boundary-element methods
discretize only the
boundaries between
homogeneous media
A Maxwell Eigenproblem
1         
 E        H i H                            First task:
c t       c                      get rid of this mess
1         0   
 H        E  J  i E
c t           c

dielectric function (x) = n2(x)

1                      2
+ constraint
    H    H
         c                                H  0
eigen-operator            eigen-value   eigen-state
Electronic & Photonic Eigenproblems
Electronic                       Photonic
 2 2               1         2

           
  V   E     H    H
 2m                         c 
nonlinear eigenproblem          simple linear eigenproblem
(V depends on e density ||2)        (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
E
TM
H
Periodic Eigenproblems
if eigen-operator is periodic, then Bloch-Floquet theorem applies:


i k x t   
can choose:    H(x ,t)  e                    Hk (x )

planewave
periodic ―envelope‖

Corollary 1: k is conserved, i.e. no scattering of Bloch wave

Corollary 2: H k given by finite unit cell,
so  are discrete n(k)
A More Familiar Eigenproblem
=1

a                             = 12                             y
x
band diagram / dispersion relation

find the normal modes                          light cone
(all non-guided modes)
of the waveguide:
i(kxt)
H(y,t)  H k (y)e
(propagation constant k
a.k.a. )
k
Solving the Maxwell Eigenproblem
1                          n 2
Finite cell  discrete eigenvalues n                  ik    ik  H n                      Hn
                          c   2

Want to solve for n(k),
& plot vs. ―all‖ k for ―all‖ n,
QuickTime™ see this
are needed toand a picture.
Graphics decompressor
constraint:     ik H  0

1

where:         H(x,y) ei(kx – t)
0.9
0.8
0.7
0.6
0.5


0.4
Photonic Band Gap
0.3
0.2
TM bands
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
Solving the Maxwell Eigenproblem: 1
1   Limit range of k: irreducible Brillouin zone
—Bloch’s theorem: solutions are periodic in k
M

2     ky
first Brillouin zone                G        X    a
kx
= minimum |k| ―primitive cell‖


irreducible Brillouin zone: reduced by symmetry

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 )               ˆ H  2 H
solve: A
m1

finite matrix problem:   Ah   Bh  2

f g   f* g                 ˆ
Am  bm A b              Bm  bm b

3   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                               Finite-element basis
constraint, boundary conditions:
H(x t )   HG e            iGxt
                                    Nédélec elements
G                                                       [ Nédélec, Numerische Math.

HG  G  k  0
35, 315 (1980) ]
constraint:
nonuniform mesh,
uniform ―grid,‖ periodic boundaries,
more arbitrary boundaries,
simple code, O(N log N)               [ figure: Peyrilloux et al.,
J. Lightwave Tech.
21, 536 (2003) ]
complex code & mesh, O(N)

              3   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:
h' Ah
―variational
theorem‖            min
2
0
minimize by preconditioned
h h' Bh
are of 2dsee
QuickTime™ Model
Band Diagram needed toand athis picture. System
Graphics decompressor

a                    (radius 0.2a rods, =12)
1
0.9

frequency  (2πc/a) = a / l
0.8
0.7
0.6
0.5
0.4
Photonic Band Gap
0.3
0.2
TM bands
0.1
0
irreducible Brillouin zone                                      G       X        M           G
M
E             gap for
k             X                                           TM
G                                                          H           n > ~1.75:1
The Iteration Scheme is Important
(minimizing function of 104–108+ variables!)
h' Ah
  min
2
0            f (h)
h h' Bh

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
QuickTime™ see this
are needed toand a picture.
Graphics decompressor

(minimizing function of ~40,000 variables)
1000000
100000
10000
1000 Ñ
J            E       EE EE
100
J
Ñ
E EEEE
EE EE
E EEEEEE
no preconditioning
% error

J
Ñ Ñ                                 EEEE
EEEE
10             J ÑÑ
J Ñ
EEEEEEE
E EEEEE
EEEEEE
E EE EE
J ÑÑ                                        EE
EEE
EEE
EE
EE
EE
1                 J J ÑÑÑÑ                                        EE
EE
EE
EE
EE
EE
ÑÑÑ
J JJ ÑÑ                                         EE
EE
EE
EE
ÑÑÑÑÑÑ
JJJJJ ÑÑÑÑÑÑ
ÑÑÑÑÑÑ
ÑÑÑÑÑ
ÑÑÑÑ                           EE
EEE
EE
EE
0.1                            JJ           ÑÑÑ
ÑÑÑ
ÑÑÑ
ÑÑÑ                         EEEE
EE
ÑÑÑ
ÑÑÑ
ÑÑÑ
ÑÑ
ÑÑ                        EE
EE
EE
EE
J                    ÑÑ
ÑÑ
ÑÑ
ÑÑ                       EE
EE
0.01                               J                    ÑÑÑÑ                       EE
EE
EE
EE
J                     Ñ
ÑÑÑ
Ñ                         EE
EE
Ñ
ÑÑ
Ñ                        EE
E
J                       ÑÑ
Ñ
Ñ                         EE
E
0.001                                                           Ñ
Ñ                         EE
preconditioned J        J                         ÑÑ
Ñ
ÑÑ
Ñ
EE
E
EE
EE
E
E
0.0001
EEE
EE
EE
E
0.00001                                    J
J
0.000001
1                     10                       100                      1000

# iterations
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 :
                 
                   E||
                 
               1 
         1       E
The -averaging is Important
QuickTime™ see this
are needed toand a picture.
Graphics decompressor

100

H
H       backwards averaging
10                H H H
H
J
B             H H H
H   H H
correct averaging
B                            H
% error

J   J           B        B                   changes order
averaging
B no B
J B    B
1                       J    J
B
B                 B     of convergence
B
J
J   J
B
from ∆x to ∆x2
tensor averaging                J   J
0.1
J J

0.01
10                                                 100     (similar effects
resolution (pixels/period)                         in other E&M
numerics & analyses)
Gap, Schmap?     QuickTime™ see this
are needed toand a picture.
Graphics decompressor

1
a                         0.9
0.8
0.7

frequency 
0.6
0.5
0.4
Photonic Band Gap
0.3
0.2
TM bands
0.1
0
G            X           M           G

But, what can we do with the gap?
Intentional ―defects‖ are good

microcavities                    waveguides (―wires‖)

3D Pho to nic C rysta l with De fe c ts
Intentional ―defects‖ in 2d
QuickTime™ see this
are needed toand a picture.
Graphics decompressor

(Same computation, with supercell = many primitive cells)
QuickTime™
GraphicsGraphics a toand a picture.
QuickTime™ see this picture.
are needed toand decompressor
are needed see this
decompressor

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                      Response is many
complex n                                               sharp peaks,
[ Mandelshtam,
J. Chem. Phys. 107, 6756 (1997) ]        one peak per mode

decay rate in time gives loss
Signal Processing is Tricky
signal processing
complex n
?
spectrum
a common approach: least-squares fit of picture.
QuickTime™ see this
are needed toand a picture.
Graphics decompressor              QuickTime™ see this
are needed toand a
Graphics decompressor

0.8                                                    450
E
0.6 E E                                               400
E
0.4   EE
E
350                  fit to:
E E
0.2 E EE    E
E
EE EEEE
E E EE E EE
300      EE
A
E EE E EE EEE EEEEEEEEEEEEEEEE
0 E E E E E EEEEEEEEEEEEEEEEEEEEEEEEEEE
E E EEE EEEEEEEEEEEEEEE
(   0 ) 2  G 2
E E E                                  250
E EE E
E E   EE
E
-0.2 E EE                                              200
-0.4
E
E
E
E
E
FFT   150
E       E
E
-0.6                                                   100       E E
E
-0.8                                                   50    E
E   E
E
E
-1 E                                                  0  EE E
E        EE E
E E E EE EEE E EE
E E EE E EEE EEE
0 1 2 3 4 5 6 7 8 9 10                               0 0.5 1 1.5 2 2.5 3 3.5 4

Decaying signal (t)                                Lorentzian peak ()
Fits and Uncertainty
completely
problem: have to run long enough to and athis picture. decay
are needed toand a picture.
QuickTime™ see this
Graphics decompressor              QuickTime™
Graphics decompressor
are needed to see

1  E                                             40000
E E
E
0.8 E E
E EE     E
E E                                  35000
E E
0.6 EE
E
E E
E E
E      E
E
E   E E
E E
E E
actual
EE    EE EE        E   E E E E
E E E E               30000
0.4                      EE           E
E       E
EE EE EE EE E    E
EE E E   E E EE                    E EE EE     25000
0.2                      EE  E E E E E E E E EE
E
0 EEEEEEEEE             EEEEEEEEEEE              20000
E E
-0.2            EE E E    E E E E E E E E EE E                           E
E EE EE                        E EE EE E
E EE EE E E E
15000
E             E E
-0.4
EE   EE EE           E E
E E E E
E E       E E E E
10000
-0.6  E EE
E   E E
E EE
E
E E
E                                 signal
E E    E E
E
E E
-0.8 E E E
E
E    E                                 5000
portion
E E
-1 E                                               0E E E E E             E E E E E
0 1 2      3    4    5   6   7   8   9   10      0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5

Portion of decaying signal (t)               Unresolved Lorentzian peak ()

There is a better way, which gets complex  to > 10 digits
Unreliability of Fitting Process
Resolving two overlapping peaks is
near-impossible 6-parameter nonlinear fit
QuickTime™ see this
are needed toand a picture.
Graphics decompressor

(too many local minima to converge reliably)
1200
EE
1000
sum of two peaks
There is a better
800                               E
E                                way, which gets
600                                                               complex 
for both peaks
400         = 1+0.033i   E        E                             to > 10 digits
200
E           E  = 1.03+0.025i
E             E
EE               EE
0 EEEEEEEEEEEEEEEEEEE                   EEEEEE
EEEEEEEEEEEE
0.5 0.6    0.7 0.8 0.9      1       1.1  1.2 1.3 1.4 1.5
Sum of two Lorentzian peaks ()
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)   ak e                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:
ˆ  i  
H                         time-∆t evolution-operator:           ˆ  eiHt /
U
ˆ

t
say:                               ˆ n  (0)
y n   (0)  (nt)   (0) U
Method
Filter-Diagonalization107, 6756 (1997) ] (FDM)
[ Mandelshtam, J. Chem. Phys.

ˆ
y n   (0)  (nt)   (0) U n  (0)        ˆ  eiHt /
U
ˆ

We want to diagonalize U: eigenvalues of U are ei∆t

…expand U in basis of |(n∆t)>:

ˆ                 ˆ ˆˆ
U m,n   (mt) U  (nt)   (0) U mUU n  (0)  y m 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,              • Arbitrary periodic(x) —
and nonlinear [(2) and (3)]       anisotropic, magneto-optic, …
(lossless, linear materials)
• ArbitraryJ(x,t)
• PML/periodic/metal bound.                            • 1d/2d/3d

• 1d/2d/3d/cylindrical               • band diagrams, group velocities
perturbation theory, …
• power spectra • eigenmodes

Free/open-source          • fully scriptable interface
software (GNU)
• built-in multivariate optimization,
• MPI parallelism         integration, root-finding, …
• exploit mirror symmetries       • field output (standard HDF5 format)
Unix Philosophy
combine small, well-designed tools, via files

Input text file       MPB/Meep           standard formats
(text + HDF5)
— have to learn several programs
Visualization / Analysis
— flexibility                          (Matlab, Mayavi [vtk],
— batch processing, shell scripting   command-line tools, …)
— ease of development
Unix Philosophy
combine small, well-designed tools, via files

Input text file      MPB/Meep           standard formats
(text + HDF5)

GNU Guile scripting interpreter
(Scheme language)
Visualization / Analysis
software
Embed a full scripting language:       (Matlab, Mayavi [vtk],
— parameter sweeps                    command-line tools, …)
— complex parameterized geometries
— optimization, integration, etc.
— programmable J(x, t), etc.
— … Turing complete
A Simple Example (MPB)
=1

a                         = 12                    y
find the normal modes n(k)               x
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, … ?)
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

10 (320 pixels)
=1

a                                     = 12                                  y
x
A File Format Made of Parentheses
Need to specify:
• computational cell size/resolution
• geometry, i.e. (y)
(set! geometry                                  (choose   units of a)
(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)                  …Then run:
• what k values                        (run)
• how many modes (n = 1, 2, … ?)
(set! num-bands 5)                      or only TM polarization:
(run-tm)

or only TM, even modes:
(run-tm-yeven)

=1

a                             = 12                    y
x
Simple Example (MPB) Results
=1

a                            = 12             y
find the normal modes n(k)       x
of the waveguide:

red = even
blue = odd
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/

```
To top