THE MASS ACCRETION RATE OF GALAXY CLUSTERS: A MEASURABLE QUANTITY

In the current model of the formation of cosmic structure, where dark matter halos form from the aggregation of smaller halos, the mass accretion of dark matter halos is a stochastic process whose average behavior can be predicted with N-body simulations and semi-analytical models (van den Bosch 2002; Zhao et al. 2003b; Sheth & Tormen 2004a, 2004b; Giocoli et al. 2007; McBride et al. 2009; Zhao et al. 2009). This process is generally investigated with the identification of the merger trees of dark matter halos, enabling the study of the mass accretion history (MAH) and the mass accretion rate (MAR) as a function of redshift z (van den Bosch 2002; Fakhouri & Ma 2008; McBride et al. 2009; Fakhouri et al. 2010; Giocoli et al. 2012).

Observationally, the exploration of the MAR of individual dark matter halos has only been attempted on the scales of galaxies with the galaxy–galaxy merger rate: one usually combines the number of observed pairs of close or disturbed galaxies with the theoretical merger probability and time scale (Lotz et al. 2011; Jian et al. 2012; Casteels et al. 2014). However, different investigations reach discrepant conclusions (Lotz et al. 2011) because the merger rate of dark matter halos is not identical to the merger rate of galaxies (Fakhouri & Ma 2008; Guo & White 2008; Moster et al. 2013) and the two rates are related by dissipative processes which are difficult to model (Hopkins et al. 2013).

In contrast, the measurement of the MAR of galaxy clusters can in principle be based on the estimated amount of mass in the cluster's surrounding regions, where we can safely neglect the dissipative processes which affect the galaxy–galaxy merger rate. Nevertheless, no measurements of the MAR of galaxy clusters have been attempted so far. This observational deficiency is due to the fact that in the large and less dense outer regions of clusters, galaxy members are difficult to distinguish from foreground and background galaxies; other probes, e.g., X-ray emission, are below the sensitivity of current instruments. In addition, the outer regions of clusters are not in dynamical equilibrium and therefore the usual mass estimation methods based on virial equilibrium are inappropriate. Considering this picture, one may think that the MAR predictions from N-body simulations are not capable of being tested.

Here we take a more optimistic perspective and explore the possibility of estimating the MAR of galaxy clusters by measuring the mass of a spherical shell surrounding the cluster. The thickness of this shell depends on the assumed infall time, on the radius at which the infall happens and on the initial velocity of the falling mass. This translates into a change of the shell thickness with redshift. Albeit rather crude when compared with the stochastic aggregation of dark matter halos in the hierarchical clustering formation model, this approach would provide a method to estimate the MAR that depends on the cluster mass profile at radii larger than the virial radius.

In theoretical investigations, the relation between the mass density profiles of galaxy clusters and their accretion history is known. For example, Ludlow et al. (2013) find that the inner part of a halo retains the information on how the halo has accreted its mass through a correlation between the mean inner density within the scale radius r_s and the critical density of the universe at the time when the mass of the main progenitor is equal to M(<r_s). Correa et al. (2015a, 2015b) confirm these findings and demonstrate that the MAH can be expressed with a general formula similar to the one originally proposed by Tasitsiomi et al. (2004) for a ΛCDM model and widely studied in McBride et al. (2009) based on N-body simulations: the MAH has an exponential evolution with redshift in the high z regime and follows a power law at low z when the accelerated expansion of the universe freezes the growth of perturbations.

Diemer & Kravtsov (2014) find that the steepness of the slope of the outer halo density profile increases with increasing MAR. In addition, the central concentration is anti-correlated with the MAR. Adhikari et al. (2014) note that the location where the steepening of the slope is observed corresponds to the radius associated to the splashback of the material that the halo has recently accreted. van den Bosch et al. (2014) show that the growth of the central potential precedes the assembly of mass and introduce a formula that can be used to compute the average MAH in any ΛCDM cosmology without running numerical simulations.

Here we derive a simple relation between the mass profile of a dark matter halo and its MAR derived from the halo merger tree extracted from N-body simulations. This result is relevant because it implies that in principle we can estimate the MAR of galaxy clusters from the estimate of the mass profile in their outer regions. This measurement can be performed with the caustic technique (Diaferio & Geller 1997; Diaferio 1999; Serra et al. 2011) which is not affected by the presence of substructures, the non-equilibrium state of the cluster outer region, and the correlated large-scale structures along the line of sight (Diaferio 1999; Serra et al. 2011; Geller et al. 2013). The caustic technique only requires a sufficiently dense galaxy redshift survey in the cluster field of view to return a mass estimate accurate to 20% on average (Serra et al. 2011).

Estimating the MAR also requires the knowledge of the infall velocity of the shell. This quantity is inaccessible to observations and remains a free parameter; we adopt the mean infall velocity of a shell surrounding dark matter halos in N-body simulations as illustrated below. Although it prevents us from estimating the MAR of individual clusters, this choice allows us to cope with the unmeasurable velocity of the falling shell. Consequently, the technique we propose here aims at estimating the mean accretion rate of a sample of galaxy clusters rather than the accretion rate of individual clusters.

