Full-wave simulation of optical waveguides via truncation in the method of moments using PML absorbing boundary conditions

G. Karagounis; D. De Zutter; D. Vande Ginste

doi:10.1364/OE.24.028326

1. Introduction

Although in the past decade advanced numerical techniques for the solution of Maxwell’s equations have been developed, the full-wave simulation of optical structures and in particular optical waveguide structures, remains very demanding in terms of memory and CPU-time resources [1]. This is mainly due to the large electrical size of the structures under study. For design and optimization purposes it is desirable to be able to concentrate on part of the structure, but this is not easily accomplished as the interactions with other parts of the structure are most often non-negligible. A possible solution is to simplify the problem at hand, e.g. by considering the coupling of a beam into a half-infinite structure and taking advantage of a ray optical approach [2], a mode-matching technique [3], or to invoke a current windowing function in combination with the Method of Moments (MoM) [4]. The use of the MoM to obtain the scattered fields without windowing, but by truncating at such a distance from the scattering edge under study that the effect of the extraneous scattering centers becomes negligible, has also been recently investigated in [5]. Although possible, the simulated structure needs to be hundreds of wavelengths long to minimize the influence of the extraneous reflections caused by the truncation.

However, it was recognized early on that, combined with mode-matching techniques, the Perfectly Matched Layer (PML) absorbing boundary condition, first used in conjunction with the Finite Difference Time Domain (FDTD) method, could be advantageously harnessed to model optical waveguides [6–8]. In this paper, we integrate the stretching formulation of a PML in a novel way, as used in Finite Element Methods (FEM) [9], into a MoM based Boundary Integral Equation (BIE) solver. The PML layer is used to remove the influence of extraneous reflection/diffraction centers by damping the interaction between them and the part of the structure that we want to focus on.

In Sec. 2 we introduce the stretching formalism of PMLs, the system of BIEs that is solved and we explain how the stretching formalism is implemented in a MoM solver. To validate the new technique, Sec. 3 discusses the classical example of diffraction by a Perfectly Electric Conducting (PEC) wedge. The novel solution for a finite wedge, whereby the effect of the wedge truncation is eliminated by the MoM-PML combination, is compared to the exact analytical solution for the infinite wedge [10], as such allowing for a critical evaluation of the truncation effect. At this point, it is important to realize that the PML stretching formalism can be implemented in several flavors. Hence, taking advantage of the canonical wedge example, we also investigate which of these flavors yields the best results when combined with the MoM. Sec. 4 presents the application of our new technique to a Polarization Beam Splitter (PBS), the design of which is based on a prototype consisting of silicon-on-insulator (SOI) nanowires excited by a Gaussian beam [11]. In Sec. 5, we provide some concluding remarks.

In this paper, we deal with two-dimensional (2-D) problems. The z-axis is chosen as the axis of invariance. A harmonic exp(jωt) time depedence is assumed, ω being the angular frequency.

2. Formalism

2.1. BIE-MoM

Consider the configuration of Fig. 1, showing a PEC or dielectric cylinder with refractive index n₁. The pertinent transverse magnetic (TM_z) BIEs for a single dielectric object are [1]

\begin{array}{l} e_{z}^{i} & - lim_{r \to C^{+}} \oint_{C} d c^{'} [e_{z} \frac{\partial G_{0}}{\partial n^{'}} - \frac{j k_{0}^{2}}{ω ∊_{0}} G_{0} h_{t}] \\ = lim_{r \to C^{-}} \oint_{C} d c^{'} [e_{z} \frac{\partial G}{\partial n^{'}} - \frac{j k_{1}^{2}}{ω ∊_{0} n_{1}^{2}} G h_{t}], \end{array}

\begin{array}{l} h_{t}^{i} & - lim_{r \to C^{+}} \oint_{C} d c^{'} [- \frac{j ω ∊_{0}}{k_{0}^{2}} e_{z} \frac{\partial^{2} G_{0}}{\partial n \partial n^{'}} - \frac{\partial G_{0}}{\partial n} h_{t}] \\ = lim_{r \to C^{-}} \oint_{C} d c^{'} [- \frac{j ω ∊_{0} n_{1}^{2}}{k_{1}^{2}} e_{z} \frac{\partial^{2} G}{\partial n \partial n^{'}} - \frac{\partial G}{\partial n} h_{t}], \end{array}

