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
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  for photonic waveguides and an efficient Sommerfeld integral representation  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
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 n1, the cladding domain by Ω2 ∈ ℝ2 and its index of refraction by n2, and the bottom substrate domain by Ω3 ∈ ℝ2 and its index of refraction by n3, respectively. The cross section of the rectangular waveguide is denoted by Ω0 with n0 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 = yt = 0 and y = yb = −hu − h − hl 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, hu and hl are the upper and lower cladding thickness, respectively.
2.2. Original PDE formulation
The electromagnetic field satisfies the Maxwell equations in each domain:
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) is12]. Let σ be a square integrable function on Γ. We define the single, double, and anti-double layer potentials by the following formulas, respectively: (5) has the following Sommerfeld integral representation : (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 = yt) have the following Sommerfeld representations: (8) decay exponentially fast as λ 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 and are two unknown surface currents on Γ. Following , we define two vector fields , by the following formulas: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 and in the Fourier domain via the Sommerfeld representation of (9)–(10). The fields in the bottom layer are generated by the unknown densities and 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  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, [10, 13]) then lead to a 12 × 12 block system Ax = 0, where . Similar arguments in  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 Nm +2 subintervals of equal length, then divide each end subinterval dyadically into Ne 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· (Nm +2Ne). 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.
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 NS 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  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  to find the propagation constant (or the effective index ne = β/kv 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/(uT M−1(β)v) via Müller’s method . 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 n0 = 1.4447, while that of the core is 2% higher, i.e., n2 = 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  shows that the waveguide supports a mode with double degeneracy with the effective index ne ≃ 1.458601414886. The result is accurate to 12 digits (see  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 hl = 4μm. The number of points in the Sommerfeld representation is fixed NS = 1000 as we found increasing NS 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, Nm = 10, Ne = 20, 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 hl 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 NS = 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)  via the formulaFigure 3 plots out the modal confinement loss L versus hl and their least squares fits.
The least squares fit of both data sets shows that the slope is about K1 ≈ −0.728μm−1. That is, the modal confinement loss L increases exponentially fast as hl is decreasing:
This can be explained as follows. The electromagnetic fields decays in the cladding region in the form . Thus, the energy loss, which is directly linked to the modal confinement loss, due to the finite lower cladding thickness is proportional to . With n2 = 1.4447 × 1.02, ne ≈ 1.4586, λ = 1.55μm, we have
We have plotted the magnitude of the electromagnetic field |Ez| and |HZ| 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.
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 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 K2 ≈ −0.431μm−1. That is, the modal confinement loss L increases exponentially fast as hl is decreasing:
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.
National Science Foundation (NSF) (DMS-1418918); Office of the Assistant Secretary of Defense for Research and Engineering and AFOSR (FA9550-10-1-0180).
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]