Galaxy cluster mass accretion rates from IllustrisTNG

We use simulated cluster member galaxies from Illustris TNG300-1 to develop a technique for measuring the galaxy cluster mass accretion rate (MAR) that can be applied directly to observations. We analyze 1318 IllustrisTNG clusters of galaxies with $M_{200c}>10^{14}$M$_\odot$ and $0.01\leq z \leq 1.04$. The MAR we derive is the ratio between the mass of a spherical shell located in the infall region and the time for the infalling shell to accrete onto the virialized region of the cluster. At fixed redshift, an $\sim 1$ order of magnitude increase in $M_{200c}$ results in a comparable increase in MAR. At fixed mass, the MAR increases by a factor of $\sim 5$ from $z=0.01$ to $z=1.04$. The MAR estimates derived from the caustic technique are unbiased and lie within 20% of the MARs based on the true mass profiles. This agreement is crucial for observational derivation of the MAR. The IllustrisTNG results are also consistent with (i) previous merger tree approaches based on N-body dark matter only simulations and with (ii) previously determined MARs of real clusters based on the caustic method. Future spectroscopic and photometric surveys will provide MARs of enormous cluster samples with mass profiles derived from both spectroscopy and weak lensing. Combined with future larger volume hydrodynamical simulations that extend to higher redshift, the MAR promises important insights into evolution of massive systems of galaxies.

Many theoretical studies investigate the MAR in ΛCDM cosmologies.The first approaches to computing the MAR in ΛCDM cosmologies were analytic models based on the extended Press-Schechter (EPS) formalism or on Monte-Carlo generated merger trees (Bower 1991;Lacey & Cole 1993).More recent approaches build merger trees from N-body dark matter only simulations (McBride et al. 2009;Fakhouri et al. 2010;Diemer et al. 2017b;Diemer 2018) or from semi-analytical models calibrated with N-body simulations (van den Bosch et al. 2014;Correa et al. 2015b).These models imply that halos with mass ≳ 10 14 M ⊙ accrete ∼ 30% − 50% of their mass over the redshift range z ∼ 0.5 to z ∼ 0. The MAR is correlated with both halo mass and redshift.Because we observe a cluster at a single redshift and thus cannot measure its history directly, the merger tree approaches are not directly applicable to cluster observations.
Soon, spectroscopic and photometric missions will provide huge samples of clusters with both dense spectroscopy and weak-lensing maps extending to large cluster-centric radius.The cluster samples will also cover redshifts ≲ 2. Both weak-lensing (Bartelmann 2010;Hoekstra et al. 2013;Umetsu 2020) and the caustic technique (Diaferio & Geller 1997;Diaferio 1999;Serra et al. 2011;Pizzardo et al. 2023) will allow estimates of cluster mass profiles without assuming dynamical equilibrium.
The use of IllustrisTNG (Pillepich et al. 2018;Springel et al. 2018;Nelson et al. 2019) enables direct application of the caustic technique to track the MAR for galaxy clusters with z ≲ 1.In contrast with a MAR recipe based on purely N-body simulations (Pizzardo et al. 2021), IllustrisTNG enables the use of galaxies as tracers of the dynamical evolution of galaxy clusters.Following Pizzardo et al. (2021) we define the MAR as the ratio between the mass of an infalling shell and the infall time.Pizzardo et al. (2021) and Pizzardo et al. (2022) compute the MAR based on mass profiles determined from the caustic technique applied to dense spectroscopic surveys of the infall regions of observed galaxy clusters (Rines & Diaferio 2006;Rines et al. 2013;Sohn et al. 2021a,b).The results are consistent with ΛCDM predictions.Pizzardo et al. (2023) use IllustrisTNG to calibrate a statistical platform for application of the caustic technique.Their approach is based on mock galaxy cluster member catalogs.The caustic technique calibrated by IllustrisTNG returns the true cluster mass profile within 10% in the radial range (0.6 − 4.2)R 200c and redshift range 0.01 − 1.04.
We build on the approach of Pizzardo et al. (2023) to develop an IllustrisTNG recipe for computing the cluster MAR.The recipe is based on simulated galaxies rather than dark matter particles in contrast with previous approaches.We demonstrate that the caustic mass profiles based on galaxy mock catalogs (Pizzardo et al. 2023) yield reliable estimates of the true MARs.They are also consistent with estimates from previous theoretical and observational investigations (McBride et al. 2009;Fakhouri et al. 2010;van den Bosch et al. 2014;Correa et al. 2015a,b;Diemer et al. 2017b;Diemer 2018;Pizzardo et al. 2021Pizzardo et al. , 2022)).
Sect. 2 describes the IllustrisTNG cluster sample.We derive the average radial velocity profiles, fundamental ingredients for the computation of the MAR.Sect. 3 describes the approach to estimating the MAR.In Sect. 4 we derive the mass of the infalling shell and the infall time.Sect. 5 compares the true MARs with the caustic MARs.Sect.6 compares these MAR results with previous models and observations.We also compare galaxy and dark matter MARs.Finally we outline future challenges and prospects.We conclude in Sect.7.

Cluster sample, mass and velocity profiles
Basic inputs to the mass accretion rate include the cluster mass and radial velocity profiles.A large sample of clusters (Pizzardo et al. 2023) from the IllustrisTNG simulations (Pillepich et al. 2018;Springel et al. 2018;Nelson et al. 2019) is the basis for the determination of the true 3D and projected caustic cumulative mass profiles.As a basis for computing the mass accretion rate, we compute the average cluster radial velocity profile at each redshift.This profile has a characteristic minimum that ultimately determines the accretion rate.
Sect. 2.1 describes the sample of clusters extracted from Il-lustrisTNG by Pizzardo et al. (2023).We summarize the derivation of the true and caustic mass profiles for each cluster.Sect.2.2 describes the determination of the radial velocity profiles.et al. (2023) extract cluster samples from the TNG300-1 run of the IllustrisTNG simulations (Pillepich et al. 2018;Springel et al. 2018;Nelson et al. 2019), a set of gravomagnetohydrodynamical simulations based on the ΛCDM model.Table 1 lists the cosmological parameters of the simulations.TNG300-1 is the baryonic run with the highest resolution among the runs with the largest simulated volumes.The simulation has a comoving box size of 302.6 Mpc.TNG300-1 contains 2500 3 dark matter particles with mass m DM = 5.88 × 10 7 M ⊙ and the same number of gas cells with average mass m b = 1.10 × 10 7 M ⊙ .0.01, 0.11, 0.21, 0.31, 0.42, 0.52, 0.62, 0.73, 0.82, 0.92, and 1.04.For ∼ 22% of the clusters sparse sampling, rich foreground or background, or the presence of significant substructures preclude robust application of the caustic technique (Pizzardo et al. 2023).We remove these clusters from the sample.

