© ESO, 2013

1. Introduction

That every galaxy of significant size harbors a supermassive black hole (SMBH) (with masses ~10⁶ − 10⁹ M_⊙) is probably one of the most remarkable discoveries of modern astrophysics. The growth of SMBHs, as manifested by the active galactic nuclei (AGN), has been observed over a broad range of electromagnetic energies, from radio to hard X-rays and gamma-rays. The observed AGN X-ray emission, believed to stem from the upscattering of softer accretion disk photons via inverse Compton mechanism by the hot electron corona, has proven to be the most effective way of selecting large samples of AGN over large cosmological volumes¹. This efficiency comes from the fact that to absorb X-rays one needs significant column densities of absorbing material, which is particularly true in the hard X-ray band, where N_H ~ 10²⁴ cm^-2 is required.

With XMM-Newton² and Chandra³ X-ray observatories, two of the most advanced X-ray instruments in existence, more than 20 deep extragalactic X-ray surveys have been performed over varying sky areas and limiting sensitivities (Brandt & Hasinger 2005). The most noticeable among these are ~2 Ms Chandra Deep Field North (Alexander et al. 2003) over ≃ 448 arcmin² sky area, ~4 Ms Chandra Deep Field South (Xue et al. 2011) over ≃ 465 arcmin², and ~3 Ms XMM-Newton deep survey covering Chandra Deep Field South (Comastri et al. 2011).

Despite the small sky areas covered by those deep surveys, a relatively large number of detected AGN has allowed reasonably good determination of the AGN luminosity function (LF) and its evolution over cosmic time (Ueda et al. 2003; Hasinger et al. 2005; Aird et al. 2010). Even though these narrow survey geometries are not very suitable for measuring the spatial clustering properties of AGN, somewhat wider and shallower surveys like X-Boötes, XMM-LSS, AEGIS, and XMM-COSMOS have allowed measuring the two-point clustering statistics (Murray et al. 2005; Gandhi et al. 2006; Coil et al. 2009; Gilli et al. 2009; Allevato et al. 2011), albeit with relatively large uncertainty on large scales. The large-scale clustering measurements are expected to improve significantly with the upcoming Spectrum-X-Gamma/eROSITA⁴⁵ space mission (Predehl et al. 2010), which is planned to cover entire sky down to the limiting sensitivity of ~10^-14 erg/s/cm² in the 0.5 − 2 keV band (Merloni et al. 2012; Kolodzig et al. 2013b). A good summary of the current status and future prospects for the X-ray-selected AGN clustering can be found in Cappelluti et al. (2012); Krumpe et al. (2013).

Even though the evolution of AGN over cosmic time is an interesting subject on its own, the remarkable discoveries of correlations between SMBH masses and host galaxy properties (see, e.g., Kormendy & Richstone 1995; Ferrarese & Ford 2005, for reviews), hinting at the underlying AGN-galaxy co-evolution mechanisms, have significantly broadened the importance and actuality of those topics in the astronomical community. Although controversial in terms of observational data (e.g., Ferrarese 2002; Kormendy & Bender 2011), in currently favored cold dark matter (CDM) cosmologies, where the structure grows in a hierarchical fashion starting from the smallest scales (i.e., the bottom-up scenario), one may expect that SMBH masses M_BH should also correlate with the masses of the host dark matter halos M_h. Indeed, this is what one commonly finds in semi-analytic AGN-galaxy co-evolution models, e.g., Fanidakis et al. (2012). There is a significant scatter in the M_BH–M_h relation, suggesting that there are other important parameters beyond M_h that determine M_BH, for example, the merger history of the dark matter halo. However, it is still reasonable to use M_h as a single proxy for M_BH, linking the latter to the quantities for which we have good theoretical predictions available, such as halo mass function (MF) and clustering bias. The black hole mass alone does not determine the AGN luminosity, since the particular process that fuels the black hole is also important. However, in some circumstances, for example in the “hot-halo” fueling model, the latter may lead to a tight correlation between AGN luminosity and the host halo mass.

Future wide-field X-ray surveys like eROSITA promise to give us very precise AGN clustering bias measurements as a function of redshift and luminosity (Kolodzig et al. 2013a). Since clustering bias provides direct information on how AGN populate dark matter halos, this also helps in determining the dominant X-ray AGN fueling mode: “hot-halo” vs. “cold-gas” accretion. As an example, semi-analytic models of Fanidakis et al. (2012) suggest that the hot-halo fueling mode has quite strong a dependence between the AGN X-ray luminosity and host halo mass, with brighter AGN populating preferably higher mass (and thus more strongly biased) halos. For the cold gas accretion, on the other hand, which is presumably responsible for most of the quasar activity, the host halo mass and AGN luminosity are only weakly dependent, with most of the activity occurring in somewhat lower mass halos.

