Abstract

A novel bidirectional operator marching method based on the Dirichlet-to-Neumann (DtN) mapping for three-dimensional optical waveguide structures is developed and implemented using iterative methods. The backward propagation wave is integrated into the classical operator marching approach which represents the forward propagating wave. The bidirectional range marching formulas are exact for each range-independent piece and a large range step is possible in both directions. The validity and effectiveness of our proposed method are verified by analyzing uniform waveguides and longitudinal waveguides with varying refractive indices.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

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 [13]. For the slowly varying waveguides or piecewise invariant structures, some efficient numerical methods, such as the beam propagation met hod (BPM) [46] and the operator marching method (OMM) [79], 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 [1013]. 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

$$\frac{\partial^2 E}{{\partial x^2}}+\frac{\partial^2 E}{{\partial y^2}}+\frac{\partial^2 E}{{\partial z^2}}+\kappa_0^2n^2(x,y,z)E=0,$$
where $\kappa _0 = 2\pi /\lambda$ is the wavenumber in free space, $n(x,y,z)$ is the complex refractive index profiles and the spatially dependent wavenumber can be given by $\kappa _0n(x, y, z)$. Here $n(x,y,z)$ is $z$-invariant for $z<0$ ($n(x,y,z) = n_{0}(x,y)$) and $z > L$ ($n(x,y,z) = n_{\infty }(x,y)$). Suppose $\textrm{Re}(\kappa _0n(x,y,z)) >0, \textrm{Im}(\kappa _0n(x,y,z)) >0$, $\textrm{Im}(\kappa _0n(x,y,z))\ll \textrm{Re}(\kappa _0n(x,y,z))$, and $\frac {1}{{\textrm{Re}(\kappa _0n(x,y,z))}}\ll 1 \ll L$, and the distinct length scale shows the wave propagates in the weakly lossy medium over a relatively large range distance. In $\Omega = \{ (x, y) \in [0, D_1] \times [0, D_2])\}$, $n(x, y, z)$ is a continuous or piecewise continuous function.

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

$$\mathscr{B}(x, y) =\frac{\partial^2}{{\partial x^2}} +\frac{\partial^2}{{\partial y^2}}+\kappa_0^2n^2(x,y,z),$$
and a Dirichlet boundary condition is imposed on the $xy$ plane:
$$\begin{aligned}&E(0^{+},y,z) = 0,\quad E(D_1^{-},y,z) = 0, \\ &E(x,0^{+},z) = 0,\quad E(x,D_2^{-},z) = 0. \end{aligned}$$

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

$$E(x, y, 0) = E_0 (x, y).$$

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

$$E_z =\textrm{i}\sqrt{ \frac{\partial^2}{{\partial x^2}} +\frac{\partial^2}{{\partial y^2}} + \kappa_0^2 n^2(x,y,L)}E, \quad z = L,$$
which represents the outgoing waves as $z \to \infty$ and $\textrm{i} = \sqrt {-1}$.

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

$$\left(\partial_z +\textrm{i}\sqrt{\mathscr{B}}\right)E_R=0,$$
and
$$\left(\partial_z -\textrm{i}\sqrt{\mathscr{B}}\right)E_L=0,$$
the right and left propagating fields in each $z$-invariant segment are, respectively, identified with the two formal solutions to Eqs. (6) and (7). The electric field at any point along the propagation direction $z$ can be represented by the right ("$E_R$") and left ("$E_L$") waves. And it follows that the wave field splitting in a transversely inhomogeneous environment,
$$\partial_z E = Q_1 E_R+ Q_2 E_L,$$
is accomplished by the DtN map, where $Q_1$ and $Q_2$ are the two operator solutions of the quadratic operator equation
$$Q^2-\mathscr{B}^2 =0.$$

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

$$\left[\frac{\textrm{d}Q_1}{{\textrm{d}z}}+Q_1^2 + \partial_x^2 + \partial_y^2 + \kappa_0^2n^2(x,y,z)\right]E_R=0,$$
and
$$\left[\frac{\textrm{d}Q_2}{{\textrm{d}z}}+Q_2^2 + \partial_x^2 + \partial_y^2 + \kappa_0^2n^2(x,y,z)\right]E_L=0.$$

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

$$\frac{dQ_1}{{dz}} ={-}\left[\partial^2 x +\partial^2 y + \kappa_0^2n^2(x,y,z)\right]-Q_1^2,$$
and
$$\frac{dQ_2}{{dz}} ={-}\left[\partial^2 x +\partial^2 y + \kappa_0^2n^2(x,y,z)\right]-Q_2^2.$$

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,

$$G_1E_R(x, y, 0) = E_R(x, y, L),\quad \frac{\textrm{d}G_1}{{dz}}={-}G_1Q_1,$$
and
$$G_2E_L(x, y, 0) = E_L(x, y, L),\quad \frac{\textrm{d}G_2}{{dz}}={-}G_2Q_2.$$

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