Pizzardo
Table 2 describes the remaining 1318 clusters we analyze here.The Table includes the number of clusters in each redshift bin, the median and the interquartile range of their masses M 3D 200c , and the minimum and maximum M 3D 200c at each redshift.Our goal is assessment of the caustic technique (Diaferio & Geller 1997;Diaferio 1999;Serra et al. 2011;Pizzardo et al. 2023) as a basis for reliable estimates of the true MAR.Estimation of the MAR requires robust knowledge of the cluster mass at large cluster-centric distances, ≳ 2R 200c , where virial equilibrium does not hold (Sects. 1 and 3).In the extended radial range (0.6 − 4.2)R 200 the caustic technique returns an unbiased estimate of the mass with better than 10% accuracy and with a relative uncertainty of 23% provided that the velocity field of the cluster outer region is sufficiently well sampled (Pizzardo et al. 2023).The caustic technique is independent of equilibrium assumptions.
The true 3D and caustic MARs rest on determination of the true and caustic mass profiles for each cluster in the sample (Table 2).Pizzardo et al. (2023) compute the true cumulative mass profile (from now on, the "true mass profile") for each cluster from the 3D distribution of matter extracted from raw snapshots.These profiles include all matter species: dark matter, gas, stars, and black holes.Pizzardo et al. (2023) compute the true mass profile for each cluster, M 3D (r), in 200 logarithmically spaced bins covering the radial range (0.1 − 10)R 3D 200c .These profiles define R 3D 200c and M 3D 200c for each cluster.The r − v los diagram, the line-of-sight velocity relative to the cluster median redshift as a function of r, v esc,los (r) (Diaferio & Geller 1997;Diaferio 1999), is the basis for estimating the caustic cumulative mass profile (caustic mass profile hereafter).In the r − v los diagram, cluster members appear within a welldefined trumpet-shaped pattern.The amplitude of the pattern decreases as r increases.The caustic technique locates the boundaries (the "caustics") of this pattern that discriminate between cluster members and the foreground/background.The caustic amplitude, A(r) is half of the vertical separation between the upper and the lower caustic at radius r.Diaferio & Geller (1997) show that the caustic amplitude approximates the escape velocity from the cluster, A(r) ≈ v esc,los (r).The square of the caustic amplitude, ⟨v 2 esc,los (r)⟩, is then a measure of the gravitational potential.
The caustic technique estimate of the mass profile is where G is the gravitational constant and F β is a constant filling factor.Pizzardo et al. (2023) use the Illustris TNG simulations to calibrate the filling factor in the redshift range 0.01 − 1.04.We use their calibration, F β = 0.41 ± 0.08, throughout.The r −v los diagram is based on catalogues of simulated cluster galaxies that include the right ascension α, declination δ, and total redshift z along the line of sight.Pizzardo et al. (2023) associate a realistic galaxy mock redshift survey with each simulated cluster.These catalogues include contaminating background and foreground galaxies.Pizzardo et al. (2023) build the catalogues by identifying galaxies with Subfind substructures that have stellar mass > 10 8 M ⊙ .This choice of mass limit mimics observable galaxies and it assures optimal performance of the caustic technique.Finally, Pizzardo et al. (2023) apply the caustic technique in the same unconstrained way to all of the mock catalogues and obtain a single calibrated caustic mass profile for each cluster.These profiles define R C 200c and M C 200c for each cluster.

