Crossed-Grating RCWA API¶

rcwa2d ¶

Differentiable RCWA for 2D-periodic (crossed) gratings of full-tensor media.

Any Bravais lattice (:class:Lattice); each patterned layer is a background and non-overlapping polygons, each with its own 3×3 ε and μ. Polygons enter through their exact Fourier transforms (a closed sum over edges), so every coefficient is differentiable in every vertex. Harmonics G = m·b₁ + n·b₂ are kept for |m| ≤ M₁, |n| ≤ M₂ (truncation="parallelogram") or |G| ≤ max(M₁, M₂)·g_min ("circular", independent of the primitive vectors chosen).

Fourier factorisation¶

factorization="laurent" uses Toeplitz matrices of ε directly: correct, but slow to converge where the normal component of E jumps at an edge.

factorization="normal" (default) is the normal-vector method generalised to full tensors. At every point the constitutive relation is written in a local frame (n, t, z) whose n is normal to the pattern's edges, and Li's anisotropic rearrangement (see :mod:hyperbolic_optics.rcwa) is applied there: (E_n, D_t, D_z) in terms of the edge-continuous (D_n, E_t, E_z), with the inverse rule on the normal component only. For isotropic media this is the textbook normal-vector method; for stripes, n = x̂ and it is the 1D solver's Li rule exactly. The normal field is the gradient direction of the pattern smoothed over a tenth of the cell.

Coefficient functions that depend on the frame are products of a shape's indicator -- whose coefficients stay exact -- with a smooth function of the frame, whose coefficients come from a fine grid; the product is a convolution in Fourier space. The assembly is averaged with its product-reversed twin, which keeps it Hermitian for lossless media: energy is conserved to rounding at every truncation.

Conventions follow :mod:hyperbolic_optics.rcwa: order (m, n) has in-plane wavevector k∥₀ + (m·b₁ + n·b₂)/k₀ (on a rectangular lattice, (kx₀ + m·λ/Λx, ky₀ + n·λ/Λy)), resolved into that order's own p and s; results are [..., order, out, in] with 0 = p and 1 = s.

Classes:

Lattice –

A 2D Bravais lattice: primitive vectors a1, a2 as (x, y) in µm.
Shape –

A polygon of one material within the unit cell.
PatternedLayer –

One layer: a background material and non-overlapping :class:Shape s.
CrossedGratingResult –

Amplitudes and efficiencies [..., order, out, in].

Functions:

uniform –

An unpatterned layer.
rectangle –

Vertices [4, 2] of an axis-aligned rectangle (differentiable in both).
ellipse –

Vertices [n, 2] of an ellipse, as a polygon of equal area.
polygon_fourier –

∫_P exp(−i(gx·x + gy·y)) dA over a polygon, exactly.
crossed_grating_smatrix –

The zeroth-order, four-port scattering matrix of a 2D-periodic stack.
solve_crossed_grating –

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

Lattice ¶

Bases: NamedTuple

A 2D Bravais lattice: primitive vectors a1, a2 as (x, y) in µm.

Use the constructors: :meth:rectangular, :meth:square, :meth:hexagonal (the triangular lattice) and :meth:oblique. A (Λx, Λy) tuple passed as a period means :meth:rectangular. Components may be JAX values (to optimise the lattice) except under truncation="circular".

Methods:

hexagonal –

Triangular lattice of nearest-neighbour distance period; a1 along x, rotated.
oblique –

|a1| = a along x and |a2| = b at angle_rad from it.

hexagonal `classmethod` ¶

hexagonal(period, rotation_rad=0.0)

Triangular lattice of nearest-neighbour distance period; a1 along x, rotated.

Source code in hyperbolic_optics/rcwa2d.py

@classmethod
def hexagonal(cls, period, rotation_rad=0.0) -> "Lattice":
    """Triangular lattice of nearest-neighbour distance ``period``; ``a1`` along x, rotated."""
    c, s = jnp.cos(rotation_rad), jnp.sin(rotation_rad)
    h = jnp.sqrt(3.0) / 2

    def turn(x, y):
        return (period * (c * x - s * y), period * (s * x + c * y))

    return cls(turn(1.0, 0.0), turn(0.5, h))

oblique `classmethod` ¶

oblique(a, b, angle_rad)

|a1| = a along x and |a2| = b at angle_rad from it.

Source code in hyperbolic_optics/rcwa2d.py

@classmethod
def oblique(cls, a, b, angle_rad) -> "Lattice":
    """``|a1| = a`` along x and ``|a2| = b`` at ``angle_rad`` from it."""
    return cls((a, 0.0), (b * jnp.cos(angle_rad), b * jnp.sin(angle_rad)))

Shape ¶

Bases: NamedTuple

