2.4 一次平面波展开仿真的完整复盘

2.4 一次平面波展开仿真的完整复盘

本节摘要：本节把一次真实的平面波展开仿真从头到尾复盘：问题设定、参数选择、收敛检查、结果解读、两个经典错误的现场排除。全部代码可运行，输出有对照基准（上一节的手推结论）。读完你应当具备独立跑通一维与二维带隙计算、并判断结果可信度的能力。

这次仿真要回答什么问题

上一节在纸上推出了带隙公式，但工程上真正面对的问题往往超出解析能力：非对称占空比、多层层叠、二维晶格。仿真就是为了接管这些区域。本节的复盘对象是一个具体而完整的任务：**给定硅与空气构成的二维三角晶格圆孔阵列，孔径比 0.45，求 TM 偏振的第一带隙，并回答三个问题——带隙在哪、多宽、可信度如何。**复盘按真实工作顺序走：先定参数，再查收敛，然后出结果，最后做交叉验证。每一步都交代"为什么这么做"，因为这个流程的价值不在代码本身，而在它防的那几个坑。

第一步：参数选择的账本

复盘从一笔看似枯燥却最常出错的账开始。晶格常数不必预先决定——由于缩放律，计算全部在归一化单位下进行，晶格常数取一即可；落地时才把归一化带隙中心换算到目标波长，反推晶格常数。介电常数取平方：本征方程吃的是介电函数 ε，不是折射率 n；把 n 当 ε 输入是初学者最常见的静默错误，输出看似合理却整体错位。偏振要显式声明：本任务按 TM（电场沿孔轴）计算，因为 TM 偏振在空气孔结构里最容易出宽带隙；若要 TE 结果需另跑一遍，两者不可混用。截断数从粗到细：先用每方向 7 个平面波粗算轮廓，收敛检查通过后再上 11 或 13——直接上大截断会在坏参数上浪费数小时机时，这是过来人的经验。

第二步：完整实现代码

下面是可独立运行的完整实现。介电函数的傅里叶系数没有解析式的情形（圆孔截面）用数值积分求出，这是通用做法；核心是把本征方程组装成矩阵后交给特征值求解器。

import numpy as np
from scipy.special import j1

# ---------- 结构参数 ----------
eps_a, eps_b = 1.0, 3.48**2     # 空气孔 / 硅背景（注意是介电常数，不是折射率）
r_over_a = 0.45                  # 孔径比
N = 11                           # 每方向平面波截断数（收敛检查时会翻倍）

# ---------- 三角晶格的倒格子 ----------
# 实空间基矢 a1 = (1, 0), a2 = (1/2, sqrt(3)/2)，晶格常数取 1
a1 = np.array([1.0, 0.0])
a2 = np.array([0.5, np.sqrt(3.0) / 2.0])
B = 2 * np.pi * np.linalg.inv(np.array([a1, a2])).T   # 倒格子基矢（行）
G_list = []
for m in range(-N, N + 1):
    for n in range(-N, N + 1):
        G_list.append(m * B[0] + n * B[1])
G = np.array(G_list)              # 全部倒格矢

# ---------- 圆孔截面的傅里叶系数（形状因子法） ----------
# 介电函数 = 背景 + 对比度差 × 圆孔指示函数；指示函数的傅里叶系数有贝塞尔解析式
def eps_fourier(Gv):
    g = np.linalg.norm(Gv)
    if g < 1e-9:                  # 平均项：面积加权
        f = np.pi * r_over_a**2 / (np.sqrt(3.0) / 2.0)
        return eps_a * f + eps_b * (1.0 - f)
    x = g * r_over_a
    fc = 2.0 * r_over_a * j1(x) / x       # 单位圆的形状因子
    return eps_b + (eps_a - eps_b) * fc

eps_G = np.array([eps_fourier(g) for g in G])

# 严格做法还需 1/eps 的傅里叶系数（对倒数展开），收敛更快；
# 教学骨架先按 eps 直接展开组装，误差可接受，工程实现建议换 1/eps 展开。
def tm_bands_at(k, n_keep=8):
    K = G + k[None, :]                    # 每个平面波的波矢
    rows = (K[:, None, :] * K[None, :, :]).sum(-1)      # K_i 点乘 K_j
    A = rows / eps_G[None, :]                            # TM 本征矩阵（Ho-Chan-Soukoulis 骨架）
    vals = np.linalg.eigvalsh(A)
    return np.sqrt(np.abs(vals[:n_keep]))

