4/27/2015
Instructor
Dr. Raymond Rumpf
(915) 747‐6958
[email protected]
EE 4395/5390 – Special Topics
Computational Electromagnetics
Lecture #5
Transfer Matrix Method Using Scattering Matrices
 These notes may contain copyrighted material obtained under fair use rules. Distribution of these materials is strictly prohibited 
Lecture 5
Slide 1
Outline
•
•
•
•
•
Alternate matrix differential equation (no mode sorting)
Scattering matrix for a single layer
Multilayer structures
Calculating reflected and transmitted power
Notes Implementation
• Advanced scattering
matrices
• Alternatives to
scattering matrices
Lecture 5
Bonus
Slide 2
1
4/27/2015
Alternate Matrix Differential Equation
Lecture 5
Slide 3
Recall Derivation from Last Lecture
Start with Maxwell’s equations from Lecture 2.
Assume LHI.
Ez E y

 k0  r H x
y
z
Ex Ez

 k0  r H y
z
x
E y Ex

 k0  r H z
x
y
H z H y

 k0  r E x
y
z
H x H z

 k0  r E y
z
x
H y H x

 k0  r E z
x
y
Lecture 5
Assume device is infinite and uniform in x and y directions.

 jk x
x
jk y Ez 

 jk y
y
dE y
dz
 k0 r H x
dE x
 jk x E z  k0 r H y
dz
jk x E y  jk y E x  k0  r H z
jk y H z 
dH y
dz
 k 0 r E x
dH x
 jk x H z  k0 r E y
dz

jk x H y  jk y H x  k0 r Ez
Normalize z and wave vectors kx, ky, and kz.
z   k0 z
k
kx  x
k0
ky
ky 
k0
Eliminate longitudinal components Ez and Hz by substitution.
k
kz  z
k0
dE y
jky Ez 
 r H x
dz 
dEx
 jkx Ez  r H y
dz 
jkx E y  jky Ex   r H z
dH y
jky H z 
  r Ex
dz 

dH x
 jkx H z   r E y
dz 
jkx H y  jky H x   r Ez

dEx kx ky
k 2 
H x   r  x  H y

dz 
r
r 


dE y  ky2
kx ky
Hy
   r  H x 

dz    r
r


dH x kx ky
k 2 
Ex    r  x  E y

dz 
r
r 


dH y  ky2
kx ky
E
    r  Ex 

dz    r
r y

Slide 4
2
4/27/2015
Derivation of Two 22 Matrix Equations
We can write our two sets of two equations in matrix form as
dEx kx ky  
k 2 

H x   r  x  H y
dz 
r
r 


dE y  ky2
k k
   r  H x  x y H y

r
dz    r

d  Ex  1
 
dz   E y   r

dH x kx ky
k 2 
Ex    r  x  E y

dz 
r
r 

2



dH y  k y
k x ky
E
    r  Ex 

dz   r
r y

d  H x  1  kx ky

 
dz   H y   r  ky2  r  r
 kx ky
 2
 ky   r  r
r  r  kx2   H x 
 
kx ky   H y 
r  r  kx2   Ex 
 
 kx ky   E y 
Note: These equations are valid regardless of the sign convention.
Lecture 5
Slide 5
Compact “PQ” Form
We can write our two matrix equations more compactly as
d  Ex  1
 
dz   E y   r
P
 kx ky
 2
 ky   r  r
1  kx ky

 r  ky2  r  r
d  H x  1  kx ky

 
dz   H y  r  ky2   r  r
r  r  kx2   H x 
 
kx ky   H y 
r  r  kx2 
 kx ky


 H x 
d  Ex 
   P  
dz   E y 
 H y 
Q
r  r  kx2   Ex 
1  kx ky

r  ky2  r  r
 
kx ky   E y 
r  r  kx2 

kx ky 
 Ex 
d  H x 
  Q 
dz   H y 
 Ey 
Note: We will see this same “PQ” form again for other methods like MoL, RCWA, and waveguide analysis. TMM, MoL, and RCWA are implemented the same after P and Q are calculated.
Lecture 5
Slide 6
3
4/27/2015
Matrix Wave Equation
Our two governing equations are
 H 
d  Ex 
 P   x  Eq. (1)


E
dz   y 
 H y 
 Ex 
d  H x 
   Q 
dz   H y 
Ey 
Eq. (2)
We can now derive a matrix wave equation. First, we differentiate Eq. (1) with respect to z.
 H 
d d  Ex  d


P   x 


E
dz  dz   y  dz   H y 

d2
dz 2
 Ex 
d  H x 
E   P    
dz  H y 
 y
Second, we substitute Eq. (2) into this result.
d2
dz 2
 Ex 
 Ex 
E   P Q E 
 y
 y
E
 Ex 
0
2  x
E   Ω E    
 y
 y  0
Ω 2  PQ
d2
dz 2
Lecture 5
Slide 7
Numerical Solution (1 of 3)
The system of equations to be solved is
d2
dz 2
E
 Ex 
0 
2  x
E   Ω E    
 y
 y  0 
Ω 2  PQ
This has the general solution of
 Ex  z 
Ωz  
 Ωz  

 e a e a

 E y  z 
a   proportionality constant of forward wave
a   proportionality constant of backward wave
Note: Here we are finally assuming a sign convention of e‐jz for forward propagation in the +z direction.
No mode sorting!  Here, we solved a second‐order differential equation where the modes we calculate are one‐way. We simply write them twice for forward and backward waves. Before we solved a first‐order differential equation that lumped forward and backward modes together.
Lecture 5
Slide 8
4
4/27/2015
Numerical Solution (2 of 3)
Recall from Lecture 4…
f  A   W  f  λ   W 1
We can use this relation to compute the matrix exponentials.
e  Ωz  We  λz W 1
W  Eigen-vector matrix of Ω 2
λ 2  Eigen-value matrix of Ω 2
eΩz  We λz W 1
e λz 
e






