The halo occupation distribution of active galactic nuclei Free

Abstract

Using a fully cosmological hydrodynamic simulation that self-consistently incorporates the growth and feedback of supermassive black holes and the physics of galaxy formation, we examine the effects of environmental factors (e.g. local gas density and black hole feedback) on the halo occupation distribution of low-luminosity active galactic nuclei (AGN). We decompose the mean occupation function into central and satellite contribution and compute the conditional luminosity functions (CLFs). The CLF of the central AGN follows a lognormal distribution with the mean increasing and scatter decreasing with increasing redshifts. We analyse the light curves of individual AGN and show that the peak luminosity of the AGN has a tighter correlation with halo mass compared to instantaneous luminosity. We also compute the CLF of satellite AGN at a given central AGN luminosity. We do not see any significant correlation between the number of satellites with the luminosity of the central AGN at a fixed halo mass. We also show that for a sample of AGN with luminosity above 10⁴² erg s⁻¹ the mean occupation function can be modelled as a softened step function for central AGN and a power law for the satellite population. The radial distribution of AGN inside haloes follows a power law at all redshifts with a mean index of −2.33 ± 0.08. Incorporating the environmental dependence of supermassive black hole accretion and feedback, our formalism provides a theoretical tool for interpreting current and future measurements of AGN clustering.

1 INTRODUCTION

Through recent observations it has been shown that active galactic nuclei (AGN) play a significant role in the evolution of galaxies. The observed correlation between the mass of the central black hole and the velocity dispersion of the bulge of the host galaxy suggests a strong connection between galaxy evolution and black hole activity (e.g. Gebhardt et al. 2000; Merritt & Ferrarese 2001; Tremaine et al. 2002; Graham et al. 2011). Several theoretical models relating AGN activity or black hole growth to galaxy evolution have been proposed (e.g. Soltan 1982; Silk & Rees 1998; Salucci et al. 1999; Kauffmann & Haehnelt 2000; Wyithe & Loeb 2003; Marconi et al. 2004; Shankar et al. 2004; Di Matteo, Springel & Hernquist 2005; Cattaneo et al. 2006; Croton et al. 2006; Hopkins et al. 2006; Lapi et al. 2006). Differentiating between these theoretical models requires several observational quantities. Measurements of the AGN luminosity function (e.g. Boyle et al. 2000; Fan et al. 2001) and the number density of black hole hosts in the present Universe can provide an estimate of the duty cycle of black holes. Alternatively, the black hole mass function (e.g. Shankar et al. 2004; Graham & Driver 2007) measured at the current epoch can provide constraints on models of black hole growth. Measurement of AGN clustering provides a unique way to study the physical characteristics of AGN through the connection with their host dark matter haloes (e.g. Croom et al. 2004; Gilli et al. 2005; Myers et al. 2006; Ross et al. 2009; Shen et al. 2009).

Theoretically, clustering properties of AGN have been mostly studied with semi-analytic models using the halo model or the black hole continuity equation approach (e.g. Lidz et al. 2006; Bonoli et al. 2009; Shankar et al. 2010a). Although semi-analytic models have provided the formalism for interpreting quasar clustering, there are certain issues that these models cannot approach. In this formalism, physical parameters (e.g.‘quasar duty cycle’) are treated as free parameters that are constrained from clustering measurements (e.g. Martini & Weinberg 2001). There also exist degeneracies between parameters in the models (e.g. Shankar, Weinberg & Shen 2010b). Moreover, in most of these models, black holes are accreting at a fixed fraction of the Eddington rate, making accretion rate dependent on black hole mass only. In reality, parameters such as accretion efficiency and the duty cycle are not fixed and depend strongly on environment (e.g. local gas density, mergers and feedback from AGN). Some recent semi-analytic studies started to consider more general cases of black hole accretion to investigate the cosmological co-evolution of black holes and the evolution of AGN luminosity functions (e.g. Malbon et al. 2007; Marulli et al. 2008, 2009; Menci et al. 2008).

Hydrodynamic simulations of galaxy formation with black hole growth can capture the environmental dependence of accretion. The growth of black holes is tied with the full dynamics of dark matter. Also, with these simulations we can model the interplay between AGN and galaxy formation in the form of feedback and self-regulated growth (e.g. Wyithe & Loeb 2003). Thus, these simulations provide an excellent platform to study the co-evolution of AGN with large-scale structures over cosmological time-scales. Studies of AGN clustering using cosmological simulations have been carried out by Thacker et al. (2009) and Degraf, Di Matteo & Springel (2011a). Degraf et al. (2011a) calculate the correlation function of black holes using a smoothed particle hydrodynamics (SPH) cosmological simulation that incorporates galaxy formation physics and self-consistent black hole growth and feedback (Di Matteo et al. 2008). They show that the black hole correlation function consists of two distinct components: contributions from intra-halo and inter-halo black hole pairs, i.e. the one-halo and two-halo terms. At small scales the one-halo term still follows a power law, which is different from that of dark matter. This boost in small-scale power is due to galaxies hosting multiple black holes.

