# Oblique Du-Fort Frankel Beam Propagation Method

Authors: Ken Chan; Philip Sewell; Ana Vukovic; Trevor Benson
Journal: Advances in OptoElectronics
 (2011)

---

## Abstract

The oblique BPM based on the Du-Fort Frankel method is presented. The paper demonstrates the accuracy and the computational improvements of the scheme compared to the oblique BPM based on Crank-Nicholson (CN) scheme.


---

## Body

## 1. Introduction

Increasingly complex optical devices demand computationally fast and memory efficient algorithms for modelling purposes. Finite difference beam propagation method (FD-BPM) is a popular numerical technique for simulating large network of optical components due to its computational advantages over classical numerical techniques such as Finite Difference Time Domain (FDTD) method. The BPM method is commonly applied in the Cartesian coordinate system. However when the boundaries of an optical component are not aligned to the Cartesian mesh, for example in the case of tilted waveguides, bends and Mach-Zehnder modulators, sampling on the Cartesian mesh introduces nonphysical staircasing noise. The noise can be minimised by using very fine mesh but that in return incurs large computational costs. To more efficiently reduce the sampling error an improved three-point formulas are used at the interface which take into account the distance between the boundary and the transverse sampling points [1–3]. Further increase in accuracy of the Cartesian BPM, particularly for strongly guided waveguides, is achieved by considering the longitudinal component of magnetic field which is commonly neglected in the standard FD-BPM method [4]. In contrast to Cartesian system, Oblique and Structure Related (SR) coordinate system offers an accurate and efficient alternative for modelling nonorthogonal structures and automatically satisfies ∂n/∂z=0. The sampling grid of the SR mesh is aligned with the component material boundary thus eliminating staircase error and allowing relaxation in mesh size. Various SR-BPM schemes have been introduced [5–11] and different schemes can be combined together to map out the optical component. Furthermore, the SR coordinate system ensures high accuracy for the simple paraxial BPM formulation even without the use of wide-angled schemes [10]. The oblique equation takes into account the propagation direction, which is usually parallel to the structure boundary. Hence the mode-mismatch error is small. One of the motivations for implementation in oblique coordinates is to remove the need for high-order wide-angle scheme which requires substantial computational resources. Wide-angle for oblique coordinate has been developed by Sujecki [11].  However the author has also confirmed that the wide-angle oblique approach should in principle only be applied to low refractive index contrast structures [12].So far SR and oblique BPM schemes have been implemented using implicit scheme such as CN. Whilst for the two-dimensional (2D) structures this is computationally fast, in the case of modelling three-dimensional (3D) structures the CN scheme uses iterative matrix solvers such as BI-CGSTAB or GMRES [13] and thus requires huge computational resources. More computationally efficient Alternate Direction Iterative (ADI) schemes [14] cannot be implemented in SR coordinates due to mixed derivatives in transverse directions and are limited to Cartesian meshes. Alternatively, Du-Fort Frankel (DFF) schemes provide larger step size and better stability condition than simple explicit schemes and can be implemented on a parallel computational platform thus providing computational efficiency for modelling realistic 3D optical components. The downside of the DFF scheme is the inherent spurious or “ghost” mode that can affect the accuracy and stability of the scheme and which can be alleviated by the right choice of parameters and initial fields [15, 16].In this paper, the oblique BPM method is implemented using the DFF scheme. Section2 outlines the formulation of the method and Section 3 presents the results for the power loss and computational efficiency of the scheme and compares it against the oblique CN BPM scheme. The results are presented for tilted 2D and 3D waveguides and scalar fields.


## 2. Formulation

