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.

Architecture

Core modules

Quick start

Add to your Cargo.toml:

[dependencies]
caustic = "0.1"

Minimal example: Plummer sphere equilibrium

use caustic::prelude::*;

fn main() -> Result<(), CausticError> {
    // Define the computational domain
    let domain = Domain::builder()
        .spatial_extent(20.0)      // [-20, 20]^3 in natural units
        .velocity_extent(3.0)      // [-3, 3]^3
        .spatial_resolution(64)    // 64^3 spatial grid
        .velocity_resolution(64)   // 64^3 velocity grid
        .boundary(Boundary::Isolated)
        .build()?;

    // Set up a Plummer sphere equilibrium
    let ic = PlummerModel::new()
        .mass(1.0)
        .scale_radius(1.0)
        .build();

    // Choose solver components
    let repr = UniformGrid6D::new(&domain);
    let poisson = FftIsolated::new(&domain);
    let advector = SemiLagrangian::new();
    let integrator = StrangSplitting::new();

    // Build the simulation
    let mut sim = Simulation::builder()
        .domain(domain)
        .representation(repr)
        .poisson_solver(poisson)
        .advector(advector)
        .integrator(integrator)
        .initial_conditions(ic)
        .time_final(50.0)          // 50 dynamical times
        .output_interval(1.0)
        .exit_on_energy_drift(1e-4)
        .build()?;

    // Run
    sim.run()?;

    Ok(())
}

Custom initial conditions

use caustic::prelude::*;

let ic = CustomIC::from_fn(|x, v| {
    // Two overlapping Gaussians — a simple merger setup
    let r1_sq = (x[0] - 3.0).powi(2) + x[1].powi(2) + x[2].powi(2);
    let r2_sq = (x[0] + 3.0).powi(2) + x[1].powi(2) + x[2].powi(2);
    let v_sq = v[0].powi(2) + v[1].powi(2) + v[2].powi(2);
    let sigma = 0.5;

    let f1 = (-r1_sq / 2.0).exp() * (-v_sq / (2.0 * sigma * sigma)).exp();
    let f2 = (-r2_sq / 2.0).exp() * (-v_sq / (2.0 * sigma * sigma)).exp();
    f1 + f2
});

Representations

The central algorithmic challenge is representing a 6D function efficiently. caustic treats this as a trait so that novel representations can be developed and benchmarked within the same framework:

pub trait PhaseSpaceRepr: Send + Sync {
    /// Integrate f over all velocities to produce the 3D density field.
    fn compute_density(&self) -> DensityField;

    /// Advect f in spatial coordinates by displacement field.
    fn advect_x(&mut self, displacement: &DisplacementField, dt: f64);

    /// Advect f in velocity coordinates by acceleration field.
    fn advect_v(&mut self, acceleration: &AccelerationField, dt: f64);

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

    /// Total mass (integral of f over all phase space).
    fn total_mass(&self) -> f64;

    /// Casimir invariant C_2 = ∫ f² dx³ dv³.
    fn casimir_c2(&self) -> f64;
}

Shipped representations

• UniformGrid6D — brute-force uniform grid. Simple, correct, useful for validation. Memory: O(N⁶).
• TensorTrain (work in progress) — low-rank tensor decomposition exploiting compressibility. Memory: O(N³ r³) where r is the rank.
• SheetTracker (work in progress) — Lagrangian sheet representation for cold (dark matter) initial conditions. Memory: O(N³).

Validation

caustic includes a built-in test suite against known analytic solutions:

cargo test --release -- --test-threads=1

Feature flags

[dependencies]
caustic = { version = "0.1", features = ["hdf5", "mpi", "cosmological"] }

Roadmap

• Uniform 6D grid representation
• FFT Poisson solver (periodic + isolated)
• Semi-Lagrangian advection with Strang splitting
• Plummer, King, Hernquist initial conditions
• Energy/momentum/Casimir conservation diagnostics
• Tensor-train representation
• Sheet-tracking representation for cold ICs
• Hybrid sheet/grid representation
• Spectral velocity-space (Hermite) representation
• Multigrid Poisson solver
• Adaptive time-stepping
• MPI domain decomposition
• Cosmological (comoving) mode
• GPU acceleration (wgpu compute shaders)

Minimum supported Rust version

caustic targets 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}
}