# Black hole spin evolution across cosmic time from the NewHorizon simulation

Ricarda Beckmann, Yohan Dubois, Marta Volonteri, Chi An Dong-Páez, Sébastien Peirani, Joanna M Piotrowska, Garreth Martin, Katarina Kraljic, Julien Devriendt, Christophe Pichon, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Ricarda Beckmann, Yohan Dubois, Marta Volonteri, Chi An Dong-Páez, Sébastien Peirani, et al. Black hole spin evolution across cosmic time from the NewHorizon simulation. Monthly Notices of the Royal Astronomical Society, 2025, 536 (2), pp.1838-1856. 10.1093/mnras/stae2595. insu-04879158

HAL Id: insu-04879158

https://insu.hal.science/insu-04879158v1

Submitted on 10 Jan 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Black hole spin evolution across cosmic time from the NEWHORIZON simulation

R. S. Beckmann<sup>1,2★</sup> Y. Dubois<sup>3</sup> M. Volonteri<sup>3</sup> C. A. Dong-Paez<sup>3</sup> S. Peirani<sup>3,4,5</sup> J. M. Piotrowska<sup>6</sup> G. Martin<sup>7</sup> K. Kraljic<sup>8</sup> J. Devriendt<sup>9</sup> C. Pichon<sup>5</sup> and S. K. Yi<sup>10</sup>

Accepted 2024 November 8. Received 2024 November 8; in original form 2024 May 22

## ABSTRACT

Astrophysical black holes (BHs) have two fundamental properties: mass and spin. While the mass-evolution of BHs has been extensively studied, much less work has been done on predicting the distribution of BH spins. In this paper, we present the spin evolution for a sample of intermediate-mass and massive BHs from the NEWHORIZON simulation, which evolved BH spin across cosmic time in a full cosmological context through gas accretion, BH-BH mergers and BH feedback including jet spindown. As BHs grow, their spin evolution alternates between being dominated by gas accretion and BH mergers. Massive BHs are generally highly spinning. Accounting for the spin energy extracted through the Blandford-Znajek mechanism increases the scatter in BH spins, especially in the mass range $10^5-10^7$ M$_\odot$, where BHs had previously been predicted to be almost universally maximally spinning. We find no evidence for spin-down through efficient chaotic accretion. As a result of their high spin values, massive BHs have an average radiative efficiency of $<\epsilon_r^{\text{thin}}> \approx 0.19$. As BHs spend much of their time at low redshift with a radiatively inefficient thick disc, BHs in our sample remain hard to observe. Different observational methods probe different sub-populations of BHs, significantly influencing the observed distribution of spins. Generally, X-ray-based methods and higher luminosity cuts increase the average observed BH spin. When taking BH spin evolution into account, BHs inject, on average, between three times (in quasar mode) and eight times (in radio mode) as much feedback energy into their host galaxy as previously assumed.

**Key words:** methods: numerical – galaxies: active – galaxies: dwarf – galaxies: jets – quasars: supermassive black holes.

## 1 INTRODUCTION

Any astrophysical black hole (BH) can be described by two fundamental properties: its mass and its spin. Mass is the simpler of the two, accumulating through gas accretion, and mergers but never decreasing (if Hawking radiation is negligible, Hawking 1975). Spin is more complex. As well as being a vector, with both a direction and magnitude evolving, BH spin can both increase and decrease throughout a BH's lifetime, through gas accretion (Bardeen 1970), spin energy extracted by feedback (Blandford & Znajek 1977), and BH-BH mergers (Rezzolla et al. 2007; Barausse & Rezzolla 2009).

The spin of a BH is tightly linked to how efficiently it converts accreted mass to radiation:

$$L_{BH} = \epsilon_r(a) \dot{M}_{BH} c^2, \qquad (1)$$

where $L_{BH}$ is the BH luminosity, $\dot{M}_{BH}$ is the BH mass accretion rate, $c$ is the speed of light, $\epsilon_r(a)$ is the radiative efficiency of the accretion disc, and $a$ is the spin magnitude. $\epsilon_r$ depends on spin as it is determined by the structure of the accretion disc close to the BH, which, in turn, depends on the BH spin magnitude and direction relative to the accretion disc spin and magnitude. How BH spin is changed during gas accretion depends on the current spin of the BH, the angular momentum of the accreted gas and the angle between the two (Bardeen & Petterson 1975; King et al. 2005; Lodato & Pringle 2006). Accretion can drive spin to a maximum value of $a = 0.998$ (Thorne 1974), but there is some more recent theoretical work suggesting that this value could be exceeded for BHs accreting at super-Eddington rates (Sądowski et al. 2011). The spin of two BHs during the merger has important consequences for the velocity kick received by the remnant following a BH-BH merger (Lousto et al. 2010, 2012).

As the BH spin evolves continuously throughout a BH's lifetime, building a nuanced understanding of the long-term evolution of

© 2024 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1">¹Institute of Astronomy and Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="2">²Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="3">³Institut d'Astrophysique de Paris, CNRS, Sorbonne Université, UMR7095, 98bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="4">⁴Université et Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, F-06304 Nice, France</footnote>

<footnote n="5">⁵Department of Physics, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan</footnote>

<footnote n="6">⁶Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="7">⁷School of Physics and Astronomy, University of Nottingham, University Park, Nottingham NG7 2RD, UK</footnote>

<footnote n="8">⁸Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="9">⁹Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="10">¹⁰Department of Astronomy and Yonsei University Observatory, Yonsei University, Seoul 03722, Republic of Korea</footnote>

<footnote>* E-mail: ricarda.beckmann@roe.ac.uk</footnote>
BHs, and their impact on their host galaxies, requires following BH mass and BH spin evolution over cosmological time-scales. Early studies using semi-analytic models based on galaxy-merger trees or analytical arguments showed that BH spin values depend strongly on the (an)isotropy of gas accreted by the BH throughout its lifetime (Shapiro 2005; Volonteri et al. 2005; Berti & Volonteri 2008; King, Pringle & Hofmann 2008; Fanidakis et al. 2011). These early works concluded that randomly oriented gas accretion (chaotic accretion) leads to low spin values ($a = 0.1–0.3$ but with significant scatter), while coherent accretion leads to maximally spinning BHs ($a > 0.9$). Berti & Volonteri (2008) concluded that gas accretion generally dominates over the contribution of BH-BH mergers, and that mergers drive spins $> 0.7$. Volonteri et al. (2013) expanded on this work and showed that at high redshift, BHs tend to spin up rapidly due to efficient gas accretion. At low redshift, BH spins for massive BHs decrease again due to dry BH-BH mergers, while lower mass BHs tend to be spun down by chaotic accretion. Dotti et al. (2013) used a parameter study to show that low-mass and high-mass BHs respond differently to anisotropy in the accreted gas: As low-mass BHs have short re-alignment time-scales, they tend to be spun up maximally but show erratically changing BH spin direction when encountering chaotic accretion. More massive BHs do not significantly change their orientation during a given accretion event, and as a result their spin value reflects the coherence of the accreted gas: high spin for coherent accretion and low spin for chaotic accretion. Using semi-analytic models, Sesana et al. (2014) confirm that while BHs in the mass range $M_{BH} = 10^6–10^7 M_\odot$ are maximally spinning, average BH spin values decrease while scatter increases for higher BH mass. Griffin et al. (2019) confirm these trends with BH mass for both chaotic and coherent accretion models and Izquierdo-Villalba et al. (2020) focus on the connection between spin and galaxy morphology, showing that BHs in elliptical galaxies have systematically lower spin than those in pseudo or classical bulges. Maio et al. (2013) demonstrated that stellar feedback insufficiently randomizes gas flows in disc-like galaxies to spin down BHs: Even in the presence of feedback, they report BH spins of $0.6–0.9$ regardless of initial BH spin.

Pioneering work using cosmological simulations was conducted by Dubois et al. (2014b). Using cosmological zoom simulations they showed that BHs spin up over cosmic time, but can have long periods of fixed spin magnitude at $z > 3.5$ as their growth is regulated by the supernovae in their host galaxy. At lower redshift, BHs tend to be highly spinning ($a > 0.7$) and spin back up efficiently after their spin is temporarily reduced during galaxy mergers. Using a larger sample of BHs and galaxies Dubois, Volonteri & Silk (2014a) showed that generally, BH spin for BHs evolved over cosmological time-scales remains high except for the most massive BHs in the local universe ($z < 0.5$) where average spin values drop by about 10 per cent. This is because most BHs are gas-accretion dominated, and BHs quickly spin up again following reductions in BH spin due to BH-BH and galaxy mergers. This result was confirmed by Bustamante & Springel (2019), who also noted that self-consistently tracking BH spin evolution increases the scatter in the $M_{star}–M_{BH}$ relation, where $M_{star}$ is the galaxy stellar mass. Beckmann et al. (2024) showed that the distribution of BH spins depends at least in part on the evolution history of its host galaxy, with BHs in merger-free galaxies having, on average, higher spin values than those in merger-dominated galaxies.

There have also recently been several efforts to study the evolution of BH spin using sub-grid accretion disc models that continuously model the evolution of the accretion disc, rather than assuming each accretion event is independent. Fiacconi, Sijacki & Pringle (2018)

present a model for a sub-grid alpha disc evolving the BH spin, and find that for BHs with mass $M_{BH} > 10^8 M_\odot$, the angular momentum of the BH generally dominates over that of the disc, which makes spin flips more likely and accretion more episodic. For lighter BHs, the opposite is generally true and accretion is more steady. Cenci et al. (2020) used a similar model and show that the long-term evolution of the BH spin depends significantly on the radius at which the gas is assumed to circularise. Koudmani et al. (2024) present a model that seamlessly transitions from an alpha-disc at high accretion efficiencies to an ADIOS-like accretion flow at low efficiencies and find that the two-disc model significantly influences the final BH spin magnitude and direction. Finally, Bollati et al. (2024) build on the model by Cenci et al. (2020) to study the impact of spin-driven feedback on BH evolution and conclude that in most situations it is comparatively small.