Radial velocity profiles
The average cluster galaxy velocity profile along the radial direction is fundamental to the evaluation of the mass accretion rate (Sect.3).We compute the set of individual cluster radial velocity profiles based on the comoving position of simulated galaxies with respect to the cluster center, r c,i , and the galaxy peculiar velocity, v p,i .From the 3D volume we select only galaxies with cluster-centric distances < 10R 3D 200c .We compute the radial velocity of each galaxy: , where H(z s ) and a(z s ) are the Hubble function and the scale factor at the redshift z s of the snapshot.We compute the mean radial velocity profile based on the galaxies within 100 linearly spaced radial bins covering the range (0, 10)R 3D 200c .For each redshift snapshot, we compute a single average radial velocity profile by averaging over all of the individual radial velocity profiles for clusters within the snapshot (Table 2).For each of the 100 radial bins of r/R 3D 200c , we compute the mean and the standard deviation of all of the velocities within all of the clusters velocities in each bin.For each redshift, we thus obtain from all the individual clusters a single average mean radial velocity profile with the dispersion around it.
The solid blue curves in Fig. 1 show mean radial velocity profiles at three example redshifts: z = 0.01, z = 0.62, and z = 1.04 (in the left, middle, and right panels, respectively).We smooth over statistical fluctuations in the profile by applying a Savitsky-Golay filter with a 10 radial bin window (Savitzky & Golay 1964).The dash-dotted curves show the resulting smoothed profiles.The shaded light blue band shows the error in the corresponding smoothed profiles.
The curves show three clear regimes.Within ∼ 1R 3D 200c , the average radial velocity has a plateau at zero or small velocity.This region is virialized, with no net infall.At larger distances from the cluster center, ∼ (1 − 4)R 3D 200c , the infall region, the radial velocity is definitely non-zero and directed toward the cluster center.Finally, at distances ≳ 4R 200 the galaxy velocity is directed out of the cluster because the galaxies are still coupled to the universal Hubble flow.
Estimation of the MAR according to Eq. ( 3) requires determination of the minimum of the radial velocity profile (see Sects. 3 and 4).At each redshift, we measure the minimum velocity of the smoothed average radial profile.The first three columns of Table 3 list the minimum radial velocity, v min , and its clustercentric location, R v min , for each redshift bin.To compute the error in R v min , we bootstrap 1000 samples of clusters at each redshift (Table 2, columns 1 and 2).
Table 2 shows that the median cluster mass generally decreases as the redshift increases because very massive systems are progressively less abundant at greater redshift.We check the impact of the changing distribution of cluster masses on R v min .For each redshift, we build homogeneous samples that include only clusters with mass M 3D 200c in the range (1.02−4.30)•10 14M ⊙ .The largest minimum and smallest maximum mass for the samples in the 11 redshift bins set this range (see last column of Table 2).According to the Kolmogorov-Smirnov test, these clipped samples share indistinguishable mass distributions; the p-values are in the range (0.4 − 0.9).For each redshift, we compute the average radial velocity profile for the clipped sample and locate the turnaround radius.The turnaround radii for these samples are within ≲ 1.7% of the turnaround radii of the full samples.They are also unbiased.Hence the R v min radii are insensitive to the difference among the distribution of cluster masses at different redshifts.
Fig. 2 shows v min as a function of its R v min , colour coded by redshift.The minimum velocity and its cluster-centric radius are strongly correlated with redshift.From z = 0.01 to z = 1.04 the v min of clusters with approximately equal mass (Table 2) increases by ∼ 100% in absolute value.The corresponding R v min decreases by ∼ 40%.
The behavior of v min as a function of redshift is in qualitative agreement with the general picture of hierarchical structure formation.For equal mass clusters, systems at higher redshift form within higher overdensity peaks.Because of their denser environment, these clusters are denser within ∼ R 200c .Thus more mass is confined within smaller cluster-centric radii than for an equally massive lower redshift cluster.The consequently deeper gravitational potential within ∼ R 200c at greater redshift is the source of the larger minimum infall velocity located at smaller cluster-centric radius.The corresponding mass accretion is also larger in these higher redshift, denser environments.Nbody simulations (e.g.van den Bosch 2002;McBride et al. 2009;Fakhouri et al. 2010;van den Bosch et al. 2014;Diemer et   2017b) show that at fixed mass, the accretion rate at z ∼ 1 exceeds the rate at z ∼ 0 by roughly an order of magnitude.
Table 3. Radial velocity minima and radii as a function of redshift.The dynamics within the infall region is independent of the turnaround radius where galaxies depart from the Hubble flow as a result of the cluster potential (Gunn & Gott 1972;Silk 1974;Schechter 1980).We identify the turnaround radius from the smoothed radial velocity profiles.Starting from the largest cluster-centric distances, the turnaround radius is the first intersection between the profile and the axis v rad = 0.The vertical solid lines in Fig. 1 show that the turnaround radii agree well with Meiksin (1985).The turnaround radii are in the range (4.52 In contrast with R v min , the turnaround radius decreases by only ∼ 3.3% as the redshift increases from z = 0.01 to z = 1.04.

Recipe for the estimation of the MAR
The average radial velocity profiles (Sect.2.2, Fig. 1) have clear minima at cluster-centric radii ∼ (2 − 3)R 200c identifying the region where clusters accrete new material.For an ideal, theoretical MAR (hereafter MAR t ), we first calculate the rate at which material falls through an infinitesimal shell at this minimum velocity.However, real and simulated systems have too few galaxies in this shell.Thus, we derive a second MAR in a larger volume that approximates MAR t .We analyze the Illustris TNG300-1 simulations to relate the observational MAR to MAR t .
Assuming spherical symmetry, the mass of an infinitesimal shell at a distance r from the cluster center is dM shell = 4πρ(r)r 2 dr, where dr is the width of the shell and ρ is the shell density.The time derivative of the shell mass is Ṁ = dM shell /dt, where dt = dr/v rad (r) is the time for infalling material with velocity v rad (r) to cross dr.We define a theoretical mass accretion rate MAR t = Ṁ(R v min ), where R v min is the minimum in the radial velocity profile (Sect.2.2, Fig. 1, and Table 3): IllustrisTNG provides an estimate of MAR t for a thin, finite shell where the radial velocity ≈ v min is nearly constant acrosss the shell.
Observations can only provide the mass of an infalling shell with a finite width, ∆R ≫ dr.We define the observational MAR of clusters of galaxies as the accretion of matter within a spherical shell onto the virialized region of the cluster: where M shell is the mass of an infalling shell with width ∆R, t inf is the time for the shell to accrete onto R 200c , and K is a scaling factor.In an infinitesimal shell (∆R = dr) located at R v min , the shell mass is 3) and (2) are equivalent.Below, we derive K for thick shells of finite width.
Deriving the MAR from Eq. (3) depends on estimates of the mass in the infalling shell and the infall time.The caustic method provides the mass of the infalling shell from observations (Diaferio & Geller 1997;Diaferio 1999;Serra et al. 2011;Pizzardo et al. 2023).To calibrate this technique, we compute the true mass of the shell by summing the mass of dark matter, gas, stars, and black holes in a shell of width ∆R.
There is no direct observational route to the infall time because v in f is currently inaccessible to observation.Although t inf depends on redshift, it is fairly independent of cluster mass (see Fig. 7 below).Thus, the combination of M shell from observations and t inf derived from simulations yield a MAR based on cluster observations.
To compute the infall time, we solve the equation of radial infall with a nonconstant acceleration derived from the true gravitational potential of the cluster.The parameters of the equation are the initial infall velocity of the infalling shell, v inf , the center of the infalling shell equivalent to the radial location of v inf , R v min , and the radius that defines the virialized region of clusters, R 200c (Sect.4.2).
We locate the infalling shell (Sect.4.1) from the simulated radial velocity profile of the cluster (Sect.2.2).The boundaries of the infalling shell are R shell,i (inner cluster-centric radius), and R shell,o (outer cluster-centric radius), where the smoothed radial velocity is a fraction A of v min .
To link MARs in Eq. 3 with MAR t , we compute MARs with a range of A = 0.45-0.95.Smaller A overlaps the virialized regions of clusters.Larger values of A approach the ideal, theoretical limit in Eq. (2).However, MAR estimates for these thin shells are less robust because the number of simulated particles in the infalling shell is small.We compare MAR C and MAR 3D , the MARs estimated according to Eq. ( 3) respectively using M C shell and M 3D shell , the caustic and 3D shell masses, for K = 1.Averaging over redshift, the median ratio between MAR 3D and MAR C , r 3c , increases monotonically from r 3c = 0.92 ± 0.05 (A = 0.45) to r 3c = 1.09 ± 0.04 (A = 0.95).The error in r 3c is the interquartile range of the 11 ratios at the redshifts we consider.We adopt A = 0.72, where r 3c = 1.00 ± 0.04, as the best value.The choice of A = 0.72 provides a large enough galaxy sample for application of the caustic technique.The shaded blue vertical regions in Fig. 1 show shells with A = 0.72 for three redshifts: z = 0.01, z = 0.62, and z = 1.04 (left to right, respectively).
Finally, we connect MARs estimated in extended shells to MAR t and estimate the appropriate value for K. Colored lines in Fig. 3 show the ratio MAR 3D (A)/MAR 3D (A = 0.72) as a function of A at six different redshifts from z = 0.01 to z = 1.04 for K = 1.The error bar shows the interquartile range of the ratios at A = 0.72.The estimated MARs decrease as A increases.
When A = 0.95, the radial velocity is nearly constant across the shell and is (0.97-0.98)v min .This thin shell provides an approximation to MAR t .Averaging over the 11 redshifts, the ratio MAR 3D (A = 0.95)/MAR 3D (A = 0.72) is K = 0.35.This value provides a plausible scaling between MAR t and the M shell s and MARs estimated from cluster observations with A = 0.72.