We investigate the feasibility of our approach by considering dark matter halos at redshift z < 2, with mass comparable to the clusters in the CIRS and HeCS catalogs (Rines & Diaferio 2006; Rines et al. 2013) whose outer mass profiles have already been measured with the caustic technique. Our interest in investigating the growth of structures at nonlinear scales of galaxy clusters from an observational perspective, which is still relatively poorly explored (e.g., Lemze et al. 2013), perfectly complements most of the current efforts that focus on constraining the growth factor in the linear and mildly nonlinear regimes with large-scale redshift surveys and weak-lensing tomography (e.g., Euclid, DES, eBOSS, DESI, PFS, LSST, and WFIRST).

In Section 2 we introduce our spherical infall prescription. In Section 3 we discuss the properties of the CoDECS set of N-body simulations (Baldi 2012). We discuss our results in Section 4.

Our spherical infall prescription assumes a shell of matter falling onto the enclosed halo. We aim to quantify $\dot{M}={dM}/{dt}$, where dM is the mass of a spherical shell of thickness ${\rm{\Delta }}r$ and proper radii R_i and ${R}_{i}+{\rm{\Delta }}r$ and dt is the time it takes to fall to R_i. R_i is the radius at which we consider the infall to happen. We choose ${\rm{\Delta }}r={\delta }_{s}{R}_{i}$ where ${\delta }_{s}$ is a free parameter. By solving the equation of motion ${d}^{2}r/{{dt}}^{2}={a}_{0}$ under the assumption of constant infall acceleration, ${a}_{0}=-{GM}{(\lt {R}_{i})/({R}_{i}+{\delta }_{s}{R}_{i}/2)}^{2}$, and initial velocity v_i, we obtain the infall time ${t}_{{\rm{inf}}}$ from the equation ${a}_{0}{t}_{{\rm{inf}}}^{2}/2+{v}_{i}{t}_{{\rm{inf}}}=-{\delta }_{s}{R}_{i}/2$, where we consider the shell to be accreted when the shell middle point, initially at ${\delta }_{s}{R}_{i}/2$, reaches R_i; we obtain

$$\begin{eqnarray}&&{t}_{{\rm{inf}}}^{2}{GM}{(\lt {R}_{i})-{t}_{{\rm{inf}}}2{R}_{i}^{2}(1+{\delta }_{s}/2)}^{2}{v}_{i}+\\ &&\quad -{R}_{i}^{3}{\delta }_{s}{(1+{\delta }_{s}/2)}^{2}=0.\end{eqnarray} \tag{ 1 }$$

We can now express the MAR obtained from our prescription as

$$\begin{eqnarray}&&\dot{M}=\displaystyle \frac{{dM}}{{dt}}\equiv \displaystyle \frac{{M}_{{\rm{shell}}}}{{t}_{{\rm{inf}}}}.\end{eqnarray} \tag{ 2 }$$

This recipe has three input parameters: the scale R_i that defines the infall radius, the initial velocity v_i, and the thickness ${\delta }_{s}$ of the shell. It is more convenient to use the infall time as an input parameter and to derive the shell thickness ${\delta }_{s}$. In this case, Equation (1) reads:

$$\begin{eqnarray}&&{\delta }_{s}^{3}\displaystyle \frac{{R}_{i}^{3}}{4}+{\delta }_{s}^{2}\left({R}_{i}^{3}+\displaystyle \frac{{R}_{i}^{2}}{2}{v}_{i}{t}_{{\rm{inf}}}\right)+{\delta }_{s}({R}_{i}^{3}+2{R}_{i}^{2}{v}_{i}{t}_{{\rm{inf}}})\\ &&\quad ={GM}(\lt {R}_{i}){t}_{{\rm{inf}}}^{2}-2{R}_{i}^{2}{v}_{i}{t}_{{\rm{inf}}}.\end{eqnarray} \tag{ 3 }$$

Therefore, by choosing ${t}_{{\rm{inf}}}$ we derive the thickness ${\delta }_{s}{R}_{i}$ of the shell centered on the cluster and we can estimate the MAR of a galaxy cluster from its mass profile by using Equation (2).

However, accretion is a stochastic process, and we need to verify that our simple approach is capable of correctly estimating the actual MAR. Below, we compare the MAR estimated with our recipe with the MAR of dark matter halos derived from their halo merger trees obtained in an N-body simulation.

CoDECS (Coupled Dark Energy Cosmological Simulations) is a suite of N-body simulations in different cosmological models (see Baldi 2012 for further details). Here we use the L-CoDECS simulation of a ΛCDM model, a collisionless N-body simulation of a flat universe, with the following cosmological parameters consistent with the WMAP7 data (Komatsu et al. 2011): cosmological matter density ${{\rm{\Omega }}}_{m0}=0.226$, cosmological constant ${{\rm{\Omega }}}_{{\rm{\Lambda }}0}=0.729$, baryonic mass density ${{\rm{\Omega }}}_{b0}=0.0451$, Hubble constant h = 0.703, power spectrum normalization ${\sigma }_{8}=0.809$, and power spectrum index n_s = 0.966. The simulated box has a comoving volume of ${(1{\rm{Gpc}}{h}^{-1})}^{3}$ containing (1024)³ dark matter (CDM) particles and the same amount of baryonic particles. The mass resolution is ${m}_{{\rm{DM}}}=5.84\times {10}^{10}\;{M}_{\odot }\;{h}^{-1}$ for CDM particles and ${m}_{b}=1.17\times {10}^{10}\;{M}_{\odot }\;{h}^{-1}$ for baryonic particles. No hydrodynamics are included in the simulation. Baryonic particles are only included to account for the different forces acting on baryonic matter in the coupled quintessence models.

