Spin evolution and feedback of supermassive black holes in cosmological simulations

TL;DR

This paper introduces a new model for supermassive black hole spin evolution in cosmological simulations, accounting for various accretion and merger processes, and examines its impact on galaxy and black hole properties.

Contribution

It presents a novel sub-grid model for black hole spin evolution integrated into the AREPO code, enhancing the physical realism of galaxy formation simulations.

Findings

High-spin BHs at lower masses due to coherent accretion

Lower-spin BHs at higher masses due to mergers and chaotic accretion

Improved galaxy-BH population agreement and galaxy morphology predictions

Abstract

It is well established that the properties of supermassive black holes and their host galaxies are correlated through scaling relations. While hydrodynamical cosmological simulations have begun to account for the co-evolution of BHs and galaxies, they typically have neglected the BH spin, even though it may play an important role in modulating the growth and feedback of BHs. Here we introduce a new sub-grid model for the BH spin evolution in the moving-mesh code {\small AREPO} in order to improve the physical faithfulness of the BH modelling in galaxy formation simulations. We account for several different channels of spin evolution, in particular gas accretion through a Shakura-Sunyaev $\alpha$-disc, chaotic accretion, and BH mergers. For BH feedback, we extend the IllustrisTNG model, which considers two different BH feedback modes, a thermal quasar mode for high accretion states and a kinetic mode for low Eddington ratios, with a self-consistent accounting of spin-dependent radiative efficiencies and thus feedback strength. We find that BHs with mass $M_{\rm{bh}}\lesssim 10^{8}\, {\rm M}_\odot$ reach high spin values as they typically evolve in the coherent gas accretion regime. On the other hand, BHs with mass $M_{\rm{bh}}\gtrsim 10^{8}\, {\rm M}_\odot$ have lower spins as BH mergers become more frequent, and their accretion discs fragment due to self-gravity, inducing chaotic accretion. We also explore the hypothesis that the transition between the quasar and kinetic feedback modes is mediated by the accretion mode of the BH disc itself, i.e.~the kinetic feedback mode is activated when the disc enters the self-gravity regime. We find excellent agreement between the galaxy and BH populations for this approach and the fiducial TNG model with no spin evolution. Furthermore, our new approach alleviates a tension in the galaxy morphology -- colour relation of the original TNG model.

Tables

Table 1. Table 1: Overview of our primary cosmologial simulation models. For all the simulations, we employ the physics model of IllustrisTNG with some modifications with respect to our new spin model and treatment of the BH feedback switch. Simulation 1 is used as a control for comparison purposes and corresponds to the unmodified IllustrisTNG model.

Simulation	Spin model	$k$	Feedback switch
1	No	-	Fiducial
2	Yes	0	Fiducial
3	Yes	2	Fiducial
4	Yes	10	Fiducial
5	Yes	0	Self-gravity

Equations

$$J_{\text{mbox}bh}=|a|\frac{G{M}_{\text{mbox}bh}^2}{c},$$

$${\hat{r}}_{\text{mbox}isco}=\frac{r_{\text{mbox}isco}}{R_g}=3+Z_2\pm [(3-Z_1)(3+Z_1+2Z_2){]}^{1/2},$$

$$R_{\text{mbox}g}=\frac{R_{\text{mbox}Schw}}{2}=\frac{G{M}_{\text{mbox}bh}}{c}.$$

$$Z_1$$

$$Z_2$$

$$\mathrm{d}{M}_{\text{mbox}bh}=\frac{{\tilde{e}}_{\text{mbox}isco}}{c^2}\mathrm{d}{M}_0,\mathrm{d}{J}_{\text{mbox}bh}={\tilde{l}}_{\text{mbox}isco}\mathrm{d}{M}_0.$$

$$\frac{\mathrm{d}a}{\mathrm{d}\ln M_{\text{mbox}bh}}=\frac{1}{M_{\text{mbox}bh}}\frac{c^3}{G}\frac{{\tilde{l}}_{\text{mbox}isco}}{{\tilde{e}}_{\text{mbox}isco}}-2a.$$

$$a^f=\frac{1}{3}{\hat{r}}_{\text{mbox}isco}^{1/2}\frac{M_{\text{mbox}bh}}{M_{\text{mbox}bh}^f}\left[4-{\left\{3{\hat{r}}_{\text{mbox}isco}{\left(\frac{M_{\text{mbox}bh}}{M_{\text{mbox}bh}^f}\right)}^2-2\right\}}^{1/2}\right],$$

$$M_{\text{mbox}d}=\frac{M_{\text{mbox}bh}^f-M_{\text{mbox}bh}}{1-{ϵ}_r}.$$

$${ϵ}_r\equiv 1-\sqrt{1-\frac{2}{3}\frac{1}{{\hat{r}}_{\text{mbox}isco}(a)}}.$$

$$\frac{\partial \mathbf{L}}{\partial t}={\mathbf{\omega }}_{\text{mbox}prec}\times \mathbf{L},$$

$$\frac{R_{\text{mbox}warp}}{R_{\text{mbox}Schw}}=3.6\times 10^3{a}^{5/8}{\left(\frac{M_{\text{mbox}bh}}{10^8{\mathrm{M}}_{\odot }}\right)}^{1/8}{f}^{-1/4}{\left(\frac{{\nu }_2}{{\nu }_1}\right)}^{-5/8}{\alpha }^{-1/2},$$

$$t_{\text{mbox}acc}=\frac{R_{\text{mbox}warp}^2}{{\nu }_1}=3\times 10^6{a}^{7/8}{\left(\frac{M_{\text{mbox}bh}}{10^8{\mathrm{M}}_{\odot }}\right)}^{11/8}$$

$$\times {\lambda }^{-3/4}{\left(\frac{{\nu }_2}{{\nu }_1}\right)}^{-7/8}{\alpha }^{-3/2}\text{mbox}yr.$$

$${\mathbf{J}}_{\text{mbox}tot}={\mathbf{J}}_{\text{mbox}bh}+{\mathbf{J}}_{\text{mbox}d}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{J}_{\text{mbox}bh}^2=0,\frac{\mathrm{d}}{\mathrm{d}t}{J}_{\text{mbox}d}\le 0,\frac{\mathrm{d}}{\mathrm{d}t}\cos {\theta }_t\ge 0,$$

$$\cos \theta <-\frac{J_{\text{mbox}d}}{2J_{\text{mbox}bh}},$$

$$J_{\text{mbox}d}\approx M_d{\Omega }_k(R_{\text{mbox}warp})R_{\text{mbox}warp}^2=M_d(G{M}_{\text{mbox}bh}{R}_{\text{mbox}warp}{)}^{1/2}.$$

$$\frac{J_{\text{mbox}d}}{2J_{\text{mbox}bh}}\approx \frac{M_{\text{mbox}d}}{a{M}_{\text{mbox}bh}}{\left(\frac{R_{\text{mbox}warp}}{R_{\text{mbox}Schw}}\right)}^{1/2}.$$

$$Q\equiv \frac{c_s{\Omega }_k}{\pi G{\Sigma }_{\text{mbox}d}},$$

$$\frac{R_{\text{mbox}sg}}{R_{\text{mbox}Schw}}=1.5\times 10^3{ϵ}^{8/27}{\left(\frac{M_{\text{mbox}bh}}{10^8\text{mbox}{M}_{\odot }}\right)}^{-26/27}{f}^{-8/27}{\alpha }^{14/27}.$$

$$M_{\text{mbox}sg}=2.13\times 10^5{ϵ}^{-5/27}{\left(\frac{M_{\text{mbox}bh}}{10^8\text{mbox}{M}_{\odot }}\right)}^{23/27}{f}^{5/27}{\alpha }^{-2/27}\text{mbox}{M}_{\odot }.$$

$$f(\theta |k)=\frac{e^{k\cos \theta }}{2\pi I_0(k)},$$

$${\mathbf{a}}^f=\frac{1}{(1+q{)}^2}({\mathbf{a}}_1+{\mathbf{a}}_2{q}^2+\mathbf{\ell }q),$$

$$\ell$$

$${\mathbf{v}}_{\text{mbox}recoil}$$

$$v_m$$

$$v_{\perp }$$

$$v_{\parallel }$$

$${\dot{M}}_{\text{mbox}Bondi}=\frac{4\pi G^2{M}_{\text{mbox}bh}^2\rho }{c_s^3},$$

$$\dot{M}=\min ({\dot{M}}_{\text{mbox}Bondi},{\dot{M}}_{\text{mbox}Edd}),$$

$$\Delta {\dot{E}}_{\text{mbox}high}={ϵ}_{\text{mbox}f,high}{ϵ}_r\dot{M}{c}^2,$$

$$\Delta {\dot{E}}_{\text{mbox}low}={ϵ}_{\text{mbox}f,kin}\dot{M}{c}^2,$$

Full text

††pagerange: Spin evolution and feedback of supermassive black holes in cosmological simulations–References††pubyear: 2018

1 Introduction

It is remarkable finding that most, quite possibly all, galaxies host a supermassive black hole (BH) at their centres (e.g. Greene et al., 2010). The origin of these black holes (i.e. the “seeds” that subsequently grew) in the high redshift universe is still unclear, but it appears relatively certain that much of the mass of the black hole population was accumulated through gas accretion (Sołtan, 1982). The accretion processes themselves are accompanied by a substantial release of energy, which is most dramatically seen in the emission of powerful quasars (Lynden-Bell, 1969).

The relatively tight scaling relations between galaxy properties and the masses of the central supermassive BH (Magorrian et al., 1998; Ferrarese & Merritt, 2000; Häring & Rix, 2004) suggest that the BHs are not merely spectators in the galaxy formation process, but rather influence it in decisive ways. This view has motivated the incorporation of black hole evolution and feedback processes in modern theories of galaxy formation (Kauffmann & Haehnelt, 2000). In particular, many simulation models explain the quenching transition of galaxies from the star forming blue cloud to the red sequence of elliptical galaxies as a result of BH feedback (Di Matteo et al., 2005; Springel et al., 2005; Sijacki et al., 2007). In large clusters of galaxies, the BHs can furthermore offset radiative cooling losses through sporadically going through phases of active galactic nuclei (AGN) activity, thereby preventing excessive cooling flows. Another importance of this AGN feedback lies in its ability to expel baryons even from deep potential wells of galaxy groups and clusters, helping to understand the observed trends of the baryon content in these objects (McCarthy et al., 2007; Puchwein et al., 2008).

The recent discovery of gravitational waves from stellar-mass black holes (Abbott et al., 2016) has additionally fuelled the interest in the merger processes of supermassive black holes, which are inevitable events during hierarchical structure formation. In the future, coalescing BH binaries in the relevant mass range should become detectable by pulsar timing arrays and the Laser Interferometer Antenna (LISA). These unprecedented experiments will shed light on the BH assembly process, but they also require theoretical predictions for the expected merger rates and a characterisation of the precursors of BH binaries. The BH mergers will emit copious amounts of gravitational waves, leading to a recoil of the merger remnant. As a result, they may be even ejected from the host galaxy in certain situations. In any case, the timescale of return of the remnant BH to the centre of the galaxy is important as during this phase the regulation of a quasi-static hot atmosphere in the halo may temporarily interrupted, allowing star formation to resume and potentially impact the evolution of massive galaxies.

