Simulations of galaxy cluster mergers with velocity-dependent, rare and frequent self-interactions

Self-interacting dark matter (SIDM) has been proposed to solve small-scale problems in $\Lambda$CDM cosmology. In previous work, constraints on the self-interaction cross-section of dark matter have been derived assuming that the self-interaction cross-section is independent of velocity. However, a velocity-dependent cross-section is more natural in most theories of SIDM. Using idealized $N$-body simulations without baryons, we study merging clusters with velocity-dependent SIDM. In addition to the usual rare scattering in the isotropic limit, we also simulate these systems with anisotropic, small-angle (frequent) scatterings. We find that the collision-less brightest cluster galaxy (BCG) has an offset from the DM peak that grows at later stages. Finally, we also extend the existing upper bounds on the velocity-independent, isotropic self-interaction cross-section to the parameter space of rare and frequent velocity-dependent self-interactions by studying the central densities of dark matter-only isolated haloes. For these upper-bound parameters, the DM-BCG offsets just after the first pericentre in the dark matter-only simulations are found to be $\leq$ 10 kpc. On the other hand, because of BCG oscillations, we speculate that the distribution of BCG offsets in a relaxed cluster is a statistically viable probe. Therefore, this motivates further studies of BCG off-centring in hydrodynamic cosmological simulations.


INTRODUCTION
Cold dark matter (CDM) is a fundamental component of the standard ΛCDM cosmology.It plays a vital role in explaining the formation of the large-scale structure of the universe and the anisotropies in the cosmic microwave background.While cosmological -body simulations within ΛCDM have successfully reproduced many observations of the large-scale structure, there seem to be discrepancies between observations and simulations on small scales (see Bullock & Boylan-Kolchin (2017) for a review of the small-scale problems).A solution to the small-scale problems was proposed by Spergel & Steinhardt (2000) via a model of DM, where DM particles can non-gravitationally scatter off each other.Constraints on the selfinteraction cross-section can be obtained by studying different astrophysical systems.In particular, relaxed galaxy clusters (e.g., Sagunski et al. 2021;Andrade et al. 2021) have provided the most stringent constraint on the cross-section.We also have constraints from galaxy cluster mergers (Randall et al. 2008;Harvey et al. 2015).For a review on astrophysical constraints on self-interacting dark matter (SIDM) see Adhikari et al. (2022).
Velocity-dependent anisotropic cross-sections are natural in most theories of SIDM (for a review of SIDM models see Tulin & ★ sabarish.venkataramani@uni-hamburg.deYu 2018).Examples of such models include light mediator models (Tulin et al. 2013;Ackerman et al. 2009), atomic DM (Cline et al. 2014), strongly interacting DM (Boddy et al. 2014).Moreover, velocity-dependent cross-sections are also well motivated observationally.Most constraints on the self-interaction cross-section per unit mass of DM particle (/  ) in the literature, have been derived assuming velocity-independent and isotropic scattering.For example, Sagunski et al. (2021) quote an upper limit of /  < 0.35 cm 2 g −1 (95%C.L.) at cluster scales and /  < 1.1 cm 2 g −1 (95%C.L.) at group scales.On the galactic scales, Ren et al. (2019) find that /  is required to be in the range 3-10 cm 2 g −1 to explain the observed diversity in the rotation curves in the SPARC data set.Similarly, Sankar Ray et al. (2022) quote an upper-bound of /  < 9.8 cm 2 g −1 (95%C.L.).At the scale of dwarf galaxies, there is no concrete upper bound on the crosssection.A cross-section with /  > 30 cm 2 g −1 is favoured by the observed central densities of Milky Way's dwarf spheroidal galaxies (Correa 2021).Elbert et al. (2015) find that /  can be as large as 50 cm 2 g −1 at these scales and still be consistent with observations.These considerations highlight the viability of velocity-dependent cross-sections.
Observationally probing the angular dependence is a daunting task.One reason being that the effects of angular dependence are not strong enough when studying the evolution of systems that do not have a pre-ferred direction.For example, Robertson et al. (2017b) and Fischer et al. (2021) simulated isolated haloes using -body simulations and found that there is no big difference in the evolution of core-sizes between isotropic and anisotropic cross-sections for a given choice of parameters.Moreover, simulating differential cross-sections which peak for tiny scattering angles, using conventional SIDM implementations (e.g.Rocha et al. 2013) is prohibitively expensive.This type of interaction will be called frequent self-interactions (as introduced in Kahlhoefer et al. 2014) as opposed to rare self-interactions corresponding to large-angle scattering.Frequently self-interacting dark matter (fSIDM) is natural in the mass-less mediator limit of light mediator models.The scattering of DM particles in this regime can be modelled as a drag force, where the drag force depends not on the total cross-section but on the momentum transfer cross-section given by Kahlhoefer et al. (2017) and Robertson et al. ( 2017b) where  cms and Ω cms are the scattering angle, and solid angle in the centre of mass frame.Frequent self-interactions have previously been studied by Kahlhoefer et al. (2014), Kummer et al. (2018Kummer et al. ( , 2019)), Fischer et al. (2021Fischer et al. ( , 2022) ) assuming velocity independence.Mergers of galaxy clusters are interesting test beds for models of SIDM since the mass distribution of the system could be measured through lensing.The gas and galaxies can be probed through their direct emission in various wavelengths.The existence of offsets between the DM component and galaxies may hint at DM selfinteractions (Randall et al. 2008).Moreover, mergers are sensitive to both velocity and angular dependence of the scattering cross-section.First, as the haloes undergo many pericentre passages, the collisional velocity changes with time.Scattering velocities are the largest at the first pericentre passage and, subsequently, the haloes slow down with every passage.The self-interactions at the pericentre passages are mainly responsible for an increase in the offset.Thus, the evolution is sensitive to the parameters of velocity-dependent cross-section.Secondly, mergers unlike isolated haloes also have a preferred direction, i.e. the merger axis.Fischer et al. (2022) find that offsets are larger for frequent self-interactions with a given  T , when compared to rare self-interactions of the same  T .They also showed that smallangle scattering can produce larger offsets than the maximal possible offset from isotropic scattering.
There have been earlier studies that have simulated mergers.For example, studies with velocity-independent isotropic cross-sections have been done by Kim et al. (2017) who simulated equal-mass mergers ; Robertson et al. (2017a) simulated a bullet cluster like system.Fischer et al. (2021Fischer et al. ( , 2022) ) studied equal and unequal-mass mergers with isotropic and anisotropic velocity-independent cross-sections.Robertson et al. (2017b) looked at mergers until just after the first pericentre passage.They used a velocity-dependent isotropic crosssection and a velocity-dependent cross-section that corresponds to Yukawa scattering under the Born approximation.It is unknown as to how the merger evolution is affected at late stages by velocitydependent self-interactions.Similarly, mergers in the fSIDM regime with velocity dependence are yet to be studied.
In this work, we aim to (i) study the qualitative differences in merger simulations between velocity-dependent and independent cross-sections, (ii) extend the upper bound on constant cross-section quoted by Sagunski et al. (2021) to the parameter space of velocitydependent cross-section, (iii) find the maximum offsets between DM and the brightest cluster galaxy (BCG) with velocity-dependent cross-section parameters that are consistent with upper bound parameters.To this end, we simulate the full evolution of galaxy cluster mergers and isolated haloes with rare and frequent self-interactions.In a companion paper (Fischer et al. 2023a), cosmological simulations are studied with velocity-dependent fSIDM.The paper is presented as follows.In section 2 we briefly describe our numerical scheme and the SIDM models that are considered.In section 3 we present our simulation results, which illustrate the qualitative differences between velocity-dependent and velocity-independent crosssection simulations.In section 4, we describe the simulations of mergers with cross-section parameters that correspond to the 95% C.L. limits provided in Sagunski et al. (2021).In section 5, we summarize our results and conclude.Additional details are provided in the appendices A and B.

