Cosmological simulations with rare and frequent dark matter self-interactions

Dark matter (DM) with self-interactions is a promising solution for the small-scale problems of the standard cosmological model. Here we perform the first cosmological simulation of frequent DM self-interactions, corresponding to small-angle DM scatterings. The focus of our analysis lies in finding and understanding differences to the traditionally assumed rare DM (large-angle) self scatterings. For this purpose, we compute the distribution of DM densities, the matter power spectrum, the two-point correlation function and the halo and subhalo mass functions. Furthermore, we investigate the density profiles of the DM haloes and their shapes. We find that overall large-angle and small-angle scatterings behave fairly similarly with a few exceptions. In particular, the number of satellites is considerably suppressed for frequent compared to rare self-interactions with the same cross-section. Overall we observe that while differences between the two cases may be difficult to establish using a single measure, the degeneracy may be broken through a combination of multiple ones. For instance, the combination of satellite counts with halo density or shape profiles could allow discriminating between rare and frequent self-interactions. As a by-product of our analysis, we provide - for the first time - upper limits on the cross-section for frequent self-interactions.


INTRODUCTION
Although many efforts have been made to uncover the nature of dark matter (DM), it remains largely unknown even after several decades of research. To narrow down the large number of models that contain potential DM candidates, a huge variety of experiments, based on direct and indirect detection, are being carried out. Moreover, forthcoming astronomical surveys with upcoming telescopes such as E 1 (Euclid Collaboration et al. 2020), Rubin Observatory 2 (Zhan & Tyson 2018) and Roman 3 (Spergel et al. 2015) promise to tighten constraints on cosmological models and to discriminate between models of DM beyond the cold collisionless DM of the standard cosmological model (Lambda cold dark matter, ΛCDM).
There exist a number of SIDM models with a set of free parameters, such as the total cross-section, the angular and velocity dependence, and the nature of the scattering (elastic or inelastic). Most studies have assumed models that are isotropic, elastic, and velocity-independent. Burkert (2000) performed the first simulations with a Monte-Carlo scheme where the numerical particles were treated analogously to physical DM particles. Since then many variants of SIDM have been studied, such as inelastic scattering (e.g. Essig et al. 2019;Huo et al. 2020;Shen et al. 2021) including multistate scattering (Schutz & Slatyer 2015;Vogelsberger et al. 2019;Chua et al. 2020) or even multicomponent DM (Todoroki & Medvedev 2018;Vogelsberger et al. 2019). Also anisotropic cross-sections have been investigated (Robertson et al. 2017b;Banerjee et al. 2020;Nadler et al. 2020).
However, if the self-interaction cross-section is strongly anisotropic, particles scatter by tiny angles. This implies a much lower momentum and energy transfer per scattering event compared to an isotropic cross-section. Hence, small-angle scattering must be more frequent to have a similar effect on the DM distribution. Models having different angular dependences might be compared by using the momentum-transfer cross-section. By frequent self-interacting DM (fSIDM), we refer to a limit where the scattering angles become infinitesimal small, while the momentum-transfer cross-section stays constant. In contrast, we refer to rare self-interacting DM (rSIDM) for less anisotropic differential cross-sections.
Frequent self-interactions have gained popularity in the context of galaxy cluster mergers because they can explain larger DM-galaxy offsets than rSIDM (Kahlhoefer et al. 2014;Fischer et al. 2021a,b). However, numerical schemes that treat numerical particles like physical ones are not capable of simulating fSIDM. Only recently, a general solution to this problem has been found. Fischer et al. (2021a) developed a new scheme that allows one to model frequent scattering within -body simulations.
However, previous fSIDM studies (Kahlhoefer et al. 2014(Kahlhoefer et al. , 2015Kummer et al. 2018Kummer et al. , 2019Fischer et al. 2021a,b) only considered idealized cases that never took the full cosmological context into account. In contrast, for rSIDM, there are a number of recently published simulations of cosmological boxes Rocha et al. 2013;Vogelsberger et al. 2016;Robertson et al. 2019Robertson et al. , 2020Banerjee et al. 2020;Stafford et al. 2020Stafford et al. , 2021Harvey et al. 2021;Ebisu et al. 2022) or using zoom-in simulations (Vogelsberger et al. 2012(Vogelsberger et al. , 2016Zavala et al. 2013Zavala et al. , 2019Vogelsberger et al. 2014;Fry et al. 2015;Robertson et al. 2018;Despali et al. 2019;Robles et al. 2019;Nadler et al. 2020;Vega-Ferrero et al. 2020;Bondarenko et al. 2021;Sameie et al. 2021;Shen et al. 2021Shen et al. , 2022Bhattacharyya et al. 2022;Silverman et al. 2022;Sirks et al. 2022). These simulations have been used to study the phenomenology of SIDM models on various mass scales, such as dwarf galaxies, Milky Way (MW)-like galaxies, and galaxy clusters. Several properties of the DM haloes such as their density profile and their shape have been measured and predictions for observations, such as gravitational lensing, have been made. This enabled constraints to be put on the total cross-section of rSIDM models, while fSIDM models have remained poorly constrained. In this study, we investigate fSIDM, for the first time using a cosmological simulation.
This paper aims to study the effects of SIDM on large scales to understand the differences between fSIDM and isotropic rSIDM. In this first effort, we assume the self-interactions to be velocityindependent and elastic. We conduct DM-only cosmological -body simulations using a full box as well as zoom-in simulations and study various properties of the DM distribution.
In Section 2 we briefly describe our numerical methods and present the setup for our cosmological simulations. The results are presented in Section 3. For example, we show the matter power spectrum, the halo mass function as well as density and shape profiles of DM haloes. A discussion of our results, their limitations, implications and further perspectives follows in Section 4. Finally, we summarize and conclude in Section 5. Additional details and plots are provided in the appendices.

NUMERICAL SETUP
In this section, we describe the numerical setup for this study. We give details of the code and algorithms that we have used as well as our simulations and their initial conditions.
In this paper, we used the cosmological -body code -3, and the predecessor -2 is described in Springel (2005). The implementation of rare and frequent self-interactions has previously been described in Fischer et al. (2021a,b). Additionally, we implemented the comoving integration for the SIDM module to perform cosmological simulations. A test problem that demonstrates that the comoving integration is working as expected can be found in Appendix A. To match rare and frequent self-interactions, we use the momentum-transfer cross-section, We simulated a full cosmological box and and also performed zoom-in simulations. All simulations are DM only and the selfinteractions are always velocity-independent and elastic. In the case of rSIDM, the differential cross-section is isotropic, while fSIDM corresponds to a very anisotropic cross-section.
To generate the initial conditions for the full box, we use N-G IC (Springel 2015). The initial conditions are similar to box4 of the Magneticum simulations 5 with a comoving side length of 48 Mpc ℎ −1 . We run simulations with different resolutions and refer to those using the naming convention of Magneticum (hr and uhr). Our highest resolution run (uhr) contains ∼ 1.9 × 10 8 simulation particles. More details on the full cosmological box are given in Tab. 1. Moreover, we performed cosmological zoom-in simulations with different resolutions of the same region. The region is selected from a large box with a comoving side length of 1 Gpc ℎ −1 . Several publications (e.g. Planelles et al. 2013;Rasia et al. 2015) have used this box for zoom-in initial conditions and it was first described in Bonafede et al. (2011). In our zoom-in region, the most massive halo has a virial mass of ∼ 8.8 × 10 11 M ⊙ ℎ −1 . Further details can be found in Tab. 2. In addition, we provide in Appendix B a convergence test of the density profile of the most massive halo.
For the analysis, we identify DM haloes using the friends-of- 4 We implicitly assume identical particles, in this case, the definition is equivalent to the one recommended by Robertson et al. (2017b) and Kahlhoefer et al. (2017). , and the mass of the numerical DM particles ( DM ) as well as the momentum-transfer cross-section per physical DM particle mass (T/ ). The non-zero cross-sections have been simulated using fSIDM and isotropic rSIDM. All simulations share the same initial conditions but with a different resolution. name high res 1x ∼ 4.51 × 10  Table 2. Properties of the zoom-in simulations. We provide the name of the simulation, the number of particles in the highly resolved region ( high res ), the mass of the high-resolution particles ( DM ) and the cross-sections we simulated (T/ ). The non-zero cross-section has been simulated using fSIDM and isotropic rSIDM. All simulations share the same initial conditions but with a different resolution. friends algorithm 6 implemented along with -3. The builtin module also identifies substructure within the haloes (Springel et al. 2001;Dolag et al. 2009). We use halo and subhalo positions, masses and radii as provided by . The virial radius, vir , and the virial mass, vir , are measured with the sphericaloverdensity approach based on the overdensity predicted by the generalized spherical top-hat collapse model (e.g. Eke et al. 1996). Here, vir is defined as the radius at which the mean density becomes larger than the one of the top-hat collapse model and vir is the mass inside vir . Every halo contains at least one subhalo, which is the primary subhalo located at the same position as the halo (determined by the location of the most gravitationally bound particle). The primary subhalo typically contains most of the particles that belong to the halo.

RESULTS
In this section, we present the results of our simulations and compare the effects of DM models. To this end, we study several statistical properties such as the matter power spectrum, the probability density function (PDF) of the DM densities, the two-point correlation function, and the halo and subhalo mass function. We then study the impact of self-interactions on the density and circular velocity profile of DM haloes. Furthermore, we investigate how the shapes of haloes change when self-interactions are present. Besides, we study qualitative differences between rSIDM and fSIDM and discuss transferring constraints on the cross-section of rare scatterings to frequent self-interactions. 6 A description of the friends-of-friends algorithm can, for example, be found in More et al. (2011).

Surface density
In Fig. 1, we show the surface density of the full cosmological box for our CDM and fSIDM (T/ = 1.0 cm 2 g −1 ) simulations. At a cosmological redshift of = 0, basically, no differences between the simulations are visible. Hence, fSIDM seems to agree well with the collisionless DM on large scales. Previous studies that examined the large-scale structure in rSIDM found that it looks like CDM, but differences arise on small scales (e.g. Rocha et al. 2013;Stafford et al. 2020). Hence, SIDM keeps the success of CDM in explaining the large-scale structure but could be capable of resolving small-scale issues. In the next sections, we investigate quantitatively the effects of fSIDM and rSIDM.

Matter power spectrum
Cosmological structure is characterized by the matter power spectrum. Stafford et al. (2021) computed it for several cosmologies including isotropic SIDM. They found a suppression of small-scale structure with increasing cross-section, while on large scales SIDM behaves like CDM.
In Fig. 2, we show the matter power spectrum of our cosmological full-box simulations and compare the various DM models. To compute the power spectrum, we use the same code as in Grossi et al. (2008). In line with Stafford et al. (2021), we find that the self-interactions affect the small scales (high -values) only and can lead here to substantial suppression of structures. The stronger the self-interactions, the stronger the suppression of structure formation on small scales. However, we only find small differences between the rSIDM and fSIDM simulations. In particular, for the larger crosssection ofT/ = 1.0 cm 2 g −1 , the effects are very similar.

Distribution of DM densities
Here, we compute the volume-weighted probability to find a given density, i.e. the PDF for a random position within the cosmological volume. First, we use an oct-tree to find neighbours of the simulation particles. For each particle, we find the radius that contains 160 neighbouring particles (for the highest resolution run, uhr). Using that radius we estimate a physical density and divide it by the particle mass to obtain a volume that is associated with the particle. From all particles, we sum their associated volumes per logarithmic density bin and divide by the simulation volume and logarithmic density bin size. The result is shown in Fig. 3. We find that the various DM models differ in the high-density regime only. Here, self-interactions suppress the highest densities compared to CDM. The suppression takes place for densities ≳ 10 7 M ⊙ kpc The densest regions are most sensitive to DM self-interactions because the effect of self-interactions depends on the density and the velocity dispersion, which tends to be high in dense regions. However, we do not find large differences between rSIDM and fSIDM, especially for the smaller cross-section ofT/

Two-point correlation function
In addition to the power spectrum, we can use the two-point correlation function to characterize the distribution of matter. Specifically, we compute the spatial two-point correlation function according to 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 [M kpc 3 h 2 ] 10 9 10 7 10 5 10 3 10 1 dp/dln CDM, 0.0 cm 2 g 1 fSIDM, 0.1 cm 2 g 1 fSIDM, 1.0 cm 2 g 1 rSIDM, 0.1 cm 2 g 1 rSIDM, 1.0 cm 2 g 1 Figure 3. The probability to find a given density per logarithmic density bin is shown as a function of density for various simulations. The colours indicate the type of self-interaction and the line style gives the strength of self-interaction as indicated in the legend. The plot is for a redshift of = 0 and produced from the high-resolution full cosmological box simulations (uhr). (Davis & Peebles 1983) We use the data points given by the simulation and draw random numbers to generate points from a uniform PDF. The number of data points is given by , and denotes the number of randomly distributed points. We compute the number of distances ( ) between data points within the interval [ , + Δ ] as well as ( ), the number of distances between data points and random points. We compute the halo-halo correlation, hh ( ), based on the halo positions and the subhalo-subhalo correlation, ss ( ), based on the subhalo positions as identified by for the full cosmological box simulations. But we do not use all substructures; instead, we introduce a mass/resolution cut to avoid our results being affected by numerical artefacts. All systems with a mass of < 9.6 × 10 10 M ⊙ ℎ −1 are excluded as they are poorly resolved. For the uhr runs, this mass corresponds to 2200 particles and for the hr runs to 220 particles. As a consequence, we have much fewer positions (very roughly 3×10 3 for hr and uhr as well as haloes and subhaloes) for the computation of the two-point correlation function. In Fig. 4, we show the results for the halo positions (left-hand panel) and the subhaloes (right-hand panel). For the haloes, we do not find much of a difference among the DM models. In contrast, the comparison of the subhalo positions reveals a difference on small scales as well as larger scales. As expected, the stronger the cross-section, the more structures on small scales are suppressed. ForT/ = 1.0 cm 2 g −1 , fSIDM simulations deviate more from CDM than the rSIDM simulations do, with the subhalo correlation function being roughly 20% lower on small scales for fSIDM than CDM. Note that these results depend on the chosen mass/resolution cut for subhaloes. The lower the cut, the larger the deviations between the CDM and SIDM runs are.
The suppression of the subhaloes compared to haloes may arise from the fact that some of them are satellites and not exclusively primary subhaloes (see Sec. 2). The SIDM satellites could dissolve faster because they have lower central densities due to DM selfinteractions. Hence, they are not as strongly bound as their CDM counterparts, which makes them more prone to tidal effects. This could be enhanced by scattering between host and satellite particles. In the next section, we turn to the halo and subhalo mass function.

Halo and subhalo mass function
Here, we study the halo and subhalo mass function as well as the abundance of satellites.
In Fig. 5, we show the halo mass function (left-hand panel) and the subhalo mass function (right-hand panel). For the computation, we used the total mass of the haloes and subhaloes calculated by . At the high-mass end, we do not find a significant difference between the CDM and SIDM models. But at lower masses, selfinteractions suppress the number of haloes and even more the number of subhaloes. However, the smallest objects we display here are not well resolved. Note that the mass cut previously used for Fig. 4 selects systems with a mass ≳ 9.6×10 10 M ⊙ ℎ −1 . At least the largest fSIDM cross-section studied here gives us a significant reduction of the number of subhaloes above this cut. But for the haloes, none of the cross-sections results in a significant reduction for masses ≳ 9.6 × 10 10 M ⊙ ℎ −1 . The difference between the halo and subhalo mass function can only arise from satellites dissolving faster, as all non-satellites are haloes. Hence, satellites, in particular low-mass satellites, appear to be an interesting test bed for SIDM models. Several SIDM studies have focused on them (e.g. Kahlhoefer et al. 2019;Kaplinghat et al. 2019;Banerjee et al. 2020;Nadler et al. 2020Nadler et al. , 2021Nishikawa et al. 2020;Sameie et al. 2020;Yang et al. 2020;Correa 2021;Zeng et al. 2022).
In Fig. 6, we study the abundance of satellites for the DM models. For each group, we compute the number of satellites per logarithmic mass as a function of the satellite mass divided by the virial mass of the group. We select all subhaloes, except the primary one, within 5 vir . The results for a selection criterion of 1 vir are shown in Appendix C. For the full-box simulations, we do this for the 100 most massive groups and take the mean (left-hand panel). In contrast, for We plot the number density of haloes/subhaloes per logarithmic mass bin as a function of the total halo/subhalo mass as identified by . The colours indicate the type of the self-interaction and the line style gives the strength of self-interaction as indicated in the legend. We display results of the higher resolution boxes (uhr, darker lines) and the lower resolution boxes (hr, fainter lines). The dash-dotted grey line indicates the mass limit of ∼ 9.6 × 10 10 M ⊙ ℎ −1 that we applied previously for the two-point correlation function. The plots are for a redshift of = 0. the zoom-in simulations, we only take the mean of the three most massive systems (right-hand panel). The comparison between SIDM and CDM in the lower panel shows that self-interactions can suppress the abundance of satellites. The stronger the self-interactions, the fewer satellites are present, and low-mass satellites are more affected than more massive ones. However, we have to note that the lowest mass satellites used in Fig. 6 are not well resolved. Although the halo mass function is not converged for the hr run at low masses (left-hand panel), the differences to CDM, in particular for fSIDM withT/ = 1.0 cm 2 g −1 , seem to be converged. This makes it plausible that the effect for the lower masses is real. In a similar manner, Stafford et al. (2020) found the difference between cosmologies to be converged for density profiles below the convergence radius proposed by Ludlow et al. (2019). We find differences for better resolved haloes at larger masses too, but not for the most massive satellites. For the larger cross-section ofT/ = 1.0 cm 2 g −1 , frequent self-interactions seem to be much more efficient in reducing the number of satellites than rare scatterings. This is in line with the suppression of the spatial two-point correlation function of subhaloes at small scales (right-hand panel of Fig. 4). It is also worth mentioning, which we found earlier, in a study of head-on collisions of unequalmass mergers (Fischer et al. 2021b), that fSIDM subhaloes dissolve faster than rSIDM subhaloes when they are matched in terms of T / . We have to note that there cannot exist a perfect matching between rare and frequent self-interactions as they are qualitatively different. TheT-matching can only provide an approximation that can be helpful in some situations, e.g. density profile (see Sec. 3.6) or shape profiles (see Sec. 3.8). In Sec. 3.9, we further investigate the qualitative difference concerning the abundance of satellites.
In addition, we study the radial dependence of the suppression of the satellite abundance in Fig. 7. Here, we find that the effect of SIDM on the number of satellites becomes stronger at smaller distances to the host. As in Figs 5 and 6, we find that fSIDM reduces the abundance of satellites stronger than rSIDM given the same value for˜/ . However, we have to note that for small radii the number of satellites we study here is low and thus the error on the ratios shown in Fig. 7 is sizeable. Thus, we can only make a qualitative statement for the smallest radii, but not quantify the difference between the DM models.
That satellites can dissolve faster in the presence of SIDM has been found earlier. The evolution of a satellite and its lifetime depends on several aspects and physical mechanisms that are at play. In CDM, satellites are torn apart by tidal forces, which act against their gravitational self-binding. The more massive and the more concentrated a satellite is, the more resistant it is against the tidal forces and can survive longer. In the context of SIDM, a DM core can form, which flattens the gravitational potential and makes the satellite more prone to tidal disruption (Yang et al. 2020). However, if the satellite is in a further evolution phase and undergoes core collapse (e.g. Balberg et al. 2002;Koda & Shapiro 2011;Essig et al. 2019), it is more protected against tidal disruption as it is even more concentrated than its CDM counterparts. Nishikawa et al. (2020) found that tidal stripping enhances the core collapse in SIDM haloes, i.e. reduces the collapse time. A recent study (Zeng et al. 2022) found that it is nearly impossible for satellites to undergo core collapse if self-interactions are velocity-independent. The evolution of a satellite is not only determined by tidal stripping and DM self-interactions between satellite particles, it is also affected by tidal heating and DM self-interactions between satellite and host particles. Whether tidal effects (stripping and heating) enhance or prevent core collapse depends on the mass concentration of the satellite. The satellite-host interactions transfer energy into the satellite and thus contribute to the core formation. For a cross-section, which is decreasing with velocity (e.g. Loeb & Weiner 2011;Tullin et al. 2012), this effect would become much less as the typical relative velocity between satellite and host particles is larger than that between satellite particles. In consequence, the survival time of satellites is expected to depend crucially on the velocity dependence of the self-interactions. Banerjee et al. (2020) studied various SIDM models and found that the suppression of the satellite abundance is weaker if the cross-section is velocity-dependent. However, they also studied an anisotropic cross-section and did not find any significant deviation to an isotropic cross-section, and this also includes satellite counts. It is worth pointing out that our fSIDM cross-section is more anisotropic. In practice, it might be infeasible to study it with an rSIDM scheme, as used by Banerjee et al. (2020), because the high scattering rate would require very small time-steps. In addition, we should mention the work by Vogelsberger et al. (2019). They studied an MW-like halo with multistate inelastic DM selfinteractions and found that this type of interaction suppresses the abundance of small structures even for small cross-sections considerably.
Overall, satellites are not only interesting objects that constrain the strength of SIDM, but may also provide bounds on the angular and velocity dependence of the differential cross-section. In Sec. 3.9, we discuss further potential strategies to constrain this angular dependence.

Density profiles
As we have seen in Fig. 3, in scenarios with SIDM the high-density regions are suppressed. In particular, density cores have been studied in the literature and used to constrain self-interactions (e.g. Correa 2021; Sagunski et al. 2021;Ray et al. 2022). Recent strong lensing observations provide further evidence for DM cores in galaxy clusters (Limousin, Marceau et al. 2022).
In cosmological simulations of rSIDM, density and circular velocity profiles of haloes have been studied previously by Rocha et al. (2013) In Fig. 8, we show median density profiles of the haloes for three different mass bins. In the outer regions, the density profiles are very similar among the various DM models. But in the central region, we observe a density core for SIDM while CDM predicts cuspy haloes. The central density is lower for a larger cross-section and frequent self-interactions lead mostly to slightly larger cores than rare self-interactions, if aT-matching (same momentum-transfer crosssection for rSIDM and fSIDM) is employed. This is in agreement with previous findings of isolated haloes (Fischer et al. 2021a).
The central density of the individual haloes is shown in Fig. 9. Specifically, we computed the mean density within 0.01 vir . As we have already seen above, a larger self-interaction cross-section leads to haloes with lower central densities. For CDM, we observe a slight decline of the central density with virial mass. But for the SIDM models, the decline is steeper. This implies an increasing difference between CDM and SIDM models with virial mass. Note that with increasing virial mass the velocity dispersion in the inner regions of the DM haloes increases and thus the self-interactions are more efficient.

Circular velocity
In addition to computing the density profile, we can also study the circular velocity, as a function of the radius. In Fig. 10 we show the circular velocity profile for the most massive subhalo of the best resolved zoom-in simulation. It is visible that in the inner regions the velocity that is needed for a circular orbit is less for SIDM compared to CDM. This is a direct consequence of the density core, i.e. the enclosed mass, (< ), in Eq. (3) is less. At larger radii, the circular velocity is the same across the DM models. In order to trace the core formation in terms of circ , we can measure it at a small radius and compare it to its maximum value. We do so in Fig. 11 and plot circ (2 kpc) vs. circ,max for the 300 most massive subhaloes of the best resolved zoom-in simulation. For the most massive subhaloes, i.e. the ones with a larger value for max , we find a difference between the DM models. The SIDM models tend to have lower circular velocities at 2 kpc than CDM. For fSIDM, circ is slightly lower than for rSIDM, which is in line with the lower densities found for fSIDM. At lower masses, we do not find a significant difference; here a measure at 2 kpc may not CDM, 0.0 cm 2 g 1 fSIDM, 0.1 cm 2 g 1 fSIDM, 1.0 cm 2 g 1 rSIDM, 0.1 cm 2 g 1 rSIDM, 1.0 cm 2 g 1 Figure 9. We display the central density of haloes as a function of their virial mass. The results of various simulations at a redshift of = 0 are shown. The central density is measured as the mean density within a sphere of 0.01 vir . Systems evolved with the smaller cross-section are marked by "+" and for the larger cross-section we use "×"; the CDM case is indicated by " ". In addition to the individual systems, we computed the mean of the distribution as a function of virial mass, indicated by the lines. The shaded region gives the standard deviation. be sensitive to the density core. Instead, a measure at smaller radii would be preferable to trace the density core. Self-interactions have turned out to be interesting to explain the diversity of observed rotation curves (e.g. Oman et al. 2015). In particular, the response of SIDM to the gravitational potential of the baryons can increase the diversity of rotation curves (e.g. Creasey et al. 2017;Kamada et al. 2017;Kaplinghat et al. 2019). However, Zentner et al. (2022) claim that it is not clear whether observed galactic rotation curves are better explained by SIDM or CDM including baryons. For late-type dwarfs Roper et al. (2022) points out that the measured rotation curves could differ significantly from their circular velocity curves. Other studies had also previously indicated that In any case, this question can only be addressed once baryonic physics is taken into account in the -body code, implying that our results cannot be directly compared to these findings. We leave a more in-depth discussion of these issues for future work.

Shapes
The shape of DM haloes can provide bounds on the strength of selfinteractions. We use the tensor, S, to compute the shape. For point masses, it is with the particle mass at position r . In order to compute the shapes, we use the eigenvalues ( 1 ≥ 2 ≥ 3 ). They are related to the semi-axes ( ≥ ≥ ). In particular, we measure the axial ratios = / = 3 / 1 and = / = 2 / 1 . Based on these two ratios, we compute the triaxiality (Franx et al. 1991), Our method to measure halo shapes is iterative, where during each iteration the computation of the tensor, S, uses all particles within an ellipsoid. The orientation and axial ratios of this ellipsoid are based on S from the previous iteration, and the volume of the ellipsoid is kept constant throughout the iterations. We iterate until the two axial ratios ( and ) are converged. Note that in the literature various methods to measure the shape of haloes have been used (e.g. Zemp et al. 2011;Peter et al. 2013;Sameie et al. 2018;Robertson et al. 2019;Banerjee et al. 2020;Chua et al. 2020;Harvey et al. 2021;Vargya et al. 2022;Shen et al. 2022). In Fig. 12, we show the median shapes for three halo mass bins. Note that for the computation we have used all particles, not only those that belong to the halo as identified by . The first mass bin contains only a few objects leading to more noise in the upper We compute the median shape as well as the median semimajor axis from ellipsoids having the same volume. The shaded regions indicated the scatter among the haloes, and the range between the 25th and 75th percentiles is displayed. For the sake of clarity, we show this only for the CDM simulation and the fSIDM run with the larger cross-section. The virial mass given in the panels indicates the median virial mass of the haloes of the corresponding mass bin. The same applies to the shown virial radius. Note that we are using the same mass bins as in Fig. 8. The shapes are computed from the highest resolved full cosmological box and for a redshift of = 0. panel. At small distances from the centre, the self-interactions lead to rounder haloes, mainly depending on the strength of the selfinteractions. This has been discovered earlier (e.g. Peter et al. 2013). The difference between fSIDM and rSIDM is small, with fSIDM leading to haloes that are a little more spherical than the rSIDM haloes, given aT-matching. At larger radii, when more particles are included, the shapes become very similar between the different DM s = c/a at r = 0.078 r vir CDM, 0.0 cm 2 g 1 fSIDM, 0.1 cm 2 g 1 fSIDM, 1.0 cm 2 g 1 rSIDM, 0.1 cm 2 g 1 rSIDM, 1.0 cm 2 g 1 Figure 14. The shape of individual haloes as a function of their virial mass at a redshift of = 0 is shown. The shapes are measured as = / within an ellipsoid that has the volume of a sphere of 0.078 vir . Systems evolved with the smaller cross-section are marked by "+" and for the larger one we use "×"; the CDM case is indicated by " ". In addition to individual systems, we computed the mean of the distribution as a function of virial mass, indicated by the lines. The shaded region gives the standard deviation.
models. These additional particles are located at lower densities and are hence hardly affected by self-interactions. However, the SIDM shapes differ from CDM at radii beyond the density core more than the density profiles do. We have evaluated this in great detail in Appendix D. Furthermore, we find that the radius where they start to differ from each other depends strongly on the strength of the selfinteractions. This has been previously pointed out by Vargya et al. (2022) and suggested as a measure of the total cross-section. In Fig. 13, we also show a complementary shape variable, the triaxiality , as given in Eq. (5) for the same mass bins as in Fig. 12. We observe that is decreasing with increasing cross-section, implying that the haloes become less prolate. For large radii, where the matter density becomes smaller, differences between the DM models vanish as expected.
We finally show the shape ( = / ) measured in ellipsoids with a volume equal to a sphere with 0.078 vir of individual systems as a function of the virial mass in Fig. 14. As already expected from the shape profiles, we observe that the haloes become more spherically symmetric with increasing cross-section. Again, fSIDM leads to rounder haloes than rSIDM. There is no significant qualitative difference between the trend of rare and frequent scattering with virial mass. At the high-mass end, we have only very few objects such that the apparent decrease of may not be significant.

fSIDM versus rSIDM
The main aim of this paper is to understand whether and how the phenomenology of rare and frequent self-interactions can help to distinguish between them. In the following, we first investigate a potential degeneracy between rare and frequent self-interactions in terms of halo shape and central densities. In Fig. 15, we plot the shape as a function of central density for individual systems. The arrows indicate how they change when increasing the cross-section. Here, we use the three most massive systems of the full-box simulation (uhr). Although the values for the same object and cross-section differ from ) of the highest resolved fullbox simulation at a redshift of = 0. How they evolve when increasing the cross-section is indicated by the arrows. each other, this is not necessarily a qualitative difference. Rather it can be interpreted as an issue of matching rare and frequent selfinteractions. As the two types of self-interaction roughly move in the same direction (in the shape-central density plane) when increasing the cross-section.
Secondly, we consider the number of satellites to find a qualitative distinction. In Sec. 3.5, we found that it can make a difference for the abundance of satellites whether DM self-interactions are frequent or rare. Compared to the differences, we found that for the central densities (Sec. 3.6) or the shapes (Sec. 3.8) the number of satellites seems to be more sensitive to the underlying form of the cross-section. In Fig. 16, we plot the central density (upper panel) and the shape (lower panel) as a function of the number of satellites for the same systems as in Fig. 15. We only count satellites that are within the virial radius and more massive than 0.0008 vir , which corresponds to a resolution of ≳ 2200 particles. If fSIDM and rSIDM behave qualitatively the same, the systems should move along the same path in the plot when increasing the cross-section, regardless of whether the scattering is rare or frequent. We use arrows to indicate how the systems change with increasing cross-section. It becomes clear that fSIDM and rSIDM are not scaled versions of each other but are qualitatively different. For a given number of satellites, we find the main halo to be less dense and more spherically symmetric for rSIDM compared to fSIDM. If this qualitative difference is still present in full-physics simulations and strong enough to be observable, it could provide an avenue to distinguish between rSIDM and fSIDM.
A possible explanation of this difference between rSIDM and fSIDM goes as follows. In addition to tidal stripping (which is present even with collisionless DM), self-interactions between host and satellite particles can reduce the mass of, or even destroy, subhaloes as they orbit their host. This is because energy is transferred from the hot (high velocity dispersion) host halo into the cold (lower velocity dispersion) subhalo by self-interactions. In the case of rSIDM, this energy transfer takes the form of a small fraction of subhalo particles receiving large velocity kicks, which unbind them from the subhalo. . We show the three most massive haloes of our full-box simulation at = 0, the same as in Fig. 15. The upper panel gives the central density within a radius of 0.01 vir as a function of the number of satellites and the lower panel gives the shape of an ellipsoid that has the same volume as a sphere with a radius of 0.078 vir as a function of the number of satellites. For the satellite count, only satellites within vir and with a mass larger than 0.0008 vir were considered. This implies that every satellite is at least resolved by ≳ 2200 particles. The arrows connect the same systems across the simulations with various cross-sections.
In contrast, with fSIDM, all particles in the subhalo will interact with host halo particles, but will not receive kicks as large as in the rSIDM case. With the rSIDM and fSIDM cross-sections matched such that the same total energy transfer from host halo to subhalo takes place, and if we consider a time where the total transferred energy is just enough to unbind the subhalo, the fSIDM halo will spread this energy amongst all its particles and become unbound, while the rSIDM halo will lose some of this energy to subhalo particles that leave the subhalo at high velocity. In consequence, fSIDM satellites dissolve faster than their rSIDM equivalents.

