Mixed Hermite-Legendre 1D Configuration Reference

This document describes how to construct a configuration file for the hermite-legendre-1d solver, which implements the mixed Hermite-Legendre spectral method for the 1D-1V electrostatic Vlasov-Poisson system (Issan, Delzanno & Roytershteyn, arXiv:2606.12322).

The electron distribution is split f = f0 + df:

• f0 (near-Maxwellian bulk) is expanded in the asymmetrically-weighted (AW) Hermite basis in velocity, with coefficients C_n(x, t), n = 0 .. Nh-1.

• df (strongly non-Maxwellian features: beams, plateaus, filamentation) is expanded in the Legendre basis on a bounded velocity window [v_a, v_b], with coefficients B_m(x, t), m = 0 .. Nl-1.

The highest Hermite coefficient C_{Nh-1} feeds the Legendre modes (one-way coupling), and both feed the self-consistent field through Poisson. The method is most accurate, at fixed total velocity DOFs, when non-Maxwellian features are localized in velocity.

Normalization (paper sec 2.1): time by 1/ω_pe, space by the Debye length λ_D, velocity by the electron thermal velocity v_the. A single electron species is evolved against an immobile neutralizing ion background of density 1.

Numerics. Space is treated spectrally (Fourier, periodic domain); both free-streaming operators are symmetric-tridiagonal in mode index and integrated exactly via prediagonalized matrix exponentials. The E-field force, the Legendre Dirichlet penalty, and the Hermite→Legendre coupling are advanced explicitly with Lawson-RK4. (The paper uses an implicit-midpoint integrator for machine-precision energy conservation; this module uses an explicit integrator — energy is then conserved to the time-integrator’s order, which converges with dt, while mass and momentum remain conserved to machine precision.)

Top-Level Structure

solver: hermite-legendre-1d
mlflow: ...
units: ...
physics: ...
grid: ...
initialization: ...
save: ...

---

physics

Field	Type	Default	Description
Lx	float	—	Domain length in x (normalized to λ_D)
alpha	float	—	AW-Hermite velocity scale parameter α (the benchmarks use √2)
u	float	0.0	AW-Hermite velocity shift parameter u
v_a, v_b	float	—	Legendre velocity-window bounds (df is resolved on [v_a, v_b])
gamma	float	0.5	Penalty coefficient γ for the weak Legendre Dirichlet BC (df(v_a)=df(v_b)=0). Applied only to modes m ≥ 3 to preserve conservation.
nu_H	float	0.0	Artificial (Lenard-Bernstein) Hermite collision rate ν_H. Keep small/zero so f0 can feed df through the last Hermite moment.
nu_L	float	0.0	Artificial Legendre collision rate ν_L. Controls filamentation/recurrence in df.
enforce_conservation	bool	true	Zero the coupling integrals J_{Nh,0}=J_{Nh,1}=J_{Nh,2}=0 so the discrete method conserves mass, momentum, and energy independent of Nh parity and of α, u (paper sec 3.4/4).
field	bool	true	Self-consistent Poisson field. Set false for the pure linear-advection test (φ = 0); the linear Hermite→Legendre closure flux still acts.

The artificial collision operator (paper sec 2.5) uses the cubic spectrum col[n] = n(n-1)(n-2) / ((N-1)(N-2)(N-3)), which is identically zero for n = 0, 1, 2 — so collisions never touch the mass/momentum/energy moments.

Choosing Nh (important). Keep the Hermite basis bulk-only — structure inside the Legendre window belongs to df. At large Nh (≳64) and saturation-scale fields, the nonlinear force ladder (pump ~|E|·√(2n)/α) outruns the cubic collision damping in the mid-n window and a spurious k=0 velocity-space cascade grows there at ~10× the physical rate, eventually destroying the run — at any practical ν_H (0, 10, and 30 were all measured to fail on the bump-on-tail benchmark). A small basis closes the window structurally: bump-on-tail with Nh=32 reproduces the Nh=128 field observables to 3 digits with machine-precision energy conservation.

---

grid

