Quantitative study of the effect of cladding thickness on modal confinement loss in photonic waveguides

Shidong Jiang; Jun Lai

doi:10.1364/OE.24.024872

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 ℝ²) shown in Fig. 1. We denote the top layer (air for dielectric waveguides) by Ω₁ and its index of refraction by n₁, the cladding domain by Ω₂ ∈ ℝ² and its index of refraction by n₂, and the bottom substrate domain by Ω₃ ∈ ℝ² and its index of refraction by n₃, respectively. The cross section of the rectangular waveguide is denoted by Ω₀ with n₀ the index of refraction of the core. The boundary of the waveguide is denoted by Γ₀ 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.

Fig. 1 Cross section of a rectangular waveguide in layered medium.

Download Full Size | PDF

2.2. Original PDE formulation

The electromagnetic field satisfies the Maxwell equations in each domain:

{\begin{matrix} \nabla \times E = - μ_{0} \frac{\partial H}{\partial t}, \nabla \cdot E = 0, \\ \nabla \times H = - ϵ_{0} n^{2} \frac{\partial E}{\partial t}, \nabla \cdot H = 0, \end{matrix}

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:

[E (x, y, z, t), H (x, y, z, t)] = [E (x, y), H (x, y)] e^{i (β z - ω t)},

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:

[Δ + (k^{2} - β^{2})] u = 0,

where k = nk_ν, and

k_{ν} = ω \sqrt{ϵ_{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:

[E_{z}] = 0, [H_{z}] = 0, [E_{τ}] = 0, [H_{τ}] = 0,

where [·] denotes the jump of the quantity across the boundary.

2.3. Layer potentials, Sommerfeld integral and representation

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

G^{k} (P, Q) = \frac{i}{4} H_{0}^{(1)} (\sqrt{k^{2} - β^{2}} ‖ P - Q ‖),

where

H_{0}^{(1)}

is the Hankel function of the first kind of order zero [12]. Let σ be a square integrable function on Γ. We define the single, double, and anti-double layer potentials by the following formulas, respectively:

{\begin{array}{l} S_{Γ}^{k} [σ] (P) = \int_{Γ} G^{k} (P, Q) σ (Q) d s_{Q}, \\ D_{Γ}^{k} [σ] (P) = \int_{Γ} \frac{\partial G^{k} (P, Q)}{\partial ν (Q)} σ (Q) d s_{Q}, \\ T_{Γ}^{k} [σ] (P) = \int_{Γ} \frac{\partial G^{k} (P, Q)}{\partial τ (Q)} σ (Q) d s_{Q} . \end{array}

Let P = (x, y) and Q = (x₀, y₀). Then the Green’s function in (5) has the following Sommerfeld integral representation [11]:

G^{k} (P, Q) = \frac{1}{4 π} \int_{- \infty}^{\infty} \frac{e^{- \sqrt{λ^{2} - k^{2} + β^{2}} | y - y_{0} |}}{\sqrt{λ^{2} - k^{2} + β^{2}}} e^{i λ (x - x_{0})} d λ .

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:

{\begin{matrix} {\hat{S}}_{Γ_{t}}^{k_{1}} [\hat{σ}] (P) = \frac{1}{4 π} \int_{- \infty}^{\infty} \frac{e^{- \sqrt{λ^{2} - {k_{1}}^{2} + β^{2}} | y - y_{t} |}}{\sqrt{λ^{2} - {k_{1}}^{2} + β^{2}}} e^{i λ x} \hat{σ} (λ) d λ, \\ {\hat{D}}_{Γ_{t}}^{k_{1}} [\hat{σ}] (P) = \frac{1}{4 π} \int_{- \infty}^{\infty} e^{- \sqrt{λ^{2} - {k_{1}}^{2} + β^{2}} | y - y_{t} |} e^{i λ x} \hat{σ} (λ) d λ, \\ {\hat{T}}_{Γ_{t}}^{k_{1}} [\hat{σ}] (P) = \frac{1}{4 π} \int_{- \infty}^{\infty} \frac{i λ e^{- \sqrt{λ^{2} - {k_{1}}^{2} + β^{2}} | y - y_{t} |}}{\sqrt{λ^{2} - {k_{1}}^{2} + β^{2}}} e^{i λ x} \hat{σ} (λ) d λ, \end{matrix}

where

\hat{σ}

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) decay only algebraically).

2.4. Integral representations of the electromagnetic fields

Suppose that $J (x, y) = J_{z} \hat{k} + J_{τ} τ$ and $M (x, y) = M_{z} \hat{k} + M_{τ} τ$ are two unknown surface currents on Γ. Following [10], we define two vector fields $ϕ_{Γ}^{k} [J, M]$ , $ψ_{Γ}^{k} [J, M]$ by the following formulas:

{\begin{array}{l} ϕ_{Γ, x}^{k} [J, M] = - \frac{1}{i k_{ν}} \frac{\partial}{\partial x} Τ_{Γ}^{k} [J_{τ}] + \frac{β}{k_{ν}} \frac{\partial}{\partial x} S_{Γ}^{k} [J_{z}] + \frac{k^{2}}{i k_{ν}} S_{Γ}^{k} [J_{τ} τ_{x}] - \frac{\partial}{\partial y} S_{Γ}^{k} [M_{z}] + i β S_{Γ}^{k} [M_{τ} τ_{y}], \\ ϕ_{Γ, y}^{k} [J, M] = - \frac{1}{i k_{ν}} \frac{\partial}{\partial y} Τ_{Γ}^{k} [J_{τ}] + \frac{β}{k_{ν}} \frac{\partial}{\partial y} S_{Γ}^{k} [J_{z}] + \frac{k^{2}}{i k_{ν}} S_{Γ}^{k} [J_{τ} τ_{y}] + \frac{\partial}{\partial x} S_{Γ}^{k} [M_{z}] - i β S_{Γ}^{k} [M_{τ} τ_{x}], \\ ϕ_{Γ, z}^{k} [J, M] = - \frac{β}{k_{ν}} Τ_{Γ}^{k} [J_{τ}] + \frac{(k^{2} - β^{2})}{i k_{ν}} S_{Γ}^{k} [J_{z}] + D_{Γ}^{k} [M_{τ}] . \end{array}

and

{\begin{array}{l} ψ_{Γ, x}^{k} [J, M] = \frac{1}{i k_{ν}} \frac{\partial}{\partial x} Τ_{Γ}^{k} [M_{τ}] - \frac{β}{k_{ν}} \frac{\partial}{\partial x} S_{Γ}^{k} [M_{z}] - \frac{k^{2}}{i k_{ν}} S_{Γ}^{k} [M_{τ} τ_{x}] - \frac{k^{2}}{k_{ν}^{2}} \frac{\partial}{\partial y} S_{Γ}^{k} [J_{z}] + i β \frac{k^{2}}{k_{ν}^{2}} S_{Γ}^{k} [J_{τ} τ_{y}], \\ ψ_{Γ, y}^{k} [J, M] = \frac{1}{i k_{ν}} \frac{\partial}{\partial y} Τ_{Γ}^{k} [M_{τ}] - \frac{β}{k_{ν}} \frac{\partial}{\partial y} S_{Γ}^{k} [M_{z}] - \frac{k^{2}}{i k_{ν}} S_{Γ}^{k} [M_{τ} τ_{y}] + \frac{k^{2}}{k_{ν}^{2}} \frac{\partial}{\partial x} S_{Γ}^{k} [J_{z}] - i β \frac{k^{2}}{k_{ν}^{2}} S_{Γ}^{k} [J_{τ} τ_{x}], \\ ψ_{Γ, z}^{k} [J, M] = \frac{β}{k_{ν}} Τ_{Γ}^{k} [M_{τ}] - \frac{(k^{2} - β^{2})}{i k_{ν}} S_{Γ}^{k} [M_{z}] + \frac{k^{2}}{k_{ν}^{2}} D_{Γ}^{k} [J_{τ}] . \end{array}

When Γ is a horizontal line, we also define similar expressions for

{\hat{ϕ}}_{Γ}^{k} [\hat{J}, \hat{M}]

,

{\hat{ψ}}_{Γ}^{k} [\hat{J}, \hat{M}]

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.

{\begin{array}{l} (E, H) = (ϕ_{Γ_{0}}^{k_{0}} [J, M], ψ_{Γ_{0}}^{k_{0}} [J, M]) in Ω_{0}, \\ (E, H) = ({\hat{ϕ}}_{Γ_{t}}^{k_{1}} [{\hat{J}}_{t}, {\hat{M}}_{t}], {\hat{ψ}}_{Γ_{t}}^{k_{1}} [{\hat{J}}_{t}, {\hat{M}}_{t}]) in Ω_{1}, \\ (E, H) = ({\hat{ϕ}}_{Γ_{t}}^{k_{2}} [{\hat{J}}_{t}, {\hat{M}}_{t}], {\hat{ψ}}_{Γ_{t}}^{k_{2}} [{\hat{J}}_{t}, {\hat{M}}_{t}]) + (ϕ_{Γ_{0}}^{k_{2}} [J, M], ψ_{Γ_{0}}^{k_{2}} [J, M]) + ({\hat{ϕ}}_{Γ_{b}}^{k_{2}} [{\hat{J}}_{b}, {\hat{M}}_{b}], {\hat{ψ}}_{Γ_{b}}^{k_{2}} [{\hat{J}}_{b}, {\hat{M}}_{b}]) in Ω_{2}, \\ (E, H) = ({\hat{ϕ}}_{Γ_{b}}^{k_{3}} [{\hat{J}}_{b}, {\hat{M}}_{b}], {\hat{ψ}}_{Γ_{b}}^{k_{3}} [{\hat{J}}_{b}, {\hat{M}}_{b}]) in Ω_{3} . \end{array}

In other words, the fields inside the core are generated by the unknown densities J and M on its boundary Γ₀ via the representations (9)–(10). The fields in the top layer are generated by the unknown densities ${\hat{J}}_{t}$ and ${\hat{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 ${\hat{J}}_{b}$ and ${\hat{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, Γ₀, 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 $x = {[{\hat{J}}_{t, x} {\hat{J}}_{t, z} {\hat{M}}_{t, x} {\hat{M}}_{t, z} J_{τ} J_{z} M_{τ}, M_{z}, {\hat{J}}_{b, x} {\hat{J}}_{b, z} {\hat{M}}_{b, x} {\hat{M}}_{b, z}]}^{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_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.

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)

λ = t - \frac{\tanh (t)}{2} i, t \in ℝ .

Fig. 2 The hyperbolic tangent contour used in the Sommerfeld representation

Download Full Size | PDF

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_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⁻¹(β)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₀ = 1.4447, while that of the core is 2% higher, i.e., n₂ = 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 [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_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.

View Table | View all tables in this article

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.

Table 2. Effective indices versus lower cladding thickness h_l for Example 1.

View Table | View all tables in this article

We now calculate the modal confinement loss (in dB/m) [15] via the formula

L = \frac{20}{\ln (10)} \cdot \frac{2 π}{λ} \cdot I m (n_{e}) \cdot 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.

Fig. 3 Modal confinement loss versus the lower cladding thickness for Example 1. The y-axis uses logarithmic scale (base 10). Dashed lines are the least squares fit to the data.

Download Full Size | PDF

The least squares fit of both data sets shows that the slope is about K₁ ≈ −0.728μm⁻¹. That is, the modal confinement loss L increases exponentially fast as h_l is decreasing:

L \propto 10^{K_{1} h_{l}} .

This can be explained as follows. The electromagnetic fields decays in the cladding region in the form $F \propto e^{- 2 π \sqrt{n_{2}^{2} - n_{e}^{2}} | y | / λ}$ . Thus, the energy loss, which is directly linked to the modal confinement loss, due to the finite lower cladding thickness is proportional to ${| F |}^{2} \propto e^{- 4 π \sqrt{n_{2}^{2} - n_{e}^{2}} h_{l} / λ} = 10^{- 4 π \sqrt{n_{2}^{2} - n_{e}^{2}} (\log_{10} e) h_{l} / λ}$ . With n₂ = 1.4447 × 1.02, n_e ≈ 1.4586, λ = 1.55μm, we have

S_{1} = - 4 π \sqrt{n_{2}^{2} - n_{e}^{2}} \cdot (\log_{10} e) / λ \approx - 0.738 μ m^{- 1},

which is in good agreement with the value of K₁ (about 1.3% relative error).

We have plotted the magnitude of the electromagnetic field |E_z| and |H_Z| 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.

Fig. 4 Magnitude of the electromagnetic field E_z and H_z of the first and second modes with h_l = 4μm in Example 1. The square represents the boundary of the waveguide. The colorbar uses logarithmic scale.

Download Full Size | PDF

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.

Table 3. Effective indices of the photonic waveguide in Example 2.

View Table | View all tables in this article

Figure 5 shows the propagation loss L with respect to the lower cladding thickness and their least squares line fit.

Fig. 5 Modal confinement loss versus the lower cladding thickness for Example 2. The y-axis uses logarithmic scale (base 10). Dashed lines are the least squares fit to the data.

Download Full Size | PDF

The least squares line fit of both data sets shows that the slope is about K₂ ≈ −0.431μm⁻¹. That is, the modal confinement loss L increases exponentially fast as h_l is decreasing:

L \propto 10^{K_{2} h_{l}} .

Similar argument as in Example 1 shows that the energy loss due to the finite lower cladding thickness is proportional to

{| F |}^{2} \propto 10^{S_{2} h_{l}}

. With n₂ = 1.4447 × 1.0075, n_e ≈ 1.449465, λ = 1.55μm, we have

S_{2} = - 4 π \sqrt{n_{2}^{2} - n_{e}^{2}} \cdot (\log_{10} e) / λ \approx - 0.468 μ m^{- 1},

which is in good agreement with the value of K₂ (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]

N	Real	Imaginary
200	1.45860122763585	2.673537E-8
300	1.45860122756971	2.498713E-8
400	1.45860122756758	2.493100E-8
500	1.45860122756751	2.492917E-8
600	1.45860122756751	2.492915E-8

	First mode		Second mode
h_l	Re(n_e)	Im(n_e)	Re(n_e)	Im(n_e)
4	1.45860122757	2.4929E-8	1.45860127384	1.26765E-7
5	1.45860138033	4.578E-9	1.45860138880	2.3309E-8
6	1.45860140847	8.47E-10	1.45860141003	4.317E-9
7	1.45860141369	1.57E-10	1.45860141398	8.04E-10
8	1.45860141466	2.9E-11	1.45860141471	1.50E-10
9	1.45860141484	5E-12	1.45860141485	2.8E-11
10	1.45860141488	1E-12	1.45860141488	5E-12

	First mode		Second mode
h_l	Re(n_e)	Im(n_e)	Re(n_e)	Im(n_e)
4	1.449463952659	1.34678E-7	1.449464113248	7.43394E-7
5	1.449465037025	4.9161E-8	1.449465094948	2.71536E-7
6	1.449465433894	1.8051E-8	1.449465454711	9.9754E-8
7	1.449465579918	6.659E-9	1.449465587215	3.6813E-8
8	1.449465633876	2.465E-9	1.449465636215	1.3635E-8
9	1.449465653883	9.15E-10	1.449465654396	5.066E-9
10	1.449465661324	3.41E-10	1.449465661161	1.887E-9
11	1.449465664098	1.27E-10	1.449465663685	7.04E-10
12	1.449465665135	4.7E-11	1.449465664628	2.63E-10

N	Real	Imaginary
200	1.45860122763585	2.673537E-8
300	1.45860122756971	2.498713E-8
400	1.45860122756758	2.493100E-8
500	1.45860122756751	2.492917E-8
600	1.45860122756751	2.492915E-8

	First mode		Second mode
h_l	Re(n_e)	Im(n_e)	Re(n_e)	Im(n_e)
4	1.45860122757	2.4929E-8	1.45860127384	1.26765E-7
5	1.45860138033	4.578E-9	1.45860138880	2.3309E-8
6	1.45860140847	8.47E-10	1.45860141003	4.317E-9
7	1.45860141369	1.57E-10	1.45860141398	8.04E-10
8	1.45860141466	2.9E-11	1.45860141471	1.50E-10
9	1.45860141484	5E-12	1.45860141485	2.8E-11
10	1.45860141488	1E-12	1.45860141488	5E-12

Quantitative study of the effect of cladding thickness on modal confinement loss in photonic waveguides

Abstract

1. Introduction

2. Integral formulation for the mode calculation

2.1. Notation

2.2. Original PDE formulation

2.3. Layer potentials, Sommerfeld integral and representation

2.4. Integral representations of the electromagnetic fields

3. Numerical algorithm

3.1. Discretization of the layer potentials

3.2. Discretization of the Sommerfeld representation

3.3. Finding propagation constant via root finding

4. Numerical results

Example 1: a high refractive index contrast silica waveguide

Example 2: a low refractive index contrast silica waveguide

5. Conclusions

Funding

Acknowledgments

References and links

Cited By

Figures (5)

Tables (3)

Equations (17)

Optics Express