## Abstract

We investigate the formation of photonic crystal waveguide (PCW) modes within the framework of perturbation theory. We derive a differential equation governing the envelope of PCW modes constructed from weak perturbations using an effective mass formulation based on the Luttinger-Kohn method from solid-state physics. The solution of this equation gives the frequency of the mode and its field. The differential equation lends itself to simple analytic approximations which reduce the problem to that of solving slab waveguide modes. By using this model, we demonstrate that the nature of the projected band structure and corresponding Bloch functions are central to the behaviour of PCW modes. With this understanding, we explain why the odd mode in a hexagonal PCW spans the entire Brillouin zone while the even mode is cut off.

©2009 Optical Society of America

## 1. Introduction

Photonic crystal waveguides (PCWs) enable control over many optical properties such as dispersion [1, 2] and slow light [3, 4]. These properties have been exploited to enable the creation of devices such as temperature-tuned slow light PCWs [5], optical delay lines [6] and more recently to enhance nonlinear effects such as third harmonic generation [7]. In such applications, the linear properties of PCWs are typically obtained using purely numerical methods. These algorithms accurately calculate the frequencies and fields of PCW modes, and are often employed because the problem is considered too complicated to admit an analytic treatment. However, in solid-state physics a range of perturbative methods have been developed to give physical insight into localised modes in bandgaps [8]. Perhaps the most well-known method is the effective mass approximation (EMA) as introduced by Luttinger and Kohn (LK) [9]. This method has been applied to investigate solid state phenomena including the behaviour of donor states in silicon [10] and the formation of superlattice band structures [11].

PCWs are typically constructed by removing a row of inclusions, creating a large perturbation, and therefore the frequencies of their modes lie deep within the bandgap. Consequently, their behaviour is typically analysed under these conditions. For example, Notomi *et al*. showed that the shape of dispersion curves of even PCW modes is due to an anti-crossing [12]. However, this only explains the shape of the dispersion curve of even modes when they are deep in the gap. Building on this, Petrov and Eich [1] qualitatively explained the shape of both odd and even modes in hexagonal lattices by the folding of the Brillouin zone (BZ) as a result of the introduction of a PCW. PCWs are perhaps experimentally most useful in the deep-gap regime where their modes are tightly bound. However, in this paper we show that the nature of PCW modes can be alternatively explained when traced back to the shallow limit.

Using the LK framework, we derive a differential equation governing the field and frequency of PCW modes, where the PCW is constructed by a weak periodic perturbation in a two-dimensional (2D) photonic crystal (PC). Crucially, the differential equation is written in terms of the envelope function of the PCW mode rather than the modal fields themselves; this results in the equation being of the same form as that governing a simple slab waveguide structure. By making further approximations, we show that the differential equation admits analytic solutions for a PCW constructed by altering a row of inclusions. The results obtained using our analytic formulation are in excellent agreement with fully numerical calculations and work particularly well when modes are near the band-edge. While the method is applied here to PCWs, the EMA is quite general and has the potential to provide insight into a variety of photonic crystal devices based on defects.

Previously, it has been shown that modes of point defects evolve from band extrema and enter the band gap; that is, their fields resemble the band-edge Bloch mode modulated by an envelope function [13]. For PCWs, periodicity along the direction of propagation is preserved and thus their modes evolve from *projected band-edges* rather than band extrema. The projected band-edge for a PCW along the *x* direction is defined by finding the frequency extremum for each value of *k _{x}* [14, 15]. This is illustrated in Fig. 1, in which a PCW is constructed in a hexagonal lattice along the Γ-K direction, which we take to be the

*k*axis. As this is at the top of the bandgap we look for frequency minima. At

_{x}*k*=0 the minimum frequency is at the BZ edge (

_{x}d*k*=2

_{y}d*π*/√3), since the global band minimum is at the M point. At

*k*~1.956 the minimum moves inside the BZ along two inequivalent trajectories, where the second trajectory is the reflection of the first about

_{x}d*k*=2

_{y}d*π*/√3 as shown Fig. 1. The path of the frequency minima within the BZ gives the

*projected band-edge*trajectory. PCWs in hexagonal lattices are typically constructed by increasing the amount of dielectric, and thus form at the top of the gap. We show in this paper that the formation of an odd and even mode in hexagonal PCWs results from coupling between modes associated with the two inequivalent trajectories shown in Fig. 1. We use this to explain why the odd mode spans the entire BZ, while the even mode is cut off.

The outline of this paper is as follows: Section 2 contains the derivation of the envelope differential equation for PCWs in the simplest (non-degenerate) case in a form amenable to numerical solution. In Section 3 we obtain analytic solutions to the differential equation for a simple instance, namely, a square PC in TM (*E _{z}*) polarisation. Section 4 deals with the PCWs in hexagonal lattices and outlines how the introduction of the PCW causes the coupling of degenerate band-edge trajectories. We discuss how this determines the formation of odd and even PCW modes in different parts of the BZ.

## 2. Formulation

We commence by writing down the polarisation independent wave equation governing an infinite 2D PC, composed of lossless dielectric, in the *x*-*y* plane:

