An analytic solver for the 1-D eigenmodes of optical multilayer step-index slab waveguides made from media with in general complex refractive index and complex permittivity. Given a potentially guiding configuration defined in terms of refractive index and/or permittivity, layer thicknesses, if applicable, and the vacuum wavelength, the script calculates the complex effective indices of the modes supported by the structure, or their phase propagation and attenuation constants, respectively, and allows to inspect the corresponding optical field patterns. Bound and leaky modes of structures with loss or gain can be identified, as well as the bound solutions related to the propagation of surface-plasmon-polaritons (SPPs) along metal-dielectric interfaces. The solver is intended as a basic tool for integrated optics design, in particular for purposes of demonstration.


For a slab structure with N intermediate layers, the input mask receives the vacuum wavelength λ, a specification of the polarization, and complex refractive index values ns (half-infinite substrate region), n1, ... , nN (inner layers 1 to N), nc (half-infinite cover region), and thicknesses t1, ... , tN of the intermediate layers. Alternatively, values εs, ε1, ... , εN, and εc for the complex permittivity can be provided. Select N = 0 to specify a single interface only. All dimensions are meant in micrometers. The figure illustrates the relevant geometry:

complex multilayer step-index slab waveguide, geometry

Cartesian coordinates x, y, z are oriented such that the x-axis is perpendicular to the layer interfaces, while the structure is constant along y and z. Light propagates in direction z; all electromagnetic fields are constant along y. The structures are supposed to be nonmagnetic with relative permeability μ = 1 at the considered frequency; permittivity ε and refractive index n are related as ε = n2. For the frequency-domain model we assume a time dependence ∼exp(i ω t) for all fields, with positive frequency ω = kc = 2πc/λ, for vacuum speed of light c, vacuum wavenumber k, and vacuum wavelength λ. This implies that an attenuating medium is characterized by a refractive index with positive real part and negative imaginary part. Likewise, the properties of an active medium with gain are given by a refractive index with positive real part and positive imaginary part.

Solver procedures

The solver relies on a heuristic procedure for generating initial guesses for the roots of the transverse resonance condition in the complex plane. A complex secant method then converges the initial guesses numerically to actual roots. Further heuristics are applied for the removal of doublets, and for the classification and ordering of these roots.

Initial guesses are generated by considering the structure with only the real part of the permittivity profile (note that this well includes the imaginary parts of the complex refractive indices, to some degree). The real propagation constants of bound modes found for these permittivity profiles then serve as starting points for tracing roots in the complex plane, for structures with gradually increasing imaginary parts of the complex permittivity. The number of intermediate configurations can be selected via the element “Im ε, #steps”. Increasing this parameter might lead to further solutions being found, if the solver must be suspected to have missed a root. This applies in particular to configurations with near-degenerate solutions (decoupled regions), and/or to configurations with balanced nonzero loss and gain.

Among the solutions with non-real effective indices, bound modes and leaky modes differ in the boundary conditions satisfied by the mode profiles at x = ± ∞, or in the choice of elementary solutions for the pieces of the profile in the half-infinite substrate and cover regions, respectively. If the option “bound & leaky modes” is selected, the solver starts with initial guesses computed for a permittivity profile with potentially reduced permittivity εext on one or both of the exterior regions. By de-selecting the “Auto”-checkbox, a value can be directly provided; otherwise, the solver uses the minimum permittivity level of all layers (N ≤ 1), or the minimum level on the interior layers (N ≥ 2) as a default. Also here the parameter “Im ε, #steps” controls the number of intermediate configurations with gradually raised outer permittivity that are considered when tracing the complex modal eigenvalues of the leaky modes.

While the found modes should be valid solutions in all cases, one must not rely too much on the scripts being able to find the most relevant modes, or to find all existing modes, for that matter. By construction, the solver is intended for configurations where the true complex eigensolutions are located close to the real modes of a “closeby” real structure. In many cases, however, these heuristics have been observed to work adequately.


A table shows, for each mode that has been found, the complex mode eigenvalue in the form of either

related as γ = k Neff, where k = 2π/λ is the vacuum wavenumber associated with the specified vacuum wavelength λ.

In case of strictly bound fields the mode identifier indicates the polarization and a mode order, the number of nodes in the real part of the principal electric component Ey of the profile of TE modes, and in the principal magnetic component Hy of the profile of TM modes, determined for the real mode that served as the initial guess (cf. the preceding section) for the true complex field. For leaky modes, that order index is replaced by two characters, one of “l” or “b”, for the type of boundary condition / elementary field type that the mode obeys in the substrate (first index) and cover regions (second index). Here “l” indicates a leaky wave that propagates outwards, accompanied by an outwards more or less strong exponential growth of the mode profile, while “b” hints at a normalizable, outwards exponentially decaying bound profile.

The solver determines forward modes, with positive propagation constant β, or positive real part of effective mode index. In line with the convention adopted for the refractive index (see first section), modal attenuation is indicated by a positive attenuation constant α, or a negative imaginary part of the effective index, respectively. Likewise, a negative α or positive imaginary part of Neff indicates modal gain. In case of an attenuated mode, the solver also shows a propagation length Lp = 1/(2 α), the distance over which the intensity (∼ absolute squared field) associated with the mode propagation diminishes by a factor of 1/e. Likewise, in case of a mode with gain, the solver shows a growth length Lg = -1/(2 α), the propagation distance over which the modal intensity increases by a factor of e.

