# SEP 23 — Far-field aperture and angular-region spectra for free-space runs | | | |---|---| | **SEP** | 23 | | **Title** | Far-field aperture and angular-region spectra for free-space runs | | **Author** | John Travers | | **Drafting agent** | Opus 5.5 (2026-10-16) | | **Status** | Final | | **Review status** | Reviewed | | **Reviewers** | John Travers (human); Opus 5.5 (AI) | | **Type** | Feature | | **Created** | 2026-10-16 | | **Resolution** | Accepted 2026-10-18 by John Travers after his final review, on the revisions of 2026-10-18 recorded in Part III; Part I is frozen from here | | **Requires** | [SEP 22](sep-0022-memory-lean-spatial-execution.md) phase 1, for the `omit_default` record rule; if this SEP lands first, its phase 1 introduces that rule | | **Supersedes** | — | | **Superseded-by** | — | | **Implementation** | phase 1: Opus 5.5 (AI) for John Travers, branch `feature/sep-0023-phase-1`; phase 2: Opus 5.5 (AI) for John Travers, branch `feature/sep-0023-phase-2`; phase 3: Opus 5.5 (AI) for John Travers, branches `feature/sep-0023-phase-3` in solprop and solprop-kokkos; phase 4: Opus 5.5 (AI) for John Travers, branches `feature/sep-0023-phase-4` in solprop and solprop-kokkos | | **Manual pages** | [Far-field apertures and angular regions](../../../guide/results/statistics.md#far-field-apertures-and-angular-regions); [far-field and Fourier-plane statistics](../../../guide/free-space/model.md#far-field-and-fourier-plane-statistics); [catalogue](../../../guide/results/statistics-catalogue.md); [definitions](../../../guide/results/statistics-definitions.md#far-field-spectra); [tutorial](../../../auto_tutorials/free_space/plot_07_far_field_statistics.rst) | ## Part I — Proposal ### 1. Abstract Add six statistics to Cartesian and line free-space runs that record what a spectrally resolved detector behind the sample would measure: the spectrum collected through each of a set of circular apertures in a plane at distance $z_m$ (a spectrometer entrance, a pinhole, the holes of a mask), the on-axis spectrum after the same aperture acts in a Fourier plane, the spectrum in boxes of transverse wavevector, the whole-beam spectrum, its axial-flux counterpart, and the power the spectral, transverse and evanescent absorbers remove per unit length. A plane-wave component of transverse wavevector $\boldsymbol\kappa$ reaches the mask plane at $\mathbf x\simeq z_mc\,\boldsymbol\kappa/\omega$, so an aperture is a *chromatic* weight on the stored transverse spectrum and needs no inverse transform. The tables travel on `StatsSpec` as leaf fields, like `energy_windows`, so the kernel schema does not change; one pass reads the state once per plane. Success is judged by closed-form profile values and Gaussian far-field fractions, identities with `energy` and `on_axis_spectrum`, Parseval and tiling to derived rounding bounds, the absorbed power against the energy decrement, three-backend parity within summation bounds, exact hard-edge membership, and thread-count agreement on `numba_cpu` within a rounding bound fixed before measurement. ### 2. Motivation and scope A spectrally resolved measurement after a thin nonlinear sample rarely sees the whole beam. The light leaves the sample, travels a distance $z_m$ (or passes a lens of focal length $f$), and part of it passes an aperture into a spectrometer. In that plane each plane-wave component has settled at a position fixed by its direction, so the aperture selects *directions*, and because a fixed direction corresponds to a transverse wavevector proportional to frequency, it selects a different part of the transverse spectrum at every frequency. Transient-grating frequency-resolved optical gating (TG-FROG) is the canonical case [1, 2]. Three beams pass through holes at three corners of a square and are focused into the sample, the folded BOXCARS geometry [4]; two of them write a transient grating that diffracts the third into the direction $\mathbf k_1-\mathbf k_2+\mathbf k_3$, which emerges towards the fourth corner. A mask passes that corner to a spectrometer, and the trace is the collected spectrum, delay by delay. The same measured quantity, a spectrum of a parametrised nonlinear process, is the input of the general retrieval algorithms of the FROG and d-scan family [3]. Simulating such a trace means propagating a state of order $768\times768$ transverse samples by 600 frequencies (5.3 GiB) through the sample for each of hundreds of delays, and keeping, at each delay and save plane, one spectrum per hole: a few kilobytes. Other uses are the angularly resolved spectrum of a filament's conical emission (concentric centred apertures, whose differences are annuli), a dark-corner diagnostic (light reaching a hole that no beam should reach: crosstalk, aliasing or scattering from the collar), and the spectrum through a slit on a [line grid](../../../guide/free-space/model.md). Today a user must save every plane's field and reduce it afterwards, which is prohibitive for gigabyte states, and the statistics catalogue has no per-frequency transverse reduction: `on_axis_spectrum` reconstructs one point, and `energy` integrates over frequency. **In scope.** Cartesian ({py:class}`~solprop.grid.FreeGrid`) and line ({py:class}`~solprop.grid.LineGrid`) free space on field and envelope grids, one or two polarisation components; circular apertures with four edge profiles in the paraxial far-field (or Fourier-plane) mapping; boxes of transverse wavevector; the beam and axial-flux spectra; the absorbed power of the channels diagonal in the stored $(\boldsymbol\kappa,\omega)$ representation; `numpy_ref`, `collect` on a saved plane, `numba_cpu` and Kokkos CPU, with GPU qualification later under the standing policy. **Out of scope.** Radial runs (the mask-plane radius is Cartesian) and the hybrid planar kind (its columns are slab modes), both refused by name; Fresnel transmission, detector response and polarisation analysers (Q1); tabulated or non-circular apertures (Q2, Q6); the exact $\tan\theta$ mapping (Q5); near-field apertures, which need an inverse transverse transform; a medium-aware normalisation, which belongs to [SEP 5](../sep-0005-solid-ionisation.md) (Q3). No new dependency, output-schema change or changed default. ### 3. Context The [statistics page](../../../guide/results/statistics.md) and the [free-space model](../../../guide/free-space/model.md) describe today's spatial statistics: `energy`, `on_axis_intensity`, `peak_intensity`, `beam_radius`, `fwhm_r` and `on_axis_spectrum`, the default set of `stats.default_spec("free")`. Their host definitions are in `solprop.spatial` (`spatial_energy`, `on_axis_spectrum`, …), using the geometry's `spectral_weights` ($\Delta x\Delta y/N$ for the raw transverse DFT, [Cartesian free space](../../../guide/free-space/cartesian.md)) and `fourier.spectral_density`. In a run, `backends.spatial.collect_spatial` evaluates them on `numpy_ref` and also on `numba_cpu`, whose spatial statistics are this same host code over a borrowed physical-space scratch; `_axis_statistics` forms the on-axis spectrum as `np.einsum("r,rpw->pw", geometry.on_axis, field)`. `stats.collect(..., space=...)` evaluates the saved-field ones on a plane after the run. In solprop-kokkos, `backend.py` maps the spatial names to the flag values 1 to 16 of the native `spatial_statistics` (`_SPATIAL_STAT_FLAGS`), whose check in `src/general.cpp` refuses any other bit (`flags & ~31`); its `axis_spectrum` pass gives one work item to each frequency and loops serially over every column. The Python layer refuses an unknown statistic by name, and a fresh run collects the entrance plane's statistics before its first step, so the refusal precedes any propagation. A `StatDefinition` records `backend` — whether only a backend can form the value because the saved field does not carry what it needs — and `between_planes`, whether the periodic channel may evaluate it ([kernel contract](../../kernel-contract.md) §8, ADR-012). The contract carries the whole request as one leaf, `stats_spec`, a frozen `StatsSpec` validated against the catalogue; `energy_windows` is a field of that leaf, not a row of the contract. [SEP 9](../sep-0009-propagation-statistics.md) proposes beam measurements on modal routes and leaves the spatial routes' statistics as they are. [SEP 5](../sep-0005-solid-ionisation.md) proposes a medium-aware axial flux for spatial routes, with admittance $Y_z=\operatorname{Re}k_z/(\mu_0\omega)$; this SEP reports that quantity's spectrum as a statistic without changing any normalisation. [SEP 22](sep-0022-memory-lean-spatial-execution.md), drafted alongside, provides statistics-only save planes, which make these spectra the only record a long scan keeps, and a transverse-wavevector absorption rate, `SeparableOperator.kperp_rate`, whose power the second channel of `absorbed_power` reports. [SEP 21](../sep-0021-absorbing-boundaries.md) proposes cumulative `absorbed_energy` statistics for the temporal and spectral collars (in J); `absorbed_power`'s spectral channel is the rate (in J/m) whose integral along $z$ is the spectral one of those. SEP 21 accumulates its energies per accepted step by trapezoidal quadrature, so the two agree within SEP 21's quadrature bound, not to rounding, when both exist (a verification item SEP 21 §7 carries once its statistics land). SEP 21 leaves the spectral collar and the transverse collars out of its own scope, so there is no overlap in what is changed. The roadmap's known issue *The vacuum flux factor in spatial functionals* applies to every spectrum here, and its *Further offline elements* item names an aperture filter that could share the profile definition (Q4). ADR-007 (absorbing boundaries as rates) and ADR-017 (the sign convention) are the decisions this builds on. ### 4. Theory SI units and the conventions of [conventions.md](../../conventions.md) throughout. $E_{c,q,k}$ is the stored internal spectrum of column $c$, polarisation $q$ and frequency bin $k$; $\omega_k$ is the bin's physical angular frequency on either grid kind. #### 4.1 Which way a stored column travels In the internal FFT convention a plane wave is $e^{i(\omega t-\mathbf k\cdot\mathbf r)}$. Column $c=i_x+n_xi_y$ carries the axis values $\mathbf k^{(s)}_c=(k_x[i_x],k_y[i_y])$ of `FreeGrid.k_x` and `k_y`, in FFT order, and the inverse transverse DFT synthesises it with $e^{+i\mathbf k^{(s)}\cdot(\mathbf x-\mathbf x_0)}$. At $\omega>0$ it is therefore the component $e^{i(\omega t+\mathbf k^{(s)}\cdot\mathbf x-k_zz)}$, which travels along $(-\mathbf k^{(s)},k_z)$. The physical transverse wavevector, in the optics sense that a beam travelling towards $+x$ has $\kappa_x>0$, is $$ \boldsymbol\kappa_c=-\mathbf k^{(s)}_c . \tag{1} $$ The [crossing-beam lesson](../../../introduction/lessons/plot_04_crossing_beams.rst) builds its $+x$ beam with the internal phase $e^{-i\kappa x}$, that is at $k_x^{(s)}=-\kappa$. Every existing statistic depends on $|\mathbf k^{(s)}|^2$ alone and never met this sign. Every position and wavevector at the interface below (aperture centres, region edges) is physical, and the sign lives in one place, the mask-plane radius (4). #### 4.2 From the transverse spectrum to the mask plane The light leaves the sample through a planar exit face normal to $z$ into vacuum (air, to the precision of the mapping). Phase matching at the face conserves $\boldsymbol\kappa$, so a component's internal transverse wavevector is also its wavevector in vacuum, where it travels at the angle $$ \sin\theta=\frac{c\,|\boldsymbol\kappa|}{\omega}, \tag{2} $$ whatever the medium of the run; only components with $|\boldsymbol\kappa|<\omega/c$ leave the face. Fresnel transmission, which depends on angle and polarisation, is not modelled, as at every end of a simulation. In the far zone of the exit plane — Fresnel number $a^2/(\lambda z_m)\ll1$ for a source of transverse size $a$ — the field in the direction $\hat{\mathbf s}$ is proportional to the angular spectrum at $\boldsymbol\kappa=(\omega/c)\hat{\mathbf s}_\perp$, by stationary phase [5]. In the plane $z=z_m$ that direction meets $$ \mathbf x=z_m\tan\theta\;\hat{\boldsymbol\kappa} =\frac{z_m\,\boldsymbol\kappa}{\sqrt{\omega^2/c^2-|\boldsymbol\kappa|^2}} . $$ The **paraxial step** replaces the vacuum $k_z$ by $\omega/c$: $$ \mathbf x\simeq s(\omega)\,\boldsymbol\kappa,\qquad s(\omega)=\frac{z_m c}{\omega} \quad[\mathrm{m^2}] . \tag{3} $$ The paraxial position is $z_m\sin\theta$, short of the exact one by the factor $\cos\theta$: a relative error $1-\cos\theta\approx\theta^2/2$, $1.25\times10^{-3}$ at 50 mrad. Within paraxial thin-lens optics the same relation is exact in the back focal plane of a lens of focal length $f$, with $z_m\to f$, whatever the distance from the sample to the lens [6]; the mask may be a far-field plane or a Fourier plane. Under (3) the components that cannot leave the face, $|\boldsymbol\kappa|\ge\omega/c$, land at $|\mathbf x|\ge z_m$. Every weight is therefore defined as zero wherever $s(\omega)|\mathbf k^{(s)}_c|\ge z_m$, whatever its profile — the `tanh` and `supergauss` tails are unbounded and would otherwise reach those columns, which carry energy inside a dense medium — so an aperture confined to $|\mathbf x|w$; for `tanh` up to its exponential tail, which leaves a relative area error of $-3.5\times10^{-3}$ at $w=D/4$, $-1.3\times10^{-5}$ at $w=D/8$ and $10^{-11}$ at $w=D/20$ by numerical quadrature, so quadrature weighting admits `tanh` only up to $w=D/8$, section 5). Evaluating the profile at the corrected diameter $$ D_q=2\sqrt{(D/2)^2-4m_2w^2},\qquad m_2=\int_0^\infty t\,W_1(t)\,dt=\frac{\sigma^2}{4w^2}, \tag{6} $$ with $W_1$ the unit-width outer tail, restores the hard hole's area exactly. The values are $m_2=\pi^2/48\approx0.2056$ for `tanh` and $m_2=1/4-2/\pi^2\approx0.04736$ for `rcos`; `hard` needs no correction. The residual for a non-uniform intensity is $\pi\sigma^2r\,\partial_\rho\bar I$, of the same order as the one-dimensional slope term, with $\bar I$ the azimuthal mean. On a line the edges are straight and $D_q=D$. These values and the exactness of (6) for a uniform intensity were checked by numerical quadrature while drafting (residuals below $2\times10^{-10}$ at $w=r/10$). Quadrature weighting is refused with `supergauss`, which has no symmetric edge, and coincides with amplitude weighting for `hard`, where $W^2=W$. #### 4.6 The on-axis spectrum through an aperture With $\alpha_c$ the geometry's existing `on_axis` weights, $N^{-1}e^{-i\mathbf k^{(s)}_c\cdot\mathbf x_0}$, $$ A_{a,q}(\omega_k)=\sum_c\alpha_c\,W_a\!\bigl(\rho_a(c,k)\bigr)E_{c,q,k},\qquad S^{\mathrm{ax}}_a(\omega_k)=\frac{\varepsilon_0c}{2}\sum_q\operatorname{sd}\bigl(A_{a,q}\bigr), \tag{7} $$ in J s/m². This is the spectral fluence at the image of the transverse origin in a unit-magnification $4f$ relay whose Fourier plane, at focal length $f=z_a$, holds the mask. $W$ enters once under both weightings, because the amplitude is linear in the field; under quadrature weighting the profile is evaluated at $D_q$, by the same edge argument. With $W\equiv1$, (7) is `on_axis_spectrum`. Mathematically $\alpha_c=\pm1/N$ on even and odd grids alike, but the stored vector carries the rounding of the exponential, so it is multiplied as stored, as a complex number. #### 4.7 Regions of transverse wavevector A region is a half-open box $R_j=[\kappa^{\rm lo}_x,\kappa^{\rm hi}_x)\times[\kappa^{\rm lo}_y,\kappa^{\rm hi}_y)$ in rad/m, with infinite edges allowed, and $$ S^R_j(\omega_k)=\frac{\varepsilon_0c}{2}\sum_{\boldsymbol\kappa_c\in R_j}w_c \sum_q\operatorname{sd}(E_{c,q,k}) . \tag{8} $$ Membership is exact: it compares $-\mathbf k^{(s)}_c$, an exact negation, with the edges. Half-open intervals that partition the real line partition the columns, so regions that tile the plane tile the columns: $\kappa=0$ belongs to $[0,\infty)$, and on an even axis the Nyquist column (stored $-\pi/\Delta x$, physical $+\pi/\Delta x$) to the non-negative side. A quadrant is one row, $(-\infty,0,-\infty,0)$. A region is not chromatic: it selects the same $\boldsymbol\kappa$ at every frequency, hence angles that shrink as $1/\omega$. Regions bounded only by zero and infinite edges (half-planes, quadrants) are the exception, being defined by the sign of $\boldsymbol\kappa$, which is the sign of the angle at every frequency. A chromatic cone is a centred aperture. #### 4.8 The beam and axial-flux spectra The beam spectrum is (5) with $W\equiv1$, $$ S^{\rm beam}(\omega_k)=\frac{\varepsilon_0c}{2}\sum_cw_c\sum_q\operatorname{sd}(E_{c,q,k}). \tag{9} $$ The axial-flux spectrum weights it by the dimensionless obliquity factor $$ O(c,k)=\frac{\operatorname{Re}k_z(c,k)}{\mu_0\omega_k\,\varepsilon_0c} =\frac{c\operatorname{Re}k_z(c,k)}{\omega_k},\qquad S^{P}(\omega_k)=\frac{\varepsilon_0c}{2}\sum_cw_c\,O(c,k)\sum_q\operatorname{sd}(E_{c,q,k}), \tag{10} $$ with $k_z$ the principal root of the compensated square $q_z=(k^2-|\mathbf k^{(s)}|^2)+(k^2_{\rm low}-|\mathbf k^{(s)}|^2_{\rm low})$ of the run's `translate.spatial.SeparableOperator`, the root the stages use, and $O=0$ on inactive and non-propagating bins ($\operatorname{Re}q_z\le0$, the stages' test) and at $\omega_k\le0$. Because $\mu_0\varepsilon_0c^2=1$, $\varepsilon_0cO$ is SEP 5's admittance $Y_z$, and $S^P$ is the spectral density of the scalar-equation axial-flux energy that SEP 5 proposes. It is written as the vacuum-normalised beam spectrum times $O$ so that it shares (9)'s arithmetic and normalisation, the ratio $S^P/S^{\rm beam}$ is directly the beam-averaged obliquity, and the paraxial vacuum limit $O\to1$ is visible; multiplying $O$ back by $\varepsilon_0c$ recovers the raw surrogate $\sum\operatorname{Re}(k_z)|E|^2/(2\mu_0\omega)$. In vacuum $O=\cos\theta\le1$, equal to one only at $\boldsymbol\kappa=0$; for a Gaussian beam of waist $w_0$, $S^P/S^{\rm beam}=1-(kw_0)^{-2}-(kw_0)^{-4}-3(kw_0)^{-6}-\dots$ with $k=\omega/c$, an asymptotic series of one-signed terms (the exact form is in section 7, item 8). In a medium, $O\approx\operatorname{Re}n(\omega)$ for a paraxial beam: this is the factor by which the vacuum functionals of (5), (8) and (9) are low (Q3). #### 4.9 Power removed by the absorbers The real part of the separable linear operator damps each stored element at a power rate that is a sum of channels. With the exact energy weights of `fourier.energy_omega`, $\nu_k=(\delta t/n_t)\times\{1,4,\dots,4,1\}$ on a field grid and $\delta t/n_t$ on an envelope grid, $$ P_j=\frac{\varepsilon_0c}{2}\sum_cw_c\sum_q\sum_k\nu_k\,r_j(c,k)\,|E_{c,q,k}|^2, \qquad j\in(\text{spectral},\ \text{transverse},\ \text{evanescent},\ \text{material}), \tag{11} $$ in J/m, where $r_{\rm spectral}=\alpha_{\rm abs}(\omega_k)$ is the operator's `spectral_rate`, applied on every bin as the stages apply it; $r_{\rm transverse}$ is SEP 22's `kperp_rate`, zero without it; and $r_{\rm evanescent}=2\min(|\operatorname{Im}k_z|,\alpha_{\rm ev})$ on active non-propagating bins, with $\alpha_{\rm ev}$ the evanescent field cap (`spatial.EVANESCENT_RATE`, or SEP 22's `evanescent_rate` setting); and $r_{\rm material}=2\operatorname{Im}k_z$ on active propagating bins, the loss (or, with the opposite sign, the gain) of the medium's complex index, exactly zero for a real index. Then $-\sum_jP_j$ is exactly the derivative of the `energy` statistic along $z$ under the linear flow restricted to these channels, which are the whole of the operator's real part. Excluded are the physical-space spatial collar and the temporal collar, which are diagonal in other representations (the collar's share is monitored by the existing collar diagnostic, `backends.spatial.spatial_collar_fraction`), and nonlinear losses. The four channels are fixed from the first release, so the statistic's shape never changes. #### 4.10 Zero frequency, end bins and the line The factor $s(\omega)$ exists only for $\omega>0$: $W=0$ and $O=0$ at and below zero frequency, so a field grid's zero-frequency bin adds nothing to (5), (7) and (10), active or not. At a field grid's zero-frequency and Nyquist bins every spectrum follows `spectral_density`'s documented rule (an overcount, as the analytic signal does not double them); both bins are zero on every grid solprop builds. On a line ($n_y=1$, $k_y=0$) the spectra (5), (8)–(10) are per metre of $y$ and the absorbed powers in J/m², while (7) stays local; an aperture's $y$ centre must be zero, and (4) reduces to $|s_ak^{(s)}_x+X_a|$, a slit of full width $D$ along $x$ with straight edges and $D_q=D$. A region selects nothing on a line unless its $\kappa_y$ interval contains zero. ### 5. Proposed interface Names are a proposal. A frozen, identity-compared dataclass in `solprop.stats`: ```text ApertureSet( centres, # (n_apertures, 2) m, mask-plane (x, y), physical axes diameters, # (n_apertures,) m, finite and positive distances, # (n_apertures,) m, exit face to mask plane, or a lens's f edges=None, # (n_apertures,): m for tanh and rcos, order >= 2 for supergauss; # None for hard profile="hard", # "hard" | "tanh" | "rcos" | "supergauss", shared by the set weighting="amplitude", # "amplitude" (W^2) | "quadrature" (W once), shared labels=None, # n_apertures unique non-empty strings; default "0", "1", ... ) ``` Validation, raising `OutputError` by name: one to 16 apertures; every array finite and of the stated shape; positive diameters and distances; for `rcos` $0pw", on_axis[rows], W[:, None, :] * field[rows])`: with $W\equiv1$ the product is the field exactly (a multiplication by one), so the result agrees with `on_axis_spectrum` within the summation bound, with one chunk or across chunks (a drafting check found NumPy's `einsum` changing at rounding level under re-chunking, so no exact identity is claimed). `stats.collect` calls the same function for the four saved-field names with `space=FreeGrid`, and refuses `poynting_spectrum` and `absorbed_power` inside `_spatial_statistic`, which every name reaches before the `backend` check when `space` is given. Because `numba_cpu` evaluates spatial statistics through `collect_spatial`, phase 1 delivers it too. **numba_cpu.** Three bodies join `define_bodies`: `mask_plane_radius` (4), `aperture_weight` (an integer profile code, with the raised cosine's early-out) and `region_member` (the half-open test on $-\mathbf k^{(s)}$). A two-stage kernel in `numba_cpu/spatial.py` replaces the host path. Stage one runs `prange` over contiguous bands of columns — whole rows on a plane, blocks of $x$ samples on a line — whose count, $\min(n_{\rm units},\lfloor B/(8n_{\rm slots}n_\omega)\rfloor)$ with a module constant $B=128$ MiB, is fixed by the shapes and never by the thread count. Each band writes its row of partials, shape `(n_bands, n_slots, n_omega)`, in the slot order `[beam, poynting, absorbed × 4, regions…, per aperture (collected, Re and Im of the on-axis amplitude per polarisation)]`, present slots only: at most $n_{\rm slots}=6+n_{\rm regions}+n_{\rm apertures}(1+2n_{\rm pol})$, plus $2n_{\rm pol}$ for `on_axis_spectrum`. Columns run in ascending order with frequency innermost; stage two sums the bands serially in ascending order (a fixed, shape-determined order), and a small finishing step applies $\varepsilon_0c/2$, the spectral prefactor and $|A|^2$. The complex `on_axis` is multiplied as complex, and $\operatorname{Re}k_z$ is recomputed per element by the shared body, never stored. A `shared_mask` path, taken when every aperture shares centre and distance, computes the radius once per element: the same strict expression on the same inputs, so its hard-edge decisions are the general path's exactly and its sums agree within a rounding bound. `on_axis_spectrum` becomes the $W\equiv1$ slot of the pass, which makes the identity of section 7 hold to rounding by construction. The partials and the scale table are allocated at construction from the spec, named in `memory_ledger()`, and reused by the periodic channel. **Kokkos CPU.** A native `aperture_statistics(flags, apertures, profile_code, exponent, regions, intensity_factor)` beside `spatial_statistics`, with its own flag mask (values from 32; 1 to 16 are `spatial_statistics`') and an argument check refusing other bits. `Problem` gains `spatial_k_x` and `spatial_k_y`, packed from `geometry.spectral_axes`; `apertures` carries the host-computed constants and scale table; the bindings expose the method as they expose `spatial_statistics`; Python packs and unpacks in a fixed order (beam, Poynting, absorbed, regions, collected, on-axis) with the existing consumed-count check, and $|A|^2$ is formed natively so only `n_omega` values per aperture cross. A team policy over the same bands, each team's threads over frequency and each thread looping its band's columns in ascending order, associates every sum as `numba_cpu` does, so agreement rests on summation order; the radius is compiled without contraction (a flag for that unit, or round-to-nearest intrinsics on a device). Folding the existing `axis_spectrum` into the pass as the $W\equiv1$ slot is recommended: one work item per frequency, each looping over all $n_xn_y$ columns, is 600 work items at the size above and does not scale on a device. The fold also adopts the host's rule at the zero-frequency and Nyquist bins, which the native pass weights by one rather than four today; the bins are zero on solprop's grids, so no value changes. On a device the slot accumulators live in registers, and a large set may run in groups that re-read the state (measured in phase 3). The partials join the native ledger beside `spatial_statistics_partials`. **Cost.** A pass reads the state once, $16n_xn_yn_{\rm pol}n_\omega$ bytes (5.66 GB at $768\times768\times600$). Each aperture adds about 40 flops per element and, for `tanh` or `supergauss`, one transcendental (none for `hard`; a sine only inside the edge for `rcos`): six smooth apertures cost a few hundred flops and $2\times10^9$ transcendentals at that size, transcendental-bound on a CPU and between memory- and compute-bound on a device. A Runge–Kutta step makes several state-sized passes and transforms per stage, so a pass is estimated at a few per cent of a step, paid at save planes and periodic rows only; the estimate is to be measured. The partials at that size, six apertures, one polarisation and no regions (24 slots, one band per row), are $768\times24\times600\times8$ B $=88.5$ MB (95.8 MB with `on_axis_spectrum`), and never exceed $B$. Benchmarks use `bench_stats.py` on a ten-step Cartesian window ([benchmarks.md](../../benchmarks.md)). | Option | For | Against | Choice | |---|---|---|---| | Save every field, reduce afterwards | No new code | Storage and transfer prohibitive for gigabyte states | Status quo | | Host reduction on every backend | One implementation | A device-to-host copy of the state per plane on Kokkos | Phase 1 only | | Fused banded pass per backend | One state read, deterministic, device-resident | Three implementations | **Recommended** | | Physical-space masks | Near-field apertures | An inverse transform per plane; a different measurement | Out of scope | No new dependency is needed. ### 7. Verification and validation Each tolerance below is set by a stated bound and asserted positive. 1. **Profiles.** Closed-form values at $\rho=0$, $D/2$ and large $\rho$; the raised cosine's compact support and monotonicity (`hypothesis`); the edge symmetry to rounding; $m_2$ against SciPy quadrature to $10^{-12}$. 2. **Direction.** A beam launched towards $+x$ as in the crossing-beam lesson is collected by an aperture at $+X$ and not by its mirror image (bounded by the Gaussian tail); an achromatic tilt's image stays at $z_m\sin\theta$ at every frequency while its spot shrinks as $1/\omega$, so a hole wider than the spot at the lowest frequency collects it fully at every frequency; a fixed-$\kappa$ launch moves as $1/\omega$. 3. **$W\equiv1$.** On a grid whose largest $|\boldsymbol\kappa|$ is below $\omega_{\min}/c$, a centred hard aperture enclosing every image gives `aperture_spectrum` equal to `beam_spectrum` and `aperture_on_axis_spectrum` equal to `on_axis_spectrum` on each backend within a rounding bound fixed before measurement (the summation bound); every weight is exactly one. 4. **Gaussian far field.** A waist-plane Gaussian beam gives the collected fraction $1-\exp[-D^2/(2w_{ff}^2)]$, $w_{ff}=2z_m/(kw_0)$, through a centred hard aperture. The closed form is exact for the statistic's own (paraxial) mapping, so the tolerance is the edge-sampling error, bounded by the energy of the columns within one lattice diagonal of the edge, plus the Poisson and window bounds of the [Cartesian page](../../../guide/free-space/cartesian.md). Smooth profiles are compared with the one-dimensional integral of the same profile (SciPy quadrature, independent of the lattice); under refinement the hard error falls at first order and erratically, the smooth error quickly once $w\ge2\Delta x_m$. The paraxial map's own physical error, $O(\theta^2)$, is documented, not tested. 5. **Quadrature unbiasedness.** A spectrum uniform across the edge region collects the hard area to lattice error, plus the `tanh` tail's residual of section 4.5, with $D_q$, misses it by $\pi\sigma^2$ without the correction, and by the first-order term (section 4.5) under $W^2$. 6. **Parseval.** $\delta\omega\sum_kS^{\rm beam}_k$ equals `energy` within $2\gamma_m$ of the total, $m$ the number of summed terms plus their roundings, with $\gamma_m=m\epsilon/(1-m\epsilon)$; a fixture with non-zero end bins checks the documented overcount instead. 7. **Regions.** Four quadrants sum to the beam spectrum within the same bound; membership of columns exactly on an edge, at $\kappa=0$ and at Nyquist is as section 4.7 states. 8. **Axial flux.** On a vacuum fixture $S^P/S^{\rm beam}$ matches the exact paraxial closed form, the mean of $\sqrt{1-u}$ over $u$ exponentially distributed with mean $2/(kw_0)^2$ (by quadrature, or its Dawson-function form), to the lattice and rounding bound; the asymptotic series $1-(kw_0)^{-2}-(kw_0)^{-4}-3(kw_0)^{-6}-\dots$ has one-signed terms, so no single term bounds its remainder (1.3 times the third term at $kw_0=5$) and it is not the test's reference. The ratio lies strictly below one for any transverse content in vacuum; in a gas at $kw_0\gg1$ it matches $\operatorname{Re}n(\omega)$ from `SpatialMedium` within the obliquity term. 9. **Absorbed power.** On linear runs with a spectral absorber, evanescent content and, after SEP 22, a transverse rate, each channel equals $-dU/dz$ of the closed-form linear decay to rounding, and on a complex-index fixture the material channel equals the propagating columns' decay (loss) and growth (gain, the sign reversed); over one fixed step $h$ at three levels, $U(z+h)-U(z)+hP(z)$ converges at first order and the trapezoid form at second, within the order tolerance of [testing.md](../../testing.md). 10. **Cross-backend.** `numba_cpu` and Kokkos CPU against `numpy_ref` within $2\gamma_m\sum|\text{terms}|$ per element, $m=4n_{\rm modes}$ for the real sums and $8n_{\rm modes}$ for the complex amplitude, plus the documented ulp error of the transcendental routines; hard-edge membership exact (a hard-mask contract) on a fixture that places columns within a few ulps of an edge. 11. **Thread counts and paths.** `numba_cpu` at 1, 2 and 4 threads within a rounding bound fixed before measurement; the `shared_mask` path's membership exactly the general one's and its sums within that bound; a native GoogleTest (`tests/native/test_aperture.cpp`) against a straightforward serial host reduction. 12. **Errors and records.** Missing tables, radial and planar geometry, a bad edge for each profile, quadrature with `supergauss`, more than 16 rows, an aperture outside $|x|= 0`; `src/solprop_kokkos/backend.py` sizes `spatial_far_field_slots` with `far_field_partials(problem, axis_folded=AXIS_SPECTRUM_ONLY)`. solprop `backends/spatial.py` (`AXIS_SPECTRUM_ONLY`, the `axis_folded` keyword of `far_field_layout` and `far_field_partials`) and `backends/memory.py` (`far_field_memory` counts Kokkos's slots without `on_axis_intensity`). - *The four-wave-mixing tutorial* of SEP 22 drops `backend="numba_cpu"`: on the default (Kokkos) backend it ran in 7.0 s against 11.9 s on `numba_cpu`, with every printed figure equal to the digits shown. - *Tests*: solprop `tests/test_far_field_statistics.py` (the warning, once, naming both wide apertures; five sets that build without it), `tests/backends/test_kokkos_far_field.py` (`on_axis_intensity` takes no native slot, ledger equal to the estimate to the byte; a beam, region and absorbed request with apertures declared against the reference); solprop-kokkos `tests/native/test_aperture.cpp` (`DeclaredAperturesCostNothingUnlessAsked`). - *Benchmark*: `benchmarks/sep23/phase4_costs.py` and the ledger entry [SEP 23 phase 4](../../benchmarks.md). - *Pages*: [statistics](../../../guide/results/statistics.md), [backends.md](../../backends.md) §3.13, [roadmap](../../roadmap.md), `CHANGELOG.md`; the solprop-kokkos CHANGELOG. - **Deviations from the proposal:** - *The tutorial's file and scope.* §8 names `docs/tutorials/free_space/plot_05_far_field_apertures.py` with two apertures (one centred, one off-axis catching a tilted beam), the closed form and the quadrature correction. `plot_05` was taken by the low-storage tutorial, so it is `plot_07`, and it covers the wider scope set for this phase: a crossing-beam run with four holes under both weightings and two profiles, the quadrants and the operator statistics. The closed form is checked on the beam 1 hole at the entrance, against the exact fraction of each profile (the hard closed form, and the profile's own integral for the raised cosines). - *What the closed form shows.* §8 expects the tutorial to "show the quadrature correction". On a lattice of six points per far-field radius the raised cosines match their own integrals to $3.5\times10^{-4}$ where the hard hole's staircase reaches $4\times10^{-2}$; quadrature weighting then departs from the hard hole by the slope term of §4.5, $2.8\times10^{-2}$ here (where the edge sits on the beam's steepest part, with $w=2\Delta x_m$), and amplitude weighting by its first-order deficit, $9\times10^{-2}$. The tutorial says so rather than claiming a smaller error for quadrature than for the hard hole on this lattice. - **Decisions:** - *Q10, as decided by the maintainer:* the warning class is a `SolpropWarning` subclass, raised once per set (one warning naming every wide aperture), at `2 <= r/w < 4` (the lower end is the existing refusal); `r/w = 4` exactly does not warn. - *The Kokkos slot count* is a keyword of solprop's layout functions rather than a Kokkos-side formula, so the payload and the estimate share one definition. - *The remaining ledger entry* of §11 is phase 4's measurement of the gate and the slot count (the passes' own entries are phases 2 and 3). - No ADR: nothing here outlasts the SEP beyond what Part I records. - **Verification** (the branch trees, on a busy Mac; every solprop command through `uv run --no-sync` with this branch's wheel installed): - `tests/test_far_field_statistics.py` (52), `tests/backends/test_far_field_statistics.py`, `tests/backends/test_numba_far_field_kernels.py`, `tests/backends/test_kokkos_far_field.py` (16), `tests/test_free_statistics.py`, `tests/backends/test_on_axis_statistics.py` together under xdist, `tests/backends/test_memory.py` (98) and `tests/driver/test_memory.py` (47) separately (the two same-named modules collected in one xdist session cannot find a conftest fixture, a collection quirk unrelated to this phase), and `tests/test_seps.py`: pass. Doctests of the four changed modules (32), basedpyright (no errors), the numpydoc lint, ruff, ruff format and codespell of the changed files. Coverage (`coverage-cpu.ini`) over the far-field modules leaves no added line or branch of `stats`, `exceptions`, `backends.spatial` or `backends.memory` uncovered. - The tutorial runs clean under `-W error`, figures inspected; its checks: `collect` against the run $4.5\times10^{-16}$ (relative), the definition against the run $2.8\times10^{-16}$ in collected fraction, the quadrants' tiling $7\times10^{-16}$, the obliquity 1.44992 against the index 1.45, the absorbed energy $1.0832\times10^{-8}$ J against the energy drop $1.0829\times10^{-8}$ J. The dark hole's reading at the exit ($6.4\times10^{-6}$ of beam 1's hole) moved by 3 % when the y spacing was halved and by $10^{-4}$ when the x spacing was refined, so it is the beams' nonlinear halo, not the grid. - solprop-kokkos `scripts/check.sh` against this solprop branch: green (format, lint, strict typing, reuse, clang-format, clang-tidy, the native tests with the eleven far-field cases among them, ASan+UBSan native tests, pytest 332 passed and 5 skipped, the distribution wheel and its smoke test), Homebrew Clang 23.1.3. Its first run failed on clang-format: two lines of `src/far_field.cpp` from `main` and the new native test, formatted with the pinned clang-format. The branch's local-profile wheel (built before that format-only commit) is installed in the solprop worktree. - Benchmarks: the ledger entry; a `beam_spectrum` sample with apertures declared cost 4.43 times a bare one before the gate and 0.87 and 1.02 times after; the band partials for `beam_spectrum` with `on_axis_intensity` fall from 1.125 MiB to 0.375 MiB. - `poe docs` (`-W -n --keep-going`): clean; the new tutorial executed in 1.3 s. - **End-of-SEP lanes**, on the branches merged with `main` after SEP 27 (solprop 61b70a8, solprop-kokkos 713abe3), on Linux (AWS, a wheel built from the branch with GCC 16): `poe check` (16 600 passed), `poe test-slow` (147 passed) and the Kokkos-selected default (16 600 passed) and slow lanes pass but for two cases of `tests/backends/test_band_plasma_parity.py` (the trapezoid rule on the mixture, on `numba_cpu` and Kokkos: the electron density 4e-8 from the reference against the test's 2e-8). The same two fail on `main` alone on Linux and pass on macOS; they arrived with SEP 27's plasma source and are not this SEP's (recorded for that work). `poe cov`, which ran out of memory beside another lane on a small instance, passed on macOS at 100 % with all 16 602 tests. - **Not run:** GPU (the device pass compiles; it is not run). - **Sources checked:** no new relation. The tutorial's closed form (the Gaussian's collected fraction through a centred hole), the coherence length of the self-diffraction signal, $\Delta k\approx-4K_0^2/k$ for $2\mathbf k_1-\mathbf k_2$ in a dispersionless medium, and the obliquity's paraxial value were re-derived; the slope residual of §4.5 accounts for the quadrature figure ($-2.7\times10^{-2}$ estimated, $-2.8\times10^{-2}$ measured). - **Limitations:** as at phase 3 (the device pass unqualified; radial and planar-hybrid geometry refused; the paraxial map only). The warning concerns amplitude `tanh` only; quadrature `tanh` is admitted at $w\le D/8$ where the centre's error is of order $2\times10^{-7}$. - **Open questions and their destinations** (every question of the SEP): - Q1 — per-polarisation outputs: decided in phase 1 (not offered); the alternative is Q9. - Q2 — tabulated radial profiles: the [roadmap](../../roadmap.md), *Far-field detection beyond SEP 23*. - Q3 — the vacuum flux factor in a dense medium: the roadmap, same item, with SEP 5. - Q4 — the offline aperture filter: the roadmap, *Further offline elements*. - Q5 — the exact mapping: decided while drafting; the option is on the roadmap, *Far-field detection beyond SEP 23*. - Q6 — other shapes: the roadmap, same item. - Q7 — a material channel: decided while drafting, built in phase 1. - Q8 — infinite energy windows: resolved at the phase 1 review. - Q9 — polarisation-resolved detection: the roadmap, *Far-field detection beyond SEP 23*. - Q10 — the `tanh` centre: decided by the maintainer, built here. - The device pass (phase 3's handoff): the roadmap, *The far-field pass on GPUs*. - No new question. - **Handoff:** the end-of-SEP lanes (full gate, `poe cov`, `poe test-slow`, the Kokkos-selected lanes, a Linux GCC build of solprop-kokkos with the strict warnings), then `Final`, the move to `docs/dev/seps/archive/` (the **Manual pages** field and the relative links follow the move) and the roadmap's accepted-list entry removed. The solprop and solprop-kokkos branches merge together: the Kokkos adapter imports `AXIS_SPECTRUM_ONLY`. ## Part III — History - 2026-10-16: Drafted. - 2026-10-16: Reviewed by Opus 5.5 (AI); seven should-fix items applied: the fourth `absorbed_power` channel carried into §6–§7, the axial-flux test against the exact paraxial form, the weight set to zero beyond the paraxial disc, the `tanh` width limit under quadrature weighting, the reference reduction placed in the physics layer, the dependence on SEP 22's record rule declared, and a refusal that names a supporting backend. - 2026-10-17: Restated for [SEP 25](sep-0025-rounding-level-reproducibility.md). Bounds replace bit identities between separate computations: `numba_cpu`'s thread-count invariance (§1, §7 item 11), the `shared_mask` path's sums against the general path's (§6, §7 item 11), the $W\equiv1$ identities with `beam_spectrum` and `on_axis_spectrum` (§6, §7 item 3), the stages after the root is factored out (§6), and `on_axis_spectrum` in the fused pass (§9), all within a rounding bound fixed before measurement; a recorded state is recaptured only after an old-against-new check at that bound. Kept exact as structural: hard-edge membership from the strict, uncontracted radius (§6, §7 items 10 and 11), the weight's exact one, and the configuration digest of a run without the new fields (§7 item 12, §9). - 2026-10-18: Revised against SEP 25 Final and SEP 21 as reviewed: the `einsum` remark restated as a rounding-level fact (§6); the agreement between `absorbed_power`'s spectral channel and SEP 21's cumulative `absorbed_energy` is to SEP 21's quadrature bound, not to rounding (§3). The Kokkos facts of §3 and §6 (flag values 1 to 16 and the serial `axis_spectrum` loop) and the host path of §3 were re-checked against the current tree and stand. - 2026-10-18: Accepted (John Travers) after his final review; Part I frozen. - 2026-10-18: In progress. Phase 1 (the reference pass, `collect` and the manual) implemented on `feature/sep-0023-phase-1`, with no kernel or output schema change; its report is in Part II. - 2026-10-19: Phase 2 (the `numba_cpu` pass) implemented on `feature/sep-0023-phase-2`, with no kernel or output schema change; its report is in Part II. - 2026-10-19: Phase 3 (the Kokkos CPU pass) implemented on `feature/sep-0023-phase-3` in solprop and solprop-kokkos, with no kernel or output schema change; its report is in Part II. - 2026-10-20: Phase 4 (the tutorial, the `tanh` warning, the Kokkos follow-ups and the open questions) implemented on `feature/sep-0023-phase-4` in solprop and solprop-kokkos, with no kernel or output schema change; its report is in Part II. - 2026-10-20: Final (John Travers, on the agreed plan). Phases 1–4 merged with green closing lanes; open questions moved to the roadmap in the phase 4 report; moved to the archive.