Supermassive black holes are therefore a key element of any modern theory of galaxy formation and evolution, and they are now routinely incorporated both in semi-analytic (e.g. Croton et al., 2006; Henriques et al., 2015) and hydrodynamical simulation models (e.g. Vogelsberger et al., 2014; Schaye et al., 2015). Despite this central importance of supermassive BHs, the physical processes governing their evolution, especially with respect to the gas accretion and the associated feedback processes, are only poorly understood. The modelling in cosmological hydrodynamic simulations is hence very sketchy, and typically encapsulated in heuristic sub-grid models. Commonly, black holes are simply treated as sink particles of a given mass, and all modelling of the physics neglects the one other important parameter characterising the black holes, namely their angular momentum. This is a significant limitation, because the BH spin affects the radiative efficiency of accretion processes decisively. In fact, while it is believed that typically of order 5% of the accreted rest mass energy of a non-spinning BH is released as radiation, this can go up to about 40% for highly spinning black holes (Blandford & Znajek, 1977; Tchekhovskoy et al., 2011). In addition, the spin of merging BHs has a strong influence on the gravitational wave emission, and thus also on the expected recoil velocity of the merger remnant (Campanelli et al., 2007).

For all these reasons, modelling BH spins appears to be a promising extension of current treatments of black hole physics. Corresponding simulations would then not only yield predictions for the BH mass function and its evolution over time, for the lifetimes and activity patterns of quasars, their clustering, and for their overall influence on the galaxy formation process, they would also be able to make predictions for the statistics of BH spin, including for the typical spin configurations expected in supermassive BH mergers. This in turn would inform attempts to forecast gravitational wave signals and their frequency, and elucidate the importance of spin for the energy release through gas accretion. It is therefore timely to go beyond the first generation of BH models and implement more detailed physical models for BH spin in the current cosmological codes for galaxy formation.

At first sight, it may seem conceptually simple to relate the BH growth to the specific angular momentum content of the accreted gas, but the huge range of scales that are in reality still unresolved in cosmological simulations can not simply be ignored easily. In fact, the accretion flows that actually develop close in to the BH can give rise to angular momentum transport, in addition to an important coupling between a preexisting BH spin and the surrounding gas, causing phenomena such as Lense-Thirring precession. This necessitates the construction of explicit sub-grid models for the BH spin evolution which are guided by analytic theory, and which respond to the larger-scale boundary conditions set by the cosmological simulation. A few such models have already been constructed. Fanidakis et al. (2011) pioneered such treatments in the context of semi-analytic modelling of galaxy formation. Dubois et al. (2014) study the BH spin evolution in hydrodynamical cosmological simulations in a post-processing fashion. Recently, and independently from our work, Fiacconi et al. (2018) presented a first effort using the moving-mesh code AREPO, which we also employ here. Our model is similar in spirit to theirs, but differs in a number of details. Importantly, we consider the back-reaction on the feedback efficiency and couple the spins to the black hole feedback mode. In particular, we conjecture to relate the transition in feedback modes, introduced in the IllustrisTNG black hole treatment in an ad-hoc fashion by Weinberger et al. (2017), to the spin evolution, and test the impact this has.

This paper is structured as follows. In Section 2, we begin by reviewing the basic problem of how BH spin evolves through gas accretion and BH binary coalescence. In Section 3, we detail our black hole growth and feedback model, as well as its coupling to the spin evolution. We then discuss our numerical implementation of these processes and our simulation set in Section 4. Our results for the statistics of BH spin in cosmological simulations are presented in Section 5, while results for the impact of our spin modelling through a novel feedback switch on the Section 6. Finally, we summarise our findings and give our conclusions in Section 7.

2 Spin evolution of Black Holes

In this section, we present a model for the evolution of the BH spin. Besides its mass, the spin of a BH is an important property as it affects the efficiency with which a BH can process parts of the matter in an accreted disc into radiation; it impacts the strength and emission direction of radio jets; it possibly regulates the radio-loudness of the host galaxy; and it determines both the detailed gravitational wave emission signal and the distribution of recoil velocities in BH-BH mergers, which is important for forecasting these events and the subsequent evolution of the BH remnant.

The spin is simply defined as the BH angular momentum with respect to the centre-of-mass frame, and is commonly expressed as:

[TABLE]

where $a$ is the dimensionless spin parameter ($0\leq|a|\leq 1$), $M_{\mbox{\scriptsize{bh}}}$ the black hole mass, $G$ the gravitational constant and $c$ the vacuum speed of light. This parametrisation results from assuming a Kerr metric for rotating bodies, in which the maximum allowed spin is $GM_{\mbox{\scriptsize{bh}}}^{2}/c$. In what follows, we assume two main channels through which BHs can be spun up or down, namely gas accretion and BH mergers.

2.1 Gas accretion

BHs in galactic centres process matter into radiation through a surrounding accretion disc that is fed by cold gas inflows from galactic scales. If a BH is accreting above $1\%$ of its Eddington rate, we assume that a classical thin, optically thick Shakura-Sunyaev $\alpha$-disc (Shakura & Sunyaev, 1973) settles around it (Fanidakis et al., 2011). On the other hand, for accretion rates below this value, the $\alpha$-disc model is no longer appropriate and advection dominated disc models (ADAF) would have to be considered. However, given that in this case the gas does not carry enough angular momentum to modify the spin of the BH significantly, we simply turn off the gas channel under these conditions.

Once an accretion disc is formed, the orbiting gas loses its angular momentum through viscous torques caused by magnetic fields (Lynden-Bell, 1969) and radiates away the binding gravitational energy, moving thus inwards and reaching the radius of the innermost stable circular orbit (isco) of the BH. After this point, further angular momentum and energy losses are no longer required as the gas will then be accreted on a very short timescale.

The radius of the isco can be expressed in terms of the BH spin as (Bardeen et al., 1972):

[TABLE]

where the negative sign is taken if the disc is counter-rotating with respect to the BH spin (i.e. $a<0$). The positive sign corresponds to the case of co-rotation. For the sake of simplicity, a normalisation with the gravitational radius $R_{g}$ is used here, which is defined as half of the Schwarzschild radius:

[TABLE]

The quantities $Z_{1}$ and $Z_{2}$ depend on the BH spin, and are defined as:

[TABLE]

Once a parcel of gas of mass ${\rm d}M_{0}$ reaches the isco, it will be eventually accreted and its energy and angular momentum per unit mass, $\tilde{e}_{\mbox{\scriptsize{isco}}}$ and $\tilde{l}_{\mbox{\scriptsize{isco}}}$, will add to the total BH mass and angular momentum, respectively:

[TABLE]

Combining these expressions with equation (1) yields the following differential equation:

[TABLE]

Integrating this equation, Bardeen (1970) obtained a solution for the evolution of the spin magnitude:

[TABLE]

where $a^{f}$ and $M_{\mbox{\scriptsize{bh}}}^{f}$ correspond to the final spin and BH mass after an accretion episode of mass

[TABLE]

Here, $\epsilon_{r}$ is the radiative efficiency of the accretion disc, which measures the fraction of gravitational binding energy released by gas spiralling in towards the BH. In the case of a Kerr BH, the radiative efficiency depends on the spin as follows (Novikov & Thorne, 1973):

[TABLE]

Note that equation (8) can be used to update the BH spin only if $M^{f}_{\mbox{\scriptsize{bh}}}/M_{\mbox{\scriptsize{bh}}}\leq\hat{r}_{\mbox{\scriptsize{isco}}}^{1/2}$. In the case of $M^{f}_{\mbox{\scriptsize{bh}}}/M_{\mbox{\scriptsize{bh}}}>\hat{r}_{\mbox{\scriptsize{isco}}}^{1/2}$, a single accretion episode is able to spin up the BH to its maximum value $a=1$; however, we cap the maximum spin to $a=0.998$ to account for the absorption of photons with angular momentum opposite to the BH spin (Thorne, 1974).

2.1.1 Accreting gas in misaligned discs

The previous analysis led us to a formula to compute the evolution of the BH spin once a mass $M_{d}$ has fallen in towards the accretion disc. However, this analysis provides no information about the accretion process and on how much mass is ultimately accreted, i.e. we do not know the value of $M_{d}$. In order to estimate this, we assume a general case in which an accretion disc settles around the BH and lies initially misaligned with respect to the BH equatorial plane. Under this configuration, the orbiting gas experiences a torque caused by the Lense-Thirring effect of the spinning BH, which is given by:

[TABLE]

where $\boldsymbol{\omega}_{\mbox{\scriptsize{prec}}}=(2G/c)\boldsymbol{J}_{\mbox{\scriptsize{bh}}}/R^{3}$ is the precession rate, $\boldsymbol{L}$ is the angular momentum per unit area of the disc and $R$ the distance from the BH (Pringle, 1992). This torque causes the gas orbits to become unstable and evolve into lower energy states, which in turn makes the plane of the disc precess around the BH spin axis at a frequency ${\omega}_{\mbox{\scriptsize{prec}}}$. If the viscosity of the gas is sufficiently high, the innermost part of the disc aligns with the equatorial plane of the spinning BH, giving rise to a warped disc.

In the $\alpha$-disc solution, the kinematic viscosity of the gas is split into two components, $\nu_{1}$ and $\nu_{2}$, which act along velocity gradients parallel and perpendicular to the disc plane, respectively. This allows us to introduce two characteristic timescales: an accretion timescale $t_{\mbox{\scriptsize{acc}}}\equiv R^{2}/\nu_{1}$, during which gas drifts inwards and angular momentum is transported outwards through radial viscous dissipation; and a warp timescale $t_{\mbox{\scriptsize{warp}}}\equiv R^{2}/\nu_{2}$, during which any vertical disturbance propagates across the disc, e.g. the warp.

Introducing the precession timescale $t_{\mbox{\scriptsize{prec}}}\equiv 2\pi/\omega_{\mbox{\scriptsize{prec}}}$, it is possible to estimate the size of the warped region by imposing the condition $t_{\mbox{\scriptsize{prec}}}(R)\lesssim t_{\mbox{\scriptsize{warp}}}(R)$, i.e. if restoring viscous forces are overcome by frame-dragging torques the disc cannot retain its original shape. Using the following parametrisation of the radial viscosity $\nu_{1}=\alpha H^{2}\Omega_{k}$, where $H$ is the disc semi-thickness and $\Omega_{k}$ the Kepler frequency, the following expression for the warp radius $R_{\mbox{\scriptsize{warp}}}$ is obtained (Volonteri et al., 2007):

[TABLE]

where $f\equiv L/L_{\mbox{\scriptsize{Edd}}}$ is the Eddington ratio, $L=\epsilon_{r}\dot{M}c^{2}$ and $L_{\mbox{\scriptsize{Edd}}}=4\pi GM_{\mbox{\scriptsize{bh}}}m_{p}/\kappa$ are the disc and Eddington luminosities, respectively, $\dot{M}$ is the total accretion rate of the system, $m_{p}$ the proton mass, and $\kappa\approx 0.3\ \mbox{cm}^{2}\mbox{g}^{-1}$ denotes the electron scattering opacity.

The warp radius is an important quantity as it provides a characteristic length scale of the accretion disc. Besides, only the material inside it can effectively transfer angular momentum to the BH (Volonteri et al., 2007). The mass of the warped region can be readily estimated as $M_{\mbox{\scriptsize{warp}}}=\dot{M}\,t_{\mbox{\scriptsize{acc}}}(R_{\mbox{\scriptsize{warp}}})$, where:

[TABLE]

One of the main uncertainties of the $\alpha$-disc solution is the relation between the two viscosities $\nu_{1}$ and $\nu_{2}$, which ultimately determines which process, warp alignment or mass accretion, occurs in a shorter timescale. In the thin disc regime, i.e. $H/R<\alpha\ll 1$, the condition $\nu_{2}/\nu_{1}\approx 1/\alpha^{2}$ has to be fulfilled (Papaloizou & Pringle, 1983). We use the relation $\nu_{2}/\nu_{1}=2(1+7\alpha)/(\alpha^{2}(4+\alpha^{2}))$ proposed by Ogilvie (1999) that satisfies this condition. Under these assumptions, we obtain that $t_{\mbox{\scriptsize{warp}}}<t_{\mbox{\scriptsize{acc}}}$, which means that the innermost disc aligns with the equatorial plane of the BH quicker than the timescale of mass and angular momentum accretion. This allows us to use $M_{\mbox{\scriptsize{warp}}}$ as a proxy of the accreted mass of a single accretion episode, i.e. $M_{\mbox{\scriptsize{d}}}\approx M_{\mbox{\scriptsize{warp}}}$. Finally, using equation (9), we calculate the final BH mass $M_{\mbox{\scriptsize{bh}}}^{f}$ needed for evolving the BH spin through equation (8).

