A Quantitative Study of the Effect of Cladding Thickness on Modal Confinement Loss in Photonic Waveguides

There has been increasing interest in making the photonic devices more and more compact in the integrated photonics industry, and one of the important questions for manufacturers and design engineers is how to quantify the effect of the finite cladding thickness on the modal confinement loss of photonic waveguides. This requires at least six to seven digits accuracy for the computation of propagation constant $\beta$ since the modal confinement loss is proportional to the imaginary part of $\beta$ that is six to seven orders of magnitude smaller than its real part by the industrial standard. In this paper, we present an accurate and efficient method to compute the propagation constant of electromagnetic modes of photonic waveguides with arbitrary number of (nonsmooth) inclusions in a layered media. The method combines a well-conditioned boundary integral equation formulation for photonic waveguides which requires the discretization of the material interface only, and efficient Sommerfeld integral representations to treat the effect of the layered medium. Our scheme is capable of calculating the propagation loss of the electromagnetic modes with high fidelity, even for waveguides with corners imbedded in a cladding material of finite thickness. The numerical results, with more than $10$-digit accuracy, show quantitatively that the modal confinement loss of the rectangular waveguide increases exponentially fast as the cladding thickness decreases.


Introduction
In many photonic devices, the input and output channels take the form of (approximately) straight waveguides with uniform cross sections. It is well known that such straight waveguide structure support only a finite number of electromagnetic modes which can propagate with little * Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, New Jersey 07102 (shidong.jiang@njit.edu).
† Courant Institute of Mathematical Sciences, New York University, New York, New York, 10012 (lai@cims.nyu.edu). or no energy loss for a given frequency. This digitizes the design of photonic devices and reduces a seemingly infinite dimensional problem to a finite dimensional one, which greatly simplifies the design process. The problem of mode calculation is concerned with characterizing the nature of the propagating waves for a waveguide of given cross-section, including their propagation constant and the structure of the associated electromagnetic fields.
Recently, as the integrated optical industry has grown, engineers have been trying to assemble more and more photonic components into a single chip. As a result, the effect of finite cladding thickness on the performance of the individual photonic components requires more care and more accurate modeling. For photonic waveguides, some obvious but important question are as follows: how does the thickness of the cladding alter/affect the propagation constant of each mode? What is the minimal cladding thickness for a given propagation loss threshold? These problems are naturally cast as nonlinear eigenvalue problems in a layered medium and present new challenges for numerical simulations, as they involves boundary conditions along the infinitely long material interface separating different layers.
Popular methods for mode calculation such as finite difference [14,3] or finite element methods [13,10] require the discretization of a finite domain, supplemented by artificial boundary conditions, such as perfectly matched layers, to simulate the effect of an infnite medium. Using these methods to study the effect of cladding thickness on propagation loss is problematic, since the computational domain has to be much larger than the waveguide cross-section in order to take the effect of layers into account. Existing boundary integral methods [1,2,4,8,12] are capable of calculating the propagation constant to high accuracy when the boundary of the waveguide consists of smooth curves in the absence of layers, but they become ill-conditioned and inaccurate for waveguides with nonsmooth geometry such as the rectangular waveguides considered in this paper and less efficient when the effect of layers has to be included. At the same time, the modal confinement loss is connected with the imaginary part of the propagation constant of the electromagnetic modes, which is very often at least six to seven orders smaller than the real part of the propagation constant. Thus, one needs at least six to seven digits accuracy in the overall computation in order to obtain a single significant digit for the modal confinement loss. This high accuracy demand presents a great challenge, requiring that the scheme be well-conditioned, high-order, and efficient, so as not to consume excessively large amounts of computational resources and for it to be of practical use in design.
In this paper, we present an integral formulation for the electromagnetic mode calculation of photonic waveguides in layered media. The formulation combines a carefully chosend boundary integral representation [6] for photonic waveguides and an efficient Sommerfeld integral representation [7] for layers. The overall numerical scheme is robust, high-order, and efficient. We demonstrate the performance of the scheme by showing that the effect of finite cladding on modal confinement loss of rectangular dielectric waveguides increases exponentially fast as the thickness decreases.

