## Abstract

There has been increasing interest in making 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 of accuracy for the computation of propagation constant *β* since the modal confinement loss is proportional to the imaginary part of *β* 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 the electromagnetic modes of photonic waveguides with an 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 embedded in a cladding material of finite thickness. The numerical results, with a more than 10-digit accuracy, show quantitatively that the modal confinement loss of the rectangular waveguide increases exponentially fast as the cladding thickness decreases.

© 2016 Optical Society of America

## 1. 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 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 [1, 2] or finite element methods [3, 4] require the discretization of a finite domain, supplemented by artificial boundary conditions, such as perfectly matched layers, to simulate the effect of an infinite 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 [5–9] 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 [10] for photonic waveguides and an efficient Sommerfeld integral representation [11] 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.

## 2. Integral formulation for the mode calculation

#### 2.1. 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 ℝ^{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} ∈ ℝ^{2} and its index of refraction by *n*_{2}, and the bottom substrate domain by Ω_{3} ∈ ℝ^{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*separating the top and bottom layers from the cladding are denoted by Γ

_{l}*and Γ*

_{t}*, respectively. Here*

_{b}*h*is the height of the rectangular waveguide,

*h*and

_{u}*h*are the upper and lower cladding thickness, respectively.

_{l}#### 2.2. Original PDE formulation

The electromagnetic field satisfies the Maxwell equations in each domain:

*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:

*β*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*, and ${k}_{\nu}=\omega \sqrt{{\mathit{\u03f5}}_{0}{\mu}_{0}}=\omega /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:

_{ν}#### 2.3. Layer potentials, Sommerfeld integral and representation

