当前位置： >>
Efficient design and optimization of photonic crystal waveguides and couplers The Interface Diffract
Efcient design and optimization of photonic crystal waveguides and couplers: The Interface Diffraction MethodAlexander A. Green, Emanuel Istrate, and Edward H. SargentDepartment of Electrical and Computer Engineering, University of Toronto, Toronto, ON M5S 3G4, Canada aa.green@utoronto.caAbstract: We present the interface diffraction method (IDM), an efcient technique to determine the response of planar photonic crystal waveguides and couplers containing arbitrary defects. Field proles in separate regions of a structure are represented through two contrasting approaches: the plane wave expansion method in the cladding and a scattering matrix method in the core. These results are combined through boundary conditions at the interface between regions to model fully a device. In the IDM, the relevant interface properties of individual device elements can be obtained from unit cell computations, stored, and later combined with other elements as needed, resulting in calculations that are over an order of magnitude faster than supercell simulation techniques. Dispersion relations for photonic crystal waveguides obtained through the IDM agree with the conventional plane wave expansion method to within 2.2% of the stopband width. 2005 Optical Society of AmericaOCIS codes: (230.3990) Mi (230.7370) Waveguides.References and links1. K. M. Ho, C. T. Chan, and C. M. Soukoulis, “Existence of a photonic gap in periodic dielectric structures,” Phys. Rev. Lett. 65,
(1990). 2. K. Busch and S. John, “Photonic band gap formation in certain self-organizing systems,” Phys. Rev. E 58,
(1998). 3. S. G. Johnson and J. D. Joannopoulos, “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Express 8, 173C190 (2001). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173. 4. S. G. Johnson, S. Fan, P. R. Villeneuve, J. D. Joannopoulos, and L. A. Kolodziejski, “Guided modes in photonic crystal slabs,” Phys. Rev. B 60,
(1999). 5. K. S. Yee, “Numerical solution of initial boundary value problems involving maxwells equations in isotropic media,” IEEE Trans. Antennas Propag. 14, 302C307 (1966). 6. A. Taove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, Norwood, MA, 2000). 7. G. Bastard, “Superlattice band structure in the envelope-function approximation,” Phys. Rev. B 24,
(1981). 8. C. M. de Sterke and J. E. Sipe, “Envelope-function approach for the electrodynamics of nonlinear periodic structures,” Phys. Rev. A 38,
(1988). 9. E. Istrate, M. Charbonneau-Lefort, and E. H. Sargent, “Theory of photonic crystal heterostructures,” Phys. Rev. B 66(075121), 075,121 (2002). 10. N. A. R. Bhat and J. E. Sipe, “Optical pulse propagation in nonlinear photonic crystals,” Phys. Rev. E 64(056604), 056,604C1C16 (2001).#8231 - $15.00 USDReceived 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 730411. J. Poon, E. Istrate, M. Allard, and E. H. Sargent, “Multiple-scales analysis of photonic crystal waveguides,” IEEE J. Quantum Electron. 39(6), 778C786 (2003). 12. J. P. Albert, D. Cassagne, and D. Bertho, “Generalized Wannier function for photonic crystals,” Phys. Rev. B 61,
(2000). 13. O. Painter, K. Srinivasan, and P. E. Barclay, “Wannier-like equation for the resonant cavity modes of locally perturbed photonic crystals,” Phys. Rev. B 68(035214), 035,214 (2003). 14. K. Busch, S. F. Mingaleev, A. Garcia-Martin, M. Schillinger, and D. Hermann, “The Wannier function approach to photonic crystal circuits,” J. Phys.: Condens. Matter 15, R (2003). 15. K. Sakoda, Optical Properties of Photonic Crystals, chap. 4 (Springer, 2001). 16. E. Istrate, A. A. Green, and E. H. Sargent, “Behavior of light at photonic crystal interfaces,” Phys. Rev. B 71(195122), 195,122C1C7 (2005). 17. E. Istrate and E. H. Sargent, “Photonic Crystal Waveguide Analysis Using Interface Boundary Conditions,” IEEE J. Quantum Electron. 41, 461C467 (2005). 18. T. Baba, N. Fukaya, and J. Yonekura, “Observation of light propagation in photonic crystal waveguides with bends,” Electron. Lett. 35, 654C655 (1999). 19. S. McNab, N. Moll, and Y. Vlasov, “Ultra-low loss photonic integrated circuit with membrane-type photonic crystal waveguides,” Opt. Express 11,
(2003). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-22-2927. 20. A. Adibi, R. K. Lee, Y. Xu, A. Yariv, and A. Scherer, “Design of photonic crystal optical waveguides with singlemode propagation in the photonic band gap,” Electron. Lett. 36,
(2000). 21. A. Chutinan and S. Noda, “Waveguides and waveguide bends in two-dimensional photonic crystal slabs,” Phys. Rev. B 62,
(2000). 22. K. Yamada, H. Morita, A. Shinya, and M. Notomi, “Improved line-defect structures for photonic-crystal waveguides with high group velocities,” Opt. Commun. 198, 395C402 (2001). 23. Z.-Y. Li and K.-M. Ho, “Light propagation in semi-innite photonic crystals and related waveguide structures,” Phys. Rev. B 68, 155,101C1 (2003). 24. Z.-Y. Li and L.-L. Lin, “Photonic band structures solved by a plane-wave-based transfer-matrix method,” Phys. Rev. E 67, 046,607 (2003). 25. J. B. Pendry, “Photonic band structures,” J. Mod. Opt. 41, 209C229 (1994). 26. L. Lifeng, “Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings,” J. Opt. Soc. Am. A 13,
(1996). 27. S. Noda, A. Chutinan, and M. Imada, “Trapping and emission of photons by a single defect in a photonic bandgap structure,” Nature (London) 407(6804), 608C610 (2000). 28. Y. Akahane, T. Asano, B. S. Song, and S. Noda, “High-Q photonic nanocavity in a two-dimensional photonic crystal,” Nature (London) 425(6961), 944C947 (2003). 29. B. S. Song, S. Noda, T. Asano, and Y. Akahane, “Ultra-high-Q photonic double-heterostructure nanocavity,” Nat. Mater. 4, 207C210 (2005). 30. R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK Users’ Guide: Solution of Large Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods (SIAM, Philadelphia, 1998). http://www.caam.rice.edu/software/ARPACK/. 31. E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. D. Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen, LAPACK Users’ Guide, 3rd ed. (SIAM, Philadelphia, 1999). http://www.netlib.org/lapack/lug/index.html. 32. S. Boscolo, M. Midrio, and C. G. Someda, “Coupling and decoupling of electromagnetic waves in parallel 2-D photonic crystal waveguides,” IEEE J. Quantum Electron. 38, 47C53 (2002). 33. D. M. Whittaker and I. S. Culshaw, “Scattering-matrix treatment of patterned multilayer photonic structures,” Phys. Rev. B 60,
(1999). 34. A. R. Cowan, P. Paddon, V. Pacradouni, and J. F. Young, “Resonant scattering and mode coupling in twodimensional textured planar waveguides,” J. Opt. Soc. Am. A 18,
(2001). 35. L. C. Andreani and M. Agio, “Photonic bands and gap maps in a photonic crystal slab,” IEEE J. Quantum Electron. 38, 891C898 (2002). 36. S. F. Mingaleev and K. Busch, “Scattering matrix approach to large-scale photonic crystal circuits,” Opt. Lett. 28, 619C621 (2003).1.IntroductionPhotonic crystals provide a exible platform for the realization of many optical components including passive, active, and nonlinear devices. Waveguides employing photonic crystals have been studied widely because they conne light strongly and guide light over sharp bends with low losses. The complexity of the interaction between electromagnetic waves and these periodic#8231 - $15.00 USD Received 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 7305structures, however, poses a large challenge for the intuitive, yet quantitative, understanding of their optical response. Because they are periodic, the behaviour of photonic crystals is completely and efciently represented using a reduced-zone representation dispersion diagram [1, 2, 3] and its accompanying Bloch modes. These are calculated from the representation of the crystal in reciprocal space, working with the Fourier series of the dielectric constant. This Fourier space representation, however, prevents the direct use of band structures and Bloch modes in the analysis and design of practical photonic crystal devices, which must necessarily be of nite size, and often incorporate deviations from periodicity. For waveguides, this problem has been overcome by the supercell method [4], where the waveguide including core and claddings is repeated periodically in order to allow a plane wave treatment. The periodic units, however, must be large enough to eliminate coupling between the parallel guides. Furthermore, the resulting modes must be inspected carefully in order to discard the unphysical modes with energy concentrated outside the core. The limitations of periodic representations are surmounted by the nite-difference timedomain (FDTD) simulation method [5, 6], one of the most common tools for the design of photonic crystal devices. FDTD simulations operate in real (direct) space and make few assumptions, if any, about the periodicity of the photonic crystal, enabling them to be applied to almost any structure. At the same time, however, these simulations are inefcient since they cannot apply results from one crystal unit cell to identical neighbouring ones. Between the extremes of complete periodicity and aperiodicity exist a variety of methods that take advantage of the results obtained efciently from the periodic parts of a device while allowing some deviation from this periodicity. Envelope approximations, similar in concept to their use with semiconductor heterostructures [7], have been used to calculate pulse propagation in a nonlinear crystal [8] and for photonic crystal heterostructures [9]. Here the periodicity of each crystal section is used to extract a set of parameters describing the crystal in the same way as effective masses are used in semiconductors. Using these parameters, an envelope equation can be written that does not include the periodicities explicitly and is therefore easy to solve. A similar approach was taken using multiple-scales techniques [10, 11]. Point and line defects have also been represented using Wannier functions forming a localized basis, similar to the tight-binding formalism [12]. Wannier functions have also been used to calculate the resonant states of graded resonant cavities [13] and to derive a set of optimally adapted functions for the simulation of waveguides and cavities [14]. The resonance of light in photonic crystals of nite size has been computed with a plane wave expansion over the entire slab [15]. All these methods replace Maxwell’s equations with a simplied set of equations to be solved in the partially periodic structure. We have recently introduced a method that uses the Bloch modes of innite photonic crystals to calculate the reection, transmission, and diffraction of light at photonic crystal interfaces [16]. These coefcients, similar to the Fresnel coefcients for dielectric interfaces, can then be used to model photonic crystal devices that include interfaces between photonic crystals and homogeneous materials as a succession of effective materials, with propagation inside each material described by its respective band structure. This method has been shown to simulate efciently the response of point defect cavities, as well as line defect waveguides and waveguide couplers [17]. So far this method has been limited by the fact that the devices could only contain large periodic sections, where the electromagnetic eld proles took the form of the bulk crystal modes, and homogeneous materials, with plane wave propagation. Most photonic crystal waveguides and defect cavities demonstrated so far are obtained by removing one cylinder, or a row of cylinders from a periodic 2D crystal [18, 19]. In the process of optimizing the properties of these waveguides and resonant cavities, it was found that sig-#8231 - $15.00 USDReceived 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 7306nicant advantages can be obtained by the introduction of more elaborate defects. The rows of cylinders adjacent to the missing rows can, for example, be modied to affect the propagation in the core [20]. Alternatively, the cylinders in the core do not have to be removed completely, as long as their diameter is changed [21, 22]. Transfer and scattering matrices have been used to model the transmission and reection from photonic crystals that are periodic in two directions but not necessarily in the third. The same matrices have been used to compute the interaction of light with semi-innite photonic crystals [23]. In this paper we demonstrate that simulations conducted over photonic crystal unit cells alone are sufcient to determine the response of photonic crystal waveguides and couplers with complex defect regions. Using our technique, which we refer to as the interface diffraction method (IDM), individual device elements, such as a row of cylindrical defects of a particular radius or a bulk crystal region, can be simulated independently and combined as needed to model a device. In contrast, the supercell waveguide simulation techniques described earlier require computational domains consisting of many unit cells and must be done on a case-bycase basis for even small changes in structure. The IDM also makes few assumptions about the crystal geometry and defect type, making it more general than other approaches. It can be applied to narrow or wide defect regions with abrupt or graded changes in structure, all within the same theoretical framework without large increases in computation time. Moreover, the IDM achieves its computational benets through conceptually simple ingredients C the plane wave expansion and scattering matrix methods C that can be obtained through a variety of different techniques. With the IDM, we take advantage of the respective strengths of reciprocal and real space methods to model different waveguide regions. For the periodic parts of the device we use Bloch modes obtained from a plane wave expansion [16]. The parts of the structure with deviations from periodicity are simulated using a plane wave based scattering matrix approach described in Ref. [24]. The photonic crystal Fresnel coefcients are used to link the two simulation methods allowing very efcient modelling of structures with arbitrary geometries. In this paper, we demonstrate the exibility of the IDM with several types of photonic crystal waveguides and a novel coupler design, at each step verifying its accuracy with comparisons to numerical simulations. 2. 2.1. Theory Bloch mode and scattering matrix interfaceOur method of simulating waveguides divides the devices into periodic cladding and aperiodic core regions connected through innitesimally thin homogeneous material layers. Inside the photonic crystal claddings, elds are described by a superposition of Bloch modes that excite a series of diffracted plane waves inside the core-cladding interface layer. These plane waves propagate into the core and produce a set of reected and transmitted waves whose amplitudes are related through a scattering matrix. In the following discussion, we employ a planar photonic crystal waveguide oriented as illustrated in Fig. 1, with the dielectric slab parallel to the xy-plane and the waveguide extending along the x-direction. The vectors x, y, and z are the Cartesian unit vectors. The Bloch modes of
the photonic crystal cladding are calculated using the plane wave expansion method introduced in Ref. [16]. In this method, we specify the angular frequency ω , the lateral wave vector k0, , and the normal vector n to a given interface to obtain a set of Bloch modes with these prop erties and their corresponding complex wave vector components along n. For a core-cladding
interface parallel to the xz-plane, k0, = k0,x x + k0,z z, n = y, and the resulting complex wave vectors are kB± = k0, + kB± y, where ± indicates the direction of propagation or decay of the j j,y #8231 - $15.00 USD Received 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 7307Bloch mode along y, j labels each of the modes at a given frequency, and kB± are the complex
j,y wave vector components. Thus crystal modes calculated from this method provide a complex band structure composed of Bloch wave vectors with real and imaginary parts. Inside the crystal passbands, the band structure has some modes with purely real wave vectors corresponding to freely propagating Bloch modes and others with complex wave vector components characterizing the exponential decay rates of evanescent Bloch modes. Inside the photonic crystal stopbands, the band structure consists entirely of evanescent modes. We represent the elds inside the cladding with a superposition of these Bloch modes given by the following equation: EB (r) = ∑ ∑ ξ jBαj αl,m,n∑ EBα exp jlmni kBα + Glmn
r j,(1)where α is either + Glmn are the reciprocal lattice vectors specied through indices l, m, EBα are the Fourier components of the B and ξ jBα are the amplitude jlmn coefcients of the Bloch modes.zyCladding
EB,LCore + EL
ELxCladding + EB,R
ER + ERnhnhInterface layersFig. 1. Schematic illustration of the IDM viewing a planar triangular photonic crystal device from above. The core and cladding areas are separated by thin interface layers. The dashed rectangles mark the simulation cells used for each device region.To model two-dimensionally patterned slab waveguides in our 3D plane wave expansion, we simulate a 3D unit cell with articial periodicity normal to the slab plane [4]. Despite the unphysical nature of this approach, extending the cell in the vertical direction can effectively eliminate the interaction of conned modes with neighbouring slabs because these modes decay exponentially above and below the slab. Although we do form a supercell in the vertical direction for these simulation cells, they are still much smaller than the supercells used in other methods, which are extended both laterally and vertically. For a planar photonic crystal with a square lattice of holes, we employ an extended unit cell of height h in the z-direction and of length a in the x- and y-directions. This crystal has a set of reciprocal lattice vectors
given by: Glmn = lb2 + mb1 + nb3 where b1 = 2π a1 x, b2 = 2π a1 y, and b3 = 2π h1 z. At the core-cladding interface, the Bloch modes excite a set of diffracted waves in the homogeneous interface layer that are periodic over the interface plane: Eh (r) = ∑ E+ exp (ikmn,y y) + E exp (ikmn,y y) exp i k0, + Gmn, mn mnm,nr ,(2)#8231 - $15.00 USDReceived 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 7308where Gmn, = (Glmn
z) = 2π ma1 x+2π nh1 z is the set of reciprocal lattice vec
tors projected onto the interface plane, r = xx + z, and kmn,y = [(nh ω /c)2
|k0, + Gmn, |2 ]1/2
z in a homogeneous layer of index nh . The application of boundary conditions over the interface ensures that coupling occurs only between Bloch modes and plane waves with the same propagation constant in the core-cladding interface plane and yields the following set of equations for each diffraction order:∑jB+ jmn,x B+ jmn,z C+ jmn,x C+ jmn,zB jmn,x B jmn,z C jmn,x C jmn,zξ jB+ ξ jB ξ jB+ ξ jB= =+ Emn,x E + mn,x , +
Emn,z Emn,z(3)+ Emn,x E
Emn,z Emn,z∑j1 kmn,ykx kz kz2
kx2 kx kz,(4)where kx = k0,x +Glmn
x = k0,x +2π ma1 and kz = k0,z +Glmn
z = k0,z +2π nh1 . The elements
and C± on the left side of Eqs. (3) and (4) reect the electric and magnetic eld B± jmn,x/z jmn,x/z amplitudes of the crystal modes at the plane y = y0 in the unit cell and are dened: B± = ∑ EB± exp iGl,⊥ y0 , jmn jlmnl l(5) (6)C± = ∑(kB± + Glmn ) × EB± exp iGl,⊥ y0 , jmn j jlmnwhere Gl,⊥ = Glmn
y = 2π la1 is the set of reciprocal lattice vectors projected onto n.
B B+ With column vectors E+ = (
, ξ jB+ , ξ j+1 ,
)T , E = (
, ξ jB , ξ j+1 ,
)T , E+ = B B + +
, Emn,x , Emn,z ,
)T , and E = (
, Emn,x , Emn,z ,
)T , we can combine Eqs. (3) and (4) to form the transfer matrix T containing the diffraction coefcients required to compute the set of plane waves excited by an arbitrary combination of Bloch modes in the cladding: E+ E =T E+ B E B = T11 T21 T12 T22 E+ B . E B (7)The scattering matrix for the interface, which describes the outgoing modes produced by an arbitrary set of incoming ones, can be readily derived from the elements of the transfer matrix [25, 26]. After calculating the transfer matrix describing the waveguide at the cladding to homogeneous material interface, we can simulate the remaining core region of the device using the plane wave based scattering matrix method described in Ref. [24]. In this approach, a dielectric structure with transverse periodicity is divided along the propagation direction into slices separated by innitesimally thin homogeneous lms. For the dielectric slices, we calculate the scattering matrices relating the diffracted plane waves excited in the two lms surrounding each slice. Applying the usual scattering matrix recursion formulae [25, 26] to each of the matrices yields the scattering matrix for the entire structure. The scattering matrix S for the waveguide determines the set of outgoing plane waves E and E+ from the core produced by an arbitrary L R set of waves E+ and E incident on the core: R L E+ R E L =S E+ L E R = S11 S21 S12 S22 E+ L . E R (8)We can now determine the properties of the entire waveguide structure using the S-matrices and T-matrices describing the behaviour of light at the homogeneous boundaries of the core and cladding regions.#8231 - $15.00 USD Received 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 73092.2.Computing waveguide dispersion relationsResonant states reecting conned waveguide modes exist when the eld prole of light inside the core is unchanged following a round-trip from one edge of the core to the other. To determine the round-trip change in the eld, we examine the elds at the core-cladding interface and divide the propagation into two distinct phases. In the rst phase, an incident set of plane-waves E in one of the interface planes strikes the cladding producing a reected set of waves E+ . L L Since the cladding is semi-innite and the frequencies of interest are inside the the photonic crystal’s stopband, any incident Bloch mode will have decayed to zero at the core-cladding interface, hence E+ = 0. This condition enables the incident and reected elds to be related B simply through E+ = T12 T1 E . In the second phase, the rst reected set of waves propaL 22 L gates through to the other side of core, reects off the cladding, and propagates back through
the core. The relationship between incident E+ and reected EL waves in this case is given by L
+ EL = [S12 + S11 T12 (T22
S21 T12 )1 S22 ]EL . Applying the resonant state condition EL = E L yields the following: E = S12 + S11 T12 (T22
S21 T12 )1 S22 T12 T1 E . L 22 L (9)We can solve Eq. (9) as an eigenvalue problem noting that resonant states for the waveguide occur when an eigenvalue of the system equals one. For waveguides that are symmetric about the centre plane of the core, guided modes can be divided into even and odd classes with respect to the mirror plane normal to the slab. For resonant states in this case, the eld proles in the two core-cladding interfaces must satisfy E+ = ±E and E+ = ±E . Applying symmetry considerations to the scattering matrix and L R R L propagating the plane waves across to the other side of the core, we arrive at the following guided mode condition: E = ± 1
S12 T12 T1 L 221S11 T12 T1 E . 22 L(10)Eq. (10) can also be solved as an eigenvalue problem with even and odd guided modes obtained for eigenvalues equal to +1 and 1, respectively. To implement the eigenvalue equations above, we obtain the T-matrices governing the reection from the claddings and the S-matrices describing propagation through the core at a given k0, and a number of frequency points inside the stopband. These matrices have no variable elements and thus Eqs. (9) and (10) can be solved using standard computational techniques. After calculating eigenvalues over a series of frequencies, we can interpolate their phase and magnitude to nd the resonant states of the system at a very ne frequency resolution. In practice, the matrices in Eqs. (9) and (10) are well-conditioned and their dimensions are not large, with several hundred elements in 3D simulations and under 100 elements in 2D, enabling exact eigenvalues to be obtained in seconds. The eigenvalues from Eq. (9) have magnitudes that typically vary little over the stopband while those of Eq. (10) tend to vary rapidly. Both sets of eigenvalues usually have phases that vary smoothly over the stopband. As a result, we prefer to employ Eq. (9) for our simulations since interpolation is simpler and faster, and use Eq. (10) to gain general information about the symmetry of the modes. It is also possible to directly combine the S- and T-matrices dened in subsection 2.1 to compute dispersion relations. In this approach, we determine the transmission through the claddingcore-cladding structure in the direction normal to the waveguide. This method, however, proves to be inefcient since Bloch modes and scattering matrices at many frequencies are required to discern the very abrupt peaks in transmission that mark resonant states.#8231 - $15.00 USDReceived 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 73102.3.Core and cladding modicationsOnce the scattering and transfer matrices for different core and cladding regions are obtained, waveguide elements can be combined arbitrarily to yield new designs. The scattering matrices for core regions consisting of several rows of defects can be formed from those of a single row defect using the scattering matrix recursion formulae [25, 26]. Another degree of freedom can be achieved by displacing the cladding and core elements along the waveguide propagation direction. For core region unit cells, the translation is performed by introducing a phase shift into the incoming and outgoing plane wave coefcients to compensate for the relative positions
of neighbouring unit cells. For a translation r0, = x0 x+z0 z in the xz-plane, the shifted scattering matrix S is related to the original matrix S through: S = eiG r0, 0 eiG r0, 0 S e+iG r0, 0 0 e+iG r0, , (11)where e±iG r0, is a diagonal matrix with diagonal elements exp(±iGmn,
r0, ). Similarly, translations of the cladding unit cells are accomplished through the following operation: T = e+iG r0, 0 0 e+iG r0, T, (12)where T and T are the initial and shifted matrices, respectively. 3. ResultsIn this section, we apply the IDM to a variety of different photonic crystal waveguide designs. The bulk crystal in each of these waveguides consists of a triangular lattice of air holes etched in a high index slab. In our rst example, we determine the response of a full 3D planar photonic crystal waveguide with a defect row of air holes. In the subsequent examples, we focus on 2D waveguides to facilitate comparison with other methods. 3.1. 3D slab waveguideThe properties of photonic crystal waveguides can be tuned using different air hole defects inside the core. In the waveguide we study here, the bulk crystal consists of holes of radius 0.29a, where a is the lattice constant, etched into a dielectric membrane of thickness 0.6a and index 3.4. This particular crystal structure has been successfully fabricated and used in waveguides [27] and high-quality factor cavities [28, 29]. It possesses a band gap ranging from 0.261c/a to 0.331c/a for even modes with respect to the in-slab mirror plane. The waveguide core is formed by reducing the radii of a row of holes to 0.2a, establishing a line defect capable of guiding light (Fig. 2(a)). To determine the properties of the waveguide, we obtain the Bloch modes and scattering matrices for its unit cell components. The bulk crystal √ Bloch modes are simulated using a computational cell of height 4a and length 3a with orthogonal axes (Fig. 1) and an expansion of 1575 plane waves. The resulting eigenvalue equation is solved using an implementation of the implicitly restarted Arnoldi method (IRAM) [30] tuned to nd the Bloch modes with the slowest decay rates. The matrices in these computations are well-conditioned and generally require fewer than 20 IRAM iterations to obtain Bloch modes converged to machine precision. The core region scattering matrices are obtained from a cell of √ the same height but of length 3/2a with an expansion of 128 diffracted waves. To compute the dispersion relation (Fig. 2(b)), we form the transfer matrix using seven to ten evanescent Bloch modes and employ Eq. (9) to nd the resonant states. In practice, it is rare to nd resonant states where an eigenvalue is identically equal to one. Instead we employ a#8231 - $15.00 USD Received 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 7311(a)0.4 IDM Even Band IDM Odd Band MPB 0.35 Stopband edge Frequency (c/a)0.3Light line0.25Stopband edge0.2 0.2 0.25 0.3 0.35 0.4 Propagation constant (2π/a) 0.45 0.5(b)Fig. 2. (a) The slab waveguide with a row of defect air holes. The structure has index 3.4 and height 0.6a with 0.29a and 0.2a radius holes in the cladding and core, respectively. (b) Dispersion relation for the slab waveguide.weaker condition requiring the eigenvalue to be purely real with a magnitude greater than 0.6. By asserting that eigenvalues are real, we ensure that the guided mode prole remains identical after a round-trip even if its amplitude may change. Using this condition, we nd that the waveguide dispersion relation agrees very closely with results obtained using the MIT Photonic Bands (MPB) software package [3], with deviations in frequency of 2% of the stopband width in the worst case and 0.9% on average. The weaker eigenvalue criterion we employ here is general and has been applied to waveguides with different core radii. To emphasize the performance advantages of the IDM, we note that the eigenvalue operation in Eq. (9) is performed on matrices of dimension 256 in our method while the analogous supercell structure in MPB requires the eigenvalues of matrices approximately two orders of magnitude larger. Given the accuracy of our results, the lower than expected magnitudes of resonant state eigenvalues are likely the result of coupling between simulation cells above and below the waveguide resulting from the articial vertical periodicity we employ. Since the magnitude of the eigenvalues is related to the power retained in the waveguide over a round-trip about the core, it could be reduced signicantly by any vertical coupling losses out of the waveguide. Further examining the results, we nd that the bands with even and odd symmetry with respect to the vertical mirror plane along the core are marked by eigenvalues with largely different magnitudes. The set of even modes with eld energy concentrated at the centre of the core have eigenvalues of#8231 - $15.00 USD Received 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 7312magnitude 0.6 to 0.7; in contrast, those with odd symmetry have magnitudes of 0.88 to 0.99 and their energy concentrated away from the centre. Since the air holes are located at the centre of the core, the index connement of the even modes must be weaker than that experienced by the odd ones. Consequently, the lower magnitudes we observe can be related to the degree of coupling to slabs above and below the core. This unwanted coupling could be reduced by further extending our simulation cells in the vertical direction. Nevertheless, our results remain very accurate even with these vertical losses. For similar simulations on 2D crystals of innite vertical extent, eigenvalues of resonant states have magnitudes much closer to one, as described in subsection 3.2. 3.2. Tuned air hole core waveguidesWhile the waveguide studied in subsection 3.1 has a very simple structure, its properties can be altered signicantly by modications to the radii and placement of the holes in the core region. We investigate these effects by obtaining the dispersion relations for a series of 2D waveguides with lattice mismatching of zero to 0.5a (Fig. 3(a)) and air hole radii in the core ranging from zero to 0.425a (Fig. 4(a)). The photonic crystal system we use for this and the remaining examples consists of a triangular lattice of air holes of radius 0.3a in a dielectric of effective index 3.4 with a band gap over 0.211c/a to 0.279c/a in the TE polarized modes. This set of simulations employs the 2D unit cell analogues of those used in the 3D slab waveguide. Bloch modes are expressed using an expansion over 512 plane waves and solved using an IRAM implementation [30]. Scattering matrices were computed for sets of 16 diffracted waves.0.29 0.28 0.27 Frequency (c/a) 0.26 0.25 0.24 0.23 0.22 0.21 0.2 0 0.1 0.2 0.3 Propagation constant (2π/a) 0.4 0.5 offset = 0.00a offset = 0.25a offset = 0.50a MPB Odd modes Stopband edgeEven modesStopband edge(a)(b)Fig. 3. (a) Schematic of the waveguides simulated illustrating the shifting of core holes from their regular lattice sites. (b) Dispersion relations for the waveguides with core radius 0.2a as the core-cladding offset is varied.Figure 3(b) shows a sample dispersion relation for a core region of radius 0.2a with holes displaced zero, 0.25a, and 0.5a from their regular lattice positions (Fig. 3(a)). The shift in core position was performed using Eq. (12). For these 2D dispersion relations with no vertical coupling losses, Eq. (9) returns modes with real eigenvalues averaging 0.98 with a standard deviation of 0.02. We attribute the deviation in eigenvalue magnitudes to numerical errors created in combining transfer and scattering matrices. Nevertheless, our results agree with MPB to within 1.9% of the stopband width in the worst case. By changing the offset of the core holes relative to the cladding, we can effectively translate the dispersion relation in frequency and#8231 - $15.00 USD Received 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 7313propagation constant space. Increasing the offset moves the core and cladding air holes closer to one another, lowering the frequencies of the bands and reducing their propagation constants. This core offset alone can change the position of the zero group velocity frequency points of the odd bands from 0.2602c/a to 0.2755c/a or 22.6% of the band gap.0.29 0.28 0.27 Frequency (c/a) 0.26 0.25 0.24 0.23 0.22 0.21 0.2 0 0.1 0.2 0.3 Propagation constant (2π/a) 0.4 0.5 radius = 0.150a radius = 0.175a radius = 0.200a MPB Even modes Odd modes Stopband edgeStopband edge(a)(b)Fig. 4. (a) Schematic of the waveguides simulated illustrating the change in core hole radius. (b) Dispersion relations for the waveguides as the core hole radius is varied with no lattice offset.The second signicant degree of freedom in this class of waveguides is obtained by varying the radii of the air holes in the core. In Fig. 4(b) we present the dispersion relations for core regions of radius 0.15a, 0.175a, and 0.2a with no lattice offset. Even for this small gradation in radii, we observe noticeable shifts in the dispersion relations. The bands in this case are translated almost directly up in frequency as the radius of the core holes increases, lowering the core’s effective index. The zero group velocity point in the odd band shifts from 0.2662c/a to 0.2755c/a in frequency with only a very small change in propagation constant. For waveguides with defect row radii ranging from zero to 0.225a, the dispersion relations share the same overall band structure with smoothly descending even bands and concave down odd bands, regardless of the offset of the core holes. For defect radii greater than this, the odd bands grow atter and are pushed out of the stopband along with the even modes. Even with this simple waveguide system, we can tailor a potential device to exhibit particular properties by varying its design over the two degrees of freedom. For instance, the location of the zero group velocity state in the odd band can be varied from 0.2513c/a to 0.2809c/a in frequency and 0.209 × 2π /a to 0.318 × 2π /a in the propagation constant as illustrated in Fig. 5. As one would expect, the effect of the lattice mismatch between core and cladding increases as the radius of the core holes increases leading to signicant changes in waveguide response at the largest radii. 3.3. Optimized waveguidesAfter building up a library of transfer and scattering matrices representing bulk crystals and core regions, the computational speed of the IDM enables us to design waveguides with desired properties simply by scanning over a large number of potential structures. We demonstrate this through the design of a 2D photonic crystal waveguide with a single-moded frequency range spanning 88.7% of the bulk crystal band gap.#8231 - $15.00 USD Received 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 7314(a)Core hole offset (a)0.5 0.4 0.3 0.2 0.1 0 0 0.05 0.1 0.15 0.2 Core hole radius (a)Frequency (c/a) 0.285 0.28 0.275 0.27 0.265 0.26 0.255 0.25(b)Core hole offset (a)0.5 0.4 0.3 0.2 0.1 0 0 0.05 0.1 0.15 Core hole radius (a) 0.2Propagation constant (2π/a) 0.32 0.3 0.28 0.26 0.24 0.22 0.2Fig. 5. (a) Frequency and (b) propagation constant of the zero group velocity point for odd waveguide modes as a function of the core hole radius and lattice offset.To shorten our search for optimal structures, we settle on a general waveguide design with an asymmetric core consisting of three rows of air holes in a rectangular lattice. The Fresnel coefcients of the bulk photonic crystal and scattering matrices for single row core regions ranging from homogeneous to holes of radius 0.425a are available from the results of subsection 3.2. These stored scattering matrices can be combined arbitrarily using the scattering matrix recursion formulae [25, 26] to represent the three row core region. In the rst step of our design process, we scan over a coarse set of nine different core radii to observe the range of dispersion relations provided by our general waveguide structure. We obtain the series of dispersion relations using LAPACK [31] subroutines to generate the eigenvalues of Eq. (9) at 64 frequency points per propagation constant, and employ an unoptimized root nding algorithm to solve for the resonant states. The acquisition of 405 dispersion relations computed at 21 propagation constants takes approximately two hours on a Pentium IV running at 2.26GHz using 650MBytes of memory, while the calculation of one such dispersion relation in MPB at the same resolution would take approximately 15 minutes or more. After analyzing the rst set of dispersion relations, we develop a qualitative understanding of the effect of different core designs on the respo in particular, we nd that positive group velocity bands spanning a large frequency range are established by core air holes with radii greater than or equal to that of the bulk crystal. With a general design rule in place, in the second scan we simulate a smaller range of core#8231 - $15.00 USD Received 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 73150.29 0.28 0.27 Frequency (c/a) 0.26 0.25 0.24 0.23 0.22 0.21 0.2 0 0.1 0.2 0.3 Propagation constant (2π/a) 0.4 0.5 Stopband edge f = 0.0598c/a Stopband edgeFig. 6. The dispersion relation for the optimized structure with a single-moded frequency range over 88.7% of the stopband width. Inset: schematic of the optimized waveguide with a three row core of radii 0.4a, 0.425a, and 0.3a.radii at smaller increments from 0.3a to 0.425a and search for dispersion relations that are single-moded over large frequency ranges. The majority of dispersion relations from this set of waveguides share a rising band over the full range of propagation vectors as expected, but also have a second band in the lower half of the stopband that limits the single-moded frequency range. By making small changes to the index prole of the core, we can expand the range of frequencies covered by the rising band while pushing the lower band out of the stopband. An analysis of the dispersion relations of the second set of waveguides reveals that a structure consisting of a centre hole of radius 0.425a bordered by holes of radius 0.4a and 0.3a produces the largest single-moded span (Fig. 6). For this asymmetric two row core structure, the rising band covers 94.5% of the band gap and the single-moded region is determined principally by the lower band, which terminates at 0.2186c/a, making the waveguide single-moded for 88.7% of the stopband. The dispersion relation for this structure obtained through our method agrees with MPB simulations to within 2.2% of the stopband width in the worst case. The width of this waveguide is particularly interesting because it is wider than most other single-moded designs and could be used to couple light from bres more efciently into photonic crystal devices. 3.4. Directional couplerDirectional couplers are important devices formed by two waveguides separated by a thin barrier region. In these structures, light travelling in one waveguide can couple to the other enabling power to be split between guides. The distance required for power to be completely transferred from one waveguide to another is known as the coupling length. Minimization of this characteristic length is crucial to the operation of microphotonic circuits. In previous work with our complex Bloch modes [17], we had been limited to simulating couplers with relatively thick barriers to ensure the elds in these regions could be accurately represented by Bloch modes. As a result, the two waveguides could not interact strongly, leading to large coupling lengths in these devices. Here we circumvent this limitation by representing the barrier region as part of a single supercore structure described by scattering matrices.#8231 - $15.00 USD Received 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 73160.29 0.28 0.27 Frequency (c/a) 0.26 0.25 0.24 0.23 0.22 0.21 0.2 0 0.1 0.2 0.3 Propagation constant (2π/a) 0.4 0.5 Stopband edge β = 0.2092π/a Single waveguide bands Even supermodes Odd supermodes Stopband edgeFig. 7. Dispersion relation for the directional coupler and the equivalent single waveguide. Inset: schematic of the coupler with a pair of identical waveguides formed by two rows of holes of radius 0.4a separated by a row of holes of radius 0.3a.We demonstrate this approach on a novel structure designed with the insights of subsection 3.3 in mind. Our coupler consists of two identical waveguides separated by a barrier made up of a single row of air holes (Fig. 7). The waveguides are formed from two parallel rows of 0.4a radius air holes such that the entire coupler consists of ve defect rows with the holes arranged in a rectangular lattice. When combined with the barrier region, these two row core waveguides are very similar in structure to the optimized waveguide described in subsection 3.3. Like the tuned structure, the coupler’s equivalent single waveguide possesses a rising singlemoded band over a relatively large frequency range, 78.0% of the band gap in this case (Fig. 7). In this coupler design, the barrier serves a dual purpose in separating the waveguides and tailoring their dispersion relations. We calculate the supermodes of the entire directional coupler with a scattering matrix representing the ve rows of the core and proceed in a manner similar to the one described in Ref. [32]. The modes of the device come in complementary pairs of even and odd modes at frequencies above and below the bands of the individual waveguide structure (Fig. 7). As a result of this mode symmetry, a superposition of modes can propagate along the coupler with different propagation constants, shifting power from one waveguide to the other as they beat in and out of phase. At a frequency of 0.25c/a, the even and odd supermodes in the structure have propagation constants βodd and βeven that differ by 0.1046 × 2π /a. The beat length LB of the directional coupler at this frequency is given by: LB = 2π = 9.56a, |βodd
βeven | (13)resulting in a coupling length of 4.78a, which agrees to within 2.7% of the value predicted by MPB. This coupler also has a fairly large potential range of operation from 0.225c/a to 0.262c/a or 54.9% of the stopband width, where only one pair of symmetrical supermodes exists. In addition, the bands have roughly constant group velocities over much of this range enabling the coupling lengths to remain fairly similar over a signicant frequency span.#8231 - $15.00 USD Received 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 73174.ConclusionsWe have demonstrated an efcient interface diffraction method of simulating planar photonic crystal waveguides and couplers with arbitrary defect regions. The IDM merges two contrasting approaches C the plane wave expansion method in reciprocal space and the scattering matrix method in real space C and takes advantage of their respective strengths to accurately determine the eld proles inside the cladding and core regions of a device. Since IDM simulations are conducted only on the unit cell elements of a structure, Bloch modes and scattering matrices can be computed in advance and used to obtain results over an order of magnitude faster than supercell techniques, in which simulations must be run on a case-by-case basis. Further, the IDM has a conceptually simple framework and is quite general, making few assumptions about device geometry and structure. It can be applied to narrow or wide defect regions with abrupt or graded changes in structure, all within the same theoretical approach without large increases in computation time. Through the IDM individual waveguide elements can be combined and positioned rapidly, facilitating the design of waveguides and couplers with elaborate defect regions. Moreover, with its high computational efciency, the IDM can be used to optimize photonic crystal structures over several degrees of freedom, ensuring they possess particular properties such as large single-moded frequency ranges, high group velocities, and small coupling lengths. In future work, the IDM could be extended by employing other techniques [33, 34] to enable mode-matching at the top and bottom surfaces of the slab to obviate the current need for simulation cells extended in the vertical direction. For waveguides etched in semiconductor slabs, the expansion could also be done in terms of the modes of the unpatterned slab [35] rather than plane waves. These extensions could further increase the method’s computational speed and eliminate any coupling between vertically adjacent unit cells. Mode-matching at the top and bottom of the slab [33, 34] could also be used together with the IDM to determine the radiation modes and radiation losses of planar photonic crystal structures, which directly determine the quality factor of resonators in such structures. The effect of material loss on the properties of waveguide modes could be incorporated into the method through complex dielectric constants. Both the plane wave expansion and scattering matrix methods employed in the IDM can accommodate such absorptive media. The guided mode eigenvalues in Eqs. (9) and (10), with magnitudes brought below one from the material loss, can be used to determine the propagation losses of the waveguide modes. Although we illustrated this method using waveguides, which are among the most common practical photonic crystal devices, the IDM is not limited to just this class of structures. We have already used the IDM to simulate the interaction between waveguides and cavities created from homogeneous defects [17]. By employing the IDM to determine the mode proles of individual devices, we could calculate a new class of scattering matrices describing the behaviour of each device. These matrices could be combined through the approach described in Ref. [36] to model the scattering between multiple devices, including that between two waveguides or a waveguide and a cavity. Waveguide bends could be simulated by dividing the device into two waveguide segments coupled to a resonator joining the segments. The IDM could be applied to the individual waveguides and to the resonator [16]. Transmission and reection at the waveguide bend could then be evaluated by calculating the coupling between these three elements. We foresee that the IDM could also be used to simulate superprisms and to optimize emission from photonic crystal waveguides. Acknowledgments The authors acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) through the Research Network programme, “Agile All-Photonic Networks (AAPN).”#8231 - $15.00 USD Received 31 July 2005; revised 29 August 2005; accepted 30 August 2005(C) 2005 OSA19 September 2005 / Vol. 13, No. 19 / OPTICS EXPRESS 7318
更多搜索：
赞助商链接
All rights reserved Powered by
文档资料库内容来自网络，如有侵犯请联系客服。