## Abstract

Applying continuous and discrete transformation techniques, new analytical expressions to calculate dispersion and stability of a Runge-Kutta Finite Difference Beam Propagation Method (RK-FDBPM) are obtained. These expressions give immediate insight about the discretization errors introduced by the numerical method in the plane-wave spectrum domain. From these expressions a novel strategy to adequately set the mesh steps sizes of the RK-FDBPM is presented. Assessment of the method is performed by studying the propagation in several linear and nonlinear photonic devices for different spatial discretizations.

© 2008 Optical Society of America

## 1. Introduction

Scalar [1,2], semivectorial [3,4] and vectorial [5] Finite Difference BPM (FDBPM) techniques are currently being used in the analysis and design of optical integrated circuits. There exists a wide variety of examples in the literature reporting the use of these techniques to analyze different types of photonic devices, including: passive waveguides (splitters and couplers in SOI technology [1], MMI couplers in low contrast technologies [2,6,7], polarization splitters [3] or reconfigurable insertion/extraction multiplexors/demultimplexors (ROADM) in AWG architectures [5]); active devices such as VCSEL lasers [6,8]; and other complex structures comprising linear and nonlinear devices such as bistable switches based on BLD lasers combined with non-linear couplers or with multimode interference devices [6,7,9]. Although FDBPM techniques are quite mature, research is still being carried out mainly concerning wide angle capabilities. The most common contributions in this field involve making improvements in the paraxial approximation, or using new techniques for the discretization of the longitudinal operator to increase the computational efficiency [10–13]. Nevertheless, in most cases there is a lack of analytical tools that can be used ‘a priori’ to set the numerical discretization mesh sizes that ensure a certain maximum amount of simulation error for a given device under analysis. Designers are thus forced to either acquire in-depth knowledge of the numerical technique’s details or to follow complicated and time-consuming trial and error iterative processes to make sure that the results of the simulations are sufficiently accurate.

In this article an analytical technique for studying the numerical dispersion and losses of the FDBPM is presented. A design criterion for the parameters of this method is introduced, so that the proposed strategy makes it possible to set up the proper mesh step sizes for any particular device under analysis. This strategy has been successfully applied to several FDBPM techniques in the analysis of linear and non-linear devices [14]. In particular, it has been applied to paraxial FDBPM techniques based on the Crank-Nicolson scheme (CN-FDBPM), the ADI technique (ADI-FDBPM) and the Runge-Kutta method (RK-FDBPM), but it can be used with any wide angle FDBPM as the one proposed in [12]. Nevertheless, to be brief, only its application to the semivectorial RK-FDBPM previously published by the authors will be discussed herein [13].

## 2. RK-FDBPM method

Under the paraxial approach, the propagation of the transverse electric field is characterized by the vectorial Fresnel equation [15],

where *ē*=*e _{x}xô*+

*e*is the slowly varying complex envelope of the electric field

_{y}ŷ*ε*, which can be written as

where *k _{0}* is the vacuum propagation constant and

*n*is the BPM reference refractive index.

_{ref}In many practical devices, coupling between different transversal components of the field can be neglected and the semivectorial BPM approximation [15] can be applied, yielding

where *M _{uu}* is the semivectorial transversal differential operator

where *u* and *v* can indistinctly be the *x* and *y* variables and *n* is the device refractive index, which can be linear or non-linear [13].

Applying the RK-FDBPM approach [13], the partial differential equations from (3) are first discretized transversally by means of a central second order finite difference scheme, providing a system of differential equations [15]

where $\overline{\Psi}$
(*z*) is a column vector containing either the e_{x} or e_{y} field components discretized in the transversal mesh points and ${\overline{\overline{M}}}_{T}$ is the system pentadiagonal matrix (becoming tridiagonal for two-dimensional problems). Thus, if *e ^{m,n}_{u}*(

*z*)=

*e*(

_{u}*m*Δ

*x*,

*n*Δ

*y*,

*z*) is the value of the field in the mesh coordinates (

*m*,

*n*), then the

*k*-th element of $\overline{\Psi}$ (z) is obtained as Ψ

_{k}(

*z*)=

*e*(

^{m,n}_{u}*z*), where

*k*=(

*m*-1)

*N*+

_{x}*n*and

*N*is the number of mesh points in the computational window on the

_{x}*x*axis. In other words, vector$\overline{\Psi}$ (z) contains the sampled electric field at position z, in a column-wise arrangement.

Finally, the system of ordinary differential equations in Eq.(5) is discretized in the longitudinal direction by means of the fourth order Runge-Kutta technique [13,14,16]. For linear and z-invariant devices, the Runge-Kutta algorithm results in

