多层介质：散射矩阵的星积与逐层场恢复

目的

多层结构可以逐层向外包入：先把末端半空间当作已知响应，再依次加入界面和传播段，直至得到整体散射矩阵。计算内部场时，沿保存下来的子结构响应逐层向内恢复，就像重新剥开已经包好的洋葱。

“包起来”消去内部相向波，留下两端口关系；“剥开来”在指定层重新求出这两个波。两步都使用同一组散射方程，不需要将 S 转换回 T。本文先推导层间连接与场恢复，再介绍多层薄板程序的核心函数及 benchmark。

仅有最终整体 S 的四个元素不能确定所有内部场；还需逐层几何和中间子结构响应。程序因此保存左右累计矩阵。整体 T 也需要相应局部关系才能恢复内部。

1. 每层场与端口的统一定义

每层正反向波

图 1：层号由左到右增加，最外为半空间。A 始终朝正 z，B 始终朝负 z；只有在最左、最右外端口，才能将它们分别直接称为外部入射或出射。内部一个 A 波，是左子结构的出射，也是右子结构的入射。

固定横向波数与一种偏振 s/p，逐通道求解。对于第 j 层，；根取出射或衰减分支。若 A、B 都在参考点定义，

$\mathbf E_j(x,z)=e^{i\kappa x}\left[A_j\mathbf e_{\sigma j}^+e^{iq_j(z-z_c)} +B_j\mathbf e_{\sigma j}^-e^{-iq_j(z-z_c)}\right].\tag{3.1}$

使用第一、二篇的偏振基矢。半空间的参考平面取最近界面；先以有限层右端为共同参考推导，再说明程序的稳定场重构接口。

一个子结构 X 的 S 定义为

$\begin{pmatrix}b_L\\b_R\end{pmatrix} =S_X\begin{pmatrix}a_L\\a_R\end{pmatrix},\qquad S_X=\begin{pmatrix}X_{11}&X_{12}\\X_{21}&X_{22}\end{pmatrix}.\tag{3.2}$

其中 a 表示进入该子结构的波，b 表示离开该子结构的波。左入射，右入射，这里大写 L 同时是最后一层的层号；外部出射。使用小写 a/b 是为了不把“端口方向”与“材料层号”混在一起。

2. T 的级联与 S 的级联，区别在哪里？

T 的定义是同一参考面的正反向波配对。若 X 在左、Y 在右，且连接处使用同一基矢、同一参考面，

$\begin{pmatrix}A_{left}\\B_{left}\end{pmatrix} =T_X\begin{pmatrix}A_{mid}\\B_{mid}\end{pmatrix},\qquad \begin{pmatrix}A_{mid}\\B_{mid}\end{pmatrix} =T_Y\begin{pmatrix}A_{right}\\B_{right}\end{pmatrix}.$

直接代入就得到。S 的中间两波则分别属于两个子结构的不同输入/输出，不能用普通矩阵乘法消去。需要重新连接端口。

3. 用连接处的 u、v 推导星积

两个散射块的连接

图 2：X 在左，Y 在右。u 向右，是 X 的右侧出射和 Y 的左侧入射；v 向左，是 Y 的左侧出射和 X 的右侧入射。中间虚线是参考切面，不是新材料界面。图中星号表示下面定义的级联运算。

两个子结构分别满足

$\begin{pmatrix}b_L\\u\end{pmatrix} =S_X\begin{pmatrix}a_L\\v\end{pmatrix},\qquad \begin{pmatrix}v\\b_R\end{pmatrix} =S_Y\begin{pmatrix}u\\a_R\end{pmatrix}.\tag{3.3}$

把四行都展开：

$\begin{aligned} b_L&=X_{11}a_L+X_{12}v,\\ u&=X_{21}a_L+X_{22}v,\\ v&=Y_{11}u+Y_{12}a_R,\\ b_R&=Y_{21}u+Y_{22}a_R. \end{aligned}\tag{3.4}$

中间两式是一对反馈方程。将 v 代入 u，令 $\delta=1-X{22}Y{11}$，得到

