1. Introduction

The beam propagation method (BPM) has become one of the most widely used techniques for the study of optical waveguide devices [1,2]. For the paraxial approximation or Fresnel equation used in BPM, this is usually accomplished by split-step methods, including the alternating-direction implicit (ADI) methods that are both fast and easy to implement [3]. However, the split-step schemes are only first-order-accurate in the step size, and a direct matrix inversion is required for acceptable accuracy if WA propagation is needed. 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]. Efforts have been taken to relax the limitations for WA simulations. Different treatments of WA-BPM based on the slowly varying envelope approximation have been developed, including the rational approximants of the square root operator [5], the one-way propagator [6], the exponential of the square root operator [7], and the Padé approximant operator [8], for rectangular coordinates as well as an oblique coordinate system [9]. The Padé-approximant-based WA-BPM is one of the most commonly used techniques for modeling optical waveguide structures. However, the method was originally limited to 2D structures due to the lack of efficient solvers. A 3D WA-BPM based on Hoekstra’s scheme was recently developed using the efficient Thomas algorithm and the splitting of the 3D Fresnel wave equation into three 2D wave equations [2]. However, this may cause splitting errors. Recently, C. Ma et al. [10] presented a new 3D WA-BPM also based on Hoekstra’s scheme that does not require the splitting of the Fresnel wave equation or the use of the ADI method. By using a technique for shifting the simulation window to reduce the dimension of the numerical equation and a threshold technique to further ensure its convergence, this approach shows accuracy and effectiveness. However, the resultant propagation scheme can be very slow if either the problem size is large or the structure or the boundary conditions are changing as the propagation proceeds, requiring frequent reinversions of the propagation matrix. Thus, it is imperative to find more efficient solution methods for 3D WA-BPM.

Recently, the complex Jacobi iterative method, a new iterative technique for solution of the indefinite Helmholtz equation, was introduced [11]. It is based on a complex generalization of the point relaxation technique proposed by Jacobi in 1845, and has been shown to converge at a rate dependent only upon the grid size and effective absorption coefficient. For beam propagation of wave profiles with a 2D cross section, the wide angle beam propagation equation can be recast in terms of a Helmholtz equation with source term, and the effective absorption coefficient appearing in this equation is high, leading to rapid convergence of the method. In this work we apply the complex Jacobi iterative (CJI) method for beam propagation and show it to be highly efficient for the solution of large problems, compared with existing techniques.

The beam propagation equation for analysis of optical propagation in waveguide structures is based on a modified Padé(1,1) approximant operator we propose here, which gives evanescent waves the desired damping. Our modified Padé approximant propagation operators allow more accurate approximation to the true Helmholtz equation and faster convergence when the CJI method is used for the modified Padé approximant-based WA-BPM than those of the standard Padé approximant approach [8]. Furthermore, since the utility of the CJI technique depends mostly upon its execution speed in comparison with the direct matrix inversion (DMI) method, we also present several speed comparisons. Numerical implementations are carried out for both 2D and 3D optical waveguide structures. The rest of this paper is structured as follows: In Section 2, we will describe the novel modified Padé approximant operators. Next, WA-BPM formulations using the CJI and DMI methods are presented in Section 3. In Section 4, the investigation of convergence rate of CJI is discussed. In Section 5, the benchmark results performed by WA-BPM using CJI and DMI are presented followed by conclusions in Section 6.

2. Modified Padé approximant operators

The scalar Helmholtz equation obtained by using the slowly varying envelope approximation is given by [8],

where $P={\nabla }_{\perp }^2+k_0^2\left(n^2-n_{\mathrm{ref}}^2\right)=\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+k_0^2\left(n^2-n_{\mathrm{ref}}^2\right)$ with k=k ₀ n_ref, n the refractive index refractive index profile, n_ref the reference refractive index, k ₀ the vacuum wavevector.

We may formally rewrite Eq. (1) in the form

Equation (2) suggests the recurrence relation

Hadley [8] proposed the rational approximation of WA beam propagation using Padé approximant operators with initial value of ${\frac{\partial }{\partial z}\mid }_0=0.\mathrm{For}{\frac{\partial }{\partial z}\mid }_2$, this gives us the well-known Padé (1,1) approximant-based WA beam propagation formula as follows:

If Eq. (4) is compared with a formal solution of Eq. (1) written in the well-known form.

where $X=\frac{P}{k^2}$, we obtain the approximation formula

