1. Introduction

Accurate calculation of electromagnetic field distributions in photonic crystals (PCs) is of fundamental importance to the design of practical micro- and nanophotonics. The motivation for the work reported below is the fact that we have coworkers who have fabricated and characterized two-dimensional (2D) PC slabs [1, 2], that often can be modeled quite accurately by a planar semiconductor slab with a 2D-periodic pattern of perfectly cylindrical air holes. Our goal (not reported here) is the simulation of the transmission and reflection of optical waves by such PC slabs.

A step towards slab simulation is modeling of a single air-PC interface, as shown in Fig. 1, where incoming plane waves from above hit the surface of a 2D PC that fills the half-space below the xy-plane. Above the PC, the electromagnetic field can be expressed as a series involving incident, reflected, diffracted and evanescent plane waves of two polarizations. Inside the PC, the field may be expressed as a series expansion of Bloch-wave solutions for the infinitely extruded PC. In the procedure called mode matching, the amplitudes of the terms are adjusted so that the electric and magnetic fields match on both sides of the boundaries.

Our contribution is the calculation of the needed series of modes (i.e. Bloch waves) of the 2D PC, needed for the mode matching calculation. We have found that to obtain this series of modes, we need knowledge about a part of the band structure of the 2D PC that in general has not been investigated before. We present new results on this part of the band structure below.

 

Fig. 1. Plane waves in three dimensions in air hitting the surface of a 2D PC. We would like to calculate the waves that are reflected and diffracted from the interface between air above the xy-plane and the PC. To do so, modes for a 2D PC extruding infinitely in the z-direction are needed.

Download Full Size | PDF

When the series of modes is known, the scattering matrix method reported in [3] can be used for the mode matching, to compute the reflected and transmitted fields of a PC slab. There is also commercial software [4] that may possibly be adapted to the purpose. The actual mode matching calculation is, however, beyond the scope of the work presented here.

The geometry that we treat is a 2D-periodic pattern of cylindrical rods or air holes, as in Fig. 1 but extruding infinitely in the z-direction. We model this geometry without approximations, and describe analytically the E-field discontinuity at the semiconductor surface around the cylindrical holes. We find that, with the exception of a few propagating low-order modes, all terms in the Bloch-mode series are evanescent. We have calculated band diagrams for real values of the z-component of the wave vectors (low-order modes) as well as imaginary values, corresponding to evanescent modes, modes that decay exponentially away from the boundary (high-order modes). In addition to modes with real and purely imaginary z-components of the wave vector, we also find complex modes (modes with a complex z-component). The existence of complex modes is intimately related to the existence of bandgaps in PCs, since bandgaps are frequency intervals where no propagating Bloch modes are allowed. When the Bloch-mode expansion in the extruded 2D PC is made to match the plane-wave expansion in the air above the PC, the complex modes are needed for completeness of the series.

The need for complex modes in mode matching was recognized already 20 years ago for a class of dielectric microwave guides [5]. Complex modes are in general needed to describe PCs with absorption losses [6–8], but we find that they are also needed to describe lossless PCs. The importance of complex modes has been noted for calculations of reflectivity and transmittivity of semi-infinite non-absorptive PCs, e.g. one-dimensional (1D) gratings [9], 2D PCs that are terminated in a plane perpendicular to the plane of periodicity [10, 11] and three-dimensional PCs [12]. However, the case depicted in Fig. 1, where a semi-infinite 2D PC is terminated in a plane parallel to the plane of periodicity, has not been extensively explored in earlier work.

There are large differences in refractive index between the air, oxide and semiconductor materials that the PC is made from, which represent a challenge for numerical simulations. Current methods for analyzing PCs include plane-wave expansion [10–13],finite-element methods [14], mode matching [15] and rigorous coupled-wave analysis [16, 17]. One can also use techniques that have been developed for the analysis of PC fibers [18–22], for instance scattering-matrix method [23–25]. Today, Bloch-wave calculations for 2D PCs can be performed with a number of software packages, available commercially [17, 26] and as open source [27]. Many of the commonly used methods are reviewed in [13].

