Solving the Helmholtz equation in conformal mapped ARROW structures using homotopy perturbation method

Kasper Reck; Erik V. Thomsen; Ole Hansen

doi:10.1364/OE.19.001808

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

\nabla^{2} ϕ (r) + (n_{r}^{2} k_{0}^{2} - β^{2}) ϕ (r) = 0,

where ϕ is the electromagnetic field, n_r is the refractive index, k₀ 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 ARROWs 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.

2. Theory

The mode field ϕ(x,y) for an ARROW structure is described by the Helmholtz equation

\frac{\partial^{2} ϕ (x, y)}{\partial x^{2}} + \frac{\partial^{2} ϕ (x, y)}{\partial y^{2}} + λ^{2} ϕ (x, y) = 0,

where

λ^{2} = n_{r}^{2} k_{0}^{2} - β^{2}

. Due to the antiresonance condition we shall assume that the boundary condition is

ϕ (Γ) = 0,

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 full-vectorial 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.

Download Full Size | PDF

2.1. Conformal mapping

The conformal mapping function for mapping a half-plane onto a polygon according to Schwarz-Christoffel is obtained from the integral

z (w) = A \int \prod_{j = 1}^{n} {(w - a_{j})}^{(α_{j} / π) - 1} d w,

where A is a scaling n is factor, 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]

\begin{array}{l} χ (w) - χ (0) & = & C_{sq} \int \frac{d w}{\sqrt{1 - w^{2}} \sqrt{1 - k_{r}^{2} w^{2}}} \\ = & \frac{1}{2 K (k_{r})} \int_{0}^{w} \frac{d ϑ}{\sqrt{1 - ϑ^{2}} \sqrt{1 - k_{r}^{2} ϑ^{2}}} = \frac{arcsn (w, k_{r})}{2 K (k_{r})}, \end{array}

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 since χ_E – χ_A = C_sq × iK′(k_r) = iK′(k_r) / [2K(k_r)] = i, where

K^{'} (k_{r}) = K (k_{r}^{'}) = K (\sqrt{1 - k_{r}^{2}})

. It follows that k_r ≃ 0.17157 is required for an aspect ratio of 1. The inverse mapping function is then simply

w (χ) = sn ([2 χ - 1] K (k_{r}), k_{r}) .

Returning now to the mapping of the rib waveguide in the z-plane onto the upper half w-plane, we see from Figure 1 that the mapping function is the integral [13]

\begin{array}{l} z (w) & = & C \int \frac{\sqrt{1 - k^{2} w^{2}}}{(1 - k_{1}^{2} w^{2}) \sqrt{1 - w^{2}}} d w \\ = & s \frac{k_{1}}{π} \frac{\sqrt{1 - k_{1}^{2}}}{\sqrt{k^{2} - k_{1}^{2}}} \int_{0}^{w} \frac{\sqrt{1 - k^{2} ϑ^{2}}}{(1 - k_{1}^{2} ϑ^{2}) \sqrt{1 - ϑ^{2}}} d ϑ, \end{array}

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₁ 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₁ = ksn(α), and using a transform of the integration variable according to τ = sn(ϑ) and thus dτ = cn(ϑ)dn(ϑ)dϑ the integral becomes

\begin{array}{l} z (ζ) & = & \frac{s}{π} \frac{sn (α) dn (α)}{cn α} \int_{0}^{ζ} \frac{1 - k^{2} {sn}^{2} τ}{1 - k_{1}^{2} {sn}^{2} τ} d τ \\ = & \frac{s}{π} (\frac{sn (α) dn (α)}{cn α} ζ - Π_{J} (ζ, α)), \end{array}

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 separating real and imaginary parts [13]

\frac{g}{s} = \frac{K (k)}{π} (\frac{sn (α) dn (α)}{cn α} - Z (α, k))

