Binary Grating — where the diffraction orders go

Download .py All tutorials

A plane wave at normal incidence hits a periodic slab — a high-index ridge filling half of each unit cell. The grating equation decides which orders escape as propagating waves:

\[\sin\theta_m = \sin\theta_i + m\,\frac{\lambda}{\Lambda} \qquad\Longrightarrow\qquad \text{order } m \text{ propagates} \iff \left|m\,\frac{\lambda}{\Lambda}\right| \le 1 .\]

Everything else is evanescent. So sweeping the period $\Lambda$ through $\lambda$ switches the $\pm1$ orders on at $\Lambda=\lambda$ — the Rayleigh anomaly. Below it the structure is sub-wavelength: only $m=0$ survives and the slab behaves as a homogeneous effective medium, which is the regime metasurfaces work in. Above it the ±1 orders open, and the very same slab sends most of what it transmits into them.

The term list is the whole problem

It is the scalar-Helmholtz weak form you would hand jno.fem — nothing RCWA-specific in it:

sol = jno.rcwa([
    ui.x * vi.x + ui.y * vi.y + ui.z * vi.z - K0**2 * eps * (u * vi),  # ∇u·∇v − k₀²·ε·u·v
    -(1j * K0 * ut) * vt,                    # outgoing radiation (top ambient)
    -(1j * K0 * ub - 2j * K0) * vb,          # incident plane wave + radiation (bottom)
    ul - ur,                                 # Floquet periodicity in x
    uf - ubk,                                # Floquet periodicity in y
], orders=160, grid=256).solve()

jno.rcwa infers the period from the Floquet ties, the ambients from the two z-normal radiation faces, the layer stack and $\varepsilon$ from the volume coefficient, and the wavelength from $\varepsilon$ in the vacuum superstrate. It infers the polarization from the equation: the ridges run along $y$, and for such a grating the scalar equation is exactly Maxwell's for $E_y$, so it is solved TE — the same equation jno.fem would solve (see Polarization). Only the truncation is genuinely yours: orders Fourier terms and a grid for sampling $\varepsilon$. At these values the reflectance at $\Lambda=1.5$ is within $2\times10^{-3}$ of an independent RCWA with exact Fourier coefficients (0.8441). Inspect what was inferred with rc.spec, which works without the [rcwa] backend:

RcwaSpec(period=(1.5, 0.3), layers=3, wavelength=1.0, k_in=(0.0, 0.0), source_face='bottom')
rc.spec.polarization       # (0j, (1+0j)) -- E along y
rc.spec.polarization_from  # 'y: every layer is invariant along y and the incidence has no y component, ...'

The mesh does not have to resolve the pattern

Because eps is an analytic jno.fn, it is evaluated directly on the RCWA grid rather than interpolated off the mesh. The tetrahedral mesh only carries the tags and periodicity, so size=0.25 is plenty — the grating profile stays exact.

Two regimes, one structure

Λ = 0.6  (λ/Λ = 1.667, sub-wavelength) — grating equation allows m ∈ [0]
  R = 0.00187   T = 0.99813   R + T = 1.000000
  T(-1) = 0.00000   T(+0) = 0.99813   T(+1) = 0.00000

Λ = 1.5  (λ/Λ = 0.667, diffractive)    — grating equation allows m ∈ [-1, 0, 1]
  R = 0.84594   T = 0.15406   R + T = 1.000000
  T(-1) = 0.06879   T(+0) = 0.01648   T(+1) = 0.06879

At $\Lambda=0.6$ the slab is nearly transparent: as an effective medium it is close to a half-wave layer, and only the straight-through order exists. $\Lambda=1.5$ happens to sit on the flank of a narrow reflection resonance of this slab — an independent RCWA puts $R$ at 0.35 for $\Lambda=1.44$ and 0.87 for $1.49$ — so there it reflects 85 %, and of what it transmits, 89 % goes into $\pm1$. At $\Lambda=0.6$ those orders are exactly dark — not small, zero — because they are not propagating solutions at all.

What is verified

There is no closed-form field here, but four statements are exact and are asserted in the script:

Check	Holds to
Energy conservation $R+T=1$ (lossless)	$10^{-15}$
Grating-equation cutoff: orders outside $\lvert m\lambda/\Lambda\rvert\le1$ carry zero power	exactly `0.0`
Mirror symmetry $T(+1)=T(-1)$ at normal incidence	$\approx 10^{-10}$
Completeness: propagating orders sum to $T$	$10^{-6}$

The $\pm1$ symmetry floor does not improve with orders

$\lvert T(+1)-T(-1)\rvert$ measures $7.5\times10^{-11}$ at orders=160 and $6.9\times10^{-11}$ at orders=320 — it is not Fourier truncation. The $\pm1$ pair is degenerate at normal incidence and the modal eigensolve splits it at its own floor. Energy conservation, computed from the same solution, holds to machine precision. Assert against the measured floor, not against zero.

$\Lambda = m\lambda$ is singular

There order $m$ is exactly grazing — a Wood–Rayleigh anomaly — and the modal problem is degenerate. jno.rcwa raises with the period and wavelength named rather than returning a NaN efficiency, so the sweep below deliberately steps around $\Lambda=1$ and $\Lambda=2$.

Result

$Diffraction efficiency against grating period. Through the shaded sub-wavelength region the ±1 orders are flat at exactly zero; at period equal to wavelength they switch on and take most of the transmitted power, while R+T stays pinned at one across the whole sweep.$

The $\pm1$ curve is flat at exactly zero through the shaded region and turns on at $\Lambda=\lambda$, matching the cutoff to the sweep's resolution; $R+T$ stays at $1$ across all 13 periods. Beyond the cutoff the $\pm1$ orders carry up to about 0.6 of the incident power; the sweep's coarse steps pass over the resonance near $\Lambda=1.5$ described above.

Going further

Every efficiency is a differentiable JAX scalar, so this is already an inverse-design objective — put a jno.np.parameter in eps and jax.grad a target order fraction through the modal solve. See RCWA for anisotropic $\hat\varepsilon$, in-plane PML, internal sources, sol.jones() for polarization, and sol.aerial() for imaging.

Full script

"""Binary dielectric grating: where the diffraction orders go, via ``jno.rcwa``.

A plane wave at normal incidence hits a periodic slab -- a high-index ridge filling half of each
unit cell. The **grating equation** decides which orders escape as propagating waves,

    sin(θ_m) = sin(θ_i) + m·λ/Λ ,   so order m propagates  ⟺  |m·λ/Λ| ≤ 1  (θ_i = 0 here),

and everything else is evanescent. Sweeping the period Λ through λ therefore switches the ±1
orders on at Λ = λ -- the **Rayleigh anomaly** (Rayleigh, *Proc. R. Soc. A* 79, 399, 1907). Below
it the structure is *sub-wavelength*: only m=0 survives, all the transmitted light goes straight
ahead, and the slab acts as a homogeneous effective medium -- the regime metasurfaces work in.
Above it the same slab sends what it transmits into the ±1 orders.

The whole problem is the scalar-Helmholtz term list you would hand ``jno.fem``:

    ∇u·∇v − k₀²·ε·u·v = 0        volume (ε carries the pattern)
    −i k₀ u·v                     outgoing radiation, top ambient
    −(i k₀ u − 2i k₀)·v           incident plane wave + radiation, bottom ambient
    u(left) − u(right) = 0        Floquet periodicity in x
    u(front) − u(back) = 0        Floquet periodicity in y

``jno.rcwa`` infers the period from the Floquet ties, the ambients from the two z-normal radiation
faces, the layer stack and ε from the volume coefficient, the wavelength from ε in the vacuum
superstrate, and the polarization from the equation itself: for a grating whose ridges run along y the
scalar equation is exactly Maxwell's for E_y, so it is solved TE -- the same equation ``jno.fem`` solves.
Only the truncation is genuinely ours to choose: ``orders`` Fourier terms and a ``grid`` for sampling ε,
which converge the reflectance here to about 2e-3 of an exact-coefficient reference (0.8441 at Λ = 1.5). Because ε
is an analytic ``jno.fn``, it is sampled directly on the RCWA grid, so the tetrahedral mesh only
has to carry the tags and can stay coarse. The modal method itself is Moharam & Gaylord,
*J. Opt. Soc. Am.* 71, 811 (1981).

What is verified (no analytic field solution needed -- these are exact statements):
  1. energy conservation, R + T = 1, on a lossless structure;
  2. the grating-equation cutoff -- orders outside |m·λ/Λ| ≤ 1 carry exactly zero power;
  3. mirror symmetry -- T(+1) = T(−1) for a symmetric grating at normal incidence (to ~1e-10: the
     ±1 pair is degenerate here, and the modal eigensolve splits it at a floor that more Fourier
     orders do not lower -- unlike the energy balance, which holds to machine precision);
  4. completeness -- the propagating order efficiencies sum to the total transmission.

Λ = m·λ is excluded throughout: there order m is exactly grazing (a Wood-Rayleigh anomaly), the
modal problem is singular, and ``jno.rcwa`` raises rather than handing back a NaN efficiency.
"""