Our numerical method is the obvious method to try to use for the particular geometry with circular cylinders. We use Fourier-Bessel series expansion for the electromagnetic fields inside the unit cell, matching the fields analytically at the core-cladding boundary, just like in the standard theory of single-mode optical fibers [14, 28]. The fields are then matched point-wise on opposing sides of the unit cell. The method was implemented in MATLAB, and we show that this series expansion is well behaved and converges very rapidly. The main results that we present do not, however, depend on the choice of numerical method for the simulations.

2. Theory

We consider time-harmonic electromagnetic fields, with the angular frequency ω, and the structure shown in Fig. 1, where a 2D PC is extruding infinitely in the positive z-direction and terminated at z = 0 so that we have air for z < 0. The angular frequency, ω, relates to the wave number for light in vacuum, k ₀, and the speed of light, c, through

With the standard assumption that the PC material is linear, isotropic, lossless and nonmagnetic, the displacement field D and the electric field E are related by a scalar relative permittivity ε(r) [13]. The structure is periodic in the xy-plane, with periods Λ_x and Λ_y in the x- and y-directions respectively.

For this 2D-periodic structure, by using Bloch’s theorem, we can search for solutions for the time-harmonic electric field on the form

where k = (k_x,k_y,k_z) is the Bloch wave vector. e(x,y,z) has the same periodicity as the dielectric structure, i.e.

k_x and k_y are chosen in the first Brillouin zone, i.e. ∣k_x∣ < π/Λ_x and ∣k_y ∣ < π/Λ_y. Note that the z-component k_z of the Bloch wave vector can be nonzero only if the PC is perfectly periodic or uniform in the z-direction. For the rest of this paper the time dependence will be omitted, and the electric field is represented by the vector phasor

that depends on space coordinates but not on time. Similar notation is used for the magnetic field, H(x,y,z).

2.1. Homogeneous region

The doubly periodic function e(x,y,z) in Eq. (4) may be expressed as a double Fourier series expansion with the integers q_x and q_y:

Using this series expansion in Eq. (2), each term will satisfy Maxwell’s equations in the homogeneous region (z < 0) if

so that each term in Eq. (5) yields a plane wave. The frequency and the z-component k_z,q_x,q_y of the wave vector are related to each other through the equation

where ε is the relative permittivity. The general solution can thus be written

For each pair of values q_x and q_y there are two possible polarizations for the plane waves, so an index p have been introduced in the sum above to include for the two polarizations. We observe that each term in the series (8) is a Bloch wave.

For plane waves in air, the electric and magnetic field vectors are orthogonal both to each other and to the wave vector, also if the wave vector has an imaginary z-component. The smaller the frequency ω is, the fewer allowed values of q_x and q_y exist that yield a real k_z,q_x,q_y. Higher-order plane waves, corresponding to larger q_x and q_y, have imaginary k_z,q_x,q_y, corresponding to evanescent plane waves, i.e. waves that decay exponentially with distance from the PC surface in the xy-plane.

2.2. Photonic crystal region

The Bloch waves propagating in the PC region must have the same transverse Bloch vector components k_x and k_y as those in the air region, since the Bloch form in Eq. (4) is valid in every period of the PC, whereas mode matching is performed for a single period only. If we consider incoming plane waves with the wave vector k _inc with components k _x,inc, k _y,inc and k _z,inc, we therefore get the discrete set of plane waves needed for mode matching by choosing k_x = k _x,inc - 2πq_x/Λ_x and k_y = k _y,inc - 2πq_y/Λ_y in Eq. (8). Furthermore, there always exist integers q_x and q_y such that the Bloch vector with components k_x and k_y belongs to the first Brillouin zone in two dimensions.

Since the PC is homogeneous in the z-direction for z > 0, we search for solutions e_q(x,y,z) of the form

The most general solution is thus a series expansion of Bloch waves

