The environment of bright QSOs at z ∼ 6: star-forming galaxies and X-ray emission

Abstract

INTRODUCTION

There is solid observational evidence that most galaxies host a supermassive black hole in their centre with properties that correlate with those of the galactic bulge (Dressler 1989; Kormendy & Richstone 1995; Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000). Currently inactive supermassive black holes are thought to have undergone a quasar (QSO) phase in the past, when they grew by accretion to their current sizes, emitting enormous amounts of radiation into their surroundings (Schmidt 1963; Salpeter 1964; Lynden-Bell 1969). Important clues to the formation of these enigmatic objects may thereby be provided by the very luminous QSOs at redshifts ≳6, revealed by deep observations such as the Sloan Digital Sky Survey (SDSS; Fan et al. 2001, 2003, 2004, 2006). These bright high-redshift QSOs are believed to be powered by supermassive black holes with masses of ∼10⁹ M_⊙. Some of these supermassive black holes appear to be even in place as early as z ≈ 7 (Mortlock et al. 2011). The existence of such massive black holes less than a Gyr after the big bang poses strong constraints on models in which supermassive black holes grow out of stellar mass seeds via Eddington-limited accretion (Volonteri, Haardt & Madau 2003; Volonteri & Rees 2005; Haiman 2006; Lodato & Natarajan 2006; Trenti & Stiavelli 2007). Models in which black holes grow from initially massive seeds are therefore perhaps the most promising in this regard. In this class of models, seed black holes of mass ∼10^{4 − 6} M_⊙ form via direct collapse of metal free gas on to protogalactic haloes at z > 10 (Haehnelt & Rees 1993; Umemura, Loeb & Turner 1993; Kauffmann & Haehnelt 2000; Oh & Haiman 2002; Bromm & Loeb 2003; Koushiappas, Bullock & Dekel 2004; Begelman, Volonteri & Rees 2006; Lodato & Natarajan 2006; Volonteri, Lodato & Natarajan 2008; Tanaka & Haiman 2009). Numerical simulations focusing on protogalactic haloes with virial temperatures just above the atomic cooling threshold at z ≈ 15 suggest that this scenario is indeed viable (Wise, Turk & Abel 2008; Regan & Haehnelt 2009a,b; Wolcott-Green, Haiman & Bryan 2011; Choi, Shlosman & Begelman 2013b; Latif et al. 2013; Prieto, Jimenez & Haiman 2013). Growth of supermassive black holes starting from massive seeds has also been modelled with full cosmological hydrodynamical simulations, where black hole sink particles of mass 10⁵–10⁶ h⁻¹ M_⊙ placed at the centre of haloes above a threshold mass of 10⁹–10¹¹ h⁻¹ M_⊙ evolve into a population of supermassive black holes with masses and accretion luminosities in line with observational estimates by z ∼ 6 (Sijacki, Springel & Haehnelt 2009; Khandai et al. 2012).

The clustering analysis of bright QSOs at intermediate redshifts, i.e. z ∼ 2–5, suggests that the black holes with mass ≳ 10⁹ M_⊙ are located at the centre of dark matter haloes with masses of ∼3 × 10¹²–10¹³ M_⊙ or higher (Porciani, Magliocchetti & Norberg 2004; Croom et al. 2005; Coil et al. 2007; Myers et al. 2007; Shen et al. 2007; da Ângela et al. 2008; Padmanabhan et al. 2009; Ross et al. 2009, but see also Trainor & Steidel 2012). At z ∼ 6 haloes that massive are still very rare and thus highly biased with respect to the overall matter distribution. Unless the fraction of the baryons finding their way into supermassive black holes is much higher than at lower redshifts, black holes with masses of ∼10⁹ M_⊙ observed at z ∼ 6 are expected to be surrounded by a wealth of other massive structures and should appear strongly clustered when compared to more ‘average regions’ of the cosmic density field. How such overdensities of massive haloes translate into an overdensity of bright galaxies surrounding the luminous QSOs at z ≳ 6 is, however, less clear (Kauffmann & Haehnelt 2002; Muñoz & Loeb 2008; Overzier et al. 2009; Romano-Diaz et al. 2011; Fanidakis et al. 2013).

Accurate direct mass measurements of the host galaxy/halo have so far proven difficult if not impossible. Some measurements of the linewidths of molecular gas have returned surprisingly low values (Walter et al. 2004; Wang et al. 2010, 2013; Willott, Omont & Bergeron 2013). It is, however, not clear to what extent these are actually representative of the full depth of the gravitational potential of the host halo and whether low-velocity widths should be interpreted as evidence that these supermassive black holes are located in significantly less massive haloes (Willott et al. 2005).

Studies of the environment of the observed high-redshift QSOs are therefore probably our best bet to find out what is the typical mass of the dark matter host haloes. None the less, observational studies aimed at verifying whether the number of galaxies surrounding high-redshift QSOs are consistent with the presence of a strong overdensity have been inconclusive so far. Significant overdensities were detected around two different SDSS QSOs at z ≈ 6 (Stiavelli et al. 2005; Zheng et al. 2006), taking the The Great Observatories Origins Deep Survey (GOODS) field (Giavalisco et al. 2004) as a reference average field. Kim et al. (2009) analysed a sample of i₇₇₅-dropout galaxy candidates identified in five Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) QSO fields which includes Stiavelli et al. (2005) sample, revealing two overdense and two underdense fields and one field with similar density compared to the mean density of the GOODS survey. Note, however, that with an angular surface area of 202 × 202 arcsec² Kim et al. (2009) only probed a region of ≈ 6 × 6 comoving Mpc² around the QSOs. Using the Subaru ‘Suprime-Cam’, with the larger field of view of 34 × 27 arcmin², Utsumi et al. (2010) (see also recent studies by Bañados et al. 2013; Husband et al. 2013) reported an overdensity of Lyman-break galaxies (LBGs) around another QSO at z ≈ 6.4. The quoted excess of LBGs occurs at a projected distance greater than 2 Mpc, suggesting the QSO that appears to be located in an underdensity on small scales might still reside in an overdensity, when larger scales are considered. Taken at face value, all of these observational results appear at odds with expectations from the Λ cold dark matter (ΛCDM) paradigm of galaxy formation, where pronounced overdensities are predicted around the massive haloes supposed to host the bright z ∼ 6 QSOs (Romano-Diaz et al. 2011). Using a mock survey of LBGs at z ≈ 6 constructed using the semi-analytic prescriptions applied to the Millennium simulation (Springel et al. 2005), Overzier et al. (2009) concluded that the 10¹²–10¹³ M_⊙ QSO host halo hypothesis for the z ∼ 6 QSOs is marginally consistent, but poorly constrained by available observational data due to the expected large scatter in the number of neighbouring LBGs.

Deeper data reaching to fainter magnitudes as well as improved modelling including relevant baryonic/feedback physics are clearly required for further progress. One of the main aims of this paper is to investigate this question for the first time with realistic mock images based on hydrodynamical cosmological simulations which include detailed supernovae and active galactic nucleus (AGN) feedback physics.

Observations of the environment of bright QSOs are also an important diagnostic tool for the study of the large-scale outflows driven by the momentum and energy released by the central engine. Strong coupling between these outflows and the interstellar medium of the host galaxy is the basis for most of the explanations of the correlations identified between the mass of supermassive black holes and the bulge properties (Haehnelt, Natarajan & Rees 1998; Silk & Rees 1998; Fabian 1999; King 2003; Fabian 2012). By expelling baryons and suppressing star formation, AGN feedback may also help explain the properties of local massive galaxies, such as their dearth in young stars and low gas content. Observational evidence for AGN feedback includes the detection of massive outflows in ultra-luminous infrared galaxies (ULIRGs; Sturm et al. 2011) and QSOs (Feruglio et al. 2010; Fischer et al. 2010; Rupke & Veilleux 2011; Cicone et al. 2012; Greene, Zakamska & Smith 2012). AGN-driven outflows have also been detected in high-redshift QSOs (Alexander et al. 2010; Nesvadba et al. 2011; Maiolino et al. 2012), for which there is some evidence of quenched star formation (Cano-Díaz et al. 2012). Typical outflow rates range from hundreds to thousands solar masses per year with wind velocities of the order of 1000 km s⁻¹. In the case of the z ≈ 6.4 SDSS QSO J1148 + 5251 (Maiolino et al. 2012), the inferred outflow rates are |${>} 3500\,\mathrm{M_{{\odot }}\,\mathrm{yr}^{-1}}$|⁠, even higher than the already large star formation rate (⁠|${\approx } 3000\,\mathrm{M_{{\odot }} \, \mathrm{yr}^{-1}}$|⁠) of the QSO host. The observed wind speeds are of the order of 1300 km s⁻¹ and are observed at a distance from the supermassive black hole as large as ≈16 kpc. A theoretical understanding of the driving mechanisms and properties of AGN-driven winds is starting to build up. A number of analytical and numerical works predicts AGN winds with speeds and outflow rates comparable to the observational findings (King, Zubovas & Power 2011; Debuhr, Quataert & Ma 2012; Faucher-Giguère & Quataert 2012; Roth et al. 2012; Sharma & Nath 2012). The spatial scale at which the energy of the wind is radiated away should depend strongly on the physical details of the driving of the wind. As the wind is launched from the AGN, it will collide with the surrounding medium driving both a shock through it and a reverse shock into the wind itself. Depending on the local gas properties, the reverse shock might cool on a time-scale which is shorter that the characteristic outflow time-scale, thus leading to a ‘momentum-driven’ flow. Alternatively, if the wind energy is not lost to the radiation but rather transferred to the surrounding medium, an ‘energy-driven’ flow develops. The presence and size of an extended cooling emission is therefore expected to be tightly linked to the detailed physics of AGN feedback as well as to the properties of gas surrounding the black hole. An additional aim of this paper is to investigate the predicted extended cooling emission around bright QSOs at z ∼ 6.

This paper is organized as follows. In Section 2, we describe the setup of our numerical simulations and lay down our method for constructing realistic mock galaxy catalogues. We provide the results of the analysis of our simulations in Section 3. Finally, we present our conclusions in Section 4.

METHOD

Numerical setup

In order to study the cosmological environments expected to host bright high-redshift QSOs, we selected the six most massive haloes at z = 6.2 from the dark matter-only Millennium simulation (Springel et al. 2005). This simulation follows cosmological structure formation in a box of 500 h⁻¹ comoving Mpc on a side. This is about a factor of 100 smaller volume than probed by high-redshift QSO surveys at z ∼ 6 and should thus be just sufficiently large to contain a few bright QSOs given their estimated spatial density (Fan et al. 2004). The descendants of these haloes all lie in rich clusters at z = 0 with masses in the range 10¹⁴ M_⊙–10¹⁵ M_⊙ (see Trenti, Santos & Stiavelli 2008; Angulo et al. 2012, for a detailed discussion of the descendants of massive high-redshift haloes). For comparison, we also selected six haloes at z = 6.2 that reside in large-scale regions with a matter density very close to the average cosmic density as well as six haloes in regions with a more moderate overdensity. As our average regions, we randomly selected six regions for which the density smoothed with a Gaussian filter with radius 5 h⁻¹ and 10 h⁻¹ comoving Mpc was within σ/2 at z = 6.2, where σ is the standard deviation of the respective Gaussian distribution. For our intermediate regions, we selected six haloes with mass in the range 4–8 × 10¹¹ M_⊙ at z = 6.2 located in regions which are intermediate in overdensity between our overdense and average regions.

We then generated high-resolution initial conditions at z = 127 centred on each of these 18 haloes and populated them with dark matter as well as gas. Beyond the high-resolution region, mass resolution was gradually degraded with distance, in order to maintain computational efficiency, while modelling large-scale tidal fields acting on the regions of interest correctly (for further details on re-simulation technique adopted see Sijacki et al. 2009). Numerical parameters of the simulations are summarized in Table 1.

The numerical setup of our simulations is similar to that of Sijacki et al. (2009), so we refer the reader to that paper for a detailed description. Here, we present a brief overview of the code and physical modules adopted and highlight the differences with respect to the Sijacki et al. (2009) setup. For most of our simulations, we employed the massively parallel treepm-sph code gadget-3 (Springel 2005). Dark matter is treated as a purely collisionless fluid and its evolution is completely determined by gravitational interactions, which are followed with a treepm algorithm. The code also tracks gravity, hydrodynamics and radiative cooling of a gaseous component, assumed to behave like an ideal optically thin mixture of hydrogen and helium. The gas is subjected to a spatially constant but time-dependent UV background with reionization occurring at z ≈ 6. Tabulated UV background intensities have been taken from Haardt & Madau (1996) as well as from Faucher-Giguère et al. (2009) (updated version from 2011 December). Additionally to the gadget simulations we have performed a number of simulations with the moving mesh code arepo (Springel 2010). arepo employs the same gravity solver as the gadget code and can be started from identical initial conditions. However, in order to evolve gas hydrodynamics, arepo uses a second-order accurate finite-volume method on an unstructured Voronoi mesh that is allowed to move together with the fluid. gadget and arepo share the same implementation for gas cooling and sub-grid star formation. The comparison of simulation results performed with these two codes therefore allows us to identify possible numerical inaccuracies stemming from the hydro solver (see e.g. Kereš et al. 2012; Sijacki et al. 2012; Vogelsberger et al. 2012) and to gauge their magnitude with respect to the uncertainties in the relevant physics that we are trying to simulate (for further details see Section 3.4.2).