2.1.2 Alignment/Counter-alignment condition

We have previously concluded that before the content of the innermost disc is accreted, it settles on the equatorial plane of the BH. Additionally, King et al. (2005) established that the warped region, once settled, can align or anti-aligned with respect to the BH spin depending on the initial configuration of the disc angular momentum and the BH spin. In order to highlight the importance of determining the sense of rotation of the gas at the isco, we analyse equation (2). For a disc counter-rotating with respect to a maximum spinning BH, the radius of the isco is $\hat{r}_{\mbox{\scriptsize{isco}}}=9$. In the case of co-rotation, $\hat{r}_{\mbox{\scriptsize{isco}}}=1$. This means that counter-rotating gas needs to lose less energy and angular momentum than co-rotating gas before being ultimately accreted and, therefore, adds up more to the BH mass and spin. On the other hand, co-rotating gas also increases the energy output of the BH, favouring feedback processes.

In what follows, we inspect what condition determines alignment/anti-alignment. We denote the angular momentum of the innermost warped disc as $\boldsymbol{J}_{d}$. The total angular momentum is then given by:

[TABLE]

This vector is constant in magnitude and orientation during the whole accretion episode due to conservation of angular momentum in the BH-disc system. In Figure 1 we illustrate the relative orientations of the vectors for different phases of an accretion episode. King et al. (2005) derived a set of equations that describes the evolution of $\boldsymbol{J}_{\mbox{\scriptsize{bh}}}$ and $\boldsymbol{J}_{\mbox{\scriptsize{d}}}$ during the precession phase:

[TABLE]

where $\theta_{t}$ is the angle subtended between $\boldsymbol{J}_{\mbox{\scriptsize{tot}}}$ and $\boldsymbol{J}_{\mbox{\scriptsize{bh}}}$. This shows that during this phase, the magnitude of the BH spin remains constant, while the disc angular momentum is continuously decreasing. The direction of the BH spin eventually aligns with the total angular momentum, while the disc will either counter-align or align depending on whether or not the following condition is satisfied (King et al., 2005):

[TABLE]

with $\theta$ the angle subtended between $\boldsymbol{J}_{\mbox{\scriptsize{bh}}}$ and $\boldsymbol{J}_{\mbox{\scriptsize{d}}}$.

Considering that the angular momentum of the innermost disc is advected inwards, passing trough $R_{\mbox{\scriptsize{warp}}}$, its magnitude can readily be computed as (King et al., 2008):

[TABLE]

We have assumed that the gas has circularised on its way to the centre. This leads to the following expression for the ratio of disc to BH angular momentum:

[TABLE]

By evaluating this expression at the beginning of every accretion episode, we can determine whether alignment or counter-alignment occurs before the disc is accreted.

The last aspect to be determined is the orientation of $\boldsymbol{J}_{\mbox{\scriptsize{d}}}$. To do so, we assume that it coincides with that of the angular momentum of the neighbouring gas that ultimately feeds the BH-disc system. Although there is a large gap between the physical scale of the gas reservoir in the galactic centre and the scale of the accretion disc, it is reasonable to assume that during the feeding process, when the gas is on its way towards the disc, its angular momentum orientation is preserved.

2.1.3 Self-gravitating discs

When the mass accretion rate of the system is sufficiently large, the accretion disc becomes very massive and the effects of self-gravity are no longer negligible. This makes the disc become unstable and fragment into small gas clumps. In order to assess this more quantitatively, we use Toomre’s stability criterion which balances rotational support of the disc against its own gravity:

[TABLE]

where $Q$ is the Toomre parameter and $c_{s}$ is the speed of sound in the gas. For $Q\leq 1$ the disc is unstable, which leads to the definition of a characteristic self-gravity radius $R_{\mbox{\scriptsize{sg}}}$ at which $Q(R_{\mbox{\scriptsize{sg}}})=1$. For the $\alpha$-disc solution, the following expression is derived (Fanidakis et al., 2011):

[TABLE]

As the parameter $Q$ is a monotonically decreasing function of radius, the disc outside $R_{\mbox{\scriptsize{sg}}}$ is subject to fragmentation whereas the inner disc remains stable. Nevertheless, as only the material within the warp radius can effectively transfer its angular momentum to the BH (Volonteri et al., 2007), we include self-gravity effects only when $R_{\mbox{\scriptsize{sg}}}\leq R_{\mbox{\scriptsize{warp}}}$ (Dubois et al., 2014). This is the defining condition of what we call the self-gravity regime. As opposed to this, we define the coherent regime when self-gravity is negligible. In Figure 1 we illustrate both regimes and the evolution steps in every accretion episode.

Note that during the self-gravity regime, the innermost disc still aligns with the equatorial plane of the BH and its mass is now given by (Fanidakis et al., 2011):

[TABLE]

Using $M_{\mbox{\scriptsize{sg}}}$ and $R_{\mbox{\scriptsize{sg}}}$ instead of $M_{\mbox{\scriptsize{d}}}$ and $R_{\mbox{\scriptsize{warp}}}$ in equations (9) and (18), we can still evolve the magnitude of the BH spin.

Due to the intrinsic complexity of the self-gravity regime, the orientation of the disc angular momentum, or equivalently the material feeding the disc, has been a matter of debate. For instance, King et al. (2008) argued that low cooling times of the gas beyond the self-gravity radius provides an ideal scenario for star formation, which can occur in about a dynamical timescale. This would have disturbing effects on the fragmented disc that most likely will induce random motion of the clumps of gas. King et al. (2008) assumed that the direction of the angular momentum of the gas clumps feeding the inner disc follows an isotropic distribution over $4\pi$ steradians. This corresponds to what is dubbed chaotic accretion.

On the other hand, Volonteri et al. (2007) proposed the so-called coherent accretion, in which the angular momentum of consecutive accretion episodes is aligned with one another. In the same fashion, Dubois et al. (2014) always align, in every accretion episode, the disc angular momentum with that of the surrounding gas. We follow this approach only in the coherent regime.

Dotti et al. (2013) assumed a somewhat intermediate approach, in which they assigned two probabilities, $F$ and $1-F$, for an accretion event to have its angular momentum pointing towards the northern or the southern hemisphere as defined by the galactic reference frame, respectively. Following this line of reasoning, we assume that the angular momentum imparted by the gas accreted from large scales is not completely erased. Nevertheless, considering that self-gravity is indeed a stochastic and complex process, the disc angular momentum should not be exactly aligned with the surrounding gas either. This is why we propose a novel approach in which the degree of anisotropy of the fuelling process can be varied, allowing its effects to be studied. In order to do so, we use the von Mises distribution function, a circular analogous of the normal distribution:

[TABLE]

where $\theta$ is the polar angle, $k$ is the concentration parameter that quantifies the anisotropy degree and $I_{0}(k)$ gives the Bessel function of zeroth-order. The $z$-axis is still defined by the direction of the angular momentum of the surrounding gas.

In Figure 2, we show different parameter choices for $k$. For instance, $k=0$ corresponds to an isotropic distribution that recovers chaotic accretion. For $k=10$, the distribution is very concentrated within the cone $\theta<\pi/4$ around the direction of the surrounding gas. For $k$ sufficiently large, we recover the coherent accretion as well.

2.2 Coalescence of BH binaries

The second channel of evolution for BH spin corresponds to mergers of BHs. This is a common event in the hierarchical structure formation process, in which two merging galaxies ultimately form a deeper common potential well. After this, the two central BHs approach each through dynamical friction, form a binary, and finally coalescence. During the final inspiral phase, when the orbital angular momentum and energy are driven away by gravitational waves, we track the spin of BH remnants by using the analytical fit of Rezzolla et al. (2008), which is obtained from relativistic numerical simulations of coalescing binary Kerr BHs. The spin of the BH remnant can be expressed as:

[TABLE]

where $\boldsymbol{a}_{1}$ and $\boldsymbol{a}_{2}$ correspond to the spin parameter vectors111This is a vector pointing in the direction of the BH spin, but with a magnitude given by the dimensionless spin paramenter. of the more massive and less massive BHs, respectively; $q=M_{2}/M_{1}\leq 1$ is the mass ratio of the BH binary; $\boldsymbol{\ell}\equiv\boldsymbol{\ell}^{\prime}/(M_{1}M_{2})$, with $\boldsymbol{\ell}^{\prime}$ corresponding to the difference between the orbital angular momentum $\boldsymbol{l}$ when both BHs are widely separated from each other; and $\boldsymbol{j}_{\mbox{\scriptsize{rad}}}$ is the angular momentum carried away by gravitational waves before coalescence, i.e. $\boldsymbol{\ell}^{\prime}=\boldsymbol{l}-\boldsymbol{j}_{\mbox{\scriptsize{rad}}}$.

We use the analytical expression derived by Rezzolla et al. (2008) for the norm of $\boldsymbol{\ell}$:

[TABLE]

where $\phi_{1}$ and $\phi_{2}$ are the angles subtended by the vectors $\boldsymbol{a}_{1}$ and $\boldsymbol{a}_{2}$ with respect to $\boldsymbol{\ell}$, respectively. The adopted numerical parameters are: $s_{4}=-0.129$, $s_{5}=-0.384$, $t_{0}=-2.686$, $t_{2}=-3.454$ and $t_{3}=2.353$.

In the same fashion as Rezzolla et al. (2008), we assume the orientation of $\boldsymbol{\ell}$ to coincide with that of the total angular momentum of the BH binary before coalescence, i.e. $\boldsymbol{l}=\boldsymbol{l}_{1}+\boldsymbol{l}_{2}$. This is computed directly from the positions and velocities of the two BHs, i.e. $\boldsymbol{l}_{i}=M_{\mbox{\scriptsize{BH,i}}}(\boldsymbol{r}_{1}-\boldsymbol{r}_{\mbox{\scriptsize{CM}}})\times(\boldsymbol{v}_{i}-\boldsymbol{v}_{\mbox{\scriptsize{CM}}})$, with $i=1,2$ and CM refering to the centre of mass of the binary.

Additionally, we also compute the recoil kick velocity imparted to the remnant due to the asymmetrical emission of gravitational waves at coalescence. We use the numerical fitting formula proposed by Campanelli et al. (2007):

[TABLE]

where $\boldsymbol{\hat{e}}_{1}$, $\boldsymbol{\hat{e}}_{2}$ are orthogonal unit vectors lying on the orbital plane, and $\boldsymbol{\hat{e}}_{z}$ is perpendicular to it, i.e. in the direction of the orbital angular momentum (see Figure 3). $\perp$ and $\parallel$ refer to vector components perpendicular and parallel to $\boldsymbol{\hat{e}}_{z}$, $\Theta$ measures the angle between the vector $\boldsymbol{a}_{1}-q\boldsymbol{a}_{2}$ and $\boldsymbol{\hat{e}}_{1}$, and $\eta\equiv q/(1+q)^{2}$ is the symmetric mass ratio. The adopted numerical parameters in these equations are: $A=1.2\times 10^{4}$ km/s, $B=-0.93$, $H=7.3\times 10^{3}$ km/s, $K=6\times 10^{4}$ km/s, $\xi=1.536\ (88^{\circ})$ and $\Theta_{0}=0.184$.

