Efficient marching solver for a class of three-dimensional optical waveguides

Jianxin Zhu; Hao Yang

doi:10.1364/OE.26.015110

1. Introduction

Numerical simulation of wave propagation is important for the optimization design of optical waveguide devices [1–3]. Applying the Fourier transform with respect to time to the wave equation, we can obtain an other equation in the frequency domain, which is called as the well-known Helmholtz equation. If the medium is homogeneous, then the refractive-index or wavenumber is constant [4, 5]. However, in many cases the medium is non-homogeneous, in which case the refractive-index or the wavenumber usually depends on the position, namely the refractive-index or wavenumber profile is a varying function.

Most Helmholtz equations with the varying refractive-index or wavenumber solvers were initially developed for the two-dimensional (2D) case, such as the finite difference, finite element and spectral methods, which are used to lead to the systems of equations with extremely large sizes of unknowns and were not very practical. There is an understandable desire to incorporate the marching (one-way) method into the finding solution algorithm for the inherently elliptic, frequency domain, a scalar Helmholtz formulation [6–8]. For medium that has gradual variations in the main propagation direction, instead of directly discretizing the equation on a compact stencil, Lu in [9,10] proposed an useful one-way operator marching method (OMM) for the 2D Helmholtz waveguides based on the Dirichlet-to-Neumann (DtN) map. Zhu in [11] generalized the OMM to solve the problem of the optical propagation in micro-waveguides with loss that is described by a 2D Helmholtz equation with complex refractive index or wavenumber. In [12], a 2D modified OMM was provided to accurately compute the optical propagation in the inhomogeneous waveguide terminated by a perfectly matched layer. In [13], the wave propagation in a curved interface waveguide was computed by a modified OMM, firstly a local orthogonal coordinate transform and an equation transform were applied, then a perfectly matched layer was used to terminate the open waveguide to a bounded one.

In this paper, we generalize the OMM algorithm to deal with the 3D case and solve the Helmholtz equation in a 3D rectangular box. The reflection operator is firstly determined, subsequently the marching method is constructed to recover the refractive-index filed from the DtN map, and finally the local base transformation is implemented by the marching method. The efficiency of the algorithm is especially vital for 3D problems which usually require heavy computations. The method that is presented here enjoys the properties of the 2D algorithm: large range step (i.e. using large marching step to achieve the prescribed accuracy) and comparatively small number of operations for each step. If the Helmholtz equation (1) is directly solved by the 3D OMM, the required memory space is only O(N⁴), with the operator T beign approximated by an N² × N² matrix, where N is the number of grid points in each one of the two perpendicular directions. If m marching steps are used over the range direction, the total operations required for the OMM is O(N⁵m). Moreover, during the local base transformation in OMM, an approach of local subtraction of the eigenmodes is implemented, which means it is only necessary to focus on the first n eigenmodes. In truth, we generally have N ≈ 10n [14]. Therefore, each integration step for T requires only O(n⁵) operations.

The structure of the paper is as follows. Sections 2,3 and 4 explain the formulation of the 3D OMM in detail. Section 5 presents the numerical simulation results for two typical boundary conditions. Whenever possible, comparisons have been made with exact solutions and larger marching steps h to assess the accuracy and the fast computation of the proposed method. Finally, Section 6 briefly explains the main conclusions of the work.

2. Basic equations

The 3D scalar Helmholtz equation, which is the basis of the OMM, is expressed by

\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} + \frac{\partial^{2} u}{\partial z^{2}} + κ^{2} (x, y, z) u = 0,

where κ(x, y, z) is called the wavenumber and

κ^{2} (x, y, z) = κ_{0}^{2} n^{2} (x, y, z);

here κ₀ is a free space wavenumber, and n(x, y, z) is a refractive-index function. If the medium is homogeneous, then κ is constant. Here κ is the varying wavenumber, which is x-invariant for x < 0 (κ(x, y, z) = κ₀(y, z)) and x > L(κ(x, y, z) = κ_∞(y, z)). Suppose

