Probing DE models with EVS of galaxy clusters from the DEUS-FUR simulations

# Probing dark energy models with extreme pairwise velocities of galaxy clusters from the DEUS-FUR simulations

Vincent R. Bouillot, Jean-Michel Alimi, Pier-Stefano Corasaniti and Yann Rasera
Centre for Astrophysics, Cosmology & Gravitation, Department of Mathematics & Applied Mathematics,
University of Cape Town, Cape Town 7701, South Africa
CNRS Laboratoire Univers et Théories (LUTh), UMR 8102 CNRS, Observatoire de Paris,
Université Paris Diderot, 5 Place Jules Janssen, 92190 Meudon, France
Institut d’Astrophysique de Paris, UMR 7095 CNRS, Université Pierre et Marie Curie, 98bis Blvd Arago, 75014 Paris, France
vincent.bouillot@uct.ac.zajean-michel.alimi@obspm.frpier-stefano.corasaniti@obspm.fryann.rasera@obspm.fr
July 1, 2019
###### Abstract

Observations of colliding galaxy clusters with high relative velocity probe the tail of the halo pairwise velocity distribution with the potential of providing a powerful test of cosmology. As an example it has been argued that the discovery of the Bullet Cluster challenges standard CDM model predictions. Halo catalogs from N-body simulations have been used to estimate the probability of Bullet-like clusters. However, due to simulation volume effects previous studies had to rely on a Gaussian extrapolation of the pairwise velocity distribution to high velocities. Here, we perform a detail analysis using the halo catalogs from the Dark Energy Universe Simulation Full Universe Runs (DEUS-FUR), which enables us to resolve the high-velocity tail of the distribution and study its dependence on the halo mass definition, redshift and cosmology. Building upon these results we estimate the probability of Bullet-like systems in the framework of Extreme Value Statistics. We show that the tail of extreme pairwise velocities significantly deviates from that of a Gaussian, moreover it carries an imprint of the underlying cosmology. We find the Bullet Cluster probability to be two orders of magnitude larger than previous estimates, thus easing the tension with the CDM model. Finally, the comparison of the inferred probabilities for the different DEUS-FUR cosmologies suggests that observations of extreme interacting clusters can provide constraints on dark energy models complementary to standard cosmological tests.

###### keywords:
cosmology: theory — galaxies: clusters: general — methods: numerical — methods: statistical

## 1 Introduction

The advent of large scale cluster survey programs (Planck Collaboration, 2011; Reichardt et al., 2012; Pierre et al., 2012; Menanteau et al., 2013) will soon provide large catalogs of galaxy clusters to test the standard cosmological scenario based on the Cold Dark Matter paradigm with Cosmological Constant (CDM).

Galaxy clusters are the largest observable structures in the universe which reside inside massive Dark Matter (DM) halos. These are gravitationally bound objects that result from the gravitational collapse of small matter density fluctuations which were present in the early universe. Because of this the number density of DM halos carries complementary information on both the statistics of the primordial matter density field and the growth rate of structures. Measurements of the abundance of cluster of galaxies aim to probe such features, but their effectiveness to constrain cosmological models depends upon the precise understanding of the survey selection function (Pacaud et al., 2006; Pierre et al., 2012; Clerc et al., 2006) as well as the availability of accurate measurements of cluster masses (see e.g. Majumdar & Mohr, 2004; Lima & Hu, 2004; Cunha, 2009). However, measuring the mass of several hundreds of clusters is a challenging task that requires costly follow-up observations of individual objects, thus inevitably spanning these projects over long periods of time. While completing these programs, cosmological information can still be inferred from reduced datasets consisting of the most massive objects. These probe the high-mass end of the halo mass function and thus have the potential to rule out entire classes of cosmological models. The advantage is that such systems, being also the most luminous, are easier to detect and being limited in number makes their follow-up observations more readily accessible.

In recent years this complementary approach to cluster cosmology has received lots of attention due to the discovery of very massive clusters at high-redshift (Jee et al., 2009; Rosati et al., 2009; Foley et al., 2011; Menanteau et al., 2012; Stalder et al., 2013). These detections have lead several authors to question the basic assumptions of the concordance CDM model (Jimenez & Verde, 2009; Hoyle, Jimenez & Verde, 2011; Holz & Perlmutter, 2012). However, estimating the probability that such extreme clusters occur is far from being trivial. As pointed out by Hotchkiss (2011) assessing the rareness of clusters in terms of the probability of observing at least one cluster of mass larger than that observed and/or at a higher redshift can lead to biased conclusions.

A natural framework to address these questions is given by the Extreme Value Statistics (EVS). This can be used to predict the probability distribution of the mass of the most massive halo in the sample from prior knowledge of the halo mass function in a given cosmological model. This has been the subject of several studies (see e.g. Davis et al., 2011; Waizmann, Ettori & Moscardini, 2011; Harrison & Coles, 2011) which have focused on the mass as measure of cluster extremeness. However, a careful analysis of the most massive systems so far observed suggests that there are other characteristics that are indicative of the extremeness of these objects. For instance, 1E0657-56 (Markevitch et al., 2002, 2006), MACS J0025.4-1222 (Bradac et al., 2008), ACT-CL J0102-4915 (Menanteau et al., 2012) and AS1063 (Gomez et al., 2012) are merging clusters with high relative velocities. Among these 1E0657-56 is one of the most well documented. Known as the “Bullet Cluster”, it is composed of two clusters which have undergone a nearly head on collision. The main cluster has a mass h M, while the smaller one has mass h M (Clowe, Gonzalez & Markevitch, 2004). The system is located at (Clowe et al., 2006; Bradac et al., 2006) and the clusters are separated by a distance of h Mpc. X-ray observations have shown that the collision has ripped off the clusters of their intra-cluster gas which is trapped in a shock with Mach number close to . This implies a velocity of the shock front of km s (Markevitch et al., 2006), which has been usually interpreted as being also the relative velocity of the colliding clusters. Under this hypothesis Hayashi & White (2006) analysed the velocities of the sub-halo distribution from the Millenium Simulation (Springel et al., 2005) and concluded that the existence of the Bullet Cluster is consistent with the standard CDM cosmology. However, due to the limited volume of Millenium Simulation their result do not rely on direct measurement but rather on extrapolating the sub-halo velocity probability distribution to host halos with mass h M. Furthermore, the relative velocity of the colliding clusters may well be different from that of the gas. As shown by Milosavljevic et al. (2007) using 2-D hydrodynamical simulations this can be up to smaller (see also Springel & Farrar, 2007).

