A refined sub-grid model for black hole accretion and AGN feedback in large cosmological simulations

Abstract

In large-scale cosmological hydrodynamic simulations simplified sub-grid models for gas accretion on to black holes and AGN feedback are commonly used. Such models typically depend on various free parameters, which are not well constrained. We present a new advanced model containing a more detailed description of AGN feedback, where those parameters reflect the results of recent observations. The model takes the dependence of these parameters on the black hole properties into account and describes a continuous transition between the feedback processes acting in the so-called radio-mode and quasar-mode. In addition, we implement a more detailed description of the accretion of gas on to black holes by distinguishing between hot and cold gas accretion. Our new implementations prevent black holes from gaining too much mass, particularly at low redshifts, so that our simulations are successful in reproducing the observed present-day black hole mass function. Our new model also suppresses star formation in massive galaxies slightly more efficiently than many state-of-the-art models. Therefore, the simulations that include our new implementations produce a more realistic population of quiescent and star-forming galaxies compared to recent observations, even if some discrepancies remain. In addition, the baryon conversion efficiencies in our simulation are – except for the high-mass end – consistent with observations presented in the literature over the mass range resolved by our simulations. Finally, we discuss the significant impact of the feedback model on the low-luminous end of the AGN luminosity function.

black hole physics, methods: numerical, galaxies: active, galaxies: evolution, galaxies: nuclei, quasars: supermassive black holes

1 INTRODUCTION

Black holes (BHs) play an essential role in the formation and evolution of galaxies. They can even influence galaxy clusters and the intra cluster medium (ICM). However, observations of active galactic nuclei (AGN) indicate that gas accretion on to BHs and AGN feedback are complex processes, which are not yet fully understood (e.g. Merloni & Heinz 2007; McNamara, Rohanizadegan & Nulsen 2011; Ma, McNamara & Nulsen 2013). There is evidence for two distinct phases of AGN activity and feedback: the radio-mode and the quasar-mode. The radio-mode is characterized by large radio jets generating hot X-ray cavities (Russell et al. 2013; Mezcua & Prieto 2014), whereas in the quasar-mode the emission is dominated by the accretion disc, which is visible as the so-called blue bump in the spectrum of quasars and Seyfert galaxies (e.g. Elvis et al. 1994; Prieto et al. 2010).

Churazov et al. (2005) characterized this distinction in a theoretical model by describing AGN feedback with two components: radiation and mechanical outflow. In their model, the amount of energy associated with each component depends on the Eddington ratio |$f_{\mathrm{Edd}}=\dot{M_\bullet }/\dot{M}_{\mathrm{Edd}}$|⁠. When a BH accretes with the Eddington accretion rate |$\dot{M}_{\mathrm{Edd}}$|⁠, gas cooling and AGN feedback are in equilibrium. Churazov et al. (2005) also took advection-dominated accretion flows (ADAFs) into account, although a jet contribution can successfully replace an ADAF (Falcke, Körding & Markoff 2004; Fernández-Ontiveros et al. 2011).

To constrain this model and to really understand the origin of different types of AGN and how they influence their environment, large cosmological simulations play a key role. They have two major advantages: first, they provide a statistically large sample of BHs. This allows us to compare the simulations to the newest and currently most complete observations of the M_•–M_* relation (e.g. McConnell & Ma 2013) or BH mass functions (e.g. Marconi et al. 2004; Shankar et al. 2004; Shankar, Weinberg & Miralda-Escudé 2009) and stellar mass functions (e.g. Bernardi et al. 2013; Muzzin et al. 2013), in particular the very massive end. Secondly, having large enough cosmological boxes where also massive galaxy clusters form, allows us to probe the influence of BHs across all scales of cosmic environment.

There already exist a number of studies discussing large cosmological simulations that include BHs (e.g. Di Matteo, Springel & Hernquist 2005; Robertson et al. 2006; Di Matteo et al. 2008; Booth & Schaye 2009; Degraf, Di Matteo & Springel 2011; Teyssier et al. 2011; Rosas-Guevara et al. 2013; Hirschmann et al. 2014; Khandai et al. 2014; Vogelsberger et al. 2014b; Schaye et al. 2015). Those simulations mostly use the BH model implemented by Springel, Di Matteo & Hernquist (2005) or are based on it. In these models – in contrast to some more simplified BH models (e.g. Battaglia et al. 2010) – BHs are typically described as sink particles which have fundamental properties like mass and accretion rate, which can be linked directly to observables. Hence, we can study BH growth and the co-evolution between BHs and their host galaxies to constrain and improve the parametrization of the underlying model. In the model from Springel et al. (2005), the gas accretion on to BHs is calculated according to the Bondi formula (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952), multiplied by a so-called boost factor α. This factor was introduced to account for the limited resolution in simulations leading to smaller densities and larger temperatures near the BH (Booth & Schaye 2009). To estimate the AGN feedback, a constant value for the radiative efficiency is typically used (Shakura & Sunyaev 1973).

For low resolutions this model works reasonably well. However, to study not only the origin of the observed fundamental relations between BHs and their host galaxies (Tremaine et al. 2002; Häring & Rix 2004; McConnell & Ma 2013), but also the impact of gas accretion and AGN feedback on the morphology of the galaxy, simulations with higher resolution are needed. Until now, this was only studied in simulations of isolated galaxies and mergers of galaxies (e.g. by Hopkins et al. 2008; Debuhr, Quataert & Ma 2011; Capelo et al. 2015; Van Wassenhove et al. 2014) as well as in cosmological zoom simulations (e.g. by Anglés-Alcázar, Özel & Davé 2013; Dubois et al. 2013; Choi et al. 2014b; Marinacci, Pakmor & Springel 2014). To reproduce both statistical BH and galaxy properties within a fully cosmological context and across various environments in a statistically relevant sample size, large cosmological boxes with high resolution are needed. This is still a challenge, but thanks to increasing computational power it now becomes feasible. However, despite of this success, new challenges arise as simulations typically overestimate the high-mass end of the BH and stellar mass function (e.g. Genel et al. 2014; Hirschmann et al. 2014; Khandai et al. 2014; Sijacki et al. 2014; Vogelsberger et al. 2014a). Therefore, a more detailed BH model is necessary.

In this work, we extend the model by Springel et al. (2005) by improving the treatment of the two modes of AGN feedback: radiation and mechanical outflows. Following theoretical predictions (White & Frenk 1991; Narayan & Yi 1995; Churazov et al. 2005) as well as recent observational results (Davis & Laor 2011; Chelouche 2013; Russell et al. 2013) gives us estimates for the corresponding two efficiencies depending on the BH mass and the accretion rate, which outreaches the simplified BH model commonly used in simulations.

Following Sijacki et al. (2007), a steep transition between radio-mode and quasar-mode is often used in current simulations (e.g. Fabjan et al. 2010; Hirschmann et al. 2014). This is only a rough approximation to the smooth transition which is observed and also theoretically expected. Adopting the model by Churazov et al. (2005) – which was already constrained by observations, e.g. Russell et al. (2013) – allows us to get a smooth transition between the two modes. This was used by Hirschmann et al. (2014) to calculate AGN luminosities, but it was never implemented into simulations. Such modifications were also suggested by a recent paper of Sijacki et al. (2014), who studied the AGN luminosity function within a cosmological simulation using a constant radiative efficiency. They concluded that in the radio-mode radiative efficiencies might depend on the accretion rate and on average should be lower than the value 0.1 used in the original BH model from Springel et al. (2005). Furthermore, Davis & Laor (2011) and Chelouche (2013) found that the radiative efficiency not only correlates with the accretion rate, but also with the BH mass.

Another deficiency in current implementations of BHs in cosmological simulations is that the (original) Bondi model predicts far too low accretion rates during the quasar-mode so that BHs do not reach the observed masses for a given bulge mass. Therefore, a so-called boost factor is commonly used to artificially raise the accretion rates. This results in realistic accretion rates for the accretion of cold gas. However, it has the disadvantage that it also raises the accretion rate when the hot gas content is large enough to fulfill the assumptions of the Bondi model, namely when the gas is distributed in an isotropic sphere. This typically is the case in old quiescent galaxies. Consequently, BHs become too massive at low redshifts. Hence, accretion rates have to be lower in the radio-mode (Li, Ostriker & Sunyaev 2013).

Indeed, several studies adapt the BH model for higher resolution simulations by using a boost factor which depends on the resolution (Choi et al. 2012; Choi et al. 2014a), density (Booth & Schaye 2009), pressure (Vogelsberger et al. 2013) or angular momentum (Rosas-Guevara et al. 2013), although none of them contains a direct distinction between the accretion of cold and hot gas, even if the existence of such two distinct accretion modes has been shown by observations (e.g. Hlavacek-Larrondo et al. 2013) and predicted by high-resolution simulations of BH accretion on sub-kpc scales (Gaspari, Ruszkowski & Oh 2013; Bourne, Nayakshin & Hobbs 2014) as well as semi-analytical models (e.g. Somerville et al. 2008; Fanidakis et al. 2011, 2013; Hirschmann et al. 2012). A distinction between accretion of cold and hot gas based on the multiphase model from Springel & Hernquist (2003) was implemented in the simulations from Pelupessy, Di Matteo & Ciardi (2007). In their study, the molecular gas of the star-forming particles was evaluated from a multiphase model, in which the accretion of this cold gas was evaluated separately without any boost factor, assuming the corresponding temperature as fixed in the underlying multiphase model.

A BH mainly grows in the quasar-mode, where cold gas forms an accretion disc around the BH which leads to higher accretion rates. During that period, BHs grow until the AGN feedback and gas cooling are in equilibrium. At that point, they reach the M_•–σ relation (Churazov et al. 2005) and thus, the M_•–M_* relation. Consequently, the accretion rate drops until the BH crosses the threshold towards the radio-mode. As reviewed by several authors (e.g. Heckman & Best 2014; Yuan & Narayan 2014), the accretion in the radio-mode, sometimes also called jet-mode, can be described with ADAFs containing hot gas (Yuan et al. 2009). Alternatively, the accretion of hot adiabatic gas can be described with the Bondi model (Gaspari et al. 2013). Therefore, we distinguish between hot and cold gas and estimate the accretion rate separately for both gas phases. This allows us to use different boost factors for hot and cold gas and thus, to account for both observed accretion modes.

The outline of this paper is as follows: in Section 2 we describe our BH model. The set-up of the cosmological simulations is presented in Section 3. In Section 4, adopting different models for BH accretion and AGN feedback, we show the results for our simulations, in particular the evolution of the BH mass, the stellar mass and the star formation rate. In Section 5, we discuss the radiative efficiency in the radio-mode and its influence on to the AGN luminosity functions. Furthermore, we compare our results with other cosmological simulations. Finally, in Section 6, we summarize our main results.

2 THEORETICAL MODEL

2.1 BH accretion

The Bondi model is commonly used in simulations to estimate the BH accretion rate. The Bondi accretion rate |$\dot{M}_\mathrm{B}$| (Bondi 1952; Shima et al. 1985) is given by

\begin{equation} \dot{M}_\mathrm{B} = \frac{4 \pi G^2 M_\bullet ^2 \rho _{\infty }}{(v^2 + c_{\mathrm{s}}^2)^{3/2}}, \end{equation}

(1)

where M_• is the BH mass, ρ is the density, c_s is the sound speed of the accreted gas and v is the velocity of the gas relative to that of the BH. Since Bondi (1952) assumed an isotropic and isothermal sphere of gas for his estimation, it is not straightforward to adopt this Bondi accretion model for hydrodynamic, cosmological simulations aiming to follow a self-consistent accretion history of BHs. For the implementations based on Springel et al. (2005), the accretion rate of the BH is estimated by

\begin{equation} \dot{M}_{\mathrm{B}} = \frac{4\pi \alpha G^2 M_{\bullet }^2 \left\langle \rho \right\rangle }{(\left\langle c_\mathrm{s} \right\rangle ^2 + \left\langle v \right\rangle ^2)^{3/2}}, \end{equation}

(2)

where 〈ρ〉, 〈v〉 and 〈c_s〉 are computed using kernel-weighted smoothed particle hydrodynamics (SPH) estimations. Due to limited numerical resolution in such simulations, the original equation (1) is multiplied by a boost factor α, which in Springel et al. (2005) is set to a value of α = 100. Note that the SPH estimates also depend on the type of SPH kernel and the number of neighbours. To make this estimation less sensitive to the actual structure of the multiphase media in the vicinity of the BH and therefore the algorithm less dependent on resolution and on the actual choice of numerical parameters for the kernel-weighted interpolation, Choi et al. (2012) suggested to use a different way of building the averages:

