## Abstract

We develop a systematic approach for simultaneous extraction of the dispersion relations and profiles of multiple modes in periodic waveguides though a special global optimization procedure applied to near-field electric field measurements in the waveguide plane. We apply this method to perform in-depth analysis of experimental data on wave propagation close to an interface between waveguide sections with different dispersion characteristics, and we successfully identify several modes contributing to the experimentally measured fields. We find clear evidence that when the group velocity is reduced across the interface, evanescent modes that facilitate the excitation of propagating slow-light waves appear, confirming previous theoretical predictions.

© 2011 Optical Society of America

## 1. INTRODUCTION

Periodically modulated optical waveguides offer new possibilities for controlling the propagation of light. Resonant scattering from periodic modulations can be used to tailor the dispersion, enabling, in particular, a dramatic modification of the group velocity and realization of slow-light propagation. Such fundamental effects can be directly visualized in experiments with near-field measurements, which can be used to recover the amplitude and phase of the electric field at all locations in the plane of the waveguide [1]. This information can then be used to determine the dispersion characteristics and spatial profiles of the guided modes.

A commonly used approach to dispersion extraction is through spatial Fourier transform (SFT) of the field profiles, because peaks in the Fourier spectra correspond to the wavenumbers of guided modes [2]. Such analysis is valuable for studying complex dispersion features such as turning, inflection, and anticrossing points [3, 4]. However, there exists a fundamental limitation on results obtained with SFT: $\mathrm{\Delta}k\ge 2\pi /L$, where $\mathrm{\Delta}k$ is the resolution of the wavenumber and *L* is the structure length. Therefore, accurate dispersion results can only be obtained for long waveguides, extending over many periods of the underlying photonic structure. Another limitation of the SFT method is that it is difficult to extract information on the dispersion of evanescent waves, which may play an important role close to the structure boundaries or interfaces between different waveguides. For example, evanescent waves enable efficient excitation of slow-light modes in photonic crystal waveguides without a transition region [5]. Moreover, all waves have decaying amplitudes in lossy media such as metal–dielectric metamaterial and plasmonic structures.

Alternative methods for dispersion extraction have been developed to overcome the shortcomings of the SFT method. An interference pattern of two counterpropagating modes in photonic crystals was used to extract their wavenumbers [6]. It was shown that in metamaterial structures, the effective refractive index can be determined through the extracted phase velocity of a single propagating or evanescent wave [7, 8]. However, these techniques are not applicable in the presence of multiple propagating or evanescent modes.

Recently, it was demonstrated that dispersion extraction in multimode waveguides with, in principle, unbounded resolution is possible even for short waveguide sections [9, 10], using approaches based on high-resolution spectral methods previously developed for the analysis of temporal dynamics [11, 12]. We have since introduced an important generalization of such methods by taking into account the spatial symmetry properties of modes in periodic waveguides [13]. It was shown that in addition to the dispersion relations, it is possible to extract the spatial profiles of all guided modes, including arbitrary combinations of propagating and evanescent waves. The method was demonstrated by application to numerical simulation data, where artificial noise was introduced to simulate experimental conditions.

In this paper, we further extend the method by formulating the global optimization procedure, which is most suitable for extracting the dispersion characteristics of multiple modes simultaneously for a range of frequencies. This enables more accurate determination of the group velocity and higher order dispersion, compared to the previously considered [13] local optimization for individual frequencies.

We apply the extraction method to experimental near-field measurement data, performed for a set of photonic crystal waveguides. The various photonic crystal waveguides were engineered to exhibit different dispersion properties in the slow-light regime. For each sample, we successfully extract the dispersion through a global optimization procedure, and we show that the values of mismatch between the measured and reconstructed fields are of the order of 10% or less. The extracted dispersion curves are also found to be in excellent agreement with direct numerical modeling of different sample designs. Additionally, we extract the dispersion and profiles of evanescent waves. We find that the amplitudes and decay rates of evanescent waves excited at the waveguide boundary increase when the propagating mode is slowed down, in agreement with the theoretical predictions [5].

