Guided Modes API

modes

Guided modes of planar stacks of anisotropic layers.

Full-vector: any 3×3 ε and μ in any layer, so the modes of a rotated crystal film are hybrid (neither TE nor TM), and absorbing layers give complex effective indices.

A mode is a field that exists with nothing incident. At a plane inside the stack, let A be the scattering matrix of everything above it and B of everything below. Then an eigenvalue of A₀₁·B₁₀ equals 1 at a mode. The claddings' waves are continued analytically to their outgoing branch, so the round trip is analytic in n_eff. The roots are found by Beyn's contour method and refined by the secant method.

Classes:

• GuidedMode –

  A guided, possibly lossy mode of a planar stack, propagating along x.

Functions:

• slab_modes –

  Guided modes of a planar stack of anisotropic layers, propagating along x.

• slab_mode_fields –

  The field profile of a guided mode, normalised to unit power along x.

GuidedMode

Bases: NamedTuple

A guided, possibly lossy mode of a planar stack, propagating along x.

Attributes:

• effective_index (complex) –

  Complex n_eff = β/k₀; Im > 0 for a mode that decays as it propagates through an absorbing stack.

• loss_db_per_cm (float) –

  Power attenuation, 20·log₁₀(e)·k₀·Im(n_eff) per cm.

• multiplicity (int) –

  How many eigenvalues of the round trip reach 1 (2 for a degenerate pair).

• residual (float) –

  |1 − μ| at the root.

slab_modes

slab_modes(
    eps,
    thickness_um,
    wavelength_um,
    eps_top,
    eps_bottom,
    index_range=None,
    max_loss=None,
    mu=None,
)

Guided modes of a planar stack of anisotropic layers, propagating along x.

Full-vector: any 3×3 ε and μ in any layer, so the modes of a rotated crystal film are hybrid (neither TE nor TM), and absorbing layers give complex effective indices. Rotate the tensors to change the propagation direction in the plane.

Parameters:

• eps –

  Layer tensors [L, 3, 3] from the top down (lab frame, x the propagation direction).

• thickness_um –

  [L].

• wavelength_um –

  Vacuum wavelength (real).

• eps_top, eps_bottom –

  Cladding permittivities (isotropic, real).

• index_range (tuple[float, float] | None, default: None ) –

  Real n_eff to search; by default from just above the higher cladding index to a bound from ε and μ in the core. Indefinite tensors require an explicit range. Keep it above the claddings' indices: the branch point at each is where a mode stops being guided.

• max_loss (float | None, default: None ) –

  Largest Im n_eff searched; by default 5 % of the range.

• mu –

  Layer permeabilities; identity if omitted.

Returns: :class:GuidedMode s, highest effective index first (fundamental first).

Source code in hyperbolic_optics/modes.py

def slab_modes(
    eps,
    thickness_um,
    wavelength_um,
    eps_top,
    eps_bottom,
    index_range: tuple[float, float] | None = None,
    max_loss: float | None = None,
    mu=None,
) -> list[GuidedMode]:
    """Guided modes of a planar stack of anisotropic layers, propagating along x.

    Full-vector: any 3×3 ε and μ in any layer, so the modes of a rotated
    crystal film are hybrid (neither TE nor TM), and absorbing layers give
    complex effective indices. Rotate the tensors to change the propagation
    direction in the plane.

    Args:
        eps: Layer tensors ``[L, 3, 3]`` from the top down (lab frame, x the
            propagation direction).
        thickness_um: ``[L]``.
        wavelength_um: Vacuum wavelength (real).
        eps_top, eps_bottom: Cladding permittivities (isotropic, real).
        index_range: Real ``n_eff`` to search; by default from just above the
            higher cladding index to a bound from ε and μ in the core. Indefinite
            tensors require an explicit range. Keep it
            above the claddings' indices: the branch point at each is where a
            mode stops being guided.
        max_loss: Largest ``Im n_eff`` searched; by default 5 % of the range.
        mu: Layer permeabilities; identity if omitted.
    Returns:
        :class:`GuidedMode` s, highest effective index first (fundamental first).
    """
    bcm._check_x64()
    eps = jnp.asarray(eps, dtype=complex)
    thickness_um = jnp.asarray(thickness_um, dtype=float)
    if not len(eps):
        return []
    if index_range is None:
        cladding, core = _index_bounds(eps, eps_top, eps_bottom, mu)
        low, high = cladding * (1 + 1e-6), core
        if high <= low:
            return []
    else:
        low, high = index_range
        if low <= max(np.sqrt(eps_top), np.sqrt(eps_bottom)):
            raise ValueError("A guided index_range must lie above both cladding indices.")
    depth = 0.05 * (high - low) if max_loss is None else float(max_loss)
    if not high > low or not depth > 0:
        raise ValueError("Need an increasing index_range and a positive max_loss.")
    split = _split(eps)
    round_trip = _planar_round_trip(
        eps, thickness_um, wavelength_um, eps_top, eps_bottom, mu, split, _DOUBLINGS
    )
    branch = float(np.sqrt(max(np.real(eps_top), np.real(eps_bottom))))
    roots = _contour_roots(
        round_trip, (low, high), depth, _QUADRATURE, _PROBES, _TOLERANCE, 0, branch
    )
    k0 = 2 * np.pi / float(wavelength_um)
    return [
        GuidedMode(
            effective_index=n,
            loss_db_per_cm=float(20 * np.log10(np.e) * k0 * n.imag * 1e4),
            multiplicity=multiplicity,
            residual=residual,
        )
        for n, residual, multiplicity in sorted(roots, key=lambda r: -r[0].real)
    ]