$$E_R(x, y, L) = G_1(0)E_R(x, y, 0),\quad E_L(x, y, L) =G_2(0)E_L(x, y, 0).$$

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

$$-\infty<0 = z_0 < z_1<\cdots<z_N = L <{+}\infty.$$

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

 figure: Fig. 1.

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

Download Full Size | PPT Slide | PDF

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

$$E_{xx} + E_{yy} + E_{zz} + \kappa^2n^2(x, y, z_{i+1/2})E =0,$$
where $z_{i+1/2} = z_i + h/2 = (z_i + z_{i+1})/2$. Exact solutions in each segment can be written down and continuity conditions at $z_{i+1/2}$ are used to match the solutions. We take $i=0$, for example, to describe the bidirectional marching form $z_0$ to $z_1$.

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

$$Q_1' ={-}Q_1^2 - \mathscr{B}_{1/2},$$
and
$$Q_2'={-}Q_2^2 -\mathscr{B}_{1/2},$$
where $\mathscr {B}_{1/2}=\mathscr {B}({z_{1/2}})$. Given the initial conditions $Q_1^{(1)}\approx Q_1(z_1)$ and $Q_2^{(1)}\approx Q_2(z_1)$, the exact solutions of $Q_1^{(0)}$ and $Q_2^{(0)}$ are used to obtain $Q_1^{(0)}\approx Q_1(z_0)$ and $Q_2^{(0)}\approx Q_2(z_0)$.

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

$$\left(\partial z + \textrm{i}\sqrt{\mathscr{B}_{1/2}}\right)E_R =0,$$
and
$$\left(\partial z - \textrm{i}\sqrt{\mathscr{B}_{1/2}}\right)E_L =0.$$

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

$$E_L(x, y, z_i) = \mathscr{R}_i E_R(x, y, z_i).$$

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

$$\mathscr{R}_1 = \left[\textrm{i}\sqrt{\mathscr{B}_{1/2}}+ Q_2^{(1)}\right]^{{-}1}\left[\textrm{i}\sqrt{\mathscr{B}_{1/2}}- Q_1^{(1)}\right].$$

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

$$Q_1^{(0)} =\textrm{i}\sqrt{\mathscr{B}_{1/2}} - \left(\textrm{i}\sqrt{\mathscr{B}_{1/2}} + Q_2^{(0)}\right)\mathscr{R}_0.$$

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

$$E_R(x, y, z_1) = e^{\textrm{i}h\sqrt{\mathscr{B}_{1/2}}} E_R(x, y, z_0), \quad E_L(x, y, z_1) = e^{-\textrm{i}h\sqrt{\mathscr{B}_{1/2}}}E_L(x, y, z_0),$$
we can have
$$\mathscr{R}_0 = e^{\textrm{i}h \sqrt{\mathscr{B}_{1/2}}}\mathscr{R}_1 e^{\textrm{i}h \sqrt{\mathscr{B}_{1/2}}}.$$

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

$$G_1^{(0)} = G_1^{(1)} (I + \mathscr{R}_1) e^{\textrm{i}h\sqrt{\mathscr{B}_{1/2}}}\left(I +\mathscr{R }_0\right)^{{-}1},$$
and
$$G_2^{(0)} = G_2^{(1)} (I + \mathscr{R}_1) e^{\textrm{i}h\sqrt{\mathscr{B}_{1/2}}}\left(I +\mathscr{R}_0\right)^{{-}1},$$
where the fundamental solution operators $G_1^{(0)}\approx G_1(z_0)$, $G_2^{(0)}\approx G_2(z_0)$, $G_1^{(1)}\approx G_1(z_1)$ and $G_2^{(1)} \approx G_2(z_1)$.

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

$$\mathscr{B}(x, y) \nu(x, y) = \lambda \nu(x, y), \quad 0 < x < D_1, 0 < y < D_2,$$
corresponding to the boundary condition (3). The square root operator,
$$\mathscr{B}^{(1/2)}(x, y) = \sqrt{\frac{\partial^2}{{\partial x^2}} +\frac{\partial^2}{{\partial y^2}}+\kappa_0^2n^2(x,y,z)},$$
is defined to satisfy
$$\mathscr{B}^{(1/2)}\nu(x, y) = \sqrt{\lambda}\nu(x, y),\quad 0 < x < D_1, 0 < y < D_2,$$
where $v(x, y)$ and $\lambda$ are the same as above, and the branch of the square root is chosen as $(-\pi /2, \pi /2]$.

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

$$\mathscr{B}(x,y,z_{1/2}) \simeq \sum^{n_x}_{i=0}\sum^{n_y}_{j=0}\mathscr{B}(x(\xi_i),y(\xi_j), z_{1/2})T_i(\psi(y))T_j(\phi(x)),$$
where $\{(x(\xi _i),y(\xi _j)),0\leq i \leq n, 0\leq j \leq m\}$ denote the nodes in a tensor product grid of $n_x\times n_y$ nodes in two transverse directions on $\Omega$, $T_i(x) = \cos (i\cos ^{-1}x)$ for $x \in [-1, 1]$ is the degree $i$ Chebyshev polynomial ( of the first kind), and $\phi (x) = 2(x/D_1-1)$ and $\psi (y) = 2(y /D_2-1)$ are affine transformations from $[0, D_1]$ and $[0, D_2]$ to $[-1, 1]$, respectively [22,23].

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