3.1. CoDECS Merger Trees

3.2. Cluster Profiles

We plan to apply our recipe for the estimate of the MAR to redshift surveys of galaxy clusters by estimating their mass profile with the caustic method (Diaferio & Geller 1997; Diaferio 1999; Serra et al. 2011). The caustic method returns mass profiles that are affected by 20%–50% uncertainties when applied to clusters with M₂₀₀ of ${10}^{14}\;{M}_{\odot }\;{h}^{-1}$ or larger. When applied to less massive clusters, the systematic errors introduced by the caustic technique can become substantially larger (Rines & Diaferio 2010; Serra et al. 2011). Therefore, with this observational perspective in mind we concentrate on halos in the M₂₀₀ mass range ${10}^{14}\mbox{--}{10}^{15}\;{M}_{\odot }\;{h}^{-1}$, which corresponds to the most common mass of the clusters in catalogs like CIRS (Rines & Diaferio 2006) and HeCS (Rines et al. 2013). We thus define two mass bins at z = 0 and follow the evolution with redshift of the clusters assigned to these bins. We select halos with median ${M}_{200}={10}^{14}$ and ${10}^{15}\;{M}_{\odot }\;{h}^{-1}$ at z = 0. The low-mass bin contains 2000 objects at z = 0, while the high-mass bin is limited to 50 objects: there are 36 clusters with ${M}_{200}\gt {10}^{15}\;{M}_{\odot }\;{h}^{-1}$ in the simulated box, but we removed the 11 most massive clusters to have a mean mass similar to the median mass in the bin. Table 1 lists the mean, standard deviation, median, 68%, and 90% of each mass bin at z = 0.

Table 1.  Mean, Median, Standard Deviation and Percentile Ranges of M₂₀₀ and Number of Halos in the Two Mass Bins at z = 0

Mean	σ	Median	68%	90%	Number of Halos
${M}_{200}\;[{M}_{\odot }\;{h}^{-1}]$
1.00 × 1014	2.90 × 1012	1.00 × 1014	(0.97–1.04) × 1014	(0.96–1.05) × 1014	2000
1.04 × 1015	1.26 × 1014	1.00 × 1015	(0.91–1.19) × 1015	(0.88–1.24) × 1015	50

Download table as:  ASCIITypeset image

For each halo in the two mass bins and for each progenitor at higher redshift (as defined in Section 3.1), we evaluate the mass profile and the profile of the radial velocity ${v}_{{\rm{rad}}}$ up to 10R₂₀₀. We show the radial velocity profile in the two bins at different redshifts in Figures 1 and 2.

Zoom In Zoom Out Reset image size

With the radial velocity profile we can identify three regions: an internal region with ${v}_{{\rm{rad}}}\simeq 0$, where matter is orbiting around the center of the cluster; an infall region, where ${v}_{{\rm{rad}}}$ becomes negative and indicates an actual infall of matter toward the center of the cluster; and a Hubble region at very large radii, where ${v}_{{\rm{rad}}}$ becomes positive and the Hubble flow dominates. Broadly speaking, the infall radius ${R}_{{\rm{inf}}}$, i.e., the radius where the minimum of ${v}_{{\rm{rad}}}$ occurs, is between 2R₂₀₀ and 3R₂₀₀, independent of mass and redshift.

More et al. (2015) suggest the use of the splashback radius as the physical halo boundary because it separates the infall region from the region where the matter has already been accreted. The splashback radius is defined as the outermost radius reached by accreted material in its first orbit around the cluster center. Its exact location depends both on redshift and MAR, but it is in general larger than R₂₀₀. Noticeably, they also show that the splashback radius is not affected by the evolution of the critical density of the universe, unlike the usual R₂₀₀, whose definition is based on the average overdensity of the halo at a given redshift. By adopting R₂₀₀, part of the evolution of the halo properties with redshift simply is a consequence of the evolution of the critical density of the universe; this effect generates a so-called pseudo-evolution of the halo properties (see Diemer et al. 2013). This pseudo-evolution substantially disappears when we adopt the splashback radius. This use of this radius might thus be more preferable in the investigations of the redshift evolution of the halo properties.

Our results in Figures 1 and 2 are consistent with the results of More et al. (2015). For massive objects at z = 0, the splashback radius ${R}_{{\rm{sp}}}$ is close to 2R₂₀₀ and the infall radius is ${R}_{{\rm{inf}}}\sim 1.4{R}_{{\rm{sp}}}$. Given that  $M(\lt {R}_{{\rm{sp}}})$ and $M(\lt {R}_{{\rm{inf}}})$ are not affected by the pseudo-evolution mentioned above, we use ${R}_{i}=2{R}_{200}\sim {R}_{{\rm{sp}}}$ as the radius at which we consider the infall to happen in our spherical infall prescription and ${v}_{i}={v}_{{\rm{shell}}}\sim {v}_{{\rm{rad}}}({R}_{{\rm{inf}}})$ (see Equation (3)).

Zoom In Zoom Out Reset image size