Constraints on fSIDM
In this paper, we do not aim to compare our simulations to observations but previous studies of rSIDM may allow placing constraints on the velocity-independent fSIDM cross-section. So far we have found that fSIDM and rSIDM behave similarly in many aspects. In particular, if we consider the density (Fig. 8) or shape (Fig. 12) profiles, we find that the momentum-transfer cross-section can roughly match fSIDM and rSIDM. This can allow transferring constraints of the cross-section for rSIDM to fSIDM. In many situations, we found that fSIDM has a slightly larger effect than rSIDM, when the same momentum-transfer cross-section is considered. This implies that upper limits on the momentum-transfer cross-section may be even more stringent for the case of fSIDM. Given this, the established upper limits on rSIDM could be viewed as conservative limits for the case of fSIDM. In particular, this applies to our shape measurements for the whole mass range we have studied and to the density profiles for the two most massive mass bins (≳ 10 13 M ⊙ , see also Fig. 9). In this context, the work of Sagunski et al. (2021) is relevant. They studied the core sizes of galaxy groups and clusters to derive limits on the self-interaction cross-section. Their upper limits at a confidence level of 95% on the total cross-section are / for clusters. As discussed above, these limits should hold for fSIDM too.
Another relevant study is Peter et al. (2013). They examined the shapes of galaxy clusters and consider a velocity-independent crosssection of / = 1.0 cm and should apply to fSIDM too. However, we have to note that Peter et al. (2013) used DM-only simulations. Robertson et al. (2019) showed that the shape of the intracluster medium (ICM) is hardly affected by the DM physics, but that the DM component becomes considerably rounder due to self-interactions, even if baryons are present. Furthermore, the presence of baryons makes the DM shapes become more round too, which rules in favour of excluding a cross-section of˜/ ≥ 0.5 cm 2 g −1 for rSIDM and fSIDM.