A BH's spin influences its ability to grow in mass in two ways: directly, as

$$ \dot{M}_{\text{growth}} = (1 - \varepsilon_r) \dot{M}_{\text{BH}}, \qquad (2) $$

i.e. any mass not converted to radiation directly contributes into the BH mass growth rate $\dot{M}_{growth}$, and indirectly, in the sense that higher radiative efficiencies imply higher amounts of feedback energy that lead to stronger self-regulation of the BH. Zubovas & King (2019) showed that both effects combine to allow slowly spinning BHs to grow up to 20 per cent more than highly spinning BHs in the same galaxy, which could have important consequences for the assembly of the first quasars (Pacucci & Loeb 2021; Zubovas & King 2021). Super-Eddington growth could play a dual role, delivering large amounts of mass growth while spinning down the BH through the jets expected to be present for the thick accretion discs that power super-Eddington growth. However, in Massonneau et al. (2023a), we showed that such Super-Eddington growth can delay but ultimately not prevent the spinup of the BH. As a result, self-regulation remains too efficient and long-term BH growth remains stifled even for an initially non-spinning BH. The length of the delay depends strongly on the adopted jet feedback efficiency (Massonneau et al. 2023b).

Observationally, BH spin is much more difficult to measure than BH mass (see Reynolds 2021, for a recent review). In recent years, observations of massive BH spins have come predominantly from a variety of methods. The most direct probe is X-ray reflection spectroscopy (e.g. Brenneman et al. 2011; Reynolds 2013; Mallick et al. 2022; Sisk-Reynés et al. 2022), but this requires long exposure times. For larger samples of active galactic nuclei (AGN), BH spin has been estimated more coarsely from the jet luminosity (Daly 2011, 2019), the radiative efficiency (Trakhtenbrot 2014) and continuum fitting to the hard X-ray spectrum (You, Cao & Yuan 2016). While this body of work gives an exciting glimpse into the distribution of BH spins in the Universe, the overall distribution of BH spins remains poorly understood as observational sample sizes are small, current observational methods are sensitive to assumptions on the structure of the accretion disc (Riaz et al. 2020) and remain biased towards high spin values (Vasudevan et al. 2016).

In this paper, we present a study of the long-term spin evolution of intermediate-mass and massive BHs in the NEWHORIZON simulation. The paper is structured as follows: After introducing the simulation in Section 2, we present results on the spin-evolution of BHs in Section 4.1, with a specific focus on the impact of BH-BH mergers and gas accretion in Sections 4.2 and 4.3, respectively. The impact of BH spin on the luminosity function is discussed in Section 4.4, and the potential observability of BH spin in Section 5. Conclusions can be found in Section 6.
## 2 SIMULATION

NEWHORIZON is a high-resolution resimulation of an average sub-volume of the HORIZON-AGN simulation (Dubois et al. 2014a). NEWHORIZON has been presented in detail in Dubois et al. (2021).

NEWHORIZON consists of a high-resolution region with a radius of 10 comoving Mpc, with a DM mass resolution of $1.2 \times 10^6 M_\odot$ which is embedded within the 142 comoving Mpc box of HORIZON-AGN. It uses a Lambda cold dark matter cosmology consistent with WMAP-7 data (Komatsu et al. 2011) with a total matter density $\Omega_m = 0.272$, a dark energy density $\Omega_\Lambda = 0.728$, a baryon density $\Omega_b = 0.045$ and a Hubble constant of $H_0 = 70.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$. The amplitude of the matter power spectrum and power-law index of the primordial power spectrum are $\sigma_8 = 0.81$ and $n_s = 0.967$, respectively.

NEWHORIZON was performed with RAMSES (Teyssier 2002), using a second-order unsplit Godunov scheme for solving the Euler equations, and an HLLC Riemann solver with a MinMod Total Variation Diminishing scheme to reconstruct interpolated variables. Refinement follows a quasi-Lagrangian scheme where a cell is refined if its mass exceeds eight times the initial mass distribution, up to a maximum resolution of 34 comoving pc. The minimum cell size is kept approximately constant throughout by adding an extra level of refinement at expansion factor $a_{\text{exp}} = 0.1, 0.2, 0.4$, and $0.8$. We supplement this quasi-Lagrangian scheme with a super-Lagrangian refinement criterion in cells with a gas number density larger than $5 \text{ H cm}^{-3}$, which enforces refinement of cells whose size is smaller than one Jeans length.

Gas follows an ideal monoatomic equation of state with an adiabatic index of $\gamma_{\text{ad}} = 5/3$, and cooling is modelled assuming equilibrium chemistry using cooling curves from Sutherland & Dopita (1993) down to $10^4 \text{ K}$. Heating from a uniform ultraviolet (UV) background takes place after redshift $z_{\text{reion}} = 10$ following Haardt & Madau (1996). Stars are formed in cells with a gas number density above $n_0 = 10 \text{ H cm}^{-3}$, following a Schmidt relation (Kimm et al. 2017; Trebitsch et al. 2017, 2021) and using a star formation efficiency that depends on the properties of the interstellar medium (ISM) such as the Mach number and virial parametres. This leads to a stellar mass resolution of $1.3 \times 10^4 M_\odot$. Stars are assumed to have a Chabrier (Chabrier 2005) initial mass function with cutoffs at 0.1 and 150 $M_\odot$, and stellar feedback separately tracks the momentum and energy-conserving phase of the explosion following Kimm et al. (2015).

### 2.1 Black hole formation, dynamics, and mergers

BHs are formed in cells whose gas and stellar density exceeds the threshold for star formation if the cell also has a stellar velocity dispersion of more than $20 \text{ km s}^{-1}$ and is located at least 50 kpc from any existing BH. BHs form with a mass of $10^4 M_\odot$. To avoid spurious motions of BHs due to finite force resolution effects, we include an explicit drag force of the gas on to the BH, following the analytic description from Ostriker (1999). BHs are merged when their relative velocity is smaller than the escape velocity of the binary, and when they approach closer than $4\Delta x$ (~150 pc). The resulting distribution of BHs mergers in NEWHORIZON is analysed in Volonteri et al. (2020). We do not model post-merger velocity kicks here.

### 2.2 Black hole accretion model

To model the evolution of BH spin through gas accretion and feedback, we use a two-disc model, tied to the two-mode feedback model:

BHs accrete gas following un-boosted Bondi–Hoyle–Lyttleton accretion $\dot{M}_{\text{BHL}}$, based on local mass-weighted, kernel-weighted gas quantities. During each accretion event, the angular momentum of the accreted gas is measured, using the same kernel-weighting as for the accretion. Accretion is capped at the Eddington rate $\dot{M}_{\text{Edd}}$, which is computed using the spin-dependent radiative efficiency $\epsilon_r^{\text{thin}}$ assuming a thin accretion disc following Bardeen (1970) such that the accretion rate on to the BH is $M_{\text{BH}} = \min(\dot{M}_{\text{BHL}}, \dot{M}_{\text{Edd}})$. For efficiently growing BHs (Eddington ratio $f_{\text{Edd}} = \frac{M_{\text{BH}}}{M_{\text{Edd}}} > 0.01$), the effective radiative efficiency $\epsilon_r$ used to determine the distribution of accreted mass between BH mass growth and feedback (see equation 2) and the bolometric luminosity (see equation 1), is taken to be that of the thin disc: $\epsilon_r = \epsilon_r^{\text{thin}}$, following Shakura & Sunyaev (1973).

To model the impact of BHs transitioning to a thick accretion disc for low Eddington ratios, we follow Benson & Babul (2009) and attenuate $\epsilon_r^{\text{thin}}$ by a factor $f_{\text{att}} = f_{\text{Edd}}/0.01$ for an effective radiative efficiency $\epsilon_r = f_{\text{att}}\epsilon_r^{\text{thin}}$ for BHs with an Eddington ratio $f_{\text{Edd}} \le 0.01$. During each accretion event, we measure the total mass accreted on to the BH, $M_{\text{acc,BH}}$ as well as the angular momentum unit vector of the accreted gas $j_{\text{gas}}$.

### 2.3 Black hole spin evolution model

In NEWHORIZON, BH spin is modelled on-the-fly and updated according to gas accretion, BH-BH mergers and BH spin-down during feedback. All BHs are seeded with zero spin, i.e. spin parameter $a = 0$. During each accretion event, we update the spin of the BH according to the angular momentum of the accreted gas, following the model first presented in Dubois et al. (2014b) and updated in Dubois et al. (2021), which we summarize here.

To update the BH spin during each accretion event, we first assume that the angular momentum direction of unresolved accretion disc $j_D$ is the same as that of the inflowing gas measured on simulation scales, i.e. $j_D = j_{\text{gas}}$. We then measure the angle $\theta$ between $j_D$ and the current BH spin unit vector $j_{\text{BH}}$. As we cannot resolve the structure of the accretion disc and its evolution explicitly, we iteratively update the spin of the BH as follows:

First, we compute the warp radius of the disc $R_{\text{warp}}$ and the mass of the disc within the warp radius, $M_d = M_d(R_{\text{warp}})$. Then we check the radius at which the disc fragments due to self-gravity $R_{\text{sg}}$. If $R_{\text{sg}} < R_{\text{warp}}$, we limit the disc mass to $M_d = M_d(R_{\text{sg}})$. Finally, we ensure mass conservation with the accreted mass by setting $M_d = \min(M_d, M_{\text{acc,BH,remaining}})$ where $M_{\text{acc,BH,remaining}}$ is the total remaining mass to be accreted on to the BH during the current accretion event.

