1. Introduction

In the early years of the study of photonic band gap systems, the principal aim usually was to create structures having total band gaps, i.e., ranges of frequency in which electromagnetic waves were non-propagating irrespective of their direction. With this achieved, researchers went on to study the creation of defects in photonic crystals [1] that form the foundation of useful devices such as cavities, waveguides, etc. In recent years, a number of workers have turned their attention to unusual properties of photonic crystals in frequency ranges for which electromagnetic waves can propagate, e.g., high dispersion and slow light. The latter case has drawn attention [2, 3] to the band-edge properties of photonic crystals, where manipulation of the photonic dispersion relation for frequencies very close to a band gap edge can lead to low group velocities over extended frequency ranges.

Here, we direct our attention to the gap-edge properties of photonic crystals, and in particular to the properties of defect modes in two-dimensional (2D) photonic crystals, created by making a small perturbation of small, finite extent in an otherwise infinite, periodic photonic crystal (PC). Such defect modes are well understood for electrons [4], and since the wave equation and the boundary conditions are similar for the 2D electromagnetic case in which the electric field is polarized parallel to the cylinder axes (i.e., E _∥ or TM polarization), one might be confident that the behaviour for electrons should carry over to photons. However, for an electric field polarized perpendicular to the axes (i.e, H _∥ or TE polarization), the boundary conditions for electrons and photons differ, and so it is an open issue whether the asymptotics for TE polarization are similar to, or quite different from, that for electrons.

The approach we adopt here combines both analytic and numerical studies which yield results that are analogous to those obtained by Economou [4] in quantum mechanics. In particular, the scalar results obtained in Ref. [4] by the tight binding method carry over to the vector case for photons. Our method requires us to express the Green function as the superposition of quasiperiodic field problems [5]. We show that in the vicinity of a band edge, the Green function in this 2D system is dominated by the contribution of the edge mode, which has a logarithmic dependence on the frequency difference between the given frequency and the band edge frequency. The domination of the Green function by the contribution of a single mode is a key requirement of the method, allowing us to undertake a first-order perturbation treatment, and leading to the exponential dependence of the difference between the defect mode frequency and the frequency of the band edge on the reciprocal of the defect perturbation parameter.

The numerical studies discussed in Secs 3 and 4 exploit the Fictitious Source Superposition (FSS) method [6,7] that we developed especially for the purpose of studying the characteristics of highly extended defect modes that occur near band edges. In Sec. 4, we consider a range of examples for both TM and TE polarizations, for square and hexagonal lattices, and for the first and second band gaps, and show that the analytically derived exponential dependence holds quite generally for 2D photonic crystals and accurately fits the evolution of the defect mode with increasing perturbation.

2. Analytic method

We consider a two-dimensional photonic crystal (PC), specified by periodically varying functions n(r) and ε(r) = n ²(r), respectively giving the refractive index and dielectric constant (permittivity) in the array. In our treatment, we assume that ε(r) is a real function of position so that the system is lossless. The Green function for the PC is expanded as a superposition of quasiperiodic Green functions. In turn, each of these is expanded in the basis of Bloch functions {ψ_m(k ₀,r)} which obey the Helmholtz equation, the Bloch condition associated with the Bloch vector k ₀, and polarization dependent boundary conditions on the interfaces between cylindrical inclusions and the matrix material in which they are placed.

We introduce functions characterizing the field (ψ _∥) in the PC: for TM (E _∥) polarization, we define ψ = E_z, p(r) = 1, s(r) = ε(r), and, for TE (H _∥) polarization, we introduce ψ = H_z, p(r) = 1/ε(r),s(r) = 1. In terms of these, the Helmholtz equation, for either polarization, can be written as

 

Fig. 1. Band surface for TE polarization for a hole-type hexagonal lattice with holes of normalized radius a/d = 0.2 and refractive index n_c = 1, in a background of refractive index n_b = 3. The frequency range corresponds to the lower edge of the first band gap shown in Fig. 3 (b), and shows the maxima (in red) at the six equivalent K points that characterize the band edge.

Download Full Size | PDF

where ℒ is the operator