Despite their advantages, studies using hydrodynamical simulations are limited due to their computational costs. For example, simulations are limited to small boxes (Di Matteo et al. 2008), restricting the analysis to small scales and low-luminosity AGN. The computational cost of running these simulations to z= 0 limits the possibility to compare the outputs with observations in the local universe. Also, since every modification of the recipes adopted to describe the ‘sub-grid’ physics requires rerun of the simulations, it is not feasible to do systematic parameter studies. Semi-analytic studies are advantageous over hydro simulations in this context. We now employ the SPH simulation of Di Matteo et al. (2008) to investigate the relation between AGN and their host haloes within the halo occupation distribution (HOD) framework, a useful analytic formalism for modelling and interpreting AGN clustering. In our approach, the expensive hydrodynamic simulation is used to obtain insights on a general analytic technique for studying AGN clustering. The HOD (e.g. Ma & Fry 2000; Seljak 2000; Berlind & Weinberg 2002) provides an analytic formalism for understanding clustering properties of galaxies or, in this case, AGN. It is characterized by the probability P(N|M) that a halo of mass M contains N AGN of a given type, together with their spatial and velocity distribution inside haloes. As long as the HOD at a fixed halo mass is statistically independent of the large-scale environments of haloes (e.g. Bond et al. 1991), it provides a complete description of the relation between AGN and haloes, allowing the calculation of any clustering statistics (e.g. two-point correlation function, three-point correlation function, void probability distribution and pairwise velocity distribution) at all scales (small, intermediate and large) for a given cosmological model. This implies that if we can constrain the HOD empirically, we will have knowledge of everything that the measured clustering properties have to tell us about AGN formation and evolution models. The HOD allows us to distinguish between the background cosmology and AGN evolution models. The cosmological parameters are encoded in the distribution of haloes, whereas the bias between mass and AGN is fully described in the probability distribution P(N|M) (see Berlind & Weinberg 2002 and Zheng & Weinberg 2007 for the strength of HOD). The HOD formalism has been widely used in interpreting galaxy clustering (see Zehavi et al. 2005).

In this paper we perform a theoretical study on the relation between AGN and dark matter haloes using cosmological simulations. We assess the validity of several simplifying assumptions in the HOD modelling of AGN and its predictive power and limitations for modelling AGN clustering data. Our paper is organized as follows: In section 2 we describe our simulations and give a brief description of our methodology. In section 3, we present the conditional luminosity function (CLF) of AGN, the mean occupation distribution, and the radial profiles of AGN. Finally, in section 4, we summarize our main results and provide future directions to this work.

2 SIMULATION

2.1 Numerical code

The numerical code uses a standard Λ cold dark matter (ΛCDM) cosmological model with cosmological parameters from the first-year Wilkinson Microwave Anisotropy Probe (WMAP) results (Spergel et al. 2003). The simulation uses an extended version of the parallel cosmological tree particle mesh-SPH code gadget2 (Springel 2005). Gas dynamics is modelled with Lagrangian SPH (Monaghan 1992); radiative cooling and heating processes are computed from the prescription given by Katz, Weinberg & Hernquist (1996). The relevant physics of star formation and the associated supernova feedback have been approximated based on a sub-resolution multiphase model for the interstellar medium developed by Springel & Hernquist (2003). The size of the simulation box is 33.75 h⁻¹ comoving Mpc with periodic boundary conditions. A detailed description of the implementation of black hole accretion and the associated feedback model is given in Di Matteo et al. (2008). Black holes are represented as collisionless ‘sink’ particles that can grow in mass by accreting gas or by mergers. The Bondi–Hoyle relation (Bondi 1952) is used to model the accretion rate of gas on to a black hole and capture the environmental dependence of black hole accretion. The accretion rate is given by ⁠, where ρ and c_s are density and speed of sound of the local gas, respectively, v is the velocity of the black hole with respect to the gas, and G is gravitational constant. Although the Bondi parametrization assumes spherical accretion, observations show that the Bondi relation adequately captures the physical state of the black hole and the Bondi scaling holds at scales much larger than the Bondi radius (e.g. Allen et al. 2006).

The bolometric luminosities of the AGN are ⁠, where η= 0.1 is the canonical efficiency for thin disc accretion. It is assumed that a small fraction of the bolometric luminosity couples to the surrounding gas as feedback energy E_f such that ⁠, with the feedback efficiency ε_f taken to be 5 per cent. Based on previous studies of galaxy mergers, the value of the feedback efficiency parameter ε_f is chosen to be 5 per cent in the simulation, a necessary fraction to reproduce the normalization of the observed relation at z= 0 (Di Matteo et al. 2005). This feedback energy is put directly into the gas-smoothing kernel at the position of the black hole (Di Matteo et al. 2008). The feedback energy is assumed to be distributed isotropically for the sake of simplicity; however, the response of the gas can be anisotropic (e.g. Chatterjee et al. 2008; Di Matteo et al. 2008). The formation mechanism for the seed black holes which evolve into the observed supermassive black holes today is not known. The simulation creates seed black holes in haloes that cross a specified mass threshold. At a given redshift, haloes are defined by a friends-of-friends group finder algorithm run on the fly. For any halo with mass M≥ 5 × 10¹⁰ h⁻¹ M_⊙ that does not contain a black hole, the densest gas particle is converted to a black hole of mass M_BH= 5 × 10⁵ h⁻¹ M_⊙. The choice of the seed mass in the present simulation was based on current models of seed black hole formation in the Universe (e.g. Bromm & Loeb 2003; Bromm & Larson 2004). The black hole then grows via the accretion prescription given above and by mergers with other black holes (Di Matteo et al. 2008). The simulation parameters are shown in Table 1.

