© The Authors 2024

1. 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 2020, 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 N-body simulations, and they showed that R_spl is usually located at ∼2R_200c¹ 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. 2014, 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. 2021; O’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. 2021, 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.

2. Simulations and samples of halos

2.1. 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 fR4 runs; the latter is an f(R) run with scalaron f_R0 = −1 × 10⁻⁴ 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₀ = 100 h km s⁻¹ Mpc⁻¹, with h = 0.6731, an initial perturbation amplitude 𝒜_s = 2.199 × 10⁻⁹, 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, σ₈, generally differs among the different simulations. In the two runs relevant for this work, ΛCDM and fR4, σ₈ = 0.847 and σ₈ = 0.967, respectively². The box size was 750 h⁻¹ Mpc on a side in comoving coordinates. The simulation contained 768³ collisionless dark matter particles with mass m_DM = 8.1 × 10¹⁰ h⁻¹ M_⊙ each. Because hydrodynamics is absent, the ΛCDM and the fR4 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 $\lambda =0.16\overline{d}$, where $\overline{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.

2.2. 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¹⁴ h⁻¹ M_⊙, and the intermediate- and high-mass samples contained all the halos with mass M_200c ∈ (4.0, 6.0)×10¹⁴ h⁻¹ M_⊙ and M_200c ∈ (9.0, 11.0)×10¹⁴ h⁻¹ 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 low-mass (L) halos, 11 samples of intermediate-mass (I) halos, and 4 samples of high-mass (H) halos. In total, the 26 ΛCDM (fR4) 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 fR4 is the result of the higher σ₈ 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 fR4, respectively), the intermediate-mass samples have ∼200 halos (154 and 251, averaging over all redshifts, for ΛCDM and fR4, respectively), and the high-mass samples have ∼35 halos (25 and 46, averaging over all redshifts, for ΛCDM and fR4, 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 intermediate-mass samples, this ratio is ∼28 − 38.

3. 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.

3.1. 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 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 ρ_bkg(z) = 3H²(z)Ω_M(z)/8πG, where Ω_M(z) = Ω_M(1 + z)³/[Ω_M(1+z)³+Ω_Λ], 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.

Fig. 1. Logarithmic derivative of the average density profiles, and Rspl 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 Table 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 (fR4) model. The vertical solid (dash-dotted) thick black lines show the location of Rspl at z = 0 in the ΛCDM (fR4) model. The pair of vertical blue (green) solid and dash-dotted lines shows Rspl at z = 1.0 (z = 0.3) for the ΛCDM and fR4 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 (fR4) 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 dlog ρshell/d log r ≃ −2, as expected from Eq. (1).

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 fR4: From z = 0 to z = 0.7, the relative differences between R_spl in fR4 and ΛCDM is in the range ∼2 − 15% in the low-mass bin, and ∼2 − 35% in the intermediate-mass 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 low-mass bin confirms this behavior: R_spl in fR4 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 fR4.

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.

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

$$\begin{array}{c}\hfill M^{″}\simeq \frac{\mathrm{d}}{\mathrm{d}r}\left(r^2{\rho }_{\mathrm{shell}}\right)=0⟺\frac{\mathrm{d}log{\rho }_{\mathrm{shell}}}{\mathrm{d}logr}=-2.\end{array}$$(1)

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²ξ_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 fR4 are larger than in ΛCDM by ∼10%, except in the low-mass sample at z = 1.0, where the ς_avg in fR4 is ∼2% smaller than in ΛCDM. Similarly, in all the high-mass samples, the ς_avg in fR4 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 fR4 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 fR4.

We assessed how well ς can distinguish between the ΛCDM and fR4 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 (fR4) 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).

Fig. 2. Median differential radial mass profiles of the samples of the simulated halos in Table 1. 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 (fR4) 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 fR4 model are comparable. The vertical thick solid (dash-dotted) black line shows the median radius ς at z = 0 in the ΛCDM (fR4) 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 fR4 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.

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 (fR4) 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 fR4. Figure 3 shows ς as a function of redshift and mass in ΛCDM and fR4.