For each halo at z = 0 in the two mass bins we build the MAH at 2R₂₀₀ (Figure 3). In Table 2 we list the mean, standard deviation, median, 68%, and 90% of M(<2R₂₀₀) at z = 0 for the two mass bins. In Figure 3 we also show the MAH model by Giocoli et al. (2013) (see Appendix) rescaled to 2R₂₀₀. We obtain this rescaling by extending the mass density profile to 2R₂₀₀ using the NFW functional and adopting the redshift-dependent relation between the concentration and mass of Zhao et al. (2003a) modified for M₂₀₀ by Giocoli et al. (2013). For a direct comparison with previous results in the literature we also estimate the MAH at the standard R₂₀₀. We report this analysis and the comparison with previous work in the Appendix.

Zoom In Zoom Out Reset image size

Table 2.  Mean, Median, Standard Deviation, and Percentile Ranges of M(<2R₂₀₀) in Each Mass Bin at z = 0

Mean	σ	Median	68%	90%
$M(\lt 2{R}_{200})\;[{M}_{\odot }\;{h}^{-1}]$
1.60 × 1014	1.73 × 1013	1.58 × 1014	(1.45–1.74) × 1014	(1.37–1.90) × 1014
1.55 × 1015	2.41 × 1014	1.56 × 1015	(1.32–1.82) × 1015	(1.19–2.00) × 1015

Download table as:  ASCIITypeset image

We find the ratio $M(\lt 2{R}_{200})/{M}_{200}$ to be ∼1.6, similar to the result of More et al. (2015) for ${M}_{{\rm{sp}}}/{M}_{200}$ in massive objects at z = 0. In real observations measuring the infall radius of a cluster, namely where the infall velocity reaches its minimum, is currently unfeasible so we keep fixed ${R}_{i}=2{R}_{200}$ for all masses at all redshifts. This choice clearly is an oversimplification, but Figures 1 and 2 indicate that this assumption is reasonable. Indeed ${R}_{{\rm{inf}}}$ lies within 2R₂₀₀ and 3R₂₀₀ for a wide range of masses and redshifts.

When we observe a cluster we look at a particular instant of its evolution and we can trace the evolution of its mass neither backward nor forward in time. Therefore, the MAH is not a measurable quantity. In contrast, the MAR can in principle be estimated as the ratio of the mass that is being accreted by the cluster at a given time and its infall time as described by our spherical infall prescription in Section 2, Equation (2). This estimate depends on the external mass profile and on the mean infall velocity. The mass profile is measurable with the caustic method, whereas the initial velocity v_i (see Equation (3)) remains a free parameter which can be inferred from N-body simulations.

As already stated in Section 3.2 we choose ${R}_{i}=2{R}_{200}\sim {R}_{{\rm{sp}}}$ as the radius where the infall occurs. As initial velocity v_i for the infall, for each redshift in each mass bin we adopt the mean velocity in our first radial bin [2–2.5] ${R}_{200}\sim {R}_{{\rm{inf}}}$ (see Figures 1 and 2); this velocity ranges from −200 to $-250\;\mathrm{km}\;{{\rm{s}}}^{-1}$. We call this velocity ${v}_{{\rm{shell}}}$. The value of ${v}_{{\rm{rad}}}$ in the radial bin [2.5–3] R₂₀₀ roughly remains in the same range and adopting this radial bin instead of the first bin has negligible effects on the final results. Although our prescription for the MAR estimate clearly depends on the choice of both R_i and ${v}_{{\rm{shell}}}$, using the mean infall velocity instead of the infall velocity of each halo enables the design of a feasible procedure for the observational estimate of the MAR of clusters.

The spread of the distribution of ${v}_{{\rm{shell}}}$, estimated as $\sigma ={\sigma }_{{v}_{i}}/\sqrt{N}$ where N is the number of halos used to evaluate ${v}_{{\rm{shell}}}$, is smaller than the spread of the distribution of the individual v_i's, ${\sigma }_{{v}_{i}}$ shown in Figures 1 and 2. For the ${10}^{14}{M}_{\odot }{h}^{-1}$ mass bin, the $1\sigma$ relative spread of ${v}_{{\rm{shell}}}$ is smaller than 2% and propagates into a relative spread on the MAR of 5% at most, which is well within 68% of the MAR distribution obtained with ${v}_{{\rm{shell}}}$. Similarly, for the ${10}^{15}{M}_{\odot }{h}^{-1}$ mass bin, where N drops from 2000 to 50, the $1\sigma$ relative spread of ${v}_{{\rm{shell}}}$ is $\sim 14\%$, implying a ∼35% relative spread on the MAR. Incidentally, even by using the $1\sigma$ spread of the distribution of the individual v_i's, the resulting MAR's are within 68% of the distribution of the MAR's obtained with ${v}_{{\rm{shell}}}$ for both mass bins. The different spread of ${v}_{{\rm{shell}}}$ in the two mass shells deriving from the different size of the two cluster samples indicates that the expected uncertainty on the MAR estimated with real data will vary substantially depending on the number of observed clusters, as one can expect. We plan to fully quantify this uncertainty in future work by estimating the MAR of simulated clusters from mock redshift surveys.