\frac{1}{κ} ≪ 1 ≪ L

, the restriction on κ(x, y, z) shows that the wave propagates in the medium over a long distance. To this end, we introduce the reference 3D optical waveguide with the reference-index function n(x, y, z) indicated in Fig. 1, where there are the outgoing waves at x = L.

Fig. 1 3D optical waveguide with outgoing waves at x = L.

Download Full Size | PDF

In the above expressions, x is represented as the principal propagating variable or range coordinate. Therefore, the transversal discretization associated with the perpendicular, or cross-range, coordinates {x_⊥} must be performed in the y and z directions. Corresponding to the pressure release surfaces at y = 0, z = 0 and total reflection rigid bottoms at y = D₁, z = D₂, we use the boundary conditions either Dirichlet boundary conditions,

\begin{array}{l} u (x, 0^{+}, z) = 0, u (x, D_{1}^{-}, z) = 0, \\ u (x, y, 0^{+}) = 0, u (x, y, D_{2}^{-}) = 0, \end{array}

or Neumann boundary conditions are imposed [4,15,16],

\begin{array}{l} \frac{\partial u}{\partial x} (x, 0^{+}, z) = 0, \frac{\partial u}{\partial x} (x, D_{1}^{-}, z) = 0, \\ \frac{\partial u}{\partial x} (x, y, 0^{+}) = 0, \frac{\partial u}{\partial x} (x, y, D_{2}^{-}) = 0 . \end{array}

Mixed Dirichlet - Neumann boundary conditions can be treated by a similar approach where each interface can be associated with one of the two types of boundary conditions. First we will introduce the algorithm for the Dirichlet problem, and later we will describe the changes necessary to accommodate the Neumann problem.

In the x direction, we also restrict ourself to a finite interval [0, L]. At x = 0, we temporarily assume that a starting field is specified

u (0, y, z) = u_{0} (y, z) .

A radiation condition is imposed at infinity such that u · e^−iωt represents the out-going waves as x → ∞. Under the assumption that the medium is x independent for x ≥ L, the radiation condition can be written down at x = L and the Helmholtz equation needs only to be solved in the cuboid domain

Ω_{e} = {(x, y, z) | 0 < x < L, 0 < y < D_{1}, 0 < z < D_{2}} .

Let κ(x, y, z) = κ_∞(y, z) for x ≥ L. We first consider the eigenvalue problem of the transverse operator:

[\nabla_{⊥}^{2} + κ_{\infty}^{2} (x_{⊥})] ϕ_{j} = λ_{j} ϕ_{j},

where

\nabla_{⊥}^{2} = \frac{\partial^{2}}{\partial y^{2}} + \frac{\partial^{2}}{\partial z^{2}}

and 0 < y < D₁, 0 < z < D₂.

The boundary conditions of ϕ_j are similarly as above (3) or (4). Since the medium is x independent for x ≥ L, the solution u satisfying the radiation condition can be written down based on u(L, y, z) [8]. In fact,

u (x, y, z) = \sum_{j = 1}^{\infty} α_{j} e^{i \sqrt{λ_{j} (x - L)}} ϕ_{j} (y, z) for x \geq L,

where {α_j} are the eigenfunction expansion coefficients of u(L, y, z). Namely,

u (L, y, z) = \sum_{j = 1}^{\infty} α_{j} ϕ_{j} (y, z), α_{j} = \frac{\int_{0}^{D_{1}} \int_{0}^{D_{2}} u (L, y, z) ϕ_{j} (y, z) d y d z}{\int_{0}^{D_{1}} \int_{0}^{D_{2}} ϕ_{j}^{2} (y, z) d y d z} .

The radiation condition at infinity is satisfied, since only out-going waves and exponentially decaying modes are retained in the solution. The positive eigenvalues correspond to the propagating modes. The negative eigenvalues correspond to the evanescent modes whose waves decay exponentially for increasing x.