$u=\frac{X_{21}a_L+X_{22}Y_{12}a_R}\delta,\qquad v=\frac{Y_{11}X_{21}a_L+Y_{12}a_R}\delta.\tag{3.5}$

再代入 bL、bR，整理成，于是

$\boxed{S_X\star S_Y= \begin{pmatrix} X_{11}+\dfrac{X_{12}Y_{11}X_{21}}\delta&\dfrac{X_{12}Y_{12}}\delta\\ \dfrac{Y_{21}X_{21}}\delta&Y_{22}+\dfrac{Y_{21}X_{22}Y_{12}}\delta \end{pmatrix},\quad\delta=1-X_{22}Y_{11}.}\tag{3.6}$

这就是本文定义顺序下的 Redheffer 星积：左边的 X 先经过，右边的 Y 后经过。有些文献把它写成带圈乘号等其他符号，也有相反的记号顺序。符号本身不能决定四个元素怎样排列，必须先明确式 (3.3) 的端口定义。本文使用，避免与通常表示 Kronecker 积的符号混用。

3.1 为什么多出一个分母？

$X{22} $是从连接处入射到左子结构的反射，$ Y{11}$ 是从连接处入射到右子结构的反射。它们的乘积表示在两个子结构之间往返一次。因此分母与薄板求和中的反馈分母是同一个来源。星积把这些往返全部计入，不是只保留一次穿过两个结构的路径。

3.2 星积的几条运算规则

在端口兼容、相关反馈逆存在的条件下，星积满足结合律，可以从左向右累积，也可以从右向左累积；一般不满足交换律。它不是逐元素相乘，也不是矩阵相加。

空结构的波直接穿过，因此星积单位元是

$J=\begin{pmatrix}0&1\\1&0\end{pmatrix},\qquad J\star S=S\star J=S.\tag{3.7}$

最简单的区别是，而普通乘法。因此即使两个块都“没有反射”，也不能把 S 直接普通相乘。

对于当前各向同性问题，每种偏振只有一个通道，式中元素是复数。若元素本身变成多通道矩阵，不能继续交换乘法次序或把分母当成一个数，见文末说明。

4. 可以一层一层加入吗？具体按什么顺序？

可以。将每个界面和每个有限层传播段各看成一个 S 块：

$S^{(j,j+1)}=\begin{pmatrix}r_{j,j+1}&t_{j+1,j}\\t_{j,j+1}&r_{j+1,j}\end{pmatrix},\qquad C_j=\begin{pmatrix}0&f_j\\f_j&0\end{pmatrix},\quad f_j=e^{iq_jd_j}.\tag{3.8}$

这里 Cj 是第 j 层的传播 S，不是传输矩阵。最外半空间不插入无限厚度的传播段。按几何顺序排出

$G_1=S^{(1,2)},\ G_2=C_2,\ G_3=S^{(2,3)},\ G_4=C_3,\ldots, G_M=S^{(L-1,L)},\quad M=2L-3.\tag{3.9}$

从开始递推

$\boxed{F_m=F_{m-1}\star G_m,\quad m=1,\ldots,M,\qquad S_{total}=F_M.}\tag{3.10}$

这就是“一层层加入”的准确含义：加入界面、加入传播段，并执行星积，不是把各层 S 相加。若左侧已经是一整个薄膜堆栈，它的整体 S 也能继续与新块作星积。

5. 包好以后保留什么：各层右端的子结构响应

先采用 CoeAB_layer_SMM 使用的共同右参考面。令第 j 个有限层范围为，厚度为。原代码规定

$z_{ref,j}=z_j^+,\qquad \mathrm{ABCoe}\{j\}=\begin{pmatrix}A_j\\B_j\end{pmatrix}. \tag{3.11}$

A、B 都在同一右端参考面，不是“一个在左端、一个在右端”。第 1 层的参考面是入口界面 z=0；最后一层没有有限右端，参考面取出口界面。

原程序在每层右端切开并恢复两个波

图 3：虚线切在第 j 层右端、尚未跨过下一材料界面的位置。左边 Pj 包括此前所有结构及第 j 层的传播段；右边 Qj 从该切面开始，包括下一个材料界面和其后结构。两条内部箭头是同一平面的 Aj、Bj。