DISCUSSION
In this section, we discuss the limitations and implications of our results. Our findings depend partially on , as we use the halo and subhalo positions as well as their mass and spatial extension. Only a few of our results are independent of , and these are the matter power spectrum and the density PDF. We want to point out that there are a number of codes for identifying substructure (e.g. Knollmann & Knebe 2009;Maciejewski et al. 2009;Tweed, D. et al. 2009;Behroozi et al. 2012;Han et al. 2017;Elahi et al. 2019). They employ different algorithms, which may give somewhat different results .
We found that the abundance of satellites is sensitive to differences between rare and frequent DM scatterings. In particular, for low-mass satellites, the difference between the DM models becomes larger. As such it would be interesting to study those low-mass satellites with a higher resolution. In this context, a hybrid approach as presented in Zeng et al. (2022) is an efficient technique that could help to improve our results.
In general, we would expect to find the largest differences between rSIDM and fSIDM in systems that are far from equilibrium. However, there are systems that may be similar or even more sensitive to the shape of the differential cross-section. Mergers have received attention in the context of SIDM (Randall et al. 2008;Kahlhoefer et al. 2014;Harvey et al. 2015;Kim et al. 2017;Robertson et al. 2017a,b;Fischer et al. 2021a,b). It has been shown that the DM-galaxy offsets can produce discernible differences between DM models. In addition to offsets, the distribution of a collisionless component (such as galaxies or stars) can give a handle to constrain DM properties.
The main simplification of our simulations is the neglect of baryonic matter and the associated physics. While DM-only simulations may provide reasonably accurate results on large scales, they fail on scales where galaxies form. In particular, the inner regions of haloes can be strongly affected by baryons. It is well known that feedback mechanisms from supernovae can create DM density cores too (Read & Gilmore 2005;Governato et al. 2012;Pontzen & Governato 2012;Di Cintio et al. 2013;Brooks & Zolotov 2014;Cintio et al. 2014;Pontzen & Governato 2014;Oñorbe et al. 2015;Tollet et al. 2016;Benítez-Llambay et al. 2019), but also black holes (e.g. Martizzi et al. 2013;Peirani et al. 2017;Silk 2017). Recently, Burger et al. (2022) studied degeneracies between cores produced by supernova feedback and SIDM cores. They found that the velocity dispersion profile produced by self-interactions is closer to isothermal profiles than in cores that are generated by supernova feedback in systems with bursty star formation. Nevertheless, there are properties that are less affected by baryons. We found that at large radii the halo shapes are more affected by self-interactions than the density profiles (see Appendix D). At these radii, baryons play a less important role than in the core region. Vargya et al. (2022) studied an MW-like galaxy and proposed to compare the shape of the stellar and gas component to the shape of the total matter distribution at a radial range of 2-20 kpc. Thus, the shapes at large radii (but still small enough to be affected by self-interactions) are more sensitive to DM physics than the central density or cores size because they are less influenced by baryonic physics. The same could be true for the abundance of satellites, at least as far as DM-rich satellites are concerned.
Given a specific DM model, a single property can help constrain the momentum-transfer cross-section, but is quite limited in providing bounds on the angular dependence of the differential crosssection. An observation, e.g. of the abundance of satellites, could be explained by rSIDM as well as by fSIDM with a different momentumtransfer cross-section. In combination with another property, such as the shape, it can become possible to derive bounds on the typical scattering angle as demonstrated in Section 3.9. For a given DM model, multiple measurements could lead to bounds onT/ , which are not compatible with each other and thus exclude the considered model.
In this context, we want to stress that matching cross-sections of various DM models is difficult because they typically behave qualitatively different. Strictly speaking, it is impossible for rare and frequent self-interactions and probably for further model variations too.
Several studies (Kaplinghat et al. 2016;Correa 2021;Gilman et al. 2021;Sagunski et al. 2021) have suggested that the self-interaction cross-section should be velocity-dependent, i.e. decrease for higher velocities. We have not studied this, but a velocity dependence should lead to qualitatively different results and is a reasonable step to extend our study. Moreover, velocity-dependent scattering is well motivated from a particle physics perspective. This is, in particular, the case for light mediator models, which interact frequently (e.g. Buckley & Fox 2010;Loeb & Weiner 2011;Bringmann et al. 2017).
Another important step forward would be to include baryonic matter and its feedback processes. At small radii, the mass distribution of haloes is typically dominated by the baryonic component. As explained above, feedback mechanisms such as outflows from supernovae can alter the DM distribution and also produce cored profiles. By this, they can mitigate the core-cusp problem (e.g. Brooks et al. 2013;Brooks & Zolotov 2014;Chan et al. 2015;Oñorbe et al. 2015;El-Badry et al. 2016) and the diversity problem (e.g. El-Badry et al. 2017). Thus, it is important to take baryonic effects into account when constraining the properties of DM models using these smallscale problems.

