Grating RCWA API

rcwa

Differentiable RCWA for 1D-periodic gratings of full-tensor anisotropic media.

Periodic in x, uniform in y, layered in z; each segment has its own 3×3 ε and μ in any orientation, and incidence may be conical. Layers use the cascade of :mod:hyperbolic_optics.bcm on 4N × 4N Fourier-space operators, so everything, the geometry's Fourier coefficients included, is differentiable.

Fourier factorisation

Across an x-discontinuity the components Dx, Ey, Ez (and Bx, Hy, Hz) are continuous. The constitutive relation is therefore first solved for (Ex, Dy, Dz) in terms of (Dx, Ey, Ez): every product in that form pairs a discontinuous coefficient with a continuous field, for which Laurent's rule (a Toeplitz matrix of Fourier coefficients) converges. The result is then rearranged back into D = [[ε]]·E with the inverse taken only in Fourier space. This is Li's factorisation for anisotropic lamellar gratings (L. Li, J. Mod. Opt. 45, 1313 (1998)); for an isotropic medium it reduces to the familiar inverse rule for Dx and Laurent's rule for Dy and Dz. factorization="laurent" uses plain Toeplitz matrices of ε (slow for TM light on high-contrast gratings).

Conventions

Fields and time dependence follow :mod:hyperbolic_optics.bcm. Harmonic m has in-plane wavevector (kx₀ + m·λ/Λ, ky) (normalised by k₀). Every order's amplitudes are resolved into that order's own p and s: s is perpendicular to its plane of diffraction, along ẑ × û with û the unit in-plane wavevector, and p completes the triad with the sign convention of the planar solver. (At exactly zero in-plane wavevector û is taken along the plane of incidence.) Results are indexed [..., order, out, in] with 0 = p and 1 = s; R and T are diffraction efficiencies, zero for evanescent orders.

Classes:

• GratingLayer –

  One layer: S segments of constant material across one period.

• GratingResult –

  Complex amplitudes and efficiencies, each [..., order, out, in].

Functions:

• uniform_layer –

  A homogeneous layer, eps [..., 3, 3].

• lamellar_layer –

  Ridge of width fill (fraction of the period) starting at offset.

• sampled_layer –

  A profile sampled at S equal-width segments, [..., S, 3, 3].

• fourier_coefficients –

  Exact Fourier coefficients c_m, |m| ≤ max_index, of a step function.

• default_doublings –

  Doublings that keep the slice error at the planar level as orders are added.

• solve_grating –

  Diffraction of a plane wave by a stack of 1D-periodic anisotropic layers.

GratingLayer

Bases: NamedTuple

One layer: S segments of constant material across one period.

Attributes:

• thickness_um (Array) –

  Layer thickness in µm, [...].

• boundaries (Array) –

  Segment edges as fractions of the period, [..., S + 1], increasing, with boundaries[S] − boundaries[0] = 1.

• eps (Array) –

  Segment permittivity tensors in the lab frame, [..., S, 3, 3].

• mu (Array | None) –

  Segment permeability tensors, same layout; None for μ = 1.

• edge_width (Array | float) –

  Optional smoothing, as a fraction of the period: the profile is convolved with a box this wide, which multiplies Fourier coefficient m by sinc(m · edge_width). 0 keeps hard edges. A few percent removes the truncation ripple that makes geometry gradients of hyperbolic (high-k) ridges unusable; see docs/user-guide/gratings.md.

GratingResult

Bases: NamedTuple

Complex amplitudes and efficiencies, each [..., order, out, in].

uniform_layer

uniform_layer(eps, thickness_um, mu=None)

A homogeneous layer, eps [..., 3, 3].

Source code in hyperbolic_optics/rcwa.py

def uniform_layer(eps, thickness_um, mu=None) -> GratingLayer:
    """A homogeneous layer, ``eps`` ``[..., 3, 3]``."""
    eps = jnp.asarray(eps)
    return GratingLayer(
        jnp.asarray(thickness_um),
        jnp.array([0.0, 1.0]),
        eps[..., None, :, :],
        None if mu is None else jnp.asarray(mu)[..., None, :, :],
    )

lamellar_layer

lamellar_layer(
    eps_ridge,
    eps_groove,
    fill,
    thickness_um,
    offset=0.0,
    mu_ridge=None,
    mu_groove=None,
    edge_width=0.0,
)

Ridge of width fill (fraction of the period) starting at offset.

fill and offset enter the Fourier coefficients analytically, so the layer is differentiable in both. edge_width smooths the edges (see :class:GratingLayer).

Source code in hyperbolic_optics/rcwa.py

