Solving the Helmholtz equation in conformal mapped ARROWstructures using homotopy perturbation method

: The scalar wave equation, or Helmholtz equation, describes within a certain approximation the electromagnetic ﬁeld distribution in a given system. In this paper we show how to solve the Helmholtz equation in complex geometries using conformal mapping and the homotopy perturbation method. The solution of the mapped Helmholtz equation is found by solving an inﬁnite series of Poisson equations using two dimensional Fourier series. The solution is entirely based on analytical expressions and is not mesh dependent. The analytical results are compared to a numerical (ﬁnite element method) solution.


Introduction
During the process of optical system design it is often required to investigate the mode shape of the propagating electromagnetic field. The mode shape can be obtained by directly solving the full vectorial Maxwell's equations. However, it is either very complicated or impossible to find analytical solutions for other than the most simple geometries, hence numerical methods are employed for all practical purposes. If the refractive index contrast is small and the field is independent of the field polarization, analytical results can be obtained for a wider range of problems using the scalar wave equation where φ is the electromagnetic field, n r is the refractive index, k 0 is the vacuum wave number, β is the propagation constant and r is the position vector. The scalar wave equation has the form of the well known Helmholtz equation. In order to solve the scalar wave equation several methods have been applied within the optics community, including iterative Lanczos reduction [1], Green's functions [2], the discrete spectral-index method [3] and the beam propagation method [4]. In this paper we will apply a combination of conformal mapping and homotopy perturbation in order to directly obtain solutions of the Dirichlet Helmholtz equation. These solutions can be used for describing the modes of Antiresonant Reflecting Optical Waveguides (ARROWs) [5,6], which are leaky waveguides based on an antiresonance Fabry-Perot reflector. Since AR-ROWs rely on anti-resonans rather than total internal reflection, as conventional waveguides does, they can guide electromagnetic waves in a medium with refractive index lower than that of its surroundings. This makes them especially useful for several sensing applications, e.g. fluid optical sensing [7,8]. Since the sensing mechanism is based on an interaction with the guided electromagnetic field, it can be very practical to compute the field distribution in order to understand and optimize the sensor design. While ARROWs can have different cross-sectional geometries, most field distribution calculations are based on a relatively simple one dimensional slab waveguide model. While such a model is useful for describing fundamental properties of the ARROW guiding mechanisms (cladding layer thickness, propagation loss etc), it cannot be used for an accurate description of the electromagnetic field distribution in two dimensionally confined ARROWs, e.g. with square, half circular [9] and rib [10] cross-sections.

Theory
The mode field φ (x, y) for an ARROW structure is described by the Helmholtz equation where λ 2 = n 2 r k 2 0 − β 2 . Due to the antiresonance condition we shall assume that the boundary condition is where Γ is the boundary; this boundary condition, however, is only strictly correct at a single wavelength, but for well confined modes at wavelengths in proximity of antiresonance it is a very good approximation [11,12]. The use of the Helmholtz equation rather than the fullvectorial wave equations is valid as long as the medium is homogeneous, i.e. the refractive index is constant throughout the waveguide core, and the cladding layers are designed properly, i.e. the field is suppressed at the boundary.
In case the refractive index is constant throughout the domain and the domain boundary is not too complicated (e.g. rectangular or circular) it is straightforward to obtain an analytical solution. For non-constant refractive index or complicated domain boundaries, numerical approaches have to be applied. In the case of a rib ARROW waveguide, the domain boundary is somewhat complicated while the refractive index is constant. In order to simplify the domain boundary a Schwarz-Christoffel conformal map is applied to transform the domain from that of the rib waveguide (z-plane) to the upper complex half-plane (w-plane), see Figure 1, top panel. Since the Helmholz equation is not easily solved in the upper half-plane, a second Schwarz-Christoffel conformal map is used to map the upper half-plane onto a square in the complex χ-plane as illustrated in the lower panel of Figure 1. Fig. 1. Sketch of the rib ARROW waveguide cross-section in the z-plane and the upper half-plane w to which the rib ARROW structure is mapped conformally. Below, the square in the χ-plane to which the upper half-plane w is mapped in a second conformal mapping step.

