© ESO 2021

1. Introduction

In the current cold dark matter (CDM) model of the formation and evolution of cosmic structures, smaller dark matter halos hierarchically aggregate into larger and more massive halos. The accretion occurs through mergers with halos of comparable (major mergers) or lower (minor mergers) mass and with the capture of diffuse dark matter particles (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).

In principle, the mass accretion of dark matter halos is a valuable tool for testing different models of structure formation. Specifically, the mass evolution M(z) of a dark matter halo, which describes its mass assembling history, or its time derivative, the mass accretion rate Ṁ(z) can be used to determine the halo formation redshift (Lacey & Cole 1993; van den Bosch 2002; Ragone-Figueroa et al. 2010; Giocoli et al. 2012). The rates are correlated with halo properties including concentration (Wechsler et al. 2002; Tasitsiomi et al. 2004; Zhao et al. 2009; Giocoli et al. 2012; Ludlow et al. 2013), shape (Kasun & Evrard 2005; Allgood et al. 2006; Bett et al. 2007; Ragone-Figueroa et al. 2010), spin (Vitvitska et al. 2002; Bett et al. 2007), degree of internal relaxation (Power et al. 2011), and fraction of substructures (Gao et al. 2004; van den Bosch et al. 2005; Ludlow et al. 2013). The assembly of dark matter halos can trace the accretion rate of baryons from the cosmic web onto the dark matter halo (van den Bosch 2002; McBride et al. 2009; Fakhouri et al. 2010; Wright et al. 2020).

In theoretical studies, the estimates of the mass accretion history (MAH) M(z) and the mass accretion rate (MAR) Ṁ(z) of a dark matter halo evolved at the present time, z = 0, are usually tackled by tracing the merger tree of the halo, either with numerical simulations (e.g. Genel et al. 2008; Kuhlen et al. 2012; Baldi 2012) or with Monte Carlo methods (Lacey & Cole 1993; Kauffmann et al. 1993; Somerville & Kolatt 1999; Parkinson et al. 2007; Jiang & van den Bosch 2014). In hierarchical clustering scenarios, when we follow the growth history of the dark matter halo backwards in time, we see that the halo separates into two or more halos; at each epoch, the growing halo, known as the descendant, has a main progenitor in the form of the most massive halo, which merges with smaller halos and thereby generates the descendant. By identifying the main progenitor of the halo at each epoch, we can trace the formation history, or MAH, of a simulated halo by tracing the main branch of its merger tree.

Within this framework, numerous studies investigate the MAH in ΛCDM cosmologies. Analytical approximations can describe the MAH as a function of the final mass of the halo and other additional parameters. For example, McBride et al. (2009), Fakhouri et al. (2010), and Correa et al. (2015) adopt the relation M(z) = M₀(1 + z)^βe^−γz for the growth in mass with redshift z, whereas van den Bosch (2002) and van den Bosch et al. (2014) adopt . In the latter case, the formation redshift, z_f, and ν are left as free parameters, whereas in the former expression, β and γ are left either free or fixed by the linear power spectrum of matter (Correa et al. 2015). In both formulae, M₀ is the halo mass at redshift z = 0.

These studies point out that the MAH can be separated into two regimes. In the first regime, at early times, the mass accretion is relatively large and the growth is nearly exponential in redshift; here, major mergers are frequent and keep the system unrelaxed. In the second regime, at later times, the accretion slows down and the growth is governed by a power-law in redshift, thus enabling the halo to reach virial equilibrium. These studies also show that more massive halos, which form at relatively low redshifts, have a greater MAR than less massive halos. This correlation is supported by (1) the correlation between the age and the concentration and (2) the anti-correlation between the mass and the concentration of dark matter halos (Zhao et al. 2009). In other words, old, low-mass, and highly concentrated dark matter halos are expected to have lower MAR than young, high-mass and less concentrated halos.

The MAH can be a probe of the cosmological parameters. Hurier (2019) uses the thermal Sunyaev-Zel’dovich (SZ) effect as a proxy for the mass of the clusters from the second Planck SZ catalogue (Ade et al. 2016) and the fit by Correa et al. (2015) to the MAH of dark matter halos in simulations to derive values for the power spectrum normalisation σ₈, the cosmic mass density Ω_m, and the Hubble parameter H₀, which yield σ₈(Ω_m/0.3)^−0.3(H₀/70)^−0.2 = 0.75 ± 0.06. This value is in rough agreement with other analyses of galaxy cluster samples and of the power spectrum of the cosmic microwave background (de Haan et al. 2016; Planck Collaboration VI 2020; Zubeldia & Challinor 2019).

The investigation of the MAH and the MAR is also connected to the splashback radius, which roughly corresponds to the first apocentre of the orbits of recently accreted material. This radius is larger than the radii R₂₀₀ or R_vir¹ that are usually adopted to quantify the halo size and is close to ∼2R₂₀₀, on average (More et al. 2015). In their simulations, Diemer & Kravtsov (2014) find that the steepness of the slope of the halo mass profile at large radii increases with increasing MAR. Moreover, the cluster-centric radius of this change of slope decreases with increasing MAR. Adhikari et al. (2014) associate the location of this feature with the splashback radius.

A change of slope that is consistent with the expectations from the simulations is indeed present in the profile of the surface number density of galaxies from the Dark Energy Survey (DES) cross-correlated with the SZ clusters from the South Pole Telescope (SPT) and the Atacama Cosmology Telescope (ACT) (Shin et al. 2019), as well as in the deprojected cross-correlation of the SZ clusters from the Planck Survey with galaxies detected photometrically in the PanSTARRS survey (Zürcher & More 2019). Similarly, the splashback radius is detected in the inferred dark matter density profiles of the redMaPPer clusters (More et al. 2016) and in clusters from either the Sloan Digital Sky Survey (SDSS; Baxter et al. 2017) or DES (Chang et al. 2018).

Although interlopers along the line of sight affect the inference of this feature both from optically selected clusters (Busch & White 2017; Shin et al. 2019; Sunayama & More 2019) and from weak lensing analyses of X-ray selected clusters (Umetsu & Diemer 2017; Contigiani et al. 2019), dense galaxy redshift surveys (Geller et al. 2011; Serra & Diaferio 2013; Sohn et al. 2018; Rines et al. 2018), and upcoming lensing surveys might potentially overcome the effects of this contamination and constrain the relation between the splashback radius and the accretion rate (Xhakaj et al. 2020).

All of these studies highlight the relevance of an observational estimate of the MAR of galaxy clusters. Unfortunately, only a handful of direct measurements have been attempted thus far. In fact, we can observe a real cluster only at a specific time and we cannot clearly identify its merger tree to quantify the MAR, as we would usually do when N-body simulations are available. A viable approach that has been pursued with real clusters is to identify galaxy groups that surround the cluster and appear to fall onto it (e.g. Rines et al. 2001, 2002); unfortunately, the estimate of their masses does not provide an estimate of the MAR, offering only a lower limit instead (e.g. Lemze et al. 2013; Haines et al. 2018).

More importantly, the cluster outer region needs to be chosen properly. For example, Lemze et al. (2013) investigate the region slightly beyond R₂₀₀ in the X-ray and optical bands, whereas Tchernin et al. (2016) detect infalling gas clumps of A2142 in X-ray and SZ out to ∼1.3R₂₀₀. Similarly, Haines et al. (2018) identify the infalling groups in the range (0.28−1.35)R₂₀₀. According to studies of the splashback radius, these radii may be too small to return a full estimate of the MAR: in fact, we expect that these regions contain matter with rather different dynamical histories: matter that is falling onto the cluster for the first time, matter that is moving outwards, and matter that is falling back again.

To avoid this complex dynamical structure and to return a proper estimate of the MAR, we must consider a region that is further out, beyond, at least, the splash-back radius of ∼2R₂₀₀. We expect that this region mostly contains matter that will fall onto the cluster in the near future (De Boni et al. 2016; Haines et al. 2018) and that the fraction of matter that is actually moving out from the central region is limited (Ludlow 2009; Diemer et al. 2017; Xhakaj et al. 2020; Bakels et al. 2021; Aung et al. 2021).

In this paper, we pursue this idea that was originally suggested by De Boni et al. (2016). We estimate the MAR from the amount of mass in the cluster’s outer region beyond ∼2R₂₀₀. Unfortunately, compared to the cluster central region, the cluster outskirts are a large and low-density region, where the system is not dynamically relaxed. The methods used to estimate the mass cannot rely on the hypothesis of virial equilibrium. The caustic method and weak gravitational lensing are two complementary methods that do not rely on this hypothesis (Geller et al. 2013) and are thus appropriate to estimate the amount of mass in these outer regions.

Here, we present the first estimation of the MAR of real clusters based on the spherical accretion model developed by De Boni et al. (2016). To estimate the cluster mass profiles in their outer region, we use the caustic technique (Diaferio & Geller 1997; Diaferio 1999). The caustic technique is known to return an unbiased mass estimate with a relative uncertainty of 50%, or better, in the regions where accretion is taking place (Serra et al. 2011). Unlike methods based on weak gravitational lensing, where the signal is stronger at intermediate redshift and rapidly drops at low and high redshift (Hoekstra 2003; Hoekstra et al. 2011), the caustic technique can be applied to clusters at any redshift, provided that the number of cluster galaxies with spectroscopic redshift is large enough to sample the velocity field properly.

In Sect. 2, we briefly summarise the spherical infall method introduced in De Boni et al. (2016). In Sect. 3, we use N-body simulations to test our method of estimating the MAR. Section 4 describes the CIRS and HeCS catalogues of real clusters that span the redshift range 0−0.3 and the mass range ∼10¹⁴ − 10¹⁵ h⁻¹ M_⊙. In Sect. 5, we illustrate and discuss the estimates of the MAR of individual clusters and the mean MAR of the cluster samples as a function of mass and redshift. We discuss our results in Sect. 6 and we present our conclusions in Sect. 7. We adopt H₀ = 100 h km s⁻¹Mpc⁻¹ throughout.

2. Spherical accretion recipe

In this section, we briefly review the model proposed in De Boni et al. (2016) to evaluate the MAR of clusters from spectroscopic redshift surveys.

De Boni et al. (2016) estimate the MAR from the merger trees of dark matter halos, of a mass of ∼10¹⁴ h⁻¹ M_⊙, extracted from N-body simulations. They find that in the redshift range of z = [0, 2], a spherical accretion recipe returns an unbiased MAR within ∼20% of the average MAR from the merger trees. For the more massive halos of a mass of ∼10¹⁵ h⁻¹ M_⊙, even though the statistics of De Boni et al. (2016) are rather poor and the MAR estimated with the spherical accretion recipe is ∼10−40% biased towards the low end, the recipe still returns a MAR within the 1-σ spread of the MAR derived from the merger trees.

In the spherical accretion recipe, the MAR is estimated by assuming that the mass within a spherical shell of a proper thickness, centred on the cluster, will fully accrete onto the cluster within a given time interval, t_inf. We assume that the massive shell falls onto the cluster with constant acceleration and a given initial proper velocity v_i. By solving the equation of motion, we end up with an equation whose unknown is the thickness δ_sR_i of the shell, where R_i is the physical inner radius of the shell:

(1)

where M(< R_i) is the mass of the cluster within the radius R_i, and G is the gravitational constant.

