Discrete exterior calculus for photonic crystal waveguides

The discrete exterior calculus (DEC) is very promising, though not yet widely used, discretization method for photonic crystal (PC) waveguides. It can be seen as a generalization of the finite difference time domain (FDTD) method. The DEC enables efficient time evolution by construction and fits well for nonhomogeneous computational domains and obstacles of curved surfaces. These properties are typically present in applications of PC waveguides that are constructed as periodic structures of inhomogeneities in a computational domain. We present a two‐dimensional DEC discretization for PC waveguides and demonstrate it with a selection of numerical experiments typical in the application area. We also make a numerical comparison of the method with the FDTD method that is a mainstream method for simulating PC structures. Numerical results demonstrate the advantages of the DEC method.

domain computations give information about the PBG 2,3 and with time domain simulations it is possible to analyze the time evolution of light waves before constructing physical PC products. 4 In order to be of practical use, the computational approaches are required to give accurate results in a reasonable computing time. The finite difference time domain (FDTD) method has become the standard tool to run computer simulations in the field of study. However, the complex geometry of the PC models is challenging to simulate accurately and efficiently with the conventional FDTD method. That is why the boundary element method (BEM) 5 or hybrid methods relying on the geometric flexibility of the finite element methods (FEM) have been applied. [6][7][8] The main drawback of the FEM-based simulation is that in time domain it typically requires the inversion or factorization of a sparse matrix at each time step. 9 The discrete exterior calculus (DEC) can be seen as a differential form based generalization of finite difference methods, including the FDTD. The DEC enables such grid structures that can handle complex domains. Since only inversion of a diagonal matrix is required at each timestep, it also provides efficient time-stepping. 10 The method is pioneered in electromagnetic simulations by Bossavit and Kettunen,11 while the groundwork presented by Marsden's group was applied in computer vision and image processing 12 and delivered the PyDEC software library. 13 During the last decade, the DEC has been applied to, for example, elastostatics, 14 electromagnetics, 15 quantum mechanics, 16,17 and flow problems. 18 Chen and Chew 2 applied the method for PC band structure computations in frequency domain. They also have shown that the method provides self-consistent and self-contained discrete electromagnetic theory. 19 Recently, Kettunen et al. 20 have advanced the method with higher order Whitney forms and Toshniwal and Hughes 21 have introduced isogeometric discrete differential forms. Despite the wide number of considered application areas and numerical experiments there exist only a few examples of wave propagation in heterogeneous domain in which the time-dependent problem is discretized with the DEC. These include electromagnetic scattering by objects such as hexagonal ice crystals or spherical particles randomly distributed into the surrounding media. 10,22 In the ice crystal simulations, the object has been considered to be large compared to the wavelength, while spherical particles were small compared to the wavelength.
In this article, we simulate a PC structure of an array of cross-sections of dielectric rods in the air substrate. Although our examples concentrate on demonstrating the feasibility of the DEC as a discretization scheme for PC waveguides, it is worth mentioning a few application areas benefiting from such computational tools. The models considered in our numerical studies form a basis in the design of optical logic gates and photonic circuits. 23,24 More sophisticated, hydrogel-based, PCs are important applications in the development of noninvasive optical glucose sensors. 25 PCs also provide future prospects in photovoltaic solar cell technologies. 26 The rest of this article is organized as follows. In Section 2, we present the mathematical model, first in the vector field notation and then in the language of differential forms. In Section 3, we discretize the differential form presentation in space by the DEC and in time by the conventional leap-frog scheme and consider the boundary conditions at the discrete stage. The numerical tests on transmission conditions, material properties, geometry approximation, mesh refinements strategies, and band gaps of PCs are based on software implemented at the University of Jyväskylä and considered in Section 4. The conclusions are presented in Section 5.

