1. Introduction

Photonic crystals (PC), which are arguably amongst the most exciting of optical structures, are at the forefront of theoretical and experimental research [1]. Their significance is based on properties which resemble those of semiconductors, and which enable them to tailor the propagation of light on the scale of optical wavelengths with minimal losses. Their ability to guide light on the wavelength scale [2–4] and to control the radiation dynamics of embedded sources [5] make them ideally suited to the miniaturization of optical components and to the development of compact devices that incorporate a large number of elements. Already demonstrated are the bending of light around corners [6], interconnections such as T- and Yjunctions [1,7], channel-drop filters [8], superprisms [9] and Mach-Zehnder interferometers [10]. Attention is now being turned to the modeling and fabrication of ultra-compact optical devices (such as the coupler in Fig. 1) which comprise various elements connected through complex waveguide structures embedded in photonic crystals [11].

The exploitation of the technological potential of PC based structures requires a thorough understanding of the mechanisms by which light is coupled into, and guided through, these devices. In order to understand the performance of such structures it is important to develop rigorous analytical techniques that can provide physical insight into the scattering and diffraction processes, facilitating the conceptualization of the device using generalizations of familiar ideas derived from other areas of optics (e.g. thin film optics). Until now, much of the modeling has been undertaken using finite difference time domain (FDTD) methods or other computational schemes. Though these methods produce accurate results, they do not easily reveal insight since they neither exploit the underlying physics to any significant extent, nor take advantage of the structure of the problem or its geometry in order to accelerate calculations. For a device such as that in Fig. 1, comprising a sequence of components and waveguides, there exists the real possibility that semi-analytic techniques, which exploit the geometry and the mode structure of the individual devices, can be superior to exclusively numerical simulations.

This paper outlines an approach based on Bloch mode expansions, which appears to be attracting increasing interest amongst the modeling community [12–15]. For these methods to be rigorous, all Bloch modes, propagating and evanescent, need to be included in the calculations. In practice, of course, the field expansions are truncated, depending on the accuracy that is required. However, often the physics is described by the propagating modes only, and a good asymptotic approximation can then be obtained by dropping all evanescent modes. The resulting techniques are flexible and efficient and have been applied by us to the study of propagation in finite PC waveguides [13], the design of compact resonant filters [11], and apodisation [16]. The computational efficiency of the technique follows from the use of the natural Bloch basis which comes to the fore in modelling propagation in extended regions (i.e., regions with many identical layers) due to the dominant contribution made by the propagating modes (i.e., the modes which can carry energy to infinity). In what follows, we outline the complete rigorous method in Section 2, focusing on Fabry-Perot style (FP) devices, and demonstrating how closely this formulation for PC devices mirrors the Airy formulae for FP interferometers. This also includes the study of important modal conservation properties. Then, in Section 3, we outline the asymptotics of the system, and demonstrate the utility of the approach in generating simple analytic expressions for reflectance and transmittance for a photonic crystal device

2. Theoretical formulation

We consider a two-dimensional (2D) structure, as in Fig. 1, that consists of three segments, M₁, M₂, and M₃, each of which comprises a set of identical layers which may be conceptualized as diffraction gratings with a transverse period D_x . While this derivation is restricted to a three segment device, the approach extends naturally to N -segment structures via recursion. Because of the gratings’ periodicity, the functional elements of the structure are contained within a supercell (see Fig. 1), the dimension D_x of which is chosen large enough to ensure effective isolation from neighbouring supercells when operated within a band gap of the bulk crystal. The periodicity imposed by the diffraction grating model means that individual layers are coupled together by plane wave diffraction orders, the directions of which are given by the grating equation. Between each layer, the field can be represented in the plane wave basis, by an expansion over all plane wave orders. The action of the individual gratings, all of which must have a common transverse period D_x , is handled by plane wave scattering matrices (R and T) that characterize the reflection and transmission of plane waves incident on a grating layer. For example, the element R_pq specifies the reflected amplitude in order p due to unit amplitude incidence in order q. For simplicity here, we assume gratings are up-down symmetric and arranged in a square or rectangular lattice.

 

Fig. 1. A typical three segment photonic crystal device showing the component regions M₁, M₂ and M₃, three lateral supercells of the model and a constituent grating bounded by dashed lines.

