Why this page exists

This is a single, high-density reference designed to make an AI coding assistant (Claude, ChatGPT, Copilot, …) an expert Ikarus user in one read — every feature, convention and sharp edge, condensed. The same text ships inside the installed package, so any chat can load it without this site:

python -c "import ikarus; print(ikarus.ai_guide())"

It is the single source of truth (the page below is generated from ikarus/AI_GUIDE.md), so it never drifts from the code. Point your assistant here, or paste the command output into a fresh conversation.

Ikarus — Reference for AI Assistants¶

You are working with Ikarus, a 2-D RCWA (Rigorous Coupled-Wave Analysis / Fourier Modal Method) solver for periodic photonics. This page is a dense, authoritative cheat-sheet so you can use it like an expert without rediscovering its API or its sharp edges. When in doubt, prefer these facts over guesses.

PyPI: ikarus-rcwa · import: import ikarus (the names differ on purpose).
Check what you are actually running: a __version__ ending in .devN (e.g. 1.2.0.dev0) means an unreleased checkout, not the published package — it has fixes PyPI does not. python -c "import ikarus; print(ikarus.__file__)" from outside the repo tells you which copy you have (inside it, the working directory shadows the installed package and the answer is always the repo).
Docs: https://cavitytechnologies.github.io/ikarus/ · GitHub: CAVITYtechnologies/ikarus.
Extras: pip install "ikarus-rcwa[inverse]" (pymoo, GA optimize), [grad] (JAX+optax, adjoint optimize + ikarus.grad), [io] (h5py), [progress] (tqdm), [all].

What it is for — and what it is not¶

Frequency-domain electromagnetics of structures periodic in x and y and layered along z: gratings, metasurfaces, photonic-crystal slabs, thin-film stacks, Bragg mirrors, metamirrors. You get per-order efficiencies, complex amplitudes, phase, exit angles and real-space fields.

Not for: isolated (non-periodic) scatterers, continuously-varying-in-z shapes that resist layer slicing, or time-domain/pulse physics. Those are not RCWA's domain. (Anisotropic/birefringent media are supported — see Materials.)

Conventions — the mistakes everyone makes¶

simulate() returns (T, R, result) — transmission FIRST. The single most common bug is unpacking it as (R, T, ...).
SI units, always. Meters for length, degrees for angles. 1550e-9, not 1550.
Time/sign convention is physics exp(-iωt): absorbers have k > 0, Im(ε) > 0. If you feed gain-signed (negative-k) data the energy balance exceeds 1 — that is the tell.
n_orders ≠ resolution. n_orders is the number of Fourier harmonics that resolve the field — the accuracy/cost dial. resolution only draws the geometry (real-space pixels) and is auto-raised to ≥ 4*M+1. You almost always tune only n_orders.
Stack order is cover → … → substrate (sky to ground). The first and last layers must be uniform and semi-infinite (height=np.inf); the cover's index sets the incident wavevector. Interior layers are finite.
Cost is O(M⁶) in 2-D (M = orders per axis). Doubling M ≈ 64× the time. A 1-D grating should use n_orders=(M, 0) with an (Nx, 2) topology — linear, not quadratic. A uniform thin-film stack runs at n_orders=0 (instant).
Energy balance is the universal smoke test. R_total + T_total (result.energy_balance): ≈1 lossless & converged; <1 absorption (physics); slightly >1 unconverged → raise n_orders; wildly >1 numerical breakdown → lower n_orders or raise resolution. Note the default factorization ("auto", the normal-vector method) is not strictly energy-conserving at finite order — a small >1 imbalance on curved/oblique structures is normal and shrinks to 0 as n_orders grows; it is a more honest convergence signal than the classic inverse rule, which can report R_total+T_total==1 while still being a percent off the true answer on curved boundaries.
Factorization is automatic — don't set it. The default factorization="auto" applies the normal-vector (Fast Fourier Factorization) method to every patterned layer: full accuracy on curved/oblique high-contrast boundaries, reducing exactly to the classic inverse rule on axis-aligned geometry. Users never need to choose. "li", "laurent", "normal" remain as explicit overrides for benchmarking.

What the code checks for you¶

