A quasianalytic mode solver for three-layer anisotropic thin film lithium niobate (TFLN) or lithium niobate on insulator (LNOI) slabs in X-cut or Z-cut configuration. The script calculates potentially direction-dependent effective indices and electromagnetic profiles of guided modes. Facilities for in-depth inspection of the vectorial mode fields, and for the convenient evaluation of dispersion curves, are provided. Most results in a comprehensive report on the topic [1] can be directly reproduced by the solver.

TFLN configuration

Lithium niobate (LN) thin films are typically prepared in either X-cut or Z-cut configuration, depending on whether the film has been cut from the birefringent bulk LN crystal in a direction perpendicular to one of its ordinary axes (X-cut) or perpendicular to the extraordinary crystal axis (Z-cut). This leads to slab waveguides with different orientations of the crystal axes {X, Y, Z} relative to the slab coordinates {x', y', z'}.

TFLN slab, X-cut configuration, coordinates
X-cut TFLN waveguide, figure adopted from [1].

TFLN slab, Z-cut configuration, coordinates
Z-cut TFLN waveguide, figure adopted from [1].

Beyond the crystal coordinates {X,Y,Z}, and the slab / device coordinates {x',y',z'}, the figures introduce a third system {x, y, z} of wave coordinates, such that the z-axis indicates the direction of guided wave propagation. All references to positions or vector components in the solver are meant with respect to that wave coordinate system. For given thin film configuration, the relative orientation of the crystal and slab coordinates is fixed. In the first place, the permittivity of the anisotropic LN layer is expressed in the slab system, then transformed to the wave system. For X-cut slabs, the outcome depends on the in-plane angle θ between the slab z' / crystal Y axes, and the direction z of wave propagation. Correspondingly, a propagation angle θ = 0° indicates a TFLN slab operated in “Y-propagation” configuration.


The solver addresses a three-layer slab waveguide, where half-infinite, isotropic substrate and cover layers with refractive indices ns and nc enclose the anisotropic core layer of thickness d, with tensor permittivity ϵf. For given choice of an X-cut or Z-cut slab, the elements of this tensor, expressed in the wave coordinate system, are determined by the ordinary and extraordinary refractive indices no, ne of the LN core medium, and by the propagation angle θ (X-cut only).

three-layer TFLN slab waveguide, geometry
Three layer slab waveguide, parameters and coordinates, adopted from [1].

Frequently, TFLN waveguides are fabricated with the LN thin film carried by a silica buffer layer on top of the actual silicon substrate. To avoid leakage, the buffer must be of a thickness such that the guided wave in the LN core is not influenced by the actual substrate. Likewise, often a protective cladding layer is applied on top of the core. If one can assume that these layers are sufficiently thick, such that the guided fields are small at their outer interfaces, the present three-layer model should accurately characterize the modal properties of the actual structure. In that situation, the buffer and cladding refractive indices should be entered in the input fields for the substrate and cover values.

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 λ. The default values for the refractive indices and vacuum wavelength have been adapted from Ref. [2].

Solver procedures

The solver operates along quasianalytical procedures as outlined in some detail in Ref. [1], relying largely on a somewhat heuristic numerical procedure for searching minima / roots of a complex expression for transverse resonance. A single computational parameter (#Steps) determines the number of cells, into which the interval for potential effective indices is divided initially. Here a higher number increases the computational load, but reduces the chance that the solver misses a mode. This concerns in particular configurations with degenerate or nearly degenerate effective indices, as occur frequently with the present TFLN slabs. Once a mode has been identified, the parameter does not influence the final accuracy of the effective index eigenvalues.

The #Steps-parameter applies also to the scan procedures. Here one could start with a low setting for a faster tentative initial overview of the parameter dependence, then select a larger #Steps for final results.


A table shows, for each mode that has been found, the modal eigenvalue in the form of the effective index Neff, and the propagation constant β, related as β = k Neff where k = 2π/λ is the vacuum wavenumber associated with the specified vacuum wavelength λ. Each mode is assigned an identifier that indicates its dominant polarization, TE or TM, and a mode order. These mode orders are distributed by ordering fields according to decreasing effective indices, then counting modes with the respective polarization, starting with zero for the fundamental modes. In this process, only the modes enter that have actually been identified by the solver (cf.\ the former remarks on the refinement of the search procedure).

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 wave 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, 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.

Modes of X-cut TFLN slabs can be pronouncedly hybrid. For a quantitative characterization of the polarization, the solver computes a ratio Π for each mode as

Π = -Re ∫ Ey* Hx dx / Re ∫ (Ex* Hy - Ey* Hx) dx.

A field with Π ≥1/2 is classified as TE, a mode with Π < 1/2 is denominated a TM mode. Note that the distinction TE / TM becomes largely obsolete in case of intermediate polarization ratios Π ≈ 1/2.

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).

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 the TFLN slab 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. Absolute and relative expressions (confinement factors, Γ) are evaluated piecewise for each layer with constant 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. Choose points on the planes of the plots to inspect precise local field levels. Clicks outside the actual axes close the coordinate displays.

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π/ω.

As a means to survey the entire vectorial profile of a mode, the button "(EH)" toggles a display that covers the Cartesian electric and magnetic components in separate small panels side by side. Field levels are comparable within each row. Controls for adjusting the global phase of the complex mode profile are provided.

Parameter scans

The Scan facilities concern variations of the quantities that define the slab waveguide, for fixed TFLN configuration. The list of options covers the vacuum wavelength λ, the refractive indices ns, nc of the substrate and cover regions, and the directional refractive indices no, ne and thickness d of the core layer. Choose the parameter interval of interest (tests for physical plausibility apply), and specify a number of samples. The solver then determines effective index levels of guided modes, for the specified range of the parameter. Note that data far varying vacuum wavelength neglects material dispersion, i.e. all refractive index values are assumed to be constant for the respective wavelength/frequency interval.

Markers colors indicate the polarization ratio Π of the mode with the respective effective index. Blue identifies a pure TE field (Π = 1), black relates to a pure TM mode (Π = 0), and intermediate states appear as levels of gray. Click on a parameter value for tabulated values of effective indices and polarization ratios.

Similar to the facilities for field inspection, Plot controls are offered that permit to enlarge (+) and to reduce (-) the size of the figure, to Export the curve data, to export the figure in SVG format, and to Detach the figure into a separate browser window. Click in the figure for a precise evaluation of curve levels. Accept a specific parameter for inspecting the fields of the modified configuration (this closes all scan-related controls).


[1]    M. Hammer, H. Farheen, J. Förstner
Guided modes of thin-film-lithium-niobate slabs
(manuscript in preparation, 2024)

[2]    L. Ebers, A. Ferreri, M. Hammer, M. Albert, C. Meier, J. Förstner, P.R. Sharapova
Flexible source of correlated photons based on LNOI rib waveguides
Journal of Physics: Photonics 4 (2), 025001 (2022)  External online source
Corrigendum: Flexible source of correlated photons based on LNOI rib waveguides (2022, J. Phys. Photonics 4, 025001)
Journal of Physics: Photonics 5 (2), 029501 (2023)  External online source