12 z 
e
22 z 







 2 z
e N 
So the overall solution can now be written as
 Ex  z 
 λz 
1 
1 
λz 

  We W a  We W a

E
z
 y  
Lecture 5
Slide 9
Numerical Solution (3 of 3)
So the overall solution can now be written as
 Ex  z  
1 
1 
 We  λz 
W
a  We λz  
W
a






 E y  z
c
c
The column vectors a+ and a‐ are proportionality constants that have not yet been determined.
The eigen‐vector matrix W multiplies a+ and a‐ to give another column vector of undetermined constants.
To simplify the math, we combine these products into new column vectors labeled c+ and c‐ .
 Ex  z 
 λz  
λz  

  We c  We c
 E y  z
Lecture 5
Slide 10
5
4/27/2015
Analytical Expressions for W and 
The dispersion relation with normalized wave vectors is
r  r  kx2  ky2  kz2
Using this relation, we can simplify the matrix equation for 2.
Ω 2  PQ 
1  kx ky

r  r  ky2  r  r
r  r  kx2   kx ky
kx ky
 2

  k y  r  r
r  r  kx2   kz2
kx ky

  0
0 
2
  kz I
 kz2 
1 0 
I
 identity matrix
0 1 
Since this is a diagonal matrix, we can conclude that
1 0 
WI

0 1 
λ 2  Ω2
 jk
λ z
 0
0 
  jkz I
jkz 
 e jkz z
e λz  
 0
0 


e jk z z 
We don’t actually have to solve the eigen‐value problem to obtain the eigen‐modes! 
Lecture 5
Slide 11
Solution for the Magnetic Field (1 of 2)
The magnetic field will have a similar solution, but will have its own eigen‐vector matrix V to describe its modes.
 H x  z   
 λz  
λz  

  Ve c  Ve c

H
z
 y   
We put the minus sign in the solution here so that both terms in the differentiated equation will be positive. Since the electric and magnetic fields are coupled and not independent, we should be able to compute V from W. First, we differentiate the solution with respect to z’.
d  H x  z   
 λz  
λz  

  Vλe c  Vλe c
dz  H y  z   
Lecture 5
Slide 12
6
4/27/2015
Solution for the Magnetic Field (2 of 2)
We now have
d  H x  z  
 λz  
λz  

  Vλe c  Vλe c
dz  H y  z   
Recall that
 Ex  z   
d  H x  z   

  Q

dz  H y  z  
 E y  z 
and
 Ex  z   
 λz  
λz  

  We c  We c

E
z
 y  
Combining these results leads to

Vλe  λzc  Vλe λzc  Q We  λzc   We λzc 

 QWe  λzc   QWeλzc 
Comparing the terms shows that
Vλ  QW 
V  QWλ 1
Lecture 5
Slide 13
Combined Solution for E and H
Electric Field Solution
 Ex  z  
λz  
 λz  

  We c  We c
 E y  z
1 0 
W

0 1 
c   amplitude coefficients of forward wave
c   amplitude coefficients of backward wave
W  eigen-vector matrix
Magnetic Field Solution
 H x  z  
λz  
 λz  

  Ve c  Ve c
 H y  z   
V  QWλ 1
Combined Solution
 Ex  z   


E y  z    W W  e  λz

ψ  z   
 H x  z     V V   0


 H y  z   
Lecture 5
0  c  
 
e λz  c  
Does this equation look familiar?
This is the same equation we had at the end of Lecture 4.
Slide 14
7
4/27/2015
Two Paths to Combined Solution
4×4 Matrix
   yz
 kx ky
 
  

 
jkx  yz  zy 
 kx zx 
  yx  yz zx 
  j  ky

  zz
 zz 
 zz 
   zz
  zz  zz 


2



k









y
 Ex 
xz  zx
xz
zx

 xz  k zy





jk
j
k






xx
y
x
y

 

 zz 
 zz 
  zz  zz 
  zz
  Ey  
  zz




2
H



z
k k
x
    
  
   yz   zx 
    x y   yx  yz zx    k x   yy  yz zy 
j
k
k


 y

x
 zz    zz
 zz 
 zz 
 H y    zz
  zz

  k 2
  k k
  

 
  y   xx   xz zx    x y   xy  xz zy 
jky  xz  zx 




 zz    zz
 zz 
   zz
  zz  zz 
Sort Eigen‐Modes


 


  E 
  x 

  Ey 
  H 
  yz  zy 
 x
jkx 


  H y 
  zz  zz 

   xz   zy  
 j  kx
 ky
 
 zz  
  zz
  kˆx2
 
k
  yy  yz zy
 
 zz
zz

 kx ky
 xz  zy
  xy 


 zz
  zz
 W
W   E
 WH
eλ z

 0

e λz
WE 

WH 
Ex
Ey
H x
H
y
0 

 λ  z
e

Anisotropic
Field Solution
Maxwell’s Equations


  E  k0  r  H


  H  k0  r  E
Isotropic or diagonally anisotropic
Lecture 5
 W
ψ  z    E
 VH
WE  e

VH   0
 λ  z
 W W  e  λz
ψ  z   

 V V   0
r  r  kx2 
1  kx ky
 2


 r  k y   r  r
 kx ky 
 r  r  kx2 