which is Hermitian with respect to the inner product

where the integration runs over the Wigner-Seitz (WSC) cell. The Bloch functions {ψ_m} are orthogonal with respect to this inner product and satisfy

where M_n = 〈ψ_n, ψ_n〉.

The Green function G that we require represents the electromagnetic field corresponding to a single source at a point r′, and satisfies the differential equation

that holds everywhere, also encapsulating the field continuity conditions. It is expressed as a superposition [5] of the quasiperiodic Green functions (expanded in the Bloch mode basis), with the superposition requiring an integration over the 2D Brillouin zone (BZ)

We next consider the situation where the frequency ω lies in a complete band gap, and close to a gap edge which has the frequency ω_L, using L to denote the closest band to the frequency ω. We begin with the simplest (nondegenerate) case for which the band edge occurs when the Bloch vector k ₀ = k _L. To lowest order, we will assume that the L ^th band surface is parabolic near its edge, i.e.,

where C_L characterizes the band curvature and is equivalent to the effective mass of the electron in semiconductor theory. Using this parabolic form, we can approximate the required integral in Eq. (6) to leading order. When the band edge corresponds to a single, nondegenerate Bloch vector k ₀ that is contained completely within the Brillouin zone (BZ), we have

The presence of the logarithm term [5, 8] in Eq. (8) ensures that the contribution from band edge L to the Green function (6) dominates all other contributions, provided that ω is very close to ω_L. We thus deduce

In cases of high symmetry, the band edge can correspond to several values of k ₀, each of which can contribute to the Green function. For example, for a hexagonal lattice, the lower edge of the first band gap (see Fig. 1) is characterized by maxima in the frequency surface at the six equivalent K points of the BZ, with each necessitating an integration over an internal angle of ϕ_L,j = 2π/3. Correspondingly, the upper edge of this band gap is characterized by minima at the six equivalent M points, each with internal angle ϕ_L,j = π. Accordingly, the expression on the right hand side of Eq. (9), corresponding to the single Bloch vector k_L, must thus be replaced by a sum over those k_L,j, values that correspond to the band edge, with each contribution weighted by θ_L,j = ϕ_L,j/(2π).

The final step in the derivation of the dispersion relation for the defect mode (when it is close to the band edge) is to apply Green’s Theorem to an infinite photonic crystal in which the dielectric constant distribution is perturbed in some region C ₀ to create the defect mode. The defect mode, which is localized and tends to zero exponentially with increasing distance from C ₀, satisfies the Helmholtz equation ▽ ∙ (p̃(r)▽ψ(r)) + (ω ²/c ²)s̃(r) ψ(r) = 0, in which the perturbation is characterized by the quantities p̃(r) and s̃(r).

Then, using Green’s Theorem, we may derive the exact result

Up to this point, the treatment is general, applying to both TE and TM polarizations. We now consider the special case of TE (H _∥) polarization, for which p(r′) = 1/ε(r′), p̃(r′) = 1/ε(r′), and s(r′) = s̃(r′) = 1, and substituting into Eq. (10) the leading order estimate (9) for G, we obtain

where j sums over all values of k ₀ = k_L,j that correspond to the edge of band L. The evaluation of this expression requires some delicacy, given the dependence on both the geometry of the Wigner-Seitz cell and the spatial characteristics of the mode. However, expression (11) simplifies considerably in cases of high symmetry, since the bulk modes at the band edge can be obtained from one another by an appropriate rotation transformation.

We now follow a particular example in order to outline clearly the approach, and which also allows us to infer some general conclusions. We consider the lower edge of the band gap, shown in Fig. 1, for a hexagonal lattice, which is associated with the six equivalent K points of the BZ. The ψ_L,j are almost circular symmetric near the cylinder and so may be taken to be effectively identical in the vicinity of C ₀. Furthermore, for this case, θ_L,j = 1/3, C_L,j = C_L, and M_L,j = M_L, for all j.