We assume that the inner radius of the shell is R_i = 2R₂₀₀. As anticipated in the introduction, this radius, R_i, is close to the average splashback radius of massive dark matter halos of cluster size at redshift z < 2 found in the N-body simulations (More et al. 2015). This radius also approximates the inner radius of the region containing the mass that will fall onto clusters in the near future. In addition, it is close to R_i that the absolute value of the infall radial velocity generally reaches its maximum (see Sect. 3.3) and the recipe thus includes the largest contribution to the MAR. The solution of Eq. (1) yields the shell thickness, δ_sR_i, and the mass of the shell, M_shell, if the mass profile of the cluster at radii larger than 2R₂₀₀ is known.

The MAR is thus simply estimated as:

(2)

Our analysis below shows that Eq. (1) typically yields δ_sR_i ∼ 0.5R₂₀₀ and, thus, M_shell is the mass of the shell of inner and outer radii 2R₂₀₀ and ∼2.5R₂₀₀, respectively. The N-body simulations suggest that a fraction of this M_shell has already been within R₂₀₀ and is thus not actually falling onto the central region of the cluster for the first time (e.g. Aung et al. 2021). According to Ludlow (2009) and Bakels et al. (2021), no more than 30% of the subhalos in the radial range 2−3R₂₀₀ have already been within R₂₀₀. Similarly, Diemer et al. (2017) and Xhakaj et al. (2020) find that 75% to 87% of the apocenters of the dark matter particles on their first orbit agree with the splashback radius identified by the maximum slope of the logarithmic derivative of the halo density profile (Diemer & Kravtsov 2014; More et al. 2015). This feature generally is within 2R₂₀₀ (Diemer & Kravtsov 2014; Xhakaj et al. 2020). We can thus conclude that less than ∼20−25% of the dark matter particles in the radial range ∼2−2.5R₂₀₀ are not falling onto the inner region of radius 2R₂₀₀ for the first time.

In principle, we could introduce a correction factor α_M < 1 for M_shell in Eq. (1), to take this effect into account. However, avoiding the introduction of this parameter has the advantage of keeping the recipe as simple as possible without affecting the effectiveness of our method: indeed, applying exactly the same recipe to real and simulated clusters yields a sensible comparison between the observed and the expected MAR. Alternatively, we could make the introduction of α_M irrelevant by adopting the inner radius of the shell larger than 2R₂₀₀ and thus decrease the fraction of M_shell that is not actually falling in for the first time; however, this choice would return a spherical shell that extends to radii larger than ∼2.5−3R₂₀₀. As we discuss in Sect. 4.3, in these regions, real clusters currently suffer from spectroscopic incompleteness and, consequently, the estimate of the MAR would be biased.

Our recipe for the MAR estimate still requires the values of t_inf and v_i, that are not currently measurable in real clusters. For t_inf and v_i, we resort to N-body simulations of dark matter halos of cluster size, assuming that these systems resemble real clusters. Clearly, these values can differ widely from halo to halo. To apply Eq. (2) properly, we adopt an average value for both t_inf and v_i for samples of dark matter halos within proper mass and redshift bins. If the ranges covered by these bins are sufficiently small, as we detail below, this approach can provide the MAR of a real cluster if its mass profile is known.

We measure the mass profile of real clusters to radii larger than R₂₀₀ with the caustic technique. This technique uses redshift data alone and does not require the assumption of dynamical equilibrium; however, it assumes spherical symmetry and deviations from this symmetry are responsible for most of the uncertainty in the mass profile and, consequently, in the MAR.

3. Testing the mass accretion recipe on mock redshift surveys of clusters

Before applying the MAR recipe to real clusters, we need to evaluate the reliability of our MAR estimate and its possible systematic errors. Here, we test the recipe on mock redshift surveys of clusters extracted from an N-body simulation of a ΛCDM model. We also use this simulation to provide the proper values for t_inf and v_i for clusters in different bins of mass and redshift.

In Sects. 3.1 and 3.2, we describe the N-body simulation and the construction of the mock redshift surveys, respectively. In Sect. 3.3 we discuss the radial velocity profiles of the clusters in the simulation and in Sect. 3.4, we illustrate how our estimate of the MAR from the redshift surveys compares with the MAR estimated with the full three-dimensional information.

3.1. CoDECS simulations

For our tests, we relied on the L-CoDECS simulations (Baldi 2012), which is a set of N-body numerical simulations of a ΛCDM cosmology and other quintessence models. For our purposes, we used only the ΛCDM run.

The simulations are normalised at the cosmic microwave background epoch with cosmological parameters according to the WMAP7 analysis (Komatsu et al. 2011): cosmological dark matter density Ω_m0 = 0.226, cosmological constant Ω_Λ0 = 0.729, baryonic mass density Ω_b0 = 0.0451, Hubble constant H₀ = 70.3 km s⁻¹ Mpc⁻¹, power spectrum normalisation σ₈ = 0.809, and power spectrum index n_s = 0.966. The box size is 1 h⁻¹ Gpc on a side in comoving coordinates. The simulation contains 1024³ dark matter particles with mass m_DM = 5.84 × 10¹⁰ h⁻¹ M_⊙ and the same number of baryonic particles with mass m_b = 1.17 × 10¹⁰ h⁻¹ M_⊙. Baryons are included only to check fifth-force effects in the quintessence cosmologies, but no hydrodynamics is included in the simulation. The ΛCDM run is a standard collisionless N-body simulation².

Groups and clusters in the simulations were identified with a friends-of-friends algorithm with linking length , where is the mean Lagrangian inter-particle separation. The algorithm is only run over the dark matter particles. Each baryonic particle is associated with the closest dark matter particle at the end of the procedure.

We considered only the clusters in two different mass bins with median mass, within R₂₀₀ at z = 0, M₂₀₀ ≃ 10¹⁴ h⁻¹ M_⊙, and M₂₀₀ ≃ 10¹⁵ h⁻¹ M_⊙. The two samples contain N = 2000 and N = 50 clusters, respectively. We identify the main progenitors of these clusters and consider the samples of these progenitors at each redshift. We consider the outputs of the simulation at six different redshifts: z = 0.0, 0.12, 0.19, 0.26, 0.35, 0.44. Table 1 shows the medians and the 68th percentile ranges of the two mass bins at different redshifts.

3.2. Mock catalogues and mass profiles

The first ingredient of the MAR recipe is the three-dimensional mass profile of the cluster extending to large radii. We use the caustic method to estimate this profile in real clusters. With a sufficiently dense redshift survey of the outer regions of an individual cluster, the caustic mass profile deviates from the real three-dimensional mass profile by less than 10%, on average, in the radial range [0.6, 4]R₂₀₀, with a 1σ relative uncertainty of ∼50% (Serra et al. 2011).

The uncertainty in the caustic mass profile is due mostly to the assumption of spherical symmetry and clearly propagates into our estimate of the MAR. To quantify how this uncertainty propagates, we apply our MAR recipe to synthetic observations of the simulated clusters. We assume that the dark matter halos we model in our N-body simulation are a realistic representation of galaxy clusters, and that their dark matter particles trace the same velocity field of the real galaxies. This latter assumption is supported by N-body/hydrodynamical simulations that suggest that the velocity bias, b_v, between the velocity dispersions of galaxies and dark matter particles is negligible in the outskirts of galaxy clusters (Diemand et al. 2004; Hellwing et al. 2016; Armitage et al. 2018). Specifically, Armitage et al. (2018) show that in their 30 dark matter halos with mass M₂₀₀ in the range (7 × 10¹³ − 1.7 × 10¹⁵) h⁻¹ M_⊙, b_v steadily decreases from ∼1.4 in the halo centre to ∼1 beyond 2R₂₀₀; this is mainly because dynamical friction, which dominates in the central region (e.g. Tormen et al. 1998; Taylor & Babul 2004), is less relevant in the outer regions.

To create the mock redshift survey of a simulated cluster, we extract a squared-basis truncated pyramid centred on the cluster, with the smaller basis closer to the observer (see Fig. 1); the pyramid axis is aligned along one of the three cartesian coordinates chosen as the line of sight and has a height of 2b_L = 140 h⁻¹ Mpc. The sizes, r_min and r_max, of the two bases are defined by the intercept theorem r_s : r_FOV/2 = (r_s∓b_L) : r_min, max, where r_s is the comoving distance of the cluster centre from the observer and r_FOV is the size of the field of view (FOV). We can use the previous relation that is appropriate for Euclidean geometry because we are considering the flat ΛCDM model. Similarly to Serra et al. (2011), we take r_FOV = 12 h⁻¹ Mpc. This value easily covers even the most massive clusters out to a few virial radii.

Fig. 1. Schematic figure of the truncated-pyramidal volume of the mock catalogue of a cluster extracted from the simulation. The solid dots show the positions of the dark matter particles in a slice of the simulation box centred on the cluster C. The orange (blue) dots are within (outside) the volume. This figure has only an illustrative purpose: the actual mock catalogues have bL ≫ rFOV and include a substantially larger number of particles than shown here.

For each cluster, we built three mock catalogues, one for each cartesian coordinate chosen as the line of sight. In general, the clusters are not spherically symmetric. Thus, to improve the statistics, we can consider these three mock catalogues as independent clusters even though they are not completely independent systems.

Each selected volume contains the number and of dark matter particles for the low- and high-mass bins, respectively; the ranges shown indicate the 10th − 90th percentile ranges. In the same field of view, the densest survey of a real cluster contains a few thousand galaxies (e.g. Hwang et al. 2014; Sohn et al. 2019). Therefore, to identify each dark matter particle with an individual galaxy, we randomly sample a limited number of dark matter particles within the volume. According to Serra et al. (2011), the caustic method performs better when the velocity field of the cluster is sampled by ∼200 galaxies within 3R₂₀₀ from the cluster centre. We follow this approach and sample the dark matter particles until we reach N_sample = 200 particles within 3R₂₀₀. This procedure yields the number and of particles within each mock redshift survey for the low- and the high-mass bin, respectively.

The observed redshift z of each particle is set by its cosmological redshift and its peculiar velocity v_p = u_p/(1 + z_s) in the comoving frame of the simulation box: z_s is the redshift of the simulation snapshot and u_p is the comoving peculiar velocity provided by the N-body simulation³.

The comoving distance from the observer to the centre of the simulation box, which coincides with the cluster centre, is , where E(z) = H(z)/H₀ = [(Ω_m0 + Ω_b0)(1 + z)³ + Ω_Λ]^1/2 in the flat ΛCDM model. The particle position vector in the observer reference frame is thus r_i = r_s + r_c,i, where r_c, i is the comoving position vector of the particle in the reference frame of the simulation box. This sum of vectors is derived in the Euclidean geometry of the ΛCDM model of the simulation.

The cosmological redshift, z_i, of the particle satisfies the relation and the observed component of the peculiar velocity is v_los,i = v_p · r_i/r_i. The redshift due to the peculiar velocity is z_p, i = v_los, i/c because v_los, i ≪ c. Combining z_p, i with z_i yields the observed redshift (1 + z_obs, i) = (1 + z_i)(1 + z_p, i), namely:

(3)

Standard geometrical transformations finally yield the celestial coordinates (α, δ) from the cartesian components of r_i. For all the clusters, we chose the celestial coordinates of the cluster centre (α, δ) = (6^h, 0). For the snapshot corresponding to z = 0, we located the cluster centre at cz = 32 000 km s⁻¹, similarly to the z = 0 mock catalogues described in Serra et al. (2011).

We have two samples of simulated clusters: 2000 clusters with M₂₀₀ ∼ 10¹⁴ h⁻¹ M_⊙ and 50 clusters with M₂₀₀ ∼ 10¹⁵ h⁻¹ M_⊙ at z = 0. By projecting each clusters along three lines of sight, we obtain two samples of 6000 and 150 cluster redshift surveys for the low- and high-mass bins, respectively.

