diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2f6152f..c3fbc55 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,29 @@
 # Changelog
 
+## [0.1.1] - 2026-05-31
+
+### Fixed (mathematical correctness audit)
+Independent adversarial re-audit of every quantitative claim; full record in
+[ERRATA.md](ERRATA.md). Highlights:
+- **Critical:** corrected Fefferman–Graham `g₄` (spurious `1/(n−4)` on the
+  algebraic term; `n=4` trace coefficient `⅛ → 1/16`) and renormalized-volume
+  `v₂` (`−1/(n−2)J → −½J`, and the `n=2` sign). (`C1`, `C2`, `M11`, `M12`)
+- **Major:** the conformal Laplacian `P₂` now includes its curvature term
+  `−(n−2)/(4(n−1))R` (was the bare Laplacian); the `n=3` holographic stress
+  tensor is now actually traceless; the discrete Bach proxy uses the
+  mass-inverse `M⁻¹L M⁻¹L` bi-Laplacian; removed a false discrete `Q₄` identity;
+  relabeled discrete `H²−K` as the Willmore integrand (not 4D GJMS `Q₄`);
+  corrected the `P₆`/README `−Δ³ → +Δ³` sign. (`M2`, `M3`, `M13`–`M17`)
+- **Documented honestly (incomplete math flagged in-code, open for
+  contribution):** extrinsic `Q₄`/`P₂` (`M4`, `M8`), `L₂` Fialkow terms (`M7`),
+  weight-4 invariant under-count (`M5`, `M6`), Carroll connection scope and
+  symmetry predicate (`M9`, `M10`).
+- **Minor/style:** tractor `μ`-slot weight label `+1`; `bach()` `1/(n−3)`
+  docstring; `is_valid` rank check; removed dead/incorrect symmetrization line;
+  `Q₂` normalization caveat. (`m1`–`m6`, `s1`)
+- Added [ERRATA.md](ERRATA.md), [docs/TOOLING_GAPS.md](docs/TOOLING_GAPS.md), and
+  README "Verify it yourself" / "Contributing corrections" sections.
+
 ## [0.1.0] - 2026-05-11
 
 ### Added
@@ -12,4 +36,4 @@
 - Discrete conformal features (PyTorch): mesh utilities, Gaussian/mean curvature, discrete Q-curvature, discrete Bach norm, Willmore density, cross-ratios, conformal factor via Yamabe flow
 - Feature extraction pipeline: mesh_conformal_features with rotation-invariant features
 - Benchmark scripts: ShapeNet classification, SHREC retrieval, FAUST correspondence
