SMBHs in EAGLE Universe

# Supermassive black holes in the EAGLE Universe. Revealing the observables of their growth.

Yetli Rosas-Guevara ; Richard G. Bower ; Joop Schaye ; Stuart McAlpine ; Claudio Dalla Vecchia ; Carlos S. Frenk ; Matthieu Schaller ; Tom Theuns .
DAS, University of Chile, Camino del Observatorio 1515, Las Condes, Santiago, Chile.
ICC, Physics Department, University of Durham, South Road, Durham DH1 3LE, UK
Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands
Instituto de Astrofísica de Canarias, C/ Vía Láctea s/n, 38205 La Laguna, Tenerife, Spain
Departamento de Astrofísica, Universidad de La Laguna, Av. del Astrofísico Francisco Sánchez s/n, 38206 La Laguna, Tenerife, Spain
E-mail: y.rosas@das.uchile.clE-mail: r.g.bower@dur.ac.uk
Accepted XXX. Received YYY; in original form ZZZ
###### Abstract

We investigate the evolution of supermassive black holes in the ‘Evolution and Assembly of GaLaxies and their Environments’ (EAGLE) cosmological hydrodynamic simulations. The largest of the EAGLE volumes covers a and includes state-of-the-art physical models for star formation and black hole growth that depend only on local gas properties. We focus on the black hole mass function, Eddington ratio distribution and the implied duty cycle of nuclear activity. The simulation is broadly consistent with observational constraints on these quantities. In order to make a more direct comparison with observational data, we calculate the soft and hard X-ray luminosity functions of the active galactic nuclei (AGN). Between redshifts and , the simulation is in agreement with data. At higher redshifts, the simulation tends to underpredict the luminosities of the brightest observed AGN. This may be due to the limited volume of the simulation, or a fundamental deficiency of the underlying model. It seems unlikely that additional unresolved variability can account for this difference. The simulation shows a similar ‘downsizing’ of the AGN population as seen in observational surveys.

###### keywords:
black hole physics, galaxies: formation, galaxies: active, methods: hydrodynamical simulations, quasars: general.
pubyear: 2016pagerange: Supermassive black holes in the EAGLE Universe. Revealing the observables of their growth.B

## 1 Introduction

Accreting supermassive black holes (SMBHs) are one of the most efficient sources of radiation in the Universe. During periods of strong activity, they are often prominent as optical nuclei of galaxies, and referred to as Active Galactic Nuclei (AGN). From a theoretical perspective, the vast energy outputs of AGN offer an appealing explanation for the steep cut-off of the massive end of the galaxy luminosity function (e.g. Bower et al. 2006, Croton et al. 2006) and the scaling of the X-ray properties of galaxy groups and clusters (e.g. Binney & Tabor 1995, Churazov et al. 2001, McCarthy et al. 2010).

From an observational perspective, the strong correlation between the mass of the central SMBH and the properties of the host galaxy, such as its velocity dispersion and bulge mass (see review by Kormendy & Ho 2013, also by Graham 2016), is consistent with a causal connection. One way to explore this connection, is to examine the evolution of AGN, for example by constructing luminosity functions at different cosmic epochs. Integrating the total energy radiated over the AGN lifetime then provides a method of charting the build-up of the rest-mass energy of SMBHs (Soltan, 1982).

Measurements of the luminosity distribution of AGN require large, unbiased samples selected over a wide range of redshifts and luminosities. Constructing such samples is difficult because a fraction of the emission that emerges from the SMBH is obscured by the surrounding gas and dust making an uncertain fraction of SMBH difficult to detect (Lansbury et al., 2015). Although spectroscopic optical surveys are able to scan wide areas and detect large numbers of AGN up to redshift , these surveys are biased to the brightest and most unobscured population of SMBHs. While the mid-infrared band can also be used to detect SMBHs via the reprocessed emission from dust heated by AGN activity, the emission from the SMBH is often overwhelmed by the host galaxy. X-rays therefore provide the most efficient and unbiased method of selection. Although soft X-rays are the most easily observed band, AGN selection is still biased due to gas extinction around the SMBH. This makes hard X-rays the least biased wavelength range to detect the full SMBH population as obscuration is greatly reduced. Recently, multiple large X-ray surveys have been carried out by *ueda2014, *aird2015 and *buchner2015. These studies have revealed that the AGN population evolves strongly and that their number density abruptly decreases between and today. Moreover, these studies show that there is strong ‘downsizing’ of the AGN population in the sense that the space density of higher-luminosity AGN peaks at higher redshifts.

Such deep X-ray surveys provide tests for models linking the build up of galaxies and their SMBHs. Recently, this has been explored using semi-analytic models where the growth of SMBHs and AGN feedback have been incorporated as analytic approximations. Typically, these models assume that AGN activity is triggered by major galaxy mergers or disc instabilities, and calibrate AGN feedback to reproduce the galaxy mass function (Bower et al., 2006; Croton et al., 2006; De Lucia & Blaizot, 2007). After accounting for strong dust obscuration of faint AGN, such models have been able to reproduce the observed AGN luminosity functions and AGN downsizing (Fanidakis et al., 2012; Hirschmann et al., 2012). Such studies rely on summarising complex hydrodynamic interactions by simple models, but provide important and useful approximations.

Hydrodynamic simulations offer an alternative approach, more clearly differentiating the resolved hydrodynamical interactions from the small-scale processes that cannot be directly resolved. AGN evolution has been explored by hydrodynamical simulations of isolated galaxy mergers (Springel et al., 2005a) and small cosmological volumes at high redshift. In these simulations, SMBH growth and AGN feedback are incorporated as subgrid physics (Springel et al., 2005a; Booth & Schaye, 2009). SMBH accretion is typically based on the pure Bondi-Hoyle model (Bondi & Hoyle, 1944) or on simple modifications of this (Springel et al., 2005a; Booth & Schaye, 2009). Recent studies, however, have recognised the importance of accounting for the effects of accreting high angular momentum material (e.g. Hopkins et al. 2009; Rosas-Guevara et al. 2015; Angles-Alcazar et al. 2016).

In addition, a fraction of the rest-mass energy of the accretion flow may be injected in the surrounding gas as thermal energy or a momentum driven wind or jet. Since such processes cannot be directly resolved, simulations choose to implement feedback in different ways, for example as thermal heating proportional to the mass accretion rate (e.g. Springel et al. 2005a, Booth & Schaye 2009), by explicitly distinguishing quasar and radio modes (e.g. Sijacki et al. 2007, Vogelsberger et al. 2014), or by injection of momentum into the surrounding gas (e.g. Power et al. 2010; Choi et al. 2013).

The latest generation of cosmological hydrodynamic simulations can now track the evolution of a galaxy population resolving the formation of individual galaxies with a resolution of within large cosmological volumes, typically on a side (Vogelsberger et al., 2014; Schaye et al., 2015; Khandai et al., 2015). In this paper we will focus on the Evolution and Assembly of GaLaxies and their Environments (EAGLE) simulations (Schaye et al., 2015; Crain et al., 2015). The EAGLE simulations reproduce many properties of galaxies, such as: the evolution of the galaxy mass functions (Furlong et al., 2015a), the evolution of galaxy sizes (Furlong et al., 2015b), the colour-magnitude diagram (Trayford et al., 2016) and the properties of molecular and atomic gas (Lagos et al., 2015; Bahé et al., 2016). Much of the success in reproducing the properties of the massive galaxies is due to the effects of AGN feedback (Crain et al., 2015). A key question is, therefore, whether the model reproduces the observables of SMBH evolution in a cosmological context. This test is almost entirely independent of the calibration procedure used to select model parameters, since the calibration procedure only considered the normalisation of the correlation between the present-day SMBH mass and galaxy stellar mass (Schaye et al., 2015).

In this paper we present the evolution of SMBHs in the EAGLE simulations from to . We compare the predicted X-ray observables in EAGLE to observational data from X-ray deep fields up to redshift 7. Such deep fields roughly correspond to the size of the largest EAGLE simulation, implying that we are restricted to densities of moderately luminous AGN. In section 2, we briefly outline the relevant subgrid physics used in the EAGLE project and describe how we compute the intrinsic X-ray emission from AGN using empirical corrections for the bolometric luminosity and the obscured fraction. In section 3, we present the results of the simulation. We summarize the properties of local SMBHs, such as their mass function in section 3.1. The evolution of the black hole mass function, the Eddington ratio distribution plane and the black hole mass-halo mass relation are investigated in section 3.2. In section 3.3 we compare the evolution of the AGN luminosity function in X-ray bands with the most recent observational estimates. We show that AGN in EAGLE follow a similar downsizing trend to that seen in observational data. Finally in section 4, we summarise and discuss our main results. Additional tests of simulation convergence and parameter dependencies are given in Appendix A and in Appendix B.