If we define the square root operator $\sqrt{\nabla_{⊥}^{2} + κ_{\infty}^{2} (x_{⊥})}$ to be a linear operator satisfying

\sqrt{\nabla_{⊥}^{2} + κ_{\infty}^{2} (x_{⊥})} ϕ_{j} = \sqrt{λ_{j}} ϕ_{j}

(with

\sqrt{λ_{j}} = i \sqrt{- λ_{j}}

if λ_j is negative), then the solution (8) is the exact solution of the one-way differential operation equation [9,10,17]:

u_{x} = i \sqrt{\nabla_{⊥}^{2} + κ_{\infty}^{2} (x_{⊥}) u} for x \geq L .

The solutions of the above equation satisfy the radiation condition automatically. Eq. (11) at x = L serves as the radiation condition for the Helmholtz equation now reduced to the cuboid domain Ω_e = {(x, y, z)|0 < y < D₁, 0 < z < D₂, 0 < x < L}.

3. The operator marching

A Dirichlet-to-Neumann map T(x) [17, 18] is utilized to map any solution of the Helmholtz equation to its x derivative, that is, u_x = T(x)u. The operator T acts as the function of y, z defined on the domain (0, D₁) × (0, D₂) and it depends on the parameter x. Substituting u_x = T(x)u into the Helmholtz equation, we get

[\frac{d T}{d x} + T^{2} + \partial_{y}^{2} + \partial_{z}^{2} + κ^{2} (x, y, z)] u = 0 .

Since the above formula is true for any solution of the Helmholtz equation (satisfying the boundary conditions at (3) or (4), and the radiation condition at x = L), we obtain

\frac{d T}{d x} = - [\partial_{y}^{2} + \partial_{z}^{2} + κ^{2} (x, y, z)] - T^{2} .

The above equation should be solved for decreasing x with an exact radiation condition at +∞. For the stability, a fundamental solution operator G has been introduced, which satisfies G(x)u(x, y, z) = u(L, y, z) and

\frac{d G}{d x} = - GT .

The first step is now to solve T and G together from x = L to x = 0. The initial condition at x = L is G(L) = I, where I is the identity operator. When G(0) is calculated, the solution at x = L can be generated by a simple multiplication with the starting field u(0, y, z), namely,

u (L, y, z) = G (0) u (0, y, z) .

The range is divided into uniform piece-wise x-independent segments. Let the discrete points {x_i} satisfy

- \infty < 0 = x_{0} < x_{1} < \dots < x_{N} = L < + \infty .

Each interval (x_i, x_i+1) corresponds to a range-independent segment and h = L/N. The approximated Helmholtz equation on (x_i, x_i+1) is

u_{x x} + u_{y y} + u_{z z} + κ^{2} (x_{i + 1 / 2}, y, z) u = 0 .

where x_i+1/2 = x_i + h/2 = (x_i + x_i+1)/2. We take i = 0 for example to describe the marching from x₀ to x₁.

On the interval (x₀, x₁), the Eq. (13) is approximated by