Notation
We assume that electromagnetic fields are propagated along the z-axis, and that the geometric structure of the photonic waveguides in a layered medium is completely determined by its cross section in the xy-plane (or R 2 ) shown in Fig. 1. We denote the top layer (air for dielectric waveguides) by Ω 1 and its index of refraction by n 1 , the cladding domain by Ω 2 ∈ R 2 and its index of refraction by n 2 , and the bottom substrate domain by Ω 3 ∈ R 2 and its index of refraction by n 3 , respectively. The cross section of the rectangular waveguide is denoted by Ω 0 with n 0 the index of refraction of the core. The boundary of the waveguide is denoted by Γ 0 with ν the unit outward normal vector and τ the unit tangential vector, respectively. Two horizontal lines at y = y t = 0 and y = y b = −h u − h − h l separating the top and bottom layers from the cladding are denoted by Γ t and Γ b , respectively. Here h is the height of the rectangular waveguide, h u and h l are the upper and lower cladding thickness, respectively.

Original PDE Formulation
The electromagnetic field satisfies the Maxwell equations in each domain: where n is the index of refraction of the domain. In the mode calculation, the fundamental assumption is that the electromagnetic field takes the following form: where β is the propagation constant and ω is the frequency of the incident wave. Combining this assumption with the Maxwell equations, we observe that every component of E(x, y) and H(x, y) satisfies the two dimensional Helmholtz equation in each domain: where k = nk v , and k v = ω √ ǫ 0 µ 0 = ω/c is the wave number in vacuum. On the material interface, the boundary conditions are that the tangential components of the electromagnetic fields are continuous. This leads to the following four boundary conditions on each boundary: where [·] denotes the jump of the quantity across the boundary.

