## Abstract

The optimization of PhC waveguides is a key issue for successfully designing PhC devices. Since this design task is computationally expensive, efficient methods are demanded. The available codes for computing photonic bands are also applied to PhC waveguides. They are reliable but not very efficient, which is even more pronounced for dispersive material. We present a method based on higher order finite elements with curved cells, which allows to solve for the band structure taking directly into account the dispersiveness of the materials. This is accomplished by reformulating the wave equations as a linear eigenproblem in the complex wave-vectors *k*. For this method, we demonstrate the high efficiency for the computation of guided PhC waveguide modes by a convergence analysis.

©2010 Optical Society of America

## 1. Introduction

Photonic crystals (PhC) are materials comprising a periodic modulation of the refractive index in the order of the wavelength of the light [1]. PhCs have attracted much attention due to their exceptional ability to engineer the properties of light propagation. Light is guided efficiently in waveguides which are formed by omitting one or a few rows of holes in an otherwise perfectly PhC of finite size. In the field of integrated optics, the ability to tailor the dispersion of guided PhC waveguide modes is of particular interest, due to the fact, that those modes may exhibit narrow stop bands [2] or regions of flat bands and are commonly referred to as slow light modes. A characteristic feature of slow light modes is the enhanced light intensity in the core of the PhC waveguide [3], which is beneficial for all devices based on intensity dependent phenomena (*e*. *g*. nonlinear optics) [4]. The tailoring of the dispersion of guided modes to obtain slow light modes with low group velocity dispersion by slightly modifying the PhC waveguide has been demonstrated by Li *et al*. [5] and Kubo *et al*. [6].

Such a systematic approach of designing the waveguide dispersion requires an extensive parameter optimization and hence efficient methods for the computation of guided modes are prerequisite. The efficiency issue is even more severe for active devices such as light sources or optical amplifiers. Those devices are typically operated close to the band gap energy of a direct semiconductor, where the complex permittivity is strongly frequency dependent. Most conventional methods solve for the eigenfrequencies *ω _{i}* and require intrinsically additional iterations of the permittivity

*ε*(

*ω*) to handle dispersive materials. In this article, we present a new method which naturally includes dispersive and lossy materials by the direct solution of the eigenvalue

_{i}*k*for a given

*ω*which exhibits numerically an exponential convergence. Hence the method is especially suited for optimization of PhC waveguides of both, dispersive and non-dispersive materials. Compared to the methods based on searching the eigen-frequencies

_{i}*ω*for a given wave-vector

_{i}*k*[7], we are able to directly access the lossy modes by computing the complex eigen wave-vectors

*k*. However, computing losses is not the object of this article.

_{i}In the presentation of our method we will emphasize three different aspects:

- the
*formulation*of the mathematical eigenvalue problem for*k*of a given_{i}*ω*,_{i} - the
*modelling*of guided PhC modes*in an infinite geometry*, and - the
*discretization*of the problem.

Our method differs fundamentally in the first and last aspect compared to the current methods. Before detailing the differences we will recall well-established methods for modelling of PhC wave-guide modes with an emphasis on these three categories.

**Formulation :** The standard procedure to obtain the photonic bands of a PhC is to solve an eigenvalue problem in the unit cell for a discrete set of wave-vectors *k _{i}* of the irreducible Brillouin zone. Subsequently, the lowest frequencies

*ω*(

_{i}*k*) for each wave-vector sample

_{i}*k*is computed for the unit cell (see [1, 8–11] and the references therein). We call the formulation underlying this procedure the

_{i}*ω*-formulation.

Typically, devices for tele-communication are operated only within particular frequency bands. Test structures are excited by a monochromatic light source with a frequency *ω*
_{0} and not by imposing a certain wave-vector *k*
_{0}. The question is hence ’how does the wave propagate when excited with light of frequency *ω*
_{0}?’. This real world situation motivates the *k*-formulations, where we solve for the eigen wave-vectors *k _{i}* at a given frequency

*π*

_{0}. In the

*k*-formulations the corresponding permittivity

*ε*(

*ω*

_{0}) can be used for each frequency

*ω*

_{0}turning the problem into a direct formulation for frequency-dependent materials. Smajic

*et al*. [12] already proposed in 2004 to tackle the problem by computing the complex wave-vector

*k*of a certain mode for a given frequency

*ω*

_{0}. However, they studied material with small absorption and used an iterative approach, which requires scanning of a large area of the complex wave-vector plane to find the solutions.

**Modelling in an infinite geometry :** Guided modes in the real PhC waveguide are localized in the ’hole-free’ dielectric channel and decay both in the crystal and the surrounding medium which motivates the following two approaches.

In PhC waveguides the dispersion of the guided modes are commonly computed by using the super-cell approach [1, 13, 14], where the standard band structure codes for PhCs are applied to a super-cell instead to the unit cell (cf. Fig. 2(b), 2(c) for examples). The super-cell system models an infinite array of identical line defects separated by a certain number photonic crystal rows. An alternative approach is to replace the finite PhC with a lateral periodicity by two semi-infinite PhCs. For the resulting infinite system, the existence and properties of guided modes have already been studied for the TE [15] and the TM case [16].

Both approaches cause a modelling error as they differ from the exact geometry in which a homogeneous dielectric region is attached to the crystal [see Fig. 1(a)]. Since a major aspect of this article is the comparability of the methods in terms of efficiency and accuracy for dispersive media rather than the investigation of the boundary conditions, we constrain to the super-cell approach.