Studies combining the AGN LF and their clustering strength have been carried out by several authors using mostly optically selected quasars (QSOs), such as Martini & Weinberg (2001); Haiman & Hui (2001); Wyithe & Loeb (2005); Lidz et al. (2006); Wyithe & Loeb (2009); Bonoli et al. (2010); Shankar et al. (2010b); Shankar et al. (2010a); Conroy & White (2013). The major differences between X-ray selected AGN and optically selected QSOs is that QSOs represent the brightest and the rarest end of the AGN population, while the X-ray-selected AGN, as detected in the deep X-ray surveys, are much more numerous and thus more typical examples of AGN. The other important difference is that despite QSOs being much rarer, their clustering strength is somewhat lower than the clustering amplitudes typically obtained for the X-ray selected AGN: QSO clustering is compatible with being practically independent of luminosity and having a similar strength as the dark matter halos with masses ~ 10¹² h^-1 M_⊙, while the clustering bias of X-ray selected AGN suggests that they populate group-sized halos with typical masses of ~10¹³ h^-1 M_⊙. Furthermore, the QSO activity is compatible with being driven by mergers of lower mass dark matter halos with plenty of available cold gas or by the accretion disk instabilities.

There is evidence that AGN clustering strength might scale with X-ray luminosity (Krumpe et al. 2012; Koutoulidis et al. 2012). The same trend is also seen in semianalytical galaxy–AGN co-evolution models of Fanidakis et al. (2012), where in the “hot-halo” AGN fueling mode they find quite a tight correlation between dark matter halo mass (and thus clustering bias) and AGN luminosity. It should be mentioned, however, that other authors analyzing similar samples of X-ray detected AGN, such as Coil et al. (2009) and Mountrichas et al. (2013), find no evidence of such luminosity dependence of AGN clustering. This is not surprising given the rather poor accuracy of clustering measurements in the X-ray band. Indeed, currently available samples of X-ray-selected AGN typically contain between a few hundred and a few thousand objects.

Motivated by these ideas, we assume that a relation exists between the dark matter halo mass and the AGN X-ray luminosity it harbors, and therefore, the X-ray luminosity function of AGN can be computed from the mass function of dark matter halos. In comparison to the QSO samples with typical sizes of ~100 000 objects, the currently available X-ray-selected samples are about two orders of magnitude smaller, which which suggests that these models should be kept from being too complex. Thus, we initially assume that the M_h − L_X scaling is given by a simple deterministic relation. The search for the functional form of this relation that is consistent with the concordance ΛCDM halo mass function and the observed X-ray AGN LF is presented in Sect. 2. It turns out that this way we indeed obtain an acceptable fit to the X-ray AGN LF and, as a bonus, automatically obtain a prediction about X-ray AGN clustering properties, which seems to agree with available observational data. This and other consequences of our simple model are discussed in Sect. 3. In Sect. 4 we investigate how an increased level of the model complexity influences our previous results. In particular, we extend our model by allowing a scatter in the M_h − L_X scaling relation and luminosity dependence of the AGN duty cycle. We bring our conclusions in Sect. 5.

Throughout this paper we assume a flat ΛCDM cosmology with Ω_m = 0.27, Ω_b = 0.045, h = 0.70, and σ₈ = 0.80.

2. Mapping halo MF to X-ray AGN LF

In the recent semianalytic galaxy–AGN co-evolution model of Fanidakis et al. (2012), the authors consider two main AGN fueling modes: “hot-halo” and “starburst” mode. Some of the results produced by these models, which were presented in Koutoulidis et al. (2012), are plotted in Fig. 5. Here the dark-shaded stripe represents the “hot-halo” mode, while the light-shaded region corresponds to for the sum of “starburst” and “hot-halo” modes. The starburst mode, where the SMBH fueling is induced by mergers or disk instabilities and where there is a lot of cold gas available, is arguably the one responsible for the QSO activity. On the other hand, once SMBH gets incorporated into group- or cluster-sized dark matter halo, which is filled with hot low density gas, the AGN fueling takes on a different character. In this case, in addition to metallicity, the gas cooling rate is largely dictated by the gas density and temperature, which under the assumption of hydrostatic equilibrium are directly related to the host dark matter halo mass. The gas cooling must in practice be moderated by the feedback effects from the SMBH. Thus, in the “hot-halo” fueling mode the mass of the host halo plays an important role in determining the amount of available fuel and its accretion rate, leading to the existence of a rather tight relation between the host dark matter halo mass and AGN luminosity, as illustrated by Fig. 5.

Motivated by these considerations, we start by assuming that there is a following simple scaling relation between the halo mass and AGN X-ray luminosity (1)where we take L₀ = 10⁴¹ erg/s. Having the above scaling relation, we can immediately write for the AGN X-ray LF: (2)where dn/dM_h is the concordance ΛCDM model halo MF and f_duty the “duty cycle”, in our case defined as a fraction of halos that contain AGN in its active state. For the halo MF we use the analytic form given by Sheth & Tormen (1999) and the hard band (2 − 10 keV) measurements from Aird et al. (2010) for the X-ray AGN LF. The Aird et al. (2010) LF data points in nine redshift bins are shown in Fig. 3. We assume, somewhat simplistically, that LF measurements in all luminosity bins are statistically independent with the errors following Gaussian distribution; i.e., the likelihood function gets factorized and has a simple analytic form. To sample the parametric likelihood function we used Markov chain Monte Carlo (MCMC) methods, in particular the Metropolis-Hastings algorithm (Metropolis et al. 1953; Hastings 1970). For finding the maximum likelihood point in parameter space we used the downhill simplex method as described in Press et al. (1992).

