FDTD中的光子晶体

前言

光子晶体是利用材料的空间周期性分布来操控光的一种方法。在特定条件下，光子晶体能形成光子带隙，阻止特定频率的光子传播。这种效应类似于半导体晶体对电子行为的影响，由此，人们可以通过光子晶体制造的器件来控制光子运动，进而制造光子作为信息载体的计算机。 为了计算光子晶体中的光学模式，我们需要在周期性介电介质条件下求解麦克斯韦方程组，尽管这项任务看似简单，但实际上我们无法获得二维或者三维周期性晶格的解析解。因此，数值计算方法被提了出来，例如有：平面波展开法、传输矩阵法、时域有限差分法等等。其中，时域有限差分法是电磁场数值计算的经典方法之一，易于实现且计算效率高，能够直观地获得能带结构、带隙大小和传输场等重要信息。

研究背景与基本原理

光子晶体的研究历史

虽然早在1887年，英国的Lord Rayleigh就揭示了一维周期性结构的光学特性。但光子晶体的概念直到1987年才由Eli Yablonovitch和Sajeev John两位科学家团队分别提出。此后，关于光子晶体的学术论文数量呈现出几何级数上升的趋势[1]。然而，由于制作光学尺寸的光子晶体难度较大，早期的研究主要集中在理论探索以及微波级光子晶体的制造上。直到1996年，Thomas Krauss利用现有的半导体技术成功研制出第一个光学尺寸的二维光子晶体。进入21世纪以来，伴随着凝聚态物理中拓扑绝缘体的发展，光子晶体拓扑绝缘体和光子外尔晶体相继实现[2] 。光子晶体的拓扑性质成为当前的研究热点。

自然界中也存在着天然的光子晶体，这类光子晶体结构能够呈现出各种鲜艳的颜色，与依靠色素反射或吸收光来呈现颜色的鲜花和绿叶不同，它们是依靠微观结构的周期性排列来实现。在这种方式下，只要微观结构不被破环，其颜色就能持久如新。例如，自然界中的蝴蝶翅膀[3]、孔雀的羽毛[4]等等，都展现出鲜艳的色彩，这些颜色正是源于尺寸与光的波长相近的周期性排列的微观光子晶体。

光子晶体的基本原理

根据介质材料周期排列的方式，可分为一维、二维、三维光子晶体。一维光子晶体可以认为是由两种不同介电材料交替叠加形成的多层膜结构。而布拉格光栅是其中的一种特殊示例。这种简单结构也拥有光子晶体的最重要的特征，如光子禁带，满足布拉格条件的光将被反射。在二维光子晶体中，波长处在禁带中的光在平面内任何方向入射均会被反射。引入点缺陷或者线缺陷时，可以将光局域在该缺陷处，或者将光限制在该缺陷内传播。可以利用这种缺陷态可以制作各种光波导器件。三维光子晶体则有着更为复杂的晶体结构和缺陷构造，人们依然最关心那些存在完全带隙的晶体结构。目前研究表明，拥有大完全带隙的光子晶体都和金刚石结构有着紧密的关系[5]。

综上，光子晶体本质上是一种特殊的混合电介质结构，需要结合电磁学以及固体物理学等知识来研究电磁波在其中的传播现象。基于这种理论基础，可以将麦克斯韦方程组转化为线性哈密顿本征值问题，这将使得其与薛定谔方程极为相似。

电磁波在空气传播时，其色散曲线可以简单表示为 $\omega =c|k|$，这条直线即为光锥线。普通电介质的色散曲线都位于光锥线下方，其方程为 $\omega =c|k|/\sqrt{{\epsilon }_r}$ 。而光子晶体的介电常数分布类似于传统原子、分子晶体，具有周期性和各种对称性。平移对称性使得 $\epsilon (\mathbf{r})=\epsilon (\mathbf{r}\pm \mathbf{R})$
， $\mathbf{R}$ 为晶格矢量，其用来度量晶格的空间位置，表示为 $\mathbf{R}=l{\mathbf{a}}_1+m{\mathbf{a}}_2+n{\mathbf{a}}_3$ ， $\mathbf{a}$ 为晶格基矢， $l,m,n$ 为任意整数。由布洛赫理论即可得到本征向量
$${\mathbf{H}}_{\mathbf{k}}(\mathbf{r})=e^{i\mathbf{k}\cdot \mathbf{r}}{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})$$
其中 ${\mathbf{u}}_{\mathbf{k}}(\mathbf{r})$ 为与晶格相关的周期函数， ${\mathbf{u}}_{\mathbf{k}}(\mathbf{r})={\mathbf{u}}_{\mathbf{k}}(\mathbf{r}+\mathbf{R})$ 对于任意的晶格矢量都成立。 而在波矢空间， $\mathbf{k}$ 表示的本征态和 $\mathbf{k}+\mathbf{G}$ 表示的本征态是相同的（ $\mathbf{G}$ 为倒格矢，用于度量波矢的变化，在倒格矢空间中表示为 $\mathbf{G}=k_1{\mathbf{b}}_1+k_2{\mathbf{b}}_2+k_3{\mathbf{b}}_3$ ， $\mathbf{b}$ 为倒格矢基矢）。晶格基矢与倒格矢基矢之间满足关系 ${\mathbf{a}}_{\mathbf{i}}{\mathbf{b}}_{\mathbf{j}}=2\pi {\delta }_{i,j}$ 。有关倒格矢的计算可参考固体物理相关资料[6]，在此不再详细展开。由上述可知，光子的本征态存在许多冗余项，因此只需要考虑倒格矢空间内的有限区域，即最接近 $\mathbf{k}=0$ 的区域，即（第一）布里渊区。