where ε_{PC}(**r**) denotes the permittivity of the photonic crystal. Here, *ω*
** _{k}** is the frequency of the mode with Bloch wavevector

**k**and Bloch mode

**E**(

_{k}**r**). We construct our PCW such that it is a periodic perturbation oriented along

*x*(see Fig. 2). The wave equation governing the modes of the PCW is

Here, *ω _{kx}* and

**E**

*(*

_{kx}**r**) are the frequency and fields of the PCW mode with associated Bloch wavevector

*k*indicating that the modes of the PCW are Bloch modes due to periodicity along

_{x}*x*. We have expressed the perturbation as an infinite sum over each of its periods. We show below that this reduces the problem to a single period of the PCW. We note that although the columns are of width Λ=

*d*, more general perturbations of size Λ=

*md*, where

*m*is a positive integer forming coupled-resonator optical waveguides [17] can also be made.

We choose to expand the eigenfunctions of the PCW, namely **E**
* _{kx}* (

**r**), in terms of the Bloch modes of the infinite 2D PC. We use an extended zone scheme, whereby all the PC bands are folded out such that we can use a single wavevector

**k**ranging over all of reciprocal space to specify a point in the band-structure uniquely (see Appendix of [9]). The ansatz takes the form

where the integral is over all of *k _{y}*. Our ansatz implies that we look for PCW mode solutions in the form of Bloch modes modulated by an envelope function in the direction transverse to their propagation. The function

*C*(

*k*) tells us how we superpose Bloch modes in reciprocal space to construct the PCW mode, ${\mathbf{E}}_{{k}_{x}}\left(\mathbf{r}\right)$. We now substitute Eq. (3) into Eq. (2) making use of Eq. (1). Taking the inner product with

_{y}**E**

^{*}

_{k′}and utilising the normalisation condition, ∫

_{ℝ2}

*ε*(

_{PC}**r**)

**E**

^{*}

_{k′}(

**r**) ·

**E**(

_{k}**r**)

*d*

^{2}

**r**=

*Mδ*(

**k**′-

**k**), where

*M*is a normalisation constant, gives the effective-mass-like equation:

Here, the normalisation constant is *M*=(2*π*)^{2}
𝓔/Ω, defined in terms of the electric energy over the unit cell, $\U0001d4d4={\int}_{\mathrm{cell}}{\epsilon}_{\mathit{PC}}\left(\mathbf{r}\right){\parallel {\mathbf{E}}_{{k}_{x},{k}_{{L}_{y}}}\left(\mathbf{r}\right)\parallel}^{2}{d}^{2}\mathbf{r}$, and the area of the unit cell, Ω. Eq. (4) thus shows how the introduction of a PCW couples Bloch modes.

Since we are deriving a band-edge theory, we are looking at perturbations where *δε _{p}*(

**r**) is shallow, or equivalently, the change in electric energy in each period of the PCW can be treated as a shallow perturbation. The function

*C*(

*k*) is then highly localised in reciprocal space because the PCW modes are loosely bound. The dominant Bloch contribution in Eq. (4) is then the mode with Bloch wavevector corresponding to the peak of

_{y}*C*(

*k*). If there are multiple points in the BZ which have the same frequency at a given value of

_{y}*k*(for example

_{x}d*k*>1.956 in Fig. 1), the function

_{x}d*C*(

*k*) has multiple narrow peaks in reciprocal space. Both single and multiple peak cases arise in PCs with typical parameters; however, in this section we assume that the perturbation does not cause any coupling between these peaks and leave the analysis of the coupled case for Section 4.

_{y}We note here that a similar treatment was previously reported by Charbonneau-Lefort *et al* [16], but that their study made the additional assumption that the perturbation, that is, the change in dielectric constant, did not vary rapidly over a unit cell. Since typical PC perturbations contain abrupt changes in dielectric constant, we do not make this assumption here. We however make an alternative assumption: recalling that *C*(*k _{y}*) is a highly localised function, the Bloch modes

**E**(

_{k}**r**) do not change much in the region of reciprocal space where

*C*(

*k*) is large. Therefore, for small perturbations, Eq. (4) is dominated by Bloch modes with

_{y}*k*~

_{y}*k*, where

_{Ly}*k*is the point in reciprocal space where

_{Ly}*C*(

*k*) is centered. Bloch modes have the form ${\mathbf{E}}_{\mathbf{k}}\left(\mathbf{r}\right)={\mathbf{u}}_{{k}_{x},{k}_{y}}\left(\mathbf{r}\right){e}^{i\mathbf{k}\xb7\mathbf{r}}$ (where ${\mathbf{u}}_{{k}_{x},{k}_{y}}\left(\mathbf{r}\right)$ is a periodic function with the period of the lattice) and we thus approximate ${\mathbf{u}}_{{k}_{x},{k}_{y}}\sim {\mathbf{u}}_{{k}_{x},{k}_{\mathit{Ly}}}$. We then immediately write ${\mathbf{E}}_{\mathbf{k}}\left(\mathbf{r}\right)={\mathbf{E}}_{{k}_{x},{k}_{\mathit{Ly}}}\left(\mathbf{r}\right){e}^{i\left({k}_{y}-{k}_{\mathit{Ly}}\right)y}$. The equivalent replacement for

