The journey of QSO haloes from $z\sim 6$ to the present

Our simulation, named the Millennium-XXL or MXXL, represents the matter content of a $\Lambda$CDM universe using more than $303$ billion particles ($6720^{3}$) within a comoving cube of side 4.1 Gpc. Our choice of cosmological parameters is identical to that used in the other Millennium Simulations (see Table 1, Springel et al., 2005; Boylan-Kolchin et al., 2009), implying a particle mass of $8.45\times 10^{9}M_{\odot}$. This set of parameters is inconsistent with the most recent constraints from observations of the Cosmic Microwave Background and low-redshift large-scale structure (Komatsu et al., 2011), but, as discussed in Section 2.2, our results can be scaled accurately to any nearby cosmology, including those that are now more favoured. The comoving softening length of the simulation was $\epsilon=13.7$ kpc, implying $\sim 10^{15}$ effective resolution elements in the full simulated volume (the forces are exactly Newtoninan beyond $2.8\times\epsilon$). The enormous statistical power and dynamical range that this implies is illustrated in Fig. 1, which shows the $z=0$ density field on four different scales, starting from the whole box in the background, then zooming progressively onto the most massive halo in the insets.

The initial phase-space distribution of the particles was set up at $z=63$ by perturbing a glass-like distribution (White, 1996; Baugh et al., 1995) using second order Lagrangian perturbation theory (2LPT Scoccimarro, 1998). The use of 2LPT has several advantages over the more common Zeldovich approximation (Zel’Dovich, 1970), including small and rapidly-decaying transients in the matter clustering and a better representation of the perturbations that seed extremely large dark matter haloes (Crocce et al., 2006; Knebe et al., 2009; Jenkins, 2010). The latter is particularly relevant for this paper. The amplitude of individual Fourier modes was set hierarchically as described in Jenkins (2012, in prep). This allows an efficient and consistent generation of initial conditions for the resimulation of any selected subregion at, in principle, arbitrarily high mass resolution.

Note that the MXXL is the first simulation with a large enough volume and a small enough particle mass to contain a representative sample of well-resolved ($\geq 1000$ particles, Trenti et al., 2010) haloes which could represent the hosts of the $z\sim 6$ SDSS quasars. Performing a simulation with these characteristics posed severe computational challenges with respect to raw execution time, scalability of the simulation algorithms, memory consumption and I/O performance. Thanks to a highly optimised simulation code and one of the largest supercomputers in Europe, these challenges were successfully met. Specifically, the MXXL was carried out with a special version of the GADGET-3 code (Springel, 2005), which aggressively reduced peak memory consumption at runtime, incorporated a number of analysis tools on the fly, and implemented a special compression of the output data. As a result the MXXL was completed in late summer 2010 in less than 3 million CPU hours (including postprocessing calculations), using 30 Tb of RAM, and $12,228$ cores of the Juropa cluster at the Jülich Supercomputer Center in Germany. We refer to Angulo et al. (2012) for more details about the simulation.

The MXXL was an extremely expensive numerical simulation and so could only be carried out once for a specific set of cosmological parameters. However, the rescaling method of Angulo & White (2010b) allows us to use the simulation to analyse any neighbouring cosmology with Gaussian initial fluctuations, and in this paper we will show results not only for the original MS cosmology but also for two other versions of the standard $\Lambda$CDM cosmology. The accuracy of the rescaling scheme is remarkably high if it is applied carefully; masses of individual objects are reproduced to better than 10% and positions to better than 100 kpc (Angulo & White, 2010b; Ruiz et al., 2011). In the following, we briefly recap the main features of the method.