A general approach for formulating oblique BPM is outlined in [5] and is limited to structures that do not vary with the propagation direction and the β reference is parallel to the propagating direction.In this section an oblique BPM method based on the paraxial approximation and implemented in the DFF algorithm is outlined. Figure1 shows schematic presentation of two coordinate systems u,y,z′, and x,y,z, where z and z′ form an angle θ. It is assumed that the fields propagate with respect to the z′ axis. The 3D oblique coordinate BPM equation adopted from [9] is used as(1)(∂2∂x2+∂2∂y2+∂2∂z2+k2(x,y,z))Ψ=0,
where Ψ represents the scalar field of the form(2)Ψ(u,z′)=φ(u,z′)exp(-jβ0z′).Figure 1
Oblique coordinate system.The relationship between the oblique and the Cartesian coordinate systems is given as(3)x=u+z′sinθ,z=z′cosθ.The difference equations are obtained using the chain rule as(4)∂Ψ∂u=∂Ψ∂x,(5)∂Ψ∂z′=∂Ψ∂xsinθ+∂Ψ∂zcosθ,(6)∂2Ψ∂u2=∂2Ψ∂x2,(7)∂2Ψ∂u∂z′=∂2Ψ∂x2sinθ+∂2Ψ∂x∂zcosθ,(8)∂2Ψ∂z′2=∂2Ψ∂x2sin2θ+2∂2Ψ∂x∂zsinθcosθ+∂2Ψ∂z2cos2θ.Rearranging (6)–(8) and substituting into (1) results in the oblique wave equation:(9)(∂2Ψ∂z′2-2∂2Ψ∂u∂z′sinθ+∂2Ψ∂u2+∂2∂y2)+k2cosθΨ=0.Substituting the field from (2) into (9) gives the scalar oblique BPM equation as(10)2jβ0∂φ∂z′+sinθ∂2φ∂u∂z′=2jβ0sinθ∂φ∂u+∂2φ∂u2+∂2∂y2+(k2cosθ-β02)φ.This equation can be straightforwardly implemented using the CN scheme. In the case of the 2D modelling where the∂2/∂2y= 0, CN scheme requires tridiagonal matrix solver known as Thomas algorithm is used to solve (10) which is much faster than the sparse matrix solver. However in the 3D modelling, term ∂2/∂2y in (10) introduces two additional unknown field points in each calculation step thus resulting in five unknown field points. In CN scheme this requires sparse matrix solver such as the commonly used BI-CGSTAB iterative solver [17].In order to implement DFF scheme the local field pointφl,m,n and the transverse second derivatives are discretized as
(11)φl,m,n≈φl+1,m,n+φl-1,m,n2,∂2φl,m,n∂u2≈φl,m+1,n+φl,m-1,n-(φl+1,m,n+φl-1,m,n)Δu2,
where l, m, and n are the discretised position in the z′, u, and y direction. Substituting the DFF scheme into (10), the numerical implementation for oblique DFF BPM is as follows:(12)L2(φl+1,m,n-φl-1,m,nΔz′-sinθφl,m+1,n-φl,m-1,n2Δu)+sinθ(φl+1,m+1,n-φl-1,m+1,n2ΔuΔz′-φl+1,m-1,n+φl-1,m-1,n2ΔuΔz′)=φu′′+φy′′+(k2-β02)(φl+1,m,n+φl-1,m,n)2,
where L=2jβ0, and φa′′=(φm,a+1+φm,a-1-φm+1,a-φm-1,a)/Δa2.The additional term in (10) diminishes the explicit nature of the DFF, the 3D oblique DFF scheme corresponds to separable tri-diagonal matrices on each layer of the 3D structure. This can be efficiently solved by Thomas algorithm and is algorithmically efficient for parallel computing. Solving of the tri-diagonal matrix is more computationally efficient than solving the sparse matrix thus ensuring better computational efficiency of the 3D oblique DFF-BPM compared to the oblique CN-BPM.


## 3. Results