除此之外，光子晶体还拥有旋转对称性、镜像反射对称性以及反演对称性。对于一个给定的光子晶体，只需要将求解问题限制在单个元胞内即可。在布里渊区中无法通过对称操作得到的最小部分，即不可约布里渊区。如上图中(e)所示，简单正方晶格光子晶体的不可约布里渊区为三角形，其顶点分别位于中心 $\Gamma$ 、角 $M$ 和一条边 $X$ 上，其中包含了布里渊区的所有信息。

光子晶体的能带分析

以上文中经典的二维正方圆柱光子晶体为例，以下展示两种求解光子晶体能带的方法。

使用FDTD求解光子晶体的带隙

时域有限差分法通过差分时域麦克斯韦方程组，获得稳定的电磁场时域解，同时也可用来提取频率本征值。在频域本征求解器中，可以观察结构对短脉冲的响应，响应谱中的峰值即为本征频率。这种方法还可用于识别谐振或泄漏模式，因为峰宽与损耗率相关。在实际仿真中，人们不仅仅可以在傅里叶变换中寻找峰，使用一些复杂的信号处理技术甚至比傅里叶不确定性原理得到的结果更精确。对于光子晶体来说，可以通过施加布洛赫周期边界条件计算能带结构。这种方法可以一次性获得多个本征频率，它可以预防伪解，并在计算少数几个本征频率时更快。

以下使用国产SimWork FDTD仿真软件[7]进行计算，下图为在软件中建立二维光子晶体元胞来求解能带结构。黄色箭头表示随机方向随机位置的偶极子源，用于激发所有的模式。布洛赫边界将使每次仿真后得到特定布洛赫波矢的频率响应，在足够的仿真时间后，不满足存在条件的模式光将消散。绿色的叉号表示随机位置的时间监视器，可确保捕获所有模式的电场，收集并相加，计算这些场的谐振频率，从而得到多个本征频率，使用布洛赫波矢扫描即可获得完整的能带图。

由于材料拟合，网格划分等因素影响，仿真结果可能需要进行收敛性测试，以获得稳定的解。最终扫描结果如下：

可以观察到TM模式下的能带结构有完整的光子带隙，而TE模式则仅有部分光子带隙，意味着该光子晶体能阻断处于光子禁带频率中的任意TM模式光传播，而对TE模式的光仅能反射部分角度的光。以上结果可在其官网案例库2D正方晶格能带结构中找到详细描述。

平面波展开法

平面波展开法是一种经典的求解光子能带的方法，它将周期性系统的麦克斯韦问题转变为频域本征求解问题。拥有求解完整的光子晶体能带以及直接得到光子的色散关系等优点，但它对于现实中较为复杂的光子晶体来说，通常需要较大的矩阵计算，结果才能收敛，这往往需要耗费较大的计算资源。在此，可以用来验证FDTD计算结果。 对于二维正方圆柱光子晶体，从上文中可以知道，其电场在z方向上周期性势场可以表示为：$E_z=U_{x,y}{e}^{i(k_x\cdot x+k_y\cdot y)}$，其中 $U_{x,y}$为周期函数， $U_{x,y}={\sum }_{m,n}{A}_{mn}{e}^{i(\frac{2\pi m}{a}x+\frac{2\pi n}{a}y)}$，其中 $A_{mn}(m,n\in [-N,N])$ 为傅里叶展开系数，N为展开级数。考虑E偏振，可以得到对应傅里叶展开的磁场分量 $H_x,H_y$ 。