where k _q = (k_x,k_y,k_z,q) is a three-dimensional Bloch vector for the 2D PC. Like for the plane waves, only a discrete set of values of k_z,q (eigenvalues) yields solutions to Maxwell’s equations (for a given frequency ω and Bloch vector components k_x and k_y). Each possible value of k_z,q corresponds to a mode e _q(x,y), and the integer q can be called the Bloch mode index.

Central for the analysis of PCs is the calculation of band diagrams, the relations between frequency ω and k_x, k_y and k_z,q. For a given frequency, allowed combinations of k_x, k_y and k_z,q define 2D surfaces in the (k_x,k_y,k_z)-space. In the homogeneous regions the equations defining the surfaces take the simple form (7). Also for the PC regions, in analogy with the plane waves, there exist a limited number of real solutions for k_z,q (low order modes) and an infinite number of imaginary solutions (higher order modes).

The equation for the modes e _q(x,y) is known to be self-adjoint, if it is considered to be an eigenvalue problem for ω ², given that k^x, k_y and k_z are all real [13]. But, seen as an eigenvalue problem for k_z, the equation is not self-adjoint even if the material is lossless. We find in our simulations that in addition to the purely real or imaginary values of k_z,q there can also exist modes for some complex-valued k_z,q. The complex modes must be included when doing the mode matching at z = 0, where the amplitudes of the modes e p,q_x,q_y and eq are adjusted so that the series expansions in Eq. (8) and Eq. (10) match. For the PC structure that we have studied, we have found that if there is a Bloch mode at the Γ point of the Brillouin zone with k_z = k_z,q, then there is also a Bloch mode with k_z = -k_z,q, and if k_z,q is complex, there is also a Bloch mode with k_z = k ^* _z,q.

2.3. Fourier-Bessel expansion in one period of the crystal

Our method for finding the modes in the 2D PC is divided into two steps. First, the time-harmonic Maxwell’s equations [14, 28] are solved analytically inside one period (cell) of the PC as shown in Fig. 2. Then Bloch’s theorem is used to set up point matching along the cell boundary, connecting fields at opposing sides of the cell. The analytical solution inside the unit cell is well-known [14, 23], but the way point matching is used to set up the boundary conditions is original for our work to the best of our knowledge.

 

Fig. 2. One period of the 2D PC. A core with the radius a is surrounded by a cladding with either higher or lower refractive index.

Download Full Size | PDF

We borrow the nomenclature from the fiber-optics, and call the region r < a the core and the surrounding area the cladding, even though this nomenclature is not standard for cylinders in PCs. The relative permittivity is ε ₁ in the core and ε ₂ in the cladding. We define the permittivity contrast as

We introduce cylindrical coordinates as in Fig. 2. Using calculations similar to those in [14], the longitudinal components E_z and H_z can be obtained solving the equation

U denotes either E_z or H_z. Like in [14], we express the transversal field components in terms of E_z and H_z, and solve Eq. (12) using the separation of variables technique, giving a set of solutions

where Ĵ_n and Ŷ_n are Bessel functions of first and second kind respectively, normalized as explained below, and n = 0, ±1, ±2…. Since Ŷ_n(β32r) grows infinitely as r → 0, no such term is allowed for r < a. For r > a, both Ĵ_n and Ŷ_n are needed. For convenience, the definition

has been made, where the branch Re (β_j) ≥ 0 is used.

In our MATLAB code, we use normalized Bessel functions for small function arguments, dividing the Bessel functions of first and second kind respectively with its value for r equal to the distance from the center of the core to the corner of the unit cell or with its value for r equal the cylinder radius. This normalization of the Bessel functions is designed to yield numerically well-behaved expansion coefficients A_n, B_n,C_n, D_n, F_n and G_n in the series Eq. (21) and Eq. (22), by ensuring that

