{'fraction_non_alphanumeric': 0.031995490790935166, 'fraction_numerical': 0.015334626415344573, 'mean_word_length': 4.7106882514437185, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electromagnetic modes of a disordered photonic crystal

We present a systematic numerical approach to compute the eigenmodes and the related eigenfrequencies of a disordered photonic crystal, characterized by small fluctuations of the otherwise periodic dielectric profile. The field eigenmodes are expanded on the basis of Bloch modes of the corresponding periodic structure, and the resulting eigenvalue problem is diagonalized on a truncated basis including a finite number of Bloch bands. The Bloch-mode expansion is very effective for modeling modes of very extended disordered structures in a given frequency range, as only spectrally close bands must be included in the basis set. The convergence can be easily verified by repeating the diagonalization for an increased band set. As illustrations, we apply the method to the W1 line-defect waveguide and to the L3 coupled-cavity waveguide, both based on a photonic crystal slab with a triangular lattice of circular holes. We compute and characterize the eigenfrequencies, spatial field profiles and radiation loss rates of the localized modes induced by disorder. For the W1 waveguide, we demonstrate that radiation losses, at the bottom of the spatially-even guided band, are determined by a small hybridization with Bloch modes of the spatially-odd band, induced by disorder in spite of their frequency separation. The Bloch-mode expansion method has a very broad range of applications: it can be also used to accurately compute the modes of structures with systematic modulations of the periodic dielectric constant, as in several designs of high-Q cavities.


I. INTRODUCTION
Thanks to the periodic spatial dependence of their dielectric constant, photonic crystal (PHC) structures 1-3 make it possible to engineer the confinement and propagation of the electromagnetic field, and hold great promise for applications to integrated optics.
Similarly to atomic crystals, PHCs always present some amount of disorder as a result of the fabrication process. 4 When considering PHC slabs -a periodic array of cylindric holes in a dielectric slab waveguide -disorder arises from the irregular shape of the holes and from the roughness of the dielectric interfaces. Disorder is mostly considered as a disturbance, in the context of applications. It prevents ballistic light propagation, by inducing losses due to scattering, [5][6][7][8][9][10][11][12][13][14][15][16] and possibly strong localization of light, 5,16-24 thus constituting a major physical limit to the use of PHCs as interconnects in photonic circuits, and in all applications relying on slow light propagation. [25][26][27] Disorder-induced scattering is also responsible for extrinsic radiation losses, 6,28,29 as it mixes Bloch momenta of truly guided and leaky modes of the nominally periodic structure.
More generally, it has been suggested that PHCs might be an election system for studying the fundamental properties of light scattering in disordered media. It is rather singular that the seminal work by Sayeed John, 2 considered as one of the first works to introduce the idea of a PHC, actually considers the photonic bands only as a tool to achieve Anderson localization of light. 30 The idea put forward in the work by John is that, in the vicinity of a point of high symmetry, the dispersion of a photonic band can be approximated by a parabola and Maxwell equations can consequently be mapped onto an effective Schrödinger equation in a continuous effective medium, including a random potential term derived from the actual disorder profile of the system.
Modeling the effects of disorder in PHCs is a current research topic and a very challenging computational task. In a regular structure, Bloch's theorem allows restricting the solution of Maxwell equations to a single elementary cell. Taking advantage of the crystal periodicity, several computational approaches have been developed, based on mode expansion 28,31,32 or scattering matrix [33][34][35][36] formalisms. Disorder however breaks the translational symmetry of the crystal, and Bloch's theorem can no longer be applied. Then, a full solution of Maxwell equations over an extended region of space, possibly including thousands of elementary cells would be needed both for studying fundamental properties, such as Anderson localization, and for modeling device operation. First principle calculations based on a finite-element sampling of real space are thus of very limited use, as the computational cost of simulating such extended samples would be prohibitive. Theoretical approaches to disordered PHCs have been proposed, based on various levels of approximation. Some properties, such as diffusive light propagation and the corresponding radiation losses, can be predicted already within a perturbation theory of light scattering. 6,16,29,37 More generally, a variety of coupledmode approaches have been proposed for studying disordered PHCs. 8,12,14,38 A coupled-mode theory in general accounts for light scattering between few selected modes. The reduced number of modes makes it possible to model the scattering processes in a fully nonperturbative fashion, but restricts the predictive power of the method. As an example, Anderson localization of light arises from multiple coherent scattering between many different propagating modes, 22,30,39 and typically requires a solution of Maxwell equations beyond fewmode coupling, to be modeled correctly. It should be noted however that the idea of mode coupling can be generalized to scattering theories including a full basis of propagating modes. In this case, the theory is formally equivalent to a full solution of Maxwell equations, as is e.g. the case for the quasi-normal Bloch-mode scattering matrix approach recently developed by Lecamp et al. 9,40 A general aspect of most theoretical approaches [33][34][35][36]40 is the fact that they compute directly the linear response of the periodic dielectric medium. For several reasons however, it would be beneficial to compute the actual Maxwell eigenmodes. Apart from the direct visualization of the Anderson-localized spatial profiles, knowledge of the actual eigenmodes and of their detailed Bloch-mode components can give insight into the mechanisms of slowlight propagation, 25,27,41-44 radiation 6,28,29 and disorderinduced [5][6][7][8][9][10][11][12][13][14][15][16] losses. When modeling systems with coupled photonic and electronic degrees of freedom, such as quantum dots embedded in cavities or guides, 19,[45][46][47][48][49][50][51][52][53] or PHC polaritons, 31,54 it is most natural to start with the eigenstates of both coupled subsystems, especially within a fully quantum mechanical treatment where second quantization of electromagnetic modes is needed. If the linear response function is known, then eigenmodes can be obtained by finding the poles of its analytical continuation on the complex frequency plane. 35 Conversely, if Maxwell eigenmodes are known, the linear response function can always be expressed as a function of these eigenmodes, by using the resolvent representation in Fredholm theory. 55 Both indirect ways are however computationally far from being optimal. A method for directly computing Maxwell eigenmodes of an extended PHC with a small perturbation of the periodic dielectric profile would thus be an indispensable complement to linear response theories.
Here, we develop the Bloch-mode expansion formalism, for directly computing Maxwell eigenmodes of a photonic structure characterized by a small perturbation of the otherwise periodic dielectric profile. The basic idea behind this formalism is that eigenmodes of the perturbed structure can be expanded on the basis of Bloch modes of the regular, unperturbed PHC. This expansion extends, in general, to all bands and all Bloch momenta of the regular structure. If the perturbation is small, as in the case of disorder in state-of-the-art PHC slabs, then the expansion can be truncated on a relatively small subset of Bloch bands -spectrally close to the frequency region of interest -and the Maxwell equation can be diagonalized restricting to this subspace, at a considerably reduced computational cost. The convergence of the method can be easily assessed by repeating the calculation on a progressively augmented Bloch-band subset. After laying down the general formalism, as examples of application we study the W1 line-defect waveguide 56,57 and the L3 coupled-cavity waveguide, 58 both based on a PHC slab structure with a triangular lattice of circular holes. For computing the Bloch modes of the regular PHC, we use the guided-mode expansion method recently proposed by Andreani and Gerace, 28 because it offers a good balance between predictivity and ease of implementation, although some limitations in its accuracy to predict outof-plane losses have been recently put forth. 9 This is not however the only possible choice, and the Bloch-mode expansion could equally rely on a first-principle numerical solution of Maxwell equations for computing Bloch modes of the regular structure. We show how a fully converged calculation of modes in the guided region of both waveguides can be achieved at a minimal computational cost, for structures as long as 1024 elementary PHC cells. We study the localization of the modes and their radiation loss rates, and discuss the influence of disorder and losses on slow-light propagation, in the light of the existing literature. For the W1 waveguide, we show that radiation losses at the bottom edge of the spatiallyeven guided band are determined by hybridization with the spatially-odd guided band, in spite of their mutual frequency separation. This result explains the physical limitations to the Q-factor of W1-based structures, such as "gentle confinement" cavities, [59][60][61] and shows the way to a significant optimization of these systems.
In Section II, we present the theoretical formalism. We also discuss the method for computing radiation loss rates, when Bloch modes are computed via guided-mode expansion. In Sections III and IV, we discuss the results of the application of the present method to a W1 linedefect waveguide and to a L3 coupled-cavity waveguide, respectively. Section V presents the main conclusions and an outlook of future applications.

II. THEORY
We consider a PHC made of a periodic array of cylindrical air holes etched in a dielectric slab of thickness d. The regular PHC -namely in the absence of disorderis described by the dielectric constant ǫ(r). In presence of disorder instead, the dielectric constant is defined as As the holes are present only in the slab, we assume that both ǫ(r) and ǫ ′ (r) take constant values in the upper and lower semi-infinite claddings. In the following, we will always consider an "air bridge" configuration, for which ǫ = ǫ ′ = 1 and δǫ = 0 in the claddings. Here, we are neglecting the roughness that might characterize the slab interfaces, as we are primarily interested in the disorder affecting the cylindrical holes. The starting point of our analysis are the Bloch modes E kn (r) of the regular PHC structure under study. They are the eigenmodes of the Maxwell equation where ω kn are the corresponding eigenfrequencies. These modes might have been obtained by any computational method that takes advantage of Bloch's theorem. For example, they could be computed by exact finite-element solution of Maxwell equations, or by an approximate scheme like e.g. guided-mode expansion. 28 In the examples of application that follow, we will consider a finite system size with periodic boundary conditions, thus implying a discrete set of Bloch-momentum values k. The formalism presented in this Section can however be straightforwardly generalized to a continuous kspectrum.
In presence of disorder, the actual eigenmodes of the system are denoted by E β (r), and are indexed by a global index β, as the Bloch momentum k is no longer a conserved quantity. These modes are solution of the Maxwell equation They can be formally expanded on the basis of the Bloch modes of the regular structure: By replacing (4) into (3), and using (1) and (2), we obtain We take the scalar product of Eq. (5) by E * k ′ n ′ (r) and integrate on r (we assume the modes E kn (r) to be normalized over the total volume of the system). By taking advantage of the orthogonality relation 62 we finally obtain Eq. (7) is a generalized eigenvalue problem, formally equivalent to the full solution of Maxwell equation (3). The matrix elements V kn,k ′ n ′ contain all information about the disorder profile. They are defined as It should be pointed out that, when assuming periodic boundary conditions for the Bloch modes, the Blochmode expansion (4) implies periodic boundary conditions also for the modes E β (r). Thus, boundary effects at the side terminations of the PHC (e.g. interface with a ridgewaveguide), that might give rise to additional spectral features in some structures, 10 should be computed additionally, e.g. through mode-coupling or transfer-matrix approaches that go beyond the scope of the present work.
It is interesting to recall the principle underlying the present formalism had partly been laid down already in the seminal work by S. John. 2 There, it is pointed out that, close to a band edge, the Maxwell problem can be approximated by an effective Schrödinger equation whose kinetic term is essentially determined by the curvature of the single Bloch band under study. The potential term is nonlocal and is given by the matrix element of the random potential between actual electric field modes. This exactly corresponds to our Eqs. (7) and (8) if restricted to the edge of a single Bloch band. According to the basic theory of electronic states in crystals, restricting to the edge of a single Bloch band coincides with the effectivemass approximation. The present formalism generalizes John's original idea to an expansion over several bands. Its advantage lies in the possibility to truncate the vector space defining the eigenvalue problem (7) to a subset of the bands of the regular PHC, in order to ensure optimal convergence of the properties under study. In particular, when studying the effect of disorder on the band n 0 , it should be possible to neglect all bands n whose frequencies are much further away from that of band n 0 , than the disorder contribution V . By inspection of Eq. (7), it is easily shown that bands can be neglected if the condition is fulfilled. The reliability of the method can be assessed by repeating the diagonalization of (7) for an increasing number of bands and checking for convergence of the eigenvalues and eigenvectors of interest. This perturbation criterion holds however only for n = n 0 . Within the same band -in particular for low-dimensional systems -even the slightest amount of disorder can induce a strong mixing of Bloch modes in correspondence to band extrema, possibly giving rise to Anderson localization, as originally suggested by John. 2 The typical situation in which condition (9) is particularly favorable is that of defect states in the bandgap of a twodimensional PHC, for which most other bands are naturally separated by the frequency bandgap. Among the systems falling in this class are all one-dimensional systems based on a two-dimensional PHC, such as linedefect waveguides, 56,57 coupled-cavity waveguides, 63,64 etc. These systems are of particular interest because of the potentially high quality factor and of the possibility of slow-light propagation. 25,27 Slow-light in particular can be exploited to enhance the coupling of light to the electronic states of embedded quantum dots, for applications as single photon emitters. 19,[45][46][47][48][49][50][51][52][53] In this context in particular, knowing the actual eigenmodes of the electromagnetic field is crucial for modeling their coupling to electronic states. 55 Otherwise, the method should also be useful to model the states of a disordered two-dimensional PHC at the edges of the bandgap. 17,20 Finally, it is also possible to apply the present method to special designs involving a small systematic perturbation δǫ(r) of a periodic PHC, as in the case of long defect cavities 60 or cavities obtained by a slow grading of the lattice parameters. [59][60][61] In the next section we will apply the present method to a W1 line-defect waveguide in a PHC made of a triangular lattice of air holes, and to the related system of a L3 coupled-cavity waveguide. In the case of the W1 waveguide, we will show that accounting for more than one Bloch band is essential in determining the properties of the Maxwell eigenmodes. One of the main consequences of disorder in PHCs are the extrinsic radiation losses. In a regular structure, modes E kn (r) lying outside the light cone in the (ω k , k)-space are truly guided within the slab and ideally lossless. Inside the light cone, on the other hand, the electromagnetic modes are characterized by a continuous frequency spectrum, with a spectral function that has narrow resonances in correspondence to the band dispersion ω kn . In this case, formally the index n only labels the resonance frequencies, while eigenmodes should be labelled by a continuous frequency variable. The continuous frequency spectrum is related to the lossy nature of these modes, and implies a finite rate of radiation out of the plane into the claddings. When disorder is present, Eq. (4) mixes discrete and continuous parts of the spectrum, resulting in a continuous spectrum over the whole (ω k , k)-space. Now all modes (4) have components in the continuous part of the spectrum, inside the light cone, and are therefore lossy. For these lossy modes, Eqs. (4)- (8) would not strictly hold and should be modified to account for the continuous frequency spectrum. Then, a formally exact but computationally demanding approach to evaluate the loss rates, would require solving the Maxwell eigenvalue problem (2) numerically in three dimensions, and investigating the properties of its continuous spectrum. Here we prefer to adopt a computationally less demanding approach, based on the guided mode expansion recently proposed by Andreani et al. 28 to compute the modes of the regular structure E kn (r). These modes are expanded on a subset of the full Maxwell basis, formed by the guided modes of an effective homogeneous dielectric slab where α runs over the different guided modes. The leaky modes of the effective slab are thus neglected in this expansion. This approximation turns out to be highly predictive in typical situations, at a marginal computational cost. Radiation losses of Bloch modes can then be evaluated by a perturbative approach analogous to the Fermi golde rule in quantum mechanics, involving overlap integrals between modes (10) and leaky modes of the effective slab. The continuous spectrum of leaky Bloch modes is thus accounted for by an irreversible decay of the approximately guided modes (10) into the leaky modes of the homogeneous slab. The same approach can be applied in the present context to compute the radiation loss rates of the eigenmodes E β (r) of the disordered structure. Then, to the real eigenvalues ω 2 β /c 2 is associated an imaginary part Γ β related to the loss rate. The expression for this quantity is where g = k + G, ρ j is the density of leaky modes for a given frequency and in-plane momentum, and H rad βg is the overlap integral between mode E β (r) and the leaky mode of the effective slab with momentum g. The index j runs over the two cladding layers, while the index λ runs over the polarization basis of the radiative modes. From here, the actual loss rate can be computed as The explicit expressions for the quantities entering Eq. (11) are presented in the Appendix. One limit of validity of this perturbative approach is of course that γ β ≪ ω β . In presence of disorder, an additional condition must be fulfilled. While in a regular structure light would propagate ballistically in each Bloch mode E kn , disorder is generally responsible for damping of the transmission along the plane, induced by multiple scattering between different Bloch momenta and in particular coherent backscattering. 9,12 This damping is typically described as a disorder-induced loss mechanism. The same multiple scattering can induce Anderson localization of modes E β (r), and the average attenuation length of the transmission at a given frequency is a measure of the average localization length of modes at that frequency. In this picture, a competition between radiation losses and multiple scattering arises. The destructive interference in multiple scattering, at the origin of forward attenuation and Anderson localization, requires that radiation losses are negligible over the total multiple scattering path. As we plan to compute modes E β (r) based on a guided mode expansion of Bloch modes E kn (r) -thus without accounting for radiation losses -then an additional condition of validity of the present approach is that the attenuation length corresponding to γ β be much longer than the localization length of modes E β (r). If this is not the case, then Eq. (4) would be inconsistent with the radiation loss rate γ β computed perturbatively afterwards, as the Bloch waves involved in the expansion would be damped before producing the interference pattern that gives rise to the localized mode. This condition holds particularly well for modes that in the regular structure would lie outside the light cone and thus be radiation lossless. A typical example are the guided modes of a W1 line-defect waveguide, as discussed in the next Section.
Before proceeding to applications, we point out that the present method is fully non-perturbative. If the Bloch-mode basis, on which the expansion is made, is computed exactly, then the converged expansion (4) coincides with the full solution of the Maxwell problem for the extended disordered structure. If, as we do in the examples below, the Bloch basis is computed via guided-mode expansion, 28 then modes (4) are still close to the full solution of Maxwell equations, with the only assumption that the admixture of out-of-plane propagating modes is negligible and can be taken into account perturbatively at a later stage through Eq. (11).

III. APPLICATION TO A W1 LINE-DEFECT WAVEGUIDE
As an example of application of the present method, we consider the problem of disorder in W1 line-defect waveguides. We choose this system, as it was widely investigated both theoretically and experimentally, thus representing a valid platform on which our method can be tested. The effect of disorder on one-dimensional systems is particularly dramatic, as all states are expected to be exponentially localized. 39 In the context of PHC waveguides, the role of disorder was highlighted already in early experimental studies. Two disorder-related effects are of primary importance and were extensively discussed in the literature. First, the restrictions on slow-light propagation, due to coherent multiple scattering that induces an attenuation of the transmitted intensity. [5][6][7][8][9][10][11][12][13][14][15][16]25,27,[41][42][43][44] Second, the appearance of extrinsic radiation losses. 6,28,29 In what follows, we characterize the eigenmodes of a W1 waveguide in presence of disorder and discuss their properties in the light of the existing studies. We show that the present method allows to actually compute the localized eigenmodes, to study their localization length and their radiation loss rates, and to characterize the disorder-induced transmission losses. In addition, the method provides new insight into the physical mechanisms governing radiation losses. In particular, we show that extrinsic radiation losses in W1 waveguides are mostly determined by hybridization of the spatiallyeven and spatially-odd guided Bloch modes -a fact that has great impact on the design and optimization of high-Q optical cavities based on "gentle confinement" of W1guided modes. [59][60][61] We consider a W1 line-defect waveguide, consisting in one missing row of holes in a triangular lattice of cylindrical air holes of radius r, with lattice parameter a and slab thickness d. The W1 waveguide is sketched in Fig. 1(a). The modes E kn (r) of the regular structure are computed by means of guided-mode expansion. 28 The waveguide is oriented along the x-axis, while the z-axis is perpendicular to the slab. We adopt a rectangular supercell of width a along x and height ba √ 3/2 along y, where b is an integer, as sketched in Fig. 1(a). Modes that are strongly confined into the guide core, are scarcely affected by the superperiodicity along y introduced by the supercell. 56 We can thus safely restrict to modes having k y = 0, corresponding to an effective one-dimensional system whose bands are characterized by Bloch momentum k x along the x-axis. For simplicity, in what follows we denote k x simply as k. For the guided-mode expansion, we restrict to the lowest guided mode of the effective slab, having TE polarization. We thus neglect the TE-TM mixing induced by the hole lattice -an approximation that holds very well for small d/a, when the frequency of the next guided mode falls outside the band gap of the 2-D lattice. 28 Calculations can be easily generalized to an expansion on several guided modes if needed, by a straightforward extension of the formalism. As usual, a scaling relation with respect to the lattice parameter a holds and we can express all lengths, momenta and frequencies in dimensionless units. For the calculations that follow, we set d/a = 0.5, r/a = 0.3, b = 10, and we consider an air-bridge structure with dielectric constant ǫ 1,3 = 1 in the cladding layers, and ǫ 2 = 12 for the slab material, corresponding to an effective-slab dielectric constantǭ 2 = 8.77 for the calculation of the guided modes (see Appendix). In the calculations we use Ho's method 56 for the Fourier transform of the dielectric profile, and we truncate the expansion at G max = 3 (in units of 2π/a), corresponding to a basis of 229 guided modes. Fig. 1(c) shows the band structure in the region of the guided mode. Bands labeled 1 to 4 will be included in the truncated band set used for the disorder calculations. 80 Band 1 is the main guided band, often referred to as index-guided or spatially-even band (with respect to theσ kz mirror plane), while band 2 is the gapguided, spatially-odd band.
We model disorder as a gaussian fluctuation of both hole radii r m and positions (x m , y m ), where m runs over all holes of the extended structure, as sketched in Fig.  1 m ) are respectively the radius and the positions of the holes in the regular structure. The fluctuations δr m , δx m , and δy m are Gauss-distributed with standard deviations σ r , σ x , and σ y . This "minimal" disorder model -in particular the radii fluctuations -is expected to catch the main effects of hole shape fluctuations on the spectral properties 4,9,10 , at a low numerical cost in computing the disorder matrix elements V kn,k ′ n ′ , as shown in the Appendix. More refined disorder models, possibly accounting for the roughness of the hole contour 6,12,65,66 , could be easily implemented. As discussed above, the converged expansion (4) corresponds to the exact solution of the Maxwell problem and would thus automatically include local field effects in proximity of the hole sidewalls, whose importance has recently been debated. 31,65,66 In what follows, we restrict our analysis to the situation in which σ x = σ y = σ r = σ. Although we didn't carry out a systematic study, our simulations seem to confirm that radii fluctuations are the dominant effect, while position fluctuations give a less significant contribution. We consider values of σ ranging from 0.001a to 0.008a, corresponding to state-of-the-art structures. 4,10 For the numerical simulations, we consider W1 structures of length L = 1024a, namely including 9216 holes for the supercell height b = 10 considered here. This corresponds to assuming a uniform grid of 1024 k-points in the interval [−π/a, π/a] for the eigenvalue problem (7). A single disorder realization of a system of this length requires a few hours on a desktop computer. Computing Maxwell eigenmodes of such an extended structure with other first-principle methods, like e.g. a finite-element sampling of real space, would constitute a practically untractable problem. The advantage of the present approach is therefore striking.
In order to study the spectral properties of the system, a natural choice would be to introduce the quantity where H β (k, k y = 0) is the Fourier transform, taken at k y = 0 (and z = 0) of the magnetic field H β (r) of the β-th eigenmode, and γ β is the rate of radiation losses computed according to (11). Eq. (13) is a scalar Green's function of the electromagnetic field, corresponding to the trace of the Green's tensor of Maxwell equations. 67 It is evaluated at equal initial and final Bloch momenta (forward propagation), thus providing information about the transmission amplitude and, ultimately, about the spectral density of the field. 39 It is expressed through the resolvent representation 67 based on the magnetic field modes H β (r). These modes constitute a complete orthonormal basis with the standard measure for the scalar product, as opposed to the dielectric measure ǫ(r) needed to define a scalar product for electric field modes. 55,62 They are related to the corresponding electric field modes as H β (r) = −i(c/ω β )∇×E β (r). In Eq. (13) however, the choice k y = 0 corresponds to integrating the field along y, thus missing all spectral features related to spatiallyodd modes with respect to theσ kz mirror symmetry. In order to be able to visualize these spectral contributions, we prefer to adopt a different choice and introduce the quantity where H β (k, y) is obtained by Fourier transforming the field H β (r) with respect to x only. Eq. (14) represents the Green's function for the propagation of the field along x, at a fixed position y. In this way, by evaluating (14) at y = 0, i.e. displaced from the symmetry planeσ kz , spectral signatures of both even and odd modes will appear. For all plots, we choose y = 0.067a. Fig. 2(a) is a log-scale intensity plot of the spectral density A(k, ω) = −Im(G y (k, ω)) for σ = 0.002a, obtained by diagonalizing Eq. (7) on the subspace spanned by bands 1 and 2. The dispersion curves of the two bands for the regular structure are superimposed, as well as the free-space light dispersion. This latter defines the boundary of the light cone: below it, modes of the dielectric slab are fully confined by total internal reflection and ideally lossless. The spectral width along band 1 shows an abrupt increase above the light cone, due to the onset of the intrinsic radiation loss mechanism. A similar abrupt change is visible along band 2 at about ka/2π = 0.4. On the left of that point on the plot, the frequency of band 2 becomes degenerate with points of the same band lying inside the light cone, thus maximizing the mode mixing due to disorder and the subsequent extrinsic radiation losses, as already suggested by Kuramochi et al. 7 As a consequence, the further increase of the spectral width across the light cone boundary is less abrupt for this band, but still noticeable. Figs. 2(b)-(e) are linear intensity plots of the same quantity in the vicinity of the bottom of the guided band 1, for σ = 0.001a, 0.002a, 0.004a, and 0.008a, respectively. A feature common to the four cases is the presence of sharp spectral signatures below the bottom of the band of the regular structure. These features are extended along k, and extend to increasingly lower frequencies as σ is increased. They originate from spatially localized modes induced by disorder, and contribute to form the well known Lifshitz tail in the density of states below the band of the regular crystal. [68][69][70] These features have been experimentally characterized by Le Thomas et al. 42 and the corresponding sharp features in the transmission spectra have been experimentally observed and theoretically studied by modeling the optical response. 10,12 The density of states can be computed by averaging over several statistical realizations of the disorder configuration. Fig. 3(a) shows the density of states obtained by averaging over 500 realizations, for σ = 0.002a, compared to the density of states of the regular structure (panel (b) shows the spectral density for reference). The Van Hove singularity of the regular crystal is smeared over an interval of the order of ωσ/(2πc) and the Lifshitz tail is clearly visible.
The main innovation brought by the present approach is the possibility to compute the field profile of the actual Maxwell eigenmodes. In connection to the study of disorder, this makes it possible to study the localization properties and brings insight to the question whether Anderson localization of light can occur in PHCs. 2,10,[18][19][20][21][22][23][24] We plot the electric field profile of three selected modes of a disorder realization with σ = 0.002a in Fig. 4(a)-(c) respectively. The three selected frequencies are indicated by dotted lines in Fig. 2(c). The main plot in each panel represents the value of the field intensity at the center of each elementary cell of the guide (in this way the main envelope can be plotted over the full length of the simulated structure), while the insets show the full electric field profile for selected representative regions of the structure. Fig. 4(a) shows the ground eigenmode at ωa/2πc = 0.27267, which lies well below the band edge and is localized over a few elementary cells of the structure. Fig. 4(b) shows a mode at ωa/2πc = 0.27283, lying very close to the band edge and displaying a less localized character with few main lobes. Most of the modes computed in this spectral region have a similar character, with several lobes and an overall envelope with exponen-  Fig. 2(c)). The main panels display the field intensity computed at the center of each elementary cell, while insets show the full field profile for selected portions of the W1 waveguide.
tially decaying tails, as expected for Anderson localized states. Modes with several lobes have been studied in various kinds of PHCs and are considered as precursors of the extended necklace states. 10,23,24,42,71,72 Finally, Fig.  4(c) shows a mode at ωa/2πc = 0.27518, well within the energy range of the guided band, which is delocalized over the whole length of the simulated sample while still showing irregular amplitude fluctuations induced by disorder. This is expected as in the spectral region of this mode, the localization length is much longer than the simulation length L considered here. The radiation loss rate γ β for each eigenmode β has been computed using Eqs. (11) and (12). They are plotted, for one disorder realization, in Fig. 5(a), for σ = 0.001a, 0.002a, 0.004a, and 0.008a. Panel (b) in this figure shows the band structure and the light line for reference. Let us focus first on frequencies below 0.287, for which the guided band of the regular structure lies outside the light cone. For this frequency range, radiation losses are extrinsic, as already discussed. For each given value of the disorder amplitude σ, the rates γ β behave quite regularly, except for the lowest frequencies in the Lifshitz tail of the spectrum. There, the values of γ β are significantly scattered over about one decade, as one would expect for modes that are strongly local- ized by the random fluctuations of the underlying crystal structure. A similar trend, with strong variation of the loss rates close to the band edge, was obtained in Ref. 29, although the perturbative approach to disorder adopted in that work cannot predict the fluctuations of γ β . These large fluctuations are typical of states in the localized part of the spectrum, namely where the localization length is much smaller than the simulation length L. Similar fluctuations are observed in the transmission spectrum 9,10,12 and are more generally expected for all properties of localized eigenmodes, as in the localization lengths discussed below. Above the fluctuating region, the rates γ β are approximately constant over a wide frequency range, until they dramatically increase when approaching the boundary of the light cone. The rates γ β , for the values of the disorder amplitude considered here, are in general very small, as expected for an extrinsic process. When the disorder amplitude σ is doubled, the extrinsic rates correspondingly increase by a factor of four. This dependence has already been discussed in the literature. 11 It can be explained by considering that the overlap integral entering Eq. (11) depends linearly on the disorder matrix element V to leading order. This can be inferred from Eq. (A9), if we consider a perturbation expansion of mode E β where the zeroth order term is a guided Bloch mode, that does not contribute to the overlap integral (A9) due to orthogonality. For the frequency range above the light line, the eigenmodes of the disordered structure are characterized by much higher loss rates that are now mostly of intrinsic nature and no longer depend on σ. 81 In the frequency range between 0.294 and 0.306, rates of modes related to the spatially-odd band 2 are also present. Here losses show again a partially extrinsic character, especially in the lowest frequency part where a quadratic dependence on σ appears. For the discussion below, it is also important to keep in mind that rates in this frequency range are significantly higher than those computed at the bottom of band 1.
The origin of the radiation loss rates can be inferred from the nature of the computed eigenmodes. These, according to Eq. (4) are linear combinations of Bloch modes at different k and different bands. This hybridization involves all modes of the regular structure, regardless of their frequency spacing ∆ω -the mixing fraction being quantified by V (ω/∆ω) to leading order. Fig.  6(a)-(c) shows the convergence of the eigenfrequencies of four selected eigenmodes, at frequencies corresponding to band 1, as a function of the number of bands included in the diagonalization. The convergence is excellent already when only band 1 is included in the Bloch-mode expansion. Convergence of the loss rates, shown in Fig.  6(d), requires instead that the spatially-odd band 2 is also included in the calculation. In particular, for the three lowest frequencies lying in the region of extrinsic radiation losses, the loss rates increase by as much as two decades when band 2 is included. The reason is that the main source of extrinsic losses, close to the edge of band 1, is the hybridization with Bloch modes of the spatially-odd band 2. This mixing would be forbidden by symmetry in the regular structure and is now allowed as disorder breaks theσ kz invariance. The fourth selected frequency instead corresponds to an eigenmode for which the intrinsic radiation loss mechanism is dominant, and the hybridization with band 2 doesn't thus produce any significant effect. Fig. 6(e) shows the components of the lowest computed eigenmode β = 1, |U β (k, n)|. Apart from the dominant contribution of Bloch modes close to k = ±0.5 of band 1, the main contribution comes from band 2 that shows about 3 × 10 −3 admixture at the edge of the Brillouin zone. This explains the dominant contribution of band 2 to the extrinsic loss rates, seen in Fig. 6(d). Band 2 losses are several orders of magnitude higher than the extrinsic losses computed when accounting for band 1 only. Therefore, even the slightest admixture of the two bands -here about 0.3% -can completely determine the final loss rate at the band edge. The role of the mixing between spatially-odd and even bands has been already pointed out by Kuramochi et al. 7 in the framework of perturbation theory, thus restricted to the frequency region where the two modes are degenerate. Our result generalizes the one by Kuramochi et al., by showing that the mixing determines all extrinsic radiation loss rates of the main guided band of a W1 waveguide. The importance of this finding goes well beyond the study of disorder effects. It suggests that radiation losses in all W1-based systems can be reduced by studying designs that minimize the mixing between bands 1 and 2, e.g. by displacing band 2 at higher frequencies. It can indeed be argued that the same mechanism limits the Q-factor of cavities based on the gentle confinement of the W1 guided mode. [59][60][61] Our conclusion thus shows a new direction to take for the optimization of such cavity designs.
To conclude the analysis of the W1 waveguide, we study the localization lengths of the simulated eigenmodes. The localization length in a disordered system can be defined in several different ways. 39 One natural definition, when studying eigenstates, is the inverse participation number (IPN). For a wavefunction φ(x) in 1-D, the IPN is defined as I = 1/ |φ(x)| 4 dx, and has the dimension of a length if |φ(x)| 2 is normalized to 1. The IPN measures the cumulative length over which the field intensity is significantly large. For very singular shapes of the disorder profile, the IPN can differ considerably from the exponential decay length of the tails of φ(x) - which is another measure of the localization length. 39 For an eigenmode H β (r) of our system, we thus define where, in order to extract the mode length, we take the integral along the guide core at y = 0. A second characteristic length related to disorder, typically studied in the context of PHCs, [8][9][10]12,22,23 is the exponential attenuation length of the transmission of a monochromatic beam. For this, we introduce the quantity which is proportional to the transmission of a monochromatic source from x to x ′ . The frequency-dependent exponential attenuation length is defined as the decay constant L ω in T (x, x ′ , ω) ∝ exp(−|x ′ − x|/L ω ). In the numerical calculation of both L ω and I β , we checked the convergenge of the computed values as a function of the system size L, up to the largest value of L considered, in order to detect the onset of finite-size effects. Fig. 7(a) shows the quantities I β and L ω , as computed for the W1 guide under study. Lines denote configuration averages over 500 disorder realizations, while the circles are IPN values for the eigenmodes of a given realization. The IPN is more affected by the finite size of the simulation and saturates at about I β = 250a, while the saturation of the attenuation length is not visible on the scale of the plot. This is expected, as the IPN is an extensive property of the eigenmode profile, contrarily to the attenuation length which can be extracted from the slope at a single point. When finite-size effects are negligible, the two quantities roughly coincide. We point out that the attenuation length includes the decay due to radiation losses, while the IPN is an exclusive property of the eigenmodes of Eq. (7). In the frequency range of Fig. 7, the computed losses correspond to a negligible contribution to the attenuation length of the transmission, which is thus uniquely determined by the coherent multiple scattering in the propagation along the guide.
As for other properties of localized eigenmodes, the IPN displays large fluctuations at each given frequency. This suggests that, independently of the average localization length at a given frequency, rare necklace modes can exist with a profile that extends over the whole sample length L. 10,23,24,42,71,72 Localization in one-dimensional photonic structures has been reported in several experiments. 10,18,19,21,23,24,42 The present result sheds light on the question concerning the nature of light localization in these structures.
We have already excluded the extrinsic radiation losses as a cause of attenuation in a W1 waveguide close to the bottom of band 1. It is clear that the origin of localization in our model calculations is to be found in coherent multiple scattering processes, corresponding to the Bloch-mode superposition in expansion (4). This is of course expected, as in one-dimensional disordered systems localization is always present, 30,39 and the W1 waveguide shows to be very closely one-dimensional. Indeed, the expansion on a single 1-D Bloch band was sufficient to obtain well converged eigenmodes, while the very small admixture with band 2 was only necessary to account correctly for radiation losses. A relevant question, raised by Vlasov et al., 22 is whether localization is of the Anderson kind or rather due to a random local energy shift of the band edge induced by large disorder fluctuations. This second mechanism can be pictured as the random alternance in space of pass-band and stop-band regions, and is thus enforced by the Bragg diffraction of the underlying periodic structure, rather than by scattering in an effective medium as in the Anderson localization case. Our approach gives direct access to the spectral function. It shows that, for the disorder amplitudes considered here and corresponding to state-of-the-art samples, 4,10 the band-edge shift is much smaller than the width of the stop band (see Figs. 1 and 2). According to Vlasov et al. then, at least for the disorder model considered here localization should be of the Anderson kind. Comparing the behavior of the IPN and the attenuation length below the band edge, as shown in the inset of Fig. 7(a), supports this conclusion. For decreasing frequency, the attenuation length ceases to follow the average IPN and decreases down to its smallest possible value, essentially given by the unit cell length a. For frequencies below this crossover, the attenuation length is no longer determined by the localization length of the eigenmodes, but rather by Bragg diffraction on the PHC, and a similar attenuation length would arise within the stop band also in the case of a regular crystal. The IPN, on the other hand, seems to approach a constant value of a few unit cells, that corresponds to the smallest IPN taken by Anderson localized modes for the given value of σ, as already observed when directly visualizing the mode in Fig.  4(a).