Download Full Size | PDF

The modelling of the composite device comprises two distinct elements: (a) Bloch mode propagation in a given region and (b) the scattering/diffraction of modes that occurs at interfaces, and which is modeled with generalized Fresnel (matrix) coefficients that satisfy various energy conservation, reciprocity and transitivity relations.

Each of the structural regions may be homogeneous (e.g., free space or dielectric) or a periodically modulated structure such as a photonic crystal device (e.g., with waveguides as shown), with its Bloch modes generated via a transfer matrix technique [12] that is based on supercell methods and diffraction grating theories [17]. Between the grating layers (e.g., on the dashed lines of Fig. 1) in any given region, the modes are characterized by plane wave expansions represented by vectors of plane wave amplitudes f- and f ₊ (which include entries for both propagating and evanescent waves), respectively denoting the downward and upward propagating plane wave field components. The Bloch modes are represented by the solutions of the eigenvalue problem for the transfer matrix 𝒯 that relates fields on either side of a grating layer according to 𝒻₂ =𝒯𝒻 ₁ and which is derived from relations between outgoing and incoming fields expressed in terms of the reflection and transmission matrices R and (e.g. [12]). The eigenvalue equation to be solved follows from the imposition of the Bloch condition (with Bloch factor µ), i.e., 𝒻₂ µ𝒻 ₁. That is,

where the complex eigenvalues, µ, corresponding to propagating Bloch modes have unit amplitude, and |µ|≠1 for non-propagating modes. The transfer matrix form in Eq. (1) is analogous to transfer matrices used in conventional thin-film optics. We observe that while the eigenvalue problem is stated formally in terms of the transfer matrix, numerical instabilities require the use of alternative methods for solving the eigenvalue problem [18,19]. While the scattering matrices can be computed in a variety of ways (e.g., differential methods, integral methods etc.), we use a multipole analysis [17]. The multipole formulation analytically preserves energy conservation and reciprocity within the scattering matrices [20] and leads to desirable corresponding properties within the Bloch mode analysis.

It is important to ensure that each supercell is sufficiently isolated from its neighbours. To achieve this, an isolating barrier of between 7–10 cylinders of bulk crystal is needed between one defect and its nearest neighbour in an adjacent supercell. In turn, this determines the period D_x of the grating supercells. Once effective isolation is achieved, the calculations are independent of the lateral component of the Bloch vector, k_0x . For analytic convenience in what follows, we set k_ox =0 and thus work with periodic boundary conditions.

We next consider some formal properties of the transfer matrix, and use these to demonstrate reciprocity and energy conservation relations, implicit within the modal formulation, and also establish orthogonality properties of the modes. From reciprocity considerations [22] of any layer or sequence of layers, we can show that the inner (scalar) product 𝓆^T𝓠𝒻, defined in terms of the anti-symmetric matrix Q=-Q^T of Eq. (2), must be independent of layer. Here, 𝒻=[f ^T ₋${\mathbf{f}}_{+}^T$] ^T and 𝓆=[g ^T ₋g ^T₊ ] ^T denote two arbitrary plane wave fields. The anti-symmetric inner product is a generalization of the familiar cross, or outer, product and its independence of layer implies that ${\mathcal{g}}_1^T$ 𝓠 𝒻 ₁=${\mathcal{g}}_2^T$ Q 𝒻 ₂, where the subscripts denote the two sides of the layer through which fields are related by the transfer matrix, e.g. 𝒻 ₂=𝒯𝒻 ₁. In this way, we can show that 𝒯 is symplectic with respect to 𝓠 [21], i.e.

In Eq. (2),Q is the reversing permutation (derived by reversing the rows of the unit matrix) and arises through the reciprocity-based derivation [22]. The symplectic nature of 𝒯 holds even in a system with loss, and ensures that the eigenstates are arranged into forward and backward propagating pairs, respectively associated with eigenvalues µ and 1/µ. For the structure in Fig. 1, in which each layer is up-down symmetric and the lattice is rectangular, the diagonalized form of 𝒯 is then