To gain an idea of the acceptable analytical form for Γ(z) and f_duty(z), we initially allow these quantities to have free values in all of the redshift bins: Γ(z_i), f_duty(z_i), i = 1...9. Thus, including M₀ we have 19 free parameters in total. Figure 1 shows best fit values along with 1σ regions for Γ(z_i) and f_duty(z_i) obtained from out MCMC calculations. We see that Γ(z_i) are quite well fitted by a simple linear form.

As a next step, we fix Γ(z) to have a linear redshift dependence Γ(z) = c₄ + c₅z and allow f_duty(z_i) to remain free parameters, so our free parameters are M₀, c₄, c₅, f_duty(z_i) (i = 1...9), i.e. 12 in total. The best fit values for f_duty(z_i) are shown with red crosses in the upper panel of Fig. 1. It turns out that these points can be approximated reasonably well with the analytic form, (3)i.e., a Lorentzian profile multiplied by z and a constant offset c₃ added⁶.

To check for the self-consistency of the above analytic description we now fix f_duty(z) to have form as given in Eq. (3), but allow the power law indices in all nine redshift bins to be free parameters, i.e., in this case we have 13 free parameters: M₀, c₁, c₂, c₃, Γ(z_i) (i = 1...9). The best fit values for Γ(z_i) are plotted as red crosses in the lower panel of Fig. 1, where with solid lines we also show the best linear and quadratic approximations. Even though the linear approximation for Γ(z) seems to capture the redshift dependence of the power law index quite well, the quadratic form with its one extra parameter, as it turns out, is statistically well justified.

Thus, our analytic MF-to-LF mapping is given by Eqs. (2), (1), and (3). In the following we call the models with linear and quadratic approximation for Γ(z) Model I and Model II, respectively. We considered models with f_duty constrained by the condition f_duty ≤ 1, as well as unconstrained models, with f_duty allowed to take any positive values. Obviously, constrained models correspond to the case when there is at most one AGN per dark matter halo which is turned on with a probability given by f_duty, whereas in the unconstrained models an arbitrary number of AGN per dark matter halo is allowed. It is interesting that in unconstrained models, f_duty does not reach values much above unity, as illustrated by the shaded region in the upper panel of Fig. 1.

In Table 1 we show the best fitting values for all the parameters of Models I and II, with and without an additional constraint on f_duty. The calculated χ² values along with expectations are given in Table 2, where we also show results for empirical models used to describe observed X-ray LFs, namely pure luminosity evolution (PLE), luminosity and density evolution (LADE), and luminosity-dependent density evolution (LDDE) models, with their analytic forms taken from Aird et al. (2010). To perform a fair comparison between different analytic fitting forms, we have recalculated the best fit models for the latter using the same fitting machinery as for our models. It can be seen that, except for PLE, all of those models provide statistically valid descriptions of the observational data. Even though our MF-to-LF mapping does not provide as good a fit as LADE or LDDE models, it is quite remarkable that this more physically motivated model indeed seems to work. As a bonus, once MF-to-LF mapping is fixed, our models make several other predictions, which are discussed in the next section.

To illustrate resulting parameter uncertainties, we show in Fig. 2 1σ and 2σ marginalized error regions along with 1D probability distribution functions for all the seven parameters of unconstrained Model II. The black dots and red crosses mark the maximum likelihood points for unconstrained and constrained (i.e., f_duty ≤ 1) cases, respectively. We see that all the model parameters seem to be sufficiently well determined. The strongest degeneracy is seen between parameters c₅ and c₆, which might leave some doubts about the usefulness of the quadratic term in Γ(z). However, as seen from the results shown in Table 2, the need for parameter c₆ is statistically justified, at a confidence level corresponding to ~2.5σ.

In Fig. 3 we plot the LFs in nine redshift bins for the best fitting PLE, LADE, LDDE models, and for our constrained Model II. Original data points with error bars, as determined by Aird et al. (2010), are also shown. Our model slightly but systematically overshoots the highest luminosity data point at high redshifts. This is discussed in more detail in the following section.

3. Implications and shortcomings of the model

In this section we investigate some of the consequences of the MF-to-LF mapping models presented above.

3.1. Comoving AGN density

We start by investigating the comoving number density of hard-band selected X-ray AGN and its dependence on redshift. In the upper panel of Fig. 4 we have plotted comoving number density for three X-ray luminosity intervals as specified in the legend. Here the highest (lowest) group of curves correspond to the lowest (highest) luminosities. Along with constrained Model II results, for comparison we have also plotted number densities resulting from Ueda et al. (2003) and Aird et al. (2010) LF fits. In the case of Ueda et al. (2003) we directly use their best fit LF parameters, whereas for the Aird et al. (2010) LDDE model, we used best fit parameters obtained in the previous section. The comparison of our results with Aird et al. (2010) LDDE parameterization should give some idea about the possible level of uncertainty in the so-obtained n(z), since both of the models give statistically valid fits to LF data.