Due to the inevitable lack of resolution, our simulations cannot follow the internal structure of star-forming regions and the formation of molecular hydrogen. We are thus forced to adopt a sub-grid prescription to account for the relevant physics of star formation and associated feedback. We here adopt the sub-grid model of Springel & Hernquist (2003). In this model, stellar particles interact purely gravitationally with other particles in the simulation and are spawned stochastically out of gas particles with densities exceeding a density threshold ρ₀. The dense star-forming gas follows an effective equation of state that accounts for the balance between supernova heating and associated cloud evaporation, gas cooling and cloud formation due to thermal instabilities. With respect to the Springel & Hernquist (2003) implementation, we slightly modify the star formation module in two ways: (i) for gas particles/cells which are hot and have densities above ρ₀, we estimate their new temperature as reduced due to radiative cooling losses; if such estimated temperature would fall below the effective temperature of the multiphase medium we set it equal to the effective temperature; (ii) we prevent star formation from gas that is hot (i.e. T > 10⁵ K) but above the density threshold ρ₀; this is done to prevent artificial spawning of star particles out of gas which e.g. was part of the multiphase medium but then heated to high temperature by AGN feedback. In some of the simulations, we additionally incorporate starburst-driven galactic outflows (as in Springel & Hernquist 2003), and we vary the mass loading of the winds in our various numerical experiments.

For our cosmological parameters, we assume a ΛCDM cosmology with parameters: h = 0.73, Ω_m = 0.25, Ω_Λ = 0.75, Ω_b = 0.041, identical to the parent Millennium simulation, but use a lower σ₈ value of 0.8, in better accordance with more recent results (Planck Collaboration et al. 2013). Note that this parameter sensitively affects the mass of the high-σ peaks where the first QSOs are thought to form and therefore plays an important role in this work. Since emphasis is placed on observational predictions or direct comparison with observations, all quantities are given in physical, rather than comoving, units throughout this paper, unless when explicitly stated otherwise.

Producing mock UV catalogues

One of the main goals of this work is to assess the likelihood of observing enhanced number counts of star-forming galaxies around bright z ∼ 6 QSOs compared to average regions of the Universe. We thus want to realistically ‘observe’ our simulation snapshots and compile galaxy catalogues that can be directly compared to HST observations. For this purpose, we construct simulated Hubble-like images, rather than relying on identifying dark matter haloes and sub-haloes and subsequently investigating the properties of the stellar populations residing within these sub-haloes. This direct procedure to construct projected luminosity maps has the advantage of being independent of the details of the halo and sub-halo finding algorithms, galaxy radius definitions, as well as to account for line-of-sight superposition.

We start from the snapshots of our simulations, which store information on the positions, masses, ages and metallicities of all the stellar particles. Using the Bruzual & Charlot (2003) models, we compute the rest-frame UV luminosity (at 1450 Å, corresponding to observed 1 μm at z ∼ 6) for each stellar particle that has been created, assuming a Salpeter initial mass function (IMF) for consistency with the chemical enrichment and supernova feedback implementation in our simulations. We then project the innermost cube of our high-resolution region (with a volume of 10³ h⁻³ comoving Mpc³) along the z-axis. At z ∼ 6, the resulting two-dimensional luminosity map corresponds to an observed field of view of approximately 6 × 6 arcmin² with a redshift depth of Δz = 0.032, which we construct on a grid with a cell-size of 0.1 arcsec. To convert the absolute luminosity into observed magnitude in the z band (for comparison with the Kim et al. 2009,HST observations), we adopt a distance modulus D = 48.9 and a K-correction K = −2.1 (flat f(ν) spectrum), in accordance with the cosmology we use throughout the paper. Finally, we add uncorrelated Gaussian noise in each pixel of the synthetic images corresponding to a 5σ sensitivity limit m_UV = 26.5 in an aperture of radius r = 0.2 arcsec (for comparison with Kim et al. 2009). We also consider higher S/N images which reach one magnitude deeper. Note that all magnitudes are given in the AB magnitude system (Oke & Gunn 1983).

To mimic galaxy detection in actual HST data, we then run the source finder code SExtractor (Bertin & Arnouts 1996) on our synthetic images, using detection threshold and deblending parameters typical for LBGs surveys at high redshifts (Bouwens et al. 2011; Trenti et al. 2011). As in Kim et al. (2009), we define the magnitude of a source as SExtractor's MAG_AUTO value, which represents the code's best estimate of the total luminosity and automatically includes aperture corrections.

This procedure closely resembles the construction of observational catalogues and provides us with a realistic sample of galaxies within our high-resolution computational volume. Minor discrepancies, which are, however, unlikely to significantly affect any conclusion drawn in this paper, are the use of Gaussian uncorrelated noise, while real data may include an extended noise tail (for example because of detector defects or residual cosmic rays), or correlated noise induced by the drizzling procedure (Casertano et al. 2000) commonly used to combine different HST exposures. Another difference with actual data is the fact that our projection only includes galaxies in a narrow redshift interval with Δz ∼ 0.032, so we are effectively neglecting crowding from foreground sources (which typically impacts ∼10 per cent of the area in medium-depth observations like (Kim et al. 2009), but which is accounted for in actual observations by means of artificial source recovery simulations).

In Fig. 1, we show how predictions for the apparent magnitudes of galaxies in one of our overdense regions compare with typical SFR-m_UV relations found in the literature. Assuming the SFR-m_UV relation found in Labbé et al. (2010), we construct mock images where we convert the star formation rates of each particle into a UV luminosity. We match galaxies between maps constructed in this way as well as according to our default method presented above. We then directly compare predictions of both methods. As shown in Fig. 1, the agreement is excellent, indicating that predictions for the UV luminosity of different galaxies should be robust as long as their star formation rates are numerically converged (see Section 3.4.1). Note also that UV bright z ∼ 6 galaxies are thought not to be significantly affected by dust attenuation (Bouwens et al. 2007; Labbé et al. 2010), which we therefore do not model in our mock catalogues.

Figure 1.

Open in new tabDownload slide

Apparent UV magnitude (AB magnitude system) obtained using the stellar population synthesis catalogues of Bruzual & Charlot (2003) (taking into account stellar age and metallicity) shown as a function of the apparent UV magnitude derived from a typical SFR-m_UV relation (here as in Labbé et al. 2010) for galaxies in one of our simulations of an overdense region at z = 6.2. The best-fitting line shows the agreement between the two to be excellent.

RESULTS

Black hole growth and cosmic environment

Fig. 2 shows the dark matter mass functions of all FoF haloes identified in the resimulated regions of our overdense and average samples. These were constructed by considering the number counts of FoF groups in a cube¹ of side length ≈1.9 Mpc (10 h⁻¹ comoving Mpc) per logarithmic mass bin.² In order to illustrate the variance in the mass functions from region to region, we shade the area between the highest and lowest mass functions for both average and overdense samples. The mass functions of the average regions are in excellent agreement with the theoretical mass function by Jenkins et al. (2001), confirming that our selected regions indeed provide a good representation of an average density field at this redshift. In the overdense regions, the mass functions reveal an excess of haloes at all masses when compared to the average regions and their high-mass end extends to higher values. This is in agreement with what is expected in the ΛCDM paradigm, where rare high-σ peaks become the first regions to host the most massive haloes at high redshift.

Figure 2.

Open in new tabDownload slide

Dark matter mass functions of the FoF groups for all overdense regions (upper band shaded in pink) and for all average regions (lower band shaded in blue) at z = 6.2. These were obtained by taking into account all FoF groups within a cube of side length ≈1.9 Mpc (10 h⁻¹ comoving Mpc) placed at the centre of the corresponding resimulation volumes. The dotted lines show the average mass function for the two cases. At the high-mass end, the variance in the mass functions increases due to the small volume of the high-resolution regions. The low-mass end is set by the minimum halo mass with at least 32 dark matter particles which for our simulation resolution corresponds to ≈2.96 × 10⁸ M_⊙. The mass function from Jenkins et al. (2001) at the same redshift is overplotted as the dashed curve for comparison.

A visual impression of our average density, intermediate density and highly overdense regions is shown in Figs 3, 4 and 5, respectively. The density field in the average density regions is characterized by thin filaments hosting low-mass haloes and large voids with relatively little structure. For the intermediate sample, the filaments become better defined and the voids appear less pronounced as the overdensity increases. In the overdense regions, the yet more prominent filaments intersect in central massive haloes with total masses M₂₀₀ ≈ 2–5 × 10¹² M_⊙. Note also that the gas density is increased on virtually all scales. The more diffuse gas in the highly overdense regions appears also denser than in the average or intermediate regions. The filled circles overplotted on the maps denote the positions of black hole particles, and their size is proportional to the black hole mass. The smallest circles correspond to a black hole mass of M_BH = 10⁶-10⁷ M_⊙ and the largest to M_BH ≥ 10⁹ M_⊙. The strong clustering of black holes in the overdense regions reflects the strong clustering of the underlying dark matter haloes which are massive enough to accumulate sufficient amounts of gas to facilitate efficient black hole growth. The most massive black holes, typically located in the most massive haloes at the intersection of the most prominent filaments (see also Dubois et al. 2012) of our overdense regions, have masses of ≈10⁸–10⁹ M_⊙. Even though our average density regions also contain black holes, they occur in significantly smaller numbers (in the case of the A4 region, none in the volume shown) and have much lower masses. While in the case of the A1, A4 and A6 regions there is one black hole with a mass of a few times 10⁷ M_⊙ (outside of the projected volume shown in Fig. 3), the most massive black holes forming in the average regions have grown to a typical mass of only a few million solar masses by z ∼ 6. Also note that among average, intermediate and overdense regions, there is considerable scatter in black hole number density due to cosmic variance. In the case of the O3 region, for example, the number of black holes is smaller than in the other overdense fields, while the central black hole is the most massive of our resimulated volumes at z ∼ 6. The growth of supermassive black holes appears therefore not to be only determined by the larger scale environment in which they grow, a point to which we will come back later.

Figure 3.

Open in new tabDownload slide

Mass-weighted maps of the gas density projected along a slice of thickness ≈1.9 Mpc (10 h⁻¹ comoving Mpc) in six regions of average density at z = 6.2. Black dots denote the location of black hole particles within the projected cube and their size is proportional to black hole mass with |$\boldsymbol\cdot$| representing black holes with mass in the range 10^{6 − 7} M_⊙ and •, black holes with mass ≥10⁹ M_⊙. Dark matter haloes (and hence black holes) form most efficiently along dense filaments that characterize the large-scale structure. The average density regions shown here do not, however, give rise to black holes more massive than a few times 10⁶ M_⊙. Note the substantial variations of the black hole abundance between the different regions due to cosmic variance.

Figure 4.

Open in new tabDownload slide

Mass-weighted maps of the gas density projected along a slice of thickness ≈1.9 Mpc (10 h⁻¹ comoving Mpc) in six regions of intermediate overdensity at z = 6.2. Black dots denote the location of black hole particles within the projected cube and their size is proportional to black hole mass with |$\boldsymbol\cdot$| representing black holes with mass in the range 10^6-7 M_⊙ and •, black holes with mass ≥10⁹ M_⊙. The typical black hole mass in these regions is a few times 10⁷ M_⊙, an order of magnitude higher than in the regions of average density. Note the substantial variations of the black hole abundance between the different regions due to cosmic variance.

Figure 5.

Open in new tabDownload slide

Mass-weighted maps of the gas density projected along a slice of thickness ≈1.9 Mpc (10 h⁻¹ comoving Mpc) in six overdense regions at z = 6.2. Black dots denote the location of black hole particles within the projected cube and their size is proportional to black hole mass with |$\boldsymbol\cdot$| representing black holes with mass in the range 10^{6 − 7} M_⊙ and •, black holes with mass ≥10⁹ M_⊙. The most massive black holes form in haloes with M₂₀₀ ≈ 4 × 10¹² M_⊙ at the intersection of very prominent filaments. Note the substantial variations of the black hole abundance between the different regions due to cosmic variance. The white square in the O6 region shows the field of view of HST/ACS.

In Table 2, we list the masses of the most massive black holes in each of our resimulated volumes as well as the total mass of their host dark matter halo. We also provide an estimate of the overdensity of each region, which we computed by simply taking the ratio of the average matter density to the average cosmic density inside a sphere of radius 5 h⁻¹ comoving Mpc and 8 h⁻¹ comoving Mpc with origin at the centre of the high-resolution region. Note that in the highly overdense regions the overdensities are already (mildly) non-linear even on scales of 8 h⁻¹ Mpc and that the ratio of these overdensities to the rms fluctuation amplitude σ ranges from 2.8 to 7.4.

