4、二维光子晶体能带计算：平面波展开法（PWE）实现

二维光子晶体的能带计算，说白了就是求解麦克斯韦方程组在周期性介质中的本征值问题。我个人习惯用平面波展开法（PWE）来做这件事，因为它物理图像清晰，代码实现也不复杂。嗯，咱们一步步来拆解。

4.1 平面波展开法的核心思想

为什么要用平面波展开？你想想看，光子晶体是周期性的结构，而平面波本身就是周期函数。把电磁场用平面波叠加起来，正好匹配这种周期性。我在项目中遇到过不少新手，一上来就想着用有限元法，其实对于能带计算，PWE 是最直接的选择。

核心思路就三步：

• 把介电常数展开成傅里叶级数
• 把电磁场也用平面波展开
• 代入麦克斯韦方程，得到本征值方程

说白了，就是把微分方程变成了代数方程。计算机最擅长解代数方程了。

4.2 二维光子晶体的能带结构

二维光子晶体，结构在 x-y 平面内周期性排列，z 方向是均匀的。这时候光传播可以分成两种模式：

• TM 模式：电场沿 z 方向，磁场在 x-y 平面内
• TE 模式：磁场沿 z 方向，电场在 x-y 平面内

这两种模式的本征方程不一样。我记得刚开始做的时候，总是把 TM 和 TE 搞混，后来养成了一个习惯：先看电场方向，电场沿 z 就是 TM，磁场沿 z 就是 TE。

4.3 TM/TE 模式分析

对于 TM 模式，我们只需要求解电场 E_z 分量。方程简化成：

∑_{G'} κ(G-G') (k+G)·(k+G') E_z(k+G') = (ω²/c²) E_z(k+G)

其中 κ 是介电常数的傅里叶系数，G 是倒格矢，k 是波矢。

对于 TE 模式，需要求解磁场 H_z 分量，方程稍微复杂一点：

∑_{G'} κ(G-G') (k+G)×(k+G')·ẑ H_z(k+G') = (ω²/c²) H_z(k+G)

4.4 Python 代码实现

下面是我常用的 PWE 计算代码。这个版本我优化过几次，跑起来很稳定。

import numpy as np
from scipy import linalg

def pwe_2d(epsilon, a, N, k_points):
    """
    二维光子晶体能带计算
    epsilon: 介电常数分布 (Nx x Ny)
    a: 晶格常数
    N: 平面波数量 (N x N)
    k_points: 布里渊区路径上的k点
    """
    # 构建倒格矢
    Gx = 2*np.pi/a * np.fft.fftfreq(N, d=1/N)
    Gy = 2*np.pi/a * np.fft.fftfreq(N, d=1/N)
    Gx, Gy = np.meshgrid(Gx, Gy)
    
    # 计算介电常数的傅里叶系数
    epsilon_k = np.fft.fft2(epsilon) / N**2
    
    bands = []
    for kx, ky in k_points:
        # 构建哈密顿矩阵
        H = np.zeros((N*N, N*N), dtype=complex)
        for i in range(N):
            for j in range(N):
                for m in range(N):
                    for n in range(N):
                        G_diff_x = Gx[i,j] - Gx[m,n]
                        G_diff_y = Gy[i,j] - Gy[m,n]
                        kG_dot = (kx+Gx[i,j])*(kx+Gx[m,n]) + \
                                 (ky+Gy[i,j])*(ky+Gy[m,n])
                        H[i*N+j, m*N+n] = epsilon_k[i-m, j-n] * kG_dot
        
        # 求解本征值
        eigvals = linalg.eigvalsh(H)
        bands.append(np.sort(np.real(eigvals[:10])))
    
    return np.array(bands)

# 使用示例
a = 1.0  # 晶格常数
N = 11   # 平面波数量
# 构建介电常数分布（空气孔在介质中）
epsilon = np.ones((N, N)) * 12.0  # 介质背景
r = 0.3 * a  # 空气孔半径
for i in range(N):
    for j in range(N):
        x = (i - N//2) * a/N
        y = (j - N//2) * a/N
        if x**2 + y**2 < r**2:
            epsilon[i,j] = 1.0  # 空气孔

# 定义布里渊区路径
k_path = [(0,0), (np.pi/a,0), (np.pi/a,np.pi/a), (0,0)]
bands = pwe_2d(epsilon, a, N, k_path)

4.5 结果分析与可视化

算完能带之后，怎么判断有没有带隙？我一般看两个指标：

• 相邻能带之间有没有频率间隔
• 这个间隔在整个布里渊区是否都存在

如果 TM 和 TE 模式都有带隙，而且带隙有重叠，那就是完全带隙。这种材料做波导、谐振腔效果特别好。

下面是我用 SVG 画的一张流程图，展示了 PWE 计算的整体流程：

4.6 避坑指南

做 PWE 计算，有几个坑我踩过，分享给你：

• 平面波数量太少：我曾经用 5x5 的平面波算能带，结果带隙位置完全不对。建议至少 11x11 起步。
• 介电常数突变：空气孔边界处介电常数突变，傅里叶级数收敛慢。可以用平滑过渡来改善。
• k 点采样不足：布里渊区边界处的能带容易漏掉。我一般每条边取 20 个点以上。

二维光子晶体的 PWE 计算，说白了就是这几步。代码跑通了，能带图出来了，你就能看到光子带隙在哪里。这个带隙宽度、位置，直接决定了你的器件性能。嗯，动手试试吧。

---

公众号:蓝海资料掘金营，微信deep3321