© The Authors 2025

1. Introduction

While astrophysical observations ranging from galactic to cosmological scales support the existence of dark matter (DM) via its gravitational impact, it is plausible that DM also features non-gravitational interactions as any other measured particle species does. Extensive searches for non-gravitational interactions between DM and visible matter have not found an unambiguous signal so far, and such work has pushed the limits even beyond the scale of weak interactions in many cases. In contrast, interactions of DM with itself with a sizeable strength, comparable to that known from the strong interaction, remain a viable possibility that is testable by searching for its imprints on the dynamics of galaxies and galaxy clusters.

Elastic 2 → 2 DM self-interactions have been proposed by Spergel & Steinhardt (2000) in the context of the core-cusp problem, with early N-body simulations by Burkert (2000) showing core formation followed by gravothermal collapse. Self-interacting DM (SIDM) has also been proposed for addressing other small-scale puzzles, such as the abundance and properties of satellite galaxy populations or the variability observed in galaxy rotation curves (see Bullock & Boylan-Kolchin (2017), Tulin & Yu (2018), Adhikari et al. (2022) for reviews).

While SIDM on galactic scales typically requires a cross section (per DM mass) of the order of σ/m ≳ 1 cm² g⁻¹ to have a sizeable impact, galaxy cluster observations have been argued to require upper limits reaching down to 0.1 cm² g⁻¹ (see e.g. Elbert et al. 2018; Sagunski et al. 2021; Andrade et al. 2021; Despali et al. 2022; Eckert et al. 2022; Zhang et al. 2024). Given the very different typical velocities on cluster and galactic scales, this suggests a velocity-dependent cross section (see for example Kaplinghat et al. 2016) such as occurs, for example, if the interaction is mediated by an exchange particle with mass m_ϕ that is light compared to the DM mass m_χ (Feng et al. 2009; Buckley & Fox 2010; Loeb & Weiner 2011; Tulin et al. 2013). Within this class of models, the scattering cross section also features a pronounced angular dependence, being strongly forward-dominated in the limit m_ϕ → 0, analogous to Rutherford scattering.

For isolated quasi-stationary halos, the impact of angular dependence on the density profile has been argued to be approximately captured by mapping simulations assuming isotropic scattering onto models with anisotropic differential cross sections by matching a suitable angle-averaged effective cross section, namely the viscosity cross section (e.g. Colquhoun et al. 2021; Yang & Yu 2022; Sabarish et al. 2024). However, for systems that are not quasi-stationary, such as during the evolution of satellite populations (Kahlhoefer et al. 2015; Fischer et al. 2022, 2024b; Ragagnin et al. 2024) or when approaching gravothermal collapse, and especially for strongly non-stationary processes, such as galaxy cluster mergers (e.g. Kahlhoefer et al. 2014; Robertson et al. 2017a; Fischer et al. 2021a,b, 2023; Sabarish et al. 2024), it is a priori less clear how they are affected by (strongly) anisotropic scattering.

These reasons motivated us to develop efficient numerical algorithms for incorporating differential cross sections featuring a pronounced angular dependence, such as for light mediator models, in order to assess the accuracy of (and possibly refine) commonly employed approximate ‘mapping’ prescriptions for a given observable and potentially identify new features that are specific to a given angular dependence. In this work, we develop such an algorithm, validate it, and apply it to study galaxy cluster mergers with angle-dependent cross sections that are characteristic for light mediator models.

Our work is based on a combination of two existing approaches, including one that is suitable for relatively rare individual self-scattering events, which feature large deflection angles (dubbed rare SIDM or ‘rSIDM’; see Sect. 2 and e.g. Koda & Shapiro 2011; Vogelsberger et al. 2012; Rocha et al. 2013; Fry et al. 2015; Robertson et al. 2017a; Banerjee et al. 2020 for details). In this first approach, SIDM can be modelled by collisions of N-body particles, with scattering angles sampled from the differential cross section. While being the most straightforward and direct, this method has the drawback of becoming prohibitively inefficient numerically for models in which the scattering rate is large but the momentum transfer in each individual scattering event is small.

The second approach is tailored precisely to efficiently describe such very frequent and strongly forward-dominated scatterings (dubbed frequent SIDM or ‘fSIDM’). The fSIDM approach formally corresponds to the limit in which the momentum transfer cross section is kept constant but the typical scattering angle approaches zero. Such frequent interactions are well known in the context of Coulomb interactions of energetic particles traversing a medium and can effectively be captured by a drag force complemented with transverse momentum diffusion. A drag force description has been derived by Kahlhoefer et al. (2014) and applied to merging galaxy clusters. Building on that work, an algorithm for the fSIDM limit has been proposed in the context of N-body simulations of SIDM by Fischer et al. (2021a) (respecting energy and momentum conservation). Their implementation has been developed further by Fischer et al. (2021b, 2022, 2024b) and used by Fischer et al. (2023, 2024b), Fischer & Sagunski (2024), Sabarish et al. (2024), Ragagnin et al. (2024) to study the SIDM halo evolution, galaxy cluster merger dynamics, the evolution of cluster satellites, and dynamical friction as compared to the usually considered opposite limit of isotropic rSIDM. However, as mentioned above, realistic models feature both small- and large-angle scattering with a roughly comparable impact over a wide range of parameter space, and thus neither the fSIDM nor rSIDM scheme can be considered to describe viable particle models in general.

Therefore, in this work we develop a new hybrid method that we dub the ‘hSIDM’ scheme, and it is designed to be able to accurately describe (in principle) arbitrary differential cross sections. In particular, the hSIDM method can be used to efficiently simulate models featuring a strongly forward-dominated enhancement at low deflections angles as well as a given contribution with large-angle scattering. The hSIDM method is therefore particularly suited to study light mediator models, but it can easily be extended to other scenarios that have been proposed in the literature (see e.g. Tulin & Yu 2018 for an overview). The hSIDM algorithm takes advantage of the numerically much more efficient fSIDM approach for small scattering angles below a certain critical angle θ_c and is complemented by the more conventional rSIDM treatment of large-angle scatterings.

After having developed the hSIDM scheme, we applied it to merging galaxy clusters. Systems such as the ‘Bullet Cluster’ (e.g. Springel & Farrar 2007; Mastropietro & Burkert 2008; Lage & Farrar 2014) or ‘El Gordo’ (e.g. Donnert 2014; Molnar & Broadhurst 2015; Zhang et al. 2015; Kim et al. 2021; Asencio et al. 2020, 2023) have been studied and used to constrain DM self-interactions (e.g. Valdarnini 2024). Besides other phenomena, such as the oscillations of the brightest cluster galaxy (BCG) at the late stages of the merger (e.g. Harvey et al. 2019; Cross et al. 2024), offsets between the distribution of the DM and the cluster galaxies mainly shortly after the first pericentre passage have received a lot of attention (e.g. Harvey et al. 2015; Robertson et al. 2017b; Sirks et al. 2024). Theoretical studies to model the offsets that arise from the self-interactions effectively decelerating the two DM halos when they pass through each other while the galaxies are only indirectly affected via gravity have been conducted to gain insights into their evolution and allow for constraints from observed systems to be inferred. Offsets between the galaxy and DM distributions are measured employing strong gravitational lensing to infer the total mass distribution. These measurements are challenging, and previous claims of large offsets have been questioned (e.g. Bradač et al. 2008; Dawson et al. 2012; Dawson 2013; Jee et al. 2014, 2015; Harvey et al. 2017; Peel et al. 2017; Taylor et al. 2017; Wittman et al. 2018, 2023). Until today, observed systems show no clear deviation from collisionless DM, and mainly upper bounds on the self-interaction cross section from DM-galaxy offsets exist.

This work is structured as follows. In Sect. 2, we first introduce the various exemplary differential cross sections considered in this work and review the rSIDM and fSIDM schemes for describing angle-dependent self-scattering. We introduce the new hSIDM scheme in Sect. 3 and furthermore discuss an analytically solvable test set-up that we use for validation later on. The numerical implementation is discussed in Sect. 4, and we present our validation checks in Appendix A. The hSIDM scheme is used to simulate galaxy cluster mergers for SIDM with angle dependence as predicted by light mediator models in Sect. 5. Our main results for DM-galaxy offsets and their model-dependence is discussed in Sect. 6. Finally, we summarise and conclude in Sect. 7. Additional information can be found in the Appendices.

2. Review of angle dependence in self-interacting dark matter

In this section, we first briefly review the angle dependence of typical differential cross sections relevant for light mediator models of SIDM as well as various definitions of angle-averaged cross sections considered in this context. We then review the existing rSIDM and fSIDM approaches for describing rare large-angle and frequent small-angle scatterings, respectively, on which the hybrid scheme introduced in this work is based.

2.1. Differential scattering cross sections

While the algorithm introduced in this work is capable of describing SIDM with an arbitrary differential cross section, we start by briefly reviewing some typical cases, that we also use below for illustration.

Light mediator models are characterised by scattering via a Yukawa interaction, for which DM particles of mass m_χ scatter via the exchange of a mediator with mass m_ϕ. One can further discriminate the cases in which the two incoming DM particles are identical (we refer to this case as ‘Møller scattering’ in analogy to e⁻e⁻ scattering, noting however that the mediator particle is massive here), or distinguishable (being e.g. a particle and an antiparticle, referred to as ‘Rutherford scattering’, but again representing the case of a massive mediator in this work). In the non-relativistic limit and in Born approximation, Møller scattering receives contributions from t- and u-channel diagrams, yielding

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\left(\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}\right)}_{\mathrm{M}}=\frac{1}{2}\frac{{\sigma }_0}{4\pi }& [\frac{1}{{(1+r{sin}^2\left(\frac{\theta }{2}\right))}^2}+\frac{1}{{(1+r{cos}^2\left(\frac{\theta }{2}\right))}^2}\hfill \\ \hfill & -\frac{1}{(1+r{sin}^2\left(\frac{\theta }{2}\right))(1+r{cos}^2\left(\frac{\theta }{2}\right))}],\hfill \\ \hfill \end{array}\end{array}$$(1)

with the factor of one-half accounting for the double counting of the identical particles, while only the t-channel process contributes to Rutherford scattering, giving

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\left(\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}\right)}_{\mathrm{R}}=\frac{{\sigma }_0}{4\pi }& \frac{1}{{(1+r{sin}^2\left(\frac{\theta }{2}\right))}^2}.\hfill \\ \hfill \end{array}\end{array}$$(2)

Following the notations used in Girmohanta & Shrock (2022), Robertson et al. (2017a), we introduced a parameter σ₀ that is related to the Yukawa coupling strength α_χ via σ₀ = α_χ²m_χ²/m_ϕ⁴. More importantly, the angle dependence is related to the parameter

$$\begin{array}{c}\hfill r\equiv {\left(\frac{{\beta }_{\mathrm{rel}}{m}_{\chi }}{m_{\varphi }}\right)}^2,\end{array}$$(3)

where β_rel = v_rel/c is the relative velocity of the incoming particles. The value of r characterises the dependence of the cross section on the deflection angle θ, with the isotropic limit being approached for r ≪ 1, while scattering becomes strongly forward-dominated (Coulomb-like) for large values of r. For illustration, we show the angle dependence of the Møller cross section in Fig. 1 for r = 100 and 10⁴, respectively, as well as an isotropic cross section for comparison. We emphasise the pronounced peaks towards small scattering angles (note the logarithmic y-axis). Since the Møller cross section is symmetric under θ → π − θ, it is also enhanced for θ → π, due to the fact that both scattering particles are indistinguishable.

Additionally, for validation purposes in Appendix A.2 we also consider a somewhat artificial case for which all scatterings proceed with a single, fixed deflection angle θ₀, being formally given by a differential cross section

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}\right)}_{\mathrm{D}}=\frac{{\sigma }_0}{2\pi }\delta (cos(\theta )-cos({\theta }_0)),\end{array}$$(4)

where δ(x) stands for the Dirac distribution.

2.2. Angle-averaged cross sections

In order to estimate the impact of DM self-scattering on various observables of interest, it is useful to define integrated, angle-averaged cross sections with specific weighting factors. These may also be used to map results obtained for a given angle dependence (e.g. isotropic) to other models, if the effect of interest is mainly sensitive to a particular averaged cross section. Here we briefly review the angle-averages that are often considered in the context of SIDM. We use them below to (i) characterise the relative strength of small- and large-angle scatterings in a given model, and (ii) investigate in how far our results taking the precise angle dependence into account can indeed be mapped among models featuring distinct differential cross sections, but agreeing in the value of a certain angle-averaged cross section.

The cross section averaged over the complete range of deflection angles can in general be written as

$$\begin{array}{c}\hfill {\sigma }_X\equiv c_X{\displaystyle \int \mathrm{d}}\mathrm{\Omega }{g}_X(\theta )\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}=2\pi c_X{\displaystyle {\int }_0^{\pi }}\mathrm{d}\theta g_X(\theta )\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}sin(\theta ),\end{array}$$(5)

where we assumed independence of the polar angle in the second expression (being true in all cases we consider). Here g_X(θ) is a weighting function that differs for each type of angle-average labelled by the index X, and c_X is a normalisation constant chosen such that c_X∫dΩ g_X(θ) = 4π. This ensures that all angle-averages agree with each other for an isotropic cross section¹.

We considered the following cases: (1) Total cross section, which determines the total scattering rate, where

$$\begin{array}{c}\hfill g_{\mathrm{tot}}(\theta )\equiv 1,c_{\mathrm{tot}}=1;\end{array}$$(6)

(2) Transfer cross section, which weights with momentum transfer for scattering of distinguishable particles, where

$$\begin{array}{c}\hfill g_{\mathrm{T}}(\theta )\equiv 1-cos(\theta ),c_{\mathrm{T}}=1;\end{array}$$(7)

(3) Modified transfer cross section, which weights with momentum transfer for scattering of identical particles, where

$$\begin{array}{c}\hfill g_{\stackrel{\sim }{\mathrm{T}}}(\theta )\equiv 1-|cos(\theta )|,c_{\stackrel{\sim }{\mathrm{T}}}=2;\end{array}$$(8)

(4) Viscosity cross section, which weights with energy transfer (for both distinguishable and identical particles), where

$$\begin{array}{c}\hfill g_{\mathrm{V}}(\theta )\equiv 1-{cos}^2(\theta ),c_{\mathrm{V}}=\frac{3}{2}.\end{array}$$(9)

The viscosity cross sections has been argued to be a good proxy for the impact of self-scattering on the distribution of DM in isolated halos, including core formation and collapse (e.g. Colquhoun et al. 2021; Yang & Yu 2022; Sabarish et al. 2024). The question whether an average cross section is a good description in mergers has also been investigated (e.g. Kahlhoefer et al. 2014; Robertson et al. 2017a; Fischer et al. 2021a,b). Below we scrutinise these findings for the evolution of DM and galaxy positions as well as their offset in cluster mergers, by comparing to a treatment capturing the exact angle dependence. The total cross section is mainly relevant for estimating the numerical effort (within the rSIDM scheme, see below) since it characterises the total scattering rate, irrespective of whether a significant amount of momentum is transferred between the scattering partners.

