Numerical Electromagnetic Frequency Domain Analysis with Discrete Exterior Calculus

TL;DR

This paper presents a novel frequency domain electromagnetic analysis method using discrete exterior calculus on arbitrary meshes, effectively handling complex boundary conditions and material properties.

Contribution

It introduces a DEC-based formulation for electromagnetic problems on unstructured meshes, incorporating effective signed dual volumes and natural charge continuity.

Findings

Accurate results for waveguides, fibers, and photonic crystals.

Effective handling of various boundary conditions.

Validation through multiple numerical examples.

Abstract

In this paper, we perform a numerical analysis in frequency domain for various electromagnetic problems based on discrete exterior calculus (DEC) with an arbitrary 2-D triangular or 3-D tetrahedral mesh. We formulate the governing equations in terms of DEC for 3D and 2D inhomogeneous structures, and also show that the charge continuity relation is naturally satisfied. Then we introduce effective signed dual volume to incorporate material information into Hodge star operators and take into account the case when circumcenters fall outside triangles or tetrahedrons, which may lead to negative dual volume without Delaunay triangulation. Then we demonstrate the implementation of various boundary conditions, including perfect magnetic conductor (PMC), perfect electric conductor (PEC), Dirichlet, periodic, and absorbing boundary conditions (ABC) within this method. An excellent agreement is achieved through the numerical calculation of several problems, including homogeneous waveguides, microstructured fibers, photonic crystals, scattering by a 2-D PEC, and resonant cavities.

Tables

Table 1. Table 1 : Fundamental resonant frequencies (GHz) of inhomogeneous cavities.

Model	This work	[30]	[29]	Error (%)
LES-mode resonator ($s=2$ mm)	15.628	15.65	$-$	0.14
LES-mode resonator ($s=4$ mm)	13.346	13.35	13.34	0.03
Hybrid mode resonator	8.431	8.43	8,42	0.01
Dielectric sphere*	152.48	$-$	153.3	0.47

Equations

$$dE=i\omega B,dH=-i\omega D+J,dB=0,dD=\rho ,dJ=i\omega \rho .$$

$${\int }_{D}d\omega ={\oint }_{\partial D}\omega$$

$$i\omega {\int }_{A_i}B={\int }_{A_i}dE={\oint }_{\partial A_i}E=\sum _{l_j\in \partial A_i}{d}_{i,j}^{(1)}{\int }_{l_j}E$$

$$E_j\triangleq {\int }_{l_j}E={\int }_{l_j}\left[\mathbf{E}(\mathbf{r})\cdot {\hat{l}}_j\right]dl,B_i\triangleq {\int }_{A_i}B={\int }_{A_i}\left[\mathbf{B}(\mathbf{r})\cdot {\hat{n}}_i\right]dA$$

$$i\omega \mathbf{B}={\overline{\mathbf{d}}}^{(1)}\mathbf{E}$$