METHODS
In this section, we describe the numerical setup of our simulations, and we discuss our choice for the self-interaction cross-section that is used in the simulations.

Numerical method
For our simulations, we use the cosmological -body simulation code OpenGadget3, adapted for frequent self-interactions using the implementation given in Fischer et al. (2021).In this section, we briefly describe the implementation of rare and frequent selfinteractions.The numerical scheme for the self-interactions of rarely self-interacting dark matter (rSIDM) follows Fischer et al. (2021), which is similar to the method described in Rocha et al. (2013).In this scheme, the probability that a numerical particle 1 with mass  1 scatters off another numerical particle 2 with mass  2 is given by, where  12 is the relative velocity between the numerical particles 1 and 2, Δ is the time-step used in the simulation, Λ 12 is the kernel overlap integral, ( 12 )/  is the total cross-section per unit mass of DM particle.For more details on the implementation and the choice for the kernel, see appendix A and B of Fischer et al. (2021).
A collection of kernels used in other modern implementations of SIDM within -body simulations can be found in Adhikari et al. (2022, equations 11-15).The drag force of frequent self-interactions is based on the relation derived in Kahlhoefer et al. (2014), which describes the deceleration rate experienced by a particle as it travels through a constant background density of DM, 0 is the velocity of the particle,  ∥ is the parallel component of the velocity of the particle,  0 is the background density,  T is the momentum transfer cross-section defined in equation (1).We can also see that the deceleration rate captures the rate of change of the parallel component of the velocity.Therefore, the above expression can then be cast into an expression for drag force for the physical particles.
The drag force as experienced by numerical particles, labelled 1 and 2, in the -body code is given as (Fischer et al. 2021), which is proportional to the momentum transfer cross-section  T .
From here on, we will drop the subscripts and let  denote the relative velocity between DM particles.In addition to the drag force, momentum is added to the particles in a random but perpendicular direction relative to the motion in the centre-of-mass frame.As for rSIDM, the post-scattered momenta in the centre-of-mass frame have the same absolute value as the pre-scattered ones.
In our SIDM implementation, we search for pairs of potentially interacting particles.A search for the neighbours of each particle achieves this, employing the same tree structure as used in the gravity computation.
The code employs an adaptive time-stepping scheme, where the time-step of an individual particle is obtained by the minimum of different time-step criteria, i.e. given by gravity and self-interactions.
The implementation of rSIDM and fSIDM does conserve momentum as well as energy explicitly.This is achieved by executing the scattering computations for pairs of particles in a consecutive manner such that their velocities are updated after each scattering event and used for the next one.At the same time, we allow for multiple interactions per particle per time-step.
Massively parallel computations of the self-interactions are enabled by a parallelized implementation using the message passing interface (MPI).Note, here we order the communication and computation of the processes in a manner to reduce waiting time.
The implementation of the velocity-dependent self-interactions into the SIDM module has been described in detail by Fischer et al. (2023a).To ensure numerically stable results, a novel time-step criterion has been added.This criterion is based on the velocity at which self-interactions are strongest, i.e. on the maximum of  T () , instead of the velocity distribution that an individual particle encountered in the previous time-step.In principle, such a scheme guarantees that the simulation time-step is always sufficiently small to account for  T () for any .