结合上文中的图5(b)所示，其材料的相对介电常数可表示为
$$\epsilon ={\epsilon }_b+({\epsilon }_a-{\epsilon }_b)f_{(x,y)},f_{(x,y)}=\{\begin{array}{l}1,x^2+y^2⩽R^2\\ 0,其他\end{array}$$
将其作傅里叶展开
$${\epsilon }_{(x,y)}=\sum _{m,n}{r}_{mn}{e}^{i(\frac{2\pi m}{a}x+\frac{2\pi n}{a}y)}$$
并代入麦克斯韦方程组中，整理可得：
$$\{\begin{array}{l}(\frac{2\pi n}{a}+k_y)A_{mn}=\omega {\mu }_0{B}_{mn}\\ -(\frac{2\pi m}{a}+k_x)A_{mn}=\omega {\mu }_0{C}_{mn}\\ (\frac{2\pi m}{a}+k_x)C_{mn}-(\frac{2\pi n}{a}+k_y)B_{mn}=-\omega {\epsilon }_0{\sum }_{m^{\prime },n^{\prime }}{A}_{m^{\prime }{n}^{\prime }}{r}_{m-m^{\prime },n-n^{\prime }}\end{array}$$
该式是无穷项的，需要进行平面波截断，并引入矩阵：
$$\{\begin{array}{l}{K}_x(l,l^{\prime })=(\frac{2\pi m}{a}+k_x){\delta }_{m,m^{\prime }}{\delta }_{n,n^{\prime }}\\ K_y(l,l^{\prime })=(\frac{2\pi n}{a}+k_y){\delta }_{m,m^{\prime }}{\delta }_{n,n^{\prime }}\\ R(l,l^{\prime })=r_{m-m^{\prime },n-n^{\prime }}\end{array}$$
故有：$R^{-1}(K_x^2+K_y^2)X=\frac{{\omega }_E^2}{c^2}X$

考虑H偏振：$(K_x{R}^{-1}{K}_x+K_y{R}^{-1}{K}_y)X=\frac{{\omega }_H^2}{c^2}X$

此时结果如下，该仿真计算结果与上文中FDTD算法计算结果完全一致。

FDTD在光子晶体研究中的应用

理论上说，FDTD可以利用bloch波矢扫描计算得到任意光子晶体的能带结构，同时可以进行光子晶体的模式求解以及传输场的仿真。上文当中提及的光子晶体均可以在FDTD中仿真得到。以下为一些复杂光子晶体案例展示，这些案例均可以在SimWork FDTD案例库中找到对应仿真工程。

光子晶体光纤已经发展出多种构型，但其本质依然是基于点缺陷的光子晶体，这种光子晶体拥有极佳的聚焦和束缚效果，减小了光信号的在远距离传输时的损耗。下图为一种光子晶体光纤的高阶模式场图。

基于线缺陷的光开关也已经发展地十分成熟，有多种不同的变型。以下为在自准直光子晶体中引入线缺陷后的传输场分布图。入射光沿特定方向入射到该光子晶体中，光束几乎不发生衍射沿直线传播，而引入线缺陷将导致自准直光束发生弯曲和分裂，此时入射光将在缺陷界面处发生完全反射。

三维光子晶体基本的面心立方以及体心立方晶格都没有完全带隙。而木堆晶体结构是人们研究制造的第一个拥有完全带隙的微米级三维光子晶体[8]。下图为进行波矢扫描得到的木堆光子晶体的能带结构图。

结语

随着光子学领域的迅速发展，越来越多复杂结构的光子晶体被研究发现出来。数值计算方法在光子晶体的理论研究中扮演着重要角色，人们不断追求更快、更高效的数值计算方法。FDTD算法稳定地求解电磁场时域解并提取频域本征值，为光子晶体的理论研究提供了有效工具。在这个快速发展的领域中，不断探索和优化数值计算方法，将有助于推动光子晶体研究取得更深入的进展。

参考

1. Wikipedia contributors. Photonic crystal. Wikipedia, The Free Encyclopedia, Photonic crystal, 2024.
   https://en.wikipedia.org/w/index.php?title=Photonic crystal&oldid=1224522534
2. Lu, L., Fu, L., Joannopoulos, J. et al. Weyl points and line nodes in gyroid photonic crystals[J]. Nature Photon 7, 294-299 2013.
3. Stavenga D G , Stowe S , Siebke K et al. Butterfly wing colours: scale beads make white pierid wings brighter[J]. Proceedings ofthe Royal Society B: Biological Sciences, 2004.
4. Zi J , Yu X , Li Y ,et al. Coloration strategies in peacock feathers[J]. Proceedings of the National Academy of Sciences of the United States of America, 2003.
5. Joannopoulos J , Johnson S , Winn J et al. Photonic Crystals: Molding the Flow of Light - Second Edition[M]. 2008.
6. 黄昆,韩汝琦改.固体物理学[M].高等教育出版社,1988.
7. SimWork application. https://www.simworks.com/zh-CN/customer-download
8. Lin S, Fleming J G, Hetherington D L, et al. A three-dimensional photonic crystal operating at infrared wavelengths[J]. Nature,1998.