This type of figure is often used to illustrate the antihierarchical growth of SMBHs (as manifested by the AGN activity) where the peak in number density moves toward lower redshifts as one considers intrinsically less luminous AGN. However, at least in the case of Aird et al. (2010) LF data, the conclusion about the behavior of the maximum in n(z) curve may depend on the particular choice of the functional form used to fit the X-ray LF data. Indeed, in our LF model, the n(z) maxima for the lowest luminosity bins stay almost at the same position, while for the highest luminosity bin we see gradual movement toward higher redshifts. This does not seem to be entirely incompatible with the pattern of the original LF data points in Fig. 10 in Aird et al. (2010). Furthermore, in our model the rate of change of n(z) has quite a strong dependence on luminosity, with the lowest luminosity AGN having the fastest rise in number density going from high redshifts down to z ~ 1, the AGN in the intermediate luminosity bin having only very mild increase and the highest luminosity bin showing a mild drop in the volume density with the redshift. The latter is related to the model behavior at high luminosity and redshift, discussed below.

3.2. A possible shortcoming of the model

A careful examination of Fig. 3 reveals that our model overpredicts the numbers of most luminous objects at high redshift (see the three lower panels in Fig. 3). This may point at the shortcoming of the model. Such a shortcoming could be caused by several reasons. One possible cause is the lack of proper treatment of the Eddington luminosity limit. The Eddington ratio analysis later on in this section suggests that this factor may become important in the halo mass range log M_h ≳ 13.5 − 14.5 at the redshift z ≳ 2 − 3. On the other hand, it may also be an indication of the shortcoming of the underlying halo mass function, somewhat overpredicting the abundance of massive halos at high redshifts.

Another possible reason may be a failure of the simple power-law scaling between the luminosity and halo mass, Eq. (1), at high halo masses and redshifts. Indeed, due to obvious considerations of time available for the black hole growth, the most massive halos at high redshifts may harbor somewhat less massive black holes than implied by Eq. (1). The parameters of the halo mass–luminosity scaling were determined from the global LF fit that is determined by the bulk of the LF data points and is not very sensitive to the data at the extremities because the contribution of high luminosity points to χ² is relatively small, owing to their small number and also to their relatively larger uncertainties.

Despite this shortcoming, however, the model provides a statistically valid overall description of the X-ray LF data (Table 2). Its (possible) failure at the high luminosity end only affects very few of the most massive halos at high redshifts and does not diminish the overall predictive power of the model for the bulk of typical AGN, as discussed in the following sections.

3.3. AGN clustering

Since our model has a fixed scaling relation between the halo mass and AGN X-ray luminosity, it can straightforwardly predict the AGN clustering strength.Indeed, the AGN clustering bias can be computed as (4)where the dark matter clustering bias is taken from the analytical model of Sheth et al. (2001), and the L_X integration bounds are appropriately adjusted to take the details of how the X-ray AGN sample was selected into account. In particular, for a flux-limited sample , where D_L(z) is the luminosity distance to redshift z.

In the lower panel of Fig. 4 we compare the thus-computed clustering bias as a function of redshift with the AGN clustering data from Allevato et al. (2011). Allevato et al. (2011) used soft band X-ray data from COSMOS field, which covers 2.13 deg² and contains a total of 780 AGN with available redshifts, its spectroscopic completeness being ~53%. In their analysis Allevato et al. (2011) applied a magnitude cut I_AB < 23 and redshift cut z < 4, resulting in 593 objects in their final sample, with effective spectroscopic completeness of ~65%.

Since the clustering bias in our model depends on AGN luminosity, with more luminous AGN having higher bias parameters, we have to adjust the effective flux limit appropriately in order to make a fair comparison. Since the selection of Allevato et al. (2011) sample was done in the soft band, we transformed their flux limit to the corresponding limit in the hard band. We do this in a simple way by assuming an effective population-averaged spectral index Γ⁷. The results for four different values of Γ are shown in the lower panel of Fig. 4. Because in the hard band one starts to see more and more sources that are obscured in the soft band, the effective spectral index around ~1.6 appears to be the most appropriate one. Considering this, we see that except for the highest redshift bin, our predictions for the bias parameter agree reasonably well with the measurements. The highest redshift bin is represented by one point only, and it remains to be seen how significant the deviation is. One has to keep in mind that the measurement of the linear bias parameter, and especially the proper estimation of the errors due to cosmic variance (e.g., accidental superclusters in small fields could significantly influence the results) from the fields as small as ~2 deg², may be quite a cumbersome task.

3.4. AGN duty cycle

In both constrained and unconstrained models, the AGN duty cycle (Eq. (2)) peaks at the redshift z ≃ 1 (Fig. 1). Although unconstrained models (when the duty cycle was allowed to take any positive value) result in a somewhat better description of the observed X-ray LF data, it is interesting that the peak value of the duty cycle in these models is nevertheless close to unity, with the upper limit of f_duty ≲ 1.5 (1σ confidence). This implies that on average there is at most ~one AGN per dark matter halo.