The coefficients A_n, B_n, C_n, D_n, F_n and G_n in Eq. (13)–(14) are related to each other through the boundary conditions at r = a (continuous H_z, E_z, H_φ, E_φ, H_r and εE_r). When the field components satisfy Maxwell’s equations and the boundary conditions are met for four of these six field components, the remaining two boundary conditions will automatically be satisfied. So we have four equations for the boundary conditions, and are left with two unknown coefficients that can be chosen arbitrarily so far. We choose A_n and B_n as independent variables and express the other four coefficients in terms of these two:

Common for all denominators is the factor

The factor $i\sqrt{{\mu }_0/{\epsilon }_0}$ that is included in Eq. (13) ensures that the coefficients A_n, B_n, C_n, D_n, F_n and G_n get same dimension (magnetic field strength), and avoids having the imaginary unit i explicitly in Eq. (16)–(19). We then approximate the z-components of the general E and Hfields by a sum of 2N + 1 contributions from Eq. (13)–(14), where N is the highest order of Bessel functions used:

2.4. Search for Bloch modes in the square periodic PC

The next step is to use Bloch’s theorem to match the fields at points r ₁ and r ₂, at opposite sides of the unit cell, chosen such that they differ with a lattice vector. We assume a square unit cell, i.e. Λ_x = Λ_y = Λ, to allow the points to be equidistant. Since we have 2(2N +1) unknown coefficients A_n and B_n, we need 2(2N +1) equations to determine the unknowns.

 

Fig. 3. Point matching for 12 sampling points around the unit cell of the PC, for the z components of the E and H fields (a) and for the transversal components (b). The small arrows in (b) indicate which transversal component that is matched.

Download Full Size | PDF

The z-components of the electric and magnetic fields are matched at N +1 point pairs as in the example in Fig. 3(a). We get N equations by matching both the electric and magnetic fields at point pairs with

where r ₁ is one of N/2 points at the left edge and r2 is one of N/2 points at the right edge. We get N equations by matching

where r ₁ and r ₂ are one of the points on the lower and the upper edge respectively. x̂ and ŷ are unit vectors in the x- and y-directions respectively. We also get two linearly independent equations by matching the z-components with

In total we get 2(N +1) equations for the z-components.

For the transversal components we choose N point pairs between the former ones, as shown in Fig. 3(b). For the x-components we match the electric and magnetic fields at points separated by Λ in the y-direction, giving us N equations, and for the y-component we use points separated by Λ in the x-direction, giving us additionally N equations. Thus, for the transversal components we get 2N equations, and in total 2(2N +1) equations.

Bloch’s theorem (2) yields for the electric field and for the magnetic field respectively

By using Eq. (16)–(19)in Eq. (13)–(14)and inserting into Eq. (21)–(22) the fields H_z and E_z can be expressed in terms of the unknown constants A_n and B_n. Likewise, transversal components can be expressed in terms of A_n and B_n. Finally the boundary conditions Eq. (26) and Eq. (27), for the point pairs shown in Fig. 3 and in Eq. (23)–(25), are used to get a homogeneous linear set of equations for the variables A_n and B_n:

This equation system can have nonzero solutions for discrete values k_z,q of k_z, as already stated in Sec. 2.2. Since the matrix elements are nonlinear functions of k_z, Eq. (28) represents a nonlinear eigenvalue problem for k_z. This equation system can be written more compactly as

where M(k_z) is the matrix in Eq. (28) and u is the vector containing the 2(2N +1) elements A_n and B_n. This nonlinear eigenvalue problem can be solved by solving

for an arbitrary, nonzero vector v and search for values of k_z for which ∣u′∣ diverges [15]. Actually, Eq. (28) can be considered a nonlinear eigenvalue problem for either one of the four quantities k ₀, k_x, k_y, or k_z, allowing the other three to be specified.

3. Numerical results