A look at Table 2 further shows a correlation between overdensity and black hole mass, with more massive black holes generally being found in more overdense environments. However, scatter is significant. For instance, regions O2 and O3 contain black holes with masses M_BH ≈ 6.5 × 10⁸ M_⊙ and M_BH ≈ 1.2 × 10⁹ M_⊙, respectively. These masses are significantly higher than those found for the most massive black holes in regions O5 and O6, even though the former are less overdense. For the latter, the FoF groups that host these QSOs are also less massive than those of the most massive black holes in the O2 or the O3 region. Despite less overdense than other regions, we have verified that the QSO host halo in the O3 region had already collapsed at a higher redshift, assembling half of its z = 6.2 mass at z ∼ 8.5 rather than at z ∼ 7 as is the case for the QSO host haloes in the remaining five overdense regions. In the O2 region, instead, the QSO host halo experiences a major merger at z ∼ 7 which boosts both halo and central black hole mass. Thus, while a significant large-scale overdensity is required for the formation of sufficiently massive host haloes, the growth history of individual haloes also plays a role in determining the mass of massive black holes.

In our sample of regions with intermediate overdensity, dark matter haloes with masses in the range 3–7 × 10¹¹ M_⊙ form by z ≈ 6. Black holes forming in these haloes have masses 4 × 10⁷–10⁸ M_⊙, up to an order of magnitude higher than the mass of the most massive black holes found in our sample of regions of average density. Clearly, despite residing in haloes which are relatively rare with a comoving spatial density of ≲ 10⁻³ h⁻¹ Mpc⁻³, these black holes are not fed enough gas to acquire masses comparable to those powering the observed z ∼ 6 QSOs, even when growing from massive initial seeds. In our model, we thus require QSOs to be located in very massive haloes, that form early enough to provide a sufficiently deep potential well. Only under such special circumstances can gas be transported to the central regions of the halo in sufficient quantities to grow massive black holes that can power bright QSOs at z ∼ 6.

Mass growth and emission of the black hole population

Mass assembly histories

In this section, we investigate the mass assembly history and the evolution of the bolometric luminosity of the QSOs powered by the most massive black holes in our simulations. Coloured solid curves in the left-hand panel of Fig. 6 show the mass of the most massive black holes for our default simulations of overdense regions as a function of redshift. At z ≈ 15–14, dark matter haloes have just become massive enough to acquire a seed black hole at their centre. Soon thereafter, at z ≈ 13–12, the seed black holes enter a phase of exponential growth which proceeds until z ≈ 9–8, when the black holes reach masses ≳10⁸ M_⊙. In this period, most mass is gained through smooth accretion of gas. Major mergers would result in collisions of black holes of similar mass and a corresponding jump by a factor of about 2 in mass should be seen in Fig. 6. While this occurs (note for instance O2 at z ≈ 10.5), the largest contribution is due to Eddington-limited accretion resulting from cold gas brought in to the resolution limit by the large-scale filaments. Note that due to lack of resolution, we are unable to determine whether clumps formed via gravitational disc instabilities also feed the black hole (see also Bournaud et al. 2012; Di Matteo et al. 2012; Dubois et al. 2013).

Figure 6.

Open in new tabDownload slide

(a) Mass assembly history of the most massive black holes in the six highly overdense regions. Results for the O1 region are also shown for simulations including moderate and extreme winds. (b) The evolution of the bolometric luminosity of the QSOs powered by these black holes. Since the bolometric luminosity is proportional to the black hole accretion rate, the curves trace the accretion history of the most massive black holes. The growth of the black holes proceeds in three phases: (i) an initial phase at z ≳ 13 of gas supply limited accretion, in which the host haloes do not yet accumulate enough gas to reach Eddington-limited accretion, (ii) an Eddington-limited accretion phase in the redshift interval 13 ≳ z ≳ 9 and (iii) a feedback-limited phase in which accretion occurs in short Eddington-limited bursts followed by quiescent periods that result from AGN-driven outflows. We also show observed luminosities and inferred black hole masses taken from Willott et al. (2010) (CFHQS, triangles) and De Rosa et al. (2011) (SDSS, squares). Error bars are included when provided in the literature.

At z ≲ 8, black hole growth progresses more slowly. Short Eddington-limited bursts of accretion are interdispersed with longer phases of accretion well below the Eddington rate. By z ≈ 6, the most massive black holes in our overdense regions have reached masses in the range ≈3 × 10⁸–2 × 10⁹ M_⊙.

On the same plot, black symbols with error bars denote observational estimates from Canada-France-Hawaii High-z Quasar Survey (CFHQS) and SDSS for the mass of black holes powering the known QSOs at z ≈ 6 as provided by Willott et al. (2010) (triangles) and De Rosa et al. (2011) (squares),³ respectively. Overall, the masses of the most massive black holes in our default simulations (solid lines) are in reasonable agreement with the observational estimates, even though perhaps somewhat on the low side. Our default simulations do not produce black holes with mass substantially higher than 10⁹ M_⊙ by z ≈ 6 as seems to be the case for a few of the observed QSOs. If the duty cycle of these QSOs is indeed very high, this is possibly due to the fact that the Millennium simulation has an about a factor 100 smaller volume than probed by current QSO surveys. According to the observational data, there is also a population of QSOs powered by less massive black holes with masses as low as M_BH ≈ 6 × 10⁷ M_⊙ to a few 10⁸ M_⊙. These lower mass observed black holes may arise as a consequence of cosmic variance, but could also refer to black holes hosted by less extreme haloes found in regions with less pronounced overdensities. Such a population of black holes is predicted by our simulations, as can been seen from Figs 4, 5, 7 and also from Table 2 for our sample of intermediate regions.

Figure 7.

Open in new tabDownload slide

The bolometric luminosity of all black holes with masses M_BH > 10⁷ M_⊙ in our six overdense regions at z = 6.2 (red diamonds) and at z = 5.7 (blue diamonds) as a function of the black hole mass. Observational data was taken from Willott et al. (2010) (CFHQS, triangles) and De Rosa et al. (2011) (SDSS, squares). At z = 6.2, accretion on to the massive black holes shown in this plot is limited by thermal AGN feedback and therefore proceeds in short Eddington bursts, followed by more quiescent accretion periods. Thus, at any given time at this redshift, only a small fraction of QSOs will be at their luminosity peak and bright enough to have been detected in published surveys. Error bars are included when provided in the literature.

We have also investigated whether the inclusion of galactic outflows, which could be very important already at z ∼ 6 (Vogelsberger et al. 2013), impairs the growth of the black holes by clearing gas away from its sphere of influence. For this purpose, we used two additional simulations of region O1 which we performed with a prescription for moderate galactic winds, O1+winds, and extreme winds, O1+winds(strong), respectively. The aim of the latter was to maximize the possible impact of the effect of galactic winds (described in detail in Section 3.3.2). We found that galactic outflows can have the counter-intuitive effect of aiding black hole growth. Early black hole growth, shown as dashed and dotted lines (moderate and extreme winds, respectively) in Fig. 6, is delayed with respect to the default simulation O1 while the host halo potential is relatively shallow. When the halo becomes massive enough, winds are not able to efficiently eject gas from its inner regions, which then stays trapped within the halo and falls back to the central galaxy. Moreover, in smaller galaxies where star formation efficiency was reduced due to galactic winds a significant amount of baryons is expelled out of their dark matter haloes. This gas collides with the surrounding intergalactic medium and is shock heated. With time as the most massive halo in the resimulated region is growing in mass and the whole region around it is collapsing, some of this gas will be incorporated into the main halo boosting its baryon fraction and allowing for more gas to reach the central supermassive black hole. These effects lead to a factor of about 2 higher black hole mass than in the case without galactic outflows by z ∼ 6, still in good agreement with observational mass estimates. Note also, that even if all simulated black holes had a mass higher by a factor of 2, their mass would still fall short of M_BH ≈ 5–6 × 10⁹ M_⊙ as reported in observations at z ∼ 6 (De Rosa et al. 2011). As they boost black hole growth, galactic winds also indirectly increase the energy released due to the AGN feedback, leading to more powerful AGN outflows. In the simulation with extreme winds, O1+winds(strong), the mass of the most massive black hole is ∼2 × 10¹⁰ M_⊙ by z ∼ 6, exceeding the estimated most massive black hole at this redshift by a factor of 2. We do not regard this particular version of galactic outflows as the most realistic model and the effect of galactic winds in this simulation can probably be considered a (generous) upper limit winds of this type are likely to have on black hole growth.

Bolometric luminosities

In the right-hand panel of Fig. 6, we show the redshift evolution of the bolometric luminosities of QSOs powered by supermassive black holes in our simulations of the highly overdense regions shown in the left-hand panel. Since the bolometric luminosity of a QSO is related to the black hole accretion rate via |$L_{\rm BOL} \,=\, \epsilon _r \dot{M}_{\rm BH} c^2$|⁠, this plot may be equally interpreted as the accretion rate history of the black holes. During the initial growth phase, at z ≳ 13 − 12, accretion on to the black hole seeds is limited by gas supply. Since the mass of the initial seeds is rather high and the potential well of their hosts not yet deep enough, gas cannot be accumulated at a sufficiently high rate to power the black holes to their Eddington limit. This situation changes as more massive haloes assemble in the interval 13 ≳ z ≳ 9, when accretion becomes Eddington limited for all QSOs shown. At z ∼ 9, the accretion rates become high enough for the first QSOs to reach a bolometric luminosity exceeding L_BOL = 10⁴⁶ erg s⁻¹. The energy injected via thermal feedback starts to become comparable to the binding energy of the gas in the QSO host haloes, thus pushing gas away from the QSO and interrupting temporarily the Eddington-limited growth. Further gas accretion gives then rise to a cycle in which black hole growth alternates between short bursts of Eddington-limited accretion and longer periods of supply-limited accretion at lower rates. During this phase (see also Choi et al. 2012), the bolometric luminosities of the QSOs can vary by factors of 10–100 on times scales much shorter than the Hubble time. In the simulations with galactic winds, the phase of Eddington-limited accretion lasts longer and the black holes grow to larger masses. The corresponding bolometric luminosities are also significantly higher. AGN feedback eventually also starts to shut off the accretion periodically, introducing similar variability in the black hole accretion rate and hence QSO bolometric luminosity as in the simulations without galactic winds. The flux limit of current observational surveys of z ∼ 6 QSOs corresponds to L_BOL ≈ 2 × 10⁴⁶ erg s⁻¹, and only the QSOs that happen to be undergoing a strong accretion event (shown as peaks in the right-hand panel of Fig. 6) can be picked up in the present surveys. Using the estimates for the Eddington ratio of observed z ∼ 6 QSOs given in Willott et al. (2010) and De Rosa et al. (2011), we have computed bolometric luminosities by simply calculating λL_Edd(M_BH), where λ and M_BH are the provided estimates for Eddington ratio and black hole mass. These are respectively shown as triangle and square symbols in the right-hand panel of Fig. 6 as well. Estimates for the bolometric luminosities obtained in this way match the peak bolometric luminosities of the simulated QSOs reasonably well. The very brightest observed QSOs however require sustained accretion rates on to black holes which may only be possible in haloes which are too rare to be contained within the Millennium volume.

In Fig. 7, we show the instantaneous bolometric luminosity of the simulated QSOs in our default simulations as a function of the black hole mass at z = 6.2 (red diamonds) and z = 5.7 (blue diamonds). We include all black holes from our default simulations of overdense regions that have M_BH > 10⁷ M_⊙. At any given instant, only about one or two black holes have bolometric luminosities comparable to 10⁴⁶ erg s⁻¹ due to variability in the accretion rate in the feedback-limited regime. The remaining black holes, even those with masses exceeding 10⁹ M_⊙ have luminosities at z = 6.2 in the yet observationally inaccessible range 10⁴⁵-10⁴⁶ erg s^{− 1}. A key prediction of our simulations is therefore a steep rise in the number of observed QSOs with bolometric luminosities in the range 10⁴⁵ ≲ L_BOL ≲ 10⁴⁶ erg s⁻¹ and masses spanning the wide range 5 × 10⁷–2 × 10⁹ M_⊙. Note that the existence of such a population of QSOs at z ∼ 6 is not trivial, as it depends on assumptions made for the black hole seeds at such early times. We also predict the galactic hosts of these QSOs to show signs of strong AGN-driven outflows (see Section 3.5).

Eddington ratios