The detailed parameter studies of this simulation, and comparison with several observations have previously been done by Sijacki et al. (2007), Di Matteo et al. (2008), Colberg & Di Matteo (2008), Croft et al. (2009) and Degraf, Di Matteo & Springel (2010). The model has reproduced the observed relation, total black hole mass density (Di Matteo et al. 2008) and the AGN luminosity functions and their evolution in optical, soft and hard X-ray band (Degraf et al. 2010). The black hole mass density from the simulation matches well with the constraints from the integrated X-ray luminosity function (Marconi et al. 2004; Shankar et al. 2004) and the accretion rate density is consistent with the constraints of Hopkins et al. (2007a). This model has also been used in galaxy merger simulations to investigate black hole growth and their correlation with host galaxies (Di Matteo, Springel & Hernquist 2005; Robertson et al. 2006; Hopkins, Richards & Hernquist 2007b).

The model has some intrinsic limitations. It does not treat the physics of the accretion disc in detail, which is not feasible in a simulation with cosmological volume. Nevertheless, the model is capable of representing some key aspects of black hole evolution in a cosmological context and in reproducing many observational results. Although the model has some limitations, it has been fairly successful in reproducing some of the key observational results. Several other teams (Johansson, Naab & Burkert 2008; Booth & Schaye 2009; Teyssier et al. 2011) have now implemented similar modelling for black hole accretion and feedback in hydrodynamic simulations. Booth & Schaye (2009) have independently explored the parameter space of the fiducial model of Di Matteo et al. (2008). Teyssier et al. (2011) have implemented this prescription in an adaptive mesh refinement code (different from the SPH formalism) and showed that AGN feedback can solve the overcooling problem in clusters. The large number of previous investigations make this particular sub-grid model a good choice for studying co-evolution of AGN with their host dark matter haloes. We emphasize that in this work we do not intend to explore the parameter space related to the sub-grid modelling of black holes. Instead we study the relation between AGN and dark matter haloes in one simulation box. The goal of this study is to obtain a useful description between AGN and dark matter to guide the interpretation of observations (e.g. AGN clustering) and to learn about AGN evolution using the model of Di Matteo et al. (2008).

2.2 Sample selection

For our analysis we consider haloes above a mass scale of 10¹¹ M_⊙ and black holes with masses above 10⁶ M_⊙. This choice of mass scales will minimize the effect of the seed mass in our results. We also restricted our analysis to black holes with Eddington-scaled accretion rate since we want to study systems with active accretion. This selection imposes a luminosity cut of 10³⁸ erg s⁻¹ in our sample. One of the key goals of this paper is to separate the contributions to the occupation distribution from the central and satellite AGN (as done for galaxy HOD by Zheng et al. 2005). We consider all the black holes within R₂₀₀ (defined as the radius within which the enclosed mean density is 200 times the critical density) of a halo to be associated with the host halo. The most massive black hole within R₂₀₀ is assumed to be the central AGN and the rest are designated as satellites. At z= 1.0, 3.0 and 5.0, about 79, 93 and 91 per cent of the central black holes lie within 20 per cent of R₂₀₀, respectively. We expect the central black hole to be located at the centre of the main galaxy in the halo since this is the place of highest gas density and thus maximal black hole growth. The fact that most ‘central AGN’ are found near the group centre shows that our technique for identifying the central AGN is reasonably adequate. Due to the limited volume of our simulation we could only probe the faint end of the AGN luminosity function (Degraf et al. 2010). For calculating the mean occupation, we study AGN samples with bolometric luminosities above 10⁴² erg s⁻¹, corresponding to the observational limit of X-ray-selected AGN in deep surveys (e.g. Boyle et al. 1993; Cowie et al. 2003; Ueda et al. 2003). We also show the distribution of AGN with luminosities above 10⁴⁰ erg s⁻¹ and 10³⁸ erg s⁻¹ for theoretical interests to see the luminosity dependence of HOD parameters.

In a companion paper, Degraf et al. (2011b, hereafter Paper I), we use the same SPH simulation to characterize the HOD of black holes as a function of black hole mass. Our present study is fundamentally different from that of Paper I. We aim to study the HOD of AGN in terms of AGN luminosity. By definition the HOD is always studied based on physical properties of the object. In the case of galaxy clustering, the galaxy HOD has been extensively studied in the literature based on colour, luminosity and stellar mass (see Zehavi et al. 2005; Zheng, Coil & Zehavi 2007). Black hole mass and luminosity (accretion rates) are different parameters and there does not exist an obvious one-to-one correspondence between mass and accretion rates (e.g. Fig. 1). Additional baryonic physics is needed to go from black hole mass to accretion rate (and thus luminosity). There is some correlation between mass and accretion rates, but the scatter is huge and the correlation depends strongly on redshift (Chatterjee et al. 2008; Colberg & Di Matteo 2008).