For the infall time, which is the last parameter of the model, we choose ${t}_{{\rm{inf}}}={10}^{9}\;{\rm{year}}$. This value is suggested by the redshift-independent $\sim 1\;{\rm{Gyr}}$ time step of the snapshots of the simulation used to estimate the MAR from the merger trees; it also has the advantage of being similar to the dynamical time, simply defined as ${t}_{{\rm{dyn}}}\sim R/\sigma$, for the clusters of our analysis. In fact, for the ${10}^{14}\mbox{--}{10}^{15}\;{M}_{\odot }\;{h}^{-1}$ clusters at z = 0 we consider here, $R\sim 1\;{\rm{Mpc}}$ and $\sigma \sim 1000\;\mathrm{km}\ {{\rm{s}}}^{-1}$, and ${t}_{{\rm{dyn}}}\sim {10}^{9}\;{\rm{year}}$. For the progenitors of these clusters at higher redshifts, which have masses at most a factor 10 smaller, the velocity dispersion σ is 10^1/2 ∼ 3 times smaller, and the virial radius R is smaller by roughly a similar factor. Therefore, ${t}_{{\rm{dyn}}}$ remains basically constant and equal to ${10}^{9}\;{\rm{year}}$. Finally, this equality between ${t}_{{\rm{dyn}}}$ and the snapshot time interval prevents us from assuming ${t}_{{\rm{inf}}}$ very different from $1\;{\rm{Gyr}}$; if ${t}_{{\rm{inf}}}$ departs too much from the snapshot time interval, the comparison of our estimated MAR's with the MAR's extracted from the merger trees would be inappropriate. We therefore investigate the dependence of our results on ${t}_{{\rm{inf}}}$ within ∼20% of ${10}^{9}\;{\rm{year}}$ and find that our results remain unaffected.

Once R_i, ${v}_{{\rm{shell}}}$ and ${t}_{{\rm{inf}}}$ are specified, the model is completely determined by Equation (3). For each halo in the two mass bins and for each progenitor at higher redshift, we evaluate the thickness ${\delta }_{s}$ of the infalling shell and its mass. We show the evolution with redshift of the shell thickness ${\delta }_{s}$ in Figure 4. The shell thickness increases with increasing redshift, and the intrinsic scatter of the distribution also enlarges. This fact reflects the individual evolution of M(<2R₂₀₀) as shown in Figure 3. The solid blue line in Figure 4 marks the value ${\delta }_{s}=0.5$ for which the external radius of the shell is equal to 3R₂₀₀, close to the cluster turnaround radius. This value of ${\delta }_{s}$ is reached between redshift 1 and 1.5, depending on the mass of the cluster.

Zoom In Zoom Out Reset image size

The accretion onto a cluster is a highly anisotropic process; nevertheless, we are confident that, given the thickness of the shell, we are taking into account almost all of the mass that is actually falling in the time interval ${t}_{{\rm{inf}}}.$

Figure 5 shows the MAR of all the clusters in the two mass bins estimated with Equation (2). It also shows the MAR derived from the merger trees of the halos.

Zoom In Zoom Out Reset image size

We see that the mean and median results from the merger trees lie within the region defined by the 68% range of the distribution of the MAR's obtained from Equation (2) for both mass bins. For the ${10}^{15}\;{M}_{\odot }\;{h}^{-1}$ bin, the mean and median values of our estimate are ∼40% smaller than the average MAR values from the merger trees. This systematic underestimate might be due to a statistical fluctuation because the sample only contains 50 clusters compared with the 2000 clusters of the less massive bin. In fact, a similar underestimate is observed when we compare the MAH from the merger trees with the Giocoli model in the right panel of Figure 3. In contrast, for the ${10}^{14}\;{M}_{\odot }\;{h}^{-1}$ mass bin, the mean and median MAR from our prescription recover the merger tree results within 20% in the redshift range z = [0, 2]. Our results are relevant because they show that our simple spherical infall prescription can in principle provide a method to estimate the average MAR of galaxy clusters from redshift surveys. In future work, we will apply our prescription to synthetic redshift surveys of galaxy clusters to quantify the uncertainties and possible systematic errors of our procedure.

Clearly Figure 5 only compares the average MAR obtained from the merger trees of individual halos with the average MAR provided by our spherical infall technique. Our recipe was not conceived to completely capture all the features of the MAR derived by the complex merging process of individual halos. Nevertheless, the average of the MAR of individual halos still is satisfactorily estimated by our recipe.

Figure 6 shows the distribution of the ratio between the MAR estimated with our recipe and the MAR derived from the merger tree for each individual halo of the ${10}^{14}\;{M}_{\odot }\;{h}^{-1}$ mass bin at four different redshifts. The distribution has a tail toward large values. Remarkably the median value of this ratio is close to the ratio between the average MAR from our model and the MAR from the merger trees (red line). The mean clearly is larger because of the tail of high values. The 68% ranges from 0.2 to 2.2. This result confirms that with our approach we can measure the mean MAR but not the MAR of the individual clusters, which is affected by lack of spherical symmetry and large variations of the infall velocity.

Zoom In Zoom Out Reset image size

Figures 1–3 show the intrinsic halo-by-halo scatter both in radial velocity and mass. Our choice to use the same ${v}_{{\rm{shell}}}$ for all clusters in a given mass bin at a given redshift implies that we neglect the scatter that originates the spread of the distribution shown in Figure 6. However, this choice keeps the model relatively simple and applicable to real clusters. It is worth saying that the impact of the large-value tail is reduced if we take the ratio of the averages of the MAR of each individual halo estimated with our prescription and with the merger trees, rather than the average of the ratio. The ratio of averages is shown with the red lines in Figure 6 and it corresponds to the result shown in Figure 5. The remarkable and encouraging result of our analysis is that the agreement shown in Figures 5 and 6 is obtained without requiring any input information from the merger trees.