where vector $\overline{\Psi}$
_{p} contains the discretized field at each numerical propagation step *z _{p}*=

*p*Δ

*z*i.e. $\overline{\Psi}$

_{p}=$\overline{\Psi}$ (

*p*Δ

*z*) and ${\overline{\overline{M}}}_{T}$ is z independent. Therefore, the successive field derivatives are obtained from Eq.(5) as

so that Eq.(6) describes the field at *z*
_{p+1}=(*p*+1)Δ*z* as a fourth order Taylor series expansion around *z _{p}*=

*p*Δ

*z*.

## 3. RK-FDBPM dispersion analysis

For linear and *z*-invariant problems, the semivectorial Fresnel equation, Eq.(3), could be considered to be a continuous-type linear and invariant system. The processed signal is the transverse electric field envelope, which depends on three independent continuous variables (*x*,*y*,*z*). Also, the difference equations obtained by means of the different transverse and longitudinal discretization strategies, Eq.(5) and Eq.(6), may be considered to be a discrete-type linear and invariant system in which the processed signal is the discretized electric field at each transversal mesh point, which therefore depends on three discrete independent variables (*m*,*n*,*p*), which are related to the continuous independent variable by the corresponding sampling periods, i.e. *x _{m}*=Δ

*xm*,

*y*=Δ

_{n}*yn*,

*z*=Δ

_{p}*zp*.

The dispersion analysis involves studying the behavior of both continuous and discrete systems by comparing their transfer functions, obtained from the correspondent dispersion relations [17]. What is intended is to compare the way how each component of the plane-wave spectrum, in which the input field profile can be expanded, propagates along the structure for both the continuous system and its discretized counterpart.

For linear and homogeneous media, transfer functions and dispersion relations are easily obtained for both continuous and discrete systems. These expressions will be later successfully applied to the study of inhomogeneous problems.

#### 3.1 Dispersion relations for the Fresnel equation in linear and homogenous media

The Fresnel equation for three-dimensional homogeneous media is given by

where *n* is now a constant. The spatial spectral analysis of any wave front can be carried out using the Fourier transform [17]. In the three-dimensional case the Fourier transform is

and the Fourier transform of Eq.(8) is [18]

This equation provides a solution for the field which differs from the trivial one only at the points of the complex Fourier space which are described by the dispersion relation *D _{u}*=0, so

Likewise, for the two-dimensional case, which is invariant in *x*, the dispersion equation is

#### 3.2 Dispersion relations for the RK-FDBPM in linear and homogenous media

Considering the field *e ^{m,n,p}_{u}*, which is sampled in the three spatial dimensions, as a discrete sequence, its three-dimensional Z transform can be defined as

Taking the Z transform of Eq.(6), a similar expression to Eq.(10) is obtained in the transformed domain for the RK-FDBPM discretization algorithm. To find out the corresponding dispersion relation, we substitute

where Ω_{x}, Ω_{y} y Ω_{z} are the corresponding angular frequencies expressed in rad/sample and A_{z} is the attenuation function in Np/sample. Contrary to its continuous counterpart, in this case the Z variable has been allowed to become complex to take into account possible numerical losses of the discretization algorithm. Once again, the relation between the frequencies of the continuous and discrete domains is established by means of the sampling periods

where superscript *RK* differentiates the numerical method propagation constant from that of the continuous equation. Using Eq. (13) to Z transform the difference equation of the Runge-Kutta algorithm in Eq.(6) and substituting Eq. (14) and Eq. (15) into the transformed equation, it is possible to calculate the dispersion relations of the RK-FDBPM as [14]

$$+\left(2\Delta z{\stackrel{\u033f}{M}}_{T,1}+\Delta {z}^{2}{\stackrel{\u033f}{M}}_{T,1}^{2}+\frac{\Delta {z}^{3}}{3}{\stackrel{\u033f}{M}}_{T,1}^{3}+\frac{\Delta {z}^{4}}{12}{\stackrel{\u033f}{M}}_{T,1}^{4}\right)\mathrm{cos}\left({\beta}_{x}\Delta x\right)$$

$$+\left(2\Delta z{\stackrel{\u033f}{M}}_{T,{N}_{x}}+\Delta {z}^{2}{\stackrel{\u033f}{M}}_{T,{N}_{x}}^{2}+\frac{\Delta {z}^{3}}{3}{\stackrel{\u033f}{M}}_{T,{N}_{x}}^{3}+\frac{\Delta {z}^{4}}{12}{\stackrel{\u033f}{M}}_{T,{N}_{x}}^{4}\right)\mathrm{cos}\left({\beta}_{y}\Delta y\right)$$