IV. APPLICATION TO A L3 COUPLED-CAVITY WAVEGUIDE
As an additional test of our method, we consider a chain of L3 coupled cavities, as sketched in Fig. 8(a). The system differs from the W1 waveguide by the presence of one hole every four along the waveguide. This particular kind of coupled-cavity waveguide has recently been experimentally characterized, 58 in order to determine the origin of out-of-plane radiation losses and the limiting factors to slow light propagation. More generally, coupled-cavity waveguides have been the object of various experimental 58,73,74 and theoretical 13,63,64,75,76 studies, because of their potential as effective slow-light channels, thanks to the very flat band expected for an impurity chain.
Here, we start our analysis from the Bloch modes of the regular structure, that we compute as for the W1 waveguide using guided-mode expansion 28 , by assuming the supercell sketched in Fig. 8(a). The disorder model and all other parameters are the same as for the W1 waveguide studied above. In particular, the same length L = 1024a, corresponding to 256 L3 cavities, is assumed. With these parameters, the band of interest 58 goes from ωa/2πc = 0.2766 to ωa/2πc = 0.2779, as plotted in Fig.  8(b). Calculations using our method fully converge already when including only this band in the Bloch-mode expansion. This holds also for the radiation loss rates, contrarily to the case of the W1 waveguide, because of the dominant role of the intrinsic radiation loss mechanism as discussed below. The spectral function for σ = 0.002a is plotted in Fig. 8(b), while Fig. 8(c) displays the radiation loss rates computed for the field eigenmodes, for the disorder amplitude varying from σ = 0.001a to σ = 0.008a. The Brillouin zone, four times smaller than for a W1 waveguide, now lies completely above the light cone, and radiation losses are mostly determined by the intrinsic mechanism predicted for the regular system. 13,75 In particular, the rates vary from a value lower than that of the single L3-cavity, 28,77,78 at the band bottom, to a higher one at the band top, as predicted in the absence of disorder from the 3-D solution of Maxwell equations 13 or from a tight-binding approach. 75 . The dependence on the disorder amplitude σ clearly shows that, up to σ = 0.004a, the extrinsic radiation loss mechanism takes over the intrinsic one only in the vicinity of the band bottom. For larger σ instead, the extrinsic mechanism becomes relevant over the whole band, and large fluctuations of the rates appear at each given frequency. Fig. 9 shows the electric field intensity computed for three selected eigenmodes, denoted by horizontal lines in Fig. 8(b). The ground eigenmode in Fig. 9(a) is localized over the length of a few L3-cavities. Modes at higher frequencies displayed in panels (b) and (c) are progressively more delocalized, again as expected for a one-dimensional disordered system. Notice that each mode has the shape of a slowly varying envelope function modulating the field profile of the ground mode of a L3 cavity, 77 as in tightbinding calculations. 75 As for the W1 waveguide studied in the previous Section, the modes displayed in Fig. 9 are computed without accounting for losses [which are computed at a later stage through Eq. (11)]. Contrarily to the W1 waveguide however, here the attenuation length corresponding to the computed radiation loss rate is actually shorter than the localization length of the modes, when considering the least localized modes at the band center. Hence, for these modes the result of our Bloch-mode expansion is not fully consistent with the initial assumption of lossless Bloch modes. To better understand this issue, we show in Fig. 10(a) the localization lengths computed from the IPN and the attenuation length computed from the decay of the transmission (16). Fig. 10(b) shows the spectral function for reference. The attenuation length is evaluated both by accounting for the computed loss rates γ β (black curve) and by assuming a vanishing loss rate for all modes (red curve). This latter coincides with the average IPN over a wide range of frequencies close to the band extrema, while at the band center the IPN starts being affected by the finite length of the simulated structure, similarly to the case of the W1 waveguide. The attenuation length including losses is instead dramatically suppressed, except for a narrow region at the band bottom (see inset). We argue that, in the region where the two curves differ significantly, the localized eigenmodes computed here are actually ill-defined. Radiation loss rates of Bloch modes are so high that they should be accounted for already in the multiple scattering process from which the eigenmodes of the disordered system arise. This could be achieved by replacing the guided-mode expansion with a full expansion on guided plus radiative modes of the slab, or simply with an expansion on 3-D plane waves. The inset of Fig. 10(a) shows a detail of the region below the band edge, displaying the same trend already observed for a W1 waveguide in Fig.  7(a), although now the minimal IPN value is larger than in the W1 case, as expected from a longer unit cell. We conclude by making some considerations on the influence of disorder on slow-light propagation, both in the W1 and in the coupled-cavity waveguide. Slow-light and is plotted as a guide to the eye. Blue: W1 waveguide. Red: L3 coupled-cavity waveguide neglecting radiation losses. Green: L3 coupled-cavity waveguide accounting for radiation losses.
has been one of the main driving forces of the development of one-dimensional PHC structures. 25,27,44,79 Ideally, close enough to the bottom of the guided band, the group velocity of a narrow light pulse, defined as v g = dω k /dk, can become arbitrarily small. Disorder hinders this ideal behavior by producing a strong attenuation of the transmitted signal, as discussed above. Past experimental studies 5,7,12,42 have focused on the dependence of the attenuation length on the group index, defined as n g = c/v g . All recent theoretical studies 6,9,10,12,14,16 show that the dominant functional dependence of the attenuation length close to the band bottom is on n −2 g . This behaviour is induced by disorder through the mechanism of coherent backscattering, determined by second and higher-order coherent scattering processes, which is a precursor of localization. 39 The n −2 g -dependence was confirmed in several experiments. 5,7,11 It was also pointed out 9,10,16 that, very close to the band bottom, the onset of localization makes multiple scattering beyond second order important and the n −2 g -dependence consequently breaks down for high n g , with the attenuation length that levels off to its minimal value given by the lattice parameter a. In Fig. 11, we plot the computed atten-uation length L ω for σ = 0.002a as a function of the nominal group index n g , both for the W1 and for the coupled cavity waveguides. By "nominal", it is meant the group index that would result at frequency ω from the dispersion of the ideal band. The plot thus does not include frequencies below the band bottom. The result for the W1 waveguide (blue) reproduces the result obtained by Mazoyer et al. 9,10 For the coupled cavities, we display two curves corresponding to the case including radiation losses (green), and to the unphysical but illustrative lossless situation (red) respectively. This latter coincides with the W1 result, except for the largest values of n g considered, that arise close to the band bottom. The interpretation of this result is straightforward. A small group index corresponds to a frequency well above the band bottom, where eigenmodes have a localization length much longer than the elementary cell. In this case, independently of the length of the elementary cell, the envelope function of the eigenmodes will be localized according to the John mechanism, 2 with the localization length solely determined by the curvature of the Bloch-mode dispersion, thus by the group index. When approaching the band bottom instead, the localization length becomes smaller and eventually reaches its lower bound corresponding to a small number of elementary cells, as shown from the analysis of the IPN. In this limit, the localization length of the coupled-L3 system is four times longer than that of the W1 waveguide, at corresponding n g . This limiting behavior can be inferred by comparing the upper branch of the red curve (corresponding to the band bottom) and the blue curve in Fig. 11. When accounting for intrinsic radiation losses in the L3-coupled system, the attenuation length is dramatically reduced except for the region close to the band bottom, where losses are the smallest and disorder concurs significantly in determining the attenuation. Following Mazoyer et al., 9,10 , a diagram of the kind displayed in Fig. 11 sets the maximum length scale over which light propagation can still be ballistic, and thus the slowdown predicted for the ideal system can still occur. Then, from this analysis we can conclude that, in spite of the increased loss rates, a L3-coupled waveguide can still bring an advantage over the W1 waveguide when close to the band bottom, by increasing four times this maximum length for a given (large) group index.

