# Constraining cluster masses from the stacked phase space distribution at large radii

###### Abstract

Velocity dispersions have been employed as a method to measure masses of clusters. To complement this conventional method, we explore the possibility of constraining cluster masses from the stacked phase space distribution of galaxies at larger radii, where infall velocities are expected to have a sensitivity to cluster masses. First, we construct a two component model of the three-dimensional phase space distribution of haloes surrounding clusters up to 50 Mpc from cluster centres based on -body simulations. We find that the three-dimensional phase space distribution shows a clear cluster mass dependence up to the largest scale examined. We then calculate the probability distribution function of pairwise line-of-sight velocities between clusters and haloes by projecting the three-dimensional phase space distribution along the line-of-sight with the effect of the Hubble flow. We find that this projected phase space distribution, which can directly be compared with observations, shows a complex mass dependence due to the interplay between infall velocities and the Hubble flow. Using this model, we estimate the accuracy of dynamical mass measurements from the projected phase space distribution at the transverse distance from cluster centres larger than Mpc. We estimate that, by using spectroscopic galaxies, we can constrain the mean cluster masses with an accuracy of 14.5% if we fully take account of the systematic error coming from the inaccuracy of our model. This can be improved down to 5.7% by improving the accuracy of the model.

###### keywords:

dark matter – galaxies: clusters: general – galaxies: kinematics and dynamics^{†}

^{†}pubyear: 2018

^{†}

^{†}pagerange: Constraining cluster masses from the stacked phase space distribution at large radii–A

## 1 Introduction

In the framework of hierarchical structure formation, galaxy clusters, which are the biggest self-gravitating system in the Universe, represent the current end point of the evolution of primordial density fluctuations. Thus, we can extract information on the initial density perturbation, the growth of structure, and cosmological parameters from observations of galaxy clusters. For instance, we can extract the matter density () and the amplitude of the density flctuation () from the abundance of galaxy clusters (e.g., Rozo et al. 2010). The uncertainty of estimating cluster masses, which are necessary to compare observations with theory involving dark matter, is one of the dominant sources of the uncertainty in such cosmological analyses. To constrain cosmological parameters accurately, we need to estimate masses of galaxy clusters accurately, which is difficult because masses of clusters are dominated by invisible dark matter.

There are several methods to estimate masses of galaxy clusters, including gravitational lensing (e.g., Schneider et al. 1992; Umetsu et al. 2011; Oguri et al. 2012; Newman et al. 2013; Okabe & Smith 2016), X-ray observations (e.g., Sarazin 1988; Vikhlinin et al. 2006), and the Sunyaev-Zel’dovich (SZ) effect (e.g., Sunyaev & Zeldovich 1972; Arnaud et al. 2010; Planck Collaboration et al. 2014). In addition, we can estimate masses of galaxy clusters by using the relative motion of galaxies surrounding galaxy clusters. Because the gravitational potential of galaxy clusters affects the motions of surrounding galaxies, such motions have information on cluster masses. In this paper, we refer to the mass estimated by motions of galaxies around clusters as the dynamical mass. It is of great importance to compare cluster masses derived by these different methods in order to understand systematic errors inherent to the individual methods. Different methods have different systematic errors, which can be inferred and hopefully corrected for by cross-checking the results of the individual methods.

Since the pioneering work by Smith (1936) who estimated the mass of the Virgo cluster by using motions of galaxies, many papers have studied dynamical masses of galaxy clusters (e.g., Busha et al. 2005; Rozo et al. 2015; Farahi et al. 2016). However, there is a room for improvement in several ways. For example, to estimate dynamical masses accurately, we have to understand the dynamical state of galaxies around galaxy clusters, taking account of motions of galaxies inside clusters as well as those of infalling galaxies. The caustic model is a method to extract the infall sequence to galaxy clusters based on a spherical collapse model (Diaferio & Geller 1997), which is extended to an expanding Universe in Stark et al. (2016). However, there is an ambiguity to define the caustic surface in observations due to the projection effect, which has to be calibrated against -body simulations. In most of the previous studies, motions of galaxies within the transverse distance of from centres of galaxy clusters are used to derive dynamical masses. However, these studies mainly focused on velocity dispersions and did not consider the full distribution of line-of-sight velocities, nor their radial dependence. Furthermore, it has been known that cosmological -body simulations, which are very often used to connect velocity measurements around clusters with cluster masses, contain a lot of uncertainty because they do not include baryonic effects, such as radiative cooling, gas pressure, and feedbacks. Such baryonic effects may change the mass profiles and dynamics of galaxies around galaxy clusters, because the efficiency of tidal stripping and dynamical friction can be significantly affected by baryonic effects (e.g., Suto et al. 2017; Peirani et al. 2017).

