## Abstract

We consider defect modes created in complete gaps of 2D photonic crystals by perturbing the dielectric constant in some region. We study their evolution from a band edge with increasing perturbation using an asymptotic method that approximates the Green function by its dominant component which is associated with the bulk mode at the band edge. From this, we derive a simple exponential law which links the frequency difference between the defect mode and the band edge to the relative change in the electric energy. We present numerical results which demonstrate the accuracy of the exponential law, for TE and TM polarizations, hexagonal and square arrays, and in each of the first and second band gaps.

©2007 Optical Society of America

## 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*
^{2}(**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*_{0},

*)} which obey the Helmholtz equation, the Bloch condition associated with the Bloch vector*

**r**

**k**_{0}, 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*(

*) = 1,*

**r***s*(

*) =*

**r***ε*(

*), and, for TE (*

**r***H*

_{∥}) polarization, we introduce

*ψ*=

*H*,

_{z}*p*(

*) = 1/*

**r***ε*(

*),*

**r***s*(

*) = 1. In terms of these, the Helmholtz equation, for either polarization, can be written as*

**r**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**

_{0}=

**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**_{0}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*_{0}, 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

*, must thus be replaced by a sum over those*

**k**_{L}*, values that correspond to the band edge, with each contribution weighted by*

**k**_{L,j}*θ*=

_{L,j}*ϕ*/(2

_{L,j}*π*).

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*
_{0} to create the defect mode. The defect mode, which is localized and tends to zero exponentially with increasing distance from *C*
_{0}, satisfies the Helmholtz equation ▽ ∙ (*p*̃(* r*)▽

*ψ*(

*)) + (*

**r***ω*

^{2}/

*c*

^{2})

*s*̃(

*)*

**r***ψ*(

*) = 0, in which the perturbation is characterized by the quantities*

**r***p*̃(

*) and*

**r***s*̃(

*).*

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

