Vectorial or semi-vectorial mode analysis of dielectric integrated optical waveguides with two-dimensional light confinement. Multilayer step-index waveguides with rectangular, piecewise constant refractive index profiles are addressed. For a waveguide configuration specified in terms of refractive indices, layer thicknesses, slice widths, and the vacuum wavelength, the online script calculates the propagation constants / effective indices of guided modes and allows to inspect the corresponding optical field patterns. The quasi-analytical solver relies on expansions of the modal fields into locally exact solutions of the mode equations, applied piecewise on the rectangular regions that constitute the waveguide cross section (wave matching method, WMM), as described in Refs. [1, 2, 3, 4]. Further implementational details can also be found in the documentation that accompanies the sources of the WMM solver.


This concerns lossless dielectric optical waveguides with a rectangular, piecewise constant refractive index profile. For a waveguide geometry with Nl interior layers and Ns interior slices, the input mask receives the layer thicknesses t1, ..., tNl, the slice widths w1, ..., wNs, and a matrix of refractive index values n0,0, ..., nNl+1,Ns+1. All dimensions are given in units of micrometers. The figure illustrates an example for a waveguide with Nl = 3 interior layers, and Ns = 1 interior slice:

Rectangular waveguide geometry

The entire structure is meant to be constant along the longitudinal z-direction; the light propagates along that axis. Note the orientation of the transverse coordinate axes x and y. In case a waveguide structure is mirror symmetric with respect to a central x-z-plane, this is detected automatically and taken into account for the calculations.

Further, a value for the vacuum wavelength needs to be provided, and a choice of the target polarization, or of the type of analysis, is required. Here the alternatives QTE, QTM, and VEC indicate the semi-vectorial computation of quasi-TE or quasi-TM modes, or the fully vectorial modal analysis, respectively:


Given the parameters in the primary input mask, the problem is pre-analyzed with the help of the OMS solver. Limits for potential effective indices of any guided modes supported by the 2-D waveguide are established in the following way.

  1. The largest refractive index on one of the four doubly-half-infinite corner rectangles of the cross section decomposition establishes a rough lower limit. A rough upper limit is given by the largest refractive index of one of the interior bounded rectangular regions.
  2. The layers stack formed by the regions of the half-infinite left- and rightmost slices that constitute the 2-D cross section are analyzed for guided slab modes, of either polarization. The highest effective index of some slab mode that is guided in these stacks forms an increased lower limit to the effective indices of guided modes of the 2-D waveguide.
  3. The layer stacks formed by the interior slices of the 2-D waveguide are subjected to a 1-D slab mode analysis, also for both TE and TM polarization. The highest effective index of any guided slab mode found for those structures, or alternatively the highest refractive index of the half-infinite upper- or lowermost regions (limits to the effective indices of the local continuum of radiation modes), establishes an upper limit for the effective indices of guided modes of the 2-D waveguide.
  4. Steps 2 and 3 are repeated for the layer stacks formed by the exterior and interior layers of the 2-D waveguide. This procedure results in the interval for effective mode indices Neff as shown by the solver.

Next, the symmetry properties of the structure are considered. In case there is a vertical x-z-mirror plane at a central y position, the solver offers analysis options related to modes that are symmetric (SYM) or antisymmetric (ASY) with respect to that plane. Here the interpretation depends on the polarization / the type of analysis.

Otherwise, if no mirror symmetry y↔-y is detected, the script continues with mode analysis for fields of no particular symmetry (NOS). Then the number of unknowns is twice as large as for a simulation with prescribed symmetry, requiring a much higher numerical effort.

In a final step, the 2-D waveguide structure is pre-analyzed by the EIMS solver. The script then offers a selection of targets for the further WMM analysis procedure. For the polarization / analysis type pol as selected initially, these options differ with respect to the symmetry alternatives, if applicable, and with respect to the interval of effective indices that is being searched for potential guided modes.

The WMM progress window (see the remarks in the next section) shows the approximate modal effective indices as obtained from the EIMS analysis, if any, as thin white dash-dotted lines.

Computational parameters

While details on the algorithm can be found elsewhere, a few remarks are necessary to make reasonable choices on the values of the parameters that affect the mode analysis process.

Separately for each rectangular region with constant refractive index, the mode profiles are expanded into a number of factorizing harmonics or exponentials that solve the basic wave equation. This spectral discretization depends on the trial value for the propagation constant. Each trial function included in the set is characterized by a transverse wave vector. The approximation of the mode profile improves with the width of the interval for these wave vectors and with the density of the selected points in the wavenumber space, i.e. with the number of trial functions. At the same time, the numerical effort increases significantly with the spectral discretization density.

To isolate valid propagation constants, the mode solver searches the interval of prospective effective mode indices for the minima of an error function, the least squares mismatch of the electromagnetic fields on the dielectric boundaries. This error function is shown in the progress window that opens while the actual computations are carried out. The error functions appear as series of sections of more or less narrow parabolas, each of which can be associated with a particular guided mode.