In the matrix 𝓕 of Eq. (3), the left and right partitions represent the forward and backward propagating modes, with the matrices F _∓ comprising the column vector components of the eigenvectors of 𝒯. In turn, the partitioned matrix contains ℒ the eigenvalues for the forward and backward propagating states respectively in the diagonal matrices Λ and Λ^-1, where Λ=diag{µ_j } and {µ_j } denotes the set of eigenvalues for the forward states. The columns of the matrix 𝓟 are then scaled so that physical properties such as reciprocity and energy conservation, and also modal orthogonality, are represented in a physically normalized form. For completeness, we state these results without proof, referring the reader to Ref. [22] for their derivation. Modal reciprocity, which follows from the symplectic nature of 𝒯 is characterized by

while modal orthogonality is expressed by the relation

arising as a consequence of the flux conservation relation 𝒯^H 𝒯_p 𝒯=𝒯_p that is satisfied by the transfer matrix [22], which in turn follows from the representation 𝒻^H 𝒯_p 𝒻 of the energy flux carried by a plane wave field 𝒻=[f ^T ₋${\mathbf{f}}_{+}^T$] ^T [20]. In Eq. (5), I _r and I _e are diagonal matrices with unit entries on the diagonal respectively for real (propagating) and evanescent plane waves. Correspondingly, I _m and I _m ̄ are diagonal matrices with unit entries on the diagonal respectively for propagating and non-propagating Bloch modes.

With the formal properties of the relevant matrices now stated, we turn now to the propagation of Bloch modes through photonic crystal structures formed from stacks of identical layers, as illustrated in Fig. 1. In the following treatment, we introduce reflection and transmission scattering matrices expressed in the Bloch mode basis of the photonic crystal stack in which we are interested. To differentiate between the Bloch mode and plane wave scattering matrices, from here on we express Bloch mode terms in a sans serif font. The reflection and transmission problem, from medium M ₁ to M ₃, through the L layers of M ₂, is characterized by an incident field δ, a vector of forward propagating Bloch mode amplitudes in M ₁, r=Rδ, a vector of backward propagating modal amplitudes in M ₁, and t=Tδ, a vector of forward propagating modes in M ₃. At the interface between each stack (e.g. the ij interface from M_i to M_j ), the reflection and transmission of Bloch modes can be characterized by Fresnel matrices R _ij and T _ij , which are expressed in terms of the Bloch modes of media M_i and M_j . These matrices are essentially generalizations of planar interface Fresnel coefficients. The Fresnel matrices are derived by considering a Bloch mode expansion represented by the amplitude vector ${\mathbf{c}}_i^{-}$ incident onto the ij interface, giving rise to reflected (${\mathrm{c}}_i^{+}$ =R _ij ${\mathbf{c}}_i^{-}$ ) and transmitted (${\mathbf{c}}_i^{-}$ =T _ij ${\mathbf{c}}_i^{-}$ ) Bloch modes, respectively in M_i and M_j . Field matching at the ij interface then leads to

where R _i =F ⁺_i(${\mathbf{F}}_i^{-}$ )^-1 denotes the plane wave reflection scattering matrix of a semi-infinite region of material i, with incidence from free space. With the modes suitably normalized, these coefficients satisfy reciprocity relations which take the form ${\mathbf{T}}_{\mathit{\text{ij}}}^T$ =T _ji and, and ${\mathbf{R}}_{\mathit{\text{ij}}}^T$ =R _ij energy conservation relations which take the form ${\mathbf{R}}_{\mathit{\text{ij}}}^H$ R _ij +${\mathbf{T}}_{\mathit{\text{ij}}}^H$ T _ij =1—a simplified form with R _ij and T _ij trimmed to contain only the propagating states.

We can now formulate the propagation problem of the structure in Fig. 1, expanding the field in M ₂ in terms of Bloch modes of amplitudes c _- and c ₊, with phase origins respectively located at the upper and lower boundaries of M₂. The mode matching equations are then

Eqs. (7) and (8) respectively derive from field matching at the upper and lower interfaces of M ₂ and express the reflection and transmission of the Bloch modes at these interfaces. The Λ ^L terms in Eqs. (7) and (8) correspond to mode propagation through the layers between the upper and lower interfaces of medium M ₂. Solving Eqs. (7) and (8), we derive the reflection (R) and transmission (T) scattering matrices of Eq. (9) and see that these are generalizations of the Airy formulation for a Fabry-Perot (FP) interferometer [23]