We then compute the disc angular momentum magnitude by integrating up to the warp radius, i.e. $J_d \sim M_d M_{\text{BH}}^{1/2} R_{\text{warp}}^{1/2}$ and the BH angular momentum magnitude $J_{\text{BH}} = a M_{\text{BH}}^{3/2} R_{\text{BH}}^{1/2}$ where $M_{\text{BH}}$ is the BH mass and $R_{\text{BH}}^{1/2}$ is the BH's Schwarzschild radius. This allows us to evaluate whether the following criteria from King et al. (2005) are fulfilled, in which case the disc is assumed to anti-align with the BH: $J_d < 2J_{\text{BH}}$ and $\cos(\theta) < -\frac{J_d}{2J_{\text{BH}}}$.

Finally, we update the BH spin direction by aligning it with the total angular momentum of BH and disc, $j_{\text{tot}} = J_{\text{BH}} + J_D j_D$. The BH spin magnitude is updated using $M_{\text{BH}}^{n+1} = M_{\text{BH}}^n + M_d$ according to one of two models: For a BH with an Eddington ratio of $f_{\text{Edd}} > 0.01$, we assume a thin disc accretion model following Bardeen (1970), in which BHs are generally spun up (for $j_D$ aligned with $j_{\text{BH}}$). For BHs with $f_{\text{Edd}} \le 0.01$, we assume a thick disc solution, with spin-up rates and $\eta = \eta_{\text{radio}} = \epsilon_{\text{MCAD}}$, where $\epsilon_{\text{MCAD}}$ takes the fourth-order polynomial form presented in Dubois et al. (2021) based on
the results of the simulated magnetically chocked accretion discs (MCAD) from McKinney, Tchekhovskoy & Blandford (2012), who predict that BHs are generally spun down due to the rotational energy extracted to power their jets.

We iterate this sub-grid model until a total mass of $M_{acc,BH}$ has been accreted on to the BH, i.e. until $M_{acc,BH,remaining} = M_{acc,BH} - \sum_n M_{d,n} = 0$ and the accretion event is completed. Each accretion event is treated as the formation of a new accretion disc i.e. $M_{acc,BH}$ and $j_D$ are independent between accretion events. This is unlike sub-grid accretion disc particle models, such as the ones used in Fiacconi et al. (2018), Cenci et al. (2020), and Koudmani et al. (2024) where the evolution of the disc is tracked continuously from one timestep to the next. All formulas used for $R_{warp}$, $R_{sg}$, $M_D$, and the change in spin magnitude can be found in Dubois et al. (2014b).

During BH-BH mergers, the spin of the remnant is calculated according to the spin of the initial BHs, and the angular momentum of the binary, following Rezzolla et al. (2007), but the orbital angular momentum of the binary is taken as random.

## 2.4 BH feedback

After each accretion event, BH feedback energy is released in the form $\dot{E}_{AGN} = \eta \dot{M}_{BH}c^2$, where $\eta$ is the efficiency with which accreted energy couples to the ISM. The feedback proceeds according to the same two-disc model also used for accretion (Section 2.2) and BH spin evolution (Section 2.3):

For a BH with an Eddington ratio of $f_{Edd} > 0.01$, we assume a thin disc accretion model and following Bardeen (1970) and release BH feedback energy in the form of an isotropic injection of thermal energy with a feedback efficiency $\eta = \eta_{quasar} = \varepsilon_f \varepsilon_r^{\text{thin}}$, with $\varepsilon_f = 0.15$ (calibrated on the BH-galaxy mass relation). For BHs with $f_{Edd} \le 0.01$, we assume a thick disc solution, with spin-up rates and $\eta = \eta_{\text{radio}} = \varepsilon_{\text{MCAD}}$, where $\varepsilon_{\text{MCAD}}$ takes the fourth-order polynomial form presented in Dubois et al. (2021) based on the results of the simulated MCAD from McKinney et al. (2012). In this mode, feedback energy is injected in the form of bipolar kinetic jets with a spin-dependent efficiency $\eta$ that follows the results of McKinney et al. (2012).

## 3 BLACK HOLE AND GALAXY CATALOGUE

The catalogue of BHs analysed here consists of all BHs in NEWHORIZON that are associated with host galaxies, which, in turn, have been identified with uncontaminated host haloes, at redshift $z = 0.25$. Uncontaminated haloes are those with high-resolution DM particles only within their virial radius.

The catalogue of DM haloes consists of all DM haloes identified by the structure-finding algorithm HOP that only contain DM particles from within the zoom region, i.e. that are at the maximum DM resolution. Galaxies are identified from the star particles of the simulation using the HOP algorithm. The galaxy catalogue consists of all main or central galaxies located within the central 0.1 virial radii of a DM halo which have a stellar mass above $10^8 M_\odot$.

To identify a galaxy's main BH, we identify the most massive BH to be contained within two effective radii of the galaxy's centre, cycling over galaxies from the most to the least massive. This BH is flagged as the galaxy's main BH and removed from the sample of unallocated BHs. We then loop over all galaxies from most to least massive for a second time, identifying all remaining unallocated BHs contained within two effective radii of the galaxy as secondary BHs. A galaxy can contain multiple BHs, but a BH can only be associated with a single galaxy. Therefore, the full sample of BHs discussed

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Distribution of BH spins versus BH mass for all $z = 0.25$, colour-coded by the fraction of that BH's mass that has been acquired via mergers, $f_{BH,merged}$. BH spin observations from Reynolds (2013), Sisk-Reynés et al. (2022), Mallick et al. (2022), and Soares & Nemmen (2020) are shown with open markers. Square markers denote BHs in quasar (thin disc) mode, while round markers denote BHs in radio (jetted, thick disc) mode in the simulations. All observations are for BHs in quasar mode.</figcaption></figure>

here contains all BHs associated with host galaxies. We exclude BHs from the sample that are contained in contaminated galaxies and haloes and 'wandering' BHs that are far from any galaxy. While this approach generally gives very satisfying matching between BHs and galaxies, it can artificially exclude BHs if their host galaxies have not been correctly identified by HOP. In the paper presented here, one of the most massive BHs, BH1049 was excluded by the automatic matching algorithm. After manually confirming that it is hosted in an uncontaminated galaxy, we have added BH1049 to our catalogue of BHs. The catalogue of BHs analysed here is the same as in Beckmann et al. (2023) except for the additional BH 1049 which was omitted in the previous paper.

## 4 RESULTS

### 4.1 The long-term spin evolution of black holes

Fig. 1 shows the distribution of absolute BH spin magnitudes $|a|$ for all BHs in NEWHORIZON at $z = 0.25$ (filled markers) in comparison to low-redshift observations (open markers).

Observationally, supermassive BH spin can be inferred in several ways. The majority of present-day spin constraints come from X-ray reflection spectroscopy – the analysis of high-energy emission from a hot X-ray corona, reprocessed in a thin accretion disc and ‘reflected’ back to an observer (e.g. García, Kallman & Mushotzky 2011; Reynolds 2013; Jiang et al. 2019; Walton et al. 2019; Mallick et al. 2022; Sisk-Reynés et al. 2022). In this approach, spin is inferred from the observed shape of the Fe Kα line at 6.4–6.9 keV, subject to relativistic broadening when emitted close to the innermost stable circular orbit (ISCO) (e.g. Fabian et al. 1989; Laor 1991; Dauser et al. 2010) and hence this method is only appropriate
for sources with prominent iron-line signatures. Another important caveat associated with reflection spin measurements is their potential vulnerability to confusion between Kα emission originating from the accretion disc at ISCO, and from high-velocity outflows launched from the disc surface (Parker et al. 2022). Although this degeneracy only affects soft-band X-ray observations and can be broken by including a constraint on the Compton hump feature peaking at ~30 keV (first demonstrated by Risaliti et al. 2013), in practice such measurements require long coordinated exposures on soft- and hard-X-ray observatories. For this reason, only ~50 per cent of supermassive BH spin measurements to date were obtrained with broad-band X-ray coverage (see Piotrowska et al. 2024, for further discussion).

An alternative, albeit less common, approach involves continuum fitting of thermal emission from the AGN disc and relies on the connection between BH spin and expected disc temperature (e.g. Done et al. 2013; Capellupo, Wafflard-Fernandez & Haggard 2017). Since this method does not rely on strong iron-line emission, it can, in principle, be applied to a wider range of sources. However, its requirement of generous multiwavelength coverage and, as of yet inaccessible, extreme UV measurement render this approach challenging in practice (Capellupo et al. 2017). Although different in spirit and precision, both of these methods can only be applied to bright, highly accreting AGN, as they both rely on the assumption of geometrically thin/optically thick accretion disc models of Shakura & Sunyaev (1973) and Novikov & Thorne (1973) to establish connection between BH spin and accretion disc emission down to the ISCO.

Another approach is to estimate the BH spin from the jet power while making an assumption about the BHs radiative efficiency and Eddington ratio, either directly from observed radio jets or by looking at jetted quasars (Soares & Nemmen 2020). As can be seen by the shape of the markers in Fig. 1, BHs in NEWHORIZON are predominantly in radio mode (i.e. have $f_{Edd} < 0.01$) at $z = 0.25$, with only four BHs in quasar mode. However, the variability of $f_{Edd}$ is high, as discussed in Section 4.3. We explore the impact of this variability, as well as the impact of luminosity on the observable sample of BH spins, further in Section 5.