In practice, it is barely possible to have accurate mass measurement of a large sample of black holes beyond the local Universe. Therefore, the study in Paper I is mostly of theoretical interest. By investigating the occupation properties of black holes at the subhalo level, Paper I develops an analytic mechanism to populate dark matter haloes with black holes in simulation. In contrast, AGN luminosity is more readily measured observationally and there are existing and ongoing efforts in measuring AGN/quasar clustering as a function of luminosity. Our study in this paper will provide the necessary ingredient in interpreting these observations. Similar to the mass function of black holes and luminosity function of AGN, the black hole-mass-based HOD and AGN luminosity-based HOD are complementary tools to advance our understanding of black hole evolution.

Since the goal of this paper is to provide the HOD framework to interpret observations, we adopt some technical differences in the definition of dark matter haloes and the choice of our sample from Paper I. In Paper I, we differentiate between central and satellite black holes at the subhalo level. Black holes residing in the central subgroup are called central black holes and black holes residing in satellite subgroups are called satellite black holes. In the nomenclature in Paper I there can be multiple central black holes. We emphasize that the model in Paper I will be useful for populating haloes (N-body simulations) with black holes using semi-analytic approaches and for studying the distribution of black holes in cosmological simulations. On the other hand, our approach in this paper is to perform a theoretical study on the relation between AGN (selected based on luminosity) and dark matter haloes and examine the assumptions in the HOD for modelling AGN clustering data. Our identification of central (i.e. most massive black hole within R₂₀₀) and satellite (rest of the black holes within R₂₀₀) AGN is equivalent to primary and secondary black holes in Paper I (see fig. 2 of Paper I) and is similar to the standard terminology used for galaxy HOD. We note that the friends-of-friends haloes were used in Paper I for studying black hole occupation. Here we choose R₂₀₀ to define the halo boundary, which is commonly adopted in defining galaxy clusters. Since the HOD describes bias at the level of systems near dynamical equilibrium, traditional virial methods of mass estimation can take advantage of this and provide a direct measurement of the probability distribution P(N|M) (Berlind & Weinberg 2002). We have verified that the two different halo definitions show a good agreement in assigning black holes to host haloes within 10 per cent.

3 RESULTS

In this section we present the distribution of accretion rates as a function of halo mass and black hole mass. From this we calculate the CLFs and the mean occupation distribution of central and satellite AGN. Finally we show the radial distribution of the satellite AGN and provide the best-fitting HOD parameters.

3.1 The accretion rates

Fig. 1 shows the distribution of the Eddington-scaled accretion rate (⁠⁠), as a function of halo mass for three redshifts. The black points are central AGN and the red open circles represent satellites. The mean rate is roughly independent of halo mass, which is commonly assumed in semi-analytic studies. However, the mean value of χ at a given halo mass varies with redshift. AGN at higher redshifts (left panel of Fig. 1 showing z= 5.0) tend to accrete at a higher value than their low-redshift counterparts (right panel of Fig. 1 showing z= 1.0). Fig. 2 shows the distribution of χ as a function of black hole mass. It is defined as the conditional probability distribution of χ for a given black hole mass P(χ|M_BH). The solid and the dotted lines denote central and satellite AGN, respectively. The mean rates are again roughly independent of black hole mass but vary substantially with redshift. The central and satellite distributions tend to trace each other at all redshifts.

In Fig. 3 we present the spatial distribution of the satellite AGN. The sizes of the circles are proportional to the black hole mass (left panel), gas mass (middle panel) and stellar mass (right panel), respectively. The gas and the stellar mass are computed within a spherical region of radius 10 kpc surrounding the black hole. The green triangles represent black holes with zero stellar mass within 10 kpc. The top, middle and bottom panels are for z= 1.0, 3.0 and 5.0, respectively. The AGN accretion rates do not show any dependence on radius or black hole mass. Also, we do not see any correlation between black hole mass and its location within the halo. The trends of accretion rates with gas and stellar mass are more prominent. The low values of χ at low redshift do not correlate with the mass distribution of the black holes (left panel of Fig. 3). The mass distribution of black holes is not significantly different at different redshifts. At low redshifts the gas mass around black holes decreases and so do the accretion rates. This is due to the effect of AGN feedback. Outflow from AGN is responsible for expelling gas from the vicinity of the black hole and suppressing its own growth. AGN feedback is responsible for shutting down star formation too. In Fig. 3 we also see the evidence of suppressed stellar mass in the immediate neighbourhood of the AGN. We see evidence of how a local phenomenon (feedback) affects the global distribution of AGN, regardless of the host halo mass (see Fig. 1 for the host halo masses). This picture is consistent with what has been seen in previous studies (e.g. Di Matteo et al. 2008; Sijacki et al. 2007) and semi-analytic models (e.g. Wyithe & Loeb 2003)

3.2 The conditional luminosity function

