An $hp$-Adaptive Sampling Algorithm for Dispersion Relation Reconstruction of 3D Photonic Crystals

In SI convention, the time harmonic Maxwell equations with $\mathbf{E}(\mathbf{x},t):=\mathbf{E}(\mathbf{x})\mathrm{e}^{-i\omega t}$ and $\mathbf{H}(\mathbf{x},t):=\mathbf{H}(\mathbf{x})\mathrm{e}^{-i\omega t}$ for linear, non-dispersive, and nonmagnetic media, with free charges and free currents, consist of a system of four equations [25]

where $\mathbf{E}$ is the electric field, $\mathbf{H}$ is the magnetic field, the scalar $\omega\geq 0$ is the frequency of the electromagnetic wave, $\mu_{0}$ is the vacuum permeability, $\epsilon_{0}$ is the vacuum permittivity, and $\epsilon\in L^{\infty}(\mathbb{R}^{3};\mathbb{R}^{+})$ is the relative permittivity.

Applying the curl operator to (2.1b) and using (2.1a) give the equations for the magnetic field

with $\epsilon_{0}\mu_{0}=c^{-2}$, where $c$ is the speed of light. The electric field $\mathbf{E}(\mathbf{x})$ can be derived from (2.1b):

which implies (2.1c). The above two equations form the governing system of 3D PhCs.

Next, 3D periodic PhCs possess a discrete translational symmetry in $\mathbb{R}^{3}$ [26], i.e., the relative permittivity $\epsilon(\mathbf{x})$ satisfies

The primitive lattice vectors, denoted as $\mathbf{a}_{1},\mathbf{a}_{2},\mathbf{a}_{3}$, are the shortest possible vectors that satisfy this condition which span the unit cell $\Omega$. For example, square lattice has primitive lattice vectors $\mathbf{a}_{i}=a{e}_{i}$ for $i=1,2,3$, with $a\in\mathbb{R}$ being the lattice constant and $(e_{i})_{i=1,2,3}$ being the canonical basis in $\mathbb{R}^{3}$. Additionally, the reciprocal lattice vectors, $\mathbf{b}_{i}$, are defined by the property

which generate the so-called reciprocal lattice. The elementary cell of the reciprocal lattice is the (first) Brillouin zone $\mathcal{B}$, i.e., the region closer to a certain lattice point than to any other lattice points in the reciprocal lattice. Note that if the materials in the unit cell have additional symmetry, e.g. mirror symmetry, we can further restrict $\mathbf{k}$ to the irreducible Brillouin zone (IBZ).

Based on Bloch’s theorem [28], the spectrum corresponding to a periodic structure in the full space is equivalent to the union of all spectra corresponding to the quasi-periodic problems in the unit cell. Thus, we apply Bloch’s theorem and aim to find the wave function $\mathbf{H}$ taking the form of a plane wave modulated by a periodic function. Mathematically, it can be written as $\mathbf{H}(\mathbf{x})=e^{i\mathbf{k}\cdot\mathbf{x}}\mathbf{u}(\mathbf{x})$, where $\mathbf{u}(\mathbf{x})$ is a periodic function satisfying $\mathbf{u}(\mathbf{x}+\mathbf{a_{i}})=\mathbf{u}(\mathbf{x})$ for $i=1,2,3$ and $\mathbf{k}$ is the wave vector varying in the Brillouin zone $\mathcal{B}$. Thus, by Bloch’s theorem, the eigenvalue problem (2.2) reduces to

where $\lambda=\left(\omega c^{-1}\right)^{2}$, $\mathbf{k}$ varies in the Brillouin zone $\mathcal{B}$, and $\mathbf{u}(\mathbf{x})$ satisfies the periodic boundary conditions $\mathbf{u}(\mathbf{x})=\mathbf{u}(\mathbf{x}+\mathbf{a_{j}})$ with $\mathbf{a_{j}}$ being the primitive lattice vector for $j=1,2,3$. For all $\mathbf{k}\in\mathcal{B}$, we can enumerate these eigenvalues in a non-decreasing manner (with their multiplicities counted) as

