Vlasov-1D Configuration Reference

This document describes how to construct a configuration file for the vlasov-1d solver.

Top-Level Structure

solver: vlasov-1d

units:
  # Physical unit normalizations

density:
  # Species definitions

grid:
  # Simulation grid parameters

save:
  # Output configuration

mlflow:
  # Experiment tracking

drivers:
  # External drivers

diagnostics:
  # Diagnostic outputs

terms:
  # Solver configuration

units

Physical unit normalizations for the simulation.

Field	Type	Description
normalizing_temperature	string	Reference temperature with unit, e.g., "2000eV"
normalizing_density	string	Reference density with unit, e.g., "1.5e21/cc"
reference	string	Reference species for the normalization: "electron" (default) or "ion"
A	float	Ion mass number, \(m_i = A\,m_p\) (only used with reference: ion; default 1.0)
Z	float	Ion charge state, \(q_i = Z e\) (only used with reference: ion; default 1.0)

Example:

units:
  normalizing_temperature: 2000eV
  normalizing_density: 1.5e21/cc

Normalization convention

All Vlasov-1D quantities are normalized using the following unit system, built from normalizing_density (\(n_0\)) and normalizing_temperature (\(T_0\)):

Unit	Definition	Meaning
time	\(\tau = 1/\omega_{p0}\), \(\omega_{p0} = \sqrt{n_0 e^2/(\epsilon_0 m_e)}\)	inverse plasma frequency
velocity	\(v_0 = \sqrt{T_0/m_e}\)	electron thermal speed (RMS / standard-deviation convention)
length	\(L_0 = v_0/\omega_{p0} = \lambda_{De}\)	electron Debye length
wavenumber	\(1/L_0\)	a code wavenumber is \(k \lambda_{De}\)

Consequences of the \(v_0 = \sqrt{T_0/m_e}\) (σ) convention:

• A Maxwellian at code temperature \(T = 1\) is \(f \propto e^{-v^2/2}\) (unit variance).

• The Bohm–Gross dispersion in code units is \(\omega^2 = 1 + 3 k^2\).

• The normalized speed of light is \(\hat c = c/\sqrt{T_0/m_e}\) (e.g. \(\hat c = 15.98\) at 2000 eV).

Dimensional inputs (strings with units, e.g. xmax: 100um) are converted with these units; plain numeric inputs are taken to already be in code units and pass through unchanged.

Ion reference (reference: ion)

For runs where only ions are evolved kinetically (e.g. terms.field: poisson-boltzmann), set reference: ion to build the normalization from the ion species instead:

Unit	Definition	Meaning
time	\(\tau = 1/\omega_{pi}\), \(\omega_{pi} = \sqrt{n_i Z^2 e^2/(\epsilon_0 m_i)}\)	inverse ion plasma frequency
velocity	\(v_0 = \sqrt{T_0/m_i}\)	ion thermal speed (same σ convention)
length	\(L_0 = v_0/\omega_{pi}\)	ion Debye length

with \(m_i = A\,m_p\) and charge \(Z e\); normalizing_density and normalizing_temperature are then the ion density \(n_i\) and ion temperature \(T_i\). Dimensional string inputs (e.g. xmax: 300um, tmax: 50ps) convert with these ion units, and the physical quantities logged by write_units() (wp0, tp0, x0, v0, …) refer to the ion reference. The electron-specific collision entries (nuee, nuee_norm, logLambda_ee) are replaced by their ion-ion counterparts (nuii, nuii_norm, logLambda_ii, from the NRL formulary).

Example (deuterium):

units:
  normalizing_temperature: 100eV   # T_i
  normalizing_density: 1.0e20/cc   # n_i
  reference: ion
  A: 2.0
  Z: 1.0

density

Species and density configuration. You can define multiple “species” by using keys prefixed with species-.

Field	Type	Description
quasineutrality	bool	Whether to enforce quasineutrality

Species Definition

Each species is defined with a key starting with species- (e.g., species-background, species-beam, species-electron1).

