© The Authors 2025

1 Introduction

A plethora of observations including the clustering of galaxies and their rotation curves can be explained by assuming an invisible matter component, dark matter (DM). Within the cosmological standard model of Lambda cold dark matter (ΛCDM), this matter component mainly interacts gravitationally. Numerous experiments have been conducted to discover DM via non-gravitational interactions with standard model particles. However, unambiguous evidence for such interactions is missing (e.g. Cirelli et al. 2024). Nevertheless, there are hints from astronomical observations that DM may have self-interactions, thus altering the evolution and structure of DM halos.

Self-interacting dark matter (SIDM) has historically been introduced to solve problems on galactic scales. In particular, it was meant to reduce the DM densities in the inner region of halos and to lower the abundance of satellites, i.e. suppress substructure (Spergel & Steinhardt 2000). Even if those problems can be solved with baryonic physics and improved comparison to observations (for reviews on small-scale problems of ΛCDM see Bullock & Boylan-Kolchin 2017; Sales et al. 2022), the underlying problem of understanding the nature of DM remains and poses one of the largest challenges in contemporary physics, thus making it imperative to search for signatures of DM physics.

There have been extensive studies on the impact of DM self-interactions on various astrophysical systems. This, for example, includes merging galaxy clusters (e.g. Robertson et al. 2017; Fischer et al. 2023; Sabarish et al. 2024; Valdarnini 2024) or dwarf galaxies (e.g. Vogelsberger et al. 2014; Fry et al. 2015; Ren et al. 2019; Ebisu et al. 2022; Mancera Piña et al. 2024; Kong & Yu 2025). Moreover, these systems have been used to constrain the strength of DM scattering as a function of the relative velocity of the particles (Kaplinghat et al. 2016; Yang et al. 2024). For a comprehensive overview, we refer to the review articles by Tulin & Yu (2018); Adhikari et al. (2022).

The impact of DM self-interactions on the evolution of a DM halo can be described with the help of heat conduction (Balberg et al. 2002). The scatterings among the DM particles lead to an energy transfer that follows the velocity dispersion gradient of the DM. This implies that a DM halo may experience heat conduction inwards to its centre when a positive velocity dispersion gradient is present. This is, for example, the case for a Navarro–Frenk–White (NFW) profile (Navarro et al. 1996). As a consequence of the heat flow, the central density decreases.

When heat conduction continues, the velocity dispersion gradient becomes negative at all radii. This means that heat only flows outwards. While the central region of the halo is losing energy, it becomes more compact and the density increases. This phase of the halo evolution has been named the collapse phase. Here, the self-interactions lead to a runaway process resulting in gravothermal catastrophe (Lynden-Bell & Wood 1968; Burkert 2000).

Interestingly, some observations point towards more compact substructures, which could potentially be explained by the gravothermal collapse of SIDM halos (e.g. Turner et al. 2021; Yang & Yu 2021; Yang et al. 2023a; Gad-Nasr et al. 2024; Dutra et al. 2024; Ragagnin et al. 2024; Shah & Adhikari 2024; Zeng et al. 2025b). This includes gravitational lensing (e.g. Meneghetti et al. 2020; Granata et al. 2022), in particular the observation of satellites acting as strong lens perturbers (e.g. Vegetti et al. 2010; Enzi et al. 2025; Despali et al. 2025; Cao et al. 2025; Li et al. 2025; Minor 2025). In particular, Nadler et al. (2023) found that the projected logarithmic density profile slope of their SIDM simulations for the collapse phase is consistent with the observations by Minor et al. (2021). While gravitational lensing is very powerful in detecting substructures (e.g. Gilman et al. 2021; Keeley et al. 2024; Gannon et al. 2025), especially if they are dark, there are other possibilities to learn about substructures.

One of them is via their impact on stellar streams (e.g. Lu et al. 2025). For example, the stellar stream GD-1, which is one of the coldest and longest streams of the Milky Way, has a rich morphology (e.g. Grillmair & Dionatos 2006; Tavangar & PriceWhelan 2025). It shows gaps (e.g. Carlberg & Grillmair 2013; de Boer et al. 2018; de Boer et al. 2020; Banik et al. 2021; Malhan et al. 2022) and a spur (e.g. Price-Whelan & Bonaca 2018; Bonaca & Hogg 2018; Bonaca et al. 2020). These structures can be explained by compact perturbers passing by the stream. Depending on the DM model, we expect compact dark substructures to be more or less abundant. Zhang et al. (2025) showed that SIDM can form objects compact enough to explain the gaps and spur-like feature in GD-1, while CDM might be able to only explain the gaps.

The abundance and compactness of the substructures leave their distinct impact on stellar streams. Predictions for GD-1-like streams in the scenario of collisionless DM have been derived by Adams et al. (2025). Additionally, self-interactions could directly impact the satellite from which the stellar stream originates. This is particularly relevant if the satellite is a dwarf galaxy in contrast to a star cluster. The self-interaction of the satellite’s DM halo can potentially have at least a small impact on the morphology of the stellar stream (e.g. Zeng et al. 2025a; Hainje et al. 2025). For the future, we expect the wealth of observational data on stellar streams to increase strongly, not only for streams of the Milky Way (Bonaca & Price-Whelan 2025), but also for extragalactic stellar streams (e.g. Pearson et al. 2022). This opens up an opportunity to constrain DM models, especially SIDM and eventually find evidence for non-gravitational DM interactions.

Modelling DM substructures from first principles, in particular gravothermally collapsing SIDM satellites, is challenging (e.g. Yang & Yu 2022; Zhong et al. 2023; Mace et al. 2024; Palubski et al. 2024). Problems leading to energy non-conservation in N-body simulations of the late collapse phase of isolated halos have been explained in detail in Fischer et al. (2024a). However, a small error in energy conservation does not necessarily imply that the simulation results are accurate. Checking the conservation of total energy can only serve as a diagnostic, but it does not guarantee reliable results. Instead, it is necessary to run a simulation with high resolution and parameters tuned to achieve highly accurate results to obtain a benchmark against which other simulations can be tested.

Several aspects concerning the accuracy of simulations of the gravothermal collapse, such as their convergence behaviour, have not yet been fully investigated. Therefore, for this paper, we studied the numerical properties of such SIDM simulations using the N-body code OPENGADGET3. In addition, we studied the qualitative differences arising from the velocity and angular dependence of the SIDM cross-section for a halo in isolation and as a satellite. For the latter, our set-up was motivated by the GD-1 perturber. As a result, we comment on how to choose numerical parameters for SIDM simulations. Moreover, we show that continuing an SIDM simulation with a too large time step may lead to incorrect predictions on the enclosed mass. Finally, we demonstrate that a King model provides a good description of the density profile at small radii in the collapse phase, and our results can be applied to study strong lensing perturbers.

The remainder of the paper is structured as follows. In Sect. 2 we describe the numerical set-up of this study. The results of isolated halo simulations are presented in Sect. 3, followed by Sect. 4 with a presentation of the findings when the halo is evolved in an external potential. In Sect. 5 we discuss the limitations of our simulations as well as directions to continue this research. Finally, we summarise and conclude in Sect. 6. Additional information is provided in the appendices.

2 Numerical set-up

In this section we explain the employed numerical set-up. We begin with the simulation code and continue with the initial conditions (ICs) and finally summarise the simulation parameters.

For our simulations we used the cosmological N-body code OPENGADGET3. It is a successor of the GADGET-2 code (Springel 2005). The domain decomposition and the neighbour search we used have been described by Ragagnin et al. (2016). Additional information on the code is also given by Groth et al. (2023) and the references therein.

OPENGADGET3 contains a module to simulate DM self-interactions, introduced by Fischer et al. (2021a,b, 2022, 2024b) and has been used for numerous studies of SIDM. A unique feature of the SIDM module is that it is capable of simulating very anisotropic cross-sections (see also Fischer et al. 2021a; Arido et al. 2025). Moreover, it allows one to perform the computations in parallel employing a message passing interface (MPI) parallelisation and/or an open multi-processing (OpenMP) parallelisation. At the same time, it explicitly conserves energy and linear momentum. A separate time step criterion was implemented for the self-interactions, which ensures small-enough time steps. For computing the interactions between numerical particles, a spline kernel was employed (Monaghan & Lattanzio 1985). The kernel size for each particle was chosen adaptively and set by the next N_ngb = 48 neighbours if not stated otherwise.

We chose the ICs for our simulations motivated by the perturber of the stellar stream GD-1. It should be a halo that is a reasonable progenitor evolving to a state potentially explaining the observations. In this sense, a high concentration is desirable, as it makes rapid gravothermal evolution driven by SIDM more plausible. Our ICs follow an NFW profile (Navarro et al. 1996): $$\rho (r)=\frac{{\rho }_0}{\frac{r}{r_{\mathrm{s}}}{(1+\frac{r}{r_{\mathrm{s}}})}^2}.$$(1)