A polygon of one material within the unit cell.

Attributes:

vertices (Array) –

[..., V, 2] in µm, in cell coordinates (either winding). Shapes may cross the cell edge; they are periodised. They must not overlap one another or their periodic images (overlaps add).
eps (Array) –

Permittivity tensor in the lab frame, [..., 3, 3].
mu (Array | None) –

Permeability tensor, or None for μ = 1.

PatternedLayer ¶

Bases: NamedTuple

One layer: a background material and non-overlapping :class:Shape s.

CrossedGratingResult ¶

Bases: NamedTuple

Amplitudes and efficiencies [..., order, out, in].

orders is [N, 2], a NumPy array (static, so it can index inside jit).

uniform ¶

uniform(eps, thickness_um, mu=None)

An unpatterned layer.

Source code in hyperbolic_optics/rcwa2d.py

def uniform(eps, thickness_um, mu=None) -> PatternedLayer:
    """An unpatterned layer."""
    return PatternedLayer(jnp.asarray(thickness_um), jnp.asarray(eps), (), mu)

rectangle ¶

rectangle(center, size)

Vertices [4, 2] of an axis-aligned rectangle (differentiable in both).

Source code in hyperbolic_optics/rcwa2d.py

def rectangle(center, size) -> jax.Array:
    """Vertices ``[4, 2]`` of an axis-aligned rectangle (differentiable in both)."""
    cx, cy = center
    hx, hy = size[0] / 2, size[1] / 2
    return jnp.stack(
        [
            jnp.stack([cx - hx, cy - hy]),
            jnp.stack([cx + hx, cy - hy]),
            jnp.stack([cx + hx, cy + hy]),
            jnp.stack([cx - hx, cy + hy]),
        ]
    )

ellipse ¶

ellipse(center, radii, n_vertices=64, rotation_rad=0.0)

Vertices [n, 2] of an ellipse, as a polygon of equal area.

The radius is scaled so the polygon's area equals the ellipse's, which removes the leading-order polygon error in the zeroth Fourier coefficient.

Source code in hyperbolic_optics/rcwa2d.py

def ellipse(center, radii, n_vertices: int = 64, rotation_rad=0.0) -> jax.Array:
    """Vertices ``[n, 2]`` of an ellipse, as a polygon of equal area.

    The radius is scaled so the polygon's area equals the ellipse's, which
    removes the leading-order polygon error in the zeroth Fourier coefficient.
    """
    angle = 2 * jnp.pi * jnp.arange(n_vertices) / n_vertices
    scale = jnp.sqrt(2 * jnp.pi / (n_vertices * jnp.sin(2 * jnp.pi / n_vertices)))
    x = radii[0] * scale * jnp.cos(angle)
    y = radii[1] * scale * jnp.sin(angle)
    c, s = jnp.cos(rotation_rad), jnp.sin(rotation_rad)
    return jnp.stack([center[0] + c * x - s * y, center[1] + s * x + c * y], axis=-1)

polygon_fourier ¶

polygon_fourier(vertices, gx, gy)

∫_P exp(−i(gx·x + gy·y)) dA over a polygon, exactly.

By the divergence theorem, (i/|g|²) Σ_edges (g·N_j) ∫_edge exp(−i g·r) dl with N_j the edge's outward normal scaled by its length; each edge integral is closed-form. g = 0 gives the area.

Parameters:

vertices (Array) –

[..., V, 2] (either winding).
gx, gy –

Wavevector components [K].

Returns:

Array –

[..., K].

Source code in hyperbolic_optics/rcwa2d.py

def polygon_fourier(vertices: jax.Array, gx: jax.Array, gy: jax.Array) -> jax.Array:
    """``∫_P exp(−i(gx·x + gy·y)) dA`` over a polygon, exactly.

    By the divergence theorem, ``(i/|g|²) Σ_edges (g·N_j) ∫_edge exp(−i g·r) dl``
    with ``N_j`` the edge's outward normal scaled by its length; each edge
    integral is closed-form. ``g = 0`` gives the area.

    Args:
        vertices: ``[..., V, 2]`` (either winding).
        gx, gy: Wavevector components ``[K]``.

    Returns:
        ``[..., K]``.
    """
    v = vertices
    w = jnp.roll(v, -1, axis=-2)
    e = w - v  # [..., V, 2]
    signed_area = 0.5 * jnp.sum(v[..., 0] * w[..., 1] - w[..., 0] * v[..., 1], axis=-1)
    winding = jnp.sign(signed_area)[..., None]  # +1 counter-clockwise

    phase = v[..., :, 0, None] * gx + v[..., :, 1, None] * gy  # [..., V, K]
    along = e[..., :, 0, None] * gx + e[..., :, 1, None] * gy
    small = jnp.abs(along) < 1e-6
    safe = jnp.where(small, 1.0, along)
    segment = jnp.where(
        small,
        1.0 - 0.5j * along - along**2 / 6.0,
        (1.0 - jnp.exp(-1j * along)) / (1j * safe),
    )
    normal = e[..., :, 1, None] * gx - e[..., :, 0, None] * gy  # g · (e_y, −e_x)
    total = jnp.sum(normal * jnp.exp(-1j * phase) * segment, axis=-2)
    g_sq = gx**2 + gy**2
    zero = g_sq == 0
    transform = 1j * total / jnp.where(zero, 1.0, g_sq)
    return jnp.where(zero, jnp.abs(signed_area)[..., None] + 0j, winding * transform)