1  kx ky
Q
 2


 r  k y   r  r
 kx ky 
0  c  
 

e λ z  c 
0  c  
 
e λz  c  
P
No sorting! 
PQ Method
Slide 15
Scattering Matrix for a Single Layer
R. C. Rumpf, “Improved Formulation of Scattering Matrices for Semi‐Analytical Methods That is Consistent with Convention,” PIERS B, Vol. 35, pp. 241‐261, 2011.
Lecture 5
Slide 16
8
4/27/2015
Motivation for Scattering Matrices
Scattering matrices offer several important features and benefits:
• Unconditionally stable method.
• Parameters have physical meaning.
• Parameters correspond to those measured in the lab.
• Can be used to extract dispersion.
• Very memory efficient.
• Can be used to exploit longitudinal periodicity.
• Mature and proven approach.
• Much greater wealth of literature available.
Excellent alternatives to S‐matrices do exist!
Lecture 5
Slide 17
Geometry of a Single Layer
Indicates a point that lies on an interface, but associated with a particular side.
Layer i
Medium 1
ψ1
ψ i  0 
ψ1
ψ i  0 
ψ  zi 

i
ψ i  zi 
c1 ci
Medium 2
ψ i  k0 Li 
ψ 2
ψ i  k0 Li 
ψ 2
c2
c2
c1 ci
Li
ψ i  z   field within i th layer
Lecture 5
z
ci  mode coefficients inside i th layer
ci   mode coefficients outside i th layer
Slide 18
9
4/27/2015
Field Relations
Field inside the ith layer:
 Ex ,i  zi  


E y ,i  zi    Wi
ψ i  zi    

 H x ,i  zi     Vi


 H y ,i  zi  
Wi   e  λi zi

Vi   0
0  ci 
 
e λi zi  ci 
Boundary conditions at the first interface:
ψ1  ψ i  0 
 W1
 V
 1
W1  c1   Wi
 
V1  c1   Vi
Wi  ci 
 
Vi  ci 
Boundary conditions at the second interface:
ψ i  k0 Li   ψ 2
 Wi
 V
 i
Wi  e  λ i k0 Li

Vi   0
0  ci   W2
   
eλ i k0 Li  ci   V2
W2  c2 
 
V2  c2 
Note: k0 has been incorporated to normalize Li.
Lecture 5
Slide 19
Definition of A Scattering Matrix
c1   S11 S12  c1 
   
 
c2  S 21 S 22  c2 
S 21
c

1
S11
c1
Lecture 5
S11  reflection
S 21  transmission
transmission
reflection
c2
S 22
S12
This is consistent with experimental convention.
c2
Slide 20
10
4/27/2015
Derivation of the Scattering Matrix
Solve both boundary condition equations for the intermediate mode ci
ci
coefficients and .
1
ci   Wi
   
ci   Vi
Wi   W1
Vi   V1
ci  e λ i k0 Li
   
c i   0
W1  c1 
 
V1  c1 
0
e
 λ i k0 Li
  Wi

  Vi
1
Wi   W2
Vi   V2
W2  c2 
 
V2  c2 
Both of these equations have
 Wi

 Vi
1
Wi   W j

Vi   V j
W j  1  A ij

V j  2  Bij
A ij  Wi1W j  Vi1V j
Bij 
A ij 
Bij  Wi1W j  Vi1V j
We set substitute this result into the first two equations and set them equal.
1  Ai1
2  Bi1
Bi1  c1  1 e λ i k0 Li
  
Ai1  c1  2  0
0
e
 λ i k0 Li
Bi 2  c2 
 
Ai 2  c2 
  Ai 2

  Bi 2
We write this as two matrix equations and rearrange terms until they have the form of a scattering matrix.
c1  ? ? c1 
   
 
c2  ? ? c2 
Lecture 5
Slide 21
The Scattering Matrix
The scattering matrix Si of the ith
layer is defined as:

c1 
 i  c1 
  S  
c2 
c2 
i 
S11
S  i

S 21
 
S12

S 22i  
i
i 
After some algebra, the components of the scattering matrix are computed as

S    A
S    A
S    A
i
 A i1  Xi Bi 2 A i21Xi Bi1
S11
i
12
i1
 Xi Bi 2 A i21Xi Bi1
i
21
i2
 Xi Bi1A i11Xi Bi 2
i
22
Lecture 5
1
i 2  Xi B i1 A i1 X i B i 2
 X B A X A
 X A  B A
 X A  B A
 X B A X A
1
i
i2
1
i2
i
i1
1
Si
 Bi1
i
i2
i2
1
i2
i
i1
i1
1
i1
1
1
i
i1
1
i1
i



Bij  Wi1W j  Vi1V j
Bi 2
Bi1
i 2  Bi 2
A ij  Wi1W j  Vi1V j

X i  e  λ i k0 Li
i is the layer number.
j is either 1 or 2 depending on which external medium is being referenced.
Slide 22
11
4/27/2015
Scattering Matrices in the Literature
For some reason, the computational electromagnetics community has: (1) deviated from convention, and (2) formulated scattering matrices inefficiently.
S 22 transmission
c   S11 S12   c 
 
 
 c  S 21 S 22  c 

i 1

i

i