Fig. 3. Medians and associated errors of ς as a function of redshift in ΛCDM (blue points) and fR4 (orange points) for the low- (top panel), intermediate- (middle panel), and high-mass (bottom panel) samples. We applied small redshift offsets to the fR4 points to facilitate comparison. The black (red) lines are the fits of Eq. (4) to the ΛCDM (fR4) 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 fR4 relative to the ΛCDM value.

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 fR4. As shown in Table 3, in the low- and high-mass bins, ς is systematically larger in fR4 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 fR4 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 fR4 at redshift z ∼ 0 − 0.5, where the relative difference decreases from ∼12% to ∼5%. In this mass bin, similarly to R_spl (Table 2), ΛCDM can present larger ςs than fR4 at z > 0.5.

3.2. 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 fR4 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 fR4.

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. (2021, 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 = δ_sR_i. The quantity δ_s is derived from the equation of motion of the infalling mass with initial velocity v_inf and constant acceleration −GM(< R_i)/(R_i + δ_sR_i/2)²,

$$\begin{array}{c}\hfill t_{\inf }^{2}GM(<R_{\mathrm{i}})-t_{\inf }2R_{\mathrm{i}}^2{(1+{\delta }_{\mathrm{s}}/2)}^2{v}_{\inf }-R_{\mathrm{i}}^3{\delta }_{\mathrm{s}}{(1+{\delta }_{\mathrm{s}}/2)}^2=0,\end{array}$$(2)

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. (2021, 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. 2021, 2022, for further details). When the mass of the infalling shell, M_shell, is known, the MAR is estimated with the equation

$$\begin{array}{c}\hfill \mathrm{MAR}\equiv \frac{M_{\mathrm{shell}}}{t_{\inf }}.\end{array}$$(3)

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).

Fig. 4. MAR as a function of the cluster mass M200c 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.

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

$$\begin{array}{c}\hfill \frac{\varsigma }{R_{200c}}(z,\text{MAR},M_{200c})=a[1+b{\mathrm{\Omega }}_{\mathrm{M}}(z)][1+cexp(-\frac{{\mathrm{\Gamma }}_{\mathrm{dyn}}}{d})],\end{array}$$(4)

where Ω_M(z) = Ω_M(1 + z)³/[Ω_M(1 + z)³ + Ω_Λ], Γ_dyn is the dynamical MAR ${\mathrm{\Gamma }}_{\mathrm{dyn}}\equiv \mathrm{d}log{M}_{\mathrm{vir}}/\mathrm{d}log\widehat{a}$, with $\widehat{a}={(1+z)}^{-1}$ the scale factor, and a, b, c, and d are four free parameters. Figure 4 proves that the difference between the M_virs 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 can therefore express the dynamical accretion rate as Γ_dyn ≈ [𝒦(t, Ω_M, H₀)/M_200c] (dM/dt), with dM/dt ≡ MAR³. 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 fR4 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 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 fR4 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 fR4 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 fR4 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⁻⁵ for the relation ς − z and for the relation ς–MAR.

Fig. 5. ς as a function of MAR, color-coded by redshift. The circles and crosses show the median values of the ΛCDM and fR4 samples, respectively. The solid (dash-dotted) curves are the fits of Eq. (4) to the ΛCDM (fR4) 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.

As expected, the fits show that the difference between the ΛCDM and fR4 models is larger at low redshifts, independently of the MAR. The difference steadily decreases with increasing redshift: consequently, because at low redshift ς is higher in fR4 than in ΛCDM, the fits in the (ς,z) plane at fixed MAR (Fig. 3) show steeper profiles in fR4 than in ΛCDM.

We now quantitatively assess to which extent ς can distinguish between ΛCDM and fR4. We applied the median test to each of our 26 ΛCDM–fR4 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 null-hypothesis that the samples are drawn from the same distribution. Generally, the fR4 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 fR4 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.

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 p-values 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.

We took p = 10⁻³ 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⁻³⁰). In contrast, only 3 of 11 intermediate-mass 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 fR4 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⁻³ 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 intermediate-mass bins, respectively. We therefore performed a median test in which we randomly undersampled all the samples of the low- and 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⁻³.

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 fR4 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.

