企业官网建设流程全解析

简介：面向光学工程、物理及相关电子信息专业的高年级本科生与研究生，这份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_z

TE 偏振对应 H_z，方程形式稍复杂，因为介电常数和微分算子不能交换顺序，梯度项会耦合进来。大多数能带图教学实现都从 TM 入手，矩阵更小，物理图像也更直观。需要说明的是，二维光子晶体的带隙对偏振敏感，同一个结构 TM 有带隙而 TE 可能完全没有，所以计算前必须明确是哪套偏振。

2.2 PWE 的矩阵化：TM 与 TE 的差异

把 E_z 按布洛赫定理展开成平面波叠加：

E_z(r) = Σ_G A_G e^(i(k+G)·r)

其中 G 是倒格矢，k 是第一布里渊区内的波矢。将展开代入标量方程，再乘 e^(−i(k+G)·r) 在晶胞内积分，利用平面波的正交性，就得到所谓 PWE 本征方程：

Σ_{G'} κ(G−G') |k+G'|² A_G' = (ω/c)² A_G

这里 κ(G) 是 1/ε(r) 的傅里叶系数，不是 ε(r) 的傅里叶系数。这是最常见的实现误区：有人直接展开 ε(r)，代入后方程形式就错了，得到的带隙位置和宽度会有明显偏差。原因在于原始方程中 ε 乘在 E_z 上，移项后变成 1/ε 作用在 ∇²E_z 上，周期函数展开的是倒数介电常数。

TM 与 TE 的差异可以整理成一张对照表，写代码时直接按表里的矩阵元素实现：

TE 的矩阵比 TM 更密，对角占优性也更差，同样 NG 下求本征值更慢。所以如果只是想快速验证一个结构的带隙，先跑 TM 是性价比最高的选择。

2.3 扫描布里渊区边界：能带图的横坐标该怎么取

能带图不是在整个二维布里渊区上画曲面，而是沿着不可约布里渊区的高对称边界扫一条折线路径。原因是带隙的上下边界通常出现在高对称点或高对称连线上，沿边界扫掠就能抓住带隙的主要特征，计算量也小得多。

不同晶格的路径和坐标约定不同，常见的两组需要记牢：

横坐标本身不是频率也不是 k 的模，而是沿路径的累积长度。把每一段的 k 点间距累加起来，横轴刻度标在高对称点处，图的可读性会好很多。这个路径坐标在第三章的代码里直接体现。

3. MATLAB 计算二维光子晶体能带图：最小实现与参数设置

3.1 程序结构：从周期结构参数到傅里叶系数

先确定物理参数：晶格常数 a、介质柱半径 ra、柱体介电常数 eps1、背景介电常数 eps2。这里默认正方格子、介质柱埋在背景介质中，TM 偏振。圆截面结构的傅里叶系数有解析表达式，不需要在实空间画网格做 FFT，这个细节决定了程序的速度和精度。