Field	Type	Description
noise_seed	int	Random seed for noise initialization
noise_type	string	"gaussian" or "uniform"
noise_val	float	Amplitude of noise
v0	float	Drift velocity in code units of \(\sqrt{T_0/m_e}\) (thermal-σ units). Numeric only — dimensional strings are not supported here
T0	float	Temperature in units of normalizing_temperature. The initialized distribution has velocity variance T0/mass. Numeric only
m	float	Exponent for super-Gaussian distribution \(f \propto \exp[-|v/(\alpha v_{th})|^m]\). 2.0 is Maxwellian. Note: \(\alpha\) is chosen to fix the moment ratio \(\langle v^4\rangle/\langle v^2\rangle = 3\,T_0/\mathrm{mass}\) for all \(m\); the variance equals T0/mass only at m: 2. For flat-top distributions (m > 2) the second-moment temperature diagnostic will read higher than T0 (e.g. ×1.24 at m: 3, ×1.37 at m: 4)
basis	string	Spatial profile type (see below)

Basis Types

uniform: Constant density profile

species-background:
  basis: uniform
  # No additional parameters required

sine: Sinusoidal density perturbation

species-background:
  basis: sine
  baseline: 1.0        # Base density
  amplitude: 1.0e-4    # Perturbation amplitude
  wavenumber: 0.3      # Wavenumber of perturbation

tanh: Hyperbolic tangent profile (for density gradients)

species-background:
  basis: tanh
  baseline: 0.001      # Minimum density
  bump_or_trough: bump # "bump" or "trough"
  center: 2000.0       # Profile center location
  rise: 25.0           # Steepness of transition
  bump_height: 0.999   # Height of bump/trough
  width: 3900.0        # Width of profile

linear: Linear density gradient

species-background:
  basis: linear
  center: 1000.0
  width: 500.0
  rise: 10.0
  gradient scale length: "100um"  # Gradient scale length with units
  val at center: 1.0              # Density value at center

exponential: Exponential density gradient

species-background:
  basis: exponential
  center: 1000.0
  width: 500.0
  rise: 10.0
  gradient scale length: "100um"
  val at center: 1.0

Multi-Species Example

density:
  quasineutrality: true
  species-electron1:
    noise_seed: 420
    noise_type: gaussian
    noise_val: 0.0
    v0: -1.5
    T0: 0.2
    m: 2.0
    basis: sine
    baseline: 0.5
    amplitude: 1.0e-4
    wavenumber: 0.3
  species-electron2:
    noise_seed: 420
    noise_type: gaussian
    noise_val: 0.0
    v0: 1.5
    T0: 0.2
    m: 2.0
    basis: sine
    baseline: 0.5
    amplitude: -1.0e-4
    wavenumber: 0.3

grid

Simulation grid parameters.

Field	Type	Description
dt	float	Timestep (normalized)
nv	int	Number of velocity grid points (not needed for multispecies)
nx	int	Number of spatial grid points
tmin	float	Start time
tmax	float	End time
vmax	float	Upper bound of the velocity grid (not needed for multispecies)
vmin	float	Lower bound of the velocity grid. Optional; defaults to -vmax (symmetric grid). Use to specify an asymmetric velocity extent (not needed for multispecies).
xmin	float	Domain minimum x
xmax	float	Domain maximum x
parallel	false or list of "x", "v"	Axes to parallelize across devices using jax.shard_map. "x" shards the edfdv push and collision operator along the spatial axis; "v" shards the vdfdx push along the velocity axis. Requires the corresponding dimension (nx or nv) to be divisible by the number of JAX devices. Defaults to false. See Running on Multiple GPUs for the memory model, the launch command, and the reverse-mode AD caveat on "v".

The velocity grid is uniform and cell-centered: dv = (vmax - vmin) / nv with cell centers spanning vmin + dv/2 to vmax - dv/2. An asymmetric grid (vmin != -vmax) is useful, for example, for a bump-on-tail distribution where less resolution is needed on one side. As with vmax, choose bounds wide enough that f ~ 0 at both edges so that the spectral velocity push and the zero-flux collision boundary conditions remain accurate.

