截止零根与重根退化
从完整解到稳定传播

关键结论

退化层的传播从“必须有一组独立本征向量”改为直接保留对应微分方程的完整解。平面有限层精确零根保留线性解；圆锥入射的约化模式退化按层切换到完整 2N 阶 Schur 传播；Python 近场使用同一套传播与界面振幅。

重根本身不会自动造成错误；缺失独立解、模式基病态、错误的上下行分离会同时影响场和散射振幅。衍射效率是振幅与法向通量共同决定的量，因此能量守恒或某一级效率接近零，都不足以判定复电场正确。

本页讲的是 C++/Python 求解器在零根与 Schur 传播上的实现、回归测试和 COMSOL 对拍；电场图已嵌入页面，COMSOL 原始数据没有随发布。报告日期 2026-10-10。以下 N = 2M + 1，M 为 Fourier 截断阶数，纵向坐标 ξ = k₀z。

本征值正确性：数值根、独立子空间与传播

1. 三种退化需要分别判断

2. 小本征残差不能证明模式基可用

残差检验每一对 (λⱼ,vⱼ) 是否满足本征方程；κ₂(V) 检验这些向量合起来是否能稳定表示任意场。重复返回同一个合法向量时，每列残差仍可为零，但 V 已经奇异。RCWA 后续的模式展开与界面匹配需要完整、可用的基，因此必须区分这两个问题。

例如 A = diag(0,0,1)，把前两列都选成 e₁，三对本征残差都是零，却丢失了 e₂。另一个例子 A = [[0,1],[0,ε]]：当 ε → 0，两列本征向量几乎平行，κ(V) 随 1/|ε| 增大；ε = 0 时只剩一个本征向量。增加小扰动可以分开数值根，却改变了所求问题。

3. 非厄米问题的判据

一般非厄米、非正规算符的本征向量本来就可以不正交。orthogonality_defect(V) 适用于考察正规问题的正交性，不能作为一般非厄米问题的失败条件。对于厄米重根，可在同一本征子空间内构造正交基；对不同本征值直接做 Gram–Schmidt 会破坏 Av = λv。对一般非正规或亏损问题，Schur 分解保留整个不变子空间及其耦合，是 StableRCWA 在退化层采用的方法。

Schur 向量是稳定的子空间坐标；T 的非对角项记录模式间耦合。只取 diag(T) 并把 Z 当作普通本征向量，会再次丢掉这些耦合。这里沿用通用 QR 的变换累积接口，RCWA 直接消费 Z、T；通用 eig() 的本征向量回代算法保持原实现。

4. 当前 QR 回代在零重根附近的实际行为

通用 QR 先经 Householder 化为 Hessenberg 形，再用酉 QR 迭代得到上三角 T，最后回代 (T−λⱼI)xⱼ=0 并乘回 Z。回代公式 xᵢ = Σ_k>iTᵢₖxₖ/(λⱼ−λᵢ)，问题集中在谱间隔很小的分母。

现有实现对严格为零的分母使用 max(eps,机器精度)×max|Tᵢⱼ| 替代；若该量下溢则使用最小正规数。对非零但很小的分母，则使用 max(eps,机器精度)×max(|λᵢ|₁,|λⱼ|₁) 作下限，|·|₁ = |Re|+|Im|。

零根簇中若 |λ| 本身已是 O(u‖A‖)，后一阈值可以只有 O(u²‖A‖)，而 Schur 非对角扰动仍可为 O(u‖A‖)。因此商依然可能很大，归一化后的向量依然可能几乎相关。分母保护能控制部分数值溢出，却不能保证重根子空间的秩；对于真正亏损的矩阵，独立本征向量在数学上就不存在。退化层直接保留 Schur 表示，避开了该回代需求。

5. 病态检测与可疑结果复核分别检查什么

本征对残差回答“每列是否是近似解”；σ_min 与 κ₂ 检查“这些列能否展开所有场”；Schur 残差 ‖HZ−ZT‖/‖H‖ 与酉性缺损检查“子空间表示是否正确”；界面线性系统残差、切向场连续性和 Maxwell 方程残差检查“解是否正确接入物理问题”。重复根之间应比较子空间或重建的场，逐列相位或排序的差异不构成物理误差。

当前生产路径采用约化谱的逐层触发判据、平方根残差和有限数检查；独立一阶场、界面连续性与 Maxwell 差分残差在回归测试中复核。全模式条件数和上述全部指标没有作为每次 RCWA 调用的固定开销加入。这样明确区分了自动运行的判据和测试时的交叉验证。