To date the most detailed study of the Bullet Cluster has been performed by Mastropietro & Burkert (2008) who have used 3-D non-cosmological hydrodynamical simulations to determine the initial physical configurations of the colliding clusters resulting in a system whose properties reproduce those observed in the Bullet Cluster. These include the displacement between the X-ray peaks and the mass distribution, the morphology of the shock velocity, the surface brightness and the projected temperature profiles across the shock. Mastropietro & Burkert (2008) have shown that the colliding halos must have an initial velocity v km s with a mass ratio of , an initial separation of Mpc (implying an initial redshift of ) and a collision angle close to 0 (i.e. head on collision). The identification of these parameters has provided criteria to select Bullet-like halo pairs in cosmological simulations, a crucial step to estimate the probability of finding the Bullet Cluster in a given cosmological setup.

Lee & Komatsu (2010) analysed the halo catalog from the MICE simulations (Crocce et al., 2010) of CDM cosmology at and to infer the pairwise velocities probability distribution for different mass ratios, distance separation and relative velocity. Quite remarkably they found that none of the analysed catalogues contains a system with parameters corresponding to that of the Bullet Cluster. Nevertheless, by fitting the probability density distribution to a Gaussian, they were able to extrapolate the probability to high relative velocities. They found the rate of occurrence of Bullet Cluster-like systems to be and at and respectively. Thompson & Nagamine (2012) performed a similar analysis of a set of CDM simulations with DM mass resolution varying from 9 to 5.7 h M and box sizes ranging from h Mpc to h Gpc. By extrapolating the cumulative distribution to high relative velocities these authors obtained at for masses M. These probabilities imply that in the observable cosmic volume the existence of the Bullet Cluster pair is either incompatible with the CDM scenario or as argued by Thompson & Nagamine (2012) that the initial conditions determined from the analysis of non-cosmological hydrodynamical simulations by Mastropietro & Burkert (2008) have to be revised to much lower values of v (see e.g. Lage & Farrar, 2013, for a recent study).

A critical point is that all these analyses rely on extrapolating the tail of the pairwise velocity probability distribution to high-velocities. This is a direct consequence of the limited volumes of the numerical simulations from which these results have been derived. Furthermore, the probability of observing the Bullet Cluster has been directly estimated from the tail of the probability density distribution fitted to a Gaussian which may suffer of potential biases especially if the tail of the distribution is non-Gaussian.

In the work presented here we improve these studies in several ways. We use the catalog of halos from the Dark Energy Universe Simulation - Full Universe Runs (DEUS-FUR) which cover the entire observable cosmic volume with a box-length of h Gpc and dark matter particles (Alimi et al., 2012; Rasera et al., 2014). These simulations provide an unprecedented large statistical sample to test the rareness of halo properties for different cosmological models. The large simulation volume allows us to perform a detailed analysis of the pairwise velocity especially in the high-velocity tail. Building upon this study we use the Extreme Value Statistics (EVS) to quantify the probability of observing the Bullet Cluster. This enables us for the first time to perform a cosmological model comparison of the Bullet Cluster observation.

The paper is organized as follows. In Section 2, we introduce the N-body simulation data, the halo finder algorithm and the halo pair selection method. In Section 3, we described the dependence of the pairwise velocity function on the halo finder parameter and discuss the physical implications, while in Section 4 we study the redshift evolution and cosmology dependence. In Section 5 we introduce the Extreme Value Statistics and present the results of its application to the Bullet Cluster. Finally in Section 6 we discuss our conclusions.

## 2 N-Body Simulation Dataset

### 2.1 DEUS Full Universe Runs

We use numerical data issued from the DEUS-FUR project. This consists of three N-body simulations with dark matter particles and box size of (21000 h Mpc) enclosing the observable volume of a flat CDM cosmology and two dark energy models with different expansion histories. The simulations have been realized using the application AMADEUS - A Multi-purpose Application for Dark Energy Universe Simulation - expressly developed for the DEUS-FUR project (Alimi et al., 2012). This includes the generator of Gaussian initial conditions for which we use an optimized version of the code MPGRAFIC (Prunet et al., 2008), the N-body solver which is a version of the RAMSES code (Teyssier, 2002) specifically improved to run on a large number of cores ( 40000 MPI tasks in production mode) and an optimized Friend-of-Friend halo finder (Roy, Bouillot & Rasera, 2014). RAMSES solves the Vlasov-Poisson equations using an Adaptive Mesh Refinement (AMR) Particle Mesh method with the Poisson equation solved with a multigrid technique (Guillet & Teyssier, 2011). A detailed description of the algorithms, optimization schemes and the computing challenges involved with the realization of DEUS-FUR is given in Alimi et al. (2012). All simulations share the same phase of the initial conditions. The coarse grid of the AMR hierarchy contains resolution elements and is allowed to refined six times, reaching a formal spatial resolution of 40 h kpc, while the particle mass resolution for the different models is h M. Such resolution roughly corresponds to the size and mass of the Milky Way.

The simulated cosmologies consist of a flat CDM model best-fit to the WMAP-7yr data (CDM-W7, Spergel et al., 2007), a quintessence model with Ratra-Peebles potential (RPCDM, Ratra & Peebles, 1988) and a “phantom” dark energy model with constant equation of state (CDM, Caldwell, Kamionkowski & Weinberg, 2003).