The minimum searching procedure begins with an initial survey of this function. Marching with constant stepsize backwards from the upper end of the specified interval of potential mode wavenumbers, the error function is evaluated for a usually coarse preliminary spectral discretization. Once a minimum in the error function is trapped (i.e. a smaller ordinate appears between two larger ones), the mode solver refines the spectral discretization in several steps and repeatedly fixes the location of the minimum. The process stops when a prescribed final spectral discretization density is reached; this discretization serves to compute the final approximation for the propagation constant. The initial survey then continues until the entire interval is investigated.

The following list of parameters contains quantities both for controlling the spectral discretization and for influencing the behavior of the search algorithm:

Further, in case of fully vectorial calculations, the online solver implements the variant of the WMM algorithm as specified by the parameters vform=HXHY, vnorm=NRMMH, ccomp=CCALL; see the respective remarks in the WMM documentation.

If the solver seems to have missed a guided mode, one should consider the shape of the least squares error function, as displayed in the solver progress window. If no minimum appears at all in the region of the prospective propagation constant, one could try to refine the initial spectral discretization. If a minimum is visible which was not trapped by the mode solver, an increase in the number of initial survey steps Ns and / or a slight alteration of the stepsize for the initial survey Δα may help.

Internally, the solver uses a simple inverse iteration procedure to determine smallest eigenvalues of the matrices that relate to the least-squares error of the field at the dielectric interfaces. That numerical procedure relies on a random initial vector. Hence, repeated runs of the solver for identical problem settings can well lead to different results, what concerns the last digits of the effective index values shown. Obviously, no physical significance should be attached to these.


A table shows, for each guided mode that has been identified by the solver, the propagation constant β (in µm-1) and the effective index Neff = β/k, where k = 2π/λ is the vacuum wavenumber associated with the specified vacuum wavelength λ. The mode identifiers indicate the type of analysis (QTE, QTM, VEC) that led to the solution in question, further an indicator (s, a) for the symmetry of the mode (cf. the respective previous remarks), and an index that counts modes with this polarization/analysis type and symmetry in the set of modes as shown, ordered according to effective index (highest first).

Arrow symbols hint at the polarization character of the modes, i.e. the predominant orientation of the electric field, either along the horizontal y-axis (↔), with a major electric profile component Ey, or along the vertical x-axis ( ↕ ), with a major electric profile component Ex. The polarization characters are determined by comparing ratios of integrals of the squared transverse magnetic components of the mode profiles. By construction, QTE modes are predominantly horizontally polarized, while QTM modes are predominantly vertically polarized.

For a structure that supports more than one guided 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.

Integrals of the mode profiles can serve for purposes of normalization, for an assessment of the confinement, or for the evaluation of perturbational expressions, where standard ratios concern the influence of material loss or gain, or phase shifts caused by small changes in refractive index. The script evaluates respective 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. Absolute and relative expressions (confinement factors, Γ) are evaluated piecewise for each rectangular part of the cross section with constant refractive index, and for the layers and slices the cross section is composed of.

Field inspection

Mode profile plots show the field (real or imaginary, as discussed above), its absolute value, or the squared field. Thin white lines indicate the dielectric interfaces. After selecting "Plot", the extent of the "color 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. This is to make plots comparable. Select the button labeled "↕" to adjust the range of the colormap to the functions that are actually displayed. Note that, due to the approximations involved, some of the profile components of modes of QTE or QTM type are zero by construction.

Click in the figure for a precise evaluation of field levels. A click outside the actual axis closes the coordinate display. 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 contour is removed by a click in the colorbar area outside the axes.

Referring to the coordinate system as introduced above, the profiles of the guided modes supported by the present lossless dielectric waveguides are of the form E(x, y) = (Ex, Ey, i Ez)(x,y), and H(x, y) = (Hx, Hy, i Hz)(x,y), where E and H are the electric and magnetic parts of the 2-D mode profile, respectively, depending on the transverse coordinates x and y. The overall phase of the mode can be adjusted such that the component functions Ex, Ey, Ez, Hx, Hy, and Hz are real (note the imaginary units in the expressions for E and H). Consequently, after selection of "Ez" or "Hz", the plot window shows the imaginary part of the longitudinal electric or magnetic component of the mode profile, while for all other components the real part is displayed.

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 over the x-y-plane of the longitudinal component Sz of the Poynting vector evaluates to a power of 1 W. 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).


[1]  M. Lohmeyer
Wave-matching-method for mode analysis of dielectric waveguides
Optical and Quantum Electronics 29 (9), 907-922 (1997)  Preprint (ps.gz)  Preprint (pdf)  External online source

[2]  M. Lohmeyer
Vectorial wave-matching-method mode analysis of integrated optical waveguides
Optical and Quantum Electronics 30 (5/6), 385-396 (1998)  Preprint (ps.gz)  Preprint (pdf)  External online source

[3]  M. Lohmeyer, N. Bahlmann, O. Zhuromskyy, P. Hertel
Wave-matching-simulations of integrated optical coupler structures
Proceedings of SPIE, Vol. 3620, Integrated Optics Devices III, 68-78 (1999)  Preprint (ps.gz)  Preprint (pdf)  External online source

[4]  M. Lohmeyer
Guided waves in rectangular integrated magnetooptic devices
Cuvillier Verlag, Göttingen, 1999  ps.gz document  pdf document