diff --git a/.github/workflows/axial-gauge-query.yml b/.github/workflows/axial-gauge-query.yml
new file mode 100644
index 000000000..1edaae840
--- /dev/null
+++ b/.github/workflows/axial-gauge-query.yml
@@ -0,0 +1,60 @@
+name: Axial gauge query mechanism
+
+on:
+  pull_request:
+    paths:
+      - src/prob4d/axial_gauge.py
+      - src/prob4d/axial_gauge_study.py
+      - tests/test_axial_gauge.py
+      - docs/axial-gauge-query.md
+      - .github/workflows/axial-gauge-query.yml
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+concurrency:
+  group: axial-gauge-query-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  controlled:
+    runs-on: ubuntu-latest
+    timeout-minutes: 15
+    env:
+      OPENBLAS_NUM_THREADS: "1"
+      OMP_NUM_THREADS: "1"
+      SOURCE_REVISION: ${{ github.event.pull_request.head.sha || github.sha }}
+    steps:
+      - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1
+        with:
+          ref: ${{ github.event.pull_request.head.sha || github.sha }}
+          persist-credentials: false
+      - uses: actions/setup-python@5fda3b95a4ea91299a34e894583c3862153e4b97 # v7.0.0
+        with:
+          python-version: "3.12"
+      - name: Install NumPy-only package and development checks
+        run: python -m pip install -e ".[dev]"
+      - name: Check focused source
+        run: |
+          python -m ruff check src/prob4d/axial_gauge.py src/prob4d/axial_gauge_study.py tests/test_axial_gauge.py
+          python -m mypy --python-version 3.12 src/prob4d/axial_gauge.py src/prob4d/axial_gauge_study.py
+          python -m pytest -q tests/test_axial_gauge.py
+      - name: Execute continuous-truth controlled experiment
+        run: |
+          python -m prob4d.axial_gauge_study \
+            --source-revision "$SOURCE_REVISION" \
+            --output outputs/axial-gauge-query-control-v1/result.json
+      - name: Verify numerical controls and show results
+        run: |
+          python - <<'PY'
+          import json
+          from pathlib import Path
+          report = json.loads(Path("outputs/axial-gauge-query-control-v1/result.json").read_text())
+          for row in report["records"]:
+              assert row["same_mean_max_abs_mm"] < 1e-10
+              assert row["same_covariance_max_abs_mm2"] < 1e-8
+              assert row["quadrature_refinement_max_abs_nll_nats"] < 1e-7
+              print(json.dumps(row, sort_keys=True))
+          print(report["claim_boundary"])
+          PY
diff --git a/docs/axial-gauge-query.md b/docs/axial-gauge-query.md
new file mode 100644
index 000000000..7b49b4b8e
--- /dev/null
+++ b/docs/axial-gauge-query.md
@@ -0,0 +1,189 @@
+# Nonlinear axial-gauge query laws
+
+**Experimental, conditional method and controlled study. Not a production provider update.**
+
+## Scientific purpose
+
+The observable-subspace factor avoids inventing information in a deficient
+Sim(3) gauge. The query-observability diagnostic then distinguishes direct
+geometric information from prior-mediated variance reduction. Both are local
+Gaussian constructions. A further issue remains: an unresolved rotation acts
+nonlinearly on an off-axis physical query. Retaining its tangent variance does
+not preserve the resulting distribution, and even exact Euclidean mean and
+covariance do not determine event probabilities.
+
+This module adds the following narrowly bounded method step:
+
+> Preserve the conditional angular law along an axial gauge orbit, update it
+> only through a declared likelihood, and propagate the same angular draw into
+> every coordinate of a joint physical query.
+
+It provides exact finite-quadrature posterior updates and pushforwards, a
+first/two-harmonic moment calculation, joint query covariance, mixture density,
+and halfspace probability. It changes neither stable provider-v2 export nor
+BayesianPhysTwin's complete-belief routing or exact fallback.
+
+## Conditional model and its limits
+
+Let eta denote the observable geometry/gauge and let theta denote rotation
+about its reference line. If the **entire positional likelihood**, including
+its covariance model, is invariant along that orbit, then
+
+    p(theta | eta, D) = p(theta | eta).
+
+This is a conditional statement. It does not imply that eta and theta are
+independent, that the angular prior is uniform, or that the marginal angular
+posterior is unchanged when the posterior over eta changes. For uncertain eta,
+a full belief must mix the eta-specific query laws with the posterior eta
+weights and preserve each prior conditional. This module operates at one fixed
+eta; it is not that complete-belief constructor.
+
+`AxialGaugeOrbit.from_line` checks only exact **mean geometry**, to a declared
+relative numerical tolerance. A local rank deficiency alone does not prove a
+global symmetry. Normals, appearance, anisotropic transformed source covariance,
+or other likelihood terms may carry angular information even for a collinear
+mean. Do not use the geometry check to erase that information.
+
+For a declared axis and fixed observable geometry, the optional
+`condition_on_correspondences` operation instead evaluates the actual nonlinear
+positional likelihood on the angular nodes. Weakly curved or off-axis geometry
+then updates the angle. The complete residual covariance may include
+cross-point dependence; it must be positive definite and fixed with respect to
+the angle. Unknown source coordinates, angle-dependent covariance, uncertain
+associations, and joint scale/rotation/translation inference are outside this
+likelihood's scope. Only observations allowed by the caller's causal/source
+protocol may be supplied.
+
+## Exact moment calculation
+
+For axis a, line point c, and nominal query point q, write
+
+    b = c + a a^T (q - c)
+    u = (I - a a^T) (q - c)
+    v = a cross u
+    q(theta) = b + u cos(theta) + v sin(theta).
+
+Let m1 = E[exp(i theta)] and m2 = E[exp(2 i theta)]. Then
+
+    E[q] = b + u Re(m1) + v Im(m1)
+
+and the covariance is `[u v] C [u v]^T`, where
+
+    C11 = (1 + Re(m2))/2 - Re(m1)^2
+    C22 = (1 - Re(m2))/2 - Im(m1)^2
+    C12 = Im(m2)/2 - Re(m1) Im(m1).
+
+Stacking the u and v vectors for several points gives their **joint** covariance.
+The same angular atom must be used for every point. For two antipodal probes,
+independent marginal draws would invent variance in their fixed midpoint.
+
+These are standard trigonometric-moment identities, not new circular-statistics
+theorems. `moments` is exact for the supplied finite angular law. A periodic
+quadrature approximates a continuous law; grid refinement is a numerical check,
+not a universal error bound.
+
+## Why covariance correction alone is insufficient
+
+A uniform angular law and equal masses at angles 0, 2pi/3, and -2pi/3 have the
+same first and second trigonometric moments: both vanish. For an off-axis probe
+at radius r, both therefore give the same Cartesian mean and covariance.
+Nevertheless, the probability that its radial reference coordinate is positive
+is 1/2 under the uniform law and 1/3 under the threefold law.
+
+Thus a Gaussian retaining even the exact mean and covariance cannot distinguish
+these two decisions. The threefold law's third harmonic matters. Under the
+threefold law, predicting 1/2 instead of 1/3 adds 1/36 to expected Brier loss.
+This counterexample is deterministic; the smooth controlled study adds small
+angular and readout noise and tests density and decision value separately.
+
+## Use with a verified rank-six factor
+
+```python
+import numpy as np
+from prob4d.axial_gauge import AxialGaugeOrbit, CircularQuadrature
+
+# Existing factor.chart.linearization maps the source cloud into reference units.
+reference_overlap = factor.chart.linearization.transform_points(source_overlap)
+orbit = AxialGaugeOrbit.from_line(reference_overlap)
+# Check that the full declared likelihood, not just the mean cloud, admits this
+# axial model. The caller supplies the conditional angle law in this axis convention.
+angular = CircularQuadrature(angles=registered_angles, weights=conditional_angle_mass)
+reference_queries = factor.chart.linearization.transform_points(source_query_points)
+query = orbit.pushforward(reference_queries, angular, noise_covariance=query_readout_covariance)
+
+# query.atoms contains complete point-major flattened query vectors.
+mean = query.mean
+joint_covariance = query.covariance
+probability = query.halfspace_probability(registered_normal, registered_threshold)
+```
+
+`from_line` fixes the largest-magnitude axis component positive. Angular priors
+must use that convention; changing the axis sign requires changing angle signs.
+The readout noise covariance is independent of the angular atom. The noiseless
+pushforward remains a discrete measure: `logpdf` refuses to label it a continuous
+density. A positive-definite readout covariance makes the Gaussian-mixture
+Lebesgue density well-defined.
+
+## Reproduce the controlled study
+
+```bash
+python -m pytest -q tests/test_axial_gauge.py
+python -m prob4d.axial_gauge_study \
+  --source-revision "$(git rev-parse HEAD)" \
+  --output outputs/axial-gauge-query-control-v1/result.json
+```
+
+The source-frozen protocol is embedded in the study module and copied with its
+SHA-256 into the result. It contains two independent seeds, five angular regimes,
+16,384 independent gauge/readout draws per seed and regime, a 50 mm off-axis
+probe, and 3 mm isotropic readout noise. Each draw is one statistical unit;
+coordinates and quadrature atoms are not independent units. Truth angles are
+sampled continuously by a separate simulator, never from the predictive grid.
+
+The three arms use exactly the same conditional geometry, prior information,
+and known readout noise:
+
+1. tangent Gaussian: nominal-position linearization with unwrapped angular variance;
+2. exact-moment Gaussian: the nonlinear mean and complete covariance, but Gaussian shape;
+3. axial-orbit mixture: the full shared-angle law with the same additive noise.
+
+The second and third arms have matched mean and covariance by construction. Their
+NLL and Brier differences test distribution shape, not an improved mean or a
+larger covariance. The primary comparison is their paired NLL, with a fixed
+halfspace Brier score as a separate decision diagnostic. Per-seed results,
+Monte Carlo intervals, numerical moment mismatch, and 512-to-1024-node density
+refinement are retained. No threshold is tuned to the simulated outcomes.
+
+Finalized numerical evidence and interpretation belong in
+`FlorianPfaff/BayesianPhysTwin-Paper`, not in a public provider-promotion claim.
+The local tests also cover likelihood-invariant prior preservation, genuine
+angular updating, cross-point noise, rigid-frame equivariance, geometry
+rejection, continuous wrapped-normal moment identities, and singular-density
+failure handling.
+
+## Relation to existing work and paper positioning
+
+Non-Gaussian rotational uncertainty, symmetry-aware pose inference, and circular
+moments are established topics. Relevant primary sources include:
+
+- Murphy et al., *Implicit-PDF*, ICML 2021:
+  https://proceedings.mlr.press/v139/murphy21a.html
+- Maken, Ramos, and Ott, *Estimating Motion Uncertainty with Bayesian ICP*, 2020:
+  https://arxiv.org/abs/2004.07973
+- Kurz et al., *Directional Statistics and Filtering Using libDirectional*:
+  https://arxiv.org/abs/1712.09718
+- Tuna et al., *X-ICP: Localizability-Aware LiDAR Registration*:
+  https://arxiv.org/abs/2211.16335
+
+The candidate Prob4D contribution is the narrower connection from deficient
+learned-window gauges to joint physical-query laws and measurable decision loss
+that survives exact mean/covariance matching. It is not the invention of a
+circular distribution, rotational symmetry, quadrature, or Bayesian conditioning.
+It specializes, rather than duplicates, the paper repository's general
+query-quotient lifting arguments.
+
+A stronger empirical paper still needs a separately frozen real-provider study,
+including a complete-belief integration and query-level proper-score/value
+comparison. The current PointWorld source-qualification route remains separate;
+this experiment opens no PointWorld, Flat'n'Fold, MotionCrafter, Deform360,
+BayesianPhysTwin, or Causal4D outcome and does not reopen any terminal cohort.
diff --git a/src/prob4d/axial_gauge.py b/src/prob4d/axial_gauge.py
new file mode 100644
index 000000000..2d1e496cf
--- /dev/null
+++ b/src/prob4d/axial_gauge.py
@@ -0,0 +1,311 @@
+"""Nonlinear query laws for an exactly unobserved rotation about a line.
+
+This experimental kernel is conditional on an observable gauge/geometry. It
+preserves a caller-supplied conditional angular law and can update it under an
+explicit positional likelihood. It does not complete a physical belief or
+authorize a real-provider update. Circular
+quadrature is shared across all query coordinates, retaining cross-point
+uncertainty. A nearly collinear cloud is not silently treated as an exact
+symmetry. See ``docs/axial-gauge-query.md`` for the statistical boundary.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from math import erfc
+from typing import Any, TypeAlias
+
+import numpy as np
+from numpy.typing import NDArray
+
+FloatArray: TypeAlias = NDArray[np.floating[Any]]
+
+
+def _array(value: object, *, name: str, ndim: int) -> FloatArray:
+    result = np.asarray(value, dtype=np.float64).copy()
+    if result.ndim != ndim or not np.all(np.isfinite(result)):
+        raise ValueError(f"{name} must be a finite {ndim}-dimensional array")
+    result.setflags(write=False)
+    return result
+
+
+def _mass(value: object, *, name: str) -> FloatArray:
+    result = _array(value, name=name, ndim=1).copy()
+    if not result.size or np.any(result < 0.0) or not np.any(result > 0.0):
+        raise ValueError(f"{name} must contain nonnegative mass with positive total")
+    result /= np.max(result)
+    result /= np.sum(result)
+    result.setflags(write=False)
+    return result
+
+
+def _covariance(value: object, *, dimension: int) -> FloatArray:
+    result = _array(value, name="noise_covariance", ndim=2).copy()
+    if result.shape != (dimension, dimension):
+        raise ValueError("noise_covariance has the wrong query dimension")
+    if not np.allclose(result, result.T, rtol=1e-12, atol=1e-12):
+        raise ValueError("noise_covariance must be symmetric")
+    result = (result + result.T) / 2.0
+    if np.min(np.linalg.eigvalsh(result)) < 0.0:
+        raise ValueError("noise_covariance must be positive semidefinite")
+    result.setflags(write=False)
+    return result
+
+
+@dataclass(frozen=True)
+class CircularQuadrature:
+    """A finite angular law; weights are masses, not unweighted density values."""
+
+    angles: FloatArray
+    weights: FloatArray
+
+    def __post_init__(self) -> None:
+        angles = _array(self.angles, name="angles", ndim=1)
+        weights = _mass(self.weights, name="weights")
+        if angles.shape != weights.shape:
+            raise ValueError("angles and weights must have the same nonempty shape")
+        object.__setattr__(self, "angles", angles)
+        object.__setattr__(self, "weights", weights)
+
+    def moment(self, order: int) -> complex:
+        """Return the trigonometric moment E[exp(i * order * angle)]."""
+        if isinstance(order, bool) or not isinstance(order, (int, np.integer)):
+            raise TypeError("order must be an integer")
+        return complex(np.sum(self.weights * np.exp(1j * int(order) * self.angles)))
+
+
+@dataclass(frozen=True)
+class GaussianQueryMixture:
+    """Joint query atoms with a shared additive Gaussian noise covariance.
+
+    Rows in ``atoms`` are complete query vectors, not independent point draws.
+    A zero covariance represents an exact discrete pushforward, for which
+    ``logpdf`` is deliberately unavailable. Noise is assumed independent of the
+    gauge atom; heteroscedastic or correlated gauge/readout noise is not modeled.
+    """
+
+    atoms: FloatArray
+    weights: FloatArray
+    noise_covariance: FloatArray
+
+    def __post_init__(self) -> None:
+        atoms = _array(self.atoms, name="atoms", ndim=2)
+        weights = _mass(self.weights, name="weights")
+        if atoms.shape[0] != weights.size or atoms.shape[1] < 1:
+            raise ValueError("atoms must have shape (number of masses, positive dimension)")
+        noise = _covariance(self.noise_covariance, dimension=atoms.shape[1])
+        object.__setattr__(self, "atoms", atoms)
+        object.__setattr__(self, "weights", weights)
+        object.__setattr__(self, "noise_covariance", noise)
+
+    @property
+    def mean(self) -> FloatArray:
+        return np.asarray(self.weights @ self.atoms, dtype=np.float64)
+
+    @property
+    def covariance(self) -> FloatArray:
+        centered = self.atoms - self.mean
+        result = centered.T @ (self.weights[:, None] * centered) + self.noise_covariance
+        return np.asarray((result + result.T) / 2.0, dtype=np.float64)
+
+    def logpdf(self, observations: FloatArray, *, batch_size: int = 256) -> FloatArray:
+        """Evaluate the continuous mixture density in bounded-memory batches."""
+        values = _array(observations, name="observations", ndim=2)
+        if values.shape[1] != self.atoms.shape[1]:
+            raise ValueError("observations have the wrong query dimension")
+        if isinstance(batch_size, bool) or not isinstance(batch_size, (int, np.integer)):
+            raise TypeError("batch_size must be an integer")
+        if batch_size < 1:
+            raise ValueError("batch_size must be positive")
+        try:
+            cholesky = np.linalg.cholesky(self.noise_covariance)
+        except np.linalg.LinAlgError as exc:
+            raise ValueError("logpdf requires positive-definite readout noise") from exc
+        inverse_root = np.linalg.solve(cholesky, np.eye(cholesky.shape[0]))
+        positive = self.weights > 0.0
+        atoms = self.atoms[positive] @ inverse_root.T
+        logweights = np.log(self.weights[positive])
+        constant = (
+            self.atoms.shape[1] * np.log(2.0 * np.pi)
+            + 2.0 * np.sum(np.log(np.diag(cholesky)))
+        ) / 2.0
+        result = np.empty(values.shape[0], dtype=np.float64)
+        for start in range(0, values.shape[0], batch_size):
+            whitened = values[start : start + batch_size] @ inverse_root.T
+            delta = whitened[:, None, :] - atoms[None, :, :]
+            logcomponents = logweights - np.sum(delta * delta, axis=2) / 2.0 - constant
+            maximum = np.max(logcomponents, axis=1)
+            result[start : start + batch_size] = maximum + np.log(
+                np.sum(np.exp(logcomponents - maximum[:, None]), axis=1)
+            )
+        return result
+
+    def halfspace_probability(self, normal: FloatArray, threshold: float) -> float:
+        """Return P(normal @ query > threshold), without Gaussianizing the law."""
+        direction = _array(normal, name="normal", ndim=1)
+        if direction.shape != (self.atoms.shape[1],) or not np.isfinite(threshold):
+            raise ValueError("normal/threshold must define a finite query-dimensional halfspace")
+        means = self.atoms @ direction
+        variance = float(direction @ self.noise_covariance @ direction)
+        if variance == 0.0:
+            return float(self.weights @ (means > threshold))
+        if variance < 0.0:
+            raise ValueError("negative projected noise variance")
+        standardized = (float(threshold) - means) / np.sqrt(2.0 * variance)
+        probabilities = np.array([0.5 * erfc(float(value)) for value in standardized])
+        return float(self.weights @ probabilities)
+
+
+@dataclass(frozen=True)
+class AxialGaugeOrbit:
+    """Rotations about a fixed line in the reference/world coordinate frame."""
+
+    center: FloatArray
+    axis: FloatArray
+
+    def __post_init__(self) -> None:
+        center = _array(self.center, name="center", ndim=1)
+        axis = _array(self.axis, name="axis", ndim=1)
+        if center.shape != (3,) or axis.shape != (3,):
+            raise ValueError("center and axis must have shape (3,)")
+        length = float(np.linalg.norm(axis))
+        if not np.isfinite(length) or length == 0.0:
+            raise ValueError("axis must have a finite positive length")
+        axis = axis / length
+        axis.setflags(write=False)
+        object.__setattr__(self, "center", center)
+        object.__setattr__(self, "axis", axis)
+
+    @classmethod
+    def from_line(
+        cls, reference_points: FloatArray, *, tolerance: float = 1e-10
+    ) -> AxialGaugeOrbit:
+        """Verify a nonzero exact line to a declared relative numerical tolerance.
+
+        ``tolerance`` bounds maximum transverse residual / RMS cloud radius.
+        It is not an observability threshold for a weakly curved object. Freeze
+        it before outcomes and do not relax it to turn weak geometry into symmetry.
+        The largest-magnitude axis coordinate is positive (deterministic sign).
+        """
+        points = _array(reference_points, name="reference_points", ndim=2)
+        if points.shape[0] < 2 or points.shape[1] != 3:
+            raise ValueError("reference_points must have shape (N>=2, 3)")
+        if not np.isfinite(tolerance) or tolerance <= 0.0 or tolerance >= 1.0:
+            raise ValueError("tolerance must be finite and lie strictly between zero and one")
+        center = np.mean(points, axis=0)
+        centered = points - center
+        radius = float(np.sqrt(np.mean(np.sum(centered * centered, axis=1))))
+        if not np.isfinite(radius) or radius == 0.0:
+            raise ValueError("a zero-extent cloud does not identify an axial stabilizer")
+        _, _, right = np.linalg.svd(centered / radius, full_matrices=False)
+        axis = right[0].copy()
+        if axis[int(np.argmax(np.abs(axis)))] < 0.0:
+            axis = -axis
+        transverse = centered - np.outer(centered @ axis, axis)
+        if float(np.max(np.linalg.norm(transverse, axis=1))) > tolerance * radius:
+            raise ValueError("reference geometry is not an exact line at the declared tolerance")
+        return cls(center=center, axis=axis)
+
+    def _components(
+        self, reference_queries: FloatArray
+    ) -> tuple[FloatArray, FloatArray, FloatArray]:
+        points = _array(reference_queries, name="reference_queries", ndim=2)
+        if points.shape[0] < 1 or points.shape[1] != 3:
+            raise ValueError("reference_queries must have shape (N>=1, 3)")
+        relative = points - self.center
+        parallel = np.outer(relative @ self.axis, self.axis)
+        cosine = relative - parallel
+        sine = np.cross(self.axis, cosine)
+        base = self.center + parallel
+        return base, cosine, sine
+
+    def positions(self, reference_queries: FloatArray, angles: FloatArray) -> FloatArray:
+        """Return (angle, point, coordinate) atoms with one shared rotation."""
+        theta = _array(angles, name="angles", ndim=1)
+        base, cosine, sine = self._components(reference_queries)
+        return (
+            base[None, :, :]
+            + np.cos(theta)[:, None, None] * cosine[None, :, :]
+            + np.sin(theta)[:, None, None] * sine[None, :, :]
+        )
+
+    def condition_on_correspondences(
+        self,
+        reference_points: FloatArray,
+        observed_points: FloatArray,
+        angular_prior: CircularQuadrature,
+        *,
+        noise_covariance: FloatArray,
+    ) -> CircularQuadrature:
+        """Bayes-update the angle using a complete, fixed positional likelihood.
+
+        Geometry is conditional/fixed. The covariance is for the point-major
+        stacked observation residual and may include cross-point dependence.
+        Exact-line observations leave the prior unchanged. Off-axis or weakly
+        curved observations may inform the angle and must not be discarded by
+        declaring a symmetry from a local rank threshold. Source uncertainty or
+        angle-dependent covariance requires a different likelihood.
+        """
+        observed = _array(observed_points, name="observed_points", ndim=2)
+        predicted = self.positions(reference_points, angular_prior.angles)
+        if observed.shape != predicted.shape[1:]:
+            raise ValueError("observed_points must match the reference point shape")
+        covariance = _covariance(noise_covariance, dimension=observed.size)
+        try:
+            root = np.linalg.cholesky(covariance)
+        except np.linalg.LinAlgError as exc:
+            raise ValueError("conditioning requires positive-definite residual noise") from exc
+        residuals = observed.reshape(1, -1) - predicted.reshape(predicted.shape[0], -1)
+        whitened = np.linalg.solve(root, residuals.T).T
+        log_likelihood = -0.5 * np.sum(whitened * whitened, axis=1)
+        positive = angular_prior.weights > 0.0
+        log_mass = np.full(angular_prior.weights.shape, -np.inf)
+        log_mass[positive] = np.log(angular_prior.weights[positive]) + log_likelihood[positive]
+        maximum = float(np.max(log_mass))
+        if not np.isfinite(maximum):
+            raise ValueError("correspondence likelihood has no finite prior-supported mass")
+        return CircularQuadrature(angular_prior.angles, np.exp(log_mass - maximum))
+
+    def moments(
+        self, reference_queries: FloatArray, angular_law: CircularQuadrature
+    ) -> tuple[FloatArray, FloatArray]:
+        """Exact finite-law joint moments from only the first two harmonics.
+
+        The mean is (N, 3); covariance uses point-major flattened coordinates.
+        This is a moment calculation, not a claim that the query is Gaussian.
+        """
+        base, cosine, sine = self._components(reference_queries)
+        first, second = angular_law.moment(1), angular_law.moment(2)
+        mean = base + first.real * cosine + first.imag * sine
+        coefficient_covariance = np.array(
+            [
+                [
+                    (1.0 + second.real) / 2.0 - first.real**2,
+                    second.imag / 2.0 - first.real * first.imag,
+                ],
+                [
+                    second.imag / 2.0 - first.real * first.imag,
+                    (1.0 - second.real) / 2.0 - first.imag**2,
+                ],
+            ]
+        )
+        basis = np.column_stack((cosine.reshape(-1), sine.reshape(-1)))
+        covariance = basis @ coefficient_covariance @ basis.T
+        return mean, (covariance + covariance.T) / 2.0
+
+    def pushforward(
+        self,
+        reference_queries: FloatArray,
+        angular_law: CircularQuadrature,
+        *,
+        noise_covariance: FloatArray | None = None,
+    ) -> GaussianQueryMixture:
+        """Retain the full shared-angle query law, optionally with readout noise."""
+        positions = self.positions(reference_queries, angular_law.angles)
+        dimension = 3 * positions.shape[1]
+        noise = np.zeros((dimension, dimension)) if noise_covariance is None else noise_covariance
+        return GaussianQueryMixture(
+            atoms=positions.reshape(positions.shape[0], dimension),
+            weights=angular_law.weights,
+            noise_covariance=noise,
+        )
diff --git a/src/prob4d/axial_gauge_study.py b/src/prob4d/axial_gauge_study.py
new file mode 100644
index 000000000..b3b329e82
--- /dev/null
+++ b/src/prob4d/axial_gauge_study.py
@@ -0,0 +1,220 @@
+"""Controlled nonlinear-gauge query study with continuous independent truth.
+
+Run with ``python -m prob4d.axial_gauge_study --output result.json
+--source-revision <40-hex-commit>``. No dataset, checkpoint, or downstream outcome
+is read. The fixed simulation protocol is included in the output with its hash.
+This is a conditional mechanism study, not real-provider calibration evidence.
+"""
+
+from __future__ import annotations
+
+import argparse
+import hashlib
+import json
+import platform
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+from .axial_gauge import AxialGaugeOrbit, CircularQuadrature, GaussianQueryMixture
+
+PROTOCOL: dict[str, Any] = {
+    "schema": "prob4d.axial-gauge-query-control-v1",
+    "evidence_kind": "controlled-conditional-mechanism",
+    "seeds": [20260830, 20260831],
+    "independent_draws_per_regime_per_seed": 16384,
+    "angular_nodes": 512,
+    "refinement_nodes": 1024,
+    "refinement_observations": 512,
+    "radius_mm": 50.0,
+    "readout_std_mm": 3.0,
+    "event": "noisy reference-frame query coordinate y > 0 mm",
+    "regimes": [
+        {"name": "narrow", "kind": "wrapped-normal", "std_rad": 0.05},
+        {"name": "moderate", "kind": "wrapped-normal", "std_rad": 0.6},
+        {"name": "broad", "kind": "wrapped-normal", "std_rad": 1.2},
+        {"name": "uniform", "kind": "uniform"},
+        {"name": "threefold", "kind": "threefold", "std_rad": 0.12},
+    ],
+    "arms": ["tangent-gaussian", "exact-moment-gaussian", "axial-orbit-mixture"],
+    "truth": "continuous angular draw plus independent isotropic Gaussian readout noise",
+    "statistical_unit": "one independent gauge-and-readout draw; not coordinates or atoms",
+    "paired_interval": "mean +/- 1.96 sample-standard-error; Monte Carlo interval only",
+    "closed_boundaries": ["real provider", "BayesianPhysTwin", "Causal4D", "protected cohorts"],
+}
+
+
+def canonical_json(value: object) -> str:
+    return json.dumps(value, sort_keys=True, separators=(",", ":"), allow_nan=False)
+
+
+def angular_rule(regime: dict[str, Any], nodes: int) -> CircularQuadrature:
+    """Periodic trapezoidal quadrature, independent of all simulated outcomes."""
+    theta = -np.pi + (np.arange(nodes, dtype=np.float64) + 0.5) * (2.0 * np.pi / nodes)
+    if regime["kind"] == "uniform":
+        weights = np.ones(nodes)
+    else:
+        std = float(regime["std_rad"])
+        centers = (
+            np.array([0.0, 2.0 * np.pi / 3.0, -2.0 * np.pi / 3.0])
+            if regime["kind"] == "threefold"
+            else np.array([0.0])
+        )
+        weights = np.zeros(nodes)
+        # Five images cover more than 8 standard deviations in the broadest regime.
+        for center in centers:
+            for winding in range(-2, 3):
+                weights += np.exp(-0.5 * ((theta - center + winding * 2.0 * np.pi) / std) ** 2)
+    return CircularQuadrature(theta, weights)
+
+
+def continuous_truth(
+    regime: dict[str, Any], rng: np.random.Generator, count: int
+) -> np.ndarray:
+    """Independent continuous simulator; no quadrature atoms or core rotations used."""
+    if regime["kind"] == "uniform":
+        theta = rng.uniform(-np.pi, np.pi, count)
+    elif regime["kind"] == "threefold":
+        centers = np.array([0.0, 2.0 * np.pi / 3.0, -2.0 * np.pi / 3.0])
+        theta = centers[rng.integers(0, 3, count)] + rng.normal(0, regime["std_rad"], count)
+    else:
+        theta = rng.normal(0, regime["std_rad"], count)
+    radius = PROTOCOL["radius_mm"]
+    positions = np.column_stack((np.zeros(count), radius * np.cos(theta), radius * np.sin(theta)))
+    return positions + rng.normal(0, PROTOCOL["readout_std_mm"], (count, 3))
+
+
+def normal_interval(values: np.ndarray) -> dict[str, float]:
+    mean = float(np.mean(values))
+    error = 1.96 * float(np.std(values, ddof=1)) / np.sqrt(values.size)
+    return {"mean": mean, "lower": mean - error, "upper": mean + error}
+
+
+def run_study(source_revision: str) -> dict[str, Any]:
+    if len(source_revision) != 40 or any(c not in "0123456789abcdef" for c in source_revision):
+        raise ValueError("source_revision must be a full lowercase hexadecimal Git commit")
+    orbit = AxialGaugeOrbit(np.zeros(3), np.array([1.0, 0.0, 0.0]))
+    query = np.array([[0.0, PROTOCOL["radius_mm"], 0.0]])
+    noise = np.eye(3) * PROTOCOL["readout_std_mm"] ** 2
+    event_normal = np.array([0.0, 1.0, 0.0])
+    records: list[dict[str, Any]] = []
+    for regime_index, regime in enumerate(PROTOCOL["regimes"]):
+        rule = angular_rule(regime, PROTOCOL["angular_nodes"])
+        nonlinear = orbit.pushforward(query, rule, noise_covariance=noise)
+        mean, covariance = orbit.moments(query, rule)
+        moment_gaussian = GaussianQueryMixture(mean.reshape(1, 3), np.ones(1), covariance + noise)
+        if regime["kind"] == "uniform":
+            angle_variance = np.pi**2 / 3.0
+        elif regime["kind"] == "threefold":
+            angle_variance = regime["std_rad"] ** 2 + 8.0 * np.pi**2 / 27.0
+        else:
+            angle_variance = regime["std_rad"] ** 2
+        tangent_covariance = noise.copy()
+        tangent_covariance[2, 2] += PROTOCOL["radius_mm"] ** 2 * angle_variance
+        tangent = GaussianQueryMixture(query, np.ones(1), tangent_covariance)
+        methods = {
+            "tangent-gaussian": tangent,
+            "exact-moment-gaussian": moment_gaussian,
+            "axial-orbit-mixture": nonlinear,
+        }
+        refinement = orbit.pushforward(
+            query, angular_rule(regime, PROTOCOL["refinement_nodes"]), noise_covariance=noise
+        )
+        log_scores: dict[str, list[np.ndarray]] = {name: [] for name in methods}
+        brier_scores: dict[str, list[np.ndarray]] = {name: [] for name in methods}
+        errors: dict[str, list[np.ndarray]] = {name: [] for name in methods}
+        event_values: list[np.ndarray] = []
+        refinements: list[float] = []
+        seed_records: list[dict[str, Any]] = []
+        for seed in PROTOCOL["seeds"]:
+            rng = np.random.default_rng(np.random.SeedSequence([seed, regime_index]))
+            truth = continuous_truth(regime, rng, PROTOCOL["independent_draws_per_regime_per_seed"])
+            event = (truth[:, 1] > 0.0).astype(np.float64)
+            event_values.append(event)
+            seed_metrics: dict[str, Any] = {"seed": seed, "arms": {}}
+            for name, method in methods.items():
+                nll = -method.logpdf(truth)
+                probability = method.halfspace_probability(event_normal, 0.0)
+                brier = (probability - event) ** 2
+                squared_error = np.sum((truth - method.mean) ** 2, axis=1)
+                log_scores[name].append(nll)
+                brier_scores[name].append(brier)
+                errors[name].append(squared_error)
+                seed_metrics["arms"][name] = {
+                    "nll_nats": float(np.mean(nll)),
+                    "brier": float(np.mean(brier)),
+                    "position_rmse_mm": float(np.sqrt(np.mean(squared_error))),
+                }
+            prefix = truth[: PROTOCOL["refinement_observations"]]
+            refinements.append(
+                float(np.max(np.abs(nonlinear.logpdf(prefix) - refinement.logpdf(prefix))))
+            )
+            seed_records.append(seed_metrics)
+        combined_nll = {name: np.concatenate(values) for name, values in log_scores.items()}
+        combined_brier = {name: np.concatenate(values) for name, values in brier_scores.items()}
+        record: dict[str, Any] = {
+            "regime": regime["name"],
+            "independent_draws": int(sum(values.size for values in event_values)),
+            "observed_event_frequency": float(np.mean(np.concatenate(event_values))),
+            "quadrature_refinement_max_abs_nll_nats": max(refinements),
+            "same_mean_max_abs_mm": float(np.max(np.abs(nonlinear.mean - moment_gaussian.mean))),
+            "same_covariance_max_abs_mm2": float(
+                np.max(np.abs(nonlinear.covariance - moment_gaussian.covariance))
+            ),
+            "arms": {},
+            "orbit_minus_moment_gaussian": {
+                "paired_nll_nats": normal_interval(
+                    combined_nll["axial-orbit-mixture"] - combined_nll["exact-moment-gaussian"]
+                ),
+                "paired_brier": normal_interval(
+                    combined_brier["axial-orbit-mixture"] - combined_brier["exact-moment-gaussian"]
+                ),
+            },
+            "per_seed": seed_records,
+        }
+        for name, method in methods.items():
+            record["arms"][name] = {
+                "nll_nats": float(np.mean(combined_nll[name])),
+                "brier": float(np.mean(combined_brier[name])),
+                "event_probability": method.halfspace_probability(event_normal, 0.0),
+                "position_rmse_mm": float(np.sqrt(np.mean(np.concatenate(errors[name])))),
+            }
+        records.append(record)
+    source_files = [Path(__file__), Path(__file__).with_name("axial_gauge.py")]
+    return {
+        "protocol": PROTOCOL,
+        "protocol_sha256": hashlib.sha256(canonical_json(PROTOCOL).encode()).hexdigest(),
+        "source_revision": source_revision,
+        "source_files_sha256": {
+            path.name: hashlib.sha256(path.read_bytes()).hexdigest() for path in source_files
+        },
+        "runtime": {"python": platform.python_version(), "numpy": np.__version__},
+        "records": records,
+        "claim_boundary": (
+            "Conditional controlled mechanism evidence with known angular laws and readout noise. "
+            "Not an end-to-end gauge estimator, real-provider result, "
+            "empirical calibration guarantee, "
+            "physical-state correction, intervention benefit, or deployment-safety result."
+        ),
+    }
+
+
+def main() -> int:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--output", type=Path, required=True)
+    parser.add_argument("--source-revision", required=True)
+    args = parser.parse_args()
+    if args.output.exists():
+        parser.error("refusing to overwrite existing evidence")
+    result = run_study(args.source_revision)
+    args.output.parent.mkdir(parents=True, exist_ok=True)
+    with args.output.open("x", encoding="utf-8") as handle:
+        json.dump(result, handle, indent=2, sort_keys=True, allow_nan=False)
+        handle.write("\n")
+    print(json.dumps({"output": str(args.output), "protocol_sha256": result["protocol_sha256"]}))
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/tests/test_axial_gauge.py b/tests/test_axial_gauge.py
new file mode 100644
index 000000000..4729b630a
--- /dev/null
+++ b/tests/test_axial_gauge.py
@@ -0,0 +1,221 @@
+"""Independent geometry, analytic-moment, and decision-law controls."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from prob4d.axial_gauge import AxialGaugeOrbit, CircularQuadrature, GaussianQueryMixture
+from prob4d.axial_gauge_study import PROTOCOL, angular_rule
+
+
+def uniform_rule(nodes: int = 256) -> CircularQuadrature:
+    return CircularQuadrature((np.arange(nodes) + 0.5) * 2 * np.pi / nodes, np.ones(nodes))
+
+
+def test_exact_line_is_fixed_but_off_axis_probe_moves() -> None:
+    center = np.array([4.0, 5.0, 6.0])
+    axis = np.array([1.0, 2.0, 3.0]) / np.sqrt(14.0)
+    line = center + np.linspace(-2.0, 2.0, 9)[:, None] * axis
+    orbit = AxialGaugeOrbit.from_line(line)
+    atoms = orbit.positions(line, uniform_rule().angles)
+    np.testing.assert_allclose(atoms, np.broadcast_to(line, atoms.shape), atol=1e-13)
+    probe = orbit.positions(np.array([[4.0, 6.0, 6.0]]), uniform_rule().angles)
+    assert np.ptp(probe[:, 0, 0]) > 0.1
+
+
+@pytest.mark.parametrize("scale", [1e-6, 1.0, 1e6])
+def test_relative_geometry_check_and_zero_extent(scale: float) -> None:
+    line = scale * np.array([[-1.0, 0.0, 0.0], [0.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
+    assert np.allclose(AxialGaugeOrbit.from_line(line).axis, [1.0, 0.0, 0.0])
+    line[1, 1] = scale * 1e-5
+    with pytest.raises(ValueError, match="not an exact line"):
+        AxialGaugeOrbit.from_line(line)
+    with pytest.raises(ValueError, match="zero-extent"):
+        AxialGaugeOrbit.from_line(np.zeros((3, 3)))
+
+
+def test_rigid_frame_equivariance() -> None:
+    rng = np.random.default_rng(7)
+    rotation, _ = np.linalg.qr(rng.normal(size=(3, 3)))
+    rotation[:, 0] *= np.linalg.det(rotation)
+    translation = rng.normal(size=3)
+    orbit = AxialGaugeOrbit(np.array([1.0, 2.0, 3.0]), np.array([1.0, 0.0, 0.0]))
+    transformed = AxialGaugeOrbit(rotation @ orbit.center + translation, rotation @ orbit.axis)
+    points = rng.normal(size=(4, 3))
+    angles = uniform_rule().angles
+    expected = orbit.positions(points, angles) @ rotation.T + translation
+    actual = transformed.positions(points @ rotation.T + translation, angles)
+    np.testing.assert_allclose(actual, expected, atol=1e-13)
+
+
+def test_shared_angle_preserves_distance_and_joint_covariance() -> None:
+    orbit = AxialGaugeOrbit(np.zeros(3), np.array([1.0, 0.0, 0.0]))
+    points = np.array([[0.0, 2.0, 0.0], [0.0, -2.0, 0.0]])
+    law = orbit.pushforward(points, uniform_rule())
+    positions = law.atoms.reshape(-1, 2, 3)
+    np.testing.assert_allclose(np.linalg.norm(positions[:, 0] - positions[:, 1], axis=1), 4.0)
+    np.testing.assert_allclose(law.covariance[1, 4], -2.0, atol=1e-13)
+    sum_y = np.array([0.0, 1.0, 0.0, 0.0, 1.0, 0.0])
+    assert abs(float(sum_y @ law.covariance @ sum_y)) < 1e-12
+    # Discarding cross-point covariance would invent variance four here.
+    assert float(sum_y @ np.diag(np.diag(law.covariance)) @ sum_y) == pytest.approx(4.0)
+
+
+def test_harmonic_moments_equal_direct_weighted_atoms() -> None:
+    rng = np.random.default_rng(8)
+    orbit = AxialGaugeOrbit(rng.normal(size=3), rng.normal(size=3))
+    rule = CircularQuadrature(rng.uniform(-np.pi, np.pi, 61), rng.uniform(size=61))
+    queries = rng.normal(size=(5, 3))
+    mean, covariance = orbit.moments(queries, rule)
+    mixture = orbit.pushforward(queries, rule)
+    np.testing.assert_allclose(mean.ravel(), mixture.mean, atol=1e-13)
+    np.testing.assert_allclose(covariance, mixture.covariance, atol=1e-13)
+
+
+def test_same_first_two_moments_do_not_determine_query_probability() -> None:
+    orbit = AxialGaugeOrbit(np.zeros(3), np.array([1.0, 0.0, 0.0]))
+    query = np.array([[0.0, 1.0, 0.0]])
+    uniform = orbit.pushforward(query, uniform_rule())
+    threefold = orbit.pushforward(
+        query, CircularQuadrature(np.array([0.0, 2 * np.pi / 3, -2 * np.pi / 3]), np.ones(3))
+    )
+    np.testing.assert_allclose(uniform.mean, threefold.mean, atol=1e-14)
+    np.testing.assert_allclose(uniform.covariance, threefold.covariance, atol=1e-14)
+    normal = np.array([0.0, 1.0, 0.0])
+    assert uniform.halfspace_probability(normal, 0.0) == pytest.approx(0.5)
+    assert threefold.halfspace_probability(normal, 0.0) == pytest.approx(1 / 3)
+
+
+@pytest.mark.parametrize("std", [0.05, 0.6, 1.2])
+def test_wrapped_normal_moments_against_continuous_analytic_formula(std: float) -> None:
+    rule = angular_rule({"kind": "wrapped-normal", "std_rad": std}, 512)
+    for order in [1, 2, 3]:
+        assert rule.moment(order) == pytest.approx(np.exp(-0.5 * order**2 * std**2), abs=1e-12)
+
+
+def test_threefold_and_uniform_continuous_moments_are_matched() -> None:
+    uniform = angular_rule({"kind": "uniform"}, 512)
+    threefold = angular_rule({"kind": "threefold", "std_rad": 0.12}, 512)
+    for order in [1, 2]:
+        assert abs(uniform.moment(order) - threefold.moment(order)) < 1e-13
+    assert abs(threefold.moment(3)) > 0.9
+    assert abs(uniform.moment(3)) < 1e-13
+
+
+def test_single_component_logpdf_matches_gaussian_and_batches() -> None:
+    rng = np.random.default_rng(9)
+    matrix = rng.normal(size=(3, 3))
+    covariance = matrix @ matrix.T + np.eye(3)
+    mean = np.array([1.0, 2.0, 3.0])
+    query = GaussianQueryMixture(mean[None, :], np.ones(1), covariance)
+    observations = rng.normal(size=(23, 3))
+    delta = observations - mean
+    expected = -0.5 * (
+        3 * np.log(2 * np.pi) + np.linalg.slogdet(covariance)[1]
+        + np.einsum("ij,jk,ik->i", delta, np.linalg.inv(covariance), delta)
+    )
+    np.testing.assert_allclose(query.logpdf(observations, batch_size=7), expected, atol=1e-13)
+    np.testing.assert_array_equal(
+        query.logpdf(observations, batch_size=1), query.logpdf(observations, batch_size=7)
+    )
+
+
+def test_zero_weight_components_and_tail_logpdf_stay_finite() -> None:
+    query = GaussianQueryMixture(np.array([[0.0], [1e6]]), np.array([1.0, 0.0]), np.eye(1))
+    result = query.logpdf(np.array([[1e3]]))
+    assert np.isfinite(result[0])
+    assert result[0] == pytest.approx(-0.5e6 - 0.5 * np.log(2 * np.pi))
+    assert query.halfspace_probability(np.array([1.0]), 0.0) == pytest.approx(0.5)
+    assert query.halfspace_probability(np.array([1.0]), 1.6448536269514722) == pytest.approx(0.05)
+
+
+def test_discrete_law_refuses_lebesgue_log_density() -> None:
+    law = GaussianQueryMixture(np.array([[0.0], [1.0]]), np.ones(2), np.zeros((1, 1)))
+    assert law.halfspace_probability(np.ones(1), 0.0) == pytest.approx(0.5)
+    with pytest.raises(ValueError, match="positive-definite"):
+        law.logpdf(np.zeros((1, 1)))
+
+
+@pytest.mark.parametrize("weights", [[0, 0], [-1, 2], [np.nan, 1], [np.inf, 1]])
+def test_bad_mass_is_rejected(weights: list[float]) -> None:
+    with pytest.raises(ValueError):
+        CircularQuadrature(np.array([0.0, 1.0]), np.asarray(weights))
+
+
+def test_mass_normalization_is_overflow_safe_and_copies_inputs() -> None:
+    weights = np.array([1e308, 1e308])
+    rule = CircularQuadrature(np.array([0.0, np.pi]), weights)
+    weights[:] = 0.0
+    np.testing.assert_array_equal(rule.weights, [0.5, 0.5])
+    assert not rule.weights.flags.writeable
+    assert rule.moment(0) == pytest.approx(1.0)
+    assert rule.moment(-1) == pytest.approx(rule.moment(1).conjugate())
+    with pytest.raises(TypeError):
+        rule.moment(True)
+
+
+@pytest.mark.parametrize("covariance", [np.array([[-1.0]]), np.eye(2), np.array([[np.nan]])])
+def test_invalid_noise_is_rejected(covariance: np.ndarray) -> None:
+    with pytest.raises(ValueError):
+        GaussianQueryMixture(np.zeros((2, 1)), np.ones(2), covariance)
+
+
+def test_query_dimension_and_argument_checks() -> None:
+    with pytest.raises(ValueError):
+        AxialGaugeOrbit(np.zeros(3), np.zeros(3))
+    with pytest.raises(ValueError):
+        CircularQuadrature(np.zeros(3), np.ones(2))
+    law = GaussianQueryMixture(np.zeros((1, 3)), np.ones(1), np.eye(3))
+    with pytest.raises(ValueError):
+        law.logpdf(np.zeros((2, 2)))
+    with pytest.raises(ValueError):
+        law.halfspace_probability(np.ones(2), 0.0)
+    with pytest.raises(ValueError):
+        law.halfspace_probability(np.ones(3), np.inf)
+    with pytest.raises(TypeError):
+        law.logpdf(np.zeros((2, 3)), batch_size=True)
+    assert PROTOCOL["statistical_unit"].startswith("one independent")
+
+
+def test_exact_line_likelihood_preserves_nonuniform_angular_prior() -> None:
+    points = np.array([[-1.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
+    orbit = AxialGaugeOrbit.from_line(points)
+    prior = angular_rule({"kind": "threefold", "std_rad": 0.12}, 512)
+    observed = points + np.array([[0.01, 0.02, -0.01], [-0.02, 0.01, 0.03]])
+    posterior = orbit.condition_on_correspondences(
+        points, observed, prior, noise_covariance=0.01 * np.eye(6)
+    )
+    np.testing.assert_allclose(posterior.weights, prior.weights, atol=1e-15)
+
+
+def test_off_axis_evidence_can_resolve_twist_instead_of_erasing_it() -> None:
+    orbit = AxialGaugeOrbit(np.zeros(3), np.array([1.0, 0.0, 0.0]))
+    points = np.array([[-1.0, 0.0, 0.0], [1.0, 0.03, 0.0]])
+    truth = 1.0
+    observed = np.array([[-1.0, 0.0, 0.0], [1.0, 0.03 * np.cos(truth), 0.03 * np.sin(truth)]])
+    posterior = orbit.condition_on_correspondences(
+        points, observed, uniform_rule(1024), noise_covariance=1e-6 * np.eye(6)
+    )
+    assert np.angle(posterior.moment(1)) == pytest.approx(truth, abs=1e-8)
+    assert abs(posterior.moment(1)) > 0.999
+    with pytest.raises(ValueError, match="not an exact line"):
+        AxialGaugeOrbit.from_line(np.vstack((points, np.zeros(3))))
+
+
+def test_correspondence_likelihood_handles_cross_point_covariance() -> None:
+    orbit = AxialGaugeOrbit(np.zeros(3), np.array([1.0, 0.0, 0.0]))
+    points = np.array([[0.0, 1.0, 0.0], [0.0, 0.8, 0.1]])
+    observed = np.array([[0.01, 0.9, 0.2], [-0.01, 0.7, 0.3]])
+    prior = uniform_rule(64)
+    covariance = np.eye(6) * 0.1 + np.ones((6, 6)) * 0.02
+    posterior = orbit.condition_on_correspondences(
+        points, observed, prior, noise_covariance=covariance
+    )
+    residual = observed.ravel() - orbit.positions(points, prior.angles).reshape(64, 6)
+    log_likelihood = -0.5 * np.einsum(
+        "ij,jk,ik->i", residual, np.linalg.inv(covariance), residual
+    )
+    expected = np.exp(log_likelihood - np.max(log_likelihood))
+    expected /= np.sum(expected)
+    np.testing.assert_allclose(posterior.weights, expected, atol=1e-14)