CONCLUSIONS
In this paper, we have presented the first cosmological simulations of DM with frequent self-interactions. We have compared DM-only simulations of CDM, rSIDM, and fSIDM in terms of various measures, such as the density PDF, the matter power spectrum, the twopoint correlation function of haloes and subhaloes, and the halo and subhalo mass functions. In addition, we have investigated the density and circular velocity profiles of the DM haloes as well as their shapes. Finally, we have examined qualitative differences between fSIDM and rSIDM. Our main results from these simulations are as follows: • On large scales, rSIDM and fSIDM are very similar to CDM, but deviate on small scales.
• Regarding the suppression of small-scale structures, rSIDM and fSIDM behave very similarly. This includes the power spectrum, the density PDF, and the two-point correlation function as well as the density core formation and shapes of haloes.
• We found an interestingly large suppression of the abundance of satellites in fSIDM compared to rSIDM.
• It may be possible to distinguish observationally between rSIDM and fSIDM using a combination of measurements. One promising avenue is the combination of shape or density profile measurements with the abundance of satellites. Further investigations, such as full-physics simulations, are needed to find out whether observations can discriminate between these DM models.
• Rare and frequent self-interactions behave similarly in many aspects. This often allows transferring upper limits on the crosssections of rSIDM to fSIDM.
We have conducted cosmological DM-only simulations to understand phenomenological differences between large-and small-angle DM scattering. Our results may prove helpful for more sophisticated studies that compare simulations to observations with the aim to discriminate between rSIDM and fSIDM. Such studies will include baryonic matter and baryonic physics, such as gas cooling, star formation, active galactic nuclei, and associated feedback mechanisms. This is the subject of forthcoming work.

