The Signatures of Self-Interacting Dark Matter and Subhalo Disruption on Cluster Substructure

The abundance, distribution and inner structure of satellites of galaxy clusters can be sensitive probes of the properties of dark matter. We run 30 cosmological zoom-in simulations with self-interacting dark matter (SIDM), with a velocity-dependent cross-section, to study the properties of subhalos within cluster-mass hosts. We find that the abundance of subhalos that survive in the SIDM simulations are suppressed relative to their cold dark matter (CDM) counterparts. Once the population of disrupted subhalos -- which may host orphan galaxies -- are taken into account, satellite galaxy populations in CDM and SIDM models can be reconciled. However, even in this case, the inner structure of subhalos are significantly different in the two dark matter models. We study the feasibility of using the weak lensing signal from the subhalo density profiles to distinguish between the cold and self-interacting dark matter while accounting for the potential contribution of orphan galaxies. We find that the effects of self-interactions on the density profile of subhalos can appear degenerate with subhalo disruption in CDM, when orphans are accounted for. With current error bars from the Subaru Hyper Suprime-Cam Strategic Program, we find that subhalos in the outskirts of clusters (where disruption is less prevalent) can be used to constrain dark matter physics. In the future, the Vera C. Rubin Observatory Legacy Survey of Space and Time will give precise measurements of the weak lensing profile and can be used to constrain $\sigma_T/m$ at the $\sim 1$ cm$^2$ g$^{-1}$ level at $v\sim 2000$ km s$^{-1}$.


INTRODUCTION
The cold dark matter (CDM) paradigm, which treats the microphysical constituents of dark matter as collisionless, has been very successful in explaining the large scale structure of the Universe.Within this paradigm, N-body simulations have been used extensively to model the formation of nonlinear structure in the Universe (e.g., Kuhlen et al. 2012).However, certain predictions derived from N-body simulations of CDM were thought to be in tension with observations of small-scale structure in the Universe, as inferred from dwarf galaxies (Bullock & Boylan-Kolchin 2017), particularly the core-cusp issue (e.g., Gentile et al. 2004;de Blok 2010) or the missing satellite problem (e.g., Klypin bhattacharyya.37@buckeyemail.osu.eduet al. 1999;Moore et al. 1999).Self-interacting dark matter (SIDM) was originally proposed as a viable candidate to mitigate some of these small-scale problems.In its most basic form this model allows dark matter particles to undergo elastic scattering, allowing for the exchange of momentum and energy via non-gravitational processes (e.g., Spergel & Steinhardt 2000;Burkert 2000;Moore et al. 2000).SIDM models have a range of phenomenological signatures; namely allowing both for the formation of cores instead of cusps at the center of collapsed halos, and also a possible suppression of substructure, within massive hosts (Tulin & Yu 2018).
With recent advances in observational precision, combined with more sophisticated numerical techniques, especially the ability to realistically model the effects of baryonic physics on these scales, some of these problems like the missing satellite problem has mainly been attributed to observational incompleteness (Kim et al. 2018;Nadler et al. 2020b) and the effects of reionization (Gnedin 2000;Somerville 2002;Benson et al. 2002).However, to constrain the microphysics of dark matter, it remains important to study the effects of these properties, like self-interactions, on different observables from cosmological surveys.Depending on the particulars of the underlying model, these self-interactions can have a velocity-dependent differential cross section.This implies that such models are best constrained by combining information from systems which have different natural velocity scales.On one end of this spectrum are dwarf galaxies, which can probe interactions at very low relative velocities.At the other end of the spectrum are galaxy clusters inhabiting massive dark matter halos, which serve as laboratories to study some of the most energetic processes in the universe (Markevitch et al. 2004;Kravtsov & Borgani 2012).These clusters are naturally the systems best suited to study effects of selfinteractions at very high velocities.
The gravitational potential of the dark matter halos around galaxy clusters can be studied in detail, both, through gravitational lensing observations and the spatial and velocity distribution of the population of satellite galaxies within them (e.g., Natarajan & Springel 2004;Sand et al. 2004;Natarajan et al. 2007;Newman et al. 2013).Historically, galaxy clusters have been used to place constraints on dark matter cross-sections at the velocity dispersion scales set by the clusters' gravitational potential, i.e., ∼ 1000 km s −1 (e.g., Gnedin & Ostriker 2001;Miralda-Escudé 2002;Markevitch et al. 2004;Peter et al. 2013;Brinckmann et al. 2018;Harvey et al. 2019).Some of the recent constraints on σ/m are 1 cm 2 g −1 at ∼ 3200 km s −1 , as derived from merging clusters (Kim et al. 2017) or < 0.5 cm 2 g −1 at ∼ 1150 km s −1 estimated from lensing and stellar kinematics (Sagunski et al. 2021).While most early studies of SIDM focused on velocity-independent differential cross-sections, various models with explicit velocity-dependence in the differential cross-section have been proposed in order to reconcile cluster constraints with lower bounds on the cross-section from smaller galaxies.Such velocity-dependent cross sections naturally arise in models where self-interactions are mediated by light particles (e.g., Kaplinghat et al. 2016).
In a universe where dark matter has a velocity-dependent interaction cross-section, the evolution of satellite galaxies, that live in their own dark matter halos, depend both on the interaction cross-section at the high velocity scale of the cluster's velocity dispersion and also the lower velocity scale of the subhalos' own internal velocity dispersion.Young objects like galaxy clusters, where a large fraction of their most massive subhalos and their satellite galaxies have not had enough time to get tidally disrupted, provide a unique opportunity to constrain the shape of the velocity dependence of the cross-section using a single system.
In this paper, we study the population of subhalos in zoomin simulations of 30 cluster-mass objects for a velocity dependent SIDM model.We use a relatively high normalization for the cross-section so that the effects of the self-interactions are prominent compared to the noise due to scatter in the masses and other properties of the zoom-in systems.The zoom-in method in particular, allows us to robustly simulate a wide range of scales that simultaneously encompass the massive host halo and its lower mass substructures.
A host of ongoing and future surveys will provide us with large samples of galaxy clusters, allowing us to carry out statistical studies of the population of satellite galaxies in these systems.In particular, surveys like the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST) (Abell et al. 2009), Dark Energy Survey (DES) (Abbott et al. 2005) and the Hyper-Suprime Cam Subaru Strategic Program (HSC-SSP) (Aihara et al. 2018) will allow us to observe satellite galaxies that are, at least, a hundred times fainter than the Milky-Way, if not more.Meanwhile, the simulations allow us to make reliable predictions for the spatial and velocity distribution of the subhalos and their detailed internal structures.Understanding the evolution of these systems can therefore significantly boost our knowledge of dark matter microphysics in the near future (Bechtol et al. 2019;Drlica-Wagner et al. 2019).
The internal structures of satellites, in particular, can also help disentangle degeneracies between dark matter microphysics and baryonic effects.Comparing the matter distribution around galaxies with similar optical properties in clusters with those in fields can help factor out the baryonic effects on the galaxy's dark matter halo to a large extent.A widely used method to probe the structure of a galaxy's dark matter is through weak lensing, the ongoing HSC-SSP being a deep survey, is ideally suited for weak lensing studies of satellite galaxies.Using mock satellite distributions from our simulations and the covariance estimates of weak lensing measurements from the HSC survey (Kumar & More in prep), we explore the limits that a HSC-like survey can place on SIDM cross-sections using these observables.By rescaling the covariance matrix to acccount for various improvements, we also make predictions for a LSST-like survey.
Importantly, this is the first study that accounts for the "orphan" galaxy population when making predictions for observations from simulations, both in CDM and SIDM.While comparing observed galaxy distributions to subhalo distributions from dark matter only N-body simulations, one must account for subhalos that have been artificially disrupted in the simulations or failed to be tracked due to numerical effects (Springel et al. 2001;Wang et al. 2006;Campbell et al. 2018;van den Bosch et al. 2018;van den Bosch & Ogiya 2018;Delfino et al. 2021).The galaxies harbored by these subhalos, being more compact, may survive even after their corresponding subhalos disrupt.This effect has been known to bias subhalo and satellite population predictions if not properly accounted for in CDM simulations (Gao et al. 2004).In principle, these effects can be even more severe in SIDM, due to the additional evaporation of particles from subhalos due to self-interactions.Orphan modeling is therefore an important systematic to consider while forward modelling observations of satellite populations.
This paper is organized as follows.In Section 2 we discuss the details of the simulations we use.In Section 3 we describe the subhalo distributions in clusters and in Section 4 we describe the weak lensing analysis.The usage of r will imply 3D distance compared to R which will stand for the 2D projection of r.

