An online solver for the guided modes supported by circular dielectric optical fibers with radially piecewise constant refractive index profiles. For a fiber defined in terms of (outer) core radius, refractive indices, and thicknesses of intermediate layers (if applicable), and a specific vacuum wavelength, the script calculates the effective mode indices and propagation constants of the vectorial, hybrid modes supported by the fiber. Facilities for detailed inspection of the mode profiles, and for exporting data and figures, are provided.


For a circular fiber with N intermediate layers, the input mask receives the vacuum wavelength λ, the radius R of the outer rim of the radial layer system, refractive index values ni (interior core region), n1, ... , nN (intermediate layers 1 to N), ne (exterior cladding region), and thicknesses t1, ... , tN of the intermediate layers. Select N=0 to specify a uniform circular dielectric core. All dimensions are meant in micrometers. The figure illustrates the relevant geometry:

Circular dielectric multi-step-index optical fiber, cross section

Polar coordinates r, θ span the Cartesian x-y-plane; the center of the fiber cross section is located at the origin of both coordinate systems. The refractive index profile is independent of θ and piecewise constant in the radial direction r; further, the refractive index distribution is assumed to be constant along the z-axis (perpendicular to the x-y- and r-θ-plane, not shown).

The vectorial modal eigenproblem is independent of the azimuthal angle θ. In the formulation employed here, all modes are thus characterized by an integer angular mode order l, which corresponds to a θ-dependence ∼exp(-i l θ) of all electromagnetic field components. In certain contexts, this number l is called the topological charge of the vortex field associated with the orbital-angular-momentum (OAM) -modes. The solver accepts a range [lmin, lmax] of positive integer angular mode orders. Supply an interval [l, l] to calculate modes of one specific order l only.

Solver details

For the given vacuum wavelength, the solver examines the range [ne, nmax] for potential effective mode indices, where ne is the refractive index of the outermost cladding region, and nmax is the maximum of the refractive indices of the central core region and of the interior layers. Valid effective indices are identified as roots in a transverse resonance condition using a numerical search procedure with a heuristically defined initial stepsize. While modes of different angular orders belong to independent eigenvalue problems, and are treated as such, the search procedure might fail to identify guided modes of the same angular order with too close effective indices.

The formalism as implemented becomes singular in case a trial effective index is too close to one of the physical refractive index levels involved. Small intervals around the values ni, n1, …, nN, ne are deliberately excluded from the numerical search for roots in the transverse resonance condition. Consequently, if the fiber structure in question supports a guided mode with an effective index very close to these levels, the solver might miss that mode.

In case a range [lmin, lmax] of angular mode orders is provided, the solver starts with the lowest value lmin and proceeds with eigenproblems for larger l. Assuming that the number of solutions decays with increasing angular order, the solver stops if no modes could be identified for some order l, even if the maximum order lmax has not been examined (exception: the order l = 1, if in the specified interval, is always examined). Therefore, if a mode of high angular order lh is suspected to have been overlooked, one might try a direct solver run for the angular mode order specification [lh, lh].

Hence, while the found solutions should be valid guided modes (check: radial profile shapes, interface conditions), one must not rely too much on the solver being able to find all guided modes supported by the fiber structure, nor on the correct classification (second index) of the modes. In many cases, however, these heuristics have been observed to work adequately.


A table shows, for each guided mode, 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 λ. B is a normalized effective permittivity, the ratio B = (Nneff2-ne2)/( nmax2-ne2), with the maximum refractive index nmax of all regions. ne denotes the refractive index of the outermost cladding layer.

The mode identifiers comprise two indices. A first index indicates the angular order l of the mode (cf. the discussion below of the functional dependence of the fields on the cylindrical coordinates). The second index counts the position of the mode in question in an ordered list of modal solutions found for this particular angular order, starting with index 1 (not 0) for the mode with highest effective index (see also the remarks on the limitations of the modesolver, on potentially missed solutions, in a paragraph above). Hence, the second index should not be interpreted as a (e.g. radial) mode order.

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 pairs 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 axial component Sz of the Poynting vector. Absolute and relative expressions (confinement factors, Γ) are evaluated piecewise for each ring-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 radial position and thickness with equal refractive index.