$$\frac{\partial^p \mathscr{B}}{{\partial x^p}}|_{x=x(\xi_i),y=y(\xi_j), z=z_{1/2}}= \sum^{n_x}_{i=0}\sum^{n_y}_{j=0}u(x(\xi_i),y(\xi_j),z_{1/2})T_j^{(p)}(x)T_i(y).$$

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:

$$\mathscr{D}^2=\mathscr{D}^2_{n_x-1} = \tilde{\mathscr{D}}^2_{n_x}(1:n_x-1, 1:n_x-1),$$
where $\tilde {\mathscr {D}}$ is the first derivative Chebyshev differentiation matrix of size $(n_x +1) \times (n_x +1)$. Similarly
$$\mathscr{D}^2_{n_y-1} =\tilde{ \mathscr{D}}^2_{n_y}(1:n_y-1, 1:n_y-1).$$
Thus we can get the matrix form of the transverse operator $\mathscr {B}(z_{1/2})$ is
$$\mathscr{M}_{1/2} = I \otimes \mathscr{D}^2 + \mathscr{D}^2 \otimes I + K_{1/2} ,$$
where $n_x = n_y$, $\otimes$ is the kronecker product, $K_{1/2}=\kappa _0^2n^2(x(\xi _k),y(\xi _t),z_{1/2})I \otimes I$ and $I$ is the identity matrix.

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

$$\left\langle {\nu_n, \hat{\nu}_m} \right\rangle = \sigma_{n, m},$$
where $\left \langle {*,*}\right \rangle$ is the Hilbert space $L^2(\Omega )$ inner product $\left \langle {\phi ,\psi }\right \rangle = \int _\Omega \phi \bar {\psi }d\Omega$, and the overline indicates complex conjugation.

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

$$\mathscr{M}_{1/2}\nu_{1/2} = \nu_{1/2} \Lambda_{1/2},\quad \mathscr{M}^{\textrm{H}}_{1/2}\hat{\nu}_{1/2} = \hat{\nu}_{1/2}\Lambda_{1/2}^\textrm{H}.$$
Denote the discrete matrices of $\mathscr {B}_{3/2}$ as $\mathscr {M}_{3/2}$. And similarly we have
$$\mathscr{M}_{3/2}\nu_{3/2} = \nu_{3/2} \Lambda_{3/2},\quad \mathscr{M}^{\textrm{H}}_{3/2}\hat{\nu}_{3/2} = \hat{\nu}_{3/2}\Lambda_{3/2}^\textrm{H}.$$
Then, by Eq. (38), we have
$$\hat{\nu}_{1/2}\nu_{1/2} = I, \quad \hat{\nu}_{3/2}\nu_{3/2} = I.$$
Meanwhile, a locally truncated technique for local eigenfunction expansion can be used in the marching solver, which means only use the first $k$ eigenvectors with the largest real parts of $\nu _{1/2}$ and $\nu _{3/2}$to increase the computing efficiency. Then the local base expansion can be implemented as
$$\mathscr{N} = \hat{\nu}_{1/2, k}^{\textrm{H}} \nu_{3/2, k},\quad {\textrm{if}} \textrm{ }\nu_{3/2, k} =\nu_{1/2, k}\mathscr{N}.$$
The bidirectional one-way formulation above allows us to evaluate $E(x, y, L)$, given $E(x, y, 0)$. However, typically we are given $E_R^{(0)} \approx E_R(x, y, 0)$, and must determine $E_L^{(0)} \approx 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.

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:

$$\varepsilon : = \frac{||E_r - E||_2}{{||E_r||_2}}= \left( \frac {1}{{\mathscr{I}}}\sum_{k=1}^{\mathscr{I}}\left(\frac{E_r -E}{{E_r}}\right)^2\right)^{1/2},$$
where $\mathscr {I}$ is a discrete set of boundary or domain error test points, $E_r$ is the reference solution and $E$ is the numerical solution. All of the simulations were on a desktop with an Intel Core CPU i7-8550 at 1.80 GHz with 8 GB of RAM.

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.

 figure: Fig. 2.

Fig. 2. Buried waveguide.

Download Full Size | PPT Slide | 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

$$E(x, y, L) = e^{\textrm{i}\sqrt{\kappa^2_0- 8\pi^2}L}\sin(2\pi x)\sin(2\pi y).$$
The exact fields in Eq. (44) are used for reference.

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.

 figure: Fig. 3.

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

Download Full Size | PPT Slide | PDF

 figure: Fig. 4.

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 | PPT Slide | PDF

Tables Icon

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

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].

Tables Icon

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

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}.$

 figure: Fig. 5.

Fig. 5. Surface plasmon stripe waveguide.