记界面 S 为，有限层传播 S 为，那么

$P_1=J,\quad P_j=I_1\star C_2\star I_2\star\cdots\star I_{j-1}\star C_j \quad(2\leq j\leq L-1),\quad P_L=S_{total}. \tag{3.12}$

右侧子结构是

$Q_j=I_j\star C_{j+1}\star I_{j+1}\star\cdots\star C_{L-1}\star I_{L-1}, \quad Q_{L-1}=I_{L-1},\quad Q_L=J. \tag{3.13}$

P、Q 在这里是 2×2 散射矩阵。尤其 不再包含第 j 层传播，因为切面已经移到了它的右端。这与“在层左端切开”的示意不同，不能把两个切面的反射系数混用。

6. 从外层包入：左右子结构的逐层递推

左向右递推的初值为 P1=J。先跨过界面 j，再传播完整个下一有限层：

$P_{j+1}=P_j\star I_j\star C_{j+1}\quad(j+1<L),\qquad P_L=P_{L-1}\star I_{L-1}. \tag{3.14}$

这正是原函数先 cascade_S(curr_SL,S_int{i})，再在下一层有限时接 S_prop{i} 的原因。数组 S_left{j} 保存的就是 Pj。

右向左递推从 QL=J 开始。为了从第 j+1 层右端回到第 j 层右端，需要将“第 j+1 层传播”和“界面 j”依照几何顺序接到现有结构前面：

$Q_j=I_j\star C_{j+1}\star Q_{j+1}\quad(j+1<L),\qquad Q_{L-1}=I_{L-1}. \tag{3.15}$

所以原代码先执行 cascade_S(S_prop{i},curr_SR)，再执行 cascade_S(S_int{i},curr_SR)，不是把两个因子顺序写反。最终 S_right{j} 保存 Qj。

例如四层结构，完整链为：

切面	左侧 `S_left`	右侧 `S_right`
第 1 层右端，即入口面	J	全部结构
第 2 层右端
第 3 层右端
第 4 层左端，即出口面	全部结构	J

因而任意切面都满足

$\boxed{S_{total}=P_j\star Q_j.}\tag{3.16}$

这张表给出了左右累计矩阵的保存范围，可直接对应程序中的两次循环。

原代码变量	本文解析量	下标的具体含义
`S_int{i}`		介质 i 与 i+1 的界面
`S_prop{i}`		第 i+1 个材料层的传播，厚度 `d_layer(i)`
`S_left{i}`		从入口到第 i 层右参考面
`S_right{i}`		从第 i 层右参考面到出口
`ABCoe{i}`		同一参考面的两个电场波系数

所以传播数组的第 i 项与材料层的第 i+1 层对应，不能把这两个 i 当成相同层号。

7. 逐层剥开：从左右响应恢复 Aj、Bj

整个结构的外部入射是、。对图 3 的左、右子结构分别写端口关系：

$\begin{pmatrix}b_L\\A_j\end{pmatrix} =P_j\begin{pmatrix}a_L\\B_j\end{pmatrix},\qquad \begin{pmatrix}B_j\\b_R\end{pmatrix} =Q_j\begin{pmatrix}A_j\\a_R\end{pmatrix}. \tag{3.17}$

从第一张矩阵的第二行、第二张矩阵的第一行，得到

$A_j=(P_j)_{21}a_L+(P_j)_{22}B_j,\qquad B_j=(Q_j)_{11}A_j+(Q_j)_{12}a_R.$

P21 是外部左入射传到切面的透射；P22 是从切面朝左看见的反射；Q11 是从切面朝右看见的反射；Q12 是外部右入射传到切面的透射。这四个量都是整段子结构的响应，不是某一个邻近界面的 Fresnel 系数。

联立后得到

$\boxed{\begin{aligned} A_j&=\frac{(P_j)_{21}a_L+(P_j)_{22}(Q_j)_{12}a_R} {1-(P_j)_{22}(Q_j)_{11}},\\ B_j&=\frac{(Q_j)_{11}(P_j)_{21}a_L+(Q_j)_{12}a_R} {1-(P_j)_{22}(Q_j)_{11}}. \end{aligned}}\tag{3.18}$