6. 为什么这种误差会传到衍射效率

若完整界面系统写为 Bc = b，数值扰动满足一阶关系 δc ≈ B⁻¹(δb − δB·c)。模式基相关会放大 B 的条件数，进而放大内部振幅和出射振幅误差。场 E = Fc 同时受到 δF 和 δc 影响；某个通道 η = f|J|²，在固定通量权重 f 下有：

因此不存在一个只由“有重根”就能确定的误差百分比。误差大小取决于丢失的解、激发方向、层厚、界面条件数和观测位置。临界外部阶 f → 0 时，其自身法向功率趋小，复场仍可有限，并通过界面耦合改变其他传播级。后文同时给出复振幅、效率和复电场误差。

平面有限层精确零根：完整推导与实现

1. q = 0 不是把指数设成 1 就结束

每个可解耦模式满足 f″(ξ) = q²f(ξ)。非零 q 的解可由 exp(±qξ) 表示；零根处必须保留它们的合流极限：

传播状态矩阵含非零的 h 项；将 q 换成一个小数或者让正反两列都变成常量，都会改变这一传播关系。

下面的图解把 q 一步步减到零，可以看到标准写法的两列怎样并成同一个函数。

2. 辅助界面基与有限散射块

物理 q 保持为零，仅在有限层的界面基中令该列 basis_q = 1。TE 的配对基为 V = W·diag(basis_q)，TM 为 V = ⟦1/ε⟧W·diag(basis_q)。这些辅助振幅表示场及导数，并非两束可独立辨认的行波。

设上表面入射辅助振幅为 a，下表面入射辅助振幅为 b；上表面出射为 u，下表面出射为 v。零模式的场与归一化导数在上表面为 a+u、−a+u，在下表面为 v+b、−v+b。由导数恒定、场跨层增量为 h 倍导数：

这个非零“反射块”是辅助坐标下的传播关系。与两侧真实材料界面组合后，得到物理反射和透射。d = 0 时 T = 1、R = 0，恰好退化为恒等传播。

下面的图解可以拖动层厚 h，看层内的直线场与 T、R 怎样随之变化。

3. 层内场和触发标准

Python 用距上下两表面的距离分别计算系数，避免 h−ξ 的消减误差。C++ 的 t₀ 写成 1/(1+0.5*h)，避免直接计算 2+h 时溢出；非有限厚度参数进入数值失败分支。

此分支对每个具有完整空间模式的零根，精确实现对应二阶方程的基本解。它解决的是有限层上下行合并；空间本征问题本身的 Jordan 链属于另一种结构。

4. 可核算的解析对照

n_I = n_II = 2，θ = 30°，λ = 1 μm，n_film = 1，d = 1/(π√3) μm。严格浮点分支取 n_film = 2 sin(30°) 的同一机器数。结果 TE：R = T = 1/2；TM：R = 1/17、T = 16/17。

剥除横向相位后，TE 电场在 z/d = 0、1/4、1/2、3/4、1 处依次为 1.5−0.5i、1.25−0.25i、1、0.75+0.25i、0.5+0.5i。非零线性斜率直接验证了第二个解被保留。

圆锥入射：按层切换完整 2N 阶 Schur 传播

1. 约化本征值 λ = 0 与真实 q = 0 的区别

原算法做两次 N 阶本征分解，再构造 2N 个完整电磁模式。构造式中第一组有 λ₁q₁，第二组有 W₂λ₂ 等因子。kᵧ ≠ 0 且 λ 接近零时，这些分量消失，两组列可能趋于相关；此时 q 可以仍等于 |kᵧ|。均匀介质在这种位置可以拥有完整半单本征空间，因此“构造式退化”与“完整 Maxwell 算符亏损”必须区分。

下面的图解用 COMSOL 对拍里的圆锥结构（M = 8）显示：越靠近约化零根，拼出的两列越接近平行，而每列的本征残差一直停在舍入误差水平。

2. 完整算符及其与一阶 Maxwell 系统的关系

在直角切向坐标下，令 E = ⟦ε⟧，原始 P、Q 为下式；实现随后用逐级次的电场、配对磁场旋转矩阵变换到 (s,t) 坐标。

旋转后的切向电场 F、配对磁场 U 满足 F′ = PU、U′ = QF，因此 F″ = PQF。设 HZ = ZT，G² = T，令：