For reference, we give the angle-averages for the two extreme cases of isotropic scattering and strongly forward-dominated scattering, respectively,

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{tot}}& ={\sigma }_{\mathrm{T}}={\sigma }_{\stackrel{\sim }{\mathrm{T}}}={\sigma }_{\mathrm{V}},(\mathrm{isotropic})\hfill \\ \hfill {\sigma }_{\mathrm{tot}}& \gg {\sigma }_{\mathrm{T}}=\frac{1}{3}{\sigma }_{\mathrm{V}},(\mathrm{forward}-\mathrm{dom}.,\mathrm{distinguishable})\hfill \\ \hfill {\sigma }_{\mathrm{tot}}& \gg {\sigma }_{\stackrel{\sim }{\mathrm{T}}}=\frac{2}{3}{\sigma }_{\mathrm{V}},(\mathrm{forward}-\mathrm{dom}.,\mathrm{identical})\hfill \end{array}$$(10)

where the last two lines refer to the cases of either distinguishable or identical particles, with ‘forward-dominated’ meaning a strong enhancement for θ → 0 as well as θ → π in the latter case.

2.3. Existing approaches for angle-dependent self-interacting dark matter

Here we briefly review the physical basis of the approaches that have been used to incorporate angle dependence into SIDM simulations. Technical details on the numerical implementation are provided in Sect. 4.

As mentioned above, we consider two existing approaches, both of which are combined in this work in order to be able to describe models with arbitrary differential cross sections. However, when used by themselves, each of the two methods is limited to a particular class of models.

The first method treats angle dependence explicitly for each scattering event, by sampling the deflection angle from a random distribution that follows from the differential cross section (see e.g. Robertson et al. 2017a; Banerjee et al. 2020). This strategy is in practice numerically feasible for models for which an 𝒪(1) amount of energy and momentum is transferred in each scattering, thus being limited to large-angle scatterings, see Robertson et al. (2017a). For these types of models the total and (modified) transfer cross sections are of comparable size, ${\sigma }_{\mathrm{tot}}/{\sigma }_{\stackrel{\sim }{\mathrm{T}}}\sim \mathcal{O}(1)$, which implies that even relatively ‘rare’ scattering events can have a sizeable impact. We thus refer to this case as rare self-interacting dark matter (rSIDM).

The second method employs an effective treatment that is applicable if DM scatters very frequently, but the momentum transfer in each scattering is small, dubbed frequent self-interacting dark matter (fSIDM)following Fischer et al. (2021a). By itself, this method becomes exact for models approaching the formal limit ${\sigma }_{\mathrm{tot}}/{\sigma }_{\stackrel{\sim }{\mathrm{T}}}\to \infty$, that is, very strongly forward-dominated scattering only.

In this work we combine both methods, using them within the range of angles for which they work best. Before discussing this hybrid approach, we briefly review the physical basis of the fSIDM scheme (see Kahlhoefer et al. 2014 and Kummer et al. 2018 for more details). In analogy to the seminal description of Brownian motion by Einstein (1905) and von Smoluchowski (1906), the collective effect of frequent small-angle scattering can be described by a drag force that acts along the direction of motion (more precisely the relative velocity vector) as well as a stochastic force leading to diffusive motion. The drag force can be related to the (modified) transfer cross section, $F_{\mathrm{drag}}\propto {\sigma }_{\stackrel{\sim }{\mathrm{T}}}$ (see Sect. 4 for details), while the stochastic force can be included as a kick in a random direction transverse to the relative velocity vector, with magnitude consistent with momentum and energy conservation.

This set-up has been implemented by Fischer et al. (2021a) and verified with multiple test problems including one for which the solution is known from Molière (1948), being essentially a Gaussian spreading of an initially collimated beam traversing a medium of target particles. In this work we generalise this test problem to models with arbitrary differential cross section in order to validate the hybrid scheme based on combining the rSIDM and fSIDM methods, which we turn to next.

3. Hybrid scheme for angle dependence

In this section, we first present an algorithm that allows us to simulate DM self-scattering with arbitrary differential cross sections efficiently, dubbed hybrid self-interacting dark matter or hSIDM. Next, we review a test set-up featuring an exact analytical solution for any given differential cross section and then use it to assess under which conditions and for which choices of various technical parameters the hSIDM approach is expected to be valid.

3.1. Hybrid approach to angle-dependent self-interacting dark matter

Here we introduce a new hybrid approach for taking the dependence of the DM self-interaction cross section on the deflection angle into account, that we refer to as hSIDM as stated above. The central idea is to combine the methods used previously to describe models with either only large (rSIDM) or only forward-dominated (fSIDM) scatterings, allowing for both regimes to be described efficiently within a single set-up. Since typical (e.g. light mediator) models can feature a very strong enhancement of the total scattering rate at small angles, (see e.g. Fig. 1) but also predict non-negligible large-angle scattering, it is necessary to capture both contributions in order to obtain accurate predictions for a given set of model parameters.

The central idea of the hybrid scheme is to introduce a critical angle, called θ_c, and use the effective drag force and transverse momentum diffusion (i.e. the fSIDM approach) for all scatterings with θ ≤ θ_c or θ ≥ π − θ_c. For angles within the interval θ_c < θ < π − θ_c, the rSIDM approach is used instead, i.e. a random sampling of deflection angles from the differential cross section. Thus, we write

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\left(\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}\right)}_{\mathrm{hSIDM}}=& {\left(\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}\right)}_{\mathrm{fSIDM}}[\mathrm{\Theta }({\theta }_{\mathrm{c}}-\theta )+\mathrm{\Theta }(\theta -(\pi -{\theta }_{\mathrm{c}}))]\hfill \\ \hfill & +{\left(\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}\right)}_{\mathrm{rSIDM}}[1-\mathrm{\Theta }({\theta }_{\mathrm{c}}-\theta )-\mathrm{\Theta }(\theta -(\pi -{\theta }_{\mathrm{c}}))],\hfill \end{array}\end{array}$$(11)

where Θ(θ) is the Heaviside function. The splitting is also illustrated by the dark and light shaded regions in Fig. 1 for two typical values θ_c = 0.1 and 0.3. In the following we refer to the range θ_c < θ < π − θ_c as ‘large-angle’ and to the ranges θ ≤ θ_c or θ ≥ π − θ_c as ‘small-angle’ regimes, noting that deflection angles close to π are also considered as being ‘small’ for identical particles².

Since the critical angle θ_c is a purely technical parameter, all physical quantities should be independent of its choice, and we check this explicitly for all our results below. Theoretically, we expect this is the case provided θ_c is chosen small enough such that the small-angle effective description underlying the fSIDM approach can be applied. As we see later on, the appropriate size of θ_c depends on characteristics of the set-up, such as the characteristic dynamical time and the considered time interval, in addition to the angle dependence of the cross section itself. We quantify the validity range for the choice of θ_c within the hSIDM scheme in Sect. 3.2, employing a test set-up that can be described analytically. In addition, we point to Appendix A.3 for a discussion of requirements constraining the choice of θ_c within typical simulation set-ups, and to Appendix A.4 for a demonstration of the performance gain of the hybrid scheme compared to a direct implementation of angular dependence.

The main reason for considering the hSIDM scheme is that the effective small-angle treatment is numerically much more efficient to capture the forward-dominated regime. We provide an explicit comparison of the speed-up when using hSIDM (as compared to a direct sampling of scattering angles over the entire range) in Appendix A.4 below. However, it is also possible to estimate the expected speed-up. For that purpose, it is convenient to consider the contributions to the various angle-averaged cross sections from Sect. 2.2 within the ‘small’ and ‘large’-angle regimes depending on the choice of θ_c. We thus generalise Eq. (5) and define

$$\begin{array}{cc}\hfill {\sigma }_{X,<}({\theta }_{\mathrm{c}})& \equiv c_X{\displaystyle \underset{\theta \in I_{<}}{\int }}\mathrm{d}\mathrm{\Omega }{g}_X(\theta )\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }},I_{<}=[0,{\theta }_{\mathrm{c}}]\cup [\pi -{\theta }_{\mathrm{c}},\pi ],\hfill \\ \hfill {\sigma }_{X,>}({\theta }_{\mathrm{c}})& \equiv c_X{\displaystyle \underset{\theta \in I_{>}}{\int }}\mathrm{d}\mathrm{\Omega }{g}_X(\theta )\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }},I_{>}=({\theta }_{\mathrm{c}},\pi -{\theta }_{\mathrm{c}}),\hfill \end{array}$$(12)

allowing us to decompose the total and (modified) transfer cross sections into the contributions treated with the fSIDM (< ) and rSIDM (> ) method within hSIDM, respectively,

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{tot}}& ={\sigma }_{\mathrm{tot},<}({\theta }_{\mathrm{c}})+{\sigma }_{\mathrm{tot},>}({\theta }_{\mathrm{c}}),\hfill \\ \hfill {\sigma }_{\stackrel{\sim }{\mathrm{T}}}& ={\sigma }_{\stackrel{\sim }{\mathrm{T}},<}({\theta }_{\mathrm{c}})+{\sigma }_{\stackrel{\sim }{\mathrm{T}},>}({\theta }_{\mathrm{c}}).\hfill \end{array}$$(13)

The drag force used for the small-angle regime is related to the (modified) transfer cross section. Within hSIDM, only scatterings with θ ∈ I_< are treated in this way, and thus (see Sect. 4 for details)

$$\begin{array}{c}\hfill F_{\mathrm{drag}}({\theta }_{\mathrm{c}})\propto {\sigma }_{\stackrel{\sim }{\mathrm{T}},<}({\theta }_{\mathrm{c}})(\mathrm{small}-\mathrm{angle}\mathrm{regime}).\end{array}$$(14)

The numerical costs depend on the size of the numerical time step. Within the fSIDM scheme, the time step is proportional to the drag force, F_drag, and hence also proportional to the (modified) transfer cross section, ${\sigma }_{\stackrel{\sim }{\mathrm{T}}}$. In contrast, the rSIDM scheme, for which individual collisions are modelled numerically as scattering of N-body particles (see Sect. 4), has a time step proportional to the total scattering rate, R, which is proportional to the total cross section of large-angle scattering, i.e. within the hSIDM scheme it is

$$\begin{array}{c}\hfill R({\theta }_{\mathrm{c}})\propto {\sigma }_{\mathrm{tot},>}({\theta }_{\mathrm{c}})(\mathrm{large}-\mathrm{angle}\mathrm{regime}).\end{array}$$(15)

We note that for light mediator models (i.e. ‘Møller scattering’), the strong enhancement of the forward-scattering rate implies that for typical values of θ_c ∼ 0.1 − 0.3 one has σ_tot, > (θ_c)≪σ_tot. The expected speed-up, S, within hSIDM as compared to a naive direct treatment of the angle dependence can thus be estimated as

$$\begin{array}{c}\hfill S\propto \frac{{\sigma }_{\mathrm{tot}}}{{\sigma }_{\mathrm{tot},>}({\theta }_{\mathrm{c}})}.\end{array}$$(16)

The dependence of this ratio on θ_c for the Møller scattering cross section Eq. (1) is shown in Fig. 2 (dashed lines), for various values of the anisotropy parameter r from Eq. (3). Even for moderate choices θ_c = 0.1 − 0.3, one has σ_tot, > (θ_c)/σ_tot ≪ 1, indicating a large expected speed-up within hSIDM. In practice, we find that the speed-up is even slightly larger than this estimate (see Appendix A.4).

On the other hand, the physical relevance of small- versus large-angle scatterings can be estimated from comparing the (modified) transfer cross sections ${\sigma }_{\stackrel{\sim }{\mathrm{T}},<}({\theta }_{\mathrm{c}})$ and ${\sigma }_{\stackrel{\sim }{\mathrm{T}},>}({\theta }_{\mathrm{c}})$. The fraction ${\sigma }_{\stackrel{\sim }{\mathrm{T}},<}({\theta }_{\mathrm{c}})/{\sigma }_{\stackrel{\sim }{\mathrm{T}}}$ related to the small-angle regime is shown by the solid lines in Fig. 2. Thus, for typical values of θ_c both small- and large-angle scatterings are expected to have a roughly comparable physical impact.

In summary, the main advantage of the hSIDM regime is that it allows us to treat very frequent small-angle scatterings much more efficiently as compared to a naive sampling of deflection angles, while in addition explicitly including the angle dependence of the physically roughly equally important contribution from large-angle scattering. With the hSIDM scheme simulation of fairly anisotropic cross sections with significant contributions from both small and large angles that have not been feasible before, become possible thanks to the efficient treatment of small-angle scattering.

3.2. Applicability of effective description for small-angle scattering: Deflection test set-up

In order to gain insight on the expected range of validity of the hSIDM scheme, and especially on the choice of the critical angle θ_c, we investigate the applicability of the effective treatment of small-angle scattering. For that purpose, it is instructive to consider a simple test set-up that admits an analytical treatment. The purpose of considering this set-up is (i) to derive analytical conditions on the angle dependence of the cross section, as well as the typical local dark matter density and velocity and dynamical timescale for which the effective treatment of small-angle scattering is valid (see Sect. 3.2.2), (ii) translate these into validity conditions for the admissible range of values of θ_c within the hSIDM scheme (see Sect. 3.2.3), (iii) demonstrate the resulting speed-up of the numerical algorithm based on hSIDM, and (iv) provide a stringent validation test for our numerical implementation (see Sect. A below). Moreover, the test set-up can also be seen as a proxy how momentum and energy is transferred in a collision of two “patches” of the dark matter population within the two halos that collide in the initial stages of a galaxy cluster merger, provided the relative velocity is large compared to the velocity dispersion with each of the patches.

3.2.1. Deflection test set-up

Following Bethe (1953), we consider a classic deflection problem of a beam of test particles traversing a homogeneous medium of target particles with number density n. All beam particles move initially in the same direction (say along the z-axis) with a given velocity v. After some time t, scatterings of the beam on target particles lead to a ‘broadening’ of the beam. This can be described by a distribution function f(t, θ), such that f(t, θ)⋅sin(θ) dθ measures the fraction of beam particles with azimuthal angle θ to the z-axis within the interval (θ, θ + dθ) at time t (assuming rotational symmetry around the z-axis). The initial condition at t = 0 is thus

$$\begin{array}{c}\hfill f(t=0,\theta )=\delta (cos(\theta )-1).\end{array}$$(17)

We study this deflection set-up in two variants, (i) in the limit for which no recoil is transmitted to the target particles (i.e. formally infinitely heavy target particles), and (ii) when including recoil and assuming equal masses for beam and target particles. We note that in this validation set-up scatterings occur only between beam and target particles, and that no gravitational force is included.

In the remainder of this subsection, we focus on case (i), for which the magnitude of the velocity is unchanged and only the direction of the test particles changes in each scattering³. In this case, a full analytical solution for the Boltzmann equation determining the distribution function is known, going back to Goudsmit & Saunderson (1940a,b). Following the more explicit derivation by Bethe (1953), one finds

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill f(t,\theta )=& {\displaystyle \sum _{\ell =0}^{\infty }}{P}_{\ell }(cos(\theta ))(\ell +\frac{1}{2})\hfill \\ \hfill & exp\{-nvt{\displaystyle \int \mathrm{d}}\mathrm{\Omega }\prime \frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }\prime }[1-P_{\ell }(cos(\theta \prime ))]\},\hfill \end{array}\end{array}$$(18)

where P_ℓ(x) are Legendre polynomials, and dσ/dΩ is the differential scattering cross section for beam-target particle scattering. We note that the total number of scatterings each beam particle experiences on average between t = 0 and time t is given by the total opacity

$$\begin{array}{c}\hfill {\tau }_{\mathrm{tot}}(t)\equiv nvt{\sigma }_{\mathrm{tot}}=nvt{\displaystyle \int \mathrm{d}}\mathrm{\Omega }\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }},\end{array}$$(19)

