Supermassive black holes in cosmological simulations
Observations suggest that almost every galaxy hosts a supermassive ($10^8 \div 10^{10}$ M$_{\odot}$) black hole (SMBH) in its innermost regions (e.g., [231] ,[232] ,[233] ,[234] ,[235] ,[236] ,[237] ,[238] ,[239] ). A fraction of these BHs exhibits ongoing activity and is called AGN (active galactic nuclei). Evidence of past and ongoing activity is observed for instance in the X-ray images of elliptical galaxies and galaxy groups and clusters, where AGN imprints often appear as depressions and ripples. Among many, a well-known example is represented by the composite (X-ray, radio, and visual) image of the MS 0735+7421 galaxy cluster (e.g. [240] ): giant X-ray cavities are filled with radio emission, and surrounded by a cocoon shock clearly visible in the Chandra image as an elliptical edge.
Before discussing AGN feedback (see §3.16), we are now recalling the main features of BHs in cosmological simulations, focusing on BH seeding and repositioning, BH-BH mergers, and AGN feeding.
The majority of cosmological simulations of galaxy and galaxy cluster formation that also include BHs and the associated AGN feedback are based on the seminal BH model by [241] ,[242] , or follow its spirit. Therefore, this model offers the base to discuss a few key aspects of the treatment of BHs in cosmological simulations, and to highlight recent improvements.
In cosmological simulations, BHs are usually described as sink particles: these are collisionless particles that can absorb neighbour elements (originally implemented to allow the removal of surrounding particles in dense regions, e.g. [243] ) and that have fundamental properties like the accretion rate which can be linked directly to observables.
BH particles are introduced in massive haloes at relatively high-redshift in cosmological simulations, and they are then allowed to grow and increase their initial or {seed} mass. The aforementioned, commonly quoted mass $M_{\bullet}$ is the theoretical mass of the BH, modelled at the sub-resolution level, as opposite to its dynamical mass, i.e. the actual gravitational mass of the BH particle. As we are still lacking a solid understanding of the formation of first SMBHs (see e.g., [244] ,[245] ,[246] ,[247] ,[248] for possible pathways) and the resolution needed to take the physics of any seed formation scenario into account, BHs are first inserted according to seeding prescriptions (see also [249] ).
Seeding prescriptions usually assume that new BHs (massive seeds of $\sim 10^4 \div 10^6$ M$_{\odot}$) are introduced in haloes which meet some criteria and do not have already BHs. New BHs are seeded if: $(i)$ the halo mass -- commonly estimated by means of a FOF (Friend-Of-Friend, [250] ) algorithm -- is larger than a threshold (e.g., [40] ,[41] ,[139] ,[6] ); $(ii)$ the stellar mass exceeds a given value (e.g., [44] ); $(iii)$ the stellar mass and gas to stellar mass fraction are larger than given thresholds; $(iv)$ based on gas properties, i.e. gas density, velocity dispersion and/or metallicity (e.g., [251] ,[252] ,[253] ,[142] ). BH mass at seeding can be either constant (e.g., [40] ,[41] ,[139] ,[6] ,[44] ) or scaled according to e.g., the M$_{BH}$/$\sigma$ or M$_{BH}$/M$_{\ast}$ scaling relations (e.g., [138] ). Additional details involve the type of particles on which the FOF algorithm is performed, whcih can be for instance DM particles only (e.g., [254] ) or stellar particles only (e.g., [138] ). Besides, BH can be seeded at the position of the densest gas particle in the halo (e.g., [254] ), of the star particle with the largest binding energy (e.g., [138] ), or of the star particle closest to the centre of mass of the structure (e.g., [44] ).
In cosmological simulations, each BH undergoes a specific evolution as an individual particle, at variance with the coarse-grained representation of collisionless fluids through particles in N-body simulations. As a result, the dynamics of BHs fails to be accurately captured and numerical artefacts can affect BH motion (see e.g., [255] ). Spurious displacements of BHs, which often occur due to scattering between particles in high-density environments and numerical heating, have dramatic consequences, such as e.g., artificial presence of wandering BHs, incorrect description of BH-BH mergers, BHs producing feedback off-centre with respect to their galaxy hosts. To avoid these artefacts, BH dynamics is commonly controlled through different approaches: re-positioning via pinning on e.g. minimum potential (e.g., [40] ,[41] ,[139] ,[44] ), boosted dynamical mass (e.g., [144] ), boosted dynamical mass and dynamical friction (e.g., [138] ,[145] ). We refer the reader to [254] ,[256] ,[35] ,[257] ,[143] ,[258] ,[249] ,[259] ,[260] for details and recent improvements.
BHs grow because of gas accretion and mergers with other BHs. As for the latter channel, BHs are expected to merge when their host galaxies and their haloes merge to form a single structure. The commonly pursued approach [242] ,[241] in cosmological simulations assumes that two BH particles merge quickly if their distance approaches the spatial resolution of the simulation (or a small multiple of it). The force resolution set by the gravitational softening determines indeed the minimum scale above which gravitational interactions can be properly followed. Additional conditions involving e.g. the merging BH relative speed have been implemented [254] . The two BHs are eventually merged into a single BH particle, with their masses combined.
As for AGN feeding, gas accretion onto a BH of mass $M_{\bullet}$ is calculated according to the Bondi formula [261] ,[262] ,[263] , multiplied by a so-called boost factor $\alpha$: $$ \dot{M}_{B} = \frac{4 \pi \, \alpha \, G^2\, M_{\bullet}^2 \, \langle \rho \rangle}{(\langle c_s\rangle^2 +\langle v\rangle ^2)^{3/2}} \,\, .
$$ Here, $G$ is the gravitational constant, $\langle\rho\rangle$, $\langle v\rangle$, and $\langle c_s\rangle$ are mean values at the scale resolved by the hydrodynamical simulation: for example, they are computed using kernel weighted estimates in the case of SPH. The BH accretion rate $\dot{M}_{B}$ is commonly capped to the Eddington acccretion rate $\dot{M}_{Edd}$ [Fn: $\dot{M}_{Edd}=(4\pi\:G\;m_p\; M_{\bullet})/(\sigma_T\;c\;\epsilon_r)$, with $m_p$ the proton mass, $\sigma_T$ the Thompson cross section, $c$ the speed of light and $\epsilon_r$ the radiative efficiency, typically assumed to be $\approx0.1$.] (but see [264] ,[44] for examples of $\dot{M}_{\bullet}$ which breaches the Eddington limit).
The boost factor $\alpha$ has been originally introduced [241] to account for the limited resolution in simulations, which leads to smaller densities and larger temperatures near the BH (and thus to an underestimate of $\dot{M}_{B}$). A typical value is $\alpha =100$. Several studies adapt the BH model by using a boost factor which depends on resolution [265] ,[266] , density [267] , or pressure [268] . Other simulations instead limit the BH accretion rate by taking into account the angular momentum of accreting gas [269] ,[270] . High-resolution simulations of BH accretion on sub-kpc scales [271] found that a boost factor of order of 100 is suitable when including cooling and turbulence, while pure adiabatic accretion suggests boost factors smaller by an order of magnitude. Hence, advanced models distinguish between hot and cold gas accretion and use different boost factors for the two components [272] , or can even get rid of fudge factors (e.g., [139] ,[270] ).
By exploting the BH accretion rate $\dot{M}_{\bullet}$, the bolometric luminosity $L_{\mathrm bol}$ can be associated to each BH in the simulation. A common assumption consists in following [273] to distinguish between high and low BH accretion state in terms of the Eddington ratio $f_{Edd} = \dot{M}_{\bullet}/ \dot{M}_{Edd}$ (i.e. the ratio between the BH and the Eddington accretion rates) and compute: $$ L_{\mathrm bol} = \left\{ \begin{array}{ll} 10 \, ( \epsilon_r \; c)^2 \, \dot{M}_{\bullet} & f_{\mathrm Edd}>0.1 \\ (10 \, \epsilon_r\; c)^2 \, f_{\mathrm Edd}\; \dot{M}_{\bullet}\;\; & f_{\mathrm Edd}\leq0.1 \end{array} \right. $$ (see [274] ,[275] ). In contrast to the original model [242] ,[241] , modern implementations often correct the accretion rate $\dot{M}_{\bullet}$ of the BH by a factor $(1-\epsilon_r)$. As a consequence, $\,L_{\mathrm bol} = \epsilon_r / (1-\epsilon_r) \, \dot{M}_{\bullet} \,c^2 \,$, and the energy radiated away during the accretion process is taken into account. Different prescriptions for the accretion rate, in combination with the change of the ISM/IGM properties produced by the combined action of all the sub-grid processes, typically leads to significant variations in the predicted evolution of the AGN luminosity function among the various state-of-the-art simulations. This is summarized in Fig. 11, where we can appreciate how predictions from simulations can strongly deviate from the observed AGN luminosity function across cosmic time (see also discussion in [275] ).

评论