**Discretization :** Many different methods and codes are proposed to compute photonic bands. The most prominent code is MIT Photonic Bands (MPB), which has been developed by S. Johnson *et al*. [17]. MPB is a method based on a global plane wave expansion (PWM) for periodic systems. Whereas the method gives exponential convergence for smooth potentials *e*. *g*. present in quantum-mechanics, it approximates the eigenmodes in PhCs only with a linear rate of convergence due to their discontinuous dielectric spatial distribution. On the other hand, the finite difference time domain method (FDTD) is ubiquitously used for electromagnetic scattering problems. Most FDTD codes [18] allow to solve for eigenmodes of periodic systems. A complete set of eigenmodes can be obtained with a single computation for a range of frequencies if the location and pulse of the excitation is properly chosen [19].

Both, FDTD and PWM are based on a uniform rectangular grid, which makes it difficult to resolve curved features and features not aligned to the grid. Opposed to FDTD and PWM, the finite element method (FEM) uses polynomial basis functions on every cell of an *irregular* – triangular or quadrilateral – mesh. Furthermore, the FEMs can be classified by their refinement strategy [20, 21]. For the *h*-version the mesh is refined for a fixed polynomial degree. On the other hand, the *p*-version keeps the mesh while enlarging the polynomial degree. The *hp*-version is a combination of both. Similar to the PWM one gets a matrix eigenvalue problem which can be solved by standard linear algebra packages for sparse matrices like Arpack. For the computation of the PhC band structure it was shown that the *h*-version of the FEM with polynomial degree 2 already achieves double the convergence rate of MPB [22]. For material interfaces having continuous contours (continuous also in their derivatives such as *e*. *g*. circles) the *p*-version of FEM with curved cells leads to an *exponential convergence*, *i*. *e*. the convergence rate increases continuously. And a limited, but high convergence rate (algebraic convergence) can be achieved for the “worst case” scenario, *i*. *e*. geometries with contours with sharp corners [11].

We pursue the following approach:

**Formulation:***Quadratic eigenvalue problem in k for each ω*._{i}Using the Bloch transform we obtain a formulation on the unit-cell with periodic boundary conditions, which exhibits a quadratic dependence on

*k*and a linear dependence on in*ω*^{2}. This formulation can be interpreted as*ω*- as well as*k*-formulation. In the following, we focus on the*k*-formulation. The quadratic equation in*k*can be linearized, such that efficient and well understood algorithms for solving eigenproblems with linear sparse matrices can be applied. The same approach was recently presented for the band structure calculation of TM modes in infinite crystals in [23, 24] and in [25] for frequency-dependent materials which we applied in this project to PhC waveguides.**Modelling in an infinite geometry:***Super-cell approach*.For comparison we pursue the most often used approach – the super-cell approach – with periodic boundary conditions at all boundaries. This approach does not influence the formulation to be quadratic in wave-vector

*k*. But it introduces a small modelling error which will be studied in a forthcoming article in comparison with two other approaches which also lead to quadratic eigenproblems in*k*.**Discretization:***High-order FEM on a coarse mesh with curved cells*.To be able to approximate the waveguide modes to a high precision with moderate computational costs we use the

*p*-version of the FEM on a coarse quadrilateral mesh with curved cells [11, 25]. For this system we implemented both, the*k*-formulation and the*ω*-formulation for strictly periodic boundary conditions using the numerical C++ library*Concepts*[26]. For practical geometries with smooth interfaces we will show the expected exponential convergence of the method and will compare it to the results of MPB and the FEM software COMSOL.

The article is organized by first formulating the eigenvalue problem in the super-cell for the TE and TM modes (Section 2). Section 3 is dedicated to the numerical solution of the eigenvalue problem which includes the discretization by the *p*-FEM and the new *k*-formulation. Finally, in Section 4 we compare the convergence rate of the presented method for the conventional *ω*-formulation to those of other methods and to that of the *k*-formulation. Finally we study the modes in a PhC waveguide with dispersive materials which shows the importance of considering the frequency-dependence.

## 2. The eigenvalue problem in the super-cell

The three-dimensional time-harmonic Maxwell’s equations (time dependence as e-^{iωt}) for a linear, isotropic, non-magnetic and dispersive material is written as:

in which the electric and magnetic fields are denoted by **E**(**x**) and **H**(**x**), the vacuum permeability and permittivity are *μ*
_{0} and *ε*
_{0} and *ε*(**x**, *ω*) stands for the frequency-dependent complex relative permittivity.

We consider 2D PhC waveguides as depicted in Fig. 2(a), 2(b) with one-dimensional periodicity *a* in propagation direction **e**
_{1} and constant permittivity in one direction let it w.l.o.g. be **e**
_{3}. The guiding core layer of the waveguide is laterally surrounded by a finite length *L* of PhC material followed by a homogeneous dielectric material which extends to infinity. The unit-cells for this configuration are infinite strips.

However, we will follow the commonly used super-cell approach, where the described geometry is replaced by a lateral finite crystal with additional periodicity conditions connecting top and bottom boundaries. One has to be aware of the fact that by this simplification a modelling error is introduced, which we ignore in this article.

We will use the Bloch transform 𝓣 [27] formally defined by

