Gross-Pitaevskii-Poisson Simulations Overview

• Gross-Pitaevskii-Poisson simulations are numerical frameworks that couple the Gross-Pitaevskii and Poisson equations to model self-gravitating quantum fluids in contexts like dark matter halos and ultracold plasmas.
• They employ advanced discretization and time integration methods—such as Crank-Nicolson and ADI schemes—to ensure stability, convergence, and physical fidelity in dynamic evolution.
• These simulations enable practical studies of soliton dynamics, collapse thresholds, and stability analyses, driving insights into both astrophysical and laboratory quantum fluid phenomena.

Gross-Pitaevskii-Poisson (GPP) simulations comprise a numerical framework for modeling self-gravitating quantum fluids whose dynamics are governed by the Gross-Pitaevskii equation (GPE) coupled to the Poisson equation for the gravitational, electrostatic, or trapping potential. These systems are broadly utilized to paper dark matter halos as Bose-Einstein condensates, ultracold plasmas, and ground-state nonlinear structures in quantum matter. The central technical challenge is solving the nonlinear, nonlocal system—where the condensate wavefunction evolves under both self-interaction and a nontrivial mean-field potential sourced by its density—using explicit numerical time integration and relaxation methods for the potential, with attention to stability, convergence, and physical fidelity (Madarassy et al., 2012, Tena-Contreras et al., 4 Apr 2025, Sakaguchi et al., 2020).

1. Fundamental Equations and Physical Regimes

The prototypical GPP system consists of the nonlinear Schrödinger equation for a condensate wavefunction $\psi(\mathbf{x}, t)$ and the Poisson equation for the potential:

• Gross–Pitaevskii equation (dimensionless form):

$$(i-\gamma)\,\partial_t \psi = \left[-\frac{1}{2}\nabla^2 + V\right]\psi$$

where $V = \phi + \kappa |\psi|^2 - \mu - \Omega L_z$ encodes gravitational (or electrostatic) potential $\phi$, nonlinear self-interaction $\kappa$, chemical potential $\mu$, and rotation $\Omega$ ($L_z = i(x\partial_y - y\partial_x)$).

• Poisson equation:

$\nabla^2 \phi = 4\pi G\,|\psi|^2$

with $G$ the (dimensionless) coupling constant. For ultracold plasma models, the system generalizes to two wavefunctions ($\psi_+, \psi_-$ for cations and anions) coupled by charge and interspecies interactions (Sakaguchi et al., 2020).

The physical regime (e.g., quantum dark matter, trapped bosonic ions) dictates the coupling terms ($\kappa$, $G$), dimensionality, and symmetry (spherical, Cartesian, or multi-component). Nondimensionalization sets units via choices such as $\hbar=1$, $m=1$, and a characteristic length, enabling direct simulation on computational grids (Madarassy et al., 2012).

2. Spatial and Temporal Discretization Schemes

Spatial discretization employs uniform Cartesian (or occasionally radial) grids, with the Laplacian operator $\nabla^2$ approximated by second-order central differences. Grid sizes span $N_x\times N_y\times N_z$, with typical values from $8^3$ up to $120^3$ in three dimensions; for 1D/2D plasma models, $1024$ points in 1D and $256^2$ in 2D are standard (Madarassy et al., 2012, Sakaguchi et al., 2020). Boundary conditions are problem-dependent: Dirichlet ($\phi=0$ on all box faces) for gravitational potential, open/implicit for the quantum field.

Time integration of the GPE leverages the unconditionally stable Crank–Nicolson (CN) method in Cayley form:

$$(1 + \tfrac{1}{2}i\Delta t\,\hat{H})\,\psi^{n+1} = (1 - \tfrac{1}{2}i\Delta t\,\hat{H})\,\psi^n$$

where $\hat{H}$ is split via alternating-direction implicit (ADI) steps across Cartesian axes, each requiring solvable tridiagonal matrices for each line/plane. Nonlinear terms $\kappa|\psi|^2$ are treated using rapid fixed-point iterations at each substep (Madarassy et al., 2012).

Time evolution codes for spherically symmetric problems employ explicit fourth-order Runge–Kutta integration in time and second-order centered differences in space, with the Poisson equation solved at each RK substep by outward integration from the origin (Tena-Contreras et al., 4 Apr 2025).

3. Efficient Potential Solvers and Relaxation

Solving the Poisson equation for $\phi$ forms a computational bottleneck in GPP simulations. In three dimensions, relaxation methods akin to Gauss-Seidel update the potential via local matrix inversion:

$$\phi_{i,j,k}^{(k+1)} = \ldots \left[\mathrm{neighbor~sums} - 4\pi G |\psi|^2_{i,j,k}\right]$$

with iterative cycling until $\|\phi^{(k+1)}-\phi^{(k)}\|_\infty<\epsilon$ (typically $\epsilon=10^{-6}$). A key efficiency strategy is initializing each time-step’s relaxation with the previous potential field, exploiting slowly varying densities for rapid convergence (Madarassy et al., 2012).

In lower dimensions, Poisson solvers employ direct Fourier inversion (1D, 2D) or Green’s function convolution for the source terms—well-suited for periodic or Neumann boundary conditions and adaptable to multi-component charge distributions (Sakaguchi et al., 2020).