Download Full Size | PPT Slide | 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.

 figure: Fig. 6.

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

Download Full Size | PPT Slide | PDF

 figure: Fig. 7.

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 | PPT Slide | PDF

Tables Icon

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

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.

Tables Icon

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

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.

Funding

Zhejiang Provincial Natural Science Foundation of China (LQ19A010012).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. A. W. Snyder and J. D. Love, Optical Waveguide Theory (Chapman and Hall, 1983).

2. R. Scarmozzino, A. Gopinath, R. Pregla, and S. Helfert, “Numerical technique for modeling guided-wave photonic devices,” IEEE J. Sel. Top. Quantum Electron. 6(1), 150–162 (2000). [CrossRef]  

3. M. Levy, Parabolic Equation Methods for Electromagnetic Wave Propagation (IET, 2009).

4. S. Obayya, Computational Photonics (Wiley, 2011).

5. S. Kawai, A. Iguchi, and Y. Tsuji, “Study on high precision and stable finite element beam propagation method based on incomplete third order hybrid edge/nodal element,” J. Lightwave Technol. 36(11), 2278–2285 (2018). [CrossRef]  

6. X. Ma, S. Zhu, L. Li, J. Li, X. Shao, and A. Galvanauskas, “Finite-difference beam-propagation method for anisotropic waveguides with torsional birefringence,” Opt. Express 26(4), 3995–4003 (2018). [CrossRef]  

7. Y. Y. Lu, “One - way large range step methods for Helmholtz waveguides,” J. Comput. Phys. 152(1), 231–250 (1999). [CrossRef]  

8. Y. Y. Lu, “Some techniques for computing wave propagation in optical waveguides,” Commun. Comput. Phys. 1(6), 1056–1075 (2006).

9. S. Zhao and J. Zhu, “Efficient marching treatment for optical propagation in the waveguide with a discontinuous interface,” J. Opt. Soc. Am. B 35(9), 2326–2333 (2018). [CrossRef]  

10. H. Rao, R. Scarmozzino, and R. M. Osgood, “A bidirectional beam propagation method for multiple dielectric interfaces,” IEEE Photonics Technol. Lett. 11(7), 830–832 (1999). [CrossRef]  

11. Y. Y. Lu and S. H. Wei, “A new iterative bidirectional beam propagation method,” IEEE Photonics Technol. Lett. 14(11), 1533–1535 (2002). [CrossRef]  

12. S. Wu and J. Xiao, “An efficient semivectorial bidirectional beam propagation method for 3-D optical waveguide structures,” J. Lightwave Technol. 34(4), 1313–1321 (2016). [CrossRef]  

13. A. Hofstrand, P. Jakobsen, and J. V. Moloney, “Bidirectional shooting method for extreme nonlinear optics,” Phys. Rev. A 100(5), 053818 (2019). [CrossRef]  

14. J. Zhu and H. Yang, “Efficient marching solver for a class of three-dimensional optical waveguides,” Opt. Express 26(12), 15110–15123 (2018). [CrossRef]  

15. H. Yang and Z. Tang, “3D adjoint-based marching scheme for optical propagation in inhomogeneous waveguides,” J. Opt. Soc. Am. B 37(5), 1298–1306 (2020). [CrossRef]  

16. L. Fisherman, “One-way wave propagation methods in direct and inverse scalar wave propagation modeling,” Radio Sci. 28(5), 865–876 (1993). [CrossRef]  

17. L. Yuan and Y. Y. Lu, “An efficient bidirectional propagation method based on Dirichlet-to-Neumann maps,” IEEE Photonics Technol. Lett. 18(18), 1967–1969 (2006). [CrossRef]  

18. J. Zhu and H. Yang, “Adaptive mode solver for varying refractive index profile’s waveguides with modified spectral element method,” J. Opt. Soc. Am. B 34(3), 538–545 (2017). [CrossRef]  

19. C. C. Huang, C. C. Huang, and J. Yang, “An efficient method for computing optical waveguides with discontinuous refractive index profiles using spectral collocation method with domain decomposition,” J. Lightwave Technol. 21(10), 2284–2296 (2003). [CrossRef]  

20. H. H. Liu and H. C. Chang, “High-resolution analysis of leaky modes in surface plasmon stripe waveguide,” J. Lightwave Technol. 34(11), 2752–2757 (2016). [CrossRef]  

21. J. Zhu and G. Wang, “High-precision computation of optical propagation in inhomogeneous waveguides,” J. Opt. Soc. Am. A 32(9), 1653–1660 (2015). [CrossRef]  

22. L. N. Trefethen, Approximation Theory and Approximation Practice (SIAM, 2013).

23. A. Abdrabou, A. M. Heikal, and S. S. A. Obayya, “Efficient rational Chebyshev pseudo-spectral method with domain decomposition for optical waveguides modal analysis,” Opt. Express 24(10), 10495–10511 (2016). [CrossRef]  

