Cluster halo shapes in CDM and SIDM models: Unveiling the DM particle nature using a weak lensing approach.

Self-interacting dark matter (SIDM) is an alternative to the standard collisionless cold dark matter model (CDM), allowing for interactions between the dark matter particles through the introduction of a self-scattering cross-section. However, the observable effects between these two scenarios are hard to detect. In this work we present a detailed analysis of an application of galaxy-galaxy lensing to measure with high precision the shapes of cluster halos and how this approach can be used to obtain information regarding the nature of the dark matter particle. Using two sets of simulated data, SIDM and CDM simulations, we compute stacked shear maps centred on several subsets of halos with masses ≳ 10 13 . 5 𝑀 ⊙ . From these maps, we obtain the quadrupole profiles related to the mean projected elongation of the particle distribution from which the shape parameters are derived. Accounting for a radial shape variation, this technique provides an enhancement of the observed differences between the simulated data-sets. In particular, we obtain a higher slope of the power law for the shape-radial relation for the halos identified in the SIDM simulation, which are rounder towards the centre. Also, as approaching to the mean virial radius, the projected semi-axis ratios converge to similar values than in the CDM simulation. Moreover, we account for the impact of the neighbouring mass, where more strongly elongated distributions are found for the halos in the SIDM simulation, indicating that under dark matter self interaction, the large scale structure imprints a more coherent accretion process.


INTRODUCTION
According to the current cosmological paradigm, ΛCDM, most of the matter content in the Universe is in the form of cold dark matter (CDM), i.e. non-relativistic and collisionless dark-matter (DM) particles.The success of this model is evidenced in accounting for most of the current cosmological observations with great precision (e.g.Weinberg et al. 2013;Alam et al. 2017;Planck Collaboration et al. 2020).However, after more than five decades after measuring the galactic rotation curves (Rubin & Ford 1970), which is one of the most substantial indicators of this kind of matter, very little is known regarding the nature of these particles.In this context, alternative models that postulate other particle characteristics in addition to those presented by the standard model, become relevant and are interesting subjects to be explored.A strategy to address the lack of knowledge of the DM particle nature is to explore the ranges in which the current cosmological model seems to be less accurate.In this sense, current observations point out at halo scales, as those ★ also at Port d'Informació Científica (PIC), Campus UAB, C. Albareda s/n, 08193 Barcelona (Barcelona), Spain.
Self-interacting dark matter (SIDM) is a theoretical candidate of dark matter that allows strong self interactions, i.e. of an amplitude comparable to the strong force.It was proposed as an explanation for certain observed astrophysical phenomena that cannot be explained by traditional cold dark matter models (Spergel & Steinhardt 2000), like the shallowed core of galaxies (Flores & Primack 1994;Simon et al. 2005;Walker & Peñarrubia 2011) and clusters (Sand et al. 2004;Newman et al. 2009;Umetsu et al. 2012).Phenomenological, SIDM is an attractive alternative to CDM since it modifies the density distribution at halo scales while conserving the distribution at larger radii, thus leaving intact the successful description of the CDM on large scales.An extensive review of the SIDM model and its predictions can be found in Tulin & Yu (2018).
One of the predictions of the SIDM is the modification of the halo shapes.Within the CDM paradigm halos are triaxial structures, while when considering collisions, the orbits of dark matter particles tend to become more isotropic, resulting in an overall rounder mass distribution.The impact of the self-interactions in shaping the halos has been proposed as a test for determining the DM cross-section (Miralda-Escudé 2002; Peter et al. 2013;Brinckmann et al. 2018;Despali et al. 2022).However, assessing to the halo shapes through observational techniques is challenging, mainly because of projection effects.A widely employed method for constraining halo shapes is strong lensing.Nevertheless, this approach is mostly limited to the inner radial regions, where the impact of baryons is significant and may lead to potential confusion with the imprints attributed to the SIDM.
In this scenario, weak-lensing offers a useful approach, since this effect is sensitive to the halo outskirts and can be used to constrain the total halo shape.In particular, weak-lensing measurements have been successfully used to measure the halo shapes of galaxy clusters, both using observations and simulated data (Evans & Bridle 2009;Oguri et al. 2010;Clampitt & Jain 2016;van Uitert et al. 2017;Shin et al. 2018).However, given the lower densities in these regions, the imprints of SIDM are expected to be almost indistinguishable from the CDM predictions.This fact plus the projection effects, make this methodology almost useless for distinguishing between the two proposed scenarios (Robertson et al. 2023).
In this work, we propose and evaluate using weak-lensing techniques to determine the projected shapes of galaxy cluster halos, with the aim of distinguishing between the predictions of the CDM and SIDM models.In particular, we propose to apply weak-lensing stacking techniques, also known as galaxy-galaxy lensing, to estimate the halo shapes.We expect that the combination of the stacking power with the projection of the lensing signal along the main-axis direction of the halo, can enhance the signal-to-noise ratio, thereby enabling a more precise shape estimation.In order to test the proposed methodology, we produce two sets of simulated data with the same initial conditions: one following the CDM scenario and the other accounting for dark matter self-interaction.Then, we create lensing maps from both simulations by taking into account the projected density distributions of the particles stacked at the centre of each halo.From these maps we compute the lensing radial profiles and fit them to obtain the aligned mean projected elongation of the stacked halos.Finally, we compare the measurements between the two simulations.
The work is organised as follows, in Sec. 2 we describe the produced simulations and how the halos are identified.We detail the shape computation using the particle distribution and the halo characterisation in Sec. 3. The formalism of the lensing analysis is presented in Sec. 4. We present the results of the maps produced and the lensing measurements in Sec. 5. Finally, we discuss our findings and conclude in Sec.6 and 7.