When looking at the distribution of NEWHORIZON BH spins, a clear trend emerges with both the current mass of the BH, $M_{BH}$ and with the fraction of this mass gained through BH-BH mergers, $f_{BH,merged}$. We remind the reader that all BHs in NEWHORIZON are seeded at $M_{BH,0} = 10^4 M_{\odot}$ with $|a| = 0$. In Fig. 1, we identify three dominant regimes for low-redshift BH spins:

(i) A gas accretion-driven spin-up for BHs close to their seed mass: after being seeded, BHs are spun up to a maximum spin ($a = 0.988$) while doubling to tripling their mass to $2 \times M_{BH,0}$ (here $2 \times 10^4 M_{\odot}$) through gas accretion. This suggests highly coherent accretion during this early spin-up phase as an initially non-spinning BH can be maximally spun up by coherently accreting approximately 1.5 times its original mass (Bardeen 1970). As a consequence, the long timescales modelled here mean that our results depend little on the choice of initial spin value for the BHs, unlike shorter studies such as, for example, Cenci et al. (2020). Some BHs continue to grow purely through gas accretion from this point onward (i.e. the fraction of their mass that they gain through mergers, $f_{BH,merged}$ remains zero) until they roughly double their mass again to $\sim 4 \times M_{BH,0}$.

(ii) A merger-induced scattering of spins for low-mass BHs: in the mass range $2-5 \times M_{BH,0}$ (here $2 \times 10^4 - 5 \times 10^4 M_{\odot}$), the distribution in BH spin becomes much broader due to BH-BH mergers: BHs that have undergone a merger (i.e. that have $f_{BH,merged} > 0$)

have $0.4 < |a| < 0.988$. Some BHs in this mass range will have acquired all their spin through the BH-BH merger (due to the orbital angular momentum of the merger, the merger-remnant BH of a merger between two non-spinning BHs will have $a = 0.69$, Berti & Volonteri 2008), while others will have at least partially spun up through gas accretion and then had their spin reduced by the merger: highly or maximally spinning BHs statistically decrease their spin during mergers due to the misalignment between the spin of the two merging BHs (Berti & Volonteri 2008). Peirani et al. (2024) showed for the same sample of BHs analysed here that during this phase gas accretion is inefficient as the repeated mergers of the host galaxy lead to rapid realignment of the stellar and gas angular momentum with respect to that of the BH.

(iii) An accretion-driven evolution for massive BHs: As BHs grow beyond $10 \times M_{BH,0}$ (here $10^5 M_{\odot}$), accretion again dominates the BH mass budget ($f_{BH,merged}$ decreases). BH-BH mergers continue to take place during this late evolution phase (see Fig. 2 below) but the mass ratio of mergers decreases as the likelihood of encountering an equally massive BHs becomes smaller. As a result, BHs assemble their mass predominantly through gas accretion. Early on, the scatter remains high, while for the most massive BHs in NEWHORIZON, spin magnitudes again increase towards maximally spinning.

While our sample of massive BHs is too small to draw firm conclusions on the spin distribution of massive BHs in the present-day Universe, those that lie within the same BH mass range as the BHs observed in the local Universe also show absolute spin magnitudes, $|a|$, in good agreement with the observations as can be seen in Fig. 1.

Fig. 2 shows how the spin of the five most massive BHs in NEWHORIZON evolves over time. The same three trends as discussed in Fig. 1 qualitatively emerge from these timeseries: First, there is an early spin-up of newly formed BHs to intermediate spin magnitudes ($0.3 < |a| < 0.9$). For the next few Gyr of BH evolution, BH-BH mergers can produce large jumps in $|a|$, often but not consistently reducing it. At late times ($z < 2$), most mergers have a modest impact on $|a|$ and gradual change through gas accretion dominates the evolution of the BH spin magnitude. This is discussed further in Section 4.2.

The time evolution of these five BHs is explored further in Figs 3–7, which show (from the top to bottom panels): BH mass $M_{BH}$, absolute and intrinsic spin magnitude $a$, the angle between accreted gas and BH spin $\theta_{spin,gas}$, the Eddington ratio $f_{Edd}$, the spin-dependent radiative efficiency $\epsilon_r$ and the feedback coupling efficiency $\eta$ for the five most massive BHs in NEWHORIZON at $z = 0.25$. The angular momentum of the accreted gas is measured as described in Section 2.2. A few quantitative trends emerge: BHs tend to have a chaotic, inefficient accretion phase early on in their evolution, which lasts anywhere from 2 to 6 Gyr. After this, accretion becomes more coherent and efficient, growing BHs in mass and spinning them up. As discussed in more detail in Section 4.3, this time evolution is at least partially driven by the mass evolution of the BH’s host galaxy.

At late times, two separate trends can be seen, linked to the accretion states of BHs: BHs predominantly in ‘quasar’ mode (BH 166 in Fig. 3, 455 in Fig. 4), which have an average $f_{Edd} = 0.03-0.6$ over the last 3 Gyr of simulation are spun up by gas accretion, even after temporary spin-down following mergers (see Fig. 2). BHs growing less efficiently (here BH 936 in Fig. 6 and 549 in Fig. 7, which have an average $f_{Edd} \sim 0.004$ over the last 3 Gyr. We note that BH 936 spins up again at late times.) are spun-down on Gyr timescales as spin energy is extracted to power the jets. BH 1049 shows an interesting mix of both behaviours, with a long spin-down phase followed by a rapid spin-up phase, followed by more spin-down. We
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Timeseries for the absolute spin magnitude $|a|$ for the five most-massive BHs in NEWHORIZON (equivalent to a mass cut of $M_{BH} > 2 \times 10^5 M_{\odot}$). Circular markers denote the timing of major BH–BH mergers (mass ratio > 1 : 4), while crosses denote minor mergers.</figcaption></figure>

discuss the impact of accretion mode on spin evolution further in Section 4.3.

Overall, our results agree with previous simulation studies on the long-term spin evolution of massive BHs using cosmological simulations by Dubois et al. (2014a) and Bustamante & Springel (2019) that showed that the spin of massive BHs is generally high and can remain roughly constant on Gyr time-scales. Unlike these works, we find examples of non-maximally spinning BHs in the mass-range $10^6-10^8 M_{\odot}$, which are spun down by coherent but inefficient accretion, in which BH spin energy is extracted to drive BH jets. This effect was not modelled in the previous two works. Unfortunately, our sample is too small to draw firm conclusions on how this affects the total spin distribution, but the two examples reported here provide intriguing hints, together with recent observations by Mallick et al. (2022), that BH spin-down by jets could play an important role in broadening the distribution of BH spins in the mass range $10^6-10^8 M_{\odot}$. Sala et al. (2024) reported that BH maximally spin up as they double their seed mass, but then report an intermediate period of high spin in the mass range $10^6-10^7 M_{\odot}$ before a mix of BH mergers and spin-down due to chaotic accretion increases scatter in spin at late times for BHs in the mass range $10^7-10^{10} M_{\odot}$ even without jet spindown, as also suggested by Fiacconi et al. (2018).

## 4.2 The impact of BH mergers on BH spin

The decreasing importance of mergers with BH mass in the mass range explored here can be seen more quantitatively in Fig. 8, which explicitly shows the distribution of the changes in spin magnitude $\Delta a$ as a function of the mass of the primary BH for all BH–BH mergers in NEWHORIZON. This plot shows a sample of stacked BH mergers that occur at different redshifts throughout the simulation. The same three phases as in Fig. 1 can be broadly identified:

(i) First BH mergers, during which the mass of the primary, $M_{BH,primary} < 2 \times M_{BH,0}$ (blue background). Such BHs are still in the process of spinning up through early gas accretion. BHs very close to the seed mass have very low $|a|$ pre-merger (see Fig. 1) and are therefore statistically spun up by mergers ($\Delta a > 0$). The opposite is true for BHs that have grown through gas accretion up to twice their original mass: they tend to be highly spinning pre-merger (see Fig. 1) and are therefore likely to be spun down by the first merger ($\Delta a < 0$). However, mergers of non-grown BHs are more common, so the average merger during this early phase will spin BHs up (see right-hand panel of Fig. 8). By definition, BH mergers for BHs close to the seed mass tend to have mass ratios close to unity.

(ii) An intermediate phase, where the primary is in the mass range $2-10 \times M_{BH,0}$ (green background). These are BHs that have already undergone either significant gas accretion or previous BH–BH mergers and will therefore have a range of spin values and orientations pre-merger. Spin changes during this period can be both positive and negative, but, on average, reduce the BH spin (see right-hand panel of Fig. 8). Mergers during this phase tend to remain major mergers (mass ratio $q > 1 : 4$) but decrease with increasing primary mass.

(iii) The long-term evolution of massive BHs (red background): Mergers continue to take place and continue to, on average, decrease the spin magnitude, but the average spin change per merger becomes significantly smaller as the average mass ratio $q$ between merging BHs decreases. The two major mergers that occur for massive BHs in NEWHORIZON also show that larger spin-changes are possible for massive BHs as long as the merger ratio is sufficiently high.

We note that having such clean mass transitions between the regimes is partially due to the uniform seed mass of $M_{BH,0} = 10^4 M_{\odot}$.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Timeseries for BH 166 of (from the top to bottom panels) BH mass $M_{BH}$, absolute and intrinsic BH spin magnitude $|a|$ and $a$, angle between accreted gas and BH spin $\theta_{spin.gas}$, Eddington ratio $f_{Edd}$, spin-dependent radiative efficiency $\epsilon_r$ and BH feedback coupling efficiency $\eta$. Dotted vertical lines mark the time of BH–BH mergers with a minimum mass ratio of $q = 0.1$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Timeseries for BH 455. See Fig. 3 for details.</figcaption></figure>