which is with the exception of the factor e^{-ikx1} equivalent to the Floquet transform 𝓕 [28, 29]. The parameter *k* ∈ *B* := [-*π*/*a*, *π*/*a* is the real wave-vector for loss-less materials, **u** is an arbitrary function and *B* the (one-dimensional) Brillouin zone.

After applying Eq. (2) to the time-harmonic Maxwell equations [Eq. (1)]

and the definition of the *shifted curl operator*

which is characterized by the property **curl**
_{k} 𝓣 ○ = 𝓣 **curl**○, we obtain

The transformed fields **E**̃ and **H**̃ are periodic in *a*
**e**
_{1}. The term *ik*
**e**
_{1} × **u**̃(*k*, **x**) in the shifted curl operator [Eq. (3)] results from applying **curl** to e^{-ikx1}. Here the derivatives in **e**
_{3} are discarded due to the assumption of a 2D configuration.

The first two components of the operator **curl**
_{k}
**u** act only on the third component of *u* and the last component **curl**
_{k}
**u** acts only on the first two components. Hence, the modes decouple in the transverse electric (TE) mode with *H*̃_{1}(**x**) ≡ *H*̃_{2}(**x**) ≡ *E*̃_{3}(**x**) ≡ 0 and the transverse magnetic (TM) mode with *E*̃_{1}(**x**) ≡ *E*̃_{2}(**x**) ≡ *H*̃_{3}(**x**) ≡ 0. Each mode is specified by a scalar problem in the out of (**e**
_{1}, **e**
_{2})-plane component and a vector valued problem for the two in-plane components. Both, the respective scalar and vector valued problem provide the same spectrum. Hence, we only consider the scalar valued problem of the out-of-plane field component which we rename according to *E*̃(**x**) := *E*̃_{3}(**x**) and *H*̃(**x**) := *H*̃_{3}(**x**) for simplicity. Therefore we need the last component of **curl**
_{k}
**curl**
_{k}
**u** for the TM mode and **curl**
_{k}
*ε*
^{-1}
**curl**
_{k}
**u** for the TE mode. It can easily be verified that they simplify to

The Bloch transform [Eq. (2)] guarantees the *a*
**e**
_{1}-periodicity of the fields *E*̃ and *H*̃ as well as their derivatives. Hence, we can formulate the problem by the following set of equations for the TM mode

and the TE mode

respectively. The two sets of equations [Eq. (5) and Eq. (6)] are stated on the two-dimensional strip Ω. [see Fig. 2(a) and 2(b)] and completed by following the super-cell approach by using periodicity conditions on the upper and lower boundaries. We write the general equation [Eq. (6c)], which results from the Bloch transform for the case of a super-cell with discontinuities of the permittivity at its borders.

## 3. The numerical solution of the eigenvalue problem

#### 3.1. The general eigenvalue problem

Now, we will state the general eigenvalue problems [Eq. (5) and Eq. (6)] for the TM and the TE mode in their weak form, *i*. *e*. we search for eigenmodes in *H*
^{1}
_{per}(Ω), the *H*
^{1}-Sobolev space of weak solutions with periodic boundary conditions on all sides of the super-cell Ω. The weak formulations are as follows:

Seek triples (*ω*,*k*,*E*̃),(*ω*,*k*,*H*̃) ∈ ∝^{+} × *B* × *H*
^{1}
_{per}(Ω) such that for all test functions *v* ∊ *H*
^{1}
_{per}(Ω)

By introducing the index ℓ ∊ {-1,0} and powers of the dielectricity *ε*
^{-1}(**x**,*ω*)= 1/*ε*(**x**, *ω*), *ε*
^{0}(**x**,*ω*) = 1 and *ε*
^{1}(**x**,*ω*) = *ε*(**x**,*ω*) we can combine the weak formulations for TM and TE mode into one equation:

$$\text{}+k\phantom{\rule{.2em}{0ex}}i\underset{{c}_{\ell}\left(u,v\right)}{\underbrace{{\int}_{\Omega}{\epsilon}^{\ell}(u\left(x\right).\frac{\partial \overline{v\left(\mathbf{x}\right)}}{\partial {x}_{1}}-\frac{\partial u\left(\mathbf{x}\right)}{\partial {x}_{1}}.\overline{v\left(\mathbf{x}\right)})d\mathbf{x}}}+{k}^{2}\underset{{b}_{\ell}\left(u,v\right)}{\underbrace{{\int}_{\Omega}{\epsilon}^{\ell}u\left(\mathbf{x}\right).v\left(\mathbf{x}\right)d\mathbf{x}}}=0.$$

For ℓ = 0, the equation provides the TM modes with *E*̃ = *u* and for ℓ = -1 it provides the TE modes with *H*̃ = *u*. With the above defined bilinear forms *a*
_{ℓ}, *b*
_{ℓ} and *c*
_{ℓ} we can write compactly

#### 3.2. Computational frameworks for the calculation of the modes

The obtained equations [Eq. (9)] for ℓ = -1,0 can be solved by two different methods:

- The search for the eigenfrequencies
*ω*for particular values of*k*which we call the_{i}*ω*-*formulation*. - The search for the eigen wave-vector
*k*for particular frequencies*ω*which we call the_{i}*k*-*formulation*.

The *ω*-formulation is quadratic in *ω* only for strictly non-dispersive permittivities for which the problem can be reduced to a linear eigenproblem in *λ* = *ω*
^{2}.

The *k*-formulation turns out to be quadratic in *k* also for dispersive materials, *i*. *e*. *ε*(**x**, *ω*) depends actually on *ω*. Hence, we benefit from two advantages of the *k*-formulations [Eq. (5) or Eq. (6)]: the sampling scheme over frequency values and the quadratic – and so simply to linearize – eigenvalue problem also for frequency-dependent material.

Note, that with a *k*-formulation derived from the Floquet transform instead from the Bloch transform, one would search directly for the quasi-periodic modes. However, due to the quasi-periodicity factor exp(*ikx*
_{1}) on the boundaries of the strip, the system can not be turned to a quadratic system in *k* anymore. Hence, the derivation of the *k*-formulation from the Bloch transform is essential.

#### 3.3. The matrix eigenvalue problem

Now, we search for eigenmodes *u*
_{FE} in some FE space (refer Sec. 3.4) with its basis functions *ϕ _{j}*,

*j*= 1,…,

*N*. This means, that we search for eigenmodes

*u*

_{FE}of the form ${u}_{\mathrm{FE}}={\mathbf{\Sigma}}_{j=1}^{N}{u}_{j}{\varphi}_{j}.$ Inserting this expression in the eigenproblem [Eq. (9)] with a frequency dependent permittivity

*ε*(

*ω*,

**x**) and each basis function

*ϕ*as test function we obtain the

_{j}*quadratic matrix eigenvalue problem*

Here *u*⃗ = (*u*
_{1}, *u*
_{2},…,*u _{N}*)

^{⊤}is the coefficient vector of the solution and the matrices are defined by

$${\mathbf{M}}_{2}\left(\omega \right):={\left({c}_{\ell}({\varphi}_{j},{\varphi}_{i})\right)}_{i,j=1}^{N},\phantom{\rule{1.2em}{0ex}}{\mathbf{M}}_{3}\left(\omega \right):={\left({b}_{\ell}({\varphi}_{j},{\varphi}_{i})\right)}_{i,j=1}^{N}.$$

With *x*⃗ := (*k*
*u*⃗,*u*⃗) we can simply linearize [30] equation [Eq. (10)] resulting in the *linear matrix eigenvalue problem*

The matrix **M**
_{3} (*ω*) and consequently **M**(*ω*) is symmetric and positive definite as the permittivity *ε*(**x**, *ω*) ≥ 1 and sup_{x∊Ω}
*ε*(**x**, *ω*) < ∞. The matrix eigenvalue problem [Eq. (11)] provides a complex spectrum *k*(*ω*). The eigenvalues *k* of the guided and loss-less modes are identified by searching the obtained spectrum for wave-vectors having a negligible imaginary part. Note, that **M**
^{2}, **M**
_{3} and consequently **M** are not frequency dependent for TM modes (ℓ = 0) and for TE modes in case of strictly frequency independent permittivities.

The linear matrix eigenvalue problem [Eq. (11)] was used to compute the guided modes for two frequencies of the following test configuration with a non-dispersive loss-less material which are shown in Fig. 3(a)–3(c). For this computation we used the *p*-version of the finite element method which we will explain in the sequel.

**The test PhC waveguide** We study the TE modes of a PhC W1 waveguide based on a hexagonal lattice structure with five lateral rows of air holes in an otherwise homogeneous and isotropic dielectric media with *ε* = 11.4.

#### 3.4. Discretization by the p-version of the FEM

In the previous section we have defined the matrix eigenvalue problem [Eq. (11)] for a FE space with basis functions *ϕ _{j}*. To define these basis functions, a computational mesh for the super-cell has to be created first. A proper mesh resolves the material interfaces as accurately as possible. The FE spaces for the continuous

*H*

^{1}-Sobolev space consist of basis functions being continuous across the boundaries of the cells of the mesh. The FEM of lowest order uses the so called “hut functions” as basis functions which take the value 1 at one node of the mesh and which are linearly decaying in the neighboring cells. For the FEM of higher orders, polynomials of an order larger than

*p*=1 are used. We distinguish between ”edge-functions”, defined on an edge, and ”bubble-functions” which are zero on the edges and on all other cells and hence, they account for the contribution of a single cell only. In this way, basis functions can be defined for triangular and quadrilateral cells for any polynomial degree

*p*, and very similarly also for curved cells. For further details see the monograph of Karniadakis and Sherwin [31].

When the meshes, which are composed of only straight cells, are refined, the (curved) material interfaces are better resolved. But when only the polynomial degree is increased the material interfaces should be sufficiently resolved by the mesh. By using circular cell boundaries in the meshes, we are able to perfectly model the circular interfaces in the geometries with round holes, *i*. *e*. no additional error is introduced by the mesh, also for the coarse grid.

The accuracy of the solution can be improved by refining the mesh (*h*-FEM), by increasing the polynomial order of the basis functions in the cells (*p*-FEM) or by a combination of both (*hp*-FEM).

The mentioned refinement methods have different efficiencies, *i*. *e*. the achieved accuracy of the solution using a particular number *N* of basis functions. The *h*-version of the FEM achieves an algebraic convergence like *N*
^{-α} where the convergence rate *α* is restricted by the polynomial degree *p* and by singular points in the material contours. By doubling the number of basis functions the error decreases by a factor 2^{α}, *e*. *g*. by 4 for a method of order 2.

The *P*-version of the FEM can approximate eigenmodes with exponential convergence like exp(-*βN*
^{½}) = *B*
^{√N} for some *β* > 0 and *B* := *exp*(-*β*). For non-smooth contours the *p*-FEM would achieve only algebraic convergence where, however, the convergence rate would be still larger than for *h*-FEM. For this “worst case scenario” a proper combination of mesh refinement and polynomial degree enhancement [11, 21] retrieves the exponential convergence.

**The implementation** The arguments above make it evident, why we have chosen the *p*-FEM as our method of choice, and therefore we use the *h*-FEM only for comparison.

For the computation of the guided modes we implemented the *ω*- and the *k*-formulation in the numerical C++ library *Concepts* [32] using a quadrilateral mesh with curved cells. A band calculation example is shown in Fig. 4(a) for the *ω*-formulation and Fig. 4(b) for the *k*-formulation. The used mesh is the one of Fig. 3(d) and some modes are depicted in Fig. 3(a)–3(c), which are obtained with the *k*-formulation using a polynomial degree of *p* = 8.

## 4. Numerical results

In this section we will study numerically the efficiency of the proposed method. Even though the focus of the article is on the *k*-formulation, we start with the convergence analysis of the *ω*-formulation for non-dispersive material (in Section 4.1), because this allows a direct comparison of the efficiency of our *hp*-FEM code to MPB and COMSOL. As the other two codes only support the *ω*-formulation, the comparison is restricted to the computation of eigenfrequencies. Then, in Section 4.2 we study the convergence behavior of the *k*-formulation with our *p*-FEM and the *h*-FEM implementation, which is very similar to the one of the *ω*-formulation. Finally, we investigate the influence of the material-dispersion on the PhC waveguide modes in Section 4.3.

*4.1. Comparison of the* p-*FEM with other methods for the ω-formulation*

The convergence behavior of the best-approximation to a solution of a partial differential equation by means of *p*-FEM discretization is proven in [21] and can be transfered to the FEM solution of many equations. This has recently been verified for the computation of the PhC band structure and the *ω*-formulation [11]. In an example, whose results are shown in Figure 5, we investigate the efficiency of the *p*-FEM within our library *Concepts* for computing the waveguide modes to the efficiency of MPB and the FEM software COMSOL as well as to our *h*-FEM.

We investigated the test PhC waveguide depicted in Section 3.3. The mesh for the super-cell used for *p*-FEM in *Concepts* is presented in Fig. 3(d). For the *h*-FEM experiments with *Concepts* we used the same mesh at the coarsest level which is then further refined. The computations with MPB and COMSOL have been performed with their respective meshes. The FEM software COMSOL uses triangular meshes and polynomial degree 2. For all the methods the same one-dimensional Brillouin zone is sampled at 51 equally spaced values. For every wave-vector *k _{i}* the 30 smallest eigenvalues

*ω*(

_{j}*k*) are computed. Subsequently, only the two guided modes in the PhC band gap are selected for the convergence analysis. These eigenvalues are compared to those of the reference solution (

_{i}*Concepts*with polynomial degree

*p*= 19) by evaluating the average error. This error is obtained by summing up the differences between the approximated and the exact eigenfrequency solutions, which is then divided by the total number of compared values. The error is plotted in Fig. 5 once versus the degree of freedom (

*left*) and once versus the computing time (

*right*). All simulations have been performed on a Debian GNU/Linux SMP system running on a Sun Fire V40z having 4 Dual Core AMD Opteron Processors at 2591 MHz and equipped with 32GB RAM. The specified time in Fig. 5 (

*right*) correspond to the time required by the software on a single CPU core to complete the program.

Figure 5 shows a quadratic convergence of COMSOL and *Concepts* with *h*-refinement (uniformly *p* = 2), and almost quadratic convergence for MPB up to an relative error of 10^{-3}. For further refinement MPB achieves approximately linear convergence. Even if the slope of the convergence of COMSOL and *Concepts* with *h*-refinement is similar, it is interesting, that the convergence curve of COMSOL has an offset of roughly one order fo magnitude to *Concepts*. We attribute that offset to a higher quality of the mesh used by COMSOL. However, for *p*-refinement the mesh quality does not contribute a lot and a high convergence rate of about 7.5 is achieved. Hence, it is the most efficient method of the presented ones for all relative error levels below 10^{-4}. With *p* = 15 and 9737 degrees of freedom the averaged relative error of the dispersion of the guided modes is about 10^{-8}. Note, that the exponential convergence of the *p*-refinement, with its property that the convergence rate should increase further, is difficult to distinguish from a high algebraic one (cf. [11]) at this refinement level. We refer to the example in next section, where the exponential convergence is better observable.

Similar behaviors of the convergence curves are observed for the error versus time: MPB is converging linearly, COMSOL and *Concepts* with *h*-refinement and *p* = 2 have roughly a quadratic convergence rate and *Concepts* with *p*-refinement converges with order 4. For the *p*-refinement the computational effort increases faster with increasing number of degrees of freedom *N* due to our current implementation of the curved cells. These results may be improved. Nevertheless, the *p*-FEM implementation in *Concepts* is superior to other methods, for all error levels below 10^{-3} for the presented example.

#### 4.2. Convergence analysis for the k-formulation

The second convergence analysis is performed with the *k*-formulation proposed in this work. A direct comparison neither to MPB nor to COMSOL is no longer adequate.

We show in Fig. 6 a similar convergence plot as in the previous section. The band structure of the same test PhC waveguide was computed using the *k*-formulation. Hence, the frequency axis *ω* was sampled at 51 equally spaced samples from 0 to $0.5=\frac{2\pi c}{a}$, and similar to the previous experiment, we computed the 30 eigenvalues of lowest magnitude were computed. The resulting complex eigenvalue spectrum contains many complex conjugate pairs of eigenvalues. The band diagram can be found by selecting the eigenvalues having an imaginary part close to zero (imag(*λ*) < 10^{-7} and simultaneously having a positive real part (real(*λ*) ≥ 0). Analogously to Fig. 5, the averaged sum of the difference between the approximated and the exact wave-vectors is plotted in Fig. 6. The convergence plots in Fig. 6 can not directly be compared to those in Fig. 5 since Fig. 5 shows the error in the frequency and Fig. 6 that in the wave-vector. Nevertheless, the posed problems are *comparable* in terms of computational resources and also in behavior of the slope of the convergence curve. We observe also quadratic convergence for *h*-refinement with *p* = 2 and a convergence rate of about 7.2 for the *p*-refinement which comes apparent to be exponential w.r.t. √*N*, where *N* again denotes the number of degrees of freedom.

The dash-dotted line above the *p*-refinement curve in Fig. 6 (*left*) is the exponential fit with parameters *β* = 0.16 and *B* = 0.85. For this example, the error decreases by a factor 0.85^{3} ≈ 0.61 if the number of degrees of freedom is increased from *N* to (*N*
^{1/2} + 3)^{2} ≈ *N*(1 + 6*N*
^{-1/2}). Thus, to decrease the error by this factor (*e*. *g*. 0.61 in our example), the corresponding relative increase of number of degrees of freedom (*N*
_{i+1} - *N _{i}*)/

*N*, becomes smaller for increasing

_{i}*N*and tends to zero. Note, that the convergence curve fits locally also to an algebraic law.

_{i}The results shown in this section verifies the superiority of the *p*-FEM for the computation of PhC waveguide modes for a single or a bunch of frequencies with the *k*-formulation. As already mentioned, the *k*-formulation is especially suited for dispersive material, as the permittivity *ε*(*ω*
_{0}) at the investigated frequency *ω*
_{0} can be taken directly. In the following section we present the numerical results for a dispersive InP-PhC waveguide.

#### 4.3. Computing Band Diagrams including Dispersive Material Models

The dispersion of semiconductors in the transparent region is commonly neglected in integrated optics, because of two reasons: first due to the small propagation distances of the light on-chip and secondly due to the relatively small dispersions for the frequency range (a few GHz) of a certain communication band. As long as a PhC device is operated in a narrow band within the transparency region of a material, one may use a constant permittivity for the center wavelength. However, in case of a photonic band gap calculation, we cannot neglect the dispersion anymore. Hence, we have to restrict the computation to a particular lattice constant *a* and to the frequency range of validity of the permittivity model. The first is a consequence of the fact that the scaling property of PhC [1] is lost in the case of a dispersive permittivity *ε*(*ω*). In the previous computations we used units, which are relative to the lattice constant *a*. Then, the photonic bands transform to any frequency by choosing the same pattern and filling factor, but a scaled lattice constant. For example the operation regime can be transformed to the doubled frequency *ω*
^{′} = 2*ω* by adjusting the lattice constant to *a*
^{′} = *a*/2. For dispersive materials the permittivity for the doubled frequency *ω*
^{′} is, not the same and we cannot simply transform the results.

**Test configuration** We study the TE modes of a PhC W1 waveguide based on a hexagonal lattice structure with lattice constant *a* = 400nm and five lateral rows of air holes in InP, the substrate material for our devices. Apart from the fixed lattice constant *a*, the geometrical settings are identic to the test PhC waveguide in Sec. 3.3. For the permittivity of InP, we use the model of S. Adachi [33], which accurately agrees with measurements by spectroscopic ellipsometry in the near infra-red regime. For wave-lengths above 920nm the model is given by:

where *α _{i}*,

*β*,

_{i}*C*and

*γ*denote positive material constants and

*δ*

_{i=1}stands for the Kronecker symbol which is 1 for

*i*= 1 and 0 otherwise. The energy levels

*E*

_{0},

*E*

_{0}+ ∆

_{0},

*E*

_{1}+ ∆

_{1}and

*E*

^{′}

_{0}correspond to different interband transitions. In this experiment we neglected the imaginary part of permittivity, thus

*ε*(

*ω*) is purely real. Note, that the complexity of the used permittivity model is irrelevant when using the

*k*-formulation, where we directly take the permittivity for the chosen frequency.

In Fig. 7 we compare the band diagram of the dispersive InP (red) to the one with a permittivity fixed to the frequency *ω*
_{0} (blue) for which $\frac{{\omega}_{0}a}{2\pi c}=0.26$ for *a* = 400nm holds. The corresponding wavelength is *λ*
_{0} = 1538nm. On the right side of the plot *ε*(*ω*), the used real part of the permittivity model with material constants according to [33], is plotted. From this illustration, we can observe a ‘down’ shift of the higher band gap frequency from 0.304 to 0.302 $\left[\frac{\omega a}{2\pi c}\right]$ and a small ‘up’ shift of the lower band gap frequency from 0.2252 to 0.2264 $\left[\frac{\omega a}{2\pi c}\right]$ – a reduction of the band gap of 4%. The deviations in this example are rather small. However, if small features in the dispersion are investigated, *e*. *g*. flat dispersion curves, resonances, etc …, a noticeable function deviation may result.

To further quantify the error of the wave-vector *k*, which is introduced by a deviation of the permittivity *ε*
^{′} = *ε*(*ω*
_{0})-∆*ε* = 0.26 from that for *ω*
_{0}, the permittivity is varied with ∆*ε* ∊ [0,0.5] while holding *ω* = *ω*
_{0}. A linear relation (∆*k* = 0.0486∆*ε*) between ∆*ε* and the deviation in the wave-vector ∆*k* is obtained and depicted in Fig. 8.

The issue of dispersive materials is more relevant for designing active PhC devices. In those designs, metallic layers for electrical contacts, active semiconductor layers for optical gain and highly doped semiconductor layers, which allow for carrier transport, are used. These layers all share a strong frequency dependence and a strong absorption. A detailed discussion, of strongly dispersive and also lossy materials in conjunction with PhC is out-of-scope of this article and will be published elsewhere.

## 5. Conclusion

We have presented a method for the direct computation of PhC waveguide modes for a given frequency *ω*
_{0} by solving an eigenproblem, whose eigenvalues are the corresponding wave-vectors *k _{i}*. The approach, which we call

*k*-formulation, allows to naturally account for frequency dependent permittivities. The numerical properties are comparable to the conventional

*ω*-search for a given wave-vector

*k*

_{0}in terms of convergence rates and numerical cost (CPU time). Using the

*p*-FEM implementation in the numerical C++ library

*Concepts*, we achieved an algebraic convergence rate of up to 7.5 for the

*ω*-formulation and up to 7.2 for the

*k*-formulation, respectively. The calculation was stopped when a relative error of 10

^{-9}was reached. Hence, the method outperforms both, the commercial FE code COMSOL and even more MPB, the standard code for band calculations. For the studied configurations with circular holes the proposed discretization and refinement scheme achieves even an exponential convergence,

*i*.

*e*. the convergence rate increases continuously.

The proposed *k*-formulation is particularly interesting for active PhC devices, where strongly dispersive materials are involved. On the other hand, the computation of slow light modes using the *k*-formulation is possibly not as favorable as with the *ω*-formalism. To detect slow light modes several sampling points in the flat region are needed which are located in a much narrower frequency than wave-vector range, and are therefore more difficult to find with the *k*-formulation. Anyhow, if the rough frequency range is already known, the eigenvalues *k _{i}* and their corresponding modes can be computed accurately at rather low computational costs.

The method can also be interesting for passive devices. The broad band power splitter investigated in [34] is an example, where a dispersive code is favored. The coupling constant of the splitter is given by the difference in the wave-vector of two adjacent guided modes. A tiny shift of one of the bands has a rather large influence on the coupling constant. Within almost the whole band gap range, a modelling error of a few percent is expected only due to the negligence of the dispersive material.

The high accuracy of the method for comparably low computational costs also for dispersive materials turns it into a tool to delve deeply into the study of photonic bands of dispersive periodic structures.

## Acknowledgements

We would like to thank Peter Kaspar for processing the W1 waveguide shown in Fig. 1(a) and Holger Brandsmeier for his very useful Matlab scripts. Further, we would like to thank Prof. Daniel Kressner and Prof. Heinz Jäckel for the fruitful discussions.

## References and links

**1. **J. D. Joannopoulos and S. G. Johnson, *Photonic Crystals*, 2nd ed. (Princeton University Press, 2008).

**2. **M. Agio and C. M. Soukoulis, “Ministop bands in single-defect photonic crystal waveguides,” Phys. Rev. E **64**(5), 055603 (2001). [CrossRef]

**3. **M. Soljacic and J. D. Joannopoulos, “Enhancement of nonlinear effects using photonic crystals,” Nat. Mater. **3**(4), 211–219 (2004). URL http://dx.doi.org/10.1038/nmat1097. [CrossRef] [PubMed]

**4. **T. F. Krauss, “Why do we need slow light,” Nat. Photon. **2**(8), 448–450 (2008). URL http://www.nature.com/nphoton/journal/v2/n8/full/nphoton.2008.139.html. [CrossRef]

**5. **J. Li, T. P. White, L. O’Faolain, A. Gomez-Iglesias, and T. F. Krauss, “Systematic design of flat band slow light in photonic crystal waveguides,” Opt. Express **16**(9), 6227–6232 (2008). URL http://www.opticsexpress.org/abstract.cfm?URI=oe-16-9-6227. [CrossRef] [PubMed]

**6. **S. Kubo, D. Mori, and T. Baba, “Low-group-velocity and low-dispersion slow light in photonic crystal waveguides,” Opt. Lett. **32**(20), 2981–2983 (2007). URL http://ol.osa.org/abstract.cfm?URI=ol-32-20-2981. [CrossRef] [PubMed]

**7. **N. Kono and M. Koshiba, “Three-dimensional finite element analysis of nonreciprocal phase shifts in magneto-photonic crystal waveguides,” Opt. Express **13**(23), 9155–9166 (2005). URL http://www.opticsexpress.org/abstract.cfm?URI=oe-13-23-9155. [CrossRef] [PubMed]

**8. **K. Busch, “Photonic band structure theory: assessment and perspectives,” C. R. Physique **3**(53), 53–66 (2002). [CrossRef]

**9. **D. C. Dobson, J. Gopalakrishnan, and J. E. Pasciak, “An Efficient Method for Band Structure Calculations in 3D Photonic Crystals,” J. Comput. Phys. **161**(2), 668–679 (200). [CrossRef]

**10. **D. Boffi, M. Conforti, and L. Gastaldi, “Modified edge finite elements for photonic crystals,” Numer. Math. **105**(2), 249–266 (2006). [CrossRef]

**11. **K. Schmidt and P. Kauf, “Computation of the band structure of two-dimensional Photonic Crystals with *hp* Finite Elements,” Comp. Meth. App. Mech. Engr. **198**, 1249–1259 (2009). [CrossRef]

**12. **J. Smajic, C. Hafner, K. Rauscher, and D. Erni, “Computation of Radiation Leakage in Photonic Crystal Waveguides,” in *Proc. Progr. Electromagn. Res. Symp. 2004, Pisa, Italy*, pp. 21–24 (2004).

**13. **K. Sakoda, *Optical Properties of Photonic Crystals*, Springer Series in Optical Sciences, 2nd ed. (Springer, Berlin, 2005).

**14. **S. Soussi, “Convergence of the supercell method for defect modes calculations in photonic crystals,” SIAM J. Numer. Anal. **43**, 1175–1201 (2005). [CrossRef]

**15. **P. Kuchment and B. Ong, “On guided waves in photonic crystal waveguides,” in *Waves in Periodic and Random Media, vol. 339 of Contemp. Math*., pp. 105–115 (AMS, Providence, USA, 2004).

**16. **H. Ammari and F. Santosa, “Guided Waves in a Photonic Bandgap Structure with a Line Defect,” SIAM J. Appl. Math. **64**(6), 2018–2033 (2004). URL http://link.aip.org/link/?SMM/64/2018/1. [CrossRef]

**17. **S. G. Johnson and J. D. Joannopoulos, “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Express **8**(3), 173–190 (2001). URL http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173. [CrossRef] [PubMed]

**18. **A. Taflove, *Computational Electrodynamics: The Finite-Difference Time-Domain Method* (Artech House, London, 1995).

**19. **G. Stark, M. Mishrikey, F. Robin, H. Jackel, C. Hafner, R. Vahdieck, and D. Erni, “Position dependence of FDTD mode detection in photonic crystal systems,” Int. J. Numer. Model. **29**, 201–218 (2009). [CrossRef]

**20. **I. Babuška and B. Q. Guo, “The h, p and h-p version of the finite-element method - basic theory and applications,” Adv. Eng. Software **15**, 159–174 (1992). [CrossRef]

**21. **C. Schwab, *p- and hp- Finite Element Methods - Theory and Applications in Solid and Fluid Mechanics* (Oxfold Science Publications, 1998).

**22. **A. von Rhein, S. Greulich-Weber, and R. B. Wehrspohn, “Berechnung von photonischen Kristallen mit Hilfe des Comsol-Elektromagnetikmoduls,” in *Proc. of the COMSOL Users Conf. 2006, Frankfurt, Germany*, pp. 91–96 (2006).

**23. **W Jiang, R. Chen, and X. Lu, “Theory of light refraction at the surface of a photonic crystal,” Phys. Rev. B **71**, 245115 (2005). [CrossRef]

**24. **C. Engstrom and M. Richter, “On the Spectrum of an Operator Pencil with Applications to Wave Propagation in Periodic and Frequency Dependent Materials,” SIAM J. Appl. Math. **70**(1), 231–247 (2009). URL http://link.aip.org/link/1SMM/70/231/1. [CrossRef]

**25. **C. Engstrom, C. Hafner, and K. Schmidt, “Computations of lossy Bloch waves in two-dimensional photonic crystals,” J. Comput. Theor. Nanosci. **6**, 775–783 (2009). [CrossRef]

**26. **P. Frauenfelder and C. Lage, “Concepts - An Object-Oriented Software Package for Partial Differential Equations,” Math. Model. Numer. Anal. **36**(5), 937–951 (2002). URL http://www.edpsciences.org/articles/m2an/abs/2002/05/m2annsl2/m2annsl2.html. [CrossRef]

**27. **K. Busch, G. Schneider, L. Tkeshelashvili, and H. Uecker, “Justification of the nonlinear Schrödinger equation in spatially periodic media,” Z. Angew. Math. Phys. **57**(6), 905–939 (2006). [CrossRef]

**28. **P. Kuchment, *Floquet Theory for Partial Differential Equations* (Birkhäuser Verlag, Basel, 1993). [CrossRef]

**29. **P. Kuchment, “The mathematics of photonic crystals,” in *Mathematical modeling in optical science*, vol. 22 of *Frontiers Appl. Math*., pp. 207–272 (SIAM, Philadelphia, USA, 2001). [CrossRef]

**30. **F. Tisseur and K. Meerbergen, “The Quadratic Eigenvalue Problem,” SIAM Review **43**(2), 235–286 (2001). [CrossRef]

**31. **G. E. Karniadakis and S. J. Sherwin, *Spectral/hp Element Methods for Computational Fluid Dynamics* (Oxford University Press, Oxford, 2005).

**32. ***Webpage of Numerical C*++ *Library Concepts 2* (2009). URL http://www.concepts.math.ethz.ch.

**33. **S. Adachi, “Optical properties of In_{1-x}Ga_{x}As_{y}P_{1-y} alloys,” Phys. Rev. B **39**(17), 12,612–12,621 (1989). [CrossRef]

**34. **P. Strasser, R. Flückiger, R. Wüest, F. Robin, and H. Jackel, “InP-based compact photonic crystal directional coupler with large operation range,” Opt. Express **15**(13), 8472–8478 (2007). URL http://www.opticsexpress.org/abstract.cfm?URI=oe-15-13-8472. [CrossRef] [PubMed]