crossed_grating_smatrix ¶

crossed_grating_smatrix(
    layers,
    period_um,
    wavelength_um,
    incidence_rad,
    eps_top,
    eps_bottom,
    max_order,
    azimuth_rad=0.0,
    factorization="normal",
    doublings=None,
    grid=None,
    truncation="parallelogram",
    *,
    kx=None,
)

The zeroth-order, four-port scattering matrix of a 2D-periodic stack.

A metasurface as a circuit element: light from either side, in either polarization, into the zeroth order on either side, power-normalised over :data:hyperbolic_optics.bcm.PORTS as :func:hyperbolic_optics.bcm.stack_smatrix is. The higher orders are solved for but are not ports. Unitarity requires lossless media, propagating zeroth orders and evanescent higher orders. Arguments as in :func:solve_crossed_grating, with incidence_rad in the top medium. Alternatively, kx supplies the conserved in-plane wavevector over k₀ and overrides incidence_rad, including evanescent top ports. Evanescent ports retain BCM's complex √kz normalization; |s|² is not power there. A crystal eps_bottom (3×3) makes the bottom ports its power-normalised modes.

Returns:

–

class:hyperbolic_optics.bcm.PlanarSMatrix.

Source code in hyperbolic_optics/rcwa2d.py

def crossed_grating_smatrix(
    layers: list[PatternedLayer],
    period_um,
    wavelength_um: jax.Array,
    incidence_rad: jax.Array,
    eps_top: jax.Array,
    eps_bottom: jax.Array,
    max_order: tuple[int, int],
    azimuth_rad: jax.Array = 0.0,
    factorization: str = "normal",
    doublings: int | None = None,
    grid: tuple[int, int] | None = None,
    truncation: str = "parallelogram",
    *,
    kx: jax.Array | None = None,
):
    """The zeroth-order, four-port scattering matrix of a 2D-periodic stack.

    A metasurface as a circuit element: light from either side, in either
    polarization, into the zeroth order on either side, power-normalised over
    :data:`hyperbolic_optics.bcm.PORTS` as :func:`hyperbolic_optics.bcm.stack_smatrix`
    is. The higher orders are solved for but are not ports. Unitarity requires
    lossless media, propagating zeroth orders and evanescent higher orders.
    Arguments as in :func:`solve_crossed_grating`, with ``incidence_rad`` in the
    top medium. Alternatively, ``kx`` supplies the conserved in-plane wavevector
    over k₀ and overrides ``incidence_rad``, including evanescent top ports.
    Evanescent ports retain BCM's complex √kz normalization; |s|² is not power there.
    A crystal ``eps_bottom`` (3×3) makes the bottom ports its power-normalised modes.

    Returns:
        :class:`hyperbolic_optics.bcm.PlanarSMatrix`.
    """
    _check_x64()
    orders, zeroth, scattering, _, flux_out = _scattering(
        layers, period_um, wavelength_um, incidence_rad, eps_top, eps_bottom, max_order,
        azimuth_rad, factorization, doublings, grid, truncation, kx,
    )  # fmt: skip
    size = orders.shape[0]
    index = jnp.array([zeroth, size + zeroth])  # p and s of the zeroth order

    blocks = tuple(block[..., index, :][..., :, index] for block in scattering)
    k_par = (
        jnp.sqrt(jnp.asarray(eps_top)) * jnp.sin(jnp.asarray(incidence_rad))
        if kx is None
        else jnp.asarray(kx)
    )
    bottom = flux_out[..., index] > 0 if _is_tensor(eps_bottom) else None
    return _normalize_smatrix(blocks, k_par, eps_top, eps_bottom, bottom)

solve_crossed_grating ¶

solve_crossed_grating(
    layers,
    period_um,
    wavelength_um,
    incidence_rad,
    eps_incident,
    eps_exit,
    max_order,
    azimuth_rad=0.0,
    factorization="normal",
    doublings=None,
    grid=None,
    truncation="parallelogram",
)

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

