Mass accretion rates of the HectoMAP clusters of galaxies

We estimate the mass accretion rate (MAR) of the 346 clusters of galaxies in the HectoMAP Cluster Survey. The clusters span the redshift range $0.17-0.42$ and the $M_{200}$ mass range $\approx (0.5 - 3.5)\cdot 10^{14} M_\odot$. The MAR estimate is based on the caustic technique along with a spherical infall model. Our analysis extends the measurement of MARs for 129 clusters at $z<0.3$ from the Cluster Infall Regions in the Sloan Digital Sky Survey (CIRS) and the Hectospec Cluster Survey (HeCS) to redshift $z \sim 0.42$. Averaging over redshift, low-mass clusters with $M_{200}\sim 0.7\cdot 10^{14} M_\odot$ accrete $\sim 3\cdot 10^4 M_\odot$yr$^{-1}$; more massive clusters with $M_{200}\sim 2.8\cdot 10^{14} M_\odot$ accrete $\sim 1\cdot 10^5 M_\odot$yr$^{-1}$. Low- and high-mass clusters increase their MAR by $\sim 46\%$ and $\sim 84\%$, respectively, as the redshift increases from $z\sim 0.17-0.29$ to $z\sim 0.34-0.42$. The MARs at fixed redshift increase with mass and MARs at fixed mass increase with redshift in agreement with $\Lambda$CDM cosmological model for hierarchical structure formation. We consider the extension of MAR measurements to $z \sim 1$.

