[rogeriojorge]

Stacked on #66 (highk/field-closure, moving-basis stack #56–#60). Stacked PRs get no CI here; the full suite was run locally (197 passed before the TOML test was added, 11/11 in the new file).

What

A physical collision operator, the conservative Dougherty (Lenard–Bernstein with self-consistent moments) operator

C_ss = nu_ss d_v · [(v − U_s) f_s + theta_s d_v f_s], theta_s = tr P_s / (d m_s n_s)

where U_s and theta_s are spatially local moments of the same f_s (Hakim, Francisquez, Juno & Hammett, arXiv:1903.08062, eq. 3–8; d = velocity dimension as in Gkeyll).

Hermite projection (exact, any moving basis)

With xi = (v − u_s)/a_s per axis, delta_i = (U_i − u_i)/a_i, tau_i = theta/a_i², and psi_n = H_n e^{−xi²}/sqrt(pi 2^n n!):
d psi_n = −sqrt(2(n+1)) psi_{n+1}, xi psi_n = sqrt((n+1)/2) psi_{n+1} + sqrt(n/2) psi_{n−1}. Composing the stencils before any truncation gives

(C f)n = nu sum_i [ −n_i C_n + sqrt(2 n_i) delta_i C{n−e_i} + (2 tau_i − 1) sqrt(n_i(n_i−1)) C_{n−2e_i} ]

Every source index is <= n (lower-triangular), so the projection onto the retained modes is exact at any N; no intermediate is truncated and no closure enters. It holds for any (u_s, a_s), including the time-dependent basis of #56–#60.

• Moments use only n <= 2 per axis: n = a_x a_y a_z C_0, U_i = u_i + a_i C_{e_i}/(sqrt2 C_0), P_ii/(mn) = a_i²[1/2 + C_{2e_i}/(sqrt2 C_0) − delta_i²].
• Conservation at finite N: rows 0 and 1 vanish identically; the row-2 energy terms sum to 2 C_0 (d theta − sum_i P_ii/mn) = 0. So density, momentum and energy are conserved to round-off whenever each active axis keeps n <= 2 (N_i >= 3). Axes with N_i <= 2 carry no second moment; on them the projected operator is identically zero, so they are skipped and not counted in d.
• Equilibrium: any Maxwellian with matching U, theta is in the kernel, truncated or not (lower-triangular rows of a kernel element).
• Shifted-centred basis (u = U, a² = 2 theta): every mode decays at exactly n nu.
• Spatial dependence: U(x), theta(x) and the products are formed in real space from the masked state, transformed and multiplied by the same strict 2/3 mask as the rest of the RHS.

Pair collisions (second part of the same commit)

m_s n_s nu_sr = m_r n_r nu_rs (pointwise), u_sr = u_rs = (U_s + U_r)/2, common
T_sr = (m_r T_s + m_s T_r)/(m_s + m_r) + m_s m_r |U_s − U_r|²/(2 d (m_s + m_r)). This is derived here from R_sr + R_rs = 0 and Q_sr + Q_rs = 0 and is Greene's form with beta = 0. Masses are m_s = 1/Omega_cs. Note: arXiv:1903.08062 leaves the cross-species coefficients to later work, and the Gkeyll collision-model docs page could not be reached, so the beta convention was not checked against Francisquez et al.; the conservation itself is tested.

Interface

• simulation(..., collision_model="dougherty") (static) with input_parameters["nu_dougherty"]: a scalar (self collisions of every species) or an (Ns, Ns) matrix (diagonal: nu_ss; upper triangle [s, r]: nu_sr; the lower triangle is ignored, nu_rs follows from momentum conservation).
• TOML: [collisions] with model = "dougherty" and nu = ....
• Separate from the numerical hypercollision nu/collision_matrix. Default collision_model=None: the RHS is unchanged (no term is added; bitwise-equality test).

Tests (tests/test_dougherty.py, CPU, ~17 s)

• projection of the operator on an anisotropic, drifting, skewed 2V distribution against direct evaluation (autodiff of the flux, quadrature projection): agreement to 1e-9 relative; moments against direct quadrature
• density / momentum / energy conservation to round-off at finite N, per x, with spatially varying moments, Ns = 1 and Ns = 3 with pair collisions
• R_sr + R_rs = 0, Q_sr + Q_rs = 0 for the pair moments
• drifting Maxwellian narrower than and off the basis: exact equilibrium (< 1e-13); anisotropic Maxwellian is not
• n·nu decay rates in the shifted-centred basis
• time-integrated: T_x − T_y decays like exp(−2 nu t) (d = 2), 1e-7; two-species drift difference decays like exp(−(nu_sr + nu_rs) t / 2) with total momentum and energy conserved to 1e-11
• off is bitwise unchanged; nu = 0 matches
• jax.grad through a collisional run vs central differences
• TOML [collisions] table

Cost per RHS (laptop CPU, 4 threads, Ns = 2 with pair collisions)

grid	off	dougherty	ratio
Nx=129, N=(64,1,1)	0.386 ms	0.585 ms	1.51
Nx=33, N=(16,16,1)	0.378 ms	0.495 ms	1.31
Nx=17, N=(8,8,8)	0.416 ms	0.603 ms	1.45

The extra cost is one rFFT of the coefficient array plus the moment products.

Limits

• nu is constant per pair (no n- or T-dependence of nu yet).
• U and theta divide by the local C_0, so a state whose density approaches zero somewhere is not supported.
• Explicit stepping: the stiffness is about N nu.