In this paper, static and dynamic thermo-optical numerical analysis of a Digital Optical Switch (DOS), based on amorphous silicon waveguide and operating at the infrared communications wavelength of 1550 nm, are presented. The aim of our design is to achieve good performances in terms of cross talk and switching time, considering relaxed requirements for the realization of device: large cross section single mode waveguides and an angle between the output branches not too small. Using a low temperature difference between the two output branches, an optical switching with a crosstalk of 25 dB and a response time of the order of ten microseconds are obtained. The device, designed for low-cost photonic applications, could be easily integrated in silicon optoelectronic circuits.
©2006 Optical Society of America
Digital Optical Switch (DOS) has become, since its invention, a very attractive component for space switching in multi-wavelength optical communication system applications. For reliable operation, most of the other types of optical switches, which rely on effects such as interference, reflection, absorption, or transparency, need sophisticated feed-back loops and tight wavelength or temperature control. On the contrary, its step-like switch response makes the DOS highly insensitive to wavelength, polarization, and other physical parameters that may affect the device operation [1–4]. On the other hand, usually, DOS has quite small transverse dimensions and branch separation angles; this means that the realization process is complex and the length of device is big. Therefore device operation is affected by high coupling and propagation losses. A variety of switch arrays consisting of DOS’s have already realized using LiNbO3 [5–7], III–V semiconductor , polymers [9–12], silicon resin  and silica on silicon . An electro-optic DOS based on LiNbO3 is attractive due to its high performance. It is capable of changing its state extremely rapidly, typically in less than a nanosecond, being the switching time limited by the capacitance of electrode configuration. Electro-optic switches are reliable, but they pay the price of high insertion loss and possible polarization dependence. Polarization independence is possible but at the cost of a higher driving voltage, which in turn limits the switching speed. Moreover, they suffer from high fabrication cost. Thermo-optic switches are generally small in size but have the drawbacks of high-driving-power, limited integration density and high-power dissipation. On the positive side, optical performance parameters, such as crosstalk and insertion loss, are acceptable for many applications. Thermo-optic effect in polymer is one of the most attractive approaches, being characterized by low power dissipation, due to its large TO coefficient, and low-cost fabrication. The main drawback is poor response time.
In recent years, hydrogenated amorphous silicon (a-Si:H), grown by Plasma Enhanced Chemical Vapor Deposition (PECVD) has gained considerable attention due to its unique characteristics of transparency at the infrared wavelengths, the refractive index tunability, and finally technological compatibility with all microelectronic processes. As a consequence, several a-Si:H based optoelectronic devices, such as light emitting diodes , photodetectors  and optical modulators  have been successfully fabricated. It is worth to note that thin films of a-Si:H can be deposited by chemical vapor deposition techniques (CVD). In this case the process takes place at low temperature, usually below 250 °C, and therefore layers of thin material can be virtually deposited on any substrate. Moreover, this technique allows the incorporation of hydrogen atoms in the film, which saturate the silicon dangling bonds and clean the semiconductor forbidden band from a conspicuous amount of undesired states. The resulting a-Si:H alloy can easily show energy gap spanning from 1.6 to 1.9 eV, depending on the H atom concentration, and therefore on the degree of saturation of the silicon bonds. Moreover, lower absorption coefficients, especially at the infrared wavelengths in the range of 1.3 and 1.55 µm, are generally shown by a-Si:H thin films deposited by PECVD . In our previous work  we already discussed the a-Si:H-based DOS structure, together with its static characteristics. In this paper, taking advantage of optical and thermal coupled analyses, both for steady-state and dynamic regime, the performances of an optimized version of the DOS are analyzed. The aim of our design is to achieve good performances in terms of crosstalk and switching time, but with relaxed requirement for the realization and operation of device. We consider large cross section single mode waveguides, so that coupling losses should be reduced, and a relatively large angle between the output branches, in order to simplify the realization and reduce the length of device. In table 1, the performances of recent devices presented in literature and regarding some of the most promising approach (namely LiNbO3 and polymers) are reported together with data coming from the analysis of our proposal.
We note that LiNbO3 based DOS gets the best performances, but it presents also the worst drawback, the high fabrication cost. Polymers and amorphous silicon are both low cost, they get comparable cross talk, but in the case of polymer the further advantage is the low power dissipation and the drawback is the response time in the order of milliseconds. For amorphous silicon based devices we get a better response time, in the order of microseconds, but the power dissipation is higher.
2. Device design and working principle
The DOS is essentially a Y-junction consisting of an input tapered waveguide, which adiabatically adapts the launched fundamental mode in the bimodal DOS input region, followed by two monomode symmetric output branches (see schematic reported in Fig 1).
The proposed device consists of a rib waveguide made by an a-Si:H core deposited onto a Zinc Oxide (ZnO) lower cladding. The substrate is standard crystalline silicon (c-Si) wafer. Crystalline silicon wafers are used as substrates and due to its high thermal conductivity it can be treated as an almost perfect heat sink. ZnO lower cladding is preferred to standard SiO2 layer because its higher thermal conductivity, which ensures faster switching operation. Moreover, the high refractive index contrast between a-Si:H and ZnO allows the infrared radiation to be effectively confined in the low-loss a-Si:H guiding layer. The high confinement of the optical radiation into the a-Si:H film reduces also the scattering losses due to the ZnO-a:Si:H interface, which roughness can be estimated in a few tens of nm.
The operation of the device is based on the modal effective index variation induced by waveguide heating, which can modify the beam propagation pattern inside the structure itself. The heating can be induced by means of a metallic (Al) film, acting as resistor on the top of the waveguide, whose length is L heater, as reported in Fig. 1. At the splitter end two straight waveguide are considered and at their output section are evaluated optical fields and power. The basic principle of DOS is the following: at the branching point, equal amounts of light are launched into each waveguide. If the two branches have the same temperature, the DOS is geometrically symmetric and acts as a -3 dB power divider. Therefore, the optical input power has to be evenly divided into the two output ports. On increasing the temperature of one arm, the branch is made asymmetric; the light is guided by adiabatically evolving the input mode to the mode of destination arm with the increased refractive index. Accurately: at the branching point, where the gap between the waveguide is small, equal amounts of lights are launched into each monomode waveguide in phase to excite the local normal mode of the branching waveguides. At the end of the branching structure, where the waveguide gap is large, the most of power of the zero order normal mode is in the waveguide in which the refractive index has been enhanced. The field of the zero order mode changes its shape as it propagates along the branch structure; this effect is called modal evolution. Such ideal behavior is encountered as long as the geometric transition represented by the branch is sufficiently adiabatic, so that mode coupling between the normal modes does not occur. Switching between the two output arms is provided by reversing the sense of asymmetry. The division of the modal power over the two branches is related to the DOS angle, to the effective indices of the output branch and to the difference between the effective index and the index of the background.
3. Thermo-optical numerical analysis
The performances of thermo-optic DOS have been analyzed using commercial simulation software packages. OlympIOs — Integrated Optics Software produced by C2V  was employed in order to analyze the optical response and a finite element model has been implemented with the FEM code Ansys Inc. 9.0 , in order to analyze the thermal response.
3.1 Modal Analysis
In order to determine the values of propagation constant and modal pattern for each desired mode supported by the layered structure, we have utilized a Finite Difference (FD) mode solver [21, 22]. Because of the scalar approximation does not discriminate between transverse electric (TE) and transverse magnetic (TM) modes and being accurate only when the refractive index difference is small, the semi-vectorial finite difference (SVFD) method is one of the most commonly used techniques for analysis of integrated optical waveguides. However SVFD method does not provide a full vectorial description of the guiding process, but given that, in our device, the hybrid nature of the modes is negligible because of the small depth of etched part we have performed only a semi-vectorial calculations. Moreover, in order to calculate leaky modes in straight waveguides, perfectly matched layers (PML) as the boundary condition, have been implemented . In order to ensure single-mode propagation in large cross section semiconductor rib waveguides, the simultaneous fulfillment of following circumstances is required :
From our simulations, we get the following results (see Fig. 2): the bimodal channel waveguide has a rib width W=12 µm and a core layer thickness H=4 µm, the rib channel waveguide is defined by a lateral etch (H - h)=0.4 µm. The monomode waveguide arms are 6-µm-wide and the thickness of the ZnO layer is 0.3 micron. The aluminum stripe on the top of the waveguide is 0.5-µm-thick and 2.0-µm-wide.
We note that dimensions have been chosen in order to ensure single-or bimodal propagation in large cross section semiconductor rib waveguides according to the rules of ref. 24, and the etching depth is relatively small for two reasons. First of all we would like to ensure low coupling losses with standard optical fiber and a small lateral perturbation in the waveguide geometry induces 2D confinement with small asymmetry of the supported mode, better matching the Gaussian field of an optical fiber. Computation of overlap factor between the fundamental mode of input waveguide and a 2D Gaussian field with a 6-µm-wide waist returns a coupling loss of about 2.5 dB. Moreover, a weak lateral confinement in the single-mode waveguide prevents bi-modal behavior during heating.
From the optical point of view the aluminum stripes on the top of the rib waveguides induce high polarization depending losses (PDL). In fact, even if the silicon is an isotropic optical material and the waveguide geometry, with its aspect ratio W/H=1.5, should not induce high birefringence in the device, the introduction of the heater (a metallic electrode) in the vicinity of the propagating field will inevitably cause a strong effect on the TM polarization. A simple 2D modeling of the structure indicates a loss of about 16 dB/cm for TM mode, against ~1dB/cm for TE mode: this will produce a large PDL for quite long electrodes. In order to make the proposed device operation polarization insensitive, we imagine to introduce a thin SiO2 layer between the Al heater and the a-Si:H waveguide acting as upper cladding. We analyze, see Fig. 3, the effect of the introduction of this SiO2 upper cladding; if such layer is thick enough the Polarization Depending Losses (PDL) are about 0.1 dB. The only price to pay is the introduction of an additional small stiffness in the thermal model of the device, whose influence on the global operation is negligible. Therefore, we add in the device design a SiO2 buffer layer which is 0.2-µm-thick and 4-µm-wide (see Fig. 2).
The optical parameters, at λ=1550 nm, utilized in the modal and propagation analysis are reported in table 2 for all the materials . In particular, together with the refractive index n and the extinction coefficient k, is reported also the thermo-optic coefficient ∂n/∂T, which takes in account the refractive index variation induced by a temperature change.
3.2 Thermal analysis
In order to evaluate the thermal transient and steady-state response of the device, the equation to be solved is the heat transfer Fourier’s equation in transient condition with constant thermal conductivity:
where ρ is the material density, c p is the specific heat, k is the thermal conductivity, and Q(x,y,z,t) is the heat generation rate per unit volume. The initial condition is:
and the boundary conditions are:
Equation 4 establishes natural convection as the heat transfer mechanism between the device and the air; s in the surface outward normal, h is the natural heat transfer coefficient , T S is the surface temperature and T A in the air temperature. The boundary condition in Eq. (5) states that the lateral surfaces are adiabatic, while Eq. (6) assigns a fixed temperature at the bottom of the device, which can be realized, for instance, with a Peltier cell. In table 3 are reported the thermal parameters used in the simulations.
We have restricted the analysis to a 2-D model in the x-y plane normal to the light propagation direction; the active region cross section used for the thermal simulations is sketched in Fig. 4, together with the finite element mesh pattern.
Disregarding the minimal border effects at the waveguide ends, a temperature gradient exist only in this plane. This hypothesis can be justified by the observation that, for the proposed device, the aspect ratio between its length and the greater transverse dimension is about 15. The structure has been meshed with 8271 iso-parametric linear elements of the type Plane55, in which we have predicted the opportunity to suggest an expected device length . In fact, this element can also have a z-depth specified, that is, in this case, L heater=7500 µm (see Fig. 1). This is the distance from the initial part of the two-branch monomode waveguide to the end of the optimized active device. In practice, metallic heater is realized on the top of two single-mode branching arms starting from the point where the bimodal waveguide splits in two arms. In the following we will discuss the heater length optimization which provides us the useful metallic strip dimension for performances maximization. The FEM model 2D geometry reported in Fig. 4 is related to the section of splitter where the arms are at a distance of 10 µm, that is a length which ensures the thermal separation between the two DOS branches. From this point to the end of the device the power dissipated on one arm won’t diffuse to the other one, while in the initial region of the active splitter the heat generate will diffuse even to the cold side of the device. This fact doesn’t limit the DOS performances but, on the contrary, allows to reduce the crosstalk.
A detailed analysis was made to evaluate the accuracy and the operation speed of the thermo-optic device. The finite element modeling (FEM) software Ansys was used to perform a steady-state and a transient thermal analysis of the system to evaluate temperature profile, power dissipation and transient response of the DOS. For both static and dynamic analysis, our target is to reach an average heating ΔT between left and right waveguide surfaces of about 30 °C. This value, as discussed in the following, guarantees the total switching of optical power between two arms. In order to simulate a real heating element and reach the desired temperature increase, in the thermal model we added two aluminum region 0.5-µm-thick and 2.0-µm-wide on the top of both single-mode waveguide arms, shown as purple regions in the meshed structure reported in Fig. 4. A heat generation density Q(x,y,z,t)=115 µW/µm3 is simulated in the right metallic element, for a total power dissipation of about 862.5 mW. As reported in Fig. 5, these driving conditions are adequate to increase the waveguide surface temperature on the hot arm of 30 °C, while in the center of the guiding layer and below the electrode the temperature increase is about 14 °C (see the A-A′ 1D temperature profile reported in Fig. 5 and the CC′ profile in Fig. 6). Moreover, the temperature of the top metal electrode reaches about 55 °C.
Steady-state simulations confirm the good heat confinement in the guiding film under the biased electrode, thanks to a-Si:H thermal properties, namely its low thermal conductivity, which allows a strongly localized heating. For this reason, even in the simulated section where the distance between the splitter arms is only 10 µm or less, heat does not diffuse laterally so much, permitting a clean temperature difference between two branches, free from undesired thermal coupling. As expected, the low a-Si:H thermal conductivity causes high heating efficiency, low thermal crosstalk but, as tradeoff, high switching times, in particular with reference to the cooling phase. In our device the presence of c-Si substrate, characterized by very high thermal conductivity, reduces the cooling transient. The use of crystalline silicon was limited to the substrate because its thermal characteristics are not useful, in our opinion, for an efficient DOS operation. First of all, c-Si films are typically available onto silicon dioxide layers in the so-called SOI structures. The presence of the lower SiO2 layer, characterized by a very low thermal conductivity, would induce higher switching times for the device, because the high thermal RC constant introduced between active region and the bottom heat sink (thermal ground) by this layer in an equivalent electrical model. Moreover, the c-Si high thermal conductivity implies that the heat generated on the top of one DOS arm will diffuse laterally in the guiding film very quickly and with high uniformity, especially when an insulating layer is present at its bottom, so destroying the switching effect. Finally, SOI substrates are very expensive compared to the standard substrate we propose to use. For these reasons we believe that the structure proposed in the paper has some advantages respect to a c-Si based devices. As we can see in Fig. 7, where are reported the temperature dynamic in the center of the waveguide under the electrode and at the point just below the electrode in the middle of the electrode itself, heating and cooling times are in the order of 10 µs, value that makes the proposed device able to operate in the range of several tens of kilohertz.
3.3 BPM Analysis
In our analysis we have utilized a Finite-Difference BPM numerical schemes (FD-BPM) based on a fourth order Padè method, in order to accurately simulate propagations with wide angles of divergence . The optical behavior of the proposed DOS was verified by means of BPM simulations. After the analysis of the device operation in OFF-state, that is with no bias applied to the electrodes, we performed an optimization with respect to the variation of some geometrical parameters and of driving temperature. For each configuration we have exported the 2D static temperature distribution from Ansys simulator, converted this into a refractive index distribution and calculated from this the propagation pattern inside the structure. Some geometrical parameters were fixed, like the branching angle α=0.5 deg, or the waveguide transverse geometry, which was discussed in the previous section. The splitter aperture angle was fixed to a value greater the typical ones reported in literature (0.1–0.2 deg) in order to minimize the total length of the device but preserving the adiabatic transition between the bi-modal input waveguide and the two single-mode DOS arms. As expected, we observe that if the two branches have the same temperature the DOS is geometrically symmetric and acts as a power divider. Therefore, the optical input power has to be evenly divided into the two output ports and the global insertion loss, due to the non-ideal power splitting but excluding the fiber coupling losses, is about 0.4 dB.
The optimization was carried out changing simultaneously the hot arm waveguide surface temperature increase from 0 to 50 °C and the heater length from 1000 to 9000 µm. Of course, on increasing the temperature of one arm, the branch is made asymmetric; the light is guided by adiabatically evolving the input mode to the mode of destination arm with the increased refractive index. Figure 8 shows the results of this analysis. The global parameter evaluated for this optimization is the optical output crosstalk, defined as -10×log10(P off/P on), where P off and P on are, respectively the output optical power at the end of the cold and hot DOS arm. The higher value of the crosstalk, indicating that more optical power is switched to the ON arm, is about 25 dB and is obtained for a heater length of 7500 µm and an induced surface heating of 30 °C. For this geometry the output separation G between the DOS arms is about 65 µm. The corresponding fundamental mode effective index imbalance, induced between the output waveguides, is Δn eff=2.47×10-3. We note that the crosstalk decreases to about 20 dB if the temperature is further increased, because the biased DOS arm becomes bi-modal and, consequently, at its output, in conjunction with the passive collecting waveguide, an additional modal mismatch loss is present. In Fig. 9 is reported the output optical intensity pattern for the optimum configuration; it is clear the optical power switching from the OFF arm, on the left, to the ON arm, on the right. The shape of the switched optical filed is symmetrical and its intensity is, as expected, 25 dB above the residual field in the cold DOS branch.
Finally, we have also carried out the analysis of the DOS performances for different lengths L wgout of the output collecting waveguides shown in Fig. 1. The additional loss in the switched arm, which implies a reduction of the crosstalk, is due to the different modal distribution between the heated waveguide and the passive structure. Even if these waveguides are geometrically identical, except for the presence of the SiO2-Al upper cover, the slight relative angle and the perturbation induced by the heating change the modal characteristics.
For that reason the optical mode emerging from the DOS section needs to propagate for a few thousands of microns before to stabilize in the output external optical circuit. The additional reduction of the crosstalk can be evaluated from Fig. 10, and reaches its regime value of about -1.7 dB after 2000 µm from the DOS output section.
In this paper, due to its unique combination of properties, i.e. a high thermooptic coefficient and a low thermal conductivity, hydrogenated amorphous silicon (a-Si:H) has been considered as active material in order to design a DOS switch. The effort of our investigation was to achieve a good compromise between performances and relaxed requirement for realization and characterization of devices. In this investigation, a device having a length of 7500 µm and branching angle of 0.5 deg, presents optical switching with crosstalk of about 25 dB, and response time in the order of a ten microseconds. With respect to a polymers approach, the advantage is that a faster device of about 3 orders of magnitude is obtained; the drawback is the higher dissipated power. However, while the dissipation of power could be a technological issue and could be solved in some way, the response time is a limitation inherent to the physics of materials. Therefore, in all the cases in which response time of polymers are not acceptable for applications, our device could be a possible candidate.
References and Links
1. Y. Silberberg, P. Perlmutter, and J. E. Baran, “Digital optical switch,” Appl. Phys. Lett. 51, 1230–1232 (1978). [CrossRef]
2. W. K. Burns, “Normal mode analysis of waveguide devices. Part I: theory,” IEEE J. Lightwave Technol. 6, 1051–1057 (1988). [CrossRef]
3. W. K. Burns, “Normal mode analysis of waveguide devices. Part II: Device output and crosstalk: theory,” IEEE J. Lightwave Technol. 6, 1058–1068 (1988). [CrossRef]
4. G J. M. Krijnen, H. J. W. M. Hoekstra, P. V. Lambeck, and T. J. M. A. Pompa, “Simple analytical description of performance of Y junctions,” Electron. Lett. 28, 2072–2074 (1992). [CrossRef]
5. R. Krahenbuhl, M. M. Howerton, J. Dubinger, and A. S. Greenblatt, “Performance and modeling of advanced Ti:LiNbO3,” IEEE J. Lightwave Technol. 20, 92–99 (2002). [CrossRef]
6. R. Krahenbuhl, M. M. Howerton, J. Dubinger, A. S. Greenblatt, and S. T. Vohra, “Reflective digital optical switch (RDOS) for DWDM optical network applications,” IEEE Photon. Technol. Lett. 13, 34–36 (2001). [CrossRef]
7. E. L. Wooten, K. M. Kissa, A. Yi-Yan, E. J. Murphy, D. A. Lafaw, P. F. Hallemeier, D. Maack, D. V. Attanansio, D. J. Fritz, G. J. McBrien, and D. E. Bossi. “A review of Lithium niobate modulators for fiber optic communications systems,” IEEE J. Sel. Top. Quantum Electron. 6, 69–80 (2000). [CrossRef]
8. M. N. Kahn, B. I. Miller, E.C. Burrows, and C. A. Burrus, “High-speed digital Y-branch switch/modulator with integrated passive tapers for fibre pigatailing,” Electron. Lett. 35, 894–895 (1999). [CrossRef]
9. R. Moosburger and K. Petermann, “4×4 digital optical matrix switch using polymeric oversized rib waveguides,” IEEE Photon. Technol. Lett. 10, 684–686 (1998). [CrossRef]
11. C. Jang and R. T. Chen, “Polymer-Based 1x6 Thermo-optic switch incorporating an Elliptic TIR Waveguide Mirror,” IEEE J. Lightwave Technol. 21, 1053–1058 (2003). [CrossRef]
12. W. H. G. Horsthuis and M. B. J. Diemeer, “Components for multiple wavelength systems in polymer optoboard technology,” IEEE LEOS 8th Annual Meeting Conference Proceedings 2, San Francisco, USA, 251–252 (1995). [CrossRef]
13. N. Ooba, S. Toyoda, and T. Kurihara, “Low crosstalk and low loss 1×8 digital optical switch using silicone resin waveguides,” Electron. Lett. 35, 1364–1365 (1999). [CrossRef]
14. M. Hoffmann, P. Kopka, and E. Voges, “Thermooptical digital switch arrays in silica-on-silicon with defined zero-voltage state,” IEEE J. Lightwave Technol. 16, 395–400 (1998). [CrossRef]
15. O. B. Usev, A. M. Kuznetsov, E. I. Terukov, M. S. Bresler, V. K. Kudoyarova, I. N. Yassievich, B. P. ZaKharchenya, and W. Fuhs, “Room-temperature electroluminescence of erbium-doped amorphous hydrogenated silicon,” Appl. Phys. Lett. 70, 240–242 (1997). [CrossRef]
16. M. Okamura and S. Suzuki, “Infrared photodetection using a-Si:H photodiode,” IEEE Photon. Technol. Lett. 6, 412–414 (1994). [CrossRef]
17. G. Cocorullo, F. G. Della Corte, R. De Rosa, I. Rendina, A. Rubino, and E. Terzini, “Amorphous silicon based guided-wave passive and active devices for silicon integrated optoelectronics,” IEEE J. Sel. Top. Quantum Electron. 4, 997–1002 (1998). [CrossRef]
18. L. Sirleto, M. Iodice, F. G. Della Corte, and I. Rendina, “Digital optical switch based on amorphous silicon waveguide,” Opt. Eng. 42, 3417–3418 (2003). [CrossRef]
19. OlympiOs Integrated Optics Software-Manual, http://www.c2v.nl.
20. Ansys 9.0 Release Guide, 2005, http://www.ansys.com.
21. S. Sujecki, T. M. Benson, P. Sewell, and P. Kendall, “Novel vectorial analysis of optical waveguides,” IEEE J. Lightwave Technol. 16, 1329–1335 (1998). [CrossRef]
22. W. P. Huang, C. L. Xu, W. Lui, and K. Yokoyama, “The perfectly matched layer (PML) boundary condition for the beam propagation method,” IEEE Photon. Technol. Lett. 8, 649–651 (1996). [CrossRef]
23. R. Scarmozzino, A. Gopinath, R. Pregla, and S. Helfert, “Numerical techniques for modeling guided-wave photonic devices,” IEEE J. Sel. Top. Quantum Electron. 6, 150–161 (2000). [CrossRef]
24. S. P. Pogossian, L. Vescan, and A. Vonsovici, “The single mode condition for semiconductor rib waveguides with large cross section,” IEEE J. Lightwave Technol. 16, 1851–1853 (1998). [CrossRef]
25. D. Y. Smith, E. Shiles, and M. InokutiE. D. Palik, in Handbook of Optical Constants of Solids, ed. (Academic Press, San Francisco, CA, 1998).
26. J. C. Sturm, W. Wilson, and M. Iodice, “Thermal effects and scaling in organic light emitting flat panel displays,” IEEE J. Sel. Top. Quantum Electron. 4, 75–82 (1998). [CrossRef]
27. I. Ilic, R. Scarmozzino, and R. M. Osgood, “Investigation of the Pade approximant-based wide-angle beam propagation method for accurate modeling of waveguiding circuits,” IEEE J. Lightwave Technol. 14, 2813–2822 (1996). [CrossRef]