\frac{t}{s} = \frac{2 K^{'} (k)}{π} (\frac{sn (α) dn (α)}{cn α} - Z (α, k)) - \frac{α}{K (k)},

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]

\nabla_{x, y}^{2} ϕ (x, y) = \nabla_{υ, ω}^{2} ψ (υ, ω) {| \frac{d χ}{d z} |}^{2},

hence the partial differential equation to be solved in the mapped domain, where χ = υ + iω, is

\frac{\partial^{2} ψ (υ, ω)}{\partial υ^{2}} + \frac{\partial^{2} ψ (υ, ω)}{\partial ω^{2}} + {| \frac{d χ}{d z} |}^{- 2} λ^{2} ψ (υ, ω) = 0 .

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

\frac{d χ}{d z} = \frac{d χ}{d w} \frac{d w}{d z} = \frac{1}{2 K (k_{r})} \frac{π}{s} \frac{cn α}{sn (α) dn (α)} \frac{1 - k_{1}^{2} w^{2}}{\sqrt{1 - w^{2} k_{r}^{2}} \sqrt{1 - k^{2} w^{2}}},

which is indeed non trivial. Note, using Equation 6 the derivative may be expressed as a function of χ.

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

𝒜 (u) - g (r) = 0, r \in Ω,

with boundary conditions

𝒝 (u, \frac{\partial u}{\partial n}) = 0, r \in Γ,

where 𝒜 is a general differential operator, 𝒝 is a boundary operator, g(r) is an analytic function and Γ is the boundary of the domain Ω. Since 𝒜 in general can be separated in a linear and nonlinear operator, ℒ and 𝒩 respectively, Equation 14 can be written

ℒ (u) + 𝒩 (u) - g (r) = 0 .

Using the concept of homotopy from topology, i.e. a continuous transformation of one function to another, one can setup a homotopy equation

ℋ (v, p) = (1 - p) {ℒ (v) - ℒ (u_{0})} + p {A (v) - g (r)} = 0,

where p ∈ [0, 1] is an embedding parameter, u₀ is an initial approximation satisfying the boundary conditions and

v = \sum_{n = 0}^{\infty} v_{n} p^{n} .

The solution of Equation 14 then is

u = lim_{p \to 1} v = \sum_{n = 0}^{\infty} v_{n} .

For Equation 12 the following homotopy can be constructed

ℋ = (1 - p) (\frac{\partial^{2} ψ}{\partial υ^{2}} + \frac{\partial^{2} ψ}{\partial ω^{2}} - \frac{\partial^{2} ψ_{0}}{\partial υ^{2}} - \frac{\partial^{2} ψ_{0}}{\partial ω^{2}}) + p (\frac{\partial^{2} ψ}{\partial υ^{2}} + \frac{\partial^{2} ψ}{\partial ω^{2}} + {| \frac{d χ}{d z} |}^{- 2} λ^{2} ψ) = 0,

which is equal to Equation 12 for p = 1. By identifying terms of identical powers of p we get the following set of equations

\begin{array}{l} p^{0} : & ψ_{0} \\ p^{1} : & \frac{\partial^{2} ψ_{0}}{\partial υ^{2}} + \frac{\partial^{2} ψ_{0}}{\partial ω^{2}} + {| \frac{d χ}{d z} |}^{- 2} λ^{2} ψ_{0} + \frac{\partial^{2} ψ_{1}}{\partial υ^{2}} + \frac{\partial^{2} ψ_{1}}{\partial ω^{2}} = 0 \\ p^{2} : & {| \frac{d χ}{d z} |}^{- 2} λ^{2} ψ_{1} + \frac{\partial^{2} ψ_{2}}{\partial υ^{2}} + \frac{\partial^{2} ψ_{2}}{\partial ω^{2}} = 0 \\ p^{3} : & {| \frac{d χ}{d z} |}^{- 2} λ^{2} ψ_{2} + \frac{\partial^{2} ψ_{3}}{\partial υ^{2}} + \frac{\partial^{2} ψ_{3}}{\partial ω^{2}} = 0 \\ ⋮ \\ p^{n} : & {| \frac{d χ}{d z} |}^{- 2} λ^{2} ψ_{n - 1} + \frac{\partial^{2} ψ_{n}}{\partial υ^{2}} + \frac{\partial^{2} ψ_{n}}{\partial ω^{2}} = 0 \end{array}

where the solution is

ψ = \sum_{n = 0}^{\infty} ψ_{n} p^{n}

for p → 1. Since ψ₀ is the initial guess and thus known what remains is to solve an infinite series of partial differential equations of the form

\nabla_{υ, ω}^{2} ψ_{n} (υ, ω) = h_{n} (υ, ω),

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.

ψ_{n} (υ, ω) = \sum_{j = 1}^{\infty} \sum_{m = 1}^{\infty} E_{m j} sin (\frac{m π}{a} υ) sin (\frac{j π}{b} ω),

where the expansion coefficients, E_mj, may be determined by inserting Equation 24 into the Poisson equation and use the orthogonality relations

\int_{0}^{a} sin (\frac{m π υ}{a}) sin (\frac{q π υ}{a}) d υ = \frac{a}{2} δ_{m q}

\int_{0}^{b} sin (\frac{j π ω}{b}) sin (\frac{r π ω}{b}) d ω = \frac{b}{2} δ_{j r} .

It follows that

E_{m j} = \frac{- 4}{a b κ_{m j}} \int_{0}^{b} \int_{0}^{a} h_{n} (υ, ω) sin (\frac{m π}{a} υ) sin (\frac{j π}{b} ω) d υ d ω

where the coefficient κ_mj is given by

κ_{m j} = {(\frac{m π}{a})}^{2} + {(\frac{j π}{b})}^{2}, m, j \in [1, 2, \dots] .

We now have a solution of the Helmholtz equation in the form of an infinite series of solutions to the Poisson equation. All that remains is to choose an initial guess. The initial guess has to fulfill the boundary conditions (Equation 3) and a reasonable choice would be the first eigenfunction to the constant coefficient Helmholtz eigenvalue problem (Equation 2) i.e.

ψ_{0} (υ, ω) = sin (\frac{π}{a} υ) sin (\frac{π}{b} ω) .

2.3. 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¹) is obtained with the source function h₁ (υ,ω) = −(∇²ψ₀ + |dχ/dz|⁻²λ²ψ₀) = −(|dχ/dz|⁻²λ² – 2π²) ψ₀ and the result is