M shell and t inf from 3D and caustic mass profiles
Computation of the MAR of clusters requires an estimate of the mass of the infalling shell, M shell , and of the infall time for the cluster to accrete the shell, t inf .Section 4.1 describes the determination of M shell .Section 4.2 describes the derivation of the infall time, t inf .

The infalling shell and its mass M shell
The average radial velocity profiles (Sect.2.2; Fig. 1) are the basis for determining the cluster-centric radius of the infalling shell.The radial velocity profiles in Sect.2.2 show a clear infall pattern at radii ∼ (2 − 3)R 200c where the minimum radial velocity occurs.Table 3 lists the average minimum radial velocity, v min , and its radial location, R v min , at each redshift.
Based on the discussion of Section 3 we identify the inner (R shell,i ) and outer (R shell,o ) boundaries of the infalling shell where the smoothed radial velocity is 0.72v min (A = 0.72; shaded blue vertical regions in Fig. 1).Table 3 (fifth column) shows the radial location of the infalling shell at each redshift.The clustercentric distance of the shells decreases as redshift increases: at low redshifts, z ≲ 0.4, the shells are located in the range ∼ (1.6 − 3.4)R 200c ; at high redshift, z ≳ 0.73, the shells are in the range ∼ (1.3 − 2.6)R 200c .This variation results from the decrease of R v min with increasing redshift (Sec.2.2, see Fig. 2).The width of the infalling shell, ∼ (1.2 − 1.3)R 200c , is nearly redshift independent.
For each simulated cluster, we compute the true mass, M 3D shell (A = 0.72), and caustic observable mass, M C shell (A = 0.72), of the infalling shell.For the appropriate radial range (Table 3) we compute two true and caustic masses that measure M(< R shell,i ; A = 0.72), the total mass within the inner boundary of the shell, and M(< R shell,o ; A = 0.72), the total mass within the outer boundary of the shell.The observable mass of the infalling shell is then the difference in these masses M shell (A = 0.72) = M(< R shell,o ; A = 0.72) − M(< R shell,i ; A = 0.72).
Fig. 4 shows the true M 3D shell = K M 3D shell (A = 0.72) as a function of the true M 3D 200c at three different redshifts: z = 0.01, z = 0.62, and z = 1.04 (cyan, violet, and magenta, respectively).The squares with error bars show the medians and interquartile ranges of the simulated data in eight logarithmic mass bins cov- ering the range (1 − 12.6) • 10 14 M ⊙ .The lines show a power law fit to the data.M 3D shell and M 3D 200c are strongly correlated at every redshift: a change of one order of magnitude in mass leads to a comparable change in the mass of the infalling shell.Kendall's test gives a correlation index of ∼ 0.44 with associated p-values in the range ∼ 10 −5 − 10 −28 .This correlation is expected in the hierarchical cluster formation model, because at fixed redshift higher mass clusters reside within greater overdensities and thus there is more mass in their accreting shells.Conversely, M shell is not strongly correlated with redshift.Although the fits show that M 3D shell generally increases by ∼ 20 − 60% from z = 0.01 to z = 1.04, the error bars show that this correlation is not statistically significant.
Coloured points in Fig. 5   The dependence of M shell and redshift is expected in the standard model of structure formation and evolution.In Appendix A we explain how a complex interplay between the density and mass profiles as a function of redshift along with the slightly different cluster mass distribution at various redshifts tends to suppress any correlation between M shell and redshift.
The caustic technique provides robust estimates of the mass profile of clusters in the infalling shell (Pizzardo et al. 2023).Points with error bars in Fig. 6 show the median and the interquartile range of individual ratios between the caustic and true mass M C shell /M 3D shell as a function of redshift.Analogously to M 3D shell , we define M C shell = K M C shell (A = 0.72).Because the infalling region is closer to the cluster center as the redshift increases (Table 3, fifth column), the median ratio slowly rises from 0.9 at z = 0 to 1.1 at z = 1.04. 1 The typical dispersion in the ratio is ∼ 34%.Thus the caustic technique slightly overestimates or underestimates the mass at smaller and larger cluster-centric distances, respectively.The trend is not statistically significant: Kendall's τ correlation coefficient is ≲ 0.05.The caustic mass M C shell is an unbiased estimate of M 3D shell .Furthermore M C shell has the same dependence on mass, redshift, and R v min as M 3D shell on average.