The paper is organized as follows. We formulate the mathematical extraction procedure in Section 2. Then in Section 3 we apply the method to extract the dispersion and determine the individual mode profiles using experimentally measured near-field profiles in a set of photonic crystal waveguides with varying dispersion characteristics in the slow-light regime. Finally, we formulate conclusions in Section 4.

## 2. DISPERSION EXTRACTION METHOD

In this section, we formulate the method for the extraction of mode dispersion and the profiles based on near-field measurements of the electric field in periodic waveguides. In Subsection 2A, we review the translational symmetry properties of modes in periodic waveguides, which are relevant for our analysis. Then in Subsection 2B we review the extraction procedure based on local optimization for individual frequencies. Finally, in Subsection 2C we formulate the global optimization approach for the simultaneous dispersion extraction over a range of frequencies.

#### 2A. Symmetries of Modes in Periodic Waveguides

Let us consider a periodic waveguide section, where the light propagation in a particular frequency range is primarily determined by a finite total number of guided modes (*M*). The value of *M* can be established based on numerical modeling, taking into account both propagating and evanescent waves. Because each of the modes of a periodic waveguide satisfies the Bloch theorem [14], the complex electric field envelope of a waveguide mode with the index *m* at the frequency *ω* can be expressed as

*x*and

*y*are the orthogonal directions transverse to the waveguide and

*z*is the direction of periodicity;

*d*is the waveguide period; and ${\psi}_{m}$ are the periodic Bloch-wave envelope functions: ${\psi}_{m}(z)={\psi}_{m}(z+d)$. Then the total field inside the waveguide can be presented as a linear superposition of

*M*guided modes with amplitudes ${a}_{m}$:

#### 2B. Locally Optimal Extraction at Individual Frequencies

To begin the procedure for the simultaneous extraction of the wavenumbers and profiles of the guided modes at a particular frequency, we follow Ref. [13] and separate the spatial domain into a number of unit cells, ${z}_{\mathrm{min}}+(n-1)d\le z<{z}_{\mathrm{min}}+nd$, where $n=1\text{:}N$ and *N* is the number of periods in the waveguide section. Let us denote with

**r**belongs to the first unit cell.

If one considers this relation only for a single point **r** in the unit cell, it becomes mathematically equivalent to the problems considered in the spectral analysis of temporal series [11, 12], and high-resolution spectral methods can be used to extract the mode wavenumbers [10]. However, the special property of periodic waveguides is that Eq. (6) needs to be satisfied simultaneously for all spatial locations **r** in the unit cell. This allows us to determine the values of ${k}_{m}$ and ${A}_{m}(\mathbf{r};\omega )$, provided that the number of measurements exceeds the number of unknowns: $N\times {N}_{p}\ge M\times {N}_{p}+{M}_{k}$, where ${N}_{p}$ is the number of measurement points per unit cell and ${M}_{k}$ is the number of independent wavenumber values [13].

For each frequency *ω*, we seek the values of ${k}_{m}$ that provide the most accurate description of the whole measured field, and we apply the least-squares method to find a minimum of the functional

It follows that for each point **r** in a unit cell, the optimal amplitudes satisfy the linear matrix equation

*C*are ${C}_{np}=\mathrm{exp}[i{k}_{p}(n-1)]$, and vector $\tilde{U}(\mathbf{r};\omega )$ components are ${U}_{n}(\mathbf{r};\omega )$ for $p=1\text{:}M$ and $n=1\text{:}N$.

The minimization procedure is illustrated in Fig. 1. It starts with an initial guess for ${k}_{m}$, from which ${A}_{m}$ is obtained by using Eq. (9). Then ${\mathrm{w}}_{n}$ can be calculated from Eq. (6), and the process is repeated for every point **r** in a unit cell to determine *W* from Eq. (7). The entire process can be combined into a single expression as

