Complex coupled mode theory electromagnetic mode solver

Timothy DeWolf; Reuven Gordon

doi:10.1364/OE.25.028337

1. Introduction

We are interested in the propagation of electromagnetic radiation in dielectric and metal waveguides [1–4], and specifically in calculating the electromagnetic modes in these systems. Analytic solutions are available for simple systems such as step-index fibers [5], but numerical methods are often necessary when the structure of the waveguide becomes more complex. Existing numerical methods include finite-difference eigenmode (FDE) solvers [6–13], finite-element solvers [14–18], and solvers based on series expansions [19–33]; review articles summarize some of the extensive work done in this field [34,35]. A non-vectorial version of (non-complex) coupled mode theory has been used in the case of lossless graded-index fibers [36,37].

We introduce a method for solving waveguide modes based on complex coupled mode theory (CCMT) [38]. CCMT is traditionally used to describe how waveguide modes couple in the presence of a refractive index perturbation. The introduced coupling allows energy transfer between the modes, and CCMT is an efficient tool for describing this energy transfer [39]. CCMT, applied in this way, is rigorous–the solutions thus found are exact solutions. The method presented here uses CCMT in a different way.

Instead of coupling the true waveguide modes, the CCMT eigenmode solver uses a known set of orthogonal electromagnetic modes as a basis. The full permittivity map of the waveguide (e.g. a fiber) is used to compute the coupling perturbation of CCMT. Being based on the complex variant of coupled mode theory means that lossy materials and metals are supported. The output set of coupled modes computed by the CCMT mode solver are the modes of the waveguide. We compare our results with analytical results and numerical finite difference eigenmode solutions.

As well as being based directly on CCMT, a difference between this method and series methods based on the Fourier transform is the ability to choose an arbitrary basis set that fits well with the geometry of the waveguide under investigation. An example basis set that can be used in the CCMT mode solver is the set of electromagnetic modes associated with a metallic pipe waveguide [40], where radial dependence is expressed in terms of Bessel functions. Another example is the set of electromagnetic modes of an elliptical waveguide [41,42], expressed in terms of Mathieu functions. It is also possible to use the modes of a rectangular metal waveguide [40], although this is similar to Fourier transform methods as it involves a sinusoidal basis.

2. Formulation of complex coupled mode theory

In this section, we present CCMT, introduce a convenient set of basis modes, and present the matrix eigenvalue form of CCMT used to numerically compute waveguide eigenmodes.

2.1. Complex coupled mode theory

The modes of a waveguide with scalar (isotropic) permittivity ∊ are e_k and h_k. The e^iωt time dependence has been factored out of these modal fields; ω is the frequency of the field. Perturbing the waveguide permittivity by Δ∊ such that the new waveguide permittivity is ∊̃ = ∊ + Δ∊ causes the modes e_k and h_k to couple. The coupled electric and magnetic fields E and H in the waveguide are,

E_{t} (x, y, z) = (a_{k} (z) + b_{k} (z)) e_{k t} (x, y)

H_{t} (x, y, z) = (a_{k} (z) - b_{k} (z)) h_{k t} (x, y) .

The subscript t denotes the vector field components in the transverse plane (x and y components). Summation over k is implied (Einstein summation convention). Expressions for the longitudinal (z) component of the coupled field are

E_{z} (x, y, z) = (a_{k} (z) - b_{k} (z)) \frac{∊}{\tilde{∊}} e_{k z} (x, y)

H_{z} (x, y, z) = (a_{k} (z) + b_{k} (z)) h_{k z} (x, y) .

CCMT specifies how the coefficients a_k(z) and b_k(z), describing the amplitude of forwards and backwards propagating modes, are computed. The equations that determine these coefficients are [38]

N_{j} (\frac{\partial a_{j} (z)}{\partial z} + i β_{j} a_{j} (z)) = - i κ_{j k} a_{k} (z) - i χ_{j k} b_{k} (z)

N_{j} (\frac{\partial b_{j} (z)}{\partial z} - i β_{j} b_{j} (z)) = i κ_{j k} a_{k} (z) + i χ_{j k} b_{k} (z) .

These equations give the time dependence of the system. The β_j are the wavevectors associated with each mode j of the uncoupled basis. The entries of the matrices κ_jk and χ_jk are the CCMT overlap integrals

κ_{j k} = \frac{ω}{2} \int d A Δ ∊ (e_{j t} \cdot e_{k t} - \frac{∊}{\tilde{∊}} e_{j z} \cdot e_{k z})

χ_{j k} = \frac{ω}{2} \int d A Δ ∊ (e_{j t} \cdot e_{k t} + \frac{∊}{\tilde{∊}} e_{j z} \cdot e_{k z})

The term e_jt · e_kt means e_jx · e_kx + e_jy · e_ky. The factors N_j,

N_{j} = \iint d A \hat{z} \cdot (e_{j t} \times h_{j t}),

are normalization factors; for real fields, N_j is power.

In the traditional application of CCMT, e_k and h_k are eigenmodes of the system of interest (e.g., modes of an optical fiber). The spatial evolution or transfer of power between modes in the presence of the permittivity perturbation Δ∊ (e.g., a bend in the fiber) is then given by Eq. (3).

2.2. CCMT eigenmode solutions

In the CCMT solver, the coupled permittivity ∊̃ is given by the refractive index map n(ρ) or n(x, y) of the waveguide being solved, ∊̃ = ∊₀n². We consider the basis in vacuum, such that ∊ = ∊₀ (the vacuum permittivity). In order to compute the waveguide eigenmodes using CCMT, we impose harmonic phase dependence in the CCMT coefficients along the propagation axis z:

