Convergence tests of self-interacting dark matter simulations

Self-interacting dark matter (SIDM) theory predicts that dark matter halos experience core collapse, a process where the halo's inner region rapidly increases in density and decreases in size. The N-body simulations used to study this process can suffer from numerical errors when simulation parameters are selected incorrectly. Optimal choices for simulation parameters are well studied for cold dark matter, but are not deeply understood when self-interactions are included. In order to perform reliable N-body simulations and model core-collapse accurately we must understand the potential numerical errors, how to diagnose them, and what parameter selections must be made to reduce them. We use the arepo N-body code to perform convergence tests of core-collapsing SIDM halos across a range of halo concentrations and SIDM cross sections, and quantify potential numerical issues related to mass resolution, time step size, and gravitational softening length. Our tests discover that halos with fewer than ${10}^{5}$ simulation particles, a resolution typically not met by subhalos in N-body simulations, suffer from significant discreteness noise that leads to variation and extreme outliers in the collapse rate. At our lowest resolution of $N={10}^{4}$ particles, this collapse time variation can reach as high as 20%. At this low resolution we also find a bias in collapse times and a small number of extreme outliers. Additionally, we find that simulations which run far beyond the age of the Universe, which have been used to calibrate SIDM gravothermal fluid models in previous work, have a sensitivity to the time step size that is not present in shorter simulations or simulations using only cold dark matter. Our work shows that choices of simulation parameters that yield converged results for some halo masses and SIDM models do not necessarily yield convergence for others.