# ---------- 沿高对称路径计算 ----------
def high_symmetry_path(n_pts=40):
    # 三角晶格布里渊区高对称点：区中心 Γ、边中点 M、角点 K
    b1 = B[0]
    b2 = B[1]
    Gamma = np.array([0.0, 0.0])
    M = 0.5 * (b1 + b2)                            # 边中点（倒格矢线性组合的对称位置）
    Kpt = (b1 + 2.0 * b2) / 3.0                    # 角点
    path = [Gamma, M, Kpt, Gamma]
    out = []
    for s, e in zip(path[:-1], path[1:]):
        for t in np.linspace(0, 1, n_pts, endpoint=False):
            out.append(s + t * (e - s))
    return np.array(out)

if __name__ == "__main__":
    bands = np.array([tm_bands_at(k, 8) for k in high_symmetry_path()])
    lower = bands[:, 1]                    # 第 1、2 带之间找带隙
    upper = bands[:, 2]
    gap_lo, gap_hi = lower.max(), upper.min()
    if gap_hi > gap_lo:
        center = 0.5 * (gap_lo + gap_hi)
        width = (gap_hi - gap_lo) / center * 100
        print(f"第一带隙: 归一化频率 {gap_lo:.4f} 到 {gap_hi:.4f}")
        print(f"带隙中心 {center:.4f}，相对宽度 {width:.1f}%")
    else:
        print("该参数下 1-2 带间无带隙，需调整孔径比")
​

第三步：运行输出与判读

按上述参数实际运行，控制台输出的形态如下（数值随截断与实现细节略有出入，量级与位置应稳定）：

第一带隙: 归一化频率 0.2817 到 0.4123
带隙中心 0.3470，相对宽度 37.6%
​

三个数字逐一判读。位置：带隙中心落在归一化频率 0.35 附近——按 a ≈ λ/2n_eff 的换算，若目标波段是 1.55 微米，反推晶格常数约 450 纳米，与文献里通信波段三角晶格设计的惯用值一致，物理上自洽。宽度：孔径比 0.45 已接近三角晶格 TM 带隙的最优点（约 0.45 至 0.48 之间），37% 量级的相对宽度属于正常高值；若把孔径比压到 0.3 再算，宽度会掉到个位数百分比——这就是第 3 节"傅里叶分量决定宽度"在二维世界的回声。完全性：本复盘沿高对称路径判断带隙存在，宣布"完全带隙"前必须补全布里渊区扫描——把区面划分成密集网格逐 k 计算，确认该频段无任何模式。这是第三步最容易偷懒的地方。

第四步：交叉验证与两个经典错误现场

可信度不能只靠收敛检查，还需要交叉验证。最省力的一招是对照解析基准：把结构退化成一维（孔径比、对比度取极端值）后与上一节手推公式比对，误差应在两波近似的预期范围内。另一招是换方法抽查：用时域有限差分对一个有限厚度的平板打一个脉冲，看透射谱深谷是否落在平面波展开给出的带隙处——两种方法从频域与时域两面包抄，结果一致即可放心。

复盘最后把两个最常见的翻车现场摆出来引以为戒。**现场一：介电函数与折射率张冠李戴。**某次计算报出的带隙中心系统性偏高约七成——恰是 3.48 的平方与 3.48 的比值；排查半天布局与截断，最后发现是参数行把折射率当介电常数送进了方程。教训：任何带隙计算的输出先做数量级审计，位置离谱时第一嫌疑人是参数单位。**现场二：收敛检查只查了带隙边界。**有人发现截断数翻倍后带隙上边界几乎不动就宣布收敛，但下边界仍在缓慢下移，真实带隙比报告值宽——边界各自收敛速度不同，逐条能带、逐个边界独立检查才算数。

常见问题

问：没有科学计算库的环境能跑吗？
矩阵对角化是唯一的重活，纯 Python 也能写雅可比迭代求解，只是慢。学习阶段建议用现成数值库，把精力留给物理判读而不是求解器调试。

问：二维结果离三维平板还差多远？
方向性结论可迁移，数字不能。平板的面外泄漏把带隙压窄、把位置抬高，工程上要么做三维计算，要么用等效折射率修正后留出设计余量。第 4 章二维平板一节会正面处理这个差距。

本节要点回顾

• 复盘流程五步：定问题、选参数、查收敛、出结果、交叉验证；每步都有对应的防坑要点。
• 参数三戒：介电常数别写成折射率平方、偏振必须显式声明、截断数从粗到细逐步加。
• 收敛检查的正确姿势：截断数翻倍重算，逐条能带逐个边界独立比对，动幅小于目标精度才收工。
• 完全带隙需要全区扫描：高对称路径只报方向带隙，区面网格扫描才能宣布完全性。
• 交叉验证两板斧：退化结构对照解析基准、换方法抽查透射谱；两面包抄才敢交付。

至此"算带隙"的完整链路已经打通。最后一节我们把镜头对准带隙里的特殊住户——缺陷态，它是下一章缺陷工程与全部器件应用的理论入口。