An approximate, semi-analytic mode solver for dielectric integrated optical waveguides with two-dimensional light confinement and weak lateral guiding. Following the waveguide definition in terms of refractive indices, layer thicknesses, slice widths, and the vacuum wavelength, the script calculates the propagation constants / effective indices of guided modes and allows to inspect the corresponding optical field patterns. It is intended as a basic tool for integrated optics design, in particular for purposes of demonstration. Note the remarks on the variant of the effective-index approximation the mode solver relies on.

Input

This concerns lossless dielectric optical waveguides with a rectangular, piecewise constant refractive index profile.
The input mask receives the vacuum wavelength, a specification of the target polarization, and,
for a waveguide geometry with N_{l} interior layers and N_{s} interior slices,
the layer thicknesses
t_{1}, ...,
t_{Nl},
the slice widths
w_{1}, ...,
w_{Ns},
and a matrix of refractive index values
n_{0,0}, ...,
n_{Nl+1,Ns+1}. All dimensions are given in units of micrometers.
The figure illustrates an example for a waveguide with
N_{l} = 3 interior layers, and
N_{s} = 3 interior slices:

The entire structure is meant to be constant along the longitudinal z-direction; the light propagates along this axis. Note the orientation of the transverse coordinate axes x and y.

Output

A table shows, for each guided mode,
the propagation constant
β (in µm^{-1}) and the effective index
N_{eff} = β/k, where k = 2π/λ
is the vacuum wavenumber associated with the specified vacuum wavelength
λ.
The mode identifiers (orders) indicate the number of horizontal and vertical nodal lines in the profile of the principal electric
component E_{y} of TE modes, and in the principal magnetic component
H_{y} of TM modes. Note that this simple classification is only possible in the framework of the approximations that are applied here.
For a structure that supports more than one guided mode, the program calculates
the coupling lengths or half beat lengths
L_{c} = π/|β_{0} - β_{1}|
that correspond to the interference pattern of each pair of modes with
different propagation constants β_{0} and β_{1}.

In certain cases (frequently ...) the solver identifies 2-D modes with propagation constants below the level of the largest propagation constant of a 1-D slab mode supported by one of the exterior half-infinite slices. Any hybrid physical modes with a propagation constant below this level must be suspected to be subject to more or less pronounced lateral leakage (i.e. no truly guided physical mode exists). The mode solver indicates this with an entry (L) in the table of effective indices / propagation constants.

The components of the present 2-D profiles are built from 1-D slab modes associated with the refractive index distribution of a specific inner slice, and with some lateral 1-D effective permittivity profile. Symbols (x) or (y) signify that either the vertical or the horizontal constituting slab mode profile belongs to only one of the parts of a numerically decoupled structure (cf. the respective remarks for the OMS solver).

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) = (E_{x}, E_{y}, i E_{z})(x,y), and
**H**(x, y) = (H_{x}, H_{y}, i H_{z})(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 (note the imaginary unit).
The overall phase of the mode can be adjusted such that the component functions
E_{x}, E_{y}, E_{z}, H_{x}, H_{y}, and H_{z}
are real (note the imaginary units in the expressions for **E** and **H**). Consequently, after selection of
"E_{z}" or "H_{z}", 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 S_{z}
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/µm^{2} = W/µm^{2}, and the electromagnetic energy density
w is measured in W·fs/µm^{3}. In this context the vacuum permittivity and
permeability, respectively, are
ε_{0} = 8.85·10^{-3} A·fs/(V·µm) and
µ_{0} = 1.25·10^{3} V·fs/(A·µm).

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 S_{z} of the Poynting vector.
Absolute and relative expressions (confinement factors, Γ) are evalutated 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.
In case of a local perturbation in only part of a region, you might wish to split that region
by introducing additional layers and/or slices, such that the respective rectangle can be distinguished.

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 approximation involved, some of the mode components are actually zero.

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.

Limitations

The solver relies on a variational variant of the effective index method [1]. More specifically, the script implements what is
called VEIM5_{1,0,0,0} (TE) and VEIM5_{0,0,1,0} (TM)
in Section 9.2 of Ref. [1]; the expressions of Sections 7.1 and 7.2 of Ref. [1] apply.

Initially, the solver examines the slab waveguides associated with the refractive index profiles of the interior slices.
The slice that supports the guided slab mode with the largest propagation constant is selected as the reference slice
(a frame indicates this slice in the EIMS-input mask). The guided slab modes of this reference slice
serve as potential vertical profiles for the actual 2-D modes. For each of these, the solver computes a lateral effective permittivity
profile and determines the 1-D modes supported by this effective permittivity. Those are then used as horizontal profiles
for the separable 2-D mode components. One arrives at semi-vectorial 5-component approximations of the true
vectorial modes. These profiles are of the form
**E**(x, y) = (0 , E_{y}, i E_{z})(x,y),
**H**(x, y) = (H_{x}, H_{y}, i H_{z})(x,y) for (quasi-) TE modes, and
**E**(x, y) = (E_{x}, E_{y}, i E_{z})(x,y),
**H**(x, y) = (0, H_{y}, i H_{z})(x,y), for (quasi-) TM modes,
where the polarization is already determined by the polarization of the constituting vertical profile.

Obviously, in the framework of this approximation, the two transverse coordinates are not treated alike. In case of waveguides with substantial (local) lateral refractive index contrast, the restriction to separable mode components cannot cope with the singular behaviour of certain physical field components close to any corners in the refractive index profile. Trying to approximate the true hybrid physical fields as good as possible with the few available degrees of freedom, the solver may respond with mode profiles that exhibit unphysical discontinuities at (parts of) the horizontal or vertical dielectric interfaces.

Still, in many cases we've so far observed quite reasonable approximations of effective indices (that e.g. can serve nicely as initial guesses for further rigorous mode analysis with some more sophisticated numerical solver), even for specific waveguides with "large" refractive index contrast. The present mode profile plots give at least some qualitative first impression of the true fields. Nevertheless, you need to be aware that the input mask accepts values for a far larger range of structures than what can reasonably be handled by the present algorithm.

**[1]** O.V. Ivanova, R. Stoffer, M. Hammer

*A variational mode solver for optical waveguides based on quasi-analytical vectorial slab mode expansion*

arXiv:1307.1315 [physics.optics]
(University of Twente, technical report, 19 pages (2009)
)