二维三角晶格光子晶体湿度特性研究

Study on the Humidity Characteristics of Two-Dimensional Triangular Lattice Photonic Crysta

摘要: 本文采用时域有限差分方法对二维三角晶格光子晶体的光子禁带湿度特性进行研究，数值模拟了环境湿度对光子禁带的起始波长、终止波长及禁带宽度的影响，并得出了在环境湿度变化的情况下，光子禁带的起始波长、终止波长及禁带宽度随湿度变化的关系。

Abstract: In this paper, the transmission matrix method is used to study the humidity characteristics of the photonic band gap of two-dimensional triangular lattice photonic crystal. The influence of the ambient humidity on the start wavelength, the end wavelength and the band gap width of the photonic band gap is numerically simulated. The relationship between the start wavelength, the end wavelength and the band gap width of the photonic band gap with the change of the ambient humidity is obtained.

1. 引言

湿度对人体、设备以及建筑的围护结构都有非常重要影响。空气湿度大小会影响人体的新陈代谢，影响设备的运行状态，影响建筑围护结构的保温性能等。因此，对湿度的监测是一项重要的工作。目前，使用较多的是电阻式、露点式氯化锂湿度计。这种传统意义上的湿度传感器在使用过程中有很多局限性。如能量损耗大、灵敏度低等。光学传感器是将光子作为信息测量载体，相对于普通的电子传感器，具有更好的抗干扰能力，以及承载信息的能力。

二十世纪八十年代，光子晶体的概念首先被Yablonovich [1] 和S. John [2] 提出。对于目前来说，光子晶体的光子带隙是一个可以将光子作为信息载体的重要研究方向。光子带隙是由于光子晶体内部折射率变化形成的一个频率区间。这个频率区间能够阻隔频率范围内的所有光子的传播。因此，光子在光子晶体中的运动与电子在半导体中的运动是极为相似的。本文涉及的光子带隙即光子禁带。根据光子晶体的这个特性，被应用在很多方面：光子晶体光纤 [3] [4] 、光子晶体滤波器 [5] 、光子晶体偏振光分束器 [6] 、光子晶体激光器 [7] 、光子晶体全光开关 [8] 、光子晶体传感器。本文主要介绍基于三角晶格结构的二维光子晶体，通过建立光子晶体内部结构模型，仿真分析环境中湿度变化，光子晶体光子禁带的变化，分析其感湿原理，验证其作为湿敏元件的可行性。

2. 数学模型

2.1. 二维光子晶体物理特性

本文选择的二维光子晶体具有较大的光子禁带，并且在一定波长条件下会出现完全禁带的效果 [9]。根据二维光子晶体中介质柱的排列方式，大致可以分为正方形晶格、三角形晶格和蜂窝状晶格。其中，三角形晶格的带隙宽度较大且制造难度也相对较低，可以比较容易的获得较大的带隙。二维光子晶体有两种模式的禁带，即在横电模(TE)下的禁带和在横磁模(TM)下的禁带。本文采用三角形晶格结构的二维光子晶体模型如图1所示，对光子禁带进行研究。

Figure 1. Schematic diagram of two-dimensional photonic crystal with triangular lattice

图1. 三角晶格二维光子晶体示意图

三角晶格的晶格周期为a = 850 nm，介质柱的半径为a = 350 nm。该二维光子晶体由折射率不同的介质材料构成，介质柱选择锐钛型TiO₂，介电常数为48，选择气孔率为30%的SiO₂作为背景介质，按照一定的周期排列构成 [10]。因为环境湿度不同，背景介质SiO₂的折射率发生变化，和介质柱组成的光子晶体则会构成起始波长、终止波长和禁带宽度变化的二维光子晶体。本论文通过数据模拟分析不同湿度下，二维光子晶体起始波长、终止波长和禁带宽度变化的趋势。采用不同的环境湿度进行软件数值模拟，得到SiO₂折射率和相对湿度的关系曲线，如图2所示 [10]。

Figure 2. Relationship between refractive index and relative humidity

图2. 折射率与相对湿度关系图

SiO₂薄膜的折射率随着相对湿度增大而增大。环境的湿度变化可以通过材料的折射率变化来反映。并且，由上图可知在30%~100%的相对湿度范围内，折射率的变化较大，有利于获得折射率与相对湿度的定量关系。折射率和介电常数的数量关系如式(1)。

$n^2=\epsilon$ (1)

由式(1)可知，环境的湿度变化可以通过光子晶体的介电常数来反应。目前分析光子晶体能带特性的传统的方法有：平面波展开法 [11] 、传播矩阵法 [12] [13] [14] 、时域有限差分法(Finite difference time domain method，FDTD) [15] 等。平面波展开法不适用与此三角晶格的研究，传播矩阵法不能从物理概念方面理解 [16]。因此，本文选用FDTD法进行研究。

2.2. Maxwell方程及其FDTD形式

Maxwell方程为

$\nabla \times H=\frac{\partial D}{\partial t}+J$ (2)

$\nabla \times E=-\frac{\partial B}{\partial t}-J_m$ (3)

其中， $E$ 为电场强度(V/m)；

$D$ 为电通量密度(C/m²)；

$H$ 为磁场强度(A/m)；

$B$ 为磁通量密度(Wb/m²)；

$J$ 为电流密度(A/m²)；

$J_m$ 为磁流密度(V/m²)。

各向同性线性介质中的本构关系为

$\begin{array}{l}D=\epsilon E\\ B=\mu H\\ J=\sigma E\\ J_m={\sigma }_{m}H\end{array}\}$

其中， $\epsilon$ 表示介电示介电(F/m)；

$\mu$ 表示磁导示磁导(H/m)；

$\sigma$ 表示电示电导(S/m)；

${\sigma }_m$ 表示导示导磁(Ω/m)。

$\sigma$ 和 ${\sigma }_m$ 分别为介质的电损耗和磁损耗，真空中 $\sigma =\text{0}$， ${\sigma }_m=\text{0}$ 以及

$\{\begin{array}{l}\epsilon ={\epsilon }_{\text{0}}=\text{8}\text{.85}\times {\text{10}}^{-12}\text{F}/\text{m}\\ \mu ={\mu }_0=4\pi \times {10}^{-7}\text{H}/\text{m}\end{array}$

对式(1)、(2)进行直角坐标系转化得到式(3)、(4)

$\begin{array}{l}\frac{\partial H_z}{\partial y}-\frac{\partial H_y}{\partial z}=\epsilon \frac{\partial E_x}{\partial t}+\sigma E_x\\ \frac{\partial H_x}{\partial z}-\frac{\partial H_z}{\partial x}=\epsilon \frac{\partial E_y}{\partial t}+\sigma E_y\\ \frac{\partial H_y}{\partial x}-\frac{\partial H_x}{\partial y}=\epsilon \frac{\partial E_z}{\partial t}+\sigma E_z\end{array}\}$ (4)

$\begin{array}{l}\frac{\partial E_z}{\partial y}-\frac{\partial E_y}{\partial z}=-\mu \frac{\partial H_x}{\partial t}-{\sigma }_m{H}_x\\ \frac{\partial E_x}{\partial z}-\frac{\partial E_z}{\partial x}=-\mu \frac{\partial H_y}{\partial t}-{\sigma }_m{H}_y\\ \frac{\partial E_y}{\partial x}-\frac{\partial E_x}{\partial y}=-\mu \frac{\partial H_z}{\partial t}-{\sigma }_m{H}_z\end{array}\}$ (5)

对式(3)、(4)进行FDTD差分离散，令 $f\left(x,y,z,t\right)$ 代表 $E$ 或 $H$ 在直角坐标系中某一分量，在时间或空间域中的离散表示为：

$f\left(x,y,z,t\right)=f\left(i\Delta x,j\Delta y,k\Delta z,n\Delta t\right)=f^n\left(i,j,k\right)$ (6)

对 $f\left(x,y,z,t\right)$ 关于时间和空间的以接偏导数取中心差分近似得式(6)。

$\begin{array}{l}{\frac{\partial f\left(x,y,z,t\right)}{\partial x}|}_{x=i\Delta x}\approx \frac{f^n\left(i+\frac{1}{2},j,k\right)-f^n\left(i-\frac{1}{2},j,k\right)}{\Delta x}\\ {\frac{\partial f\left(x,y,z,t\right)}{\partial y}|}_{y=j\Delta y}\approx \frac{f^n\left(i,j+\frac{1}{2},k\right)-f^n\left(i,j-\frac{1}{2},k\right)}{\Delta y}\\ {\frac{\partial f\left(x,y,z,t\right)}{\partial z}|}_{z=k\Delta z}\approx \frac{f^n\left(i,j,k+\frac{1}{2}\right)-f^n\left(i,j,k-\frac{1}{2}\right)}{\Delta z}\\ {\frac{\partial f\left(x,y,z,t\right)}{\partial t}|}_{t=n\Delta t}\approx \frac{f^{n+1/2}\left(i,j,k\right)-f^{n-1/2}\left(i,j,k\right)}{\Delta t}\end{array}\}$ (7)

在二维情况下，假设与z坐标无关，即 $\partial /\partial z=0$。对于在TE模式下 $H_x=H_y=E_Z=0$，FDTD公式为

$\begin{array}{c}{E}_x^{n+1}\left(i+\frac{1}{2},j\right)=CA\left(m\right)\cdot E_x^n\left(i+\frac{1}{2},j\right)\\ +CB\left(m\right)\cdot \frac{H_z^{n+1/2}\left(i+\frac{1}{2},j+\frac{1}{2}\right)-H_z^{n+1/2}\left(i+\frac{1}{2},j-\frac{1}{2}\right)}{\Delta y}\end{array}$ (8)

$\begin{array}{c}{E}_y^{n+1}\left(i,j+\frac{1}{2}\right)=CA\left(m\right)\cdot E_y^n\left(i,j+\frac{1}{2}\right)\\ +CB\left(m\right)\cdot \frac{H_z^{n+1/2}\left(i+\frac{1}{2},j+\frac{1}{2}\right)-H_z^{n+1/2}\left(i-\frac{1}{2},j+\frac{1}{2}\right)}{\Delta x}\end{array}$ (9)

$\begin{array}{c}{H}_z^{n+1/2}\left(i+\frac{1}{2},j+\frac{1}{2}\right)=CP\left(m\right)\cdot H_z^{n-1/2}\left(i+\frac{1}{2},j+\frac{1}{2}\right)\\ -CQ\left(m\right)\cdot [\frac{E_y^n\left(i+1,j+\frac{1}{2}\right)-E_y^n\left(i,j+\frac{1}{2}\right)}{\Delta x}\\ -\frac{E_y^n\left(i+\frac{1}{2},j+1\right)-E_y^n\left(i+\frac{1}{2},j\right)}{\Delta y}]\end{array}$ (10)

对于在TM模式下， $E_x=E_y=H_z=0$，FDTD公式为

$H_x^{n+1/2}\left(i,j+\frac{1}{2}\right)=CP\left(m\right)\cdot H_x^{n-1/2}\left(i,j+\frac{1}{2}\right)-CQ\left(m\right)\cdot \frac{E_z^n\left(i,j+1\right)-E_z^n\left(i,j\right)}{\Delta y}$ (11)

$H_y^{n+1/2}\left(i+\frac{1}{2},j\right)=CP\left(m\right)\cdot H_y^{n-1/2}\left(i+\frac{1}{2},j\right)+CQ\left(m\right)\cdot \frac{E_z^n\left(i+1,j\right)-E_z^n\left(i,j\right)}{\Delta y}$ (12)

$\begin{array}{c}{E}_z^{n+1}\left(i,j\right)=CA\left(m\right)\cdot E_z^n\left(i,j\right)+CB\left(m\right)\cdot [\frac{H_y^{n+1/2}\left(i+\frac{1}{2},j\right)-H_y^{n+1/2}\left(i-\frac{1}{2},j\right)}{\Delta x}\\ -\frac{H_x^{n+1/2}\left(i,j+\frac{1}{2}\right)-H_x^{n+1/2}\left(i,j-\frac{1}{2}\right)}{\Delta y}]\end{array}$ (13)

上述方程的系数分别为

$CA\left(m\right)=\frac{1-\frac{\sigma \left(m\right)\Delta t}{2\epsilon \left(m\right)}}{1+\frac{\sigma \left(m\right)\Delta t}{2\epsilon \left(m\right)}}$， $CB\left(m\right)=\frac{\frac{\Delta t}{\epsilon \left(m\right)}}{1+\frac{\sigma \left(m\right)\Delta t}{2\epsilon (m)}}$

$CP\left(m\right)=\frac{1-\frac{{\sigma }_m\left(m\right)\Delta t}{2\mu \left(m\right)}}{1+\frac{{\sigma }_m\left(m\right)\Delta t}{2\mu \left(m\right)}}$， $CP\left(m\right)=\frac{\frac{\Delta t}{\mu \left(m\right)}}{1+\frac{{\sigma }_m\left(m\right)\Delta t}{2\mu (m)}}$

上述式子中的m随着FDTD方程变化。

3. 数据分析

理论模型进行数据仿真实验，分析二维三角晶格光子晶体的带隙特征。以本文采用的二维光子晶体的人背景介质是锐钛型TiO₂，将SiO₂作为介质柱，嵌入其中。通过模拟环境湿度变化，测量光子晶体产生的光子禁带的起始波长、终止波长及禁带宽度的大小和变化趋势。并分别分析在TM与TE模式下，湿度与起始波长、终止波长及禁带宽度的定量函数关系，为二维光子晶体作为湿度传感器的应用做理论分析。

根据FDTD模型，模拟湿度在0%与100%时的能带结构图，在TM模式下的能带图如图3所示，在TE模式下的能带图如图4所示。图3、图4的纵坐标为归一化频率，可以根据归一化转化公式(14)转化为波长：

$\frac{a\cdot 2\cdot \pi }{w}=\frac{\lambda }{a}$ (14)

式(14)中a代表周期；w代表角频率； $\lambda$ 代表波长。

(a) 相对湿度0% (b) 相对湿度100%

Figure 3. Band structure of two-dimensional photonic crystal for TM mode

图3. TM模式二维光子晶体能带结构

通过公式(14)求出二维光子晶体的起始波长、终止波长和禁带宽度，并分析环境湿度与它们的关系。

根据FDTD模型，利用图2理论曲线相对湿度和折射率的关系以及公式(2)，分别建立二维光子晶体在TM与TE模式下的起始波长与相对湿度关系曲线。如图5所示，起始波长在湿度0%~100%范围内的变化趋势。

如图5所示，相对湿度越大，三角晶格二维光子晶体的起始波长也越大。在TE模式下光子禁带的起始波长从6218.99 nm增加到6243.85 nm；在TM模式下光子禁带的起始波长从4502.316 nm增加到4508.299 nm。在TE模式下的光子禁带起始波长比TM模式长，并且TE模式下的波长变化也更大一些。在40%~100%之间，相对湿度和起始波长趋于线性关系，经过数据拟合，其拟合公式为：

$fx=p1x+p2$ (15)

TE模式：p1 = 0.3751；p2 = 6206。

TM模式：p1 = 0.09016；p2 = 4499。

(a) 相对湿度0% (b) 相对湿度100%

Figure 4. Band structure of two-dimensional photonic crystal for TE mode

图4. TE模式二维光子晶体能带结构

(a) TM模式 (b) TE模式

Figure 5. Relationship between initial wavelength and relative humidity

图5. 起始波长与相对湿度之间的关系

光子禁带终止波长的测量计算方式与上述起始波长相同，三角晶格二维光子晶体的终止波长随相对湿度的变化趋势如图6所示。

如图6所示，当相对湿度增加时，光子禁带的终止波长也随之增大。在TM模式下，光子禁带的终止波长从4161.799 nm增加到4183.296 nm；在TE模式下，禁带的终止波长从3193.287 nm增加到3251.396 nm。与起始波长相反，光子禁带的终止波长在TM模式下比在TE模式的波长要长，但是TE模式下的终止波长变化更大。在40%~100%之间，相对湿度和终止波长趋于线性关系，经过数据拟合，其拟合公式为：

$fx=p1x+p2$ (16)

TE模式：p1 = 0.88；p2 = 3162。

TM模式：p1 = 0.3248 ；p2 = 4150。

(a) TM模式 (b) TE模式

Figure 6. Relationship between termination wavelength and relative humidity

图6. 终止波长与相对湿度之间的关系

光子禁带的起始波长与终止波长都与相对湿度建立了一次线性函数关系，禁带宽度与相对湿度的函数关系建立方法与之类似，通过理论模型测算三角晶格二维光子晶体的禁带宽度，绘制其与相对湿度的关系曲线图如图7所示，相对湿度与光子禁带宽度的函数关系如式(17)。

(a) TM模式 (b) TE模式

Figure 7. Relationship between band width and relative humidity

图7. 禁带宽度与相对湿度之间的关系

图7表明，在TE模式和TM模式下，三角晶格二维光子晶体的禁带宽度都随着环境的相对湿度的增加而减小，在TE模式下，晶体的禁带宽度从3025.702 nm减少到2992.456 nm；在TM模式下，晶体的禁带宽度从340.516 nm减少到325.003 nm。显而易见，在TM模式下的禁带宽度要比TE模式下小的多，并且在TE模式下，禁带宽度的变化也更加明显。在40%~100%之间，相对湿度和禁带宽度趋于线性关系，经过数据拟合，其拟合公式为：

$fx=p1x+p2$ (17)

TE模式：p1 = −0.5049 ；p2 = 3044。

TM模式：p1 = −0.2346；p2 = 348.8。

4. 结论

本文研究了三角晶格二维光子晶体带隙的起始波长、终止波长和禁带宽度与相对湿度的关系，并通过软件仿真建立了曲线图与线性关系表达式。通过实验数据和曲线图的对比，发现在TE模式下光子禁带的终止波长比TM模式要短；禁带的起始波长和禁带宽度在TE模式下比TM模式下要长；TE模式下起始波长、终止波长和禁带宽度的变化都要比TM模式下大。并且，TM所产生的禁带落在TE模式的禁带中，说明了TM模式产生的禁带区域属于完全光子禁带。以上研究结果可以表明基于三角晶格二维光子晶体的湿度传感器在原理上是可行的，它相对于其他原理的传感器具有非常大的优势，光子晶体湿度传感器是未来湿度传感器的发展方向。

基金项目

宁波工程学院王伟明基金资助项目(2017010)。