The form of the reflection matrix in Eq. (9) can be simplified further with the introduction of transitivity relations which are implicit in this formulation. These relations evolve from (9) by setting the width of M ₂ to L=0. In doing so, we derive a pair of relations that are analytically satisfied by the Fresnel matrices of Eq. (6). We now consider the special case when the media M ₁ and M ₃ are identical, for which T ₁₁=I and R ₁₁=0, and the corresponding problem of the transition from M ₂ to M ₂ via a L=0 length layer of M ₁. From these we derive the four transitivity relations that are summarized in the notable result

denotes the Bloch mode S-matrix for the M ₁-M ₂ interface. With these relations, which are independent of the energy and reciprocity relations, we arrive at a simplified form for R,

With the modes in their normalized form (i.e. satisfying Eqs. (4) and (5)), the reflected and transmitted fluxes may be computed directly by taking the square magnitude of the relevant elements in the reflection and transmission scattering matrices. In fact, the generalized form of the energy conservation relations [22] are best expressed in terms of the matrix. In the case of the propagating modes, we may trim the Bloch mode scattering matrices so they contain only entries for the propagating states and show, using an approach analogous to that in Ref. [20], that

where I ₁and I ₃ are identity matrices, the dimensions of which are given by the number of propagating states in each of regions M ₁ and M ₃ respectively. When Eq. (12) is expanded, the diagonal partitions contain the energy conservation relations (e.g. ${\mathbf{R}}_{13}^H$ R ₁₃+${\mathbf{T}}_{13}^H$ T ₁₃=l), while the off-diagonal partitions contain a host of phase relationships (e.g. ${\mathbf{R}}_{13}^H$ T ₁₃+${\mathbf{T}}_{13}^H$ R ₁₃=0), the origins of which lie in the application of time reversibility for lossless systems. The generalization of Eq. (12) to include both propagating and evanescent modes is given in Ref. [22].

We stress that the energy conservation relations Eq. (12), the transitivity relations Eq. (10), and the reciprocity relations which correspond to the symmetry of the S matrix all hold analytically within the multipole formulation and are verified numerically to within effectively machine procession (14 significant figures), a consequence of the use of the multipole formulation for the computation of the scattering matrices. This therefore rules out the use of energy conservation and reciprocity as valid physical tests of the accuracy of the Bloch mode formulation in the case of an implementation in which the grating scattering matrices already preserve these properties analytically. In the case of the multipole formulation that we use, such test are merely tests of the efficacy of the coding of the algorithm. We note, however, that for alternative implementations, in which the grating scattering matrices are not imbued with these analytic properties, energy conservation and reciprocity are indeed valid physical tests of the scattering matrix calculation and have flowon effects in the Bloch mode method. In our treatment, however, convergence of the method, which is dependent on the truncation dimensions of plane wave and modal fields (i.e. the number of evanescent terms included), is the only real means of validating results and comparing them against those obtained by entirely different means. We have confirmed the accuracy of this method using results obtained from a recently developed Wannier function method [24], demonstrating agreement of results to better than 1 part in 1000.

Finally, we observe that the forms in Eqs. (9) and (11) are precise generalizations of the Airy formulae for a Fabry-Perot interferometer [23], with the Bloch mode theory having enabled the formulation of the propagation problem in a stratified photonic crystal device to be cast into a form that is structurally identical to that of thin film optics in uniform media.

3. Applications

3.1 Asymptotic expansions and the folded directional coupler

Recently, we reported upon the design of the folded directional coupler (FDC) [11], a novel compact resonant filter which combines the characteristics of both a directional coupler and a Fabry-Perot interferometer, and which is the basis of a high-Q notch rejection filter. This device (Fig. 1) is simply two semi-infinite waveguides with a common coupling region of length L where the guides run parallel to each other, separated by N_c columns of cylinders. As is discussed in Ref. [11], the properties of the FDC are governed by only the propagating modes, since the coupling region is sufficiently long to ensure substantial decay of the evanescent states. The behaviour of the coupler is thus described by its two propagating Bloch modes—the supermodes which, for well separated guides, are well approximated by the symmetric and anti-symmetric superpositions of the modes of a single guide. We now form asymptotic approximations of the reflection and transmission R ₁₃ matrices T ₁₃ and in these limits. First, we introduce notation that will enable us to take projections that extract the relevant rows and columns of matrices. The projection matrices comprise columns of the identity matrix, the number of which is identical to the number of propagating modes in the given region. In the case of the FDC, we have

