## Higher-order wide-angle split-step spectral method for non-paraxial beam propagation |

Optics Express, Vol. 21, Issue 13, pp. 15815-15825 (2013)

http://dx.doi.org/10.1364/OE.21.015815

Acrobat PDF (4812 KB)

### Abstract

We develop a higher-order method for non-paraxial beam propagation based on the wide-angle split-step spectral (WASSS) method previously reported [**41**, 849 (2010)

© 2013 OSA

## 1. Introduction

1. G. R. Hadley, “Multistep method for wide-angle beam propagation,” Opt. Lett. **17**, 1743–1745 (1992) [CrossRef] [PubMed] .

2. K. Q. Le, R. Godoy-Rubio, P. Bienstman, and G. R. Hadley, “The complex Jacobi iterative method for three-dimensional wide-angle beam propagation,” Opt. Express **16**, 17021–17030 (2008) [CrossRef] [PubMed] .

3. Y. Y. Lu and P. L. Ho, “Beam propagation method using a [(*p*−1)/*p*] Padé approximant of the propagator,” Opt. Lett. **27**, 683–685 (2002) [CrossRef] .

4. A. Sharma and A. Agrawal, “New method for nonparaxial beam propagation,” J. Opt. Soc. Am. B **21**, 1082–1087 (2004) [CrossRef] .

4. A. Sharma and A. Agrawal, “New method for nonparaxial beam propagation,” J. Opt. Soc. Am. B **21**, 1082–1087 (2004) [CrossRef] .

5. A. Sharma and A. Agrawal, “Non-paraxial split-step finite-difference method for beam propagation,” Opt. Quantum. Electron. **38**, 19–34 (2006) [CrossRef] .

6. C. D. Clark and R. Thomas, “Wide-angle split-step spectral method for 2D or 3D beam propagation,” Opt. Quantum. Electron. **41**, 849–857 (2010) [CrossRef] .

4. A. Sharma and A. Agrawal, “New method for nonparaxial beam propagation,” J. Opt. Soc. Am. B **21**, 1082–1087 (2004) [CrossRef] .

## 2. Formulation

*ψ*(

*z*,

**r**) is the complex scalar electric field,

*z*is a Cartesian coordinate in the direction of propagation,

**r**are the transverse coordinates, and

*k*

_{0}is the free space wavenumber,

*n*(

*z*,

**r**) is the refractive index distribution, and

*n̄*

^{2}= min [

*n*

^{2}(

*z*,

*x*)]. For simplicity we will only consider the two-dimensional (2D) case, however generalization to three dimensions (3D) follows quite trivially. In two dimensions, the electric field is denoted

*ψ*(

*z*,

*x*). Note that

*x*is not necessarily a Cartesian coordinate but could, for example, be a radial coordinate. Expanding

*ψ*(

*z*,

*x*) in terms of the eigenfunctions of the transverse Laplace operator we write where Multiplying the 2D form of Eq. (1) by

*N*and

_{x}*N*be the number of discrete grid points in x and z respectively, and let

_{z}*i*,

*j*∈ [1,

*N*], and

_{x}*l*∈ [1,

*N*]. We let

_{z}*x*=

_{i}*i*Δ

*x*+

*x*

_{0}and

*z*=

_{l}*l*Δ

*z*+

*z*

_{0}, where Δ

*x*and Δ

*z*are the corresponding step sizes. The matrix

**N**(

*z*) is defined by The discrete eigentransform matrix,

**S**[6

6. C. D. Clark and R. Thomas, “Wide-angle split-step spectral method for 2D or 3D beam propagation,” Opt. Quantum. Electron. **41**, 849–857 (2010) [CrossRef] .

**S**

^{−1}transforms vectors back into the spatial basis. In the case of Cartesian coordinates with hard boundary conditions this matrix becomes the discrete Fourier transform matrix, and in cylindrical coordinates with azimuthal symmetry and hard boundary conditions,

**S**takes the form of a discrete Hankel transform [7

7. M. Guizar-Sicairos and J. C. Gutiérrez-Vega, “Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields,” J. Opt. Soc. Am. A **21**, 53–58 (2004) [CrossRef] .

**M**with elements Notice that if

*λ*>

_{i}*k*

_{0}

*n̄*then

*M*becomes imaginary. If we allow eigenvalues such that

_{i,j}*λ*>

_{i}*k*

_{0}

*n̄*, then the trigonometric functions appearing in

**P**will become hyperbolic (see Eq. (25)), making the method numerically unstable. For this reason, we will limit

*N*to satisfy

_{x}*λ*<

_{i}*k*

_{0}

*n̄*. However, we must include at least enough eigenfrequencies to sufficiently represent the initial conditions, otherwise large errors will be incurred. Using these definitions Eq. (5) becomes Now if we define we can write the Helmholtz equation as a first order vector differential equation

8. W. Magnus, “On the exponential solution of differential equations for a linear operator,” Comm. Pure Appl. Math. **7**, 649–673 (1954) [CrossRef] .

**H**(

*t*

_{1}),

**H**(

*t*

_{2})] is the usual commutator. Note that Eq. (11) is essentially the time-dependent Schrödinger equation, in which context the Magnus expansion is more well-known as the time-ordered exponential operator [9

9. M. Bauer, R. Chetrite, K. Ebrahimi-Fard, and F. Patras, “Time-ordering and a generalized Magnus expansion,” Lett. Math. Phys. **103**, 331–350 (2012) [CrossRef] .

**H**(

*z*) can be written as the sum of two matrices, one that is constant and one that depends on

*z*, Using this and the trapezoid rule to approximate the integral in Eq. (13) we obtain, To approximate Ω

_{2}(

*z*

_{l}_{+1},

*z*) we first need to compute the commutator. After making use of Eq. (15) and simplifying we find Note that [

_{l}**H**

_{2}(

*z*

_{1}),

**H**

_{2}(

*z*

_{2})] = 0. We apply the trapezoid method twice to approximate the double integral in Eq. (14) and obtain The next term in the Magnus series will contribute a factor of (Δ

*z*)

^{3}as it involves a triple integral. Keeping terms up to (Δ

*z*)

^{3}Because this contains the exponential of a dense matrix, which is difficult to handle numerically, we will split this exponential into a form that will be easier to work with via the symmetric exponential splitting technique First, we split the Δ

*z*terms in the exponential from the (Δ

*z*)

^{2}terms. Then we split the terms involving

**H**

_{1}from those involving

**H**

_{2}, and we are left with where Each operator

**P**,

**Q**(

*z*), and

**C**(

*z*

_{1},

*z*

_{2}) can be expressed by a 2 × 2 block matrix where the matrices within each operator are diagonal. This allows the exponential operators to be simply evaluated by expanding the exponential into its power series Now all matrices inside of functions are diagonal making them quite simple to numerically evaluate. In this form the physical interpretation of these operators is most transparent.

**P**is the operator corresponding to a propagation of the pulse through a homogenous medium with an index of refraction of

*n̄*.

**Q**(

*z*) takes into account the difference

*n*(

*z,x*) −

*n̄*, while

*C*(

*z*

_{1},

*z*

_{2}) depends on the change in the refractive index over the step taken. Thus, in this sense we would like to draw the analogy that

*Q*(

*z*) acts like the constant term of a Taylor’s series expansion, and

*C*(

*z*

_{1},

*z*

_{2}) acts in a manner similar to the first derivative term in such a series. We can save some additional matrix-vector multiplications by combining

**Q**(

*z*

_{l}_{+1})

**Q**(

*z*)

_{l}## 3. Numerical example

*θ*from the positive

*z*axis [5

5. A. Sharma and A. Agrawal, “Non-paraxial split-step finite-difference method for beam propagation,” Opt. Quantum. Electron. **38**, 19–34 (2006) [CrossRef] .

*ψ*(

*z*,

*x*

_{0}) = 0 and

*ψ*(

*z*,

*x*) = 0 are assumed. The eigenfunctions of the transverse Laplace operator are then given by The resulting eigentransform matrix is given by the discrete Fourier matrix This particular Fourier matrix was chosen because the hard boundary conditions at

_{f}*x*

_{0}and

*x*are automatically satisfied. We do not need to include these points in our grid, defined by Δ

_{f}*x*= (

*x*−

_{f}*x*

_{0})/(

*N*+ 1) and Δ

_{x}*z*= (

*z*−

_{f}*z*

_{0})/

*N*. The refractive index profile is where

_{z}*x̃*is a shifted coordinate to align the waveguide in the center of our computational region to make the hard boundary conditions as negligible as possible, Δ

*n*is the height of the refractive index shift,

*w*is the width of the waveguide. The shifted coordinate is given by

*x̃*=

*x*− (1/2)(

*x*−

_{f}*x*

_{0}) + (1/2)(

*z*−

_{f}*z*

_{0}) tan(

*θ*). The initial electric field is given by the ze-roth order mode of the Epstein-layer waveguide Here,

*i*is the unit imaginary number, not to be confused with the counting index used elsewhere.

*W*and

*K*

_{0}are given by The exact solution for this tilted Epstein-layer waveguide is [10] To measure the error we use the correlation factor, as this provides a measure for both the profile shape and amplitude of the beam [5

5. A. Sharma and A. Agrawal, “Non-paraxial split-step finite-difference method for beam propagation,” Opt. Quantum. Electron. **38**, 19–34 (2006) [CrossRef] .

6. C. D. Clark and R. Thomas, “Wide-angle split-step spectral method for 2D or 3D beam propagation,” Opt. Quantum. Electron. **41**, 849–857 (2010) [CrossRef] .

*k*

_{0}was selected to ensure that we could include 1000 transverse modes while still keeping

**M**real-valued. When the waveguide was aligned (

*θ*= 0°) roughly an order of magnitude decrease in the error was observed with the HOWASSS method over the WASSS method (see Fig. 1). However, the improvement was more pronounced when we tilted the waveguide at an angle of 50 degrees, especially for large step sizes (see Fig. 2). Figure 3 shows the calculated electric field next to the exact solution to further illustrate the accuracy of the HOWASSS method. The simulations shown in Figs. 1, 2, and 3 were run using double precision accuracy numbers on the GPU. Both methods were capable of producing accurate results using single precision arithmetic, until about Δ

*z*= 0.05

*μ*m. At this point, round-off error began to affect the results of the HOWASSS method due to the increased number of matrix-vector multiplications required in each time step. The WASSS method does not suffer as quickly from this problem because it is computationally more simple. The limit where round-off error begins to affect the double precision code was not observed for either method.

## 4. Discussion

### 4.1. Stability

*n*, while keeping the width of the waveguide fixed. We see that at a 50° waveguide tilt the WASSS method actually becomes unstable with a large stepsize, while the HOWASSS is able to remain stable for at least an additional order of magnitude increase in Δ

*n*. However, both methods do loose accuracy as the change in the refractive index becomes steeper, but by decreasing Δ

*z*accurate results can still be obtained in a reasonable run time. This same analysis was done for a 0° tilted waveguide and the results were quite similar and so are not shown here.

### 4.2. Speed

*N*so that the matrix

_{x}**M**is real then, for the example in Section 3, all the matrices will be real, however the vector

**A**will be complex. In general, the eigenfunctions and the refractive index could be complex. However, for simplicity we will assume these to be real.

*N*×

_{x}*N*diagonal matrix by an

_{x}*N*element complex vector requires 2

_{x}*N*operations. Multiplying an

_{x}*N*×

_{x}*N*dense matrix by an

_{x}*N*element complex vector is equivalent to performing 2

_{x}*N*dot products, each requiring

_{x}*N*multiplications and

_{x}*N*− 1 additions, for a total of

_{x}**P**requires a total of 4 diagonal matrix-vector multiplications for a total of 8

*N*operations. Applying

_{x}**Q**(

*z*

_{l}_{+1})

**Q**(

*z*) requires 2 diagonal matrix-vector multiplications, 2 dense matrix-vector multiplications, 2 vector-vector additions, and 1 scalar-vector multiplication for a total of

_{l}**C**(

*z*

_{l}_{+1},

*z*) requires 4 dense matrix-vector multiplications, 1 vector-vector addition, and 2 element-by-element exponential operations (assumed to only be 1 operation per element) for a total of

_{l}**P**is applied four times,

**Q**(

*z*

_{l}_{+1})

**Q**(

*z*) is applied twice, and

_{l}**C**(

*z*

_{l}_{+1},

*z*) is applied once each step, giving a total of

_{l}**41**, 849–857 (2010) [CrossRef] .

*N*, the HOWASSS method is approximately a factor of 4 slower for a given Δ

_{x}*z*(see Fig. 5). However, the initial hypothesis was that a significant improvement in accuracy may lead to a more efficient algorithm. To test this hypothesis, we executed the HOWASSS and WASSS methods using the GPU implementation with various step sizes, holding all other aspects of the problem constant. Figure 6 summarizes the result. For a propagation length of 100 micrometers, and a waveguide tilt angle of 50 degrees, we show the compute time per micron of propagation as a function of the maximum absolute error. For two different values of

*N*, the efficiency at constant error is better for the HOWASSS method, as indicated by shorter compute times. In fact, for error values smaller than about 10

_{x}^{−4}, the HOWASSS method is much better, with efficiency rapidly exceeding an order of magnitude for absolute error values of less than 10

^{−6}.

*N*is large, around 2000 for the specific example considered in this paper. If the application does not require extreme accuracy, then a substantial speed-up can be achieved by using single-precision arithmetic on the GPU.

_{z}### 4.3. Parallelization

### 4.4. Generalizations to other coordinate systems, boundary conditions, and higher dimensions

**41**, 849–857 (2010) [CrossRef] .

## 5. Conclusion

*z*)

^{3}in the Magnus expansion we use a symmetric operator splitting technique in order to analytically reduce the exponential matrices into a more simple form. The solution to the Helmholtz equation is then approximated via straightforward matrix multiplication. We have demonstrated the method in a simple geometry with 2D Cartesian coordinates with hard boundary conditions. The results obtained show our higher-order approach significantly improves the overall efficiency, when measured as compute time per distance of propagation, in cases where high accuracy results are required. The method can be easily extended to more generalized coordinate systems, higher dimensions, and various boundary conditions, provided that the eigenfunctions and eigenvalues of the transverse Laplace operator can be found for that geometry.

## Acknowledgments

## References and links

1. | G. R. Hadley, “Multistep method for wide-angle beam propagation,” Opt. Lett. |

2. | K. Q. Le, R. Godoy-Rubio, P. Bienstman, and G. R. Hadley, “The complex Jacobi iterative method for three-dimensional wide-angle beam propagation,” Opt. Express |

3. | Y. Y. Lu and P. L. Ho, “Beam propagation method using a [( |

4. | A. Sharma and A. Agrawal, “New method for nonparaxial beam propagation,” J. Opt. Soc. Am. B |

5. | A. Sharma and A. Agrawal, “Non-paraxial split-step finite-difference method for beam propagation,” Opt. Quantum. Electron. |

6. | C. D. Clark and R. Thomas, “Wide-angle split-step spectral method for 2D or 3D beam propagation,” Opt. Quantum. Electron. |

7. | M. Guizar-Sicairos and J. C. Gutiérrez-Vega, “Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields,” J. Opt. Soc. Am. A |

8. | W. Magnus, “On the exponential solution of differential equations for a linear operator,” Comm. Pure Appl. Math. |

9. | M. Bauer, R. Chetrite, K. Ebrahimi-Fard, and F. Patras, “Time-ordering and a generalized Magnus expansion,” Lett. Math. Phys. |

10. | M. J. Adams, |

**OCIS Codes**

(000.4430) General : Numerical approximation and analysis

(080.1510) Geometric optics : Propagation methods

(350.5500) Other areas of optics : Propagation

(080.1753) Geometric optics : Computation methods

**ToC Category:**

Physical Optics

**History**

Original Manuscript: March 25, 2013

Revised Manuscript: May 24, 2013

Manuscript Accepted: June 14, 2013

Published: June 25, 2013

**Citation**

Brett H. Hokr, C. D. Clark, Rachel E. Grotheer, and Robert J. Thomas, "Higher-order wide-angle split-step spectral method for non-paraxial beam propagation," Opt. Express **21**, 15815-15825 (2013)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-21-13-15815

Sort: Year | Journal | Reset

### References

- G. R. Hadley, “Multistep method for wide-angle beam propagation,” Opt. Lett.17, 1743–1745 (1992). [CrossRef] [PubMed]
- K. Q. Le, R. Godoy-Rubio, P. Bienstman, and G. R. Hadley, “The complex Jacobi iterative method for three-dimensional wide-angle beam propagation,” Opt. Express16, 17021–17030 (2008). [CrossRef] [PubMed]
- Y. Y. Lu and P. L. Ho, “Beam propagation method using a [(p−1)/p] Padé approximant of the propagator,” Opt. Lett.27, 683–685 (2002). [CrossRef]
- A. Sharma and A. Agrawal, “New method for nonparaxial beam propagation,” J. Opt. Soc. Am. B21, 1082–1087 (2004). [CrossRef]
- A. Sharma and A. Agrawal, “Non-paraxial split-step finite-difference method for beam propagation,” Opt. Quantum. Electron.38, 19–34 (2006). [CrossRef]
- C. D. Clark and R. Thomas, “Wide-angle split-step spectral method for 2D or 3D beam propagation,” Opt. Quantum. Electron.41, 849–857 (2010). [CrossRef]
- M. Guizar-Sicairos and J. C. Gutiérrez-Vega, “Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields,” J. Opt. Soc. Am. A21, 53–58 (2004). [CrossRef]
- W. Magnus, “On the exponential solution of differential equations for a linear operator,” Comm. Pure Appl. Math.7, 649–673 (1954). [CrossRef]
- M. Bauer, R. Chetrite, K. Ebrahimi-Fard, and F. Patras, “Time-ordering and a generalized Magnus expansion,” Lett. Math. Phys.103, 331–350 (2012). [CrossRef]
- M. J. Adams, An Introduction to Optical Waveguides (Wiley, 1981).

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.