## 2 Code and Simulations

### 2.1 Code

In this study we use simulations from the EAGLE project. This consists of a large number of cosmological simulations, with variations in parameters, galaxy formation subgrid models and numerical resolutions, as well as a large, volume reference calculation. Full details of the EAGLE simulations can be found in Schaye et al. (2015) and Crain et al. (2015) (hereafter S15 and C15); here we give only a brief overview.

The EAGLE simulations were performed with a modified version of the parallel hydrodynamic code Gadget-3 which is a computationally efficient version of the public code Gadget-2 (Springel, 2005b). The improvements to the hydrodynamics solver, which are collectively referred to as Anarchy, aim to better model hydrodynamical instabilities, as described in Dalla Vecchia (in prep) (see also S15 and Schaller et al. 2015). Here, we concentrate on the reference model denoted as, Ref-L100N1504, which corresponds to a cubic volume of comoving Mpc () on a side. Initially, it employs particles. In order to study numerical weak convergence, we also use the simulations AGNdT9-L050N0752 and Recal-L025N0752 with box sizes and respectively, containing particles per simulation. Numerical weak convergence is defined in S15 and reflects the need of recalibrating the subgrid parameters to model more faithfully the physical processes at increasing resolution. Further simulation variations are considered in Appendix A.

The EAGLE simulations start from cosmological initial conditions at . The transfer function for the linear matter power spectrum was generated with CAMB (Lewis et al., 2000), adopting the Planck Cosmology parameters (Planck collaboration et al., 2013). The Gaussian initial conditions were generated using the linear matter power spectrum and the random phases from the public multiscale white noise Panphasia field (Jenkins, 2013). Particle displacements and velocities are calculated using second-order perturbation theory (Jenkins, 2010).

The setup of these simulations gives a mass resolution of for dark matter ( and for baryonic) particles. The gravitational interaction between particles is calculated using a Plummer potential with a softening length of comoving limited to a maximum physical size of . The box sizes, particle numbers and mass and spaced resolutions are summarised in Table 1.

### 2.2 Subgrid physics

The galaxy formation subgrid physics included in these simulations is largely based on that used for the OWLS project (Schaye et al. 2010, see also Crain et al. 2009). Many improvements have been implemented, in particular in the modelling of stellar feedback and black hole growth. We provide a brief overview below. Further details can be found in S15 and an extensive comparison of the effects of varying the subgrid physics parameters is given in C15. The values of the parameters that differ between the simulations can be found in Table 2.

• Radiative cooling and photoheating, star formation and stellar feedback.

Radiative cooling and photoheating are as described in *wiersma2009a. The radiative rates are computed element by element in the presence of the cosmic microwave background (CMB) and the UV and X-ray background radiation from quasars and galaxies (model of Haardt & Madau. 2001). Eleven elements are tracked. The radiative cooling and heating rates are computed with the software Cloudy (Ferland et al., 2013). Prior to reionization, the gas is in collisional ionization equilibrium and no ionizing background is present.

Star formation is implemented following the model of *scha08, including a metallicity dependent density threshold, (Schaye, 2004) above which gas particles are allowed to form stars. The model parameters are chosen to reproduces the empirical Schmidt-Kennicutt law which is encoded in terms of a pressure law. A temperature floor is imposed as a function of density, , for gas with . This value of leads to the Jeans mass, and the ratio of the Jeans length to the SPH kernel length, being independent of density, avoiding spurious fragmentation due to a lack of resolution. Gas particles are stochastically selected for star formation and converted to collisionless star particles. Each star particle represents a simple stellar population formed with a *chabrier2003 IMF.

Stellar evolution is implemented as described in S15 and *wiersma2009b. The stellar mass loss and consequent metal enrichment of 11 elements are modelled via three channels: (1) AGB stars, (2) Supernova (SNe) type Ia and (3) Massive stars and core collapse SNe. The mass loss of the stellar population, including metals, is added to the gas particles that are within an SPH kernel of the star particle.

Feedback from star formation is treated stochastically, using the thermal injection method described in *dallavecchia_schaye12. The total energy available to inject into the ISM per core SN is assumed to be . This amount of energy is injected 30 Myr after the birth of the star particle. Neighbouring gas particles are selected to be heated stochastically based on the available energy, and then heated by a fixed temperature difference of . The stochastic heating distributes the energy such that the cooling time relative to the sound crossing time across a resolution element allows the thermal energy to be converted to kinetic energy, limiting spurious losses. The fraction of the available energy injected into the ISM depends on the local gas metallicity and density.

• Black hole growth and AGN feedback.

Halos that become more massive than are seeded with black holes of ( for the simulation Small-seeds-L050N0752  presented in Appendix A) using the method of Springel et al. (2005a). In order to mimic dynamical friction, at each timestep the black holes less massive than 100 times the initial mass of the gas particles are relocated to the minimum of its local gravitational potential. SMBHs can then grow via gas accretion, where the accretion rates are calculated by the modified Bondi-Hoyle model presented in Rosas-Guevara et al. (2015):

 ˙Maccr=min(˙MBondi[C−1visc(cs/VΦ)3],˙MBondi), (1)

where is the sound speed, is the SPH-average circular speed of the gas around the black hole and is a viscosity parameter that controls the degree of modulation of the Bondi-Hoyle rate, , in high circulation flows. SMBH accretion rates are also Eddington limited. In contrast to Rosas-Guevara et al. (2015) and Booth & Schaye (2009), EAGLE accretion rates do not include an additional ‘-factor’ to boost the accretion rates when the surrounding gas density is high. This parameter is largely degenerate with the parameter. The values of in the simulations are found in Table 2. SMBHs also grow via mergers when they are within their smoothing length and have sufficiently small radial velocity. Further details are given in S15. Following Springel et al. (2005a), two masses are adopted for BH particles: a sub grid mass that is applied to the computation of the gas accretion rates and AGN feedback, and a particle mass that is used in the gravitational calculations. Initially, the sub-grid mass is smaller than the particle mass. Once the subgrid mass exceeds the particle mass, the SMBH accretes stochastically gas particles in its vicinity so both masses grow a step. This method ensures that sub-grid masses can be smaller than particle mass whilst conserving the gravitating mass.

AGN feedback is implemented following the stochastic model of *booth_schaye2009. Thermal energy is injected into the surrounding gas as a fraction of the rest mass energy of the gas accreted by the SMBH. Neighbouring gas particles of the SMBH are stochastically selected and heated by a temperature difference of for the simulation Ref-L100N1504 and for the simulation AGNdT9-L050N0752. The scheme is similar to that used to implement feedback from star formation, but uses a significantly higher heating temperature for the energy injection events. It is important to emphasise the simplicity of the feedback scheme that we adopt: a single mode of AGN feedback is implemented throughout using a fixed efficiency of , from which, a fraction of 0.15 is coupled to the surrounding gas.

### 2.3 Simulation outputs

Most published papers by the EAGLE collaboration are based on the analysis of 29 “snapshot” outputs, containing the full information on all particles at a particular redshift. These provide a good census of the masses of SMBHs at one particular time. Since we are interested in the dominant black hole of each dark matter halo, we do not use the quantities tabulated in the public database (McAlpine et al., 2015) as these correspond to summed quantities of all SMBHs within the halo. As discussed in that paper, these can differ significantly in the case of low-mass black holes.

Although “snapshot” outputs can be used to construct the AGN luminosity function, the strong variability of AGN in the simulation means that the statistics of luminous AGN are poorly sampled. To obtain a better determination of the AGN luminosity functions, we make use of the more frequent “snipshot” outputs. These are partial copies of the particle state of the simulation, which are output in order to track critical simulation quantities with higher time resolution. There are 406 (400, 406) ‘snipshots’ outputs for the Ref-L100N1504 (AGNdT9-L050N0752, Recal-L025N0752) simulation, with a temporal separation between and Myr. In the simulation, AGN are highly variable on significantly shorter timescales, and we average the luminosity functions in ranges of snipshots in order to improve the statistical sampling of luminous outbursts. A detailed analysis of AGN variability will be presented in McAlpine et al. (2016, in preparation). Although this procedure allows us to reduce sampling uncertainties due to variability, it does not allow us to include rare objects that are not present in the simulation volume. In appendix A, we compare simulations with different volumes, but the same parameters. The analysis presented then suggests that variability is the dominant uncertainty and that the procedure we use does not appear to cause a significant underestimate of the abundance of luminous AGN.