Initial conditions and simulation parameters
In this work, we simulate both DM-only isolated haloes and galaxy cluster merger with two different merger mass ratios (MMR).Firstly for the isolated DM-only haloes, they have a virial mass of 10 15 M ⊙ and are initialized with an NFW density profile Navarro et al. (1996), i.e.

𝜌(𝑟)
(5) DM positions are sampled from a probability density function such that the density profile follows the NFW profile.Given the virial mass of the halo and the critical density  crit of the universe, the parameters of the NFW profile are determined as follows: (i) the concentration parameter  vir is determined using the concentration-mass relation given in Dutton & Macciò (2014), (ii) characteristic density   is computed using  vir from the earlier step (using Eq.2 of Navarro et al. 1996), (iii) scale radius   is computed using its definition   =  vir / vir .The NFW parameters thus obtained are given in Table 1.
Using the resulting density profile, the initial velocity dispersion  2 ini () is obtained by integrating the Jeans equation (Binney & Tremaine 2008).Then, initial velocities in each radial bin are drawn from a Gaussian distribution with the variance  2 ini ().The isolated DM halo simulations use the same NFW parameters as the main halo of the merger for initial conditions (ICs).We simulate mergers with MMR ∈ {1, 5}.The barycentre of the clusters are initially 4000 kpc apart and they are put on course towards each other with a relative velocity of 1000 km s −1 as summarized in Table 2.In the galaxy cluster used in the merger simulations, the cluster has three particle species: DM, galaxy and BCG.Galaxies and BCG are approximated to be collision-less, while DM is collisional with self-interaction characterized by its cross-section.An equal number of DM and galaxy particles is used in the simulation.A sufficient number of galaxy particles are chosen to ensure that it is easier to find the peak position of the galaxies.The main halo has a virial mass of 10 15 M ⊙ and both species initially follow an NFW profile.The particle masses are as follows: for DM,  DM = 2 × 10 9 M ⊙ , for galaxies,  Gal = 4×10 7 M ⊙ .In addition, the brightest cluster galaxy (BCG) is represented by a single particle at the centre of the halo with a mass  BCG = 7×10 11 M ⊙ .As the BCG is approximated to be a point particle, the effects of gravitational scattering become strong if the mass is large, therefore the BCG particle is taken to be less massive than the observed BCGs.This choice is adopted from earlier studies (e.g.Fischer et al. 2022;Kim et al. 2017).The mass resolution chosen works reasonably well for measuring offsets in simulation since the simulation is run only for a few dynamical time-scales.The detailed effects of such a treatment on measurements other than offsets are yet to be studied.We use a fixed gravitational softening length of  = 1.2 kpc for all particles.For both mergers and isolated haloes we use an adaptive kernel size for the DM self-interactions, such that the number of neighbours within each particles' kernel,  ngb is equal to 64.This choice follows from Fischer et al. (2021).
All simulations have been performed with a resolution of O (10 6 ) particles.For certain cross-section parameters, the simulations were rerun at a higher resolution of O (10 7 ) particles to validate the lower resolution runs (see Appendix A).

Dark matter cross-section model
We assume a fairly generic form for the self-interaction momentum transfer cross-section for rare and frequent self-interactions, parametrized as follows (Gilman et al. 2021;Yang et al. 2023), where  is the relative velocity of the DM particles.We furthermore assume that the total cross-section has the same velocity dependence as the momentum transfer cross-section, such that This assumption is automatically satisfied if the differential crosssection is a separable function, i.e. it is of the following form, where  is a normalization constant, Θ() captures the angular dependence and (, ) = (1 +  2 / 2 ) −2 .However, even for nonseparable differential cross-sections, such as Rutherford scattering, the assumption in equation ( 7) gives a useful approximation.
Depending on the choice of  and Θ(), the above differential cross-section can correspond to either frequent or rare selfinteractions.In this work, we consider isotropic scattering for the case of rare self-interactions, such that Θ() is simply a number.The total cross-section to calculate the scattering probability in equation ( 2) is then given by () = 2 T (), see equation ( 1).For frequent self-interactions, on the other hand, Θ() is strongly peaked for small angles.The normalization, , is then chosen such that the momentum transfer cross-section given in equation ( 6) is recovered, which is used in the simulations to compute the drag force experienced by the particles, see equation (4).
To run the simulations, the parameters  0 and  must be chosen, where  0 :=  0 /  .In this paper, we are interested in studying the qualitative differences in the evolution of galaxy cluster mergers between the different regions of  0 −  parameter space.A priori, it is conceivable that there are degeneracies in this parameter space, such that different parameter combinations lead to very similar merger observables.Such a degeneracy was found by Yang et al. ( 2023) (see their fig.8) when analysing the rotation curve of LSB dwarf galaxy UGC 128.In other words, it was possible to compensate a change in  by an appropriate change in  0 .We refer to a prescription to determine the cross-section normalization  0 for a given velocity dependence (determined by ) as matching.In the following, we will explore different matching procedures.
The choice for  follows from the typical relative velocity scales.These scales can be estimated by running CDM simulations.The observed values are displayed as dashed lines in Figure 1.The largest observed scales are the relative velocities of the infalling BCGs shown as the last two lines corresponding to MMR:5 and MMR:1 respectively.The average scattering velocity with which DM particles scatter off each other within 100 kpc around the barycentre at the first pericentre passage are shown.Finally, the relative velocity dispersion within 100 kpc of a 10 15 M ⊙ halo is displayed as the leftmost vertical dashed line and it is approximately 1400 km s −1 .This implies that for any value of  larger than 1400 km s −1 , the self-interactions within the halo will be in the weakly velocitydependent regime.Therefore, the following choice for  is made,  ∈ {1000, 2000, 3000, 4000}km s −1 .
We have analysed cluster mergers with  0 chosen according to two matching procedures: The first is to choose the same value of  0 for all chosen values of , the second is to choose  0 such that the evolution of the central density of the isolated haloes for different values of  is similar.These procedures will be explained further in the following sections.