References

  • View by:

  1. A. W. Snyder and J. D. Love, Optical Waveguide Theory (Chapman and Hall, 1983).
  2. R. Scarmozzino, A. Gopinath, R. Pregla, and S. Helfert, “Numerical technique for modeling guided-wave photonic devices,” IEEE J. Sel. Top. Quantum Electron. 6(1), 150–162 (2000).
    [Crossref]
  3. M. Levy, Parabolic Equation Methods for Electromagnetic Wave Propagation (IET, 2009).
  4. S. Obayya, Computational Photonics (Wiley, 2011).
  5. S. Kawai, A. Iguchi, and Y. Tsuji, “Study on high precision and stable finite element beam propagation method based on incomplete third order hybrid edge/nodal element,” J. Lightwave Technol. 36(11), 2278–2285 (2018).
    [Crossref]
  6. X. Ma, S. Zhu, L. Li, J. Li, X. Shao, and A. Galvanauskas, “Finite-difference beam-propagation method for anisotropic waveguides with torsional birefringence,” Opt. Express 26(4), 3995–4003 (2018).
    [Crossref]
  7. Y. Y. Lu, “One - way large range step methods for Helmholtz waveguides,” J. Comput. Phys. 152(1), 231–250 (1999).
    [Crossref]
  8. Y. Y. Lu, “Some techniques for computing wave propagation in optical waveguides,” Commun. Comput. Phys. 1(6), 1056–1075 (2006).
  9. S. Zhao and J. Zhu, “Efficient marching treatment for optical propagation in the waveguide with a discontinuous interface,” J. Opt. Soc. Am. B 35(9), 2326–2333 (2018).
    [Crossref]
  10. H. Rao, R. Scarmozzino, and R. M. Osgood, “A bidirectional beam propagation method for multiple dielectric interfaces,” IEEE Photonics Technol. Lett. 11(7), 830–832 (1999).
    [Crossref]
  11. Y. Y. Lu and S. H. Wei, “A new iterative bidirectional beam propagation method,” IEEE Photonics Technol. Lett. 14(11), 1533–1535 (2002).
    [Crossref]
  12. S. Wu and J. Xiao, “An efficient semivectorial bidirectional beam propagation method for 3-D optical waveguide structures,” J. Lightwave Technol. 34(4), 1313–1321 (2016).
    [Crossref]
  13. A. Hofstrand, P. Jakobsen, and J. V. Moloney, “Bidirectional shooting method for extreme nonlinear optics,” Phys. Rev. A 100(5), 053818 (2019).
    [Crossref]
  14. J. Zhu and H. Yang, “Efficient marching solver for a class of three-dimensional optical waveguides,” Opt. Express 26(12), 15110–15123 (2018).
    [Crossref]
  15. H. Yang and Z. Tang, “3D adjoint-based marching scheme for optical propagation in inhomogeneous waveguides,” J. Opt. Soc. Am. B 37(5), 1298–1306 (2020).
    [Crossref]
  16. L. Fisherman, “One-way wave propagation methods in direct and inverse scalar wave propagation modeling,” Radio Sci. 28(5), 865–876 (1993).
    [Crossref]
  17. L. Yuan and Y. Y. Lu, “An efficient bidirectional propagation method based on Dirichlet-to-Neumann maps,” IEEE Photonics Technol. Lett. 18(18), 1967–1969 (2006).
    [Crossref]
  18. J. Zhu and H. Yang, “Adaptive mode solver for varying refractive index profile’s waveguides with modified spectral element method,” J. Opt. Soc. Am. B 34(3), 538–545 (2017).
    [Crossref]
  19. C. C. Huang, C. C. Huang, and J. Yang, “An efficient method for computing optical waveguides with discontinuous refractive index profiles using spectral collocation method with domain decomposition,” J. Lightwave Technol. 21(10), 2284–2296 (2003).
    [Crossref]
  20. H. H. Liu and H. C. Chang, “High-resolution analysis of leaky modes in surface plasmon stripe waveguide,” J. Lightwave Technol. 34(11), 2752–2757 (2016).
    [Crossref]
  21. J. Zhu and G. Wang, “High-precision computation of optical propagation in inhomogeneous waveguides,” J. Opt. Soc. Am. A 32(9), 1653–1660 (2015).
    [Crossref]
  22. L. N. Trefethen, Approximation Theory and Approximation Practice (SIAM, 2013).
  23. A. Abdrabou, A. M. Heikal, and S. S. A. Obayya, “Efficient rational Chebyshev pseudo-spectral method with domain decomposition for optical waveguides modal analysis,” Opt. Express 24(10), 10495–10511 (2016).
    [Crossref]

2020 (1)

2019 (1)

A. Hofstrand, P. Jakobsen, and J. V. Moloney, “Bidirectional shooting method for extreme nonlinear optics,” Phys. Rev. A 100(5), 053818 (2019).
[Crossref]

2018 (4)

2017 (1)

2016 (3)

2015 (1)

2006 (2)

