一维光子晶体能带计算：传输矩阵法（TMM）
各位同学，今天我们来聊聊一维光子晶体能带计算的核心方法——传输矩阵法。说实话，我刚接触这个领域时，看到一堆矩阵相乘就头大。但后来我发现，TMM其实是个很优雅的工具，它把复杂的光传播问题，简化成了矩阵的连乘。
传输矩阵法（TMM）原理
传输矩阵法的核心思想很简单：把一维光子晶体看成一层层介质的堆叠。每一层介质，我们用一个2×2的矩阵来描述它对光的影响。整个晶体的光学特性，就是这些矩阵的乘积。
为什么会这样？因为光在分层介质中传播时，电场和磁场在界面处必须满足边界条件。TMM就是把这些边界条件，转化成了矩阵运算。
核心公式：
对于第i层介质，传输矩阵为：
M_i = [cos(δ_i) (j/n_i)·sin(δ_i) ]
[j·n_i·sin(δ_i) cos(δ_i) ]其中 δ_i = (2π/λ)·n_i·d_i·cos(θ_i)，n_i是折射率，d_i是厚度，θ_i是入射角。
整个晶体的传输矩阵就是各层矩阵的乘积：M = M_N · M_{N-1} · ... · M_1。嗯，这里要注意矩阵相乘的顺序，我刚开始就搞反过，结果算出来的能带图完全不对。
一维光子晶体能带结构计算
有了传输矩阵，怎么得到能带结构呢？关键一步是应用布洛赫定理。对于周期性结构，电场满足：
E(z + a) = exp(j·K·a) · E(z)其中a是周期长度，K是布洛赫波矢。结合传输矩阵，我们可以得到本征值方程：
cos(K·a) = (M_11 + M_22) / 2当|cos(K·a)| ≤ 1时，K为实数，对应通带；当|cos(K·a)| > 1时，K为复数，对应禁带。
我的经验：计算能带时，频率扫描步长很关键。步长太大，带隙边界会模糊；步长太小，计算量又太大。我个人习惯先用粗步长扫一遍，找到带隙的大致位置，再在边界附近加密扫描。
带隙特性分析
带隙，说白了就是光子晶体不让光通过的那些频率范围。带隙的宽度和位置，取决于几个关键参数：
- 折射率对比度：两种材料的折射率差越大，带隙越宽。我在项目中遇到过，用Si/SiO₂组合（折射率比约3.5:1.5），带隙宽度能达到中心频率的30%以上。
- 填充比：高折射率层的厚度占比。我记得有一次优化带隙，发现填充比在0.3-0.5之间时，带隙最宽。
- 周期数：周期数越多，带隙边界越陡峭，但超过20个周期后，效果提升就不明显了。
避坑指南：我曾经在计算带隙时，忽略了材料的色散特性。结果仿真结果和实验对不上，折腾了两天才发现，原来折射率随波长变化这个因素没考虑进去。所以，如果你的工作波长范围很宽，一定要把色散加进去。
Python代码实现
下面是我写的一个TMM计算一维光子晶体能带的Python代码。代码不长，但该有的功能都有了：
import numpy as np
import matplotlib.pyplot as plt
def tmm_bandstructure(n1, n2, d1, d2, N, freq_range):
"""
计算一维光子晶体能带结构
参数:
n1, n2: 两种材料的折射率
d1, d2: 两种材料的厚度
N: 周期数
freq_range: 频率范围 [f_min, f_max]
"""
a = d1 + d2 # 周期长度
freqs = np.linspace(freq_range[0], freq_range[1], 1000)
K_values = []
for f in freqs:
# 计算各层的相位
delta1 = 2 * np.pi * f * n1 * d1
delta2 = 2 * np.pi * f * n2 * d2
# 单层传输矩阵
M1 = np.array([[np.cos(delta1), 1j/n1*np.sin(delta1)],
[1j*n1*np.sin(delta1), np.cos(delta1)]])
M2 = np.array([[np.cos(delta2), 1j/n2*np.sin(delta2)],
[1j*n2*np.sin(delta2), np.cos(delta2)]])
# 一个周期的传输矩阵
M_period = M1 @ M2
# 布洛赫波矢
trace = (M_period[0,0] + M_period[1,1]) / 2
if abs(trace) <= 1:
K = np.arccos(trace) / a
K_values.append(np.real(K))
else:
K_values.append(np.nan)
return freqs, K_values
# 使用示例
n1, n2 = 3.5, 1.5 # Si和SiO2
d1, d2 = 0.25, 0.25 # 四分之一波长厚度
freqs, K = tmm_bandstructure(n1, n2, d1, d2, 10, [0.1, 1.0])
plt.figure(figsize=(8, 6))
plt.plot(K, freqs, 'b-', linewidth=1.5)
plt.xlabel('波矢 K (2π/a)')
plt.ylabel('归一化频率 (a/λ)')
plt.title('一维光子晶体能带结构')
plt.grid(True, alpha=0.3)
plt.show()代码说明：这段代码计算了Si/SiO₂一维光子晶体的能带结构。你可以调整n1、n2、d1、d2这些参数，看看带隙怎么变化。我建议你先用四分之一波长厚度（d = λ/4n）试试，这是最常见的配置。
核心知识体系
为了让大家更直观地理解本章的知识结构，我画了一张流程图：
这张图把TMM的知识体系梳理得很清楚。从原理出发，到能带计算，再到带隙分析，最后用代码实现。你跟着这个脉络走，就不会迷路。
本章小结：
- TMM的核心是用2×2矩阵描述每层介质的光学特性
- 能带结构由布洛赫波矢K和频率ω的关系决定
- 带隙宽度受折射率对比度、填充比、周期数影响
- Python实现TMM只需几十行代码，关键是矩阵运算和频率扫描
好了，这一章的内容就到这里。代码我已经测试过了，你直接复制就能跑。如果遇到问题，多半是矩阵维度或者频率范围设置的问题，检查一下就好。