CEM计算电磁学 -- Lecture 2 学习笔记 (1) ---TMM 传输矩阵法(1)

一、电磁场中的一维结构

将三维结构降为一维结构：

二、传递矩阵法

1、模型结构

模型结构的建立是基于上一节的内容，也就是通过将三维结构降为一维结构来进行分析的，类似于切片的方法。如下图所示，假设有一束光打在由四层介质板构成的器件上。
将Layer1最靠前一侧的界面0（z=0处）的电磁场定义为E0和H0，将Layer1靠后一侧的界面1的电磁场定义为E1和H1。当发射电磁波时，对于这个器件而言，我们假设电磁波是从界面0出发，经过Layer1后，到达界面1的，因此若要建立起两个界面电磁场之间的关系，关键就在于对Layer1层的分析。

而传输矩阵（Transfer Matrix，通常称为T矩阵），就是我们用来构建起两个界面电磁场之间关系的一个矩阵，即E0和H0经过T矩阵的变换后变为了E1和H1，也由此我们可以认为：事实上，T矩阵是表征了Layer1自身的一些属性（如：电导率、磁导率等）对电磁波产生的影响。
基于此，我们可以用一个简单的公式来表示两个界面电磁场间的关系，即：[E1, H1]=T1 • [E0, H0]，T为矩阵；同样地，对于界面2的电场E2和界面1的电场E1，我们也可以获得关系式：[E2, H2]=T2 • [E1, H1]=T2 • T1 • [E0, H0]，即如上图所示。
按照以上的方法，对整个器件涉及到的界面关系均进行推导，我们即可以获得输入界面和输出界面之间的电磁场关系，如下图所示。对于该关系中所有涉及到的T矩阵的累乘，我们将之称为全局传输矩阵（Global Transfer Matrix）。

此时问题的关键，就变为了要如何推导出T矩阵？ 这里，我们先给出结果（下图），再对过程进行分析。由下图可以看出，有两种方法可以获取最终的解，一种方法是需要 sort（分类）的，在这个方法过程中需要使用到4×4矩阵；另一个方法是不需要 sort 的，其使用到的是2×2矩阵。

2、4×4矩阵方法（需要sort）

（1）前提内容

在均匀介质中，假设电磁波的入射方向是 +z 方向，则：
同理，在均匀介质中，假设电磁波的入射方向是 k⃗\vec{k}k 方向，则：

对于器件模型而言，z轴通常表示了多层的结构，x,y方向则表示的是面。通常来说，对于z轴某一点所在面，其沿x和y方向延伸出去的部分都是均匀的，而沿z轴方向延伸出去的部分则是非均匀的（因为层与层的性质不同）。因此，我们可以对x和y方向的电场进行求导，如下图所示：

但由于沿z方向结构的非均匀性，故而其不能直接进行求导：

（2）4×4矩阵方程式（用于求解电磁场方程）

在 Lecture 1 学习笔记 (1) 中，通过对麦克斯韦方程组的变形和展开，得到了一组关于 E 和 H 在三个方向的等式，它们是通常用于计算电磁学中的麦克斯韦方程组，我们将之称为Starting Point for CEM，即：
需要注意的是：这些等式对应的是LHI的情况（ [ε\varepsilonε]=ε\varepsilonε，[μ\muμ]=μ\muμ ）

利用之前推导的公消去关于Ex、Ey和Hx、Hy的偏导，由于此时只剩下z变量，因此可将偏分改为微分：

将以上方程进行归一化处理，并对各自的第三个式子进行变形以得到Ez和Hz：

将Ez和Hz代入前两个式子，此处之所以要消去Ez和Hz，是因为我们将三维问题降为了一维，因此我们现在的处理对象是切片后的面，即x和y轴构成的面。然后再对式子进行化简，即可得到Exy和Hxy关于z‘变化的方程：

将之整理成矩阵形式（需要注意的是：这些等式对应的是LHI的情况（ [ε\varepsilonε]=ε\varepsilonε，[μ\muμ]=μ\muμ ）），即：

至次，我们获得了一个4×4的矩阵方程。

（3）LHI 情况下的解

将上式的表达形式进行简化：