$$+\left(\Delta {z}^{2}{\stackrel{\u033f}{M}}_{T,2}^{2}+\frac{\Delta {z}^{3}}{3}{\stackrel{\u033f}{M}}_{T,2}^{3}+\frac{\Delta {z}^{4}}{12}{\stackrel{\u033f}{M}}_{T,2}^{4}\right)\mathrm{cos}\left(2{\beta}_{x}\Delta x\right)$$

$$+\left(\Delta {z}^{2}{\stackrel{\u033f}{M}}_{T,2{N}_{x}}^{2}+\frac{\Delta {z}^{3}}{3}{\stackrel{\u033f}{M}}_{T,2{N}_{x}}^{3}+\frac{\Delta {z}^{4}}{12}{\stackrel{\u033f}{M}}_{T,2{N}_{x}}^{4}\right)\mathrm{cos}\left(2{\beta}_{y}\mathrm{\Delta y}\right)$$

$$\phantom{\rule{.2em}{0ex}}+\left(\frac{\Delta {z}^{3}}{3}{\stackrel{\u033f}{M}}_{T,3}^{3}+\frac{\Delta {z}^{4}}{12}{\stackrel{\u033f}{M}}_{T,3}^{4}\right)\mathrm{cos}\left(3{\beta}_{x}\Delta x\right)+\left(\frac{\Delta {z}^{3}}{3}{\stackrel{\u033f}{M}}_{T,3{N}_{x}}^{3}+\frac{\Delta {z}^{4}}{12}{\stackrel{\u033f}{M}}_{T,3{N}_{x}}^{4}\right)\mathrm{cos}\left(3{\beta}_{y}\Delta y\right)$$

$$+\frac{\Delta {z}^{4}}{12}{\stackrel{\u033f}{M}}_{T,4}^{4}\mathrm{cos}\left(4{\beta}_{x}\Delta x\right)+\frac{\Delta {z}^{4}}{12}{\stackrel{\u033f}{M}}_{T,4{N}_{x}}^{4}\mathrm{cos}\left(4{\beta}_{y}\Delta y\right)$$

$$+\left(2\Delta {z}^{2}{\stackrel{\u033f}{M}}_{T,{N}_{x}+1}^{2}+\frac{2\Delta {z}^{3}}{3}{\stackrel{\u033f}{M}}_{T,{N}_{x}+1}^{3}+\frac{\Delta {z}^{4}}{6}{\stackrel{\u033f}{M}}_{T,{N}_{x}+1}^{4}\right)\mathrm{cos}\left({\beta}_{x}\Delta x\right)\mathrm{cos}\left({\beta}_{y}\Delta y\right)$$

$$+\left(\frac{2\Delta {z}^{3}}{3}{\stackrel{\u033f}{M}}_{T,{N}_{x}+2}^{3}+\frac{\Delta {z}^{4}}{6}{\stackrel{\u033f}{M}}_{T,{N}_{x}+2}^{4}\right)\mathrm{cos}\left(2{\beta}_{x}\Delta x\right)\mathrm{cos}\left({\beta}_{y}\Delta y\right)$$

$$\phantom{\rule{.2em}{0ex}}+\left(\frac{2\Delta {z}^{3}}{3}{\stackrel{\u033f}{M}}_{T,{2N}_{x}+1}^{3}+\frac{\Delta {z}^{4}}{6}{\stackrel{\u033f}{M}}_{T,2{N}_{x}+1}^{4}\right)\mathrm{cos}\left({\beta}_{x}\Delta x\right)\mathrm{cos}\left(2{\beta}_{y}\Delta y\right)+\text{}$$

$$+\frac{\Delta {z}^{4}}{6}{\stackrel{\u033f}{M}}_{T,{N}_{x}+3}^{4}\mathrm{cos}\left(3{\beta}_{x}\Delta x\right)\mathrm{cos}\left({\beta}_{y}\Delta y\right)+\frac{\Delta {z}^{4}}{6}{\stackrel{\u033f}{M}}_{T,3{N}_{x}+1}^{4}\mathrm{cos}\left({\beta}_{x}\Delta x\right)\mathrm{cos}\left(3{\beta}_{y}\Delta y\right)$$

$$+\frac{\Delta {z}^{4}}{6}{\stackrel{\u033f}{M}}_{T,2{N}_{x}+2}^{4}\mathrm{cos}\left(4{\beta}_{x}\Delta x\right)\mathrm{cos}\left(4{\beta}_{y}\Delta y\right)],$$