THE MAXWELL EQUATIONS
Time-dependent propagation of electromagnetic waves is presented by the hyperbolic system of the Maxwell equations where E = (E 1 , E 2 , E 3 ) T and H = (H 1 , H 2 , H 3 ) T are the electric and magnetic field strengths, is the electric permittivity, is the magnetic permeability, and J E = (J E1 , J E2 , J E3 ) T is the current density of the electric charge which is assumed to be divergence free, that is, ∇ ⋅ J E = 0. Respectively, we define the magnetic current density J H = (J H1 , J H2 , J H3 ) T . Electric and magnetic charges are presented by E and H , respectively. The electric flux density D, electric field strength E, magnetic flux density B, and magnetic field strength H are related by the constitutive equations D = E and B = H. The electromagnetic wave propagates with the angular frequency at speed c(x) = 1∕ √ (x) (x). The wavenumber (x) = ∕c(x) describes how many waves there are for a 2 unit. The corresponding wavelength is given by (x) = 2 ∕ (x) = c(x)∕f , where f = ∕2 is the frequency. The corresponding time-harmonic problem is considered with To truncate the domain Ω by an absorbing boundary Γ ext , we can use the Silver-Müller boundary condition where n is the outward pointing normal vector, to complete the circulation along missing dual edges of a boundary vertex. Another choice to restrict ourselves to a finite computational domain is to use a perfectly matched layer (PML). 27 To simulate dielectric rods in air, we concentrate on transverse magnetic (TM) modes. That is, Under the circumstances, the differential form representation corresponding to problem (1)-(4) is dB =̃H, where d is the exterior derivative,Ẽ = E 3 is a differential 0-form,H = H 1 dx 1 + H 2 dx 2 ,B = B 1 dx 2 − B 2 dx 1 , andJ H = J H1 dx 2 − J H2 dx 1 are differential 1-forms, andD = D 3 dx 1 ∧ dx 2 ,J E = J E3 dx 1 ∧ dx 2 , and̃H = H dx 1 ∧ dx 2 are differential 2-forms, and ∧ is the exterior product (wedge product) operator. The constitutive relations areD = ⋆Ẽ andB = ⋆H, where ⋆ and ⋆ are permittivity and permeability related Hodge operators mapping a differential p-form to a differential n − p form, where n is the spatial dimension. Respectively, we would have a 0-formH, 1-formsẼ,D, andJ E , and 2-forms B and̃E for considering transverse electric (TE) modes. 28 By replacing the time derivative operator t by a multiplier −i , we get the corresponding time-harmonic problem. 19

DISCRETIZATION
To discretize the computational domain Ω we construct a computational mesh. It is a collection of N Ω polygonal surface elements Ω k , k = 1, … , N Ω , such that Ω= ⋃ N Ω k=1 Ω k . There are N e edges e j , j = 1, … , N e , and N v vertices v i , i = 1, … , N v , in the mesh. In addition to this (primal) mesh we construct its dual mesh by connecting the nearest circumcenters (see Figure 1). The elements of both the primal and dual mesh can be seen as surfaces bounded by edges bounded by vertices. Thus, we can see a hierarchy of cells such that a vertex is a 0-cell, an edge between two 0-cells is a 1-cell, and a surface surrounded by edges is a 2-cell. For each p-cell in the primal mesh there is a (2 − p)-cell in the dual mesh. That is, the surface elements Ω * i of the dual mesh are dual to the primal vertices v i and the edges e * j of the dual mesh are dual to the primal edges e j . The dual mesh of a primal mesh constructed of squares is constructed of squares and the dual mesh of a primal mesh constructed of triangles is constructed of hexagons (see Figure 1). The spatial mesh step size, or edge length, is denoted by h. In a nonuniform mesh, we mark the shortest edge length of the primal mesh as h min and the largest edge length of the primal mesh as h max .