SIMULATIONS
To simulate self-interacting dark matter (SIDM) we use the method adopted in Rocha et al. (2013), which modifies the evolution of dark matter particles in the GADGET-2 (Springel 2005) simulations, by introducing a Monte-Carlo particle scattering scheme.Banerjee et al. (2020) extrapolates this scheme to include a velocity dependent scattering cross-section.In this paper, we extend this method to run zoom-in simulations of massive galaxy clusters that have been selected from a parent box of 1 Gpc h −1 and 1024 3 particles.We resimulate 30 clusters with virial mass between (2.0 − 4.5) × 10 14 M h −1 at redshift, z = 0.The particle mass resolution in the zoom-in region is m p = 1.5 × 10 8 M h −1 .We set the search radius for self-interactions equal to the gravitational softening scale, h SI = 1 kpc h −1 .These simulations were run with the cosmological parameters set to Ω m = 0.3, Ω Λ = 0.7, A s = 2.2×10 −9 , n s = 0.96, H 0 = 70 km s −1 Mpc −1 and σ 8 = 0.85 (Banerjee et al. 2020).
In this work we focus on a velocity dependent simulation cross-section where the differential cross-section also has an angular dependence (Ibe & Yu 2010;Robertson et al. 2017).This is natural in models where dark matter interacts via dark photon mediators, or in other words a Yukawa type of potential (Kaplinghat et al. 2016), This model is parameterized by w, a characteristic velocity, below which the cross-section is isotropic as σ ∼ σ 0 but above which the cross-section not only decreases with increasing velocity (∝ v −4 ) but also becomes anisotropic, favouring scatterings by small angles (Kahlhoefer et al. 2014).We use a value of w = 1500 km s −1 in this work as this corresponds to the upper limit of velocity scale of particles in a typical cluster sized halo.The normalization σ 0 is chosen such that the momentum transfer cross-section σ T /m at v = 2000 km s −1 (the typical relative velocity between particles in the cluster) is 1 cm 2 g −1 (Robertson et al. 2017;Markevitch et al. 2004).This particular scale is chosen to approximately correspond to the constraints from the Bullet cluster which is σ/m = 2 cm 2 /g (note σ/m ∼ 2σ T ).The momentum transfer cross-section σ T /m shown in Fig. 1 , is defined as Kahlhoefer et al. (2014), The velocity dispersion scales of the subhalos and their hosts as encountered in the simulations are shown as the green and pink bands respectively.We refer the reader to Banerjee et al. (2020) for a detailed description of the simulation method.
Subhalos Hosts Figure 1.The velocity dependence of the momentum transfer crosssection σT /m for the SIDM model used in this paper.The green and pink vertical bars represent the distribution of velocity dispersion of the subhalos and their hosts respectively.Darker bands have higher numbers of objects.The relevant interaction cross-section between host particles and subhalo particles, which have relative velocities of the order of the host velocity dispersion is 2 cm 2 g −1 and decreases in the Rutherford scattering limit.The relevant cross-section for the interaction between subhalo particles is at the subhalo velocity dispersion scale.
Terminology: The dark matter halo corresponding to the largest structure of the cluster will be referred to as the "cluster host" or simply as the "host."Following the terminology adopted in Nadler et al. (2020a), smaller halos that have been accreted on to such a host will be divided into two categories: • Surviving subhalos represent halos identified by ROCKSTAR in the zoom-in simulations at z = 0. We refer to galaxies hosted by surviving subhalos as satellite galaxies.
• Disrupted subhalos are halos once identified by ROCKSTAR at z > 0 that no longer exists at z = 0 because they deposit the majority of their mass onto the host halo in the interim redshifts.We refer to galaxies hosted by disrupted subhalos as orphan galaxies.

Surviving Subhalos
We use halo catalogs generated using the ROCKSTAR halofinder (Behroozi et al. 2013a) and merger trees generated using the CONSISTENT-TREES merger algorithm throughout this work (Behroozi et al. 2013b).Here we mainly focus on the halo sample at z = 0.
We extract all particles and halos around the clusters in CDM and SIDM within a radius of 15 Mpc h −1 of the center.To study the distribution of halos around the clusters, we select halos based on their V peak , which is the peak value of the the maximum circular velocity within the halo as recorded in the merger-tree catalogs.Peak quantities like V peak are known to best correlate with luminosities of ob-served galaxies within subhalos (e.g., Wechsler & Tinker 2018;Hadzhiyska et al. 2021), as current virial masses of subhalos are often affected significantly by tidal stripping in the cluster potential (Merritt 1983;Niemiec et al. 2019).Galaxies, being more compact, are less likely to be affected by tides and their luminosity traces the original mass of the subhalo before infall (Conroy et al. 2006).
In order to select a well resolved population of subhalos, i.e., each having not less than 1500 particles when V peak is achieved, a condition of V peak > 136.5 km s −1 is placed on the CDM population of subhalos.The precise value of the cut is based on the abundance matching of SDSS galaxies (Blanton et al. 2003) and corresponds to all galaxies with rband magnitude, M r < −19.43 (see e.g., Appendix B in More et al. 2016) 1 .For the SIDM sample, we find that a lower threshold with V peak > 130 km s −1 is appropriate to match the abundance of subhalos (surviving & disrupted) in CDM.This is partly due to the enhanced disruption in SIDM and also the fact that cored halos in SIDM have smaller maximum circular velocities compared to their CDM counterparts.

Disrupted Subhalos
Since we use dark matter-only N-body simulations, we use the subhalos as proxies for galaxies in the observed clusters.This is a reasonable assumption in the dark matter paradigm considering all galaxies are expected to reside within their own halos.However, while the overall dynamics of galaxies and subhalos are expected to be similar, there are some subtleties that arise while mapping galaxies to simulated subhalos.Subhalos, being extended objects, are subject to stronger tidal forces and lose mass from their outskirts more easily compared to the galaxies within them.While tidal stripping is the primary mechanism of mass loss for infalling satellites (Merritt 1983), enhanced by effects like dynamical friction (Chandrasekhar 1949;Binney & Tremaine 1987) in CDM, subhalos in SIDM can additionally experience evaporation of particles due to self-interactions with both their own particles and host halo particles (Gnedin & Ostriker 2001;Markevitch et al. 2004).The time-integrated effect of scattering can be approximated as a net pressure-force given by ∼ ρv 2 σ/m (Dooley et al. 2016;Kummer et al. 2018) and is often referred to as SIDM ram-pressure.
While the aforesaid mechanisms constitute modes of physical disruption, artificial disruption can arise due to numerical discreteness effects (e.g., van den Bosch et al. 2018;van den Bosch & Ogiya 2018) and due to the mass resolution threshold inherent to the simulation.In particular, subhalos are not tracked once their mass passes below a resolution threshold, and the rate of mass loss preceding this may be artificially enhanced.Therefore, even if a subhalo ceases to exist in a halo catalog, this does not necessarily imply that the orphan galaxy within it gets disrupted.To alleviate these issues, orphan modeling is often necessary to understand the full dis-tribution of observed galaxies in a given dataset.In this section, we describe briefly our disrupted subhalo tracking and orphan modeling methods.