Field	Type	Default	Description
Nx	int	—	Number of Fourier modes in x
Nh	int	—	Number of AW-Hermite modes for f0 (closure by truncation: C_{Nh}=0)
Nl	int	—	Number of Legendre modes for df (closure by truncation: B_{Nl}=0)
tmax	float	—	Final simulation time (normalized). Snapped to an exact multiple of dt.
dt	float	0.01	Timestep
integrator	str	"lawson"	Time integrator: "lawson" (explicit Lawson-RK4), "imex" (Lawson-RK4 + implicit Lorentz substep), or "implicit" (implicit midpoint, AD-JFNK) — see below.
newton_iters	int	3	(implicit) Newton iterations per step.
gmres_restart, gmres_maxiter, gmres_tol	int/int/float	20/4/1e-8	(implicit) matrix-free GMRES controls for the Newton linear solves.
precondition	bool	true	(implicit) use the streaming+collision operator as a physics-based GMRES preconditioner (see below).

integrator: implicit (implicit midpoint via AD-JFNK)

Advances the full RHS with the implicit-midpoint rule y1 = y0 + dt·F((y0+y1)/2), solved by Jacobian-free Newton-Krylov: each Newton linear system uses a matrix-free GMRES whose Jacobian-vector products are exact autodiff JVPs (jax.linearize) — the Jacobian is never assembled (memory is the state plus a few Krylov vectors). Implicit midpoint is A-stable (no CFL at all) and conserves quadratic invariants, so it conserves mass exactly and energy to the solve tolerance, and stays stable into the saturated / long-time regime where both lawson and imex blow up (e.g. bump-on-tail). Cost: each step does newton_iters × (GMRES iterations) × (RHS evals), so it is the most expensive per step — use it for the hard cases, not the cheap ones.

Preconditioning (precondition: true, default). The implicit operator’s stiffness is dominated by the skew streaming term, whose eigenvalues smear along the imaginary axis (~dt/2·α·k_max·√(2Nh)) — the worst case for unpreconditioned GMRES, which then needs many iterations and can fail to converge at large dt/Nx (Newton then injects energy). The preconditioner M = I − dt/2·(L_streaming + L_collision) is block-diagonal in k and tridiagonal in mode index, so M⁻¹ is a cheap per-k tridiagonal solve that captures exactly that stiff spectrum; GMRES on M⁻¹A then converges in a handful of iterations.

Two measured caveats. (1) Do not precondition the Lorentz-force block: its triangular factors are strongly non-normal (nilpotent ladders, norm ~Nl²/width·|E|), and applying their inverse stalls GMRES entirely at saturation-scale fields. (2) At large dt (≳0.05) combined with large |E|, the stream-preconditioned solve itself degrades while unpreconditioned GMRES still converges — set precondition: false for large-dt experiments in strongly nonlinear regimes.

integrator: imex

The stiffness that limits the explicit step is the E·∂_v f Lorentz force: in the spectral velocity bases it is strictly lower-triangular (nilpotent for Hermite, lower-triangular + a rank-2 penalty for Legendre) with operator norm ~Nl²/width·|E| — so explicit RK4’s |dt·‖L‖|≲2.8 limit tightens as modes/field grow. Setting integrator: imex keeps free-streaming, collisions, and the Hermite→Legendre closure flux in the explicit Lawson step, and advances the Lorentz force with an unconditionally stable frozen-E Backward-Euler substep (a per-x triangular/dense linear solve; first-order Lie split). This removes the CFL limit, letting two-stream run at dt ≈ 0.02 instead of 0.002. Trade-offs: Backward Euler is mildly dissipative and the split is first-order in dt, so for high-accuracy/conservation studies prefer small-dt lawson; for robustness at large mode counts or large Nx, prefer imex.

