solprop architecture

This is the current architecture of solprop, the reference that changes to the design are checked against. Companion documents: conventions.md (units, signs, shapes), stepper.md (the solver), backends.md (execution), kernel-contract.md (the seam), decisions.md (why) and roadmap.md (what is next). What this page describes in prose, structure.md draws: six figures of the layers, the data flow, the contract, one stage, the backends and the output side.

1. What solprop is

solprop models ultrafast nonlinear pulse propagation: the unidirectional pulse propagation equation (UPPE) for the full electric field or for a complex envelope, in a single spatial mode, in many modes of a waveguide, and in free space (radially symmetric and full 3-D). It began as a Python rewrite of Luna.jl: it keeps Luna’s scope and its most valuable design choices, and is designed around the following (the user-facing differences are listed in Luna and solprop):

• exponential treatment of the represented diagonal linear flow, with adaptive estimates for the integrated source (stepper.md);

• absorbing boundaries defined as rates, with explicit and split-step treatments;

• a specified one-photon input-noise model normalised for dispersive media;

• field and envelope propagation with explicit supported combinations and limitations;

• polarisation as an ordinary array axis, at no cost when unused;

• a hot loop designed for threads and GPUs from the start, behind one seam;

• lean memory: bounded buffers, node chunking, no materialised stages;

• an in-process batch mode for many small or several large problems;

• a clean, fast-loading Python package that plays well with NumPy, SciPy, matplotlib and h5py.

Priorities, in order: accuracy and reliability; then as much CPU performance as can be had; then GPU performance (ADR-020). Double precision by default: realistic modelling always carries noise, so the stepper’s job is efficiency at sufficient accuracy, not extended precision. Single precision is an explicit per-run choice, precision="single", under SEP 33’s contract (ADR-054): the state, stages, transforms and responses in binary32, the controller, host scalars and propagator exponents in binary64, every narrowed table with its own power-of-two exponent and the state at an exact scale. numpy_ref’s single-precision path is its oracle; numba_cpu and native compute in binary64 only.

2. The model

In the frequency domain the propagation equation for every mode m (and polarisation component p) is

∂z Ẽ_m(ω, z) = D_m(ω, z) Ẽ_m(ω, z) + N_m[Ẽ](ω, z)