def lamellar_layer(
    eps_ridge,
    eps_groove,
    fill,
    thickness_um,
    offset=0.0,
    mu_ridge=None,
    mu_groove=None,
    edge_width=0.0,
) -> GratingLayer:
    """Ridge of width ``fill`` (fraction of the period) starting at ``offset``.

    ``fill`` and ``offset`` enter the Fourier coefficients analytically, so the
    layer is differentiable in both. ``edge_width`` smooths the edges (see
    :class:`GratingLayer`).
    """
    fill, offset = jnp.asarray(fill), jnp.asarray(offset)
    boundaries = jnp.stack([offset, offset + fill, offset + 1.0], axis=-1)
    eps_ridge, eps_groove = jnp.broadcast_arrays(jnp.asarray(eps_ridge), jnp.asarray(eps_groove))
    eps = jnp.stack([eps_ridge, eps_groove], axis=-3)
    mu = None
    if mu_ridge is not None or mu_groove is not None:
        eye = jnp.eye(3)
        mu_ridge = eye if mu_ridge is None else jnp.asarray(mu_ridge)
        mu_groove = eye if mu_groove is None else jnp.asarray(mu_groove)
        mu = jnp.stack(jnp.broadcast_arrays(mu_ridge, mu_groove), axis=-3)
    return GratingLayer(jnp.asarray(thickness_um), boundaries, eps, mu, edge_width)

sampled_layer

sampled_layer(eps_samples, thickness_um, mu_samples=None)

A profile sampled at S equal-width segments, [..., S, 3, 3].

Source code in hyperbolic_optics/rcwa.py

def sampled_layer(eps_samples, thickness_um, mu_samples=None) -> GratingLayer:
    """A profile sampled at ``S`` equal-width segments, ``[..., S, 3, 3]``."""
    eps_samples = jnp.asarray(eps_samples)
    boundaries = jnp.linspace(0.0, 1.0, eps_samples.shape[-3] + 1)
    mu = None if mu_samples is None else jnp.asarray(mu_samples)
    return GratingLayer(jnp.asarray(thickness_um), boundaries, eps_samples, mu)

fourier_coefficients

fourier_coefficients(values, boundaries, max_index)

Exact Fourier coefficients c_m, |m| ≤ max_index, of a step function.

f(x) = Σ_m c_m exp(2πi·m·x) over one period, with f = values[s] on [boundaries[s], boundaries[s + 1]). Analytic in the boundaries.

Parameters:

• values (Array) –

  [..., S].

• boundaries (Array) –

  [..., S + 1].

• max_index (int) –

  Highest |m|.

Returns:

• Array –

  [..., 2·max_index + 1], index m + max_index.

Source code in hyperbolic_optics/rcwa.py

def fourier_coefficients(values: jax.Array, boundaries: jax.Array, max_index: int) -> jax.Array:
    """Exact Fourier coefficients ``c_m``, ``|m| ≤ max_index``, of a step function.

    ``f(x) = Σ_m c_m exp(2πi·m·x)`` over one period, with ``f = values[s]`` on
    ``[boundaries[s], boundaries[s + 1])``. Analytic in the boundaries.

    Args:
        values: ``[..., S]``.
        boundaries: ``[..., S + 1]``.
        max_index: Highest ``|m|``.

    Returns:
        ``[..., 2·max_index + 1]``, index ``m + max_index``.
    """
    m = jnp.arange(-max_index, max_index + 1)
    phase = jnp.exp(-2j * jnp.pi * m * boundaries[..., :, None])  # [..., S + 1, M]
    zero = m == 0
    denominator = jnp.where(zero, 1.0, -2j * jnp.pi * m)  # safe at m = 0
    integrals = jnp.where(
        zero,
        (boundaries[..., 1:] - boundaries[..., :-1])[..., None],
        (phase[..., 1:, :] - phase[..., :-1, :]) / denominator,
    )  # [..., S, M]
    return jnp.sum(values[..., None] * integrals, axis=-2)

default_doublings

default_doublings(max_order)

Doublings that keep the slice error at the planar level as orders are added.

The operator's norm grows as kx_max² ∝ max_order² and the fourth-order start leaves an error ∝ (norm · k₀d)⁵ / 2⁴ᴺ, so holding it fixed takes 2.5 · log₂(max_order / 10) extra doublings beyond ten orders. At 161 orders this is 28 doublings, which brings a lossless grating's energy error from ~1e-7 (at 20) down to its ~1e-9 rounding floor.

Source code in hyperbolic_optics/rcwa.py

def default_doublings(max_order: int) -> int:
    """Doublings that keep the slice error at the planar level as orders are added.

    The operator's norm grows as ``kx_max² ∝ max_order²`` and the fourth-order
    start leaves an error ``∝ (norm · k₀d)⁵ / 2⁴ᴺ``, so holding it fixed takes
    ``2.5 · log₂(max_order / 10)`` extra doublings beyond ten orders. At 161
    orders this is 28 doublings, which brings a lossless grating's energy error
    from ~1e-7 (at 20) down to its ~1e-9 rounding floor.
    """
    if max_order <= 10:
        return DEFAULT_DOUBLINGS
    return DEFAULT_DOUBLINGS + math.ceil(2.5 * math.log2(max_order / 10))

solve_grating

solve_grating(
    layers,
    period_um,
    wavelength_um,
    incidence_rad,
    eps_incident,
    eps_exit,
    max_order,
    azimuth_rad=0.0,
    factorization="li",
    doublings=None,
)

Diffraction of a plane wave by a stack of 1D-periodic anisotropic layers.

Parameters:

• layers (list[GratingLayer]) –

  :class:GratingLayer s from the incident side down; build them with :func:uniform_layer, :func:lamellar_layer or :func:sampled_layer.