and increases linearly with time t.

Let us first discuss some general features. Using the properties ∫₀^πdθ sin(θ)P₀(cos(θ)) = 2 and ∫₀^πdθ sin(θ)P_ℓ(cos(θ)) = 0 for ℓ > 0 one sees that

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{\pi }}\mathrm{d}\theta sin(\theta )f(t,\theta )=1,\end{array}$$(20)

for all times t, i.e. the distribution is properly normalised at all times. Furthermore, the relation

$$\begin{array}{c}\hfill {\displaystyle \sum _{\ell }}{P}_{\ell }(x)(\ell +1/2)=\delta (x-1),\end{array}$$(21)

with x ≡ cos(θ), ensures that the correct initial condition is reached for t → 0. Additionally, since P_ℓ(cos(θ)) ≤ 1 (and strictly less than unity for θ > 0 and ℓ > 0), all summands in Eq. (18) are exponentially damped for t → ∞ except for ℓ = 0, implying

$$\begin{array}{c}\hfill f(t\to \infty ,\theta )\to 1/2,\end{array}$$(22)

approaches an isotropic distribution at late times, as expected after a large number of multiple scatterings.

The distribution thus evolves from a peaked to a flat shape, and the differential cross section determines how the intermediate evolution occurs. For example, for an isotropic cross section one has

$$\begin{array}{cc}\hfill f(t,\theta ){|}_{\mathrm{iso}}& =\frac{1}{2}+{\displaystyle \sum _{\ell >0}}{P}_{\ell }(cos(\theta ))(\ell +\frac{1}{2})e^{-{\tau }_{\mathrm{tot}}(t)}\hfill \\ \hfill & =\frac{1}{2}(1-e^{-{\tau }_{\mathrm{tot}}(t)})+\delta (cos(\theta )-1)e^{-{\tau }_{\mathrm{tot}}(t)}.\hfill \end{array}$$(23)

Thus for isotropic scattering the distribution is given by a superposition of a flat and a peaked component at any time, with time-dependent relative weighting factors 1 − e^−τ_tot and e^−τ_tot, respectively. We note that the exponential damping factor in Eq. (18) is identical for all ℓ > 0 in the isotropic case.

This changes completely when the differential cross section is strongly forward-dominated. In this case the Legendre polynomial inside the integral over dσ/dΩ′ in Eq. (18) can be approximated for small deflection angle as

$$\begin{array}{c}\hfill 1-P_{\ell }(x)=\frac{1}{2}\ell (\ell +1)(1-x)+\mathcal{O}{(1-x)}^2,\end{array}$$(24)

which yields

$$\begin{array}{c}\hfill f(t,\theta ){|}_{\frac{\mathrm{forward}}{\mathrm{dominated}}}={\displaystyle \sum _{\ell \ge 0}}{P}_{\ell }(cos(\theta ))(\ell +\frac{1}{2})exp[-\frac{1}{2}\ell (\ell +1){\tau }_{\mathrm{T}}(t)],\end{array}$$(25)

where higher ℓ become more and more suppressed. The relevant quantity determining the evolution is the transfer cross section σ_T, see Sect. 2.2, that we use to define a transfer opacity via

$$\begin{array}{c}\hfill {\tau }_{\mathrm{T}}(t)\equiv nvt{\sigma }_{\mathrm{T}}=nvt{\displaystyle \int \mathrm{d}}\mathrm{\Omega }(1-cos(\theta ))\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}.\end{array}$$(26)

This already indicates that for purely forward-dominated scatterings, the transfer and not the total cross section is the relevant quantity for the dynamical evolution (see Kahlhoefer et al. 2014)⁴.

The approximation given by Eq. (24) is the essential assumption underlying the effective treatment of small-angle scattering entering the fSIDM approach. This can be seen more directly when assuming in addition that τ_T is small enough such that sufficiently many ℓ contribute in the summation in Eq. (25). In this case the summation over ℓ can be replaced by a continuous integration using the Euler formula (see e.g. Bethe 1953), and Eq. (25) yields the Molière approximation (compare with Molière 1948)

$$\begin{array}{c}\hfill f(t,\theta )\approx \frac{1}{{\tau }_{\mathrm{T}}(t)}exp[-\frac{{\theta }^2}{2{\tau }_{\mathrm{T}}(t)}].\end{array}$$(27)

The effect of multiple small-angle scattering thus yields a Gaussian broadening of the beam, as long as τ_T(t) ≲ 𝒪(1). This is in direct correspondence to the effective fSIDM treatment of forward-dominated scatterings in terms of a transverse momentum diffusion process, which also leads to a Gaussian broadening when applied to the deflection set-up (with identical width, see Fischer et al. 2021a)⁵.

3.2.2. Validity of the fSIDM approach

Importantly, the analytically solvable deflection set-up allowed us to derive a general criterion for the validity of the effective forward-dominated approximation. In this subsection, we start by discussing the validity in the case of allowing for the entire range of possible scattering angles, i.e. the fSIDM approximation. Subsequently, we generalise the discussion to the hSIDM case, for which only angles below a critical angle are included in the effective small-angle approximation.

To obtain a quantitative estimate of the corrections to the forward-dominated approximation, we keep the next term in the Taylor expansion of Eq. (24), of order (1 − x)² (where x ≡ cos(θ)), given by

$$\begin{array}{cc}\hfill 1-P_{\ell }(x)& =\frac{\ell (\ell +1)}{2}(1-x)-\frac{(\ell -1)\ell (\ell +1)(\ell +2)}{16}{(1-x)}^2\hfill \\ \hfill & +\mathcal{O}{(1-x)}^3.\hfill \end{array}$$(28)

The fSIDM approach is valid if the contribution of the (1 − x)² term to the distribution function f(t, θ) (entering via the exponential in Eq. (18)) is small. This suggests to consider a ‘transfer-squared’ cross section and a corresponding opacity given by

$$\begin{array}{c}\hfill {\tau }_{{\mathrm{T}}^2}(t)\equiv nvt{\sigma }_{{\mathrm{T}}^2},{\sigma }_{{\mathrm{T}}^2}\equiv {\displaystyle \int \mathrm{d}}\mathrm{\Omega }{(1-cos(\theta ))}^2\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}.\end{array}$$(29)

Including the correction term, the exponential factor in the forward-dominated approximation of Eq. (25) becomes

$$\begin{array}{c}\hfill exp[-\frac{1}{2}\ell (\ell +1){\tau }_{\mathrm{T}}(t)+\frac{(\ell -1)\ell (\ell +1)(\ell +2)}{16}{\tau }_{{\mathrm{T}}^2}(t)(t)+\cdots ],\end{array}$$(30)

where the ellipsis stand for even higher-order terms. The effective small-angle approach is thus valid if the second summand in the exponential is small compared to the first one, i.e. τ_T² ≪ 8τ_T/ℓ², for all values of ℓ that yield a sizeable contribution in the sum over ℓ in Eq. (25). Since the exponential factor becomes strongly suppressed when ℓ(ℓ + 1)τ_T/2 ≫ 1, the relevant range is up to ℓ² ∼ 2/τ_T. This yields the validity condition

$$\begin{array}{c}\hfill {\tau }_{{\mathrm{T}}^2}(t)\ll 4{({\tau }_{\mathrm{T}}(t))}^2.\end{array}$$(31)

In practice, this means the small-angle treatment becomes valid only after a finite time t_min, for which sufficient scatterings have taken place, given by

$$\begin{array}{c}\hfill t\gg t_{\min }=\frac{1}{\rho v\frac{{\sigma }_{\mathrm{T}}}{m_{\chi }}}\frac{{\sigma }_{{\mathrm{T}}^2}}{4{\sigma }_{\mathrm{T}}}=t_{\mathrm{dyn}}\frac{{\sigma }_{{\mathrm{T}}^2}}{4{\sigma }_{\mathrm{T}}},\end{array}$$(32)

where we expressed the number density in terms of the energy density using ρ = m_χn. For the effective small-angle approach to be useful, t_min should be much smaller than the timescale on which the distribution reaches its late-time isotropic limit, which is of order of the typical dynamical evolution timescale $t_{\mathrm{dyn}}\equiv 1/(\rho v\frac{{\sigma }_{\mathrm{T}}}{m_{\chi }})$, where v is the beam velocity. Thus, we arrive at the condition for the effective small-angle treatment to be applicable and useful, given by

$$\begin{array}{c}\hfill {\sigma }_{{\mathrm{T}}^2}\ll {\sigma }_{\mathrm{T}}.\end{array}$$(33)

We note that this is the condition that is relevant when using the effective fSIDM approach for all possible deflection angles, being in practice limited to models that exclusively feature forward-dominated scattering. For realistic models, this is typically not the case. As an example, Møller scattering with r = 1, 10, 100, 10³, 10⁴ yields σ_T²/σ_T = 0.32, 0.23, 0.14, 8.5 × 10⁻², 5.9 × 10⁻² respectively. Analogously Rutherford scattering yields 1.2, 0.82, 0.51, 0.33, 0.24. This motivates the use of the hybrid approach instead. We discuss its validity next.

3.2.3. Validity of the hSIDM approach

Let us now discuss how the validity of the effective treatment of forward-dominated scattering is changed within the hSIDM scheme, i.e. when only including deflection angles below a certain critical angle θ_c (and above π − θ_c for indistinguishable particles) in the small-angle approximation. We start the discussion for the test set-up with distinguishable particles and no recoil (i.e. infinitely heavy target particles), for which the analytical solution of the deflection set-up has been given in Eq. (18) above, and then generalise to the case of identical particles and with recoil (i.e. identical masses for beam and target particles).

To obtain the hSIDM validity conditions, we replace the cross section integrated over all angles by the cross section integrated only over the small-angle regime, see Eq. (12), in the validity conditions obtained above. Thus, hSIDM validity requires

$$\begin{array}{c}\hfill {\tau }_{{\mathrm{T}}^2,<}(t,{\theta }_{\mathrm{c}})\ll 4{({\tau }_{\mathrm{T},<}(t,{\theta }_{\mathrm{c}}))}^2,\end{array}$$(34)

where τ_X, < (t, θ_c)≡nvtσ_X, < (θ_c) for X = T (transfer cross section) and X = T² (transfer-squared cross section, defined analogously to Eq. (29) but integrated only over the angle-interval I_< ). As above, this implies a minimal timescale after which the hSIDM approach is valid, given by

$$\begin{array}{c}\hfill t\gg t_{\min }({\theta }_{\mathrm{c}})=\frac{1}{\rho v\frac{{\sigma }_{\mathrm{T},<}}{m_{\chi }}}\frac{{\sigma }_{{\mathrm{T}}^2,<}}{4{\sigma }_{\mathrm{T},<}}=t_{\mathrm{dyn}}\frac{{\sigma }_{{\mathrm{T}}^2,<}}{4{\sigma }_{\mathrm{T},<}}.\end{array}$$(35)

Thus the hSIDM approach is valid and advantageous if

$$\begin{array}{c}\hfill {\sigma }_{{\mathrm{T}}^2,<}({\theta }_{\mathrm{c}})\ll {\sigma }_{\mathrm{T},<}({\theta }_{\mathrm{c}}).\end{array}$$(36)

The conditions from above can be easily generalised to the case of identical particles with equal mass by using the appropriate definition of the cross sections σ_X, < (θ_c) for that case, (a) integrated over the interval corresponding I_<  = [0, θ_c]∪[π − θ_c, π] instead of I_<  = [0, θ_c], and (b) replacing the transfer (X = T) by the modified transfer ($X=\stackrel{\sim }{\mathrm{T}}$) cross section, see Sect. 2.2, and analogously for X = T².

As an illustrative example, for Møller (Rutherford) scattering with r = 10³, one obtains σ_T², < (θ_c)/σ_T, < (θ_c) = 1.3 × 10⁻³ (2.6 × 10⁻³) for θ_c = 0.1. We provide a detailed comparison of hSIDM simulations to a full, explicit treatment of the angle dependence for various models in Sect. A, where we also cross-check the theoretically expected validity conditions derived here.

We expect that the condition on the cross section Eq. (36) is valid also beyond the deflection set-up. We note that it is independent of the details of the test problem and only involves the angle dependence of the cross section itself, and thus Eq. (36) depends only on the properties of the particle physics model underlying the angle dependence of the self-scattering cross section. Moreover, as discussed above, the deflection test can be viewed as a proxy for the dynamics of two colliding patches of dark matter in a galaxy cluster merger, such that the main physical mechanism of self-interactions on the dark matter distribution is broadly similar. This supports the applicability of Eq. (36) in cluster merger simulation set-ups. We also expect that Eq. (35) can be generalised beyond the test problem when replacing t_dyn on the right-hand side by the appropriate dynamical timescale on which the system of interest evolves.

4. Implementation

We implement the hSIDM scheme for efficiently describing DM self-interactions with arbitrary differential cross section in OPENGADGET3 (see Groth et al. 2023, and the references therein), an upgraded version of the N-body code GADGET-2 (Springel 2005), following existing implementations for either purely forward-dominated (fSIDM) or large-angle scattering (rSIDM), respectively. In Sect. 4.1 we briefly review existing implementations, and then introduce the numerical algorithm underlying this work in Sect. 4.2. We discuss criteria for choosing the associated numerical time steps in Sect. 4.3.

4.1. Existing implementations

In this section, we describe the numerical algorithm for the two existing methods reviewed in Sect. 2.3. We start with rSIDM for large-angle scattering and then turn to the purely forward-dominated fSIDM case.

4.1.1. Algorithm for rare large-angle scattering: rSIDM

N-body simulations of DM self-interactions via large-angle scatterings, dubbed rSIDM, are well-established. Most studies assume the simplest case of isotropic scattering, but also an explicit sampling of deflection angles from a given distribution has been considered (see e.g. Robertson et al. 2017b; Banerjee et al. 2020). Microscopic two-body interactions of DM particles are modelled by collisions between numerical N-body particles within rSIDM. Various approaches along these lines exist, differing in the computation of scattering probabilities (e.g. Koda & Shapiro 2011; Vogelsberger et al. 2012; Rocha et al. 2013; Robertson et al. 2017b). Starting point is the interaction between two phase-space patches, represented by two numerical particles assumed to have equal mass. To simulate isotropic rSIDM, Fischer et al. (2021a) utilised the method presented by Rocha et al. (2013). The kernel function, W, characterises the density distribution of DM for a numerical particle in configuration space. Specifically, the scattering probability of a numerical particle i to scatter with a particle j having the mass m_j can be expressed as

$$\begin{array}{c}\hfill P_i=\frac{{\sigma }_{\mathrm{tot}}}{m_{\chi }}{m}_j|\mathrm{\Delta }{\mathit{v}}_{\mathit{ij}}|\mathrm{\Delta }t{\mathrm{\Lambda }}_{\mathit{ij}},\end{array}$$(37)

where Δv_ij is the relative velocity of the two numerical particles, Δt denotes the time step, and Λ_ij = ∫dVW(|x − x_i|,h_i) W(|x − x_j|,h_j) gives the kernel overlap, with h being the kernel size. Following Fischer et al. (2021a), we employ a spline kernel (Monaghan & Lattanzio 1985) and choose h adaptively such that it includes the N_ngb next neighbours. Based on the probability P_i, we determine whether two numerical particles scatter during the time step. If a randomly selected number x from the interval [0, 1] satisfies x ≤ P_i, the particles scatter. Within the rSIDM scheme, the scattering angle is chosen randomly according to the differential cross section. For that purpose, it is convenient to consider the cumulative probability (as e.g. by Robertson et al. 2017b)