For each mock catalogue, we construct the R − v_los diagram, the line-of-sight velocity relative to the cluster mean as a function of the projected distance from the cluster centre. The caustic method returns the mass profile estimated from the amplitude of the caustics (see Diaferio 1999; Serra et al. 2011, for further details on the caustic technique). Figures 2 and 3 show some examples of the procedure for the low- and high-mass bins. In both figures, the left and right columns show simulated clusters at z = 0.12 and z = 0.19, respectively; the upper and lower panels show the R − v_los diagrams and the mass profiles.

Fig. 2. R − vlos diagram (top panels) and mass profile (bottom panels) of two simulated clusters in the low-mass bin. Left (right) column: cluster at z = 0.12 (z = 0.19). In the bottom panels, the red (blue) solid curves show the caustic (real) mass profile, whereas the red (blue) dot-dashed curves show the NFW fits to the caustic (real) mass profile. In these two examples, the two NFW fits are indistinguishable. The shaded areas show the 50% confidence level of the caustic location and of the caustic mass profile according to the caustic technique recipe. In the R − vlos diagrams the shaded areas are present, but very thin.

Fig. 3. Same as Fig. 2 but for two simulated clusters in the high-mass bin.

The mass profiles estimated from the caustic amplitude are within the expected uncertainty from the real mass profile within ∼4R₂₀₀. At larger radii, the caustic mass profiles generally flattens, unlike the real mass profile. These regions can include large nearby groups and clusters that would increase the caustic amplitude with increasing radius; however, the caustic technique is designed to favour decreasing amplitudes with increasing radius, to avoid the inclusion of background and foreground systems, as detailed in Diaferio (1999) and Serra et al. (2011). Thus, at radii where the amplitude would increase by an anomalous amount, the algorithm sets the caustic amplitude to zero and the cumulative mass profile gets flattened.

We fit both the real and caustic mass profiles with the Navarro et al. (1997, NFW) model. We fit the caustic mass profiles up to 4R₂₀₀, unless the caustic amplitude shrinks to zero at a smaller radius. Similarly, we fit the real mass profiles up to 4R₂₀₀. In the bottom panels of Figs. 2 and 3, the red and blue dot-dashed lines show the NFW fits to the caustic and the real mass profiles, respectively. When averaged over the entire radial range, the NFW fits overestimate, by ∼5%, on average, the caustic and the real mass profile. This analysis shows that the NFW profile properly models the mass profiles of the clusters in our simulation to radii substantially larger than R₂₀₀. This result is also observed in real galaxy clusters in the CAIRNS (Rines et al. 2003), CIRS (Rines & Diaferio 2006), and HeCS catalogues (Rines et al. 2013). Indeed, by applying the caustic technique to a dense redshift survey of the Coma cluster, Geller et al. (1999) were the first to prove that the NFW model describes the mass profile of real clusters to these large radii; a few years later, this result was confirmed by weak gravitational lensing observations of other clusters (Clowe & Schneider 2001; Kneib et al. 2003; Lemze et al. 2009).

The concentration parameters derived from the NFW fits to the real mass profiles are in the range between 2−6 for both mass bins, which is consistent with the concentration parameters of dark matter halos with similar mass and redshift in other N-body simulations (e.g. Ludlow et al. 2014, 2016). The concentrations derived from the fits to the caustic mass profiles are, on average, 15% larger than those derived from the real mass profiles. This overestimation confirms the results of Serra et al. (2011) and is derived from the fact that the caustic method generally overestimates the mass of the cluster within ∼0.6R₂₀₀.

3.3. Radial velocity profiles and the accretion time, t_inf

The second ingredient of our MAR recipe is the radial velocity profile of the cluster extending to large radii. This piece of information cannot be derived from observations. Therefore, we need to rely on the modelling of clusters in the N-body simulations.

For each cluster, we construct the radial velocity profile by computing the mean radial velocity v_rad of the particles within individual radial bins. The radial velocity of each particle is v_i = [v_p + H(z_s)a(z_s)r_c,i] · r_c,i/r_c,i, where, as in the previous section, v_p = u_p/(1 + z_s) is the proper peculiar velocity and r_c,i is the comoving position vector of the particle from the cluster centre; H(z_s) and a(z_s) are the Hubble parameter and the cosmic scale factor at the snapshot redshift z_s, respectively.

Figure 4 shows the mean and median profiles of the radial velocity of the simulated clusters in our low- and high-mass samples at z = 0 and z = 0.44. We also show the 68th percentile range of the profile distribution. In the estimate of the velocity profiles, we include both the dark matter and the baryonic particles. The simulated clusters at other redshifts show qualitatively similar results.

Fig. 4. Mean (blue) and median (orange) profiles of the radial velocity of the particles within the dark matter halos extracted from the simulation for the two mass bins at z = 0 and z = 0.44, as indicated in the panels. 68% of the profiles of the individual halos lie within the light blue areas.

By inspecting the radial velocity profile, we can identify three regions: an internal region with radial velocity v_rad ≃ 0, where the matter is moving close to the centre within an isotropic velocity field; an infall region, where v_rad becomes negative and indicates an actual infall of matter towards the centre of the cluster; and a Hubble region at radii r ≳ 4R₂₀₀, where v_rad becomes positive and the Hubble flow dominates. Broadly speaking, the infall radius R_inf, defined as the radius where the minimum of v_rad occurs, is between 2R₂₀₀ and 3R₂₀₀, independently of the cluster mass and redshift.

We estimate the mean radial velocity profile at discrete values of the radius r. For the infall velocity v_i of the shell (Eq. (1)), we take the radial velocity associated with the bin [2−2.5]R₂₀₀. The choice of this velocity is consistent with the shell adopted for the estimate of the MAR: the internal radius of the shell is R_i = 2R₂₀₀, comparable with the splashback radius (More et al. 2015; Adhikari et al. 2014), and the thickness typically returned by Eq. (1) is ∼0.5R₂₀₀. For both bins of cluster mass, this radial bin corresponds to the radial range where the radial velocity profile has its minimum thus capturing the largest contribution to the MAR.

Our prescription for the MAR estimate clearly depends on the choice of both R_i and v_i. We will discuss the effect of the value of v_i on the MAR in Sect. 5.2 below. Here, we show that using the mean infall velocity rather than the infall velocity of each cluster has very little impact on the estimate of the MAR. We compare the spreads of the velocity profiles shown in Fig. 4 with the spreads of the MAR estimates obtained with the mean v_i and either the three-dimensional mass profiles or the caustic mass profiles that we compute in Sect. 3.4 below (and shown in Fig. 5).

Fig. 5. MAR of simulated clusters in the low- (lower set of points) and high- (upper set of points) mass bins. The blue squares and the red triangles show the median MAR based on the three-dimensional and the caustic mass profiles, respectively. The blue shaded areas show the 68th percentile range of the distribution of the MAR derived from the three-dimensional mass profiles; the red error bars show the 68th percentile ranges of the estimates obtained with the caustic mass profiles.

We use the 1σ relative standard deviation of the mean , where is the mean of a sample of N measurements of the quantity X. For the low-mass bin, the infall velocity standard deviation ϵ_vi ≈ 1.2% propagates into the standard deviation of the three-dimensional MAR, ϵ_MAR3D ≈ 0.95%; this spread is well within the standard deviation ϵ_MARcaustic ≈ 1.9% of the MAR distribution estimated from the caustic mass profiles. For the high-mass bin, where the number of clusters decreases from 2000 to 50, the velocity standard deviation becomes ϵ_vi ≈ 4.0%, which implies an ϵ_MAR3D ≈ 6.6% standard deviation of the real three-dimensional MAR; this ϵ_MAR3D is smaller than the ϵ_MARcaustic ≈ 10.9% standard deviation of the MAR estimated from the caustic mass profiles. We thus conclude that assuming the same mean radial velocity profile for every cluster in a given mass and redshift bin does not introduce any systematic bias in the MAR estimate. The uncertainty in the MAR is actually dominated by the uncertainty in the mass profile.

To estimate the MAR of galaxy clusters, we also need the accretion time t_inf (see Eq. (2)). Following De Boni et al. (2016), we adopt t_inf = 10⁹ yr, which is comparable with the dynamical time t_dyn = R/σ for the clusters we consider here, where R and σ are their size and one-dimensional velocity dispersion, respectively. To see that this value of t_inf is a reasonable choice, consider a homogeneous spherical system which is α time denser than the critical density and thus has mass . When the system is in virial equilibrium, its potential energy W = −3GM²/5R and kinetic energy K = 3Mσ²/2 can be combined in the virial relation K = |W|/2, which yields GM = 5Rσ²/4. Therefore, yr when α ∼ 200. The size and the velocity dispersion of the cluster progenitors at earlier times decrease by comparable factors; therefore, t_dyn remains roughly constant thus justifying the choice t_inf = 10⁹ yr for any cluster at any redshift.

In principle, we could adopt a different value for t_inf. However, the thickness of the shell derived to estimate the MAR is proportional to t_inf, and adopting t_inf larger than 10⁹ yr would generally return a shell thicker than ∼0.5R₂₀₀, which is the typical value we estimate in our analysis below. Thicker shells would thus include regions beyond ∼2.5−3R₂₀₀. In these regions, real clusters suffer from spectroscopic incompleteness (see Sect. 4.3) and t_inf > 10⁹ yr would thus return a biased estimate of the MAR. Adopting t_inf smaller than 10⁹ yr would instead return thinner shells; in this case, the shell would contain a handful number of galaxies in the redshift diagram and the MAR estimate would be affected by relevant shot noise.

3.4. Estimates of the MAR of simulated clusters

We now apply the caustic technique to each of our mock catalogues to derive the caustic mass profile and to estimate the MAR with our recipe (Eq. (2), Sect. 2). We bin the mock catalogues according to the cluster redshift and mass.

As we show with real clusters (Sect. 5), some of the individual R − v_los diagrams do not support an estimate of the MAR. Our recipe requires that we estimate the mass of the shell of radii R_i = 2R₂₀₀ and (1 + δ_s)R_i ∼ 2.5R₂₀₀. The caustic method estimates this mass from

(4)

where ℱ_β is the filling factor (Diaferio 1999; Serra et al. 2011) and 𝒜(R) is the amplitude of the caustics, namely, the vertical separation of the upper and lower caustics at radius R. At such large radii, a system can return a caustic amplitude 𝒜(R) = 0, either because of poor sampling (especially in real systems) or because of galaxy-rich background or foreground structures that inhibit the caustic technique from properly identifying the caustic location. In these cases, the mass of the shell, and thus the MAR, cannot be estimated.

In the samples of real clusters, we visually inspect the R − v_los diagrams to identify systems where the caustic technique fails. With the mock catalogues, we adopt an automatic procedure: to be conservative we remove the R − v_los diagrams where, in the range [R_i; (1 + δ_s)R_i], the caustic amplitude 𝒜(R) < 200 km s⁻¹ and 𝒜(R) < 400 km s⁻¹, for the low- and high-mass bin, respectively. The caustic technique algorithm also prohibits unphysical increases of 𝒜(R) with increasing R. This feature of the algorithm is more effective in real clusters than in mock clusters; in fact, the simulated clusters tend to have less sharp separation between the cluster members and the foreground and background galaxies (Diaferio 1999) and the caustic amplitude 𝒜(R) can remain artificially large at large R. Therefore, we also remove the R − v_los diagrams where the caustic amplitude 𝒜(R) > 2000 km s⁻¹ in the range [R_i; (1 + δ_s)R_i]. This procedure removes 31% and 16% of the R − v_los diagrams for the low- and high-mass bin, respectively. In Fig. 5, the red triangles show the median of the caustic MARs of the clusters at each redshift bin. The upper and lower sets of points are relative to the two cluster mass bins, as indicated in the figure. The error bars also show the 68th percentile range of the distribution of the estimated caustic MARs.

