B modes

NOTE_lowl_bb_excess.md

@@ -0,0 +1,132 @@
+# low-ℓ raw tensor-BB ~0.4% excess vs CLASS — ROOT (half) FOUND 2026-06-16
+
+**Not** the §5 ℓ=477 item (resolved: CubicSpline interpolation; see bottom).
+
+## ROOT CAUSE (≈half the excess): tensor integration STARTS TOO LATE.
+
+`TensorPerturbationEvolver.evolution_one_k` (tensors.py:376) caps the start at
+lna=-10 (z~22000) — too late for the photon polarization quadrupole Π to settle
+before recomb. Fixed-start sweep (`diag_bb_starttime.py`): BB is start-converged
+only for t0 ≤ -13 (ℓ=10: cap -10 → +4.20e-3; -11 → +2.98e-3; -13/-16/-22 →
++2.00e-3). Fix = cap at -14 (`proposed_diff_tensor_starttime.patch`).
+End-to-end verify (`diag_bb_starttime_verify.py`): improves EVERY node ℓ=2-450,
+**zero wall-clock cost** (6.84→6.87s); low-ℓ max 4.20e-3→3.09e-3, recomb max
+1.59e-3→0.85e-3. Earlier start is the more-correct answer (closer to CLASS,
+which is start/TCA-independent). The SAME -10 cap is in the scalar evolver
+(perturbations.py:279) — NOT changed (scalars pass + are far less sensitive;
+re-validate before touching).
+
+**After the fix, two residuals remain:** (a) ℓ=2-4 ~+3.0e-3 (reion-bump region,
+start-insensitive — likely the ~0.5% HyRex-vs-C-HYREC-2 visibility-wing diff,
+see #9); (b) a converged ~1-1.7e-3 recomb-region residual (the genuine
+inter-code Π difference below, now ~halved). Both are within inter-code scatter
+for raw tensor BB; (a) is the next target if pursued.
+
+## The symptom
+
+ABCMB raw tensor BB (lensing off, r=0.1, fiducial LCDM) is systematically
+**~0.4% HIGH at low ℓ** vs high-precision tensor-only CLASS
+(`class_tensor_hp_reference()`), decaying smoothly to the recomb-tail floor:
+
+| ℓ | 2 | 10 (peak) | 64 | 100 | 152 | 237 | 331 | 490 |
+|---|---|---|---|---|---|---|---|---|
+| ABCMB/CLASS−1 | +3.3e-3 | **+4.2e-3** | +2.6e-3 | +1.8e-3 | +9e-4 | +1.3e-3 | +5e-4 | ~0 |
+
+At the spline **nodes** (computed exactly) — a real transfer-level difference,
+not interpolation. Smooth, single-signed (ABCMB high), broad, peaks at ℓ≈10,
+→0 by ℓ≈490 (so it is k-dependent, largest for low-k / super-horizon-at-recomb
+modes).
+
+## CONCLUSION (this session): genuine inter-code difference in the tensor
+## moment evolution, established in the tight-coupling epoch. NOT a TCA effect,
+## NOT recombination, NOT a localizable bug.
+
+The excess lives in the **polarization source** `source_E = √6·g·Π`. Perturbation-
+level ABCMB-vs-CLASS at recomb (`diag_bb_pi_pert_xcheck.py`,
+`diag_bb_tca_origin.py`): GW amplitude **h is exact (+0.004% at ℓ=10)**, but the
+**polarization quadrupole Π is +0.1% high**, which squares up to +0.4% in
+`Cl_BB ∝ |∫√6·g·Π·radB|²`. EE carries the identical excess (shared source); TT
+carries ~⅔ (its `−hdot·e^{−κ}` term is cleaner). The Π residual is established
+by z≈1400 (tight-coupling era) and *decreases* toward recomb (z=1400 +0.16% →
+z=950 +0.08%) — the flavor of an early-time IC transient. It is in the larger
+moments (δ_g, G0, G2), NOT the tiny frozen shear_g.
+
+**Who is right (h is exact, so this is the polarization quadrupole only):**
+undetermined, but it is NOT a TCA artifact on either side. Pushing CLASS's TCA
+triggers 1.5e-3 → 1e-7 (CLASS runs full-hierarchy for ~all the evolution)
+moves CLASS BB by **+0.006%** — i.e. TCA accounts for only ~1.4% of the 0.4%
+gap (`diag_class_minimal_tca.py`). CLASS is TCA-independent/converged; ABCMB's
+full-hierarchy result differs by +0.4% for a non-TCA reason. For raw low-ℓ
+tensor BB (a tiny observable) a 0.4% inter-code difference is within normal
+CLASS-vs-CAMB-class scatter.
+
+## Ruled OUT this session (with evidence; do NOT re-investigate)
+
+1. **Reionization.** `diag_bb_srcgrid_reion.py` part C: reion fraction of BB
+   swings 99%→0% from ℓ=2→12, excess stays flat +4e-3. Decisive.
+2. **k-transfer grid** (`diag_bb_kgrid_conv.py` A): x1/x2/x4 → ℓ=10 stable 1e-5.
+3. **k-source/perturbation grid** (`diag_bb_srcgrid_reion.py` B): source
+   x1/x2/x4 re-solve, no trend.
+4. **lna/time grid** (`diag_bb_lna_conv.py`): Nlna 500→4000 changes ℓ=10 by <1e-8.
+5. **CLASS under-convergence at low ℓ** (`diag_class_lowl_ladder.py`): CLASS
+   converged to <0.01% across k/q/tol/time/TCA; converging it *widens* the gap.
+6. **Equations / ICs / radials / source-formula / GW-eq / truncations.**
+   Line-by-line vs CLASS perturbations.c + transfer.c, flat limit: all verbatim
+   (P2, F/G/U hierarchies, l_max closures, radB=½(j'+2j/x), GW eq, IC=gw_ini/√6).
+7. **GW-source densities** (`diag_gw_source_norm.py`): ρ_unit·ρ_g, ρ_unit·ρ_ν,
+   and ρ_ν/ρ_g all match CLASS to ~2e-6 (ppm). Not a normalization bug.
+8. **GW amplitude h**: exact (+0.004% at ℓ=10, perturbation level).
+9. **Recombination is NOT the driver.** Both ABCMB (HyRex) and default CLASS
+   use HYREC-2 (CLASS 3.3.4 default `recombination=hyrec`, input.c:5772). ABCMB
+   xe(z) is only +0.07% off C-HYREC-2 (`diag_recomb_hyrec_recfast.py`); the
+   HyRec↔RecFast spread moves BB only ~0.07%; ABCMB BB sits +0.42% above HyRec
+   AND +0.35% above RecFast (above *both*). Recomb explains ≲0.07% of the 0.4%.
+   (The visibility g wings differ ~0.5% antisymmetrically but that mostly
+   cancels in the BB integral — proven by the small HyRec↔RecFast BB spread.)
+10. **Solver under-resolution of the tight-coupling quadrupole**
+    (`diag_bb_dtmax_test.py`): shear_g is frozen at ~3e-12 (below atol_ten=1e-9)
+    through tight coupling, BUT forcing resolution (dtmax 0.005, atol 1e-13)
+    changes BB by <2e-5. shear_g is real-but-irrelevant (negligible Π weight).
+11. **Tensor TCA (CLASS-side) / gw_source photon-drop during TCA.** See
+    conclusion: CLASS TCA-independent → not it. The photon gw_source term CLASS
+    drops during TCA is negligible (δ_g, shear_g suppressed in tight coupling).
+12. **Scalar EE cross-check** (`diag_scalar_ee_xcheck.py`): scalar EE is clean
+    (mean −3e-4 over ℓ=100–1500), so this is tensor-specific, not a shared
+    visibility/recomb bug at the 0.4% level.
+
+## Only untested structural difference left
+
+**Early-time start / IC transient.** ABCMB starts every tensor mode at
+lna ≤ −15 (z > 3e6) with exactly-zero moments (h=1/√6); CLASS starts later,
+per-mode. The Π bias being largest early and decaying toward recomb is
+consistent with an IC transient that depends on the start epoch. Test: sweep
+ABCMB's tensor t0 (force earlier/later, check BB start-time convergence). If
+insensitive → irreducible accumulated inter-code difference; if sensitive →
+the start criterion. Not yet run.
+
+## How the test handles it (`pytests/accuracy_test_bb.py::test_bb_raw`)
+
+- recomb nodes ℓ≥100: held 2.5e-3 (actual ~1.6e-3).
+- low-ℓ nodes 3≤ℓ<100: held 5e-3 (actual 4.2e-3 — this excess).
+- full band: 1e-2 (interp-limited ℓ=477 residual). ℓ=2 exempt; 491–500 excluded.
+
+Since the excess is NOT a fixable bug (verbatim-correct tensor code; within
+inter-code scatter for raw tensor BB), the 5e-3 low-ℓ bar should STAY, with the
+comment updated to "converged tensor-quadrupole inter-code difference, not
+reion/grid/recomb/TCA" rather than "unknown structural diff." Do NOT tighten to
+2.5e-3.
+
+Diagnostics in worktree root: diag_bb_kgrid_conv.py, diag_bb_srcgrid_reion.py,
+diag_bb_eett_xcheck.py, diag_class_lowl_ladder.py, diag_bb_lna_conv.py,
+diag_bb_visibility_xcheck.py, diag_bb_pi_pert_xcheck.py, diag_scalar_ee_xcheck.py,
+diag_recomb_hyrec_recfast.py, diag_gw_source_norm.py, diag_bb_tca_origin.py,
+diag_bb_dtmax_test.py, diag_class_minimal_tca.py (+ logs).
+
+---
+
+### Resolved earlier (2026-06-15): §5 ℓ=477 interpolation residual
+
+raw-BB 6.9e-3 at ℓ=477 = CubicSpline interpolation error over the 40-wide node
+gap (450→490), not transfer physics. Proven by `diag_bb_dense_hp.py`
+(CLASS-internal self-interp reproduces ~7e-3 at ℓ=477, 0 at nodes). Nodes are
+maximally dense (bessel_l_tab); finer needs new Bessel tables.

SCOPE_bessel_merge.md

@@ -0,0 +1,151 @@
+# Scope: bring the curvature-branch Bessel recurrence into `bmodes`
+
+**Goal (user):** integrate the curvature branch's changes so the B-mode path
+uses the every-ℓ hyperspherical-Bessel recurrence instead of the phi-tables +
+sparse-ℓ CubicSpline. This kills the §5 ℓ=477 raw-BB residual (pure spline
+error over the 40-wide node gap) and drops the tensor solver's table
+dependence. **`origin/curvature` (== local `worktree-curvature`, @ `3174b25`)
+is authoritative and read-only; only `bmodes` is modified.**
+
+Both branches are siblings off `5eabbab`. Verified `origin/curvature` and
+`worktree-curvature` are the *same commit*.
+
+---
+
+## 1. What each branch changed (vs base `5eabbab`, `abcmb/` only)
+
+| File | curvature | bmodes | overlap? |
+|------|-----------|--------|----------|
+| `spectrum.py` | rewrote `get_Cl` → recurrence; **deleted** `Cl_one_ell`, sparse-ℓ spline, `ells_indices`, module-level phi-tables + `j`, the ℓ-cap | added `tensor_cls` arg to `get_Cl`; `lensed_Cls` → 4-tuple EE↔BB mixing | **YES — the real merge** |
+| `main.py` | `+curvature=` kwarg to SS; `omega_k`/`K` derived params; `omega_Lambda −= omega_k`; `expected_keys += omega_k,K` | TPE/TSS fields+init; tensor wiring in `_run_post_recomb`; `r`/`n_t` derived; `Output.ClBB`; `expected_keys += r,n_t` | trivial (only `expected_keys` literal collides) |
+| `model_specs.py` | curvature specs; `species.Curvature`; curved k-grids (`omega_k_ref`, `closed_integer_nu`, `K_ref`) | tensor specs block | none (disjoint regions) |
+| `ABCMBTools.py` | +109: `sin_K`, `cot_K`, `_curv_f`, `_curv_g`, `_curv_g_diff` | — | clean (curvature-only) |
+| `species.py` | +179: `Curvature` species + curved ρ/ICs | — | clean |
+| `background.py`, `perturbations.py` | curved distances / scalar ICs | — | clean |
+| `bessel_tab/*.txt` | **deleted** (79 MB) | — (still imported by `tensors.py`!) | see §3 |
+| `tensors.py` | — (does not exist on curvature) | +750 new | — |
+
+**Net:** a `git merge origin/curvature` into `bmodes` applies cleanly
+everywhere except `spectrum.py` (substantial, expected) and a one-line
+`expected_keys` collision in `main.py` (keep both: `omega_k,K` + `r,n_t`).
+The bulk of the curvature feature (new species, curved background, curved
+k-grids, ABCMBTools helpers, table deletions) lands without conflict.
+
+---
+
+## 2. The recurrence is a clean fit for the tensor radials
+
+`spectrum.py:_Cl_all_ells_curved` (curvature) walks ℓ via the three-term
+hyperspherical recurrence, and its `step` already computes, per (Nlna, Nk):
+
+- `Phi_l`            = Φ_ℓ           → at K=0 = **jₗ(kχ)**
+- `dPhi`  = dΦ/dχ    → `dPhi/k`      → at K=0 = **jₗ′(kχ)**
+- `d2Phi` = d²Φ/dχ²  → `d2Phi/k²`    → at K=0 = **jₗ″(kχ)**
+- `x_eff` = q·S_K(χ) → at K=0 = **kχ = x**
+
+The tensor radial functions (`tensors.py:Cl_one_ell`, CLASS `transfer.c`) are
+*exactly* combinations of those quantities (flat limit):
+
+```
+radT = sqrt(3/8 (l+2)(l+1)l(l-1)) · jₗ / x²
+radE = 1/4 [ jₗ″ + 4 jₗ′/x − (1 − 2/x²) jₗ ]
+radB = 1/2 [ jₗ′ + 2 jₗ/x ]
+```
+
+So the port is: inside a recurrence `step`, form `radT/radE/radB` from
+`Phi_l, dPhi/k, d2Phi/k², x_eff`, contract against the tensor sources
+`source_T2, source_E` (and primordial `P_h = r A_s (k/k_pivot)^{n_t}`),
+accumulate `ClTT_t/ClTE_t/ClEE_t/ClBB_t` at **every** integer ℓ in
+`2..l_tensor_max`, zero-pad to `out_ells`. No tables, no spline → the ℓ=477
+residual is gone by construction.
+
+**Tensors stay flat (K=0).** The scalar Φ recurrence is *not* the correct
+spin-2 radial basis in curved space, so the tensor port is validated only at
+`omega_k = 0`. `omega_k ≠ 0` AND `tensors=True` must be guarded/documented as
+unsupported (both are opt-in, default off, so no regression). Carrying `K`
+through the tensor recurrence (reusing `tools.sin_K/cot_K`) is fine because it
+is exact at K=0; just don't claim the curved-tensor result.
+
+---
+
+## 3. The hard dependency that forces the choice
+
+`tensors.py:11` imports `bessel_l_tab, xphi0_tab, phi0_tab, … , j` from
+`spectrum.py`. The curvature branch **deleted all of those**. So after merging
+curvature, `tensors.py` will not import until the recurrence port is done —
+the port is *mandatory*, not optional, once curvature is merged.
+
+---
+
+## 4. Two ways to scope it
+
+### Option A — tensor-only recurrence, NO curvature merge (minimal)
+Write a flat (K=0) spherical-Bessel recurrence self-contained in `tensors.py`;
+replace `TensorSpectrumSolver.{Cl_one_ell, get_Cl}` sparse-node+spline with the
+every-ℓ recurrence; drop the table imports. Leave bmodes' scalar `spectrum.py`
+and the tables untouched (scalars keep tables+spline; they pass at 1%).
+- **Touches:** `tensors.py` only (+ maybe copy `sin_K`/`cot_K` or hardcode flat forms).
+- **Fixes:** ℓ=477 raw-BB residual; tensor table dependence.
+- **Does NOT:** reconcile bmodes with curvature (bmodes keeps tables for
+  scalars; curvature deleted them → the two still conflict on a future main
+  merge). Does NOT bring omega_k.
+- **Risk:** low. Smallest, most surgical.
+
+### Option C — full `git merge origin/curvature` into bmodes, then port (alignment)
+Merge curvature wholesale (resolving spectrum.py + the trivial main.py
+collision), then port the tensor recurrence in `tensors.py`. Scalars **and**
+tensors use the recurrence; tables deleted; bmodes gains omega_k; the two
+branches are reconciled (what "bring bmodes into alignment with curvature"
+literally asks).
+- **Touches:** `spectrum.py` (merge: adopt curvature's recurrence `get_Cl`,
+  graft bmodes' `tensor_cls` add-in + 4-tuple `lensed_Cls`), `main.py`
+  (trivial), `tensors.py` (port), + accept the rest of curvature clean.
+- **Fixes:** ℓ=477 for BB *and* removes the scalar sparse-ℓ spline; full
+  table deletion; mergeable-with-curvature.
+- **Risk:** medium. Biggest piece is the `spectrum.py` 3-way reconciliation
+  and re-validating the scalar path (curvature's `get_Cl` must still produce
+  the bmodes 4-tuple incl. tensor add-in + lensed BB).
+
+**Recommendation:** Option C — it is what the user asked ("integrate the
+changes in the curvature branch … bring it into alignment"), and once
+curvature is merged the tensor port is required anyway (§3). Option A is the
+fallback if the spectrum.py reconciliation proves larger than expected.
+
+---
+
+## 5. Execution plan for Option C
+
+1. **Baseline first:** confirm reverse-AD through `tensors=True` is finite on
+   the *current* bmodes (in flight) and record scalar `accuracy_test.py` +
+   `accuracy_test_bb.py` + tensors=True warm time as the pre-merge reference.
+2. `git merge origin/curvature` on `bmodes`. Resolve:
+   - `main.py`: union `expected_keys` (`omega_k,K` + `r,n_t`); keep both
+     add_derived blocks; keep both SS-init kwargs (curvature= AND the TPE/TSS
+     init is separate).
+   - `model_specs.py`: should auto-merge (disjoint); verify tensor specs +
+     curvature specs both present.
+   - `spectrum.py`: take curvature's recurrence `get_Cl`; re-introduce the
+     `tensor_cls=None` arg, the "add tensor to unlensed totals + bb_unlensed"
+     block, and switch `get_lensed_Cls` to the 4-tuple `lensed_Cls`; keep
+     bmodes' 4-tuple `lensed_Cls` body (curvature didn't touch it).
+3. **Port `tensors.py:TensorSpectrumSolver`** to the recurrence (§2): new
+   `_Cl_all_ells_tensor`, drop table imports + spline + `Cl_one_ell`. Keep
+   `TensorPerturbationEvolver` untouched (it is table-independent).
+4. **Validate:** `accuracy_test.py` (scalar path unchanged), `accuracy_test_bb.py`
+   (ℓ=477 band bar can tighten now — every ℓ computed; nodes already
+   sub-permille), reverse-AD `tensors=True` finite both lensings, tensors=True
+   warm time. All on GPU.
+5. The §5 / band-tolerance test assertions in `accuracy_test_bb.py` can be
+   tightened from 1e-2 once the spline is gone (the interp floor was the only
+   reason for the loose 1e-2). The low-ℓ 4e-3 bar STAYS (inter-code scatter,
+   not interp).
+
+## 6. Open risks
+- **curved + tensors** unsupported (flat-only tensor radials) — guard it.
+- **reverse-AD**: curvature's recurrence is revAD-validated at K=0 for scalars;
+  the tensor port reuses the same `lax.scan`/`jax.checkpoint` chunked pattern,
+  so it should inherit that — but must be re-tested (ties to task #1).
+- **k-grid**: tensor transfer currently truncates the scalar grid at the tensor
+  k_max (`get_tensor_k_axes`). The recurrence consumes `self.k_axis_transfer`;
+  confirm the tensor k-axis still feeds the recurrence correctly (the
+  contraction is over the tensor k-grid, not the full scalar one).

abcmb/ABCMBTools.py

@@ -360,3 +360,112 @@ def bilinear_interp(x, y, z, xq, yq):
     return (1 - tx) * (1 - ty) * z00 + tx * (1 - ty) * z01 + (1 - tx) * ty * z10 + tx * ty * z11
 
 
+
+
+### CURVED-GEOMETRY HELPERS ###
+
+# Generalized trig functions for curved FLRW as entire functions of u = K chi^2,
+# smooth through K = 0 (exact flat limit). K > 0 closed, K < 0 open. Series for
+# |u| < _CURV_U0; clipped sqrt args keep the unused branch finite under reverse AD.
+
+_CURV_U0 = 1.e-2
+
+
+def _curv_f(u):
+    """
+    f(u) = sin(sqrt(u))/sqrt(u), continued through u=0 (sinh for u<0, 1 at 0).
+
+    Parameters:
+    -----------
+    u : array
+        K * chi^2.
+
+    Returns:
+    --------
+    array
+        f(u); used as sin_K(chi) = chi * f(K chi^2).
+    """
+    series = 1. - u/6. + u**2/120. - u**3/5040. + u**4/362880.
+    sp = jnp.sqrt(jnp.clip(u, _CURV_U0, None))
+    sm = jnp.sqrt(jnp.clip(-u, _CURV_U0, None))
+    return jnp.where(u > _CURV_U0, jnp.sin(sp)/sp,
+                     jnp.where(u < -_CURV_U0, jnp.sinh(sm)/sm, series))
+
+
+def _curv_g(u):
+    """
+    g(u) = sqrt(u)*cot(sqrt(u)), continued through u=0 (coth for u<0, 1 at 0).
+
+    Parameters:
+    -----------
+    u : array
+        K * chi^2.
+
+    Returns:
+    --------
+    array
+        g(u); used as cot_K(chi) = g(K chi^2)/chi.
+    """
+    series = 1. - u/3. - u**2/45. - 2.*u**3/945. - u**4/4725.
+    sp = jnp.sqrt(jnp.clip(u, _CURV_U0, None))
+    sm = jnp.sqrt(jnp.clip(-u, _CURV_U0, None))
+    return jnp.where(u > _CURV_U0, sp/jnp.tan(sp),
+                     jnp.where(u < -_CURV_U0, sm/jnp.tanh(sm), series))
+
+
+def _curv_g_diff(a, b):
+    """
+    Cancellation-safe g(a) - g(b) for the Phi_1 seed.
+
+    Parameters:
+    -----------
+    a, b : array
+        Arguments of g (K chi^2 and q^2 chi^2 at the call site).
+
+    Returns:
+    --------
+    array
+        g(a) - g(b); factored series when both |a|, |b| < _CURV_U0.
+    """
+    both_small = jnp.maximum(jnp.abs(a), jnp.abs(b)) < _CURV_U0
+    series = (a - b) * (-1./3. - (a + b)/45. - 2.*(a**2 + a*b + b**2)/945.
+                        - (a**3 + a**2*b + a*b**2 + b**3)/4725.)
+    return jnp.where(both_small, series, _curv_g(a) - _curv_g(b))
+
+
+def sin_K(chi, K):
+    """
+    Curved-space comoving distance S_K(chi) = chi * f(K chi^2).
+
+    Parameters:
+    -----------
+    chi : array
+        Comoving radial distance, Mpc.
+    K : array
+        Curvature constant, Mpc^-2 (K > 0 closed, K < 0 open).
+
+    Returns:
+    --------
+    array
+        S_K(chi), Mpc: sin(sqrt(K) chi)/sqrt(K) (closed), sinh form (open), chi (flat).
+    """
+    return chi * _curv_f(K * chi**2)
+
+
+def cot_K(chi, K):
+    """
+    Generalized cotangent cot_K(chi) = g(K chi^2)/chi.
+
+    Parameters:
+    -----------
+    chi : array
+        Comoving radial distance, Mpc (kept > 0; diverges as 1/chi at 0).
+    K : array
+        Curvature constant, Mpc^-2.
+
+    Returns:
+    --------
+    array
+        cot_K(chi), Mpc^-1: sqrt(K) cot(sqrt(K) chi) (closed), coth form (open), 1/chi (flat).
+    """
+    return _curv_g(K * chi**2) / chi