Struct SpectralV 

pub struct SpectralV {
    pub coefficients: Vec<f64>,
    pub spatial_shape: [usize; 3],
    pub n_modes: usize,
    pub velocity_scale: f64,
    pub domain: Domain,
}

Spectral-in-velocity representation using orthonormal Hermite function basis.

The distribution function is represented as: f(x,v,t) = sum_{n1,n2,n3} a_{n1,n2,n3}(x,t) * psi_{n1}(v1/sigma) * psi_{n2}(v2/sigma) * psi_{n3}(v3/sigma)

where psi_n(u) = H_n(u) * exp(-u^2/2) / sqrt(2^n * n! * sqrt(pi)) are the orthonormal Hermite functions (eigenfunctions of the Fourier transform and harmonic oscillator).

Layout

Coefficients are stored as a flat array. For each spatial cell (ix, iy, iz), there are n_modes^3 coefficients. The flat index is: spatial_flat * n_modes^3 + m0 * n_modes^2 + m1 * n_modes + m2

where spatial_flat = ix * ny * nz + iy * nz + iz.

Fields

coefficients: Vec<f64>

Hermite expansion coefficients: length = n_spatial * n_modes^3.

spatial_shape: [usize; 3]

Spatial grid dimensions [nx, ny, nz].

n_modes: usize

Number of Hermite modes per velocity dimension.

velocity_scale: f64

Velocity scaling parameter sigma (controls Hermite basis width). Chosen so that the Gaussian weight exp(-v^2/(2*sigma^2)) covers the velocity domain.

domain: Domain

Domain specification.

Implementations

impl SpectralV

pub fn new(domain: Domain, n_modes: usize) -> Self

Create a new SpectralV with zero coefficients.

The velocity scale sigma is set to Lv / 3 so that 3-sigma covers the velocity half-extent, placing most of the Gaussian weight inside the domain.

pub fn hermite_raw(n: usize, x: f64) -> f64

Evaluate the raw (unnormalized) physicist’s Hermite polynomial H_n(x).

Uses the three-term recurrence: H_0(x) = 1 H_1(x) = 2x H_n(x) = 2x * H_{n-1}(x) - 2(n-1) * H_{n-2}(x)

pub fn from_snapshot( snap: &PhaseSpaceSnapshot, n_modes: usize, domain: &Domain, ) -> Self

Project a 6D snapshot onto the Hermite basis.

For each spatial cell (ix, iy, iz), the coefficients are: a_{n1,n2,n3}(x) = integral f(x,v) * psi_{n1}(v1/sigma) * psi_{n2}(v2/sigma) * psi_{n3}(v3/sigma) dv^3 / sigma^3

This integral is approximated by quadrature over the velocity grid in the snapshot.

Trait Implementations

impl PhaseSpaceRepr for SpectralV

fn compute_density(&self) -> DensityField

Compute density rho(x) = integral f(x,v) dv^3.

In the Hermite representation: rho(x) = a_{0,0,0}(x) * sigma^3 * (integral psi_0(u) du)^3

Only the zeroth mode contributes because all higher modes integrate to zero (they are orthogonal to the constant = psi_0 * norm_0, and integral H_n * w du = 0 for n > 0).

fn advect_x(&mut self, _displacement: &DisplacementField, dt: f64)

Spatial drift: advect each coefficient grid a_{n1,n2,n3}(x) by displacement v*dt.

In the Hermite representation, the spatial advection operates on each coefficient grid independently: the Hermite modes do not couple during spatial transport. Each a_{n1,n2,n3}(x) grid is shifted using semi-Lagrangian interpolation.

For a given mode, the “velocity” associated with the drift is the mean velocity of that Hermite function. However, since the Vlasov spatial drift is v * df/dx, and v is not a single value but a distribution, the correct approach in SpectralV is to use the recurrence relation: d/dt a_n = -sigma * (sqrt((n+1)/2) * da_{n+1}/dx + sqrt(n/2) * da_{n-1}/dx)

For implementation simplicity and because this avoids CFL constraints, we use the mode-coupling approach with finite differences for the spatial gradient, but applied dimension-by-dimension.

