1. technology and computing
2. consumer electronics
3. tv and video equipment
4. projectors

5/11/2015
Instructor
Dr. Raymond Rumpf
(915) 747‐6958
[email protected]
EE 4395/5390 – Special Topics
Computational Electromagnetics
Lecture #4
Transfer Matrix Method
 These notes may contain copyrighted material obtained under fair use rules. Distribution of these materials is strictly prohibited 
Lecture 4
Slide 1
Outline
• Maxwell’s equations for 1D structures
• Solution to Maxwell’s equations in a homogeneous layer
• Multilayer structures
• TMM is an inherently unstable method
Lecture 4
Slide 2
1
5/11/2015
Maxwell’s Equations for 1D Structures
Lecture 4
Slide 3
1D Structures
Sometimes it is possible to describe a physical device using just one dimension. Doing so dramatically reduces the numerical complexity of the problem and is ALWAYS GOOD PRACTICE.
x
Region I
Reflection Region
y
z
Region II
Transmission Region
Lecture 4
Slide 4
2
5/11/2015
3D  1D Using Homogenization
Many times it is possible to approximate a 3D device in one dimension. It is very good practice to at least perform the initial simulations in 1D and only moving to 3D to verify the final design.
Physical Device
Effective Medium Approximation in 1D
Lecture 4
Slide 5
3D  1D Using Circuit‐Wave Equivalence
d2
d1
d3
d5
d4
Z trn  Z L
Z ref



N   r r
Lecture 4
d6
Z1
Z2
Z3
Z4
Z5
Z6
in
1
2
3
4
5
6
L
N in
N1
N2
N3
N4
N5
N6
NL
d1
d2
d3
d4
d5
d6
Slide 6
3
5/11/2015
Starting Point
We start with Maxwell’s equations in the following form. Here we have assumed isotropic materials.
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

 k 0 r E x
y
z
H x H z

 k 0 r E y
z
x
H y H x

 k 0 r E z
x
y
Lecture 4
Slide 7
Calculation of the Wave Vector Components
The components kx and ky are determined by the incident wave and are continuous throughout the 1D device. The kz component is different in each layer and calculated from the dispersion relation in that layer.
k x  k0 r ,inc r ,inc sin  cos 
k y  k0 r ,inc r ,inc sin  sin 
k z ,i  k02 r ,i r ,i  k x2  k y2
i  Layer #
Lecture 4
Slide 8
4
5/11/2015
Waves in Homogeneous Media
A wave propagating in a homogeneous layer is a plane wave. It has the following mathematical form.
 
  
jk y
E  r   E0 e jk r  E0 e jk x x e y e jk z z
 
  
jk y
H  r   H 0 e jk r  H 0 e jk x x e y e jk z z
Note: e+jkz sign convention was used for propagation in +z direction.
When we take derivatives of these solutions, we see that
 jk y
 
  
  jk x x jk y y jkz z
E r  
E0 e e e
 jk x E0 e y e jkz z e jkx x  jk x E  r 
x
x


 jk x
x
 jk y
 
  
  jk x x jk y y jk z z
E r  
E0 e e e
 jk y E0 e y e jk z z e jk x x  jk y E  r 
y
y


 jk y
y
We cannot say that because the structure is not  z  jk z
homogeneous in the z direction.

 jk z
z




Lecture 4
Slide 9
Reduction of Maxwell’s Eqs. to 1D
Given that

 jk x
x

 jk y
y
Maxwell’s equations become
jk y Ez 
dE y
 k0 r H x
dz
dEx
 jk x Ez  k0 r H y
dz
jk x E y  jk y Ex  k0 r H z
dH y
jk y H z 
 k 0 r E x
dz
dH x
 jk x H z  k0 r E y
dz
jk x H y  jk y H x  k0 r Ez
Note: z is the only independent variable left so its derivative can be ordinary.
Lecture 4
Slide 10
5
5/11/2015
Normalize the Parameters
We normalize the parameters according to
z   k0 z
k
kx  x
k0
ky
ky 
k0
k
kz  z
k0
Using the normalized parameters, Maxwell’s equations become
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
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
Lecture 4
Slide 11
Solve for the Longitudinal Components Ez and Hz
We solve the third and sixth equations for the longitudinal field components Hz and Ez.
dE
jky Ez  y   r H x
dz 
dEx
 jkx Ez   r H y
