1. Introduction

The classical finite-difference beam propagation method (FD-BPM) [1-5] is a technique for calculating the electromagnetic field transmitted through materials with inhomogeneous refractive index profiles. The electric field is discretized on a planar grid perpendicular to the direction of propagation. The field on the following plane is calculated based on a numerical solution of the scalar Helmholtz equation using the Crank-Nicholson scheme [6], and solved using the highly efficient Thomas algorithm [6] along with transparent boundary conditions (TBC) [7]. BPM simulations have been used extensively to simulate propagation in 2-D, but 3-D methods are more difficult to implement. 3-D algorithms for the paraxial case have been simplified using the alternating-direction implicit (ADI) method [6] by splitting the 3-D equation into two 2-D steps. Besides the limitation to paraxial beams, these methods also restrict the simulations to a low refractive index contrast ratio between the core and cladding of the waveguide [4-5, 8-9]. Efforts have been taken to relax these limitations for wide-angle (WA) simulations [10-19], of which the 2-D wide-angle BPM based on Hoekstra’s scheme [18, 19] is simple and practical, because the wide-angle term is included in the formulation by simply replacing the second-order derivative using the derivative of the first-order derivatives. Hadley has also developed WA-BPMs using the Padé approximant operators and the multistep method, which are widely used for 2-D structures, which, however, the author has also claimed to be restricted for 3-D structures due to the lack of efficient solvers [16, 17]. A 3-D wide-angle BPM was recently developed using the efficient Thomas algorithm, Hoekstra’s scheme and the splitting of the 3-D Fresnel wave equation into three 2-D wave equations [20]. However, both this splitting method and the ADI method may cause splitting errors. In this work, we demonstrate a new 3-D wide-angle BPM based on Hoekstra’s scheme that does not require the splitting of the Fresnel wave equation or the use of the ADI method. As a result, the 3-D equation is approximated more strictly by only one numerical equation. The coefficients matrix of the algebraic equation of this method is not tridiagonal but large and sparse, and its dimension (number of elements) is usually on the order of >10⁸. In this case, use of direct matrix equation solvers is prohibitively expensive and in some cases impossible even with the best available computing power. This large sparse matrix equation, however, can be solved using iterative methods, such as the bi-conjugate gradients method (BICG) and the bi-conjugate gradients stabilized method (BICGSTAB) [6].

While a different iterative method, the Orthogonal Minimum Residual (ORTHOMIN), has been proposed for 3-D structures [21], there are no published efforts demonstrating how the well-known convergence problem can be relaxed for the simulations of 3-D optical waveguide structures. When solving a large sparse matrix equation Ax = b , iterative methods have convergence problems if either the coefficient matrix A is ill-conditioned or the elements of the vector b have a wide range of values [6]. A large simulation window is needed for wide-angle simulations, which may result in an ill-conditioned large size algebraic equation with a corresponding convergence problem. In order to overcome this problem, we introduce a simulation window shifting scheme to decrease the dimension of the equation and a threshold technique to reduce the volatility of field distribution b.

2. Formulation

The 3-D scalar Helmholtz equation is given by

where k₀ 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 (SVEA), 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₀n₀z) (the rapidly varying phase term). The field is then expressed as

where n₀ is the refractive index in the cladding. Substituting Eq. (2) into Eq. (1) yields

with a = 2k ₀ n ₀ and b = k ² ₀(n ² - n ²). In the usual paraxial formulation, the second-order derivative ∣∂² ψ/∂z ²∣ in Eq. (3) is dropped, giving

Rather than dropping this second order derivative, we can discretize Eq. (3) in the z direction into the following

and use Eq. (4) at the different z positions l and l+1 to find expressions for the first derivative with respect to z

Substituting Eq. (6a) and Eq. (6b) into Eq. (5), and discretizing the equation, we find

where

Interestingly, Eq. (7) with its coefficients in Eq. (8) is same as the result when applying Padé (1, 1) approximant operator to the 3-D wave equation [16]. The second-order truncation error of the algorithm is maintained as in tridiagonal BPMs. However, the numerical algebraic equation (7) is not tridiagonal and cannot be solved using the Thomas method. Equation (7) is an M ² by M ² matrix equation for an M by M mesh grid, which requires a large amount of computer resources, further resulting in the inapplicability of direct matrix equation solvers, if the matrix is fully stored. For example, a relatively small simulation window with 200 by 200 mesh grids will have a 40,000 by 40,000 matrix equation, which requires 25.6 GB dynamic storage space only for the coefficient matrix of Eq. (7). Additional matrices of the same size are typically needed in direct matrix equation solvers. As a result, direct methods would be prohibitively expensive and in some cases impossible even with the best available computing power. Some direct solvers can store the matrix sparsely, i.e. only the non-zero elements. However, Gaussian elimination method or LU decomposition method in these direct solvers might generate non-sparse triangular matrices, which will require a large amount of computer memory as the matrix is fully stored. Whether the generated triangular matrices are sparse or not is quite problem dependent. Fortunately, each row of the coefficient matrix of Eq. (7) has no more than five non-zero values. As a result, this sparse matrix equation can be efficiently solved using iterative methods [6].