\begin{array}{l} ψ_{1} (υ, ω) & = & \sum_{j = 1}^{\infty} \sum_{m = 1}^{\infty} \frac{4 sin (m π υ) sin (j π ω)}{κ_{m j}} \\ \times & \int_{0}^{1} \int_{0}^{1} ({| \frac{d χ}{d z} |}^{- 2} λ^{2} - 2 π^{2}) sin (π υ) sin (π ω) sin (m π υ) sin (j π ω) d υ d ω . \end{array}

All following terms are obtained by repeated application of Equation 24 on the next equations in 21, thus

\begin{array}{l} ψ_{n} (υ, ω) & = & \sum_{j = 1}^{\infty} \sum_{m = 1}^{\infty} \frac{4 sin (m π υ) sin (j π ω)}{κ_{m j}} \\ \times & \int_{0}^{1} \int_{0}^{1} ({| \frac{d χ}{d z} |}^{- 2} λ^{2} ψ_{n - 1} (υ, ω)) sin (m π υ) sin (j π ω) d υ d ω, \end{array}

for n > 1.

By use of Equation 12 the eigenvalue λ² is found from a ratio of integrals over the model domain

λ^{2} = - \frac{\int ψ \nabla_{ν, ω}^{2} ψ d Ω}{\int {| \frac{d χ}{d z} |}^{- 2} ψ^{2} d Ω},

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₀, corresponding to the physical domain z = r exp (iθ) with r_a ≥ r ≥ r₀ and π ≥ θ ≥ 0. Here the mapping functions $χ (z) = ln \frac{z}{r_{0}} = ln \frac{r}{r_{0}} + i θ$ , or z(χ) = r₀ exp χ maps the physical domain on the rectangle $ln \frac{r_{a}}{r_{0}} \geq υ \geq 0$ and π ≥ ω ≥ 0. The Jacobian for this problem becomes ${| d χ / d z |}^{- 2} = r_{0}^{2} 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.

3. Results

The Jacobian |dχ/dz|⁻² 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.

Fig. 2 The Jacobian of the square to rib waveguide conformal mapping. The vertical axis of the plot has been truncated at 15.

Download Full Size | PDF

While the complete solution can be built from Equation 21, the solutions become rather intractable 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 λ² 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.

Fig. 3 Eigenvalue λ² 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 with the number of Fourier terms as parameter. Calculations are shown for 2, 4, 6 and 12 Fourier terms and compared to the result from FEM.

Download Full Size | PDF

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.

Fig. 4 The initial guess and first three solutions to the homotopy equations (Equation 21). The effect of the two peaks of the Jacobian close to the real axis is clearly seen in the HPM solutions.

Download Full Size | PDF

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 Jaco-bian 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.

Fig. 5 Sixth order homotopy solution in model domain (a), physical domain using HPM (b) and using FEM (c) for t/s = 1 and g/s = 1.

Download Full Size | PDF

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.

Fig. 6 Sixth order homotopy solution in model domain (a), physical domain using HPM (b) and using FEM (c) for t/s = 2 and g/s = 1.

Download Full Size | PDF

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.