One important application of dynamical mass measurements is to test General Relativity by comparing cluster masses estimated by gravitational lensing effect with dynamical masses (e.g., Schmidt 2010; Lam et al. 2012; Lam et al. 2013). Because the dynamical mass and the lensing mass contain different information on the spacetime metric, we can test General Relativity by comparing these two masses. However, in some classes of modified gravity, such as gravity (e.g., Nojiri & Odintsov 2011), the deviation from General Relativity is suppressed at the regions close to cluster centres where the gravity is relatively strong. Hence, to test General Relativity efficiency, it is important to study the outer part of cluster mass profiles (Lombriser 2014; Terukina et al. 2015).

In this paper, we study the phase space distribution of galaxies around galaxy clusters on large separations well beyond the conventional radii used for dynamical measurements, up to several tens of Mpc, by using -body simulations, and propose to use the stacked phase space distribution at these large radii to constrain cluster masses. For this purpose, we construct a new analytical model of the phase space distribution around clusters. While our model of the phase space distribution is similar to those shown in Zu & Weinberg (2013) and Lam et al. (2013) in several ways, we explicitly consider the dependence on the cluster mass so that it can be used to estimate the dynamical mass. After we construct a model of the full three-dimensional phase space distribution, we take account the impact of the Hubble flow, which is indistinguishable from peculiar velocities of galaxies on an individual basis. We discuss how we can treat this effect in projecting the phase space model to finally obtain the projected pairwise line-of-sight velocity distribution between galaxies and clusters.

In this paper, we demonstrate how dynamical masses can be estimated by the stacked phase space distribution at large transverse distances beyond , which is highly complementary to traditional methods to estimate dynamical masses from motions of galaxies within the transverse distance of from cluster centres. A possible advantage of using the outer phase space distribution is that it is expected to be much less sensitive to baryonic effects compared with the phase space distribution within Mpc as discussed above. This suggests that the observed phase space distribution may be compared with -body simulations more directly at such large distances from cluster centres.

This paper is organised as follows. In Section 2, we show our model of the three-dimensional phase space distribution around clusters. In Section 3, we project the three-dimensional distribution along the line-of-sight with the effect of the Hubble flow. In Section 4, we discuss how the dynamical mass is measured using our model. Finally, we conclude in Section 5.

## 2 Three-Dimensional Phase Space Distribution

### 2.1 Simulations

We perform four random realizations of cosmological -body simulations.
Each simulation is performed with a TreePM code - (Springel, 2005), which runs from to in a box of comoving on a side with the periodic boundary condition.
The number of dark matter particles is , corresponding to the mass of each particle of .
The gravitational softening length is fixed at comoving .
The initial condition is generated by a code based on the second order Lagrangian perturbation theory to compute the displacement from the regular lattice preinitial configuration developed in Nishimichi
et al. (2009) and parallelised in Valageas &
Nishimichi (2011).
The transfer function at is computed by the linear Boltzman code (Lewis
et al., 2000).
We adopt , , , , and following the WMAP nine year result (Hinshaw
et al., 2013).
To identify haloes^{1}^{1}1Rockstar identifies both haloes and subhaloes from -body simulations. In this paper, we do not distinguish haloes and subhaloes, and collectively call them haloes. in our -body simulations, we use six-dimensional friend of friend (FoF) algorithm implemented in (Behroozi
et al., 2013).

### 2.2 Stacked Phase Space Distribution

We use these simulations to obtain the phase space distribution of dark matter haloes around galaxy clusters. Because we are interested in average features of dynamics of dark matter haloes, we stack simulated haloes around a large number of galaxy clusters to derive accurate three-dimensional phase space distributions. We show an example of stacked phase space distribution from -body simulations in Fig, 1.

We use haloes with masses in this paper. Note that we define as the virial mass of halo. The virial mass is the mass within , which is the radius within which the average density is times the background density at , which is calculated following Bryan & Norman (1998). We also define cluster mass () as the virial mass of galaxy cluster. We define galaxy clusters as haloes with masses . To mimic observations, we remove galaxy clusters from our analysis if there are any other clusters with larger masses within from those clusters. To analyse the cluster mass dependence of the phase space distribution, we divide galaxy clusters into 24 log-equal mass bins. We regard the geometric mean mass of the lower and upper limits of each mass bin as the mass of the bin. Note that galaxy clusters can become haloes when we focus on other galaxy clusters. We use only snapshots at in this paper for simplicity.