Wide-angle simulations, such as those for a tilted waveguide or a branch of a Y-switch, need a large simulation window, which can also be computationally demanding. To ease this problem, we introduce a simulation window shifting technique, in which the window is shifted in the transverse direction every few steps, maintaining the main part of the pulse around the center of the window and avoiding computations involving large numbers of edge cells. In other words, the window is numerically shifted along the flow of the light energy. This is performed by dropping cells at the two edges closer to the origin in transverse directions x and y and generating the same number of cells at the two opposite edges. These new cells are initialized with zero field values. The number of steps along the direction of propagation before the window is shifted depends on the waveguide geometry, but is such that the peak position is kept within several grid points of the origin. This technique can be directly applied to a tilted waveguide structure. It can also be applied to a Y-branch by splitting the simulation window into two after the two branches are decoupled.

Iterative methods may have convergence problems when the coefficient matrix of the equation A is ill-conditioned or the vector b is widely distributed. In our simulations, all nonzero elements of the coefficient matrix are in a narrow range, indicating that the matrix is not ill-conditioned. However, the values of the elements of the vector b are widely distributed due to the wide range of the field [Ψ] in the equation. For example, a normalized Gaussian beam with a waist radius of 3.0 μm drops exponentially to the order of 10^-16 roughly 18 μm transversely away from the center of the pulse. In order to implement the simulations using iterative methods, we narrow the field distribution by setting the elements of the field profile to zero if they are below some threshold. This reduces the volatility of the field profile to a large extent. Simulations performed without using the threshold technique confirm that iterative solvers will not converge if the threshold technique is not used.

3. Simulations and discussions

In order to show the effectiveness of our method, we carried out simulations on several situations for which we have both analytic solutions and results from our previous method [20]. The first example is a 3-D Gaussian beam with a waist radius w₀ = 3.0 μm 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. We also analyze the propagation in a single-mode channel waveguide, with the same structure and parameters as the second example in reference 20, namely, a 4 μm by 4 μm channel waveguide with the refractive indices of the core and the cladding of 1.55 and 1.52 respectively, again tilted at 30-degrees with respect to z axis and 45 degrees with respect to the x and y axes. The wavelength is 0.85 μm for both cases. Table 1 shows the comparison of the relative L² norm errors and relative position shift with respect to the analytical profile for both cases between the output profiles at z = 60 μm the simulation results obtained using the classical BPM (hereafter referred to as method 1) [1-6], the 3-D method from [20] (hereafter referred to as method 2), and our new iterative method (hereafter referred to as method 3). As can be seen, both methods 2 and 3 give more accurate results than method 1, and method 3 gives both better relative L² norm errors and better relative position shifts. Simulations were also performed using different step sizes (Δz) in the range from 0.05 μm to 0.25 μm, only resulting in slight or even indistinguishable differences in output. The dependence of relative L² norm errors and relative position shifts on the values of the threshold were tested in the range of 10^-4-10^-16 for the Gaussian beam propagation case, and no significant change was found. This is also validated in the simulations of the channel waveguide simulation, in which the refractive index contrast between the core and the cladding may change the range of the field profile. However, no matter how the range of the field profile is changed, a different threshold can be selected to ensure the convergence. Based on this technique, zeros can be directly inserted to the edge of the new simulation windows if the threshold is satisfied when shifting the simulation window.

Table 1. The relative L2 norm errors and relative position shift

View Table

^a Δy is set equal to Δx, and Δz = 0.2 μm for all simulations. ^ball errors and shifts are in percentage.

The incorporation of the wide-angle effects enables the more accurate simulation for devices, while the use of the window shift significantly reduces the computational time required. For example, a simulation of the propagation in the tilted channel waveguide described above, with a step size of 0.1 μm, required ∼380 seconds with our new method, compared to ∼1000 seconds with a traditional BPM, and ∼950 seconds with our non-iterative 3D wide angle method, described earlier [20]. These simulations were all run on an Intel P4 based PC with a 3.0 GHz CPU and 2 GB DDR400 RAM. In addition, our simulations do not seem to show any problems of convergence at longer path lengths. Changing the total path length from 60 to 300 μm, while keeping the step size constant at 0.2 μm, had no effect on the number of iterations needed at each step (between 10 and 20) to reach a fixed tolerance of 10^-6.

In summary, we derived a new form of 3-D wide-angle BPM based on Hoekstra’s scheme and an iterative matrix solution, 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. This method yields good convergence and low errors. The simulation window shifting technique is useful for wide-angle simulations of large photonic structures. The threshold technique should be adopted to ensure convergence for any relatively large simulation windows which may result in widely distributed field profile. The threshold technique, however, cannot be applied if there is a considerable amount of radiation leaving the simulation window, from which no widely distributed field profile results. We are currently trying to adapt the new algorithm to a semi-vectorial formulation and investigating other large sparse matrix equation solving methods. This method also enables us to develop higher order 3-D wide-angle algorithms using Padé approximant operators and multistep method as was developed by Hadley for 2-D simulations [16, 17].