Note: For multispecies simulations, nv, vmax, and vmin are defined per-species in terms.species and the global values are not used.

Example:

grid:
  dt: 0.1
  nv: 256
  nx: 32
  tmin: 0.0
  tmax: 100.0
  vmax: 6.4
  xmin: 0.0
  xmax: 20.94
  parallel: ["x", "v"]   # shard edfdv+collisions over x, vdfdx over v

save

Configures what data to save and at what times.

Structure

save:
  fields:
    t:
      tmin: 0.0
      tmax: 100.0
      nt: 601
  electron:
    main:
      t:
        tmin: 0.0
        tmax: 100.0
        nt: 11

Save Key	Description
fields	Electric and magnetic field data
<species>	Dict of named distribution function saves for that species
diag-vlasov-dfdt	Time derivative from Vlasov operator (optional)
diag-fp-dfdt	Time derivative from Fokker-Planck operator (optional)

Multiple distribution saves per species

Each species section is a dict of named saves (main, full, monitor, etc.). Multiple saves with different resolutions and time cadences can be configured under the same species key. Each produces a dist-<species>/<label>.nc file in binary/ and is accessible as result.ys["<species>/<label>"].

save:
  fields:
    t:
      tmin: 0.0
      tmax: 100.0
      nt: 601
  electron:
    # Full-resolution dump, saved infrequently
    full:
      t:
        tmin: 0.0
        tmax: 100.0
        nt: 5
      x:
        xmin: 0.0
        xmax: 20.94
        nx: 32
      v:
        vmin: -6.4
        vmax: 6.4
        nv: 512
    # Coarser dump at higher frequency for monitoring
    monitor:
      t:
        tmin: 0.0
        tmax: 100.0
        nt: 101
      x:
        xmin: 0.0
        xmax: 20.94
        nx: 16
      v:
        vmin: -3.0
        vmax: 3.0
        nv: 128

Each named save contains a t sub-key with:

Field	Type	Description
tmin	float	Start time for saving
tmax	float	End time for saving
nt	int	Number of save points

mlflow

Experiment tracking configuration.

Field	Type	Description
experiment	string	MLflow experiment name
run	string	MLflow run name

Example:

mlflow:
  experiment: basic-epw-for-plots
  run: nonlinear

drivers

External electromagnetic drivers. Both ex and ey are dictionaries of named drivers.

Field	Type	Description
ex	dict	Longitudinal electric field drivers
ey	dict	Transverse electric field drivers (for electromagnetic simulations)
ex_stochastic	object	(Optional) Band-limited, time-correlated stochastic forcing of the longitudinal field (see below)

Each driver is identified by a string key (e.g., "0", "1") and has these parameters:

Field	Type	Description
a0	float	Normalized vector potential amplitude of the wave (see below)
k0	float	Wavenumber
w0	float	Frequency
dw0	float	Frequency offset/chirp
t_center	float	Temporal envelope center
t_rise	float	Temporal envelope rise time
t_width	float	Temporal envelope width
x_center	float	Spatial envelope center
x_rise	float	Spatial envelope rise distance
x_width	float	Spatial envelope width
source_type	string	"extended" (default) or "point". Only affects ey drivers. See below.

Driver normalization and a0

The parameter a0 has the same meaning for both ex and ey drivers: the normalized vector potential amplitude of the wave. The difference is how a0 determines the forcing magnitude in each context:

• ex drivers (longitudinal electric field): The driver produces an electric field with amplitude ω · a0, derived from E = -∂A/∂t. This field enters the Vlasov equation directly.