and initial spin ($a = 0$) for all BHs in NEWHORIZON. One consequence of a uniform initial mass function of BH seed masses is that all early mergers are by necessity major mergers. For a more diverse initial mass function of $M_{BH,0}$, we would expect early mergers to

continue to play an important role in increasing the scatter in BH spins for BHs of $10^5 M_{\odot}$, but the mass transitions would become less clear as early mergers would encompass a range of minor and major mergers with different mass ratios.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Timeseries for BH 1049. See Fig. 3 for details.</figcaption></figure>

**4.3 The impact of gas accretion on black hole spin**

Understanding how gas accretion influences the long-term spin evolution of BHs is more complex because the efficiency with which

angular momentum is transferred from the gas to the BH (i.e. with which the BH is spun up or down by gas accretion) depends on (1) the existing spin of the BH (2) the long-term alignment between the accreted gas angular momentum vector and the spin vector of the BH,
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Timeseries for BH 936. See Fig. 3 for details.</figcaption></figure>

here called $\theta_{gas,BH}$, and (3) the Eddington ratio $f_{Edd}$ as it determines the BH feedback mode.

BHs with $f_{Edd} > 0.01$ are assumed to be in ‘quasar’ mode, during which gas accretion will always spin BHs up if the angular

momentum of the accreted gas is aligned with the BH spin, with the spin-up rate depending on the current spin of the BH. BHs with $f_{Edd} < 0.01$ are assumed to be in ‘jet’ mode where the combined action of angular momentum transferred to the BH and BH spin
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Timeseries for BH 549. See Fig. 3 for details.</figcaption></figure>

energy extracted to power the jet means that $|a| \to 0$ over time (see section 2.5.2 in Dubois et al. 2021, for details).

Here, we assume that the angular momentum direction of the accretion disc during any accretion event is equal to the accreted

gas' angular momentum direction. Long periods of coherent gas flows on galactic scales allow BH spin vectors to realign themselves with the angular momentum vector of the inflowing gas, which then creates ideal conditions for efficiently spinning BHs up. By
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Left panel: distribution of BH spin changes $\Delta a$ during BH–BH mergers versus the primary BH mass at the time of merger, $M_{BH,primary}$ (left panel). Right panel: Gaussian fit to the distribution for the total sample (black, solid) and the three subsamples (coloured, dashed and dotted).</figcaption></figure>

contrast, if the angular momentum of the inflowing gas frequently changes, due to chaotic clumpy accretion on to the BH, we expect the spin magnitude $a$ to decrease if accretion is sufficiently efficient as retrograde accretion carries more specific angular momentum than pro-grade accretion.

The interplay between $\theta_{gas,BH}$, $f_{Edd}$, and $a$ is complex, as can be seen by the evolutionary histories of the most massive BHs shown in Section 4.1. Qualitatively, all BHs undergo an early period where $f_{Edd}$ is generally low and $\theta_{gas,BH}$ varies rapidly. Early in their evolution, when galaxies are regulated by stellar feedback (Dubois et al. 2015; Bower et al. 2016; Anglés-Alcázar et al. 2017; Habouzit, Volonteri & Dubois 2017; Trebitsch et al. 2017), they have highly disturbed morphologies and lack coherent gas inflows into the centre. As a result, gas falling on to the BHs has rapidly changing angular momentum, and $a$ oscillates between being positive and negative, as the sign of $a$ is measured with respect to the angular momentum of gas that is currently being accreted. This stochasticity could partially be due to the resolution of the simulation. It has been shown that resolving the gas flows around BHs in more detail can turn accretion from chaotic to episodic as a nuclear disc around the BH becomes resolved (Beckmann, Devriendt & Slyz 2019; Hopkins et al. 2024), or make accretion more stochastic as the inner disc fragments (Levine, Gnedin & Hamilton 2010) depending on the conditions of the disc.

This phase can last anywhere from 1 to 3 Gyr, depending on the BH, and is also a time of frequent BH–BH mergers with sufficiently high mass ratios to significantly boost the BH mass and influence the absolute spin magnitude $|a|$. For most BHs, eventually $\theta_{gas,BH}$ settles down and $f_{Edd}$ increases as the BH enters an accretion-driven regime where it grows steadily in mass. Then follows a period of coherent accretion on to the BH where BH spin is well aligned with that of the gas in the centre of the galaxy, which we assume flows coherently all the way to the accretion disc. Mergers continue to occur but as merger ratios drop, the impact on the BH mass and spin evolution diminishes. We note that $\theta_{spin,gas}$ only measures the relative angle between the BH spin and the accreted gas angular momentum. It does not quantify the direction of either gas or BH with respect to the host galaxy or large-scale cosmological environment. $\theta_{spin,gas} \sim 0$ means that the angular momentum of accreted gas is changing slowly

enough for the BH to continue aligning with it, not that the BH spin vector (and accreted gas angular momentum vector) is fixed in 3D space during this period of time.

In Fig. 9, we explore the impact of gas accretion on BH spin evolution more quantitatively. In this plot, we compute the average absolute spin magnitude $|a|$ and the average angle between gas angular momentum and BH spin vector every 0.5 Gyr over intervals of 100 Myr for the five most massive BHs in NEWHORIZON. Markers show the average value of $|a|$ and $\theta_{gas,BH}$ over each interval of 100 Myr, while grey error bars denote the variance in $\theta_{gas,BH}$ over the same time interval. Colour-coding indicates the average $f_{Edd}$ for each interval, and marker size denotes cosmic time, with smaller markers presenting earlier sampling times. As can be seen in Fig. 9, early on all BHs shown here undergo chaotic accretion (error bars in $\theta_{gas,BH}$ are large, which shows a lack of coherence in the accreted gas). During this time, BH accretion is very inefficient ($f_{Edd}$ is low). Peirani et al. (2024) showed that this chaotic phase is driven by a rapid reorientation of the stellar angular momentum, rather than by a reorientation of the BH spin. Eventually either the BH undergoes a period of efficient, more coherent accretion (small, reddish markers) or the BH is spun up by an early merger (see Section 4.2).

At later times, as the host galaxy grows in mass and settles, the coherence of the accreted gas settles down (i.e. error bars in $\theta_{gas,BH}$ become small) and all BHs studied here undergo a period of efficient accretion. Peirani et al. (2024) showed that this late transition occurs around a host galaxy mass of $5 \times 10^9 - 5 \times 10^{10} M_{\odot}$.

During this time, gas accretion keeps BH spins high, as the high $f_{Edd}$ put BHs in quasar feedback mode, allowing them to be spun even when BH mergers briefly decrease their spin magnitudes. Examples of such BHs are BH 166 and BH 455, which both show persistent high spin magnitudes at late times (see also Fig. 2). For other BHs, $f_{Edd}$ eventually drops, pushing BHs into jet feedback mode. As a result, such BHs are slowly spun down over Gyr time-scales, as can be seen for BHs 936 and 549 in Fig. 2). Interestingly, BH 1049 undergoes a period of persistent anti-aligned accretion, also evident in its timeseries in Fig. 5 around $z = 0.5$ that causes slow spindown.

We investigate the impact of (in)coherence further in Fig. 10 which is based on the same data points as Fig. 9. We now show the variance
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Distribution of the absolute spin magnitude $|a|$ versus the angle between accreted gas angular momentum and BH spin vector, $\theta_{gas,BH}$, for the five most massive BHs in NEWHORIZON (the same ones as shown in Fig. 2). BHs are sampled every 0.5 Gyr, with markers denoting the average $|a|$ and $\theta_{gas,BH}$ over $\Delta t = 100$ Myr. Each marker is colour-coded by the average $f_{edd}$ over the same period. Grey error bars show the variance of $\theta_{gas,BH}$ over each 100 Myr time interval. Small error bars indicate a period of high coherence in the accreted gas, while large error bars are a sign of chaotic accretion. Marker size encodes cosmic time, with smaller markers presenting earlier points in a BHs evolution history.</figcaption></figure>

$\Delta\theta_{BH,gas}$ on the x-axis (shown as the length of the error bars in Fig. 9) versus the Eddington ratio $f_{edd}$. As can be seen in the top panel of Fig. 10, many BHs undergoing highly incoherent accretion (high $\Delta\theta_{BH,gas}$) show little spin change $\Delta a$ because they are growing inefficiently (low $f_{edd}$). Contrary to expectation, we find that BHs that are undergoing efficient incoherent accretion (high $f_{edd}$) are actually spun up ($\Delta a > 0$). We quantify the anisotropy in the accreted gas using the anisotropy parameter $F$, following the analysis in Dotti et al. (2013). $F$ measures the fraction of accretion events for which the dot product between BH spin and gas angular momentum is positive. Truly random accretion has $F = 0.5$ while ($F < 0.5$) $F > 0.5$ suggests, on average, (anti-)aligned accretion. We find that for our BHs, even during their highly incoherent accretion phases, $F$ remains sufficiently high to cause long-term spinup. The few BHs that show high $f_{edd}$ and low $\Delta a$ are almost maximally spinning already (black outline), so cannot be spun up further. Note the anti-aligned period of accretion of BH 1049 again which leads to a net spin-down during this time.

Overall, we find little evidence for an early suppression of spin-
up due to chaotic accretion in the early Universe or the late spin-
down of massive BHs due to chaotic accretion (Volonteri et al.
2013; Bustamante & Springel 2019; Griffin et al. 2019). This is
unlike previous results reported by Izquierdo-Villalba et al. (2020),
Bustamante & Springel (2019), and Sesana et al. (2014), who
stress the importance of spin-down through chaotic accretion to
produce the diverse distribution of BH spins observed in the local
Universe. Instead, in NEWHORIZON a qualitatively similar spin-down
is produced by extracting BH spin to drive jets for some massive
BHs, an effect that was not included in the works cited above.

