二维色散介质光子晶体的能带结构求解问题研究

In this study, the characteristics of two-dimensional dispersive dielectric photonic crystals are studied via a numerical approach that combines the finite-element method (FEM) with Newton's iterative method. For a dispersive photonic crystal whose permittivity is dependent on the frequency, the band structure problem is formulated as a nonlinear eigenvalue problem. In this study, first, a discrete variational formulation is derived from Maxwell's equations based on the FEM. Thereafter, by selecting approximate solutions as the initial values for the iteration, this nonlinear problem is solved for different values of the wave vector k based on Newton's iterative method. Consequently, the band structure of dispersive dielectric photonic crystals is obtained numerically. Several dispersive photonic crystals in the transverse electric (TE) and transverse magnetic (TM) modes are investigated. The numerical results reveal that the proposed method is effective for dispersive photonic crystals.

1　引言

光子晶体是不同折射率材料呈周期排列而形成的人工微结构，其出现为在亚波长尺度上实现光传输操控以及光与物质的相互作用调控提供了强大而有效的手段^［1］。光子晶体的概念于1987年被提出，其因特殊的控光能力被广泛应用在很多领域，如光子芯片、高品质光学谐振腔、滤波器等^［2］。

色散介质光子晶体是光子晶体众多种类中的一种，其是由色散介质和一般介质组成的周期性空间分布的人工结构^［3］。研究色散介质光子晶体主要考虑二维色散介质光子晶体的能带结构，对开发利用光子晶体，制造特殊的光学器件等都具有十分重要的意义。区别于普通光子晶体，色散介质光子晶体的介电常数或磁导率是关于频率的函数。此时，光子晶体能带结构的求解是一个非线性特征值问题。当前求解此类问题的方法有很多，例如，改进的平面波展开法、时域有限差分与频域有限差分法、光谱指示方法（SIM）、Dirichlet-to-Neuman映射法等，这些计算方法也各有特点。在改进的平面波展开法^［4］中，特征模是否收敛是需要着重考虑的一个问题，此方法要保证有较高精度的初始值；频域有限差分法^［5］在求解色散问题时，需要对网格进行剖分，这时通常会用到矩形网格剖分，但是这样很容易引入较大的误差，使得求解的准确度下降。时域有限差分法^［6］在求解色散光子晶体能带结构时，需分步进行，在计算出特征模态后要对其作进一步分析，求解过程较为繁琐。近年来，还有一些新的方法用来求解光子晶体色散问题，如SIM方法^［7］，该方法在求解此问题时，首先利用有限元方法将问题离散，随后利用SIM-H算法求解特征值，该方法为求解Hermitian矩阵的非线性特征值问题提供了一个新思路。此外，Dirichlet-to-Neuman映射法^［8-9］在求解色散关系时，将能带结构求解问题转化为一个二次非线性特征值方程，通过给定频率 $ω$ 求解特征值方程，以波矢 $k$ 作为特征值，得到关系 $ω (k)$ ，从而获得能带结构分布。

本文给出了一种计算二维色散介质光子晶体能带结构的方法。首先，基于电磁波传播的控制方程组推导出变分形式，接着将其转化为一个非线性特征值问题，利用有限元方法将问题离散。离散之后，选择多个合适的解作为初始值，通过扫描不可约布里渊区域的边界获得 $k$ 值，基于牛顿法对不同的 $k$ 值求解该非线性问题，得到色散关系 $ω (k)$ ，从而获得色散介质光子晶体的能带结构。利用COMSOL软件，通过其自带的参数化扫描，对不同的 $k$ 值计算出可能存在能带的频率 $ω$ ，并将其作为牛顿法求解的初始值。将有限元离散与牛顿迭代相结合，利用有限元方法得到问题的离散形式，随后基于牛顿迭代法求解非线性问题，最终求得色散介质光子晶体的能带结构。

2　光子晶体原胞与第一布里渊区域

