A rigorous full-vector analysis based on the finite-difference mode-matching method is presented for three-dimensional optical wave propagation problems. The computation model is facilitated by a perfectly matched layer (PML) terminated with a perfectly reflecting boundary condition (PRB). The complex modes including both the guided and the radiation fields of the three-dimensional waveguide with arbitrary index profiles are computed by a finite-difference scheme. The method is applied to and validated by the analysis of the facet reflectivity of a buried waveguide and the power exchange of a periodically loaded dielectric waveguide polarization rotator.
©2008 Optical Society of America
The advances in nanofabrication technologies of high index contrast optical waveguide devices hold great promise for large-scale, high-density photonic integrated circuits for a wide range of applications [1–4]. Miniaturized photonic devices and high-density integrated photonic chips lead to cost reduction of optical systems and make optical functionality more affordable. There have been a large number of such high-performance photonic devices demonstrated. Examples include Bragg reflection filters , silicon optical modulators , silicon Mach-Zehnder modulators , just to name a few. One of distinct features of the high index contrast optical waveguides and high-density integrated optical circuits is the strong interaction between the optical field and the structure. Consequently, polarization effects and hybrid features of the vector optical field and the successive reflections due to the waveguide discontinuities can no longer be neglected. Among the widely used modeling tools, coupled mode theory (CMT) provides a physically intuitive picture and is considered as mathematically simple and accurate for relatively weakly guided optical waveguides [8–10]. In practice, numerical methods are more widely employed to deal with the field propagation problems. The application of finite difference time domain (FDTD) method has been encouraged to model electromagnetic problems by its comparative accuracy and versatile features but limited by intensive computations and storage requirements. Bi-directional beam propagation method (BiBPM) based on the scattering operators can correctly model the evanescent waves by using the complex propagation operators, and the robustness and effectiveness are enhanced by the pre-conditioners. However, the phase errors associated with BiBPM depend critically on the transverse discretization schemes . Using relative fine mesh may alleviate the problem but leads to a high computation effort and memory requirements.
Mode-matching method (MMM) has been known as an efficient and rigorous method in dealing with waveguide structures with longitudinal discontinuities especially for periodical structures [13–18]. In MMM approach, the structure is divided into multiple homogenous sections. For each section, the total propagating field is expressed as the eigen-modes superposition. By applying the continuity conditions of tangential components of electric and magnetic field at the interfaces and the mode orthogonality relations, the different sections can be related by a scattering matrix. The entire reflection and transmission are calculated by cascading the scattering matrices. The theoretical formulations of the mode-matching method have been developed in the literatures. The area of continuing interests and research activities has mainly concerned the treatment of the radiation fields for simulation of open waveguide structures. One of the widely used approaches was by enclosing the object with a sufficient large metal box so as to discretize the continuous radiation fields. The shortcoming associated with the computation model is that a large number of the modes in the box have to be considered to ensure adequate accuracy. Further, a detrimental oscillation error has been observed with change of the size of the computation window that is slow to diminish. Another scheme for implementation of the MMM is to represent the radiation field in terms of Fourier expansion. In this treatment, a periodic boundary condition is implied and spurious reflections from the boundaries need to be dealt with .
Recently, a new scheme for the mode-matching method was proposed and demonstrated by applying perfectly matching layers (PML) terminated with a perfectly reflecting boundary (PRB) condition [19, 20]. In the new MMM scheme, the radiation fields are represented in terms of a set of complex modes, some of which resemble the conventional leaky modes and others associated with the interaction between the PML media and the reflecting numerical boundaries. The mode spectrum is therefore split into the guided modes and complex modes which possess the normal mode features such as normalization and modal orthogonality. The seemingly paradoxical application of both the PML and PRB in the mode-matching method has in fact overcome one of the main challenges associated with this traditional method, i.e., the desire for discrete, orthogonal, and normalized modes to represent radiation fields and the need for elimination and reduction of spurious reflections from the edges of the finite computation window. So far, the finite-difference based numerical implementation of MMM for the analysis of three-dimensional structures has been reported under the semi-vector approximation . It is true that the semi-vector approximation is sufficient for many practical optical waveguides in which the modal fields are predominantly linearly polarized and the coupling between the field components is negligible. However, the modal hybridness in high index contrast waveguides is usually pronounced and can not be ignored. In order to accurately model such devices—polarization rotator, for example, the semi-vector approximation is deficient and the full vector approach is mandatory [22–24].
In this work, we develop the finite-difference based full-vector MMM for the three-dimensional propagation problems. The theoretical formulations are briefly described in Section 2. The applications in the three-dimensional propagation problems are discussed in Section 3, followed by a summary in Section 4.
Considering the propagation of optical waves in a three-dimensional waveguide structure, the propagation direction is assumed to be z direction and the transverse cross section is in x,y plane. The waveguide is characterized by M distinct z-invariant index profiles. n(x,y,z=nj(x,y) for zj-z<z<zj where j=1, 2,⋯,M. Further we suppose that the waveguide is z-invariant for z<0 and z>l, where l is the length of region including the multiple discontinuities. The permittivity and permeability of the vacuum are denoted as ε 0 and µ0, respectively. Suppressing the time variation factor ejωt, the vector wave equation for each piecewise uniform waveguide can be expressed in terms of transverse electrical field as 
where the differential operators are defined as
Here k0 is the free space wave number; β is the propagation constant and n=n(x,y,z). The full-vector equations contain all the vector properties of the electromagnetic field. The items Pxx≠Pyy will cause the polarization dependence whereas Pxy≠0 and Pyx ≠0 induce the polarization coupling. Equation (1) is solved by finite-difference schemes. The complex coordinate stretching formulation of PML boundary condition is obtained by replacing x,y in the original wave equations with the complex coordinates x̃,ỹ :
Here ζ0 denotes the initial PML interface, d is the PML thickness, and m is the order of the polynomial profile which determines the evolution of σ from minimum value of zero at the initial PML interface to a maximum value at the computation boundary. Typically m is chosen as 2 to generate a parabolic PML profile. The decay factor σζmax is related to the reflection coefficient R at the computation boundary by
where c is the speed of light in the free space.
The transverse magnetic fields in terms of the transverse electric fields components are:
where is the characteristic admittance of the free space, Neff=β/k0 is the modal index. For general media, the time average power for each mode can be normalized as:
where t denotes the transverse components, and k denotes the mode number. The expression is reduced to
for lossless media.
All the modes including the guided modes and the complex modes are orthogonal which can be expressed as
for lossless media, and
for general media. Here δmk is defined by
where m,k stand for the mode index and Pn is referred to as Petermann factor.
Considering a wave propagating from z<0, at the first discontinuity (z=0), part of the wave will be reflected into the region z<0 and part of the wave will be transmitted into the region z>0. Assuming that N modes are significant in both sides, the total transverse fields in terms of the forward waves (denoted by “+”) and backward waves(denoted by “-”) can be expressed as:
Here, the superscripts A and B denote the waveguide sections z<0 and the first section as z>0, respectively; subscript t denotes the transverse component; the subscripts i and k denote mode index in waveguide A and waveguide B, respectively. e⃗t and h⃖t are the transverse electric and magnetic fields, and a±(b±) are the amplitudes of forward and backward waves in waveguide A (B). The continuity of the tangential electric and magnetic fields at z=0 gives rise to
where the inner product of the field vectors is defined by
In the matrix form, we have
where p± and q± are defined as mode amplitude vectors (b ± 1,b ± 2,⋯b ± N) and (a ± 1,a ± 2,⋯a ± N), respectively. TAB denotes the transfer matrix across the discontinuous interface z=0. We may use scattering matrix to eliminate the numerical instability resulting from the conventional transfer matrix algorithm [27, 28],
The elements of scattering matrix SAB are given by
The complete scattering matrix which characterizes a unit cell can be acquired by cascading the scattering matrix across the junction and propagation matrix within the z -invariant sections. Taking the example of section B, the propagation matrix is given by,
where Λb is the section length. The elements of scattering matrix of a unit cell are given by
where u denotes the unit cell. The entire scattering matrix for the structures with multiple discontinuities can be obtained by cascading the scattering matrices layer by layer. Assuming that the scattering matrix linking the first interface with the (m-1)th interface is given by
and the scattering matrix linking the (m-1)th with m th interfaces is given by
the scattering matrix linking the first interface with the m th interface can be easily obtained by
Here I is the identity matrix.
It is noted that the computation efforts could be significantly reduced when the structure contains periodic cells such as Bragg reflectors. The key idea is to find the equivalent scattering matrix of s cells. The scattering matrix of 2 s periods is obtained by cascading the two identical scattering matrices of s cells. Taking the advantage of symmetric features, for any M periods, only log2 M+1 or log2 (M+1) steps (depends on whether M is an integer power of 2) needed, therefore, this method can speed up the calculation significantly.
3. Simulation results and discussion
The full-vector mode-matching method developed in the previous section will be applied to analyze the facet reflectivity of a buried waveguide structure (Fig. 1). The transmission medium is assumed to be air and the operation wavelength is 1.30 µm.
The modal facet reflectivities of the fundamental quasi-TE and quasi-TM modes against the core width wx have been calculated. The refractive indexes of the core and the cladding media are 3.54 and 3.17, respectively. The core thickness t is 0.4 µm. The simulation parameters are summarized as follows: the computational window size is 1.4µm×2.4 µm (thickness ×width) including a perfectly matched layer (PML) of 0.2µm thickness on each side; the horizontal and vertical mesh sizes are 50 nm. Three hundred modes have been employed to increase the accuracy.
The modal reflectivities of the quasi-TE and quasi-TM modes against the core width of the buried waveguide are shown in Fig. 2. Comparing with the reported results from full-vector and semi-vector approaches [30, 31], it is observed that the results from full-vector mode-matching method are in good agreement with that from the full-vector Bi-BPM. However, it is noted that there is a prominent discrepancy between the semi-vector and the full-vector approaches as the width and thickness of the structure becomes comparable. In other words, the semi-vector method is valid for large structure aspect ratio r (defined by max (wx,t)/min (wx,t)).
We also investigate the convergence behavior of the facet modal reflectivity of the buried waveguide as a function of mode numbers. The simulation parameters are summarized as follows: the width and thickness of the cross section are fixed at 0.4 µm; the computation window is set to 1.4µm ×1.4 µm (thickness ×width) including a perfectly matched layer (PML) of 0.2 µm thickness on each side. The facet reflectivity of the quasi-TE mode as a function of mode numbers is shown in Fig. 3. It is observed the modal reflectivity oscillates then converges to stable.
To validate the full-vector finite-difference mode-matching method for multiple discontinuities, the polarization rotator which has been widely used in optical systems is simulated. The structure consists of a buried waveguide (λ=20mm) with periodic perturbations (Fig. 4). The structure parameters are summarized as follows: The refractive indexes of buried waveguide and the loadings are εr=2.8; the width (W) and the height (H) of the buried waveguide are 13mm and 6.5 mm, respectively. The width of the loading (w) is chosen to be half of the waveguide structure, i.e., w=6.5mm. The height of the loaded strip (h) is chosen to be 2mm or 3mm. The phase matching condition between the fundamental quasi-TE and quasi-TM modes in the loaded sections gives rise to the loading length d=π/(βTE-βTM), where βTE and βTM are the propagation constants of the quasi-TE and quasi-TM modes, respectively.
The parameters used in the simulations are summarized as follows: the computational window size is 35×35mm including a perfectly matched layer (PML) of 6 mm thickness on each side; the mesh sizes are fixed asΔx=Δy=0.5mm, 60 modes are employed in the mode matching scheme. The power exchange as the function of loadings for h=2mm and h=3mm are shown in Fig.5 and Fi.g.6, respectively. It is observed that the simulation results agree well with the experimental results from . The overall loss is proportional to the loading numbers which results from the radiation loss at the junctions.
We have developed a highly accurate and efficient full-vector mode-matching method based on finite-difference method. The perfectly matched layer is employed to reduce the parasite boundary reflections. The average power and the orthogonality of the guided, as well as complex modes are discussed in detail. The scattering matrix for multi-discontinuities is derived. The developed method is applied to analyze the facet modal reflectivity of a buried waveguide. Comparison with the reported results shows that the semi-vector method is only valid for large structure aspect ratio. We also investigate the modal reflection behavior for varied mode numbers. It is noted the reflectivity oscillates first and converges as the mode numbers increased. The capability and accuracy of the developed method in multiple three-dimensional discontinuities have been validated by investigating the power exchange in periodically loaded waveguide polarization rotators.
The authors would like to express their appreciation to Dr. Xun Li and Mr. Kai Jiang for their insightful discussions.
References and links
1. S. Janz, P. Cheben, D. Dalacu, A. Delage, A. Densmore, B. Lamontagne, M. J. Picard, E. Post, J. H. Schmid, P. Waldron, D. X. Xu, K. P. Yap, and W. N. Ye, “Microphotonic elements for integration on the silicon-on-insulator waveguide platform,” IEEE J. Sel. Top. Quantum Electron. 12, 1402–1415 (2006). [CrossRef]
2. M. Lipson, “Guiding, modulating, and emitting light on silicon - Challenges and opportunities,” J. Lightwave Technol. 23, 4222–4238 (2005). [CrossRef]
4. C. Manolatou, S. G. Johnson, S. H. Fan, P. R. Villeneuve, H. A. Haus, and J. D. Joannopoulos, “High-density integrated optics,” J. Lightwave Technol. 17, 1682–1692 (1999). [CrossRef]
5. T. E. Murphy, J. T. Hastings, and H. I. Smith, “Fabrication and characterization of narrow-band Bragg-reflection filters in silicon-on-insulator ridge waveguides,” J. Lightwave Technol. 19, 1938–1942 (2001). [CrossRef]
6. O. Limon, A. Rudnitsky, Z. Zalevsky, M. Nathan, L. Businaro, D. Cojoc, and A. Gerardino, “All-optical nano modulator on a silicon chip,” Opt. Express 15, 9029–9039 (2007), http://www.opticsexpress.org/abstract.cfm?uri=OE-15-14-9029. [CrossRef] [PubMed]
7. L. Liao, D. Samara-Rubio, M. Morse, A. S. Liu, D. Hodge, D. Rubin, U. D. Keil, and T. Franck, “High speed silicon Mach-Zehnder modulator,” Opt. Express 13, 3129–3135 (2005), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-13-8-3129. [CrossRef] [PubMed]
8. W. Streifer, D. Scifres, and R. Burnham, “Coupled wave analysis of DFB and DBR lasers,” IEEE J. Quantum Electron. 13, 134–141 (1977). [CrossRef]
9. A. Yariv, “Coupled-mode theory for guided-wave optics,” IEEE J. Quantum Electron. 9, 919–933 (1973). [CrossRef]
10. A. M. Shams-Zadeh-Amiri, J. Hong, X. Li, and W. P. Huang, “Second- and higher order resonant gratings with gain or loss - Part I: Green’s function analysis,” IEEE J. Quantum Electron. 36, 1421–1430 (2000). [CrossRef]
11. A. Taflove and S. C. Hagness, Computational Electrodynamics : the finite-difference time-domain method (Artech House, Boston, 2005).
12. P. L. Ho and Y. Y. Lu, “A stable bidirectional propagation method based on scattering operators,” IEEE Photon. Technol. 13, 1316–1318 (2001). [CrossRef]
13. P. Bienstman and R. Baets, “Optical modelling of photonic crystals and VCSELs using eigenmode expansion and perfectly matched layers,” Opt. Quantum Electron. 33, 327–341 (2001). [CrossRef]
14. S. F. Helfert, A. Barcz, and R. Pregla, “Three-dimensional vectorial analysis of waveguide structures with the method of lines,” Opt. Quantum Electron. 35, 381–394 (2003). [CrossRef]
15. T. E. Rozzi, “Rigorous Analysis of Step Discontinuity in a Planar Dielectric Waveguide,” IEEE Trans. Microwave Theory Tech. 26, 738–746 (1978). [CrossRef]
16. H. Shigesawa and M. Tsuji, “Mode Propagation through a Step Discontinuity in Dielectric Planar Wave-Guide,” IEEE Trans. Microwave Theory Tech. 34, 205–212 (1986). [CrossRef]
17. E. Silberstein, P. Lalanne, J. P. Hugonin, and Q. Cao, “Use of grating theories in integrated optics,” J. Opt. Soc. Am. A 18, 2865–2875 (2001). [CrossRef]
18. C. Vassallo, “Finite-difference derivation of the reflectivity at output facets of dielectric waveguides with a highly diverging output beam,” J. Opt. Soc. Am. A 15, 717–726 (1998). [CrossRef]
19. H. Derudder, D. De Zutter, and F. Olyslager, “Analysis of waveguide discontinuities using perfectly matched layers,” Electron. Lett. 34, 2138–2140 (1998). [CrossRef]
20. H. Derudder, F. Olyslager, D. De Zutter, and S. Van den Berghe, “Efficient mode-matching analysis of discontinuities in finite planar substrates using perfectly matched layers,” IEEE Trans. Antennas Propag. 49, 185–195 (2001). [CrossRef]
21. K. Jiang and W. P. Huang, “Finite-difference-based mode-matching method for 3-D waveguide structures under semivectorial approximation,” J. Lightwave Technol. 23, 4239–4248 (2005). [CrossRef]
22. Y. Shani, R. Alferness, T. Koch, U. Koren, M. Oron, B. I. Miller, and M. G. Young, “Polarization Rotation in Asymmetric Periodic Loaded Rib Wave-Guides,” Appl. Phys. Lett. 59, 1278–1280 (1991). [CrossRef]
23. T. Ando, T. Murata, H. Nakayama, J. Yamauchi, and H. Nakano, “Analysis and measurement of polarization conversion in a periodically loaded dielectric waveguide,” IEEE Photon. Technol. Lett. 14, 1288–1290 (2002). [CrossRef]
24. F. J. Mustieles, E. Ballesteros, and F. Hernandezgil, “Multimodal Analysis Method for the Design of Passive Te/Tm Converters in Integrated Wave-Guides,” IEEE Photon. Technol. Lett. 5, 809–811 (1993). [CrossRef]
25. W. P. Huang and C. L. Xu, “Simulation of 3-Dimensional Optical Wave-Guides by a Full-Vector Beam-Propagation Method,” IEEE J. Quantum Electon. 29, 2639–2649 (1993). [CrossRef]
26. W. C. Chew, J. M. Jin, and E. Michielssen, “Complex coordinate stretching as a generalized absorbing boundary condition,” Microwave Opt. Technol. Lett. 15, 363–369 (1997). [CrossRef]
27. L. F. Li, “Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings,” J. Opt. Soc. Am. A 13, 1024–1035 (1996). [CrossRef]
28. L. F. Li, “Multilayer Modal Method for Diffraction Gratings of Arbitrary Profile, Depth, and Permittivity,” J. Opt. Soc. Am. A 10, 2581–2591 (1993). [CrossRef]
29. P. L. Ho and Y. Y. Lu, “A bidirectional beam propagation method for periodic waveguides,” IEEE Photon. Technol. Lett. 14, 325–327 (2002). [CrossRef]
30. H. A. Jamid and M. Z. M. Khan, “3-D full-vectorial analysis of strong optical waveguide discontinuities using Pade approximants,” IEEE J. Quantum Electron. 43, 343–349 (2007). [CrossRef]
31. K. Kawano, T. Kitoh, M. Kohtoku, T. Takeshita, and Y. Hasumi, “3-D semivectorial analysis to calculate facet reflectivities of semiconductor optical waveguides based on the bi-directional method of line BPM (MoL-BPM),” IEEE Photon.Technol. Lett. 10, 108–110(1998). [CrossRef]