We set the density parameter to ρ₀ = 4.42 × 10⁷ M_⊙ kpc⁻³ and the scale radius to r_s = 1.28 kpc. The halo has a mass of about ≈2.8 × 10⁹ M_⊙ and its concentration parameter is c ≈ 22, which is slightly higher than the typical concentration at this halo mass (Dutton & Macciò 2014). For our simulations, we generated the ICs using SPHERIC (Garrison-Kimmel et al. 2013) and sampled the halo up to a radius r_cut = 10 r_s or r_cut = 15 r_s. The ICs are all without baryons, i.e. DM only.

We set a spherical boundary condition at 200 kpc to avoid problems arising from particles travelling very far from the halo. The boundary condition reflects particles that are travelling outwards, i.e. it does not directly affect energy conservation (as done by Fischer et al. 2024a). The self-interactions locally thermalise the velocity distribution and consequently populate the high-velocity tail of the Maxwell-Boltzmann distribution. This leads to a few particles escaping the halo. When those particles travel to large radii, this can lead to significant round-off errors in the gravity calculations, given that we store the positions of the particles using double-precision floating-point numbers. In our case, the precision of the particles may not be of much concern, given that we do not run our simulations for a very long time, and thus the particles may not travel to very large distances. However, to ensure that no problem arises from a reduced precision in the position of the particles, we introduced the boundary condition.

Moreover, we simulated the halo with three different mass resolutions and also varied the gravitational softening length. The most highly resolved simulations contain N = 5 × 10⁷ particles. Furthermore, we employed different values for the gravitational time step parameter η (Eq. (34) by Springel 2005) or the SIDM time step parameter τ (Eq. (5) by Fischer et al. 2024b) and employed for some simulations a minimum allowed time step Δt_min. In addition, we set the opening angle for the gravitational force computations by setting the value of the parameter α = 5 × 10⁻⁴ (Eq. (18) by Springel 2005). An overview with the details of our simulations is provided in Table A.1.

For most of our simulations we assumed an isotropic velocity-independent cross-section. Simulations with a different angular or velocity dependence are described in Sect. 3.8. In addition to the isolated halo simulations, we placed the halo in an external potential. This allowed us to study its evolution when being subject to tidal forces. The details are described in Sect. 4.1. An overview of the 29 simulations that we show in this paper is given in Appendix A. The choices of the various numerical parameters, such as the softening length or the values for the time step criteria, are highlighted.

3 Isolated halo evolution

In this section, we present our simulation results and discuss the relevance of various numerical parameters for achieving accurate simulations. This includes various aspects such as the role of variable time steps, the radius up to which the ICs are sampled, the role of mass resolution, the gravitational softening length and the SIDM kernel size, as well as setting a minimal time step. We also present results for different angular and velocity dependences and comment on their qualitative differences. Furthermore, we explore the evolution of the velocity anisotropy. In addition, we discuss how much the evolution of an SIDM halo simulation depends on the choice for the gravitational time step criterion in Appendix B. We also demonstrate in Appendix D that the inner region during the collapse phase can be well described by a King model (King 1962).

For all simulations, we used the peak find algorithm introduced by Fischer et al. (2021b) to find the centre of the halo. It corresponds to the minimum of the gravitational potential. We computed physical quantities such as the enclosed mass with respect to this centre.

Fig. 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 A.1; the corresponding names are A, B, C, D, E, and F.

3.1 Relation between softening length and variable time steps

Here, we investigate the size of the gravitational softening length and its relation to the use of an adaptive time-stepping scheme. For this purpose, we consider only CDM simulations.

In Fig. 1, we can see that a too small gravitational softening length can lead to a loss of energy over the course of the simulation (compare light green and light blue). However, when fixing the adaptive time steps to a constant value for all particles over the whole simulation (yellow curve), the error in energy conservation vanishes. This can be understood by a too small softening length, which causes more frequent changes in time steps, thus implying a larger error in total energy. The energy error related to the time steps is a consequence of choosing the time steps in a non-time-symmetric way (Dehnen 2017).

We note that we have chosen a fairly small time step for the simulation with the fixed time step. This ensures that all particles are evolved on a sufficiently small time step, but at the same time reduces integration errors that are not directly related to the change in time steps as well. If we instead choose a much larger constant time step for all particles, this leads to an increase in total energy by amplifying the errors that arise from the asymmetric tree evaluation that we do in OPENGADGET3, implying that it decreases when shrinking the opening angle for the tree nodes (see also the explanation given by Fischer et al. 2024a, at the end of their Sect. 3.2.1). In addition, we also checked that the energy conservation error for the simulation with the small softening length and the variable time step (light blue) does not arise from the error related to the asymmetric tree evaluation, i.e. we explicitly checked that decreasing the opening angle does not affect the error in energy conservation. In conclusion, the statement from the previous paragraph about the energy conservation error arising from the asymmetric change of time steps is valid.

Fig. 2 Variation of the maximum sampling radius for the ICs. Analogously to Fig. 1, 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 rs (orange and light blue) and 15 rs (purple and dark blue) are shown. We note that the lower resolution simulations (orange and purple) employ a softening length of ϵ = 6.4 pc, whereas the higher resolution simulations (light and dark blue) employ a softening length of ϵ = 1.0 pc. All parameters for the shown simulations (G, J, M, and O) are given in Table A.1.

3.2 Sampling radius of ICs

As a next step, we investigated the role of the ICs. In particular, we made an approximation that is commonly made when sampling a halo model such as the NFW profile.

Sampling the ICs for a density profile that implies an infinite mass, such as the NFW profile, requires truncating the halo at a radius r_cut. The truncation of the halo can influence its gravothermal evolution and, for example, can speed up the collapse for satellite galaxies, where the outer regions of the halo are removed due to tidal force (e.g. Nishikawa et al. 2020; Sameie et al. 2020). In contrast to this physical truncation, we investigated the role of the artificially set truncation radius for SIDM simulations.

In Fig. 2 we test the influence of r_cut on the evolution of the halo. A comparison of simulations with ICs sampled up to 10 r_s and 15 r_s reveals a significant deviation in the evolution. Although the simulations agree with each other during the core formation and early collapse phase within the estimated Poisson noise based on the number of simulation particles within 10 pc, they start to deviate at later times. Here the enclosed mass grows super-exponentially, and the relative uncertainty on the inferred mass goes down. At the same time, the difference between the simulations in terms of the enclosed mass becomes relatively large, although the deviation in the collapse time is small. Because of the super-exponential nature of the enclosed mass growth in the deep core-collapse phase, a small deviation in the collapse time results in a large difference in the mass. Depending on how accurately the evolution should be predicted, r_cut = 10 r_s could be too small, and a higher value would be preferable.

Finally, we note that we chose r_cut in units of r_s compared to the virial radius. This is because the only relevant length scale for the NFW profile is r_s, whereas the virial radius is only defined by making additional assumptions about cosmology, namely the critical density.

Fig. 3 Variation of resolution and gravitational softening length. Similar to Fig. 1, 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 A.1 are shown.

3.3 The role of resolution and softening length

We now look at the role of the numerical resolution, i.e. the number of numerical particles that are used to resolve the DM halo. At the same time, we consider two different choices for the gravitational softening length.

The results for four of our SIDM simulations are given in Fig. 3. Interestingly, from a comparison of the purple and black curves, we can see that the collapse time is increasing slightly with higher resolution. An effect that could contribute to this is gravitational two-body relaxation being stronger in simulations with lower resolution. An increased effective heat conduction due to a large SIDM kernel size compared to the mean free path can also contribute to this, as we discuss in Sect. 3.5. We note that the difference between these two runs may not be caused by inaccurate energy conservation. In fact, the simulation indicated by the purple curve is slightly gaining energy, which would slow down the evolution (see also Sect. 3.4). However, this is in contrast to what we find when we compare it with the black curve.

Moreover, we do not find clear evidence that the gravitational softening length has a significant impact on the evolution for the values that we have tested. This suggests that the results are not very sensitive to gravitational softening as long as the softening length is chosen within a reasonable range. The very slight difference between the green and blue curves for the enclosed mass might be related to differences in energy conservation.

A common choice is to set the gravitational softening length according to Power et al. (2003) or van den Bosch & Ogiya (2018). However, the optimal gravitational softening length for an SIDM halo in the collapse phase might be smaller than that for the pre-collapsing halo (e.g. an NFW halo). This is especially the case if one is interested in the inner regions of the collapsing object, i.e. fairly small length scales. As a consequence, as a compromise one may eventually want to choose a softening length that is somewhat smaller, closer to undersoftening than recommended according to the criteria based on CDM simulations.

Fig. 4 Impact of energy conservation on the evolution time of the halo. Following Fig. 1, 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. 3.1). All parameters for the displayed simulations J and L are given in Table A.1.

3.4 The role of energy conservation

Several numerical difficulties can affect the conservation of the total energy over time. Next, we comment on its role in simulations of SIDM halos.

The gravothermal evolution leads to a collapse because the SIDM halo is losing energy in its central region, making DM sink further inside. Similarly, the gravothermal evolution can be impacted by numerical errors in the conservation of total energy. This may imply that a loss of total energy speeds up the collapse while an increase slows it down.