$${\left[{\overline{\mathbf{d}}}^{(1)}\right]}_{i,j}=d_{i,j}^{(1)}=\{\begin{array}{ll}\pm 1& \text{if lj is an edge of Ai,}\\ 0& \text{otherwise.}\end{array}$$

$${\left[{\overline{\mathbf{d}}}^{(2)}\right]}_{i,j}$$

$${\left[{\overline{\mathbf{d}}}^{(0)}\right]}_{i,j}$$

$${\int }_{V_i}dB={\oint }_{\partial V_i}B=\sum _{A_j\in \partial V_i}{d}_{i,j}^{(2)}{\int }_{A_j}B=0$$

$${\overline{\mathbf{d}}}^{(2)}\mathbf{B}=0.$$

$$\mathbf{D}$$

$$\mathbf{J}$$

$$D_i$$

$$J_i$$

$${\overline{\mathbf{d}}}_{\text{dual}}^{(n-k)}=(-1{)}^k({\overline{\mathbf{d}}}^{(k-1)}{)}^T$$

$${\overline{\mathbf{d}}}_{\text{dual}}^{(1)}=({\overline{\mathbf{d}}}^{(1)}{)}^T,{\overline{\mathbf{d}}}_{\text{dual}}^{(2)}=-({\overline{\mathbf{d}}}^{(0)}{)}^T,{\overline{\mathbf{d}}}_{\text{dual}}^{(0)}=({\overline{\mathbf{d}}}^{(2)}{)}^T.$$

$${\overline{\mathbf{d}}}_{\text{dual}}^{(1)}\mathbf{H}=({\overline{\mathbf{d}}}^{(1)}{)}^T\mathbf{H}=-i\omega \mathbf{D}+\mathbf{J},$$

$${\overline{\mathbf{d}}}_{\text{dual}}^{(2)}\mathbf{D}=-({\overline{\mathbf{d}}}^{(0)}{)}^T\mathbf{D}=\mathbf{\rho }.$$

$${\overline{\mathbf{d}}}_{\text{dual}}^{(2)}\mathbf{J}=-({\overline{\mathbf{d}}}^{(0)}{)}^T\mathbf{J}=i\omega \mathbf{\rho },$$

$$\mathbf{D}={ϵ}_0{\star }_{ϵ}^{(1)}\mathbf{E},\mathbf{H}={\mu }_0^{-1}{\star }_{{\mu }^{-1}}^{(2)}\mathbf{B}.$$

$$\left[({\overline{\mathbf{d}}}^{(1)}{)}^T{\star }_{{\mu }^{-1}}^{(2)}{\overline{\mathbf{d}}}^{(1)}\right]\mathbf{E}=k_0^2\left[{\star }_{ϵ}^{(1)}\right]\mathbf{E}+i\omega \mathbf{J},$$

$${\mathbf{e}}^T\left[({\overline{\mathbf{d}}}^{(1)}{)}^T{\star }_{{\mu }^{-1}}^{(2)}{\overline{\mathbf{d}}}^{(1)}\right]\mathbf{e}={\mathbf{b}}^T{\star }_{{\mu }^{-1}}^{(2)}\mathbf{b}\ge 0$$

$$B_{z,i}$$

$$D_{z,i}$$

$$E_{s,i}$$

$$B_{s,i}$$

$$D_{s,i}$$

$$H_{s,i}$$

$${\overline{\mathbf{d}}}_{\text{dual}}^{(0)}=({\overline{\mathbf{d}}}^{(1)}{)}^T,\text{and}{\overline{\mathbf{d}}}_{\text{dual}}^{(1)}=-({\overline{\mathbf{d}}}^{(0)}{)}^T.$$

$$\overline{\mathbf{ϵ}}=\left[\begin{array}{cc}{\overline{\mathbf{ϵ}}}_s& 0\\ 0& {ϵ}_{zz}\end{array}\right],\overline{\mathbf{\mu }}=\left[\begin{array}{cc}{\overline{\mathbf{\mu }}}_s& 0\\ 0& {\mu }_{zz}\end{array}\right],$$

$${\mathbf{D}}_z$$

$${\mathbf{H}}_z$$

$${\nabla }_s\cdot {\nabla }_s{E}_z(\mathbf{r})+k_s^2{E}_z(\mathbf{r})$$

$${\nabla }_s\cdot {\nabla }_s{H}_z(\mathbf{r})+k_s^2{H}_z(\mathbf{r})$$

$$\left[({\overline{\mathbf{d}}}^{(0)}{)}^T{\star }^{(1)}{\overline{\mathbf{d}}}^{(0)}\right]{\mathbf{E}}_z-k_s^2\left[{\star }^{(0)}\right]{\mathbf{E}}_z=0.$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Finite difference time domain method (FDTD) and finite element methods (FEM) have been widely applied to solve electromagnetic problems [1]. And they are both based on the vectorial version of Maxwell’s equations. As is well known, differential forms can be used to recast Maxwell’s equations in a more succinct fashion, which completely separate metric-free and metric-dependent parts, see [2, 3]. Maxwell’s equations and charge continuity condition in terms of differential forms within frequency domain are written as:

[TABLE]

In this set of equations, $E$ and $H$ are $1$-forms; $D$, $B$ and $J$ are $2$-forms; $\rho$ is the only $3$-form. Electric potential $\phi$, which is not shown above, is a [math]-form. Operator $d$ is the exterior derivative, and it takes $k$-form to $(k+1)$-form. Intuitively, $k$-form is an integrand, which can be integrated over $k$-D space. For example, a line integral of $1$-form $E$111Generally speaking, with a $k$-form $\omega$ and $k$-D oriented domain ${{\mathcal{D}}}$, we use $\int_{{{\mathcal{D}}}}\omega$ to denote the integral of $\omega$ over ${{\mathcal{D}}}$., $\int_{a}^{b}E=\int_{a}^{b}{{\mathbf{E}}}\cdot d{{\mathbf{l}}}$, leads to the potential difference from point $a$ to $b$. Indeed, the calculus of differential forms has significant advantages in illustrating the theory of EM theory compared to traditional vector analysis, see [2, 3, 4].

Discrete exterior calculus (DEC) provides a numerical treatment of differential forms [5, 6], which means that Equation (1) can be solved directly. In fact, since we are only able to obtain discrete values in numerical calculation, instead of solving for exact forms, the integral of unknown forms on finite line, area or volume are formulated and solved in DEC. Compared to FDTD, DEC can be implemented on unstructured simplicial mesh as in FEM, e.g. triangular mesh in $2$-D and tetrahedral mesh in $3$-D. It should be noted that DEC can be implemented on any mesh as long as circumcenter dual exists, such as regular Yee grid and hexahedral (layered triangular) grid. Therefore, this method is more adaptable over complex structures. In fact, the FDTD method can also be viewed as DEC method on Yee grid. In FDTD, the vectorial fields in fact are average values over edges, or square faces. These average values can also be thought of as small integrals, which is the origin of finite integration technique (FIT) or finite volume technique [7, 8]. Moreover, in contrast to FEM, this method exactly preserves important structural features of Maxwell’s equations, e.g. Gauss’s law $\nabla\cdot{{\mathbf{D}}}=\rho$ and ${{\mathbf{E}}}=-\nabla\phi$. Besides, since $\nabla\times\nabla=0$ and $\nabla\cdot\nabla\times=0$ are naturally and exactly preserved in DEC, which means that this method will not give rise to spurious solutions due to spurious charge [9, 10].

It should be pointed out that, the numerical work with unstructured grid based on FIT mainly dealt with $2$-D problems [11, 12]. More importantly, the material matrix (Hodge star) needs special treatment to be symmetric [13]. The reason is that FIT uses barycentric dual; then neighboring field values are needed in interpolation for building constitutive relation. In contrast, DEC adopts circumcenter dual, or Voronoi dual, which leads to the orthogonality of primal and dual elements. Therefore, the Hodge star operators, which represent constitutive relations, are diagonal matrices (for isotropic material or anisotropic material with diagonal ${\boldsymbol{\overline{\epsilon}}}$ and ${\boldsymbol{\overline{\mu}}}$). Since the Hodge star operators, the metric dependent parts, are closely related to the value of ${\boldsymbol{\overline{\epsilon}}}$ and ${\boldsymbol{\overline{\mu}}}$, on the other hand, a novel design of ${\boldsymbol{\overline{\epsilon}}}$ and ${\boldsymbol{\overline{\mu}}}$ can equivalently change the metric of system, which is merit of transformation optics [14]. There have been some efforts to apply DEC to computational electromagnetics. However, some limitations, mainly in two aspects, still remains in these work. First, the mesh is not totally unstructured. For example, there has been work based on polyhedral mesh [15, 16, 17] (layered triangular mesh for $3$-D) and partly structured nonuniform grids [18]. These special meshes largely limit the use of DEC method. Second, the previous calculations are mainly in time domain with homogeneous medium.

In our present work, we formulate and numerically calculate several kinds of electromagnetic problems in frequency domain, such as inhomogeneous waveguides (microstructured fibers), photonic crystals, and inhomogeneous resonant cavities. We also show that the charge continuity relation is exactly preserved with DEC, which prevents the generation of spurious charge and is nontrivial in FEM simulation [19]. Effective signed dual volume is introduced to incorporate material information to construct Hodge star operators with arbitrary simplicial mesh. It should be pointed out that if Delaunay triangulation is performed, the positivity of dual volumes can be guaranteed [20], which leads to a positive definite Laplacian operator in free space. But for frequency domain analysis, this is not a forced requirement. We also illustrate the construction and implementation of various boundary conditions in frequency domain, including perfect magnetic conductor (PMC), perfect electric conductor (PEC), first-order and second-order absorbing boundary conditions (ABCs), and periodic boundary conditions.

In Section 2, the framework of DEC is introduced and the governing equation for $3$-D and $2$-D problems are also formulated. Then in Section 3, we present in detail how we construct the Hodge star operators. In Section 4, various boundary conditions in frequency domain are discussed in detail. After every element of Maxwell’s equations is discussed and presented with the language of DEC, some numerical examples are shown in the last section. We apply this method to solve the decoupled or coupled modes in waveguides or optical fibers, the band diagram of a photonic crystal, and resonant frequencies of $3$-D inhomogeneous resonators with both closed and open boundary.

2 Maxwell’s Equations with DEC

As mentioned above, the FDTD method can also be understood with DEC by viewing the average field values as a series of finite integrals. Just as in the Yee grid, it has been shown that the Maxwell’s equations can also be expressed with DEC based on a simplicial mesh [5, 6, 17]. In this section, after a brief introduction to the framework of DEC, the governing equations for $3$-D and $2$-D problems in frequency domain are both formulated.

2.1 DEC approach to three-dimensional problems

As mentioned above, DEC is based on a simplicial mesh. Mathematically, $k$-simplex refers to the convex hull of $(k+1)$ vertices222The convex hull of a finite point set $S=\{x_{i}\}_{i=1}^{N}\subseteq\mathbb{R}^{d}$ refers to a region defined as Conv(S)=\left\{\sum_{i=1}^{N}c_{i}x_{i}\bigg{|}c_{i}\geq 0,\forall i;\sum_{i=1}^{N}c_{i}=1\right\}. in $k$-D space, e.g. points in [math]-D, lines in $1$-D, triangles in $2$-D and tetrahedrons in $3$-D. Figure 1 shows the comparison between a regular Yee grid and a general simplicial mesh in $2$-D. In DEC, the circumcenters are adopted to form the dual mesh. Therefore, a dual edge, denoted with ${{\mathcal{L}}}_{i}$, is orthogonal to its related primal edge $l_{i}$, e.g. $\overline{d_{1}d_{2}}$ and $\overline{p_{0}p_{1}}$333Here, $\overline{v_{1}v_{2}}$ is adopted to denote the line segment between vertices $v_{1}$ and $v_{2}$. in Figure 1. The dual face enclosed by $6$ dual edges and centered at $p_{0}$ is denoted with ${{\mathcal{A}}}_{p_{0}}$.

In $3$-D, the dual mesh is constructed by connecting circumcenters of nearest tetrahedrons. The dual edge, face and volume element in each tetrahedron is illustrated in Figure 2. By definition, the orthogonality between primal and dual elements are also satisfied.

Exterior calculus is induced by the Stokes’ theorem

[TABLE]

where $D$ and $\partial D$ represent any $k$-D domain and its $(k-1)$-D boundary, and $\omega$ is a $(k-1)$-form which can be integrated over any $(k-1)$-D space. Exterior derivative operator $d$, which transforms $(k-1)$-form $\omega$ to a $k$-form, acting on [math]-forms, $1$-forms and $2$-forms corresponds to gradient, curl and divergence operation in $3$-D space of vector calculus, respectively. In fact, Stokes’ theorem is a generalization of Newton’s, Green’s and Gauss’ theorem. However, in practice, both domain $D$ and its boundary $\partial D$ are discretized. And DEC is based on the discretized version of this theorem. In the following, we take Faraday’s law $dE=i\omega B$ as an example to illustrate this. Suppose primal face element $A_{i}$ (triangle) is enclosed by three edges, then the Stokes’ theorem on $A_{i}$ is written as

[TABLE]

where $d_{i,j}^{(1)}\in\{+1,-1\}$ depends on if the direction of $l_{j}$ coincides with the direction of $A_{i}$, which can be clockwise or anti-clockwise. The superscript $(1)$ means this relation is about edges.

Then we can introduce the discrete counterpart of differential forms, namely cochains444Chain, on the other hand, is used to denote a linear combination of simplices. In practice, a $k$-chain refers to a discrete domain a $k$-cochain can be integrated on., as column vectors, e.g. $1$-cochain $\bm{E}=(E_{1},\dots,E_{N_{1}})^{T}$ and $2$-cochain $\bm{B}=(B_{1},\dots,B_{N_{2}})^{T}$. Noted that $N_{1}$ and $N_{2}$ are the number of primal edge and face elements. The elements of these two cochains are defined as

[TABLE]

where $\hat{l}_{j}$ is the unit vector along $l_{j}$, and $\hat{n}_{i}$ is the unit vector normal to $A_{i}$. It should be pointed out that although traditionally cochain $\bm{E}$ and $\bm{B}$ are defined on primal mesh and are called primal cochains, they can also be defined on edges and faces of dual mesh. Then the relations described in Equation (3) can be represented as

[TABLE]

where $\overline{\mathbf{d}}^{(1)}$, as an $N_{2}\times N_{1}$ matrix, is the discrete version of exterior derivative or co-boundary operator, and the superscript $(1)$ means that the operator only acts on primal $1$-cochains. The $(i,j)$ element of this matrix is defined as

[TABLE]

Similarly, we can also define other co-boundary operators $\overline{\mathbf{d}}^{(2)}$ (an $N_{3}\times N_{2}$ matrix) acting on primal $2$-cochains and $\overline{\mathbf{d}}^{(0)}$ (an $N_{1}\times N_{0}$ matrix) acting on primal [math]-cochains as

[TABLE]

where $A_{j}$ and $V_{i}$ are the $j$-th primal face and $i$-th primal volume (a tetrahedron), and $p_{j}$ is the $j$-th primal vertex.

Therefore, since the integral form of the divergence relation $dB=0$ can be represented as

[TABLE]

where $d_{i,j}^{(2)}\in\{+1,-1\}$ depends on if the normal direction of each face $A_{j}$ is pointing outside the volume $V_{i}$. Then this divergence relation can be described as

[TABLE]

Different from primal cochains $\bm{E}$ and $\bm{B}$, cochains $\bm{D}$, $\bm{H}$, $\bm{J}$, $\bm{\rho}$ are defined on dual mesh and are called dual cochains. They are defined as

[TABLE]

where

[TABLE]

where ${{\mathcal{L}}}_{i},{{\mathcal{A}}}_{i}$, and ${{\mathcal{V}}}_{i}$ are the $i$-th dual edge, face and volume element, as illustrated in Figure 2. Noted that due to orthogonality between primal and dual elements, $\hat{l}_{i}$, the unit vector along $i$-th primal edge $l_{i}$, is also the normal vector of dual face ${{\mathcal{A}}}_{i}$, and $\hat{n}_{i}$, the unit normal vector through $i$-th primal face $A_{i}$, is parallel to ${{\mathcal{L}}}_{i}$. Also noted that the number of dual edges is the same with the number of primal faces $N_{2}$, and the number of dual faces is the same with the number of primal edges $N_{1}$. And this explains the length of these dual cochains.

To build relations between dual cochains as in (5) and (10), co-boundary operators $\overline{\mathbf{d}}_{\text{dual}}^{(k)}$ on dual mesh needs to be introduced. In fact there is a relation between co-boundary operators on $n$-dimensional primal and dual mesh [5, 6]

[TABLE]

where superscript $T$ refers to transpose. Specifically in $3$-D,

[TABLE]

Therefore, we can represent the discrete version of $dH=-i\omega D+J$ and $dD=\rho$ as

[TABLE]

Moreover, since $\overline{\mathbf{d}}^{(1)}\overline{\mathbf{d}}^{(0)}=0$ and $\overline{\mathbf{d}}^{(2)}\overline{\mathbf{d}}^{(1)}=0$ exactly holds for both primal and dual mesh, another $\overline{\mathbf{d}}^{(2)}_{\text{dual}}$ left multiplying to (15) leads to the charge continuity relation

[TABLE]

which is naturally satisfied. An intuitive reason is that a $(k-2)$ simplex always appears twice with different signs in the boundary of the boundary of a $k$-simplex [6]. In fact, the exact preservation of this relation will prevent the generation of spurious charge in electromagnetic simulation. Therefore, only Equations (5) and (15) are independent.

Moreover, the constitutive relation can be illustrated by two characterized Hodge star operators, $N_{1}\times N_{1}$ matrix $\star^{(1)}_{\epsilon}$ and $N_{2}\times N_{2}$ matrix $\star^{(2)}_{\mu^{-1}}$. The construction scheme will be shown in the next section. These two operators map primal cochains $\bm{E}$ and $\bm{B}$ to dual cochains $\bm{D}$ and $\bm{H}$ as

[TABLE]

Then combining Equations (5), (15) and (18), we can obtain the equation for primal cochain $\bm{E}$ as

[TABLE]

where $k_{0}=\omega\sqrt{\mu_{0}\epsilon_{0}}$ is the vacuum wave number. Then Equation (19) can be used to solve eigen modes in sourceless $3$-D space or the excitation field by a given current source.

It should be pointed out that if the diagonal elements of $\star^{(2)}_{\mu^{-1}}$ are all positive, the semi-positiveness of $(\overline{\mathbf{d}}^{(1)})^{T}\star^{(2)}_{\mu^{-1}}\overline{\mathbf{d}}^{(1)}$ is guaranteed. For example, with arbitrary $1$-cochain $\bm{e}$,

[TABLE]

where $\bm{b}=\overline{\mathbf{d}}^{(1)}\bm{e}$. But if there is a [math]-cochain $\bm{\phi}$, such that $\bm{e}=\overline{\mathbf{d}}^{(0)}\bm{\phi}$, cochain $\bm{b}$ is an exact zero cochain. Therefore, $(\overline{\mathbf{d}}^{(1)})^{T}\star^{(2)}_{\mu^{-1}}\overline{\mathbf{d}}^{(1)}$ has a large null space, just like $\nabla\times\nabla\times$ operator.

2.2 Two-dimensional case

Two-dimensional mesh, assumed in $x$-$y$ plane, can be thought of a degenerated $3$-D layered triangular mesh (extended along $z$-direction), or hexahedral mesh. In this $3$-D mesh, there are edges both on $x$-$y$ plane and along $z$-direction, and there are faces both on $x$-$y$ plane and parallel to $z$-direction. After degeneration to a $2$-D mesh, those edges and faces on $x$-$y$ plane remain the same, while the edges along $z$-direction become points and the faces parallel to $z$-direction become edges. Therefore, previous $1$-cochains associated with edges along $z$-direction become [math]-cochains in $2$-D mesh, and $2$-cochains associated with faces parallel to $z$-direction become $1$-cochains. In $2$-D problems, the fields are separated into $z$-component and transverse component, and each component corresponds to one cochain. By convention, we still place cochains associated with field ${{\mathbf{E}}}({{\mathbf{r}}})$ and ${{\mathbf{B}}}({{\mathbf{r}}})$ on primal mesh, and cochains related with field ${{\mathbf{D}}}({{\mathbf{r}}})$ and ${{\mathbf{H}}}({{\mathbf{r}}})$ on dual mesh. More specifically, the corresponding relations between fields and cochains are listed below:

$z$-component:

1. (a)

$E_{z}({{\mathbf{r}}})\hat{z}$ $\Rightarrow$ primal [math]-cochains $\bm{E}_{z}=(E_{z,1},\dots,E_{z,N_{0}})^{T}$; 2. (b)

$B_{z}({{\mathbf{r}}})\hat{z}$ $\Rightarrow$ primal $2$-cochain $\bm{B}_{z}=(B_{z,1},\dots,B_{z,N_{2}})^{T}$; 3. (c)

$D_{z}({{\mathbf{r}}})\hat{z}$ $\Rightarrow$ dual $2$-cochain $\bm{D}_{z}=(E_{z,1},\dots,E_{z,N_{0}})^{T}$; 4. (d)

$H_{z}({{\mathbf{r}}})\hat{z}$ $\Rightarrow$ dual [math]-cochains $\bm{H}_{z}=(H_{z,1},\dots,H_{z,N_{2}})^{T}$. 2. 2.

transverse component:

1. (a)

${{\mathbf{E}}}_{s}({{\mathbf{r}}})$ $\Rightarrow$ primal $1$-cochain $\bm{E}_{s}=(E_{s,1},\dots,E_{s,N_{1}})^{T}$; 2. (b)

${{\mathbf{B}}}_{s}({{\mathbf{r}}})$ $\Rightarrow$ primal $1$-cochain $\bm{B}_{s}=(B_{s,1},\dots,B_{s,N_{1}})^{T}$; 3. (c)

${{\mathbf{D}}}_{s}({{\mathbf{r}}})$ $\Rightarrow$ dual $1$-cochain $\bm{D}_{s}=(D_{s,1},\dots,D_{s,N_{1}})^{T}$; 4. (d)

${{\mathbf{H}}}_{s}({{\mathbf{r}}})$ $\Rightarrow$ dual $1$-cochain $\bm{H}_{s}=(H_{s,1},\dots,H_{s,N_{1}})^{T}$.

For [math]-cochain $\bm{E}_{z}$ or $\bm{H}_{z}$, each element is just the field value at corresponding primal or dual vertex. For $2$-cochain $\bm{B}_{z}$ or $\bm{D}_{z}$, each element is the field integral on one primal face $A_{i}$ or dual face ${{\mathcal{A}}}_{i}$

[TABLE]

Also noted that, although both $\bm{E}_{s}$ and $\bm{B}_{s}$ are $1$-cochains on primal mesh, they are defined in different ways

[TABLE]

This is because $E_{s,i}$ denotes the potential change along $l_{i}$, while $B_{s,i}$ is the magnetic flux through $l_{i}$. In fact, $\bm{B}_{s}$ represents a $2$-cochain mesh in $3$-D defined on faces parallel to $z$-direction. Similarly, dual $1$-cochain $\bm{D}_{s}$ and $\bm{H}_{s}$ are also defined differently

[TABLE]

The operator $\overline{\mathbf{d}}^{(1)}$ and $\overline{\mathbf{d}}^{(0)}$ in $2$-D mesh is defined the same way as in Equation (6) and (8). For dual mesh, according to (13), we can obtain relations

[TABLE]

It should be pointed out that with $2$-D gradient operator $\nabla_{s}=\hat{x}\frac{\partial}{\partial x}+\hat{y}\frac{\partial}{\partial y}$, there are two possible operations on a transverse vector field: curl $\nabla_{s}\times$ and divergence $\nabla_{s}\cdot$. They are both represented by $\overline{\mathbf{d}}^{(1)}$ or $\overline{\mathbf{d}}^{(1)}_{\text{dual}}$, because they both acts on primal or dual $1$-cochains. However, observation of the $2$-D mesh in Figure 1 shows that primal edge $l_{j}$ enclosing primal face $A_{i}$ clockwise contradicts with $\hat{z}\times\hat{l}_{j}$ pointing outside $A_{i}$, and dual edge ${{\mathcal{L}}}_{j}$ (with direction $\hat{z}\times\hat{l}_{j}$) enclosing dual face ${{\mathcal{A}}}_{i}$ clockwise coincide with $\hat{l}_{j}$ pointing outside ${{\mathcal{A}}}_{i}$. Therefore, divergence operation $\nabla_{s}\cdot$ acting on primal $1$-cochains induces an extra negative sign to $\overline{\mathbf{d}}^{(1)}$.

The constitutive relations in $2$-D involve two components. Assuming the permittivity and permeability matrices are

[TABLE]

where ${\boldsymbol{\overline{\epsilon}}}_{s}$ and ${\boldsymbol{\overline{\mu}}}_{s}$ are $2\times 2$ tensors functions and their components are along the transverse directions. For simplicity, we assume ${\boldsymbol{\overline{\epsilon}}}_{s}=\epsilon_{s}({{\mathbf{r}}})\overline{\mathbf{I}}$ and ${\boldsymbol{\overline{\mu}}}_{s}=\mu_{s}({{\mathbf{r}}})\overline{\mathbf{I}}$, and we use $\epsilon_{z}$ and $\mu_{z}$ to replace $\epsilon_{zz}$ and $\mu_{zz}$. Then four Hodge star operators are needed, $\star^{(0)}_{\epsilon_{z}}$, $\star^{(1)}_{\epsilon_{s}}$, $\star^{(2)}_{\mu^{-1}_{z}}$, and $\star^{(1)}_{\mu^{-1}_{s}}$. The mapping relations are

[TABLE]

For homogeneous waveguides or inhomogeneous ones with $k_{z}=0$, e.g. $2$-D photonic crystals, TM and TE modes are decoupled and can be analyzed independently. The reduced wave equations of $E_{z}({{\mathbf{r}}})$ or $H_{z}({{\mathbf{r}}})$ for homogeneous waveguides can be written as:

[TABLE]

where $k_{s}^{2}=\omega^{2}\mu\epsilon-k_{z}^{2}$.

For TM modes, $\nabla_{s}E_{z}({{\mathbf{r}}})$ leads to a $1$-cochain on primal edges, which is $\overline{\mathbf{d}}^{(0)}\bm{E}_{z}$. For the next operator $\nabla_{s}\cdot$ to function appropriately, a Hodge star operator $\star^{(1)}$ (no subscript due to homogeneity) needs to be inserted to transform $\overline{\mathbf{d}}^{(0)}\bm{E}_{z}$ to it dual $1$-cochain. Therefore, $\nabla_{s}\cdot\nabla_{s}$ can be replaced with $[\overline{\mathbf{d}}_{\text{dual}}^{(1)}\star^{(1)}\overline{\mathbf{d}}^{(0)}]$, or $-[(\overline{\mathbf{d}}^{(0)})^{T}\star^{(1)}\overline{\mathbf{d}}^{(0)}]$. Finally the TM mode governing equation of primal [math]-cochain $\bm{E}_{z}$ can be written with DEC as:

[TABLE]

For TE modes, the equation for dual [math]-cochain $\bm{H}_{z}$ can be similarly rewritten as:

[TABLE]

For inhomogeneous waveguide with $k_{z}=0$, Equation (32) and (33) need to be modified accordingly.

[TABLE]

where $k_{0}^{2}=\omega^{2}\mu_{0}\epsilon_{0}$.

For inhomogeneous waveguide with nonzero $k_{z}$, TE and TM modes are coupled to each other, and the resulting modes are hybrid. The governing equation of the transverse field ${{\mathbf{E}}}_{s}({{\mathbf{r}}})$ and ${{\mathbf{H}}}_{s}({{\mathbf{r}}})$ can be written as [1]:

[TABLE]

and

[TABLE]

For example, we rewrite Equation (36) with DEC one by one. As mentioned, transverse field ${{\mathbf{E}}}_{s}({{\mathbf{r}}})$ is represented as a primal $1$-cochain $\bm{E}_{s}$. And ${\boldsymbol{\overline{\epsilon}}}_{s}$, ${\boldsymbol{\overline{\mu}}}_{s}^{-1}$, $\epsilon_{z}$, and $\mu_{z}^{-1}$ can be replaced by four Hodge star operators. It should also be pointed out that although $\hat{z}\times$ rotates the vector field by $90$ degree, the corresponding cochain remains the same. For differential operators, in the first term of (36), first $\nabla_{s}\times$ acting on primal $1$-cochains $\bm{E}_{s}$ refers to $\overline{\mathbf{d}}^{(1)}$, and results in a primal $2$-cochain. Then subsequent ${\boldsymbol{\overline{\mu}}}_{s}\cdot\hat{z}\times\nabla_{s}\times\mu_{z}^{-1}$ acting on this primal $2$-cochain555For compactness, we only defined Hodge star operators mapping primal to dual cochains. If a reverse mapping is needed, we use the inverse of the corresponding Hodge star operator. leads to $(\star_{\mu_{s}^{-1}}^{(1)})^{-1}\overline{\mathbf{d}}^{(0)}_{\text{dual}}\star_{\mu_{z}^{-1}}^{(2)}$. In the second term of (36), $\nabla_{s}\cdot{\boldsymbol{\overline{\epsilon}}}_{s}$ acting on $\bm{E}_{s}$ refers to $\overline{\mathbf{d}}^{(1)}_{\text{dual}}\star^{(1)}_{\epsilon_{s}}$, and results in a dual $2$-cochain. Then the subsequent gradient operator $\nabla_{s}\epsilon_{z}^{-1}$ acting on this dual $2$-cochains leads to $\overline{\mathbf{d}}^{(0)}(\star_{\epsilon_{z}}^{(0)})^{-1}$. The third term of (36) only involves two Hodge star operators. Therefore, (36) can be transformed to

[TABLE]

Therefore, Equations (32) to (35) and Equation (38) can be applied with certain boundary conditions to solve eigen mode problems for $2$-D structures.

3 Hodge Star Operators

Hodge star operators defined in Section 2 map a primal cochain to its corresponding dual cochain. It should be pointed out the “orthogonal” dual and corresponding constitutive relation was first proposed in [21], then fully developed in the frame of DEC. From the definition for $E_{i}$ and $D_{i}$ in Equation (4) and (12), with locally constant field assumption, we can infer that there should be a one to one relation between them

[TABLE]

Here, $(\star^{(1)}_{\epsilon})_{i}$ is the $i$-th diagonal matrix element of $\star^{(1)}_{\epsilon}$, and $\epsilon_{i}$ is the average permittivity. Then $(\star^{(1)}_{\epsilon})_{i}$ is defined to be the ratio of dual face element’s area to primal edge element’s length, and multiplied by $\epsilon_{i}$. Note that, from now on, we use volume as a general term for length, area and volume in $1$-D, $2$-D and $3$-D, respectively. In this case, since the local value of $\epsilon$ needs to be included to obtain a characteristic Hodge star, and this is why $\epsilon$ is in the subscript. If we adopt this definition, (39) can be written as

[TABLE]

Here, $|{{\mathcal{A}}}_{i}|$ and $|l_{i}|$ denote the volume of dual face ${{\mathcal{A}}}_{i}$ and primal edge $l_{i}$. Similarly the relation between $H_{i}$ and $B_{i}$ defined in (4) and (12) reads

[TABLE]

where $|{{\mathcal{L}}}_{i}|$ and $|A_{i}|$ are the volume of dual edge ${{\mathcal{L}}}_{i}$ and primal face $A_{i}$, and $\mu_{i}$ is the local value of $\mu({{\mathbf{r}}})$. Then the Hodge star operator $\star^{(1)}_{\epsilon}$ and $\star^{(2)}_{\mu^{-1}}$ in $3$-D are constructed as below

[TABLE]

In a $2$-D case, from the definition in Equation (20)-(25), there are four constitutive relations involved

[TABLE]

where $E_{z,i}$ and $H_{z,i}$ are the field value at $i$-th primal and dual vertex.

Noted that, in practice, $\epsilon({{\mathbf{r}}})$ and $\mu({{\mathbf{r}}})$ are assumed to be constant in one tetrahedron (triangle in $2$-D). However, since generally dual edge elements or dual face elements do not belong to a unique tetrahedron, a weighted average needs to be performed, which will be discussed next. It should be pointed out that even with Delaunay triangulation, circumcenter may falls outside the tetrahedron (triangle in $2$-D). However, there are some problems for which Delaunay triangulation is not a good idea. In this volume calculation, we also put this special case into consideration and show that our scheme is suitable for a general triangular mesh (or tetrahedral mesh). We will start from a $2$-D case.

3.1 Volume of dual cells in 2D

Then in a simple $2$-D mesh, as shown in Figure 3, dual edge elements, such as ${{\mathcal{L}}}_{d_{1}d_{6}}$, are composed by two components from two neighboring triangles ($\bigtriangleup_{p_{0},p_{1},p_{2}}$ and $\bigtriangleup_{p_{0},p_{6},p_{1}}$), while dual face elements are composed by several parts from triangles sharing the same primal vertex. For example, dual face ${{\mathcal{A}}}_{p_{0}}$ associated with $p_{0}$ is composed by 6 parts $A_{1},\dots,A_{6}$. Then we can obtain the length of dual edge ${{\mathcal{L}}}_{d_{1}d_{6}}$ by adding its two components, and the area of gray shaded $2$-cell by summing over all its six parts:

[TABLE]

Here the second equality sign in (46) is because that the dual face ${{\mathcal{A}}}_{p_{0}}$ can also be decomposed into $6$ triangles, and each of them has one dual edge as bottom edge and $p_{0}$ as another vertex. In this way, the volume of dual edge and face elements can be calculated systematically.

With inhomogeneous material information, a weighted average needs to be performed to obtain effective dual volume. For example, if each region $A_{i}$ in Figure 3 is associated with a distinct permittivity value $\epsilon_{i}$, then effective volume of ${{\mathcal{L}}}^{*}_{d_{1}d_{6},\epsilon}$ and ${{\mathcal{A}}}^{*}_{p_{0},\epsilon}$ can be obtained as

[TABLE]

Therefore, Hodge star operators $\star^{(0)}_{\epsilon_{z}}$ and $\star^{(1)}_{\epsilon_{s}}$ can be constructed as follows

[TABLE]

Noted that, $\epsilon_{z}$ and $\epsilon_{s}$ can be chosen differently for anisotropic material.

Similarly, we can also obtain effective dual volume based on localized value of $\mu^{-1}$

[TABLE]

Here $A_{d_{i}}$ is the primal triangle with circumcenter $d_{i}$.

Then Hodge star operators $\star^{(2)}_{\mu_{z}^{-1}}$ and $\star^{(1)}_{\mu_{s}^{-1}}$ can be constructed as

[TABLE]

In the above, $\mu_{z}$ and $\mu_{s}$ can also be chosen differently.

If an inverse mapping is needed, such as $\bm{D}$ to $\bm{E}$ and $\bm{H}$ to $\bm{B}$, we can just use the direct inverse of these Hodge star operators since they are all diagonal.

3.1.1 Special case: dual vertex is outside the triangle

In Figure 4(a), a special mesh with one circumcenter ($d_{1}$) falling outside the associated triangle is shown, which is frequently encountered even with Delaunay triangulation. Then in (45), $s_{1,1}$ is negative, while the calculation for ${{\mathcal{A}}}_{p_{0}}$ based on the second equality in (46) stays unchanged. It should be noted that since $s_{1}$, the volume of dual edge $\overline{d_{6}d_{1}}$, is still positive, this mesh is still Delaunay triangulation ($p_{2}$ is outside the circumcircle of $\bigtriangleup_{p_{0},p_{1},p_{6}}.$).

The scheme described by Equation (46) also works for an even more “twisted” case, a non-Delaunay triangulation, shown in Figure 4(b). In fact, the total area of the dual face element can be found to be the difference between left cyan region and right gray region. From the arrow directions, we can see that the left cyan region has an anti-clockwise orientation, while the right gray region has a clockwise orientation. This means that they are signed volume. Therefore, we need to subtract the clockwise one to obtain a corrected dual volume for this dual face.

3.2 Volume of dual cells in 3D

The calculation for the volume of dual elements is far more complicated in $3$-D. We can calculate the volume of dual faces one tetrahedron by one tetrahedron, shown in Figure 5. For example, the area for the shaded dual face in Figure 5(a) can be obtained as:

[TABLE]

Here the subscript $\overline{p_{0}p_{1}}$ is used because this dual face is associated with primal edge $\overline{p_{0}p_{1}}$. Then the total volume for dual edge and face elements can be obtained by adding the values in a single tetrahedron appropriately. For example, as shown in Figure 5(b), the total area of shaded dual $2$-cell is composed by 6 components from 6 tetrahedrons.

As long as the volume of dual faces are obtained, the volume for dual $3$-cells can be obtained simply. Figure 6 shows one component of the dual $3$-cell centered at vertex $p_{0}$. We can see that each of this component is a generalized cone associated with one dual face, or primal edge (the red line segment). Then the total volume for this dual $3$-cell ${{\mathcal{V}}}_{p_{0}}$ can be obtained by summing over all its cone components as:

[TABLE]

Here we have applied the volume formula of a cone shape structure: $V=\frac{1}{3}S\cdot h$, with bottom surface area $|{{\mathcal{A}}}_{i}|$ and height $\frac{1}{2}|l_{i}|$.

Incorporating the material information, like localized values of $\epsilon$ and $\mu^{-1}$, the effective volume can be obtained as

[TABLE]

Therefore, the Hodge star operators $\star^{(1)}_{\epsilon}$ and $\star^{(2)}_{\mu^{-1}}$ can be constructed from Equation (42).

3.2.1 Special case: dual vertex is outside the tetrahedron

As shown in Figure 7, since some dual vertices are outside their associated tetrahedrons, so the dual $2$-cell seems to be very “twisted”. To obtain the corrected volume of dual cells for this special case, extra negative signs need to be assigned to some dual edges as we did in $2$-D,

Unlike in $2$-D, the complex geometry of dual face shown in Figure 7(b) is originated from two reasons: dual vertices fall outside tetrahedrons and circumcenters of primal faces may also fall outside. This is why the red primal edge does not even go through its dual face. Two signs, instead of just one, are needed to take these two factors into account simultaneously. For example, as in Figure 5(a), the area of triangle $\bigtriangleup_{d,f_{2},e_{0}}$, $|\bigtriangleup_{d,f_{2},e_{0}}|$, is written as:

[TABLE]

But if $d$ is outside the tetrahedron, or $f_{2}$ is outside the triangle $\bigtriangleup_{p_{0},p_{1},p_{3}}$, extra signs need to be multiplied:

[TABLE]

Then the following calculation for the volume of dual $3$-cells in Equation (54) stays the same.

4 Boundary Conditions

In this section, we will show how to implement various of boundary conditions. The direct reason for a boundary condition is that the dual mesh is incomplete and truncated by primal mesh at the boundary, as illustrated in Figure 8. In fact, PEC, PMC, and first-order ABC in time domain analysis have been investigated in [22]. We will show a different formulation for these boundary conditions in frequency domain, and in addition, we will also formulate second-order ABC and periodic boundary condition. For illustration, we use Equation (32) for $2$-D and Equation (19) for $3$-D. In fact, there is already one boundary condition embedded in $(\overline{\mathbf{d}}^{(0)})^{T}$ for (32) and $(\overline{\mathbf{d}}^{(1)})^{T}$ for (19). We first examine which kind this default boundary condition is.

4.1 Default boundary condition: PMC

From (32), $\left[\star^{(0)}\overline{\mathbf{d}}^{(0)}\right]\bm{E}_{z}$ generates a dual $1$-cochain, and we denote it with $\bm{H}_{s}$. Then $\left[(\overline{\mathbf{d}}^{(0)})^{T}\right]\bm{H}_{s}$ leads to a dual $2$-cochain, which we represent it with $\bm{D}_{z}$. Next we use $D_{z,p_{j}}$ and $H_{s,\overline{v_{i}v_{j}}}$ to denote cochains’ value on dual face ${{\mathcal{A}}}_{p_{j}}$ and its surrounding dual edge $\overline{v_{i}v_{j}}$ ($v_{i},v_{j}$ are two vertices) in Figure 8. Without discretization, the relation between field $D_{z}({{\mathbf{r}}})\hat{z}$ and ${{\mathbf{H}}}_{s}({{\mathbf{r}}})$ is

[TABLE]

Then we can conclude the relation below from the integral form of this relation (Ampere’s law):

[TABLE]

Here, $\partial{{\mathcal{A}}}_{p_{2}}$ denotes the edges of dual face ${{\mathcal{A}}}_{p_{2}}$. However, with discrete relation $\bm{D}_{z}=(\overline{\mathbf{d}}^{(0)})^{T}\bm{H}_{s}$, only the first three terms in Equation (60) are included. The last two terms $H_{s,\overline{e_{1}p_{2}}}$, and $H_{s,\overline{p_{2}e_{2}}}$ need to be determined by certain boundary conditions. Therefore, without extra boundary condition implemented, the last two terms of (60) are set to be zero implicitly. And this is the perfect magnetic conductor (PMC) boundary condition with zero tangential magnetic field. For a $3$-D equation (19), this is also true. In the language of partial differential equation (PDE), this is the homogeneous natural (Neumann) boundary condition.

4.2 PEC and Dirichlet boundary condition

Perfect electric conductor (PEC) boundary condition is a little different, and it refers to essential boundary condition in PDE. For TM modes, $\bm{E}_{z}$ defined on the boundary vertices all have zero value. Next we use subscript I to denote field or value strictly inside the boundary, and subscript $B$ to denote field or value on the boundary. To consider the vertices strictly inside only and keep the constructed derivative operator, we can use a projection operator $\overline{\mathbf{P}}_{E_{z}}$. Suppose there are $N_{0,\text{I}}$ vertices strictly inside out of total $N_{0}$ vertices. Then, $\overline{\mathbf{P}}_{E_{z}}$ is a $N_{0}\times N_{0,\text{I}}$ matrix, and it can insert zero boundary values to map a $\bm{E}_{z,\text{I}}$ cochain to the total $\bm{E}_{z}=\{\bm{E}_{z,\text{I}};\bm{E}_{z,\text{B}}\}$ cochain.

[TABLE]

where $\overline{\mathbf{I}}$ is $N_{0,\text{I}}\times N_{0,\text{I}}$ identity matrix. Then (32) with PEC boundary condition should be adjusted as:

[TABLE]

For a Dirichlet boundary condition, e.g. in scattering problem, $E_{z,\text{scat}}=-E_{z,\text{inc}}$ must be satisfied for scattered field such that $E_{z,\text{total}}=0$ on a PEC boundary. Therefore, by using operator $\overline{\mathbf{P}}_{E_{z}}$, the scattered field cochain is represented as

[TABLE]

Then Equation (62) for a $2$-D scattering problem is adjusted as

[TABLE]

For Equation (38), PEC boundary condition can be implemented similarly by introducing projection operator $\overline{\mathbf{P}}_{E_{s}}$. Similarly, $\overline{\mathbf{P}}_{E_{s}}$ is a $N_{1}\times N_{1,I}$ matrix with $N_{1}$ and $N_{1,I}$ representing the number of all edges and edges strictly inside. Then with $\bm{E}=\{\bm{E}_{s,I};\bm{E}_{s,B}\}$, (38) with PEC is written as

[TABLE]

with

[TABLE]

where $\overline{\mathbf{I}}$ is $N_{1,I}\times N_{1,I}$ identity matrix.

However, this is not the only way to implement PEC boundary condition. Inspired by the embedded PMC boundary condition discussed above, we can conclude that if we switch primal and dual cochains, PEC is implied instead of PMC. More specifically, we can set $\bm{E}$, $\bm{B}$ as dual cochains and $\bm{H}$, $\bm{D}$, $\bm{J}$ as primal cochains instead. Then equation of cochain $\bm{H}$ without extra boundary condition will imply PEC boundary condition.

Especially for $3$-D equation (19), although PEC boundary condition can be implemented simply by ignoring the boundary elements, it will result in considerable error. Instead, we can write equivalent equation for cochain $\bm{H}$ as

[TABLE]

Here, $\star^{(2)}_{\epsilon^{-1}}$ and $\star^{(1)}_{\mu}$ are defined in a similar way with $\star^{(2)}_{\mu^{-1}}$ and $\star^{(1)}_{\epsilon}$. Then, without extra boundary condition inserted on (67), PEC is implicitly implemented.

This approach can also applies to $2$-D problems. Instead of solving Equation (33) of dual [math]-cochain $\bm{H}_{z}$ with PEC boundary condition for TE modes, we can set $\bm{H}_{z}$ as primal [math]-cochain and rewrite Equation (31) with DEC as

[TABLE]

In the language of PDE, by switching primal and dual cochains, PEC boundary condition is changed from an essential boundary condition to natural boundary condition.

4.3 Periodic boundary condition

Periodic boundary condition can also be implemented in our method with a procedure similar to FEM [23]. For a periodic structure with square unit cell, shown in Figure 9.

A projection matrix $\overline{\mathbf{P}}$ is needed to reduce unknowns on the boundary. If we use subscript $I$ to denote vertices strictly inside, subscripts $L$, $R$, $B$, $T$ for four boundaries, and subscripts $C_{i}$ ($i=1,2,3,4$) for four corners, then the cochain $\bm{E}_{z}$ can be divided into nine components

[TABLE]

This can be reduced to a vector only containing independent unknowns by applying periodic boundary conditions:

[TABLE]

where

[TABLE]

and

[TABLE]

where $\Psi_{x}=k_{x}L$ and $\Psi_{y}=k_{y}L$ are the phase shift between adjacent unit cells along $x$-axis and $y$-axis. It should be noted that operator $\overline{\mathbf{P}}$ will have different form for lattice with different unit cell. Then with periodic boundary condition, Equation (34) for TM modes can be written as

[TABLE]

For TE modes, using the same technique introduced in PEC boundary condition, we can place [math]-cochain $\bm{H}_{z}$ on primal mesh, and defined $\bm{H}_{z,\text{reduced}}$ in the same way as above. Then equation for TE modes with periodic boundary condition is written as

[TABLE]

Therefore, TM and TE band diagrams in photonic crystals can be obtained by solving Equations (73) and (74) with different $k_{x}$ and $k_{y}$.

4.4 Two-dimensional ABCs

Absorbing boundary conditions (ABCs) can also be implemented easily, and they represent inhomogeneous natural (Neumann) boundary conditions. For first-order ABC, we should have a proportional relation between tangential component of ${{\mathbf{H}}}_{s}({{\mathbf{r}}})$ and $E_{z}({{\mathbf{r}}})\hat{z}$ [9, 23].

[TABLE]

Here $\eta=\sqrt{\frac{\epsilon}{\mu}}$ is the impedance, and $\kappa$ is the curvature at the boundary ([math] for flat boundary). Then the last two terms of Equation (60) can be approximated as

[TABLE]

Here, $|{{\mathcal{L}}}_{p_{2}}^{\partial}|=|\overline{e_{1}p_{2}}|+|\overline{p_{2}e_{2}}|$ is the total length of dashed lines in Figure 8, and superscript $\partial$ stands for boundary. Then with first-order ABC, (32) should be adjusted as:

[TABLE]

where

[TABLE]

Here, $\star^{\kappa}$ and $\star^{\partial}$ are constructed as $N_{0}\times N_{0}$ matrices with nonzero diagonal elements only for boundary vertices, and subscript $i$ refers to a vertex on the boundary. Therefore, they are highly sparse matrices and only act on boundary elements.

Second-order ABC is derived in [1, 23] as

[TABLE]

Here $\frac{\partial^{2}E_{z}}{\partial s^{2}}$ refers to second-order derivative on the boundary. The first term of (79) can be implemented the same way as in first-order ABC. To interpret the second term of (79) with DEC, we introduce the derivative operator $\overline{\mathbf{d}}^{(0,\partial)}$ and Hodge star operator $\star^{1,\partial}$ confined on the boundary. Here $\overline{\mathbf{d}}^{(0,\partial)}$ is only nonzero for boundary primal edge and point elements, and is defined the same way as $\overline{\mathbf{d}}^{(0)}$, while $\star^{1,\partial}$ is defined as

[TABLE]

Here subscript $l_{j}$ refers to a primal edge on the boundary. Then $\frac{\partial^{2}E_{z}}{\partial s^{2}}$ corresponds to a dual $1$-cochain only on the boundary. More specifically, as shown in Figure 8, we can derive

[TABLE]

Therefore, $\frac{\partial^{2}E_{z}}{\partial s^{2}}$ can be represented with dual $1$-cochain $[(\overline{\mathbf{d}}^{(0,\partial)})^{T}\star^{1,\partial}\overline{\mathbf{d}}^{(0,\partial)}]\bm{E}_{z}$. Then with second-order ABC, (32) is adjusted as

[TABLE]

4.5 Three-dimensional ABCs

For a $3$-D case, the simplest ABC is the Sommerfeld radiation condition [23]

[TABLE]

This relation can be implemented in a similar way as in above $2$-D case. With this first-order ABC, (19) is adjusted as

[TABLE]

with

[TABLE]

Similarly, $\star^{(1,\partial)}$ is a $N_{1}\times N_{1}$ matrix only containing diagonal elements for boundary edges, and ${{\mathcal{L}}}_{i}^{\partial}$ is the surface dual edge shown in Figure 10.

The second-order ABC [24] requires the electric field on the boundary to satisfy

[TABLE]

where $\nabla_{t}$ is the surface tangential gradient operator, the subscript $n$ represents the normal component on the surface, and parameter $\beta$ is defined as

[TABLE]

where $\kappa$ is the curvature of the boundary surface (zero for flat surfaces).

While in the second term of (86), $\hat{n}\hat{n}\cdot(\nabla\times{{\mathbf{E}}})$, or $\hat{n}(\nabla\times{{\mathbf{E}}})_{n}$ can be simply represented by $\overline{\mathbf{d}}^{(1,\partial)}\bm{E}$, where $\overline{\mathbf{d}}^{(1,\partial)}$ is defined the same way with $\overline{\mathbf{d}}^{(1)}$ only containing relation between boundary surface faces and edges. It should be noted that on the $2$-D surface $\partial\Omega$, $\hat{n}(\nabla\times{{\mathbf{E}}})_{n}$ leads to a primal $2$-cochain $\overline{\mathbf{d}}^{(1,\partial)}\bm{E}$. Therefore, for the second curl operator $\nabla\times$ to operate on this primal $2$-cochain, a Hodge star operator $\star^{2,\partial}$ needs to be implemented first. Here $\star^{(2,\partial)}$ is defined as

[TABLE]

where surface primal face $A_{i}^{\partial}$ represents a triangle on boundary surface. Then $\nabla\times$ corresponds to $(\overline{\mathbf{d}}^{(1,\partial)})^{T}\star^{(2,\partial)}$, because

[TABLE]

In the third term of (86), with the language of DEC, surface divergence $\nabla_{t}\cdot$ is denoted by the transpose of $\overline{\mathbf{d}}^{(0,\partial)}$ and a surface Hodge star operator $\star^{(1,\partial)}$ together acting on primal $1$-cochain $\bm{E}$. Then $[(\overline{\mathbf{d}}^{(0,\partial)})^{T}\star^{(1,\partial)}]\bm{E}$ represents a dual $2$-cochain on boundary surface. For the surface gradient $\nabla_{t}$ to operate appropriately, another Hodge star operator $(\star^{(0,\partial)})^{-1}$ needs to be inserted to map this dual $2$-cochain to primal [math]-cochain. Hodge star $\star^{(0,\partial)}$ is defined only for surface points as

[TABLE]

where ${{\mathcal{A}}}_{p_{0}}^{\partial}$ is illustrated in Figure 10. Then surface gradient $\nabla_{t}$ can be represented by $\overline{\mathbf{d}}^{(0,\partial)}$, and $[\overline{\mathbf{d}}^{(0,\partial)}(\star^{(0,\partial)})^{-1}(\overline{\mathbf{d}}^{(0,\partial)})^{T}\star^{(1,\partial)}]\bm{E}$ is a primal $1$-cochain. However, since $\hat{n}\times(\nabla\times{{\mathbf{E}}})$ (tangential component of ${{\mathbf{H}}}$ field), the left hand side of Equation (86), normally is denoted by a dual $1$-cochain, Hodge star $\star^{(1,\partial)}$ needs to be implemented again in this third term.

Therefore, Equation (19) with second-order ABC is written as

[TABLE]

where we assume that the curvature is a constant on boundary surface for simplicity.

5 Numerical Examples

5.1 Validation: homogeneous circular waveguide

Then we can solve for the TM and TE modes in a hollow circular waveguide by using (62) and (68). The analytical value of $k_{s}$ for the TM and TE modes are roots of Bessel functions and roots of derivatives of Bessel functions. More specifically, $k_{s}=2.40482555769577$ for TM01 mode, and $k_{s}=1.84118378134065$ for TE11 mode. Comparing mesh with different fineness, the relation between relative error and maximum edge length $\Delta$ can be plotted as in Figure 11. A fitting shows that the convergence order is $2.0932$ for TM01 mode and $2.0689$ for TE11 mode. 666This second order convergence will be proved in our future publication.

The profiles of $E_{z}({{\mathbf{r}}})$ field for first six TM modes in a circular waveguide are plotted in Figure 12.

5.2 Microstructured optical fibers

Using Equation (65), we investigate the fundamental effective index of a step-index optical fiber and a air-hole assisted optical fiber (AHAOF), shown in Figure 13. The results are compared to both analytical solution (step-index fiber) and numerical solution by finite difference method [25].

The step-index optical fiber has parameters shown in Figure 13(a). It is surrounded with air (refractive index $=1$). The fundamental mode index is defined as

[TABLE]

Then, the fundamental mode index can be calculated analytically as $n_{\text{eff}}=1.438604$. This step-index optical fiber can also be viewed as an inhomogeneous waveguide, and its modes can be solved by Equation (65) with PEC boundary condition. Our numerical solution gives the fundamental mode index $n_{\text{eff}}=1.4386043519$ with 2,972 triangles and 1,561 vertices, which agrees very well with the analytical value. The transverse electric field intensity of first $9$ non-degenerate modes are plotted in Figure 14. There are two main methods to map a $1$-cochain to a vector field. One is to assume the vector field is constant in each triangle patch [26]; the other is to expand the vector field with Whitney forms [5, 6, 27].

The structure of AHAOF we have considered is shown in Figure 13(b). The guiding core is surrounded by air-holes. The advantage of this structure compared to normal optical fibers is that its dispersion is easily tailorable. Literature [25] used a $120\times 120$ grid for a quadrant window with 28,800 unknowns and solved the fundamental mode index as $n_{\text{eff}}=1.4353602$. We adopted a triangular mesh for the entire squire domain with 4,240 triangles and 6,400 edges (number of unknowns), and obtained $n_{\text{eff}}=1.4353696$. The intensity of $|E_{x}|$ and $|E_{y}|$ are plotted in Figure 15, and the profiles reflect the position of surrounding air-holes.

5.3 Photonic crystals

Photonic crystal, a periodic optical nanostructure, has very broad applications. Among these, two dimensional photonic crystals is not only used to produce commercial photonic-crystal fibers, they are also applied to form nano-cavities in quantum optics [28].

The band structure of photonic crystals can also be investigated in our framework with periodic boundary condition. By solving an eigen problem in Equations (73) and (74) with different $k_{x}$ and $k_{y}$, we can obtain a band diagram of a periodic structure. We considered two structures, shown in Figure 16. The structure in Figure 16(a) is dielectric cylinders positioned in squared lattice, and the structure in Figure 16(b) is a dielectric slab with air-holes placed on a triangular lattice. The band structure shown in Figure 16 agrees really well with results in reference [23] which uses FEM. Observation shows that the left structure admits a TM photonic band gap, while the right structure has a TE photonic band gap.

5.4 2-D scattering problems

Here we investigate the problem of open region scattering by a $2$-D perfect electrical conductor. We have already formulated the Dirichlet boundary condition at the conductor surface with an known incident field ${{\mathbf{E}}}_{\text{inc}}=\hat{z}E_{0}e^{ik_{0}x}$ as in Equation (64). We can adopt first-order or second-order ABC, as in Equations (77) and (82), for the outer truncation boundary. The structure we consider is NACA0012 airfoil, shown in Figure 17.

Then the scattered field and total field with second-order ABC can be plotted as in Figure 18.

5.5 Resonant cavities

Modal analysis can be used to determine the natural resonant frequencies and mode shapes of a structure in a broad field, such as structural mechanics, acoustics and electromagnetics [29].

By solving the eigen problem shown in Equation (67) with $\bm{J}=0$, we can obtain resonant frequencies of a resonator with conducting enclosures. We first examine an LES-mode resonator and a hybrid-mode resonator with structures introduced in [30] as shown in Figure 19. They are both inhomogeneous cavities. The comparison between current work and results from [30] used FEM is presented in Table 1.

With ABC boundary condition introduced, we can also calculate the resonant frequency and quality factor of a open structure. However, we need to reformulate the Equation (84) as

[TABLE]

We applied Equation (91) to find the resonant frequency and the $Q$ factor of the lowest TE mode for a dielectric sphere with radius $r=160$ $\mu$m and $\epsilon_{r}=36$. We placed this sphere in a $320\times 320\times 320$ $\mu\text{m}^{3}$ cube and discretized with $23,669$ tetrahedrons. The result is also summarized in Table 1.

6 Discussion and Conclusions

In this work, we have adopted discrete exterior calculus (DEC) to formulate and numerically solve various electromagnetic problems in frequency domain. In other words, we have provided an alternative method for computational electromagnetics analysis based on an arbitrary simplicial mesh.

Due to the nature of electromagnetics, the unknown fields are separated into primal cochains and dual cochains. But in practice, we always prefer to solve for the primal cochains, because the error introduced by the boundary can be minimized and the results can be interpolated with Whitney forms. And this is the reason why we treat $H_{z}({{\mathbf{r}}})$ and ${{\mathbf{H}}}({{\mathbf{r}}})$ field as primal cochains to solve for TE modes and closed $3$-D problems.

Since DEC keeps the structure and terseness of differential form description of Maxwell’s equations, charge continuity relation is exactly preserved, which leads to a great potential in problems involving motions of charged particles [31]. Another important feature is that all operators acting on cochains are naturally symmetric due to the diagonal Hodge stars.

In fact, since DEC is a tool to solve all kinds of partial differential equations, this method can also be applied to solve equations in many other fields, such as Navier-Stokes equations in fluid dynamics [32], Boltzmann equation in statistical mechanics, and Schrödinger equation in quantum mechanics.