The model parameters have been calibrated to fit the Cosmic Microwave Background anisotropy power spectra from WMAP 7-year observations (Spergel et al., 2007) and luminosity distances measurements to Supernova Type Ia from the UNION dataset (Kowalski et al., 2008). In particular, the mean cosmic matter density, , has been chosen within the marginalized contour in the plane along the degeneracy line of the SN Ia data (see left panel in Fig. 1); while the values of the root-mean-square fluctuation amplitude of the density contrast at h Mpc, , has been chosen within the marginalized confidence contours in the plane nearly parallel to the degeneracy line of the CMB data (see right panel in Fig. 1). This particular choice is motivated by the fact that through the analysis of the DEUS-FUR simulations we aim to test whether observables of the non-linear clustering of matter can break the degeneracies affecting current cosmological parameter constraints. A summary of the cosmological model parameter values and the simulation characteristics are reported in Table 1.

### 2.2 Halo Finder and Halo Pair Selection

The detection of halos in the DEUS-FUR simulations is performed with a highly-scalable parallelized version of the Friend-of-Friend halo finder algorithm (Roy, Bouillot & Rasera, 2014) implemented in the AMADEUS application. This algorithm detects halos as groups of particles characterized by an interparticle distance smaller than a given linking length (in units of the mean interparticle distance) or percolation parameter .

In the study presented in Section 3 we consider halos detected with and respectively. In order to limit the effect of numerical artifacts due to low number of particles we only consider halos with particles, which corresponds to halo masses h M.

In Table 2 we report the total number of FoF(b=0.2) halos detected in the DEUS-FUR simulations at different redshifts. These are vast halo catalogues for which the selection of halo pairs and the calculation of relevant quantities poses a challenging computational problem. In fact, the complexity of a standard brute force computation of the relative velocities for all pair separations grows with the square of the number of halos. This leads to a prohibitive computational time as soon as the number of halos exceeds . However, we are interested only in halo pairs characterized by a small distance separation such as the Bullet Cluster initial configuration found in Mastropietro & Burkert (2008). Therefore, we can significantly reduce the number of computations by using an octree space decomposition which enables the computation of pairs up to a maximum distance . To be conservative we set , where is the radius of the most massive halo in the catalogues. This gives a maximal distance separation h Mpc. Thanks to such a space decomposition we compute all relevant quantities for all pairs within the pruning radius, thus avoiding the most time consuming long-range computations. The implementation of this algorithm allows us to compute the relative velocities of about 100 million pairs in less than 2 minutes on a 64 core (Intel Xenon CPU X75502.00Ghz) local machine. The number of pairs within 15 h Mpc for the three cosmological catalogs at three redshifts of interest is given in Tabel 2. As we can see at the number of pairs is a significant fraction of the total number of halos detected in the simulations which is indicative of the high level of clustering of such objects compared to that of the average density field with a mean inter-halo separation of 40-70 h Mpc.

## 3 Pairwise Velocities and Friend-of-Friend Halos

In this section we study the dependence of the pairwise velocity function on the halo definition specified by the value of the percolation parameter . It is usually understood that for a given value of the FoF algorithm selects halos whose boundary has approximatively a fixed isodensity surface. For instance, the local surface overdensity with respect to the mean matter density of two particles within a sphere of radius is (Summers, Davis & Evrard, 1995; Audit, Teyssier & Alimi, 1998). In the case this gives which assuming an isothermal density profile, , corresponds to an enclosed overdensity with respect to the mean matter density of . This is close to the value of the virial overdensity predicted by the spherical collapse model in the Einstein-De Sitter cosmology. That is why the percolation parameter is commonly set to . However, More et al. (2011) have shown that the boundary of FoF halos is not associated to a unique local surface overdensity but is distributed around a characteristic value (for this is ). In addition the profile is not isothermal. As a consequence the enclosed overdensity is much higher than 180 and is found to vary in the range to at . It is important to keep this in mind when comparing the properties of N-body halos to observations. In fact, the mass of clusters is usually estimated in terms of the enclosed spherical overdensity with respect to the critical cosmic density that, depending on redshift and cosmology, may correspond to different percolation parameter values.

In Fig. 2 we plot the pairwise velocity distribution for halo pairs with a separation hMpc from the DEUS-FUR CDM-W7 simulation at detected with and respectively (corresponding to a typical variation of the enclosed overdensity by a factor of 8). We can see that for increasing values of the peak of the distribution increases while its width decreases. This trend is a direct consequence of the halo definition corresponding to the different values of . Such a difference implies a change of the halo mass function that manifests in the pairwise velocity distribution.

We may gain a better insight from Fig. 3 where we plot isocontours of the multivariate pairwise velocity distribution as function of the relative velocity and the arithmetic mean mass M of halo pairs detected with and (panels from left to right) and distance separation and h Mpc (panels from top to bottom) respectively.

In the case of pairs with h Mpc (top panel) the size of the isocontours increases for decreasing values of . This is because FoF(b=0.1) not only detects individual groups of particles but also sub-structures within halos detected by FoF(b=0.2), thus the former selects a greater number of halo pairs. This population of pairs is distributed in a lobe structure (Lobe 1) which extends from small average masses with low relative velocities to large average masses with moderate relative velocities. The spread of Lobe 1 increases for distance separations h Mpc, simply because FoF detects a greater number of halo pairs. Notice that the tails of the distribution does not exceeds km s. This indicates that halo pairs in the tail of Lobe 1 correspond to low velocity massive mergers. This could be the case of some observed systems such as MACS J0025.4-1222 (Bradac et al., 2008) or ACT-CL J0102-4915 (Menanteau et al., 2012) which are examples of pairs of massive interacting clusters with Bullet-like baryonic features.