Subhalo Tracking Methodology
To track disrupted subhalos we look up the merger histories of each subhalo generated by the CONSISTENT-TREES algorithm.To find these halos, we track any subhalo that enters the virial radius of its host halo at any point in the simulation and subsequently disrupts.We only study disrupted subhalos with V peak above the threshold value of 136.5 km s −1 for CDM and 130 km s −1 for SIDM.
In this work, crucially, we use the most bound dark matter particle (MBP) of the disrupted subhalo as a tracer of the location of the associated orphan at z = 0.This is a standard choice in CDM simulation analysis because galaxies are expected to be located at the minimum of their halos' potential wells (Gao et al. 2004;De Lucia et al. 2006;Guo et al. 2011;Han et al. 2016).This choice also anticipates our weak-lensing studies in Section 4, which explore the distribution of matter around predicted galaxy locations.Thus, instead of using orbit modeling (Tollet et al. 2017;Nadler et al. 2019;Behroozi et al. 2019), where the orbit of the disrupted subhalo is predicted to infer its z = 0 position, we use actual particles associated with the disrupted subhalo to trace orphan locations.As our subhalos have ∼ 1000 particles at peak mass we do not expect that the results will change significantly if we used a different modeling method.We treat our SIDM subhalos in a similar fashion.However, we note that using a single particle as a tracer for disrupted systems in SIDM, makes it susceptible to being scattered out from the minimum of the potential well (Dooley et al. 2016) due to interactions with the particles of the host.We test the robustness of MBP as a tracer for a disrupted subhalo in Appendix A by looking at the z = 0 positions of few of the particles located around the MBP at the time when it is selected.
To isolate the MBP for each of the disrupted subhalos, we look at the snapshot corresponding to the redshift at which its V peak was attained, following which we find the particle with the minimum mechanical energy E in the reference frame fixed at the center of the subhalo.We model the potential energy, V of the particles using the functional form of that expected for an NFW halo (Binney & Tremaine 2008) for both CDM and SIDM, using the scale radius r S that has been calculated by ROCKSTAR using only the constituent particles of the subhalo: We assume that the true potential in a SIDM halo does not deviate much from our 'model' NFW potential and hence does not affect our conclusions significantly.We accept only those disrupted subhalos that have >10 particles within 0.25r 200 , with r 200 representing the virial radius of the subhalo, at the redshift of V peak .2Having identified the MBP for each of the disrupted candidates, we find their position in the snapshot at z = 0.

Orphan Contribution Model
When we trace the disrupted subhalos using their MBPs, the mode of their disruption-i.e., whether it is physical or artificial-is not specified by our model.Furthermore, the dominant mode of disruption may differ between the CDM and SIDM models.While artificial disruption due to the mass resolution limit is plausibly the dominant source of disruption in CDM simulations (e.g., Green et al. 2021), SIDM subhalos may experience both increased amounts of physical disruption due to ram-pressure stripping as well as numerical effects due to potential biases in halo-finding algorithms (Nadler et al. 2020a).In light of these uncertainties, (Pujol et al. 2017;Jiang et al. 2021), we apply a simple model to associate disrupted subhalos with orphan galaxies, the details of which are as follows: We allow a fraction of the disrupted subhalos selected using the V peak threshold to host orphan galaxies within them.Denoting the number of surviving and disrupted subhalos in a given sample as N surv and N dis respectively, we define the orphan galaxy fraction f orp as a function of a free parameter η that can be continuously varied between 0 and 1, The free parameter η encapsulates our ignorance in assigning orphan galaxies to the disrupted subhalos.It can be also thought of as the probability of a disrupted subhalo hosting an orphan (P sat,dis = η); in contrast, we assume that all surviving subhalos host satellites (P sat,surv = 1).The extreme values of η = 1 and η = 0 correspond to all and none of the candidates hosting orphans respectively.In detail, η depends on the accretion and disruption history of the disrupted progenitor, and it is interesting to explore this dependence in a future work.We emphasize that the interpretation of η depends on the mass resolution limit of the simulations.
First, we study the distribution of dark matter particles in the cluster halos in zoom-in CDM and SIDM simulations.Fig. 2 depicts the density profiles of these clusters as a function of radius normalized by the virial radii of the hosts, r 200m as determined by ROCKSTAR.We confirm on average that the mean profile shows a core within the scale radius for the SIDM clusters.This is expected since the number of scatterings per Hubble time is significantly high for particles within this radius.Furthermore, the density profile in SIDM is higher and steeper than its CDM counterpart right outside the core of the cluster (r ∼ 0.1r 200 ).Particles which get scattered to higher energies near the center of the halo end up getting transferred to larger apocentric orbits (Rocha et al. 2013).These contribute to a slight increase in density at such radial distances.We also note that the scatter in the profiles for the SIDM cores is larger than the CDM cusps.

V peak distribution of Subhalos
Fig. 3 depicts N(> V peak ), the complementary cumulative distribution of V peak of the separate and combined populations of surviving and disrupted subhalos in our simulations.The numbers of surviving subhalos are significantly suppressed in SIDM relative to CDM.As a result this gives rise to a larger number of disrupted subhalos in SIDM.When comparing the total population of subhalos we note that V peak distributions agree well with each other.However, we also observe that the total number of low-mass subhalos is slightly suppressed in SIDM, presumably as these subhalos are disrupted before they fall into the massive cluster; this can happen as low mass subhalos can often enter their current hosts as parts of groups, and group environments may have disrupted low-mass subhalos in SIDM even before infall (Nadler et al. in prep).

