The mass distribution in the outskirts of clusters of galaxies as a probe of the theory of gravity

We show that $\varsigma$, the radial location of the minimum in the differential radial mass profile $M^\prime(r)$ of a galaxy cluster, can probe the theory of gravity. We derived $M^\prime(r)$ of the dark matter halos of galaxy clusters from N-body cosmological simulations that implement two different theories of gravity: standard gravity in the $\Lambda$CDM model, and $f(R)$. We extracted 49169 dark matter halos in 11 redshift bins in the range $0\leq z\leq 1$ and in three different mass bins in the range $0.9<M_{200c}/10^{14}h^{-1}$M$_\odot<11$. We investigated the correlation of $\varsigma$ with the redshift and the mass accretion rate (MAR) of the halos. We show that $\varsigma$ decreases from $\sim 3R_{200c}$ to $\sim 2R_{200c}$ when $z$ increases from 0 to $1$ in the $\Lambda$CDM model. At $z\sim 0.1$, $\varsigma$ decreases from $2.8R_{200c}$ to $\sim 2.5R_{200c}$ when the MAR increases from $\sim 10^4h^{-1}$M$_\odot$~yr$^{-1}$ to $\sim 2\times 10^5h^{-1}$M$_\odot$~yr$^{-1}$. In the $f(R)$ model, $\varsigma$ is $\sim 15$% larger than in $\Lambda$CDM. The median test shows that for samples of $\gtrsim 400$ dark matter halos at $z\leq 0.8$, $\varsigma$ is able to distinguish between the two theories of gravity with a $p$-value $\lesssim 10^{-5}$. Upcoming advanced spectroscopic and photometric programs will allow a robust estimation of the mass profile of enormous samples of clusters up to large clustercentric distances. These samples will allow us to statistically exploit $\varsigma$ as probe of the theory of gravity, which complements other large-scale probes.


Introduction
According to the standard ΛCDM model, the gravitational dynamics of the Universe and its content is described by General Relativity in a flat space-time with a positive cosmological constant, Λ.The matter content of the Universe mainly consists of cold dark matter (CDM), and baryons are a minor component.In this model, cosmological structures form hierarchically through consecutive mergers of smaller halos into larger halos (e.g., Press & Schechter 1974;White & Rees 1978;Lacey & Cole 1993;Bower 1991;Sheth & Tormen 2002;Zhang et al. 2008;De Simone et al. 2011;Corasaniti & Achitouv 2011;Achitouv et al. 2014;Musso et al. 2018).
For a long time, the outskirts of clusters of galaxies have been regarded as a potentially powerful probe of the model of structure formation and evolution (e.g., Diaferio 2004;Reiprich et al. 2013;Diemer & Kravtsov 2014;Lau et al. 2015;Walker et al. 2019;Rost et al. 2021).At radii larger than the virial radius, this region is linked to the cosmic large-scale structure, which is sensitive to the model of gravity and to the nature of dark matter.Furthermore, unlike the central regions of clusters, comparisons of the dynamical properties of these regions with simulations are more robust because the relevant physical processes are dominated by gravity and not by baryonic physics (e.g., Diemand et al. 2004;Hellwing et al. 2016;Armitage et al. 2018;Velliscig et al. 2014;van Daalen et al. 2014;Shirasaki et al. 2018).
The most widely investigated feature of cluster outskirts is the splashback radius, R spl , which roughly corresponds to the first apocenter of the orbits of recently accreted material (Adhikari et al. 2014).This radius has an important physical meaning: It separates particles that are orbiting in the gravitational potential of the halo from matter falling onto it (Diemer & Kravtsov 2014;Adhikari et al. 2014;More et al. 2015;Diemer 2020Diemer , 2021)).For this reason, More et al. (2015) proposed it as the physical boundary of clusters.Diemer & Kravtsov (2014) and More et al. (2015) showed that this feature can be detected as a sudden drop in the logarithmic slope of the halo density profile.These works employed Nbody simulations, and they showed that R spl is usually located at ∼2R 200c1 and that it decreases with increasing mass accretion rate (MAR) and redshift.
The splashback radius is a genuine prediction of the hierachical model of structure formation.It is therefore expected to be sensitive to the model of gravity and the nature of dark matter and dark energy (Diemer & Kravtsov 2014;Adhikari et al. 2014Adhikari et al. , 2018;;Contigiani et al. 2019;Banerjee et al. 2020;Diemer 2020;Umetsu 2020).Adhikari et al. (2018) proposed the exploitation of R spl as a probe of cosmic expansion and the model of gravity for the first time.They used analytical models and N-body simulations to study the impact of dark energy and modifications in gravity on the value of R spl .They showed that this feature is sensitive both to the w parameter of the equation of state of the dark energy fluid and to two screened modified gravity models, f (R) and nDGP.They pointed out that the outskirts of galaxy clusters are the transition region between the screened and unscreened regimes of modified gravity models, and thus they are the ideal region in which to detect signatures of modified gravity.
In the data, several works detected a change in the slope of the mass density profile at a radius compatible with the expected R spl (e.g.,More et al. 2016;Baxter et al. 2017;Umetsu & Diemer 2017;Chang et al. 2018;Zürcher & More 2019;Shin et al. 2019;Murata et al. 2020;Bianconi et al. 2021;Gonzalez et al. 2021).However, although these results are encouraging, observational estimates are affected by both projection effects and by the fact that we observe galaxies rather than dark matter particles.Dark matter particles and dark matter subhalos, with which galaxies might be associated, do not necessarily share the same apocenter (Diemer et al. 2017).Recent investigations based on hydrodynamical simulations showed that R spl derived from galaxies or gas particles can substantially differ from R spl derived from dark matter (Baxter et al. 2021;Deason et al. 2021;O'Neil et al. 2021O'Neil et al. , 2022;;Dacunha et al. 2022).These discrepancies behave in a complex way that depends on the halo mass and selection criteria.Therefore, further observational and methodological efforts are required.
High-precision measurements require high-quality spectroscopic and wide-field surveys.These will be available in the near future by ongoing or planned experiments; in addition, an improved understanding of the systematics that affect the methods that are used to estimate the mass profile of clusters is necessary.Currently, weak gravitational lensing (Bartelmann 2010;Umetsu 2020) and the caustic technique (Diaferio 1999;Serra et al. 2011) are the only methods that estimate the cluster mass profile at large cluster-centric distances without relying on dynamical equilibrium.Weak-lensing is steadily improving with the acquisition of large combined spectroscopic and photometric datasets and methodological advancements (Umetsu 2013;Dell'Antonio et al. 2020;Umetsu 2020).Synergies with the caustic technique (Diaferio 1999;Serra et al. 2011) promise unbiased and accurate mass profiles that extend to large radii.
On the theoretical side, we need to thoroughly characterize the dynamics of clusters in their outer region.This plan requires investigating other characteristic radii of cluster outskirts in addition to R spl .Beyond R spl , clusters present an extended region in which matter is falling onto them.The radius at which the radial velocity profile has its minimum traces the region of clusters in which the dynamics is a genuine infall (De Boni et al. 2016;Vallés-Pérez et al. 2020;Pizzardo et al. 2021Pizzardo et al. , 2023b)).Characteristic radii also arise from the halo-matter bias (Fong & Han 2021;García et al. 2021).These radii depend on the cluster mass and redshift.Because these radii are tightly connected with the cluster MAR, they are probably extremely sensitive to the gravity model.Studies aiming at the quantitative determination of the dependence of these radii on the gravity model are vital for the exploitation of cluster outskirts as a cosmological probe.
Here, we focus on a radius that is different from both R spl and the radius of the minimum of the radial velocity profile.We consider the radial location, ς, of the minimum of the radial differential mass profile of galaxy clusters, M (r).We investigate whether ς can distinguish between standard gravity and f (R) gravity based on galaxy cluster halos in wide redshift and mass ranges.The radius ς complements R spl as a probe of the external region of clusters of galaxies: ς reduces the large errors in R spl that originate from the logarithmic slope of the density profile, and it can thus in principle be located in individual systems.We prove that although ς and R spl are not equivalent, they share similar trends with redshift, mass, and MAR.
The paper is organized as follows.In Sect. 2 we present the simulations (Sect.2.1) and the samples of simulated halos (Sect.2.2).In Sect. 3 we present our theoretical results: in Sect.3.1 we compare ς with R spl , and in Sect.3.2 we study whether ς can distinguish between two different gravity models.We discuss our results in Sect. 4 and conclude in Sect. 5.