In Fig. 4, we show two simulations that differ in energy conservation. The one with the smaller softening length (red) loses energy in line with the explanation given in Sect. 3.1. It loses approximately 0.6% of its total energy, resulting in a speed-up of the evolution of about 4.2%, when compared to the simulation with the larger softening length (purple). We note that the larger softening length is not directly responsible for a significant slowdown of the collapse for the corresponding simulation (purple), as we found in Sect. 3.3. Overall, this illustrates that an accurate prediction of the collapse time requires a fairly accurate conservation of total energy.

Fig. 5 Evolution of an isolated halo as a function of time for two different choices of Nngb. The upper panel gives the mass enclosed within 10 pc and the lower panel displays the energy conservation as a function of time. The choice of Nngb = 384 corresponds to twice the kernel size compared to the SIDM computations for Nngb = 48. Table A.1 gives the parameters employed for the displayed simulations J and K.

3.5 SIDM kernel size

It may have been Koda & Shapiro (2011) who first pointed out that SIDM N-body simulations should resolve the mean free path set by the local density and cross-section. Fischer & Sagunski (2024) found for a set-up rather different from a halo simulation that the ratio of the SIDM kernel size, h, to the mean free path, l, actually matters. Here we investigate the relevance of this for simulating the collapse phase of SIDM halos.

For this purpose, we reran one of our simulations, increasing the neighbour number from N_ngb = 48 to N_ngb = 384. This implies an increase in kernel size by a factor of two. In Fig. 5, we show the results for the mass enclosed within 10 pc and the energy conservation. For almost all of the evolution, there is basically no difference between the two simulations. Only at very late stages, where h/l has increased a great deal at small radii due to the collapse, do we find that the simulation with the larger value for N_ngb reaches high densities earlier. This could be explained by a larger SIDM kernel size leading to an artificially strong effective heat conduction.

We computed the ratio of kernel size, h, to mean free path, l, following Fischer et al. (2024a): $$\frac{h}{l}=\frac{{\sigma }_{\mathrm{e}\mathrm{f}\mathrm{f}}}{m}{\rho }^{2/3}\sqrt[3]{\frac{3m_{\mathrm{n}}{N}_{\mathrm{n}\mathrm{g}\mathrm{b}}}{\pi \sqrt{2}}}.$$(2)

Here, σ_eff is defined as the effective cross-section given by Eq. (C.1). In the case of a velocity-independent isotropic cross-section, σ_eff = σ. The result is given for various simulations and times as a function of radius in Fig. 6. We find that a kernel size that is even a few multiples larger than the mean free path appears to be unproblematic (compare also Fig. 3). In contrast, ratios exceeding an order of magnitude appear to be problematic. However, in general, it depends on the simulation set-up how problematic a specific value for h/l is; in particular, it should depend on the velocity dispersion gradient. In this sense, stages of the gravothermal collapse later than those that we have simulated could be even more problematic. Furthermore, we want to add that there might be a mild dependence on how the geometric factor Λ for the interaction probability and the drag force is computed (see Eqs. (9), (13) and (B4) by Fischer et al. 2021a); for example, the chosen kernel function may matter.

Finally, we note that the error arising from a too large value for h/l does not show up in the energy conservation of the simulation. Hence, this is a good example of where energy conservation alone is not enough to assess the simulation’s accuracy.

Fig. 6 SIDM kernel size divided by the mean free path as a function of radius for several times. The results, computed according to Eq. (2), are shown for several simulations with different resolutions. All simulations were run with Nngb = 48. More details for the shown simulations J, O, and Q are given in Table A.1.

3.6 Using a minimal time step

Next, we discuss how imposing a minimum time step impacts the simulation results in the collapse phase. The time step criterion for SIDM, but also gravity, requires smaller and smaller time steps while the halo is collapsing. In particular, the SIDM time step criterion becomes prohibitively small. This poses a challenge for cosmological simulations as one aims to run a simulation until a specific redshift (e.g. z = 0). A strategy to complete such a simulation, although some halos are collapsing, is to set a minimum allowed time step despite the time step criteria, i.e. evolving some particles on too large time steps.

Following this idea, we set a minimum time step for the late evolution phase. In practice, we continued some of our simulations with this minimum time step. In our case, we chose Δt_min = 2.8 × 10⁻⁵ Gyr. The results of the simulations for which we tested this are displayed in Fig. 7. Here we show the mass enclosed within 10 pc and the energy conservation error. For the late collapse phase, we find that the simulations with the minimum time step agree well with those that can contain particles on even smaller time steps. Having a minimum time step in place, we were able to continue the simulations to even later stages with higher central densities. However, at some point, the central density stopped growing, and the total energy was no longer conserved, but increased drastically.

Interestingly, the increase in total energy seems to be sensitive to the gravitational softening length and number of particles. For the three runs with the fixed minimum time step, the one with the highest number of particles and the smallest softening length actually has the most severe energy violation. The cause of this surprising behaviour is not entirely clear. Although we impose a minimum time step, the size of the time step at which a particle is evolved can still change. As a consequence, we expect a contribution from the non-time-symmetric change of time steps to the energy conservation error. We note that this error can depend on the choice of the gravitational softening length (see Sect. 3.1). Another source of the energy conservation error comes from the asymmetry of the one-sided oct-tree used in the gravitational force computations (Fischer et al. 2024a). In addition, how quickly the numerical errors grow depends on the size of the employed time step. In particular, the ratio of the actually used time step to the one implied by the time step criteria can matter, and are independent of the adopted minimum time step. Thus, although the minimum time step is fixed to a constant, the time step required by the time step criteria depends effectively on the number of particles and the softening length. The interplay of these factors complicates the evolution of the energy conservation error.

We would like to emphasise that although the runs with relatively low resolution have better energy conservation, it does not mean they are accurate. For example, the simulations shown reach densities for which the SIDM kernel sizes might be relatively large compared to the mean free path, leading to an underestimation of the collapse time, as discussed in Sect. 3.5. To efficiently and accurately simulate such late stages, more work is needed, especially in developing measures to evaluate the accuracy of the simulations, as well as alternative models to describe the inner region of the halo, as we discuss in Sect. 5.

Studying the mass enclosed only within 10 pc may give an incomplete picture of the impact of imposing a minimum time step. Instead, we consider the density, the velocity dispersion and the Knudsen number, $$Kn=\sqrt{\frac{4\pi }{\rho {\nu }^2}}{\left(\frac{{\sigma }_{\mathrm{e}\mathrm{f}\mathrm{f}}}{m}\right)}^{-1}.$$(3)

We use the effective cross-section, σ_eff, given by Eq. (C.1). Moreover, for a velocity-independent isotropic cross-section, σ_eff = σ. We computed those quantities as a function of radius and show them for two of our continued simulations in Fig. 8. For the first time that we show (4.7 Gyr), the two simulations give fairly similar results. However, at later times they deviate from each other, sometimes leading to huge differences (e.g. for the velocity dispersion). This demonstrates that measuring quantities from simulations of collapsing halos where a minimum time step had been imposed can be far off.

Nevertheless, the simulation Jt shown in purple in Fig. 8 might be more accurate as the energy conservation is much better. The results for the very late collapse phase, where we find an increase in the central density compared to the maximum core formation by a factor of 6 × 10⁵, also appear to be qualitatively in agreement with the findings in Fig. 3 by Balberg et al. (2002) for the gravothermal fluid model.

In Sect. 3.7, we comment further on the errors arising from imposing a minimum time step. In detail, we look at the enclosed mass in projection, which is relevant for gravitational lensing.

Fig. 7 Employing a minimum time step. As in Fig. 1, we show the enclosed mass within 10 pc and the energy conservation for various simulations. This time, we investigate how much the simulations are affected when employing a minimum time step. The brighter lines indicate simulations without a minimum time step, i.e. the time step can decrease further. In contrast the darker lines feature a minimum time step of Δtmin = 2.8 × 10−5 Gyr. The displayed simulations J, N, O, Jt, Nt, and Ot are described in Table A.1.

Fig. 8 Density, velocity dispersion, and Knudsen number for simulations employing a minimum time step. The purple lines give the results for the simulation with N = 2 × 106 and ϵ = 6.4 pc (grey) and Δtmin = 2.8 × 10−5 Gyr and the blue lines are for the simulation with N = 1 × 107 and ϵ = 1.0 pc and Δtmin = 2.8 × 10−5 Gyr. These are the same simulations as in Fig. 7. The top panel gives the density as a function of radius for different times. We note that these times are during the collapse phase, and for most of them the energy conservation error is sizeable. The middle panel gives the velocity dispersion as a function of radius, and the bottom panel displays the Knudsen number (Eq. (3)). All parameters for the shown simulations Jt and Ot are given in Table A.1.

Fig. 9 Time evolution of the enclosed mass for our most highly resolved simulation (simulation Q in Table A.1). The relative change of the enclosed mass within various radii is shown as a function of time. The enclosed mass is computed in projection, i.e. it is the mass within a cylinder of radius r2D.

3.7 Evolution of projected enclosed mass