Y. Y. Lu, “Some techniques for computing wave propagation in optical waveguides,” Commun. Comput. Phys. 1(6), 1056–1075 (2006).

L. Yuan and Y. Y. Lu, “An efficient bidirectional propagation method based on Dirichlet-to-Neumann maps,” IEEE Photonics Technol. Lett. 18(18), 1967–1969 (2006).
[Crossref]

2003 (1)

2002 (1)

Y. Y. Lu and S. H. Wei, “A new iterative bidirectional beam propagation method,” IEEE Photonics Technol. Lett. 14(11), 1533–1535 (2002).
[Crossref]

2000 (1)

R. Scarmozzino, A. Gopinath, R. Pregla, and S. Helfert, “Numerical technique for modeling guided-wave photonic devices,” IEEE J. Sel. Top. Quantum Electron. 6(1), 150–162 (2000).
[Crossref]

1999 (2)

Y. Y. Lu, “One - way large range step methods for Helmholtz waveguides,” J. Comput. Phys. 152(1), 231–250 (1999).
[Crossref]

H. Rao, R. Scarmozzino, and R. M. Osgood, “A bidirectional beam propagation method for multiple dielectric interfaces,” IEEE Photonics Technol. Lett. 11(7), 830–832 (1999).
[Crossref]

1993 (1)

L. Fisherman, “One-way wave propagation methods in direct and inverse scalar wave propagation modeling,” Radio Sci. 28(5), 865–876 (1993).
[Crossref]

Abdrabou, A.

Chang, H. C.

Fisherman, L.

L. Fisherman, “One-way wave propagation methods in direct and inverse scalar wave propagation modeling,” Radio Sci. 28(5), 865–876 (1993).
[Crossref]

Galvanauskas, A.

Gopinath, A.

R. Scarmozzino, A. Gopinath, R. Pregla, and S. Helfert, “Numerical technique for modeling guided-wave photonic devices,” IEEE J. Sel. Top. Quantum Electron. 6(1), 150–162 (2000).
[Crossref]

Heikal, A. M.

Helfert, S.

R. Scarmozzino, A. Gopinath, R. Pregla, and S. Helfert, “Numerical technique for modeling guided-wave photonic devices,” IEEE J. Sel. Top. Quantum Electron. 6(1), 150–162 (2000).
[Crossref]

Hofstrand, A.

A. Hofstrand, P. Jakobsen, and J. V. Moloney, “Bidirectional shooting method for extreme nonlinear optics,” Phys. Rev. A 100(5), 053818 (2019).
[Crossref]

Huang, C. C.

Iguchi, A.

Jakobsen, P.

A. Hofstrand, P. Jakobsen, and J. V. Moloney, “Bidirectional shooting method for extreme nonlinear optics,” Phys. Rev. A 100(5), 053818 (2019).
[Crossref]

Kawai, S.

Levy, M.

M. Levy, Parabolic Equation Methods for Electromagnetic Wave Propagation (IET, 2009).

Li, J.

Li, L.

Liu, H. H.

Love, J. D.

A. W. Snyder and J. D. Love, Optical Waveguide Theory (Chapman and Hall, 1983).

Lu, Y. Y.

Y. Y. Lu, “Some techniques for computing wave propagation in optical waveguides,” Commun. Comput. Phys. 1(6), 1056–1075 (2006).

L. Yuan and Y. Y. Lu, “An efficient bidirectional propagation method based on Dirichlet-to-Neumann maps,” IEEE Photonics Technol. Lett. 18(18), 1967–1969 (2006).
[Crossref]

Y. Y. Lu and S. H. Wei, “A new iterative bidirectional beam propagation method,” IEEE Photonics Technol. Lett. 14(11), 1533–1535 (2002).
[Crossref]

Y. Y. Lu, “One - way large range step methods for Helmholtz waveguides,” J. Comput. Phys. 152(1), 231–250 (1999).
[Crossref]

Ma, X.

Moloney, J. V.

A. Hofstrand, P. Jakobsen, and J. V. Moloney, “Bidirectional shooting method for extreme nonlinear optics,” Phys. Rev. A 100(5), 053818 (2019).
[Crossref]

Obayya, S.

S. Obayya, Computational Photonics (Wiley, 2011).

Obayya, S. S. A.

Osgood, R. M.

H. Rao, R. Scarmozzino, and R. M. Osgood, “A bidirectional beam propagation method for multiple dielectric interfaces,” IEEE Photonics Technol. Lett. 11(7), 830–832 (1999).
[Crossref]

Pregla, R.

R. Scarmozzino, A. Gopinath, R. Pregla, and S. Helfert, “Numerical technique for modeling guided-wave photonic devices,” IEEE J. Sel. Top. Quantum Electron. 6(1), 150–162 (2000).
[Crossref]

Rao, H.

H. Rao, R. Scarmozzino, and R. M. Osgood, “A bidirectional beam propagation method for multiple dielectric interfaces,” IEEE Photonics Technol. Lett. 11(7), 830–832 (1999).
[Crossref]

Scarmozzino, R.