First, consider a “target” cosmological model at $z=z_{B}$ which we seek to match using the results of an “original” cosmological simulation of side-length $L_{A}$ (in units of $h^{-1}$ Mpc). The heart of the method is to find a length transformation, $L_{A}\rightarrow L_{B}=s~L_{A}$, and a relabelling of the time variable $z_{A}\rightarrow z_{B}$, by requiring that the variance of the linear density field in the target cosmology, $\sigma_{B}^{2}(R,z_{B})$, over the range $[R_{1},R_{2}]$ is as close as possible to that of the original cosmology, $\sigma_{A}^{2}(R,z_{A})$, over the range $[s^{-1}R_{1},s^{-1}R_{2}]$ at redshift $z_{A}$. Thus, we minimise:

In the Press-Schechter theory (Press & Schechter, 1974) the halo mass function is determined by the linear variance of the underlying dark matter field, thus Eq. (1) also minimises the difference between the halo mass functions in the target and original cosmologies over the mass range $[M(R_{2}),M(R_{1})]$. We usually take $M(R_{2})$ to be the mass of the largest halo in the simulation at the lowest redshift of interest ($z=0$ here) and $M(R_{1})$ to be that of the least massive resolved halo.

As a result, the original box size will be expanded by a factor $s$, and the output at redshift $z_{A}$ will represent redshift $z_{B}$ in the target cosmology. Redshifts in the target cosmology ($z^{\prime}$) corresponding to the redshifts ($z$) of stored data in the original simulation are then determined implicitly by $D_{B}(z^{\prime})=[D_{A}(z)/D_{A}(z_{A})]D_{B}(z_{B})$. The mass of a simulation particle (in units of $M_{\odot}$) in the target cosmology is $m_{B}=(\Omega_{m,B}H_{B}^{2}/\Omega_{m,A}H_{A}^{2})\,s^{3}\,m_{A}$, where $\Omega_{m,X}$ and $H_{X}$ ($X=A$ or $B$) are the dimensionless total matter densities and Hubble parameters of the two cosmologies.

The second part of the algorithm of Angulo & White (2010b) corrects differences in power spectrum shape on large scales between the original and target cosmologies by altering the amplitude of quasi-linear modes using the Zel’dovich approximation. In this paper we do not make this correction since we are interested in the internal structure, abundance and evolution of massive haloes rather than in their spatial distribution, and we use only the original cosmology when looking at the overdensity around haloes (a quantity that is slightly affected by our scaling).

With this technique we have created two additional halo catalogues in alternative cosmologies which we denote M3 and M7. These have cosmological parameters motivated by the 3-year (Spergel et al., 2007) and 7-year (Komatsu et al., 2011) analyses of data from the WMAP satellite. The main feature of these models are values for $\sigma_{8}$ which are lower than in the MXXL and the other Millennium Simulations, and different values for $\Omega_{m}$ (see Table 1). The corresponding length scalings are $s=1.222$ and $s=1.072$ for M3 and M7, respectively. The $z=0.623$ and $z=0.319$ outputs of the MXXL represent $z^{\prime}=0$ in the M3 and M7 cosmologies.

In this paper we will assume that the QSO luminosity is a monotonically-increasing function of the FoF host halo mass at any given time, and that there is a duty cycle that is independent of the halo mass. Models with these characteristics appear to be preferred by clustering analyses (e.g White et al., 2008; Bonoli et al., 2010), but we note that they are not the only possibility. In physically motivated models the QSO host depends on the details of the galaxy formation model and in particular on AGN feedback (Marulli et al., 2008; Bonoli et al., 2009; Fanidakis et al., 2011, 2012).

Therefore, we consider halo samples limited by FoF mass at two different thresholds corresponding to two different abundances at $z\sim 6$ which we keep constant across our three cosmologies, “Long-lived QSO haloes” (LLQ haloes) have a comoving number density of $n=0.4\,{\rm Gpc}^{-3}$, requiring FoF masses above $[15.4,\,9,\,5]\times 10^{12}\,M_{\odot}$ for the M1, M7 and M3 cosmologies, respectively. The number density of this sample (which is well below the limit that could be probed by the MS) matches that of the extremely bright quasars observed at $z>5$ in the SDSS ($n=0.6\,{\rm Gpc}^{-3}$, Fan et al., 2003, 2006). Thus, they mimic the assembly history of SDSS QSOs if they have a $100\%$ duty cycle, i.e. if they shine constantly at their full brightness.