Then $\{\lambda_{n}(\mathbf{k})\}_{n=1}^{\infty}$ is an infinite sequence with $\lambda_{n}(\mathbf{k})$ being a continuous function with respect to the wave vector $\mathbf{k}$ and $\lambda_{n}(\mathbf{k})\to\infty$ when $n\to\infty$ [21].

Below we state the regularity of band functions, which can be proved in a similar manner as the 2D case [51]. Here, $\mathrm{Lip}(\mathcal{B})$ is the space of Lipschitz continuous functions in the Brillouin zone $\mathcal{B}$ and $\mathring{A}(\mathcal{B})$ denotes the space of piecewise analytic functions with a zero Lebesgue measure of singular point set composed of points of degeneracy (i.e., the set of points at which certain band functions intersect, with zero Lebesgue measure). Throughout, the set of singular points is denoted by $\mathbf{S}$.

Next, we introduce the weak formulation to (2.5). We first introduce useful notation for function spaces. For an open, bounded and simply-connected Lipschitz domain $\Omega\subset\mathbb{R}^{3}$, we define the space of square integrable functions $L^{2}(\Omega)$ equipped with the inner product $(u,v)=\int_{\Omega}u(\mathbf{x})\cdot\bar{v}(\mathbf{x})\,\mathrm{d}\mathbf{x}$ and the corresponding norm by $\|u\|=\left(\int_{\Omega}|u(\mathbf{x})|^{2}\,\mathrm{d}\mathbf{x}\right)^{\frac{1}{2}}$, where $\mathbf{\bar{\cdot}}$ denotes the complex conjugate. We denote by $H^{1}(\Omega)$ the usual Sobolev space with norm

We also write $(\cdot,\cdot)$, $\|\cdot\|$ and $(\cdot,\cdot)_{\mathbf{H}^{1}}$, $\|\cdot\|_{\mathbf{H}^{1}}$ for the inner products and norms of vector-valued functions spaces $L^{2}(\Omega)^{3}$, $H^{1}(\Omega)^{3}$, respectively. Let $(\mathbf{u},\mathbf{v})_{\mathbf{H}(\text{curl};\Omega)}=(\mathbf{u},\mathbf{v})+(\nabla\times\mathbf{u},\nabla\times\mathbf{v})$ be the standard inner product on the Sobolev space $\mathbf{H}(\text{curl};\Omega):=\{\mathbf{v}\in L^{2}(\Omega)^{3}:\nabla\times\mathbf{v}\in L^{2}(\Omega)^{3}\}$ and $(\mathbf{u},\mathbf{v})_{\mathbf{H}(\text{div};\Omega)}=(\mathbf{u},\overline{\mathbf{v}})+(\nabla\cdot\mathbf{u},\nabla\cdot\mathbf{v})$ be the inner product on $\mathbf{H}(\text{div};\Omega):=\{\mathbf{v}\in L^{2}(\Omega)^{3}:\nabla\cdot\mathbf{v}\in L^{2}(\Omega)^{3}\}$. Also we will use the following spaces for the weak formulation:

Here $\mathbf{n}$ is the unit outward normal vector to $\partial\Omega$. The periodic Sobolev spaces $H_{\text{per}}^{1}(\Omega)$, $\mathbf{H}_{\text{per}}(\text{div};\Omega)$, $\mathbf{H}_{\text{per}}(\text{curl};\Omega)$ are the closures of the spaces $C^{\infty}_{\text{per}}({\Omega})$, $C^{\infty}_{\text{div}}({\Omega})$, $C^{\infty}_{\text{curl}}({\Omega})$ with respect to the Sobolev space norms $\|\cdot\|_{H^{1}}$,$\|\cdot\|_{\mathbf{H}(\text{div};\Omega)}$, $\|\cdot\|_{\mathbf{H}(\text{curl};\Omega)}$ induced by the inner product mentioned above respectively. The periodic space $L^{2}_{\text{per}}(\Omega)$ is the closure of the space $C^{\infty}_{\text{per}}({\Omega})$ with respect to the norm $\|\cdot\|$.