• period_um (Array) –

  Grating period in µm, along x.

• wavelength_um (Array) –

  Vacuum wavelength in µm.

• incidence_rad (Array) –

  Polar angle of incidence in the incident medium.

• eps_incident (Array) –

  Real, positive permittivity of the incident half-space.

• eps_exit (Array) –

  Real, positive permittivity of the exit half-space.

• max_order (int) –

  Keep diffraction orders −max_order … max_order (static).

• azimuth_rad (Array, default: 0.0 ) –

  Angle of the plane of incidence from the grating vector (x); nonzero gives conical diffraction.

• factorization (str, default: 'li' ) –

  "li" (default) or "laurent" (static).

• doublings (int | None, default: None ) –

  Doublings per layer (static), as in :mod:hyperbolic_optics.bcm. Defaults to :func:default_doublings for max_order.

Returns:

• GratingResult –

  class:GratingResult; orders lists the orders, and r, t,

• GratingResult –

  R, T are [..., order, out, in] with 0 = p and 1 = s. Only the

• GratingResult –

  zeroth order is incident.

Source code in hyperbolic_optics/rcwa.py

def solve_grating(
    layers: list[GratingLayer],
    period_um: jax.Array,
    wavelength_um: jax.Array,
    incidence_rad: jax.Array,
    eps_incident: jax.Array,
    eps_exit: jax.Array,
    max_order: int,
    azimuth_rad: jax.Array = 0.0,
    factorization: str = "li",
    doublings: int | None = None,
) -> GratingResult:
    """Diffraction of a plane wave by a stack of 1D-periodic anisotropic layers.

    Args:
        layers: :class:`GratingLayer` s from the incident side down; build them
            with :func:`uniform_layer`, :func:`lamellar_layer` or
            :func:`sampled_layer`.
        period_um: Grating period in µm, along x.
        wavelength_um: Vacuum wavelength in µm.
        incidence_rad: Polar angle of incidence in the incident medium.
        eps_incident: Real, positive permittivity of the incident half-space.
        eps_exit: Real, positive permittivity of the exit half-space.
        max_order: Keep diffraction orders ``−max_order … max_order`` (static).
        azimuth_rad: Angle of the plane of incidence from the grating vector
            (x); nonzero gives conical diffraction.
        factorization: ``"li"`` (default) or ``"laurent"`` (static).
        doublings: Doublings per layer (static), as in :mod:`hyperbolic_optics.bcm`.
            Defaults to :func:`default_doublings` for ``max_order``.

    Returns:
        :class:`GratingResult`; ``orders`` lists the orders, and ``r``, ``t``,
        ``R``, ``T`` are ``[..., order, out, in]`` with 0 = p and 1 = s. Only the
        zeroth order is incident.
    """
    _check_x64()
    if doublings is None:
        doublings = default_doublings(max_order)
    size = 2 * max_order + 1
    orders = jnp.arange(-max_order, max_order + 1)
    eps_incident, eps_exit = jnp.asarray(eps_incident), jnp.asarray(eps_exit)
    wavelength_um, azimuth_rad = jnp.asarray(wavelength_um), jnp.asarray(azimuth_rad)

    k_par = jnp.sqrt(eps_incident) * jnp.sin(jnp.asarray(incidence_rad))
    kx = (k_par * jnp.cos(azimuth_rad))[..., None] + orders * (
        wavelength_um / jnp.asarray(period_um)
    )[..., None]
    ky = k_par * jnp.sin(azimuth_rad)
    direction = (jnp.cos(azimuth_rad)[..., None], jnp.sin(azimuth_rad)[..., None])
    k0 = 2 * jnp.pi / wavelength_um

    ky = jnp.asarray(ky)[..., None]  # shared by every harmonic of a 1D grating
    w_ref = _reference_basis(kx, ky, direction)
    w_in, flux_in = _harmonic_modes(eps_incident[..., None], kx, ky, direction)
    w_out, flux_out = _harmonic_modes(eps_exit[..., None], kx, ky, direction)

    scattering = _interface(w_in, w_ref)
    for layer in layers:

        def fourier(tensor, layer=layer):
            tensor = jnp.asarray(tensor) + 0j
            return _factorize(tensor, layer.boundaries, max_order, factorization, layer.edge_width)

        eps = fourier(layer.eps)
        if layer.mu is None:  # μ = 1
            mu = jnp.einsum(
                "ij,mn->ijmn", jnp.eye(3, dtype=eps.dtype), jnp.eye(size, dtype=eps.dtype)
            )
        else:
            mu = fourier(layer.mu)
        delta = _fourier_system_matrix(eps, mu, kx, ky)
        delta_ref = _solve(w_ref, delta @ w_ref)
        k0d = k0 * jnp.asarray(layer.thickness_um)
        scattering = _star(scattering, _layer(delta_ref, k0d, doublings))
    scattering = _star(scattering, _interface(w_ref, w_out))

    r, t, R, T = _order_results(scattering, flux_in, flux_out, max_order, size)
    return GratingResult(orders=orders, r=r, t=t, R=R, T=T)