To estimate cluster masses from the line-of-sight velocity () histogram, we first construct a model of the three-dimensional phase space distribution, and then project it along the line-of-sight to obtain the histogram of . The three-dimensional phase space distribution shown in Fig. 1 clearly indicates that the phase space distribution is quite complicated. In what follows, we construct a two component model to take account of the complexity of the phase space distribution.

### 2.3 Fitting with a Two Component Model

We construct a model of the three-dimensional phase space distribution of haloes surrounding galaxy clusters based on the stacked phase space distributions obtained from the -body simulations. Because we stack a large number of galaxy clusters without aligning their orientations, the spherical asymmetry of the phase space distributions for individual clusters should be averaged out. Hence, throughout the paper we assume a spherically symmetric phase space distribution.

We divide a velocity of a halo into three orthogonal components, the radial velocity () and two tangential velocities () as shown in Fig. 2. At this point we only consider pairwise peculiar velocities between clusters and haloes and do not consider the Hubble flow. Since one of the two components of the tangential velocities do not contribute to , we ignore in this paper, and denote as .

Under the assumption of the spherically symmetric phase space distribution, we can describe the probability distribution function (PDF) of the phase space as

(1) |

We then assume that the PDF of the phase space distribution can be divided into two components, the infall component and the splashback component. The infall component corresponds to haloes that are now falling into galaxy clusters, and the splashback component corresponds to haloes that are on their first (or more) orbit after falling into galaxy clusters. Such two component model is also proposed in Zu & Weinberg (2013), although they consider a virial component for which the average value of the radial distribution is fixed to zero. In contrast, in our splashback component we allow a non-zero average radial velocity. In the two component model, equation (1) is described as

(2) |

where both and are properly normalised, and thus denotes the fraction of the splashback component at given . For simplicity, we also assume that there is no correlation between radial and tangential velocities in both of two components. With this assumption, we can rewrite equation (2) as

(3) |

We discuss each component one by one in what follows.

#### 2.3.1 Radial Velocity Distribution

We start with the distribution of the radial velocities. We find from the simulations that the radial velocity distributions at large radii, where the distributions are dominated by the infall component, have non-negligible skewness and kurtosis, as was already shown in Scoccimarro (2004). To incorporate the skewness and kurtosis, we adopt the Johnson’s SU-distribution (Johnson, 1949) as the model function for the radial velocity distribution of the infall component,

(4) |

where

(5) |

The Johnson’s SU-distribution has four free parameters, and ths is sufficiently flexible to model the phase space distributions from -body simulations including the skewness and kurtosis. We note that these four parameters, , and are functions of the radius .

As for the splashback component, on the other hand, we adopt the Gaussian distribution

(6) |

since we do not find a clear evidence for higher order cumulants. Again, we note that the free parameters, and , are determined as functions of .

The model function for the radial velocity distribution is

(7) |

For a given radial and mass bin, we fit the histogram of radial velocities obtained from our -body simulations with equation (7) that contains seven parameters.

Fig. 3 shows examples of radial velocity distributions. Here the error bars are Poisson errors from the number of haloes in each bin. The Figure indicates that our model of the radial velocity distribution is in good agreement with the histogram from our -body simulations. Since the splashback component is negligibly small at large , in fitting we always fix at larger than .

#### 2.3.2 Tangential Velocity Distribution

Next we model the tangential velocities. We adopt the Johnson’s SU-distribution (Johnson, 1949) again as a model function for the tangential velocity distribution of the infall component, and the Gaussian distribution for the splashback component. Since we assume the spherically symmetric phase space distribution, the tangential velocity distribution must be symmetric aruond . Hence, for the tangential velocity distribution, we fix , and in the Johnson’s SU-distribution and for the Gaussian distribution to ensure the symmetry

(8) |

where

(9) |

(10) |

Again, we note that , , and are functions of .

The model function for the tangential velocity distribution is

(11) |

In equation (11), is always fixed to the value obtained from the analysis of the radial velocity distribution, otherwise it is usually very difficult to separate the two components from the tangential velocities alone. Hence, equation (11) contains additional three free parameters.

Fig. 4 shows examples of the tangential velocity distributions. We can see that our model function for the tangential velocity distribution is also in good agreement with the histogram from our -body simulations.

### 2.4 Radial and Cluster Mass Dependences

By fitting histograms of radial and tangential velocity distributions from -body simulations with equations (7) and (11), we can obtain 10 parameters in these equations for a given and . In order to construct a model of the three-dimensional phase space distribution as a smooth function of and , we fit these ten parameters with analytic functions. Because the four parameters of Johnson’s SU-distribution are highly degenerate, we impose some additional assumptions.