_{y}**E**′ (

_{k}**r**) in Eq. (4) is also made, as both sides of the equality are of appreciable magnitude only when

*k*

^{′}

_{y}is near

*k*, and therefore these points contribute most to the PCW mode.

_{Ly}Our goal is to obtain a differential equation satisfied by the envelope function in real space so we must take the inverse Fourier transform of our equation. We thus multiply by exp[*i*(*k*′_{y}-*k*
_{Ly})y′] and integrate over *k*
^{′}
* _{y}*. The right hand side (RHS) of Eq. (4) is then

We now introduce the function

where the integral is again over all *k _{y}*. We show below that

*f*(

*y*) is the envelope function of the PCW mode in the direction perpendicular to its propagation. Using this definition we obtain for Expression (5),

Using Eq. (6), the left hand side (LHS) of Eq. (4), after multiplication by exp[*i*(*k*′_{y}-*k*
_{Ly})y′] and integration over all *k*
^{′}
* _{y}*, reads

Here, ${\omega}_{{k}_{x},-\mathit{id}\u2044\mathit{dy}\prime}^{2}$ is an operator, expanded about the frequency extremum in a Taylor series. This operates on *f* (*y*
^{′}) and its mixed subscript is mathematically a result of inverse Fourier transforming with respect to *k _{y}* but not

*k*. Physically, this is because the PCW is periodic in the

_{x}*x*direction, but not in

*y*, and therefore the modes propagate with specified kx. The confinement of the modes in the

*y*direction means a mode is associated with a range of

*k*and hence the differential equation is for the envelope function in the

_{y}*y*direction. The operator thus expands about a frequency extremum at a given

*k*in powers of -

_{x}*id*/

*dy*

^{′}. This frequency extremum is simply the projected band-edge. The expansion in powers of -

*id*/

*dy*

^{′}is related to how the band surface curves in the

*k*direction along the projected band-edge trajectory. Since

_{y}*f*(

*y*) is a slowly varying envelope function, we can take only the leading order term in -

*id*/

*dy*

^{′}. Typically, this means the operator is of the form -

*d*

^{2}/

*d*

_{y}^{′}

^{2}, however in Section 4 we show that points exist where higher order derivatives need to be included.

Throughout, we have kept the quasi-periodicity of the Bloch modes in the *x* direction, which we now use to simplify the summation in Eq. (7). We proceed by replacing the integral over all space with an integral of the Bloch modes over a single period of the perturbation multiplied by an infinite sum over the phases of the Bloch modes over each period of the perturbation. Using the Poisson summation formula, the integral in Eq. (7) becomes

Dropping the primes, we equate the LHS and RHS of our equation giving a differential equation governing the PCW mode,

where again we note that *ω _{kx}* is the frequency of the PCW mode. Here the integral is over those regions for which

*δε*(

**r**) is non-zero. We recall that ${\mathbf{E}}_{{k}_{x},{k}_{\mathit{Ly}}}\left(\mathbf{r}\right)$ is the electric field of the Bloch mode at

*k*along the projected band-edge.

_{x}Finally the PCW mode is obtained from our ansatz Eq. (3), by again making the approximation ${u}_{\mathbf{k}}\left(\mathbf{r}\right)\sim {u}_{{k}_{x},{k}_{\mathit{Ly}}}\left(\mathbf{r}\right)$, giving ${\mathbf{E}}_{{k}_{x}}\left(\mathbf{r}\right)=f\left(y\right){\mathbf{E}}_{{k}_{x},{k}_{\mathit{Ly}}}\left(\mathbf{r}\right)$. Thus, *f*(*y*) is the envelope function of the PCW mode, which is expressed as a band-edge Bloch mode modulated by *f*(*y*).

## 3. Solution for TM polarisation in a square lattice

Equation (10) can be solved using a variety of numerical methods for arbitrarily shaped shallow PCWs. This has several advantages over solving the wave equation directly using finite element or finite difference methods. In our approach, calculating the modes of the periodic structure is separated from finding the modes of the perturbed structure. This means that we first calculate the Bloch modes of the 2D PC, and with these, we solve Eq. (10) to obtain modes of the PCW. We also note that Eq. (10) can be solved in a similar fashion to that for a symmetric graded index slab waveguide. This is less numerically intensive than solving the full wave equation directly.

In this section, we solve Eq. (10) analytically for a PCW constructed by uniformly altering the refractive index of a row of cylindrical inclusions in a square lattice PC along the *x* direction as shown in Fig. 2(a). This leads to a uniform change in permittivity *δε* in a row of cylinders. The parameters of the PC we use for numerical examples are: cylinder index *n _{c}*=3, background index

*n*=1 and cylinder radius

_{b}*a*=0.3

*d*in TM (

*E*) polarisation. We begin by writing down our operator along a projected band-edge trajectory, which, to first order, is quadratic in the

_{z}*k*direction:

_{y}$${\omega}_{{k}_{x},-i\genfrac{}{}{0.1ex}{}{d}{\mathit{dy}}}^{2}=\left[{\omega}^{2}({k}_{x},{k}_{\mathit{Ly}})-\genfrac{}{}{0.1ex}{}{\omega ({k}_{x},{k}_{\mathit{Ly}})}{{D}_{\mathit{Ly}}\left({k}_{x}\right)}\genfrac{}{}{0.1ex}{}{{d}^{2}}{{\mathit{dy}}^{2}}+\dots \right].$$

Here, *ω*(*k*
_{x},*k*
_{Ly}) denotes the frequency along the projected band-edge and *D*
_{Ly}(*k*
_{x}) is the curvature in the *k _{y}* direction along the projected band-edge trajectory. We construct the PCW by decreasing the index of a row of inclusions and therefore its mode appears at the bottom of the bandgap. The projected band-edge trajectory for the first projected band-edge ranges from the Y to the M point in the BZ and hence is non-degenerate. Our perturbations are shallow and therefore the envelope function

*f*(

*y*) is slowly varying with respect to the lattice period. Thus, when operating with ${\omega}_{{k}_{x},-\mathit{id}\u2044\mathit{dy}}^{2}$ we keep only lowest order terms in the derivative. Equation (10) now becomes

where *C* indicates the integral is over the domain of the perturbation, namely the region where *δε* is non-zero and uniform. We cannot yet solve this equation analytically because of the term containing Bloch modes. The typical shape of this term is shown by the red curve in Fig. 3. Since we are pursuing an analytic solution of this equation, we wish to replace this term with one which is piecewise continuous, as shown by the dashed line in Fig. 3. Given that our formulation is based on first order perturbation theory, the quantity governing the frequency of the PCW modes is the total change in electric energy resulting from the introduction of the PCW. Therefore, if we alter the shape of the change in electric energy as a function of y while preserving its total magnitude, we can solve for the frequency of the PCW mode exactly. We then make the following replacement: ${\int}_{C}\delta \epsilon {\mid \mid {\mathbf{E}}_{{k}_{x},{k}_{\mathit{Ly}}}\left(\mathbf{r}\right)\mid \mid}^{2}\mathit{dx}\to 1\u2044\left(2a\right){\int}_{C\prime}\delta \epsilon {\mid \mid {\mathbf{E}}_{{k}_{x},{k}_{\mathit{Ly}}}\left(\mathbf{r}\right)\mid \mid}^{2}\mathit{dx}\mathit{dy}$, where *C*′ now indicates the perturbed domain in two dimensions. The inclusion of the integral over *dy* and division by 2*a* ensures that the change in magnitude of the perturbation is unaltered. The two curves in Fig. 3 thus show the form of this term before and after the alteration. As the area under the two curves is the same they lead to the same total change in electric energy. We can now solve for the frequency analytically. The cost of this simplification is that the resultant envelope function and corresponding fields inside the perturbation are slightly altered. This effect is almost irrelevant for weakly perturbed PCWs since most of the energy of modes is outside the perturbation. Equation (12) then becomes

where $\delta {\U0001d4d4}_{{k}_{x}}={\int}_{C\prime}\delta \epsilon \mathit{dx}\mathit{dy}{\parallel {\mathbf{E}}_{{k}_{x},{k}_{\mathit{Ly}}}\left(\mathbf{r}\right)\parallel}^{2}$. This equation has an internal solution (where $\delta {\U0001d4d4}_{{k}_{x}}\ne 0$) and an external solution (where $\delta {\U0001d4d4}_{{k}_{x}}=0$). In both cases, Eq. (13) is of the form of a one-dimensional Helmholtz equation and is solved in the same fashion as a slab waveguide or potential well. We begin with the external solution, where the equation simplifies to

with solution f(y)=Ae-α|y|.

Here, *α*
^{2} is always positive, implying that guided solutions are obtained for modes evolving from both the top and bottom gap-edges; that is, bound PCW modes exist when increasing or decreasing the index of a row of inclusions. The internal solution is

with solutions *f*(*y*)=*B*cos(*β _{y}*) (even solution).

For shallow perturbations there exists only a single PCW mode with an even envelope. Odd envelope functions are associated with higher order modes due to strong perturbations and are thus not discussed here. The appropriate boundary condition for both polarisations requires that the envelope function and its derivative be continuous at the boundary of the perturbation, giving the transcendental equation