a_{k} (z) = a_{k} e^{- i λ z} and b_{k} (z) = b_{k} e^{- i λ z} .

The wavevector λ gives the effective index of a given mode via λ = n_effk, where k = c/ω is the wavevector of the field. Using the Eq. (6), the CCMT equations (Eq. (3)) are algebraically rearranged to an eigenvalue equation,

(\begin{matrix} β + κ_{N} & χ_{N} \\ - χ_{N} & - β - κ_{N} \end{matrix}) (\begin{matrix} a \\ b \end{matrix}) = λ (\begin{matrix} a \\ b \end{matrix}) .

This is our main equation. In this matrix equation, β is a vector constructed out of the entries β_j (11), and χ_N and κ_N are matrices constructed out of the elements of κ_jk and χ_jk (Eq. (4)); each row j has been normalized (divided) by N_j as indicated by the subscript N.

Numerical diagonalization of this matrix yields eigenvectors, which give the CCMT coefficients a_k and b_k, and eigenvalues, which give the effective index for each coupled mode. The vector fields of coupled eigenmodes are obtained by substituting these coefficients into Eq. (1) and (2); the squared electric field magnitude I of a given CCMT eigenmode is given by

I \propto \sum_{t = (x, y)} {| \sum_{k} (a_{k} (z) + b_{k} (z)) e_{k t} (x, y) |}^{2} + {| \sum_{k} (a_{k} (z) - b_{k} (z)) \frac{∊}{\tilde{∊}} e_{k z} (x, y) |}^{2} .

3. Eigensolver implementation details

In this section, the finite difference eigenmode (FDE) method and analytic expressions describing a step-index fiber and a plasmonic wire will be introduced, as these will be useful in validating results obtained using CCMT. Two methods of computing the CCMT overlap integrals are discussed. Details related to our GPU implementation of the CCMT mode solver are given.

3.1. Basis modes

In the CCMT mode solver we use arbitrary basis modes e_k and h_k, which are not eigenmodes of the system being solved. We use the electromagnetic modes of a dielectric-filled metal pipe waveguide [40]. This basis set is complete and orthogonal; the circular geometry is also convenient for many applications when Cartesian coordinates are not ideal. The TE_nm modes are