To give an impression of the size of the SMBH population, the first SMBH appears in the Ref-L100N1504 simulation at , the most massive black hole has a mass of and is located in the most massive halo (which has a mass of ) at . The total number of SMBHs at with mass larger than is of which have masses , have masses and have masses .

### 2.4 Post-processing and definition of accretion regimes

As we have previously stressed, the subgrid models of SMBH accretion and feedback used do not make any distinction between different regimes of SMBH accretion. In post-processing, we distinguish between the activity levels of SMBHs based on their Eddington ratio,

 λEdd≡˙Maccr/˙MEdd, (2)

where and are the SMBH mass accretion rate and the Eddington limit respectively. We define two ‘active’ accretion regimes. For Eddington ratios larger than , we assume that the nuclear disc around the SMBH is thin and radiative cooling is efficient. We therefore assume that the luminosity of the disc can be described by the standard Shakura-Sunyaev disc model (Shakura & Syunyaev 1973). Such sources will be highly luminous in X-rays. We will refer to SMBHs in this regime as Shakura-Sunyaev discs (SSDs). For in the range , we assume that the nuclear accretion disc is thick and radiatively inefficient. We will refer to these SMBH as Advection Dominated Accretion Flows (ADAFs) (Rees et al., 1982; Narayan & Yi, 1994; Abramowicz et al., 1995). By default, we assume that these sources make negligible contributions to the X-ray luminosity function. Such sources are, however, expected to dominate radio source counts. Finally, we classify those with as inactive and assume that such sources are essentially undetectable against the emission of the host galaxy.

Note that the choice of a threshold in the Eddington ratio to define both active states SSDs and ADAFs does not have a significant effect on the X-ray AGN luminosity functions as shown in Appendix B.

### 2.5 Predicting X-ray observables

In this section, we describe the method used to predict X-ray luminosities from the SMBH accretion rates. We consider only AGN in the SSD regime (). For such sources, the bolometric luminosity is

 Lbol=ϵr1−ϵr˙MBHc2=ϵr˙Maccrc2, (3)

where is the radiative efficiency and is set to as suggested by Shakura & Syunyaev (1973).

We use the redshift independent bolometric corrections of Marconi et al. (2004) to convert the bolometric luminosity into intrinsic hard (2–10 keV) and soft (0.5-2 keV) X-ray band luminosities. The bolometric corrections are third degree polynomial relations defined as follows:

 log10(LHXLbol) = −1.54−0.24L−0.012L2+0.0015L3 log10(LSXLbol) = −1.64−0.22L−0.012L2+0.0015L3, (4)

where . The bolometric corrections are computed with a template spectrum that is truncated at to exclude the IR bump (Marconi et al., 2004) produced by reprocessed UV radiation. The correction is assumed to be independent of redshift. We note that Vasudevan & Fabian (2009) have suggested that the bolometric corrections may be a function of the Eddington ratio, but the differences are not significant except in AGN with . Lusso et al. (2012) suggested that the bolometric corrections could be lower than those of Marconi et al. (2004) at high bolometric luminosities, but the offset is small in the context of this work. Hopkins et al. (2007) also proposed expressions for the bolometric corrections based on dust absorbed luminosities. For us, this is inappropriate since we base our analysis on intrinsic X-ray luminosities. Thus, we opt for the relations of Marconi et al. (2004) which has the benefit of being consistent with previous studies. Hereafter, we always refer to luminosity () as the intrinsic luminosity in the 2–10 kev (0.5–2 keV) rest-frame energy range.

The emission of AGN may be absorbed if the circumnuclear environment is rich in gas and dust. The absorption is likely to be highly anisotropic, making a fraction of sources undetectable in the soft X-ray band. From the observational data, it is unclear whether the obscured fraction is a function of redshift and other sample properties. For example, early studies e.g. Ueda et al. 2003; Steffen et al. 2003 did not find clear evidence for a redshift dependence, but recent studies have established that the obscured AGN fraction increases with increasing redshift e.g. Hasinger 2008;Treister et al. 2009;Ueda et al. 2014. Because of these uncertainties, we prefer to compare the simulations to observed soft X-ray luminosity functions for which the obscured fraction has already been taken into account by simultaneously fitting to both hard and soft X-ray data (Aird et al., 2015). In future work, we will investigate the obscuration of AGN due to gas and dust, taking into account the properties of the host galaxy.

## 3 Results

### 3.1 Properties of nearby SMBHs

The SMBH mass function provides a useful overview of the SMBH population at low redshift. To determine the average SMBH mass function with reduced sampling noise, we combine snipshot outputs as explained in § 2.3. Fig. 1 shows the SMBH mass function for the Ref-L100N1504 (dark blue line), AGNdT9-L050N0752 (light blue line) and Recal-L025N0752 (green line) simulations. In order to facilitate comparison with later plots and observational data, we include only the central black hole of each galaxy. The level of agreement between simulations is good (better than 0.2 dex) for SMBHs with mass . This level of agreement is encouraging, but note that the EAGLE simulations were calibrated to reproduce the normalisation of - relation at (see S15, Fig. 10). It is interesting to note the similarity of the high-mass end of the Ref-L100N1504 and AGNdT9-L050N0752 simulations. Given the more effective AGN feedback in the AGNdT9-L050N0752 model, we might have expected a divergence at the massive end due to its greater gas mass loss from galaxy groups. Clearly this is not the case. In Appendix A, we show that for SMBHs with mass the mass function depends strongly on the seed black hole mass and we indicate this region by dotted lines.

It is interesting to compare the simulation SMBH mass functions to the estimates based on observational data. Because SMBH masses can be directly determined only in an incomplete sample of galaxies (Kelly & Shen, 2013), it is important to note that the observational estimates of the SMBH mass function are indirect and must be inferred from the correlations between SMBHs and the properties of their host galaxy bulge. In the figure 1, we show estimate from Marconi et al. (2004) (grey circles), Shankar et al. (2004) (grey squares) and more recent data from Shankar et al. (2013) (red and grey regions). Shankar et al. (2013) and Shankar et al. (2004) use the - correlation, while Marconi et al. (2004) use the relation between SMBH mass and bulge luminosity. The simulated SMBH mass function is in reasonable agreement with the observational estimates from Shankar et al. (2004) and Marconi et al. (2004) over a wide mass range, but underestimates the abundance of the high-mass SMBHs when compared to Shankar et al. (2013). This discrepancy is somewhat surprising since both the simulation and Shankar et al. (2013) are calibrated to SMBH masses from McConnell & Ma (2013). The discrepancy, and the variance between observational estimates, illustrates the uncertainty in deriving the SMBH mass function from observational data. Indeed, *shankar2013 note that adopting different variations of the SMBH scaling relations leads to large variations in the inferred SMBH mass function, since the scatter and mass range covered by the data must be taken into account. For example, Shankar et al. (2013) show that the low-mass end of their mass function depends strongly on how the - correlation is applied to galaxies of different morphological types. The red region assumes that the relation can be applied regardless of morphology, while the grey region assumes that Sa and late-type galaxies do not have a SMBH. Thus, although there are some differences between the simulated SMBH mass function and the more recent observational estimates, these depend heavily on how the observational calibration data is extrapolated. For this reason, it is far better to validate the simulation by comparing to the black hole mass-stellar mass relation directly and in *schaye2015 we show that the simulation reproduces the observational data within their uncertainties.

Integrating the mass function, we obtain the predicted black hole mass density at . In the Ref-L100N1504 simulation, we find it to be , closely matching the observational value estimated by Yu & Tremaine (2002) (, adjusted to Planck Cosmology parameters), whose calculations are based on the velocity dispersion of early-types galaxies in the Sloan Digital Survey (SDSS). This value is also consistent with the result from Aird et al. (2010) who used hard X-ray luminosities to compute the total energy released by the SMBH population through time (), although is lower than the estimate of *marconi2004 ().