此时，我们要做的就是求出Ψ\PsiΨ的解，即可获得该器件关于z′z'z′的电磁场方程。该方程是一个一阶微分方程，其解的形式为：
其中eΩz′e^{^{\Omega z'}}eΩz′为矩阵，因此可利用矩阵的性质来对其计算进行简化。
根据矩阵的运算性质，对于一个矩阵函数 f(A)f(A)f(A) 而言，A为矩阵，用V表示A的特征矢量矩阵，用D表示A的特征值矩阵，则f(A)=V⋅f(D)⋅V−1f(A)=V\cdot f(D) \cdot V^{^{-1}}f(A)=V⋅f(D)⋅V−1，其中f(D)f(D)f(D)为关于f(D1)f(D1)f(D1)，f(D2)f(D2)f(D2)…的对角线矩阵。（ Matlab中的语句就是：[V,D]=eig(A) ）
因此，对于f(Ω)=eΩz′f(\Omega)=e^{^{\Omega z'}}f(Ω)=eΩz′来说，求出其特征矢量矩阵W，和特征值矩阵λ\lambdaλ，则可得：

将该式代入Ψ\PsiΨ的通解中，获得关于Ψ\PsiΨ的解，其中，用ccc来表示W−1Ψ(0)W^{^{-1}} \Psi(0)W−1Ψ(0)，即：

至此，我们获得了电磁场关于z′z'z′的表达式。

（4）计算传输矩阵T

在获取到电场和磁场的表达式后，需要引入边界条件来对器件边界情况进行一个约束。
根据器件建立的模型（如下图所示），以Layer i 为例，其存在左右两个界面，因此有两个边界条件，再将Ψ(z′)\Psi(z')Ψ(z′)表达式代入，可以得到：
① Ψi−1(k0Li−1)=Ψi(0)\Psi_{i-1}(k_{0}L_{i-1})=\Psi_{i}(0)Ψi−1(k0Li−1)=Ψi(0) -------> Wi−1(eλik0Li−1)Ci−1=WiCiW_{i-1}(e^{\lambda_{i} k_{0} L_{i-1}}) C_{i-1}=W_{i}C_{i}Wi−1(eλik0Li−1)Ci−1=WiCi
① Ψi(k0Li)=Ψi+1(0)\Psi_{i}(k_{0}L_{i})=\Psi_{i+1}(0)Ψi(k0Li)=Ψi+1(0) -------> Wi(eλik0Li)Ci=Wi+1Ci+1W_{i}(e^{\lambda_{i} k_{0} L_{i}}) C_{i}=W_{i+1}C_{i+1}Wi(eλik0Li)Ci=Wi+1Ci+1
由于Ci+1=TiCiC_{i+1}=T_{i} C_{i}Ci+1=TiCi，因此可以得到变换矩阵TiT_{i}Ti：Ti=(Wi+1)−1Wi(eλik0Li)T_{i}=(W_{i+1})^{-1} W_{i} (e^{\lambda_{i} k_{0} L_{i}})Ti=(Wi+1)−1Wi(eλik0Li)

（5）稳定性问题

我们从麦克斯韦方程组开始，推导出了用来表征电磁场随z′z'z′变化情况的场方程等式；而后使用一阶微分方程通解公式和矩阵性质，求出了该场方程的表达式；随后，利用各层之间的边界条件对各界面进行约束，求解出了传输矩阵 T。虽然看起来已经求得了我们需要的结果，但事实上，我们通过该方法得出的结果是不稳定的，或者说这个模型是不稳定的。
在一开始，我们写平面波的电场和磁场方程的时候，我们将其入射方向设为的是k⃗\vec{k}k，但是k波矢通常是一个复数。对于前向传播（forward propagation）的电磁波来说，其可以展开为k+=k′+jk′′k^{+}=k'+jk''k+=k′+jk′′；对于反向传播（backward propagation）的电磁波来说，其可以展开为k−=k′+jk′′k^{-}=k'+jk''k−=k′+jk′′，其中k′k'k′表现为波的振荡，k′′k''k′′表现为波的衰减或增强，这两种波都会存在于实际的器件层中，如下图。

但是在我们的模型中，是默认所有波为前向传播的，因为我们的方程是：E⃗(r⃗)=E0⃗⋅ejk⃗r⃗\vec{E}(\vec{r} )=\vec{E_{0}} \cdot e^{j\vec{k}\vec{r}}E(r)=E0⋅ejkr（前向），而没有考虑到E⃗(r⃗)=E0⃗⋅e−jk⃗r⃗\vec{E}(\vec{r} )=\vec{E_{0}} \cdot e^{-j\vec{k}\vec{r}}E(r)=E0⋅e−jkr（反向）。所以，为了让反向波能够像前向波一样被处理，我们将k−k^{-}k−的表达式变为：k−=−k′−jk′′k^{-}=-k'-jk''k−=−k′−jk′′，如下图所示。如果我们假设原来的k′′k''k′′是衰减的，那么−k′′-k''−k′′在指数上则是成增长的态势，是非常不稳定的。

Sort！
因此，我们在稳定性这一步要做的就是分离出前向波和反向波。此处，要引入坡印廷矢量，当S>0时，为前向波，反之为后向波。

举个栗子，假设 εr=9.0\varepsilon_{r} = 9.0εr=9.0 且 μr=9.0\mu_{r} = 9.0μr=9.0，则：

将WWW各项代入坡印廷矢量的公式，可以求出WWW的第一列和第三列结果为正（前向波），第二和第四列结果为负（后向波），此时我们将为正的放在一起，为负的放在一起，需要注意的是变换位置后的λ\lambdaλ矩阵要重新调整为单位矩阵的形式，如下图所示：

那么，我们的WWW和λ\lambdaλ矩阵可以简化为以下形式：

最终Ψ\PsiΨ的解可以写为如下形式：

然后再用和之前相同的方式求解出传输矩阵T，那么此时我们获得的解才是稳定的。

以上就是使用4×4矩阵来求解T矩阵，该方法的关键在于最后要考虑一个稳定性的问题。如前所述，除了该方法，我们还可以使用2×2矩阵来求解，该部分内容将写在下一篇中。

如有错误，欢迎指正…