For h Mpc we can see the emergence of a second lobe (Lobe 2) that causes the multivariate probability density distribution to become bimodal. This second lobe consists of small average mass pairs with velocities extending up to km s. This is the population that seems to better correspond to the characteristics of 1E0657-56, the Bullet Custer. Notice that the velocity tail of Lobe 2 shifts to larger values for decreasing values of . This is because FoF detects a greater number of satellite halos that translates into an increase of the number of small mass high-velocity pairs as decreases.

Although the differences of the pairwise velocity probability density function between the case and might look minor in the high-velocity tail, these may have an important impact on the evaluation of the probability of high-velocity colliding clusters. Hence, when comparing to observations an important point concern the choice of the percolation parameter which has to be as consistent as possible with the observational mass definition. For instance in Mastropietro & Burkert (2008) the colliding halos are spherical objects with a Navarro-Franck-White (NFW) profile (Navarro et al., 1996, 1997) characterized by a concentration parameter for which the mass is defined in terms of the virial mass , which is the mass contained in a spherical region of radius enclosing an overdensity with respect to the critical density of the universe at the redshift of the halo. These halos are sampled with  particles, thus following More et al. (2011) an enclosed overdensity of at roughly corresponds to applying FoF with a linking-length .

In the following, we therefore adopt a linking-length and limit our analysis to halo pairs with distance separation h Mpc which is the minimum distance for which two massive merging halos with virial radius h Mpc can be detected as distinct objects.

## 4 Statistics of pairwise velocities

### 4.1 Redshift evolution and cosmology dependence

We now focus on evaluating the dependence of the pairwise velocity function on redshift and cosmology.

In Fig. 4 we plot the probability density function associated with for halo pairs from the DEUS-FUR CDM-W7 simulation with distance separation h Mpc at and respectively. First, we may notice that the amplitude remains constant with redshift, this is the normalization of each curve is different as the number of halo pairs grows with cosmic time. Another effect concerns the tail of the velocity function which tends to slightly increase towards higher velocities from to . For instance we find the maximal relative velocities to slowly evolve with redshift with and km s at and respectively. Because of this, we can expect that in a given cosmological model the probability (defined as the ratio of the velocity function to the total number of pairs) of finding a halo pair with a large relative velocity at redshift and is not significantly different.

The advantage of using the halo catalogues from the DEUS-FUR simulations can be appreciated from Fig. 4. Despite statistical scatter, even at , the high-velocity tail of the pairwise velocity function is resolved to km s. In this range previous analyses had to strongly rely on extrapolation from fitting functions calibrated to lower relative velocities. For instance, Thompson & Nagamine (2012) approximated the cumulative pairwise velocity distribution of FoF(b=0.2) halo pairs in CDM model simulations at with a quadratic fit which we plot in the top panel of Fig. 5 against the DEUS-FUR CDM-W7 results. Quite remarkably we can see that it provides a very good approximation up to intermediate velocities.111The quadratic fit from Thompson & Nagamine (2012) has been derived from the analysis of FoF(b=0.2) halo pairs in simulations with mass and spatial resolution different from those of the DEUS-FUR simulations. Such differences may affect the halo mass function and introduce a systematic bias in the cumulative pairwise velocity distribution. However, for halo pairs with average masses this effect remains negligible. In contrast, large discrepancies occurs in the high-velocity tail as can be seen from the bottom panel of Fig. 5. In the same figure we also plot the cumulative distribution from our reference catalogue of halo pairs detected with FoF(b=0.15). As expected from the discussion in Section 3 this is characterised by both a greater number of halo pairs and a longer tail at high-velocity compared to .

Let us now turn to the cosmological dependence of the velocity function. In Fig. 6 we plot at for a distance separation h Mpc for RPCDM, CDM-W7 and CDM. The first noticeable difference is the overall amplitude of the various curves which is essentially caused by the difference of the mass function of the DEUS-FUR simulated cosmologies. We can also notice that in the high-velocity interval the velocity functions have slightly different slopes, with the CDM and CDM-W7 models showing a heavier tail than RPCDM. This is indicative of the fact that halo pairs with extreme relative velocities are a sensitive probe of the underlying cosmological model. Understanding the mechanisms responsible for this dependence requires a physical analysis that is beyond the scope of this paper.

Even in the lack of such study we can have an idea of the dominant cosmological parameter dependence by evaluating the average pairwise velocity . This can be directly inferred from the pairwise probability distribution function obtained from and confronted with predictions from the following relation derived in the context of stable clustering (Juszkiewicz, Springel & Durrer, 1999; Caldwell et al., 2001):

 −¯v12(x,a)Hr=a3[1+ξ(x,a)]∂¯ξ(x,a)∂a , (1)

where is the two-point correlation function of the density field, is the expansion factor, is the proper separation, is the Hubble rate and

 ¯ξ(x,a)=3x3∫x0ξ(y,a)y2dy (2)

is the two-point correlation function averaged over a sphere of radius . Evaluating Eqs. (1) and (2) using the correlation function of the density field of each of the DEUS-FUR simulations, we obtain the following average pairwise velocities: km s for the RPCDM model, km s for the CDM-W7 and km s for the CDM model. These values are within of those directly estimated from the numerical data. Their variation is essentially due to the different values of of the DEUS-FUR cosmologies. This can be understood by taking the ratio of the average pairwise velocity given by Eq. (1) for a given with respect to a reference one :

 ¯v12 (σ8)¯v12 (σ8,ref)=σ8σ8,ref1+ξ(r)1+σ8σ8,refξ(r) . (3)

Assuming that the shape of the power spectrum does not change over the range of scales where the stable clustering regime occurs (Smith et al., 2003) and considering as a reference case the CDM-W7 values of and , we recover from Eq. 3 the estimated values of the average pairwise velocity of both RPCDM and CDM to better than a few per cent.