In Fig. 2 we dissect the SMBH mass function according to accretion regime. We distinguish inactive SMBHs (light blue), ADAFs (green) and SSDs (red). The three SMBH populations differ greatly in normalization. Most SMBHs are inactive, corresponding to of the total SMBH population with , while ADAFs (green line) correspond to and SSDs correspond to only of the population. The figure shows results from the Ref-L100N1504 simulation, but the breakdown is similar in the other simulations. Since SMBHs frequently switch between states between output times, we can view the differences in normalization as an average duty cycle, and interpret the relative normalisation of the mass functions as the probability of finding a SMBH in any given state.

The duty cycle, that is the fraction of the time a SMBH is active, shortens with increasing SMBH mass, with the probability of classifying a SMBH as an active SSD varying from 0.05 for to 0.01 . At higher masses, the probability of finding a present-day SMBH in the SSD state becomes extremely small. Restricting the comparison to the SMBH population with , the SSD fraction is on average. This is consistent with the observational estimates of the average AGN lifetime for that corresponds to years ( e.g. Marconi et al. 2004, Yu & Tremaine 2002).

These trends are not particularly sensitive to the choice of the threshold used to define SSD systems; however, an observational survey will only detect black holes that exceed an X-ray luminosity (or flux) limit. Purples lines in Fig 2 show the effect this has on the fraction of black holes that are detectable in an ideal hard (2–10 keV) X-ray survey. Our estimates account for bolometric corrections (as described in section 2.5) but do not account for additional selection effects, such as the difficulty of distinguishing faint AGN from emission associated with star formation. We focus on the results for an observational survey with a luminosity limit of . For the highest mass black holes, all the SSD systems are detected, but below a black hole mass of , the observable population become increasingly biased. For black holes of mass , the detected population accounts for only 6% of the SSD population and 0.3% of all black holes of that mass. As black hole masses drop below , the population becomes undetectable because of the Eddington limit. Fortunately, in practice, observational surveys are flux limited so that a range of luminosity limits can be probed; nevertheless, this exercise highlights the difficulty of constructing a complete census of the black hole population. We examine the X-ray emission of the simulation’s black hole population in more detail in section 3.3.1.

### 3.2 Properties of high-redshift SMBHs

Having examined the properties of SMBHs at low redshift, we now investigate the evolution of SMBHs. We will look at the evolution of the SMBH mass function, the - distribution and the SMBH mass-halo mass relation.

#### 3.2.1 The SMBH mass function

Fig. 3 investigates the evolution of the SMBH mass function. Between and , the normalization of the SMBH mass function rapidly evolves by an order of magnitude. While the overall normalisation changes little at lower redshifts (), the abundance of the most massive objects increases as the break in the mass function becomes shallower with cosmic time, and the dip seen at low masses at intermediate redshifts is filled in. The evolution of the SMBH mass function is similar to the evolution of the galaxy stellar mass function (see Furlong et al., 2015a, Fig. 2).

In Fig. 5, we show the evolution of the different accretion regimes. The contributions of ADAFs (green lines), SSDs (red lines) and inactive SMBHs (light blue) evolve relative to each other. At the abundance of the SMBH mass function is dominated by ADAFs and inactive SMBHs, preserving a similar shape to the total SMBH mass function. In contrast, the SSD population evolves rapidly in normalization from to , also showing a rapid decrease in characteristic mass. This is a feature of AGN ‘downsizing’ which we will return in the following sections. At the inactive SMBH population declines and the dominant populations are SSDs and ADAFs. The results described above show that there is switch in the nature of black hole accretion with cosmic time. Below , the population is dominated by inactive SMBHs or ADAFs and only a tiny fraction is undergoing strong accretion. At high redshift, SMBHs undergo much more frequent high Eddington-rate accretion events.

#### 3.2.2 The λEdd-MBH plane

In Fig. 5 we show Eddington ratio, , as a function of black hole mass in the Ref-L100N1504 simulation from redshifts 0 to 5. Overall, the median Eddington ratio decreases as a function of the black hole mass, with the SMBHs with mass the most active population in the simulation through cosmic time. However, there is large scatter ( dex) in the distribution of Eddington ratios for any given SMBH mass due to the high variability of the mass accretion rates.

For a given SMBH mass, the median value of moves towards lower values as redshift decreases. This trend is more evident in SMBHs with , where the median of log declines from at to at . For SMBHs with higher mass, the median changes less dramatically, consistent with the difference in the evolution of the active and inactive SMBH populations shown in Fig. 5. The figure also highlights that SMBHs of mass have an increasing tendency to be limited by the Eddington accretion rate at higher redshift, making it possible to build quasar-mass black holes early in the history of the Universe. In general, however, the SMBHs in the simulation accrete well below their Eddington limit.

#### 3.2.3 The MBH-M200 relation

Fig. 6 shows the evolution of the the central SMBH mass – halo mass relation for different redshifts from to .

We also include the estimates at from observations of strong gravitational lenses by Bandara et al. (2009) as grey diamonds and the grey line. In comparison to the simulations, Bandara et al. (2009) find a steeper relation, but this calculation was based on the assumed form of the - relation and not on direct observations of the mass of the central SMBHs.

The two simulations agree closely, and both show rapid BH growth in the halo mass range - demonstrating that this rapid growth phase does not depend on the details of the SMBH feedback scheme (or indeed the SMBH seed mass, as we show in Appendix A). The origin of this transition will be investigated in Bower et al.( 2016, in preparation), showing that it emerges as a result of a change in the hot gas content of the halo (see also Bower et al. 2006). In halo masses below the transition, the gas content of galaxies is regulated by stellar feedback; however, in more massive haloes, supernova-driven outflows stall as the hot halo becomes established and the gas content of the galaxy is regulated by the more energetic SMBH driven feedback. The median - relation evolves very little with redshift, so that at many massive SMBHs have already been assembled and a population of halos with already host SMBHs with . At higher redshifts the transition between the SMBH mass regimes becomes more abrupt, and SMBHs in this regime must grow rapidly as their halo mass increases which is consistent with the increasing median Eddington ratio seen in Fig. 5. SMBHs in higher-mass haloes () have released enough energy into the host halo to ensure that the cooling time becomes long, and the galaxy is starved of further fuel for star formation. This later process results in a self-regulated growth as noted by *booth_schaye2010 leading, together with the BH growth due to mergers, to the small scattering well-defined slope seen in Fig. 6.

### 3.3 Observable diagnostics of SMBH growth

In this section we investigate observables related to gas accretion onto SMBHs. We will focus on the hard and soft X-ray AGN luminosity functions and the evolution of the space density of AGN in hard X-rays through cosmic time.

#### 3.3.1 Hard X-ray luminosity functions

Fig. 7 shows the predicted hard (2–10 keV) X-ray luminosity function (HXRLF) over the redshift range to . Intrinsic X-ray luminosities have been derived using the bolometric corrections described in Sec. 2.5, and include only SMBHs in the SSD regime. We show two simulations, Ref-L100N1504 (dark blue lines) and AGNdT9-L050N0752 (light blue lines). The HXRLF shows strong evolution in shape and normalization in both simulations. Below the simulations agree with each other within dex in normalisation, with AGNdT9-L050N0752 slightly above the HXRLF of Ref-L100N1504. At the simulations are still similar but present higher discrepancies in the bright end of the HXRLF. The bright end of the HXRLF may be affected the limited statistics available in our volume or by variability on short timescales that is not resolved.

We also compare the predicted HXRLF to recent observational estimates based on the deep X-ray fields from Miyaji et al. (2015) (green circles), Aird et al. (2015) (red pentagons) and Buchner et al. (2015) (bright blue hexagons). Overall, the comparison between the observations and simulations is encouraging, given the simplicity of the subgrid model used. We stress that the parameters of SMBH growth and AGN feedback were calibrated to match the stellar mass function at and the normalization of the - relation, not to match the evolution of AGN. Looking in detail, the observations are in good agreement (within the observational error bars) out to ; however, at and higher luminosities (), the simulation HXRLF appears to decline with more quickly than seen in the observations. The discrepancy does not appear to be due to the sampling statistics but we cannot entirely rule out the possibility that it is due to the finite volume of the simulation (see Appendix A) because we sample a relatively small number of massive galaxies in the simulations. Above , there is an overabundance of low luminosity () AGN in EAGLE. In this regime, however, the observational constraints are quite uncertain, as can be seen by comparing different observational datasets.

