通过场分布得到光子晶体的色散光子晶体那彩虹般的色散特性总让人着迷但真正上手计算时总有种「知道原理却不知怎么操作」的尴尬。今天咱们来点硬核实操直接通过电磁场分布数据倒推色散关系——这个思路在缺陷态分析里尤其好用。先看核心逻辑当特定频率的光波在光子晶体中传播时其场分布会呈现周期性特征。假设我们通过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准则。比如晶格常数a0.5μm时dx至少得小于a/4否则高频模式会混叠。接着做时空傅里叶变换# 沿时间轴做FFT Ez_fft np.fft.fft(Ez, axis0) 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, height0.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决定了频率分辨率Δf1/T。如果仿真只跑了1ps频率分辨率顶多到1THz这对可见光波段显然不够用。这种方法的妙处在于能处理非理想结构比如存在缺陷或无序的情况。传统平面波展开法算不出的复杂结构用场分布反推反而更直接。不过要注意边界反射的影响最好在仿真时加PML层吸收边界。