Efficient Computation of Photonic Crystal Waveguide Modes with Dispersive Material

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 OCIS codes: (000.4430) Numerical approximation and analysis; (130.5296) Photonic crystal waveguides 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), 055,603 (2001). 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. 4. T. F. Krauss, “Why do we need slow light,” Nat. Photonics 2(8), 448–450 (2008). URL http://www.nature.com/nphoton/journal/v2/n8/full/nphoton.2008.139.html. 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. 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. 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. 8. K. Busch, “Photonic band structure theory: assessment and perspectives,” C. R. Physique 3(53), 53–66 (2002). 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). 10. D. Boffi, M. Conforti, and L. Gastaldi, “Modified edge finite elements for photonic crystals,” Numer. Math. 105(2), 249–266 (2006). 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). 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). 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. 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. 18. A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, London, 1995). 19. G. Stark, M. Mishrikey, F. Robin, H. Jäckel, 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). 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. In Engineering Software 15, 159–174 (1992). 21. C. Schwab, pand hpFinite 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, 245,115 (2005). 24. C. Engström 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/?SMM/70/231/1. 25. C. Engström, C. Hafner, and K. Schmidt, “Computations of lossy Bloch waves in two-dimensional photonic crystals,” J. Comput. Theor. Nanosci. 6, 775–783 (2009). 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/m2anns12/m2anns12.html. 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). 28. P. Kuchment, Floquet Theory for Partial Differential Equations (Birkhäuser Verlag, Basel, 1993). 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). 30. F. Tisseur and K. Meerbergen, “The Quadratic Eigenvalue Problem,” SIAM Review 43(2), 235–286 (2001). 31. G. E. Karniadakis and S. J. Sherwin, Spectral/hp Element Methods for Computational Fluid Dynamics (Oxford University Press, Oxford, 2005). 32. Concepts Development Team, Webpage of Numerical C++ Library Concepts 2 (2009). URL http://www.concepts.math.ethz.ch. 33. S. Adachi, “Optical properties of In1−xGaxAsyP1−y alloys,” Phys. Rev. B 39(17), 12,612–12,621 (1989). 34. P. Strasser, R. Flückiger, R. Wüest, F. Robin, and H. Jäckel, “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.


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 ε(ω i ) 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 k for a given ω i 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 k [7], we are able to directly access the lossy modes by computing the complex eigen wave-vectors k i .However, computing losses is not the object of this article.In the presentation of our method we will emphasize three different aspects: 1. the formulation of the mathematical eigenvalue problem for k i of a given ω i , 2. the modelling of guided PhC modes in an infinite geometry, and 3. 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 i ) for each wave-vector sample k i is computed for the unit cell (see [1,8,9,10,11] and the references therein).We call the formulation underlying this procedure the ω-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 absorp- tion 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),(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 semiinfinite 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: 1. 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.2. 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.

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),(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 T [27] formally defined by which is with the exception of the factor e −ikx 1 equivalent to the Floquet transform F [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 (2) to the time-harmonic Maxwell equations ( 1) and the definition of the shifted curl operator which is characterized by the property The transformed fields E and H are periodic in ae 1 .The term ike 1 × u(k, x) in the shifted curl operator (3) results from applying curl to e −ikx 1 .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 and The Bloch transform (2) guarantees the ae 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 ( 5) and ( 6) are stated on the two-dimensional strip Ω (see Fig. 2(a) and (b)) and completed by following the super-cell approach by using periodicity conditions on the upper and lower boundaries.We write the general equation (6c), which results from the Bloch transform for the case of a super-cell with discontinuities of the permittivity at its borders.

The general eigenvalue problem
Now, we will state the general eigenvalue problems ( 5) and ( 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: 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: 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

Computational frameworks for the calculation of the modes
The obtained equations ( 9) for ℓ = −1, 0 can be solved by two different methods: • The search for the eigenfrequencies ω for particular values of k i which we call the ω-formulation.
• The search for the eigen wave-vector k for particular frequencies ω i which we call the 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 (5) or ( 6): the sampling scheme over frequency values and the quadratic -and so simply to linearizeeigenvalue 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 quasiperiodicity 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.

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 FE = ∑ N j=1 u j φ j .In- serting this expression in the eigenproblem (9) with a frequency dependent permittivity ε(ω, x) and each basis function φ j as test function we obtain the 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  With x := (k u, u) we can simply linearize [30] equation (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 (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 (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)-(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.

Discretization by the p-version of the FEM
In the previous section we have defined the matrix eigenvalue problem (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 1/2 ) = 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 [21,11] 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 kformulation.The used mesh is the one of Fig. 3(d) and some modes are depicted in Fig. 3(a)-(c), which are obtained with the k-formulation using a polynomial degree of p = 8.

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.

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 supercell 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 i ) 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 (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 prefinement 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  5: A comparison of the convergence for various codes, each using an ω-formulation.We plot the average error of all eigenfrequency values, which are part of the dispersion of the guided PhC modes versus the degree of freedom (left) and computing time (right).MPB (red) shows a algebraic convergence, COMSOL (blue) and Concepts with h-refinement (green) converge quadratically.The best convergence is achieved with Concepts and refinement of the polynomial degree of the basis functions (black curve).The configuration is that of the test PhC waveguide of Sec.3.3.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.

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 2πc a , and similar to the previous The test configuration is that of Sec.3.3.The shown error is the averaged error of the wave-vector eigenvalues of the guided modes for frequency values between ωa 2πc = 0.23 and ωa 2πc = 0.28 which is shown in dependence on the degrees of freedom (left) and the computational time (right).The curve for h-refinement and uniform polynomial degree p = 2 (dashed line) reveals coarsely the behavior of a potential COMSOL implementation of the k-formulation.The p-FEM curves were fitted both to an algebraic convergence (convergence rate of 7.2) as well as to an exponential law.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 wavevectors 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 + 6N −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 i , becomes smaller for increasing N i and tends to zero.Note, that the convergence curve fits locally also to an algebraic law.
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.

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 920 nm 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 ω 0 a 2π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 [ ωa 2πc ] and a small 'up' shift of the lower band gap frequency from 0.2252 to 0.2264 [ ωa 2πc ] -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.
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 Fig. 7: On the left side, a comparison of the band diagrams for the case of a constant dielectric constant (blue) and for the case of a dispersive permittivity (red) are shown.Both have been computed using the k-formulation.Despite their non-scalability we labelled the frequency and the wave-vector for comparison purpose in dimensionless units.The absolute wave-length is labelled on the right axis.On the right side the real part of the permittivity model of Adachi [33] for InP is plotted.
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.

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 Fig. 8: The error ∆k in the wave-vector k for a given frequency ω 0 ( ω 0 a 2πc = 0.26 and a = 400nm) as a function of the deviation ∆ε of the permittivity ε(ω 0 ) = 10.007.The frequency is kept to ω 0 .For the k-formulation for each frequency the right permittivity value can be used directly (corresponding to ∆ε = 0) whereas with the ω-formulation an iteration procedure has to be used.
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.

Fig. 1 :
Fig. 1: (a) Scanning electron microscope (SEM) picture showing the top view of a PhC W1 waveguide fabricated by etching deep holes in a InP/InGaAsP/InP layer structure with finite PhC-width.(b) 3D sketch of a supercell including the vertical layer structure.

Fig. 2 :
Fig. 2: (a) Illustration of a PhC W1 waveguide in 2D based on a hexagonal lattice with three rows of holes on both sides of the channel and a possible super-cell Ω.(b) Analogous illustration of PhC W1 waveguide based on a square lattice with a super-cell Ω.

Fig. 3 :
Fig. 3: In the sub-figures (a)-(c) the real part of the magnetic field component for three typical guided PhC waveguides modes at different frequencies are shown: (a) an even mode where the slope of the dispersion is relatively steep, (b) another even mode where the slope is flat and (c) a typical odd mode.In sub-figure (d) the mesh of the super-cell is illustrated.The configuration is that of Sec.3.3.
(a) Band diagram with equally spaced wave-vectors.(b) Band diagram with equally spaced frequencies.

Fig. 4 :
Fig. 4: The same band diagram, once computed with the ω-formulation and uniformly distributed wave-vectors k i shown in (a), and once -see (b) -with the k-formulation and uniformly distributed frequencies ω i .The PhC configuration is that introduced in Sec. 3.3.
Fig.5: A comparison of the convergence for various codes, each using an ω-formulation.We plot the average error of all eigenfrequency values, which are part of the dispersion of the guided PhC modes versus the degree of freedom (left) and computing time (right).MPB (red) shows a algebraic convergence, COMSOL (blue) and Concepts with h-refinement (green) converge quadratically.The best convergence is achieved with Concepts and refinement of the polynomial degree of the basis functions (black curve).The configuration is that of the test PhC waveguide of Sec.3.3.

Fig. 6 :
Fig.6: This figure shows the convergence rate of p-FEM and h-FEM based on the kformulation.The test configuration is that of Sec.3.3.The shown error is the averaged error of the wave-vector eigenvalues of the guided modes for frequency values between ωa 2πc = 0.23 and ωa 2πc = 0.28 which is shown in dependence on the degrees of freedom (left) and the computational time (right).The curve for h-refinement and uniform polynomial degree p = 2 (dashed line) reveals coarsely the behavior of a potential COMSOL implementation of the k-formulation.The p-FEM curves were fitted both to an algebraic convergence (convergence rate of 7.2) as well as to an exponential law.