3 Black Hole model

Modelling BHs in a cosmological context is a challenging task due to many fundamental uncertainties, including the nature of the non-linear coupling between galaxy formation and black hole properties; poorly understood and constrained models of mass accretion onto BHs and AGN-gas interactions; the inherent multi-scale physics of the problem, in which a comparatively small region around the BH impacts significantly largers scales in the host galaxy. These difficulties make an ab-initio treatment of small-scale BH physics within computational models of galaxy formation impossible with current numerical techniques and resources. Consequently, sub-grid approximations have to be adopted. Specifically, we use the fiducial BH model developed by Weinberger et al. (2017) and references therein. This BH model was adopted in the IllustrisTNG project as well.

3.1 Fiducial model

Our fiducial default model (based on Weinberger et al., 2017) comprises different aspects of BH physics that are crucial for modelling BH-galaxy co-evolution, namely BH growth, seeding of BHs in newly formed proro-galaxies, and feedback associated with BH accretion. In what follows, we describe the modelling of each of these aspects in turn.

3.1.1 Black hole growth

BHs gain mass and angular momentum through gas accretion from the ISM and mergers with other BHs. For gas accretion, we adopt the standard Bondi-Hoyle-Lyttleton formula (Hoyle & Lyttleton, 1939; Bondi & Hoyle, 1944; Bondi, 1952):

[TABLE]

where $\rho$ and $c_{s}$ are the density and sound speed of the gas near the BH, respectively. When the accretion rate of the BH-disc system exceeds the Eddington rate, its value is capped as follows:

[TABLE]

where $\dot{M}_{\mbox{\scriptsize{Edd}}}\equiv L_{\mbox{\scriptsize{Edd}}}/\epsilon_{r}c^{2}$ and $L_{\mbox{\scriptsize{Edd}}}$ is the Eddington luminosity as defined in the spin evolution model. For BH mergers, the masses of the two BHs in the coalescing BH binary are simply added up in the remnant, neglecting mass loss due to gravitational wave emission.

3.1.2 Black hole orbits

Due to our inability to resolve the local environment of each BH, it is not possible to accurately capture dynamical friction exerted by the background dark matter, gas and stars on the BHs. However, accounting for this process would be important for the correct computation of the BH orbits and the BH binary coalescence timescale. To circumvent this problem, the position of each BH is simply fixed to the local gravitational potential minimum of the host galaxy and its velocity is set to the mean mass-weighted velocity of the region. This has as a consequence that BH binary coalescence cannot be treated in a realistic fashion and a “numerical” coalescence must be implemented instead, i.e. two BHs merge once their relative separation is smaller than two softening lengths. Although this approach is simplistic and leads to underestimated BH coalescence timescales and fractions of free floating BHs, it is also numerically robust and prevents BHs to unrealistically migrate out of the centres of their host galaxies.

Note that we have provided an estimate of the recoil velocity in a BH binary coalescence in equation (25), however, this has no real impact on the BH population in the present set of simulations due to our special BH positioning scheme, i.e. BHs cannot recoil by construction. Therefore, we adopt the same approach as Sijacki et al. (2009), i.e. we compare the recoil velocity at coalescence time with the escape velocity of the host halo in order to assess the impact of our spin evolution model on the fraction of escaping BHs, which are the result of strong gravitational recoils. Mild gravitational recoils might not be powerful enough to expel a BH from its galaxy, but they can lead to off-centred BHs. As highlighted by Sijacki et al. (2011) and Blecha et al. (2011), BHs wandering off the galactic centre have a sizeable impact on the BH growth, increasing thus the scatter in the BH mass – host galaxy relationships. In future work we intend to include a sub-grid model of dynamical friction in order to study the impact of off-centred BHs more realistically (see also Tremmel et al., 2015).

3.1.3 Black hole seeding

Although virtually every galaxy hosts a supermassive BH, there is a significant theoretical uncertainty about the origin and early growth of these objects (e.g. Volonteri 2010). This makes it difficult to physically motivate models of BH seeding in simulations. Many recent hydrodynamical simulations of galaxy formation employ heuristic approaches, where a halo is seeded with an already massive BH upon fulfilment of a given condition, circumventing thus the phase of early growth. For instance, Bellovary et al. (2011) and Tremmel et al. (2017) seed BHs in locations where the local gas becomes dense, cold and has a low metal content. In the fiducial model of Weinberger et al. (2017) and in other previous works, a BH of mass $M_{\mbox{\scriptsize{seed}}}$ is placed at the position of the potential minimum of any halo that exceeds a total mass $M_{\mbox{\scriptsize{FOF}}}$ and does not contain a BH yet. Despite its simplicity, this method is relatively insensitive to the details of the early BH growth physics. It only assumes that this process is efficient enough to populate all halos with a low-mass supermassive black hole, whether or not this then grows further is a matter of the local conditions around it.

3.1.4 Black hole feedback

There is widespread theoretical consensus that there exist at least two types of accretion flows onto BHs (e.g. Begelman 2014, and references therein). One type operates in BHs with high accretion rates and is associated with the formation of classical accretion discs (Shakura & Sunyaev, 1973). The second type occurs only in BHs with low accretion rates and produces a more spherical and hotter accretion flow (Shapiro et al., 1976; Ichimaru, 1977). Weinberger et al. (2017) associate these types of accretion with two different BH feedback modes, namely a thermal quasar feedback mode that is activated during the high accretion state, and a kinetic feedback mode that operates in BHs with low-accreting BHs.

In the quasar mode, the energy liberated by the BH is injected thermally into the surrounding gas and is given by (Springel et al., 2005):

[TABLE]

where, as previously defined, $\dot{M}$ is the accretion rate of the BH-disc system and $\epsilon_{r}$ the radiative efficiency, which is spin-dependent in our case whereas in the fiducial model a fixed value $\epsilon_{r}=0.1-0.2$ is commonly adopted. Finally, the parameter $\epsilon_{\mbox{\scriptsize{f, high}}}$ quantifies the efficiency of the energy coupling process to the local environment.

For the kinetic mode, the liberated energy is

[TABLE]

where $\epsilon_{\mbox{\scriptsize{f, kin}}}$ is the coupling efficiency. Unlike the thermal quasar mode, the coupling process is mechanical, i.e. pure kinetic energy is injected into the gas without changing its temperature. This injection is isotropic when averaged in time, which yields a net conservation of linear momentum. Note that there is no dependence on the BH spin through $\epsilon_{r}$ as the low accretion state is thought to be radiatively inefficient. A direct spin dependence might, however, be present in reality, depending on the mechanism responsible of delivering the feedback energy in the kinetic mode. For instance, if small-scale jets produced by the Blandford-Znajek mechanism take place, the jet power has a direct dependence on the BH spin (Blandford & Znajek, 1977; Yuan & Narayan, 2014):

[TABLE]

where $B_{\mbox{\scriptsize{pol}}}$ is the poloidal magnetic field at the BH horizon. This scenario might also provide a third spin evolution channel as the rotational energy of the BH is drained by the jet. Nevertheless, in the fiducial model the kinetic feedback mode is assumed to be independent on the BH spin. In future work, we intend to return to this point and explore the consequences of including such a dependence, and also of injecting the kinetic feedback in the direction of the BH spin instead of using a random direction.

A remaining issue is the transition between the high and the low accretion modes, or equivalently, between the thermal quasar mode and the kinetic mode. Given the inability to resolve the small-scale accretion onto a BH, such transition cannot be self-consistently accounted for in cosmological simulations. Therefore, more effective approaches must be adopted instead. For instance, Sijacki et al. (2015) and other previous works have used a fixed threshold value of the fraction $\chi\equiv\dot{M}/\dot{M}_{Edd}$, i.e. when $\chi\lesssim\chi_{\mbox{\scriptsize{th}}}$ the low accretion mode is adopted, otherwise the BH is assumed to be in a high accretion state. This idea is extended in the fiducial model, and the threshold value is scaled with BH mass (Weinberger et al., 2017):

[TABLE]

with $\chi_{0}$ and $\beta$ free parameters that are calibrated to reproduce observational trends.

3.2 Self gravity in the accretion disc as a switch of BH kinetic feedback

The implementation of a kinetic BH feedback mode in the IllustrisTNG simulations has proven to be extremely important to simultaneously account for quenching of star formation and correct gas abundances in high mass galaxies (Weinberger et al., 2017, 2018). Nevertheless, apart from the notion that the kinetic BH feedback mode occurs in massive galaxies, the exact conditions under which it operates remain unknown. This uncertainty is addressed in the fiducial model by adopting an ad hoc feedback switch (equation 31) that is calibrated to reproduce observational trends. An important example for the outcome of this is the strong galaxy colour bimodality in IllustrisTNG, which is remarkably consistent with observations (Nelson et al., 2018). This represents a major improvement with respect to the first generation of Illustris simulations, in which the galaxy colour bimodality was much weaker (Nelson et al., 2015).

In spite of many encouraging successes, there are still some tensions between IllustrisTNG and observations. For example, Rodriguez-Gomez et al. (2019) report a galaxy morphology – colour relation that is inconsistent with observations, i.e. there is an excess of red discs and blue spheroids in the simulations. This discrepancy might be caused by the limited mass resolution of the simulations, which seems to be insufficient to produce starbursts (Sparre & Springel, 2016). Given that galaxy mergers trigger starbursts and drive morphological changes, and that post-merger galaxies tend to have a quiescent star forming activity, the inability to produce starbursts might naturally lead to a lacking connection between galaxy morphology and colour. Rodriguez-Gomez et al. (2019) suggest that the missing link is BH gas accretion during galaxy mergers. They argue that the limited mass resolution likely causes an underestimation of merger-induced gas inflows that feed central BHs. This would ultimately lead to an under-powered BH feedback and a potential excess of blue spheroids due to poor quenching. Nevertheless, we note that the BH mass – stellar mass relation is well reproduced in IllustrisTNG (Weinberger et al., 2018), which means that BHs are not significantly underfed. Therefore, gas accretion is not likely the source of the problem.

An alternative explanation of the tension in the morphology – colour relation could lie in the use of an inadequate BH feedback switch. Recall that the colour bimodality is well reproduced in IllustrisTNG (Nelson et al., 2018), demonstrating that the kinetic BH feedback mode is an efficient quenching mechanism. However, the inconsistent morphology means that the BH feedback switch is selecting the wrong galaxy population to quench. Thus, spheroids that should be quenched but are not, become blue, and discs that are quenched but should not, become red. The origin of this problem might be related to the ad-hoc nature of the fiducial BH feedback switch, and its lack of a physical connection to the galaxy morphology. In this case it would be very difficult, if not impossible, to calibrate equation (31) to also reproduce the morphology – colour relation in the simulations.

In order to construct a more physically motivated BH feedback switch, we start from the notion that there are two different BH accretion modes (Begelman, 2014). A high accretion mode, in which a radiatively efficient thin disc develops; and a low accretion mode, in which the accretion flow is spherical, hot, and radiatively inefficient. We note that these modes closely resemble the coherent and self-gravity regimes of the accretion disc in the BH spin evolution model (see Section 2). The coherent accretion regime is compatible with the formation of a classical Shakura-Sunyaev disc in the high accretion mode, whereas the fragmentation of the disc in the self-gravity regime provides an ideal scenario for the spherical hot accretion flow in the low accretion mode. Following Weinberger et al. (2017), we associate the quasar thermal feedback mode with the coherent regime (high accretion mode), and the kinetic feedback mode with the self-gravity regime (low accretion mode). Therefore, we can use the self-gravity condition $R_{\mathrm{sg}}\leq R_{\mathrm{warp}}$ as a BH feedback switch. A similar scenario is proposed by Garofalo et al. (2010), in which merger-induced retrograde accretion onto a BH is much more efficient in producing non-thermal jet outflows. Note that our argument is a mere working hypothesis, meant to explore the still poorly understood physical connection between BH accretion and BH feedback. Nonetheless, our approach does not require calibration, and self-consistently predicts the activation of the kinetic feedback mode based only on the physical properties of BHs, which represents a big improvement with respect to previous ad-hoc schemes.

