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 ε(ω_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:

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–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 ω ₀ and not by imposing a certain wave-vector k ₀. The question is hence ’how does the wave propagate when excited with light of frequency ω ₀?’. This real world situation motivates the k-formulations, where we solve for the eigen wave-vectors k_i at a given frequency π ₀. In the k-formulations the corresponding permittivity ε(ω ₀) can be used for each frequency ω ₀ 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 ω ₀. 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:

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 μ ₀ and ε ₀ 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 ₁ and constant permittivity in one direction let it w.l.o.g. be e ₃. 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^-ikx₁ 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 Ẽ and H̃ are periodic in a e ₁. The term ik e ₁ × ũ(k, x) in the shifted curl operator [Eq. (3)] results from applying curl to e^-ikx₁. Here the derivatives in e ₃ 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̃₁(x) ≡ H̃₂(x) ≡ Ẽ₃(x) ≡ 0 and the transverse magnetic (TM) mode with Ẽ₁(x) ≡ Ẽ₂(x) ≡ H̃₃(x) ≡ 0. Each mode is specified by a scalar problem in the out of (e ₁, e ₂)-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 Ẽ(x) := Ẽ₃(x) and H̃(x) := H̃₃(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 ₁-periodicity of the fields Ẽ and H̃ as well as their derivatives. Hence, we can formulate the problem by the following set of equations for the TM mode

 

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 Ω.

Download Full Size | PDF

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 ¹ _per(Ω), the H ¹-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,Ẽ),(ω,k,H̃) ∈ ∝⁺ × B × H ¹ _per(Ω) such that for all test functions v ∊ H ¹ _per(Ω)

By introducing the index ℓ ∊ {-1,0} and powers of the dielectricity ε ^-1(x,ω)= 1/ε(x, ω), ε ⁰(x,ω) = 1 and ε ¹(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 Ẽ = 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 ω-formulation is quadratic in ω only for strictly non-dispersive permittivities for which the problem can be reduced to a linear eigenproblem in λ = ω ².

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 ₁) 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 ϕ_j as test function we obtain the quadratic matrix eigenvalue problem

 

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.

Download Full Size | PDF

Here u⃗ = (u ₁, u ₂,…,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 [Eq. (10)] resulting in the linear matrix eigenvalue problem

The matrix M ₃ (ω) 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 ², M ₃ 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.

 

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.

Download Full Size | PDF

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 ¹ -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_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 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.

 

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.

Download Full Size | PDF

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.

 

Fig. 6. This figure shows the convergence rate of p-FEM and h-FEM based on the k-formulation. 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 $\frac{\omega a}{2\pi c}=0.23$ and $\frac{\omega a}{2\pi 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.

Download Full Size | PDF

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³ ≈ 0.61 if the number of degrees of freedom is increased from N to (N ^1/2 + 3)² ≈ 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 ε(ω ₀) at the investigated frequency ω ₀ 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 ₀, E ₀ + ∆₀, E ₁ + ∆₁ and E ^′ ₀ 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 ω ₀ (blue) for which $\frac{{\omega }_{0}a}{2\pi c}=0.26$ for a = 400nm holds. The corresponding wavelength is λ ₀ = 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.

 

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.

Download Full Size | PDF

To further quantify the error of the wave-vector k, which is introduced by a deviation of the permittivity ε ^′ = ε(ω ₀)-∆ε = 0.26 from that for ω ₀, the permittivity is varied with ∆ε ∊ [0,0.5] while holding ω = ω ₀. 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 ω ₀ 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 ₀ 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.

 

Fig. 8. The error ∆k in the wave-vector k for a given frequency ω ₀ ($\frac{{\omega }_{0}a}{2\pi c}=0.26$ and a = 400nm) as a function of the deviation ∆ε of the permittivity ε(ω ₀) = 10.007. The frequency is kept to ω ₀. 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.

Download Full Size | PDF

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.