1. 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 refractive-index 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 z-invariant 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-radiation-leaky 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–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–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–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 non-conservative 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.

2. 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 z-coordinate. 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].

 

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.

Download Full Size | PDF

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–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²=-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 ω₁ and ω₂ for the same permittivity and permeability distributions,

Applying the Green-Ostrogradski formula to the vector E ₂×H ₁ on a closed surface S of volume V enclosing the sources J ₁ and J ₂, 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 ₁ and A ₂ defined by z=z ₁ and z=z ₂>z ₁, 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 ₁ and A ₂. Additionally, noting that E ^T ₂ εE ₁=E ^T ₁εE ₂ and that H ^T ₂µH ₁=H ^T ₁µH ₂ since the materials are reciprocal, we obtain from Eq. (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 Φ₁=|E ₁,H ₁> and by Φ ₂=|E ₂,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 Φ ₁ and Φ ₂ that are respectively anti-symmetric and symmetric. With these notations, Eq. (5) is conveniently rewritten

Importantly, let us note that for ω₁=ω₂ and in the absence of any source (J ₁=J ₂=0), the term on the right side of Eq. (7) is null and thus F_z(Φ ₁, Φ ₂) is independent of z. In this case, it will be simply denoted F(Φ ₁,Φ ₂).

3. 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).

3.1 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 ⁽¹⁾,H ⁽¹⁾>exp(jk₁(ω)z) with Im(k₁)=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 pseudo-periodic 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.

Let us consider two QNBMs, Φ ^(p,ω’) and Φ ^(-q,ω), that satisfy Maxwell’s equations in the absence of source for the same periodic-waveguide geometry, at two frequencies ω’ and ω. Equation (8) implies that F_z+a(Φ ^(p,ω’),Φ ^(-q,ω))=exp{j[k_p(ω’)-k_q(ω)]a} F_z(Φ ^(p,ω’),Φ ^(-q,ω)). The Lorentz reciprocity relation, Eq. (7), applied in the absence of source to a volume V denoted Cell(z) and delimited by the cross-section planes z and z+a, becomes

We now consider two important properties of the QNBMs, which are straightforwardly derived from Eq. (9). Note that similar approaches are classically used for z-invariant waveguides [8,10].

3.2 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.

3.3 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 v⁽¹⁾ _g=dω/dk₁. 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\frac{\partial k_p}{\partial \omega }(\omega -\omega ’)$ we have

where E ^(p,ω)=E_Cell(z)(Φ ^(p,ω), Φ ^(-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, v^(p) _g=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,ω) and Φ ^(-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₁z) and Φ ^(-1,r)=|E ^(-1,r),H ^(-1,r)>exp(-jk₁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 k1 are all real, it is easily shown from the conjugate form of the curl Maxwell’s equations that

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

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/v⁽¹⁾g)^1/2.

4. 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 symmetry-property 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.

4.1 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₁), labelled by “L” for “left”, are denoted by Φ ^(p,L) and those on the right side (z>L₂), 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₂. The whole system is assumed to be surrounded by PMLs (not shown) everywhere in the transverse directions.

Download Full Size | PDF

We consider two arbitrary solutions, Φ ₁=|E ₁,H ₁> and Φ ₂=|E ₂,H ₂> (represented in Fig. 2) of Maxwell’s equations at frequency ω with different dipole sources, J ₁ and J ₂, respectively. The sources are located at positions r ₁ and r ₂ 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₁ and z=L₂, we obtain

Assuming again that the QNBMs form a complete set in the periodic-waveguide sections, the fields Φ ₁ and Φ ₂ can be written as a superposition of backward- and forward-propagating QNBMs in the left periodic-waveguide section. At z=L₁, we have

where the I^(p,L) ₁, D^(p,L) ₁, I^(p,L) ₂ and D^(p,L) ₂ 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₂, we have

A combination of Eqs. (16), (18) and (19) leads to

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 ₁=J ₂=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.

 

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₀) located at point r=r₀ is equal to the scalar product between the source J and the field E(r₀) scattered at the dipole location by exciting the same geometry with the reciprocal QNBM.

Download Full Size | PDF

4.2 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 ₀ 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 ₀) can be expanded in terms of the complete set of outgoing QNBMs

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 ₁=r ₂=r ₀. For solution 1, we consider Φ ₁=Φ [i.e. I^(p,R) ₁=I^(p,L) ₁=0, J ₁=J and r ₁=r ₀, in Eq. (20)], while solution 2 is successively prescribed by Φ ₂=Φ ^(m) [i.e.I^(p,R)=0,I^(p,L) ₂=δ_p,m and J ₂=0, in Eq. (20)] and then by Φ ₂=Φ ^(-m) [i.e. I^(p,R) ₂=δ_p,m,I^(p,L)=0 and J ₂=0, in Eq. (20)]. We obtain

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.

 

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

Download Full Size | PDF

Truly guided Bloch modes deserve a particular attention. By virtue of Eq. (14), we have E ^(-1)(r ₀)=[E ⁽¹⁾(r ₀)]*, 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=I u, I being a complex number and u a 3x1 real vector. From Eq. (23) and from E ^(-1)(r ₀)=[E ⁽¹⁾(r ₀)]*, 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₁ or P-₁, coupled into either Φ ^(1,ω) or Φ ^(-1,ω), are equal to |D^(1,R)|2 and |D^(1,L)|², 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/v⁽¹⁾g)^1/2, the power emitted diverges as 1/v⁽¹⁾ _g 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.

5. 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.

5.1 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 [23]

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.

 

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.

Download Full Size | PDF

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 W(p)n and W ^(p)-_n (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 ⁽ⁱ⁾ and f ⁽ⁱ⁾, 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⁽ⁱ⁾ and f ^(t)=ρf⁽ⁱ⁾, 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 [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 $\left[\begin{array}{c}{\mathbf{b}}_m^{\mathrm{QNBM}}\\ {\mathbf{f}}_m^{\mathrm{QNBM}}\end{array}\right]$, with m=1,…2N. The b^QNBM _m and f^QNBM _m are Nx1 vectors, 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.

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 fPML 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 TE2 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 leaky-mode 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].

5.2 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(a)]. To solve these problems, the outgoing wave conditions have 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 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 z-invariant waveguides, elementary algebraic manipulations lead to

where B- and B ⁺ denote the NxN matrices whose column vectors are the Nx1 eigenvectors, b ^QNBM _-m and b ^QNBM _m and F- and F+ are related matrices formed with the vectors f ^QNBM _-m and f ^QNBM _m. The S _T scattering matrix defined in Eq. (30) and (26) can be used with standard S-matrix products [38] to handle intricate scattering geometries.

 

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

Download Full Size | PDF

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 mx 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 my 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.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 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(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. 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. (b) Convergence of the a-FMM for the β-factor defined as the power emitted into Φ⁽¹⁾ normalized to the total power emitted. The calculation is performed for a/λ=0.255, point A in the inset of Fig. 5(a).

Download Full Size | PDF

Let us now consider the emission of a dipole source J δ(r-r ₀) 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. For the calculation, one may first compute the scattering matrix ST 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=z-₀

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₀. 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. (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=z+₀ and by b ^W- and f ^W- the similar quantities defined at plane z=z-₀, we have

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 well-converged 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 cross-checking 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.

6. Conclusion

A rigorous modal formalism for modeling light propagation and light emission in three-dimensional 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 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].