Radial distribution of Subhalos
Here we compare the number density of subhalos in concentric shells around the cluster center in the CDM and SIDM hosts.We evaluate the 3D number density, n(r), in twelve logarithmically spaced radial bins in units of r 200 of their hosts in the range, 0.08 < r/r 200 < 1.5.When dealing with the surviving subhalos, we do not use the subhalo tags from ROCKSTAR, instead we compute the number density of all resolved CDM (SIDM) halos above a V peak threshold of 136.5 km s −1 (130 km s −1 ) around the hosts.In Fig. 4 we show the stacked radial distribution of the surviving subhalos at z = 0 in the upper panel with the thick solid lines.The number of surviving halos is significantly reduced in the SIDM simulations.This effect persists out to nearly the virial radius of the clusters.
The distributions of the disrupted subhalos, traced by their MBPs, as a function of clustercentric radius, are shown with the dot-dashed lines in the same figure.The fractional difference between the CDM and SIDM clusters increases as we move outwards.In other words, the drop-off in the number of disrupted subhalos is steeper in CDM compared to SIDM, implying enhanced disruption in SIDM through the halo's interior compared to CDM.
In the lower panel of Fig. 4 we compare the subhalo densities to the particle densities in the simulation.We normalize the the number density at a given radius by the number density at r 200m .The dotted lines correspond to the normalized number density of dark matter particles and the shaded regions correspond to the subhalo distributions.The upper and lower envelopes of the shaded region represent the radial distribution of subhalos in the scenarios when all and none of the disrupted subhalos are taken into account.The shaded region is meant to demonstrate how adding in different fraction of the disrupted subhalos to the surviving population changes the radial distribution.If we assume orphan galaxies populate only a fraction of the disrupted subhalos, the total number density of substructure in any radial bin will be a weighted mean of the contribution from the surviving and disrupted population, 5 130) km s −1 with the surviving subhalos n sat (thick solid) and the disrupted subhalos selected using the most bound particle n orp (thick dashed).Bottom: The stacked radial distribution of dark matter particles (dashed) and satellites (shaded regions), including orphan galaxies for orphan fractions η ∈ [0, 1].All distributions are plotted as a function of clustercentric radial distances normalized with respect to the cluster virial radius r200.
where n surv and n dis are the stacked number density profiles of surviving and disrupted subhalos respectively.The shaded region corresponds to the range η ∈ {0, 1}.We observe that the full radial distribution of subhalos, including surviving and disrupted (η = 1, upper envelope), agrees quite well with the dark matter distribution both in CDM and SIDM.This is consistent with the results of Han et al. (2016); Bose et al. (2019); Green et al. (2021), and here we confirm that it holds for SIDM as well to a large extent.However, we note that the while the CDM subhalo distribution can be as steep as the dark matter distribution (in fact it can even be slightly steeper than dark matter) the SIDM subhalo profile never becomes quite as steep as the dark mat- ter profile within 0.8 Mpc h −1 .This can be partly attributed to the fact that the dark matter profile itself becomes steeper in SIDM right outside the core, as dark matter particles are pushed outside from the center region.This phenomena is specific to particles, and does not necessarily effect the profiles of subhalos.A comparison between matter and subhalo distribution can therefore be a probe of the dark matter physics in this regime.
In Fig. 5 we show the radial distribution of the redshift of disruption for subhalos in the SIDM and CDM simulations.In general, we find that subhalos tend to disrupt earlier in SIDM compared to CDM. 4. LENSING AROUND SUBHALOS Current and future lensing surveys like HSC-SSP (Aihara et al. 2019), DES (Abbott et al. 2018) and the LSST (Ivezić et al. 2019) will give us an unprecedented sample of clusters and member galaxies allowing us to measure the detailed mass distribution around subhalos in clusters using satellite galaxy-galaxy weak lensing (Li et al. 2014;Sifón et al. 2015).In this section we measure the matter distribution around surviving subhalos and evaluate the projected excess surface density profile, ∆Σ, that is the relevant weak lensing observable.We compare the stacked profiles around surviving subhalos at different projected clustercentric distances, to those around isolated centrals with the same V peak with the aim of probing the effects of dark matter self-interaction between host and subhalo particles, i.e., at the scale of the cluster velocity dispersion.To avoid the biases incurred in not accounting for the disrupted subhalos (Han et al. 2016) we also examine the mass distribution around their tracers, i.e. the most bound particles, and study the stacked profiles as a function of the orphan fraction.

Stacked 3D Density Profiles
Before we study the lensing signal around subhalos, it is instructive to look at the 3D distribution of matter around the subhalos in the simulation directly.In Fig. 6 we plot the stacked density profile of dark matter particles around surviving and disrupted subhalos along with the profiles around isolated halos of the same mass.We study these profiles for subhalos at different clustercentric distances r sub and in two bins of V peak .The solid and dashed curves depict the dark matter density around the surviving subhalos and the isolated centrals respectively.In general the core and cusp features are prominent at the centers of the respective SIDM and CDM subhalos, particularly in the massive ones (V peak > 200 km s −1 ).Two features distinguish the subhalos from the isolated halos: a major contribution from the density peak of the cluster hosts, and a relatively minor suppression in density associated with tidal stripping in the subhalos' inner structure.
The density profile around the MBP tracer for the disrupted subhalos close to the cluster center (r sub < 0.5 Mpc h −1 ) is drowned out by the core or cusp features of their respective hosts.On the other hand, the ones beyond 0.5 Mpc h −1 show signs of feeble remnant of their now disrupted cores or cusps.It is seen that their central densities are suppressed by typically 2 dex relative to the surviving subhalo sample with the same V peak .The suppression is greater in SIDM as cores are known to be more vulnerable to tidal stripping compared to cusps (Peñarrubia et al. 2010;Errani & Navarro 2021).
In the following subsections we study the lensing signal around these objects in detail in finer radial bins.

Stacked Excess Surface Density Profiles
Weak lensing measures the distorted ellipticities of background galaxies behind a lensing source.For any mass distribution the shear field is determined by the excess surface density ∆Σ of mass (Mandelbaum et al. 2005;Schneider 2005).This is connected to the azimuthally averaged tangential shear field as where Σ(R) is the azimuthally averaged projected mass density in a narrow annulus at R and Σ(< R) is the average projected mass density integrated within R. Σ crit is the critical density given by which depends on the angular diameter distances to the lens (D l ), the source galaxies (D s ), and between the lens and source (D ls ).
The value of ∆Σ can be measured directly from the dark matter particles in the simulations.We compute the ∆Σ profile around the subhalos in our simulations as a function of Figure 6.The stacked density profiles around surviving subhalos (solid) and isolated halos (dashed) with Vpeak> 136.5 km s −1 in CDM (blue) and SIDM (red) as a function of 3D radial distances from the halo centers.Besides this, the stacked density profiles around the disrupted subhalo (dotted) around the position of the MBP is shown.While the cores and cusps are present in both SIDM and CDM halos, the surviving subhalos exhibit a contribution from their hosts and tidal stripping throughout their radial extent.There are remnants of cores/cusps of disrupted subhalos is visible in the cluster outskirts.
projected radius R from the center of each subhalo.As the subhalo is embedded in the massive cluster potential there are two separate contributions to the ∆Σ profile-one from the enhanced density around the subhalo and the other from the host cluster-mass distribution at the location of the subhalo.
We measure ∆Σ in 20 logarithmically spaced bins between 0.01 − 5 Mpc h −1 centered around the subhalo centers.We count the number of particles in the 2D projected annuli or radius R around every subhalo belonging to the clusters.We assume the z direction as the line of sight direction and compute the projected quantities in the x − y plane.We project over the whole length of the simulation box.For the subhalo population comprising members from all the 30 clusters, we split them into 4 bins according to their projected clustercentric distances, R sub ∈ {0.1 − 0.3, 0.3 − 0.5, 0.5 − 0.7, 0.7 − 0.9} Mpc h −1 and study the stacked ∆Σ profile in each such bin.
We compare the stacked ∆Σ profiles of subhalos with that of centrals or isolated centrals with the same V peak threshold.These centrals were selected using the condition that they are not within 5 Mpc h −1 of any host with mass greater than 10 13 M .In a given zoom-in box, the number of centrals was found to be typically an order of magnitude larger compared to the satellites in the cluster.
In Fig. 8 we show the stacked ESD profiles around the surviving subhalos (solid lines) and isolated centrals (dashed lines) in our simulations.The top and bottom panels show the ESD profiles for two different bins of V peak .The differences between isolated centrals and subhalos as observed in the 3D density profiles are reflected here.We note that the small radii (R 0.2 Mpc h −1 ) are dominated by the (sub)halo's own overdensity and at the location of the host center (R ∼ R sub ), ∆Σ which essentially traces the slope of the density profile, changes sign.The overall amplitude of the density within surviving subhalos is suppressed compared to the isolated halos of the same V peak due to stripping of mass.This effect is significantly more severe for the lower V peak subhalos, both in CDM and SIDM.
The stripping of CDM subhalos throughout their extent is surprising; tides are expected to strip material from their outskirts, but we observe a depletion throughout subhalos' cusps.A possible explanation for this effect can be that some of the particles in the outskirts have radial or plunging orbits within the subhalo, and they get stripped when they are near the outskirts of their orbits.Moreover, subhalos can become significantly aspherical inside a cluster, which my invaldiate assumptions about their cusps not being disrupted in a spherical potential. 3.The stacked ∆Σ profiles of the disrupted subhalos (dotted lines) are very flat close to their centers because they have lost most of their mass and the projection makes the profiles appear shallower than in 3D.The net effect of adding these disrupted subhalos to the stacked ∆Σ profile of the surviving subhalos, is to suppress the average density profiles inside their virial radius.As both dark matter interactions and disruption lead to a suppression of the density profile around a galaxy, it is important to consistently account for the orphan contribution when making inferences about dark matter microphysics.
The sample of satellite galaxies that are observed in clusters may correspond in part to "orphan" galaxies in our simulations, i.e. they may exist in disrupted subhalos.Using the formalism of Section 2.2 we can identify the orphan candi-3 A peculiar feature is noticed in the innermost bin of R sub for V peak > 200 km s −1 that has CDM satellites being more dense than CDM centrals, this is due to the cluster particles themselves elevating the number density of particles near the center dates, and measure the ∆Σ profile around their MBP tracers.Because we do not know what fraction of the observed satellites are orphans a priori, we allow the fraction to vary before adding their contribution to the stacked subhalo profiles.In order to obtain a meaningful observable we appropriately combine them as a weighted average, where f orp is calculated using Eq. 4. For context, η = 1 corresponds to the case where all disrupted subhalos host orphans and η = 0 corresponds to the case where none of them do.Note that the fraction of disrupted subhalos itself varies as a function of the projected cluster-centric radius R sub as can be seen in the bottom panel of Fig. 7.
In the following section, we explore the degeneracies due to orphan modelling and self-interactions on the ∆Σ profiles, and study how well we can distinguish between them two by constructing mock lensing observables for an HSC like data set around satellite galaxies.