Simulated data
Multiple analyses of massive groups have constrained the crosssection of the dark matter particle to values of / ≲ 1.25, cm 2 /g (see, for instance: Randall et al. 2008;Jee et al. 2014;Wittman et al. 2018).Taking this into account, we produce two simulations with the same initial conditions and two different dark matter particle cross-sections, / = 0  2 / and / = 1  2 /, named CDM and SIDM simulations, respectively.Each simulation consists of two combined boxes with 120 ℎ −1 Mpc on side and 1024 3 particles, reaching a resolution of ∼ 10 8 ℎ −1 M ⊙ .The cosmological parameters used are Ω  = 0.3, Ω Λ = 0.7 and  0 = 70 km s/Mpc.The initial conditions were generated with the MUSIC code (Hahn & Abel 2011) at  = 50 and the simulations were run with the GIZMO code (Hopkins 2015).
The implementations of elastics self-interactions in GIZMO are described in detail in Rocha et al. (2013).In a nutshell, the probability of scattering in a time step is calculated as: (1) In the equation, / is the transfer cross section per unit mass,   the mass of the macroparticle,    is the relative velocity between these particles and    is a number density factor that accounts for the overlap of two particles smoothing kernels.For each pair of particles with an interaction probability greater than zero, the interaction is determined by drawing a random number.If interactions occur, the elastic scattering determines the exit velocity.These velocities are calculated as a function of the masses, the relative velocity and the velocity of the mass centre of the two particles.We use the isotropic scattering implementation, so the direction of scattering is chosen randomly.(Meskhidze et al. 2022).

Halo identification and properties
From both simulated data-sets snapshots at  = 0, CDM and SIDM, we identify dark matter halos using the rockstar phase-space halo finder code (Behroozi et al. 2013).In particular, for our analysis, we consider only host halos, i.e. halos not contained within another larger halo.Through this work, we use several halo properties calculated by rockstar code.All of these calculations are performed on the complete set of particles classified as gravitationally bound within each respective halo.The viral radius, denoted as  vir , is determined through a profile fitting technique and the viral mass,  vir , corresponds to the mass enclosed within this radius.The position of each halo is established using the innermost particles of the halo, while the halo-offset, referred to as  off , is defined as the spatial separation between the halo location and its centre of mass.This parameter can be used to establish the halo relaxation state (3.1).
In this study, we investigate halos with masses 10 13.5 ℎ −1  ⊙ <  vir < 10 15 ℎ −1  ⊙ .Taking into account the position centres of the selected halos in the CDM simulation, we select those in the SIDM simulation with a distance difference between the centres lower than 0.5  vir and differences between the masses up to 20%, obtaining in total 501 halos in each simulation.This choice was adopted after varying the distance difference and checking the correlation between the halos according to the measured masses.With the adopted cuts, roughly 97% of the halos identified in the considered mass range, have their counterpart in both simulated data-sets.The proximity and mass similarity criteria used in our matching process, ensure robust comparisons between the two simulated scenarios We are also interested in testing the variations of the halo shapes according to their local environment.With this aim, we characterise the environment of each halo by computing the distance to the 5 th neighbour in each simulation,  5 , considering as neighbours those halos with  vir > 10 13 ℎ −1  ⊙ .This parameter is a common observable proxy to characterise the local density of the environment (e.g.Lackner & Gunn 2012;Muldrew et al. 2012;Ching et al. 2017;Santucci et al. 2023).We inspect the variations of this parameter with the halo shapes in 3.1.
In this section, we present the characterisation of halo shapes based on the particle distribution.Additionally, we explore the connection between the observed differences in shape within both simulated data-sets and halo properties, while also defining the subsets of halos to be used in the stacking procedure.

Individual halo shape estimates
In order to perform the stacking procedure, it is necessary to first compute the main orientations of the halos before their combination.With this aim, we characterise the shape of the particle distribution for each halo using only the identified bound particles given by rockstar.Since we are interested in the projected shapes, we obtain the projected particle positions by defining a tangential plane perpendicular to the line-of-sight vector pointing towards each halo centre.We consider all three primary axes of each simulation for the projection, thereby increasing the number of stacked halos.
We calculate the orientation and determine the 3D and projected semi-axis ratios of the particle distribution according to the eigenvectors and eigenvalues of the shape tensors.We use a common definition of the shape tensor, the reduced tensor: where the summation runs over the   bound particles,   ,   are the coordinates of the −particle,   is the distance to the halo centre and the sub-indexes ,  ∈ [0, 1, 2] and ,  ∈ [0, 1] for the 3D and projected positions, respectively.We choose to use the reduced tensor definition since it accounts for the radial distance, being more sensitive to the inner particle distribution than the other commonly used definition, the standard one (in which   is neglected, Eq. 5 in Brinckmann et al. 2018).Given that the particle density is higher towards the halo centre, the orbits of these particles are more affected by differences in the cross-section.Therefore, a shape definition based on this tensor highlights the differences between the considered simulations.

HALO CHARACTERISATION
By diagonalising the tensors, we obtain the semi-axis values ( >  >  in 3D and  >  in projection) corresponding to the square root of the eigenvalues, while their respective eigenvectors ( ê , ê and ê in 3D and ê , ê in projection) define the principal axis directions associated with their eigenvalues.From these quantities, we can obtain the semi-axis ratio of the projected particle distribution  = / and the position angle of the main semi-axis, , computed from ê .We also obtain for the 3D inertia tensor the sphericity  = /, which would be used to inspect the differences in the elongation related to the halo properties (see 3.1).Additionally, we consider an alternative definition of the shape (see e.g.Katz 1991;Dubinski & Carlberg 1991;Zemp et al. 2011) calculated by using all the particles contained within an ellipsoidal radius ( ell ) given by: ell = .ê ,  ell = .ê and  ell = .ê , where  is the position vector of the particle in 3D with respect to the halo centre.Analogously, in 2D, we take into account the particles within an ellipse with radius: The procedure for calculating this shape is iterative because the value of the ellipsoidal radius is not known in advance and depends on the semi-axes of the ellipsoid, which need to be determined during the process.In practice, we first consider all the bound particles and calculate the initial set of eigenvectors and eigenvalues.Then, we compute the new set of shape parameters excluding those particles outside the obtained ellipsoidal/ellipse region.This calculation is repeated considering 10 iterations or until the difference between the newly and the previously computed shape parameter ( in 3D or  in projection) represents less than the 1%.We obtain the shape parameters using both iterative and noniterative approaches.In general, the semi-axis ratios obtained iteratively exhibit a strong correlation with those calculated using the non-iterative approach, with the latter being approximately ∼ 20% lower.For the stacking procedure, we use the orientations given by the non-iterative method.This last choice is motivated by the fact that in observational studies the distribution of the galaxy members is commonly used as tracers of the underlying dark matter distribution, in order to estimate the cluster orientations (Huang et al. 2016;van Uitert et al. 2017;Shin et al. 2018;Gonzalez et al. 2021a).Therefore, due to the typically small number of identified members (≲ 100), discarding the tracers in each iteration can introduce bias into the shape parameters.This is because the introduced shot noise tends to predict more elongated shapes.
On the other hand, to characterise the individual halo shapes based on the dark matter particle distribution, we use the semi-axis ratios obtained from the iterative method.This is done since this methodology is more commonly used given that leads to less biased results (Behroozi et al. 2013;Robertson et al. 2023).Shape definition is a matter of importance specially for this kind of analysis, since halo shapes are particularly difficult to estimate in cored density profiles that get monotonically more spherical towards the centre, which are expected mainly in the SIDM scenario.A detailed discussion on several shape definitions can be found in Zemp et al. (2011) and in particular for the comparison between SIDM and CDM halo shapes in Peter et al. (2013); Brinckmann et al. (2018).It is important to highlight that these shape parameters are only used for comparative purposes and do not impact the lensing determinations, which constitute the focus of this work.