Here, w ₁ and w ₃, are used to handle projections in regions M ₁ and M ₃ in which there is only a single propagating mode, while w ₂ is used in the coupler (region 2 w M ₂) which has two propagating modes. We then observe that for L sufficiently large, Λ ^L ≈w ₂Λ ^L ̃${\mathbf{w}}_2^T$ 2, where $\tilde{\Lambda }$ is a 2×2 diagonal matrix containing the eigenvalues of the two supermodes of the coupler. Some manipulation of Eq. (9) then reveals

where T̃ _ij =${\mathrm{w}}_j^T$ T _ij w_i and R̃ _ij =$w_i^T$ R _ij w_i , together with an analogous expression for R̃₁₃. The derivation of these follows from the Woodbury formula [25] for the inversion of a rank-p matrix perturbation. In the case of the FDC, both R̃₁₃ and T̃₁₃ are scalars, the value of which can be obtained in a variety of ways including direct computation. In our earlier paper [11], we derived, in a heuristic manner,

where β̄=(β ₁+β ₂)/2,Δβ=(β ₁-β ₂)/2, with β_j denoting the modal propagation constants of the propagating supermodes in M ₂, for which µ_j =exp(iβ_j ). However, these results Eq. (15) may also be derived [22] as a limiting case of the general theory by assigning particular, idealized values to the Fresnel matrices R̃ _ij and T̃ _ij in Eq. (14) and its reflection matrix equivalent.

3.2 Serpentine waveguide

We now build upon the treatment of the previous section to investigate the serpentine waveguide, a coupled waveguide device formed by cascading multiple FDC structures to form the linear periodic array of resonant waveguides shown in Fig. 2. While coupled resonators such as coupled cavity waveguides [26] and coupled ring resonators [27] have been studied previously, this serpentine structure is of interest because of its potential as a slow-wave photonic crystal structure [28].

The unit cell of the serpentine waveguide is a pair of coupled FDC devices joined input-to- output by a single waveguide of length L ₂ as shown in Fig 2(a). Each of the individual components of this structure is essentially the same as those of the isolated FDC device. In the sections where there are two parallel waveguides (length L ₁), a pair of propagating modes is supported, and light is coupled from one guide to the other. The sections of single guide (of length L ₂) act as the input/output guides to the next/previous coupling region. If L ₂ is sufficiently long that evanescent tunnelling between the blocked waveguide ends is negligible, the only connection between the coupling regions is via the single (even-symmetry) mode of the single guide. In this limit, the structures shown in Figs. 2(a) and 2(b) are equivalent. This property means that a single FDC, rather than a complete unit cell, is the minimum segment of the serpentine required to determine the band structure.

 

Fig. 2. Two equivalent serpentine waveguide geometries (assuming no tunneling through the guide ends). Both are characterized by the period D_y and the double and single guide lengths, L ₁ and L ₂.

Download Full Size | PDF

We now proceed to analyze the structure using the monomodal approximations developed in Section 3.1 for the FDC, valid when L ₂ provides a few lattice periods of separation between each cavity region. The complex Fresnel reflection and transmission coefficients, ρ_f and τ_f of Eq. (15), for the single propagating mode through the double guide cavity of the FDC, are calculated using either the full Bloch method or the approximate analytic form given in Eq. (15). Since ρ_f and τ_f are the reflection and transmission amplitudes of the cavity itself, they must be padded to provide for propagation through the input and output guides to give the properties of the minimum serpentine segment. Padding is included as an extra phase term due to propagation through the single guide, exp(iβL ₂), where β is the modal propagation constant in the single mode guide that connects the coupler regions of length. L ₂

Analysis of the periodic structure begins with the transfer matrix 𝒯_s through a complete period of the serpentine