Conformal mapping
The conformal mapping function for mapping a half-plane onto a polygon according to Schwarz-Christoffel is obtained from the integral where A is a scaling factor, n is the number of sides in the polygon, α j are the internal angles of the polygon and a j are the coordinate points on the real axis in the w-plane corresponding to which the polygon vertices in the z-plane are transformed. The second conformal map of a unit square in the χ-plane on the upper half-plane (w-plane), as sketched in the lower panel of Figure 1, is obtained by the elliptic integral of the first kind in Jacobi form [13,14] where χ(0) = 1/2, K(k r ) is the complete elliptic integral of the first kind of modulus k r , arcsn(w, k r ) is the inverse of Jacobi's elliptic sine amplitude sn(z, k r ) both of modulus k r . The scaling factor C sq = [2K(k r )] −1 is determined by the length of the base-line (A-B) of the square, since χ B − χ A = C sq × 2K(k r ) = 1. The modulus k r controls the aspect ratio of the rectangle It follows that k r ≃ 0.17157 is required for an aspect ratio of 1. The inverse mapping function is then simply Returning now to the mapping of the rib waveguide in the z-plane onto the upper half wplane, we see from Figure 1 that the mapping function is the integral [13] where we in the definite integral use that origo in the w-plane should be mapped on origo in the z-plane. The scaling factor C is determined by the requirement that the integral must increase by ∆z = is/2 when w pass 1/k 1 corresponding to the points H and B in Figure 1. As suggested by Gibbs [13] this integral can be rewritten in terms of Jacobi's incomplete elliptic integral of the third kind, Π J (ζ , α, k) (see appendix for definition), by introducing the two parameters ζ and α. Using the Jacobi elliptic functions sn, cn and dn, and defining ζ from w = sn(ζ ) and α from k 1 = ksn(α), and using a transform of the integration variable according to τ = sn(ϑ ) and thus dτ = cn(ϑ ) dn(ϑ ) dϑ the integral becomes where all Jacobi elliptic functions are to modulus k. The constants k and α are determined by the aspect ratios of the rib structure by considering mapping of the point G in Figure 1 and where Z(α, k) is Jacobi's zeta function and all elliptic function are to modulus k. Thus by using Equations 8, 9 and 10 one can easily map rib waveguides of different dimensions t, g, and s conformally to the upper complex half plane, which then may be mapped onto a unit square using Equation 5, and thus a much simpler boundary results. The price paid for the simplification of the boundary is that the Helmholtz equation becomes nonlinear as described below.
Assume that the potential in the physical z-plane φ (z) = φ (x, y) and the potential in the model χ-plane ψ (χ) = ψ (υ, ω) are related such that ψ (υ (x, y) , ω (x, y)) = φ (x, y), then the Laplace operator is affected as [14] hence the partial differential equation to be solved in the mapped domain, where Obviously, Equation 12 can be quite involved to solve, since the derivative of the conformal map except for few simple cases is very difficult to handle. The most often encountered approach in literature is therefore to apply numerical methods [15]. For the rib to square transformation considered here, the derivative is which is indeed non trivial. Note, using Equation 6 the derivative may be expressed as a function of χ.

Homotopy perturbation
A relatively new analytical method for solving nonlinear differential equations is the Homotopy Perturbation Method (HPM) [16,17]. HPM does not rely on a small parameter, as conventional perturbation methods do and has been successfully applied to a number of classic nonlinear differential equations. We consider the nonlinear differential equation with boundary conditions where A is a general differential operator, B is a boundary operator, g(r) is an analytic function and Γ is the boundary of the domain Ω. Since A in general can be separated in a linear and nonlinear operator, L and N respectively, Equation 14 can be written Using the concept of homotopy from topology, i.e. a continuous transformation of one function to another, one can setup a homotopy equation where p ∈ [0, 1] is an embedding parameter, u 0 is an initial approximation satisfying the boundary conditions and The solution of Equation 14 then is For Equation 12 the following homotopy can be constructed which is equal to Equation 12 for p = 1. By identifying terms of identical powers of p we get the following set of equations where the solution is for p → 1. Since ψ 0 is the initial guess and thus known what remains is to solve an infinite series of partial differential equations of the form that is Poisson equations where h n (υ, ω) is the source term. Since the conformal map transformed the rib waveguide into an a × b rectangle (here a unit square), the solution to the Poisson equations can be expressed as a two dimensional Fourier series, i.e.
where the expansion coefficients, E m j , may be determined by inserting Equation 24 into the Poisson equation and use the orthogonality relations It follows that where the coefficient κ m j is given by

Special cases
Hence the solutions to Equations 21 for a = b = 1 are found by substituting the Equations 29 and 24 into 21. The first term of the homotopy solution (p 1 ) is obtained with the source function h 1 (υ, ω) = − ∇ 2 ψ 0 + |dχ/dz| −2 λ 2 ψ 0 = − |dχ/dz| −2 λ 2 − 2π 2 ψ 0 and the result is for n > 1. By use of Equation 12 the eigenvalue λ 2 is found from a ratio of integrals over the model domain which can be evaluated both using the initial guess as well as using solutions including higher order terms. Equations 30 and 31 are quite general and may be applied to any conformal mapping of the Helmholtz equation to a rectangle. The rib waveguide considered here is just one such example, another simple example could be the half co-axial waveguide with outer radius r a and inner radius r 0 , corresponding to the physical domain z = r exp (iθ ) with r a ≥ r ≥ r 0 and π ≥ θ ≥ 0. Here the mapping functions χ (z) = ln z r 0 = ln r r 0 + iθ , or z (χ) = r 0 exp χ maps the physical domain on the rectangle ln r a r 0 ≥ υ ≥ 0 and π ≥ ω ≥ 0. The Jacobian for this problem becomes |dχ/dz| −2 = r 2 0 exp (2υ). For this rather simple mapping, infinite sum HPM solutions can easily be found using the method described in this paper. Other useful analytic mappings could include those of triangular or circular cross-sections, while arbitrarily shaped cross-sections could be analyzed using numerically approximated conformal transformations.

