SyncMoments

Documentation: https://syncmoments.readthedocs.io. Releases: https://github.com/zzhang0123/syncmoments/releases.

Differentiable JAX/Equinox synchrotron kernels, finite statistical contractions, and polarised radiative transfer for A statistical framework for synchrotron emission (Zhang & Chluba).

This research implementation keeps four operations explicit: physical kernel approximation, finite kernel Taylor expansion, statistical closure, and numerical solution. Passing a kernel or regression test does not certify a Galactic foreground model at 21-cm precision. All relevant discrepancies and uncertainties must be propagated to the chosen scientific quantity.

The positive finite-cumulant PDF reconstruction and the independent 21-cm assessment live in the manuscript repository's validation/ directory. They do not import this package. CumulantExpansion here is the historical name of a quadratic kernel average; it is not that PDF reconstruction.

The main paper covers radiation kernels, population averages, moment/cumulant representations and their error control, including internally mixed emission and pure Faraday rotation. External screens are a special case. General ordered transfer, Magnus, absorption/conversion and the 21-cm examples remain exploratory modules beyond this main model.

Statistical interfaces

QuadraticTaylorExpansion is the descriptive alias for CumulantExpansion; existing imports and Equinox trees retain their identity. Its low-level call continues to contract supplied tensors. For a validated entry point, use model.checked_average(mean, covariance, absolute_error=...): it checks finite real inputs, shapes, symmetry and positive semidefiniteness and returns (prediction, envelope). The check uses correlation coordinates so a large independent variance cannot hide a small invalid block. Its numerical tolerance is 32 * P * eps; no covariance projection or eigenvalue clipping is performed.

The required absolute_error is a caller-established componentwise envelope in the prediction's units, broadcastable to its shape. The routine validates and returns it; it does not derive a remainder from two moments, certify physical support or bound floating-point error. For a downstream linear response A, propagate it as abs(A) @ envelope (flattening the Stokes/harmonic axes as needed). For example, zero is justified for this exactly quadratic synthetic kernel:

import syncmoments
import jax.numpy as jnp

model = syncmoments.QuadraticTaylorExpansion(
    harmonics=(1,), S0=jnp.zeros((1, 3)),
    dS=jnp.zeros((1, 3, 1)), ddS=2*jnp.ones((1, 3, 1, 1)),
)
prediction, error = model.checked_average([0.3], [[0.25]], absolute_error=0.)
# Each synthetic output is E[x**2] = 0.3**2 + 0.25 = 0.34.

rm_moments is the descriptive alias for gaussian_rm_cumulants. Neither name establishes Gaussianity. Use burn_depolarisation for a declared Gaussian screen, or screen_polarisation for a discrete weighted screen without that closure:

from syncmoments.rm import screen_polarisation, rm_moments

mean_rm, var_rm = rm_moments([-1., 1.])
P = screen_polarisation([1., -1.], [-1., 1.], jnp.sqrt(jnp.pi / 4))
# P is -1j: incident polarisation and RM are correlated across the two rays.

RM is in rad/m² and wavelength in metres. Incident complex P0 = Q + iU is a scalar or a vector matching the 1D RM samples; wavelengths have any shape, which is also the output shape. Optional weights are nonnegative relative masses. They are normalised through an exact power of two in float64 (float16, bfloat16 and float32 weights included), so a largest weight anywhere from 2^-1022 (about 2.2e-308) to 1.79e308 gives w / max(w); a subnormal largest weight raises. On XLA CPU a weight below 2^-1022 times the largest is flushed to zero, which drops at most n 2^-1022 of the mass. Second derivatives with respect to the weights use a closed-form second tangent of w / sum(w) (rm._sum_tangent, summed before one exact power-of-two scaling), and none is refused. Hessians in log-weights (w = exp(theta)) are exact in every mode down to max(w) of about 1e-307. For max(w) >= 2^-958, Hessians in w are exact where representable and inf with the exact sign where not, except that jacfwd over jacfwd returns NaN or inf once sum(w) < 2^-512. For max(w) < 2^-958 the exact Hessians in w exceed float64; they come out NaN, and jacrev over jacrev can give inf with the wrong sign (CHANGELOG, "Known limitations"). Rescaled weights give the Hessian exactly: at w it is lambda^2 times that at lambda w. Forward evaluation scans rays, avoiding a samples-by-wavelengths phase array; reverse-mode differentiation can still retain per-ray intermediates. JIT and gradients are supported. Invalid inputs and arithmetic overflow raise errors. Finite phase values alone do not guarantee accurate argument reduction for arbitrarily large phases; floating-point and sampling accuracy need separate checks in such regimes.

For fixed normalized weights, perturbing incident polarisation and RM gives |delta P| <= sum(w*|delta P0|) + 2*lambda²*sum(w*|P0|*|delta RM|). Use consistent intermediate values when also changing weights; their additional bound is max(|P0|)*sum(|delta w|). Sampling/quadrature uncertainty remains an external input. For internally distributed emission and frequency-dependent emitter spectra, use syncmoments.faraday.emission_polarisation below. Absorption and conversion require the separate transfer interfaces.

Independent symbolic derivation

derivation/ contains local Wolfram scripts that derive and check the helical orbit, retarded radiation, absolute harmonic powers and Stokes basis conventions from the classical equations. Run wolframscript -file derivation/verify.wls to regenerate the report. It also compares the complex harmonic amplitudes with direct Fourier integration of the original electric field, independently of this package's Python implementation.

Install and verify

Python 3.12 is the tested runtime. Install the release from PyPI:

python -m pip install syncmoments

To run the tests, install from a checkout (or the source distribution, which carries the tests):

python -m pip install -e '.[validation]'
python -m pytest tests -q

A plain run skips the 1661 tests marked slow (1573 of them in the automatic-differentiation transform matrix, 24 third-order derivative checks of the transfer, the others mostly the full-resolution benchmark); run them with python -m pytest tests -m slow or SYNCMOMENTS_RUN_SLOW=1.

Two environments are tested: JAX 0.10.0, Equinox 0.13.7, NumPy 2.3.5, SciPy 1.16.3, Matplotlib 3.10.8 and pytest 9.0.2; and JAX 0.10.2, Equinox 0.13.8, NumPy 2.5.3, SciPy 1.18.1 without mpmath. The high-precision references of the tests are literals (most generated by scripts/reference_constants.py); two transfer tests recompute theirs with mpmath on random slabs and are skipped without it. The benchmark tests use a verbatim, hash-checked copy of the manuscript's reference implementation and saved results (tests/model/reference), so they run without the manuscript repository; the source distribution includes the tests, scripts and documentation. The dependency ranges in pyproject.toml do not imply every version combination has been tested. SciPy is confined to validation; it supplies independent special-function, integration and matrix-exponential references.

Importing syncmoments enables JAX float64 globally for scientific accuracy. Import it before creating arrays; this is a documented compatibility side effect. For a source checkout without installation, PYTHONPATH=. remains supported.

PYTHONPATH=. python scripts/run_all.py                 # pytest, diagnostics, figures
PYTHONPATH=. python scripts/run_all.py --no-figures    # pytest and diagnostics

Printed diagnostics alone are not pass/fail physics tests. The runner first executes pytest, then reports script completion. Generated companion figures are package demonstrations; the manuscript uses its separate independent figure generators and captions. Do not substitute one set without checking assumptions.

Harmonic radiation and derivatives

from syncmoments.stokes import stokes_harmonic
from syncmoments.derivatives import derivative_spectra

I, Q, V = stokes_harmonic(10, 5.0, 0.785, 1.047, B=5e-6)
value, gradient, hessian = derivative_spectra(10, 5.0, 0.785, 1.047, B=5e-6)

B is in Gauss; harmonic power is per source time in erg/s/sr. The reference is vacuum radiation from a prescribed helical orbit, not a plasma-corrected or self-consistent particle trajectory. Without B, the API preserves its legacy dimensionless normalization. Differentiating that normalized response is a different operation: physical fixed-B energy derivatives must include the Lorentz-factor dependence of the gyrofrequency. These are fixed-harmonic quantities. Fixed-frequency predictions also need the moving harmonic line positions and a declared frequency/channel response.