We have mentioned previously that one possible factor that affects the HXRLF is the variability of AGN that is not resolved in the simulation. Short time-scale variations could originate from fluctuations in accretion disk viscosity, for example. In galactic binary systems such as stellar remnant compact objects, order of magnitude variations in the accretion rate arise from the ionisation instability in accretion disks, and similar instabilities may be present in SMBHs (Done et al., 2007). Such ‘flickering’ cannot be resolved in our simulations which only attempt to model variations in the gas supply rate on pc scales.

An illustration of the effect of short-timescale variability is shown in Fig. 8. Here, we show the effect of convolving source luminosities with a log-normal distribution with between 0.3 and 0.5 dex per luminosity.

We have chosen the log-normal distribution as a simple way to illustrate the potential impact of flickering since we are interested in the effect of order of magnitude variations in source luminosity. It is important to note that the central value of the convolution kernel has been set so that the average (expectation) luminosity is independent of . We have explored relatively high values in order to assess the maximum impact of unresolved variability in the simulation. A source with is, instantaneously, an order of magnitude brighter or fainter than the mean luminosity for 5% of the time, and a factor of 3 brighter or fainter for 32% of the time. Values higher than 0.5 would imply that the instantaneous is almost unrelated to the SMBH accretion rate.

Solid lines reproduce the Ref-L100N1504 simulation and observational data from Fig. 7. The effect of including unresolved variability is shown as red dot-dashed and red dashed lines in Fig. 8. As the width of the convolving Gaussian is increased intrinsically low luminosity sources scatter in the high luminosity bins, and the simulation tends to agree better with the observational data. The overall effect is relatively small, however, and does not seem able to reconcile the simulation with the observational data. The convolution has little effect on the fainter luminosities () and so cannot account for the overabundance of faint sources in the simulation at . We have already stressed that the observational measurements are uncertain in this regime.

Although the volume of the simulation is rather small for the characterisation of extreme high redshift events, for completeness we show the predicted hard X-ray luminosity function in EAGLE at in Fig. 9. We see that the HXRLF amplitude decreases with redshift and evolves rapidly in shape. Above , the simulation suffers from particularly poor sampling for AGNs with (indicated by the dashed line). Comparing to observational estimates from Buchner et al. (2015), we find reasonable agreement between and over the luminosity range that we are able to probe. Observational data is not available at higher redshifts in this luminosity range, and the figure presents the model predictions.

#### 3.3.2 Soft X-ray luminosity function

Soft (0.5–2 keV) X-ray measurements provide a useful complement to hard X-ray surveys. The evolution of the soft X-ray AGN luminosity function (SXRLF) has been investigated by e.g. Miyaji et al. (2000); Hasinger et al. (2005); Aird et al. (2015). A complication, however, is that an uncertain fraction of sources will be obscured by the gas column along the line of sight through the host galaxy to the AGN. Rather than trying to model the effects of obscuration in the simulation, we compare to the results of *aird2015 in which the effects of obscuration have been empirically corrected.

Fig. 10 shows the SXRLF for the simulation Ref-L100N1504 (solid lines). To illustrate the importance of obscuration, we have included blue dot-dashed lines to show the expected abundance of detectable (i.e. unobscured) objects using the prescription of Hasinger (2008). Comparing the SXRLF with obscuration to the one without, the fraction of obscured AGN can vary between and with the largest values found at low luminosities () and the smallest at high luminosities (). Note, however, that these ratios are not a theoretical prediction of the simulations, but the effect of corrections derived from

We compare the predicted SXRLF (solid lines) to the observational estimates from *aird2015 (red pentagons). The observed counts are corrected for the effects of obscuration by comparing the hard- and soft-band X-ray data. The data points should therefore be compared to the solid lines derived from the simulations. The simulation broadly reproduces the observed evolution across cosmic time, particularly for the faintest part of the SXRLF () at . Even at low redshift, the brightest part of the SXRLF is steeper than observed, but the discrepancy is only greater than the observational uncertainties in the region where we have fewer than 10 sources per bin (shown as dashed lines).

As we discussed above, we would expect some additional contribution from unresolved variability, and we show the effect of log-normal luminosity variations between width of 0.3 and 0.5 dex as red dot-dashed and red dashed lines respectively. As we found in Fig. 8, this has relatively little impact. Towards high redshifts (), the simulations predict a higher amplitude of the SXRLF (particularly at luminosities of ). This might indicate an overabundance of the faint AGN in the simulations, but may also be due to a greater redshift dependence of obscuration than accounted for by Aird et al. (2015).

However, a similar over-abundance is seen in both soft and hard X-ray luminosity functions, and in general in both X-ray bands the luminosity functions follow a similar evolution. It seems therefore that the offset between the simulation and the observations must either be real (in the sense that the numerical implementation of SMBH accretion used in the simulations generates an excess of low luminosity sources) or be due to observational selection effects (for example, we have not attempted to model observational selection effects such as the difficulty of detecting faint AGN against the galaxy’s nuclear star formation).

In general terms, the evolution of the SXRLF in the simulation evolves in broad agreement with the observational constraints, opening a window to explore more deeply the connection between BH accretion rates, obscuration and the gas and star formation properties of galaxies.

#### 3.3.3 Evolution of the comoving number density of AGN

Fig. 11 shows the evolution of the comoving number density of AGN in the simulation Ref-L100N1504 (solid lines). As we discussed above, some AGN variability may not be accounted for in the simulation and we illustrate the effect of convolving source luminosities with a log-normal distribution of width 0.3 to 0.5 dex. We split the AGN population into three luminosity bands as indicated in the legend. This figure reproduces the information already seen in Fig. 8, but allows us to investigate the evidence for AGN ‘downsizing’ more clearly. In the observational data, higher luminosity sources peak in abundance at progressively higher redshifts. A similar trend is seen in the simulations: brighter AGN (; green lines) peak at redshifts greater than 3 while fainter AGN (; dark blue lines) peak at .

For comparison, we show recent estimates from *ueda2014 and *aird2015. We also show values obtained by integrating the luminosity functions from *buchner2015. There is reasonable agreement between simulations and observations for the lower luminosity bins, (dark blue line) and (light blue line), out to . Moving towards higher redshifts, the simulations predict too many faint AGN in comparison to observations, the discrepancy becoming almost 0.8 dex in the comoving number density at . We note, however, similar discrepancies are seen in observational data due to uncertainties in detecting faint AGN at high redshifts. For example, the comoving number density of AGN of from Ueda et al. (2014) (blue circles, ) is higher by 0.8 dex compared to Aird et al. (2015) (blue pentagons, ) at . Moreover, *giallongo2015 recently reported a more abundant population of faint AGNs at - by studying AGN candidates from the multiwavelength CANDELS deep surveys (Grogin et al., 2011; Koekemoer et al., 2011). The sample was initially selected by the UV rest frame of the parent galaxy and thus is able to account for sources with marginal X-ray nuclear detections. This population might extend to even higher redshift (). *madau2015, based on this result, suggest that such contribution of active galaxies drives the reionization of hydrogen and helium satisfying several observational constrains such as the observed flatness of HI photoionization rate between and and the estimated integrated Thompson scattering optical depth () found in the Lymann opacity of the intergalactic medium and cosmic microwave (CMB) polarization.

The comoving number density of the brightest AGN is low in the simulations compared to the observational estimates. However, the comoving number density of the brightest AGN can be affected by additional AGN variability combined with the low numbers of bright AGN in the finite simulation volume. As we show in the figure, convolution with log-normal flickering of 0.5 dex goes some way to account for the high abundance of bright AGN seen in the observations.

Overall, while the ‘downsizing’ trend is present in the simulations, it is not as clear as suggested by the observational data. In particular, the abundance of AGN of a particular luminosity has a broader plateau than suggested by observations, principally because the rapid decline observed in the abundance of faint AGN at high redshift seen in observations is shallower in the simulations. At almost all redshifts, however, the simulations tend to underpredict the abundance of the brightest AGN, and their peak of their abundance occurs at too high redshift. We should be careful not to over-interpret this apparent discrepancy, however, since the abundance of these objects is poorly sampled in the simulation volume.

## 4 Conclusions