Mock lensing profiles for satellite galaxies
In the upper row of Fig. 9 we show the stacked ∆Σ profiles for subhalos with V peak > 195 km s −1 in SIDM and V peak > 200 km s −1 in CDM simulations when different fractions of disrupted subhalos are populated with orphans.From light to dark, the curves correspond to assigning successively higher fraction of orphans f orp .The four panels correspond to different bins of projected distances of the subhalo lens from the cluster center, R sub .We note that the disparity between the ∆Σ profiles of SIDM and CDM subhalos decrease with their proximity to the cluster center.This is because when f orp is sufficiently high, the cusp of the stacked CDM profiles is damped out and resembles cored like the SIDM profiles, confounding inferences about dark matter interactions.
In the lower panel of Fig. 9 we demonstrate a typical scenario that we will encounter while measuring the stacked weak lensing signal in cluster satellites in a universe with SIDM.The red points correspond to mock measurement of the stacked lensing profile in the SIDM simulations when 100% of orphans are assumed to have galaxies, i.e. η = 1; to this profile we add the error bars from the HSC survey.The errors are determined from a cross-correlation of the positions of satellite galaxies in SDSS redMaPPer clusters (Rykoff et al. 2016) with the shear obtained from the first year shape catalog release of the HSC-SSP (Mandelbaum et al. 2018b).These errors shown here reflect the diagonal component of the shape noise covariance C HSC which has been calculated by 320 different realizations of randomly rotated shapes of HSC galaxies around redMaPPer satellite galaxies (Kumar & More in prep).We also overlay the CDM profiles with different orphan fractions to demonstrate that given the current error bars, a typical scenario in an SIDM universe can be degenerate with a CDM model that has a high orphan fraction, particularly for satellites near the cluster center.
With the aim of studying the joint effects of selfinteractions and subhalo disruption on the subhalo ∆Σ pro-  files, we conduct a mock observation of an SIDM universe and try to fit it with a CDM model.The fraction of orphans in both the mock and model samples are varied and the goodness-of-fit is checked based on a χ 2 measure.We study the stacked ∆Σ profiles for each of the 4 bins of R sub .We construct the set of mock observations Σ SIDM from the SIDM simulations by varying the orphan fraction through 100 bins of η SIDM ∈ {0, 1} using Eq. 4. For our model space, we use the CDM simulated ∆Σ profiles Σ CDM for the same range of the parameter η CDM .For each observation in a given bin of R sub we compute the χ 2 with the model CDM profiles using the weak lensing covariance C HSC .Therefore the χ 2 at each point in our 2D space of mock and model orphan fractions is given by 9) where i and j iterate over the bins of η CDM and η SIDM .The number of degrees of freedom in our model is 19, with 20 radial bins around each lens and the free parameter, η CDM .
We show the 2D distribution of log(χ 2 /d.o.f) in Fig. 10.The x-axis and the y-axis are η SIDM and η CDM respectively.We note that in the innermost regions around the cluster, low values of η SIDM in an SIDM universe will be inferred as a high η CDM in CDM.This is because the density deficit arising due to the cored nature of the SIDM satellites can also be compensated by CDM satellites with a larger contribution of orphans.But as we move to the outer most bins there are no good fits to the assumed model for the observed set of curves.In the outer regions (R sub > 0.7 Mpc h −1 ), an orphan fraction of > 0.2 in the CDM model would not be permissible to explain the data in an SIDM cosmology, because we find that the fraction of CDM orphans can be at most 0.2 (see Fig. 7).Therefore a way to distinguish between SIDM and CDM in a HSC-like survey will be to observe the lensing profiles in the cluster outskirts.
For each value of η SIDM the minimum of the reduced chisquare χ 2 min and with it the η CDM at which the minima is obtained are plotted against each other as the solid navy blue line in the lower row of Fig. 10.The solid and dotted black lines represent χ 2 /d.o.f = 1 and the 95% confidence interval of a χ 2 distribution with d.o.f = 19 respectively.If all of the minima, χ 2 min falls outside this interval, the probability of the ∆Σ profile arising from CDM rather than SIDM substructure can be rejected at > 95% confidence independent of the underlying orphan fraction.For bins of R sub other than the outermost one, the SIDM and CDM ∆Σ profiles tend to give χ 2 min that fall within the interval.This implies that it is very difficult to constrain σ/m using the weak lensing signal from satellites at projected distances < 0.7 Mpc h −1 .However, in the outermost bin, the abundance of orphans decrease enough for their effect on ∆Σ to be insufficient in compensating for the reduced density due to self-interactions.As a result χ 2 min remains outside the interval for the full range of the best-fit η CDM which implies that this can provide a possible way to place an upper-limit on σ/m.

