## Abstract

The development of three dimensional (3-D) waveguide structures for chip scale planar lightwave circuits (PLCs) is hampered by the lack of effective 3-D wide-angle (WA) beam propagation methods (BPMs). We present a simple 3-D wide-angle beam propagation method (WA-BPM) using Hoekstra’s scheme along with a new 3-D wave equation splitting method. The applicability, accuracy and effectiveness of our method are demonstrated by applying it to simulations of wide-angle beam propagation and comparing them with analytical solutions.

©2006 Optical Society of America

## 1. Introduction

The finite-difference beam propagation method (FD-BPM) is an effective and widely used technique for the study of optical waveguide devices [1–4]. It is based on a numerical solution of the scalar Helmholtz equation using the Crank-Nicholson scheme [5] which is unconditionally stable. In addition, the FD-BPM makes use of the highly efficient Thomas algorithm [5], the transparent boundary condition (TBC) [6], which efficiently suppresses the reflections at the edges of the finite computational window and can be easily implemented, the alternating-direction implicit finite difference method (ADIFDM) [5] which can greatly simplify the simulations for 3-D structures by splitting 3-D equation into two steps of 2-D simulation, and the slowly varying envelope approximation (SVEA). The SVEA substantially simplifies the algorithm by neglecting the second-order derivative with respect to *z* in the wave equation. However, it limits the simulations to paraxial beams along the *z* axis and low refractive index contrast ratio between the core and cladding of the waveguide [3, 7–9]. Fortunately, several schemes have been proposed to relax these limitations for wide-angle (WA) simulations, such as the method of lines BPM [10], the Lanczos reduction method [11], the Eigenmode expansion BPM [12] and the matrix expansion BPM [12], the power series expansion method [13], the collocation method [14], the BPM based on alternating direction implicit preconditioner [15], the Padé approximant scheme and the multistep method [16, 17], and the BPM using Hoekstra’s scheme [18]. Each of them has some advantages over the others, but also has drawbacks. For example, the Lanczos method has convergence problems, the method of lines requires that Eigenfunction expansions be calculated for each different cross section structure, the Padé approximant scheme increases the matrix bandwidth in higher-order Padé operators, etc. All these published methods can be applied to 2-D wide-angle simulations, but only a few have been proposed explicitly for 3-D structures. Even in commercialized software packages, such as BeamProp^{™}, 3-D wide-angle calculation can only be performed using up to Padé (1,1) approximant, and for more accurate calculation, only smaller step sizes is suggested. Because of the utility of 3-D waveguide structures, especially in the development of chip scale planar lightwave circuits, effective 3-D algorithms need to be developed. In this article, we propose a simple and practical 3-D BPM using Hoekstra’s scheme along with a new 3-D wave equation splitting method.

## 2. Formulation

The 3-D scalar Helmholtz equation is given by

where *k*
_{0} denotes the free-space wavenumber, *E*(*x*, *y*, *z*)exp(*iωt*) is the electric field component with angular frequency *ω*, and *n*(*x*, *y*, *z*) is the refractive index profile. If we assume the slowly varying envelope approximation (SEVA), then *E*(*x*, *y*, *z*) can be separated into two parts: the complex field amplitude *ψ*(*x*, *y*, *z*) (the axially slowly varying envelope term) and a propagation factor *exp*(-*ik*
_{0}
*n*
_{0}
*z*) (the rapidly varying phase term). The field is then expressed as

where *n*
_{0}
is the refractive index in the cladding. Substituting Eq. (2) into Eq. (1) yields

With *a*=2*k*
_{0}
*n*
_{0} and *b*=${k}_{0}^{2}$(*n*
^{2}-${n}_{0}^{2}$). In the paraxial approximation one assumes |*ia*∂*ψ*/∂*z*≫∂^{2}
*ψ*/∂*z*
^{2}| and so neglects the second order derivative |∂^{2}
*ψ*/∂*z*
^{2}| in Eq. (3), yielding

The FD-BPM calculates how an input optical field distribution propagates through a waveguide by first discretizing the field into a grid in the plane perpendicular to the z (propagation) direction. The field is then numerically propagated along the z-direction to the next section. This enables the wave equation to be solved numerically for waveguide structures which cannot be solved analytically. The electric field at the grid point of *x*=*mΔx*, *y*=*jΔy* and *z*=*lΔz* is expressed by *ψ*(*x*,*y*,*z*)=${\psi}_{m,j}^{\mathit{1}}$. First, we discretize Eq. (3) in the *z* direction as