import jax

jax.config.update("jax_enable_x64", True)  # modal eigen-decomposition of a high-contrast stack

import jax.numpy as jnp  # noqa: E402
import numpy as np  # noqa: E402

import jno  # noqa: E402

WL = 1.0  # vacuum wavelength -- the yardstick every length below is measured against
EPS_RIDGE = 11.0  # ridge permittivity (n ≈ 3.32, silicon-like)
Z0, Z1 = 0.4, 0.6  # slab occupies z ∈ (0.4, 0.6); vacuum ambients above and below
PY = 0.3  # y-period, deliberately sub-wavelength -> a 1-D grating
DUTY = 0.5  # ridge fills the left half of each cell
K0 = 2 * jnp.pi / WL
ORDERS, GRID = 160, 256


def grating(period, orders=ORDERS, grid=GRID):
    """Solve one binary grating and return its order-resolved efficiencies."""
    d = jno.shape.box(0, 0, 0, period, PY, 1.0, size=0.25).domain()  # coarse: ε is analytic
    e = 1e-6
    d.tag("bottom", lambda x, y, z: z < e)
    d.tag("top", lambda x, y, z: z > 1.0 - e)
    d.tag("left", lambda x, y, z: x < e)
    d.tag("right", lambda x, y, z: x > period - e)
    d.tag("front", lambda x, y, z: y < e)
    d.tag("back", lambda x, y, z: y > PY - e)

    u, v = d.fem_symbols()
    xi, yi, zi, _ = d.variable("interior", split=True)
    ui, vi = u.bind(x=xi, y=yi, z=zi), v.bind(x=xi, y=yi, z=zi)

    def on(tag):  # bind u, v on a named face
        s = d.variable(tag, split=True)
        return u.bind(x=s[0], y=s[1], z=s[2]), v.bind(x=s[0], y=s[1], z=s[2])

    (ut, vt), (ub, vb) = on("top"), on("bottom")
    ul, ur, uf, ubk = on("left")[0], on("right")[0], on("front")[0], on("back")[0]

    ridge = jno.fn(lambda x, y, z: jnp.where((Z0 < z) & (z < Z1) & (x < DUTY * period), 1.0, 0.0), [xi, yi, zi])
    eps = 1.0 + (EPS_RIDGE - 1.0) * ridge  # 1 in the ambients, EPS_RIDGE in the ridge

    sol = jno.rcwa(
        [
            ui.x * vi.x + ui.y * vi.y + ui.z * vi.z - K0**2 * eps * (u * vi),  # ∇u·∇v − k₀²εuv
            -(1j * K0 * ut) * vt,  # outgoing (top)
            -(1j * K0 * ub - 2j * K0) * vb,  # incident + outgoing (bottom)
            ul - ur,  # Floquet, x
            uf - ubk,  # Floquet, y
        ],
        orders=orders,
        grid=grid,
    ).solve()

    real = lambda a: float(jnp.real(jnp.asarray(a)))  # noqa: E731
    return dict(
        T=real(sol.efficiency("T")),
        R=real(sol.efficiency("R")),
        orders={m: real(sol.order(m, 0)) for m in (-2, -1, 0, 1, 2)},
    )