Halo subset definitions
The stacking procedure is commonly applied by combining the halos that share a common property, such as a mass proxy.In order to take advantage of this combination technique, we initially investigate how the nature of the particles influences the shaping of cluster halos based on certain intrinsic properties.The properties inspected are the halo relaxation state according to the scaled central position offset,  off / vir , the halo mass,  vir and the halo environment according to  5 .We characterise the shape differences between the two simulated data sets using the computed 3D semi-axis ratios and the sphericity  according to the iterative method presented in the last section.
In Fig. 1 we show how the ratio between the shape estimates in the SIDM and CDM simulations ( CDM / SIDM − 1), are related with the considered halo properties.As it can be noticed, relaxed halos classified as those with  off / vir < 0.07 following Neto et al. (2007), tend to show larger differences in the shapes.This can be related with the fact that less-relaxed halos yield less precise shape measurements due to a wrong centre definition.However, we do not obtain a clear relation between the  CDM / SIDM − 1 and the halos masses or the environment characterisation,  5 .
In view of these results, we consider six halo subsets for the stacked analysis, classified considering their relaxation state, mass and environment.Although we do not obtain a clear relation between the shape differences and the environment, we decide to inspect samples selected based on  5 in order to assess the influence of the -2.45 -1.Relation between the halo properties and the shape differences observed in the CDM and SIDM simulations. CDM and  SIDM refers to the halos sphericities (/) measured in the CDM and SIDM simulations, respectively.In the left panel, we show the relation between  CDM / SIDM − 1 and the halo relaxation state characterised using  off / vir .The vertical line indicates the criteria used to select relaxed halos ( off / vir < 0.07), those halos selected as relaxed are at the left of the solid line in darker colours.The colour code indicates the halo virial masses.In the next panel, we show the normalised  CDM / SIDM − 1 distributions of the selected halo subsets, Lower-mass and Massive (see Table 1).The next panel shows the relation between the shape differences and halos masses.The colour code indicates the  5 parameter, lower values are related to halos in denser environments.Finally, in the last panel, we show the normalised  CDM / SIDM − 1 distributions for the Clustered and Isolated selected subsamples (see Table 1).neighbouring mass component on the analysis.In Table 1 we show the criteria used to define each sample and the number of halos included.Since the analysis is performed by projecting the particles in the three main directions of the simulations, these numbers are tripled in the stacking.In particular, we focus our analysis on relaxed halos since their shapes are better constrained and they are also less affected by systematic effects, such as the centre definition.
The samples selected according to the mass and environment are obtained according to the quartiles of the distributions of  vir and  5 to enhance the differences.

Stacked halo shapes
For a more direct comparison between the lensing estimates, we perform a stacked analysis of the particle distribution.In this approach, we only consider the bound particles since we intend to characterise the halo shapes.With this aim, we compute the projected shapes of the combined particle distribution for the different halo subsets defined in Table 1, previously rotating them according to the position angle, .We characterise the shapes of the combined halos by computing the shape tensor defined in Eq. 2, taking into account the positions of the rotated particles.We perform this procedure by considering the particles within varying radii to get an insight into the shape radial variation.Results are shown in Fig. 2. We highlight that the stacking procedure dilutes the scatter introduced by the shape variance in the halo sample and represents the shape radial variation of the mean particle distribution, instead of the mean shapes of the individual halos.Therefore, the shaded areas obtained according to the bootstrap resampling, do not represent the variance of the  distribution, for which we measure a standard deviation of ∼ 0.12 when considering all bound particles.The observed radial variation of the semi-axis ratios can be approximately modelled using a power law relation: (5) We show in Fig. 2 the fitted parameters,  0 and , for each halo subset.As can be noticed, halos in SIDM simulations tend to be more spherical across the entire range of radii considered.Halos in both simulations tend to be rounder towards the centre, especially those in the SIDM simulation.This is evidenced by the fitted slopes, , being in general ∼ 20% larger for the halos identified in the SIDM simulation.

