{'fraction_non_alphanumeric': 0.03604717201750576, 'fraction_numerical': 0.018107105869775833, 'mean_word_length': 4.4864369501466275, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Second Kind Integral Equation Formulation for the Mode Calculation of Optical Waveguides

We present a second kind integral equation (SKIE) formulation for calculating the electromagnetic modes of optical waveguides, where the unknowns are only on material interfaces. The resulting numerical algorithm can handle optical waveguides with a large number of inclusions of arbitrary irregular cross section. It is capable of finding the bound, leaky, and complex modes for optical fibers and waveguides including photonic crystal fibers (PCF), dielectric fibers and waveguides. Most importantly, the formulation is well conditioned even in the case of nonsmooth geometries. Our method is highly accurate and thus can be used to calculate the propagation loss of the electromagnetic modes accurately, which provides the photonics industry a reliable tool for the design of more compact and efficient photonic devices. We illustrate and validate the performance of our method through extensive numerical studies and by comparison with semi-analytical results and previously published results.


Introduction
Optical fibers and waveguides are important building blocks of many photonic devices and systems in telecommunication, data transfer and processing, and optical computing. Indeed, most photonic devices consist of approximately straight waveguides as input and output channels with complicated functional structures between the two. Two main mechanisms by which the electromagnetic wave can be confined in optical fibers or waveguides are total internal reflection and photonic band gap guidance [42,23]. Generally speaking, when the refractive index of the core is greater than that of the surrounding material, the light is confined in the core by total internal reflection; when the (hollow) core has a smaller refractive index, confinement can be achieved through photonic band gap guidance. In both cases, the propagating electromagnetic modes of optical fibers and waveguides depend on physical parameters such as the input light wavelength, refractive indices, and the geometry of the cross section of fibers and waveguides. To reduce the cost of designing new photonic devices, accurate and efficient simulation tools are in high demand in integrated photonics industry. The first step in the photonics simulation is to compute a complete set of propagating modes accurately and efficiently for optical fibers or waveguides.
There has been extensive research on the mode calculation of optical fibers and waveguides and various numerical methods have been developed. These include the effective index method [38,10], the plane wave expansion method [17,48,24], the multipole expansion method [53,54,32,12,11,52,30], finite difference methods [19,50,14,56], finite element methods [4,16,26,5,45,41,44], boundary integral methods [51,18,33,34,13,15,3,35,36,49,43], etc. Here we do not intend to review these methods in great detail, but note that the effective index method is generally of low order making it difficult to calculate the propagation constant to high accuracy; the plane wave expansion method implies an infinite periodic medium; the multipole expansion method requires that each core be of circular shape and that the cores be well separated from each other; finite difference and finite element methods requires a volume discretization of the cross section in a truncated computational domain with some artificial boundary conditions or perfectly matched layers imposed on or near the boundary of the truncated domain. When optical fibers and waveguides consist of many cores of arbitrary shape, these methods often need excessively large amount of computing resource in order to accurately calculate the imaginary part of the propagation constant, which is related to the propagation loss of the electromagnetic modes and thus of fundamental importance for the design purpose.
On the other hand, boundary integral methods represent the electromagnetic fields via layer potentials which satisfy the underlying partial differential equations automatically. One then derives a set of integral equations through the matching of boundary conditions with the unknowns only on the material interfaces. Thus the dimension of the problem is reduced by one and complex geometries can be handled relatively easily. Among the aforementioned work on boundary integral methods, [13,3,35,43] present numerical examples with high accuracy for smooth geometries. In [3] and [43], the field components E z and H z (with zaxis the longitudinal direction of the waveguide) are represented via four distinct single layer potentials and the resulting system is a mixture of first kind and singular integral equations; both authors apply the circular case as a preconditioner to obtain a well conditioned system for smooth boundaries. In [13], E z and H z are represented via a proper linear combination of single and double layer potentials in such a way that the hypersingular terms are cancelled out. The resulting system still contains the tangential derivatives of the unknown densities and layer potentials and thus is not of the second kind. In [35], Dirichlet-to-Neumann (DtN) maps for H x and H y are used to construct a system of two integral equations, where each DtN map is in turn computed by a boundary integral equation with a hypersingular integral operator and a method in [29] is applied to evaluate the DtN map to high accuracy for smooth boundaries. While these methods are all capable of computing the propagation constant to high accuracy for smooth cases, it is not straightforward to extend them to treat nonsmooth cases such as standard dielectric rectangular waveguides in integrated optics.

Remark 1
We would like to remark here that [3] has a subsection titled "Buried Rectangular Dielectric Waveguide". In that subsection, the authors approximate the rectangular waveguide via a smooth super-ellipse and compute the propagation constant for the super-ellipse. Though Fig. 2 in [3] achieves about 9 digit accuracy for the super-ellipse, Fig. 3 in [3] shows only about two digit accuracy for the propagation constant of the genuine rectangular waveguide which is regarded as a limit of the super-ellipse.
In this paper, we construct a system of SKIEs formulation for the mode calculation of optical waveguides. Our starting point is the dual Müller's formulation [39] for the time-harmonic Maxwell's equations in three dimensions. We then reduce the dimension of the integration domain by one using the key assumption of the mode calculation that the dependence on z of the electromagnetic fields and the unknown densities is of the form e iβz , where β is the propagation constant of the mode. Accordingly, the layer potentials defined on the surface of the waveguide are reduced to the layer potentials defined only on the boundary of the cross section. And the boundary conditions lead to a system of four SKIEs where each integral equation consists of a sum of a constant multiple of the identity operator and several compact operators. All involved integral operators have logarithmically singular kernels and thus are straightforward to discretize with high order quadratures. Hence, our formulation leads to a numerical algorithm which can handle the mode calculation of optical waveguides of arbitrary geometries (smooth and nonsmooth) accurately and efficiently.
The paper is organized as follows. In Section 2, we develop the SKIE formulation. We discuss briefly the discretization scheme and numerical algorithm for the mode finding in Section 3. In Section 4, we illustrate the performance of our scheme numerically and compare our results with semi-analytical results or previously published results in [13,35,43] for PCFs with smooth holes. High accuracy results for dielectric rectangular waveguides are presented in Section 5. And concluding remarks are contained in Section 6.

SKIE formulation for the electromagnetic mode calculation
In this section, we first derive the SKIE formulation for the electromagnetic mode calculation when the waveguide consists of only one core. The extension to the case of multiple cores or holes is straightforward.

Notation
We follow standard conventions in the mode calculation literature and assume that electromagnetic fields are propagated along the z-axis, and that the geometric structure of the optical waveguides is completely determined by its cross section in the xy-plane (see Fig. 1 for an illustration). The waveguide consists of several cores (or holes) denoted by Ω 1 , · · · , Ω N with refractive indices n 1 , · · · , n N , respectively. The surroundings of these cores or the cladding is denoted by Ω 0 with refractive index n 0 . The boundary of each core is denoted by Γ i with ν the unit outward normal vector and τ the unit tangential vector, respectively. We denote points in R 2 by P and Q.

PDE formulation
The source free Maxwell equations in each homogeneous region are given by: where n is the index of refraction of the region. Assuming the time dependence of the electromagnetic field is e −iωt , i.e., E(x, y, , t) = E(x, y, z)e −iωt and H(x, y, , t) = H(x, y, z)e −iωt , and rescaling the magnetic field H by µ 0 ǫ 0 , equations (1)-(2) become where k = k v n is the wave number in the region, k v = ω/c is the wave number in vacuum, and c = 1 √ ǫ 0 µ 0 is the speed of light in vacuum. In the mode calculation, a further important assumption is that the electromagnetic field takes the following form: H(x, y, z) = H(x, y)e iβz .
We observe that with this assumption each component of the electromagnetic fields of the mode satisfies the Helmholtz equation in two dimensions: On the material interface, the boundary conditions are that the tangential components of the electromagnetic fields are continuous. That is, where [·] denotes the jump across the material interface. It is clear that z and τ are two tangential directions for the waveguide geometry. Thus the boundary conditions can be written explicitly as follows: Finally, we would like to remark that combining the assumptions (7), (8) and Maxwell's equations, it is straightforward to verify that there are only two independent components for the electromagnetic fields in each region. For instance, the components E x , E y , H x , H y are completely determined by E z and H z via the following relations: In fact, [13,3,43] start from the integral representation of E z and H z to develop the integral equation formulation for the mode calculation.

SKIE formulation
We first provide an informal description about our construction of the SKIE formulation for the mode calculation of optical waveguides. Previous integral equation formulations in [13,3,43,35] start from two scalar variables, i.e., two components of the electromagnetic fields (E z and H z in [13,3,43] and H x , H y in [35]) and then set up the integral equations through the boundary conditions. Here we start from the dual Müller's representation in [39] for the time harmonic electromagnetic fields which leads to an SKIE formulation for the dielectric interface problems in three dimensions. As an integral representation for three dimensional problem, the integration domain is over the boundary ∂Ω and the unknowns are the surface currents J and M on ∂Ω. For the mode calculation of an optical waveguide, the boundary is an infinitely long cylinder, i.e., ∂Ω = Γ × (−∞, ∞) where Γ is the boundary of the cross section of the waveguide. Since all electromagnetic field components have e iβz dependence on z, it is natural to assume that the unknown surface currents J and M have the same dependence on z as well.
Combining these two factors, we are able to reduce the dimension of the representation by one and derive an integral representation of E(x, y) and H(x, y) using layer potentials defined only on Γ. It is readily to verify that the boundary conditions lead to a system of SKIEs for the mode calculation, with the propagation constant β appearing as a nonlinear parameter of the system.

Reduction of layer potentials in the mode calculation
The dual Müller's representation [39] assumes that in each region the time harmonic electromagnetic field E(x, y, z) and H(x, y, z) have the following representation: Here ∂Ω is the boundary of the three dimensional domain, J and M are the unknown surface electric and magnetic currents, and is the single layer potential for the Helmholtz equation in 3D with similar expression for S k For the waveguide geometry, ∂Ω = Γ × (−∞, ∞). Furthermore, since both E(x, y, z) and H(x, y, z) depend on z in the form of e iβz , it is natural to assume that the surface currents J(x, y, z) and M(x, y, z) have the same z dependence as well, that is, We now introduce some notation to be used subsequently. We denote k β = k 2 − β 2 and the Green's function for the Helmholtz equation (9) in 2D by G(P, Q) 0 is the zeroth order Hankel function of the first kind [1]. We further denote by S the single layer potential operator defined by the formula D the double layer potential operator defined by the formula and T the anti-double layer potential operator (for the lack of standard terminology) defined by the formula The following lemma shows that the Fourier transform of the 3D Helmholtz Green's function is the 2D Helmholtz Green's function.
where P = (x, y) and Q = (x ′ , y ′ ) are the projections of r = (x, y, z) and r ′ = (x ′ , y ′ , z ′ ) onto the xy-plane, respectively, and the wavenumber k β in G is chosen to have nonnegative imaginary part.
We now note that the well-known Sommerfeld integral identity (for y > 0) [47] is given by: (21) follows from the substitution λ = it, y = k, x = iβ and k s = a and appropriate contour deformation.
The following corollary reduces the single layer potential S k ∂Ω in 3D to a single layer potential in 2D in the mode calculation of an optical waveguide.
where the second equality follows from (21).

Integral representations of the electromagnetic field in the mode calculation
Let (î,ĵ,k) be the basis of the Cartesian coordinate withk pointing along the positive z direction. We write the unit tangent vector in its component form τ = τ 1î + τ 2ĵ . Then the outward unit normal vector ν is given by ν = ν 1î + ν 2ĵ = τ 2î − τ 1ĵ . Since J and M are unknown surface currents and τ ,k are two locally orthonormal tangential vectors at a point (x, y, z) ∈ ∂Ω, we may write J(x, y) and M(x, y) in (17) for (x, y) ∈ Γ as follows: We now write down our integral representations for the electromagnetic field in each region in the mode calculation of optical waveguides.

Theorem 2.3 (Integral representations of the electromagnetic field)
Suppose that the unknown surface currents J(x, y, z) and M(x, y, z) are given in the form of (17) and (26).
where E x , E y and E z are defined via the formulas: where H x , H y and H z are defined via the formulas: Proof By Corollary 2.2, we have Thus, and similar expression holds for S k ∂Ω [M]. Since and we have Using the identity and Thus, Substituting (30), (33) (with J replaced by M), (34), (36) into (14), we obtain (27). And (28) can be obtained similarly.
Remark 2 Though β is assumed to be real in Lemme 2.1, we will use (27)-(28) to represent the electromagnetic field even when β is complex for the mode calculation. In fact, it is straightforward (though a little tedious) to verify (27)-(28) satisfy the relations (12)-(13) even when β is complex. The hard part is to verify the following identities: and other similar identities. In other words, we might start directly from the representation (27)-(28) for the mode calculation without even mentioning the dual Müller's representation, which only serves as a formal derivation tool.

Main theoretical result
We are now in a position to derive the SKIE formulation for the mode calculation of photonic waveguides. Since the boundary conditions for dielectric problems are that the tangential components of E and H must be continuous across the boundary. We first combine the first two equations in (27) and (28) to obtain E τ and H τ . Together with E z and H z , we list all four tangential components as follows: We now scale and nondimensionalize the quantities in the above equation by multiplying all length quantities with k v . We havẽ where n e = β/k v is the effective index.
For each region Ω i (i = 0, 1 for now), we denote the single, double, and anti-double layer potential operators in that region by S i , D i , and T i , respectively. We now summarize our main theoretical result in the following theorem.
Here the nonzero entries of A are given by the following formulas: The effective index n e = β/k v of the electromagnetic mode is a complex number for which the above linear system has a nontrivial solution.
Proof We first obtain the nondimensionalized tangential components in Ω i (i = 0, 1) by replacing n, S, D, T in (39) with n i , S i , D i , T i , respectively. We then substitute the resulting expressions into the boundary conditions (11) and observe that the principal part of the linear system is exactly the matrix A defined in (41)- (42). It is easy to see now that only D and S ν have nonzero jumps when the target point P approaches the boundary and their jump relations lead to the diagonal matrix D in (41) (see, for example, [27]). Furthermore, all entries A ij are compact operators due to following reasons. First, S, and the principal part of D and S ν are compact from standard potential theory [27]. Second, T 0 − T 1 and S 0,τ − S 1,τ are compact since their kernels are only logarithmically singular due to the cancellation of more singular terms. Third, T 0,τ − T 1,τ is also compact since the hypersingular terms in difference kernel cancel out. Thus, the resulting system is of the second kind.

Extension to the waveguide with multiple holes
When the waveguide consists of multiple holes, say, Ω 1 , · · · , Ω N , we simply represent E(x, y) and H(x, y) in each region Ω i (i = 1, · · · , N ) via (27)-(28) with k, S, D, T , Γ replaced by k i , S i , D i , T i , Γ i , respectively. For the exterior domain Ω 0 , we have similar representations for E(x, y) and H(x, y), except that the boundary Γ is replaced by Γ 0 = Γ 1 ∪ · · · ∪ Γ N . The boundary conditions on each Γ i lead to a 4N × 4N block system. We note that all 4 × 4 diagonal blocks are similar to (40)- (42) and that the entries in off-diagonal blocks are all compact operators since their kernels are smooth due to the fact that the target point P is bounded away from the source curve. Thus, the system is of the second kind.

Discretization and numerical algorithms
As pointed out in Theorem 2.4, all the integral operators in (40)- (42) have kernels with logarithmic singularity. When the boundary curves are smooth, there are many high-order quadrature rules based on local modifications of the trapezoidal rule [2,25] or kernel splitting method [46,28]. When the boundary curves contain corners, recent developments in [21,20,22,6,8,9] treat such cases with high order quadratures very efficiently. In order to achieve optimal efficiency in discretization, one should treat smooth and nonsmooth cases separately and implement the aforementioned schemes (say, [2,9]) accordingly.
Here we adopt a simpler scheme to discretize both smooth and nonsmooth cases. We divide the boundary curve into N c chunks as follows. If the curve is smooth, the chunks are of equal length in the parameter space. If the curve contains corners, we will have finer and finer dyadic chunks toward the corner. On each chunk, the unknown densities are approximated via a p − 1th order polynomial and the set of collocation points are the collection of the images of p shifted and scaled Gauss-Legendre nodes on each chunk. When the target points and the source points lie on the same chunk, we apply a precomputed generalized Gaussian quadrature (see, for example, [7,37,55]) to discretize the associated (logarithmically) singular integrals; and when the target points lie outside the integration chunk, we simply use an adaptive Gaussian integrator to compute the corresponding matrix entries. The total number of discretization points is N = N c p, where N c is the number of chunks, and the size of the resulting matrix M (n e ) is 4N × 4N .
We apply the same method as in [13] to find the effective index n e such that M (n e ) has a nontrivial nullspace. That is, we use Müller's method [40] to find the root of the function where u and v are two fixed random column vector of length 4N .
4 Numerical examples for smooth boundaries 4.1 Example 1: optical fiber with a circular core For our first example, we consider the silica optical fiber consisting of a single circular core of radius 50µm. At the incident wavelength λ = 1500nm, the refractive index of the core is n 1 = 1.4475 while that of the cladding is n 0 = 1.444. The effective index can be calculated independently by the multipole method [31] and our SKIE formulation. For this simple case, the multipole method is semi-analytical and one only needs to solve a 4 × 4 system to find the effective index every time. Thus its results can be taken as the reference values for comparison purpose. Table 1 presents the first five modes obtained by our formulation with 50 discretization points of the boundary and by the multipole method, which shows IEEE double precision agreement of these two results. We would like to remark here that our algorithm already achieves 13 digit accuracy with only 20 discretization points.

Example 2: PCFs with a hexagon ring and its perturbation
In this example, we consider two PCFs studied in [13] and [35]. The first PCF consists of a hexagon ring of circular air holes as shown in Figure 2(a). The six holes are equally spaced along the hexagon with hole diameter d = 5µm and hole pitch Λ = 6.75µm. The refractive index of the glass matrix is assumed to be n 0 = 1. 45   wavelength is 1450nm. The second PCF shown in Fig. 2(b) is the perturbation of the first one, where the holes have the parametrization: Here r(θ) is the boundary of the i-th perturbed circle centered at c i and the perturbation level h goes from 1% to 6%. We discretize each boundary curve by 100 points. Table 2 lists the results obtained by our algorithm and the corresponding results in [13] for the first PCF. We observe that our algorithm is able to recover the first 10 digits in [13] for all modes. Table 3 lists the results for the perturbed PCF. Compared with Table 2 in [13], we observe that the two results again have at least 10 digit agreement with each other. We have also studied the convergence rates for these two PCFs. The results are presented in Fig. 3, which shows the rapid convergence of our method.   Fig. 2(a).  Table 3: Effective index of the PCF shown in Fig. 2(b).

Example 3: PCFs with multiple-layer hexagon rings
We now consider PCFs with multiple layers. The examples are taken from [35,43]. The first PCF has five layers of hexagon rings surrounding a circular core as shown in Figure 4(a). The small air holes are equally spaced with the hole pitch Λ = 2.74µm and the diameter d = 0.95Λ. The diameter of the core is l = 2.5d. The refractive indices of the glass surroundings and the air hole are 1.45 and 1, respectively. The wavelength of the incident field is 1510nm. The second PCF has an elliptic core with minor axis 2a = 2.3µm and major axis 2b = 4.6µm surrounded by a three layers of hexagon rings. The hole pitch of the small air holes is Λ = 2µm and the diameter of each hole is d = 0.9Λ. The incident wavelength for the second PCF is λ = 1420nm. These two PCFs are shown in Fig. 4.
The results of our computation for these two PCFs are shown in Table 4. For comparison, we also list the corresponding results in [35] in columns 4 and 5. The results in [35] are seen to be very accurate and there is a 13 digit agreement between our results and the results in [35]. For the first PCF shown in Fig. 4(a), we discretize each small hole with 100 points and the center hole with 200 points; while for the second PCF shown in Fig. 4(b), we discretize each circular hole with 120 points and the center ellipse with 240 points.

Rectangular waveguides
In this section, we present a detailed numerical study on a high refractive index contrast silica waveguide. The cross section of the waveguide is of the square shape with the side length equal to 3.4µm. The refractive index of the cladding is n 0 = 1.4447, while that of the core is 2% higher, i.e., n 1 = 1.4447 × 1.02. The wavelength of the incident field is 1550nm.  Table 4: Effective index of the PCFs shown in Fig. 4. The second and third columns list the values obtained by our method and the last two columns list the values in [35].

Accuracy and conditioning of the SKIE formulation
We first check the accuracy and conditioning of the SKIE formulation. For this, we will solve the linear system M (n e )x = b with a nonzero right hand side for some n e ∈ [n 0 , n 1 ]. We set n e = 1.451 for testing purpose and we construct artificial electromagnetic field by placing a point source outside the square for the interior field and a point source inside the square for the exterior field. We then obtain the right hand side vector b by simply computing the difference of the tangential components of these artificial electromagnetic fields. After we have solved the linear system, we may use the representation (27)- (28) to evaluate the electromagnetic field both inside and outside and compare the numerical result with the exact solution (which is the artificial electromagnetic field by construction). We use GMRES to solve the linear system and GMRES terminates when the relative residual is less than 10 −14 . In Table 5  the number of discretization points on each side of the square, while the second row lists the number of iterations needed in GMRES. We observe that the number of iterations in GMRES is almost independent of the size of the linear system, which is a characteristics of the SKIE formulation. x-axis shows the number of discretization points on each side of the square; y-axis shows the relative error (in logarithmic scale) of the computed electromagnetic field against the exact value. Figure 5 shows the relative error with various number of discretization points for each side of the square. We observe that the order of convergence is about 10, which is in agreement with the theoretical value since the number of collocation points on each chunk is set to p = 10.

Computation of the effective index
It is known that this waveguide admits a single propagation constant (or effective index) with double degeneracy. The method in [38] produces an effective index approximately equal to 1.458 · · · . For comparison purpose, we have also implemented the formulation in [3,43] Table 6: Effective index of the rectangular waveguide. The first column lists the number of discretization points on each side of the square; the second and third columns list the real and imaginary parts of the effective index found via the formulation in [3,43]; the fourth and fifth columns list the real and imaginary parts of the effective index found via our SKIE formulation. Table 6 shows the effective index n e found by the formulation in [3,43] (denoted by SREP in the table) and by our SKIE formulation for various number of discretization points. We observe that while the SKIE formulation exhibits a consistent convergence behavior to about 13 digit accuracy as N increases, the single layer representation behaves much more erratically and achieves about only 4 digit accuracy. This is because that the resulting matrix from the discretization of the single layer representation becomes more and more ill-conditioned as N increases even when n e is sufficiently far away from the root of the function defined in (43).
N 150 300 450 600 κ SREP 4.8E+5 3.4E+7 2.0E+9 1.0E+10 Table 7: Condition number of the formulation in [3,43] for the rectangular waveguide. The first row lists the number of discretization points on each side of the square, while the second row lists the condition number of the resulting matrix. n e is set to 1.451 as in Table 5. Table 7 lists the condition numbers of the matrix M SREP for various number of discretization points, which shows that the condition number of the matrix increases very rapidly even when n e is sufficiently far away from the root of the function in (43).
We have also carried out the singular value decomposition (SVD) for the matrices M SREP and M SKIE at the effective indices listed in Table 6 for N = 300. The singular values are plotted out in Figure 6. We observe that the SKIE formulation has much cleaner singular value distribution. Indeed, the smallest two singular values are both about 2. while the third smallest singular value is about 1.6 × 10 −9 . To sum up, the SKIE formulation enables us to find the effective index accurately and robustly; while non-SKIE formulation will either give low accuracy or lead to spurious propagation mode due to ill-conditioning.

Conclusions
We have constructed a second kind integral equation formulation for the mode calculation of optical waveguides or fibers. The resulting numerical algorithm is capable of finding the propagation modes of optical waveguides or fibers with an arbitrary number of cores or holes of arbitrary shape (smooth or with corners). The algorithm is high order accurate so that it is capable of computing the propagation constant (including its imaginary part which is related to the propagation loss of the mode) with high fidelity. The algorithm is robust and well-conditioned due to its SKIE formulation. The algorithm is efficient since one only needs to discretize the material interfaces. This enables practitioners in the integrated photonics industry to have a reliable simulation tool for designing more compact and efficient optical components or devices.