The value for the duty cycle close to unity is a consequence of the fact that at the redshift z ~ 1, the volume density of X-ray AGN is comparable to the volume density of the dark matter halos. Indeed, a power-law M_h − L_X scaling relation does not introduce any new scale, therefore the only way to obtain agreement with the observed LF is by matching the MF characteristic mass M_∗ with a characteristic AGN X-ray luminosity L_∗. In this MF to LF mapping, the role of the f_duty is to tally their volume densities. As it turns out, the comoving number densities of M_∗ (see Appendix A) and L_∗ objects are comparable at z ~ 1: ; i.e., there is approximately as many dark matter halos as AGN. Correspondingly, f_duty ~ 1. The duty cycle drops to f_duty ~ 0.2 at the redshift z = 0 (Fig. 1), a manifestation of the well known paucity of AGN in the local Universe. It is interesting to note that the duty cycle at z = 0 is not as low as predicted by the simple continuity equation models with a fixed Eddington ratio, but is closer to the ones in which the median of the Eddington ratio distribution decreases with redshift.

3.5. M_h–L_X scaling relation

As emphasized above, the crucial component in our model is the assumption about the existence of a scaling relation between the halo mass and AGN X-ray luminosity. Koutoulidis et al. (2012) have recently attempted to infer this type of scaling by using AGN clustering data from several deep X-ray fields: Chandra Deep Field South and North, the AEGIS, the extended Chandra Deep Field South, and the COSMOS field. In Fig. 5 we show their measurements with error bars along with our constrained Model II predictions for redshifts 0.5, 1.0, and 1.5, i.e. approximately covering the redshift range of AGN used in their analysis. The dotted lines with surrounding gray bands display the results from the semianalytic galaxy–AGN co-evolution model of Fanidakis et al. (2012) (as presented in Fig. 9 of Koutoulidis et al. 2012), showing their “hot-halo” (dark grey region) and “all AGN” case, i.e., the sum of “hot-halo” and “starburst” modes (light gray region). Although the uncertainties are rather large, we see that our model is able to capture the major trend where the more luminous AGN tend to populate more massive halos. The quantitative agreement between our M_h − L_X relation and results of Koutoulidis et al. (2012) is truly remarkable since they are based on totally independent arguments.

The long-dashed line in Fig. 5 shows the M_h − L_X scaling one obtains when converting the M_h–M_BH relation of Bandara et al. (2009) by assuming a constant Eddington ratio f_Edd ≡ L_X/L_Edd(M_BH) of 3 × 10^-4. In reality, because the Bandara et al. (2009)M_BH–M_h scaling relation used strong lensing masses derived from the SLACS lens sample with mean redshift z ≃ 0.2, we have adjusted the amplitude of the scaling relation to the redshift z ~ 1, using results from Croton (2009) semi-analytic models. In Bandara et al. (2009) the authors obtain with γ = 1.55 ± 0.31, which agrees quite well with the analytic and semi-analytic models of Wyithe & Loeb (2003) and Croton (2009), which in turn provide values γ = 5/3 ≃ 1.67 and 1.39, respectively. In Croton (2009) and Wyithe & Loeb (2003) models the amplitude of the scaling relation increases by factors of ~2 and ~4, respectively, while going from z = 0 to z = 1. However, the power law index γ is independent of redshift in both models.

3.6. AGN Black hole masses, Eddington ratios, and the growth of SMBHs

By assuming the M_BH–M_h scaling of Bandara et al. (2009) with amplitude adjusted according to Croton (2009) semi-analytic model, we can convert our L_X–M_h relation to f_Edd–M_h scaling relation. This is detailed in Fig. 6 where we have plotted various scalings for five different redshifts. The upper panel gives our results for L_X–M_h with 1σ uncertainties as determined in Sect. 2 for the constrained Model II. The middle panel shows the amplitude-adjusted Bandara et al. (2009)M_BH–M_h relation with 1σ uncertainties. In the lower panel we show the Eddington ratios as determined by combining the results from the two uppermost panels. The X-ray luminosities and Eddington ratios refer to the 2–10 keV band and no bolometric corrections were applied. For a typical Type I AGN spectrum, the bolometric correction factor should be in the ~5−20 range (Lusso et al. 2012).

As we can see, the Eddington ratios for high redshift (z ≳ 2) AGN are predicted to be steeply rising functions of halo mass, while at redshifts z ~ 1, i.e. close to the peak in X-ray AGN number density, the dependence of f_Edd on M_h flattens out, and at lower redshifts the trend is mildly reversed; i.e., the AGN in less massive halos have somewhat higher Eddington ratios. For z ~ 3, the 2–10 keV band f_Edd reaches values of ~0.1 for halo masses ~2 × 10¹³ M_⊙; i.e. the bolometric Eddington ratio will approach ~1 at this mass. At higher masses, our M_h − L_X scaling relation would imply f_Edd ≳ 1. Obviously, this is because our simple M_h − L_X scaling relation does not impose the Eddington luminosity limit. This may be considered as one of the shortcomings of the model, as discussed in Sect. 3.2. One has to note, however, that this does not affect the overall performance of the model, in the major part of the parameter space. Indeed, such massive halos at redshifts z ~ 3 are extremely rare, corresponding to ~4−5σ peaks in the initial fluctuation field. The bulk of the objects are shining at much lower Eddington ratios, ≪1.