We now work in the framework of first order perturbation theory, in which the eigenfunctions are not perturbed while the eigenvalues may vary. Accordingly, we make the approximation ψ(r) = ψ(k _L,1,r), corresponding to a particular mode. The sum in Eq. (11) now contains 6 identical terms each of which are weighted by θ_L,j = 1/3. This leads us to define a geometry factor f_L, for which here f_L, = 6 × 1/3 = 2. Of course, for a singularity that is completely enclosed within the BZ, f_L = 1.

We can further simplify the combination of the integral of Eq. (11) and the normalization term M_L by noting that the transverse resolute of the electric field is E_t = -i(z̃ × ▽ψ)/(ωε). Thus, for all j, this leads to

after recasting the normalization factor as

which follows from the differential equation and an application of Green’s Theorem. Apart from the factor c ²/ω ² _L, the result in Eq. (12) is the ratio of the change in the electric energy (∫E∙D d ² r) resulting from the defect, to the electric energy of the unperturbed mode within the Wigner-Seitz cell. Since the perturbed dielectric constant $\tilde{\epsilon }$ (r) differs from ε(r) only in the region C ₀, it follows that the integral in the numerator can be calculated over the entire Wigner-Seitz cell. Thus,

in which 𝓔 denotes the electric energy contained within the region specified by the subscript and

where N_L = |C_L|f_LA_WSC (2π) is the Density of States at the band edge L[4, 5].

Finally, we find the most general form of the result for the frequency shift generated by the introduction of a defect in terms of the relative change in the electric energy as

where 𝓐 is a constant which cannot be determined from first order perturbation theory. The results in the form presented in Eq. (16) are entirely general, applying also for TM polarization. We briefly outline the derivation in this case which is somewhat simpler because it involves directly the electric field component ψ = E_z. We now have p(r′) = p̃(r̃) = 1, s(r′) = ε(r′) and s̃(r′) = $\tilde{\epsilon }$ (r′) and, as for TE polarization, the substitution of these relations into Eq. (10) leads to

for the TM polarization; the relations (14)–(16) then follow from the fact that E = (0,0, ψ).

From the general result (16), we may derive particular forms. In Section 4 we consider PCs in which the each unit cell has a constant refractive index n_b, except for a circular inclusion which has the refractive index n_c. For a defect caused by changing the refractive index of a single cylinder to n_d, we have δ𝓔_C0/𝓔_C0 = δε_C0/𝓔_C0 and so we derive the result

which is applicable to either polarization, and in which

denotes the sensitivity parameter.

3. Numerical details

To demonstrate the method and verify the theoretical analysis, we compute the evolution of the defect mode using the FSS method, and then proceed to fit the exponential form (18), thus yielding a numerical estimate of the sensitivity parameter. We then compare these with values obtained from the expressions derived in the previous section. The evaluation of the parameters occurring in the expressions for the sensitivity factor requires an accurate knowledge of the relevant Bloch function at the band edge, together with its frequency and band curvature. However, this is a routine calculation and such quantities may be determined to high accuracy using either the multipole method, as described in Ref. [5], or other techniques, including the finite element method [9] which has been used for the computations here.

The testing of the exponential dependence of the defect frequency difference on δε_c, when both are small, is computationally difficult, since the defect mode becomes arbitrarily extended, the more closely we search for modes near to the band edge. While ultimately any numerical technique fails in such a search, the FSS method [6, 7] is better positioned than any alternative technique to locate accurately modes that are very close to the band edge. The FSS method exploits a fictitious source to tailor the existence of a defect, and then constructs the mode as a superposition of solutions of quasiperiodic field problems. The superposition entails an integration over the 2D Brillouin zone and thus the computational complexity of the method is comparable to that of 2D supercell methods which become overwhelmed by the demands of handling a supercell that is sufficiently large to enclose the defect mode. However, in our implementation of the FSS, we model the defect structure as a diffraction grating surrounded by two semi-infinite PCs, thus reducing the superposition task to require an integration over only a 1D Brillouin zone [6, 7], and permitting us to locate modes that are extremely close to the band edge.

 

