Photonic crystal devices modelled as grating stacks : matrix generalizations of thin film optics

A rigorous semi-analytic approach to the modelling of coupling, guiding and propagation in complex microstructures embedded in twodimensional photonic crystals is presented. The method, which is based on Bloch mode expansions and generalized Fresnel coefficients, is shown to be able to treat photonic crystal devices in ways which are analogous to those used in thin film optics with uniform media. Asymptotic methods are developed and exemplified through the study of a serpentine waveguide, a potential slow wave device. ©2004 Optical Society of America OCIS codes: (050.1950) Diffraction Theory; (290.4210) Multiple Scattering References and links 1. 2.


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][3][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][13][14][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

Theoretical formulation
We consider a two-dimensional (2D) structure, as in Fig. 1, that consists of three segments, M 1 , M 2 , and M 3 , each of which comprises a set of identical layers which may be conceptualized as diffraction gratings with a transverse period x D .While this derivation is restricted to a three segment device, the approach extends naturally to -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 N x D 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 x D , is handled by plane wave scattering matrices ( and T ) that characterize the reflection and transmission of plane waves incident on a grating layer.For example, the element R pq R specifies the reflected amplitude in order due to unit amplitude incidence in order .For simplicity here, we assume gratings are up-down symmetric and arranged in a square or rectangular lattice.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 T 1 = 2 f T f 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 T µ ), i.e., 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 x D of the grating supercells.Once effective isolation is achieved, the calculations are independent of the lateral component of the Bloch vector, 0 x k .For analytic convenience in what follows, we set and thus work with periodic boundary conditions.0 0 x k = 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 g Q f , defined in terms of the anti-symmetric matrix of Eq. (2), must be independent of layer.Here, and = g g g 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 1 1 2 where the subscripts denote the two sides of the layer through which fields are related by the transfer matrix, e.g.

=
1 f T f .In this way, we can show that is symplectic with respect to [21], i.e.
T Q where = .In Eq. ( 2), 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 Q T µ 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 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 , where 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

I
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 1 M to 3 M , through the layers of L ) and transmitted ( T c ) Bloch modes, respectively in i M and j M .Field matching at the ij interface then leads to where denotes the plane wave reflection scattering matrix of a semi-infinite region of material , with incidence from free space.With the modes suitably normalized, these coefficients satisfy reciprocity relations which take the form and , and conservation relations which take the form ( ) , 7) and ( 8) respectively derive from field matching at the upper and lower interfaces of , .
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 2 M to 0 L = .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 , .
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] , where , , where and are identity matrices, the dimensions of which are given by the number of propagating states in each of regions M respectively.When Eq. ( 12) is expanded, the diagonal partitions contain the energy conservation relations (e.g. ), while the off-diagonal partitions contain a host of phase relationships (e.g. ), 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].

T T I
13 31 13 31 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 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.

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

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 where the guides run parallel to each other, separated by 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 matrices 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 propagating modes.We then observe that for sufficiently large, , where is a diagonal matrix containing the eigenvalues of the two supermodes of the coupler.Some manipulation of Eq. ( 9) then reveals where and , together with an analogous expression for .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 and 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 , ( ) , with j β denoting the modal propagation constants of the propagating supermodes in 2 M , for which exp( ) 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 and in Eq. ( 14) and its reflection matrix equivalent.~ij R ~ij T

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 inputto-output by a single waveguide of length 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 ), a pair of propagating modes is supported, and light is coupled from one guide to the other.The sections of single guide (of length ) act as the input/output guides to the next/previous coupling region.
2 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.L We now proceed to analyze the structure using the monomodal approximations developed in Section 3.1 for the FDC, valid when provides a few lattice periods of separation between each cavity region.The complex Fresnel reflection and transmission coefficients, 2 L 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, where s ρ and s τ are the appropriately padded Fresnel coefficients for a full period, double FDC structure of length .As observed above, the minimum segment of the serpentine is the single FDC, which leads to 1 2

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

/ L
The dispersion relation for the periodic serpentine structure is then obtained by solving the eigenvalue problem s µ = T f f where the eigenvalues are exp( ) and is the Bloch factor in the longitudinal direction of the serpentine.Since T T , a consequence of the serpentine period being regarded as a pair of FDC devices in series, the eigenvalue equation can be written as , 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 ' 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 , we can express Eq. ( 18) in the form ' 2 exp( ) Transmission bands exist for frequencies where the right hand side of Eq. ( 19), ( ) has magnitude less than unity.Correspondingly, band gaps, which lie between the transmission bands, occur when ( ) 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 , also for odd .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., the double guide cavity only supports a single, odd mode, and thus the analytic result of ( 19) does not apply.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.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.Another band structure of interest occurs when two or more consecutive band gaps are resonator gaps.  .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.

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.

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

(
C) 2004 OSA 19 April 2004 / Vol. 12, No. 8 / OPTICS EXPRESS 1595 that is satisfied by the transfer matrix[22], which in turn follows from the representation H p f I f of the energy flux carried by a plane wave field[20].In Eq. (5), and are diagonal matrices with unit entries on the diagonal respectively for real (propagating) and evanescent plane waves.Correspondingly, and are diagonal matrices with unit entries on the diagonal respectively for propagating and non-propagating Bloch modes.
an incident field δ , a vector of forward propagating Bloch mode amplitudes in 1 M , , a vector of backward propagating modal amplitudes in = the reflection and transmission of Bloch modes can be characterized by Fresnel matrices and , which are expressed in terms of the Bloch ij 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 i − c incident onto the ij interface, giving rise to reflected ( T T -a simplified form with and trimmed to contain only the propagating states.ij R ij T We can now formulate the propagation problem of the structure in Fig. 1, expanding the field in 2 M in terms of Bloch modes of amplitudes − c and + c , with phase origins respectively located at the upper and lower boundaries of M 2 .The mode matching equations are then 12 21

2M
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 L 2 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] these relations, which are independent of the energy and reciprocity relations, we arrive at a simplified form for R , 2004 OSA  19 April 2004 / Vol. 12, No. 8 / OPTICS EXPRESS 1598

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

LT
the modal propagation constant in the single mode guide that connects the coupler regions of length . 2 Analysis of the periodic structure begins with the transfer matrix s 2004 OSA 19 April 2004 / Vol. 12, No. 8 / OPTICS EXPRESS 1600 #3950 -$15.00US Received 1 March 2004; revised 5 April 2004; accepted 5 April 2004 Fig. 3. (a) Band diagram for a serpentine waveguide with 1 2 5 L L d = =, where is the Bloch

Figure 3
Figure 3(a) shows the band diagram of a serpentine waveguide with 1 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 L is apparent that resonator gaps of type RG1 occur less frequently than RG2-type gaps.
[11]  it was shown that the width of a resonance of type RG2 is largely a 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.

Fig
Fig. 4. (a) Band diagram for a serpentine waveguide with 1 2 7 L L = = d.The solid curve is calculated with the full numerical simulation while the dashed curve is calculated using the approximation (19) 1 2 7.5 , 6.7 L d L d = = (b) Transmission through a FDC with 7 L d = (dashed), 2 periods of the serpentine guide (dotted) and 3 periods of the serpentine guide.In Fig.4(a), the band structure is shown for a serpentine waveguide with .The constituent waveguides are identical to those used previously.The band gaps at centre frequencies1 2 7 L L d = =