1 Introduction

Over the last years, tremendous progress in microfabricating two- and three-dimensional compounds with a spatially periodic dielectric function has lead to a dramatically increased interest in studying these so-called photonic crystals [1]. In such structures, Bragg-scattering of electromagnetic waves leads to the formation of a photonic band-structure, the engineering of which offers enormous potential for applications both in basic science as well as in telecommunication. Due to the advanced state in various planar lithography techniques, the fabrication of two-dimensional photonic crystals has matured to a point which suggests that the first practical applications are imminent [2]. The theoretical description of photonic crystals is generally based on methods of electronic bandstructure theory and is playing a key role in characterizing existing as well as designing novel structures. For instance, the seminal paper of Ho et al. [3] described the first structure that exhibits a complete photonic band gap (PBG), i.e., a range of frequencies for which propagation is forbidden irrespective of the direction of propagation. This work suggested that the formation of a complete PBGis favored by structures exhibiting diamond or hexagonal symmetry and has subsequently lead to the fabrication of the so-called “Yablonovite” and Lincoln-log structures.

In this paper we describe in detail the calculation of photonic band and associated electromagnetic mode structures of two-dimensional photonic crystals. In Section 2, we introduce a novel algorithm for efficient computation of the electromagnetic mode structures. In addition, we outline the calculation of group velocities and effective photon masses in section 3.

2 Two-dimensional photonic crystals