• ey drivers (transverse wave equation source): The driver produces a source term for the wave equation. The source type controls the spatial profile:

  • "extended" (default): Source is S = -ω² · a0 · envelope(x,t) · sin(kx - ωt), spread over the spatial envelope width. The resulting wave amplitude depends on the source geometry.

  • "point": Source is concentrated at a single grid cell (nearest to x_center). The amplitude is calibrated using the vacuum Green’s function so that the outgoing wave has amplitude a0. Point sources radiate equally in both directions; place the source near a boundary with absorbing BCs to get a unidirectional wave.

Example: Single ex driver

drivers:
  ex:
    '0':
      a0: 1.e-2
      k0: 0.3
      w0: 1.1598
      dw0: 0.0
      t_center: 40.0
      t_rise: 5.0
      t_width: 30.0
      x_center: 0.0
      x_rise: 10.0
      x_width: 4000000.0
  ey: {}

Stochastic forcing: ex_stochastic

Drives the plasma with an external longitudinal electric field built from a set of periodic-box Fourier modes whose complex amplitudes evolve as independent Ornstein-Uhlenbeck (OU) processes. This is the standard “compressive, time-correlated” forcing used for driven-turbulence simulations (e.g. ion-acoustic turbulence with Boltzmann electrons):

dE(x, t) = sum_m Re[ a_m(t) exp(i k_m x) ],   k_m = 2 pi m / L

Each a_m(t) has correlation time tau and stationary RMS amplitude. The realization is drawn once from seed on a uniform time grid of spacing dt_update and linearly interpolated in time, so a run is exactly reproducible and independent of the solver timestep.

Field	Type	Default	Description
modes	list[int]	[1]	Box mode numbers m to force (k_m = 2πm/L)
amplitude	float	—	Stationary RMS of each mode’s complex amplitude (code-unit E field)
tau	float	—	OU correlation time (code units)
seed	int	42	RNG seed for the forcing realization
dt_update	float	tau/10	Time resolution of the precomputed OU series (clamped to at most tau/2)

Example (critically balanced box-scale driving — the correlation time is the ion thermal box-crossing time tau = L/v_ti and the amplitude is v_ti^2/L, so the velocity kick per correlation time is ~v_ti; here v_ti = 1, L = 458):

drivers:
  ex: {}
  ey: {}
  ex_stochastic:
    modes: [1]
    amplitude: 2.1834e-3   # v_ti^2 / L
    tau: 458.0             # L / v_ti
    seed: 42

Example: Multiple ey drivers (SRS simulation)

drivers:
  ex: {}
  ey:
    '0':
      a0: 1.0
      k0: 1.0
      w0: 2.79
      dw0: 0.0
      t_center: 4000.0
      t_rise: 20.0
      t_width: 7900.0
      x_center: 50.0
      x_rise: 0.1
      x_width: 1.0
    '1':
      a0: 1.e-6
      k0: 1.0
      w0: 1.63
      dw0: 0.0
      t_center: 4000.0
      t_rise: 20.0
      t_width: 7900.0
      x_center: 3950.0
      x_rise: 0.1
      x_width: 1.0

diagnostics

Enable/disable diagnostic outputs.

Field	Type	Description
diag-vlasov-dfdt	bool	Save df/dt from Vlasov operator
diag-fp-dfdt	bool	Save df/dt from Fokker-Planck operator

Example:

diagnostics:
  diag-vlasov-dfdt: False
  diag-fp-dfdt: False

terms

Solver algorithm configuration.

Field	Type	Description
field	string	Electric field solver: "poisson", "poisson-boltzmann", "ampere", or "hampere"
edfdv	string	Velocity advection scheme: "exponential" or "cubic-spline"
time	string	Time integrator: "sixth" (6th order Hamiltonian) or "leapfrog"
fokker_planck	object	Fokker-Planck collision operator configuration
krook	object	Krook collision operator configuration
hou_li_filter	object	Hou-Li spectral filter (optional, default off)
species	list	(Optional) List of species configurations for multispecies simulations
boltzmann_electrons	object	(Optional) Linearized Boltzmann electron closure parameters, used with field: poisson-boltzmann

species (Multispecies Configuration)