DATA AVAILABILITY
The data underlying this paper will be shared on reasonable request to the corresponding author. (2021a), but not the second step, which re-adds the energy lost in the first one (described in section 2.2). This test problem is conducted without any other physics, i.e. gravity is not present. In Fig. A1, we show the cosmic deceleration problem by plotting the canonical momentum of the test particle as a function of the scale factor. Note that in the absence of self-interactions the canonical momentum would stay constant over the cosmic expansion. We can see that the simulation result matches the analytical prediction. Hence, we assume that the comoving integration is properly implemented.

APPENDIX B: CONVERGENCE OF DENSITY PROFILES
In Fig. B1, we study the convergence of the density profile of the most massive subhalo in the zoom-in simulation. We show the density for various resolutions and DM models as given in Tab. 2. For all three DM models, we find that the density profiles converge. However, it seems that they converge at a different speed. Comparing the two best resolved runs, CDM seems to converge the fastest, followed by fSIDM and rSIDM is the slowest. The difference in the convergence speeds might be caused by the use of random numbers to model SIDM. For fSIDM, they have a smaller influence on the particle trajectories than in rSIDM, which eventually could explain the deviation.

APPENDIX C: SUBHALO MASS FUNCTION
In Section 3.5, we studied the abundance of satellites as a function of mass. Here, we computed Fig. 6 with a smaller selection radius for the satellites. In Fig. C1, we only consider satellites within 1 vir as we have done previously in Sec. 3.9.  Figure C1. We display the number of satellites per logarithmic mass as a function of their total mass relative to the virial mass of their host. This is the same as in Fig. 6, but with a selection radius of 1 vir .