这与原代码末尾 Ai、Bi 的两行表达式逐项相同。分母表示切面两侧子结构之间的往返反馈，外部双侧入射也已包括。

左入射 1、右入射 0 时，才简化为

$A_j=\frac{(P_j)_{21}}{1-(P_j)_{22}(Q_j)_{11}},\qquad B_j=(Q_j)_{11}A_j.\tag{3.19}$

因此“通过 S 恢复场”的实质，是保留被总 S 消去的中间切面，再解其内部两波。原 CoeAB_layer_SMM 已经实现，不需要再造一个求解器。

8. 对应 Get_Field_From_ABCoe 的场重构

令 $\zeta=z-z{ref,j} $，其中有限层$ z{ref,j}=z_j^+$。每层的矢量场是

$\boxed{\mathbf E_j(x,z)=e^{i\kappa x} \left[\mathbf e_{\sigma j}^+A_je^{iq_j\zeta} +\mathbf e_{\sigma j}^-B_je^{-iq_j\zeta}\right].}\tag{3.20}$

原 Get_Field_From_ABCoe 按同样的右端参考平面，计算标量。对于 TE，这正是 Ey。在 TM 下它是什么物理分量，要结合第 10 节的偏振约定判断，不能都叫 Cartesian Ez。

第二篇曾将板内 A2、B2 都定义在左端。两套系数并不矛盾，换参考面即可：

$\boxed{A_j^{right}=A_j^{left}e^{iq_jd_j},\qquad B_j^{right}=B_j^{left}e^{-iq_jd_j}.}\tag{3.21}$

在同一个位置代入式 (3.20)，结果应相同。原代码返回的是 right 这组，不要拿它直接与左参考系数比较。

9. 用第二篇薄板核对原程序

空气 | n=1.5、100 nm | 空气，600 nm、正入射 TE，，。原程序在第 2 层右端切开，所以

$P_2=\begin{pmatrix}-0.2&1.2i\\0.8i&-0.2\end{pmatrix},\qquad Q_2=\begin{pmatrix}0.2&0.8\\1.2&-0.2\end{pmatrix}.\tag{3.22}$

P2 已包含一次板内传播，Q2 只是右界面。恢复公式给出

$\boxed{A_2=\frac{0.8i}{1-(-0.2)(0.2)}=\frac{10i}{13},\qquad B_2=0.2A_2=\frac{2i}{13}.}\tag{3.23}$

换回左参考面，得到第二篇的。因此原程序的标量 TE 场为

$E_{2y}(z)=\frac{10i}{13}e^{iq_2(z-d)}+\frac{2i}{13}e^{-iq_2(z-d)} =\frac{10}{13}e^{iq_2z}-\frac2{13}e^{-iq_2z}.\tag{3.24}$

可以直接调用已有原函数核对，不需要新增函数：

% 从发布目录根部，使用完整保留的原仓库函数。
addpath('modules/01_PlaneWaveMultilayer/examples/original_transfermatrix/Functions');
ABCoe=CoeAB_layer_SMM(100e-9,[1 1.5 1],600e-9,0,"TE",1,0);
disp(ABCoe{2});                 % [10i/13; 2i/13]，都在右界面。
kz=2*pi/600e-9*[1 1.5 1];
z=linspace(-100e-9,200e-9,601);
Ey=Get_Field_From_ABCoe(z,ABCoe,100e-9,kz);

这里直接运行共同右参考面的接口；ABCoe、S_left、S_right 的层号一致，可逐项检查第 2 层的解析结果。

10. 核心函数与三个程序的分工

整个项目只组织三部分核心程序。第一部分计算无源多层薄板的反射、透射与内部场；第二部分把偶极子展开成平面波，调用第一部分，再组合 Green 张量、Purcell 因子与远场；第三部分处理球外偶极子的球多极子辐射。

本篇介绍第一部分。原始接口保留在 TransferMatrix 的快照中，整理后的公共入口位于模块 1：

计算步骤	原函数	整理后的入口
单界面 s/p 散射	interface_S	mlpw.interface
散射块星积	cascade_S	mlpw.cascade
整体响应、左右累计及层内系数	CoeAB_layer_SMM/TMM	mlpw.solve
空间场分布	Get_Field_From_ABCoe	mlpw.field