Let Ω be a bounded domain in ℝ^{2} with boundary Γ. We denote the points in ℝ^{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 ℝ^{2}. The Green’s function for (3) is

*σ*be a square integrable function on Γ. We define the single, double, and anti-double layer potentials by the following formulas, respectively:

*P*= (

*x*,

*y*) and

*Q*= (

*x*

_{0},

*y*

_{0}). Then the Green’s function in (5) has the following Sommerfeld integral representation [11]:

*(*

_{t}*y*=

*y*) have the following Sommerfeld representations:

_{t}*σ*on Γ

*. The advantage of the Sommerfeld representation is that the kernels in (8) decay exponentially fast as*

_{t}*λ*approaches infinity (while the kernels using the Green’s function in (6) decay only algebraically).

#### 2.4. Integral representations of the electromagnetic fields

Suppose that $\mathbf{J}\left(x,y\right)={\mathbf{J}}_{z}\widehat{k}+{\mathbf{J}}_{\tau}\tau $ and $\mathbf{M}\left(x,y\right)={\mathbf{M}}_{z}\widehat{k}+{\mathbf{M}}_{\tau}\tau $ are two unknown surface currents on Γ. Following [10], we define two vector fields ${\varphi}_{\mathrm{\Gamma}}^{k}\left[\mathbf{J},\mathbf{M}\right]$, ${\psi}_{\mathrm{\Gamma}}^{k}\left[\mathbf{J},\mathbf{M}\right]$ by the following formulas:

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
${\widehat{\mathbf{J}}}_{t}$ and
${\widehat{\mathbf{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
${\widehat{\mathbf{J}}}_{b}$ and
${\widehat{\mathbf{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 [10] that the above representation satisfies the corresponding Helmholtz equation and also the Maxwell equations in each region. The boundary conditions on Γ* _{t}*, Γ

_{0}, and Γ

*together with the well-known jump relations of the layer potentials (see, for example, [10, 13]) then lead to a 12 × 12 block system*

_{b}*Ax*= 0, where $x={\left[{\widehat{J}}_{t,x}\phantom{\rule{0.2em}{0ex}}{\widehat{J}}_{t,z}\phantom{\rule{0.2em}{0ex}}{\widehat{M}}_{t,x}\phantom{\rule{0.2em}{0ex}}{\widehat{M}}_{t,z}\phantom{\rule{0.2em}{0ex}}{J}_{\tau}\phantom{\rule{0.2em}{0ex}}{J}_{z}\phantom{\rule{0.2em}{0ex}}{M}_{\tau},\phantom{\rule{0.2em}{0ex}}{M}_{z},\phantom{\rule{0.2em}{0ex}}{\widehat{J}}_{b,x}\phantom{\rule{0.2em}{0ex}}{\widehat{J}}_{b,z}\phantom{\rule{0.2em}{0ex}}{\widehat{M}}_{b,x}\phantom{\rule{0.2em}{0ex}}{\widehat{M}}_{b,z}\right]}^{T}$. Similar arguments in [10] show that

*A*is a second kind Fredholm integral operator for smooth boundaries. And the propagation constant

*β*of the electromagnetic 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 thus 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 [5–9] due to the intrinsic ill-conditioning of their formulations.

## 3. Numerical algorithm

#### 3.1. 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*smaller and smaller subintervals. On each subinterval, we place

_{e}*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*+2

_{m}*N*). 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.

_{e}#### 3.2. 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 Fig. 2)

We then apply a truncated trapezoidal rule to discretize the Sommerfeld integrals, which achieves spectral accuracy due to the smoothness and exponential decay of the integrand. Specifically, we truncate *t* at [−*T*, *T*] and discretize *t* uniformly by *N _{S}* points on that interval. 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 [11] for a more detailed discussions of the advantages of the Sommerfeld representation as compared with the original representation using Green’s function directly.

#### 3.3. Finding propagation constant via root finding

We apply the method in [6] to find the propagation constant (or the effective index *n _{e}* =

*β*/

*k*in the actual implementation). Suppose that

_{v}*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 [14]. Here

*u*and

*v*are two fixed random column vectors.

## 4. 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 1550*nm*. In the simplified model where the top and bottom layers are absent, our calculation in [10] 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 [10] 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*= 1000 as we found increasing

_{S}*N*will give about the same values under double precision computation. We observe that 12 digit accuracy is achieved with

_{S}*N*= 500 (to be more precise,

*N*= 10,

_{m}*N*= 20,

_{e}*p*= 10) points. Thus we set

*N*= 500 for our subsequent calculation.

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

*λ*the vacuum wavelength measured in nanometers. Figure 3 plots out the modal confinement loss

*L*versus

*h*and their least squares fits.

_{l}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
$\mathbf{F}\propto {e}^{-2\pi \sqrt{{n}_{2}^{2}-{n}_{e}^{2}}\left|y\right|/\lambda}$. Thus, the energy loss, which is directly linked to the modal confinement loss, due to the finite lower cladding thickness is proportional to
${\left|\mathbf{F}\right|}^{2}\propto {e}^{-4\pi \sqrt{{n}_{2}^{2}-{n}_{e}^{2}}{h}_{l}/\lambda}={10}^{-4\pi \sqrt{{n}_{2}^{2}-{n}_{e}^{2}}\left({\mathrm{log}}_{10}e\right){h}_{l}/\lambda}$. With *n*_{2} = 1.4447 × 1.02, *n _{e}* ≈ 1.4586,

*λ*= 1.55

*μm*, we have

*K*

_{1}(about 1.3% relative error).

We have plotted the magnitude of the electromagnetic field |*E _{z}*| and |

*H*| of the two modes in Fig. 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.

_{Z}#### 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 1550*nm*. 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 lists the effective indices of two modes for various lower cladding thickness.

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:

*n*

_{2}= 1.4447 × 1.0075,

*n*≈ 1.449465,

_{e}*λ*= 1.55

*μm*, we have

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

## 5. 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 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.

## Funding

National Science Foundation (NSF) (DMS-1418918); Office of the Assistant Secretary of Defense for Research and Engineering and AFOSR (FA9550-10-1-0180).

## Acknowledgments

S. J. is supported by the National Science Foundation (NSF) under grant DMS-1418918. J. L. was supported in part by the Office of the Assistant Secretary of Defense for Research and Engineering and AFOSR under grant FA9550-10-1-0180. The authors would like to thank Prof. Leslie Greengard at Courant Institute for useful suggestions.

## References and links

**1. **N. Thomas, P. Sewell, and T. M. Benson, “A new full-vectorial higher order finite-difference scheme for the modal analysis of rectangular dielectric waveguides,” J. Lightwave Technol. **25**, 2563–2570 (2007). [CrossRef]

**2. **Y. P. Chiou, Y. C. Chiang, C. H. Lai, C. H. Du, and H. C. Chang, “Finite difference modeling of dielectric waveguides with corners and slanted facets,” J. Lightwave Technol. **27**, 2077–2086 (2009). [CrossRef]

**3. **S. Selleri, L. Vincetti L, A. Cucinotta, and M. Zoboli, “Complex FEM modal solver of optical waveguides with PML boundary conditions,” Opt. Quant. Electron. **33**, 359–371 (2001). [CrossRef]

**4. **S. S. A. Obayya, B. M. A. Rahman, K. T. V. Grattan, and H. A. El-Mikati, “Full vectorial finite-element-based imaginary distance beam propagation solution of complex modes in optical waveguides,” J. Lightwave Technol. **20**, 1054–1060 (2002). [CrossRef]

**5. **S. V. Boriskina, T. M. Benson, P. Sewell, and A. I. Nosich, “Highly efficient full-vectorial integral equations method solution for the bound, leaky, and complex modes of dielectric waveguides,” IEEE J. Selected Topics in Quantum Electron. **8**, 1225–1231 (2002). [CrossRef]

**6. **H. Cheng, W. Y. Crutchfield, M. Doery, and L. Greengard, “Fast, accurate integral equation methods for the analysis of photonic crystal fibers I: Theory,” Opt. Express **12**(16), 3791–3805 (2004). [CrossRef] [PubMed]

**7. **A. Ferrando, E. Silvestre, J.J. Miret, P. Andres, and M.V. Andres, “Full vector analysis of a realistic photonic crystal fiber,” Opt. Lett. **24**, 276–278 (1999). [CrossRef]

**8. **W. Lu and Y. Y. Lu, “Waveguide mode solver based on Neumann-to-Dirichlet operators and boundary integral equations,” J. Comput. Phys. **231**, 1360–1371 (2012). [CrossRef]

**9. **E. Pone, A. Hassani, S. Lacroix, A. Kabashin, and M. Skorobogatiy, “Boundary integral method for the challenging problems in bandgap guiding, plasmonics and sensing,” Opt. Express **15**, 10231–10246 (2007). [CrossRef] [PubMed]

**10. **J. Lai and S. Jiang, “Second kind integral equation formulation for the mode calculation of optical waveguides,” Appl. Comput. Harmon. Anal. in press, doi:http://dx.doi.org/10.1016/j.acha.2016.06.009.

**11. **J. Lai, M. Kobayashi, and L. Greengard, “A fast solver for multi-particle scattering in a layered medium,” Opt. Express **14**, 302–307 (2014).

**12. **F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, eds. *NIST Handbook of Mathematical Functions* (Cambridge University Press, 2010).

**13. **R. Kress, *Linear Integral Equations* (Springer–Verlag, 1989). [CrossRef]

**14. **D. E. Müller, “A method for solving algebraic equations using an automatic computer,” Math. Tables Aids Comput. **10**, 208–215 (1956). [CrossRef]

**15. **T. P. White, B. T. Kuhlmey, R. C. McPhedran, D. Maystre, G. Renversez, C. Martijn de Sterke, and L. C. Botten, “Multipole method for microstructured optical fibers. I. formulation,” J. Opt. Soc. Am. B **19**(10), 2322–2330 (2002). [CrossRef]