Referring to the polar coordinate system as introduced above, this concerns optical electric fields E and magnetic fields H that depend on the spatial cross sectional coordinates r, θ (these span the x-y-cross sectional plane), on the axial propagation coordinate z, and on time t, as

E(r, θ, z, t) = Re{E(r) exp(i ω t-i l θ-i β z)},   H(r, θ, z, t) = Re{H(r) exp(i ω t-i l θ-i β z)},

with angular frequency ω = k c = 2 π c/λ, integer angular mode order l, and real propagation constant β = k neff.

For nonzero angular mode order |l| ≥ 1, these are hybrid vectorial modes; in this case, all six electromagnetic components of the radial mode shape E, H are nonzero. These orbital-angular-momentum (OAM) modes correspond to vortex fields with a topological charge l. OAM modes of angular orders l and -l are degenerate. Plots of the physical field components Er, …, Hz on the cross sectional plane show profiles that rotate anticlockwise (l ≥ 1) or clockwise (l ≤ -1) for progressing time. Superpositions of the directional variants of an OAM mode with angular order ±l constitute further valid modal eigensolutions for the same effective index eigenvalue. Relative amplitudes of unit absolute value lead to the well-known hybrid modal solutions of HE- or EH-type, where the relative phase determines the orientation of the symmetry plane (if any) and the polarization of the HE/EH modes.

Solutions for zero angular mode order l = 0 are independent of the coordinate θ. Two classes of modes with different polarization are distinguished. The radial shapes of transverse electric (TE) fiber modes are of the form E(r) = (0, Eθ, 0)(r), and H(r) = (Hr, 0, Hz)(r), i.e. are characterized by a purely transverse electric field, with zero Ez component. Likewise, the radial profiles of transverse magnetic (TM) modes can be written as E(r) = (Er, 0, Ez)(r), and H(r) = (0, Hθ, 0)(r) with a purely transverse magnetic field and zero Hz component. Note that it depends on the refractive index layering whether a TE/TM or an OAM mode is the fundamental mode of the fiber structure.

In many cases, modes of potentially all of the above types appear in groups with close-by effective indices. This can be observed in particular for fiber configurations with small refractive index contrast. Superpostions of these nearly-degenerate modes (not directly supported by the solver) then constitute the (approximate) linear polarized (LP)-modes observed in weakly-guiding fibers.

Mode profile inspection

Being solutions of eigenvalue problems, the fiber modes are determined up to some complex constant only. Somewhat arbitrarily, the solver chooses a global phase for the electromagnetic components of the mode profile, such that Ez, Eθ, and Hr are real, while Hz, Hθ, and Er are imaginary (note that this does not hold for superpositions of OAM modes). No units are shown for the electric or magnetic fields. Still, units of V/µm can be assumed for all electric fields, 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). Then the values given for the electromagnetic mode profile components correspond to a normalization of individual modes to unit power flow: The integral over the x-y-cross sectional plane of the longitudinal component Sz of the Poynting vector evaluates to 1 W (but note the increased power associated with a superposition of OAM±l,m-modes).

Plots of the radial mode shapes show the real and imaginary parts of the complex field, its absolute value, or the absolute squared profile. The background shading indicates the dielectric structure, where darker shading means higher refractive index. 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 radial 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 ψ ∈ {Er, …, Hθ} of the complex-valued radial electromagnetic mode profile relate to time-varying physical fields Ψ(r, t) = Re ψ(r) exp(i ω t). The animations show the respective field at recurring points in time.

Alternatively, plots of the mode profile on the cross sectional plane can be displayed. These show a component of the total physical electromagnetic field E, H on a rectangular domain in the Cartesian x-y-plane. The global phase of the solution can be adjusted; the animations show the respective component at recurring points in time.

Referring to the coordinates as introduced in the figure, the OAM modes as calculated by the solver, with positive angular mode order l ≥1, represent fields that propagate anticlockwise in positive θ-direction. For each mode there exists an eigenfunction with the same effective index and negative angular mode order -l, that corresponds to a clockwise propagating wave. The solver permits to inspect the respective fields (Rev-checkbox). Note that the profiles of the reversed modes differ from the profiles of the mode with positive angular orders in the signs of certain field components. Superpositions of those two degenerate modes, with arbitrary complex coefficients, constitute further valid solutions for that particular effective index. Options for inspecting the fields of those superimposed states are available for the cross sectional profile plots.

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.