$$\begin{array}{c}\hfill P(\theta )={\displaystyle {\int }_0^{\theta }}\mathrm{d}\theta \prime \frac{2\pi sin(\theta \prime )}{{\sigma }_{\mathrm{tot}}}\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }\prime }\in [0,1]\mathrm{for}\theta \in [0,\pi ],\end{array}$$(38)

which gives the probability of scattering by an angle less than θ. The deflection angle can then be determined by drawing another random number y ∈ [0, 1] with uniform distribution, and solving P(θ) = y. We note that the deflection angle refers to the centre of mass frame, such that the computation of the vectorial velocities after the collision involves (i) a transformation of the initial velocity vectors into the centre of mass frame, (ii) application of the azimuthal deflection angle θ along with a random, uniform polar angle ϕ ∈ [0, 2π], and (iii) an inverse transformation of the modified velocity vectors back to the original ‘laboratory’ frame. For the case of isotropic scattering, the velocity directions after the scattering can simply be chosen randomly in the centre of mass frame. Details on how the vectorial velocities are calculated in practice, are given in Appendix F.

4.1.2. Algorithm for frequent small-angle scattering: fSIDM

The effective description of frequent, purely forward dominated scattering has been implemented by Fischer et al. (2021a). The self-interaction is described by a drag force and a diffusive random ‘kick’. As opposed to rSIDM, all pairs of numerical particles that are close enough to each other are affected by the interaction within a given time step. Here ‘close enough’ means in practice that the kernel overlap is positive (i.e. Λ_ij > 0). The individual ‘interaction’ of two numerical particles i and j within fSIDM is then composed of two steps. Firstly, as already discussed in Sect. 2.3, a drag force proportional to the (modified) transfer cross section is applied, with magnitude given by (see Eq. (9) by Fischer et al. 2021a)⁶

$$\begin{array}{c}\hfill F_{\mathrm{drag}}=\frac{1}{4}{|\mathrm{\Delta }{\mathit{v}}_{\mathit{ij}}|}^2\frac{{\sigma }_{\stackrel{\sim }{\mathrm{T}}}}{m_{\chi }}{m}_i{m}_j{\mathrm{\Lambda }}_{\mathit{ij}}.\end{array}$$(39)

The drag force acts along the direction of the relative velocity vector Δv_ij, and decelerates the numerical particles, reducing the velocity of particle i by an amount Δv_drag = (F_drag/m_i)Δt. The second step involves a random ‘kick’ in the direction transverse to the relative velocity vector, in a random direction within the transverse plane. This step corresponds to the diffusive part of the effective description of small-angle scattering. The magnitude of the transverse velocity kick is determined such that total energy is conserved, in other words, it makes up for the energy lost due to the drag force, given by

$$\begin{array}{c}\hfill \frac{2\mathrm{\Delta }E}{m}=|\mathrm{\Delta }{\mathit{v}}_{\mathrm{drag}}|\left(|\mathrm{\Delta }{\mathit{v}}_{\mathit{ij}}|-|\mathrm{\Delta }{\mathit{v}}_{\mathrm{drag}}|\right).\end{array}$$(40)

This two-step algorithm ensures that energy as well as three-momentum is explicitly conserved. We note that the deceleration due to the drag force and the transverse momentum kick are computed in the centre of mass frame, to determined the updated velocity vectors in the ‘laboratory’ frame.

4.2. Implementation of the hybrid approach: hSIDM

In this work we implement the hybrid approach for DM self-interaction as described by an arbitrary differential cross section. For that purpose, we combine the approaches describe above, using the drag force and momentum diffusion method (fSIDM) for scattering with ‘small angle’, and the explicit sampling technique (rSIDM) for ‘large-angle’ scatterings. Both are separated by a critical angle θ_c.

In each time step, for every pair with a positive kernel overlap Λ_ij, we simulate the small-angle scattering, before we model large scattering angles. Moreover, the pairwise computation of the particle interactions is done in a consecutive manner to conserve energy explicitly. This means that a numerical particle cannot interact with multiple particles at the same time. Only after small- and large-angle scattering have been computed for the particles of a given pair, do further pairwise interactions with other particles follow. When the scatterings for all relevant pairs of numerical particles have been computed, the time step is complete and the next one can follow.

To model the small-angle scatters, we apply the drag force and transverse momentum kick to all relevant pairs of numerical particles. Compared to purely forward-dominated scattering (i.e. fSIDM), the cross section entering the drag force is replaced in Eq. (39) as ${\sigma }_{\stackrel{\sim }{\mathrm{T}}}\mapsto {\sigma }_{\stackrel{\sim }{\mathrm{T}},<}({\theta }_{\mathrm{c}})$, see Eq. (12). This takes into account the range of scattering angles treated with the effective ‘small-angle’ approach. The analytical expression of the modified transfer cross section for the hSIDM scheme can be found in Appendix E.

Next, we draw uniform random numbers x, y ∈ [0, 1] for each pair of numerical particles with Λ_ij > 0 to determine whether they undergo ‘large-angle’ scattering (if x ≤ P_i) and, if yes, with which value of the scattering angle (P(θ) = y). Compared to the pure rSIDM implementation, we replace σ_tot ↦ σ_tot, > (θ_c) in Eq. (37) and compute the cumulative probability as⁷

$$\begin{array}{c}\hfill P(\theta )={\displaystyle {\int }_{{\theta }_{\mathrm{c}}}^{\theta }}\mathrm{d}\theta \prime \frac{2\pi sin(\theta \prime )}{{\sigma }_{\mathrm{tot},>}}\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }\prime }\in [0,1]\mathrm{for}\theta \in [{\theta }_{\mathrm{c}},\pi -{\theta }_{\mathrm{c}}].\end{array}$$(41)

As before, the update of particle velocity vectors in the ‘laboratory frame’ is computed analogously as described below Eq. (38). In Appendix D we provide the analytical expression of the cumulative probability for models considered in this work.

We note that within hSIDM both the drag force F_drag = F_drag(θ_c) describing small-angle scattering as well as the probability P_i = P_i(θ_c) and cumulative distribution P(θ) = P(θ; θ_c) for large-angle scattering depend on the critical angle θ_c that is used as a technical parameter to split between small- and large-angle regimes. We check below that all our results are independent of θ_c within an expected range of validity. Finally, we stress that the potentially large (and numerically challenging) cross section σ_tot, < (θ_c) does not enter in hSIDM, but only the much smaller ${\sigma }_{\stackrel{\sim }{\mathrm{T}},<}({\theta }_{\mathrm{c}})$ as well as σ_tot, > (θ_c) and the differential angle-distribution dσ/dΩ for θ_c ≤ θ ≤ π − θ_c.

4.3. Choice of the time step

Compared to conventional N-body simulations taking only gravity into account, the maximal possible size of a single time step should be further constrained for SIDM to ensure accurate results, see for example Vogelsberger et al. (2012) and Fischer et al. (2021b, 2024b). Next, we explain how we choose the time step.

When simulating hSIDM and gravity, the upper limit on the time step related to particle i is determined as

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_i=\min (\mathrm{\Delta }{t}_{\mathrm{grav}.,i},\mathrm{\Delta }{t}_{\mathrm{large}-\mathrm{angle},i},\mathrm{\Delta }{t}_{\mathrm{small}-\mathrm{angle},i}),\end{array}$$(42)

where Δt_grav., i is set by the gravitational acceleration following the criterion given by Springel (2005, see Eq. (34)). The second one (Δt_{large − angle, i}) arises from the requirement that the scattering probability for large-angle scattering within a single time step needs to be sufficiently smaller than unity, P_i ≪ 1⁸, leading to

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{large}-\mathrm{angle},i}\approx \frac{\kappa }{{max}_j\left(|\mathrm{\Delta }{\mathit{v}}_{\mathit{ij}}|\right)}\frac{1}{m_j{\mathrm{\Lambda }}_{\mathit{ii}}}{\left(\frac{{\sigma }_{\mathrm{tot},>}({\theta }_{\mathrm{c}})}{m_{\chi }}\right)}^{-1},\end{array}$$(43)

where κ is a dimensionless accuracy parameter (we use κ = 0.01), m_j is the numerical particle mass, and Λ_ii quantifies the maximum kernel overlap, determined from Λ_ij with j = i. The third criterion (Δt_{large − angle, i}) arises from the effective treatment of frequent small-angle scattering. The description via a drag force and momentum diffusion requires that the typical distribution of deflection angles generated after a single time step is still sufficiently narrow. In practice, this means that the contribution to the (modified) ‘transfer opacity’ (compare to Eq. (26)) from scatterings between t and t + Δt, given by $\mathrm{\Delta }{\tau }_{\stackrel{\sim }{\mathrm{T}},<}({\theta }_{\mathrm{c}})=\rho v\mathrm{\Delta }t{\sigma }_{\stackrel{\sim }{\mathrm{T}},<}({\theta }_{\mathrm{c}})/m_{\chi }$ for a particle moving though a medium with mass density ρ and relative velocity v, needs to be sufficiently small, $\mathrm{\Delta }{\tau }_{\stackrel{\sim }{\mathrm{T}},<}({\theta }_{\mathrm{c}})\ll \mathcal{O}(1)$. The origin of this criterion can be illustrated within the deflection test set-up discussed above, by considering the Molière approximation from Eq. (27), which corresponds to the result obtained within the effective small-angle approach. Here the transfer opacity determines the amount of Gaussian broadening. Hence, within hSIDM the condition $\mathrm{\Delta }{\tau }_{\stackrel{\sim }{\mathrm{T}},<}({\theta }_{\mathrm{c}})\ll \mathcal{O}(1)$ ensures that the contribution to the Gaussian width of the angular distribution from small-angle scatterings, and within a single time step, is sufficiently small (compared to the full range of the azimuthal angle). This leads to

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{small}-\mathrm{angle},i}\approx \frac{\kappa }{{max}_j\left(|\mathrm{\Delta }{\mathbf{v}}_{\mathit{ij}}|\right)}\frac{1}{m_j{\mathrm{\Lambda }}_{\mathit{ii}}}{\left(\frac{{\sigma }_{\stackrel{\sim }{\mathrm{T}},<}({\theta }_{\mathrm{c}})}{m_{\chi }}\right)}^{-1},\end{array}$$(44)

where κ is again an accuracy parameter (we use κ = 0.1). This condition is analogous to the one considered in Fischer et al. (2024b) in the context of purely forward dominated scattering (i.e. fSIDM, for which ${\sigma }_{\stackrel{\sim }{\mathrm{T}},<}({\theta }_{\mathrm{c}})$ is replaced by ${\sigma }_{\stackrel{\sim }{\mathrm{T}}}$), see in particular Sect. 2.3 and Appendix B therein.

As already emphasised, the virtue of the hSIDM scheme is that for models with an enhancement of the differential cross section in the forward direction, the effective treatment of such scatterings is much more efficient as compared to the case when applying the sampling method (rSIDM approach) to the entire range of deflection angles. Indeed, within pure rSIDM, the time step in Eq. (43) would contain σ_tot, instead of σ_tot, > (θ_c) as for hSIDM. Due to σ_tot ≫ σ_tot, > (θ_c), this would lead to the requirement of excessively small time steps in order to be able to treat also small-angle scatterings as individual collisions of numerical particles within a given time step. Within hSIDM, only large-angle scatterings are treated in this way, while small-angle scatterings are treated separately by the effective drag force method. This relaxes the time step requirement by one to two orders of magnitude for typical light mediator models.

We provide a detailed discussion of validation checks of the hSIDM approach in Appendix A. This includes a demonstration of the substantial numerical speed-up by 1–2 orders of magnitude compared to a naive implementation of angle dependence.

5. Simulation of merging galaxy clusters

In this section, we apply the hSIDM method for efficiently treating dark matter self-interaction with angular dependence characteristic for light mediator models. Specifically, we simulate a merger of two galaxy clusters of equal mass, and compare the dynamical evolution of spatial dark matter and galaxy distributions under the influence of gravity as well as dark matter self-interactions with a differential cross section given by Eq. (1). This corresponds to self-interaction via exchange of a light mediator in the Born approximation for indistinguishable particles, to which we refer as ‘Møller’-scattering. In order to demonstrate the impact of the angular dependence, we fix the anisotropy parameter r entering in Eq. (1) for all of our galaxy cluster simulations⁹.

For the remainder of this work, we consider the following configurations:

• (A)

  hSIDM: Angle-dependent self-interaction with a ‘Møller’ cross section (Eq. (1)) for various values of r in the range 0.1 ≤ r ≤ 10⁵, using the hybrid scheme characterised by a dedicated treatment of small- and large-angle scatterings separated by a critical angle θ_c. We verify that all our results are independent of θ_c, by checking agreement between runs with θ_c = 0.1 and θ_c = 0.3.

• (B)

  isotropic rSIDM: Isotropic self-interaction (‘billiard ball scattering’). This commonly considered case can be viewed as arising from taking the ‘heavy mediator’ limit r → 0 in Eq. (1), for which ${\sigma }_{\mathrm{tot}}/{\sigma }_{\stackrel{\sim }{\mathrm{T}}}\to 1$. We often simply refer to this case as rSIDM in this section¹⁰.

• (C)

  fSIDM: Purely forward-dominated scattering (see e.g. Fischer et al. 2021a). Here all dark matter self-interactions are modelled by the effective drag force and momentum diffusion technique. It can be considered as the formal ‘massless mediator’ limit obtained for r → ∞ in Eq. (1) (while keeping the modified transfer cross section ${\sigma }_{\stackrel{\sim }{\mathrm{T}}}$ fixed), for which ${\sigma }_{\mathrm{tot}}/{\sigma }_{\stackrel{\sim }{\mathrm{T}}}\to \infty$.

• (D)

  CDM: For comparison, we also consider the case of collisionless dark matter.

The parameters for each run considered in this work are summarised in Table B.1.

Next, we describe the set-up of the simulations, explain their analysis and demonstrate that the results are insensitive to the choice of the critical angle (Sect. 5.1). It follows a general description of the merger evolution (Sect. 5.2) as well as a discussion how well models of different angular dependencies can be mapped onto each other in the case of merging galaxy clusters (Sect. 5.3). A detailed investigation and discussion of the DM-galaxy offsets follow in Sect. 6.

5.1. Merger simulation set-up, analysis, and validation

Following Fischer et al. (2021b), we consider a merger of two galaxy clusters of equal mass, moving towards each other with an initial relative velocity of 1000 km s⁻¹ along the merger axis (head-on collision), separated by 4000 kpc. The initial dark matter content of each cluster is modelled by a radially symmetric Navarro-Frenk-White (NFW, Navarro et al. 1996) profile, with a concentration parameter c = 5.4 and a scale radius r_s = 3.9 × 10² kpc, corresponding to a virial mass of M_vir = 1 × 10¹⁵ M_⊙. We represent each halo using 10.1 × 10⁶ DM particles with a particle mass of m_DM = 2 × 10⁸ M_⊙. The halo is sampled up to a radius of 7.8 Mpc. We furthermore set the gravitational softening length to ϵ = 1.2 kpc and use N_ngb = 48 neighbouring particles to determine the kernel overlap entering the implementation of dark matter self-interactions, see Sect. 4.