倒格矢按 (2NG+1)×(2NG+1) 截断，总平面波数为 N=(2NG+1)²。NG=6 时 N=169，矩阵是 169×169，MATLAB 里一次 eig 求解在零点几秒量级；NG=10 时 N=441，仍然可以接受。解析傅里叶系数的公式为：G=0 时等于填充率加权平均倒数介电常数，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*NG+1)^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(G1==dG(1) & G2==dG(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(s+1,:), 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构造。原始方程中 |k+G′|² 只乘在列索引上，矩阵不对称；这里把 |k+G| 和 |k+G′| 各分一半乘到 Kappa 两侧，构成相似变换，特征值不变，但矩阵变成对称矩阵，数值稳定性更好。E(E > 1e-10)是为了滤掉 Γ 点处 k+G=0 引入的零特征值，这些零模不是物理模式，不滤掉会占用能带序号。

调用脚本如下：

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'); % 手动标出高对称点，npts=200 时每段分界索引为 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=5 扫一遍参数空间找趋势，确定感兴趣的结构后，再用 NG=8 精算最终能带图。不要在参数扫描阶段用大 NG，否则一次扫几十个半径值会等很久。

矩阵组装还有一个可以立刻优化的点：Kmat 只依赖倒格矢差值，与 k 无关，所以放在 k 循环外只算一次；(kg * kg.') .* Kmat本身就是向量化操作，NG=8 时单 k 点也很快。如果还想继续提速，把 Kmat 改成稀疏存储，再用 eigs 求最低若干条带，效率能再上一个台阶，这个放在第 5 章展开。

4. 场分布：本征矢回代、实空间成像与超胞法

4.1 能带图上取一个点，怎么把它还原成实空间场

能带图给出的是色散关系，同一个频率可能对应多个模式，不画场图就无法知道电磁能局域在哪里。PWE 的优势在于，本征值对应的本征矢本身就是平面波展开系数 A_G，频率解出来之后，场图几乎是免费的。

具体做法是：在目标 k 点和目标能带序号处取本征矢，按 E_z(r) = Σ_G A_G e^(i(k+G)·r) 叠加到实空间网格上。这里有一个物理细节要分清：如果只叠加 G 的周期项，得到的是布洛赫函数的周期部分 u_k(r)，它反映晶胞内部的场调制；如果把 k+G 一起放进相位，得到的是完整波函数 E_z(r)，能看到波长远小于晶胞时的快速振荡。画带边模式时两者都能用，但解释方式不同。

4.2 对称化本征矢回代：一个容易错的换算

第三章用对称化矩阵求解，特征向量不能直接当 A_G 用。对称化后的本征矢 x 与原方程振幅 A_G 之间差一个对角变换：A_G = x / |k+G|。如果跳过这一步，场图在 k+G 接近零的位置会出现错误的幅度放大。

% 假设已在某 k 点取得对称化本征矢 V(:, b) % kg 是该 k 点下的 |k+G| 向量，长度为 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 叠加平面波，等效于给每个平面波分量乘了 |k+G|，场分布节点位置可能不明显变，但幅度分布会偏向高倒格矢分量。逻辑上，还原到原方程之后，A_G 的量纲和物理意义才一致。

实空间网格数 nx 与 NG 要匹配。NG=6 时倒格矢最大到 6，一个晶胞内最高空间频率对应 12 个振荡周期，奈奎斯特条件要求 nx 至少 24，实用中取 96 或 128 足够平滑。网格过大不会增加物理信息，只会拖慢 surf 渲染。

4.3 超胞近似：缺陷模与带隙内平带

完整光子晶体器件通常在周期结构里引入缺陷，比如拿掉一根柱子或改变某根柱子的半径。缺陷破坏了平移对称性，严格来说不能直接用原晶格的 k 点扫描；工程上常用超胞近似，把若干个原胞拼成一个超胞，让缺陷位于超胞中心，然后在超胞的 Γ 点计算能带。

超胞方法的代价是倒格子缩小。3×3 超胞第一布里渊区缩小到原来的 1/3，倒格矢步长变小，原来在 k 路径上的模式全部折叠到 Γ 点附近；为了保持傅里叶级数收敛，NG 通常要同步增加，矩阵维数上升很快。实际中 3×3 或 5×5 超胞配上 NG=4～6 是比较常见的折中。缺陷模式在超胞能带图里表现为带隙中的一条平带，这条带的频率对缺陷半径和位置很敏感，是微腔设计的主要依据。

超胞法只适合单个局域缺陷；缺陷之间的间距小于超胞尺寸时，相邻镜像耦合会污染结果。当需要计算波导或慢光模式时，超胞尺寸还要进一步加大，但这时我会直接换用 FDTD 或有限元，PWE 的矩阵规模容易失控。判断标准很简单：看缺陷场分布是否在超胞边界处衰减到可忽略水平。

5. 能带图可靠性的三个把关点：收敛性、归一化与 Γ 点零模

5.1 NG 截断收敛性：用表格而不是感觉判断

能带图画出来之后，第一个要回答的问题是“这个图可信吗”。最直观的做法是固定所有物理参数，只增加 NG，观察带隙边界的变化。下面是 r/a=0.2、ε=12:1、正方格子 TM 偏振下带隙上下边频的典型收敛趋势，数值本身是示意性的，关键是看变化率随 NG 的下降：

判断标准可以这样定：当 NG 增大两档后，带隙宽度变化小于 1%，认为收敛。论文级别的图我会跑到 NG=9 或 10，日常工程用 NG=6 足够。收敛性检查还有一个附带作用：如果带隙随 NG 剧烈震荡，说明结构里有尖锐的介电界面，这时优先检查 1/ε 的傅里叶系数是否正确，而不是盲目加大 NG。

5.2 单位约定与 Γ 点零频率模

与文献对比时最常踩的坑是单位不一致。PWE 本征值解出来的是 (ωa/2πc)²，开方就是 a/λ，这也是我输出的单位。另一些文献把纵轴标成 ωa/2πc，两者数值完全相同；但如果对方用的是归一化到ωa/c，数值会差 2π 倍。对比前先看清图纵轴标签，比对着曲线形状硬猜可靠得多。

Γ 点零频率模也是必踩项。k=0 且 G=0 时，矩阵对应行全部为零，本征值 0 会被算出来。第三章代码里用E > 1e-10滤掉了它，但如果改用 eigs 求部分本征值，这个方法要保留。正确做法是求 nbands+1 条带，去掉零模后取前 nbands 条。否则能带图第一条带会是一条恒为零的直线，带隙判断直接出错。

5.3 eigs 与小批量扫掠技巧

大 NG 下全矩阵 eig 求解 N 个本征值，浪费在不需要的高频带上。改用 eigs 只求最低若干条带，是提速最明显的一步：

opts.tol = 1e-8; opts.issym = true; [V, D] = eigs(M, nbands+1, '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 基本可以确认整套傅里叶系数和回代换算没有手滑。

本文还有配套的精品资源，点击获取