The infall time t inf
Based on the full 3D cluster properties at each redshift, we compute t inf , the time for an infalling shell to accrete onto the virialized region of the cluster.We derive t inf by tracking the radial infall with nonconstant acceleration driven by the cluster gravitational potential.The locations of the minima in the average radial velocity profiles (Sect.2.2) set the limit of the radial range where we compute the infall time.We ultimately employ the infall times based only on the true full 3D data because the un-1 (Pizzardo et al. 2023) show that on average in the radial range (0.6 − 4.2)R 200c the caustic mass is within the 10% of the true mass.However, the caustic to true mass ratio slowly decreases from the lower to the upper end of the calibrated range.certainty in these measures derived by application of the caustic method are too large.
We compute the time for the center of the infalling shell to reach R 200c as a measure of t inf .We measure the infall time starting from the radius of the minimum infall velocity, R v min (Figure 1 and Table 3).In other words, t inf is the time required for the center of the infalling shell initially located at R v min to reach R 200c .
Over the radial range R v min to R 200c , the gravitational acceleration changes substantially.To compute t inf of each cluster in each redshift bin we use an iterative procedure.
We divide the radial range (R 200c , R v min ) into N +1 = 101 bins.The radial step between two contiguous bins is ∆r At each step n < N, we calculate ∆t n , the time for the center of the shell to move from r n = R v min − n∆r to r n+1 = R v min − (n + 1)∆r.Within this small radial interval with ∆r ∼ 0.01R 200c ≈ 0.01Mpc the gravitational acceleration is nearly constant.From this acceleration and the initial infall velocity of the shell, a n and v n respectively, we compute ∆t n as the positive (physical) solution of the equation that is2 We obtain the acceleration of the infalling shell by computing the cluster gravitational potential ϕ(r) based on the true cluster shell density profile ρ(r).The Poisson equation for an isolated spherical system with shell density profile ρ I (r) is: The uniform cosmological background density, ⟨ρ(z)⟩ = Ω M (z)ρ c (z), exerts no net gravitational effect.Thus we compute the gravitational potential based on the mass density fluctuations by replacing ρ I (r) with ρ(r) − ⟨ρ⟩.The second integral is finite; at sufficiently large cluster-centric distances, ∼ 10R 200c (see Appendix A), the correlated cluster density is ∼ 0. We replace the upper limit of the second integral of Eq. ( 6) with 10R 200c .From the true cluster potential ϕ(r) we compute the gravitational acceleration induced by the cluster, a(r), by simple differentiation: a(r) = − d dr ϕ(r).At each step n, we set a n (Eq.( 4) to the value of a(r) at the position of the center of the infalling shell, a n = a(R v min − n∆r).
At the first step n = 0, the initial velocity is v 0 = v inf , where v inf is the average cluster radial velocity of the infalling shell at that redshift (fourth column of Table 3).At each succeeding step n we increment the initial infall velocity of the shell by taking a constant acceleration in the small time step: Application of Eq. ( 5) from n = 0 to n = N − 1 yields a set of N time steps, ∆t n={0,...,N−1} .The sum of these time steps is our estimate of the cluster infall time: Figure 7 shows the resulting infall times for individual clusters as a function of M 3D 200c in three redshift bins: z = 0.01, 0.62, and 1.04 (cyan, violet, and magenta, respectively).The squares with error bars show the medians and interquartile ranges of the simulated data in eight logarithmic mass bins covering the range (1 − 12.6) • 10 14 M ⊙ (as in Fig. 4).The lines show a power law fit.
The infall t inf is correlated with redshift.The increased cluster radial acceleration resulting from the larger cluster density at high redshift produces this correlation (Eq.6).Equation 5shows that the increased acceleration produces a corresponding decrease in the infall time, t inf .
The squares in Fig. 7 show the absence of correlation between t inf and the cluster mass M 3D 200c at each redshift.The higher density of more massive clusters generates a larger acceleration decreasing t inf .However, more massive clusters are also more extended thus increasing the radial range (R 200c , R v min ) and correspondingly increasing t inf .These effects result in a minimal dependence of t inf on M 3D 200c .

The Mass Accretion Rate
We apply Eq. 3 to estimate the MARs at different redshifts.We compute MAR 3D , the true MAR based on M 3D shell and t inf , and MAR C , the caustic MAR based on M C shell .As described earlier, we use 3D data to model the infall time.We compare MAR C with MAR 3D to assess the caustic technique as a robust method for estimating the true MAR.
According to Eq. 3, the MAR 3D of each individual cluster is the ratio between its respective M 3D shell (see Sect. 4.1) and t inf (see Sect. 4.2).Fig. 8 shows the true MARs of individual clusters as a function of M 3D 200c in three different redshift bins: z = 0.01, 0.62, and 1.04 (cyan, violet, and magenta, respectively).The squares with error bars show the medians and interquartile ranges of the simulated data in eight logarithmic mass bins covering the range (1 − 12.6) • 10 14 M ⊙ .The lines show a power law fit.
The MARs increase both with increasing mass and with increasing redshift.A correlation between MAR and M 200c is expected because more massive clusters tend to be surrounded by larger amounts of mass. Figure 8 shows that at fixed redshift  and 4 show that the increase of MAR 3D with M 3D 200c is consistent with the increase of M 3D shell with M 3D 200c .Fig. 7 shows that the infall time does not play a significant role in this correlation.
In the hierarchical clustering paradigm, clusters of fixed mass at higher redshift reside within denser regions and thus accrete faster than clusters of the same mass at lower redshift.Figure 8 shows that at fixed mass MAR 3D increases by a factor ∼ 2.2 from z = 0.01 to z = 0.62, and by a factor ∼ 5 from z = 0.01 to z = 1.04.This effect originates from the anticorrelation between t inf and redshift (Sect.4.2): Figs. 8 and 7 show that the increase of MAR 3D with redshift is consistent with this decrease of t inf with redshift.
We fit the individual MAR 3D s and M 3D 200c s at each redshift to the relation MAR = a M 3D 200c /10 14 M ⊙ b (Table 4).The fits show  We fit the analytic relation proposed by Fakhouri et al. (2010) to the MAR as a function of redshift: Because of the limited mass range we sample at greater redshifts, we fix V = b and we fit the median MAR 3D as a function of the median M 3D 200c and the redshift.The resulting coefficients are: U = 128 ± 13, V = b = 0.90 ± 0.18, and W = 1.69 ± 0.34.We compare these results with earlier work in Sect.6.
The MAR C of individual clusters is the ratio between M C shell and the infall time t fit inf (M C 200c ) (Sect.4.2), where we base the computation on the cluster caustic mass M C 200c .We compute the cluster infall time t fit inf (M C 200c ) by evaluating the power law fit of the individual infall times at the cluster caustic mass M C 200c as a function of mass derived by true profiles at the given redshift (Fig. 7).The upper panel of Fig. 9 shows the median and interquartile range of the true (blue squares and solid error bars) and caustic (red triangles and dotted error bars) MARs as a function of redshift.The black curve is a fit of Eq. ( 8) to the MARs as a function of redshift.We fix the mass to the median mass of the cluster sample.The residuals around the fit are, on average, only ∼ −0.6%.Points in the lower panel show the ratios between the median values of the caustic and true MARs.
The figure shows that MAR C is a robust platform for estimating the true MAR 3D at every redshift.The median MAR C s are within 15% of the median MAR 3D s at each redshift.The caustic MARs are also unbiased relative to the true MARs.Thus the caustic technique provides accurate and robust estimates of the MARs of clusters in the redshift range 0.01 − 1.04.

Discussion
We use the Illustris TNG300-1 simulation to estimate the MAR of galaxy clusters based on the radial velocity profile of cluster galaxies and on the cluster total mass profile.The caustic technique provides robust and unbiased estimates of the true MARs derived from 3D data over the redshift range 0.01 − 1.04.
We next (Sect.6.1) compare the Illustris results with previous work on simulated and observed clusters.In Sect.6.2 we assess the bias between the galaxy MARs and MARs derived from the dark matter halos for the same sample of clusters drawn from the IllustrisTNG simulations.In Sect.6.3 we discuss future simulations and observational applications of the MAR.

Comparison with Previous Results
The dynamically motivated MAR recipe we develop differs significantly from previous merger tree approaches (e.g., McBride et al. 2009;Fakhouri et al. 2010;van den Bosch et al. 2014;Correa et al. 2015b;Diemer et al. 2017b;Diemer 2018).In contrast with a merger tree procedure that is not directly observable, the approach we outline allows the estimation of the MAR of real clusters and comparison with the true MARs of comparable simulated systems.
Most previous theoretical investigations of the MAR are based on N-body dark matter only simulations.These studies employ merger trees that trace the mass accretion of a halo at z = 0 back in time.The merger trees follow the change in mass between the halo descendant on the main branch identified at z i > 0 and its most massive progenitor for z i + ∆z where ∆z is the time step.The details of the simulation including the halo fragmentation algorithm and the choice of mass and time step may affect the results (Diemer 2018;Xhakaj et al. 2019).Some studies (e.g., van  Merger-tree MARs depend on the subhalo finder, the merger tree builder, and the definition of MAR.The difference between the two merger tree models in the upper panel of Fig. 10 shows the impact of these underlying differences.At each time step a merger tree MAR generally results from the difference between the mass of the descendant and the mass of the most massive progenitor.The most massive progenitor may not be the main branch progenitor (see Fakhouri et al. 2010, for more details) and thus the MAR is a lower limit.The halo mass definition may also vary.In the Fakhouri et al. (2010) model the mass is the sum of the masses of the subhalos; in Diemer et al. (2017b) the mass definition is M 200m , the mass enclosed within a sphere centered on the halo center with matter density equal to 200 times the background matter density.The choice of time step can also affect the MARs (Xhakaj et al. 2019).In Sect.6.2, we demonstrate that the use of dark matter only simulations may also bias the MARs toward lower values, but the bias is small.In all models, the MAR increases with cluster mass.Because Fakhouri et al. (2010) and Diemer et al. (2017b) do not report errors in their fits we assume fractional errors comparable with ours (shadowed areas).The slope of the IllustrisTNG MAR then agrees with Fakhouri et al. (2010) and Diemer & Kravtsov (2014).When averaged over the entire redshift range, the slope of the Illustris MARs as a function of M 3D 200c is 0.90 ± 0.18 (Sect.5) in agreement with the slopes of 1.1-1.2obtained by Fakhouri et al. (2010) and Diemer et al. (2017b).
The IllustrisTNG MARs are consistent with merger tree MARs at every redshift (upper panel of Fig. 10).Assuming similar fractional errors for Fakhouri et al. (2010) and Diemer et al. (2017b) and IllustrisTNG, the 0-50% difference in the rates is usually within the 1σ error and always within the 2σ error.The shaded regions in Fig. 11 indicate the general consistency of the results.The difference between the IllustrisTNG and merger tree MARs increases as the redshift increases, but so does the error (Table 4).The parameter W of Eq. ( 8) characterizes the redshift-dependence of the MAR.The IllustrisTNG simulations yield W = 1.69 ± 0.34 (Sect.5); Fakhouri et al. (2010) obtain W = 1.17.The difference is ≲ 2σ.The qualitative agreement between the IllustrisTNG and merger tree approaches is reassuring given the substantial differences in the algorithms.
Based on the spherical accretion prescription proposed by De Boni et al. (2016), Pizzardo et al. (2021) develop the first systematic approach to estimating the MARs of real clusters.Pizzardo et al. (2021) compute the MAR as the ratio between the mass of an infalling shell and the infall time.Their approach is similar to the IllustrisTNG approach we follow, but they use a ΛCDM N-body dark matter only simulation.Pizzardo et al. (2022) apply the Pizzardo et al. (2021) recipe to ten stacked clusters from the HectoMAP redshift survey (Sohn et al. 2021a,b).They also derive MARs based on the ΛCDM Nbody simulation L-CoDECS (Baldi 2012).The bottom panel of Fig. 10 shows the simulated MARs (triangles) and the observed MARs of HectoMAP stacked clusters (stars) as a function of mass and colour coded by redshift compared with TNG300-1 (squares).
The HectoMAP MARs are consistent with the TNG300-1 MARs.The agreement between MARs of observed clusters and TNG300-1 may reflect the use of simulated galaxies to calibrate the MAR derived from IllustrisTNG.

The Dark matter MAR
Previous theoretical work on the MAR is based on N-body dark matter only simulations.Galaxies are generally biased tracers of the underlying distribution of dark matter (e.g., Kaiser 1984;Davis et al. 1985;White et al. 1987).With Illustris TNG300-1 we can measure the bias by estimating the MAR for both galaxies and dark matter for the identical sample of clusters.We use 3D data from the simulation for this test.
We begin by locating the infalling shell based on the dark matter.We follow the procedure outlined in Sect.2.2 using the average radial velocity profiles of the dark matter particles.For each redshift bin, we compute the dark matter radial velocity profile of the individual clusters in 200 logarithmically spaced bins covering the radial range (0 − 10)R 3D 200c .We choose narrower binning for the dark matter profiles than we did for the galaxies because the number of dark matter particles is much larger than the number of galaxies.
We compute a single mean radial velocity profile and smooth it as we did in Sect.2.2.We identify the minimum radial velocity of the average profile, v dm min , and its cluster-centric location, R dm v min .As in 2.2, the boundaries of the infalling are the cluster-centric distances where the average velocity is 0.75v dm min .The average radial velocity profiles of dark matter and galaxies agree at every redshift.The radial location of v dm min , R dm v min , is on average ∼ 0.04% (∼ 7.9% at most) smaller than R v min based on the galaxies.The infall velocity based on the dark matter, v dm inf , is on average ∼ 2.1% (∼ 5.9% at most) less than v inf determined from the galaxies.
We compute the mass of the infalling shell, M 3D,dm shell , and the shell infall time, t dm inf , from the dark matter field following Sects.4.1 and 4.2 but based only on the mass profile of the dark matter component.To compute M 3D,dm shell and t dm inf we multiply the dark matter profile by (1 + Ω b0 /Ω m0 ) to account for the baryonic fraction.The resulting MAR 3D,dm are then directly comparable with the MAR 3D computed based on the total mass profile.
Figure 12 compares the dark matter and galaxy MARs.Blue squares and orange points in the upper panel show the median true MARs derived from galaxies and dark matter, respectively.The corresponding coloured error bars show the interquartile ranges of the MARs.Points in the lower panel show the median ratio between the dark matter and galaxies 3D MARs.The dash-dotted line show the global median.
On average, the dark matter MAR 3D,dm s are ∼ 6.5% below the MARs based on the galaxies.The scatter between the two MARs is < 30%, generally less than the uncertainty in the determination of the respective MARs.Thus galaxies are indeed biased tracers, but the bias is small on these scales.
The MAR derived from galaxies should exceed the dark matter MAR because the clustering amplitude of galaxies relative to dark matter is larger at smaller scales and at higher redshift (Davis & Peebles 1983;Davis et al. 1985;Jenkins et al. 1998).On the scale of the accretion region of galaxy clusters with redshift z ≲ 1, the galaxy clustering excess is, however, small (Springel et al. 2018).Thus the bias between the galaxy and dark matter MARs derived from IllustrisTNG is also small.
Next generation wide-field spectroscopic surveys will observe the infall regions around large numbers of galaxy clusters with high sampling rates.These dense and deep spectroscopic surveys will provide the necessary observational baseline for measuring the MARs of thousands of clusters extending to higher redshifts.
The multi-object William Herschel Telescope Enhanced Area Velocity Explorer spectrograph on WHT (WEAVE, Balcells et al. 2010;Dalton et al. 2012;Dalton et al. 2014;Dalton 2016) will explore the infall regions of galaxy clusters and their connections to the cosmic web.The Weave Wide Field Cluster survey (WWFC, Cornwell et al. 2022Cornwell et al. , 2023) ) will measure thousands of galaxy spectra in and around 20 clusters with 0.04 < z < 0.07 out to radii ≲ 5R 200c .A deeper cluster survey will provide dense spectroscopic surveys of 100 clusters for redshift ≲ 0.5, a new baseline for measuring the MAR in this redshift range.Planned observations with the Prime Focus Spectrograph on Subaru (PFS, Tamura et al. 2016) and the Maunakea Spectroscopic Explorer on CFHT (MSE, Marshall et al. 2019) will provide spectroscopic redshifts of hundreds to thousands of galaxy cluster members for thousands of individual clusters with z ≲ 0.6.
The caustic technique (Diaferio & Geller 1997;Diaferio 1999) and weak gravitational lensing (Bartelmann 2010;Hoekstra et al. 2013; Umetsu 2020) will provide two independent measurements of cluster mass profiles extending to large radii.Neither the caustic technique nor weak lensing relies on the assumption of dynamical equilibrium.These techniques can thus be applied throughout the accretion region where dynamical equilibrium does not hold (Ludlow 2009;Bakels et al. 2021).Present weak-lensing maps from HST and Subaru already provide cluster mass profiles up to 5.7 Mpc for ∼ 20 systems (Umetsu et al. 2011(Umetsu et al. , 2016)).Future facilities will extend these measurements to thousands of clusters.The VRO (Ivezić et al. 2019) and the Euclid mission (Laureijs et al. 2011;Sartoris et al. 2016) will provide extended weak-lensing mass profiles for combination with extensive spectroscopy samples of clusters with z ≲ 2.
The Illustris TNG300-1 MARs are based on < 100 clusters at z > 0.52.The mass distribution of TNG300-1 mostly samples clusters with M 3D 200c ∼ (1.2 − 2) • 10 14 M ⊙ .High redshift massive clusters have larger MARs and place tight constraints on models of structure formation and evolution (e.g., Bardeen et al. 1986;White 2002).Larger volume hydrodynamical simulations, including MillenniumTNG (Pakmor et al. 2023) with its 740 Mpc comoving size, will provide larger samples of the most massive systems up to higher redshift.The next generation of simulations should enable tracing of the MAR to redshifts ≳ 1. Extension of MAR determination to early epochs in cluster history will provide new insights into the astrophysics of cluster formation and evolution.

Conclusion
We use the Illustris TNG300-1 simulation (Pillepich et al. 2018;Springel et al. 2018;Nelson et al. 2019) to compute the MAR of clusters of galaxies.The recipe, based on the dynamics of cluster galaxies, computes the MAR as the ratio between the mass inside a spherical shell within the cluster infall region and the time for the shell to reach the virialized region of the cluster.The method builds on the approach by De Boni et al. ( 2016) and Pizzardo et al. (2021) and incorporates the caustic technique (Diaferio & Geller 1997;Diaferio 1999;Serra et al. 2011) that provides robust, unbiased estimates of the true MARs.IllustrisTNG refines estimation of clusters MARs.It permits derivation of true MARs in a narrow shell, an approach that provides improved scaling between observational and true MARs.A major goal of the approach is direct application to cluster observations.We use 1318 clusters extracted from TNG300-1 (Pizzardo et al. 2023).This sample includes both the 3D and caustic mass profiles of each cluster.The clusters have median mass M 3D 200c ∼ (1.3 − 1.6) • 10 14 M ⊙ and cover the redshift range 0.01 − 1.04.We locate the infalling shell based on the average radial motion of cluster galaxies as a function of cluster-centric distance and redshift.We compute the infall time by solving the equation for radial infall of the infalling shell to R 200c with nonconstant acceleration derived from the true cluster gravitational potential.
The MARs increase with increasing cluster mass and redshift.At fixed redshift, a change of ∼ 1 order of magnitude in M 200c yields a comparable increase in the MAR.This dependence tracks the increase of the mass of the infalling shell as a function of M 200c .At fixed mass, the MAR increases by a factor of ∼ 5 from z = 0.01 to z = 1.04 because of the anticorrelation of the infall time with redshift.The correlations between the MAR and cluster mass and redshift are predicted by hierarchical structure formation scenarios.
The MARs from IllustrisTNG build on similar approaches based on N-body simulations (e.g., Pizzardo et al. 2021Pizzardo et al. , 2022)).In Illustris TNG300-1 we can test the dark matter MARs against the galaxy MARs for the identical set of simulated systems.The dark matter MARs are ∼ 6.5% lower than the galaxy MARs reflecting the relative amplitudes of the clustering of galaxies and dark matter as a function of scale and redshift.
The IllustrisTNG MARs complement approaches based on merger trees which cannot be linked as directly to the observations (e.g., McBride et al. 2009;Fakhouri et al. 2010;Correa et al. 2015b;van den Bosch et al. 2014;Diemer et al. 2017b;Diemer 2018).The IllustrisTNG MARs lie within 2σ of the merger tree results.On average, the IllustrisTNG MARs exceed the merger tree MARs by ∼ 50 − 70%; the difference increases with redshift.At fixed redshift, the dependence of the merger tree and Illustris MARs on cluster mass agree well.The IllustrisTNG MARs are remarkably consistent with available observations of the MAR as a function of redshift (Pizzardo et al. 2022).
IllustrisTNG enables exploration of accretion by galaxy clusters with simulated galaxies.The approach provides a framework for obtaining robust MARs based on observed spectroscopic samples with ≳ 200 cluster members.Future spectroscopic (e.g. with instruments like WEAVE, PFS, and MSE) and weak lensing (e.g. with facilities like Euclid and VRO) surveys will provide these samples.The next generation of large volume hydrodynamical simulations including MillenniumTNG will guide the interpretation of observations of the MAR for higher redshift and and an extended cluster mass range.Determination of MARs from the combined large datasets and enhanced simulations will test models of formation and evolution of cosmic structure.