slab_mode_fields

slab_mode_fields(
    eps,
    thickness_um,
    wavelength_um,
    eps_top,
    eps_bottom,
    effective_index,
    mu=None,
    samples_per_layer=50,
    cladding_um=(1.0, 1.0),
    cladding_samples=50,
)

The field profile of a guided mode, normalised to unit power along x.

Parameters:

• effective_index –

  A root from :func:slab_modes.

• samples_per_layer (int, default: 50 ) –

  Samples through each layer, both faces included.

• cladding_um (tuple[float, float], default: (1.0, 1.0) ) –

  How far into the top and bottom claddings to sample.

• cladding_samples (int, default: 50 ) –

  Samples in each cladding.

• Others –

  As in :func:slab_modes.

Returns:

• –

  (z_um, fields): depth below the top of the first layer [Z] and

• –

  (Ex, Ey, Ez, Hx, Hy, Hz) there [Z, 6] (H scaled by the free-space

• –

  impedance), with ½ Re ∫ (Ey·Hz* − Ez·Hy*) dz = 1: unit power along x.

The fields are marched from the internal plane with bounded scattering matrices, so thick claddings or strongly evanescent layers stay stable.

Source code in hyperbolic_optics/modes.py

def slab_mode_fields(
    eps,
    thickness_um,
    wavelength_um,
    eps_top,
    eps_bottom,
    effective_index,
    mu=None,
    samples_per_layer: int = 50,
    cladding_um: tuple[float, float] = (1.0, 1.0),
    cladding_samples: int = 50,
):
    """The field profile of a guided mode, normalised to unit power along x.

    Args:
        effective_index: A root from :func:`slab_modes`.
        samples_per_layer: Samples through each layer, both faces included.
        cladding_um: How far into the top and bottom claddings to sample.
        cladding_samples: Samples in each cladding.
        Others: As in :func:`slab_modes`.

    Returns:
        ``(z_um, fields)``: depth below the top of the first layer ``[Z]`` and
        ``(Ex, Ey, Ez, Hx, Hy, Hz)`` there ``[Z, 6]`` (H scaled by the free-space
        impedance), with ``½ Re ∫ (Ey·Hz* − Ez·Hy*) dz = 1``: unit power along x.

    The fields are marched from the internal plane with bounded scattering
    matrices, so thick claddings or strongly evanescent layers stay stable.
    """
    bcm._check_x64()
    eps = jnp.asarray(eps, dtype=complex)
    thickness_um = jnp.asarray(thickness_um, dtype=float)
    n_eff = jnp.asarray(effective_index, dtype=complex)
    w_ref, w_top, w_bottom, deltas, k0d, mu_full, k0 = _planar_parts(
        eps, thickness_um, wavelength_um, eps_top, eps_bottom, mu, n_eff
    )
    count = len(deltas)
    fractions = np.linspace(0.0, 1.0, samples_per_layer)
    step = [_layer(deltas[i], k0d[i] / (samples_per_layer - 1), _DOUBLINGS) for i in range(count)]

    # planes: every sample of every layer, top down; slices between consecutive ones
    planes = [(i, j) for i in range(count) for j in range(samples_per_layer)]
    identity = (jnp.eye(2, dtype=complex), jnp.zeros((2, 2), complex),
                jnp.zeros((2, 2), complex), jnp.eye(2, dtype=complex))  # fmt: skip

    def between(a, b):  # the slice from plane a to plane b = a + 1
        return identity if planes[b][0] != planes[a][0] else step[planes[a][0]]

    above = [_interface(w_top, w_ref)]
    for k in range(1, len(planes)):
        above.append(_star(above[-1], between(k - 1, k)))
    below = [_interface(w_ref, w_bottom)]
    for k in range(len(planes) - 2, -1, -1):
        below.insert(0, _star(between(k, k + 1), below[0]))

    split = _split(eps)
    s = split[0] * samples_per_layer + int(round(split[1] * (samples_per_layer - 1)))
    mu_values, vectors = np.linalg.eig(np.asarray(above[s][1] @ below[s][2]))
    down = {s: jnp.asarray(vectors[:, np.argmin(np.abs(1 - mu_values))])}
    up = {s: below[s][2] @ down[s]}
    eye = jnp.eye(2, dtype=complex)
    for k in range(s + 1, len(planes)):  # march down: f_k = (I − L01·B10)⁻¹ L00 f_{k−1}
        l00, l01, _, _ = between(k - 1, k)
        down[k] = jnp.linalg.solve(eye - l01 @ below[k][2], l00 @ down[k - 1])
        up[k] = below[k][2] @ down[k]
    for k in range(s - 1, -1, -1):  # march up: b_k = (I − L10·A01)⁻¹ L11 b_{k+1}
        _, _, l10, l11 = between(k, k + 1)
        up[k] = jnp.linalg.solve(eye - l10 @ above[k][1], l11 @ up[k + 1])
        down[k] = above[k][1] @ up[k]

    tops = np.concatenate([[0.0], np.cumsum(np.asarray(thickness_um))])
    z_list, field_list = [], []

    def full(tangential, eps_here, mu_here):
        a, b = _longitudinal(eps_here, mu_here, n_eff)
        ez, hz = a @ tangential, b @ tangential
        return jnp.stack([tangential[0], tangential[1], ez, tangential[2], tangential[3], hz])

    # the top cladding: outgoing (upward, decaying) waves only
    eye3 = jnp.eye(3, dtype=complex)
    top_amplitudes = jnp.linalg.solve(w_top, w_ref @ jnp.concatenate([down[0], up[0]]))
    kz_top = 1j * jnp.sqrt(n_eff**2 - eps_top + 0j)
    for z in np.linspace(-cladding_um[0], 0.0, cladding_samples, endpoint=False):
        tangential = w_top[:, 2:] @ (top_amplitudes[2:] * jnp.exp(-1j * k0 * kz_top * z))
        z_list.append(z)
        field_list.append(full(tangential, jnp.asarray(eps_top) * eye3, eye3))
    for k, (i, j) in enumerate(planes):
        tangential = w_ref @ jnp.concatenate([down[k], up[k]])
        z_list.append(tops[i] + fractions[j] * float(thickness_um[i]))
        field_list.append(full(tangential, eps[i], mu_full[i]))
    bottom_amplitudes = jnp.linalg.solve(
        w_bottom, w_ref @ jnp.concatenate([down[len(planes) - 1], up[len(planes) - 1]])
    )
    kz_bottom = 1j * jnp.sqrt(n_eff**2 - eps_bottom + 0j)
    for z in np.linspace(0.0, cladding_um[1], cladding_samples + 1)[1:]:
        tangential = w_bottom[:, :2] @ (bottom_amplitudes[:2] * jnp.exp(1j * k0 * kz_bottom * z))
        z_list.append(tops[-1] + z)
        field_list.append(full(tangential, jnp.asarray(eps_bottom) * eye3, eye3))

    z_um = np.asarray(z_list)
    fields = np.asarray(jnp.stack(field_list))
    flux = 0.5 * np.real(
        fields[:, 1] * np.conj(fields[:, 5]) - fields[:, 2] * np.conj(fields[:, 4])
    )
    order = np.argsort(z_um, kind="stable")
    power = np.trapezoid(flux[order], z_um[order])
    return z_um, fields / np.sqrt(power)