In this section  the accuracy and computational stability of the oblique DFF-BPM method is analyzed and compared with the oblique CN-BPM method. Both 2D and 3D tilted waveguides are analyzed.In an oblique coordinate system a tilted waveguide is essentially a straight waveguide and ideally the power transmitted in a tilted waveguide in the oblique system is identical to the power transmitted in a straight waveguide in the Cartesian mesh.  Figure2 compares power loss of the 2D tilted waveguide analyzed using the oblique DFF-BPM with that analyzed using the DFF scheme in Cartesian mesh for different tilted angles and different sampling meshes. Perfectly matched layer (PML) is used for absorbing any leakage from the simulation window. The implementation of PML for the DFF method is described in [16]. The PML layer is set as 1.0 μm and strength, σ/ω is set as 10. Guided mode at the wavelength of 1.15 μm is launched in the 1 μm wide slab waveguide with core refractive index n=1.1 surrounded by air. Longitudinal sampling Δz′ is fixed at 0.1 μm and the waveguide length is 409.6 μm. It can be seen that by reducing the sampling of the mesh the power loss in the Cartesian BPM is reduced whilst for the oblique DFF-BPM the change of the mesh size does not significantly change the power loss. This indicates that the mesh size can be more relaxed in the case of the oblique DFF-BPM ensuring faster run time.Figure 2
Power loss for a range of tilted angles for the oblique DFF-BPM and the Cartesian DFF-BPM method.Figure3 compares the computational runtime of the 3D oblique and Cartesian DFF- and CN-BPM methods for different number of total mesh points N. A simple 1 μm by 1 μm square metal-air waveguide is chosen for a fair speed test. Longitudinal step Δz′ is fixed at 0.05 μm, Δu=Δy, and the waveguide is 100 μm long. For the oblique simulation, the same waveguide is used but tilted at an angle of 5°. Figure 3 shows that the Cartesian DFF excels in speed even without any parallelization involved. Figure 3 also shows that oblique DFF-BPM is substantially slower than Cartesian DFF-BPM due to the implementation of the Thomas algorithm. However, when compared with the Cartesian and oblique CN-BPM methods, the oblique DFF-BPM method is much faster, especially for large mesh sizes. However, it should be noted that DFF-BPM requires smaller longitudinal step size than CN-BPM to achieve the same level of accuracy and maintain stability. The 3D oblique DFF-BPM is also suitable for parallel computing platform allowing for more computationally efficient simulations.Figure 3
Comparison of the computational time for the oblique and Cartesian DFF-BPM and CN-BPM w.r.t. the total number of mesh pointsN.Figure4 analyses the stability of the oblique DFF-BPM method. It is well known that the main weakness of the DFF algorithm is the appearance of the spurious solution [15]. The position of the spurious mode can be controlled by appropriate choice of the mesh size and the excitation [16]. It is desirable that the spurious mode is not too close to the waveguide true mode so that the stability of the algorithm is not affected. Generally increasing the transverse mesh size and reducing the propagation size step will place the spurious mode further away from the true mode [16]. Figure 4 investigates the impact of the tilted angle on the position of the spurious mode for a fixed mesh size. Figure 4 gives the Fourier transform of the field overlap between the incident field and the field evolved along the waveguide for different tilted angles. The 3D waveguide is a rectangular metal-air waveguide with dimensions of 1 μm by 0.5 μm and 204.8 μm long.  Transverse mesh size is Δu=Δy=0.1μm and longitudinal step sized is Δz′=0.05μm. Half sine wave of 1.15 μm wavelength is launched at the input. Figure 4 shows that the increase of the tilted angle brings the spurious and the true mode closer together. This will have implications on the maximum tilted angle that can be considered using the oblique DFF-BPM method.Figure 4
Fourier transform of field overlaps along the waveguide to identify the true mode (neff=0.82) and the appearance of spurious modes at various tilted angles.The stability of oblique DFF is determined by the oblique angle and the mesh size. The effect of the mesh size is examined in Figure5 for Δu=0.01μm and different tilted angles. The waveguide parameters are as in Figure 2. It shows that for small transverse mesh size, it is necessary to keep the propagation step size small to maintain stability. It can be seen that when the oblique angle is small, the oblique DFF behaves similar to the Cartesian DFF. As the oblique angle increases, the oblique DFF requires smaller propagation step or larger transverse mesh size to maintain the stability.Figure 5
Power gain for oblique DFF. It shows the necessity for keeping the propagation step size small to maintain stability when transverse mesh size is small (Δu=0.01μm).Figure6 examines the parameter choice and instability on the calculation of the effective index. The waveguide parameters are the same as in Figure 2. The obtained effective index is plotted for Δz′=0.025μm with various mesh sizes and compared between the Cartesian DFF applied to the straight waveguides and the oblique DFF applied to the 10° tilted waveguide. Figure 6 shows the Cartesian DFF and the 10° tilted oblique DFF agree very well for small mesh sizes but have significant discrepancy when the transverse mesh size is increased. It can be concluded that the stability condition has restrained the use of very small transverse mesh size in DFF. However, using large transverse mesh size would risk losing accuracy. However, it should be noted that the difference in coordinate system makes it difficult to compare results directly. A slice in the oblique coordinate is equivalent to a diagonal cross-section through multiple slices on a Cartesian coordinate system.Figure 6
Effective index of a slab waveguide using Cartesian DFF and 10 degree tilted oblique DFF method. Accuracy deteriorates as it reaches the unstable region.Figures7(a)–7(c) shows the field profiles of the 2D tilted waveguide obtained using the oblique and Cartesian DFF-BPM method. The slab waveguide is as in Figure 2 but 10 μm long and tilted by 10°. Figure 7(a) shows the field profile of the 10° tilted waveguide mapped in the Cartesian mesh and Figure 7(b) shows the same waveguide analysed using the oblique DFF-BPM method. It can be seen that the Cartesian mesh introduces large staircasing noise even for very fine sampling mesh. For comparison, the field profile of a standard nontilted waveguide modelled on the Cartesian mesh in Figure 7(c) behaves similarly to the tilted waveguide modelled on oblique coordinate. Perfectly Matched Layer (PML) boundary condition is used to absorb the wave leakage from the waveguide and the nonphysical staircase noise. The slight leakage of the guided mode near the input of the waveguide is due to the mismatch of the analytical mode used to excite the waveguide and the actual numerical field solution. This leakage disappears as the propagating mode is settled in the waveguide.(a) Field surface plot for10° tilted waveguide on Cartesian mesh. Scale has been capped to show the stair case noise. (b) Field surface plot for 10° tilted waveguide on oblique coordinate. (c) Field surface plot for straight waveguide on Cartesian mesh.
(a)(b)(c)


## 4. Conclusion

The paper demonstrates the implementation of the DFF algorithm into the 3D scalar oblique BPM method. The accuracy and stability of the oblique DFF is investigated on the 2D and 3D tilted waveguides and compared against the oblique CN-BPM method. The resulting method is computationally faster than oblique CN-BPM method and is suitable for parallel computing for further computational savings. It is noted that the appearance of the spurious mode can potentially limit the application of the DFF to very large titled angles.



---

*Source: 196707-2010-09-16.xml*