Faraday rotation and screens

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.

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.