Theoretical and computational concepts for periodic optical waveguides

: We present a general, rigorous, modal formalism for modeling light propagation and light emission in three-dimensional (3D) periodic waveguides and in aggregates of them. In essence, the formalism is a generalization of well-known modal concepts for translation-invariant waveguides to situations involving stacks of periodic waveguides. By surrounding the actual stack by perfectly-matched layers (PMLs) in the transverse directions, reciprocity considerations lead to the derivation of Bloch-mode orthogonality relations in the sense of E × H products, to the normalization of these modes, and to the proof of the symmetrical property of the scattering matrix linking the Bloch modes. The general formalism, which rigorously takes into account radiation losses resulting from the excitation of radiation Bloch modes, is implemented with a Fourier numerical approach. Basic examples of light scattering like reflection, transmission and emission in periodic-waveguides are accurately resolved.


Introduction
In recent years, the race for faster optical telecommunication and data processing has motivated research towards solutions to minimize the involvement of electronics in signal manipulation.In the quest for ultimate miniaturization, circuits relying on high refractiveindex contrast waveguides, like semiconductor-air/oxide photonic wires, appear as a promising approach [1].Alternatively to these systems relying on a conventional guiding in a high index core surrounded by a low index cladding, one can use a periodic structure, photonic crystals (PhCs) for instance, with high refractive index contrasts and a period of the wavelength order [2].These strong contrast and a well-chosen periodicity can create a photonic bandgap, and when a line defect is created in the crystal structure, monomode guidance is achieved.The PhC waveguides share many common features with their zinvariant photonic-wire counterparts [3] but, in addition, due to the periodicity along the propagation axis, they may welcome guided modes with dramatically slow group velocities [4] and anomalous dispersion, with potential applications for true all-optical signal processing [5,6].Thus, photonic-crystal waveguides and aggregates of them, as shown in Fig. 1(a), offer a robust platform for light-matter interaction phenomena, which is becoming more and more important in advanced integrated circuit designs.
In the analysis and design of integrated-optics circuits with classical z-invariant waveguides, modal concepts play a central role.The normal mode theory of bound-radiationleaky modes [7,8] provides a powerful formalism and allows a physical description of many important phenomena, such as light propagation in z-invariant sections, mode-coupling at discontinuities, or light emission in waveguides.The exploitation of this formalism has given birth to a large variety of numerical tools for z-invariant waveguides, see [9][10][11][12][13][14] and references therein.By contrast, a similar modal approach with Bloch mode is often ignored in the design and analysis of integrated components relying on periodic waveguides, like photonic-crystal microcavities or tapers [15][16][17].
Of course, semi-analytical methods relying on Bloch-mode transfer techniques [18-21] have already been developed but they are mostly restricted to the analysis of lossless 2D problems and thus do not tackle the important problem of radiation into the cladding.For 3D problems, softwares are available for calculating Bloch modes [22][23][24][25][26], and yet numerical tools that relies on a full sampling in space and time are often preferred to handle with aggregates of periodic waveguides.Although these methods may provide accurate results, they are mainly numerical and do not provide physical insight into the underlying physics of light propagation.Recent works have tackled the important problem of handling nonconservative periodic-waveguide aggregates with Bloch mode concepts in the frame of light emission [27] and extrinsic radiation losses [28,29] computations.To handle with the radiation losses, they use Born-like approximations and restrict the analysis to the bound Bloch modes without explicitly considering radiation Bloch modes.Finally and despite its importance, the works on computation and theory of Bloch modes still remains scattered and to our knowledge, no complete description of the formalism used to calculate the Bloch-mode modal reflection or transmission coefficients has already been given.
In this work, we present a general theoretical formalism for modeling light emission and propagation in periodic waveguides and in aggregates of them.Conceptually, the formalism can be seen as a generalization of the classical formalism developed for z-invariant waveguides.It allows to handle both the bound and radiation Bloch modes of periodic waveguides and by combining independent Bloch mode calculations of intermediate sections, it also allows a semi-analytical treatment of light propagation in aggregates of periodic waveguides.In Section 2, we establish the Lorentz reciprocity relation for periodic waveguide aggregates surrounded by perfectly-matched absorbers (PMLs).This relation is frequently referred to in Sections 3 and 4, where important properties, like the Bloch-mode orthogonality and the symmetrical property of the scattering matrix are derived.These derivations provide a theoretical background to the numerical aspects that are presented in Section 5. Finally, we illustrate several problems related to the calculation of light emission and propagation in periodic waveguides, especially focussing on the accuracy of the computational results.Section 6 summarizes the main results.