DUSTGRAIN simulations
We extracted our samples of simulated halos from DUSTGRAIN (Giocoli et al. 2018), which is a suite of N-body cosmological simulations implementing the standard ΛCDM model and various extensions of General Relativity in the form of f (R).The simulations might include a non-negligible fraction of massive neutrinos.Appendix A reports a basic review of the f (R) model considered in these simulations.For our purposes, we used only the ΛCDM and the f R4 runs; the latter is an f (R) run with scalaron f R0 = −1 × 10 −4 and without massive neutrinos (Ω ν = 0).
The DUSTGRAIN simulations were normalized at the cosmic microwave background (CMB) epoch with cosmological parameters consistent with the Planck 2015 constraints (Ade et al. 2016).The cosmological matter density is Ω M ≡ Ω CDM + Ω b + Ω ν = 0.31345, where Ω CDM , Ω b , and Ω ν are the current CDM, baryonic, and neutrino density parameters, respectively.The DUSTGRAIN simulations had Ω b = 0. Furthermore, they had a cosmological constant Ω Λ = 0.68655, a Hubble constant H 0 = 100 h km s −1 Mpc −1 , with h = 0.6731, an initial perturbation amplitude A s = 2.199 × 10 −9 , and a power spectrum index n s = 0.9658.The normalization of the power spectrum of the perturbations at the CMB epoch implies that the value of the power spectrum normalization at z = 0, σ 8 , generally differs among the different simulations.In the two runs relevant for this work, ΛCDM and f R4, σ 8 = 0.847 and σ 8 = 0.967, respectively2 .The box size was 750 h −1 Mpc on a side in comoving coordinates.The simulation contained 768 3 collisionless dark matter particles with mass m DM = 8.1 × 10 10 h −1 M each.Because hydrodynamics is absent, the ΛCDM and the f R4 runs are two standard collisionless N-body simulations with Ω M = Ω CDM = 0.31345 and Ω b = Ω ν = 0.
The CDM halos were identified with a friends-of-friends (FoF) algorithm (Davis et al. 1985) with linking length λ = 0.16 d, where d is the mean Lagrangian interparticle separation.The SUBFIND algorithm (Springel et al. 2001) was run on the FoF groups to identify their gravitationally bound substructures.A synthetic cluster corresponds to the most massive detected substructure, whose most tightly bound particle is the cluster center.

Samples of simulated halos
For the theoretical investigations presented in this work, we considered the three-dimensional synthetic clusters of DUSTGRAIN (Sect.2.1) in 11 redshift bins in the range 0 ≤ z ≤ 1. Table 1 lists the main properties of the cluster samples at each redshift: the number of clusters, the median mass, and the median MAR with the 68th percentile ranges of these quantities.
The clusters were separated into low-, intermediate-, and high-mass samples.The low-mass samples contained all the simulated halos with mass M 200c ∈ (0.9, 1.1) × 10 14 h −1 M , and the intermediate-and high-mass samples contained all the halos with mass M 200c ∈ (4.0, 6.0) × 10 14 h −1 M and M 200c ∈ (9.0, 11.0) × 10 14 h −1 M , respectively.The high-mass samples are missing at z > 0.3 because the number of halos in this mass   range is ∼0−10 and these small samples would not provide statistically robust results.We had 26 independent samples of three-dimensional halos for each of the two models we considered: 11 samples of lowmass (L) halos, 11 samples of intermediate-mass (I) halos, and 4 samples of high-mass (H) halos.In total, the 26 ΛCDM ( f R4) halo samples are provided with the six phase-space coordinates of the dark matter particles in 20 887 (28 282) synthetic clusters; the higher number of halos for f R4 is the result of the higher σ 8 of this model.
As an effect of the cluster mass function, the low-mass samples have ∼2000 halos at each redshift (1736 and 2303, averaging over all redshifts, for ΛCDM and f R4, respectively), the intermediate-mass samples have ∼200 halos (154 and 251, averaging over all redshifts, for ΛCDM and f R4, respectively), and the high-mass samples have ∼35 halos (25 and 46, averaging over all redshifts, for ΛCDM and f R4, respectively).The redshift dependence of the mass function has clear consequences for the size of the samples in the same mass bin but different redshift: for the low-mass samples, the number of halos at z = 0.00 is ∼3−4 times higher than at z = 1.00; for the intermediatemass samples, this ratio is ∼28−38.