where ρ_s and τ_s are the appropriately padded Fresnel coefficients for a full period, double FDC structure of length D_Y =(L ₁+L ₂). As observed above, the minimum segment of the serpentine is the single FDC, which leads to

is the transfer matrix through a single FDC with input and output guides of length L ₂/2, and the superscripted prime denotes the padded Fresnel coefficients described above.

The dispersion relation for the periodic serpentine structure is then obtained by solving the eigenvalue problem 𝒯_s𝒻=µ𝒻 where the eigenvalues are µ=exp(iqD_y ) and q is the Bloch factor in the longitudinal direction of the serpentine. Since 𝒯_s =${\mathcal{t}}_f^2$, a consequence of the serpentine period being regarded as a pair of FDC devices in series, the eigenvalue equation can be written as 𝒯_f𝒻=±µ ^1/2 𝒻, the solution of which leads to the dispersion relation

where the simplification exploits the energy conservation properties satisfied by the Fresnel coefficients Eq. (12), guaranteeing that |ρ′_f |2+|τ′_f |²=1 and arg (τ′_f )-arg (ρ′ _f )=±π/2.

The right hand side of Eq. (18) can be calculated using the full Bloch mode method, or alternatively, from the analytic expressions Eq. (15) obtained above, with the appropriate padding terms included. The latter method results in an analytic form of the right hand side, which can provide a better understanding of the serpentine waveguide properties. With τ′=τ_f exp(iβL ₂), we can express Eq. (18) in the form

Transmission bands exist for frequencies where the right hand side of Eq. (19), ℜ(1/τ′_f ), has magnitude less than unity. Correspondingly, band gaps, which lie between the transmission bands, occur when |ℜ(1/τ′_f )|>1. There are three distinct types of band gap, each originating through a distinct physical effect. Following Eq. (19), or from the analysis of a single FDC, there are two conditions for which τ_f may vanish. The first occurs when sin(ΔβL ₁)=0, i.e., for cavity lengths L ₁=nL_B /2=nπ(2Δβ), for an odd integer n, where L_B is the beat length of the double waveguide. The second condition occurs when, cos($\overline{\beta }$ L ₁)=0, L ₁=nπ/(2$\overline{\beta }$), also for odd n. When either of these conditions is satisfied, the right hand side of Eq. (19) diverges, and a resonator band gap [27] appears in the serpentine band diagram. We label these two types of resonator band gap RG1 and RG2 respectively. The third class of band gap that appears in the serpentine band diagram is a Bragg gap [27] that occurs as a result of the overall periodicity of the serpentine structure, rather than a single feature of the FDC.

 

Fig. 3. (a) Band diagram for a serpentine tine waveguide with L ₁=L ₂=5d, where q is the Bloch coefficient along the waveguide, and D_y =L ₁+L ₂. (b) Transmission through an FDC with same parameters (dashed), one period of the serpentine waveguide (dotted) and two periods (solid). Note that for frequencies below ωd/(2πc)=0.3064, the double guide cavity only supports a single, odd mode, and thus the analytic result of (19) does not apply.

Download Full Size | PDF

Figure 3(a) shows the band diagram of a serpentine waveguide with L ₁=L ₂=5d, formed in a square symmetric photonic crystal, with cylinders of normalized radius 0.3d and refractive index 3.0, operated in TM polarized light with N_c =1. The FDC structures making up this device have the same parameters as those in Ref. [11], and indeed, the very narrow RG2-type band gap at ωd/(2πc)=0.33107 corresponds to the sharp resonance of the FDC filter. In Fig. 3(b), the calculated intensity transmission is shown for the single FDC (half a period), a single whole period and two whole periods of the serpentine waveguide. The band diagram was calculated using τ_f obtained from the full Bloch mode method and the transmission spectra were also calculated using the full numerical approach. Results for the approximation Eq. (19) do not agree exactly because of the choice of and the phase change on reflection from the guide ends, which is not included in the FDC model. Counting from the top of Fig 3(a), we see that top two resonator gaps are both of type RG2. Given the resonance conditions for RG1 and RG2, and since β̄>Δβ, it is apparent that resonator gaps of type RG1 occur less frequently than RG2-type gaps.