Then the spaces for the weak formulation are given by

For $\mathbf{u},\mathbf{v}\in\mathbf{V}$ and $q\in Q$, we define the sesquilinear forms

Then the variational formulation to (2.5) can be written in a mixed form: for a given $\mathbf{k}\in\mathcal{B}$, find non-trivial eigentriplets $\{(\lambda_{n},\mathbf{u}_{n},p_{n})\}_{n=1}^{\infty}\subset(\mathbb{R},\mathbf{V},Q)$ such that

Finally, we introduce a regularity assumption for the finite element analysis [4, Assumption $\operatorname{H}_{0}$],

In this section we describe the discretization of problem (2.7) for each wave vector $\mathbf{k}\in\mathcal{B}$ using the $\mathbf{H}(\text{curl})$- conforming $\mathbf{k}$-modified Nédélec edge element method, which is stable with respect to $\mathbf{k}$. Since the divergence of the curl of any vector is zero, taking the divergence on both sides of the first equation in (2.5) gives its second equation. Furthermore, the zero eigenvalue of (2.5) has an infinite multiplicity since the curl operator has an infinite-dimensional kernel. Thus, a discretization of the eigenvalue problem (2.5) solely using the scalar node-based FEM may contain non-physical solutions (i.e., spurious modes [48]). Indeed, most eigenvalue solvers are inefficient when dealing with such eigenvalue problems [2].

To eliminate these spurious modes, the standard approach is to use the edge FEM and add the constraint (2.5) explicitly [31, 48]. The Lagrange multiplier in the mixed form in (2.7) is one common approach to add the divergence constraint to the system. Also, several methods have been proposed to eliminate this large null space [14, 11, 24, 38, 36, 35], including the special $\mathbf{H}(\text{curl})$-conforming $\mathbf{k}$-modified Nédélec edge elements [15, 4] and the mixed discontinuous Galerkin (DG) formulation with modified Nédélec basis functions [35]. The $\mathbf{H}(\text{curl})$- conforming $\mathbf{k}$-modified Nédélec edge element method was first proposed and analyzed in [14, 15] under strong regularity restrictions and in the case of uniform mesh sequences, with an improved analysis under minimal assumptions on the regularity of the eigensolutions and on the mesh sequences [4]. We recap the Nédélec edge element method below.

Let $\mathcal{T}_{h}$ be a periodic, shape-regular, conformal tetrahedral mesh for $\Omega$ aligned with the possible discontinuities of $\alpha,\beta$, where the phrase “periodic” means that the meshes are the same on each pair of the boundary. Let $h=\max_{K\in\mathcal{T}_{h}}h_{K}$, where $h_{K}$ is the diameter of $K$. We denote the set of all edges by $\mathcal{E}$ and the set of all vertexes by $\mathcal{V}$. The conforming finite element spaces of (2.6) are [4]

where $P_{k}(K)$ is the space of local polynomials of total degree less than or equal to $k$ on $K$, and $S_{k}(K)$ is the space of vector polynomials with the form $\mathbf{a}(\mathbf{x})+\mathbf{b}(\mathbf{x})\times\mathbf{x}$ with $\mathbf{a},\mathbf{b}\in P_{k}(K)^{3}$.

Given $\mathbf{k}\in\mathcal{B}$, the discretization of problem (2.7) reads: find a non-trivial eigentriplet $(\lambda_{n,h},\mathbf{u}_{n,h},p_{n,h})\in(\mathbb{R},\mathbf{V}_{h}^{\mathbf{k}},Q_{h}^{\mathbf{k}})$ with $(\mathbf{u}_{n,h},\lambda_{n,h})\neq(\mathbf{0},0)$, such that

For $\mathbf{k}=\mathbf{0}$, one should make the spaces $\mathbf{V}$, $Q$ and their corresponding FEM spaces $\mathbf{V}^{\mathbf{0}}_{h}$, $Q^{\mathbf{0}}_{h}$ orthogonal to constant vectors or constants in order to guarantee the stability of (2.7) and (2.9) [5].