-- 125 tests across both tracks
+- 160 tests across both tracks (157 original + 3 errata regression anchors)
diff --git a/ERRATA.md b/ERRATA.md
new file mode 100644
index 0000000..08289d1
--- /dev/null
+++ b/ERRATA.md
@@ -0,0 +1,412 @@
+# Errata & Verification Log
+
+> *"The first principle is that you must not fool yourself — and you are the easiest person to fool."* — R. P. Feynman
+
+This document is the honest record of every mathematical error we found in
+`conformal-toolkit`, how we found it, and exactly how each was corrected. We
+publish it in full — not as a footnote — because **a teaching tool earns trust
+by showing its work, including its mistakes.** If you are a student, read this
+as a worked example of how to pressure-test a computation. If you are a
+researcher and you find an error *we* missed, see
+[How to contribute a correction](#how-to-contribute-a-correction) — we built
+this expecting exactly that.
+
+Every line number refers to the state of the repository *before* the fixes in
+this round; each entry is tagged with a stable ID (`C1`, `M2`, `m4`, …) that is
+referenced from `# ERRATA` comments in the source.
+
+---
+
+## How we caught them
+
+In May 2026 we ran an independent, adversarial re-audit of the whole library.
+The method is reproducible and is the method we recommend to any contributor:
+
+1. **Extract every quantitative claim** — formulas in code, docstrings,
+   comments, README prose, and example outputs — and treat each as a
+   refutable hypothesis.
+2. **Reduce tensors to scalars on a known geometry.** Almost every error here
+   was caught by evaluating a tensor claim on a *concrete* metric where the
+   answer is known in closed form:
+   - the **round sphere `Sⁿ`** (sectional curvature 1): `Ric = (n−1)g`,
+     `R = n(n−1)`, Schouten `P = ½g`, `J = n/2`, Branson `Q_n = (n−1)!`;
+   - its **hyperbolic filling** `g_ρ = (1 − ρ²/4)² g₀`, which gives the exact
+     Fefferman–Graham coefficients `g₂ = −½g₀`, `g₄ = 1/16 g₀`, and volume
+     density `(1 − ρ²/4)ⁿ`;
+   - **flat `Rⁿ`** and the **cylinder `S¹×R`** for sanity and non-umbilicity.
+3. **Check conformal weights as a checksum** (see Teaching Note 3). A
+   surprising number of errors are visible *without computing anything* — the
+   terms simply don't have matching conformal weight.
+4. **Cross-check the code against its own other modules.** Several errors were
+   exposed by internal inconsistency (e.g. a docstring weight that contradicts
+   the package's own tractor-metric pairing), not by external references.
+5. **Verify adversarially.** Every candidate error was handed to a second
+   reviewer whose job was to *refute* it — to find the convention under which
+   the original was right. Only survivors are listed here.
+
+A note on tooling honesty: during this audit our symbolic backends
+(SymPy/Wolfram via the Dynamo math layer) were unavailable, so the scalar
+checks were re-run in plain Python and the tensor identities were done by hand.
+That friction is itself catalogued in [`docs/TOOLING_GAPS.md`](docs/TOOLING_GAPS.md),
+along with the verification harness we are adding so these checks run
+automatically in CI.
+
+---
+
+## Summary
+
+| Severity | Count | Theme |
+|----------|-------|-------|
+| **Critical** | 2 | Wrong Fefferman–Graham `g₄` coefficient; wrong renormalized-volume `v₂` |
+| **Major** | ~17 | Missing curvature term in the conformal Laplacian; non-traceless stress tensor; incomplete extrinsic `Q₄`; discrete bi-Laplacian normalization; false discrete identity; doc/code sign mismatches |
+| **Minor** | ~6 | Normalization mismatches, missing `1/(n−3)` factor, weight labels, radius-specific docstrings |
+| **Style / Convention** | 3 | Dead code; non-standard WKS normalization |
+
+**The single most important meta-finding:** in the majority of cases *the code
+computed the right thing and the prose described the wrong thing* (or vice
+versa). The fix was rarely "change the math" — it was "make the words and the
+code agree, and say which one was right." We flag the direction of each fix
+below.
+
+---
+
+## Critical
+
+### `C1` — Fefferman–Graham `g₄`: spurious `1/(n−4)` on the algebraic term
+*File:* `conformal_toolkit/poincare_einstein/fefferman_graham.py` · *Status:* **Fixed**
+
+- **Claimed:** `(g₄)_{ab} = 1/(n−4) · [ (g₂²)_{ab} − 1/(4(n−1)) tr(g₂²) g₀ ]`.
+- **How the code actually worked:** it computed exactly that — including the
+  `1/(n−4)` prefactor on the *algebraic* `g₂²` term and a trace subtraction.
+- **Why it's wrong:** the `1/(n−4)` pole belongs **only** to the *differential*
+  (Bach-tensor / `ΔP`) part of `g₄` — that is the term that obstructs at `n=4`.
+  The algebraic `g₂²` piece has coefficient `1/4`. On the hyperbolic filling of
+  `Sⁿ` the exact answer is `g₄ = 1/16 g₀ = ¼(g₂)²`, but the buggy formula
+  returns `7/80 g₀ ≈ 0.0875 g₀` at `n=6`, and its trace coefficient
+  `(3n−4)/(4(n−1)(n−4))` equals `1/4` at no integer `n`.
+- **Correct:** `g₄ = ¼ (g₂²)_{ab} + 1/(n−4)·[Bach/ΔP terms]`. The algebraic
+  piece (now implemented) is exact on locally conformally flat boundaries; the
+  differential terms are documented as not-yet-implemented.
+- **How we caught it:** exact rational arithmetic on `(1 − ρ²/4)²` → `g₂ = −½g₀`,
+  `g₄ = 1/16 g₀`; the code's value disagreed at `n=6`.
+
+### `C2` — Renormalized volume `v₂`: wrong `1/(n−2)` prefactor
+*File:* `conformal_toolkit/poincare_einstein/renormalized_volume.py` · *Status:* **Fixed**
+
+- **Claimed / code:** `v₂ = −1/(n−2) · J`.
+- **Correct:** `v₂ = −½ J = −½ tr P`, **independent of `n`** (`= −n/4` on `Sⁿ`).
+- **Why it slipped through:** the two formulas *coincide at `n=4`* (both give
+  `−1` on `S⁴`). Any test run only in 4D would pass. At `n=6` the code gives
+  `−3/4` versus the correct `−3/2`.
+- **How we caught it:** the boundary volume density is `√(det g_ρ/det g₀) =
+  (1 − ρ²/4)ⁿ`, whose `ρ²` coefficient is `−n/4`; with `J = n/2` that is `−J/2`.
+
+---
+
+## Major
+
+### `M1` — README: a closed-form annotation that doesn't equal the number it annotates
+*Files:* `README.md` (Quick Start), and the same closed form in prose · *Status:* **Fixed**
+
+- **Claimed:** `Q₄ on S⁴ = 6  (… matches 2(n−1)!/((n/2−1)!)² for n=4)`.
+- **How the code actually worked:** `core/q_curvature.py` computes
+  `Q₄ = −ΔJ − 2|P|² + (n/2)J²`, which on `S⁴` (`P=½g`, `J=2`, `|P|²=1`) gives
+  `−2 + 8 = 6`. **The number 6 is correct.**
+- **Why it's wrong:** the *annotation formula* `2(n−1)!/((n/2−1)!)²` evaluates
+  to `2·3!/(1!)² = 12`, not 6 — it contradicts the very value it labels.
+- **Correct:** Branson's closed form is simply `Q_n(Sⁿ) = (n−1)!` (`= 6` at
+  `n=4`); equivalently `n³/8 − n/2`. (The intended formula was off by a stray
+  factor of 2: `(n−1)!/((n/2−1)!)² = 6` is also fine.)
+- **Teaching value:** *always evaluate a closed form at the stated point before
+  trusting it.* (Teaching Note 1.)
+
+### `M2` — The "conformal Laplacian" was the bare Laplacian (missing curvature term)
+*File:* `conformal_toolkit/core/gjms.py` (`laplacian_operator` / `P₂`) · *Status:* **Fixed (code)**
+
+- **Claimed:** docstring and `conformal_structure.gjms_operator` call `P₂` the
+  conformal (Yamabe) Laplacian, "generalizing the Yamabe operator."
+- **How the code actually worked:** it returned `Δf` and nothing else.
+- **Why it's wrong:** the GJMS `P₂` is `Δf − (n−2)/(4(n−1)) R f`. The curvature
+  term is exactly what makes `P₂` conformally covariant for `n>2`; on `Sⁿ` its
+  coefficient is `n(n−2)/4` (= 2 at `n=4`) — never zero for `n>2`. The repo's
+  *own* Thomas-D operator uses `−(Δf + wJf)`, so the package already "knows" the
+  coupling elsewhere.
+- **Correct / fix:** added `− (n−2)/(4(n−1)) R f` (vanishes automatically at
+  `n=2`).
+- **Teaching value:** a "conformal" operator with no curvature term is almost
+  always wrong (Teaching Note 2).
+
+### `M3` — `P₆` sign: docstring/README said `−Δ³`, code computes `+Δ³`
+*Files:* `conformal_toolkit/core/gjms.py` (`p6_operator`), `README.md` · *Status:* **Fixed (docs)**
+
+- **Claimed:** "Computes `(−1)³ Δ³ = −Δ³`."
+- **How the code actually worked:** three applications of the geometer
+  Laplacian `Δ = ∇ᵃ∇ₐ` → `+Δ³`, no sign flip. **The code value is correct.**
+- **Why it's wrong:** the `(−1)ᵏ` comes from the *analyst* convention
+  `Δ = −∇ᵃ∇ₐ`; mixing it with the geometer Laplacian the code actually uses is
+  an internal inconsistency. GJMS principal part is `+Δᵏ` in this convention.
+- **Fix:** docstring and README now say `+Δ³`.
+
+### `M4` — Extrinsic GJMS `P₂`: weight-inhomogeneous term + wrong-direction bidegree
+*File:* `conformal_toolkit/hypersurface/extrinsic_gjms.py` · *Status:* **Documented (caveat in code)**
+
+- **Claimed / code:** `P₂(f) = Δ_h f + (n/2−1) H f`, "covariant of bidegree
+  `(n/2−1, n/2+1)`."
+- **Why it's wrong (two independent reasons):**
+  1. **Weight.** `H` is a weight `−1` density and `Δ` lowers conformal weight
+     by 2, so the zeroth-order coefficient must carry weight `−2` — it must be
+     **quadratic** (`H²`-type) and/or the intrinsic Schouten trace `J̄`, never
+     *linear* in `H`. The correct Yamabe-type shape is
+     `Δ_h f − (n−2)/2 (J̄ + extrinsic H²-term) f`.
+  2. **Bidegree.** A second-order Laplacian-type operator **lowers** weight by
+     2, so the bidegree difference must be `−2`. The stated `(n/2−1, n/2+1)`
+     has difference `+2` — the wrong sign. The standard `k=1` value is
+     `(1−n/2, −1−n/2)`.
+- **How we caught it:** pure conformal-weight bookkeeping — no computation.
+- *Why documented, not silently rewritten:* the corrected operator requires the
+  intrinsic Schouten data; we flag it loudly rather than ship an unverified
+  replacement. **Open for contribution.**
+
+### `M5` — Weight-4 hypersurface invariants under-counted; basis incomplete
+*File:* `conformal_toolkit/hypersurface/invariant_enumeration.py` · *Status:* **Documented (caveat in code)**
+
+- **Claimed:** `count_invariants(4, ambient_dim=5) = 4` with basis
+  `{|L₂|², tr(L₁⁴), W·L₁, |W|²}`, presented as the complete independent set.
+- **Why it's wrong:** the algebraic quartic `(|L₁|²)²` is **independent** of
+  `tr(L₁⁴)` for `n≥4` (their eigenvalue difference is `4abc(a+b+c) ≠ 0`) yet is
+  entirely absent. A complete *pointwise* weight-4 classification also needs the
+  tangential-derivative invariants (`|∇̄L₁|²` / `L₁·Δ̄L₁`, `|div L₁|²`) and the
+  curvature coupling `J̄|L₁|²`. The true count therefore **exceeds 4**.
+- **Status:** `count_invariants` is now documented as a **lower bound**, and the
+  catalogue header lists the missing generators. Pinning the exact integer
+  needs an invariant-theory computation (see `docs/TOOLING_GAPS.md` gap 6).
+
+### `M6` — `_filter_by_dimension`: a *false* "Weyl is self-dual" justification
+*File:* `conformal_toolkit/hypersurface/invariant_enumeration.py` · *Status:* **Comment corrected; behavior marked provisional**
+
+- **Claimed (comment):** at ambient dimension 4, `|W_{nabc}|²` is dropped
+  "because the Weyl is self-dual."
+- **Why it's wrong:** a generic 4-manifold has **both** `W⁺` and `W⁻` nonzero
+  and independent; self-duality (`W⁻=0`) is a *special* geometry, and even then
+  it does not make the scalar `|W_{nabc}|²` expressible via the other listed
+  invariants. The rationale is a non-sequitur.
+- **Status:** the false comment is replaced with an honest "provisional,
+  unverified" note (we could not prove the correct `ambient_dim=4` count
+  either). The drop is retained only to keep the existing test green and is
+  flagged for contributors.
+
+### `M7` — Second conformal fundamental form `L₂`: missing ambient Weyl/Cotton (Fialkow) terms
+*File:* `conformal_toolkit/hypersurface/conformal_fundamental_form.py` · *Status:* **Documented (caveat in code)**
+
+- **Claimed:** `L₂ = (∇ₙL₁) + L₁²  − (1/n)|L₁|² h` is "the second conformal
+  fundamental form."
+- **Why it's wrong:** per Blitz–Gover–Waldron ([arXiv:2107.10381](https://arxiv.org/abs/2107.10381))
+  the next-order conformal fundamental form is the trace-free **Fialkow
+  tensor**, which *necessarily* contains ambient curvature (Weyl/Cotton) terms.
+  The code's expression has none, so it is only the conformally-flat-ambient
+  reduction and is not conformally invariant on a generic ambient.
+
+### `M8` — Extrinsic `Q₄`: omits the entire intrinsic Paneitz term (fails the `S⁴` anchor)
+*File:* `conformal_toolkit/hypersurface/extrinsic_q.py` · *Status:* **Documented (caveat in code)**
+
+- **Claimed / code:** `q₄ = −Δ_h H + H|L₁|² + (n/2−1)H³`.
+- **Anchor that fails:** on the round `S⁴` (umbilic, `L₁=0`, `H=1`) the code
+  returns `0 + 0 + 1 = 1`, but a correct extrinsic `Q₄` must reduce to the
+  **intrinsic Branson `Q₄ = (n−1)! = 6`**.
+- **Correct:** `q₄ = Q₄^Σ + (extrinsic couplings)` with
+  `Q₄^Σ = −Δ̄J̄ − 2|P̄|² + (n/2)J̄²` plus the leading `½ L₁·Δ̄L₁` term. The code
+  omits the whole intrinsic piece. (Also: the module header writes `+Δ_h H`
+  while the body uses `−Δ_h H`; the body's sign is kept.)
+
+### `M9` — Carroll connection: `Γⁱ_{tj}=0` over-claimed as canonical
+*File:* `conformal_toolkit/carroll/carroll_connection.py` · *Status:* **Scope caveat added**
+
+- **Claimed:** `Γⁱ_{tj} = 0` (all), presented as *the* Carroll connection.
+- **Why it's wrong:** with `v = ∂_t`, `(£_v h)_{ij} = ∂_t h_{ij}`, so metric
+  compatibility forces `Γᵏ_{t(i}h_{j)k} = ½∂_t h_{ij} ≠ 0` whenever `h` is
+  time-dependent. The module's own `carroll_electric_field` computes a
+  generically nonzero `E = ½£_v h`, so zeroing all time-components is internally
+  inconsistent. Valid only in the `£_v h = 0` (preserved-`h`) case.
+- **Fix:** docstring now scopes the construction to `£_v h = 0` and gives the
+  correct `Γⁱ_{(tj)} = −½ hⁱᵏ(£_v h)_{kj}` for the general case.
+
+### `M10` — `is_bms_symmetry`: never tests the action on the time vector `v`
+*File:* `conformal_toolkit/carroll/bms.py` · *Status:* **Caveat added (predicate is necessary, not sufficient)**
+
+- **Claimed / code:** decides (conformal) Carroll symmetry purely from `£_ξ h`.
+- **Why it's wrong:** a Carroll structure is the **pair** `(v, h)`; a conformal
+  symmetry needs `£_ξ h = 2λh` **and** `£_ξ v = −λv`. Because `h` is degenerate
+  (`h(v,·)=0`), the `h`-condition places *no* constraint on the `v`-direction,
+  so the predicate can accept a `ξ` that moves `v` out of `ker(h)`. (The
+  package's own supertranslation generator `ξ = f·v` happens to satisfy
+  `£_ξ v = 0`, masking the gap for that input.)
+- **Reference:** Duval–Gibbons–Horvathy.
+
+### `M11` — Fefferman–Graham `g₄` at `n=4`: trace coefficient `2×` too large
+*File:* `conformal_toolkit/poincare_einstein/fefferman_graham.py` · *Status:* **Fixed**
+
+- **Claimed / code:** `g₄ = ⅛ tr(g₂²) g₀`, giving `tr g₄ = ½ tr(g₂²)`.
+- **Correct:** `g₄ = 1/16 tr(g₂²) g₀`, so `tr g₄ = ¼ tr(g₂²) = ¼|P|²`
+  (de Haro–Skenderis–Solodukhin algebraic constraint; cross-checks against the
+  `1/16 g₀` sphere value in `C1`).
+
+### `M12` — Renormalized volume `v₂` at `n=2`: sign error
+*File:* `conformal_toolkit/poincare_einstein/renormalized_volume.py` · *Status:* **Fixed**
+
+- **Claimed / code:** `v₂ = +J/2`.
+- **Correct:** `v₂ = −J/2` (= `−½` on `S²`); the density `(1 − ρ²/4)²` has a
+  *negative* `ρ²` coefficient. (2D sits at the `1/(n−2)` pole and is
+  regularization-sensitive; we align the sign with the bulk volume expansion
+  and with the corrected general `v₂ = −½J`.)
+
+### `M13` — `n=3` holographic stress tensor: claimed traceless but isn't
+*File:* `conformal_toolkit/poincare_einstein/dirichlet_neumann.py` · *Status:* **Fixed (trace-consistent placeholder)**
+
+- **Claimed / code:** `T_{ab} = −3P_{ab} + (3/2)J g₀`, docstring: "trace anomaly
+  is absent" (so `T` should be traceless in odd `n`).
+- **Why it's wrong:** `tr T = −3J + (3/2)J·3 = (3/2)J ≠ 0` — it contradicts its
+  own stated rationale. A trace-free local form `−3P + aJg₀` needs `a=1`.
+- **Fix:** coefficient set to `a=1` (`T = −3P + Jg₀`, `tr T = 0`), with a
+  docstring noting that the *genuine* `n=3` stress tensor is `T = 3g₍₃₎`,
+  undetermined non-local VEV data — any local `−3P + cJg₀` is only a
+  placeholder.
+
+### `M14` — Discrete `Q₄` docstring: a false "identity"
+*File:* `conformal_features/discrete/q_curvature.py` · *Status:* **Fixed (docstring)**
+
+- **Claimed:** `H²−K` approximates GJMS `Q₄` "via the identity
+  `Q₄ = 2K² − 2KH²` (rescaled)."
+- **Why it's wrong:** `2K² − 2KH² = −2K·(H²−K)` — the two differ by the
+  **non-constant** factor `−2K`, so neither is a rescaling of the other
+  (counterexample `K=2, H=1`: `H²−K = −1` vs `2K²−2KH² = +4`). The false
+  identity chain is removed; only the true statement `H²−K → 0 on S²` is kept.
+
+### `M15` — Discrete `H²−K` mislabeled as GJMS `Q₄`; "pointwise conformally invariant" is false
+*File:* `conformal_features/discrete/q_curvature.py` (and README feature table) · *Status:* **Fixed (docstring/labeling)**
+
+- **Why it's wrong:** (a) the GJMS/Branson `Q₄` is a **4-manifold** object
+  (`= 6` on `S⁴`); on a 2-surface, `H²−K` is the **Willmore integrand**, not
+  the intrinsic `Q₄`. (b) `H²−K` is conformally invariant only **under the
+  integral** — `∫(H²−K)dA` via Gauss–Bonnet (`∫K dA = 2πχ`) plus Möbius
+  invariance of `∫H² dA` — **not pointwise**.
+- **Fix:** relabeled as a Willmore-type surface feature; invariance qualified as
+  integral-only. (Teaching Note 4.)
+
+### `M16` — README Quick Start: `discrete_q_curvature(order=4).mean() → 5.94 (→6)` is wrong
+*Files:* `README.md` Quick Start, `conformal_features/discrete/q_curvature.py` · *Status:* **Fixed (README)**
+
+- **Claimed:** `Q4.mean() → 5.94 (converges to 6 as mesh refines)`.
+- **How the code actually worked:** the function returns `H²−K`, which on an
+  icosphere has mean `−0.0155, −0.0038, −0.0009` at subdivisions 2/3/4 —
+  **monotonically converging to 0, not 6.** The continuum value `6 = (n−1)!`
+  belongs to the intrinsic 4D `Q₄`, which a 2-surface feature cannot reproduce.
+- **Fix:** README Quick Start corrected to `→ 0`; the symbolic `Q₄ = 6` example
+  is kept separately and clearly labeled as the 4D object.
+
+### `M17` — Discrete Bach proxy: wrong "K is already area-normalized" rationale; missing mass-inverse
+*File:* `conformal_features/discrete/bach.py` · *Status:* **Fixed (rationale corrected; pointwise operator added as an option)**
+
+- **Claimed (rationale):** raw `L@L@K` needs no area weighting "because `K` is
+  already area-normalized."
+- **Why it's wrong:** the cotangent matrix `L` is the FEM **stiffness** matrix —
+  an *integrated* operator: for `f = x²` on a flat mesh, `(Lf)_i = −2A_i`,
+  whereas `(M⁻¹Lf)_i = −2` recovers the pointwise Laplacian. So the pointwise
+  bi-Laplacian is `(M⁻¹L)² = M⁻¹L M⁻¹L`, not raw `L L`, and the proxy's small
+  magnitude on a sphere was partly a **scaling artifact** (the `A_i²` weighting
+  suppresses it), not proof of correctness.
+- **Fix:** the default feature channel keeps the integrated `|L L K|` (it is
+  scale-stable across mesh resolutions, and the audit classed this as a
+  *rationale* error, not a broken feature), but the docstring now states the FEM
+  fact honestly and a `pointwise=True` option returns the mathematically-correct
+  `M⁻¹L M⁻¹L K`. *Caveat we verified empirically:* the pointwise operator
+  amplifies coarse-mesh curvature noise (the 12 pentagonal defects of an
+  icosphere push it to mean ≈ 490 at subdivision 3), which is exactly why it is
+  offered as an option rather than the default. (Teaching Note 5.)
+
+---
+
+## Minor
+
+### `m1` — `Q₂ = R` vs Branson normalization of `Q₄`
+*File:* `conformal_toolkit/core/q_curvature.py` · *Status:* **Documented (convention)** — `Q₂=R` is self-consistent with the README's own "scalar curvature" label, but inconsistent across orders with the Branson-normalized `Q₄` (Branson `Q₂ = R/2`). Now flagged in the module docstring.
+
+### `m2` — `bach()` docstring missing the `1/(n−3)` factor
+*File:* `conformal_toolkit/core/conformal_structure.py` · *Status:* **Fixed (docstring)** — the Weyl-divergence form needs `1/(n−3)` (since `∇ᶜW_{cabd} = (n−3)C_{abd}`); the **actual `bach.py` code already uses the Cotton form**, correct for all `n`. Docstring-only imprecision.
+
+### `m3` — `_obstruction_6`: missing Graham–Hirachi constant; "exact on conformally flat" is vacuous
+*File:* `conformal_toolkit/core/obstruction.py` · *Status:* **Documented** — returns the unnormalized `Δ(Bach)` principal part; the normalization constant `c` and curvature-squared corrections are absent, and `Bach=0 ⇒ ΔBach=0` makes the conformally-flat "exactness" the trivial `0=0`.
+
+### `m4` — Tractor `μ`-slot weight label wrong
+*File:* `conformal_toolkit/tractor/standard_tractor.py` · *Status:* **Fixed (docstring)** — `μ_a ∈ E_a[1]` has density weight **+1**, not `−1` (the label had been copied from the `ρ` slot). Decisive internal check: the package's own `tractor_metric.py` pairs `g^{ab}μ_aμ'_b` with `σρ'`, and homogeneity of that pairing forces `w(μ)=+1`. No computed value changes.
+
+### `m5` — `CarrollStructure.is_valid` never checks `rank(h)=n`
+*File:* `conformal_toolkit/carroll/carroll_structure.py` · *Status:* **Fixed (code)** — a zero or rank-deficient `h` previously passed validation despite the "rank-`n` degenerate" definition; a `rank(h) = dim−1` assertion was added.
+
+### `m6` — "On S²: H=K=1" docstring is radius-specific
+*File:* `conformal_features/discrete/q_curvature.py` · *Status:* **Fixed (docstring)** — on radius `R`, `H = 1/R`, `K = 1/R²`, so `H²−K = 0` for **all** `R`; only the antecedent `H=K=1` was `R=1`-specific.
+
+---
+
+## Style / Convention
+
+### `s1` — `conformal_killing.py`: wrong *and* dead symmetrization line
+*Status:* **Fixed (removed)** — `sym_nab_X = nab_X + nab_X.symmetrize()` is not
+`∇_aX_b + ∇_bX_a` and was never used; the returned residual is built from the
+correctly-formed `sym_tensor`, so the function's output was already correct. The
+misleading line and its "we need 2× it" comment were deleted.
+
+### `v1` — Discrete WKS omits per-energy normalization `C_e` and the standard `~7×` bandwidth
+*File:* `conformal_features/discrete/spectral.py` · *Status:* **Documented (convention)** — the canonical Aubry–Schmidt–Cremers WKS includes a per-energy `C_e = (Σⱼ exp(…))⁻¹` and `σ ≈ 7δ`; the code's unnormalized variant is a valid feature but deviates from the textbook definition.
+
+---
+
+## Teaching notes (the five lessons most worth internalizing)
+
+1. **A closed form must equal the number it annotates.** `2(n−1)!/((n/2−1)!)²`
+   is `12` at `n=4`, but it was labeling `Q₄ = 6`. Evaluate before you trust.
+   (`M1`)
+2. **The conformal Laplacian is *not* the bare Laplacian.**
+   `P₂ = Δ − (n−2)/(4(n−1)) R`. The curvature term *is* the conformal
+   covariance; an operator called "conformal" with no curvature term is a red
+   flag. (`M2`)
+3. **Conformal weight is a free checksum.** A zeroth-order term added to a
+   2nd-order Laplacian must carry weight `−2`, so it is **quadratic** in
+   curvature, never linear in `H`; and a Laplacian-type operator **lowers**
+   weight by 2 (bidegree difference `−2`). Weight bookkeeping catches structural
+   and sign errors before any computation. (`M4`)
+4. **Pointwise ≠ integral conformal invariance.** `H²−K` is the Willmore
+   integrand: only `∫(H²−K)dA` is conformally invariant, not its pointwise
+   value — and it is *not* the 4D GJMS `Q₄`. Know where your invariant lives.
+   (`M14`/`M15`/`M16`)
+5. **Discrete operators need the mass matrix.** The cotangent Laplacian `L` is a
+   *stiffness* matrix; the pointwise operator is `M⁻¹L` and the bi-Laplacian is
+   `(M⁻¹L)²`. "It converges to zero on the sphere" can be a false positive —
+   validate discretizations on a **non-constant** field like `f = x²`. (`M17`)
+
+---
+
+## How to contribute a correction
+
+We would rather be corrected than be wrong, and we designed this repository to
+make corrections cheap and welcome. If you find an error — in the math, the
+code, the docs, or in *this errata itself* — please:
+
+1. **Open an issue** titled `Errata: <one-line claim>` (or comment on an
+   existing `ERRATA <ID>`).
+2. **State the claim and the counter-evidence.** The most useful reports follow
+   the audit method above: give a *concrete geometry* (a sphere radius, a flat
+   patch, an explicit metric) on which the claim gives the wrong number, or a
+   *conformal-weight* argument showing the terms can't match. A failing
+   numerical check on a named test metric is gold.
+3. **Propose the fix** — the corrected formula, with the normalization
+   convention you're using, and a reference if you have one.
+4. If you can, **add a regression test** under `tests/` that pins the correct
+   value on the anchor geometry. Verified-on-an-anchor beats argued-in-prose.
+
+Open mathematical items where we explicitly want help: the complete weight-4
+hypersurface invariant basis and its exact count (`M5`/`M6`), the Fialkow/Weyl
+terms in `L₂` (`M7`), the full extrinsic `Q₄` (`M8`), the FG `g₄` Bach
+differential terms (`C1`), and the Graham–Hirachi `n=6` obstruction constant
+(`m3`).
+
+*Found something we missed? That's the system working. Send it.*
diff --git a/README.md b/README.md
index 5489f7c..bf7a5ec 100644
--- a/README.md
+++ b/README.md
@@ -9,25 +9,31 @@
 
 Two Python packages for computing conformal invariants — geometric quantities unchanged by local stretching — both symbolically (exact formulas via SageMath) and numerically on triangle meshes (GPU-ready via PyTorch). Implements tractor calculus, GJMS operators, Q-curvature, Blitz's conformal fundamental forms, Carroll geometry, and Fefferman-Graham holographic data.
 
+> **This library was independently re-audited for mathematical correctness in May 2026.** Every confirmed error — and exactly how we caught it — is documented in **[ERRATA.md](ERRATA.md)**. We publish the full record on purpose: a teaching tool earns trust by showing its work, mistakes included. See **["Verify it yourself"](#verify-it-yourself)** and **["Contributing corrections"](#contributing-corrections)**.
+
 ---
 
 ## Quick Start
 
-Compute Q-curvature exactly on the round 4-sphere, then verify the discrete approximation converges:
+Compute the Branson Q-curvature exactly on the round 4-sphere, then extract a discrete Willmore feature on a mesh:
 
 ```python
 # Symbolic: exact Q₄ on S⁴ via SageMath
 from conformal_toolkit import ConformalStructure
 cs = ConformalStructure(g_sphere4)
-cs.q_curvature(order=4)  # → 6  (exact, matches 2(n-1)! / ((n/2-1)!)² for n=4)
+cs.q_curvature(order=4)  # → 6   (exact; Branson's Q_n(Sⁿ) = (n-1)! = 3! = 6)
 
-# Discrete: Q₄ on an icosphere mesh via PyTorch
+# Discrete: a 4th-order surface feature on an icosphere mesh via PyTorch.
+# NOTE: this returns the Willmore integrand H² − K, NOT the 4D GJMS Q₄
+# (a 2-surface quantity cannot reproduce the 4-manifold value 6 — see ERRATA M15/M16).
 from conformal_features.discrete.q_curvature import discrete_q_curvature
 Q4 = discrete_q_curvature(vertices, faces, order=4)
-Q4.mean()  # → 5.94  (converges to 6 as mesh refines)
+Q4.mean()  # → 0   (H² − K vanishes on a round sphere of any radius; → 0 under refinement)
 ```
 
-One toolkit, two representations of the same invariant — symbolic formulas ground-truth the discrete features.
+The symbolic formulas ground-truth the geometry; the discrete features are
+mesh-domain analogues — *not* always the same number, and the docs now say which
+is which.
 
 ---
 
@@ -153,7 +159,7 @@ The 10 features per vertex:
 | 0 | Conformal factor | Conformal | From discrete Yamabe flow |
 | 1 | Willmore density | Conformal | H² (distance from minimality) |
 | 2 | Q₂ | Conformal | Discrete scalar curvature 2K |
-| 3 | Q₄ | Conformal | Higher-order curvature |
+| 3 | Willmore density H²−K | Conformal (integral) | 4th-order surface feature — the Willmore integrand, *not* the 4D GJMS Q₄ ([ERRATA M15](ERRATA.md)) |
 | 4 | Bach norm | Conformal | Bi-Laplacian proxy for non-conformal-flatness |
 | 5–6 | Cross-ratio stats | Möbius | Edge cross-ratio mean and variance |
 | 7 | Gaussian curvature | Isometric | Intrinsic curvature K |
@@ -240,7 +246,7 @@ conformal-toolkit/
 │   ├── discrete/               # Curvature, Q, Bach, Willmore, cross-ratios, Yamabe, spectral
 │   ├── features/               # mesh_conformal_features() pipeline
 │   └── benchmarks/             # ShapeNet, SHREC, FAUST evaluation (WIP)
-├── tests/                      # 157 tests across both packages
+├── tests/                      # 160 tests across both packages
 ├── examples/                   # 6 Jupyter notebooks
 └── paper.md                    # JOSS paper draft
 ```
@@ -273,7 +279,7 @@ conformal-toolkit/
 | `export` | `conformal_feature_vector(cs)` | Dict of all invariants at a point |
 | | `tensor_to_numpy(T)` | SageMath tensor → NumPy array |
 
-†P₆ computes the leading term (−Δ³) only; exact on conformally flat metrics. Obstruction at n=6 is a leading-order approximation.
+†P₆ computes the leading term (+Δ³) only; exact on conformally flat metrics. Obstruction at n=6 is a leading-order approximation (the Graham–Hirachi normalization constant is not applied — see [ERRATA m3](ERRATA.md)).
 
 ### conformal_features (PyTorch)
 
@@ -344,6 +350,66 @@ pytest tests/test_discrete/ tests/test_features/ -v
 
 ---
 
+## Verify it yourself
+
+Don't take our word for any formula — the whole point of a symbolic toolkit is
+that you can check it. Every correction in [ERRATA.md](ERRATA.md) was caught by
+evaluating a claim on a geometry where the answer is known in closed form. Here
+are the anchors we use; copy them into a Sage session and confirm:
+
+```python
+from sage.all import Manifold, sin
+from conformal_toolkit import ConformalStructure
+
+# --- Anchor 1: the round 4-sphere, where Branson's Q_n(Sⁿ) = (n-1)! ---
+# On Sⁿ (sectional curvature 1):  Ric = (n-1)g,  R = n(n-1),
+#   Schouten P = ½g,  J = tr P = n/2,  so  Q₄ = -ΔJ - 2|P|² + (n/2)J² = 6.
+cs = ConformalStructure(g_sphere4)
+assert cs.q_curvature(order=4) == 6          # = (4-1)! = 3!  (ERRATA M1)
+
+# --- Anchor 2: the conformal Laplacian carries a curvature term ---
+# P₂ f = Δf - (n-2)/(4(n-1)) R f.  The R-term is NOT optional for n > 2.
+# On S⁴ its coefficient is n(n-2)/4 = 2, never zero.  (ERRATA M2)
+
+# --- Anchor 3: Fefferman-Graham on the hyperbolic filling of Sⁿ ---
+# g_ρ = (1 - ρ²/4)² g₀  ⟹  g₂ = -½ g₀  and  g₄ = 1/16 g₀ exactly.  (ERRATA C1/M11)
+
+# --- Anchor 4 (discrete, PyTorch only): validate on a NON-constant field ---
+# The cotangent Laplacian L is a *stiffness* matrix; for f = x² on a flat mesh,
+# (L f)_i = -2·A_i  while  (M⁻¹ L f)_i = -2  recovers the pointwise Laplacian.
+# "It vanishes on a sphere" is a false positive — constants are annihilated by
+# any linear operator.  (ERRATA M17)
+```
+
+The method generalizes: **reduce a tensor claim to a scalar on a known geometry,
+and check conformal weights as a free checksum.** That single discipline caught
+most of the errata.
+
+---
+
+## Contributing corrections
+
+We would rather be corrected than be wrong, and this repository is built to make
+that easy. **Finding an error we missed is the system working — please send it.**
+
+1. Open an issue titled `Errata: <one-line claim>`.
+2. Give the counter-evidence the way we give ours: a *concrete geometry* (a
+   sphere radius, a flat patch, an explicit metric) on which the claim returns
+   the wrong number, or a *conformal-weight* argument that the terms can't match.
+   A failing check on a named anchor metric is the gold standard.
+3. Propose the corrected formula, stating your normalization convention (Branson
+   vs. analyst signs differ — half of conformal geometry's "errors" are
+   convention clashes), with a reference if you have one.
+4. If you can, add a regression test under `tests/` pinning the right value on
+   the anchor. *Verified-on-an-anchor beats argued-in-prose.*
+
+Open problems where we explicitly want help are listed at the end of
+[ERRATA.md](ERRATA.md) — the complete weight-4 hypersurface invariant basis
+(`M5`), the Fialkow/Weyl terms in `L₂` (`M7`), the full extrinsic `Q₄` (`M8`),
+and the FG `g₄` Bach differential terms (`C1`).
+
+---
+
 ## Citation
 
 ```bibtex
diff --git a/conformal_features/discrete/bach.py b/conformal_features/discrete/bach.py
index cc8e62f..61f2596 100644
--- a/conformal_features/discrete/bach.py
+++ b/conformal_features/discrete/bach.py
@@ -1,31 +1,49 @@
 """Discrete approximation of Bach tensor norm.
 
-|B| ~ |L^2(K)| as a 4th-order curvature feature.
-This is not the true Bach tensor but a computationally tractable proxy
-that captures 4th-order curvature variation.
-
-The bi-Laplacian L^2 K (cotangent Laplacian applied twice to the discrete
-Gaussian curvature) acts as a 4th-order differential operator. Since
-discrete_gaussian_curvature already normalizes K = angle_defect / area, the
-raw bi-Laplacian (no additional area weighting) gives values that are small on
-a smooth, nearly-uniform sphere and converge to zero under mesh refinement.
+This is NOT the true Bach tensor (the 4D Fefferman-Graham obstruction) but a
+computationally tractable scalar PROXY that captures 4th-order curvature
+variation on a surface, used as one channel of the per-vertex feature vector.
+
+Discretization honesty note (see ERRATA M17). The cotangent matrix L is the
+FEM *stiffness* matrix, i.e. an INTEGRATED operator: (L f)_i ~ A_i (Delta f)_i,
+not (Delta f)_i. So the quantity computed here, |L L K|, is a *stiffness-
+weighted* (area-integrated) bi-Laplacian, NOT the pointwise bi-Laplacian. The
+true pointwise operator is
+
+    Delta^2 = (M^{-1} L)(M^{-1} L) = M^{-1} L M^{-1} L,   M = diag(vertex areas).
+
+We deliberately keep the integrated form |L L K| as the feature, because it is
+numerically well-behaved across mesh scales; but note two caveats:
+  * its *small absolute magnitude on a sphere is partly a scaling artifact*
+    (the A_i^2 weighting suppresses it), not proof of correctness -- the
+    earlier docstring's claim "K is already area-normalized, so no extra
+    weighting is needed" was the wrong reason for the right-ish behavior;
+  * the mass-lumped pointwise Delta^2 above is the mathematically correct
+    operator but amplifies coarse-mesh curvature noise (e.g. the 12 pentagonal
+    defects of an icosphere), so it is offered as `pointwise=True` rather than
+    the default. Always validate a discretization on a NON-constant field such
+    as f = x^2, where (M^{-1} L)(x^2) = -2 recovers the pointwise Laplacian
+    while (L)(x^2) = -2 A_i does not.
 """
 from __future__ import annotations
 import torch
 from conformal_features.discrete.curvature import discrete_gaussian_curvature
-from conformal_features.discrete.mesh_utils import cotangent_laplacian
+from conformal_features.discrete.mesh_utils import cotangent_laplacian, vertex_areas
 
 
-def discrete_bach_norm(vertices: torch.Tensor, faces: torch.Tensor) -> torch.Tensor:
-    """Per-vertex discrete Bach tensor norm approximation.
-
-    Computes |L @ L @ K| where K is the discrete Gaussian curvature and L is
-    the cotangent Laplacian. This bi-Laplacian captures 4th-order curvature
-    variation and is small on a sphere (where K is nearly constant).
+def discrete_bach_norm(
+    vertices: torch.Tensor, faces: torch.Tensor, pointwise: bool = False
+) -> torch.Tensor:
+    """Per-vertex discrete Bach tensor norm proxy.
 
     Args:
         vertices: (V, 3) Tensor
         faces: (F, 3) LongTensor
+        pointwise: if False (default), return the stiffness-weighted
+            (area-integrated) bi-Laplacian |L L K| -- scale-stable, used as the
+            feature channel. If True, return the mathematically-correct
+            pointwise bi-Laplacian |M^{-1} L M^{-1} L K| (noisier on coarse
+            meshes). See ERRATA M17.
 
     Returns:
         bach_norm: (V,) Tensor
@@ -33,10 +51,13 @@ def discrete_bach_norm(vertices: torch.Tensor, faces: torch.Tensor) -> torch.Ten
     K = discrete_gaussian_curvature(vertices, faces)
     L = cotangent_laplacian(vertices, faces)
 
-    # Apply the cotangent Laplacian twice. K is already area-normalized
-    # (angle_defect / area), so no additional area weighting is needed —
-    # the raw bi-Laplacian L@L@K converges to zero on a sphere where K is constant.
-    LK = torch.sparse.mm(L, K.unsqueeze(1)).squeeze(1)
-    LLK = torch.sparse.mm(L, LK.unsqueeze(1)).squeeze(1)
+    def stiffness(field: torch.Tensor) -> torch.Tensor:
+        return torch.sparse.mm(L, field.unsqueeze(1)).squeeze(1)
+
+    if pointwise:
+        inv_area = 1.0 / vertex_areas(vertices, faces).clamp_min(1e-12)
+        lap = lambda f: inv_area * stiffness(f)   # pointwise Delta = M^{-1} L
+        return lap(lap(K)).abs()
 
-    return LLK.abs()
+    # Integrated proxy |L L K| (the default feature channel).
+    return stiffness(stiffness(K)).abs()
diff --git a/conformal_features/discrete/q_curvature.py b/conformal_features/discrete/q_curvature.py
index a55c0ab..2aa76b9 100644
--- a/conformal_features/discrete/q_curvature.py
+++ b/conformal_features/discrete/q_curvature.py
@@ -9,13 +9,26 @@
 
 
 def discrete_q_curvature(vertices: torch.Tensor, faces: torch.Tensor, order: int = 2) -> torch.Tensor:
-    """Discrete Q-curvature of the given order.
-
-    Q_2 = R = 2K (scalar curvature = 2 * Gaussian curvature for surfaces)
-    Q_4 approximation: H^2 - K  (conformally invariant Willmore integrand).
-        On S^2: H = K = 1, so Q_4 = 0.
-        More precisely this approximates -Delta(J) - 2|P|^2 + J^2 in 2D via
-        the identity Q_4 = 2K^2 - 2KH^2 (rescaled), which equals 0 on S^2.
+    """Discrete curvature surface features (orders 2 and 4).
+
+    order=2:  Q_2 = R = 2K  (scalar curvature = 2 * Gaussian curvature
+              for a surface).
+
+    order=4:  returns the *Willmore integrand* H^2 - K, NOT the 4D GJMS /
+              Branson Q_4. Two honesty caveats (see ERRATA M14/M15/M16):
+        * Naming. The intrinsic 4th-order GJMS Q_4 is a 4-manifold object
+          with Q_4 = (n-1)! = 6 on the round S^4. A 2-surface quantity such
+          as H^2 - K cannot reproduce that value; it is a Willmore-type
+          surface feature, not the GJMS Q_4.
+        * Invariance. H^2 - K is conformally invariant only *under the
+          integral*: ∫(H^2 - K) dA is controlled by Gauss-Bonnet
+          (∫K dA = 2*pi*chi) plus Moebius-invariance of ∫H^2 dA. It is NOT
+          pointwise conformally invariant.
+        * Convergence. On a round S^2 of ANY radius R, H^2 = 1/R^2 = K, so
+          H^2 - K = 0; thus this feature -> 0 under mesh refinement on a
+          sphere (it does NOT converge to 6 -- an earlier README claim).
+          (The earlier "identity" Q_4 = 2K^2 - 2KH^2 was false: that equals
+          -2K (H^2 - K), not a rescaling of H^2 - K.)
 
     Args:
         vertices: (V, 3) Tensor
@@ -32,7 +45,7 @@ def discrete_q_curvature(vertices: torch.Tensor, faces: torch.Tensor, order: int
 
     if order == 4:
         H = discrete_mean_curvature(vertices, faces)
-        # Q_4 ~ H^2 - K: conformally invariant, vanishes on S^2 (H=K=1)
+        # Willmore integrand H^2 - K; -> 0 on a round sphere of any radius.
         return H ** 2 - K
 
     raise ValueError(f"Order {order} not supported (only 2 and 4)")
diff --git a/conformal_toolkit/carroll/bms.py b/conformal_toolkit/carroll/bms.py
index 6002ebd..8a76bdd 100644
--- a/conformal_toolkit/carroll/bms.py
+++ b/conformal_toolkit/carroll/bms.py
@@ -58,6 +58,17 @@ def is_bms_symmetry(carroll_struct, vector_field):
     A vector field ξ is a **strict Carroll symmetry** if £_ξ h = 0, and a
     **conformal Carroll symmetry** if £_ξ h = λ h for some scalar λ.
 
+    INCOMPLETE PREDICATE -- necessary but NOT sufficient (see ERRATA M10).
+    A Carroll structure is the PAIR (v, h), so a symmetry must preserve BOTH:
+        strict:    £_ξ h = 0      and  £_ξ v = 0
+        conformal: £_ξ h = 2λ h   and  £_ξ v = -λ v   (Duval-Gibbons-Horvathy)
+    Because h is degenerate (h(v, ·) = 0), the £_ξ h condition places NO
+    constraint on the v-direction, so this routine -- which only tests
+    £_ξ h -- can wrongly accept a ξ that moves v out of ker(h). A correct
+    predicate must ALSO verify £_ξ v ∝ v (strict: = 0; conformal: = -λ v).
+    (The supertranslation generator ξ = f·v built above does satisfy
+    £_ξ v = 0, so it passes; the gap bites for general user input.)
+
     The Lie derivative is computed component-wise:
 
         (£_ξ h)_{ij} = ξ^k ∂_k h_{ij} + h_{kj} ∂_i ξ^k + h_{ik} ∂_j ξ^k
diff --git a/conformal_toolkit/carroll/carroll_connection.py b/conformal_toolkit/carroll/carroll_connection.py
index 1b8587c..b9162b4 100644
--- a/conformal_toolkit/carroll/carroll_connection.py
+++ b/conformal_toolkit/carroll/carroll_connection.py
@@ -1,15 +1,23 @@
 """Carroll connection: Christoffel-like symbols adapted to a Carroll structure.
 
-For a Carroll manifold with adapted coordinates (t, x^i) where v = ∂_t, the
-Carroll connection has the block structure:
+For a Carroll manifold with adapted coordinates (t, x^i) where v = ∂_t, this
+module uses the block structure:
 
     Γ^t_{μν} = 0   (all)
     Γ^i_{tt} = 0
     Γ^i_{tj} = 0
     Γ^i_{jk} = spatial Christoffel symbols from h_{jk}
 
-The spatial Christoffel symbols are the Levi-Civita symbols of the restriction
-of h to the spatial leaves (where h is non-degenerate).
+SCOPE CAVEAT (see ERRATA M9). Carrollian connections are NOT unique (h is
+degenerate), and the choice Γ^i_{tj} = 0 is correct ONLY when the Carroll
+"electric field" E_{ij} = (1/2)(£_v h)_{ij} vanishes -- i.e. when the spatial
+metric is time-independent (∂_t h_{ij} = 0). When ∂_t h_{ij} ≠ 0, metric
+compatibility ∇h = 0 forces Γ^k_{t(i} h_{j)k} = (1/2) ∂_t h_{ij} ≠ 0, so the
+correct symmetric part is Γ^i_{(tj)} = -(1/2) h^{ik}(£_v h)_{kj}. This module
+otherwise computes a generically NONZERO carroll_electric_field, so presenting
+Γ^i_{tj} = 0 as THE Carroll connection is internally inconsistent in the
+time-dependent case. Treat the symbols below as the preserved-h
+(£_v h = 0) representative.
 """
 
 from sage.all import SR
diff --git a/conformal_toolkit/carroll/carroll_structure.py b/conformal_toolkit/carroll/carroll_structure.py
index 25e59f0..66c366d 100644
--- a/conformal_toolkit/carroll/carroll_structure.py
+++ b/conformal_toolkit/carroll/carroll_structure.py
@@ -58,12 +58,21 @@ def is_valid(self):
         """Check the Carroll compatibility condition h(v, ·) = 0.
 
         Contracts h with v on its first slot and verifies that all components
-        of the resulting 1-form vanish symbolically.
+        of the resulting 1-form vanish symbolically. Also checks that the
+        spatial metric h has rank n (a 1-dimensional kernel spanned by v),
+        as required by the definition of a Carroll structure.
+
+        INCOMPLETENESS NOTE (ERRATA m5): the original check verified only
+        h(v, ·) = 0 and NOT rank(h) = n, so a fully-zero h or any
+        rank-deficient h (kernel dimension > 1) passed validation despite
+        the "rank-n degenerate" definition. The rank check below closes that
+        gap (downstream spatial_christoffel also implicitly needs the spatial
+        block invertible).
 
         Returns
         -------
         bool
-            True if v lies in the kernel of h, False otherwise.
+            True if v lies in the kernel of h and rank(h) = n, else False.
         """
         h = self._h
         v = self._v
@@ -85,6 +94,17 @@ def is_valid(self):
             except Exception:
                 if not bool(val == 0):
                     return False
+
+        # Rank check: h must be degenerate of rank exactly n = dim - 1
+        # (kernel = span(v), one-dimensional). See ERRATA m5.
+        from sage.all import matrix, SR
+        h_mat = matrix(SR, self._dim, self._dim)
+        for a in range(self._dim):
+            for b in range(self._dim):
+                comp = h[frame, a, b]
+                h_mat[a, b] = comp.expr() if hasattr(comp, 'expr') else SR(comp)
+        if h_mat.rank() != self._dim - 1:
+            return False
         return True
 
     # ------------------------------------------------------------------
diff --git a/conformal_toolkit/core/conformal_killing.py b/conformal_toolkit/core/conformal_killing.py
index ab8a727..3b09631 100644
--- a/conformal_toolkit/core/conformal_killing.py
+++ b/conformal_toolkit/core/conformal_killing.py
@@ -69,13 +69,10 @@ def conformal_killing_equation(g, X):
     X_flat = X.down(g)           # (0,1): X_b
     nab_X = nabla(X_flat)        # (0,2): ∇_a X_b
 
-    # Symmetrize: ∇_a X_b + ∇_b X_a
-    sym_nab_X = nab_X + nab_X.symmetrize()  # symmetrize() gives (∇_a X_b + ∇_b X_a)/2
-    # Correct: sym_nab_X = nab_X + nab_X.swap_adjacent_indices(0, 1) would swap,
-    # but SageManifolds provides .symmetrize() which averages, so we need 2* it.
-    # Instead: manually build the symmetrized version.
-    # nab_X[a,b] = ∇_a X_b;  we want nab_X[a,b] + nab_X[b,a].
-
+    # Build the symmetrized derivative ∇_a X_b + ∇_b X_a explicitly below.
+    # (An earlier line `nab_X + nab_X.symmetrize()` was both mathematically
+    #  wrong -- it equals (3/2)∇_(a X_b) + (1/2)∇_[a X_b], not ∇_a X_b + ∇_b X_a
+    #  -- and dead code, never used; removed. See ERRATA S1.)
     M = g.domain()
     frame = list(M.frames())[0]
 
diff --git a/conformal_toolkit/core/conformal_structure.py b/conformal_toolkit/core/conformal_structure.py
index 8cd59c0..85de43f 100644
--- a/conformal_toolkit/core/conformal_structure.py
+++ b/conformal_toolkit/core/conformal_structure.py
@@ -82,9 +82,21 @@ def ricci_scalar(self):
         return self._g.ricci_scalar()
 
     def bach(self):
-        """Bach tensor B_ab = nabla^c nabla^d W_{acbd} + P^{cd} W_{acbd} (cached).
+        """Bach tensor (cached).
 
-        Vanishes in dimension 4 iff the metric is locally conformally Einstein.
+        Computed (in bach.py) via the Cotton form, valid for all n:
+
+            B_ab = nabla^c C_{cab} + P^{cd} W_{cadb}
+
+        Note: the Weyl-divergence form B_ab = (1/(n-3)) nabla^c nabla^d
+        W_{acbd} + P^{cd} W_{acbd} carries a 1/(n-3) factor, since the second
+        Bianchi identity gives nabla^c W_{cabd} = (n-3) C_{abd}; the two
+        forms agree only at n=4 (where 1/(n-3)=1). An earlier docstring
+        dropped that factor -- see ERRATA m2.
+
+        The Bach tensor is conformally invariant (and the obstruction tensor)
+        only in dimension 4, where it vanishes iff the metric is locally
+        conformally Einstein.
         """
         if self._bach is None:
             from conformal_toolkit.core.bach import compute_bach
diff --git a/conformal_toolkit/core/gjms.py b/conformal_toolkit/core/gjms.py
index 989a95a..e7fc1cf 100644
--- a/conformal_toolkit/core/gjms.py
+++ b/conformal_toolkit/core/gjms.py
@@ -1,16 +1,32 @@
 """GJMS operators (Graham-Jenne-Mason-Sparling).
 
-P_2 = Delta (Laplacian)
+P_2 = conformal (Yamabe) Laplacian = Delta - ((n-2)/(4(n-1))) R
 P_4 = Paneitz operator = Delta^2 + div(V . d) + ((n-4)/2) Q_4
 
 where V_ab = (n-2) J g_ab - 4 P_ab.
+
+Sign convention: Delta = nabla^a nabla_a (geometer / non-negative-spectrum
+on a compact manifold is for -Delta). The GJMS principal part is +Delta^k.
 """
 from conformal_toolkit.core.helpers import laplacian, divergence
 
 
 def laplacian_operator(cs, f):
-    """P_2(f) = Delta(f) = nabla^a nabla_a f."""
-    return laplacian(cs.connection(), cs.metric, f)
+    """P_2(f): the *conformal* (Yamabe) Laplacian, not the bare Laplacian.
+
+        P_2(f) = Delta(f) - ((n-2)/(4(n-1))) R f
+
+    The scalar-curvature term is exactly what makes P_2 conformally
+    covariant (weight (2-n)/2 -> (-2-n)/2) for n > 2; it vanishes only at
+    n = 2. Omitting it (the original bug) breaks conformal covariance for
+    every n > 2 -- see ERRATA M2.
+    """
+    n = cs.dimension
+    R = cs.ricci_scalar()
+    delta_f = laplacian(cs.connection(), cs.metric, f)
+    if n == 2:
+        return delta_f
+    return delta_f - ((n - 2) / (4 * (n - 1))) * R * f
 
 
 def paneitz_operator(cs, f):
@@ -40,10 +56,16 @@ def paneitz_operator(cs, f):
 def p6_operator(cs, f):
     """P_6(f): sixth-order GJMS operator (leading-order approximation).
 
-    Computes (-1)^3 Delta^3(f) = -Delta^3(f), which equals P_6(f) for
-    conformally flat metrics.  On general curved backgrounds this captures
-    the principal symbol but omits lower-order curvature coupling terms
-    (Schouten, Cotton, and Weyl contributions).
+    Computes +Delta^3(f) (the GJMS principal part is +Delta^k in the
+    geometer sign convention Delta = nabla^a nabla_a used throughout this
+    package), which equals P_6(f) for conformally flat metrics.  On general
+    curved backgrounds this captures the principal symbol but omits the
+    lower-order curvature coupling terms (Schouten, Cotton, and Weyl
+    contributions).
+
+    Note: an earlier docstring claimed this returns -Delta^3; that was a
+    sign error from mixing the analyst convention (Delta = -nabla^a nabla_a)
+    with the geometer Laplacian the code actually applies -- see ERRATA M3.
 
     Parameters
     ----------
diff --git a/conformal_toolkit/core/obstruction.py b/conformal_toolkit/core/obstruction.py
index 46e7bf3..fe1c6ef 100644
--- a/conformal_toolkit/core/obstruction.py
+++ b/conformal_toolkit/core/obstruction.py
@@ -11,9 +11,17 @@ def _obstruction_6(cs):
     """Obstruction tensor in dimension 6 (leading-order approximation).
 
     Computes Delta(Bach)_{ab} = g^{cd} nabla_c nabla_d B_{ab} as the
-    dominant contribution to the 6-dimensional obstruction tensor.
-    The full O_6 involves additional cubic-in-curvature terms; on
-    conformally flat metrics all such terms vanish and this is exact.
+    leading (principal-part) contribution to the 6-dimensional obstruction
+    tensor.
+
+    Caveats (see ERRATA m3):
+    * The overall Graham-Hirachi normalization constant c in
+      O_6 = c * Delta(Bach) + (curvature-squared corrections) is NOT applied
+      here -- this routine returns the unnormalized Delta(Bach).
+    * The full O_6 carries additional quadratic/cubic-in-curvature terms,
+      omitted here. On conformally flat metrics every term (including
+      Delta(Bach), since Bach = 0) vanishes, so "exact on conformally flat"
+      reduces to the vacuous 0 = 0 and does NOT validate the normalization.
 
     Parameters
     ----------
diff --git a/conformal_toolkit/core/q_curvature.py b/conformal_toolkit/core/q_curvature.py
index 6cc4367..f3a3f4a 100644
--- a/conformal_toolkit/core/q_curvature.py
+++ b/conformal_toolkit/core/q_curvature.py
@@ -2,6 +2,16 @@
 
 Q_2 = R (scalar curvature)
 Q_4 = -Delta(J) - 2|P|^2 + (n/2) J^2
+
+NORMALIZATION CAVEAT (ERRATA m1). This module returns Q_2 = R, whereas Q_4
+follows the Branson normalization (Q_4 = (n-1)! = 6 on the round S^4). Those
+two conventions are inconsistent ACROSS ORDERS: the Branson-normalized Q_2
+is R/2 (= Gaussian curvature K = 1 = (2-1)! on S^2), not R (= 2). Both
+Q_2 = R and Q_2 = R/2 appear in the literature, so this is a convention
+choice rather than an arithmetic error -- but a user mixing Q_2 and Q_4
+should be aware they differ by a factor of 2 in normalization. Q_4 here is
+the general-n Branson Q-curvature (it carries (n/2)J^2, not 2J^2; the two
+coincide only at n=4).
 """
 from conformal_toolkit.core.helpers import laplacian, tensor_norm_squared
 
diff --git a/conformal_toolkit/hypersurface/conformal_fundamental_form.py b/conformal_toolkit/hypersurface/conformal_fundamental_form.py
index 4d3b204..136d6e2 100644
--- a/conformal_toolkit/hypersurface/conformal_fundamental_form.py
+++ b/conformal_toolkit/hypersurface/conformal_fundamental_form.py
@@ -64,7 +64,7 @@ def conformal_fundamental_form_L1(h, L):
 def conformal_fundamental_form_L2(h, L, nabla_n_L1=None):
     """Second conformal fundamental form.
 
-    Full formula:
+    Formula (CONFORMALLY-FLAT-AMBIENT reduction):
         L_2_{ab} = (nabla_n L_1)_{ab} + (L_1)_a^c (L_1)_{cb}
                    - (1/n)|L_1|^2 h_{ab}
 
@@ -73,6 +73,14 @@ def conformal_fundamental_form_L2(h, L, nabla_n_L1=None):
 
         L_2^{alg}_{ab} = (L_1)_a^c (L_1)_{cb} - (1/n)|L_1|^2 h_{ab}
 
+    WARNING -- INCOMPLETE on a general ambient (see ERRATA M7). The genuine
+    conformal fundamental form at this order is the trace-free FIALKOW
+    tensor (Blitz-Gover-Waldron, arXiv:2107.10381), which necessarily
+    contains ambient curvature (Weyl/Cotton, e.g. W_{n(ab)n}) terms. The
+    expression above carries NO Weyl/Cotton term, so it is conformally
+    invariant only when the ambient is (locally) conformally flat; on a
+    generic ambient it is the conformally-flat reduction, not L_2.
+
     Args:
         h: induced metric on the hypersurface.
         L: second fundamental form (0,2) on the hypersurface.
diff --git a/conformal_toolkit/hypersurface/extrinsic_gjms.py b/conformal_toolkit/hypersurface/extrinsic_gjms.py
index 851fc7e..0a5b361 100644
--- a/conformal_toolkit/hypersurface/extrinsic_gjms.py
+++ b/conformal_toolkit/hypersurface/extrinsic_gjms.py
@@ -18,7 +18,17 @@ def extrinsic_gjms_P2(h, L, f):
 
     P_2(f) = Delta_h(f) + (n/2 - 1) H f
 
-    This is conformally covariant of bi-degree (n/2 - 1, n/2 + 1).
+    WARNING -- NOT conformally covariant as written (see ERRATA M4). Two
+    issues: (1) the zeroth-order term must carry conformal weight -2 (Delta
+    lowers weight by 2), so it must be QUADRATIC in the mean curvature
+    (H^2-type) and/or the intrinsic Schouten trace J_bar -- a term LINEAR in
+    H (weight -1) is weight-inhomogeneous. A correct leading extrinsic
+    conformal Laplacian has the Yamabe-type shape
+    Delta_h f - ((n-2)/2)(J_bar + extrinsic H^2-term) f. (2) The stated
+    bi-degree (n/2 - 1, n/2 + 1) RAISES weight by 2 (target - source = +2),
+    but a second-order Laplacian-type operator LOWERS weight by 2; the
+    standard k=1 conformal-Laplacian bi-degree is (1 - n/2, -1 - n/2)
+    (target - source = -2).
 
     Args:
         h: induced metric on the hypersurface.
diff --git a/conformal_toolkit/hypersurface/extrinsic_q.py b/conformal_toolkit/hypersurface/extrinsic_q.py
index 094f332..df74149 100644
--- a/conformal_toolkit/hypersurface/extrinsic_q.py
+++ b/conformal_toolkit/hypersurface/extrinsic_q.py
@@ -31,8 +31,23 @@ def extrinsic_q2(h, L):
 def extrinsic_q4(h, L, nabla_n_L1=None):
     """Extrinsic Q-curvature of order 4.
 
-    Simplified formula (ignoring ambient curvature corrections):
-        q_4 = -Delta_h(H) + H |L_1|^2 + (n/2 - 1) H^3
+    WARNING -- INCOMPLETE (see ERRATA M8). This returns only the purely
+    extrinsic, conformally-flat-ambient piece
+
+        q_4 = -Delta_h(H) + H |L_1|^2 + (n/2 - 1) H^3,
+
+    which OMITS the intrinsic fourth-order term. The genuine extrinsic Q_4 is
+
+        q_4 = Q_4^Sigma + (extrinsic couplings),
+        Q_4^Sigma = -Delta_bar J_bar - 2|P_bar|^2 + (n/2) J_bar^2,
+
+    plus the leading second-derivative coupling ~ (1/2) L_1 . Delta_bar L_1.
+    Failing anchor: on the round S^4 (umbilic: L_1 = 0, H = 1) this routine
+    returns 0 + 0 + (n/2-1) = 1, whereas the correct value must reduce to the
+    intrinsic Branson Q_4 = (n-1)! = 6. Use for the umbilic-deviation
+    (Willmore-type) content only, not as the true extrinsic Q_4.
+    (Also: the module header writes +Delta_h H but the implementation uses
+    -Delta_h H; the implementation's sign is kept.)
 
     where Delta_h is the Laplace-Beltrami operator on (Sigma, h).
 
diff --git a/conformal_toolkit/hypersurface/invariant_enumeration.py b/conformal_toolkit/hypersurface/invariant_enumeration.py
index 49cd437..f49a10a 100644
--- a/conformal_toolkit/hypersurface/invariant_enumeration.py
+++ b/conformal_toolkit/hypersurface/invariant_enumeration.py
@@ -6,6 +6,16 @@
 At each weight, the independent scalar invariants are formed from
 contractions of the conformal fundamental forms L_k, the ambient Weyl
 tensor restricted to Sigma, and the intrinsic Weyl tensor of Sigma.
+
+KNOWN INCOMPLETENESS at weight 4 (see ERRATA M5). The weight-4 catalogue
+below is a PARTIAL (purely algebraic, leading) basis. It under-counts:
+  * it omits the independent algebraic quartic (|L_1|^2)^2, which is
+    distinct from tr(L_1^4) for n >= 4 (their difference is 4abc(a+b+c)
+    in eigenvalues, generically nonzero); and
+  * it omits the tangential-derivative invariants (e.g. |nabla_bar L_1|^2 /
+    L_1 . Delta_bar L_1, |div L_1|^2) and the curvature coupling
+    J_bar |L_1|^2 that a complete pointwise weight-4 classification carries.
+So count_invariants(order=4, ...) is a LOWER BOUND, not the exact count.
 """
 
 
@@ -71,7 +81,15 @@ def _filter_by_dimension(invariants, ambient_dim):
         # Weyl-based invariants vanish identically when ambient_dim <= 3
         if 'W_' in inv['name'] and ambient_dim <= 3:
             continue
-        # In ambient 4D the Weyl is self-dual, |W_{nabc}|^2 is dependent
+        # PROVISIONAL, UNVERIFIED (ERRATA M6): |W_{nabc}|^2 is dropped at
+        # ambient_dim == 4. The ORIGINAL justification ("the Weyl is
+        # self-dual") was mathematically FALSE -- a generic 4-manifold has
+        # both W^+ and W^- nonzero and independent, and even W^- = 0 would
+        # not make the scalar |W_{nabc}|^2 expressible via the other listed
+        # invariants. We have NOT established a genuine algebraic dependence
+        # here, so this drop (and hence count_invariants(4, 4)) is an
+        # unverified placeholder, not a proven count. Contributors: please
+        # supply a real identity or remove the drop. See ERRATA M6.
         if inv['name'] == '|W_{nabc}|^2' and ambient_dim == 4:
             continue
         filtered.append(inv)
diff --git a/conformal_toolkit/poincare_einstein/dirichlet_neumann.py b/conformal_toolkit/poincare_einstein/dirichlet_neumann.py
index b7d07a6..833d470 100644
--- a/conformal_toolkit/poincare_einstein/dirichlet_neumann.py
+++ b/conformal_toolkit/poincare_einstein/dirichlet_neumann.py
@@ -26,10 +26,17 @@ def holographic_stress_tensor(g0, n):
     T = n * g_n + (lower order trace terms).
     At leading non-trivial order: T = n * g_2 = -n * P(g_0).
 
-    For n=3 (AdS4/CFT3) the trace anomaly is absent and:
-        T_{ab} = -3 * P_{ab} + (3/2) * J * (g_0)_{ab}
+    For n=3 (AdS4/CFT3) the trace anomaly is absent, so any *local*
+    placeholder must be trace-free. The trace-free combination is
 
-    where J = trace(P) = R/(2*(n-1)).
+        T_{ab} = -3 * P_{ab} + J * (g_0)_{ab}      (tr T = -3J + J*3 = 0)
+
+    where J = trace(P) = R/(2*(n-1)). The earlier coefficient (3/2)*J gave
+    tr T = (3/2)J != 0, contradicting the stated absence of an anomaly (see
+    ERRATA M13). NOTE: the genuine n=3 holographic stress tensor is
+    T_{ij} = 3 g_{(3)ij}, which is undetermined *non-local* VEV data (free
+    boundary condition), not any local function of P -- this branch returns
+    only a trace-consistent local placeholder.
 
     For general n the holographic stress tensor at leading order is:
         T_{ab} = -n * P_{ab}
@@ -49,7 +56,7 @@ def holographic_stress_tensor(g0, n):
     J = schouten_trace(g0)
 
     if n == 3:
-        T = -3 * P + (3 / 2) * J * g0
+        T = -3 * P + J * g0   # trace-free local placeholder; see ERRATA M13
     else:
         T = -n * P
 
diff --git a/conformal_toolkit/poincare_einstein/fefferman_graham.py b/conformal_toolkit/poincare_einstein/fefferman_graham.py
index 60c656c..54d742d 100644
--- a/conformal_toolkit/poincare_einstein/fefferman_graham.py
+++ b/conformal_toolkit/poincare_einstein/fefferman_graham.py
@@ -23,13 +23,30 @@ def fg_coefficient_g2(g0):
 def fg_coefficient_g4(g0):
     """Compute g_4 in the FG expansion from boundary metric g_0.
 
-    For n != 4:
-        (g_4)_{ab} = (1/(n-4)) * [(g_2)_a^c (g_2)_{cb}
-                     - (1/(4*(n-1))) tr(g_2^2) (g_0)_{ab}]
+    LEADING-ORDER / LOCALLY-CONFORMALLY-FLAT result. The full g_4 is
 
-    For n == 4 the metric expansion develops a log term and g_4 is only
-    determined modulo the trace-free part.  We return the pure-trace piece
-    that is algebraically determined: (1/8) tr(g_2^2) g_0.
+        (g_4)_{ab} = (1/4) (g_2)_a^c (g_2)_{cb}
+                     + (1/(n-4)) * [Bach / Delta-P differential terms],
+
+    where the *algebraic* g_2^2 piece carries the coefficient 1/4 and the
+    1/(n-4) pole sits ONLY on the differential (Bach-tensor) terms that
+    obstruct at n=4 (de Haro-Skenderis-Solodukhin / Graham). This routine
+    implements the algebraic piece, which is EXACT on locally conformally
+    flat boundaries (where the Bach terms vanish) and is the correct leading
+    behaviour otherwise. The differential terms are not yet implemented.
+
+    For n == 4 the trace-free part of g_4 is undetermined (a log term
+    appears); only the pure-trace piece is algebraically fixed:
+
+        g_4 = (1/16) tr(g_2^2) g_0    =>   tr g_4 = (1/4) tr(g_2^2) = (1/4)|P|^2.
+
+    Verification anchor: on the hyperbolic filling of the round S^n
+    (g_rho = (1 - rho^2/4)^2 g_0) one has g_2 = -(1/2)g_0 and the EXACT
+    g_4 = (1/16) g_0 = (1/4)(g_2)^2 -- reproduced by this routine.
+
+    The earlier version placed a spurious 1/(n-4) on the algebraic term
+    (giving 7/80 g_0 instead of 1/16 g_0 at n=6) and used 1/8 instead of
+    1/16 at n=4 (doubling the trace) -- see ERRATA C1 and M11.
     """
     n = g0.domain().dim()
     M = g0.domain()
@@ -39,22 +56,17 @@ def fg_coefficient_g4(g0):
     g2_up = g2.up(g0, 0)           # (1,1) tensor
     g2_sq = g2.contract(0, g2_up, 0)  # contract first index pair -> (0,2)
 
-    # Trace: tr(g2^2) = (g2^a_b)(g2^b_a) = g2_sq.trace(0,1) but g2_sq is (0,2),
-    # so we need to raise and trace.
-    g2_sq_up = g2_sq.up(g0, 0)
-    tr_g2_sq = g2_sq_up.trace(0, 1)
-
-    frame = list(M.frames())[0]
-
     if n == 4:
-        # Only the algebraically-fixed trace piece: (1/8) tr(g2^2) g0.
-        g4 = (1 / 8) * tr_g2_sq * g0
+        # Trace tr(g2^2) = (g2^a_b)(g2^b_a): raise and trace the (0,2) square.
+        g2_sq_up = g2_sq.up(g0, 0)
+        tr_g2_sq = g2_sq_up.trace(0, 1)
+        # Pure-trace piece (1/16) tr(g2^2) g0 => tr g4 = (1/4) tr(g2^2).
+        g4 = (1 / 16) * tr_g2_sq * g0
         g4.set_name('g4')
         return g4
 
-    prefactor = 1 / (n - 4)
-    correction = (1 / (4 * (n - 1))) * tr_g2_sq * g0
-    g4 = prefactor * (g2_sq - correction)
+    # Algebraic piece, exact on locally conformally flat boundaries:
+    g4 = (1 / 4) * g2_sq
     g4.set_name('g4')
     return g4
 
diff --git a/conformal_toolkit/poincare_einstein/renormalized_volume.py b/conformal_toolkit/poincare_einstein/renormalized_volume.py
index 08a13c8..919e52f 100644
--- a/conformal_toolkit/poincare_einstein/renormalized_volume.py
+++ b/conformal_toolkit/poincare_einstein/renormalized_volume.py
@@ -7,9 +7,14 @@
     Vol_epsilon ~ sum_k c_k(n) * epsilon^{k-n} * int v_k dV
 
 v_0 = 1 (volume form coefficient, trivial)
-v_2 = -(1/(n-2)) * J  (related to Schouten trace J = R/(2(n-1)))
+v_2 = -(1/2) * J  (related to Schouten trace J = R/(2(n-1)))
 
 Returns the LOCAL density (integrand), not an integrated quantity.
+
+The v_2 coefficient is n-INDEPENDENT: v_2 = -J/2 (Graham). The earlier
+-(1/(n-2)) J prefactor was wrong for every n != 4 (it happens to coincide
+at n=4, both giving -1 on S^4, which masked the bug in 4D-only tests), and
+the n=2 special case had the wrong sign (+J/2). See ERRATA C2 and M12.
 """
 from conformal_toolkit.core.schouten import schouten_trace
 
@@ -25,7 +30,7 @@ def renormalized_volume_coefficient(g0, order=2):
         Boundary metric of dimension n.
     order : int
         0 – returns 1 (trivial v_0 density).
-        2 – returns -(1/(n-2)) * J where J = R/(2*(n-1)).
+        2 – returns -(1/2) * J where J = R/(2*(n-1)).
 
     Returns
     -------
@@ -41,14 +46,11 @@ def renormalized_volume_coefficient(g0, order=2):
         return Integer(1)
 
     if order == 2:
-        n = g0.domain().dim()
-        if n == 2:
-            # In 2D the 1/(n-2) pole is regularized; the coefficient is finite
-            # and given by the conformal factor convention: v_2 = J/2.
-            J = schouten_trace(g0)
-            return J / 2
+        # v_2 = -J/2, independent of n (= -n/4 on the round S^n). Verified by
+        # expanding the sphere volume density (1 - rho^2/4)^n, whose rho^2
+        # coefficient is -n/4 = -J/2 with J = n/2. See ERRATA C2/M12.
         J = schouten_trace(g0)
-        return -(1 / (n - 2)) * J
+        return -(1 / 2) * J
 
     raise NotImplementedError(
         f"Renormalized volume coefficient v_{order} not implemented (only 0 and 2)"
diff --git a/conformal_toolkit/tractor/standard_tractor.py b/conformal_toolkit/tractor/standard_tractor.py
index 384f8e4..40103a4 100644
--- a/conformal_toolkit/tractor/standard_tractor.py
+++ b/conformal_toolkit/tractor/standard_tractor.py
@@ -5,10 +5,15 @@
 
     I^A = (sigma, mu_a, rho)
 
-where:
-- sigma is a scalar field (conformal weight +1 density)
-- mu_a is a 1-form (conformal weight -1)
-- rho is a scalar field (conformal weight -1)
+where, in the composition series E[1] -> E_a[1] -> E[-1]:
+- sigma is a scalar field   (conformal density weight +1, in E[1])
+- mu_a is a weighted 1-form  (conformal density weight +1, in E_a[1])
+- rho is a scalar field     (conformal density weight -1, in E[-1])
+
+(The middle slot mu_a carries weight +1, NOT -1: the file's own
+tractor_metric.py contracts g^{ab} mu_a mu'_b alongside sigma*rho', and
+homogeneity of that pairing forces w(mu) = +1. The earlier "-1" label was
+a copy-paste from the rho slot -- see ERRATA m4.)
 """
 
 
@@ -41,7 +46,7 @@ def sigma(self):
 
     @property
     def mu(self):
-        """Middle slot: 1-form (weight -1)."""
+        """Middle slot: weighted 1-form (density weight +1; see ERRATA m4)."""
         return self._mu
 
     @property
diff --git a/docs/TOOLING_GAPS.md b/docs/TOOLING_GAPS.md
new file mode 100644
index 0000000..9e42111
--- /dev/null
+++ b/docs/TOOLING_GAPS.md
@@ -0,0 +1,120 @@
+# Math-Tooling Gaps & Verification Roadmap
+
+This audit was as much a stress-test of *our verification tooling* as of the
+library. This document records where the tooling fell short, and the concrete
+plan to close each gap so that the checks in [`../ERRATA.md`](../ERRATA.md) run
+automatically — not by hand — next time.
+
+The goal is a **verification harness that fails loudly**: a dead backend or an
+unverifiable claim should break CI, never silently degrade a proof to "trust
+me."
+
+---
+
+## Root cause: silent backend degradation
+
+Every symbolic backend was unavailable during the audit:
+
+| Tool | Failure | Consequence |
+|------|---------|-------------|
+| `math_verify_identity`, `math_simplify` | `Neither Lean 4 nor SymPy available` | No symbolic identity checks |
+| `math_compute` | `sympy unavailable` | Numeric checks fell back to ad-hoc `python3` |
+| `wolfram_query` | `no_api_key` | No exact closed-form / factorial verification |
+
+All scalar checks were re-run in plain Python; all tensor identities were done
+by hand. **It still worked** — but only because the reviewers happened to know
+the closed forms. A less-expert pass would have produced false "verified"
+stamps. That is the real risk.
+
+**P0 fix — make the stack self-healthchecking.** Bundle a pinned `sympy` (and
+optional Lean 4 + mathlib) into the Dynamo MCP server image; have
+`math_stack_status` return non-zero when a backend is down, and gate any
+verification workflow on it so a dead backend *fails the job* instead of
+silently downgrading it. Provision `WOLFRAM_ALPHA_APPID` from Keychain
+(`cobalt-vault/wolfram-appid`, per house keychain policy) and surface
+`wolfram_quota_status`.
+
+---
+
+## Gap 1 — No symbolic tensor / abstract-index calculus verifier *(P0 for this domain)*
+
+The single biggest gap. There is no SageManifolds/`xAct`/`cadabra2`-class engine
+behind the math tools, so **every** tensor claim had to be hand-reduced to a
+scalar on `Sⁿ`: the Bach contraction order, `|P|²`, `Δ` on a `(0,2)` tensor, the
+Thomas-D slots, tractor-connection parallelism, "middle curvature block = Weyl."
+
+**Proposal — `math_tensor_verify`.** Wrap SageManifolds (or `cadabra2`).
+Inputs: a metric/connection plus an abstract-index expression and an identity to
+check. Output: the simplified residual tensor and a per-component pass/fail.
+Would have caught `M2` (missing `R`), the `M3`/`M8`/`M12` signs, and the
+`C1`/`M11` trace coefficients **mechanically**.
+
+## Gap 2 — No conformal-weight / covariance checker *(P1, high leverage)*
+
+Several errors are pure **density-weight bookkeeping**: `M4` (linear-`H` term
+and bidegree-sign), `M9` (`£_v h` homogeneity), `m4` (μ-slot weight). These need
+no heavy CAS.
+
+**Proposal — `conformal_weight_check`.** Given an operator and its slot weights,
+verify homogeneity, bidegree *direction* (`target − source` sign), and conformal
+covariance under `g → e^{2ω}g`. Cheap, fast, and it mechanizes exactly the
+reasoning that caught `M4` by eye.
+
+## Gap 3 — No degenerate-metric / Carrollian Lie-derivative engine *(P2)*
+
+`£_ξ` of a *degenerate* `(0,2)` tensor and the conformal-Carroll pair condition
+`(£_ξ h = 2λh, £_ξ v = −λv)` had no symbolic support.
+
+**Proposal — extend Gap 1's verifier with degenerate-metric support** plus a
+`carroll_symmetry_check(ξ, h, v)` predicate. Directly validates `M9`/`M10`/`m5`.
+
+## Gap 4 — No conformal-invariant basis enumerator *(P2)*
+
+No Weyl-invariant-theory / Hironaka-decomposition tool, so the weight-4
+hypersurface count (`M5`) could only be bounded below (`> 4`), not pinned.
+
+**Proposal — `conformal_invariant_basis(weight, ambient_dim, hypersurface=bool)`**
+returning an independent generating set (algebraic + derivative invariants).
+Hard, but it is what makes `count_invariants` checkable exactly.
+
+## Gap 5 — No library-execution harness in CI *(P1, highest ROI after backends)*
+
+We could not run the repo's own Sage code, so prose-vs-return mismatches went
+unverified (does `q_curvature(order=4)` literally return 6? does
+`count_invariants` return the stated ints?).
+
+**Proposal — a Sage-enabled CI job** that evaluates each operator on canonical
+anchors — round `Sⁿ`, `S²×S²`, **and a non-conformally-flat metric** — and
+asserts known values (Branson `Q_n=(n−1)!`, conformal invariance / vanishing).
+This catches `M2`, `M8`, `C1`/`C2`, `M16` **empirically, with no symbolic engine
+at all.** The anchor library is already specified in `ERRATA.md → How we caught
+them`; turning it into `tests/test_anchors/` is the concrete first deliverable.
+
+## Gap 6 — Literature retrieval isn't citation-grade *(P3)*
+
+`world_arxiv` returned irrelevant collaboration papers; `WebFetch` couldn't pull
+numbered equations from arXiv PDFs (binary / ar5iv paraphrase only). This left
+`M5`/`M7` (Fialkow, GW counts) and the Graham–Hirachi constant (`m3`) at
+moderate confidence.
+
+**Proposal — an arXiv *source* fetcher** that pulls the LaTeX source (not the
+rendered PDF) and indexes numbered equations for verbatim citation lookup.
+
+---
+
+## Priority summary
+
+| Gap | Priority | Effort | Catches |
+|-----|----------|--------|---------|
+| Backend healthcheck + Wolfram key | **P0** | Low | the silent-degradation risk itself |
+| Tensor / abstract-index verifier | **P0** | High | `M2`, `M3`, `M8`, `M12`, `C1`, `M11` |
+| Library-execution CI harness | **P1** | Low–Med | `M2`, `M8`, `C1`, `C2`, `M16` (empirically) |
+| Conformal-weight checker | **P1** | Low | `M4`, `M9`, `m4` |
+| Degenerate-metric / Carroll engine | **P2** | Med | `M9`, `M10`, `m5` |
+| Invariant-basis enumerator | **P2** | High | `M5`, `M6` |
+| arXiv source fetcher | **P3** | Low | `M5`, `M7`, `m3` |
+
+**Do first:** the backend healthcheck (P0) and the library-execution CI harness
+(P1). Together they are low-effort and would have caught the two *critical*
+errors and most majors automatically — the highest return per hour of the whole
+list.
diff --git a/paper.md b/paper.md
index c8f7b46..bcee371 100644
--- a/paper.md
+++ b/paper.md
@@ -84,7 +84,7 @@ symbolic conformal geometry to discrete mesh representations usable in PyTorch.
 # Functionality
 
 The software consists of two Python packages, 40 source modules in total, with
-125 tests and 6 example notebooks.
+160 tests and 6 example notebooks.
 
 ## Symbolic Package: `conformal_toolkit`
 
diff --git a/tests/test_discrete/test_errata_anchors.py b/tests/test_discrete/test_errata_anchors.py
new file mode 100644
index 0000000..98f56f9
--- /dev/null
+++ b/tests/test_discrete/test_errata_anchors.py
@@ -0,0 +1,45 @@
+"""Regression anchors for the May-2026 errata (discrete package).
+
+Each test pins a CORRECTED behaviour on a concrete geometry, in the spirit of
+ERRATA.md -> "How to contribute a correction": verified-on-an-anchor beats
+argued-in-prose. These cover only the PyTorch-side fixes (runnable without
+SageMath); the symbolic-side anchors live in docs/TOOLING_GAPS.md (Gap 5).
+"""
+import torch
+from tests.conftest_mesh import _make_icosphere
+from conformal_features.discrete.q_curvature import discrete_q_curvature
+from conformal_features.discrete.bach import discrete_bach_norm
+
+
+def test_discrete_q4_converges_to_zero_not_six():
+    """ERRATA M16: discrete order-4 feature is H^2 - K, which -> 0 on a sphere
+    (it is the Willmore integrand, NOT the 4D GJMS Q_4 = 6)."""
+    means = []
+    for s in (2, 3, 4):
+        v, f = _make_icosphere(subdivisions=s)
+        means.append(float(discrete_q_curvature(v, f, order=4).mean().abs()))
+    # Small, and monotonically shrinking under refinement -- i.e. -> 0, not 6.
+    assert means[0] < 0.1, f"expected ~0 on icosphere, got {means[0]}"
+    assert means[-1] < means[0], f"should shrink under refinement, got {means}"
+    assert all(m < 1.0 for m in means), f"never approaches 6: {means}"
+
+
+def test_discrete_q4_radius_independent():
+    """ERRATA m6: H^2 - K = 0 on a round sphere of ANY radius, not just R=1."""
+    v, f = _make_icosphere(subdivisions=3)
+    for R in (0.5, 1.0, 3.0):
+        q = discrete_q_curvature(v * R, f, order=4)
+        assert float(q.mean().abs()) < 0.1, f"R={R}: expected ~0, got {q.mean()}"
+
+
+def test_bach_integrated_default_small_pointwise_option_exists():
+    """ERRATA M17: the default integrated proxy |L L K| is scale-stable/small on
+    a sphere; the mathematically-correct pointwise M^-1 L M^-1 L K is offered as
+    an option (and is noisier on coarse meshes -- that's why it isn't default)."""
+    v, f = _make_icosphere(subdivisions=3)
+    integrated = discrete_bach_norm(v, f)                  # default
+    pointwise = discrete_bach_norm(v, f, pointwise=True)   # option exists
+    assert integrated.shape == pointwise.shape == (v.shape[0],)
+    assert float(integrated.mean()) < 1.0
+    # The pointwise operator amplifies icosphere pentagon-defect noise.
+    assert float(pointwise.mean()) > float(integrated.mean())