#### 2C. Globally Optimal Extraction across a Frequency Range

The procedure described above in Subsection 2B enables one to extract the wavenumbers and amplitudes of guided modes for each *ω*. However, in experimental data, the noise level can vary widely between different frequencies, especially if the group velocity has a strong frequency dependence, as is the case in the examples presented in Section 3 below. This can dramatically reduce the accuracy for determining group velocity and higher order dispersion characteristics, which are defined through the derivatives of the dispersion curves. In order to overcome this issue, we suggest here a globally optimal extraction procedure, which recovers dispersion properties simultaneously across a range of frequencies.

To perform the global extraction, we first identify the expected general shape of the dispersion curves. Suppose that for the considered frequency range, dispersion of mode number *m* can be approximated by Taylor expansion, i.e.,

*ω*(a particular branch will need to be chosen based on physical considerations). The significance of this transformation is that the modal dispersion ${k}_{m}(\omega ,\{Q\})$ can be fully described by the

*same*set of parameters $\{Q\}$ for

*all values*of

*ω*in the data. Thus, instead of extracting ${k}_{m}$ for each

*ω*one at a time, we can perform a single global optimization and extract ${k}_{m}$ for all

*ω*:

*ω*.

## 3. DISPERSION AND MODE PROFILE EXTRACTION FROM NEAR-FIELD EXPERIMENTAL MEASUREMENTS

In this section we apply the Bloch mode extraction techniques described in the previous section to experimental near-field optical measurements of silicon photonic crystal waveguides.

#### 3A. Dispersion-Engineered Photonic Crystal Waveguides

In order to demonstrate the ability of the general methods introduced in Section 2 above to extract the dispersion properties and mode profiles of multiple propagating and evanescent modes, we have specially designed the photonic crystal waveguides by locally modifying the lattice geometry on either side of the waveguide [15, 16, 17]. Specifically, following the previous analysis [5], we have engineered waveguides supporting slow-light modes away from the edge of the Brillouin zone, because in such structures both the propagating and evanescent modes can be simultaneously excited at the waveguide boundary [5].

We present results for three different geometries where the waveguide dispersion was engineered by shifting two rows of holes on each side of the waveguide parallel to the waveguide axis, as illustrated in Figs. 2d, 3d, 4d. The first rows were moved by ${p}_{1}=0.3\text{\hspace{0.17em}}d$, and second rows by ${p}_{2}=0.375\text{\hspace{0.17em}}d$, $0.4\text{\hspace{0.17em}}d$ and $0.425\text{\hspace{0.17em}}d$, respectively, in the same direction. These values were chosen based on 3D plane-wave dispersion calculations using the MIT photonic bands software package [18] in order to obtain progressively more “extreme” dispersion features in the slow-light region; see the numerically calculated dispersion curves shown with dotted lines labeled “*p*” in Figs. 2a, 3a, 4a, where the slow-light regime is realized at wavenumbers around $kd/2\pi \simeq \pm 0.6$. The first two waveguides (Figs. 2, 3) support a single propagating TE-polarized mode at all wavelengths of interest, and a narrow slow-light region corresponding to an inflection point in the dispersion curve. The third waveguide (Fig. 4) has a more complex dispersion curve: the slow-light region has been further distorted to introduce two turning points where the group velocity is zero and a region between them that supports three propagating modes.

The dispersion curves in Figs. 2a, 3a are similar to those studied in Ref. [5], where it was shown that very efficient excitation of slow modes at dispersion inflection points is possible without a transition region. Such efficient excitation is enabled by weakly evanescent modes that assist in mode matching at the interface, but do not transmit any power. To provide the required “fast” mode from which to couple to the slow mode, the fabricated structures were designed to include on either end of the waveguide a short 10-period section where the lattice period is increased by $30\text{\hspace{0.17em}}\mathrm{nm}$ (to $450\text{\hspace{0.17em}}\mathrm{nm}$) in the direction of the waveguide. This waveguide section supports a mode with a group-velocity of approximately $c/5$ in the studied wavelength range. The complex shape of the dispersion curves and the presence of the evanescent modes make it an ideal case to test the advanced analysis techniques presented in this paper.