Most of the rules above are enforced at construction time, so a slip fails loudly instead of returning a confident wrong number. You still have to get #1 right yourself: the first two elements are both complex scalars of the same kind, so no check can tell a swapped unpack from a correct one.

And note what they are. T and R from the tuple are the zero-order complex amplitude coefficients — the same objects as result.T / result.R, not power. For transmitted/reflected power use result.T_total / result.R_total (or abs(T)**2 for the zero order alone). Writing power = T is wrong and quietly optimistic: on a bare Air/SiO2 interface T = 0.9827 while the actual transmittance is 0.9657, and the two diverge further on lossy or strongly patterned structures.

you write	what happens
`wavelength=1550`, `period_x=500`, `height=200`	warns: looks like nm passed as m (`>= 1 mm`). Legitimate for THz/mm-wave, so it warns rather than raises. `height=np.inf` never warns.
`n_orders=-5`, `resolution=0`	`ValueError`. Note `n_orders=(M, 0)` is valid — that is a 1-D grating.
`n_orders=(25, 25)`, or a bare `n_orders=200`	warns: too expensive. A scalar applies to both axes, so `200` means 160801 harmonics, not 200.
`theta=90` or more	`ValueError` — `theta` is from the normal, so `\\|theta\\| >= 90` is not an incident wave.
non-integer or negative `topology`	`ValueError`. A negative index would silently wrap and select the wrong material; a greyscale map belongs in `ikarus.inverse.pixels`.
cover/substrate not `np.inf`, or `np.inf` in the middle	`ValueError` — stack must be cover → interior → substrate.
unknown material name	`KeyError` listing every built-in.

Warnings are real signals, not noise. If one fires, fix the input rather than suppressing it — each marks a case that otherwise produces a plausible wrong answer.

Never report a number without its confidence¶

The checks above stop bad input. They cannot tell you whether a perfectly valid run is converged. A truncation that is too low returns a clean, confident, entirely wrong efficiency, and nothing about the number itself reveals it.

So when you report any result to the scientist, report these alongside it:

T, R, result = rcwa.simulate()
print(f"R = {result.R_total:.4f}  "
      f"(n_orders={rcwa.n_orders}, energy balance {result.energy_balance:.4f})")