where *N _{x}* and

*N*are the number of samples of the computational window in the

_{y}*x*and

*y*direction, respectively. The size of ${\overline{\overline{M}}}_{T}$ is

*N*×

_{x}N_{y}*N*, for linear homogeneous media it is symmetric and constant, meaning that all the elements in each subdiagonal are identical, so ${\overline{\overline{M}}}_{T}$(

_{x}N_{y}*k*,

*k*+

*p*)=${\overline{\overline{M}}}_{T}$(

*k*,

*k*-

*p*), ∀

*k*∈[1,

*N*]. Each coefficient ${\overline{\overline{M}}}_{T,i}^{n}$ in Eq. (16) is any of the i-th subdiagonal elements in(

_{x}N_{y}*M̿*)

_{T}^{n}, where

*M̿*

^{n}_{T,0}represents the diagonal elements.

In the two-dimensional *x*-invariant case, the obtained dispersion equation is,

$$+\left(2\Delta z{\stackrel{\u033f}{M}}_{T,1}+\Delta {z}^{2}{\stackrel{\u033f}{M}}_{T,1}^{2}+\frac{\Delta {z}^{3}}{3}{\stackrel{\u033f}{M}}_{T,1}^{3}+\frac{\Delta {z}^{4}}{12}{\stackrel{\u033f}{M}}_{T,1}^{4}\right)\mathrm{cos}\left({\beta}_{y}\mathrm{\Delta y}\right)$$

$$+\left(\Delta {z}^{2}{\stackrel{\u033f}{M}}_{T,2}^{2}+\frac{\Delta {z}^{3}}{3}{\stackrel{\u033f}{M}}_{T,2}^{3}+\frac{\Delta {z}^{4}}{12}{\stackrel{\u033f}{M}}_{T,2}^{4}\right)\mathrm{cos}\left(2{\beta}_{y}\mathrm{\Delta y}\right)$$

$$+\left(\frac{\Delta {z}^{3}}{3}{\stackrel{\u033f}{M}}_{T,3}^{3}+\frac{\Delta {z}^{4}}{12}{\stackrel{\u033f}{M}}_{T,3}^{4}\right)\mathrm{cos}\left(3{\beta}_{y}\Delta y\right)+\left(\frac{\Delta {z}^{4}}{12}{\stackrel{\u033f}{M}}_{T,4}^{4}\right)\mathrm{cos}\left(4{\beta}_{y}\Delta y\right)].$$

## 4. Selecting RK-FDBPM discretization parameters

For homogeneous media, Eq. (12) and Eq. (17) make it possible to calculate the discretization phase error for a two-dimensional problem. In the three-dimensional case, Eq. (11) and Eq.(16) must be used. In both cases, the cumulative phase error along a propagating distance *L* can be calculated as

The proposed technique for selecting the RK-FDBPM discretization parameters involves estimating the bandwidth of the plane-wave spectrum of the propagating field profile. Then the mesh step sizes will be chosen which are small enough to yield an accumulated phase error which is much smaller than 2π at the superior end of the plane-wave spectrum band.

Homogeneous media are assumed in the calculation of the phase error curves, taking the value of the analyzed structure fundamental mode effective index for the refractive index n. Due to intrinsic paraxial limitations of the BPM technique, in practice the homogenous medium approach is usually sufficient to limit the dispersion. Once the optic envelope bandwidth has been estimated and the dispersion curves for the homogeneous problem have been calculated, the following condition is enforced in the selected bandwidth

Finally, calculating the real part of Eq. (17) and Eq. (16) yields the numerical losses for the selected mesh and provides information about the method’s stability.

In this paper, the proposed discretization criterion for linear problems is based on calculating the envelope’s spatial bandwidth, where 90% of the plane wave spectrum energy is located; then the mesh steps are set up to achieve an accumulated phase error which is lower than 0.1×2π rad. For nonlinear problems it is more complicated to propose a universal criterion, since operating principles may be very different among different devices. For the soliton propagation analyzed in this paper, the most important issue is to preserve the equilibrium between nonlinear and dispersion effects, so the mesh steps will be selected to completely eliminate the cumulative phase error. The goal is considered to have been achieved if the total phase error is smaller than 0.01×2π rad/µm.

As previously commented, the proposed technique may be applied to any FDBPM. In case of wide angle FDBPM, the procedure described in section *3.1* can be readily applied to the wide angle continuous partial derivatives equation to calculate the continuous dispersion equation. Also, provided wide angle FDBPM algorithms have the form [12]

