Accurately simulating core-collapse self-interacting dark matter halos

TL;DR

This work tackles the numerical challenges of simulating gravothermal collapse in SIDM halos, focusing on both isolated and satellite halos within a Milky Way–like external potential. Using the OpenGadget3 SIDM module with adaptive kernel sizing and careful time-stepping, the authors quantify how energy conservation, kernel size, resolution, and IC sampling influence deep-collapse dynamics, and they show that a King model can describe inner density profiles during collapse. They provide a high-resolution benchmark at $N=5\times 10^7$ particles and derive practical guidelines (e.g., keep energy errors $<1\%$, use $r_{\mathrm{cut}}=15 r_s$, constrain $h/l$, and be cautious with minimum time steps) to reliably simulate SIDM halos. The results illuminate how SIDM physics can generate compact substructures relevant to GD-1 perturbations and strong lensing, and they offer concrete benchmarks and modeling approaches for future studies of gravothermal collapse in SIDM systems.

Abstract

The properties of satellite halos provide a promising probe for dark matter (DM) physics. Observations have motivated current efforts to explain surprisingly compact DM halos. If DM is not collisionless, but has strong self-interactions, halos can undergo gravothermal collapse, leading to higher densities in the central region of the halo. However, it is challenging to model this collapse phase from first principles. To improve on this, we sought to better understand the numerical challenges and convergence properties of self-interacting dark matter (SIDM) N-body simulations in the collapse phase. Especially, our aim was to better understand the evolution of satellite halos. To do so, we ran SIDM N-body simulations of a low-mass halo in isolation and within an external gravitational potential. The simulation set-up was motivated by the perturber of the stellar stream GD-1. We find that the halo evolution is very sensitive to energy conservation errors, and a SIDM kernel size that is too large can artificially speed up the collapse. Moreover, we demonstrate that the King model can describe the density profile at small radii for the late stages that we have simulated. Furthermore, for our most highly resolved simulation (N = 5x10^7) we have made the data public. It can serve as a benchmark. Overall, we find that the current numerical methods do not suffer from convergence problems in the late collapse phase and provide guidance on how to choose numerical parameters, for example that the energy conservation error is better kept well below 1%. This allows simulations to be run of halos that become concentrated enough to explain observations of GD-1-like stellar streams or strong gravitational lensing systems.

Figures (25)

• Figure 1: Interplay of softening length and adaptive time steps for energy conservation. We show the time evolution of the mass enclosed within 10 pc (upper panel) and the energy conservation (lower panel) for simulations of collisionless DM. The shaded regions indicate the uncertainties estimated based on shot noise. The simulations employ an adaptive time step, except for the results given by the yellow curve. More details can be found in Table \ref{['tab:sim_para']}; the corresponding names are A, B, C, D, E, and F.
• Figure 2: Variation of the maximum sampling radius for the ICs. Analogously to Fig. \ref{['fig:soft_tstep']}, the mass enclosed within 10 pc (upper panel) and the energy conservation (lower panel) are shown as a function of time. The results for ICs sampled up to $10 \, r_\mathrm{s}$ (orange and light blue) and $15 \, r_\mathrm{s}$ (purple and dark blue) are shown. We note that the lower resolution simulations (orange and purple) employ a softening length of $\epsilon=6.4$ pc, whereas the higher resolution simulations (light and dark blue) employ a softening length of $\epsilon=1.0$ pc. All parameters for the shown simulations (G, J, M, and O) are given in Table \ref{['tab:sim_para']}.
• Figure 3: Variation of resolution and gravitational softening length. Similar to Fig. \ref{['fig:soft_tstep']}, we show the enclosed mass within 10 pc and the energy conservation as a function of time. The simulations J, N, O, and Q as given in Table \ref{['tab:sim_para']} are shown.
• Figure 4: Impact of energy conservation on the evolution time of the halo. Following Fig. \ref{['fig:soft_tstep']}, we show the enclosed mass within 10 pc (upper panel) and the energy conservation (lower panel) as a function of time. The shown simulations differ in softening length, causing differences in the accuracy of energy conservation (see explanations in Sect. \ref{['sec:soft_tstep']}). All parameters for the displayed simulations J and L are given in Table \ref{['tab:sim_para']}.
• Figure 5: Evolution of an isolated halo as a function of time for two different choices of $N_\mathrm{ngb}$. The upper panel gives the mass enclosed within $10 \, \mathrm{pc}$ and the lower panel displays the energy conservation as a function of time. The choice of $N_\mathrm{ngb} = 384$ corresponds to twice the kernel size compared to the SIDM computations for $N_\mathrm{ngb} = 48$. Table \ref{['tab:sim_para']} gives the parameters employed for the displayed simulations J and K.