光子晶体中的电磁场计算在倒易空间中进行，只需计算第一布里渊区域中的电磁场即可获悉整个光子晶体的场分布。给出两种二维光子晶体原胞和其对应的第一布里渊区域，如图1所示。图1（a）为正方晶格光子晶体原胞，晶格基矢 $a_{1} = a e_{x}, a_{2} = a e_{y}$ ，其中 $a$ 为晶格常数， $e_{x} 、 e_{y}$ 为单位向量，第一布里渊区域为 ${[- π / a, π / a]}^{2}$ ，如图2（a）所示。以顶点 $Γ = 2 π / a (0,0)$ ， $X = 2 π / a (1 / 2,0)$ ， $M = 2 π / a (1 / 2,1 / 2)$ 围成的三角形区域 $Γ X M$ 表示其不可约布里渊区域。图1（b）为三角晶格光子晶体原胞，其晶格基矢为 $a_{1} = (1 / 2, \sqrt[]{3} / 2) a e_{x},$ $a_{2} = (1 / 2, - \sqrt[]{3} / 2) a e_{y}$ ，如图2（b）所示。以顶点 $Γ = (0,0) 2 π / a$ ， $M =$ $(0,1 / \sqrt[]{3}) 2 π / a$ ， $X = (1 / 3,1 / \sqrt[]{3}) 2 π / a$ 围成的三角形区域 $Γ X M$ 是三角晶格的不可约布里渊区域。

图 1. 光子晶体原胞。（a）正方晶格；（b）三角晶格

Fig. 1. Photonic crystal unit cell. (a) Square lattice; (b) triangular lattice

下载图片查看所有图片

图 2. 第一布里渊区域。（a）正方晶格；（b）三角晶格

Fig. 2. First Brillouin zone. (a) Square lattice; (b) triangular lattice

下载图片查看所有图片

3　理论模型

以正方晶格为例（原胞区域 $D_{0}$ 如图1（a）所示），推导光子晶体特征值问题。电磁波在二维光子晶体结构中存在不同的传播模式，若只考虑电磁波传播的横电（TE）模式和横磁（TM）模式，二维光子晶体的标量波动方程为

\nabla \cdot p \nabla ϕ + (ω^{2} / c^{2}) q ϕ = 0

，（1）

式中： $\nabla$ 为梯度算子；c为真空光速。TE模式下， $ϕ = H_{z}$ ， $p = 1 / ε (r, ω)$ ， $q = 1$ 。TM模式下， $ϕ = E_{z}$ ， $p = 1$ ， $q = ε (r, ω)$ 。引入波矢 $k \in Κ$ ， $Κ$ 为正方晶格的第一布里渊区域。根据Bloch原理，函数 $ϕ (r)$ 可描述为 $ϕ (r) = e x p (i k \cdot r) u (r)$ ， $u (r)$ 为 $ϕ$ 的Floquet变换，在 $x$ 方向、 $y$ 方向都满足周期性。

由 $\nabla ϕ (r) = \nabla (e x p (i k \cdot r) u (r)) = e x p (i k \cdot r) (\nabla + i k) u (r)$ ，代入式（1），可得TM模式下：

(\nabla + i k) \cdot (\nabla + i k) u + (ω^{2} / c^{2}) ε (r, ω) u = 0

，（2）

式中： $ε$ 为介电常数；u为待求未知量，表示各点处的场值函数。

TE模式下：

(\nabla + i k) \cdot [\frac{1}{ε (r, ω)} (\nabla + i k) u] + \frac{ω^{2}}{c^{2}} u = 0

，（3）

式（2）和式（3）综合可得：