Layer Potentials, Sommerfeld Integral and Representation
Let Ω be a bounded domain in R 2 with boundary Γ. We denote the points in R 2 by P and Q. Here Q is usually a point (the source point) on the boundary Γ, and P is in general an arbitrary point (the target point) in R 2 . The Green's function for (3) is where H (1) 0 is the Hankel function of the first kind of order zero [11]. Let σ be a square integrable function on Γ. We define the single, double, and anti-double layer potentials by the following formulas, respectively: Let P = (x, y) and Q = (x 0 , y 0 ). Then the Green's function in (5) has the following Sommerfeld integral representation [7]: Using the above representation, the layer potentials defined in (6) and their derivatives have an alternative representation for sources lying on the interface of a layered medium. For example, the single, double, and anti-double layer potentials for the top layer with sources on Γ t (y = y t ) have the following Sommerfeld representations: where σ is the Fourier transform of the unknown density σ on Γ t . The advantage of the Sommerfeld representation is that the kernels in (8) decay exponentially fast as λ approaches infinity (while the kernels using the Green's function in (6) decays only algebraically).

Integral Representations of the Electromagnetic Fields
Suppose that J(x, y) = J z k +J τ τ and M(x, y) = M z k +M τ τ are two unknown surface currents on Γ. Following [6], we define two vector fields φ k and When Γ is a horizontal line, we also define similar expressions for φ with the associate layer potentials replaced by their corresponding Sommerfeld representations. We now propose the following representations for the electromagnetic fields in various regions for the mode calculation of photonic waveguide in a layered medium shown in Fig. 1.
In other words, the fields inside the core are generated by the unknown densities J and M on its boundary Γ 0 via the representations (9)- (10). The fields in the top layer are generated by the unknown densities J t and M t in the Fourier domain via the Sommerfeld representation of (9)- (10). The fields in the bottom layer are generated by the unknown densities J b and M b in the Fourier domain. Finally, the fields in the cladding region are generated by all six unknown vector densities via suitable representations.
It is tedious, yet straightforward to show [6] that the above representation satisfies the corresponding Helmholtz equation and also the Maxwell equations in each region. The boundary conditions on Γ t , Γ 0 , and Γ b together with the well-known jump relations of the layer potentials (see, for example, [5,6]) then lead to a 12 × 12 block system Ax = 0, where Similar arguments in [6] show that A is a second kind Fredholm integral operator for smooth boundaries. And the propagation constant β of the eletromagnetic mode is a complex number for which the integral operator A has a nontrivial nullspace. We would like to emphasize again that our formulation is well-conditioned and thuse capable of achieving high accuracy even in the case of nonsmooth geometries, say, waveguides with corners, while it is difficult to obtain high accuracy for nonsmooth cases using existing boundary integral methods [1,2,4,8,12] due to the intrinsic ill-conditioning of their formulations.

Discretization of the Layer Potentials
The layer potentials involve weakly singular integrals and many photonic waveguides in integrated optics are of rectangular shape. Thus we need to deal with the singularities in the kernel and the densities (induced by the corner singularity in the geometry). Here we use a collocation Nyström method with dyadic refinement toward the corners to discretize the layer potentials.
We first divide each side of the rectangle into N m +2 subintervals of equal length, then divide each end subinterval dyadically into N e smaller and smaller subintervals. On each subinterval, we place p shifted and scaled Gauss-Legendre nodes and the solution is approximated by a polynomial of degree less than p. So the total number of discretization points N on each side is p · (N m + 2N e ). For each collocation point, the integrals in the layer potentials are discretized via either a precomputed generalized Gaussian quadrature when they are weakly or nearly singular or regular Gaussian quadrature when they are smooth.

Discretization of the Sommerfeld Representation
To numerically evaluate the Sommerfeld integral to high order, we need to avoid the square root singularity in (7). This can be achieved by contour deformation. In particular, we choose the following hyperbolic tangent contour (see Figure 2) We then apply a truncated trapezoidal rule to discretize the Sommerfeld integrals, which The value of T depends on the wave numbers and the distance between the core of the waveguide and the interface separating the layers. In our numerical experiments, we set T = 13, after which the contribution from Sommerfeld integral is exponentially small. We refer the reader to [7] for a more detailed discussions of the advantages of the Sommerfeld representation as compared with the original representation using Green's function directly.

Finding Propagation Constant via Root Finding
We apply the method in [2] to find the propagation constant (or the effective index n e = β/k v in the actual implementation). Suppose that M (β) is the resulting matrix from the discretization. The propagation constant is obtained by finding the roots of a scalar function f (β) = 1/ u T M −1 (β)v via Müller's method [9]. Here u and v are two fixed random column vectors.

Numerical Results
In this section, we present some benchmark calculation on the effect of lower cladding thickness on the modal confinement loss of rectangular dielectric waveguides. Example 1: a high refractive index contrast silica waveguide. In this example, the cross section of the waveguide is of 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 2 = 1.4447 × 1.02. The refractive indices of the silicon base and the air are 3.476, and 1.0003, respectively. The wavelength of the incident field is 1550nm. In the simplified model where the top and bottom layers are absent, our calculation in [6] shows that the waveguide supports a mode with double degeneracy with the effective index n e ≃ 1.458601414886. The result is accurate to 12 digits (see [6] for details).
We first check the convergence rate of our numerical scheme. Table 1 shows the effective index of the first mode found by our scheme for various number of discretization points on each side of the square when h l = 4µm. The number of points in the Sommerfeld representation is fixed N S = 1000 as we found increasing N S will give about the same values under double precision computation. We observe that 12 digit accuracy is achieved with N = 500 (to be more precise, N m = 10, N e = 20, p = 10) points. Thus we set N = 500 for our subsequent calculation.  Table 1: Convergence study for the effective index of the second mode when the lower cladding thickness is 4µm. 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 of the second mode.
We now consider the effect of the finite thickness of the cladding. To simplify our discussion, the upper thickness of the cladding is fixed at 15µm, while the thickness of the lower cladding h l is varied from 10µm to 4µm. The presence of the layers destroys the symmetry of the square waveguide, and the doubly degenerate mode is split into two single modes. Table 2 lists the effective indices of two modes for various lower cladding thickness. The results are obtained by setting N S = 1000 for Sommerfeld integrals and N = 500 for each side of the waveguide, which achieves about 12-digit accuracy by the aforementioned convergence study.
We now calculate the modal confinement loss (in dB/m) [15] via the formula L = 20 ln (10) · 2π λ · Im(n e ) · 10 9  with λ the vacuum wavelength measured in nanometers. Figure 3 plots out the modal confinement loss L versus h l and their least squares fits.
Lower cladding thickness(µm) The least squares fit of both data sets shows that the slope is about K 1 ≈ −0.728µm −1 . That is, the modal confinement loss L increases exponentially fast as h l is decreasing: This can be explained as follows. The electromagnetic fields decays in the cladding region in the form F ∝ e −2π √ n 2 2 −n 2 e |y|/λ . Thus, the energy loss, which is directly linked to the modal confinement loss, due to the finite lower cladding thickness is proportional to e h l /λ = 10 −4π √ n 2 2 −n 2 e (log 10 e)h l /λ . With n 2 = 1.4447 × 1.02, n e ≈ 1.4586, λ = 1.55µm, we have which is in good agreement with the value of K 1 (about 1.3% relative error).
We have plotted the magnitude of the electromagnetic field |E z | and |H Z | of the two modes in Figure 4. We observe that the field is concentrated near the core and decays exponentially fast away from the core, which confirms our previous claim about the behavior of the field in the cladding region.  Example 2: a low refractive index contrast silica waveguide. In this example, the cross section of the waveguide is of square shape with the side length equal to 5.2µm. The refractive indices of the cladding, core, silicon base, and the air are 1.4447, 1.4447 × 1.0075, 3.476, and 1.0003, respectively. The wave length of the incident field is 1550nm. The upper thickness of the cladding is fixed at 15µm. And we vary the thickness of the lower cladding from 4µm to 12µm. The discretization is the same as Example 1. Table 3 Figure 5 shows the propagation loss L with respect to the lower cladding thickness and their least squares line fit.
The least squares line fit of both data sets shows that the slope is about K 2 ≈ −0.431µm −1 . That is, the modal confinement loss L increases exponentially fast as h l is decreasing: Similar argument as in Example 1 shows that the energy loss due to the finite lower cladding thickness is proportional to |F| 2 ∝ 10 S 2 h l . With n 2 = 1.4447 × 1.0075, n e ≈ 1.449465, λ = 1.55µm, we have S 2 = −4π n 2 2 − n 2 e · (log 10 e)/λ ≈ −0.468µm −1 , which is in good agreement with the value of K 2 (about 7.8% relative error). A comparison of these two examples show that high contrast silica waveguides not only have smaller core size, but also need much thinner cladding for the same modal confinement loss. Hence, more compact photonic devices can be made using high contrast silica waveguides.

Conclusions
We have presented a well-conditioned integral formulation for the calculation of electromagnetic modes of photonic waveguides in layered media. Unlike finite difference or finite element methods which requires the volume discretization and the imposition of artificial boundary conditions. Our numerical scheme discretizes the material interface only and treats the effect of layers efficiently via the Sommerfeld representation, leading to a discrete linear system with Lower cladding thickness(µm) very small number of unknowns. The scheme is as well conditioned as the underlying physical problem and thus capable of achieving high accuracy even for waveguides having complex geometries and corners. Our benchmark computation on rectangular waveguides imbedded in a cladding of finite thickness provides more than 10 significant digits for the propagation constants and shows quantitatively the relationship between the modal confinement loss and the cladding thickness. The scheme, being highly accurate and efficient for the mode calculation, provides a powerful design and simulation tool for the photonics industry.