The natural basis uses Q=parallel-perpendicular to the projected magnetic field and V=-2 Im(E_parallel E_perpendicular*). A fixed sky azimuth phi gives Q_sky=Q*cos(2phi), U_sky=Q*sin(2phi). The recurrence implementation includes the viewing-axis limit. Isotropic pitch angles alone do not force V to vanish.

Integer Bessel J uses a periodic integral with a contour shift to resolve exponentially small values. A concrete harmonic sets its static resolution; for traced orders the default is 2048 nodes. bessel_jn(..., n_nodes=...) exposes that static choice. bessel_jn returns NaN instead of an aliased answer for unresolved n+abs(x)>n_nodes/2, and for negative or non-integral n. bessel_jn_neighbours(n, x) returns J_{n-1}, J_{n+1} and J_n' from one quadrature on the order-n contour. It needs J_{n+1}, so it returns NaN for n+1+abs(x)>n_nodes/2, for n<1 and for non-integral n (at n_nodes=128, n=10, x=53.5 it gives NaN while bessel_jn is finite). Automatic derivatives through these quadratures lose accuracy at abs(x) << n: the third x-derivative of J_0 from bessel_jn_neighbours(1, x) at x=1e-6 (128 nodes) is -83.2 against the true 3.75e-7 (see the bessel_jn_neighbours docstring). HarmonicKernel differentiates the Bessel triple by the order recurrence instead. These guards are not accuracy certificates: validate resolution and derivatives over the actual domain. Modified K and continuum F/G use continuous hyperbolic integrals. Their public tail-bound helpers bound only the finite integral tail, excluding quadrature and roundoff. F(0)=G(0)=0; the continuum slopes diverge there, so derivative tests concern positive x.

sed.power_law_emissivity_abs(..., theta=...) gives directional ordered-field continuum emissivity. Its theta=None default additionally averages viewing angles/random field axes with the normalized sin(theta)/2 measure. Isotropic electron pitch angles alone do not perform that viewing-direction average.

Finite statistical response

import equinox as eqx
import jax.numpy as jnp
from syncmoments.expansion import build_expansion

response = build_expansion([1, 2, 5, 10], 5.0, 0.785, 1.047, B=5e-6)
mean_deviation = jnp.zeros(3)                    # gamma, alpha, theta
covariance = jnp.diag(jnp.array([0.01, 1e-4, 0.0]))
stokes = eqx.filter_jit(response)(mean_deviation, covariance)

build_expansion precomputes values, gradients and Hessians; repeated calls contract the same arrays with the supplied mean and covariance. Construction choices are fixed harmonic orders/reference/units; model arrays are ordinary Equinox leaves, so callers decide which parameters to differentiate. Runtime array shapes should remain stable when reusing a compiled response.

At the same Taylor degree and with the same statistics, moment and cumulant contractions are equal. A Gaussian has no cumulants beyond second order, but still has fourth and higher moments: a quadratic kernel average leaves those kernel terms out. The general scalar Bell helper can construct additional moments by setting unprovided cumulants to zero; that is an explicit closure, not a guarantee of a legal PDF. Vector contractions support at most four orders and reject higher orders. No universal pitch or energy width is safe: assess the kernel remainder at the actual reference, frequency and science precision.

apply_B uses the exact B-squared factor at fixed harmonic only when the remaining population parameters are fixed or B is independent of them. Correlated populations require conditional or mixed statistics. It does not supply physical units to a normalized spectrum or account for frequency shifts.

Correlated field, energy and orientation

An exact conditional B average leaves a function of the other parameters. apply_B does not infer that function: multiplying an already averaged spectrum by a global field factor assumes independence. The correlated alternative is response.mixed_average(field_moments, phase_moments, absolute_error=...). It expands only the known natural-basis kernel at fixed B0 and contracts its derivatives with unnormalised mixed moments. No unknown conditional function is differentiated. mixed_moments computes these statistics from paired samples:

from syncmoments import mixed_moments

def average_joint_population(response, offsets, B_over_B0, phi, weights, envelope):
    # response was built at B0>0; offsets use the SAME (gamma, alpha, theta)
    # reference. Angles are in radians; alpha is not the pitch cosine mu.
    W = B_over_B0**2
    field = mixed_moments(offsets, W, weights)
    phase = mixed_moments(offsets, W*jnp.exp(2j*phi), weights)
    return response.mixed_average(field, phase, absolute_error=envelope)

offsets has shape (samples,P), while B ratios, phi and optional relative electron-number weights have shape (samples,). Keep the samples paired to retain correlations. The helper returns (E[W], E[W*dq], E[W*dq*dq]), with shapes (), (P,), (P,P). Weights are normalized, W is not. The moments can instead come from a specified joint model. A complex zeroth moment may be zero while its higher mixed moments still contribute. Ordinary covariance of q and marginal B moments cannot in general determine these inputs.

mixed_average returns sky (I,Q,U,V) and a supplied componentwise envelope, both (harmonics,4); the older call and checked_average still return natural (I,Q,V). It keeps B-squared and the sky rotation exact, but truncates the q kernel at degree two. Thus E[B²*dq_i*dq_j] is retained even though it has total degree four in B and q. Moments of q must use alpha for this builder; converting only the mean/covariance from mu=cos(alpha) is not generally sufficient.

For a pointwise kernel remainder abs(R_s(q)) <= rho_s(q), use the weighted bound E[(B/B0)²*rho_s]; the natural Q bound controls each sky Q and U component. Uncertainty in the three supplied moment tensors adds abs(S0)*epsilon0 + abs(dS)@epsilon1 + abs(ddS):epsilon2/2 componentwise. These are mathematical input contracts. The method requires and returns an externally justified absolute_error in output units; it does not derive it, validate joint-moment realizability or physical support, or bound numerical and sample/quadrature errors. Zero is appropriate for an exactly quadratic kernel, not automatically for synchrotron radiation. For a subsequent linear observation map propagate the envelope with the elementwise absolute response matrix.

The new contraction is for fixed harmonic numbers. A frequency channel must also include B-dependent line positions and cannot use this simple B-squared factorisation. The power-law SED functions remain restricted continuum examples; they are not assumptions of the general kernel or mixed-moment method.

Finite statistics for a spectral fit

At a fixed expansion reference the response contracts a finite vector of joint statistics with known spectral coefficients. The same statistics apply to all harmonics; correlations do not require a separate unknown function for each harmonic. In the current three-coordinate quadratic builder, each moment block contains ten distinct monomials (one constant, three linear and six symmetric quadratic terms). The real field block and complex phase block thus contain thirty real statistical components before constraints and degeneracies.

These are not thirty unrestricted fit parameters. They must be consistent with one nonnegative joint electron-number distribution on the declared support. For example, abs(E[W*exp(2j*phi)]) <= E[W] for nonnegative W; analogous cross-moment constraints couple the two blocks. mixed_average checks neither their full realizability nor whether the spectral response identifies them. The method also does not establish the required remainder envelope from the retained moments. If higher moments are unknown, they must not silently be set to independent or Gaussian values to supply that envelope.

A frequency-channel fit needs coefficients built from the channel kernel, including moving line positions, with a single fixed population measure. Chromatic observing weights belong in those coefficients. The fitting layer must additionally specify the admissible joint statistics, number amplitude, measurement likelihood and propagated model discrepancy. An absolute error envelope is not automatically an independent Gaussian noise covariance. Weak or degenerate spectral responses constrain combinations of moments rather than each physical statistic separately; a small residual is not a remainder certificate or a reconstruction of the full population PDF.

The Faraday extension uses the joint distribution of emitting parameters and the intervening depth from each emitter to the observer. It includes mixed emission/rotation and retains spectrum-depth correlations. joint_faraday_average contracts a finite, frequency-independent moment matrix with intrinsic response coefficients and a finite phase series. The response and its truncation envelope are implemented; the application still supplies the intrinsic/channel basis, physical support, admissible moment set and likelihood. The Gaussian burn_depolarisation helper remains a declared factorised foreground model.