T^{'} = - T^{2} - [\partial_{y}^{2} + \partial_{z}^{2} + κ^{2} (x_{1 / 2}, y, z)] .

For a given initial condition T₁ ≈ T(x₁), the exact solution of T₀ is used to obtain T₀ ≈ T(x₀). For this purpose, we explore the relationship with the Helmholtz equation and consider the associated equation

u_{x x} + u_{y y} + u_{z z} + κ^{2} (x_{1 / 2}, y, z) u = 0

on (x₀, x₁). Eq. (19) can be factored as [19,20]:

(\partial_{x} + i β_{1 / 2}) (\partial_{x} - i β_{1 / 2}) u = 0,

where

β_{1 / 2} = \sqrt{\partial_{y}^{2} + \partial_{z}^{2} + κ^{2} (x_{1 / 2}, y, z)}

. Then the wavefield can be decomposed into the right and the left going waves, denoted by

u_{R}^{+}

and

u_{L}^{-}

respectively, satisfying

(\partial_{x} + i β_{1 / 2}) u_{R}^{+} = 0,

and

(\partial_{x} - i β_{1 / 2}) u_{L}^{-} = 0 .

Then we can denote $ℜ_{R}^{x_{i}}$ , $ℜ_{L}^{x_{i}}$ and $ℑ_{R}^{x_{i}}$ , $ℑ_{L}^{x_{i}}$ are the reflection and transmission operators for the transition region at x_i [6], which are depicted as in Fig. 2. And we can obtain

u_{R}^{+} (x_{1}) = ℑ_{R}^{x_{1}} u_{R}^{-} (x_{0}) + ℜ_{L}^{x_{1}} u_{L}^{+} (x_{1}),

and

u_{L}^{-} (x_{0}) = ℜ_{R}^{x_{1}} u_{R}^{-} (x_{0}) + ℑ_{L}^{x_{1}} u_{L}^{+} (x_{1}) .

Fig. 2 The reflection and transmission operators for N transverse transition regions.

Download Full Size | PDF

It is easy to obtain from the above relation [7,10,17] that

ℜ_{L}^{x_{1}} = {[i β_{1 / 2} + T_{1}]}^{- 1} [i β_{1 / 2} - T_{1}],

and obtain

T_{0} = i β_{1 / 2} [I - ℜ_{L}^{x_{0}}] {[I + ℜ_{L}^{x_{0}}]}^{- 1} .

Consider the marching equations,

u^{+} (x_{1}, y, z) = e^{i h β_{1 / 2}} u^{+} (x_{0}, y, z), u^{-} (x_{1}, y, z) = e^{- i h β_{1 / 2}} u^{-} (x_{0}, y, z),

we can have

ℜ_{L}^{x_{0}} = e^{i h β_{1 / 2}} ℜ_{L}^{x_{1}} e^{i h β_{1 / 2}} .

On the interval (x₀, x₁), supposing that the Helmholtz equation is approximated by (17), we can get

u (x_{1}, y, z) = (I + ℜ_{L}^{x_{1}}) e^{i h β_{1 / 2}} {(I + ℜ_{L}^{x_{0}})}^{- 1} u (x_{0}, y, z) .

Therefore, we have

G_{0} = G_{1} (I + ℜ_{L}^{x_{1}}) e^{i h β_{1 / 2}} {(I + ℜ_{L}^{x_{0}})}^{- 1} .

When the DtN map Q and the fundamental solution operator G are solved from x = L to x = 0, the formulas (25), (28), (29), and (30) are used in the step from x_i+1 to x_i, where 1 ≤ i ≤ N. Moreover, the impedance transport across range-independent segments are the same, i.e., the impedance of the free-sapce [21].

4. Local base transformation

In this section, we consider the numerical implementation of the above marching formulas (21), (23), (24) and (26). A direct approach is to approximate the operators by matrices.

The transverse operator of Helmholtz equation is

T (x) = \frac{\partial^{2}}{\partial y^{2}} + \frac{\partial^{2}}{\partial z^{2}} + κ^{2} (x, y, z),

where x is fixed. The eigenvalue problem of T(x) is important. It is defined as

T (x) ϕ (y, z) = λ ϕ (y, z), 0 < y < D_{1}, 0 < z < D_{2},

corresponding to the boundary conditions (3) or (4).

The square root operator,

T^{1 / 2} (x) = \sqrt{\frac{\partial^{2}}{\partial y^{2}} + \frac{\partial^{2}}{\partial z^{2}} + κ^{2} (x, y, z)},

is definded to satisfy

T^{1 / 2} (x) ϕ (y, z) = \sqrt{λ} ϕ (y, z), 0 < y < D_{1}, 0 < z < D_{2},

where ϕ(y, z) and λ are the same as above. The branch of the square root is chosen as (−π/2, π/2].

We also take the segment (x₀, x₁) for example. The marching formulas indicate the relation of operators from x₁ to x₀. For Eqs. (2) and (3), the y or z is discretized by z_j = y_j = jδ, δ = (x₁ − x₀)/n, j = 1, 2, · · · , n. The second order approximation to $\partial_{y}^{2}$ or $\partial_{z}^{2}$ is

T_{n} = \frac{1}{δ^{2}} (\begin{matrix} - 2 & 1 \\ 1 & ⋱ & ⋱ \\ ⋱ & - 2 & 1 \\ 1 & - 2 \end{matrix}),

For Eqs. (2) and (4), the y or z axis is discretized by y_j = z_j = jδ, δ = (x₁ − x₀)/(n + 1/2), j = 1, 2, · · · , n. The second order approximation to $\partial_{y}^{2}$ or $\partial_{z}^{2}$ is [10]

T_{n} = \frac{1}{δ^{2}} (\begin{matrix} - 1 & 1 \\ 1 & - 2 & 1 \\ ⋱ & ⋱ & ⋱ \\ 1 & - 2 & 1 \\ 1 & - 1 \end{matrix}) .

So the operator

T (x_{1 / 2}) = \partial_{y}^{2} + \partial_{z}^{2} + κ^{2} (x_{1 / 2}, y, z)

can be discretized as

\bar{M} (x_{1 / 2}) = I \otimes T_{n} + T_{n} \otimes I + k^{2} (x_{1 / 2}, y, z) I \otimes I,

where ⊗ is the kronecker product [22,23], T_n is the n × n matrix in Eq. (31) (or (32)) and I is the n × n identity matrix.

To implement formulas (21), (23), (24) and (26), it is necessary to find the eigendecomposition for the matrix that approximates M̄(x_1/2).

Let T_n = ZΛZ be the eigendecompositon of T_n, with Z = [z₁, · · · , z_n] the orthogonal matrix whose columns are eigenvectors, and Λ = diag(λ₁, · · · , λ_n) is the eigenvalue matrix. It can be easily proved that [23] the eigendecomposition of M̄(x_1/2) is

\bar{M} (x_{1 / 2}) = V A_{1 / 2} V^{T} = (Z \otimes Z) \cdot (I \otimes Λ + Λ \otimes I + k^{2} (x_{1 / 2}, y, z) I \otimes I) \cdot {(Z \otimes Z)}^{T},

where A_1/2 = I ⊗ Λ + Λ ⊗ I + k²(x_1/2, y, z)I ⊗ I is a diagonal matrix. V = Z ⊗ Z is an orthogonal matrix whose (i · n + j) th column, the corresponding eigenvector, is z_i ⊗ z_j, i, j = 1, · · · , n. The resulting matrix eigenvalue and eigenfunctions problems could be both solved by the MATLAB eigenvalue solver.

Denote the discrete matrices of M̄(x_1/2) and M̄(x_3/2) as M̄₀ and M̄₁, respectively. Their truncated local eigenvector decompositions are

{\bar{M}}_{0} = V_{0} A_{0} V_{0}^{T}, {\bar{M}}_{1} = V_{1} A_{1} V_{1}^{T},

where V₀, V₁ are the orthogonal matrices of the eigenvectors, and A₀, A₁ are the diagonal matrices of the eigenvalues. As the so-called coupled mode method and operator marching method for 2D problem, here we only consider the first m = ratio · n ² largest eigenvalues (Usually, it is sufficient to choose m slightly larger than the number of propagating modes, including all the positive eigenvalues and few negative eigenvalues of the operator

\partial_{y}^{2} + \partial_{z}^{2} + κ^{2}

[10].) and eigenfunctions of M̄₀ and M̄₁. So we can replace (36) by

{\bar{M}}_{0} V_{0, m} = V_{0, m} A_{0, m}, {\bar{M}}_{1} V_{1, m} = V_{1, m} A_{1, m},

where V_0,m, V_1,m are the n² × m eigenvextor matrices and A_0,m, A_1,m are the m × m diagonal eigenvalue matrices.

Here the local base transform can be done by

N = V_{0, m}^{T} V_{1, m} .

Then the marching formulas from x₁ to x₀ (19), (20), (21) and (22) can be transformed to [10,11]

\begin{array}{l} O = {NO}_{1} N^{- 1}, \\ P_{1} = {(i \sqrt{A_{0, m}} + O)}^{- 1} (i \sqrt{A_{0, m}} - O), \\ P_{0} = e^{i h_{1} \sqrt{A_{0, m}}} P_{1} e^{i h_{1} \sqrt{A_{0, m}}}, \\ S_{0} = N S_{1} N^{- 1} (I + P_{1}) e^{i h_{1} \sqrt{A_{0, m}}} {(I + P_{0})}^{- 1}, \end{array}

where O₁ and S₁ are known at x₁, and h₁ = x₁ − x₀. The wave field at x = L is solved by

u L ≜ u (L, y, z) = V_{0, m} S_{0} V_{0, m}^{T} u (0, y, z) .

5. Numerical examples and discussion

In this section, several examples are presented to show the performance of our method. The problems (2) with both homogeneous and inhomogeneous media are considered, which are denoted as Cases 1 and 2, respectively. To verify the validity of our method, we firstly consider the homogeneous medium, for the reason that exact solutions can be obtained.

For the homogeneous medium, let

n^{2} (x, y, z) = 1, κ^{2} (x, y, z) = κ_{0}^{2};

and for the inhomogeneous media, let

n^{2} (x, y, z) = [1 + 0.05 e^{- 20 {(x / L - 0.5)}^{2}} {sin}^{2} (π y) {sin}^{2} (π z)], κ^{2} (x, y, z) = κ_{0}^{2} n^{2} (x, y, z);

where the wavelength is assumed to be λ₀ = 0.6328μm [11], k₀ = 2π/λ₀, L = 10, n = 400 and ratio = 0.4 for both cases. The fields u(L, y, z) are computed by different step lengths to show the merit of large steps for the proposed method.

5.1. Problems (2) with boundary conditions (3), (5) and (8)

Firstly, we consider the homogeneous medium. The starting field at x = 0 is u₀(y, z) = sin(2πy) sin(2πz), which is the third propagating mode. Then Case 1 has an exact solution

u (x, y, z) = e^{i \sqrt{κ_{0}^{2} - {(2 π)}^{2} - {(2 π)}^{2} x}} sin (2 π y) sin (2 π z),

which should be used for reference.

In Figs. 3 and 4, the exact field (47) and the numerical field u(L, y, z) with the step length h = 1 are plotted separately. One can see little difference between the numerical fields and the exact solutions except the oscillations around the y, z coordinates due to the chosen local base. It indicates that reasonably good solutions have been obtained with h = 1. Since κ is constant here, the relative errors almostly don’t change with h, which can be indicated in Fig. 5 for Case 1. In Fig. 5, $E (u) = \frac{{‖ u_{r} - u ‖}_{2}}{{‖ u_{r} ‖}_{2}}$ , where u_r is the exact field and u is the numerical field with the same step size h.

Fig. 3 Exact solution of fields uL for Case 1: homogeneous medium.

Download Full Size | PDF

Fig. 4 Numerical fields uL for Case 1: homogeneous medium.

Download Full Size | PDF

Fig. 5 Relative error curves of u(L, y, z) for the two cases in homogeneous media.

Download Full Size | PDF

Secondly, the inhomogeneous media is considered. Its starting field at x = 0 is

u_{0} (y, z) = \sum_{j = 1}^{5} sin (m_{j} y_{0}) sin (m_{j} z_{0}) sin (m_{j} y) sin (m_{j} z) / \sqrt{k_{0}^{2} - m_{j}^{2} - m_{j}^{2}},

for m_j = (j − 1/2)π, y₀ = z₀ = 0.65.

This problem doesn’t have an exact solution. In this case, the numerical fields with step length h = 1/128 are calculated for reference, which are plotted in Fig. 6. After that, we calculate the solution with much larger range steps and compare them with the more accurate reference solution [10]. The relative errors are shown in Tab. 1, where L₂-err (the relative L₂ error) is $E (u) = \frac{{‖ u_{r} - u ‖}_{2}}{{‖ u_{r} ‖}_{2}}$ , where u_r is the reference field with the step length 1/128 and u is the numerical field with the step size h. From Tab. 1, it indicates a reasonably accurate solution is already got with h = 1, which is also plotted in Fig. 7.

Fig. 6 Numerical fields uL of h = 1/128 for Case 2: inhomogeneous media.

Download Full Size | PDF

Table 1. Relative errors E(u) for different range step sizes.

View Table | View all tables in this article

Fig. 7 Numerical fields uL of h = 1 for Case 2: inhomogeneous media.

Download Full Size | PDF

5.2. Problems (2) with boundary conditions (4), (5) and (8)

As above, we firstly consider the homogeneous medium. The starting field at x = 0 is u₀(y, z) = cos(2πy) cos(2πz) corresponding to the third propagating mode. Then the problem has an exact solution

u (x, y, z) = e^{i \sqrt{κ_{0}^{2} - {(2 π)}^{2} - {(2 π)}^{2} x}} cos (2 π y) cos (2 π z) .

In Figs. 8 and 9, the exact field (49) and the numerical field u(L, y, z) with the step size h = 1 are plotted separately. The numerical fields in Fig. 9 coincide well with the exact solutions in Fig. 8 except a little difference along the x coordinate, which shows good solutions could be got with h = 1. The relative errors also don’t change with h, for the reason that κ is constant here, which can be shown in Fig. 5 for Case 2.

Fig. 8 Exact solution of fields uL for Case 2: homogeneous medium.

Download Full Size | PDF

Fig. 9 Numerical fields uL for Case 2: homogeneous medium.

Download Full Size | PDF

In the inhomogeneous media, it is reasonable to assume the starting field at x = 0 is Eq. (48) too. The numerical field uL are calculated with h = 1/128 as a reference that are plotted in Fig. 10. Meanwhile, the numerical fields uL with h = 1 are also plotted in Fig. 11. The relative numerical errors as above are presented in Tab. 2, which also indicates reasonably accurate solutions could be obtained with h = 1.

Fig. 10 Numerical fields uL of h = 1/128 for Case 2: inhomogeneous media.

Download Full Size | PDF

Table 2. Relative errors E(u) for different range step sizes.

View Table | View all tables in this article

Fig. 11 Numerical fields uL of h = 1 for Case 2: inhomogeneous media.

Download Full Size | PDF

6. Conclusion

An efficient marching solver for 3D Helmholtz equation with the varying wavenumber or refractive-index is given in this study. Compared to other standard numerical methods, such as the finite difference, finite element and spectral methods, our algorithm requires the less order of operations, that is, O(n⁶) for each range step, where n is assumed as the number of grid points in each perpendicular direction [10]. Furthermore, if the truncated technique is used, the amount of the operations may be reduced to O(m³) for each range step where m = ratio · n² (0 < ratio < 1). Numerical experiments has proved a second order accuracy of the piecewise approximation. It can be further improved if the fourth order method [24] is used. If the transport of polarization and coupling in electric field are taken into account, the fast marching method for the eikonal [25] could be combined to construct a vectorial operator marching method.

Funding

Ministry of Science and Technology of the People’s Republic of China (2018YFB1107601); National Natural Science Foundation of China (11371319).

References and links

1. O. K. Ersoy, Diffraction, Fourier Optics and Imaging (Wiley, 2007). [CrossRef]

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

3. K. Gallo and G. Assanto, “All-optical diode based on second-harmonic generation in an asymmetric waveguide,” J. Opt. Soc. Am. B 16(2), 267–269 (1999). [CrossRef]

4. E. Braverman, M. Israeli, and A. Averbuch, “A fast spectral solver for a 3D Helmholtz equation,” SIAM J. Sci. Comput. 20(6), 2237–2260 (1998). [CrossRef]

5. E. Turkel, D. Gordon, R. Gordon, and S. Tsynkov, “Compact 2D and 3D sixth order schemes for the Helmholtz equation with variabel wave number,” J. Comput. Phys. 232(1), 272–287 (2013). [CrossRef]

6. J. McCoy and L. N. Frazer, “Propagation modelling based on wavefield factorization and invariant imbedding,” Geophys. J. R. Astron. Soc. 86(3), 703–717(1986). [CrossRef]

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

8. L. Fishman and J. J. McCoy, “A new class of propagation models based on a factorization of the Helmholtz equaiton,” Geophys. J. Roy. Astron. Soc. 50(2), 439–461 (1985).

9. Y. Y. Lu, “Exact one-way methods for acoustic waveguides,” Math. Comput. Simulat. 50 (5–6), 377–391 (1999). [CrossRef]

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

11. J. Zhu and R. Song, “Fast and stable computatoin of optical propagation in micro-waveguides with loss,” Microelectron. Reliab. 49(12), 1529–1536 (2009). [CrossRef]

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

13. J. Zhu and G. Wang, “Fast computation of wave propagation in the open acoustical waveguide with a curved interface,” Wave Motion 57, 171–181 (2015). [CrossRef]

14. Y. Y. Lu and J. R. McLaughlin, “The Riccati method for the Helmholtz equation,” J. Acoust. Soc. Am. 100(3), 1432–1446 (1996). [CrossRef]

15. A. Ortega-Moñux, I. Molina-fernández, and J. G. Wangüemert-pérez, “3D-scalar fourier eigenvector expansion method (Fourier-EEM) for analyzing optical waveguide discontinuities,” Opt. Quant. Electron. 37(1), 213–228 (2005). [CrossRef]

16. G. Ciraolo, F. Gargano, and V. Sciacca, “A spectral approach to a constrained optimization problem for the Helmholtz equation in unbounded domains,” Comput. Appl. Math. 34(3), 1035–1055 (2015). [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. D. Givoli, I. Patlashenko, and J. B. Keller, “Discrete Dirichlet-to-Neumann maps for unbounded domains,” Comput. Methods Appl. Mech. Engrg. 164(1–2), 173–185 (1998). [CrossRef]

19. D. A. Angus, “The one-way wave equation: A full-waveform tool for modeling seismic body wave phenomena,” Surv. Geophys. 35(2), 359–393 (2014). [CrossRef]

20. D. Lee and W. L. Siegmann, “A mathematical model for the 3-dimensional ocean sound propagation,” Math. Model. 7(2), 143–162 (1986). [CrossRef]

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

22. S. H. Lui, “Kron’s method for symmetric eigenvalue problems,” J. Comput. Appl. Math. 98(1), 35–48 (1998). [CrossRef]

23. J. W. Demmel, Applied Numerical Linear Algebra (SIAM, 1997). [CrossRef]

24. Y. Y. Lu, “A Fourth order derivative-free operator marching method for the Helmholtz equation in waveguides,” J. Comput. Math. 25(6), 719–729 (2007).

25. A. Capozzoli, C. Curcio, A. Liseno, and S. Savarese, “Two-dimensional fast marching for geometrical optics,” Opt. Express , 22(22), 26680–26695 (2014). [CrossRef] [PubMed]

Efficient marching solver for a class of three-dimensional optical waveguides

Abstract

1. Introduction

2. Basic equations

3. The operator marching

4. Local base transformation

5. Numerical examples and discussion

5.1. Problems (2) with boundary conditions (3), (5) and (8)

5.2. Problems (2) with boundary conditions (4), (5) and (8)

6. Conclusion

Funding

References and links

Cited By

Figures (11)

Tables (2)

Equations (49)

Optics Express