Starting from the condition $R_{\mathrm{sg}}\leq R_{\mathrm{warp}}$, we can derive a parametric equation for the mass threshold $M_{\mathrm{bh,th}}(\dot{M},a)$ at which the transition from the coherent regime to the self-gravity regime occurs. In order to do so, we use equations (12) and (20) for the warp radius $R_{\mathrm{warp}}(M_{\mathrm{bh}},\dot{M},a)$ and the self-gravity radius $R_{\mathrm{sg}}(M_{\mathrm{bh}},\dot{M},a)$, respectively. In Figure 4, we show the BH mass – Eddington ratio diagram with two different mass thresholds for two extreme values of the spin parameter, i.e. $M_{\mathrm{bh,th}}(\dot{M},a=0.998)$ and $M_{\mathrm{bh,th}}(\dot{M},a=0.001)$. For comparison purposes, we also show the mass threshold associated to the BH feedback switch of the fiducial model (we use $\chi_{0}=0.002$ and $\beta=2.0$ as in Weinberger et al., 2017).

In this figure, we also sketch the evolutionary path of a maximum spinning BH of mass $\sim 10^{8}\ M_{\odot}$222According to the BH mass – stellar mass relation (Kormendy & Ho, 2013), a BH mass of $\sim 10^{8}\ M_{\odot}$ corresponds to a galaxy stellar mass of $\sim 10^{10}\,{\rm M}_{\odot}$, which is roughly the mass scale above which quenching becomes efficient and galaxies migrate to the red sequence. during a galaxy merger. The grey arrow represents what happens if the fiducial BH feedback switch is adopted. Likewise, the black arrow represents the BH evolution if the self-gravity switch is used instead. Considering that, during a galaxy merger, the central BH experiences a boosted growth due to gas accretion and/or coalescence with a second BH, both quantities, BH mass and Eddington ratio, will unavoidably increase. This sets a specific direction of evolution in the diagram. Adopting the self-gravity switch will clearly favour quenching of post-merger galaxies compared to the fiducial switch. This, in conjunction with merger-induced morphological changes, might help recovering the colour-morphology relation in the simulation. We also note that galaxies hosting BHs with lower spin values will be quenched at higher stellar masses.

Some high redshift quasars exhibit extremely large BH masses ($M_{\mathrm{bh}}\gtrsim 10^{9}\,{\rm M}_{\odot}$) and high accretion rates. The kinetic feedback mode should clearly not be operating under these conditions as it would expel most of gas content of the host galaxy, so the thermal quasar mode must be operating instead. The fiducial model solves this by specifying an upper limit for the Eddington rate ($\dot{M}/\dot{M}_{\mathrm{Edd}}=0.1$) for the kinetic feedback mode to operate.

For the self-gravity switch, a self-consistent solution might instead emerge out of the spin evolution model. Once the self-gravity regime is reached, the kinetic feedback mode is activated and most of the central gas is heated and/or pushed away. This results in a significant decrease of the BH accretion rate and a suppression of the gas evolution channel. After this, only BH binary coalescence can alter the spin. Nevertheless, if a gas-rich event manages to substantially increase the BH accretion rate again, the BH will quickly be spun down, thereby increasing the BH mass threshold of the self-gravity switch and activating again the thermal quasar mode. Additionally, if a spin dependent kinetic feedback based on the Blandford-Znajek mechanism were operating (equation 30), a BH with a decreasing spin would weaken the feedback energy as $L_{\mbox{\scriptsize{jet}}}\propto a^{2}$, thereby aiding further accretion. The feasibility of this scenario is yet to be tested in full; the low abundance of these objects demands simulations of large cosmological volumes at at least moderate resolutions, which is out of the scope of this work.

4 Numerical Set-up

4.1 Simulation Code

We run cosmological simulations to self-consistently track the spin evolution of supermassive black holes along with other quantities of interest such as BH mass accretion rate, BH feedback and properties of the host galaxies. We use the $N$-body magneto-hydrodynamics moving-mesh code AREPO (Springel, 2010) to solve the non-linear system of differential equations describing gas hydrodynamics and gravitational interactions of dark matter and baryons in such simulations. One of the unique features of AREPO is that the mesh continuously transforms and flows along with the fluid, unlike in standard grid-based methods. This is achieved through a dynamic unstructured mesh that is constructed as the Voronoi tessellation of a set of mesh-generating points that move with the local fluid velocity. This ensures that every cell approximately maintains a given target mass, and in this fashion, regions of high density are resolved automatically with smaller cells than low density regions, which is crucial for modelling processes such as gas accretion onto BHs residing in high density galaxy centres. Furthermore, because the mesh moves with the fluid, the advection errors inherent in Eulerian mesh codes are strongly reduced. In fact, a discretization that is manifestly Galilean invariant is achieved in this way. Also, there are no preferred coordinate directions in this scheme, avoiding potential spurious alignments of gas structures with the grid, as it can happen in Cartesian mesh codes when the resolution is low. This helps in accurately tracking the co-evolution of the spin of BHs and the angular momentum of the gas surrounding BHs in an isotropic way.

Gravitational forces in AREPO are calculated using the Tree-PM scheme (Xu, 1995), in which long-range forces are computed with a particle-mesh method and short-range forces are followed via a hierarchichal octree algorithm (Barnes & Hut, 1986). Apart from gravitational interactions and hydrodynamics, the code incorporates a complete physics model that includes: gas cooling of primordial and metal lines; a sub-resolution ISM model based on the two phase medium (Springel & Hernquist, 2003), in which stars are formed in a stochastic fashion above a density threshold of $0.13\,$cm*-3*; stellar evolution; chemical enrichment; gas recycling; kinetic stellar feedback driven by SNe and black hole physics (section 3).

4.2 Numerical implementation of spin evolution

As discussed in Section 2, the BH spin parameter $a$ can be computed in an iterative fashion through equation (8); however, an initial condition $a_{0}$ is still required. Hence, we assign an initial spin value $a_{0}$ to every BH seed. Once initialised, if a BH seed accretes gas above $1\%$ of its Eddington rate, the gas evolution channel is activated. In what follows, we discuss how to integrate BH spin evolution through gas accretion into the AREPO code.

The density and the direction of the angular momentum of the gas around BHs are needed to compute the mass accretion rate (equation 26) and the direction of the BH spin (equation 14). In order to estimate these properties, we use a spherical SPH kernel $w(r;h)$ around each BH that encloses a prescribed number of gas cells $n_{\mathrm{ngb}}$ in a radius $h$:

[TABLE]

Here, $m_{\mathrm{baryon}}$ is the target mass of a gas cell and $i$ runs over the cells within the kernel radius.

One important issue that needs to be addressed concerns the time integration of the spin evolution. On the one hand, gravitational interactions and gas hydrodynamics fix the evolution time step $\Delta t$ for every particle/gas cell of the simulation (Springel, 2010). For BHs, this time step determines how often properties such as $M_{\mathrm{bh}}$, $\dot{M}$ and the direction of the angular momentum of neighbouring gas are updated. On the other hand, the BH spin is updated after an accretion timescale $t_{\mathrm{acc}}$ has elapsed in the coherent regime, or after the disc has accreted a mass $M_{\mathrm{sg}}$ in the self-gravity regime. This implies that processes such as BH growth and BH feedback, which are governed by $\Delta t$ in the code, are not necessarily synchronised with the evolution of BH spin. This leads to four different cases:

• (i)

BH in coherent regime, $t_{\mathrm{acc}}<\Delta t$: in this case, several accretion episodes occur within the same time step. At the beginning of the time step, we store the BH properties needed for the spin model, and they are used for every accretion episode occurring within that time step. We also set a cumulative time counter that keeps track of how many accretion episodes have been completed. We stop once the counter equals $\Delta t$ and proceed to the next time step.

• (ii)

BH in coherent regime, $t_{\mathrm{acc}}>\Delta t$: in this case, several times steps elapse in the same accretion episode. We set a time counter at the beginning of the accretion episode that keeps track of how many time steps have been completed. We stop once the counter equals $t_{\mathrm{acc}}$ and proceed to the next accretion episode. The input BH properties for the spin model are averaged over the elapsed time steps.

• (iii)

BH in self-gravity regime, $M_{\mathrm{sg}}<\dot{M}\Delta t$: this case is similar to case (i), i.e. several accretion episodes occur within the same time step. Given that the BH spin must be updated once a mass $M_{\mathrm{sg}}$ is swallowed by the BH, we set, at the beginning of the time step, a cumulative mass counter that keeps track of how many times this has happened. We stop once the mass counter equals $\dot{M}\Delta t$ and proceed to the next time step. The BH properties are assumed constant for all the accretion episodes in this time step.

• (iv)

BH in self-gravity regime, $M_{\mathrm{sg}}>\dot{M}\Delta t$: this case is similar to case (ii). At the beginning of the accretion episode, we set a cumulative mass counter that keeps track of how much mass is accreted in every time step. We stop once the mass counter equals $M_{\mathrm{sg}}$ and proceed to the next accretion episode. The BH properties are averaged over the elapsed time steps.

At the beginning of every accretion episode, we compute $R_{\mathrm{warp}}$, $R_{\mathrm{sg}}$, $t_{\mathrm{acc}}$ and $M_{\mathrm{sg}}$. This allows us to assess which of the previous cases applies. At the end of an accretion episode, we update the BH spin.

The second evolution channel, BH binary coalescence, can be easily integrated into the code as BH mergers are treated in an instantaneous fashion. This makes it straightforward to apply equation (23) to compute the spin of the BH remnants.

4.3 Simulation set

We run a set of 5 cosmological simulations in a periodic box with side length of $25\,h^{-1}\ \mathrm{Mpc}$. This is a relatively small size compared to state-of-art cosmological simulations like Illustris (Vogelsberger et al., 2014), IllustrisTNG (Pillepich et al., 2018; Springel et al., 2018) or Eagle (Schaye et al., 2015); however, we still form around $1200$ BHs, with roughly $60$ of them being more massive than $10^{8}\,\mathrm{M}_{\odot}$. This makes our simulations suitable to sample the evolution of BH spin in the coherent and self-gravity regimes. For simulation 1, we do not use our BH spin evolution model as this is the control simulation used for comparison purposes. For simulations 2, 3 and 4, we use our self-consistent BH spin evolution model with different values of the parameter $k$. BH feedback and BH growth are affected by BH spin owing to the spin-dependent radiative efficiency of the accretion disc. Besides, we use the fiducial feedback switch as we want to isolate the effects of BH spin evolution. We show the results of these three simulations in Section 5. Finally, simulation 5 is run with the BH spin model as well, but additionally, we use our BH feedback switch based on the self-gravity of the accretion disc. The results of this simulation are discussed in Section 6. We summarise the model choices for our simulations in Table 1.