The method described in Sec. 2 was implemented in MATLAB, and was used to calculate the modes for the geometry in Fig. 2 with the parameters Λ_x = Λ_y = Λ and a = 0.2Λ. This geometry is the same as one in reference [13], and our method has been verified against the band diagram presented there in Fig. 5.2, where k_z is set to zero and ω is calculated as a function of k_x and k_y. An independent verification of our code is the fact that it yields results that converges to the analytical plane-wave result (7) as the permittivity contrast Δ approaches zero.

For mode matching, we need to find the set of allowed values of k_z for specified ω, k_x and k_y. We confine our attention to k_x = k_y = 0, at the Γ point of the 2D Brillouin zone, believing that the numerical method that finds k_z at the Γ point will find k_z anywhere in the 2D Brillouin zone. The band diagrams that we calculate correspond to Fig. 5.11 in reference [13], but are extended to include modes for imaginary k_z as well.

Our algorithm for calculating the band diagrams starts with finding eigenvalues for k_z = 0, considering Eq. (28) as an eigenvalue problem for k ₀ instead of k_z, finding zeros of 1/∣u′∣. For each eigenvalue, k ₀ is then either iteratively increased or decreased in small steps. In each step the previous value of k ² _z is used as a starting guess in the calculation of the new one (now with Eq. (28) considered as an eigenvalue problem for k_z). If no real solution for k ² _z is found, Muller’s method [29] is used to find complex zeros of 1/∣u′∣.

 

Fig. 4. Automatically computed band diagram for a high-contrast structure. ε ₁ = 8.9, ε ₁ = 1.0. Λ_x = Λ_y = Λ and a = 0.2Λ. k_x = k_y = 0. The highest Bessel function order is N = 10. Two occurrences of complex-conjugated pairs of complex modes are indicated with dotted lines, where the black dotted lines represent Re (k ² _z)/k ² _Λ and the blue dotted lines represent [Re (k ² _z )+Im (k ² _z)]/k ² _Λ. Eigenvalues that are used in Fig. 9 are here marked with small red circles.

Download Full Size | PDF

In Fig. 4 the band diagram for a structure with large permittivity contrast (ε ₁ = 8.9 and ε ₂ = 1.0, giving Δ = 1.6) is shown. The calculations were done with N = 10. k_z and k ₀ are normalized with respect to k _Λ, defined as

One of the important features in the band diagram is that the dispersion curves for k ² ₀ as a function of k ² _z appear to have minima and maxima, leaving intervals (bandgaps) of k ² ₀ where pairs of modes with complex k_z exist, as indicated with dotted lines in Fig. 4. To compute the band diagram in Fig. 4, a preliminary computation with low resolution in k ₀ was first done. The resolution was then manually increased at some intervals of k ₀ that were troublesome, and the whole band diagram was recomputed. This was then repeated a few times. The last computation, which produced the band diagram in Fig. 4, took approximately 10 minutes to compute on a standard personal computer (including a separate computation for complex-valued k_z).

The dispersion relation for a structure with lower permittivity contrast (ε ₁ = 4.0 and ε ₂ = 1.7, giving Δ = 0.81) is shown in Fig. 5 (black lines). We have also plotted the straight lines (green) that Eq. (7) yields for plane waves in a medium with a relative permittivity that equals the average relative permittivity

 

Fig. 5. Automatically computed band diagram (black lines) for a low-contrast structure, ε ₁ = 4.0 and ε ₂ = 1.7. The unit cell geometry is as in Fig. 4. Green lines show plane waves for ε _avg = 2.0 (the average dielectric constant for the low-contrast structure). Figure (b) shows a magnification of the marked area in Fig. (a) where two of the bands form a complex-conjugated pair of k ² _z (marked with dotted lines).

Download Full Size | PDF

The relative permittivities ε ₁ and ε ₂ are chosen such that the average relative permittivity ε _avg is the same in both Fig. 4 and Fig. 5. We observe that many of the bands are nearly straight lines lying close to the plane-wave results.