their dispersion equation may be easily obtained by means of the procedure described in section *3.2*. In fact, Eq. (20) has the same form as the ones obtained from paraxial Crank-Nicolson or ADI FDBPM, and the proposed dispersion characterization technique has already been applied in both cases [14].

## 5. Results

To show the validity of the proposed method, several devices (including linear and nonlinear structures, two dimensional and three dimensional cases) were analyzed by means of the semivectorial RK-FDBPM with complex coordinate PML [14]. Dispersion analysis was performed by means of the expressions in Eq.(16) and Eq.(17), and discretization parameters were varied to show the effects of the accumulated phase errors on simulated propagation.

#### 5.1. 2D linear devices

As a first example, Eq. (12) and Eq. (17) were used to select the proper mesh steps to ensure the correct calculation of fundamental mode propagation in a tilted slab.

Slab refractive indexes are n_{f}=1.571 in the film, n_{s}=1.57 in the substrate and n_{c}=1.569 in the cover. The film’s width is 10 µm, the slab propagation axis is tilted θ=10° from *z* axis and the excitation wavelength is λ_{0}=633 nm.

In Fig. 1 phase error curves for the different tested meshes are shown. In all cases the propagation length is 250 µm. The medium refractive index and the BPM reference index were both set equal to the fundamental mode effective index (n_{eff}) at λ_{0}. Mesh steps are Δy=0.5 µm, Δz=0.1 λ_{0} in the case of (a)-mesh, Δy=0.1 µm, Δz=0.1 λ_{0} for (b)-mesh, and Δy=0.025 µm, Δz=0.0125 λ_{0} in the case of (c)-mesh.

The field profile at *z*=*0* is the tilted slab fundamental mode, so the excitation spectrum is centered at β_{y}=-k_{0}
*n _{eff}* sinθ=-2.7 rad/µm, and 90% of the spatial spectrum energy is located in the range from -3 rad/µm to -2.47 rad/µm, so it will be considered a spatial bandwidth of -3 rad/µm. For this βy value, one can see in Fig. 1 how the phase error is not much smaller than 2π when the (a)-mesh is used. In the case of the (b)-mesh, the phase error is approximately 0.085×2π, and for the (c)-mesh it is roughly 0.01×2π. In the latter two cases, the numerical dispersion is considered to be sufficiently small.

In Fig. 2 the propagation of TE_{0} excitation is shown. Fig. 2(a) shows the solution obtained using the (a)-mesh. It can be seen how the high numerical dispersion makes the BPM method unable to properly characterize the mode guidance along the structure. Fig. 2(b) shows both the (b)-mesh and the (c)-mesh solutions, which are indistinguishable. It can be seen how in these cases the numerical dispersion remains limited and the achieved solution is sufficiently accurate.

In order to accurately quantify the effect of the numerical dispersion, the scattering parameter S21 has been taken as a figure of merit. In the S_{21} calculation, the cross-section overlap integral of the propagating field and the fundamental mode has been solved at each propagation step, as defined in [19].

Figure 3 shows the *z* evolution of S_{21}, where the curve’s ripple is caused by the lateral misalignment between the exact theoretical position of the slab longitudinal axis and the position of the transverse mesh points at each propagation step. In Fig. 3(a) it can be seen how the S_{21} modulus remains near unity for the (b)-mesh and (c)-mesh. However, for the (a)-mesh, not far from the excitation the field profile begins to differ greatly from the fundamental mode, so the S_{21} modulus dramatically decreases. In Fig. 3(b) one can see how the S_{21} phase error agrees well with the values calculated in the dispersion analysis.

#### 5.2. 3D linear devices

Next, the application of the dispersion analysis to a three-dimensional problem is shown. Expressions in Eq.(11) and Eq.(16) were used for mesh sizing in the analysis of the MMI coupler shown in Fig. 4.

The coupler is a device with refractive index n_{1}=1.4549 embedded in a medium with index n_{2}=1.4440. At the design wavelength (λ_{0}=1550nm) the multi mode section supports 10 guided modes. The device is designed to work as a 120° 2×3 coupler with S parameters as in [20],

To perform the device analysis the phase error inside the bandwidth of the highest order MMI mode spectrum has to be limited. The field profile of the highest propagating mode in the MMI section is obtained by modal analysis. Application of Fast Fourier Transform (FFT) to this field profile yields the plane-wave spectrum which exhibits a spatial bandwidth with 90% of its energy concentrated in the band *β _{x}*,

*β*∈[-1,1](rad/µm). Therefore, this is the band in which the RK-FDBPM phase error is intended to be limited, considering a total propagation length L=3570 µm.