Projections for the LSST Survey
We also conduct the same experiment by estimating the lensing covariance for a future survey like LSST.Assuming that the magnitude of the covariance matrix for galaxy- galaxy tangential shear is inversely proportional to the survey coverage area, Ω, and the effective number density of background galaxies, n eff (Shirasaki 2015), C HSC is scaled by the appropriate ratios such that Here, we have used n eff,HSC = 21.6 arcmin −2 (Mandelbaum et al. 2018a) and n eff,LSST = 37 arcmin −2 (Chang et al. 2013).
For the sky-coverages of LSST and HSC we use Ω LSST 20000.The results are summarized in the bottom panel of Fig. 9 with the dark green solid line.Based on our estimation we note that the LSST will shrink the error bars significantly making it easy to rule out the effects of subhalo disruption in an SIDM universe when we measure the ∆Σ profiles around satellites.
In summary this demonstrates the challenges that we will face in using the weak lensing observable ∆Σ to infer the effects of self-interactions.Irrespective of the nature of dark matter, the orphan fraction should be allowed to be a free parameter when ∆Σ is measured (Kumar & More in prep).Ideally if one can constrain the orphan fraction as a function of radius, independently using the radial distribution of subhalos (Fig. 4) and then use it as a prior for the weak lensing analysis, the nature of dark matter interactions can be inferred.Alternatively, the weak lensing profiles around the satellites in the cluster outskirts may be used because in these regions the disrupted remnants are less abundant.
In the final analysis section we discuss some possible systematics that need to be considered for a realistic model.

Stellar Mass Contribution and Miscentering
Here we explore how the contribution of the baryonic mass in the satellite galaxies affects the weak lensing profiles and therefore our inferences.Stellar masses M are assigned to the subhalos using a M -V peak relation from Campbell et al. (2018), from which we derive ∆Σ profile by assuming the galaxy to be a centrally located point mass.The details of our method are described in Appendix B. Furthermore, we also test the effects of cluster miscentering (Mandelbaum et al. 2008) on the subhalo profiles.Miscentering may be significant source of systematic uncertainty as it will cause subhalo distances from the cluster center to be mislabelled.Our method of introducing miscentering in the simulated ∆Σ profiles is described in Appendix C.
The impact of including stellar mass and miscentering on the χ 2 min estimate is depicted by the dotted and dashed lines respectively in the lower panel of Fig. 10.Although the addition of a stellar component may appear to drastically change the dark matter only ∆Σ profile compared to the feeble effect of miscentering (see Fig. 12), these two have opposite effects on the χ 2 estimate.This is because the covariance for the inner radial bins (R < 0.01 Mpc h −1 ) is much larger compared to the radial bins at R ∼ 1 Mpc h −1 , e.g., the errorbars reflect this in Fig. 9. Therefore, systematics at large scale like miscentering contribute to the χ 2 more than bary- onic effects at small-scales.Nonetheless, Fig. 10 show that, satellites with R sub >0.7 Mpc h −1 could be still used to probe SIDM as miscentering leads to poorer rather than better fits to CDM.Our miscentering estimates are derived using the miscentering fractions for the RedMaPPer cluster finder (Melchior et al. 2017).Optical cluster finders like RedMaPPer, that assign centers to bright central galaxies have larger fractions of miscentered clusters than X-Ray or SZ selected clusters (Zhang et al. 2019), therefore the effects of miscentering can be mitigated by using alternative cluster finders.

DISCUSSION & OUTLOOK
In this section we highlight key takeaways from our analysis and some caveats that must be accounted for to make robust inferences for dark matter physics based on comparisons of cluster satellites to cosmological simulations.

Artificial Disruption and Satellite Abundances
While van den Bosch et al. (2018); van den Bosch & Ogiya (2018); Errani & Navarro (2021) discuss forms of numerical disruptions inherent in CDM simulations, an equivalent study with the same detailed analysis for SIDM is yet to occur.By focusing on satellites hosted by relatively well-resolved subhalos (for example, we only studied SIDM subhalos with V peak > 195 km s −1 for the lensing analysis), we have mitigated the impact of artificial disruption; however, this effect may become severe for less massive substructure in SIDM due to evaporation (Gnedin & Ostriker 2001).
In Rocha et al. (2013) the subhalo counts for cluster hosts are suppressed in a σ/m = 1 cm 2 g −1 realization of SIDM by a few percent relative to the CDM equivalent, especially in the inner region of the halo (r < 0.5r vir ).Our results for surviving subhalos seem consistent with this result, however when disrupted subhalos are taken into account, both SIDM and CDM subhalos are found to be equally abundant within clusters.Therefore, SIDM is comparable with CDM in being able to explain the abundance of massive substructure in clusters (Moore et al. 1999;Natarajan et al. 2017).This is expected if all halos that form in CDM also form in SIDM, given that the primordial matter power spectrum is the same for both models down to the mass scales we are observing in this paper (Vogelsberger et al. 2016;Huo et al. 2018).However, there may still be small differences in the SIDM subhalo abundance in clusters from disruption within larger groups before infall into a larger host halo (Nadler et al. in prep).
To use the abundance of satellites in clusters for constraining models of dark matter physics, whether using studies of the spatial distribution of bright satellites (Budzynski et al. 2012;Shin et al. 2021) or lensing mass maps (Natarajan et al. 2017), it is therefore imperative to understand how galaxies populate their subhalos and how subhalo disruption is related to the disruption of the galaxy within them.While central dark matter cusps are never completely disrupted in CDM (Errani & Navarro 2021), the existence of cores in SIDM subhalos makes these systems more susceptible to disruption (Peñarrubia et al. 2010).In addition, the central cores may lead to more extended satellite galaxies with shorter tidal disruption timescales compared to CDM.A detailed study of this effect is required to interpret observed satellite populations in the context of SIDM.

Baryonic Contribution and Galaxy-halo Connection
In principle, to get a complete picture of the evolution of satellite galaxies in clusters, we need to robustly account for their baryonic component and its influence on dark matter (e.g., Gnedin 2000;Benson et al. 2002;Brooks & Zolotov 2014;Schaller et al. 2015;Sawala et al. 2016;Kim et al. 2018;Nadler et al. 2018).In this paper we model the baryonic contribution to the ∆Σ profiles in an ad hoc manner (see Appendix B) and therefore do not account for co-evolution of the galaxy's dark matter halo and its baryonic component.
There are two primary ways in which baryons affect their halos.Firstly, adiabatic contraction of dark matter orbits in the presence of a central galaxy can make the centers of halos appear cuspy (Gnedin et al. 2004); it has been shown that this effect is enhanced in SIDM halos making them appear as cuspy as CDM halos, if not more (Despali et al. 2019;Sameie et al. 2021;Kaplinghat et al. 2014).Secondly, feedback processes from supernovae (Pontzen & Governato 2012) and AGN (Peirani et al. 2017) remove dark matter particles from the center creating a core.Both these competing effects can, in principle, make the density profiles of subhalos in SIDM appear degenerate with a scenario that has both CDM and baryons, particularly within 0.01 Mpc h −1 .However, this is unlikely to change our inferences from the stacked ∆Σ profiles as most of our signal comes from the enhanced stripping of SIDM subhalos due to interaction between host and subhalo particles, which impacts the profiles throughout the subhalos entire extent.
Several studies have been conducted on the evolution of Milky Way satellites in the presence of SIDM (Fry et al. 2015;Robles et al. 2017Robles et al. , 2019;;Fitts et al. 2018;Sameie et al. 2018;Elbert et al. 2018;Dooley et al. 2016).On cluster mass scales Robertson et al. (2017) and Robertson et al. (2021) have studied the evolution of cluster in the presence of baryons, but have focussed mostly on host profiles.We note that the presence of a massive central galaxy in the cluster can also affect the survivability of satellites that are on highly radial orbits-this effect is more severe in disk centrals due to its axis-symmetric potential.On Milky Way scales, it has been shown that the disruption of satellites due to the central disk can be significant in shaping the abundance of satellites and their radial distribution.On cluster scales, however, this effect is less severe considering cluster centrals are mostly elliptical and generally the bright satellite galaxies have been on fewer orbits within a cluster compared to satellites of the Milky Way.
Finally, it has recently been pointed out that the timescale for gravothermal core-collapse can be shortened for satellite galaxies at large interaction cross-sections ( 10 cm 2 /g), due to tidal stripping of their outer profiles (Nishikawa et al. 2019;Kaplinghat et al. 2019;Correa 2021).This can further accelerate in the presence of baryons that can aid adiabatic contraction and generally the satellite survival probability (Haggar et al. 2021).However, we note that the cross-section at the satellite mass-scales investigated in this paper are unlikely to be much higher than a few cm 2 g −1 , making these systems unlikely candidates for core-collapse (Kaplinghat et al. 2016).Moreover, because clusters are relatively young objects, a large fraction of satellites are unlikely to have been within the cluster potential for longer than a few orbital timescales, making them less susceptible to enter the collapse phase.
Given these subtleties and open questions, it is important to explore the galaxy-halo connection in an SIDM universe to account for the interplay of baryonic and dark matter physics in these models.Our work is a step towards disentangling some of the degeneracies in the domain of dark matter-only simulations.Our findings imply that while the radial distribution of subhalos, including both surviving and disrupted systems can be quite similar in CDM and SIDM, the inner structure of these objects can be significantly different.Furthermore, these differences can significantly affect the stellar properties of galaxies that live within them, and vice versa (Dooley et al. 2016).