We investigate the feasibility of directly measuring the mean MAR of a sample of galaxy clusters from their mass profile. To measure the mean MAR we suggest a prescription based on the infall of a spherical shell of matter onto the halo, with constant acceleration and initial velocity derived from the average infall velocity of matter around dark matter halos in N-body simulations. Once we fix the scale R_i, that defines the halo radius at which the infall occurs, the initial velocity ${v}_{{\rm{shell}}}$, and the infall time ${t}_{{\rm{inf}}}$ our method only depends on the mass profile at large radii, beyond R_i.

We consider dark matter halos from the CoDECS set of N-body simulations and compare their MAR estimated with our prescription with the MAR estimated from the merger trees extracted from the simulations. We focus on two sets of halos with mass M₂₀₀ around 10¹⁴ and ${10}^{15}\;{M}_{\odot }\;{h}^{-1}$ at z = 0.

We recover the mean and median MAR obtained from the merger trees without bias and with 20% accuracy in the redshift range z = [0, 2] for the ${10}^{14}\;{M}_{\odot }\;{h}^{-1}$ mass bin. The accuracy is about 40%, with a systematic underestimation of ∼40%, for the ${10}^{15}\;{M}_{\odot }\;{h}^{-1}$ mass bin. This result is impressive given the simple assumptions of our prescription and the fact that no input parameter of the model is taken from the merger trees. Our result does show that measuring the mean MAR of a sample of real galaxy clusters is in principle possible.

A fundamental step to assess the feasibility of our approach is to apply the caustic technique to realistic mock redshift surveys of galaxy clusters extracted from N-body simulations and quantify the accuracy of the estimated MAR. This investigation remains for further studies. Similarly the analysis of the dependence of the method parameters and its results on the cosmological model and on the theory of gravity remain to be investigated. Specifically, ${v}_{{\rm{shell}}}$ might turn out to vary substantially with the assumed cosmology and thus to be a crucial parameter of the method.

We might expect that the accuracy will be better than 20% on the average MAR if we estimate the MAR of individual clusters in a given mass bin and estimate their average MAR; this approach would minimize the systematic errors due to projection effects which dominate the estimate of the mass profile with the caustic technique. Rines et al. (2013) have already applied this approach to measure the total mass of clusters within their turnaround radius, the ultimate mass ${M}_{{\rm{ta}}}$. By combining 50 CIRS clusters with 58 HeCS clusters, they find ${M}_{{\rm{ta}}}/{M}_{200}=1.99\pm 0.11$, a measure accurate to 5% and in agreement with the ΛCDM prediction where ${M}_{{\rm{ta}}}/{M}_{200}$ has a log-normal distribution with a peak at mass ratio 2.2 and dispersion 0.38 (Busha et al. 2005).

Our measurements of the average MAR may provide an additional tool to discriminate among different cosmological models if deviations from the ΛCDM model generate measurable differences in the MAR.

We are deeply grateful to the anonymous referee whose very accurate comments and suggestions urged us to substantially revise our analysis which led to more solid and valuable results. We also sincerely thank Fabio Fontanot for providing us with the merger trees of the H-CoDECS simulation. We thank Aaron Ludlow, Margherita Ghezzi, Giulio Falcioni, and Andrea Vittino for useful discussions. C.D.B., A.L.S., and A.D. acknowledge partial support from the grant Progetti di Ateneo/CSP TO_Call2_2012_0011 "Marco Polo" of the University of Torino, the INFN grant InDark, and the grant PRIN_2012 "Fisica Astroparticellare Teorica" of the Italian Ministry of University and Research. A.L.S. also acknowledges the support of Dipartimento di Fisica, University of Torino, where most of this project was carried out. C.G.'s research is part of the project GLENCO, funded under the European Seventh Framework Programme, Ideas, Grant Agreement n. 259349. G.C. also thanks CNES for financial support. M.B. is supported by the Marie Curie Intra-European Fellowship "SIDUN" within the 7th Framework Programme of the European Commission. This work was partially supported by grants from Région Ile-de-France.

In this appendix, we investigate the MAH of M₂₀₀ of dark matter halos: we compare the merger tree results with two known fitting formulae and a theoretical model. We calculate the mean and median MAH for the objects in four mass bins (Table 3). The two most massive bins are the same used in the main body of the paper. The four bins are centered on ${M}_{200}={10}^{15},{10}^{14},{10}^{13},5\times {10}^{12}\;{M}_{\odot }\;{h}^{-1}$. Each bin contains 2000 objects at z = 0 with the exception of the largest-mass bin whose sample is limited to 50 objects. The lowest-mass bin is centered on ${M}_{200}=5\times {10}^{12}\;{M}_{\odot }\;{h}^{-1}$ because ${M}_{200}={10}^{12}\;{M}_{\odot }\;{h}^{-1}$ is below the resolution limit of about 100 particles per subhalo set by SUBFIND. Table 3 also lists the 90% mass range of each mass bin.

Table 3.  Mean, Median, Standard Deviation, and Percentile Ranges of M₂₀₀ and Number of Halos in Each Mass Bin at z = 0