We wanted to get closer to the observations, in particular gravitational lensing, and to measure the enclosed mass in projection. In Fig. 9, the enclosed mass within several radii is shown as a function of time for our most highly resolved simulation (N = 5 × 10⁷). For small radii, we find that the enclosed mass first decreases during the core expansion phase and later increases when the gravothermal collapse dominates, as we also found in our three-dimensional analysis. This is in contrast to the bound mass of the system. The self-interactions create a few particles that have velocities exceeding the escape velocity, leading to a decreasing mass of the halo. Even more particles may gain enough kinetic energy to travel to fairly large radii, but still remain bound. In line with this, we find that for the largest radius shown in Fig. 9, the mass decreases slightly. Related to this, the velocity anisotropy at large radii increases, i.e. radial motion becomes more important, as we show in Appendix G.

If one defines a radius r_eq as the one at which the mass inflow and the mass outflow have the same size, i.e. $\dot{M}(<r_{\text{eq}})=0$. Then r_eq decreases with time when the halo is undergoing the gravothermal collapse. This implies that we do not expect the mass within a fixed radius (e.g. 1 kpc) to keep continuously increasing. Instead, it starts to decrease as we can see for a few of the radii that we have chosen for Fig. 9. The smaller the radius, the later the point in time where the enclosed mass starts to decrease¹.

The projected logarithmic density profile slope, γ_2D, shows a similar behaviour. In Fig. 10, we give γ_2D as a function of the enclosed mass M_<r2D. The two quantities can be inferred from observations via gravitational lensing analysis. Here we use our most highly resolved simulation, which models the halo in isolation. The figure shows how the system evolves in the γ_2D–M_<r2D plane considering several radii, r_2D. It can be seen that not only does the enclosed mass make a turn during the collapse phase, as discussed above, but the density slope also starts to become flatter at a specific point in time. Similar to the enclosed mass, this turn occurs earlier at larger radii. Interestingly, γ_2D turns even earlier to becoming flatter than the projected enclosed mass starts to decrease.

In addition to the results above, we also show how the enclosed mass evolves for our simulations employing a minimum time step. The results for all three of them are given in Fig. 11. The figure shows that for the late phase that we would hardly be able to simulate without imposing a minimum time step; the values for the enclosed masses differ significantly between the different runs. As a result, cosmological simulations employing a minimum for the allowed time step may overestimate the collapse time and underestimate the density at the halo centre. But importantly, they may also overestimate the mass of the halo within the scale radius of the initial NFW profile or larger radii. As a consequence, the ability of SIDM simulations to explain gravitational lensing signals supposedly arising from fairly concentrated objects can be easily overestimated.

Fig. 10 Projected logarithmic density profile slope as a function of projected enclosed mass. For our most highly resolved simulation evolved in isolation (simulation Q from Table A.1), we show how the projected logarithmic density profile slope, γ2D, and the projected enclose mass, M<r2D, evolve with time. We display γ2D and M<r2D for several radii r2D (see inset for legend). The larger circles mark the value at the beginning of the simulation, i.e. for an NFW profile. Moreover, the small white dots are placed equidistant in time, every 0.98 Gyr.

Fig. 11 Enclosed mass within several radii as a function of time. The results are shown for simulations with the standard time step criterion (darker curves) and those employing a minimum time step (lighter curves). The three panels give the enclosed mass within different radii for the same simulations, as shown in Fig. 7. Moreover, we show the enclosed mass for the initial NFW halo as a reference (grey).

Fig. 12 Velocity and angular dependent cross-sections. Following Fig. 1, the enclosed mass within 10 pc (upper panel) and the energy conservation (lower panel) are shown. In detail, a velocity-independent isotropic cross-section (blue), a velocity-independent forward-dominated cross-section (light blue), a velocity-dependent isotropic cross-section (red), and a velocity-dependent forward-dominated cross-section (light red) are displayed. They all have the same resolution of N = 107 and a gravitational softening length of ϵ = 1 pc. More details on the simulations O, P, R, and S can be found in Table A.1.

3.8 Velocity and angular dependence

So far, we have discussed the simulation results for an isotropic velocity-independent cross-section only. Now we focus on the angular and velocity dependence. In contrast to the previous part, we are mainly interested in qualitative differences arising from particle physics. Given that all the simulations we discuss here are based on similar numerical schemes, they also have similar numerical properties.

We simulated cross-sections with different velocity and angular dependences. For the velocity-dependent cross-sections we used $$\frac{{\sigma }_{\mathrm{V}}}{m}=\frac{{\sigma }_0}{m}{[1+{\left(\frac{v}{w}\right)}^2]}^{-2}.$$(4)

Here, σ_V is the viscosity cross-section as given by Eq. (A.1). The velocity-dependent cross-sections we simulated are given by σ₀/m = 6593.89 cm² g⁻¹ (in terms of the normalised viscosity cross-section, Eq. (A.1)) and w = 20 km s⁻¹. All values used for our simulations are given in Table A.1. Moreover, we simulated two different angular dependences, isotropic scattering and a forward-dominated cross-section, specifically the limit of keeping the momentum transfer constant while the scattering angle approaches zero (e.g. Kahlhoefer et al. 2014; Fischer et al. 2021a).

In Fig. 12, we display the evolution of the isolated halo with different velocity and angular dependences. The angular dependences were matched using the viscosity cross-section (Eq. (A.1)). As we can see, this results in a similar time evolution for the velocity-independent and dependent cases. Moreover, the energy error is similar as well. This demonstrates that the viscosity cross-section works very well for matching different angular dependences, as previously stated by Yang & Yu (2022); Sabarish et al. (2024). In contrast to the angular dependence, the velocity dependence (Eq. (4)) leads to a qualitatively different evolution of the halo, with a larger maximum core size and a longer collapse time compared to the core expansion phase. This qualitative behaviour arising from the velocity dependence is discussed in greater detail in Fischer et al. (2024b) (Sect. 3.2.3). We provide further discussion in Appendix C.

Fig. 13 Velocity anisotropy as a function of radius. We show β (Eq. (5)) for our most highly resolved simulation (simulation Q of Table A.1) at different times.

3.9 Velocity anisotropy

As a last part on the evolution of the isolated halo, we present the velocity anisotropy, and compute it as $$\beta \equiv 1-\frac{{\sigma }_{\theta }^2+{\sigma }_{\varphi }^2}{2{\sigma }_r^2},$$(5)

with the polar velocity dispersion ${\sigma }_{\theta }^2$, the azimuthal velocity dispersion ${\sigma }_{\varphi }^2$, and the radial velocity dispersion ${\sigma }_r^2$. The results are displayed in Fig. 13. For t = 0 the velocity distribution is isotropic at all radii (β = 0), but increases at later times at the larger radii beyond the inner density core. Here, the radial motion becomes more pronounced (see also Gurian & May 2025). We note that at the largest radii shown in Fig. 13, the anisotropy increases. This might be a consequence of the particles moving radially outwards becoming more dominant as the density decreases strongly with radius (ρ ∝ r⁻³). In addition, the collision rate is low at those radii, and the self-interactions do not efficiently isotropise the velocity distribution.

4 Satellite evolution

In this section we discuss the evolution of the SIDM halo under the influence of tidal forces arising from a host galaxy. We first describe the changes to our numerical set-up in Sect. 4.1. An illustration of the satellite orbit and the tidal evolution can be found in Appendix E. We show the difference between a satellite halo and an isolated system (Sect. 4.2) as well as the difference arising from the angular and velocity dependence (Sect. 4.3). Moreover, we demonstrate in Sect. 4.4 that the density profile during the collapse phase can be well described by a King model. Overall, we mainly focus on the physics, in contrast to the numerical aspects of the previous section. Given that we tested various numerical parameters for the isolated case, we now assume that the same choice also works when evolving the halo in an external potential. We note that we used the peak find algorithm by Fischer et al. (2021b) to determine the position of the satellite. As in the previous section, we compute physical quantities with respect to this position for all of this section.

4.1 Numerical set-up

To model the tidal forces of the host system acting on the satellite at low computational costs, we described the gravitational potential of the host analytically. We focused only on the role of the tidal forces and neglected the dynamical friction that acts on the satellite and causes its orbit to decay. Moreover, we also did not take the scattering of DM particles of the satellite halo on the host’s DM into account, which can play an important role especially for velocity-independent cross-sections (e.g. Zeng et al. 2022).

The external potential is chosen to mimic the gravitational potential of the Milky Way and consists of six components. The details of our description follow that employed by Zhang et al. (2025). We describe the DM halo with an NFW profile (Navarro et al. 1996). The corresponding potential is $${\mathrm{\Phi }}_{\mathrm{N}\mathrm{F}\mathrm{W}}(r)=-4\pi \mathrm{G}\frac{{\rho }_{\mathrm{s}}{r}_{\mathrm{s}}^3}{r}\ln (1+\frac{r}{r_{\mathrm{s}}}).$$(6)