In Fig. 7 we plot the normalized probability density distribution of the pairwise velocity of the three simulated cosmologies. The dependence described above can be seen here on the fact that the distributions have very similar average and overall amplitude. This is because their respective normalizations encode differences of the mass function of the underlying cosmological models that are mostly due to the different values. On the other hand we can see that the high-velocity tail of the distributions is where the cosmological models differ the most. This indicates that the high-velocity tail carry information not only on , but also on the cosmic matter density and the properties of the Dark Energy which characterize the simulated models.

### 4.2 Bullet-like halo pairs in DEUS-FUR cosmologies

In order to identify extreme halo pairs in the DEUS-FUR simulation it is instructive to consider the redshift and cosmology dependence of the multivariate pairwise velocity distribution (see Appendix B). This is shown in Fig. 8 for halo pairs with distance separation h Mpc at and (panels top to bottom) for DEUS-FUR RPCDM, CDM-W7 and CDM simulations (panels left to right) respectively. As expected from Section 4, for a given cosmological model the range of relative velocities does not vary significantly as function of the redshift compared to interval variation of the average mass of the pairs. Indeed, this is due to the different redshift evolution of the halo mass function compared to that of the pairwise velocity function. Furthermore, we can see that in the RPCDM case the tail of velocity distribution at all redshift remains confined to low velocities. This is not the case of CDM-W7 and CDM for which we have a large number of pairs with average mass larger than M and relative velocity v km s. For these models the bimodality of the multivariate pairwise distribution is more pronounced than in RPCDM. In particular, the halo pairs in the tail of Lobe 1 and Lobe 2 of the multivariate distribution clearly indicate that observations of high-mass moderate-velocity merging clusters and low-mass high-velocity ones can provide powerful cosmological probes.

We are especially interested in extremal halos belonging to the latter category. In Appendix A we summarize the characteristics of the highest velocity pairs and the most massive ones detected in the DEUS-FUR simulations. We find that the properties of such extremal pairs vary with the cosmological model. As an example, from the values quoted in Table 6 we notice that even at the RPCDM simulation has no pairs with relative velocities exceeding km s and an average mass  hM. More generally from Tables 6 and 7 we can infer that the deficiency of high-velocity halo pairs with mass above h M and that of massive pairs with relative velocities above km s tend to disfavor such a cosmological model. The extremal halo pairs in RPCDM are low-velocity massive mergers with large distance separation. In contrast, in CDM the highest velocity pairs all exceed km s, though their average masses remain small.

For comparison, in the top panels of Fig. 8 we plot the average mass and relative velocity of the Bullet Cluster (blue solid lines) as inferred from the analysis of Mastropietro & Burkert (2008). Using this as a reference of the extreme halo pairs we can see that RPCDM has no candidate pairs reproducing the Bullet Cluster characteristics. In the CDM-W7 case the best candidate pair has a main halo with mass h M and a smaller halo with mass h M (corresponding to a mass ratio ) separated by a distance of 8.4 h Mpc, which experience an head on interaction with a relative velocity of 3011 km s. Notice that such an object is absent from the list of extreme halo pairs shown in Appendix since it has neither an extreme velocity nor a very high mass. In the CDM model the best candidate is characterised by a main halo with a mass h M and a lighter halo of mass h M (corresponding to a mass ratio of ) separated by a distance of 8.8 h Mpc and a relative velocity 2839 km s. This candidate is no better than that of CDM. This seems contradictory given the cosmological dependence of the high-velocity tail of the multivariate distribution shown in Fig. 8. However, such discrepancy can be simply a consequence of the specific realization of the simulation run, such as the phase of the initial conditions. Hence, an object by object comparison is no meaningful in assessing the extremeness of the Bullet Cluster. Such estimation can only be performed through a statistical analysis of extreme halo pairs. In the next Section we will discuss the use of these pairwise velocity catalogs to infer the probability of observing the Bullet Cluster.

## 5 Extreme Value Statistics of Pairwise Velocities

### 5.1 Methodology

Extreme Value Statistics, originally pioneered by Fréchet (1927), Fisher & Tippett (1928), Gumbel (1935) and Gnedenko (1943), has been applied to a wide variety of problems to model the probability of extreme events. Here, we briefly review the basic formalism.

Consider a set of independent identically distributed random variates drawn from a cumulative distribution and . It is easy to show that in such case the cumulative distribution function of the maximum of the first observations is given by

 P(Xmax≤x)=FN(x). (4)

This is the so called “exact” EVS approach in which is well known. In the large limit, it is possible to show that the cumulative distribution function of extreme observations tends to the Generalized Extreme Value (GEV) distribution:

 P[μ,σ,ξ](x)=exp{−[1+ξ(x−μσ)]−1/ξ}, (5)

defined for , where is the location parameter, is the scale parameter and is the tail index (or shape parameter), which generalizes the central-limit theorem to extremal subset of data. Depending on the value of , Eq. (5) reduces to three possible functional forms: the Gumbel (or type I) distribution (), the Fréchet (or type II) distribution () and the Reversed Weibull (or type III) distribution ().

Contrary to the exact EVS approach, the use of the Generalized Extreme Value distribution does not require prior knowledge of the underlying cumulative function of the random variates. Instead, it uses these observations to infer the GEV distribution parameters. This is done by classifying the data into blocks of arbitrary size, determining the maxima in each block and inferring the GEV parameters by best fitting the GEV function to the distribution of maxima. A potential disadvantage of this block maxima method is the fact that data need to be sampled. This may cause some loss of valuable rare information or inclusion of non-extremal events.

A complementary approach that is more suited to our purposes consists in using the Generalized Pareto distribution (GPD). This corresponds to Taylor expanding the tail of Eq. (5) to obtain the cumulative distribution function of observing extreme events above a fixed threshold . This reads as

 P[μ,σ,ξ](x)=1−[1+ξ(x−μσ)]−1/ξ, (6)