\{\begin{array}{l} (\nabla + i k) \cdot [ρ (\nabla + i k) u] + ω^{2} κ u = 0, i n D_{0} \\ u |_{Γ_{1}} = u |_{Γ_{3}}, u |_{Γ_{2}} = u |_{Γ 4}, o n \cup_{i = 1}^{4} Γ_{i} \end{array}

，（4）

式中： $ρ = \{\begin{array}{l} 1, T M \\ 1 / ε (r, ω), T E \end{array}; κ = \{\begin{array}{l} {\tilde{k_{0}}}^{2} ε (r, ω), T M \\ {\tilde{k_{0}}}^{2}, T E \end{array}$ ， $\tilde{k_{0}} = 1 / c$ 。

定义Sobolev空间 $H^{1} (D_{0})$ 上满足周期边界条件的函数子空间 $H_{p}^{1} (D_{0})$ ：

H_{p}^{1} (D_{0}) = \{v \in H^{1} (D_{0}) | v |_{Γ_{1}} = v |_{Γ_{3}}, v |_{Γ_{2}} = v |_{Γ_{4}}\}

，（5）

任取 $v \in H_{p}^{1} (D_{0})$ ，对式（4）中第一个方程两边同乘以测试函数 $v$ ，得到变分方程：

\int_{D_{0}} ρ (\nabla + i k) u \cdot \bar{(\nabla + i k) v} d x d y = \int_{D_{0}} ω^{2} κ u \bar{v} d x d y, \forall v \in H_{p}^{1} (D_{0})

。（6）

定义算子函数 $T (ω) : S \to ℒ (H_{p}^{1} (D_{0}), H_{p}^{1} (D_{0}))$ ，频率 $ω$ 为复数， $S \subset C$ 为复平面上的一个有界闭集， $ℒ (H_{p}^{1} (D_{0}), H_{p}^{1} (D_{0}))$ 为从 $H_{p}^{1} (D_{0})$ 到 $H_{p}^{1} (D_{0})$ 上的有界线性算子空间。进一步 $T (ω)$ 满足：

(T (ω) u, v) = \int_{D_{0}} ρ (\nabla u \cdot \nabla \bar{v} + i k u \cdot \nabla \bar{v} - i \nabla u \cdot k \bar{v} + {|k|}^{2} u \bar{v}) - ω^{2} κ u \bar{v} d x d y, \forall v \in V

，（7）

从而将式（4）转化为下列特征值问题，求 $ω \in S, u \in H_{p}^{1} (D_{0})$ ，使得：

T (ω) u = 0

。（8）

4　数值解法

类似文献［7］，利用有限元方法离散 $T (ω)$ ，其中 $ω \in S$ 。对区域 $D_{0}$ 进行三角网格剖分得到 ${#x1D4AF}_{h}$ ，且 ${#x1D4AF}_{h}$ 是对称剖分，即 $Γ_{1}$ 和 $Γ_{3}$ 上的节点、 $Γ_{2}$ 和 $Γ_{4}$ 上的节点分别关于区域中心轴对称。构造试探函数空间 $H_{p}^{1} (D_{0})$ 的一个有限元空间 $V_{h}$ ，然后求 $ω$ 以及 $u_{h} \in V_{h}$ ，使得：

\int_{D_{0}} ρ (\nabla u_{h} \cdot \nabla {\bar{v}}_{h} + i k u_{h} \cdot \nabla {\bar{v}}_{h} - i \nabla u_{h} \cdot k {\bar{v}}_{h} + {|k|}^{2} u_{h} {\bar{v}}_{h}) d x d y = \int_{D_{0}} ω^{2} κ u_{h} {\bar{v}}_{h} d x d y, \forall v_{h} \in V_{h}

。（9）

设 $ϕ_{i}, i = 1,2, \dots, N$ 是空间 $V_{h}$ 的基函数， $u_{h}$ 可以表示为

u_{h} (x, y) = \sum_{j = 1}^{N} u_{h, j} ϕ_{j} (x, y),

（10）

式中： $u_{h, j}$ 为待定的标量； $ϕ_{j} (x, y)$ 为第 $j$ 个节点对应的线性基函数。用 $(x_{j}, y_{j})$ 表示第 $j$ 个节点的坐标，且基函数 $ϕ_{j} (x, y)$ 的性质：

ϕ_{j} (x_{i}, y_{i}) = \{\begin{array}{l} 1, i = j \\ 0, i \neq j \end{array}

。（11）

因此，得到有限元变分的离散形式：

\begin{matrix} \int_{D_{0}} ρ [(\nabla (u_{h, j} ϕ_{j}) \cdot \nabla {\bar{v}}_{h}) + i k u_{h, j} ϕ_{j} \cdot \\ \nabla {\bar{v}}_{h} - i \nabla (u_{h, j} ϕ_{j}) \cdot k {\bar{v}}_{h} + {|k|}^{2} u_{h, j} ϕ_{j} \cdot {\bar{v}}_{h}] d x d y = \\ \int_{D_{0}} ω^{2} κ \sum_{j = 1}^{N} u_{h, j} ϕ_{j} \cdot {\bar{v}}_{h} d x d y, \forall v_{h} \in V_{h} \end{matrix}

。（12）

对给定的波矢 $k$ ， $k = {(α, β)}^{T}$ ，将式（12）中的试探函数 $v_{h}$ 选取为基函数，采用N个基函数，即可得到离散方程组系统：

\begin{matrix} ρ u_{h, j} \int_{D_{0}} [(\nabla ϕ_{j} \cdot \nabla \bar{ϕ_{i}}) + i k ϕ_{j} \cdot \nabla \bar{ϕ_{i}} - i \nabla ϕ_{j} k \\ \bar{ϕ_{i}} + {|k|}^{2} ϕ_{j} \bar{ϕ_{i}}] d x d y = ω^{2} κ \sum_{j = 1}^{N} u_{h, j} \int_{D_{0}} ϕ_{j} \cdot \bar{ϕ_{i}} d x d y, \\ i = 1,2, \dots, N \end{matrix}

。（13）

设 $M_{ρ}^{h}, M^{h}, S^{h}, A_{1}^{h}, A_{2}^{h}$ 都是 $N \times N$ 矩阵：

M_{ρ (i j)}^{h} = ρ \int_{D_{0}} ϕ_{j} \bar{ϕ_{i}} d x d y, M_{(i j)}^{h} = \int_{D_{0}} ϕ_{j} \bar{ϕ_{i}} d x d y,

（14）

S_{(i j)}^{h} = ρ \int_{D_{0}} \nabla ϕ_{j} \nabla \bar{ϕ_{i}} d x d y, A_{1 (i j)}^{h} = ρ κ \int_{D_{0}} k ϕ_{j} \bar{ϕ_{i}} d x d y,

（15）

A_{2 (i j)}^{h} = ρ κ \int_{D_{0}} \nabla ϕ_{j} \cdot k \bar{ϕ_{i}} d x d y, i, j = 1,2, \dots, N

。（16）

故原问题经有限元离散可写成矩阵形式的求解问题：给定波矢 $k$ ，求解 $ω$ 以及N维向量 $u$ ，满足：

T^{h} (ω) u = 0

，（17）

式中： $T^{h} (ω) = S^{h} + i A_{1}^{h} - i A_{2}^{h} + {|k|}^{2} M_{ρ}^{h} - ω^{2} κ M^{h}$ 。

基于牛顿法在 $S (S \subset C)$ 上求解 $T^{h}$ 的所有特征值。将式（17）中的 $T^{h} (ω)$ 简记作 $T (ω)$ 。类似文献［10］中的推导，对给定的波矢 $k$ ，设 $(u^{*}, ω^{*})$ 是式（17）的一组真实解，取 $(u^{i}, ω^{i})$ 作为初始值，由泰勒定理得：

\begin{matrix} 0 = T (ω^{*}) u^{*} = T (ω^{*}) u^{i} + \\ T (ω^{*}) (u^{*} - u^{i}) = T (ω^{i}) u^{i} + T (ω^{i}) (u^{*} - u^{i}) + (ω^{*} - ω^{i}) T_{ω} (ω^{i}) u^{i} + \dots, \end{matrix}

（18）

忽略高阶部分，用 $(u^{i + 1}, ω^{i + 1})$ 代替 $(u^{*}, ω^{*})$ 得到：

T (ω^{i}) y^{i + 1} = - T_{ω} (ω^{i}) u^{i}

，（19）

式中： $y^{i + 1} = u^{i + 1} / (ω^{i + 1} - ω^{i}),$ 将 $T (ω^{*})$ 在 $ω^{i}$ 处泰勒展开，得到 $0 = T (ω^{i}) = T (ω^{i}) +$ $(ω^{*} - ω^{i}) T_{ω} (ω^{i}) + \dots,$ 忽略高阶部分，并用 $ω^{i + 1}$ 代替 $ω^{*}$ 得到：

T (ω^{i}) + (ω^{i + 1} - ω^{i}) T_{ω} (ω^{i}) = 0

。（20）

式（20）两端同乘以 $u^{i + 1}$ ，得：

ω^{i + 1} = ω^{i} - \frac{(u^{i + 1})^{H} T (ω^{i}) u^{i + 1}}{(u^{i + 1})^{H} T_{ω} (ω^{i}) u^{i + 1}}

。（21）

式（19）、式（21）是算法的基础，下面给出牛顿法求解式（17）的算法1流程图，如图3所示。

图 3. 算法1流程图

Fig. 3. Algorithm 1 flow chart

下载图片查看所有图片

需要注意到，式（21）中第2项的分母不为0。算法在满足非退化条件式（22）时具有二次收敛性，证明过程可见文献［11］。

{u^{*}}^{H} T_{ω} (ω) u^{*} \neq 0

。（22）

使用基于有限元方法的COMSOL软件与牛顿迭代相结合来求解晶体的能带结构，进一步描述COMSOL与牛顿法结合共同计算的实现过程：

1）建立几何模型。进入COMSOL软件选取空间维度为2D，然后选择电磁波、频域，随后进入特征频率模块。添加表示周期性单元几何结构的参数，绘制原胞区域。

2）添加表示倒易晶格基矢和波矢 $k$ 的变量。

3）设置材料参数。在COMSOL中将色散关系定义为解析函数，赋给相应的原胞区域。

4）设置周期性边界条件并选择计算模式（TE/TM模式）。

5）进行网格剖分。为满足周期性边界条件使对边上的网格相同。

6）COMSOL初始特征频率分析。选取1.55 $μ m$ 作为工作波长，该波长对应的光波频率作为工作频率（约193 THz，记作freq），继而得到该频率处的介电常数。利用COMSOL在193 THz附近查找零波矢 $k = 0$ 时的所有特征频率，将接下来迭代求解用到的频率（即freq1）的初始值取自本次初始特征频率分析得到的所有特征值。

7）求解及数据后处理。将上一步初始特征频率分析得到的解导出，作为牛顿迭代的初始值，随后利用算法1完成后续求解，画出光子晶体的能带分布图。

通常情况下，研究光子晶体的能带情况需要选择某一中心频率，在此频率附近探究光子带隙和能带分布。上述步骤6）使用COMSOL软件确定初始特征频率分析得到的解作为初始值，物理上可以这样理解：为研究光子晶体的能带结构提供多个实验所用的估算值，在此基础上经过后续的算法研究和验证，围绕多个初始频率作为中心频率，进一步计算真实的能带分布情况。

5　数值算例

5.1　算例1

考虑空气柱在色散材料背景中周期排列构成三角晶格色散光子晶体。利用Drude模型来表示背景材料的介电函数：

ε (ω) = ε_{\infty} [1 - (ω_{p}^{2} / ω^{2})],

（23）

式中： $ω_{p}$ 为材料的等离子体频率； $ε_{\infty}$ 为无限频率的介电常数。以 $ε_{\infty} = 1, ω_{p} = (2 π c) / a, r = 0.48 a$ （a为三角晶格的晶格常数，r为晶体原胞的半径）的三角形色散光子晶体为例，分别求解其TE、TM模式下的能带结构。给出此材料的色散关系曲线，如图4所示。

图 4. 算例1色散曲线

Fig. 4. Dispersion relation curve of example 1

下载图片查看所有图片

TM、TE模式下的能带结构分别如图5、图6所示。

图 5. TM模式下的算例1能带图

Fig. 5. Computed band structure for example 1:TM mode

下载图片查看所有图片

图 6. TE模式下的算例1能带图

Fig. 6. Computed band structure for example 1: TE mode

下载图片查看所有图片

图5为三角形晶格色散光子晶体在TM模式下的能带图，可以看出存在一条带隙，光子带隙频段范围为 $0.6826 ω_{p}$ ~ $0.7323 ω_{p}$ 。计算结果与文献［12-13］中结果一致。图6为TE模式下的能带结构，可以看出当频率大于 $0.5 ω_{p}$ ，开始出现带隙。图6中填充部分为2条带隙：第一条带隙频率范围为 $[0.5164 ω_{p}, 0.5369 ω_{p}]$ ，带宽 $0.0205 ω_{p}$ ；第二条带隙频率范围 $[0.5546 ω_{p},$ $0.5706 ω_{p}]$ ，带宽为 $0.0160 ω_{p}$ 。

在其他参数保持不变的前提下，调整晶柱半径r，晶格填充比 $f (f = r / a)$ 从0到0.5，画出TM模式下光子晶体带隙变化图，结果如图7所示。通过观察可以明确该模式下色散光子晶体的带隙随填充比的变化情况。

图 7. TM模式下的算例1带隙

Fig. 7. Gap map of example 1 in TM mode

下载图片查看所有图片

5.2　算例2

研究正方晶格有耗色散光子晶体的能带结构，选择中心晶柱为色散材料，基体部分是空气的色散光子晶体。其中，色散材料介电函数：

ε (ω) = 1 - \frac{ω_{p}^{2}}{ω^{2} - i ω ω_{τ}},

（24）

式中： $ω_{p}$ 为等离子频率； $ω_{τ}$ 为衰减频率，当 $ω_{τ} \neq 0$ 时称为有耗色散光子晶体。 $ω_{p} = 2 π \times 1914 T H z$ ， $ω_{τ} = 2 π \times 8.34 T H z$ ，填充比 $f = 0.2$ 。此材料的色散关系曲线如图8所示。

图 8. 算例2色散曲线

Fig. 8. Dispersion relation curve of example 2

下载图片查看所有图片

为进一步说明本算法的创新性，以此例为研究对象，进一步探讨TM模式下本算法与未使用本算法（仅选用频率 $2 \times 10^{15}$ Hz时对应的介电常数，此时 $ε = 0.09$ ，如图8所示）所得的光子能带，并将两种方法的计算结果进行对比。

可以看出两种算法求得的光子能带分布差异较大，与文献［7］中结果进行对比，可验证本算法所得是正确的，如图9、图10所示。仅使用一个频率处的介电常数计算得到的光子能带与真实情况有很大差距。

图 9. 本算法结果图

Fig. 9. Results chart obtained by algorithm in this paper

下载图片查看所有图片

图 10. 另一种方法结果图

Fig. 10. Results chart obtained by relative algorithm

下载图片查看所有图片

5.3　算例3

求解无损有耗介电材料在空气中均匀分布所构成的正方晶格光子晶体的能带结构分布情况。材料的介电函数：

ε (ω) = 1 - \frac{ω_{p}^{2}}{ω (ω + i γ)},

（25）

式中： $γ$ 为反电子弛豫时间， $ω_{p} = 2 π c / a$ ， $γ = 0.01 ω_{p}$ 。此材料的色散关系曲线如图11所示。

图 11. 算例3色散曲线

Fig. 11. Dispersion relation curve of example 3

下载图片查看所有图片

由文献［14］可知归一化频率 $(ω a) / (2 π c)$ 在 $[0,2]$ 范围内，TE模式下的光子晶体是没有带隙的。以TM模式为例，分别计算填充比 $f = 0.0564 、 f = 0.4720$ 时的能带结构分布，TM模式下两种不同填充比的光子晶体能带结构如图12、图13所示。

图 12. f=0.0564时TM模式下的算例3能带图

Fig. 12. Computed band structure for example 3: TM mode, f=0.0546

下载图片查看所有图片

图 13. f=0.4720时TM模式下的算例3能带图

Fig. 13. Computed band structure for example 3: TM mode, f=0.4720

下载图片查看所有图片

文献［14］中填充率 $(π r^{2}) / a^{2} = 0.01$ 对应本文 $f = 0.0564$ ， $(π r^{2}) / a^{2} = 0.7$ 对应 $f = 0.4720$ ，计算结果与文献［14］中方法所得结果相一致。

6　结论

基于有限元方法和牛顿迭代法研究了二维色散介质光子晶体的带隙结构，根据给定的数值算法，计算了3个不同的二维色散介质光子晶体的能带分布。首先，利用有限元方法将特征值问题离散为矩阵形式；然后，选择合适的初值基于牛顿法迭代求解，从而得到色散关系，获得光子晶体的能带结构，分析TE、TM两种模式下的色散介质光子晶体的能带结构。此外，通过调整光子晶体填充比，求得了光子晶体带隙变化图。有限元离散与牛顿迭代相结合方法为二维色散介质能带结构分布问题提供了有效的求解途径，具有一定的理论意义和实用价值。

参考文献

[1] 陈剑锋, 梁文耀, 李志远. 磁光光子晶体中拓扑光子态研究进展[J]. 光学学报, 2021, 41(8): 0823015.

Chen J F, Liang W Y, Li Z Y. Progress of topological photonic state in magneto-optical photonic crystal[J]. Acta Optica Sinica, 2021, 41(8): 0823015.

[2] 马锡英. 光子晶体原理及应用[M]. 北京: 科学出版社, 2010: 6-12.

MaX Y. Principle and application of photonic crystal[M]. Beijing: Science Press, 2010: 6-12.

[3] 李春早. 色散介质光子晶体电磁特性的研究[D]. 南京: 南京航空航天大学, 2012: 2.

LiC Z. Study on electromagnetic characteristics of dielectric photonic crystals[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2012: 2.

[4] Jiang B, Zhou W J, Chen W, et al. Improved plane-wave expansion method for band structure calculation of metal photonic crystal[J]. Chinese Physics Letters, 2011, 28(3): 034209.

[5] Fietz C, Urzhumov Y, Shvets G. Complex k band diagrams of 3D metamaterial/photonic crystals[J]. Optics Express, 2011, 19(20): 19027-19041.

[6] 魏源, 肖峰, 吴博, 等. 基于FDFD方法周期结构的薄膜太阳能电池特性[J]. 光子学报, 2014, 43(1): 131001.

Wei Y, Xiao F, Wu B, et al. Thin-film solar cell’s characteristic with periodic structure based on FDFD method[J]. Acta Photonica Sinica, 2014, 43(1): 131001.

[7] Xiao W Q, Sun J G. Band structure calculation of photonic crystals with frequency-dependent permittivities[J]. Journal of the Optical Society of America A, 2021, 38(5): 628-633.

[8] Yuan J H, Lu Y Y. Photonic bandgap calculations with Dirichlet-to-Neumann maps[J]. Journal of the Optical Society of America A, 2006, 23(12): 3217-3222.

[9] Yuan J H, Lu Y Y. Computing photonic band structures by Dirichlet-to-Neumann maps: the triangular lattice[J]. Optics Communications, 2007, 273(1): 114-120.

[10] SpenceA, PoultonC G. Inverse iteration for nonlinear eigenvalue problems in electromagnetic scattering[M]∥Movchan A B. IUTAM symposium on asymptotics, singularities and homogenisation in problems of mechanics. Solid mechanics and its applications. Dordrecht: Kluwer Academic Publishers, 2006, 113: 585-594.

[11] Spence A, Poulton C. Photonic band structure calculations using nonlinear eigenvalue techniques[J]. Journal of Computational Physics, 2005, 204(1): 65-81.

[12] 牛凯坤, 王丽华, 黄志祥, 等. 三角晶格有耗色散光子晶体的能带结构分析[J]. 光子学报, 2016, 45(3): 0319002.

Niu K K, Wang L H, Huang Z X, et al. Band structure of triangular lattices photonic crystals with lossy and dispersive materials[J]. Acta Photonica Sinica, 2016, 45(3): 0319002.

[13] 王瀚林. 有耗色散光子晶体能带结构研究[D]. 合肥: 安徽大学, 2015: 30-31.

WangH L. Research on band structure of photonic crystal with lossy and dispersive materials[D]. Hefei: Anhui University, 2015: 30-31.

[14] Xiao W Q, Gong B, Sun J G, et al. Finite element calculation of photonic band structures for frequency dependent materials[J]. Journal of Scientific Computing, 2021, 87(1): 27.

钟相辉, 袁健华. 二维色散介质光子晶体的能带结构求解问题研究[J]. 激光与光电子学进展, 2023, 60(9): 0926001. Xianghui Zhong, Jianhua Yuan. Numerical Analysis of the Band Structure of Two-Dimensional Dispersive Dielectric Photonic Crystals[J]. Laser & Optoelectronics Progress, 2023, 60(9): 0926001.