For all our simulations, we adopt a $\Lambda$CDM cosmology consistent with the Planck intermediate results (Planck Collaboration et al., 2016), given by $\Omega_{M}=0.3089$, $\Omega_{\Lambda}=0.6911$, $\Omega_{b}=0.0486$, $h=0.6774$ and $\sigma_{8}=0.8159$. Initial conditions are generated by using an Eisenstein & Hu (1998) power spectrum starting at redshift $z=127$. The initial conditions contain $512^{3}$ dark matter particles and the same number of gas cells. This corresponds to a mass resolution of $8.4\times 10^{6}h^{-1}\,\mathrm{M}_{\odot}$ and $1.5\times 10^{6}h^{-1}\,\mathrm{M}_{\odot}$ for dark matter particles and gas cells, respectively, which is similar to the highest resolution of Illustris TNG100 (Pillepich et al., 2018). The softening length for gas cells, dark matter and star particles is $1$ comoving kpc with a maximum value of $0.5$ proper kpc. For the black hole model, we adopt the following parameters: $M_{\mathrm{FOF}}=5\times 10^{10}h^{-1}\,\mathrm{M}_{\odot}$, $M_{\mathrm{seed}}=8\times 10^{5}h^{-1}\,\mathrm{M}_{\odot}$, $\epsilon_{\mathrm{f},\,\mathrm{high}}=0.1$, $\epsilon_{\mathrm{f},\,\mathrm{kin}}=0.2$, $\epsilon_{r}=0.2$ (only for fiducial run), $\chi_{0}=0.002$, $\beta=2.0$ and $n_{\mathrm{ngb}}=128$. For the spin evolution model, we use $\alpha=0.1$ and $a_{0}=0$.

5 Black hole spin in cosmological simulations

5.1 Individual histories

Before discussing the main demographic results of the BH population in our simulations, it is instructive to examine individual BH histories and how they relate to the spin evolution. In Figures 5 and 6, we show two representative cases of the BH evolution. The first example is a very massive BH of mass $10^{8.8}-10^{9}\,{\rm M}_{\odot}$333We report a range of BH masses because the anisotropy degree $k$ of the feeding process is varied in each simulation, which yields slightly different mass values even though we show the same black hole in each case. residing in a galaxy of stellar mass $1.1\times 10^{11}\,{\rm M}_{\odot}$. This BH is seeded at a redshift $z=8.3$ (lookback time of $13.12\ \mathrm{Gyr}$) when its host halo reaches the threshold mass $M_{\mathrm{FOF}}$. The BH survives until $z=0$, i.e. it does not merge with any more massive BH throughout its history.

The evolution of the spin parameter is shown in the top panel of Fig. 5. Once the BH is seeded, we find that gas accretion takes about $500\ \mathrm{Myr}$ to spin up the BH to the maximum spin value. Especially in the last $100\ \mathrm{Myr}$ of this period, from lookback time $12.8\ \mathrm{Gyr}$ to $12.7\ \mathrm{Gyr}$, the spin parameter increases very quickly, which is due to the accretion disc being in the coherent regime and the BH reaching the Eddington accretion rate (bottom panel). This is also accompanied by a rapid increment of the radiative efficiency from $\sim 5\%$ to $\sim 30\%$ (third panel), which represents an enhancement of a factor of six in the supplied feedback energy per unit of accreted mass.

For this BH, the coherent accretion regime extends up to $z=3.6$ (blue shade). During this time, every accretion episode is by construction aligned with the surrounding gas of the BH. Although this type of accretion is more stable than the chaotic counterpart, the orientation of the gas angular momentum can still fluctuate considerably, which is seen in the evolution of the Cartesian $z$-component of the angular momentum $j_{\mathrm{gas},z}$ (fifth panel). Owing to the high mass accretion rate of the disc in the coherent regime, the angular momentum transferred to the BH in every accretion episode is sufficient to keep the BH spin mostly oriented with the gas. This explains why the $z$-projection of the BH spin, $j_{\mathrm{bh},z}$, mimics the global behaviour of $j_{\mathrm{gas},z}$ (fourth panel), and why the angle $\theta$ is small (sixth panel). Small fluctuations in $\boldsymbol{j}_{\mathrm{gas}}$ cannot be followed up by the BH spin as they occur on timescales smaller than the accretion timescale $t_{\mathrm{acc}}$, which in turn has a lower limit imposed by the Eddington rate at a given BH mass and spin, i.e. corresponding to $\lambda=1$ in equation (13).

Between $z=4.3$ and $z=4$, the host galaxy experiences a merger with a smaller galaxy (stellar mass ratio 1:4). This event disturbs the central gas considerably and is reflected in the rapid and violent fluctuations in $j_{\mathrm{gas},z}$ (Mayer et al., 2007). Although the BH is already accreting at its Eddington rate during this time, in this particular case the BH spin is not reoriented fast enough to keep up with the gas fluctuations. This creates an interesting situation in which the BH spin is suddenly misaligned, and counter-rotating accretion takes place even in the coherent regime, i.e. for $\cos\theta<0$. The spin parameter drops to lower values as the BH is spun down, and thus the radiative efficiency decreases to $5\%$ as well. Due to the self-consistent nature of our model, the BH mass is also impacted during this time and its growth is temporarily boosted (second panel). Note that this behaviour cannot be accounted for in a static model with constant $\epsilon_{r}$. We remark that these fluctuations of the gas angular momentum are resolved in the simulation and driven by local turbulence (Dubois et al., 2014). Therefore, they are different than those produced by internal processes in the accretion disc during the self-gravity regime. Once the merger has completed, the central gas becomes stable again and the BH spin rapidly reaches the maximum value.

In general, during the coherent regime the orientation of the BH spin tends to be more stable over time than the angular momentum of the central gas, which allows to interpret the former as a time average of the latter. This might have interesting consequences for the alignment of the BH spin with the host galaxy and the cosmological environment as several processes funnel gas towards the galactic centre from galactic and extra-galactic reservoirs. We defer an analysis of this point to future work.

The self-gravity regime kicks in at $z=4.6$ once the condition $R_{\mathrm{sg}}\leq R_{\mathrm{warp}}$ is satisfied. At this point, the concentration parameter $k$ starts to considerably impact the evolution of the BH. For a fully chaotic accretion ($k=0$), the BH spin does not follow the gas angular momentum anymore, and thus its magnitude rapidly decreases as counter-rotating accretion episodes become very frequent. This reduces the radiative efficiency of the disc and enhances mass growth. For a very coherent accretion ($k=10$), the BH behaves similarly as during the coherent regime. Medium concentrations (e.g. $k=2$) yield a somewhat intermediate behaviour, i.e. the BH is fairly aligned with the gas angular momentum, but counter-rotating accretion episodes can still occur occasionally. At $z=3.2$, the kinetic feedback mode is activated and the central gas is quickly heated up and swept away, hence the gas accretion channel is suppressed and only BH binary mergers can significantly alter the BH mass and spin. We note that during the self-gravity regime, $j_{\mathrm{gas},z}$ appears to display a discontinuous behaviour. This is, however, a spurious effect due our recording of all relevant BH properties only when the BH spin is modified, which is done for computational reasons.

Finally, we show in Figure 6 a second example of a more quiescent BH of mass $10^{7.8}\,{\rm M}_{\odot}$ formed at $z=6.05$, residing in a galaxy of stellar mass $1.2\times 10^{10}\,{\rm M}_{\odot}$. The evolution of this BH is in part very similar to the previous one, e.g., during the coherent regime, the spin parameter eventually reaches the maximum allowed value while the BH spin is, with a few exceptions, mostly oriented with the gas angular momentum. Nevertheless, as this BH almost never reaches its Eddington rate, its whole evolution is considerably slower. A very interesting aspect is that a major galaxy merger at $z=0.7$ (stellar mass ratio 1:2) triggers both the self-gravity regime and the kinetic feedback mode at the same time, hinting thus a possible connection between these two processes. We explore this hypothesis more thoroughly in Section 6.

5.2 Black Hole demographics

5.2.1 Black hole spin distribution

Having shown how the spin model works in individual BHs, we proceed to study the demographics of the BH population in our simulations. We start by analysing the distribution of the absolute value of the BH spin as a function of BH mass. We plot in Figure 7 the results for simulations 2, 3 and 4, which include the spin model, but which use different values for the concentration parameter $k$.

One common result is the initial monotonic growth of the spin for BHs with masses of $10^{6}-10^{6.5}\,{\rm M}_{\odot}$, i.e. right after they are seeded. The small dispersion during this first phase indicates that all BHs follow the same path, i.e., at the beginning, BH spin evolution depends mainly on BH mass, no matter how fast the mass is accreted. This behaviour is explained by the explicit mass dependence in equation (8), and the ease with which the spin of a low-mass BH can be reoriented during an accretion episode, i.e. counter-rotating accretion is very unlikely to occur and the BH always gets spun up. For masses larger than $10^{7.5}\,{\rm M}_{\odot}$, the spin decreases to intermediate values and the distribution becomes more homogeneous. This is caused by BH binary coalescence and BHs entering the self-gravity regime, in which the frequency of counter-rotating accretion episodes is controlled by the concentration parameter. Approximately at this mass scale, the kinetic feedback mode kicks in as well, which limits gas accretion and leaves BH coalescence as the only effective spin evolution channel. This also explains why the different choices of the concentration parameter do not affect the spin distribution considerably, i.e. the gas evolution channel is briefly active during the self-gravity regime. From now on, we adopt for definiteness in all our analysis the chaotic accretion in the self-gravity regime as the default mode, i.e. we use only the simulations with $k=0$.

Although direct observational measurements of BH spins are difficult and challenging to obtain, there have been several attempts made through X-ray spectroscopy of the iron fluorescence line of the accretion disc. Reynolds (2013) made a compilation and a quality assessment of the measurements, which is necessary due to several contradictory results. We include these data in Figure 7 (orange diamonds) and compare them with our BH spin distribution. We find a very good agreement, i.e. spins tend to have a high value at intermediate BH masses, but they decrease for more massive BHs. Our results are also consistent with previous numerical studies of spin evolution in a cosmological context (e.g. Dotti et al., 2013; Dubois et al., 2014).

In order to test the robustness of our results, we have run simulations with different initial spin values, namely $a_{0}=0.5$ and $a_{0}=0.9$. In all cases, the spin distributions for masses larger than $10^{6.5}\,{\rm M}_{\odot}$ were almost identical, thereby demonstrating the insensitivity of the model on the initial conditions. This can be easily explained by the fact that BH binary mergers become common for massive BHs, and thus they rapidly erase any trace of the initial conditions and introduce a large dispersion in the spin distribution. Counter-rotating accretion also aids this as it becomes increasingly more likely to occur and is a relatively stochastic process.

5.2.2 Black hole mass – stellar mass relation

Amongst the achievements of the IllustrisTNG project and its kinetic BH feedback mode is the reproduction of the BH mass – stellar mass relation (Weinberger et al., 2017, 2018). Given that our spin model has a direct impact on BH accretion and the associated feedback through the spin-dependent radiative efficiency, we investigate how it affects the BH mass – stellar mass relation and the BH mass function. These are shown in Figure 8, top and bottom panels, respectively. In order to more easily understand the effects of spin evolution alone, we do not consider the full model with the self-gravity feedback switch in this discussion (rightmost panels), i.e. we exclude simulation 5.

Using the fiducial model as a reference, we find that the BH mass – stellar mass function in our model follows the same trend and is consistent with the observational fit of Kormendy & Ho (2013). A pronounced change in the accretion rate in the BH population at BH masses of about $10^{8}\,{\rm M}_{\odot}$ appears in both simulations as consequence of the fiducial feedback switch (equation 31). Given that the radiative efficiency has a strong non-linear dependence on the spin parameter (equation 10), periods of chaotic accretion can alter the BH mass accretion history significantly. This is evinced in the larger dispersion in the relation for our model. However, we remark that gas accretion is still self-regulated correctly by feedback for BH masses smaller than $10^{8}\,{\rm M}_{\odot}$ since the global slope of the relation is preserved. For BH masses larger than $10^{8}\,{\rm M}_{\odot}$, gas accretion is suppressed and mass growth is driven by mergers with lower mass BHs, which naturally leads to a scaling relation.