Our second sample, “Short lived QSO haloes” (SLQ haloes) is selected at 30 times higher abundance, i.e. $n=11.6\,{\rm Gpc}^{-3}$. This corresponds to minimum FoF halo masses of $[7.0,\,4.2,\,2.0]\times 10^{12}\,M_{\odot}$ for the M1, M7 and M3 cosmologies. This sample could be regarded as representing the hosts of the high-redshift SDSS quasars if these have a 1/30 duty cycle, i.e. each object shines at full brightness only 1/30 of the time.

Due to our length scaling, the total number of haloes in our two samples varies by factors of $s^{3}$ between the three cosmologies; the SLQ samples contain $810$, $997$ and $1478$ objects for the M1, M7 and M3 cases, respectively (note that we expect four objects above this mass threshold in the MS which, in fact, contained only two), whereas the LLQ samples contain $27$, $34$, and $50$ haloes. The redshift of the samples also differs between cosmologies because the MXXL data were stored only at a discrete set of times: we use $z=6.18$ (M1), $z=6.19$ (M3) and $z=5.96$ (M7) for the three cases.

Finally, we note that alternative sample definitions, for example, using virial masses or circular velocities rather than FoF masses, do not produce significant differences in our results. $70\%$ of the $1000$ most massive FoF haloes rank within the $1800$ most massive haloes according to $M_{200}$, and within the $4000$ according to peak circular velocity.

Since we are assuming that high-redshift quasars live in the most massive objects present at that time, one might naively expect that their descendants will lie at the centres of the most massive haloes at any later time, in particular, in the central galaxies of large galaxy clusters today. However, in this section we confirm the results of previous studies showing this is not generally true (Trenti et al., 2008; Overzier et al., 2009). The descendants of QSOs can be found in haloes spanning a factor of $30$ in mass, and in a few cases, they are not even located in the central object of this halo, but in an orbiting subhalo.

The upper panel of Fig. 2 shows differential mass functions for the two $z\sim 6$ QSO host halo samples (defined in Section 2.3) and for the three cosmological models of Table 1, which differ primarily in their value of $\sigma_{8}$ and $\Omega_{m}$. By construction, these samples consist of very massive haloes. Their mean mass ranges from $2.6$ to $8\times 10^{12}\,M_{\odot}$ for the SLQ samples, and from $5.2$ to $15.2\times 10^{12}\,M_{\odot}$ for the sparser LLQ samples – for both samples these mean masses increase smoothly with $\sigma_{8}$, as expected. Within each sample more than than $99\%$ of all haloes lie within a factor of two in mass of the threshold for inclusion. This is a consequence of the exponentially falling high-mass tail of the halo mass function.

The magnitude of the shifts between the three cosmologies are well described by the Jenkins et al. (2001) and Angulo et al. (2012) fitting formulae, which we display as solid or dashed lines in Fig. 2. However, at this redshift these formulae overpredict the number of haloes of a given mass by a factor of two to three. In part, this reflects the fact that these formulae were calibrated using simulations and redshifts where the most massive haloes were much less extreme than those considered here, but the disagreement might also be in part a consequence of the non-universal behaviour of the halo mass function (see e.g Tinker et al., 2008).

We identify the $z=0$ descendant of each halo in our samples as the object that contains the majority of its 15 most bound particles. Although in principle it is possible that a different structure contains most of the mass of our high-$z$ haloes, following the innermost particles should represent well the fate of a hypothetical black hole sitting at the centre of the halo. The bottom panel of Fig. 2 displays the FoF mass distributions of these descendant haloes. For comparison, we also present predictions made using the analytic merger tree algorithm of Parkinson et al. (2008) which is based on the Extended Press-Schechter (EPS) formalism and calibrated using merger trees extracted from the Millennium Simulation.