Fig. 7 Second eigenfunction for sixth order homotopy solution in model domain (a), physical domain using HPM (b) and using FEM (c) for t/s = 1 and g/s = 1.

Download Full Size | PDF

The relative deviation ɛ between the HPM solution ψ_HPM and the FEM solution ψ_FEM is quantified using the L₂ norm as

ɛ = \frac{{‖ ψ_{HPM} - ψ_{FEM} ‖}_{L_{2}}}{{‖ ψ_{FEM} ‖}_{L_{2}}} = {[\frac{\int {(ψ_{HPM} - ψ_{FEM})}^{2} d Ω}{\int ψ_{FEM}^{2} d Ω}]}^{\frac{1}{2}} .

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₂ 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 $\frac{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.

Table 1. Eigenvalues for 6^th order HPM $λ_{HPM}^{2}$ compared to eigenvalues from FEM $λ_{FEM}^{2}$ for two different rib waveguides and the two first modes. The 2^nd mode eigen-values are based on a higher order initial guess. The relative L₂ norm deviations ɛ of 6^th order HPM solutions from the FEM solutions are also listed. In the calculations 12 Fourier terms were used.

View Table | View all tables in this article

The convergence of the eigenvalue λ² 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.

Fig. 8 The eigenvalue λ² and the HPM field amplitude ψ(χ_ref) for increasing number of homotopy terms, where χ_ref is a specific reference point in the square domain (here χ_ref = 1/2 + i3/5). The changes in eigenvalue and field amplitude are plotted on the right axis.

Download Full Size | PDF

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.

Table 2. Comparison of the performance (memory used and CPU time) of FEM and conformal mapping/homotopy (HPM) solutions.

View Table | View all tables in this article

4. Conclusion

We have shown that the general Helmholtz equation may be solved using a combination of conformal mapping and homotopy perturbation. The method was applied to a general rib waveguide structure and field distributions of fundamental and higher order modes as well as eigenvalues were calculated. The results were verified by comparison to finite element method solutions and by convergence analysis. While the mathematical framework was applied to the case of a rib waveguide, the method is very general and may easily be applied to several other interesting waveguide geometries by appropriate choice of conformal mapping functions. The method may be applied to other waveguide structures than ARROW’s, e.g. metal-cladding waveguides, as long as a null boundary condition is a valid approximation for these waveguides.

A. Appendix: Special functions

The notation used with elliptic integrals and elliptic functions differs somewhat in literature [13, 18, 19, 20, 21], thus we briefly list the notation used here.

Jacobi’s elliptic functions

Jacobi’s elliptic sine amplitude function sn(u, k) = sn(u) and Jacobi’s amplitude function am(u, k) = am(u) of modulus k are implicitly defined from an elliptic integral of the first kind; let

u = arcsn (z, k) = \int_{0}^{z} \frac{d t}{\sqrt{(1 - t^{2}) (1 - k^{2} t^{2})}} = F (ϕ, k), with z = sin ϕ

then

am (u, k) = am (u) = ϕ,

sn (u, k) = sn (u) = z = sin ϕ = sin (am (u)),

where the modulus k may be omitted if no misunderstanding is possible. Note, Equation 34 also defines the inverse of the Jacobi elliptic sine amplitude, arcsn(z, k). Jacobi’s elliptic cosine amplitude function cn(u, k) = cn(u) and delta amplitude dn(u,k) = dn(u) may be defined from

cn (u, k) = cn (u) = cos (am (u)),

dn (u, k) = dn (u) = \frac{\partial am (u)}{\partial u} .

The Jacobi elliptic functions fulfil

{sn}^{2} (u) + {cn}^{2} (u) = 1, and k^{2} {sn}^{2} (u) + {dn}^{2} (u) = 1 .

Elliptic integrals

Legendres elliptic integral of the first kind F(ϕ,k) is the definite integral

F (ϕ, k) = \int_{0}^{ϕ} \frac{d θ}{\sqrt{1 - k^{2} {sin}^{2} θ}} = \int_{0}^{sin ϕ} \frac{d t}{\sqrt{(1 - t^{2}) (1 - k^{2} t^{2})}},

F_{1} (z, k) = \int_{0}^{z} \frac{d t}{\sqrt{(1 - t^{2}) (1 - k^{2} t^{2})}}, with z = sin ϕ,

where k is the modulus, ϕ the amplitude and z = sinϕ, and F₁(z, k) is the elliptic integral in Jacobi form. The complete elliptic integral of the first kind K(k) equals the elliptic integral of the first kind at an amplitude of ϕ = π/2, thus

K = K (k) = \int_{0}^{π / 2} \frac{d θ}{\sqrt{1 - k^{2} {sin}^{2} θ}} = \int_{0}^{1} \frac{d t}{\sqrt{(1 - t^{2}) (1 - k^{2} t^{2})}},

where the modulus may be omitted, if no misunderstanding is possible. With the complementary modulus

k^{'} \equiv \sqrt{1 - k^{2}}

the complementary complete elliptic integral of the first kind comes

K^{'} = K^{'} (k) = K (k^{'}) = K (\sqrt{1 - k^{2}}) .

Legendres elliptic integral of the second kind E(ϕ,k) is defined from the definite integral

E (ϕ, k) = \int_{0}^{ϕ} \sqrt{1 - k^{2} {sin}^{2} θ} d θ = \int_{0}^{sin ϕ} \frac{\sqrt{1 - k^{2} t^{2}}}{\sqrt{1 - t^{2}}} d t

E_{1} (z, k) = \int_{0}^{z} \frac{\sqrt{1 - k^{2} t^{2}}}{\sqrt{1 - t^{2}}} d t, with z = sin ϕ,

again the complete elliptic integral E(k) of the second kind is obtained at an amplitude of ϕ = π/2

E = E (k) = E (\frac{π}{2}, k) = \int_{0}^{π / 2} \sqrt{1 - k^{2} {sin}^{2} θ} d θ = \int_{0}^{1} \frac{\sqrt{1 - k^{2} t^{2}}}{\sqrt{1 - t^{2}}} d t .

Jacobi’s elliptic integral of the third kind Π_J (z, α, k) is the definite integral

Π_{J} (z, α, k) = k^{2} sn (α, k) cn (α, k) dn (α, k) \int_{0}^{z} \frac{{sn}^{2} (u, k)}{1 - k^{2} {sn}^{2} (α, k) {sn}^{2} (u, k)} d u

which differs from Legendres elliptic integral of the third kind Π(z,α,k), since that integral in Jacobi’s form is defined as

\begin{array}{l} Π (z, α, k) & = & \int_{0}^{z} \frac{d u}{1 - k^{2} {sn}^{2} (α, k) {sn}^{2} (u, k)} \\ = & \int_{0}^{z} \frac{d t}{(1 - k^{2} {sn}^{2} (α, k) t^{2}) \sqrt{1 - t^{2}} \sqrt{1 - k^{2} t^{2}}} . \end{array}

It follows that the two elliptic integrals of the third kind are related

Π (z, α, k) = z + \frac{sn (α, k)}{cn (α, k) dn (α, k)} Π_{J} (z, α, k) .

Jacobi’s zeta function Z(u, k) = Z(u) is related to the incomplete elliptic integrals of the first and second kind by

Z (u, k) = Z (u) = E_{1} (u, k) - F_{1} (u, k) \frac{E (k)}{K (k)} .

References and links

1. R. P. Ratowsky, J. Fleck, and M. D. Feit, “Helmholtz beam propagation in rib waveguides and couplers by iterative Lanczos reduction,” J. Opt. Soc. Am. A 9, 265–273 (1992). [CrossRef]

2. M. Balagangadhar, T. Sarkar, I. Rejeb, and R. Boix, “Solution of the general Helmholtz equation in homogeneously filled waveguides using a static Green’s function,” IEEE Trans. Microwave Theory Tech. 46, 302–307 (1998). [CrossRef]

3. W. Ng and M. Stern, “Analysis of multiple-rib waveguide structures by the discrete-spectral-index method,” in Proceedings of IEEE Conference on Optoelectronics (IEEE, 1998), 365–371 (1998). [CrossRef]

4. C. T. Shih and S. Chao, “Simplified numerical method for analyzing TE-like modes in a three-dimensional circularly bent dielectric rib waveguide by solving two one-dimensional eigenvalue equations,” J. Opt. Soc. Am. B 25, 1031–1037 (2008). [CrossRef]

5. M. A. Duguay, Y. Kokubun, T. L. Koch, and L. Pfeiffer, “Antiresonant reflecting optical waveguides in SiO₂-Si multilayer structures,” Appl. Phys. Lett. 49, 13–15 (1986). [CrossRef]

6. T. Baba and Y. Kokubun, “Dispersion and radiation loss characteristics of antiresonant reflecting optical waveguides—Numerical results and analytical expressions,” IEEE J. Quantum Electron. 28, 1689–1700 (1992). [CrossRef]

7. D. Yin, J. P. Barber, A. R. Hawkins, and H. Schmidt, “Low-loss integrated optical sensors based on hollow-core ARROW waveguide,” Proc. SPIE 5730, 218–225, (2005). [CrossRef]

8. D. Yin, D. W. Deamer, H. Schmidt, J. P. Barber, and A. R. Hawkins, “Integrated optical waveguides with liquid cores,” Appl. Phys. Lett. 85, 3477–3479 (2004). [CrossRef]

9. D. Yin, H. Schmidt, J. P. Barber, E. J. Lunt, and A. R. Hawkins, “Optical characterization of arch-shaped ARROW waveguides with liquid cores,” Opt. Express 13, 10564–10570 (2005). [CrossRef] [PubMed]

10. A. M. Young, C. L. Xu, W. Huang, and S. D. Senturia, “Design and analysis of an ARROW-waveguide-based silicon pressure transducer,” Proc. SPIE 1793, 42–53 (1993). [CrossRef]

11. K. J. Rowland, S. V. Afshar, and T. M. Monro, “Bandgaps and antiresonances in integrated-ARROWs and Bragg fibers; a simple model,” Opt. Express 16, 17935–17951 (2008). [CrossRef] [PubMed]

12. J.-L. Archambault, R. J. Black, S. Lacroix, and J. Bures, “Loss calculations for antiresonant waveguides,” J. Lightwave Technol. 11, 416–423 (1993). [CrossRef]

13. W. J. Gibbs, Conformal Transformations in Electrical Engineering (Chapman & Hall, 1958).

14. R. Schinzinger and P. A. A. Laura, Conformal Mapping: Methods and Applications (Elsevier, 1991).

15. C. Lee, M. Wu, and J. Hsu, “Beam propagation analysis for tapered waveguides: taking account of the curved phase-front effect in paraxial approximation,” J. Lightwave Technol. 15, 2183–2189 (1997). [CrossRef]

16. S. Liao, “An approximate solution technique not depending on small parameters: a special example,” Int. J. Non-Linear Mech. 30, 371–380 (1995). [CrossRef]

17. J. He, “Homotopy perturbation technique,” Comput. Methods Appl. Mech. Eng. 178, 257–262 (1999). [CrossRef]

18. I. S. Gradshteyn and I. M. Ryzhik, Tables of Integrals, Series, and Products, Corrected and Enlarged Edition (Academic, 1980).

19. M. Abramowitz and I. A. Stegun, eds., Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Dover, 1972).

20. M. Hazewinkel, ed., Encyclopaedia of Mathematics, Springer online Reference Works, http://eom.springer.de/default.htm.

21. NIST Digital Library of Mathematical Functions, http://dlmf.nist.gov/.

Mode	$λ_{HPM}^{2}$	$λ_{FEM}^{2}$	ɛ	t/s	g/s
1	12.16	12.13	0.29%	1	1
1	6.788	6.789	0.14%	2	1
2	19.06	18.84	1.78%	1	1
2	13.99	13.96	0.75%	2	1

Points	Memory [MB]		CPU time [s]
Points	FEM	HPM	FEM	HPM
12417	40	7	0.9	1.1
56000	65	30	4.4	4.9
109857	173	45	13.3	9.9
143265	208	59	16.5	12.8
172257	254	73	22.8	16.0

Mode	$λ_{HPM}^{2}$	$λ_{FEM}^{2}$	ɛ	t/s	g/s
1	12.16	12.13	0.29%	1	1
1	6.788	6.789	0.14%	2	1
2	19.06	18.84	1.78%	1	1
2	13.99	13.96	0.75%	2	1

Points	Memory [MB]		CPU time [s]
Points	FEM	HPM	FEM	HPM
12417	40	7	0.9	1.1
56000	65	30	4.4	4.9
109857	173	45	13.3	9.9
143265	208	59	16.5	12.8
172257	254	73	22.8	16.0

Solving the Helmholtz equation in conformal mapped ARROW structures using homotopy perturbation method

Abstract

1. Introduction

2. Theory

2.1. Conformal mapping

2.2. Homotopy perturbation

2.3. Special cases

3. Results

4. Conclusion

A. Appendix: Special functions

Jacobi’s elliptic functions

Elliptic integrals

References and links

Cited By

Figures (8)

Tables (2)

Equations (50)

Optics Express