The outskirts of clusters of galaxies are a potentially powerful probe of the model of structure formation and evolution (e.g., Diaferio 2004;Reiprich et al. 2013;Diemer & Kravtsov 2014;Lau et al. 2015;Walker et al. 2019;Rost et al. 2021): Far from the center (at radii 2R vir 1 ), hydrodynamic simulations show that baryonic processes play a minor role in cluster dynamics (e.g., Diemand et al. 2004;van Daalen et al. 2014;Velliscig et al. 2014;Hellwing et al. 2016;Armitage et al. 2018;Shirasaki et al. 2018). This feature makes the comparison with simulations more robust because the relevant physical processes are dominated by gravity.
The mass accretion rate (MAR) of galaxy clusters offers a valuable tool for probing the outskirts of clusters. N -body simulations show that the infall region of clusters, where new material is falling onto the clusters for the first time, is located beyond ∼ 2R 200 (Ludlow 2009;Diemer & Kravtsov 2014;More et al. 2015;Diemer et al. 2017;Bakels et al. 2021;Xhakaj et al. 2020). The MAR is tightly linked to the cluster properties (Vitvitska et al. 2002;Gao et al. 2004;van den Bosch et al. 2005;Kasun & Evrard 2005;Allgood et al. 2006;Bett et al. 2007;Ragone-Figueroa et al. 2010;Ludlow et al. 2013;Diemer & Kravtsov 2014;More et al. 2015;M. Pizzardo et al 2022, in preparation). For example, consider the splashback radius, the average location of the first apocentre of infalling material (Adhikari et al. 2014). At the splashback radius, there is a sudden decrease of the radial logarithmic derivative of the density profile, d log ρ/d log r (Diemer & Kravtsov 2014;More et al. 2015). N -body simulations show that the location of the splashback radius depends on the cluster mass and on the MAR: larger masses and accretion rates are associated with smaller splashback radii (Diemer & Kravtsov 2014;More et al. 2015;M. Pizzardo et al 2022, in preparation).
Despite the importance of the cluster outskirts, observations are limited compared with the intensive study of clusters within their approximate virial radii, R 200 . There are two reasons for this contrast. Because these regions cover a large region on the sky and the number density of galaxies in the outer regions of clusters is much smaller than in the inner regions, the contrast with the foreground/background is much lower, and it is challenging to obtain substantial samples of cluster members at these large radii. Other observational techniques for estimating mass profiles at large cluster-centric distances are also limited. The X-ray surface brightness drops at large radii. Sunyaev-Zel'dovich (SZ) effect measurements are less sensitive at large radii because the gas pressure is low. The weak-gravitational-lensing signal is more sensitive to projection effects at large radius (Serra et al. 2011;Ettori et al. 2013;Reiprich et al. 2013).
Most techniques for estimating the mass profiles of clusters assume dynamical equilibrium (Zwicky 1937;The & White 1986;Merritt 1987;Sarazin 1988;Pierpaoli et al. 2003;Rasia et al. 2006;Serra et al. 2011;Ettori et al. 2013;Reiprich et al. 2013). Only weak gravitational lensing (Bartelmann 2010;Hoekstra et al. 2013;Umetsu 2020) and the caustic technique (Diaferio & Geller 1997;Diaferio 1999;Serra et al. 2011) avoid the assumption of dynamical equilibrium. Both techniques can be applied at larger radii. Weak-lensing estimates are also redshift dependent, with a signal that peaks at intermediate redshift (Hoekstra 2003;Hoekstra et al. 2011). Both methods are steadily improving with the acquisition of large combined spectroscopic and photometric datasets. There are also advances in lensing-mass reconstruction methods. Umetsu et al. (2011) andUmetsu (2013) improve the precision of cluster-mass measurements from weak lensing by ∼ 30% with the joint use of lensing distortion and magnification. Using this method on the CLASH survey (Postman et al. 2012), Umetsu et al. (2014 estimate the clusters' mass profiles for radii 3 Mpc. Umetsu et al. (2016) include strong gravitational lensing shear and magnification and extend the estimates to even larger distances, ∼ 5.7 Mpc.
Although the caustic technique is also subject to projection effects, it returns an unbiased mass estimate with a relative uncertainty of 50% at large radii for sufficiently densely sampled systems. The caustic technique is independent of redshift (Serra et al. 2011), but obtaining large samples of cluster members at higher redshift is more demanding. In principle, weak lensing and the caustic technique offer complementary approaches that combine to elucidate the physics of the outer regions of clusters. Here we focus our attention on the application of the caustic technique.
Taking advantage of the Cluster Infall Regions in the Sloan Digital Sky Survey (CIRS) (Rines & Diaferio 2006) and the Hectospec Cluster Survey (HeCS) (Rines et al. 2013) spectroscopic surveys (now included in the HeCS-omnibus; Sohn et al. (2020)), Pizzardo et al. (2021) first estimated the MARs of 129 clusters in the redshift range 0.01 < z < 0.3 and mass range M 200 ∼ (0.1 − 25) · 10 14 M . They adopt a procedure based on the caustic technique to estimate the mass profile of the clusters, and they apply a spherical infall model (De Boni et al. 2016) to estimate the cluster MARs. The MAR estimation method is informed by N -body simulations.
The MAR estimates of Pizzardo et al. (2021) are limited to redshifts z 0.3. Extension to higher redshift offers the opportunity for more sensitive tests of models for the growth of structure in the universe. We thus explore the additional constraints obtained from an independent cluster sample reaching z ∼ 0.42. Sohn et al. (2021a) identified 346 clusters with redshift 0.17 z 0.42 by applying a Friends-of-Friends (FoF) algorithm to the full HectoMAP redshift survey (Geller et al. 2011;Hwang et al. 2016). The HectoMAP survey includes ∼ 110, 000 galaxies with spectroscopic redshifts with z < 0.6 and with an average redshiftẑ ∼ 0.31 (see Sohn et al. 2021b, for the first data release).
In Section 2 we review the HectoMAP cluster sample. In Section 3 we briefly introduce the recipe we use for estimating of the MAR. We then measure the MARs of the HectoMAP clusters. In Section 4 we discuss the results. We conclude in Section 5. We use the standard ΛCDM parameters Ω m0 = 0.27, Ω Λ0 = 0.73, and H 0 = 70 km s −1 Mpc −1 .

HECTOMAP CLUSTER SAMPLE
HectoMAP is a dense redshift survey designed to study galaxy clustering at redshifts 0.2 < z < 0.6 (Geller et al. 2011;Hwang et al. 2016;Sohn et al. 2021b). Sohn et al. (2021a) use an FoF algorithm to identify galaxy clusters in HectoMAP. Sohn et al. (2021b) and Sohn et al. (2021a) describe the details of the HectoMAP redshift survey and the HectoMAP FoF cluster catalog. We briefly review the HectoMAP survey in Section 2.1 and the HectoMAP FoF cluster catalog in Section 2.2.

The HectoMAP Redshift Survey
HectoMAP is a large-scale redshift survey covering 54.64 deg 2 over a narrow strip of the sky at 200 < R.A. (deg) < 250 and 42.5 < Decl. (deg) < 44.0. HectoMAP is based on Sloan Digital Sky Survey (SDSS) Data Release (DR) 16 photometry (Ahumada et al. 2020). The HectoMAP redshift survey includes a modest number of redshifts from SDSS DR16. Sohn et al. (2021b) measured most of the redshifts with the Hectospec instrument mounted on the MMT 6.5 m telescope (Fabricant et al. 1998(Fabricant et al. , 2005. The primary targets of HectoMAP are galaxies with r < 20.5 and (g − r) > 1, and galaxies with 20.5 ≤ r < 21.3, g − r > 1, and r − i > 0.5. The full survey is > 80% complete at r = 20.5 and 62% complete at r = 21.3. The completeness of the survey is much less for bluer galaxies outside the target color range (Sohn et al. 2021b). For red objects within the selection limits, the HectoMAP survey has remarkably uniform completeness on the sky. This uniformity results from the strategy of revisiting every position in the field typically ∼ 9 times (see Section 2.4 of Sohn et al. 2021b for details). The uniformity of HectoMAP makes the survey a robust basis for the examination of properties of dense systems and their evolution.
HectoMAP includes a total of ∼ 110, 000 galaxies with spectroscopic redshifts; the average galaxy number density is ∼ 2000 galaxies deg −2 . The typical uncertainty of HectoMAP redshifts is ∼ 40 km s −1 . The high density and uniform completeness of the survey enable the identification and study of HectoMAP galaxy clusters based on spectroscopy (Sohn et al. 2021a).

HectoMAP FoF Cluster Catalog
Sohn et al. (2021a) built a cluster catalog by applying an FoF algorithm to the HectoMAP redshift survey. The FoF algorithm bundles sets of neighboring galaxies within given linking lengths into candidate clusters. Sohn et al. (2021a) apply the FoF algorithm in redshift space. The standard FoF algorithm requires two linking lengths: one in the projected direction (∆D) and one in the radial direction (∆V ). Sohn et al. (2021a) choose the optimal linking lengths empirically using photometrically identified redMaPPer clusters (Rykoff et al. 2014(Rykoff et al. , 2016. The redMaPPer cluster catalog includes 104 redMaPPer clusters in Hec-toMAP. Sohn et al. (2021a) demonstrate that the FoF algorithm with ∆D = 900 kpc and ∆V = 500 km s −1 identifies more than 90% of the HectoMAP redMaPPer clusters. With these linking lengths, the FoF algorithm also identifies all of the 15 known X-ray clusters in Hec-toMAP (Sohn et al. 2018).
The HectoMAP FoF cluster catalog includes 346 systems with 10 or more spectroscopic members. Sohn et al. (2021a) report properties of the FoF clusters including the number of FoF members, the cluster center, and the cluster velocity dispersion. The HectoMAP FoF clusters typically consist of 17 members.
We use the center of each cluster as listed in Sohn et al. (2021a). The centers of the FoF clusters are determined based on the center-of-light method (Robotham et al. 2011), which weights the positions of FoF members by their luminosity. This method identifies a bright central galaxy. Sohn et al. (2021a) identify the position and redshift of this central galaxy as the FoF cluster center.
Finally, Sohn et al. (2021a) use the biweight technique (Beers et al. 1990) to compute the line-of-sight velocity dispersion of the cluster.
We select spectroscopic galaxies within R cl < 5 Mpc and |∆cz/(1 + z cl )| < 5000 km s −1 of each cluster center as a basis for measuring the mass accretion rate. The line-of-sight velocity boundaries are 5 times larger than the velocity dispersions of HectoMAP clusters. The radial boundaries are much larger than the region where we measure the mass accretion rate (see Figure 1, and Tables 1 and 2). Thus, the selection of the survey boundaries does not impact our analysis. In Section 3.1 we briefly review the procedure for estimating the MAR from spectroscopic catalogs (Pizzardo et al. 2021). In Section 3.2 we describe the sample of stacked HectoMAP clusters and their mass profiles. In Section 3.3 we estimate the MAR of the HectoMAP clusters and compare the results with previous measurements and with ΛCDM predictions. . This approach assumes that a spherical infall model describes the accretion of new material onto a cluster. During the infall time t inf , a cluster accretes all of the material within a spherical shell centered on the cluster center and with an inner radius R i and thickness δ s .

The MAR Recipe
We derive the shell thickness from the classical equation of motion for accretion with constant acceleration and initial velocity v inf , is the mass of the cluster within a radius R i , and G is the gravitational constant.
Given the mass of the infalling shell, M shell , the estimate of the MAR is where v p is the proper peculiar velocity and r c,i is the comoving position vector of the particle relative to the cluster center. At the snapshot redshift z s , H(z s ) and a(z s ) are the Hubble parameter and the cosmic scale factor. The minimum of the velocity profile is generally within the radial range [2, 2.5]R 200 . This range thus provides the most reliable estimate of the MAR. We adopt v inf as the value of the infall velocity for the median profile in the radial range [2, 2.5] R 200 .
The radial position of v inf , and thus of the infall region of a cluster, is coincident with the approximate location of the splashback radius (Diemer & Kravtsov 2014;More et al. 2015). N -body simulations show that the splashback radius separates the inner region where accretion is nearly complete from the outer region where accretion is ongoing (Diemer et al. 2017;Xhakaj et al. 2020).
Defining the observable MAR relies on the v inf derived from N -body simulations. This approach is necessary because measurement of the radial velocity profiles of observed cluster galaxies is not currently feasible. The results are insensitive to the exact choice of v inf . For example, Pizzardo et al. (2021) show that varying v inf within ±40% of the true value produces results within the typical ∼ 40% spread of the resulting MARs.

The Caustic Technique
Use of Eqs. (1) and (2) requires an estimate of the cluster mass at large cluster-centric distances where virial equilibrium does not hold. Pizzardo et al. (2021) use the caustic method (Diaferio & Geller 1997;Diaferio 1999;Serra et al. 2011). At large distances, (0.6−4)R 200 , the caustic method returns an unbiased estimate of the mass with better than 10% accuracy and with a relative uncertainty of 50% (Serra et al. 2011) provided that the velocity field of the cluster outer region is sufficiently well sampled. According to Serra et al. (2011), ∼ 200 spectroscopic measurements within a three-dimensional distance ∼ 3R 200 provide a sufficient sample. Increased sampling generally improves performance. Significantly sparser samples lead to underestimation of the mass. For example, a sample of only ∼ 100 spectroscopically identified members results in a mass underestimated by ∼ 35%. Pizzardo et al. (2021) demonstrate that, in ΛCDM simulations, the MARs estimated by applying the spherical accretion recipe based on caustic mass profiles are unbiased. The MARs are within ∼ 19% of the MARs computed by applying the same accretion recipe to the true mass profiles of large samples of simulated clusters.
The caustic technique estimates the three-dimensional mass profile of a cluster from the line-of-sight escape velocity of the cluster galaxies as a function of clustercentric distance R, v 2 esc,los (R). The R − v los diagram, the line-of-sight velocity relative to the cluster median as a function of R, is the basis for the extraction of the caustic mass profile. In this diagram, the cluster galaxies appear in a well-defined trumpet-shaped pattern; the amplitude decreases as R increases. The vertical separation between the upper and lower caustics at radius R is the caustic amplitude, A(R), which approximates the average line-of-sight escape velocity profile. The caustic technique locates the caustics from the R − v los diagram. The square of the caustic amplitude, A 2 (R), estimates v 2 esc,los (R). With a form (or filling) factor F β accounting for the cluster velocity anisotropy, 2 the caustic technique provides an estimate of the threedimensional escape velocity profile, v 2 esc (R), which is related to the three-dimensional cluster gravitational potential φ, v 2 esc (R) = −2φ. The resulting estimate of the mass profile is where G is the gravitational constant and F β = 0.7 (see Diaferio & Geller 1997;Diaferio 1999;Serra et al. 2011, for further details).

Stacking the Clusters
A sufficiently dense spectroscopic catalog is a fundamental requirement for the application of the caustic technique. The HectoMAP cluster catalogs include a median of 102 member galaxies within a projected cluster-centric distance of 5 Mpc, with a 68th percentile range (65-157). For these clusters, R 200 ∼ 1 Mpc and thus the HectoMAP catalogs include a median of ∼ 50 galaxies within a projected distance of 3R 200 . This population is an upper limit to the number of galaxies within the three-dimensional cluster-centric distance, ∼ 3R 200 . In general, individual HectoMAP systems are not sufficiently well sampled for optimal performance of the caustic technique.
Given the limitations of redshift sampling for individual HectoMAP clusters, we estimate the MAR of HectoMAP clusters as a function of their velocity dispersion and redshift by stacking the observed clusters. We construct 10 stacked samples in five redshift bins, in the range [0, 0.6], and two velocity dispersion bins, [200, 400) km s −1 and [400, 1020] km s −1 . We select the FoF clusters to build these samples based on the velocity dispersion of each cluster (Sohn et al. 2021a, see Section 2.2). Our procedure removes 25 FoF groups with velocity dispersion < 200 km s −1 . These systems have masses 10 13 M . The stacked clusters include 321 individual FoF clusters.
We compute the velocity dispersion of each ensemble cluster by applying the biweight technique (Beers et al. 1990) to the ensemble catalog of relative line-of-sight velocities of galaxies with respect to the central galaxy of the individual FoF cluster. In this computation, we identify the mean line-of-sight velocity of a constituent FoF cluster with the redshift of its central galaxy. We use the bootstrap technique to derive the uncertainty in the ensemble velocity dispersion. Table 1 lists the total number of galaxies, the median redshift and its 68th percentile range, and the velocity dispersion for each of the stacked clusters.
We then apply the caustic technique to the 10 stacked clusters in Table 1. For each ensemble cluster, we build an R − v los diagram containing all of the galaxies with projected distance from their respective host FoF cluster center < 5 Mpc and with absolute line-of-sight velocity relative to the host cluster mean < 5000 km s −1 . We use centers of the individual FoF clusters listed in Sohn et al. (2021a) who identify the mean FoF cluster redshift with the redshift of the central galaxy. Finally, we derive the caustics for each stacked R − v los diagram and compute the related mass profiles according to Eq. (3).
Following Pizzardo et al. (2021), the R − v los diagrams of the stacked clusters are based on the projected cluster-centric distances and line-of-sight velocities of the galaxies in the reference frame of each of the individual systems included in the stack. In other words, there is no additional normalization. This approach is robust if the bins in velocity dispersion are small enough to avoid the presence of systems with widely different masses. The procedure avoids introducing systematics from the uncertain individual constituent cluster estimate of R 200 and velocity dispersion. Adopting a massvelocity dispersion relation (e.g., Evrard et al. 2008), the 90th percentile range of the clusters in the low-and highvelocity dispersion bins contains systems with masses in the range (0.1 − 0.8) · 10 14 M and (1 − 4) · 10 14 M , respectively. Figure 1 shows the R − v los diagrams of the 10 stacked clusters of Table 1 and the caustics. Table  Table 1 This stacking method, adopted by Rines & Diaferio (2006) and Serra et al. (2011), is particularly suitable with the HectoMAP FoF catalog because the R − v los diagram of a stacked cluster is generally dense enough to guarantee good performance of the caustic technique. Furthermore, the caustic technique assumes spherical symmetry for the cluster member distribution, but individual clusters are often triaxial (Frenk et al. 1988;Dubinski & Carlberg 1991;Warren et al. 1992). Stacking the member galaxies within multiple clusters averages out the asphericity of individual clusters (Rines & Diaferio 2006;Serra et al. 2011). Pizzardo et al. (2021) demonstrate the validity of this procedure by showing that, at comparable mass and redshift, the MARs of the stacked clusters agree with the average MARs of individual systems. They also show that the method is very robust against the overrepresentation of the richest individual systems and that it is insensitive to the impact of fore-and background galaxies. The red selection of HectoMAP clusters members (Sect. 2.1) does not affect our measurements: Pizzardo et al. (2021) show that missing ∼ 35% of the blue galaxies has no impact on the MARs. Rines et al. (2013) also show that blue galaxies have undetectable effects on dynamical mass estimates.

Estimates of the MAR
In addition to the caustic mass profiles of the 10 stacked HectoMAP clusters, we also require the initial radial velocities of the infalling shells to estimate the MAR. To obtain these velocities, we use numerical simulations (Sect. 3.1).
Following Pizzardo et al. (2021), we use the ΛCDM run of the L-CoDECS N -body simulations (Baldi 2012). The simulation has a box size of 1.43 Gpc in comov-ing coordinates. The simulation includes two collisionless fluids of 1024 3 particles each, with masses m DM = 8.34 × 10 10 M and m b = 1.67 × 10 10 M , respectively. The simulation is normalized at the cosmic microwave background epoch, with cosmological dark matter density Ω m0 = 0.226, cosmological constant Ω Λ0 = 0.729, baryonic mass density Ω b0 = 0.0451, Hubble constant H 0 = 70.3 km s −1 Mpc −1 , power spectrum normalization σ 8 = 0.809, and power spectrum index n s = 0.966. An FoF algorithm identifies groups and clusters in the simulations. The system centers are identified as the most bound particle within the system.
The ΛCDM run of L-CoDECS is a standard collisionless N -body simulation. For measuring the MAR, the simulation must trace the galaxy velocity field reliably in the outer regions of clusters, beyond 2R 200 : Nbody/hydrodynamical simulations (e.g., Diemand et al. 2004;Hellwing et al. 2016;Armitage et al. 2018) show that indeed the velocity bias between the velocity dispersions of galaxies and of dark matter particles is negligible in the outskirts of galaxy clusters. Hence, collisionless N -body simulations provide an unbiased measure of the MAR.
Here, we consider the same 12 samples of simulated clusters used by Pizzardo et al. (2021, see their Sect. 3.1 and Table 1). At z = 0, they include two samples of 2000 and 50 halos with median masses M 200 1.43×10 14 M and M 200 1.43×10 15 M , respectively. They trace the main progenitors of these halos at five higher redshifts with z ≤ 0.44. Over this redshift range, six low-mass bin samples cover the mass range ∼ (0.6 − 1.6) · 10 14 M and six high-mass bin samples cover the range ∼ (4 − 16) · 10 14 M .
We derive the median radial velocity profiles of the 12 samples of simulated clusters. We then use these radial profiles to compute the radial infall velocity, v inf , for   each stacked cluster according to its redshift and mass. The procedure involves three interpolations in velocityredshift, mass-redshift, and velocity-mass space. We also compute the uncertainty in v inf , but we ignore the uncertainty in v inf when estimating the MAR with Eq.
(2) because this uncertainty is negligible compared with the uncertainty in the mass profile. Sect. 3.3 of Pizzardo et al. (2021) contains a more detailed description of this procedure. For each simulated cluster, Table 2 lists the infall velocities and their uncertainty. The solid diamonds in Fig. 2 show the resulting MARs of the 10 stacked HectoMAP clusters. Table 2 lists the MARs. For comparison, the four crosses with error bars show the MARs of the four stacked clusters Pizzardo et al. (2021) constructed from the CIRS and HeCS surveys with the same stacking method we apply to Hec-toMAP. To compare the observations with the simulations, we separate each of the 12 simulated samples into five mass bins. We then have 30 low-mass and 30 highmass subsamples of halos. For each of these subsamples we compute the median of all of its individual halos MARs obtained by applying the recipe for the MAR to the true mass profiles. Fig. 2 (open circles) shows the results.
The MARs of the stacked clusters from HectoMAP are fully consistent with the MARs from the independent sample of stacked clusters from the CIRS and HeCS surveys. The MARs also appear consistent with ΛCDM expectations. In particular, the results underscore the positive correlation between the MAR and mass at fixed redshift and between the MAR and redshift at fixed mass (Pizzardo et al. 2021).
The uncertainties in the HectoMAP MAR estimates differ from the uncertainties for the CIRS and HeCS results. The mean relative uncertainty for the HectoMAP estimates is ∼ 44%; the CIRS-HeCS uncertainty is typically ∼ 64%. This difference originates from the difference in the number of galaxies within the R − v los plane of the stacked clusters: there are ∼ 2000 − 6000 galaxies for HectoMAP (Table 1) and ∼ 12000 − 27000 for CIRS-HeCS (see the 'Total' column of Table 7 in Pizzardo et al. 2021). Pizzardo et al. (2021) explain this effect. For a larger number of galaxies within the R − v los plane of the stacked cluster, the density contrast between the galaxy cluster members and the foreand background is smaller. The distribution of foreground and background galaxies also tends to be more uniform. These effects increase the uncertainty in the caustic mass profile because the boundary between the cluster and its surroundings is less well defined. Consequently the MAR is also more uncertain. However, the larger number density in the R − v los plane affects only the uncertainties in the mass profiles. The mass profiles themselves and the MARs are substantially unaffected. In other words, the MAR estimates based on the caustic mass profiles of stacked clusters are robust. Pizzardo et al. (2021) argue that the robust correlations among the MAR, M 200 , and redshift are linked directly to the correlations between M 200 , redshift, and the mass of the infalling shell with fixed radial thickness. Unlike the MAR, the mass of the infalling shell is unrelated to v inf . In the analysis of Pizzardo et al. (2021), the spherical infall model returns an average thickness δ s R i ≈ 0.5R 200 . Thus in addition to considering the individual values of the thickness of the infalling shells, Pizzardo et al. (2021) also consider the mass M 2−2.5 of the spherical shell with inner and outer radii 2R 200 and 2.5R 200 .
We also compute M 2−2.5 for each HectoMAP stacked cluster. We compare the results with the CIRS and HeCS stacked clusters and with the ΛCDM predictions. Figure 3 shows the relation between M 2−2.5 , M 200 , and redshift. As in Fig. 2, the M 2−2.5 s for the HectoMAP stacked clusters are consistent with those from CIRS and HeCS and with the estimates based on true mass profiles from ΛCDM simulations. Comparison between Fig. 3 and Fig. 2 confirms the correlation between MAR and M 2−2.5 for z 0.3, but M 2−2.5 cannot be used as a proxy for the MAR at higher redshift. Strikingly, the simulated data (empty circles) in Fig. 3 have a distribution that differs from the one in Fig. 2 for the z ∼ 0.35 − 0.44 range that exceeds the largest redshifts investigated by Pizzardo et al. (2021). At these larger redshifts, the clear increase of the MAR with increasing redshift at fixed mass is absent for M 2−2.5 . The dependence of the MAR on redshift revealed by using shells with increasing width for larger redshifts is an expected feature of hierarchical structure formation. Starting with two halos of the same mass at different redshifts, the halo at greater redshift grows faster than the halo at lower redshift; this difference occurs because the higher-redshift halo is embedded within a higher-overdensity region. This result indicates that, in general, the radial velocity of the infalling region cannot be ignored in estimating the MAR.
To underscore this point, Fig. 4 shows the thickness of the infalling shell, δ s R i , as a function of halo redshift and mass. The thickness, δ s R i , increases by ∼ 200% and ∼ 75% from z = 0 to z = 0.44, for low-and highmass halos, respectively. At fixed mass, the shell thickness depends mainly on the infall velocity. Fig. 4 thus shows that the absolute value of v inf steadily increases with redshift. The radial velocity profiles of halos confirm this behavior. This behavior of v inf also explains why the correlation between MAR and M 2−2.5 weakens as the redshift increases. The MAR increases with red- shift because the infalling shells of high-redshift halos are consistently thicker than those for the low-redshift counterparts.

DISCUSSION
The independent sample of HectoMAP clusters reproduces the observed MARs of the HeCS and CIRS clusters (Pizzardo et al. 2021). The HectoMAP clusters extend the redshift range of observed MARs to z ∼ 0.42. These observed MARs are consistent with ΛCDM expectations.
Here we discuss the potential for reducing the error in the MAR and thus increasing its power to constrain models for structure formation (Sect. 4.1). We also highlight basic model predictions for higher accretion at greater redshift. Observations will soon access these redshifts, again increasing in the power of the MAR as a discriminant among models (Sect. 4.2). Future approaches to computing the MAR will also benefit from enhanced hydrodynamical simulations covering large volumes (Sect. 4.2).

Observational Challenges
Combining weak gravitational lensing with the caustic technique holds promise for reducing the uncertainty in the MAR estimates. This ambitious project requires deep and dense spectroscopic surveys along with extensive imaging surveys. Imaging surveys already exist and dense spectroscopic surveys of clusters over a large redshift range will be possible with, e.g., the Prime Focus Spectrograph (PFS) on Subaru (Tamura et al. 2016).
Deep, dense spectroscopic surveys of clusters can substantially reduce both systematic and statistical errors affecting weak-lensing mass estimates. The spectroscopic surveys can robustly identify potentially contaminating foreground and background structures, thus enabling correction for their effects. Very large, deep spectroscopic surveys also promise to limit the systematic error in the source redshift currently based on photometric redshifts (von der Linden et al. 2014).
At low redshifts, 0.1, spectroscopic surveys of massive clusters and the weakly lensed background galaxies are feasible. These surveys enable spectroscopic tomographic weak lensing. Dell'Antonio et al. (2020), using spectroscopic redshifts of A2029 (Sohn et al. 2019), at z = 0.078, extracted the tangential ellipticity of the background galaxies with z 0.8. The weak-lensing mass agrees with the X-ray mass. This approach avoids calibration issues inherent in the use of photometric redshifts for background sources, and it effectively removes contamination of the lensing signal by faint cluster members.
The signal-to-noise (S/N) of this method is currently limited by small samples of redshifts for faint galaxies. Powerful spectrographs such as the PFS on Subaru will provide spectroscopic catalogs of clusters with ∼ 6000 cluster members and ∼ 15, 000 background objects; such catalogs will allow spectrotomographic weaklensing measurements with S/Ns comparable to current photometric-redshift-based weak-lensing measurements for hundreds of galaxy clusters (Dell'Antonio et al. 2020). This method holds promise for improving the weak-lensing precision of cluster-mass profiles at large radius where the MAR can be derived.
At higher redshift, a synergy between weak-lensing and caustic mass profiles could result in a substantial decrease of the uncertainty in the MARs. With dense spectroscopic redshift catalogs of clusters, the caustic mass profiles are unbiased because projection effects decrease. However, they have large uncertainties that do not decrease significantly with better sampling. In contrast, weak-lensing mass profiles are increasingly precise (e.g., Umetsu et al. 2011;Umetsu 2013), and some of the inherent biases will be reduced with large spectroscopic surveys. Taking advantage of the respective strengths of the two methods promises much more robust, unbiased, and accurate mass profiles extending to a large enough radius for the determination of the MAR.
Powerful instruments such as the PFS on Subaru (Takada et al. 2014) have the potential to return large sets of members of individual galaxy clusters. These surveys would allow the determination of MARs for individual clusters. For HectoMAP, however, we must stack the clusters to obtain large samples (Rines & Diaferio 2006;Biviano & Poggianti 2010;Serra et al. 2011;Biviano 2020;Pizzardo et al. 2021) and to obtain the most robust possible results (Sect. 3.2). Stacking does have the advantage that it averages over spatial and kinematic anisotropies. Extensive redshift surveys like DESI (Dey et al. 2019) will yield large samples of sparsely sampled clusters, potentially yielding stacked systems similar to those we construct from the HectoMAP. All of these future surveys will access systems at higher redshift.

Theoretical Predictions and Challenges
Deeper redshift surveys allow estimation of the MAR at higher redshifts where the MAR at fixed cluster mass should be larger (van den Bosch 2002;McBride et al. 2009;Fakhouri et al. 2010;van den Bosch et al. 2014;Correa et al. 2015;Diemer et al. 2017). The HectoMAP sample extends to z ∼ 0.42. Here we use simulations to briefly explore the MAR at higher redshift. Fig. 5 shows the MARs of synthetic clusters in the redshift range z = [0, 1]. In the low-and high-mass samples (bottom left and top right, respectively), the six circles at lower redshifts, 0 − 0.44, use the same simulation data as in Fig. 2 (Sect. 3.3), but at each redshift we use only one mass bin for each sample. The upper and lower dark blue circles indicate the median MARs as a function of the median M 200 of the progenitors at z = 1 of the halos in the high-and low-mass bin, respectively. The error bars show the 68th percentile ranges of mass and MAR. The MARs of high-redshift halos are much greater than the MARs of present epoch halos with similar mass: at z = 1 the progenitors of the high-mass halos have masses ∼ (1 − 4) · 10 14 M . This mass range overlaps the low-mass bin at lower redshifts. However, the MARs of the z = 1 halos are an order of magnitude larger than their counterparts at z = 0.
At higher redshift, the MAR may be more sensitive to the nature of dark matter, dark energy, and gravity (Cimatti et al. 2008;Guzzo et al. 2008;Baldi 2012;Candlish 2016). New-generation simulations including the hydrodynamical N -body Illustris TNG simulation suite (Springel 2010;Weinberger et al. 2017;Pillepich et al. 2018) are a crucial platform for exploring the entire range of physical insights possible with the determination of the MAR at greater redshift.
The Illustris TNG data release (Nelson et al. 2019) provides catalogs of galaxies for cleaner comparison with the data. These galaxy catalogs enable the computation of the MAR of simulated clusters based on galaxies rather than the underlying dark matter field. The comparison between the models and the data includes the following improvements: (i) simulated galaxies are a better proxy for observed galaxies, (ii) matter ) and velocity (Ye et al. 2017;Kuruvilla et al. 2020) biases that could affect dark-matter-only estimators are minimized, and (iii) the center of a cluster can be identified with a bright galaxy (Sandage & Hardy 1973;Dressler 1979;Dubinski 1998) or by applying the center-of-light method to cluster member galaxies (Robotham et al. 2011) in analogy with observed systems.
Larger simulated volumes will permit increases in the statistical significance of theoretical constraints. Large volumes are important for obtaining adequate samples of clusters at the upper end of the mass function. These massive systems place tight constraints on structure formation models (Bardeen et al. 1986;White 2002;Courtin et al. 2011;Cui et al. 2012;Giocoli et al. 2018). These systems are also the candidate systems most likely to be accessible observationally.

CONCLUSION
We derive MARs that characterize 321 clusters in the HectoMAP Cluster Survey (Sohn et al. 2021a). The cluster sample covers the redshift range 0.17 z 0.42 and mass range M 200 ≈ (0.5 − 3.5) · 10 14 M .
To estimate the MARs, we adopt the approach of Pizzardo et al. (2021). They apply the caustic technique (Diaferio 1999;Serra et al. 2011) to estimate the mass profiles of the clusters at ∼ 2−3R 200 , where accretion occurs. They adopt a spherical accretion prescription (De Boni et al. 2016) to evaluate the MAR. N -body simulations show that this recipe returns MARs within 19% of the MARs obtained by applying the same recipe to the true mass profiles of synthetic clusters. The technique is robust against the typical photometric and spectroscopic incompleteness of spectroscopic redshift surveys (Pizzardo et al. 2021).
The HectoMAP cluster MARs agree well with the MARs derived (Pizzardo et al. 2021) for independent samples of lower-redshift clusters, CIRS (Rines & Diaferio 2006) and HeCS (Rines et al. 2013) at z 0.3. In the low-velocity dispersion bin, with M 200 ∼ 0.7 · 10 14 M , the average MAR is ∼ 3 · 10 4 M yr −1 , whereas in the high-velocity dispersion bin, with M 200 ∼ 2.8 · 10 14 M , the average MAR is ∼ 1 · 10 5 M yr −1 . The uncertainty in the HectoMAP cluster MAR measurements is ∼ 44%, roughly one-third smaller than that of Pizzardo et al. (2021). We show that an increase in the width of the mass shell with redshift is critical for obtaining the correct MAR at larger redshift.
We compare the HectoMAP MARs with ΛCDM using the L-CoDECS N -body simulation suite (Baldi 2012  Larger cluster catalogs with more densely sampled systems reaching greater redshift and a more extensive cluster-mass range will provide more sensitive measures of the MAR. The dependence of the MAR on both redshift and cluster mass can provide new insights into the development of structure in the universe. Higherredshift catalogs extend the MAR probes to epochs where the mass accretion is larger and the sensitivity to the details of the ΛCDM paradigm is more pronounced (Cimatti et al. 2008;Guzzo et al. 2008;Baldi 2012;Candlish 2016). These measures can complement approaches based on large-scale statistical analyses.
Large dense surveys obtained with heavily multiplexed spectrographs on large telescopes (e.g. the PFS on Subaru; see Tamura et al. 2016) are crucial for reducing the uncertainty in the measurement of the MAR. The exploitation of weak-gravitational-lensing mass profiles estimated with sophisticated mass reconstruction methods (e.g., Umetsu et al. 2011;Umetsu 2013) combined with caustic profiles will improve the accuracy of mass profiles at large distances from the cluster center.
Finally, the new generation of hydrodynamical simulations like the Illustris TNG suite (Springel 2010;Weinberger et al. 2017;Pillepich et al. 2018) improves theoretical constraints on cluster accretion, including the high-redshift regime. These simulations allow the computation of the MAR directly from galaxies, thus improving consistency with results from observations. Larger simulated volumes will allow the identification of larger sets of the most massive halos, paving the way to tighter testing of the structure formation model.
We thank the referee for suggestions that improved the clarity of the paper. We thank Marco Baldi for providing snapshots of the ΛCDM L-CoDECS N -body simulation. We acknowledge Benedikt Diemer for useful discussions. The graduate-student fellowship of Michele Pizzardo is supported by the Italian Ministry of Education, University and Research (MIUR) under the Departments of Excellence grant L.232/2016. We also acknowledge partial support from the INFN grant InDark. The Smithsonian Institution supports the research of Margaret Geller and Jubee Sohn. This research has made use of NASA's Astrophysics Data System Bibliographic Services.