Scattering API

scattering

Numerically-stable scattering-matrix backend.

The default engine multiplies per-layer transfer matrices Mᵢ = Vᵢ·Pᵢ·Vᵢ⁻¹ with Pᵢ = diag(exp(-i·kz·k₀·d)) (waves.get_matrix). For thick / lossy / strongly-evanescent layers the propagation diagonal mixes exponentially growing and decaying terms, so the products overflow and lose precision (energy conservation drifts, eventually NaN).

This module reformulates the same 4×4 anisotropic problem as a scattering matrix cascaded with the Redheffer star product, where every propagation factor is a decaying exponential (exp(-|Im(kz)|·k₀·d) ≤ 1) — the growing matrices are never materialised — so it is stable in all cases. It is the opt-in backend behind Structure.execute(payload, backend="scattering") and produces the same reflection/transmission coefficients as the transfer method where both are well-conditioned, and correct coefficients where the transfer method fails.

Algorithm and conventions follow PyLlama (Bay, Lafait & Lequime / Bay et al., J. Opt. Soc. Am. A 39, 1431 (2022); arXiv:2012.05945), eqs 18-28, adapted to this codebase's field-vector ordering [Ex, Ey, Hx, Hy] and partial-wave sort (WaveProfile.tangential_modes: columns [t0, t1, r0, r1] = forward, forward, backward, backward). Per-medium eigenvectors come from each layer's profile (interior layers and a semi-infinite crystal exit) or from the ambient dynamical matrices A_inc / A_exit (prism, isotropic exit).

Functions:

• scattering_interface_fields –

  Tangential field [Ex, Ey, Hx, Hy] at the top of each layer, stably.

• scattering_coefficients –

  Reflection/transmission coefficients via the stable scattering-matrix method.

scattering_interface_fields

scattering_interface_fields(layers, k_0, a_p, a_s)

Tangential field [Ex, Ey, Hx, Hy] at the top of each layer, stably.

The transfer route to these fields propagates Gᵢ = Mᵢ·G_{i+1} through the per-layer matrices, which carry the growing exponentials this module exists to avoid. Here the mode amplitudes come out of the cascade instead and the field is rebuilt from each medium's own eigenvectors, so no growing term is ever formed.

For layer i, split the stack at its top interface into the cascade above it (L) and the cascade below (R). The forward amplitude entering the layer and the backward amplitude returning into it satisfy

c_f = L₀₀·a + L₀₁·c_b ,   c_b = R₁₀·c_f

-- the layer is illuminated from above by a and from below by whatever the rest of the stack sends back -- giving c_f = (I − L₀₁·R₁₀)⁻¹·L₀₀·a. The field is then P·[c_f, c_b].

Parameters:

• layers (list) –

  The executed structure's layers.

• k_0 (ndarray) –

  Canonical free-space wavenumber.

• a_p (complex) –

  Incident forward p amplitude.

• a_s (complex) –

  Incident forward s amplitude.

Returns:

• list[ndarray] –

  One [..., 4] field array per layer, ordered prism-first, matching the

• list[ndarray] –

  convention of :mod:hyperbolic_optics.fields.

Source code in hyperbolic_optics/scattering.py