i 1
S12 reflection
S 21
S11
ci1
ci
ci1
ci
Lecture 5
• Here s11 is not reflection. Instead. it is backward transmission!
• Here s21 is not transmission. Instead, it is a reflection parameter!
• Scattering matrices can not be interchanged.
• Scattering matrices are not symmetric so they take twice the memory to store and are more time‐consuming to calculate.
Slide 23
Limitation of Conventional S‐Matrix Formulation
Note that the elements of a scattering matrix are a function of materials outside of the layer.
This makes it difficult to interchange scattering matrices arbitrarily.
For example, there are only three unique layers in the multilayer structure below, yet 20 separate computations of scattering matrices are needed.
Three unique layers
20 layer stack
Lecture 5
Slide 24
12
4/27/2015
Solution
To get around this, we will surround each layer with an external regions of zero thickness. This lets us connect the scattering matrices in any order because they all calculate fields that exist outside of the layers in the same medium. This will have no effect electromagnetically as long as we make the external regions have zero thickness.
Gap Medium
Layer i
Gap Medium
Li
Lecture 5
Slide 25
Visualization of the Technique
We calculate the scattering matrices for just the unique layers.
Three unique layers
Then we just manipulate these same three scattering matrices to “build” the global scattering matrix.
Gaps between the layers are made to have zero thickness so they have no effect electromagnetically.
Faster! Simpler! Less memory needed!
Lecture 5
Slide 26
13
4/27/2015
Revised Geometry of a Single Layer
Gap Medium
Gap Medium
Layer i
ψ1
ψ i  0 

1
ψ 0
ψ i  zi 
ψ i  k0 Li  ψ 2
r ,h
r ,h
ψ
 r ,h

i
ψ i  zi 
ψ

i
 k0 Li 
c1 ci
ψ

2
 r ,h
c2
c2
c1 ci
Li
Lecture 5
Slide 27
Calculating Revised Scattering Matrices
The scattering matrix Si of the ith
layer is still defined as:

c1 
 i  c1 

S
 
 
c2 
c2 
S  i 
S  i    11i
S 21
i  
S12

S 22i  
But the elements are calculated as

S    A
i
12
i
i 
S21i   S12
i 
i 
S 22  S11
Lecture 5
 Xi Bi A i1
Si
 X B A X A  B 
X B  X A  B A B 
i 
S11
 A i  Xi Bi A i1Xi Bi
1
i
1
i
i
i
1
i
i
i
i
i
i
1
i
i
i
• Layers are symmetric so the scattering matrix elements have redundancy.
• Scattering matrix equations are simplified.
• Fewer calculations.
• Less memory storage.
A i  Wi1Wh  Vi1Vh
Bi  Wi1Wh  Vi1Vh
X i  e  λ i k0 Li
Slide 28
14
4/27/2015
Scattering Matrices of Lossless Media
If a scattering matrix is composed of materials that have no loss and no gain, the scattering matrix must conserve energy. That is, all incident energy must either reflect or transmit.
This implies that the scattering matrix is unitary.
If the scattering matrix is unitary, it must obey the following rules:
S H  S 1
S H S  SS H  S 1S  SS 1  I
Lecture 5
Slide 29
Hints About Stability in These Formulations
• Diagonal elements S11 and S22 tend to be the largest numbers. Divide by these instead of any off‐diagonal elements for best numerical stability.
 S11 S12 
S

S 21 S 22 
• X describes propagation through an entire layer. Don’t divide by X or your code can become unstable.
Lecture 5
Slide 30
15
4/27/2015
Multilayer Structures
Lecture 5
Slide 31
Solution Using Scattering Matrices
The scattering matrix method consists of working through the device one layer at a time and calculating an overall scattering matrix.
S 1
S 2
S  3
S 4
S 5
S  device   S1  S 2   S  3  S 4   S  5
Redheffer star product.
NOT matrix multiplication!
Lecture 5
Slide 32
16
4/27/2015
Redheffer Star Product
Two scattering matrices may be combined into a single scattering matrix using Redheffer’s star product.
S  A 
S A   11A
S21 
S AB   S A  S B 
A 
S12

S22A  
S  B
S B    11B
S21 
 B 
S12

S22B 
The combined scattering matrix is then

S11
AB 
S
 AB 
 AB
S11
  AB
S21 
 AB 
S12

S22AB 
1
 
 
   
   
 S11
 S12
I  S11 S 22  S11 S 21
A
A
B
A
B
A
1
 AB 
A 
 B  A  
 B
 S12
S12
I  S11 S 22  S12
1
 
 
S21   S21  I  S22 S11
 S 21
AB
B
A
B
A
1
 B 
 A   B
S22AB  S22B  S21B I  S 22A S11
 S 22 S12
R. Redheffer, “Difference equations and functional equations in transmission-line theory,”
Modern Mathematics for the Engineer, Vol. 12, pp. 282-337, McGraw-Hill, New York, 1961.
Lecture 5
Slide 33
Derivation of the Redheffer Star Product
We start with the equations for the two adjacent scattering matrices.
A
c1  S11


 
A
c2  S 21
A 
S12
c1 

 
A
S22   c2 
 B
c2  S11


 
 B
c3  S 21
 B 
S12
c2 

 
B
S 22   c3 
We expand these into four matrix equations.
A 
A 
 B 
 B 
Eq. 1
Eq.  3
c1  S11
c1  S12
c2
c2  S11
c2  S12
c3
A 
A 
 B 
 B 


Eq.  2 
Eq.  4 
c2  S 21 c1  S 22 c2
c3  S 21 c2  S 22 c3
c2
We substitute Eq. (2) into Eq. (3) to get an equation with only .
c2
We substitute Eq. (3) into Eq. (2) to get an equation with only .
 I  S S   c
 I  S   S    c
B
11
A
22

2
 B  A  
 B 
 S11
