简介这是一份面向光子晶体与微纳光学研究者的COMSOL仿真与理论分析资料包聚焦连续域束缚态BIC的远场偏振特性及Q值能带计算同时覆盖k空间模拟与Matlab脚本实现适合具备一定电磁仿真基础、正在开展BIC相关课题的高校师生或工程师使用。资源包共15个文件大小约2.83MB类型包括docx文献笔记、html技术分析页面、jpg结果截图与txt说明文档docx和html用于梳理理论推导与仿真思路jpg直观展示远场偏振和能带分布整体结构清晰便于按需查阅。目前已有41人学习下载。内容既包含基于COMSOL的连续域束缚态远场偏振仿真模型也附有Matlab脚本用于k空间与Q值提取同时整理了相关文献与深度探讨可帮助读者快速复现计算流程、理解BIC物理图像并进一步拓展到更复杂的光子晶体能带结构研究。1. 光子晶体的 BIC为什么远场偏振比色散关系更重要看到“连续域束缚态”这个词大部分人的第一反应是能带曲线上那个尖锐的 Fano 共振峰或者一个理论上可以做到无穷大 Q 值的模式。但真正做器件的人会告诉你BIC 的价值从来不在 Q 值本身而在于它强迫远场辐射通道发生拓扑性关闭这种关闭会在 k 空间留下一个明确的偏振奇点。如果你只把 Q 值从 COMSOL 里拉出来画个趋势线那等于只拿到了结果的一半。另一半也就是远场偏振的演化轨迹才是决定这个模式能否应用在涡旋激光器或非线性频率转换里的关键。实际上从近场模式到远场偏振中间隔着一次不平凡的矢量场投影投影的相位参考没定对后面算偏振椭圆率、拓扑荷都是错上加错。下面我把基于 COMSOL 光子晶体仿真和 Matlab 脚本的完整流程拆开讲先解决能带和 Q 值怎么算再给出远场偏振的具体提取方式。2. COMSOL 建模篇构建单元、布洛赫边界与特征频率搜索2.1 周期性晶格与仿真单元矩形设置要进行能带计算首先要确定模型维度。很多初学者习惯直接建一个三维全波模型把单元结构、衬底、上层空气一次全建出来。这样虽然能求但代价是每求解一个特征值都要花费数十秒扫描一条完整的高对称路径通常需要二十到四十个点累计时间难以接受。实际做光子晶体仿真时更标准做法是用二维模型配合面外波矢来等效或者建立三维薄板模型但把 x 和 y 方向的边界设置成周期性边界条件。这里以最常见的介质柱型光子晶体板为例正方晶格晶格常数a 800 nm介质柱半径r 120 nm板厚度t 220 nm。在 COMSOL 的射频模块或波动光学模块下使用特征频率研究。几何上只需建立单个电介质柱和周围空气域柱体上方留出空气层底部可以暂时不设衬底。材料先设置成无损耗介质折射率n 3.5。必须强调使用无损耗材料是能带计算的一个核心前提因为只有损耗为零特征频率的虚部才能纯粹反映辐射损耗从而代入 Q 值公式。如果一开始就加入吸收损耗虚部里混杂了材料吸收你提取出的 Q 值就不是 BIC 本身的辐射 Q后面的分析全部失去意义。几何构建完成后需要设定周期性条件。COMSOL 中有两种常见方式一是直接使用“周期性条件”功能中的 Floquet 周期边界二是通过“端口”配合扫频。能带计算应当使用前者因为特征频率研究和 Floquet 边界能直接给出复数特征值。这里要特别检查一下源边界和目标边界的对应关系以及相位因子的方向变量是否一致。周期边界条件的相位因子本质是exp(i*k*(r - r_source))如果源和目标设反了k 的符号就会反转最终画出来的能带沿 k 方向是镜像错误的。这个问题在视觉上不容易发现因为能带看起来仍是连续的只有当你拿 Γ 点模式与已知文献对比时才会发觉频率对不上。2.2 设置 Floquet 周期边界与 k 空间扫描在 COMSOL 参数表中定义扫描变量这一步直接决定你能拿到什么样的 k 空间能带数据。建议建立参数kx和ky它们作为晶格倒空间中的坐标。在特征频率研究的参数化扫描中把kx和ky同时作为扫描维度。# COMSOL 参数表配置说明 # lambda_0 1.55[um] # 工作波长 # a 0.8[um] # 晶格周期 # r 0.12[um] # 介质柱半径 # h 0.22[um] # 介质板厚度 # kx 0 # 倒空间 x 分量由扫描定义 # ky 0 # 倒空间 y 分量此处kx的单位在 COMSOL 中需要根据模型几何使用的单位保持一致通常直接给1/m或者无量纲归一化值。扫面路径建议沿高对称线进行也就是Γ - X - M - Γ。在Γ - X段固定ky 0逐步增加kx每步约取a / 40作为步长。到了X - M段kx固定在倒空间边界值ky从 0 扫到边界值。最后M - Γ段则让kx和ky同时减小。这里每个点的特征频率求解都是在上一步的解基础上作初始值所以连续性很好不会出现模式跳变。一个比较关键的设置是特征频率研究里的“期望模式数”。光子晶体板在目标频率附近通常会存在多个模式包括横电模式和横磁模式。如果只填一个求解器会随机返回一个填太多又会增加大量无效模式求解时间。经验上在目标频段 0.7 到 1.2 倍之间模式数 6 到 8 个就足够覆盖感兴趣的两个能带。求解后按实部频率排序通过查看电场分布图判断模式是类横电还是类横磁。此外还应开启“搜索该值附近”选项填入一个初始猜测频率如c / (a*n_eff)附近的值能显著加快收敛。2.3 提取 Q 值与能带数据特征频率扫描得到的是一个复数频率f Re(f) j*Im(f)。其中实部是模式频率虚部代表辐射损耗和材料损耗的总和。在无耗材料下Q 值的计算公式为Q Re(f) / (2 * abs(Im(f)))当扫描点非常靠近 Γ 点或某些高对称点时BIC 的辐射损耗理论上趋向于 0此时虚部基本由数值噪声决定频率虚部可能在1e-4到1e-6之间无规律跳动。把这个真实值直接代入计算出的 Q 值会在1e5到1e8之间剧烈震荡绘制在对数坐标轴上就是一条竖线乱刺的曲线完全没法看。在实际项目处理中我一般是在后处理时对 Q 做一次截断当Q 10^7时将其值赋为10^7并加一个标记。这个处理并非掩盖物理事实而是绘图需要因为真正无限大的 Q 值在数值结果里不可能出现截断上限能保留趋势又不破坏标度。如果要做收敛性验证或发表图完整虚部数据还是要保留截断只用于趋势图。数据导出建议使用 COMSOL 的“全局计算”功能把kx、ky、Re(f)、Im(f)一次性导出到文本文件。导出时需要勾选“每个扫描参数作为单独列”这样后续用读取工具解析时k 路径的重建会方便很多。3. 远场偏振计算球面上坐标定向投影与相位奇点3.1 为什么要从近场算远场近场分布看起来是电磁场在结构周围的振荡形态比如介质柱内部电场增强或者表面电场束缚。但BIC的远场特性不同于普通模式它的偏振态并不是直接看几个切面电场分量就能确定的。从近场到远场需要借助 Stratton-Chu 公式做一次外推把边界上的切向电场和磁场积分到无穷远球面上。COMSOL 的远场计算逻辑正是如此它把计算域边界上的场投影到设定的远场方向上得到一组复振幅。关键点是投影球面的半径必须处于远场条件之下同时模型外部必须添加完美匹配层。否则边界反射会污染远场幅值导致偏振椭圆取向完全失真的。我在第一次做一个类似结构时省了完美匹配层只用散射边界结果远场Ex和Ey分量的幅值出现明显的干涉波纹但当时误以为是物理上的偏振不均匀。后来加了完美匹配层并留出至少lambda / 2的间隔结果立刻干净了。这里建议使用 COMSOL 内置的球形完美匹配层域并将远场边界设置在完美匹配层内表面外一个波长处。还有一个容易被忽略的问题周期结构的远场不是连续的球面分布而是分成离散的衍射级次。在计算远场之前需要去查阅 COMSOL 的远场设置中是否有针对周期性结构的衍射级次选项。对于亚波长晶格我们研究的就是零阶衍射高于零阶的都处于倏逝波范围。如果设置了错误的衍射级次得到的远场偏振就不再是 BIC 对应辐射通道的真实表现。3.2 偏振椭圆率与偏振角计算在拿到远场复振幅之后计算偏振可以按照以下步骤进行。假设导出的两个正交切向分量为Ex_cos、Ex_sin、Ey_cos、Ey_sin重组为复数形式import numpy as np # 从COMSOL导出的远场分量分别为实部和虚部 Ex Ex_cos 1j * Ex_sin Ey Ey_cos 1j * Ey_sin # 计算两个分量之间的相位差 delta np.angle(Ey) - np.angle(Ex) # 计算偏振椭圆长轴方位角 psi # 需要处理 cos(delta) 符号和分母为零的情况 denom np.abs(Ex)**2 - np.abs(Ey)**2 psi 0.5 * np.arctan2(2 * np.abs(Ex) * np.abs(Ey) * np.cos(delta), denom) # 计算椭圆率角 chi chi 0.5 * np.arcsin(np.clip(2 * np.abs(Ex) * np.abs(Ey) * np.sin(delta) / (np.abs(Ex)**2 np.abs(Ey)**2), -1, 1))这里arctan2函数是为了避免在分母接近零时产生角度跳变。np.clip则是防止数值溢出导致反正弦函数输入超出定义域。计算完成后psi对应的就是偏振长轴相对于实验室坐标系的旋转角。围绕 BIC 在 k 空间绕一圈psi会旋转整数倍的 π而这个整数就是我们常说的拓扑荷。如果仅仅看某个 k 点的偏振椭圆计算速度很快但若要观察整个 k 空间网格上的偏振分布必须要保证远场投影方向的一致性。一般是设置观测方向固定在 z 轴上也就是theta 0, phi 0这样所有 k 点的远场才具有可比性。如果选取了不同的theta、phi偏振椭圆的参考系会跟着旋转数据分析就会变成一场灾难。我在项目处理中习惯写一个脚本把所有 k 点的远场分量的实部和虚部存入一个三维数组再一次性矩阵化算完所有偏振角效率远高于逐点循环。4. Matlab 脚本实战能带绘图与 Q 值筛选4.1 解析 COMSOL 导出文件从 COMSOL 导出的文本文件表头会包含模型信息和解算器信息需要在数据处理中跳过。常见格式是逗号分隔或制表符分隔。使用readmatrix时要注意指定NumHeaderLines否则第一行被误认成数据。以下脚本给出实际操作中读取数据并生成 Q 值的完整逻辑% 跳过表头读取数值部分 데이터 data readmatrix(band_data.csv, NumHeaderLines, 9); % 列顺序一般按照导出时的设置 kx data(:,1); ky data(:,2); freq_real data(:,3); % 单位可能为 Hz freq_imag abs(data(:,4)); % 转换单位到 THz便于绘图 freq_real freq_real * 1e-12; freq_imag freq_imag * 1e-12; % 计算品质因子 Qfactor freq_real ./ (2 * freq_imag); % 针对 BIC 点虚部消失导致的超大 Q 值设置绘图截断上限 Qfactor_plot min(Qfactor, 1e6);脚本执行逻辑是先把行数据切割成独立列然后对频率进行单位换算。Q 值计算时分母中2 * freq_imag的系数来自品质因子的标准定义。处理完成后Qfactor_plot用于绘图而Qfactor则保留原始数值供后续寻找 BIC 精确位置或分析虚部收敛趋势时使用。截断上限设置为1e6是因为在这个量级继续分开显示已经不会影响能带图中尖峰形状的辨识度。4.2 生成能带图和 Q 值趋势线能带图需要把倒空间路径映射到一维距离轴。在高对称路径Γ-X-M-Γ上距离按倒空间坐标差值累加。例如Γ点为[0,0]X点为[pi/a, 0]那么这段距离就是pi/a之后再累加下一段。重新整理数据之后就能绘制出横轴均匀的能带图。% 计算 k 空间路径距离 k_dist zeros(size(kx)); for i 2:length(kx) if i 1 dk sqrt((kx(i)-kx(i-1))^2 (ky(i)-ky(i-1))^2); k_dist(i) k_dist(i-1) dk; end end % 按距离排序避免连线错乱 [k_dist_sorted, idx] sort(k_dist); freq_sorted freq_real(idx); Q_sorted Qfactor_plot(idx); figure(Position, [100, 100, 800, 450]); yyaxis left; plot(k_dist_sorted, freq_sorted, b-o, LineWidth, 1.2, MarkerSize, 5); ylabel(Frequency (THz)); ylim([min(freq_sorted)*0.98, max(freq_sorted)*1.02]); yyaxis right; semilogy(k_dist_sorted, Q_sorted, r-s, LineWidth, 1.2, MarkerSize, 5); ylabel(Q factor (log scale)); xlabel(k path (Γ-X-M-Γ));这段代码使用双纵轴是因为频率和 Q 值的量纲和数量级差异过大。左轴展示色散关系右轴以对数坐标展示 Q 值。在绘制 BIC 时右轴会在对应位置出现一个尖锐的单点凸起这个位置旁边的频率值就是 BIC 的束缚频率。通过对 Q 值排序后取前几个峰值的索引可以反向定位到对应的 k 点坐标进而明确 BIC 位于哪条高对称线上。4.3 远场偏振的矢量场图绘製远场偏振图展示的是 k 空间网格上的偏振方向分布可以直接用quiver绘制每个箭头代表该 k 点的偏振长轴方向。对 BIC 而言远场偏振在奇点周围形成涡旋分布这种涡旋结构直接证明 BIC 的拓扑性质。% 假设已经从多个 k 点提取远场分量 Ex, Ey % 这里使用二维网格表示 k 空间网格 kx_grid linspace(-1, 1, 31); ky_grid linspace(-1, 1, 31); [kx2, ky2] meshgrid(kx_grid, ky_grid); % Ex_data 和 Ey_data 是 31x31 复数矩阵 Ex Ex_data; % 复数远场 x 分量 Ey Ey_data; % 复数远场 y 分量 % 计算偏振长轴角度 delta angle(Ey) - angle(Ex); psi 0.5 * atan2(2 * abs(Ex) .* abs(Ey) .* cos(delta), abs(Ex).^2 - abs(Ey).^2); % 绘制矢量场图 figure; quiver(kx2, ky2, cos(psi), sin(psi), AutoScale, off); axis square; xlabel(kx (2π/a)); ylabel(ky (2π/a));这里cos(psi)和sin(psi)构成偏振矢量的两个分量。由于atan2返回的角度范围在-π/2到π/2之间绘制时箭头方向不会出现突然反转。若观察到箭头在某个点周围绕圈就说明这里存在偏振奇点可进一步使用后续章节中的相位积分验证拓扑荷。5. 排错与进阶验证数值噪声、网格收敛与拓扑荷5.1 如何判断 Q 值虚部收敛特征频率虚部收敛是光子晶体仿真中最容易翻车的地方。理论上 BIC 的虚部为零但实际网格离散化会引入额外辐射路径相当于给模式加了一个假的损耗通道。网格越粗这个假损耗越大。初学者最典型的错误是用一套标准网格算出一个Q 1e5的结果然后认为 BIC 就是有限 Q 值其实完全是网格限制。要判断收敛性做法是对同一结构依次剖分 3 至 4 组网格从粗糙网格到极细网格每个方向加密一倍。记录每个网格下的 Q 值绘制log(Q)对网格尺寸的曲线。真正的 BIC 点Q 值会随着网格加密持续上升并且斜率趋于稳定也就是说每加密一倍Q 值大致提升一个固定的比例。如果 Q 值在某个网格尺寸后开始波动不变说明已经收敛到了数值极限可以接受。如果不加密就无法继续上升那说明仿真结果没有稳定在 BIC 的物理行为上需要检查边界条件或模式识别是否正确。5.2 偏振涡旋拓扑荷的判定偏振涡旋的拓扑荷是描述 BIC 性质的重要指标它定义为在 k 空间绕行偏振奇点一圈时偏振长轴方向角psi的总变化量除以2π。实际计算时可以选取以奇点为中心的圆周路径依次提取路径上各点的远场分量然后计算展开后的相位% 在闭合路径上计算拓扑荷 % 假设已经沿圆周取了 100 个采样点对应的 Ex, Ey theta_path atan2(ky_points, kx_points); % 极角 psi_wrapped 0.5 * atan2(2 * abs(Ex_p) .* abs(Ey_p) .* cos(delta_p), abs(Ex_p).^2 - abs(Ey_p).^2); % 由于 psi 在 (-pi/2, pi/2) 之间缠绕这里用 unwrap 恢复连续变化 psi_unwrapped unwrap(psi_wrapped * 2) / 2; % 展开时乘以 2 再除以 2 % 拓扑荷 charge round((psi_unwrapped(end) - psi_unwrapped(1)) / (2 * pi));这里隐藏细节在于unwrap的操作对象。由于psi的真实物理变化范围是-π/2到π/2而绕 BIC 一圈时psi应连续变化 π 的整数倍所以先把psi乘以 2使变化范围扩展到-π到π再执行相位展开最后把结果除以 2恢复实际角度。取round是为了消除数值噪声带来的微小偏差。如果计算出的charge为1、-1或更大整数可以判定该奇点对应的拓扑荷。结构对称性会直接决定拓扑荷的数值。四方晶格 Γ 点上的对称性保护型 BIC 通常对应charge 1的涡旋。如果网格剖分时没有保持结构完整对称可能得到非整数的charge此时应先检查网格剖分而不是怀疑物理模型。实际处理中我常用两个不同旋转角度的模型互为验证只有当两次计算均给出相同整数结果时才认为拓扑荷计算可靠。5.3 结合理论文献的边界设定与参数复核资料包中的文献和文档虽然看起来只是背景阅读材料但它们其实是排除模型错误的重要依据。和文献结果对照不顺时核心检查点有三个衬底折射率是否设成了无限衬底、边界条件是否错用成理想导体边界、以及远场投影的积分面是否放置在了过于靠近结构的位置。先看衬底很多论文的计算模型是悬浮膜结构也就是上下都是空气。如果你用有衬底结构模式的泄漏通道和远场辐射方向都会发生重大变化能带图上表现为 Q 值整体下调、BIC 频率移动。再看边界周期性结构要求在 k 空间区分布洛赫周期边界但有些初学者把上下边界也设置成了周期边界就相当于造了一个纵向也无限周期的结构和物理事实完全不符。正确的做法是上下方向使用完美匹配层加散射边界水平方向使用 Floquet 周期边界。最后看远场积分面积分半径不宜过小至少应在结构上方一个波长以上同时保持积分面与完美匹配层内边界平行。只要这三个位置设置得物理自洽COMSOL 算出的远场偏振结果与文献里的远场偏振图就能在拓扑结构上高度吻合。每当我拿到一个新模型总是先按这三个维度检查一遍通常五分钟内就能定位绝大多数仿真数据对不上的原因。本文还有配套的精品资源点击获取