In order to split the equation to the stage where the tridiagonal matrix equations may be formed, we apply the approximation *ψ*
^{l+1}+*ψ*
^{l}
=2*ψ*
^{l+1/2}, where *ψ*
^{l+1/2} is the intermediate field between the steps of *l* and *l*+1 in the propagation direction *z*. Eq. (5) can then be rewritten as

$$=\frac{{\partial}^{2}}{\partial {x}^{2}}\left(\frac{2{\psi}^{l+1\u20442}}{2}\right)+\frac{{\partial}^{2}}{\partial {y}^{2}}\left(\frac{{\psi}^{l+1}+{\psi}^{l}}{2}\right)+\frac{{\left(b\psi \right)}^{l+1}+2{\left(b\psi \right)}^{l+1\u20442}+{\left(b\psi \right)}^{l}}{4}$$

Using the alternating-direction implicit finite difference method (ADIFDM) [5], this can be further split into two equations:

In the 2-D algorithm [18, 19], ∂*ψ*/∂*z* can be replaced by the expression obtained from the paraxial wave equation. However, in the 3-D algorithm, if ∂*ψ*/∂*z* in Eqs. (7a) and (7b) is directly replaced by the paraxial wave Eq. (4), both *x* and *y* will be present in all steps at z direction, which results in the equations becoming intractable with the efficient Thomas method. Therefore, the paraxial equation needs to be separated. Equation. (4) can be rewritten as

Using the same approximation for the half-step field (*ψ*
^{l+1}+*ψ*
^{l}
=2*ψ*
^{l+1/2}), Eq. (8) becomes

$$+\frac{{\partial}^{2}}{\partial {y}^{2}}\left(\frac{{\psi}^{l+1}+{\psi}^{l}}{2}\right)+\frac{{\left(b\psi \right)}^{l+1}+2{\left(b\psi \right)}^{l+1\u20442}+{\left(b\psi \right)}^{l}}{4}$$

Again using the ADIFDM, Eq. (9) can be split into two equations

Finally, Eq. (10) is separated into the following three equations

Eqs. (11) were also derived by J. Shibayama in a different way and successfully applied to improve the accuracy of BPM [20]. The physical significance is that the effects from the second order derivatives in x and y are separated into two distinct half-step propagations. Replacing Eq. (7) by Eq. (11), we get the following two wide-angle equations

where

${P}_{\xi ,m,j}^{\pm}=\frac{1}{\Delta {\xi}^{2}}\left(\mp 1+\frac{4i}{a\Delta z}\right)$ and ${Q}_{\xi ,m,j}^{\pm}=\pm \frac{2}{\Delta {\xi}^{2}}+\frac{2\mathit{ia}}{\Delta z}-\frac{8i}{a\Delta z\Delta {\xi}^{2}}\pm \left(\frac{2i}{a\Delta z}+\frac{1}{2}\right){b}_{m,j}^{h},$

with ξ=*x* or *y* and *h*=*l* for *Q*^{-}
_{y}
, *l*+*1/2* for *Q*
_{x}
, and *l*+*1* for *Q*^{+}
_{y}
, respectively. In this algorithm, the second-order derivative with respect to *z* is shown as the difference between two first-order derivatives with respect to *z*, which are finally replaced by the second-order derivatives in transverse directions. As a result, the truncation error of this algorithm is maintained at *O*(*Δx*
^{2})+*O*(*Δy*
^{2}
) in transverse directions, This, however, is not ensured at the optical interfaces, as is in the classical BPM. The electric field at the intermediate step *l*+*1/2* is first calculated by using Eq. (12a), and then the field at *l*+*1* can be obtained by using the results of Eq. (12a) and Eq. (12b). Since Eqs. (12a) and (12b) are tridiagonal, they can be solved efficiently using the Thomas method, combined with Transparent Boundary Conditions (TBC). Higher accuracy is expected because the second-order derivative with respect to *z* is included in the new scheme.

## 3. Simulations and discussions

In order to show the applicability, accuracy and effectiveness of our method, we first compare the simulations of a 3-D Gaussian beam propagating in free space (unity refractive index) with a 30-degree tilt with respect to *z* axis and along a direction mid way between the x and y axes. Both our technique and the classical (paraxial) BPM are compared with the analytical results. The free space wavelength for this simulation is 0.85 µm, and the input Gaussian beam profile has a waist radius *w*
_{0}=3.0 µm at (*x*, *y*)=(0, 0). The simulations are performed in a window of *x*=[-20 µm, 55 µm] and *y*=[-20 µm, 55 µm], with the step sizes of *Δx*=*Δy*=0.15 µm and *Δz*=0.2 µm. Figure 1 shows the input at *z*=0 and the outputs at *z*=60 µm in the same simulation window, in order to facilitate the comparison of their relative positions and moduli. As is shown, our result (in green) has the closer modulus to the exact solution (in blue), with a relative error in L^{2} norm of 1.44%, but with a slight shift of 3.66% in position, while the paraxial BPM calculation (in red) propagates at the wrong angle and incorrectly preserves the symmetric Gaussian shape, rather than the spread out, asymmetric profile. In our method, the paraxial approximation is still implicit at each *z* step of the simulation, even though the second-order derivative with respect to *z* is included in the derivation, which causes the slight mismatch between our result and the exact solution. It would