Here G denotes the gravitational constant. For our simulations, we set the NFW halo parameters to ρ_s = 8.54 × 10⁶ M_⊙ kpc⁻³ and r_s = 19.6 kpc. To mimic the gravitational potential of the stellar bulge, we employed a Hernquist profile (Hernquist 1990). $${\mathrm{\Phi }}_{\mathrm{H}\mathrm{e}\mathrm{r}\mathrm{n}}(r)=-\frac{\mathrm{G}{M}_{\mathrm{H}}}{a_{\mathrm{H}}}{(1+\frac{r}{a_{\mathrm{H}}})}^{-1}.$$(7)

For our set-up, we used M_H = 9.23 × 10⁹ M_⊙ and a_H = 1.3 kpc. In addition, we added four disk components described by an axisymmetric Miyamoto–Nagai profile (Miyamoto & Nagai 1975), $${\mathrm{\Phi }}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}(r)=-\mathrm{G}{M}_{\mathrm{d}}{[R^2+{(a_{\mathrm{d}}+\sqrt{z^2+b_{\mathrm{d}}^2})}^2]}^{-1/2},$$(8)

with $R=\sqrt{x^2+y^2}$. The components are as follows: a thin stellar disk with M_d = 3.52 × 10¹⁰ M_⊙, a_d = 2.5 kpc, and b_d = 0.3 kpc; a thick stellar disk with M_d = 1.05 × 10¹⁰ M_⊙, a_d = 3.02 kpc, and b_d = 0.9 kpc; a thin gas disk with M_d = 1.2 × 10⁹ M_⊙, a_d = 1.5 kpc, and b_d = 0.045 kpc; and last a thick gas disk with M_d = 1.1 × 10¹⁰ M_⊙, a_d = 7.0 kpc, and b_d = 0.085 kpc.

The satellite halo is initially placed at a distance of r = 123.6 kpc from the centre of the external potential. The exact coordinates are [58.16, 40.12, 101.41] (kpc). The position is relative to the origin of the external potential. The initial velocity of the satellite is v = 109 km s⁻¹. The velocity vector reads [45.94, 83.42, 53.12] (km s⁻¹). With this choice, we also follow Zhang et al. (2025). Moreover, we illustrate the orbit of the satellite halo in Appendix E.

Fig. 14 Evolution of an isolated halo (solid line) and a satellite halo (dashed line). The simulations are for a velocity-independent cross-section, employ N = 107 particles, and use a gravitational softening length of ϵ = 1 pc. All parameters for the displayed simulations O and W are given in Table A.1.

4.2 Isolated versus satellite

Next, we discuss how the evolution of the density and velocity dispersion profiles changes for a satellite halo compared to a system evolved in isolation. In Fig. 14, we show the evolution of the density, the one-dimensional velocity dispersion, and the Knudsen number, as a function of radius for a halo in isolation (solid lines) and for a satellite system (dashed lines).

By and large, the evolution of the satellite system is similar to that of the isolated halo. First, the two undergo a core expansion phase and later collapse. However, the velocity dispersion profile (middle panels of Fig. 14) reveals a significant difference. For the satellite halo, the velocity dispersion is strongly affected by tidal forces. Close to the pericentre passage (e.g. at 1.47 Gyr) tidal heating injects energy into the halo and the velocity dispersion increases. Subsequently, particles become unbound and leave the system, and the velocity dispersion decreases over all relevant radii. This is in contrast to the isolated halo, where the velocity dispersion is only affected by the self-interactions. The speed-up of the halo evolution due to tidal stripping has been found by various authors (e.g. Kahlhoefer et al. 2019; Sameie et al. 2020). This makes the energy outflow of the satellite more efficient and speeds up the halo evolution, i.e. reduces the collapse time. Later, in the collapse phase, the velocity dispersion in the centre reaches unprecedentedly high values.

We note that the satellite system reaches a central density of ρ = 3 × 10¹¹ M_⊙ kpc⁻³ in the last snapshot at t = 4.6 Gyr. At sufficiently small radii, the density gradient becomes flat and the Knudsen number reaches a value of about 0.04, i.e. the central region is deep in the short-mean-free-path regime. In contrast, the lowest central density we find at the maximum core expansion is about ρ = 6 × 10⁷ M_⊙ kpc⁻³. It is lower compared to the isolated case due to tidal interactions. Moreover, the tidal interactions also greatly reduce the density in the outer regions of the halo, which leads to a large loss in gravitationally bound mass.

4.3 Angular and velocity dependence of the cross-section

In this section, we investigate the role of the velocity and angular dependence of the self-interaction for the evolution of the satellite halo. As for the isolated halo, we employed isotropic and forward-dominated cross-sections that either follow the velocity dependence given by Eq. (4) or are velocity-independent. In detail, we employed the same cross-section as in Sect. 3.8. The results for the enclosed mass within 10 pc are shown in Fig. 15.

The tidal forces acting on the satellite halo only slightly speed up the gravothermal evolution in the case of a velocity-independent cross-section. Interestingly, we find that this is very different for our velocity-dependent cross-sections. When the halo experiences tidal forces, the collapse time reduces to about half its value compared to the evolution in isolation (compare Figs. 12 and 15).

This implies that a velocity-dependent cross-section not only speeds up the collapse by suppressing interactions between DM particles of the host and the satellite, but there is already a contribution coming from within the satellite halo. We note that our cross-section has a fairly strong velocity dependence, i.e. w is smaller than the relative velocities within the halo. This implies an increase in the effective cross-section relative to the velocity-independent case when the velocity dispersion decreases due to tidal stripping. We show this in detail in Appendix F.

Interestingly, we find that the angular dependence starts to play a role in the velocity-independent case, although the scatterings between the DM of the host and the satellite are neglected. The small-angle scattering leads here to a longer collapse time compared to the isotropic cross-section (Fig. 15). For the isolated halo, the two cross-sections resulted in the same collapse time (Fig. 12). Fischer et al. (2021a) show that the thermalisation process depends on the angular dependence of the cross-section, i.e. large-angle scattering populates the high-velocity tail of the Maxwell-Boltzmann distribution more quickly compared to small-angle scattering. In combination with the tidal forces, this allows us to explain the different collapse times. Moreover, the difference arising from the angular dependence is minor in the velocity-dependent case where scattering with a high relative velocity is strongly suppressed (bottom panel of Fig. 15). We also checked that both angular-dependent and -independent cases give rise to almost the same total gravitationally bound mass of the satellite halo regardless of the velocity dependence. For the runs shown here, the angular dependence of the scattering does not affect the global properties of the halo.

Fig. 15 Mass enclosed within 10 pc as a function of time for the satellite halos. The upper panel gives the results for a velocity-independent cross-section (simulations U, V, W, and X), and the lower panel is for velocity-dependent cross-sections (simulations Y and Z). More details on the simulations can be found in Table A.1. In addition, the grey dotted line indicates the preferred mass range for the perturber of the stellar stream GD-1 (Bonaca et al. 2019).

4.4 Fit with King model

In this last part about the satellite halo simulations, we fit the density profile of the late collapse phase with a King model (King 1962), $$\rho (r)={\rho }_0{[1+{\left(\frac{r}{r_{\mathrm{c}}}\right)}^2]}^{-3/2}.$$(9)

With this, we follow Zhang et al. (2025) who also used a King profile to describe the collapse phase of a potential GD-1 perturber, motivated by research on star clusters. To determine the density and radial parameters (ρ₀, r_c), we used a limited radial range with r < 0.5 kpc only. For the fit, we maximised a likelihood based on Poisson statistics analogously to the description in Sect. 4 by Fischer et al. (2021a).

In Fig. 16, we show the fit of the King model for the collapse phase using our simulation with a velocity-dependent isotropic cross-section (simulation Y according to Table A.1). It becomes visible that the King model provides a reasonable fit to the inner region of a satellite halo that undergoes core collapse. This is in line with the findings of Zhang et al. (2025) and may work well enough to fit the observational data, as discussed by the same authors. In addition, we fit the King mode to the isolated halo and show the density profile as well as the time evolution of the fitted parameter in Appendix D.

Fig. 16 Density profile fitted by a King model for a satellite halo. The results are for the simulation with a velocity-dependent isotropic cross-section. All simulation parameters are given in Table A.1 under simulation Y.

5 Discussion

In this section, we discuss the limitations of our work and of the employed numerical scheme. We mention the physical processes that we neglected and that would be worth including in the modelling, and we highlight directions to improve the modelling of collapsing SIDM halos.

In our simulations, we neglected the scatterings between satellite and host particles. However, they can play a crucial role in shaping the evolution of the satellite. How strongly the DM of the satellite is affected depends on how strong the cross-section is on the relevant velocity scale, which is typically larger than for the velocities within the satellite. Hence, for a cross-section that is decreasing as a function of velocity, this is less important than for a velocity-independent cross-section (e.g. Silverman et al. 2022; Zeng et al. 2022, 2025b). In addition, the angular dependence of the cross-section can play an important role. Satellite halos are strongly affected if the cross-section is forward-dominated (Fischer et al. 2022) and the effect becomes weaker for velocity-dependent cross-sections (e.g. Nadler et al. 2020; Fischer et al. 2024b). The DM self-interactions between host and satellite could be incorporated in simulations such as the ones we have presented following the approach by Zeng et al. (2022) or within the more general framework of Klemmer et al. (in prep.).