\begin{equation} \dot{M}_{\mathrm{B}} = \left\langle \frac{4\pi \alpha G^2 M_{\bullet }^2 \rho }{(c_\mathrm{s}^2 + v^2)^{3/2}}\right\rangle . \end{equation}

(3)

Still, choosing the correct value for the boost factor α is not trivial. Since due to the limited resolution the density in the not resolved vicinity of BHs is large, it will be underestimated and – in turn – the temperature (and thus the sound speed) will be overestimated. Following this argument, Booth & Schaye (2009) parametrize α, which is chosen to be α = 1 as long as the density is below the critical value where one can assume the gas to be in the hot phase. For larger densities, when gas is accreted mainly in a cold phase, α increases with density. Alternatively, Vogelsberger et al. (2013) have presented a recipe for modelling α based on the equilibrium between cooling losses and AGN feedback. However, both models do not directly account for the different accretion modes of hot and cold gas phase, where cold gas usually is accreted in turbulent streams, whereas hot gas indeed can be assumed to be isotropic and isothermal.

In our model, we use a sixth-order Wendland kernel (Dehnen & Aly 2012) with 295 neighbours, building the mean values according to equation (2) and directly distinguishing between the accretion of hot and cold gas. In this way, we can safely use the original estimate of building the averages, which has the advantage to be more sensitive to density structures close to the BH. In general, we assume that hot gas has temperatures above T ≈ 10⁶ K, whereas cold gas has temperatures below T ≈ 10⁵ K (Gaspari et al. 2013). Since we do not account for a third warm phase, we choose T = 5 × 10⁵ K as threshold between hot and cold gas. In contrast to Pelupessy et al. (2007), who use the molecular fraction of the gas for star-forming particles from the multiphase model (Springel & Hernquist 2003) to account for cold gas accretion, we also assign gas with a temperature below our threshold in addition to the star-forming gas to the cold phase. For both gas phases, the accretion rate is calculated separately according to equation (2), but with different values for α according to the result by Gaspari et al. (2013), who argue that due to turbulence the assumptions of the Bondi model are not fulfilled for the cold gas. When they include cooling and turbulence in their simulation, they find an accretion rate which is around 100 times larger than the Bondi accretion rate. Interestingly, this is the same value which is used as boost factor α in the original model from Springel et al. (2005). But for adiabatic accretion, the difference, Gaspari et al. (2013) find, is about one order of magnitude smaller. Hence, we use α = 10 for hot gas and α = 100 for cold gas.

Furthermore, the BH accretion rate |$\dot{M_\bullet }$| is limited to the Eddington accretion rate

\begin{equation} \dot{M}_{\mathrm{Edd}} = \frac{4 \pi G M_{\bullet } m_\mathrm{p}}{\eta _{\mathrm{Edd}} \sigma _\mathrm{T} c}, \end{equation}

(4)

where m_p is the proton mass, σ_T the Thompson scattering cross-section and η_Edd the feedback efficiency if the BH would accrete with |$\dot{M}_{\mathrm{Edd}}$|⁠. Then the accretion rate is given by

\begin{equation} \dot{M_\bullet } = \mathrm{min}(\dot{M}_\mathrm{B, hot} + \dot{M}_\mathrm{B, cold}, \dot{M}_{\mathrm{Edd}}). \end{equation}

(5)

The distinction between hot and cold gas accretion leads to a faster BH growth in the quasar-mode, because when calculating the mean value of the sound speed 〈c_s〉 and the gas velocity 〈v〉 only for cold gas, the accretion rate estimated with equation (2) is higher than calculating the mean values of both cold and hot gas together. This solves the well-known problem of too low gas accretion, which was addressed in other simulations by increasing the maximum accretion rate to a few times |$\dot{M}_{\mathrm{Edd}}$| (e.g. Di Matteo et al. 2012), which is not needed in our simulations.

2.2 AGN feedback

In the commonly used BH model by Springel et al. (2005), the feedback energy per unit time is calculated as

\begin{equation} \dot{E} = \epsilon _\mathrm{f} \epsilon _\mathrm{r} \dot{M_\bullet } c^2, \end{equation}

(6)

where ϵ_f is the efficiency with which the energy radiated from the BH is coupled to the ISM (Springel et al. 2005; Booth & Schaye 2009) and ϵ_r is the radiative efficiency.

The original model as used in Hirschmann et al. (2014) is simplified, since it uses a constant radiative efficiency and thus does not allow for a smooth transition between quasar- and radio-mode. Furthermore, it neglects mechanical feedback, which was already implemented in other simulations as AGN-driven winds (i.e. Choi et al. 2014b). To account for both mechanical and radiative feedback, we adopt a new feedback scheme based on Churazov et al. (2005). In this study, they propose that AGN feedback can be split up into two components.

• Outflow. The outflow component is a mechanical feedback, which dominates at accretion rates below |$\sim\! 0.01\dot{M}_{\mathrm{Edd}}$| and diminishes at accretion rates above |$\sim\! 0.1\dot{M}_{\mathrm{Edd}}$|⁠. The corresponding gas heating power is given by

  \begin{equation} P_\mathrm{o} = \epsilon _\mathrm{o} \dot{M_\bullet } c^2, \end{equation}

  (7)

  where ϵ_o is the outflow efficiency.
• Radiation. The radiative component dominates near the Eddington limit (f_Edd > 0.1) and has the luminosity

  \begin{equation} L=\epsilon _\mathrm{r} \dot{M_\bullet } c^2. \end{equation}

  (8)

We implement both radiative and mechanical AGN feedback as thermal feedback due to the inability to resolve the sub-kpc scales, where the jets provide the mechanical feedback. The feedback energy per unit time in this model is then the sum of P_o and the fraction ϵ_f of the luminosity:

\begin{equation} \dot{E} = (\epsilon _\mathrm{o} + \epsilon _\mathrm{f} \epsilon _\mathrm{r}) \dot{M_\bullet } c^2. \end{equation}

(9)

The effect of accreted matter can be split into outflow and radiation components:

\begin{equation} \frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}} = \frac{P_\mathrm{o}}{L_{\mathrm{Edd}}} + \frac{L}{L_{\mathrm{Edd}}}, \end{equation}

(10)

where the Eddington accretion rate

\begin{equation} \dot{M}_{\mathrm{Edd}} = \frac{L_{\mathrm{Edd}}}{\eta _{\mathrm{Edd}} c^2} \end{equation}

(11)

depends on the total efficiency

\begin{equation} \eta := \epsilon _\mathrm{o} + \epsilon _\mathrm{r}. \end{equation}

(12)

This model is shown as solid lines (blue corresponds to mechanical outflow and red to radiation) in Fig. 1, which were adopted from Churazov et al. (2005). For the outflow-dominated regime they assume

\begin{equation} \frac{L}{L_{\mathrm{Edd}}} = 10 \times \ \left(\frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}}\right)^2 \end{equation}

(13)

as a lower limit for the radiation, which is a consequence of ADAFs (Narayan & Yi 1995). In the radiation-dominated regime, the outflow decreases with the Eddington ratio:

\begin{equation} \frac{P_\mathrm{o}}{L_{\mathrm{Edd}}} = 10^{-4} \times \left(\frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}}\right)^{-1.8431}. \end{equation}

(14)

This guarantees that the minimum value for the outflow efficiency is ϵ_o = 10⁻⁵, which was calculated by Churazov et al. (2005) assuming that gas cooling and AGN feedback balance each other at the Eddington limit. We choose |$\frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}} = 0.05$| as the threshold between radio- and quasar-mode. The value for the outflow at |$\frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}} = 1$| follows the calculations of Churazov et al. (2005), who find ϵ_o ≈ 10⁻⁵ for BHs accreting with the Eddington accretion rate.

Figure 1.

The lines show the predictions by Churazov et al. (2005, C05) for the power of the radiation (red line), the mechanical outflow (blue line) and the sum of both (black dashed line). Observations of jet powers (blue error bars and edges) and luminosities (red error bars and edges) constrain the difference between both components. This figure includes two different observations: the big stars and squares show recent observations by Mezcua & Prieto (2014, MP14) and the data with blue and black error bars are observations by Russell et al. (2013, R13). Black triangles mark upper limits. Furthermore, the BH masses are indicated by the colours of the symbols. Since the masses used by R13 are based on K-band magnitudes, which are known to be inaccurate, we used the dynamical masses by McConnell & Ma (2013) for the sources included in both samples.

The feedback model of Churazov et al. (2005) was recently confirmed by observations (see also Russell et al. 2013) measuring luminosities and cavity powers of a large sample of unresolved nuclear X-ray sources. Most of the selected brightest cluster galaxies have large X-ray cavities. The data from Russell et al. (2013) show a large scattering of the luminosities in the radio regime illustrated by round filled circles with black error bars in Fig. 1, implying that a secondary quantity influences the luminosity. A few data points are below the theoretical lower limit, albeit the uncertainties in the observations are relatively high. Uncertainties can occur, for example, when measuring the cavity volume due to projection effects. In Fig. 1, the BH masses are colour-coded as indicated by the colourbar. The masses from Russell et al. (2013) are based on K-band magnitudes, which are known to be problematic. Therefore, we use the dynamical masses from McConnell & Ma (2013) for the sources included in both samples. Nearly all BHs that lie below the prediction are very massive (>10⁹ M_⊙). For lower masses, the observations are in better agreement with the predictions. We will discuss the uncertainties in Section 5.2 in more detail.

Recently, Mezcua & Prieto (2014) presented measurements of luminosities of a much smaller sample of AGN, but with sufficiently larger angular resolution and sensitivity. Their estimations for L_bol are more reliable than those presented in Russell et al. (2013), because they measure L_bol after integrating the radio to X-ray spectral energy distribution. Furthermore they explicitly provide values for X-ray cavity powers. For CenA, M87 and NGC 1052, they used X-ray cavities of maser emission from the literature (Prieto et al. 2010; Fernández-Ontiveros et al. 2012; Russell et al. 2013). All other values were estimated using the correlation between core radio luminosity at 5 GHz and P_o of Merloni & Heinz (2007). The data from Mezcua & Prieto (2014) is also included in Fig. 1, where the filled stars represent the luminosities and the squares the cavity powers. Since equation (13) is a lower limit, their luminosities are in very good agreement with the predictions. The cavity powers do not always match the blue line, but as described by Mezcua & Prieto (2014), they are expected to be lower limits, because the estimations of P_o do not take into account the energy which is used to compress the gas when the jet advances the ISM/ICM.

In simulations, the theoretical and observational results shown in Fig. 1 can be used to calculate the efficiencies ϵ_o and ϵ_r. To estimate the radiative and outflow efficiencies, we first have to assume a value for the total efficiency η and then use the predictions from Churazov et al. (2005) to separate the AGN feedback into radiation and mechanical outflow. In theoretical studies, the total efficiency is often assumed to be 0.1 (e.g. Churazov et al. 2005), however, observations of Davis & Laor (2011) and Chelouche (2013) suggest a mass dependence of this parameter. In the model from Churazov et al. (2005), both ϵ_o and ϵ_r depend on the accretion rate and the total efficiency. For |$\dot{M_\bullet }/\dot{M}_{\mathrm{Edd}} < 0.05$|⁠, the lower limit for ϵ_r can be calculated with equations (8) and (13), i.e.

\begin{equation} \epsilon _\mathrm{r, min} = 10 \eta \frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}} . \end{equation}

(15)

Since this is only a lower limit, all solutions between ϵ_{r, min} and ϵ_{r, max} = η are possible. Therefore, we introduce the slope β, which is in the range between 0 and 1, to get a general expression for ϵ_r:

\begin{equation} \epsilon _\mathrm{r} = A \eta \left(\frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}}\right)^{\beta }, \end{equation}

(16)

where A = 10⁻⁴ × 0.05^{−2.8431 − β}. The outflow efficiency is calculated with equations (16) and (12).

For |$\dot{M}/\dot{M}_{\mathrm{Edd}} > 0.05,$| the radiation dominates. The origin of the blue line in Fig. 1 in this regime is the analytical calculation by Churazov et al. (2005), which is based on the equilibrium between gas cooling and heating of gas due to AGN feedback. Hence, it is not only a lower limit and it is not necessary to introduce a slope as in the radio regime. In that respect from equations (7) and (14) follows

\begin{equation} \epsilon _\mathrm{o} = 10^{-4} \eta \left(\frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}}\right)^{-2.8431} \end{equation}