Results
The Jacobian |dχ/dz| −2 for a rib waveguide with aspect ratios t/s = g/s = 1 is shown in Figure 2 as a function of the square model domain coordinates; the Jacobian is seen to have two distinct peaks at points in the χ-plane corresponding to z → ±∞. Since the Jacobian may be considered an effective refractive index in the model domain, the mode field is expected to be attracted towards these points as is also seen below. tractable for the conformal mapping to the rib waveguide. It is therefore practical to make some simplifications rather than using the complete solution. First of all, the Homotopy perturbation method is known to converge rapidly, hence only few terms are needed for any practical application, secondly, higher order terms in the Fourier expansion solution in Equations 30 and 31 may be neglected. Figure 3 illustrates the consequences of such simplifications by showing the calculated eigenvalue λ 2 for the second mode in a rib waveguide with aspect ratios t/s = 2 and g/s = 1 as a function of the HPM order (0-6) with the number of Fourier terms (2, 4, 6, or 12) as parameter. The calculated eigenvalues are compared to the results of a MATLAB Finite Element Model (FEM) solution, where again zero field boundary conditions were assumed. Obviously, in this case the HPM order should be 5-6 and the number of Fourier terms 6-12; for the first mode 6 Fourier terms may be sufficient since the spatial spectral requirements here are less.
In Figure 4 the initial guess and the three first homotopy solutions are shown in the model domain for a rib waveguide with aspect ratios t/s = g/s = 1, and with the Fourier expansion limited to 6 th order. The effect of the Jacobian peaks is clearly visible in the HPM solutions where the mode field is shifted towards the peaks. However, it should also be noted that while the effective refractive index is approaching infinity at these peaks, the field is approaching zero due to the boundary conditions, causing finite field amplitude and a guided wave in the rib structure. Figure 5a shows a 6 th order homotopy solution in the unit square model domain for a rib waveguide with aspect ratios t/s = g/s = 1; the mode field is clearly shifted towards the Jacobian peaks and thus the mode is asymmetric in the ω-direction. In Figure 5b the HPM solution is shown mapped to the physical domain, while Figure 5c shows a Finite Element Method (FEM) solution in the physical domain for the same waveguide. Obviously, the overall field distribution of the HPM solution matches that of the FEM solution, with the maximum amplitude at the center of the waveguide (z = i/2). Since the transformation is general, it can be concluded that for all rib waveguides where the rib height is equal to the gap height (t/s = 1), the maximum amplitude is located at the center. Figure 6a shows the HPM solution in the model domain for a rib waveguide with aspect ratios t/s = 2 and g/s = 1, again the mode field is attracted towards the Jacobian peaks but the mode is less asymmetric in the ω-direction than that of Figure 5a. Figure 6b shows the HPM solution mapped to the physical domain; the solution is seen to agree well with the FEM solution in Figure 6c. The modes are seen to be shifted downwards in the rib when compared to the mode in Figure 5. Figure 7a shows the HPM solution in the model domain for the second mode of a rib waveguide with aspect ratios t/s = g/s = 1; here a higher order initial guess was used. The asymmetry of the mode in the ω-direction is similar to that in Figure 5a as expected. Figure 7b shows the HPM solution mapped to the physical domain while Figure 7c shows the FEM solution. Due to the finite width of the FEM solution domain, the field is shifted towards ±∞, compared to the HPM solution.
The relative deviation ε between the HPM solution ψ HPM and the FEM solution ψ FEM is quantified using the L 2 norm as Approximate eigenvalues found from HPM solutions using Equation 32 and eigenvalues found by solving the Helmholtz equation using the Finite Element Method (FEM) are listed in Table 1; values calculated for the first two modes for two different aspect ratios t/s = 1 or 2 and g/s = 1 are shown. The two methods are seen to agree well, especially for the first mode, while the second mode values differ more. Table 1 also lists calculated values of relative L 2 norm deviations ε, Equation 33, for the four cases. Even though the relative deviation for the second modes is larger than the relative deviation of the first modes, as also observed in the field distribution plots, the HPM solution is still within less than 2% of the FEM solution and the agreement could probably be further improved by an improved initial guess or by including higher order terms; a perfect match between the two methods would however require an infinite FEM domain. A relative deviation of less than 1 2 % for first mode is sufficient for mode-overlap calculations, while more demanding tasks related to e.g. waveguide dispersion would require higher accuracy.
The convergence of the eigenvalue λ 2 and field amplitude ψ(χ ref ) = ψ(1/2 + i3/5) of the HPM solution for a rib waveguide with t/s = 2 and g/s = 1 for increasing number of homotopy terms (iteration) is shown in Figure 8a and b, respectively. Clearly, both converge rapidly towards a finite value. In Table 2 we show some statistics (memory used and CPU time) comparing the performance of FEM and conformal mapping/homotopy solutions for 6 homotopy terms and 6 fourier terms. The homotopy solution clearly requires far less memory and is faster except for low resolution solutions. For the homotopy solutions both CPU time and RAM used seem to increase linearly with the number of points, while a super-linear trend is seen for the FEM solution.