Lorentz reciprocity theorem in periodic waveguides aggregates
To illustrate our purpose, let us consider the geometry shown in Fig. 1(a).The geometry is composed of aggregates of periodic waveguides (possibly incorporating z-invariant waveguides), for which the claddings are assumed to extend toward infinity in the transverse x-y directions.If the material is lossless, at any frequency, every periodic section (as that shown in Fig. 1(c) supports a finite set of bound modes and a continuum of radiation modes.These are the so-called normal Bloch modes that are pseudo-periodic functions of the zcoordinate.As for classical z-invariant waveguides, there are two main difficulties related to normal-Bloch-mode concepts.The first one is the fact that the radiation modes form a continuum which is difficult to take into account in numerical methods.The second one is much more problematic: indeed for open periodic systems that extend toward infinity in the transverse directions, the radiation Bloch modes of the waveguide carry an infinite power and only linear combination of these modes may represent physical quantities, see the related discussion for radiation modes in z-invariant waveguides in [7,8].To remove these theoretical difficulties, mathematicians consider analytical continuations by making complex some variables that are usually real, like the wave frequency for instance [30].The same approach is adopted hereafter and instead of considering the scattering problem of the actual geometry in Fig. 1(a), we will consider that of the associated geometry shown in Fig. 1(b).The new computational system is obtained by surrounding the actual geometry with PMLs.We adopt here the approach in [31,32] and consider that the PML can be viewed as a linear complex coordinate stretching.Within this picture, the solution of Maxwell's equations remain unchanged except in the PMLs where the cladding permittivity and permeability distributions are renormalized by constant complex numbers.These complex numbers actually define the analytical continuation used to satisfy the outgoing wave conditions in the claddings.Because all the electromagnetic-field components fall off exponentially in the PMLs, we will assume that, after propagating over a given finite propagation distance, these components are negligible.This allows to define an outer boundary for the PMLs and to deal with finite transverse grid volume.The PML thickness, not represented in Fig. 1(b), depends on the specific stretching coefficients chosen for the numerical implementation.The new computational system of Fig. 1(b) largely differs from that of Fig. 1(a).We are now facing a closed system and not an open one and even if the system in Fig. 1(a) is composed of lossless materials, the new system of Fig. 1(c) is dissipative because of the complex continuation so that each radiation Bloch mode carries on a finite power.Moreover, it implies that the modes of any periodic section of the computational system, see Fig. 1(d), are quantized (there is an enumerable set of modes, not a continuum) and, except for the truly guided modes, they depend on the specific analytical continuation, i.e. on the choice of the PMLs used for the calculation.Hereafter, they will be referred as quasi-normal Bloch modes (QNBMs) to emphasize that they are not identical to the normal Bloch modes of the actual periodic waveguide.They are conceptually similar to the quasi-normal modes used for solving Maxwell's equations in aggregates of z-invariant waveguides [11][12][13][14].
We now derive the explicit form of the Lorentz reciprocity theorem [7] for the computational system of Fig. 1(b).This theorem will be repeatedly used in the remaining.In the orthogonal cartesian coordinate (x, y, z) system and in the presence of a dipole source J, the curl Maxwell equations in the Gaussian system of units are where H and E are the magnetic and electric fields at a given frequency ω, ε and µ are the relative permittivity and permeability tensors, δ is the Dirac distribution and j 2 = -1.In the following, we assume that the materials of the waveguide and of the cladding are reciprocal and possibly anisotropic.Thus everywhere we have µ = µ T and ε = ε T , where the upper-script T denotes matrix-transposition.We consider two solutions (labeled by the subscripts 1 and 2) at two frequencies ω 1 and ω 2 for the same permittivity and permeability distributions, Applying the Green-Ostrogradski formula to the vector E 2 × H 1 on a closed surface S of volume V enclosing the sources J 1 and J 2 , one gets By subtracting to Eq. ( 3) the related relation obtained by exchanging the indices 1 and 2, we further obtain In the following, we will consider closed surfaces formed by two waveguide cross-sections, A 1 and A 2 defined by z = z 1 and z = z 2 > z 1 , and by the outer-boundary-PML surface S' that closes the two cross-sections.Because of the PML attenuation, the field can be assumed to be negligible on the outer-boundary-PML surface so that the integral on S' is null and the surface integral in the left-side of Eq. ( 4) reduces to two surface integrals over A 1 and A 2 .
#  4) Equation ( 5) represents the Lorentz reciprocity relation applied to a waveguide (or more generally to an aggregate of waveguides) surrounded by PMLs.Since it will play a central role hereafter, we conveniently introduce a few notations for the sake of simplicity.Denoting by Φ 1 = |E 1 ,H 1 > and by Φ 2 = |E 2 ,H 2 > the vector formed by the 6-components of the electric and magnetic fields, we define F z and E V define two continuous bilinear forms of Φ 1 and Φ 2 that are respectively antisymmetric and symmetric.With these notations, Eq. ( 5) is conveniently rewritten Importantly, let us note that for ω 1 = ω 2 and in the absence of any source (J 1 = J 2 = 0), the term on the right side of Eq. ( 7) is null and thus F z (Φ 1 , Φ 2 ) is independent of z.In this case, it will be simply denoted F(Φ 1 ,Φ 2 ).

Quasi-normal Bloch mode properties
In the analysis of light scattering in aggregates of periodic waveguides, the QNBMs of every periodic sections play a key role.In this Section, we aim at deriving the main properties of these modes like their orthogonality and their group velocity.The whole set of radiation and bound modes is studied but a special attention is paid to bound QNBMs as they govern most of the physics in waveguide sections sufficiently far from any source of excitation (spatial steady state).

Quasi-normal Bloch modes
For a given frequency ω, the forward-propagating QNBMs are potentially constituted of a finite number of truly guided modes for lossless materials and of an enumerable set of other modes.For the sake of simplicity, we assume that the waveguide is monomode and we denote the forward-propagating guided mode by Φ (1,ω) =|E (1) ,H (1) >exp(jk 1 (ω)z) with Im(k 1 )=0 .
Similarly we denote by Φ (m,ω) =|E (m) ,H (m) >exp[jk m (ω)z], with m > 1 and Im(k m ) > 0, the other modes.Apart from the bound mode, the set of forward-propagating QNBMs (m > 1) potentially encompasses leaky QNBMs that are modes operating below the cutoff [23,33] and also a finite number of evanescent truly guided modes with a complex propagation constant k, Real(k) = π/a (modulo 2π/a) [17].Although the latter do not carry energy, they have to be taken into account in any calculation.Because of the Bloch theorem, the QNBMs are pseudoperiodic functions of the z-coordinate with a being the waveguide period.In the Annex 1, we show that it is possible to associate a backward-propagating mode Φ (-m,ω) =|E (-m) , H (-m) >exp[-jk m (ω)z], with E (-m) and H (-m) periodic functions of the z-coordinate, to every forward-propagating mode Φ (m,ω) .Without proof, we will assume that the forward-and backward-propagating QNBMs form a complete set and that, in any intermediate section of periodic-waveguide aggregates, the electromagnetic solution of a given scattering problem for a fixed frequency ω can be written as a linear combination of the section QNBMs.

Quasi-normal Bloch mode orthogonality
To derive the orthogonality of QNBMs, we apply Eq. ( 9) for ω = ω'.Since the term on the right side of Eq. ( 9) is null, we deduce that F z (Φ (p,ω) ,Φ (-q,ω) ) = 0, provided that exp[j{k p (ω)k q (ω)}a] ≠ 1.Therefore, we obtain the following orthogonality relation where δ q,p is equal to 1 for p = q, and 0 otherwise.F (p,ω) is a complex constant, which can be used for normalizing the QNBMs.In the next subsection, it will be expressed as a function of the QNBM group velocity and mode volume.Note that because of the presence of PMLs, the orthogonality relation of Eq. ( 10) is not defined with usual Poynting E×H* products and therefore the energy carried in any periodic section of the system cannot be decomposed as a sum of energies carried by individual QNBMs, even if the waveguide materials are lossless.Nevertheless, the orthogonality relation remains central in the analysis of light propagation into aggregates of periodic waveguides since it allows for closed-form expressions of the scattering reflection and transmission coefficients at the waveguide-section interfaces.

Group velocity of quasi-normal Bloch modes
As mentionned in the introduction, one of the most promising property of periodic waveguides is their ability to offer very slow group velocities.Classically, the group velocity of a bound mode Φ (1,ω) is defined as In the following, this definition is generalized to radiation modes, even if the group velocity looses its physical meaning: the propagation speed of the energy.
For p=q and for two different frequencies ω and ω', Eq. ( 9) becomes As the frequency ω' tends towards the frequency ω, Φ (p,ω') → Φ (p,ω) and since the bracketed term of the left-hand side of Eq. ( 11) tends towards ja where ω) , Φ (-p,ω) ) is the complex modal volume that is actually independent of the z origin of the unitary cell and v g can be complex for radiation QNBMs.In a general sense, ) p ( g v =dω/dk p is understood as a proportionality factor between the two most important quantities of QNBMs, F (p,ω) and E (p,ω) .Furthermore, by applying Eq. (3) to Φ (p,ω)  onto a unit cell in the absence of source, we obtain 0 = ∫∫∫ Cell [E (-p)T ε E (p) + H (p)T µ H (-p) ] dV, by noting that the surface integral on the left-hand side of the Eq. ( 3) is null because (E (p) × H (-p) ) is fully periodic with period a.Therefore, we have for any cell From a numerical point of view, it is important to note that a better accuracy for computing E (p,ω) is obtained with the integral on the right-hand side of Eq. ( 13) rather than with the one on the left-hand side.Indeed as one is generally dealing with non-magnetic materials (µ = 1), the magnetic fields H (p) and H (-p) , contrary to the electric ones, are continuous everywhere except at the inner PML boundary.Finally, note that since we are dealing with E × H products and not with E × H* ones because of the PMLs, the complex mode volume and the complex group velocity are not related to energy flow, except for the important case of truly guided QNBMs that we now consider.
To deal with this case, it is convenient to start with the "real" Bloch modes of the actual system of Fig. 1(c), i.e. without PMLs.We will denote by Φ (1,r) = |E (1,r) ,H (1,r) > exp(jk 1 z) and Φ (-1,r) = |E (-1,r) ,H (-1,r) > exp(-jk 1 z) the forward-and backward-guided Bloch modes propagating in the periodic waveguide at frequency ω, the upper-script 'r' being used to differentiate these modes from the corresponding truly guided QNBM Φ (1,ω) and Φ (-1,ω) of Fig. 1(d).Truly guided Bloch modes exist only if the waveguide and cladding materials are lossless.Since ε(r), µ(r) and k 1 are all real, it is easily shown from the conjugate form of the curl Maxwell's equations that | E (1,r) ,H (1,r) at least if the modes are non-degenerate.The truly guided QNBMs of the computational system of Fig. 1(d), Φ (1,ω) and Φ (-1,ω) , are identical to the real Bloch modes, Φ (1,r) and Φ (-1,r) , inside the inner PML boundaries, but differ in the PML regions.In practical situations, it is necessary to normalize the power flow or the volume energy of the truly-guided QNBMs.Indeed, since the latter are only known in a finite volume (inside the inner PML boundaries), it is not easy to reconstruct the missing field in real space outside the inner PML boundaries.
Actually, this problem can be fully released.In the annex 2, we show that PMLs allow for some conservation laws, and more specifically, that F (p,ω) and E (p,ω) are invariant quantities that do not depend on the specific choice of the PML.In particular, this implies that F (1,ω) and E (1,ω) , which are defined for a given PML implementation, are equal to their associated quantities defined for the actual Bloch mode obtained for a unitary PML stretching coefficient.In other words, we have ∫∫∫ Cell [E (1,r)* εE (1,r) + H (1,r)* µH (1,r) ]dV = E (1,ω) , ( where the integrals over the "real" Bloch modes on the left-hand sides extend over infinite transverse cross-sections, while those related to F (1,ω) and E (1,ω) on the right-hand side run over finite transverse grid sizes bounded by the PML outer boundary.In practice, Eqs.(15a) and (15b) can be used to calculate the power flow or the cell energy of a bound QNBM.They can also be used for normalization, a bound QNBM with a unitary power flux being obtained with F (1,ω) =4.Furthermore, Eq. (12) shows that the electric and magnetic fields of a normalized bound QNBM (or of a bound Bloch mode) scale as (a/ )

Scattering-matrices of periodic waveguides aggregates
In this Section, we assume that the QNBMs are all normalized, i.e.F (p,ω) =4 for any p.For this specific normalization, we show that important reciprocity relations hold for light scattering or light emission in aggregates of periodic waveguides.We especially focus on the symmetryproperty of the scattering matrix related to periodic-waveguide sections.Since the analysis is fully performed at a single frequency ω, we drop the ω-dependence in the notation of the QNBMs.

Symmetry property of the scattering matrix
Figure 2 shows a general scattering problem between two periodic waveguides.Hereafter, we will define a scattering matrix S that relates the QNBM amplitudes of the left and right periodic-waveguide sections, potentially incorporating coupling with dipole sources and we show that S = S T like in classical z-invariant waveguides.
Fig. 2. The two solutions used for deriving the S-matrix reciprocity relation.At a given frequency ω, the QNBMs of the left periodic-waveguide section (z<L 1 ), labelled by "L" for "left", are denoted by Φ (p,L) and those on the right side (z>L 2 ), labelled by "R" for "right", by Φ (p,R) , p being a relative integer.The geometries are arbitrary and may contain sources for the left and right ends of the structure and for L <z<L 2 .The whole system is assumed to be by PMLs (not shown) everywhere in the transverse directions.
We consider two arbitrary solutions, Φ 1 =|E 1 ,H 1 > and Φ 2 =|E 2 ,H 2 > (represented in Fig. 2) of Maxwell's equations at frequency ω with different dipole sources, J 1 and J 2 , respectively.The sources are located at positions r 1 and r 2 in the arbitrary geometry section.By applying the Lorentz reciprocity theorem of Eq. ( 7) to a volume delimited by the transverse-section plane z=L 1 and z=L 2 , we obtain Assuming again that the QNBMs form a complete set in the periodic-waveguide sections, the fields Φ 1 and Φ 2 can be written as a superposition of backward-and forward-propagating QNBMs in the left periodic-waveguide section.At z=L where the are the complex modal QNBM-expansion coefficients for the solutions 1 and 2 in the left periodic waveguide.Since the form F z is bilinear, the application of the QNBM orthogonality relation straightforwardly leads to Note that for deriving Eq. ( 18), we have used the fact that the QNBMs are normalized (F (p,ω) =4 for any p).With similar notations and for z=L 2 , we have A combination of Eqs. ( 16), ( 18) and ( 19) leads to 4 ( ) By introducing the scattering matrix S that relates the vector D = |D (p,L) ,D (q,R) > (formed by the association of the diffracted modal QNBM-expansion coefficients) to the vector I = |I (p,L) ,I (q,R) > (formed by the incident modal QNBM-expansion coefficients), Eq. ( 20) can be rewritten in a compact form In particular for J 1 =J 2 =0, we have S=S T , like for classical z-invariant waveguides.Note that the symmetry property holds the appropriate QNBMs normalization defined in Section 3. We emphasize that Eq. ( 21) is a direct consequence of the use of materials with symmetrical constitutive parameters.Some of its important implications are summarized in Fig. 3.

Excitation of QNBMs by localized currents: local density of states
Given a source of excitation prescribed by a distribution of currents, either at the periodic waveguide endface or within the waveguide, the modal amplitude of the radiation and bound QNBMs can be straightforwardly calculated.In this subsection, we determine the amplitudes of the bound and leaky QNBMs excited by a dipole source located at r=r 0 outside the PML region.In the electromagnetic description adopted here and in the weak coupling regime, this amounts to determine the classical photonic local density of states (LDOS) [34-35] of periodic waveguides.Referring to Fig. 4, the total electromagnetic field, Φ=|E, H>, radiated by a Dirac dipole source J δ(r-r 0 ) can be expanded in terms of the complete set of outgoing QNBMs and for z<z 0 , where the D (p,R) and D (p,L) are the modal amplitude coefficients of the forward and backward QNBMs, respectively.To calculate these coefficients, we apply Eq. ( 20) for r 1 =r 2 =r 0 .For solution 1, we consider , J 1 =J and r 1 =r 0 , in Eq. ( 20)], while solution 2 is successively prescribed by Φ 2 =Φ (m) [i.e.
These equations show that the total electromagnetic field radiated by the Dirac source is analytically determined by the sole knowledge of the QNBM electric fields on the source.Note that the excitation coefficient D (m,R) of a forward-propagating QNBM Φ (m) is directly proportional to the field amplitude of the backward-propagating QNBM Φ (-m) at the dipole source location, and vice versa.Truly guided Bloch modes deserve a particular attention.By virtue of Eq. ( 14), we have E (-1) (r 0 )=[E (1) (r 0 )] * , since the truly guided QNBMs of the computational system, Φ (1,ω) and Φ (-  1,ω)   , are identical to the real Bloch modes, Φ (1,r) and Φ (-1,r) , inside the inner PML boundaries.For a linearly polarized source, J = Iu, I being a complex number and u a 3x1 real vector.From Eq. ( 23) and from E (-1) (r 0 )=[E (1) (r 0 )] * , it is easily shown that the excitation coefficients D (1,R) and D (1,L) of the forward-and backward-propagating fundamental Bloch modes are related by the simple relation : D (1,R) I * =D (1,L)* I. Remembering that the QNBM are all normalized in this Section, Φ (1,ω) and Φ (-1,ω) have a unitary Poynting flux in the z direction.Thus the powers P 1 or P -1 , coupled into either Φ (1,ω) or Φ (-1,ω) , are equal to |D (1,R) | 2 and |D (1,L) | 2 , respectively.They are found to be equal and are given by Perhaps counter intuitively, Eq. ( 24) shows that the coupled powers into Φ (1,ω) or Φ (-1,ω) are equal, independently of the location of the source in the unit cell of the periodic waveguide, or of the existence of a symmetry axis for the waveguide.Additionally, since the electric fields of a normalized bound QNBM scale as (a/ ) 1 ( g v ) 1/2 , the power emitted diverges as 1/ ) 1 ( g v in the slow light regime.Note that this divergence holds for lossless and infinitely long periodic waveguides, and that any waveguide termination or any residual material absorption would prevent such a divergence.

Numerical implementation
In this Section, we discuss practical issues in relation with QNBM calculations and with their use for the analysis of light scattering in periodic waveguides.Firstly, the main steps to implement the aperiodic Fourier Modal Method (a-FMM) are described.The principles of this method originate from a generalization [12,36] of classical grating methods known as the Rigorous-Coupled-Wave-Analysis [37], to handle non-periodic geometries through an artificial periodization and the use of PML.A particular attention is paid to the influence of the PML on the QNBM calculation and also to the scattering matrix formalism.Then, to illustrate our purpose and to test the accuracy of the a-FMM, two basic 3D scattering problems related to light reflection and emission in periodic waveguides are considered.
In the following, all numerical results are obtained for a PhC silicon (n=3.55)slab waveguide, composed of a line-defect in a 2D triangular lattice of air holes, see Fig. 5(a) for a schematic view of the structure along with definitions of the different parameters.This waveguide supports a single guided TE-like Bloch mode in the gap, as shown by the dispersion curve in the inset.

Quasi-normal Bloch mode calculation
The method used for calculating the QNBMs has been described in a previous work [23] and hereafter, we simply summarize it for the sake of consistency.In general, the computational system associated to the actual periodic waveguide geometry incorporates PMLs in both the x-and y-directions as in Fig. 1(d).Nevertheless, for the PhC-waveguide of Fig. 5(a), we use purely periodic boundary conditions in the y-direction.The validity of this assumption has been checked a posteriori by verifying that the calculated data do not depend on the number P of hole rows surrounding the line defect.All the results provided hereafter are obtained for P = 6.Because of the periodic boundary in the y-direction and of the PML in the x-direction, the field solution of the scattering problem is null on the boundaries of the computational domain.Moreover, the artificial periodization along the x-and y-directions makes the field periodic, which allows to expand the electric E and magnetic H fields in a Fourier basis [12], where G x =2π/Λ x and G y =2π/Λ y , Λ x and Λ y being the lengths of the unit cell of the transverse section.In Eqs.(25a) and (25b), the S αlm and U αlm (α=x, y or z) are unknown z-dependent coefficients.In practice, the Fourier series have to be truncated, we denote by m x and m y the truncation ranks, -m x <p<m x -m y <q<m y .By incorporating expressions (25a) and (25b) into the curl Maxwell's equations, and by expanding the permittivity and the permeability in the Fourier basis, we obtain after elimination of the z-components In Eq. ( 26), Ψ is equal to [S x S y U x U y ], a column-vector formed by the electric-and magnetic-field coefficients and Ω is a matrix formed by the Fourier coefficients of the permittivity and of the permeability.The matrix Ω in Eq. ( 26) depends on the variable z since the permittivity function is z-dependent.The calculation of the QNBMs requires the integration of Eq. ( 26) over a unit cell from z to z + a.For the integration, we approximate the real profile of the circular holes by a stack of slices with locally z-invariant permittivities [23].We have used N s = 9 slices per holes, i.e. a z-discretisation step of ≈ λ/35.Indeed, the accuracy of the computational results increases as N s increases, but calculations performed for N s =19 have revealed that the discretisation error is much smaller than the truncation errors discussed hereafter.Within this approximation, the integration along the z-direction can be performed analytically.The N local modes of each slice (p) can be called quasi-normal modes (QNMs) as they correspond to the mode of a z-invariant waveguide surrounded by PMLs [12].The QNMs, denoted by the vectors (p) n W and (p) n − W (n=1,…N) in the Fourier basis, are calculated in every slice (p) as the eigenvector of a local matrix Ω (p) .Denoting by λ (p) the corresponding eigenvalue (the QNM propagation constant), the electromagnetic fields Ψ (p) can be written as a superposition of QNMs in very slice (p) where b (p) and f (p) are column vectors whose elements are the amplitudes of the QNMs propagating backward (in the negative z-direction) and forward (in the positive z-direction), respectively.The linearity of Maxwell's equations assures the existence of a linear relationship between the mode amplitudes of the slice (i), b (i) and f (i) , in the input z-plane and those of the slice (t), b (t) and f (t) , in the output (z+a)-plane.To avoid any numerical problems, a S-matrix approach is used to relate these amplitudes.In a compact form, we have where the matrix on the right-hand side of the equation is simply the S-matrix of a unit cell.Details concerning the recursive calculation of S can be found in [38].
Since the QNBMs are pseudo-periodic, their mode amplitudes at planes z and z + a are proportional.Denoting by ρ the proportionality factor, we have b (t) =ρb (i) and f (t) =ρf (i) , and incorporating Eq. ( 28), we obtain Eq. ( 29) is solved as a generalized eigenproblem without numerical instabilities [24,39].Note that for periodic waveguides possessing a mirror-symmetry plane in the z-direction, enhanced performance is achieved by a related simplified eigenproblem scheme Figure 5(b) shows the 300-first QNBM propagation constants k p obtained for the geometrical parameters defined in the caption and for m x =15 and m y =26.The calculation has been performed for two different PMLs, implemented as a linear complex coordinate stretching [36].Blue dots and red squares respectively correspond to (f PML ) -1 =(1+i) and (f PML ) -1 =5(1+i), where f PML is the stretching parameter defined in Annex 2. With the exception of the fundamental QNBM located on the real axis near Real(k p )=9, the distribution of k p depends on the specific PML implementation.This is understood by considering that the actual radiation Bloch modes of the periodic waveguide in Fig. 1(a) are modified by the PML and that the modification depends on the PML themselves.In other words, a PML represents a cut in the complex plane that allows to satisfy outgoing-wave conditions.By changing the PML parameters (thickness, stretching, gradual variation …), the cut is modified and consequently the QNBMs distribution varies.Yet, some of the QNBMs weakly depend on the PML parameters, they deserve a special attention.An inspection of their electromagnetic field has revealed that some of them, especially those with small Im(k p ) values (attenuation) and with 10<Re(k p )<18, possess two zeros in the x-direction of the membrane, indicating that they are related to the leaky TE 2 mode of the membrane without holes.More generally, we believe that these modes which are almost insensitive to the PMLs' choice, are leaky QNBMs but a quantitative discussion would deserve further studies.Note that the relation between leakymode expansion and PML with varying absorption for a simple slab waveguide has already been discussed in the literature, see Ref [10].A thorough discussion of the accuracy of the a-FMM method to calculate the attenuation of leaky Bloch modes operating above the cladding light lines can be found in [41].

Light scattering in periodic waveguides
To illustrate the potential of the approach, we consider two basic 3D scattering problems in PhC-waveguide aggregates, the reflection of light onto a semi-infinite PhC mirror [Fig.6(a)] and the emission of a dipole source into a semi-infinite PhC waveguide closed by a mirror at one of its extremity [Fig.7 to be satisfied in the QNBM basis to take into account terminations by semi-infinite periodic waveguides.For that purpose, one needs to derive the scattering matrices which link the forward-and backward-modal coefficients, f QNBM and b QNBM , in the QNBM basis to the forward and backward modal coefficients, f QNM and b QNM , in the QNM basis.In General, two cases have to be considered.We simply consider a termination with a semi-infinite periodic waveguide that extends to the positive z-direction hereafter, but similar considerations can be straightforwardly applied for termination with semi-infinite periodic waveguides extending towards the negative z-direction.For this termination, the S-termination matrix S T reads as .The S T scattering matrix defined in Eq. ( 30) and ( 26) can be used with standard Smatrix products [38] to handle intricate scattering geometries.The S-termination matrices approach has been implemented for solving the scattering problem of Fig. 6(a) by using the products of S T matrices associated to the two different QNBMs encountered for z < 0 (mirror) and z >0 (waveguide).Figure 6(b) shows the convergence performance for the modal reflectivity R of the fundamental PhC bound mode as a function of m x for several values of m y .The calculation is performed for a/λ = 0.255.For the sake of convergence performance, the cladding PML are implemented as complex nonlinear coordinate transforms [36].As expected, both m x and m y impact the computational accuracy.For all curves, we note that a plateau is obtained for m x > 30, the peak-to-peak residualoscillation amplitude being smaller than 0. m x values (m x = 45).On overall, R is estimated to be equal to 0.9740 ± 0.0001.Note that similar calculations performed for a related scattering geometry have been previously reported [42] and that the agreement between the experimental and numerical data has indicated an absolute error below 0.0002 for the calculated modal reflectivity.Typical CPU times for the computation of R in Fig. 6  Let us now consider the emission of a dipole source J δ(r-r 0 ) in the geometry shown in Fig. 7(a).As in the previous example, two periodic sections are involved.They will be labelled by the subscripts "M" and "W" hereafter, "M" referring to the PhC mirror on the left side and "W" to the single-row defect waveguide on the right side of the figure.
In Section 4.2, we have solved the emission of a dipole source in a periodic waveguide by providing an analytic experession for the QNBM excitation coefficients D (m,L) and D (m,R) .In fact, as shown by Eq. ( 20), the presence of a dipole source in a periodic waveguide results in a field discontinuity, which is described in the QNBM basis as a step discontinuity of the QNBM excitation coefficients at z=z 0 .Because of the linearity, the step discontinuity is simply that obtained in Section 4.2 in the absence of incident QNBM illumination, and is equal to the D (m,L) and D (m,R) coefficients given by Eqs.where the Nx1 vectors, D (L) and D (R) , are formed by the D (m,L) and D (m,R) coefficients.Since these vectors are known analytically, Eqs. ( 32) and (33) can be solved for the unknown f W+ , b W-, b M and f W-vectors with b W+ = f M = 0.The electromagnetic-field distribution is then recovered everywhere.Figure 7(b) shows the convergence performance of the spontaneous emission β-factor as a function of m x for different values of m y .The β-factor is defined as the fraction of light emitted into the fundamental forward-propagating QNBM of the waveguide W. It is an important parameters which may drastically impact many optoelectronic devices, like lasers or single-photon sources, if it may be made close to unity [43].Again the data shows a wellconverged situation: a plateau with a negligible residual fluctuation is obtained for m x values as small as 20 and the plateau height is weakly dependent on m y .In addition, cross-checking tests obtained by varying the PML thickness and the number of hole rows surrounding the defect in the y-direction have further confirmed the accuracy of the results.This crosschecking tests are important, because the QNBM basis changes as one varies the PML implementation.But accordingly to Eqs. ( 32) and (33), the excitation coefficients D (m,L) and D (m,R) also change.Obtaining nearly identical β-factor predictions for totally different extension basis thus represents a strong test for the theoretical and computational aspects developed in this work.For this Green-function problem, CPU times are almost identical to those obtained for the previous scattering problem.Other computational data obtained with this approach for the broadband and directive emission of light in single-row-defect PhC waveguides can be found in [44] for the two in-plane dipole orientations.

Conclusion
A rigorous modal formalism for modeling light propagation and light emission in threedimensional periodic waveguides and in aggregates of them has been presented.In essence, this work is a generalization of known modal concepts for translation-invariant waveguides to situations involving aggregates of periodic waveguides.By surrounding the actual stack by PMLs in the transverse directions, we have shown that both radiation and bound modes can be handled and that reciprocity considerations lead to the derivation of Bloch-mode orthogonality relations in the sense of E×H products, to normalization of these modes, and to the proof of the symmetrical property of the scattering matrix linking the Bloch modes.The general formalism which rigorously takes into account radiation losses resulting form the excitation of radiation Bloch modes has been implemented with a Fourier numerical approach.Basic examples of light scattering like the emission of a dipole source in periodic-waveguide aggregates are accurately resolved.
All the above results (orthogonality, reciprocity ...) can be straightforwardly generalized for other systems involving additional periodicities.For instance, for 2D-periodic systems in thin film stack, like semiconductor membranes perforated by triangular hole arrays, the QNBMs possess an in-plane parallel wave-vector and PMLs have to be introduced only above and below the stack.When using the Lorentz reciprocity theorem, the integral over every cross-section planes of the unit cell is not null, but due to the in-plane pseudo-periodicity, the integral over two parallel cross-section planes of a unit cell are identical and their contributions exactly cancel.For 3D periodic systems, similar results hold and there is no need to introduce PMLs.Furthermore, if the media of the fully-periodic system are lossless, the orthogonality relation could be used with E×H * products.
By surrounding the cross-section with PMLs, complex permittivity and permeability are introduced even if the ε and µ of the physical system are real.As a consequence, orthogonality in the sense of the Poynting vector no longer holds and the system energy cannot be reduced to a sum of QNBMs energies even for lossless materials.At first sight, this might appear a drawback of the approach, but we emphasize that inside the PML inner boundaries, since the solution for the total electromagnetic field is virtually exact, the total energy flow on a closed surface can be fully predicted, as shown by the analysis of light emission in Section 5.2.In addition, we note that lossy materials, like metals at optical , where X'=dX/dx, Y'=dY/dy and Z'=dZ/dz.
In general, X', Y' and Z' are piecewise-constant functions with a gradual variation [46].For the numerical results shown in Fig. 5, the X' function extends over a finite thickness d PML and is constant: X'=f PML ; Y' and Z' are both equal to 1.After a few elementary algebraic operations, it is shown that, if E and H satisfy the Maxwell's equations of Eqs.(2a) and (2b) in the (x,y,z) coordinate systems, the new electromagnetic fields H ˆ=L -1 H and E ˆ=L -1 E also satisfy the same Maxwell's equations in the new coordinate system provided that µ ˆ= LµL/Det(L), ε ˆ= LεL/Det(L), J ˆ= LJ, (A2.2) where µ ˆ, ε ˆand J ˆ represent the permeability, permittivity and dipole source in the new coordinate system and Det(L) denotes the determinant of the operator L. With this formalism, it is easily found that important quantities like F z and E V are conserved quantities that are independent of the complex transformation.Let us start with the mode volume ∫∫∫ V (E T εE - #80752 -$15.00USD Received 6 Mar 2007; revised 16 May 2007; accepted 16 May 2007; published 20 Aug 2007 (C) 2007 OSA 3 September 2007 / Vol. 15, No. 18 / OPTICS EXPRESS 11044

Fig. 1 .
Fig. 1.Actual and computational geometries considered in this work.(a) Sketch of a geometry composed of aggregate of different periodic-waveguide sections.The geometry is assumed to be surrounded by infinite uniform claddings in the x-and y-transverse directions.(b) Associated computational system obtained by bounding the actual waveguide with PMLs (in blue) in the transverse directions.(c) A periodic-waveguide section of (a).(d) Associated periodic waveguide bound with PMLs.In (b) and (d), the PMLs have a finite thickness that is not represented for the sake of clarity.

Fig. 3 .
Fig. 3. Some implications of Eq. 21.(a) The complex modal transmission-coefficients do not depend on the propagation sense.(b) Same property for the cross modal reflection-coefficients.(c)The excitation amplitude of a QNBM by a dipole source J δ(r-r 0 ) located at point r = r 0 is equal to the scalar product between the source J and the field E(r 0 ) scattered at the dipole location by exciting the same geometry with the reciprocal QNBM.

Fig. 4 .
Fig. 4. Excitation of QNBMs by a Dirac dipole source J δ(r-r 0 ) located at point r = r 0 .The D (p,R) and D (p,L) coefficients represent the modal amplitude coefficients of the excited forwardand backward-QNBMs, respectively.The periodic waveguide is not necessarily symmetric for the study, as shown by the échelette profile.

Fig. 5 .
Fig. 5. QNBM calculation of a PhC waveguide.(a) Schematic view of the PhC waveguide formed by removing a line defect in the ΓK direction of a 2D PhC structure composed of a triangular lattice of air holes (lattice constant a=0.24 µm) etched into a silicon slab (n=3.55).The slab thickness is 0.6a and the air holes radii 0.29a.The inset shows the dispersion relation of the fundamental guided QNBM Φ (-1) .(b) Display of the 300-first normalized propagation constants of the QNBMs for a frequency a/λ=0.255, point A in the inset.Blue dots and red squares are obtained for (f PML ) -1 =(1+i) and (f PML ) -1 =5(1+i), respectively.
[40].Because of the truncation, a total of 2N=4(2m x +1)(2m y +1) QNBMs are calculated.One half corresponds to forward-propagating modes and the other half to backward-propagating ones.The QNBMs are denoted , whose coefficients represent the components of the m th QNBM in the local QNM basis associated to the slice at planes z and z+a.From the eigenvalues ρ, the propagation constants of the QNBMs are obtained from ρ=exp[ik a], where the vector k is formed by the propagation constants k p of the QNBMs.
(a)].To solve these problems, the outgoing wave conditions have #80752 -$15.00USD Received 6 Mar 2007; revised 16 May 2007; accepted 16 May 2007; published 20 Aug 2007 (C) 2007 OSA 3 September 2007 / Vol. 15, No. 18 / OPTICS EXPRESS 11055 matrix is simply related to a basis transformation.By matching the continuous tangential field components at the termination in the Fourier basis, like for classical zinvariant waveguides, elementary algebraic manipulations lead to S 22 = (F + ) -1 , S 12 = B + S 22 , S 21 = -S 22 F -and S 11 = B --S 12 F -, (31) where B -and B + denote the NxN matrices whose column vectors are the Nx1 eigenvectors, and F + are related matrices formed with the vectors QNBM m − f and QNBM m f

Fig. 6 .
Fig. 6.Scattering at the interface between two periodic sections.(a) Schematic top view of the 3D scattering problem.The PhC parameters are the same as in the caption of Fig. 5.(b) Convergence of the a-FMM for the modal reflectivity R of the fundamental guided QNBM Φ (-1) .The calculation is performed for a/λ = 0.255 , point A in the inset of Fig. 5(a).
0001.The curves obtained for different m y are vertically shifted from each other by an offset value that rapidly decreases as m y increases.This property allows for an accurate interpolation of R. Starting from the calculated data obtained for small m x values (m x = 15) and for large m y values (m y = 30 for instance), the data can be corrected from the deviation observed on the curve with m y = 15 and m x = 15 to large #80752 -$15.00USD Received 6 Mar 2007; revised 16 May 2007; accepted 16 May 2007; published 20 Aug 2007 (C) 2007 OSA 3 September 2007 / Vol. 15, No. 18 / OPTICS EXPRESS 11056 (a) with (m y /m x )=(15/15), (25/20) and (20/35) are approximately 30 min, 300 min and 900 min on a PhC computer equipped with a 3-GHz Intel Pentium 4 processor and with Matlab.

Fig. 7 .
Fig. 7. Dipole emission into PhC waveguides closed at one extremity by a PhC mirror.(a) Schematic top view of the 3D problem.The PhC parameters are the same as in the caption of Fig. 5.The dipole is parallel to the x-axis and is located in the central plane of the membrane at z 0 = 0. (b) Convergence of the a-FMM for the β-factor defined as the power emitted into Φ (1) normalized to the total power emitted.The calculation is performed for a/λ = 0.255 , point A in the inset of Fig. 5(a).
For the calculation, one may first compute the scattering matrix S T which links the forward-and backward-modal coefficients, b M and f M , in the QNBM basis of the PhC mirror to those (denoted b W-and f W-) of the single-row defect waveguide at plane z = − (23a) and (23b) for the QNBMs of the single-row defect waveguide.Thus if we denote by b W+ and f W+ the forward-and backwardmodal of the single-row defect waveguide at plane z = + 0 z and by b W-and f W-the similar quantities defined at plane z = $15.00USD Received 6 Mar 2007; revised 16 May 2007; accepted 16 May 2007; published 20 Aug 2007 (C) 2007 OSA 3 September 2007 / Vol. 15, No. 18 / OPTICS EXPRESS 11057 #80752 -$15.00USD Received 6 Mar 2007; revised 16 May 2007; accepted 16 May 2007; published 20 Aug 2007 (C) 2007 OSA 3 September 2007 / Vol. 15, No. 18 / OPTICS EXPRESS 11058 wavelengths, have by essence complex permittivities.Therefore the formalism developed here fully applies, as shown by the calculation of PhC surface-plasmon-polariton reflectance and transmittance reported in Ref. [45].by L the 3x3 diagonal matrix L =