(17)

and thus ϵ_r = η − ϵ_o. This is shown in Fig. 2 for different BH masses. The filled circles and diamonds in Fig. 2 are the observations from Davis & Laor (2011) and Chelouche (2013) illustrating that they are not consistent with the model for η = 0.1 (green lines). Therefore, we account for the observed spin of BHs by following the observations of Davis & Laor (2011) for quasars and of Chelouche (2013) for Seyfert 1 AGN, who both find a correlation between the radiative efficiency and the BH mass. Hence, we use the relation found by Davis & Laor (2011) to estimate the total efficiency at the Eddington limit, which is approximately the same as the radiative efficiency at the Eddington limit:

\begin{equation} \eta _{\mathrm{Edd}}(M_{\bullet }) \approx \epsilon _\mathrm{r, Edd}(M_{\bullet }) = 0.089 \left(\frac{M_{\bullet }}{10^8 \,\mathrm{M}_{{\odot }}}\right)^{0.52}. \end{equation}

(18)

We limit η_Edd(M_•) by the value 0.42, which is the theoretical maximum efficiency of a rotating BH. To calculate the outflow efficiency, the constant value of η = 0.1 is used as it is currently difficult to estimate outflow efficiencies with observations (see Section 5.2 for further discussion). Equations (12), (16) and (17) then lead to the following set of equations:

\begin{eqnarray} \epsilon _{\mathrm{r}} =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}A \eta _{\mathrm{Edd}}(M_{\bullet }) \left(\frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}}\right)^{\beta },\quad {\rm if }\frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}} < 0.05,\\ \eta _{\mathrm{Edd}}(M_{\bullet }) - 10^{-4} \eta _{\mathrm{Edd}}(M_{\bullet }) \left(\frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}}\right)^{-2.8431}, \\ {\rm otherwise} \end{array}\right. \end{eqnarray}

(19)

and

\begin{eqnarray} \epsilon _{\mathrm{o}} =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}0.1 - A 0.1\left(\frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}}\right)^{\beta },&\quad {\rm if }\;\frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}} < 0.05, \\ 10^{-5} \left(\frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}}\right)^{-2.8431},&\quad {\rm otherwise.} \end{array}\right. \end{eqnarray}

(20)

In our simulations both radiative and mechanical feedback are implemented as thermal feedback, since we do not resolve jets.

Figure 2.

Our NFM includes both outflow (dotted line) and radiation (dashed lines) as described by Churazov et al. (2005) as well as a mass-dependent radiative efficiency following Davis & Laor (2011). The solid lines show the sum of ϵ_o and ϵ_r. The small dots and diamonds are observations by Davis & Laor (2011, D11) and Chelouche (2013, Ch13), who both estimated radiative efficiencies. In the radio regime we assume η = 0.1. The large stars and squares correspond to recent observations by Mezcua & Prieto (2014, MP14) of the outflow and radiation. From left to right the observed galaxies are M87, NGC 4594, NGC 1097, NGC 3169, NGC 1386, NGC 2911, NGC 1052 and Cen A. Small stars and squares correspond to observations by Russell et al. (2013, R13). The BH masses are colour-coded as indicated by the colourbar.

The three coloured lines in Fig. 2 show the model from Churazov et al. (2005) for β = 0.5 (thick dashed lines) and β = 1 (thin dashed lines) and different BH masses. The red lines correspond to M_• = 10¹⁰ M_⊙, the green ones to M_• = 10⁸ M_⊙ and the blue ones to M_• = 10⁶ M_⊙. This is in much better agreement with the observations than choosing a constant total efficiency. In the radio regime, we included observations by Russell et al. (2013) and Mezcua & Prieto (2014), who measured the power of the radiation and outflow as well as L_Edd. With equation (10) they calculated |$\dot{M}/\dot{M}_{\mathrm{Edd}}$|⁠. Using the equations (7), (8) and (11) we can derive the efficiencies

\begin{equation} \epsilon _\mathrm{o} = \eta \frac{P_0 / L_{\mathrm{Edd}}}{\dot{M_\bullet }/\dot{M}_{\mathrm{Edd}}} \end{equation}

(21)

and

\begin{equation} \epsilon _\mathrm{r} = \eta \frac{L / L_{\mathrm{Edd}}}{\dot{M_\bullet }/\dot{M}_{\mathrm{Edd}}}. \end{equation}

(22)

In the radio regime, it is justified to use η = 0.1. As can be seen in Fig. 2, the data points for the radiative efficiency do not show the simple trend as assumed in Churazov et al. (2005). In fact, they seem to be consistent with randomly scattering between 10⁻¹ and 10⁻⁵. There also seems to be no mass dependence in the radio regime.¹

For NGC 1097 and NGC 1386, the radiation dominates. The observations by Mezcua & Prieto (2014) show that these sources have small jets, whereas the other sources have larger jets. Interestingly, both NGC 1097 and NGC 1386 have a bar at large scales, but they show no evidence of a bar on small scales. They both also have a ring of star-forming regions. This indicates that the morphology of the galaxies will play a key role for future studies. For simulations, this implies that the resolution has to be high enough to resolve the morphology of galaxies. Note that this is not the case for the simulations performed in this work, but will be the aim for forthcoming studies.

3 THE SIMULATIONS

This work is based on a set of cosmological simulations called the Magneticum Pathfinder Simulations² (Dolag et al. in preparation). The simulations are performed with an updated version of the TreePM-SPH code p-gadget3 (Springel 2005).

We adopt a Λ cold dark matter cosmology with parameters according to the seven year results of the Wilkinson Microwave Anisotropy Probe with Ω_m = 0.272, Ω_Λ = 0.728, Ω_b = 0.0456 and h = 0.704 (Komatsu et al. 2011). We follow the hydrodynamics of the gas using the smoothed particle hydrodynamics method (see Price 2012 for a recent review on the SPH method). We use an entropy conserving formulation (Springel & Hernquist 2002), where star formation is based on a multiphase sub-resolution model by Springel & Hernquist (2003). Additionally, we include complex treatment for a wide range of physical processes such as isotropic thermal conduction (Dolag et al. 2004) with an efficiency of κ = 1/20 of the classical Spitzer value, stellar evolution, metal enrichment and supernova feedback (Tornatore et al. 2003; Tornatore et al. 2007), a cooling function which depends on the individual metal species following Wiersma, Schaye & Smith (2009) as well as the treatment of BHs and their associated feedback based on the model implemented by Springel et al. (2005). We improve the accuracy, stability and reliability of our hydrodynamical method with several state-of-the-art improvements of the SPH method. This includes the higher order Wendland kernel functions (Dehnen & Aly 2012) as well as time-dependent artificial viscosity to properly track turbulence within galaxy clusters (Dolag et al. 2005; Donnert et al. 2013).

Regarding the BH physics we use the modifications as described by Fabjan et al. (2010), in contrast to the original model implemented by Springel et al. (2005), and made changes to the seeding and further treatment of BHs as described in detail by Hirschmann et al. (2014). The most important one of these changes is that we do not pin the BHs to the most bound particles anymore. This ‘pinning’ is used in other simulations to keep the BHs in the centre of their host galaxy, but it also has the side effect that BHs ‘jump’ from the less massive galaxy to the more massive one during merger events. To avoid that the BH particles are wandering away from the centre of galaxies by numerical effects, we first implemented the conservation of momentum and centre of mass when two BH particles are merging. Secondly, we enforce momentum conservation for the smooth accretion of gas and therefore do not model any momentum transfer when swallowing gas. Without pinning, we have BHs not only in central galaxies, but also keep them in satellite systems until they fully merge. Thus, we are able to track BH growth much better, in particular in massive galaxy clusters (following all the BHs in satellite galaxies).

Hirschmann et al. (2014) already presented a detailed analysis of BH growth in the Magneticum Pathfinder Simulations particularly focusing on the origin of the antihierarchical growth of BHs within a hierarchical structure formation scenario. Various observational trends can be already explained using the simplified BH model described by Springel et al. (2005). However, implementing the more detailed description of AGN feedback and BH accretion as described in Section 2 leads to further improvements in predicting a more realistic population of BHs and AGN in our hydrodynamic simulations.

We performed six simulation runs with the same resolution as in the large (500 Mpc)³ box with an initial particle number of 2 × 1564³ analysed by Hirschmann et al. (2014). In the context of the set of Magneticum Pathfinder Simulations from Dolag et al. (in preparation), we refer to this resolution as hr (‘high resolution’). The particle masses are M_dm = 6.9 × 10⁸ M_⊙ h⁻¹, M_gas = 1.4 × 10⁸ M_⊙ h⁻¹ and M_stars = 3.5 × 10⁷ M_⊙ h⁻¹ and the softening length is 3.75 kpc h⁻¹ for dark matter and gas and 2.0 kpc h⁻¹ for stars. BHs are represented as collisionless sink particles. They are seeded in galaxies with stellar masses above 2.3 × 10¹⁰ M_⊙ with an initial mass of 4.6 × 10⁵ M_⊙.

Four of our simulations are ‘test’ runs with a smaller box size of (68 Mpc)³, which were performed to be able to test the effect of the new BH accretion and AGN feedback model separately. The first run adopts the ‘original’ BH model as described in Hirschmann et al. (2014) to which we refer as the fiducial model. The second run adopts only the new accretion model (NAM), the third run only adopts the new feedback model (NFM), and finally, our fourth run combines both new implementations (NFAM).

The other two simulations have the same resolution but a larger box size of (182 Mpc)³ to achieve a larger statistical sample of galaxies and BHs. The first box uses the original implementation of BH growth and the second box adopts the NFAM model, enabling us to statistically see the effects of the new model, in particular on the more massive galaxy and BH population.

As described in Section 2 in detail, the NAM, NFM and NAFM models contain improvements of the BH model regarding the calculation of the accretion rate and/or the feedback energy of BHs:

• NAM: for the estimation of the BH accretion rate we use different boost factors for cold (α = 100) and hot (α = 10) gas. For this run, we use the fiducial feedback model.

• NFM: for the calculation of the energy of the AGN feedback we consider not only radiative, but also mechanical feedback. The two different feedback mechanisms have different efficiencies. The radiative efficiency ϵ_r depends on the BH mass and the Eddington ratio, whereas the outflow efficiency ϵ_o depends only on the Eddington ratio. Like in the fiducial model only a fraction ϵ_f of the radiation couples to the surrounding medium. Both kinds of feedback are implemented as thermal feedback. Hence, the total feedback energy is computed with equation (9). We use the old accretion model for this simulation.

• NFAM: our final run contains both the new feedback and the NAM.

The NFM as shown in Fig. 2 was implemented into the code using equations (19) and (20). In reality the slope β can be between 0 and 1. However, the choice of β does not play a significant role for the simulations, as the mechanical outflow dominates over the radiation in the radio regime. Furthermore, the AGN luminosities are not calculated during the simulation but only for the analysis afterwards. Thus, we choose the fixed value of β = 0.5 for all simulations.

For the NAM run and the two fiducial runs, we use the standard feedback model with ϵ_f = 0.15 and a constant radiative efficiency ϵ_r = 0.2 (Hirschmann et al. 2014). In the other runs we use ϵ_f = 0.2. The parameters of the simulations used in this work are summarized in Table 1.

Table 1.

General settings of the simulations performed in this study. Variable values of ϵ_r and ϵ_o are calculated with equations (19) and (20).

	Box size	Initial particle number	ϵf	ϵr	ϵo
	[(Mpc/h)3]
68 Mpc/hr fiducial model	483	2 × 2163	0.15	0.2	–
68 Mpc/hr NFM	483	2 × 2163	0.2	Variable	Variable
68 Mpc/hr NAM	483	2 × 2163	0.15	0.2	–
68 Mpc/hr NFAM	483	2 × 2163	0.2	Variable	Variable
182 Mpc/hr fiducial model	1283	2 × 5763	0.15	0.2	–
182 Mpc/hr NFAM	1283	2 × 5763	0.2	Variable	Variable

Note that we identify the dark matter haloes and the corresponding galaxies in the simulation using the friends-of-friends and then the subfind algorithm (Springel et al. 2001; Dolag et al. 2009).

4 RESULTS

4.1 BH growth

4.1.1 BH–galaxy mass scaling relations at z = 0