Another physical effect that we neglect here is dynamical friction. However, given the low mass of our satellite halo compared to the host mass, it is expected to be hardly affected by dynamical friction. Nevertheless, it would be possible to include the effect of dynamical friction without resolving the host system, but by approximating the deceleration experienced by the subhalo (e.g. Petts et al. 2015).

Our numerical scheme has the advantage that the pairwise interactions do not harm energy conservation, although we employ shared and distributed memory parallelisation. Many of the codes used today do not have this property. In our case, it implies that the energy conservation is independent of the chosen kernel function for the SIDM computations, by construction. This is because the value of the interaction probability or the drag force does not have a direct impact on energy conservation. However, the chosen kernel function could eventually make a difference in the context of angular momentum conservation and for the effective heat conduction in the case of large values for h/l (see Sect. 3.5).

Ways to improve the accuracy of SIDM simulations have already been discussed by Fischer et al. (2024b). We want to highlight that using a symmetric evaluation of the oct-tree used for the gravitational force computation would be helpful (e.g. as done by Appel 1985). Another aspect is the time-irreversible choice of the time step. Approaches that make the time-stepping function more time symmetric can help to improve the accuracy, for example by reducing the error on energy conservation (Dehnen 2017).

We also note that the above-mentioned strategies may not fully solve the challenges in simulating the late collapse phase. For example, it remains problematic to resolve smaller and smaller length scales while the halo is collapsing. Even if this were resolved, and all numerical issues directly related to simulating SIDM as well, the modelling of gravity would still require a time step that continues shrinking while the halo is collapsing. The time step would become so small that the simulation could no longer be reasonably fast advanced in time. One way to overcome the challenges posed by gravothermal collapse could be the introduction of a subgrid model that describes the very inner part of the halo. It might be a sink particle that accounts for the DM mass in the centre of the halo. In the post-collapse state, it may describe the black hole that finally formed and the DM spike that it is accreting. Idealised simulations of such a system can help build such a model (Sabarish et al. 2025).

6 Conclusion

In this work, we used the N-body simulation code OPENGADGET3 to investigate the gravothermal collapse of SIDM halos. We studied the role of various numerical parameters and simulated the system in isolation and as a satellite halo. Additionally, we studied the qualitative differences arising from the velocity and angular dependence of the self-interaction cross-section. Furthermore, the data for our most highly resolved simulation (N = 5 × 10⁷) are available² and may serve as a benchmark. Our main results and recommendations for SIDM simulations are as follows:

1. A very small softening length can lead to larger errors in the total energy because the time steps change more often in an asymmetric way. However, the softening length chosen should be small enough to resolve the relevant length scales (e.g. the size of the constant density core during collapse);

2. Overall, halo evolution is very sensitive to energy conservation errors;

3. We did not find that a stricter gravity time step criterion is needed for our SIDM simulations compared to what is commonly used for CDM simulations, i.e. η = 0.025 is fine. At the same time, we may use a more stringent SIDM time step criterion than some other codes. It should be noted that simulations that evolve an SIDM halo over many more gravitational relaxation timescales than we did may need a smaller value for η (see also Sect. 4 by Fischer et al. 2024a);

4. Moreover, we recommend using a value for the opening criterion as small as α = 5 × 10⁻⁴, when using a one-sided oct-tree.

5. A SIDM kernel size that is too large effectively leads to heat conduction that is too strong, which causes the halo to collapse faster. We find that h/l should be below 10, but exceeding 1 by a few multiples appears to be unproblematic for the range that we have tested. This has implications on the number of particles required to resolve the collapse phase, i.e. it depends on how deep one wants to simulate in the gravothermal collapse;

6. We also note that the maximum sampling radius of ICs has an influence on how fast the halo collapses. For a controlled simulation of an isolated NFW halo, we recommend using 15 r_s instead of 10 r_s for the truncation radius;

7. In our tests that enforce a minimum time step, we found that the central density increases further and stays constant at the same time that the total energy artificially increases. Although imposing a minimum time step in cosmological simulations may allow correct predictions of the fraction of collapsed halos to be made, we do not find evidence that the enclosed mass within a specific radius can be accurately predicted. The mass within larger radii (≳r_s) might be overestimated;

8. During the late collapse stages the velocity distribution is no longer isotropic, but becomes progressively more radial outside the inner density core.

9. The speed-up in halo evolution for a satellite compared to an isolated system is much greater for a velocity-dependent cross-section compared to a velocity-independent one, even when ignoring the scattering between the DM particles of the satellite with the host’s DM;

10. We find that a King model provides a good fit for the density profile in the inner regions of the halo undergoing gravothermal collapse at the stage that we have simulated. This is the case for the isolated halo and the satellite.

Despite the numerical challenges, we have shown that it is possible to accurately simulate SIDM halos in the deep-collapse phase to compare with observations. Taking the dense perturber of the GD-1 stellar stream as an example, our simulations can resolve the simulated halo within 10 pc, which is necessary to match the observations. It becomes prohibitively expensive to simulate extremely deep into the collapse phase using the N-body method, such as collapsing into a black hole. Nevertheless, in that regime, simulations based on schemes derived from first principles exploiting the symmetries of the problem can provide a complementary method (e.g. Gurian & May 2025; Kamionkowski et al. 2025).

Acknowledgements

We thank all participants of the Darkium SIDM Journal Club for discussion, in particular Antonio Ragagnin. MSF is also grateful to Andrew Robertson for his comments and suggestions. This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 “Origins” – 390783311. KD and MSF acknowledge support by the COMPLEX project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679. HBY acknowledges support by the U.S. Department of Energy under grant No. de-sc0008541. Software: NumPy (Harris et al. 2020), Matplotlib (Hunter 2007).

References

Appendix A Overview of simulations

With Table A.1, we provide an overview of the simulations shown in this paper. Moreover, we use a consistent colour scheme for the figures with the colours as indicated in the table. The strength of the self-interactions is specified in terms of the viscosity cross-section. It is given by $${\sigma }_{\mathrm{V}}=3\pi {\int }_{-1}^1\frac{\mathrm{d}\sigma }{\mathrm{d}{\mathrm{\Omega }}_{\mathrm{c}\mathrm{m}\mathrm{s}}}{\sin }^2{\theta }_{\mathrm{c}\mathrm{m}\mathrm{s}}\mathrm{d}\cos {\theta }_{\mathrm{c}\mathrm{m}\mathrm{s}},$$(A.1)

and normalised to match the total cross-section in the case of isotropic scattering.

Appendix B The size of the time step

Previous studies (Mace et al. 2024; Palubski et al. 2024) indicated that one may need to be careful in choosing the time step for simulations of SIDM halos to ensure accurate results. It has been claimed that depending on the set-up, simulations of the SIDM halo collapse require more stringent settings for the time step criteria than in common CDM simulations. This has been particularly expressed in terms of choosing smaller values for the time step parameter of the gravitational time step criterion.

In our simulations, we find that the value of the time step parameter η has no significant influence on the results for the range that we have tested as shown in Fig. B.1. We note that $\mathrm{\Delta }{t}_{\text{grav}}\propto \sqrt{\eta }$ (Eq. 34 by Springel 2005). In conclusion, we do not need a stricter time step criterion for gravity compared to the CDM simulations. However, we also note that we use at the same time our time step criterion for SIDM (Fischer et al. 2024b). It becomes relevant at the late stage of the evolution (t ≳ 4.4 Gyr) for the particles in the inner region. Here, it results in smaller time steps than solely based on the gravitational time step criterion. Moreover, we want to point out that it might be important how exactly the SIDM time step criterion is formulated (see appendix B by Fischer et al. 2024b).

Fig. B.1 Variation of the gravitational time step parameter η. In the upper panel, we show the mass enclosed within the inner 10 pc as a function of time. The lower panel gives the evolution of the total energy relative to the absolute initial value. The shown simulations are for the low resolution of N = 2 × 106 and share the same parameters except η, which is varied as indicated in the legend. All parameters are given in Table A.1 under the names H, I, and J. The shaded regions in the top panel indicate the uncertainty estimated based on shot noise.

Furthermore, we want to point out that we have simulated a set-up for which the impact of the self-interactions is much stronger than numerical errors arising from solving for gravity. When simulating a DM halo over many more gravitational timescales, it could be that the desired accuracy can only be achieved when choosing a smaller time step, i.e. a smaller value for η (see also the discussion about the gravitational relaxation in Sect. 4 by Fischer et al. 2024a). However, if the numerical integration errors become significantly large, the collapse time of an SIDM halo simulation may no longer be inversely proportional to the cross-section. This might be in line with the results Mace et al. (2024) obtained for their simulations of small cross-sections that they have run for long time spans.

Appendix C Effective cross-section of the isolated halo