Apart from the DM particles, we include a population of collisionless particles as tracers representing the galaxy distribution, with each halo hosting 10.1 × 10⁶ particles as well, but with a particle mass of m_Gal = 4 × 10⁶ M_⊙. We use an equivalent number of particles for galaxies as for DM to sample their distributions with approximately equal resolution. These particles, which are more abundant than galaxies in clusters, do not represent individual galaxies but rather represent a smoothed-out galaxy distribution. It is noteworthy that while ‘galaxy’ particles in our simulation are approximated to be collisionless, real galaxies may not behave entirely like that (see Kummer et al. 2018). We also include an additional particle with mass m_BCG = 7 × 10¹⁰ M_⊙ positioned at the centre of each halo to simulate the BCG. It is acknowledged that this representation entails an idealised treatment of BCGs, neglecting their spatial extension. The BCG mass is selected to be much lower than typical BCG masses to reduce numerical artefacts. This is in line with previous studies of merging galaxy clusters (e.g. Kim et al. 2017; Fischer et al. 2023; Sabarish et al. 2024); see also Valdarnini (2024) for an alternative approach.

A characteristic signature of dark matter self-interactions are offsets between the dark matter and galaxy populations of the individual galaxy clusters, that are dynamically generated during the merger process. To quantify such offsets, it is necessary to identify a suitable ‘peak’ location of the dark matter and galaxy distributions, respectively, for which various methodologies have been explored (e.g. Power et al. 2003; Kim et al. 2017; Robertson et al. 2017b). Here we adopt the peak identification algorithm proposed by Fischer et al. (2021b), which tracks for each particle whether it belonged originally to the first or the second halo¹¹, and utilises a peak search approach based on the gravitational potential. To compute the peak position for a component, for example, the galaxies of the first halo (association based on the ICs), the gravitational potential sourced by these particles only is computed and the peak position is set by the potential minimum. For the head-on collision set-up studied in this work, the peaks are always located along the merger axis. We denote the cluster moving initially in the negative (positive) direction of the merger axis as ‘halo 1’ (‘halo 2’). The peak positions along the merger axis are denoted by $d_{\mathrm{DM}}^{\mathrm{halo}n}$ and $d_{\mathrm{Galaxies}}^{\mathrm{halo}n}$ for the DM and galaxy populations and the two halos (n = 1, 2), respectively. In addition, we consider the location $d_{\mathrm{BCG}}^{\mathrm{halo}n}$ of the BCG for each cluster. We define the offsets as¹²

$$\begin{array}{c}\hfill d_{\mathrm{DM}}-d_i\equiv \frac{-(d_{\mathrm{halo}1}^{\mathrm{DM}}-d_{\mathrm{halo}1}^i)+d_{\mathrm{halo}2}^{\mathrm{DM}}-d_{\mathrm{halo}2}^i}{2},\end{array}$$(45)

for i = Galaxies or i = BCG, respectively. The sign convention ensures that the offsets for both halos have the same sign. We quote the mean value over the two halos unless stated otherwise.

The hSIDM set-up involves the critical θ_c, that determines the separation between the regimes treated using the dedicated small-angle and large-angle approaches, respectively. As stated above, we checked that all our results are independent of the choice of θ_c within numerical uncertainties. While we provide more details for specific quantities of interest below, we give a first example in Fig. 3, showing the one-dimensional density profile within a small region¹³ around the merger axis, for various times, and for the DM density (solid lines) as well as the galaxy density (dashed lines). Here we used a strongly anisotropic Møller cross section with r = 10⁴. The results obtained from the hSIDM runs for θ_c = 0.1 and θ_c = 0.3 are shown by different colours, and are hardly distinguishable from each other in Fig. 3. This already indicates that peak positions and the resulting offsets between DM and galaxy distributions are independent of θ_c, as we also verify later on.

We also stress that the hSIDM scheme allow us to simulate the strongly forward-dominated model for self-interactions described by r = 10⁴. It would be numerically prohibitively expensive when using an algorithm that does not employ an effective treatment of small-angle scatterings (i.e. for θ_c = 0).

5.2. Peak positions

In Fig. ¹⁴, we present the time evolution of the peak positions of DM and galaxy distributions, as well as the BCG position, along the merger axis. The two sets of lines within each panel correspond to the first and second halo, respectively. For comparison, the upper left panel shows the purely collisionless DM case, that is, without self-interactions. The peaks of DM and galaxies as well as the BCG position are then indistinguishable from each other for all times, as expected. Since we consider an equal-mass merger that occurs head-on, the peaks oscillate along the merger axis (around the centre of mass ) in a symmetrical way, with an amplitude that is damped due to dynamical friction. The other three panels in Fig. show the evolution when assuming various models for DM self-interactions, specifically purely forward-dominated scattering (fSIDM), isotropic scattering (isotropic rSIDM), as well as an anisotropic Møller cross section with r = 10 (hSIDM), respectively. We fixed the absolute self-interaction strength in all these cases such that they correspond to the same viscosity cross section, specifically σ_V/m_χ = 1 cm² g⁻¹. While we discuss the rationale for this choice in more detail below (see Sect. 5.3), we already notice that it leads to an evolution of peak positions that is broadly similar in all cases. Nevertheless, significant differences between the models remain, see Sect. 6.2.

Before going in details, we briefly review some general features for convenience. An initial observation is that, due to the presence of self-interactions, the DM component merges faster, i.e. coalesces on a shorter timescale compared to the galaxy component. Consistent with previous findings by for example Kim et al. (2017) and Fischer et al. (2021a), we observe that while the DM component experiences a rapid merging process, the galactic components undergo relatively stable and prolonged oscillations. These can be attributed to lower central DM densities and thus reduced dynamical friction as compared to the collisionless DM case. Famously, self-interactions also give rise to a spatial separation between DM and galaxy/BCG positions at earlier times, for example shortly after the first pericentre passage at t ≃ 1.8 Gyr or around the first apocentre at t ≃ 2.3 Gyr. The insets in Fig. show the evolution for these times in more detail. Shortly after the first pericentre, namely for t ≲ 2 Gyr, the DM peaks are behind those for galaxies (left inset in each panel), providing a classic signature originating from the deceleration of the DM halos due to self-interaction. Subsequently, the galaxies are decelerated in their outward motion due to the gravitational pull of the DM. When approaching the first apocentre at t ≃ 2.3 Gyr (right inset), the galaxies therefore turn around at a somewhat smaller distance from the centre of mass than the DM peaks. Thus, during this stage, the offset is inverted, and the DM peaks are more outwards directed. A detailed comparison of the right insets for the three models of SIDM also indicate quantitative differences in the relative offsets between DM and galaxies/BCG, see Sect. 6.2 for details.

5.3. Matching of models with different angular dependence

The hSIDM approach allows us to simulate SIDM models with strongly anisotropic differential cross section. Using these results, we are able to scrutinise the question whether certain characteristic features of SIDM can be related among models with distinct differential cross sections by considering a suitable angle-averaged cross sections. Specifically, in this context, we are interested in the question whether the time evolution of DM and galaxy peak positions during the merger process, as well as their offset, can be ‘mapped’ across different models. If such a mapping exists, it would imply that models with different differential cross section but identical angle-averaged cross section yield indistinguishable results for one particular quantity (or ideally even several) of interest.

To investigate this question we consider two well-known candidates to define such a mapping, being the modified transfer cross section ${\sigma }_{\stackrel{\sim }{\mathrm{T}}}$ as well as the viscosity cross section σ_V, see Sect. 2.2. We note that it is clear that the total cross section σ_tot is not a good candidate for identifying a viable ‘mapping’, since it strongly increases for forward-dominated scattering, while the physical effect of self-interactions does not increase correspondingly due to the tiny momentum transfer incurred by each individual scattering event. The viscosity cross section has for example been proposed to provide a ‘mapping’ in the context of gravothermal collapse dynamics of isolated halos in Yang & Yu (2022), see also Sabarish et al. (2024). In the context of galaxy cluster mergers, Fischer et al. (2021a) showed that when matching the modified transfer cross section, offsets for purely forward-dominated scattering are typically larger than for isotropic scattering, indicating that ${\sigma }_{\stackrel{\sim }{\mathrm{T}}}$ may not provide a viable mapping for the DM-galaxy or DM-BCG offset in general. Here we come back to this question, and provide an in-depth investigation.

To be concrete, we considered three specific models, and we simulated all of them with the hSIDM scheme:

• (1)

  Differential Møller cross section Eq. (1) with r = 1, being almost isotropic, with ${\sigma }_{\stackrel{\sim }{\mathrm{T}}}/m_{\chi }\simeq {\sigma }_{\mathrm{V}}/m_{\chi }=1{\mathrm{cm}}^2{\mathrm{g}}^{-1}$.

• (2)

  Differential Møller cross section Eq. (1) with r = 10⁴, being strongly anisotropic, with modified transfer cross section being matched to (1), i.e. ${\sigma }_{\stackrel{\sim }{\mathrm{T}}}/m_{\chi }=1{\mathrm{cm}}^2{\mathrm{g}}^{-1}$. This implies σ_V/m_χ ≃ 1.4 cm² g⁻¹.

• (3)

  As (2), but with viscosity cross section being matched to (1), i.e. σ_V/m_χ = 1 cm² g⁻¹. This implies ${\sigma }_{\stackrel{\sim }{\mathrm{T}}}/m_{\chi }\simeq 0.7{\mathrm{cm}}^2{\mathrm{g}}^{-1}$.

Let us now turn to the DM and galaxy peak positions, as well as their offset. The time evolution is shown in the three columns of Fig. 4, respectively, with each panel including lines for model (1) in black, model (2) in orange, and model (3) in blue. In addition, the lower panels show the difference (2)–(1) (orange) and (3)–(1) (blue). From the left panel in Fig. 4, we see that the DM peak position for model (1) and (2) is similar at all times. This indicates that matching via the transfer cross section is a reasonable approximation for the evolution of the DM peaks. On the other hand, for the galaxy peak position (middle panel in Fig. 4) neither model (2) nor (3) are close to (1) at all times. However, during the most relevant stage from the first pericentre (t ≃ 1.8 Gyr) until around the first apocentre (t ≃ 2.3 Gyr), the galaxy peak position of model (3) is closer to model (1). This suggests that galaxy positions can be better matched using the viscosity cross section, but only until around the first apocentre. Afterwards, no clear matching seems possible.

Finally, we turn to the offset between the DM and galaxy peak positions, shown in the right panel of Fig. 4. We observe that the offsets for both the anisotropic models (2) and (3) are significantly larger than for the isotropic model (1). In particular, this is the case shortly after the first pericentre passage at t ≃ 2 Gyr (when the DM peaks of the two clusters are closer to each other than the galaxy peaks, corresponding to negative offset) as well as around the first apocentre (t ≃ 2.3 Gyr, positive offset). The maximal absolute offsets in the two stages differ by a factor of around two between model (3) and (1), and even by a factor of around four for model (2) and (1). The finding of larger absolute values for the offsets in anisotropic models versus the isotropic case is in line with Fischer et al. (2021a). Thus we find that neither the transfer nor the viscosity cross section yield a useful matching between models when regarding the DM-galaxy offset during the most relevant stages of the merger evolution.

In conclusion, while an approximate mapping of DM peak positions with the transfer cross section is reasonable, the galaxy positions as well as the DM-galaxy offset cannot be mapped with this prescription, with offsets differing by a factor of four between isotropic and strongly anisotropic cases. When comparing models with different angular dependence but identical viscosity cross section, we observe that the galaxy peak position can be roughly mapped until about the first apocentre. However, in that case the DM peak positions differ for the isotropic and anisotropic models, such that the DM-galaxy offsets differ by a factor of around two. Altogether, this implies that the complex merger dynamics requires to include the actual angle dependence for a given model of self-interactions in order to achieve reliable predictions.

6. Model dependence of dark matter-galaxy offsets

In this section we discuss implications of the results obtained from the simulation of galaxy cluster mergers in the previous section. In particular, a main quantity of interest in this context is the offset between the DM halo and either the galaxy population or the brightest galaxy of the cluster (i.e. the BCG). Here we discuss our findings on how these offsets depend on the underlying model describing DM self-interactions. Specifically, our focus is on the characteristic angle dependence of the differential cross section as predicted in realistic (e.g. light mediator) models, that we take into account in the simulations based on the hSIDM scheme developed in this work.

6.1. Dependence on viscosity cross section for fixed anisotropy

As a first step, we verify the expectation that the offsets increase when increasing the overall DM self-interaction rate. To demonstrate this dependence, we simulate galaxy cluster mergers for models with various values of the viscosity cross section in the range 0.2 cm² g⁻¹ ≤ σ_V/m_χ ≤ 3 cm² g⁻¹, but keeping the angle dependence of the cross section fixed (we choose a ‘Møller’ cross section Eq. (1) with r = 10⁴). The time evolution of the DM peak positions for the two halos is shown in Fig. 5 (left panel). Compared to collisionless DM (CDM, black dashed lines), the peaks remain the more close to each other the larger σ_V, as expected. We note that, following Fischer et al. (2021b), we use a linearly re-scaled time variable, t_in, here, defined such that the first (second) pericentre occurs at t_in = 0 (1), respectively. The maximal absolute DM-galaxy and DM-BCG offsets between the first and second pericentre are shown in the right panels of Fig. 5, versus σ_V/m_χ. We observe that the offsets increase roughly linearly within the range of shown viscosity cross sections. To illustrate the magnitude, we extract an approximate estimate for both DM-galaxy and DM-BCG offsets given by

$$\begin{array}{c}\hfill \mathrm{maximal}\mathrm{offset}\approx 80\mathrm{kpc}·\frac{{\sigma }_{\mathrm{V}}/m_{\chi }}{1{\mathrm{cm}}^2{\mathrm{g}}^{-1}}(\mathrm{for}r={10}^4),\end{array}$$(46)

for the equal-mass merger considered in this work.

6.2. Dependence on anisotropy for fixed viscosity cross section

The differential cross section of DM self-scattering within light mediator models is the more forward-dominated the lighter the mediator mass m_ϕ, while it approaches the isotropic limit for a sufficiently heavy mediator. The level of anisotropy is captured by the parameter r defined in Eq. (3), with r ≲ O(1) (r ≫ 1) corresponding to the almost isotropic (strongly forward-dominated) limits, respectively. To assess the dependence of the DM-galaxy offset on the anisotropy, we simulate galaxy cluster mergers as described above using the hSIDM scheme with an angle dependence described by the ‘Møller’ cross section Eq. (1), and for r ∈ {0.1,  1,  10,  10²,  10³,  10⁴,  10⁵}. For comparison, we also consider the purely forward-dominated case (fSIDM) as well as isotropic scattering (isotropic rSIDM). We fix the overall size of the differential cross section such that the angle-averaged viscosity cross section is identical in all cases, with σ_V/m_χ = 1 cm² g⁻¹.

The left panels of Fig. 6 show the time evolution of the DM-galaxy (top) and DM-BCG (bottom) offset, respectively. As expected from Sect. 5, the offsets are negative shortly after the first pericentre passage, since the DM halos were slowed down when passing through each other, while the offsets are positive around the first apocentre, since the galaxy populations are subsequently slowed down by the gravitational force exerted by the DM mass distribution. The overall size of the offsets tend to be the larger the more anisotropic the scattering is, as can be seen by comparing the lines showing the evolution for r = 1, 10, 10², respectively. This is in line with the findings of Fischer et al. (2021a), where offsets for fSIDM were compared to those for isotropic rSIDM¹⁵. Indeed, we find that our results for r = 1, 10, 10² yield offsets that typically lie between those for fSIDM and isotropic rSIDM within the numerical accuracy. To highlight this point, we show the maximal absolute offset within the time interval t ∈ [1.5 Gyr, 3 Gyr] in the right panels of Fig. 6, versus the model parameter r. These results allow us to infer the relevance of model-dependence, being one of the main results to be discussed next.

