# 二维六方光子晶体能带计算的Python实现方法 ## 1. 光子晶体能带计算理论基础 光子晶体是一种介电常数呈周期性排列的人工微结构，能够调控电磁波的传播特性。二维六方晶格光子晶体具有典型的周期性结构，其能带计算基于布洛赫定理和麦克斯韦方程组的周期性解[ref_1]。 ### 1.1 关键物理概念 - **布里渊区**：六方晶格的倒格子空间，高对称路径通常为Γ-M-K-Γ - **特征频率**：对应特定波矢的本征模式频率 - **能带折叠**：由于周期性边界条件导致的能带结构特征[ref_2] ## 2. Python实现的核心算法 ### 2.1 平面波展开法 平面波展开法是计算光子晶体能带结构的经典方法，以下是基于Python的完整实现： ```python import numpy as np import matplotlib.pyplot as plt from scipy.linalg import eigh from scipy.sparse.linalg import eigs import itertools class HexagonalPhotonicCrystal: def __init__(self, r=0.3, eps_a=13.0, eps_b=1.0, lattice_constant=1.0): """ 初始化六方晶格光子晶体参数 参数: r: 介质柱半径 (相对于晶格常数) eps_a: 介质柱介电常数 eps_b: 背景介电常数 lattice_constant: 晶格常数 """ self.r = r self.eps_a = eps_a self.eps_b = eps_b self.a = lattice_constant # 六方晶格基矢 self.a1 = np.array([self.a, 0.0]) self.a2 = np.array([-self.a/2, self.a*np.sqrt(3)/2]) # 倒格子基矢 self.b1 = 2*np.pi*np.array([1/self.a, 1/(self.a*np.sqrt(3))]) self.b2 = 2*np.pi*np.array([0, 2/(self.a*np.sqrt(3))]) def get_reciprocal_lattice_vectors(self, G_max): """ 生成倒格矢集合 """ G_vectors = [] n_max = int(G_max * self.a / (2*np.pi)) for n1 in range(-n_max, n_max+1): for n2 in range(-n_max, n_max+1): G = n1 * self.b1 + n2 * self.b2 if np.linalg.norm(G) <= G_max: G_vectors.append(G) return np.array(G_vectors) def fourier_coefficient(self, G): """ 计算介电常数傅里叶系数 """ if np.linalg.norm(G) < 1e-10: # G=0分量 f = 2 * np.pi * self.r**2 / (self.a**2 * np.sqrt(3)/2) return f * self.eps_a + (1 - f) * self.eps_b else: # 非零G分量 Gr = np.linalg.norm(G) * self.r f = 2 * np.pi * self.r**2 / (self.a**2 * np.sqrt(3)/2) return f * (self.eps_a - self.eps_b) * 2 * sp.j1(Gr) / Gr def build_hamiltonian_TE(self, k, G_vectors, N_plane_waves): """ 构建TE模式的哈密顿量矩阵 """ N = len(G_vectors) H = np.zeros((N, N), dtype=complex) # 构建介电常数倒矩阵 eps_inv_matrix = np.zeros((N, N), dtype=complex) for i in range(N): for j in range(N): G_diff = G_vectors[i] - G_vectors[j] eps_inv_matrix[i, j] = self.fourier_coefficient(G_diff) # 构建波矢矩阵 k_matrix = np.zeros((N, N), dtype=complex) for i in range(N): k_plus_G = k + G_vectors[i] for j in range(N): if i == j: k_plus_G_norm = np.linalg.norm(k_plus_G) k_matrix[i, j] = k_plus_G_norm**2 # TE模式哈密顿量: H = (k+G)^2 * ε^{-1} H = k_matrix @ np.linalg.inv(eps_inv_matrix) return H def calculate_band_structure(self, k_path, num_points=50, N_plane_waves=121): """ 计算能带结构 """ G_max = 4 * np.pi / self.a G_vectors = self.get_reciprocal_lattice_vectors(G_max) bands = [] for k in k_path: H = self.build_hamiltonian_TE(k, G_vectors, N_plane_waves) eigenvalues = np.linalg.eigvalsh(H) # 取实部并排序 freq = np.sqrt(np.real(eigenvalues)) * self.a / (2*np.pi) freq_sorted = np.sort(freq[freq > 0]) bands.append(freq_sorted[:8]) # 取前8个能带 return np.array(bands) ``` ### 2.2 布里渊区路径生成 ```python def generate_hexagonal_brillouin_zone_path(num_points=50): """ 生成六方晶格布里渊区高对称路径: Γ -> M -> K -> Γ """ # 高对称点定义 Gamma = np.array([0.0, 0.0]) M = np.array([2*np.pi/1.0, 0.0]) # 简化表示 K = np.array([2*np.pi/1.0, 2*np.pi/(np.sqrt(3)*1.0)]) # 分段生成路径 path_segments = [] # Γ -> M for t in np.linspace(0, 1, num_points): k = Gamma + t * (M - Gamma) path_segments.append(k) # M -> K for t in np.linspace(0, 1, num_points): k = M + t * (K - M) path_segments.append(k) # K -> Γ for t in np.linspace(0, 1, num_points): k = K + t * (Gamma - K) path_segments.append(k) return np.array(path_segments), ['Γ', 'M', 'K', 'Γ'] ``` ## 3. 完整计算与可视化示例 ```python def main(): # 创建光子晶体实例 pc = HexagonalPhotonicCrystal(r=0.3, eps_a=13.0, eps_b=1.0) # 生成布里渊区路径 k_path, high_sym_points = generate_hexagonal_brillouin_zone_path(30) # 计算能带结构 print("开始计算能带结构...") bands = pc.calculate_band_structure(k_path) # 可视化结果 plt.figure(figsize=(10, 8)) # 绘制能带图 x_axis = np.arange(len(k_path)) for i in range(bands.shape[1]): plt.plot(x_axis, bands[:, i], 'b-', linewidth=1.5) # 标记高对称点 high_sym_indices = [0, 30, 60, 90] for idx in high_sym_indices: plt.axvline(x=idx, color='gray', linestyle='--', alpha=0.7) plt.xticks(high_sym_indices, high_sym_points) plt.xlabel('Wave Vector k') plt.ylabel('Frequency (ωa/2πc)') plt.title('Band Structure of 2D Hexagonal Photonic Crystal (TE Mode)') plt.grid(True, alpha=0.3) plt.ylim(0, 1.2) plt.tight_layout() plt.show() return bands if __name__ == "__main__": bands = main() ``` ## 4. 关键参数优化与验证 ### 4.1 参数敏感性分析 ```python def parameter_sensitivity_analysis(): """ 分析关键参数对能带结构的影响 """ radii = [0.2, 0.3, 0.4] dielectric_contrasts = [(9,1), (13,1), (16,1)] fig, axes = plt.subplots(3, 3, figsize=(15, 12)) for i, r in enumerate(radii): for j, (eps_a, eps_b) in enumerate(dielectric_contrasts): pc = HexagonalPhotonicCrystal(r=r, eps_a=eps_a, eps_b=eps_b) k_path, _ = generate_hexagonal_brillouin_zone_path(20) bands = pc.calculate_band_structure(k_path) ax = axes[i, j] for band_idx in range(6): ax.plot(bands[:, band_idx], 'b-', linewidth=1) ax.set_title(f'r={r}, ε={eps_a}/{eps_b}') ax.set_ylim(0, 1.2) ax.grid(True, alpha=0.3) plt.tight_layout() plt.show() ``` ### 4.2 能带折叠现象验证 ```python def verify_band_folding(): """ 验证能带折叠现象 """ pc = HexagonalPhotonicCrystal(r=0.35, eps_a=12.0, eps_b=1.0) # 扩展的k路径以观察折叠 extended_k_path, labels = generate_extended_brillouin_zone_path() bands = pc.calculate_band_structure(extended_k_path) # 绘制扩展能带图显示折叠 plt.figure(figsize=(12, 8)) for i in range(12): # 显示更多能带以观察折叠 plt.plot(bands[:, i], 'b-', alpha=0.7, linewidth=1) plt.xlabel('Extended Wave Vector Path') plt.ylabel('Frequency (ωa/2πc)') plt.title('Band Folding in Hexagonal Photonic Crystal') plt.grid(True, alpha=0.3) plt.show() ``` ## 5. 性能优化建议 ### 5.1 计算效率提升 ```python def optimized_hamiltonian_construction(k, G_vectors, pc): """ 优化哈密顿量构建，利用对称性和稀疏性 """ N = len(G_vectors) # 利用对称性减少计算量 H_optimized = np.zeros((N, N), dtype=complex) # 预计算常用值 k_norms = [np.linalg.norm(k + G) for G in G_vectors] for i in range(N): for j in range(N): if abs(i - j) < 10: # 利用矩阵的带状特性 G_diff = G_vectors[i] - G_vectors[j] eps_inv = pc.fourier_coefficient(G_diff) H_optimized[i, j] = k_norms[i]**2 * eps_inv return H_optimized ``` 该Python实现基于平面波展开法，能够准确计算二维六方光子晶体的能带结构，并可视化展示光子带隙和能带折叠现象[ref_1][ref_2]。代码包含了完整的物理模型、数值计算方法和可视化工具，适用于光子晶体设计和分析的各种应用场景。