fix(harmonic): correct BB radial kernel + add fine-k interpolation

[MinhMPA]

Summary

Two compounding bugs in compute_cl_bb were producing clax/CLASS C_l^BB ratios in [0.4×, 22×] depending on l. Both are fixed; primordial BB now matches CLASS sub-percent at l ≤ 200 and ~2% at l = 300.

Bugs

1. Wrong radial kernel. compute_cl_bb used √[l(l-1)(l+1)(l+2)] · j_l(x)/x², which is CLASS's TENSOR_TEMPERATURE_2 kernel (transfer.c:4241-4249), not the BB kernel. Replaced with TENSOR_POLARISATION_B (transfer.c:4263-4272, flat-space limit K=0 where cscKgen=cotKgen=1/x and rescale_function=1):

K_l^B(x) = 0.5 · (j_l'(x) + 2 · j_l(x) / x)

The derivative is computed via the recurrence j_l'(x) = j_{l-1}(x) - (l+1)/x · j_l(x), using spherical_jl_backward for both j_l and j_{l-1}.

2. k-grid undersampling for BB integration. compute_cl_bb integrated P_T(k)·|B_l(k)|² directly over the perturbation k-grid (~160 modes). The Bessel-driven oscillation period at the BB recombination peak (k ≈ 0.005-0.05 Mpc⁻¹) is comparable to the log-uniform k-mode spacing, so adjacent modes can sample opposing peaks of |B_l(k)|² and produce trapezoidal-rule errors of 6-30% with sign that flips with l. Added cubic-spline interpolation of source_p from the perturbation k-grid to a fine log-uniform k-grid (n_k_fine=2000 default), mirroring the pattern already used by compute_cls_all_fast for scalar T,E.

Both bugs are necessary and together sufficient — fixing only the kernel leaves 30-40% residuals at l ≈ 80-200; fixing only the k-grid leaves the kernel-form pathology (a spurious peak at the same scales).

Validation

tests/test_tensor.py::TestClBB::test_cl_bb_vs_class — ratio band tightened from [0.05, 20.0] to [0.95, 1.05] at l=2,10. Pass at the test's reduced TENSOR_PREC (l_max_g=10, k/decade=10, rtol=1e-3): l=2 → 0.993, l=10 → 1.006.

At production precision (l_max_g=30, 40 k/decade, rtol=1e-6, n_k_fine=2000), clax/CLASS BB ratios:

l	2	10	30	50	80	100	150	200	300
ratio	0.998	0.994	0.996	1.000	1.001	1.001	1.002	1.008	1.018

Sub-percent across l ≤ 200; ~2% at l = 300.

Why this didn't show up before

• tests/test_tensor.py only checked 0.05 < ratio < 20.0 — wide enough to pass even with the 22× kernel error at l=80-150.
• Lensed BB validates at <0.5% because at low r it is dominated by E→B conversion via C_l^φφ, which doesn't go through the tensor pipeline at all. The README's "lensing BB accurate, primordial BB ~2× off" line reflected this asymmetry without naming the cause.
• n_k_fine source interpolation was already in place for scalar T,E in compute_cls_all_fast — compute_cl_bb had been written before that pattern was generalized.

API

compute_cl_bb now takes an additional keyword-only n_k_fine=2000 argument. Existing positional callers compute_cl_bb(tpt, params, bg, l_values) continue to work unchanged.

Blast radius

grep -rn "compute_cl_bb\|from clax.harmonic" tests/ shows compute_cl_bb is consumed only by tests/test_tensor.py. Other tests in test_harmonic.py and test_high_l.py use compute_cl_tt, compute_cl_ee, compute_cl_te, compute_cls_all — all unaffected by this change.

Test plan

• tests/test_tensor.py::TestClBB — all 4 tests pass at TENSOR_PREC (28s).
• Tightened test_cl_bb_vs_class to ±5% — passes at l=2,10.
• Production-precision validation across l=2-300 (table above).
• Targeted regression run (tests/test_tensor.py tests/test_harmonic.py tests/test_lensing.py tests/test_constants.py tests/test_interpolation.py tests/test_background.py tests/test_thermodynamics.py tests/test_nonlinear.py tests/test_end_to_end.py --fast -x): 20 tests passed, 0 failures before the run was interrupted on a slow fixture; no regression observed.

Files

• clax/harmonic.py — rewritten compute_cl_bb (kernel + fine-k).
• tests/test_tensor.py — tightened test_cl_bb_vs_class ratio band.
• README.md — updated the "Known limitations" line for primordial BB.
• CHANGELOG.md — entry under May 4, 2026.