To quantify the systematic errors introduced by the projection effects in realistic observations, the blue squares in Fig. 5 show the MAR estimated with the recipe of Sect. 2 but with the correct three-dimensional mass profile. The procedure applied to obtain these estimates coincides with the procedure described in De Boni et al. (2016). The two estimates of the average MAR agree at all redshifts and for both mass bins: the average relative difference is ≲19%.

The difference between the three-dimensional MAR and the caustic MAR appears in their relative spreads. The relative spreads σ_MAR/MAR of the MAR obtained with the three-dimensional mass profiles are ∼42% and ∼46% for the low- and the high-mass bin, respectively. In contrast, using the mass profiles estimated with the caustic method, the relative spreads are ∼86% and ∼75%, respectively. These spreads are consistent with the spread of the caustic mass profile around the true mass profile estimated in N-body simulations (Serra et al. 2011, Fig. 12). According to Serra et al. (2011), this spread mainly originates from projection effects. Consequently, projection effects are also mainly responsible for the spreads of the MAR estimated with the caustic method shown in Fig. 5.

In addition, the larger spread in the low-mass bin originates from the fact that in N-body simulations, less massive systems have less well defined structure in redshift space compared to more massive systems (Diaferio 1999). Consequently, the identification of the caustics is prone to larger random errors. The algorithm takes this effect into account by associating a larger uncertainty with the caustic location and with the mass profile. Redshift space structures are sharper in the real universe (Schmalzing & Diaferio 2000; Casagrande & Diaferio 2006) and applications of the caustic technique to numerous real clusters have indeed shown that the errors in the caustic mass profiles estimated with N-body simulations are probably upper limits (Geller et al. 1999; Diaferio et al. 2005; Geller et al. 2013).

Figure 5 demonstrates that the caustic technique should return an unbiased estimate of the median three-dimensional MAR of real clusters at any redshift. In addition, the uncertainties overestimate the spread in the MAR based on three-dimensional mass profiles of a sample of clusters of comparable mass. We conclude that applying our recipe to real clusters should return a robust estimate of their MAR.

4. Catalogues of real clusters

Estimating the MAR of real clusters requires a dense redshift survey of the cluster outer regions. The largest catalogues currently available that satisfy this condition are the Cluster Infall Regions in the Sloan Digital Sky Survey (CIRS; Rines & Diaferio 2006) and the Hectospec Cluster Survey (HeCS; Rines et al. 2013). The former catalogue contains clusters at z < 0.1 and the latter contains clusters in the redshift range 0.1 < z < 0.3. These redshift ranges enable us to measure the MAR as a function of redshift. In Sects. 4.1 and 4.2, we review the main features of the two catalogues; Sect. 4.3 discusses some systematic effects in the selection of the galaxy samples and Sect. 5 describes the estimated MARs.

4.1. CIRS

The CIRS project extended the analysis of the CAIRNS survey (Rines et al. 2003), which pioneered the study of the infall region of clusters; CAIRNS used nine nearby galaxy clusters observed by the 2MASS survey (Jarrett 2004), exploiting extensive spectroscopy and near-infrared photometry.

CIRS is based on the Fourth Data Release (DR4) of SDSS, a photometric and spectroscopic wide-area survey at high galactic latitudes and low redshifts (Stoughton et al. 2002). The DR4 includes 6670 deg² of imaging data and 4783 deg² of spectroscopic data (Adelman-McCarthy et al. 2006). It was thus possible to extend the study of infall patterns around clusters initiated by CAIRNS to a larger number of clusters. By matching four X-ray cluster catalogues derived from the ROSAT All Sky Survey (RASS; Voges et al. 1999) to the spectroscopic area covered by the DR4, Rines & Diaferio (2006) obtained the CIRS catalogue, a sample of 74 clusters at z < 0.1, perfectly suited to the study of infall regions with the caustic technique.

In our analysis, we used the updated catalogues of the CIRS dataset obtained by compiling the SDSS DR14 spectroscopic sample; these new catalogues are now part of the HeCS-omnibus survey (Sohn et al. 2020). Here, three clusters of the original sample are removed: NGC4636, NGC5846 and Virgo. These clusters have redshifts < 0.01 and are poorly sampled. The remaining 71 clusters constitute our CIRS sample. Figure 6 shows the selection function of the CIRS clusters: their X-ray flux limit in the 0.1–2.4 keV band is 3 × 10⁻¹² erg s⁻¹ cm⁻² and their redshift range is [0, 0.1].

Fig. 6. Luminosities of the CIRS (blue points) and HeCS (orange points) clusters as a function of redshift. The superimposed curves show the flux and redshift selection functions of the two catalogues. Filled points refer to the clusters where we can compute an individual MAR.

Table 2 lists the celestial coordinates, redshift, size R₂₀₀, and the corresponding mass M₂₀₀. It also lists the individual MARs estimated with the caustic method and discussed in Sect. 5.1⁴. The blue histogram in Fig. 7 shows the mass distribution of the CIRS clusters. Table 3 lists the medians and percentile ranges of the redshift and mass distributions for both the complete catalogue and the subset of the CIRS clusters for which we estimate the individual MAR.

Fig. 7. Mass distribution of the CIRS (blue histogram) and HeCS (orange histogram) clusters. The black and red dashed lines show the median masses for the CIRS and HeCS catalogues, respectively. The area under each histogram is normalised to unity.

4.2. HeCS

HeCS is the first systematic and extensive spectroscopic survey of the infall regions of clusters at z ≥ 0.1 (Rines et al. 2013). HeCS takes advantage of the SDSS and RASS surveys. In particular, existing X-ray cluster catalogues based on RASS were used to define flux-limited cluster samples that were then matched to the imaging footprint of the SDSS DR6 (Adelman-McCarthy et al. 2008). Figure 6 shows the selection function of the HeCS clusters: their X-ray flux limit in the 0.1–2.4 keV band is 5 × 10⁻¹² erg s⁻¹ cm⁻² and their redshift range is [0.1, 0.3]. We include four additional clusters below the flux limit but with fluxes > 3 × 10⁻¹² erg s⁻¹ cm⁻².

The imaging footprint of the SDSS DR6 includes 8417 deg² of imaging data. The multicolour photometry enabled the selection of candidate cluster members using the red-sequence technique. At z ≥ 0.1, the SDSS spectroscopic survey is not dense enough for accurate measurements of the cluster masses with the caustic technique; therefore, the MMT/Hectospec instrument (Fabricant et al. 2005) was used to obtain spectroscopic data for the candidate members. Recently, Sohn et al. (2020) updated the HeCS dataset by using the spectroscopic data of SDSS DR14 and incorporated the new catalogues in the HeCS-omnibus survey.

The HeCS survey contains 58 clusters in the redshift range 0.1 < z < 0.3, for a total amount of 22 680 observed galaxy redshifts, of which 10,145 are cluster members. Each cluster survey typically includes ∼400−550 redshifts; in general, roughly half of these galaxies are cluster members and the remaining galaxies are foreground or background objects.

The orange histogram of Fig. 7 shows the mass distribution of the HeCS cluster sample. This sample includes fewer low-mass clusters than CIRS because it covers a deeper redshift range [0.1, 0.3]. The HeCS sample also contains more high-mass clusters as a result of the larger survey volume. Due to the extended mass range of the HeCS cluster sample (see Fig. 7), we separate the 58 clusters sorted by mass into two subsamples of 29 clusters each. The low-mass and high-mass samples have a median mass of M₂₀₀ = 1.86 × 10¹⁴ h⁻¹ M_⊙ and M₂₀₀ = 5.61 × 10¹⁴ h⁻¹ M_⊙, respectively.

Table 4 lists the HeCS clusters separated into the two subsamples⁵. Table 3 lists the medians and the percentile ranges of the redshift and mass distributions of the entire HeCS catalogue and of its subsamples, including the subsets of clusters for which we estimate the individual MARs in Sect. 5.1.

4.3. Effects of the selection of the galaxy samples

4.3.1. Photometric completeness

The galaxies in the HeCS clusters are selected according to their red colours, whereas the galaxies in CIRS are from a magnitude-limited survey. Therefore, in principle, unlike CIRS, a substantial number of blue galaxies could be missing in HeCS.

To quantify the impact of these different selections, Fig. 8 shows the colour-magnitude diagram of the two catalogues. We only show the galaxies at projected distance smaller than 3R₂₀₀ from the cluster centre and with a line-of-sight velocity |v_los − v_cl|< 3000 km s⁻¹, where v_cl is the line-of-sight velocity of the cluster. In the right panel, the black solid line shows the fit of Rines et al. (2013) to the red sequence for all the member galaxies in HeCS; in the left panel, the solid line shows the red sequence fit of Rines et al. (2013) shifted by -0.03 mag. The black dashed, dotted, and dash-dotted lines show the boundaries of the stripes used for sample selection obtained by changing the intercept of the fit by ±0.30, ±0.40, ±0.50 mag.

Fig. 8. (g − r) colour-magnitude diagrams of the CIRS (left) and HeCS (right) galaxies, including k-corrections. The galaxies have a cluster-centric distance smaller than 3R200 and line-of-sight velocity |vlos − vcl|< 3000 km s−1. The black solid curve in the right panel shows the fit by Rines et al. (2013) derived from the HeCS galaxies. The black dashed, dotted, and dash-dotted lines show the ±0.30, ±0.40, ±0.50 shifts of the same fit. In the left panel, the black solid line shows the Rines et al.’s fit with an offset of −0.03 mag. The black dashed, dotted, and dash-dotted lines show the ±0.30, ±0.40, ±0.50 mag shifts of this line.

Some of the galaxies outside the selection stripes are extremely red. However, they are a tiny minority and they are almost never cluster members (Rines et al. 2013). The impact of red objects can thus be ignored here. The spatial distribution of the red galaxies peaks within R₂₀₀; the blue galaxy distribution peaks at significantly larger radius. The distributions of the line-of-sight velocities of both red and blue galaxies are centred on zero and are approximately Gaussian; the width of these distributions is smaller for the red galaxies than for the blue galaxies, which is in agreement with previous observations of galaxies in clusters (e.g. Dickens & Moss 1976; Carlberg et al. 1997; Geller et al. 2014; Barsanti et al. 2016) and with the expectations from simulations of galaxy formation (e.g. Diaferio et al. 2001).

Table 5 lists the ratio Ξ between the number of the blue galaxies outside each of the stripes shown in Fig. 8 and the total number of galaxies in CIRS and HeCS. The fraction of blue galaxies is ≲10% in both catalogues, for any stripe. The relative difference in blue galaxies between CIRS and HeCS, δ_Ξ = 1 − Ξ_HeCS/Ξ_CIRS, decreases from 29% to 2% from the ±0.30 mag-stripe to the ±0.50 mag-stripe.

These numbers suggest that HeCS might roughly miss, at most, one third of blue galaxies compared to CIRS. The effect on the caustic location should thus be mild. In addition, the caustic technique locates the caustics by adopting a threshold of the number density distribution of the galaxies in the R − v_los diagram that is set by the galaxies within ∼R₂₀₀. The spatial distribution of the red galaxies peaks within this radius, whereas the blue galaxy distribution peaks further out, as mentioned above.