where

e_{z}^{i}

and

h_{t}^{i}

are the incoming tangential electric and magnetic fields on the boundary of the cylinder, respectively, ∊₀ = 8.85 10⁻¹² F/m, k₁ = n₁ω/c₀ and c₀ = 299792458 m/s. The BIEs are solved for e_z and h_t, the unknown total tangential electric and magnetic fields on the boundary, respectively. The Green’s function for a source at r′ = x′x̂ + y′ŷ and an observation point at r = xx̂ + yŷ is given by

G (r; r^{'}) = \frac{j}{4} H_{0}^{(2)} (k_{1} d (r, r^{'})),

and similarly for k₀ and G₀. d(r, r′) is the distance between the source and observation points, i.e.

\sqrt{{(x - x^{'})}^{2} + {(y - y^{'})}^{2}}

. The contour C denotes the boundary of the dielectric object, C⁻ and C⁺ being its interior and exterior boundary, respectively.

At the surface of a PEC body, only Eq. (1) and the unknown h_t at C⁺ remain. So, e_z and the entire right-hand side of the equation should be replaced by zero.

Similar expressions exist for the transverse electric (TE_z) problem and the set of equations can readily be extended to include several scatterers [1].

To solve the BIEs Eq. (1) and Eq. (2) with the well-known Galerkin MoM, pulse basis functions are used to expand h_t along the contour C and as test functions for Eq. (1). Hat functions are used to expand e_z and as test functions for Eq. (2) [12].

Fig. 1 Object illuminated by an electromagnetic wave.

Download Full Size | PDF

2.2. The coordinate stretching formalism

In order to damp the influence of subparts of the objects, the PML’s coordinate stretching formalism is employed, in a similar way as in FEM [9]. Ideally, a PML layer (a) exponentially damps impinging fields and (b) does not introduce reflections at the free space/PML interface. For simplicity, we focus on stretching in the x-direction. A coordinate transformation x → x̃ is introduced as follows (see also Fig. 2 for an illustration):

\tilde{x} = x_{r} + j x_{i} = \int_{0}^{x} d x^{'} \tilde{χ} (x^{'}),

with χ̃(x) a complex function of x. x_r and x_i denote the real and imaginary part of x̃, respectively. This is equivalent to the introduction of a new medium. In this new medium, the derivative in the x-direction is stretched in the following way:

\frac{\partial}{\partial x} \to \frac{\partial}{\partial \tilde{x}} = \frac{1}{\tilde{χ} (x)} \frac{\partial}{\partial x} .

It can be shown [9, 13] that by introducing these changes into Maxwell’s equations, the characteristic impedance of the new medium remains unchanged w.r.t. to the background medium, irrespective of χ̃(x). Moreover, the wave equation

(\frac{\partial^{2}}{\partial {\tilde{x}}^{2}} + k_{x}^{2}) \tilde{ϕ} (x) = 0,

is satisfied by

\tilde{ϕ} (x) = A exp (- j k_{x} \tilde{x}) + B exp (j k_{x} \tilde{x}),

which further indicates that no reflections occur at the interface between free space and the new medium, irrespective of frequency or angle of incidence, if x̃ is a continuous function of x.

Assuming that the interface between the free-space background medium and the new medium is situated at x = x₀, as shown in Fig. 2, x̃ becomes

\tilde{x} = {\begin{array}{l} x, & x \leq x_{0} \\ x_{r} + j x_{i} (x_{i} \neq 0), & x > x_{0} \end{array}

Introducing Eq. (8) into the plane wave solution propagating to the right in Eq. (7) leads to

A exp (- j k_{x} x_{r}) exp (k_{x} x_{i}),

such that the wave is exponentially damped if x_i is negative for x > x₀. Thus, the new medium acts as a PML.

Fig. 2 Illustration of the PML’s coordinate stretching formalism.

Download Full Size | PDF

2.3. PML layers in the BIE-MoM

After discretization of the objects in the MoM, i.e. segmentation of their boundaries, the begin-and endpoints of each segment are transformed according to the stretching formalism of Sec. 2.2. Consequently, a re-interpretation of geometrical quantities such as segment length, tangent vector, etc. imposes itself. First, the distance function appearing in Eq. (3) has to be adjusted. For complex vectors r̃ and r̃′, the “distance” function becomes

\begin{array}{l} d (\tilde{r}, {\tilde{r}}^{'}) & = \sqrt{{(\tilde{x} - {\tilde{x}}^{'})}^{2} + {(y - y^{'})}^{2}} \\ = [{(x_{r} - x_{r}^{'})}^{2} - {(x_{i} - x_{i}^{'})}^{2} + {(y - y^{'})}^{2} \\ {+ 2 j (x_{r} - x_{r}^{'}) (x_{i} - x_{i}^{'})]}^{1 / 2} . \end{array}