For the BH mass function, both simulations exhibit similar trends that are consistent with the observational constraints compiled by Shankar (2009) from X-ray and optical measurements. The only obvious discrepancy with observations is the overabundance of BHs of masses of $10^{7.8}-10^{8}\,{\rm M}_{\odot}$ present in both the fiducial model and our model. In order to check the validity of this finding, we recall the observational galaxy colour – stellar mass relation and the redshift independent transition from the blue to the red sequence at a stellar mass of about $10^{10}\,{\rm M}_{\odot}$ (Baldry et al., 2006; Peng et al., 2010). This stellar mass scale corresponds to a BH mass of $10^{8}\,{\rm M}_{\odot}$ according to the BH mass – stellar mass relation (Kormendy & Ho, 2013), which implies that the activation of an efficient quenching mechanism (e.g. our kinetic feedback mode) must occur precisely at the BH mass scale at which we find the BH overabundance in our simulations. We hypothesise that the origin of this feature is the suppression of gas accretion induced by kinetic feedback at this specific mass scale, which causes a “bottleneck” effect in BH growth as BH binary coalescence is abruptly left as the only active mass accretion channel.

We also compute separately the mass function of BHs that have been in the kinetic feedback mode at least once during their lifetime (red histograms) and BHs that have been always in the quasar feedback mode (blue histograms). We obtain again almost an identical result as in the fiducial model. For instance, we find that the red histograms extend over almost the entire BH mass range, which means that the kinetic feedback mode plays always a role in BH growth. However, for low mass BHs, this role is subdominant as this mode only fully kicks in for BHs more massive than $10^{8}\,{\rm M}_{\odot}$. We return to this result in Section 6 as it represents one of the more significant differences when our full model with the self-gravity feedback switch is used.

5.3 Galaxy mergers and BH binary coalescence

In this subsection, we study the effects of our spin model on the population of coalescing BH binaries. We start by analysing the evolution of different BH properties during galaxy mergers. In order to do so, we centre every BH binary history at the BH coalescence time, with negative values indicating pre-merger stages. We finally stack all individual profiles to obtain general trends. For pre-merger stages, we restrict our analysis to the most massive BH of each binary, with the additional constraint of being in the quasar feedback mode before coalescence. This guarantees that gas accretion is still active. We show the results in Figure 9.

We compute the fraction of BHs with counter-rotating accretion in the first panel. This is done by counting how many BHs experience counter-rotating accretion in every time bin. At about $20\,\mathrm{Myr}$ before coalescence, we find that the fraction monotonically increases until a maximum of $20\%$, which is reached at about $40\,\mathrm{Myr}$ after coalescence. The fraction drops to pre-merger values after $300\,\mathrm{Myr}$. Counter-rotating accretion becomes more likely to occur due to the rapid fluctuations in the central gas angular momentum induced by gas inflows and turbulence during galaxy mergers (Mayer et al., 2007). This situation was already encountered in the individual BH history analysed in Figure 5. In the second, third and fourth panels we plot the median trends for the magnitude of the spin parameter, the radiative efficiency and the angle between the gas angular momentum and the BH spin, respectively. There, we also find evidence of a rapid spin misalignment, i.e. the spin parameter drops to a median value of about $0.8$, and with it, the radiative efficiency decreases from $20\%$ to about $10\%$. In some cases, galaxy mergers are also able to trigger the self-gravity regime, which by construction leads to incoherent accretion episodes.

In the last panel, we follow the fractional change of the BH mass accretion rate with respect to the fiducal model, $\Delta\log\dot{M}_{\mathrm{bh}}\equiv\log\dot{M}_{\mathrm{bh}}-\log\dot{M}_{\mathrm{bh,fid}}$. To do so, we match every BH binary with its counterpart in the fiducial simulation; then, we follow the merger in both simulations and compare the mass accretion rate histories centred on the respective coalescence times. For pre-coalescence times, we find that BH accretion in our simulations follows the same trend as in the fiducial model as the median radiative efficiency coincides with the fixed fiducial value ($\epsilon_{r}=0.2$). For post-merger times, the drop in radiative efficiency has a twofold effect on BH accretion. On one hand, less radiated energy means more matter is available to the BH. Besides, a temporarily diminished BH feedback aids the global accretion onto the BH-disc system. On the other hand, some BHs are accreting at Eddington rates as we select only those in the quasar feedback mode before coalescence. The inverse proportionality of Eddington rate on radiative efficiency implies that mass accretion in those BHs is capped with a higher value. All this results in enhancements of up to $0.25\,\mathrm{dex}$ in BH accretion a few tens of Myr after coalescence with respect to the fiducial case.

We proceed by analysing the alignment of the two BH spins in BH binaries prior to coalescence as a function of binary mass. The results are shown in Figure 10. For almost the entire mass range, the spins in BH binaries tend to be slightly aligned as the trend goes towards positive values. The only exception appears in BHs of masses $10^{8}-10^{8.5}\,{\rm M}_{\odot}$, in which we find a more isotropic distribution. In this mass range, the self-gravity regime also kicks in, and given that we adopted the chaotic accretion scenario (i.e. $k=0$), the BH spin orientation is randomised. For less massive BH binaries, their host galaxies are mostly gas rich, and therefore, galaxy mergers are frequently wet. This implies that the two BHs in a binary shared a common gas reservoir before coalescence, and through gas accretion, both spins can be orientated in a similar direction. For more massive binaries, gas accretion is already suppressed and most galaxy mergers are dry. Spin alignment in this case is somewhat unexpected and therefore arguably more interesting. We hypothesise that galaxy alignment with the cosmic web might play a role here. For instance, Mesa et al. (2018) find that the axis connecting SDSS galaxy pairs is aligned with the host cosmic filament, showing thus that the galaxy merger process is not isotropic and depends on environment. The positive spin alignment in BH binaries might then be an imprint of a large-scale galaxy alignment.

Finally, we plot in Figure 11 the distribution of gravitational recoil velocities for remnants of coalescing BH binaries as a function of BH mass. We normalise with the local escape velocity of each host galaxy as derived from the central gravitational potential; therefore, individual systems with $v_{\mathrm{recoil}}/v_{\mathrm{esc}}>1$ are BHs that can potentially escape their host and form a population of free floating BHs (Favata et al., 2004; González et al., 2007). This is indeed the case for a few BHs with masses below $10^{8}\,{\rm M}_{\odot}$. For higher mass BHs, their host galaxies are massive enough to keep them gravitationally bound. The recoil produced by a gravitational wave burst depends on the mass asymmetry of the BH binary and the BH spins (see equation 25). The term depending on the mass asymmetry only accounts for a small fraction of the recoil velocity and is not enough to eject the BH remnant. Consequently, the configuration of the spins in a BH binary is what most likely drives the strong recoil that ultimately produces a free floating BH (Campanelli et al., 2007). Milder recoil kicks, which are more ubiquitous in almost the entire mass range, are also extremely important as they yield a population of off-centred BHs, which can impact lifetimes of AGNs, enhance central star formation, ease the formation of massive stellar cusp, and increase the intrinsic scatter of the $M_{\mbox{\scriptsize{bh}}}-\sigma_{*}$ relation (Sijacki et al., 2011; Blecha et al., 2011). Unfortunately, we cannot assess these effects self-consistently as BHs are artificially glued to the galaxy potential minima in our numerical model and thus the recoil velocity is not applied. Our result therefore merely demonstrates a capability of our model that has yet to be exploited self-consistently.

As a caveat, we also remind the reader that our treatment of BH coalescence is simplistic and does not account for processes such as dynamical friction, stellar hardening and circumbinary discs, which can exert an important influence on BH binaries. For example, Dotti et al. (2010) show that large-scale nuclear discs, which are likely to form in gas-rich mergers, can orientate both BH spins perpendicular to the binary orbital plane in a relatively short time scale ($\lesssim 10\,\mathrm{Myr}$), which ultimately results in weaker gravitational recoil kicks compared to an isotropic spin distribution. The positive alignment of BH binaries in our simulation (Figure 10) reinforces this conclusion as it shows that spins are already mildly oriented, even without including a circumnuclear disc. In the light of this discussion, our estimates of the velocity recoil of BH binary remnants in gas rich mergers are likely overestimated, and we consequently interpret them as an upper limit case. A proper treatment in which dynamical friction and other dissipative processes are included along with a self-consistent model of spin evolution is essential to properly address the effects of off-centred and free floating BHs in a cosmological context. We plan to address this in future work.

6 Black hole spin and feedback

In this section, we explore the impact of our novel BH feedback switch, which is based on the self-gravity of the accretion disc, on the galaxy population. To this end, we plot in Figure 12 the morphology – colour relation (right panels), the colour – stellar mass relation (central panels), and the morphology – size relation (left panels), for the fiducial run (top panels) and for our full model with BH spin evolution and the self-gravity feedback switch (bottom panels). We quantify galaxy morphology with the disc-to-total mass ratio (D/T), which is computed from the total mass of stellar particles within one tenth of the virial radius ($0.1\,R_{200,c}$) that have relatively high circularity ($\epsilon>0.7$) and do not belong to the spheroidal component. For galaxy colour we use the $\mathrm{g-r}$ index. Galaxy stellar mass is computed from all star particles within twice the stellar half mass radius ($2\,R_{1/2,\star}$), as in Figure 8. Finally, the stellar half mass radius is used as a proxy of galaxy size.

We recall that a strong galaxy colour bimodality consistent with observations is found in IllustrisTNG (Nelson et al., 2018). It is thus not surprising that in our fiducial run, which implements the same IllustrisTNG physics model, we also find a strong bimodality with a population transition at stellar masses of about $10^{10}\,{\rm M}_{\odot}$. Interestingly, our new feedback switch produces an even stronger colour bimodality with a population transition also at $10^{10}\,{\rm M}_{\odot}$. Specifically, our model yields a larger population of red sequence galaxies, but with a somewhat lower stellar mass, i.e. the average stellar mass of galaxies with $\mathrm{g-r}>0.6$ is $10^{10.3}\,{\rm M}_{\odot}$ and $10^{10.6}\,{\rm M}_{\odot}$ in our model and in the fiducial mode, respectively. Galaxies also exhibit smaller sizes in our model.

An analogous result appears for the BH mass – stellar mass relation and the BH mass function in our full model (see right panels of Figure 8), in which BHs more massive than $10^{8}\,{\rm M}_{\odot}$ seem to have a somewhat stagnant growth. The “bottleneck” effect observed in the fiducial model and the model with BH spin and fiducial feedback switch is also present in the full model, and is even stronger. The fact that the BH mass – stellar mass relation is still followed demonstrates that the underlying mechanism responsible for this effect is also keeping host galaxies from growing their stellar content. Looking at the Eddington ratio vs BH mass diagram (Figure 4), we can deduce that a BH of mass $10^{8}\,{\rm M}_{\odot}$ and accretion rate $\dot{M}=0.1\,\dot{M}_{\mbox{\scriptsize{Edd}}}$, which is a common scenario at high redshift ($z\sim 2$), would be promptly put into the kinetic feedback mode in our model, whereas it would remain in the quasar mode if the fiducial switch was adopted instead. This means that stellar quenching and suppression of gas accretion by AGN feedback in massive galaxies are maintained for a longer time in our model.