The photonic crystal structures were fabricated using electron-beam lithography and dry etching, as described in Ref. [15]. The samples consisted of $80\text{\hspace{0.17em}}\mathrm{\mu m}$ long silicon membrane photonic crystal waveguides with a lattice period $d=420\text{\hspace{0.17em}}\mathrm{nm}$ and a hole radius $r=128\text{\hspace{0.17em}}\mathrm{nm}$ $(0.305\text{\hspace{0.17em}}d)$, connected to silicon ridge access waveguides at either end. To study the evanescent modes directly, we perform phase- and polarization-sensitive near-field optical measurements on the waveguide [19, 20]. The near-field electric field profiles were measured along the center of the waveguide core.

#### 3B. Spatial Fourier Transform Spectra

Before carrying out the extraction using the methods presented in Section 2, the SFT spectra of the near-field electric field profiles were first investigated because they give useful information about the modes. However, the resolution of the SFT is poor, as it is limited by an inverse of the waveguide length; see Figs. 2c, 3c, 4c. For all experimental samples, we can see the modes appearing at $kd/2\pi \simeq \pm 0.6$ (and their Bloch harmonics at $kd/2\pi \simeq \pm 0.4$), as expected, based on numerical simulations shown with the dotted dispersion curves in Figs. 2a, 3a, 4a, respectively, although the slow-light regime is hardly distinguishable in the SFT spectra.

Also visible in the SFT spectra is an additional mode at $kd/2\pi \simeq 0.345$, which was observed in all three waveguides [Figs. 2c, 3c, 4c]. It is consistent with the fundamental TM-like mode of the waveguide [the numerical dis persion of which is also shown with dotted lines labeled “TM” in Figs. 2a, 3a, 4a, respectively], which is tightly confined to the waveguide core and thus is almost unaffected by the ${p}_{1}$ and ${p}_{2}$ shifts. Its presence is most likely due to scattering at the access to the silicon ridge waveguide.

#### 3C. Application of Optimal Extraction Method

Based on the SFT spectra and numerical modeling of the dispersion curves, we see that the structure supports two types of propagating modes: a TE-like mode, which exhibits slow-light propagation for certain frequencies, and a TM-like mode with simple linear dispersion. Following the approach of global optimization formulated above in Subsection 2C, we now represent the mode dispersion through Taylor series expansions with coefficients that need to be determined.

The dispersion of the TE mode can feature an inflection point, and such dependencies can be described by a third- degree polynomial, i.e., we put ${j}_{1}=3$ for the first mode, the dispersion of which is then approximated as

As we invert this relation, we obtain the dependence ${k}_{1}(\omega )$. Because this is a third-degree polynomial, for each*ω*there are three solutions for ${k}_{1}$, ${k}_{2}$, and ${k}_{3}$. Because the dependence $\omega (k)$ is an analytic function in lossless dielectric photonic structures, all three roots correspond to actual modes that can be excited in the structure. In particular, in the vicinity of an inflection point in the dispersion dependence, there will be one real ${k}_{1}$ corresponding to a propagating mode and a pair of complex wavenumbers ${k}_{2}={k}_{3}^{*}$ corresponding to the evanescent modes. The presence of evanescent modes is actually of significant physical importance, as they can facilitate efficient light coupling to the slow-light mode [5, 21], and the amplitude of evanescent waves can be very large at the waveguide boundary. It is therefore critically important to include the evanescent modes in the extraction procedure. Accordingly, we include in the analysis all three roots ${k}_{1,2,3}(\omega )$ of the cubic polynomial Eq. (13). As we are considering a dielectric structure, the dispersion has a symmetry between forward and backward modes, $k\to -k$. Accordingly, we also need to take into account the modes with ${k}_{4,5,6}(\omega )=-{k}_{1,2,3}(\omega )$.