Choosing dt. Free-streaming and collisions are integrated exactly, but the explicit Lawson-RK4 treatment of the E-field force has a stability (CFL) limit that tightens as the self-consistent field grows. For small-amplitude/linear runs (e.g. driven Landau damping) dt = 0.05 is fine; for nonlinear instabilities that saturate to a large field (two-stream) a smaller step is needed — dt ≈ 0.002 is stable and converged for the two-stream benchmark. (The paper’s dt = 0.01 relies on its unconditionally stable implicit-midpoint integrator; this explicit module trades that for a smaller step and a much smaller memory footprint.) A run that goes NaN partway through is the signature of dt above the CFL limit — halve it.

---

initialization

Selects how the initial C_n(x) and B_m(x) coefficients are built.

type	Parameters	Description
linear-advection	eps, mode	f0 = (1 + eps·cos(k x))/√(2π)·exp(-v²/2); df = 0. (C_0 = n(x)/α.)
two-stream	eps, mode	f0 ∝ (1 + eps·cos(k x))·v²·exp(-v²/2): C_0 = n(x)/α, C_2 = √2·C_0; df = 0.
bump-on-tail	eps, mode, n_beam, v_drift, v_th	Bulk Maxwellian in f0; a drifting Gaussian beam n_beam/(√(2π) v_th)·exp(-(v-v_drift)²/2v_th²) projected onto Legendre as df.
custom	hermite: {n: {base, eps, mode}}, df: {beams: [{amp, v_drift, v_th}], eps, mode}	Generic Hermite coefficient profiles plus a beam/sum-of-Gaussians df projected onto Legendre.

Here k = 2π·mode/Lx. The Legendre projection uses Gauss-Legendre quadrature.

---

drivers (optional)

An external longitudinal field ex can be applied to the velocity-space force (it never enters the Poisson solve), e.g. to drive a resonant EPW for a Landau-damping measurement — the analogue of the Vlasov-1D ex driver. Omit the drivers block for self-consistent runs.

drivers:
  ex:
    '0':                 # one entry per pulse
      k0: 0.4            # wavenumber
      w0: 1.285          # angular frequency (e.g. Re(omega) from the dispersion relation)
      dw0: 0.0           # frequency offset (added to w0)
      a0: 1.0e-3         # amplitude
      t_center: 20.0     # pulse: center / full width / rise(+fall) time
      t_width: 20.0
      t_rise: 5.0
      x_center: 7.85     # spatial envelope: center / width / rise (defaults span the box)
      x_width: 1.0e6
      x_rise: 1.0

The driver field is E_drive(x,t) = Σ env(x,t)·(w0+dw0)·a0·sin(k0 x − (w0+dw0) t) and is saved as de in the fields group.

---

save

Standard ADEPT save block with t: {nt: ...} (or tmin/tmax/nt) sub-axes.

Key	Contents
fields	Electric field e(x,t), potential phi(x,t), and external driver field de(x,t)
hermite	AW-Hermite-Fourier coefficient timeseries Ck (shape nt × Nh × Nx)
legendre	Legendre-Fourier coefficient timeseries Bk (shape nt × Nl × Nx)
default	Scalar invariants mass, momentum, energy (paper eqns 26, 28, 30-31) plus field energy and density extrema. Always added; the primary correctness gate.

post_process writes netCDF binaries and spacetime/scalar plots, and reports the relative drift of each invariant as the metrics reldrift_{mass,momentum,energy}.

---

Example: two-stream instability

solver: hermite-legendre-1d
mlflow: {experiment: hermite-legendre-1d, run: two-stream}
units: {normalizing_density: 1e20/cc, normalizing_temperature: 1keV}
physics:
  Lx: 12.566370614359172   # 4π
  alpha: 1.4142135623730951
  u: 0.0
  v_a: -2.5
  v_b: 2.5
  gamma: 0.5
  nu_H: 0.0
  nu_L: 1.0
  enforce_conservation: true
  field: true
grid: {Nx: 64, Nh: 85, Nl: 171, tmax: 35.0, dt: 0.01}
initialization: {type: two-stream, eps: 0.01, mode: 1}
save:
  fields: {t: {nt: 351}}
  legendre: {t: {nt: 71}}

See configs/hermite-legendre-1d/ for the linear-advection, two-stream, and bump-on-tail benchmark configurations.