4. Stationary Solutions, Genetic Algorithms, and Empirical Ansätze

Construction of stationary (ground-state) solutions involves a single-parameter shooting or, more robustly, genetic algorithm (GA)-based optimization of eigenvalue problems. For spherically symmetric systems:

$$-i\frac{1}{2r^2}\frac{d}{dr}\left(r^2\psi'\right) + V\psi + g\psi^3 = \omega\psi,\quad \frac{1}{r^2}\frac{d}{dr}(r^2 V') = \psi^2$$

The ground-state is mapped via the ansatz $\Psi(t,r) = \psi(r)e^{-i\omega t}$, with boundary values tuned by a GA that optimizes $$f(\omega,r_\mathrm{max}) = 1/[\psi(r_\mathrm{max})^2 + \psi'(r_\mathrm{max})^2]$$. Population-based crossover, local/differential mutation, and elitist selection drive the search, converging when errors reach $\sim 10^{-8}$ (Tena-Contreras et al., 4 Apr 2025).

Empirical core profile formulae for the density,

$$\rho_\mathrm{core}(r) = \rho_c\left[1+(2^{1/8}-1)\left(\frac{r}{r_c}\right)^{2+\beta}\right]^{-8}$$

are parameterized by core amplitude $\rho_c$, scale $r_c$, and exponent $\beta$, fitted by GA to match numerical solutions (with best-fit coefficients $a_1$, $a_2$, $b_1$, $b_2$, $b_3$ for different interaction strengths $g$ and central amplitudes $[1,2]$). The fit is accurate to $\sim1\%$ over the relevant parameter range (Tena-Contreras et al., 4 Apr 2025).

5. Stability Analyses and Dynamic Evolution

GPP simulations validate mass, energy, and dynamical stability through long-run mass-conservation (within a few percent after $5\times10^4$ iterations on $100^3$ grids) and by monitoring unitary evolution under nonlinearity. Operator-splitting errors from ADI are $O(\Delta t^3)$ and controlled by small time steps ($\Delta t=0.01$ in code units $\simeq 2.2\times10^4$ yr) (Madarassy et al., 2012). Grid robustness is substantiated by consistent behavior across coarse (as small as $8^3$) and fine grids.

Empirical stability criteria relate collapse to parameters:

$$\alpha = g\,\rho_c^{1/2},\quad \alpha_c\approx -0.7140$$

where $\alpha<\alpha_c$ signals instability (collapse) in the attractive regime, and $\alpha>\alpha_c$ ensures oscillatory stability (Tena-Contreras et al., 4 Apr 2025). This criterion is consistent for both numerical and empirical initial conditions and ties collapse thresholds directly to model parameters.

In multi-component GPP plasmas, the modulational instability (MI) analysis of uniform mixed states yields explicit thresholds:

$u_0^2 > \frac{8\pi q^2}{G - g}$

where density waves emerge as ground states, predicted by both MI and variational approximation (VA) theory (Sakaguchi et al., 2020).

6. Applications, Model Extensions, and Soliton Phenomenology

GPP simulations are widely used in dark matter galactic halo studies (self-gravitating BEC), ultracold plasma modeling, and nonlinear wave propagation. Sakaguchi & Malomed demonstrated stable spatially periodic density waves and a rich taxonomy of bright solitons (neutral, dipole, quadrupole) in 1D/2D plasmas, with existence and energy minima verified via both analytical VA and numerical imaginary-time propagation (Sakaguchi et al., 2020).

Three distinct soliton types arise:

• Neutral: Fully mixed, exact sech profile, stable if $|g|>G$.
• Dipole: Density splitting under interspecies repulsion, with interaction-driven separation quantified analytically.
• Quadrupole: Further symmetry breaking with central and side-lobe structures, stabilized at higher $G$.

A (G,N) phase diagram delineates the ground-state type, and collision simulations reveal elastic dipole-dipole scattering and quadrupole formation via dipole-antidipole merging, consistent with the underlying field theory.

7. Numerical Challenges, Remedies, and Benchmarks

The relaxation method for Poisson's equation can be inefficient for abrupt density changes; use of previous solutions as initial guesses expedites convergence. Fixed-point iterations for nonlinear terms must be sufficiently converged at each ADI split to guarantee unitarity. Consistent implementation of boundary conditions for all fields suppresses artificial reflections and preserves physical fidelity. Time-step, grid spacing, and operator-splitting parameters must be judiciously chosen to balance computational cost, truncation error, and dynamic accuracy, as demonstrated by stability on both large ($100^3$–$120^3$) and small ($8^3$–$20^3$) grids (Madarassy et al., 2012).

Published codes based on these methods are capable of advancing 3D self-gravitating BEC simulations over timescales $\sim10^9$ years in under a day on contemporary desktops, evidencing both robustness and efficiency for astrophysical and condensed matter contexts.

---

These methodologies and results collectively establish Gross-Pitaevskii-Poisson simulations as a technically mature framework for exploring nonlinear, self-consistent quantum fluid systems, with proven fidelity to both analytic structures and empirical benchmarks across astrophysical and laboratory regimes (Madarassy et al., 2012, Tena-Contreras et al., 4 Apr 2025, Sakaguchi et al., 2020).