We prove that missing this fraction of blue galaxies in HeCS does not affect our estimates of the MAR by taking ten CIRS clusters and randomly removing 35% of the galaxies outside the ±0.30 mag-stripe, a slightly larger fraction of the upper limit of δ_Ξ listed in Table 5. The differences between these caustic mass profiles and the original profiles are within the uncertainties of the caustic technique and they are thus statistically indistinguishable.

Our conclusion is further supported by Rines et al. (2013). They present a case study of three HeCS clusters: with additional spectroscopic observations including blue galaxies, they can quantify how these galaxies can affect the estimate of the velocity dispersion and the dynamical mass. They find that a small fraction of blue galaxies are actual cluster members and their inclusion only increases the velocity dispersion by 0.3%. Rines et al. (2013) conclude that targeting the galaxies on the red-sequence alone does not produce any bias on the velocity dispersion and the mass estimates.

4.3.2. Spectroscopic completeness

In general, the spectroscopic completeness, namely, the ratio between the number of galaxies with spectroscopic measurements and the number of galaxies with photometric measurements, can decrease with increasing distance from the cluster centre. In CIRS, the spectroscopic completeness could be affected by the edges of the footprint of SDSS, whereas in HeCS, spectroscopic redshifts are only available within 30 arcmin from the cluster centre. Hence, the spectroscopic measurements can be incomplete at large radii. For the caustic method, this incompleteness can cause an underestimate of the caustic amplitude and thus an underestimation of the cluster mass and of the MAR. However, when the sample is too sparse, the caustic method refrains from locating the caustics and returns no mass estimates. In contrast, if the incompleteness is not too severe to prevent the location of the caustics, the underestimate of the mass might be within the caustic mass uncertainty. Towards the end of this section, we confirm that this is indeed the case for the CIRS and HeCS clusters.

We are first concerned with the fact that the spectroscopic incompleteness could affect the two catalogues differently and, thus, bias our MAR estimates in different ways. We show that this case does not occur and that the spectroscopic incompleteness is comparable in the two CIRS and HeCS samples.

The left panel of Fig. 9 shows the ratio 𝒫 = 𝒫_]2, 2.5]/𝒫_]0, 1] between the numbers of galaxies with photometric measurements in the range of r/R₂₀₀ ∈  ]2, 2.5] and in the range r/R₂₀₀ ∈  ]0, 1]. The right panel shows the ratio 𝒮 = 𝒮_]2, 2.5]/𝒮_]0, 1] of the numbers of galaxies with spectroscopic redshift in the same two regions. Blue and orange points show the CIRS and HeCS clusters, respectively. The ratios, 𝒫 and 𝒮, of each cluster are computed with the galaxies brighter than M_r = −20. This limit is , where is the characteristic red-band magnitude of the Schechter luminosity function, as indicated by photometric studies of A2029 and Coma (Sohn et al. 2017). The left panel shows that, on average, the photometric samplings of the CIRS and HeCS surveys in the central and outer regions of the clusters are comparable: ⟨𝒫⟩ = 0.99 ± 0.61 and ⟨𝒫⟩ = 0.97 ± 0.27 for CIRS and HeCS, respectively. A similar result holds for the spectroscopic samplings: ⟨𝒮⟩ = 0.65 ± 0.30 and ⟨𝒮⟩ = 0.73 ± 0.22 for CIRS and HeCS, respectively.

Fig. 9. Left panel: ratio of the numbers of galaxies with photometric data in the range r/R200 ∈ ]2, 2.5] and within r/R200 = 1 within each cluster, against the cluster redshift. Right panel: ratio between the numbers of galaxies with spectroscopic redshifts within the same regions. Blue and orange points show the CIRS and HeCS clusters, respectively. Solid points show the clusters for which we estimate the individual MAR.

The values of the ratio 𝒮 < 1 suggest that the spectroscopic completeness might decrease with distance from the cluster centre. In fact, Fig. 10 shows the spectroscopic completenesses 𝒮_]0, 1]/𝒫_]0, 1] and 𝒮_]2, 2.5]/𝒫_]2, 2.5] in the central and outer regions of the clusters. It also confirms a spectroscopic incompleteness at large radii. Nevertheless, the spectroscopic completenesses in the central and outer regions are comparable in the two catalogues. In the centre, CIRS and HeCS have mean spectroscopic completeness 0.75 ± 0.12 and 0.65 ± 0.15, respectively. In the outer region, the mean completenesses are 0.53 ± 0.19 and 0.49 ± 0.15. We obtain similar results by considering the subsample of clusters for which, as we illustrate in the next section, we estimate the MAR individually. The similarity of the ratios 𝒫 and 𝒮 of CIRS and HeCS in Fig. 9 shows that the comparable completenesses are not a fluke originating from different photometric samplings in the central and outer regions of the two catalogues.

Fig. 10. Spectroscopic completeness of the CIRS (blue points) and HeCS (orange points) clusters, against their redshifts, within r/R200 = 1 from the cluster centre (upper panel) and within the range r/R200 ∈ ]2, 2.5] (lower panel). Solid points show the clusters for which we estimate the individual MAR.

Based on these results, we conclude that the spectroscopic incompleteness as a function of radius is present in our samples, but it might similarly bias our MAR estimates of the CIRS and HeCS clusters. This conclusion is further supported by the fact that the fraction of clusters with null amplitude beyond 2R₂₀₀ is ∼50% in both catalogues, as we show in the next section. In turn, these fractions are larger than the fractions ∼13−18% found for the simulated clusters (Sect. 3.4).

We now quantify the possible systematic error resulting from the spectroscopic incompleteness as a function of radius. We create mock catalogues by stacking the simulated clusters at z = 0.12, the intermediate redshift between the average redshifts of the CIRS and HeCS clusters. We stack all the 141 clusters of the high-mass bin, whereas we stack a random sample of only 35 clusters of the low-mass bins, to limit the computational time. The R − v_los diagrams of these two stacked clusters contain ∼50,000 particles, and simulate catalogues that are spectroscopically complete. We also create mock catalogues by randomly removing some particles within r/R₂₀₀ = 1 and within the range of r/R₂₀₀ ∈  ]2, 2.5] to simulate the spectroscopic incompleteness as a function of radius. We consider two undersampled catalogues: in the first catalogue, we remove 25% and 47% of the particles in the central and outer regions, respectively; in the second catalogue, we remove 35% and 51% of the particles in the two regions. These fractions are chosen accordingly to the mean spectroscopic incompleteness estimated above for the CIRS and HeCS clusters.

The undersampled mock catalogues return mass and MAR estimates consistent with the estimates obtained from the complete catalogues: with the undersampled catalogues, the mass M₂₀₀ is underestimated by ∼13%, whereas the MARs is overestimated by ∼22%, on average. Both values are within the uncertainties of the estimates obtained from the complete catalogues and are thus statistically indistinguishable. The MAR is, on average, overestimated rather than underestimated, as one might have naively expected; this overestimation confirms that the spectroscopic incompleteness as a function of radius in the CIRS and HeCS clusters can generate statistical fluctuations on the MAR estimates, rather than a systematic error.

5. Measure of the mass accretion rate of real clusters

Here, we discuss our estimates of the MAR of individual clusters (Sect. 5.1) and the average MAR of the cluster samples (Sect. 5.2) as a function of mass and redshift.

5.1. Individual MARs

As anticipated in Sects. 4.1 and 4.2, we only estimate the individual MAR of a subset of clusters, selected by visually inspecting their R − v_los diagrams. In fact, we remove (i) the clusters whose caustic amplitude shrinks to zero within the infalling mass shell and (ii) the clusters whose caustics have unphysical spikes. These cases usually occur in the presence of galaxy-rich background or foreground structures that prevent the caustic technique algorithm from properly identifying the caustic location. This procedure is similar to the automatic procedure adopted for the mock catalogues in Sect. 3.4. Visual inspection selects 30 clusters (out of 71) for the CIRS sample. We selected 18 out of 29 and 16 out of 29 clusters, respectively, for the low- and high-mass HeCS samples.

The initial infall velocity v_i entering Eq. (1) depends on the cluster mass and redshift. Rather than consider the mass and redshift of each cluster, we consider the same v_i for all the clusters within each sample. We thus estimate the value of v_i appropriate for the median redshift and the median mass of each cluster sample.

The median mass and redshift of the CIRS subsample are M₂₀₀ = 1.1 × 10¹⁴ h⁻¹ M_⊙ and z = 0.064 (Table 3). The closest redshifts of the simulated snapshot are z = 0 and z = 0.12. We thus estimate v_i appropriate for the CIRS median mass with three linear interpolations on the simulation information. For each of the two simulated samples of mass 10¹⁴ and 10¹⁵ h⁻¹ M_⊙, we first interpolate between the two median masses at redshifts z = 0 and z = 0.12 listed in Table 1 to estimate the appropriate median masses at z = 0.064, and . The second interpolation returns the velocities appropriate for z = 0.064, and , for each of the two simulated samples: for each sample, we consider the radial velocity profiles at the two redshifts z = 0 and z = 0.12 and consider the value of the velocity at the single radius lying in the range [2−2.5]R₂₀₀; and are derived from the interpolation between the two values at redshifts z = 0 and z = 0.12 for each sample. Finally, to obtain v_i appropriate for the median mass M₂₀₀ = 1.1 × 10¹⁴ h⁻¹ M_⊙ at z = 0.064, we interpolate between the two median masses, and , and the two velocities, and . We find v_i = −170 ± 3 km s⁻¹.

For the sake of completeness, the value of v_i above also includes the uncertainty derived with an analogous interpolation based on the profile of the standard deviation of the mean radial velocity profile; the standard deviation is comparable to the 68th percentile range shown by the shaded bands in Fig. 4. As discussed in Sect. 3.3, the spread on v_i is roughly a factor three smaller than the uncertainty on the MAR derived from the caustic technique and does not generate any bias on the MAR itself. Therefore, we ignore the uncertainty on v_i when computing the shell thickness (Eq. (1)).

For the two HeCS subsamples, we adopt the same procedure. The median mass and redshift of the low-mass HeCS subsample are M₂₀₀ = 2.2 × 10¹⁴ h⁻¹ M_⊙ and z = 0.14 (Table 3), and we use the snapshots of the simulation at redshifts z = 0.12 and z = 0.19. Finally, we find that v_i = −288 ± 8 km s⁻¹.

For the high-mass HeCS subsample, the median mass is M₂₀₀ = 5.2 × 10¹⁴ h⁻¹ M_⊙ and the median redshift is z = 0.21 (Table 3). The closest redshifts of the simulated snapshots are z = 0.19 and z = 0.26. We thus find v_i = −566 ± 19 km s⁻¹.

The uncertainty in the MAR only depends on the uncertainty, σ_Mshell, on the mass M_shell of the infalling shell because we adopt a value for t_inf without an uncertainty. The caustic method estimates the mass profile from the caustic amplitudes, 𝒜_i, measured at a set of radii, r_i. The mass, M_i, of the ith shell is proportional to . According to Diaferio (1999) and Serra et al. (2011), we thus estimate the uncertainty on the mass of the infalling shell as:

(5)

where the sum extends over the shells of the caustic mass profile within the infalling shell, and σ_𝒜i is the uncertainty on the caustic amplitude of the ith shell; this uncertainty increases with decreasing ratio between the number of galaxies within the caustics and the number of galaxies within the R − v_los diagram at each radius r_i. Below, we find that the relative uncertainties of the individual MARs are 17% on average; for only 7 out of 64 clusters we measure a MAR with a relative uncertainty larger than 35%. The uncertainties mostly come from (1) the assumption of spherical symmetry and (2) the presence of very populated R − v_los diagrams that make it challenging to locate the caustics.