Actually, the cleaner approach: since displacement already encodes v*dt per spatial cell, we can treat each mode’s coefficient grid as being advected by the mode-weighted displacement. For zeroth mode (Gaussian), the effective velocity is zero; for first mode it is proportional to sigma, etc.

Following the standard SpectralV literature (Schumer & Holloway 1998), we implement the x-advection as mode coupling: The spatial derivative couples modes: v * psi_n(v/sigma) = sigma * (sqrt((n+1)/2) * psi_{n+1} + sqrt(n/2) * psi_{n-1})

For simplicity and correctness with the existing semi-Lagrangian infrastructure, we shift each coefficient grid by a uniform displacement derived from the displacement field (which encodes v*dt averaged over the domain). This is exact for uniform advection; for non-uniform v fields, we use the displacement field values directly.

fn advect_v(&mut self, acceleration: &AccelerationField, dt: f64)

Velocity kick: acceleration couples adjacent Hermite modes.

The recurrence relation for Hermite functions under velocity translation gives: da_n/dt = (g/sigma) * [sqrt(n/2) * a_{n-1} - sqrt((n+1)/2) * a_{n+1}]

This is applied per velocity dimension, per spatial cell. The coupling is tridiagonal in mode space and exact for constant acceleration within a cell. We use forward Euler time integration of the mode coupling.

fn moment(&self, position: &[f64; 3], order: usize) -> Tensor

Compute velocity moment at a given spatial position.

• Order 0: density rho = a_{0,0,0} * normalization
• Order 1: mean velocity <v_i> = a_{e_i} * sigma * sqrt(2) * norm / rho where e_i is the unit vector with 1 in dimension i (mode (1,0,0), (0,1,0), or (0,0,1))
• Order 2: velocity dispersion tensor (reconstructed from coefficients)

fn total_mass(&self) -> f64

Total mass M = integral f dx^3 dv^3.

Since only the zeroth Hermite mode contributes to the velocity integral, M = sum_x a_{0,0,0}(x) * norm * dx^3.

fn casimir_c2(&self) -> f64

Casimir invariant C_2 = integral f^2 dx^3 dv^3.

By Parseval’s theorem for the Hermite expansion: C_2 = sum_x sum_{n1,n2,n3} |a_{n1,n2,n3}(x)|^2 * dx^3 * sigma^3

The sigma^3 factor comes from the change of variable v -> u = v/sigma and the fact that the psi_n are orthonormal in u-space.

fn entropy(&self) -> f64

Entropy S = -integral f * ln(f) dx^3 dv^3.

No closed form exists in Hermite space; we reconstruct f on a velocity quadrature grid and integrate numerically.

fn stream_count(&self) -> StreamCountField

Stream count: number of peaks in the marginal f(v_x | x) at each spatial cell.

Reconstructs f on a 1D velocity grid (marginalised over v_y, v_z) and counts peaks.

fn velocity_distribution(&self, position: &[f64; 3]) -> Vec<f64>

Extract the local velocity distribution f(v | x) at a given spatial position.

Returns f on the velocity grid defined by the domain’s velocity resolution.

fn total_kinetic_energy(&self) -> f64

Total kinetic energy T = (1/2) integral f * v^2 dx^3 dv^3.

Using the Hermite recurrence for v^2: v_d^2 = sigma^2 * (1/2 + u_d^2) u^2 * psi_n(u) = sqrt(n(n-1)/4) * psi_{n-2} + (n + 1/2) * psi_n + sqrt((n+1)(n+2)/4) * psi_{n+2}

So integral v_d^2 * psi_{m_d}(u_d) * psi_0(u_other) du^3 is nonzero only for specific mode indices. For simplicity, we compute via quadrature.

fn to_snapshot(&self, time: f64) -> PhaseSpaceSnapshot

Reconstruct the full 6D phase-space snapshot from Hermite coefficients.

fn as_any(&self) -> &dyn Any

Downcast to concrete type for implementation-specific queries (e.g. HT rank data).

fn load_snapshot(&mut self, snap: PhaseSpaceSnapshot)

Replace the current state with data from a dense 6D snapshot.