In this appendix, we discuss the differences between the velocity-independent and dependent runs for the isolated halo, as shown in Sect. 3.8. In addition to Fig. 12, we display the mass enclosed within 0.15 kpc as a function of time in Fig. C.1. We choose a larger radius to reduce the uncertainty in the measured mass. Importantly, we rescale the time axis of the velocity-independent runs to match the collapse phase of the velocity-dependent simulations. It becomes visible that there is a qualitative difference between the simulations, given that the core-expansion phase cannot be matched simultaneously. We note that the collapse timescale of the halo does not scale inversely proportional with the cross-section. This is because the halo is partially in the short-mean-free-path regime. If we nevertheless assume this scaling, the rescaled version of the velocity-independent cases would correspond to a cross-section of σ_V/m = 49.6 cm² g⁻¹.

Fig. C.1 Enclosed mass and energy conservation as a function of time. We show the same simulations as in Fig. 12 (O, P, R, and S in Table A.1). The upper panel gives the mass within a radius of 0.15 kpc and the lower panels shows how well the total energy is conserved. The time for the simulations with a velocity-independent cross-section is rescaled to match the collapse phase of the simulations with velocity-dependent interactions.

The effective cross-section (Yang & Yu 2022), σ_eff (see Eq. (C.1)) is a common way to map the effect of SIDM models with different velocity and angular dependences on each other (see also Yang et al. 2023b). It can be expressed as $${\sigma }_{\mathrm{e}\mathrm{f}\mathrm{f}}=\frac{⟨v^5{\sigma }_{\mathrm{V}}(v)⟩}{⟨v^5⟩}.$$(C.1)

The average is computed by integrating over a Maxwell-Boltzmann distribution set by a characteristic velocity dispersion. For the σ₀/m value as in Eq. (4), we employed a value of 90% of the maximum velocity dispersion of the initial NFW halo for the characteristic velocity-dispersion, which is 17.28 km s⁻¹, and found the corresponding effective cross-section is 80 cm² g⁻¹, which is 61% larger than the value from the rescaling in Fig. C.1, as discussed above.³ For comparison, w = 20 km s⁻¹ in the velocity-dependent models. Moreover, σ_V denotes the viscosity cross-section normalised to match the total cross-section for isotropic scattering (Eq. (A.1)).

To allow for a better understanding of the difference between velocity-dependent and velocity-independent cases, we compute the effective cross-section (Eq. (C.1)) as a function of radius. Here, we do not use a single characteristic velocity dispersion, but instead use the value of the velocity dispersion at each radius from the simulation. In Fig. C.2, we show the results for the simulation with isotropic scattering (run R according to Table A.1). We can see that the effective cross-section initially is large at small and large radii but low for the intermediate range. At the early stage of core expansion, the constant cross-section 80 cm² g⁻¹ is a reasonable approximation for the effective cross-section. However, in later stages, the effective cross-section decreases a lot, which is a consequence of the increasing velocity dispersion and is in line with the long collapse time that we find in Fig. 12. This continuous decrease during the collapse may limit the usability of the effective cross-section using a single characteristic velocity dispersion to map the effect that different SIDM models have on the evolution of an isolated halo on each other. For the velocity-dependent models, we consider w = 20 km s⁻¹, which is comparable to the characteristic velocity dispersion of the initial halo 17.28 km s⁻¹. Thus, the effective cross-section is sensitive to the change of the velocity dispersion in the halo during the core collapse. For models where w is larger than the characteristic velocity dispersion, we expect that the mapping approach should work better. We will leave detailed studies on this topic for future work.

Fig. C.2 Effective cross-section as a function of radius. For the isolated halo evolved with an isotropic velocity-independent cross-section (simulation R according to Table A.1), we show the effective cross-section for several times as indicated in the legend. In addition, we display the cross-section for the velocity-independent cases (black line).

Appendix D Fit with King model

In this appendix, we fit our highest-resolved simulation for the halo in isolation (simulation Q according to Table A.1) with the King model (Eq. (9)). As for the satellite halo in Sect. 4.4, we use a limited radial range only, r < 0.5 kpc.

In Fig. D.1, we give the density profile together with the fitted King model. The King model appears to provide a reasonable description of the central region of a collapsing SIDM halo. However, it does not describe the outer regions well and may also eventually fail for the inner region at even later stages. Nevertheless, it may work well enough to fit observational data as discussed by Zhang et al. (2025). We note that the central density of the fit, ρ₀, rises super-exponentially with time, and at the same time, the mass within r_c is decreasing as a function of time as we show in more detail below.

We display the parameters of the King model (Eq. (9)) that we obtained by fitting to our highest-resolved simulation in Fig. D.2. It is visible that the central density quickly rises while the core size decreases. Moreover, we show the enclosed mass within r_c and the maximum circular velocity ($v_{\text{vcirc\hspace{0.17em}}}=\sqrt{\mathrm{G}M/r}$) as a function of time in the bottom panel. Here, we can see that the core mass is quickly decreasing during the gravothermal collapse.

Fig. D.1 Density as a function of radius. Simulation results and a fit using the King model (Eq. (9)) are shown for different times of the late collapse phase. The results are for the isolated simulation with a resolution of N = 5 × 107 particles (simulation Q of Table A.1).

Fig. D.2 Parameters of fitted King model and further quantities as a function of time. The parameters that we determine by fitting a King model (Eq. (9)) to our highest-resolved simulation (Q in Table A.1) are shown in the top panel as a function of time. The bottom panel gives the mass enclosed in rc according to the fitted model as well as the maximum circular velocity directly computed from the simulation data.

Appendix E Orbit of the satellite and tidal evolution

To give an overview, we show the trajectory of the satellite in Fig. E.1. The results for a CDM and an SIDM simulation are shown. The two orbits are almost identical. This is also related to the fact that we neglect dynamical friction and the non-gravitational interaction between the satellite’s DM particle with the host’s DM.

Fig. E.1 Trajectory of the satellite. The three-dimensional position of the subhalo for a CDM and an SIDM simulation is shown. The corresponding simulations are T and W according to Table A.1, where their parameters are given.

Fig. E.2 Tidal radius and its enclosed mass as a function of time. For the satellite halo, evolved with different DM models, the tidal radius (solid lines) according to Eq. (E.1) and the mass within it (dashed lines) are shown. All simulations shown here share a resolution of N = 107 particles and a gravitational softening length of ϵ = 1.0 kpc. Further information on the simulations T, W, X, Y, and Z is provided in Table A.1.

Next, we investigate the tidal evolution of the satellite for the DM models that we have simulated. Here, we compute the tidal radius, r_tidal, of the satellite halo. For this, we implicitly define it as $$\frac{\mathrm{G}{M}_{\text{sat}}(<r_{\text{tidal}})}{r_{\text{tidal\hspace{0.17em}}}^2}=|\mathrm{\nabla }{\mathrm{\Phi }}_{\text{host\hspace{0.17em}}}(\mathbf{x}+r_{\text{tidal\hspace{0.17em}}}\mathbf{n})-\mathrm{\nabla }{\mathrm{\Phi }}_{\text{host}}(\mathbf{x}-r_{\text{tidal\hspace{0.17em}}}\mathbf{n})|.$$(E.1)

Here, M_sat(< r_tidal) denotes the mass of the satellite halo that is enclosed within r_tidal and Φ_host is the gravitational potential of the host, which in our case is given analytically as described in Sect. 4.1.

In Fig. E.2, we display the tidal radius as a function of time for various DM models and show its enclosed mass. It is visible that all DM models evolve fairly similarly, especially the SIDM models. This is related to the fact that the tidal radius is typically relatively large compared to the initial scale radius of the satellite halo (r_s = 1.28 kpc).

Appendix F Effective cross-section of the satellite

In this appendix, we study the effective cross-section given by Eq. (C.1) for our simulation of a satellite halo begin subject to an isotropic velocity-dependent cross-section. We compute the results given in Fig. F.1 using simulation Y as specified in Table A.1. The effective cross-section requires a characteristic velocity dispersion. Here, we use the average velocity dispersion within two different radial ranges, namely r < 0.05 kpc and 0.7 kpc < r < 1.5 kpc. In addition, we show the case of the constant cross-section of 80 cm² g⁻¹ as we use it for the simulations for a velocity-independent cross-section.

It is visible in Fig. F.1, that the effective cross-section for the velocity-dependent case is larger in the centre of the halo (red, solid) than for the constant case (blue). This leads to an enhanced core formation. However, later on, the effective cross-section drops below 80 cm² g⁻¹. At about the first pericentre passage (≈ 1.5 Gyr) tidal stripping is effective and reduces the velocity dispersion of the satellite halo, as discussed in Sect. 4.3. As a consequence, the effective cross-section sharply rises for the velocity-dependent case. This is the case at small (dotted) and large (dashed) radii. Subsequently, the effective cross-section at small radii starts to decline because of the increase in velocity dispersion due to the gravothermal collapse. At larger radii, the impact of the second pericentre passage is visible, increasing the effective cross-section for later stages.

Fig. F.1 Effective cross-section as a function of time. We compute the effective cross-section for different radial ranges from the average velocity dispersion using the satellite simulation for an isotropic velocity-dependent cross-section (simulation Y of Table A.1). Moreover, we show the velocity-independent case, which does not have any radial dependence.