Tables 2 and 4 list the individual MARs and Fig. 11 shows the distributions of the individual MARs of our three cluster subsamples. The density distributions are normalised to unity for an easier comparison. The blue, orange, and green histograms refer to the CIRS, low-mass HeCS, and high-mass HeCS clusters, respectively; the same colours are used for the vertical dashed lines showing their median MAR. The distributions generally show a peak at small MARs and an extended tail of large MARs. In general, the small and large MARs are associated with low- and high-mass clusters, respectively.

Fig. 11. Distributions of the individual MARs. The blue, orange and green histograms refer to the CIRS, low-mass HeCS, and high-mass HeCS clusters, respectively. The dashed lines with the same colours show the median MAR of each sample. The area under each histogram is normalised to unity.

Figures 12 and 13 show R − v_los diagrams and cumulative mass profiles for six sample clusters. In all the panels, the black dashed line shows the inner radius of the infalling shell; the red dashed line shows the outer radius of the shell (Eq. (1)). In the two figures, the mass of the cluster, M₂₀₀, increases from top to bottom. More massive systems tend to be surrounded by larger amounts of mass and the caustic amplitude decreases more slowly with radius. Consequently, we expect that the mass of the infalling shell, and thus the MAR, is correlated with M₂₀₀.

Fig. 12. Three examples of the estimation of the MAR for individual clusters in CIRS. The left column shows the R − vlos diagram for each cluster. The right column shows the corresponding caustic mass profile (solid line) and the NFW fit (dot-dashed line). The black and red dashed lines show the inner and outer radius of the infalling spherical shell. The M200 of the clusters increases from top to bottom.

Fig. 13. Same as Fig. 12 for three HeCS clusters. Zw3179 and A655 belong to the low-mass sample, whereas A963 is in the high-mass sample.

Figure 14 shows the correlation between M₂₀₀ and the MAR. There is also, as expected, a correlation with redshift. The positive correlation with redshift and M₂₀₀ is clearly expected in the hierarchical clustering scenario where more massive halos lie in higher density regions and they are surrounded by larger amounts of mass (Bardeen et al. 1985; Lacey & Cole 1993; McBride et al. 2009; Fakhouri et al. 2010; van den Bosch et al. 2014). Similarly, halos with comparable masses are expected to have larger MARs at larger redshifts. The correlations are statistically significant according to Kendall’s test: the coefficients, τ, are 0.516 and 0.528 for the MAR versus M₂₀₀ and versus z, respectively, with corresponding significance levels of p_M200 = 1.7 ⋅ 10⁻⁹ and p_z = 7.6 ⋅ 10⁻¹⁰.

Fig. 14. MAR of individual clusters as a function of their mass M200; the colour code shows the dependence on redshift. Filled circles, squares and triangles refer to the CIRS, low- and high-mass HeCS sample, respectively. The open circles show the median MARs of the N-body clusters of our two simulated samples; for these clusters, we estimate the MAR from their three-dimensional mass profiles.

To compare these measurements with our simulated clusters, we consider each simulated cluster sample at the four redshifts z = 0.0, 0.12, 0.19, and 0.26, and separate each sample into four mass bins; these splittings yield 16 subsamples for each of the 10¹⁴ and 10¹⁵ h⁻¹ M_⊙ sample of simulated clusters. The open circles in Fig. 14 show the median MARs of these 32 subsamples. The ΛCDM expectations appear fully consistent with our measurements.

We close this section with a brief comment about two specific clusters, namely: A750 and A1758. In agreement with Rines et al. (2013), we find that the galaxies in the redshift catalogue of A750 belong to two different clusters at different redshifts: A750 at z = 0.164, and MS0906 at z = 0.177. The caustic technique successfully identifies the two clusters and estimates the two MARs separately. Their individual MARs are listed in Table 4. A750 belongs to the low-mass HeCS sample, whereas MS0906 is in the high-mass sample.

Similarly, a weak lensing analysis, based on B and V passbands images of A1758 at redshift z = 0.28, by Ragozzine et al. (2011), suggests that this cluster is a system of four gravitationally bound substructures currently undergoing two separate mergers. The mass estimate derived by the caustic technique is not affected by the presence of substructures (Diaferio 1999) and we can thus estimate the MAR of A1758 without any particular precaution (Table 4).

5.2. Average MAR

Figure 15 shows the median of the individual MARs as a function of redshift. The green filled circle at redshift z = 0.064 shows the median MAR of the CIRS clusters; the green filled circles at redshift z = 0.14 and z = 0.21 show the median MARs of the two HeCS subsamples. We plot each green circle at the median redshift of the subsample. The number close to each of them shows the median mass of each sample in units of h⁻¹ M_⊙. The error bars show the 68th percentile ranges of the distributions of the MAR and of the redshift of the clusters.

Fig. 15. Median MARs of real (green) and simulated (blue and red) clusters. The green error bars show the 68th percentile ranges of the MAR and redshift distributions of the real clusters. The value close to each green filled circle is the median mass M200 of the real clusters in units of h−1 M⊙. The blue squares and red triangles are for the MAR based on the three-dimensional and the caustic mass profiles of simulated clusters, respectively (from Fig. 5). The light-blue shaded areas and the red error bars are the 68th percentile range of the distribution of the individual MAR of the simulated clusters. The open circles show the median MARs of real clusters for different values of the infall velocity. We show the results where we changed the initial infall velocity that we adopted by: −40% (blue), −20% (cyan), +20% (orange), and then +40% (red). The black circles show the MAR estimated with an initial infall velocity of vi = 0.

To compare these measured MARs with the expectations of the ΛCDM model, Fig. 15 shows the results of Fig. 5 for the simulated clusters extracted from the L-CoDECS simulation (Sect. 3.2). The blue squares show the medians of the MARs derived from the three-dimensional mass profiles. The red triangles show the medians of the MARs derived from the caustic mass profiles estimated from the mock redshift surveys.

This figure confirms the result of Fig. 14: the medians of the MARs of real clusters fall within the range of the MAR of simulated clusters. The three median masses of the real cluster samples, 1.1, 2.2, and 5.2 × 10¹⁴ h⁻¹ M_⊙, are in between the two median masses, 10¹⁴ and 10¹⁵ h⁻¹ M_⊙, of the two samples of the simulated clusters.

Figure 15 also shows that the spreads in the median MARs of the real clusters are comparable with or even smaller than the spreads of the mock catalogues. This result supports the conclusion that the caustic technique returns robust estimates of the average MAR of clusters and that mock catalogues tend to overestimate the expected uncertainties because the caustics are usually less well-defined in the simulations than in the data (Sect. 3.4).

We now quantify the effect of the value of the initial velocity v_i on the estimate of the MAR. The open circles in Fig. 15 show the median MARs of each cluster sample when we decrease v_i by 20 or 40% (cyan and blue circles) or increase v_i by 20 or 40% (orange and red circles) with respect to the adopted v_i (green dots with error bars). The estimated MAR does indeed depend on v_i, but the resulting MARs remain within the 68th percentile range of the MAR distribution. The extreme and unrealistic choice of v_i = 0 makes the estimated MAR (black circles) decrease substantially: the MAR then disagrees significantly with the ΛCDM model expectations. Nevertheless, the correlations between the MAR, redshift, and cluster mass persist regardless of the choice of v_i.

A simple analysis explains qualitatively why the correlations between MAR, redshift, and cluster mass are so robust. Around each cluster, we consider the mass M_2−2.5 of the shell with inner and outer radius 2R₂₀₀ and 2.5R₂₀₀, respectively. This shell is comparable to the shell usually identified by the solution of Eq. (1). Figure 16 shows M_2−2.5 as a function of the cluster mass M₂₀₀. The colour and symbols are the same as in Fig. 14. The two figures are qualitatively similar and show, as expected, that the correlations we see for the MAR derives directly from the correlations between M_2−2.5, M₂₀₀, and z. We also show the shell mass M_2−2.5 derived from the NFW fits to the three-dimensional mass profiles of our simulated samples. On average, for the clusters in the low-mass bin, the shell mass derived from the NFW fits is ∼22% larger than the actual mass. This average overestimate increases to ∼28% for the clusters in the high-mass bin. Because this shell mass strongly correlates with the MAR, we expect a comparable overestimation of the MAR if we adopt the NFW fits rather than the mass profiles estimated with the caustic method. In Sect. 3.2, we show that, over the entire radial range [0−4]R₂₀₀, the NFW mass profile overestimates the true mass profile by ∼5%, on average. The mass of the shell estimated with the NFW mass profile overestimates the true shell mass by a larger factor, ∼22−28%, because the shell mass derives from the difference of the masses estimated at two different radii, 2R₂₀₀ and ∼2.5R₂₀₀.

Fig. 16. Mass M2−2.5 of the shell with radii 2R200 and 2.5R200 as a function of M200. The colour code shows the dependence on redshift. Circles, squares and triangles refer to the CIRS, low- and high-mass HeCS sample, respectively. The open circles show the median M2−2.5 of the N-body clusters of our two simulated samples, estimated from their three-dimensional mass profiles. The crosses show the median M2−2.5 obtained from the NFW fits to the three-dimensional mass profiles: on average, they overestimate the individual true mass of the shell by ∼22−28%, depending on the mass bin.

A different strategy for estimating the mean MAR of real clusters is to apply our recipe to an average cluster obtained by stacking the real clusters (Rines & Diaferio 2006; Serra et al. 2011). This approach has the advantage of considering all the clusters of each sample, including those clusters where the individual MAR cannot be estimated. In addition, stacking the clusters averages out the deviations from spherical symmetry of the individual clusters.

Similarly to HeCS, we separated the CIRS sample into a low-mass and a high-mass subsample, according to the cluster median mass M₂₀₀ = 1.78 × 10¹⁴ h⁻¹ M_⊙. Table 6 lists the median masses and the 68th percentile ranges of the mass and redshift distributions of the four subsamples. For each cluster subsample, we built a R − v_los diagram containing all the galaxies in the field of view with line-of-sight velocity within 4000 km s⁻¹ from the cluster velocity and within 10 h⁻¹ Mpc from the cluster centre. The ‘Total’ column of Table 7 shows the number of galaxies within these four stacked clusters.

In the R − v_los diagram of the stacked cluster, the velocities and the radial distances are expressed in the reference frame of each cluster without any additional normalisation. This simple approach is justified by the homogeneity of mass and size for the clusters within each subsample. This stacking procedure avoids introducing errors from the estimate of the size and velocity dispersion of the clusters. If we stacked clusters of widely different mass we would require an additional normalisation. The R − v_los diagrams of the stacked clusters are shown in the first and third panels of Figs. 17 and 18. The mass profiles estimated with the caustic technique applied to these R − v_los diagrams are shown in the second and fourth panels of Figs. 17 and 18. The dot-dashed lines in these panels show the NFW fits to these profiles performed over the radial range within 4R₂₀₀ or within the radius where the caustic amplitude shrinks to zero, if this radius is smaller than 4R₂₀₀. The masses M₂₀₀ of the stacked clusters are listed in Table 6.