Having specified M_BH–M_h scaling relation, we can also calculate how SMBHs with different masses grow throughout the cosmic time. Taking the duty cycle into account, the population-averaged e-folding time of the SMBH growth is (5)where η is the accretion efficiency, defined as (6)L_X refers to the X-ray luminosity in the 2–10 keV band used throughout this paper, and f_bol ~ 5−20 is the bolometric correction factor for the 2–10 keV band. The e-folding time can be expressed further though the Eddington ratio: (7)In the lower panel of Fig. 7 we plot the e-folding time in units of the Hubble time, i.e. τH(z), vs redshift for several values of the black hole mass. Shown in the upper panel is the evolution of the population average 2 − 10 keV band X-ray luminosity with redshift. AGN harboring black holes with masses less than a ~few× 10⁸ M_⊙ are too faint at high redshifts and fall outside the luminosity range of the X-ray LF data of Aird et al. (2010) used in this study. Therefore they are not directly constrained by the X-ray data. At these redshifts we only show with dotted lines the extrapolations of the central model curves and omit the 1σ uncertainty regions.

We see that the X-ray luminosity of more massive SMBHs declines faster with the redshift than for less massive ones. For ~10⁸ M_⊙ black holes, L_X stays almost constant, while for even lower mass SMBHs, the evolutionary trend is slightly reversed; i.e., the objects are getting somewhat brighter as one moves toward lower redshifts. As one might have expected from the results plotted in the bottom panel of Fig. 6, for the typical X-ray-selected AGN, the corresponding SMBH e-folding times are long in comparison to the Hubble time. With the exception of most massive SMBHs at high redshifts, the typical X-ray-selected AGN are not in the rapid mass growth regime in the z ~ 0−3 redshift range.

The crossing of the model curves in the lower panel of Fig. 7 at the redshift of z ~ 0.85 is artificial and does not seem to have any physical meaning. It is driven by the equality of the power law indices of the M_BH − M_h and M_h − L_X scaling relations, which leads to the τH being independent of M_BH at that particular z.

Thus, with the particular form of the M_BH-M_h scaling relation from Bandara et al. (2009) with amplitude evolution from Croton (2009), our model implies relatively low Eddington ratios, f_Edd ~ 10^-4−10^-3. In the redshift range z ~ 0−3, for a typical X-ray-selected AGN, accretion proceeds in the “hot-halo mode” characterized by a relatively low mass accretion rate and long e-folding time for SMBH mass growth. Obviously, this conclusion critically depends on the assumed M_BH − M_h scaling relation. Indeed, as X-ray luminosity for the given halo mass is directly determined by our model from the L_X − M_h mapping, the Eddington ratio scales inversely with the black hole mass, while e-folding time is linearly proportional to the black hole mass.

4. Extensions of the model

The minimalistic model considered in the previous sections has two major simplifications: a deterministic M_h − L_X relation and a luminosity-independent duty cycle. In this section we lift these assumptions and explore a class of more sophisticated models.

We start with introducing scatter in the M_h − L_X relation. We assume that scatter can be described by the lognormal probability distribution⁸: (8)whose width σ is assumed to be independent of M_h and z. Eq. (2) for the luminosity function is now replaced by (9)There are several obvious possibilities for μ(M_h,z), depending on whether one chose the (i) mode; (ii) median; or (iii) mean of the lognormal probability distribution function to obey our previous deterministic M_h − L_X relation defined by Eq. (1). This gives for the mode, median, and mean, respectively. Since this analysis is largely inspired by comparison with observations, which pick up the most probable objects, it is natural to choose the mode of the distribution, Eq. (10).

It is easy to see that in case of small scatter, all the above descriptions are equivalent, and in the limit σ → 0 and we recover the deterministic relation given in Eq. (1).

To calculate the AGN clustering bias, Eq. (4) now generalizes to (13)Similar to Eq. (4), the integration limits for L_X have to be adjusted appropriately, taking the ways AGN sample was selected into account.

In Fig. 8, which is the analog of the previous Fig. 2, we show the results of introducing scatter in the M_h − L_X mapping. The scatter is parametrized via a single parameter – the width of the log-normal distribution σ ≡ c₈. All other conventions are exactly the same as in Fig. 2. The best fitting values for all the parameters are given in the first row of Table 3 for the case of unconstrained f_duty. As is obvious from Fig. 8 and Table 3, no scatter is required by the X-ray LF data, since the best-fit value of σ is close to zero. Correspondingly, the χ² and the best fit values of other parameters are nearly identical to the ones obtained in the deterministic case (cf. Table 1). By comparing with the results shown in Fig. 2, it is interesting to note that after allowing for the scatter in the M_h − L_X relation, no new strong parameter degeneracies appear and all parameters are as well determined as before.