Note that, by the conventional definition of the square root function for complex arguments, the sign of the imaginary part of d(r̃, r̃′) and the sign of the imaginary part inside the square brackets of the right-hand side of Eq. (10) are the same. Thus, to obtain a negative imaginary part for d(r̃, r̃′), which is required given Eq. (9), x_i must be monotonically decreasing and x_r is logically chosen to be monotonically increasing as a function of x. Naturally, the proposed square root function does not represent a distance measure, as it is complex-valued. With a slight abuse of terminology, we will continue to use the term “distance function”, as previously done in literature [14].

The concept of tangential and normal vectors must also be generalized to complex coordinates. The concept of a line between two complex vectors p̃₀ and p̃₁ is interpreted as the collection of complex vectors p̃ fulfilling

\tilde{p} = (1 - s) {\tilde{p}}_{0} + s {\tilde{p}}_{1} \equiv {\tilde{p}}_{0} + s \tilde{l} \tilde{t},

with s ∈ [0, 1], l̃ = d(p̃₀, p̃₁) and t̃ = 1/l̃(p̃₁ − p̃₀). t̃ is the extension of the definition of a tangential vector along a line. The extension of the definition of a vector ñ normal to t̃ follows naturally: ñ = ±(ẑ × t̃).

Using Eq. (5), the normal derivatives appearing in Eq. (1) and Eq. (2) are now replaced in the following way:

\frac{\partial}{\partial n} = n \cdot \nabla \to \tilde{n} \cdot (\frac{1}{\tilde{χ} (x)} \frac{\partial}{\partial x} \hat{x} + \frac{\partial}{\partial y} \hat{y}) .

The coupling integrals appearing in classical Galerkin-MoM are computed after substitution of the pertinent quantities. Eq. (11) ensures that integration limits over the domain of test and basis functions are real, such that a traditional Gaussian quadrature rule may be employed. For singular coupling integrals, the analytical expressions given in [1] remain valid. Lastly, in a classical Galerkin-MoM, pulse and hat functions are normalized such that the area below the function measures unity. Here, the normalization factors are adjusted by a formal integration over complex limits, leading to identical expressions as traditional results, where segment lengths are replaced by their complex counterparts.

Equipped with these changes to the MoM solver, the novel MoM-PML layer is easily integrated into already existing simulation software. The influence of the part of the objects that resides inside the PML will be damped, thus creating the possibility to efficiently focus on only part of the structure.

2.4. Coordinate stretching functions

There are several types of stretching functions that fulfill the requirements posed by Eq. (9) and Eq. (10). We investigate three different types, inspired by previously published results [9,15,16].

Table 1 summarizes the χ̃(x) function and the resulting stretched x̃ for each type inside the PML layer (x > x₀). For x ≤ x₀, no stretching is applied (χ̃(x) = 1 and x̃ = x).

Table 1. Three different types of coordinate stretching for x > x₀. Here, $σ (x) = σ_{\max} {(\frac{x - x_{0}}{D})}^{m}$ .

View Table | View all tables in this article

3. Accuracy and PML optimization

For reasons of efficiency, the computational domain inside the PML is kept as small as possible. By making the attenuation parameters (type I: α, type II: θ, type III: σ_max) large, interactions inside the PML are damped faster. It is however known that discretization introduces reflections at the free space/PML interface. In [17], it is shown, for a FEM solver and employing a type I stretching, that a reflection at the free space/PML interface arises with a magnitude proportional to α²Δ², with Δ the discretization length. Consequently, the attenuation parameter cannot be chosen arbitrarily large.