6.3. Discussion

Models of dark matter self-interactions typically feature angle-dependent scattering. Here we use our results to assess in how far the limits of purely forward-dominated (fSIDM) or isotropic (isotropic rSIDM) scattering yield an accurate approximation to a true angle dependence described by the ‘Møller’ cross section Eq. (1). Its anisotropy is determined by the parameter r (see Eq. (3)), with fSIDM (isotropic rSIDM) being approached in the limiting cases r → ∞ (r → 0), respectively. Here, we discuss for which finite values of r these limits are approached in practice. In particular, we focus on the maximal DM-galaxy and DM-BCG offsets shown in the right column of Fig. 6. We notice several features:

• The maximal DM-galaxy and DM-BCG offsets that we obtain by simulating self-interactions described by the ‘Møller’ cross section Eq. (1) with the hSIDM scheme converge towards those obtained for purely forward-dominated scattering (fSIDM, shown by the light grey horizontal band in Fig. 6) for r ≳ 10². In terms of the model parameters, this corresponds to a mass ratio of mediator and dark matter mass

  $$\begin{array}{c}\hfill m_{\varphi }/m_{\chi }\lesssim 0.1v/c,\end{array}$$(47)

  where v stands for the typical velocity scale of the considered system. Equivalently, the condition for approaching the fSIDM limit can be expressed as ${\sigma }_{\mathrm{tot}}/{\sigma }_{\stackrel{\sim }{\mathrm{T}}}\gtrsim 10$. We note that this ratio is formally assumed to be infinite for fSIDM. The inequality from above can thus be viewed as a criterion for when this limiting case is applicable in practice for realistic models of self-interacting dark matter, and constitutes one of the main results of this work.

• On the other hand, for r ≲ 1 the maximal offsets are compatible with those obtained for isotropic scattering (isotropic rSIDM, shown by the dark grey horizontal band in Fig. 6). This corresponds to

  $$\begin{array}{c}\hfill m_{\varphi }/m_{\chi }\gtrsim v/c,\end{array}$$(48)

  or equivalently ${\sigma }_{\mathrm{tot}}/{\sigma }_{\stackrel{\sim }{\mathrm{T}}}\lesssim 1.1$. We note that the ratio is unity for the limit of exactly isotropic scattering. Thus, our result shows by ‘how much’ a model may quantitatively deviate from isotropy, while still leading to offsets that are compatible with those obtained when assuming isotropic scattering.

• Consequently, for a mediator mass in the range 0.1v/c ≲ m_ϕ/m_χ ≲ v/c, the precise angular dependence of DM self-scattering needs to be taken into account to obtain accurate predictions for DM-galaxy and DM-BCG offsets. This corresponds to $1.1\lesssim {\sigma }_{\mathrm{tot}}/{\sigma }_{\stackrel{\sim }{\mathrm{T}}}\lesssim 10$.

• Finally, coming back to the discussion of validity checks from Appendix A.3 and Sect. 5.1, we note that (like all other considered quantities) the offsets shown in Fig. 6 are independent of the critical angle θ_c entering the hSIDM scheme (within numerical errors), as demonstrated by showing results for θ_c = 0.1 and 0.3, respectively.

Let us stress that the angle dependence enters the merger dynamics in general in a complex way, affecting galaxy and DM distributions differently, see Sect. 5.3. Nevertheless, when considering only the DM-galaxy offset, one may ask how large the error on the inferred self-interaction cross section would be when assuming an incorrect model for the angle dependence. Given that the offsets differ by about a factor of two between the limits of small and large r, see Fig. 6, and that their dependence on the viscosity cross section is approximately linear, see Sect. 6.1, we estimate that the deviation in σ_V/m_χ would amount to about a factor of two as well. For the modified transfer cross section, this number would go up to a factor of approximately three (see Sect. 5.3). However, we stress that these estimates should be taken as indicative only. Fischer et al. (2021b) found that the difference in the size of offsets between isotropic and forward-dominated scattering is larger in unequal-mass mergers compared to equal-mass mergers. Accordingly, this would imply a larger estimate for the error than the factor of two for the viscosity cross section in the case of merging galaxy clusters with a significantly unequal mass ratio. Modelling an observed galaxy cluster merger with detailed simulations of various models for dark matter self-interactions and astrophysical settings could reveal more information than when just considering (maximal) offsets based on a peak identification method that is not applicable to observations.

We note that our study considers head-on equal mass mergers only and does not model the intra-cluster medium. It would be interesting to extend this analysis to unequal-mass merges and allow for a non-central collision. This as well as modelling the intra-cluster medium and including the velocity-dependence of the cross section is left to future work. Nevertheless, our results can give a valuable hint for whether a given model for DM self-interactions can be approximately described by the purely forward-dominated fSIDM or the isotropic rSIDM limits, respectively.

7. Conclusions

In this work, we have proposed, validated, and employed a numerical scheme to efficiently perform N-body simulations of SIDM while taking into account the angle dependence predicted in realistic particle physics models. Most previous works considered either the case of relatively rare scatterings by large angles (‘rSIDM’, typically assuming an isotropic distribution) or very frequent and purely forward-dominated scatterings (‘fSIDM’). The hybrid scheme developed in this work, dubbed ‘hSIDM’, combines the methods underlying both of these approaches, yielding an efficient algorithm for incorporating DM self-interactions described by an (in principle) arbitrary differential cross section dσ/dΩ.

As a first step, we extensively validated the approach based on an analytically solvable test set-up. In this context, we also derived a quantitative validity criterion for the effective treatment of small-angle scattering by a drag force combined with a diffusive transverse motion. The insights were used to show under which conditions the hSIDM scheme is expected to improve over a naive sampling of scattering angles in each collision of N-body particles, and we demonstrated a speed-up by several orders of magnitude in computation time. The theoretically expected speed-up was estimated in terms of a ratio of angle-weighted cross sections in Eq. (16). We checked the independence of our results of the auxiliary ‘critical angle’ entering the hSIDM scheme within a validity regime derived from the condition given in Eq. (36).

Next, we applied hSIDM to simulate equal-mass galaxy cluster mergers including DM self-interactions with the characteristic angle dependence predicted by light mediator models. Concretely, we assumed a dependence of the scattering probability on the deflection angle given by Eq. (1). In this parameterisation, the model parameter r (related to the ratio of mediator and DM mass; see Eq. (3)) controls the amount of anisotropy, with r ≫ 1 corresponding to strongly forward-dominated scattering, and r ≪ 1 corresponds to the isotropic limit. We first investigated whether simulations with a different r but an identical transfer or viscosity cross section can be mapped onto each other, and we found that this is in general not the case, especially when jointly considering the time evolution of the DM distribution, the galaxy population, the position of the BCG, and their relative offsets. Nevertheless, the dependence of the offset on r is somewhat weaker when fixing the viscosity cross section as compared to the case when keeping the (modified) transfer cross section fixed.

Finally, we investigated the dependence of the DM-galaxy and DM-BCG offsets on the properties of the model for DM self-interactions. We find that the maximal offset approaches the (larger) value predicted by purely forward-dominated scattering for r ≳ 10². This corresponds to a mediator to DM mass ratio of m_ϕ/m_χ ≲ 0.1v/c, where v is the typical velocity scale. Equivalently, this implies ${\sigma }_{\mathrm{tot}}/{\sigma }_{\stackrel{\sim }{\mathrm{T}}}>10$ for the ratio of the total and modified momentum transfer cross sections. Conversely, the maximal offset converges to the (smaller) value obtained for isotropic self-scattering when r ≲ 1 (i.e. m_ϕ/m_χ ≳ v/c or ${\sigma }_{\mathrm{tot}}/{\sigma }_{\stackrel{\sim }{\mathrm{T}}}\lesssim 1.1$). Our results thus yield a quantitative criterion for when the fSIDM and rSIDM limits are applicable for extracting offsets. Beyond that, we find that the angle dependence of DM self-scattering influences the time evolution of DM and galaxy populations in a distinct way, suggesting that it should be taken into account when modelling observed galaxy cluster mergers with simulations of self-interacting DM.

Our work can be extended by considering broader classes of galaxy cluster merger set-ups (unequal mass, non-central) with a more realistic modelling (e.g. including the intracluster medium) and by also taking the interplay of the angle- and velocity-dependence of the differential cross section into account. Furthermore, the hSIDM approach could be applied to other situations, such as the evolution of satellite populations, for a precision study of the impact of self-interactions with pronounced angle dependence on isolated halos or for studying the approach of gravothermal collapse.

Movie

Movie 1 associated with Fig. 3 Access here

Data availability

Movie associated to Fig. 3 is available at https://www.aanda.org

---

Acknowledgments

The authors thank Tobias Binder, Felix Kahlhoefer, Sowmia Balan, Marc Wiertel and all participants of the Darkium SIDM Journal Club for discussion. CA thanks Christoph Selbmann for helpful discussion. This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 ‘Origins’ – 390783311. Software: NumPy (Harris et al. 2020), Matplotlib (Hunter 2007).

References

Appendix A: Validation based on deflection test set-up

In this section we validate the hSIDM implementation described in Sect. 4. Moreover, we demonstrate the numerical gain in efficiency of this scheme as compared to a naive sampling of deflection angles.

A.1. Simulation set-up and basic validation

We consider the deflection test set-up discussed in Sect. 3.2, consisting of a beam of test particles moving through a target. So we compare our numerical results to the analytical solution Eq. (18) for the time evolution of the distribution function f(t, θ) of azimuthal angles of beam particles relative to the initial beam axis.

We set up an initial configuration with 8 000 test particles, initially moving with velocity v = 2 km s⁻¹ in a common beam axis direction, through a medium described by 92 000 target particles, with total mass 1 × 10¹⁰ M_⊙ in a cube with side length of 14 kpc, corresponding to a background density ρ = 3.353 × 10⁶ M_⊙ kpc⁻³. Here we use N_ngb = 64 neighbouring particles. Beam particles can scatter off target particles with a given differential cross section dσ/dΩ. We did not allow for scattering of the beam with beam or target with target particles, and we switched off gravity for the deflection test set-up. However, each beam particle can undergo multiple scatterings, with in principle any of the target particles.

As mentioned already in Sect. 3.2, we consider two configurations:

• (i)

  no recoil is transmitted from beam to target particles. Physically this can be realised if beam and target particles are distinguishable, and target particles are formally considered to have infinite inertia. In this case the magnitude of the initial velocity is conserved for beam particles. This case is valuable due to the known exact analytical solution Eq. (18) of the angular distribution function.

• (ii)

  recoil is transmitted from beam to target particles, assuming they have equal inertial mass. This is the physically more relevant case, realised if beam and target particles are identical particle species.

Furthermore, for the test set-up, we consider the three following methods:

• (a)

  hSIDM scheme: The hybrid approach introduced above, for various choices of the critical angle θ_c to divide between small and large-angle scattering regimes.

• (b)

  anisotropic rSIDM scheme: The case for which all scatterings are sampled explicitly from the differential cross section. We often refer to this case simply as rSIDM in this section. It corresponds to the (computationally expensive) limit of hSIDM for θ_c → 0.¹⁶

• (c)

  fSIDM scheme: The case for which all scatterings are treated according to the effective drag force and momentum diffusion approach. We refer to this case as fSIDM. It corresponds to hSIDM for θ_c → π/2 (π) for indistinguishable (distinguishable) particles.

We furthermore consider the various examples for differential cross sections as given in Sect. 2.1.

Since the parameters of the test set-up are chosen somewhat arbitrarily, we normalise all timescales to a given reference time t₁, defined to be the dynamical time t_dyn on which the angular distribution is expected to evolve significantly, see Sect. 3.2. Concretely, we define t₁ such that the opacity τ(t)≡ρvtσ/m_χ satisfies τ(t₁) = 1 for a generic cross section σ/m_χ = 1 cm² g⁻¹. We note that this definition is insensitive to which type of cross section σ stands for.

As a basic check of the implementation, we also consider an artificial case where each beam particle is allowed to scatter only once during the entire simulation. In this case the angular distribution simply follows the differential cross section, see Appendix G. After this check, we allowed each beam particle to scatter multiple times, as determined according to the chosen scheme.

As a first non-trivial validity check, we consider the no-recoil case (i),¹⁷ for which the exact analytical solution is known, see Eq. (18). In Fig. A.1 we show the angular distribution function f(t, θ) of beam particles at three different times t (solid lines) obtained numerically within the (anisotropic) rSIDM scheme (b), and compare them to the exact analytical solution (dashed), observing good agreement.

In addition, we show the angular distribution obtained from the fSIDM scheme (c) within the no-recoil case (i) for three different times in Fig. A.2 (solid lines). In this case, we compare to the Gaussian Molière distribution from Eq. (27) (dashed lines). Agreement between fSIDM and Eq. (27) is expected as long as the width of the Gaussian distribution is much smaller than π, which is the case for the times t shown in Fig. A.2. As an additional check, the two panels in Fig. A.2 show results for two values of the transfer cross section that differ by a factor of ten. Due to the absence of gravitational forces in the test set-up, the solutions for these two cases are expected to match when rescaling the simulation time by a factor of one-tenth. We find that this is indeed the case, serving as an additional validation of the implementation of the time-step criterion. We note that this check is similar to the one performed by Fischer et al. (2021a), except that we consider the set-up without recoil (case (i)) here. This check validates that the effective treatment of small-angle scattering via a drag force and momentum diffusion leads to the expected Gaussian beam broadening. However, it does not provide a check of whether this treatment is a valid approximation for a given differential cross section and choice of critical angle. We turn to these checks next.

A.2. Validity check of effective small-angle approach

The effective treatment of small-angle scattering via a drag force and momentum diffusion is valid provided the scattering is sufficiently forward dominated, and sufficiently frequent. In this sub-section we provide a quantitative check of these intuitive validity requirements. In order to expose the validity criteria most clearly, we consider the artificial differential cross section introduced in Eq. (4), for which all individual scatterings occur with a fixed scattering angle θ₀ in the centre of mass frame. Details of this fixed scattering angle case are described in Appendix I. We solve the deflection test set-up for various fixed values of θ₀. Since all scatterings occur with a single angle, we expect the effective small-angle approach to work the better the smaller θ₀.

In order to check this expectation, we run two simulations for each value of θ₀: one using the effective small-angle approach (fSIDM, scheme (b)), and one using the explicit sampling technique (rSIDM, scheme (c)). The latter is numerically more expensive, but yields the exact angular distribution function f(t, θ) at all times t (within the numerical accuracy). This allowed us to check whether the effective small-angle approach is valid by comparing the angular distribution obtained from the fSIDM simulation to the rSIDM result.

We find that both angular distributions agree with each other (within numerical errors) after a finite simulation time t_min(θ₀) which depends on the choice of θ₀. By agreement we mean, that the deviation between the distribution functions is strictly smaller than a certain threshold for all angles. Numerically, we carry this out by subtracting both angular distributions and normalising it to the square root of the number of particles in each bin. In our case, we use an error threshold of 0.5. This minimal time is shown by the cross symbols in Fig. A.3 (relative to the common reference time t₁ introduced above), for three different values of the transfer cross section, and for the case (i) without recoil as well as the case (ii) with recoil, respectively.