Fig. 17. R − vlos diagrams and corresponding mass profiles of the two stacked clusters from the CIRS catalogue. The upper (lower) panels refer to the low-mass (high-mass) subsample. The black and red dashed vertical lines show the inner and outer radius of the spherical shell used to estimate the MAR. The shaded areas show the 50% confidence level of the caustic location and of the caustic mass profile according to the caustic technique recipe. The dot-dashed lines show the NFW fits to the caustic mass profiles; the green vertical dotted lines show the upper limit of the radial range over which the NFW fits have been performed.

Fig. 18. Same as Fig. 17 for the HeCS catalogue.

To estimate the MAR, we need to set the value of the initial infall velocity v_i. With the same procedure described in Sect. 5.1, we use the L-CoDECS simulations to find v_i = −164 ± 3 km s⁻¹, and v_i = −303 ± 12 km s⁻¹ for the low- and high-mass CIRS subsample, and v_i = −324 ± 10 km s⁻¹ and v_i = −735 ± 27 km s⁻¹, for the low- and high-mass HeCS subsample, respectively.

Figures 17 and 18 show the inner (black dashed line) and outer radius (red dashed line) of the infalling shell, according to Eq. (1). Averaging out the deviations from spherical symmetry and increasing the number of member galaxies make the caustic location and the mass profiles appear smoother than those of individual clusters (e.g. Figs. 12 and 13).

Figure 19 shows the MAR of the stacked clusters (orange circles), superimposed on the results obtained from the individual clusters (Fig. 15, green circles). The MARs of the stacked clusters are consistent with the median MARs of the individual clusters; they confirm the dependence of the MAR on cluster mass and redshift and agree with the predictions of the ΛCDM model.

Fig. 19. MARs of the stacked clusters of CIRS (the two orange circles at low redshift) and HeCS (the two orange circles at high redshift). The median MARs of the individual clusters from Fig. 15 are shown in green. The value close to each circle is the mass M200 in units of h−1 M⊙. The blue squares with the shaded bands are from Fig. 5 and show the median MARs and their spreads of the simulated clusters estimated with the three-dimensional mass profiles.

The uncertainty in the MARs of the stacked clusters are also consistent with the spread of the distributions of the individual MARs. The relative uncertainty in the MARs of the stacked clusters are ∼64%, an uncertainty that is larger than the relative uncertainties of ∼17% of the individual MARs. Although the stacked clusters more closely satisfy the assumption of spherical symmetry, the increase of the number of background and foreground galaxies in the R − v_los diagram makes the location of the caustics more uncertain and, consequently, increases the uncertainty on the mass profile⁶. Nevertheless, as expected, most members of the stacked clusters are members of the individual clusters: 75%, 81%, 80%, and 94% of the members of CIRS-STK1, CIRS-STK2, HeCS-STK1, and HeCS-STK2, respectively, are members of the individual clusters. Hence, the location of the caustic curves appears to be robust.

We can reduce the uncertainty on the mass profile by adopting a different strategy to build the stacked clusters. For each cluster, we then include only the member galaxies identified by the caustic technique, namely, the galaxies within the caustics in each R − v_los diagram. Again, we include all the clusters, irrespective of the fact that their caustic amplitude can be zero beyond 2R₂₀₀. The ‘Members-only’ column of Table 7 shows the number of galaxies within these four stacked clusters. Clearly, the number of galaxies in the R − v_los diagrams substantially drops compared to the previous stacked clusters: for the two mass bins of CIRS, the number of galaxies is reduced by 70% and 55%, respectively; for the HeCS subsamples, these numbers decrease by 55% and 43%. The R − v_los diagrams of these new stacked clusters have now well-defined trumpets and the caustics identify substantially fewer interlopers. As expected, the mean relative uncertainty in the estimated MAR drops to ∼36%, approximately half of the relative uncertainty of the previous stacked clusters. Although the new MAR estimates are ∼12% smaller, on average, than the previous MARs, because the R − v_los diagrams are now slightly less sampled than the original stacked clusters, the difference between the two MAR estimates is within the uncertainties and the estimates are thus statistically consistent with each other.

In the stacking procedure, for each cluster, we include all the galaxies in the field of view within 4000 km s⁻¹ along the line of sight and within 10 h⁻¹ Mpc from the cluster centre without any further normalisation of the line-of-sight velocity and projected distance of the cluster galaxies. The properties of the stacked cluster might thus mirror the properties of the clusters that contribute most galaxies and are presumably more massive. To quantify this effect, we then create new stacked clusters. For each cluster, we include the same number N of galaxies randomly sampled from the field of view: for each mass bin, we set N close to the number of galaxies of the cluster with the smallest number of galaxies. For the two mass bins of CIRS, we have N = 263 and N = 571; for HeCS, we have N = 402 and N = 408. These values of N preserve the optimal number of galaxies ≳200 within 3R₂₀₀ (Serra et al. 2011). Only 5 out of 36 clusters in the low-mass subsample of CIRS have fewer galaxies than N = 263; the clusters with fewer galaxies than N are 2 out of 35, 1 out of 26, and 2 out of 26 for the high-mass subsample of CIRS, and for the low- and high-mass subsamples of HeCS, respectively. The ‘Equally-weighted’ column of Table 7 shows the number of galaxies within these four stacked clusters. These equally weighted stacked clusters have a R − v_los diagram that is slightly less sampled than the original stacked clusters. As a consequence, their MARs are, on average, 9% lower than the previous MARs and their relative uncertainties are, on average, ∼25% larger. Nevertheless, their MAR estimates are still consistent with the results of Fig. 19.

6. Discussion

Our result is the first attempt to measure the MAR of clusters of galaxies by estimating the amount of mass in the outer regions. Our approach is crucially different from the approach usually adopted in the theoretical investigations of the accretion of dark matter halos. Most attention is focused on the connection between the accretion history of halos and their global properties at a particular epoch (e.g. Rey et al. 2019). These properties include the concentration parameter of the NFW profile (Navarro et al. 1997) and the amount of substructure (Jing 2000). More recently, Ludlow et al. (2013) show that the ubiquity of the NFW density profile in CDM models results from the mass independence of the accretion histories of dark matter halos.

The accretion histories of dark matter halos can, in principle, be a very sensitive probe of the nature of the dark matter particles and of the theory of gravity. For example, both self-interacting dark matter particles (Brinckmann et al. 2018) and warm dark matter particles (Ludlow et al. 2016) predict smaller concentrations for similar halo masses, which, in turn, affect the determination of the formation redshift of the dark matter halos. On the other hand, Oleśkiewicz et al. (2019) show that in the Hu & Sawicki (2007) form of f(R) gravity, the linear relation between concentration and formation time that is seen in CDM models becomes more complicated due to gravitational screening.

Unfortunately, when moving from dark matter halos to real cluster of galaxies, the differences detected in the theoretical relations might appear out of reach because the formation redshift is an unobservable quantity and the concentration parameter is prone to various systematic errors and sample biases (Mandelbaum et al. 2008; Bartelmann 2010; Oguri et al. 2012). Moreover, the properties of the central region of clusters, say at ≲0.1R₂₀₀, are complicated by the physics of baryons, as they can play a relevant role in the dynamics of the cluster centre (Blumenthal et al. 1986; Kravtsov et al. 2005).

The N-body simulations suggest that the mass that is falling onto the cluster for the first time is located beyond the splash-back radius ∼2R₂₀₀ (Diemer & Kravtsov 2014; More et al. 2015). However, to date most investigations of the accretion rate of real clusters focus on regions inside this radius (e.g. Lemze et al. 2013; Haines et al. 2018). Therefore, comparing their results with the MAR derived from N-body simulations is not straightforward.

On the contrary, the estimation of the mass at distances ≳2R₂₀₀ is a more straightforward approach because most of this mass is falling onto the cluster for the first time and is thus a direct estimate of the accretion rate. In addition, in this region, baryons are expected to have a limited effect because gravity is the main driver of the properties of the density and velocity fields. Estimating the mass in these outer regions is inaccessible to strong lensing measurements and to current X-ray observations, but this can be tackled with the caustic technique or weak gravitational lensing. Thanks to the caustic technique, estimating the mass at these large radii, where dynamical equilibrium does not hold, becomes feasible when appropriately dense spectroscopic surveys are available (Geller et al. 2011; Serra & Diaferio 2013; Sohn et al. 2018; Rines et al. 2018).

According to N-body simulations, at radii larger than ∼2R₂₀₀, not all the mass is falling onto the cluster for the first time (e.g. Aung et al. 2021), but ∼20−25% of this mass has already been within R₂₀₀ (Ludlow 2009; Diemer et al. 2017; Bakels et al. 2021; Xhakaj et al. 2020). When estimating the MAR of real clusters from the amount of mass in this outer region, this evidence can, in principle, prove problematic. In our approach, we could avoid this problem by adopting an inner radius of the shell greater than 2R₂₀₀; this choice would decrease the fraction of the mass that has already been within R₂₀₀ but it would yield a shell extending to radii larger than ∼2.5−3R₂₀₀; here, real clusters currently suffer from spectroscopic incompleteness and the estimate of the MAR would be biased. Adopting 2R₂₀₀ as the inner radius of the shell is thus a trade-off between the maximisation the fraction of infalling mass and the current observational limitations.

With our approach, the comparison between the observed MAR and the expectations derived from the N-body simulations is immediate – once exactly the same procedure is applied to the real and mock clusters. Moreover, the identification of the systematic errors is straightforward. On the contrary, making a clean comparison of our estimated MAR with the expectations available in the literature is not possible: the theoretical modelling is based on the mass accretion history of dark matter halos derived from their merger trees (e.g. McBride et al. 2009; Fakhouri et al. 2010; van den Bosch 2002; van den Bosch et al. 2014; Ludlow et al. 2013, 2014, 2016) or from the application of the extended Press-Schechter theory (Bond et al. 1991; Lacey & Cole 1993; Correa et al. 2015). All these models relate the MAR of a cluster of mass M(z) at redshift z to the mass M(z_i) of the same cluster at redshift z_i. For real clusters, the knowledge of the cluster mass at two different redshifts clearly is unavailable.

The technique we present here has two potential observational limitations: the use of N-body simulations for setting the value of the infall velocity and the removal of the real clusters with a sparse spectroscopic sampling when computing the individual MAR. At present, the former issue cannot be solved because the measurement of the radial velocity of cluster galaxies is not currently feasible. We can, in principle, remove the requirement of including the infall velocity from N-body simulations, as well as the assumption of the infall time, by adopting the mass within the shell [2, 2.5]R₂₀₀ as an appropriate prior of the MAR, as Figs. 16 and 14 suggest.

The issue of removing clusters from the samples of clusters with individual MAR concerns the distribution of the galaxies within the R − v_los diagram. The caustics of ∼65% of the removed clusters have unphysical spikes that lead the caustic technique to refrain from estimating the mass at large radii. Usually, a spike emerges when there is an asymmetric lack of galaxies with respect to the line v_los = v_cl in the R − v_los diagram. This lack of galaxies can originate either from a poor spectroscopic sampling or from a real underdense region around the cluster: clusters are generally not spherically symmetric, as assumed by the caustic technique, and asymmetric galaxy distributions on the sky can generate inhomogeneous galaxy distributions in the R − v_los diagrams. For clusters with a spike resulting from poor sampling rather than an actual underdense region, we estimate that the caustic technique could properly locate the caustics if the spectroscopic sampling were increased by ∼15% in a stripe, of the R − v_los diagram, of thickness ∼0.5 h⁻¹ Mpc centred on the spike. In principle, we could thus be able to recover ∼55% of the clusters removed in both CIRS and HeCS. The remaining clusters are either: (1) less massive and thus have a poorly populated R − v_los diagram or (2) embedded in a particularly dense region and the caustic algorithm is unable, at large radii, to properly distinguish between galaxy members and interlopers and, thus, to locate the caustics.