We explored the dependence of one of the band gap regions in the band diagram on permittivity contrast, expecting the band gap to close with decreasing contrast. We observe the plane-wave approximation (7) to be better, the smaller the permittivity contrast is. Figure 6(a) shows the maximum of the imaginary part of k ² _z in the bandgap. This value clearly decreases as the contrast decreases. We fail, however, to observe a significant reduction in the band gap size (Fig. 6(b)), and the two modes with a complex k_z seem to exist at a significant range of frequencies, also when the two modes are nearly degenerate and k_z is well approximated by the plane-wave result for both modes. To verify that this behavior is not an artifact of using too few orders of Bessel functions N, we have also computed Fig. 6 with N = 18, with no significant changes.

 

Fig. 6. (a) The largest value of Im k ² _z in the bandgap as a function of the permittivity contrast Δ. (b) The upper and lower boundaries (of k ² ₀) of the bandgap and the bandgap size, as a function of the permittivity contrast.

Download Full Size | PDF

To visualize the electromagnetic mode patterns for k ² ₀ = 0.36k ² _Λ and for three of the eigenvalues marked with red circles in Fig. 4, we calculate the quantity

over one period of the PC. Taking the real part of S_z(x,y,z) we get the z-component of the time-average Poynting vector, i.e. the power per unit area crossing the xy-plane.

In Fig. 7 and Fig. 8 the real respectively the imaginary values of S_z(x,y,z = 0) is shown, representing the average power and the reactive power respectively at the plane z = 0. The calculation was done with N = 22 as the highest order of the Bessel functions. For each eigenvalue k_z, S_z(x,y,z) has been normalized so that its largest absolute value inside the unit cell is one.

When k_z is real, S_z(x,y,z = 0) is purely real too, which is seen in Fig. 8(a), where the imaginary part of S_z(x,y,z = 0) is negligible. So there is an average power flow in the z-direction, concentrated in the core of the PC. For purely imaginary k_z (Fig. (b)) the situation is the opposite, with a negligible real part of S_z(x,y,z = 0), meaning that even if the fields are oscillating, on average no power is transported through the PC by these modes. The mode shown is a higher order mode, with larger spatial variation than in Fig. (a). The modes with complex k_z have a complex S_z(x,y,z = 0) (Fig. 7(c) and Fig. 8(c)), implying that these modes in principle carry power, but that the power exhibit exponential decay or growth as a function of z. Such a power flow, decreasing in the z-direction, is unphysical in a lossless PC. Hopefully, the two modes with k_z = k_z,q and k_z = -k ^* _z,q combine to provide a net zero power.

 

Fig. 7. Average power per unit area perpendicular to the xy-plane, Re [S_z(x,y,z = 0)], within one period of the PC, for some allowed k_z.

Download Full Size | PDF

 

Fig. 8. Reactive power, Im [S_z(x,y,z = 0)], within one period of the PC, for some allowed k_z.

Download Full Size | PDF

In the band diagrams in Fig. (4)–(5) the highest order of the Bessel functions used is N = 10, giving a 42 × 42 matrix in Eq. (29). To estimate how the truncation of the sums in Eq. (21)–(22) affects the accuracy of the calculation of k_z, we gradually increase the number of terms used in the calculation of k_z, and see how the value changes. The result for a few selected eigenvalues (marked with red circles in Fig. 4) is shown in Fig. 9, where the values obtained for N = 46 are used as a reference.

 

Fig. 9. An estimate of the error in the calculation of k_z/k _Λ, as a function of the highest order for the Bessel functions, N. The value obtained for N = 46 is used as a reference.

Download Full Size | PDF

The general trend is that the truncation error decreases exponentially with the number of terms used in the Fourier-Bessel series, and that more Bessel function orders are needed to achieve a given accuracy, the larger the imaginary part of k_z is. That can be understood via the observation that modes with a large imaginary k_z (high-order modes) in general have fields that vary more rapidly as a function of position in space than low-order modes. This trend does not prevent some low-order results from being accurate by coincidence, like for k ² _z = -14.6k ² _Λ and N = 22.