S 21 c1  S12
c3
Eq.  5 
A
22
B
11

2
 B 
 S 21A c1  S 22A S12
c3
Eq.  6 
c
c
We eliminate and by substituting these equations into Eq. (1) and (4). We then rearrange terms into the form of a scattering matrix.

2
c1  ? ? c1 
   
 
c3  ? ? c3 
Lecture 5

2
Overall, this is just algebra. We start with 4 equations and 6 unknowns and reduce it to 2 equations with 4 unknowns.
Slide 34
17
4/27/2015
Putting it All Together
We have one remaining problem. In the revised framework, the global scattering matrix places the device in free space. In many applications, we may want something other than free space outside the device.
We connect the global scattering matrix to the external materials by surrounding it by “connection” scattering matrices that have zero‐thickness
Device exists within gap medium
S
global 
 S
ref 
 S    S      S     S 





1
2
N
trn 
Device in gap medium
S device 
Lecture 5
Slide 35
Reflection/Transmission Side Scattering Matrices
The reflection‐side scattering matrix is
 ref 
1
ref
 ref 
1
ref
S11   A B ref
S12  2 A
s ref
A ref  Wh1Wref  Vh1Vref

S 21ref   0.5 A ref  B ref A ref1 B ref

1
h
r ,I
r ,h
 r ,I
 r ,h
1
h
B ref  W Wref  V Vref
S 22ref   B ref A ref1
lim
L0
The transmission‐side scattering matrix is
 trn 
S11  B trn A
 trn 

1
trn
1
trn
S12  0.5 A trn  B trn A B trn
S21   2A trn1
trn
S22trn    A trn1 B trn

s trn
A trn  Wh1Wtrn  Vh1Vtrn
1
h
r ,h
r ,II
 r ,h
 r ,II
1
h
B trn  W Wtrn  V Vtrn
lim
L0
Lecture 5
Slide 36
18
4/27/2015
Calculating Transmitted and Reflected Power
Lecture 5
Slide 37
Recall How to Calculate Source Parameters Incident Wave Vector
sin  cos  

kinc  k0 ninc  sin  sin  
 cos  
Right‐handed
coordinate system
x
y
z
Surface Normal
0
nˆ  0
1 
Unit Vectors Along Polarizations
 aˆ y


aˆTE   nˆ  kinc

 nˆ  k
inc


aˆTE  kinc
aˆTM 

aˆTE  kinc
  0
  0
Composite Polarization Vector

P  pTE aˆTE  pTM aˆTM
In CEM, we usually make

P 1
Lecture 5
Slide 38
19
4/27/2015
Solution Using Scattering Matrices
The external fields (i.e. incident wave, reflected wave, transmitted wave) are related through the global transfer matrix.
c ref 
 global  cinc 
c   S
 0 
 
 trn 
E
1  x ,inc 
cinc  Wref
E 
 y ,inc 
We get Ex,inc and Ey,inc from the polarization vector P.  E   P 
x ,inc
x
E   P 
 y ,inc   y 
This matrix equation can be solved to calculate the mode coefficients of the reflected and transmitted fields.
 global 
c ref  S11
 c     global
 trn  S 21
 global  
S12
cinc 

 global   0 
S 22   
 global 
c ref  S11
cinc
c trn  S21globalcinc
cinc 
right
cinc
not typically used
Lecture 5
Slide 39
Calculation of Transmitted and Reflected Fields
The procedure described thus far calculated cref and ctrn. The transmitted and reflected fields are then
inc
 Exref 
 global 
 global  1  Ex 
 ref   Wref c ref  Wref S11 cinc  Wref S11 Wref  inc 
 E y 
 E y 
inc
 Extrn 
 global 
 global  1  Ex 
 trn   Wtrn c trn  Wtrn S 21 cinc  Wtrn S 21 Wref  inc 
 E y 
 E y 
Lecture 5
Slide 40
20
4/27/2015
Calculation of the Longitudinal Components
We are still missing the longitudinal field component Ez on either size of the layer stack.
These are calculated using Maxwell’s divergence equation. 
E  0


E0, x e
x
 
 jk  r
  y  E
0, y e
 
 jk  r
 
  z  E
0, z e
 
 
 jk  r
Note:

   E  0 reduces to

  E  0 when  is homogeneous.
 
0
 
 jk x E0, x e  jk r  jk y E0, y e  jk r  jk z E0, z e jk r  0
k x E0, x  k y E0, y  k z E0, z  0
k z E0, z   k x E0, x  k y E0, y
E0, z  
k x E0, x  k y E0, y
kz
Ezref  
kx E xref  ky E yref
k ref
z
Eztrn
kx Extrn  ky E ytrn

k trn
z
Lecture 5
Slide 41
Calculation of Power Flow
Reflectance is defined as the fraction of power reflected from a device.
 2
Eref
2
2
2
2
R  2
E  Ex  E y  Ez
Einc
Transmittance is defined as the fraction of power transmitted through a device.
 2
Etrn

k 
T   2 Re  r ,ref z ,trn 
  r ,trn k z ,inc 
Einc
Note: We will derive these formulas in Lecture 7.
It is always good practice to check for conservation of energy.
 1 materials have loss
R  T   1 materials have no loss and no gain
 1 materials have gain
