Hermite-Poisson: AW skew-symmetric closure (Dai & Després 2405.07811) for quadratic stability

[joglekara]

Summary

Adopt the structure-preserving skew-symmetric closure of Dai & Després (arXiv:2405.07811, "On the quadratic stability of asymmetric Hermite basis with application to plasma physics with oscillating electric field") for the asymmetrically-weighted (AW) Hermite acceleration operator in _hermite_poisson_1d (and, by extension, the shared _spectrax1d coupling).

We commit to the AW (Parker–Dellar) convention rather than switching to symmetrically-weighted (SW) Hermite, because (a) it keeps us aligned with the Hermite-spectral plasma community whose work we want to connect to, and (b) AW retains the fluid-moment interpretability (mass/momentum/energy as exact low-order moments).

Background: the instability mechanism

The paper analyses exactly our setting. In the AW basis the acceleration operator e ∂_v is lower-bidiagonal — a single index-raising shift C_{n-1} → C_n — and the natural conserved quantity is the Gram-weighted norm ⟨C, A C⟩, not the Euclidean Σ|C_n|². The operator is skew-symmetric w.r.t. the Gram matrix A in the infinite system, but truncation at the top mode N drops the C_N → C_{N+1} coupling and destroys that skew-symmetry, allowing the quadratic norm to grow without bound — even for a constant E and with zero time-stepping error.

This matches our code exactly. In _hermite_e_coupling the term is

dCk_n/dt|E = (q/m)·√(2n)/α · F(x) · C[n-1](x)

i.e. the paper's operator D with d_{mn} = -e√(2m/T) δ_{m-1,n}, a one-sided shift (C_up[n] = C[n-1]) with no closure at the top mode. Free-streaming (FreeStreamingExp1D) is built symmetric and exponentiated unitarily, so it stays bounded; the unbounded-growth channel is concentrated in the explicit RK4 acceleration term — precisely the operator the paper isolates.

This is a strong candidate root cause for #262. The symptoms logged there and in the spectrax instability debugging are consistent with a spatial-truncation effect rather than a time-integration one:

• refining dt doesn't help (onset comes earlier), and
• adding a Hou–Li filter made things worse (it breaks the conservation/energy balance) — exactly the failure mode the paper warns about when stabilising with ad-hoc dissipation.

The fix (Eq. 29)

A cheap closure on the highest Hermite mode that restores exact skew-symmetry:

D̄^{11,N} = D^{11,N} + (A^{11,N})^{-1} A^{12,N} D^{21,N}

Computed once at setup by solving a single small linear system against the (constant) Gram matrix; it modifies only the last row/column of the acceleration operator. Fully JIT-compatible and differentiable.

The Gram matrix is built from the stable recurrences in their Thm 4.4:

• a_00 = (2T)^{-1/2}, a_{m+1,m+1} = (2m+1)/(2m+2) · a_{mm}
• off-diagonal: a_{m-l-1,m+l+1} = √((m-l)/(m+l+1)) · a_{m-l,m+l}, with a_{mn}=0 for m+n odd.

The paper explicitly prefers this over a Fokker–Planck/collision perturbation or a dynamically reweighted norm, because it is structure-preserving — which is what we want for a differentiable solver where artificial filters inject nonphysical gradients.

Implementation plan

• Build the AW Gram matrix A (Thm 4.4 recurrences) at module setup in _hermite_poisson_1d.
• Apply the Eq. 29 top-mode correction inside _hermite_e_coupling (precompute z from A·z = first column of A^{12,N} once; apply as a fixed correction each step).
• Add a Gram-weighted energy diagnostic ⟨C, A C⟩ to _hermite_poisson_1d/storage.py so we can measure norm conservation.
• Once validated, dial the hyper-diffusion / hypercollision / Hou–Li knobs in DiagonalExp1D back to physical-only levels (or off) and confirm stability holds without them.
• Mirror the fix into the shared _spectrax1d coupling (same operator shape) — likely resolves Hermite solver has trouble with the nonlinear term #262.
• Update docs/source/solvers/hermite_poisson_1d/config.md if any knobs change.

Validation

• Constant-E / free-streaming limit: confirm ⟨C, A C⟩ is conserved to round-off with the closure on and all dissipation off (this directly confirms the mechanism).
• Landau damping (tests/test_hermite_poisson_1d/test_landau_damping.py): rate and frequency unchanged with dissipation off.
• Strongly driven EPW (the Hermite solver has trouble with the nonlinear term #262 MWE): check whether the closure alone removes the blowup; reference the successful adaptive-DoPri8 runs (mlflow experiment 188596) for ground truth, and post the new experiment number.

Caveat

The paper's stability proof is for constant E only — despite the title, time-dependent/oscillating E(t) is not rigorously covered, and our driven (EPW/SRS) cases are exactly that regime. So the closure removes one rigorous source of spurious growth but is necessary-not-sufficient for the driven runs; we still validate empirically on the EPW MWE.

References

• Dai & Després, arXiv:2405.07811 (Eq. 9, Lemma 3.5, Eqs. 29/32, Thm 4.1/4.4, Lemmas 5.1–5.3).
• Parker & Dellar (2015), AW Hermite force term (eq. 3.11) — already cited in _hermite_e_coupling.
• Related: Hermite solver has trouble with the nonlinear term #262 (Hermite solver nonlinear-term instability).