R. Scarmozzino, A. Gopinath, R. Pregla, and S. Helfert, “Numerical technique for modeling guided-wave photonic devices,” IEEE J. Sel. Top. Quantum Electron. 6(1), 150–162 (2000).
[Crossref]

H. Rao, R. Scarmozzino, and R. M. Osgood, “A bidirectional beam propagation method for multiple dielectric interfaces,” IEEE Photonics Technol. Lett. 11(7), 830–832 (1999).
[Crossref]

Shao, X.

Snyder, A. W.

A. W. Snyder and J. D. Love, Optical Waveguide Theory (Chapman and Hall, 1983).

Tang, Z.

Trefethen, L. N.

L. N. Trefethen, Approximation Theory and Approximation Practice (SIAM, 2013).

Tsuji, Y.

Wang, G.

Wei, S. H.

Y. Y. Lu and S. H. Wei, “A new iterative bidirectional beam propagation method,” IEEE Photonics Technol. Lett. 14(11), 1533–1535 (2002).
[Crossref]

Wu, S.

Xiao, J.

Yang, H.

Yang, J.

Yuan, L.

L. Yuan and Y. Y. Lu, “An efficient bidirectional propagation method based on Dirichlet-to-Neumann maps,” IEEE Photonics Technol. Lett. 18(18), 1967–1969 (2006).
[Crossref]

Zhao, S.

Zhu, J.

Zhu, S.

Commun. Comput. Phys. (1)

Y. Y. Lu, “Some techniques for computing wave propagation in optical waveguides,” Commun. Comput. Phys. 1(6), 1056–1075 (2006).

IEEE J. Sel. Top. Quantum Electron. (1)

R. Scarmozzino, A. Gopinath, R. Pregla, and S. Helfert, “Numerical technique for modeling guided-wave photonic devices,” IEEE J. Sel. Top. Quantum Electron. 6(1), 150–162 (2000).
[Crossref]

IEEE Photonics Technol. Lett. (3)

H. Rao, R. Scarmozzino, and R. M. Osgood, “A bidirectional beam propagation method for multiple dielectric interfaces,” IEEE Photonics Technol. Lett. 11(7), 830–832 (1999).
[Crossref]

Y. Y. Lu and S. H. Wei, “A new iterative bidirectional beam propagation method,” IEEE Photonics Technol. Lett. 14(11), 1533–1535 (2002).
[Crossref]

L. Yuan and Y. Y. Lu, “An efficient bidirectional propagation method based on Dirichlet-to-Neumann maps,” IEEE Photonics Technol. Lett. 18(18), 1967–1969 (2006).
[Crossref]

J. Comput. Phys. (1)

Y. Y. Lu, “One - way large range step methods for Helmholtz waveguides,” J. Comput. Phys. 152(1), 231–250 (1999).
[Crossref]

J. Lightwave Technol. (4)

J. Opt. Soc. Am. A (1)

J. Opt. Soc. Am. B (3)

Opt. Express (3)

Phys. Rev. A (1)

A. Hofstrand, P. Jakobsen, and J. V. Moloney, “Bidirectional shooting method for extreme nonlinear optics,” Phys. Rev. A 100(5), 053818 (2019).
[Crossref]

Radio Sci. (1)

L. Fisherman, “One-way wave propagation methods in direct and inverse scalar wave propagation modeling,” Radio Sci. 28(5), 865–876 (1993).
[Crossref]

Other (4)

M. Levy, Parabolic Equation Methods for Electromagnetic Wave Propagation (IET, 2009).

S. Obayya, Computational Photonics (Wiley, 2011).

A. W. Snyder and J. D. Love, Optical Waveguide Theory (Chapman and Hall, 1983).

L. N. Trefethen, Approximation Theory and Approximation Practice (SIAM, 2013).

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (7)

Fig. 1.
Fig. 1. Propagation along the transition region of $N$ transverse inhomogeneous layers.
Fig. 2.
Fig. 2. Buried waveguide.
Fig. 3.
Fig. 3. Numerical fields $EL_R$ for a buried waveguide with $h=1$ $\mu$ m.
Fig. 4.
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.
Fig. 5.
Fig. 5. Surface plasmon stripe waveguide.
Fig. 6.
Fig. 6. Numerical fields $EL_R$ for a surface plasmon stripe waveguide with $h = 60$ $\mu$ m.
Fig. 7.
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.

Tables (4)

Tables Icon

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

Tables Icon

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

Tables Icon

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

Tables Icon

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

Equations (44)

Equations on this page are rendered with MathJax. Learn more.