It is also clear from Fig. 8 (bottom row) that moderate scatter, up to σ ~ 0.5 is still consistent with the LF data, within ~1σ confidence. To investigate models with even greater scatter, we also performed two MCMC runs where the scatter was fixed to 0.25 and 0.50 dex (i.e., σ = 0.25ln10 and σ = 0.50ln10, respectively). The best parameter values for the case of unconstrained f_duty for these two runs are shown in the second and third rows in Table 3. It turns out that the data and model are able to comfortably accommodate M_h − L_X scatter up to ~0.50 dex.

We now investigate the effect of the luminosity dependence of the duty cycle. One possibility of such dependence was suggested by observations (Bundy et al. 2008), as well as by semi-empirical models (Merloni & Heinz 2008; Shankar et al. 2013). We consider L_X-dependence in the power law form: (14)where f_duty(z) is given by Eq. (3) and is fixed at erg/s. As before, we consider the case of unconstrained f_duty.

Results of these calculations are presented in the last two rows of Table 3, for cases without and with scatter in the M_h − L_X relation. The confidence contours for the model parameters are shown in Fig. 9, where for comparison we have also plotted our contours from Fig. 8. Introduction of the L_X −dependence of the duty cycle dramatically improves the quality of the fit with the χ² reducing by more than ~30−40%. Including the scatter in the M_h − L_X relation does not lead to significant further improvement, although statistically it is justified at the ~1.5−2σ level. Even though the best fit value of the scatter σ is quite high (~0.6 dex), its 2σ contour includes the null value.

Contrary to what one might have guessed, the best fitting values for α turn out to be positive, i.e., f_duty grows with L_X⁹. It is interesting to note that the L_X-independent part of the duty cycle describing its redshift dependence, f_duty(z), (Eq. (3)), is sufficiently robust and insensitive to the details of the model; it is nearly identical in shape and normalization to the f_duty(z) in our default model in previous sections. In particular, it is insensitive to the choice of , with changes in mostly absorbed by corresponding changes in other parameters, mostly by the Γ(z) in Eq. (1). Changes in the latter lead to modifying the M_h − L_X relation, such that the characteristic halo mass M_∗ is no longer mapped into the characteristic luminosity L_∗. For this reason, the overall duty cycle f_duty(L_X,z) can now take values over unity, for example, f_duty(L_X = 10⁴⁶,z = 1) ~ 6. This should be interpreted as meaning that high mass halos may contain several SMBHs that are in the active state. The robustness and potential implications of this result still need to be explored.

The statements of the previous section regarding the low Eddington ratios of X-ray-selected AGN and long mass growth times hold for the extended models, as shown in Fig. 10. Here the solid lines and the shaded areas show distribution medians and regions between 16% and 84% quantiles, respectively. The strong level of asymmetry in these probability distributions can be seen clearly by looking at the dotted lines, which represent the distribution modes.

The shift of the M_h − L_X mapping toward higher halo masses in the models with luminosity-dependent f_duty(L_X,z), leads to an increase in its predictions for the clustering bias at z ≳ 2. This places the new models in somewhat better agreement with the observational data, as demonstrated in Fig. 11. The figure compares predictions of the extended models with our default model and with results of clustering bias measurements. The lowest stripe shows the results from Fig. 4, while the middle and upper ones correspond to the models with luminosity-dependent f_duty(L_X,z), discussed in this section. While at lower redshifts, z ≲ 1.5 − 2, the behavior of different models is nearly identical and equally consistent with the data. At z ≳ 1.5 they start to diverge, with the extended models in better agreement with the highest redshift data point. This suggests a possible way to distinguish between different models, using the clustering bias measurements at redshifts z ≳ 2.

Comparing the χ² values shown in Tables 2 and 3, we note that models with a luminosity-dependent duty cycle give a better fit to the hard band X-ray AGN LF than do the specifically tailored and commonly used mathematical forms like LADE or LDDE. Therefore these models can be used as an easily calculable analytic representation of the LF. To assist in computing the luminosity function, in Appendixes A and B we give a detailed description of the algorithm and provide easy-to-use approximation formulae.

5. Conclusions

Even though dark matter halos are highly nonlinear objects, with the aid of cosmological N-body simulations along with heuristic analytic models, we now have good knowledge about their redshift dependent MF and clustering bias. In CDM type cosmologies, where structures form from the bottom up, one may expect that the buildup of more massive objects has also been accompanied by the formation of more massive SMBH in their centers. This should also be manifested on average by more luminous AGN. Although there certainly are other key parameters beyond the halo mass, which influence the central SMBH mass, one may expect that M_h is one of the main contributors.

In this paper, we assumed that a scaling relation between X-ray AGN luminosity and its host halo mass does exist. Assuming further that the dark matter halo mass function is described by the concordance ΛCDM model, we can predict X-ray AGN luminosity function. Comparing the latter with the redshift-dependent AGN X-ray luminosity function known from observations we can determine the shape and parameters of the M_h − L_X scaling relation.