The most time-consuming part of the calculations can be divided into two steps; the calculation of the matrix M(k_z), and to solve for u′ in the equation system (30). For large matrices the time needed to compute all the matrix elements is proportional to the number of elements, i.e. to (4N+ 2)², and the time for solving the equation system is known to be proportional to (4N + 2)³ [29]. It is interesting to see which step that is most time-consuming in our calculations. In Fig. 10 computation times for the matrix calculation and the equation system solving (matrix inversion) are compared. Together with the data a quadratic fit, t ~ (4N+ 2)², is shown for the matrix computation and a cubic fit, t ~ (4N + 2)³, for the matrix inversion.

Even if the matrix inversion time asymptotically grows more rapid as N increases, it is clear from Fig. 10 that the calculation of the matrix itself dominates for reasonable values of N (N < 100). Due to the large amount of memory needed, no calculations have been done for very large N to determine where the matrix inversion time start to dominate, but one can still tell that the crossing point occurs for a very large value of N – much larger than anything needed to get a reasonable accuracy, as can be seen in Fig. 9. So, for the range of N values that is needed for our calculations, most of the computing time is spent on evaluation of the elements of the matrix M(k_z) in Eq. (30), and not on solving the equation for the unknown u′.

 

Fig. 10. Computation times for calculating the matrix M(k_z) and taking the inverse, on a standard personal computer.

Download Full Size | PDF

4. Discussion

The aim in this work is to calculate the Bloch modes in a 2D PC of circular cylinders that extends infinitely in the z-direction, to be able to expand the electromagnetic field in a series expansion as in Eq. (10), for a given frequency. For high-order modes, the order of a mode is mainly determined by the imaginary part of k_z, the z-component of the Bloch vector. We need to find all the modes: the low-order propagating modes with a real k_z, the high-order evanescent modes with a purely imaginary k_z, and the many modes with complex k_z that are found in PC structures even if they are lossless with a very small permittivity contrast.

4.1. Comparison with other methods

Cylindrical coordinates and Fourier-Bessel expansions have been used by others, i.e. scattering-matrix method [23–25] and an early method for multicore fibers [18]. In the analysis of photonic crystal fibers (PCFs) [19–22], geometries with cylinders in a periodic pattern appear frequently, but mostly in models using a large but finite number of objects We, on the other hand, are interested in an infinite, perfectly periodic array of cylinders, a 2D PC with a single cylinder per unit cell. We have not found anyone else describing our application of the Fourier-Bessel expansions for Bloch-mode calculations in this PC, with point matching along the edge of the unit cell, and solving a nonlinear matrix eigenvalue problem for the out-of-plane component of the Bloch wave vector.

Our analysis differs from that of PCFs in what 2D modes are needed and how they are used. The evanescent higher-order modes that decay along the direction of the cylinders are usually not relevant for the analysis of PCFs, but those are the modes we are interested in.

Other methods exist that like ours have the advantage of using 1D expansions instead of 2D-expansions, like mode matching in rectangular coordinates [15] and rigorous coupled-wave analysis [16, 17]. Here the structure is approximated by a set of homogeneous rectangular boxes, which has the drawback that, analytically, the electric field perpendicular to the edge diverges at the edge [30], and more so the larger the permittivity contrast is. Therefore a relative large number of orders are needed to obtain accurate results, which have been noted for 1D gratings [31]. To retain an accurate description of the field around the cylinder, independently of the number of terms used in the field expansion, cylindrical coordinates are better suited than rectangular coordinates.

One benefit of our method is that the computed Bloch-wave vector depends continuously on the cylinder radius, also when very few terms are included in the Fourier-Bessel series expansion. A continuous dependence allows the Bloch-wave vector to be easily tracked of as a function of cylinder radius.