An explicit form is obtained by taking a small argument approximation of tan(*βa*) and realising that, typically in perturbation theory, the magnitude of the perturbation is larger than the corresponding shift in modal frequency and thus the second term in *β*
^{2} is of leading order. With the additional approximation *ω*
^{2}
_{kx}/*ω*(*k*
_{x},*k*
_{Ly})~*ω*(*k*
_{x},*k*
_{Ly}), this gives

We have thus obtained a simple expression for the frequency of a PCW mode, analytically expressing how the mode evolves from the band-edge. The quadratic dependence of the frequency on the perturbation is typical of 1D perturbation theory. The PCW has a fundamental mode when $\delta {\U0001d4d4}_{{k}_{x}}\ne 0$ and the mode emerges from the projected band-edge at all *k _{x}* values. This is shown for a PCW created by decreasing the index of a row of cylinders in Fig. 4(a). Here we observe that the PCW mode “peels” off the projected band structure – by this we mean that the dispersion curve of the PCW mode has a similar shape to the projected band-edge but is shifted upward in frequency. The magnitude of this shift is different for each

*k*and is governed by Eq. (17). The agreement between fully numerical calculations and our analytic theory is shown in Fig. 4(b). The numerical method used throughout this paper employs Bloch modes to calculate the reflection matrix of a half plane above and below the PCW as shown by de Sterke

_{x}*et al*[20]. As expected the agreement is best when the change in frequency from the projected band-edge is smallest. Fig. 5 shows a comparison of numerical and analytic calculations for the

*E*field of a PCW mode constructed by altering the refractive index of a row of cylinders from

_{z}*n*=3.0 to

_{c}*n*=2.9. Fig. 5(a) shows a contour plot of the waveguide mode constructed by modulating a Bloch mode by the envelope function calculated from Eqs. (14) and (15). The result was indistinguishable from a contour plot calculated numerically. This is highlighted in Figs. 5(b) and 5(c) where we compare |

_{w}*E*| calculated numerically and using our analytic method at

_{z}*d*/

*λ*=0.245 and

*d*/

*λ*=0.266. The agreement is better for the former, since from Fig. 4(b) we observe that asymptotic model is more accurate when Δ(

*d*/

*λ*) is smaller.

## 4. Degeneracy and solutions in hexagonal lattices for TE polarisation

In Sections 1 and 2 we mentioned that the introduction of a PCW breaks the symmetry of the PC and can lead to coupling between degenerate projected band-edge trajectories. In this section we illustrate that these types of degeneracies are central to understanding the evolution of PCWs in hexagonal PCs. Using our formulation we show that the coupling of these degenerate points leads to the formation of two modes typically observed in hexagonal PCWs. In particular we explain why the odd mode spans the entire BZ while the even mode is cut off.

We now consider the modes of PCWs constructed in hexagonal PCs, as shown in Fig. 2(b). Figure 1 shows the projected band-edge trajectory for the second band of a hexagonal PC with, cylinder index *n _{c}*=1, background index

*n*=3 and cylinder radius

_{b}*a*=0.3

*d*for TE (

*H*) polarisation. The projected band-edge trajectory moves along the BZ-edge and at

_{z}*k*~1.956 the trajectory moves into the BZ. This point is a supercritical pitchfork bifurcation due to a single minimum turning into two minima above and below the BZ-edge as well as a maximum on the BZ edge. When

_{x}d*k*<1.956 there is a single non-degenerate band-edge trajectory and thus the PCW modes are obtained as in Section 3. The solution for the frequency of PCW modes where the index of a row of air holes is changed to

_{x}d*n*=1.05 is shown in Fig. 6(a). Here we have plotted the frequency difference between the PCW mode and the projected band-edge. The frequency difference is small and thus our asymptotic theory works well. When

_{w}*k*>1.956, there are two projected band trajectories which have the same frequency and are thus degenerate. The existence of a degeneracy means the degenerate modes have the potential to couple under a perturbation. Therefore, the introduction of the PCW can cause coupling of the Bloch modes along the two projected band-edge trajectories. We note that although the band-structure calculated in Fig. 1 is for a particular set of parameters, the occurrence of the bifurcation point along M-K is a common feature of the second band in hexagonal PCs for TE polarisation. The band surface for any typical parameter set has a band minimum at the M point, while the K point is a local maximum as a function of

_{x}d*k*([14] pg. 76 and pg. 137). With these two conditions, a bifurcation point of the form shown in Fig. 1 must exist along M-K.

_{y}The effect of the coupling of degenerate points in the BZ has previously been analysed for point defects in PCs [18] and is typical of degenerate perturbation theory [19]. We now show that the treatment derived in Section 2 can be adapted to the degenerate case.

We commence our formulation from Eq. (4) in Section 2 outlining the coupling of points in reciprocal space due to the introduction of a PCW. As there are two projected band-edge trajectories for a given *k _{x}*, the Fourier transform of the envelope function has narrow peaks centred on each of the trajectories. Provided we are not too close to the bifurcation point we may write

*C*(

*k*

_{y})=

*C*

_{1}(

*k*

_{y}-

*k*

^{(1)}

_{Ly})+

*C*

_{2}(

*k*

_{y}-

*k*

^{(2)}

_{Ly}), where

*C*

_{1}(

*k*

_{y}-

*k*

^{(1)}

_{Ly}) and

*C*

_{2}(

*k*

_{y}-

*k*

^{(2)}

_{Ly}) are the functions associated with the two peaks centered on the upper and lower trajectory respectively. Re-writing Eq. (4) in this form gives

$$={\omega}_{{k}_{x}}^{2}\int d{k}_{y}\left({C}_{1}\left({k}_{y}-{k}_{\mathit{Ly}}^{\left(1\right)}\right)+{C}_{2}\left({k}_{y}-{k}_{\mathit{Ly}}^{\left(2\right)}\right)\right)\int \sum _{p}\delta {\epsilon}_{p}{\mathbf{E}}_{\mathbf{k}\prime}^{*}\left(\mathbf{r}\right)\xb7{\mathbf{E}}_{k}\left(\mathbf{r}\right){d}^{2}\mathbf{r}.$$