• D_m(ω, z) = −i(β_m(ω, z) − β₀(z) − β₁(z)(ω − ω₀)) − α_m(ω, z)/2 − α_abs(ω)/2 is the linear operator: dispersion and loss of the mode in a co-moving frame, plus the spectral absorbing collar as a rate. It is diagonal in ω, which is the one property the whole solver design rests on: D(ω, z₁) and D(ω, z₂) commute, so exp(∫_{z₁}^{z₂} D dz') is the linear flow, exactly. Which frame terms apply depends on the grid: an envelope grid removes β₀ unless a third-harmonic term is present; a field grid keeps it, because a constant phase on the half-spectrum of a real field rotates its carrier-envelope phase, and measures the β₁ term from zero frequency (docs/guide/numerics/linear-operator.md, “Which terms the frame removes”).

• N_m is the nonlinear term: the projection onto mode m of the nonlinear polarisation evaluated in the time domain at points of the fibre cross-section (Kerr, Raman, plasma, and later χ⁽²⁾ and spectral-domain responses such as two-photon absorption), multiplied by the UPPE prefactor −iω/(2ε₀c Re n_eff(ω)) (field) or its envelope equivalent, with the temporal absorbing collar as a rate. Only the real part of the index enters the prefactor: the loss belongs to D, and carrying it here would give the nonlinear term a spurious gain component.

Two grids: a field grid (real E(t), half-spectrum state, r2c/c2r transforms, oversampled fine grid for the cubic response) and an envelope grid (complex A(t) about ω₀, full-spectrum state, c2c transforms, optional oversampling for third-harmonic terms). Both are first class, with parity tests for shared physics on narrowband pulses. Supported combinations are explicit: plasma supports scalar and Cartesian vector fields on both grids, using a separately derived cycle-averaged envelope model. Explicitly declared tensor Raman supports both grids on NumPy, Numba CPU and native; unlabelled custom kernels and the default silica model remain scalar-only. Envelope plasma with THG is refused. The native CPU/CUDA/HIP paths have their own qualification. The running guide and backends.md record the qualified scope.

Optional atomic response (SEP 12). A scalar helium TDSEGasResponse owns the gas nonlinearity. Its immutable specification becomes a TDSEDescriptor; the driver lazily owns native plans and stable radial chunks and passes a narrow host binding explicitly to NumPy. The optical march uses Adams. Atomic device execution is independent of the optical backend. Numba and the native backend refuse this response by name. The kernel schema includes separate polarisation and current factors: current enters with -1/(2 epsilon_0 c Re(n_eff)), without division by frequency. The provider’s raw discrete shadows are subtracted before density and calibration. Native retention subtracts only the first-order current and applies its explicitly declared positive scale to the whole nonlinear current; its polarisation contribution is zero. Its source is always gated and edge-free (SEP 26): the backend multiplies the already subtracted remainder by its gate, grids.t_window_fine or the atomic window’s, before source projection/FFT, and subtracts the gate’s edge term F[g' S]/(i omega) on the active bins. The raw provider result, atomic drive and separate in-stage field absorber are unchanged. This is a recorded finite-window approximation, not a model of absorption, walk-off or phase mismatch. Cubic-replacement policies also subtract the third-order current and supply their separately specified polarisation; their sources are never gated.

The TDSE binding uses the shared Adams hooks and the step and plane policies of the optical march. Every descriptor has one layout, solprop.tdse-descriptor.v5, which binds the closure, model, window and grid collar geometry into its identity. The saved output records the absolute launch, the TDSE identity, the outcome and any terminal checkpoint beside the step limits and the dense and step-plane policies. The full-band sensitivity to the window and the atomic cut error are unresolved (TDSE accuracy and limitations). The current versions are kernel schema 36 and saved-output schema 25 (kernel contract).

haCC helium (SEP 20). The same route carries a second atomic model, the two-electron haCC atom (TDSEGasResponse.model="hacc" with a HaccBasis), at a fixed current scale of one. Its descriptor may carry an atomic window whose gate replaces the grid collar as the source gate; the binding slices and gates the drive to the window and pads the currents back to the full lattice, so the seam’s arrays keep their shapes. The weak-field consistency check lives in the driver layer (solprop.driver.tdse_consistency), runs before the first propagation call on the run’s open binding, and is cached on disk; its policy never enters an identity. An optional raw response archive (output/tdse_archive.py) records every provider evaluation in a separate append-only HDF5 file: the driver opens its session before any atomic work, the binding writes each evaluation’s input before its first provider call and each raw result before validation, and the reference backend reports each Adams trial’s outcome to the binding from its step record. It changes no numerical step and adds no call; the optical file records the archive’s run UUID and relative path.

Plan lifetime, allocation admission and attempted-call budgets belong to the driver/binding; modal synthesis, P/J projection and spectral assembly belong to the reference backend. See the TDSE manual.

Geometry: mode-averaged propagation is the single-node case of multimode propagation (the effective-area normalisation replaces the quadrature); the GNLSE is the mode-averaged envelope equation with a dispersion-specified mode (β from Taylor coefficients or a tabulated β(ω), γ or n₂ + A_eff, Raman fraction), so it is an ordinary problem rather than a hack. Free-space propagation (radial and 3-D) uses the same equation with a transverse spectral axis in place of the mode index; nothing in the layering below prevents it, and the kernel layout carries the transverse columns in its mode axis (see §12).

3. Layers and their rules

layer	contents	may import	performance rules
P physics / convenience	grids, materials, modes, profiles, pulses, noise, linear-operator assembly, response definitions, ionisation and Raman models, boundaries, normalisation, tools	numpy, scipy, standard library, other P modules	none: clarity first; never in a hot loop
T translation	turns P objects into validated kernel-ready arrays and the KernelProblem contract	P, numpy, scipy	may be expensive once per problem; results cached and validated
B backends	state buffers, transform plans, kernels, the inner z-loop, statistics kernels	numpy, numba, pyfftw; the native package in native/ (Kokkos, Kokkos-FFT and nanobind); from P only the types and constants the seam names — never a physics object	all hot-path rules of backends.md; in-place and out= allowed here only
D driver	Simulation, run loop over save planes, batch runners, prop_capillary/prop_gnlse, output handlers, statistics specs	P, T, B (through the seam only), h5py	none
X processing / plotting	analysis of outputs, standard plots, design tools	numpy, scipy, matplotlib, D (outputs), T (reference tables and their identities, for processing.modal; the recorded plane measure, for processing.mixtures and processing.planes) — never B	none; may call backend statistics kernels through the D layer when a computation belongs on the device

Transverse rules are built in P. modes/quadrature.py constructs its rules numerically – Laurie’s Kronrod recursion, and the field-scale rules of SEP 19 by the Lanczos process and the Golub–Welsch eigenvalue method, certified on construction – once per problem, and caches them. It is the one numerical construction kept in P: it yields nodes and weights for T to fold into the transform, and no compiled code, per-stage evaluation or hot loop enters P with it.

The X → T edges. solprop.processing.modal imports solprop.translate.modal_tables.ModalTransformTables and solprop.translate.general.guide_identity, and that is inherent to the output schema rather than incidental. A general modal file records what its E/H reference is – the table identity, the guide identity – and not the reference’s tens of megabytes, so the reference is rebuilt in T by general_modal_route and read in X, which checks it against the stored identity, path by path, before it synthesises anything. X takes the record type and the identity functions, never an execution map, never a KernelProblem and never a backend. Two further modules read what a file records the same way. processing.mixtures and processing.planes rebuild a spatial plane’s measure from the record with translate.flux.flux_from_record and translate.spatial.axial_admittance (SEP 34), and processing.planes reuses processing.modal’s checked tables and evaluates a general basis Gram with translate.coupled.metric_at. These edges are deliberate (SEP 31): the measure a file records is rebuilt by the code that built it for the run, not by a second copy in X. No other module of X imports T.

Coupling crosses the layers without breaking them (SEP 31). The field plane, the input port, the remap and couple are P; the driver’s Simulation.input_port() asks T (translate.coupling) to fill the port from the translated problem – its metric, reference index, the scalar profiles of its modes (built in P from the modes, coupling.models.ScalarModeBasis) or a general guide’s E/H and contrast – and extraction is X (processing.field_plane, which builds the same profiles from a source= simulation’s modes), because P must not import X. A coupling computed in P therefore never imports T or reads a file. The driver’s incident-beam launch calls couple and installs its launch through with_input.

Rules: imports flow P → T → B; B never builds a physics object — it does not know how to make a mode or a gas, and imports P only for the grid, mode and profile types its statistics are handed at construction, the constants, and the Raman drive form; X never imports B directly. Every public module declares __all__. No utils modules: organise by domain concept. Hand-rolled numerics live in T (translate.chebyshev) and B only; everything else uses NumPy/SciPy/h5py/matplotlib. A value type in P may carry what T designed and apply it with NumPy alone: the plasma column’s nonlinear.ionisation.RunningIntegral holds taps that translate.plasma_integration fits and certifies, and filters and sums with them (SEP 18).

4. Package layout

The public package is organised by domain. structure.md figure 1 lists the modules individually and is checked against src/solprop by tests/test_structure_figures.py.

Location under src/solprop/	Responsibility	Layer
__init__.py, exceptions.py, constants.py	Lazy public imports, package logging, exception hierarchy and SI constants	package / P
grid.py, fourier.py, hankel.py, spatial.py, adiabatic.py	Temporal and spatial grids (radial, the Cartesian plane, and LineGrid, its x axis uniform in y per unit width), Fourier conventions, raw Hankel transforms, spatial media, the axial-flux measure of a plane (SpatialFlux, SEP 34) and the adiabatic path certificate of a z-varying medium	P
profiles.py, gasflow.py, pulses.py, noise.py, polarisation.py, optics.py, tools.py	Longitudinal paths, steady gas flow in a hollow fibre (the taper- and temperature-aware, choked, regime-checked fill, its flow report and calculators), temporal/spatial inputs, noise, polarisation algebra, linear optical elements as transfer functions, and the soliton and phase-matching design calculators	P
coupling/	Coupling simulations through a coherent field at a reference plane (SEP 31): the field plane with its reference, geometry, metric, basis identity and provenance (plane.py), the input port of a receiving simulation and its bound launch (port.py), temporal remapping (remap.py), the coupling models, their scalar mode profiles and the recorded pair envelopes (models.py, envelopes.py), couple with its admission budget and report (coupler.py), and the synthesis of a guide’s exit field (outcouple.py), and the transport core: the linear elements between planes (elements.py) and the ordered transport of a plane with its ledger (transport.py); NumPy, SciPy and other P modules only	P
physdata/	Species, mixtures, equations of state, transport properties (viscosity, sound speed, mean free path), optical material data, metals, nonlinear and Raman data	P
modes/	Mode protocols and derived quantities; capillary, antiresonant, tabulated, step-index, rectangular and planar (slab, uniform in y, per unit width) modes, with PlanarGrid for the planar guide resolved in y; shared gas fill, quadrature and collections; full-vector modal records, material laws and regions, bound step-index and radial graded-index providers, and the guide that pairs a provider with a mode collection, and the connection of a basis that moves along z	P
linear.py, boundaries.py, nonlinear/	Linear assembly, absorbing rates, response definitions, ionisation, Raman, normalisation and the GNLSE mapping	P
translate/	Validated KernelProblem; parameter expansions and z tables; modal/radial/Cartesian transforms and the hybrid planar one (slab modes in x × FFT in y, translate/planar.py); dense tracked modal tables, the general modal route and the coupled-flow and metric records (translate/general.py, translate/coupled.py); the spatial plane measure (translate/flux.py, the FluxPlan of schema 35); the persisted frame integral and the filled input port of a problem’s entrance (translate/coupling.py); response tables, masks and stage bundles	T
backends/	Registry and protocols, shared step and batch drivers, scalar bodies, spatial and population statistics	B
backends/numpy_ref/, backends/numba_cpu/	Reference and compiled CPU execution	B
stepper/, log.py, stats.py	Tableaux, controllers, schedules, progress reporters and statistics specifications	D
driver/, driver/executors/, scan.py	Simulation and convenience interfaces, scan/batch orchestration, Local/Pool/Queue/Slurm execution and scan CLI	D
output/	Memory and HDF5 outputs, provenance/resume, scan storage and artifacts	D
processing/, plotting/, sol.py	Field, spectrum, time, coherence, compression, mixture and scan analysis; a saved plane read back as a field plane (processing/planes.py); propagation/statistics plots and spectrograms; the interactive session namespace, the one module whose import has side effects	X

solprop-native is the backend’s second distribution, built from native/ in this repository at solprop’s own version and installed through the native extra. It depends on the versioned seam; solprop does not depend on its native build.

driver/scaffold.py is an independent fixed-step IP-RK4 reference used by tests and its tutorial; it is not the production stepper. It propagates a GnlseProblem directly, without a contract, a backend or an output handler, which is what makes it a cross-check of the reference backend. processing/windows.py is a candidate module, drawn dashed in figure 1, to be written only for an agreed use case. _display.py is private: the formatter behind every class’s printed block, imported by any layer and importing none. _numba_cache.py is private in the same way: numba’s njit with a disk cache that concurrent processes cannot cross-wire, through which every kernel in T, D and B compiles (backends.md §3.2).

The layout does not mirror Luna’s module tree: what NonlinearRHS.jl does spans response normalisation, translation and backend passes; Maths.jl mostly becomes SciPy calls; the cumulative-operator stepper takes the place of RK45.jl; scans use the driver and executor layers above.

5. Data flow of one run

user script / prop_capillary(...)                                  (D)
   │  builds physics objects: Grid, Modes, ZProfiles, Responses, Pulses, Noise,
   │  Boundaries, StepperSettings, Output
   ▼
Simulation(...)                                                    (D)
   │  .translate() → KernelProblem  (validated arrays, schema_version)   (T)
   │  .backend = backends.select(name, problem)                     (B)
   ▼
run loop over save planes and stats periods                        (D)
   │  backend.set_initial(field)                                    host → device once
   │  for z_save in planes: backend.propagate_to(z_save)            shared step loop
   │                        obs = backend.observe()                 small host arrays
   │                        stats = backend.stats(spec)             one fused pass
   │                        output.save(z_save, obs, stats)
   ▼
MemoryOutput / HDF5Output  ──►  processing / plotting               (X)

The inner loop (many adaptive or scheduled steps between two save planes) is the shared Python StepDriver. NumPy evaluates reference passes, Numba launches compiled kernels and library transforms, and the native backend executes each trial through one native call. The outer driver visits save planes. Device arrays never cross the seam; observe() returns host arrays of the field at the plane and the step history; stats() returns a small dictionary of host scalars/vectors.

The outer loop is solprop.driver.run.run_problem, which takes a translated problem, a backend that already holds its initial state, and an output handler. Nothing in it is specific to a Simulation; it is the piece a batch runner or a scan calls. What the driver does not do is as much a part of the design as what it does:

• It never conjugates a spectrum. The backend works in the FFT convention, the output layer converts on the way to disk, and those are the only two places a sign changes.

• It never touches the state between planes. Everything between two save planes happens inside the backend’s march.

• It never computes physics. The mapping from gamma to a susceptibility is solprop.nonlinear.gnlse, the statistics are solprop.stats, the storage is solprop.output; the driver is the wiring.

A general modal guide adds a chain of host-side records in front of the KernelProblem (vector fields and projection):

ModalProvider  ->  DenseModalTable  ->  ModalTransformTables  ->  KernelProblem
 E, H, beta at      tracked knots in     S and Q on the run's      the three
 nodes, one         frequency, gated     coarse grid, the dense    general-modal
 (omega, z)         on its sampled       E/H and the node rule     leaves and the
                    Gram                                           diagonal operator

6. The physics layer

6.1 Grids (grid.py, fourier.py)

FieldGrid and EnvelopeGrid implement one protocol: time axis t (n_t points, spacing δt), coarse frequency axis omega (n_omega: n_t/2+1 for the field grid, n_t for the envelope grid in FFT order), the oversampled fine axes t_fine/omega_fine (factor ffac: 6 for a cubic field response with third-harmonic content, 4 for |E|²E, 1 or 2 for envelopes), the in-band mask (physical frequencies inside the wavelength limits), the spectral and temporal collar profiles, and the reference frequency ω₀. Sizes are powers of two or 2-3-5-smooth on request. The invariants between coarse and fine grids (Luna Grid.jl:107–108, 244–250) are asserted in the constructors and tested. fourier.py owns the sign convention (conventions.md), the r2c/c2c helpers, the analytic-signal and oversampling operations used outside the hot loop, and the crop/pad scaling factors that the translation layer hands to kernels.

6.2 Modes (modes/) and profiles (profiles.py)

The Mode protocol: neff(omega, z) (complex), dneff_domega(omega, z) (analytic where the model allows; a validated finite-difference fallback otherwise, ADR-006), field(points, z) (transverse field components at points), dimlimits(z), aeff(z), and the traits z_constant, scale_invariant, azimuthal_order. Derived quantities (β, α, β₁, β₂, ZDW, loss length) are protocol-level functions. MarcatiliMode (hollow capillary, full and reduced models, loss as a Boolean, scale factor or callable of ω) shares gas-fill machinery with the tabulated and rectangular hollow modes; the core index is a ZProfile-aware object so gradients and tapers are ordinary cases. n_eff − 1 is computed directly, never as (1 + x) − 1.

ZProfile (profiles.py) represents a scalar function of z — gas density, core radius, pre-ionisation fraction — with its breakpoints (stage boundaries are hard kinks), analytic derivative where available, and evaluation on arrays. It also declares its jumps (a subset of the breakpoints) and encloses its variation over any interval as a PathVariation: pieces with one-sided limits, value bounds and a variation bound, exact on pieces certified monotone. An arbitrary callable must declare its turning points to have one; sampling is not an enclosure. The adiabatic certificate (adiabatic.py, SEP 34 §4.5) builds on it. Pressure gradients use the steady-flow profile p(z) = sqrt(p₀² + (z/L)(p₁² − p₀²)); fills must not start or end at exactly zero pressure (a vacuum-terminated fill caps every scheme at order ~3, stepper.md §10; the interface substitutes a small floor and logs it). A steady-flow fill of a hollow fibre (gasflow.py) follows the bore and chokes at its vacuum end (ADR-032).

modes/quadrature.py builds the fixed transverse rules (Gauss–Legendre radial with Kronrod extension for the embedded estimate; periodic trapezoid in θ, exact for trigonometric polynomials of degree ≤ n_θ − 1, with the 4·h_max + 1 warning; Cartesian Gauss for rectangular modes; a one-dimensional Gauss–Kronrod rule in x with unit y weight for slab modes; for one Marcatili HE₁₁ mode, Gauss and Gauss–Radau rules in the field scale \(\eta = J_0(u_{01}r/a)\) that carry radial area weights) and the synthesis/projection matrices with weights folded in (Luna TransModalFixed §3.3–3.4; ADR-010). radial_rule selects the family; a one-node rule is an ordinary modal transform, and only mode_averaged_transform makes Transform.kind == "mode_averaged".

6.3 Materials (physdata/)

solprop keeps no material database of its own (SEP 3). Every built-in material is optmat’s, adapted at the boundary: a gas’s polarizability (OptmatSusceptibility), equation of state and transport (OptmatGas) and electronic Kerr model at the carrier (OptmatKerr), its first ionisation energy, and the solids’ complex indices (OptmatGlass) and fused silica’s nonlinear bundle (SilicaBundle), each with its exact derivative, its qualified intervals, its status and an identity that the run’s material manifest records. A MaterialSelection names optmat’s models, and the solid factories’ model= names a solid’s. What physdata/ keeps is solprop’s physics around the data — linear density scaling, the mixture rule (mixtures.py), the χ⁽³⁾/n₂ conversions — and the types a caller’s own data take: the gas formula types (sellmeier.py), a custom Species with TabulatedSusceptibility or IdealGas, SellmeierGlass, TabulatedGlass and CompositeGlass, and the molecular Raman constants of raman_types.py. Air’s molar mass is read from CoolProp through pyfluids, the one direct use left (see dependencies.md).

6.4 Pulses, noise, optics

pulses.py builds initial fields on either grid (Gaussian, sech, custom shape, measured data, CW), scaled to energy or peak power, with Taylor spectral phase; noise.py draws shot-noise realisations for the input model and the modified model (one photon per mode per frequency bin, amplitudes per conventions.md, masked to the in-band region, seeded with np.random.default_rng); optics.py applies linear optical elements (Taylor phase, bulk material, chirped mirrors, grating and prism pairs, and chains of them) as spectral transfer functions to a field outside the propagation.

6.5 Linear operator assembly (linear.py) and boundaries (boundaries.py)

linear.py evaluates D_m(ω, z) from the modes and the frame (β₀(z), β₁(z) from analytic derivatives at ω₀) and folds in the spectral absorber rate; it knows nothing about Chebyshev tables. For Marcatili modes with ZProfile core indices it exposes the three-scalar structure (a(z), ρ(z), β₁(z)) that lets the translation layer sample D cheaply at many z (the idea behind Luna’s MarcatiliLinop). boundaries.py defines the collar geometry (spectral edges beyond the wavelength limits; temporal collar as a fraction of the window, guaranteed non-degenerate) and converts profiles to rates with a reference length ℓ: α = −2 log(W)/ℓ for the spectral collar and the bounded 12.6·(1 − W)/ℓ for the temporal one (ADR-007).

6.6 Responses (nonlinear/)

There is one definition per response, not one per layout or per device: a frozen dataclass that carries its coefficients, declares four traits (ResponseTraits, ADR-011) and provides a NumPy evaluate that is the reference implementation — the oracle every backend kernel is tested against, written for clarity rather than speed.

trait	values	Kerr	dispersive Kerr	plasma	envelope plasma	Raman	silica tensor Raman
domain	pointwise, memory (sequential in t), convolution, spectral	pointwise	convolution	memory	memory	convolution	convolution
grids	field, envelope, both	both	both	field	envelope	both	field
polarisation	scalar, vector_isotropic, tensor	vector_isotropic	vector_isotropic	vector_isotropic	vector_isotropic	tensor, or scalar for a scalar-only declaration	tensor
z_profiles	names of the scalar paths it consumes	density	density	density, preionisation	density, preionisation	density	none

A response bound to a mixture component names its paths with the label, density:N2 (labelled mixture responses). The translation layer reads the traits and compiles the declared set into one fused time-domain traversal per column. Batched variants do not exist as types: batching is the kernel layout. Adding a response means adding a definition, a scalar body and a test. Envelope plasma uses the same traits but exposes sources rather than one evaluate polarisation: its temporal population scan produces a Drude drive and a directional loss current, which keep their own frequency factors through projection and the temporal transforms (envelope plasma sources).

ionisation.py holds the rate protocol (IonisationRate, EvaluableRate), its optical-cycle conventions, the energy binding of the plasma loss, the rate ceiling MAX_RATE (clipping the rate so a diverged step gives no NaN) and the Keldysh and barrier-suppression diagnostics, and no rate formula; optmat_rates.py adapts optmat’s rate models to the protocol (OptmatRate); raman.py the response models (damped oscillator, rotational/vibrational models of a supplied species, silica’s sampled bundle shape and tensor response); normalisation.py the UPPE and envelope prefactors and the effective-area normalisation for mode-averaged runs; gnlse.py the mapping from (γ or n₂/A_eff, β coefficients, f_R, shock) to a mode-averaged envelope problem. optmat_raman.py adapts optmat’s gas Raman responses onto the same terms: one isotropic and one rank-two row per damping group, its decay law optmat’s own, and a mixture’s one response driven by the summed density of its labels.

nonlinear/envelope_plasma.py defines an independent Cartesian-cycle oracle for envelope plasma. Translation compiles instantaneous-rate tables, validated angular rules and separate Drude/loss frequency factors. Backends recompute each channel’s cycle moments and temporal population for every RHS input, then project and transform both current streams before adding their equivalent polarisations to the ordinary source. The envelope plasma manual derives the directional loss current, source signs and averaging limits; the kernel contract defines the numeric data crossing into a backend.

RamanResponse stores immutable RamanTensorTerm declarations. Translation preserves each term’s density and decay law, then supplies numeric trace and traceless weights. The local Cartesian response retains delayed off-diagonal coherence. Scalar and isotropic-only vector responses use one real history; rank-two vector responses use three histories and two effective kernel spectra. The Raman manual derives the contraction and its scalar and circular limits. Circular launch still uses the existing real modal basis.

7. The translation layer (translate/)

The translation layer produces the KernelProblem (kernel-contract.md): plain arrays with fixed dtypes, shapes and axis names, validated (finite, contiguous, sizes consistent), immutable, versioned. What it computes, once per problem:

product	from	representation
Φ table	linear.py + ZProfiles	piecewise Chebyshev in z, coefficients (n_intervals, n_deg, n_batch_or_1, n_modes, n_pol, n_omega) as in kernel-contract.md §3, breakpoints, hard kinks; constant-D fast path (degree 1: one interval, or one per stretch between the kinks of a response path)
unitary phase for noise	Im Φ	same table, imaginary part
quadrature	modes/quadrature.py	nodes, weights, coarse weights, synthesis and projection matrices; taper scaling rule
response spec	nonlinear/responses.py	compiled ordering of the fused pass, coefficients, ionisation-rate spline tables (log-spaced, converged, clip limits), Raman kernels on the doubled grid in the frequency domain, normalisation prefactors
masks and rates	grid.py, boundaries.py	in-band mask, spectral rate (already in D), temporal rate, hard band mask, crop/pad scaling; analytic fine-grid collar derivative and steady-current certificate for field plasma
initial state and noise realisation	pulses.py, noise.py	canonical layout
stepper settings	stepper/	tableau, tolerances, controller parameters, schedule, save planes, caps
statistics requests	stats.py	which statistics, which windows
spatial plane measure	spatial.py, adiabatic.py, translate/flux.py	the FluxPlan: separable coefficients whose root gives each column’s axial admittance, the material square’s z table for a varying radial medium, and its adiabatic admission (SEP 34)
batch parameters	driver/batch.py	per-member scalars or per-member table indices

Validation is not optional: a backend calls problem.validate() and refuses anything inconsistent. Tables are cached on the Simulation and rebuilt only when a physics object changes.

8. The solver, in one page

Details and evidence are in stepper.md. The anatomy of one stage i of the interaction-picture DOPRI5 step from anchor z_n, for a multimode field-grid problem, as the backend executes it (each numbered line is one pass over its data; the FFT counts are the cost metric):

1. stage combine + to lab frame + noise + pad: u = exp(ΔΦ_i)·(y + h Σ_j a_ij k_j) + N(z_n + c_i h), written zero-padded into the fine-grid spectral buffer — one fused pass over (n_modes, n_pol, n_omega).

2. batched inverse FFT (fine grid), n_modes·n_pol columns.

3. synthesis GEMM onto a chunk of quadrature nodes: (nodes·n_pol) × n_modes by n_modes × n_to.

4. fused response pass per column: |E|² → ionisation rate (table) → two plasma running integrals (exposure and free current) → physical current J. Kerr and Raman accumulate separate non-plasma polarisation. Raman adds one or three real histories per physical node, transformed in a doubled-grid batch. The column writes three complete streams: collared non-plasma polarisation gP, gJ and g′J.

5. projection GEMM for each stream back onto modes, accumulated over chunks.

6. batched forward FFT (fine grid) for each stream, followed by DFT[gP] + DFT[gJ]/(iω) − DFT[g′J]/(iω)² on active bins. Inactive bins are exactly zero and no division is evaluated there.

7. crop + prefactor + spectral window + band mask + back to the interaction frame + accumulate: k_i = exp(−ΔΦ_i)·p(ω)·W(ω)·crop(source), acc_sol += b_i k_i, acc_err += e_i k_i — one fused pass; on the last stage the weak error norm is reduced in the same pass (a two-stage reduction, its order free within its rounding bound).

Without field plasma there is one projected source stream. The held-current source adds two temporal transforms to that modal route. Radial propagation instead transforms and combines the three streams per node, crops, then projects once in frequency. Cartesian and line propagation combine before the transverse FFT; planar propagation projects in x before the temporal transforms and combines before the y transform. These orders and their buffer budgets are part of SEP 27.

These are passes over data, not launches.

Per step: six new nonlinear evaluations with FSAL (seven without it), with the modal batched transforms above plus the doubled-grid pairs for Raman, 6 exp(ΔΦ_i) evaluations (or 5 with FSAL and memoisation), one controller decision (lagged, with rejection allowed by default; reject_free is explicit), one accept pass (y ← exp(ΔΦ)·(y + h·acc_sol) — the stages live in the interaction variable anchored at the step start, so the whole sum is carried by one exponential, and the embedded difference takes the same factor — rebase the FSAL stage, apply the temporal absorber per ADR-025 when active), one hard-mask multiply. Steps land on save planes and hard kinks (ADR-009). Where the linear flow couples channels (§16), each exp(±ΔΦ_i) multiply becomes the action of a small per-bin matrix on the mode axis, and a frozen-residual section adds its explicit residual to the stage source (stepper.md §2).

Memory budget (field-sized unit F = n_batch·n_modes·n_pol·n_omega·16 B): state y, stages k₁…k₆, two accumulators, up to five propagator factors (tabulated) or none (recomputed) — 9–14 F by design, more when both propagator signs are tabulated (recompute halves that share); fine-grid spectral and modal time-domain buffers — 2–3 F·ffac; node-chunk buffers (chunk, n_pol, n_to) for E and P, one or three doubled real Raman histories per node and their spectra/convolved values — bounded by the chunk size, not by the node count. Field plasma also charges its current and product streams, projected streams, transforms and provider plans. Planar node admission reserves at least one complete x row. Effective kernels and provider workspace are counted separately in the backend ledger. Examples:

case	F	solver registers	node buffers (chunk 16)	Luna’s node buffers (all nodes)
RDW VUV: n_t 4096, 4 modes	0.13 MB	~2 MB	~3 MB	~8 MB
CtC: n_t 32768, 6 modes, Raman	1.6 MB	~20 MB	~45 MB	~210 MB

The low-storage march (SEP 24, ADR-050) is RK4 in the midpoint interaction picture, replay-only and accepted in place: each slope is folded into a running sum as it is produced, so a step holds the state, an accumulator and a stage register beside the evaluation’s scratch, plus the half-step flow exp(hD/2) when tabulated — four or five F where the Runge–Kutta march’s ledgers hold 13 to 16. It makes four source evaluations a step and keeps no estimate.

A backend exposes a memory_ledger() naming every buffer and its size; the ledger is tested. The table shows the design’s scale, not a per-configuration allocation promise: read the constructed backend’s ledger for its actual buffers.

9. Backends, in one page

Details in backends.md. The seam is at the level of operations on the propagated state, never at the array-namespace level (solprop.backends.protocol.Backend; signatures abbreviated):

class Backend(Protocol):
    name: str
    max_batch: int

    def set_initial(self, field, noise, *, z: float = 0.0) -> None: ...
    def set_progress(self, reporter) -> None: ...  # an observer, never an input
    def propagate_to(self, z: float) -> StepSummary: ...  # shared step loop
    def observe(self, *, field: bool = True) -> Observation: ...  # host arrays only
    def stats(self, spec) -> dict[str, np.ndarray]: ...  # one fused pass
    def metadata(self) -> dict[str, object]: ...  # what the execution was
    def restore_call_ledgers(self, logical, atomic=None) -> None: ...  # after a resume
    def available_device_memory(self) -> int | None: ...
    def memory_ledger(self) -> dict[str, int]: ...
    def close(self) -> None: ...

• numpy_ref: the readable oracle. Plain NumPy + scipy.fft, no fusion, no in-place tricks, one function per pass. Every other backend is validated against it.

• numba_cpu: fused stage kernels (@njit(cache=True, nogil=True) through solprop._numba_cache.njit), pyfftw plans executed in place on preallocated aligned buffers (FFTW_MEASURE + wisdom on a quiet machine), node-parallel kernels and controlled FFTW/BLAS threads, with two-stage reductions and the shared Python step loop. Its floating-point licence is SEP 25’s: contraction by default, strict opt-outs for the protected helpers (stable increments, small-argument functions, defects, compensated sums), reassociation only on the loops of a kernel that asks for it, and no other fast-math flag (backends.md §3.2), applied by the njit wrapper. Fixed fibre batches have a compiled accept update over members; other numerical work remains member-owned. An all-compiled march is a deferred experiment, subject to a measured dispatch-cost case.

• native (the solprop-native distribution, native/): C++20/Kokkos kernels and persistent Kokkos-FFT plans behind nanobind. It runs mode-averaged and multimode scalar/two-component problems, Kerr, plasma and declared tensor Raman, weighted control, embedded quadrature and backend-owned statistics with one native call per trial, and on its CPU build also the general modal route and the coupled flow. Raman uses doubled-grid plans, actual-width tails and memoisation by decay rates and mixture amplitude densities. KokkosKernels supplies the portable GEMM seam; measured local builds may dispatch eligible whole transforms to Accelerate, OpenBLAS or cuBLAS. The CUDA build uses time-parallel response teams, deterministic plasma scans and device reductions and is qualified on data-centre NVIDIA GPUs; the HIP build is qualified on AMD Instinct GPUs, with the HIP sanitizer run still outstanding. Fixed-schedule fibre batching is implemented through select_batch. The private state is exposed to the shared driver only through its protected state type and the _observed_field host-copy hook.

• JAX: evaluated in a measured comparison (ADR-002); not a planned backend unless the data reverse the decision.

backends.md §9 records which route each backend admits and where each is qualified, and §11 its current limits.

10. Output and statistics

MemoryOutput (everything in RAM: fast, small problems, batch scans) and HDF5Output (streaming to disk, chunked per saved plane, optional compression, cache/resume) implement one save protocol and one metadata schema (output/schema.py): grid, modes, simulation type, the full configuration that built the Simulation (serialised physics objects), package and dependency versions, git provenance of the running code when available, the noise seed and realisation, the stepper settings and the step schedule actually taken (z, h, error norm, rejections, wall time), the statistics, and the rows a periodic subset of them recorded between planes. Fields are stored in the canonical layout (n_z, n_modes, n_pol, n_omega) (batch: leading n_batch), in the optics sign convention with the attribute spectral_convention = "optics" (ADR-017: the output layer conjugates on save and on resume); processing accessors squeeze trivial axes and accept convention="fft". A TDSE run may also write its raw provider evaluations to a separate archive (output/tdse_archive.py, schema solprop.tdse-archive.v1), which the optical record references and never needs (the archive page).

output/snapshots.py coordinates optional immutable disk snapshots at the existing accepted-step commit boundary. It persists retained march state and execution facts, suspends the optical and optional archive handles while their advisory writer locks remain held, copies and verifies the pair, publishes a manifest last, then reopens the same writers. Backend and atomic resources are untouched; no additional observe, landing plane or response evaluation occurs. The store bounds completed and incomplete generations and leaves transport/retirement to an external controller. The manifest has its own version, outside the scientific configuration and kernel schema. See the manual for admission, restoration, ownership and interruption rules.

Statistics (energy per mode and window, peak power and intensity, FWHM, electron density and peak ionisation rate, density and pressure, core radius, ZDW, the quadrature error estimate, ω₀) belong to saved planes (ADR-012). They are computed by the backend in one fused pass with two-stage reductions and returned as small host arrays; an optional periodic mode computes a subset every N accepted steps for step-level diagnostics.

The periodic mode is the shared step loop’s, not a backend’s: after every N-th accepted step the march calls the backend’s own stats with the subset and appends one row, which observe drains into Observation.step_stats at each plane, so there is one implementation of every statistic (backends.md §1, “The statistics channel between planes”). The catalogue’s between_planes trait says which statistics may be in a subset, and StatsSpec refuses the rest by name — a root search or a mode-profile rebuild per accepted step is what cost Luna a fifth of its step. The native backend reduces and packs the statistics without observing the complete field. The rows are stored beside the step history and in the same ragged shape.

11. Scans and batch mode

Scans run above the backend seam. Scan gives each point a stable global index, parameters and seed. run_scan selects an executor; Local is the API default, Pool uses spawned processes for independent CPU runs, and Queue/Slurm provide durable distributed execution. A controller owns the collected HDF5 scan file.

Fixed-schedule fibre batches use BatchBackend and select_batch. A batch adopts independently translated and loaded member backends, each with its own tables, state and noise; it advances them under one compatible explicit schedule. The accepted update is fused on compiled backends, while response and transform work remains member-owned. observe_members preserves each member’s history and statistics. Admission uses memory ledgers and a free-device-memory query with headroom.

Every backend has exact member-versus-singleton identity checks under the same schedule. That does not imply a throughput improvement on every workload: whether a batch beats independent runs depends on the problem size and the device. Adaptive lockstep is not implemented. A group-level failure fails every member of that group; recovering the surviving members is an open question (roadmap).

Independent simulations use spawned processes, normally one scientific thread per worker, rather than threads, because thread pools, FFTW plans and the OpenMP runtime are process-global (ADR-013). Each member may use its own adaptive schedule outside a fixed batch.

12. Parallelism and the GPU model

Where parallelism lives at each scale. Construction-time pool choices use measured size rules with explicit overrides (backends.md §3.4); independent process counts belong to the scan executor:

problem	inside one step	across steps	across problems
small single	serial kernels and small transform pools by measured size rules	sequential	spawned processes
large single	prange over node chunks; FFTW/BLAS threads on batched plans	sequential	processes
compatible fibre problems on GPU	member-owned stage work; fused accepted updates	fixed lockstep	backend batch
large on GPU	one thread per element	sequential	one problem per device

The plasma cumulative integral is sequential in value but not necessarily in execution. On the CPU it is the inner loop of the fused pass. The CUDA path uses one team per node: rate evaluation is time-parallel and a fixed deterministic scan tree forms exposure, fraction and currents before the response finish. Response teams are 256 below 8,192 time samples and 512 at and above it; deterministic reduction blocks are 512. Fixed batches are implemented; their measured crossover depends on workload (backends.md). Public selection prefers the installed native backend, otherwise Numba CPU (ADR-029); the installed native build determines its execution space.

13. Memory discipline

• Every buffer is named, preallocated at backend construction and listed in memory_ledger().

• The Runge–Kutta march holds every stage row its tableau and dense output read, and running accumulators for the solution and the error estimate; the low-storage march folds each slope into its registers as it is produced (ADR-009 as amended by ADR-050).

• Space-domain work runs over node chunks; the chunk size is chosen from a memory budget (default: the whole node set when it fits a per-problem budget, else the largest chunk that does).

• HDF5Output never holds more than the current plane, in datasets preallocated for the planned planes with an n_saved count; MemoryOutput is allocated for the planned planes, grows by doubling past them, and is charged by the memory admission.

• A run keeps its fields on every plane, the named planes or none (FieldRetention, SEP 22): a plane without its field is not copied, converted or stored, the stores allocate field storage for the kept planes only, and the admission charges those alone (its output group). The backend is asked for the field only on the kept planes: StepDriver.observe(field=False) transfers nothing from a device, scans nothing for finiteness, skips the collar diagnostic (counted in the execution record) and forms no noise-phase volume, so a field-free plane costs its statistics alone. Under the low-storage march a field-free plane holds one state volume so that a run stopped in flight can resume (held_plane, in the ledger and as its own admission group; on the native backend in the memory space).

• The Cartesian and planar collar diagnostic reconstructs the physical field by block of frequency, at most a node chunk’s worth of samples at a time, never the whole volume.

• A SimulationRunner keeps one backend’s buffers and plans across many inputs, and is admitted once; StepDriver.set_initial restores everything a march reads, so a reused backend repeats a fresh one’s run. An input array of more than 2²⁰ elements is recorded by its SHA-256 digest, not inline, so the run record stays small.

• Φ tables are shared where possible; the constant-D fast path stores one array.

• Fine-grid buffers use the smallest ffac the declared responses need.

14. Logging and progress

logging.getLogger(__name__) per module; the package installs a NullHandler; solprop.configure_logging(level, file) for scripts. Debug logs are plentiful in P, T and D (translation summaries, table sizes, schedule statistics, memory ledger) and guarded with logger.isEnabledFor(logging.DEBUG) when formatting costs anything. Nothing logs inside compiled kernels or per-stage code. Progress is reported by the driver through one ProgressReporter protocol (solprop/log.py): started once, step from inside the step loop on every attempted step, plane per save plane and finished at the end. The default LogReporter turns those into three kinds of INFO line – a periodic one on a wall-clock period (progress_period, 10 s, 0 to disable) carrying position, counts, current step, elapsed and an ETA; one per plane with the step count and mean h; and a closing one with the achieved-against-requested tolerance ratio and how far off the run’s own ETAs were. The period is wall time, not a step count, so the log volume of a run does not depend on how fast its steps are. Where a person is watching – standard error a terminal, or a Jupyter kernel – the default is instead a BarReporter: a LogReporter whose periodic event redraws a tqdm bar carrying the same content and the same measured ETA, with the plane line demoted to DEBUG; redirected or captured output always gets the lines. A reporter that implements interrupted (InterruptibleReporter) is told when the loop raises, before the error propagates, so a bar can give the terminal back.

15. Extension points

• Modes: implement the Mode protocol and declare the profile geometry, frequency support, derivatives and z dependence. Capillary, antiresonant, tabulated, step-index and rectangular implementations provide current examples. A full-vector E/H source implements the separate ModalProvider protocol instead, declaring its material regions, its domains and a serializable identity; the bound step-index and radial graded-index providers are the current examples. A cross-section is polar, cartesian or slab; a new one needs a quadrature_rule branch, taper_scale’s scaling law and on_axis_rule’s node. A mode uniform in one transverse direction declares transverse_measure = "per_unit_width", which the run records and processing and plotting follow for units; the kernel never sees it. A grid uniform in y (LineGrid) is the same measure, recorded from the free geometry’s single y row. A spatial kind whose columns are not its nodes ("slab") says so through SpatialGeometry.n_columns; its transform carries the x half of the map in synthesis and projection and the y half as FFT axes.

• Materials and responses: species and mixtures carry optical and nonlinear data; the translation layer compiles supported Kerr, scalar/tensor Raman and field/envelope plasma terms. A new tensor or spectral response needs a written design with its physical derivation, grid/backend scope and contract changes, and independent tests.

• Spatial propagation: radial and Cartesian transforms already share the seam. Additional geometry or varying Cartesian media are separate extensions to qualify against the existing spatial contracts.

• Execution: scans separate point identity and output collection from the executor. New execution policies must preserve that identity and be measured on complete jobs.

The roadmap records selected work and optional extensions.

16. Coupled linear flows, and non-goals

The linear operator need not be diagonal in the modes. Where channels couple linearly — a capillary’s coherent wall, the taper coupling of a moving multimode guide (the connection of a basis that moves along z), a declared unitary basis of a split block on the general modal route — the problem carries the optional linear.coupled_flow leaf (translate/coupled.py) beside the diagonal Φ table. Its generator does not commute with itself at different z, so the flow is built per section of the operator table’s breakpoints, in one of three ways:

• "constant_schur": a constant block per section, whose exponential is formed from its Schur decomposition and is exact over any step inside the section;

• "frozen_residual": the block frozen per section inside the exponential, and the rest of the generator — its z variation and the connection — as an explicit residual that the Runge–Kutta source carries, so the embedded pair sees its error;

• "cf4": a fourth-order commutator-free exponential of the whole generator at Gauss nodes of each stage interval, whose flow error the embedded pair does not see (stepper.md §2); admitted on numpy_ref only and never selected by default.

The coupled flow runs on numpy_ref, numba_cpu and the native CPU build, which form the per-bin flow actions in their own kernels and agree with the reference to a derived rounding bound; a certificate records each flow’s conditioning. A multimode capillary taper carries the taper coupling by default; taper_coupling=False is the explicit opt-out to the adiabatic route (the taper coupling). The general modal route and its moving guides are described in general modal propagation.

Non-goals. Extended precision; adaptive spatial quadrature; Julia interoperability. Queue and Slurm executors are implemented; more general distributed task graphs are outside the scan model.