Figure 3 indicates that if we do the field matching in 2N points (where N = 6 in Fig. 3) around the unit cell of the 2D PC, we get a spatial resolution in the mode field that roughly corresponds to (N/2)² sampling points inside the unit cell. The analysis of Sec. 2.4 then shows that we need matrices M of dimension 2 (2N + 1) × 2(2N+1) inEq. (29)–(30)to obtain the (N/2)² samples of the Bloch mode field inside the unit cell. A large number of very clever methods have been developed for efficient calculation of the two lowest-order Bloch modes, or a few of the lowest-order modes, of a perfect 2D PC, and many of the methods are reviewed in Appendix D of [13]. As pointed out on p. 257 of [13], if we want to find as many as (N/2)² modes, however, the time taken to perform the matrix operations needed to find all these modes will in general be proportional to N ⁶. Similar reasoning allows us to conclude that the time needed to find (N/2)² modes using Eq. (29) is proportional to N ⁴ if the computation is dominated by matrix element evaluation, and proportional to N ⁵ if the computation is dominated by matrix inversions. So we should get a gentler scaling of computing time with N, by exploiting the knowledge that the unit cell has cylinder geometry.

4.2. Considerations on the method used

There is a wide range of conceivable choices for the point matching around the unit cell boundary. Care must be taken to ensure that the matching procedure always yields a set of equations (28) that are linearly independent, except when a solution for k_z is found. Not all of them will work, so in our work we have tested quite a few different options, and we would like to point out that our paper provides a prescription for setting up point matching equations that fulfills the requirement above. Among the six field components in Eq. (26)–(27), we have chosen to match the components that are parallel to the edges of the unit cell, in points that are evenly spread around the unit cell boundary. The point matching is always done along a curve (a square) where all electric and magnetic field components are analytically well-behaved, resulting in series expansions that converge very rapidly, as shown in Fig. 9.

It should be noted that with the scheme for setting up the equations as in Fig. 3, N must be an even number. Furthermore, for the symmetric case where the core is centered in the unit cell, and k_x = k_y = 0, N cannot be a multiple of 4, for the scheme to yield a set of linearly independent equations. The option to use path integrals instead of point matching to set up the boundary conditions has not been tested, since point matching works well and no problems in converging have been encountered as the number of points is increased. The time needed for the calculations would not scale differently with the number of Bessel orders that is used, if path integrals were used instead of point matching.

With our method an eigenvalue to Eq. (28) appears at k_z = k ₀√ε ₁. An eigenvalue also appears for a value of k_z > k ₀√ε ₁. This eigenvalue approaches k_z = k ₀√ε ₁ as N increases. These eigenvalues are obviously unphysical, and have been discarded.

5. Conclusion

We have studied propagating and evanescent Bloch waves inside one of the canonical two-dimensional photonic crystals (2D PCs); identical circular cylinders in a square pattern, with the interface perpendicular to the cylinders. For a given frequency, all waves are evanescent, except for a few low-order propagating waves. We have presented computed out-of-plane band diagrams for propagating as well as evanescent Bloch waves in the 2D PC. We have found that many of the evanescent Bloch modes do not have a purely imaginary z-component of the Bloch wave vector, but a complex z-component. We believe that the existence of modes with a complex z-component is a property of PCs in general, intimately related to the existence of band gaps in PCs.

To compute the band diagrams, we have used a numerical method based on Fourier-Bessel expansions of the electromagnetic fields, well suited to treat the fields inside and surrounding a cylinder. The method is numerically very well behaved, to the extent that the accuracy of the calculated Bloch vector improves exponentially with the number of terms included in the series expansions. The method involves manipulation of matrices with small dimension, making fast exploration of complex dispersion maps feasible on a desktop computer. To provide assurance that our computer code can be trusted, we have reproduced the band diagrams published in [13] for a structure with high permittivity contrast, and we have observed the computed band diagrams to converge to the analytic plane-wave results when we let the contrast approach zero.

Finally, we believe that our results represent a contribution to the development of a modal expansion method to calculate transmission and reflection coefficients of 2D PC slabs. The existence of modes with a complex z-component of the Bloch wave vector represents a challenge for modal expansions methods in general, regardless of the method used to compute the 2D modes.