Since the operator X has a real spectrum, it is useful to consider the approximation of $\sqrt{1+X}-1$ by the Padé approximant propagation operator. Figure 1 shows the absolute value of $\sqrt{1+X}-1$ and its first-order Padé approximant called Hadley (1,1) as a function of X.

However, as the denominator of Hadley (1,1) gradually approaches zero, its absolute value approaches ∞ as can clearly be seen in Fig. 2. Physically this means that the standard Padé approximant incorrectly propagates the evanescent modes.

To circumvent this problem we propose a modified Padé approximant operator. First of all, following [12], by multiplying both sides of Eq. (2) with $-\frac{i}{k}$, we obtain

We may rewrite Eq. (7) as follows

where $f\left(X\right)=-\frac{i}{k}\frac{\partial }{\partial z}$.

Eq. (8) suggests the recurrence relation

Y. Y. Lu [12] has proved that Eq. (9) can provide a good approximation to $\sqrt{1+X}-1$ with the initial value of

Subsequently, we use this fact to go back to the original recurrence relation (3) and construct modified Padé approximant operators by using a different initial value of ${\frac{\partial }{\partial z}\mid }_0=-k\beta$. For ${\frac{\partial }{\partial z}\mid }_2$ the first-order modified Padé(1,1) approximant operator is given as follows:

The absolute value of the modified Padé(1,1) approximant of $\sqrt{1+X}-1$ is also depicted in Fig. 1. It is seen that our modified Padé approximant operator (with β=2) allows more accurate approximations to the true Helmholtz equation than the standard Padé approximant operator. Furthermore, the standard rational Padé approximant incorrectly propagates the evanescent modes as their denominator gradually approaches zero while the modified Padé approximant gives the waves propagating in the evanescent region the desired damping as clearly seen in Fig. 2.

 

Fig. 1. The absolute values of (1+X)^1/2-1 (solid line), its first-order standard Padé approximant (X/2)/(1+X/4) (solid line with circles) and modified Padé approximant (X/2)/(1+X/{4(1+ibeta/2)}) (dotted line).

Download Full Size | PDF

 

Fig. 2. The absolute value of (1+X)^1/2-1 (solid line), the first-order standard (X/2)/(1+X/4) (solid line with circles) and modified (X/2)/(1+X/{4(1+ibeta/2)}) (dotted line) Padé approximant of (1+X)^1/2-1 with respect to X.

Download Full Size | PDF

3. WA-beam propagation formulation

3.1 Basic equation

By using the modified Padé(1,1) approximant, the 3D semivectorial WA beam propagation equation can be written as follows [13]:

where $\xi =\frac{1}{4k^2(1+i\beta /2)}-\frac{i\Delta z}{4k}$, ξ*the complex conjugate of ξ and Δz the propagation step.

3.2. WA-BPM using CJI

By dividing both sides of Eq. (12) by ξ, it may be written as an inhomogeneous Helmholtz equation

or

Thus the beam propagation can be cast as a 2D Helmholtz equation with source term in an effective medium with loss of $4k_0^2{n}_{\mathrm{ref}}^2\frac{k\Delta z}{{(1+k_0^2{n}_{\mathrm{ref}}^2\beta \Delta z/2)}^2+k^2\Delta z^2}$. This loss is high for a typical choice of k Δz. This is a condition that favors rapid convergence for the CJI method.

3.3 WA-BPM using DMI

By discretizing Eq. (13), we find that

where

Actually, the coefficients in Eq. (16) are the same as the result when applying Hoekstra’s scheme to the 3D wave equation [10]. Eq. (15) is an M ² by M ²matrix equation for an M by M mesh grid. However, each row of the coefficient matrix has no more than five non-zero values. As a result, this sparse matrix equation can be efficiently solved using various methods [3]. In our calculations, the sparse matrix solver-UMFPACK package has been used [14].

4. Convergence studies of CJI

In this section, we investigate the convergence rate of the CJI method for 3D WA-BPM based on the modified Padé(1,1) approximant and Hadley(1,1). We simulate the propagation of a 3D Gaussian beam through a symmetric Y-branch waveguide. The structure parameters of the Y-branch waveguide are the same as in [17] with d ₁=2µm and d ₂=1µm. In Fig. 3 we show that the CJI method converges faster with WA-BPM based on the modified Padé(1,1) approximant than those of Hadley(1,1). However, it suffers from the fact that the iteration count between two successive 2D cross sections increases throughout the propagation direction. To overcome this problem we propose the use of a perfectly matched layer (PML), which can absorb incident radiation without any additional parasitic reflections, regardless of wavelength, incident angle or polarization as boundary conditions [15–16]. It is seen that the CJI technique becomes more stable with the use of a PML as shown in Fig. 3.

 