We also analyze the propagation in a single-mode channel waveguide, again tilted at 30-degrees with respect to *z* axis and 45 degrees with respect to the *x* and *y* axes. The dimensions of the channel are 4 µm by 4 µm, and the refractive index of the core and the cladding is 1.55 and 1.546 at the wavelength of 0.85 µm, respectively. The simulations are performed using the same step sizes and the same simulation window as in the case of the Gaussian beam propagation above. Figure 2 shows the input at *z*=0 and the outputs at *z*=60 µm in the same simulation window. Our result (in green) has a similar profile to the analytical calculation (in blue) with a relative error in L^{2} norm of 5.21%, with again a slight shift of 4.91% in position. The red profile is the result of the paraxial BPM. Similar to the case with the free space propagation, it is not propagating at the right direction and its profile is spreading out more quickly even inside a waveguide channel.

Simulations were performed using different *Δx* (*Δy* is set to equal *Δx*) and *Δz*. For a fixed *Δx* in the examples above, changes in *Δz* in the range from 0.05 µm to 0.25µm, only result in slight or even indistinguishable differences in output. Figure 3 shows the dependences of the relative L2 norm errors on the grid spacing, *Δx* (*Δy*), of the output field moduli and the relative position shift of the output field profile, with respect to the exact solutions, at *z*=60 µm of the Gaussian beam propagation and mode propagation in the waveguide. The simulation calculated using our wide-angle BPM has much better relative L^{2} norm errors and relative position shifts than the classical BPM for the Gaussian beam propagation in uniform media. It also has better results for the waveguide beam propagation. The errors calculated for the waveguide beam propagation using the classical BPM change irregularly because the classical BPM does not properly account for the wide-angle term (the second-order derivative with respect to *z*), which causes an asymmetrical output, rather than the symmetric analytical profile. The truncation error in this case is overwhelmed and the L^{2} norm error no longer varies in a predictable fashion. The plots also suggest small transverse step sizes are preferred. Both the errors and position shifts obtained using our wide-angle BPM for the beam propagation in the waveguide are higher than for the Gaussian propagation in uniform media, and become worse for higher refractive index contrasts between the core and cladding. This is because the truncation error *O*(*Δx*
^{2})+*O*(*Δy*
^{2}
) of the algorithm is not ensured at the transverse discontinuities. Our algorithm can be directly applied to a semi-vectorial scheme to overcome this problem by starting with the semi-vectorial Fresnel wave equation following a similar derivation as above. The extension to a full-vectorial treatment is not trivial and is currently under investigation.

A comparison with other variations of the BPM, such as the examples simulated in Ref. [15] and references therein, is difficult, because these examples have z-variant structures which cannot be accurately calculated only by using wide-angle BPMs, including our method. It is commonly known that classical BPMs and wide-angle BPMs cannot suppress the numerical reflections and the staircasing approximation for z-variant structures such as tapers and Y-branches, which require special BPM algorithms [21–23]. These algorithms typically involve nonorthogonal coordinate systems. It is likely that a combined use of these structure related algorithms and wide-angle algorithms would improve the accuracy even for wide-angle z-variant structures.

The wide-angle BPMs based on Padé approximant operators and the multistep method developed by G. R. Hadley [16, 17] are the most commonly employed techniques to improve the numerical accuracy for wide-angle simulations of 2-D structures. However, when they are applied to 3-D structures, large sparse matrices will be involved. The matrix algebraic equation is no longer tridiagonal and needs to be solved using iterative methods, such as the bi-conjugate gradients method [5]. Iterative methods have convergence problems for ill-conditioned sparse matrices and cannot be applied to very wide-angle structures, which usually result in ill-conditioned sparse matrices.

In summary, we derived a simple and practical 3-D wide-angle BPM, and demonstrated its utility for two cases: the wide-angle propagation of a Gaussian beam with a 30 degree tilt in free space, and the wide-angle propagation in a single-mode waveguide channel 30 degree tilted with respect to *z* axis. Our results have much better matches in both position and modulus with the exact solutions than the paraxial calculation. It can be easily implemented because its tridiagonal form is same as the classical 3-D paraxial BPM except for the coefficients.

