Statistics of extreme objects in the Juropa Hubble Volume simulation††thanks: http://jubilee-project.org
We present the first results from the JUropa huBbLE volumE (Jubilee) project, based a large N-body, dark matter-only cosmological simulation with a volume of , containing 6000 particles, performed within the concordance CDM cosmological model. The simulation volume is sufficient to probe extremely large length scales in the universe, whilst at the same time the particle count is high enough so that dark matter haloes down to can be resolved. At we identify over 400 million haloes. The cluster mass function is derived using three different halofinders and compared to fitting functions in the literature. The distribution of clusters of maximal mass across redshifts agrees well with predicted masses of extreme objects, and we explicitly confirm that the Poisson distribution is very good at describing the distribution of rare clusters. The Poisson distribution also matches well the level to which cosmic variance can be expected to affect number counts of high mass clusters. We find that objects like the Bullet cluster exist in the far-tail of the distribution of mergers in terms of relative collisional speed. We also derive the number counts of voids in the simulation box for , and .
keywords:cosmology: large-scale structure of Universe—galaxies: haloes—galaxies: clusters—methods: numerical
Surveys mapping a substantial portion of the observable Universe (e.g. BOSS - Dawson et al. 2013, WiggleZ - Drinkwater et al. 2010, BigBoss - Schlegel et al. 2009, PanSTARRS - Magnier et al. 2013, DES - Mohr et al. 2012, PAU - Benítez et al. 2009, LSST - LSST Dark Energy Science Collaboration 2012, Euclid - Amiaux et al. 2012, etc.) aim to constrain the cosmological model to unprecedented accuracy. As they will be able to capture very faint objects they will have a shot-noise level low enough to be close to sampling variance-limited. This requires impressive handling of every step of the observational pipeline in order to limit the possibility of systematic errors that may degrade the information contained in them. Aside from these observational efforts, there will also be a similarly high demand placed on our ability to generate theoretical predictions that are equally accurate. This undoubtedly calls for numerical simulations of cosmic structure formation that resolve galactic scales in volumes comparable to the ones covered by these surveys. This is a non-trivial task. A simulation must cover a wide dynamic-range in order to accurately sample large-scale structure (LSS) in the universe. In particular, simulations need to resolve dark matter haloes, which are believed to host the observed galaxies, groups and clusters of galaxies, and accurately model the physics of galaxy formation and other non-linear physics, whilst adequately sampling large-scale matter fluctuations. Only recently have such simulations become feasible and nowadays full-box simulations of considerable fractions of the observable Universe are being conducted utilising close to a trillion particles (for a review of dark matter N-body simulations see Kuhlen et al. 2012).
Whilst a careful comparison of the statistical clustering properties of objects, in particular galaxies, will put tighter constraints on the parameters of any cosmological model, it is worth noting that the mere existence of individual outliers might pose challenges. Following observations of a series of apparently extreme objects such as high-mass clusters (for example XMMU J2235.3-2557 – a cluster with mass at redshift – Mullis et al. 2005; Rosati et al. 2009), or high-velocity collisions (for example the Bullet Cluster – see § 3.4) some authors have claimed that such objects are highly unlikely to exist in a concordance CDM cosmology and hence pose a challenge to its validity (Jimenez & Verde 2009; Lee & Komatsu 2010; Cayón et al. 2011; Enqvist et al. 2011; Hoyle et al. 2011; Holz & Perlmutter 2012). However, others have found no tension in their analyses (Harrison & Coles 2011; Mortonson et al. 2011; Waizmann et al. 2012b) and it has been recently pointed out by Hotchkiss (2011) that many studies into high-mass cluster observations have used biased statistical methods, a result that has been corroborated by other studies (Waizmann et al. 2012a; Hoyle et al. 2012; Stalder et al. 2013). The case of the Bullet cluster is less clear and we discuss it in § 3.4. Debate of this nature highlights the need to clarify our understanding of the statistics of such rare objects, for example by using large cosmological simulations. But this again requires simulations of large enough volumes and sufficient resolution to properly capture the likelihood of the formation of such rare clusters. Theoretical (as opposed to numerical) studies of such objects are challenging too due to their rarity and highly non-linear nature.
In this work we present one of the largest cosmological dark matter only simulations to date, the so-called Jubilee Universe, consisting of 6000 particles in a cubical volume of side-length 6. In this paper we focus on a presentation of the simulation itself and its general properties with respect to the cosmic web, clustering properties, and halo statistics. Subsequent papers will focus on topics including the generation of mock catalogues of LRGs and Quasars, a calculation of the Integrated Sachs-Wolfe effect signal and its cross-correlation to LSS, the weak lensing signal, and the SZ effect (see Watson et al. 2013b, for initial results in some of these areas). This paper is laid out as follows. In § 2 we outline our methodology for running the simulation and deriving from it results including halo and void catalogues. In § 3 we present our main results and in § 4 briefly discuss their potential implications.
The results presented in this work are based on two large-scale structure N-body simulations, whose parameters are listed in Table 1. Our main simulation has (216 billion) particles in a volume of . The particle mass is , yielding a minimum resolved halo mass (with 20 particles) of , corresponding to galaxies slightly more massive than the Milky Way. Luminous Red Galaxies (LRGs; ) are resolved with 100 particles, and galaxy clusters () are resolved with particles or more. This main simulation is accompanied by a second, smaller ‘control’ one with (29 billion) particles in a volume of and exactly the same minimum resolution. We used the CUBEPM N-body code, a PM (particle-particle-particle-mesh) code (Harnois-Deraps et al. 2012). It calculates the long-range gravity forces on a 2-level mesh and short-range forces exactly, by direct summation over local particles. The code is massively-parallel, using hybrid (combining MPI and OpenMP) parallelization and has been shown to scale well up to tens of thousands of computing cores (see Harnois-Deraps et al. 2012, for complete code description and tests). Both simulations and most analyses were performed on the Juropa supercomputer at the Jülich Supercomputing Centre in Germany (17,664 cores, 53 TB RAM, 207 TFlops peak performance) and required approximately 70,000 and 1.5 million core-hours for the and boxes respectively. The larger simulation was run on 8,000 computing cores (1,000 MPI processes, each with 8 OpenMP threads), and the smaller one on 2,048 cores.
We base our simulation on the 5-year WMAP results (Dunkley et al. 2009; Komatsu et al. 2009). The cosmology used was the ‘Union’ combination from Komatsu et al. (2009), based on results from WMAP, baryonic acoustic oscillations and high-redshift supernovae; i.e. , , , , , . These parameters are similar to the recent cosmology results of the Planck collaboration (Planck Collaboration 2013), where, considering a combination of data from Planck, WMAP, and LSS surveys (showing baryon acoustic oscillations) the parameters were calculated to be: , , , , and . The power spectrum and transfer function used for setting initial conditions was generated using CAMB (Lewis et al. 2000). The initial condition generator employed in the run uses first-order Lagrangian perturbation theory (1LPT), i.e. the Zel’dovich approximation (Zel’dovich 1970), to place particles in their starting positions. The initial redshift when this step takes place was . For a more detailed commentary on the choice of starting redshift for this simulation see Watson et al. (2013a).
We use two complementary definitions of haloes in this study. The first is the Spherical Overdensity (SO) definition of Lacey & Cole (1993). In this approach haloes are taken to be spheres that have overdensities that are above a chosen threshold, . The mass enclosed in these spheres is then given by
where is the radius of the halo and is the background matter density in the universe. We choose the overdensity threshold to be , i.e. an overdensity of 178 times the background matter density. This is a common choice motivated from the top-hat model of non-linear collapse in an Einstein de-Sitter (EdS) universe (Gunn & Gott 1972).
The second halo definition we adopt is that of the Friends-of-Friends (FOF) algorithm, first proposed by Davis et al. (1985). Haloes defined by this algorithm are identified within a simulation volume as agglomerations of particles that lie within a certain parameterised distance from one another. This distance is typically defined as the ‘linking-length’ the mean interparticle separation of particles in the simulation. Groups of particles within this distance of each other are identified as individual dark matter haloes. For our FOF haloes we follow various previous authors (Jenkins et al. 2001; Reed et al. 2003, 2007; Crocce et al. 2010; Courtin et al. 2011; Angulo et al. 2012) and use a linking length of 0.2. For further analysis on how the choice of halofinding parameters affect the mass function see Watson et al. (2013a) and references therein.
We employ three halofinding codes in our analysis: CUBEPM’s own on-the-fly SO halofinder (hereafter ‘CPMSO’) (Harnois-Deraps et al. 2012), the Amiga Halo Finder (hereafter ‘AHF’, Gill et al. 2004; Knollmann & Knebe 2009), and the FOF halofinder from the Gadget-3 N-Body cosmological code (Springel 2005, an update to the publicly available Gadget-2 code).
The CPMSO halofinder utilises a fine mesh from the CUBEPM code (with spacing twice as fine as the mean interparticle separation) and an interpolation scheme to identify local peaks in the density field. The code first builds the fine-mesh density using either Cloud-In-Cell (CIC) or Nearest-Grid-Point (NGP) interpolation. It then proceeds to search for and record all local density maxima above a certain threshold (typically set to 100 above the mean density) within the physical volume. It then uses quadratic interpolation on the density field to determine more precisely the location of the maximum within the densest cell. The halo centre determined this way agrees closely with the centre-of-mass of the halo particles. Each of the halo candidates is inspected independently, starting with the highest peak. The grid mass is accumulated in spherical shells of fine grid cells surrounding the maximum until the mean overdensity within the halo drops below . While the mass is accumulated it is removed from the mesh, so that no mass element is double-counted. For further details on the CPMSO method see Harnois-Deraps et al. (2012).
The halofinder ahf111ahf is freely available from http://www.popia.ft.uam.es/AMIGA (amiga Halo Finder) is a spherical overdensity finder that identifies (isolated and sub-)haloes as described in Gill, Knebe & Gibson (2004) and Knollmann & Knebe (2009). It employs a recursively refined grid to locate local overdensities in the density field. The identified density peaks are then treated as centres of prospective haloes. The resulting grid hierarchy is further utilised to generate a halo tree containing the information of which halo is a (prospective) host and subhalo, respectively. Halo properties are calculated based on the list of particles asserted to be gravitationally bound to the respective density peak. For a comparison of its performance to other finders in the field we refer the reader to Knebe et al. (2011); Onions et al. (2012); Knebe et al. (2013).
The specifics of the FOF halofinder packaged in with the Gadget-3 code currently have not been detailed in any publication but the algorithm itself is outlined in Davis et al. (1985). The main difference in the algorithm that exists in the Gadget-3 version is that the code is parallelised for distributed-memory machines. Specifically, haloes are found in local subvolumes of the simulation assigned to individual MPI tasks (created using the Gadget-3 domain decomposition which utilises a space-filling Peano-Hilbert curve – for details see the Gadget-2 paper, Springel 2005) and then haloes that extend spatially beyond the edges of the subvolumes are linked together in a final MPI communication step. We have altered the Gadget-3 code to read CUBEPM’s particle output format and significantly reduced its memory footprint by stripping away extraneous data structures.
Due to limitations in the scaling of the codes with processor numbers and the large memory footprint of the Jubilee simulation it was necessary to split the simulation time-slices into 27 subvolumes and run the halofinding algorithms on each subvolume independently. Each subvolume included a buffer zone which overlapped with the neighbouring ones, for correct handling of haloes straddling two or more sub-regions. We then stitched the subvolumes back together to create the final AHF and FOF halo catalogues, removing any duplicated structures in the overlapping buffers. This approach allows the handling of much larger amounts of data than otherwise possible and provides additional flexibility in terms of computational resources needed for post-processing.
2.3 Void finding
The formation of structure in the universe is a hierarchical process: small objects form, grow by accretion and merging and form more and more massive objects up to clusters of galaxies. Between the clusters large filaments can be seen both in observational data as well as in numerical simulations (Figure 1). These filaments surround large regions of low density which do not contain objects as massive as the ones found in the filaments or the knots at the end of filaments. These low density regions – voids – are the most extended objects in the universe. There are many different ways to define voids and correspondingly there are many different void-finding algorithms (for a review, see Colberg et al. 2008). In the following we are interested in the largest spherical regions of the universe which do not contain any object above a certain threshold in mass. In principle, one could extend this definition of voids also to non-spherical regions, however in this case one can get arbitrary volumes depending on the shape allowed. Since we are interested in the void function we restrict ourselves to spherical voids which are described only by one parameter, their radius. We identify voids in a sample of point-like objects distributed in space. Here these objects are our AHF haloes above a certain mass, but one could also use galaxies above a certain luminosity. Thus our voids are characterised by a threshold mass. If one decreases this mass threshold the number of objects increases and the size of the void decreases. In fact, a given void defined with objects at a higher mass becomes decomposed into many smaller voids defined in the distribution of lower mass objects (Gottlöber et al. 2003). This reflects the scale-free nature of structure formation. The algorithm searches first for the largest empty sphere then repeats taking into account the previously found voids so that no region is double-counted.
2.4 Online databases
It is our intention to make the data from the Jubilee simulation publicly available. This data will consist of three complementary halo catalogues of CPMSO, AHF and FOF haloes, in addition to LRG catalogues derived from the halo data and a catalogue of voids. The CPMSO will be available across a wide number of redshifts ( from , whereas the AHF and FOF data will be initially available only for . Further datasets will include smoothed density fields and maps of weak lensing and ISW signals. An SQL database has been set up so that the data can be queried to suit the requirements of individual users. Further information can be found at the Jubilee project website: http://jubilee-project.org.
3.1 Large-scale structure and the Cosmic Web
In Figure 1 we show a slice of the Cosmic Web at extracted from our 6 simulation. Perhaps most striking is the homogeneity of the matter distribution at large scales. This is expected from the cosmological principle which states that, on large enough scales, the universe is homogeneous and isotropic. On smaller scales significant features in the density field can be observed including voids, walls, filaments and clusters. For any simulated observations (for example of ISW and weak lensing signals, Watson et al. 2013b) we place virtual observers inside the simulation volume at a given location. As can be seen in Figure 1 the full sky as observed by an observer will show a highly homogeneous distribution of galaxies past a proper distance of a few hundred Mpc (i.e. a redshift of around ).
In Fig. 2 we show the evolution of the power spectra of the density field, , from redshift to . The particles were interpolated onto a regular grid of cells using the cloud-in-cell (CIC) interpolation scheme. From these data we then applied a correction for aliasing and the CIC window function and another for the effect of Poisson noise, all based on the prescription laid out in Jing (2005).
The baryonic acoustic oscillation (BAO) scale, , is well within the simulation box size and the BAO are visibile in the power spectra. At high redshift, , the power spectrum is largely linear, except at the smallest scales (), where the power grows faster than the linear growth factor predicts. As the hierarchical structure formation proceeds this non-linearity scale propagates to ever larger scales, reaching at redshift , and thereby affecting the BAO scale.
We calculate halo mass functions using the three halofinding algorithms outlined in § 2.2. Figure 3 shows the residuals between our haloes and two fits from the literature at . We compare our CPMSO and AHF haloes to the Tinker et al. (2008) mass function, noting that the Tinker et al. (2008) fit was calibrated to haloes with an overdensity criteria of versus the background matter density and our haloes were calculated using a value of . Despite this we see a good correspondence to within between the Tinker et al. (2008) fit and our AHF data for haloes with particle counts greater than . For the very largest haloes there is evidence that the Tinker et al. (2008) fit may be overpredicting the mass function, although this is where shot noise begins to severely affect number counts of objects. The CPMSO data follows a similar trend to the AHF data. We overlay on these plots two of the fits from Watson et al. (2013a) – mass function results calibrated to data that included the Jubilee haloes presented here. The fit used for the left-hand panel is the redshift-dependent fit based on the CPMSO halofinder, the fit for the central panel a fit based on AHF results for (see Watson et al. 2013a, for further details).
We compare the FOF results to those of the Millennium, Millennium-2 and Millennium-XXL simulations (Angulo et al. 2012), the latter containing particles in a box with length 3Gpc. The FOF halo data shows agreement to within with the Angulo et al. (2012) fit for haloes with 300 or more particles. The FOF haloes are being compared with a linking length of 0.2, which makes the similarity between the mass functions a good test of the validity of the Jubilee halo distribution as this was the same choice made in Angulo et al. (2012). For a more detailed study of the mass function across a broad redshift range, including results from the 6Gpc simulation, see Watson et al. (2013a).
3.2 Cosmic variance
Due to the large size of our simulated volume we are able to quantify cosmic variance on scales smaller than our box size in terms of the number counts of objects one expects to find in a given volume. To that end we have compared halo counts in different mass bins in different sized subvolumes. We chose the subvolumes such that they filled the entire full-box with no overlap. The results are shown in Figure 4. We show the 1 standard-deviation error in the number counts of haloes by mass bin relative to our entire volume for sub-box lengths of 3, 2, 1 and 0.5 Gpc. These choices directly compare to the box lengths of some contemporary simulations (Millennium-XXL - Angulo et al. 2012, Horizon - Teyssier et al. 2009, MultiDark - Prada et al. 2012, and Millennium - Springel et al. 2005 respectively). We also show a prediction for this error calculated by assuming that the halo number counts follow the CPMSO redshift parameterised mass function from Watson et al. (2013a), and by assuming that the observed error in number counts follows a Poisson distribution. This theoretical prediction matches the high-mass data very well. For lower masses the error becomes dominated by sample variance, as discussed in Smith & Marian (2011), and the Poisson prediction presented here begins to break down. We discuss how well the Poisson distribution matches the counts of high mass clusters in § 3.3 below.
As expected, the error is minimal for lower mass haloes and increases for rarer objects. At the 0.5 Gpc box has an error of under 10% up until haloes of mass around , while for box lengths of 1, 2 and 3 Gpc, a 10% error in number counts per mass interval is realised at around and respectively. The errors are exacerbated at higher redshifts, due to the haloes of a fixed size growing more rare at earlier times.
One subtlety should be mentioned: because the sub-volumes considered in this analysis were derived from a larger simulation volume they include the effect of matter fluctuations that exist on scales larger than their box lengths. We stress here that this is not the case for simulations with equivalent volumes to these sub-volumes, as modes of power in the density field that are larger than the box length of a simulation are typically set to zero. This implies that the variation in number counts presented here is slightly different from one that occurs due to a lack of appropriate large-scale power in a simulation volume. This mis-representation of reality (by all simulations, including the Jubilee despite its large volume) leads to an additional set of errors but is, fortuitously, only an issue for very small volume simulations with box lengths of the order of up to a few tens of Mpc (for example see Yoshida et al. 2003; Barkana & Loeb 2004; Sirko 2005; Power & Knebe 2006; Bagla & Prasad 2006; Lukic et al. 2007). Observational volumes, sampling the universe, do not suffer from this effect and the results presented here can be expected to translate reasonably well into counts of high-mass objects in LSS surveys.
3.3 Statistics of rare objects
A current topic in cosmology that relates to the number counts of very high mass objects is that of whether large observed clusters are in conflict with the standard CDM model. In Figure 5 we show a theoretical prediction for the expected distribution of maximal mass clusters. This prediction was created using the Extreme Value Statistics (EVS) prescription of Harrison & Coles (2011) and the redshift parameterised redshift-dependent mass function, based on CPMSO haloes, from Watson et al. (2013a). It is comparable to the plot in Figure 1 of Harrison & Coles (2012) which was created using the mass function from Tinker et al. (2008), except for the fact that the mass of the haloes in this version is taken to be set by the overdensity criterion, rather than the criterion used in Harrison & Coles (2012). The black data points shown on Figure 5 correspond to the largest clusters observed in the Jubilee simulation by a central observer, also based on a , as per the configuration of the CPMSO halofinder. The redshift shells for both the EVS contours and the maximal mass clusters are identical. The data all lies within the 3 range showing the expected result that the there is no tension between objects observed in a CDM cosmological simulation and the theoretical expectation from Extreme Value Statistics. Of interest is whether there are observed clusters in the Universe that have masses that are in tension with the CDM model. To date observations have shown this to not be the case, as shown in a systematic review by Harrison & Hotchkiss (2012).
The question of how well the Poisson distribution fits our rare cluster number counts is addressed in Figure 6. The simulation volume at was split up into 5438 independent subvolumes. For each subvolume we calculated the number of objects above a given threshold mass (, and for the panels in Figure 6, left-to-right respectively) found in each subvolume. The mass thresholds were chosen so that only a very small number (around 0–2) of objects were found in each subvolume, which represents the regime where we expect Poisson statistics to be dominant. We then compared the histogram of the measured distribution of the objects in the simulation to that predicted by a Poisson distribution with a mean set by the average across all the subvolumes. The correspondence between the two is very close. This is an interesting result as it validates the common choice of Poisson statistics for describing the expected distribution of these objects, and is the first time it has been validated using a simulation of this scale (for a detailed investigation of the applicability of the Poisson distribution in cluster counts across different masses see Smith & Marian 2011, who used simulations of box length for their study).
3.4 High Mergers and the Bullet Cluster
There has been recent debate regarding whether the Bullet Cluster (1E0657-56, which resides at a redshift of ) poses a challenge to the CDM model. 1E0657-56 consists of a large cluster of mass and a sub-cluster – the ‘bullet’ – of mass that has traversed through the larger cluster, creating a substantial bow shock along the way (Markevitch et al. 2002; Barrena et al. 2002; Clowe et al. 2004, 2006; Bradač et al. 2006). Tension with CDM arises from the calculated value for the speed of the shock of (Markevitch et al. 2002, 2004; Springel & Farrar 2007), which was originally calculated to be too high for a CDM universe (Farrar & Rosen 2007) – whereas it might be better accommodated in alternative cosmologies (e.g. Llinares et al. 2009). Other studies have concluded that the velocity is not in tension with CDM (Hayashi & White 2006). An important clarification of this issue was presented by groups working on simulations of Bullet-like systems (Takizawa 2005, 2006; Milosavljević et al. 2007; Springel & Farrar 2007; Mastropietro & Burkert 2008) where, in general, it was found that the shock speed was substantially higher than the speed of the mass centroid of the infalling subcluster. For example Springel & Farrar (2007) found that a Bullet-like system in their simulations had a shock speed of whereas the sub-cluster had a speed of only . Milosavljević et al. (2007) found that in an illustrative simulation the sub-cluster CDM halo had a speed that was 16% lower than that of the shock.
Even given this moderation of the extreme sub-cluster speed in 1E0657-56 there have still been claims in the literature that the CDM model may be incapable of creating such a system (Lee & Komatsu 2010; Thompson & Nagamine 2012). This is not wholly unexpected as a) Mastropietro & Burkert (2008) have shown that the properties of the bow shock are not well described by simulations and b) even with a moderation in sub-cluster speed along the lines of Springel & Farrar (2007) or Milosavljević et al. (2007), the speed may still be too high for the CDM model to accommodate. These studies have relied on numerical simulations to observe the distribution of relative velocities in colliding clusters. From these distributions 1E0657-56 can be assessed and deemed to be either rare for a CDM universe or so rare that it puts the whole model in doubt.
Alternative approaches have also been taken in addressing this question. Forero-Romero et al. (2010) looked in 2-D-projected position-space for Bullet-like systems in the MareNostrum Universe, a large hydrodynamical cosmological simulation. The characteristic distribution of gas and dark matter in 1E0657-56, as projected on the sky – with a large displacement between the cluster’s gas and dark matter – was found to be expected in 1% - 2% of clusters with masses larger than . Nusser (2008) performed a ‘back in time’ analysis to place bounds on the relative overdensity the system resides in the universe, concluding that for a relative speed of the system would need to have a mass of and exist in a local overdensity of 10 times the background density of the universe.
Here we use the huge number counts of clusters in the Jubilee simulation to add to the debate. We consider AHF (sub-)haloes with mass greater than that are colliding with (host) haloes of mass greater or equal to at . Our results are shown in Figure 7, along with the original Bullet speed presented in Markevitch et al. (2004), and the moderated result from Springel & Farrar (2007), which represents the lowest value from the literature to-date. We show in blue haloes that are colliding pairs and in red haloes that are a colliding subhalo and halo pair. We have added a normally distributed random scatter to our velocities with a width given by the error in the observed value for . This is to mimic the effect of Eddington bias in our simulated data. This observational bias arises when observing extreme measurements in a distribution of measurements all with some associated scatter. As there are many more data points that exist with less extreme velocities than the one in question it is likely that an extreme data point is an upscattered less extreme one. As we know very precisely the pairwise velocities of halo pairs in our simulation adding random scatter to this distribution serves to create the effect in our measurement.
As can be seen from the distribution the Bullet cluster is an extreme object, but only when the radial separation of the halo pairs is considered. We find many candidate mergers in our volume with a collision speed that equals or exceeds the more conservative speed estimate for the cluster, and a few objects that are not far from the higher velocity estimate of Markevitch et al. (2004). However, we find no objects that, at a closer separation, give rise to a large enough merging velocity. This is likely to be due to the effect of only considering one simulation output in our analysis. At any given output redshift only a handful of haloes will be undergoing a major merger event of the sort we are interested in and this is reflected by the paucity of data points that lie at a separation of less than 0.6 Mpc. It is likely, therefore, that over the course of the Jubilee simulation run, high velocity mergers of the type observed in the Bullet Cluster do occur.
This result is in line with previous attempts to use large cosmological simulations to address this issue where the bullet was not found to be extreme (Hayashi & White 2006; Thompson & Nagamine 2012). Interestingly, Thompson & Nagamine (2012) extrapolated their results from smaller simulation volumes and concluded that a volume of would be required in order to observe a Bullet-like cluster.
The conclusion that we put forward based on this result is that there is at present no tension between our data and the standard cosmological model. This conclusion is tentative, however, and there would appear to be a need for careful further research into this question based on a number of points. Firstly, the results are very sensitive to the mass cuts imposed on the candidate search. It would be very difficult to find a precise analogue to the Bullet Cluster in terms of the masses, velocities, and spatial separation of the haloes. Here we have taken a cut-off in mass that allows us to search for Bullet-like systems rather than a precise Bullet Cluster. Secondly, we have placed no restrictions on the directions of the relative velocities of the halo pairs. The bow shock observed in the Bullet Cluster has arisen from the Bullet subhalo having passed through the parent halo (it is this occurrence, that fortuitously lies almost in the plane of the sky as we observe it, that has allowed us to identify the relative pairwise velocities of the halo and subhalo in the system). In our analysis we plot all the pairwise velocities of the haloes, making no distinction between haloes that are infalling and haloes that have already undergone a collision and not considering how the orientations of the collisions might appear to a specific observer. This is a fair way to assess the data as the actual collision in a Bullet Cluster-like system is expected to take a few hundred Megayears so is a relatively short event. Canvasing all our haloes in this manner assesses whether there is likely to be or whether there has been a Bullet Cluster-like collision in the simulation around . Further studies should include how random observers would observe these events. Lastly, halofinding algorithms are notoriously sub-optimal when trying to find and separate haloes that are merging (this is discussed in detail in Knebe et al. 2011, we draw the reader’s attention to Figure 10 from that paper in particular). However, as the separation in question of the two haloes is relatively large, this is unlikely to be affecting our results.
3.5 The Jubilee void function
The distribution of voids for a given threshold is characterised by the void function, the number of spheres with radii larger than per volume. We have studied the void distribution at redshifts , 0.5, and 1. At we have identified the voids in the distribution of haloes more massive than , , , and . At redshift we identified 244,989, 1,753,982, 5,596,627, 91,615,821 haloes more massive than these thresholds, respectively. Thus the mean distance between them (i.e. the box length divided by the cube root of the number counts) is about 96Mpc, 50Mpc, 34Mpc, and 13Mpc. Nevertheless, we found huge volumes which do not contain any of these objects.
In Figure 8, top panel, we show the void functions at for four different threshold masses. For the largest threshold we find a few very large spheres with radii of 150Mpc which do not contain any cluster more massive than . For smaller thresholds the void function is very steep, i.e. there are a number of voids with a volume almost as large as the volume of the largest voids defined by the threshold. This means that the voids are almost uniformly distributed, as there are so many of a similar size. At higher redshifts (middle and bottom panel of Figure 8) we observe similar behaviour but, due to the evolution of the mass function, only with lower threshold masses. Note, that at the lowest threshold () the maximum void radius is almost redshift independent between and occurs at a void radius of about 40Mpc. This may seem in contradiction to the fact low density regions expand slightly faster than the mean expansion rate of the universe. However, since the tracers of the voids are also evolving the number of objects above the threshold evolves. For mass haloes, the number counts rise from 38,994,056 at to 91,615,821 at . Therefore, the mean distance shrinks from 18Mpc to 13Mpc and using this threshold mass we see the interesting result that the maximum void radius remains almost constant in time.
4 Summary and Discussion
In this paper we have presented a broad range of results from a CDM-based simulation. The results have focussed on predictions on very large scales, such as extremely massive clusters and large void regions. The simulation itself represents one of the largest undertaken to-date, with a volume of , and haloes resolved down to , a resolution that allows the creation of mock LRG and cluster catalogues.
The distribution of dark matter haloes in the Jubilee was found to be well described by fitting functions from the literature, and the dark matter haloes from the Jubilee have been used in a separate paper (Watson et al. 2013a) to construct mass function fits across a broad range of redshifts and volumes. For the rare tail of the mass function we have confirmed that the Poisson distribution describes well the number counts of objects. The masses of clusters with extremal masses in the Jubilee simulation were investigated across a range of redshifts and were found to agree well with both observation and theory, in particular the expected masses of the very largest objects found when using Extreme Value Statistics.
4.1 Implications for precision cosmology
An important prediction from this simulation is the expected effect of cosmic variance on the counts of massive clusters. This result can be used to gauge number-count errors in survey and simulation data. Understanding this is a vital component of the drive towards high precision cosmology. We showed in Figure 4 how the expected number of clusters in given volumes are likely to vary. In general quantifying the effect of cosmic variance in simulations is notoriously tricky due to the requirement for either multiple repeats of a simulation or for a large simulated volume (or preferably both of these). Due to the large scale of the Jubilee volume we are able to use the latter and do so in a manner that includes the long-wavelength modes of the matter distribution. The variation in cluster counts for smaller boxes or surveys is highly significant if one is investigating the distribution of high-mass objects such as galaxy clusters, which form an important cosmological probe.
4.2 The largest voids
Our largest void, defined using a threshold mass of , is across. To put this void in context it is around one fifth of the volume of the Millennium simulation and it contains no clusters with mass greater than . The probability, based on volume-occupation alone, of finding yourself within this void in the universe represented by the Jubilee simulation is . There have been investigations as to whether our occupying a local underdensity might explain the apparent existence of an accelerated expansion in the late-time universe (Ellis 1979; Mustapha et al. 1997; Zehavi et al. 1998; Tomita 2001; Iguchi et al. 2002; Barausse et al. 2005; Wiltshire 2005; Moffat 2005; Alexander et al. 2009; February et al. 2010; Marra & Pääkkönen 2010; Nadathur & Sarkar 2011). The void in question would need to have very specific characteristics that include its radius, sphericity, density and density profile. Predictions for these void parameters vary but have typically required the void to be of at least a few hundred Mpc in radius and, importantly, close to spherically symmetric, with us as observers very near its centre. This latter requirement is due to the type Ia supernovae data implying that dark energy is close to isotropic across the sky. We see from our void functions in Figure 8 that there are a few hundred voids in the Jubilee volume with radii , for the mass threshold. We estimate the proportion of the entire simulation volume taken up by voids with a radius of to be 0.04%. Adding an additional requirement that an observer occupy the central of the void volumes in question we arrive at the total spatial volume in the Jubilee box that would contain observers in the centre of voids of radii greater than to be . This is a rough statistical estimate and ignores the fact that observers might be better considered to only exist at the locations of galaxies in the simulation. In addition the simulation contains a dark energy component so has already modelled the effect of late-time accelerated expansion on structure formation. This latter point does not alter the order of magnitude of the result as void sizes in universes without dark energy are comparable to void sizes in CDM (Müller et al. 2000). We intend to look more closely into putting a probability on this figure for voids in the Jubilee simulation in a follow-up paper.
4.3 The CDM model versus observations
The distribution of most-massive clusters in the Jubilee was found to be in line with current theoretical predictions based on Extreme Value Statistics. The nature of Extreme Value Statistics is that it lacks predictive power in terms of constraining models, but it is a powerful method for ruling out models based on only a handful of extreme data points (various authors have previously implemented studies in cosmology based on it, for example Coles 1988; Antal et al. 2009; Colombi et al. 2011; Davis et al. 2011; Harrison & Coles 2011). Had the masses of observed clusters in Harrison & Hotchkiss (2012) lain significantly away from the expected EVS prediction then the CDM model would be immediately placed in doubt. One result in this paper, that of the extreme nature of the Bullet cluster, is suggestive of a possible tension with CDM. An Extreme Value Statistics approach is likely to cast this result in a more comprehensive context, but it is beyond the scope of this paper to attempt analysis along these lines.
In follow-up work we intend to investigate in more detail the existing discrepancy between observations of the ISW signal and the expected CDM signal. For example Hunt & Sarkar (2010) claim that the observed voids from SDSS data are too large for a CDM universe. This result was based largely on analysis of the ISW signal in Granett et al. (2008) and Granett et al. (2009). This highlights the intimate link between the void distribution and the ISW signal – underdense regions imprint themselves on the CMB via the ISW effect – and represents a current challenge to the CDM model.
The simulation was performed on the Juropa supercomputer of the Jülich Supercomputing Centre (JSC). We would like to give special thanks to Alexander Schnurpfeil as well the whole the Jülich support team for their help during the phase of code implementation as well as during the production runs. We thank Peter Coles, Aurel Schneider and Leonidas Christodoulou for their helpful comments and feedback on this work. We also thank Steffen Knollman for his help in applying the AHF halofinder to our data. We thank Volker Springel for allowing us to use the Gadget-3 code with its FOF halofinder. WW thanks The Southeast Physics Network (SEPNet) for providing funding for his research. ITI was supported by The SEPNet and the Science and Technology Facilities Council grants ST/F002858/1 and ST/I000976/1. AK is supported by the Spanish Ministerio de Ciencia e Innovación (MICINN) in Spain through the Ramón y Cajal programme as well as the grants AYA 2009-13875-C03-02, AYA2009-12792-C03-03, CSD2009-00064, CAM S2009/ESP-1496 and the Ministerio de Economa y Competitividad (MINECO) through grant AYA2012-31101. GY acknowledges support from MINECO (Spain) under research grants AYA2009-13875-C03-02, FPA2009-08958, AYA2012-31101 and Consolider Ingenio SyeC CSD2007-0050. He also thanks Comunidad de Madrid for support under PRICIT project ASTROMADRID (S2009/ESP-146). JMD and EMG acknowledge support of the consolider project CSD2010-00064 funded by the Ministerio de Economia y Competitividad and AYA2012-39475-C02-01.
- Alexander et al. (2009) Alexander S., Biswas T., Notari A., Vaid D., 2009, JCAP, 9, 025
- Amiaux et al. (2012) Amiaux J. et al., 2012. Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8442
- Angulo et al. (2012) Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046
- Antal et al. (2009) Antal T., Sylos Labini F., Vasilyev N. L., Baryshev Y. V., 2009, EPL (Europhysics Letters), 88, 59001
- Bagla & Prasad (2006) Bagla J. S., Prasad J., 2006, Mon. Not. R. Astron. Soc. , 370, 993
- Barausse et al. (2005) Barausse E., Matarrese S., Riotto A., 2005, Phys Rev D, 71, 063537
- Barkana & Loeb (2004) Barkana R., Loeb A., 2004, ApJ, 609, 474
- Barrena et al. (2002) Barrena R., Biviano A., Ramella M., Falco E. E., Seitz S., 2002, A.&A, 386, 816
- Benítez et al. (2009) Benítez N. et al., 2009, ApJ, 691, 241
- Bradač et al. (2006) Bradač M. et al., 2006, ApJ, 652, 937
- Cayón et al. (2011) Cayón L., Gordon C., Silk J., 2011, MNRAS, 415, 849
- Clowe et al. (2004) Clowe D., Gonzalez A., Markevitch M., 2004, ApJ, 604, 596
- Clowe et al. (2006) Clowe D., Bradač M., Gonzalez A. H., Markevitch M., Randall S. W., Jones C., Zaritsky D., 2006, ApJL, 648, L109
- Colberg et al. (2008) Colberg J. M. et al., 2008, MNRAS, 387, 933
- Coles (1988) Coles P., 1988, MNRAS, 231, 125
- Colombi et al. (2011) Colombi S., Davis O., Devriendt J., Prunet S., Silk J., 2011, MNRAS, 414, 2436
- Courtin et al. (2011) Courtin J., Rasera Y., Alimi J. M., Corasaniti P. S., Boucher V. et al., 2011, Mon. Not. R. Astron. Soc. , 410, 1911
- Crocce et al. (2010) Crocce M., Fosalba P., Castander F. J., Gaztañaga E., 2010, Mon. Not. R. Astron. Soc. , 403, 1353
- Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
- Davis et al. (2011) Davis O., Devriendt J., Colombi S., Silk J., Pichon C., 2011, MNRAS, 413, 2087
- Dawson et al. (2013) Dawson K. S. et al., 2013, AJ, 145, 10
- Drinkwater et al. (2010) Drinkwater M. J. et al., 2010, MNRAS, 401, 1429
- Dunkley et al. (2009) Dunkley J. et al., 2009, ApJS, 180, 306
- Ellis (1979) Ellis G. F. R., 1979, General Relativity and Gravitation, 11, 281
- Enqvist et al. (2011) Enqvist K., Hotchkiss S., Taanila O., 2011, JCAP, 4, 017
- Farrar & Rosen (2007) Farrar G. R., Rosen R. A., 2007, Physical Review Letters, 98, 171302
- February et al. (2010) February S., Larena J., Smith M., Clarkson C., 2010, MNRAS, 405, 2231
- Forero-Romero et al. (2010) Forero-Romero J. E., Gottlöber S., Yepes G., 2010, ApJ, 725, 598
- Gill et al. (2004) Gill S. P., Knebe A., Gibson B. K., 2004, Mon. Not. R. Astron. Soc. , 351, 399
- Gottlöber et al. (2003) Gottlöber S., Łokas E. L., Klypin A., Hoffman Y., 2003, MNRAS, 344, 715
- Granett et al. (2008) Granett B. R., Neyrinck M. C., Szapudi I., 2008, ApJL, 683, L99
- Granett et al. (2009) Granett B. R., Neyrinck M. C., Szapudi I., 2009, ApJ, 701, 414
- Gunn & Gott (1972) Gunn J. E., Gott III J. R., 1972, ApJ, 176, 1
- Harnois-Deraps et al. (2012) Harnois-Deraps J., Pen U. L., Iliev I. T., Merz H., Emberson J. et al., 2012, arXiv:1208.5098
- Harrison & Coles (2011) Harrison I., Coles P., 2011, MNRAS, 418, L20
- Harrison & Coles (2012) Harrison I., Coles P., 2012, MNRAS, 421, L19
- Harrison & Hotchkiss (2012) Harrison I., Hotchkiss S., 2012, arXiv:1210.4369
- Hayashi & White (2006) Hayashi E., White S. D. M., 2006, MNRAS, 370, L38
- Holz & Perlmutter (2012) Holz D. E., Perlmutter S., 2012, ApJL, 755, L36
- Hotchkiss (2011) Hotchkiss S., 2011, JCAP, 7, 004
- Hoyle et al. (2011) Hoyle B., Jimenez R., Verde L., 2011, Phys Rev D, 83, 103502
- Hoyle et al. (2012) Hoyle B., Jimenez R., Verde L., Hotchkiss S., 2012, JCAP, 2, 009
- Hunt & Sarkar (2010) Hunt P., Sarkar S., 2010, MNRAS, 401, 547
- Iguchi et al. (2002) Iguchi H., Nakamura T., Nakao K., 2002, Progress of Theoretical Physics, 108, 809
- Jenkins et al. (2001) Jenkins A., Frenk C. S., White S. D. M., Colberg J. M., Cole S., Evrard A. E., Couchman H. M. P., Yoshida N., 2001, Mon. Not. R. Astron. Soc. , 321, 372
- Jimenez & Verde (2009) Jimenez R., Verde L., 2009, Phys Rev D, 80, 127302
- Jing (2005) Jing Y. P., 2005, ApJ, 620, 559
- Knebe et al. (2011) Knebe A., Knollmann S. R., Muldrew S. I., Pearce F. R., Aragon-Calvo M. A. et al., 2011, arXiv:1104.0949
- Knebe et al. (2013) Knebe A. et al., 2013, MNRAS, 428, 2039
- Knollmann & Knebe (2009) Knollmann S. R., Knebe A., 2009, ApJS, 182, 608
- Komatsu et al. (2009) Komatsu E. et al., 2009, ApJS, 180, 330
- Kuhlen et al. (2012) Kuhlen M., Vogelsberger M., Angulo R., 2012, Physics of the Dark Universe, 1, 50
- Lacey & Cole (1993) Lacey C., Cole S., 1993, Mon. Not. R. Astron. Soc. , 262, 627
- Lee & Komatsu (2010) Lee J., Komatsu E., 2010, ApJ, 718, 60
- Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, Astrophys. J., 538, 473
- Llinares et al. (2009) Llinares C., Zhao H. S., Knebe A., 2009, ApJL, 695, L145
- LSST Dark Energy Science Collaboration (2012) LSST Dark Energy Science Collaboration, 2012, arXiv:1211.0310
- Lukic et al. (2007) Lukic Z., Heitmann K., Habib S., Bashinsky S., Ricker P. M., 2007, ApJ, 671, 1160
- Magnier et al. (2013) Magnier E. A. et al., 2013, ApJS, 205, 20
- Markevitch et al. (2002) Markevitch M., Gonzalez A. H., David L., Vikhlinin A., Murray S., Forman W., Jones C., Tucker W., 2002, ApJL, 567, L27
- Markevitch et al. (2004) Markevitch M., Gonzalez A. H., Clowe D., Vikhlinin A., Forman W., Jones C., Murray S., Tucker W., 2004, ApJ, 606, 819
- Marra & Pääkkönen (2010) Marra V., Pääkkönen M., 2010, JCAP, 12, 021
- Mastropietro & Burkert (2008) Mastropietro C., Burkert A., 2008, MNRAS, 389, 967
- Milosavljević et al. (2007) Milosavljević M., Koda J., Nagai D., Nakar E., Shapiro P. R., 2007, ApJL, 661, L131
- Moffat (2005) Moffat J. W., 2005, JCAP, 10, 012
- Mohr et al. (2012) Mohr J. J. et al., 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series. Vol. 8451
- Mortonson et al. (2011) Mortonson M. J., Hu W., Huterer D., 2011, Phys Rev D, 83, 023015
- Müller et al. (2000) Müller V., Arbabi-Bidgoli S., Einasto J., Tucker D., 2000, MNRAS, 318, 280
- Mullis et al. (2005) Mullis C. R., Rosati P., Lamer G., Böhringer H., Schwope A., Schuecker P., Fassbender R., 2005, ApJL, 623, L85
- Mustapha et al. (1997) Mustapha N., Hellaby C., Ellis G. F. R., 1997, MNRAS, 292, 817
- Nadathur & Sarkar (2011) Nadathur S., Sarkar S., 2011, Phys Rev D, 83, 063506
- Nusser (2008) Nusser A., 2008, MNRAS, 384, 343
- Onions et al. (2012) Onions J. et al., 2012, MNRAS, 423, 1200
- Planck Collaboration (2013) Planck Collaboration, 2013, arXiv:1303.5076
- Power & Knebe (2006) Power C., Knebe A., 2006, Mon. Not. R. Astron. Soc. , 370, 691
- Prada et al. (2012) Prada F., Klypin A. A., Cuesta A. J., Betancort-Rijo J. E., Primack J., 2012, MNRAS, 423, 3018
- Reed et al. (2003) Reed D., Gardner J., Quinn T. R., Stadel J., Fardal M. et al., 2003, Mon. Not. R. Astron. Soc. , 346, 565
- Reed et al. (2007) Reed D., Bower R., Frenk C., Jenkins A., Theuns T., 2007, Mon. Not. R. Astron. Soc. , 374, 2
- Rosati et al. (2009) Rosati P. et al., 2009, A.&A, 508, 583
- Schlegel et al. (2009) Schlegel D. J. et al., 2009, arXiv:0904.0468
- Sirko (2005) Sirko E., 2005, ApJ, 634, 728
- Smith & Marian (2011) Smith R. E., Marian L., 2011, arXiv:1106.1665
- Springel (2005) Springel V., 2005, Mon. Not. R. Astron. Soc. , 364, 1105
- Springel & Farrar (2007) Springel V., Farrar G. R., 2007, MNRAS, 380, 911
- Springel et al. (2005) Springel V. et al., 2005, Nature, 435, 629
- Stalder et al. (2013) Stalder B. et al., 2013, ApJ, 763, 93
- Takizawa (2005) Takizawa M., 2005, ApJ, 629, 791
- Takizawa (2006) Takizawa M., 2006, PASJ, 58, 925
- Teyssier et al. (2009) Teyssier R. et al., 2009, A.&A, 497, 335
- Thompson & Nagamine (2012) Thompson R., Nagamine K., 2012, MNRAS, 419, 3560
- Tinker et al. (2008) Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709
- Tomita (2001) Tomita K., 2001, MNRAS, 326, 287
- Waizmann et al. (2012a) Waizmann J. C., Ettori S., Moscardini L., 2012a, MNRAS, 420, 1754
- Waizmann et al. (2012b) Waizmann J. C., Redlich M., Bartelmann M., 2012b, A.&A, 547, A67
- Watson et al. (2013a) Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R., Yepes G., 2013a, MNRAS, 433, 1230
- Watson et al. (2013b) Watson W. A. et al., 2013b, arXiv:1307.1712
- Wiltshire (2005) Wiltshire D. L., 2005, arXiv:gr-qc/0503099
- Yoshida et al. (2003) Yoshida N., Sokasian A., Hernquist L., Springel V., 2003, ApJ, 598, 73
- Zehavi et al. (1998) Zehavi I., Riess A. G., Kirshner R. P., Dekel A., 1998, ApJ, 503, 483
- Zel’dovich (1970) Zel’dovich Y. B., 1970, A.&A, 5, 84