This could mean that previous works overestimated the impact of
chaotic accretion to compensate for the lack of jet-induced spin-
down. Alternatively, it might show that our assumption of the gas
at BH accretion scales inheriting the angular momentum from gas
at resolution scales makes efficient accretion artificially coherent,
reducing spin-down via chaotic accretion. A combination of both
effects is also possible.

In conclusion, through a combination of efficient thin disc accre-
tion, non-negligible periods of anti-aligned accretion and jet-induced
spin-down, NEWHORIZON predicts a diverse population of BHs at
low redshift with the potential for a large dispersion in spin values.
Unfortunately, our sample is too small, and does not reach very high
BH masses, to robustly quantify the distribution of massive BH spins
in the local Universe or study the impact of environment driving this
diversity.

We find little evidence for spin-down through chaotic accretion
even at early times, as the chaotic accretion episodes we observe
early-on in the evolution of most of our massive BHs is either too
inefficient to significantly influence BH spin or leads to a new spin-up
due to persistent coherence during such episodes.

4.4 The impact of black hole spin on feedback efficiencies

One important consequence of following BH spin is that it changes
both the radiative efficiency $\epsilon_r$, which controls mass growth and
luminosity of the BH, and the feedback efficiency $\eta$, which controls
how much of the feedback energy of the BH couples to the ISM
close to the BH (see Section 2.4 for details) Both potentially have
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Distribution of the variance of $\theta_{\text{gas,BH}}$ versus the average Eddington ratio $f_{\text{edd}}$ for the five most massive BHs in NEWHORIZON sampled every 0.5 Gyr over 100 Myr intervals, as in Fig. 9. Black outlines show BHs with average spin magnitude $|a| > 0.95$ i.e. almost maximally spinning. Data points are colour-coded by the total spin change due to accretion $\Delta a_{\text{acc}}$ during each interval (top panel) and by the anisotropy parameter $F$ (bottom panel, see Dotti et al. 2013). The impact of mergers has been removed from $\Delta a_{\text{acc}}$.</figcaption></figure>

important consequences for BH and galaxy coevolution, as well as their observability.

In the absence of spin information, most existing galaxy evolution simulations use a fixed efficiency of $\epsilon_r = 0.1$ as radiative efficiency, often attenuated by a factor calibrated on low-redshift observation for the feedback efficiency. In our model, the radiative efficiency $\epsilon_r$, and feedback efficiency $\eta$ are both based on analytic models that take spin evolution into account (except for a factor of 0.15 in the quasar feedback efficiency used for calibration, see Section 2.4 for details).

In Fig. 11, we show the spin-based distribution of $\epsilon_r^{\text{thin}}$ and $\epsilon_r$ (where $\epsilon_r = f_{\text{att}}\epsilon_r^{\text{thin}}$ includes the attenuation for BHs in radio mode, see Section 2.4) for all BHs at four different redshifts (left panel), and the feedback efficiency $\eta$ (right panel). Both are highly variable (see the bottom two panels of Figs 3–7 for examples). Fig. 11 shows that BHs close to the seed mass have a wide range of $\epsilon_r^{\text{thin}}$. The bimodal nature of the distribution is driven by BHs that are aligned ($a > 0$, upper branch) and anti-aligned ($a < 0$, lower branch) with the angular momentum of gas that they are accreting. While some BHs in the mass range $10^5-10^6 M_{\odot}$ also have thin discs radiative efficiency $\epsilon_r^{\text{thin}}$ below the classic value (horizontal line), most massive

BHs have $\epsilon_r^{\text{thin}} > 0.1$. The average value for all BHs with masses $M_{\text{BH}} > 10^5 M_{\odot}$ is $< \epsilon_r^{\text{thin}} > = 0.19$. This means that massive BHs in quasar mode in NEWHORIZON will shine two to three times as brightly as at a given accretion rate as previously assumed, which could also have important consequences when using the luminosity to derive the spin value of BHs (see e.g. Soares & Nemmen 2020, who assume an efficiency of 0.15). However, attenuating the radiative efficiency by a factor of $f_{\text{att}}$ for BHs in radio mode, to account for the fact that such BHs have radiatively inefficient accretion discs, has a significantly bigger impact on the average radiative efficiency than using a spin-driven $\epsilon_r^{\text{thin}}$, as can be seen by the fact that most AGN in jet mode, even those with high masses, end up with an effective radiative efficiency $\epsilon_r < 10^{-3}$ (shown as grey upper limits). As a result, the average $\epsilon_r$ for all AGN above $10^5 M_{\odot}$ is only $< \epsilon_r > = 0.1$, but the distribution is very double peaked between those in quasar mode (which typically have $\epsilon_r > 0.25$) and those in radio mode (which typically have $\epsilon_r < 0.1$). From this, we conclude that the high spin values of massive BHs will make those with thin accretion discs easier to observe but delay their mass growth. As most BH accumulate most of their accreted mass in quasar mode, this could be the reason why central BHs in massive galaxies in NEWHORIZON tend to lie towards the lower end of the observed distribution of BH masses (see Beckmann et al. 2023). For BHs with thick discs, observing their disc emission remains challenging, even with high spins, but mass growth is more efficient for such BH as almost all the mass flowing on to the accretion disc from large scales can be accreted due to their low $\epsilon_r$.

Following the spin evolution of the BH also allows us to more accurately estimate the amount of feedback energy from the BH that couples to the ISM. This quantity is encoded in the feedback efficiency $\eta$, shown in the right-hand panel of Fig. 11 as a function of BH mass. Here, the impact of spin on the distribution is even more striking. When taking the spin into account, BHs in the mass range $> 10^5 M_{\odot}$ now have an average $< \eta > = 0.04$, roughly three times higher than the fiducial value, and $< \eta > = 0.8$ in jet mode which is eight times higher than the previous assumptions based on a fixed spin and hence fixed radiative efficiency of $\epsilon_r^{\text{thin}} = 0.1$. So tracking spin evolution means that all galaxies hosting massive BHs experience significantly more feedback energy for a fixed accretion rate than previously assumed.

The overall expected impact on galaxy evolution from an increased feedback efficiency is difficult to estimate. First, the effect is not uniform across a galaxy’s evolution: only as BHs spin up over time does their feedback efficiency increase. Secondly, the total amount of feedback energy experienced by a galaxy remains determined predominantly by variations in the accretion rate (see Figs 3–7), rather than the feedback efficiency. As BH feedback self-regulates, the overall amount of energy experienced by a given galaxy might not be so different for a fixed-spin and a spin-driven model.

In Fig. 12, we test whether using a spin-based $\epsilon_r$ has a significant impact on the bolometric AGN luminosity function by comparing the spin-based luminosity function (solid line) to the one that assumes $\epsilon_r = 0.1$ for all BHs, as in previous comparable work. As can be seen by comparing the two lines, the difference in bolometric luminosity is negligible, both at high ($z = 2$) and low ($z = 0.25$) redshift. The fact that NEWHORIZON generally underpredicts the luminosity function at low redshift is explored further in Beckmann et al. (2023). Two further lines in Fig. 12 show that even under extreme assumptions (all BHs are maximally radiatively efficient, dash-dotted line and all BHs are minimally radiatively efficient, dotted line, assuming their current accretion disc structure) the impact of a spin-dependent $\epsilon_r$ on the luminosity function remains small (within a factor of 3) for the sample of BHs in NEWHORIZON.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Distribution of the spin-dependent radiative efficiency of a thin disc $\epsilon_r^{\text{thin}}$ and the effective radiative efficiency $\epsilon_r$ (left panel) and the feedback efficiency $\eta$ (right panel) as a function of BH mass for BHs in NEWHORIZON at redshifts $z = 3, 2, 1$, and $0.25$. For BHs in quasar mode, where $\epsilon_r = \epsilon_r^{\text{thin}}$, only one data point is shown. For BHs in radio mode, where $\epsilon_r = f_{\text{att}}\epsilon_r^{\text{thin}}$, any BHs with $\epsilon_r < 10^{-3}$ are shown as limits. The dashed lines in the left plot show the average $\epsilon_r^{\text{thin}}$ and $\epsilon_r$ as a function of mass across all redshifts.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Bolometric BH luminosity function at $z = 0.25$ (left panel) and $z = 2$ (right panel) for all BHs in NEWHORIZON. The luminosity function is shown for three different assumptions: BHs have the spin-dependent radiative efficiency $\epsilon_r$ as extracted from the simulation (solid), all BHs have a fixed radiative efficiency of $\epsilon_r = 0.1$ (dashed), all BHs have a spin-dependent $\epsilon_r = 0.1$ and are assumed to be maximally (dash-dotted) and non-spinning (dotted), respectively. Observations from (Soares & Nemmen 2020) (S20) and (Hopkins, Richards & Hernquist 2007, H07) are shown for comparison.</figcaption></figure>

We caution that the sample of BHs in NEWHORIZON is dominated by intermediate-mass, low-luminosity AGN. Conclusions on the impact of a spin-dependent radiative efficiency on the luminosity function might change significantly for a more representative sample. The consequences of increasing the radiative efficiency for massive BHs might require re-calibration of models elsewhere. For example,

Sijacki et al. (2015) found that using an average radiative efficiency of $\epsilon_r = 0.2$ over-predicted the AGN luminosity function, with $\epsilon_r = 0.1$ in better agreement with observations (but also readily acknowledge that in their model $\epsilon_r$ is degenerate with a numerical feedback efficiency similar to $\epsilon_f$). Overall, a careful study of a large sample of massive BHs with on-the-fly spin evolution will be needed to
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Fraction of time $f_{\text{time}}$ individual BHs spend in quasar (solid) or radio (dashed) mode with a bolometric luminosity above a given luminosity cut. Fractions are measured across cosmic history (left panel) and in the range $z = 0.5–0.25$ (right panel).</figcaption></figure>