Analysis methods
In order to find the peak position of any component, i.e.DM or galaxies, we use the peak finding method based on the gravitational potential, see Fischer et al. (2022).In the simulations, all particles have a unique particle ID assigned to them.Using the ID, particles belonging to a given halo can be identified.Then, the gravitational potential in each cell in a grid is computed.The cell with the lowest potential corresponds to the position of the peak.Fischer et al. (2022) also propose the isodensity-based peak finding algorithm.In this algorithm, the peaks are identified as the cell in the merger plane with the highest projected density.This method is closer to observations where, for example, gravitational shear measurements can be used to infer mass densities.We find that gravitational potential based peak finding is more reliable when the simulation is run with low resolution.Hence, this is our choice for finding peaks.To find the errors on the peak position, we bootstrap the particle distribution 20 times and determine 20 such peaks.Then we estimate the error, by finding the standard deviation in the obtained peak positions.
We define offsets as the distance between the peaks of two different species of the same cluster.For example,  DM−BCG is the distance between the DM and BCG peak of a given cluster.

VARYING 𝑤 ONLY
In this scheme, the same value of  0 is used for different values of the parameter .Even though it is a very simple choice, it is easier to observe the qualitative difference introduced by velocity-dependent cross-sections.To make the differences stand out, we use a large value of  0 = 5.0 cm 2 g −1 .We confirmed the results to be qualitatively similar but less pronounced for smaller values of  0 .The names for the individual runs shown in the plots are tabulated in Table 3. Table 3. Simulation labels and the corresponding cross-section parameters.
The generic form of labels is XwDsNpM, where X is F(frequent), R(rare), or C(constant).D that follows w is the value of the parameter  to be read as D km s −1 and NpM following the s is the value of  0 to be read as N.M cm 2 g −1 .
2.0 2.5 3.0  3. Simulation results presented in this figure corresponds to varying values of  and a fixed  0 . .

DM peak position
The plots of the peak positions of DM against time for equal and unequal-mass mergers are given in Figure 2. The peaks of the main halo are marked by dotted lines, while those of subhalo are indicated by solid lines.First we observe that for constant   the drag force from self-interaction is strong enough to stop the DM component in their tracks and they coalesce at the first pericentre passage.For the same  0 , the DM peak positions in the velocity-dependent SIDM runs are closer to the CDM run for smaller values of .This observation matches our expectation that, for fixed  0 , increasing  increases the effective self-interaction strength.The separation of the peaks at the first apocentre is largest for CDM and decreases with increasing self-interaction strength.This can be attributed to the fact that DM particles experience a drag force and thus stay closer to the centre-of-mass of the system.
The first pericentre passage occurs at different times for different parameter combinations.Albeit a small effect, self-interactions tend to increase the time taken to reach the first pericentre passage.Then the second passage occurs earlier for increasing self-interaction strength.For example, we see that the first pericentre passage occurs slightly later in the simulation with constant   .This could be understood as self-interactions generating pressure that resists the infall, thus slowing the halo down.The difference is large for constant   because the effective self-interaction strength of the velocitydependent ones is much smaller than the one of constant   at early stages, as  T () is small for  > .
At later stages of the evolution, the scenario changes.For example, at the third apocentre, we see that the oscillation in the DM peak for  = 2000 km s −1 has an amplitude that is greater than that of CDM.During these stages, the central density of the haloes are lower for SIDM simulations (see Figure C2).Therefore, the oscillations are less dampened by dynamical friction for SIDM simulations.
For the unequal-mass mergers, the subhalo dissolves faster, making it difficult to identify the DM peaks during later stages of the evolution.Both equal and unequal mass mergers have identical initial separation and initial relative velocity.Therefore, in unequal mass merger, the less massive cluster traverses more distance than the massive one and this leads to fewer oscillations.For instance, in Figure 2, we see that within 5 billion years, the subhalo in MMR:5 system has undergone fewer pericentre passages than the equal-mass merger.

