Complex coupled mode theory electromagnetic mode solver

We present a method for computing waveguide eigenmodes based on complex coupled mode theory (CCMT). This approach generalizes Fourier transform methods by allowing an arbitrary but convenient basis set to be used for the expansion. In the presented approach, one is free to choose an arbitrary basis representation; for example, we show the use of electromagnetic modes of a cylindrical metal waveguide. CCMT-computed modes are compared with modes computed using analytic expressions and results obtained using a finite difference solver. In cases where the basis set is small, the method can efficiently re-compute modes after structural refinements are made, and can efficiently compute dispersion. The parallel nature of the algorithm makes it well suited to a graphics processing unit implementation, as demonstrated here. © 2017 Optical Society of America under the terms of the OSA Open Access Publishing Agreement OCIS codes: (050.1755) Computational electromagnetic methods; (060.2280) Fiber design and fabrication; (130.5296) Photonic crystal waveguides. References and links 1. P. Russell, “Photonic Crystal Fibers,” Science 299(5605), 358–362 (2003). 2. P. Russell, “Photonic-Crystal Fibers,” J. Lightwave Technol. 24(12), 4729–4749 (2006). 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). 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). 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). 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). 8. Z. Zhu and T. G. Brown, “Full-vectorial finite-difference analysis of microstructured optical fibers,” Opt. Express 10(17), 853–864 (2002). 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). 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). 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). 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). 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) 14. J. A. M. Svedin, “A Numerically Efficient FiniteElement Formulation for the General Waveguide Problem Without Spurious Modes," IEEE Trans. Microw. Theory Techn. 37(11), 1708–1715 (1989). 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) 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). 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). 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). Vol. 25, No. 23 | 13 Nov 2017 | OPTICS EXPRESS 28337 #307072 https://doi.org/10.1364/OE.25.028337 Journal © 2017 Received 13 Sep 2017; revised 24 Oct 2017; accepted 25 Oct 2017; published 31 Oct 2017 19. J. E. Goell, “A Circular-Harmonic Computer Analysis Of Rectangular Dielectric Waveguide,” Bell Syst. Tech. J. 48(7), 2133–2160 (1969). 20. Y.-H. Wang and C. Vassallo, “Circular Fourier analysis of arbitrarily shaped optical fibers," Opt. Lett. 14(24), 1377–1379 (1989). 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). 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). 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). 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). 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). 26. Z. Zhu and T. G. Brown, “Multipole analysis of hole-assisted optical fibers,” Opt. Commun. 206(4), 333–339 (2002). 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). 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). 29. P. Lalanne and E. Silberstein, “Fourier-modal methods applied to waveguide computational problems,” Opt. Lett. 25(15), 1092–1094 (2000). 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). 31. G. Sztefka and H. P. Nolting, “Bidirectional eigenmode propagation for large refractive index steps,” IEEE Photon. Technol. Lett. 5(5), 554–557 (1993). 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). 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). 35. C. Vassallo, “1993—1995 Optical mode solvers,” Opt. Quant. Electron. 29(2), 95–114 (1997). 36. J.-R. Qian and W.-P. Huang, “Coupled-mode theory for LP modes,” J. Lightwave Technol. 4(6), 619–625 (1986). 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). 38. W.-P. Huang and J. Mu, “Complex coupled-mode theory for optical waveguides,” Opt. Express 17(21), 19134–19152 (2009). 39. A. Yariv, “Coupled-mode theory for guided-wave optics,” IEEE J. Quant. Electron. 9(9), 919–933 (1973). 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). 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). 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). 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). 45. A. Mohammadi, H. Nadgaran, and M. Agio, “Contour-path effective permittivities for the two-dimensional finitedifference time-domain method,” Opt. Express 13(25), 10367–10381 (2005). 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). 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). 49. R. Gordon, “Reflection of Cylindrical Surface Waves,” Opt. Express 17(21), 18621–18629 (2009).

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.

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.

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, 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 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 ∂a j (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 The term e jt • e kt means e j x • e k x + e jy • e ky .The factors N j , 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).

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, ˜ = 0 n 2 .We consider the basis in vacuum, such that = 0 (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: The wavevector λ gives the effective index of a given mode via λ = n eff k, 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, 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

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.

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 e z = 0, and, The TM nm modes are e z = (A sin nφ + B cos nφ)J n (k c ρ), and, 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, µ = µ 0 is the vacuum magnetic permeability and 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 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.

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.

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), The transverse terms in I are I 1 + I 2 = e j x • e k x + e jy • e ky .But the transformation mapping the transverse part of the basis modes from polar coordinates to Cartesian coordinates is Ones sees that e j x • e k x + e jy • e ky = e jρ • e kρ + e jφ • e kφ .In this way, I 1 and I 2 are each separable into products of functions, R i (ρ)Φ i (φ).The longitudinal or z term in I, I 3 , 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, 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.

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

Parallelization
With N basis modes, CCMT requires the computation of N normalization factors N j and approximately 3N 2 /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.

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

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][45][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.

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 In the cladding it is: (17) 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 The values u and k 1 are in the core; w and k 2 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] J ν (ub) uJ ν (ub) This characteristic equation can be solved numerically, yielding the set of modal wavevectors β.
The quantity b is the radius of the fiber.

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]: The functions I are modified Bessel functions of the first kind.The quantities p m,d are given by p m,d = β 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.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.

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

A lossy dielectric shell
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.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.

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 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.33iand the analytic result is 5.81+0.34i.
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 β.
Table 4. Computed effective indexes for the plasmonic nanowire shown in Fig. 6.

Method
Effective Index, n eff CCMT 5.93+0.35iFDE 5.92+0.33iAnalytic 5.81+0.34i 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.

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][32][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.

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 Step-index fiber; n f = 1.6.

Fig. 1 .
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.

Fig. 2 .
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.

Figure 1 (
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.Figure4shows the CCMT eigenmodes.In this thin shell, the staircasing of the circular fiber

Fig. 3 .
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.

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

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

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

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