于是下行分量 F = W exp(−ξG)a、U = −V exp(−ξG)a，和上行分量 F = W exp[−(h−ξ)G]b、U = V exp[−(h−ξ)G]b，均逐项满足一阶 Maxwell 方程。这给出了 Schur 传播与原物理系统的直接等价关系。代码通过右上三角求解 VG = QZ 得到 V，避免显式求逆。

3. 上三角平方根与 Jordan 项

平方根按原求解器分支选择，使传播模保持相位约定、衰减模向外衰减。C++ 与 Python 均验证 ‖G²−T‖₁ ≤ 10⁻¹⁰ max(‖T‖₁,tiny)。精确零对角根、零除数或非有限中间量会触发数值失败。

若 T = q₀²I + J，J² = 0、q₀ ≠ 0，则：

hJ 项就是广义模式的传播。删除上三角项，会得到错误的场和错误的层散射块。测试用 T = [[4,3],[0,4]] 核对 G₁₂ = 0.75 及矩阵指数的非对角项。

下面的图解在这个 T 的右下角加 ε，比较逐模式写法与 Schur 写法算出的矩阵指数。

4. 与散射递推的接合

非零 G 的层传播块为 T_dd = T_uu = X，R_du = R_ud = 0。界面继续由 W、V 匹配切向场，随后使用既有 S 矩阵递推。每个方向的振幅锚定在它进入该层的表面，传播只经过非负距离，厚吸收层也沿这一约定处理。

更换基 W → WD、V → VD 时必须同步令 G → D⁻¹GD、a → D⁻¹a。原算法的逐列独立归一化只适用于对角传播；Schur 分支保留 Z 的酉基，不再独立缩放 W、V 的列。

5. 数值实现与可观察结果

C++ 的矩阵指数使用缩放 Taylor 展开及平方恢复，每次平方后用标量指数刷新对角线；这避免 diag(−10²⁰,−1) 中较慢衰减的分量在大量平方中被误算。空矩阵、10^±200 平方根尺度、Jordan 耦合和溢出失败均有专门测试。

n_schur 记录实际切换层数，每层额外一次 QR 计入 n_eig。min_abs_lambda 和 degenerate_mode_warning 保留近退化诊断含义。自定义本征后端继续负责原两次 N 阶分解，Schur 回退使用内置 QR。这里“回退”由层内判据提前选择，不依赖整次求解产生 NaN 后重新猜测原因。

Python 近场：与层散射保持相同的解

1. 一个共享传播入口

新增 gui/mode_propagation.py。平面 rcwa1d.py 与圆锥入射的 modes_conical.py、rcwa_conical.py、fields_conical.py 共同使用它。ModeSet.schur_generator 保存 G，普通层保存 q。求解层散射块、恢复层内振幅、单点场和场图都读取同一模式表示。

零根分支使用上一节的线性系数，普通层使用逐模式指数。之后按原材料卷积规则恢复纵向分量并做 Fourier 合成，输出 Ex、Ey、Ez、Hx、Hy、Hz。层内模式振幅属于所选基，跨实现比较应在物理复场或外部 Jones 振幅上进行。

2. 计算场图时保留耦合，同时控制内存

Python 的层 S 块用 scipy.linalg.expm；多深度场图使用 expm_multiply 直接计算指数对振幅的作用，场采样输出的存储量为 O(N·n_z)。对恰好等间距的深度采用区间接口；任意、重复、含界面的深度通过 unique/inverse 映射恢复原顺序。这样每一个请求深度都对应正确位置。

3. 同时修复未激发倏逝入射分量的溢出

原实现先为所有入射阶计算指数，在离结构较远的位置，未激发的倏逝阶也可能得到 Inf，随后 0×Inf 污染整个场。现在仅对非零入射振幅计算入射指数。平面和圆锥的单点场、场图均已修改；z = −100 μm 的回归测试在浮点溢出与 invalid 均设为抛错时通过。

4. GUI 与调用约定

GUI 返回 python_n_schur 与 python_zero_modes，显示实际使用的近场方法。mode_method="schur" 可强制非均匀层采用 Schur 进行诊断；默认 reduced 模式按层自动检测。C++ 的 JSON 和示例同步输出 Schur 层统计。

验证：从矩阵函数到 Maxwell 方程

验证覆盖解析解、独立一阶系统、界面连续性、微分方程残差、C++/Python 一致性和外部 FEM。各项针对不同错误来源。

独立一阶参照的构造