If the interior stack contains (thick) layers with, when compared to other guiding regions, low permittivity, it can happen that the multilayer structure supports modes that are, for all practical purposes, restricted to some part of the structure only. A typical example would be a directional coupler (here 1-D), consisting of two high index cores with a thick buffer layer of low index material in between. With increasing distance between the cores, the “supermodes” of the full structure become more and more degenerate, i.e. “numerically decoupled”. The solver detects these scenarios, and replaces the true supermodes by modes of the individual cores, with (near) zero amplitudes in the respective other core. Modes determined in this way show up with an additional entry (x) in the list of effective indices and propagation constants. Note that these fields are true eigenfunctions for all numerical and practical purposes.

For a structure that supports more than one mode, the program calculates the coupling lengths or half beat lengths Lc = π/|β0 - β1| that correspond to the interference pattern of each pair of modes with different propagation constants β0 and β1. Note that also any (pronounced) modal attenuation or gain, if present, influences the interference pattern.

Referring to the Cartesian coordinate system as introduced above, the present frequency-domain model concerns optical electric fields E and magnetic fields H that depend on the spatial coordinates x, z (these span the plane of interest), and on time t, with angular frequency ω, as

E(x, z, t) = Re{E(x) exp(i ω t-i γ z)},   H(x, z, t) = Re{H(x) exp(i ω t-i γ z)}.

All fields are constant along the y-direction. The profiles of TE slab modes are of the form E(x) = (0, Ey, 0)(x), and H(x) = (Hx, 0, Hz)(x), where E and H are the electric and magnetic parts of the mode profile, respectively, depending on the normal coordinate x only. Likewise, the profiles of TM slab modes can be written as E(x) = (Ex, 0, Ez)(x), and H(x) = (0, Hy, 0)(x). In general, the profile components Ex, …, Hz are complex functions of x, with a phase that may also vary with x.

Being solutions of eigenvalue problems, the mode profiles are determined up to some complex constant only. No units are shown for their electric or magnetic fields. Still the given values correspond to a normalization of the modes to unit power flow: The integral along the x-axis of the longitudinal component Sz of the Poynting vector evaluates to 1 W/µm (power per lateral (y) unit length). Correspondingly, all electric fields are given in units of V/µm, magnetic fields are measured in A/µm, the components of the Poynting vector S have units of V·A/µm2 = W/µm2, and the electromagnetic energy density w is measured in W·fs/µm3. In this context the vacuum permittivity and permeability, respectively, are ε0 = 8.85·10-3 A·fs/(V·µm) and µ0 = 1.25·103 V·fs/(A·µm). In case of (partly) leaky modes with outwards exponentially growing profiles on one or both of the exterior regions, the normalization is computed by excluding the respective region(s) from the power integration. The global phase of the mode profiles is chosen such that the profile becomes real at the position of the maximum of the principal field component (Ey for TE, Hy for TM polarization) on the domain of the interior layers (or at the single interface, for N = 0).

Integrals of the mode profiles can serve for purposes of normalization, for an assessment of the confinement, or for the evaluation of perturbational expressions. The script evaluates respective integrals of the squared principal components |Ey|2, |Hy|2 of TE and TM modes, integrals of the squared absolute values |E|2, |H|2 of the vectorial electric and magnetic fields, and integrals of the longitudinal component Sz of the Poynting vector. Where possible (leakage, outwards unbounded fields), absolute and relative expressions (confinement factors, Γ) are evaluated piecewise for each layer with constant refractive index. In case of a local perturbation in only part of a layer, you might wish to split the respective region by introducing additional layers of proper location and thickness with equal refractive index.

Mode profile plots show the real and imaginary parts of the complex field profile, its absolute value, or the absolute squared profile. The background shading indicates the medium properties, where darker shading means higher permittivity. After selecting "Plot", the extent of the vertical axis is being adjusted such that it covers the maximum values, determined separately for the electric field strength, magnetic field strength, Poynting vector, and the energy density, over all modes (and all their field components) that have been identified by the solver, on a default x range. This is to make the plots comparable. Select the button labeled "↕" to adjust the vertical plot range to the functions that are actually displayed.

Single components ψ ∈ {Ex, …, Hz} of the complex-valued electromagnetic mode profile relate to time-varying physical fields Ψ(x, t) = Re ψ(x) exp(i ω t) at z = 0. The animations show the respective component at recurring points in time, equally distributed over the period of 2π/ω.

Alternatively, propagation plots can be displayed. These show a component of the total physical electromagnetic field E, H on a rectangular domain in the x-z-plane, for a unit (initial) amplitude of the mode at z = 0. The global phase of the solution can be adjusted; the animations show the respective component at recurring points in time, equally distributed over the period of 2π/ω.

Select points on the planes of the plots to inspect precise local field levels. Clicks outside the actual axes close the coordinate displays. In case of a propagation plot, the "▯"-button toggles a colorbar. Select a field level on the colorbar to superimpose the field plot with a pseudo-contour at that level. Also here the contours are removed by a click in the colorbar area outside the axes.