caustic

A 6D Vlasov–Poisson solver framework for collisionless gravitational dynamics.

caustic is a modular, general-purpose library for solving the Vlasov–Poisson equations in full 6D phase space (3 spatial + 3 velocity dimensions). It targets astrophysical problems — dark matter halo formation, galaxy dynamics, tidal streams, stellar system stability — that are traditionally handled by N-body methods but suffer from artificial collisionality and loss of fine-grained phase-space structure.

The library provides a pluggable architecture where the phase-space representation, Poisson solver, time integrator, and initial condition generator can be swapped independently.

Why not N-body?

N-body simulations sample the distribution function with discrete particles. This introduces noise and artificial two-body relaxation that destroys exactly the structures a collisionless solver should preserve: caustic surfaces, thin phase-space streams, and the true velocity distribution at any point. caustic solves the governing equation directly — no particles, no sampling noise, no artificial collisionality.

Quick start

Add to your Cargo.toml:

[dependencies]
caustic = "0.0.4"

Minimal example: Plummer sphere equilibrium

use caustic::prelude::*;
use caustic::{
    FftPoisson, PlummerIC, SemiLagrangian, StrangSplitting,
    SpatialBoundType, VelocityBoundType, sample_on_grid,
};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let domain = Domain::builder()
        .spatial_extent(20.0)      // [-20, 20]^3 in natural units
        .velocity_extent(3.0)      // [-3, 3]^3
        .spatial_resolution(32)    // 32^3 spatial grid
        .velocity_resolution(32)   // 32^3 velocity grid
        .t_final(50.0)
        .spatial_bc(SpatialBoundType::Periodic)
        .velocity_bc(VelocityBoundType::Open)
        .build()?;

    // Set up a Plummer sphere: mass=1, scale_radius=1, G=1
    let ic = PlummerIC::new(1.0, 1.0, 1.0);
    let snap = sample_on_grid(&ic, &domain);

    let poisson = FftPoisson::new(&domain);
    let mut sim = Simulation::builder()
        .domain(domain)
        .poisson_solver(poisson)
        .advector(SemiLagrangian::new())
        .integrator(StrangSplitting::new(1.0))
        .initial_conditions(snap)
        .time_final(50.0)
        .build()?;

    let exit = sim.run()?;
    exit.print_summary();
    Ok(())
}

Architecture

Each solver component is a Rust trait; implementations are swapped independently:

The PhaseSpaceRepr trait

The central abstraction. All phase-space storage strategies implement this interface:

pub trait PhaseSpaceRepr: Send + Sync {
    /// Integrate f over all velocities: rho(x) = integral f dv^3.
    fn compute_density(&self) -> DensityField;

    /// Drift sub-step: advect f in spatial coordinates by dx = v*dt.
    fn advect_x(&mut self, displacement: &DisplacementField, dt: f64);

    /// Kick sub-step: advect f in velocity coordinates by dv = g*dt.
    fn advect_v(&mut self, acceleration: &AccelerationField, dt: f64);

    /// Compute velocity moment of order n at given spatial position.
    fn moment(&self, position: &[f64; 3], order: usize) -> Tensor;

    /// Total mass M = integral f dx^3 dv^3.
    fn total_mass(&self) -> f64;

    /// Casimir invariant C_2 = integral f^2 dx^3 dv^3.
    fn casimir_c2(&self) -> f64;

    /// Boltzmann entropy S = -integral f ln f dx^3 dv^3.
    fn entropy(&self) -> f64;

    /// Number of distinct velocity streams at each spatial point.
    fn stream_count(&self) -> StreamCountField;

    /// Extract the local velocity distribution f(v|x) at a given position.
    fn velocity_distribution(&self, position: &[f64; 3]) -> Vec<f64>;

    /// Total kinetic energy T = 1/2 integral f v^2 dx^3 dv^3.
    fn total_kinetic_energy(&self) -> f64;

    /// Extract a full 6D snapshot of the current state.
    fn to_snapshot(&self, time: f64) -> PhaseSpaceSnapshot;
}

Initial conditions

All implemented ICs satisfy the IsolatedEquilibrium trait and can be sampled onto a grid with sample_on_grid():

• PlummerIC — Plummer sphere via analytic distribution function f(E)
• KingIC — King model via Poisson-Boltzmann ODE (RK4 integration)
• HernquistIC — Hernquist profile via closed-form f(E)
• NfwIC — NFW profile via numerical Eddington inversion
• ZeldovichSingleMode — single-mode Zel'dovich pancake (cosmological)
• MergerIC — two-body superposition f = f_1 + f_2 with offsets
• TidalIC — progenitor equilibrium model in an external host potential
• CustomIC / CustomICArray — user-provided callable or pre-computed array

Diagnostics

Conserved quantities monitored each timestep via GlobalDiagnostics:

• Total energy (kinetic + potential), momentum, angular momentum
• Casimir C_2, Boltzmann entropy
• Virial ratio, total mass in box
• Density profile (radial binning)

Additional output modules: VelocityMoments (surface density, J-factor), PhaseSpaceDiagnostics (power spectrum, growth rates), CausticDetector (caustic surface detection, first caustic time).

Validation suite

Run with cargo test --release -- --test-threads=1:

Plus 2 integration tests (smoke_test, end_to_end_run) exercising the full pipeline from Domain through Simulation::run() to ExitPackage.

Feature flags

[dependencies]
caustic = { version = "0.0.4", features = ["jemalloc"] }

Performance

• Parallelism: rayon data parallelism across all hot paths (compute_density, advect_x, advect_v, FFT axes)
• Release profile: fat LTO, codegen-units = 1, target-cpu=native (via .cargo/config.toml)
• Benchmarks: criterion benchmarks (cargo bench), benchmark binary: solver_kernels
• Instrumentation: tracing::info_span! on all hot methods (zero overhead without a subscriber)
• Profiling profile: [profile.profiling] inherits release with debug symbols for perf/samply

Roadmap

• Uniform 6D grid with rayon parallelism
• FFT Poisson (periodic + Hockney-Eastwood isolated)
• Semi-Lagrangian advection (Catmull-Rom) + Strang/Lie/Yoshida splitting
• Isolated equilibrium ICs (Plummer, King, Hernquist, NFW)
• Cosmological, merger, tidal, and custom ICs
• Conservation diagnostics + 11-test validation suite
• Criterion benchmarks + tracing instrumentation
• Binary snapshot I/O, CSV diagnostics, JSON checkpoints
• Tensor-train low-rank representation
• Lagrangian sheet tracker for cold dark matter
• Multigrid / spherical harmonics Poisson solvers
• Adaptive mesh refinement
• GPU acceleration
• MPI domain decomposition

Companion: phasma

phasma is a ratatui-based terminal UI that consumes caustic as a library dependency. It provides interactive parameter editing, live diagnostics rendering, density/phase-space heatmaps, energy conservation plots, and radial profile charts — all from the terminal. phasma contains no solver logic; it delegates entirely to caustic.

Minimum supported Rust version

Rust edition 2024, targeting stable Rust 1.75+.

License

This project is licensed under the GNU General Public License v3.0. See LICENSE for details.

Citation

If you use caustic in academic work, please cite:

@software{caustic,
  title  = {caustic: A 6D Vlasov--Poisson solver framework},
  url    = {https://github.com/resonant-jovian/caustic},
  year   = {2026}
}