二维光子晶体能带图计算：平面波展开法与MATLAB实现

简介面向光学工程、物理及相关电子信息专业的高年级本科生与研究生这份Matlab源码包聚焦二维光子晶体的能带结构与电磁场分布计算借助平面波展开法PWE呈现正方晶格、六角晶格与DFB分布反馈结构下的TM/TE模式仿真结果帮助读者直观理解光子禁带、色散关系与模式场分布等核心概念适用于课程设计、毕业设计或科研入门。压缩包共20个文件核心为5个.m脚本分别对应正方形、六角形和DFB结构的主程序与求解函数14张PNG图按晶格类型分类存放展示各结构的能带曲线与电场分布快照另有1份README.md说明文件便于对照运行。整个资源仅2.19MB轻量紧凑便于快速下载与复现。目前已有270人学习可作为熟悉Matlab光学仿真的实用参考。通过运行示例、查看注释和比对结果图读者不仅能获得可直接输出的能带图和场图还能掌握从结构参数定义、PWE求解到后处理绘图的完整流程并迁移至其他光子晶体器件设计中。1. 二维光子晶体能带图到底在算什么一张图背后的本征值问题第一次做二维光子晶体能带图的人很容易把它当成一个时域仿真问题。实际上稳态线性光学下它是个标准的矩阵本征值问题结构在 xy 平面周期排布z 方向均匀电磁场按布洛赫定理展开成平面波叠加把介电常数展开到倒格矢空间后每一个 k 点都变成一个代数方程等着解码。用 MATLAB 写一个平面波展开法PWE主程序几十行就能同时得到能带图和带边模式的场分布不需要 FDTD也不需要商业电磁软件。能带图解决的是“哪些频率能在这个周期结构里存在”场图解决的是“这些频率的电磁能到底局域在哪个位置、长成什么样”。两者合在一起才能回答带隙、缺陷模、波导耦合、微腔设计这些实际问题。这篇内容按“理论推导 → MATLAB 最小实现 → 场分布回代 → 可靠性检查”的顺序展开适合正在做光学课程项目、光子晶体器件仿真或者刚接触平面波展开法想少走弯路的人。2. 平面波展开法二维光子晶体能带图从 Maxwell 方程到矩阵方程的推导2.1 为什么二维问题可以降成 TE/TM 两个标量方程二维光子晶体指介电常数只在一个平面内周期变化另一个方向完全均匀。设周期方向为 xy 平面z 方向不变那么 Maxwell 方程组在无源、无磁响应介质中会解耦成两套独立的偏振电场沿 z 方向的 TM 偏振以及磁场沿 z 方向的 TE 偏振。这个解耦是二维问题比三维问题简单一个量级的根本原因。TM 偏振下只有 E_z 这一个电场分量控制方程从矢量方程坍缩成标量亥姆霍兹方程−∇²E_z (ω/c)² ε(x, y) E_zTE 偏振对应 H_z方程形式稍复杂因为介电常数和微分算子不能交换顺序梯度项会耦合进来。大多数能带图教学实现都从 TM 入手矩阵更小物理图像也更直观。需要说明的是二维光子晶体的带隙对偏振敏感同一个结构 TM 有带隙而 TE 可能完全没有所以计算前必须明确是哪套偏振。2.2 PWE 的矩阵化TM 与 TE 的差异把 E_z 按布洛赫定理展开成平面波叠加E_z(r) Σ_G A_G e^(i(kG)·r)其中 G 是倒格矢k 是第一布里渊区内的波矢。将展开代入标量方程再乘 e^(−i(kG)·r) 在晶胞内积分利用平面波的正交性就得到所谓 PWE 本征方程Σ_{G} κ(G−G) |kG|² A_G (ω/c)² A_G这里 κ(G) 是 1/ε(r) 的傅里叶系数不是 ε(r) 的傅里叶系数。这是最常见的实现误区有人直接展开 ε(r)代入后方程形式就错了得到的带隙位置和宽度会有明显偏差。原因在于原始方程中 ε 乘在 E_z 上移项后变成 1/ε 作用在 ∇²E_z 上周期函数展开的是倒数介电常数。TM 与 TE 的差异可以整理成一张对照表写代码时直接按表里的矩阵元素实现偏振标量方程PWE 矩阵元素展开系数TM−∇²E_z (ω/c)² ε E_zκ(G−G)·kGTE−∇·(1/ε ∇H_z) (ω/c)² H_z含 (kG)·(kG′) 的完整卷积同样用 FT(1/ε)但梯度算符引入额外耦合TE 的矩阵比 TM 更密对角占优性也更差同样 NG 下求本征值更慢。所以如果只是想快速验证一个结构的带隙先跑 TM 是性价比最高的选择。2.3 扫描布里渊区边界能带图的横坐标该怎么取能带图不是在整个二维布里渊区上画曲面而是沿着不可约布里渊区的高对称边界扫一条折线路径。原因是带隙的上下边界通常出现在高对称点或高对称连线上沿边界扫掠就能抓住带隙的主要特征计算量也小得多。不同晶格的路径和坐标约定不同常见的两组需要记牢晶格类型扫描路径高对称点坐标单位 2π/a正方格子Γ → X → M → ΓΓ(0,0)X(0.5,0)M(0.5,0.5)三角/六角晶格Γ → K → M → ΓΓ(0,0)K(1/3,1/3)M(0,0.5)横坐标本身不是频率也不是 k 的模而是沿路径的累积长度。把每一段的 k 点间距累加起来横轴刻度标在高对称点处图的可读性会好很多。这个路径坐标在第三章的代码里直接体现。3. MATLAB 计算二维光子晶体能带图最小实现与参数设置3.1 程序结构从周期结构参数到傅里叶系数先确定物理参数晶格常数 a、介质柱半径 ra、柱体介电常数 eps1、背景介电常数 eps2。这里默认正方格子、介质柱埋在背景介质中TM 偏振。圆截面结构的傅里叶系数有解析表达式不需要在实空间画网格做 FFT这个细节决定了程序的速度和精度。倒格矢按 (2NG1)×(2NG1) 截断总平面波数为 N(2NG1)²。NG6 时 N169矩阵是 169×169MATLAB 里一次 eig 求解在零点几秒量级NG10 时 N441仍然可以接受。解析傅里叶系数的公式为G0 时等于填充率加权平均倒数介电常数G≠0 时用 2f·Δ(1/ε)·J1(x)/x其中 x|G|·R。3.2 可运行的 MATLAB 主程序TM 偏振下面是最小可运行的 TM 偏振能带图函数。为了可读性Kappa 矩阵的组装用了显式查找NG≤8 时速度完全够用。function [kv, freq] pwe_2d_tm(a, ra, eps1, eps2, NG, kpath, nbands) % 正方格子二维光子晶体 TM 偏振能带结构平面波展开法 % 输入: % a 晶格常数长度与 ra 保持同一单位即可 % ra 介质柱半径 % eps1 介质柱相对介电常数 % eps2 背景相对介电常数 % NG 每个方向的倒格矢截断数总平面波数 (2*NG1)^2 % kpath 高对称点路径坐标按 2*pi/a 归一化 % nbands 需要输出的能带数 % 输出: % kv 沿路径累积长度用作横坐标 % freq 归一化频率 a/lambda多行 nbands 列 g (-NG:NG).; [G1, G2] meshgrid(g, g); G1 G1(:); G2 G2(:); N length(G1); % 1/eps 的傅里叶展开系数圆截面解析式 f pi * ra^2 / a^2; x 2*pi*ra/a * sqrt(G1.^2 G2.^2); inv_eps zeros(N, 1); for i 1:N if x(i) 1e-12 inv_eps(i) f/eps1 (1-f)/eps2; else inv_eps(i) 2*f*(1/eps1 - 1/eps2) * besselj(1, x(i)) / x(i); end end % 预组装 Kappa 矩阵元素为 kappa(G_i - G_j) Kmat zeros(N, N); for i 1:N for j 1:N dG [G1(i)-G1(j), G2(i)-G2(j)]; idx find(G1dG(1) G2dG(2), 1); Kmat(i, j) inv_eps(idx); end end % 沿高对称路径生成 k 点 nseg size(kpath, 1) - 1; npts 200; pts []; for s 1:nseg pts [pts; linspace(kpath(s,:), kpath(s1,:), npts)]; end nk size(pts, 1); % 每个 k 点组装矩阵并求解 freq zeros(nk, nbands); for ik 1:nk kg sqrt((pts(ik,1)G1).^2 (pts(ik,2)G2).^2); M (kg * kg.) .* Kmat; % 对称化特征值与原方程一致 E eig(M); E sort(real(E(E 1e-10))); freq(ik,:) sqrt(E(1:nbands)); end % 横坐标路径累积长度 dk sqrt(sum(diff(pts).^2, 2)); kv [0; cumsum(dk)].; end代码里有一个容易看漏的细节M 矩阵用(kg * kg.) .* Kmat构造。原始方程中 |kG′|² 只乘在列索引上矩阵不对称这里把 |kG| 和 |kG′| 各分一半乘到 Kappa 两侧构成相似变换特征值不变但矩阵变成对称矩阵数值稳定性更好。E(E 1e-10)是为了滤掉 Γ 点处 kG0 引入的零特征值这些零模不是物理模式不滤掉会占用能带序号。调用脚本如下a 1; ra 0.2*a; eps1 12; % 硅 eps2 1; % 空气 NG 6; kpath [0 0; 0.5 0; 0.5 0.5; 0 0]; % Gamma-X-M-Gamma [kv, freq] pwe_2d_tm(a, ra, eps1, eps2, NG, kpath, 8); plot(kv, freq, LineWidth, 1.2); axis([0 kv(end) 0 1]); xlabel(波矢路径); ylabel(归一化频率 a/\lambda); % 手动标出高对称点npts200 时每段分界索引为 200 和 400 set(gca, XTick, [0 kv(200) kv(400) kv(600)]); set(gca, XTickLabel, {\Gamma, X, M, \Gamma});输出的 freq 单位是 a/λ这是光子晶体文献里最常见的归一化方式。画图时纵轴取 0 到 1 就够用更高的带通常不是关注对象。跑通这个脚本后能带图里应当能看到低频段近似直线、在高对称点出现能带折叠如果介质柱与背景折射率对比足够大比如 12:1会在某个频率区间看到明显的空白带隙。3.3 倒格矢截断数 NG矩阵大小、耗时与精度的取舍NG 是整个计算里最重要的收敛参数。矩阵维度随 NG 平方增长但能带频率的收敛速度约为一阶NG 从 4 加到 8带边频率的变化通常从几个百分点降到零点几个百分点。不同场景下的推荐取值如下NG平面波总数单 k 点 eig 参考耗时适用场景3490.02 秒快速验证、教学演示51210.1 秒初步扫描结构参数61690.3 秒常规能带图、论文插图82891.5 秒带边频率精算、场分布回代这里的耗时是个人电脑上的相对量级只用于选型参考。实际项目中我一般先用 NG5 扫一遍参数空间找趋势确定感兴趣的结构后再用 NG8 精算最终能带图。不要在参数扫描阶段用大 NG否则一次扫几十个半径值会等很久。矩阵组装还有一个可以立刻优化的点Kmat 只依赖倒格矢差值与 k 无关所以放在 k 循环外只算一次(kg * kg.) .* Kmat本身就是向量化操作NG8 时单 k 点也很快。如果还想继续提速把 Kmat 改成稀疏存储再用 eigs 求最低若干条带效率能再上一个台阶这个放在第 5 章展开。4. 场分布本征矢回代、实空间成像与超胞法4.1 能带图上取一个点怎么把它还原成实空间场能带图给出的是色散关系同一个频率可能对应多个模式不画场图就无法知道电磁能局域在哪里。PWE 的优势在于本征值对应的本征矢本身就是平面波展开系数 A_G频率解出来之后场图几乎是免费的。具体做法是在目标 k 点和目标能带序号处取本征矢按 E_z(r) Σ_G A_G e^(i(kG)·r) 叠加到实空间网格上。这里有一个物理细节要分清如果只叠加 G 的周期项得到的是布洛赫函数的周期部分 u_k(r)它反映晶胞内部的场调制如果把 kG 一起放进相位得到的是完整波函数 E_z(r)能看到波长远小于晶胞时的快速振荡。画带边模式时两者都能用但解释方式不同。4.2 对称化本征矢回代一个容易错的换算第三章用对称化矩阵求解特征向量不能直接当 A_G 用。对称化后的本征矢 x 与原方程振幅 A_G 之间差一个对角变换A_G x / |kG|。如果跳过这一步场图在 kG 接近零的位置会出现错误的幅度放大。% 假设已在某 k 点取得对称化本征矢 V(:, b) % kg 是该 k 点下的 |kG| 向量长度为 N kg_safe kg; kg_safe(kg_safe 1e-8) 1; % 避免除零 A_G V(:, b) ./ kg_safe; % 还原为原方程平面波振幅 % 在单个晶胞内画布洛赫周期部分 u_k nx 96; [XX, YY] meshgrid(linspace(0, 1, nx), linspace(0, 1, nx)); u_k zeros(nx, nx); for ig 1:N u_k u_k A_G(ig) * exp(1i*2*pi*(G1(ig)*XX G2(ig)*YY)); end % 画全波 E_z 时把上面的 G1(ig) 换成 k(1)G1(ig)G2 同理 figure; surf(XX, YY, real(u_k), EdgeColor, none); view(2); axis equal tight;这个换算很多人会在第一次实现时漏掉。原因在于对称化矩阵的特征向量本身满足的是另一种归一化直接用 x 叠加平面波等效于给每个平面波分量乘了 |kG|场分布节点位置可能不明显变但幅度分布会偏向高倒格矢分量。逻辑上还原到原方程之后A_G 的量纲和物理意义才一致。实空间网格数 nx 与 NG 要匹配。NG6 时倒格矢最大到 6一个晶胞内最高空间频率对应 12 个振荡周期奈奎斯特条件要求 nx 至少 24实用中取 96 或 128 足够平滑。网格过大不会增加物理信息只会拖慢 surf 渲染。4.3 超胞近似缺陷模与带隙内平带完整光子晶体器件通常在周期结构里引入缺陷比如拿掉一根柱子或改变某根柱子的半径。缺陷破坏了平移对称性严格来说不能直接用原晶格的 k 点扫描工程上常用超胞近似把若干个原胞拼成一个超胞让缺陷位于超胞中心然后在超胞的 Γ 点计算能带。超胞方法的代价是倒格子缩小。3×3 超胞第一布里渊区缩小到原来的 1/3倒格矢步长变小原来在 k 路径上的模式全部折叠到 Γ 点附近为了保持傅里叶级数收敛NG 通常要同步增加矩阵维数上升很快。实际中 3×3 或 5×5 超胞配上 NG46 是比较常见的折中。缺陷模式在超胞能带图里表现为带隙中的一条平带这条带的频率对缺陷半径和位置很敏感是微腔设计的主要依据。超胞大小倒格矢范围变化推荐 NG矩阵维数参考1×1原倒格矢61693×3倒格矢缩小 3 倍51215×5倒格矢缩小 5 倍481超胞法只适合单个局域缺陷缺陷之间的间距小于超胞尺寸时相邻镜像耦合会污染结果。当需要计算波导或慢光模式时超胞尺寸还要进一步加大但这时我会直接换用 FDTD 或有限元PWE 的矩阵规模容易失控。判断标准很简单看缺陷场分布是否在超胞边界处衰减到可忽略水平。5. 能带图可靠性的三个把关点收敛性、归一化与 Γ 点零模5.1 NG 截断收敛性用表格而不是感觉判断能带图画出来之后第一个要回答的问题是“这个图可信吗”。最直观的做法是固定所有物理参数只增加 NG观察带隙边界的变化。下面是 r/a0.2、ε12:1、正方格子 TM 偏振下带隙上下边频的典型收敛趋势数值本身是示意性的关键是看变化率随 NG 的下降NG归一化带隙 a/λ与 NG9 的偏差30.0209.5%50.02182.1%70.02220.5%90.0223基准判断标准可以这样定当 NG 增大两档后带隙宽度变化小于 1%认为收敛。论文级别的图我会跑到 NG9 或 10日常工程用 NG6 足够。收敛性检查还有一个附带作用如果带隙随 NG 剧烈震荡说明结构里有尖锐的介电界面这时优先检查 1/ε 的傅里叶系数是否正确而不是盲目加大 NG。5.2 单位约定与 Γ 点零频率模与文献对比时最常踩的坑是单位不一致。PWE 本征值解出来的是 (ωa/2πc)²开方就是 a/λ这也是我输出的单位。另一些文献把纵轴标成 ωa/2πc两者数值完全相同但如果对方用的是归一化到ωa/c数值会差 2π 倍。对比前先看清图纵轴标签比对着曲线形状硬猜可靠得多。Γ 点零频率模也是必踩项。k0 且 G0 时矩阵对应行全部为零本征值 0 会被算出来。第三章代码里用E 1e-10滤掉了它但如果改用 eigs 求部分本征值这个方法要保留。正确做法是求 nbands1 条带去掉零模后取前 nbands 条。否则能带图第一条带会是一条恒为零的直线带隙判断直接出错。5.3 eigs 与小批量扫掠技巧大 NG 下全矩阵 eig 求解 N 个本征值浪费在不需要的高频带上。改用 eigs 只求最低若干条带是提速最明显的一步opts.tol 1e-8; opts.issym true; [V, D] eigs(M, nbands1, smallestabs, opts); E diag(D); E sort(real(E(E 1e-10))); % 去掉 Γ 点零模smallestabs在 MATLAB 对称特征值问题里按绝对值找最小的若干本征值issym必须设成 true因为 M 是对称化之后的矩阵。需要回代场图时eigs 返回的 V 直接接上第四章的 A_G 换算整条流程闭环。最后一招是保留参数化调用。把半径、NG、k 路径都写成函数参数用一个外层循环扫半径批量输出带隙随 r/a 的变化曲线。这样能带图不再是单张图而是一条设计曲线光子晶体器件优化的第一个步骤就落地了。场图完成后检查 u_k 在晶胞边界上的连续性正方晶格下旋转 90 度后场分布应与原图一致归一化误差小于 1e-3 基本可以确认整套傅里叶系数和回代换算没有手滑。本文还有配套的精品资源点击获取