Fig. 1 .Fig. 2 .
Fig. 1.Galaxy radial velocity profiles, infalling shells, and turnaround radii.The solid blue curves show the average galaxy radial velocity profiles for redshifts z = 0.01, z = 0.62, and z = 1.04 (left, middle and right panels, respectively).Dashed blue curves show the Savitzky-Golay smoothed profiles based on a ten-bin window (Savitzky & Golay 1964).The shadowed region shows the error in the smoothed profile.The vertical black line shows the average turnaround radius based on the velocity profiles.In each panel the blue vertical band indicates the infalling shell with boundaries at cluster-centric radii where the velocity is 0.72v min (Sect.4).

Fig. 3 .
Fig. 3. Relation between MAR and MAR t .Curves show MAR 3D (A)/MAR 3D (A = 0.72) as a function of A in the range 0.45 − 0.95 for six redshifts from z = 0.01 to z = 1.04.The error bar shows the median interquartile range of the ratio at A = 0.72.
Fig. 4. M 3D shell as a function of M 3D 200c .Points show M 3D shell as a function of M 3D 200c in the three redshift bins z = 0.01, 0.62, and 1.04 (cyan, violet, and magenta, respectively).Squares with error bars indicate the median and interquartile range in eight fixed logarithmic mass bins.Lines show power law fits to the data.