We also need to include in the analysis the TM-like modes, which dispersion is a simple linear dependence (${j}_{7}=1$):

such that ${k}_{7}=(\omega -{Q}_{7}^{(0)})/{Q}_{7}^{(1)}$. There also can be a backward propagating mode with ${k}_{8}(\omega )=-{k}_{7}(\omega )$.To summarize, we consider eight modes: six TE (three forward and three backward, including propagating and evanescent modes) and two propagating (forward and backward) TM modes. The dispersion of all these modes is defined through six independent coefficients: ${Q}_{1}^{(0)}$, ${Q}_{1}^{(1)}$, ${Q}_{1}^{(2)}$, ${Q}_{1}^{(3)}$, ${Q}_{7}^{(0)}$, and ${Q}_{7}^{(1)}$. Their values are determined through the optimization procedure, as described in Subsection 2C, where we estimate the initial guesses for $\{Q\}$ based on the SFTs. As a result, we extract the dispersion for all the dominant propagating and evanescent waveguide modes.

#### 3D. Results of Dispersion Extraction

The extraction was carried out for the slow waveguide section of $\sim 50$ periods long where the evanescent modes were strong, i.e., near the input interface between the fast waveguide and slow waveguide. Because the exponentially growing evanescent modes in the propagating direction (*z*) should be very weak [5], they were excluded from the extraction.

The extracted wavenumbers of the propagating and evanescent waves (labeled “*p*” and “*e*,” respectively) for different photonic structures are shown in Figs. 2a, 3a, 4a, and the corresponding mismatch values ${W}_{\mathrm{min}}$ are shown in Figs. 2b, 3b, 4b. We see that, in all cases, the shape of the dispersion curves matches very closely the numerically designed types. For the structures with ${p}_{2}=0.375\text{\hspace{0.17em}}d$ and ${p}_{2}=0.4\text{\hspace{0.17em}}d$ shown in Figs. 2, 3 respectively, the mismatch values for wavelengths in the slow-light region are below 15%. The increase of the mismatch in the slow-light region is possibly due to the increased sensitivity of the disorder of the slow-light modes [22], which would break the structure of the eigenmodes, as they would no longer satisfy the Bloch periodicity condition, which is assumed in the extraction procedure.

For the structure with ${p}_{2}=0.425\text{\hspace{0.17em}}d$ shown in Fig. 4, the dispersion curve has two turning points instead of an inflection point, and the extracted wavenumbers in Fig. 4 clearly demonstrate this feature. In the wavelength region between the turning points, the complex wavenumbers of the evanescent modes become purely real and they all turn into propagating modes. Despite this difference in the dispersion characteristics, the modes were successfully extracted in the same way as before, although the mismatch values in this region are larger (up to 30%). We expect that this increase in the mismatch can be attributed to even larger sensitivity to disorder of the slow-light modes in regions close to turning points in dispersion curves where $\partial \omega /\partial k=0$ and ${\partial}^{2}\omega /\partial {k}^{2}\ne 0$ compared to slow-light near inflection points [22]. Conversely, the large mismatch values could serve as an indication of the increased disorder effects.

In addition to dispersion extraction, our procedure also recovers the spatial profiles of all the individual modes. To illustrate the effect of the evanescent modes, we present here the mode profiles for the case of ${p}_{2}=0.4\text{\hspace{0.17em}}d$ (corresponding to Fig. 3) at $\lambda =1537\text{\hspace{0.17em}}\mathrm{nm}$ and $1534\text{\hspace{0.17em}}\mathrm{nm}$ in Figs. 5, 6, respectively. In both cases, the sum of the mode profiles is in good agreement with the measured field, see plots (c) and (d) in these figures, which confirms that the experimental profile can indeed be decomposed into a sum of Bloch waves. For the wavelength of $\lambda =1537\text{\hspace{0.17em}}\mathrm{nm}$, which is detuned away from the slow-light region, the evanescent modes decay quickly away from the waveguide boundary; see Fig. 5b. On the other hand, the evanescent modes penetrate deeper into the structure and play an essential role in the total field distribution when the wavelength is tuned to the slow-light region at $1534\text{\hspace{0.17em}}\mathrm{nm}$; see Fig. 6b. These differences are consistent with the previous theoretical predictions on the role of evanescent modes facilitating light coupling to the slow-light waveguides [5, 13, 21].