For multispecies simulations, you can define multiple particle species (e.g., electrons and ions) with independent velocity grids and charge-to-mass ratios. Each species is defined as an entry in the terms.species list.

Field	Type	Description
name	string	Species identifier (e.g., "electron", "ion")
charge	float	Charge in units of fundamental charge (e.g., -1.0 for electrons, 10.0 for Z=10 ions)
mass	float	Mass in units of electron mass (e.g., 1.0 for electrons, 1836.0 for protons)
vmax	float	Upper bound of the velocity grid for this species
vmin	float	Lower bound of the velocity grid for this species. Optional; defaults to -vmax (symmetric grid)
nv	int	Number of velocity grid points for this species
density_components	list	List of density component names (from density section) that belong to this species

Important notes:

• Each species can have its own velocity grid parameters (vmax, vmin, nv), allowing heavier species to use smaller velocity extents or asymmetric velocity bounds

• The density_components field maps species to their density profiles defined in the density section

• For single-species simulations, terms.species can be omitted; a default electron species will be auto-generated from the grid.vmax, grid.nv, and all species-* density components

• When using multispecies, grid.vmax and grid.nv are not used (each species defines its own)

Multispecies Example: Ion Acoustic Wave

density:
  quasineutrality: true
  species-electron-background:
    noise_seed: 420
    noise_type: gaussian
    noise_val: 0.0
    v0: 0.0
    T0: 1.0              # Electron temperature
    m: 2.0
    basis: sine
    baseline: 1.0
    amplitude: 1.0e-3
    wavenumber: 0.1
  species-ion-background:
    noise_seed: 421
    noise_type: gaussian
    noise_val: 0.0
    v0: 0.0
    T0: 0.01             # Ion temperature (cooler than electrons)
    m: 2.0
    basis: sine
    baseline: 1.0
    amplitude: 1.0e-3
    wavenumber: 0.1

terms:
  field: poisson
  edfdv: exponential
  time: sixth
  species:
    - name: electron
      charge: -1.0
      mass: 1.0
      vmax: 6.4           # Larger velocity range for light electrons
      nv: 512
      density_components:
        - species-electron-background
    - name: ion
      charge: 10.0        # Z=10 ions
      mass: 18360.0       # Ion mass relative to electron mass
      vmax: 0.15          # Smaller velocity range for heavy ions
      nv: 256
      density_components:
        - species-ion-background
  fokker_planck:
    is_on: False
    # ...

See tests/test_vlasov1d/configs/multispecies_ion_acoustic.yaml for a complete working example.

boltzmann_electrons (kinetic ions with adiabatic electrons)

With field: poisson-boltzmann, only ions are evolved kinetically; the electrons respond adiabatically through the linearized Boltzmann (screened Poisson) closure

e phi / T_e = delta n_i / n_0 + lambda_De^2 d^2/dx^2 (e phi / T_e)

solved spectrally as e phi_k / T_e = (delta n_k / n_0) / (1 + k^2 lambda_De^2). The k=0 force vanishes identically, so no static neutralizing background is needed. This is the standard setup for ion-acoustic wave and IAW-turbulence simulations: it removes all electron timescales from the problem, so the timestep is limited only by ion dynamics.

Field	Type	Default	Description
Te	float	1.0	Electron temperature in code units (units of normalizing_temperature)
lambda_De	float	sqrt(Te / n_0)	Screening length in code units. Omit for the physically consistent value; set explicitly to move the screening scale (e.g. to the grid spacing, mimicking PIC references where the cell size takes the place of the Debye length); 0.0 gives the quasineutral closure e phi / T_e = delta n / n_0

Requirements: terms.species must be given explicitly and all kinetic species must have positive charge (electrons are the closure, not a species).