BCG peak position
The BCG positions for subhaloes are given in Figure 3, the upper panel corresponds to the equal-mass merger system, while the lower panel corresponds to the unequal-mass merger.In CDM simulations, the BCG oscillations are damped faster compared to SIDM simulations.This general feature has already been observed in earlier work (Fischer et al. 2022(Fischer et al. , 2021;;Kim et al. 2017).This could be explained by noting that the merger remnants have a cored density distribution at the centre owing to the self-interactions.On the other hand, merger remnants in CDM simulations have larger central densities.As a result, the oscillations dampen out faster in CDM due to dynamical friction.
In the equal-mass merger, we observe that the peak positions of subhalo BCGs are closer to the CDM for smaller values of  at the early stages of the merger evolution.At later stages, around 5 billion years, the BCG oscillations in the CDM simulation have dampened considerably.On the other hand, the oscillations approximately stay constant for the constant  T simulation for the period 1 − 8 billion years, as shown in Figure 3.The position of the BCG in velocity-dependent simulations start deviating from CDM with time.For example, let us consider the  = 2000 km s −1 simulation : (i) we see that the curve is initially close to the CDM simulation (ii) during the period 4 − 7 billion years, the oscillations have approximately a constant amplitude, and they have significantly deviated from the CDM simulation.The typical velocities around the DM peak at various stages of the evolution are given in Figure D2.The central relative velocity dispersion around the DM peak is large at the first pericentre owing to the active merging.At later stages, when the haloes have slowed down, they have a very slowly rising relative velocity dispersion for the rest of the simulated time period.Therefore, from Figure 1 we see that for  = 2000 km s −1 and  > 2000 km s −1 , the  T is less than 1 cm 2 g −1 , whereas for  < 2000 km s −1 the crosssection is larger than 1 cm 2 g −1 .Thus, merger remnants experience larger self-interactions at later stages, leading to more cored distributions and lesser dynamical friction.Cumulatively, this leads to steady oscillations at later stages.
The distance travelled by the BCGs at the first apocentre is observed to become smaller with increasing values of .This can be understood by looking at the DM peaks.For example, the DM haloes coalesce at the first pericentre passage for the constant   simulation.As a result, the BCG experiences a larger gravitational force due to the coalesced DM distribution at the barycentre.This accumulation of DM at the barycentre reduces with decreasing  since the average interaction strength reduces with .This effectively leads to smaller amplitudes at the first apocentre in both equal and unequal-mass mergers.Immediately after the second pericentre passage, the amplitude of the velocity-dependent SIDM and CDM curves have decreased significantly due to dynamic friction.In the constant   simulation , the DM peaks have come to rest, and the merger remnants start coring.This leads to the persistence of the BCG oscillation amplitude.

Morphology
In Figure 4, we show the physical density of DM within the slice || < 100 kpc projected on the merger plane.There are three columns, each correspond to a DM model and each row corresponds to a particular simulation time.First we look at time  = 1.4 Gyr, i.e. before the first pericentre passage.Both haloes in the constant   simulation (column 2) have lower central densities than in the CDM (column 1) and velocity-dependent simulations (column 3).The second row shows results at a time  = 2.1 Gyr, which is just after the first pericentre passage.Owing to the large self-interaction strength, the haloes in the constant   simulation have coalesced, while for CDM and velocity-dependent self-interactions the DM haloes pass through each other.At later stages,  > 5 Gyr, the merger remnant in the constant   simulation has a lower central density.While for the velocity-dependent simulations, the effects of velocity-dependence slowly becomes relevant as the system slows down.Thus, this leads to more cored distribution at the centre of the merger remnant when compared to the one from CDM simulation.This feature essentially leads to the persistent BCG oscillations at late stages.In Figure 5 we have a similar plot displaying only the subhalo's projected density for MMR:5.Independent of the choice for the   (), the subhaloes are observed to evaporate with time.With self-interactions, the evaporation is more pronounced.At  = 4.5 Gyr, the subhalo experiencing constant   has its core dissolved significantly and comes to rest.For CDM, the core has remained relatively intact.Finally, in the velocitydependent simulation, the core has dissolved.Although the merger remnant seems to be oscillating even at these stages.In addition, in the CDM simulation, we see shell-like features.These features  3.
are missing in the constant   simulation, since the haloes have coalesced.

CENTRAL DENSITY MATCHED CROSS-SECTION
In this section, we explore a more refined matching procedure based on simulations of isolated haloes.In this matching scheme, parameters { 0 , } are chosen such that different parameter combinations lead to similar central density evolution.We will refer to parameters matched according to this scheme as CD-matched.To avoid performing multiple simulations to find the matched parameter set, we make use of the self-similar nature of the gravothermal fluid equations of an isolated halo (Balberg et al. 2002;Essig et al. 2019).This allows us to obtain the central density evolution without running a suite of simulations.In Balberg et al. (2002), they assume that the total crosssection is velocity independent.In order to illustrate the rescaling, consider two constant cross-section parameters   0 ,   0 .Then, the central density evolution obeys the following scaling relation: (5, ∞), (5,2000).The rows correspond to different times.The first row being before the pericentre passage, the second, just after the first pericentre passage and the third being at later stages.The time stamps in the images are in units of Gyr.
where   ,   correspond to the evolution time of the isolated haloes simulated with parameters   0 ,   0 .Velocity-dependent cross-sections contain two parameters, and it is not immediately clear how the central densities can be rescaled.Yang & Yu (2022) propose an effective cross-section  eff to model the halo evolution.For a differential cross-section d/d cos , the effective cross-section is given by In the expression given above,  is the relative velocity of DM particles,  1D is the characteristic velocity dispersion of the halo.Yang & Yu (2022) show that the evolution of central density in a simulation with the differential cross-section can be mimicked by a constant cross-section simulation with the same  eff .They also note that the equivalence holds well when the halo is in short-mean-free-path regime.In the long-mean-free-path regime,  eff does not capture the effects of self-scattering accurately.However, it provides a reasonable approximation to the halo evolution.Therefore, we extend the rescaling procedure given in equation ( 9) to any differential cross-section by using the corresponding effective cross-section  eff .
Integrating the angular part of equation ( 10), we get where, is the viscosity cross-section.Thus, for a given ,   eff /  eff =   0 /  0 .In other words, rescaling by  eff is equivalent to rescaling by  0 for the same .We verify this method by performing tests with some parameter combinations, see Appendix B.
In order to find the value  0 given a , such that the evolution of the central density matches that of a target simulation with parameter set,  = {  0 ,   } we follow the procedure given below.
(i) Simulate an isolated halo with the target parameter set.Find the evolution of the central density    (  ) from the simulation snapshots.
(ii) Simulate an isolated halo with the parameter set  = {  0 , }, followed by the estimation of the evolution of central density    (  ) from the simulation data.
(iii) To obtain the evolution    (  ) corresponding to the parameter set  = {  0 , }, rescale the time axis of the simulation A, i.e., Thus, we have obtained the CD-matched   0 for the given .To make a guess for the value of   0 for one of our chosen values of , we solve, To this end, we need a value for  1D .In Yang & Yu (2022) show that the 1D velocity dispersion is 0.64 max at the maximal core expansion phase.Using  max ≈ 1.65  √︁   0 for our NFW parameters, we find  1D ≈ 980 km s −1 .On the other hand, we observe  1D ≈ 1000 km s −1 in our simulations.Thus, we find consistency between our simulations and the semi-analytic result.
We calculate the initial guess   0 for rare self-interactions and use the same for frequent self-interactions.For this calculation, we use the differential cross-section given in subsection 2.3.
We choose the target set  to correspond to the values quoted in Sagunski et al. (2021).They quoted a 95% upper limit on /  of 0.35 cm 2 g −1 , assuming isotropy and velocity independence.Therefore, this value translates to  0 = 0.175 cm 2 g −1 for rare self-interactions (because  = 2 T ).In other words, we choose the target set  = {0.175cm 2 g −1 , ∞}.In Figure 6, we show an example for the matching procedure in detail.The black band corresponds to the target central density of constant   with  0 = 0.175 cm 2 g −1 .The orange band corresponds to a simulation with {2.72 cm 2 g −1 , 2000 km s −1 }.This value for  0 is our initial guess calculated using equation ( 14).After rescaling by trial and error, the desired value of  0 is found to be 1.2 cm 2 g −1 .This procedure can be extended to all chosen values of  and the obtained results are tabulated in Table 5.The inferred values of  0 at different values of  can then be used to calculate the corresponding viscosity cross-section   .This is shown in Figure 7.The orange triangles and blue stars represent the values of   obtained using the inferred values of  0 from -body simulations for rare and frequent self-interactions, respectively.Similarly, the solid line corresponds to the   calculated from the values of  0 inferred by solving equation ( 14).The fact that the results obtained from  eff and -body simulations are different can be attributed to the fact that the isolated halo is in the long-mean-free-path regime.This was already noted in Yang & Yu (2022).
The evolution of an isolated halo has a feature that, as long as it is in the core expansion phase, at any given time the central density is monotonically decreasing with  eff .This implies that for a given 1000 2000 3000 4000 5000 Rare -NBody Frequent -NBody  2022) constrain  0 using the observed core-size in clusters.Hence, for every value of , there is a value of  0 that produces core sizes, or central densities, similar to the current upper bound.Increasing  0 any further would increase the core size to values larger than what is observed.Thus, this matching procedure can be used to extend the bounds from constant   to different values of .
When matched using  eff or central density evolution, we observe that the ratio between the  0 's of rare and frequent is approximately 0.64 at every chosen values of .This can be understood from the definition of  eff .From equation (11), we have for any , for rare self-interactions, Here,  = 1/(512 8 1D ).Now from equation ( 6), for rare selfinteractions we have  T =  0 (, ), which implies that For frequent self-interactions, we consider a differential cross-section of the form given in equation ( 8) with the function Θ() having support only for values of  close to 0. A simple choice for Θ() is a step-function that is non-zero in the interval [0, ], where  is some small number.Therefore,  eff is given as Table 4.This table contains labels and cross-section parameters matched to constant cross-section,  0 = 1.0 cm 2 g −1 of frequent self-interactions using central density matching scheme.
where in the second equality we have Taylor-expanded and retained only the leading order in  in the angular integrand.Similarly, for the momentum transfer cross-section we have ≈  0 (, ) 4 /8.