defined for . Again depending on the value of we have different probabilities of the extreme events. In particular () corresponds to a long (short) tail distribution. Instead the case corresponds to the distribution of events in the tail of a Gaussian distribution. Hence, studies that have extrapolated the probability of the relative velocity of Bullet Cluster in terms of a Gaussian pairwise velocity probability distribution can be seen as a limiting case of the EVS approach described here, with fixed to zero. Notice also that since the tail of the pairwise velocity function depends on the underlying cosmology, we can expect the GPD parameters to carry a strong cosmological dependence.

In the GPD approach the issue of sampling the data is replaced by the problem of choosing a suitable value of the threshold. A high threshold would result in a drastic reduction of the data sample, whereas a low threshold may include non-extremal data and thus bias the results toward a gaussian behaviour. This can be seen in Fig. 9 for a subset of the pairwise velocities in the DEUS-FUR CDM-W7 catalog at where we have classified the halo pairs according to three different velocity thresholds. Using pairs above the lowest threshold would lead to a gaussian biased estimation of the GPD parameters, while using points above the highest threshold provides a too small sample of extremal events to determine the GPD.

Several statistical diagnostics have been considered to estimate a suitable threshold that segregates common events from extreme ones (for a review see e.g. Scarrott & MacDonald, 2012). Here, we focus on the mean residual life method and the threshold stability plot that have been developed in relation with EVS data analysis problems.

The mean residual life method consists in plotting the mean excess, defined as the mean of the exceedances of the data minus the threshold, as function of the threshold itself. An optimal choice is then given by the lowest threshold value for which all higher thresholds give a sequence of mean excesses that is consistent with a straight line (Scarrott & MacDonald, 2012).

A mean residual life plot is shown in Fig. 10 for the halo pairs of the full DEUS-FUR CDM-W7 simulation volume (blue solid line) and three subvolumes of boxlength h Mpc (red dashed line), h Mpc (orange dot-dashed line) and h Mpc (yellow squared solid line) respectively. We can see a characteristic trend with the mean of the exceedances rapidly decreasing at low threshold values, while increasing linearly at intermediate thresholds and then sharply decreasing at large values. The main difference among the various volume catalogs is the interval extent and statistical uncertainty of the linear trend. In the case of the full DEUS-FUR volume this interval has the maximum extent implying a precise selection of the GPD threshold which also guarantee a stability of GPD inferred results. For smaller volumes the interval ranges is much smaller and more uncertain such that it becomes impossible to reliably select a threshold value.

The threshold stability plot is joint diagnostic that consists in plotting the GPD shape and scale parameter values best fitting the data as function of the threshold. Then, an optimal threshold value is chosen such that for higher values the GPD parameters remain stable (Scarrott & MacDonald, 2012). We show such a plot in Fig. 11.

From Fig. 10 and Fig. 11 we can see that the curves are nearly constant straight lines in the threshold range to km s, thus we set km s. This guarantees the stability of the results with respect to the choice of the threshold. Performing a similar analysis for the other DEUS-FUR cosmological simulations we set km s for the RPCDM case and km s for the CDM model respectively.

In Fig. 12 we plot the Generalized Pareto density distributions best fitting the probability density distribution functions of halo pairs with distance separation h Mpc and average pair mass 10 h M from the DEUS-FUR cosmological simulation catalogs at . The selected thresholds for the different cosmological models are indicated by vertical dashed lines, while the dot-dashed lines correspond to the Gaussian tails (with ) for the same threshold and scale parameter values.

The best-fit values and the confidence interval of and have been determined through a Monte Carlo Markov Chain likelihood analyses of the binned numerical data assuming Poisson errors. The inferred values are quoted in Table 3. Notice that in all cases a Gaussian tail () is excluded at more than confidence level. Since the best-fit value of the shape parameter is positive this implies the probability distribution of extreme pairwise velocities is slightly heavy tailed, which increases the probability of finding high relative velocity pairs compared to previous studies that have simply assumed a Gaussian distribution (Lee & Komatsu, 2010; Thompson & Nagamine, 2012).

### 5.2 Application to the Bullet Cluster

Having determined the GPD parameters from each of the DEUS-FUR halo pair catalogs, we are now able to estimate the probability of observing the Bullet Cluster for different DEUS-FUR cosmological models. This is obtained by integrating the probability density functions shown in Fig. 12 from to infinity. This probability has to be interpreted as the rate of occurrence of Bullet Cluster-like systems in comoving space.

In the CDM-W7 case we find , which is two orders of magnitude larger than previous estimates (Lee & Komatsu, 2010; Thompson & Nagamine, 2012). For the RPCDM model we obtain , while we find for the CDM case.

As shown in Section 4 the pairwise velocity distribution carries information on cosmological model parameters such as , as well as and which differentiate the DEUS-FUR cosmologies. The value of these parameters have been selected along the (and ) degeneracy line of the CMB (and SN Ia) data. Henceforth, the Bullet Cluster inferred probabilities can be used to provide us with some qualitative constraints on these class of models. In particular, these suggest that the observation of the Bullet Cluster strongly disfavors Dark Energy models, such as RPCDM, which have an equation of state for which CMB data enforce smaller values to compensate for the greater amplitude of the integrated Sachs-Wolfe (Sachs & Wolfe, 1967) effect on the CMB temperature anisotropy power spectrum (Kunz et al., 2004) while SN Ia data enforce a lower value of to compensate for the shorter luminosity distance. These models are characterized by a lower level of matter clustering with respect to the standard CDM-W7 model. In contrast, the probability of finding the Bullet Cluster increases in the case of Dark Energy models with more negative values of the equation of state for which CMB data enforces larger values while the SN Ia data requires larger values of . These models are characterized by a higher level of matter clustering compared to the CDM-W7 case. Henceforth, it is plausible that the statistical measurements of the rate of occurrence of bullet cluster-like systems which sample the tail of the pairwise velocity distribution have the potential probe Dark Energy and break degeneracy lines of the underlying cosmological parameters.

## 6 Conclusion

