## Abstract

An efficient marching method of the wave propagation computation in the three-dimensional (3D) rectangular box is developed by solving the scalar Helmholtz equation with the weakly range-dependent refractive-index or wavenumber profile. The method is based on the application of the one-way formulation accompanied by a truncated technique that allows for the reduction of the number of eigenmodes and for quickly achieving reliable numerical results. The merits of the high accuracy and fast computation with a larger range step are shown to solve the equation. This treatment can be also applied to other complicated 3D waveguides.

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

## 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*^{4}), with the operator *T* beign approximated by an *N*^{2} × *N*^{2} 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*^{5}*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* ≈ 10*n* [14]. Therefore, each integration step for *T* requires only *O*(*n*^{5}) 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

*κ*(

*x*,

*y*,

*z*) is called the wavenumber and here

*κ*

_{0}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*) =

*κ*

_{0}(

*y*,

*z*)) and

*x*>

*L*(

*κ*(

*x*,

*y*,

*z*) =

*κ*

_{∞}(

*y*,

*z*)). Suppose $\frac{1}{\kappa}\ll 1\ll 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*.

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*_{1}, *z* = *D*_{2}, we use the boundary conditions either Dirichlet boundary conditions,

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*· 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

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

*y*<

*D*

_{1}, 0 <

*z*<

*D*

_{2}.

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,

*α*} are the eigenfunction expansion coefficients of

_{j}*u*(

*L*,

*y*,

*z*). Namely,

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}_{\perp}^{2}+{\kappa}_{\infty}^{2}({x}_{\perp})}$ to be a linear operator satisfying

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

_{j}*x*=

*L*serves as the radiation condition for the Helmholtz equation now reduced to the cuboid domain Ω

*= {(*

_{e}*x*,

*y*,

*z*)|0 <

*y*<

*D*

_{1}, 0 <

*z*<

*D*

_{2}, 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*

_{1}) × (0,

*D*

_{2}) and it depends on the parameter

*x*. Substituting

*u*=

_{x}*T*(

*x*)

*u*into the Helmholtz equation, we get

*x*=

*L*), we obtain

*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

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,

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

*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 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*

_{0}to

*x*

_{1}.

On the interval (*x*_{0}, *x*_{1}), the Eq. (13) is approximated by

*T*

_{1}≈

*T*(

*x*

_{1}), the exact solution of

*T*

_{0}is used to obtain

*T*

_{0}≈

*T*(

*x*

_{0}). For this purpose, we explore the relationship with the Helmholtz equation and consider the associated equation on (

*x*

_{0},

*x*

_{1}). Eq. (19) can be factored as [19,20]:

Then we can denote ${\Re}_{R}^{{x}_{i}}$, ${\Re}_{L}^{{x}_{i}}$ and ${\Im}_{R}^{{x}_{i}}$, ${\Im}_{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

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

Consider the marching equations,

On the interval (*x*_{0}, *x*_{1}), supposing that the Helmholtz equation is approximated by (17), we can get

Therefore, we have

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

*x*is fixed. The eigenvalue problem of

*T*(

*x*) is important. It is defined as

The square root operator,

*ϕ*(

*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*_{0}, *x*_{1}) for example. The marching formulas indicate the relation of operators from *x*_{1} to *x*_{0}. For Eqs. (2) and (3), the *y* or *z* is discretized by *z _{j}* =

*y*=

_{j}*jδ*,

*δ*= (

*x*

_{1}−

*x*

_{0})/

*n*,

*j*= 1, 2, · · · ,

*n.*The second order approximation to ${\partial}_{y}^{2}$ or ${\partial}_{z}^{2}$ is

For Eqs. (2) and (4), the *y* or *z* axis is discretized by *y _{j}* =

*z*=

_{j}*jδ*,

*δ*= (

*x*

_{1}−

*x*

_{0})/(

*n*+ 1/2),

*j*= 1, 2, · · · ,

*n.*The second order approximation to ${\partial}_{y}^{2}$ or ${\partial}_{z}^{2}$ is [10]

So the operator

can be discretized as where ⊗ is the kronecker product [22,23],*T*is the

_{n}*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*, with

_{n}*Z*= [

*z*

_{1}, · · · ,

*z*] the orthogonal matrix whose columns are eigenvectors, and Λ = diag(

_{n}*λ*

_{1}, · · · ,

*λ*) is the eigenvalue matrix. It can be easily proved that [23] the eigendecomposition of

_{n}*M̄*(

*x*

_{1/2}) is

*A*

_{1/2}=

*I*⊗ Λ + Λ ⊗

*I*+

*k*

^{2}(

*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̄*_{0} and *M̄*_{1}, respectively. Their truncated local eigenvector decompositions are

*V*

_{0},

*V*

_{1}are the orthogonal matrices of the eigenvectors, and

*A*

_{0},

*A*

_{1}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*

^{2}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}+{\kappa}^{2}$[10].) and eigenfunctions of

*M̄*

_{0}and

*M̄*

_{1}. So we can replace (36) by

*V*

_{0,m},

*V*

_{1,m}are the

*n*

^{2}×

*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

Then the marching formulas from *x*_{1} to *x*_{0} (19), (20), (21) and (22) can be transformed to [10,11]

*O*

_{1}and

*S*

_{1}are known at

*x*

_{1}, and

*h*

_{1}=

*x*

_{1}−

*x*

_{0}. The wave field at

*x*=

*L*is solved by

## 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

*λ*

_{0}= 0.6328

*μ*m [11],

*k*

_{0}= 2

*π*/

*λ*

_{0},

*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*_{0}(*y*, *z*) = sin(2*πy*) sin(2*πz*), which is the third propagating mode. Then Case 1 has an exact solution

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{{\Vert {u}_{r}-u\Vert}_{2}}{{\Vert {u}_{r}\Vert}_{2}}$, where *u _{r}* is the exact field and

*u*is the numerical field with the same step size

*h*.

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

*m*= (

_{j}*j*− 1/2)

*π*,

*y*

_{0}=

*z*

_{0}= 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*_{2}-err (the relative *L*_{2} error) is $E(u)=\frac{{\Vert {u}_{r}-u\Vert}_{2}}{{\Vert {u}_{r}\Vert}_{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.

#### 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*_{0}(*y*, *z*) = cos(2*πy*) cos(2*πz*) corresponding to the third propagating mode. Then the problem has an exact solution

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.

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.

## 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*^{6}) 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*^{3}) for each range step where *m* = *ratio* · *n*^{2} (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]