def scattering_interface_fields(
    layers: list, k_0: np.ndarray, a_p: complex, a_s: complex
) -> list[np.ndarray]:
    """Tangential field ``[Ex, Ey, Hx, Hy]`` at the top of each layer, stably.

    The transfer route to these fields propagates ``Gᵢ = Mᵢ·G_{i+1}`` through the
    per-layer matrices, which carry the growing exponentials this module exists
    to avoid. Here the mode amplitudes come out of the cascade instead and the
    field is rebuilt from each medium's own eigenvectors, so no growing term is
    ever formed.

    For layer ``i``, split the stack at its top interface into the cascade above
    it (``L``) and the cascade below (``R``). The forward amplitude entering the
    layer and the backward amplitude returning into it satisfy

        c_f = L₀₀·a + L₀₁·c_b ,   c_b = R₁₀·c_f

    -- the layer is illuminated from above by ``a`` and from below by whatever
    the rest of the stack sends back -- giving
    ``c_f = (I − L₀₁·R₁₀)⁻¹·L₀₀·a``. The field is then ``P·[c_f, c_b]``.

    Args:
        layers: The executed structure's layers.
        k_0: Canonical free-space wavenumber.
        a_p: Incident forward p amplitude.
        a_s: Incident forward s amplitude.

    Returns:
        One ``[..., 4]`` field array per layer, ordered prism-first, matching the
        convention of :mod:`hyperbolic_optics.fields`.
    """
    media = [_medium(layer) for layer in layers]
    count = len(media)
    interfaces = [
        _interface_scattering(media[i][0], media[i][1], media[i][2], media[i + 1][0], k_0)
        for i in range(count - 1)
    ]

    # below[i] cascades interfaces i .. end; above[i] cascades 0 .. i-1.
    below: list[np.ndarray | None] = [None] * count
    accumulated = None
    for i in range(count - 2, -1, -1):
        accumulated = interfaces[i] if accumulated is None else _star(interfaces[i], accumulated)
        below[i] = accumulated

    above: list[np.ndarray | None] = [None] * count
    accumulated = None
    for i in range(1, count):
        accumulated = (
            interfaces[i - 1] if accumulated is None else _star(accumulated, interfaces[i - 1])
        )
        above[i] = accumulated

    # Ambient columns are reordered to [p_fwd, s_fwd, ...], so the incident
    # forward pair is (p, s) in that order.
    reference = interfaces[0][..., 0, 0]
    incident = np.stack(
        [np.full_like(reference, a_p), np.full_like(reference, a_s)], axis=-1
    )  # [..., 2]

    identity = np.eye(2, dtype=np.complex128)
    fields = []
    for i in range(count):
        if above[i] is None:  # the prism: nothing above it to correct for
            forward = incident
        else:
            l00, l01 = above[i][..., :2, :2], above[i][..., :2, 2:]
            forward = l00 @ incident[..., np.newaxis]
            if below[i] is not None:
                r10 = below[i][..., 2:, :2]
                forward = np.linalg.inv(identity - l01 @ r10) @ forward
            forward = forward[..., 0]

        if below[i] is None:  # the exit half-space carries no backward wave
            backward = np.zeros_like(forward)
        else:
            backward = (below[i][..., 2:, :2] @ forward[..., np.newaxis])[..., 0]

        amplitudes = np.concatenate([forward, backward], axis=-1)  # [..., 4]
        fields.append((media[i][0] @ amplitudes[..., np.newaxis])[..., 0])

    return fields

scattering_coefficients

scattering_coefficients(layers, k_0)

Reflection/transmission coefficients via the stable scattering-matrix method.

Parameters:

• layers (list) –

  The executed structure's layers (structure.layers), each exposing eigenvectors via profile or an ambient dynamical matrix.

• k_0 (ndarray) –

  Canonical free-space wavenumber [1, 1, F, 1].

Returns:

• dict[str, ndarray] –

  Dict of the eight complex coefficients ``r_pp, r_ss, r_ps, r_sp, t_pp,

• dict[str, ndarray] –

  t_ss, t_ps, t_spin canonical[A, B, F, T]`` layout (un-presented).

Note

Reflection coefficients are in the prism's s/p basis. For an isotropic exit the transmission coefficients are too; for a semi-infinite anisotropic exit they name that crystal's eigenmodes, ordered p-like first to match the prism -- the same convention the transfer backend uses, so the two agree coefficient for coefficient.

Source code in hyperbolic_optics/scattering.py

def scattering_coefficients(layers: list, k_0: np.ndarray) -> dict[str, np.ndarray]:
    """Reflection/transmission coefficients via the stable scattering-matrix method.

    Args:
        layers: The executed structure's layers (``structure.layers``), each
            exposing eigenvectors via ``profile`` or an ambient dynamical matrix.
        k_0: Canonical free-space wavenumber ``[1, 1, F, 1]``.

    Returns:
        Dict of the eight complex coefficients ``r_pp, r_ss, r_ps, r_sp, t_pp,
        t_ss, t_ps, t_sp`` in canonical ``[A, B, F, T]`` layout (un-presented).

    Note:
        Reflection coefficients are in the prism's s/p basis. For an isotropic
        exit the transmission coefficients are too; for a semi-infinite
        anisotropic exit they name that crystal's eigenmodes, ordered p-like
        first to match the prism -- the same convention the transfer backend
        uses, so the two agree coefficient for coefficient.
    """
    media = [_medium(layer) for layer in layers]
    scattering = None
    for (p_left, q_left, d_left), (p_right, _, _) in zip(media[:-1], media[1:], strict=True):
        interface = _interface_scattering(p_left, q_left, d_left, p_right, k_0)
        scattering = interface if scattering is None else _star(scattering, interface)

    assert_canonical(scattering, matrix_ndim=2, name="scattering_matrix")
    s = scattering
    # PyLlama r_kj is (out k, in j); this codebase uses r_{in->out}, so the cross
    # terms transpose (verified against the transfer backend).
    return {
        "r_pp": s[..., 2, 0],
        "r_ss": s[..., 3, 1],
        "r_ps": s[..., 3, 0],
        "r_sp": s[..., 2, 1],
        "t_pp": s[..., 0, 0],
        "t_ss": s[..., 1, 1],
        "t_ps": s[..., 1, 0],
        "t_sp": s[..., 0, 1],
    }