The median FoF mass of the descendants of our SLQ haloes varies by a factor of just $1.7$ across our three cosmologies (see the lower set of coloured arrows in Fig. 2). This is significantly smaller than the spread of a factor of $3.1$ in the the initial median masses. The typical descendant mass is $M\sim 5.6\times 10^{14}\,M_{\odot}$, corresponding to a moderately rich $z=0$ cluster. The median masses of the descendants of the sparser LLQ haloes are larger by a factor of about $1.5$ and also vary slightly more with $\sigma_{8}$ (see the upper set of coloured arrows in Fig. 2). The median initial masses of the LLQ haloes are typically a factor of $1.9$ larger than those of the SLQ haloes, so both within a single cosmology and between cosmologies the evolution is convergent in the sense that descendant masses are more similar than those of the original $z\sim 6$ objects. This is probably a result of QSO haloes, in all cases, descending into much less extreme peaks for which the differences among cosmologies are smaller. For example, the mass function for the M1 and M7 cosmologies are almost identical at $z=0$ for masses below $M\sim 10^{14}\,M_{\odot}$.

In all cases, there is a large spread in the masses of the descendants, as shown by the bottom panel of Fig. 2. This distribution at $z=0$ can be well described by a log normal with $\sigma=0.36$ dex. The most massive tail indeed corresponds to rare, high-mass clusters, but the least massive one corresponds to galaxy groups. The scatter is slightly smaller for the LLQ sample. We also would like to note that every descendant of an LLQ halo is the dominant structure of its $z=0$ group, but in the M1 cosmology $39$ out of the $810$ descendants of SLQ haloes are satellite subhaloes. For the M3 and M7 cases the corresponding numbers are similar. As we will see in the next subsection, all this reflects the variety of formation histories among dark matter haloes of similar mass.

Our distributions are in good qualitative agreement with the EPS-based calculations of Trenti et al. (2008). Using the same M1 cosmology, these authors found $68\%$ of the descendants of a sample analogous to our LLQ to have masses in the range $2.5$ to $12.2\times 10^{14}\,M_{\odot}$ with a median of $5.6\times 10^{14}\,M_{\odot}$. For this cosmology, our results show a similar scatter, but with a median which is $35\%$ larger. As can be seen in Fig. 2 for all cosmologies, similar scatter and median masses are predicted by the merger tree algorithm of Parkinson et al. (2008) if we use a halo sample that matches the number density of the SLQ haloes. This confirms that our scaling algorithm captures the main features of structure growth in different background cosmologies.

We now look in more detail at the paths connecting our samples of $z\sim 6$ quasar host haloes to their present-day descendants, concentrating on results in the unscaled MXXL, i.e. in the original M1 cosmology. In Fig. 3 we plot mass growth, defined as the ratio of descendant mass at each redshift to initial mass at $z\sim 6$. Light blue and brown regions indicate $86\%$ of the trajectories for SLQ and LLQ samples, respectively. Darker regions of each colour outline the regions containing 68% of the trajectories. Growth histories seem to be remarkably similar in the two samples and to be faster than exponential, described approximately by $\log_{10}M(z)/M(z=6.19)=0.21\times(6.19-z)^{1.2}$. This seems independent of initial mass at $z=6$, at least over the relatively restricted range of (high) masses considered here. Of course, this formula only describes the typical behaviour, and very different growth histories occur for haloes of similar initial mass. Some massive $z=6$ haloes have grown only by a factor of 20 to 30 by $z=0$ – much less than the change of the characteristic nonlinear mass-scale $M_{*}$ over the same redshift range – while others have increased their mass by factors of 200 to 300.