Noting that *C*
_{1}(*k*
_{y}-*k*
^{(1)}
_{Ly}) and *C*
_{2}(*k*
_{y}-*k*
^{(2)}
_{Ly}) are narrow features in reciprocal space, when one is large the other is small. Therefore we separate Eq. (18) into two equations

$$={\omega}_{{k}_{x}}^{2}\int \left({C}_{1}\left({k}_{y}-{k}_{\mathit{Ly}}^{\left(1\right)}\right)+{C}_{2}\left({k}_{y}-{k}_{\mathit{Ly}}^{\left(2\right)}\right)\right){\mathit{dk}}_{y}\int \sum _{p}\delta {\epsilon}_{p}{\mathbf{E}}_{\mathbf{k}\prime}^{*}\left(\mathbf{r}\right).{\mathbf{E}}_{\mathbf{k}}\left(\mathbf{r}\right){d}^{2}\mathbf{r},$$

$$M\left[{\omega}_{{k}_{x},{k\prime}_{y}}^{2}-{\omega}_{{k}_{x}}^{2}\right]{C}_{2}\left({k\prime}_{y}-{k}_{\mathit{Ly}}^{\left(2\right)}\right)\delta \left({k}_{x}-{k\prime}_{x}\right)$$

$$={\omega}_{{k}_{x}}^{2}\int \left({C}_{1}\left({k}_{y}-{k}_{\mathit{Ly}}^{\left(1\right)}\right)+{C}_{2}\left({k}_{y}-{k}_{\mathit{Ly}}^{\left(2\right)}\right)\right){\mathit{dk}}_{y}\int \sum _{p}\delta {\epsilon}_{p}{\mathbf{E}}_{\mathbf{k}\prime}^{*}\left(\mathbf{r}\right)\xb7{\mathbf{E}}_{\mathbf{k}}\left(\mathbf{r}\right){d}^{2}\mathbf{r}.$$

We then proceed with each equation as in Sections 2 and 3, obtaining the matrix equation

where we have used the notation in Section 3, and where *δ𝓔*
_{12} refers to the overlap integral over the perturbation of Bloch modes with wavevectors corresponding to points *k*
^{(1)}
_{Ly} and *k*
^{(2)}
_{Ly} respectively. Here, the operators and the integrals in the diagonal elements are identical due to reflection symmetry about the BZ edge and the PCW centre, that is, *δ𝓔*
_{11}=*δ𝓔*
_{22}. We assume the perturbation is lossless and thus the off-diagonal terms are in general non-zero and are complex conjugates of each other, *δ𝓔*
_{12}=*δ𝓔*
^{*}
_{21}. The matrix is then Hermitian and is diagonalised to give two equations

$$\left[\left({\omega}_{{k}_{x}}^{2}-{\omega}_{{k}_{x},-i\genfrac{}{}{0.1ex}{}{d}{\mathit{dy}}}^{2}\right)+\genfrac{}{}{0.1ex}{}{{\omega}_{{k}_{x}}^{2}\mathrm{\Omega}}{\U0001d4d4\mathrm{\Lambda}}\delta {\U0001d4d4}_{\mathit{BB}}\right]{f}_{B}\left(y\right)=0.$$

Here, matrix diagonalisation is equivalent to choosing the Bloch mode superpositions such that they are orthogonal to each other over the perturbation. Therefore, *δ𝓔 _{AA}* and

*δ𝓔*correspond to the change in electric energy associated with the two Bloch modes. Equations (21) governing

_{BB}*f*(

_{A}*y*) and

*f*(

_{B}*y*) are now uncoupled giving envelope functions and frequencies of the two PCW modes. The frequency difference between the PCW mode and the projected band-edge is shown in Fig. 6(b). The asymptotic theory agrees better with numerical results when Δ(

*d*/

*λ*) is smaller.

We have shown that the introduction of the PCW breaks the degeneracy associated with the projected band-edge trajectory and that two modes form. Equation (21) gives the envelope functions associated with the two PCW modes. Recalling from Section 3 that the envelope functions are even, we conclude that the underlying Bloch modes determine the symmetry of the PCW modes. In diagonalising the matrix in Eq. (20) to obtain (21) we have chosen the Bloch mode superpositions such that the off-diagonal terms vanish. Since the PCW is symmetric about its centre, its modes must be either even or odd, and thus the superposition of Bloch modes that diagonalise the matrix in Eq. (20) must be those that are even and odd about the waveguide centre. Such linear combinations are shown in Figs. 7(b) and 7(c) for *k _{x}d*=

*π*. For 1.956 <

*k*<

_{x}d*π*the projected band-edge trajectory is inside the BZ, the Bloch modes are complex, and hence the appropriate superposition of Bloch modes changes. However, the diagonalisation condition remains unchanged and the modes remain either odd and even. We have thus established that the two envelope equations associated with the degenerate section of the BZ (

*k*> 1.956) yield an even and an odd PCW mode. Fig. 7(a) shows the Bloch mode at