Fig. 8 shows the mean black hole accretion rates relative to the Eddington rate as a function of black hole mass. The top panels are for the six overdense regions (left), the six intermediate regions (middle) and the six average regions (right) at z = 7.9, while the bottom panels are for the same regions at z = 6.2. Note that results for O1+winds are also shown as the dashed orange histogram in the left-hand column. At z = 7.9, the highest Eddington ratios occur mainly for black holes with masses 10⁶ ≲ M_⊙ ≲ 10⁷ M_⊙. For overdense and intermediate regions, the accretion rate on to black holes in this mass range is Eddington limited. In the average regions, the accretion rates on to black holes in this mass range do not quite reach the Eddington limit due to a lower supply of gas. At z = 6.2, the spatial density of haloes hosting black holes with masses 10⁶ ≲ M_⊙ ≲ 10⁷ M_⊙ is higher and these haloes therefore no longer provide the rare potential wells which efficiently collect gas from a large region. As a consequence, typical accretion rates become lower. At both redshifts shown, efficient AGN feedback leads to a decrease in the typical Eddington ratios for black holes with mass ≳ 10⁸ M_⊙. For the overdense regions, where several very massive black holes can be found, the accretion rates are variable due to feedback-limited episodic accretion and the most massive black holes accrete at only a few per cent of their Eddington rate. Note how in our run with moderate galactic winds the accretion rates are suppressed for all the black holes but the most massive. The latter are the only black holes that reside in haloes which are massive enough to trap sufficient quantities of gas for efficient growth in their potential wells despite the influence of the galactic winds.

Figure 8.

Open in new tabDownload slide

Mean accretion rates relative to the Eddington rate as a function of black hole mass for overdense regions (left), intermediate regions (middle) and average regions (right) at redshifts: z = 7.9 (top) and z = 6.2 (bottom). The Eddington ratios are large for black holes with masses in the range 10⁶–10⁷ M_⊙ for all regions at both redshifts. At the higher redshift, accretion on to black holes in this mass range occurs very close to the Eddington limit. Note further that for the regions of average density, the lower gas supply means that on average black holes never quite reach Eddington-limited accretion. For the most massive black holes, with M_BH ≳ 10⁷ M_⊙, accretion is limited by strong thermal feedback and alternates between quiescent and Eddington-limited accretion episodes. The accretion rates of such massive black holes therefore vary considerably from region to region.

Star-forming galaxies

Number counts from mock images

Having discussed the evolution of the (most massive) black holes in our simulations, we now focus on the number counts of bright star-forming galaxies at z = 6.2. In order to identify galaxies in the different fields in a way as close to the observations as possible, we adopt the procedure outlined in Section 2.

An example of the spatial distribution of simulated mock galaxies is shown in Fig. 9, where the large-scale field is shown for the A1 average region (left), for the O1 overdense region (middle) and for the O1 overdense region with moderate galactic winds, O1+winds (right). We consider two magnitude limits for galaxies to be included in our catalogues: (i) m_UV ≤ 26.5 (as in Kim et al. 2009) and (ii) a selection going one magnitude deeper, m_UV ≤ 27.5. The galaxies identified by SExtractor as included in our m_UV ≤ 26.5 (m_UV ≤ 27.5) sample are marked by closed (open) black circles in Fig. 9, and in the second row they are visible as concentrated clumps against the noisy background. In the average region, galaxies appear evenly distributed within the field of view of ACS, with the exception of a small group of clustered galaxies. In the overdense region, galaxies clearly trace the dense filaments, where high gas density favours galaxy formation (see O1 in Fig. 5 for comparison). In the second row of Fig. 9, we zoom into the regions marked by the red rectangles, to show how they appear in our UV mock observational images. Apart from containing both brighter and more numerous sources, O1 (and O1+winds in the right-hand panel) also contain diffuse regions with significant UV emission, particularly surrounding the QSO host in the centre. Fainter galaxies whose apparent magnitude satisfies m_UV > 27.5 are also visible as lighter coloured clumps. The absence of such faint clumps in A1 indicates that the galaxy overdensity persists at lower magnitudes than considered with our magnitude selection, in agreement with the enhancement of the mass function at all resolved scales shown in Fig. 2.

Figure 9.

Open in new tabDownload slide

The first row shows the spatial distribution of mock galaxies within a cube with side length of 1.9 Mpc (10  h⁻¹ comoving Mpc) from one of our simulations of a region of average density, A1 (left), compared to one of our overdense regions, O1 (middle) and to a simulation of the same overdense region with galactic winds, O1+winds (right). We denote the galaxies, as identified by SExtractor with m_UV ≤ 26.5 (solid circles) and m_UV ≤ 27.5 (open circles). The large black square delimits the field of view of HST/ACS. The regions denoted by red rectangles – in the average region focusing on a small group of clustered galaxies, and in the overdense region centred on the QSO host – are expanded in the second row and shown as they appear in our mock images. As above, open black circles denote the identified galaxies with m_UV ≤ 27.5. In the overdense region, sources surrounding the QSO are brighter and far more numerous compared to the average region. Regions with diffuse UV emission are also present around the QSO host in the overdense region. The mock image of O1+winds shows that the number density of galaxies is very sensitive to galactic winds and it decreases by a factor of 3.7 (2.6) within the field of view of HST/ACS for m_UV ≤ 26.5 (m_UV ≤ 27.5). Note also the overdensity of galaxies drops quickly beyond the HST/ACS field of view.

In Table 3, we list the number of bright galaxies identified in the mock images of each of our simulations and for both of our magnitude selections. The excess of bright galaxies with m_UV ≤ 26.5 around our bright QSOs is quantified in the left-hand panel of Fig. 10. Here, we show the cumulative number of galaxies within a projected radius r for our different simulations. We shade the area enclosed between the minimum and maximum curves for both average and overdense regions in order to illustrate cosmic variance. Within a projected radius of 200 arcsec, the number of bright galaxies in our six average regions, as indicated by the blue band, lies within the range 1–5. For the six overdense regions, this number lies in the range 40–95 (pink band).

Figure 10.

Open in new tabDownload slide

(a) Cumulative number of galaxies within a projected radius r for galaxies of apparent brightness m_UV ≤ 26.5. The pink shaded band shows the cumulative number counts obtained for our overdense regions, and the dotted line passing through it is the mean of the six overdense regions. We also show the results for two simulations with modified UV backgrounds (green continuous and dashed lines; see Section 3.3.3 for further description) and for the simulations with supernova-driven winds of different strength (red continuous and dashed lines; see Section 3.3.2 for further details). The blue shaded band refers to the average regions, and the dotted line passing through it is the mean of the six average regions. Since for A1+winds, there is only one source in the region of interest, we mark the radius at which this source is located with a yellow square. Overdense regions typically contain an order of magnitude more galaxies than average regions. Galactic winds decrease the number counts in both overdense and average regions. (b) Luminosity function of detected galaxies down to m_UV = 27.5 (M_UV = −19.3) in the overdense regions (shaded pink) and in the average regions (shaded blue). Dotted lines passing through the pink and blue shaded regions of the plot denote the mean luminosity function of overdense and average regions, respectively. The luminosity functions for the simulations with a modified UV background is given by the green circles. Red circles denote the luminosity functions for the simulations including supernova-driven galactic outflows. Error bars illustrate the effect of Poissonian fluctuations in number counts in each bin shown. The inclusion of galactic outflows brings the luminosity function of one of our average regions in reasonable agreement with the Bouwens et al. (2007) luminosity function (shown as black line, together with 1σ error level as a light grey shade) and strongly decrease the expected overdensity of galaxies around the luminous QSOs in our simulations.

The right-hand panel of Fig. 10 compares the luminosity functions obtained from our mock images down to a magnitude m_UV = 27.5 to the observed (field) luminosity function by Bouwens et al. (2007). Note that despite matching the dark matter mass functions of average regions to Jenkins et al. (2001), our average regions contain significantly more bright galaxies than observed. This mismatch is not surprising given that we have not yet considered any strong stellar feedback processes e.g. due to reionization and/or supernova-driven winds. Especially the latter will be very efficient given the high star formation rates of the galaxies we are investigating.

Galactic outflows

The high star formation rates of bright high-redshift galaxies suggests that energy output from numerous supernova events is an important source of feedback. Supernova explosions drive galactic outflows that can drastically reduce the gas content of galaxies, especially from host haloes with shallow potential wells. To test this in our simulations, we performed two runs for which we switch on the Springel & Hernquist (2003) prescription for supernova-driven galactic winds. In this model, the gas outflow rate |$\dot{M}_{\rm out}$| is assumed to be proportional to the star formation rate |$\dot{M}_{\rm SFR}$| such that |$\dot{M}_{\rm out} = \eta \dot{M}_{\rm SFR}$|⁠, where the mass loading factor η is assumed constant for every halo. In these two runs, we left the wind speed constant at ≈480 km s⁻¹ and we only varied η, considering in turn the typical case where η = 2 (O1+winds) and a more extreme case where η = 10 (O1+winds(strong)), which should bracket the possible suppression of star formation due to supernova-driven winds.

In the right-hand panel of Fig. 11, we show how the spatial entropy distribution changes when supernova-driven outflows are included in our simulations. Additional entropy is generated due to dissipation of kinetic energy of the outflow leading to overall higher levels of entropy around star-forming regions. The mock image of O1+winds in the right-hand panel of Fig. 9 shows that the number density of galaxies is very sensitive to galactic winds and it decreases by a factor of 3.7 (2.6) within the field of view of HST/ACS for m_UV ≤ 26.5 (m_UV ≤ 27.5). Note also that the overdensity of galaxies drops quickly beyond the HST/ACS field of view. Fig. 10 shows the corresponding number counts of bright galaxies. Galactic winds affect the number counts in both overdense and average regions by a similar factor. Within 6.86 Mpc³ this factor is ∼5–6 for η = 2 and ∼17 for η = 10. In the O1+winds simulation, the number of galaxies falls from 50 (110) to 11 (43) for m_UV ≤ 26.5 (≤ 27.5). In O1+winds(strong), the winds are so efficient at ejecting gas out of dark matter haloes that only four bright galaxies remain within the projected volume. For the A1+winds simulation, the observed number counts drop from 6 (16) to one (four) and to zero for A1+winds(strong). We however caution that in our simulations with η = 10, we are intentionally considering an extreme scenario that represents very strong supernova-driven winds which are probably too strong to be realistic. In Fig. 12, we show that the inclusion of galactic winds suppresses star formation across the entire range of dark matter halo masses. As expected, this effect is more pronounced for O1+winds(strong). Note thereby that the star formation rate is suppressed comparatively mildly for the FoF group mass range 10¹²–10¹³ M_⊙. This is a consequence of extremely efficient suppression of star formation in lower mass haloes. The surplus gas is instead incorporated in the most massive haloes at a later time and accumulates there. With moderate galactic winds with η = 2, the luminosity function of the A1 region is only slightly higher than the observed luminosity function of Bouwens et al. (2007) (see the right-hand panel of Fig. 10). As A1 is the average region with the highest number of sources of our sample of average regions, it is likely that other average regions would have luminosity functions similar or slightly lower than that of Bouwens et al. (2007), leading to overall good agreement with observations (see also numerical studies by Finlator, Davé & Özel 2011; Jaacks et al. 2012). Note however that whereas good agreement with observations is achieved for global properties of simulated galaxies, higher resolution is necessary in order to investigate whether currently adopted sub-grid models lead to internal properties of simulated galaxies in line with observations.

Figure 11.

Open in new tabDownload slide

Mass-weighted maps of the entropy projected along a slice of thickness ≈1.9 Mpc (10 h⁻¹ comoving Mpc) at z = 6.2 in O1 (left), O1+UV (middle) and O1+winds (right) simulations. The entropy was estimated as S = T/ρ^2/3, where T and ρ are the gas temperature and density, respectively. The entropy distribution is strongly dependent on the input physics and is mainly generated in the virial shocks that form in very massive haloes (seen in all figures), but also by photoionization of diffuse IGM gas (middle panel) and by dissipation of kinetic energy in galactic outflows (right-hand panel).

Figure 12.

Open in new tabDownload slide

(a) The total star formation rate as a function of FoF group dark matter mass for the simulations of region O1 including different input physics at z = 6.2. Histograms shaded in pink give the star formation rates for our default simulation. Results from our simulations using different UV backgrounds are given by the green histograms; the solid line shows results for the default Faucher-Giguère et al. (2009) tables and the dashed-shaded histogram shows results for the same tables, but with the ionizing flux increased by a factor of 100. The solid and dashed orange histograms are for our simulations of O1 including moderate and strong galactic outflows, respectively. The biggest contribution to star formation in the default O1 simulation is provided by haloes with masses in the range 10^10-11 M_⊙. The onset of earlier reionization in the simulation with the Faucher-Giguère et al. (2009) UV background suppresses star formation in haloes with M_FoF < 10¹⁰ M_⊙. This effect is more pronounced for O1+UV(strong), where the additional gas that does not collapse on to lower mass haloes falls at lower redshifts into more massive haloes instead leading to a moderate increase in star formation for these haloes. Inclusion of galactic outflows leads to suppression of star formation in haloes of any mass and also shifts the mass range of the biggest contribution to star formation to haloes with masses in the range 10⁹⁻¹⁰ M_⊙. (b) Number of identified sources in different magnitude bins for the same simulations. Earlier reionization slightly alters the predicted magnitude of detected galaxies without changing their number counts. The inclusion of galactic outflows on the other hand clearly reduces the predicted number of galaxy detections.