where ϕ(L) is the AGN luminosity function and dn/dM is the halo mass function. The CLF constitutes the differential form of the HOD. The CLF provides a tool to examine the distribution of halo mass for a given luminosity and study luminosity-dependent clustering. This formalism has been widely used in modelling galaxy clustering with data from 2dF Galaxy Redshift Survey (2dFGRS) and DEEP2 redshift surveys (e.g. van den Bosch, Yang & Mo 2003). The CLF is useful in generating mock catalogues, which can then be used to test the cosmological model (e.g. Mo et al. 2004; Yan, White & Coil 2004). Fig. 4 shows the CLF in bins of halo mass running from log M= 11.25 ± 0.25 (leftmost panel) to log M≥ 12.5 (rightmost panel). The dashed lines show the distribution of bolometric luminosities for the central AGN while the dotted lines show the distributions for the satellite AGN. In both cases the distribution has been normalized by the total number of haloes in the given mass bin. The top, middle and bottom panels show results for redshifts 1.0, 3.0 and 5.0, respectively.

The luminosity distributions reveal several features. The central AGN luminosity distribution closely traces a lognormal distribution (e.g. Martini & Weinberg 2001). Fig. 5 captures the features we observe in the CLF. The top panel of Fig. 5 shows the evolution of the mean log luminosity of the central AGN as a function of halo mass. The bottom panel of Fig. 5 shows the 1σ scatter in the log of the luminosity. The mean luminosity at a given halo mass evolves with redshift, with higher mean luminosity at higher redshifts. The scatter in the luminosity distribution of central AGN increases towards lower redshift, indicating that the luminosities for a given halo mass are more uniform at high redshift than at low redshift. As a result of this, at low redshifts, AGN samples based on luminosity will have a wider range of host halo masses and hence clustering will depend weakly on luminosity. However, at high redshifts luminosity-dependent clustering will be more prominent. Similar redshift evolution of luminosity-dependent clustering has been observed with SDSS quasars (Shen et al. 2009, 2007). We note that the clustering properties of AGN will be largely related to their host halo mass. However we find that accretion rates and hence luminosity depends strongly on local properties particularly at low redshifts (e.g. feedback discussed in Section 3.1), which erases some of the dependences of bolometric luminosity on halo mass and hence the luminosity-dependent clustering at low redshifts. This will not be the case in black hole mass-selected sample since black hole mass is more tightly correlated with halo mass regardless of redshifts (Colberg & Di Matteo 2008). The shape of the satellite luminosity function cannot be identified definitively due to lack of statistics. We note that the lower end of the satellite luminosity function is affected by the resolution of the simulation and the seed black hole mass, which is manifested in the cut-off of the luminosity function. The satellite distribution also shows some variation with halo mass. For a given halo mass the peak of the satellite distribution tends to be at a lower luminosity than the central AGN.

To examine the effect of the central AGN on the number distribution of satellites within a halo, we compare the CLFs of satellite AGN for haloes differing in central AGN luminosities. This is defined as the conditional distribution of satellite AGN luminosities for a fixed M_halo and ⁠. We used a large halo mass bin to increase the statistics of our sample. We note that there is a correlation between central AGN luminosity and halo mass. To eliminate the effect of mass-dependent central AGN luminosity, we divide haloes according to the central AGN luminosity as follows. For each redshift, at each halo mass, we tag haloes as ‘high Lcen’ (‘low Lcen’) if their central AGN luminosities are above (below) the mean central AGN luminosity at that mass (Fig. 5). Then, the CLFs of satellite AGN are computed for the ‘high Lcen’ and ‘low Lcen’ haloes, respectively. The results are shown in Fig. 6. The solid lines show the distribution of satellite luminosities for ‘high Lcen’ sample and the dashed lines represent the distribution of satellite luminosities for the ‘low Lcen’ sample.

The two groups of curves are for two halo mass bins. The open and the filled black circles connected by the dashed and solid lines respectively represent the halo mass bin 11.0 ≤ log(M_halo) ≤ 12.0. Similarly, the blue open triangles and the open squares connected by the dashed and solid lines respectively represent the halo mass bin 12.0 ≤ log(M_halo) ≤ 13.0. The curves are normalized to the total number of haloes hosting the central AGN above and below the mean luminosity. We have shown the result for z= 3.0. Our results are similar for z= 1.0 and z= 5.0. The correlation between central and satellite AGN luminosity at a fixed halo mass will have important implications on the small-scale clustering strength. If there exists a strong large-scale feedback effect from the central AGN, it can possibly alter the local gas distribution around satellite AGN and suppress the growth of the satellite black holes. This in effect will decrease the number of luminous satellites in a halo with a higher central AGN luminosity and we would observe an anti-correlation between satellite number and central AGN luminosity.

In Fig. 6 we do not see any correlation between the number of satellites and the luminosity of the central AGN. It has been shown in Chatterjee et al. (2008) that in group-scale haloes (the most massive haloes in our simulation), feedback from the central AGN can extend up to a few hundred kiloparsec and thus can potentially affect the satellite distribution. However, in higher mass haloes, satellite AGN will also be residing in more massive subhaloes and so they will be relatively unaffected from the feedback effects of the central AGN. We note that in the self-regulatory growth paradigm, the AGN shuts down its own growth by blowing up the gas around it (as seen in Fig. 3) but the feedback effects do not affect the gas distribution around neighbouring black holes.

3.3 Peak versus instantaneous luminosity