_{y}Figure 5(a) shows the accumulated phase error for the higher order propagating mode, after a propagation length L=3570 µm and for discretization steps Δx=Δy=1µm and Δz=0.1 µm. In this case an excessive phase error can be seen, so one would expect the numerical technique to fail to achieve a proper device characterization. Reducing the mesh sizes to Δx=Δy=0.1µm and Δz=0.01 µm, the phase error shown in Fig. 5(b) is obtained. In this case, it can be verified that the total phase error is much smaller than 2π, so it is assumed that this mesh setup leads to a proper device characterization.

In Fig. 6 the results of the MMI RK-FDBPM simulation are presented. The propagation shown in Fig. 6(a) was calculated with the mesh size considered in Fig. 5(a). It can be seen how the field solution is highly contaminated by the phase error shown in Fig. 5(a). The phase error introduced in each propagating mode prevents the image formation at the proper places in the device’s cross-section at the end of the multimode zone. On the other hand, the field in Fig. 6(b) was calculated with the mesh used in Fig. 5(b). In this case the result of the dispersion analysis shown in Fig. 5(b) predicts proper numerical technique behavior. In fact, the results shown in Fig. 6(b) agree well with the ones obtained in [20] by means of a spectral BPM technique.

#### 5.3. Non-linear devices

To validate the use of the dispersion curves in the analysis of non-linear devices, Eq. (12) and Eq. (17) were used to set the mesh sizes to analyze simultaneous propagation of TE_{0} and TM_{0} modes in a non-linear straight slab. The slab’s refractive index profile is the same as that of the linear slab analyzed in section 5.*1*, but in this case the substrate index has an additional non-linear component [13,21]. The substrate non-linear refractive index is *n*
^{2}
_{NL}=*n*
^{2}
_{s}+*α _{s}*(|

*e*|

_{x}^{2}+|

*e*|

_{y}^{2}+|

*e*|

_{z}^{2}), and the non-linear coefficien0t is

*α*=10

_{s}^{-10}

*V*

^{2}/

*m*

^{2}. The guiding structure is excited simultaneously with the linear slab TE

_{0}and TM

_{0}modes at 633nm wavelength. The excitation power is 2.92 mW/mm, equally divided among both modes. The propagation length is

*L*=2400

*λ*

_{0}.

The analyzed device is a weakly guiding structure, so the field profiles of both propagating modes are very similar and indistinguishable. In both cases, with the previously defined criterion, their transverse spectra extend into an approximate 0.3 rad/µm bandwidth centered at the origin. The phase errors for the three different meshes used are shown in Fig. 7(a). The first tested mesh is Δ*y*=0.5*µm* and Δ*z*=8*λ*
_{0}. At the upper band end (0.3 rad/µm), the normalized phase error is 1.5×10^{-3}, so from previous considerations, negligible numerical phase dispersion is expected. However, in the device simulation results shown in Fig. 8 for this specific mesh size, a couple of solitons generated at z=400λ_{0} and z=1200λ_{0} are observed. These propagate in the non-linear substrate with a propagation angle smaller than 1,5° from the longitudinal axis. It is obvious that this behavior will make the plane-wave spectrum to be shifted in frequency, making it now be centered at 0.35 rad/µm, and thus its bandwidth is shifted up to 0.65 rad/µm.

As can be seen in Fig. 7(a), for this new value of plane-wave spectrum bandwidth the accumulated phase error is no longer negligible (normalized phase error higher than 0.02) and it will cause a phase distortion which will modify the field behavior. This effect will be more noticeable in focusing Kerr non-linear devices, like the one analyzed in this example. Since soliton generation in Kerr focusing media obeys a delicate balance between dispersion and non-linear effect, this balance may be easily altered by the numerical dispersion.

To minimize the numerical dispersion effects, the mesh size is redefined. Dispersion curves show that, to diminish phase error, the simultaneous reduction of the transversal and longitudinal mesh sizes is required. In principle, the mesh steps are divided by two, setting up a mesh with Δ*y*=0.25µm and Δ*z*=4*λ*
_{0}. However, the corresponding numerical attenuation curve (calculated as the real part of Eq. (17) and shown in blue line in Fig. 7(b)) is positive in the higher frequencies, meaning that the numerical technique is unstable for this mesh setup. In order to encourage the finite difference technique’s stability, the propagation step is reduced to Δ*z*=2.5*λ*
_{0}, resulting in the Fig. 7(b) red curve, which still yields instability problems. Finally, with Δ*z*=2*λ*
_{0} the green curve in Fig. 7(b) is obtained which, being negative for all frequencies, ensures the RK-FDBPM’s stability.

