1. Introduction

An efficient computational method for simulation of optical waveguide structures is of fundamental importance in the design of new optical devices. Since a certain number of grid points are needed for each wavelength, a direct numerical approach by the finite difference (FD) method or the finite element (FE) method that discretize the whole wave-guiding structure is prohibitively expensive [1–3]. For the slowly varying waveguides or piecewise invariant structures, some efficient numerical methods, such as the beam propagation met hod (BPM) [4–6] and the operator marching method (OMM) [7–9], have been extremely successful in computation of the forward propagating optical fields.

However, many practical optical devices, such as laser facets, grating structures, multiple dielectric layers or junctions of different waveguides, contain refractive index discontinuities along the direction of propagation of the optical field. Therefore, a bidirectional technique which considers coupling of the forward and backward waves is needed for accurate modeling of reflection or transmission [10–13]. The OMM, based on the Dirichlet-to-Neumann (DtN) mapping with the exact reformulations for pairs of operators, is used to solve the propagation problem in three dimensional optical waveguides that is formulated in a domain with just one direction (i.e. $z$) having a particularly large length, due to its excellent simulation speed with sufficient accuracy. In [14], a three-dimensional OMM of the scalar Helmholtz equation in the three-dimensional rectangular box with the radiation boundary was developed. In [15], a three-dimensional adjoint-based OMM (AOMM) of the complex Helmholtz-PML equation for long-range inhomogeneous optical waveguide devices was proposed. Hence, there is an understandable desire to incorporate the bidirectional technique and OMM into developing the so-called bidirectional OMM (Bi-OMM) in three dimensional optical waveguides since it can simultaneously consider the forward and backward waves.

In this paper, we propose a new Bi-OMM based on the bidirectional DtN mapping for three-dimensional optical waveguide structures. The bidirectional reflection operator is firstly determined, subsequently the exact one-way re-formulation is constructed through the bidirectional DtN map [7,16,17], and finally the local base transformation is implemented by the series-expansion technique [18,19]. To validate the present approach, a typical buried waveguide and the surface plasmon stripe waveguide with complex varying refractive indices are considered as the numerical examples [5,20]. Results show that the proposed Bi-OMM can deal with evanescent waves properly with accuracy comparable to that of the well-established bidirectional beam propagation method (Bi-BPM) [4,10]. And in the inhomogeneous waveguides, our proposed Bi-OMM also shows more effective with comparison of the AOMM.

The rest of this paper is organized as follows. Section 2 presents the basic mathematical model in unbounded waveguides. In Section 3 and 4 explain the formulation of the bidirectional operator marching scheme in detail, which is the main goal of this paper. To show the effectiveness and usefulness of the proposed approach, some numerical results are shown in Section 5. The conclusion is drawn in Section 6.

2. Basic problem

Considering the optical waves propagating in a three-dimensional waveguide structure, the propagation direction is supposed to be $z$-direction and the transverse cross-section is set to be $xy$-plane. In the Bi-OMM, the longitudinally varying integrated optical waveguide is first discretized and replaced by a suitable number of longitudinally invariant inhomogeneous segments. Then, by utilizing the governing wave equations in frequency-domain for a monochromatic wave, the electric field satisfies the following three-dimensional complex Helmholtz equation

The transverse operator on the $xy$ plane is defined as

At $z =0$, the incident wave is given as

And an outgoing condition can be written down at $z=L$,

Therefore, the propagation problem is formulated in a domain with just the main propagation direction (i.e., $z$ direction), having a particularly large length, and a relatively small transverse region.

3. Bidirectional operator marching

The total wave field can be split into right and left propagating components, denoted by $E_R$ and $E_L$, each of which satisfies a one-way re-formulation based on the DtN map. We write

Substituting Eq. (8) into the Eq. (1), we have

Since the above formulas Eqs. (10) and (11) are true for any solution of Eq. (1) (satisfying the boundary conditions at (3), and the outgoing condition at $z = L$), we obtain