The upper panel in Fig. 3 shows the predictions for the present-day M_•–M_* relation for the 68 Mpc/hr NFAM simulation. In our simulations M_* is the total stellar mass of a galaxy and not only the stellar mass of the bulge, because our resolution is not high enough to resolve the internal structures of the individual galaxies. Hence, all galaxies consist mainly of a spheroidal component. The solid black lines in Fig. 3 indicate the observations of McConnell & Ma (2013) and the dashed line is the fit for all BHs in our simulations with M_• > 5 × 10⁷. This threshold is necessary to exclude newly seeded BHs, as they are seeded far below the relation and need time to grow on to the relation. BHs with masses above M_• > 5 × 10⁷ are close enough to the M_•–M_* relation to exclude seeding effects. The figure shows the excellent agreement of our NFAM model with observations, in particular in comparison to other simulations, i.e. the Illustris simulation (Sijacki et al. 2014) and the MassiveBlack-II simulation (Khandai et al. 2014). The dark grey shaded area marks the 1σ scatter of the observations and the light grey shaded area the 1σ scatter for our simulation. For a quantitative comparison with the observations, Table 2 shows the best-fitting parameters a and b corresponding to the fit function log(M_•/M_⊙) = a + b · log(M_*/10¹¹M_⊙) for all six runs. It also contains the 1σ scatter of McConnell & Ma (2013) and our simulations. For the 182 Mpc/hr runs, we consider only stellar masses below 10¹² M_⊙ to exclude the central galaxies of very massive clusters (see discussion in Section 4.2).

Figure 3.

Upper panel: present-day relation between the BH mass and the host galaxy stellar mass for the 68 Mpc/hr NFAM run. The dots represent the BHs in the simulations at z = 0. The solid black line shows the fit to the observations by McConnell & Ma (2013, M&M13) and the dark shaded area the corresponding 1σ error. The dashed lines illustrate the fit to our simulation for M_• > 5 × 10⁷ M_⊙ (to exclude seeding effects) and the light shaded area the corresponding 1σ error. For comparison with other simulations, we also show the results from Sijacki et al. (2014, S14) and Khandai et al. (2014, K14) as dotted and dot–dashed lines. Lower panel: ratio of the simulated BH mass in all different models (fiducial: dark blue, NFM: light blue, NAM: green, NFAM: red) to the observed BH mass M_•obs (McConnell & Ma 2013, black solid line and grey shaded area) versus the galaxy stellar mass.

Table 2.

Best-fitting parameters and standard deviation for our runs in comparison to the observations by McConnell & Ma (2013). All BHs with masses smaller than 5 × 10⁷ M_⊙ have been excluded for the fit. For the 182 Mpc/hr runs, we took only stellar masses below 10¹² M_⊙ into account to exclude clusters.

	a	b	σ
McConnell & Ma (2013)	8.46 ± 0.08	1.05 ± 0.11	0.45
68 Mpc/hr fiducial model	8.53	1.28	0.17
68 Mpc/hr NFM	8.52	1.03	0.16
68 Mpc/hr NAM	8.44	1.24	0.19
68 Mpc/hr NFAM	8.51	1.00	0.16
182 Mpc/hr fiducial model	8.46	0.93	0.15
182 Mpc/hr NFAM	8.40	1.09	0.14

While the slope of the M_•–M_* relation turns out to be relatively insensitive to the values of ϵ_r and ϵ_f, the normalization depends strongly on these parameters as already shown by Di Matteo et al. (2005), because the final BH mass follows the proportionality M_• ∝ (ϵ_fϵ_r)⁻¹. Hence, many recent simulations which include BHs (e.g. Di Matteo et al. 2005; Robertson et al. 2006; Degraf et al. 2011; Hirschmann et al. 2014) tuned these parameters in order to reproduce the normalization of the observed M_•–M_* relation. In addition, the normalization depends on the cooling function (Churazov et al. 2005), i.e. the values of ϵ_r and ϵ_f must be larger to get the same normalization if the cooling is more effective. Since ϵ_r is not a constant parameter in our new AGN feedback model, the slope of the M_•–M_* relation changes. This is shown in the lower panel of Fig. 3. Here, we show the ratio of the simulated to the observed BH mass (from McConnell & Ma 2013) versus the galaxy stellar mass for all different models, i.e. the fiducial, NFM, NAM and NFAM runs (coloured dashed lines), as well as for the results from Sijacki et al. (2014) and Khandai et al. (2014, black dotted and dot–dashed lines, respectively). Since they use a constant radiative efficiency, their slopes are similar to our fiducial simulation. In our NFM, however, ϵ_r is not a free parameter anymore. Therefore, it is encouraging that both the slope and the normalization of the M_•–M_* relation are self-consistently predicted with less free parameters than in the standard model.

However, even in our new model one free parameter remains, i.e. the fraction of radiation coupling to the surrounding medium ϵ_f, for which we choose a value of ϵ_f = 0.2 (to be consistent with the observed relation).³ For lower efficiencies the feedback would be higher and the BHs would grow too much. We would like to remark that the normalization of the M_•–M_* relation in simulations always depends on the observations used for the calibration of ϵ_f. However, there are discrepancies in observational estimations of the M_•–M_* relation. For example, Scott, Graham & Schombert (2013) find a slightly higher normalization, but a similar slope as McConnell & Ma (2013), which would change the calibration of ϵ_f.

In our simulations, the NFAM model reproduces the observed slope better than the fiducial model, in which the BHs accrete slightly too much gas, resulting in too large masses, particularly at low redshifts and in the most massive galaxies. The new AGN feedback model is more efficient in preventing gas accretion on to massive BHs. Thus, the gas in the vicinity of the BH has a higher thermal and kinetic energy, which results in lower accretion rates. Consequently, as can be seen in Fig. 3, the massive end of the M_•–M_* relation is now in excellent agreement with the observations from McConnell & Ma (2013).

Our second implementation is the separation of hot and cold gas (NAM). For an increasing amount of hot gas in the vicinity of the BH, this results in slightly lower accretion rates due to the smaller boost factor. Even if the NAM by itself cannot prevent the most massive BHs from growing too much, it can decrease the BH masses slightly. Consequently, a combination of both modifications results in the best match with the observed M_•–M_* relation.

The best-fitting parameters in Table 2 summarize the excellent agreement of the NFM-run and the NFAM-run with the observations. Particularly, the slope b is in better agreement with the observations than in the other runs and also in the analysis of the Illustris simulation shown by Sijacki et al. (2014). Note that in the simulations, the 1σ scatter is significantly smaller than in the observations. As the typical measurement errors in the observations are still substantial, future observations are needed to distinguish, whether this relation indeed has such a small scatter as seen in the simulations, or if there are still additional processes missing in the simulations which influence the growth and evolution of the BHs.

Furthermore, the scatter in the BH mass in the simulations decreases with increasing BH mass. This is most likely a consequence of statistical merging (Peng 2007; Hirschmann et al. 2010; Jahnke & Macciò 2011) and is also visible in the Illustis simulation (Sijacki et al. 2014). Nevertheless, the relative role of AGN feedback and statistical merging in establishing the M_•–M_* relation and producing the observed slope still remains a matter of debate.

To explore BH growth in our simulations in more detail, Fig. 4 shows the cosmic evolution of four BHs selected due to their different present-day mass (different colours) on the M_•–M_* relation.⁴ When BHs are merging, the most massive progenitor is followed back in time. As can be seen in this figure, we can distinguish between two different phases of BH growth: during the first phase, they grow rapidly until they reach the M_•–M_* relation and thus the Eddington limit. In this phase BH accretion is primarily triggered by smooth accretion of cold gas, because below the Eddington limit AGN feedback is not strong enough to suppress gas cooling. Hence, the cold gas reservoir is large enough to trigger BH growth. In our simulations, this phase is a consequence of the small BH seeding mass. However, recent observations seem to indicate that the slope of the M_•–M_* relation is steeper for BHs with masses below 10⁸ M_⊙ (Graham & Scott 2013; Scott et al. 2013). Therefore, we can speculate that the phase of rapid BH growth is actually present and that simulations in which BHs are seeded on or above the M_•–M_* relation might miss the first phase of BH growth.

Figure 4.

Evolution of the total BH mass and the corresponding host galaxy stellar mass of four haloes (diamonds in different colours) in the 68 Mpc/hr NFAM simulation. The black line shows the fit from McConnell & Ma (2013).

In the second phase BHs grow along the M_•–M_* relation. In this phase, gas cooling and AGN feedback are in equilibrium and hence both star formation and BH growth are suppressed. Only the in-fall of cold gas either in the form of streams or clumps as well as merger events can trigger star formation and BH growth during this period.

To demonstrate that at low redshifts BHs grow faster compared to the growth of the stellar mass than at high redshifts, we show exemplarily the results for four typical objects, where we verified that they reflect the typical growth of BHs with the chosen final mass. For example, the stellar mass of the host galaxy corresponding to the red diamonds grows very little, whereas the BH mass increases by more than two orders of magnitude. This galaxy reaches the M_•–M_* relation already 1.08 Gyr after the seeding. In contrast, the stellar mass of the host galaxy corresponding to the black and blue diamonds grows much more during the first phase of BH growth. Here, the object reaches⁵ the M_•–M_* relation after 2.29 Gyr. This trend is also visible in Fig. 6, which shows the M_•–M_* relation at different redshifts, in particular when looking at the data points corresponding to the lowest stellar masses. The figure will be discussed later in more detail. Hence, we suspect that the BH mass at the threshold between the two phases – namely when the M_•–M_* relation is reached – depends on the seeding redshift. We suggest, that these differences might be a consequence of the star formation rate, which decreases with time (see Section 4.3).

Furthermore, since BHs are seeded upon a certain galaxy mass, they are seeded earlier in a dense environment and can thus become more massive. We plan to study the evolution of BHs and their host galaxies in a forthcoming study in more detail, performing a simulation with resolution high enough to resolve the internal structure of galaxies. In particular, we are interested in the effect of merger events on BH growth and star formation, because the BH and stellar masses in Fig. 4 seem to grow mainly in steps after reaching the M_•–M_* relation. These steps also explain the scatter around the M_•–M_* relation in our simulations. It furthermore indicates that BH growth and star formation are both triggered by merger events. However, for this study it is more important to increase the box size instead of the resolution, in particular to extend our simulation results towards more massive galaxies and BHs.

4.1.2 Evolution of the BH mass function

Fig. 5 shows the BH mass function of both the fiducial and the NFAM 182 Mpc/hr run. We compare our simulations to observed BH mass functions of the local Universe by Marconi et al. (2004), Shankar et al. (2004), Shankar et al. (2009) and Shankar (2013). We would like to remark that the uncertainties in these relations are large, in particular because the BH masses are estimated using different scaling relations as recently discussed by Shankar (2013) and therefore, we also show the BH mass functions derived from the best-fitting velocity dispersion function and stellar mass function from Bernardi et al. (2010) using different scaling relations, i.e. from McConnell & Ma (2013, dotted grey lines) and Kormendy & Ho (2013, dashed grey lines). Since the high-mass end of all of these curves is lower than in Shankar (2013), we take – following their discussion – the two data points at the high-mass end of Shankar (2013) as upper limits. One should also keep in mind that as discussed in Tundo et al. (2007), the different BH scaling relations are not necessarily consistent with each other or with the M_•–M_* relation from McConnell & Ma (2013), which we use in this work to calibrate the value of the free parameter ϵ_f. The uncertainties in the scaling relations are also reviewed and discussed in Kormendy & Ho (2013).

Figure 5.

BH mass function of the fiducial (dashed coloured lines) and the NFAM (solid coloured lines) 182 Mpc/hr simulation at different redshifts. For comparison, we show observations from Marconi et al. (2004, black solid line), Shankar et al. (2004, black diamonds and lines with grey shaded areas), Shankar et al. (2009, dark grey shaded area) and Shankar (2013, black dots). To show the uncertainties in deriving BH mass functions from observations, we show as dotted and dashed grey curves the BH mass functions derived from the best-fitting velocity dispersion function and stellar mass function from Bernardi et al. (2010) using different scaling relations, i.e. from McConnell & Ma (2013, MM) and Kormendy & Ho (2013, KK).

The high-mass end of the fiducial simulation is just in agreement with the upper limits of Shankar (2013), but the NFAM simulation matches previously published BH mass functions much better, because the new accretion and feedback models suppress the growth of massive BHs more efficiently. As already shown in Fig. 3, the smaller masses of the most massive BHs are mainly caused by the new feedback scheme, where the mass dependence of the radiative efficiency for the model is taken from Davis & Laor (2011), which is quite similar to the results presented in Trakhtenbrot (2014). From a theoretical point of view, this relation is motivated by the fact that the spin of the BH should increase with mass. However, the slope of this relation might actually be flatter than in Davis & Laor (2011) due to selection effects (see discussion in Raimundo et al. 2012 and Laor & Davis 2011). Thus, the massive end of the BH mass function of the NFAM simulation could be a lower limit. Furthermore, we already mentioned that it is uncertain whether in general the normalization of the M_•–M_* relation could be larger than in McConnell & Ma (2013).

For less massive galaxies, the effects of the seeding become dominant which cause the deviation from the observed BH mass function at small masses. However, especially at low masses, observations are uncertain and only give an upper limit (Shankar 2013), in particular because pseudo-bulges do probably not follow the observed scaling relations like the M_•–σ relation or the M_•–M_* relation as reviewed by Kormendy & Ho (2013).

4.1.3 Evolution of the BH–galaxy mass scaling relations

Fig. 6 shows the relation between the BH mass and the stellar mass of the host galaxy for our NFAM 182 Mpc/hr run at different redshifts, again in comparison to the observations by McConnell & Ma (2013) and the simulations from Sijacki et al. (2014) and Khandai et al. (2014). Again, we only show BHs with masses above 5 × 10⁷ M_⊙. Below this limit BHs generally grow fast, while M_* stays relatively constant until they reach the M_•–M_* relation. The reason is the equilibrium between AGN feedback and gas cooling, when BHs accrete with |$\dot{M}_{\mathrm{Edd}}$| as described by Churazov et al. (2005). Afterwards BHs can only grow along the M_•–M_* relation together with their host galaxy through smooth accretion or merging.

Figure 6.

Evolution of the relation between the BH mass and the host galaxy stellar mass for the NFAM 182 Mpc/hr run (red dots). The dashed lines are fits for both 182 Mpc/hr runs including all BHs with masses larger than 5 × 10⁷ M_⊙ and stellar masses with masses smaller than 10¹² M_⊙ to exclude clusters. The light grey shaded area marks the corresponding 1σ error of the NFAM run. The black line with the dark grey shaded area represents the fit through the observations from McConnell & Ma (2013) with the 1σ error. The dotted and dot–dashed lines show the results from other simulations, i.e. from Sijacki et al. (2014) and Khandai et al. (2014).

In the NFAM run, the M_•–M_* relation is much earlier in place than in the original run, namely already at z = 3. Furthermore, the panels at z = 2 and 1 show that in the fiducial simulation the slope of the M_•–M_* relation is larger than at z = 0, where it is in agreement with the observed M_•–M_* relation.

In our very massive galaxies (M_* ≈ 10¹³ M_⊙), i.e. the central galaxies of galaxy clusters, most BHs are lying slightly below the M_•–M_* relation. This is most likely caused by a still too large stellar mass in these very massive galaxies, also visible in the high mass excess of the stellar mass function and the still too large baryon conversion efficiency for large haloes as discussed later on. The reason for the overestimation of stellar masses of cluster galaxies might be the purely thermal feedback in our model, which fails to reproduce the mechanical feedback in such massive systems, visible as large X-ray cavities in observed clusters. Hence, an implementation of mechanical jets (e.g. Ostriker et al. 2010; Dubois et al. 2013; Choi et al. 2014b) might play an important role for future simulations, in which both the resolution and the size of the cosmological boxes will get larger and larger. Furthermore, in our analysis we do not distinguish between the stars belonging to the central galaxy and the ones which would be related to the intra cluster light, which can be substantial for such massive systems. It is also possible that some merging systems are identified as one galaxy. Thus, the predicted stellar mass for cluster galaxies might actually be slightly larger than in observations.

For comparison, Fig. 6 also includes the fit to the data points of the fiducial model, where BHs in galaxy clusters are substantially more massive compared to the stellar mass, especially at redshifts around z = 1. Although the fit at z = 0 is in agreement with the fit from McConnell & Ma (2013), it is evident from the BH mass function that the BH masses are too large at the high-mass end implying that the galaxy stellar masses must be too large (compensating for the large BH masses) which will be investigated in more detail in Section 4.2.

4.1.4 Eddington ratio distribution

The modifications in our NFAM simulations are also expected to significantly affect the Eddington ratios of the BHs. Therefore, in Fig. 7 we present the Eddington ratio distributions of both 182 Mpc/hr simulations at different redshifts. The black dotted vertical line shows the threshold between radio-mode and quasar-mode and the vertical lines in the top mark the mean values. For redshifts below z = 3, the Eddington ratios are clearly smaller in the NFAM run than in the fiducial simulation. For higher redshifts the Eddington ratios in the NFAM run are larger than in the fiducial simulation. We suggest that the wide range of values for the feedback efficiency leads to broader distributions. Especially the range of very low accretion rates is represented much better in the NFAM simulation than in the fiducial run.

Figure 7.

Eddington ratio distributions for the two 182 Mpc/hr simulations at different redshifts. The black dotted vertical line marks the threshold between radio-mode and quasar-mode. The vertical lines in the top show the mean values.

In contrast to the recent study from Sijacki et al. (2014) our simulations – in particular the NFAM run – show two peaks in the Eddington ratio distribution for z < 4, one in the radio-mode and a second peak either in the radio-mode or in the quasar-mode. This indicates that we have a clear separation between two accretion modes. In the fiducial model, where a step function was used to distinguish between radio-mode and quasar-mode (Hirschmann et al. 2014), the two peaks are only visible at z = 1. In the NFAM simulation, the second peak appears at z = 3 in the quasar-mode. For smaller redshifts it is much more distinct. Interestingly, at z = 1 and 2, which is the redshift range where most quasars are observed, a very clear second peak is visible in the quasar-mode. For z = 4 the Eddington ratios are even higher, because here the first phase of BH growth is dominant. At z = 0 both peaks are in the radio-mode and even a third peak is visible at very low Eddington ratios.

4.2 Evolution of the stellar mass function

Fig. 8 shows the evolution of the stellar mass function in the simulations (blue: fiducial model, red: NFAM model) and observations (black symbols from Panter et al. 2004, Cole et al. 2001, Bell et al. 2003, Pérez-González et al. 2008, Borch et al. 2006, Bundy et al. 2005, Drory et al. 2004, Fontana et al. 2006 and Marchesini et al. 2007 and black lines from Muzzin et al. 2013 and Bernardi et al. 2013). The figure illustrates that the new feedback scheme can slightly suppress late star formation at the high-mass end, mainly because the radiative efficiency now depends on the BH mass. Hence, compared to the fiducial model, the modifications in the NFAM model lower the amount of massive galaxies resulting in an overall better match with the massive end of the observed stellar mass function (SMF), at least down to z = 0.2.

Figure 8.

Stellar mass functions in different redshift ranges for the fiducial (blue lines) and the NFAM (red lines) 182 Mpc/hr runs. The solid black lines with the shaded areas show the observed stellar mass functions presented by Muzzin et al. (2013, M13) and their Poisson errors. The black diamonds are observations from Panter, Heavens & Jimenez (2004), Cole et al. (2001), Bell et al. (2003), Pérez-González et al. (2008), Borch et al. (2006), Bundy, Ellis & Conselice (2005), Drory et al. (2004), Fontana et al. (2006) and Marchesini et al. (2007). The black dashed and dot–dashed lines show the result from Bernardi et al. (2013, B13) using a Sérsic model and a Sérsic-bulge + exponential-disc model.

For the entire redshift range, a small peak in the SMFs is visible at stellar masses of about 2 × 10¹⁰ M_⊙. The origin of this peak is caused by a subtle effect of our BH seeding. Since BHs are seeded below the M_•–M_* relation, the AGN feedback is efficient during the first phase of BH growth and hence suppresses star formation until the equilibrium between cooling and AGN feedback is reached. During that phase, the stellar mass stops growing and consequently, there are more galaxies with a certain stellar mass. The peak moves towards higher stellar masses at higher redshifts because of the effect seen in Fig. 4, namely that BHs which are seeded earlier have larger stellar masses when they reach the M_•–M_* relation.

The overestimation of the low-mass end of the stellar mass function at high redshifts happens most likely due to the chosen wind model (constant winds as in Springel & Hernquist 2003) as described by Hirschmann et al. (2014) in more detail. Apart from that, our simulations – especially the NFAM run – are in good agreement with observations at high redshifts.

For z < 0.2, the high-mass end is still overestimated. However, we have to keep in mind that observations in this mass range contain also relatively large uncertainties. Bernardi et al. (2013) showed that different measurements of stellar masses differ from each other significantly, especially at the high-mass end. They demonstrate that the stellar masses are higher using a Sérsic model instead of standard models. Their fits using a single Sérsic and a Sérsic-bulge + exponential-disc model are shown as black dashed and dot–dashed line in the upper-left panel of Fig. 8. In comparison to other observational estimates this is in better agreement with our simulations. Nevertheless, the high-mass end still appears to be slightly overestimated in our simulations as also indicated by the massive end of the M_•–M_* relation (see lower-right panel of Fig. 6).

To study the effect of our new accretion and feedback models on the stellar masses in more detail, Fig. 9 shows the stellar mass functions separately for quiescent and star-forming galaxies in our simulations – again in comparison to the observations from Muzzin et al. (2013). Following Franx et al. (2008), we use a specific star formation rate of 0.3/t_Hubble as threshold to distinguish between quiescent and star-forming galaxies. We would like to mention that this is a different selection criterion than in the observations, where a threshold in the UVJ diagram is used (Muzzin et al. 2013). Hence, this criterion might lead to discrepancies with the observations, which may e.g. falsely identify metal-rich, star-forming galaxies to be red and thus quiescent.

Figure 9.

Stellar mass functions of quiescent (dashed lines) and star-forming (solid lines) galaxies in different redshift ranges for the fiducial (blue lines) and the NFAM (red lines) 182 Mpc/hr runs. For the threshold between quiescent and star-forming galaxies, we use the specific star formation rate of 0.3/t_Hubble following Franx et al. (2008). The black lines with the shaded areas (light grey for star forming and dark grey for quiescent galaxies) show the observations from Muzzin et al. (2013, M13) and their Poisson errors.

Fig. 9 illustrates that our new implementations increase the amount of quiescent galaxies at z > 1.5. Consequently, for this redshift range, the discrepancies between simulated and observed SMFs are much smaller for the NFAM simulation than for the fiducial run. Star formation is suppressed, when cooling and AGN feedback are in equilibrium (Churazov et al. 2005) and the gas in the vicinity of the AGN cannot cool enough to form stars. Hence, the increase of the amount of quiescent galaxies can be explained with the upper-left panel in Fig. 6, which shows that the M_•–M_* relation – and thus the phase of equilibrium – is earlier in place for the NFAM run. This is due to higher BH accretion rates during the phase of rapid BH growth as a consequence of both new implementations: first, the NAM leads to higher accretion rates when cold gas dominates. Secondly, the NFM results in less AGN feedback for low BH masses and thus to lower gas temperatures.

In contrast to the equilibrium phase, which can be associated with the radio-mode, the phase of star formation and rapid BH growth is not much affected by our new implementations. We conclude that the overestimation of the high-mass end is mainly due to star-forming galaxies. At z < 1 the amount of star-forming galaxies is too low for 2 × 10¹⁰ M_⊙ < M_* < 2 × 10¹¹ M_⊙. First, this is an effect of the low seeding mass of BHs, which also leads to the overproduction of quiescent galaxies. Secondly, it is a consequence of the overestimation of the high-mass end.

For both runs, Fig. 9 shows an artefact at low redshifts, namely that the amount of star-forming galaxies decreases rapidly after the seeding of BHs. We speculate that this decrease might be due to our very low BH seeding mass, which leads to artificially high accretion rates. This also explains why the number of star-forming galaxies is reduced in the NFAM model compared to the fiducial one. Fig. 4 illustrates why this artefact becomes even larger with decreasing redshift: for BHs that are seeded later, the evolutionary track during the first phase of BH growth is steeper than for early BH seeds. All in all Fig. 9 shows that our new implementations cannot significantly improve the stellar mass functions at low redshifts, but at high redshifts they predict a larger amount of quiescent galaxies, which is in better agreement with observations.

To quantify how efficient baryons are converted into stars for a given halo mass, we calculate the mean baryon conversion efficiencies, which are defined as M_*/(f_barM_halo), where f_bar = 0.17 is the baryon fraction of the Universe, for different redshifts. To be comparable to other studies we do not use M_vir for the halo mass, but M_200c, which is the mass inside the radius where the density is 200 times larger than the critical density of the Universe. Fig. 10 shows the conversion efficiencies versus halo mass for our two 182 Mpc/hr runs (different panels illustrate z = 0, 1, 2). The black vertical line shows the resolution limit for the baryon content as estimated by Vazza et al. (2011), which is given by 500 dark matter particles. Furthermore, the dashed and solid red vertical lines mark the minimum and mean value of M_200c, respectively, in the NFAM simulation corresponding to the minimum stellar mass for BH seeds. Below the mean seeding limit our resolution does not allow reliable predictions (dashed lines). The figure clearly shows, that the new implementations lower the stellar content in a halo for a given mass above this limit, which is also reflected by the reduced high-mass end of the stellar mass functions (see Fig. 8). At z = 2 and 1, this effect is even stronger than at z = 0. The dotted and dot–dashed black lines show the predictions of the abundance matching models by Moster et al. (2013) and Behroozi et al. (2013). The peak at M_halo ≈ 10¹² M_⊙ is in agreement with these models, which also find a maximum baryon conversion efficiency of around 20 per cent. At larger halo masses, the stellar content decreases due to AGN feedback and because the gas is consumed by star formation. Although the baryon conversion efficiencies in the NFAM simulation are smaller than in the fiducial run, they are still higher than in the abundance matching models of Moster et al. (2013) and Behroozi et al. (2013) for M_200c > 10¹³ M_⊙ galaxies.

Figure 10.

Mean baryon conversion efficiencies versus halo mass at different redshifts for the two 182 Mpc/hr runs. The grey shaded area shows the 1σ error of the NFAM run. The dashed and solid red vertical lines mark the minimum and mean value of M_200c in the NFAM simulation corresponding to the minimum stellar mass for BH seeds. Below the mean seeding limit our resolution does not allow reliable predictions (dashed lines). The black vertical line shows the resolution limit for the baryon content as estimated by Vazza et al. (2011), which is given by 500 dark matter particles. We compare our simulation with abundance matching models (Behroozi, Wechsler & Conroy 2013; Moster, Naab & White 2013) and with observations estimating the halo mass with weak lensing (Mandelbaum et al. 2006; Reyes et al. 2012; Hudson et al. 2013) or X-ray temperatures (Kravtsov, Vikhlinin & Meshscheryakov 2014).

For the NFAM simulation, at low redshifts a slight ‘upturn’ of the baryon conversion efficiencies occurs for stellar masses above 10¹⁴ M_⊙ corresponding to galaxy clusters due to too inefficient AGN feedback. This might indicate that other AGN feedback processes like mechanical jets should be included in future simulations. Since the most massive BHs accrete less in the NFAM model, we suspect that there is more cold gas left to form stars than in the fiducial run. Therefore, the upturn is only visible in the NFAM simulation. However, except for the high-mass end, our simulations – in particular the NFAM run – are in agreement with observations using weak lensing (Mandelbaum et al. 2006; Reyes et al. 2012; Hudson et al. 2013) or X-ray temperatures Kravtsov et al. (2014) to estimate the total halo mass.⁶

4.3 Evolution of the star formation rate

Fig. 11 shows the SFR–stellar mass plane (number density is colour-coded) for our two 182 Mpc/hr runs at different redshifts. The panels illustrate all galaxies classified as subhaloes using the subfind algorithm (Springel et al. 2001; Dolag et al. 2009). For comparison with observations, we also show the main sequence for star-forming galaxies estimated by Steinhardt et al. (2014) for 4 < z < 6 (red line), by Daddi et al. (2007) for z = 2 (orange line) and by Elbaz et al. (2007) for z = 1 and 0 (yellow line). At z = 2 and 1, the simulated SFRs at a given stellar mass are slightly below the observations. This trend is also visible in the recently published analysis of the Illustris simulation by Sparre et al. (2014). At z = 0 and at redshifts above z = 4 our simulation results are in very good agreement with the observed main sequence, independent of the adopted BH model. The redshift evolution of the SFR–stellar mass plane nicely demonstrates that the most massive galaxies become more and more quiescent with cosmic time. Furthermore, in the NFAM simulation star formation is suppressed earlier than in the fiducial one. This is consistent with Fig. 9, where we demonstrated that in the NFAM run the amount of quiescent galaxies is larger at earlier times. In the NFAM simulation, the SFRs of the most massive galaxies decrease already at redshifts above z = 4.8 such that they lie below the observed main sequence of star-forming galaxies. In the fiducial simulation, this decrease starts at redshifts below z = 4. This may be unrealistic, because – as shown in Fig. 9 – Muzzin et al. (2013) observe much more quiescent galaxies at high redshifts (z > 3) than in our fiducial simulation. Looking at the star formation main sequence of the Illustris simulation (Sparre et al. 2014) shows that this is not only a problem in our fiducial run, but seems to be a general issue. Therefore, it is encouraging that in the NFAM run galaxies become quiescent much earlier due to both of our new implementations, even if there are still discrepancies between the observed and simulated SMFs for star-forming and quiescent galaxies. The NFM leads to a lower feedback energy for low BH masses, whereas for large BH masses the AGN feedback is stronger as long as the BHs are accreting in the quasar-mode and star formation is suppressed.

Figure 11.

Comparison of the star formation rates of all galaxies in the two 182 Mpc/hr runs at different redshifts. The solid lines represent the observed main sequence of galaxies derived by Steinhardt et al. (2014, S14), Daddi et al. (2007, D07) and Elbaz et al. (2007, E07).

The NAM leads to lower accretion rates when the hot gas phase dominates. Hence, BHs grow less strongly and the SFR decreases already in less massive galaxies as can be seen in the panels corresponding to z = 1. From the earlier and more rapid decrease of the SFR follows that at z = 1 star-forming galaxies with stellar masses above 2 × 10¹⁰ M_⊙ are more concentrated along the observed main sequence in the NFAM simulation than in the fiducial one. At z = 0 there are only very few star-forming galaxies above log(M_*/M_⊙) = 10.5, which is the mass at which AGN feedback becomes important. At that redshift both runs predict galaxies with similar SFRs at a given stellar mass. Hence, our modifications mainly affect the evolution of high-redshift galaxies.

Fig. 12 depicts the redshift evolution of the mean specific SFR for our two 182 Mpc/hr runs. As in Biffi & Maio (2013) – who studied early proto-galaxies at z > 9 – we compare our simulations with other theoretical models (i.e. Davé, Oppenheimer & Finlator 2011; Biffi & Maio 2013; Dayal et al. 2013) and observations (i.e. Daddi et al. 2007; Noeske et al. 2007; Dunne et al. 2009; Pannella et al. 2009; Stark et al. 2009, 2013; Yabe et al. 2009; Michałowski, Hjorth & Watson 2010; Schiminovich et al. 2010; González et al. 2012; Reddy et al. 2012; Zheng et al. 2012; Coe et al. 2013; Bouwens et al. 2014). Irrespectively of the assumed accretion and feedback models, our simulations are both in better agreement with observations than many other theoretical models, especially at low redshifts (where the observational constraints are tighter). Fig. 12 also demonstrates that our new implementations have no effect on the specific SFR. Hence, the changes in the SFR and in the stellar mass are the same. However, star formation is certainly not only regulated by AGN feedback. Recent studies (e.g. Aumer et al. 2013; Hirschmann et al. 2013; Hopkins et al. 2014; Kannan et al. 2014) showed that stellar feedback also plays an important role, particularly for low-mass galaxies. Fig. 13 provides further evidence that our model is still not sufficient for reproducing galaxies with realistic SFRs. It illustrates the history of the star formation and the BH accretion rate densities as shown by Hirschmann et al. (2014) for our two 182 Mpc/hr runs compared to observations of the SFR density (squares) by Hopkins & Beacom (2006). In comparison to the fiducial model, the star formation rate density in the NFAM model is slightly lower above z ≈ 1.5, although it is still too high in comparison to the observations except for very high redshifts, which are, however, affected by resolution.

Figure 12.

History of the specific star formation rate in our 182 Mpc/hr runs in comparison to different observations and other theoretical predictions.

Figure 13.

History of the star formation (orange lines) and BH accretion rate (red lines) density in both 182 Mpc/hr runs (fiducial model: dashed lines, NFAM: solid lines) in comparison to observations from Hopkins & Beacom (2006, squares).

As expected due to the lower BH masses in the NFAM model, the BH accretion rate density is significantly lower at z < 4.5 than in the fiducial model. For higher redshifts, it is larger than in the fiducial model, which leads to a much shallower increase up to the maximum. Fig. 13 demonstrates that in the NFAM simulation the SFR and the BH accretion rate evolve very similar with redshift. The reason is that both depend on the amount of cold gas. With our NAM the analogy between SFR and BH accretion is even stronger, because the accretion factor for hot gas is smaller than for cold gas. Thus, in the NFAM simulation, hot gas results not only in less star formation, but also in smaller BH accretion rates. This shows that the gas temperature plays a key role in both galaxy formation and BH growth. A similar accordance between the history of the star formation and BH accretion rate density was also found by Zheng et al. (2009), who adopted the luminosity functions from Hopkins et al. (2007) to estimate the BH accretion rate densities.

5 DISCUSSION

5.1 The effect of the feedback model on to the luminosity functions

As already mentioned before, the choice of the slope β of the feedback model should not have a significant influence on the resulting galaxy and BH properties in the simulations since ϵ_r is much smaller than ϵ_o. However, it has an influence on the AGN luminosity functions, which are calculated during post-processing using the accretion rates calculated by the simulation and the radiative efficiencies, which can be varied.

In that way we can test the effect of the parameter β on the AGN luminosity function. We calculate the bolometric AGN luminosities of the NFAM simulation for different values of β using equations (8) and (19). Fig. 14 shows the resulting luminosity functions in comparison to the observational compilation of Hopkins et al. (2007). For a comparison of moderately luminous AGN, particularly at high redshifts, one has to keep in mind that simulations are affected by resolution (see discussion of Hirschmann et al. 2014). In addition, dust obscuration effects in observational data typically result in an underestimation of their number density (e.g. Hasinger 2008; Merloni et al. 2014) which complicates a comparison between simulations and observations. Even if luminosity-dependent obscuration effects on a torus level are already considered in Hopkins et al. (2007), an additional redshift dependence (of X-ray luminosities, as suggested by e.g. Hasinger 2008 and Merloni et al. 2014) may change the low-luminous end at high redshifts.

Figure 14.

AGN luminosity function of our 182 Mpc/hr NFAM run at different redshifts for different values for the slope β in comparison to the observational compilation by Hopkins et al. (2007).

Fig. 14 shows that the effect of the choice of β on the AGN luminosity functions is not significant, especially at high redshifts, because β changes only the efficiencies in the radio-mode and not in the quasar-mode. For lower redshifts, when more BHs accrete with low Eddington ratios, it has an influence on the amount of AGN with luminosities smaller than 10⁴⁵ erg s⁻¹ in the sense that with decreasing β the radiative efficiency and thus the amount of moderately luminous AGN is increasing and thus the result is in better agreement with the observational constraints. However, due to the fact that observations constrain very low values of ϵ_r, we suspect that the accretion rates in the quasar-mode are slightly underestimated in our simulations.

As shown in Fig. 2, the actual value of ϵ_r is entirely unconstrained in the radio regime. It might depend on many properties like the morphology of the host galaxy or the merger history of an individual BH. For that reason, calculating a more realistic value of ϵ_r is beyond the current feasibility.

Nevertheless, according to the observations by Russell et al. (2013), one should consider different models to estimate ϵ_r in the radio-mode. Fig. 15 shows the AGN luminosity functions in comparison to observational compilation by Hopkins et al. (2007) for four models adopting different values for ϵ_r in the radio regime.

• The commonly used value ϵ_r = 0.1 (green lines) seems to match the observations reasonably well, although such a value is unlikely according to the results from Russell et al. (2013) and Mezcua & Prieto (2014).

• ϵ_r = 10⁻³ is the mean value of the data points from Russell et al. (2013). Because we change only values in the radio regime, the high-luminosity end is not affected. At lower luminosities, the AGN number densities are significantly underestimated as AGN become way too faint⁷ (blue lines).

• We choose random values in log space in the range 10⁻⁵ < ϵ_r < 0.4. This is approximately the range of the data points from Russell et al. (2013) with a maximum value equal to the theoretical maximum efficiency of a rotating BH (since we assumed η = 0.1). It leads to a reasonably good match (magenta line) with the observational constraints, even if the low-luminous end is slightly lower than when adopting the commonly used value (green lines). Since we may speculate that the curve will probably be shifted upwards when choosing a higher resolution (Hirschmann et al. 2014), the concordance with the observations might be even better.

• Now we exclude very low values for ϵ_r and hence choose random values in the range 10⁻³ < ϵ_r < 0.4. This leads to a slightly, but not significantly larger number density of moderately luminous AGN (red lines) and hence to a better agreement with observations.

Figure 15.

AGN luminosity function of our 182 Mpc/hr NFAM run at different redshifts for different values of ϵ_r in the radio regime in comparison to the observational compilation by Hopkins et al. (2007). The green and blue curves show the result for two constant values of ϵ_r. For the purple and red curve we took random values in two different intervals.

In comparison to the AGN luminosity functions from the Illustris simulation (Sijacki et al. 2014), we have less luminous AGN for redshifts below z = 1, although our cosmological box is larger. Nevertheless, to investigate the high-mass end in more detail larger cosmological boxes are needed. Hirschmann et al. (2014) already presented luminosity functions of a larger box from the set of Magneticum Pathfinder Simulations, which are in good agreement with the observations from Hopkins et al. (2007). Furthermore, our simulation matches better with the observed amount of AGN with luminosities below L ≈ 10⁴⁵ erg s⁻¹ than in Sijacki et al. (2014). This confirms the conclusion from Sijacki et al. (2014) that the radiative efficiency is not constant and might actually be very low in the radio regime.

This analysis shows that the efficiency of the radiative component in the radio regime is indeed not yet understood because the theoretical lower limit is not captured by observations. Interestingly, choosing random values for the radiative efficiency in the range of the observed values leads to a good agreement with observed AGN luminosity functions. This may indicate that in the radio regime the radiative efficiency depends neither on the mass of the BH, nor on its accretion rate. It also implies that – as we are matching the observed luminosity function by randomly choosing the radiative efficiency within the observed values – the distribution of the accretion rates as predicted by the simulations are similar to the observed ones. We conclude, that it is theoretically not fully understood how efficient AGN radiate and we suspect that the morphology of the galaxy, but also turbulence or even magnetic fields might play an important role. Since jets dominate in the radio-mode, they can also prevent efficient accretion. The similar morphologies of the two radiation dominated sources from Mezcua & Prieto (2014), i.e. NGC 1097 and NGC 1386, give a first evidence for these speculations, because they both have a ring of star-forming regions and a bar on large scales, but no bar on small scales. However, a better understanding of BH accretion and AGN feedback processes is a great challenge for the future, because more accurate observations are needed to learn in which cases ADAF/Bondi models are a good estimate and in which cases we have to include additional physical processes.

5.2 The unconstrained total efficiency in the radio regime

Besides the radiative efficiency, the total efficiency η in the radio regime is also unconstrained. Throughout this study, we always assumed η = 0.1 to calculate ϵ_r and ϵ_o, making, thus, our conclusions for the radio regime rather uncertain. The reason for this assumption are missing or unconstrained estimations of |$\dot{M_\bullet }$|⁠. According to equation (11), η is given by

\begin{equation} \eta = \frac{L_{\mathrm{Edd}}}{\dot{M}_{\mathrm{Edd}} c^2} = \frac{L_{\mathrm{bol}} \frac{\dot{M_\bullet }}{\dot{M}_{\mathrm{Edd}}}}{\frac{L_{\mathrm{bol}}}{L_{\mathrm{Edd}}} \dot{M_\bullet } c^2}. \end{equation}

(23)

In observations, however, usually only the AGN luminosity, the jet power and the BH mass are measured. Using the BH mass, one can calculate L_Edd. Equation (10) is then used to calculate |$\dot{M_\bullet }/\dot{M}_{\mathrm{Edd}}$|⁠. Hence |$\dot{M_\bullet }$| is the parameter which is typically missing. Nevertheless, for some of the sources from Russell et al. (2013) and Mezcua & Prieto (2014), |$\dot{M_\bullet }$| has been estimated. We use these estimations to calculate the corresponding total efficiencies with equation (23). With these values and equations (21) and (22), we then compute ϵ_o and ϵ_r.

Before we calculate the efficiencies for the selected sources, we want to focus on the nearest SMBH, namely Sagittarius A* (Sgr A*). For the luminosity we adopt L_bol = 2.1 × 10³⁶  erg s⁻¹ (Narayan et al. 1998) and for the power of the mechanical outflow we assume P_o = 1.2 × 10⁴¹ erg s⁻¹ (Yusef-Zadeh et al. 2012). With these values and the mass M_SgrA* = 4 × 10⁶ M_⊙, we calculate the Eddington ratio using equation (10). Although Sgr A* is the nearest SMBH, there are different estimates for the accretion rate. Quataert, Narayan & Reid (1999) estimated a Bondi accretion rate of ∼3 × 10⁻⁵ M_⊙ yr⁻¹. However, there are other models suggesting the actual accretion rate might be much lower than the Bondi accretion rate (e.g. Quataert & Gruzinov 2000). Cuadra et al. (2006) derived |$\dot{M} \approx 3 \times \ 10^{-6} \,\mathrm{M}_{{\odot }}\,\mathrm{yr}^{-1}$| from stellar winds. We calculated the efficiencies corresponding to both values using equations (21) and (22). They are shown in Fig. 16. The upper data points belong to |$\dot{M} \approx 3 \times \ 10^{-6} \,\mathrm{M}_{{\odot }}\,\mathrm{yr}^{-1}$| and the lower ones to |$\dot{M} \approx 3 \times \ 10^{-5} \,\mathrm{M}_{{\odot }}\,\mathrm{yr}^{-1}$|⁠. Assuming that the ADAF model really provides a lower limit, this illustrates that |$\dot{M} \approx 3 \times \ 10^{-6} \,\mathrm{M}_{{\odot }}\,\mathrm{yr}^{-1}$| is in good agreement with our model for the radiative efficiency. It also indicates that it is necessary to choose different lower limits for different BH masses, because the dashed green line – which corresponds to η ≈ 0.1 – is far above the data point. However, the corresponding value for ϵ_o is larger than the commonly used value 0.1. This indicates, that the outflow efficiencies might differ significantly from this value, which is not well constrained. For the second estimation of the accretion rate, i.e. |$\dot{M} \approx 3 \times \ 10^{-5} \,\mathrm{M}_{{\odot }}\,\mathrm{yr}^{-1}$|⁠, the radiative efficiency is clearly below the prediction, although ϵ_o is near 0.1. This implies that Bondi estimations of the accretion rate indeed tend to be too high.

Figure 16.

Same as in Fig. 2, but with efficiencies calculated using values for |$\dot{M_\bullet }$| from Russell et al. (2013, R13) and from other authors, i.e. Evans et al. (2004), Allen et al. (2006) and Li et al. (2011) (R13*, MP14). The three data points from Mezcua & Prieto (2014), for which we know estimations of |$\dot{M_\bullet }$| are from left to right M87, NGC 4594 and CenA. We also included values for Sgr A*, which have been calculated using different estimations of |$\dot{M_\bullet }$|⁠, i.e. |$\dot{M} \approx 3 \times \ 10^{-6} \,\mathrm{M}_{{\odot }}\,\mathrm{yr}^{-1}$| from Cuadra et al. (2006, upper symbols) and |$\dot{M} \approx 3 \times 10^{-5} \,\mathrm{M}_{{\odot }}\,\mathrm{yr}^{-1}$| from Quataert & Gruzinov (2000, lower symbols).

Now, we consider the sources from Russell et al. (2013) and Mezcua & Prieto (2014), for which |$\dot{M_\bullet }$| has been estimated using the Bondi model. Russell et al. (2013) investigated a subsample of 13 objects for which they estimated |$\dot{M_\bullet }$|⁠. The efficiencies corresponding to these sources are plotted in Fig. 16 (R13). Other authors also estimated |$\dot{M_\bullet }$|⁠: for Centaurus A and NGC 4216, we use the result from Evans et al. (2004) and for the Sombrero galaxy (NGC 4594) we take |$\dot{M_\bullet }$| from Li et al. (2011). For M87, M84, M89, NGC 4636, NGC 4472, NGC 407 and NGC 5846 we take values from Allen et al. (2006). The efficiencies calculated with these values and the data from Russell et al. (2013) are marked with grey symbols (R13*). Most of these sources are also in the selected sample from Russell et al. (2013). We can, thus, directly compare the results of two independent measurements. This shows a clear discrepancy between different estimations of |$\dot{M_\bullet }$|⁠. Overall, the efficiencies are larger using the |$\dot{M_\bullet }$| from Russell et al. (2013). In contrast to Fig. 2, the lowest values of the radiative efficiency now tend to increase with increasing Eddington ratio as predicted by theory. Nevertheless, the observations are in better agreement with theory using only the 13 objects of the selected subsample. Furthermore, Fig. 16 indicates that the value ϵ_o = 0.1 is indeed a reasonable assumption for the mean value of the observed values, although the observations can be nearly two dex lower.

However, all these estimations are highly uncertain and very speculative. On the one hand, all data points are upper limits due to the approximation of using the Bondi model. On the other hand, there are studies showing that accretion rates can also be much smaller than |$\dot{M}_\mathrm{B}$| (i.e. Quataert & Gruzinov 2000; Baganoff et al. 2003; Li et al. 2013). Moreover, values for L_bol might be underestimated when the jet is emitting in the plane of the sky. In that case, the measured flux is smaller than if the jet were located close to the line of sight. This would lead to higher radiative efficiencies and to an even better agreement with our model. Furthermore, uncertainties in the determination of BH masses make it almost impossible to investigate whether the lower limit for the radiative efficiency splits up for different BH masses as seen in the quasar regime (Davis & Laor 2011; Chelouche 2013).

Nevertheless, the data shown in Fig. 16 is one of the best constrained samples. The comparison between Figs 16 and 2 shows that we need more accurate measurements to learn more about the feedback of radio jets and the corresponding efficiencies. Due to the fact that knowing the efficiencies is (at least with the currently available computational power) essential for performing large-scale cosmological simulations, it is worth and necessary spending more effort on observational estimates of BH accretion rates.

5.3 Comparison with other simulations

During the last couple of years, several other groups have also been working on large cosmological simulations including baryons and BHs. As our simulations, some of these simulations, for example the MassiveBlack-II simulation (Khandai et al. 2014), earlier simulations from Di Matteo et al. (2008) and the new EAGLE simulation (e.g. Schaye et al. 2015), are based on the SPH code gadget-3, but differ in their physical sub-resolution models, including the model for BH growth. In contrast, the recent Illustris simulation (e.g. Genel et al. 2014; Vogelsberger et al. 2014a) has been performed with a different hydrodynamic scheme, the moving mesh code arepo (Springel 2010), and also slightly different sub-resolution models. A comparison between these models can help to understand which effects the different sub-resolution models for BH growth and AGN feedback may have on basic galaxy and BH properties.

Fig. 17 shows the stellar mass function in the NFAM model below z = 0.2 in comparison to other simulations. As for the BH mass function, the number density of massive galaxies in the Illustris simulation (Genel et al. 2014, green lines) is by half an order of magnitude larger than the one in the Magneticum simulation. For stellar masses below 4 × 10¹¹ M_⊙, the galaxy number densities in the Illustris simulation are in reasonably good agreement with the observations, while our simulations produce slightly too few low-mass galaxies. Since the difference between the SMFs of the fiducial model and the NFAM model are very small at z = 0, we suggest that other physical processes (e.g. stellar feedback or cooling) or the lower resolution might be the reason for the lower stellar masses. The prediction from the MassiveBlack-II simulation (Khandai et al. 2014, orange line) has no pronounced exponential cut-off with the consequence that they overestimate the low- and the high-mass end, but slightly underestimate the number density of galaxies around the exponential cut-off. In contrast, the stellar mass function obtained by the EAGLE simulation (Schaye et al. 2015, blue line), where the feedback is especially calibrated to match the stellar mass functions, is in good agreement with observations for the entire stellar mass range.

Figure 17.

Comparison of the SMF in the NFAM model (red line) at z = 0 with the Illustris simulation (Genel et al. 2014, G14, green line), the MassiveBlack-II simulation (Khandai et al. 2014, K14, orange line) and the EAGLE simulation (Schaye et al. 2015, Sch14, blue line). The observations shown are the same as in Fig. 8.

Compared to our results – the BHs in the Illustris simulation are much more massive than in the Magneticum simulation (as shown in Fig. 18). This discrepancy might have several reasons, for example the different implementations of radiative AGN feedback. Furthermore, given that there may still be resolution-dependent details of the BH feedback model (e.g. the estimation of the Bondi accretion rate or the distribution of the feedback) the higher resolution of the Illustris simulation could contribute to these differences. In addition, there could be differences due to the different numerical techniques, namely SPH and moving mesh, especially in the way the feedback gets transported away from the centre of the galaxies. In addition, a more efficient gas cooling in arepo (Nelson et al. 2013) might lead to higher BH accretion rates. Furthermore, the underlying physics referring to the energy transport might influence how much gas is driven outwards and which fraction of this gas is recycled as for instance discussed by Nelson et al. (2014).

Figure 18.

Comparison of the BH mass function in the 182 Mpc/hr NFAM run with that in the Illustris simulation (Sijacki et al. 2014) at z = 0 and with that in the MassiveBlack-II simulation (Di Matteo et al. 2008) at z = 1. The observations shown are the same as in Fig. 5.

Due to the large uncertainties in different observational estimates, it is not clear which simulation matches the observations of the local Universe best. At z = 1, we also compare the BH mass function of our NFAM model to the predictions of Di Matteo et al. (2008). This simulation produces slightly more massive BHs than the Magneticum simulation, which might be due to a more inefficient AGN feedback of massive BHs in Di Matteo et al. (2008).

Obviously, the other simulations shown here capture BHs down to smaller BH masses. First, this is due to the higher resolutions. Secondly, they use the so-called ‘pinning’ to keep the BHs at the potential minimum and therefore in the centre of the galaxies. Hence, they can seed the BHs in less massive galaxies. In our simulations this is not possible, because the BHs in less well-defined galaxies would not be able to stay in the centre of their host galaxy due to numerical effects. However, not using the so-called ‘pinning’ avoids other drawbacks of this method as discussed in Hirschmann et al. (2014). As discussed by Shankar (2013), also the low-mass end of the BH mass function is relatively uncertain and depends on the BH scaling relations. For example, the low-mass end could be significantly smaller when excluding galaxies with pseudo-bulges. Therefore, it will be quite challenging to compare observed BH mass functions to any simulation at the low-mass end.

Fig. 19 shows a comparison of the AGN luminosity function in our NFAM run (purple line) with the predictions from the Illustris simulation (Sijacki et al. 2014, green solid line) and from the MassiveBlack-II simulation (Khandai et al. 2014, orange solid line). The luminosity function of the Illustris simulation matches both the observations and our simulation, whereas the MassiveBlack-II simulation widely fails to reproduce the observed shape of the observed luminosity functions of Hopkins et al. (2007). Since the latter simulation contains the original model from Springel et al. (2005) with only one mode of AGN feedback, we can speculate that this might be one possible reason for the discrepancies. The Illustris simulation uses a so-called ‘radiative’ efficiency, which is implemented as a change in the net cooling rate and is most efficient in the quasar-mode (Sijacki et al. 2014). This seems to have a similar effect as our variable radiative efficiency, which increases for large BH masses in the quasar mode.

Figure 19.

Comparison of the AGN luminosity function in the NFAM model (using random radiative efficiencies in the radio regime in the range 10⁻⁵ < ϵ_r < 0.4) with the predictions of the Illustris simulation (Sijacki et al. 2014, S14, green solid lines) and of the MassiveBlack-II simulation (Khandai et al. 2014, K14, orange solid lines).

Nevertheless, we want to emphasize that despite of the general importance for understanding the (physical or numerical) origin of different simulation predictions, such a comparison must remain speculative: besides different models for BH growth and AGN feedback, many other physical details (e.g. models for star formation, stellar feedback) or different hydrodynamic schemes may cause more fundamental changes in basic galaxy properties. Such an investigation is, however, clearly beyond the scope of this work.

6 SUMMARY AND CONCLUSIONS

In this paper, we presented an improved implementation of the BH model originally introduced by Springel et al. (2005). We combined theoretical predictions of Churazov et al. 2005, Narayan & Yi 1995 and Gaspari et al. (2013) with observations from Russell et al. 2013, Mezcua & Prieto 2014, Davis & Laor 2011 and Chelouche 2013 in order to model the underlying sub-grid processes more realistically.

The new model includes a combination of mechanical outflow and radiation, which we both implemented as thermal feedback due to the inability of resolving sub-kpc scales, where jets provide the mechanical feedback. Both feedback processes are modelled as a function of the actual accretion rate with respect to the Eddington rate, which leads to a smooth transition between the outflow-dominated radio-mode and the radiation-dominated quasar-mode. In addition, our model includes a mass-dependent radiative efficiency to account for the observed spin of the BHs.

Furthermore, we distinguish between the hot and the cold gas component within the environment of the BHs and calculate the accretion rate for these two components separately. This allows us to model the Bondi accretion differently for the two phases, where we use two different boost factors (α = 10 for the hot and α = 100 for the cold gas) according to the results of small-scale simulations of Gaspari et al. (2013).

Besides that, free parameters of the model (like the various efficiencies) are now more strictly linked to values inferred from observations. Compared to the fiducial model, our new implementations predict a more realistic population of BHs and their host galaxies, when compared to fundamental observational constraints, in several aspects.

• The slope and normalization of the produced M_•–M_* relation are in much better agreement with observations over a larger range of galaxy masses and redshifts than in the fiducial model. In particular, these improvements are due to the faster BH growth at large redshifts and the lower BH masses at the massive end for redshifts below z ≈ 2.

• Our new feedback scheme is also able to efficiently suppress the late growth of massive black holes. Hence, the resulting present-day BH mass function provides an excellent match to the observed one.

• In the NFAM simulations, the equilibrium between gas cooling and AGN feedback within the galaxies is reached earlier. Consequently, star formation starts to be suppressed at earlier times. This leads to a better agreement with observed stellar mass functions than before. In particular, in the NFAM simulation there are much more quiescent galaxies at high redshifts than in the fiducial simulation, in which galaxies become quiescent far too late. However, some inconsistencies between observed and simulated SMFs for quiescent and star-forming galaxies remain.

• The baryon conversion efficiencies are more consistent with observations and abundance matching predictions than before, although they are still too high by a factor of 2–3 at very high stellar masses.

A comparison with other large cosmological simulations (e.g. Illustris, MassiveBlack-II) illustrates that the original BH model from Springel et al. (2005) needs to be extended to be able to reproduce observations. In particular, we find that

• our NFAM simulation successfully matches the observed M_•–M_* relation. As our fiducial model, the simulations from Sijacki et al. (2014) and Khandai et al. (2014) do not manage to entirely reproduce the observed slope. This may be due to the constant values adopted for their radiative efficiencies.

• In contrast to the MassiveBlack-II simulation, both our NFAM simulation and the Illustris simulation are able to reproduce the observed luminosity functions. We suggest that this might be due to the distinction between quasar-mode and radio-mode.

• our model predicts a lower high mass end of the BH mass functions than other simulations (i.e. Di Matteo et al. 2008, Sijacki et al. 2014), because the new AGN feedback model is more efficient in limiting BH growth at higher masses. Although all simulations are compatible with the upper limits of the BH mass function estimated from observations by Shankar (2013), our model is in excellent agreement with the observational data from Marconi et al. (2004), Shankar et al. (2004) and Shankar et al. (2009).

• We predict lower stellar masses than Genel et al. (2014) and Khandai et al. (2014). Since our new implementations do not change the SMFs at z = 0 significantly, we suggest that other physical processes like stellar feedback or cooling might be the reason for the differences. In addition, we find that improvements in the model for star formation and stellar feedback like in Schaye et al. (2015) might be necessary to better reproduce the observed shape of the SMFs.

Despite of the overall success of the NFAM model, open questions regarding the actual values of the feedback efficiencies remain. In contrast to the quasar-mode, the radiative efficiency in the radio-mode does not show clear trends in observations, which generally have large uncertainties, especially due to the difficulties in accurately determining the accretion rate. At high redshifts, the quasar luminosity function predicted by the simulations is quite insensitive to the choice of the radiative efficiency in the radio-mode. However, the best match between simulated and observed quasar luminosity functions – especially at low redshifts – is obtained when applying a random radiative efficiency to the simulated AGN in the radio-mode with no dependence on BH mass or actual accretion rate.

Studying the growth of BHs in more detail (i.e. for individual objects) provides evidence for a two phase process controlling the evolution of the accretion on to the BH and the associated feedback.

• As long as BHs have masses below the M_•–M_* relation, they grow mainly due to continuous gas accretion. This phase is primarily driven by cold gas accretion with an accretion rate that increases up to the Eddington limit. In this phase, AGN are observed as luminous X-ray sources. This means that the most luminous AGN are not necessarily driven by merger events as long as they are below the M_•–M_* relation.

• When the M_•–M_* relation is reached, gas cooling and AGN feedback are in equilibrium. Consequently, hot gas accretion begins to dominate. This means that the accretion rate, compared to the original implementation, is lowered since we correctly reduce the boost factor for the hot phase. In this phase, AGN feedback is mostly visible as radio jets. This low accretion phase can be disturbed by mergers or other processes driving cold gas into the centre of the galaxy. In a forthcoming study of the most luminous AGN in a simulation with higher resolution we will investigate in more detail whether those objects are mainly triggered by major mergers.

Regarding the latter point, more detailed studies are needed to better differentiate the AGN triggering mechanisms (as galaxy major and minor mergers) and their correlation with the BH accretion processes within a cosmological context. The next generation of simulations will also allow us to distinguish between morphological types of galaxies in more detail and thus, to investigate the connection between AGN luminosities and the host galaxy morphologies, hopefully shedding more light on the main trigger mechanisms for AGN activity in different redshift and luminosity regimes. Such future simulations will also help to understand the dependence of the AGN driving mechanisms on the large-scale environment.

In addition, we plan to further improve the current implementations by taking the angular momentum of the accreted material into account, which in turn would allow us to better model the direction of the feedback. This would especially have an important effect on the spatial distribution of the feedback energy in the surroundings of the AGN. Indeed, current BH accretion and feedback models are purely empirically motivated and have the major drawback that they do not capture the underlying small-scale physical processes, which is, within the framework of large-scale cosmological simulation, currently not feasible due to limited computational power. Nevertheless, despite of the rather crude approximations, the BH model, in particular with our new modifications, seems to capture the essence of how BHs grow and how feedback affects the host galaxies in reality.

Future observations will improve our understanding of the different accretion modes and their relation to the multiphase nature of the ICM/IGM. In particular, studies of Seyfert galaxies (Mezcua et al. in preparation) will allow an investigation of the role of warm H₂ gas (with temperatures of ∼10³ K). In combination with X-ray observations, this will shed more light on the complicated interplay between the various accretion modes of AGN.

We thank the referee for a careful and constructive reading of our paper. We also thank Alexander Beck, Veronica Biffi, Andreas Burkert, Massimo Gaspari, Mar Mezcua, David Schlachtberger, Francesco Shankar and Adelheid Teklu for many fruitful discussions. Particularly, we thank Madhura Killedar for carefully reading the text and editing. Furthermore, we would like to thank Shane W. Davis and Helen Russell for providing us with observational data and Veronica Biffi and Umberto Maio for providing us with the data for Fig. 12.

We are especially grateful for the support by M. Petkova through the Computational Center for Particle and Astrophysics (C2PAP). Computations have been performed at the ‘Leibniz-Rechenzentrum’ with CPU time assigned to the Project ‘pr86re’.

This research was supported by the DFG Cluster of Excellence ‘Origin and structure of the Universe’ and the SFB-Tansregio TR33 ‘The Dark Universe’.

MH acknowledges financial support from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007- 2013)/ERC grant agreement no. 202781 and support from the European Research Council via an Advanced Grant under grant agreement no. 321323NEOGAL.

MAP thanks the hospitality of the Universitäts-Sternwarte München and the Max-Planck-Institut für extraterrestrische Physik where she stayed as visitor.

1

For the data from Russell et al. (2013), the dynamical masses from McConnell & Ma (2013) were taken if available. If not, the same masses were taken which Russell et al. (2013) used to calculate L_Edd.

2

www.magneticum.org

3

Note that this value depends on the resolution, because at lower resolutions the feedback energy is spread further away from the BH. Hence, for our simulations, this value is comparatively high.

4

The two outliers (black and red diamond with M_• ≈ 2 × 10⁸M_*) are due to temporary attributions to different haloes.

5

We excluded the outlier (black diamond on the left with M_• ≈ 2 × 10⁸M_*).

6

For the observations we computed M_200c out of M_500c using the NFW profile.

7

The amount of BHs does not change, because we use the same simulations for all different feedback models. Consequently, lower number densities of AGN with L > 10⁴² erg s⁻¹ are equivalent to higher number densities for fainter AGN.