The growth of cosmic structures on linear scales can be different for different dark energy models and modified gravity models (e.g. Burrage et al. 2017; Giocoli et al. 2018; Perenon et al. 2019). In principle, the MAR, which measures the growth on non-linear scales, can be a tool for discriminating among different theories, similarly to the cluster mass function (Kopp et al. 2013; Lombriser et al. 2013; Cataneo et al. 2016; Barreira et al. 2013). However, the sensitivity of the MAR on the theory of gravity remains to be investigated in detail (e.g. Zentner 2007). Similarly, it remains to be seen whether the MAR estimated with the caustic method for clusters at different redshift and with different mass are accurate enough to distinguish among the models. Specifically, we need to assess the optimal balance between the number of observed clusters and the spectroscopy sampling required to reach the necessary sensitivity. However, estimating the mass in the cluster outer regions might not even be necessary for probing the theory of gravity and the velocity field may suffice. For example, in his cosmological N-body simulations in the MOND framework, Candlish (2016) shows that the velocity field in the outskirts of clusters is significantly enhanced compared to standard gravity. We plan to tackle these issues in future work.

7. Conclusion

We use the dense redshift surveys of the CIRS (Rines & Diaferio 2006) and HeCS (Rines et al. 2013) clusters to make the first measurement of the MAR of 129 clusters in the redshift range of 0.01 < z < 0.3 and in the mass range of 10¹⁴ − 10¹⁵ h⁻¹ M_⊙.

Our measurement is based on the spherical accretion recipe suggested by De Boni et al. (2016), where the MAR is estimated from the mass of a shell, of appropriate thickness, surrounding the cluster. We set the inner radius of the shell at 2R₂₀₀. The shell thickness is typically ∼0.5R₂₀₀. To estimate the shell mass at these large radii, we use the caustic technique. When applied to mock redshift surveys of simulated clusters in the same range of redshift and mass, our procedure returns an unbiased MAR within ∼19% of the true MAR, on average.

Our estimates of the MARs of the CIRS and HeCS clusters agree with the MARs of the dark matter halos extracted from a ΛCDM N-body model (Baldi 2012) in the same mass and redshift range. As in the simulations, the observed cluster MARs increase with cluster mass and redshift.

The observed correlations of the MAR with cluster mass and redshift are shared by the mass of the infalling shell. We thus suggest that the mass of the shell can be a proxy of the MAR.

The MAR estimates of the individual clusters have typical uncertainties of ∼17%; for the stacked clusters, when we only stack the members of each individual cluster identified with the caustic technique, the MAR uncertainty is ∼36%. Systematic errors resulting from sample selection could impact these results. We show that neither CIRS nor HeCS have significant photometric incompleteness as a function of radius. However, HeCS misses at most ∼30% of the blue galaxies compared to CIRS; the spectroscopic completeness of both CIRS and HeCS drops to ∼50% at radii larger than ∼2R₂₀₀. Nevertheless, neither the photometric nor the spectroscopic incompletenesses lead to a bias in the MAR estimates. The systematics are within the random errors of the caustic technique.

The measurement of the MAR holds promise for more stringent tests of the standard ΛCDM paradigm. With much larger cluster samples, it may also be able to distinguish between different models of dark matter and of modified gravity. These tests require dense, complete spectroscopic surveys that extend to large cluster-centric radius.

The extension of the cluster survey to a significantly higher redshift provides a broader baseline for these tests. For example, N-body simulations predict that the accretion rate at z = 0.7 is ∼2.6 times larger than at the current epoch (McBride et al. 2009). Both the caustic technique and the weak lensing method can readily be applied to this redshift. Deep photometric surveys, along with such multi-object spectrographs as the Prime Focus Spectrograph on Subaru (Takada et al. 2014) will enable this extension on the observational side. Advances in simulations (Giocoli et al. 2018) will be important for refining the model testing.

Acknowledgments

We sincerely thank the referee whose sensible comments helped us to substantially clarify some relevant issues of our analysis. An early version of this work was the Master thesis in Physics at the University of Torino of SDG. CDB thanks Benedikt Diemer for useful discussion on the splashback radius. The graduate-student fellowship of MP is supported by the Italian Ministry of Education, University and Research (MIUR) under the Departments of Excellence Grant L.232/2016. We also acknowledge partial support from the INFN grant InDark. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

References

All Tables

All Figures

Fig. 1. Schematic figure of the truncated-pyramidal volume of the mock catalogue of a cluster extracted from the simulation. The solid dots show the positions of the dark matter particles in a slice of the simulation box centred on the cluster C. The orange (blue) dots are within (outside) the volume. This figure has only an illustrative purpose: the actual mock catalogues have bL ≫ rFOV and include a substantially larger number of particles than shown here.

In the text

Fig. 2. R − vlos diagram (top panels) and mass profile (bottom panels) of two simulated clusters in the low-mass bin. Left (right) column: cluster at z = 0.12 (z = 0.19). In the bottom panels, the red (blue) solid curves show the caustic (real) mass profile, whereas the red (blue) dot-dashed curves show the NFW fits to the caustic (real) mass profile. In these two examples, the two NFW fits are indistinguishable. The shaded areas show the 50% confidence level of the caustic location and of the caustic mass profile according to the caustic technique recipe. In the R − vlos diagrams the shaded areas are present, but very thin.

In the text

	Fig. 3. Same as Fig. 2 but for two simulated clusters in the high-mass bin.
In the text

	Fig. 4. Mean (blue) and median (orange) profiles of the radial velocity of the particles within the dark matter halos extracted from the simulation for the two mass bins at z = 0 and z = 0.44, as indicated in the panels. 68% of the profiles of the individual halos lie within the light blue areas.
In the text

Fig. 5. MAR of simulated clusters in the low- (lower set of points) and high- (upper set of points) mass bins. The blue squares and the red triangles show the median MAR based on the three-dimensional and the caustic mass profiles, respectively. The blue shaded areas show the 68th percentile range of the distribution of the MAR derived from the three-dimensional mass profiles; the red error bars show the 68th percentile ranges of the estimates obtained with the caustic mass profiles.

In the text

	Fig. 6. Luminosities of the CIRS (blue points) and HeCS (orange points) clusters as a function of redshift. The superimposed curves show the flux and redshift selection functions of the two catalogues. Filled points refer to the clusters where we can compute an individual MAR.
In the text

	Fig. 7. Mass distribution of the CIRS (blue histogram) and HeCS (orange histogram) clusters. The black and red dashed lines show the median masses for the CIRS and HeCS catalogues, respectively. The area under each histogram is normalised to unity.
In the text

Fig. 8. (g − r) colour-magnitude diagrams of the CIRS (left) and HeCS (right) galaxies, including k-corrections. The galaxies have a cluster-centric distance smaller than 3R200 and line-of-sight velocity |vlos − vcl|< 3000 km s−1. The black solid curve in the right panel shows the fit by Rines et al. (2013) derived from the HeCS galaxies. The black dashed, dotted, and dash-dotted lines show the ±0.30, ±0.40, ±0.50 shifts of the same fit. In the left panel, the black solid line shows the Rines et al.’s fit with an offset of −0.03 mag. The black dashed, dotted, and dash-dotted lines show the ±0.30, ±0.40, ±0.50 mag shifts of this line.

In the text

Fig. 9. Left panel: ratio of the numbers of galaxies with photometric data in the range r/R200 ∈ ]2, 2.5] and within r/R200 = 1 within each cluster, against the cluster redshift. Right panel: ratio between the numbers of galaxies with spectroscopic redshifts within the same regions. Blue and orange points show the CIRS and HeCS clusters, respectively. Solid points show the clusters for which we estimate the individual MAR.

In the text

	Fig. 10. Spectroscopic completeness of the CIRS (blue points) and HeCS (orange points) clusters, against their redshifts, within r/R200 = 1 from the cluster centre (upper panel) and within the range r/R200 ∈ ]2, 2.5] (lower panel). Solid points show the clusters for which we estimate the individual MAR.
In the text

	Fig. 11. Distributions of the individual MARs. The blue, orange and green histograms refer to the CIRS, low-mass HeCS, and high-mass HeCS clusters, respectively. The dashed lines with the same colours show the median MAR of each sample. The area under each histogram is normalised to unity.
In the text

Fig. 12. Three examples of the estimation of the MAR for individual clusters in CIRS. The left column shows the R − vlos diagram for each cluster. The right column shows the corresponding caustic mass profile (solid line) and the NFW fit (dot-dashed line). The black and red dashed lines show the inner and outer radius of the infalling spherical shell. The M200 of the clusters increases from top to bottom.

In the text

	Fig. 13. Same as Fig. 12 for three HeCS clusters. Zw3179 and A655 belong to the low-mass sample, whereas A963 is in the high-mass sample.
In the text

	Fig. 14. MAR of individual clusters as a function of their mass M200; the colour code shows the dependence on redshift. Filled circles, squares and triangles refer to the CIRS, low- and high-mass HeCS sample, respectively. The open circles show the median MARs of the N-body clusters of our two simulated samples; for these clusters, we estimate the MAR from their three-dimensional mass profiles.
In the text

Fig. 15. Median MARs of real (green) and simulated (blue and red) clusters. The green error bars show the 68th percentile ranges of the MAR and redshift distributions of the real clusters. The value close to each green filled circle is the median mass M200 of the real clusters in units of h−1 M⊙. The blue squares and red triangles are for the MAR based on the three-dimensional and the caustic mass profiles of simulated clusters, respectively (from Fig. 5). The light-blue shaded areas and the red error bars are the 68th percentile range of the distribution of the individual MAR of the simulated clusters. The open circles show the median MARs of real clusters for different values of the infall velocity. We show the results where we changed the initial infall velocity that we adopted by: −40% (blue), −20% (cyan), +20% (orange), and then +40% (red). The black circles show the MAR estimated with an initial infall velocity of vi = 0.

In the text

Fig. 16. Mass M2−2.5 of the shell with radii 2R200 and 2.5R200 as a function of M200. The colour code shows the dependence on redshift. Circles, squares and triangles refer to the CIRS, low- and high-mass HeCS sample, respectively. The open circles show the median M2−2.5 of the N-body clusters of our two simulated samples, estimated from their three-dimensional mass profiles. The crosses show the median M2−2.5 obtained from the NFW fits to the three-dimensional mass profiles: on average, they overestimate the individual true mass of the shell by ∼22−28%, depending on the mass bin.

In the text

Fig. 17. R − vlos diagrams and corresponding mass profiles of the two stacked clusters from the CIRS catalogue. The upper (lower) panels refer to the low-mass (high-mass) subsample. The black and red dashed vertical lines show the inner and outer radius of the spherical shell used to estimate the MAR. The shaded areas show the 50% confidence level of the caustic location and of the caustic mass profile according to the caustic technique recipe. The dot-dashed lines show the NFW fits to the caustic mass profiles; the green vertical dotted lines show the upper limit of the radial range over which the NFW fits have been performed.

In the text

	Fig. 18. Same as Fig. 17 for the HeCS catalogue.
In the text

Fig. 19. MARs of the stacked clusters of CIRS (the two orange circles at low redshift) and HeCS (the two orange circles at high redshift). The median MARs of the individual clusters from Fig. 15 are shown in green. The value close to each circle is the mass M200 in units of h−1 M⊙. The blue squares with the shaded bands are from Fig. 5 and show the median MARs and their spreads of the simulated clusters estimated with the three-dimensional mass profiles.

In the text