evaluate the impact of spin-dependent radiative efficiency on BH-galaxy coevolution.

# 5 BLACK HOLE SPIN OBSERVABILITY

To be observable, BHs have to be sufficiently luminous to meet observational luminosity cuts. In this section, we explore how the observable distribution of BH spins can vary over time due to a combination of evolving BH luminosity and BH accretion state. As discussed briefly in Section 4.1, different observational methods probe BHs in different accretion states.

Looking at the final redshift of the simulation ($z = 0.25$), the potentially observable sample is small: out of a sample of 583 BHs, only six exceed a bolometric luminosity of $10^{39}$ erg s$^{-1}$, with none brighter than $3 \times 10^{43}$ erg s$^{-1}$. This is because our sample of BHs is dominated by intermediate-mass BHs (see fig. 1 and Beckmann et al. 2023 for a full analysis) and the sample of massive BHs is small. All six bright BHs are in jet mode, which is not unexpected as there are only four BHs in quasar mode within the whole sample at this point (see Fig. 1).

However, BHs show strong variability over time, so when exactly their luminosity and accretion state is measured can make a significant difference to their observability, as can be seen in the example timeseries in Section 4.1. We quantify this variability for the five most massive BHs in NEWHORIZON in Fig. 13, which shows the percentage of time a BH would be observable in either quasar mode or radio mode for a given luminosity cut. Fractions are measured for all time (left panel) and since $z < 0.5$ only (right panel). At high luminosity, the probability of observing any of the BHs in NEWHORIZON drops significantly, but all observable BHs will be in quasar mode, so only observable with X-ray-based methods. At lower luminosity ($< 10^{42}$ erg s$^{-1}$), some BHs are much more likely to be in a jetted radio-mode, especially at $z < 0.5$ (e.g. BH 936 and 549) while others have a constant probability of being observed in quasar mode (here around 20 per cent for BHs 166 and 455) but a small probability of being found in radio mode at any reasonable luminosity cut.

It has recently been reported that all massive galaxies with stellar masses $M_{\text{star}} > 10^{11} M_{\odot}$ in the local Universe show radio emission

associated with the AGN (Sabater et al. 2019). At first glance, we struggle to reproduce this result as our most massive BHs (BHs 166 and 455), which are also hosted in our most massive galaxies, spend more time in quasar mode than radio mode at low redshift. At any given time, the fraction of AGN in radio mode in galaxies with stellar masses $> 10^{11} M_{\odot}$ in NEWHORIZON is therefore less than 100 per cent. This could be a sign that BH growth in NEWHORIZON is artificially delayed or reduced, as was also suggested by the fact that all massive BHs in NEWHORIZON lie in the lower region of the $M_{\text{bh}}-\sigma$ relation, and that NEWHORIZON generally fails to reproduce the observed fraction of AGN in dwarf galaxies the local Universe (see Beckmann et al. 2024, for more information). It could also mean that our binary disc model is too abrupt, and a smoother transition between quasar mode and radio mode via a truncated disc, such as in Koudmani et al. (2024) is required, or that radio emitting features driven by AGN are more long-lived than the short-scale variability in their accretion mode. Given that our sample of such massive galaxies, and their BHs is extremely small, we note the importance of using such population-wide observations to constrain BH accretion but post-pone a careful study of its impact to future work with a sample better suited to do population-wide studies.

Shown in Fig. 14 is the population of observable BHs in both radio and quasar mode for a stacked sample of BHs with different minimum bolometric luminosity cuts. All BHs in NEWHORIZON have been sampled 30 times at even time intervals from $z = 0.5$ to 0.25 to account for variability in both accretion state and luminosity. As can be seen in Fig. 14, not surprisingly the distribution of observable BHs shifts towards higher masses with increasing minimum luminosity. At the same time, higher luminosity cuts mean a shift towards BHs in quasar mode, and a decreasing sample size for BHs observable via their radio jets. It also shifts the observed BH spins to higher spin values, with low-spin BHs dropping out of the observable sample. This includes both low-mass BHs with masses $M_{\text{BH}} < 10^5 M_{\odot}$, which can have spins as low as $|a| < 0.2$, and higher mass BHs with masses $M_{\text{BH}} > 10^6 M_{\odot}$ in radio mode that tend to have spins in the range $|a| = 0.4–0.8$. In comparison to observations, we find fewer quasar-mode BHs in this mass and spin range than reported by X-ray-based observations, but we caution that our sample is too small to draw firm conclusions.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Distribution of observable spins from both quasar mode (red triangles) and jet mode (green circles) for different luminosity cuts. All NEWHORIZON BHs are sampled 30 times in the redshift range $z = 0.5–0.25$ to account for variability in both luminosity and accretion state. Shown for comparison are observed BH spins in the local Universe for X-ray-based (empty triangles, Reynolds 2021; Mallick et al. 2022; Sisk-Reynés et al. 2022) and radio-based techniques (empty circles, Soares & Nemmen 2020). No luminosity cut has been applied to observations beyond that of the original observational sample.</figcaption></figure>

Overall, we conclude that the required high luminosity cuts (e.g. Reynolds (2013) only consider sources with a minimum of 0.01 the Eddington luminosity, equivalent to $L > 1.3 \times 10^{42}$ erg s$^{-1}$ for a $10^6$ M$_\odot$ BH) for observations of BH spin bias the distribution of spins towards higher spin values than the full sample of BHs (see also a similar discussion in Brenneman et al. 2011). We caution that the sample of BHs in NEWHORIZON is too small for statistical analysis, and the numbers quoted are highly specific to our individual BHs. This discussion is therefore only intended to highlight the variety of BH evolution histories and their potential impact on the observability of BH spin. Further work is required to understand how the observable distribution of BH spin reflects the underlying BH spin distribution for different observational methods, luminosity cuts and galaxy types surveyed (see e.g. Beckmann et al. 2024). We leave a robust statistical analysis to future work based on a larger volume sample.

## 6 CONCLUSIONS

In this paper, we studied the spin evolution of a sample of BHs in the mass range $10^4–10^8$ M$_\odot$ from their formation to $z = 0.25$ using the cosmological simulation NEWHORIZON. We tracked how BH spin changes through gas accretion, BH-BH merger and BH feedback, and concluded that

(i) BH spin evolution goes through three phases: an early gas-driven spin-up to maximum spin following formation, a merger-induced scattering phase up to a BH mass of $10 \times M_{BH,0}$, where $M_{BH,0}$ is the BH seed mass (here $M_{BH,0} = 10^4$ M$_\odot$). This is followed by a second accretion-driven phase for BHs in the mass range $10^5–10^8$ M$_\odot$.

(ii) During the second accretion-driven phase, when BHs are in the mass range of $10^5–10^8$ M$_\odot$, they can be both spun up (due to efficient
accretion through a thin disc) or spun down (through inefficient accretion through a thick disc, where BH spin energy is extracted to drive BH jets). Due to the high variability of the accretion rate, all BHs alternate between both modes, significantly increasing the scatter in the spin distribution in this mass range where BHs had been reported as predominantly maximally spinning when jet spindown was not taken into account (Dubois et al. 2014b; Bustamante & Springel 2019).

(iii) BH-BH mergers for BHs close to their seed mass (here $10^4 M_\odot$) predominantly spin-up BHs, while spin typically decreases during mergers between more massive BHs, in agreement with work by Dong-Páez et al. (2023). As the BH mass increases the average merger ratio decreases, making BH mergers less important to the long-term spin evolution of the BH. We cannot comment on the previously reported second merger-driven spin-down phase for BHs with masses $M_{BH} > 10^8 M_\odot$ (see Dubois et al. 2014b and Bustamante & Springel 2019) as our BHs are insufficiently massive by the end of the simulation.

(iv) After an early accretion burst after formation, gas accretion on to BHs undergoes a highly chaotic, inefficient phase during which BH spin (and mass) evolves predominantly through BH-BH mergers. At late times, BH spin is typically well aligned with the angular momentum of the accreted gas, allowing for significant accretion-driven spin up (or down, depending on the accretion mode and alignment) over Gyr time periods.

(v) We find no evidence for spin-down through efficient chaotic accretion on to BHs, either in the early Universe or today. At early times, chaotic accretion is either too inefficient to allow for significant BH growth, or too coherent to reduce BH spin. At late times, we find little evidence for persistent chaotic accretion on to massive BHs at all. We do, however, report on one episode of persistent anti-aligned accretion over ~2 Gyr that reduces BH spin.

(vi) At $<\epsilon_r^{\text{thin}}>=0.19$, the average radiative thin disc efficiency for BHs with masses $M_{BH} > 10^5 M_\odot$ is, on average, almost twice as high as the commonly used value of 0.1. For BHs in radio mode, even high spin is unable to compensate for the fact that their accretion discs are inherently radiatively inefficient, making such BHs hard to observe.

(vii) Due to their high spin, massive BHs deposit, on average, between ~ three times (in quasar mode) and ~eight times (in radio mode) more feedback energy into their host galaxy than previously assumed in fixed-spin models.

(viii) The impact of the spin-based radiative efficiency on the luminosity function is negligible for the sample of BHs presented here.

(ix) X-ray-based and radio-based observations of BH spin probe BHs in different feedback modes, with X-ray-based methods probing, on average, a more highly spinning population of BHs than radio-based approaches.

(x) When observing BH spins, higher luminosity cuts lead to higher average observed BH spin values.