测试内部采用辅助状态 [F;−U]，所以实际指数生成矩阵写作 [[0,−P],[−Q,0]]，返回物理场时恢复 U 的符号。对这个等价的 4N 阶系统计算一个小厚度的矩阵指数，将传输关系换成 S 矩阵，再通过倍增组合到目标厚度。场点上下各自组装散射系统，从该点匹配得到 F、U。这一路径与生产实现的 2N 阶 Schur 平方根有不同的计算流程，能检出符号、坐标旋转、G 耦合、层厚锚定和界面接合错误。

另外的混合层数值复核中，与一阶参照的最大 Jones 差为约 1.21×10⁻¹¹；M = 16 的复核约为 5.88×10⁻¹²。自动测试的源文件见页面末尾；COMSOL 对拍脚本没有随发布。

实现中的失败与诊断边界

检测到 QR 不收敛、奇异线性系统、无效矩阵函数或非有限数时，C++ 按原接口返回 NaN 结果，Python 抛出相应数值异常。圆锥 Schur 平方根的对角根若精确为零，仍触发此失败路径；它使用 G⁻¹ 的场配对关系，因此与平面零根辅助基是两种不同分支。当前外部精确掠射的部分界面路径也仍含零分母；对应条件在头文件和零厚度/外部截止回归中明确保留。

COMSOL 对拍：关键数字与电场

Windows COMSOL 6.4 实算。原批次 27 个模型；PML 主批次 90 个，另有 2 个场图与 2 个探针模型。每模型分别计算 S/P 激励。复场按已核实的坐标、相量和入射单位电场归一化对齐，保留绝对相位。

圆锥结构：周期 0.574 μm，光栅厚 0.38 μm，齿宽 0.1722 μm，齿/槽折射率 1.46/1，空气入射、衬底 1.46，θ = 40°、φ = 25°。内部约化零根定位于 M8 的 λ* = 0.980484673140678 μm，此时真实 q ≈ 0.271654，衬底 m=+1 的 k_z/k₀ ≈ 0.889288。真正衬底 m=+1 截止在 λ_II = 1.1577971963084606 μm；以下 δ = λ/λ_II−1。

底部 Port 与 PML 对复电场的影响

δ = +10⁻⁷、z = 0.88 μm、RCWA M200。原底部 Port 配置遗漏了已进入倏逝区的临界级；改为底部真实 PML 后，线采样误差如下。PML 厚 2 μm、factor = 2000、128 层，物理衬底缓冲 1.6 μm。

此处 E = (Ex,Ey,Ez) 为完整复电场；峰值误差为 max_x√Σ_α|E_COMSOL,α−E_RCWA,α|²，相对 L₂ 为 ‖ΔE‖_F/‖E_RCWA‖_F。效率差采用 0–1 的功率比例。

二维场图：结构内部的局部误差

δ = −10⁻⁵，PML factor = 200、96 层，96×96 物理域采样。下图共享强度色标，右列为完整复电场误差。M120 → M200 时，全图相对 L₂ 误差随截断提高而下降；介质不连续处的点态峰值仍受 Fourier 截断影响。

已确认的边界敏感性

临界阶的直接复 E/Z₀H 分量误差随 PML 尺度变化明显。下表相对 Python M120，取两采样面、两种入射的最大临界阶分量误差；它与上面的空间电场向量范数是不同指标。

名义 δ = 0 的实际浮点 k_z/k₀ 约为 1.5–2.1×10⁻⁸，最大尺度下直接场误差仍约 1.16×10⁻⁴。数学上的严格 k_z = 0 常量模在有限 PML 内不产生指数衰减。另一方面，近截止处前后向谐波分离包含 1/k_z，会放大场插值误差；直接场与尺度扫描揭示了这一影响。能量 R+T 接近 1 与较大的局部场误差在这些算例中同时出现。

代码、证据与复现

页面内图片全部嵌入，可离线阅读。COMSOL 的建模脚本、求解日志与复场原始数据没有随发布；上表的代码链接指向 v0.1.0 版本的源码。

cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build
cd build && ctest --output-on-failure && cd ..

# 使用已安装 NumPy / SciPy 的 Python；从项目根目录执行
cd gui
RCWA_JSON_EXE=../build/rcwa_json python3 -m unittest \
  validation.test_degenerate_fields \
  validation.test_cpp_backend.Comparison \
  validation.test_cpp_backend.Routing \
  validation.test_cpp_backend.WorkerUsesCpp \
  validation.test_cpp_backend.DegenerateModes
cd ..