As described in Section 4, we adopt the same set of numerical parameters used in IllustrisTNG, which includes the coupling efficiency of the kinetic feedback mode ($\epsilon_{\mbox{\scriptsize{f, kin}}}=0.2$). Because this parameter quantifies how much energy is ultimately delivered to the gas in the form of kinetic energy, its specific value can in principle be set by tuning the simulation to reproduce specific observational trends related to BH-galaxy co-evolution. Given that we inject kinetic feedback energy to the gas at the same rate, but for more prolonged time periods than the fiducial model, we are likely overdoing AGN galaxy quenching in massive galaxies, which results in an atrophied BH and galaxy population. A potential solution to this problem is to scale down the value of $\epsilon_{\mbox{\scriptsize{f, kin}}}$ in order to deliver a similar amount of total kinetic energy as in the fiducial model, but in a more steady and prolonged fashion. A second, and perhaps more realistic alternative is to consider a spin dependent kinetic feedback (see equation 30). The feedback energy will be initially very powerful as a BH just entering the self-gravity regime should have a high spin value from the coherent regime. Nevertheless, once stochastic accretion in the self-gravity regime sets in, the BH will be rapidly spun down, and thus the feedback energy will diminish in a few hundreds of Myr. In this way, we could potentially produce a more consistent galaxy population at $z=0$, while at the same time account for high redshift galaxy quenching.

Although no parameter re-calibration was carried out in our simulations for the purpose of this initial study, our results are in remarkably good agreement with the fiducial model. An interesting discrepancy arises in the morphology – colour relation, however. As discussed in Section 3, galaxies in IllustrisTNG exhibit morphologies inconsistent with their colours, i.e. there is an excess of red discs and blue spheroids (Rodriguez-Gomez et al., 2019). This is indeed seen in our fiducial simulation, too, in which bluer galaxies ($\mathrm{g-r}<0.6$) do not exhibit a clear morphological trend, i.e. there is an almost equal fraction of spheroids and disc-like galaxies that are blue. In our model, however, a clearer trend towards disc-like morphologies appears for blue galaxies, which is more consistent with observations. The population of red galaxies is also clearly dominated by spheroids. This result is very interesting because, in spite of the somewhat stronger kinetic feedback in our simulation, our self-gravity feedback switch still manages to select the correct morphological type of galaxies to quench, even if done slightly too strongly.

With the aim of further investigating the previous result, we compute in Figure 13 the evolution of the fraction of BHs in the self-gravity regime that have masses between $10^{7.5}\,{\rm M}_{\odot}$ and $10^{8.5}\,{\rm M}_{\odot}$ and are hosted by merging galaxies. This specific mass range is where most BHs enter the self-gravity regime and, therefore, hosts galaxies which become quenched and red. We notice a significant increase in the BH fraction from $500\,\mathrm{Myr}$ before to $500\,\mathrm{Myr}$ after the merger time. This demonstrates that galaxy mergers speed up the onset of the self-gravity regime, because a boosted BH growth, either through gas accretion or BH coalescence or both, facilitates the transition from the coherent regime (see Figure 4). For our model, this has the consequence that the kinetic feedback mode is activated predominantly in merging galaxies, in which major morphological changes take place as well. This explains why we recover the galaxy morphology – colour relation in our simulations, while still reproducing the colour bimodality and the BH mass – stellar mass relation.

7 Summary and Conclusions

We have implemented a new sub-grid model for BH spin evolution in the $N$-body magneto-hydrodynamical moving-mesh code AREPO. We account for two different channels of spin evolution, namely gas accretion and BH binary coalescence. Additionally, we test a hypothesis in which the self-gravity of the BH accretion disc regulates the transition from a quasar BH feedback mode to a kinetic BH feedback mode. We have run several cosmological simulations to explore the new models, including a fiducial run with the original IllustrisTNG physics that is used as a comparison reference. Our main findings can be summarized as follows:

• •

Gas accretion in low and intermediate mass BHs ($M_{\mbox{\scriptsize{bh}}}\lesssim 10^{8}\,{\rm M}_{\odot}$) occurs in a coherent fashion as small-scale turbulence induced by a self-gravitating disc is negligible. For higher redshift ($z\gtrsim 2$), the angular momentum of the central gas reservoir around the BH can fluctuate considerably. However, low-mass BHs can be easily reoriented, and thus, although the BH spin direction can fluctuate, the net effect is an increase of the spin magnitude. An exception occurs in galaxy mergers, in which strong gas inflows can cause rapid fluctuations in the gas angular momentum that cannot be followed up by the BH, even if it is accreting already at its Eddington rate. This results in a decrease of both the spin parameter and the radiative efficiency, because counter-rotating accretion becomes more likely. A consequence of this is a boost in BH growth of up to $0.25\,\mathrm{dex}$ during the merger phase. A few hundreds of Myr after the merger, accretion becomes coherent again, and the BH is quickly spun up. In some intermediate mass BHs, a galaxy merger can trigger the self-gravity regime due to the boosted BH growth.

• •

More massive BHs ($M_{\mbox{\scriptsize{bh}}}\gtrsim 10^{8}\,{\rm M}_{\odot}$) can develop massive accretion discs in which self-gravity effects become important and counter-rotating accretion occurs frequently. This leads to a drop of the BH spin magnitude that depends on the degree of anisotropy of the incoming directions of accretion episodes that replenish the inner disc. Nevertheless, these BHs usually reside in massive, early-type galaxies in which star formation activity must be promptly quenched. Therefore, most of them are in the kinetic BH feedback mode, and thus, gas accretion is drastically suppressed. This implies that spin evolution is mainly driven by BH binary coalescence, while gas accretion in the self-gravity regime plays only a weak role in shaping the BH spin evolution. Our predicted BH spin distribution is in very good agreement with the observational data compiled by Reynolds (2013), and with previous studies of BH spin evolution in a cosmological context (e.g. Dotti et al., 2013; Dubois et al., 2014).

• •

Due to the non-linear dependence of the radiative efficiency on the BH spin, even short periods of incoherent accretion can significantly impact BH growth. This is seen in the larger scatter of the BH mass – stellar mass relation in our simulation compared to the fiducial TNG model, in which a calibrated fixed value for the radiative efficiency is adopted. The fact that our relation still follows the same trend as the fiducial model, and is consistent with the observational fit of Kormendy & Ho (2013), reveals however that self-regulated gas accretion still operates in a broadly similar fashion if BH spin evolution is accounted for.

• •

The spins of BH binaries are mildly aligned prior to coalescence, with the only exception occurring in BH binaries of masses $10^{8}\,{\rm M}_{\odot}$ to $10^{8.5}\,{\rm M}_{\odot}$, which exhibit a more isotropic distribution due to incoherent accretion induced by the disc self-gravity. For low-mass BH binaries, their coalescence occurs usually in gas-rich environments, in which a common gas reservoir can orient, through gas accretion, both spins in a similar direction. Further alignment can be reached by the action of a circumnuclear disc (Dotti et al., 2010). For more massive BHs, gas accretion is already suppressed by the strong BH feedback. Therefore, spin alignment in this case might presumably reflect a large-scale alignment of galaxies with their environment. An alignment of BH spins can weaken gravitational wave bursts from BH binary coalescence and the subsequent recoil kick of the BH remnant, and has thus potentially direct consequences for the population of off-centred and free floating BHs, as well as for the predicted gravitational wave background.

• •

We propose a novel scheme for activating the kinetic BH feedback mode that is based on the self-gravity of the accretion disc. In this scenario, massive accretion discs, in which self-gravity becomes significant, get fragmented outside their self-gravity radius $R_{\mbox{\scriptsize{sf}}}$. As the disc gets more massive, the self-gravity radius shrinks, and so does the radiatively efficient innermost disc. This leads to a natural shut down of the quasar feedback mode. On the other hand, incoherent accretion induced by disc fragmentation provides favourable conditions for a spherical hot accretion flow. Furthermore, counter-rotating accretion becomes more frequent, in which case the production of non-thermal jet outflows is highly efficient (Garofalo et al., 2010). Our conjecture is therefore to activate the kinetic feedback mode once $R_{\mbox{\scriptsize{sg}}}\leq R_{\mbox{\scriptsize{warp}}}$. We compare with the warp radius $R_{\mbox{\scriptsize{warp}}}$ as only the material inside the warped disc can effectively transfer angular momentum to the BH (Volonteri et al., 2007), i.e. self-gravity does not significantly impact mass and angular momentum accretion if $R_{\mbox{\scriptsize{sg}}}>R_{\mbox{\scriptsize{warp}}}$. For a highly spinning BH ($a=0.998$), the feedback transition occurs at a BH mass of about $10^{8}\,{\rm M}_{\odot}$, which, according to the BH mass – stellar mass relation (Kormendy & Ho, 2013), corresponds to a host galaxy of stellar mass $\sim 10^{10}\,{\rm M}_{\odot}$. This stellar mass coincides surprisingly well with the mass scale at which galaxies enter the red sequence (Peng et al., 2010), thereby making our proposed scenario a plausible explanation for the onset of star formation quenching in massive galaxies.

• •

In order to test our new feedback switch, we analyse different properties of the galaxy population and compare them with the fiducial TNG model. We find a bimodality in the galaxy colour relation that is consistent with the fiducial model, although being somewhat more pronounced. Our high mass galaxies are in general less massive and have smaller sizes, which indicates a stronger kinetic feedback in our simulations. Our proposed feedback switch activates the kinetic feedback mode in massive BHs earlier than in the fiducial model. Considering that we adopted the same set of parameters as in IllustrisTNG, including the coupling efficiency of the kinetic feedback mode $\epsilon_{\mbox{\scriptsize{f, kin}}}$, we are consequently delivering more kinetic energy to the central gas. A possible solution to this problem would be to scale down the efficiency, either by simply adopting a lower value, or by implementing a spin-dependent kinetic feedback (e.g. equation 30), in such a way that the kinetic feedback energy is injected to the gas in a more prolonged fashion.

• •

Rodriguez-Gomez et al. (2019) reported a tension between IllustrisTNG and observations with respect to the galaxy morphology – colour relation, with red discs and blue spheroids appearing to be in excess in the simulations. As expected, our fiducial run displays a similar trend, as it is based on the original IllustrisTNG physics model. When we use our new model with BH spin evolution and the kinetic feedback switch, galaxies exhibit morphologies more consistent with their colours, i.e. blue galaxies are predominantly disc types, while red galaxies are mostly spheroids. We find that galaxy mergers are efficient triggers of self-gravity in the accretion discs as a boosted BH growth facilitates the transition from the coherent accretion regime (see Figure 4). Considering that we activate the kinetic feedback mode at the onset of the self-gravity regime, and that merging galaxies undergo significant morphological changes, this explains why our approach succeeds in recovering the missing relation between morphology and colour.

Our implementation of a self-consistent BH spin evolution model into realistic hydrodynamical cosmological simulations is an important advance with respect to the first generation of BH models in which BHs were simply treated as massive sink particles, ignoring the spin. Our results demonstrate the importance of coupling the spin to the physics of BH accretion and feedback, especially during galaxy mergers. The potential role of the BH spin, through the self-gravity of the accretion disc, in shaping the galaxy colour bimodality and a correct morphology – colour relation is a further tantalising finding, suggesting that future simulation campaigns should include refined models of BH spin evolution in order to yield more realistic results.

Acknowledgements

SB acknowledges support from the International Max-Planck Research School for Astronomy and Cosmic Physics of Heidelberg (IMPRS-HD) and financial support from the Deutscher Akademischer Austauschdienst (DAAD) through the program Research Grants - Doctoral Programmes in Germany (57129429). VS acknowledges the European Research Council through ERC-StG grant EXAGAL-308037, and the SFB-881 ‘The Milky Way System’ of the German Science Foundation. Authors like to thank the Klaus Tschira Foundation. We have used NASAs ADS Bibliographic Services. Most of the computations that made possible this work were performed with Python 2.7 and their related tools and libraries, iPython (Perez & Granger, 2007), Matplotlib (Hunter, 2007), scipy and numpy (Van Der Walt et al., 2011).