def propagating(period):
    """The transmitted orders the grating equation lets escape at normal incidence."""
    return [m for m in (-2, -1, 0, 1, 2) if abs(m * WL / period) <= 1.0]


# --- two regimes, side by side ------------------------------------------------------------------
SUB, DIF = 0.6, 1.5  # sub-wavelength (Λ < λ) and diffractive (Λ > λ)
res = {P: grating(P) for P in (SUB, DIF)}

print("\nBinary dielectric grating -- where the orders go   (λ = 1, ε_ridge = 11, 50% duty)")
for P, r in res.items():
    prop = propagating(P)
    tag = "sub-wavelength" if P < WL else "diffractive"
    print(f"\n  Λ = {P}  (λ/Λ = {WL / P:.3f}, {tag}) -- grating equation allows m ∈ {prop}")
    print(f"    R = {r['R']:.5f}   T = {r['T']:.5f}   R + T = {r['R'] + r['T']:.6f}")
    print("    " + "   ".join(f"T({m:+d}) = {r['orders'][m]:.5f}" for m in (-1, 0, 1)))

for P, r in res.items():
    prop, o = propagating(P), r["orders"]
    assert abs(r["R"] + r["T"] - 1.0) < 1e-6, f"lossless structure must conserve energy (Λ={P})"
    for m in o:
        if m not in prop:
            assert o[m] < 1e-9, f"order {m} is evanescent at Λ={P} but carries {o[m]:.2e}"
    # The ±1 pair is degenerate at normal incidence, so the modal eigensolve splits it at its own
    # floor: measured |T(+1) − T(−1)| ≈ 7e-11, and it does NOT shrink with `orders` (7.5e-11 at 160,
    # 6.9e-11 at 320) — it is not Fourier truncation. Energy conservation, by contrast, holds to 1e-15.
    assert abs(o[+1] - o[-1]) < 1e-6, f"symmetric grating at normal incidence: T(+1) must equal T(-1) (Λ={P})"
    assert abs(sum(o[m] for m in prop) - r["T"]) < 1e-6, f"propagating orders must sum to T (Λ={P})"

# the physical punchline, stated as a test rather than prose
assert res[SUB]["orders"][+1] == 0.0, "below the Rayleigh cutoff the ±1 orders must be exactly dark"
assert res[DIF]["orders"][+1] + res[DIF]["orders"][-1] > 0.8 * res[DIF]["T"], (
    "above the cutoff this grating should throw most of its light into ±1"
)

# --- sweep the period across the Rayleigh anomaly at Λ = λ ---------------------------------------
# Λ = m·λ puts order m exactly grazing (a Wood-Rayleigh anomaly): the modal problem is singular
# there and `jno.rcwa` raises rather than returning a NaN efficiency. Step around Λ = 1 and Λ = 2.
periods = np.concatenate([np.linspace(0.55, 0.95, 5), np.linspace(1.05, 1.95, 8)])
sweep = [grating(float(P)) for P in periods]
print(f"\n  swept {len(periods)} periods across the Rayleigh anomaly at Λ = λ = {WL}")

first_on = next(P for P, r in zip(periods, sweep) if r["orders"][+1] > 1e-6)
assert first_on > WL, f"the ±1 orders must not switch on below Λ = λ (got {first_on})"
assert all(abs(r["R"] + r["T"] - 1.0) < 1e-6 for r in sweep), "energy conserved across the whole sweep"
print(f"  ±1 orders first carry power at Λ = {first_on:.2f}  (cutoff is Λ = λ = {WL})")

Check	Holds to
Energy conservation \(R+T=1\) (lossless)	\(10^{-15}\)
Grating-equation cutoff: orders outside \(\lvert m\lambda/\Lambda\rvert\le1\) carry zero power	exactly `0.0`
Mirror symmetry \(T(+1)=T(-1)\) at normal incidence	\(\approx 10^{-10}\)
Completeness: propagating orders sum to \(T\)	\(10^{-6}\)