energy_balance ≈ 1 for a lossless stack means converged; < 1 means absorption; > 1 means not converged (see convention #7 for the detail).
Quoting n_orders makes the result reproducible and shows what it cost.
If you changed n_orders and the answer moved, it was not converged — say so, rather than reporting the last value you happened to compute.

A number without this context is not a result, it is a guess with four decimal places. Reporting it bare is the single easiest way to mislead someone who is trusting you to have checked.

Two points can agree by luck. optimize()'s built-in check, and any hand-rolled "compare M against M+4", is a smoke test, not proof. Silence from it does not mean converged — it means those two truncations happened to land close. A field test found a 62-pixel freeform grating whose specular efficiency oscillated between 0.78 and 0.89 out to n_orders=210, never settling; any two-point check that happened to straddle the oscillation would have blessed it. The honest number (78.4 %) only emerged from sweeping eight truncations and requiring the curve to be visibly flat.

So: a two-point check is enough for a simple, well-conditioned structure. A freeform or high-contrast design needs a real ladder — sweep several n_orders and look at the curve (ikarus.tools.convergence_curve). This is convention #7 one level up: energy balance is not convergence, and neither is a single pair of agreeing points.

Minimal forward simulation¶

import numpy as np
from ikarus import RCWA, shapes

rcwa = RCWA(period_x=500e-9, period_y=500e-9, resolution=(96, 96), n_orders=(8, 8))
rcwa.add_uniform_layer(np.inf, "Air")                                   # cover
rcwa.add_layer(220e-9, shapes.circle(radius=0.3, grid_shape=(96, 96)),
               ["Air", "Si"])                                           # patterned
rcwa.add_uniform_layer(np.inf, "SiO2")                                  # substrate

rcwa.set_source(wavelength=600e-9, theta=0, polarization="linear")
T, R, result = rcwa.simulate()                                         # T FIRST
print(result.R_total, result.T_total, result.energy_balance)
print(result.R_phase)                                                  # zero-order phase (rad)

The API surface¶

RCWA(period_x, period_y, resolution=32, n_orders=25, materials=None, convergence_tol=1e-4, factorization="auto") - .add_uniform_layer(height, material, name="") — height=np.inf for cover/substrate. - .add_layer(height, topology, materials, resolution=None, name="") — topology is an integer array, a parametric Shape, or any object exposing an .img array. materials[i] dresses topology index i. - .set_source(**kwargs) — remembers unspecified fields between calls (so a sweep mutates one parameter at a time). First call needs wavelength. - .simulate(auto_converge="never"|"once"|"always") -> (T, R, result). - .get_fields(plane="xy"|"xz"|"yz", ...) → field maps; .visualize_structure(plane="xz", layer_index=None, savefig=...). - .save_results(path, include=("T","R","metadata"[, "fields"]), result=None) and RCWA.load_results(path) (needs [io]).

Source (via set_source): wavelength, theta=0 (from +z, normal), phi=0 (azimuth from +x), polarization="linear"|"RCP"|"LCP", linear_pol_angle=0 (degrees from TE: 0=TE/s, 90=TM/p).

SimulationResult (the 3rd return value): - T, R — zero-order complex amplitudes (a scalar for linear pol, or a {"co","cross"} dict for circular). T_total, R_total — summed over all propagating orders. - T_orders, R_orders — per-order efficiencies; orders = (p, q) integer arrays; order_index(p, q) → the array index of one order. - T_phase, R_phase — zero-order phase in radians (np.angle of the coeff). - theta_out_ref/trn, phi_out_ref/trn — exit angles in degrees (NaN = evanescent). theta_out_* is an unsigned polar magnitude; the ± direction lives in phi_out_* (compare a grating-equation angle against abs(...)). theta_out_trn is the angle inside the substrate, not in air. The substrate is semi-infinite, so there is no back surface to refract at. A grating steering +1 to 50.8° in air reports 32.5 on an n=1.444 substrate — both look plausible. Use res.theta_out_trn_in(n_exit=1.0) rather than doing the Snell yourself: it returns NaN for orders past the critical angle, which are totally internally reflected and never leave the chip. Those read as perfectly ordinary deflection angles in theta_out_trn (e.g. 53.6°), and a hand-rolled arcsin turns them into a domain error or a nonsense number instead of telling you the design is dead. theta_out_ref needs no correction: the cover is the reflected light's medium. - energy_balance — R_total + T_total. solution — raw modal solution for fields.

Materials: a shared default_library ships Air, Ag, Au, GaN, GaP, Si, Si₃N₄, SiO₂, TiO₂, aSi. Anywhere a material is wanted you may pass a name ("Si"), a bare number (1.5, 2.4+0.01j = constant index), a JSON path, or a Material. default_library.get("Si", 1.55e-6) returns the complex index; default_library.available() lists names. Custom materials: default_library.add_from_file(path, name=..., persist=False) imports a CSV (lambda_nm, n[, k]) or JSON (persist=True writes it into the library for future sessions); Material.constant(2.4+0.01j) for non-dispersive; a JSON lorentz block for an oscillator model.

Anisotropic (birefringent) materials: anywhere a material goes you may also pass a 3-tuple (n_x, n_y, n_z) (diagonal tensor; each element any scalar spec, so dispersive components work), an ikarus.AnisotropicMaterial(n_x, n_y, n_z, angle=deg) (in-plane rotation of the principal axes → eps_xy off-diagonals), or ikarus.uniaxial(n_o, n_e, axis="z"|"x"|"y"|angle_deg). Works in uniform and patterned layers, at every factorization. Scope: the in-plane tensor plus a distinct z response — NO tilted optic axis (eps_xz/eps_yz) and NO magneto-optic gyrotropy; the cover and substrate must be isotropic (a clear ValueError otherwise). (n, n, n) reduces exactly to the scalar material. Convention reminder: at normal incidence linear_pol_angle=0 is E along y (sees eps_yy), 90 is E along x (sees eps_xx).

Shapes (from ikarus import shapes): functions in fractional [0,1) unit-cell coords — circle, rectangle, ring, ellipse, cross, polygon, all taking grid_shape=(Nx,Ny); combine, rotate. Parametric classes (ikarus.shapes): Circle, Ellipse, Rectangle, Ring, Cross, SplitRing — any parameter may be free(lo,hi), plus a rotation angle; .to_grid(shape) rasterizes.

Fields and real-space visualization¶

rcwa.get_fields(plane="xz"|"yz"|"xy", nx=…, ny=…, y_position=…, z_positions=[…]) → a dict of FieldMap keyed by plane (xy returns one map per z). Each FieldMap carries .E/.H (complex, last axis = x,y,z components), .intensity (|E|²), .coords, and .eps (the structure on the same grid). Plot with from ikarus.visualization import plot_field: plot_field(fmap, component="intensity"|"Ex".."Hz"|"Eyphase", overlay=True) — a component magnitude, or "<comp>phase" for its phase; overlay outlines the material boundaries. Orientation (a real gotcha): xz/yz cross-sections are drawn stack-vertical — cover on top, light entering from the top (z increases downward); axes auto-scale to nm/µm. rcwa.visualize_structure(plane="xz") draws the layer stack; plane="xy" draws the first patterned layer's topology. Shortcuts: rcwa.visualize_fields(component=…) does get_fields+plot_field in one call; plot_field_xy(field_dict, component=…) tiles an xy multi-z dict.

Sweeps and progress bars¶

from ikarus import Sweep, progress

sw = Sweep(rcwa).over(wavelength=np.linspace(1.4e-6, 1.7e-6, 31)).run(progress=True)
sw.R_total           # ndarray shaped like the axes; also T_total, energy_balance
sw.order(0, 0, "R")  # per-order efficiency across the grid; sw.axes, sw.results

Sweep.over(**axes) sweeps source parameters only (wavelength/theta/ polarization/…); geometry changes rebuild, so keep those in a manual loop wrapped in progress(iterable, enable=True, desc=...). Reuse one RCWA across a sweep for convenience/correctness — not speed: there is no eigenmode caching, so every set_source re-solve is a full solve (wavelength, angle, polarization alike). Budget one solve per sweep point (see the cost table below).

Inverse design (`from ikarus.inverse import …`)¶

Two engines behind ONE call — the user never chooses: - Adjoint (gradient-based), via the optional differentiable JAX solver ([grad] extra): picked automatically for freeform pixels(...) maps and free heights/periods with a single (possibly multi-wavelength / worst-case) target. Scales to thousands of pixel DOFs (the whole gradient costs ~1 extra solve); relax-and-project binarization with a min_feature=<meters> fab filter; the final design is hard-thresholded and re-simulated with the NumPy core at the optimization n_orders (so .achieved is honest for that truncation, but NOT necessarily converged -- optimize() runs a convergence check and warns if the metric still moves with n_orders). Optimize at a modest n_orders, then confirm/report higher with optimize(..., verify_n_orders=16); energy balance ~1 does NOT prove convergence here. Its edge GROWS with DOF count (measured: ~40 pixels -> GA competitive or better; ~2000 pixels -> adjoint ~2x the GA's result in half the time). For steering/deflection objectives use init="random" AND restarts=3+ (multi-modal -- single-seed results vary a lot; not optional). Beware trivial attractors (maximize-total-R is satisfied by a uniform slab) and keep n_orders faithful to the pixel size. Adjoint-only knobs (all defaulted): steps=150, learning_rate=0.05, min_feature, beta=(8, 256), init, restarts. - GA / NSGA-III (gradient-free), via pymoo ([inverse] extra): picked for parametric-Shape DOFs, discrete/anisotropic materials, circular-pol coefficient metrics, and whenever a full Pareto front (>= 2 targets) is wanted. algorithm='adjoint'|'ga'|'nsga2'|'nsga3' forces an engine.

Anything exposing variables() + build(params, n_orders) is GA-optimizable — MetaAtom, Structure, or your own class (adjoint currently supports MetaAtom).

Roll your own differentiable objective ([grad] extra) when optimize doesn't fit your figure of merit: from ikarus.grad import solve mirrors solve_stack (pinned to ~1e-13; gradients cross-checked against FMMax to ~4e-7), so jax.grad of any FoM built on it is the adjoint method — the gradient w.r.t. every pixel costs ~one extra solve. Importing ikarus.grad auto-enables JAX x64 (RCWA needs double precision). Also public: GradSolution, tangent_fields_for, and eig (the differentiable non-Hermitian eigensolver, custom VJP). Building blocks for a hand-rolled topology-opt loop live in ikarus.grad.topology: conic_filter (min-feature), tanh_projection (binarization), beta_schedule.

MetaAtom(period, cover, substrate, polarization="linear", pol_angle=0.0) then .add_pattern(topology, materials, height): - period — a float (square), a (period_x, period_y) tuple (rectangular), or free(lo, hi) (square DOF). height — a float or free(lo, hi). - topology — a fixed integer array, pixels(nx, ny, symmetry=...) (freeform binary map), or a parametric Shape (its free(...) params become DOFs). - pixels symmetry: None, 'mirror_x', 'mirror_y', 'mirror_xy', 'c2', 'c4', 'c4v' (c4/c4v need a square grid). Symmetry cuts the DOF count.

Target — the figure of merit. Metrics: - 'R' / 'T' — efficiency into order (default specular (0,0); order=None → total). - 'r_co' / 't_co' — the complex zero-order coefficient (for linear pol, the co-pol one). - 'r_phase' / 't_phase' — its phase in radians (matched modulo 2π).

Constructors (all take at=/band=, order=(0,0), weight=, worst_case=): - Target.maximize(metric, at=1550e-9) · Target.minimize(metric, band=(300e-9,600e-9,6)) - Target.match(metric, value, at=...) — drive metric to value. - Wavelengths: at=1550e-9, at=[1064e-9, 1550e-9], or band=(lo, hi[, n]); multiple are aggregated by the mean, or the worst case if worst_case=True.

optimize(atom, targets, n_orders=8, algorithm="auto", pop=100, n_gen=60, seed=0, verbose=True, progress=False, verify_n_orders=None, restarts=1, **adjoint_options) -> OptimizeResult. algorithm="auto" picks adjoint or GA per the rules above (one differentiable Target → adjoint; a list → NSGA-III; pop/n_gen are GA-only). Result: .achieved (metric units — the number to quote; .F is the internal minimization loss: for maximize targets F = 1 − achieved), .plot() (one-line convergence curve in metric units, both engines), .algorithm (which engine ran: 'adjoint'/'ga'/'nsga3'), .params (best dict), .metaatom/.rcwa (a ready-to-simulate RCWA), .report(), .X, .F, .history.

Multi-objective / Pareto (a list of ≥2 Targets → NSGA-III): result.multi is True; .X is the Pareto set (a list of param dicts) and .F the objective matrix (designs × objectives). .achieved/.report() give the best of each metric across the front, while .params/.rcwa are just the first front point — materialize any other design i with atom.build(result.X[i], result.n_orders). .plot() is single-objective only; for the trade-off curve scatter np.asarray(result.F) yourself.

Structure — multi-layer / shared-parameter inverse design. Subclass it, declare params as class attributes (free(...) = DOF, plain value = fixed; a period attribute is required), and implement define(self, p) calling self.add_layer(height, topology, materials) (p is the resolved-parameter namespace). cover, substrate, resolution, polarization, pol_angle are configuration, not DOFs. It plugs into optimize() unchanged.

from ikarus.inverse import Structure, free
from ikarus.shapes import Circle

class ARStack(Structure):
    cover, substrate, resolution = "Air", "SiO2", 64
    period = free(0.20e-6, 0.40e-6)
    h1 = free(0.05e-6, 0.30e-6)
    def define(self, p):
        self.add_layer(p.h1, Circle(radius=0.3), ["Air", "Si3N4"])

Which construct? Forward sweep already answers it → use Sweep, skip the optimizer. One patterned layer with a few meaningful knobs → MetaAtom + a parametric Shape. Freeform topology, no shape prior → MetaAtom + pixels. Multiple layers or parameters shared/derived across layers → Structure.

Idiom — efficiency and a target phase (metamirror or metalens). An atom that is both efficient (|r| → 1 or |t| → 1) and imparts a chosen phase is the core metasurface move. Do not fold it into a single unit-modulus match (Target.match("r_co"/"t_co", np.exp(1j*phi))): that minimizes |coeff - e^(i*phi)|, so when the phase is out of reach at high efficiency (the reachable high-|coeff| phase is more than ~90° from phi), the lowest-loss point is |coeff| → 0 -- the optimizer switches the atom off and the run collapses to R ≈ 0 / T ≈ 0. It looks like an optimizer/DoF failure but is really the objective preferring nothing to a wrong-phase design. Keep efficiency and phase as separate objectives so efficiency is never traded to zero. Choose by degrees of freedom:

One knob (e.g. a pillar whose radius sweeps 0..2*pi at high T): skip the optimizer -- Sweep the parameter and look up the geometry per target phase. This is the standard metalens/metamirror-library route and cannot collapse.
A few parametric knobs: pass a Pareto pair and pick the knee -> optimize(atom, [Target.maximize("R", at=lam), Target.match("r_phase", value=phi, at=lam)]) (use T / t_phase for a lens). With efficiency as its own objective, R = 0 sits at a dominated corner of the front and never masks a good design.
Freeform / many pixels (where only the gradient engine scales): write the scalar figure of merit yourself on the differentiable core, loss = (1 - R) + w*|wrap(arg(r) - phi)|/pi (phase weight w), via ikarus.grad.solve + jax.grad. R = 0 costs a full 1.0 there and can never win.

Performance and accuracy¶

Budget by this rough 2-D cost table (one 96² patterned layer, 1 BLAS thread): n_orders=(4,4)≈0.06 s, (6,6)≈0.4 s, (9,9)≈3 s, (12,12)≈15 s — steeply O(M⁶). A 40×60 design map at (9,9) is ~2 h; coarsen n_orders/grid while exploring, refine once. (A 1-D (M,0) grating is far cheaper — linear.)
Pin BLAS to one thread for GA/sweep loops (many small solves): set OMP_NUM_THREADS/OPENBLAS_NUM_THREADS/MKL_NUM_THREADS/VECLIB_MAXIMUM_THREADS to "1" before importing numpy — often ~10× faster on many-core machines. For one big solve (M ≳ 20) do the opposite: let BLAS thread.
Converge n_orders at the worst-case wavelength/polarization. Prefer simulate(auto_converge="once") — it raises n_orders until the complex zero-order R/T (magnitude and phase) stop changing, and caches the result. For a one-off solve, simulate(check_convergence=True) re-solves once higher and warns if it's still moving. Do not trust R+T≈1 as convergence — a lossless structure conserves energy at every n_orders while R/phase still drift (the classic high-contrast-TM trap that has cost real optimization runs). Manual study: from ikarus.tools.convergence import convergence_curve; orders, vals = convergence_curve(rcwa, range(4,21,2), metric="R"|"R_phase"|"T_phase") — its phase metrics are in DEGREES, unlike result.T_phase/R_phase and Target.match("r_phase", ...) which are radians, so a tolerance carried across without converting is off by 57.3x. An unrecognized metric raises and lists the valid names. (it restores your n_orders afterward).
Fourier factorization: the default factorization="auto" applies the normal-vector (Fast Fourier Factorization) method, giving fast TM / high-contrast convergence on sharp and curved boundaries — high-contrast TM gratings converge by n_orders≈10–15 instead of drifting. On axis-aligned geometry it reduces exactly to Li's inverse rule ("li", the default through v0.7). It works automatically for any topology and any number of materials (it acts on the rendered ε(x,y) grid); "li", "laurent", "normal" are explicit overrides for benchmarking. Caveat: R+T≈1 does not prove convergence for high-contrast TM — always check that R/phase have stopped moving with n_orders.

Known gaps (do not promise these)¶

Tilted-optic-axis anisotropy (eps_xz/eps_yz) and magneto-optic gyrotropy (eps_xy != eps_yx) — anisotropy covers the in-plane tensor + z only, and the cover/substrate must be isotropic.
No GPU for the NumPy core (the optional ikarus.grad JAX solver runs on GPU where JAX does; the default engine is CPU NumPy/SciPy).

Runnable examples¶

Canonical, self-contained scripts ship in the package — run any with python -m ikarus.examples.<name>, or read the source for a copy-paste idiom: feature_tour (a guided tour of most features), validation_fresnel (checks Ikarus against analytic Fresnel — the credibility demo for a skeptic), grating_diffraction (1-D efficiencies vs wavelength), metasurface_spectrum (resonant dielectric-pillar spectrum), inverse_metamirror (inverse design in a few lines), save_load (HDF5 round-trip). A fuller gallery is in the docs.