We have also verified the impact of AGN feedback on the predicted number counts of galaxies by comparing to the simulation O1nobh, identical in setup to O1 but without black holes. Including thermal AGN feedback leads to a drop from 56 to 50 galaxies for the m_UV ≤ 26.5 mag limit. The number counts are almost unchanged when considering the number of galaxies down to m_UV = 27.5, dropping from 111 to 110. The small impact of AGN feedback on number counts can be attributed to the fact that only the most massive haloes (with masses above ∼ 10¹¹ M_⊙) hosting the brightest galaxies experience efficient enough black hole growth. Feedback-driven outflows in these galaxies reduce the star formation rates and consequently lower their UV luminosity. However, galaxies with mass ≲ 10¹¹ M_⊙ contain black holes with masses of at most ∼ 4 × 10⁶ M_⊙ and have therefore not released enough energy to expel gas from their interior. Since detectable galaxies reside in haloes with masses down to about 10¹⁰ M_⊙ and AGN feedback affects only the most massive objects, the number of galaxies in our sample is not significantly affected by this process.

Reionization history

We now investigate the effect of changing the reionization history. We have run two additional simulations which are identical to our default simulation setup apart from the UV background, which we have adopted from Faucher-Giguère et al. (2009) instead of Haardt & Madau (1996). The purpose of the two additional simulations ran with the Faucher-Giguère et al. (2009) UV tables was twofold: to test if a different reionization history, with the UV background becoming important at a higher redshift, or if a stronger ionizing background due to the higher number star-forming galaxies and QSOs expected in overdensities, can lead to lower number counts of bright galaxies. To obtain a lower limit for the number of galaxies in our overdense regions, we selected the region O1 since this is the field with the lowest count of galaxies to begin with. In the first test run, which we shall refer to as O1+UV, we replaced the UV background of Haardt & Madau (1996) with that of Faucher-Giguère et al. (2009). In the second test run, O1+UV(strong), we increase the ionizing flux by a (large) factor of 100, in order to reproduce the enhanced UV background due to the presence not only of a bright QSO but also brighter and more numerous galaxies. These very different choices for the reionization history should bracket the possible effects of local, UV bright sources on the star formation rates of simulated galaxies. Fig. 10 shows that the number counts of observed galaxies is not significantly affected by the choice of the reionization history. In order to understand this we have examined the time evolution of the neutral hydrogen fraction. As expected, the neutral hydrogen fraction drops at higher redshift in the O1+UV and O1+UV(strong) simulations, indicating that reionization occurs at higher redshift with the Faucher-Giguère et al. (2009) UV background, especially in the simulation where we boost the UV flux by a large factor. The projected mass-weighted entropy maps shown in Fig. 11 for the O1 (left) and O1+UV (middle) simulations at z = 6.2 show that the earlier onset of photo-heating alters the spatial entropy distribution considerably. In the O1 simulation, reionization occurs at the slightly lower redshift of z ≈ 6 and entropy is mostly produced within the virial shock of the central QSO halo. In the O1+UV simulation, however, the spatially uniform heat injection due to reionization leads to higher entropy in the entire volume.

The left-hand panel of Fig. 12 shows that an earlier reionization redshift leads to efficient suppression of star formation only for haloes residing in FoF groups with masses in the range 10⁸-10⁹ M_⊙. For more massive haloes, the effect of the heat injection during reionization is marginal. Galaxies hosted by these small mass haloes are not bright enough to meet our magnitude selection criterion and our predictions for the number counts are therefore not very sensitive to the heat injection during reionization. The right-hand panel of Fig. 12 shows that the suppression of star formation in the lowest mass haloes only somewhat lowers the prediction of the observed number of galaxies with magnitude in the range 27–27.5 for the O1+UV simulation. In the O1+UV(strong) run, star formation is suppressed by about an order of magnitude in haloes in the mass range 10⁸–10⁹ M_⊙. The excess gas which fails to collapse on to these low-mass haloes falls later on into the deeper potential wells of more massive haloes. The resulting boost in star formation rates leads to more UV emission and hence overall brighter galaxies, thus leading to the opposite effect than needed to suppress the number of galaxies. Thus, for the magnitude limits investigated here even a very strong ionizing UV background has very little effect on the predicted number counts of galaxies in QSO fields.

Predicted (over) densities of star-forming galaxies in ACS/HST-fields

We are now in a position to estimate the number of star-forming galaxies in the environment of bright z ∼ 6 QSOs as well as in average density regions. In order to estimate how many sources an HST/ACS observation employing the usual drop-out technique would be expected to detect, we consider a pencil beam with an area given by the field of view of ACS, i.e. 202 × 202 arcsec². Following Romano-Diaz et al. (2011), i.e. assuming a redshift range of Δz = 1, we get V ≈ 87.1 Mpc³ (11628 h⁻³ comoving Mpc³) for the volume probed by the Kim et al. (2009) observations. Due to incompleteness as well as contamination by foregrounds, the effective volume probed by these observations is likely to be equal to half of V (Romano-Diaz et al. 2011). From Table 3, we first compute the mean number of sources for our simulations without strong stellar feedback. Within 6.86 Mpc³ and for m_UV ≤ 26.5, we obtain N_avg = 3.5 ± 0.7 and N_QSO = 68.8 ± 10.6 for our samples of average density and highly overdense regions, respectively. The space density of galaxies is strongly enhanced by a factor of about 20 around the most massive haloes in the Millennium simulation. The expected number of sources in an average field as seen by HST/ACS, can be estimated as |$N_{\mathrm{ACS}} \,=\, 0.5 \times \frac{N_{\mathrm{avg}}}{V_{\mathrm{cube}}} \times V \, = \, 20.3 \pm 4.1$|⁠, where V_cube is the volume of the cube used for producing the mock images. As expected, this is significantly larger than the number of observed drop-out galaxies to the same magnitude limit in a field of this size (according to Kim et al. 2009, 4–8 per ACS field depending on the exact colour selection criterion), as we have used our simulations without strong stellar feedback for the calculation.

Repeating now the calculation for our simulations with galactic winds, we estimate for Kim et al. (2009)-like observations a number of sources of N_ACS = 5.8 for A1+winds in an ACS-sized field. Assuming that winds reduce the number of sources by a similar factor in all of our simulated average density regions, the predicted number of galaxies in an ACS-size field is N = 3.4 ± 0.7. This is in good agreement with the observed number which is quite sensitive to the details of the drop-out selection. To calculate the expected number of sources in fields containing one of our simulated bright z ∼ 6 QSOs, we assume that ACS observations probe a pencil beam volume made up of one of our highly overdense regions complemented with our average density regions so as to fill a redshift range of Δz = 1. As the excess of sources around the bright QSOs falls off rapidly with radius (see again the upper right-hand panel of Fig. 9) this should be a good approximation. Based on the A1+winds and O1+winds simulations, we estimate N_ACS = 10.1 for the number of star-forming galaxies in an ACS-size field. If we assume again that winds reduce the number of sources by a similar factor in the different simulated average density and highly overdense regions, we get N = 9.1 ± 2.1. This corresponds to a difference of about five sources in ACS – fields with bright z ∼ 6 QSOs compared to average density fields, in good agreement with the two observed overdense fields in Kim et al. (2009). Kim et al. (2009), however, also observed two underdense fields. We have not noted any significant reduction in the number of galaxies due to AGN feedback in any of our simulations and the occurrence of such underdense field would thus have to be attributed either to physical effects not included in our simulations or perhaps more likely to subtleties in the colour-selection criteria (physical or observational) or low-number statistics. These observational limitations arise primarily because of the relatively shallow depth of the Kim et al. (2009) observations. Since our simulations predict that the galaxy overdensity around the QSO haloes extends to fainter luminosities with an increasing number of sources, deeper imaging of the underdense fields identified by Kim et al. (2009) would discriminate between small number statistics fluctuations and/or scatter in the colour selection for LBGs. In addition, surveys with redshift information for the star-forming galaxies should – at least in principle – be able to provide stronger constraints. Obtaining redshifts at z ∼ 6 is, however, very challenging. Most interesting in this regard is perhaps a survey by Husband et al. (2013) at somewhat lower redshift which finds significant number of LBGs with redshifts within Δz = 0.05-0.1 of three bright z ∼ 5 QSOs. There have also been several attempts using Lyα emitters which, unfortunately, have been rather inconclusive, not least because of the difficulty of interpreting the presence or absence of Lyα emission at these redshifts (Bañados et al. 2013). In summary, our predictions for the overdensities around bright QSOs appear – apart from the reported occurrence of underdense regions in Kim et al. (2009) – to be consistent with observations, but there is no strong observational evidence yet that observed overdensities require the bright z ∼ 6 QSOs to be hosted in very massive haloes. The possibility that the host dark matter haloes have considerably lower masses and hence a lower number count of galaxies has been recently discussed by Fanidakis et al. (2013). In our simulations however, we have verified that haloes with viral masses in the range 4–8 × 10¹¹ M_⊙ only host M_BH = 3 × 10⁷–10⁸ M_⊙ instead of the required ∼10⁹ M_⊙. The formation of QSOs in lower mass haloes would therefore require a picture significantly different from the one in which supermassive black holes grow from massive seeds by essentially Eddington-limited accretion.

Numerical issues

Numerical resolution and comparison to semi-analytical models based on the Millennium simulation

In this section, we assess the dependence of the predicted number counts and star formation rates of bright galaxies on the resolution of our simulations and compare our results to semi-analytic models based on the Millennium simulation. Semi-analytic models are computationally cheaper and can therefore be used to obtain global properties of a sample of galaxies in very large volumes. Cosmological hydrodynamical simulations cannot presently simulate such large volumes at sufficient resolution. However, they treat gas hydrodynamics more self-consistently together with the interaction between baryons and dark matter. They also can more realistically follow gas inflows, outflows and non-linear processes occurring in phenomena such as mergers, turbulence and fragmentation. A more detailed representation of such processes then enables closer comparison with observations and subsequent tests of sub-grid physics. For this purpose, we compare our default O1 simulation, with a dark matter particle resolution of 6.75 × 10⁶ h⁻¹ M_⊙ and spatial resolution 1 h⁻¹ comoving kpc to a lower resolution simulation of O1lowres with dark matter particle mass m_DM = 8.44 × 10⁸ h^{− 1} M_⊙ and softening length ϵ = 5 h⁻¹ comoving kpc, parameters very similar to those of the Millennium simulation.

We ran subfind (Springel et al. 2001) in post-processing for snapshots at z = 6.2 for the O1 and O1lowres simulations in order to identify the gravitationally bound dark matter structures to which galaxies are more directly related. The resulting catalogue of sub-haloes in each of our simulations makes it possible to draw a fair comparison with galactic properties from the semi-analytic models of De Lucia & Blaizot (2007), where prescriptions for galaxy formation were applied to the dark matter sub-haloes identified in the same parent Millennium simulation and are available in an online data base.⁴ We focus here on resolved sub-haloes located within a spherical volume with radius ≈1.3 Mpc (7 h⁻¹ comoving Mpc), which is well inside our high-resolution region. In the left-hand panel of Fig. 13, we show the average star formation rate of sub-haloes for different sub-halo (dark matter) mass bins. The average star formation rates in both our low- and high-resolution simulations are in close agreement for M_sub-halo ≳ 10¹¹ M_⊙, where sub-haloes are well resolved. The agreement is poor for the lower mass sub-haloes, M_sub-halo ≲ 10¹¹ M_⊙. Due to low resolution, compact dense regions where star formation is particularly efficient are not accounted for, leading to star formation rates lower by a factor of about 4. Note that spurious angular momentum transfer, as can arise in cosmological simulations due to low resolution, (Okamoto et al. 2005; Kaufmann et al. 2007; Guedes et al. 2011) would normally lead to higher (and not lower) star formation rates as gas flows into the centre of haloes more efficiently. Moreover, as shown in Fig. 13, star formation rates are converged for haloes with mass ≳ 10¹¹ M_⊙ in the low-resolution simulations. Since our high-resolution simulations resolve haloes with a mass ∼ 10⁸ M_⊙, star formation rates in haloes hosting detectable galaxies, i.e. with mass ≳ 10¹⁰ M_⊙, should also be well within numerical convergence, indicating that numerical angular momentum transfer does not significantly affect our results. We have further verified that simulated galaxies have very flat rotation curves, with a ratio of peak circular velocities to circular velocity evaluated at the virial radius of ≲ 1.05. A detailed study of internal properties of simulated galaxies however requires even higher resolution simulations and is therefore well beyond the scope of this study.

Figure 13.

Open in new tabDownload slide

(a) The average star formation rate of sub-haloes as a function of sub-halo mass (dark matter only) for our O1, O1+winds, O1lowres and O1lowres+winds simulations as well as according to the semi-analytic catalogues of De Lucia & Blaizot (2007). Agreement between these three models is very good for sub-haloes with mass >10¹¹ M_⊙, but becomes poor for lower mass sub-haloes in O1lowres, O1lowres+winds simulations and in the semi-analytic catalogue. For low-mass haloes star formation is underestimated when compared to our high-resolution simulation, O1. (b) Cumulative plot of number counts as a function of apparent magnitude in O1, O1+winds, O1lowres and O1lowres+winds simulations. In O1lowres, haloes hosting observable galaxies are not well resolved leading to a drastic underprediction of sources with respect to O1. Incorporating galactic winds does not affect this conclusion.