Finite joint response and spectral fits

syncmoments.model predicts channel-integrated Stokes spectra from a finite set of joint population moments (eq: finite joint response, eq: channel derivative coefficients) and fits those moments to channel data (extra eq: finite fit model). The population is never factorised unless an assumption object is declared, and every output carries an ErrorBudget whose terms are bound, estimate, measured, unbounded or not_applicable; a missing input is unbounded, never zero. docs/DESIGN.md describes the object model, the index layout, the free-parameter counts and the acceptance run. The four examples below are the printed output of scripts/model_examples.py (PYTHONPATH=. python scripts/model_examples.py, 1 min 43 s wall time with JAX 0.10.0 and 1 min 50 s with JAX 0.10.2 on an otherwise idle machine, most of it example (a); examples (b) to (d) live in scripts/model_examples_continuum.py and scripts/model_examples_fit.py); rerun it after any change and replace the quoted numbers. Example (a) compares with the manuscript's saved results in the verbatim copy tests/model/reference/validation, which the tests use too.

Example (a): harmonic kernel, the manuscript's smooth-channel benchmark

The setting of main-text Section 5.3.1: gamma0 = 20, B0 = 1 G, three bump channels at y_j = 2, 4, 8 (y = nu / nu_*, nu_* = e B0 / (2 pi gamma0 m_e c)) with 65 % support half-widths, the correlated population of eq: channel toy population on its full 16 x 16 latent and 64 x 64 angular grid (1,048,576 samples), harmonics through m = 40, and the truncation N = 2, L = 8 (1306 real moments).

from syncmoments.model.assumptions import independent_screen
from syncmoments.model.basis import build_basis
from syncmoments.model.channels import Channels
from syncmoments.model.harmonic import HarmonicKernel
from syncmoments.model.index import Truncation
from syncmoments.model.moments import JointMoments, Reference, Support
from syncmoments.model.predict import predict

y = np.array([2.0, 4.0, 8.0])
channels = Channels.bump(centres_hz=y * NU_STAR, widths_hz=0.65 * y * NU_STAR)
reference = Reference(20.0, 1.0, depth_ref=4 * S_DEPTH, scales=(20.0, 1.0, S_DEPTH))
support = Support(gamma=(16, 24), B=(0.8, 1.2), depth=(0, 10 * S_DEPTH), truncated=False)
basis = build_basis(HarmonicKernel(m_max=40), channels, Truncation(8, 8, 2), reference,
                    support=support)                       # 18 s, convergence="angular"
joint = JointMoments.from_samples(samples, basis.index, reference)
screened = JointMoments.from_samples(samples, basis.index, reference,
                                     parameter_map=independent_screen(basis.index),
                                     discrepancy="measured")
pred = predict(basis, screened, amplitude=1.0)
print(pred.summary()); print(pred.budget.total().kind); print(pred.budget.unbounded())

S_DEPTH = 1 / (2 (c / nu_*)^2) converts the manuscript's depth coordinate zeta to rad/m^2. With the joint moments the prediction reproduces the manuscript's saved finite response and its direct average:

required_m_max(support, channels) = 40
index: n0=486 n2=410 n_real=1306  build time 18.2 s
numerics: {'n_nodes': 256, 'm_range': (1, 26), 'width_ratio_min': 1.2999999999999998, 'cells': 64, 'route': 'product'}
population size: 1048576
joint moments:  max |pred - manuscript finite| / I_direct per channel: [1.52748222e-06 5.14411821e-07 1.39817674e-06]
joint moments:  max |pred - manuscript direct| / I_direct per channel: [0.00175259 0.00053129 0.00021887]

The second line is the manuscript's 1.75e-3 (Figure fig: full channel response, width 1.0, N = 2, L = 8). The build time includes the one convergence="angular" rebuild (the same route at doubled nodes for the numerical envelope). The product-versus-tensor route check is opt-in (build_basis(..., cross_route=True) or convergence="full") and is recorded as a named finite check, not a budget term. With independent_screen declared and discrepancy="measured" the moments are the factorised tensor and the budget carries the measured |m_joint - m_fac| propagated as N_src |C| Delta:

Prediction (eq: finite joint response): 3 channels, units per-electron channel Stokes per unit z-moment: erg/s/sr (unit_peak) or erg/s/sr/Hz (unit_integral); z = ((gamma-gamma0)/s_gamma, (B-B0)/s_B, (depth-depth_ref)/s_depth), Gauss, Hz, rad/m^2, amplitude 1
  channel 0: I=3.586096e-19 Q=1.882597e-20 U=-9.944719e-20 V=8.295205e-20
  channel 1: I=1.199577e-18 Q=-1.814062e-19 U=-3.997255e-19 V=3.017935e-19
  channel 2: I=3.461261e-18 Q=-7.788530e-19 U=-1.138940e-18 V=7.987727e-19
  error envelope unbounded; unbounded terms: basis_remainder, statistical_input, physical_kernel, depth_model, amplitude
    basis_remainder: unbounded -
    statistical_input: unbounded -
    physical_kernel: unbounded -
    harmonic_truncation: bound 0.000e+00
    excluded_tail: bound 0.000e+00
    numerical: estimate 7.826e-24
    screen_exponent: not_applicable -
    depth_model: unbounded -
    amplitude: unbounded -
    assumption:independent_screen: measured 2.776e-20
  assumptions: independent_screen
  note: incident polarisation not modelled
  note: moment-level discrepancy of kind 'measured' propagated as N_src |C| Delta (assumption)
budget.total().kind = unbounded
budget.unbounded() = ('basis_remainder', 'statistical_input', 'physical_kernel', 'depth_model', 'amplitude')
assumption 'independent_screen': kind measured; max |pred_screen - pred_joint| / I = 7.566e-03; max term / I = 7.731e-02; covered: True

The screen assumption moves the prediction by 7.6e-3 of channel I (the toy population correlates depth with energy and field); the measured term covers that shift. The numerical term is the per-column envelope |C_2x - C| of the product route against itself at doubled nodes, contracted with |m|. The total is unbounded because five inputs are missing. Declaring them for this synthetic population (E_phys zero because the harmonic reference is the truth, statistical_input zero because the moments are exact sums over the measure, amplitude_uncertainty=0, depth_model zero because the depths are the sample values) and probing the basis remainder gives the output below. The remainder inputs come from two measures. The order-3 derivative envelope H (at the reference point only) and the angular residual rho_ang are probed with RemainderInputs.from_samples on a separate 2-node (16-atom) Gauss-Legendre discretisation of the same toy model, which is not a subset of the population; this takes 35 s (147 s in 0.2.0, before the probe used HarmonicKernel.angular_taylor). The absolute moments <|P_l P_k| ||z||^3> of eq: local response remainder are recomputed on the predicted population itself (weighted sums, under a second); the 16-atom probe underestimates their row sums by the factors printed:

remainder probe (H at the reference, rho_ang on the 16-atom discretisation, absolute moments on the 1048576 population samples): 35.1 s, kind estimate
absolute-moment row sums, population / 16-atom probe (z_2, z_3): [1.732 1.589]
Prediction (eq: finite joint response): 3 channels, units per-electron channel Stokes per unit z-moment: erg/s/sr (unit_peak) or erg/s/sr/Hz (unit_integral); z = ((gamma-gamma0)/s_gamma, (B-B0)/s_B, (depth-depth_ref)/s_depth), Gauss, Hz, rad/m^2, amplitude 1
  channel 0: I=3.586096e-19 Q=1.882597e-20 U=-9.944719e-20 V=8.295205e-20
  channel 1: I=1.199577e-18 Q=-1.814062e-19 U=-3.997255e-19 V=3.017935e-19
  channel 2: I=3.461261e-18 Q=-7.788530e-19 U=-1.138940e-18 V=7.987727e-19
  error envelope (estimate): max 1.330e-15, 3.843e+02 of the largest channel I
    basis_remainder: estimate 1.330e-15
    statistical_input: bound 0.000e+00
    physical_kernel: bound 0.000e+00
    harmonic_truncation: bound 0.000e+00
    excluded_tail: bound 0.000e+00
    numerical: estimate 7.826e-24
    screen_exponent: not_applicable -
    depth_model: bound 0.000e+00
    amplitude: bound 0.000e+00
    assumption:independent_screen: measured 2.776e-20
  assumptions: independent_screen
  note: incident polarisation not modelled
  note: moment-level discrepancy of kind 'measured' propagated as N_src |C| Delta (assumption)