The spread in these trajectories increases substantially with decreasing redshift. At $z=3$ the rms scatter in the log of the fractional growth is $0.168$, at $z=1$ it is $0.271$ and at $z=0$ it is $0.33$, slightly larger than the initial separation in median mass between our SLQ and LLQ samples. Although we do not display them, the mass accretion histories of haloes in other cosmologies are very similar. This is, of course, expected.

Fig. 3 also shows the mass growth for the most massive (dotted) and least massive (dashed) haloes in our SLQ sample. By chance, the least massive halo is one of the fastest growing, with a fractional growth rate faster than that of the most massive halo at all redshifts. Indeed, this “low-mass” halo grows into a $2\times 10^{15}\,M_{\odot}$ object by $z=0$ , whereas the present-day descendant of the initially most massive halo is actually slightly smaller ($1.8\times 10^{15}\,M_{\odot}$). Neither of these haloes is in the extreme tail at $z=0$, ranking 2,224-th and 3,459-th in mass among MXXL haloes, and being four times less massive than the most extreme object. Curiously, none of the 20 most massive $z=0$ haloes in the MXXL has a progenitor in the SLQ sample.

Finally, the thick solid line in Fig. 3 indicates the mass growth that the most massive halo at $z=6$ would have to have in order to be the most massive halo at all redshifts. This is close to the actual trajectory of the largest halo until $z\sim 3$, but at later times, other objects take over the top spots.

All these examples are not exclusive to extremely massive haloes at high redshifts, but are a generic illustration that the most massive haloes at any given epoch will no longer be the most massive haloes at a later time. This is a consequence of the diverse assembly histories of haloes in a hierarchical universe.

We explore this behaviour directly in Fig. 4. The top panel indicates the fraction of SLQ halo descendants whose masses place them above three abundance thresholds. The bottom panel shows a complementary picture, showing the fraction of SLQ haloes that rise above various mass thresholds at later times.

At redshift 3, all descendants of SLQ haloes have masses above $2\times 10^{13}\,M_{\odot}$ and rank among the $3000$ most massive haloes per cubic Gigaparsec. (Recall that initially these objects are defined as the high-mass tail of the halo distribution with an abundance of 30 per cubic Gpc.) In contrast, only $5\%$ are in haloes with $M>1.4\times 10^{14}\,M_{\odot}$ and only $20\%$ still rank among the $30$ most massive haloes per cubic Gigaparsec. With time, the SLQ descendants fall further and further behind the most massive haloes in the MXXL. By the present day, only $2-5\%$ (depending on cosmology) would still be included in a mass-limited sample with $n=11.6\,{\rm Gpc^{-3}}$ and about $50-70\%$ would be included in a sample with 100 times greater abundance.

Why do some high-redshift haloes keep growing rapidly until the present, whereas others appear to shut down their accretion? Is this a random process or can it be related to some halo property at $z\sim 6$?

Fig 5 shows how the masses of the $z=0$ descendants of our SLQ and LLQ haloes correlate with their environment densities. The overdensity surrounding each high-redshift halo is computed on a variety of scales by mapping the underlying dark matter distribution onto a $2048^{3}$ grid and then convolving this density field with a Gaussian filter of different sizes.

Clearly, most of the QSO haloes live in overdense regions, but there is considerable diversity in how extreme these regions are. On small scales, all live in at least a $3\sigma$ region, but $7\sigma$ is the typical overdensity of an SLQ halo and $10\sigma$ that of an LLQ halo. With increasing smoothing, our quasar host haloes are found in progressively less extreme regions and the separation between the SLQ and LLQ haloes decreases. For a smoothing of $28$ Mpc the typical SLQ halo lives in a region of overdensity $1.4\sigma$, and 10% are located in regions of below-average density. Thus, our results suggest that it might be possible to find a quasar in the middle of a $14-28$ Mpc underdense region, even if quasars reside in the most massive haloes at $z\sim 6$.