APPENDIX D: DENSITY AND SHAPE AT LARGER RADII
In this appendix, we quantitatively evaluate how much the density and shape profiles differ at larger radii between the DM models. We compute the shape of particles within elliptical shells using the tensor as given in Eq. (4). As in Sec. 3.8 we keep the volume of the shells during the iteration constant. Furthermore, we use the shells to compute the density profile too. This is in contrast to the profiles shown in Fig. 8, which were computed from spherical shells. For the computation, all particles are considered, not only those that belong to the halo as identified by . This implies that we also take the satellites into account, but it has almost no influence on the results shown here. We use the scatter between the individual objects (16th and 84th percentile) to estimate how significant the deviation between the DM models is. In Fig. D1, we show the results for the lowest mass bin of our uhr simulations. We find that the difference between the models for radii beyond the density core is somewhat larger for the shapes than the density profiles (see the bottom row). That seems to be true for ≲ 0.25 vir in the case of the smaller cross-section or ≲ 0.4 vir for the larger cross-section. For most of the radial range we covered here, the SIDM shapes are rounder than the CDM shapes. In contrast, the picture for the density profiles is less clear. The density ratios are smaller (middle row) and the differences are noisier (bottom row). For studying DM physics, this may make the measurement of shapes at larger radii preferable compared to densities. This paper has been typeset from a T E X/L A T E X file prepared by the author. . In the middle row, we show the ratio between the SIDM models and CDM. We further compute the difference of the SIDM models and CDM and divide it by the scatter among the individual systems. The result is displayed in the bottom row. The -axis is in terms of the radius that a spherical shell with the same volume has. Here, we use the mean of the two shell boundaries.