Overall, we conclude that BH spin evolution changes significantly over cosmic time, with both accretion and mergers dominating during different phases of the BH's evolution. This has important consequences for their radiative efficiency, and therefore the feedback energy experienced by the BH's host galaxy. Massive galaxies, in particular, have the potential to experience feedback up to twice as strong as canonically modelled if their SMBH are in quasar mode. As has been explored here, the distribution of observed BH spins depends significantly on the observational method. As the sample studied here is dominated by intermediate-mass BHs, with only a

small number of massive BHs up to ~$10^8 M_\odot$, we leave a robust exploration of this effect to future work.

## ACKNOWLEDGEMENTS

For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

This work was granted access to the HPC resources of CINES under the allocations c2016047637, A0020407637, and A0070402192 by Genci, KSC-2017-G2-0003 by KISTI, and as a 'Grand Challenge' project granted by GENCI on the AMD Rome extension of the Joliot Curie supercomputer at TGCC. This research is part of the Spin(e) ANR-13-BS05-0005 (http://cosmicorigin.org), Segal ANR-19-CE31-0017 (http://secular-evolution.org) and Horizon-UK projects. This work has made use of the Infinity cluster on which the simulation was post-processed, hosted by the Institut d'Astrophysique de Paris. We warmly thank S. Rouberol for running it smoothly. The large data transfer was supported by KREONET which is managed and operated by KISTI. RSB would like to thank Newnham College, Cambridge, for financial support. SKY acknowledges support from the Korean National Research Foundation (NRF-2020R1A2C3003769 and NRF-2022R1A6A1A03053472).

## DATA AVAILABILITY

Data is available upon request to the corresponding author.

## REFERENCES

Anglés-Alcázar D., Faucher-Giguère C.-A., Quataert E., Hopkins P. F., Feldmann R., Torrey P., Wetzel A., Kereš D., 2017, MNRAS, 472, L109

Barausse E., Rezzolla L., 2009, ApJ, 704, L40

Bardeen J., 1970, Nature, 226, 64

Bardeen J., Petterson J., 1975, ApJ, 195, L65

Beckmann R. S., Devriendt J., Slyz A., 2019, MNRAS, 483, 3488

Beckmann R. S. et al., 2023, MNRAS, 523, 5610

Beckmann R. S. et al., 2024, MNRAS, 527, 10867

Benson A. J., Babul A., 2009, MNRAS, 397, 1302

Berti E., Volonteri M., 2008, ApJ, 684, 822

Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433

Bollati F., Lupi A., Dotti M., Haardt F., 2024, A&A, 690, A194

Bower R., Schaye J., Frenk C. S., Theuns T., Schaller M., Crain R. A., McAlpine S., 2016, MNRAS, 465, 32

Brenneman L. W. et al., 2011, ApJ, 736, 103

Bustamante S., Springel V., 2019, MNRAS, 490, 4133

Capellupo D. M., Wafflard-Fernandez G., Haggard D., 2017, ApJ, 836, L8

Cenci E., Sala L., Lupi A., Capelo P. R., Dotti M., 2020, MNRAS, 500, 3719

Chabrier G., 2005, in Corbelli E., Palla F., Zinnecker H., eds, Astrophysics and Space Science Library Vol. 327, The Initial Mass Function 50 Years Later. Springer-Verlag, Berlin. p. 41

Daly R. A., 2011, MNRAS, 414, 1253

Daly R. A., 2019, ApJ, 886, 37

Dauser T., Wilms J., Reynolds C. S., Brenneman L. W., 2010, MNRAS, 409, 1534

Done C., Jin C., Middleton M., Ward M., 2013, MNRAS, 434, 1955

Dong-Páez C. A., Volonteri M., Beckmann R. S., Dubois Y., Trebitsch M., Mangiagli A., Vergani S. D., Webb N. A., 2023, A&A, 673, A120

Dotti M., Colpi M., Pallini S., Perego A., Volonteri M., 2013, ApJ, 762, 68

Dubois Y., Volonteri M., Silk J., 2014a, MNRAS, 440, 1590

Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., 2014b, MNRAS, 440, 2333

Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., Teyssier R., 2015, MNRAS, 452, 1502

Dubois Y. et al., 2021, A&A, 651, A109
Fabian A. C., Rees M. J., Stella L., White N. E., 1989, MNRAS, 238, 729

Fanidakis N., Baugh C. M., Benson A. J., Bower R. G., Cole S., Done C., Frenk C. S., 2011, MNRAS, 410, 53

Fiacconi D., Sijacki D., Pringle J. E., 2018, MNRAS, 477, 3807

García J., Kallman T. R., Mushotzky R. F., 2011, ApJ, 731, 131

Griffin A. J., Lacey C. G., Gonzalez-Perez V., Lagos C. d. P., Baugh C. M., Fanidakis N., 2019, MNRAS, 487, 198

Haardt F., Madau P., 1996, ApJ, 461, 20

Habouzit M., Volonteri M., Dubois Y., 2017, MNRAS, 468, 3935

Hawking S. W., 1975, Commun. Math. Phys., 43, 199

Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731 (H07)

Hopkins P. F. et al., 2024, Open J. Astrophys., 7, 19

Izquierdo-Villalba D., Bonoli S., Dotti M., Sesana A., Rosas-Guevara Y., Spinoso D., 2020, MNRAS, 495, 4681

Jiang J. et al., 2019, MNRAS, 489, 3436

Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, MNRAS, 451, 2900

Kimm T., Katz H., Haehnelt M., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 466, stx052

King A. R., Lubow S. H., Ogilvie G. I., Pringle J. E., 2005, MNRAS, 363, 49

King A. R., Pringle J. E., Hofmann J. A., 2008, MNRAS, 385, 1621

Komatsu E. et al., 2011, ApJS, 192, 18

Koudmani S., Somerville R. S., Sijacki D., Bourne M. A., Jiang Y.-F., Profit K., 2024, MNRAS, 532, 60

Laor A., 1991, ApJ, 376, 90

Levine R., Gnedin N. Y., Hamilton A. J. S., 2010, ApJ, 716, 1386

Lodato G., Pringle J. E., 2006, MNRAS, 368, 1196

Lousto C. O., Nakano H., Zlochower Y., Campanelli M., 2010, Phys. Rev. D, 81, 084023

Lousto C. O., Zlochower Y., Dotti M., Volonteri M., 2012, Phys. Rev. D, 85, 084015

Maio U., Dotti M., Petkova M., Perego A., Volonteri M., 2013, ApJ, 767, 37

Mallick L. et al., 2022, MNRAS, 513, 4361

McKinney J. C., Tchekhovskoy A., Blandford R. D., 2012, MNRAS, 423, 3083

Massonneau W., Dubois Y., Volonteri M., Beckmann R. S., 2023a, A&A, 669, A143

Massonneau W., Volonteri M., Dubois Y., Beckmann R. S., 2023b, A&A, 670, A180

Novikov I. D., Thorne K. S., Gordon and Breach, Science Publishers, 1973, Black Holes (Les Astres Occlus). p. 343

Ostriker E. C., 1999, ApJ, 513, 252

Pacucci F., Loeb A., 2021, MNRAS, 509, 1885

Parker M. L., Matzeu G. A., Matthews J. H., Middleton M. J., Dauser T., Jiang J., Joyce A. M., 2022, MNRAS, 513, 551

Peirani S. et al., 2024, A&A, 686, A233

Piotrowska J. M. et al., 2024, Front. Astron. Space Sci., 11, 1324796

Reynolds C. S., 2013, Class. Quantum Gravity, 30, 244004

Reynolds C. S., 2021, ARA&A, 59, 117

Rezzolla L., Barausse E., Dorband E. N., Pollney D., Reisswig C., Seiler J., Husa S., 2007, Phys. Rev. D, 78, 044002

Riaz S., Ayzenberg D., Bambi C., Nampalliwar S., 2020, ApJ, 895, 61

Risaliti G. et al., 2013, Nature, 494, 449

Sabater J. et al., 2019, A&A, 622, A17

Sala L., Valentini M., Biffi V., Dolag K., 2024, A&A, 685, A92

Sesana A., Barausse E., Dotti M., Rossi E. M., 2014, ApJ, 794, 104

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Shapiro S. L., 2005, ApJ, 620, 59

Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G. F., Nelson D., Hernquist L., 2015, MNRAS, 452, 575

Sisk-Reynés J., Reynolds C. S., Matthews J. H., Smith R. N., 2022, MNRAS, 514, 2568

Sądowski A., Bursa M., Abramowicz M., Kluźniak W., Lasota J.-P., Moderski R., Safarzadeh M., 2011, A&A, 532, A41

Soares G., Nemmen R., 2020, MNRAS, 495, 981

Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253

Teyssier R., 2002, A&A, 385, 337

Thorne K., 1974, ApJ, 191, 507

Trakhtenbrot B., 2014, ApJ, 789, L9

Trebitsch M., Blaizot J., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 470, 224

Trebitsch M. et al., 2021, A&A, 653, A154

Vasudevan R. V., Fabian A. C., Reynolds C. S., Aird J., Dauser T., Gallo L. C., 2016, MNRAS, 458, 2012

Volonteri M., Madau P., Quataert E., Rees M. J., 2005, ApJ, 620, 69

Volonteri M., Sikora M., Lasota J.-P., Merloni A., 2013, ApJ, 775, 94

Volonteri M. et al., 2020, MNRAS, 498, 2219

Walton D. J. et al., 2019, MNRAS, 484, 2544

You B., Cao X.-W., Yuan Y.-F., 2016, Res. Astron. Astrophys., 16, 003

Zubovas K., King A., 2019, MNRAS, 489, 1373

Zubovas K., King A., 2021, MNRAS, 501, 4289

This paper has been typeset from a TEX/LATEX file prepared by the author.