V. CONCLUSIONS
The main result of the present analysis is the demonstration that Boch-mode expansion can be an extremely effective and reliable tool for studying the properties of the eigenmodes of Maxwell equations, in PHCs where disorder or small systematic variations of the otherwise periodic structure are present.
The approach relies on a method to compute the Bloch modes of the regular PHC. Ideally, this first elementary step should be accomplished by exact diagonaliza-tion of the Maxwell problem within a single elementary cell, for example by using a finite element approach. For the scope of the present work however, we preferred to adopt the guided-mode expansion method, 28 because of its reliability at a minimal computational cost. In this case, however, out-of-plane propagating modes are not included in the expansion, and the radiation losses of the final Maxwell eigenmodes must be computed non selfconsistently within a perturbative approximation at a later stage.
The general advantage of Bloch-mode expansion lies in the optimal choice of the Bloch-mode basis, for situations where the variations from the periodic structure are small. Then, an expansion restricted to a very limited number of spectrally close Bloch bands, can produce excellent convergence, thus practically reproducing the exact solution of Maxwell equations at a reasonable computational cost.
As an example, we have applied the method to study the effect of disorder in a W1 line-defect waveguide and a L3 coupled-cavity waveguide. The method gives direct access to the shape and spectrum of the eigenmodes, from which the localization properties, and the mechanism underlying radiation losses have been characterized.
One very important result of this work -apart from the assessment of the Bloch-mode expansion method -is the explanation of the origin of radiation losses of the spatially even guided modes of a W1 waveguide close to the band bottom. Radiation losses are mainly determined by the hybridization of these modes with those of the guided band of opposite spatial parity. A very small admixture of the two bands is sufficient to increase the radiation loss rate close to the band edge by two decades. Hence, the effect is dominant in spite of the frequency separation of the two bands that hybridize -a feature overlooked by previous analyses. This result was only possible thanks to the separate analysis of the contribution of different Bloch bands. It unveils the physical mechanism limiting the quality factor of all structures based on a W1 waveguide. In particular, all designs of high-Q cavities based on slow modulation of a W1 waveguide should be affected. We suggest that a significant increase in the quality factor of such structures can be obtained by studying a waveguide design that reduces the hybridization, e.g. by maximizing the frequency separation between the two guided bands.
The Bloch-mode expansion has a very broad spectrum of potential applications. These range from the realistic modeling of high-Q cavities, to the study of radiationmatter coupling and cavity quantum electrodynamics effects for systems of several quantum dots coupled to the modes of PHC cavities or waveguides. All these application should take advantage of the possibility of computing the actual eigenmodes of Maxwell equations, thus making the Bloch-mode expansion an election method for the modeling and design of photonic systems.
In this work we give explicit expressions only for guided modes having TE polarization, as we assume quasi-TE PHC modes in the results that we present in Sections III and IV. More general expressions, for modes where TE and TM polarizations are mixed, can be obtained in a straightforward way by following the guidelines in Ref. 28.
For the calculation of the matrix elements V kn,k ′ n ′ , only the expression of E kn (r) for r in the slab is needed, as we assumed δǫ(r) = 0 in the claddings. By further assuming a structure symmetric with respect to a mirror reflectionσ xy , namely equal upper and lower claddings, the guided-mode expansion reads where µ = (k + G, α) is the global index for the guided mode, ρ is the in-plane position vector and z the coordinate orthogonal to the plane, q µ = ǭ 2 ω 2 µ /c 2 − g 2 , and S is a normalization area that will disappear from the final result. The complex coefficient A 2µ is defined in Ref. 28. We have introduced the unit polarization vector ǫ g = z × g of the guided mode, where z is the unit vector along the z-axis. By replacing (A2) into (8), we obtain V kn,k ′ n ′ = µ,µ ′ c * n (µ)c n ′ (µ ′ )δǫ(g ′ − g)M µµ ′ ,