做光子晶体或者光力学方向的人应该都见过这两年冒出来的各种摩尔光子学工作——两层一样的介质柱阵列叠在一起转一个几度的小角度原本平平无奇的能带会突然变密集出现平带、局域模式甚至拓扑性质大变。我最近折腾的项目就是这个方向再往前推一步不仅看双层扭曲光子晶体里辐射场的拓扑性质还专门研究这个场的分布怎么决定光力。项目名可以概括为“双层扭曲光子晶体辐射场拓扑光力”。这篇就把整个项目的物理图像、建模思路、计算方法和踩过的坑都拆开写一遍给想入坑这个方向的人一个可以直接参考的路线图。1. 项目定位与整体设计思路1.1 从双层扭曲到摩尔光子学为什么又双叒是“转角度”扭曲双层结构的灵感最早来自石墨烯。两个六角晶格叠在一起相对转过一个很小的角度比如 1.1 度左右就会形成一个波长极大的摩尔超晶格。电子在这个超周期势里感受到的相互作用被强烈放大和完全没转角的双层比起来电子关联效应完全不同。这个物理在电子系统里被反复验证之后光学方向的人很快就意识到光子晶体里的介质柱阵列本质上也是一个周期介电函数完全可以把“扭转”这个自由度搬进光子晶体。两层二维光子晶体只要材料是低损耗介质电场传播基本由介电常数分布决定不像电子系统那么复杂。但正因为它不需要考虑复杂的库仑相互作用光子晶体里的摩尔结构往往更容易在实验上搭出来也更容易用数值方法算清楚。实际构建时通常的做法是取两片同参数的晶格一片固定作为底层另一片绕中心旋转一个小角度再把它们放在近场距离内。两个周期性介电函数叠加在一起之后整体的介电调制的包络周期远大于原晶格常数于是布洛赫定理描述的快速振荡波包上多了一个缓慢变化的“超胞势场”。这个超胞势场正是整个项目的核心舞台。它可以产生很多单层甚至双层对齐结构里看不到的现象能带平化、慢光、模式局域、角动量重整化等。我们做拓扑光力看中的就是摩尔超晶格能把普通的能带结构打散重组把原本弥散的辐射场压缩在某些空间区域而辐射场的压缩、涡旋、单向传播都会直接改变光力在结构上的分布。1.2 为什么盯着“辐射场”不放光子晶体里通常最吸引人的是束缚模式也就是局域在结构内部、不会往外跑的那些本征模。束缚模式用来做波导、谐振腔、滤波都很好但它们和外界环境交换能量很弱做传感或者机械驱动时效率反而受限。辐射场则相反它是开放系统里能往外泄漏电磁能量的模式即使一个结构没有任何输入也存在自然的辐射通道可以耦合成出射波。那辐射场和光力有什么关系宏观上光压是光照射物体时动量转移的结果辐射场既然不断有光子跑进跑出动量转移在方向和大小上就丰富得多。尤其是带拓扑结构的辐射场它的动量流可以做得很复杂包括横向自旋霍尔效应、涡旋相位、结构手性等。这些性质投射到力学量上就会产生普通平面波照射无法实现的力分量比如局域扭矩、指向性的散射力、鲁棒性很弱的横向力都可以设计出来。我看重辐射场的另一个原因是实验可测度高。束缚模式要耦合进探针往往需要额外加棱镜或者近场探针过程繁琐且容易扰动原结构。辐射场模式天然在外部空间有可探测的光功率用远场探测器、焦平面测量或者数值孔径收集就能反演内部场的对称性和相位信息。对拓扑光力来说能用实验验证模型是后续能否形成完整闭环的关键。1.3 把“拓扑”和“光力”捏在一起的价值单谈拓扑光子学大家第一印象是边缘态、单向传播、背散射抑制这是通信的视角。单谈光力大家第一印象是光镊、悬浮光力学、引力波探测这是精密测量和量子操控的视角。双层扭曲光子晶体辐射场拓扑光力的思路是把两块拼图合成一个完整的研究对象用拓扑结构决定了辐射场的相位与空间对称性再把这些性质转化为可控的力学响应的阵列分布。好处很明显。拓扑保护的传播模式不害怕结构微扰如果结构表面有一点粗糙、有一点误差光场模式和对应的光力分布依然能保持形状普通光子晶体只要碰上一个杂质或者缺陷局域场就变了光力也跟着乱套。拓扑结构还把自由度分得很清楚比如某个能带携带非零拓扑荷那么它辐射出来的场在动量空间就会有涡旋转矩方向就锁定了。这种输出“带有鲁棒性”的光力是以后做集成光机械逻辑、粒子筛选、甚至微型马达驱动都很值得借力的点。2. 双层扭曲光子晶体的物理图像与数值建模2.1 扭转角的核心角色从失配晶格到摩尔周期摩尔超晶格的形成本质是失配。两层晶格的晶格常数一样但坐标方向相对转了角度于是两套倒格子方向不一致。真实空间的超胞周期是典型摩尔干涉的结果理论计算常用公式A_M a / (2 sin(θ/2))其中 a 是原晶格常数θ 是层间扭转角。小角度时这个公式可以近似成 A_M ≈ a / θ。简单说1 度的扭转角配上 500 nm 的晶格常数摩尔周期大约 28.65 μm这个尺度远远超出了单独一层光子晶体的周期也让超胞内可以容纳非常丰富的衍射阶次。需要小心的是这个公式只在两层结构都严格是无损周期晶格时才精确成立。实际建模时两层之间往往有真空层或低折射率间隔层摩尔周期会轻微偏移有限尺寸的样品还会引入边缘效应周期不再是无限延展。更麻烦的是光子晶体的介质柱如果是有限高度扭转面位置本身也会改变局部模式的耦合强度。因此数值模型里不能只设一个“完美至极的旋转复制”最好先扫描层间距和扭转角两个参数让超胞周期与目标摩尔周期对齐之后再谈能带和光力。2.2 材料与结构选型介质柱还是空气孔怎么选双层扭曲光子晶体在材料选择上没有标准答案但要同时做辐射场和光力研究我有几个建议方向。第一折射率对比度要足够大。氮化硅、硅、砷化镓这些常规光子晶体材料都行硅的折射率约 3.48对比空气衬底带隙宽、模式局域能力强对辐射场的束缚也更好。不过硅在近红外有吸收长波通讯窗口虽然低损耗但做光力实验时吸收热不可避免如果要做高灵敏度的机械响应反而要考虑色散低、损耗更小的碳化硅或者氮化铝。第二结构类型上介质柱型光子晶体的模式场主要集中在低折射率空气区域带隙里的慢光模式很适合做梯度力增强空气孔型光子晶体的高折射率背景连续和平板波导耦合方便但辐射损耗相对难控制。我第一版模型用的是硅介质柱三角晶格双层之间用一层二氧化硅做间隔。这样做的好处是两层晶格可以通过同一套沉积和刻蚀工艺完成转角度也可以通过电子束光刻对准实现数值模型的周期边界也容易设置不必费劲处理悬浮结构。2.3 数值模型的建立要点数值模型里最关键的是摩尔超胞到底怎么建。由于整体几何不再是单一周期结构常见做法是把两层晶格同时放进同一个“超胞”里超胞尺寸等于摩尔周期 A_M再对整个区域施加布洛赫周期性边界条件。这个过程说起来容易真操作时有个隐藏麻烦如果扭转角不是能被 60 度整除的某个特殊角度两层结构的相对位置不会在有限超胞内自闭合边界的原子/柱子分布会产生截断伪影。所以很多人的做法不是直接扭转两个大晶格而是先构造一个周期性的摩尔超胞内部包含第一层完整晶格和第二层完整晶格然后用扫描参数微调扭转角使得超胞边界处的柱子位置几乎对齐。这本质上是一个几何优化问题需要你在脚本里对两个原胞做旋转平移的数值匹配。我用的方案是把扭转角固定在某几个“魔角”附近比如 1.05 度、1.47 度、2.01 度分别对应整数个摩尔周期可以比较干净地跑出能带。网格分辨率也需要注意。摩尔超胞很大而原晶格尺度又很小这意味着计算区域内必须同时分辨亚波长的柱子细节和远超波长的整体周期。COMSOL 或者 FDTD 软件里面如果直接用统一网格内存需求会恐怖到无法启动。我的经验是采用局部加密网格在介质柱边界和层间界面区域加密其余背景区域用梢稀网格。辐射场分析还要在结构外设置足够大的完美匹配层吸收外行波避免边界反射污染远场相位。3. 辐射场拓扑性质的解构手法3.1 能带结构与拓扑不变量先判平庸不平庸判断辐射场是不是拓扑非平庸经典做法还是先算光子能带图。如果你能在带隙里看到一条连接上下能带的边缘态色散说明体系具有拓扑保护的边缘传播通道如果这个边缘态直接从辐射连续谱中间穿过会产生强方向性的表面波辐射这种场当然是非平庸的。更数学一点每一条能带都可以算布洛赫态的贝里联络沿着二维布里渊区闭合路径积分得到贝里相位。若整个价带群累积的陈数不为零就意味着能带之间存在能级交叉系统至少在最低能级意义上携带了非零拓扑荷。这个拓扑荷和电子系统里的陈数意义完全一样只是自变量从晶体动量变成了光子动量。但辐射场的分析不能只看能带。开放系统中本征频率是复数虚部对应损耗拓扑不变量必须推广到复波矢空间或者扩展布里渊区才能处理“带隙里的实频模式”与“连续谱里的泄漏模式”交错的问题。这一步很容易把初学者搞晕实际上你不需要真的去推导完整的非厄米拓扑理论只需要在 FDTD 结果里提取复频率的实部和虚部看带隙处的边缘态频率虚部是否显著小于体态泄漏率就能定性判断可用的拓扑保护窗口。3.2 近远场变换与涡旋辐射场拓扑的“指纹”辐射场最直接的拓扑表达方式是远场偏振或相位的涡旋结构。做一个近场到远场的变换比如用 k 空间的角谱展开你能看到每个瑞利衍射阶次上光场的振幅和相位都有确定分布。拓扑光力真正关心的不是单一方向出射平面波而是出射波在动量空间携载的相位奇点。你想象一下在动量空间中某个偏振分量绕着一个中心点转一圈相位改变了 2π 的整数倍。这个中心点就是相位涡旋它的绕数就是拓扑荷。如果围绕这个涡旋的等相位线是逆时针转拓扑荷为 1顺时针转则为 -1。这种涡旋的中心往往是辐射场完全暗点也就是远场强度为零的位置。物理上它也对应一个奇点围绕它的角动量密度不为零产生的光力会出现切向分量给结构一个旋转倾向。这是普通平面波或者高斯光束照射光子晶体完全没有的现象。分析这类涡旋我建议直接在 k 空间画两个偏振分量的相位图用强度掩膜过滤掉噪声区域之后用相位跳跃检测算法定位奇点。计算时注意相位缠绕unwrap 之前对矩阵做至少 2 倍冗余采样否则奇点识别会因空间分辨率不足而失败。3.3 BIC 与拓扑保护的绑定效应光子晶体辐射场里还有一类特殊现象叫连续谱束缚态BIC就是模式虽然和外部辐射通道频率重叠但因为对称性或者参数调制的保护不会泄漏能量。BIC 在能带图上表现为无限大 Q 因子实际对应远场涡旋的拓扑荷会落在贝里丘上。严格说每个 BIC 都是一个动量的相位奇点源/汇取决于拓扑荷的正负号。在双层扭曲结构里BIC 的存在与六角晶格的 C6v 对称性保护密切相关。两层相对转角只要保持一致对称性BIC 就能稳定存在于整个布里渊区的高对称点上。一旦引入不对称扰动BIC 变成准 BICQ 因子断崖式下降但在拓扑荷上仍保留涡旋痕迹。对光力来说BIC 附近是理想工作点场局域强度极高但和自由空间的耦合又可以通过改变扰动大小精细调控。这意味着你不用像普通光镊那样靠强功率来获取大光力而是在低输入功率下利用高 Q 值增强的局域场来放大梯度和力密度。这是双层扭曲体系里值得死死抓住的一条路线也是设计光力传感器可以发力的核心窗口。4. 拓扑光力的计算与调控4.1 光力计算的物理基础光力在宏观上可以描述为电磁场对物质的力学作用微观上则是洛伦兹力或麦克斯韦应力张量在物体表面的积分。要计算辐射场的光力分布经典公式是取麦克斯韦应力张量T_ij ε0 ε(E_i E_j - 1/2 δ_ij E^2) μ0 μ(H_i H_j - 1/2 δ_ij H^2)然后对包含目标的封闭面积分F_i ∮ T_ij n_j dA这里 n_j 是封闭面积外法向量。时间平均之后实部对应平均力虚部对应周期振荡的力分量。如果你把积分面积从目标表面外移到离结构几个波长的位置远场的辐射反力也自动包在里面这就是辐射场对结构反作用力的来源。在这个积分表达式里电磁场的相位分布非常重要。只有强度而丢失相位的情况下算不出稳态光力因此仿真后处理里必须保留 E 和 H 场的复数形式不能用幅值图当输入。我的实际脚本是先从 FDTD 保存的 monitor 点读入复数 E、H 场然后在后处理里构造 T 张量再用数值面积分算力。积分面的形状不用特别讲究长方体或者球形都行但必须完全包围目标结构且避开源区和吸收边界。4.2 拓扑结构如何塑造光力分布不同拓扑模式对光力的贡献机制差异巨大我按常见的几个类型梳理一下。一是边缘态主导的光力。拓扑边缘态在真实空间里集中在结构边界它的动量流沿边界定向传播。应力张量积分之后会得到一个横向驱动分量作用于边界媒质。这类光力的价值在于哪怕边界有缺陷回环光场绕开缺陷继续保持传播方向力的分布也基本锁定在边界从而得到一个鲁棒性的机械推动。二是平带模式增强的梯度力。摩尔超晶格很容易在布里渊区中心制造出局域平带群的传播速度近零场局域强烈。光子在这里几乎没有群速但场的梯度非常大光力密度可比普通传播模式高两个数量级以上。这种设计适合做低输入功率的光机械谐振器用很弱的信号就能驱动纳米机械结构振动。三是拓扑荷产生的自旋角动量和轨道角动量。辐射场的涡旋结构对应非零角动量密度结构表面会产生扭矩项。扭矩不等于力但它在积分里对应应力张量的非对称部分在扭矩传感和粒子旋转操控场景很敏感。四是手性耦合带来的方向选择性。如果体系带非零陈数边缘态的自旋动量锁定向左还是向右传播由光场的圆偏振方向决定。对光力而言这代表你可以用偏振控制机械运动方向不用重调结构参数。4.3 关键参数对光力的调控窗口调控光力最直接的两个旋钮是扭转角和层间距。扭转角改变摩尔超胞周期进而改变能带平带位置与 BIC 频率层间距则直接调控两层之间的模式耦合强度间距减小层间有效折射率变化剧烈场会向下层挤压梯度力增大。举个例子我用 COMSOL 扫过两层硅介质柱晶格、层间距从 50 nm 变到 400 nm 的情况。间距 400 nm 时层间耦合几乎可忽略能带结构基本上就是两个独立能带的直接叠加光力分布也毫无亮点间距降到 120 nm 以下时能带出现明显的劈裂和平化边缘态出现BIC 频率附近的 Q 值也从几百跳升至几十万。这时光力分布出现两个高峰一个位于两层之间的高折射率区域一个位于边缘通道最大值增大了接近两个数量级。光力的频率选择性也不能忽视。BIC 附近频点的力和边带上的力差了两个量级因此在实验上锁定到超高 Q 值频率点需要配合激光器的波长扫描先用实际结构的 Q 因子实测定位再锁定光力测试。这个实验流程里的技巧是不要只做单频激励最好扫一小段频率把光力响应谱记录下来找峰值和旋向反转的位置那里往往就是拓扑涡旋的奇点位置。5. 实操记录、常见问题与排查表格5.1 搭建仿真流程的实操步骤我习惯用两步走的流程第一步用平面波展开法快速扫能带结构和拓扑荷第二步用 FDTD 或者 FEM 算近场场分布和光力。平面波展开那步非常快推荐用 Python 的 scikit-finite-difference 或者 MIT Photonic BandsMPB扫一个二维模型。把两层结构建模为同一个超胞晶格常数设为摩尔周期设置相对介电常数的空间分布扫描波矢路径 G-M-K-G输出能带和电磁场分量。粗扫时每个 k 点用 64×64 平面波就可以看到平带轮廓精细计算再提到 128×128 以上注意周期性边界条件必须包含两层晶格的旋转信息。FDTD 的部分用 Lumerical 或者 MEEP 都可以。有一点大家容易忽略扭曲结构对网格的敏感度远高于普通光子晶体。因为两层的柱子边界在超胞里形成丰富的干涉图案网格稍粗一点介电常数边界就会在数值上变模糊BIC 频率漂移几百个纳米都是可能的。我最后固定用 20 nm 网格在柱边界局部加密到 3 nm模拟时间设到 4000 飞秒以上确保泄漏模式的场衰减到噪声平台以下。存场数据时至少存两个监听面一个在结构内部截取双层之间的近场一个在结构外 5 个波长处截取远场。远场监听面要用复场数据这样才可以做相位展开和涡旋定位。如果软件只存了 rms 值后面拓扑荷算法就用不了那整个项目基本白做。5.2 常见问题与排查速查表我把这段时间踩过的坑整理成了一张表按项目优先级排列问题现象可能原因排查与解决办法能带图里没有预期平带扭转角公式给错超胞周期失配用 a/sin(θ/2) 验证临界角直接打印两层柱坐标检查摩尔轮廓BIC 处的 Q 值只有几千没有上升趋势PML 吸收不够强或网格太脆增大 PML 层数网格细化到预期分辨率的一半重新验证光力计算值为零或方向混乱场数据缺相位信息积分面太靠近源区域切换为复数场输出把积分面外移两个波长以上拓扑荷扫描结果出现成对正负电荷网格采样不足或相位缠绕未解开包裹相位图做冗余采样检查采样步长是否小于 π/6边缘态光力实测方向反转偏振锁定方向预测与实际结构手性相反检查结构镜像是否把陈数符号翻转用相反圆偏振光重测双层间距参数实际为 0 后被建模软件合并成一个实体几何布尔操作没有正确设置层间隔在介质柱之间加一层 1 nm 的低折射率间隔层做强制区分5.3 后处理脚本的通用骨架拓扑光力计算的后处理脚本并不复杂但要注意索引和单位换算。我下面给一个简化的 Python 伪代码说明主要思路import numpy as np # 假设已经从 FDTD 导出了 grid_x, grid_y, E_complex, H_complex # E_complex 形状: (Nx, Ny, 3, Nfreq) def stress_tensor_time_average(E, H, eps, mu): 输入复数场返回时间平均麦克斯韦应力张量 (Nx, Ny, 3, 3) T np.zeros((E.shape[0], E.shape[1], 3, 3), dtypenp.float64) ex, ey, ez E[..., 0], E[..., 1], E[..., 2] hx, hy, hz H[..., 0], H[..., 1], H[..., 2] E2 0.5 * (np.abs(ex)**2 np.abs(ey)**2 np.abs(ez)**2) H2 0.5 * (np.abs(hx)**2 np.abs(hy)**2 np.abs(hz)**2) # 对角线项简化为电场和磁场的贡献 for i in range(3): T[..., i, i] 0.5 * eps * E2 0.5 * mu * H2 # 非对角线项来自不同场分量共轭积 # 省略展开实际需按 Re(E_i * conj(E_j)) 计算 return T def compute_force_from_surface(T, normals, area_weights): 输入应力张量和法向量沿封闭面积分得到力 F np.zeros(3) for i in range(3): for j in range(3): F[i] np.sum(T[..., i, j] * normals[..., j] * area_weights) return F这段代码的核心是用复数场而不是幅值场同时把应力张量积分落到网格离散点上。实际操作时注意力的单位是牛顿每单位深度二维模型要找实验对标的话得乘上器件厚度。5.4 给后来者的三个建议第一不要一上来就搞大扭转角。角度越小摩尔周期越大超胞越大计算量爆炸式增长。建议先用 1.2 度或者 1.5 度这种适中角度跑通整个流程确认能带、拓扑荷、光力三者的对应关系之后再转向小角度极限。第二拓扑荷计算要随着角度变化反复验证。小扭转角的摩尔超胞非常大布里渊区很小相位采样点只要少一点涡旋定位就会出现假正负电荷你在那调半天光力方向其实问题在网格采样上。第三如果实验平台允许建议在同一个样品上先做线性光学表征透射谱、远场强度分布再上光力测试。这样拓扑荷的远场指纹和光力峰值频率就能互相印证不用靠仿真预测硬猜实验窗口。这个方向目前还在快速成长期仿真先行是对的但结构建模和光力计算之间很容易脱节——很多论文只给能带和场图光力分析一笔带过真正把两者做到定量关联的很少。我在实操中最大的体会是拓扑光力不是一个靠“算一个数”就能出结果的课题它要求你对电磁场相位、动量流分布和材料机械响应都有直觉而这种直觉只能靠大量扫参数、看场图、盯积分值慢慢攒出来。建议先从单层体系把光力算法跑得滚瓜烂熟再上双层扭曲结构一步步来比直接挑战最强复杂度稳得多。