通过场分布得到光子晶体的色散
光子晶体那彩虹般的色散特性总让人着迷，但真正上手计算时总有种「知道原理却不知怎么操作」的尴尬。今天咱们来点硬核实操，直接通过电磁场分布数据倒推色散关系——这个思路在缺陷态分析里尤其好用。
先看核心逻辑：当特定频率的光波在光子晶体中传播时，其场分布会呈现周期性特征。假设我们通过FDTD仿真得到某时刻的电场分布（比如Ez分量），可以用二维傅里叶变换把空间分布转成波矢k的信息。不同频率对应的模式会在(k,ω)空间形成亮线，这就是色散曲线。
上代码！假设我们有个方格子光子晶体，先读取场数据：
import numpy as np import matplotlib.pyplot as plt with h5py.File('Ez_snapshot.h5', 'r') as hf: Ez = hf['Ez'][:] # 形状为(nt, nx, ny) x = hf['x'][:] y = hf['y'][:] t = hf['t'][:] dt = t[1] - t[0] # 时间步长 dx = x[1] - x[0] # 空间步长这里有个坑：空间采样间隔必须满足Nyquist准则。比如晶格常数a=0.5μm时，dx至少得小于a/4，否则高频模式会混叠。接着做时空傅里叶变换：
# 沿时间轴做FFT Ez_fft = np.fft.fft(Ez, axis=0) freqs = np.fft.fftfreq(Ez.shape[0], dt) pos_freq = freqs > 0 # 取正频率部分 # 空间二维FFT kx = np.fft.fftshift(np.fft.fftfreq(Ez.shape[1], dx)) ky = np.fft.fftshift(np.fft.fftfreq(Ez.shape[2], dx)) Ez_k = np.fft.fftshift(np.fft.fft2(Ez_fft, axes=(1,2)), axes=(1,2))注意fftshift的用法——把低频分量移到中间。这时候Ez_k的形状是（频率点数，kx点数，ky点数）。对于正方晶格，只需关注Γ-X方向（kx从0到π/a），所以切片处理：
# 提取Γ-X路径 k_line = kx[kx >= 0] Ez_line = Ez_k[:, len(kx)//2:, len(ky)//2]现在需要找出每个k对应的主导频率。一个骚操作是用局部极大值检测：
from scipy.signal import find_peaks k_points = [] omega_points = [] for ik in range(Ez_line.shape[1]): spectrum = np.abs(Ez_line[:,ik]) peaks, _ = find_peaks(spectrum, height=0.2*np.max(spectrum)) if len(peaks) > 0: main_peak = peaks[np.argmax(spectrum[peaks])] omega_points.append(freqs[main_peak]) k_points.append(k_line[ik]) # 转换为归一化频率 a = 0.5e-6 # 晶格常数 c = 3e8 omega_norm = np.array(omega_points)*a/(2*np.pi*c)最后画出的色散曲线可能会发现某些频段缺失，这时候得检查仿真时间是否够长——时间窗口T决定了频率分辨率Δf=1/T。如果仿真只跑了1ps，频率分辨率顶多到1THz，这对可见光波段显然不够用。
这种方法的妙处在于能处理非理想结构，比如存在缺陷或无序的情况。传统平面波展开法算不出的复杂结构，用场分布反推反而更直接。不过要注意边界反射的影响，最好在仿真时加PML层吸收边界。