2 E x 2 + 2 E y 2 + 2 E z 2 + κ 0 2 n 2 ( x , y , z ) E = 0 ,
B ( x , y ) = 2 x 2 + 2 y 2 + κ 0 2 n 2 ( x , y , z ) ,
E ( 0 + , y , z ) = 0 , E ( D 1 , y , z ) = 0 , E ( x , 0 + , z ) = 0 , E ( x , D 2 , z ) = 0.
E ( x , y , 0 ) = E 0 ( x , y ) .
E z = i 2 x 2 + 2 y 2 + κ 0 2 n 2 ( x , y , L ) E , z = L ,
( z + i B ) E R = 0 ,
( z i B ) E L = 0 ,
z E = Q 1 E R + Q 2 E L ,
Q 2 B 2 = 0.
[ d Q 1 d z + Q 1 2 + x 2 + y 2 + κ 0 2 n 2 ( x , y , z ) ] E R = 0 ,
[ d Q 2 d z + Q 2 2 + x 2 + y 2 + κ 0 2 n 2 ( x , y , z ) ] E L = 0.
d Q 1 d z = [ 2 x + 2 y + κ 0 2 n 2 ( x , y , z ) ] Q 1 2 ,
d Q 2 d z = [ 2 x + 2 y + κ 0 2 n 2 ( x , y , z ) ] Q 2 2 .
G 1 E R ( x , y , 0 ) = E R ( x , y , L ) , d G 1 d z = G 1 Q 1 ,
G 2 E L ( x , y , 0 ) = E L ( x , y , L ) , d G 2 d z = G 2 Q 2 .
E R ( x , y , L ) = G 1 ( 0 ) E R ( x , y , 0 ) , E L ( x , y , L ) = G 2 ( 0 ) E L ( x , y , 0 ) .
< 0 = z 0 < z 1 < < z N = L < + .
E x x + E y y + E z z + κ 2 n 2 ( x , y , z i + 1 / 2 ) E = 0 ,
Q 1 = Q 1 2 B 1 / 2 ,
Q 2 = Q 2 2 B 1 / 2 ,
( z + i B 1 / 2 ) E R = 0 ,
( z i B 1 / 2 ) E L = 0.
E L ( x , y , z i ) = R i E R ( x , y , z i ) .
R 1 = [ i B 1 / 2 + Q 2 ( 1 ) ] 1 [ i B 1 / 2 Q 1 ( 1 ) ] .
Q 1 ( 0 ) = i B 1 / 2 ( i B 1 / 2 + Q 2 ( 0 ) ) R 0 .
E R ( x , y , z 1 ) = e i h B 1 / 2 E R ( x , y , z 0 ) , E L ( x , y , z 1 ) = e i h B 1 / 2 E L ( x , y , z 0 ) ,
R 0 = e i h B 1 / 2 R 1 e i h B 1 / 2 .
G 1 ( 0 ) = G 1 ( 1 ) ( I + R 1 ) e i h B 1 / 2 ( I + R 0 ) 1 ,
G 2 ( 0 ) = G 2 ( 1 ) ( I + R 1 ) e i h B 1 / 2 ( I + R 0 ) 1 ,
B ( x , y ) ν ( x , y ) = λ ν ( x , y ) , 0 < x < D 1 , 0 < y < D 2 ,
B ( 1 / 2 ) ( x , y ) = 2 x 2 + 2 y 2 + κ 0 2 n 2 ( x , y , z ) ,
B ( 1 / 2 ) ν ( x , y ) = λ ν ( x , y ) , 0 < x < D 1 , 0 < y < D 2 ,
B ( x , y , z 1 / 2 ) i = 0 n x j = 0 n y B ( x ( ξ i ) , y ( ξ j ) , z 1 / 2 ) T i ( ψ ( y ) ) T j ( ϕ ( x ) ) ,
p B x p | x = x ( ξ i ) , y = y ( ξ j ) , z = z 1 / 2 = i = 0 n x j = 0 n y u ( x ( ξ i ) , y ( ξ j ) , z 1 / 2 ) T j ( p ) ( x ) T i ( y ) .
D 2 = D n x 1 2 = D ~ n x 2 ( 1 : n x 1 , 1 : n x 1 ) ,
D n y 1 2 = D ~ n y 2 ( 1 : n y 1 , 1 : n y 1 ) .
M 1 / 2 = I D 2 + D 2 I + K 1 / 2 ,
ν n , ν ^ m = σ n , m ,
M 1 / 2 ν 1 / 2 = ν 1 / 2 Λ 1 / 2 , M 1 / 2 H ν ^ 1 / 2 = ν ^ 1 / 2 Λ 1 / 2 H .
M 3 / 2 ν 3 / 2 = ν 3 / 2 Λ 3 / 2 , M 3 / 2 H ν ^ 3 / 2 = ν ^ 3 / 2 Λ 3 / 2 H .
ν ^ 1 / 2 ν 1 / 2 = I , ν ^ 3 / 2 ν 3 / 2 = I .
N = ν ^ 1 / 2 , k H ν 3 / 2 , k , if   ν 3 / 2 , k = ν 1 / 2 , k N .
ε := | | E r E | | 2 | | E r | | 2 = ( 1 I k = 1 I ( E r E E r ) 2 ) 1 / 2 ,
E ( x , y , L ) = e i κ 0 2 8 π 2 L sin ( 2 π x ) sin ( 2 π y ) .

Metrics