We have examined the evolution of supermassive black holes (SMBHs) across cosmic time predicted by the EAGLE simulations (\al@schaye2015,crain2015; \al@schaye2015,crain2015). The EAGLE project consists of a suite of hydrodynamical simulations with state-of-the-art subgrid models of galaxy formation including radiative cooling, star formation, reionization, abundance evolution, stellar evolution and mass loss, feedback from star formation, and SMBH growth and AGN feedback. The parameters of these subgrid models were calibrated to reproduce the observed galaxy mass function and sizes at . In particular, the efficiency of AGN accretion and feedback were set to reproduce the break in the stellar mass function at and the normalization of the SMBH mass-stellar mass relation at . It is important to emphasise that the subgrid models of SMBH growth and AGN feedback do not make any explicit distinction between quasar and radio modes, and that we only distinguish sources with high and low Eddington ratios during the analysis. The main findings are summarised as follows:

• The main properties of nearby SMBHs are reproduced well in the EAGLE simulations. Within the observational uncertainties, the SMBH mass function is similar to estimates from *shankar2004, *marconi2004 and *shankar2013, and the density of SMBHs in the local Universe is also comparable to that observed. This agreement is partly the result of the calibration strategy (see Fig. 1). As a post-processing step, we divide the present-day black hole mass function as a function of Eddington ratio, . We associate sources with with X-ray luminous Shakura-Sunyaev discs (SSDs), sources with as Advection Dominated Accretion Flows (ADAFs) and classify that sources with as inactive. At low redshift the mass function is dominated by inactive and ADAF black holes (Fig. 2). Assuming that SMBHs cycle between the SSD, ADAF and inactive states, we estimate that the duty cycle for SMBHs in the SSD state is , which is comparable to the observational estimates (e.g. Soltan 1982, Marconi et al. 2004, Yu & Tremaine 2002).

• The mass function of SMBHs in the EAGLE simulation evolves rapidly in amplitude from to . At redshift a large fraction of the population has already been formed (Fig. 3). Between and the present-day the mass function evolves more gradually in normalisation. When we break this evolution down by accretion state, we find that luminous SSD systems, while a minor contribution at the present-day, become increasingly dominant at high redshift (Fig. 5). This trend can also be clearly seen by examining the evolution of the Eddington ratio distribution directly (Fig. 5).

• We examine the dependence of black hole mass on the dark matter halo mass, (Fig. 6). The relation has a characteristic shape, with SMBH masses growing little above the seed mass in haloes less massive than , but showing a sharp rise in more massive haloes. The fast growth of the SMBH ends when its mass exceeds . SMBHs follow an almost linear trend with in larger haloes. The characteristic shape of this relation evolves little with redshift, with a suggestion that the steep rise in mass becomes more abrupt as redshift increases.

• The black hole mass function, the Eddington ratio distribution and the SMBH dependence on halo mass cannot be directly observed and must be inferred by combining observational surveys with, for example, a calibration between black hole mass and stellar mass. To facilitate a more direct comparison between the model and observational data, we compute the X-ray luminosity function in the rest frame for SMBHs in the SSD state. We use bolometric corrections from *marconi2004 to convert the bolometric output predicted by the model into the intrinsic hard and soft X-ray luminosities. We compare the hard-band X-ray luminosity functions with the observational measurements from *miyaji2015, *aird2015 and *buchner2015 (Fig. 7). The finite volume of the simulation limits the comparison to hard X-ray AGN luminosities lower than . At low redshifts, the simulations agree extremely well with the observational data. At higher redshift () the simulations tend to underpredict the abundance of high luminosity sources with , although we cannot rule out the possibility that this is a result of the simulation’s limited volume. At and above, the amplitude of the predicted luminosity function appears higher than observed, particularly around , although some caution is required since we have not attempted to include observational selection effects (other than bolometric corrections) in our predictions. We find a similar result when we compare to obscuration-corrected soft-band X-ray luminosity functions from *aird2015.

• The hard X-ray luminosity functions we derive include the effect of variability captured by the model due to gas flows on kpc scales, but unresolved variability (for example due to flickering in the accretion disk) that may cause additional fluctuations in luminosity. Given the steepness of the X-ray luminosity function, this could have an important impact. However, we present a simple model based on additional log-normal distributed flickering to show that this has only a limited impact on the comparison with observational data (Fig. 8).

• We investigate AGN downsizing in the simulation (Fig. 11). The observed trend seen in observational data is qualitatively reproduced: the comoving number density of higher-luminosity AGN peaks at higher redshift, and the simulations are in good quantitative agreement with the observational data for at . At higher redshifts, the simulations produce more active SMBHs than observed, resulting in a shallower roll-over of the AGN abundance. The finite volume of the simulations and the possible effects of flickering make difficult to reliably compare the abundance of the more luminous AGN. Taken at face value, the simulations do not predict the rapid rise in the abundance of the brightest () objects seen in some observational surveys between and . Larger volume simulations, and a better understood model for AGN flickering, are required to determine if this is due to a real discrepancy between the hydrodynamical model and the observational data.

The results we find are broadly consistent with other simulations and semi-analytic calculations. For example, *fanidakis12 used a version of the semi-analytic GALFORM code, similar to that of Bower et al. (2006), in which galaxy formation is approximated as a network of analytic differential equations that are applied to halos that grow in a dark matter N-body simulation. It is assumed that SMBHs grow either by accretion from the diffuse gas halo, if this is stable against cooling, or as a result of gas flows produced during merger and disk instability driven starbursts. The semi-analytic model is able to probe large volumes and hence more luminous sources, and the model indeed generates a population of very luminous sources at high redshift, improving the match to the observational data. Similarly to EAGLE, theses calculations do not show the strong ‘down-sizing’ trend inferred from the observational data and the authors conclude that the perceived evolution is largely the result of obscuration.

Recently, *hirschmann2014 and *sijacki2015 have presented an analysis of black hole properties in large volume cosmological simulations. *hirschmann2014 combined simulations of a Mpc region at low resolution (a factor of higher than the particle mass of the Ref-L100N1504 simulation) run to , with the results obtained from a 68 Mpc region with a resolution similar to that of the EAGLE simulation, but run only to . Their prescription of AGN feedback is extensively based on Springel et al. (2005a) and two explicit modes of AGN feedback are assumed, namely, a quasar and a radio mode with a switchover depending on the source Eddington ratio. The efficiency of feedback in radio mode is times larger than that in quasar mode. These models generally fit the observed AGN luminosity functions reasonably well, although they also find that the abundance of sources tends to be overestimated when the resolution of their large volume calculation is increased. Accounting for the differences in resolution and volume, their results appear compatible with our own.

*sijacki2015 presented an analysis of SMBHs in the Illustris simulation (Vogelsberger et al., 2014) which, while similar in volume and resolution of our work, differs greatly in its implementation of AGN feedback and accretion onto SMBHs. In particular, the Illustris simulation employs different schemes for feedback in high and low Eddington ratio sources. In low accretion states radio feedback is implemented by depositing energy in thermal ‘bubbles’ at some distance from the central galaxy (Sijacki et al., 2007). In this high accretion mode feedback energy is dumped at the location of the BH at every timestep, a procedure that is expected to result in significant radiative losses at this resolution (e.g. Dalla Vecchia & Schaye 2012) The black hole mass function derived from the Illustris simulations differs significantly in shape from that in EAGLE, being essentially a power-law that increases in steepness as a function of mass. At low redshift, rare SMBHs (more massive than ) are more abundant, but most of these SMBHs accrete at rates less than of the Eddington limit. The results for the luminosity functions are broadly similar to ours (within the uncertainty of the data that are shown). In particular, they also find that the model tends to overpredict the abundance of moderate luminosity () AGN at . In terms of AGN downsizing, the model does appear to capture the rapid decline of the abundance of high-luminosity sources, although it is unclear whether this is largely affected by the selection of sources based on an Eddington ratio criterion of , given the significant difference in the black hole mass functions of the simulations.

Although most of our qualitative results seem compatible with these earlier works. It is important to stress the greater simplicity of the AGN feedback model used in EAGLE, and the fact that basic galaxy properties like stellar masses and sizes are better reproduced. It is therefore encouraging that the model, which uses a single mode of AGN feedback and in which AGN feedback energy is a fixed fraction of the accretion rate, captures so many of the trends seen in observational data. The results we have presented from the EAGLE simulations open a new window to investigate the co-evolution of the SMBH growth and galaxy evolution. In future work, we will investigate more consistently the obscuration of AGN due to gas and dust by including the properties of the host galaxy. We will also investigate the effects of AGN feedback on the host galaxies and how this evolves through cosmic time.