Ion-unit convention: with a single ion species at charge: 1.0, mass: 1.0 and T0: 1.0, code time is 1/omega_pi, velocity is the ion thermal speed v_ti, and length is v_ti/omega_pi. Then Te is the electron-to-ion temperature ratio Te/Ti, the sound speed is c_s = sqrt(Te), and the electron Debye length is lambda_De = sqrt(Te) in code units. Set units.reference: ion (see units) so that dimensional string inputs and the physical quantities reported by write_units() follow this same ion convention — with the default electron reference they are electron-referenced and therefore nominal for such runs.

Example (ion-acoustic turbulence, Te/Ti = 0.05):

terms:
  field: poisson-boltzmann
  boltzmann_electrons:
    Te: 0.05
  edfdv: cubic-spline
  time: sixth
  species:
    - name: ion
      charge: 1.0
      mass: 1.0
      vmax: 6.4
      nv: 512
      density_components:
        - species-ion-background

See configs/vlasov-1d/iaw-turbulence.yaml for a complete working example with box-scale stochastic forcing, and tests/test_vlasov1d/configs/boltzmann_iaw.yaml for an ion-acoustic dispersion test setup.

solver: vlasov-1d-iaw (IAW turbulence module)

solver: vlasov-1d-iaw selects IAWTurbulence1D, a problem-specific module for driven ion-acoustic turbulence that extends the base Vlasov-1D module with

• an nk save stream: the low-|k| complex spectrum of the total kinetic charge density, sampled densely in time (each sample is only 2 * (nk_modes + 1) floats, so it can be saved far more often than full moment profiles). Written to binary/nk.nc with real/imag parts and P = |n_k|^2.

• IAW post-processing: plots/iaw/density_spectrum.png (late-window-averaged P(m) with a k*lambda_De axis and m^-2 / m^-3/2 guides), plots/iaw/nk_spectrogram.png (P(m, t)), and phase_space_dfx.png (f - _x panels) for every configured distribution save.

Options live in an optional top-level iaw_diagnostics block:

Field	Type	Default	Description
nk_modes	int	1024	Number of box modes in the nk stream (clipped to nx // 2)
nk_nt	int	2001	Number of nk time samples
spectrum_window	[float, float]	[0.5, 1.0]	Averaging window for the spectrum plot, as fractions of tmax

solver: vlasov-1d-iaw

iaw_diagnostics:
  nk_modes: 2048
  nk_nt: 4001
  spectrum_window: [0.4, 1.0]

hou_li_filter

Hou-Li exponential spectral filter applied after each timestep. Damps high-wavenumber modes to suppress numerical oscillations without significantly affecting well-resolved physics. Can be applied in position space (x), velocity space (v), or both.

The filter kernel in Fourier space is:

sigma(j) = exp(-alpha * (j / N)^(2*order))

where j is the mode index, N is the Nyquist mode, and modes near j = N are damped strongly while low modes are left nearly unchanged.

Field	Type	Default	Description
is_on	bool	—	Enable/disable the filter
alpha	float	36.0	Filter strength. The default 36.0 ≈ -log(float64 machine epsilon) ensures the Nyquist mode is zeroed to machine precision
order	int	36	Filter order; higher values give a sharper roll-off that preserves more low-wavenumber content
dimensions	list	["x", "v"]	Dimensions to filter. Can be [], ["x"], ["v"], or ["x", "v"]

Example:

terms:
  hou_li_filter:
    is_on: True
    alpha: 36.0
    order: 36
    dimensions: ["x", "v"]

If hou_li_filter is omitted entirely, it defaults to is_on: False.

Field Solvers

• "poisson": Spectral Poisson solver (works with any time integrator)

• "poisson-boltzmann": Spectral screened-Poisson solver for kinetic ions with linearized Boltzmann electrons (works with any time integrator; see boltzmann_electrons above)

• "ampere": Ampere solver (requires "leapfrog" time integrator)

• "hampere": Hamiltonian Ampere solver (requires "leapfrog" time integrator)

fokker_planck

Fokker-Planck collision operator.

Field	Type	Description
is_on	bool	Enable/disable
type	string	Collision type (see table below)
m	float	Super-Gaussian exponent of the equilibrium (only used by type: super_gaussian; default 2.0 = Maxwellian)
self_consistent_beta	object	Optional Newton refinement of the equilibrium shape parameter (see below)
time	object	Temporal profile
space	object	Spatial profile

Available type values (case-insensitive):

type	Model	Scheme	Equilibrium
lenard_bernstein	Lenard-Bernstein	central differencing	Maxwellian at v=0
chang_cooper	Lenard-Bernstein	Chang-Cooper	Maxwellian at v=0
dougherty	Dougherty	central differencing	Maxwellian at \(\bar v\)
chang_cooper_dougherty	Dougherty	Chang-Cooper	Maxwellian at \(\bar v\)
dougherty_nodrag	Dougherty (diffusion of the deviation)	central differencing	Maxwellian at \(\bar v\)
super_gaussian	Super-Gaussian Dougherty	Chang-Cooper	Super-Gaussian of order m at \(\bar v\)

The Chang-Cooper scheme is positivity-preserving and conserves density exactly; central differencing is provided for comparison only.

super_gaussian

Drift-diffusion operator whose equilibrium is a super-Gaussian \(f_0 \propto \exp(-\beta\,|v-\bar v|^m)\) instead of a Maxwellian (\(m=2\) recovers chang_cooper_dougherty). Use this to maintain a prescribed non-Maxwellian order — e.g. a Langdon/DLM inverse-bremsstrahlung-heated distribution — against collisional relaxation during Landau-damping or LPI runs. Note it does not generate super-Gaussians dynamically (there is no competition between heating and e-e collisions as in the VFP-1D inverse-bremsstrahlung operator); it relaxes toward the prescribed order m.

The drift coefficient is the exact finite difference of the equilibrium potential \(\varphi = \beta|v-\bar v|^m\), so the sampled super-Gaussian is the exact fixed point of the Chang-Cooper discretization. The shape parameter \(\beta\) is set by the energy-conservation closure \(\beta = n/(m\,\langle|v-\bar v|^m\rangle)\) (the generalization of the Lenard-Bernstein/Dougherty \(\beta = 1/(2T)\)), refined by a Newton solve of the discrete energy-flux condition when self_consistent_beta is enabled.

Conservation properties:

• Density: exact (zero-flux boundary conditions).

• Energy: no secular drift at equilibrium when self_consistent_beta is enabled (recommended for long collisional runs; without it the continuum closure leaves an O(dv²) quadrature residual that drifts T by ~5e-7 per collision time at nv=128 for m=3 — already negligible unless many collision times elapse). A single Newton step (max_steps: 1) reduces the drift to machine level. Off-equilibrium transients carry a one-time O(nu·dt) offset from operator splitting, negligible at production nu·dt.

• Momentum: exact for distributions symmetric about \(\bar v\); skewed transients exchange momentum at O(skewness), unlike \(m=2\) where \(C \propto (v-\bar v)\) makes it exact.

Note on initialization: the super-Gaussian initializer (m in the species config) uses a width convention that fixes \(\langle v^4\rangle/\langle v^2\rangle\), while this operator preserves the super-Gaussian carrying the distribution’s own variance. Any super-Gaussian of order m is (approximately) a fixed point regardless of width convention, so the two compose without drift.

terms:
  fokker_planck:
    is_on: True
    type: super_gaussian
    m: 3.0
    self_consistent_beta:
      enabled: True
      max_steps: 1   # one Newton step suffices (the closure lands within O(dv²) of the root)
    time:
      baseline: 1.0e-5
      # ...
    space:
      baseline: 1.0
      # ...

dougherty_nodrag

Diffusion-of-the-deviation operator \(C[f] = \nu\, T\, \partial_v^2 \left(f - f_M[n, \bar v, T]\right)\). It uses the same moments machinery as dougherty, but the implicit solve runs with the drag coefficient zeroed (pure diffusion) and the diffusion of the self-consistent Maxwellian is subtracted explicitly using the same discrete zero-flux stencil, so \(f_M\) is a discrete fixed point and \(n\), \(P\), and \(E\) are conserved. Because the operator carries no drag acting on the deviation, it is parity-preserving in \(v - v_\phi\) — useful as a control for isolating the role of collisional drag (e.g. in ratchet/pawl-style transport experiments) while retaining Maxwellian-preserving velocity diffusion.

self_consistent_beta

Optional sub-object controlling the Newton refinement of the equilibrium shape parameter β. For the Maxwellian operators it matches the discrete temperature of the equilibrium to that of f (eliminating equilibrium drift in Chang-Cooper schemes); for super_gaussian it solves the discrete energy-flux condition.

Field	Type	Default	Description
enabled	bool	False	Enable the Newton solve
max_steps	int	3	Maximum Newton iterations
rtol	float	1e-8	Relative tolerance
atol	float	1e-12	Absolute tolerance

krook

Krook collision operator (relaxation towards Maxwellian).

Field	Type	Description
is_on	bool	Enable/disable
time	object	Temporal profile
space	object	Spatial profile

Profile Objects (time/space)

Both fokker_planck and krook use profile objects to define spatiotemporal variation of collision frequency:

Field	Type	Description
baseline	float	Base collision frequency in code units of \(\omega_{p0}\) (i.e. \(\hat\nu = \nu[\mathrm{s}^{-1}]/\omega_{p0}\), a true rate — no \(2\pi\)). For reference, the NRL e–e rate in these units is logged as nuee_norm in units.yaml. Note the Dougherty/Lenard–Bernstein \(\nu\) agrees with the NRL \(\nu_{ee}\) only up to an O(1) factor
bump_or_trough	string	"bump" or "trough"
center	float	Profile center
rise	float	Transition steepness
slope	float	Linear slope (typically 0.0)
bump_height	float	Height of bump/depth of trough
width	float	Profile width

Example: With Fokker-Planck collisions

terms:
  field: poisson
  edfdv: cubic-spline
  time: sixth
  fokker_planck:
    is_on: True
    type: Dougherty
    time:
      baseline: 1.0e-5
      bump_or_trough: bump
      center: 0.0
      rise: 25.0
      slope: 0.0
      bump_height: 0.0
      width: 100000.0
    space:
      baseline: 1.0
      bump_or_trough: bump
      center: 0.0
      rise: 25.0
      slope: 0.0
      bump_height: 0.0
      width: 100000.0
  krook:
    is_on: False
    time:
      baseline: 1.0
      bump_or_trough: bump
      center: 0.0
      rise: 25.0
      slope: 0.0
      bump_height: 0.0
      width: 100000.0
    space:
      baseline: 1.0
      bump_or_trough: bump
      center: 0.0
      rise: 25.0
      slope: 0.0
      bump_height: 0.0
      width: 100000.0

Complete Examples

Electron Plasma Wave (EPW)

See configs/vlasov-1d/epw.yaml - driven electron plasma wave with Fokker-Planck collisions.

Bump-on-Tail Instability

See configs/vlasov-1d/bump-on-tail.yaml - two-species configuration with beam distribution.

Two-Stream Instability

See configs/vlasov-1d/twostream.yaml - counter-streaming electron beams.

Stimulated Raman Scattering (SRS)

See configs/vlasov-1d/srs.yaml - electromagnetic simulation with ey drivers.

Nonlinear EPW with Initial Condition

See configs/vlasov-1d/nlepw-ic.yaml - large-amplitude initial perturbation without external driver.

Ion Acoustic Wave (Multispecies)

See tests/test_vlasov1d/configs/multispecies_ion_acoustic.yaml - two-species (electron + ion) simulation demonstrating multispecies support.

Ion-Acoustic Turbulence (Boltzmann Electrons)

See configs/vlasov-1d/iaw-turbulence.yaml - kinetic ions with the linearized Boltzmann electron closure (field: poisson-boltzmann), driven by box-scale time-correlated stochastic forcing (drivers.ex_stochastic).