The corresponding phase error curve is also plotted in red in Fig. 7(a). In this case, for a transverse frequency of 0.6 rad/µm, the total phase error is roughly 5×10^{-3}, so numerical dispersion will not affect simulation results. To prove this assumption, a third finer mesh was used. For this case, mesh sizes were set to Δ*y*=0.125*µm* and Δ*z*=*λ*
_{0}/2, which also ensure stability. The green curve in Fig. 7(a) is the phase error for this third mesh setup, showing a normalized value of 1×10^{-3} at 0.6 rad/µm.

Field profiles at z=2400λ_{0} for each of the proposed discretization meshes are plotted in Fig. 9. It can be seen that the simulation results obtained with the two last mesh setups are formally identical.

Additional insight can be gained into understanding the effects of phase errors by analyzing the results of Fig. 9. Three different pulses can be clearly observed in Fig. 9 (corresponding to different propagation phenomena observed in Fig. 8): the pulse located at *y*=0 µm corresponds to the field guided by the linear slab; the pulse located at *y*=- 35 µm comes from the first soliton emitted at approximately *z*=400 λ_{0} (as seen in Fig. 8); while the pulse located around *y*=- 20 µm corresponds to the second soliton formed at *z*=1200 λ_{0}. The first pulse is not affected by numerical dispersion, as its plane-wave spectrum is centered at 0 where, as shown by Fig.7(a), numerical dispersion is negligible. This is the reason why all the results in Fig. 9 are almost the same for this pulse. On the other hand, the plane-wave spectrum of the two solitons is passband-type as a consequence of their propagation angle, but as their propagation angle is similar, they will suffer from the same amount of dispersion phase error per unit length. However, as the first soliton formed at *z*=400 λ_{0} suffers an accumulated phase error which is much higher than the second one generated at *z*=1200 λ_{0}, the total error for the first pulse is much more sensitive to mesh size settings. This explains why the error in Fig. 9 simulations with the coarser discretization mesh is higher for the first emitted soliton than for the second one.

Numerical attenuation curves calculated in the interest frequency band show that total losses are under er 2×10^{-4} dB, so it can be stated that losses do not affect the simulation results.

## 4. Conclusions

In this paper, analytical expressions to calculate numerical dispersion and losses of the RK-FDBPM technique have been presented for the first time. These expressions make it possible to calculate the phase error and attenuation per unit length due to the numerical discretization as a function of mesh sizes. This way, a powerful method for analyzing the numerical technique performance as a function of the numerical parameters is provided. Phase error curves provide numerical phase distortion dependence. Attenuation constant curves provide stability characterization. This analysis made it possible to determine a criterion for the limitation of the cumulative phase error, which has shown to be useful for the numerical technique’s parameters setup. The criterion involves first estimating the spatial bandwidth of the optical envelope where 90% of the spectrum energy is located, and second making use of dispersion relations to set up the mesh steps to have a negative attenuation curve and an accumulated phase error which is significantly less than 2π.

This strategy has shown to be a very useful tool for the use of the RK-FDBPM method in the design of integrated optical devices, as it provides a good starting point for setting up the numerical method’s parameters, while avoiding the use of slow and inefficient trial and error procedures. The proposed tools have been tested in the analysis of linear and non-linear optical devices, in 2D and in 3D cases, with excellent results. It has been demonstrated that the numerical phase error can significantly modify field behavior in linear devices, where it can inhibit the guidance phenomena, and also in non-linear devices, where the equilibrium between chromatic dispersion and non-linear effect can be altered.

## Acknowledgments

This work was supported by Spanish National Research Project TEC2006-02868, the Andalusian Regional Ministry of Science, Innovation and Business (CICYE) Excellence Research Project TIC-02946 and the “Juan de la Cierva” National Fellowship program.

## References and Links

**1. **A. Koster, E. Cassan, S. Laval, L. Vivien, and D. Pascal, “Ultracompact splitter for submicrometer silicon-on-insulator rib waveguides,” J. Opt. Soc. Am. A **21**, 2180–2185 (2004). [CrossRef]

**2. **Z. Chen, Z. Li, and B. Li, “A 2-to-4 decoder switch in SiGe/Si multimode interference,” Opt. Express **14**, 2671–2678 (2006). [CrossRef] [PubMed]

**3. **J. Yamauchi, K. Sumida, and H. Nakano, “Analysis of a polarization splitter with a multilayer filter using a padé-operator-based power-conserving fourth-order accurate beam-propagation method,” IEEE Photon. Technol. Lett. **18**, 1858–1860 (2006). [CrossRef]