Fig. 6 .
Fig. 6.Median and interquartile ranges of the ratio between the caustic and true M shell as a function of redshift.

Figure 5
Figure 5 shows the same trend observed in Fig. 4 between M 3D shell and M 3D 200c .At fixed redshift M 3D shell increases with M 3D 200c : a change of ∼ 1 order of magnitude in M 3D 200c corresponds to an analogous change in M 3D shell .The median M 3D shell s, denoted by the black circles, are uncorrelated with redshift.The dependence of M shell and redshift is expected in the standard model of structure formation and evolution.In Appendix A we explain how a complex interplay between the density and mass profiles as a function of redshift along with the slightly different cluster mass distribution at various redshifts tends to suppress any correlation between M shell and redshift.The caustic technique provides robust estimates of the mass profile of clusters in the infalling shell(Pizzardo et al. 2023).Points with error bars in Fig.6show the median and the interquartile range of individual ratios between the caustic and true mass M C shell /M 3D shell as a function of redshift.Analogously to M 3D shell , we define M C shell = K M C shell (A = 0.72).Because the infalling region is closer to the cluster center as the redshift increases (Table3, fifth column), the median ratio slowly rises from 0.9 at z = 0 to 1.1 at z = 1.04.1 The typical dispersion in the ratio is ∼ 34%.Thus the caustic technique slightly overestimates or underestimates the mass at smaller and larger cluster-centric distances, respectively.The trend is not statistically significant: Kendall's τ correlation coefficient is ≲ 0.05.The caustic mass M Cshell is an unbiased estimate of M 3D shell .Furthermore M C shell has the same dependence on mass, redshift, and R v min as M 3D shell on average.