budget.total().kind = estimate
budget.unbounded() = ()
envelope covers |pred_screen - manuscript direct|: True

Every slot is now valued and the total is an estimate, but the envelope is 384 times channel I while the measured finite-vs-direct error is 1.75e-3 of I: the eq: local response remainder envelope sums absolute contributions over the 81 Legendre pairs with the Frobenius norm of the third-order tensor and the absolute moments (the depth displacement reaches |z_depth| = 3 in the manuscript's scale), so at L = 8 it exceeds the actual error by five orders of magnitude. H is evaluated at one point and rho_ang on the 16-atom measure, not the predicted population; neither is a supremum over the support. The amplitude slot is a zero bound: amplitude_uncertainty=0 declares the amplitude exact. E_phys and the excluded tail are declared by the caller. predict(samples=population, kernel=kernel) would run the probe itself, on the basis phase route and with the Support checks of the harmonic tail and of the declared-complete excluded tail, but it would probe H at every sample; at N >= 3 it leaves basis_remainder unbounded (order N + 1 is not among the validated orders) and the caller passes errors=RemainderInputs(...) with a supplied derivative_envelope.

Example (b): continuum kernel in the Galactic regime

At B0 = 5 uG, gamma0 = 3000 and channels between 0.1 and 3 GHz the harmonic index needed to reach the channels is of order 1e12, so the harmonic kernel is refused and ContinuumKernel (eq: directional continuum) is the only route. It ignores the pitch cosine mu, which forces isotropic_pitch (L_mu = 0), and it does not model V. The Faraday screen is a Burn screen (GaussianScreen(mean=30, sigma=5) rad/m^2, exact in tau for that screen, depth_degree=0), the population is a product measure local to the reference (energies 2700..3300 with N(gamma) ~ gamma^-2.5, field within 10 %, uniform pitch, one viewing angle of 60 degrees, one sky azimuth), and the remainder is probed on eight corner samples.

from syncmoments.model.kernels import ContinuumKernel, required_m_max
from syncmoments.model.phase import GaussianScreen

centres = np.geomspace(0.1e9, 3.0e9, 8)
channels = Channels.bump(centres_hz=centres, widths_hz=0.3 * centres)
reference = Reference(gamma0=3000.0, B0=5e-6, depth_ref=30.0, scales=(300.0, 5e-7, 5.0))
support = Support(gamma=(2700.0, 3300.0), B=(4.5e-6, 5.5e-6), depth=(0.0, 60.0))
screen = GaussianScreen(mean=30.0, sigma=5.0)
basis = build_basis(ContinuumKernel(), channels, Truncation(0, 2, 2, depth_degree=0),
                    reference, support=support, phase=screen)
moments = JointMoments.from_samples(samples, basis.index, reference,
                                    parameter_map=isotropic_pitch(basis.index),
                                    discrepancy="measured")
remainder = RemainderInputs.from_samples(probe, basis, kernel=ContinuumKernel(), phase=screen,
                                         segment_points=(1.0,), angular_residual="probe")
pred = predict(basis, moments, amplitude=1e20, errors=remainder,
               statistical_input=np.zeros(basis.index.n_real))

channels (GHz): [0.1   0.163 0.264 0.43  0.698 1.135 1.845 3.   ]
required_m_max(support, channels) = 2043405178518
critical frequency a_B gamma^2 (GHz) at gamma = 2700, 3000, 3300: [0.133 0.164 0.198]
Burn factor exp(-2 sigma^2 lambda^4) at the channel centres:
  [0.0000e+000 6.7602e-252 1.0801e-036 7.0752e-006 1.8302e-001 7.8414e-001
 9.6578e-001 9.9503e-001]
HarmonicKernel refused: required_m_max=2043405178518 lies outside the harmonic design regime (m <~ 16368); use ContinuumKernel with isotropic_pitch
index: components=('I', 'Q') n0=12 n2=12 n_real=36  build time 3.0 s
forced assumptions: ['independent_screen', 'gaussian_screen', 'isotropic_pitch']
notes: ('V not modelled [ContinuumKernel]', 'incident polarisation not modelled')
L_mu > 0 refused: ContinuumKernel requires the uniform_mu closure (isotropic pitch): use L_mu = 0
population size: 4096
remainder probe on 8 samples: 8.5 s, kind estimate
Prediction (eq: finite joint response): 8 channels, units per-electron channel Stokes per unit z-moment: erg/s/sr (unit_peak) or erg/s/sr/Hz (unit_integral); z = ((gamma-gamma0)/s_gamma, (B-B0)/s_B, (depth-depth_ref)/s_depth), Gauss, Hz, rad/m^2, amplitude 1e+20
  channel 0: I=2.321612e-01 Q=0.000000e+00 U=0.000000e+00 V=0.000000e+00
  channel 1: I=2.990974e-01 Q=-7.074508e-100 U=2.470894e-99 V=0.000000e+00
  channel 2: I=3.120780e-01 Q=-1.633457e-18 U=1.954271e-18 V=0.000000e+00
  channel 3: I=2.349223e-01 Q=8.577717e-06 U=-1.834397e-05 V=0.000000e+00
  channel 4: I=1.075459e-01 Q=2.940888e-03 U=1.145815e-03 V=0.000000e+00
  channel 5: I=2.343829e-02 Q=-5.012653e-03 U=5.761419e-03 V=0.000000e+00
  channel 6: I=1.700234e-03 Q=1.211610e-03 U=-6.457802e-04 V=0.000000e+00
  channel 7: I=2.302571e-05 Q=-2.624099e-06 U=-2.138403e-05 V=0.000000e+00
  error envelope unbounded; unbounded terms: physical_kernel, excluded_tail, depth_model, amplitude, assumption:isotropic_pitch, assumption:independent_screen, assumption:gaussian_screen
    basis_remainder: estimate 1.374e-02
    statistical_input: bound 0.000e+00
    physical_kernel: unbounded -
    harmonic_truncation: not_applicable -
    excluded_tail: unbounded -
    numerical: estimate 2.314e-08
    screen_exponent: bound 0.000e+00
    depth_model: unbounded -
    amplitude: unbounded -
    assumption:isotropic_pitch: unbounded -
    assumption:independent_screen: unbounded -
    assumption:gaussian_screen: unbounded -
  assumptions: independent_screen, gaussian_screen, isotropic_pitch
  note: V not modelled [ContinuumKernel]
  note: incident polarisation not modelled
  note: moment-level discrepancy of kind 'measured' not propagated: assumption 'isotropic_pitch' is forced by the basis and stays unbounded
budget.total().kind = unbounded
budget.unbounded() = ('physical_kernel', 'excluded_tail', 'depth_model', 'amplitude', 'assumption:isotropic_pitch', 'assumption:independent_screen', 'assumption:gaussian_screen')
assumption 'isotropic_pitch': kind unbounded, max value -
assumption 'independent_screen': kind unbounded, max value -
assumption 'gaussian_screen': kind unbounded, max value -
polarisation fraction |P|/I per channel: [0.000e+00 0.000e+00 0.000e+00 1.000e-04 2.930e-02 3.258e-01 8.075e-01
 9.357e-01]
max |pred - direct| / I per channel: [0.03 0.02 0.01 0.04 0.12 0.25 0.42 0.59]
sum of the valued terms / I per channel: [0.03 0.03 0.01 0.04 0.13 0.32 1.   4.37]
valued terms cover |pred - direct|: True

Four things to read off. physical_kernel stays unbounded: the continuum replacement is not bounded against the harmonic reference, and no channel comparison can be run at m ~ 1e12. All three assumptions are unbounded because no allowance is passed (example (d) shows assumption_allowances). The screen route's independent_screen and gaussian_screen have no declared discrepancy. isotropic_pitch is forced by the continuum kernel, and the retained L_mu = 0 moments cannot test it: the measured |m_joint - m_fac| on those rows is zero by construction, so predict keeps the term unbounded and notes that the moment-level discrepancy was not propagated. The Burn factor exp(-2 sigma^2 lambda^4) removes the polarisation below 0.4 GHz and is 0.995 at 3 GHz; the 94 % |P|/I printed for the top channel is the continuum kernel's intrinsic polarisation fraction on its tail. The upper channels sit on the exponential tail of the kernel (the critical frequency a_B gamma^2 of eq: directional continuum is 0.164 GHz at gamma0, 0.133 to 0.198 GHz over the gamma support), where the quadratic Taylor expansion in z_gamma is poor: the direct average differs by up to 59 % of I there, and the probed basis_remainder says so (the valued terms cover the difference in every channel, at the price of an envelope larger than the signal in the top two channels). Reducing that difference needs a narrower gamma support or a higher N.

Non-Gaussian independent screens

LaplaceScreen(mean, sigma) and GammaScreen(mean, sigma, shape=4.0, sign=1) (syncmoments.model.screens) are the illustrated screens of the manuscript (eq: illustrated screens). Each is an exact characteristic function with the same depth mean and standard deviation as GaussianScreen, applied line by line inside the channel sum, with depth_degree=0. At mean=0, sigma=1 and t = tau*sigma, the three weights are exp(-t**2/2), 1/(1 + t**2/2) and exp(-2j*t)*(1 - 1j*t/2)**-4. The Gamma screen adds the residual position angle -t + 2*arctan(t/2). Each screen records independent_screen plus laplace_screen or gamma_screen, all with unbounded discrepancy: the screen shape is a declared input, and zero covariance between depth and emission does not justify the factorisation. In predict these forced assumptions stay unbounded even when the moments carry a measured discrepancy, unless the caller passes an allowance through assumption_allowances (example (d)); without one a screen-route total is unbounded. tests/model/test_screens.py checks the weights against SciPy integrals of the screen densities, and the build_basis/predict route and direct_channel_average against a per-atom NumPy channel sum.

Example (c): fit_linear on synthetic continuum data

Synthetic data from the configuration of example (b) with twelve channels: d = A C m* + eta with Gaussian noise at 0.1 % of channel I, V masked (not modelled), isotropic_pitch declared (affine, so fit_linear applies), the amplitude fitted. StokesData takes the channel Stokes (n_ch, 4), per-row variances (channel-major) and a mask of rows to keep.

from syncmoments.model.fit.diagnostics import feasibility_checks, identifiability
from syncmoments.model.fit.linear import fit_linear
from syncmoments.model.fit.observation import StokesData

data = StokesData(stokes=stokes, noise=sigma.ravel() ** 2,
                  mask=jnp.array([True, True, True, False]))
result = fit_linear(basis, data, isotropic_pitch(basis.index))
report = identifiability(basis, data, isotropic_pitch(basis.index))
feasibility = feasibility_checks(result.moments, basis.index, support)

-- Truncation(L_mu=0, L_eta=0, N=0, depth_degree=0): n_real=3, free=2, data rows kept 36 of 48
rank=3 chi2=22.099 dof=33 converged=True amplitude=9.995439e+19
max |m_fit - m_true| = 6.164e-04; covariance is None: False
identifiability: rank=3 null_dim=0 weak(0.5)=()
feasibility: Necessary conditions only (eq: joint moment feasible set): 5 passed, 0 failed; 4 entries not computable from the supplied moments. Passing does not certify membership of conv{psi(x): x in D}.
  failed: ()
bias_bound kind: unbounded | prediction statistical_input kind: unbounded
-- Truncation(L_mu=0, L_eta=0, N=1, depth_degree=0): n_real=9, free=8, data rows kept 36 of 48
rank=6 chi2=28.914 dof=30 converged=True amplitude=9.974895e+19
max |m_fit - m_true| = 3.131e-02; covariance is None: True
identifiability: rank=6 null_dim=3 weak(0.5)=('<z_B>', 'Re <z_B e^{2i phi}>', 'Im <z_B e^{2i phi}>')
  null direction 0: (('amplitude', -0.08805817557044662), ('Re <e^{2i phi}>', -0.008213623149586166), ('<z_B>', 0.6412648032524726))
  null direction 1: (('amplitude', -0.048810333453748564), ('Re <e^{2i phi}>', -0.004821811624159282), ('<z_B>', 0.3554507991579424))
feasibility: Necessary conditions only (eq: joint moment feasible set): 4 passed, 3 failed: hermitian_psd[1,e^{2i phi}], support_bound[M2], abs_bound[f=1]; 6 entries not computable from the supplied moments. Passing does not certify membership of conv{psi(x): x in D}.
  failed: ('hermitian_psd[1,e^{2i phi}]', 'support_bound[M2]', 'abs_bound[f=1]')
bias_bound kind: unbounded | prediction statistical_input kind: unbounded
total time 5.6 s

The two printed null directions belong to an orthonormal basis of a three-dimensional null space whose singular values are at roundoff (1e-16). The SVD's choice of basis inside that space, and so the printed components, change with roundoff between runs, versions and environments (this output is from one JAX 0.10.0 run; later runs of the same code and JAX 0.10.2 print other vectors). The null space itself is the same: with JAX 0.10.0 the projectors onto it from 0.2.0 and 0.3.0 differ by 1.3e-15 in the 2-norm.

At N = 0 (three unknowns: the amplitude and the complex <e^{2i phi}>) the design has full rank, the moments are recovered to 6e-4 in z units and every computable necessary condition passes. At N = 1 the design is rank-deficient by three: the continuum kernel depends on (gamma, B) only through B_perp and B_perp gamma^2, so the energy and field displacements enter through one combination each in I and in P. The rank is decided on the column-equilibrated design G D, D = diag(1/||G_i||), so it does not depend on the Reference scales; weak(0.5) names <z_B> and <z_B e^{2i phi}>, and the null-space vectors (their first three components are printed) are in those equilibrated coordinates (report.column_scales holds D). The minimum equilibrated-norm solution then fails three necessary conditions (|<e^{2i phi}>| <= 1 among them), covariance is None, and the report says which combinations the data do constrain. bias_bound is unbounded at both orders because the data carry no declared discrepancy (at N = 1 the rank deficiency alone would also make it unbounded), and so is the fitted prediction's statistical_input.

Example (d): a caller-supplied screen allowance

assumption_allowances={name: ErrorTerm} (in predict and direct_channel_average) replaces the budget term of a named assumption, including one forced by the basis, with an allowance in Stokes units. The setting is that of example (b) with two rays of equal mass: the field is 4.7 uG with Faraday depth 25 rad/m^2 on one ray and 5.3 uG with 35 rad/m^2 on the other, so emission and screen are correlated. The screen is the exact EmpiricalScreen of the two depths, which leaves independent_screen as the only screen assumption. Its allowance is bounds.screen_factorisation_bound (detail-eq: screen factorisation error) in its per-line form, with the channel quadrature nodes nu_n as the lines: sigma_{j,n} = w_n |R_j(nu_n)| sigma_P(nu_n), where sigma_P(nu_n) is the ray standard deviation of the incident polarisation <K_Q e^{2i phi}> of each ray, and Phi_{j,n} = |<exp(i tau_n depth)>| over the rays. Each ray is uniform in mu by construction, so the continuum kernel's forced isotropic_pitch gets a declared zero.

from syncmoments.model.bounds import screen_factorisation_bound
from syncmoments.model.errors import ErrorTerm
from syncmoments.model.phase import EmpiricalScreen

screen = EmpiricalScreen(jnp.asarray([25.0, 35.0]), jnp.asarray([0.5, 0.5]))
basis = build_basis(ContinuumKernel(), channels, Truncation(0, 2, 2, depth_degree=0),
                    reference, support=support, phase=screen)
moments = JointMoments.from_samples(samples, basis.index, reference)
allowances = {
    "independent_screen": screen_factorisation_bound(sigma, phi, amplitude=1e20),
    "isotropic_pitch": ErrorTerm.declared_zero("each ray is uniform in mu by construction",
                                               shape=(channels.n_ch, 4)),
}
pred = predict(basis, moments, amplitude=1e20, statistical_input=0.0,
               amplitude_uncertainty=0.0, assumption_allowances=allowances)

forced assumptions: ['independent_screen', 'isotropic_pitch']
without allowances: budget.unbounded() = ('basis_remainder', 'physical_kernel', 'excluded_tail', 'depth_model', 'assumption:independent_screen', 'assumption:isotropic_pitch')
with allowances:    budget.unbounded() = ('basis_remainder', 'physical_kernel', 'excluded_tail', 'depth_model')
  assumption:independent_screen: bound; note: allowance supplied by caller: detail-eq: screen factorisatio...
  assumption:isotropic_pitch: bound; note: allowance supplied by caller: declared zero: each ray is uni...
screen error |direct exact - direct screen| / I: [0.007 0.011 0.    0.001 0.015 0.115 0.141 0.093]
screen_factorisation_bound / I: [0.03  0.045 0.064 0.104 0.181 0.252 0.165 0.094]
bound covers the screen error: True
direct screen route: assumption:independent_screen is bound
total time 3.7 s

Both assumption terms are now bounds with the note "allowance supplied by caller", and they leave the unbounded list. The total stays unbounded because basis_remainder (no RemainderInputs; example (b) shows the probe), physical_kernel, excluded_tail (Support.truncated is None) and depth_model are not supplied. The screen error is measured as the difference of two direct averages, the exact per-emitter phase against the factorised screen route; the bound covers it in every channel and is nearly attained in the top channel, since two rays with equal masses attain the Cauchy-Schwarz inequality at each node. The allowance holds for the stated sigma_P and Phi only; predict does not check them.

Per-variable caps and symmetry assumptions

Truncation(L_mu, L_eta, N, depth_degree=None, *, max_orders=(N_gamma, N_B, N_depth)) [extension] keeps only the rows with r <= N_gamma, s <= N_B and b <= N_depth in addition to the total cutoff (each cap an int or None). At L = 8, N = 2 the default keeps 1306 real moments, max_orders=(None, 1, None) keeps 1143 and (None, 0, None) keeps 735, and build_basis differentiates only in the directions those rows need. Caps that remove no row are dropped, so max_orders=(5, 5, None) at N = 2 is the default truncation. The retained multi-indices form a lower (downward-closed) set Lambda, and basis_remainder of a capped basis uses the margin S(Lambda) returned by syncmoments.model.lower_set_margin:

|f(z) - T_Lambda f(z)| <= sum_{beta in S} |z^beta| / beta! sup_{0<=t<=1} |d^beta f(p_k(t))|
p_k(t) = (z_1, ..., z_{k-1}, t z_k, 0, ..., 0),  k = first nonzero coordinate of beta

It follows from one-dimensional Taylor steps, one variable at a time; for the total-degree set, S is the shell |beta| = N + 1 of eq: local response remainder. S contains the minimal elements of the complement of Lambda and can contain more, and the minimal elements alone do not bound the remainder: for Lambda = {0,1}^2 and f = (x - x^2)(y - y^2) at (1, 1) the remainder is 1 while the sum over (2,0) and (0,2) with the supremum over the box is 0.5 (S adds (2,1)). RemainderInputs.from_samples probes the per-multi-index envelopes margin_H at finitely many path points, so the probed term is an estimate; with both envelopes supplied as arrays it is a bound for those inputs.

pitch_symmetric(index) and field_reversal_symmetric(index) [extension] declare symmetries of the population as ParameterMaps, alone or combined with assume:

from syncmoments.model import field_reversal_symmetric, pitch_symmetric

both = pitch_symmetric(basis.index).assume(field_reversal_symmetric(basis.index))
moments = JointMoments.from_samples(samples, basis.index, reference,
                                    parameter_map=both, discrepancy="measured")

Pitch symmetry is invariance under mu -> -mu with every other variable fixed; P_l(-mu) = (-1)^l P_l(mu), so every odd-l moment vanishes and its row is removed. Field reversal is invariance under B -> -B with the electrons fixed, which maps (mu, eta, phi) to (-mu, -eta, phi + pi) at fixed gamma, B and Faraday depth; e^{2i phi} is unchanged and P_l P_k changes sign for odd l + k, so those moments vanish and the finite V response is exactly zero. The reversal is local: reversing the field along the whole line of sight would also reverse the sign of the Faraday depth, which this declaration does not model. Free parameters at (2, 2, 2) are 153 without an assumption, 115 with pitch symmetry, 129 with field reversal and 103 with both; at (8, 8, 2), 1305, 769, 1065 and 649. Neither symmetry assumes independence between variables. The recorded discrepancy is unbounded unless it is measured or supplied; for a population without the symmetry the measured term equals |m_joint| on the removed rows. Combined with a factorisation, the symmetry is applied to each group, because a product measure is invariant exactly when every group marginal is: with mu and eta in different groups, field reversal drops the rows with odd l or odd k, so fully_independent with field reversal has 10 free parameters at (2, 2, 2) (12 without). Two consequences for closures combined with a symmetry:

• A delta or fixed_table closure whose own marginal is not symmetric is accepted. Only its even entries enter, so the map represents the symmetrised marginal, not the declared one. The assumption records list both declarations; nothing checks that they are compatible.
• A fixed_table combined with a symmetry takes the full table, the shape it has without the symmetry, both in one ParameterMap.build call and through assume; the two routes give the same map. The entries the symmetry drops (for example <P_1(mu)> on a mu singleton under field_reversal_symmetric with mu and eta in different groups) are ignored, whatever their values. The record keeps the full table and names the dropped entries in closure_kind, e.g. fixed_table (full table, symmetrised: dropped <P_1(mu)>)+field_reversal_symmetric. A table of the kept entries only raises ValueError naming the full shape.

No route averages the kernel over a declared distribution of B independent of the other variables: in the target populations the field strength is correlated with the electron energy, the Faraday depth and the field direction. max_orders with a small N_B reduces the cost of the field direction without that assumption.

HarmonicKernel(m_max, ..., chunk_budget=2**20, derivatives="analytic") bounds the working set of one lax.map step of the projections and the tail probe by chunk_budget Bessel integrand values and computes the (z_gamma, z_B) derivative tensors without tangents through the Bessel quadrature. derivatives="autodiff" nests jacfwd through the quadrature, the 0.2.0 structure on the new shared-contour Bessel rule (not bit-identical to 0.2.0); the two agree within 1e-12 per Stokes block. direct_channel_average caps the samples per step by the kernel's samples_per_step (102 for HarmonicKernel(40)), so batch_size is an upper cap. This trades time for memory: for the benchmark channels with S = 4096 samples and batch_size = 1024 the average takes 1.9 s instead of 1.3 s at the default budget, and a larger chunk_budget restores the old speed at the old memory. ContinuumKernel(..., chunk_budget=2**20) blocks its eta nodes under the same kind of budget. Forward tangents of an outer jacfwd are not counted in a budget unless the caller divides it with for_tangents, as the basis build does, and the nested-jacfwd remainder probe for the continuum and "autodiff" kernels (by 4^q at order q; the analytic harmonic kernel keeps its budget, because its tangents do not pass through the Bessel quadrature). The whole example (b) now peaks at 1.2 GB with JAX 0.10.0 and 1.8 GB with JAX 0.10.2 (6.7 GB and 27.1 GB before the budgets), example (d) at 0.85 GB and 1.1 GB (4.6 GB and 12.2 GB before). Peak RSS and time of build_basis for the benchmark basis of example (a) with convergence=False, on one shared machine:

Identifiable combinations of a continuum SED

reduce_response and fit_combinations (syncmoments.model.fit, since 0.4.0) [extension] keep three statements apart that a fit of the full moment vector mixes: exact redundancy of the response columns, the numerical rank of the design, and which combinations the data determine at the declared noise.

from syncmoments.model.fit import fit_combinations, reduce_response

reduction = reduce_response(basis, relations="continuum")
fit = fit_combinations(basis, data, reduction=reduction, max_sigma=0.1)
comparison = fit.against_truth(a_true, discrepancy=direct - forward, noise=noise)

reduce_response takes no data. For ContinuumKernel, whose kernel B Phi(B gamma^2) obeys B dK/dB - (gamma/2) dK/dgamma - K = 0, it groups the columns into C = H T with q = T a. For the other kernels relations="auto" keeps structural zeros only and relations="continuum" raises. fit_combinations fits beta = B x in a declared coefficient metric (Euclidean in the full coordinates by default), keeps the modes whose noise standard deviation 1/s_i is at most max_sigma, and reports the estimator, the resolution operator Pi, the covariance, the retained, weak, numerical-null and analytic-null directions, and four separate error terms: declared discrepancy, unresolved directions, approximate reduction and, for synthetic data, measured bias. An unconstrained term is unbounded. The fitted representative a_hat has only as many independent numbers as retained modes; it is not a moment vector and is not checked for realisability.

scripts/sed_reconstruction_example.py runs the continuum example of main.tex Section 5.3 (about 25 s): 60 full coordinates, 52 active columns and 30 groups before any data; numerical rank 26, whose smallest counted singular value is 1.058 times the cutoff; 9 combinations retained at max_sigma = 0.1; chi2 = 55.24 on 63 degrees of freedom. The mock data come from an independent NumPy/SciPy population integration. The rank, the retained set and beta depend on the channels, noise, masks, metric and reference scales. The guide page docs/guide/reconstruction.md ("SED fits of identifiable combinations" in the documentation) explains the example, its outputs and its comparison with the historical research run; docs/DESIGN.md Section 10.1 gives the design and Section 12.8 the limitations.

What the layer does not certify

• physical_kernel (E_phys) and excluded_tail (E_tail) are inputs.
• The numerical term is the per-column envelope |C' - C| of the basis against the same angular route at doubled nodes, contracted with |m|; it is an estimate, not a bound. The product-versus-tensor comparison is opt-in (cross_route=True, or convergence="full") and is a named finite check in provenance.finite_checks, not a budget term.
• The remainder probe evaluates envelopes at finitely many points (estimate) and computes absolute moments on the samples it is given; those must be the predicted population. from_samples(phase=None) probes the basis phase route. The probe is refused for N + 1 (for a capped truncation, any margin order |beta|) outside the orders validated against finite differences (N <= 2; the field certified_orders keeps its name for compatibility, and these are finite checks, not certificates); predict(samples=...) then leaves basis_remainder unbounded instead of raising.
• The zero harmonic_truncation bound assumes the population lies inside the declared Support. HarmonicKernel.truncation_error, predict(samples=..., kernel=...) and direct_channel_average check concrete samples: the term, and a declared-complete excluded_tail, become unbounded when an emitting sample lies outside. Without samples (or with traced samples) the hypothesis stays unchecked; the harmonic-tail note says so, and with traced samples so does the note of a declared-zero excluded_tail.
• A measured assumption term holds for the supplied population only. Assumptions forced by the basis (screen routes, the continuum kernel's isotropic_pitch) stay unbounded whatever the moments carry, because the retained moments cannot test them, unless the caller supplies an allowance through assumption_allowances. An allowance keeps its kind and value (not_applicable is refused); nothing checks that it holds. A moment-level discrepancy is attributed to the first assumption record only; assume() keeps it only when the other maps add no constraint.
• The amplitude slot is delta N_src (|C m| + e) with e the valued per-electron error slots (allowances divided by N_src), in predict and in direct_channel_average. A missing amplitude_uncertainty leaves it unbounded; a concrete 0 is a zero bound (amplitude declared exact); a negative value raises.
• A fit's bias_bound is |K| |delta| for the discrepancy the caller supplied, with delta = |R| E + |delta R| (|S_hat| + E): over the unknowns u for fit_linear, and over z for fit_bfgs, where it is an estimate linearised at the optimum (Gauss-Newton K = J^+ L^-1, prior and log-Jacobian curvature excluded, linearisation error not bounded). It is unbounded without a declared discrepancy, with an observing response but no response_uncertainty, and for a rank-deficient design; a linear fit's is an estimate when the plug-in |delta R| |S_hat| term is nonzero. The fit's statistical_input adds |dm/dx| bias_bound to the one-sigma noise and is unbounded when the bias bound is. A nodal fit has no bias bound and an unbounded statistical_input: the softmax gauge leaves its Fisher matrix singular. Rank decisions use the column-equilibrated design, so they do not depend on the Reference scales. Feasibility checks are necessary conditions, not a proof that the fitted moments come from a nonnegative population.
• A map with fitted hyper-parameters records NaN placeholders (closure kind "<kind> (fitted)"), so the moments' pytree structure does not depend on theta; the fitted values are in theta and in a provenance note of the FitResult and of its prediction. Every to_dict() (FitResult, Prediction, JointMoments, Provenance) writes non-finite floats as null, so its output is strict JSON.

Mixed emission and Faraday rotation

In the pure-rotation model, with no incident background, absorption, scattering or conversion, the observed polarisation is P(nu) = N_src * E[K_P(nu, p) * exp(2j*lambda**2*depth)]. The positive measure counts source electrons. Depth is in rad/m² and wavelength in metres; K_P is complex Q+iU in the observer's sky basis. Depth, source energy, field and orientation may all be correlated. Different emitters can have different spectral shapes; no common spectral factor is imposed.

emission_polarisation(emission, depths, lam, weights=None, source_column=1.) computes this discrete average. Depths have shape (n,). Emission is scalar, (n,) for wavelength-independent per-emitter values, or exactly (n, *lam.shape) for per-emitter spectra. The result has lam.shape. Relative nonnegative weights are normalized internally; the physical source_column is a separate finite nonnegative scalar. Complex emissivity is never normalized as a probability weight. A zero source column gives zero emission; the supplied measure must still be well defined. The original screen_polarisation retains its existing scalar/per-ray input contract and foreground interpretation.

faraday_depth_practical(n_e_cm3, B_par_uG, s_pc) in syncmoments.rm integrates from each node to the last node, using the rounded 0.812 coefficient and trapezoids. Positions increase towards the observer. Field reversals are allowed; depth need not be monotone. Add any exterior foreground depth to all nodes. For spatial quadrature, use weights proportional to source density times path-quadrature weights and set N_src to their total. Across unresolved rays, include each ray's column as well as its fixed nonnegative ray weight before normalizing. Chromatic instrumental weights belong in the observing response. This NumPy preprocessing routine is not differentiable; gradients with respect to supplied depths in the JAX average are supported. Path quadrature error is not inferred from the grid spacing.

For a finite spectral response, represent the intrinsic per-source kernel as sum_a c_a(lam)*psi_a(p) + r. The following synthetic exact-basis example keeps its intrinsic source error zero and bounds only the phase truncation:

import jax.numpy as jnp
from syncmoments.faraday import (
    emission_polarisation, joint_faraday_moments, joint_faraday_average,
)

lam = jnp.linspace(0.0, 0.8, 9)
x = jnp.array([-0.8, 0.2, 1.0])
depth = jnp.array([-0.4, 0.1, 0.6])
weights = jnp.array([1.0, 2.0, 4.0])
psi = jnp.stack([jnp.ones_like(x), x, jnp.exp(2j*x)], axis=-1)
c = jnp.stack([1 + lam, 0.3j*lam, -0.2 + lam**2], axis=-1)
emission = psi @ c.T
reference = emission_polarisation(emission, depth, lam, weights, source_column=5.)
M, absolute_next = joint_faraday_moments(
    psi, depth, 8, weights, reference_depth=0.1,
)
prediction, envelope = joint_faraday_average(
    c, M, lam, reference_depth=0.1, source_column=5.,
    absolute_next=absolute_next, source_error=0.,
)
assert jnp.all(jnp.abs(prediction-reference) <= envelope + 1e-13)

joint_faraday_moments takes basis values (n,n_basis) and a static nonnegative phase degree L, and returns M[a,b]=E[psi_a*(depth-reference_depth)**b] for b=0..L, plus A[a]=E[abs(psi_a)*abs(depth-reference_depth)**(L+1)]. The moment matrix can instead be fitted directly; its shape fixes the degree in joint_faraday_average. With coefficients (*lam.shape,n_basis), that routine returns the prediction and absolute envelope, both lam.shape:

prediction = N_src*exp(it*reference_depth) * sum_ab c_a*(it)^b/b!*M_ab
error      = N_src*(source_error + |t|^(L+1)/(L+1)! * sum_a |c_a|*A_a)
t          = 2*lam^2

source_error bounds E|r| per source and is required, as is absolute_next. For continuous populations they must have independent support/derivative or distributional justification; a sampled estimate alone is not a certificate. Retained moments do not determine these additional inputs. A complex zeroth moment can vanish while mixed moments still contribute. Joint cumulants can parametrise the same finite statistics; marginal depth cumulants alone cannot. Moment realizability and uncertainty, excluded tails, numerical error and physical-model discrepancy remain separate requirements. Large phase ranges can need more modes or direct quadrature; no universal low-order accuracy is claimed. Strong phases also require control of floating-point argument reduction.

The real-space source/depth pairing is sufficient for pure rotation; no recovery of a unique spatial geometry is implied. At fixed source measure and column, paired perturbations give |delta P| <= N_src*(E|delta K_P| + 2*lambda²*E[|K_P|*|delta depth|]), using the reference K_P in the second term. Different columns or weights add their normalization errors. For a frequency channel, integrate the combined emitted/rotated spectrum and propagate the envelope with the absolute observing weights; rotating an already integrated channel generally gives a different answer. These functions support JIT and gradients (phase degree is static when collecting moments). Test cases include the sinc slab, an independent transfer matrix exponential, swapped layers, correlated spectra and analytic gradients.

Physical transfer and reduced demonstrations

from syncmoments.los_moments import moment_driven_slab_cgs

# Local electron-number moments M_k=int N(gamma)(gamma-gamma0)^k dgamma.
# This narrow-population illustration does not certify an energy PDF closure.
gamma0, number_density = 2500.0, 1e-12
S = moment_driven_slab_cgs(
    1e8, gamma0, 5e-6,
    number_density, 0.0, number_density*50.0**2,
    number_density/gamma0,                 # illustrative inverse moment approximation
    3.0856775814913673e21,
    n_e=0.03, B_par=2e-6, phi=0.2,
)

The physical interface requires frequency Hz, fields Gauss, length cm and thermal/electron densities cm^-3. Output is erg/s/cm²/Hz/sr. It restores the separate emission and absorption prefactors; they do not cancel in a physical source function. Supply a measured or independently computed inverse moment when precision requires it. For finite support, pass the three endpoint terms [N(gamma)*(gamma-gamma0)^k]_lower^upper through boundary_terms; zero defaults assert that these vanish or are separately represented. Distributional jumps at physical hard cutoffs must be treated consistently.

The wrapper assumes isotropic ultra-relativistic electrons and a second-order local energy-kernel Taylor approximation. It uses h=gamma² for absorption; derivative_weighted(..., relativistic=True) evaluates the exact radial h=gamma*sqrt(gamma²-1) on a NumPy grid, with its own discretization error. The wrapper omits intrinsic V and alphaV; their propagated error is not bounded by this illustration. Generic mueller_matrix retains all Stokes absorption terms.

Cold-plasma mueller_rotation is the signed Stokes rV, twice the position-angle rate returned by rotation_coefficient. mueller_conversion and conversion_coefficient both return signed natural rQ, with no further factor two. Natural conversion mixes U and V; sky azimuth rotates rQ/rU. These leading high-frequency cold-dielectric coefficients need distribution-dependent replacements for hot or non-thermal plasma. The Gaussian Burn formula is an external-screen average, not internal emission and rotation; a uniform emitting slab has sinc depolarisation at finite Faraday depth.

moment_driven_slab remains a reduced-unit teaching interface and rejects physical Faraday arguments. Use the CGS wrapper to combine those effects. A uniform slab uses one augmented exponential. transfer_los handles ordered, nonuniform slabs; varying-medium discretization must be checked separately. A source-column coordinate rescaling by an exact power of two preserves zero-emissivity derivatives and keeps values correct up to max|eps ds| = 1.79e308. Derivatives use the linearity of the slab, out = Phi S + G eps: Phi, G and their derivatives in K and ds come from one 8x8 exponential, so forward and reverse Jacobians do not depend on the scale of eps (tested against SciPy from 1e297 to 1e307 for non-diagonal K and to 1.79e308 for K = a I). That exponential is a fixed Pade-13 scaling and squaring with ceil(log2(|A|_1 / theta_13)) squarings and no branch on a batch member's value, so reverse mode works through vmap and nested vmap over batched K or ds. Since 0.4.0 the value's 5x5 exponential uses the same method; against mpmath its relative error is at most 4.13 u max(|K ds|_1, 1) (u = 2^-53) on 214 test slabs, a measured constant, not a proven bound. 0.3.0 used jax.scipy.linalg.expm, whose floor count applies Pade 13 up to 2 theta_13. Its error grows with |K ds|_1 (measured up to 2.85e6 u max(|K ds|_1, 1)), and values differ from 0.3.0 by that error: 1.7e-9 relative on the sampled Faraday-dominated slabs, 7.2e-10 on a realistic moment_driven_slab_cgs spectrum, and 2.8e-5 of max|out| (3.4e-4 of Q) on an independent check slab at 36.3 MHz with L = 600 pc and |K ds|_1 = 8.7e4, where 0.4.0 is within 1.4e-12 of mpmath. The difference is not bounded by the small-depth values at large Faraday depth. d out / d ds is computed as Phi (eps - K S), with the error bound stated in the transfer_slab docstring. Second derivatives in every tested nesting, including derivatives of the primal output of jvp or vjp, chained slabs and transfer_los under vjp, agree with plain autodiff of Phi S + G eps and are symmetric up to max|eps ds| = 1e306: where the binary exponent of the source scale, plus terms for small ds and large |K ds|, reaches 896 (always in float32), the value carries the derivatives of Phi S + G eps from one extra 8x8 exponential, which about doubles the primal cost there (for a whole vmap batch or transfer_los path when any member needs it). Reverse Jacobians in K and ds are finite near 1.7e308. A derivative whose own value exceeds float64 overflows to inf without an error; only the value is checked. max_squarings=32 is the static exponential budget; tests include scalar optical depth 10^6, but do not guarantee arbitrary matrix depth, conditioning or gain. A non-finite slab is refused with an equinox error naming its cause: non-finite inputs, eps ds beyond float64, an exhausted max_squarings budget, or a result beyond float64 (accumulated intensity or extreme gain; 0.2.0 returned NaN or inf). The cores are compiled once per shape (equinox.filter_jit), so repeated eager calls do not recompile, and eager refusals raise EquinoxRuntimeError; under an outer jax.jit the error arrives as a JaxRuntimeError.

Magnus order two uses a running-prefix commutator, O(N) matrix work instead of an explicit O(N²) pair sum. It supports unequal slab widths and rejects orders other than 1 or 2. The commutator is deterministic ordering information, not a statistical cumulant. Ensemble transfer also needs source–propagator dependence.

Module boundaries

Compatibility and scientific scope

This review corrects signed Q, the conversion sign/factor/axis, fixed-B physical derivatives, and absolute isotropically averaged emissivity (an omitted angular probability factor made the legacy result twice too large). Physical predictions using those APIs will change. Dimensionless harmonic defaults are preserved. Invalid PDF/moment assumptions are not made valid by the corrected numerics.

The accompanying manuscript supplies a total-error propagation interface and bounded independent PDF benchmarks. A complete Galactic angular/field/energy-tail, plasma and instrument error budget is still an input to scientific sufficiency. fit_combinations reports which linear combinations of the truncated response a given data set constrains; identifiability of the population itself remains future work.

Citation and license

Please cite Zhang & Chluba, A statistical framework for synchrotron emission. MIT license; see LICENSE.