An optimization test for the attenuation parameters is performed using the geometry shown in Fig. 3. A finite PEC wedge CBAED of opening angle 45° is truncated at x = 10 m. A PML layer is added at x = 7 m to damp the influence of the right part BAED and the two right tips A and E. A unit-amplitude, TM_z-polarized, plane wave with wavelength λ = 1 m, impinging at an angle of 100° with the x-axis, illuminates the wedge.

Fig. 3 Geometry used for the optimization tests of the PML layer.

Download Full Size | PDF

The dominant scattering mechanisms are the reflection at the upper face CA and the diffraction at the left tip C of the wedge. As a reference, we employ the analytical solution for the scattering at an infinite PEC wedge with external opening angle νπ, illuminated by a unit-amplitude, TM_z-polarized plane wave, coming in at angle ϕ₀ with the positive x-axis. The z-oriented current on the upper face of the wedge, at distance r from the tip, is given by [10]

j_{z} (r) = - \frac{4 j ∊_{0} c_{0}^{2}}{ν r ω} \sum_{l = 1}^{\infty} \frac{l}{ν} exp (j \frac{l π}{2 ν}) J_{l / ν} (k r) sin (\frac{l}{ν} ϕ_{0}),

and the currents at the lower face by

j_{z} (r) = - \frac{4 j ∊_{0} c_{0}^{2}}{ν r ω} \sum_{l = 1}^{\infty} {(- 1)}^{l + 1} \frac{l}{ν} exp (j \frac{l π}{2 ν}) J_{l / ν} (k r) sin (\frac{l}{ν} ϕ_{0}) .

We truncated this series at 200 terms.

Figures 4a–4c show the current amplitude on the two faces of the wedge for each stretching type and for different values of the attenuation parameter. The absciccae correspond to the pertinent positions shown in Fig. 3. For all three stretching types, increasing the value of the attenuation parameter leads to better damping of the current density near the undesired tips A and E at x = 10 m. For high values of the attenuation parameter of stretching types I and III, a peak in the current density is observed on the segments just inside the PML layer, close to the interface B.

Fig. 4 Results of the optimization tests. The points A–E correspond to the pertinent points in Fig. 3. For type III stretching, D = 4.

Download Full Size | PDF

To assess the fitness of the different stretching formalisms, we investigate in more detail the trade-off between these two apparently counteracting demands, i.e. the damping of the PML versus the accuracy of the solution at the free space/PML interface. In Fig. 5, the current at the tip A, which is a measure for the damping of the PML layer, is plotted against the relative error (RE) of the solution at B, just to the left of the free space/PML interface. Compared to the analytical reference result Eq. (13), this RE is defined as follows:

RE = | \frac{solution - reference}{reference} | .

The two competing demands plotted in Fig. 5 form Pareto fronts. Type I and type II stretching are both linear stretching formalisms, i.e. the imaginary part of the stretched coordinate increases linearly with increasing x. As such, they perform similarly over a certain range of their respective PML parameters. Damping by types I and II is negligible over the range of parameters that lead to an RE<10%. More damping can only be achieved at the cost of a worse RE. The nonlinear stretching of type III provides gradual damping, leading to the much better results, i.e. a high level of damping can be achieved while maintaining a low RE. Nonlinear stretching has been used extensively in the past due to its better trade-off of damping and accuracy [9,16,17].

Fig. 5 Pareto fronts for demands on accuracy and damping when adding a PML layer.

Download Full Size | PDF

In Fig. 6, the RE is shown for the z-component of the total electric field in the shadow zone of the wedge, obtained by using a type III stretching PML layer, with m = 3, σ_max = 10 and D = 4, when compared to the analytical solution [10]. These parameters correspond to the stretching that provides the best combination of high damping and low RE at the free space/PML interface in Fig. 5. In the shadow zone, where the only field contribution is due to the diffraction at the edge C under study, the RE remains below 3%, even very close to the diffracting tip and also near the free space/PML interface. In the other zones, the dominant contributions stem from the incident and/or reflected waves, leading to even better accuracy.

Fig. 6 RE on the z-oriented electric field compared to the analytical solution.

Download Full Size | PDF

4. Polarization beam splitter