Predictive Power of Weak Lensing
Measurements of the weak lensing profile around satellite galaxies have previously been used to infer the dark matter distribution around them (Sifón et al. 2015(Sifón et al. , 2018)).The shape of the weak lensing profile around the satellites is a sensitive probe of the velocity dependent interaction cross-section as it depends both on the cross-section at the low mass subhalo scale and the cluster mass scale.However, we show that inferences about dark matter can be complicated by degeneracies between the galaxy occupation in disrupted subhalos and the SIDM model.The contribution to the total weak lensing profile around satellite galaxies that comes from disrupted subhalos, which have dramatically stripped dark matter profiles needs to be accounted for to accurately infer the effect of particle interactions.
To forward model the weak lensing signal around satellite galaxies from a given simulation, one can first attempt to constrain the contribution from disrupted subhalos ( orphan galaxies ) by measuring the number of observed galaxies as a function of radius.This will allow us to use the orphan fraction as a prior to estimate differences from CDM subhalo profiles.Alternatively, the weak lensing signal can be measured around satellite galaxies near the outskirts of the cluster (R sub > 0.5r 200 ) where there are fewer disrupted systems (Fig 7).However, we expect the exact value of the inner radial cut-off region to depend on sample of galaxies that we observe.
In this work, we have selected subhalo samples with V peak thresholds to match the total number of subhalos in CDM and SIDM.This choice was made envisioning a measurement where the satellite galaxy sample is selected based on a luminosity threshold in observations and its corresponding subhalo sample is selected based on abundance matching.We set the abundance using subhalos from CDM with a V peak threshold that is matched to the abundance of typical galaxies used for cluster studies in galaxy surveys.However, as noted previously, we find that the same abundance is obtained in SIDM at a smaller V peak (see Sec. 2), implying that the SMHM relation in SIDM can potentially be different.Alternatively, we can also obtain the appropriate V peak selection by directly measuring weak lensing profiles around isolated galaxies at a given luminosity and use it for our simulation comparisons between CDM and SIDM.We find that using the latter method, as expected, the V peak inferred from the stacked ∆Σ profiles of SIDM, isolated subhalos is smaller than CDM (due to coring of the profiles).When subhalos are selected using the V peak inferred from isolated profiles, we find that the degeneracies between SIDM and CDM become generally more severe.However, we find that our overall inferences do not change significantly, i.e. subhalos in the cluster outskirts are still the most promising probes of dark matter and observations from the LSST survey should capture the subtle changes throughout the subhalo profiles and help constrain deviations from CDM.