The above Eqs. (12) and (13) should be solved for decreasing $z$ with an exact radiation condition at $+\infty$. Then, the fundamental solution operators $G_1$ and $G_2$ have been introduced for stability reasons,

The initial condition at $z=L$ is $G_1(L) = I$ and $G_2(L) = I$, where $I$ is the identity operator. When $G_1(0)$ and $G_2(0)$ are calculated, the solution at $z =L$ can be generated by the multiplications with the starting field $E_R(x, y, 0)$ and $E_L(x, y, 0)$, namely

The formulation above allows us to evaluate $E(x, y, L)$, given $E(x, y, 0)$. However, typically we are given $E_R(x, y, 0)$, and should determine $E_L(x, y, 0)$ such that the boundary condition that no backward field will be present at the output, i.e., $E_L(x, y, L) = 0$, is satisfied.

Now, we consider the multiple-layered structure shown in Fig. 1, which consists of $N$ interfaces, separating $N-1$ transversely inhomogeneous layers. Let the discrete points $z_i$ ($1\leq i < N$) satisfy

Each interval $(z_i, z_{i+1})$ corresponds to a range-independent segment and $h = L/N$.

 

Fig. 1. Propagation along the transition region of $N$ transverse inhomogeneous layers.

Download Full Size | PDF

The approximated Helmholtz equation by piecewise range-independent segments using the medium values at $(z_i, z_{i+1})$ is

On the interval $(z_0, z_1)$, the Eqs. (12) and (13) are approximated by

Then the total wave field can be decomposed into the right and left going waves, satisfying

The solution to this general problem can be written symbolically, in terms of the introduced reflection operators $\mathscr {R}$ for right and left propagating incident fields. At the point $z_i$, we define the reflection operator $\mathscr {R}_i$ by

And it is easy to obtain from $E_z(x, y, z_1) = Q_1^{(1)} E_R + Q_2^{(1)} E_L$ and the above relationship that

Similarly, at the point $z_0$, we can obtain

Compare the definition of $\mathscr {R}_1$ and $\mathscr {R}_0$ and notice that

The formulas (24), (25) and (27) reveal the exact relationship between $Q_1^{(0)}$, $Q_2^{(0)}$ and $Q_1^{(1)}$, $Q_2^{(1)}$ for Eqs. (19) and (20). In the DtN formulation above, the continuity condition is implicit in the requirement that $Q_1^{(1)}$ and $Q_2^{(1)}$ obtained from the previous calculation in $(z_1, z_2)$ be the same $Q_1^{(1)}$ and $Q_2^{(1)}$ used to calculate $Q_1^{(0)}$ and $Q_2^{(0)}$ in $(z_0, z_1)$.

Therefore, on the interval $(z_0, z_1)$, we have

At this point, we know how to construct the bidirectional DtN operators $Q_1$ and $Q_2$, and the fundamental solution operators $G_1$ and $G_2$. Meanwhile, solving such operators from $z=L$ to $z =0$ with the formulas (24), (25) and (27) are used in the bidirectional marching step from $z_{i+1}$ to $z_i$, where $1\leq i\leq N$.

4. Local base transformation and iteration schemes

In this section, a new treatment for local base transformation is introduced [21], then the above bidirectional operator marching formulas are considered into numerical implementation through the matrices.

The eigenvalue problem of the transverse operator $\mathscr {B}$ is important, and it is defined as

We will also take the segment $(z_0, z_1)$ for example. For integers $n_x$ and $n_y$, the transverse operator $\mathscr {B}$ are expanded by bivariate Chebyshev coefficients such that

The derivatives of two transverse directions $x$ and $y$ can be expressed in terms of the grid values $\mathscr {B}(x(\xi _i),y(\xi _j),z_{1/2})$. That is, the $p$th derivatives

In consideration of the boundary conditions (3), we can make use of the $(n_x-1) \times (n_x-1)$ matrix ${\mathscr {D}}^2_{n_x-1}$ obtained by stripping $\tilde {\mathscr {D}}_{n_x}^2$ of its first and last rows and columns. In MATLAB notation:

Suppose $\nu (x, y)_m$ and $\hat {\nu }(x,y)_m$ are the normalized eigenvectors of $\mathscr {M}$ and $\mathscr {M}^{\textrm{H}}$, respectively, and assume their corresponding eigenvalues are $\lambda _m$ and $\bar {\lambda }_m$, where $\mathscr {M}^{\textrm{H}}$ is the conjugate transverse of $\mathscr {M}$. Then

Otherwise, the eigenvector decomposition of $\mathscr {M}_{1/2}$ and $\mathscr {M}^{\textrm{H}}_{1/2}$ are

Thus the bidirectional operator marching algorithm (Bi-OMM) can be outlined as follows:

(1) $E_R^{(0)}$ is given, $E_L^{(0)}$ is guessed (e.g., set to 0) and choose a step size $h >0$;

(2) for $i = 0, 1, 2, \ldots$ do:

a) $\mathscr {O}_1 = \mathscr {N}\hat {\mathscr {O}}_1\mathscr {N}^{-1}$,

$\mathscr {O}_2 = \mathscr {N}\hat {\mathscr {O}}_2\mathscr {N}^{-1},$

$\mathscr {P}_1 := ( \textrm{i} \sqrt { \Lambda _{1/2, k} } + \mathscr {O}_2 )^{-1}(\textrm{i}\sqrt {\Lambda _{1/2, k}}-\mathscr {O}_1),$

$\mathscr {P}_0 = e^{\textrm{i}h\sqrt {\Lambda _{1/2, k}}}\mathscr {P}_1e^{\textrm{i}h\sqrt {\Lambda _{1/2, k}}},$

$\mathscr {S}_0 := \mathscr {N}\mathscr {S}_1\mathscr {N}^{-1}(I + \mathscr {P}_1)e^{\textrm{i}h\sqrt {\Lambda _{1/2, k}}}(I + \mathscr {P}_0)^{-1},$

$\left ( {\begin {array}{c} {{E_R}^{(i)}}\\ {{E_L}^{(i)}} \end {array}} \right ) = \hat {\nu }_{1/2}\mathscr {S}_0\hat {\nu }^{\textrm{H}}_{1/2}\left ( {\begin {array}{c} {{E_R}^{(0)}}\\ {{E_L}^{(i)}} \end {array}} \right )$;

b) if $\left | {{E_L}^{(k)}} \right | < 1e - 6$, stop; otherwise,

c) $\hat {\mathscr {O}}_1 = \mathscr {N}^{-1}\mathscr {O}_1\mathscr {N}$,

$\hat {\mathscr {O}}_2 = \mathscr {N}^{-1}\mathscr {O}_2\mathscr {N}$

$\mathscr {P}_0 := ( \textrm{i} \sqrt { \Lambda _{1/2, k} } + \hat {\mathscr {O}}_2 )^{-1}(\textrm{i}\sqrt {\Lambda _{1/2, k}}-\hat {\mathscr {O}}_1),$

$\mathscr {P}_1 = e^{-\textrm{i}h\sqrt {\Lambda _{1/2, k}}}\mathscr {P}_0e^{-\textrm{i}h\sqrt {\Lambda _{1/2, k}}},$

$\mathscr {S}_1 := \mathscr {N}^{-1}\mathscr {S}_0\mathscr {N}(I + \mathscr {P}_0)e^{-\textrm{i}h\sqrt {\Lambda _{1/2, k}}}(I + \mathscr {P}_1)^{-1},$

$\left ( {\begin {array}{c} {*}\\ {{E_L}^{(i+1)}} \end {array}} \right ) = \hat {\nu }^{\textrm{H}}_{3/2}\mathscr {S}_1\hat {\nu }_{3/2}\left ( {\begin {array}{c} {{E_R}^{(i)}}\\ {\alpha \cdot {E_L}^{(i)}} \end {array}} \right )$, where $0 < \alpha <1$ and * indicates an unused component.