Results
We now describe our results.In Sect.3.1, we show the relation between the splashback radius R spl and the radial location ς of the minimum of the derivative of the cumulative radial mass profile.In Sect.3.2 we investigate the dependence of ς on redshift, MAR, and the theory of gravity.

Splashback radius and the radial location ς of the minimum of the differential mass profile
In this section, we compute the radial logarithmic derivative of the density profiles of the synthetic clusters listed in Table 1 and identify R spl .We then compare R spl with the radial location ς of the minimum of the derivative of the cumulative radial mass profiles.
For each redshift and mass bin, we computed the average radial shell mass density profile, ρ shell (r), of the halos as follows.For each halo, we computed the mass within 200 adjacent spherical shells progressively from the center.The shells had outer cluster-centric radii spanning the range (0.1, 10)R 200c , and they were logarithmically spaced.We then computed the A80, page 3 of 13  1.The left, middle, and right panels show the average profiles for the low-, intermediate-, and high-mass samples, respectively, at the lowest and highest redshifts of each sample, as listed in the insets.The solid (dash-dotted) curves are the profiles in the ΛCDM ( f R4) model.The vertical solid (dash-dotted) thick black lines show the location of R spl at z = 0 in the ΛCDM ( f R4) model.The pair of vertical blue (green) solid and dash-dotted lines shows R spl at z = 1.0 (z = 0.3) for the ΛCDM and f R4 models, respectively.Bottom panels: the thin solid (dash-dotted) black curves show the radial derivative of the average mass profile in ΛCDM (fR4) at z = 0; the vertical thin solid (dash-dotted) black lines show ς avg of the low-, intermediate-, and high-mass samples at z = 0 in the ΛCDM ( f R4) models (Table 2, first row, values in brackets).When extended to the upper panels, these thin solid and dash-dotted black lines cross the average density profiles at d log ρ shell /d log r −2, as expected from Eq. (1).
average shell mass profile of each halo sample: The value of the mass at each radial bin r i of this profile is the mean of the individual shell masses of all the halos in the sample at r i .Each value of these average profiles was determined by a number of data equal to the size of the corresponding sample (columns "size" in Table 1).From these average profiles ρ shell (r), we also computed the logarithmic derivative, d log ρ shell /d log r.This procedure automatically includes the mean background matter density and H(z) and G are the Hubble parameter and the gravitational constant, respectively.This procedure is the same as the approach adopted for real clusters.
In Fig. 1, the curves in the upper left, middle, and right panels show examples of the logarithmic derivative of the average shell density profile in the low-, intermediate-, and high-mass bins, respectively.Each profile is smoothed with a Savitzky-Golay filter (Savitzky & Golay 1964) with a window length of 25 bins.Colors and line styles show the dependence of the profiles on redshift and theory of gravity, respectively.
All the simulated average profiles of d log ρ shell /d log r show an absolute minimum.Adhikari et al. (2014) associated this feature with R spl , which roughly corresponds to the first apocenter of the orbits of recently accreted material.We identify R spl of each halo sample as the point of absolute minimum of the logarithmic derivative of the associated smoothed average density profile ρ shell (r).In the upper panels of Fig. 1, each vertical line shows the point of minimum, R spl , of the logarithmic derivatives of ρ shell (r) for the corresponding redshift (color) and theory of gravity (style).
As shown in Table 2, R spl decreases with increasing redshift at fixed mass, and with increasing mass at fixed redshift.At z ≤ 0.7, at fixed redshift and mass, ΛCDM systematically predicts a smaller R spl than f R4: From z = 0 to z = 0.7, the relative differences between R spl in f R4 and ΛCDM is in the range ∼2−15% in the low-mass bin, and ∼2−35% in the intermediatemass bin.In the high-mass bin, the same difference increases from ∼5% at z = 0 to ∼55% at z = 0.3.At z ≥ 0.8, the lowmass bin confirms this behavior: R spl in f R4 is larger than in ΛCDM by ∼2−12%.In contrast, in the intermediate-mass bin, at z = 0.8−1.0,R spl in ΛCDM can be ∼22% larger than in f R4.
A major issue in the use of the d log ρ shell /d log r profile is that when it is computed for real clusters from the estimated mass profile, the uncertainty in the ith bin is proportional to r i /(r i+1 − r i ), where r i and r i+1 are the radii used to compute the numerical derivative of the shell mass density profile, ρ shell (r).The uncertainty is thus increasingly large at increasing radius.Furthermore, the dependence of R spl on ρ shell (r) hinders the robust estimation of R spl due to the scatter caused by the high sensitivity of the logarithmic derivative to rapid changes of ρ shell (r) in two adjacent shells.We can overcome these problems by considering, instead of ρ shell (r) and R spl , the radial differential mass profile dM/dr = M (r) and the radial location, ς, of the minimum of this profile.Notes.The values reported in brackets refer to the point of the minimum of the radial derivative of the average mass profile, ς avg .
R spl and ς are not equivalent: The condition for R spl is d[d log ρ shell /d log r]/dr = 0, whereas the condition for ς is The values in brackets in Table 2 show the points of minimum ς avg computed from the radial derivative of the averaged cumulative mass profiles of the halos in each bin.The ς avg appear larger than R spl : By averaging over all the masses and redshifts, ς avg exceeds R spl by 43% and 40% in the ΛCDM and fR4 model, respectively.This result agrees with the recent study by García et al. (2021), who reported a ratio of the minimum of r 2 ξ hm , where ξ hm is the halo-matter correlation function, and R spl of ∼1.3.In Fig. 1, the thin solid (dash-dotted) black curves in the bottom panels show the M average profiles in ΛCDM (fR4) at z = 0; the thin black vertical solid (dash-dotted) lines show the points of minimum of these curves, namely ς avg in ΛCDM (fR4) at z = 0. Despite the obvious difference between ς avg and R spl , these two quantities share remarkable similarities in their dependence on redshift, mass, and gravity model.Table 2 shows that ς avg , similarly to R spl , decreases with increasing redshift at fixed mass, and with mass at fixed redshift.Moreover, analogously to R spl , in all the low-mass samples the ς avg in f R4 are larger than in ΛCDM by ∼10%, except in the low-mass sample at z = 1.0,where the ς avg in f R4 is ∼2% smaller than in ΛCDM.Similarly, in all the high-mass samples, the ς avg in f R4 are larger than in ΛCDM by 5−55%.In the intermediate-mass bin, the similarity with R spl also holds globally: At z ≤ 0.4, the ς avg in f R4 are larger than in ΛCDM by 5−35%.At higher redshifts, this hierarchy can be inverted, with ς avg in ΛCDM larger up to 29% than in f R4.
We assessed how well ς can distinguish between the ΛCDM and f R4 models through its dependence on redshift and mass.To do this, in addition to ς avg extracted from the average mass profiles and discussed above, we used a different estimate of ς extracted from the individual mass profiles.In the simulations, the M (r) profiles are less sensitive to sudden mass differences caused by random variations in the number of dark matter particles in two adjacent shells.We are thus able to exploit the ςs of the individual halos in our analysis.More precisely, for each redshift and mass bin, we computed the cumulative radial mass profiles of all the halos.The profiles were obtained by progressively adding from the center the mass of 200 adjacent spherical shells with logarithmically spaced outer cluster-centric radii spanning the range (0.1, 10)R 200c .For each halo, we computed the corresponding radial derivative of the mass profile.
In Fig. 2, the upper, middle, and bottom panels show the median radial differential mass profiles M (r) of the low-, intermediate-, and high-mass samples of Table 1, respectively, with their 68th percentile ranges.The solid (dash-dotted) curves are for the ΛCDM ( f R4) samples: for each sample, the value of the corresponding curve at each radial bin r i is the median of the set of the individual M (r i ) of all the halos in the sample; each point of these curves is determined by a number of data equal to the size of the corresponding sample (columns "size" in Table 1).
We identified the final ς of each mass and redshift sample as the median of the set of all the ςs of the individual halos in the sample.The ς of each individual halo is the minimum of the corresponding profile M (r) beyond 0.5R 200c .This choice avoids possible inconsistencies caused by high fluctuations in the inner region of the halo.In Fig. 2, the solid (dash-dotted) vertical lines show the value of ς in ΛCDM ( f R4) at the lowest and highest redshift of the three mass bins.In general, the position of ς does not coincide with the minimum of the corresponding median profile in Fig. 2 because we did not select ς at the minimum of the median profile of the individual halos, but as the median of the ςs of the individual halos.Table 3 lists the median ς of the individual differential mass profiles for the samples listed in Table 1.We also provide an estimate of the error of the median with the difference between the median and the corresponding 16th percentile, and the difference between the median and the 84th percentile value divided by the square root of the sample size.We also display the relative differences of ς in ΛCDM and f R4. Figure 3 shows ς as a function of redshift and mass in ΛCDM and f R4.
These new estimates of ς confirm the previous trends shown by ς avg , but the results are less prone to large scatter: ς decreases with increasing redshift and mass; furthermore, at fixed redshift and mass, ΛCDM generally predicts a smaller ς than f R4.As shown in Table 3, in the low-and high-mass bins, ς is systematically larger in f R4 than in ΛCDM: In the low-mass bin, the relative difference decreases with redshift from ∼12% at z = 0 to 0% at z = 1; in the high-mass bin, the ςs in f R4 are increasingly larger than in ΛCDM as redshift increases, from ∼19% at z = 0 to ∼25% at z = 0.3.In the intermediate-mass bin, ς is systematically larger in f R4 at redshift z ∼ 0 − 0.5, where the relative difference decreases from ∼12% to ∼5%.In this mass A80, page 5 of 13  z Low-mass samples Intermediate-mass samples High-mass samples  2), ΛCDM can present larger ςs than f R4 at z > 0.5.