RESULT & CONCLUSION
We run 30 N-body, zoom-in simulations of cluster-mass (>10 14 Mpc h −1 ) dark matter halos with a velocity dependent SIDM cross-section to conduct a detailed comparison of the distribution and properties of the massive subhalo population with peak velocity V peak 130 km s −1 .Throughout, we have aimed to consistently account for the population of disrupted subhalos in the dark matter simulations and understand their impact on the statistics of various observables that can potentially help understand the nature of dark matter.In particular, we have focused on the subhalo radial distribution and the weak-lensing profile around subhalos in observed galaxy clusters.The principle findings of our work are as follows, -The potential remnants of disrupted subhalo and satellite galaxies must be accounted for in order to generate robust predictions for subhalo and satellite populations from SIDM simulations.
-The V peak and radial distributions of the combined sample of disrupted and surviving subhalos in the CDM and SIDM scenarios agree with each other reasonably well.
-While the radial number density profile of subhalos can be as steep as the dark matter density profile in CDM, the subhalo profile remains shallower than dark matter in SIDM even when disrupted systems are accounted for.
-The coring and enhanced stripping of subhalos prevalent in SIDM can be degenerate with the parameters that control the galaxy occupation of disrupted subhalos, e.g., the effects of surviving subhalos' cored profiles in SIDM can be mimicked by a CDM model with an enhanced orphan fraction.
-The degeneracy between coring in SIDM and orphan modelling can be broken by studying the weak lensing signal around satellite galaxies in cluster outskirts (> 0.5r 200 ), where disrupted subhalos are rarer, especially in CDM.
-Given state-of-the-art weak lensing covariances, largescale systematics like cluster miscentering are important compared to baryonic effects at the galaxy center.Nonetheless, we forecast that LSST will able to constrain σ/m at the ∼ 1 cm 2 g −1 level from satellite galaxy-galaxy weak lensing measurements.
N-body simulations have helped us study the evolution substructure in massive halos to a great detail.The dynamics of subhalos in clusters can potentially be complicated; the precise orbit of the subhalo determining the evolution of its internal structure through time.We find that the density profiles of subhalos, particularly those that have been disrupted, can be significantly altered from their isolated counterparts and we attempt for the first time in this work to consistently incorporate the effect of these disrupted systems on the observables both in CDM and SIDM.We also note that while the central cusp in disrupted systems survives in CDM, the central regions of SIDM subhalos do not (Appendix A); this poses a challenge to assign galaxies to SIDM subhalos by traditional methods using bound particles and in future it will be important to explore alternative methods.As massive clusters are relatively easy to observe, the rich diversity of substructure allows us to study hierarchical structure formation and understand the nature of dark matter.Developing semi-analytical treatments of galaxy evolution using N-body simulations of SIDM is therefore essential to exploit the large statistical samples that will be available to us in the near future.
While in this work we have focussed primarily on bright satellites of cluster halos, in principle many of our inferences can be extended to lower mass subhalos.Recent observations help us probe the lower end of the halo mass function allowing us to study fainter systems like ultra-diffuse galaxies (UDGs; van Dokkum et al. 2015;Koda et al. 2015;Mihos et al. 2015;Tanoglidis et al. 2021.Intriguing deviations from CDM have been pointed out in (Meneghetti et al. 2020) where they find a potential excess of small-scale lenses in strong lensing studies, while we do not see an excess of substructure in SIDM systems in the mass range that we have explored, at lower mass scales, core collapse can potentially make substructure robust to disruption near the central region of clusters.Extending SIDM zoom-in simulations into the dwarf galaxy regime will therefore be important to study substructure in the low surface brightness regime.
Future surveys that probe much larger volumes will give us the opportunity to probe thousands of galaxy clusters significantly reducing statistical uncertainty in lensing measurements.eROSITA (Pillepich et al. 2012) will provide X-ray samples that will allow us to study objects that extend down to group mass, surveys like CMB S4 (Abazajian et al. 2016) and Simons Observatory (Ade et al. 2019) will add to the already existing catalog of Sunyaev-Zeldovich (SZ) selected clusters from Planck (Ade et al. 2014), Atacama Cosmology Telescope (Hilton et al. 2021) and the South Pole Telescope (Williamson et al. 2011).Moreover, the sample of optical clusters will also largely increase with the advent of the Rubin Observatory.While accurate substructure modeling is an ongoing, open problem, our set of zoom-in cluster simulations has enabled us to develop an understanding of the differences in the distributions in surviving and disrupted subhalos when self-interactions are introduced.A further understanding of the co-evolution of baryons and dark matter in SIDM simulations to constrain the detailed galaxy-halo connection should lead to more precise predictions, which, coupled with future surveys, can be used to elucidate the nature of dark matter.According to Campbell et al. (2018) a good fit is provided by, log(M 0 ) = 9.95 ± 0.01, log(V 0 ) = 2.177 ± 0.005, α = −5.9± 0.1 and β = −0.25 ± 0.02 .The contribution of stars to the stacked subhalo ∆Σ profile ∆Σ is obtained by assuming that the distribution of the stars can be treated as a centrally located point with respect to the extent of the dark matter halo.
The mean value of M is calculated for the relevant bin of R sub and V peak and the ∆Σ thus derived is added to the existing dark matter only ∆Σ profile.∆Σ(R, R sub ) = ∆Σ sub (R) + ∆Σ host (R, R sub ) + ∆Σ (R) (B3) C. MISCENTERING The positions of the centers of the simulated halos are determined in ROCKSTAR as the mean position of a confined set of particles around the density peak of the halo (Behroozi et al. 2013a).We artificially miscenter the ∆Σ profiles we obtain from the simulations in a Monte Carlo fashion to ascertain its effect on the mock lensing observable analysis.
We apply the method of Melchior et al. (2017), assuming that a fraction f mis of the cluster centers are miscentered and the stacked ∆Σ profile around their subhalos is ∆Σ mis .Similarly the fraction of clusters that are well centered have a stacked subhalo ∆Σ profile ∆Σ 0 , ∆Σ = (1 − f mis ) ∆Σ 0 + f mis ∆Σ mis (C4) The miscentered clusters are chosen by randomly sampling a fraction f mis of all the simulated clusters.The profile of each surviving and disrupted subhalo belonging to each of them is miscentered by reassigning a new value of R sub to them.We obtain the miscentered profiles by stacking in bins of R sub,mis , the projected distances between the host center and subhalo centers that have been modified by using the cosine law, Here we assume that the uncertainties on the position of the subhalo centers are negligible compared to the miscentering offset of their hosts.A value of the offset radius R off is sampled from a Rayleigh distribution (Johnston et al. 2007) that is parameterized with σ off which represents the mode of the distribution.
We choose a value of σ off = 0.2 Mpc h −1 and f mis = 0.22.Likewise, for each subhalo, cosθ is drawn from a uniform distribution in the interval (−1, 1)

Figure 2 .Figure 3 .
Figure2.The dark matter density profile of the cluster hosts (faint dotted) and their mean (dashed).Prominent core and cusp features are visible at the centers of the SIDM and CDM hosts, respectively.

Figure 4 .
Figure 4. Top: The stacked radial distribution of CDM (SIDM) subhalos with Vpeak > 136.5 (130) km s −1 with the surviving subhalos n sat (thick solid) and the disrupted subhalos selected using the most bound particle n orp (thick dashed).Bottom: The stacked radial distribution of dark matter particles (dashed) and satellites (shaded regions), including orphan galaxies for orphan fractions η ∈ [0, 1].All distributions are plotted as a function of clustercentric radial distances normalized with respect to the cluster virial radius r200.

Figure 5 .
Figure5.The dependence of zdisrupt, the scale at which the disrupted subhalos were last resolved as halos, on the z = 0 clustercentric distances of the orphan tracers.Subhalos tend to disrupt earlier in SIDM compared to CDM.In both models, the oldest disrupted subhalos are concentrated near the host center.

Figure 7 .
Figure 7. Top: The histograms of Rsub, the projected clustercentric distances of satellites (solid) and orphans (dotted).Bottom: The maximum possible value of the orphan fraction forp, i.e., with η = 1, as a function of Rsub.The vertical black dotted line represents the average value of 0.5r200 of the cluster hosts.

Figure 8 .
Figure 8.The ∆Σ profiles for the two bins of Vpeak corresponding to the upper and lower rows and the four columns from the left to right stand for the increasing bins of Rsub.The solid and dashed lines represents the signal from surviving subhalos and isolated centrals respectively, with the shaded band implying the 1σ uncertainty.The faint dotted curves are the ∆Σ profile around the positions of the MBP tracers representing the disrupted subhalos.

Figure 9 .
Figure 9.The stacked ∆Σ profiles for satellites+orphans with Vpeak > 200 km s −1 , created using Eq. 8 at four different distances from the cluster center (in Mpc h −1 ).The upper panel depicts the effect of varying η for the SIDM and CDM models (colorbar on extreme right), with a darker shades implying a larger orphan fractions.The bottom panel shows a mock measurement created from SIDM satellites and orphans with η = 1.The errorbars are derived from the data covariance matrices of the weak lensing shear measured in the HSC-SSP survey (Kumar & More in prep).The χ 2 red from fitting the mock profiles with the CDM model profiles is shown in the inset plots in the upper corners.

Figure 10 .
Figure10.Upper row: The logarithmic reduced χ 2 obtained by fitting the mock ∆Σ profiles around SIDM subhalos using HSC error bars, with CDM profiles by varying the orphan contribution, η.The x-axis shows the different values of ηSIDM used to create the observable, the y-axis shows the model ηCDM.The loci of minima for each ηSIDM is plotted as a dark blue line.The interplay of the effects of inclusion of orphans and self-interactions makes it difficult to disentangle the two for Rsub < 0.7 Mpc h −1 .Lower row: The values of the minimum χ 2 /d.o.f. are plotted against the position of the minima, ηCDM.The horizontal black solid line represents χ 2 /d.o.f.= 1 and the horizontal black dotted line is the 2σ upper limit for this χ 2 /d.o.f.distribution.A larger value of the minimum χ 2 /d.o.f., implies a poor fit of the CDM model to the SIDM mock observation, thereby represents a possible way to constrain a SIDM model.The dotted lines shows the result of incorporating a stellar mass contribution to the ∆Σ profile and the dashed lines reflect the effect of miscentering the host.The dark green lines shows the expected constraining power of the LSST survey.

Figure 12 .
Figure12.The dark matter only profile with an orphan fraction forp = 0.38 is depicted by the solid line and the addition of the stellar component is shown with a dotted line.When the hosts are miscentered using forp = 0.22 and σoff = 0.2 Mpc h −1 followed by stacking, the result is the dashed line.