We have explored the possibility of testing cosmological models through observations of extreme pairwise velocities of interacting galaxy clusters. To this purpose we have studied the properties of pairwise velocities from the halo catalogs of the DEUS-FUR cosmological simulations. Thanks to the large simulation volume we have been able to resolve the high-velocity tail of the pairwise velocity distribution. We have studied its dependence on the percolation parameter of the FoF halo finder, the distance separation, the redshift evolution and cosmology. We have shown that a particular attention has to be paid to the halo mass definition, since the choice of the percolation parameter especially alter the tail of pairwise velocity function. In the redshift range to the latter show minor evolution, while it significantly varies with cosmology. To have an idea of the cosmological model parameter dependence we have estimated the average pairwise velocity of the DEUS-FUR cosmologies using a model based on stable clustering. From the comparison with the mean value inferred from the DEUS-FUR halo pairs catalogs we have shown that most of the average cosmological dependence is driven by the value of while the tail of the distribution carries information on , and which differentiate the DEUS-FUR cosmologies. As such observations of extreme relative velocities can be used as a different probe to measure the equation of state of Dark Energy and test cosmological models. In particular, the analysis of the multivariate pairwise velocity distribution indicates that observations of low-mass high-velocity interacting clusters (e.g. Bullet Cluster) as well as massive systems with moderate relative velocities are most sensitive to the underlying cosmology.

Focusing on the Bullet Cluster system, we have found a number of halo pairs candidates in CDM-W7 and CDM catalogs respectively, while we have found none in the RPCDM case within the simulated volume of the observable universe. Built upon these results we have quantified the probability of observing the Bullet Cluster in the context of Extreme Value Statistics. To this end we have used Generalized Pareto distributions to model the probability distribution of the pairwise velocities. We find the probability of observing a halo pair with average mass 10 h M, distance separation h Mpc and relative velocity km s to strongly vary across the DEUS-FUR simulated cosmologies with probabilities of , for RPCDM, CDM-W7 and CDM respectively. Thus, we can deduce that the observation of the Bullet Cluster strongly disfavours cosmologies with low value of the (low level of matter clustering). In the case of Dark Energy models calibrated against CMB observations this occurs for since CMB data enforce low values primarily to compensate for the enhanced amplitude of the ISW effect. In contrast, the probability of the Bullet Cluster increases for models with larger and thus a higher level of matter clustering.

The study presented here suggests that observations of extreme interactive clusters sampling the tail of the pairwise velocity distribution can provide complementary information on Dark Energy models potentially capable of breaking standard cosmological parameter degeneracies. A first step in this direction will be the development of an accurate theoretical model of the pairwise velocity distribution density function along the line of the original work by (Sheth, 1996; Diaferio & Sheth, 2001). This will help elucidating the cosmological dependence of the high-velocity tail and provide an estimate of the number of bullet-like systems in the sky that need to be observed to improve current parameter constraints on Dark Energy.

## Acknowledgements

We warmly thank Vincent Reverdy and Irène Balmes for fruitful discussions about EVS and their deep implication in the achievement of the DEUS-FUR project. This work was granted access to HPC resources of TGCC through allocations made by GENCI (Grand Équipement National de Calcul Intensif) in the framework of the “Grand Challenge” DEUS. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 279954. We acknowledge support from the DIM ACAV of the Region Ile de France. V. R. Bouillot is supported financially (AW) by the National Research Foundation of South Africa. Any opinion, findings and conclusions or recommendations expressed in this material are those of the authors and therefore the NRF does not accept any liability in regard thereto.

## References