dz 
jkx E y  jky Ex   r H z
j 
 H z 
k x E y  ky Ex




r
dH y
  r Ex
jky H z 
dz 
dH x
 jkx H z   r E y
dz 
jkx H y  jky H x   r Ez
Lecture 4
 Ez 
j  
k x H y  ky H x
r
Slide 12
6
5/11/2015
Eliminate the Longitudinal Components
We eliminate the longitudinal field terms by substituting them back into the remaining equations.
dE
jky Ez  y  r H x
dz 
dEx
 jkx Ez  r H y
dz 
j 
H z 
k x E y  ky Ex
r


dH y
jky H z 
  r Ex
dz 
dH x
 jkx H z   r E y
dz 
j  
Ez 
k x H y  ky H x
r

dE y
  r  r H x
ky2 H x  kx ky H y   r
dz 
dE
 r x  kx2 H y  kx ky H x  r  r H y
dz 

dH y
ky2 Ex  kx ky E y   r
 r  r Ex
dz 
dH x  2
 k x E y  kx ky Ex   r  r E y
r
dz 
Lecture 4
Slide 13
Rearrange Maxwell’s Equations
Here we simply change the order that we our previous equations.
dE y
  r  r H x
ky2 H x  kx ky H y   r
dz 
dE
 r x  kx2 H y  kx ky H x  r  r H y
dz 
dH y
  r  r Ex
ky2 Ex  kx ky E y   r
dz 
dH x  2
r
 k x E y  kx ky Ex   r  r E y
dz 
Lecture 4

dEx kx ky
k 2 

H x   r  x  H y
r
r 
dz 


dE y  ky2
kx ky
Hy
   r  H x 

r
dz    r


dH x kx ky
k 2 
Ex    r  x  E y

r
r 
dz 


dH y  ky2
kx ky
E
    r  Ex 

r y
dz   r

Slide 14
7
5/11/2015
Matrix Form of Maxwell’s Equations
The remaining four equations can be written in matrix form as
dE x 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

dz    r
r





dH x k x k y
k2 
E    x  E

dz 
r x  r r  y

dH y  ky2
k k
    r  Ex  x y E y

dz   r
r


 0


 Ex  
0
 
d  Ey  

dz   H x   kx ky
  
 H y   r
 ky2
  r
  r
kx ky
0
r
k
0
r 

2
y
r
kx2
r
kx ky
r
Lecture 4
 r
0
0
kx2 

r 
kx ky   Ex 


 r   E y 
 
 Hx 
0   H 
 y
 

0 

r 
Slide 15
BTW…for Anisotropic Materials
   yz
 
 kx zx 
  j  ky
 zz 
   zz

 Ex   jk   xz   zx 

y
  
  zz  zz 
  Ey  

z   H x   kx ky
 yz zx 
   
  yx 

 zz 
 H y    zz

  k 2
  
  y   xx  xz zx 
 zz 
   zz
  yz  zy 
jkx 


  zz  zz 
 zy 
 
 j  kx xz  ky


 zz 
zz

 yz zy 
 kx2
  yy 



 zz 
zz

 kx ky
 xz  zy 
  xy 


  zz
 zz 

 kx ky
 yz  zx 
  yx 



 zz 
  zz
 ky2
  
  xx  xz zx 

  zz
 zz 

  yz   zx 
 j  ky
 kx

 zz 
  zz

 
jky  xz  zx 
  zz  zz 
 kˆx2
 yz  zy  
  yy 
 

 zz  
  zz

 kx ky
 xz  zy    Ex 
  xy 

  
  zz
 zz   E y

 
  H 
x
  yz  zy 


 

jk x 

H y 





 zz
zz 

   xz   zy  
 j  kx
 ky
 
 zz  
  zz
e j z
Note: This is for the sign convention.
Lecture 4
Slide 16
8
5/11/2015
Solution to Maxwell’s Equations in a Homogeneous Layer
Lecture 4
Slide 17
Matrix Differential Equation
Maxwell’s equations can now be written as a single matrix differential equation.
dψ
 Ωψ  0
dz 
 Ex  z   


E y  z 
ψ  z   
 H x  z  


 H y  z   
Lecture 4

 0


 0

Ω
 kx ky
 
r

 ky2
  r
 r
kx ky
0
r
k
0
r 

2
y
r
kx2
r
kx ky
r
 r
0
0
kx2 

r 
k k 
 x y 
r 


0 


0 

r 
Slide 18
9
5/11/2015
Solution of the Differential Equation (1 of 3)
The matrix differential equation is
dψ
 Ωψ  0
dz
This is actually a set of four coupled differential equations. The system of four equations can be solved as a single matrix equation as follows.
ψ  z    eΩz ψ  0 
This is easy to write, but how do we compute the exponential of a matrix?
Lecture 4
Slide 19
Functions of Matrices
It is sometimes necessary to evaluate the function of a matrix.
f A  ?
It is NOT correct to calculate the function of every element in the matrix A individually. A different technique must be used.
To do this, we first calculate the eigen‐vectors and eigen‐values of the matrix A.
f  A  ?
A

W eigen-vector matrix of A
λ eigen-value matrix of A
1 0
0 
2
λ
0 0

0 0
0
0
3
0
0
0 
0

4 
The function of the matrix is then evaluated as
f  A   W  f  λ   W 1
Lecture 4
This is very easy to evaluate because  is a diagonal matrix so the function only has to be performed individually on the diagonal elements.
Slide 20
10
5/11/2015
Solution of the Differential Equation (1 of 2)
We had the following matrix differential equation and general solution
dψ
 Ωψ  0
dz

ψ  z   eΩz ψ  0 
We can now evaluate the matrix exponential using the eigen‐values and eigen‐vectors of the matrix .
Ω

W eigen-vector matrix
λ eigen-value matrix
eλz
eΩz  Weλz W 1
e1z

0

 0

 0
0
e
0
0
2 z 
0
0
e3 z
0





e4 z 
0
0
0
Lecture 4
Slide 21
Solution of the Differential Equation (2 of 2)
The solution to the matrix differential equation is therefore
dψ
 Ωψ  0
dz
ψ  z    eΩz ψ  0 

ψ  z    We λz W 1ψ  0 

c
We can combine the unknown initial values (0) with W-1 because that product just leads to another column vector of unknown constants.
Our final solution is then
dψ
 Ωψ  0
dz 
Lecture 4

ψ  z    We λzc
c  W 1ψ  0 
Slide 22
11
5/11/2015
Interpretation of the Solution
ψ  z    Weλzc
(z’) – Overall solution which is the sum of all the modes at plane z’.
W – Square matrix who’s column vectors describe the “modes” that can exist in the material. These are essentially pictures of the modes which quantify the relative amplitudes of Ex, Ey, Hx, and Hy. c – Column vector containing the amplitude coefficient of each of the modes. This quantifies how much energy is in each mode.
ez’ – Diagonal matrix describing how the modes propagate. This includes accumulation of phase as well as decaying (loss) or growing (gain) amplitude.
Lecture 4
Slide 23
Getting a Feel for the Numbers (1 of 2)
For a layer with r=9.0 and r=1.0 (i.e. n=3.0) and a wave at normal incidence, we will have
0
0
Ω
0

 9
0 0 1
0 1 0 
9 0 0

0 0 0
This has the following eigen‐vectors and eigen‐values.
  j 0.32
 0
W
 0

 0.95
Lecture 4
j 0.32
0
0
0.95
0
0 
j 0.32  j 0.32 
0.95
0.95 

0
0 
0
 j 3.0
 0
 j 3.0
λ
 0
0

0
 0
0
0
j 3.0
0





 j 3.0 
0
0
0
Slide 24
12
5/11/2015
Getting a Feel for the Numbers (2 of 2)
We see that the modes occur as either an Ex‐Hy or Ey‐Hx pair. This is consistent with plane waves. Due to the normalization, they are 90° out of phase. A sign difference indicates forward and backward waves. Only the relative amplitude difference between E and H is important here.
  j 0.32
 0
W
 0

 0.95
j 0.32
0
0
0.95
0
0 
j 0.32  j 0.32 
0.95
0.95 

0
0 
r
E
   0
r
H
E

H
r 1

r 3
The modes in W only contain information about the relative amplitudes of the field components.
We know the refractive index (n = 3.0), so the eigen‐values are consistent with what we would expect. The signs correspond to forward and backward waves.
0
 j 3.0
 0
 j 3.0
λ
 0
0

0
 0
0 
0 
j 3.0
0 

 j 3.0 
0
0
0
e z  e jn cosinc z
  jn cos inc
n  r  r  3
The numbers in  describe how the modes accumulate phase in the z direction. This is essentially just the refractive index of the material.
Lecture 4
Slide 25
Visualizing the Modes
  j 0.32
 0
W
 0

 0.95
0.95
0
0 
j 0.32  j 0.32 
0.95
0.95 

0
0 
0
 j 3.0
 0
j
3.0

λ
 0
0

0
 0
0 
0
0 
j3.0
0 

0
 j 3.0
j 0.32
0
0
Mode 1
0
Mode 2
0.95
0.95
j0.32
Mode 4
Mode 3
Mode 2
Mode 1
-j0.32
j0.32
Mode 3
Mode 4
0.95
0.95
-j0.32
Lecture 4
Slide 26
13
5/11/2015
Multilayer Structures
Lecture 4
Slide 27
Geometry of an Intermediate Layer
+z
Layer i-1
ψ i 1  k0 Li 1 
Layer i
ψ i  zi 
ci 1
ψ i  k0 Li 
ψ i1  0 
ψi  0
Li 1
Layer i+1
Li
ci
Li 1
ci 1
zi is a local z‐coordinate inside the ith layer that starts at zero at the layer’s left side.
Lecture 4
Slide 28
14
5/11/2015
Field Relations
Field inside the ith layer:
 Ex ,i  zi  


E y ,i  zi  
 Wi e λi zi ci
ψ i  zi    
 H x ,i  zi  


 H y ,i  zi  
Boundary conditions at the first interface:
ψ i 1  k0 Li 1   ψ i  0 
Wi 1eλi1k0 Li1 ci 1  Wi ci
We need to include k0 in the exponential to normalize Li-1 because The parameter i-1 expects to multiply a normalized coordinate.
Boundary conditions at the second interface:
ψ i  k0 Li   ψ i 1  0 
Wi eλi k0 Li ci  Wi 1ci 1
Note: We must equate the field  on either side of the interfaces, not the mode coefficients c.
Lecture 4
Slide 29
The Transfer Matrix
The transfer matrix Ti of the ith
layer is defined as:
ci 1  Ti  ci
Ti
After some algebra, the transfer matrix is computed as
Ti  Wi11Wi eλi k0 Li
Lecture 4
Slide 30
15
5/11/2015
The Transfer Matrix Method
The transfer matrix method (TMM) consists of working through the device one layer at a time and calculating an overall (global) transfer matrix.
T1
T2
T3
T4
T5
Transmission
Region
Reflection
Region
  T5  T4  T3  T2  T1
Tglobal
This is standard matrix multiplication.
The order of multiplication may seem backwards here, but it is not. Recall the definition of the transfer matrix to have this make sense.
Lecture 4
Slide 31
The Global Transfer Matrix
The transfer matrix so far is not yet the “true” global transfer matrix because it does not connect the reflection region to the transmission region. It only connects the amplitude coefficients of Layer 1 to the amplitude coefficients in the transmission region. This is a result of how we defined the transfer matrix.
 c1
c trn  Tglobal
The global transfer matrix must connect the amplitude coefficients in the reflection region to the amplitude coefficients in the transmission region. Boundary conditions at the first interface require
Wref c ref  W1c1
Solving this for c1 yields
c1  W11Wref c ref
The global transfer matrix is derived by substituting this result into the first equation.
 W11Wref c ref
c trn  Tglobal
 W11Wref
Tglobal  Tglobal
Lecture 4
Tglobal  T5  T4  T3  T2  T1  W11Wref
Slide 32
16
5/11/2015
TMM is an Inherently Unstable Method
Lecture 4
Slide 33
The Multi‐Layer Problem
The figure below is focused on an arbitrary layer in a stack of multiple layers. We will be examining the wave solutions in the ith layer.
Lecture 4
Slide 34
17
5/11/2015
Wave Solutions in ith Layer
Recall that the wave vector can be purely real (pure oscillation), purely imaginary (pure exponential decay), or complex (decaying oscillation).
k  k
k  k   jk 
k  jk 
Lecture 4
Slide 35
Backward Waves in ith Layer
Due to reflections at the interfaces, there will also be backward traveling waves in each of the layers. These can also have wave vectors that are real, imaginary or complex, so they can oscillate, decay/grow, or both.
Lecture 4
Slide 36
18
5/11/2015
All Waves are Treated as Forward Waves
The true transfer matrix method treats all waves as if they are forward propagating. Decaying fields associated with backward waves become exponentially growing fields and quickly become numerically unstable.
 Ex ,i  zi  


E y ,i  zi  
 Wi eλi zi ci
ψ i  zi    
 H x ,i  zi  


 H y ,i  zi  
Lecture 4
Slide 37
TMM is Inherently Unstable
Our wave solution was
 Ex  z   


E y  z 
 We λzc
ψ  z   
 H x  z  


 H y  z   
This treats all energy as forward propagating.
We know that backward waves exist. We also know that decaying fields exist when a wave is evanescent or propagating in a lossy
material.
When backward waves are decaying and treated as forward propagating waves, they grow exponentially. This leads to numerically instability.
The TMM is inherently an unstable method because it treats everything as forward propagating.
Lecture 4
Slide 38
19
5/11/2015
The Fix
We are treating all energy as forward propagating because we did not distinguish between forward and backward waves.
Clearly, the first part of the fix is to distinguish between forward and backward propagating waves.
This can be accomplished by calculating the Poynting vector associated with the modes and looking at the sign of the z component. Be careful! We are using a normalized magnetic field.
  
 E  H
z  Ex H y  E y H x
 H y 
 H 
z  Ex  j
 Ey  j x 
  
0 
 0 

j
z   Ex H y  E y H x 
0
0 
 i 0.32 i 0.32
 0

i
i
0
0.32
0.32


W
 0
0
0.95 0.95 


0.95
0
0 
 0.95
0
Lecture 4
Slide 39
Rearrange Eigen Modes
Now that we know which eigen‐modes are forward and backward propagating, we can rearrange the eigen‐vector and eigen‐value matrices to group them together.
0
0 
 i 0.32 i 0.32
 0

i

i
0
0.32
0.32

W
 0
0
0.95 0.95 


0.95
0
0 
 0.95
i 0.32
0
0 
 i 0.32
 0

i
i
0.32
0
0.32


W  
 0
0.95
0
0.95 


0
0.95
0 
 0.95
0
0
0 
i3.0
 0 i3.0 0
0 
λ
 0
0
i3.0
0 


0
0 i3.0 
 0
0
0 
i3.0 0
 0 i3.0
0
0 
λ  
 0
0 i3.0
0 


0
0
i3.0 
 0
You will also need to adjust the vertical positions of the eigen‐values so that ’ remains a diagonal matrix.
Lecture 4
Slide 40
20
5/11/2015
New Interpretation of the Matrices
0
i 0.32
0 
 i 0.32
 0

i
0.32
0

i
0.32

W  
 0
0.95
0
0.95 


0
0.95
0 
 0.95
Ex
Ey
H x
H
y
0
0 
i3.0 0
 0 i3.0
0
0 
λ  
 0
0 i3.0
0 


0
0
i3.0 
 0
We have now partitioned our matrices into forward and backward propagating elements.
 W
W   E
 WH
e λ z

 0

e
Note: For anisotropic materials, all the eigen‐vectors and eigen‐
values are in general unique.
λz 
i3.0 0 
λ  

 0 i3.0 
WE 

WH 
Ex
Ey
H x
H
y
0 

 λ  z

e
0 
 i3.0
λ  
i3.0 
 0
Lecture 4
Slide 41
Revised Solution to Differential Equation
The matrix differential equation and its original solution was
dψ
 Ωψ  0
dz

ψ  z   We λzc
After distinguishing between forward and backward propagating waves and grouping them in the matrices, we can write our solution as
 WE

ψz    
 WH
WE  e

WH   0
λ  z
0  c  
 

e  λ z  c 
We now have separate mode coefficients c+ and c- for forward and backward propagating modes, respectively.
We will pick up here next lecture.
Lecture 4
Slide 42
21