In brief, at each iteration, the guess for $E_L^{(i)}$ is propagated along with $E_R^{(0)}$ forward through the system. The fields $E_R^{(i)}$ and $E_L^{(i)}$ are propagated backward through the structure to obtain a new guess for $E_R^{(i+1)}$. Then, $E_L^{(i)}$ should be zero when we reach a solution. The damping factor, $\alpha$, is selected to optimize convergence. Compared with the Bi-BPM, our method is developed from the exact formulation of pairs of operators, and no paraxial or slowly varying envelope approximation has been used as in the BPM.

5. Numerical results

Two numerical examples will be presented here to illustrate the validity of the proposed Bi-OMM. The relative error (RE) $\varepsilon$ using the discrete $L^2$ norm is defined:

5.1 Buried waveguide

To validate the proposed Bi-OMM, we first apply it to analyze a typical Buried waveguide, as shown in Fig. 2 [4,15]. The core of the buried waveguide is assumed to have the reflective index $n_2 = 3.54$. The reflective index of the cladding material is $n_1 = 3.17$. The core thickness is fixed to $d = 0.4$ $\mu$m, while the core width is set to be $\omega = 0.4$ $\mu$m. The operating wavelength is assumed to be $\lambda = 1.30$ $\mu$m. The propagation distance is $L = 10$ $\mu$m. The width of the transverse computational window is $D_1 = D_2 = 1$ $\mu$m. Meanwhile, the transverse operator is approximated with a grid of $20\times 20$ and $k = 8$ in the local subtraction of the eigenmodes.

 

Fig. 2. Buried waveguide.

Download Full Size | PDF

The starting field at $z=0$ is $E(x, y, 0) =\sin (2\pi x)\sin (2\pi y)$, which is the third propagating mode, and this case has an analytical solution

In Fig. 3, the numerical fields $EL_R\approx E_R(x, y, L)$ with the step length $h = 1$ $\mu$m are plotted. Meanwhile, in Fig. 4, the RE is shown for the numerical fields $EL_R$ with $h = 1$ $\mu$m compared to the exact fields for a buried waveguide in the $XY$ plane. In contrast, the RE for the numerical fields and the exact solution is small enough, and the accuracy of the proposed Bi-OMM is pretty good in the buried waveguide. Thus good solutions are already obtained with $h =1$ $\mu$m. In Table 1, we list the relative errors ( compared with the "exact" solution) of the numerical solutions of $E(x, y, L)$ with different choices of the range step size.

 

Fig. 3. Numerical fields $EL_R$ for a buried waveguide with $h=1$ $\mu$m.

Download Full Size | PDF

 

Fig. 4. RE on the numerical fields $EL_R$ with $h = 1$ $\mu$m compared to the exact fields for a buried waveguide in the $XY$ plane.

Download Full Size | PDF

Table 1. RE for different step sizes by Bi-OMM in the buried waveguide.

View Table | View all tables in this article

For this case, we also compare the efficiency of the Bi-OMM, FD-Bi-BPM, and AOMM for the buried waveguide by finding the number of nodes, CPU time and the required computing memory for two predetermined $\varepsilon$ error goals of $0.01$ and $0.001$. The results in Table 2 show the obvious superiority of the proposed Bi-OMM with the much lower computational cost over the FD-Bi-BPM and the results are especially impressive for the higher accuracy goal.For the buried waveguide, our proposed Bi-OMM can achieve similar accuracy as the adjoint-based OMM (AOMM), meanwhile, avoiding the construction of the adjoint operator for the transverse operator that usually is complicated. When the transverse operator is approximated with a grid of $40 \times 40$, and similarly $k=8$ for the eigenmode local subtraction, the effectiveness of the proposed Bi-OMM in the buried waveguide has no significant change. It is mainly because the differentiation matrices Eq. (37) are dense, and the series-expansion method is a powerful mode solver, with its low memory storage requirements for a small matrix, could obtain highly accurate results [18,22].

Table 2. Comparison of the computational cost of different methods for a fixed error goal with $h/4$ in the buried waveguide.

View Table | View all tables in this article