Our main conclusions are the following:

• We have shown that a simple scaling relation in the form of Eq. (1), where Γ(z) is a linear or quadratic polynomial, allows one to obtain X-ray AGN LFs consistent with observational data. The best fitting model parameters along with χ² values are summarized in Tables 1 and 2.

• Since the M_h − L_X scaling establishes a link to the halo MF, our model is certainly more predictive than the usual multi-parameter fitting forms tailored to fit the X-ray LF data. In particular, our model predicts the redshift and luminosity dependence of the X-ray AGN clustering bias, which agrees reasonably with observational data.

  	Fig. A.1 Comparison of the LF fit from Appendix B with the “new model” and LDDE model of Fig. 3. In addition, with dotted lines we have shown the model with the lowest χ2 value obtained throughout this paper, i.e., the model with scatter in Mh − LX scaling relation along with LX-dependent duty cycle (the last entry of Table 3).
  Open with DEXTER

• The obtained M_h–L_X scaling relation is in good agreement with the data presented in Koutoulidis et al. (2012) and with semianalytic galaxy–AGN co-evolution models of Fanidakis et al. (2012). Comparison with semianalytical models suggests that for X-ray AGN, the dominant accretion mode is the “hot-halo mode”, in contrast to the “starburst mode” which is compatible with being a dominant accretion mode for somewhat less clustered but more energetic quasars.

• The AGN duty cycle peaks at the redshift z ≃ 1, and its value at the peak is consistent with unity, implying that there is at most one AGN per dark matter halo. The value for the duty cycle close to unity is a consequence of the fact that at the redshift z ~ 1, the volume density of X-ray AGN is comparable to the volume density of the dark matter halos: .

• We further combined our M_h-L_X relation with the M_h–M_BH relation from Bandara et al. (2009), along with its redshift evolution taken from Croton (2009), and obtained Eddington ratios f_Edd ≡ L_X/L_Edd ~ 10^-4 − 10^-3 (2–10 keV band, no bolometric correction applied). This result presents another evidence that at the redshifts below z ≲ 3, SMBHs in X-ray selected AGN typically grow through “hot-halo” accretion mode characterized by low mass accretion rate and a long e-folding time. It is also consistent with the fact that according to their clustering strengths, X-ray selected AGN populate group-sized dark matter halos and thus are typically embedded in halos of hot low density gas.

• We also investigated a class of extended models that include a scatter in the M_h–L_X scaling relation and a luminosity-dependent f_duty. While the scatter alone is not required by the X-ray LF data, including the luminosity dependence of the duty cycle leads to significant improvement of the fit quality. In the new models, the statements regarding the low Eddington ratios of X-ray-selected AGN and their long SMBH growth times hold. However, they differ in their predictions regarding the behavior of the clustering bias at redshifts z ≳ 2, with extended models in better agreement with the most recent bias measurements at z ~ 2.2 by Allevato et al. (2011).

• In Appendix B, we summarize the formulae that can be used to construct an analytical model of the hard-band AGN LF. From the statistical point of view, it provides a better description of the X-ray LF data as specially tailored analytical forms commonly used in the literature.

NOTE: during the refereeing process, a new paper by Fanidakis et al. (2013) appeared, which touches similar topics to the ones discussed in the current study.

---

Acknowledgments

We thank James Aird for providing the tables of hard-band X-ray luminosity function measurements and our referee for suggestions that greatly helped to improve the paper.

References

Appendix A: Analytic fits for concordance ΛCDM halo MF and bias parameters

To facilitate a fast check of our results for those who do not have numerical routines for the halo MF and clustering bias available, we provide here analytic fits for those quantities, valid for the mass and redshift ranges relevant for this study: M_h ≳ 10¹² M_⊙, z < 3 − 4.

Halo MF can be approximated as (A.1)where Similarly, for the halo clustering bias, as a function of halo mass and redshift, one approximately obtains

where Under the mass and redshift constraints as given above, those analytic forms provide fits to the Sheth & Tormen (1999) MF and Sheth et al. (2001) bias parameter with accuracies ~ 10% and better than 5%, respectively. In the above formulae, we have assumed that MF is measured in units of h³Mpc^-3/h^-1M_⊙ and halo mass in h^-1M_⊙, while in the main part of our paper, we fixed the reduced Hubble parameter to h = 0.7.

Appendix B: A simple physically motivated fit to the observed hard band X-ray AGN LF

The hard band X-ray AGN LF can be approximated as (B.1)where f_duty(L_X,z), M_h(L_X,z), and are given by Eqs. (14,3), (1), and (A.1), respectively. Using the approximation for the MF from Appendix A leads to slightly different best fitting parameter values compared to the ones shown in the last row of Table 3: {c_i} _{i  =  1...9} = {0.3654 ,0.8243,0.2402,1.002, –0.2997, 0.03216, 26.97, 0., 0.3198}, which results in χ² = 47.5 for 71 degrees of freedom.

A detailed comparison of the above LF fit with the “new model” and LDDE model of Fig. 3 is shown in Fig. A.1.

All Tables

All Figures