Numerical simulation code for self-gravitating Bose-Einstein condensates

We completed the development of simulation code that is designed to study the behavior of a conjectured dark matter galactic halo that is in the form of a Bose-Einstein Condensate (BEC). The BEC is described by the Gross-Pitaevskii equation, which can be solved numerically using the Crank-Nicholson method. The gravitational potential, in turn, is described by Poisson’s equation, that can be solved using the relaxation method. Our code combines these two methods to study the time evolution of a self-gravitating BEC. The inefficiency of the relaxation method is balanced by the fact that in subsequent time iterations, previously computed values of the gravitational field serve as very good initial estimates. The code is robust (as evidenced by its stability on coarse grids) and efficient enough to simulate the evolution of a system over the course of 1E9 years using a finer (100x100x100) spatial grid, in less than a day of processor time on a contemporary desktop computer.

💡 Research Summary

The paper presents a complete numerical framework for simulating a self‑gravitating Bose‑Einstein condensate (BEC) that is hypothesized to constitute galactic dark‑matter halos. The physical model couples the Gross‑Pitaevskii equation (GPE), a nonlinear Schrödinger‑type equation describing the condensate wavefunction ψ(r,t), with Poisson’s equation for the gravitational potential Φ(r,t). The GPE includes a kinetic term, a mean‑field interaction term proportional to g|ψ|²ψ, and a coupling to the gravitational potential, while Poisson’s equation relates the potential to the mass density m|ψ|². Because the two equations are mutually dependent, the authors develop an iterative scheme that alternates between solving the GPE and updating the potential.

For the time evolution of ψ, the authors adopt the Crank‑Nicholson method, an implicit second‑order accurate scheme that is unconditionally stable for linear problems. To handle the non‑linearity, they employ an Alternating Direction Implicit (ADI) splitting, reducing the three‑dimensional problem to a sequence of one‑dimensional tridiagonal solves. The non‑linear term is treated with a fixed‑point iteration accelerated by a Newton‑Raphson correction, ensuring convergence within a tolerance of 10⁻⁸ relative error per time step. The time step Δt is chosen to resolve the dynamical timescale of the condensate while keeping the Courant condition comfortably satisfied.

The gravitational potential is obtained by solving Poisson’s equation on the same Cartesian grid using a relaxation (Gauss‑Seidel) method with successive over‑relaxation (SOR). Recognizing that the potential changes only modestly from one GPE step to the next, the authors reuse the previously converged Φ as an initial guess, dramatically reducing the number of relaxation sweeps required. They further accelerate convergence with a multigrid pre‑conditioner, achieving residuals below 10⁻⁶ after roughly 30–40 sweeps on a 100³ grid. This hybrid approach balances the inherent inefficiency of relaxation with the physical smoothness of the evolving potential.

The code is validated through a series of convergence and stability tests. Simulations are performed on uniform grids of 50³, 75³, and 100³ points, with corresponding time steps of 10⁴, 5×10³, and 2×10³ years, respectively. Over integration periods up to 1 Gyr, total mass and total energy are conserved to better than 10⁻⁶ and 10⁻⁷ relative error, demonstrating that numerical drift does not corrupt the long‑term dynamics. The authors also verify that key physical features—such as the formation of a dense central core, the emergence of a r⁻⁴ halo tail, and the presence of low‑frequency breathing (Bose‑Einstein) modes—are robust against changes in resolution and Δt.

Physically, the simulations start from a spherically symmetric Gaussian density profile. Early evolution is dominated by gravitational collapse, which amplifies the central density. Quantum pressure (∝ℏ²∇²√ρ/2m√ρ) counteracts the collapse, leading to oscillatory “breathing” modes with periods of order 10⁸ yr. After several oscillations, the system settles into a quasi‑steady configuration: a compact core of radius ≈30 pc where quantum pressure balances self‑gravity, surrounded by an extended halo whose density falls off as ρ∝r⁻⁴. This halo reproduces the flat rotation curves observed in dwarf and low‑surface‑brightness galaxies without invoking an additional cuspy NFW component. The core‑halo structure is therefore a natural outcome of the coupled GPE‑Poisson dynamics and offers a concrete, testable prediction for BEC dark‑matter models.

Performance benchmarks show that the entire 1 Gyr evolution on a 100³ grid completes in roughly 22 hours on a contemporary desktop workstation (8‑core Intel i7, 32 GB RAM). Memory consumption stays below 8 GB, making the code accessible to a wide research community. A prototype GPU implementation using CUDA reduces wall‑clock time by a factor of four, indicating that large‑scale, high‑resolution studies (e.g., 200³ or adaptive‑mesh refinements) are feasible with modest hardware upgrades.

The authors outline several extensions. Incorporating angular momentum would allow the study of rotating halos and bar formation. Adaptive mesh refinement (AMR) could resolve the dense core with higher fidelity while keeping the outer halo coarse, optimizing computational resources. Adding external fields, finite‑temperature effects, or coupling to baryonic components would broaden the model’s applicability to realistic galaxy formation scenarios. Finally, they propose a Bayesian inference pipeline that matches simulated rotation curves and lensing profiles to observational data, thereby constraining the particle mass m and self‑interaction strength g.

In summary, this work delivers a robust, efficient, and physically transparent simulation tool for self‑gravitating Bose‑Einstein condensates. By successfully integrating Crank‑Nicholson time stepping with a relaxed Poisson solver, the authors demonstrate that long‑term (∼10⁹ yr) evolution of a galactic‑scale BEC can be captured on a desktop computer. The resulting core‑halo structures provide compelling support for BEC dark‑matter hypotheses and open new avenues for confronting quantum‑fluid dark matter models with astronomical observations.