## Acknowledgements

Authors thank Johannes Buchner, James Aird, Francesco Shankar and Takamitsu Miyaji for providing their observational data and the referee for the useful comments. This work would have not be possible without Lydia Heck’s technical support. YRG thanks Johannes Buchner for his useful comments about observations. YRG gratefully acknowledges financial support from the Mexican Council for Science and Technology (CONACYT) (Studentship No. 213183) and partial support from Center of Excellence in Astrophysics and Associated Technologies (PFB 06). This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/L00075X/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. We also gratefully acknowledge PRACE for awarding us access to the resource Curie based in France atTres Grand Centre de Calcul. This work was sponsored by the Dutch National Computing Facilities Foundation (NCF) for the use of supercomputer facilities, with Financial support from the Netherlands Organization for Scientific Research (NWO). The research was supported in part by the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) / ERC Grant agreements 278594-GasAroundGalaxies, GA 267291 Cosmiway, and 321334 dustygal, the Interuniversity Attraction Poles Programme initiated by the Belgian Science Policy OWNce ([AP P7/08 CHARM]), the National Science Foundation under Grant No. NSF PHY11-25915, the UK Science and Technology Facilities Council (grant numbers ST/F001166/1 and ST/I000976/1. We thank contributors to SciPy , Matplotlib , and the Python programming language .

## References

• Abramowicz et al. (1995) Abramowicz, M. A. , Chen, X. , Taam, R. E., 1995, ApJ, 452, 379.
• Aird et al. (2010) Aird, J. , Nandra, K. , Laird, E. S. and Georgakakis, A. , Ashby, M. L. N. , Barmby, P. , Coil, A. L. , Huang, J.-S. , Koekemoer, A. M. , Steidel, C. C. and Willmer, C. N. A. 2010, MNRAS, 401, 2531.
• Aird et al. (2015) Aird, J., Coil, A. L., Georgakakis, A., et al. 2015, arXiv:1503.01120
• Angles-Alcazar et al. (2016) Anglés-Alcázar, D., Davé, R., Faucher-Giguère, C.-A., Özel, F., & Hopkins, P. F,, 2016, arXiv:1603.08007
• Bahé et al. (2016) Bahé, Y. M., Crain, R. A., Kauffmann, G., et al. 2016, MNRAS, 456, 1115
• Bandara et al. (2009) Bandara, K., Crampton, D., Simard,  L., 2009, ApJ, 704, 1135.
• Binney & Tabor (1995) Binney J., Tabor G., 1995, MNRAS, 276, 663.
• Bondi & Hoyle (1944) Bondi  H., Hoyle  F., 1944, MNRAS, 104, 273.
• Booth & Schaye (2009) Booth C. M., Schaye  J., 2009, MNRAS, 398, 53.
• Booth & Schaye (2010) Booth C. M., Schaye  J., 2010, MNRAS, 405, L1.
• Bower et al. (2006) Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645.
• Buchner et al. (2015) Buchner J., Georgakakis A., Nandra K., et al. 2015, ApJ, 802, 89.
• Chabrier (2003) Chabrier, G., 2003, PASP, 115, 763.
• Choi et al. (2013) Choi, E., Naab, T., Ostriker, J. P., Johansson, P. H., & Moster, B. P., 2014, MNRAS, 442, 440.
• Churazov et al. (2001) Churazov E., Brüggen M., Kaiser C. R., Böhringer H., Forman W., 2001, ApJ, 554, 261.
• Crain et al. (2009) Crain R. A., Theuns  T., Dalla Vecchia  C., Eke V. R. , Frenk S. C., Jenkins  A., Kay S. T. and Peacock J. A. et al., 2009, MNRAS, 399, 1773.
• Crain et al. (2015) Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, arXiv:1501.01311
• Croton et al. (2006) Croton D.  J., Springel V., White S.  D.  M., De Lucia  G., Frenk S. C., Gao L., Jenkins  A. and Kauffmann  G. et al., 2006, MNRAS, 365, 11.
• Cullen & Dehnen (2010) Cullen, L. and Dehnen, W., 2010, MNRAS, 408, 669.
• Dalla Vecchia & Schaye (2012) Dalla Vecchia  C. and Schaye  J., 2012, MNRAS, 426, 140.
• De Lucia & Blaizot (2007) De Lucia G., Blaizot J., 2007, MNRAS, 375, 2
• Di Matteo et al. (2008) Di Matteo, T. and Colberg, J. and Springel, V. and Hernquist, L. and Sijacki, D., 2008, ApJ, 676, 33.
• Done et al. (2007) Done, C., Gierliński, M., & Kubota, A. 2007, A&ARv, 15, 1
• Durier & Dalla Vecchia (2012) Durier, F. and Vecchia D. C., 2012, MNRAS, 419, 465.
• Fanidakis et al. (2012) Fanidakis N. and Baugh C. M. and Benson A. J. and Bower R. G. and Cole S. and Done C. and Frenk C. S. and Hickox R. C. and Lacey C. and Del P. Lagos C., 2012, MNRAS, 419, 2797.
• Ferland et al. (2013) G. J. Ferland, R. L. Porter, P. A. M. van Hoof, R. J. R. Williams, N. P. Abel, M. L. Lykins, Gargi Shaw, W. J. Henney, and P. C. Stancil, 2013, Rev. Mex. Soc., 49, 1.
• Furlong et al. (2015a) Furlong, M., Bower, R. G., Theuns, T., et al. 2015, MNRAS, 450, 4486
• Furlong et al. (2015b) Furlong, M., Bower, R. G., Crain, R. A., et al. 2015, arXiv:1510.05645
• Giallongo et al. (2015) Giallongo, E., Grazian, A., Fiore, F., et al., 2015, AAP, 578, A83
• Graham (2016) Graham, A. W., 2016, Galactic Bulges, 418, 263
• Grogin et al. (2011) Grogin, N. A., Kocevski, D., Faber, S. et al. 2011, ApJS, 197, 37
• Haardt & Madau. (2001) Haardt, F. and Madau, P., 2001, Clusters of Galaxies and the High Redshift Universe Observed in X-rays, astro-ph/0106018.
• Hasinger et al. (2005) Hasinger, G. and Miyaji, T. and Schmidt, M., 2005, A&A, 441, 417.
• Hasinger (2008) Hasinger, G., 2008, A& A, 490, 905.
• Hirschmann et al. (2014) Hirschmann, M., Dolag, K., Saro, A., et al. 2014, MNRAS, 442, 2304
• Hirschmann et al. (2012) Hirschmann M., Somerville R. S., Naab T., Burkert A., 2012, MNRAS, 426, 237
• Hopkins et al. (2007) Hopkins, P. F., Richards, G. T., & Hernquist, L., 2007, ApJ, 654, 731
• Hopkins et al. (2009) Hopkins P. F., and Hernquist L., Cox T. J. , Keres D. and Wuyts S., 2009, ApJ, 691, 1424.
• Jenkins (2010) Jenkins  A., 2010, MNRAS, 403, 1859.
• Jenkins (2013) Jenkins  A., 2013, MNRAS, 434, 2094.
• Kelly & Shen (2013) Kelly B. C., Shen Y., 2013, ApJ, 764, 45.
• Kennicutt (1998) Kennicutt, Jr., R. C., 1998, ARA&A, 36, 189.
• Khandai et al. (2015) Khandai, N. and Di Matteo, T. and Croft, R. and Wilkins, S. M. and Feng, Y. and Tucker, E. and DeGraf, C. and Liu, M.-S., 2015, MNRAS, 450, 1349
• Koekemoer et al. (2011) Koekemoer, A. M., Faber, S., Ferguson, H. et al. 2011, ApJS, 197, 36
• Kormendy & Ho (2013) Kormendy, J. and Ho, L. C., 2013, ARA&A ,51,511.
• Lagos et al. (2015) Lagos, C. d. P., Crain, R. A., Schaye, J., et al., 2015, arXiv:1503.04807
• Lansbury et al. (2015) Lansbury, G. B., Gandhi, P., Alexander, D. M., et al., 2015, ApJ, 809, 115
• Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
• Lusso et al. (2012) Lusso, E., Comastri, A., Simmons, B. D., et al., 2012, MNRAS, 425, 623
• Madau & Haardt. (2015) Madau, P., & Haardt, F., 2015, ApJ, 813, L8
• Marconi et al. (2004) Marconi, A. and Risaliti, G. and Gilli, R. and Hunt, L. K. and Maiolino, R. and Salvati, M., 2004, MNRAS, 351, 169.
• McAlpine et al. (2015) McAlpine S., Helly J.C., Schaller M., et al. 2015, arXiv:1510.01320
• McCarthy et al. (2010) McCarthy I. G., Schaye  J., Ponman T. J., Bower R. G., Booth C. M., Dalla Vecchia C., Crain R. A. and Springel V., 2010, MNRAS, 406, 822.
• McConnell & Ma (2013) McConnell, N. J., & Ma, C.-P., 2013, ApJ, 764, 184
• Miyaji et al. (2000) Miyaji T., Hasinger G., Schmidt M., 2000, A&A, 353, 25.
• Miyaji et al. (2015) Miyaji, T., Hasinger, G., Salvato, M., et al., 2015, ApJ, 804, 104
• Narayan & Yi (1994) Narayan, R. and Yi, I., 1994, ApJ, 428, L13.
• Planck collaboration et al. (2013) Planck Collaboration and Ade, P. A. R. and Aghanim, N. and Alves, M. I. R. and Armitage-Caplan, C. and Arnaud, M. and Ashdown, M. and Atrio-Barandela, F. and Aumont, J. and Aussel, H. and et al., 2013, arxiv e-prints, 1303.5062.
• Portinari et al. (1998) Portinari L., Chiosi C., & Bressan A. 1998, A& A, 334, 505.
• Power et al. (2010) Power, C. and Nayakshin, S. and King, A., 2011, MNRAS, 412, 269.
• Rees et al. (1982) Rees, M. J. and Begelman, M. C. and Blandford, R. D. and Phinney, E. S., 1982, NATURE, 295, 17.
• Rosas-Guevara et al. (2015) Rosas-Guevara Y. M., Bower R.G., Schaye J., et al. 2015, MNRAS, 454, 1038
• Shakura & Syunyaev (1973) Shakura N.I ., Syunyaev R. A., 1973, A&A , 24, 337.
• Schaller et al. (2015) Schaller, M., Dalla Vecchia, C., Schaye, J., et al., 2015, MNRAS, 454, 2277
• Schaye (2004) Schaye, J., 2004, ApJ, 609, 667
• Schaye & Dalla Vecchia (2008) Schaye J., Dalla Vechia C., 2008, MNRAS, 383, 1210.
• Schaye et al. (2010) Schaye  J., Dalla Vecchia  C., Booth C. M., Wiersma R. P. C., Theuns  T., Haas M. R., Bertone  S. and Duffy A. R.et al., 2010, MNRAS, 402, 1536.
• Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521
• Schmidt (1959) Schmidt, M., 1959, ApJ, 129, 243.
• Shankar et al. (2004) Shankar, F. and Salucci, P. and Granato, G. L. and De Zotti, G. and Danese, L., 2004, MNRAS, 354, 1020.
• Shankar et al. (2013) Shankar, F., 2013, Classical and Quantum Gravity, 30, 244001
• Sijacki et al. (2007) Sijacki  D., Springel  V., Di Matteo  T. and Hernquist  L., 2007 MNRAS, 380, 877.
• Sijacki et al. (2015) Sijacki, D., Vogelsberger, M., Genel, S., et al. 2015, MNRAS, 452, 575
• Soltan (1982) Soltan  A., 1982, MNRAS, 200, 115.
• Springel et al. (2005a) Springel  V., Di Matteo  T. and Hernquist  L., 2005, MNRAS, 361, 776.
• Springel (2005b) Springel,  V., 2005b, MNRAS, 364, 1105.
• Steffen et al. (2003) Steffen, A. T. and Barger, A. J. and Cowie, L. L. and Mushotzky, R. F. and Yang, Y., 2003, ApJ, 596, L23.
• Trayford et al. (2016) Trayford, J. W., Theuns, T., Bower, R. G., et al. 2016, arXiv:1601.07907
• Tremaine et al. (2002) Tremaine  S., Gebhardt  K., Bender  R., Bower  G., Dressler  A., Faber S.  M., Filippenko A.  V. and Green  R. et al., 2002, ApJ, 574, 740.
• Treister et al. (2009) Treister, E., Urry, C. M., & Virani, S., 2009, ApJ, 696, 110
• Ueda et al. (2003) Ueda, Y. and Akiyama, M. and Ohta, K. and Miyaji, T., 2003,ApJ , 598, 886.
• Ueda et al. (2014) Ueda, Y., Akiyama, M., Hasinger, G., Miyaji, T., & Watson, M. G. 2014, ApJ, 786, 104
• Vasudevan & Fabian (2009) Vasudevan R. V., Fabian A. C., 2009, MNRAS, 392, 1124
• Vogelsberger et al. (2014) Vogelsberger, M., Genel, S., Springel, V., et al., 2014, MNRAS, 444, 1518
• Yu & Tremaine (2002) Yu, Q. and Tremaine, S., 2002, MNRAS, 335, 965.
• Wiersma et al. (2009a) Wiersma R. P. C., Schaye  J., Theuns  T., Smith B. D, 2009, MNRAS, 393, 99.
• Wiersma et al. (2009b) Wiersma R. P. C., Schaye  J., Theuns  T., Dalla Vecchia  C. and Tornatore  L., 2009, MNRAS, 399, 574.

## Appendix A Convergence Tests

Following the discussion of the numerical convergence in S15, we use this appendix to investigate the impact of the simulation volume on AGN observables. We also present the effects of varying the initial seed SMBH mass and simulation resolution on the AGN luminosity functions. The simulations that we consider are described in Table 3 and Table 1.

Fig. 12 investigates the sensitivity to the volume of the simulation of the hard X-ray luminosity function (HXRLF) at , the median of the SMBH mass-halo mass relation and the SMBH mass function at . Dash lines show the limit of our sampling statistics. The left panel of Fig. 12 shows a good convergence in the HXRLF. The discrepancies for the HXRLF are smaller than 0.2 dex between simulations in the brightest part of the HXRLF. The plots in the figure also show excellent agreement within these limitations.

Fig. 13 shows the same panels, but compares the simulations Ref-L100N1504 (blue), Small-seeds-L050N0752 (pink) and Recal-L025N0752 (green). Although we compare different box sizes, the previous figures show that this is not a concern providing that the sampling statistics are appropriately accounted for. The left panel shows that the HXRLF is insensitive to resolution and to a change in the model for AGN feedback. Decreasing the SMBH seed mass by an order of magnitude has only a small effect (compared to observational uncertainties) at the faintest luminosities shown. The middle panel of Fig. 13 shows the median of the distribution of the - relation at . The median and scatter of each simulation shows a similar shape, however, the simulation Small-seeds-L050N0752 presents a sharper rise in halos of in comparison to the other simulations. This results from the SMBHs in Small-seeds-L050N0752 having to grow faster to reach the self-regulated - relation. Note however, that this steep rise persists in the three simulations, indicating that this fast growth of SMBHs in such haloes are independent of the subgrid parameters of the simulation. The right panel presents the SMBH mass function at in simulations Ref-L100N1504, Small-seeds-L050N0752 and Recal-L025N0752. The agreement between Ref-L100N1504 and Recal-L025N0752 is better than 0.2 dex, comparable to the differences in the galaxy mass functions of these simulations. In contrast, the lower mass end of the SMBH mass function ( ) is strongly affected by the SMBH seed mass. The simulation Small-seeds-L050N0752 predicts lower values for the SMBH mass function by dex for SMBHs with mass smaller than . Nevertheless, the massive end of the SMBH mass function present an impressive level of agreement between simulations.

## Appendix B Choice of Accretion Regimes

In section 2.4, we define two active accretion regimes that depend on the value of Eddington ratio. We assume that black holes with that are higher than are luminous sources of X-rays (since the nuclear disc is thin and radiative cooling efficient) and consider them as SSDs. Lower luminosity sources are assumed not to contribute to the X-ray luminosity functions we show in the main text. In this appendix, we explore the impact of varying this limit. Fig. 14 shows the HXRLF considering SMBH to be X-ray luminous when is larger than a minimum that varies from to . We show redshifts from to . The dependence is weak, especially for the bright end of HXRLF. For the faint end of the HXRLF, (AGN with lower than ) the difference becomes of dex, comparing and . This discrepancy becomes smaller with increasing redshift and is smaller than the observational error bars.

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

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