Lecture 5
Slide 42
21
4/27/2015
Reflectance and Transmittance on a Decibel Scale
Decibel Scale
PdB  20 log10  A 
PdB  10 log10  P 
How to calculate decibels from an amplitude quantity.
How to calculate decibels from a power quantity.
P  A2
PdB  10 log10  A2   20 log10  A 
Reflectance and Transmittance
Reflectance and transmittance are power quantities, so
RdB  10 log10  R 
TdB  10 log10 T 
Lecture 5
Slide 43
Notes on
Implementation
Lecture 5
Slide 44
22
4/27/2015
Storing the Problem
How is a the device described and stored for TMM?
We don’t use a grid for this method!
Store the permittivity for each layer in a 1D array.
Store the permeability for each layer in a 1D array.
Store the thickness of each layer in a 1D array.
ER = [ 2.50 , 3.50 , 2.00 ];
UR = [ 1.00 , 1.00 , 1.00 ];
L = [ 0.25 , 0.75 , 0.89 ];
Input arrays for three layers
We will also need the external materials, and source parameters.
er1, er2, ur1, ur2, theta, phi, pte, ptm, and lam0
Lecture 5
Slide 45
Storing Scattering Matrices
We often talk about the scattering matrix S as a single matrix.
S12 
S
S   11

S 21 S 22 
However, we never actually use the scattering matrix S this way. We only ever use the individual terms S11, S12, S21, and S22.
So, scattering matrices are actually stored as the four separate components of the scattering matrix.
S12 
S
S   11

S 21 S 22 
Lecture 5

S11 , S12 , S 21 , and S 22
Slide 46
23
4/27/2015
Calculating Xi = exp(-ik0Li)
Recall the correct answer:
Xi  e
 Ωi k0 Li
 e  jkz k0 Li

 0



e  jk z k0 Li 
0
It is incorrect to use the function exp() because this calculates a point‐by‐point exponential, not a matrix exponential.
X =
0.0135 + 0.9999i
1.0000
X = exp(-OMEGA*k0*L);
1.0000
0.0135 + 0.9999i
Approach #2: diag()
Approach #1: expm()
X = expm(-OMEGA*k0*L);
X = diag(exp(-diag(OMEGA)*k0*L));
X =
0.0135 + 0.9999i
0
X =
0.0135 + 0.9999i
0
0
0.0135 + 0.9999i
0
0.0135 + 0.9999i
Lecture 5
Slide 47
Efficient Star Product
After observing the equations to implement the Redheffer star product, we see there are some common terms. Calculating these multiple times is inefficient so we calculate them only once using intermediate parameters.
S

S11
AB 
 AB 
S
 A
S
 B
1
 
 
   
   
 S11
 S12
 I  S11 S 22  S11 S 21
A
A
B
A
B
1
 B
 AB
 B
S 21  S 21
1
 A   B 
A

I  S 22 S11  S 21
 B
 A   B  1  A   B 
S 22  S 22  S 21  I  S 22 S11  S 22 S12
Lecture 5
1
 B 
F  S21B I  S22A S11

1
A
 AB
A 
 B  A  
 B
S12
 S12
I  S11 S 22  S12
 AB
A 
 B  A  
D  S12
I  S11 S 22 
 AB 
A
 B  A 
S11
S 21
 S11
 DS11
 AB 
 B
S12
 DS12
S21AB   FS21A 
 B
S22AB   S22B  FS22A S12
Slide 48
24
4/27/2015
Using the Star Product as an Update
Very often we update our global scattering matrix using a star product.
When we use this equation as an update, we MUST pay close attention to the order that we implement the equations so that we don’t accidentally overwrite a value that we need.
S global  Si   S global 
global 
i
F  S21
global 
global 
reverse order
S22
 global 
S 21
 S22
global 
S22 
i
1
 global  
I  S22i S11


 global  
 i   global  
D  S12
I  S11 S 22 
1
 FS22S1 2
i
global 
i 
 FS 21
 global 
 global 
S12
 DS12
 global 
i 
 i
S11
 S11
 DS1 global
S 21
1
global  i 
i

F  S21 I  S22 S11

standard order


D  S12
I  S11
S global  S global  S i 
1
1
 global 
 global 
 i   global 
 S11
 DS11
S11
S 21
 global 
i
 DS12
S12
S21global  FS21global 
i 
S22global  S22i   FS22global S12
Lecture 5
Slide 49
Simplifications for TMM in LHI Media
In LHI media,
1 0 
Wi  I  

0 1 
and
Ωi  jkz ,i I
1 0 
I
 identity matrix
0 1 
Now we do not actually have to calculate  because
λ i  Ωi
Given all of this, the eigen‐vectors for the magnetic fields can be calculated as
Vi  Qi Wi λ i1  Qi Ωi1
When calculating scattering matrices, the intermediate matrices Ai and Bi are
A i  Wi1Wh  Vi1Vh  I  Vi1Vh
Bi  Wi1Wh  Vi1Vh  I  Vi1Vh
The fields and mode coefficients are now related through
P
 Px 
1  x 
cinc  Wref
P   P 
 y  y
Lecture 5
 Exref 
 ref   Wref S11cinc  S11cinc
 E y 
 Extrn 
 trn   Wtrn S 21cinc  S 21cinc
 E y 
Slide 50
25
4/27/2015
Initializing the Global Scattering Matrix
Before we iterate through all the layers, we must initialize the global scattering matrix as the scattering matrix of “nothing.”
What are the ideal properties of nothing?
1.
Transmits 100% of power with no phase change.
 global 
S12
 S21global   I
2.
Does not reflect.
 global 
S11
 S 22global  0
We therefore initialize our global scattering matrix as
0 I 
S global   

 I 0
This is NOT an identity matrix!
Look at the position of the 0’s and I’s.
Lecture 5
Slide 51
Calculating the Parameters of the Homogeneous Gaps
Our analytical solution for a homogeneous layer is
1  kx ky
Qh 