Discrete exterior calculus
The discrete counterpart of a differential p-form is a discrete differential p-form, or p-cochain. 29 The cochains defined on the primal mesh are primal cochains, and the cochains defined on the dual mesh are dual cochains. The primal 0-cochain E, associated with the primal vertices, and the primal 1-cochains B and J H , associated with the primal edges, are presented as column vectors with components Respectively, the dual 1-cochain H, associated with the dual edges, and dual 2-cochains D and J E , associated with the dual surfaces, are presented as column vectors with components The discrete analog of the Hodge star operator presents the constitutive relations at the discrete level. At this point the choice of circumcentric duality between primal and dual mesh is crucial since it provides the diagonality of the discrete Hodge operator that is also a key element of the efficiency of the time-stepping. We present the discrete Hodge star operator, that may vary with the properties of the medium, as ⋆ and clarify the associated material properties with a subscript. The discrete Hodge star operator ⋆ maps a primal 0-cochain E to the dual 2-cochain D taking into account permittivity . Since D = ⋆ E should hold, we can construct its inverse, ⋆ −1 , as an where d is the differential vector element of surface area  normal to the ith dual mesh surface Ω * i , |Ω * i | is the surface area of the ith dual mesh surface, and x v i gives the coordinates of the vertex v i . Respectively, we set H = ⋆ B. The discrete Hodge star operator ⋆ maps a primal 1-cochain B to the dual 1-cochain H taking into account permeability , such that = −1 . It is an N e × N e diagonal matrix with the diagonal components where |e * j | is the length of the jth dual mesh edge and |e j | is the length of the jth primal mesh edge.
Discrete exterior derivative d 0 is a boundary operator acting on 0-cochains to produce 1-cochains. It is a sparse N e × N v matrix presenting the incidence number between each edge e j , j = 1, … , N e , and each vertex v i , i = 1, … , N v . The incidence number is 0 if v i is not a point of e j . If an edge e j has a vertex v i , the value at the jth row and the ith column of the incidence matrix is ±1, depending on the relative orientation of the edge, that is, whether vertex v i is the end or start point of edge e j . Respectively, the discrete exterior derivative d 1 is a N Ω × N e matrix and it operates on 1-cochains to produce 2-cochains. If the edge e j is a boundary of face Ω k , the value at the kth row and the jth column of the incidence matrix is ±1, depending on the relative orientation between the face and the edge element (see Figures 2 and 3). The other values of the discrete exterior derivative, for the cases edge e j is not a boundary of face Ω k , are zeros.
We can summarize the discrete calculus presented above as a diagram F I G U R E 2 An example orientation of the edges e i , i = 1, … , 12, and the corresponding discrete exterior derivative matrices d 0 and d 1 for a square primal grid with vertices v j , j = 1, … , 9. The counterclockwise orientation of faces Ω k , k = 1, … , 4, is assumed to be positive.

F I G U R E 3
An example orientation of the edges e i , i = 1, … , 19, and the corresponding discrete exterior derivative matrices d 0 and d 1 for a triangular primal grid with vertices v j , j = 1, … , 10. The counterclockwise orientation of faces Ω k , k = 1, … , 10, is assumed to be positive. and present the spatially discretized counterpart of system (5)-(8) as