Central density matched simulations -qualitative features
In this subsection, we look at the qualitative differences in mergers when parameters are chosen according to the central density matching procedure.Again, as in section 3, we only simulate frequent selfinteractions.In section 3, we investigated the effects arising from changing the value of .To ensure that the effects of self-interaction are not negligible for  = 1000 km s −1 , we used a large value for  0 of 5.0 cm 2 g −1 .On the other hand, in this section the crosssection parameters are CD-matched for which we choose the target set to be  = {1.0cm 2 g −1 , ∞}.The matched parameters are given in Table 4.We have chosen a value for  0 such that the effects of self-interaction are observable, but not as large as the previous value 5.0 cm 2 g −1 .
For brevity, we show only the BCG oscillations in Figure 8.We see that the BCG oscillations of velocity-dependent simulations all have similar amplitudes at early stages.Only at a later stages they start to deviate.This similarity stems from the fact that the parameters are CD-matched.Similar to what was explained in subsection 3.2, at early times the DM particles have a large relative velocity.As a result, most of the velocity-dependent cross-sections have a small effective self-interaction strength.At later stages, the system slows down and  ≲  and the effective self-interaction strength increases.In addition, the internal evolution around a DM peak is similar for all the cross-section parameters chosen, since they are CD-matched.See appendix C for the evolution of the central density of the main halo.