原接口的 ABCoe 使用共同右参考面，TE 的标量场是 Ey。通用 Cartesian 场应按偏振基矢恢复，不能把所有偏振都写成一个标量 A+B。整理后的 mlpw.field 返回 E、H 三个分量，内部使用 Aleft、Bright 两端参考来避免增长指数；两端参考和共同参考之间仅相差已知传播因子，边界物理完全相同。

模块 2 的 mldipole.rt 只负责上下结构的层序、厚度和源平面的相位转换。它直接调用模块 1 的 mlpw.solve 获取每个角谱分量的反射、透射，不另写 Fresnel 系数或多层递推。对角各向异性使用的平面波本征模与场重构也归模块 1；模块 2 只完成源投影及互易远场组合。

11. 数值实现与 benchmark

从发布目录运行：

1
2
3

setup;
demo_original_transfermatrix; % 四个原结构的复算与比较
run_tests;                   % 三个核心程序的验证

首先用第二篇的单块玻璃板核对复振幅与板内系数；再检查多层 TMM/SMM、切向 E/H 连续与功率守恒。复场对比须统一 TM 基矢和参考面，不能只比较未转换的系数符号。

基准	检查量	已执行结果
100 nm 玻璃板、600 nm、TE	r、t 与解析分数	-5/13、12i/13
普通双层薄膜	SMM/TMM 复系数差	2.31e-16
空间界面检查	切向 E/H 采样残差	6.61e-9
五周期四分之一波长镜	新旧实现反射率差	9.13e-15
原布拉格场分布例子	内部 TE 复场差	约 3.7e-15
20 片高折射率层	新旧实现功率差	不超过 8.78e-12

这些指标分别检查归一化、矩阵级联和内部场恢复。界面采样残差包含跨越界面的有限采样距离，不是纯线性方程残差；不能将不同误差定义当成精度排名。极厚吸收/倏逝层和窄共振需要另外检查数值范围。

四分之一波长镜的反射透射谱

图 4：空气 | (HL)^5 | 空气，nH=3.5、nL=1.45，设计波长 1300 nm，各层光学厚度为四分之一波长，正入射 TE。

布拉格结构内部场

图 5：空气 | (HG)^4H | 空气，nH=2、nG=1，dH=2L、dG=L，omega L/c=3。图 4、5 是两个不同配置。另两个原例子——20 片层频谱与吸收光刻胶驻波——及全部参数、脚本、CSV 和比较记录见完整 benchmark 集合。

下一篇在这里已经建立的多层散射程序上引入点偶极子的角谱与 Green 张量，第五篇再单独展示 Purcell 因子和远场图的数值结果。

补充：元素变成矩阵时的星积

以上 s/p 独立通道的四个元素是复数。若以后处理耦合偏振或多个模式，X、Y 的每个元素块会变成矩阵。定义 $Du=I-X{22}Y{11} $、$ D_v=I-Y{11}X_{22}$，相同消元给出

$\begin{aligned} S_{11}&=X_{11}+X_{12}D_v^{-1}Y_{11}X_{21},&S_{12}&=X_{12}D_v^{-1}Y_{12},\\ S_{21}&=Y_{21}D_u^{-1}X_{21},&S_{22}&=Y_{22}+Y_{21}D_u^{-1}X_{22}Y_{12}. \end{aligned}\tag{3.28}$

这时次序不能交换，也不能把 Du 与 Dv 当成同一个分母。若以后扩展程序，应解线性方程，例如 Du\X21，而不是显式求逆。cascade_S 实现的是本篇的标量通道，不能直接把每个元素换成矩阵后照抄标量除法。

关于散射连接及其数值背景，可参阅 R. C. Rumpf, “Improved Formulation of Scattering Matrices for Semi-Analytical Methods That Is Consistent with Convention,” Progress In Electromagnetics Research B 35, 241–261 (2011), DOI:10.2528/PIERB11083107。本文的实际公式均由式 (3.3) 的端口关系直接消元，不依赖猜测文献里的星号方向。