$$+\frac{{\omega}^{2}}{{c}^{2}}{\int}_{{C}_{0}}\left(s\left(\mathrm{r\prime}\right)-\tilde{s}\left(\mathrm{r\prime}\right)\right)\psi \left(\mathrm{r\prime}\right)G(r,\mathrm{r\prime};\omega ){d}^{2}\mathrm{r\prime}.$$

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*̃(

*′) = 1/*

**r***ε*(

*′), and*

**r***s*(

*′) =*

**r***s*̃(

*′) = 1, and substituting into Eq. (10) the leading order estimate (9) for*

**r***G*, we obtain

$$\times {\int}_{{C}_{0}}\frac{\tilde{\epsilon}\left(r\prime \right)-\epsilon \left(r\prime \right)}{\tilde{\epsilon}\left(r\prime \right)\epsilon \left(r\prime \right)}\nabla \psi \left(r\prime \right)\bullet \nabla {\psi}_{L}^{*}({k}_{L,j}r\prime ){d}^{2}r\prime ,$$

where *j* sums over all values of **k**_{0} = * 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*

_{0}. Furthermore, for this case,

*θ*= 1/3,

_{L,j}*C*=

_{L,j}*C*, and

_{L}*M*=

_{L,j}*M*, for all

_{L}*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},

*), corresponding to a particular mode. The sum in Eq. (11) now contains 6 identical terms each of which are weighted by*

**r***θ*= 1/3. This leads us to define a geometry factor

_{L,j}*f*, for which here

_{L}*f*, = 6 × 1/3 = 2. Of course, for a singularity that is completely enclosed within the BZ,

_{L}*f*= 1.

_{L}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*
^{2}/*ω*
^{2}
_{L}, the result in Eq. (12) is the ratio of the change in the electric energy (∫* E*∙

**D***d*

^{2}

*) 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*

**r***C*

_{0}, 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*(2

_{L}A_{WSC}*π*) 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*̃(

*̃) = 1,*

**r***s*(

*′) = ε(*

**r***′) and*

**r***s*̃(

*′) = $\tilde{\epsilon}$ (*

**r***′) and, as for TE polarization, the substitution of these relations into Eq. (10) leads to*

**r**$$\phantom{\rule{.2em}{0ex}}\times {\int}_{{C}_{0}}\left(\tilde{\epsilon}\left(r\prime \right)-\epsilon \left(r\prime \right)\right)\psi \left(r\prime \right){\psi}_{L}^{*}({k}_{L,j},r\prime ){d}^{2}r\prime $$

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*. For a defect caused by changing the refractive index of a single cylinder to

_{c}*n*, we have δ𝓔

_{d}_{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.

## 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*= 1. Fig. 2 (a) shows the band diagram which has band gaps for

_{b}*ω*̃ ≡

*ω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*) of the single defect cylinder. For

_{d}*n*>

_{d}*n*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.

_{c}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.

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*= 1, while for the TE polarization, we invert the lattice so that

_{b}*n*= 1,

_{c}*n*= 3, preserving the value of

_{b}*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*) 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

_{c}*n*>

_{d}*n*since for

_{c}*n*<

_{d}*n*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 𝓢

_{c}_{FSS}and 𝓢

_{asymp}is much better than 1% for all cases in Table 2.

Table 2 shows why we could not find defect modes in TE polarization for *n _{d}* <

*n*using the FSS method: the computed value of the sensitivity ε

_{c}_{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*(

*) = ε∥*

**r**the bulk modes for the hexagonal lattice*E*(

*)∥*

**r**^{2}/

*𝓔*

_{WSC}, so that ∫

_{WSC}

*f*(

*)*

**r***d*

^{2}

*= 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).*

**r**## 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.

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*
_{0} where the perturbation is applied, is of finite spatial extent.

## Acknowledgments

The support of the Australian Research Council through its Centres of Excellence and Discovery Grants programs is acknowledged. We also acknowledge discussions with Prof. Costas Soukoulis that motivated this work, and with Prof. Harald Giessen.

## References and links

**1. **J. D. Joannopoulos, R. D. Meade, and J. N. Winn, *Photonic crystals: molding the flow of light* (Princeton University Press, Princeton, 1995).

**2. **M. Ibanescu, S. G. Johnson, D. Roundy, C. Luo, Y. Fink, and J. D. Joannopoulos, “Anomalous dispersion relations by symmetry breaking in axially uniform waveguides,” Phys. Rev. Lett. **92**, 063903 (2004). [CrossRef] [PubMed]

**3. **A. Figotin and I. Vitebskiy, “Slow light in photonic crystals,” Waves in Random and Complex Media **16**, 293–382 (2006). [CrossRef]

**4. **E. N. Economou, *Green’s functions in quantum physics*, 2nd ed. (Springer-Verlag, Berlin, 1983).

**5. **R. C. McPhedran, L. C. Botten, J. McOrist, A. A. Asatryan, C. M. de Sterke, and N. A. Nicorovici, “Density of states functions for photonic crystals,” Phys. Rev. E **69**, 016609 (2004). [CrossRef]

**6. **S. Wilcox, L. C. Botten, R. C. McPhedran, C. G. Poulton, and C. M. de Sterke, “Modeling of defect modes in photonic crystals using the fictitious source superposition method,” Phys. Rev. E **71**, 056606 (2005). [CrossRef]

**7. **L. C. Botten, K. B. Dossou, S. Wilcox, R. C. McPhedran, C. M. de Sterke, N. A. Nicorovici, and A. A. Asatryan, “Highly accurate modelling of generalized defect modes in photonic crystals using the FSS method,” Int. J. Microwave and Optical Technology **1**, 133–145 (2006).

**8. **D. P. Fussell, R. C. McPhedran, and C. M. de Sterke, “Two-dimensional treatment of the level shift and decay rate in photonic crystals,” Phys. Rev. E **72**, 046605 (2005). [CrossRef]

**9. **K. Dossou, M. A. Byrne, and L. C. Botten, “Finite element computation of grating scattering matrices and application to photonic crystal band calculations,” J. Comput. Phys. **219**, 120–143 (2006). [CrossRef]

**10. **W. Kohn and J. Luttinger, “Theory of donor states in silicon,” Phys. Rev. **98**, 915 – 922 (1955). [CrossRef]