Fig. 7 .
Fig. 7. t inf as a function of M 3D 200c .Points show t inf of individual clusters as a function of their mass M 3D 200c in the three redshift bins z = 0.01, 0.62, and 1.04 (cyan, violet, and magenta, respectively).Squares with error bars indicate the median and the interquartile range in eight fixed logarithmic bins of mass.Lines show power law fits to the data.

Fig. 8 .
Fig. 8. MAR 3D as a of M 3D 200c .Points show MAR 3D as a function of M 3D 200c in the three redshift bins z = 0.01, 0.62, and 1.04 (cyan, violet, and magenta, respectively).Squares with error bars show the median and interquartile range in eight fixed mass bins.Lines show power law fits to the data.

Fig. 9 .
Fig. 9. Comparison between MAR C and MAR 3D .Upper panel: Median true (blue squares) and caustic (red triangles) MARs as a function of redshift.The solid blue and dotted red error bars show the interquartile ranges of the true and caustic MARs, respectively.The black curve shows a fit to Eq. (8).Lower panel: Points show the median ratio between the caustic and true MARs as a function of redshift.
den Bosch et al. 2014; Correa et al. 2015a,b) are based on the analytic or semi-analytic extended Press-Schechter formalism (Bond et al. 1991) calibrated with N-body simulations.The upper panel of Fig. 10 compares MARs from Fig. 8 (squares) with two merger tree calculations based on N-body simulations by Fakhouri et al. (2010, Eq. 2) 3 (dotted lines) and Diemer et al. (2017b, Eqs.9-10) (solid lines).Fakhouri et al. (2010) extend the McBride et al. (2009) investigation of the Millennium simulation (Springel et al. 2005) to Millennium II (Boylan-Kolchin et al. 2009).For the mass of each FoF halo, they take the sum of the masses of their subhalos.They then compute the MARs from one snapshot to the next.The results by Fakhouri et al. (2010) are in excellent agreement with the earlier results of McBride et al. (2009) and the later results of van den Bosch et al. (2014, see their Fig.12) and of Correa et al. (2015a,b).Diemer et al. (2017b) use a large set of ΛCDM N-body simulations run with GADGET2 (Springel et al. 2005).They use M 200m as a halo mass proxy and compute the MARs at time steps comparable with the halo dynamical time (Diemer 2017a) 4 .

Fig. 10 .
Fig. 10.Comparison between our MARs and previously published results.Upper panel: Squares show MAR 3D s from Fig. 8 as a function of M 3D 200c and colour coded by redshift from dark blue to dark red as the redshift increases from bin to bin:, z = 0.01, 0.21, 0.52, 0.82, and z = 1.04.Lines show MARs from merger trees (Fakhouri et al. (2010) (dotted line) and Diemer et al. (2017b) (solid line).Lower panel: Squares again show the true MARs from TNG300-1 as a function of M 3D 200c .Triangles show the MAR 3D s from Pizzardo et al. (2022) as a function of M 3D 200c in four redshift bins: z = 0.01, 0.19, 0.44, and z = 1.00.Stars show the MARs of ten observed stacked clusters from the HectoMAP survey (Pizzardo et al. 2022).

Fig. 11 .
Fig. 11.Comparison between fits to our MARs and previously published fits.The solid, dotted, and dashed lines show the fit of the analytic relation proposed byFakhouri et al. (2010) to the Illustris MARs (Eq.(8), see Sect.5),Fakhouri et al. (2010), andDiemer et al. (2017b), respectively, for z = 0.31.The shadowed turquoise band shows the typical 1σ IllustrisTNG scatter.The gray shadowing indicates a similar assumed scatter for the merger tree models (grey bands).

Table 4 .
MAR 3D as a function of M 3D coefficient a that measures the MAR at fixed mass increases with redshift by a factor ∼ 7 from z = 0.01 to z = 1.04.The slope of the power law is essentially redshift independent.The average of the 11 slopes b's yields a mean slope b = 0.90 ± 0.18 in the redshift range 0.01 − 1.04.
Comparison between MARs from galaxies and dark matter.Upper panel: Blue squares (orange dots) and blue (orange) error bars show the median and the interquartile range of the true MARs based on galaxies (dark matter) as a function of redshift.Lower panel: Points show the median ratio between the dark matter and galaxy MARs as a function of redshift.The dash-dotted line shows the median ratio.