Dependence of ς on redshift, MAR, and the theory of gravity
In the previous section, we showed that ς, analogously to R spl , is sensitive to the redshift and the mass of the cluster, and that it can be up to 25% larger in f R4 than in ΛCDM (Table 3).Mass and MAR are tightly correlated (Diemer & Kravtsov 2014;Pizzardo et al. 2021).Our findings therefore agree with previous results showing that R spl is correlated with redshift and MAR (More et al. 2015).In this section, we model the behavior of ς with z and MAR as independent variables, and we investigate whether ς can distinguish between ΛCDM and f R4.
Our definition of the MAR can be easily extended to real clusters because it does not rely on the merger history of the halos, as is usually adopted in investigations based on N-body simulations.We considered the procedure by Pizzardo et al. (2021Pizzardo et al. ( , 2022)), who measured the MAR of ∼500 real clusters of galaxies.The recipe is based on the assumption that the infall of new material onto the clusters is described by the spherical accretion model (De Boni et al. 2016).Namely, in the infall time t inf , a cluster accretes all the material within a spherical shell, centered on the cluster center, with inner radius R i and thickness ∆ s = δ s R i .The quantity δ s is derived from the equation of motion of the infalling mass with initial velocity v inf and constant accel- where M(< R i ) is the mass of the cluster within R i , and G is the gravitational constant.Here, according to De Boni et al.
(2016) and Pizzardo et al. (2021Pizzardo et al. ( , 2022)), we fixed t inf = 1 Gyr and R i = 2R 200c , whereas for the infall velocity, v inf , we relied on the radial velocity profiles of the simulated halos of DUSTGRAIN in the two theories of gravity.For a sample of simulated clusters, v inf is the minimum velocity of the median radial velocity profile in the range [2, 2.5]R 200c of the sample, whereas for a sample of real clusters, v inf results from suitable interpolations of the simulated samples (see Pizzardo et al. 2021Pizzardo et al. , 2022, for further details).When the mass of the infalling shell, M shell , is known, the MAR is estimated with the equation We computed the MARs of all the halos in the samples of Table 1.This table lists the median MAR of each sample with their corresponding 68th percentile range.In Fig. 4, the open circles show the median MARs of the 26 ΛCDM samples as a function of the median mass of the sample, color-coded by redshift.By considering the merging history of simulated CDM halos, Diemer et al. (2017) proposed an accretion law that we computed at our 11 simulated redshifts and show with the curves in Fig. 4. The agreement between this law and our estimates shows that the MARs of galaxy clusters estimated with our recipe, which can be applied to real clusters, are consistent with the MARs derived from the merger trees, which remain inaccessible to observations.The estimated MARs of the CIRS and HeCS clusters (Pizzardo et al. 2021, squares in Fig. 4) and of the HectoMAP clusters (Pizzardo et al. 2022, diamonds in Fig. 4) agree with the theoretical expectations of ΛCDM within the uncertainties.Both M 200c and MAR are estimated from the same caustic mass profile, hence their errors are not independent.However, the error of the caustic mass profile at a given radius is determined by the number and the distribution of the galaxies within the spherical shell centered at this radius, and M 200c and MAR are estimated at well-separated radii.Therefore, we can assume that the errors on these quantities are independent, and we obtain a roughly correct estimate of the expected uncertainty ranges (error bars in Fig. 4).
For each of the two sets of simulations, we fit the 26 points in the four-dimensional space (ς, z, MAR, M 200c ) with the relation where the scale factor, and a, b, c, and d are four free parameters.Figure 4 proves that the difference between the M vir s at two different epochs is comparable with the difference between the masses of the infalling shell that we adopted in our definition of MAR.With an obvious change in the variables, we A80, page 7 of 13 We applied small redshift offsets to the f R4 points to facilitate comparison.The black (red) lines are the fits of Eq. ( 4) to the ΛCDM ( f R4) points: specifically, the curves are the result of a spline interpolation of the values of ς (11 for the low-and intermediate-mass bin, 4 for the high-mass bin) resulting from Eq. ( 4) using the redshift, mass, and mass accretion rate of each simulated bin (Table 1).In each panel, the lower plot shows the differences of the median ςs in ΛCDM and f R4 relative to the ΛCDM value.
can therefore express the dynamical accretion rate as Γ dyn ≈ [K(t, Ω M , H 0 )/M 200c ] (dM/dt), with dM/dt ≡ MAR 3 .Table 1 reports Γ dyn for each cluster sample.Equation (4) was proposed by More et al. (2015) to model R spl rather than ς.The similar behavior of these two quantities with redshift, MAR, and mass leads us to heuristically use the same functional form for ς, however.We also kept Ω M (z) in the form of ΛCDM in the fit to the f R4 data in order to direct any difference to the four fitting parameters a, b, c, and d.We fit the median points with a nonlinear least-squares Marquardt-Levenberg algorithm by weighting 3 K transforms the time coordinates from the scale factor to cosmic time: d ln â = K −1 (t, Ω M , H 0 )dt.  the points with the inverse of their variance (Table 3).The values of the coefficients are reported in Table 4. Table 4 shows that the ΛCDM fit is robust.The parameter b, which controls the dependence of ς on redshift, is constrained with a relative uncertainty of ∼10%; the parameters c and d, which control the dependence of ς on the logarithmic MAR, are constrained with uncertainties of ∼45−47%.This result shows that in ΛCDM the dependence of ς on redshift and MAR are statistically significant and are similar to those expected for R spl .The f R4 fit is less strongly constrained: Except for the parameter b with a ∼10% uncertainty, the remaining three parameters are affected by uncertainties of ∼100% or larger.However, the model of Eq. ( 4) was only proposed to fit the relation in the ΛCDM universe.The large uncertainties on the parameters for the f R4 fit can thus simply suggest that Eq. ( 4) is less adequate for this gravity model.
The absolute value of the relative differences between the fit and the actual values is ∼4.2% and ∼3.5% for the ΛCDM and the f R4 samples on average, respectively.Figure 3 shows examples of these fits in the ς−z plane.Figure 5 shows projections of the fits on the ς-MAR plane at three different redshifts.The correlation of ς with redshift and MAR is confirmed by the Kendall test, which returns a tight correlation with a p-value 10 −5 for the relation ς−z and for the relation ς-MAR.
As expected, the fits show that the difference between the ΛCDM and f R4 models is larger at low redshifts, independently A80, page 8 of 13 of the MAR.The difference steadily decreases with increasing redshift: consequently, because at low redshift ς is higher in f R4 than in ΛCDM, the fits in the (ς, z) plane at fixed MAR (Fig. 3) show steeper profiles in f R4 than in ΛCDM.We now quantitatively assess to which extent ς can distinguish between ΛCDM and f R4.We applied the median test to each of our 26 ΛCDMf R4 pairs of the distributions of the individual ς at equal redshift and mass bin.This test counts the number of elements of the two distributions that are either larger or smaller than the grand median of the two samples.Based on a Pearson χ-squared test, the test assigns a p-value to the nullhypothesis that the samples are drawn from the same distribution.Generally, the f R4 sample is more densely populated than the corresponding sample in ΛCDM (see Table 1).To perform the test on samples with the same number of halos, we therefore randomly undersampled each f R4 sample to the number of halos in the corresponding ΛCDM sample.The results of the median test are shown in Fig. 6: the blue points in the upper panel are for the low-mass samples, whereas the orange and green points in the lower panel are for the intermediate-and high-mass samples, respectively.
We took p = 10 −3 as the upper limit to conclude that the two theories of gravity can be distinguished: Values below this limit are evidence at a significance level higher than 3σ that the two samples are drawn from different populations.Figure 6 shows that the low-mass samples can generally distinguish between the two models in the redshift range z ∼ 0−0.8.The samples at the lowest redshifts show extremely strong evidence that gravity is modified (p ≈ 10 −30 ).In contrast, only 3 of 11 intermediatemass samples show some evidence of modified gravity, and none of the high-mass samples is able to do so.
The progressively weaker discerning power of the median test in the low-mass samples as redshift increases is due to the intrinsic degeneracy of f R4 and ΛCDM at high redshift.We illustrate in Appendix A that the main differences between the two models are more evident at late times.This conclusion is also supported by the steady decrease in the relative difference between the ςs of the two models as redshift increases (see Table 3).In contrast, the largest part of the intermediate-and high-mass samples shows p > 10 −3 even at low redshift, where the relative difference between the values of ς is similar to or larger than for the low-mass halos (Table 3).The reason for this behavior is the small number of halos in these mass samples and not the intrinsic similarity of the models.
To test the effect of poor sampling, we considered for the low-and intermediate-mass bins the corresponding sample with the lowest number of halos, regardless of redshift and model: Table 1 shows that the smallest samples are the ΛCDM samples at z = 1.0, with 646 and 11 halos for the low-and intermediatemass bins, respectively.We therefore performed a median test in which we randomly undersampled all the samples of the lowand intermediate-mass bins to 646 and 11 halos, respectively.
The results are shown by the gray crosses in Fig. 6.In the low-mass bin, the undersampling increases the p-value of the median test for the samples at z ∼ 0−0.8; furthermore, the increment of the p-value decreases with redshift because the higher the redshift, the fewer halos are removed.Nevertheless, at z < 0.9, although they are substantially reduced, the p-values remain below 10 −3 .A80, page 9 of 13 The important result of this test is the minimum size of the sample required to distinguish between the two theories of gravity: As long as the sample size is 400, the distinguishing power remains substantially unaffected.Therefore, in the low-mass bin, the loss of significance of the median test to distinguish the theory of gravity at high redshifts, z ∼ 0.9−1.0, is caused by the intrinsic degeneracy between f R4 and ΛCDM at these early times and not by the sample size.Conversely, the failure of the distinguishing ability of ς in the intermediate-and high-mass bins is caused by the small halo samples.As long as the sample size is ∼230−400, as in the first four redshift bins of the full intermediate-mass bin, the distinguishing power is still present.With smaller sample sizes, however, the p-values are aligned to the results of the intermediate-mass full samples at z 0.4, and of the high-mass full sample, where we have ∼10−150 halos.

