The Role of Axion Stars in the Early Universe
A Schrödinger–Poisson–hydrodynamics simulator for ultralight dark matter and early star formation.
Developed as part of a Jugend forscht research project by Ben Waldmann (17), Ulf-Merbold-Gymnasium Greiz
HaloGenesis is a numerical simulation framework designed to explore whether axion-stars stable, quantum-mechanical density configurations of ultralight axions, could have played a role in triggering the formation of the first stars in the Universe.
Motivation The standard cosmological model (ΛCDM) faces well-known challenges:
It predicts an overabundance of low-mass galaxies It struggles to accurately reproduce the inner structure of dark matter halos
Fuzzy Dark Matter (FDM), composed of ultralight axions (mₐ ≈ 10⁻²² eV), provides a potential resolution. Due to their wave-like nature, these particles generate a quantum pressure that suppresses small-scale structure formation.
Hypothesis
Axion stars may act as pre-galactic gravitational attractors for baryonic gas.
They can locally compress gas to such a degree that Jeans instability occurs earlier than in classical cold dark matter (CDM) minihalos—potentially even before hierarchical structure formation is complete.
The key methodological contribution of this project is a modified Schrödinger–Poisson-framework with asymmetric coupling between dark matter and baryonic gas.
Timescale Separation Gas response timescale: τ_gas ≈ 187 Myr (sound-crossing time across the solitonic core) Axion field oscillation timescale: τ_Φ ≈ 0.1 years Ratio: τ_gas / τ_Φ ≈ 10⁹
This extreme separation implies:
Baryonic gas cannot dynamically respond to rapid axionic oscillations.
Φ̄(x, t) = (1/ΔT) ∫ Φ(x, t') dt'
This leads to the following approximation of asymmetric coupling:
- Axion field → responds to the full gravitational potential Φ
- Baryonic gas → responds only to the averaged potential Φ̄
This is not a numerical shortcut, but a physically justified limit arising from timescale separation.
Axion field (Schrödinger equation):
iħ ∂ₜψ = [−ħ²/(2mₐ)∇² + mₐΦ] ψ
Self-gravity (Poisson equation):
∇²Φ = 4πG (ρₐ + ρ_b) with ρₐ = mₐ|ψ|²
Baryonic gas (compressible Euler equations):
∂ₜρ_b + ∇·(ρ_b v) = 0
∂ₜ(ρ_b v) + ∇·(ρ_b v⊗v + P I) = −ρ_b ∇Φ̄
∂ₜE + ∇·[(E+P)v] = −ρ_b v·∇Φ̄
P = (γ−1)(E − ½ρ_b|v|²), γ = 5/3
The baryons are driven by the adiabatically averaged potential Φ̄ - a time-averaged version of Φ that filters out the rapid quantum oscillations of the axion field at frequency ω ~ mₐc²/ħ.
The stable ground-state solution of the Schrödinger–Poisson system is the soliton (Schive et al. 2014, PRL 113, 261302):
ρ_sol(r) = ρ_c · [1 + 0.091 (r/r_c)²]⁻⁸
where ρ_c is the central density and r_c is the core radius (half-maximum point). The soliton is stabilised against gravitational collapse by the quantum pressure:
Q = −ħ²/(2mₐ) ∇²√ρₐ / √ρₐ
This is the fundamental observable that the simulation validates against.
Throughout the code, G = ħ = mₐ = 1. Physical conversions are
provided in utils/units.py.
| Quantity | Code unit | Physical value (m_a = 10⁻²² eV) |
|---|---|---|
| Length | L | ~ 1 kpc (for L=1 in box) |
| Time | t₀ | ~ 100 Myr |
| Mass | M₀ | ~ 10⁸ M_sun |
Yoshida-4 symplectic splitting (default, schroedinger_order=4):
S₄(Δt) = S₂(w₁Δt) ∘ S₂(w₀Δt) ∘ S₂(w₁Δt)
w₁ = 1/(2 − 2^(1/3)) ≈ 1.351
w₀ = 1 − 2w₁ ≈ −1.702 (backward sub-step, unitary)
Each Strang sub-step S₂(τ) = V(τ/2) K(τ) V(τ/2):
- V(τ) — pointwise phase rotation exp(−imₐΦτ/ħ) in real space
- K(τ) — kinetic propagator exp(−iħk²τ/2mₐ) in Fourier space
Properties:
- 4th-order accuracy in Δt (vs 2nd-order for plain Strang)
- Norm ‖ψ‖² preserved to machine precision (unitary operators)
- Allows Δt 3–5× larger than Strang for equal accuracy
- 6 FFTs per global step (vs 2 for Strang)
Also available: schroedinger_order=2 (Strang, 2 FFTs/step).
Spectral inversion in O(N log N):
Φ̃(k) = −4πG ρ̃(k) / k², Φ̃(k=0) = 0 (zero mean field)
Spectral accuracy, periodic boundary conditions.
MUSCL-HLLC finite volume with HLL fallback:
- MUSCL reconstruction — linear interpolation to cell faces with minmod slope limiter (2nd-order spatial accuracy)
- HLLC Riemann solver (Toro et al. 1994) — exact for contact and shear waves; 4-region upwind selection based on wave speeds S_L, S_R, S*
- HLL fallback — activated cell-by-cell when |S_R − S_L| < ε·max(|S_L|,|S_R|) (vacuum states, near-zero pressure); unconditionally entropy-satisfying
- SSP-RK2 time integration (Shu-Osher 1988)
- Spectral gravity source — ∂ᵢΦ̄ computed via FFT for consistency with the Poisson solver
Strang splitting for the coupled system:
S(Δt) = S_schr(Δt/2) → S_grav → S_hydro(Δt) → S_grav → S_schr(Δt/2)
The baryons see Φ̄, a ring-buffer mean of the last n_phi_avg potential
snapshots. This adiabatic averaging decouples the fast axion oscillation
from the baryon dynamics — physically correct because the baryon sound
crossing time >> axion oscillation period.
Two CFL conditions are enforced simultaneously:
dt_schr = cfl_schr · dx² · mₐ / ħ (Schrödinger phase velocity)
dt_hydro = cfl_hydro · dx / max(cs + |v|) (acoustic signal speed)
dt = clip(min(dt_schr, dt_hydro), dt_min, dt_max)
The active limiter (Schrödinger or hydro) is printed at each progress line. During collapse, dt_hydro tightens as the gas accelerates.
halogenesis/
├ core/
│ ├ backend.py NumPy / CuPy abstraction (GPU support)
│ ├ grid.py Cartesian grid, k-vectors, cell volumes
│ ├ schroedinger.py Strang-2 and Yoshida-4 split-step solver
│ ├ poisson.py FFT Poisson solver, Φ̄ averaging
│ ├ hydro.py MUSCL-HLLC + HLL-fallback Euler solver
│ └ coupling.py Strang operator-splitting time integrator
├ output/
│ ├ plots.py Density slices, power spectra, diagnostics
│ ├ animation.py MP4 export via FuncAnimation
│ ├ checkpoints.py NPZ / HDF5 state I/O
│ └ vtk_export.py VTK/ParaView 3D export
├ utils/
│ ├ units.py Physical constants and unit converters
│ ├ ics.py Soliton, Gaussian-random, uniform-baryon ICs
│ ├ diagnostics.py Mass, energy, norm conservation monitors
│ └ soliton_validation.py Schive-2014 fit, radial profile, validation plot
├ tests/
│ ├ test_schroedinger.py 10 tests: norm, convergence order, phase
│ ├ test_hydro.py 12 tests: conservation, Sod, HLL, CFL
│ ├ test_poisson.py 6 tests: plane-wave, linearity, symmetry
│ └ test_soliton_validation.py 8 tests: Schive fit, profile recovery
├ config.py SimConfig dataclass (all parameters + defaults)
├ simulation.py Simulation runner with adaptive CFL loop
├ run.py Entry point with argparse (soliton scenario)
└ requirements.txt Dependencies
git clone https://github.com/benwaldmann/HaloGenesis.git
cd HaloGenesis
pip install -e .Dependencies:
pip install numpy scipy matplotlib h5py pyyaml tqdm
# Optional, aber empfohlen für Performance:
pip install pyfftw numba
# Optional für GPU:
pip install cupy-cuda12x# Install dependencies
pip install numpy scipy matplotlib h5py
# Run from inside the HaloGenesis folder
cd HaloGenesis
# Quick smoke test (~5 m on CPU)
python run.py --N 32 --t-end 0.05
# Standard run: N=64 3D soliton (~10-15 min on CPU)
python run.py
# 2D run (much faster, good for development)
python run.py --N 128 --dim 2 --t-end 1.0
# Production 3D with GPU
pip install cupy-cuda12x
python run.py --N 128 --backend cupy
# Experimentation example
python /run.py --ic two_solitons --N 128 --dim 3 --L 1.0 --t-end 1.0 --m-axion 0.595 --rcore 0.1 --ic-seed 11121 --merger-separation 0.65 --merger-phase 0.3 --n-phi-avg 32 --dt-min 5e-4 --jeans --order 6 --split 4
# Run tests (stdlib only, no pytest needed)
python -m unittest discover -s tests -p "test_*.py" -vfrom halogenesis import Simulation, SimConfig
cfg = SimConfig(
N=64, L=1.0, dim=3,
t_end=0.5,
ic_type="soliton", # or "gaussian_random"
adaptive_dt=True,
schroedinger_order=4, # Yoshida-4
backend="auto", # CuPy if GPU available
output_dir="my_output",
)
sim = Simulation(cfg)
sim.setup()
sim.run()
# Post-run: inspect adaptive timestep statistics
print(sim.dt_breakdown())import numpy as np
from halogenesis import Simulation, SimConfig
from halogenesis.core.grid import Grid
from halogenesis.utils.ics import soliton_psi, uniform_baryon
cfg = SimConfig(N=64, L=1.0, dim=3, ic_type="custom")
sim = Simulation(cfg)
# Two solitons on a collision course
grid = sim.grid
psi1 = soliton_psi(grid, r_core=0.06, rho_c=1.0, center=[0.3, 0.5, 0.5])
psi2 = soliton_psi(grid, r_core=0.06, rho_c=1.0, center=[0.7, 0.5, 0.5])
# Add relative velocity via Galilean boost
x = np.linspace(0, 1, cfg.N, endpoint=False)
kboost = 20.0 # k in code units
psi1 *= np.exp(+1j * kboost * x[:, None, None])
psi2 *= np.exp(-1j * kboost * x[:, None, None])
psi = psi1 + psi2
U = uniform_baryon(grid, rho_b0=0.05, gamma=5/3)
sim.setup(psi=psi, U=U)
sim.run()All output is written to output_dir/ (default: halogenesis_output/):
| Format | Output |
|---|---|
checkpoint_NNNNNN.npz |
ψ, U, Φ, t, step (restart-capable) |
density_slice_NNNNNN.pdf |
Mid-plane ρₐ, ρ_b, Φ slices |
power_spectrum_NNNNNN.pdf |
P(k) of axion density field |
soliton_validation_NNNNNN.pdf |
Radial profile + Schive-2014 fit |
diagnostics_time.pdf |
M, ‖ψ‖², E_grav vs time |
fields_NNNNNN.vti |
VTK for ParaView/VisIt (if vtk installed) |
import numpy as np
from halogenesis import Simulation, SimConfig
cfg = SimConfig(N=64, t_end=2.0, output_dir="halogenesis_output")
sim = Simulation(cfg)
ckpt = np.load("halogenesis_output/checkpoint_000500.npz")
sim.setup(psi=ckpt["psi"], U=ckpt["U"])
sim.t = float(ckpt["t"])
sim.step = int(ckpt["step"])
sim.run()All parameters live in SimConfig (dataclass with full validation). Some are presented in the following:
| Parameter | Default | Note |
|---|---|---|
N |
128 | Grid points per dimension (must be even) |
L |
1.0 | Box side length [code units] |
dim |
3 | Spatial dimensions (2 or 3) |
| Parameter | Default | Note |
|---|---|---|
dt |
1e-3 | Initial timestep [code units] |
t_end |
1.0 | End time [code units] |
adaptive_dt |
True | Enable CFL-adaptive Δt |
cfl_hydro |
0.4 | Hydro CFL safety factor |
cfl_schr |
0.5 | Schrödinger CFL safety factor |
dt_min |
1e-7 | Hard minimum Δt |
dt_max |
0.1 | Hard maximum Δt |
dt_report_every |
100 | Console progress interval |
| Parameter | Default | Note |
|---|---|---|
hbar |
1.0 | Reduced Planck constant [code] |
m_axion |
1.0 | Axion mass [code] |
G |
1.0 | Gravitational constant [code] |
gamma_eos |
5/3 | Adiabatic index |
rho_b0 |
0.1 | Background baryon density [code] |
n_phi_avg |
10 | Adiabatic Φ̄ averaging window (steps) |
| Parameter | Default | Note |
|---|---|---|
splitting_order |
2 | Operator splitting: 1=Lie, 2=Strang |
schroedinger_order |
4 | Schrödinger: 2=Strang, 4=Yoshida-4 |
hydro_solver |
"hllc" | Riemann solver (always uses HLL fallback) |
| Parameter | Default | Note |
|---|---|---|
ic_type |
"soliton" | "soliton", "gaussian_random", "two_solitons" or "custom" |
ic_seed |
42 | Random seed for gaussian_random IC |
The test suite uses only Python's standard library (unittest) — no pytest
required.
python -m unittest discover -s tests -p "test_*.py" -v36 tests covering:
- Schrödinger (10): norm conservation to 10⁻¹², convergence orders p≥2 and p≥4, free-particle dispersion, Yoshida coefficients, quantum pressure
- Hydro (12): mass/momentum/energy conservation, positivity, Sod shock tube, HLL fallback (vacuum, NaN check), acoustic convergence p≥1.5, rest state
- Poisson (6): plane-wave solution to 10⁻¹⁰, zero mean, linearity, symmetry
- Soliton validation (8): Schive profile, half-max, core mass, radial profile, parameter recovery to <5%, graceful failure on degenerate input
The simulation is validated against three analytical benchmarks:
-
Norm conservation: ‖ψ‖² must be constant to machine precision. Measured by
Diagnostics.conservation_drift(). -
Soliton profile: After t > several dynamical times, the azimuthal radial profile ρₐ(r) should fit the Schive-2014 formula with L₂_rel < 5%. Measured automatically by
soliton_validation.plot_soliton_validation(). -
Acoustic waves: A linear acoustic wave in the baryon gas should propagate with the correct phase velocity cs = √(γP/ρ) and show temporal convergence order p ≥ 1.5 (MUSCL). Tested in
tests/test_hydro.py.
- Schive, H.-Y., Chiueh, T., Broadhurst, T. (2014). Cosmic structure as the quantum interference of a coherent dark wave. Nature Physics 10, 496–499.
- Mocz, P. et al. (2017). Galaxy formation with BECDM. MNRAS 471, 4559–4570.
- Toro, E.F. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed. Springer. Ch. 10 (HLLC solver).
- Yoshida, H. (1990). Construction of higher order symplectic integrators. Physics Letters A 150, 262–268.
- Shu, C.-W., Osher, S. (1988). Efficient implementation of essentially non-oscillatory shock capturing schemes. JCP 77, 439–471.
Ben (Class 10a, Ulf-Merbold-Gymnasium Greiz) — Jugend forscht 2025/26
HaloGenesis is an independent research project in computational astrophysics. The simulation code was written from scratch in Python using only NumPy/SciPy. In case of questions, make a pull request or write me under my e-mail waldmannben93@gmail.com! Have fun!