Fig. 3. The iteration count per propagation step for propagation through a symmetric Y-branch waveguide obtained using the CJI method for standard Padé approximant-based WA BPM without (1) and with (2) PML and the CJI method for the modified Padé without (3) and with (4) PML. (um is defined as µm)

Download Full Size | PDF

5. Benchmark results

We now employ the WA-BPM using the new CJI and the traditional DMI methods to perform benchmark tests on both 2D and 3D optical waveguide structures. The simulations were all run on a notebook PC using Matlab. For the 2D case, we consider a 10-degree tilted waveguide [7] and Y-branch waveguide [17]. In the tilted waveguide the fundamental mode for the slab of width w=0.2µm is propagated through 30µm at wavelength λ=1.55µm in a medium of refractive index n=3.4 with the propagation step size of Δz=0.1µm. With a very strict propagation error tolerance of 10^-9 the CJI method only took 84.8 seconds, whereas the DMI method took 296.5 seconds.

In the Y-branch waveguide the parameters needed for calculation are the same as in [17]. With a small propagation step size Δz=0.01µm (requiring frequent matrix inversion) the DMI method performed the propagation in 1414.8 seconds while the CJI method took only 584.1 seconds. It is obvious that for these 2D waveguide structures the CJI method is faster than DMI.

For the 3D case, we consider Gaussian beam propagation in a straight rib waveguide [18- 19] and guided-mode propagation in a Y-branch rib waveguide. The width and height of the straight rib waveguide are w=2µm and h=1.1µm, as seen in Fig. 4 of [19]. The guiding core has an index n_f=3.44 and a thickness t=0.2µm while the refractive index of substrate and cover is n_s=3.4 and n_c=1, respectively. The Gaussian beam with a waist radius w ₀=0.3µm has been injected into the rib waveguide at wavelength λ=1.55µm. Due to the large memory required for DMI, the small computational window of 2x2µm is discretized with a grid size of Δx= Δy=0.1µm, and the short path length of 2µm is discretized with a propagation step size Δz=0.1µm. The resulting runtime of DMI is 177.9 seconds while runtime for CJI is only 4.7 seconds.

For a Y-branch, the initial rib waveguide is split into two 5-degree tilted waveguides as shown in Fig. 4. The longitudinal dimension is h₁=1µm. The other structure parameters are the same as the above straight rib waveguide. The fundamental TE mode of the ridge waveguide of width w=2µm at 1.55µm wavelength is used as the excited field at z=0. The propagation step size is Δz=0.1µm. The field pattern at z=3µm calculated by DMI and CJI are depicted in Fig. 5(a) and 5(b), respectively. Due to the high effective loss in the propagation medium the complex Jacobi method performed the propagation in only 5.9 seconds while DMI required 268.9 seconds.

Table 1. Quantitative comparison of runtimes of the direct matrix inversion and the complex Jacobi iteration for WA beam propagation in waveguide (WG) structures

View Table

Table 1 shows the performance of the two methods for the optical waveguide structures chosen here. It is clearly seen that the runtime of the iterative method is substantially lower than that of the DMI method. For large problems requiring very large storage space and also for structures with a long path length with small propagation step size that require frequent matrix inversions, the DMI technique is numerically very intensive. In contrast, for typical choices of k Δz the CJI technique offers rapid convergence and shorter runtimes.

 

Fig. 4. Y-branch optical rib waveguide

Download Full Size | PDF

 

Fig. 5. Intensity contours for TE mode propagating in a 3D Y-branch rib waveguide at z=3µm calculated by (a) DMI and (b) the CJI method.

Download Full Size | PDF

4. Conclusions

A new complex Jacobi iterative method adapted for the solution of 3D WA beam propagation has been presented. We proposed the modified Padé approximant of the WA beam propagation operator that gives the waves propagating in the evanescent region the desired damping. The resulting approach allows accurate approximations to the true Helmholtz equation. Furthermore, a quantitative comparison of runtimes between the traditional direct matrix inversion and the new complex Jacobi iterative method for both 2D and 3D WA beam propagation demonstrates convincingly that the complex Jacobi iterative method is very competitive for demanding problems. This solution technique will also enable a development of higher order 3D Padé approximant-based WA-beam propagation algorithms using the multistep method [20], which will be presented in a future publication.