Leap-frog time-stepping
For time discretization we apply the staggered leap-frog type discretization familiar from the conventional FDTD method.
The one-dimensional time domain from t = t 0 to t = t N is divided into time step intervals, each of length Δt. We evaluate Equations (13) and (15) at time t = t n+ 1 2 and Equations (14) and (16) at time t = t n and approximate the time derivatives by the second order central finite difference method, applied to variables B and E staggered in time, such that That is, B is updated at half-integer timesteps, while E is updated at integer timesteps. The initial conditions are set as B − 1 2 and E 0 . To update the primal 1-cochain B, representing the magnetic flux through the mesh edges, and the primal 0-cochain E, representing the electric field located on the mesh vertices, we use the formulas The coupled system of equations (18) provides time-stepping that is explicit and conditionally stable. To consider the stability of the time evolution, we present the system of equations (18) in a modified recursive form where I N e on the upper left block of the coefficient matrix is an N e × N e identity matrix, whereas I N v on the lower right block of the coefficient matrix is an N v × N v identity matrix. By substituting, time step by time step, all the formulas of the earlier time steps to the formula with n = N − 1, we find that the closed form representation to the solution, at the end of the considered time period, is By assuming J n H = J For a conditionally stable scheme, we need to find such choices of parameters that the numerical solution remains bounded. A necessary condition for stability is that the absolute values of the all eigenvalues of G must be less or equal to one and G must have a complete set of distinct eigenvalues and eigenvectors. 30 Instead of the eigenvalues for the matrix G we search for eigenvalueŝ= − 1 for matrixĜ = G − I (N e +N v ) . We consider the eigenvalue equationĜu =̂u, where u is the eigenvector, multiply the first row from left by ⋆ −1 d T 0 ⋆ and solve it with respect to the first component of u, and substitute the result to the second row of the equation to get the characteristic equation̂2 Since ⋆ −1 d T 0 ⋆ d 0 is a diagonalizable matrix, multiplying it from right by the matrix of eigenvectors P and from left by the inverse of the matrix of eigenvectors P −1 results to the diagonal matrix with the eigenvalues of ⋆ −1 d T 0 ⋆ d 0 . Consequently, by multiplying Equation (23) from right by P and from left by P −1 results to the characteristic equation where are the eigenvalues of ⋆ −1 d T 0 ⋆ d 0 . Thus, the eigenvalues for G are If Δt 2 4 > 1, and the largest absolute value of the eigenvalues is larger than one. For Δt 2 4 = 1, there exists nondistinct eigenvalues = 1 − Δt 2 ∕2 and the stability is not satisfied. Thus, a necessary stability condition is Δt 2 ∕4 < 1, that is, Δt < 2∕ √ . In the case of square primal mesh elements, presented as an example in Figure 2, symbolic formulas for the eigenvalues are obtained by an eigenvalue algorithm for heptadiagonal matrices. 31

Computational domain and boundary conditions
PCs are periodically structured electromagnetic media. In certain cases, PCs may exhibit photonic band gaps (PBG), which are ranges of frequency in which light cannot propagate through the structure. To consider the phenomenon, we model a 2D square symmetrical PC, constructed of the cross-section of dielectric rods in the air substrate, in a rectangular domain illustrated in Figure 4. The radius of each dielectric rod is r and the lattice constant a = 1 is the distance between the centers of two nearest dielectric rods. Due to symmetry, we can restrict ourselves to simulate a rectangular domain [0, 10] × [0, 1], that is, a single row of the crystal, and terminate the upper and lower boundary with a perfect magnetic conductor (PMC) or periodic boundary condition. The update formulas (18) naturally implement a PMC boundary. 2 To enforce the periodic boundary condition we could use a cylindrical shape of the domain and a mesh with a suitable metric to derive the discrete Hodge star operators. However, here we get to the same result by modifying the update formulas derived for the rectangular domain. The computational domain is discretized, also with a triangular primal mesh, such that each point on the lower boundary corresponds to a point on the upper boundary with the same where i is the index of a boundary vertex and |e l | + |e m | is the sum of the lengths of the boundary edges linked to the ith boundary vertex v i (see Figure 5). All the terms ⋆ ii have the same sign because the simplices have the same orientation.

NUMERICAL RESULTS
In this section, we present numerical tests on transmission conditions, material properties, geometry approximation, mesh refinements strategies, and band gaps of PCs and compare the triangular DEC discretization with the square element based FDTD discretization. We have used the DEC code implemented in C++ at the University of Jyväskylä. 15 In principle, the solver is based on generalized finite differences, covering as a special case also the FDTD method. That is, it is possible to compare the DEC discretization with a computational mesh with triangular face elements with the FDTD as a special case constructed with a computational mesh with rectangular face elements. However, we used for the FDTD tests a program implemented in Matlab by the second author.
For all the simulations we set = 1.0 and modeled the material properties by varying the values of . The time domain simulations are initialized at t = 0 and finished at t = 200. With the FDTD we have used time step length Δt = 9h∕10 √ 2 and with the DEC Δt = h max ∕4, where h max is the largest edge length of the primal mesh. This implies that in the FDTD discretizations Δt∕h ≈ 0.64 and in the DEC discretizations Δt∕h min varies between 0.55 and 0.88 and the stability criteria are satisfied. The band structure graphs are obtained as eigenvalue problem solutions in frequency domain by following Chen and Chew. 2 We start with naive examples to demonstrate the simulation challenges and proceed to conquer them by adjusting material properties and optimizing and refining the mesh.