Another feature of both types of resonator band gap is the positioning of the bands above and below the gaps. Recall that in a one-dimensional Kronig-Penney model [28], the right hand side (RHS) of the dispersion relation is continuous, and hence only direct gaps exist, the RHS of the relation being unable to change sign without passing through zero. In the case of the resonator gaps however, the RHS of the dispersion relation (19) diverges when either of the resonance conditions is satisfied. When this happens, the RHS can change sign without passing back through zero, thus resulting in an indirect band gap. Although band gaps with this property occur commonly in two-dimensional photonic crystals, such features have only recently been observed in other coupled resonator waveguides [27] consisting of a linear periodic array of resonant elements.

A number of interesting band structures can be created with appropriate choices of the lengths and and waveguide separation. In Ref. [11] it was shown that the width of a resonance of type RG2 is largely a function of Δβ, being only weakly dependent on β̄. As can be seen in Fig. 3(a), very flat bands with low group velocities occur for sharp resonances, and hence there is the potential to tune the group velocity and group delay with Δβ, by changing either the guide separation or the properties of the cylinders between the guides. Another band structure of interest occurs when two or more consecutive band gaps are resonator gaps.

 

Fig. 4. (a) Band diagram for a serpentine waveguide with L ₁=L ₂=7d. The solid curve is calculated with the full numerical simulation while the dashed curve is calculated using the approximation (19) L ₁=7.5d, L ₂=6.7d (b) Transmission through a FDC with L=7d (dashed), 2 periods of the serpentine guide (dotted) and 3 periods of the serpentine guide.

Download Full Size | PDF

In Fig. 4(a), the band structure is shown for a serpentine waveguide with L ₁=L ₂=7d. The constituent waveguides are identical to those used previously. The band gaps at centre frequencies ωd/(2πc)=0.3135 and ωd/(2πc)=0.3183 are resonator gaps of type RG1 and RG2 respectively. Since these are consecutive gaps, and each is indirect, the three bands above, below and between the gaps are almost parallel to one other and all have a positive slope. With careful optimization, it may be possible to design a band structure with several equally-spaced bands. Fig. 4(b) shows the transmission spectrum of a single FDC with L=7d and a 2-period long serpentine waveguide with L ₁=L ₂=7d. Observe that, for these parameters, the resonator gaps are relatively strong even for such a small number of periods, whereas the Bragg gaps only appear as very weak perturbations in the transmission spectrum for two periods.

4. Conclusions

This paper has provided an overview to the development and application of Bloch mode methods for modelling extended photonic crystal devices. The method reveals itself to be intuitive, analytically tractable and computationally easy. Indeed, our implementation is built in a combination of both Mathematica [29] and Fortran—with the applications suite being implemented in the former, exploiting the convenient programming language and numerical linear algebra library, and the grating scattering matrices implemented in the latter, with the two linked together using the MathLink toolkit. The method demonstrates computational advantages when handling extended structures which have components with many identical gratings, thus allowing the propagation problem in a long segment to be handled by a single set of modes. However, there is no performance advantage when dealing with varying structures such as tapers, in which each layer can differ from its neighbours.

The choice of the diffraction grating paradigm naturally imposes a supercell on the structure, thus modelling the behaviour of not a single device, but of an infinite parallel array of devices. The utility of the method is thus limited to gap frequencies, which, of course, is where photonic crystals are at their most useful. That aside, the key to the successful implementation of the method is a computationally accurate and efficient method of calculating the grating scattering matrices. While these can be generated in a variety of ways, we have found the multipole method to be particularly effective as it analytically encapsulates within the implementation key physical properties, including reciprocity and energy conservation. All of these properties are inherited by the modal analysis, in turn making the method delightfully tractable from an analytic viewpoint and very easy to validate. The essential highlight of the method, however, is its ability to provide physical and theoretical insight into coupling and propagation problems for photonic crystal devices. The method achieves this through the use of the natural (Bloch) modal bases of the system and generates a solution in a form which maps naturally onto concepts of conventional optics—in this case, thin film optics in uniform media. This suggests that it may be possible to import aspects of the elegant theory and successful designs of thin film optics into the new context of photonic crystals—an exciting possibility which merits investigation.