First, we assume that , , and depend only on and have no dependence. We fit , , and with analytic functions of chosen to match the shape obtained from the fitting in Section 2.3. For , , and , we choose the model functions as

(12) |

(13) |

where and are free parameters. After calculating , , and , we refit the histograms from -body simulations with equations (7) and (11) by replacing these parameters to the smooth functions derived above. Next, we fit the remaining parameters with analytic functions of for each galaxy cluster mass bin. The model functions are

(14) |

(15) |

(16) |

(17) |

where , , and are free parameters, and runs over and , and runs over , , and . We cannot determine parameters of the splashback component where there is no splashback component, i.e. . Hence, for , , , and , we need to set the upper limit of the fitting range. We set the upper limit as the minimum radius where derived by fitting radial velocity distributions with equation (14) become smaller than 0.1. We also set the lower limit of the fitting range for computational reasons. We set the lower limit as the half of the upper limit radius. Finally, we fit these parameters with one common functional form for the mass dependence

(18) |

where is a pivot mass, are free parameters, and runs over , , , , , , and . We adopt the pivot mass .

Now, we can calculate parameters introduced in Section 2.3 as a smooth functions of and . Hence, by using these parameters, we can compute the dependence of the average of , , the standard deviations of , and , on and as

(19) |

(20) |

(21) |

Fig. 5 show the dependence of , , and on and . We can see that more massive clusters have lower and higher , and , even at larger than .

As we will discuss later, one of the key parameters that constrain cluster masses is the radial velocity , whose mass dependence can be understood by a simple analytic argument. First, when haloes are sufficiently close to clusters, based on the energy conservation, we expect that the radial (i.e., infall) velocity scales as

(22) |

Next, as shown in Fisher (1995) and Agrawal et al. (2017), at sufficiently large the radial velocity is derived by using the linear perturbation theory as

(23) |

where is the logarithmic growth rate with being the linear growth factor, is the linear bias factor (e.g., Bhattacharya et al. 2011), is the matter power spectrum, and is the spherical Bessel function of the first order. This linear perturbation theory result suggests that the radial velocity should scale as

(24) |

Fig. 6 shows the cluster mass dependence of the average of , . We determine the normalization of equations (22) and (24) so as to fit the points. We can see that at a small radius the cluster mass dependence of follows equation (22), whereas at a large radius it follows equation (24). The success of this simple analytic argument implies that we may be able to construct a fully analytic model of the average infall velocity that smoothly connects equations (23) and (24), which we leave for future work.

## 3 Projected Phase Space Distribution

### 3.1 Effect of the Hubble Flow

Next we derive the PDF of the line-of-sight velocity that can directly be compared with observations. We can do so by projecting the three-dimensional phase space distribution derived in the previous section along the line-of-sight with the effect of the Hubble flow. In Section 2.3 and Section 2.4, we obtain equation (1) as a smooth function of and . We can describe the PDF of by using equation (1) as

(25) |

where is the transverse distance from the cluster centre, is the line-of-sight distance from the cluster centre, is the number density profile of haloes, is the Dirac’s delta function, and is the three-dimensional phase space distribution defined in equation (1). We note that is defined as

(26) |

and is a normalization factor defined as

(27) |

where is the range of we calculate. Because of practical reasons, we set in this paper. The choice of is unimportant because we can also apply this cut to observational data. Also, is defined as

(28) |

where is the Hubble parameter, and corresponds to the angle between the line from the cluster centre and the line-of-sight, as shown in Fig. 7.

In what follows, we set the integration range of and as , because we find that there is almost no probability out of this range in the three-dimensional phase space distribution at any radii of our interest. We also set the integration range of as , because for our choice of , this integration range is sufficiently large.

Since we construct the three-dimensional phase space distribution () only at larger than the lower limit we set in Section 2.4, we cannot describe the PDF of at less than the lower limit. Since in this paper we focus on the PDF of at radii larger than 2 Mpc that is larger than the lower limit, this lower limit does not affect our analysis.

In the analysis of observed PDF of , one should adopt that is estimated directly from the observation. However, in comparison with the -body simulation results we use directly measured in the -body simulations. We now check the accuracy of our model calculation of PDFs of . We do so by comparing the PDF of computed from equation (25) with the PDF directly measured from our -body simulations. In order to allow accurate comparisons, here we estimate the error bars of directly from our -body simulations, as described as follows. First, we divide our -body simulations into thirty-two sub-regions. We then measure in each sub-regions and calculate the sample variances of for each sub-regions by comparing from the variance of values among the thirty two sub-regions. We find that errors estimated by this method are similar to Poisson errors, but are slightly larger than Poisson errors typically by 16%. This small difference of errors may arise from the contribution from the large-scale structure of the Universe. For simplicity, we also ignore the covariance between different bins, as we do not find any significant covariance between different bins in our sub-region analysis. In Fig. 8, we compare calculated in equation (25) with directly measured from our -body simulations. We obtain relatively good , which implies that our model is reasonably accurate and the estimated errors from the sub-regions are also reasonable.