This material is based upon work supported by the National Science Foundation under Grant No. 0348955. The authors also express gratitude to Professor Amy Liu in the Department of Physics at Georgetown University, for her helpful discussion.

## References and links

**1. **M. D. Feit and J. A. Fleck Jr., , “Light propagation in graded-index optical fibers,” Appl. Opt. **17**, 3990–3998 (1978). [CrossRef] [PubMed]

**2. **D. Yevick, “A guide to electric field propagation techniques for guided-wave optics,” Opt. Quantum Electron. **26**, S185–S197 (1994). [CrossRef]

**3. **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]

**4. **J. Van Roey, J. van der Donk, and P. E. Lagasse, “Beam-propagation method: analysis and assessment,” J. Opt. Soc. Am. **71**, 803–810 (1981). [CrossRef]

**5. **W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vettering, *Numerical recipes: The art of scientific computing*, (Cambridge University Press, New York, 1986).

**6. **G. R. Hadley, “Transparent boundary condition for beam propagation,” Opt. Lett. **16**, 624–626 (1991). [CrossRef] [PubMed]

**7. **Y. Chung and N. Dagli, “Assessment of finite difference beam propagation,” Opt. Lett. **26**, 1335–1339 (1990).

**8. **M. D. Feit and J. A. Fleck Jr.,, “Analysis of rib waveguides and couplers by the propagation method,” J. Opt. Soc. Am. A **7**, 73–79 (1990). [CrossRef]

**9. **C. Vassallo, “Reformulation for the beam-propagation method,” J. Opt. Soc. Am. A **10**, 10 2208–2216 (1993). [CrossRef]

**10. **J. Gerdes and R. Pregla, “Beam-propagation algorithm based on the method of lines,” J. Opt. Soc. Amer. A **8**, 389–394 (1991). [CrossRef]

**11. **R. P. Ratowsky and J. A. Fleck Jr., , “Accurate numerical solution of the Helmholtz equation by iterative Lanczos reduction,” Opt. Lett. **16**, 787–789 (1991). [CrossRef] [PubMed]

**12. **P.-C. Lee, D. Schulz, and E. Voges, “Three-dimensional finite difference beam propagation algorithms for photonic devices,” J. Lightwave Technol. **10**, 1832–1838 (1992). [CrossRef]

**13. **P.-C. Lee and E. Voges, “Three-dimensional semi-vectorial wave-angle beam propagation method,” J. Lightwave Technol. **12**, 215–225 (1994). [CrossRef]

**14. **A. Sharma and A. Agrawal, “New method for nonparaxial beam propagation,” J. Opt. Soc. Am. A **21**, 1082–1087 (2004). [CrossRef]

**15. **S. L. Chui and Y. Y. Lu, “Wide-angle full-vector beam propagation method based on an alternating direction implicit preconditioner,” J. Opt. Soc. Am. A **21**, 420–425 (2004). [CrossRef]

**16. **G. R. Hadley, “Wide-angle beam propagation using Padé approximant operators,” Opt. Lett. **17**, 1426–1428 (1992). [CrossRef] [PubMed]

**17. **G. R. Hadley, “Multistep method for wide-angle beam propagation,” Opt. Lett. **17**, 1743–1745 (1992). [CrossRef] [PubMed]

**18. **H. J. W. M. Hoekstra, G. J. M. Krijnen, and P. V. Lambeck, “On the accuracy of the finite difference method for applications in beam propagating techniques,” Opt. Commun. **94**, 506–508 (1992). [CrossRef]

**19. **Z. Ju, J. Fu, and E. Feng, “A simple wide-angle beam-propagation method for integrated optics,” Microwave Opt. Technol. Lett. **14**, 345–347 (1997). [CrossRef]

**20. **J. Yamauchi, J. Shibayama, O. Saito, O. Uchiyama, and H. Nakano, “Improved finite-difference beam-propagation method based on the generalized Douglas scheme and its application to semivectorial analysis,” J. Lightwave Technol. **14**, 2401–2406 (1996). [CrossRef]

**21. **S. Sujecki, P. Sewell, T. M. Benson, and P. C. Kendall, “Novel beam propagation algorithms for tapered optical structures,” J. Lightwave Technol. **17**, 2379–2388 (1999). [CrossRef]

**22. **Y. P. Chiou and H. C. Chang, “Beam-propagation method for analysis of z-dependent structures that uses a local oblique coordinate system,” Opt. Lett. **23**, 998–1000 (1998). [CrossRef]

**23. **D. Z. Djurdjevic and T. M. Benson, et. al, “Fast and accurate analysis of 3-D curved optical waveguide couplers,” J. Lightwave Technol. **22**, 2333–2340 (2004). [CrossRef]