r ,h  ky2  r ,h r ,h
r ,h r ,h  kx2 


kx ky
kz2,h  r ,h r ,h  kx2  ky2
Wh  I
λ h  jkz ,h I
We are free to choose any r and
Vh  Q h Wh λ h1
r that we wish. We also wish to avoid the case of kz,h = 0. For convenience, we choose
r ,h  1.0
and
 r ,h  1  kx2  ky2
We then have
 kx ky
Qh  
2
  1  kx

Lecture 5

1  ky2 

kx ky 

Wh  I
Vh   jQ h
Slide 52
26
4/27/2015
Block Diagram of TMM Using S‐Matrices
Initialize Global Scattering Matrix
Wh  I
0 I 
S  global   

 I 0
Calculate Parameters for Layer i
Vh   jQ h
Calculate Scattering Matrix
for Layer i
kz ,i  i i  kx2  ky2
Ai  I  Vi1Vh
i i  kx2 
1  kx ky
Qi   2

i  ky  i i
kx ky 
Ωi  jkz ,i I
Vi  Qi Ωi1
Calculate Transverse Wave Vectors
kx  ninc sin  cos 
Calculate Gap Medium Parameters
Update Global Scattering Matrix
 global 
S11
  global
S 21
B i  I  Vi1Vh
Xi  e Ωi k0 Li


i 
S11
 S 22i   A i  Xi B i A i1Xi Bi
i 
S12
 S 21i   A i  Xi B i A i1Xi Bi
i
1
i
i
i
1
i
i
i
i
1
i
 global     global 
S12
S11

S 22global   S 21global
 global     i 
S12
S11

S22global  S 21i 
i 
S12

S22i  
Start
1
 global 
 global 
 global  
 i   global  
 i   global 
S11
 S11
 S12
I  S11 S 22  S11 S 21
 X B A X A  B 
 X A  B A B 
1
ky  ninc sin  sin 
i
1
 global 
 global  
 i   global  
i
S12
 S12
I  S11 S 22  S12
i
1
i  
 global 
S 21global  S 21i  I  S22globalS11
 S 21
1
i  
 global   i 
S 22global  S 22i   S 21i  I  S 22global S11
 S 22 S12
no
yes
Loop through all layers
Done?
Connect to External Regions
S  global   S  ref   S  global 
S  global   S  global   S  trn 
Finish
Calculate Source

P  pTE aˆTE  pTM aˆTM

P 1
 Px 
cinc   
 Py 
Calculate Transmitted and Reflected Fields
Calculate Longitudinal Field Components
 E xref 
 ref   S11cinc
 E y 
E zref  
 E xtrn 
 trn   S 21cinc
 E y 
E ztrn  
kx Exref  ky E yref
kˆ ref
z
kx E xtrn  ky E ytrn
k trn
z
Calculate Transmittance and Reflectance

R  Eref
2
 2
  k trn 
T  Etrn  Re  ref zref 

  trn k z 
Lecture 5
Slide 53
How to Handle Zero Layers
Follow the block diagram!!
Setup your loop this way…
NLAY = length(L);
for nlay = 1 : NLAY
...
end
For zero layers:
ER = [];
UR = [];
L = [];
If NLAY = 0, then the loop will not execute the global scattering matrix will remain as it was initialized.
0 I 
S  global   

 I 0
Lecture 5
Slide 54
27
4/27/2015
Can TMM Fail?
Yes!
The TMM can fail to give an answer and behave numerically strange any time kz = 0. This happens at a critical angle when the transmitted wave is very near its cutoff.
We fixed this problem in the gap medium, but this can also happen in any of the layers or in the transmission region.
This happens in any medium where
 r  r  kx2  ky2
Lecture 5
Slide 55
Advanced Scattering Matrices
Lecture 5
Slide 56
28
4/27/2015
Longitudinally Periodic Devices
Suppose we just calculated the scattering matrix for the unit cell of a longitudinally periodic device.
Unit Cell
S
ABC
Unit Cell 1
Unit Cell 2
Unit Cell 3
Unit Cell 4
Unit Cell 4
1
 S   S   S
A
Unit Cell 6
B
Unit Cell 7
C
Unit Cell 8
There exists a very efficient way of calculating the global scattering matrix of a longitudinally periodic device without calculating and combining all the individual scattering matrices.
S
8 
S
 S   S   S   S   S   S   S   S   S   S   S   S   S   S   S   S   S   S   S   S   S   S   S   S
8 
A
 S
B
1
C
 S
A
1
B
 S
C
1
 S
A
1
B
 S
C
1
A
 S
B
1
 S
C
1
A
 S
B
C
A
B
C
A
B
C
A
B
C
1
Both are inefficient!!!
Lecture 5
Slide 57
Cascading and Doubling
We can quickly build an overall scattering matrix that describes hundreds and thousands of unit cells.
We start by calculating the scattering matrix for a single unit cell.
S 1  S  A  S  B   S C 
A BC
Next, we keep connecting the scattering matrix to itself to keep doubling the number of unit cells it describes.
S 2  S 1  S 1
S   4   S   2   S  2 
S   8   S  4   S   4 
Lecture 5
Slide 58
29
4/27/2015
Algorithm for Arbitrary Number of Unit Cells
Step 1 – Calculate the scattering matrix for the unit cell.
S 1  S  A  S  B   S  C 
A BC
Step 2 – Determine the binary digits for the total number of repetitions.
Chain of 10 unit cells  1010
Step 3 – Perform cascading a doubling while updating the global scattering matrix only with the scattering matrices corresponding to binary digits of ‘1’.
1. Initialize Algorithm
S  N 
0 I 


 I 0