This finding is in line with the analytical discussion in Sect. 3.2.2, based on the analysis of the exact analytical solution of the deflection problem in the no-recoil case. Indeed, our numerical results precisely follow the expected behaviour; that is, the fSIDM approach becomes valid only after some minimal time that is required to lead to sufficiently many scatterings for the effective drag force description to apply. We find that the dependence of the minimal time t_min derived analytically in Eq. (32) on θ₀ as well as on the overall size of the cross section (shown by solid lines in Fig. A.3) agrees with our results for the numerically determined minimal time t_min(θ₀) (cross symbols).¹⁸ This demonstrates that the validity criteria derived analytically in Sect. 3.2.2 are applicable in practice. Moreover, our numerical results for the deflection test cases with and without recoil suggest that the validity criteria derived for the no-recoil set-up can be carried over to the case with recoil. This is an important observation, since all our later results are based on scatterings among identical particles, including recoil.

As can be observed in Fig. A.3, the minimal time after which the fSIDM approach becomes applicable is the shorter the smaller θ₀. This agrees with the intuitive notion that the fSIDM approach works the better the smaller the (typical) scattering angles are in individual collisions (which are given by θ₀ for the differential cross section considered in this sub-section, and in general bounded from above by θ_c within the hSIDM approach). As discussed in Sect. 3.2.2, the fSIDM approach is useful if t_min ≪ t_dyn, where the latter is the dynamical timescale over which the angular distribution evolves significantly. For the deflection set-up, t_dyn = (1, 0.1, 0.01)⋅t₁ for σ_T/m_χ = (1, 10, 100) cm² g⁻¹, respectively. We find

$$\begin{array}{c}\hfill t_{\min }/t_{\mathrm{dyn}}\simeq {10}^{-2}·{\left(\frac{{\theta }_0}{0.1}\right)}^2,\end{array}$$(A.1)

as long as θ₀ ≲ 1, independently of σ_T/m_χ. The dependence on θ₀² as well as the independence of the ratio t_min/t_dyn on the overall size of the cross section is consistent with the analytical result Eq. (32), as stated above.

Altogether, this implies scatterings by angles θ ≲ 1 can be accurately treated within the effective small-angle (fSIDM) approach. This is valid on timescales that are one to two orders of magnitude shorter than the typical dynamical scale of the system.

A.3. Validity check of hSIDM approach

Let us now return to the hSIDM approach for efficiently describing angle-dependent DM self-scattering. We solve the deflection test problem using the hSIDM scheme (a) employing various critical angles θ_c, see Sect. A.1. We consider first the case (i) without recoil, for which we can compare the angular distribution f(t, θ) to the analytical solution Eq. (18). We assume a ‘Rutherford’ cross section given in Eq. (2), with three values of the anisotropy parameter (related to the mediator-to-dark matter mass ratio) r = 10², 10³, 10⁴.

The angular distribution of beam particles at three different simulation times is shown in Fig. A.4 for r = 10⁴. As expected, the hSIDM approach yields an angular distribution that agrees with the analytical solution only after some minimal time t > t_min(θ_c). The reason is that scatterings by small angles (θ < θ_c) are treated with the effective drag force and momentum diffusion method within hSIDM, and this approach requires a minimal time after which scatterings are sufficiently frequent to apply, as discussed in detail above. We observe that for the smaller critical angle θ_c = 0.3 (left panel in Fig. A.4) the agreement is reached at earlier times as compared to θ_c = 0.5 (right panel), i.e. t_min(0.3) < t_min(0.5). This is also expected, since in the former case only scatterings by smaller angles are included in the effective small-angle approach, such that it becomes valid earlier on.

We note that t_min(θ_c)≪t_dyn(=t₁) is much smaller than the dynamical timescale for both values of θ_c shown in Fig. A.4. This means the hSIDM approach yields a valid description on the relevant timescale of order t_dyn. In the upper panel of Fig. A.5 we show t_min(θ_c) as obtained from hSIDM simulations with various critical angles θ_c = 0.1, 0.3, 0.5, 0.7, 0.9, and three values of r (cross symbols). In addition, the dependence of the analytical prediction for t_min from Eq. (35) on θ_c and on r is shown (solid lines). We observe that the analytical result yields a useful estimate for the dependence of the minimal validity timescale on the model parameters and on the critical angle. We note that for very small values of θ_c < 0.1, the analytical result for t_min(θ_c) strongly increases. This is due to the fact that for such small values of the critical angle, only a very small range of scattering angles is treated via the drag force method, such that it formally leads to sufficiently frequent scatterings only after a rather long time. However, we notice that for θ_c = 0.1 the value of t_min extracted from our numerical solutions is actually still rather small (left-most crosses in Fig. A.5). We attribute this to the fact that the physical impact of scatterings with θ > θ_c is dominant for such small values of θ_c. Since those are treated based on the explicit sampling technique within hSIDM, the potentially larger relative error of the small-angle regime is not noticeable in the complete angular distribution function for θ_c = 0.1.

Next, we considered the deflection test set-up with recoil, see case (ii) described in Sect. A.1. We use a ‘Møller’ cross section given in Eq. (1) (as appropriate for beam and target particles with identical mass), again with three values r = 10², 10³, 10⁴ and for θ_c = 0.1, 0.3, 0.5, 0.7, 0.9, respectively. Since no analytical result is available in this case, we instead compare the angular distribution when using the hSIDM scheme (a) to the one obtained when using the sampling technique for scatterings by all angles (anisotropic rSIDM scheme (b)), see Sect. A.1. Similarly as before, we find that the angular distributions f(t, θ) obtained from hSIDM and rSIDM agree after some minimal time t_min(θ_c), shown by the plus symbols in the lower panel of Fig. A.5. The dependence of t_min(θ_c) on θ_c and on r is similar as for the no-recoil case (i). For comparison, we also show the dependence of the analytical prediction for t_min(θ_c) from Eq. (35) on θ_c and on r by the solid lines in the lower panel of Fig. A.5. We stress that this analytical result has been derived within the no-recoil case. Nevertheless, we find that it yields a useful analytical estimate of t_min(θ_c) also in the physically relevant case with recoil, similarly as observed before for the set-up considered in Sect. A.2.

Thus, we find that the validity criteria for applying the effective small-angle approach derived in Sect. 3.2.3 are applicable in practice also beyond the strict set-up for which they have been derived, and yield a useful quantitative estimate for when the hSIDM approach can be used. Finally, we notice that t_min(θ_c) is significantly smaller than the dynamical timescale t_dyn(≡t₁) as long as θ_c is sufficiently small for both cases with (i) and without (ii) recoil. In addition, the ratio t_min/t_dyn is the smaller the larger r. This implies the hSIDM approach is valid particularly for very strongly forward-dominated differential cross sections. Since the gain in numerical efficiency is also most pronounced for large r, this implies that hSIDM is indeed applicable in the cases for which it is expected to be most advantageous from the point of view of computational resources. We note that as an additional check, we explicitly verified that all results presented in Sect. 5 are insensitive to the choice of critical angle θ_c, being particularly relevant for moderate values of the anisotropy parameter r ≲ 𝒪(10²).

For the explanations given above, we assumed that the dynamical time scale is set by the self-interactions. This is not true in general, for example, the evolution of astrophysical systems typically involves gravity setting its own time scale. There is no positive value for θ_c that can always be safely used. Instead, accurately modelling the evolution of a system requires to choose the critical angle in respect to its dynamical time scale (i.e. t_min(θ_c)≪t_dyn). We come back to this point in case of galaxy cluster mergers in Sect. 5.

A.4. Performance gain due to the hSIDM approach

The hSIDM approach to angle-dependent SIDM allows us to simulate differential cross sections that are characteristic for realistic (e.g. light mediator) models in a numerically significantly more efficient way as compared to a naive approach for which scatterings by all angles are sampled from dσ/dΩ. Instead, the division into small- and large-angle scattering regimes and the treatment of the former by an effective drag force and momentum diffusion approach within hSIDM speeds up the numerical performance. In this section we quantify the numerical efficiency. We consider again the deflection test set-up, using the variant with recoil (case (ii) in Sect. A.1) and assuming a ‘Møller’ cross section, see Eq. (1), with r = 10⁴. We perform runs using the hSIDM approach, with various values for the critical angle θ_c used to separate small and large-angle regimes. We compare to the naive approach for which all scattering angles are treated with the sampling technique (anisotropic rSIDM), which corresponds to hSIDM in the limit θ_c → 0.

The required CPU run-time t_cpu is shown in the upper panel of Fig. A.6, versus the simulation time t. We observe that the CPU run-time significantly decreases when increasing θ_c. We note that the physical results for the distribution of beam particles are compatible for all shown cases of θ_c, i.e. hSIDM yields a valid description, but at a numerical cost that is almost a factor 100 lower for θ_c ∼ 0.3. This is comparable to (but slightly better than) the theoretically expected speed-up, see Eq. (16) and Fig. 2.

The lower panel of Fig. A.6 shows the required CPU run-time Δt_cpu per simulation time-interval Δt_sim, which is roughly t-independent. We normalise the y-values to unity for the anisotropic rSIDM scheme, and show the performance gain within hSIDM when increasing the value of the critical angle θ_c. The improvement is most drastic for values θ_c ∼ 0.1, and saturates for θ_c ≳ 0.3. The saturation indicates that the numerical cost of the scatterings treated via the effective small-angle approach becomes noticeable for these values.

In summary, we verified the applicability and accuracy of the hSIDM approach. In particular, we demonstrated that it yields an improvement in computational performance by up to a factor 100 for models of SIDM featuring a differential cross section with enhancement at small scattering angles, as is typical for light mediator models.

Appendix B: Table of simulated self-interacting dark matter models and parameters

In this Appendix, we present and summarise the parameters for each run considered in this work (see Sect. 5) in Tab. B.1.

Appendix C: Figure of time evolution of peak positions in merging galaxy clusters

In this Appendix we present the time evolution of the peak positions of DM and galaxy distributions, as well as the BCG position, along the merger axis in Fig. C.1.

Appendix D: Implementation of anisotropic large-angle scattering for hSIDM

The contribution of large-angle scattering within the hSIDM scheme is implemented by explicitly sampling deflection angles from the differential cross section, as discussed in Sect. 4.2, following the method used in Robertson et al. (2017a). However, there is a small difference. Within the hSIDM scheme, this requires the cumulative probability function for scattering between an angle θ_c and θ given in Eq. (41). We find that it can be computed analytically for the models considered in this work. For the ‘Møller’ scattering cross section Eq. (1) the cumulative distribution function reads

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P(\theta )& =\frac{{\sigma }_0}{4{\sigma }_{\mathrm{tot},>}({\theta }_{\mathrm{c}})}\frac{a(\theta )}{b(\theta )},\mathrm{with}\hfill \\ \hfill a(\theta )& =4r(2+r)c_2({\theta }_{\mathrm{c}})cos(\theta )\hfill \\ \hfill & -c_3(\theta )(8r(2+r)cos({\theta }_{\mathrm{c}})+c_2({\theta }_{\mathrm{c}})log\left(c_1(\theta )\right)),\hfill \\ \hfill b(\theta )& =r(2+r)c_3(\theta )c_3({\theta }_{\mathrm{c}}),\hfill \\ \hfill c_1(\theta )& =\frac{(2+r-rcos(\theta ))(2+r+rcos({\theta }_{\mathrm{c}}))}{(2+r+rcos(\theta ))(2+r-rcos({\theta }_{\mathrm{c}}))},\hfill \\ \hfill c_2({\theta }_{\mathrm{c}})& =-8-r(8+r)+r^{2}cos(2{\theta }_{\mathrm{c}}),\hfill \\ \hfill c_3(\theta )& =-{(2+r)}^2+r^2{cos}^2(\theta ),\hfill \end{array}\end{array}$$(D.1)

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{tot},>}({\theta }_{\mathrm{c}})={\sigma }_0(\frac{4cos({\theta }_{\mathrm{c}})}{{(2+r)}^2-r^2{cos}^2({\theta }_{\mathrm{c}})}+\frac{log(-1+\frac{2(2+r)}{2+r+rcos({\theta }_{\mathrm{c}})})}{r(2+r)}).\end{array}$$(D.2)

For the ‘Rutherford’ scattering cross section Eq. (2), the cumulative distribution function is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P(\theta )& =\frac{2{\sigma }_0}{{\sigma }_{\mathrm{tot},>}({\theta }_{\mathrm{c}})}\frac{(cos({\theta }_{\mathrm{c}})-cos(\theta ))}{(2+r-rcos({\theta }_{\mathrm{c}}))(2+r-rcos(\theta ))},\hfill \\ \hfill \end{array}\end{array}$$(D.3)

with

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{tot},>}({\theta }_{\mathrm{c}})={\sigma }_0\frac{(1+cos({\theta }_{\mathrm{c}}))}{(1+r)(2+r-rcos({\theta }_{\mathrm{c}}))}.\end{array}$$(D.4)

Appendix E: Implementation of frequent small-angle scattering for hSIDM

The contribution from small-angle scattering within the hSIDM scheme is implemented using an effective description based on a drag force as well as transverse momentum diffusion, following Fischer et al. (2021a). The drag force is determined by the (modified) transfer cross section. For the hSIDM scheme, we use the contribution to this quantity from small-angle scatterings only, see Sect. 4.2.

For the ‘Møller’ scattering cross section Eq. (1) the modified transfer cross section for scattering by angles θ ≤ θ_c or θ ≥ π − θ_c reads

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\sigma }_{\stackrel{\sim }{\mathrm{T}},<}({\theta }_{\mathrm{c}})& =2{\displaystyle \underset{\theta \le {\theta }_{\mathrm{c}}\vee \theta \ge \pi -{\theta }_{\mathrm{c}}}{\int }}\mathrm{d}\mathrm{\Omega }(1-|cos(\theta )|)\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}\hfill \\ \hfill & =2{\sigma }_0\frac{a({\theta }_{\mathrm{c}})}{b({\theta }_{\mathrm{c}})},\hfill \end{array}\end{array}$$(E.1)

where the numerator is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill a(\theta )=& 24log(c_1(\theta ))\hfill \\ \hfill & -2r[-4log(c_2(\theta ))+r[2(2+r)cos(\theta )\hfill \\ \hfill & -(7+r)log(2+r-rcos(\theta ))-(11+2r)log(c_3(\theta ))\hfill \\ \hfill & +{cos}^2(\theta )(-4+3log(c_4(\theta ))\hfill \\ \hfill & +r(-2+log(c_5(\theta )))-6log(2))+3(6+r)log(2)]],\hfill \end{array}\end{array}$$(E.2)

and the denominator reads

$$\begin{array}{c}\hfill b(\theta )=r^2(2+r)(2+r-rcos(\theta ))(2+r+rcos(\theta )).\end{array}$$(E.3)

Here we defined several terms:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill c_1(\theta )& =\frac{(2+r-rcos(\theta ))(2+r+rcos(\theta ))}{4(1+r)},\hfill \\ \hfill c_2(\theta )& =\frac{{(2+r-rcos(\theta ))}^4{(2+r+rcos(\theta ))}^5}{512{(1+r)}^5},\hfill \\ \hfill c_3(\theta )& =\frac{(2+r+rcos(\theta ))}{(1+r)},\hfill \\ \hfill c_4(\theta )& =\frac{(2+r-rcos(\theta ))(2+r+rcos(\theta ))}{(1+r)},\hfill \\ \hfill c_5(\theta )& =\frac{(2+r-rcos(\theta )){(2+r+rcos(\theta ))}^2}{8{(1+r)}^2}.\hfill \end{array}\end{array}$$(E.4)