where L_Bol is the peak bolometric luminosity. We overplot the instantaneous luminosity at z= 1.0 (red open circles). The best-fitting slope for the instantaneous luminosity is (0.95 ± 0.12). As we go to higher redshifts, the AGN have luminosities closer to their peak luminosities and hence we see a higher mean value and a tighter correlation (Fig. 5). The difference between the peak and the instantaneous luminosity is most prominent at higher mass haloes. We associate this suppression with feedback from AGN. In Paper I, we used the relation to show that black holes residing in higher mass haloes enter the feedback-dominated phase at an earlier time than black holes populating lower mass haloes. Thus, feedback effects will alter the correlation resulting in a wider distribution of host halo masses for a given AGN luminosity at lower redshift. This result is also in agreement with Lidz et al. (2006), who conclude that the peak luminosity is more correlated with halo mass than the instantaneous luminosity of AGN and hence a better indicator of clustering. However, we cannot measure the peak luminosity of the AGN in practice.

3.4 The mean occupation function of AGN

where 〈N_cen(M)〉 represents the mean occupation function of central AGN and 〈N_sat(M)〉 represents the mean occupation function of satellites (see discussions in Section 2.2 for the differences in terminology with Paper I). This formalism implicitly assumes that the AGN content of haloes of a given mass is statistically independent of the large-scale environments within which the haloes reside and thus the mean occupation 〈N(M)〉 depends only on halo mass (e.g. Bond et al. 1991; Lemson & Kauffmann 1999). 〈N(M)〉 is shown as a function of halo mass in Fig. 8. The black filled circles show the total occupation and the red open circles show the satellite occupation. The black solid, blue dot–dashed and the red dashed lines are the best-fitting models for the total, central and satellite occupation, respectively. The top, middle and bottom panels show redshifts 1.0, 3.0 and 5.0, respectively, while the left, middle and right columns denote luminosity thresholds 10³⁸, 10⁴⁰ and 10⁴² erg s⁻¹, respectively.

In this formalism there are four parameters for modelling the HOD: M_min, defining the halo mass where the occupation of central AGN of a given type (in the present case we choose AGN in terms of luminosity type) is 0.5; σ_logM, the characteristic transition width; M₁, the mass scale at which the mean number of satellites above a given luminosity threshold equals unity; and α, the power-law exponent for the satellite occupation. However, it has been observed that the mean occupation of the satellites drops off faster than a power law at lower halo mass (e.g. Kravtsov et al. 2004; Zheng et al. 2005; Conroy, Wechsler & Kravtsov 2006) and we also see this trend in our case. We use the parametrization of Tinker et al. (2005) to model this drop-off. The parameter M_cut is used to model the rolling off the power law. So we have a five-parameter HOD model. The best-fitting HOD parameters are shown in Table 2. The change in the luminosity threshold does affect the mean value of the power-law exponent but the values are consistent within 1σ. The power-law exponent also shows weak evolution (seen in Paper I) with redshift but is consistent with no evolution within statistical limits.

We compute the distribution of AGN with respect to the mean P(N|〈N〉). In Paper I, we showed that the distribution of satellite number (without any mass or luminosity cut) follows a Poisson distribution. We have verified that this is also the case for AGN in any luminosity threshold sample. We further performed a Kolmogorov–Smirnov (KS) test to check whether the satellite distribution follows a Poisson distribution in our luminosity-selected sample. The mean P-value that we obtain by performing the KS test over all redshifts and over all mass and luminosity bins is 0.92. This shows that the null hypothesis is strongly accepted and the distribution of the satellites is close to a Poisson distribution.

3.5 Luminosity versus mass

Our model of central occupation is slightly different from that of Paper I. In Paper I we do not impose any mass cut and hence the fraction of haloes (above the threshold mass) containing black holes is always unity (see equations 1 and 2 of Paper I). In this paper, the softening of the step function arises from the luminosity-based selection. We note that at lower luminosities, the duty cycle is extremely close to 1.0 at all halo mass bins and we converge to the model of Paper I. However, the drop-off from unity is evident in the 10⁴² erg s⁻¹ sample (right panel of Fig. 8). In Fig. 9, we show the differences in the mean occupation function between a mass-selected and a luminosity-selected sample. The black solid lines show the mean occupation function of AGN with bolometric luminosities greater than 10⁴⁰ erg s⁻¹ at z= 1.0 (shown in the top middle panel of Fig. 8). The triangles represent the mean occupation of central AGN and the open circles represent the mean occupation of satellite AGN. We now use the best-fitting relation between black hole mass and AGN luminosity at z= 1.0 (Chatterjee et al. 2008) to obtain the equivalent mass corresponding to our bolometric luminosity of 10⁴⁰ erg s⁻¹. The red dashed lines represent the mean occupation functions for the mass-selected sample with triangles representing central black holes and open circles representing satellite black holes.