In the same plot in Fig. 13, we also show the average star formation rates of sub-haloes for the corresponding galaxies according to De Lucia & Blaizot (2007). There is good agreement between the semi-analytic predictions and our numerical simulations, particularly with O1lowres, where the force resolution is the same as in the Millennium simulation on which De Lucia & Blaizot (2007) is based. The small differences between the O1lowres and De Lucia & Blaizot (2007) galaxy catalogues can be explained by somewhat different prescriptions for star formation and associated feedback processes. We note, however, that the semi-analytic galaxy simulations by De Lucia & Blaizot (2007) lack sufficient resolution to resolve dark matter haloes below ≈10¹⁰ M_⊙ and thus do not account for star-forming galaxies in these low-mass haloes.

In the right-hand panel of Fig. 13, we show the cumulative number count of galaxies as a function of apparent magnitude identified following the method outlined in Section 2.2 for the simulations O1 and O1lowres. A drastic reduction in galaxy numbers occurs by varying the resolution alone. By matching the IDs of the stellar particles situated within the pixels associated with each galaxy by SExtractor with their parent FoF group, we verified that the bulk of galaxies absent in O1lowres simulation resides in FoF groups of mass ∼10¹⁰ M_⊙ and in a few cases even in as low as a few times 10⁹ M_⊙. The fact that these haloes are only marginally, if at all, resolved in our low-resolution run leads to the discrepancy of one order of magnitude in the number of galaxies brighter than m_UV = 27.5. Note also a similar order of magnitude discrepancy for galaxies with magnitude m_UV ≤ 26.5. On the other hand, since our high-resolution simulations resolve haloes of mass ≈3 × 10⁸ M_⊙ which contain galaxies fainter than the 27.5 mag limit, our predictions for the number counts of bright galaxies are numerically converged.

Overzier et al. (2009) match the observed surface density of i-dropouts at z = 6 using the galaxy catalogues of De Lucia & Blaizot (2007). However, our results unambiguously indicate that haloes with masses down to ∼10¹⁰ M_⊙ can host detectable galaxies and must therefore be well resolved. Thus, similar semi-analytic studies based on higher resolution dark matter simulations would likely require stronger feedback in order to still match the observed surface density of i-dropouts at z = 6. Differences in galaxy numbers in overdense regions according to our simulations and Overzier et al. (2009) may also arise from the treatment of dust in their models (ignored in our simulations). However, the impact of dust in galaxies at z = 6 is thought to be low (Bouwens et al. 2007; Labbé et al. 2010) and likely small when compared to the order of magnitude effect of varying the resolution.

Hydro solver comparison

Different choices of hydro solvers used in cosmological simulation codes can give rise to significant differences in the properties of simulated galaxies and the intergalactic medium they are embedded in. Recent cosmological simulations using the moving mesh code arepo have highlighted a number of systematic differences in gas cooling and star formation rates, as well as in galaxy sizes and morphologies with respect to the identical simulations performed with the gadget code (Kereš et al. 2012; Sijacki et al. 2012; Torrey et al. 2012; Vogelsberger et al. 2012; Bird et al. 2013; Nelson et al. 2013). It is therefore important to verify that our predictions for the number counts of galaxies are robust against changes in the adopted hydro-solver scheme.

We have resimulated our O6 region with the moving mesh code arepo using the same setup as outlined in Section 2.1. Since our focus here is on comparing the global star formation rate and the number of observable galaxies, we have excluded black holes and AGN feedback from these simulations. Also by dropping the black hole module, we can more straightforwardly compare the two codes, given that gadget and arepo share the same implementation for gas cooling and sub-grid star formation, allowing us thus to isolate the possible differences that arise only due to the hydro solver. In the left-hand panel of Fig. 14, we show the total star formation rate within a cube of 1.9 Mpc (10 h⁻¹ comoving Mpc) on a side as a function of redshift for arepo and gadget runs (which we shall refer to as O6arepo and O6gadget, respectively). Even though at high redshifts, star formation rates are somewhat higher in arepo, from z ≲ 11 they are in very close agreement. The reason why we do not see any significant differences in the star formation rates predicted by these two codes is due to the fact that gas is in the rapid cooling regime and subject to an identical UV ionizing background. In fact, a similar result can be seen in the top panel of fig. 2 of Nelson et al. (2013), where distributions of the maximum past temperature of gas accreting on to different mass haloes are shown. Only when haloes are sufficiently massive and are at a sufficiently low redshift such that the quasi-hydrostatic hot atmosphere develops the differences between gadget's and arepo's hydro solver start to result in significant and systematic discrepancies. This finding also supports the view that at this very high redshift, our simulations should not significantly suffer from numerical overcooling or spurious angular momentum transfer.

Figure 14.

Open in new tabDownload slide

(a) The total star formation inside a cube of side length 1.9 Mpc (10 h⁻¹ comoving Mpc) centred on the most massive halo as a function of redshift for O6gadget (black filled circles) and O6arepo (red triangles). The close agreement between the predictions of both codes shows that the two differing hydro solvers converge on the global star formation rates at very high redshifts. (b) Cumulative number of sources at an apparent magnitude below a given value m_UV for O6gadget (black) and O6arepo (red). There is a systematic overproduction of galaxies at the faint end (m_UV ≳ 25.5) in the simulation performed using gadget. The excess of galaxies picked up in the gadget simulation also exist in the arepo simulation, albeit it is below the magnitude cut. This is a consequence of deficiencies in the gadget SPH solver, which lead to more compact galaxies and hence higher star formation surface density, leading to brighter objects at the faint end.

We then proceeded with the construction of the mock images for these two simulations at z = 6.2. We employed the same method as in Section 2.1, except that here for both runs we assumed the star formation rate and the UV luminosity to be related as in Labbé et al. (2010). The right-hand panel of Fig. 14 shows the cumulative number of galaxies below a given apparent magnitude in the UV band. While for the magnitude range from 23.6 to 25, O6arepo contains about five more bright galaxies, there are in total ≈1.2 more sources in O6gadget due to the differences at the faint end. The excess of bright galaxies in arepo is caused by more efficient gas cooling into haloes that already at z ∼ 6 are starting to become sufficiently massive so that the differences between the two hydro solvers start to show up. To understand why there are differences at the faint end of the luminosity function we matched galaxies between our gadget and arepo simulations. Considering the coordinates of a given object in one simulation, we checked whether a matching object exists within a radius of five pixels of the same position in the other simulation. If such a source exists, we then take it to be the same galaxy. Otherwise, we assume the given source to be missing if it is selected in one simulation but not in the other. This analysis revealed that there are no missing galaxies in arepo, but that they simply fall below the magnitude threshold of 27.5. Visually inspecting the galaxies just above this magnitude threshold in gadget and below in arepo reveals that these galaxies are generally less compact and more extended in arepo and thus appear fainter. While we attribute this difference to the deficiencies of the standard SPH solver, we note that all detected galaxies in gadget have a dark matter halo associated with them and are not artificial SPH ‘blobs’ (Sijacki et al. 2012; Torrey et al. 2012). While the discrepancies in the results produced by gadget and arepo as discussed above are of a numerical nature, we note that they are small compared to previous findings at lower redshifts. In particular, focusing on the total number of detectable galaxies around the bright QSO in our O6 region, our arepo simulation highlights the same problem we identified with gadget runs already – there is a need for very strong feedback processes to significantly suppress star formation rates in a few times 10⁹ M_⊙ to 10¹¹ M_⊙ haloes to alleviate the discrepancy between cosmological simulations and current observational findings.

Uncertainties in the sub-grid physics

Despite the high resolution of our simulations, we are still far from reaching the resolution necessary for resolving crucial baryonic processes relevant for galaxy formation ab initio. Low resolution means that we are forced to employ simplified sub-grid models to follow black hole seeding, accretion and associated AGN feedback as well as star formation and consequent outflows driven by supernova explosions. We address each of these uncertainties in this section.

In their cosmological simulations, Sijacki et al. (2009) investigated the effect of varying both the halo seeding threshold mass as well as the mass of the black hole seed (within the direct collapse model). They verified that raising the halo threshold mass for seeding by a factor of 10 over the value assumed in this paper leads to a black hole mass higher by a factor of ≈1.5 by z ∼ 6. This uncertainty however is smaller than variations in the black hole mass due to galactic winds, which can also increase black hole mass as shown in Section 3.2. Varying the black hole seed mass from 10⁵ h⁻¹ M_⊙ (our adopted value) to 10⁶ h⁻¹ M_⊙ was also shown to lead to only a slight increase in final black hole mass.

In a comparison of five popular black hole growth and feedback sub-grid models, Wurster & Thacker (2013) use simulations of major galaxy mergers to highlight the differences in the predicted black hole mass assembly history. Besides leading to differing accretion histories, the five models compared in Wurster & Thacker (2013) result in post-merger black hole masses that vary by a factor up to ∼7. The predicted black hole growth differs due to differences in the evolution of the density and temperature of the surrounding gas caused by the different AGN feedback prescriptions. In our simulations of overdense regions, gas inflows transport gas close to the black holes in sufficient quantities to power accretion at the Eddington limit until it becomes limited by feedback at z ≈ 8–9. Note, however, that the redshift where this is expected to happen will depend somewhat on the mass of the assumed black hole seeds as well as the mass of the haloes in which these are placed (Sijacki et al. 2009). Suggestions that the large gas supply reaches below our resolution limit during the exponential growth phase (Dubois et al. 2012, 2013) imply that Eddington accretion rates (as is the case in our simulations) should be robust and nearly independent on the underlying black hole model assumed. However, there remains a possibility that black hole growth progresses more slowly than suggested by our simulations, where the scales at which angular momentum of the inflows may slow down accretion are not accurately resolved. More physical recipes of black hole accretion applied to galaxy mergers, e.g. Hopkins & Quataert (2010) indeed suggest that accretion can progress less efficiently at scales below ∼10 pc. These simulations however do not capture the cosmological environment of QSOs, where large amounts of gas are brought to the vicinity of the accreting black hole, thus this question can only be settled using cosmological simulations with very high resolution. Finally, since black hole accretion rates are numerically converged (Sijacki et al. 2009), we expect spurious angular momentum transfer not to significantly affect black hole growth in our simulations.

The evolution after AGN feedback starts to shut off the accretion depends more sensitively on the details of the AGN feedback implementation. For example, Choi et al. (2012); Wurster & Thacker (2013); Choi et al. (2013a) give examples of how varying both input physics (such as X-ray heating, or winds driven by momentum rather than energy injection) and numerical implementation can lead to significant differences in accretion rate (see also Dubois et al. 2013) outside the regime in which gas is supplied very efficiently. We note that an important property of our feedback model is that it leads to self-regulated black hole mass growth. Thus, once the black hole growth becomes feedback limited, significant changes in the feedback prescription will still lead only to moderate changes in the black hole mass growth (for further details see e.g. Sijacki et al. 2009). None the less, the properties of AGN-driven outflows can be very sensitive to the specific feedback implementation, see e.g. Choi et al. (2012, 2013a) and Debuhr et al. (2010), and we will discuss this point in detail in Section 3.5.

Our cooling prescription may also be too simplistic as it relies on a cooling function that is only applicable to the primordial gas. As stellar mass builds up and heavy elements are produced and released into the surrounding interstellar medium, further cooling can occur via metal lines. Efficient metal line cooling would only exacerbate the discrepancy in galaxy counts between overdense and average regions discussed in Section 3.3. Vogelsberger et al. (2013) have recently performed cosmological simulations using the moving mesh code arepo, where the effects of metal-line cooling are investigated together with various feedback prescriptions (see also Schaye et al. 2010). Metal-line cooling and gas recycling are found to have little effect at high redshift, when gas is still relatively unpolluted by metals. We therefore expect metal line cooling not to affect any, but the most massive galaxies in our simulations of overdense regions at z ∼ 6.

