Transcript From electrons to photons: Quantum

From electrons to photons: Quantuminspired modeling in nanophotonics
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
wavevector
interacting: hard problem
photon frequency
electron energy
Bloch waves:
Band Diagram
Periodic
Medium
atoms in diamond structure
wavevector
non-interacting: easy problem
Electronic & Photonic Modelling
Electronic
Photonic
• strongly interacting
—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
Option 1: Numerical “experiments” — discretize time & space … go
Option 2: Map possible states & interactions
using symmetries and conservation laws: band diagram
Fun with Math
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)
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
D
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(nDt)   ake
i k nDt
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ˆ Dt /
ˆ
U e
yn  (0) (nDt)  (0) Uˆ n (0)
Filter-Diagonalization
Method
(FDM)
[ Mandelshtam, J. Chem. Phys. 107, 6756 (1997) ]
yn  (0) (nDt)  (0) Uˆ n (0)
iHˆ Dt /
ˆ
U e
We want to diagonalize U: eigenvalues of U are ei∆t

…expand U in basis of |(n∆t)>:
Um,n  (mDt) Uˆ (nDt)  (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
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/