We see a clear difference between the two populations at lower halo masses. For the central AGN (triangles), the solid and the dashed lines converge at a mass scale of 10^11.5 M_⊙. Also, the mass-selected occupation function falls off very steeply below this mass scale. This difference arises from lower mass black holes residing in lower mass haloes with high Eddington ratios. Also, the correspondence between black hole mass and host halo mass is tighter than luminosity and hence we see this sharp cut-off in the occupation function for the mass-selected sample. Recently, Gallo et al. (2010) found evidence of such low-mass high accretion rate black holes in local early-type galaxies in the AGN Multi-Wavelength Survey of Early Type Galaxies (AMUSE)-Virgo survey. For the satellite population (open circles) the lack of convergence in the occupation function between the two populations is prominent even at higher halo masses. The minimum halo mass for hosting satellite black holes (based on mass) is much higher than minimum halo mass for hosting satellite AGN with equivalent luminosity. We thus see that the HOD properties will be significantly different between a mass-selected sample and a luminosity-selected sample. This is because a relatively small halo can host a bright AGN if there is high-density gas. On the other hand, it is difficult to get a massive black hole in a small halo, since that would require significant amounts of dense gas over a long period of time.

3.6 Radial profiles

Fig. 10 shows the radial distribution of the number density of satellite AGN (luminosity greater than 10⁴⁰ erg s⁻¹) within host haloes of masses 11.0 ≤ log(M_halo) ≤ 13.0. The red circles, blue triangles and the green squares show the profiles at z= 1.0, 3.0 and 5.0, respectively. It has been shown in Paper I that the radial distribution of black holes follows a power law at all redshifts (equation 6 of Paper I). Our findings are similar to Paper I. The black dashed line is the best-fitting power law to the average profile over all redshifts (since we do not see any significant evolution with redshift). The best-fitting values for β (power-law index of AGN: averaged over all redshifts) and log(n₀) (normalization: averaged over all redshifts) are −2.33 ± 0.08 and −0.67 ± 0.05 for the sample shown in Fig. 10. The power-law index β does not show any significant dependence on halo mass or AGN luminosity either. For comparison, we also show the profile of dark matter and galaxies within AGN host haloes. The solid and the dotted lines show the average profiles (averaged over all redshifts) for dark matter and satellite galaxies, respectively. The minimum stellar mass that we used for obtaining the galaxy profile is 10⁸ M_⊙ (roughly 100 times the stellar mass resolution). The profiles do not show any variation if we change the threshold stellar mass for selecting galaxies. The current choice of stellar mass is an optimization between resolution elements and statistics. The profiles are normalized to the mean number density of the corresponding species (satellite AGN/dark matter/satellite galaxies) within R₂₀₀ (e.g. Nagai & Kravtsov 2005).

We also fit the Navarro, Frenk & White (1997) profile (NFW profile) to the AGN radial distribution. We fit the profiles for different concentration parameters and calculate the corresponding P-values. At all redshifts our data strongly disfavours the null hypothesis and the NFW profile is ruled out at 3σ. We note that the AGN are more centrally concentrated than dark matter and galaxies. The reason that we see an enhanced population of AGN at the centre of the halo compared to galaxies is because of the merging process of black holes. When two galaxies merge there exists a time lag between merging galaxies and the merging of AGN that reside in them. This time lag is due to the time it takes for the satellite AGN to fall into the halo centre, where it can merge with the central AGN. After the two galaxies merge, the time that it takes for the AGN to merge can then be further affected by the gas content of these galaxies. Lin & Mohr (2007) measure the radial profile of radio sources in clusters and show that it is consistent with an NFW profile with a concentration of 25. Martini, Mulchaey & Kelson (2007) study the radial distribution of X-ray-selected AGN in clusters and find that AGN with X-ray luminosities above 10⁴² erg s⁻¹ show stronger central concentration than cluster host galaxies. However, a larger sample with L_X≥ 10⁴¹ erg s⁻¹ (closer to our sample of AGN) does not show any stronger evidence of central concentration, different from what we observe in simulation.

In Fig. 11 we show the relation between stellar mass and AGN luminosity at z= 1.0. The stellar mass is computed within a spherical region of radius 30 kpc surrounding the black hole. We selected different spatial scales to compute the stellar mass. The values converge between radii of 25 and 30 kpc and hence we chose 30 kpc to be the relevant radius for computing the stellar mass content in AGN host galaxies. The black points represent central AGN and the red open circles show the satellite distribution. In the case of central AGN, the stellar mass for more luminous AGN is generally higher. AGN will have higher luminosity when the accretion rate is high, which is more likely when the gas density is high. Also, when the gas density is high, the cooling rate will be higher and hence there will be more stars. So, in general, higher density will correspond to high accretion rates of AGN and higher stellar mass. For satellite AGN, this trend is very weak. But even for satellite AGN if we consider a luminosity-selected sample we will systematically choose AGN with higher stellar mass. Also, in Fig. 3 we see that AGN with higher stellar mass content tend to reside in the central regions of the haloes. This can potentially increase the higher concentration of AGN seen in Fig. 10 for a luminosity-selected sample. To further check this we selected AGN based on their stellar mass content and computed their average radial profile. We find a steeper mean slope (β∼−2.7) for the stellar mass-selected sample compared to the luminosity-selected sample. This again is an indication of enhanced central concentration of AGN with higher stellar mass content in their host galaxies.

4 SUMMARY AND DISCUSSIONS