Fig. 2. For a square rod-type PC with n_b = 1, n_c = 3 and a/d = 0.3 operated in TM polarization: (a) Band diagram showing the first two band gaps. (b) Evolution of the defect mode in each of the gaps with changing defect refractive index n_d. Thick, red curves show the defect mode frequency computed by the FSS method. The dashed black curves indicate analytic result (18), with prefactor 𝓐, obtained as described in the text. The vertical dashed lines indicate the refractive index n_c of unperturbed cylinders, while the horizontal dashed lines show the band gap edges. The dotted blue curves are for other defect modes.

Download Full Size | PDF

 

Fig. 3. Band diagrams for the hexagonal lattices for (a) TM polarization in a rod-type PC with n_b = 1, n_c = 3, and a/d = 0.2, and (b) TE polarization in a hole-type PC with n_b = 3, n_c = 1, and a/d = 0.2.

Download Full Size | PDF

4. Defect modes

In this section, we demonstrate and validate the theory by studying the evolution of defect modes constructed by varying the refractive index of a single cylinder in both rod-type and hole-type PCs, for both square and hexagonal lattices. We compute the sensitivity parameter 𝓢_asymp from the asymptotic analysis (19) and then compare this with the estimate 𝓢_FSS obtained by fitting the evolution curve obtained from the FSS calculations to the model ln |$\tilde{\omega }$ - $\tilde{\omega }$ _L| = ln 𝓐_FSS - ε_C0𝓢_FSS/δε in the vicinity of the band edges (10^-6 <|$\tilde{\omega }$ - $\tilde{\omega }$ _L| < 10^-3). From this least squares fit, we determine values for the parameters 𝓐_FSS and 𝓢_FSS.

We begin with a square symmetric rod-type PC operated in TM polarized light, with the refractive index and radius of the unperturbed cylinders being n_c = 3.0 and a/d =0.3, and with a background refractive index of n_b = 1. Fig. 2 (a) shows the band diagram which has band gaps for ω̃ ≡ ωd/(2πc) = d/λ ∈[0.265265,0.334947] and ω̃ ∈ [0.468027,0.562353]. In Fig. 2 (b) we show the evolution of a defect mode, in each of the first and second band gaps, generated by varying the refractive index (n_d) of the single defect cylinder. For n_d > n_c in the second band gap, there are two additional defect modes indicated in blue. These originate at a higher band edge and are not considered in our analysis.

The dashed black curves in Fig. 2 (b) indicate the analytic result (18), with 𝓐 = 𝓐_FSS and 𝓢 = 𝓢_FSS while the red curves display the results of the FSS calculations. Note the excellent agreement between the numerical and analytic results for frequencies close to the band edge, where the analytic result is valid. Note also the excellent agreement between the analytic (𝓢_asymp) and fitted (𝓢_FSS) values of the sensitivity parameter, shown in Table 1, which confirm the validity of the asymptotic analytic method.

Table 1. Overview of the results for the structure in Fig. 2. For the four lowest band edges we show the normalized frequency, the model fit parameters 𝓐_FSS and 𝓢_FSS, the analytic value 𝓢_asymp, and the relative difference between the latter.

View Table | View all tables in this article

The next example is for hexagonal lattices, for which we consider both rod-type structures in TM (E _∥) polarization, and hole type structures in TE (H _∥) polarization. The geometry for TM polarization comprises rods of refractive index n_c = 3 and normalized radius a/d = 0.2 in a background of refractive index n_b = 1, while for the TE polarization, we invert the lattice so that n_c = 1, n_b = 3, preserving the value of a/d = 0.2. Fig. 3 shows band diagrams for (a) TM and (b) TE polarizations, which respectively exhibit band gaps for $\tilde{\omega }$ ∈ [0.31312,0.480299] and $\tilde{\omega }$ ∈ [0.225553,0.237806], respectively. Fig. 4(a) is similar to Fig. 2(b), but is for the hexagonal rod type lattice with TM polarization. To show the quality of the agreement between the analytic results and numerics we show in Fig. 4(b) the same results (for n_d > n_c) on a logarithmic vertical scale: the curves overlap for a normalized frequency difference (off the band edge) of more than 4 orders of magnitude. Figs 2(c) and (d) are similar, but refer to the hole-type lattice with TE polarization. We only show data for n_d > n_c since for n_d < n_c we could not accurately find the defect mode numerically (see the discussion below). Note again the excellent agreement between the analytic and numerical results. This agreement is further illustrated in Table 2, in which the data is presented as in Table 1: in all cases the agreement between 𝓢_FSS and 𝓢_asymp is much better than 1% for all cases in Table 2.