• Alimi et al. (2010) Alimi, J.-M., Füzfa, A., Boucher, V., et al. 2010, MNRAS, 401, 775
• Alimi et al. (2012) Alimi, J.-M., et al., 2012, IEEE Computer Soc. Press, CA, USA, SC2012, art. 73, arXiv:1206.2838
• Audit, Teyssier & Alimi (1998) Audit E., Teyssier R., Alimi J.-M., 1998, A&A, 333, 779
• Bradac et al. (2006) Bradac. M., et al., 2006, ApJ, 652, 937
• Bradac et al. (2008) Bradac. M., et al., 2008, ApJ, 687, 959
• Caldwell et al. (2001) Caldwell R., Juszkiewicz R., Steinhardt P.J., Bouchet F.R., 2001, ApJ, 547, L93
• Caldwell, Kamionkowski & Weinberg (2003) Caldwell R.R., Kamionkowski M., Weinberg N.N., 2003, PRL, 91, 071301
• Clerc et al. (2006) Clerc, N., Pierre, M., Pacaud, F., Sadibekova, T. 2012, MNRAS 423, 3545
• Clowe, Gonzalez & Markevitch (2004) Clowe D., Gonzalez A., Markevitch M., 2004, ApJ, 604, 596
• Clowe et al. (2006) Clowe D., et al., 2006, ApJ, 648, L109
• Crocce et al. (2010) Crocce M., et al., 2010, MNRAS, 403, 1353
• Cunha (2009) Cunha, C. 2009, PRD 79, 063009
• Davis et al. (2011) Davis O., et al., 2011, MNRAS, 413, 2087
• Diaferio & Sheth (2001) Diaferio A., Sheth R.K., 2001, MNRAS, 322, 901
• Fisher & Tippett (1928) Fisher R. A., Tippett L. H. C., 1928, Proc. Cambridge Philosophical Soc. 24, 180
• Foley et al. (2011) Foley R. J., et al., 2011, ApJ, 731, 86
• Fréchet (1927) Fréchet M., 1927, Ann. Soc. Polon. Math. 6, 93
• Gnedenko (1943) Gnedenko V. B., 1943, Annals of Mathematics 44, 423
• Gomez et al. (2012) Gomez P.L., et al., 2012, AJ, 144, 79
• Guillet & Teyssier (2011) Guillet, T., & Teyssier, R. 2011, Journal of Computational Physics, 230, 4756
• Gumbel (1935) Gumbel E. J., 1935, Ann. Inst. Henri Poincaré 5(2), 115
• Gumbel (1958) Gumbel E. J. 1958, Statistics of Extremes, Columbia Univ. Press
• Harrison & Coles (2011) Harrison I., Coles P., 2011, MNRAS, 418, L20
• Holz & Perlmutter (2012) Holz D. E., Perlmutter S., 2012, ApJL, 755, L36
• Hotchkiss (2011) Hotchkiss S., 2011, JCAP, 07, 004
• Hayashi & White (2006) Hayashi E., White S.D., 2006, MNRAS, 370, L38
• Hoyle, Jimenez & Verde (2011) Hoyle B., Jimenez R., Verde L., 2011, PRD, 83, 103502
• Jee et al. (2009) Jee G., et al., 2009, ApJ, 704, 672
• Jimenez & Verde (2009) Jimenez R., Verde L., 2009, PRD, 80, 127302
• Juszkiewicz, Springel & Durrer (1999) Juszkiewicz R., Springel V., Durrer R., 1999, ApJ, 518, L25
• Kowalski et al. (2008) Kowalski M. et al., 2008, ApJ, 686, 749
• Kravtsov & Borgani (2012) Kravtsov A.V., Borgani S., 2012, Annu. Rev. Astron. Astrophys., 50, 353
• Kunz et al. (2004) Kunz M, et al., 2004, PRD, 70, 041301
• Lacey & Cole (1994) Lacey C., Cole S., 1994, MNRAS, 271, 676
• Lage & Farrar (2013) Lage, C., & Farrar, G. 2013, arXiv:1312.0959
• Lee & Komatsu (2010) Lee J., Komatsu, E., 2010, ApJ, 718, 60
• Lima & Hu (2004) Lima, M. & Hu, W. 2004, PRD 70, 043504
• Machado & Lima Neto (2013) Machado, R. E. G., & Lima Neto, G. B. 2013, MNRAS, 430, 3249
• Majumdar & Mohr (2004) Majumdar, S. & Mohr, J. J. 2004, ApJ 613, 41
• Markevitch et al. (2002) Markevitch M., et al., 2002, ApJ, 567, L27
• Markevitch et al. (2006) Markevitch M., et al., 2006, astro-ph/0511345
• Mastropietro & Burkert (2008) Mastropietro C., Burkert A., MNRAS, 289, 967
• Menanteau et al. (2012) Menanteau F., et al., 2012, ApJ, 748, 7
• Menanteau et al. (2013) Menanteau F., et al., 2013, ApJ, 765, 67
• Milosavljevic et al. (2007) Milosavljevic M., et al., 2007, ApJ, 661, L131
• More et al. (2011) More S., Kravtsov A. V., Dalal N., & Gottlöber, S. 2011, ApJS, 195, 4
• Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
• Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
• Pacaud et al. (2006) Pacaud, F., et al. 2006, MNRAS, 372, 578
• Pierre et al. (2012) Pierre, M., et al. 2011, MNRAS 414, 1732
• Planck Collaboration (2011) Planck Collaboration 2011, A&A 536, 28
• Prunet et al. (2008) Prunet, S., Pichon, C., Aubert, D., Pogosyan, D., Teyssier, R., & Gottloeber, S., 2008, ApJs, 178, 179
• Rasera et al. (2014) Rasera, Y., Corasaniti, P.-S., Alimi, J.-M., et al. 2014, MNRAS, 440, 1420
• Ratra & Peebles (1988) Ratra B., Peebles P.J.E., 1988, PRD, 37, 3406
• Roy, Bouillot & Rasera (2014) Roy, F., Bouillot, V. R. & Rasera, Y. 2014, A&A, 564, A13
• Reichardt et al. (2012) Reichardt, C. L., et al. 2012, arXiv:1203.5775
• Rosati et al. (2009) Rosati P., et al., 2009, A&A, 508, 583
• Sachs & Wolfe (1967) Sachs R.K., Wolfe A.M., 1967, ApJ, 147, 73
• Scarrott & MacDonald (2012) Scarrott C., MacDonald A., 2012, Statistical Journal, 10, 33
• Sheth (1996) Sheth R.K., 1996, MNRAS, 279, 1310
• Smith et al. (2003) Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311
• Spergel et al. (2007) Spergel, D. N., et al., 2007, ApJs, 170, 377
• Springel et al. (2005) Springel V., et al., 2005, Nature, 435, 629
• Springel & Farrar (2007) Springel V. and Farrar G. R., 2007, MNRAS, 380, 911
• Stalder et al. (2013) Stalder, B., Ruel, J., Šuhada, R., et al. 2013, ApJ, 763, 93
• Summers, Davis & Evrard (1995) Summers F.J., Davis M., Evrad A.E., 1995, ApJ, 454, 1
• Thompson & Nagamine (2012) Thompson R., Nagamine K., 2012, MNRAS, 419, 3560
• Teyssier (2002) Teyssier R., 2002, Astron. Astrophys., 385, 337
• Waizmann, Ettori & Moscardini (2011) Waizmann J.-C., Ettori S., Moscardini L., 2011, MNRAS, 418, 456

## Appendix A Properties of Extremal Halo Pairs in DEUS-FUR simulations

In this Appendix, we present the characteristics of pairs of halos detected in the DEUS-FUR simulations at different redshifts for different cosmological models. For each pair we quote the relative velocity v, the mass of the main halo and the satellite , the mass ratio, the distance separation and the colliding angle. We only consider pairs with distance separation h Mpc.

## Appendix B Details of the multivariate distribution

In this Appendix, we present the maximum relative velocity or the maximum average mass of the pairs of halos detected in the DEUS-FUR simulations at different redshifts for different cosmological models having set a prior on the minimum mass (Table 8) or on the minimum relative velocity (Table 9). This table selects halos following the same procedure as in section 5 and aims at highlighting some Bullet Cluster candidates in the three DEUS-FUR cosmologies. For each sample we quote the maximal relative velocity v. We only consider pairs with distance separation h Mpc.

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters