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 kx x k0 ky ky k0 Eliminate longitudinal components Ez and Hz by substitution. k kz z k0 dE y jky Ez r H x dz dEx jkx Ez r H y dz jkx E y jky Ex r H z dH y jky H z r Ex dz dH x jkx H z r E y dz jkx H y jky H x r Ez dEx kx ky k 2 H x r x H y dz r r dE y ky2 kx ky Hy r H x dz r r dH x kx ky k 2 Ex r x E y dz r r dH y ky2 kx ky E r Ex dz r r y Slide 4 2 4/27/2015 Derivation of Two 22 Matrix Equations We can write our two sets of two equations in matrix form as dEx kx ky k 2 H x r x H y dz r r dE y ky2 k k r H x x y H y r dz r d Ex 1 dz E y r dH x kx ky k 2 Ex r x E y dz r r 2 dH y k y k x ky E r Ex dz r r y d H x 1 kx ky dz H y r ky2 r r kx ky 2 ky r r r r kx2 H x kx ky H y r r kx2 Ex kx ky 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 kx ky 2 ky r r 1 kx ky r ky2 r r d H x 1 kx ky dz H y r ky2 r r r r kx2 H x kx ky H y r r kx2 kx ky H x d Ex P dz E y H y Q r r kx2 Ex 1 kx ky r ky2 r r kx ky E y r r kx2 kx ky 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‐jz 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 kx2 ky2 kz2 Using this relation, we can simplify the matrix equation for 2. Ω 2 PQ 1 kx ky r r ky2 r r r r kx2 kx ky kx ky 2 k y r r r r kx2 kz2 kx ky 0 0 2 kz I kz2 1 0 I identity matrix 0 1 Since this is a diagonal matrix, we can conclude that 1 0 WI 0 1 λ 2 Ω2 jk λ z 0 0 jkz I jkz e jkz 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 λzc Vλe λzc Q We λzc We λzc QWe λzc QWeλzc 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 kx ky jkx yz zy kx zx yx yz zx j ky 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 jky xz zx zz zz zz zz zz zz Sort Eigen‐Modes E x Ey H yz zy x jkx H y zz zz xz zy j kx ky zz zz kˆx2 k yy yz zy zz zz kx ky 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 kx2 1 kx ky 2 r k y r r kx ky r r kx2 1 kx ky Q 2 r k y r r kx ky 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 c2 c2 c1 ci Li ψ i z field within i th layer Lecture 5 z ci mode coefficients inside i th layer ci 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 c2 V2 c2 Note: k0 has been incorporated to normalize Li. Lecture 5 Slide 19 Definition of A Scattering Matrix c1 S11 S12 c1 c2 S 21 S 22 c2 S 21 c 1 S11 c1 Lecture 5 S11 reflection S 21 transmission transmission reflection c2 S 22 S12 This is consistent with experimental convention. c2 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 c2 V2 c2 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 Wi1W j Vi1V j Bij A ij Bij Wi1W j Vi1V 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 c2 Ai 2 c2 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 c2 ? ? c2 Lecture 5 Slide 21 The Scattering Matrix The scattering matrix Si of the ith layer is defined as: c1 i c1 S c2 c2 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 i21Xi Bi1 S11 i 12 i1 Xi Bi 2 A i21Xi Bi1 i 21 i2 Xi Bi1A i11Xi 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 Wi1W j Vi1V j Bi 2 Bi1 i 2 Bi 2 A ij Wi1W j Vi1V 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 ci1 ci ci1 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 c2 c2 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 c2 c2 S i S i 11i S 21 i S12 S 22i But the elements are calculated as S A i 12 i i S21i S12 i i S 22 S11 Lecture 5 Xi Bi A i1 Si X B A X A B X B X A B A B i S11 A i Xi Bi A i1Xi 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 Wi1Wh Vi1Vh Bi Wi1Wh Vi1Vh 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 S1 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 S21 S AB S A S B A S12 S22A S B S B 11B S21 B S12 S22B The combined scattering matrix is then S11 AB S AB AB S11 AB S21 AB S12 S22AB 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 S21 S21 I S22 S11 S 21 AB B A B A 1 B A B S22AB S22B S21B 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 c2 S 21 A S12 c1 A S22 c2 B c2 S11 B c3 S 21 B S12 c2 B S 22 c3 We expand these into four matrix equations. A A B B Eq. 1 Eq. 3 c1 S11 c1 S12 c2 c2 S11 c2 S12 c3 A A B B Eq. 2 Eq. 4 c2 S 21 c1 S 22 c2 c3 S 21 c2 S 22 c3 c2 We substitute Eq. (2) into Eq. (3) to get an equation with only . c2 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 c3 Eq. 5 A 22 B 11 2 B S 21A c1 S 22A S12 c3 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 c3 ? ? c3 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 Wh1Wref Vh1Vref 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 L0 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 S21 2A trn1 trn S22trn A trn1 B trn s trn A trn Wh1Wtrn Vh1Vtrn 1 h r ,h r ,II r ,h r ,II 1 h B trn W Wtrn V Vtrn lim L0 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 S21globalcinc 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 kx E xref ky E yref k ref z Eztrn kx Extrn ky 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 jkz 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 S21B I S22A 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 S21AB FS21A B S22AB S22B FS22A 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 Si S global global i F S21 global global reverse order S22 global S 21 S22 global S22 i 1 global I S22i S11 global i global D S12 I S11 S 22 1 FS22S1 2 i global i FS 21 global global S12 DS12 global i i S11 S11 DS1 global S 21 1 global i i F S21 I S22 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 S21global FS21global i S22global S22i FS22global S12 Lecture 5 Slide 49 Simplifications for TMM in LHI Media In LHI media, 1 0 Wi I 0 1 and Ωi jkz ,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 λ i1 Qi Ωi1 When calculating scattering matrices, the intermediate matrices Ai and Bi are A i Wi1Wh Vi1Vh I Vi1Vh Bi Wi1Wh Vi1Vh I Vi1Vh 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 S21global 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 kx ky Qh r ,h ky2 r ,h r ,h r ,h r ,h kx2 kx ky kz2,h r ,h r ,h kx2 ky2 Wh I λ h jkz ,h I We are free to choose any r and Vh Q h Wh λ h1 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 kx2 ky2 We then have kx ky Qh 2 1 kx Lecture 5 1 ky2 kx ky 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 kz ,i i i kx2 ky2 Ai I Vi1Vh i i kx2 1 kx ky Qi 2 i ky i i kx ky Ωi jkz ,i I Vi Qi Ωi1 Calculate Transverse Wave Vectors kx ninc sin cos Calculate Gap Medium Parameters Update Global Scattering Matrix global S11 global S 21 B i I Vi1Vh Xi e Ωi k0 Li i S11 S 22i A i Xi B i A i1Xi Bi i S12 S 21i A i Xi B i A i1Xi 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 S22global S 21i i S12 S22i 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 ky 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 S22globalS11 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 kx Exref ky E yref kˆ ref z kx E xtrn ky 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 kx2 ky2 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 S1 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 S1 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. cN 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
© Copyright 2024