We consider a two-dimensional photonic crystal consisting of a periodic lattice of infinitely extended parallel rods with diameter d. The rods are made from an isotropic dielectric material (dielectric constant ∊_a ) and are embedded in a matrix of dielectric material characterized by an isotropic dielectric constant ∊_b . The direction of the rods defines the z-axis. For electromagnetic waves propagating in the xy-plane, the two transverse polarizations decouple leading to two separate scalar problems. For the TM-polarization (E-field parallel to the rods, i.e., E⃗ (r⃗)≡(0, 0, E(r⃗)), we obtain from Maxwell’s equations

while for the TE-polarization (H-field parallel to the rods, i.e., H⃗ (r⃗)≡(0, 0,H(r⃗)), we have

Here, ∂_x ≡∂/∂x,∂_y ≡∂/∂_y , and r⃗ (x, y) is a two-dimensional vector. The lattice defined through the centers of cylinders is given by the set R={R⃗=n ₁ a⃗₁+n ₂ a⃗2,n ₁, n ₂∊N} of two-dimenisonal lattice vectors R⃗ that are generated by the primitive translations a⃗₁ and a⃗₂. The corresponding reciprocal (dual) lattice is defined through G={G⃗:G⃗·R⃗=2πn, n∊N,R⃗∊R}. The dielectric function ∊(r⃗)=∊_b -(∊_a -∊_b )∑ _R ⃗θ(d/2-|r⃗-R⃗|) is a lattice periodic function which contains all the information of the photonic crystal. Here, θ(r) denotes the Heaviside step function. Eq. (1) and (2) represent differential equations with periodic coefficients. Therefore, their solutions obey the Bloch-Floquet theorem

where i=1, 2. The wave vector k⃗∊1.BZ is a vector of the first Brillouin zone (BZ) known as the crystal momentum.

A straightforward way to solving Eqs. (1) and (2) is to expand all the periodic functions into Fourier series over the reciprocal lattice G, thereby transforming the differential equations into a infinite matrix eigenvalue problem, which is then truncated and solved numerically. Details of this plane wave method (PWM) for isotropic systems can, for instance, be found in [4] and for anisotropic systems in [5].

3 Efficient computation of the mode structure

While the PWM provides a straightforward approach to computing the band structure of photonic crystals, it also exhibits a number of shortcomings. The truncation of Fourier series in the presence of discontinuous changes in the dielectric function may lead to serious convergence problems. In particular, we find that an accurate computation of the electromagnetic mode structure through PWM generally requires more than 2000 plane waves. This leads to very large matrices for which not only eigenvalues but also eigenvectors have to be computed. Finally, the Fourier series for the Bloch functions have to be evaluated. This clearly is a formidable task. In what follows, we describe the details of a new approach to efficiently compute the electromagnetic mode structure of two-dimensional photonic crystals.

The pairs of Eqs. (1) and (3) as well as Eqs. (2) and (4) comprise elliptic partial differential eigenproblems with Bloch-type boundary conditions at the boundaries of the unit cell. Upon a real-space discretization of the differential operators, both problems for any wave vector k⃗ may be cast into the matrix form of a linear system of equations

where Λ _i =${\omega }_i^2$ /c ² and U_i represents the discretized fields E _k⃗(r⃗) and H _k⃗(r⃗) in case of TM- and TE-polarization, respectively. The precise form of the operator matrix L _k⃗ and the vectors U_i depend on the type of polarization (TE or TM) and the type of discretization (finite difference or finite element method). Mathematically, the problem is then to find approximations to the first few eigenfrequencies (bands) Λ₁≤Λ₂≤Λ₃≤… and associated Bloch functions U ₁; U ₂; U ₃…. Due to the local nature of the underlying differential operators, the operator matrix L _k⃗ is sparse, which suggests that Eq. (5) should be solved iteratively by Rayleigh quotient iteration or Lanczos. Such approaches of treating the discrete problem as a purely algebraic one can result in loss of valuable information, especially concerning the smoothness of the Bloch functions. In general, the Bloch functions corresponding to the desired first bands are very smooth, so that they are well approximated on coarser grids. Certain multigrid methods take full advantage of this smoothness and are therefore very effective for solving such problems [6]. In addition, multigrid methods can be applied to nonlinear problems such as Eq. (5), where both eigenfrequency and Bloch functions are unknown. To do so, Eq. (5) has to be amended by a ortho-normalization condition [7]

where V denotes the volume of the Wigner-Seitz cell (WSC) and δ_nm is the Kronecker symbol for the band indices n and m.

The efficiency of multigrid methods results from the fact that, although iteration (relax-ation) of Eq. (5) using a simplified version of the matrix operator L _k⃗ (Jacobi- or Gauss-Seidel iteration [6]) is usually slow to converge, it is quick to reduce high-frequency error components [6]. This allows the problem to be transfered to a coarser grid where the error can be resolved with much less work. Not only is the relaxation cheaper per sweep on coarser grids, but the solution process is also much more effective. The coarse grid equation can be solved by relaxation and appeal to still coarser grids. The coarsest grid used is chosen so that solution of the problem there is inexpensive compared to the work performed on the fine grid. The number of relaxation sweeps needed to smooth the error on each grid is generally small.

We discretize Eqs. (1) and (2) using a finite volume discretization scheme on a uniform mesh and employ a simple white-black ordered Gauss-Seidel relaxation method [6]. The multigrid (MG) iteration follows a V-cycle as illustrated in Fig. 1. The relaxation steps on all but the coarsest level are performed with fixed eigenfrequencies Λ _i i=1,2,…, q. The eigenfrequencies are updated after every relaxation step on the base level using the Rayleigh quotient

 

Fig. 1. Graphical illustration of the V-cycle used in the multigrid iteration of Eq. (5) depicting four levels of grids. A movie (881.mov(0.8M)) shows the multigrid iteration using two levels for the third band at the x-point in the TM-polarized case (see Fig. 2).

Download Full Size | PDF

where the superscript 1 indicates the representation of L _k⃗ and U_i on the lowest level and 〈.,.〉 stands for the discrete version of scalar products defined in Eq. (6). This V-cycle scheme is repeated until stable values for eigenfrequencies Λ _i and Bloch-functions U_i emerge (generally two V-cycles are sufficient) [7].

As an example, we consider a square lattice of dielectric rods (∊_a =13) in air (∊_b =1). Then, the primitive translations are given by a⃗₁=(a, 0) and a⃗₂=(0, a), where a is the lattice constant. The radius r of the rods is r/a=0.45. The photonic band structure of this structure is shown in Fig. 2. In Tab. 1, we compare the band structure data

 

Fig. 2. Photonic band structure for a square lattice of dielectric rods (∊_a =13, r/a=0.45) in air (∊_b =1) for TM-polarization (blue) and TE-polarization (green).

Download Full Size | PDF

from the MG-method for the first four bands of TM-polarization using different mesh sizes with corresponding data using PWM with different numbers N of plane waves. We observe, that for large numbers of plane waves, the PWM aproaches the values of the MG-method to within 1% but requires about 100% more CPU-time. In addition, the result from the MG-computations suggest that a 128×128 mesh is sufficient to obtain converged results. It should be noted that in this example PWM was used to compute the band structure only. If an evaluation of the electromagnetic mode structure were

Table 1. Comparison of numerical data obtained from PWM (upper part) and the multigrid method (lower part) using different numbers N of plane waves and different mesh-sizes (128×128 and 256×256), respectively. The values for the defect give an estimate of the absolute error of the Bloch-functions [7]. The computations have been performed at the X-point for TM-polarization using a 466 MHz Celeron processor.

View Table

required, the computational effort would increase substantially leading to much longer CPU-times as quoted above. In contrast, the MG-method obtains the Bloch functions without any additional numerical effort.

4 Computation of group velocities and effective photon masses

In order to understand pulse propagation in photonic crystals, additional calculations have to be performed. In this section we demonstrate how to obtain group velocities and effective masses characterizing the speed of pulse propagation and the broadening of pulses. This approach closely follows the so-called kp-perturbation theory (kp-PT) of electronic band structure theory and has been applied to one-dimensional systems [8]. For simplicity, we will only discuss the TM-polarization.

We begin by writing the Bloch-Floquet theorem, Eq. (3), as

where u _k⃗(r⃗) is a lattice periodic function. The equation of motion for E _k⃗(r⃗), Eq. (1) may now be transformed into an equation of motion for u _k⃗(r⃗)

Introducing the operator Ĥ (k⃗)=Δ+2i∇·k⃗+k⃗², where Δ=${\mathit{\partial }}_x^2$ +${\mathit{\partial }}_y^2$, allows to cast the equation of motion, Eq. (10), for u _k⃗+q⃗^(r⃗) (|q⃗≪|π/a) into the form

Here, we have introduced Ω⃗=-i(∇+ik⃗). Eq. (11) together with the fact that |q⃗|≪π/a suggests that we treat the second and third part on its l.h.s as a perturbation on Ĥ(k⃗) such that the eigenfrequency ω _k⃗+q⃗ may be regarded as a perturbed eigenvalue and compared to a Taylor-expansion of ω _k⃗+q⃗ around k⃗

 

Fig. 3. Group velocity of the first three bands (first band: black; second band: red; third band: green) for a square lattice of dielectric rods (∊_a =13, r/a=0.45) in air (∊_b =1) for TM-polarized light. The corresponding band structure is displayed in Fig. 2.

Download Full Size | PDF

The computation of group velocities υ⃗_k⃗=∂ _k ⃗ω _k⃗ and the inverse effective photon mass tensors (group velocity dispersion tensors) $M_{\vec{k}}^{\mathit{\text{ij}}}$=∂_ki ∂ _kj ω _k⃗, i=1, 2 are thus reduced to second order perturbation theory [8].

In Fig. 3 we show the group velocity for the first three TM-polarized bands along the high-symmetry lines of the 1.BZ. As expected from the corresponding band structure in Fig. 2, the higher bands exhibit rather low values of the group velocity. In particular, the values of group velocity for band 3 along the Γ-X direction are about an order of magnitude smaller than the group velocity in band 1, making band 3 an excellent candidate for realizing nonlinear optical properties in this material.

5 Discussion

In summary, we have introduced a novel scheme for computing photonic band and associated electromagnetic mode structures based on multigrid methods. This approach is much more efficient than the plane wave method especially when the electromagnetic mode structure is required as, for instance, is the case when calculating group velocities and effective photon masses in photonic crystal. The MG-approach may be further refined using finite elements instead of finite differences as well as adaptive instead of uniform grids. An extension of this approach to three-dimensional photonic band structures is straightforward. In addition, we have shown how group velocities and effective photon masses may be calulated directly from photonic band structure computations through the photonic analogue of the well-known kp-perturbation theory. These results are of direct relevance to studies of pulse propagation in photonic crystals.