_{x}d*k*=0, which is clearly odd symmetric about

_{x}d*x*; hence the PCW mode at this point is also odd symmetric. Since the projected band-edge trajectory has two paths at the BZ-edge and a single path at the BZ centre (see Fig. 1), we expect one of the PCW modes to cut off at the bifurcation point. Since the Bloch mode at the BZ centre is odd symmetric then the even mode must cutoff at the bifurcation point while the dispersion curve of the odd mode spans the entire BZ. Fig. 8 shows the dispersion curves in the vicinity of the bifurcation point for waveguide cylinder index

*n*=1.05. Our theory predicts that, for shallow perturbations, the even mode cuts off at the bifurcation point as the two projected band-edge trajectories merge. Analytically we know that the electric energy term must vanish for one of the equations in (21), however we note that in the vicinity of the bifurcation point the energy term is numerically sensitive to the value of kx and is often not exactly zero when computed.

_{w}## 5. Discussion and conclusion

The theory in Section 4 makes assumptions that do not hold close to the bifurcation point (*k _{x}d*~1.956). Accordingly, we now discuss the physics of the PCW mode in the vicinity of this point. In the above, we have written

*C*(

*k*) as a superposition of

_{y}*C*

_{1}(

*k*

_{y}-

*k*

^{(1)}

_{Ly}) and

*C*

_{2}(

*k*

_{y}-

*k*

^{(2)}

_{Ly}) and have postulated that when one is large the other is small. This assumption clearly breaks down as we get sufficiently close to the bifurcation point, for which the spacing between upper and lower projected band-edge trajectories becomes comparable to the width of

*C*

_{1}(

*k*

_{y}-

*k*

^{(1)}

_{Ly}) and

*C*

_{2}(

*k*

_{y}-

*k*

^{(2)}

_{Ly}). Furthermore, a stronger assumption is that we have used an effective mass approximation, that is, we have taken the band-edge to be quadratic to first order. At the bifurcation point the band-edge is quartic to first order. This is because the curvature goes from being positive in

*k*at M to being negative at K. Since the band-edge is symmetric in

_{y}*k*about the M-K line the frequency of the band-surface must scale quartically in

_{y}*k*at the bifurcation point. This means that the curvature term

_{y}*D*

_{Ly}(

*k*

_{x}), diverges. Thus from Eq. (17), [${\omega}_{{k}_{x}}-\omega ({k}_{x},{k}_{\mathit{Ly}})$] diverges as shown by the behaviour of the dashed green and purple curves in Fig. 8. The assumption that the band-edge is quadratic to leading order is not necessary in our formulation and we can choose to keep higher order terms in

*d*/

*dy*. Analytic solutions potentially exist when keeping terms of order

*d*

^{4}/

*dy*

^{4}, but are beyond the scope of this paper. The effective mass approximation thus breaks down when approaching the bifurcation point, either from the left or from the right. However, we assert that this region accounts for only a small fraction of the BZ and from Figs. 6(a), 6(b) and 8, we observe our results retain accuracy very close to the bifurcation point. Furthermore, the symmetry arguments we have made based on odd and even modes hold at this point.

In this paper we have only shown examples of PCWs constructed by altering the index of a row of cylinders, however, our treatment works equally for PCWs created by constructing a periodic perturbation of any form, such as by altering the shape, position or size of inclusions periodically. Moreover, although we derived our formulation for PCWs in 2D PCs, the only assumptions we have made about the unperturbed structure is that it be periodic and lossless, and it admits a complete set of Bloch mode solutions. Since our formulation is polarisation independent we can apply it to photonic crystal slabs. The envelope function is then the in-plane envelope function, while out-of-plane confinement due to the slab is unaffected by small perturbations.

We complement the previous figures showing the mode evolution with an online animation (Fig. 9) demonstrating the changes in the odd and even modes as the refractive index of the PCW cylinders increases from *n _{w}*=1.1 to

*n*=3.0. From this, we see that the even mode has bigger frequency shifts than the odd mode throughout this index range, and so our asymptotic estimate for its frequency is less accurate than for the odd mode. While the estimate breaks down for the odd mode at the bifurcation point, the cutoff of the even mode is accurately predicted when the perturbation is small. Nevertheless, our asymptotic model incorrectly predicts that the cutoff point for this mode does not change as the perturbation strength increases. This is because in our analytic formulation, cutoff occurs when the overlap of the projected band-edge Bloch mode with the perturbation is zero, which, for a uniform perturbation, occurs at a fixed value of

_{w}*k*.

_{x}dOur asymptotic theory is based primarily on band-edge Bloch modes and the topology of projected band-edges. Therefore, the primary restriction on the accuracy of our formulation for larger perturbations is that we only consider the contribution of the Bloch mode along the projected band-edge trajectory. We previously argued that this approximation is valid for modes of shallow PCWs since they are highly localised in reciprocal space; however, as PCW modes move deeper into the bandgap they become more tightly bound, and thus spread in reciprocal space. Hence, rather than approximating the integral in Eq. (4) by taking only the Bloch modes along the projected band-edge, more Bloch modes need to be taken into account. By taking a sufficient number of Bloch mode components it is possible to obtain solutions for the frequency and fields of modes of strongly perturbed PCWs. However, such solutions are not available analytically.