HALO SHAPE ESTIMATES USING LENSING
The gravitational lensing effect is originated under the presence of a gravitational field which bundles the light rays of the luminous sources located behind it.In our case of interest, the gravitational field is generated mainly by the halo which hosts the galaxy cluster and the lensing effect magnifies and distorts the shapes of the background galaxies, i.e. those located behind the cluster.In the inner regions of the clusters, which correspond to the strong lensing regime, distortions are significant leading to the detection of arcs and multiple images.Conversely, in the weak lensing regime at the .Measured semi-axis ratios according to the stacked particle distribution within a radius , using the reduced tensor (see Eq. 2).The shadow regions corresponds to the errors set by a bootstrap with 50 iterations.Grey and pink correspond to the results for CDM and SIDM simulations, respectively.Fitted relations according to Eq. 5 are shown in red and black and the fitted parameters are shown in the legend.outskirts of the cluster, the produced effect is smaller and can be only assessed using an statistical approach.The magnitude of the distortion is related to the projected surface density distribution, Σ, and encoded in the shear parameter, , which is a complex quantity that can be obtained from the ellipticity components of the background galaxies.
A key limitation in weak-lensing studies is the relatively low signalto-noise ratio of the measured effect.Moreover, the observed distortions are combined with the galaxy intrinsic projected elongation,   .Hence, the observed galaxy ellipticity is a combination of both quantities, e ∼ e s + .The strategy that can be applied to isolate the lensing effect from the measured galaxy ellipticity, is to combine the measurements of many background galaxies located at roughly the same radial distance from the cluster centre.Then, by averaging the tangential ellipticity component of these galaxies, the shear estimator can be obtained as γ = ⟨e⟩.This approach assumes that the cluster is spherically symmetric so that the galaxies at the same radial distance are similarly affected by the lensing effect.Also, it assumes that the intrinsic orientations of background galaxies are random, which is valid when the galaxies combined are located in a wide redshift range.With this approach, the precision in the derived lensing estimates will depend on the number of galaxies combined.
The signal-to-noise ratio can be upgraded by combining the background galaxies of many galaxy clusters that share a similar observed property, such as the number of identified galaxy members or their X-ray luminosity.This strategy is commonly known as stacking techniques or galaxy-galaxy lensing.In this approach, the measured radial shear, is related to the mean surface density distribution of the combined clusters (e.g.Niemiec et al. 2017;McClintock et al. 2019;Pereira et al. 2020).In principle, this combination results in softening the surface distribution, by blurring the impact of substructure.However, if the main direction of the projected mass of each cluster is taken into account in the stacking procedure, this technique can be successfully applied in order to obtain the mean aligned projected elongation of the combined clusters ( In order to perform our analysis, we construct the lensing maps by stacking the particles related with the halos included in the different subsets defined in 3.1, previously rotating them according to the main axis direction of each halo, .To take into account that the lensing effect is sensitive to the whole particle distribution along the line-ofsight, all particles within a box of side-length 20ℎ −1 Mpc centred at each halo's position are taken into account.From the lensing maps we compute the radial profiles that can be modeled by considering three halo properties: mass, concentration and shape.In this section, we first summarise the formalism related to the lensing analysis and how the mean ellipticity of combined halos can be obtained using weaklensing stacking techniques.Then, we describe the computation of lensing estimators for combined halo samples and the derivation of mean halo elongation using this procedure.

Surface mass distribution and convergence
An elliptical surface mass density distribution can be modelled considering confocal elliptical isodensity contours, Σ(), where  is the elliptical radial coordinate,  2 =  2 ( cos 2 () + sin 2 ()/) (van Uitert et al. 2017).This distribution can be approximated using a multipole expansion in terms of the ellipticity defined as  := (1 − )/(1 + ) (Schneider & Weiss 1991): where we neglect the higher order terms in .In this approximation,  is the angle relative to the major semi-axis of the surface density distribution.Σ 0 and Σ 2 are the monopole and quadrupole components, respectively.Σ 0 is related to the axis-symmetrical mass distribution while the quadrupole component is defined in terms of the monopole as Σ 2 = − (Σ 0 ())/.The surface density distribution can be related with the lensing convergence parameter: where Σ crit is the critical density defined as: Here  OL ,  OS and  LS are the angular diameter distances from the observer to the lens, from the observer to the source and from the lens to the source, respectively.
In this work, we adopt a similar approach as presented in Gonzalez et al. (2022) and model the radial surface density distribution, Σ 0 , by taking into account two surface density components: the main halo host and the neighbouring mass contribution.This model assumes that the density distribution of the main halo (1−halo term) as well as the contribution of the neighbouring masses (2−halo term) are elongated by different amounts and that the main directions are correlated.Thus, although the direction taken into account to align the halos is related with the main halo component, the contribution of the neighbouring mass distribution to the quadrupole is also considered in this modelling.The expected misalignment between the two components will bias the estimated elongation of the neighbouring mass to lower values.
The main halo component is modelled assuming a spherically symmetric NFW profile (Navarro et al. 1997), which can be parametrised by the radius that encloses a mean density equal to 200 times the critical density of the Universe,  200 , and a dimensionless concentration parameter,  200 .The 3D density profile is given by: where   is the scale radius,   =  200 / 200 ,  crit is the critical density of the Universe and   is the characteristic overdensity: Therefore, to model the total surface density profile we consider an elliptical distribution for the main halo component with an elongation  1ℎ plus the term introduced by the neighbouring distribution also elongated and characterised by the aligned ellipticity component,  2ℎ , as: where Σ 1ℎ corresponds to the projected NFW profile (Eq.9) and Σ 2ℎ is the projected density of the neighbouring mass (Eq.11).
The quadrupoles are related to each monopole component with Finally,  is the position angle with respect to the major semi-axis of the halo mass distribution as in Eq. 12.The correspondent semi-axis ratios for each component are obtained as

Shear components
The tangential and cross shear components related with the observed background galaxy ellipticities, can be obtained from the deflection potential corresponding to the defined mass distribution and can also be decomposed into the monopole and quadrupole contributions: × (, ) =   ×,2 () sin(2).
These shear components are related with the surface density distribution through (van Uitert et al. 2017): , where  2 () is the quadrupole component of the lensing potential and is obtained as: By averaging the tangential shear in 13 in annular bins we obtain the usually defined density contrast ΔΣ, equivalent to the specified in Eq. 14, which is the only term observed in the case of an axis-symmetric mass distribution and it is only related with the monopole.
If we average the   and  × projections in annular bins, we can isolate the quadrupole components scaled according to the ellipticity: Here we define the distance independent quantities related with the quadrupole, Γ T and Γ × .We model the shear profiles described in Eq. 16, 17 and 18, considering the projected surface density as the sum of the surface components, Σ 0 = Σ 1ℎ + Σ 2ℎ , computed by using COLOSSUS1 astrophysics toolkit (Diemer 2018).

Profile computation and estimators
Shear profiles are computed from the surface density maps, Σ(, ), of the stacked halos, which are related to the convergence through Eq. 7. We first obtain the stacked Σ(, ) maps, by considering all the particles within a 20ℎ −1 Mpc size box, centred in each halo and projected along the three Cartesian axes of the simulation.In this way, we account for the line-of-sight contribution present in the lensing studies and we increase the number of halos considered in the analysis.Before all the particles are combined to compute the maps, we rotate their positions taking into account the main orientation, , based on the particle distribution of each halo (see 2.3).Therefore, the main direction of the stacked map is aligned with the horizontal axis, .Σ(, ) is then obtained by computing the particle density of all the considered halos within a grid of 8×8 (ℎ −1 Mpc) 2 area divided in 500×500 bins, leading to a pixel resolution of ∼ 32 × 32(ℎ −1 kpc)2 .This map is re-scaled by dividing each pixel density by the number of projected halos combined.This procedure is done for all the halo subsets defined in Table 1.Shear maps are then obtained from the density maps using lenspac 2 package, which applies the inverse of the KS-inversion (Kaiser & Squires 1993) to compute the  components.For illustration, in Fig. 3 we show the stacked maps derived for the total sample of halos analysed in SIDM and CDM simulations.It can be seen . Fitted contrast density profiles, ΔΣ, for the halo subsets analysed.In solid grey and dashed pink are the profiles obtained for the CDM and SIDM halos, respectively.Fitted relations are shown in solid and dashed lines, which are almost identical.In orange and light-green we show the corresponding 1-halo and 2-halo components, and the red lines represent the sum of these two components.The panels below each profile shows the differences between the observed contrast density and the fitted model (red solid lines).
included in annular bins: Here the sums runs over the  pix pixels which centres are included within a radial bin  ± .(Σ crit  t )  and (Σ crit  × )  are the −pixel value of tangential and cross shear maps, respectively, and   is the position angle of the −pixel with respect to the −axis.Profiles are computed considering 20 logarithmic annular bins from 100ℎ −1 kpc up to 5ℎ −1 Mpc.The respective errors for each profile,  ΔΣ ,  Γ T and  Γ × , are computed using a bootstrap resampling with 50 iterations over the total number of halos in each subset.In order to compute the errors we produce for each halo sample 50 resampled maps.

RESULTS
We perform the fitting procedure in order to obtain the parameters that characterise the surface density distribution, in the same way as it would be conducted using the observational data.First, we obtain the mass and the concentration parameters by fitting the contrast density profiles, ΔΣ, since this profile is sensitive only to the monopole component of the modelling and has a larger signal-to-noise ratio compared with the quadrupole profiles.After that, we constrain the projected halo shapes by fitting the quadrupole components.In this section, we begin by presenting the results related to the monopole component and then we present the halo shapes derived according to the quadrupole components, by considering two different models.

Fitted parameters from the monopole
We constrain fit the contrast density profiles by using the Markov chain Monte Carlo (MCMC) method, implemented through emcee python package (Foreman-Mackey et al. 2013), to optimise the loglikelihood function for the monopole profile: where ΔΣ is the profile computed from the tangential shear map according to Eq. 19,  ΔΣ its respective bootstrap error and ΔΣ is the model (Eq.16) computed considering Σ 0 as the sum of the main and second halo components, Σ 0 = Σ 1ℎ + Σ 2ℎ .To fit the data we use 15 steps and 250 steps with flat priors, 12.5 < log( 200 /(ℎ −1  ⊙ )) < 15.5, 4 <  200 < 8. Our best fit parameters are obtained after discarding the first 50 steps of each chain, according to the median of the marginalised posterior distributions and errors enclose the central 64% of the marginalised posterior.Fitted profiles are shown in Fig. 4 together with the reduced chisquares values.The adopted modelling is suitable for all the sub-sets Massive 1 3 .9 3 0 1 3 .9 4 5 1 3 .9 6 0 1 3 .9 7 5 1 3 .9 9 0 logM 200 but for the Isolated sample of halos, in which the amplitude of the 2-halo term overestimates the surface density distribution at larger radial scales.Fitted posterior density distributions for each sample are shown in Fig. 5. Mass and concentrations derived for each combined halo sample are in general in agreement within the errors for both simulated data-sets.However, profiles from the SIDM simulation tend to describe more massive (by ∼ 1%) and concentrated (by ∼ 3%) distributions.Previous studies have shown that density profiles in SIDM simulations tend to show higher concentrations which can be related with a mass displacement from the inside out, enhancing the density of SIDM halos around or beyond   , resulting in a steeper profile at such radii (Banerjee et al. 2020), thus in agreement with our findings.Although with low significance, differences in the masses might be related with the well-known interplay between mass and concentration (e.g.Bullock et al. 2001;Duffy et al. 2008;Bhattacharya et al. 2013;Okoli 2017;Ishiyama et al. 2021).
The combined relation between these parameters can be also noticed from the contour distributions.

Mean aligned elongation constrained from the quadrupoles
Once we have constrained  200 and  200 for each halo subset, we proceed to simultaneously fit the quadrupole components using a similar approach, minimising the sum of the likelihoods ln L (Γ T |, q1ℎ , q2ℎ ) + ln L (Γ × |, q1ℎ , q2ℎ ), defined as: where Γ T and Γ × are the profiles computed from the tangential and cross shear maps according to Equations 20 and 21, respectively,  Γ T and  Γ × their respective bootstrap error.Γ T and Γ × are the adopted models (Eq.17 and Eq.18) computed considering Σ 0 = Σ 1ℎ +Σ 2ℎ .To fit the data we use 15 chains for each parameter and 1000 steps in this case, considering flat priors: 0.6 < q1ℎ < 0.9 and 0.1 < q2ℎ < 0.5.The parameters are obtained after discarding the first 200 steps of each chain and the obtained posterior distributions for both fitted parameters are shown in Fig. 6.The fitted semi-axis ratios are tightly constrained and, as expected, they indicate rounder shapes for the halos in the SIDM simulation, with semi-axis ratios roughly ∼ 10% larger than their counterparts in the CDM simulation.An unexpected result is that the elongation of the neighbouring distribution of the halos tends to be larger for the SIDM simulation.Although the observed differences are not as significant as for the halo elongation, given that  2ℎ has larger errors, halo subsets from the SIDM simulation systematically show larger values, roughly around ∼ 10%.This result, together with the larger differences obtained for  1ℎ when compared to those based on the particle distribution, might reflect the deficiencies of the model in properly distinguishing between the shapes of both components.This is also accounted by the large differences between the modelling and the computed profiles quantified from reduced chi-squared values that can be observed in Fig. 7, es- pecially for the cross component and the halos included in the SIDM simulation.
A potential shortcoming in the adopted modelling approach might be linked to the parameterization used for fitting the main halo component.The impact of the adopted modelling for fitting the quadrupole was already tested in a previous work (Gonzalez et al. 2022), in which we considered three alternative approaches: (1) fitting an Einasto model (Einasto & Haud 1989;Retana-Montenegro et al. 2012) instead of the NFW, (2) fixing the NFW concentration in the analysis by using a concentration relation with mass and redshift and (3) restricting the radial range in the fitting procedure to avoid the impact of the neighbouring component.In that work, we showed that the fitted semi-axis ratios are in general agreement, regardless of the adopted modelling specially when discarding highly un-relaxed halos.This leads us to the conclusion that the constrained shapes are not expected to be significantly influenced by the particular model adopted.In this work, we extend this analysis by considering an NFW model with a core, since it is expected to better constrain the density distributions observed in SIDM halos.We adopt the model presented in (Neto et al. 2007) and fit three parameters to the contrast density profiles, mass, concentration and a parameter that accounts for the presence of a core (1/ in Eq. 2 in Neto et al. 2007).We obtain that the fitted parameters are highly correlated.This is expected since, as was already shown in the previous subsection, there is a correlation between the mass and the concentration that cannot be properly disentangled since the adopted radial range is not highly sensitive to the inner mass distribution.Despite this, we test the impact on the quadrupole fitting by using this modelling.We obtain equivalent results as for the NFW model without a core.The fitted semi-axis ratios are in agreement within 1%.The results obtained when varying the model further reinforces our findings, concluding that the chosen model for the main halo component has a limited impact on the estimated halo shapes based on the quadrupole fitting.However, it is worth noting that in contrast, the constrained masses and concentrations can be significantly influenced by the particular modelling.This result is expected given the low values obtained for the derived −squares when fitting ΔΣ, indicating that the adopted modelling accurately describes the observed density distributions with highly correlated parameters.

Radial variation modelling
According to the results presented in 3.2, the shapes of the halo particle distributions show a radial variation, which is more pronounced in the SIDM simulation, as indicated by a steeper fitted power law slope.Taking this into account and the results presented in the previous subsection, we aim to propose a model that accurately captures the shape radial variation.In order to do that, we generate a surface mass map that follows a projected NFW elliptical distribution, Σ(()), with a radial variation of the semi-axis ratio,  := (), set following Eq. 5. From this map, we can compute the expected quadrupole profiles using Equations 20 and 21.In Fig. 8, we show a comparison between the obtained profiles with and without considering a radial variation for the elongation of the density distribution.As it can be noticed, a radial variation of the semi-axis ratio can significantly affect the Figure 7. Fitted quadrupole profiles, tangential (upper panels) and cross (lower panels) components, for the halo subsets analysed.In solid grey and dashed pink are the profiles obtained for the CDM and SIDM halos, respectively.The shadow region express the obtained errors according to bootstrap resampling.Fitted relations are shown in solid and dashed lines.In orange and light-green we show the corresponding to the 1-halo and 2-halo components, respectively, and the red lines represent the sum of these two components. 1ℎ and  2ℎ are obtained according to the fitting procedure described in 5.2.The panels below each profile shows the differences between the observed quadropole profiles and the fitted model (red solid lines).
quadrupoles, especially at lower radial ranges, where the signal is suppressed due to the rounder distribution towards the centre.
In this case, we minimise the sum of the likelihoods defined in Equations 23 and 24, however, the models are modified.For the main halo component, the model for the quadrupole profiles is derived from surface density maps of an elliptical NFW model with a radial variation, Σ[()].In order to build the maps, we fix the previously obtained mass and concentrations for each sample, according to the fitted contrast density distributions, ΔΣ.We add to this model for the 1−halo term, the contribution from the neighbouring mass distribution considering the same model as in the previous subsection, i.e. using Eq. 17 and Eq.18 with Σ 0 = Σ 2ℎ .With this approach, we proceed to fit the quadrupoles considering three free parameters with uniform priors: −0.15 <  < 0.15, 0.4 <  0 < 0.8 and 0.3 <  2ℎ < 0.7.
The posterior distributions are presented in Fig. 10.In this case, it can be noticed that the semi-axis ratios for both simulations are consistent for all the halo subsets analysed.This is evidenced in both, the main halo component,  0 , and the neighbouring mass distribution,  2ℎ .On the other hand, significant differences are obtained for the power law slope that characterises the radial variation of the elongation.In particular, we obtain a detectable radial variation for the halos in the SIDM simulation, while the halos in the CDM simulation exhibit a negligible variation, compatible with  = 0.The proposed model seems to improve the fitting of the cross-quadrupole component especially for the halos in the SIDM simulation, resulting in lower chi-square values.However, according to the contours of the posterior distributions, there is also a noticeable interplay between the fitted parameters  0 and  which reflects that the signal is sensitive to a combination of these parameters.A deeper discussion of these results is presented in the next section.

DISCUSSION
Allowing strong elastic self-interactions of the DM particle produces more spherical shapes in the particle distribution of simulated halos.While this result is consistent and halos in SIDM simulations systematically exhibit rounder shapes, the actual quantitative differences are relatively small.This can be seen in the upper panel of Fig. 11 in which the differences of the semi-axis ratios of the particle distribution between CDM and SIDM simulations are lower than 5%.Due to the combination of the stacking power and the fact that the shear is sensitive to the whole inner distribution where the impact of considering a non-zero / is expected to be larger, weak-lensing estimates show significant differences between the fitted shape parameters from both simulations.
In the bottom panel of Fig. 11 we show the results when accounting for a radial variation of the elongation.In this case, the results based on the stacked particle distribution, show systematically lower values of  0 for the halos in the CDM simulation, but still in agreement within a 5%.On the other hand, as discussed in 3.2, the power law slopes differ by a ∼ 20%.Regarding the results predicted by the lensing analysis, there are no discernible variations in the fitted  0 and  2ℎ parameters.However, we found significant differences in the fitted power law slope, .Our results show a very subtle radial variation in the shapes of halos from the CDM simulation (−0.04 ≲  ≲ 0.0), whereas halos in the SIDM simulation exhibit a steeper relationship ( ≲ −0.06).
Further visualisation of these results is presented in Fig. 12, in which we show the predicted () relation for the main halo component, according to the fitted parameters,  and  0 , together with the fitted elongation for the neighbouring mass contribution,  2ℎ .In this Figure, it can be clearly seen how this modelling captures the expected variation of the elongation, in which at the inner radial ranges the halos in the SIDM simulation are rounder towards the centre compared to those in the CDM simulation.For all the samples in both simulations, the neighbouring mass distribution is rounder than the mass distribution of the main halo, except for the Clustered sample.This result indicates that the larger clustering amplitude expected for these halos, preferentially occurs in a direction aligned with the halo main orientation.We marked in darker colours the regions in which the main and the 2-halo component are expected to be more dominant, setting the limit at the mean virial radius of the stacked halos.At these regions, the semi-axis ratios of the main halo component from the SIDM and CDM simulations, converge roughly at the same value for most of the samples.However, Clustered and Lower-mass halo samples show higher differences within this region, being the halos in the SIDM simulation more elongated.
In order to interpret these results, we extend the analysis presented in 3.2 by considering all the dark matter particles used to calculate the lensing maps, as opposed to solely the bound particles, and computing the shapes up to projected radii of 5ℎ −1 Mpc.We show in Fig. 13 the ratio between the estimated shape parameters from the CDM and the SIDM simulations.When considering only the bound particles the differences between the simulated data are more significant at the inner regions and then remains constant given that we are far from the halo radial extension.On the other hand, the whole particle distribution shows a stepper variation and the semi-axis ratios converge at larger radii.Clustered and Lower-mass halo sub-sets show a similar behaviour as the modelled from the quadrupole profiles, in which the halos in the SIDM simulation show a more elongated distribution at the outskirts.This last result does not necessarily indicate that these halos are indeed more elongated at the outskirts.Instead it suggests that the neighbouring mass distribution is better aligned with the halo.We should keep in mind that the orientations considered to align the halos and perform the stacking are computed using only the bound dark-matter particles.This result is in agreement with the presented in Banerjee et al. (2020).In this work, they obtain from the phase-space of the particles in CDM and SIDM simulations, that the potential becomes more isotropic in the presence of velocity-independent self-interactions.
The observed shapes of the whole DM particle distribution clarify the obtained results presented in 5.2 and 5.3 and summarised in Fig. 11 and Fig. 12.When considering a radial variation, the fitted slopes of the power law are higher for the SIDM since the whole particle distribution shows a steeper radial variation than the predicted according only to the bound particles.On the other hand, results presented in 5.2, in which  2ℎ show significant differences between both simulations, is an expected result given that the model tends to compensate for the radial variation mainly present in the SIDM halos.In fact, fitted  2ℎ in the CDM simulation are in agreement for both models within a ∼ 3%, while the fitted values for the SIDM are 10% larger when considering a radial variation.

SUMMARY AND CONCLUSIONS
In this study, we provide a comprehensive analysis of the application of a stacked weak-lensing approach to determine the shapes of cluster halos, with the aim of gathering insights into the nature of dark matter particles.With this aim, we produce two sets of simulated data, with different cross-sections for the dark-matter particle, /, but with the same initial conditions.One of the simulated data considers a non-interactive DM particle, / = 0  2 / (CDM simulation), while the other allows for self-interactions with / = 1  2 / (SIDM simulation).In agreement with previous studies, derived shapes based on the particle distribution of the halos identified in the SIDM simulation are systematically rounder than those in the CDM simulation.However, assessing these differences through observations certainly constitutes a challenging task.
With the aim of obtaining the halo shapes to discern between the two simulated scenarios, we propose to measure the mean halo elongation by stacking the expected aligned component of the weaklensing signal.The approach considered in this work to recover the halo shapes, takes advantage of the stacking power,the sensitivity of the lens effect to the entire inner mass distribution within a projected radius and the projection of the signal in the main direction of the halo elongation.This results in radial profiles that efficiently capture the differences introduced due to the shape variation.In fact, our results show that the obtained fitted parameters show larger discrepancies between the simulated data, than when considering the particle distribution.Hence, the adopted methodology highlights the differences between the two simulated data-sets.
When applying the stacked lensing analysis, we found that the fitted shapes can vary by roughly a 10% for the halos identified in the SIDM simulation compared to those in the CDM simulation, with the latter tending to be more elongated.We have also tested a radial variation modelling that considers a relation between the mass elongation and the projected radial distance.The results show a higher slope of the power law for the shape radial variation for the halos in the SIDM simulation, describing rounder mass distributions towards the centre.The shape parameters roughly agree with the measured in CDM simulation at the mean virial radius of the stacked sample.
We also notice that all the lensing results are supported by the trends of the measured shapes based on the the whole particle distribution included in the analysis.For both models, with and without a radial dependent shape variation, we systematically obtain a more elongated distribution of the neighbouring mass component for the combined halos in the SIDM simulation.This result is more significant for the lower-mass halos and those in denser local environments.The results obtained agree with those presented in Banerjee et al. (2020) where a more coherent infall from larger scales is expected for the halos when considering self-interacting particles.
It is important to take into account that the results derived in this work are bound to the fact that we neglect the impact of the baryon physics in shaping the halos.Nevertheless, the fitted radial ranges include information well beyond the halo radial extension, where the impact of baryons is expected to be less significant.In fact, simulations that include baryon physics found that shape differences in the DM distributions persisted out to large radii and that the shape differences between CDM and SIDM were significantly larger than those between DM-only and hydrodynamical simulations (Robertson et al. 2019).Another effect that will be important to take into account is the estimate of the main-cluster orientation.In the anal-ysis presented we use the particle distribution in order to align the halos for the stacking.These orientations are in principle unknown in observational studies, which consider luminous proxies such as the galaxy member distribution in order to constrain them.Although the member distribution has been shown to offer a good prediction for the halo main orientation (Gonzalez et al. 2021b), hydrodynamical simulations will be necessary to test the impact of this misalignment in the analysis and to account for possible systematic effects.
Although other effects need to be assessed to improve the calibration of the presented methodology, its application shows a strong potential in capturing the differences introduced by dark-matter particle candidates.In view of the new upcoming large-scale surveys, such as the Legacy Survey of Space and Time (LSST, Ivezić et al. 2019) and Euclid (Laureĳs et al. 2011), that will provide high-quality measurements of lensing statistics, our findings show lensing techniques as a promising approach in order to test the nature of the dark-matter particle.11.Ratio between the fitted shapes parameters from each sample of halos tacked from the CDM and SIDM simulations.The shadow grey region corresponds to a 5% agreement.Upper panel: Red diamonds represent the mean semi-axis ratio obtained for the stacked halos, according to the bound particle distribution using the iterative method (2.3).Green and purple dots correspond to the fitted  1ℎ and  2ℎ from the quadrupole components (5.2) Bottom panel: Black and grey squares correspond to  and  0 fitted according to the radial variation of the stacked bound particle distribution (3.2).Light-blue, green and purple correspond to the fitted parameters, ,  0 and  2ℎ , according to the quadrupole profiles (5.3).
Córdoba (SeCyT-UNC, Argentina).Also, it was supported by Ministerio de Ciencia e Innovación, Agencia Estatal de Investigación and FEDER (PID2021-123012NA-C44).This work has been supported by the call for grants for Scientific and Technical Equipment 2021 of the State Program for Knowledge Generation and Scientific and Technological Strengthening of the R+D+i System, financed by MCIN/AEI/ 10.13039/501100011033 and the EU NextGeneration/PRTR (Hadoop Cluster for the comprehensive management of massive scientific data, reference EQC2021-007479-P) Figure2.Measured semi-axis ratios according to the stacked particle distribution within a radius , using the reduced tensor (see Eq. 2).The shadow regions corresponds to the errors set by a bootstrap with 50 iterations.Grey and pink correspond to the results for CDM and SIDM simulations, respectively.Fitted relations according to Eq. 5 are shown in red and black and the fitted parameters are shown in the legend.

Figure 3 .
Figure3.Maps of the stacked halos in the CDM (upper panels) and SIDM (middle panels) simulations, and their differences (lower panel).From left to right: The first map corresponds to the surface density map (Σ ( , ) = Σ  ), computed by stacking the particle distribution of the total sample of halos, taking into account all the particles included within a box of side-length 20ℎ −1 Mpc and after rotating the particles to get the halo major semi-axis aligned with the −axis.The particle positions are centred at each halo position and rotated according to the main axis direction of the halo.The second map is related to the tangential shear, Σ    ( , ), and leads to contrast density profile, ΔΣ, when averaged in annular bins.The third map corresponds to the tangential quadrupole component obtained after subtracting the contrast density distribution, ΔΣ ( , ), to the tangential shear map ( Σ   ,2 ( , )).ΔΣ ( , ) is computed with the fitted mass and concentration derived from adopted modelling.Finally, the last map corresponds to the cross component and is only related to the quadrupole signal, Σ   × ( , ).

Figure 5 .
Figure 5. Posterior density distributions for the fitted parameters from the contrast density profiles, ΔΣ, mass and concentration.In grey and pink are the results for the CDM and SIDM, respectively.The distributions are shown after discarding the first 50 steps of each chain.Vertical lines indicate the median values of the distributions adopted as log  200 and  200 .

Figure 6 .
Figure 6.Posterior density distributions for the fitted parameters ( 1ℎ and  2ℎ ) from the tangential and cross quadrupole component profiles, Γ  and Γ × .In grey and pink are the results for the CDM and SIDM, respectively.The distributions are shown after discarding the first 200 steps of each chain and the vertical lines indicate the median values.

Figure 8 .
Figure 8. Quadupole profiles for a NFW surface density distribution with  200 = 10 14  ⊙ and  200 = 5 elongated by a fix semi-axis ratio,  = 0.6 (purple solid line) and by a radial dependent semi-axis ratio  =  0   (green solid line).Bottom panels show the difference between both profiles (black dashed line).

Figure 9 .
Figure9.Same as in Fig.7. 0 ,  and  2ℎ are obtained according to the fitting methodology described in (5.3).

Figure 10 .
Figure 10.Posterior density distributions for the fitted parameters (,  0 and  2ℎ ) from the tangential and cross quadrupole component profiles, Γ  and Γ × .In grey and pink are the results for the CDM and SIDM, respectively.The distributions are shown after discarding the first 200 steps of each chain and the vertical lines indicate the median values.

Figure
Figure11.Ratio between the fitted shapes parameters from each sample of halos tacked from the CDM and SIDM simulations.The shadow grey region corresponds to a 5% agreement.Upper panel: Red diamonds represent the mean semi-axis ratio obtained for the stacked halos, according to the bound particle distribution using the iterative method (2.3).Green and purple dots correspond to the fitted  1ℎ and  2ℎ from the quadrupole components (5.2) Bottom panel: Black and grey squares correspond to  and  0 fitted according to the radial variation of the stacked bound particle distribution (3.2).Light-blue, green and purple correspond to the fitted parameters, ,  0 and  2ℎ , according to the quadrupole profiles (5.3).

Table 1 .
Criteria chosen to select the halo subsets for the stacked analysis and the number of halos included in each sub-sample.