We now investigate the PBS presented in Fig. 7. The design is based on a prototype consisting of SOI nanowires reported in [11]. Signals of both TM_x and TE_x polarization and wavelength λ = 1550 nm are coupled into the PBS at the left boundary of the upper slab. Note that in this example, in contrast to the previous sections, TM and TE refer to the axis of beam propagation instead of the axis of invariance, to conform with the convention most often used in optical waveguide theory. The widths w₁ (narrow upper and lower slabs) and w₂ (wide middle slab) are designed so that the phase matching between TM_x modes of the narrow slabs and the wider slab is optimally satisfied and this is not the case for TE_x modes. This allows transfering TM_x polarized beams from the upper slab to the lower slab, while ensuring that TE_x polarized beams remain confined in the upper slab. The refractive index of the slabs in Fig. 7 is n₁ = 3.45 and they are surrounded by free space.

Fig. 7 Polarization beam splitter with three silicon slabs.

Download Full Size | PDF

Modeling the fields inside the PBS is not possible with the traditional MoM, as back-reflections occur at the extraneous backends of the slabs and this even independently of the truncation size. With our novel, advocated method, however, we can damp the effect of these backends. Thereto, the structure is terminated by a PML (see Fig. 7) with the same parameters as in the previous example.

In order to choose an appropriate pair of slab widths w₁ and w₂, first the effective index n_eff of the propagating modes in a dielectric slab waveguide is computed as a function of the slab’s width. The result is shown in Fig. 8. Then, for a fixed w₁, w₂ (>w₁) is chosen such that n_eff of the TM_x,1 mode in the wide slab corresponds to n_eff of the TM_x,0 mode of the narrow slab (lower horizontal line in Fig. 8). Moreover, for the same pair of widths, n_eff of the TE_x,0 mode of the narrow slab does not correspond to n_eff of any of the modes of the wide slab. Therefore, the PBS shown in Fig. 7 is expected to split the propagating TE_x,0 and TM_x,0 modes in the upper slab.

Fig. 8 Effective indices of guided modes in a silicon slab waveguide.

Download Full Size | PDF

In Fig. 9, the fields are shown when the upper slab of the PBS is excited at its left side by a Gaussian beam with width equal to w₁. For TM_x excitation in Fig. 9(a), the coupling of the beam from the upper slab to the lower slab is clearly visible. For TE_x excitation, shown in Fig. 9(b), the beam remains confined in the upper slab.

Fig. 9 Field density in a polarization beam splitter with three optical waveguides (the colors are an indication for the field strength).

Download Full Size | PDF

The requirements in terms of memory and computational time for this simulation are summarized in Table 2. The simulation was run on a dual Intel Core CPU E6750 at 2.66 GHz with 8 GB of RAM.

Table 2. Computational resources for the simulation of the silicon PBS of Fig. 7.

View Table | View all tables in this article

5. Conclusion

We have integrated the coordinate stretching of PMLs, as used in FEM, into an MoM solver. This allows the full-wave study of parts of optical waveguides, e.g. the coupling region, while damping the interaction with the remaining parts of the structure. The necessary adjustments to already existing MoM solvers are minimal. Using a traditional wedge diffraction example, the technique proves to be very accurate compared to the analytical result. In particular, nonlinear stretching provides the best results in terms of accuracy. Finally, the method was applied to the coupling of a Gaussian beam into a silicon polarization beam splitter.

Acknowledgments

This research was kindly funded by the Research Foundation Flanders (FWO).

References and links

1. J. Fostier and F. Olyslager, “A GRID computer implementation of the multilevel fast multipole algorithm for full-wave analysis of optical devices,” IEICE Trans. Commun. E90B, 2430–2438 (2007). [CrossRef]

2. M. Stallein, “Coupling efficiency of Gaussian beams into step-index waveguides–an improved ray-optical approach,” J. Lightwave Technol. 26, 2937–2945 (2008). [CrossRef]

3. M. Stallein, “Coupling of a Gaussian beam into a planar slab waveguide using the mode matching method,” in “Progress in Electromagnetic Research Symposium (PIERS),” (2004), pp. 309–312.

4. G. F. Herrmann, “Numerical computation of diffraction coefficients,” IEEE Trans. Antennas Propag. 35, 53–61 (1987). [CrossRef]

5. G. Apaydin and L. Sevgi, “A novel wedge diffraction modeling using method of moments (MoM),” ACES J. 30, 1053–1058 (2015).