S  bin   S1
2. Perform modified cascading and doubling
a. If binary digit is ‘1’  S N   S  N   S  bin 
Loop # binary digits
b. S bin   S bin   S bin 
c. repeat through all binary digits
Lecture 5
Slide 59
Block Diagram for Modified Cascading and Doubling Algorithm
Example
Inputs
S
1
 S-matrix of one unit cell
N  Number of times to repeat unit cell
Convert N to binary
10110
N = 22  10110
Initialize Algorithm
1: S(N) now 2 unit cells
S(bin) now 4 unit cells
0 I 
S  N   

 I 0
Loop through all binary digits starting with the least significant digit.
Lecture 5
0: No update to S(N)
S(bin) now 2 unit cells
S bin   S 1
1: S(N) now 6 unit cells
S(bin) now 8 unit cells
Output
Done?
S  N 
no
digit = 1?
yes
0: No update to S(N)
S(bin) now 16 unit cells
1: S(N) now 22 unit cells
S(bin) now 32 unit cells
Update S(N)
S N   S N   S bin 
no
Update Doubling
S bin   S  bin   S bin 
Slide 60
30
4/27/2015
Dispersion Analysis (1 of 2)
An overall scattering matrix is calculated that describes the unit cell.
 uc 
 c0  S11
      uc 
c N 1  S11
 uc  
S11
 c0 

 uc  c  
S11
  N 1 
S uc  same as S1 on previous slide
The terms are rearranged in “almost” the form of a transfer matrix.
 uc   
 uc 
0 S12
c N 1  S11
  uc


 uc    
 
 I S 22  c N 1  S 21
I  c0 
 
0  c0 
If the device is infinitely periodic in the z direction, then the following periodic boundary condition must hold.
cN 1 
jk z  z
  e
c N 1 
c0 
 
c0 
Here kz is the “effective” propagation constant of the mode.
Lecture 5
Slide 61
Dispersion Analysis (1 of 2)
We substitute the periodic boundary condition into our rearranged equation to get
 uc 
S11
  uc 
S 21
I  c0 
jk 
   e z z
0  c0 
 uc   
0 S12
c 0 




 uc 
 I S 22  c0 
This is a generalized eigen‐value problem.
Ax   Bx
S uc 
A   11uc
 
S 21
I 

0 
c  
x   0 
c 0 
 uc  
0 S12
B

 uc 
 I S 22 
  e jk 
z
z
[V,D] = eig(A,B);
Eigen vectors
Lecture 5
Eigen values
Slide 62
31
4/27/2015
Alternatives to Scattering Matrices
Lecture 5
Slide 63
Transmittance Matrices (T‐Matrices)
The T‐matrix method is the transfer matrix method where forward and backward waves are distinguished.
left

 c trn   T11 T12  cinc

c right  T T   
22   c ref 
 inc   21
Benefits
• Much faster (5 to 10 times)
• Unconditionally stable
Drawbacks
• Less memory efficient
• Cannot exploit longitudinal periodicity
• Less popular in the literature
M. G. Moharam, Drew A. Pommet, Eric B. Grann, “Stable implementation of the rigorous coupled-wave analysis for surfacerelief gratings: enhanced transmittance matrix approach,” J. Opt. Soc. Am. A, Vol. 12, No. 5, pp. 1077-1086, 1995.
Lecture 5
Slide 64
32
4/27/2015
Hybrid Matrices (H‐Matrices)
The h‐matrix method is borrowed from electrical two‐port networks.
V1   h11
 I   h
 2   21
h12   I1 
h22  V2 
h11 
V1
I2
I
h21  2
I1
V2  0
V1
V2
I1  0
V2  0
I
h22  2
V2
I1  0
h12 
In the framework of fields, the h‐matrix is defined as
 Ex ,i 1 
 E   i 
 y ,i 1    H11
 H x ,i   H i 
    21
 H y ,i 
 H x ,i 1 

i    
H12
H y ,i 1 


H 22i    Ex ,i 


 E y ,i 
Claimed Benefits
• Improved numerical stability
• More concise formulation
• Simpler to implement
• Improved numerical efficiency (30% better than ETM)
• Unconditionally stable
Eng L. Tan, “Hybrid-matrix algorithm for rigorous coupled-wave analysis of
multilayered diffraction gratings,” J. Mod. Opt., Vol. 53, No. 4, pp. 417-428, 2006.
Lecture 5
Slide 65
R‐Matrices
The R‐matrix method is essentially the impedance matrix framework borrowed from electrical two‐port networks.
V1   z11
V    z
 2   21
z12   I1 
z22   I 2 
z11 
z21 
V1
I1
V2
I1
V1
I2
z12 
I2 0
z22 
I2 0
V2
I2
I1  0
I1  0
In the framework of fields, the h‐matrix is defined as
 Ex ,i 1 
 E   i 
 y ,i 1    R11
 E x ,i   R  i 

  21
 E y ,i 
 H x ,i 1 

i    
R12
H y ,i 1 


R 22i    H x ,i 
  
 H y ,i 
Claimed Benefits
• Unconditionally stable
• Improved numerical efficiency
Lifeng Li, “Bremmer series, R-matrix propagation algorithm, and numerical modeling of
diffraction gratings,” J. Opt. Soc. Am. A, Vol. 11, No. 11, pp. 2829-2836, 1994.
Lecture 5
Slide 66
33