Appendix G Highest-resolved simulation

In this Appendix we provide additional information for our highest-resolved simulation evolved in isolation (simulation Q in Table A.1). Figure G.1 gives the density (top panel), the one-dimensional velocity dispersion (middle panel), and the Knudsen number (bottom panel) as a function of radius for several times. Here, the gravothermal evolution driven by a velocity-independent cross-section including the core expansion and collapse phase is well visible.

Fig. G.1 Density, velocity dispersion, and Knudsen number as a function of radius. We show the results for our most highly resolved halo (simulation Q of Table A.1) at various stages of the gravothermal evolution.

All Tables

All Figures

Fig. 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 A.1; the corresponding names are A, B, C, D, E, and F.

In the text

Fig. 2 Variation of the maximum sampling radius for the ICs. Analogously to Fig. 1, 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 rs (orange and light blue) and 15 rs (purple and dark blue) are shown. We note that the lower resolution simulations (orange and purple) employ a softening length of ϵ = 6.4 pc, whereas the higher resolution simulations (light and dark blue) employ a softening length of ϵ = 1.0 pc. All parameters for the shown simulations (G, J, M, and O) are given in Table A.1.

In the text

	Fig. 3 Variation of resolution and gravitational softening length. Similar to Fig. 1, 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 A.1 are shown.
In the text

Fig. 4 Impact of energy conservation on the evolution time of the halo. Following Fig. 1, 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. 3.1). All parameters for the displayed simulations J and L are given in Table A.1.

In the text

Fig. 5 Evolution of an isolated halo as a function of time for two different choices of Nngb. The upper panel gives the mass enclosed within 10 pc and the lower panel displays the energy conservation as a function of time. The choice of Nngb = 384 corresponds to twice the kernel size compared to the SIDM computations for Nngb = 48. Table A.1 gives the parameters employed for the displayed simulations J and K.

In the text

	Fig. 6 SIDM kernel size divided by the mean free path as a function of radius for several times. The results, computed according to Eq. (2), are shown for several simulations with different resolutions. All simulations were run with Nngb = 48. More details for the shown simulations J, O, and Q are given in Table A.1.
In the text

Fig. 7 Employing a minimum time step. As in Fig. 1, we show the enclosed mass within 10 pc and the energy conservation for various simulations. This time, we investigate how much the simulations are affected when employing a minimum time step. The brighter lines indicate simulations without a minimum time step, i.e. the time step can decrease further. In contrast the darker lines feature a minimum time step of Δtmin = 2.8 × 10−5 Gyr. The displayed simulations J, N, O, Jt, Nt, and Ot are described in Table A.1.

In the text

Fig. 8 Density, velocity dispersion, and Knudsen number for simulations employing a minimum time step. The purple lines give the results for the simulation with N = 2 × 106 and ϵ = 6.4 pc (grey) and Δtmin = 2.8 × 10−5 Gyr and the blue lines are for the simulation with N = 1 × 107 and ϵ = 1.0 pc and Δtmin = 2.8 × 10−5 Gyr. These are the same simulations as in Fig. 7. The top panel gives the density as a function of radius for different times. We note that these times are during the collapse phase, and for most of them the energy conservation error is sizeable. The middle panel gives the velocity dispersion as a function of radius, and the bottom panel displays the Knudsen number (Eq. (3)). All parameters for the shown simulations Jt and Ot are given in Table A.1.

In the text

	Fig. 9 Time evolution of the enclosed mass for our most highly resolved simulation (simulation Q in Table A.1). The relative change of the enclosed mass within various radii is shown as a function of time. The enclosed mass is computed in projection, i.e. it is the mass within a cylinder of radius r2D.
In the text

Fig. 10 Projected logarithmic density profile slope as a function of projected enclosed mass. For our most highly resolved simulation evolved in isolation (simulation Q from Table A.1), we show how the projected logarithmic density profile slope, γ2D, and the projected enclose mass, M<r2D, evolve with time. We display γ2D and M<r2D for several radii r2D (see inset for legend). The larger circles mark the value at the beginning of the simulation, i.e. for an NFW profile. Moreover, the small white dots are placed equidistant in time, every 0.98 Gyr.

In the text

Fig. 11 Enclosed mass within several radii as a function of time. The results are shown for simulations with the standard time step criterion (darker curves) and those employing a minimum time step (lighter curves). The three panels give the enclosed mass within different radii for the same simulations, as shown in Fig. 7. Moreover, we show the enclosed mass for the initial NFW halo as a reference (grey).

In the text

Fig. 12 Velocity and angular dependent cross-sections. Following Fig. 1, the enclosed mass within 10 pc (upper panel) and the energy conservation (lower panel) are shown. In detail, a velocity-independent isotropic cross-section (blue), a velocity-independent forward-dominated cross-section (light blue), a velocity-dependent isotropic cross-section (red), and a velocity-dependent forward-dominated cross-section (light red) are displayed. They all have the same resolution of N = 107 and a gravitational softening length of ϵ = 1 pc. More details on the simulations O, P, R, and S can be found in Table A.1.

In the text

	Fig. 13 Velocity anisotropy as a function of radius. We show β (Eq. (5)) for our most highly resolved simulation (simulation Q of Table A.1) at different times.
In the text

	Fig. 14 Evolution of an isolated halo (solid line) and a satellite halo (dashed line). The simulations are for a velocity-independent cross-section, employ N = 107 particles, and use a gravitational softening length of ϵ = 1 pc. All parameters for the displayed simulations O and W are given in Table A.1.
In the text

Fig. 15 Mass enclosed within 10 pc as a function of time for the satellite halos. The upper panel gives the results for a velocity-independent cross-section (simulations U, V, W, and X), and the lower panel is for velocity-dependent cross-sections (simulations Y and Z). More details on the simulations can be found in Table A.1. In addition, the grey dotted line indicates the preferred mass range for the perturber of the stellar stream GD-1 (Bonaca et al. 2019).

In the text

	Fig. 16 Density profile fitted by a King model for a satellite halo. The results are for the simulation with a velocity-dependent isotropic cross-section. All simulation parameters are given in Table A.1 under simulation Y.
In the text

Fig. B.1 Variation of the gravitational time step parameter η. In the upper panel, we show the mass enclosed within the inner 10 pc as a function of time. The lower panel gives the evolution of the total energy relative to the absolute initial value. The shown simulations are for the low resolution of N = 2 × 106 and share the same parameters except η, which is varied as indicated in the legend. All parameters are given in Table A.1 under the names H, I, and J. The shaded regions in the top panel indicate the uncertainty estimated based on shot noise.

In the text

Fig. C.1 Enclosed mass and energy conservation as a function of time. We show the same simulations as in Fig. 12 (O, P, R, and S in Table A.1). The upper panel gives the mass within a radius of 0.15 kpc and the lower panels shows how well the total energy is conserved. The time for the simulations with a velocity-independent cross-section is rescaled to match the collapse phase of the simulations with velocity-dependent interactions.

In the text

	Fig. C.2 Effective cross-section as a function of radius. For the isolated halo evolved with an isotropic velocity-independent cross-section (simulation R according to Table A.1), we show the effective cross-section for several times as indicated in the legend. In addition, we display the cross-section for the velocity-independent cases (black line).
In the text

	Fig. D.1 Density as a function of radius. Simulation results and a fit using the King model (Eq. (9)) are shown for different times of the late collapse phase. The results are for the isolated simulation with a resolution of N = 5 × 107 particles (simulation Q of Table A.1).
In the text

Fig. D.2 Parameters of fitted King model and further quantities as a function of time. The parameters that we determine by fitting a King model (Eq. (9)) to our highest-resolved simulation (Q in Table A.1) are shown in the top panel as a function of time. The bottom panel gives the mass enclosed in rc according to the fitted model as well as the maximum circular velocity directly computed from the simulation data.

In the text

	Fig. E.1 Trajectory of the satellite. The three-dimensional position of the subhalo for a CDM and an SIDM simulation is shown. The corresponding simulations are T and W according to Table A.1, where their parameters are given.
In the text

Fig. E.2 Tidal radius and its enclosed mass as a function of time. For the satellite halo, evolved with different DM models, the tidal radius (solid lines) according to Eq. (E.1) and the mass within it (dashed lines) are shown. All simulations shown here share a resolution of N = 107 particles and a gravitational softening length of ϵ = 1.0 kpc. Further information on the simulations T, W, X, Y, and Z is provided in Table A.1.

In the text

	Fig. F.1 Effective cross-section as a function of time. We compute the effective cross-section for different radial ranges from the average velocity dispersion using the satellite simulation for an isotropic velocity-dependent cross-section (simulation Y of Table A.1). Moreover, we show the velocity-independent case, which does not have any radial dependence.
In the text

	Fig. G.1 Density, velocity dispersion, and Knudsen number as a function of radius. We show the results for our most highly resolved halo (simulation Q of Table A.1) at various stages of the gravothermal evolution.
In the text