## 6. Acknowledgements

S. Mahmoodian would like to acknowledge P.Y. Chen andW.H. Dupree for their useful discussions. C.G. Poulton acknowledges support from a startup grant from University of Technology, Sydney. The support of the Australian Research Council through its Centres of Excellence and Discovery Grants Programs is acknowledged.

## Footnotes

^{1} | In fact we chose to take only the periodic function ${\mathbf{u}}_{{k}_{x},{k}_{\mathit{Ly}}}\left(\mathbf{r}\right)$ at this point rather than the Bloch mode itself. |

## References and links

**1. **A. Y. Petrov and M. Eich, “Zero dispersion at small group velocities in photonic crystal waveguides,” Appl. Phys. Lett. **85**, 4866–4868 (2004).
[CrossRef]

**2. **L. H. Frandsen, A. V. Lavrinenko, J. Fage-Pedersen, and P. I. Borel, “Photonic crystal waveguides with semi-slow light and tailored dispersion properties,” Opt. Express **14**, 9444–9450 (2006).
[CrossRef] [PubMed]

**3. **T. F. Krauss, “Slow light in photonic crystal waveguides,” J. Phys. D: Appl. Phys. **40**, 2666–2670 (2007).
[CrossRef]

**4. **J. Li, T. P. White, L. O’Faolain, A. Gomez-Iglesias, and T. F. Krauss, “Systematic design of flat band slow light in photonic crystal waveguides,” Opt. Express **16**, 6227–6232 (2008).
[CrossRef] [PubMed]

**5. **Y. A. Vlasov, M. O’Boyle, H. F. Hamann, and S. J. McNab, “Active control of slow light on a chip with photonic crystal waveguides,” Nature **438**, 65–69 (2005).
[CrossRef] [PubMed]

**6. **D. Mori and T. Baba, “Dispersion-controlled optical group delay device by chirped photonic crystal waveguides,” Appl. Phys. Lett. **85**, 1101–1103 (2004).
[CrossRef]

**7. **B. Corcoran, C. Monat, C. Grillet, D. J. Moss, B. J. Eggleton, T. P. White, L. O’Faolain, and T. F. Krauss, “Green light emission in silicon through slow-light enhanced third-harmonic generation in photonic-crystal waveguides,” Nature Photonics **3**, 206–210 (2009).
[CrossRef]

**8. **C. Kittel and P. B. Kahn, *Quantum theory of solids* (John Wiley and Sons, 1987).

**9. **J. M. Luttinger and W. Kohn, “Motion of electrons and holes in perturbed periodic fields,” Phys. Rev. **97**, 869–883 (1955).
[CrossRef]

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

**11. **G. Bastard, “Superlattice band structure in the envelope-function approximation,” Phys. Rev. B **24**, 5693–5697 (1981).
[CrossRef]

**12. **M. Notomi, K. Yamada, A. Shinya, J. Takahashi, C. Takahashi, and I. Yokohama, “Extremely large group-velocity dispersion of line-defect waveguides in photonic crystal slabs,” Phys. Rev. Lett. **87**, 253902 (2001).
[CrossRef] [PubMed]

**13. **K. B. Dossou, L. C. Botten, R. C. McPhedran, C. G. Poulton, A. A. Asatryan, and C. M. de Sterke, “Shallow defect states in two-dimensional photonic crystals,” Phys. Rev. A **77**, 063839 (2008).
[CrossRef]

**14. **J. D. Joannopoulos, S. G. Johnson, and J. N. Winn, *Photonic crystals: molding the flow of light* (Princeton university press, 2008).

**15. **A. Mekis, S. Fan, and J. D. Joannopoulos, “Bound states in photonic crystal waveguides and waveguide bends,” Phys. Rev. B **58**, 4809–4817 (1998).
[CrossRef]

**16. **M. Charbonneau-Lefort, E. Istrate, M. Allard, J. Poon, and E. H. Sargent. “Photonic crystal heterostructures: Waveguiding phenomena and methods of solution in an envelope function picture,” Phys. Rev. B **65**, 125318 (2002).
[CrossRef]

**17. **A. Yariv, Y. Xu, R. K. Lee, and A. Scherer, “Coupled-resonator optical waveguide: a proposal and analysis,” Opt. Lett. **24**, 711–713 (1999).
[CrossRef]

**18. **S. Mahmoodian, R. C. McPhedran, C. M. de Sterke, K. B. Dossou, C. G. Poulton, and L. C. Botten, “Single and coupled degenerate defect modes in two-dimensional photonic crystal band gaps,” Phys. Rev. A **79**, 013814 (2009).
[CrossRef]

**19. **L. D. Landau and E. M. Lifshitz, *Quantum mechanics: non-relativistic theory* (Pergamon Press, 1958).

**20. **C. M. de Sterke, L. C. Botten, A. A. Asatryan, T. P. White, and R. C. McPhedran “Modes of coupled photonic crystal waveguides,” Opt. Lett. **29**, 1384–1386 (2004).
[CrossRef]