In this paper, we study the relation between AGN and dark matter haloes using a cosmological SPH simulation that incorporates black hole and galaxy formation physics. Specifically, we have investigated the effect of environmental factors (e.g. feedback, local gas density and host galaxy mass) on the accretion rates of AGN within dark matter haloes and how it affects the occupation distribution of AGN. We examine a number of simplifying assumptions in the HOD modelling (e.g. dependence of AGN luminosity on halo mass and dependence of AGN accretion on black hole mass) of AGN and provide the necessary tools to model AGN clustering data. We have characterized the HOD of faint AGN (the sample probed by our simulation box) which can be generalized to incorporate the bright end of the luminosity function.

We compute the CLFs of AGN and separate the contributions from the central and satellite AGN. Our key findings are as follows. (1) The central AGN luminosities follow a lognormal distribution similar to the assumption in semi-analytic studies. (2) The mean of the CLF for a given halo mass shows a strong dependence on redshift with higher luminosities at higher redshifts. (3) The scatter in central AGN luminosity is large at low redshift, but decreases with increasing redshift. This implies that the dependence of AGN clustering on luminosity is weak at low redshift, but it can become stronger at high redshift. We analyse the light curves of individual AGN and show that there exists a tighter correlation between halo mass and peak luminosity rather than instantaneous luminosity. We present the joint distribution of satellite occupation as a function of halo mass and central AGN luminosity. We do not see any significant correlation between the satellite number and the luminosity of the central AGN. We also show that the mean occupation function of the central AGN resembles a softened step function while the satellite population follows a power law with an exponential roll-off at lower mass, similar to what has been observed with the galaxy HOD. We show that low-mass black holes with high Eddington ratios residing in low-mass haloes make the luminosity-based HOD significantly different from the black hole mass-based HOD.

We will now compare our simulation results with semi-analytic models that have been widely used for clustering analysis (e.g. Martini & Weinberg 2001; Shankar et al. 2010a). We compare our results with the model described in Shankar et al. (2010a), who computed the mass function of black holes at z= 3.0 using the halo mass function and a redshift-dependent relation between halo mass and black hole mass. They parametrize the redshift dependence in the normalization of the black hole mass–halo mass relation while the slope is kept fixed at 1.52. In our simulation, we observe the average slope to be (1.04 ± 0.06) at z= 3.0. There is a large scatter from the mean slope and the 1σ scatter is 0.44. We note that this relation in our simulation is dependent on the AGN model and the ratio of the threshold halo mass that can host black holes to the mass of the seed black holes. Shankar et al. (2010a) assume the black holes to be accreting with a constant Eddington fraction. Although the mean χ in our simulation is roughly independent of halo mass as seen in Fig. 1, they do show an evolution with redshift. This has also been noted by Shankar et al. (2010a), who suggest that the assumption of constant Eddington fraction might break down at the faint end. The other important parameter is the scatter in the relation. We find a lognormal distribution as assumed in semi-analytic studies. However, the scatter strongly depends on redshift with increasing scatter at low redshifts (shown in Fig. 5).

Observationally, the most relevant physical quantity describing an AGN is luminosity. Recently, Miyaji et al. (2011) and Starikova et al. (2011) have used X-ray-selected AGN to constrain the HOD empirically. We propose to compare our HOD model with observational samples in a follow-up paper. However, an alternative approach would be selecting black holes based on their mass and measure clustering statistics based on black hole mass. Although we do not have reliable observational measurements of black hole masses, we have also provided the framework for predicting clustering properties based on the black hole mass function in Paper I. These two complimentary approaches would impose even tighter constraints on theoretical models of AGN growth and feedback. Our analysis provides the first step towards comparing semi-analytic and simulation results on AGN clustering and assessing some of the simplifying assumptions in the present interpretations of AGN clustering observations and constraining physical parameters. Because of the small simulation box, our current study is limited to low-luminosity AGN, and the results on satellite AGN suffer from small number statistics. A larger simulation box with larger statistical samples spanning both the faint and the bright ends of the luminosity function can provide tighter constraints on AGN HOD. Also, the mass of the seed black hole is dependent on the resolution of the simulation, imposing an artificial mass cut in the simulation. Future work on high-resolution simulation and more accurate modelling of accretion and feedback (e.g. Booth & Schaye 2009) will also be needed to understand the full implication of gas physics and black hole accretion on clustering studies. The HOD formalism has been successfully incorporated in galaxy evolution. We hope that our work will have the same impact for studying the co-evolution of AGN and their hosts and shed light on the AGN–galaxy connection.

We would like to thank Douglas Rudd, Laurie Shaw, Yue Shen and Frank van den Bosch for some useful discussions which helped in the analysis and physical interpretation of some results. We thank David Weinberg for suggesting the analysis for looking at the correlation between central and satellite AGN luminosities. We also thank the referee for the comments which helped in improving the draft. JR was supported by the Yale College Dean’s Science Research Fellowship through summer 2010. ZZ gratefully acknowledges support from Yale Center for Astronomy and Astrophysics through a YCAA fellowship. DN was supported in part by the NSF grant AST-1009811, by NASA ATP grant NNX11AE07G and by Yale University. The work in CMU was supported by the National Science Foundation, NSF Petapps, OCI-0749212 and NSF AST-1009781. The simulations were carried out at the NSF Teragrid Pittsburgh Supercomputing Center (PSC).

REFERENCES