feat: super-Gaussian-preserving Fokker-Planck operator for Vlasov-1D

[joglekara]

Summary

Adds a collision operator for vlasov1d whose equilibrium is a super-Gaussian $f_0 \propto \exp(-\beta|v-\bar v|^m)$ instead of a Maxwellian:

terms:
  fokker_planck:
    is_on: True
    type: super_gaussian
    m: 3.0
    self_consistent_beta: {enabled: True, max_steps: 3}   # recommended for long runs

This is the 1D-1V counterpart of maintaining a Langdon/DLM inverse-bremsstrahlung-heated distribution: the standard LB/Dougherty operators actively maxwellianize a super-Gaussian background, which is unwanted in Landau-damping / LPI studies with DLM distributions. This operator maintains a prescribed order m (it does not generate m dynamically — no heating/e-e competition as in the VFP-1D IB operator; a prescribed m(x,t) from Matte's fit is a possible follow-up).

Design

• Exact-potential drift: C_edge = D·Δφ/dv with φ = β|v-v̄|^m, rather than the midpoint derivative. Chang-Cooper's fixed point satisfies f_{i+1}/f_i = exp(-C·dv/D), so this makes the sampled super-Gaussian the exact discrete equilibrium. (The midpoint form is exact only at m=2 and causes secular temperature drift otherwise.)
• Energy-conserving β closure: β = n/(m·⟨|v-v̄|^m⟩) — the generalization of the LB/Dougherty β = 1/(2T), which it equals at m=2 — refined by a Newton solve of the δ-weighted discrete energy-flux condition Σₑ vₑ(w·f̃(w) + Δf) = 0. A sampled super-Gaussian is an exact root of this condition, so the discrete equilibrium is exactly stationary (no secular drift, verified < 1e-8 over 100 collision times). The Newton reuses the existing self_consistent_beta config knobs.
• m=2 reduces to chang_cooper_dougherty step for step (tested).
• find_self_consistent_beta is generalized to super-Gaussian targets (m parameter, default 2.0 preserves existing behavior bitwise).

Conservation properties

• Density: exact (zero-flux BC).
• Energy: exact at equilibrium (with self_consistent_beta enabled); off-equilibrium transients carry a one-time O(ν·dt) offset from operator splitting (β frozen from fⁿ during the implicit step) — first-order convergence in dt verified in tests; negligible at production ν·dt.
• Momentum: exact for f symmetric about v̄; skewed transients exchange momentum at O(skewness) (the m=2 C ∝ v-v̄ cancellation doesn't generalize). Documented in the config reference.

Testing

• 8 new tests in tests/test_vlasov1d/test_super_gaussian_fp.py: SG fixed point (m=3, 4), no secular drift at equilibrium, Maxwellian→SG relaxation with O(dt) energy-error convergence, Dougherty control (maxwellianizes the same IC), shifted-SG momentum conservation, m=2 ≡ Dougherty, self-consistent-β unit test.
• Existing suites pass: vlasov1d FP relaxation/conservation/momentum + config validation (24), vfp1d FP model tests (10).
• End-to-end smoke via ergoExo on the FP-conservation config with type: super_gaussian: kurtosis 3.00 → 2.44 toward the m=3 target 2.418, density conserved to 1e-14.

Note: test_config_regression.py has 3 pre-existing failures on main (stale baselines from the vmin feature), unrelated to this PR and being fixed separately.

🤖 Generated with Claude Code

[ManharGupta25]

Interesting. It very much seems to be in the same spirit as configuring the Krook operator to relax to the initialized super-Gaussian. But this is better physics-grounded since we have the conservation properties that Jon was talking about for vfp1d.

I am wondering for SRS sims, if the momentum leak would matter since we usually have the hot electron tail. I am not sure about this, but would this operator also end up eroding that hot electron tail to maintain the symmetric super-Gaussian?

Just asking, even though the hot tail stuff doesn't matter as much for studying threshold shifts

[joglekara]

I think it will be a minor loss of conservation so it should probably just be negligible... we'll find out!

[ManharGupta25]

I guess we'll see.

[ManharGupta25]

Just put the PR through claude. Its pointing out if the gradients through this operator have been tested to work or not?

Claude's comment with some edits:

The design is clean — the exact-Δφ Chang–Cooper construction making the sampled super-Gaussian the exact discrete fixed point is nicely done, and the conservation limitations (O(ν·dt) energy transient, momentum only for symmetric f) are documented up front rather than buried. m=2 ≡ Dougherty being tested bitwise is good insurance against regressions.

Worth confirming before merge: since this operator will likely end up inside optimization loops, is the gradient through it exercised anywhere? It now carries an optx.root_find for β in the collision step — implicit diff should handle that fine, but a quick jax.grad sanity check (finite + finite-difference-consistent) would close it, or an ack that differentiability is tracked separately.

[joglekara]

im not sure but I am a bit scared of having that inside the solve. it'd be best to avoid having to do a nonlinear solve like that, it might be very expensive.

[ManharGupta25]

Updated comment as a whole from a previous, now deleted, PR comment:

Just verified the non-linear solve and the effect of removing it with Claude. I guess we can still have the non-linear solve since its optional and only required at certain cases (very long, high nu_fp runs). For our runs, having it off wouldn't have much of an effect since the secular temperature/variance drift its trying to remove is already small (more in Claude's comment).

Claude's comment with some edits:
TL;DR: the optx.root_find for β is optional (self_consistent_beta) and unnecessary for our runs — turn it off and both the cost worry and the gradient-through-root_find worry go away, with no loss to what we actually need.

Why the solve is optional — the connection: the Newton solve exists for one purpose: it refines β so the sampled super-Gaussian is the exact discrete fixed point, removing a small secular temperature drift (a slow monotonic creep in ⟨(v−v̄)²⟩ from the O(dv²) error in the cheap closure β = n/(m·⟨|v−v̄|ᵐ⟩)).

So "do we need the solve?" collapses to "is that drift big enough to matter for us?" — and it isn't:

run	coll-times (ν·t)	drift @ nv=128 (~1e-5/ct)	drift @ nv=256 (~2.5e-6/ct)
25 ps	1.34	~1.3e-5	~3.4e-6
40 ps	2.15	~2.2e-5	~5.4e-6
60 ps	3.22	~3.2e-5	~8.1e-6

(ν_fp = 1e-5.) The drift stays ~1e-5–3e-5 across our full run range — orders of magnitude below anything that moves a threshold; it wouldn't reach ~1e-3 until ~1.9 ns (nv=128).

Key point: turning the solve off costs only that tiny drift — it does not cost shape maintenance. The super-Gaussian order m is held by the exact-Δφ Chang–Cooper drift term, which is a separate mechanism from the β refinement. So self_consistent_beta: off still prevents the erosion we care about.

---

Update (re: the max_steps: 1 commit): the corrected measured drift (~5e-7/ct at nv=128, m=3) makes the numbers above ~20× smaller — over our runs it's ~7e-7 (25 ps) to ~1.6e-6 (60 ps), even more negligible. So the takeaway holds either way: a single Newton step (max_steps: 1) gets machine-level at trivial cost, and off is also fine on drift grounds. A single step isn't the iterated nonlinear solve that prompted the cost worry, so this basically resolves it. The one place off still has a distinct edge is inside an optimization loop (keeps root_find out of the gradient graph) — not our current threshold-scan campaigns. Net: max_steps: 1 for the scans, and the gradient-through-root_find check only matters if/when this operator goes in an opt loop.