Central density matched upper bound simulations
In section 4, the conservative upper bounds were obtained using the central density matching scheme.Values are tabulated in Table 5.We run merger simulations for this set of parameters and estimate the offsets  DM−BCG .For velocity-independent cross-sections up to 0.5 cm 2 g −1 , Fischer et al. (2021) found that the DM-BCG and DM-Galaxy offsets increase with increasing values of  0 .By running merger simulations at the upper bound values of the cross-section parameters, we estimate the order of magnitude of the largest possible offsets allowed by current bounds.
In Figure 9, we show the DM-BCG offset at the first pericentre passage for the equal mass merger.The offset is O (1)kpc, while the offsets after the third pericentre passage ( ∼ 4 billion years) start increasing and is seen to be O (10)kpc.The effective self-interaction strength is not large enough to produce an observable offset at the first pericentre.Hence, it might be difficult to find such an offset in real observations.We do not show the offsets for unequal mass merger because, the offsets just after the first pericentre are smaller than the ones of equal mass merger.In addition, at later stages due to the evaporation of halo, we do not have reliable peak positions.We can also see from Figure 9 that the offsets produced by fSIDM are  larger than that of rSIDM, and this is due to the fact that mergers are sensitive to the angular dependence.

CONCLUSIONS
We first discuss the assumptions made in the paper before we conclude.The first assumption that we make is that we use idealized initial conditions.Yet, it is informative to study them to find the appropriate features to look out for in more realistic simulations and observations.One example is the amplitude of BCG oscillation at late stages, which is seen to depend on the velocity-dependent crosssection parameters.Therefore, it is instructive to simulate cosmological boxes with velocity-dependent cross-section at higher resolution in order to estimate the distribution of DM-BCG offsets.Later, this could be compared to observations (e.g.Cross et al. 2023;Lauer et al. 2014).For example, Harvey et al. (2019) studied the DM-BCG offsets in the BAHAMAS-SIDM suite of cosmological simulations and placed constraints on /  assuming velocity-independent isotropic SIDM.
In addition, we have modelled the BCG and galaxies as collisionless point particles.A more realistic treatment of BCG's and galaxies is necessary for better comparison with observations.Furthermore, we neglect the effect of galaxies having their own DM halo (Kummer et al. 2018).In reality, approximately 10% of the mass of clusters is made up by the intracluster medium (ICM), which is not included in our simulations.Merger studies that include the ICM, such as Robertson et al. (2017a), find that the DM-galaxy offset is not significantly affected by the presence of the ICM.Fischer et al. (2023b) also finds similar results.Therefore, we argue that the absence of an ICM component does not significantly affect our conclusions.We leave the study of these systems with hydrodynamic simulations to a future study.
As mentioned in subsection 2.3, we have assumed that the angular and velocity dependence of the differential cross-sections can be separated into two functions.This assumption has to be dropped when dealing with realistic SIDM models (Tulin et al. 2013;Feng et al. 2009).There are other SIDM models leading to different effects that are not included in our studies.For example, in addition to elastic scattering, inelastic scattering could be included (O'Neil et al. 2023).We leave the study of such models to future work.
On the observation side, all observed DM-BCG offsets are inferred to be just after pericentre passage and have uncertainties that make them consistent with zero (Bradač et al. 2008;Dawson et al. 2012;Dawson 2013;Jee et al. 2014Jee et al. , 2015)).Similarly, estimating DMgalaxy offsets are difficult due to shot noise arising from the smaller galaxy count.In our simulations we had 10 6 galaxy particles, but in reality we observe at most ∼ 10 3 of them.On the other hand, Lauer et al. (2014) find a median offset of ∼ 10 kpc, with the offset measured between BCG and cluster centre, the latter being identified by X-ray observations.Their sample comprised 433 BCGs that are located in Abell galaxy clusters.The DM-BCG offset distribution could be compared to predictions of cosmological simulations.Overall, the situation with observations is expected to improve with forthcoming surveys, such as SuperBIT (Romualdez et al. 2016) and Euclid (Laureĳs et al. 2011).
We have simulated idealized, isolated haloes and galaxy cluster mergers with equal and unequal merger mass ratios with velocitydependent, frequent and rare self-interactions.Mergers are interesting astrophysical probes since the system is sensitive to selfinteraction cross-sections with both angular and velocity dependencies.Therefore, we focused on understanding the qualitative effects that arise from velocity dependence in mergers.On the quantitative side, we also investigate the maximum offsets that can be observed given the current bounds on /  .
• Independent of the matching procedure used in the paper, the effects of velocity-dependent cross-sections can be observed on galaxy cluster mergers by comparing the early time and late time oscillations of BCG.In particular, the degeneracy in the cross-section parameters when studying the evolution of central density in isolated haloes is broken when studying mergers.This is due to the fact that the relative velocities of the merging clusters change with time.
• The evolution of central densities of isolated haloes are similar between rare and frequent self-interaction, when the momentum transfer cross-section  T () of fSIDM is chosen to be 2/3 T () of rSIDM.The factor 2/3 follows from matching the angular dependence of fSIDM and rSIDM with viscosity cross-section, as seen in Figure 7.
• We extend the existing upper bounds on the constant crosssection  0 to the parameter space { 0 , } of velocity-dependent, rare and frequent self-interactions.
• In the equal-mass merger simulations with upper-bound crosssection parameters, we find that the offsets after the first pericentre is approximately O (1)kpc.In particular, the offsets are the largest in the constant cross-section simulation.As the system evolves further, offsets grow.After the third pericentre passage, due to the oscillations of BCG, and the galactic component the offsets are O (10)kpc.Thus, mergers in their late stages are interesting to test and constrain SIDM.
In conclusion, we have studied the qualitative effects of velocitydependent SIDM cross-sections in galaxy cluster mergers.Our models do not have the realism required for a direct comparison with astronomical data, owing to the neglection of baryonic effects.However, they offer insights into the physical processes that govern the phenomenology of SIDM.More realistic predictions can be obtained by performing full hydrodynamical simulations that include stars, cooling and feedback effects.The significantly larger complexity of such models with additional degrees of freedom render the interpretation much harder.Clearly, this is an avenue for future work.