Table 2. Similar to Table 1, but for the lowest band gap of the hexagonal lattices from Fig. 4 for TE and TM polarization.

View Table | View all tables in this article

Table 2 shows why we could not find defect modes in TE polarization for n_d < n_c using the FSS method: the computed value of the sensitivity ε_C0 𝓢_asymp ≈ 900, implying that the defect modes are extremely close to the band edge, and thus difficult to locate. In turn, the origin of this high sensitivity can be seen from Fig. 5 which shows the normalized electric energy distribution f(rthe bulk modes for the hexagonal lattice ) = ε∥E(r)∥²/𝓔 _WSC, so that ∫_WSC f(r)d ² r = 1. Figs 5(a) and 5(b) show contour maps of the normalized electric energy distribution (on a base 10 logarithmic scale) for the bulk mode, respectively at the lower and upper edge of the first band gap, while Figs 5(c) and 5(d) respectively show horizontal and vertical sections through the centre of the defect for modes at the lower (red) and upper (blue) edge of the band gap. Clearly evident from these is the very low energy density within the cylindrical inclusion for the mode at the lower edge of the band gap, thereby explaining the very high sensitivity factor (19).

 

Fig. 4. (a) Similar to Fig. 2(b), but for a hexagonal lattice with TM polarization in a rod type structure. (b) Replotting of the data in (a) showing |ω̃ -ω̃_L| for n_d > n_c for both the numerical (red dots) and the analytic (curve) results. Panels (c) and (d) are similar to (a) and (b) respectively, but for TE polarization and a hole-type structure.

Download Full Size | PDF

5. Discussion and conclusion

The analytic results (15) and (16) indicate that for high-curvature bands (|C| small; “effective mass” m* small), defect levels tend to remain close to the band edge. This is consistent with well known behaviour in solid state physics: for example, a shallow defect in a semiconductor can be described by a hydrogen atom model [10], in which the “ionization energy,” the energy of the defect’s ground state, is proportional to m*.

The theoretical approach we have outlined has relied on ideas presented by Economou [4] for two-dimensional periodicity, but has presented them in a very general framework. Economou’s discussion was in the context of the tight-binding model for the Schrödinger equation, whereas our treatment has been for the Helmholtz equation, with boundary conditions dependent on polarization. We have shown however that there is no requirement for employing a tight-binding model, and that the results may be cast in a polarization-independent form relying on two factors: the change in electric energy caused by a perturbation, and the density of states at the band edge. We stress that, although the numerical examples have concentrated on the case of defects introduced by changing a single cylinder in a periodic two-dimensional lattice, the results apply to defects caused by perturbing multiple cylinders in a localized region, and even to photonic crystals in which the dielectric constant and its perturbations are continuous functions of position. We also note that the quality of the agreement exhibited between the perturbation formula and the numerical results based on the fictitious source superposition method highlights the value of the FSS technique in dealing with physical situations involving states which, though spatially localized, are nonetheless highly extended.

 

Fig. 5. Normalized electric energy distribution (on a base 10 logarithmic scale) of the bulk modes for the hexagonal lattice (a/d = 0.2, ν_b = 3, ν_c = 1) in TE polarization respectively at the (a) lower and (b) upper edges of the band gap. Panels (c) and (d) are horizontal and vertical section plots through the centre of the defect for the normalized energy density for modes at the lower (red) and upper (blue) edge of the band gap.

Download Full Size | PDF

As a final point, we note that even though our examples in Section 4 referred to PCs with periodically placed cylindrical inclusions in an otherwise uniform background, our analytic result (16) can apply to any periodic refractive index distribution, capable of creating a complete band gap. The only restriction is that the region C ₀ where the perturbation is applied, is of finite spatial extent.