4. 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. 2003, 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₂₀₀.

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 multi-object William Herschel Telescope Enhanced Area Velocity Explorer spectrograph on the WHT (WEAVE, Balcells et al. 2010; Dalton et al. 2012, 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 cluster-centric 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 fR4 runs of the DUSTGRAIN simulations that we used here were normalized at the CMB epoch. Therefore, their power-spectrum normalization at z = 0, σ₈ 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 σ₈, 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. (2021, 2022) in the ΛCDM L-CoDECS N-body simulation of Baldi (2012), which has σ₈ = 0.809. Despite the lower σ₈, 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, fR4, which means that if ς were inversely proportional to σ₈, we should expect ςs larger in ΛCDM than in fR4 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 σ₈, 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. 2023, 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 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¹⁵ h⁻¹ 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¹⁴ h⁻¹ M_⊙, can already provide observed samples that are sufficiently large to yield tight constraints on the theory of gravity.

5. 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¹⁴ h⁻¹ 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⁻⁴ (fR4). By means of ∼21 000 and ∼28 000 dark matter halos in ΛCDM and fR4, 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⁴ h⁻¹ M_⊙ yr⁻¹ to ∼2 × 10⁵ h⁻¹ M_⊙ yr⁻¹. We computed the MAR with the recipe presented in Pizzardo et al. (2021, 2022), which, we demonstrated, agrees with results based on the merger trees of dark matter halos (Diemer et al. 2017). The behavior of ς in fR4 is similar to that of the ΛCDM, but at low redshifts, ς in fR4 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¹⁴ h⁻¹ M_⊙ at z ≤ 0.8 are considered, ς can distinguish between the two theories of gravity with high significance, with a p-value ≲10⁻⁵. At z ≤ 0.5, the same halos give a p-value ≲10⁻¹⁵. At the highest redshift, ΛCDM and fR4 are degenerate, and ς cannot distinguish between the two models. Conversely, the number of halos with mass ∼10¹⁵ h⁻¹ M_⊙ is too small to distinguish between the two models at any redshift: We showed that obtaining a p-value ≲10⁻⁵ 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. 2012, 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.

Acknowledgments

We sincerely thank the referee for her/his insightful comments that greatly improved the quality of the paper. We thank Marco Baldi for providing snapshots of the ΛCDM and fR4 runs of the DUSTGRAIN N-body simulation. We thank Margaret Geller, Jubee Sohn, and Benedikt Diemer for numerous stimulating discussions. The graduate-student fellowship of M.P. was supported by the Italian Ministry of Education, University and Research (MIUR) under the Departments of Excellence grant L.232/2016. M.P. acknowledges the support of the Canada Research Chair Program and the Natural Sciences and Engineering Research Council of Canada (NSERC, funding reference number RGPIN-2018-05425). We also acknowledge partial support from the INFN grant InDark. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

References

Appendix A: The f(R) theory of gravity

In the f(R) model of gravity, the action S of General Relativity is modified by a scalar function, f(R), of the Ricci scalar curvature, R, as follows (Buchdahl 1970):

$$\begin{array}{c}\hfill S={\displaystyle \int }{\mathrm{d}}^{4}x\sqrt{-g}(\frac{R+f(R)}{16\pi G}+{\mathcal{L}}_m),\end{array}$$(A.1)

where g is the determinant of the metric, G is the universal gravitational constant, and ℒ_m is the matter Lagrangian. Hu & Sawicki (2007) proposed for f(R) the form

$$\begin{array}{c}\hfill f(R)=-m^2\frac{c_1{\left(\frac{R}{m^2}\right)}^n}{c_2{\left(\frac{R}{m^2}\right)}^n+1},\end{array}$$(A.2)

where c₁, c₂, and n are non-negative constant free parameters of the model, and m² is a mass scale, defined as $m^2\equiv \frac{8\pi G}{3}{\rho }_{M,0}=H_0^2\frac{{\rho }_{M,0}}{{\rho }_{crit,0}}=H_0^2{\mathrm{\Omega }}_M$. This specific form of f(R) guarantees (i) that at the recombination time, the model mimics ΛCDM, that is, lim_R → 0f(R) = 0; (ii) that at the late stages of cosmic evolution, the model encompasses cosmic acceleration through the presence of a cosmological constant-like quantity, namely lim_{R → ∞}f(R) = const; and (iii) the consistence with a large plethora of low-redshift phenomena through the tuning of the three free parameters.

In the Hu & Sawicki (2007) class of f(R) models, any departure from ΛCDM can only arise at low redshifts, or equivalently, in the high-curvature regime. This therefore is the only relevant limit for a comparison with ΛCDM. In the high-curvature limit, where $c_2{\left(\frac{R}{m^2}\right)}^n\to \infty$, we can expand Eq. (A.2) to first order in (R/m²) as

$$\begin{array}{c}\hfill f(R)\approx -m^2\frac{c_1}{c_2}+m^2\frac{c_1}{c_2^2}{\left(\frac{R}{m^2}\right)}^{-n}.\end{array}$$(A.3)

In this limit, if $c_1/c_2^2\to 0$ at fixed c₁/c₂, we see that f ≈ const, resembling the onset of a cosmological constant term.

Using Eq. (A.3), it is customary to introduce f_R,

$$\begin{array}{c}\hfill f_R=\frac{\mathit{df}}{\mathit{dR}}=-n\frac{c_1}{c_2^2}{\left(\frac{R}{m^2}\right)}^{-n-1},\end{array}$$(A.4)

a quantity that plays a central role in the description of the evolution of the space-time and of the dynamics of the matter fields in the f(R) model. f_R is useful for a simple characterization of the model: From Eqs. (A.3) and (A.4), for instance, we see that larger ns mimic ΛCDM until later times, whereas f(R) and ΛCDM become increasingly similar as $c_1/c_2^2$ decreases.

Equation (A.2) covers a wide variety of possible behaviors, which can be much different from ΛCDM. We are interested in obtaining realistic observations with small deviations from ΛCDM, however. We can pursue this aim by imposing conditions on the free parameters, in particular, by fixing a relation between c₁ and c₂.

Following Hu & Sawicki (2007), we start with the variation in the action S with respect to the metric to obtain the field equations

$$\begin{array}{c}\hfill G_{\mu \nu }+f_R{R}_{\mu \nu }-(\frac{f}{2}-\square f_R)g_{\mu \nu }-{\mathrm{\nabla }}_{\mu }{\mathrm{\nabla }}_{\nu }{f}_R=8\pi G{T}_{\mu \nu },\end{array}$$(A.5)

where $G_{\mu \nu }=R_{\mu \nu }-\frac{1}{2}R$, ∇_μ is the covariant derivative, and ▫ is the d’Alembert operator. The trace of Eq. (A.5) is

$$\begin{array}{c}\hfill 3\square f_R-R+f_{R}R-2f=-8\pi G\rho ,\end{array}$$(A.6)

which is the equation of motion of the field f_R. If ▫f_R is interpreted as a function of an effective potential V_eff,

$$\begin{array}{c}\hfill \square f_R=\frac{\partial V_{\mathrm{eff}}}{\partial f_R},\end{array}$$(A.7)

Eq. (A.6) yields

$$\begin{array}{c}\hfill \frac{\partial V_{\mathrm{eff}}}{\partial f_R}=\frac{1}{3}(R-f_{R}R+2f-8\pi G\rho ),\end{array}$$(A.8)

which states that V_eff has an extremum at

$$\begin{array}{c}\hfill R-f_{R}R+2f-8\pi G\rho =0.\end{array}$$(A.9)

We note that in the absence of any modification due to f(R), Eq. (A.9) reduces to the standard form of the Ricci scalar,

$$\begin{array}{c}\hfill R=8\pi G\rho .\end{array}$$(A.10)

At late times, when R/m² ≫ 1, Eq. (A.3) shows that the modification of gravity is approximately constant, holding $f(R)\approx -m^2\frac{c_1}{c_2}$, because the second term on the right-hand side subleads in (m²/R). Analogously, we note from Eq. (A.4) that f_RR ∝ (R/m²)⁻ⁿ → 0. This means that when it is approximated at leading order, Eq. (A.9) takes the following form:

$$\begin{array}{c}\hfill R=8\pi G\rho -2f\approx 8\pi G\rho +2\frac{c_1}{c_2}{m}^2.\end{array}$$(A.11)

In ΛCDM, when we consider the Friedmann equations with the explicit presence of the cosmological constant Λc² = 8πGρ_Λ, the ΛCDM Ricci scalar (Eq. A.10) can be rewritten in the form

$$\begin{array}{c}\hfill R=8\pi G({\rho }_M+4{\rho }_{\mathrm{\Lambda }}).\end{array}$$(A.12)

Small deviations from the ΛCDM expectations are achieved when we impose equality between Eqs. (A.11) and (A.12) because in this way, any discrepancy is generated only by the subleading terms of f(R). This equality translates into the equality between the two second terms on the right-hand side of Eqs. (A.11) and (A.12) when we set ρ ≡ ρ_M in Eq. (A.11). We obtain the following constraint between c₁ and c₂:

$$\begin{array}{c}\hfill \frac{c_1}{c_2}=6\frac{{\mathrm{\Omega }}_{\mathrm{\Lambda }}}{{\mathrm{\Omega }}_M}.\end{array}$$(A.13)

By taking this contraint into account and by considering a pure stiff-matter model, we can finally write Eq. (A.11) as

$$\begin{array}{c}\hfill R\approx 8\pi G{\rho }_M+12\frac{{\mathrm{\Omega }}_{\mathrm{\Lambda }}}{{\mathrm{\Omega }}_M}{m}^2=3m^2(a^{-3}+4\frac{{\mathrm{\Omega }}_{\mathrm{\Lambda }}}{{\mathrm{\Omega }}_M}),\end{array}$$(A.14)

whose current value is

$$\begin{array}{c}\hfill R_0=3m^2(1+4\frac{1-{\mathrm{\Omega }}_M}{{\mathrm{\Omega }}_M})=m^2(\frac{12}{{\mathrm{\Omega }}_M}-9).\end{array}$$(A.15)

When n is fixed, the model is left with only one free parameter: c₂, which is more commonly rephrased in terms of the current value of f_R, f_R0, known as the scalaron. Its expression follows from Eqs. (A.4)-(A.13)-(A.15). When we take n = 1, the chosen value in the DUSTGRAIN suite, the scalaron is

$$\begin{array}{c}\hfill f_{R0}=-\frac{1}{c_2}\frac{6{\mathrm{\Omega }}_{\mathrm{\Lambda }}}{{\mathrm{\Omega }}_M}{\left(\frac{m^2}{R_0}\right)}^2.\end{array}$$(A.16)

For the purposes of our work, we considered the DUSTGRAIN f(R) run with f_R0 = −1 ⋅ 10⁻⁴.

All Tables

All Figures

Fig. 1. Logarithmic derivative of the average density profiles, and Rspl 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 Table 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 (fR4) model. The vertical solid (dash-dotted) thick black lines show the location of Rspl at z = 0 in the ΛCDM (fR4) model. The pair of vertical blue (green) solid and dash-dotted lines shows Rspl at z = 1.0 (z = 0.3) for the ΛCDM and fR4 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 (fR4) 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 dlog ρshell/d log r ≃ −2, as expected from Eq. (1).

In the text

Fig. 2. Median differential radial mass profiles of the samples of the simulated halos in Table 1. 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 (fR4) 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 fR4 model are comparable. The vertical thick solid (dash-dotted) black line shows the median radius ς at z = 0 in the ΛCDM (fR4) 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 fR4 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.

In the text

Fig. 3. Medians and associated errors of ς as a function of redshift in ΛCDM (blue points) and fR4 (orange points) for the low- (top panel), intermediate- (middle panel), and high-mass (bottom panel) samples. We applied small redshift offsets to the fR4 points to facilitate comparison. The black (red) lines are the fits of Eq. (4) to the ΛCDM (fR4) 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 fR4 relative to the ΛCDM value.

In the text

Fig. 4. MAR as a function of the cluster mass M200c 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.

In the text

	Fig. 5. ς as a function of MAR, color-coded by redshift. The circles and crosses show the median values of the ΛCDM and fR4 samples, respectively. The solid (dash-dotted) curves are the fits of Eq. (4) to the ΛCDM (fR4) 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.
In the text

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 p-values 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.

In the text