APPENDIX E: MOMENTUM TRANSFER CROSS-SECTION OF CD-MATCHED PARAMETERS
Plots similar to Figure 1 but with CD-matched parameters are provided.Figure E1 and E2 correspond to the parameters in Table 4 and 5, respectively.3.

Figure 1 .
Figure1.Momentum transfer cross-section for  0 = 5 cm 2 g −1 , and the values of  are given in the legend.The vertical coloured dashed lines correspond to different velocity scales in the system.The left most is the relative velocity dispersion within 100 kpc of the 10 15 M ⊙ cluster.The second and third from the left correspond to the relative velocity dispersion within 100 kpc around the barycentre at the first pericentre for MMR:5 and MMR:1 respectively.Finally, the second last and the last are the relative velocity of the BCGs at the first pericentre passage of the system with MMR:5 and MMR:1, respectively.

Figure 2 .
Figure 2. The dotted and solid lines correspond to the DM peak position of the main halo and subhalo respectively.Upper panel: DM peaks in equal-mass merger, lower panel: Merger with mass ratio 5. Plot labels are described in Table3.Simulation results presented in this figure corresponds to varying values of  and a fixed  0 .

Figure 3 .
Figure 3. BCG position of the subhalo vs. time.The upper panel shows peak positions in the equal-mass merger, while the lower panel corresponds to the unequal-mass merger.For better visibility, we plot only the position of the BCG of the subhalo.Plot labels are described in Table3.

Figure 4 .
Figure 4.The density of DM haloes, accounting for the particles within |  | < 100 kpc.The first column corresponds to CDM simulations, while the second and the third columns correspond to simulations with the following momentum transfer cross-sections (  0 , ) :(5, ∞),(5, 2000).The rows correspond to different times.The first row being before the pericentre passage, the second, just after the first pericentre passage and the third being at later stages.The time stamps in the images are in units of Gyr.

Figure 5 .
Figure 5. Plot similar to Figure 4, but displaying only the subhalo of MMR:5 simulation.

Figure 6 .
Figure6.Evolution of central density of an isolated halo of virial mass 10 15 M ⊙ .There are three bands in the plot corresponing to, (i) simulation with rare self-interactions with constant   of  0 = 0.175 (ii) simulation for  = 2000 km s −1 with an initial guess for  0 = 2.72 cm 2 g −1 (iii) the same simulation as the earlier one, but rescaled to  0 = 1.2 cm 2 g −1 .The width of the band corresponds to the uncertainty in the estimation of central density, and it is proportional to 1/ √  .

Figure 7 .
Figure7.The plot contains viscosity cross-section evaluated at  = 0, for different values of .At each , the inferred value of  0 is used to compute   .The triangles and stars represent   calculated using the  0 inferred from  -body simulations for fSIDM and rSIDM, respectively.The solid line corresponds to  0 inferred by solving  eff = 0.35 cm 2 g −1 .

Figure 8 .
Figure 8. BCG positions of the subhalo against time.The upper panel corresponds to equal mass merger, while the lower one corresponds to unequal mass merger.

Figure 9 .
Figure9.DM -BCG offsets in the equal-mass merger.Upper panel corresponds to fSIDM, while bottom panel to rSIDM.The plot labels are described in Table5.The vertical dashed lines correspond to the time of the first pericentre passage.

Figure A2 .Figure B1 .
Figure A2.Comparison of peaks positions of different components between low and high resolutions for CDM and SIDM simulations in the unequal mass merger.Top, middle and bottom panels correspond DM, galaxy and BCG components.The dashed lines correspond to low resolution and solid lines correspond to high resolution.The SIDM case corresponds to the frequent self-interactions with  0 = 0.5 m 2 g −1 .

Figure C1 .Figure C2 .Figure D1 .Figure D2 .
Figure C1.Evolution of the central density measured within 100 kpc around the DM peak of main halo.The parameters are CD-matched and the labels are explained in Table4.

Table 1 .
This table contains the NFW parameters used in generating the ICs for merger simulations.The first column contains the virial mass, the second the density parameter  0 :=    crit , followed by scale radius   , concentration parameter and number of DM and galaxy particles.DM-only isolated haloes have the same NFW parameters as given in the first row.

Table 2 .
This table contains the initial separation, initial relative velocity between the two clusters and the relative velocity of the BCGs at the first pericentre passage.

Table 5 .
This table contains the conservative upper bound on  0 for different values of  given in the first column.The second column contains  0 the value for frequent self-interactions, and the third for rare self-interactions.