Fig. 5 also shows that there is clearly a strong correlation between the overdensity around a high-$z$ halo and the mass of its $z=0$ descendant. In fact, this correlation is stronger than between the mass of the $z=6$ progenitor and that of its $z=0$ descendant (or with any other property in our catalogues). For example, if we consider the environment at $14$ Mpc, haloes that live in $<1\sigma$ regions end up in haloes of $M\sim 2.8\times 10^{14}\,M_{\odot}$, whereas those found in $>4\sigma$ regions typically end up in ten times more massive haloes ($2.8\times 10^{15}\,M_{\odot}$). The scatter in descendant mass at a fixed overdensity is $\sigma_{logM}=[0.061,\,0.044,\,0.045,\,0.063]$ for the fields smoothed on $[3.5,\,7,\,14,\,28]$ Mpc scales, respectively. These figures are to be compared with a scatter of $0.36$ for the sample as a whole. Thus, if the environment density surrounding a quasar is known, then the mass of its $z=0$ descendant can be predicted $6-8$ times more accurately than if it is not.

These results are easy to understand, since the masses of the $z=0$ descendants correspond to the mass contained within a sphere of radius about $7-14$ Mpc, and an object can only form by $z=0$ if its material is already overdense by a factor of about $1.68/D_{+}(z=6)$²²2$D_{+}(z)$ is the growth factor in units of its present-day value.$=0.32$ at $z=6$. Our findings also warn against a naive connection between objects at different redshifts – the linkage depends not only on the actual properties of the objects, but also on their environment.

Our $z\sim 6$ quasar host haloes are the most massive objects present at that time and so, by definition, are the objects which had the highest mean mass growth rates averaged over previous epochs. It is these extreme growth rates which must supply the baryons needed to build up the supermassive black holes and to fuel the extraordinarily luminous quasars which assume to lie in their cores. In this final section, we examine when and how fast our LLQ and SLQ haloes achieve their extreme masses.

We define an accretion history for each halo in our samples from the growth in mass along the main branch of its assembly tree. This branch is defined by stepping back in time from the $z\sim 6$, object, selecting the most massive progenitor of the current main branch object as the main branch object at the immediately preceding step. Note that with this algorithm the main progenitor at $z=10$, say, is not necessarily the most massive of all the $z=10$ progenitors. In fact, only about $\sim 40\%$ of our SLQ haloes have an identifiable main progenitor at $z\sim 10$ (i.e. with a mass above the resolution limit, 20 particles or $1.2\times 10^{11}M_{\odot}$) even though the MXXL contains 125,000 identifiable haloes at this time. Conversely, of the $100$ most massive MXXL haloes at $z=10$, only $30$ have a descendant among our SLQ sample at $z\sim 6$.

A consequence of these statistics is that SLQ haloes typically accrete most of their mass in a relatively short period before they are identified. This is illustrated explicitly in Fig. 6, which shows histograms of the $25$ and $50\%$ growth times of haloes in each of our three SLQ samples. These are defined for each halo as the times since it had a quarter and a half of its final mass, and they are given in units of the age of the universe, $t_{\rm H}$, at the time the samples were defined. Median values are $t_{\rm FT}/t_{\rm H}=0.89$ and $0.78$, respectively for our two definitions of formation time, and they are almost independent of the cosmological model. These values correspond to roughly $100$ and $200$ Myr prior $z\sim 6$.

Recent accretion rates are clearly substantially larger than the mean value required to grow the halo in the Hubble time. Median accretion rates for the last half of halo growth are $[9.7,\,8.9,\,7.6]\times 10^{3}\,M_{\odot}\,{\rm yr^{-1}}$ for the M1, M7 and M3 samples, respectively. The growth times of our LLQ haloes are similarly distributed to those shown in Fig. 6 and, as a result, their median accretion rates are two to three times higher. In all three cosmologies, the growth rates we find appear to be comfortably large enough to fuel even quasars as bright as the SDSS objects at $z\sim 6$, provided, of course that the associated baryons are able to shed most of their angular momentum and reach the central regions despite the tremendous luminosity being generated there.