Next, we derive the formula of the gradient $\partial_{j}\lambda_{n}:=\frac{\partial\lambda_{n}}{\partial k_{j}}$ for $n\in\mathbb{N}_{+}$ and $j=1,2,3$, and provide its numerical approximation. Recall that $(e_{j})_{j=1,2,3}$ denotes the canonical basis in $\mathbb{R}^{3}$.

The use of modified basis functions (2.8) greatly simplifies the computation of $\partial_{j}\lambda_{n}$, and hence reduce the computational complexity. For example, for the lowest order Nédélec FEM space with $k=0$ in (2.8), the degrees of freedom for $\mathbf{V}_{h}^{\mathbf{k}}$ are

Here, $\mathbf{z}_{e}:=(\mathbf{x}_{e}+\mathbf{y}_{e})/2$ denotes the midpoint of edge $e\in\mathcal{E}$, identified with the ordered vertices $(\mathbf{x}_{e},\mathbf{y}_{e})$, and $\mathbf{t}_{e}:=(\mathbf{y}_{e}-\mathbf{x}_{e})/\|\mathbf{y}_{e}-\mathbf{x}_{e}\|$ represents the unit tangent vector of $e$. Hence, $\mathbf{V}_{h}^{\mathbf{k}}$ of the lowest order can be expressed as

with $\ell_{e,\mathbf{k}}(\psi_{e,\mathbf{k}})=1$ and $\ell_{e,\mathbf{0}}(\psi_{e,\mathbf{0}})=1$. Similarly, $Q_{h}^{\mathbf{k}}$ of the lowest-order can be represented by

with $\ell_{v,\mathbf{k}}(q_{v,\mathbf{k}})=1$ and $\ell_{v,\mathbf{0}}(q_{v,\mathbf{0}})=1$, which has the vertex-based degrees of freedom defined as $\ell_{v,\mathbf{k}}(q):=e^{i\mathbf{k}\cdot(x-\mathbf{x}_{v})}q(\mathbf{x}_{v})$, for all $v\in\mathcal{V}$, $\mathbf{x}_{v}$ is the coordinate of a vertex in $\mathcal{V}$.

With the basis functions defined above, we obtain

Together with (2.11), we obtain

Thus, the use of the modified basis functions separates the integral with the parameter $\mathbf{k}$, transforming it into a standalone exponential function during the computation. This separation streamlines evaluating the integrals, making it a conventional integration involving parameter-free polynomials.

We first recall the definition of a regular mesh without hanging nodes [45]. Let $\mathcal{B}\subset\mathbb{R}^{3}$ be a polygon with straight sides and a mesh $\mathcal{T}=\{T_{i}\}_{i=1}^{M}$ be a partition of $\mathcal{B}$ into $M$ open, disjoint tetrahedrons $T_{i}$ with $\overline{D}=\bigcup_{i=1}^{M}\overline{T_{i}}$. For each $T\in\mathcal{T}$, we denote by $h_{T}$ its diameter and $\rho_{T}$ the diameter of its largest inscribed ball.