### 3.2 Cluster Mass Dependence

Fig. 9 shows the cluster mass dependence of . We can see a significant difference of between different even at larger than . The cluster mass dependence of originates from the number density profile and from the three-dimensional phase space distribution . Since our main interest is the cluster mass dependence purely coming from the latter, we check the cluster mass dependence of for a fixed but including the cluster mass dependence of . The result shown in Fig. 10 indicates that the cluster mass dependence is similar to that shown in Fig. 9. This confirms that the mass dependence from Fig. 9 mostly comes from and we have a sufficiently large cluster mass dependence of even if we fix to the observed radial number density profile of galaxies.

### 3.3 Origin of Cluster Mass Dependence of the PDF of

In Figs. 9 and 10, we show that the cluster mass distribution of is quite complicated despite the fact that the three-dimensional phase space distribution shows a rather simple cluster mass dependence (Fig. 5). Here we explore how such complicated cluster mass dependence appears. To do so, we check the contribution from different line-of-sight position to . Following equation (25), we define as

(29) |

where is a function of and as shown in equation (26).

In Fig. 11, we show for four sets of parameters for a fiducial cluster mass. We can see that high segments contribute to the tail of and low segments contribute to the peak. These are explained as follow. At high segments, the Hubble flow significantly increases , whereas at low segments, the contribution of the Hubble flow is small so that we observe peculiar velocities directly. We then check the dependence of each part of on cluster masses, focusing on at .

First, we focus on a segment with . In the left panel of Fig. 12, we compare for the fiducial mass at and . We confirm that at low , contributes to the peak. In the left panel of Fig. 12, we also compare for fiducial mass and lower mass at . We can see that for the higher cluster mass has a larger width than that for the lower cluster mass, because of larger for higher mass clusters as shown in Fig. 5. Note that at segments, does not contribute to as shown in Section 3.1.

Next we focus on a segment with . In the right panel of Fig. 12, we show a plot similar to the left panel of Fig. 12, but we change from to . In this segment, , , and the Hubble flow contributes to , and this segment contributes to the edge of the peak of . While the effect of the Hubble flow is significant, also makes non-negligible contribution to the PDF. Since clusters with higher masses have lower (see Fig. 5), they shift to the lower mean more significantly, which also results in the wider width of the PDF.

## 4 Prospect for constraining cluster masses

### 4.1 Accuracy of Mass Estimation

Our goal is to estimate masses of galaxy clusters by using for a given . We assume that we obtain observed by stacking galaxies around many galaxy clusters and estimate the mean mass of galaxy clusters by using in the range . We estimate the mean masses of galaxy clusters from observed as follows. First, we calculate by using given and . We compute for different cluster masses, and calculate by comparing observed with our model . We determine the mean mass of galaxy clusters as the mass of our model that results in the minimum .

The can be described as

(30) |

where is the difference between directly measured from our -body simulations and our model purely caused by statistical fluctuations, is that caused by our model inaccuracy, is the error of directly measured from our -body simulations (see Section 3.1). When we compare directly measured from our -body simulations with our model with different cluster masses, there exists the difference arising from the difference between assumed and true cluster masses. We denote this difference as . As there is no correlation between and or between and , we rewrite equation (30) as

(31) |

where,

(32) |

Note that must obey chi-square distribution with . Assuming that statistical errors derived from our -body simulations are dominated by the Poisson errors, the dispersion depends on the number of haloes in the -th bin as

(33) |

where is the number of haloes at the -th bin.

Fig. 13 shows calculated by comparing directly measured from our -body simulations for a fixed cluster mass range with our model for different cluster masses. We fix to the value that directly measured from our -body simulations, and calculate for different cluster masses. In total we use 431,925 haloes to derive directly from our -body simulations.

Fig. 13 indicates that we can indeed constrain cluster masses with our approach, although there is an offset between input and best-fitting cluster masses, which originates from the inaccuracy of our model. We fit the derived curve with the following function

(34) |

where , , and are free parameters. The best fit line is shown in Fig. 13. We obtain , , and for Fig. 13. The parameter represents the systematic error of galaxy cluster mass estimation caused by the inaccuracy of our model, and represents the statistical accuracy of cluster mass estimation. As we see in equations (31) and (33), depends on the number of haloes we use as

(35) |

where