Guided modes of circular multi-step index optical fibers

An online solver for the guided modes supported by circular dielectric optical fibers with radially piecewise constant refractive index profiles. Following the fiber definition in terms of (outer) core radius, refractive indices, and thicknesses of intermediate layers (if applicable), and the specification of a 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.

Input

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
n_{i} (interior core region),
n_{1}, ... ,
n_{N} (intermediate layers 1 to N),
n_{e} (exterior cladding region),
and thicknesses
t_{1}, ... ,
t_{N} 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:

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 [l_{min}, l_{max}] 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 [n_{e}, n_{max}] for potential effective mode indices, where n_{e} is the refractive index of the outermost cladding region, and n_{max} 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 n_{i}, n_{1}, …, n_{N}, n_{e} 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 [l_{min}, l_{max}] of angular mode orders is provided, the solver starts with the lowest value l_{min} 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 l_{max} 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 l_{h} is suspected to have been overlooked, one might try a direct solver run for the angular mode order specification [l_{h}, l_{h}].

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.

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
λ.
B is a normalized effective permittivity, the
ratio
B = (N_{neff}^{2}-n_{e}^{2})/(
n_{max}^{2}-n_{e}^{2}), with
the maximum refractive index n_{max} of all regions.
n_{e} 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
L_{c} = π/|β_{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 S_{z} 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

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

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 E_{r}, …, H_{z} 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) = (H_{r}, 0, H_{z})(r), i.e. are characterized by a purely transverse electric field, with zero E_{z} component.
Likewise, the radial profiles of **transverse magnetic (TM) modes** can be written as
**E**(r) = (E_{r}, 0, E_{z})(r), and
**H**(r) = (0, H_{θ}, 0)(r) with a purely transverse magnetic field and zero H_{z} 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
E_{z},
E_{θ}, and
H_{r} are real, while
H_{z},
H_{θ}, and
E_{r} 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/µ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).
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 S_{z} 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 ψ ∈ {E_{r}, …, 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**,

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.