Now we develop a strategy to generate a family of nested regular tetrahedral meshes $\mathcal{T}_{n}$, $n=1,2,\cdots$, such that $\mathcal{T}_{n+1}$ is a refinement of $\mathcal{T}_{n}$. To this end, we introduce several notations for the given regular tetrahedral mesh $\mathcal{T}_{n}$. The mesh size of $\mathcal{T}_{n}$ is defined as $h_{n}:=\max_{T\in\mathcal{T}_{n}}h_{T}$. The collection of all the faces, edges, and vertices over $\mathcal{T}_{n}$ is denoted as $\mathcal{F}_{n}$, $\mathcal{E}_{n}$ and $\mathcal{V}_{n}$, respectively. For any $T\in\mathcal{T}_{n}$, we define its faces as $\mathcal{F}_{T}:=\partial T\bigcap\mathcal{F}_{n}$, and the $j$th face is denoted by $F_{T}^{j}$ for $j=1,2,3,4$. Similarly, we define its edges and vertices as $\mathcal{E}_{T}:=\partial T\bigcap\mathcal{E}_{n}$ and $\mathcal{V}_{T}:=\partial T\bigcap\mathcal{V}_{n}$, and the $j$th edge and vertex are denoted by $E_{T}^{j}$ for $j=1,\cdots,6$, and $V_{T}^{j}$ for $j=1,2,3,4$. For any $F\in\mathcal{F}_{n}$, we define the collection of elements having $F$ as a face as $\mathcal{T}_{F}:=\bigcup\{T\in\mathcal{T}_{n}:F\in\mathcal{F}_{T}\}$ and the collection of elements having $E$ as an edge as $\mathcal{T}_{E}:=\bigcup\{T\in\mathcal{T}_{n}:E\in\mathcal{E}_{T}\}$. Finally, we define $P_{m}(T)$ as the space of polynomial functions with total degrees at most $m\in\mathbb{N}^{+}$ on element $T$, and $P_{m}(F)$, $P_{m}(E)$ as the space of polynomial functions of degrees at most $m\in\mathbb{N}^{+}$ on each face $F\in\mathcal{F}_{n}$ and each edge $E\in\mathcal{E}_{n}$, respectively.

Theorem 2.1 implies that each $\lambda_{n}(\mathbf{k})$ is a piecewise analytic function of the wave vector $\mathbf{k}$, with singular points at the branch points. Our mesh design strategy aims to identify elements within $\mathcal{T}_{n}$ where $\lambda_{\ell}$ and $\lambda_{\ell+1}$ are not smooth. To accomplish this, we introduce an element-wise indicator

which uses the adjacent values on the vertices of each tetrahedron element to ascertain the regularity of the eigenvalues as functions of $\mathbf{k}$ within each element. Furthermore, the element-wise tolerance is

where $\kappa\geq 2\sqrt{3}$ is a constant and

is an upper bound of the gradients of $\lambda_{q}(\mathbf{k})$ with $\max\{1,\ell-1\}\leq q\leq\ell+1$.

Give a mesh $\mathcal{T}_{n}$, the collection of marked elements is given by

Then we apply the iterated longest edge bisection algorithm (LEB) [42] to divide each marked element $T\in\mathcal{M}_{n}$ into two children elements, which typically follows two steps. First, creating a new vertex at the midpoint of its longest edge and then dividing this tetrahedron by the plane that passes through this new vertex and the two vertices opposite the longest edge. Second, performing additional bisection to avoid any hanging nodes by means of recursive algorithm, which was proven to be finite [42].

We summarize the proposed mesh refinement procedure in Algorithm 1, where the parameters $\operatorname{tol}_{2}$ and $n_{\max}^{1}$ denote the smallest allowed element size and the maximum number of iterations, respectively.

Next, we show that all elements containing singular points are marked for refinement. For the mesh $\mathcal{T}_{n}$ at the $n$th iteration, let

be a partition of the mesh $\mathcal{T}_{n}$ based upon the singular point set $\mathbf{S}$.

In practice, we use a practical tolerance given by

as a substitute for $C_{\max,\ell}(T)$ which only requires $\partial_{j}\lambda_{q}$ for $i=1,2,3$ and $\max\{1,\ell-1\}\leq q\leq\ell+1$ over $\mathcal{V}_{T}$, which can be obtained from (2.11).

Let $\{\mathcal{T}_{n}\}_{n=1}^{n_{\operatorname{iter}}}$ be a sequence of nested regular tetrahedral meshes generated by Algorithm 1. First, we establish a finite element space of continuous element-wise polynomial functions $V_{n}$ as the ansatz space to approximate the target band functions. Then, we define element-wise interpolation operator $\Pi_{\mathcal{T}_{n}}:C(\mathcal{B})\to V_{n}$ to reconstruct the target band functions on conforming finite element spaces. First, we introduce the concept of a layer of each element $T\in\mathcal{T}_{n}$.