[kim-fp-03] Verify the off-diagonal Fokker-Planck Fourier kernel

[marjohma]

Summary

• replace the diagonal-only shortcut with the verified off-diagonal FP kernel
• stabilize Bessel/exponential evaluation and retain explicit model choices
• add oracle, symmetry, diagonal-collapse, normalization, and mutation tests

Why

Dense periodic Fourier assembly now uses the intended nonlocal FP response.

Validation

Validated cumulatively at the stack head:

• cmake --build build -j2
• ctest --test-dir build -E 'test_rhs_balance|test_periodic_convergence|test_periodic_vs_global' --output-on-failure — 35/35 passed

Known red gates are intentionally excluded from that command: the pre-existing test_rhs_balance baseline failure, the periodization-deformation acceptance gate, and the unresolved periodic-vs-global 5% agreement gate.

Stack

Base: fix/kim-fp-196-charge-kernel-oracle
Head: fix/kim-fp-187-offdiagonal-kernel

Closes #187

[marjohma]

Review (standards + spec vs #187) — keeping this in draft. The structural spec items are genuinely delivered: fp_ks_model with FP_KPERP_FULL/FP_KR_ONLY_APPROX matches the issue scaffold exactly (flr2_fourier_kernel.f90:16-56), the no-k_s form survives only as a named approximation, the test reads the committed oracle through CTest with good failure diagnostics, and the rho-Φ fixtures cover both charge signs, k_r≠k_r′, m_phi=−2..+2, finite k_s, b_x=600, zero-FLR, and Debye. What blocks it:

Blocking

1. The rho-B oracle assertion is vacuous. assert_oracle_complex uses |diff|/max(1,|expected|) > 1e-12 (test_flr2_fourier_kernel.f90:175-177), i.e. effectively an absolute 1e-12 tolerance — but all 13 committed rhoBDynamic magnitudes are ≤4.6e-15 (down to 1e-22). A fp_rho_B_harmonic stubbed to zero, a flipped sign, wrong Bessel order, dropped 4π, or wrong phase all pass. The issue's success criterion — "Mutating a sign, 4 pi, Fourier phase, k_s^2, or Bessel order must fail the test" — fails for the entire K^{ρB} half. Fix: scale the tolerance per-row to the row's magnitude (or normalize the committed rows).
2. K^{ρB} is not independently derived. The .wl addition (rhoBDynamic) transcribes the Fortran assembly (same −vT³/(λ²ω_c ν c) factorization) with no new derivation or cross-integration checks — only regenerated-vs-committed self-comparison. The issue's goal is to "replace the fused off-diagonal FP kernel's internal-consistency argument with an independent derivation"; that is currently true for rho-Φ (via the [math-04] Derive the charge-potential integral kernels #196 Weber/velocity-integral machinery) but not for rho-B.
3. ion_flr_scale_factor is silently bypassed in the new production path: fp_rho_phi_harmonic/fp_rho_B_harmonic recompute rho_L = abs(vT/omega_c) inline (flr2_fourier_kernel.f90 harmonic routines), while the stored ion rho_L includes the namelist ion_flr_scale_factor (species_mod.f90:511). The off-diagonal path now ignores a config option the diagonal path honors.
4. Production default flip is unplumbed and undocumented: the old kern_include_ks2=.false. default becomes fp_ks_model=FP_KPERP_FULL (physics change in the assembled kernels), the new knob is a bare module-save public integer — no namelist wiring, no README/nml documentation — and docs/plans/*phase1/3*.md still describe the old toggle. Also, the diagonal entry points (hatG_rho_phi/B_diag_sp) now route through flr_arg_pair_sp, so their production caller FLR2_asymptotics.f90:211-214 silently gains k_s² under the new default although the issue targets the off-diagonal kernel.

Also missing vs spec

• Negative fixture "an m_phi truncation whose estimated tail exceeds its configured tolerance: fail the validation case" — absent (only positive bound asserts at one synthetic point, not the standard q=3 case; runtime neither estimates tails nor rejects).
• "Define the near-axis treatment separately" is realized only as a non-finite-k_s guard (status=2) that can never fire on the production axis path (k_s there is finite/divergent, not NaN).

Minor

• Oracle columns phase, F0/F2, x1/x2, kpar/omega_E/Vpar/omega are read but never asserted in the Fortran test.
• Tests leave module state at FP_KR_ONLY_APPROX (the old "restore default" was removed) — fragile under test reordering.
• Collapse-test tolerances loosened 1e-12 → 5e-12 without a note.
• Line width >100 at test_flr2_fourier_kernel.f90:102,104,129; FORMULA_INDEX marks KIM-05/06 "[kim-fp-03] Verify the off-diagonal Fokker-Planck Fourier kernel #187 verified" though this PR touches neither file.

[marjohma]

Review follow-up fixes pushed (638fadb), addressing the blocking findings above. All changes developed red→green with the failure observed first; full suite (33/33 non-excluded) passes locally.

1. Vacuous rho-B oracle — fixed, and it found a real bug.

• Reproduced the vacuity first: with fp_rho_B_harmonic stubbed to return zero, the oracle subtest passed all 13 rows (only the internal-consistency collapse test caught it).
• assert_oracle_real/complex now scale per-row relative (|diff|/|expected|; exactly-zero rows must stay exactly zero), with a documented ORACLE_RELATIVE_TOLERANCE = 1e-12.
• The enforced oracle immediately failed on the pristine kernel with a scaled error of 2.6418e-8 — exactly the relative error of sol in constants_mod.f90, whose literal lacked _dp and was rounded through single precision (29979246592 vs 29979245800). Fixed (this is the same one-liner [math-06] Verify equilibrium and plasma-background formulas #218 carries, pulled forward with disclosure since the rho-B oracle requires it). Note: every c-dependent result shifts by ~2.6e-8, so the golden baseline will need a re-bless.
• After the fix the pristine max scaled error is 6.4e-14 (two orders inside the tolerance), and the issue's full mutation battery — sign, 4π, Fourier phase, Bessel order, hidden k_s² removal, each applied to the rho-B path — fails the ctest in every case (verified one mutation at a time).

2. ion_flr_scale_factor bypass — fixed. fp_rho_phi_harmonic/fp_rho_B_harmonic take rho_L as an explicit argument; hatG_* pass the stored per-species value, matching the diagonal path's flr_arg_pair_sp. Regression test added (test_harmonic_uses_stored_rho_L, watched fail before the fix).

3. fp_ks_model plumbing/docs — fixed. New &KIM_CONFIG key fp_ks_model_name ('kperp_full' default, 'kr_only' legacy), mapped in kim_read_config via fp_ks_model_from_name with a hard startup error on unknown names (mapping + rejection unit-tested). Documented in both shipped nmls and a new KIM/README.md section that states the default-behavior change relative to the old kern_include_ks2 = .false., and that the FLR arguments now honor ion_flr_scale_factor.

4. Truncation-tail negative fixture — added. New production gate fp_mphi_truncation_check (outermost-shell/total tail estimate; fails closed for sums too short to assess, i.e. mphi_max=0). The oracle's |m_phi|=2 bound now routes through it (positive arm); a slow-decaying shell sequence fails it (negative arm).

5. Smaller items. Infinite (on-axis) k_s now covered by the model-domain negative fixture alongside NaN; all kernel tests restore fp_ks_model to the production default (no more FP_KR_ONLY_APPROX leakage between tests); the three >100-col test lines wrapped; FORMULA_INDEX KIM-05/KIM-06 rows and the #187 paragraph reworded to claim only what is verified.

Still open (needs Mathematica, which I don't have here): the K^{ρB} entry in 04_charge_potential_kernels.wl remains a transcription of the Fortran assembly rather than an upstream derivation — the Fortran↔oracle agreement is now genuinely enforced at relative precision, but for rho-B that establishes consistency with the transcribed formula, not independent physics. FORMULA_INDEX now discloses this explicitly. A derivation pass (Weber/velocity-integral machinery as done for rho-Φ) plus .dat regeneration would close it; the tightened Fortran test needs no further changes to benefit.

[marjohma]

Mathematica derivation pass done (314a564) — the K^{ρB} entry in 04_charge_potential_kernels.wl is no longer a transcription.

Derivation. Starting from thesis (4.39) with drive (4.41) and gradient (13.23): the B drive u'_∥ B^r/B0 differs from the Φ drive −i c k_s Φ/B0 only in the parallel factor and coupling, so every shared factor (Fourier measures, 4π, Jacobian, Maxwellian norm, gyro phases, Weber J_⊥ integrals → F0/F2) cancels in the ratio. Per the susceptibility definition (5.12), the drive factor u' raises each source-side moment index by one and contributes one power of vT, mapping the verified rho-Φ assembly onto

K^{ρB} = −ph/(4π λ²) · vT³/(ω_c ν c) · (I01 (A1 F0 + A2 F2) + A2 I03 F0/2)

Two consequences:

1. The prefactor is vT³ — thesis (14.4) prints vT² which is a dimensional misprint (now disclosed in FORMULA_INDEX next to the known b_+/m_φ sign misprint). The Fortran was already correct.
2. The naive derivation produces the cubic source moment I03, not the I21 the code stores. These are equal: I03 = I21 is an identity of the energy-conserving FP moments, consistent with the (A.3)–(A.4) recurrence algebra (chain I^{03} = I^{12} = I^{21} via k↔l symmetry of (A.2) plus (A.4)). The WL script now extends the FP moment matrix to third order and gates every oracle case on I02 = I20, I03 = I21, and derived rho-B == production form before writing/checking rows.

Oracle unchanged. Derived and production forms agree to ~1e-49 relative, far below the 34 committed digits, so rho_phi_kernels.dat is byte-identical and the Fortran side needs nothing (no re-bless from this change; the sol-related golden re-bless noted earlier still stands independently).

Execution caveat + independent cross-check. The local Wolfram license on this machine turned out to be expired, so the extended .wl could not be executed here (make verify-mathematica-charge-kernels will exercise the new gates wherever a licensed kernel exists). To not leave the check unexecuted, the same derivation was run through an independent engine: an exact-series evaluation of the (A.2) generating integral (expand e^{−x1²s} in s=e^{−τ}; each term integrates to 1/(j+x1²−ix2) — no quadrature) in Python/mpmath at 50 digits, validated against the committed oracle intermediates:

check	max over 13 cases
replica vs committed I00/I20/I01/I21 columns	2.9e-34 (= 34-digit floor)
symmetry I02 = I20	exact (0.0)
identity I03 = I21	1.0e-49
re-assembled rho-Φ vs committed column	7.9e-34
derived rho-B vs committed rhoBDynamic	3.9e-34
derived rho-B vs production (I21) form	1.0e-49

mpmath reproduction script

#!/usr/bin/env python3
"""Independent verification of the K^{rho B} FP Fourier kernel (PR #217 / issue #187).

Derivation chain (Markl thesis, DOI 10.3217/efp2p-0x485):
  (4.39)  f_m(v_perp,u) = -Int du' G_m(u,u') df0/dr(v_perp,u') v^r_m(v_perp,u')
  (4.41)  v^r_m = (u'/B0) B^r_m - (i c k_s/B0) Phi_m
  (13.23) df0/dr_g = (A1 + (v_perp^2+u'^2)/(2 vT^2) A2) f0
  (5.12)  I^{kl}(x1,x2) = nu/(sqrt(2 pi) vT^{k+l+1}) Int du Int du' G(u,u')
                          exp(-u'^2/(2 vT^2)) u^k u'^l
  (5.13)-(5.15), (A.1)-(A.2): practical generating-integral definition.

I_D^{kl} series: with s = exp(-tau), the (A.2) integrand derivative at
alpha=beta=0 is P_kl(s; x1) * exp(x1^2) * exp(-x1^2 s) * exp((i x2 - x1^2) tau),
and Int_0^inf s^j exp((i x2 - x1^2) tau) dtau = 1/(j + x1^2 - i x2), so
  I_D^{kl} = exp(x1^2) * Sum_j d_j / (j + x1^2 - i x2),
  Sum_j d_j s^j = P_kl(s) exp(-x1^2 s)   (entire; truncated at J=80).
"""

import sympy as sp
import mpmath as mp

mp.mp.dps = 50
KMAX = 3
JMAX = 80
SOL = mp.mpf(29979245800)  # speed of light, cm/s (matches constants_mod fix)

a, b, s, x1s = sp.symbols('a b s x1')
h = (a + sp.I*x1s)*(b + sp.I*x1s)*(s - 1) + (a + b)**2/2
Pkl = {}
for k in range(KMAX + 1):
    for l in range(KMAX + 1):
        expr = sp.diff(sp.exp(h), a, k, b, l) * sp.exp(-h)
        poly = sp.Poly(sp.expand(expr.subs({a: 0, b: 0})), s)
        Pkl[(k, l)] = [(int(mono[0]), sp.lambdify(x1s, coef, 'mpmath'))
                       for mono, coef in zip(poly.monoms(), poly.coeffs())]

def i_diff_matrix(x1, x2):
    denom = [mp.mpc(j) + x1**2 - 1j*x2 for j in range(JMAX + 1)]
    e = [mp.mpf(1)]
    for n in range(1, JMAX + 1):
        e.append(e[-1] * (-x1**2) / n)
    out = {}
    pref = mp.e**(x1**2)
    for k in range(KMAX + 1):
        for l in range(KMAX + 1):
            total = mp.mpc(0)
            for r, coef_fn in Pkl[(k, l)]:
                c = coef_fn(x1)
                if c == 0:
                    continue
                acc = mp.mpc(0)
                for j in range(r, JMAX + 1):
                    acc += e[j - r] / denom[j]
                total += c * acc
            out[(k, l)] = pref * total
    return out

def i_cons(idm, k, l):
    den = 1 - idm[(0, 0)] + 2*idm[(2, 0)] - idm[(2, 2)]
    return idm[(k, l)] + (idm[(k, 0)] - idm[(k, 2)])*(idm[(l, 0)] - idm[(l, 2)])/den

ORACLE = 'verification/oracles/rho_phi_kernels.dat'
rows = []
with open(ORACLE) as fh:
    for line in fh:
        line = line.strip()
        if not line or line.startswith('#'):
            continue
        rows.append(line.split())

COLS = ('id charge lambdaD vT omegaC nu ks kr krp rg mphi A1 A2 kpar omegaE '
        'Vpar omega debye x1 x2 rhoL bplus bcross phaseRe phaseIm I00Re I00Im '
        'I20Re I20Im I01Re I01Im I21Re I21Im F0 F2 rhoPhiDynamicRe '
        'rhoPhiDynamicIm debyeRe debyeIm rhoPhiTotalRe rhoPhiTotalIm '
        'rhoBDynamicRe rhoBDynamicIm').split()
assert all(len(r) == len(COLS) for r in rows) and len(rows) == 13

def g(row, name):
    return mp.mpf(row[COLS.index(name)].replace('E', 'e'))

def gc(row, re_name, im_name):
    return mp.mpc(g(row, re_name), g(row, im_name))

def rel(x, y):
    d = abs(x - y)
    m = max(abs(x), abs(y))
    return d/m if m > 0 else d

report = {t: mp.mpf(0) for t in ('T1', 'T2', 'T3', 'T4', 'T5', 'T5b')}
for row in rows:
    x1, x2 = g(row, 'x1'), g(row, 'x2')
    idm = i_diff_matrix(x1, x2)
    I00, I20 = i_cons(idm, 0, 0), i_cons(idm, 2, 0)
    I01, I21 = i_cons(idm, 0, 1), i_cons(idm, 2, 1)
    I02, I03 = i_cons(idm, 0, 2), i_cons(idm, 0, 3)

    for val, rn, im in ((I00, 'I00Re', 'I00Im'), (I20, 'I20Re', 'I20Im'),
                        (I01, 'I01Re', 'I01Im'), (I21, 'I21Re', 'I21Im')):
        report['T1'] = max(report['T1'], rel(val, gc(row, rn, im)))
    report['T2'] = max(report['T2'], rel(I02, I20))
    report['T3'] = max(report['T3'], rel(I03, I21))

    lam, vT, om_c, nu = g(row, 'lambdaD'), g(row, 'vT'), g(row, 'omegaC'), g(row, 'nu')
    ks, kr, krp, rg = g(row, 'ks'), g(row, 'kr'), g(row, 'krp'), g(row, 'rg')
    mphi = int(row[COLS.index('mphi')])
    A1, A2 = g(row, 'A1'), g(row, 'A2')
    rho = abs(vT/om_c)
    bp = rho**2*(2*ks**2 + kr**2 + krp**2)/2
    bx = rho**2*mp.sqrt(ks**2 + kr**2)*mp.sqrt(ks**2 + krp**2)
    al, alp = mp.atan2(kr, ks), mp.atan2(krp, ks)
    ph = mp.e**(1j*((kr - krp)*rg + mphi*(alp - al)))
    F0 = mp.exp(-bp)*mp.besseli(abs(mphi), bx)
    F2 = mp.exp(-bp)*((1 - bp - mphi)*mp.besseli(abs(mphi), bx)
                      + bx*mp.besseli(abs(mphi - 1), bx))

    dyn_phi = (ph/(4*mp.pi*lam**2)) * 1j*vT**2*ks/(om_c*nu) \
        * (I00*(A1*F0 + A2*F2) + A2*I20*F0/2)
    report['T4'] = max(report['T4'], rel(dyn_phi, gc(row, 'rhoPhiDynamicRe',
                                                     'rhoPhiDynamicIm')))

    rhoB_derived = -(ph/(4*mp.pi*lam**2)) * vT**3/(om_c*nu*SOL) \
        * (I01*(A1*F0 + A2*F2) + A2*I03*F0/2)
    rhoB_transcribed = -(ph/(4*mp.pi*lam**2)) * vT**3/(om_c*nu*SOL) \
        * (I01*(A1*F0 + A2*F2) + A2*I21*F0/2)
    committed = gc(row, 'rhoBDynamicRe', 'rhoBDynamicIm')
    report['T5'] = max(report['T5'], rel(rhoB_derived, committed))
    report['T5b'] = max(report['T5b'], rel(rhoB_derived, rhoB_transcribed))

for key, desc in (('T1', 'replica vs committed I columns'),
                  ('T2', 'symmetry I02=I20'),
                  ('T3', 'identity I03=I21'),
                  ('T4', 'rho-Phi re-assembly vs committed'),
                  ('T5', 'derived rho-B vs committed'),
                  ('T5b', 'derived vs production rho-B form')):
    print(f'{desc:36s} max rel err = {mp.nstr(report[key], 3)}')

FORMULA_INDEX updated accordingly (KIM-04 coverage paragraph). With this, all four blocking findings plus the derivation follow-up from the review are addressed; the remaining known reds are the pre-existing test_rhs_balance, the periodization-deformation gate, and the periodic-vs-global 5% gate, unchanged by this stack.

[marjohma]

Execution caveat resolved — the local Wolfram license issue is fixed (kernel path configure + reactivation), so the official gate has now been executed locally on 314a564:

• make verify-mathematica-charge-kernels: 91 checks, 0 failures — including all 39 new per-case derivation gates (I02=I20, I03=I21, derived rho-B == production form; scaled error 0. throughout) and the 13-row oracle comparison.
• make verify-mathematica-charge-kernels-negative: all six mutations (four-pi, fourier-measure, phase, bessel-index, bplus-sign, hidden-ks2) rejected.
• make regenerate-rho-phi-kernel-oracle: regenerated rho_phi_kernels.dat is byte-identical to the committed file (git diff --quiet clean), confirming no re-bless is needed.

The mpmath cross-check above stands as an independent-engine corroboration; the canonical Mathematica path is now verified as well. Nothing remains open from the review on this PR.