Interestingly, Vogelsberger et al. (2013) however found gas cooling in lower mass haloes to be strongly limited by stellar feedback, which we have shown to play a major role in our simulations at high redshift as well. In this paper, such outflows have been shown to be a very promising means to prevent overcooling in galaxies already at z ∼ 6. In the adopted sub-grid model, we have assumed a very simple prescription in which the mass loading factor η is constant for every halo. There is, however, some evidence for a scaling between the mass loading and the depth of the potential well of the galaxy the wind is launched from. If the wind velocity is proportional to the galaxy's escape velocity v_w, then |$\eta \propto v_{\rm w}^{-1}$| for momentum-driven and |$\eta \propto v_{\rm w}^{-2}$| for energy-driven winds. The implementation of winds with such galaxy-dependent mass loading and wind speed has been recently tested in Puchwein & Springel (2013) (see also Vogelsberger et al. 2013), where good agreement with the low-mass end of the stellar mass function at z = 0 (and up to z = 2) is obtained using energy-driven outflows. Puchwein & Springel (2013) report star formation suppression due to galactic winds in haloes with masses up to ∼10¹² M_⊙, which at z = 6 completely encompasses the mass range of haloes of all the galaxies in the QSO fields as well as in average and intermediate regions. We therefore expect galactic winds to remain a viable mechanism to suppress the number counts of bright galaxies in both overdense and average regions at z ∼ 6. Different implementations may, however, affect our results quantitatively, by predicting a different shape particularly at the faint end of the luminosity function. Another feature of our sub-grid model for stellar feedback is that the emitted wind particles are temporarily hydrodynamically decoupled while they are still within the dense star-forming region. Thereafter, wind particles are hydrodynamically coupled again with the rest of the gas. Puchwein & Springel (2013) and Vogelsberger et al. (2013) show that such implementation of SN-driven winds leads to luminosity functions in very good agreement with observations. It is nevertheless possible that this implementation can affect the internal properties of individual galaxies such as their rotational velocity and morphology. Capturing these properties reliably would, however, require considerably higher resolution simulations.

The effects of (thermal) AGN feedback

AGN-driven outflows

In our implementation of AGN feedback, a fraction ϵ_f = 0.05 of the radiated bolometric luminosity is assumed to couple with the surrounding gas. This energy is distributed as thermal energy over gas particles located within a smoothing length of the black hole particle. For accretion at the Eddington limit, the energy is injected at a rate |$L_{\rm f} \approx 1.26 \epsilon _{\rm f} \times 10^{38} \frac{M_{\rm BH}}{{\rm M}_{{\odot }}} \, \mathrm{erg\, s^{-1}} \,=\, 6.3 \times 10^{45} \frac{M_{\rm BH}}{10^9 {\rm M}_{{\odot }}}\,\mathrm{erg\, s^{-1}}$|⁠. For one e-folding of a 10⁹ M_⊙ black hole, the total amount of energy injected corresponds to ≈9 × 10⁶⁰ erg. For comparison, the binding energy of gas within the most massive haloes in our overdense sample at z = 6.2 is typically only ∼10⁵⁹ erg.

The energy injected by the AGN is therefore more than sufficient for driving an outflow, the details of which depend on how much of this energy is dissipated by cooling. As we will see later, once the black hole mass has grown sufficiently, the energy is injected so efficiently that the gas in its environment is heated and overpressurized so rapidly that the expansion time-scale is faster than the cooling time-scale and the resulting pressure gradient leads to the acceleration of the surrounding gas, generating an ‘energy-driven outflow’ (King 2010; Ostriker et al. 2010; King et al. 2011; Choi et al. 2012; Zubovas & King 2012). The outflow transports energy to the outer regions of the halo, modifying the density and temperature distribution and therefore the spatially extended (Bremsstrahlung) emission from the surrounding gas. As we show later, the gas close to the QSO is heated by the AGN to temperatures exceeding 10⁸ K, well in excess of the virial temperature of the dark matter haloes hosting the luminous high-redshift QSOs in our simulations. Note that we have not tried to model cooling of the gas due to Compton scattering by the radiation emitted by the accreting black hole. Compton cooling due to the AGN emission should be important on scales smaller or comparable to our resolution limit and we assume here that this is accounted for by our choice of ϵ_f.

In this section, we investigate the physical characteristics of the AGN outflow and explore the properties of the spatially extended X-ray emission in our simulations. We mainly investigate the X-ray properties for all default simulations of overdense regions including the O1+winds and O1+winds(strong), but also compare to an additional simulation, in which region O6 is simulated without black holes. We shall refer to this simulation without AGN feedback as O6nobh.⁵

Temperature and density distribution

In Fig. 15, we show mass-weighted radial profiles of gas density and temperature of the most massive halo at the centre of the O6 and O6nobh simulations (black solid and black dashed lines, respectively) at z = 6.2 together with the corresponding profiles of the QSO host haloes from all other overdense regions. The coloured band indicates the softening length of the simulations. Densities, temperatures and radii are scaled to the virial gas density, the virial temperature and the virial radius of the QSO host halo of the corresponding simulation. The central gas density drops towards the centre in all simulations, except O6nobh and O1+wind(strong), where the density at the centre rises by about two orders of magnitude. Since the only difference between O6 and O6nobh is the presence of an accreting black hole and the corresponding thermal AGN feedback, this decrement is obviously attributable to the AGN heating which occurs in episodic Eddington-limited bursts (see Section 3.2). The reduction in density mirrors a corresponding increase in temperature within the same radius, also of about two orders of magnitude. As mentioned before, the AGN heats the gas to temperatures well in excess of 10⁸ K, while the typical virial temperature of these haloes is about ∼10⁷ K.

Figure 15.

Open in new tabDownload slide

Mass-weighted radial profiles of gas (a) density and (b) temperature at z = 6.2 shown for all overdense regions (coloured solid lines), O6nobh (dashed), O1+winds (dashed yellow) and O1+winds(strong) (dotted yellow). Densities, temperatures and radii are scaled to the virial density, temperature and radius of each respective simulation. The vertical yellow band marks the location of the gravitational softening length which has a fixed physical scale and thus varies somewhat when scaled to the virial radius of the different haloes. The central density is highest for the simulations O5 and O6, where the central supermassive black holes are correspondingly accreting at the highest rate in our sample of overdense regions (see fourth column in Table 4). Comparing O6 with O6nobh clearly shows that AGN feedback reduces the central density by one to two orders of magnitude. The decrease in density corresponds to the increase in temperature for all regions. For the case without AGN heating, gas remains cold in the central regions.

Outside the very central regions, the temperature profiles are pretty much flat, showing that the underlying gravitational potential is approximately isothermal. The only simulation for which the temperature decreases in the centre is O6nobh. Without AGN heating the gas at the centre of this halo cools efficiently due to its high density, leading to a core with radius ∼0.006r₂₀₀ ≈ 400 pc of cold gas. Note that in several cases, particularly for the models shown by the green (O3) and blue (O2) lines, the profiles do not extend all the way to the centre of the halo. This is also a consequence of AGN feedback which drives outflows and clears most of the gas away from the central regions of these haloes.

In order to further illustrate the differences due to thermal AGN feedback, we show a phase diagram of all gas particles with positive radial velocity in the halo rest frame, located within the virial radius of O6 in Fig. 16. We overplot contours enclosing the region of the phase diagram populated by more than 10 particles for the O6nobh (black) and O6 (white) simulations for comparison. The phase diagram can be divided into four regions, separated by dashed lines in Fig. 16 and each labelled with a different letter.

• This region of the phase diagram consists of shock-heated hot and diffuse gas. In O6, the amount of gas in this region is considerably larger and extends to both higher and lower temperatures than in O6nobh.

• The gas here is above the critical density to form stars, but is too hot to do so. It is above region D containing the star-forming gas lying on the effective equation of state. Most gas in this region was once lying along the equation of state but was heated by the central AGN. The gas in region B, and a portion of that in region A, is more abundant in the phase diagram of the O6 simulation than that of O6nobh as revealed by the contours. We expect the gas in this area of the phase diagram to constitute the outflow and we mark this region with a black rectangle.

• The gas in this region has cooled to T ≈ 10⁴ K, which is the lowest temperature to which gas can cool radiatively in our simulation.

• Gas in this region is star-forming and lies along the effective equation of state.

Figure 16.

Open in new tabDownload slide

Phase diagram of gas particles with positive radial velocity (in the central halo frame) within the virial radius of the QSO host halo at z = 6.2 for the O6 simulation (in colour). The colour encodes the number density of particles, with blue representing low and red high number densities. Contours enclose the region of the phase diagram with 10 or more particles for the O6 simulation (white) and O6nobh simulation (black). After the onset of strong AGN feedback, the crossover region between dense hot and diffuse hot gas (between regions A and B) becomes populated. This is consistent with a picture where the very dense cold gas in the vicinity of the QSO is heated by the thermal AGN feedback and moved away from the effective equation of state. As thermal energy is converted to kinetic energy, this gas escapes to regions of lower density. This region (marked with a black rectangle) therefore likely contains the bulk of the AGN-driven outflow.

Properties of the outflow

In Fig. 17, we show a histogram of the radial velocities of all outflowing particles contained within the virial radius of O6 (red histogram) and O6nobh (blue histogram). For O6nobh, the maximum outward (positive) radial velocity is about 1700 km s⁻¹. In contrast, thermal AGN feedback in O6 simulation gives rise to outflowing gas with outward radial velocities reaching values as high as 3000 km s⁻¹. The yellow histogram shown in Fig. 17 represents the SPH particles located within the black rectangle drawn on the phase diagram in Fig. 16, confirming that the bulk of the outflow is contained within this region. The AGN-driven wind thus mainly consists of hot gas that transports energy injected by the AGN away from the centre with velocities ≳1000 km s⁻¹.

Figure 17.

Open in new tabDownload slide

(a) Histograms showing the positive radial velocities of all gas particles within the virial radius of the most massive halo in the O6 simulation at z = 6.2. The red histogram shows the results for O6, while the blue histogram shows results for O6nobh. In yellow is the histogram for all particles contained in the box drawn in Fig. 16. The yellow histogram clearly occupies the tail of the O6 histograms (red). Comparison with Fig. 16 shows that the outflow consists of hot and relatively dense gas. (b) Cumulative plot of mass of all gas particles in simulation O6 moving with radial speed >v_rad enclosed by spheres with different radii as shown in the legend. The bulk of the outflow is contained within 10 per cent of the virial radius of the QSO host halo.

The right-hand panel of Fig. 17 shows a plot of the cumulative mass of outflowing gas contained within a sphere with four different radii surrounding the centre of the halo. The curves overlap at large velocities indicating that the fastest outflowing gas is located at r ≲ 0.1 r₂₀₀. At the halo's virial radius, the mass in gas outflowing with v ≥ 1000 km s⁻¹ is ≈5 × 10⁹ M_⊙.

Recently, Maiolino et al. (2012) have detected an outflow consisting of about 7 × 10⁹ M_⊙ of atomic gas from a QSO at z ≈ 6.4. The outflow is resolved at scales of ≈16 kpc and has a speed of ≈1300 km s⁻¹ as suggested by the width of the observed C ii transition line. The estimate of the gas mass is likely only a lower bound on the total mass of the outflow since some molecular gas is also expected to be present in the wind, suggesting that our simulated outflow may be somewhat weaker than the strongest AGN outflows observed so far. The feedback luminosity in our simulation O6 is, however, not particularly strong. As discussed in Section 3.2, observed luminous z ∼ 6 QSOs in the current flux-limited samples are generally caught at luminosities and accretion rates higher by a factor of about 30 or more than that predicted for the QSO in our z = 6.2 snapshot of the O6 simulations. We therefore note here that we expect observed luminous z ∼ 6 QSOs to have significantly stronger AGN-driven winds. Our simulation O1+winds(strong) which contains the most massive black holes at z = 6.2 with a mass of ≈ 1.3 × 10¹⁰ M_⊙ drives such a powerful wind with velocities of several thousand km s⁻¹ and outflowing gas reaching to hundreds of kpc. Inspection of Fig. 15 shows that the increased supply of gas driven out of smaller galaxies in O1+wind(strong) leads to a very high central gas density which efficiently radiates away the thermal AGN feedback energy resulting in a momentum- rather than energy-driven wind. This together with the increased depth of the central potential is responsible for the factor of 10 higher black hole mass required for the AGN feedback to limit further accretion. While we have argued that this exceptional mass growth may have been mediated by possibly unrealistically strong galactic winds, this may be an example of what kind of AGN-driven wind may be realized in even deeper and rarer potential wells than we have considered here.

In Fig. 18, we show mass-weighted maps of the projected radially outward velocity of outflowing gas for thin slices of thickness ≈57 kpc (0.3 h⁻¹ comoving Mpc) with origin at the centre of the halo for O6nobh (left) and O6 (right) simulations. The outflow in simulation O6 is clearly visible and strikingly non-spherical (see also Fig. 19), with gas flowing out with maximum speeds exceeding ∼1000 km s⁻¹. The white density contours show that the anisotropic outflow escapes along directions of least resistance into lower density regions. The outflow thereby avoids the dense filaments and sweeps across the voids to large distances. For O6nobh, shown on the left-hand side, positive radial velocities resulting from random motion of gas particles as well as pressure forces in the fluid are distinctly smaller than in O6.

Figure 18.

Open in new tabDownload slide

Mass-weighted maps of positive radial velocity projected along a thin slice of thickness 57 kpc (0.3 h⁻¹ comoving Mpc) for O6nobh (left) and O6 (right) simulations at z = 6.2. In O6 simulation, powerful outflows with speed ∼1000 km s⁻¹ are driven away from the supermassive black hole at the centre. For O6, we overplot plot a gas density contour level of ≈2 × 10⁵ M_⊙ kpc⁻³. The contours show that the AGN-driven outflow escapes preferentially into underdense regions, where it encounters lower ram pressure.