6. H. Derudder, F. Olyslager, and D. De Zutter, “Efficient calculation of far-field patterns of waveguide discontinuities using perfectly matched layers,” Electron. Lett. 36, 798–799 (2000). [CrossRef]

7. W. Bogaerts, P. Bienstman, D. Taillaert, R. Baets, and D. De Zutter, “Out-of-plane scattering in photonic crystal slabs,” IEEE Photon. Technol. Lett. 13, 565–567 (2001). [CrossRef]

8. W. Bogaerts, R. Baets, P. Dumon, V. Wiaux, S. Beckx, D. Taillaert, B. Luyssaert, J. Van Campenhout, P. Bienstman, and D. Van Thourhout, “Nanophotonic waveguides in silicon-on-insulator fabricated with CMOS technology,” J. Lightwave Technol. 23, 401–412 (2005). [CrossRef]

9. W. C. Chew, J. M. Jin, and E. Michielssen, “Complex coordinate stretching as a generalized absorbing boundary condition,” Microw. Opt. Technol. Lett. 15, 363–369 (1997). [CrossRef]

10. P. Y. Ufimtsev, “Fast convergent integrals for nonuniform currents on wedge faces,” Electromagn. 18, 289–313 (1998). [CrossRef]

11. D. Dai, “Silicon polarization beam splitter based on an asymmetrical evanescent coupling system with three optical waveguides,” J. Lightwave Technol. 30, 3281–3287 (2012). [CrossRef]

12. F. Olyslager, D. De Zutter, and K. Blomme, “Rigorous analysis of the propagation characteristics of general lossless and lossy multiconductor transmission lines in multilayered media,” IEEE Trans. Microw. Theory Techn. 41, 79–88 (1993). [CrossRef]

13. L. Knockaert and D. De Zutter, “On the stretching of Maxwell’s equations in general orthogonal coordinate systems and the perfectly matched layer,” Microw. Opt. Technol. Lett. 24, 31–34 (2000). [CrossRef]

14. R. Torretti, Philosophy of Geometry from Riemann to Poincaré (Springer, 1978). [CrossRef]

15. D. Pissoort, D. Vande Ginste, and F. Olyslager, “Including PML-based absorbing boundary conditions in the MLFMA,” IEEE Antennas Wireless Propagat. Lett. 4, 312–315 (2005). [CrossRef]

16. S. D. Gedney, “An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices,” IEEE Trans. Antennas Propag. 44, 1630–1639 (1996). [CrossRef]

17. F. Collino and P. B. Monk, “Optimizing the perfectly matched layer,” Comput. Methods in Appl. Mech. Eng. 164, 157–171 (1998). [CrossRef]

	type I	type II	type III
χ̃(x)	1 − jα	e^−jθ	1 − j(1 + m)σ(x)
x_r	x	x₀ + cos θ(x − x₀)	x
x_i	−α(x − x₀)	− sin θ(x − x₀)	−σ(x)(x − x₀)
	α > 0	θ ∈ ]0, π/2[	m ∈ ℕ σ_max > 0 D > sup(x − x₀)

	TM_x excitation	TE_x excitation

unknowns	12452	12452
memory	2.3 GB	2.3 GB
iterations	27	34
solution time	58 s	1 min 13 s

	type I	type II	type III
χ̃(x)	1 − jα	e^−jθ	1 − j(1 + m)σ(x)
x_r	x	x₀ + cos θ(x − x₀)	x
x_i	−α(x − x₀)	− sin θ(x − x₀)	−σ(x)(x − x₀)
	α > 0	θ ∈ ]0, π/2[	m ∈ ℕ σ_max > 0 D > sup(x − x₀)

	TM_x excitation	TE_x excitation

unknowns	12452	12452
memory	2.3 GB	2.3 GB
iterations	27	34
solution time	58 s	1 min 13 s

Full-wave simulation of optical waveguides via truncation in the method of moments using PML absorbing boundary conditions

Abstract

1. Introduction

2. Formalism

2.1. BIE-MoM

2.2. The coordinate stretching formalism

2.3. PML layers in the BIE-MoM

2.4. Coordinate stretching functions

3. Accuracy and PML optimization

4. Polarization beam splitter

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (9)

Tables (2)

Equations (15)

Optics Express