Parameters:

layers (list[PatternedLayer]) –

:class:PatternedLayer s from the incident side down.
period_um (tuple[float, float]) –

A :class:Lattice, or (Λx, Λy) in µm for a rectangular one.
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) –

The exit half-space: a real, positive permittivity, or a 3×3 tensor for a crystal substrate (anisotropic, hyperbolic or lossy, any orientation). Its orders are then the crystal's two transmitted modes, each power-normalised, in place of p and s.
max_order (tuple[int, int]) –

(M₁, M₂): orders |m| ≤ M₁, |n| ≤ M₂ along b₁, b₂ (static).
azimuth_rad (Array, default: 0.0 ) –

Plane of incidence, from the x axis.
factorization (str, default: 'normal' ) –

"normal" (default) or "laurent" (static).
doublings (int | None, default: None ) –

Doublings per layer (static); by default grows with the order.
grid (tuple[int, int] | None, default: None ) –

Sampling (G₁, G₂) of the cell (along a₁, a₂) for the normal-vector method's smooth parts.
truncation (str, default: 'parallelogram' ) –

"parallelogram" (default) or "circular" (static). Circular keeps |G| ≤ max(M₁, M₂)·g_min, with g_min the shortest reciprocal lattice vector: it respects the lattice's symmetry (hexagonal included), and the results do not depend on the primitive vectors chosen to describe it.

Returns:

CrossedGratingResult –

class:CrossedGratingResult: orders [N, 2] and r, t,
CrossedGratingResult –

R, T [..., N, out, in]. Only order (0, 0) is incident.

Source code in hyperbolic_optics/rcwa2d.py

def solve_crossed_grating(
    layers: list[PatternedLayer],
    period_um: tuple[float, float],
    wavelength_um: jax.Array,
    incidence_rad: jax.Array,
    eps_incident: jax.Array,
    eps_exit: jax.Array,
    max_order: tuple[int, int],
    azimuth_rad: jax.Array = 0.0,
    factorization: str = "normal",
    doublings: int | None = None,
    grid: tuple[int, int] | None = None,
    truncation: str = "parallelogram",
) -> CrossedGratingResult:
    """Diffraction of a plane wave by a stack of 2D-periodic anisotropic layers.

    Args:
        layers: :class:`PatternedLayer` s from the incident side down.
        period_um: A :class:`Lattice`, or ``(Λx, Λy)`` in µm for a rectangular one.
        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: The exit half-space: a real, positive permittivity, or a
            3×3 tensor for a crystal substrate (anisotropic, hyperbolic or
            lossy, any orientation). Its orders are then the crystal's two
            transmitted modes, each power-normalised, in place of p and s.
        max_order: ``(M₁, M₂)``: orders ``|m| ≤ M₁``, ``|n| ≤ M₂`` along ``b₁``, ``b₂`` (static).
        azimuth_rad: Plane of incidence, from the x axis.
        factorization: ``"normal"`` (default) or ``"laurent"`` (static).
        doublings: Doublings per layer (static); by default grows with the order.
        grid: Sampling ``(G₁, G₂)`` of the cell (along ``a₁``, ``a₂``) for the
            normal-vector method's smooth parts.
        truncation: ``"parallelogram"`` (default) or ``"circular"`` (static).
            Circular keeps ``|G| ≤ max(M₁, M₂)·g_min``, with ``g_min`` the
            shortest reciprocal lattice vector: it respects the lattice's
            symmetry (hexagonal included), and the results do not depend on the
            primitive vectors chosen to describe it.

    Returns:
        :class:`CrossedGratingResult`: ``orders`` ``[N, 2]`` and ``r``, ``t``,
        ``R``, ``T`` ``[..., N, out, in]``. Only order (0, 0) is incident.
    """
    _check_x64()
    orders, zeroth, scattering, flux_in, flux_out = _scattering(
        layers, period_um, wavelength_um, incidence_rad, eps_incident, eps_exit, max_order,
        azimuth_rad, factorization, doublings, grid, truncation,
    )  # fmt: skip
    r, t, R, T = _order_results(scattering, flux_in, flux_out, zeroth, orders.shape[0])
    return CrossedGratingResult(orders=orders, r=r, t=t, R=R, T=T)

Crossed-Grating RCWA API¶

rcwa2d ¶

Fourier factorisation¶

Lattice ¶

hexagonal classmethod ¶

oblique classmethod ¶

Shape ¶

PatternedLayer ¶

CrossedGratingResult ¶

uniform ¶

rectangle ¶

ellipse ¶

polygon_fourier ¶

crossed_grating_smatrix ¶

solve_crossed_grating ¶

hexagonal `classmethod` ¶

oblique `classmethod` ¶