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_xxô+e_yŷ is the slowly varying complex envelope of the electric field ε, which can be written as

where k₀ is the vacuum propagation constant and n_ref is the BPM reference refractive index.

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_x is the number of mesh points in the computational window on the 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]

where N_x and N_y are the number of samples of the computational window in the x and y direction, respectively. The size of ${\overline{\overline{M}}}_T$ is N_xN_y×N_xN_y, for linear homogeneous media it is symmetric and constant, meaning that all the elements in each subdiagonal are identical, so ${\overline{\overline{M}}}_T$(k,k+p)=${\overline{\overline{M}}}_T$(k,k-p), ∀k∈[1,N_xN_y]. Each coefficient ${\overline{\overline{M}}}_{T,i}^n$ in Eq. (16) is any of the i-th subdiagonal elements in(M̿_T)ⁿ, where M̿ⁿ _T,0 represents the diagonal elements.

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

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 λ₀=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 λ₀. Mesh steps are Δy=0.5 µm, Δz=0.1 λ₀ in the case of (a)-mesh, Δy=0.1 µm, Δz=0.1 λ₀ for (b)-mesh, and Δy=0.025 µm, Δz=0.0125 λ₀ 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₀ 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.

 

Fig. 1. Normalized phase error for the different meshes. (a)-mesh Δy=0.5 µm, Δz=0.1 λ₀. (b)-mesh Δy=0.1 µm, Δz=0.1 λ₀. (c)-mesh Δy=0.025 µm, Δz=0.0125 λ₀. The superimposed black line is the normalized envelope spatial spectrum.

Download Full Size | PDF

In Fig. 2 the propagation of TE₀ 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₂₁ 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].

 

Fig. 2. Propagation of the TE₀ excitation in the tilted slab. In Fig. 2(b) the (b)-mesh and the (c)-mesh solutions are indistinguishable.

Download Full Size | PDF

Figure 3 shows the z evolution of S₂₁, 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₂₁ 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₂₁ modulus dramatically decreases. In Fig. 3(b) one can see how the S₂₁ phase error agrees well with the values calculated in the dispersion analysis.

 

Fig. 3. Propagating and fundamental mode fields overlap integral: S₂₁. Fig. 3(a): S₂₁ amplitude. Fig. 3(b): S₂₁ phase error regarding fundamental mode phase at d=z/cosθ. The blue line corresponds to the (a)-mesh, the red line to the (b)-mesh and the green line to the (c)-mesh.

Download Full Size | PDF

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.4549 embedded in a medium with index n₂=1.4440. At the design wavelength (λ₀=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],

 

Fig. 4. MMI 2×3 coupler geometry

Download Full Size | PDF

 

Fig. 5. Variations in RK-FDBPM phase error depending on mesh size. n_ref=1,4487, n=1.4549, λ₀=1.55µm, L=3570µm.

Download Full Size | PDF

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,β_y∈[-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.

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.

 

Fig. 6. Longitudinal distribution of the electric field for n_ref=1,4487, λ₀=1.55µm. (a) Δx=Δy=1µm, Δz=0.1 µm. (b) Δx=Δy=0.1µm., Δz=0.01 µm.

Download Full Size | PDF

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₀ and TM₀ 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 ² _NL=n ² _s+α_s(|e_x|²+|e_y|²+|e_z|²), and the non-linear coefficien0t is α_s=10^-10 V ²/m ². The guiding structure is excited simultaneously with the linear slab TE₀ and TM₀ modes at 633nm wavelength. The excitation power is 2.92 mW/mm, equally divided among both modes. The propagation length is L=2400λ ₀.

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λ ₀. 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λ₀ and z=1200λ₀ 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.

 

Fig. 7. Numerical analysis results for the non-linear slab characterization. (a) Variation of phase error for different mesh step sizes. (b) Variation with Δz of the numerical attenuation constant for the mesh with Δy=0.25 µm. For those Δz with positive attenuation in the bandwidth ends (red and blue curves in (b)), the method becomes unstable.

Download Full Size | PDF

 

Fig. 8. Electric field in the non-linear slab with Δy=0.5 µm and Δz=8λ₀. It can be seen how two solitons are generated at z=400λ₀ y z=1200λ₀. They propagate in the non-linear substrate (y<-10µm).

Download Full Size | PDF

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λ ₀. 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λ ₀, resulting in the Fig. 7(b) red curve, which still yields instability problems. Finally, with Δz=2λ ₀ 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=λ ₀/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λ₀ 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.

 

Fig. 9. Field in the non-linear slab at z=2400λ₀ for the different proposed meshes. Field strength is normalized with respect to its maximum value in the structure.

Download Full Size | PDF

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 λ₀ (as seen in Fig. 8); while the pulse located around y=- 20 µm corresponds to the second soliton formed at z=1200 λ₀. 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 λ₀ suffers an accumulated phase error which is much higher than the second one generated at z=1200 λ₀, 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.