Discussion
We have investigated the radial location, ς, of the minimum of the radial differential mass profile, M (r), of simulated galaxy clusters halos.Although not equivalent, we show that ς shares similarities with the splashback radius, R spl .We show that ς can distinguish among different gravity models.We investigated General Relativity, with the standard ΛCDM model, and f (R) gravity.Our work strengthens the hypothesis that galaxy clusters might be used as cosmological probes.In particular, we showed that the outskirts of clusters promise to effectively complement more traditional cluster-based cosmological tests, such as cluster counts and the gas fraction (see Cataneo & Rapetti 2020, for a review of the tests of gravity with galaxy clusters).
Like R spl , ς can be located directly from the mass profiles of galaxy clusters.Most works that searched for R spl relied on the projected galaxy number density profile or on the projected luminosity profile.These quantities are tightly related and are good proxies of the three-dimensional matter density profile in general (Carlberg et al. 1997;Rines et al. 2003Rines et al. , 2004;;Rines & Diaferio 2006;Rines et al. 2013).However, many studies have pointed out that the galaxy number density and luminosity profiles are not directly proportional to each other (e.g., Proctor et al. 2015).In addition, their correlation is sensitive to the cluster mass and to the selection function (Mulroy et al. 2014).Moreover, if the galaxy density is constrained by photometrically detected clusters, projection effects and member mismatch become increasingly relevant in the outer regions (e.g.,More et al. 2016;Baxter et al. 2017;Shin et al. 2019).The use of the mass profile to detect R spl or ς is clearly preferred because the definitions of these quantities are based on the mass profile that is measured in N-body simulations.
Furthermore, the detection of ς from mass profiles is prone to smaller fluctuations than the detection of R spl because the radial derivative M (r) is much less noisy than the logarithmic slope of the matter density.For this reason, unlike R spl , ς can be located in individual clusters if the mass profile is accurately estimated up to ∼3−4R 200 .
Using ς as a probe of the theory of gravity requires an unbiased estimate of hundreds of cluster mass profiles up to large cluster-centric radii with high-precision and accuracy.Weak gravitational lensing (Bartelmann 2010;Hoekstra et al. 2013;Umetsu 2020) and the caustic technique (Diaferio & Geller 1997;Diaferio 1999;Serra et al. 2011) are the only methods that can be reliably used beyond the inner region of clusters, where dynamical equilibrium does not hold.
The current cluster samples do not allow us to robustly locate ς and to statistically use it as a probe of the theory of gravity.However, future facilities such as the VRO (Ivezić et al. 2019) and the Euclid mission (Laureijs et al. 2011;Sartoris et al. 2016) will provide weak-lensing mass measurements of thousands of clusters.Powerful spectrographs such as the multiobject William Herschel Telescope Enhanced Area Velocity Explorer spectrograph on the WHT (WEAVE, Balcells et al. 2010;Dalton et al. 2012Dalton et al. , 2014;;Dalton 2016) will explore the infall regions of galaxy clusters and provide dense spectroscopic surveys out to radii 5R 200c of more than 100 clusters at redshift 0.5.Planned observations with the Prime Focus Spectrograph on Subaru (PFS, Tamura et al. 2016) and the Maunakea Spectroscopic Explorer on the CFHT (MSE, Marshall et al. 2019) will provide spectroscopic redshifts of hundreds to thousands of galaxy cluster members for thousands of individual clusters at z 0.6.These surveys will allow a precise and unbiased determination of ς from weak lensing and the caustic technique.
These future observations will enable synergies between photometric and spectroscopic techniques.These synergies will be crucial for promoting ς as a probe of gravity.At clustercentric distances of about ς, a misidentification of background galaxies and large-scale structure can severely bias weak-lensing mass profiles obtained from photometric redshifts (Hoekstra 2003;Hoekstra et al. 2011;Becker & Kravtsov 2011;Sommer et al. 2022).The caustic technique provides unbiased mass profiles on average, but the scatter is relevant (Serra et al. 2011;Pizzardo et al. 2023a).The joint use of these two techniques could take advantage of the respective strengths, and it would be crucial to obtain exquisitely robust mass estimates up to large cluster-centric distances (e.g., Dell'Antonio et al. 2020).
As a cautionary note about our analysis presented in this work, we recall that the ΛCDM and f R4 runs of the DUSTGRAIN simulations that we used here were normalized at the CMB epoch.Therefore, their power-spectrum normalization at z = 0, σ 8 was 0.847 and 0.967 for the two models (Sect.2.1).To show that our results are not affected by these different values of σ 8 , we compared the median MARs at z = 0 of our low-and high-mass ΛCDM samples with the median MARs estimated by Pizzardo et al. (2021Pizzardo et al. ( , 2022) ) in the ΛCDM L-CoDECS N-body simulation of Baldi (2012), which has σ 8 = 0.809.Despite the lower σ 8 , we find that the MARs in the L-CoDECS simulation are ∼14% and ∼3% larger than the MARs of our DUSTGRAIN samples for the low-and the high-mass bin, respectively.Moreover, ς is ∼4.5% smaller in the DUSTGRAIN clusters than in the L-CoDECS simulations for both the low-and high-mass bin.In addition, in DUSTGRAIN, σ 8,ΛCDM < σ 8, f R4 , which means that if ς were inversely proportional to σ 8 , we should expect ςs larger in ΛCDM than in f R4 at fixed redshift and mass.This expectation is the opposite of the results from the simulations (Table 3).Hence, the behavior of ς does not appear to be correlated with σ 8 , and therefore, the differences between the models that we find here are mainly driven by the different theory of gravity.
One might expect that the outer regions of clusters smoothly approach the uniform cosmic background and that this background contributes to the mass profile of the clusters.However, the distribution of matter around both simulated (e.g., Gao et al. 2004;van den Bosch et al. 2005;Cornwell et al. 2023Cornwell et al. , 2024) ) and real clusters (e.g., Beers et al. 1983;Dressler & Shectman 1988;Martínez et al. 2016;Sarron et al. 2019) is inhomogeneous, with filaments, groups of galaxies, and low-density regions.Therefore, we did not subtract the contribution of a uniform background from the mass profile the clusters in our analysis (see Sect. 3.1).Nevertheless, we verified that the effect of subtracting a uniform background with density Ω M ρ c from the mass profiles has a negligible effect on our results: The radii ς at A80, page 10 of 13 which the minimum of the radial derivative of the cluster mass profile occurs are ∼10% larger on average than the ςs from the profiles without subtraction.More importantly, the correlations between ς, redshift, and MAR are preserved at similar significance levels.
We showed that extended ranges of mass and redshift are critical to avoid possible degeneracies between the models.Thus, our analysis pointed out the importance of further efforts toward the realization of simulations of large cosmological volumes in modified gravity.Larger samples that include the most massive clusters can help us to establish the correlation between ς and MAR with more significance because these samples will allow us to stack clusters based on mass and to probe a more extended range of MARs at fixed redshift.Furthermore, the most massive clusters, with a mass ∼10 15 h −1 M , are expected to provide the most robust evidence of possible deviations from the standard scenario because these clusters probe the upper tail of the mass function, which is extremely sensitive to the cosmological model (Bardeen et al. 1986;White 2002;Courtin et al. 2011;Cui et al. 2012;Giocoli et al. 2018), and because they are the easiest objects to detect observationally.In Sect. 3 we showed that in the high-mass samples of simulated clusters, we observe the largest discrepancy between the values of ς in the two models of gravity (Table 3).Nevertheless, the larger sizes of the low-and intermediate-mass samples, compared to the high-mass sample, enable the best distinction between the two models of gravity (Sect.3.2).On the one hand, large simulated volumes can therefore provide samples with an adequate size for a solid test bench based on which the performance of ς as a gravitational probe can be investigated.On the other hand, low-and intermediate-mass clusters, with ∼1−5 × 10 14 h −1 M , can already provide observed samples that are sufficiently large to yield tight constraints on the theory of gravity.