The ‘Rutherford’ scattering cross section Eq. (2) is used for cases with distinguishable particles, for which the transfer cross section is relevant. In addition, only scatterings by θ ≤ θ_c are included in the small-angle regime. We thus consider

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\sigma }_{\mathrm{T},<}({\theta }_{\mathrm{c}})& ={\displaystyle \underset{\theta \le {\theta }_{\mathrm{c}}}{\int }}\mathrm{d}\mathrm{\Omega }(1-cos(\theta ))\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}\hfill \\ \hfill & ={\sigma }_0\frac{c_1({\theta }_{\mathrm{c}})+2log\left(c_2({\theta }_{\mathrm{c}})\right)-2}{r^2},\hfill \end{array}\end{array}$$(E.5)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill c_1({\theta }_{\mathrm{c}})& =\frac{4}{2+r-rcos({\theta }_{\mathrm{c}})},\hfill \\ \hfill c_2({\theta }_{\mathrm{c}})& =\frac{(2+r-rcos({\theta }_{\mathrm{c}}))}{2}.\hfill \end{array}\end{array}$$(E.6)

Appendix F: Vector calculation for anisotropic scattering

For large-angle anisotropic scattering, the deflection angle is defined in the centre of mass frame of the two scattering particles. In this section, we describe an algorithm how to obtain the velocity vectors after the scattering. We consider the most general case where both particles can have different masses (relevant for distinguishable particles), denoted by m₁ and m₂, respectively.

Given two arbitrary initial velocity vectors of particles that scatter off each other u₁, u₂ ∈ ℝ³, we first transform to the rest frame of the second particle using (see illustration in Fig. F.1)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathit{V}}_1& \equiv {\mathit{u}}_1-{\mathit{u}}_2,\hfill \\ \hfill {\mathit{V}}_2& \equiv {\mathit{u}}_2-{\mathit{u}}_2=0.\hfill \end{array}\end{array}$$(F.1)

Next, we consider the centre of mass frame, for which

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathit{v}}_1={\mathit{u}}_1-{\mathit{v}}_{\mathrm{com}}=\frac{m_2}{m_1+m_2}({\mathit{u}}_1-{\mathit{u}}_2),\hfill \\ \hfill & {\mathit{v}}_2={\mathit{u}}_2-{\mathit{v}}_{\mathrm{com}}=-\frac{m_1}{m_1+m_2}({\mathit{u}}_1-{\mathit{u}}_2),\hfill \\ \hfill \end{array}\end{array}$$(F.2)

where

$$\begin{array}{c}\hfill \begin{array}{cc}& {\mathit{v}}_{\mathrm{com}}\equiv \frac{m_1{\mathit{u}}_1+m_2{\mathit{u}}_2}{m_1+m_2}.\hfill \end{array}\end{array}$$(F.3)

To simplify the calculation, we consider in addition a three-dimensional rotation, defined by the orthogonal matrix

$$\begin{array}{c}\hfill R=\left(\begin{array}{ccc}cos(\alpha )& -sin(\alpha )cos(\beta )& sin(\alpha )sin(\beta )\\ sin(\alpha )& cos(\alpha )cos(\beta )& -cos(\alpha )sin(\beta )\\ 0& sin(\beta )& cos(\beta )\end{array}\right),\end{array}$$(F.4)

with

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \beta & =arctan(-\frac{\mathrm{z}}{\mathrm{y}}),\hfill \\ \hfill \alpha & =arctan(-\frac{cos(\beta )\mathrm{y}-sin(\beta )z}{\mathrm{x}}).\hfill \end{array}\end{array}$$(F.5)

In this rotated frame, the initial velocities in the centre of mass frame are aligned along the x-axis. The velocity vectors before (un-primed) and after (primed) scattering read as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathit{v}}_1& =v_1{\widehat{e}}_{\mathrm{x}}\hfill \\ \hfill {\mathit{v}}_2& =-\frac{m_1}{m_2}{v}_1{\widehat{e}}_{\mathrm{x}}\hfill \\ \hfill {\mathit{v}}_1^{\prime }& =v_1(cos(\theta ){\widehat{e}}_{\mathrm{x}}+sin(\theta ){\widehat{e}}_{\mathrm{y}})\hfill \\ \hfill {\mathit{v}}_2^{\prime }& =\frac{m_1}{m_2}{v}_1(-cos(\theta ){\widehat{e}}_{\mathrm{x}}-sin(\theta ){\widehat{e}}_{\mathrm{y}})\hfill \\ \hfill {\mathit{V}}_1& =(1+\frac{m_1}{m_2})v_1{\widehat{e}}_{\mathrm{x}}\hfill \\ \hfill {\mathit{V}}_1^{\prime }& =v_1((cos(\theta )+\frac{m_1}{m_2}){\widehat{e}}_{\mathrm{x}}+sin(\theta ){\widehat{e}}_{\mathrm{y}})\hfill \\ \hfill {\mathit{V}}_2^{\prime }& =\frac{m_1}{m_2}{v}_1((1-cos(\theta )){\widehat{e}}_{\mathrm{x}}-sin(\theta ){\widehat{e}}_{\mathrm{y}}),\hfill \\ \hfill \end{array}\end{array}$$(F.6)

where θ is the scattering angle in the centre of mass frame, which in practice is chosen according to the cumulative distribution function in Appendix D. We note that the angles ψ and ζ indicated in Fig. F.1 are related to θ by

$$\begin{array}{c}\hfill tan(\psi )=\frac{sin(\theta )}{cos(\theta )+\frac{m_1}{m_2}},\end{array}$$(F.7)

and ζ = (π − θ)/2. For the choice from above the velocity vectors lie in the xy-plane after scattering. To capture the most general case, we rotate them around the x-axis by an arbitrary rotation angle φ ∈ [0, 2π] using

$$\begin{array}{c}\hfill R_{\mathrm{x}}=\left(\begin{array}{ccc}1& 0& 0\\ 0& cos(\phi )& -sin(\phi )\\ 0& sin(\phi )& cos(\phi )\end{array}\right).\end{array}$$(F.8)

Now we rotate the vectors back to the initial frame with the inverse matrix R⁻¹ = R^T and finally add back the initial velocity u₂ of the second particle, giving

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathit{u}}_1^{\prime }& :={\mathit{V}}_1^{\prime }+{\mathit{u}}_2,\hfill \\ \hfill {\mathit{u}}_2^{\prime }& :={\mathit{V}}_2^{\prime }+{\mathit{u}}_2.\hfill \end{array}\end{array}$$(F.9)

Appendix G: Single scattering test

As an initial sanity check of the implementation of anisotropic (large-angle) scattering, we consider the case in which each beam particle within the deflection set-up discussed in Sect. A is allowed to scatter only exactly once during the entire duration of the simulation. In this artificial case the angular distribution f(t, θ) has to approach the underlying differential cross section dσ/dΩ for t → ∞. We performed this test for various choices of the differential cross section, and setting θ_c = 0 to capture the entire range of scattering angles for this test purpose. The result is shown in Fig. G.1. We find good agreement between the simulation and the expected shape in all cases.

Appendix H: No-recoil case implementation

For some of the validation tests, and in particular to be able to compare to the exact analytical Goudsmit–Saunderson solution Eq. (18) of the deflection tests problem, we also implement scattering of distinguishable types of particles. Specifically, we consider the case in which ‘beam’ particles have some (finite) mass m_χ, while ‘target’ particles have a much larger mass M → ∞.¹⁹ In this limit, the recoil transmitted to the target vanishes, i.e. target particles remain at rest even after scattering.

Here we describe how we adapt the implementation of both the small- and large-angle scattering regimes to this situation. For large-angle scattering, we follow the steps described in Appendices D and F, setting m₁ ≡ m_χ and taking the limit m₂ ≡ M → ∞. We note that the centre of mass frame and the frame where the heavy test particle is at rest formally coincide in this limit, and that the magnitude of the velocity of the beam particle is unchanged before and after the scattering, such that energy is conserved.

For small-angle scattering, we infer the drag force by considering the deceleration rate of an individual physical DM particle travelling with velocity v₀ through a background density ρ_j, and undergoing small-angle scatterings. This is the same ansatz as by Kahlhoefer et al. (2014) in their Appendix A for the case of equal masses.

In the no recoil case the modified drag force²⁰ is given by

$$\begin{array}{c}\hfill F_{\mathrm{drag}}={|\mathrm{\Delta }{\mathit{v}}_{\mathit{ij}}|}^2\frac{{\sigma }_{\mathrm{T}}}{m_{\chi }}{m}_i{m}_j{\mathrm{\Lambda }}_{\mathit{ij}}.\end{array}$$(H.1)

Furthermore, in re-adding the energy, one should use

$$\begin{array}{c}\hfill \frac{2\mathrm{\Delta }E}{m}=|\mathrm{\Delta }{\mathit{v}}_{\mathrm{drag}}|\left(2|\mathrm{\Delta }{\mathit{v}}_{\mathit{ij}}|-|\mathrm{\Delta }{\mathit{v}}_{\mathrm{drag}}|\right).\end{array}$$(H.2)

With this, we can modify the fSIDM routine and transform it from recoil to no-recoil case.

Appendix I: Fixed-angle deflection test

In Sect. A.2 we discussed a validation test of the hSIDM approach, using the differential cross section Eq. (4) that is peaked at a single value of the deflection angle θ = θ₀. Thus, all individual scatterings occur with this fixed angle. We used this approach to check the validity criterion for the effective treatment of small-angle scattering by a drag force and momentum diffusion derived in Sect. 3.2.2. In the context of the deflection test set-up from Sect. A.2, this corresponds to the minimal time t_min(θ₀) (see Eq. (32)) after which the distribution function f(t, θ) agrees with the Gaussian Molière approximation Eq. (27). While t_min(θ₀) is shown in Fig. A.3 in the main text, we additionally provide more detailed simulations results of the deflection test problem when assuming a differential cross section given by Eq. (4). Specifically, we show the resulting angular distribution functions f(t, θ) at various times t in Fig. I.1, and compare them to the Molière approximation Eq. (27). We see that both agree for late enough times, as stated in the main text.

Appendix J: Multiple scattering and time step

The sampling of the scattering angle in the rSIDM scheme is based on the assumption that the physical particles represented by the numerical particles do not scatter more than once per pairwise interaction per time step. This assumption is implicitly made when sampling the scattering angle for the two numerical particles from the differential cross section instead of a distribution that takes multiple scatterings into account such as Eq. (18) for the no-recoil case. In the following, we test which role this assumption plays in the case of our fixed-angle deflection test, see Sect. A.2.

When considering the pairwise interaction of numerical particles, the size of the time step affects how many of the represented physical particles would scatter. It also alters how many physical particles scatter multiple times relative to only once or not at all. We note that this refers to (1) scatterings of physical particles with in a single time step and (2) scatterings of physical particles belonging to the same pair of numerical particles.

To estimate the probabilities for no and one scatter, we assume that the scattering probability is constant over time:

$$\begin{array}{c}\hfill \frac{\mathrm{d}P}{\mathrm{d}t}=\mathrm{const}>0.\end{array}$$(J.1)

The probability not to scatter if we have only one time step over an interval Δt is

$$\begin{array}{c}\hfill P_{\mathrm{ns}}=1-\frac{\mathrm{d}P}{\mathrm{d}t}\mathrm{\Delta }t.\end{array}$$(J.2)

Given that we have n equally spaced time steps of size Δt/n, it becomes (because of Bernoulli’s inequality)

$$\begin{array}{c}\hfill P_{\mathrm{nsm}}={(1-\frac{\mathrm{d}P}{\mathrm{d}t}\frac{\mathrm{\Delta }t}{n})}^n>P_{\mathrm{ns}}.\end{array}$$(J.3)

Here P_nsm > P_ns for n > 1. Hence, smaller time steps slightly increase the number of unscattered particles. The probability to scatter once can be expressed as

$$\begin{array}{c}\hfill P_{\mathrm{os}}=\frac{\mathrm{d}P}{\mathrm{d}t}\mathrm{\Delta }t.\end{array}$$(J.4)

Given that we have n equally spaced time steps, it becomes

$$\begin{array}{c}\hfill P_{\mathrm{osm}}=n\frac{\mathrm{d}P}{\mathrm{d}t}\frac{\mathrm{\Delta }t}{n}{(1-\frac{\mathrm{d}P}{\mathrm{d}t}\frac{\mathrm{\Delta }t}{n})}^{n-1}<P_{\mathrm{os}}.\end{array}$$(J.5)

Here P_osm < P_os for n > 1. Hence, smaller time steps slightly decrease the number of particles that have scattered only once. Furthermore, there is a recent study by Fischer et al. (2024b) that gives the limit for n → ∞. We note that these equations are not straightforwardly applicable to the simulations, for example, a numerical particle can interact with multiple neighbours per time step and N_ngb is not taken into account in the equations. However, the equations provide a qualitative understanding and thus provide insight into the limitations of the simulations.

To illustrate these limitations, we show in Fig. J.1 the fixed-angle deflection test simulated with a varying number of time steps for a given time t. More precisely, we simulate this deflection test with a fixed angle of θ₀ = 0.2 for a reference time step of Δt = 0.0012 Gyr and other equally spaced time steps of Δt ∈ {0.1536, 0.0768, 0.0384, 0.0192, 0.0096, 0.0048, 0.0024} Gyr. We use the same simulation set-up described in Sect. A.1, a transfer cross section of σ_T/m_χ = 100 cm² g⁻¹ and N_ngb = 16 neighbouring particles. The angle distribution is the (already normalised) mean value of 8 simulations employing different seeds for the pseudo-random number generator.

The three panels of Fig. J.1 show the distribution of the deflection angle at different times. The first peak in the plot (from left to right) corresponds to the unscattered particles, the second peak to the particles that have scattered once, and the third peak is mainly particles that have scattered twice. The height of the first peak decreases over the course of the simulation, while the height of the second peak first increase and later decreases. It is clearly visible that the height of the peaks depends on the time resolution, that is, how many time steps we use to arrive at the solution. In line with the previous estimates, we find that the number of particles that have not scattered yet increases with increasing time resolution. In contrast, for the particles that have scattered once the picture is more complicated, in an initial phase a better time resolution leads to a smaller peak while at a later stage, it leads to a higher peak. At that late time also the peak height is decreasing with time. We note that these findings are in line with Eq. (J.5).

For each of our simulation times (set by multiples of the largest time step), t ∈ {0.1536, 0.3072, 0.4608, 0.6144} Gyr we can compute the ratio of peak heights between the different time resolutions and the reference (Δt = 0.0012 Gyr) for both peaks to better see the multiple scattering effect. We show the result for the first (second) peak in the left (right) panel of Fig. J.2 as a function of time. Again, we see from the first peak that the number of time steps we use clearly influences the number of unscattered particles. For both peaks, we find that the relative difference is increasing with time at the later times we show here.

We can conclude that one needs time steps not only small enough to ensure P_os < 1 but P_os ≪ 1 to obtain accurate results, because numerically we can only describe the single scatterings but not the multiple ones of physical particles. How important this is for a simulation, depends on how crucial resolving the precise angular dependence for the evolution of the SIDM distribution is. While we have chosen a sensitive test, this problem might be less relevant for typical astrophysical simulations. In consequence, the optimal choice of the time step to achieve a good trade-off between accuracy and performance might be problem dependent. Thus the value for κ in the time step criterion (Eq. (43)), which effectively limits the interaction probability, could be chosen depending on the problem at hand.

All Tables

All Figures