最近我在调试一维光子晶体的反射相位谱时遇到一件让我重新审视“常识”的事我把入射角从5度一路往0度收按道理光束垂直打上去就该垂直弹回来但在带隙边缘的波长附近反射光束的重心居然横向漂了将近十几个波长。垂直入射却不垂直反射这听起来像悖论但在光子晶体领域“正入射光束位移”确实是真实存在的现象而且它背后牵涉到对称性、布洛赫模式和等频面倾斜这些硬核概念。这篇文章我想把这套东西从头捋一遍先讲清楚正入射位移为什么“不应该存在”再讲光子晶体怎么把它“变出来”最后给出一套完整的一维复现流程包括传输矩阵法、高斯光束角谱分解和位移计算代码。无论你是刚接触计算光学的研究生还是正在复现文献里巨大位移结果的工程师希望这条“从理论到代码”的路径能帮你少走弯路。1. 正入射位移真的存在吗先从对称性说起1.1 光束位移是什么从Goos-Hänchen效应讲起光束位移这个现象最早被系统研究是从德国物理学家Goos和Hänchen在1947年的实验开始的。当一束光以大于临界角的入射角打到界面上发生全反射时反射光束并不在几何光学预测的那个点离开界面而是沿界面方向平移一段距离。这段位移通常是波长的量级肉眼看不见但用干涉手段完全可以测到。为什么会有这种位移简单理解全反射时虽然光在宏观上被“弹”回来了但电磁场并没有完全被拒绝在界面之外它会以倏逝波的形式渗入到低折射率介质里一段深度再被“拉”回来。这段渗入和返回的过程等效于光束在界面上多走了一段路于是反射点看起来就平移了。Artmann在同年给出了一个漂亮的稳态相位公式位移 ( \Delta ) 正比于反射系数的相位 ( \phi ) 对横向波矢 ( k_x ) 的导数[ \Delta -\frac{\partial \phi}{\partial k_x} ]这个公式是整个光束位移计算的灵魂。它的物理图像是一束有限宽度的光可以分解成一系列不同 ( k_x ) 的平面波分量每个分量的反射相位不一样合在一起后等相位面发生了倾斜而能流方向垂直于等相位面于是光束整体就横向漂移了。后面我给的数值复现本质上就是把这个公式从“解析近似”推广到“全数值积分”。1.2 对称性保护为什么正入射时位移通常为零现在回到正入射。想象一个完全均匀的平面镜光束垂直打上去。入射光的 ( k_x0 )反射光当然也应该是 ( k_x0 )。光束里那些微小的角谱分量假设有个分量略微带一点 ( k_x )那一定也存在一个对称的 ( -k_x ) 分量。对于普通的平面镜这两个分量的反射系数大小相等、相位也相等它们的横向漂移会完全抵消最终反射光束的重心稳稳地停在入射点上。用Artmann公式的语言来说结构如果在横向存在镜像对称性那么反射系数满足[ r(k_x)r(-k_x) ]相位 ( \phi(k_x) ) 是偶函数在 ( k_x0 ) 处一阶导数为零所以位移严格为零。这不是数值误差而是对称性保护的结果。这一点非常关键。很多初学者在复现文献时把一维多层膜算来算去发现正入射位移归零就以为是自己代码写错了。其实不是是结构本身限制了它。所以真正要做“正入射光束位移”第一步不是急着写代码而是先问自己我打算打破哪个对称性1.3 打破对称性的三条路线文献里真正实现正入射非零位移的机制大致可以分成三类各有各的物理源头。第一类是等频面倾斜。光子晶体内部传播的不是简单的平面波而是布洛赫波。布洛赫波的群速度方向垂直于等频面如果等频面相对界面法线是倾斜的那么即使在正入射条件下横向波矢为零能流方向也会有一个横向分量光进入晶体后就会横向漂移。宏观世界里最像这个现象的是双折射晶体里的e光偏移——光垂直入射到一块光轴倾斜的晶体板上出射光束照样平移。第二类是磁光效应。给光子晶体加一个外磁场介电常数张量会出现反对称的非对角元时间反演对称性被打破。这时候圆偏振的正负模式不再简并正入射的线偏振光反射后两个圆偏振分量会分离出现类似自旋霍尔效应的横向分离。这个方向的工作通常和“光子自旋霍尔效应”放在一起讲位移量级可以做得很可观。第三类是表面模式耦合。光子晶体截断表面如果存在Tamm态或泄漏模式反射相位在共振点附近会急剧变化。配合一个横向不对称的耦合结构比如非对称光栅单元就能在正入射下获得非常大的位移。我这篇文章的复现路线走的是“一维光子晶体带隙边缘相位增强”这条最稳的路。它严格来说算“近正入射”的增强方案正入射极限下位移会被对称性压回零但理解它的数值框架之后你想往二维等频面倾斜方向扩展只需要换色散计算位移计算的底子完全一样。2. 光子晶体里的布洛赫模式与等频面位移的“方向盘”2.1 从布拉格反射到带隙为什么要用光子晶体普通介质里的光束位移哪怕是全反射量级也就是几个波长。想获得大位移需要让反射系数对 ( k_x ) 有足够大的相位梯度。光子晶体的优势在于它可以提供非常陡峭的相位响应。一维光子晶体本质上就是多层膜高低折射率层交替排列。在布拉格条件下每层反射的光相位一致叠加后反射率可以无限接近1形成一个“光子带隙”——就像半导体带隙禁止电子通过一样这个频率范围内的光不能在晶体里传播。最妙的是带隙边缘的反射率从接近0迅速爬升到接近1伴随而来的就是反射相位在很窄的频率/角度范围内发生大角度变化。相位变化越陡位移越大。我经常拿一个比喻给组里同学讲光子晶体带隙边缘就像信号处理里的“边沿滤波器”相位响应曲线陡峭而位移本质上就是相位对波矢的斜率。滤波器越陡输出信号对频率越敏感光学里就是位移越大。2.2 布洛赫模式与等频面倾斜能流方向不等于波矢方向在光子晶体里我们处理的是周期介质中的布洛赫波。布洛赫定理告诉我们本征场可以写成振幅受周期调制的平面波[ E_k(\mathbf{r}) e^{i\mathbf{k}\cdot\mathbf{r}} u_k(\mathbf{r}) ]其中 ( u_k ) 具有和晶格相同的周期性。色散关系 ( \omega(\mathbf{k}) ) 在周期介质中不再是简单的直线或圆而是被折叠成能带结构。群速度[ \mathbf{v}g \nabla{\mathbf{k}} , \omega(\mathbf{k}) ]是能流传播的方向它总是垂直于等频面。对于各向同性均匀介质等频面是圆群速度沿径向和波矢方向一致。但光子晶体的等频面经常发生形变甚至被扭曲成斜椭圆。一旦等频面法线不再沿波矢方向就会出现“斜向能流”。这就是正入射位移在物理上最本质的来源入射光的波矢是垂直于界面的但激发出来的布洛赫模式的能量流动方向可以不垂直于界面。光进到晶体里以后像进了传送带横向漂移一段再离开。2.3 截断带来的表面态相位陡变的放大器除了能带本身的色散还有一个经常被低估的机制表面态。完整周期性晶体被截断后表面处可能出现局域在界面附近的模式——一维对应Tamm态二维和三维对应表面态或表面泄漏模式。这些模式通常位于带隙内频率上靠近带隙边缘时它们会和入射光发生耦合。表面态对反射相位的影响相当大。正入射光在表面态共振频率附近反射相位会发生接近 ( 2\pi ) 的快速跳变这时候 ( \partial\phi/\partial k_x ) 可以达到很大的值位移也就被放大数倍甚至数十倍。所以做这个方向表面终止层怎么截断非常关键。同一块光子晶体最后一层是高折射率还是低折射率厚度偏了多少都会显著改变表面态的位置和耦合强度进而改变位移符号和大小。这个后面参数扫描部分还会详细说。2.4 复现选型建议先一维后二维如果你是想快速建立“理论到复现”的完整闭环我强烈建议先做一维光子晶体多层膜。原因有三传输矩阵法TMM精确且数值稳定200行Python就能搞定。带隙和相位响应有清晰的解析预期方便校验代码。位移计算框架和高维完全一致——都是角谱分解加相位响应一旦跑通换成二维光子晶体只是把“反射系数计算器”换掉而已。等你想看真正的等频面倾斜导致的正入射位移再上二维平面波展开法PWEM或者FDTD不迟。但数值方法的核心思想一维里已经全部有了。3. 用传输矩阵法把多层膜搬进代码3.1 传输矩阵的核心逻辑传输矩阵法适合处理一维分层结构。它的思路是每层介质对应一个 ( 2\times2 ) 矩阵把层一侧的电场和磁场切向分量映射到另一侧。把所有层的矩阵按顺序相乘得到整个结构的等效矩阵再结合入射介质和出射介质的边界条件就能解出反射系数和透射系数。对于TE偏振电场垂直入射面第 ( j ) 层折射率 ( n_j )、厚度 ( d_j )、入射角 ( \theta_j )的相位厚度是[ \delta_j \frac{2\pi}{\lambda} n_j d_j \cos\theta_j ]在代码里我不会真的组装矩阵然后解方程而是用等价但更不容易出错的导纳递推法从出射介质开始逐层向前递推结构的输入导纳 ( Y )最终反射系数就是[ r \frac{q_{\mathrm{inc}} - Y_{\mathrm{in}}}{q_{\mathrm{inc}} Y_{\mathrm{in}}} ]其中 ( q_{\mathrm{inc}} n_{\mathrm{inc}} \cos\theta_{\mathrm{inc}} ) 是入射介质的TE导纳。递推公式一层层套下去逻辑非常清晰。3.2 Python实现一个精简但完整的函数下面这个函数够我们后面所有计算用了。结构定义成从空气入射、以低折射率层结束的多层膜也就是 ( \text{air} | (HL)^N | \text{air} )。import numpy as np def bragg_reflection(n_H, n_L, N, lam0, lam, theta0.0, polTE): 一维光子晶体反射系数。 结构: air | (H L)^N | air n_H, n_L: 高低折射率 N: 周期数 lam0: 设计波长, 层厚取光学厚度 lam0/4 lam: 工作波长 theta: 入射角(rad) k0 2 * np.pi / lam d_H lam0 / (4 * n_H) d_L lam0 / (4 * n_L) n_inc n_sub 1.0 # 入射介质的横向波矢, 正入射时 kx 0 kx k0 * n_inc * np.sin(theta) if pol TE: q_inc n_inc * np.cos(theta) # 出射介质 sin_theta_sub np.clip(n_inc * np.sin(theta) / n_sub, -1, 1) q_sub n_sub * np.sqrt(1 - sin_theta_sub**2) Y q_sub # 从出射一侧开始 # 注意递推方向: 从靠近出射介质的最后一层 L 开始, 所以先处理 L 再 H for _ in range(N): for n, d in ((n_L, d_L), (n_H, d_H)): cos_t np.sqrt(np.clip(1 - (kx / (k0 * n))**2, 0, 1)) q n * cos_t delta k0 * n * d * cos_t Y q * (Y 1j * q * np.tan(delta)) / (q 1j * Y * np.tan(delta)) r (q_inc - Y) / (q_inc Y) return r else: # TM偏振只需把导纳换成 1/q, 这里留作扩展 raise NotImplementedError(这里只写TE, TM的同理换导纳即可)这是个最小实现。说实话这个版本为了便于阅读牺牲了一点点效率——如果后面要做几万个 ( k_x ) 点的高斯光束角谱建议用numba或者把for循环向量化。不过对于演示完全够用。3.3 先用反射谱验证代码代码写完千万别直接算位移先画个反射谱确认带隙位置是准的。我的习惯是扫描400到2000nm波长在正入射下看反射率。参数取 ( n_H2.35 )( n_L1.45 )( N8 )设计波长 ( \lambda_01.55\mu m )。此时每层光学厚度都是 ( \lambda_0/4 )带隙中心应该在1550nm附近反射率接近1。带隙边缘大概落在1400nm和1750nm附近反射率会急剧掉下来而就在这两个掉下来的地方相位谱会出现陡坡。如果你跑出来发现带隙中心不在设计波长先检查 ( d_H ) 和 ( d_L ) 的公式。四分之一波长条件是 ( d \lambda_0/(4n) )光在介质里的波长缩短了 ( n ) 倍这个最容易写错。4. 从平面波到高斯光束位移量的数值计算4.1 为什么平面波算不出位移平面波是单一 ( k_x )反射后还是单一 ( k_x )虽然反射系数的相位会变但这个相位只是整体载波相位不会让光束的能量中心在横向上移动。位移必须靠不同 ( k_x ) 分量之间的干涉才能体现。所以计算位移之前要把入射场换成有横向分布的光束最常用的是高斯光束。高斯光束在束腰处的横向电场分布是[ E(x,0) E_0 \exp\left(-\frac{x^2}{w_0^2}\right) e^{ik_{x0}x} ]其中 ( w_0 ) 是束腰半径( k_{x0} ) 是光束整体的横向中心波矢正入射时为零。4.2 角谱分解与逐分量反射我们把高斯光束做傅里叶变换得到它在 ( k_x ) 空间的角谱分布 ( A(k_x) )。高斯光束的角谱仍然是高斯形宽度约为 ( 1/w_0 )。每个 ( k_x ) 分量对应一个入射角[ \theta \arcsin\left(\frac{k_x}{k_0 n_{\mathrm{inc}}}\right) ]用bragg_reflection函数算出这个角度下的反射系数 ( r(k_x) )反射光束的角谱就是[ A_r(k_x) A(k_x) \cdot r(k_x) ]再逆傅里叶变换回实空间就得到反射光束的横向电场分布 ( E_r(x) )。有一点实操上的提醒入射光和反射光的传播方向相反但在界面处横向波矢是连续的所以这里用傅里叶变换的逆变换就可以不需要额外处理传播方向。符号细节如果写错反射光束的重心可能会算成负的但绝对值不受影响。我的经验是先用一个纯相位响应比如 ( re^{i\alpha k_x} )自检一下看位移方向是否符合预期。4.3 数值重心法得到反射光束的电场分布后位移量就直接用能量重心之差[ \Delta \frac{\int x |E_r(x)|^2 , dx}{\int |E_r(x)|^2 , dx} ]入射高斯光束的重心在 ( x0 )所以这个 ( \Delta ) 就是反射光束相对入射点的横向位移。这里有个很容易被忽略的点入射光束的重心虽然是你自己定义的但实际计算中FFT采样网格里的坐标偏移、半像素误差都会混进重心计算里。我建议先算一束“不加结构”的反射也就是直接 ( r1 )作为基准再把位移结果减去这个基准值能有效消除数值坐标系统的固有偏移。4.4 完整计算流程代码下面这段是位移计算的完整函数参数上我留了一个theta0方便你做“近正入射扫描”或者“严格正入射”两种情况。def gaussian_beam_shift(n_H, n_L, N, lam0, lam, w0, Nx4096, dxNone, theta00.0): 计算高斯光束经一维光子晶体反射后的横向位移。 if dx is None: dx w0 / 10 # 空间采样步长, 要远小于束腰 x (np.arange(Nx) - Nx / 2) * dx k0 2 * np.pi / lam # 入射高斯光束(束腰在结构表面) E_in np.exp(-(x / w0)**2) * np.exp(1j * k0 * np.sin(theta0) * x) # 角谱 A np.fft.fftshift(np.fft.fft(np.fft.ifftshift(E_in))) kx 2 * np.pi * np.fft.fftshift(np.fft.fftfreq(Nx, dx)) # 逐 kx 分量计算反射系数 r np.zeros(Nx, dtypecomplex) for i, kxi in enumerate(kx): theta np.arcsin(np.clip(kxi / (k0 * 1.0), -1, 1)) r[i] bragg_reflection(n_H, n_L, N, lam0, lam, theta) # 反射角谱 - 反射光束 A_r A * r E_refl np.fft.ifftshift(np.fft.ifft(np.fft.fftshift(A_r))) # 重心位移 I_r np.abs(E_refl)**2 shift np.sum(x * I_r) / np.sum(I_r) return shift, x, E_in, E_refl这个函数跑一遍如果参数取 ( N8 )( w_020\lambda )在带隙边缘波长附近扫描入射角会看到位移随角度迅速增大在接近正入射时却回落。这不是代码错了而是对称性在起作用。想观察增强效应可以把入射角固定在比如 ( 1^\circ )然后去扫波长在带隙边缘附近你会看到位移出现一个明显的尖峰。5. 参数扫描位移的旋钮到底在哪里5.1 周期数 ( N )相位陡度的第一旋钮光子晶体周期数越多带隙越“硬”反射率在带隙内越接近1带隙边缘的相位变化也越陡峭。我用N4, 8, 12分别跑过在 ( \theta1^\circ )、波长取带隙边缘附近时位移大约是周期数 N带隙边缘位移示意值4约 ( 2\lambda )8约 ( 12\lambda )12约 ( 35\lambda )注意这组数字是个定性趋势具体值取决于折射率对比度和波长选择。但趋势是可靠的层数越多位移越大。代价是带宽变窄可调节的波长窗口变小而且结构对加工的精度要求更高。5.2 波长位置必须“钉”在带隙边缘带隙中心附近反射率虽然高但相位变化平缓位移很小带隙深处光是强度耦合不进去的相位响应近乎平稳只有在带隙边缘反射率处于“将掉未掉”的状态相位梯度最大位移也最大。实际操作中我习惯把波长步长取到带隙宽度的1/50以内去扫描。如果步长太大很可能直接错过位移尖峰误以为这个结构效果不行。带隙边缘是这种效应最敏感的区域值得你多花点算力去“蹲守”。5.3 入射角近正入射增强的极限位移对角度的依赖很有意思。在带隙边缘位移既可以在某个非零角度下达到峰值又在接近0度时被压回零。峰值角度通常在带隙边缘的“泄漏锥”附近和我上面说的表面态耦合有关。你的结构表面终止方式变了峰值角度也会跟着变。想往严格正入射推就得换机制——比如引入等频面倾斜的二维结构让 ( \partial\phi/\partial k_x ) 在 ( k_x0 ) 处不为零。我在文末会再提一句。5.4 束腰 ( w_0 )一个最容易被忽略的陷阱光束束腰对位移结果的影响很多人第一次算都会栽跟头。束腰太小角谱铺得太宽一部分角谱分量落在带隙外位移被“掺水”摊平束腰太大角谱太窄理论上也能算出接近解析值的结果但数值上角谱采样点不够FFT的 ( k_x ) 分辨率不足计算结果会抖动很厉害。我用的经验法则是束腰不小于 ( 10\lambda )否则角谱展宽太严重。束腰不大于 ( 100\lambda )否则要非常大的 ( N_x ) 才能保证角谱分辨率。空间采样步长 ( dx ) 至少要小于 ( \lambda/(10 n_{\mathrm{max}}) )才能容纳介质内部的场变化。如果你发现位移随束腰变化超过10%那基本是角谱采样或者空间采样的问题先把网格加密再说。6. 复现过程中的实际经验与坑6.1 相位解包arctan的2π跳变会让你怀疑人生如果你直接用np.angle(r)去取相位然后试图数值求导大概率会得到一堆锯齿状的东西——因为np.angle把相位卷绕到 ( (-\pi,\pi] ) 区间每跨过 ( 2\pi ) 就跳一次。位移计算用到了 ( \partial\phi/\partial k_x )相位不连续会导致位移出现虚假的尖刺。两条路可以选一是对原始复数反射系数逐点计算位移不用显式求相位我上面给的代码就是这么干的二是用np.unwrap()先把相位解开再求导。前者更稳因为它在复数域直接做乘法不引入任何人工处理。6.2 角谱采样范围不足位移算出来是振荡的有一次我把 ( N_x ) 从1024改成4096之后位移突然从乱跳变成了平缓曲线才知道之前算的振荡全是采样不足造成的假象。判断标准很直接固定其他参数只增加 ( N_x )如果位移曲线还在明显变那就是没收敛。角谱的 ( k_x ) 范围是 ( 2\pi/(2dx) )你要保证这个范围至少覆盖高斯角谱的5个标准差以上也就是 ( 5/w_0 )。另外( k_x ) 的分辨率是 ( 2\pi/(N_x dx) )要保证在一个带隙的特征角宽度内有足够多的采样点。6.3 表面终止层厚度微小偏差改变位移符号这是我从一个“意外的负位移”里学到的。当时我微调了表面高折射率层厚度从标准的 ( \lambda/4 ) 改成 ( 0.9\lambda/4 )位移峰值不仅大小变了符号还反转了。原因在于表面终止层决定了反射相位谱的形状尤其是表面模式的位置。表面层厚度变化时表面态共振频率移动相位陡变区域移向不同角度/波长位移的峰值和方向自然跟着变。所以如果你的目标是“复现文献结果”第一件事就是把文献的结构参数特别是表面终止层的厚度精确到小数点后两位拿过来用。哪怕差5%结果都可能对不上。6.4 如何验证你的数字是靠谱的我自己的验证习惯分三步解析极限验证在小角度极限下位移应该趋近Artmann公式的预测值。跑几个不同角度把数值位移和 ( -\partial\phi/\partial k_x ) 解析值画在一起两者应该重合。结构退化验证让高低折射率相等( n_Hn_L )光子晶体退化成均匀介质位移应该归零。如果这时候不为零你的代码有bug。束腰收敛验证连续增大 ( w_0 )位移应该逐渐逼近“无穷大束腰极限”——这个极限对应的是完整角谱加权积分而不是单点稳态相位近似。三步都过了我才敢说这个位移结果是可信的。如果你手头有FDTD工具可以把一维结果和FDTD对一遍但FDTD网格色散会带来微小的频移要对齐反射谱的带隙边缘来比别直接对同一个波长。最后再分享一个小经验正入射光束位移这个题目初看像个“计算物理练习题”真正踩进去才发现它对物理直觉的要求很高。我这段时间最深的体会是计算框架本身并不难难的是知道自己算出来的这个位移对应的是哪个机制。同样是位移尖峰到底是带隙边缘相位陡变贡献的还是表面态共振贡献的还是等频面倾斜贡献的这三者的参数依赖关系完全不一样。所以建议大家做参数扫描的时候不要只盯着位移曲线看。把反射谱、相位谱、角谱分布一起画出来三张图对着看。一旦发现位移峰的位置恰好和相位梯度峰重合而反射率又处于过渡区基本就能锁定机制了。这套代码和思路接下来还能往两个方向扩展一个是把TE换成TM对比偏振行为另一个是换成二维光子晶体的平面波展开法算真正的等频面倾斜正入射位移。希望这篇从理论到复现的记录能让你在自己的项目里少踩几个坑。