Conclusion
We performed the first assessment of the ability of the radial location, ς, of the minimum of the differential radial mass profiles, M (r), of clusters of galaxies to distinguish between two different models of gravity, General Relativity, with the ΛCDM model, and f (R).We investigated the behavior of ς in a wide range of masses, M 200c ∼ (0.9−11) × 10 14 h −1 M , and redshift, 0.0 ≤ z ≤ 1.0.We explored the dependence of ς on redshift and MAR in the two models.We also showed that ς has a similar redshift, mass, and MAR dependence as the splashback radius.
We used the DUSTGRAIN N-body simulations (Giocoli et al. 2018).We considered the ΛCDM run and the Hu-Sawicki f (R) (Hu & Sawicki 2007) run with scalaron f R = −1 × 10 −4 ( f R4).By means of ∼21 000 and ∼28 000 dark matter halos in ΛCDM and f R4, respectively, arranged into three mass bins and 11 redshift bins, we showed that in both models ς decreases with increasing redshift at fixed mass and decreases with increasing mass (or equivalently, MAR) at fixed redshift.Specifically, in the ΛCDM model, ς decreases from ∼3R 200c to ∼2R 200c when z increases from 0 to 1.At z ∼ 0.1, ς decreases from 2.8R 200c to ∼2.5R 200c when the MAR increases from ∼10 4 h −1 M yr −1 to ∼2 × 10 5 h −1 M yr −1 .We computed the MAR with the recipe presented in Pizzardo et al. (2021Pizzardo et al. ( , 2022)), which, we demonstrated, agrees with results based on the merger trees of dark matter halos (Diemer et al. 2017).The behavior of ς in f R4 is similar to that of the ΛCDM, but at low redshifts, ς in f R4 is ∼15% larger on average than in ΛCDM, whereas at higher redshifts (z ∼ 0.8−1.0) the ςs of the two models overlap.In ΛCDM, the correlations of ς with redshift and MAR agree with the correlations expected for the splashback radius (Diemer & Kravtsov 2014;Adhikari et al. 2014;Deason et al. 2021;O'Neil et al. 2021;Baxter et al. 2021;Dacunha et al. 2022).
By using the median test, we find that when the dark matter halos with mass ∼10 14 h −1 M at z ≤ 0.8 are considered, ς can distinguish between the two theories of gravity with high significance, with a p-value 10 −5 .At z ≤ 0.5, the same halos give a p-value 10 −15 .At the highest redshift, ΛCDM and f R4 are degenerate, and ς cannot distinguish between the two models.Conversely, the number of halos with mass ∼10 15 h −1 M is too small to distinguish between the two models at any redshift: We showed that obtaining a p-value 10 −5 requires 400 halos.
Developing the full potential of the ς radius as a gravity probe thus requires further observational and computational efforts.New facilities such as the Vera Rubin Observatory (Ivezić et al. 2019) and Euclid (Laureijs et al. 2011;Sartoris et al. 2016) will provide extended weak-lensing mass profiles for thousands of clusters up to z ∼ 2. Large dense surveys obtained with advanced spectrographs on large telescopes such as WEAVE (Balcells et al. 2010;Dalton et al. 2012Dalton et al. , 2014;;Dalton 2016), PFS (Tamura et al. 2016), and eventually, MSE (Marshall et al. 2019) will provide dense spectroscopic surveys of a tremendous number of clusters up to z ∼ 0.5 that will complement the photometry.Weak gravitational lensing mass profiles estimated with sophisticated mass reconstruction methods (e.g., Umetsu et al. 2011;Umetsu 2013) combined with the caustic mass profiles (Diaferio & Geller 1997;Diaferio 1999;Serra et al. 2011;Pizzardo et al. 2023a) will improve the accuracy and the precision of the mass profiles at large distance from the cluster center.Larger simulated volumes will allow the identification of larger sets of the most massive halos, paving the way to tighter tests of structure formation models.

Fig. 1 .
Fig.1.Logarithmic derivative of the average density profiles, and R spl as function of mass, redshift, and theory of gravity.Upper panels: average profiles of d log ρ shell /d log r for the samples of simulated halos in Table1.The left, middle, and right panels show the average profiles for the low-, intermediate-, and high-mass samples, respectively, at the lowest and highest redshifts of each sample, as listed in the insets.The solid (dash-dotted) curves are the profiles in the ΛCDM ( f R4) model.The vertical solid (dash-dotted) thick black lines show the location of R spl at z = 0 in the ΛCDM ( f R4) model.The pair of vertical blue (green) solid and dash-dotted lines shows R spl at z = 1.0 (z = 0.3) for the ΛCDM and f R4 models, respectively.Bottom panels: the thin solid (dash-dotted) black curves show the radial derivative of the average mass profile in ΛCDM (fR4) at z = 0; the vertical thin solid (dash-dotted) black lines show ς avg of the low-, intermediate-, and high-mass samples at z = 0 in the ΛCDM ( f R4) models (Table2, first row, values in brackets).When extended to the upper panels, these thin solid and dash-dotted black lines cross the average density profiles at d log ρ shell /d log r −2, as expected from Eq. (1).

Fig. 2 .
Fig.2.Median differential radial mass profiles of the samples of the simulated halos in Table1.The top, middle, and bottom panels show the results for the low-, intermediate-, and high-mass samples, respectively.The solid (dash-dotted) curves are the median profiles for the ΛCDM ( f R4) halos.The profiles of the halos at different redshifts are shown with different colors, as listed in the insets.The shaded gray areas show the 68th percentile range of the distributions of the profiles of the individual halos at z = 0 in the ΛCDM model.The spreads of the profiles at the other redshifts and in the f R4 model are comparable.The vertical thick solid (dash-dotted) black line shows the median radius ς at z = 0 in the ΛCDM ( f R4) model.The vertical thick blue (green) lines show ς for the profiles at z = 1.0 (z = 0.3): the solid and dash-dotted lines show the ΛCDM and f R4 models, respectively.The dash-dotted blue line in the upper panel is overplotted on the solid blue line.The vertical stripes indicate the error bands of the median ςs.

PizzardoFig. 3 .
Fig. 3. Medians and associated errors of ς as a function of redshift in ΛCDM (blue points) and f R4 (orange points) for the low-(top panel), intermediate-(middle panel), and high-mass (bottom panel) samples.We applied small redshift offsets to the f R4 points to facilitate comparison.The black (red) lines are the fits of Eq. (4) to the ΛCDM ( f R4) points: specifically, the curves are the result of a spline interpolation of the values of ς (11 for the low-and intermediate-mass bin, 4 for the high-mass bin) resulting from Eq. (4) using the redshift, mass, and mass accretion rate of each simulated bin (Table1).In each panel, the lower plot shows the differences of the median ςs in ΛCDM and f R4 relative to the ΛCDM value.

Fig. 4 .
Fig. 4. MAR as a function of the cluster mass M 200c in the simulations and the MAR of the CIRS, HeCS, and HectoMAP clusters.The open circles show the median MARs of the 26 ΛCDM samples of Table 1, color-coded by redshift.The colored curves are the relation of Diemer et al. (2017) computed at the 11 redshifts of our simulated samples.The squares and the diamonds with error bars show the MAR estimates with their uncertainties of the stacked clusters of the CIRS and HeCS (Pizzardo et al. 2021) and HectoMAP clusters (Pizzardo et al. 2022), respectively.

z
Fig. 5. ς as a function of MAR, color-coded by redshift.The circles and crosses show the median values of the ΛCDM and f R4 samples, respectively.The solid (dash-dotted) curves are the fits of Eq. (4) to the ΛCDM ( f R4) data at three redshifts.From top to bottom, z = 0.0, 0.5, and 1.0.At fixed redshift, ς decreases with increasing MAR; at fixed MAR, ς decreases as the redshift increases.
Fig. 6. p-values of the median test to distinguish the two theories of gravity performed on the samples of ςs of the simulated clusters in Table 1.Top panel: blue points show the p-values of the test on the full samples of ςs in the low-mass bin, and the gray crosses show the p-values of the test performed by undersampling all the low-mass samples to 646 halos.Bottom panel: orange and green points show the pvalues of the test on the full samples in the intermediate-and high-mass bin, respectively; the gray crosses show the p-values of the test performed by undersampling the intermediate-mass samples to 11 halos.In both panels, the horizontal lines denote p = 10 −3 .

Table 1 .
Samples of simulated clusters.

Table 2 .
Radii R spl and ς avg of the simulated clusters.

Table 3 .
ς of the simulated clusters.