5.2 Surface plasmon stripe waveguide

To further confirm the applicability of the proposed Bi-OMM to plasmonic devices, we consider a surface plasmon stripe waveguide structure (as shown in Fig. 5) and discuss the numerical accuracy of Bi-OMM in lossy and leaky inhomogeneous waveguide analysis. The refractive indices of the cladding (air), core and substrate, are $n_0 = 1.0$, $n_1 = (-26.1437 - \textrm{i}1.8479) [1 + 0.05e^{-(z/L-0.5)^2}\sin ^2(\pi x)\sin ^2(\pi y)]$ and $n_2 = 1.46^2$. The core height is fixed to $d = 55$ nm, while the core width is set to be $\omega = 3.3$ $\mu$m. The propagation distance is set to be $L = 2$ mm, and $D_1 = D_2 = 10$ $\mu$m. The operating wavelength is set to be $0.8$ $\mu$m. Meanwhile, the transverse operator is discretized by a grid of $24\times 24$ and $k = 10$ in the local subtraction of the eigenmodes. The starting field at $z = 0$ is set to be $E(x, y, 0) = \sin (2\pi x)\sin (2\pi y)e^{-x-y}.$

 

Fig. 5. Surface plasmon stripe waveguide.

Download Full Size | PDF

The numerical fields calculated by our Bi-OMM with $h = 60$ $\mu$m are shown in Fig. 6. The RE on the numerical fields $EL_R$ with $h = 60$ $\mu$m compared to a much more accurate solution obtained with $h/32$ for a surface plasmon stripe waveguide in the $XY$ plane are thoroughly presented in Fig. 7. We can see from the RE of Fig. 7 that the results are still satisfactory and a reasonably good solution is already obtained with $h$. As above, we list the relative errors of the numerical solutions of $EL_R$ with different choices of the range step size in Table 3. And the numerical fields with a step length $h/128$ are used for reference, because analytic solutions are not available here.

 

Fig. 6. Numerical fields $EL_R$ for a surface plasmon stripe waveguide with $h = 60$ $\mu$m.

Download Full Size | PDF

 

Fig. 7. RE on the numerical fields $EL_R$ with $h = 60$ $\mu$m compared to the reference fields with $h/32$ for a surface plasmon stripe waveguide in the $XY$ plane.

Download Full Size | PDF

Table 3. RE for different step sizes by Bi-OMM in the surface plasmon stripe waveguide.

View Table | View all tables in this article

A comparison of the computational cost of the proposed Bi-OMM, FD-Bi-BPM and AOMM for the surface plasmon stripe waveguide with two predetermined $\varepsilon$ error goals of $0.01$ and $0.001$ is also shown in Table 4. The proposed Bi-OMM has much lower computational cost compared to the FD-Bi-BPM, the required computer memory and CPU time needed are also much smaller. For the varying refractive index profiles, it is even harder than the constant refractive index profile case of the Bi-BPM to do rational approximations in each segment. For the surface plasmon stripe waveguide with the varying refractive index profiles, our proposed Bi-OMM is more efficient than the AOMM, and what is more, the backward waves can also be considered by the Bi-OMM. When the transverse operator is approximated with a grid of $48 \times 48$, and similarly $k = 10$ for the eigenmode local subtraction, the effectiveness of the proposed Bi-OMM in a surface plasmon stripe waveguide is still satisfactory as above.

Table 4. Comparison of the computational cost of different methods for a fixed error goal with $h/16$ in the surface plasmon stripe waveguide.

View Table | View all tables in this article

6. Conclusion

In this paper, we have developed an efficient Bi-OMM based on the bidirectional DtN mapping and the large range step method for three-dimensional optical waveguide structures. Using a traditional buried waveguide example, the technique proves to be very accurate compared to the analytical result. Numerical results also show that the present Bi-OMM can deal with inhomogeneous waveguide structures properly with accuracy comparable to the well-established FD-Bi-BPM and less demanding in terms of memory and CPU-time resources. Currently, work is being done to extend the proposed method to the full-vectorial formulation.