\begin{matrix} e_{ρ} = - \frac{i ω μ n}{k_{c}^{2} ρ} (A cos n ϕ - B sin n ϕ) J_{n} (k_{c} ρ), e_{ϕ} = \frac{i ω μ}{k_{c}} (A sin n ϕ - B cos n ϕ) J_{n}^{'} (k_{c} ρ), \\ h_{ρ} = - \frac{i β}{k_{c}} (A sin n ϕ + B cos n ϕ) J_{n}^{'} (k_{c} ρ), h_{ϕ} = \frac{- i β n}{k_{c}^{2} ρ} (A cos n ϕ - B sin n ϕ) J_{n} (k_{c} ρ), \\ e_{z} = 0, and, h_{z} = (A sin n ϕ - B cos n ϕ) J_{n} (k_{c} ρ), with k_{c} = \frac{p_{n m}^{'}}{a} . \end{matrix}

The TM_nm modes are

\begin{matrix} e_{ρ} = - \frac{i β}{k_{c}} (A sin n ϕ + B cos n ϕ) J_{n} (k_{c} ρ), e_{ϕ} = \frac{- i β_{n}}{k_{c}^{2} ρ} (A cos n ϕ - B sin n ϕ) J_{n} (k_{c} ρ), \\ h_{ρ} = \frac{i ω ∊ n}{k_{c}^{2} ρ} (A cos n ϕ - B sin n ϕ) J_{n} (k_{c} ρ), h_{ϕ} = \frac{- i ω ∊}{k_{c}} (A sin n ϕ + B cos n ϕ) J_{n}^{'} (k_{c} ρ), \\ e_{z} = (A sin n ϕ + B cos n ϕ) J_{n} (k_{c} ρ), and, H_{z} = 0, with k_{c} = \frac{p_{n m}}{a} . \end{matrix}

The vector components are given in cylindrical coordinates; ρ is the radial coordinate, ϕ the angular one. The expressions are given at z = 0. In these expressions, μ = μ₀ is the vacuum magnetic permeability and

β = \sqrt{k^{2} - k_{c}^{2}}

is the modal wavevector. The value k_c is called the cutoff wavevector. Modes where k_c > k have imaginary wavevectors. They decay rather than propagate in a real waveguide, but are important in the CCMT mode solver as they can describe higher spatial frequency components. The values n ≥ 0 and m ≥ 1 are integers enumerating the set of basis modes. The functions J_n are Bessel functions of the first kind of order n; the functions J′_n are the derivatives of the Bessel functions of the first kind, which can be computed as

J_{n}^{'} (ξ) = \frac{1}{2} (J_{n - 1} (ξ) - J_{n + 1} (ξ)) = - J_{n + 1} (ξ) + \frac{n}{ξ} J_{n} (ξ) .

The latter form [5] is convenient in our implementation (Sec. 3.3). The values p_nm and p′_nm give the zeros of J_n and J′_n (the values m index the set of Bessel function zeros for a given order n). The radius a of the pipe waveguide must contain the refractive index map of the waveguide being studied as well as any penetrating field appearing in the computed eigenmodes.

The input basis set for CCMT is constructed from these values by building TE_nm or TM_nm modes with either A = 1 and B = 0, or A = 0 and B = 1. This basis set is sorted by increasing k_c, which is real; once sorted, these are the modes e_j and h_j (associated with wavevectors β_j) that are input into CCMT.

3.2. Two implementations

We present two distinct CCMT eigenmode solver implementations that differ in the way they evaluate the overlap integrals (Eq. (4)). We call them the radial solver and the Cartesian solver.

3.2.1. The radial solver

In systems where there is no angular dependence in ∊̃, ∊̃ and Δ∊ are functions of the radial variable ρ only. Consider the integral in Eq. (4a),

\int d A Δ ∊ (ρ) (e_{j t} \cdot e_{k t} - \frac{∊}{\tilde{∊}} e_{j z} \cdot e_{k z}) \equiv \int d A Δ ∊ (ρ) I .

The transverse terms in I are I₁ + I₂ = e_jx · e_kx + e_jy · e_ky. But the transformation mapping the transverse part of the basis modes from polar coordinates to Cartesian coordinates is

(\begin{matrix} e_{x} \\ e_{y} \end{matrix}) = (\begin{matrix} cos ϕ & - sin ϕ \\ sin ϕ & cos ϕ \end{matrix}) (\begin{matrix} e_{ρ} \\ e_{ϕ} \end{matrix}) .

Ones sees that e_jx · e_kx + e_jy · e_ky = e_j_ρ · e_k_ρ + e_jϕ · e_kϕ. In this way, I₁ and I₂ are each separable into products of functions, R_i(ρ)Φ_i(ϕ). The longitudinal or z term in I, I₃, is also separable into a product of functions. This means that the integral over each I_i can be evaluated as a product two 1D integrals,

\int d A Δ ∊ (ρ) I_{i} = \int d ρ ρ Δ ∊ (ρ) R_{i} (ρ) \int d ϕ Φ_{i} (ϕ) .

The angular integrations over ϕ can be done analytically. The radial integrations are best evaluated numerically (e.g. direct lattice summation) as they involve products of Bessel functions of mixed order and argument. Similar considerations (using Eq. (14)) for N_j (Eq. (5)) reveal that it is also a separable integration, and that the cylindrical basis components can also be used directly in place of the Cartesian ones in this expression. This method substantially reduces the computational cost of the integrations, as summation is over a 1D rather than a 2D lattice; as a consequence, a higher resolution lattice can be employed for the 1D integration.

3.2.2. The Cartesian solver

Where ∊̃ is a function of both ρ and ϕ (or x and y), the above method cannot be employed, and the integrations in Eq. (4) will be done on a 2D Cartesian lattice. Eq. (14) is used to rotate the basis function components specified in cylindrical coordinates to a x and y components on Cartesian lattice. In both the radial solver and this Cartesian solver, the final presentation of eigenmodes on a 2D lattice requires a Cartesian representation of the basis eigenmodes to be computed.

We have aligned the coordinate systems such that ρ is never zero. For example, in 1D, we align the lattice such that the values of ρ being sampled around the origin are −0.5 and 0.5, rather than −1, 0. Of course, there are more sophisticated methods of integration around singularities [43].

3.3. Parallelization

With N basis modes, CCMT requires the computation of N normalization factors N_j and approximately 3N²/2 integrations for the terms in κ_jk and χ_jk (there are three distinct terms in the integrands of these two integrals, and the coefficient matrices are symmetric). These quantities can be computed in parallel on a graphics processing unit (GPU) capable of general-purpose computation. The point here is not to perform detailed benchmarking nor to do extensive code optimization, but to simply report that mathematical form of CCMT form makes effective use of hardware capable of doing highly parallelized computation, and that this is a boon to this method.

We use the Nvidia CUDA parallel computing platform; CUDA libraries also provide built-in support for evaluation of Bessel functions on the GPU. C language CUDA kernels are called directly from MATLAB. The kernel that renders the basis modes is straightforward, computing individual pixels of all basis modes in parallel. For the overlap integrals, each kernel invocation (at a particular j and k) computes a complete numerical integral. Parallel reduction methods are not needed as the distinct sums run in parallel. The second form of the Bessel function derivative given in Eq. (12) is useful as the CUDA Bessel implementation does not allow the computation of Bessel functions of negative order.

In the GPU implementation, the construction of the cache of CCMT basis modes and the evaluation of all the CCMT integrations is very rapid–about an order of magnitude faster than our unoptimized FDE implementation.

4. Comparison methods

In this section we briefly review methods with which we will compare the CCMT solutions with other methods.

4.1. Finite-difference eigenmode solutions

An effective method for solving waveguide modes is to discretize Maxwell’s partial differential equations and use finite differences in place of derivatives, and then algebraically rearrange the problem to an eigenvalue equation. This yields a set of eigenvectors (the field distributions) and eigenvalues (the effective indexes) describing the modes. Mode solvers based on this technique are called finite difference eigenmode (FDE) solvers.

We implemented a fully vectorial FDE mode solver as specified in the literature [8], and use it, where possible, to make comparisons with the fully vectorial CCMT eigenmodes that we compute. In some cases, however, it is desirable or necessary to use conformal mesh technology [44–46] to “smooth” curved material interfaces that otherwise appear jagged on a Cartesian lattice. In these cases we use a commercial mode solver [47] whose base implementation of FDE is otherwise similar to ours.

4.2. Analytic results for a step-index fiber

The modes and effective indexes of a step-index fiber can be computed numerically using analytic expressions [5]. The transverse electric field in the core is

e_{ρ} = i (- \frac{A β}{u} J_{ν}^{'} (u ρ) + B \frac{μ ω}{u^{2}} \frac{ν}{ρ} J_{ν} (u ρ)) sin ν ϕ, e_{ϕ} = i (- \frac{A β}{u^{2}} \frac{ν}{ρ} J_{ν} (u ρ) + B \frac{μ ω}{u} J_{ν}^{'} (u ρ)) cos ν ϕ .

In the cladding it is:

e_{ρ} = i (\frac{C β}{w} K_{ν}^{'} (w ρ) - D \frac{μ ω}{w^{2}} \frac{ν}{ρ} K_{ν} (w ρ)) sin ν ϕ, e_{ϕ} = i (\frac{C β}{w^{2}} \frac{ν}{ρ} K_{ν} (u r) - D \frac{μ ω}{w} K_{ν}^{'} (w ρ)) cos ν ϕ .

Expressions giving the other components will not be given here. The functions K are modified Bessel functions of the second kind; the primed versions are their derivatives. The modal wavevector β appears in these expressions directly as well as in the quantities

u^{2} = k_{1}^{2} - β^{2}

and

w^{2} = β^{2} - k_{2}^{2}

, with k_i = n_ik. The values u and k₁ are in the core; w and k₂ are in the cladding. Boundary conditions in terms of the longitudinal and the ϕ part of the fields give the coefficients A, B, C and D. The characteristic equation of the 4×4 matrix associated with these four boundary condition equations and coefficients can be written as [5]

(\frac{J_{ν}^{'} (u b)}{u J_{ν} (u b)} + \frac{K_{ν}^{'} (u b)}{w K_{ν} (u b)}) (k_{1}^{2} \frac{J_{ν}^{'} (u b)}{u J_{ν} (u b)} + k_{2}^{2} \frac{K_{ν}^{'} (w b)}{w K_{ν} (w b)}) = {(\frac{β ν}{b})}^{2} {(\frac{1}{u^{2}} + \frac{1}{w^{2}})}^{2} .

This characteristic equation can be solved numerically, yielding the set of modal wavevectors β. The quantity b is the radius of the fiber.

4.3. Analytic results for a Sommerfeld modes

A cylindrical metal wire functions as a plasmonic waveguide; the propagating mode is a Sommerfeld surface wave. The propagation constant β is found by numerically solving the following expression [48,49]:

\frac{K_{0} (p_{d} b) I_{1} (p_{m} b)}{K_{1} (p_{d} b) I_{0} (p_{m} b)} = - \frac{∊_{d} p_{m}}{∊_{m} p_{d}} .

The functions I are modified Bessel functions of the first kind. The quantities p_m,d are given by $p_{m, d} = \sqrt{β^{2} - k^{2} ∊_{m, d}}$ . The wire radius is b, ∊_d = ∊ is the permittivity of the surrounding dielectric, and ∊_m is the permittivity of the metal.

5. Results

To demonstrate the CCMT mode solver, we find the eigenmodes of three dielectric waveguides, shown in Fig. 1. The first two, Figs. 1(a) and 1(b), have no angular dependence in ∊̃ and can be solved using the radial solver methods of Sec. 3.2.1. We also compute the Sommerfeld mode of a metal nanowire, using the radial solver with a modified basis.

Fig. 1 Refractive maps n(ρ) (in (a)–(b)) and n(x, y) (in (c)) of three dielectric waveguides used to illustrate the coupled mode theory (CCMT) eigenmode solver. Each figure gives the refractive index in the transverse plane. The specified n_f is the “fiber” index in the black region; the white background has an index of 1.0.

Download Full Size | PDF

5.1. A step-index fiber

The step-index fiber used here is shown in Fig. 1(a). The core has refractive index n_f = 1.6 and the cladding is air; the fiber core is 4.2 μm in radius. Modes are computed using the radial solver, as the refractive index can be specified as a function of ρ only. Step-index fiber modes are computed for radiation with vacuum wavelength λ = 1.5 μm; N = 600 basis modes are used. With this choice of λ and N the basis is a mixture of modes with real and imaginary β. A radial lattice with 7000 points is used; after the integration and diagonalization steps, the 2D eigenmode intensity is constructed on a 215 × 215 lattice.

Figure 2 shows the CCMT eigenmodes. With this size of basis set, good agreement is observed between the CCMT eigenmodes and those calculated using either FDE (Sec. 4.1) or analytic expressions (Sec 4.2) for dozens of modes starting at the fundamental. Table 1 summarizes the effective indexes n_eff computed using CCMT, FDE, and analytic expressions. By including a comparison with the analytic expressions (Eqs. (16), (17) and (18)) we see how well CCMT and FDE compare with theory, giving some calibration of the meaning of the size of errors between FDE and CCMT. For this and the other example systems, we have not shown the eigenmode field intensity from these comparison methods as they are visually identical to the CCMT eigenmodes.

Fig. 2 Electric field intensity of the eigenmodes of a step-index fiber (given in Fig. 1(a)) calculated using CCMT. The wavelength is λ = 1.5 μm, and the basis set size is N = 600. The corresponding effective indexes n_eff obtained using FDE are 1.5945 and 1.5750; using analytic expressions, they are 1.5947 and 1.5760.

Download Full Size | PDF

Table 1. Computed effective indexes for the step-index fiber modes shown in Fig. 2.

View Table | View all tables in this article

The convergence behavior of the solver for this step-index fiber structure is illustrated in Fig. 3. Timings showing the total computation time spent as a function of the number of basis modes, N. These timings were taken on a system built with an Intel Core i7–7700HQ CPU and a Nvidia GeForce GTX 1070 GPU. The discontinuity in the slope, between N = 300 and N = 400, arises as the propagation constant β of the basis modes becomes imaginary past cutoff. The introduction of imaginary modes causes the CCMT overlap matrix (Eq. (7)) to transition from purely real to complex, and the (CPU-based) numerical diagonalization routine runs slower with complex matrices. These timings reflect finding the set of all eigenvectors and eigenvalues (of the CCMT matrix) at each point.

Fig. 3 Convergence and timings for the step-index fiber (given in Fig. 1(a)). In (a), the dependence of the effective index of the fundamental mode, n_eff, on the number of basis modes, N, is given. In (b) the times spent computing the eigenmode solutions in (a) are shown, as a function of N. The sample size for each data point is ten. The reason for the kink at N = 300 is discussed in the text.

Download Full Size | PDF

5.2. A lossy dielectric shell

Figure 1(b) shows another system with no angular dependence, a lossy dielectric ring. The shell has refractive index n_f = 1.6 + 0.2i. The inner radius is 3.6 μm and the outer radius is 4.8 μm. Modes are computed using the radial solver; a vacuum wavelength of λ = 2.5 μm is used. We used N = 1300 basis modes. With this choice of λ and N the majority of basis modes are have imaginary β. The same size of radial and 2D eigenmode lattices are used as were used for the step-index fiber.

Figure 4 shows the CCMT eigenmodes. In this thin shell, the staircasing of the circular fiber boundaries led to spurious features in the modes obtained using our own FDE solver (or any solver where staircasing is present). As such, we use a commercial FDE solver here to compute modes using conformal mesh technology (using the method labelled Conformal Variant 0 in the software) (Sec. 4.1). Table 2 summarizes the effective indexes n_eff computed using CCMT and FDE (with a conformal mesh). The modal fields agree well when the FDE conformal mesh is used.

Fig. 4 Electric field intensity of the eigenmodes of a dielectric shell waveguide (given in Fig. 1(b)) calculated using CCMT. The wavelength is λ = 2.5 μm, and the basis set size is N = 1300. The corresponding effective indexes n_eff obtained using FDE (with a conformal mesh) are 1.4473+0.1970i, 1.4447+0.1972i and 1.4447+0.1972i.

Download Full Size | PDF

Table 2. Computed effective indexes for the lossy dielectric shell fiber modes shown in Fig. 4.

View Table | View all tables in this article

5.3. A microstructured fiber

This third structure, Fig. 1(c), is solved on a Cartesian lattice using the methods of Sec. 3.2.2; this is our only example where the radial solver is not employed. The eigenmodes are shown in Fig. 5. The fiber material is of index n_f = 1.5. The outside radius of the fiber is about 5.0 μm. These modes are solved using λ = 800 nm, and N = 1200 basis modes. With this choice of λ and N all of the basis modes have a real propagation constant β. The modes are computed using a 215 × 215 lattice. Table 3 summarizes the effective indexes n_eff computed using CCMT and FDE. Good agreement is obtained, even to high mode numbers (e.g. several dozens of modes), with the modes found using FDE.

Fig. 5 Electric field intensity of the eigenmodes of a microstructured dielectric waveguide (given in Fig. 1(c)) calculated using CCMT. The wavelength is λ = 800 nm, and the basis set size is N = 1200. The corresponding effective indexes n_eff obtained using FDE are 1.4940, 1.4749 and 1.4621.

Download Full Size | PDF

Table 3. Computed effective indexes for the microstructured optical fiber modes shown in Fig. 5.

View Table | View all tables in this article

5.4. Nanowire

We also test the mode solver with a metal nanowire as CCMT can handle plasmonic materials. These metal wire modes were introduced in Sec. 4.3 along with a method for analytically computing the propagation constant β. The nanowire has a radius of 6.0 nm and dielectric permittivity ∊_m = −12.95 + 1.12i. The wavelength is λ = 650 nm, and the basis set size is N = 500. In order to find this mode, we have restricted the metal pipe basis modes to only those with n = 0, matching the expected radial intensity distribution of the Sommerfeld mode. The n > 0 modes do not contribute to the solution and decrease performance. The modes have been computed using a radial lattice with 4000 points. With this choice of λ and N all of the basis modes have a real propagation constant β.

Fig. 6 Electric field intensity of the Sommerfeld mode in a gold nanowire 6 nm in radius, calculated using CCMT. The wavelength is λ = 650 nm, and the basis set size is N = 550 (this is a modified basis set; see the text for details). The effective index n_eff calculated using CCMT, using the parameters given here and in the text, is 5.93+0.35i. For comparison, the effective index n_eff found using FDE (using a conformal mesh) is 5.92+0.33i and the analytic result is 5.81+0.34i.

Download Full Size | PDF

Table 4. Computed effective indexes for the plasmonic nanowire shown in Fig. 6.

View Table | View all tables in this article

Figure 6 shows the computed Sommerfeld mode. Sinusoidal artifacts can be seen in the eigenmode solution, a manifestation of Gibbs phenomenon. This arises when recreating the sharp edge at the metal-dielectric boundary with a basis of Bessel functions. Otherwise, good agreement is seen between the effective index and modal field computed using CCMT and those found using FDE and analytic results. We again use a commercial FDE solver to compute the modes using conformal mesh technology (Sec. 4.1); in this case, a particular type of conformal mesh well suited to metal-dielectric interfaces is used (labelled Conformal Variant 2). Table 4 summarizes the effective indexes n_eff computed using CCMT, FDE, and the analytic result.

6. Discussion

As seen, CCMT can generate modal solutions for waveguides that agree well with analytic solutions (when available) and those computed using FDE. In this section we discuss a few of the strengths and weaknesses of the CCMT eigenmode solver.

One benefit of the CCMT solver is that it allows a choice of a convenient coordinate system when evaluating the overlap integrals (Eqs. 4). This was discussed in Sec. 3.2. Closely related to the choice of coordinate system is the choice of basis. CCMT allows you to choose any convenient basis. This could be a basis of fiber LP modes, modes of a metal pipe or box (as used here), or any other basis that fits well with the geometry of the problem.

In the radial solver examples, circular coordinates allowed us to separate the area integration into two 1D integrations, one of which could be done analytically. Other systems were best evaluated on a 2D Cartesian lattice. A third solver implementation where a basis of rectangular waveguide modes was used was also tested; in this case, all area integrals could be done as separate analytic integrations along x and y (when, e.g., the structure of interest can be expressed in terms of a small number of rectangular parts). The metal pipe modes used here allow relatively quick convergence is many cases because they match well with the symmetry of the waveguides studied; a basis of rectangular modes converges more slowly, despite the advantages of analytic integration. The use of radial coordinates in the dielectric shell proved valuable in obtaining accurate solutions, as staircasing artifacts along the circular boundary would have been introduced if this were done on a 2D lattice.

Being based on complex CCMT, the method allows you to handle lossy materials. This was seen in Sec. 5.2.

Another benefit of the CCMT mode solver is that the computational cost of the initial overlap integrations and basis set construction can be reduced or eliminated. It can be reduced if one desires to change only a small area of the structure and then re-compute modes—with localized changes to n(ρ) or n(x, y)), only a small fraction of the original integration lattice needs to be re-evaluated. Overlap integrations can be completely eliminated (except for the initial point) when computing dispersion curves. This is possible as all of the wavelength dependence in the chosen basis (Sec. 3.1) is contained in the constants β and ω. One factors these out, and then performs the overlaps integrations only once. (Dispersive materials where the wavelength cannot be factored out will not work.) Multiplying CCMT matrix entries by the appropriate frequency-dependent terms is then sufficient to introduce the wavelength dependence before each diagonalization. This is most beneficial if the basis set is small. If the basis is small, in this case the CCMT coefficient matrix (Eq. (7)) is relatively quick to diagonalize. (If iterative methods are used to find a subset of the eigenvalues of the CCMT coefficient matrix, one can also seed each successive eigenvalue computation with a value from previous iterations.)

The CCMT mode solver may be beneficial in the waveguide mode matching method [40]. Ordinarily, mode overlap integrals must be performed on each side of the interface, leading to a large number of integrations. When the modes have been solved using CCMT, however, modes on each side of the interface are always the same: they are the basis modes e_k and h_k weighted by the coefficients a and b. Thus only a single computation of the overlaps between the basis modes is needed, for all wavelengths and structures. This is analogous to the advantages allowed for by Fourier modal method, where a common Fourier basis throughout the structure allows for efficient matching at the boundaries [29]. Indeed, there have been several works dedicated to mode matching with different bases [31–33] that provide more efficient computation by selecting a convenient basis set for matching or expansion. The CCMT present here provides a general framework for representing a common basis for efficient matching at boundaries including losses and perfectly matched layers, and so it is a good candidate for efficient mode matching representations that may be investigated in future works.

A possible drawback of the method is the computational cost of matrix diagonalization, if the number of basis modes grows large. That said, the small dense matrix diagonalizations performed here were substantially faster (e.g. an order of magnitude) than the large sparse matrix diagonalizations used in FDE (comparing, e.g., a 215 × 215 2D lattice CCMT solution to a similarly sized 2D FDE solution in MATLAB). Numerical stability criterion and further convergence studies are important issues for follow-up studies.

7. Conclusion

We have shown that the CCMT mode solver is able to generate accurate modal solutions of Maxwell’s equations, and have compared them to analytic and other numerical results. As this work fits into a vast body of literature related to modal solutions of Maxwell’s equations in optical waveguides, we have given a brief discussion of some of the possible strengths of the method. We have also attempted to highlight its suitability to GPU implementation, as the algorithm is highly parallelizable in nature.

Appendix

In the examples of Sec. 5 we have avoided the use of a perfectly matched layer (PML) or other coordinate stretching or absorbing boundary condition by considering structures with appreciable refractive index contrast. In such structures, the modes are strongly guided; little field amplitude exists near the simulation boundaries. In the examples above, we verified this by using a FDE solver with PML support and comparing the solutions obtained with PML turned on vs. the results obtained with it turned off.

As our solver is based on complex coupled mode theory, it should be able to accommodate perfectly matched layer (PML) boundaries. As a demonstration, we find eigenmode solutions of the 1D slab waveguide system studied in [38]. This waveguide, shown in Fig. 7(a), consists of an core of index n_core = 1.458, an inner cladding with index n_cl = 1.450. The surrounding “substrate” index n_s is varied to produce either a strongly guiding or a quasi-leaky waveguide. PML layers 2.5 μm thick, on both sides, are computed using parameters as defined in [38]. Figure 7 shows the guided, quasi-leaky and PML modes thus obtained using the CCMT method. It is seen that when the refractive index contrast between the core/cladding and the substrate is weak, the obtained mode is quasi-leaky, and exhibits the expected attenuation near the boundary due to the presence of the PML.

Fig. 7 A 1D slab waveguide system illustrating the use of a perfectly matched layer (PML) in the CCMT mode solver. In (a) an example structure is given; the details are given in the text. In (b)–(d) we see the field profile for the fundamental guided mode, the lowest order quasi-leaky mode, and the lowest order PML mode as solved using CCMT. The interfaces between the core, inner cladding, substrate and PML regions are marked with vertical lines.

Download Full Size | PDF

Funding

This work has been supported by an NSERC Discovery Grant.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References and links

1. P. Russell, “Photonic Crystal Fibers,” Science 299(5605), 358–362 (2003). [CrossRef] [PubMed]

2. P. Russell, “Photonic-Crystal Fibers,” J. Lightwave Technol. 24(12), 4729–4749 (2006). [CrossRef]

3. B. R. West and A. S. Helmy, “Properties of the quarter-wave Bragg reflection waveguide: theory,” J. Opt. Soc. Am. B 23(6), 1207–1220 (2006). [CrossRef]

4. S. A. Maier and H. A. Atwater, “Plasmonics: Localization and guiding of electromagnetic energy in metal/dielectric structures,” J. Appl. Phys. 98(011101), 1–9 (2005). [CrossRef]

5. G. Keiser, Optical Fiber Communications, 3rd ed. (McGraw-Hill, 2000).

6. E. Schweig and W.B. Bridges, “Computer Analysis of Dielectric Waveguides: A Finite-Difference Method,” IEEE Trans. Microw. Theory Techn. 32(5), 531–541 (1984). [CrossRef]

7. K. Bierwirth, N. Shulz, and F. Arndt, “Finite-Difference Analysis of Rectangular Dielectric Waveguide Structures,” IEEE Trans. Microw. Theory Techn. 34(11) 1104–1114 (1986). [CrossRef]

8. Z. Zhu and T. G. Brown, “Full-vectorial finite-difference analysis of microstructured optical fibers,” Opt. Express 10(17), 853–864 (2002). [CrossRef] [PubMed]

9. C.-P. Yu and H.-C. Chang, “Yee-mesh-based finite difference eigenmode solver with PML absorbing boundary conditions for optical waveguides and photonic crystal fibers,” Opt. Express 12(25), 6165–6177 (2004). [CrossRef] [PubMed]

10. Y.-C. Chiang, Y.-P. Chiou, and H.-C. Chang, “Improved Full-Vectorial Finite-Difference Mode Solver for Optical Waveguides With Step-Index Profiles,” J. Lightwave Technol. 20(8), 1609 (2002). [CrossRef]

11. C.-P. Yu and H.-C. Chang, “Applications of the finite difference mode solution method to photonic crystal structures,” Opt. Quant. Electron. 36(1–3), 145–163 (2004). [CrossRef]

12. A. B. Fallahkhair, K. S. Li, and T. E. Murphy, “Vector Finite Difference Modesolver for Anisotropic Dielectric Waveguides,” J. Lightw. Technol. 26(11), 1423–1431 (2008). [CrossRef]

13. Y.-C. Lu, L. Yang, W.-P. Huang, and S.-S. Jian, “Improved Full-Vector Finite-Difference Complex Mode Solver for Optical Waveguides of Circular Symmetry,” J. Lightwave Technol. 26(13), 1868–1876 (2008) [CrossRef]

14. J. A. M. Svedin, “A Numerically Efficient Finite-Element Formulation for the General Waveguide Problem Without Spurious Modes,” IEEE Trans. Microw. Theory Techn. 37(11), 1708–1715 (1989). [CrossRef]

15. S. S. A. Obayya, B. M. A. Rahman, and H. A. El-Mikati, “New Full-Vectorial Numerically Efficient Propagation Algorithm Based on the Finite Element Method,” J. Lightwave Technol. 18(3), 409 (2000) [CrossRef]

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

17. A. Cucinotta, S. Selleri, L. Vincetti, and M. Zoboli, “Holey fiber analysis through the finite-element method,” IEEE Photon. Technol. Lett. 14(11), 1530–1532 (2002). [CrossRef]

18. S. S. A. Obayya, B. M. Azizur 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(6), 1054 (2002). [CrossRef]

19. J. E. Goell, “A Circular-Harmonic Computer Analysis Of Rectangular Dielectric Waveguide,” Bell Syst. Tech. J. 48(7), 2133–2160 (1969). [CrossRef]

20. Y.-H. Wang and C. Vassallo, “Circular Fourier analysis of arbitrarily shaped optical fibers,” Opt. Lett. 14(24), 1377–1379 (1989). [CrossRef] [PubMed]

21. A. Khavasi, K. Mehrany, and B. Rashidian, “Three-dimensional diffraction analysis of gratings based on Legendre expansion of electromagnetic fields,” J. Opt. Soc. Am. B 24(10), 2676–2685 (2007). [CrossRef]

22. J. P. Hugonin, P. Lalanne, I. D. Villar, and I. R. Matias, “Fourier modal methods for modeling optical dielectric waveguides,” Opt. Quant. Electron. 37(1–3), 107–119 (2005). [CrossRef]

23. G. Colas des Francs, J.-P. Hugonin, and J. Čtyroký, “Mode solvers for very thin long–range plasmonic waveguides,” Opt. Quant. Electron. 42(8), 557–570 (2011). [CrossRef]

24. J. Xiao and X. Sun, “Full-vectorial mode solver for anisotropic optical waveguides using multidomain spectral collocation method,” Opt. Commun. 283(14), 2835–2840 (2010). [CrossRef]

25. A. Ferrando, E. Silvestre, J. J. Miret, P. Andrés, and M. V. Andrés, “Full-vector analysis of a realistic photonic crystal fiber,” Opt. Lett. 24(5), 276–278 (1999). [CrossRef]

26. Z. Zhu and T. G. Brown, “Multipole analysis of hole-assisted optical fibers,” Opt. Commun. 206(4), 333–339 (2002). [CrossRef]

27. E. Silvestre, M. V. Andres, and P. Andres, “Biorthonormal-Basis Method for the Vector Description of Optical-Fiber Modes,” J. Lightwave Technol. 16(5), 923 (1998). [CrossRef]

28. A. Weisshaar, J. Li, R. L. Gallawa, and I. C. Goyal, “Vector and quasi-vector solutions for optical waveguide modes using efficient Galerkin’s method with Hermite-Gauss basis functions,”J. Lightwave Tech. 13(8), 1795–1800 (1995). [CrossRef]

29. P. Lalanne and E. Silberstein, “Fourier-modal methods applied to waveguide computational problems,” Opt. Lett. 25(15), 1092–1094 (2000). [CrossRef]

30. M. G. Moharam, E. B. Grann, D. A. Pommet, and T. K. Gaylord, “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” J. Opt. Soc. Am. A 12(5), 1068–1076 (1995). [CrossRef]

31. G. Sztefka and H. P. Nolting, “Bidirectional eigenmode propagation for large refractive index steps,” IEEE Photon. Technol. Lett. 5(5), 554–557 (1993). [CrossRef]

32. A. S. Sudbo, “Film mode matching: a versatile numerical method for vector mode field calculations in dielectric waveguides,” Pure Appl. Opt. 2(3), 211–233 (1993). [CrossRef]

33. P. Bienstman, CAMFR full-vectorial Maxwell solver. http://camfr.sourceforge.net

34. K. S. Chiang, “Review of numerical and approximate methods for the modal analysis of general optical dielectric waveguides,” Opt. Quant. Electron. 26(3), S113–S134 (1994). [CrossRef]

35. C. Vassallo, “1993—1995 Optical mode solvers,” Opt. Quant. Electron. 29(2), 95–114 (1997). [CrossRef]

36. J.-R. Qian and W.-P. Huang, “Coupled-mode theory for LP modes,” J. Lightwave Technol. 4(6), 619–625 (1986). [CrossRef]

37. Z. H. Wang, “Application of the Coupled Mode Theory to Eigenvalue Problems of Graded-Index Optical Fibers,” Microw. Opt. Technol. Lett. 12(2), 90–93 (1996). [CrossRef]

38. W.-P. Huang and J. Mu, “Complex coupled-mode theory for optical waveguides,” Opt. Express 17(21), 19134–19152 (2009). [CrossRef]

39. A. Yariv, “Coupled-mode theory for guided-wave optics,” IEEE J. Quant. Electron. 9(9), 919–933 (1973). [CrossRef]

40. D. M. Pozar, Microwave Engineering, 4th ed. (John Wiley & Sons, 2012).

41. C. Yeh, “Elliptical Dielectric Waveguides,” J. Appl. Phys. 33(11), 3235–3243 (1962). [CrossRef]

42. D. A. Goldberg, L. J. Laslett, and R. A. Rimmer, “Modes of Elliptical Waveguides: A Correction,” IEEE Trans. Microw. Theory Techn. 38(11), 1603–1608 (1990). [CrossRef]

43. S. Amari and J. Bornemann, “Efficient numerical computation of singular integrals with applications to electromagnetics,” IEEE Trans. Antennas Propag. 43(11), 1343–1348 (1995). [CrossRef]

44. W. Yu and R. Mittra, “A conformal finite difference time domain technique for modeling curved dielectric surfaces,” IEEE Microw. Wirel. Compon. Lett. 11(1), 25–27 (2001). [CrossRef]

45. A. Mohammadi, H. Nadgaran, and M. Agio, “Contour-path effective permittivities for the two-dimensional finite-difference time-domain method,” Opt. Express 13(25), 10367–10381 (2005). [CrossRef] [PubMed]

46. A. Farjadpour, D. Roundy, A. Rodriguez, M. Ibanescu, P. Bermel, J. D. Joannopoulos, S. G. Johnson, and G. W. Burr, “Improving accuracy by subpixel smoothing in the finite-difference time domain,” Opt. Lett. 31(20), 2972–2974 (2006). [CrossRef] [PubMed]

47. Lumerical Solutions, Inc.http://www.lumerical.com/tcad-products/mode/

48. A. Sommerfeld, “Mathematische theorie der diffraction,” Math. Ann. 47(2), 317–374 (1896). [CrossRef]

49. R. Gordon, “Reflection of Cylindrical Surface Waves,” Opt. Express 17(21), 18621–18629 (2009). [CrossRef]

Complex coupled mode theory electromagnetic mode solver

Abstract

1. Introduction

2. Formulation of complex coupled mode theory

2.1. Complex coupled mode theory

2.2. CCMT eigenmode solutions

3. Eigensolver implementation details

3.1. Basis modes

3.2. Two implementations

3.2.1. The radial solver

3.2.2. The Cartesian solver

3.3. Parallelization

4. Comparison methods

4.1. Finite-difference eigenmode solutions

4.2. Analytic results for a step-index fiber

4.3. Analytic results for a Sommerfeld modes

5. Results

5.1. A step-index fiber

5.2. A lossy dielectric shell

5.3. A microstructured fiber

5.4. Nanowire

6. Discussion

7. Conclusion

Appendix

Funding

Disclosures

References and links

Cited By

Figures (7)

Tables (4)

Equations (23)

Optics Express

Method	Effective Index, n_eff
Method	Mode 1	Mode 7
CCMT	1.5944	1.5748
FDE	1.5945	1.5750
Analytic	1.5947	1.5760

Method	Effective Index, n_eff
Method	Mode 1	Mode 2	Mode 3
CCMT	1.4474+0.1971i	1.4448+0.1973i	1.4448+0.1973i
FDE	1.4473+0.1970i	1.4447+0.1972i	1.4447+0.1972i

Method	Effective Index, n_eff
Method	Mode 1	Mode 19	Mode 39
CCMT	1.4938	1.4745	1.4617
FDE	1.4940	1.4749	1.4621