Figure 19.

Open in new tabDownload slide

Mass-weighted maps of positive radial velocity projected along a thin slice of thickness ≈57 kpc in four of our overdense regions at z = 6.2. White contours enclose regions where the gas density is higher than 10³ M_⊙ kpc⁻³. Powerful outflows with speeds ∼1000–1500 km s⁻¹ are driven by AGN feedback out of the centre of the most massive haloes in our overdense regions to distances up to ∼100 kpc from the QSO. The outflows are anisotropic because they preferentially propagate into lower density regions.

X-ray emission

Figure 20.

Open in new tabDownload slide

The X-Ray surface brightness around the QSO host halo at z = 6.2 in O6 simulation for a soft energy band 0.5–2 keV (left) and a hard band 2–10 keV (both in the observer's frame). Surface brightness contours of levels 10⁻²¹ and 10⁻²³ erg s⁻¹ cm⁻² arcsec⁻² are shown as solid lines for O6 and dashed lines for O6nobh. AGN feedback changes the distribution of density and temperature of gas as the AGN-driven wind flows to the outer regions of the QSO host halo, leading to more extended Bremsstrahlung emission. In the soft X-ray band, the emission is more spatially extended than it would be without AGN heating. In the hard band, the X-ray emission is dominated by rarer regions where the gas is extremely hot. The profile therefore appears more compact and its spatial extent does not vary significantly between O6 and O6nobh.

In the top row of Fig. 21, we compare radial profiles of the X-ray surface brightness and cumulative luminosities for O6 (solid lines) and O6nobh (dashed lines) for the soft and hard bands as well as the bolometric emission as shown in the legend. The bolometric luminosity L_Brems is higher in O6nobh than in O6 at almost all radii and particularly in the centre, where the gas density is higher for O6nobh due to the absence of strong outflows. In O6, almost all of the feedback energy is used to accelerate the energy-driven outflow and the strongly reduced density at the centre of the halo leads to an overall reduction of the Bremsstrahlung emission. In the observed X-ray bands, however, the thermal AGN feedback leads to a clear enhancement of the X-ray emission at small radii ∼400 pc due to the much hotter gas found in this region. The central peak in fact dominates the X-ray emission from gas in the proximity of the QSO as the cumulative plot shown at the right-hand side of the top row indicates. With an angular scale considerably smaller than 1 arcsec, this spatially extended X-ray emission would still appear as a point source even with the resolution of the Chandra X-ray observatory. In the bottom row of Fig. 21, we show the X-ray surface brightness profile and the cumulative X-ray luminosity in the observed soft band for all QSO host haloes in overdense regions (solid coloured lines), including results for simulation O1+winds (dashed yellow line). The central X-ray emission in the soft band is strongest for the QSO haloes in O1, O5 and O6 simulations. A look at Table 4, where we summarize properties of the supermassive black holes including QSO bolometric luminosity as well as properties of X-ray emission, reveals that these haloes are those where the largest amount of energy is being injected by the QSO. Comparison with the profiles shown in Fig. 15 shows these haloes also to have the highest central densities, which explains why the corresponding Eddington ratios are higher. Note that the higher central densities also promote cooling, raising the ratio of cooling losses to injected power, as shown in the last column of Table 4. In the cumulative plot at the right-hand side of the bottom row in Fig. 21, the QSO host haloes in the O1, O5 and O6 simulations therefore show the strongest central X-ray emission in the soft band and the corresponding AGN-driven winds loose significant amounts of the injected thermal energy. The emission is almost completely dominated by the central peak, especially in O5 and O6, where the Eddington ratio is ≳20 per cent. In other haloes, QSOs are experiencing comparatively quiescent accretion since gas was moved outwards due to AGN outflows, explaining the lower central densities found in these haloes (see Fig. 15). The redistribution of the hot central gas leads to moderately more extended emission. However, emission levels are up to two of orders of magnitude lower due to the absence of a dominant central peak of X-ray emission. Thus, either central densities are high, promoting high black hole accretion rates and hence high AGN energy injection, and cooling via Bremsstrahlung is boosted in the central regions or the central regions have been evacuated by an outflow and the central emission is less dominant and emission more extended. Including galactic winds further alters the X-ray emission profile, reinforcing the notion that X-ray emission is critically dependent on feedback physics. In O1+winds(strong), where we are probing the effect of extremely efficient accretion and outflows, dissipation of kinetic energy heats the dense gas at the centre so much that X-ray emission temporarily exceeds the injected power from the AGN. We note, however, that the results of O1+winds(strong) should be considered as a (generous) upper limit for the effect of thermal AGN feedback on the surrounding gas and its X-ray emission.

Figure 21.

Open in new tabDownload slide

Various properties of the X-ray emission in the QSO host haloes of our sample of overdense regions. The first row compares Bremsstrahlung surface brightness profiles (left-hand side) and cumulative Bremsstrahlung luminosity for different energy bands (observer's frame) for O6 (solid lines) and O6nobh (dashed lines) at z = 6.2. In the X-ray bands, O6 shows systematically (moderately) higher surface brightness than O6nobh at r ≳ 2 kpc. Due to the presence of a hot bubble at the centre, a large amount of radiation in the X-rays is emitted in the very central regions, in stark contrast with O6nobh, where no such feature is present. This central peak dominates the X-ray emission within the halo, as apparent from the cumulative luminosity on the right-hand side. The total Bremsstrahlung emission instead is higher for O6nobh due to the presence of large quantities of very dense gas at the centre. Due to the high redshift most of the energy released will be observed in the optical and UV bands. In the bottom row, we compare the X-ray surface brightness profile and the cumulative X-ray luminosity in the observed soft X-ray band for all different QSO host haloes, including our simulation of O1 with moderate galactic winds. The cumulative plot on the right illustrates a variation of about two orders of magnitude in total X-ray emission at any radii among the overdense sample. The brightest sources, O1, O5 and O6 are those where the QSOs are injecting the highest amounts of energy into their surroundings. For these regions, the central peak completely dominates the X-ray emission.

Inspection of Table 4 suggests that predicting observable effects of thermal AGN feedback on the X-ray emission from the host haloes of luminous z ∼ 6 QSOs and predicting prospects of discriminating this emission from the direct emission released during the accretion on to the central supermassive black hole will be difficult. AGN-driven winds generally lead to X-ray emission which is more extended than expected from hot halo gas alone. This is a non-trivial result, because it depends sensitively on where injected energy is radiated away and on the input physics of our simulations. At low accretion rates an appreciable fraction of the total emission in the soft band (0.5–2 keV) should be spatially extended, but our simulations predict flux levels too low for a detection with Chandra. With increasing accretion rate/feedback luminosity, the X-ray emission due to AGN feedback becomes more compact in our simulations and is probably detectable but not resolvable with Chandra. In this case, the spectral signature of thermal emission from hot gas with temperatures and luminosities significantly higher than expected from the cooling of gas in a halo without AGN feedback would be the diagnostic of choice. Our simulation O1+winds(strong) thereby suggests that an ultrahard spectrum due to very high temperatures may be such a signature to look out for. The predicted X-ray emission will, however, undoubtedly, depend also strongly on the details of the AGN feedback implementation. Indeed, Choi et al. (2012) show that injection of momentum instead of thermal energy can have a strong effect on both the predicted wind velocities and the X-ray emission. We therefore leave more definite statements to a future detailed study.

CONCLUSIONS

We have used a comprehensive suite of high-resolution hydrodynamical resimulations of the Millennium simulations to study the environment of luminous high-redshift (z ∼ 6) QSOs harbouring supermassive black holes with (several) billion solar masses. In our simulations, these black holes grow in the most massive dark matter haloes by Eddington-limited accretion from massive seed black holes and their growth is self-regulated by strong thermal AGN feedback. We have produced realistic mock images of star-forming galaxies and the spatially extended X-ray emission to investigate whether these probes of the environment of the observed luminous high-redshift QSOs can verify our picture for the growth of these early forming supermassive black holes.

Our simulation suite consists of cosmological hydrodynamical simulations of six highly overdense regions, six regions of moderate/intermediate overdensity and six regions of average density. We find that the central dark matter haloes in the highly overdense regions are not only more massive, but their environment is also markedly different with the presence of well-defined massive filaments intersecting the central dark matter halo. This position at the intersection of prominent filaments is instrumental for the continuous Eddington-limited fuelling of the central supermassive black hole which in turn is responsible for their very efficient growth. From z ∼ 8 black hole growth becomes AGN feedback limited and continues in short bursts of accretion at the Eddington limit which result in strong mostly energy-driven AGN winds interdispersed by longer phases of more moderate accretion. In contrast, in our simulations of average regions of the Universe, the massive seed black holes generally only grow to a few times 10⁶ M_⊙.

Our simulations predict that at the luminosities of the observed z ∼ 6 QSOs the luminosity function is very steep. The peak accretion luminosities of the supermassive black holes are in good agreement with the luminosities of observed bright high-redshift QSOs which are expected to be strongly biased towards large Eddington ratios. If AGN-driven winds indeed severely limit the gas supply already at z ∼ 8 then there should be large numbers of somewhat fainter QSOs with luminosities in the range 10^{45 − 46} erg s⁻¹ which are powered by similarly massive black holes accreting well below the Eddington limit. The masses of our supermassive black holes appear thereby to be somewhat lower than the highest values measured in observed bright QSOs at this redshift, albeit introducing (strong) galactic winds due to stellar feedback may raise the predicted masses of the most massive black holes. It is thus not clear whether the most massive black holes may have to be hosted in haloes more massive than the most massive haloes present in the Millennium simulation with its about a factor of 100 smaller volume than probed by current surveys.

We found that our predictions of the number of star-forming galaxies in the field as well as in the environment of luminous high-redshift QSOs at currently achievable flux limits are only weakly affected by photo(re-)ionization feedback, but are very sensitive to the presence and strength of supernova-driven galactic winds. We have further shown that previous studies of this kind based both on numerical simulations and semi-analytic modelling have suffered significantly from insufficient resolution. For efficient galactic wind implementations consistent with observation of galaxy properties at lower redshift, the predicted number of star-forming galaxies with small angular separation from luminous high-redshift QSOs is moderate and consistent with observations and the apparent failure to detect the large excess of star-forming galaxies which had been predicted by simple analytical models of QSOs that are hosted by the most massive dark matter haloes. The predicted number of excess sources increases, however, strongly for fainter flux limits. Deeper observations should thus be key for making progress with constraining the host halo masses of bright z ∼ 6 QSOs.

Once the black hole masses in the most massive haloes in our simulation have grown to a few billion solar masses, the (thermal) AGN feedback in our simulations strongly affects the overall spatial distribution and the thermal state of the gas in the haloes. The injection of 5 per cent of the bolometric luminosity of the AGN powered by the accreting supermassive black hole as thermal energy leads to ‘blast-waves’ which efficiently launch (mostly) energy-driven winds which move several billion solar masses of gas at speeds ≳1000 km s⁻¹ to distances of ≳10 kpc similar to what has been recently observed for AGN outflows, albeit mostly at lower redshift. We have created mock images of the predicted spatially extended X-ray emission from the environment of the luminous high-redshift QSOs. At low accretion rates/luminosities, our simulations predict a significant fraction of the X-ray emission to be spatially extended, but unfortunately at flux levels too low for a detection with Chandra. Spectral signatures of hot gas in the more compact and probably unresolvable emission at higher accretion rates/luminosities are therefore the most promising diagnostic of the effect of AGN feedback on the gas surrounding luminous high-redshift QSOs.

In summary, our simulations show that current observations of the environment of bright z ∼ 6 QSOs are fully consistent with the picture of Eddington-limited early growth of supermassive black holes from massive seed black holes with masses in the range 10⁵–10⁶ M_⊙. If massive seed black holes can indeed form at z ∼ 15 and bright z ∼ 6 QSOs are hosted in haloes as massive as in our simulations there appears to be no need for super-Eddington accretion. The observed number of star-forming galaxies suggests that efficient galactic winds are already operating at this high redshift. There is also the exciting prospect to detect not only the direct X-ray emission from the QSOs powered by accretion on to these supermassive black holes but also the effect of the thermal AGN feedback on the spatial extent and the spectrum of the X-ray emission. At these high redshifts, the detection of resolved spatially extended X-ray emission will be very challenging. At lower redshifts however, spatially extended X-ray emission around bright QSOs should be an important tool to study the mechanism of how the feedback from the energy released during the accretion on to supermassive black holes drives observed AGN outflows.

We thank George Becker for useful discussions and suggestions and Roberto Maiolino, Roderik Overzier and the anonymous referee for comments on the manuscript. We are also grateful to Simon White for granting access to Millennium Simulation data and Stefan Hilbert for assistance. All simulations presented in this paper were run using the dirac Complexity cluster based at the University of Leicester and the COSMOS Supercomputer based at the University of Cambridge funded by STFC. TC acknowledges an STFC studentship. MH acknowledges support by the FP7 ERC Grant Emergence-320596.

REFERENCES