Mean	σ	Median	68%	90%	Number of Halos
${M}_{200}\;[{M}_{\odot }\;{h}^{-1}]$
1.04 × 1015	1.26 × 1014	1.00 × 1015	(0.91–1.19) × 1015	(0.88–1.24) × 1015	50
1.00 × 1014	2.90 × 1012	1.00 × 1014	(0.97–1.04) × 1014	(0.96–1.05) × 1014	2000
1.00 × 1013	2.36 × 1010	1.00 × 1013	(1.00–1.00) × 1013	(1.00–1.00) × 1013	2000
5.00 × 1012	6.04 × 109	5.01 × 1012	(5.00–5.01) × 1012	(5.00–5.01) × 1012	2000

Download table as:  ASCIITypeset image

Figure 7 shows the MAHs of the four mass bins. We limit our study to the low redshift range z = [0, 2] because we are interested in the observational relevance of our analysis. As we can see in Figure 7, for the two largest-mass bins the mean MAH agrees with the median MAH within 20%. In the two smallest-mass bins the difference between the mean and the median MAH is never larger than 5%. In all four cases the standard deviation and the 68% are comparable. The results from the largest-mass bin are noisier because of the low-number statistics. The number of objects ${N}_{{\rm{hal}}}$ at each z decreases with increasing z due to the resolution limit: not all the objects selected at z = 0 have merger trees that reach z = 2. Indeed, the decrease is larger for less massive objects which are already closer to the mass resolution limit at z = 0.

Zoom In Zoom Out Reset image size

A.1. Fitting Formulae

Different fitting formulae for the MAH shown in Figure 7 exist in the literature. We focus on two of them. By using the extended Press-Schechter (EPS) formalism (Bond et al. 1991; Lacey & Cole 1993), van den Bosch (2002) proposed

$$\begin{eqnarray}&&M(z)={M}_{0}\mathrm{exp}\left\{\mathrm{ln}(1/2){\left[\displaystyle \frac{\mathrm{ln}(1+z)}{\mathrm{ln}(1+{z}_{f})}\right]}^{\nu }\right\},\end{eqnarray} \tag{ 4 }$$

where z_f and ν are free parameters. The redshift formation z_f indicates the redshift when the halo has accreted half of its final mass. We fit both the median and the mean MAH with the equation above by assuming Poisson errors weighted by the number of halos in each redshift bin. We quantify the deviation from this analytic description of the MAH with the rms of the fit

$$\begin{eqnarray}&&{(\mathrm{rms})}^{2}=\displaystyle \frac{1}{N}\displaystyle \sum _{N}\displaystyle \frac{{[M({z}_{i})-{M}_{{\rm{model}}}({z}_{i})]}^{2}}{M{({z}_{i})}^{2}}.\end{eqnarray} \tag{ 5 }$$

We list the best-fit parameters of Equation (4) along with the rms of the fits in Table 4. As expected in hierarchical clustering scenarios, the value of the best-fit parameter z_f increases with decreasing mass because more massive objects tend to form later than less massive ones.

Table 4.  Best-fit Parameters of the Formulae by Van Den Bosch, Equation (4), and McBride, Equation (6)

	Median				Mean
	van den Bosch
${M}_{200}[{10}^{14}{M}_{\odot }\;{h}^{-1}]$	zf	ν		rms	zf	ν		rms
10	0.358 ± 0.011	1.227 ± 0.035		0.044	0.381 ± 0.012	1.252 ± 0.037		0.043
1	0.711 ± 0.003	1.702 ± 0.012		0.006	0.709 ± 0.003	1.582 ± 0.012		0.022
0.1	0.944 ± 0.005	1.675 ± 0.019		0.025	0.925 ± 0.005	1.519 ± 0.018		0.015
0.05	1.033 ± 0.010	1.317 ± 0.024		0.032	1.036 ± 0.011	1.242 ± 0.022		0.024
	McBride
${M}_{200}[{10}^{14}{M}_{\odot }\;{h}^{-1}]$	β	γ	$\beta -\gamma$	rms	β	γ	$\beta -\gamma$	rms
10	−1.000 ± 0.270	1.150 ± 0.180	−2.160	0.048	−0.690 ± 0.270	1.280 ± 0.180	−1.970	0.028
1	+1.283 ± 0.044	1.920 ± 0.029	−0.637	0.030	+0.812 ± 0.043	1.566 ± 0.028	−0.754	0.010
0.1	+0.898 ± 0.051	1.342 ± 0.035	−0.444	0.039	+0.480 ± 0.051	1.065 ± 0.035	−0.585	0.029
0.05	−0.005 ± 0.070	0.658 ± 0.051	−0.664	0.053	−0.217 ± 0.069	0.513 ± 0.050	−0.730	0.037

Download table as:  ASCIITypeset image

The second formula we considered was first proposed by Tasitsiomi et al. (2004) and widely studied by McBride et al. (2009):

$$\begin{eqnarray}&&M{(z)={M}_{0}(1+z)}^{\beta }{e}^{-\gamma z},\end{eqnarray} \tag{ 6 }$$