Transmission
Here we examine a method for finding the transmission coefficients for different wavelengths of TM-polarized light through a 2D square symmetrical PC. We send a incident wave from all of the left side boundary vertices by adding to their values at time instance t the modified Gaussian function where d = 3∕ f m defines the half-bandwidth of the pulse and f m is the center of the frequency area of interest. 32 The wavelength inside the surrounding medium is one half of the lattice constant a. Furthermore, the time shift t 0 = 3d reduces the effect of the initial discontinuity when the source is turned on at t = 0. Equation (28) is discretized as g n = g(nΔt), leading to an additional update formula E n | lhs = E n | lhs + g n on the left-hand side boundary vertices. 33 The transmission coefficient is computed by repeating the simulation twice, one simulation measuring transmitted waves with the dielectric rods and one as incident wave in the homogeneous domain without dielectric rods. 34 During the simulations, we obtain in time domain the field histories at a reference point on the right-hand side of the most right-hand rod as E n rods and E n no rods , n = 0, … , N − 1, respectively. Further, the time domain values are transformed to the frequency domain values E rods ( ) and E no rods ( ), respectively, with the discrete Fourier transform to get the transmission coefficient in decibels for frequency as 20 ⋅ log 10 ( |E rods ( )| In the first example, we discretize the computational domain by a triangle mesh with 10 edges (11 mesh points) per unit length on the boundary (see Figure 6). We set = 1, radius of each dielectric rod r = 0.3 , and permittivity of rods max = 5.9, while for air = 1.0. Simulating the setup with f m = 0.45, d = 6a∕ and = 6∕ leads to graphs presented in Figure 7 showing clear gaps in the transmission coefficient around frequencies 0.3 and 0.6.
The results presented in Figure 8 are in good agreement with the band structure considerations presented by Chen and Chew. 2 However, in our simulations, the PBG seems to extend below the frequency 0.3 unlike in the band structure diagram. This is likely because these modes in the band structure diagram correspond to diagonally traveling waves which our simulations do not address. Modeling of the circular cross-sections of the rods together with the coarseness of the mesh is a source of inaccuracies.

Adjusted material and geometry approximation for the dielectric rods
Since geometric modeling of the circular cross-sections of the rods is an important source of error, we observe this error source by numerical tests and propose remedies by adjusting both the mesh and the material properties. With both the conventional FDTD and the DEC it is common to model the rods or other inhomogeneities by adjusting the permittivity of a mesh point based on whether the point is within the radius of the cross-section of the rod or not. We tested how a strict rectangular approximation limits the choice of mesh resolution, that is, the number of mesh points per unit length, by setting the radius of cross-sections of the rods r = 0.2 and electric permittivity of the rods max = 8.9 and using several mesh resolutions. The effective area of the cross-section of the dielectric rod was set as where the sum goes over all the vertices inside a periodic cell, i is the permittivity assigned to the vertex i, and |Ω * i | is the area of the dual face of the vertex i. Effective radius r eff is the radius of a circle with an area of A eff . The wave field strength is measured at reference point (9.0, 0.5).
In Figure 9, we see that effective radius r eff = 0.2 provides the transmission curve (left bottom figure) without significant loss of accuracy for all tested mesh resolutions of 22, 32, 34, and 39 points per unit length (marked by red stars in the top figure). The deviation in the effective radius of the cross-section of the rod is shown to decrease the accuracy of the transmission simulation, as is presented for the mesh resolutions of 29, 33, 37, and 40 points per unit length (marked by blue stars in the top figure) in the right bottom figure in Figure 9. Thus, for the FDTD, that is, rectangular face elements, the geometry approximation of the rods is poor, unless a very fine mesh is used, and can be significantly improved by scaling the permittivity of a point based on how much of the corresponding dual area is inside the circle.
F I G U R E 9 Approximating the circular cross-sections of a dielectric rod with rectangles, the accuracy of the approximation and simulation results depend strongly on the number of mesh points per unit length.
In Figure 10, we demonstrate the effect of scaling the permittivity. The horizontal axis represents the modeling accuracy used to calculate the scaled permittivity for points near the boundary of a rod, and the value 1 corresponds to no scaling. The accuracy of the solution (red line in Figure 10) is directly related to the difference in effective radius (blue line in Figure 10) when the difference between radiuses is greater than 0.001. The radius of cross-sections of the rods was r = 0.2, so this corresponds to an error of 0.5%. For values smaller than this, it is likely still an important factor.
With the DEC, utilizing also triangular elements, we can model the mesh so that it approximates the cross-section of a dielectric rod more accurately. First we model the rod with a regular n-gon so that the effective area is exactly equal to the area of the rod cross-section. Since the edges of the n-gon need to be edges of the dual mesh, we draw appropriate sized circles to the n-gon vertices and construct tetragons based on the intersection points of these circles. The remaining areas are filled with a triangular mesh. The resulting mesh, for a low level discretization, of 11 points per unit length on the boundary is illustrated in Figure 11. In what follows, we call the mesh obtained by this construction as optimal or optimized mesh. This method can be straightforwardly generalized to accurately model arbitrary shaped material boundaries.
We consider the accuracy of the space discretization methods by comparing the transmission graphs, and measure the difference between the methods X and Y as the root mean square deviation (RMSD) calculated as where the sum goes over all the discrete frequencies of interest F, x f and y f are the transmissions in decibels corresponding to the frequency f , obtained by methods X and Y , respectively, and |F| is the number of discrete frequencies. The compared models are the FDTD with no permittivity scaling, the FDTD with permittivity scaling, and the DEC with an optimized mesh. Each of the compared simulations is setup so that they have exactly the same effective radius. For the first comparison we use the radius of the dielectric rod r = 0.2, electric permittivity max = 8.9, and 36 mesh points per unit length on the boundary. For the second comparison r = 0.3, max = 5.9, and there were 31 mesh points per unit length on the boundary. The frequency-transmission graphs through 6 layers and comparison between the models are shown in Figures 12 and 13, where the difference between the compared solutions is presented as a black graph, the transmission of the first mentioned method as a blue graph and the last mentioned method as a red graph. The deviation seems to be F I G U R E 10 Difference between transmission graphs compared to the difference between the effective radius of the cross-sections of dielectric rods. The compared simulations used a discretization of 23 and 24 mesh points per unit length. largest in areas where there is a sudden change in transmission. This is natural since these are also the most challenging parts to model accurately. The deviation seems to also increase with frequency, which is likely due to the decreasing wavelength and consecutively the number of mesh points per wavelength getting smaller. Some of the spikes in the error graphs could also be attributed to a slight frequency shift. The increase in the edge length, or mesh stepsize, raised all of the RMSDs as could be expected. In all the test cases, the obtained transmission graph was qualitatively valid. Permittivity scaling and mesh optimization have the benefit that they can be applied with high accuracy to any radius and discretization.

Mesh refinement
We study numerically the convergence rate and accuracy of the FDTD with parameter scaling and the DEC with an optimal mesh. The FDTD setup is truncated with a CFS-PML absorbing boundary layer 27 on the left and right side of the computational domain. The top and bottom sides are joined to effectively make the simulation area a cylinder. For the first test case, we study transmission through 6 layers, r = 0.2, and max = 8.9. The reference solution is calculated by linear extrapolation of three FDTD simulations using discretizations of 80, 90, and 100 points per unit length. The error for each simulation result is calculated as the RMSD to this reference solution. The results shown in Figure 14 show clear convergence in terms of mesh edge length. For the DEC, we have reported both the largest edge length (red stars in Figure 14, "DEC max dx") and the smallest edge length (black stars in Figure 14, "DEC min dx"). The convergence of the DEC in terms of the largest edge length has the same trend as the convergence of the FDTD (blue stars in Figure 14). Thus, the maximum edge length is the limiting factor for the accuracy and the use of a mesh that optimizes the material boundary modeling gains no significant improvement in the overall accuracy with this choice of parameters, geometry, and mesh refinement. Repeating the experiment for r = 0.3 and max = 5.9 shows (see Figure 15) a slight increase in accuracy for the DEC. By comparing the transmission graphs in Figures 12 and 13, we also observe that the transmission graph with this test's parameters is a bit smoother in the sense that it has gentler slopes.
As earlier mentioned, the error is related to the rate of change in the transmission graph and increases with frequency. Both of the methods tended to overestimate transmission in places where the slope of the transmission graph is positive and underestimate if the slope was negative, which could be the result of a slight frequency shift. This is demonstrated in Figure 16. Based on the above mentioned observations, we choose three points (marked as black stars  Figure 16) in the transmission graph and take a closer look of the transmission with respect to the edge length. In the bottom figure in Figure 16, we present the results with the FDTD with scaling (blue stars) and the DEC in terms of the largest edge length (red stars) and the reference solution (black star) for each of the chosen frequencies.
Since the compared methods had different boundary conditions (Silver-Müller and PMC for DEC, CFS-PML and periodic for FDTD) we should investigate if the choice of boundary condition has an affect on the results. Using the same rectangular mesh for both methods and measuring the RMSD of the transmission graphs, we get the results presented in Figure 17. Comparing this to Figure 14 or Figure 15, we conclude that the deviation is nearly an order of magnitude smaller and conclude that the use of different boundary conditions is not a significant error factor in the comparisons presented above.

Photonic band gap
A photonic band gap is a range of frequencies in which there does not exist propagating waves for any wavenumber, but there are propagating waves above and below the band gap. To find the band gap frequencies we compute, for r = 0.2, band structure graphs with two different permittivities. We present the graphs along the points G, X, and M (see Figure 18) of the Brillouin zone 35 for the TM modes. We are considering a normally incident wave, and thus we can concentrate on the section [G, X] of the band structure. In Figures 19 and 20, we illustrate the band structures, side-by-side with the transmission graphs, with max = 8.9 and max = 5.9, respectively. We observe two band gaps bounded by the blue dashed lines in both figures. The photonic band gaps are clearly positioned in the frequencies related to the pits in the transmission graphs. By introducing a line defect to the PC, frequencies within the bandgap are strongly confined inside the line defect. In the last numerical tests, we use the DEC space discretization with r = 0.2, max = 8.9, h max = 1∕32, and h min ≈ 0.4h max in an area of 15 × 7 periodic cells with the rods in the middle row removed. We set a narrow width modified Gaussian pulse, first centered around the frequency 0.4 and in another test run centered around the frequency 0.5, as a source in the middle of the left-hand side boundary and simulate with respect to time with Δt ≈ 0.008. As expected from the band structure and transmission graph 19, only the pulse corresponding to the frequency 0.4 is confined within the line defect (see Figure 21).

F I G U R E 18
The band structure is calculated by solving an eigenvalue problem at the points of the irreducible Brillouin zone highlighted in red F I G U R E 19 Comparison of the band structure and transmission graphs for r = 0.2, max = 8.9

F I G U R E 20
Comparison of the band structure and transmission graphs r = 0.3, max = 5.9 F I G U R E 21 Pulses centered around frequencies 0.4 (left) and 0.5 (right) in a photonic crystal array

CONCLUSIONS
We presented two-dimensional DEC discretization for PCs and demonstrated it with numerical experiments on material properties, geometry approximation, and mesh refinement strategies and demonstrated how the band gap frequencies are related to the transmission conditions. Modeling the circular cross-sections of the dielectric rods accurately and running the simulations computationally efficiently is a challenging task especially if only rectangular mesh elements can be used. However, the accuracy can be increased by scaling the permittivity on the rod boundary and adjusting the radius to better match the area. With the DEC, utilizing more generalized or mixed types of polygonal elements, there is also the option of designing the mesh to better match the nonpolygonal material boundaries. Hence, the DEC is well-suited for PC waveguide problems and with the diagonal discrete Hodge star operator it is computationally efficient especially in time-domain simulations.

ACKNOWLEDGMENTS
The authors would like to thank Jukka Räbinä for the discrete exterior calculus implementation and Tytti Saksa for discussions on discrete schemes.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.