In order to further confirm the presence of evanescent waves, we have repeated the dispersion extraction calculations by only taking into account the propagating waves. The corresponding results presented in Fig. 7 show that the mismatch becomes very high in the slow-light region, increasing by over three times compared to analysis including evanescent waves presented in Fig. 3 above. This provides additional confirmation of the excitation of evanescent waves and their key role in light reshaping at the interfaces of slow-light waveguide sections.

#### 3E. Comparison with Local Extraction Results

Finally, we perform a comparison of the extraction results based on the global optimization result presented in Subsection 3D above, with the results obtained through a local optimization procedure introduced in Subsection 2B. For comparison, we show in Fig. 8 the results of the locally optimal dispersion extraction using the same parameters as in Fig. 3, where a globally optimal extraction procedure was used. By comparing these figures, we see that the local optimization provides only very small improvement to the mismatch values in the slow-light region, to the detriment of the highly irregular extracted dispersion shape. Also, the locally optimized mismatch values even turn out to be larger than for the global optimization procedure at certain frequencies, and the average mismatch appears to be smaller in the case of global optimization. This happened because it is numerically challenging to locate the minimal mismatch through local multidimensional optimization, whereas in global optimization knowledge of the generic dispersion properties leads to extraction results that more accurately recover the mode characteristics.

## 4. CONCLUSIONS

In conclusion, we have presented the extraction of wave numbers and amplitude profiles of multiple propagating and evanescent Bloch modes in periodic waveguides. Most importantly, we were able to recover dispersion curves from experimental near-field data by using the knowledge of the overall shapes of the dispersion dependencies and performing a single global extraction for all frequency points. This approach enables us to correctly extract and identify different types of slow-light dispersion in the vicinity of inflection points and confirm the presence of evanescent modes close to the structure boundary.

## ACKNOWLEDGMENTS

This work was supported by the Australian Research Council.

**1. **R. J. P. Engelen, Y. Sugimoto, H. Gersen, N. Ikeda, K. Asakawa, and L. Kuipers, “Ultrafast evolution of photonic eigenstates in *k*-space,” Nat. Phys. **3**, 401–405 (2007). [CrossRef]

**2. **H. Gersen, T. J. Karle, R. J. P. Engelen, W. Bogaerts, J. P. Korterik, N. F. van Hulst, T. F. Krauss, and L. Kuipers, “Direct observation of Bloch harmonics and negative phase velocity in photonic crystal waveguides,” Phys. Rev. Lett. **94**, 123901 (2005). [CrossRef] [PubMed]

**3. **H. Gersen, T. J. Karle, R. J. P. Engelen, W. Bogaerts, J. P. Korterik, N. F. van Hulst, T. F. Krauss, and L. Kuipers, “Real-space observation of ultraslow light in photonic crystal waveguides,” Phys. Rev. Lett. **94**, 073903 (2005). [CrossRef] [PubMed]

**4. **N. Le Thomas, V. Zabelin, R. Houdre, M. V. Kotlyar, and T. F. Krauss, “Influence of residual disorder on the anticrossing of Bloch modes probed in *k* space,” Phys. Rev. B **78**, 125301 (2008). [CrossRef]

**5. **T. P. White, L. C. Botten, C. M. de Sterke, K. B. Dossou, and R. C. McPhedran, “Efficient slow-light coupling in a photonic crystal waveguide without transition region,” Opt. Lett. **33**, 2644–2646 (2008). [CrossRef] [PubMed]