**4. **D. Dai, J. He, and S. He, “Compact silicon-on-insulator-based multimode interference coupler with bilevel taper structure,” Appl. Opt. **44**, 5036–5041 (2005). [CrossRef] [PubMed]

**5. **S. T. Lee, C. E. Png, F. Y. Gardes, and G. T. Reed “Optically switched arrayed waveguide gratings using phase modulation,” IEEE J. Lightwave Technol. **18**, 1858–1860 (2006).

**6. **M. Takenaka and Y. Nakano, “Multimode interference bistable laser diode,” IEEE Photon. Technol. Lett. **15**, 1035–1037, (2003). [CrossRef]

**7. **M. Raburn, M. Takenaka, K. Takeda, X. Song, J. S. Barton, and Y. Nakano, “Integrable multimode interference distributed Bragg reflector laser all-optical flip flop,” IEEE Photon. Technol. Lett. **18**, 1421–1423 (2006). [CrossRef]

**8. **N. N. Elkin, A. P. Napartovich, V. N. Troschieva, D. V. Vysotsky, T. Lee, S. C. Hagness, N. Kim, L. Bao, and J.L Mawst, “Antiresonant reflecting optical waveguide-type vertical-cavity surface emitting lasers: Comparison of full-vectorial finite difference time domain and 3D bidirectional beam propagation methods,” IEEE J. Lightwave Technol. **24**1834–1842 (2006). [CrossRef]

**9. **M. Takenaka and Y. Nakano, “Simulation of all-optical flip-flops based on bistable laser diodes with nonlinear couplers,” in Proceedings of the 4^{th} Int. Conf. on Numerical Simulation of Optoelectronic Devices (NUSOD ‘04),15–18 (2004).

**10. **C. Ma and E. Van Keuren, “A three-dimensional wide-angle BPM for optical waveguide structures,” Opt. Express **15**, 402–407 (2007). [CrossRef] [PubMed]

**11. **J. Shibayama, T. Takahashi, J. Yamauchi, and H. Nakano, “A three-dimensional multistep horizontally wide-angle beam-propagation method based on the generalized Douglas scheme,” IEEE Photon. Technol. Lett. **18**, 2535–2537 (2006). [CrossRef]

**12. **S. Sujecki “Wide-angle, finite-difference beam propagation in oblique coordinate system,” J. Opt. Soc. Am. A **25**, 138–145 (2008). [CrossRef]

**13. **J. de-Oliva-Rubio and I. Molina-Fernández, “Fast semivectorial non-linear finite-difference beampropagation method,” Microwave Opt. Technol. Lett. **40**, 73–77 (2004). [CrossRef]

**14. **J. de-Oliva-Rubio, *Desarrollo y validación de técnicas de diferencias finitas para el análisis de dispositivos ópticos lineales y no-lineales*, Ph. D. Thesis (in Spanish), ISBN: 84-690-3321-2, (Universidad de Málaga, 2006).

**15. **W. P. Huang and C. L. Xu
, “Simulation of three-dimensional optical waveguides by a full-vector beam propagation method,” IEEE J. Quantum Electron. **29**, 2639–2649 (1993). [CrossRef]

**16. **M. Matsuhara, “A novel beam propagation method based on the Galerkin method,” Electron. Commun. Jpn. Pt. II **73**, 41–47 (1990). [CrossRef]

**17. **B. E. A. Saleh and M. C. Teich, *Fundamentals of photonics, 2 ^{nd} Ed.* (John Wiley & Sons, 2007), Chap. 4, pp. 105–112.

**18. **I. Molina Fernández, C. Camacho Peñalosa, and J. I. Ramos, “Application of the two-dimensional Fourier transform to nonlinear wave propagation phenomena,” IEEE Trans. Microwave Theory Tech. **42**, 1079–1085 (1994). [CrossRef]

**19. **C. Vassallo, *Optical waveguides concepts, ser. Optical waveguides sciences and technology* (Elsevier, Amsterdam, 1991), Chap. 1, pp. 18–26.

**20. **I. Molina-Fernández, J. G. Wangüemert-Pérez, A. Ortega-Moñux, R. G. Bosisio, and KeWu, “Planar lightwave circuit six-port technique for optical measurements and characterizations,” IEEE J. Lightwave Technol. **23**, 2148–2157 (2005). [CrossRef]

**21. **K. Hayata, A. Misawa, and M. Koshiba, “Spatial polarization instabilities due to transverse effects in nonlinear guided-wave systems,” J Opt Soc Am B **7**, 1268–1280 (1990). [CrossRef]