where β and $\gamma (\geqslant 0)$ are free parameters. This formula represents an exponential growth with a redshift-dependent correction. It is a revision of the simple one-parameter exponential form $M(z)={M}_{0}{e}^{-\alpha z}$ (Wechsler et al. 2002), where $\alpha =\mathrm{ln}(2)/{z}_{f}$. By using the EPS formalism, Correa et al. (2015a) showed that in a ΛCDM model the exponential growth is a good description of the MAH at high z, while the power-law behavior at low z is necessary because the accelerated expansion of the universe slows down the accretion. For this reason Equation (6) appears to be a general description of the MAH of dark matter halos in a ΛCDM model independently of the halo mass. The parameters β and γ are related to the power spectrum (Correa et al. 2015a). The value of $\beta -\gamma$ is a parameter describing the mass growth rate at low redshift. We use Equation (6) to perform a fit with Poisson errors weighted by the number of halos in each redshift bin and evaluate the rms as in Equation (5). We list the best-fit parameters of Equation (6) and the rms of the fits in Table 4.

A.2. Comparison with a Theoretical Model

By following and generalizing the formalism of Lacey & Cole (1993) and Nusser & Sheth (1999), Giocoli et al. (2012) introduced a new theoretical model to describe the MAH of dark matter halos. This simple model, which enables the derivation of a generalized redshift formation distribution, has already been applied to the CoDECS simulations in Giocoli et al. (2013).

Here we summarize the relevant definitions and refer to Giocoli et al. (2012) for further details. The model defines the redshift formation z_f of a halo of a given mass M₀ at a given redshift z₀ as the redshift when the object has accreted a fraction f of its final mass M₀, for any fraction 0 < f < 1. The variance of the linear fluctuation field is

$$\begin{eqnarray}&&S(M)=\displaystyle \frac{1}{2{\pi }^{2}}{\displaystyle \int }_{0}^{\infty }{W}^{2}({kR}){P}_{{\rm{lin}}}(k){k}^{2}{dk},\end{eqnarray} \tag{ 7 }$$

where W(kR) is a top-hat window function of scale $R={(3M/4\pi {\rho }_{m})}^{1/3}$, ${\rho }_{m}$ is the comoving background density of the universe, and ${P}_{{\rm{lin}}}(k)$ is the linear power spectrum.

The initial threshold overdensity for spherical collapse is

$$\begin{eqnarray}&&{\delta }_{c}(z)=\displaystyle \frac{{\delta }_{c,\mathrm{DM}}}{{D}_{+}(z)},\end{eqnarray} \tag{ 8 }$$

where ${\delta }_{c,\mathrm{DM}}$ is the linear overdensity at redshift z and ${D}_{+}(z)$ is the growth factor normalized to unity at z = 0.

The cumulative formation redshift distribution, in terms of the scaled variable

$$\begin{eqnarray}&&{w}_{f}=\displaystyle \frac{{\delta }_{c}({z}_{f})-{\delta }_{c}({z}_{0})}{\sqrt{S({fM})-S(M)}},\end{eqnarray} \tag{ 9 }$$

is given by

$$\begin{eqnarray}&&P(\gt {w}_{f})=\displaystyle \frac{{\alpha }_{f}}{{e}^{{w}_{f}^{2}/2}+{\alpha }_{f}-1}.\end{eqnarray} \tag{ 10 }$$

The model has a single free parameter ${\alpha }_{f}$ which depends on the fraction f of the final mass that is used to define the formation redshift z_f. For the same set of simulations we use here, Giocoli et al. (2013) find

$$\begin{eqnarray}&&{\alpha }_{f}=\displaystyle \frac{1.365}{{f}^{0.65}}{e}^{-2{f}^{3}}.\end{eqnarray} \tag{ 11 }$$

Since Equation (10) can be inverted, it is possible to evaluate the median redshift z_f when a halo accretes a fraction f of its final mass M with the relation

$$\begin{eqnarray}&&{\delta }_{c}({z}_{f})={\delta }_{c}({z}_{0})+\tilde{{w}_{f}}\sqrt{S({fM})-S(M)},\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&\tilde{{w}_{f}}=\sqrt{2\mathrm{ln}({\alpha }_{f}+1)}\end{eqnarray} \tag{ 13 }$$

is the median value of w_f defined by the usual relation $P(\gt {\tilde{w}}_{f})=1/2$. Equation (12) can be translated into a MAH for a given final mass M₀.

We compare the simulation results and the Giocoli model in Figure 7. The model of Giocoli et al. (2012) is built by evaluating the median redshift at which the halo has accreted a fixed fraction of its final mass whereas we evaluate the mean and the median M₂₀₀ for all the objects in a given mass bin at each redshift. We list the rms defined in Equation (5) in Table 5. In general, the global agreement between the theoretical model of Giocoli et al. (2012) and the simulation in each mass bin is similar for both the mean and the median, as it can be seen from the rms values in Table 5. In the largest-mass bin the model overestimates the MAH obtained from the simulation, while in the other three bins the agreement is within a few percent up to z ∼ 1. Toward higher z, the model starts underestimating the simulation MAH. This discrepancy is more pronounced and appears at decreasing redshifts for decreasing halo mass. This behavior originates from the mass resolution and the consequent decrease of the number of halos ${N}_{{\rm{hal}}}$ at a given redshift.

Table 5.  MAH rms for the Giocoli Model

	Median	Mean
${M}_{200}[{10}^{14}{M}_{\odot }\;{h}^{-1}]$	rms
10	0.144	0.136
1	0.058	0.052
0.1	0.059	0.058
0.05	0.101	0.103

Download table as:  ASCIITypeset image