**6. **S. H. Fan, I. Appelbaum, and J. D. Joannopoulos, “Near-field scanning optical microscopy as a simultaneous probe of fields and band structure of photonic crystals: a computational study,” Appl. Phys. Lett. **75**, 3461–3463 (1999). [CrossRef]

**7. **B. I. Popa and S. A. Cummer, “Determining the effective electromagnetic properties of negative-refractive-index metamaterials from internal fields,” Phys. Rev. B **72**, 165102 (2005). [CrossRef]

**8. **A. Andryieuski, R. Malureanu, and A. V. Lavrinenko, “Wave propagation retrieval method for metamaterials: unambiguous restoration of effective parameters,” Phys. Rev. B **80**, 193101 (2009). [CrossRef]

**9. **B. Dastmalchi, A. Mohtashami, K. Hingerl, and J. Zarbakhsh, “Method of calculating local dispersion in arbitrary photonic crystal waveguides,” Opt. Lett. **32**, 2915–2917 (2007). [CrossRef] [PubMed]

**10. **A. A. Sukhorukov, S. Ha, I. V. Shadrivov, D. A. Powell, and Yu. S. Kivshar, “Dispersion extraction with near-field measurements in periodic waveguides,” Opt. Express **17**, 3716–3721 (2009). [CrossRef] [PubMed]

**11. **R. Roy, B. G. Sumpter, G. A. Pfeffer, S. K. Gray, and D. W. Noid, “Novel methods for spectral analysis,” Phys. Rep. **205**, 109–152 (1991). [CrossRef]

**12. **V. A. Mandelshtam, “FDM: the filter diagonalization method for data processing in NMR experiments,” Prog. Nucl. Magn. Reson. Spectrosc. **38**, 159–196 (2001). [CrossRef]

**13. **S. Ha, A. A. Sukhorukov, K. B. Dossou, L. C. Botten, C. M. de Sterke, and Yu. S. Kivshar, “Bloch-mode extraction from near-field data in periodic waveguides,” Opt. Lett. **34**, 3776–3778 (2009). [CrossRef] [PubMed]

**14. **J. D. Joannopoulos, R. D. Meade, and J. N. Winn, *Photonic Crystals: Molding the Flow of Light* (Princeton University Press, 1995).

**15. **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]

**16. **S. Kubo, D. Mori, and T. Baba, “Low-group-velocity and low-dispersion slow light in photonic crystal waveguides,” Opt. Lett. **32**, 2981–2983 (2007). [CrossRef] [PubMed]

**17. **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]

**18. **S. G. Johnson and J. D. Joannopoulos, “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Express **8**, 173–190 (2001). [CrossRef] [PubMed]

**19. **M. L. M. Balistreri, J. P. Korterik, L. Kuipers, and N. F. van Hulst, “Local observations of phase singularities in optical fields in waveguide structures,” Phys. Rev. Lett. **85**, 294–297 (2000). [CrossRef] [PubMed]

**20. **M. Burresi, R. J. P. Engelen, A. Opheij, D. Oosten, D. van Mori, T. Baba, and L. Kuipers, “Observation of polarization singularities at the nanoscale,” Phys. Rev. Lett. **102**, 033902 (2009). [CrossRef] [PubMed]

**21. **C. M. de Sterke, K. B. Dossou, T. P. White, L. C. Botten, and R. C. McPhedran, “Efficient coupling into slow light photonic crystal waveguide without transition region: role of evanescent modes,” Opt. Express **17**, 17338–17343 (2009). [CrossRef]

**22. **M. Patterson, S. Hughes, S. Schulz, D. M. Beggs, T. P. White, L. O’Faolain, and T. F. Krauss, “Disorder-induced incoherent scattering losses in photonic crystal waveguides: Bloch mode reshaping, multiple scattering, and breakdown of the Beer– Lambert law,” Phys. Rev. B **80**, 195305 (2009). [CrossRef]