11institutetext: Universitäts-Sternwarte, Fakultät für Physik, Ludwig-
Maximilians-Universität München, Scheinerstr. 1, D-81679 München, Germany
22institutetext: Astronomy Unit, Department of Physics, University of Trieste,
via Tiepolo 11, I-34131 Trieste, Italy 33institutetext: INAF – Osservatorio
Astronomico di Trieste, via Tiepolo 11, I-34131 Trieste, Italy
44institutetext: ICSC - Italian Research Center on High Performance Computing,
Big Data and Quantum Computing 55institutetext: Max-Planck-Institut für
Astrophysik, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany

# Supermassive black hole spin evolution in cosmological simulations with
OpenGadget3

Luca Sala email<EMAIL_ADDRESS>Milena Valentini,,, 22331144 Veronica
Biffi 33 Klaus Dolag, 1155

(Received ; accepted .)

###### Abstract

Context. Mass and spin of massive black holes (BHs) at the centre of galaxies
evolve due to gas accretion and mergers with other BHs. Besides affecting e.g.
the evolution of relativistic jets, the BH spin determines the efficiency with
which the BH radiates energy.

Aims. Using cosmological, hydrodynamical simulations, we investigate the
evolution of the BH spin across cosmic time and its role in controlling the
joint growth of supermassive BHs and their host galaxies.

Methods. We implement a sub-resolution prescription that models the BH spin,
accounting for both BH coalescence and misaligned accretion through a
geometrically thin, optically thick disc. We investigate how BH spin evolves
in two idealised setups, in zoomed-in simulations, and in a cosmological
volume. The latter simulation allows us to retrieve statistically robust
results as for the evolution and distribution of BH spins as a function of BH
properties.

Results. We find that BHs with $M_{\rm BH}\lesssim 2\times
10^{7}\leavevmode\nobreak\ {\rm M}_{\sun}$ grow through gas accretion,
occurring mostly in a coherent fashion that favours spin-up. Above $M_{\rm
BH}\gtrsim 2\times 10^{7}\leavevmode\nobreak\ {\rm M}_{\sun}$ the gas angular
momentum directions of subsequent accretion episodes are often uncorrelated
with each other. The probability of counter-rotating accretion and hence spin-
down increases with BH mass. In the latter mass regime, BH coalescence plays
an important role. The spin magnitude displays a wide variety of histories,
depending on the dynamical state of the gas feeding the BH and the relative
contribution of mergers and gas accretion. As a result of their combined
effect, we observe a broad range of values of the spin magnitude at the high-
mass end. Reorientation of the BH spin direction occurs on short timescales
($\lesssim 10$ Myr) only during highly-accreting phases ($f_{\rm Edd}\gtrsim
0.1$). Our predictions for the distributions of BH spin and spin-dependent
radiative efficiency as a function of BH mass are in very good agreement with
observations.

###### Key Words.:

accretion, accretion discs – black hole physics – galaxies: nuclei – methods:
numerical

##  1 Introduction

Supermassive black holes (BHs) ($10^{6}\lesssim M_{\rm
BH}/\leavevmode\nobreak\ {\rm M}_{\sun}\lesssim 10^{10}$) at the centre of
massive galaxies are observed as active galactic nuclei (AGN), if they accrete
material from and release energy in their proximity. This process, dubbed AGN
feedback, is able to significantly affect the surroundings of the BHs and it
is thought to play a significant role in the evolution of massive galaxies
(e.g. Benson et al. 2003; Di Matteo et al. 2005). Observationally, feedback
has been broadly categorized in two modes, depending on the main channel of
energy release (Churazov et al. 2005; Sijacki et al. 2007; Fabian 2012). In
the so-called quasar mode, operating close to the Eddington limit, most of the
energy is released in the form of radiation and winds (Silk & Rees 1998;
Fabian 1999; Harrison et al. 2018). This mode is prevalent at earlier epochs
(e.g. Croton et al. 2006). A second mode occurs at highly sub-Eddington
accretion rates, when energy is channelled mostly into powerful relativistic
jets (McNamara & Nulsen 2007; Cattaneo et al. 2009; McNamara & Nulsen 2012).
This mode is often referred to as maintenance or radio mode.

Rotating BHs have been advocated to have a fundamental role in the formation
of jets (Davis & Tchekhovskoy 2020) and related feedback (McNamara & Nulsen
2012). The spin $J_{\rm BH}$ of a BH with mass $M_{\rm BH}$ can range from 0
to the maximal value $J_{\rm max}=GM_{\rm BH}^{2}/c$. Thus, a BH is
characterised by the dimensionless spin parameter $a=J_{\rm BH}/J_{\rm max}$.
A maximally-spinning $10^{9}\leavevmode\nobreak\ {\rm M}_{\sun}$ BH can
potentially provide $\sim 10^{62}$ erg of rotational energy (McNamara & Nulsen
2012) that could be injected into the surroundings and offset cooling, once a
mechanism to extract this energy is in action. Theoretically, a viable
mechanism has been suggested by Blandford & Znajek (1977) (see also Lasota et
al. 2014). This process requires poloidal magnetic fields anchored to an
accretion flow, surrounding a rotating BH. The accreting matter drags the
poloidal field into the BH ergosphere, where the BH rotation twists the field
lines, generates a toroidal component of the field and exerts an effective
pressure that accelerates the plasma (Tchekhovskoy 2015). Thus, the power of
the jet is affected by the BH spin magnitude and the magnetic field flux
threading the BH event horizon. The latter can be maximised in the so-called
magnetically-arrested disc state. In this configuration, magnetic flux is
accumulated near the event horizon, until it becomes dynamically relevant. In
this state, the ratio between jet power and accretion power (i.e. $\eta=P_{\rm
jet}/(\dot{M}c^{2})$) can be larger than one (Tchekhovskoy et al. 2011;
McKinney et al. 2012; Narayan et al. 2022). This is possible by extracting
rotational energy from the BH. Launching powerful jets requires the presence
of an accretion flow that acts as an anchor for the magnetic field and drags
it inwards. If a jet is created, it is possible to tap into the rotational
energy of the BH. Conversely, it is possible to power an AGN by accretion
alone (without extracting spin energy) but in this case the maximum power
output can only reach few tens percent of the accretion power.

Observationally, a connection between the BH spin and the jet power has been
proposed to explain the observed radio-loudness (radio-to-optical flux density
ratio, Sikora et al. 2007) in samples of AGN (e.g. Sikora et al. 2007;
Martínez-Sansigre & Rawlings 2011). Ghisellini et al. (2010) find that the
power of the most powerful jets in their sample of blazars (see Giommi et al.
2012, for a definition) can exceed the luminosity of the accretion disc by a
factor of $\sim$10\. They argue that extraction of rotational energy may
provide the additional power required. They also argue that if this is the
case, the observed correlation between $P_{\rm jet}$ and accretion disc
luminosity may arise from a link between the accretion and rotational energy
extraction processes (as discussed above). Finally, even though instrumental
capabilities prevent a definite confirmation so far, future polarimetry
measurements with the Event Horizon Telescope (EHT) could provide
observational support to the process of energy extraction by jets (Chael et
al. 2023).

The direction of the BH spin also plays a key role in jet formation. Spinning
BHs create an outward Poynting flux in the direction of the BH spin, through
the Blandford-Znajek mechanism (Hawley & Krolik 2006; Nakamura et al. 2008).
Several observations of nearby clusters exhibit the presence of multiple pairs
of jet-inflated cavities, that have distinct angular orientations relative to
the central galaxy (e.g. Forman et al. 2007; David et al. 2009; Sanders et al.
2009; Fabian et al. 2011). This suggests multiple generations of jets
propagating along different directions. The observations may be interpreted as
spin re-orientation, provided that processes close to the BH manage to modify
the BH spin and thus the launching direction of the jet (e.g. accretion
proceeding on different planes across subsequent events). A few numerical
works demostrated that it is possible to recover the observed morphologies by
assuming a re-orienting jet (Cielo et al. 2018; Horton et al. 2020; Lalakos et
al. 2022).

BH spins evolve across cosmic time due to accretion and mergers. Observational
estimates of the BH spin parameter are therefore crucial to gain insight into
its evolution. Several methods are adopted in literature to infer the BH spin
of AGN (e.g. X-ray reflection methods or thermal continuum fitting, see
Reynolds 2021, for a review). Most of the measurements published to date are
based on X-ray reflection spectroscopy based on Fe K lines. Recent works
developed a new class of high-density disc models (Jiang et al. 2019, 2022),
used to estimate spin in low-mass BHs by modelling the X-ray relativistic
reflection continuum (Mallick et al. 2022). The observed jet powers can also
be used to determine the spin (Daly 2009, 2011, 2019, 2021). This method is
particularly useful because it can be applied to large samples of AGN
exhibiting jets. Lastly, interferometric observations with the EHT might
provide constraints on the spin of M87∗, from the circularity of the shadow
(Broderick et al. 2022) or from the phase-twisting of light propagating near
the BH (Tamburini et al. 2020). The former method was only able to infer a
clockwise rotation with the current observational capabilities (i.e. spin
vector pointing away from Earth). The latter provided an estimate of
$a=0.90\pm 0.05$ and an angle with respect to the line of sight
$i=163^{\circ}\pm 2^{\circ}$.

From the theoretical standpoint, several analytical studies have been focused
on investigating the mechanism driving spin evolution, that is likely linked
to the coupling between the accretion disc and the BH spin (Papaloizou &
Pringle 1983; Pringle 1981; Scheuer & Feiler 1996; Natarajan & Pringle 1998;
Natarajan & Armitage 1999; King et al. 2005; Martin et al. 2007; Perego et al.
2009). If a geometrically thin, optically thick disc is misaligned with
respect to the BH spin rotation axis, the innermost part aligns with the BH
spin and is connected to the unperturbed outermost part through a smooth warp
(Bardeen-Petterson effect). The flow through this region exerts a torque on
the BH spin, modifying its direction. Recently, full general-relativistic
magneto-hydrodynamic simulations performed by Liska et al. (2019) were able to
reproduce this effect in magnetised thin discs. The BH spin magnitude changes
when matter is accreted at the innermost stable circular orbit (ISCO).

Further works have focused on building upon the analytical theory to
understand the evolution of BH spin as a response of accretion and mergers, as
well as the relative contribution between the two channels. A number of works
adopted a semi-analytical treatment in a hierarchical cosmological context,
with various recipes for the drivers of gas accretion and its effect on spin
evolution (Volonteri et al. 2005; Volonteri & Rees 2005; Berti & Volonteri
2008; Lagos et al. 2009; Fanidakis et al. 2011; Izquierdo-Villalba et al.
2020; Griffin et al. 2019). Other studies focused on simulating evolutionary
histories of BHs under the assumption of a fixed fuelling angular momentum
distribution and analysed how it affects the distribution of BH spins
(Volonteri et al. 2007; King et al. 2008; Dotti et al. 2013; Zhang & Lu 2019).
Going a step further, some works developed models suitable to be used on-the-
fly in hydrodynamical simulations (Maio et al. 2013; Fiacconi et al. 2018), to
study isolated systems at high resolution. The model by Fiacconi et al. (2018)
was also coupled to novel feedback recipes for winds (Cenci et al. 2020; Sala
et al. 2021) and jets (Talbot et al. 2021). Recently, using the same model,
Bollati et al. (2023) studied the dynamics of massive BH binaries in the
presence of spin-dependent radiative feedback, whereas Talbot et al. (2023)
focused on the effect of jets in gas-rich galaxy mergers. While all the above
works assumed a thin accretion disc, Huško et al. (2022) developed a spin
evolution model that assumes a thick disc and applied it to study jet feedback
in galaxy clusters. Dubois et al. (2014a) made a further step forward and
developed a model suited to cosmological, hydrodynamical simulations, built
upon Volonteri et al. (2007) and King et al. (2005, 2008). The model was used
to perform zoom-in simulations (Dubois et al. 2014a) and later updated to
include jets, with BH spin-dependent direction and power (Beckmann et al.
2019; Dubois et al. 2021; Dong-Páez et al. 2023).

Comparatively fewer works have been dedicated to the study of spin evolution
in a cosmological box and using a large statistical sample of BHs. Dubois et
al. (2014b) run their spin evolution model, although only in post-processing,
on the output of the Horizon-AGN simulation. Bustamante & Springel (2019)
performed a set of simulations of a cosmological volume, evolving the spin on-
the-fly with a similar model, although with a few differences at the sub-
resolution level with respect to our implementation (see Sec. 2.2 and 6). In
this work, we focus on studying supermassive BH spin evolution with a sub-
resolution model that follows the latter two works. Our aim is to take full
advantage from the statistically significant simulated population of BHs in a
cosmological box. Using information provided by cosmological, hydrodynamical
simulations, we aim to analyse in detail the coupled evolution of the spin
direction and magnitude, the effect of the dynamical state of the gas fuelling
the BHs and the impact of mergers.

This paper is organised as follows. Sec. 2 describes the theoretical basis of
our model and its implementation. Sec. 3 presents the suite of simulations we
used to test our model and statistically study BH spins. Sec. 4 presents our
simulation results, whereas in Sec. 5 we discuss their implications. In Sec. 6
we compare our work with previous studies in literature. In Sec. 7 we
summarise and conclude.

##  2 The model

In this Section, we present our sub-resolution model for spin evolution due to
gas accretion onto BHs and mergers and its coupling to the resolved scales. We
implement our model in the TreePM code OpenGadget3 (see Groth et al. 2023),
descendant of a non-public evolution of the Gadget-3 code (originally from
Springel 2005). The code features a modern smoothed particle hydrodynamics
(SPH) description (Beck et al. 2016) and adopts a bias-corrected, sixth-order
Wendland kernel (Dehnen & Aly 2012) with 295 neighbours.

###  2.1 BH treatment

Our sub-resolution model is built within the broader BH physics module for
cosmological simulations originally described in Springel et al. (2005b), with
modifications as in Hirschmann et al. (2014) and Steinborn et al. (2015).

Each BH in our simulations is treated as a collisionless sink particle,
coupled to other particles only by gravity. We associate an SPH kernel to
every BH (with smoothing length $h_{\rm BH}$) in a similar fashion as for the
gas particles, with the same number of neighbours. The gravitational force on
the BH is softened using a Plummer-equivalent softening length
$\varepsilon_{\rm BH}$. Its value is specific for each simulation and stated
in the corresponding subsections and in Table 2. To select haloes where to
seed new BH particles, we apply a FoF (Friends-of-Friends) algorithm to the
stellar particles alone. In this way, we identify stellar bulges. We adopt a
linking length of $l=0.05$ and require that a halo has a stellar mass
corresponding to $M_{\rm*,seed}$ to host a new BH seed. We also demand that
the halo has a stellar over dark matter (DM) mass fraction $f_{\rm*}>0.05$ and
a gas over stellar mass fraction $f_{\rm gas}>0.1$, to ensure that the
identified stellar bodies are newly-forming galaxies and not tidally stripped
stellar remains. The seed is created at the position of the star particle with
the largest binding energy within the group, with an initial mass $M_{\rm
BH,seed}$. The values of $M_{\rm BH,seed}$ and $M_{\rm*,seed}$ for each
simulation are stated in the corresponding subsections and in Table 2.
Alternatively, BHs with specific properties can be inserted in the initial
conditions (ICs) to study a specific configuration of an idealised system
(Sec. 3.1 and 3.2).

We do not apply any pinning of the BHs to the position of the potential
minimum within $h_{\rm BH}$. Besides, we do not apply any drag force from the
continuous gas accretion onto the BH and keep position and velocity of the
central BH when BHs merge. We allow BHs to merge if the following criteria are
met: (i) their relative velocity is smaller than 0.5$c_{s}$, where $c_{s}$ is
the sound speed of the gas, averaged in a kernel-weighted fashion; (ii) their
distance is smaller than 5 times $\varepsilon_{\rm BH}$; (iii) the binary
binding energy is smaller than $0.5c_{s}^{2}$ (see Hirschmann et al. 2014, for
further details).

Besides mergers, BHs also increase their mass111In our code BHs are
characterised by two masses, a dynamical mass used to compute gravitational
interactions and a physical mass used in all the BH sub-resolution models. by
accreting gas. The mass inflow rate onto the BH is computed using the
Bondi–Hoyle–Lyttleton prescription (hereafter Bondi; Hoyle & Lyttleton 1939;
Bondi & Hoyle 1944; Bondi 1952):

$\dot{M}_{\rm B}=\frac{4\pi\alpha_{\rm acc}G^{2}M^{2}_{\rm
BH}\langle\rho\rangle}{(\langle c_{\rm s}\rangle^{2}+\langle
v\rangle^{2})^{3/2}},$ (1)

where $\rho$ is the gas density, $v$ is the velocity of the BH with respect to
the gas and $G$ is the gravitational constant. We compute the properties
marked with $\langle\cdot\rangle$ as a kernel-weighted average over the BH
neighbouring particles. We also enforce a maximum radius on $h_{\rm BH}$ for
the computation, $r_{\rm acc}=100$ kpc. $\alpha_{\rm acc}$ is a dimensionless
parameter that is often adopted to account for the detailed structure of the
interstellar medium (ISM) in the vicinity of the BH (e.g. Springel et al.
(2005a); Booth & Schaye (2009); but see e.g. Valentini et al. (2020)),
typically unresolved in cosmological simulations because of the limited
spatial resolution ($\sim$kpc). In our code we compute two different Bondi
rates using $\alpha_{\rm acc,hot}=10$ and $\alpha_{\rm acc,cold}=100$ (unless
stated otherwise) for the hot and cold gas phases222The hot phase is composed
by particles with a temperature $T>2\times 10^{5}$ K. The cold phase is
composed by star-forming particles (see Sec. 3) and gas particles that have a
temperature $T<2\times 10^{5}$ K. respectively. We assume that the accretion
rate onto the BH is given by

$\dot{M}_{\mathrm{BH}}=\min\left(\dot{M}_{\mathrm{B},\mathrm{hot}}+\dot{M}_{\mathrm{B},\mathrm{cold}},\dot{M}_{\mathrm{Edd}}\right)$
(2)

where $\dot{M}_{\rm Edd}=4\pi GM_{\rm BH}m_{\rm p}/(\sigma_{\rm
T}\epsilon_{\rm r}c)$ is the Eddington (1916) accretion rate, $m_{\rm p}$ the
proton mass, $\sigma_{\rm T}$ the Thomson cross section, $c$ is the speed of
light and $\epsilon_{\rm r}$ is the radiative efficiency. The latter depends
on the BH dimensionless spin parameter in our simulations, as explained in
Sec. 2.2.1.

###  2.2 Spin evolution algorithm

Our model is designed for simulations whose spatial resolution ranges from few
tens of pc to a few kpc (see Sec. 3), several orders of magnitudes larger than
the physical size of an accretion disc around a BH ($<1$ pc, Frank et al.
2002). In the approach often adopted in cosmological simulations, an accretion
disc is not included in the modelling (e.g. Springel et al. 2005b; Booth &
Schaye 2009). In our implementation, we introduce a sub-grid accretion disc as
an intermediate step in the mass transfer between the resolved scales and the
BH. We then assume that the mass transfer rate from the resolved scales onto
the accretion disc is equal to the mass rate from the accretion disc onto the
BH. The mass rate is defined by Eqs. (1) and (2). We also assume that the gas
maintains the angular momentum direction it has at the resolved scales. The
inclusion of an accretion disc is necessary to model the physical effects that
modify the spin due to gas accretion, as explained in the following.

####  2.2.1 Accretion disc and spin evolution

We define the BH angular momentum vector $\bm{J}_{\rm BH}$ as

$\bm{J}_{\rm BH}=J_{\rm BH}\cdot\bm{j}_{\rm BH}=aJ_{\rm max}\bm{j}_{\rm BH},$
(3)

where $\bm{j}_{\rm BH}$ is the unit vector encoding its direction and $J_{\rm
BH}$ is its magnitude. We define $0\leq a\leq 1$. In what follows, a negative
sign encodes counter-rotating conditions on a BH with spin parameter $a$.
$\bm{J}_{\rm BH}$ evolves because of its interaction with the distribution of
matter in a surrounding accretion disc, whose angular momentum is misaligned
with respect to $\bm{J}_{\rm BH}$. Our model assumes an optically thick,
geometrically thin Shakura & Sunyaev (1973) accretion disc.

If we start from a completely flat disc surrounding a spinning BH, Lense-
Thirring precession forces the fluid elements close to the BH to precess and
induces them to rotate in the BH equatorial plane. If viscosity is
sufficiently high, the shear stresses propagate the perturbation diffusively
until a warped steady state of the disc is reached (Bardeen & Petterson 1975;
Martin et al. 2007). The largest deviation from a flat profile – i.e. disc
annuli where the gas angular momentum is misaligned with respect to both the
unperturbed outer region of the disc and the BH spin – occurs at the warp
radius. The latter is defined as the radius $R_{\rm w}$ at which the Lense-
Thirring precession timescale $t_{\rm LT}=R_{\rm w}^{3}c^{2}/(2GJ_{\rm BH})$
(Wilkins 1972) is equal to the warp propagation timescale $t_{\rm w}\sim
R_{\rm w}^{2}/\nu_{2}$ (Pringle 1981; Perego et al. 2009). $\nu_{2}$ is the
vertical shear viscosity governing the propagation of vertical perturbations
(e.g. Volonteri et al. 2007; Perego et al. 2009; Dotti et al. 2013; Dubois et
al. 2014b). Following Dubois et al. (2014b) we compute

$\displaystyle\frac{R_{\text{w}}}{R_{\mathrm{BH}}}\simeq 4\times
10^{2}a^{5/8}M_{\mathrm{BH},8}^{1/8}\left(\frac{f_{\rm
Edd}}{\epsilon_{\mathrm{r},01}}\right)^{-1/4}\left(\frac{\nu_{2}/\nu_{1}}{85}\right)^{-5/8}\alpha_{\nu_{1},01}^{-1/2}.$
(4)

where $M_{\mathrm{BH},8}=M_{\rm BH}/(10^{8}\leavevmode\nobreak\ {\rm
M}_{\sun})$, $\epsilon_{\mathrm{r},01}=\epsilon_{\mathrm{r}}/0.1$ and we
define the Eddington ratio $f_{\rm Edd}=\dot{M}_{\rm BH}/\dot{M}_{\rm Edd}$.
$\nu_{1}$ is the radial shear viscosity, responsible for the radial inward
drift of gas across the disc. We also assume that the viscosity
$\alpha$-parameter (Shakura & Sunyaev 1973) is
$\alpha_{\nu_{1},01}=\alpha_{\nu_{1}}/0.1\equiv 1$, following King et al.
(2005), and $\nu_{\rm 2}/\nu_{\rm
1}=2(1+7\alpha_{\rm\nu_{1}})/(4+\alpha^{2}_{\rm\nu_{1}})/\alpha^{2}_{\rm\nu_{1}}$,
following Ogilvie (1999), leading to a fiducial value for the ratio $\nu_{\rm
2}/\nu_{\rm 1}=85$.

The flow of matter through the warp is responsible to exert a torque that
modifies only the direction of the BH spin. Indeed, matter flowing from the
unperturbed outermost region to the aligned innermost region through the warp
changes its angular momentum direction, therefore forcing the BH spin
orientation to change and ensure conservation of the total angular momentum
(King et al. 2005; Dotti et al. 2013). Since the torque acting to modify the
direction of the BH spin is produced by matter flowing through the warped
region at $R_{\rm w}$, we adopt the same assumption as Volonteri et al. (2007)
and Dubois et al. (2014b) and define $\bm{J}_{\rm d}=J_{\rm d}\bm{j}_{\rm d}$
as the angular momentum of the disc within $R_{\rm w}$. King et al. (2005)
showed that starting from a configuration where $\bm{J}_{\rm BH}$ and
$\bm{J}_{\rm d}$ are initially misaligned, the torque resulting from the
warped configuration always leads the BH spin to align with the total angular
momentum $\bm{J}_{\rm tot}$ of the system disc+BH

$\displaystyle\bm{J}_{\rm tot}$ $\displaystyle=\bm{J}_{\rm BH}+\bm{J}_{\rm
d}.$ (5)

Moreover, the torque acts dissipatively on the disc, whose angular momentum
ends up either aligned or counter-aligned with the BH spin, depending on the
ratio between the angular momentum magnitudes. Namely, the counter-aligned
configuration is possible if and only if

$\cos{\theta_{\rm BH-d}}<-\frac{J_{\rm d}}{2J_{\rm BH}}$ (6)

where $\theta_{\rm BH-d}$ is the angle subtended by the initial angular
momentum directions $\bm{j}_{\rm BH}$ and $\bm{j}_{\rm d}$, i.e.
$\cos{\theta_{\rm BH-d}}=\bm{j}_{\rm BH}\cdot\bm{j}_{\rm d}$. Only the
(counter-) aligned innermost – i.e. within $R_{\rm w}$ – region can
effectively transfer its angular momentum to the BH, once the matter enclosed
within this region is eventually accreted (Volonteri et al. 2007). Therefore,
$R_{\rm w}$ is also the relevant radial scale to estimate the variation of the
BH spin magnitude as a result of the accretion of this innermost part. Since
the radial shear viscosity $\nu_{1}$ is responsible for the inward drift of
matter across the disc, the region within $R_{\rm w}$ is accreted on a
timescale $t_{\rm\nu_{1}}(R_{\rm w})\sim R_{\rm w}^{2}/\nu_{1}$ (King et al.
2005; Perego et al. 2009; Dubois et al. 2014b). We follow Dubois et al.
(2014b) and compute

$\displaystyle t_{\rm\nu_{1}}(R_{\rm w})\sim 3.4\times
10^{5}a^{7/8}M_{\mathrm{BH},8}^{11/8}\left(\frac{f_{\rm
Edd}}{\epsilon_{\mathrm{r},01}}\right)^{-3/4}\left(\frac{\nu_{2}/\nu_{1}}{85}\right)^{-7/8}\alpha_{\mathrm{\nu_{1}{}},01}^{-3/2}\mathrm{yr}.$
(7)

Note that $\nu_{\rm 1}\ll\nu_{\rm 2}$ and the warp propagation timescale
$t_{\rm w}$ is thus much shorter than the radial drift timescale
$t_{\rm\nu_{1}}$. As a consequence, the timescale on which the warp forms is
much shorter than that over which the spin magnitude and direction change. We
then estimate the mass enclosed within the warped region $M_{\rm d}$ as

$M_{\rm d}\simeq\dot{M}_{\rm BH}t_{\rm\nu_{1}}(R_{\rm w}).$ (8)

where $\dot{M}_{\rm BH}$ is given by Eq. (2). When $M_{\rm d}$ is accreted
onto the BH, the BH spin magnitude $a$ changes due to the accretion of its
ISCO angular momentum, according to the expression derived in Bardeen &
Petterson (1975)

$a_{\rm f}=\frac{1}{3}\frac{r^{1/2}_{\rm isco}}{M_{\rm
ratio}}\left[4-\left(3\frac{r_{\rm isco}}{M^{2}_{\rm
ratio}}-2\right)^{1/2}\right],$ (9)

where

$M_{\rm ratio}=\frac{M_{\rm BH}^{i}+M_{\rm d}(1-\epsilon_{\rm r})}{M_{\rm
BH}^{i}}.$ (10)

and, normalising to the gravitational radius $R_{\rm g}=R_{\rm BH}/2=GM_{\rm
BH}/c^{2}$,

$r_{\rm isco}=R_{\rm isco}/R_{\rm
g}=3+Z_{2}\pm[(3-Z_{1})(3+Z_{1}+2Z_{2})]^{1/2},$ (11)

where the positive (negative) sign is for counter(co)-rotating orbits. $Z_{1}$
and $Z_{2}$ are functions of the BH dimensionless spin parameter only:

$Z_{1}=1+(1-a^{2})^{1/3}[(1+a)^{1/3}+(1-a)^{1/3}],$ (12) $\displaystyle
Z_{2}=(3a^{2}+Z_{1}^{2})^{1/2}.$ (13)

$R_{\rm isco}$ ranges from 1 to 9 $R_{\rm g}$, for co- and counter-rotating
orbits on a maximally spinning BH respectively, as illustrated in the bottom
panel of Fig. 1. Physical sizes are therefore between $\sim 5\times
10^{-6}M_{\rm BH,8}\text{ pc }$ and $\sim 5\times 10^{-5}M_{\rm BH,8}\text{ pc
}$. The mass accreted is also corrected for the radiated energy, where

$\epsilon_{\rm r}=1-\sqrt{1-\frac{2}{3r_{\rm isco}}}.$ (14)

The top panel of Fig. 1 shows the value of $\epsilon_{\rm r}$ as a function of
$a$.

Figure 1: Radiative efficiency – Eq. (14) – (top panel) and radius of the
innermost stable circular orbit – Eq. (11) – (bottom panel), as a function of
the BH dimensionless spin parameter. Some reference values of these quantities
are also highlighted: $a=-1$, for a counter-rotating orbit around a maximally
spinning BH; $a=0$, for a non-spinning BH; $a=0.998$, for the maximum spin
allowed in our simulations.

Finally, as mentioned before, we estimate the angular momentum $J_{\rm d}$ of
the accreted disc region that determines the spin variation in magnitude and
direction (and is also used to evaluate condition 6) as

$J_{\rm d}\simeq M_{\rm d}(R_{\rm w})\omega_{\rm k}(R_{\rm w})R_{\rm
w}^{2}=M_{\rm d}(R_{\rm w})(GM_{\rm BH}R_{\rm w})^{1/2},$ (15)

where $\omega_{\rm k}=\sqrt{(GM_{\rm BH}/R^{3})}$ is the keplerian angular
frequency. While $R_{\rm isco}\sim 1-9R_{\rm g}$, the warped region is at
hundreds of gravitational radii (see Eq. 4), thus a few orders of magnitude
larger. Therefore, the BH spin direction changes on a shorter timescale than
its magnitude (Scheuer & Feiler 1996; Perego et al. 2009; Dotti et al. 2013).

####  2.2.2 BH spin update iteration

Figure 2: Schematic of the steps that compose a single accretion episode. The
vector schemes in the upper part of the figure represent the initial and final
configurations of the angular momenta, in a case similar to the one shown in
Fig. 1b by King et al. (2005). Left: an accretion disc settles around the BH,
in a misaligned configuration. Center: a warp develops, and the innermost part
is forced to rotate in the BH equatorial plane and either co- or counter-
align. Right: the BH spin changes in magnitude when gas is accreted at the
innermost stable orbit.

The algorithm that we adopt to update BH mass and spin in our simulations
models the above-mentioned processes, and makes sure that the BH mass grows
consistently with the rate given by Eq. (2). Following Dubois et al. (2014b),
BH mass and spin are updated in iterations (hereafter accretion episodes)
composed of the following steps:

  1. 1.

we assume that an accretion disc forms around the BH, characterised by the
accretion rate given by Eq. 2. We further assume that the initial angular
momentum direction of the disc $\bm{j}_{\rm d}$ is set by the resolved scales
of the simulation, i.e. parallel to the direction of the angular momentum of
the gas within the BH kernel $\bm{j}_{\rm g}$:

$\bm{j}_{\rm d}=\bm{j}_{\rm g}\equiv\bm{L}_{\rm BH,kernel}/|\bm{L}_{\rm
BH,kernel}|,$ (16)

where

$\bm{L}_{\rm BH,kernel}=\sum_{j}m_{\rm j}(\bm{r}_{\rm j}-\bm{r}_{\rm
BH})\times(\bm{v}_{\rm j}-\bm{v}_{\rm BH})w(\bm{r}_{\rm j}-\bm{r}_{\rm
BH},h_{\rm BH}).$ (17)

$\bm{L}_{\rm BH,kernel}$ is computed considering only the cold gas particles –
i.e. the component that is able to settle into an accretion disc. $w$ is the
dimensionless SPH kernel function. This is the initial misaligned
configuration illustrated on the left in Fig. 2.

  2. 2.

we compute $R_{\rm w}$ using Eq. 4. This defines the region of the disc
(marked in blue in Fig. 2) that exerts the alignment torque on the BH spin and
whose gas is eventually accreted by the end of the accretion episode.

  3. 3.

we compute $M_{\rm d}$ and $J_{\rm d}$, using Eq. 8 and 15 respectively. This
defines the total angular momentum of the accretion episode $\bm{J}_{\rm
tot}=\bm{J}_{\rm BH}+\bm{J}_{\rm d}=\bm{J}_{\rm BH}+{J}_{\rm d}\bm{j}_{\rm d}$
(black solid vector in Fig. 2). Note that we are assuming that the warped
distribution that defines the innermost aligned region of the disc (central
panel in Fig. 2) develops on a timescale shorter than those over which the BH
spin direction and magnitude change, as explained in Sec. 2.2.1.

  4. 4.

we establish whether the innermost part of the disc is co- or counter-rotating
using Eq. (6), by computing

$\displaystyle\frac{J_{\mathrm{d}}}{2J_{\mathrm{BH}}}$
$\displaystyle\simeq\frac{M_{\mathrm{d}}\left(R_{\text{w
}}\right)}{aM_{\mathrm{BH}}}\left(\frac{R_{\text{w
}}}{R_{\mathrm{g}}}\right)^{1/2}$ (18) $\displaystyle\sim 6.8\times
10^{-2}a^{3/16}M_{\mathrm{BH},8}^{23/16}\left(\frac{f_{\rm
Edd}}{\epsilon_{\mathrm{r},01}}\right)^{1/8}\left(\frac{\nu_{2}/\nu_{1}}{85}\right)^{-19/16}\alpha_{\nu_{\rm
1},01}^{-7/4}.$ (19)

In case of counter-rotating conditions, $r_{\rm isco}$ and $\epsilon_{\rm r}$
are computed with a negative sign in front of $a$.

  5. 5.

the final BH spin direction changes as a result of the alignment torque and
ends up parallel to the total angular momentum (vector scheme on the right in
Fig. 2), i.e. $\bm{J}_{\rm BH}^{f}\parallel\bm{J}_{\rm tot}$.

  6. 6.

the disc within $R_{\rm w}$ is consumed and the BH spin magnitude changes
according to Eq. (9) (right panel in Fig. 2), with $M_{\rm d}$ computed in
step 3. Counter-alignment is taken into consideration as described in step 4.
We also cap the spin parameter to 0.998, which is the maximum spin allowed if
photon trapping is assumed (Thorne 1974). The BH mass is increased by $\Delta
M_{\rm BH}=M_{\rm d}(1-\epsilon_{\rm r})$.

We stress that it is possible to update the BH magnitude and direction
separately because the latter changes over a shorter timescale than the
former, i.e. first the BH spin aligns with the total angular momentum, then
the magnitude changes because of accretion at the ISCO. Moreover, Eq. (9) does
not depend on the direction, because it models angular momentum accretion from
the innermost disc region, that is (counter-)aligned with the equatorial plane
of the spinning BH due to the Bardeen-Petterson effect.

We also note that the mass per accretion episode $M_{\rm d}$ might be smaller
or larger than the amount of mass required to be accreted during the
simulation timestep $\Delta t$. Therefore, we follow Bustamante & Springel
(2019) and adopt the following strategy:

  * •

if $M_{d}<\dot{M}_{\rm BH}\Delta t$, multiple accretion episodes occur over
one BH simulation timestep $\Delta t$. We allow for timestep sub-cycles
indexed with a counter variable, therefore executing $N=\dot{M}_{\rm BH}\Delta
t/M_{d}$ accretion episodes. All the sub-cycles share the same accretion rate,
computed at the beginning of the timestep. At the end of each sub-cycle (i.e.
accretion episode) we update BH spin and mass.

  * •

if $M_{d}>\dot{M}_{\rm BH}\Delta t$, we adopt the same strategy extending the
counter across multiple timesteps. The code executes $N=M_{d}/(\dot{M}_{\rm
BH}\Delta t)$ timesteps, then the accretion episode ends and the BH mass and
spin are updated using averages for the accretion rate necessary in Eqs. (4),
(8), and (18), i.e.

$\langle\dot{M}_{\rm
BH}\rangle_{t}=\sum_{i=1}^{N}\dot{M}_{\mathrm{BH},i}\Delta
t_{i}\Big{/}\sum_{i=1}^{N}\Delta t_{i},$ (20)

Fig. 3 illustrates how the final value of $a$ after a single accretion episode
depends on $M_{\rm ratio}$, for a few example values of initial $a$, assuming
no misalignment is present. Note that the plotted lines do not correspond to
evolutionary tracks of the BH in our simulations, the latter ones being
defined instead by a succession of multiple accretion episodes, each one
characterised by different initial directions and different values for $M_{\rm
ratio}$ and $r_{\rm isco}$. The solid and dotted lines correspond to counter-
rotating events. Fig. 3 shows that a single counter-rotating accretion episode
on a maximally spinning BH (solid line) would be able to spin it down to
$a=0$, if the accreted mass per episode were $\simeq 25$% of the BH initial
mass. It would require an accreted mass of $\simeq 2$ times the BH mass to
spin it up to $a=1$, in a direction opposite to the initial one. Similarly, a
counter-rotating accretion episode on a BH with $a=0.5$ (dotted line) would
require $\simeq 0.1M_{\rm BH}$ and $1.75M_{\rm BH}$ to achieve the same
results, respectively. An accretion episode on a non-spinning BH (dashed line)
needs to be $1.5M_{\rm BH}$ to spin it up to its maximal value, whereas it
needs to be equal to $\simeq M_{\rm BH}$ to obtain the same result on a BH
with $a=0.5$ in co-rotating conditions (dash-dotted line).

####  2.2.3 Self-gravity regime

Figure 3: Final BH spin dimensionless parameter after a single accretion
episode as a function of $M_{\rm ratio}$ as defined in Eq. (10). The solid,
dotted, dashed and dot-dashed lines represent $a^{\rm f}$ for the initial spin
values -1, -0.5, 0, and 0.5 respectively. The solid and dotted lines represent
the event of an accretion episode in which the accretion disc is counter-
rotating with respect to the BH.

Depending on the physical conditions, parts of the accretion disc could become
unstable because of their own self-gravity (Dotti et al. 2013). This occurs at
radii that are beyond

$\frac{R_{\mathrm{sg}}}{R_{\mathrm{BH}}}\simeq 5\times
10^{2}M_{\mathrm{BH},8}^{-52/45}\left(\frac{f_{\rm
Edd}}{\epsilon_{\mathrm{r},01}}\right)^{-22/45}\alpha_{\nu_{\rm
1},01}^{28/45}.$ (21)

The mass stable against fragmentation (see Dotti et al. 2013) is then

$M_{\mathrm{sg}}\simeq 6\times
10^{5}M_{\mathrm{BH},8}^{34/45}\left(\frac{f_{\rm
Edd}}{\epsilon_{\mathrm{r},01}}\right)^{4/45}\alpha_{\mathrm{\nu_{1}},01}^{-1/45}\mathrm{M}_{\odot}.$
(22)

Under this condition, only the region of the disc within $R_{\rm sg}$ can be
accreted. Therefore, if $R_{\rm sg}<R_{\rm w}$ we set $M_{\rm d}=M_{\rm sg}$
and substitute $R_{\rm w}$ with $R_{\rm sg}$ in Eqs. (15) and (18). In this
case,

$\displaystyle\frac{J_{\mathrm{d}}}{2J_{\mathrm{BH}}}$ $\displaystyle\simeq
9.4\times 10^{-2}a^{-1}M_{\mathrm{BH},8}^{-37/45}\left(\frac{f_{\rm
Edd}}{\epsilon_{\mathrm{r},01}}\right)^{-7/45}\alpha_{\nu_{\rm
1},01}^{13/45}.$ (23)

Whenever a BH is seeded, its spin is set to zero. As soon as $\dot{M}_{\rm
BH}>0$, a disc is initialised with an angular momentum equal to

$J_{\rm d}=M_{\rm sg}(GM_{\rm BH}R_{\rm sg})^{1/2},$ (24)

since we assume that only the portion of the disc that is stable against
fragmentation can eventually accrete on the BH. The rest of the algorithm
proceeds as before. $M_{\rm sg}$ and $R_{\rm sg}$ are computed using Eqs. (21)
and (22) respectively.

###  2.3 BH mergers

We implement a prescription to account for spin evolution in the case of BH
mergers. In what follows, we adopt the same strategy as Dubois et al. (2014b)
and Bustamante & Springel (2019). We use equations retrieved in a full
general-relativistic framework by Rezzolla et al. (2008a) to compute the final
spin after a merger event. Following their notation, we define
$\bm{a}=a\bm{j}_{\rm BH}$. Once two BHs characterised by masses $M_{\rm 1}$,
$M_{\rm 2}$ and spins $a_{\rm 1}$, $a_{\rm 2}$ have been selected to merge
according to the criteria described in Sec. 2.1, the final spin vector of the
remnant BH $\bm{a}^{f}$ is given by

$\bm{a}^{f}=\frac{1}{(1+q)^{2}}\left(\bm{a}_{1}+\bm{a}_{2}q^{2}+\bm{\ell}q\right),$
(25)

where $q=M_{2}/M_{1}$, with $M_{1}\geq M_{2}$. In this formula, $M_{\rm 1}$
and $M_{\rm 2}$ are the BH physical masses.
$\bm{\ell}=\bm{\ell}^{\prime}/(M_{1}M_{2})$ where $\bm{\ell}^{\prime}$ is the
binary orbital angular momentum that cannot be radiated away in gravitational
waves before coalescence. The magnitude of $\bm{\ell}$ is provided by Rezzolla
et al. (2008a) and reads

$\displaystyle\ell=$
$\displaystyle\frac{s_{4}}{\left(1+q^{2}\right)^{2}}\left(a_{1}^{2}+a_{2}^{2}q^{4}+2\bm{a}_{1}\cdot\bm{a}_{2}q^{2}\right)$
(26)
$\displaystyle+\left(\frac{s_{5}\mu+t_{0}+2}{1+q^{2}}\right)\left(a_{1}\cos\phi_{1}+a_{2}q^{2}\cos\phi_{2}\right)$
$\displaystyle+2\sqrt{3}+t_{2}\mu+t_{3}\mu^{2}.$

Here, $\cos{\phi}=\bm{a}\cdot\bm{\ell}/(a\ell)$ depicts the angle subtended by
the spin vector of BH with $\bm{\ell}$, $\mu=q/(1+q)^{2}$.
$s_{4}=-0.129,s_{5}=-0.384,t_{0}=-2.686,t_{2}=-3.454$, $t_{3}=2.353$ are the
parameters of the fit to their numerical results. We assume that $\bm{\ell}$
is parallel to the binary angular momentum $\bm{L}=\bm{L}_{1}+\bm{L}_{2}$
(Rezzolla et al. 2008a). $\bm{L}_{i=1,2}$ is the angular momentum vector of BH
$i$ with respect to the binary centre of mass (CM), computed as
$\bm{L}_{i}=M_{i}(\bm{r}_{i}-\bm{r}_{\rm CM})\times(\bm{v}_{i}-\bm{v}_{\rm
CM})$. The values for the radii $\bm{r}_{i}$ and velocities $\bm{v}_{i}$ are
retrieved when the BHs are selected to merge.

###  2.4 AGN feedback

Our AGN feedback model assumes that a fraction $\epsilon_{\rm f}$ of the
radiated energy is coupled to the local ISM (refer to Table 2 and following
sections for values), as often adopted in cosmological simulations (see e.g.
Springel et al. 2005a; Booth & Schaye 2009). The total rate of coupled energy
is thus:

$\dot{E}=\epsilon_{\rm f}\epsilon_{\rm r}\dot{M}_{\rm BH}c^{2},$ (27)

where $\epsilon_{\rm f}$ is the coupling efficiency of the feedback energy to
the surrounding ISM. An amount of energy $\Delta E=\dot{E}\Delta t$ is
distributed among the particles surrounding each BH in a weighted fashion,
using the SPH kernel. Moreover, we follow Hirschmann et al. (2014) and
implement a transition from quasar- to maintenance-mode feedback by assuming a
coupling efficiency $\epsilon_{\rm f}$ larger by a factor of 4 when $f_{\rm
Edd}<0.01$. In our model, the radiative efficiency $\epsilon_{\rm r}$ depends
on $a$ through Eq. (14). This is at variance with other implementations where
$\epsilon_{\rm r}$ is either kept fixed (e.g., Springel et al. 2005b; Booth &
Schaye 2009; Hirschmann et al. 2014) or depends on the accretion rate and BH
mass using a phenomenological prescription (e.g., Steinborn et al. 2015).

##  3 The suite of simulations

In the following sections we present our suite of simulations. The simulations
are carried out with an advanced formulation of SPH (Beck et al. 2016),
including a low-viscosity scheme (Dolag et al. 2005; Donnert et al. 2013) and
isotropic thermal conduction (Dolag et al. 2004). The code also features
metallicity-dependent radiative cooling, using the procedure described in
Wiersma et al. (2009), and accounts for the presence of a uniform, time-
dependent ionising background (Haardt & Madau 2001). We adopt the model by
Springel & Hernquist (2003) to include a stochastic treatment of star
formation and to describe the multiphase ISM, so that each gas particle (with
hydrogen number density above a threshold $n_{\rm H}=0.5\;\rm cm^{-3}$) is
made up of a cold, star-forming phase in pressure equilibrium with a hot
phase. A detailed chemical evolution model (Tornatore et al. 2007) allows us
to account for the metal enrichment of the ISM by ageing and exploding stars.
We assume the lifetime function by Padovani & Matteucci (1993), and stellar
yields for supernovae type Ia by Thielemann et al. (2003), supernovae type II
by Nomoto et al. (2013) and asymptotic giant branch stars by Karakas &
Lattanzio (2007). Each star particle represents a simple stellar population
described by a Chabrier (2003) initial mass function.

We consider two idealised cases, a galaxy in isolation (Sec. 3.1) and a galaxy
merger (Sec. 3.2), to test the model in a simple setup and study the evolution
of the BH spin in a well-controlled environment. We also consider three setups
in a full cosmological context (Sec. 3.3 and 3.4), which is key to capture the
effects that the complex interplay between accretion and feedback can have on
BH spin evolution. All the tests are performed including cooling, as well as
star formation and evolution. The cosmological simulations include AGN and
stellar feedback, whereas the idealised tests are performed with stellar
feedback but no AGN feedback.

Table 1: Parameter summary for the idealised tests. $M_{\rm BH,0}$ is the BH
mass in the ICs, $\theta_{z,0}$ is the angle subtended by the initial BH spin
direction and the positive z axis, $f_{\rm Edd}=\dot{M}_{\rm BH}/\dot{M}_{\rm
Edd}$ is the Eddington ratio. Idealised Milky Way galaxy
---
| $\theta_{z,0}\;\rm[^{\circ}]$ | $M_{\rm BH,0}\;\rm[M_{\odot}]$ | $f_{\rm Edd}$
IdealGal-fid | 170 | $5\times 10^{6}$ | 1
IdealGal-120 | 120 | $5\times 10^{6}$ | 1
IdealGal-30 | 30 | $5\times 10^{6}$ | 1
IdealGal-5e5 | 170 | $5\times 10^{5}$ | 1
IdealGal-5e7 | 170 | $5\times 10^{7}$ | 1
IdealGal-0.3 | 170 | $5\times 10^{6}$ | 0.3

###  3.1 Idealised Milky Way galaxy

The first test we present is based upon the setup described in Steinwandel et
al. (2019) aimed at modelling a Milky Way-like galaxy with total halo mass of
$M_{200}=10^{12}\;\rm M_{\rm\odot}$ in isolation with a BH at the potential
minimum. $M_{200}$ is the mass enclosed within the spherical region whose
average density is 200 times the critical density of the Universe. We defer
the reader to their paper for the detailed procedure used to generate the ICs.

Table 1 summarises the initial properties of the central BH in our tests. We
consider a fiducial case (IdealGal-fid), initialised with $M_{\rm BH,0}=5\cdot
10^{6}\leavevmode\nobreak\ {\rm M}_{\sun}$, $a_{\rm 0}=0.5$, $\theta_{\rm
z,0}=170^{\circ}$ and $f_{\rm Edd}=1$. $\theta_{\rm z}$ is the angle subtended
by the BH spin and the z-axis. We perform further tests modifying one
parameter at a time with respect to the fiducial run. Each simulation is
labelled according to the value of the modified parameter, as in Table 1. Each
test assumes a fixed Eddington fraction $f_{\rm Edd}$ for the duration of the
simulation.

The resolution is the same for all the tests. Masses of DM, gas and star
particles are: $m_{\rm DM}=9.6\times 10^{5}\;\rm M_{\rm\odot}$, $m_{\rm
g}=m_{\rm*}=4.8\times 10^{4}\;\rm M_{\rm\odot}$, respectively; the softening
lengths are: $\varepsilon_{\rm DM}=218\;\rm pc$ and $\varepsilon_{\rm
g}=\varepsilon_{\rm*}=\varepsilon_{\rm BH}=50\;\rm pc$.

###  3.2 Idealised galaxy merger

The second set of ICs describes a galaxy merger with a mass ratio of 1:1,
initialised following the procedure described in Karademir et al. (2019).

The setup consists of two identical Milky Way-like galaxies, each with mass
$M_{\rm 200}=1.89\cdot 10^{12}\;\rm M_{\odot}$. Each of them is embedded in a
DM halo and has a spherical stellar bulge component, as well as an exponential
stellar and gas discs. The initial separation of their CMs is 80 kpc and they
rotate one around the other in the same direction in which they complete their
revolution. The angle between the line connecting the two galaxy centres and
their initial velocity vectors is equal to $40^{\circ}$. Fig. 4 shows the
setup. The BHs are initialised with $a=0$ and $M_{\rm BH}=2\times 10^{5}\;\rm
M_{\odot}$.

Figure 4: Gas density map of the ICs for the galaxy merger. The black arrows
trace the local gas projected velocity field, while the white arrows indicate
the initial direction of the CM of each galaxy. The arrows are scaled
arbitrarily for visualisation purposes.

In this setup, $m_{\rm DM}=3.22\times 10^{6}\;\rm M_{\rm\odot}$, $m_{\rm
g}=m_{\rm*}=1.86\times 10^{6}\;\rm M_{\rm\odot}$, $\varepsilon_{\rm
DM}=83\;\rm pc$ and $\varepsilon_{\rm g}=\varepsilon_{\rm*}=\varepsilon_{\rm
BH}=20\;\rm pc$.

Table 2:  Parameter summary of the cosmological simulations. $a_{\rm seed}$ is the BH spin parameter at seeding; $\epsilon_{f}$ is the feedback coupling efficiency (Sec. 2.4); $M_{\rm BH,seed}$ is the BH mass at seeding; $M_{\rm*,seed}$ is the stellar mass considered for BH seeding (Sec. 2.1); $m_{\rm DM}$ is the DM particle mass; $m_{\rm g}$ is the initial gas particle mass; $m_{\rm*}$ is the star particle mass; $\varepsilon_{\rm DM}$, $\varepsilon_{\rm g}$, and $\varepsilon_{\rm*}$ are the DM, gas, and stellar and BH softening length, respectively.  | $a_{\rm seed}$ | $\epsilon_{f}$ | $M_{\rm BH,seed}$ | $M_{\rm*,seed}$ | $m_{\rm DM}$ | $m_{\rm g}$ | $m_{\rm*}$ | $\varepsilon_{\rm DM}$ | $\varepsilon_{\rm g}$ | $\varepsilon_{\rm*,BH}$
---|---|---|---|---|---|---|---|---|---|---
|  |  | $\rm[M_{\odot}/h]$ | $\rm[M_{\odot}/h]$ | $\rm[M_{\odot}/h]$ | $\rm[M_{\odot}/h]$ | $\rm[M_{\odot}/h]$ | $\rm[kpc/h]$ | $\rm[kpc/h]$ | $\rm[kpc/h]$
Zoom-in simulations
asin/dfrogin | 0 | 0.05 | $2\times 10^{5}$ | $1\times 10^{9}$ | $3.3\times 10^{7}$ | $6.25\times 10^{7}$ | $1.6\times 10^{7}$ | 1 | 1 | 0.25
Cosmological box
Box4 | 0 | 0.0775 | $4\times 10^{5}$ | $1.6\times 10^{10}$ | $6.9\times 10^{8}$ | $1.4\times 10^{8}$ | $3.5\times 10^{7}$ | 3.75 | 3.75 | 2.0

###  3.3 Zoom-in simulations

We perform our zoom-in simulations starting from ICs generated with the
procedure described in Bonafede et al. (2011). We consider two halos: a region
with target mass $M_{\rm 200}\simeq 10^{12}\;h^{-1}\;\rm M_{\odot}$, dubbed
asin; a region with target mass $M_{\rm 200}\simeq 10^{13}\;h^{-1}\;\rm
M_{\odot}$, dubbed dfrogin. They are selected from a parent DM-only simulation
of a periodic box with side length $L_{\rm box}=1\;h^{-1}\;\rm{cGpc}$ and
resolution $m_{\rm DM}=10^{9}\;h^{-1}\;\leavevmode\nobreak\ {\rm M}_{\sun}$.
The final resolution is quoted in Table 2, together with the properties and
parameters of each simulation. This set of simulations adopts a flat
$\Lambda$CDM cosmological model, with a matter density parameter $\Omega_{\rm
m}=0.24$, a baryon density parameter $\Omega_{\rm b}=0.04$ and $h=0.72$. The
initial power spectrum follows $n=0.96$ and is normalised to $\sigma_{\rm
8}=0.8$.

###  3.4 Cosmological box

We carry out a simulation starting from the ICs of the cosmological volume
with $L_{\rm box}=48\;h^{-1}\;\rm cMpc$ (hereafter Box4) of the
Magneticum333http://www.magneticum.org simulations suite. Such a simulation,
that contains $\sim 10^{3}$ BHs at redshift $z=0$, allows us to analyse the
properties of a statistically significant BH population, rather than focusing
only on the specific behaviour of a few BHs.

Table 2 summarises the properties and parameters of the simulation. It assumes
a flat $\Lambda$CDM model, with $\Omega_{\rm m}=0.272$, $\Omega_{\rm
b}=0.0456$ and $h=0.704$, $n=0.963$ and $\sigma_{\rm 8}=0.809$.

##  4 Results

###  4.1 Idealised Milky Way galaxy

Figure 5: Spin alignment process in the idealised Milky-Way galaxy. Top
panels: angle subtended by the BH angular momentum and the z-axis $\theta_{\rm
z}$. Bottom panels: BH spin parameter $a$. Quantities are plotted as a
function of time (left panels) and of the ratio between accreted and initial
BH mass $M_{\rm BH}/M_{\rm BH,0}$ (right). Each line represents one simulation
of the set summarised in Table 1. Dotted lines represent the instant after
which condition (6) is no longer satisfied.

Fig. 5 presents the BH spin evolution in the tests listed in Table 1. Top left
panel shows the time evolution of $\theta_{\rm z}$, the angle between
$\bm{j}_{\rm BH}$ and the unit vector indicating the positive z-axis
$\bm{j}_{\rm z}$. The isolated galaxy is initialised with the angular momentum
of the accreting gas along $\bm{j}_{\rm z}$ (i.e. $\theta_{\rm z}=\pi$ refers
to the BH spin direction being anti-parallel to the galaxy angular momentum,
while $\theta_{\rm z}=0$ refers to complete alignment). According to condition
(6), counter-rotating accretion requires both $\bm{j}_{\rm d}\cdot\bm{j}_{\rm
BH}>\pi/2$ and $J_{\rm d}<2J_{\rm BH}$. In IdealGal-fid, IdealGal-120,
IdealGal-5e5, IdealGal-5e7 and IdealGal-0.3 both conditions are satisfied at
the beginning. However, the accreting gas in the BH surroundings keeps the
same direction across the entire simulation, i.e. $\bm{j}_{\rm g}$ remains
constant. As a result, the initial sub-grid disc direction $\bm{j}_{\rm d}$
for every accretion episode is also constant (Eq. 16, $\bm{j}_{\rm
g}=\bm{j}_{\rm d}$). Therefore, each of them induces the BH spin to tilt
towards the angular momentum direction of the large-scale external reservoir.
Note that also $J_{\rm d}/J_{\rm BH}$ is approximately constant in these runs
($a$ and $\epsilon_{\rm r}$ vary slightly, but $J_{\rm d}/J_{\rm BH}$ depends
weakly on them, see Eq. 19). Whether or not condition (6) is verified in these
tests depends only on $\bm{j}_{\rm d}\cdot\bm{j}_{\rm BH}$ (actually only on
$\bm{j}_{\rm BH}$, since $\bm{j}_{\rm d}$ is constant). Each vertical dotted
line marks the first instant in which accretion becomes co-rotating (from
counter-rotating; i.e., condition (6) is no longer verified after this
instant). In IdealGal-fid the BH spin completes its alignment with the
external reservoir in $\sim 2$ Myr. If we keep the same Eddington ratio but
change the BH mass (IdealGal-5e5 and IdealGal-5e7, respectively), alignment
takes the same time. On the other hand, alignment takes longer ($\sim 5$ Myr)
to complete with the same mass but a lower Eddington ratio (IdealGal-0.3). At
fixed BH mass and Eddington ratio (IdealGal-fid, -120, and -30), the alignment
timescale depends weakly on the initial misalignment angle. A smaller
$\theta_{\rm z,0}$ leads to a faster alignment.

The bottom left panel of Fig. 5 illustrates the time evolution of the BH spin
parameter $a$. Since IdealGal-fid, IdealGal-120, IdealGal-5e5, IdealGal-5e7
and IdealGal-0.3 start with counter-rotating accretion conditions, $a$
decreases. At $t\simeq 1$ Myr co-rotating conditions take over and $a$ starts
increasing. The turnaround point occurs in correspondence of the dotted line,
after which condition (6) is no longer verified. In IdealGal-30, $a$
monotonically increases as condition (6) is not verified from the very
beginning.

The top right panel of Fig. 5 shows $\theta_{\rm z}$ as a function of the
ratio $M_{\rm BH}/M_{\rm BH,0}$. At fixed $\theta_{\rm z,0}$ (IdealGal-fid,
-5e5, -5e7, and -0.3), the same ratio of accreted mass over BH mass ($M_{\rm
BH}\sim 5\%\,M_{\rm BH,0}$) is needed for complete alignment to occur,
regardless of other parameters. The actual time it takes to align depends on
how fast the BH accretion proceeds. BHs whose spin starts with smaller
$\theta_{\rm z,0}$ (IdealGal-120 and -30) require less time.

The bottom right panel of Fig. 5 presents $a$ as a function of the ratio
$M_{\rm BH}/M_{\rm BH,0}$. Similarly to the bottom left panel, the turnaround
point occurs in correspondence of the dotted line, when accretion switches
from counter- to co-rotating. This plot further shows that it occurs at about
the same value of $M_{\rm BH}/M_{\rm BH,0}$ at fixed $\theta_{\rm z,0}$.
Finally, for smaller $\theta_{\rm z,0}$ a smaller ratio is required to reach
the turnaround point.

The behaviour of both $\theta_{\rm z}$ and $a$ as a function of $M_{\rm
BH}/M_{\rm BH,0}$ is similar for simulations assuming the same $\theta_{\rm
z,0}$ because $f_{\rm Edd}$ is fixed. Indeed, the BH grows according to the
following differential equation

$\frac{\mathrm{d}M_{\rm BH}}{\mathrm{d}t}=(1-\epsilon_{\rm r})f_{\rm
Edd}\dot{M}_{\rm Edd}=f_{\rm Edd}\frac{1-\epsilon_{\rm r}}{\epsilon_{\rm
r}}\frac{M_{\rm BH}}{\tau_{\rm S}},$ (28)

where $\tau_{\rm S}=\sigma_{\rm T}c/(4\pi Gm_{\rm p})\sim 4.5\times 10^{8}$ yr
is the Salpeter time. As a result,

$\frac{M_{\rm BH}}{M_{\rm BH,0}}=\exp{\left(f_{\rm Edd}\frac{1-\epsilon_{\rm
r}}{\epsilon_{\rm r}}\frac{t}{\tau_{\rm S}}\right)}.$ (29)

Therefore, the evolution of $\theta_{\rm z}$ and $a$ as a function of time is
different depending on $f_{\rm Edd}$ (compare IdealGal-fid and IdealGal-0.3).
On the other hand, the curves overlap if we consider $M_{\rm BH}/M_{\rm BH,0}$
as the independent variable (and they are independent of $M_{\rm BH,0}$ – see
IdealGal-fid, IdealGal-5e5, IdealGal-5e7).

###  4.2 Idealised galaxy merger

Figure 6: Time evolution of a few properties of the two BHs (each of them is
depicted by either blue solid or orange dotted lines) in the idealised galaxy
merger simulation. From top to bottom: separation distance between the two
BHs; BH dimensionless spin parameter; BH mass; x, y, z component of the unit
vector of the BH angular momentum. Left panel: entire simulated timespan.
Right panel: timespan of 2 Myr, centered on the time of the merger (gray,
dotted vertical line). The vertical green dash-dotted and purple dashed lines
mark the instants in time illustrated in the top and bottom panel of Fig. 7,
respectively.

Fig. 6 presents the evolution of the spins of the two BHs (one identified with
the solid blue lines and the other with the dashed orange lines) in the
idealised merger simulation. The left panel shows the evolution of quantities
for the entire simulated timespan, whereas the right panel shows a 2 Myr
window, centred on the merger instant, marked in both panels by the vertical
grey dotted line.

The first row from the top of Fig. 6 shows the separation distance between the
BHs. The galaxies undergo 3 pericenter passages before the BHs merge (left
panel), at $\sim 336.15$ Myr (right panel).

The second row of Fig. 6 follows the evolution of the spin parameter $a$. The
left panel shows that both BHs follow an identical path from zero spin to
maximal within 50 Myr (the two galaxies of the pair are identical). They are
kept at maximal spin by accretion occurring consistently on the same plane,
until the third pericenter passage. At this stage accretion becomes less
coherent – i.e. the direction of the local angular momentum has larger
variability – and $a$ decreases slightly due to counter-rotating accretion
conditions. In the right panel, second row, just before the merger instant we
identify two progenitors with $a\sim 0.9$. The final remnant has $a\sim 0.9$.

The third row of Fig. 6 shows the evolution of $M_{\rm BH}$. On the left, we
see that the mass increases steadily and identically for the two BHs. The
right panel allows us to identify two equal-mass progenitors with $M_{\rm
BH}=1.2\times 10^{6}\leavevmode\nobreak\ {\rm M}_{\sun}$ merging into a
remnant of $M_{\rm BH}=2.4\times 10^{6}\leavevmode\nobreak\ {\rm M}_{\sun}$.

Rows 4 to 6 of Fig. 6 show the three components of the BH spin direction. In
the left panel we see that the two BHs consistently accrete from the same
plane, hence their spins are kept very well aligned with the angular momentum
direction of the host galaxies, which is along $\bm{j}_{\rm z}$. The
pericenter passages create some disturbance in the gas close to the BH,
temporarily inducing the BH spin to tilt. This effect is more prominent close
to the merger. In the right panel, we see that the BH spin components
immediately before the merger are close to $\bm{j}_{\rm z}$, although with
some disturbances induced by the local environment. Furthermore, the spin of
the remnant is also aligned with $\bm{j}_{\rm z}$.

Fig. 7 represents visually the process just described, using a volume
rendering of the gas density in the simulation, for two simulation snapshots.
The field of view depicted in the figure spans a width of 55 kpc and a height
of 31 kpc, as measured along a plane that intersects the centre of the
rendering volume, offering a perspective view of the merging galaxies. The
arrows represent the instantaneous directions of the BH spins. The top panel
of Fig. 7 corresponds to the instant marked by the green dash-dotted line in
the right panel of Fig. 6, showing the pre-merger configuration. The bottom
panel of Fig. 7 corresponds to the time marked by the purple dashed line,
displaying the post-merger configuration.

Figure 7: 3D perspective volumetric rendering of the gas density in the
idealised galaxy merger simulation. The field of view depicted in the figure
spans a width of 55 kpc and a height of 31 kpc, as measured along a plane that
intersects the centre of the rendering volume. Top panel: last snapshot before
the merger. Bottom panel: first snapshot after merger. The arrows mark the
instantaneous direction of the BH spin vector. (Movie online)

The idealised galaxy merger test enables a simplified computation of the
expected BH spin of the merger remnant to validate the result obtained. The
pre-merger configuration (as shown in the right panels of Fig. 6 and in the
top panel of Fig. 7) can be approximated to the idealised configuration of an
equal-mass BH binary with aligned spins treated in Rezzolla et al. (2008b). In
this configuration, the final spin magnitude $a_{\rm f}$ depends only on the
initial spin magnitudes $a_{\rm 1}$ and $a_{\rm 2}$ (and the final direction
is equal to the initial one). Two BHs with $a\sim 0.9$ and spins in the same
direction are indeed expected to end up in a merger remnant with $a\sim 0.9$
and spin along the original, common rotation axis, as observed in our test.

###  4.3 Zoom-in simulations

Figure 8: Gas surface density maps at $z=0$ for the asin (top) and dfrogin
(bottom) runs. The white circle indicates a spherical region of radius
$5R_{\rm 200}$, centred on the target halo of the zoom-in region (i.e. the
most massive at $z=0$). The white dots mark the positions of the BHs in the
simulation.

Fig. 8 presents projected gas density maps at $z=0$ of asin (top) and dfrogin
(bottom). The white dots mark the positions of the BHs in the simulation. In
what follows, we consider only the BHs that lie within a spherical region of
radius $5R_{\rm 200}$ at $z=0$, centred on the target halo of interest of the
zoom-in region (i.e. the most massive one at $z=0$). $R_{\rm 200}$ is the
radius of the spherical volume, centred on the subhalo, whose average density
is 200 times the critical density of the Universe. The region is marked by the
white circle in Fig. 8. The main purpose of this selection is to exclude BHs
that are close the border of the high-resolution region.

Figure 9: $M_{\rm BH}$ as a function of stellar mass $M_{\rm*}$, at $z=0$.
Diamonds and stars represent BHs in asin and dfrogin, respectively. The dashed
line shows the experimental fit by McConnell & Ma (2013), while crosses with
associated uncertainties show data from Kormendy & Ho (2013).

Fig. 9 shows the BH sample at $z=0$ in asin (diamonds) and dfrogin (stars) in
the BH mass-stellar mass ($M_{\rm BH}-M_{\rm*}$) plane. The BHs are selected
according to the criterion explained above; the sample is composed by 4 BHs in
asin with $10^{6.5}\lesssim M_{\rm BH}/\leavevmode\nobreak\ {\rm
M}_{\sun}\lesssim 10^{8}$ and 21 BHs in dfrogin, with $10^{6}\lesssim M_{\rm
BH}/\leavevmode\nobreak\ {\rm M}_{\sun}\lesssim 10^{9}$. Each BH is associated
to a subhalo identified with the sub-structure finder algorithm SubFind
(Springel et al. 2001; Dolag et al. 2009) based on particle ID matching.
Whenever more than one BH is associated to the same subhalo, the closest to
the subhalo center is chosen. $M_{\rm*}$ is the stellar mass as computed by
SubFind. The dashed line shows the experimental fit by McConnell & Ma (2013),
while crosses with associated uncertainties show observations from Kormendy &
Ho (2013).

Figure 10: BH spin parameters $a$ of the selected BH sample in asin (diamonds)
and dfrogin (stars) as a function of $M_{\rm BH}$, at $z=0$. The squares and
pentagons show the collection of observational measurements of the BH spin
parameter by Reynolds (2021) and Mallick et al. (2022), respectively.

Fig. 10 presents the spin parameters of the BHs in the simulated regions as a
function of $M_{\rm BH}$ at $z=0$, with the same symbols used in Fig. 9. We
compare our simulation output to the available observations that provide
estimates of $a$ and BH mass (Reynolds 2021; Mallick et al. 2022), with
associated uncertainties. We observe a range of BH masses ($10^{6}\lesssim
M_{\rm BH}/\leavevmode\nobreak\ {\rm M}_{\sun}\lesssim 5\times 10^{7}$) where
BH spins tend to be larger than $a\sim 0.70$, in both the simulated regions.
Several of them are maximally spinning. For masses above $5\times
10^{7}\leavevmode\nobreak\ {\rm M}_{\sun}$, there are no maximally spinning
BHs, whereas the most massive BH in each region systematically has a lower
spin, $a\simeq 0.25$ in asin and $a\simeq 0.1$ in dfrogin.

Figure 11:  Summary of properties related to the most massive BH in the
dfrogin run, as they evolve across time. From top to bottom: z-component of
the BH spin direction $j_{\rm BH,z}$ (black, solid line) and direction of the
angular momentum of the gas in the BH kernel $j_{\rm g,z}$ (orange, dashed
line); cosine of the angle between the accretion disc and the BH spin at the
beginning of each accretion episode $\cos{\theta_{\rm BH-d}}$ (i.e. left-hand
side of Eq. (6), black solid line) and $-J_{\rm d}/(2J_{\rm BH})$ (i.e. right-
hand side of Eq. (6), orange dashed line); ratio of the magnitudes of the disc
and BH angular momenta $J_{\rm d}/J_{\rm BH}$; BH dimensionless spin parameter
$a$; radiative efficiency of the accretion disc $\epsilon_{\rm r}$; Eddington
ratio $f_{\rm Edd}$; BH mass $M_{\rm BH}$.

Fig. 11 shows the time evolution of a few BH properties predicted by the sub-
resolution model. We restrict our analysis to the most massive BH of the
sample in dfrogin, and focus on the relationship between BH spin and gas
accretion as they evolve with time. The evolutionary tracks of the other BHs
in the sample display features which are similar to those of our reference,
most massive BH as for the interplay between the BH spin and the fuelling gas.
The upper axis of Fig. 11 marks the time since the Big Bang, while the lower
axis shows the redshift $z$. Across the panels, the vertical dashed lines mark
the occurrence of mergers.

The top panel of Fig. 11 shows the z-components of $\bm{j}_{\rm BH}$ (black,
solid line) and $\bm{j}_{\rm g}$ (orange, dashed line). Before 1.7 Gyr,
$\bm{j}_{\rm g}$ varies gradually with time, with a clear, average trend and
limited scatter. From $\simeq 1.7$ to $\simeq 2.6$ Gyr a trend in $\bm{j}_{\rm
g}$ is still present but the scatter increases. $\bm{j}_{\rm BH}$ follows the
average evolution of $\bm{j}_{\rm g}$. From $\simeq 2.6$ to $\simeq 6$ Gyr,
$\bm{j}_{\rm g}$ exhibits drastic and sudden changes, whereas $\bm{j}_{\rm
BH}$ remains stable with minimal variations. After $\simeq 6$ Gyr,
$\bm{j}_{\rm BH}$ often undergoes abrupt variations due to the erratic
behaviour of $\bm{j}_{\rm g}$.

The second row of Fig. 11 shows $\cos\theta_{\rm BH-d}$ (black solid line),
i.e. the left-hand-side term of Eq. (6), as a function of time. The orange
dashed line shows $-J_{\rm d}/(2J_{\rm BH})$, i.e. the right-hand side of Eq.
(6). Before $\simeq 1.7$ Gyr each accretion episode is characterised by
consistently small misalignment. Between $\simeq 1.7$ and $\simeq 2.6$ Gyr
several accretion episodes display $\theta_{\rm BH-d}$ close to $\pi/2$. After
$\simeq 2.6$ Gyr accretion episodes are mostly misaligned. This panel also
allows us to easily identify counter-rotating accretion episodes, when
$\cos{(\theta_{\rm BH-d})}<-J_{\rm d}/(2J_{\rm BH})$ (e.g. around 2.2 Gyr and
4.5 Gyr).

The third row of Fig. 11 shows $J_{\rm d}/J_{\rm BH}$, the ratio between the
accretion disc and the BH angular momentum per accretion episode. This
quantity controls how much a single accretion episode is able to modify the
direction of the BH spin and ultimately determines to which degree the BH spin
is able to follow the variations of $\bm{j}_{\rm g}$. Being the BH spin
direction $\bm{j}_{\rm BH}^{f}$ after an accretion episode parallel to
$\bm{J}_{\rm tot}=\bm{J}_{\rm d}+\bm{J}_{\rm BH}$, if $J_{\rm d}\ll J_{\rm
BH}$, then $\bm{j}_{\rm BH}^{f}\sim\bm{j}_{\rm tot}\sim\bm{j}_{\rm BH}$ (i.e.
the direction change is negligible; Sec. 2.2.2). Conversely, if $J_{\rm d}\gg
J_{\rm BH}$, then $\bm{j}_{\rm BH}^{f}\sim\bm{j}_{\rm d}$, i.e. $\bm{j}_{\rm
BH}$ aligns with the direction imparted by $\bm{j}_{\rm d}$ (and hence
$\bm{j}_{\rm g}$). In an intermediate configuration $\bm{j}_{\rm BH}$ only
partially aligns with $\bm{j}_{\rm d}$. The latter situation is expected to
occur before $\simeq 6$ Gyr, when $J_{\rm d}/J_{\rm BH}\lesssim 0.1$. However,
at $t\gtrsim 6$ Gyr $J_{\rm d}/J_{\rm BH}$ sometimes approaches 1; as a
result, $\bm{j}_{\rm BH}$ exhibits larger variations as a response to large
changes in $\bm{j}_{\rm g}$. Although here we inspect the evolution of $J_{\rm
d}/J_{\rm BH}$ and ${\theta_{\rm BH-d}}$ for one BH as an example, we properly
quantify how the change $\bm{j}_{\rm BH}$ depends on these two quantities by
analysing accretion episodes statistically in the cosmological box (Sec. 4.4).

The fourth row of Fig. 11 shows the BH dimensionless spin parameter $a$. Its
evolution is characterised by a maximally-spinning phase until $t\simeq 1.95$
Gyr ($z\gtrsim 3.3$). The largest spin-downs occur because of mergers,
although we also observe counter-rotating phases that decrease the spin due to
gas accretion (e.g. at $t\simeq 2.2$ and $\simeq 2.35$ Gyr, as well as around
$4.5$ Gyr).

The fifth row of Fig. 11 plots the radiative efficiency, which depends on $a$
(see Eq. 14). The dash-dotted, dashed and dotted lines mark the values of
efficiency for $a=0.998$, 0 and for a counter-rotating episode on a BH with
$a=1$, respectively. The maximally-spinning period before $z\sim 3.5$
corresponds to a maximal efficiency of 0.32, whereas later times are
characterised by a lower efficiency. Counter-rotating accretion conditions are
clearly visible as a drastic decrease in efficiency, close to the minimum
theoretical value marked by the dotted line.

The sixth row of Fig. 11 illustrates the Eddington ratio $f_{\rm Edd}$.
Accretion is Eddington-limited soon after the BH is seeded ($z\sim 4.9$).
$f_{\rm Edd}$ is typically $\gtrsim 10^{-1}$ until $z\sim 2.6$. Between
$0.9\lesssim z\lesssim 2.6$, $f_{\rm Edd}\ll 10^{-1}$ for most of the time. A
few highly accreting ($f_{\rm Edd}\simeq 10^{-1}$) episodes occur in proximity
of the mergers at $z\sim 0.9$ and $z\sim 0.7$, and lead to significant
reorientation of the BH spin (see first row of Fig. 11). After $z\sim 0.7$
$f_{\rm Edd}$ is mostly $\lesssim 10^{-3}$. Comparing with the evolution of
$a$ (fourth row of Fig. 11), we observe that the largest $f_{\rm Edd}$, the
largest is the change induced in $a$.

The last row of Fig. 11 plots $M_{\rm BH}$. The BH increases its mass by three
orders of magnitude during the highly-accreting phase at $z\gtrsim 2.6$. After
$z\sim 2.6$, the BH gets more massive (by about one order of magnitude) mostly
due to mergers.

###  4.4 Cosmological box

For the analysis of Box4 we consider all the BHs in the cosmological volume
and we select them at $z=0$. We associate each BH to a subhalo based on
particle ID matching using SubFind and whenever more than one BH is associated
to the same subhalo, the closest to the subhalo center is chosen. The sample
includes 1790 BHs, in a mass range between $6\times 10^{6}$ and $2\times
10^{10}\leavevmode\nobreak\ {\rm M}_{\sun}$.

Figure 12: $M_{\rm BH}$ as a function of stellar mass $M_{\rm*}$ for the Box4
run. Circles represent the simulated sample at $z=0$, while the observational
points are as in Fig. 9. The points are colour-coded according to the number
of mergers BH have undergone.

Fig. 12 frames the Box4 BH sample at $z\sim 0$ in the $M_{\rm BH}-M_{\rm*}$
plane. $M_{\rm*}$ corresponds to the subhalo stellar mass as computed by
SubFind. We compare our sample to observations by Reynolds (2021) and Mallick
et al. (2022) (as in Fig. 9). Each point is colour-coded according to the
number of mergers of the BH (including their progenitors). The number
increases with increasing BH mass.

Figure 13: Distribution of BH spin parameters as a function of the BH mass in
Box4 (circles), at $z=0$. The squares and pentagons show the collection of
observational measurements of the BH spin parameter by Reynolds (2021) and
Mallick et al. (2022), respectively.

In Fig. 13 we show how BH spin parameters change as a function of $M_{\rm BH}$
in Box4 (similar to Fig. 10). The simulated BHs are indicated by circles. We
compare our simulation output to the observational data by Reynolds (2021) and
Mallick et al. (2022). The BH sample is much more numerous than in the zoom-in
regions and the mass range extends further on the high-end, up to $2\sim
10^{10}\leavevmode\nobreak\ {\rm M}_{\sun}$. We identify three mass ranges in
which we observe different distributions of $a$. Close to the seeding mass
($\sim 5.5\times 10^{5}\leavevmode\nobreak\ {\rm M}_{\sun}$), we can see a
steep increase of $a$ with $M_{\rm BH}$. The BH population characterised by
$10^{6}\lesssim M_{\rm BH}/\leavevmode\nobreak\ {\rm M}_{\sun}\lesssim 2\times
10^{7}$ shows a systematic tendency for highly spinning BHs ($a\gtrsim 0.85$),
with most of them close to the maximal value. BHs with masses above $2\times
10^{7}\leavevmode\nobreak\ {\rm M}_{\sun}$ display a wider range of $a$,
extending as low as $a\sim 0.1$. The robustness of this result is consolidated
by the larger number of BHs probing this high-mass regime with respect to the
zoom-in simulations. We also note that for $M_{\rm BH}\gtrsim 5\times
10^{8}\leavevmode\nobreak\ {\rm M}_{\sun}$ there are no maximally spinning
BHs.

Figure 14: BH dimensionless spin parameter $a$ as a function of $M_{\rm BH}$,
for the five most massive BHs at $z=0$, for the Box4 run. Each line
corresponds to one BH and is colour-coded by redshift. The symbols mark the
position of each BH in the plot at a few specific instants in time,
corresponding to $z={0,1,2,3,4}$, colour-coded accordingly.

Fig. 14 shows the evolution of $a$ of the five most massive BHs in the sample
as a function of mass. The redshift is encoded in the colour gradient of each
curve. The symbols mark the position of each BH in the plot at a few specific
instants in time, corresponding to $z={0,1,2,3,4}$, colour-coded accordingly.
As soon as the BHs are seeded, they undergo a phase of rapid increase of $a$,
then they reach and maintain a maximally spinning state. We also notice that
when $M_{\rm BH}$ is between $10^{6}$ and $2\times 10^{7}\leavevmode\nobreak\
{\rm M}_{\sun}$, BHs undergo short transitory periods of spin-down due to
counter-rotating accretion or mergers. On the other hand, $a$ reaches again
the maximal state afterwards, highlighting that co-rotating accretion is
dominant. After the BHs overcome the $2\times 10^{7}\leavevmode\nobreak\ {\rm
M}_{\sun}$ mass scale, $a$ exhibits a tendency to decrease. In addition, we
observe large and sudden changes due to mergers (when also the mass increases
significantly at the same time). A binary with BH spins oriented towards
opposite directions results in a severe spin-down of the remnant compared to
the state immediately before the merger. The wider distribution of $a$ at the
highest masses in Fig. 13 suggests that spin-down due to counter-rotating
accretion and mergers occurs frequently.

Figure 15: Statistical analysis of a few key properties of the accretion
episodes occurred in the Box4 run. The entire set of accretion episodes
occurred during the simulation has been considered (i.e. for every BH and at
every redshift). The top (middle) panel shows a 2D histogram of the values of
$J_{\rm d}/J_{\rm BH}$ ($\theta_{\rm BH-d}$) as a function of mass. In the top
panel, the dashed line pinpoints $J_{\rm d}/J_{\rm BH}\propto M_{\rm
BH}^{-37/45}$ (for the self-gravitating case, Eq. 23); the dotted line shows
$J_{\rm d}/J_{\rm BH}\propto M_{\rm BH}^{23/16}$ (non self-gravitating case,
Eq. 19). $a=0.998$ and $f_{\rm Edd}=1$ are assumed to plot these reference
lines. The bottom panel shows the fraction of counter-rotating accretion
episodes over the total, per BH mass bin.

In Fig. 15 we carry out a statistical analysis of a few key properties of the
accretion episodes occurred in the Box4 run, to gain insight on the mechanisms
with which gas accretion drives the trends observed in Fig. 13 and 14. For the
analysis, the entire set of accretion episodes occurred during the simulation
is considered (i.e. for every BH at every redshift). The top panel of Fig. 15
shows a 2D histogram where accretion episodes are binned according to their
values of $J_{\rm d}/J_{\rm BH}$ and $M_{\rm BH}$. Each 2D bin is colour-coded
by the number of accretion episodes in that bin. The dashed and dotted lines
represent the dependence of $J_{\rm d}/J_{\rm BH}$ on mass from Eq. 18 (i.e.
$\propto M_{\rm BH}^{-37/45}$ for the self-gravitating case and $\propto
M_{\rm BH}^{23/16}$ for the standard case), assuming $a=0.998$ and $f_{\rm
Edd}=1$. The distribution observed in $J_{\rm d}/J_{\rm BH}$ at fixed mass bin
is due to different values of $a$ and $f_{\rm Edd}$. However, $J_{\rm
d}/J_{\rm BH}$ depends weakly on the latter two quantities, while it depends
on $M_{\rm BH}$ quite strongly (see Eq. 19). Above $M_{\rm BH}\sim
10^{8}\leavevmode\nobreak\ {\rm M}_{\sun}$ accretion occurs mostly in the
self-gravitating regime (i.e. $R_{\rm sg}<R_{\rm w}$). The middle panel of
Fig. 15 shows a 2D histogram where accretion episodes are binned according to
their values of $\theta_{\rm BH-d}$ and $M_{\rm BH}$. The solid white line
indicates the median value of $\theta_{\rm BH-d}$ per BH mass bin, whereas the
shaded region represents the 25th and 75th percentile of the distribution of
$\theta_{\rm BH-d}$ in that BH mass bin. We observe that accretion episodes
are characterised predominantly by small misalignment (75% of them has
$\theta_{\rm BH-d}\lesssim 0.25\text{ rad}\sim 15^{\circ}$) below $M_{\rm
BH}\sim 10^{8}\leavevmode\nobreak\ {\rm M}_{\sun}$. Above this mass threshold,
the distribution broadens with increasing mass, showing that large
misalignment is increasingly more probable. However, we caution that in this
regime less accretion episodes occur, therefore the significance of this trend
is limited by low-number statistics. Note that when $J_{\rm d}\ll J_{\rm BH}$,
the minimum angle required for an episode to satisfy condition (6) is
$\theta_{\rm BH-d}=\pi/2$. The minimum angle is larger for larger $J_{\rm
d}/J_{\rm BH}$, whereas for $J_{\rm d}\geq 2J_{\rm BH}$ accretion will be
always co-rotating regardless the initial misalignment. Since a large
misalignment is more probable at high BH masses (middle panel of Fig. 15)
whereas $J_{\rm d}/J_{\rm BH}$ decreases with mass above $M_{\rm BH}\sim
10^{8}\leavevmode\nobreak\ {\rm M}_{\sun}$ (top panel of Fig. 15), counter-
rotating accretion episodes are more likely. In the bottom panel of Fig. 15 we
plot the fraction of accretion episodes that are counter-rotating with respect
to the total number per BH mass bin. This fraction increases with increasing
mass and it is as high as 0.5 for the highest mass bins.

It is now possible to assess the contribution of gas accretion to the trends
discussed in Fig. 13 and 14 in light of the results shown in Fig. 15. Below
$M_{\rm BH}\sim 10^{8}\leavevmode\nobreak\ {\rm M}_{\sun}$, most of the
accretion episodes are co-rotating, therefore spin-up is favoured. Counter-
rotating accretion episodes do occur, but the decrease in spin is transitory
and co-rotating accretion leads the spin back to maximal. Conversely, above
$M_{\rm BH}\sim 10^{8}\leavevmode\nobreak\ {\rm M}_{\sun}$ the probability of
counter-rotating accretion (and hence spin-down) increases as a function of
mass.

Figure 16: BH spin direction variation per accretion episode (i.e.
$\Delta\theta_{\rm BH}=\bm{j}_{\rm BH}\cdot\bm{j}_{\rm BH}^{f}$) as a function
of $J_{\rm d}/J_{\rm BH}$, colour-coded by $\theta_{\rm BH-d}$, the
misalignment angle between disc and BH angular momenta at the beginning of the
episode. The dashed lines illustrate the analytical dependence of $\bm{j}_{\rm
BH}\cdot\bm{j}_{\rm BH}^{f}$ on $J_{\rm d}/J_{\rm BH}$, computed using Eq. 5.

In Fig. 16 we quantify the direction variation imparted to the BH spin as a
function of the properties of each accretion episode. The plot shows the angle
$\Delta\theta_{\rm BH}$ between the direction of the BH spin before and after
each accretion episode (i.e. $\bm{j}_{\rm BH}\cdot\bm{j}_{\rm BH}^{f}$) in the
Box4 run, as a function of $J_{\rm d}/J_{\rm BH}$. The accretion episodes
(circles) are colour-coded by $\theta_{\rm BH-d}$, i.e. the BH spin-disc
misalignment at the beginning of the episode. The dashed lines indicate the
analytical dependence computed using Eq. (5) and expressing $\bm{j}_{\rm
BH}\cdot\bm{j}_{\rm BH}^{f}$ as a function of $J_{\rm d}/J_{\rm BH}$, at fixed
$\theta_{\rm BH-d}$. If $J_{\rm d}\ll J_{\rm BH}$ then $\Delta\theta_{\rm
BH}\sim 0$, regardless $\theta_{\rm BH-d}$. If $J_{\rm d}\gg J_{\rm BH}$, then
$\Delta\theta_{\rm BH}\sim\theta_{\rm BH-d}$. Increasing values of $J_{\rm
d}/J_{\rm BH}$ induce larger alignment of $\bm{j}_{\rm BH}$ with $\bm{j}_{\rm
d}$ per accretion episode. Fig. 16 also shows that, overall, $J_{\rm d}/J_{\rm
BH}\gtrsim 1$ is rare. Therefore complete alignment of $\bm{j}_{\rm BH}$ with
$\bm{j}_{\rm d}$ ($\Delta\theta_{\rm BH}=\theta_{\rm BH-d}$) never occurs in a
single accretion episode, i.e. accretion episodes lead at most to partial
alignment with the instantaneous direction of the accreting gas.

Figure 17: Radiative efficiencies of the BH populations at redshift $z=0$, as
a function of mass, for the Box4 run. The triangles show the collection of
empirical estimates of the radiative efficiency by Daly (2021). The dotted,
dashed and dash-dotted lines mark the values of the efficiency corresponding
to $a=-1,0,0.998$, respectively, for reference. The population in the bottom
part of the figure represents the BH that are accreting in counter-rotating
conditions.

Fig. 17 shows the radiative efficiency $\epsilon_{\rm r}$ across the BH sample
at $z=0$, as a function of $M_{\rm BH}$. The trends shown in Fig. 13 are
reflected in the distribution of $\epsilon_{\rm r}$. Indeed, each BH has its
own value of $\epsilon_{\rm r}$ at each instant, dependent on $a$ (see Fig.
1). We observe predominantly high efficiency ($\sim 0.32$) at intermediate BH
masses ($10^{6}\lesssim M_{\rm BH}/\leavevmode\nobreak\ {\rm M}_{\sun}\lesssim
2\times 10^{7}$). A lower efficiency value becomes more likely at higher
masses (above $4\times 10^{7}\leavevmode\nobreak\ {\rm M}_{\sun}$), whereas at
the highest masses (above $5\times 10^{8}\leavevmode\nobreak\ {\rm M}_{\sun}$)
efficiencies have systematically lower values ($\sim 0.06-0.1$). The radiative
efficiency also depends on whether, at a given instant, accretion is
proceeding in co- or counter-rotating accretion conditions. Points with
$\epsilon_{\rm r}$ below the dashed line correspond to BHs that are accreting
in counter-rotating conditions. Indeed, the probability of having counter-
rotating conditions increases with mass above $4\times
10^{7}\leavevmode\nobreak\ {\rm M}_{\sun}$ (bottom panel of Fig. 15). We also
compare our simulated sample with a collection of radiative efficiency factors
provided in Daly (2021). Since the simulation points refer to the sample at
$z=0$, we consider all the sources in the observational catalogue that have
$z<0.2$444Note that we exclude the sources catalogued as LINERs..

##  5 Discussion

###  5.1 Evolution of the BH spin direction

Our algorithm for spin evolution proceeds through accretion episodes that
modify $\bm{j}_{\rm BH}$ as an effect of the external change in $\bm{j}_{\rm
g}$. As shown in Fig. 16, only if $J_{\rm d}\gg J_{\rm BH}$ a single accretion
episode is able to induce complete alignment with $\bm{j}_{\rm d}$ and
therefore with $\bm{j}_{\rm g}$. Conversely, if $J_{\rm d}\ll J_{\rm BH}$,
then $J_{\rm BH}$ is insensitive to external change. We observe that in our
simulations accretion episodes are in general characterised by $J_{\rm
d}/J_{\rm BH}\lesssim 1$ (Fig. 16), therefore accretion episodes induce only
partial alignment of $\bm{j}_{\rm BH}$ with $\bm{j}_{\rm d}$. Larger values of
$J_{\rm d}/J_{\rm BH}$ reduce the misalignment between $\bm{j}_{\rm BH}$ and
$\bm{j}_{\rm g}$ to a larger degree. The relation between $\bm{j}_{\rm BH}$
and $\bm{j}_{\rm d}$ (and hence $\bm{j}_{\rm g}$) on timescales longer than a
single accretion episode depends on $J_{\rm d}/J_{\rm BH}$ and on the
variability of $\bm{j}_{\rm g}$. For the reference BH considered in dfrogin,
the top panel of Fig. 11 shows that before $t\sim 2.6$ Gyr $\bm{j}_{\rm g}$
varies gradually with time, although with some variability on short (i.e.
$\lesssim$ 1 Myr) timescales. On the other hand, $\bm{j}_{\rm g}$ changes
erratically after $t\sim 2.6$ Gyr. In the former case $\bm{j}_{\rm BH}$
manages to follow the average evolution of $\bm{j}_{\rm g}$, whereas the two
directions are decoupled in the latter. The middle panel of Fig. 15 shows that
statistically such large misalignment is more probable at the highest masses.
Combined with low values of $J_{\rm d}/J_{\rm BH}$, it results in counter-
rotating accretion conditions to be more frequent with increasing mass (bottom
panel of Fig. 15).

###  5.2 Evolution of the BH spin magnitude

The BH spin magnitude evolution via gas accretion is driven by the amount of
accreted mass and by the radius of the ISCO (Eq. 9). Therefore, it evolves
more rapidly in high accretion rate phases. Furthermore, counter-rotating
accreting gas is characterised by a larger ISCO (bottom panel of Fig. 1), thus
a larger angular momentum per unit mass. The same amount of accreted mass
induces a larger (negative) change in $a$ than if it were accreted in co-
rotating conditions. At fixed mass, if co- and counter-rotating episodes
occurred in equal number, the net effect would actually be a decrease in $a$
(Dotti et al. 2013). Whether there is a trend to increase or decrease the spin
magnitude depends on the accretion rate and on how frequent co- or counter-
rotating accretion is. Moreover, mergers also contributes to influence $a$. In
Figs. 10 and 13 we analyse the trends of $a$ with BH mass and observe that BHs
with $M_{\rm BH}\lesssim 2\times 10^{7}\leavevmode\nobreak\ {\rm M}_{\sun}$
show systematic spin-up. This effect is attributed to the prevalence of co-
rotating accretion conditions (bottom panel of Fig. 15). Mergers do not occur
in this range of masses (see Fig. 12), hence the trend is exclusively due to
gas accretion. For $M_{\rm BH}\gtrsim 2\times 10^{7}\leavevmode\nobreak\ {\rm
M}_{\sun}$ we observe a wide distribution of $a$, indicating that the
behaviour strongly depends on the detailed history, and there is no systematic
behaviour. In this mass regime, several elements contribute to the trend: $i)$
the fraction of counter-rotating episodes increases with mass (bottom panel of
Fig. 15) but it remains generally lower than 50%; $ii)$ co-rotating conditions
and hence spin-up still occur; $iii)$ only the BHs in the highest mass bins
($\gtrsim 5\times 10^{9}\leavevmode\nobreak\ {\rm M}_{\sun}$) reach a fraction
$\gtrsim 0.5$ and exhibit systematically low spins ($a\simeq 0.2-0.3$). We
also note that above $M_{\rm BH}\gtrsim 10^{8}$ BHs undergo several mergers
(Fig. 12) and Eddington ratios are generally highly sub-Eddington (e.g. second
to last panel of Fig. 11). Therefore, mergers significantly contribute to the
variation of the spin. Moreover, mergers with misaligned directions tend to
significantly decrease $a$ (e.g. fourth panel of Fig. 11). From Fig. 13, we
infer that mergers and a more likely counter-rotating accretion are
responsible for the wider distribution of values of $a$.

In Fig. 13 we also note that BHs close to the seeding mass ($M_{\rm BH}\simeq
5.5\times 10^{5}\leavevmode\nobreak\ {\rm M}_{\sun}$) are characterised by a
steep increase of $a$ with mass. The region in the $a-M_{\rm BH}$ plane that
these BHs occupy is mostly determined by $M_{\rm BH,seed}$. However, the
initial value of $a$ does not affect the following evolution, since any
initial spin value is quickly evolved to maximal due to large accretion rates.

Overall, the trends in $a$ as a function of $M_{\rm BH}$ are consistent
between the zoom-in simulations (Fig. 10) and the Box4 (Fig.13), which even
have different resolutions. The distribution of $a$ is compatible with the
observations within the uncertainties, across the entire mass range.

###  5.3 Radiative efficiency

The radiative efficiency plays an important role. First of all, it enters the
computation of the Eddington accretion rate. Since the BH accretion rate
cannot exceed this rate in our model, the efficiency directly affects the
Eddington-limited growth phases. According to Eq. 29, the e-folding timescale

$\tau_{\rm M_{\rm BH}}=\frac{\epsilon_{\rm r}\tau_{\rm S}}{f_{\rm
Edd}(1-\epsilon_{\rm r})}$ (30)

depends on the efficiency. Eddington-limited co-rotating accretion on a
maximally spinning BH (i.e. $\epsilon_{\rm r}\sim 0.32$) implies $\tau_{\rm
M_{\rm BH}}\sim 210$ Myr, whereas in counter-rotating conditions
($\epsilon_{\rm r}\sim 0.038$) $\tau_{\rm M_{\rm BH}}\sim 18$ Myr. A lower
efficiency leads to a significantly faster BH growth. In our simulations,
Eddington-limited phases are characterised by maximally-spinning BHs and co-
rotating accretion, thus growth proceeds with $\epsilon_{\rm r}\sim 0.32$ (see
e.g. $\epsilon_{\rm r}$ panel in Fig. 11).

The efficiency also controls the amount of feedback energy released to the
surroundings, at a given $\dot{M}_{\rm BH}$ (see Eq. 27). A lower efficiency
implies less released energy and subsequent increased accretion. It also means
a faster increase in $M_{\rm BH}$ (because of the factor $1-\epsilon_{\rm r}$
in Eq. 10). This in turn boosts $\dot{M}_{\rm BH}$ (due to the $\propto M_{\rm
BH}$ dependence in Eddington-limited phases or $\propto M_{\rm BH}^{2}$
otherwise, see Eq. 1), leading to stronger feedback outbursts that hinder
accretion. In addition, since the Eddington accretion rate depends on the
efficiency, the switch to maintenance mode feedback is also affected. Overall,
the effects just discussed contribute in a non-trivial way to modify the
evolutionary path of a BH through the feedback loop. Our simulations – where
all these processes are self-consistently taken into account – show that BHs
are on the observed correlation between BH mass and stellar mass (Figs. 12 and
9). While the detailed evolutionary path along the plane $M_{\rm BH}-M_{\rm*}$
changes, each BH is eventually able to approach the correlation, implying that
the BHs still grow in an overall self-regulated scenario.

Fig. 17 highlights that $\epsilon_{\rm r}$ tends to decrease with increasing
mass at the high-mass end. We also find that the distribution of our simulated
sample is compatible with the distribution of empirical estimates by Daly
(2021), within the uncertainties. Note that the empirical sample is obtained
with a method that does not rely on a specific accretion disc model.
Therefore, while our sample has an upper and lower limit for the efficiency
determined by our theoretical assumption on the disc structure, the
interpretation of the empirical estimates is not bound to such an assumption.

##  6 Comparison with previous works

Dubois et al. (2014b) and Bustamante & Springel (2019) perform a statistical
study similar to ours. The former simulate a cosmological volume with size and
resolution comparable to Box4 and the latter a smaller box with higher
resolution ($L_{\rm box}=25\;h^{-1}\;\rm cMpc$ and $m_{\rm DM}=8.4\times
10^{6}\leavevmode\nobreak\ {\rm M}_{\sun}$). Our spin evolution model follows
Dubois et al. (2014b), although we track the spin on-the-fly whereas in their
work the spin is tracked in post-processing. $\epsilon_{\rm r}$ is assumed to
be fixed to 0.1. They include a thermal feedback channel for $f_{\rm
Edd}>0.01$ and a bipolar outflow launched in the direction of the local gas
otherwise, but both do not depend on $a$. Bustamante & Springel (2019) adopt
an on-the-fly spin update algorithm similar to ours, with a variable
$\epsilon_{\rm r}$ affecting the thermal feedback mode. However, below a mass-
dependent $f_{\rm Edd}$ threshold, BHs inject purely kinetic energy, with a
fixed efficiency and random direction, isotropic on average. Moreover, when
the sub-grid accretion disc is affected by self-gravity the angular momentum
direction of each accretion episode is extracted from a chosen angular
distribution. In our implementation we assume it is fixed and equal to the gas
angular momentum direction at all times. Interestingly, the trends we find in
Fig. 13 are in agreement with both works, regardless of whether or not
$\epsilon_{\rm r}$ depends on $a$ or whether the spin is evolved in post-
processing rather than on-the-fly. This indicates that the variability of the
radiative efficiency does not play a significant role in setting these trends.
On the other hand, $\epsilon_{\rm r}$ does affect AGN feedback and BH growth
(Sec. 2.4). As a result, a different prescription for $\epsilon_{\rm r}$
changes the detailed time evolution of BH mass and stellar mass, leading to a
different path towards the $M_{\rm BH}-M_{\rm*}$ relation and more scatter
around it (Bustamante & Springel 2019). Nonetheless, Fig. 12 shows that the
self-regulated scenario is still present and eventually leads the BHs on the
relation.

We find that the dynamical state of the feeding material, combined with the
key parameter $J_{\rm d}/J_{\rm BH}$, plays a significant role in the
evolution of the BH spin due to gas accretion (Fig. 15 and 16). Dotti et al.
(2013) conclude that when $J_{\rm d}/J_{\rm BH}\ll 1$ (as it occurs in our
simulations, see Fig. 15), the feeding angular momentum distribution
determines whether $a$ increases or decreases. Our findings regarding the
behaviour of the spin in response to the feeding conditions are in agreement
with theirs, as well as with Dubois et al. (2014b) and Bustamante & Springel
(2019). When the gas angular momentum direction varies slowly and its
misalignment with respect to the BH spin is small (indicating a preferential
direction), the BHs are maximally spinning. This generally occurs at $M_{\rm
BH}\lesssim 10^{8}\leavevmode\nobreak\ {\rm M}_{\sun}$ (see Fig. 11, 14, 15).
Conversely, uncorrelated gas angular momentum directions lead to increasingly
probable counter-rotating accretion and BH spin-down (bottom panel of Fig. 15,
Sec. 4.4). We note that the large scatter in $a$ at high masses (Fig. 13) is
found also by Dubois et al. (2014a) and Bustamante & Springel (2019), despite
the different prescriptions for AGN feedback. This might indicate that its
effect is to generally induce loss of coherence in the angular momentum
distribution (middle panel of Fig. 15), regardless of the specific channel.

In contrast to our model, Bustamante & Springel (2019) introduce stochasticity
in $\bm{j}_{\rm d}$ (Eq. 16) at the sub-resolution level, on top of the
resolved variability in the simulation. They assume that in the self-gravity
regime $\bm{j}_{\rm d}$ is extracted from a distribution that ranges from
random (isotropic) to concentrated around the preferential axis set by the
local gas angular momentum (anisotropic). As a result, they observe a more
pronounced decoupling between $\bm{j}_{\rm d}$ and $\bm{j}_{\rm BH}$ compared
to us, in case of an isotropic distribution. On the other hand, they find the
same widening of the distribution of $a$ at high BH masses, regardless of the
degree of anisotropy. Dubois et al. (2014a) also explore the effect of varying
the distribution of $\bm{j}_{\rm d}$ at the sub-resolution level. However,
they introduce stochasticity in all accretion regimes. They find that even a
slight anisotropy leads to results that are very similar to the completely
coherent case (i.e. when the gas preserves the angular momentum direction as
measured by the simulation). Only if the gas angular momentum is random
oriented at all masses, then BHs settle on $a\sim 0.2-0.3$ (King et al. 2008)
and an increasing trend of spin with mass is found (Fanidakis et al. 2011) (in
those conditions the trend is produced by mergers, that tend to bring slowly-
spinning BHs to values around $0.7$ for $M_{\rm BH}\gtrsim
10^{9}\leavevmode\nobreak\ {\rm M}_{\sun}$). One possibility to explain these
results is that when gas accretion is sub-dominant, the trend is mainly driven
by mergers. However, we observe that above $M_{\rm BH}\gtrsim
10^{8}\leavevmode\nobreak\ {\rm M}_{\sun}$ gas accretion does contribute to
spin evolution, inducing both spin-up and spin-down (Fig. 14). The net effect
depends crucially on the accretion rate, the BH environment and the level of
anisotropy of the feeding gas. This level varies widely over time and across
the BH population. Mergers also contribute to modify the spin value, but the
relative contribution of gas accretion and mergers to spin evolution is
different depending on the detailed cosmological history. We postpone a
systematic study to a future work.

Finally, we note that we do not integrate the full differential equation that
describes the precession and alignment process, due to the coupling between
the gas distribution in the accretion disc and the BH spin (at variance with,
e.g., Fiacconi et al. 2018). Rather, we assume that the process is discretised
in accretion episodes and is globally taken into account, while the total
angular momentum is conserved (Sec. 2.2), following King et al. (2005).
Fiacconi et al. (2018) developed an algorithm that includes the full detailed
treatment, although it requires high temporal and spatial resolution. In fact,
timesteps as low as $10^{-3}$ Myr are required in some cases to meaningfully
integrate the differential equation, which are prohibitive in a full
cosmological context. Moreover, their model measures directly the inflow
properties as resolved by the simulation at the sub-resolution boundary,
rather than using an effective prescription such as the Bondi parametrisation.
Therefore, it is suitable for high-resolution simulations (e.g. $\sim$ pc
scale). In contrast, we assume that the mass rate through the sub-resolution
accretion disc is equal to the mass accretion rate onto the BH and it is
identical to the Bondi accretion rate at all times. We further note that
Fiacconi et al. (2018) model spin evolution only due to gas accretion,
although mergers are expected to contribute to spin evolution for $M_{\rm
BH}\gtrsim 10^{8}\leavevmode\nobreak\ {\rm M}_{\sun}$ (Fanidakis et al. 2011;
Dubois et al. 2014a). Although they do not produce a statistical sample of
BHs, they perform a suite of simulations aimed at mimicking a range of
realistic conditions. Despite the approach is different in a number of
numerical aspects, the expected effect of the gas accretion channel on spin
evolution is in line with ours: systematic spin-up for $M_{\rm BH}\lesssim
10^{7}\leavevmode\nobreak\ {\rm M}_{\sun}$ and a wider distribution of $a$ at
higher masses.

##  7 Conclusions

We implement a sub-resolution model to track the evolution of BH spins due to
gas accretion and mergers in large scale cosmological simulations. The model
assumes the presence of a misaligned thin accretion disc perturbed by the
metric of a spinning BH, which in turn experiences a torque that modifies its
spin direction. The BH radiative efficiency and Eddington accretion rate are
dependent on the BH spin and thus variable across cosmic time. Their impact on
accretion and feedback is therefore captured self-consistently. We design a
simulation suite featuring idealised, isolated systems (Sec. 4.1 and 4.2) to
validate the model and cosmological setups (Sec. 4.3 and 4.4) to investigate
statistical properties of the BH population. We summarise our findings as
follows.

  * •

The ability of a single accretion episode to modify the BH spin depends on the
amount of mass and angular momentum it carries with respect to the BH ($J_{\rm
d}/J_{\rm BH}$, Fig. 16). An accretion episode with larger $J_{\rm d}/J_{\rm
BH}$ induces the BH spin direction $\bm{j}_{\rm BH}$ to tilt more towards the
gas angular momentum direction $\bm{j}_{\rm g}$.

  * •

The evolution of the direction and the magnitude of the spin are tightly
coupled. When $J_{\rm d}/J_{\rm BH}\lesssim 1$ (Fig. 15), the feeding
distribution of the gas angular momentum directions determines whether $a$
preferably increases or decreases. If accretion occurs consistently along the
same plane (e.g. Fig. 5), spin-up is expected. Conversely, if $\bm{j}_{\rm g}$
changes direction erratically, counter-rotating conditions and spin-down can
occur (Fig. 11 and 15).

  * •

In a cosmological context, we identify two regimes, depending on the
distribution of $a$ with BH mass (Fig. 13). BHs with $M_{\rm BH}\lesssim
2\times 10^{7}\leavevmode\nobreak\ {\rm M}_{\sun}$ tend to be highly spinning
($a\gtrsim 0.85$). At the high-mass range ($M_{\rm BH}\gtrsim 2\times
10^{7}\leavevmode\nobreak\ {\rm M}_{\sun}$), $a$ exhibits a broad range of
values.

  * •

We observe a wide variety of evolutionary histories of $a$ (Fig. 14),
depending on the dynamical state of the gas feeding the BH and the occurrence
of coalescences. This indicates that the level of anisotropy of $\bm{j}_{\rm
g}$ and the relative contribution of mergers and accretion varies across the
BH population.

  * •

When $\bm{j}_{\rm g}$ exhibits some degree of coherence and varies slowly
(generally at $z\gtrsim 2$ and $M_{\rm BH}\lesssim 2\times
10^{7}\leavevmode\nobreak\ {\rm M}_{\sun}$), $\bm{j}_{\rm BH}$ follows the
average evolution of $\bm{j}_{\rm g}$ and small misalignment is observed (Fig.
15). At late times ($z\lesssim 2$), $\bm{j}_{\rm g}$ shows large and abrupt
changes, while $\bm{j}_{\rm BH}$ is stable over long periods (hundreds of Myr,
e.g. Fig. 11). Indeed, since $J_{\rm d}/J_{\rm BH}\lesssim 1$, accretion
episodes are not able to modify significantly $\bm{j}_{\rm BH}$. This
frequently leads to large misalignment.

  * •

The tendency for maximal spin in the low-mass range is due to accretion of co-
rotating gas with small misalignment (Fig. 15). The wide range of values of
$a$ in the high-mass range is due to mergers and more isotropically
distributed $\bm{j}_{\rm g}$, this results in a probability of counter-
rotating accretion that increases with mass (Figs. 14 and 15). The
distribution of $\bm{j}_{\rm g}$ arises self-consistently, as measured from
the simulation.

  * •

Including a self-consistent $\epsilon_{\rm r}$ has an important effect in
determining the BH growth rate in the Eddington-limited phases. A higher
efficiency during these phases implies a lower growth rate. Our statistical
sample shows that BHs with $M_{\rm BH}\lesssim 4\times
10^{7}\leavevmode\nobreak\ {\rm M}_{\sun}$ have efficiencies always around
$0.32$, whereas the most massive BHs have generally lower efficiencies (Fig.
17).

  * •

Although the variable $\epsilon_{\rm r}$ modifies the detailed path on the
$M_{\rm BH}-M_{\rm*}$ plane, BHs eventually approach the observed correlation,
indicating self-regulated growth.

  * •

The spin with which BHs are initialised is erased quickly after they are
seeded. Massive BHs can retain some information on the dynamical state of the
gas they recently accreted (Figs. 14 and 15).

We caution that so far we have coupled the spin evolution model to a single
channel of energy injection, namely purely thermal, aimed at reproducing the
radiative feedback. An additional mechanism, in which large-scale jets
(10s-100s of kpc) are the key actors, is thought to be specifically relevant
at the high-mass end. We plan to include this feature in future works, by
coupling the spin evolution model to a feedback channel from jets, dependent
on the spin magnitude and direction.

###### Acknowledgements.

LS acknowledges support from ‘BiD4BEST’ - European Innovative Training Network
(ITN) funded by the Marie Skłodowska-Curie Actions (860744) in Horizon 2020.
MV is supported by the Fondazione ICSC, Spoke 3 Astrophysics and Cosmos
Observations, National Recovery and Resilience Plan (PNRR) Project ID
CN_00000013 ”Italian Research Center on High-Performance Computing, Big Data
and Quantum Computing” funded by MUR (Missione 4, Componente 2, Investimento
1.4: Potenziamento strutture di ricerca e creazione di campioni nazionali di
R&S, M4C2-19) - Next Generation EU. MV also acknowledges partial financial
support from the INFN Indark Grant. This work was supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy - EXC-2094 - 390783311 and we are especially grateful for
the support through the Computational Center for Particle and Astrophysics
(C2PAP). Some of the analyses reported in this work have been performed using
pynbody (Pontzen et al. 2013)

## References

  * Bardeen & Petterson (1975) Bardeen, J. M. & Petterson, J. A. 1975, ApJ, 195, L65
  * Beck et al. (2016) Beck, A. M., Murante, G., Arth, A., et al. 2016, MNRAS, 455, 2110
  * Beckmann et al. (2019) Beckmann, R. S., Dubois, Y., Guillard, P., et al. 2019, A&A, 631, A60
  * Benson et al. (2003) Benson, A. J., Bower, R. G., Frenk, C. S., et al. 2003, ApJ, 599, 38
  * Berti & Volonteri (2008) Berti, E. & Volonteri, M. 2008, ApJ, 684, 822
  * Blandford & Znajek (1977) Blandford, R. D. & Znajek, R. L. 1977, MNRAS, 179, 433
  * Bollati et al. (2023) Bollati, F., Lupi, A., Dotti, M., & Haardt, F. 2023, MNRAS, 520, 3696
  * Bonafede et al. (2011) Bonafede, A., Dolag, K., Stasyszyn, F., Murante, G., & Borgani, S. 2011, MNRAS, 418, 2234
  * Bondi (1952) Bondi, H. 1952, MNRAS, 112, 195
  * Bondi & Hoyle (1944) Bondi, H. & Hoyle, F. 1944, MNRAS, 104, 273
  * Booth & Schaye (2009) Booth, C. M. & Schaye, J. 2009, MNRAS, 398, 53
  * Broderick et al. (2022) Broderick, A. E., Tiede, P., Pesce, D. W., & Gold, R. 2022, ApJ, 927, 6
  * Bustamante & Springel (2019) Bustamante, S. & Springel, V. 2019, MNRAS, 490, 4133
  * Cattaneo et al. (2009) Cattaneo, A., Faber, S. M., Binney, J., et al. 2009, Nature, 460, 213
  * Cenci et al. (2020) Cenci, E., Sala, L., Lupi, A., Capelo, P. R., & Dotti, M. 2020, MNRAS, 500, 3719
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Chael et al. (2023) Chael, A., Lupsasca, A., Wong, G. N., & Quataert, E. 2023, ApJ, 958, 65
  * Churazov et al. (2005) Churazov, E., Sazonov, S., Sunyaev, R., et al. 2005, MNRAS, 363, L91
  * Cielo et al. (2018) Cielo, S., Babul, A., Antonuccio-Delogu, V., Silk, J., & Volonteri, M. 2018, A&A, 617, A58
  * Croton et al. (2006) Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11
  * Daly (2009) Daly, R. A. 2009, ApJ, 696, L32
  * Daly (2011) Daly, R. A. 2011, MNRAS, 414, 1253
  * Daly (2019) Daly, R. A. 2019, ApJ, 886, 37
  * Daly (2021) Daly, R. A. 2021, MNRAS, 500, 215
  * David et al. (2009) David, L. P., Jones, C., Forman, W., et al. 2009, ApJ, 705, 624
  * Davis & Tchekhovskoy (2020) Davis, S. W. & Tchekhovskoy, A. 2020, ARA&A, 58, 407
  * Dehnen & Aly (2012) Dehnen, W. & Aly, H. 2012, MNRAS, 425, 1068
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604
  * Dolag et al. (2009) Dolag, K., Borgani, S., Murante, G., & Springel, V. 2009, MNRAS, 399, 497
  * Dolag et al. (2004) Dolag, K., Jubelgas, M., Springel, V., Borgani, S., & Rasia, E. 2004, ApJ, 606, L97
  * Dolag et al. (2005) Dolag, K., Vazza, F., Brunetti, G., & Tormen, G. 2005, MNRAS, 364, 753
  * Dong-Páez et al. (2023) Dong-Páez, C. A., Volonteri, M., Beckmann, R. S., et al. 2023, A&A, 673, A120
  * Donnert et al. (2013) Donnert, J., Dolag, K., Brunetti, G., & Cassano, R. 2013, MNRAS, 429, 3564
  * Dotti et al. (2013) Dotti, M., Colpi, M., Pallini, S., Perego, A., & Volonteri, M. 2013, ApJ, 762, 68
  * Dubois et al. (2021) Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109
  * Dubois et al. (2014a) Dubois, Y., Volonteri, M., & Silk, J. 2014a, MNRAS, 440, 1590
  * Dubois et al. (2014b) Dubois, Y., Volonteri, M., Silk, J., Devriendt, J., & Slyz, A. 2014b, MNRAS, 440, 2333
  * Eddington (1916) Eddington, A. S. 1916, MNRAS, 77, 16
  * Fabian (1999) Fabian, A. C. 1999, MNRAS, 308, L39
  * Fabian (2012) Fabian, A. C. 2012, ARA&A, 50, 455
  * Fabian et al. (2011) Fabian, A. C., Sanders, J. S., Allen, S. W., et al. 2011, MNRAS, 418, 2154
  * Fanidakis et al. (2011) Fanidakis, N., Baugh, C. M., Benson, A. J., et al. 2011, MNRAS, 410, 53
  * Fiacconi et al. (2018) Fiacconi, D., Sijacki, D., & Pringle, J. E. 2018, MNRAS, 477, 3807
  * Forman et al. (2007) Forman, W., Jones, C., Churazov, E., et al. 2007, ApJ, 665, 1057
  * Frank et al. (2002) Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics: Third Edition
  * Ghisellini et al. (2010) Ghisellini, G., Tavecchio, F., Foschini, L., et al. 2010, MNRAS, 402, 497
  * Giommi et al. (2012) Giommi, P., Padovani, P., Polenta, G., et al. 2012, MNRAS, 420, 2899
  * Griffin et al. (2019) Griffin, A. J., Lacey, C. G., Gonzalez-Perez, V., & Lagos, C. d. P. 2019, arXiv e-prints, arXiv:1912.09490
  * Groth et al. (2023) Groth, F., Steinwandel, U. P., Valentini, M., & Dolag, K. 2023, MNRAS, 526, 616
  * Haardt & Madau (2001) Haardt, F. & Madau, P. 2001, in Clusters of Galaxies and the High Redshift Universe Observed in X-rays, ed. D. M. Neumann & J. T. V. Tran, 64
  * Harrison et al. (2018) Harrison, C. M., Costa, T., Tadhunter, C. N., et al. 2018, Nature Astronomy, 2, 198
  * Hawley & Krolik (2006) Hawley, J. F. & Krolik, J. H. 2006, ApJ, 641, 103
  * Hirschmann et al. (2014) Hirschmann, M., Dolag, K., Saro, A., et al. 2014, MNRAS, 442, 2304
  * Horton et al. (2020) Horton, M., Krause, M., & Hardcastle, M. 2020, MNRAS, 499, 5765
  * Hoyle & Lyttleton (1939) Hoyle, F. & Lyttleton, R. A. 1939, Math. Proc. Camb. Phil. Soc., 35, 405
  * Huško et al. (2022) Huško, F., Lacey, C. G., Schaye, J., Schaller, M., & Nobels, F. S. J. 2022, Spin-Driven Jet Feedback in Idealised Simulations of Galaxy Groups and Clusters
  * Izquierdo-Villalba et al. (2020) Izquierdo-Villalba, D., Bonoli, S., Dotti, M., et al. 2020, MNRAS, 495, 4681
  * Jiang et al. (2022) Jiang, J., Abdikamalov, A. B., Bambi, C., & Reynolds, C. S. 2022, MNRAS, 514, 3246
  * Jiang et al. (2019) Jiang, J., Fabian, A. C., Dauser, T., et al. 2019, MNRAS, 489, 3436
  * Karademir et al. (2019) Karademir, G. S., Remus, R.-S., Burkert, A., et al. 2019, MNRAS, 487, 318
  * Karakas & Lattanzio (2007) Karakas, A. & Lattanzio, J. C. 2007, PASA, 24, 103
  * King et al. (2005) King, A. R., Lubow, S. H., Ogilvie, G. I., & Pringle, J. E. 2005, MNRAS, 363, 49
  * King et al. (2008) King, A. R., Pringle, J. E., & Hofmann, J. A. 2008, MNRAS, 385, 1621
  * Kormendy & Ho (2013) Kormendy, J. & Ho, L. C. 2013, ARA&A, 51, 511
  * Lagos et al. (2009) Lagos, C. D. P., Padilla, N. D., & Cora, S. A. 2009, MNRAS, 395, 625
  * Lalakos et al. (2022) Lalakos, A., Gottlieb, O., Kaaz, N., et al. 2022, ApJL, 936, L5
  * Lasota et al. (2014) Lasota, J. P., Gourgoulhon, E., Abramowicz, M., Tchekhovskoy, A., & Narayan, R. 2014, Phys. Rev. D, 89, 024041
  * Liska et al. (2019) Liska, M., Tchekhovskoy, A., Ingram, A., & van der Klis, M. 2019, MNRAS, 487, 550
  * Maio et al. (2013) Maio, U., Dotti, M., Petkova, M., Perego, A., & Volonteri, M. 2013, ApJ, 767, 37
  * Mallick et al. (2022) Mallick, L., Fabian, A. C., García, J. A., et al. 2022, MNRAS, 513, 4361
  * Martin et al. (2007) Martin, R. G., Pringle, J. E., & Tout, C. A. 2007, MNRAS, 381, 1617
  * Martínez-Sansigre & Rawlings (2011) Martínez-Sansigre, A. & Rawlings, S. 2011, MNRAS, 414, 1937
  * McConnell & Ma (2013) McConnell, N. J. & Ma, C.-P. 2013, ApJ, 764, 184
  * McKinney et al. (2012) McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083
  * McNamara & Nulsen (2007) McNamara, B. R. & Nulsen, P. E. J. 2007, ARA&A, 45, 117
  * McNamara & Nulsen (2012) McNamara, B. R. & Nulsen, P. E. J. 2012, New Journal of Physics, 14, 055023
  * Nakamura et al. (2008) Nakamura, M., Tregillis, I. L., Li, H., & Li, S. 2008, ApJ, 686, 843
  * Narayan et al. (2022) Narayan, R., Chael, A., Chatterjee, K., Ricarte, A., & Curd, B. 2022, MNRAS, 511, 3795
  * Natarajan & Armitage (1999) Natarajan, P. & Armitage, P. J. 1999, MNRAS, 309, 961
  * Natarajan & Pringle (1998) Natarajan, P. & Pringle, J. E. 1998, ApJ, 506, L97
  * Nomoto et al. (2013) Nomoto, K., Kobayashi, C., & Tominaga, N. 2013, ARA&A, 51, 457
  * Ogilvie (1999) Ogilvie, G. I. 1999, MNRAS, 304, 557
  * Padovani & Matteucci (1993) Padovani, P. & Matteucci, F. 1993, ApJ, 416, 26
  * Papaloizou & Pringle (1983) Papaloizou, J. C. B. & Pringle, J. E. 1983, MNRAS, 202, 1181
  * Perego et al. (2009) Perego, A., Dotti, M., Colpi, M., & Volonteri, M. 2009, MNRAS, 399, 2249
  * Pontzen et al. (2013) Pontzen, A., Roškar, R., Stinson, G., & Woods, R. 2013, pynbody: N-Body/SPH analysis for python, Astrophysics Source Code Library, record ascl:1305.002
  * Pringle (1981) Pringle, J. E. 1981, Annu. Rev. Astron. Astrophys., 19, 137
  * Reynolds (2021) Reynolds, C. S. 2021, Annu. Rev. Astron. Astrophys., 59, 117
  * Rezzolla et al. (2008a) Rezzolla, L., Barausse, E., Dorband, E. N., et al. 2008a, Phys. Rev. D, 78, 044002
  * Rezzolla et al. (2008b) Rezzolla, L., Dorband, E. N., Reisswig, C., et al. 2008b, ApJ, 679, 1422
  * Sala et al. (2021) Sala, L., Cenci, E., Capelo, P. R., Lupi, A., & Dotti, M. 2021, MNRAS, 500, 4788
  * Sanders et al. (2009) Sanders, J. S., Fabian, A. C., & Taylor, G. B. 2009, MNRAS, 396, 1449
  * Scheuer & Feiler (1996) Scheuer, P. & Feiler, R. 1996, MNRAS, 282, 291
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
  * Sijacki et al. (2007) Sijacki, D., Springel, V., Di Matteo, T., & Hernquist, L. 2007, MNRAS, 380, 877
  * Sikora et al. (2007) Sikora, M., Stawarz, Ł., & Lasota, J.-P. 2007, ApJ, 658, 815
  * Silk & Rees (1998) Silk, J. & Rees, M. J. 1998, A&A, 331, L1
  * Springel (2005) Springel, V. 2005, MNRAS, 364, 1105
  * Springel et al. (2005a) Springel, V., Di Matteo, T., & Hernquist, L. 2005a, MNRAS, 361, 776
  * Springel & Hernquist (2003) Springel, V. & Hernquist, L. 2003, MNRAS, 339, 289
  * Springel et al. (2005b) Springel, V., White, S. D. M., Jenkins, A., et al. 2005b, Nature, 435, 629
  * Springel et al. (2001) Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726
  * Steinborn et al. (2015) Steinborn, L. K., Dolag, K., Hirschmann, M., Prieto, M. A., & Remus, R.-S. 2015, MNRAS, 448, 1504
  * Steinwandel et al. (2019) Steinwandel, U. P., Beck, M. C., Arth, A., et al. 2019, MNRAS, 483, 1008
  * Talbot et al. (2021) Talbot, R. Y., Bourne, M. A., & Sijacki, D. 2021, MNRAS, 504, 3619
  * Talbot et al. (2023) Talbot, R. Y., Sijacki, D., & Bourne, M. A. 2023, arXiv e-prints, arXiv:2306.07316
  * Tamburini et al. (2020) Tamburini, F., Thidé, B., & Della Valle, M. 2020, MNRAS, 492, L22
  * Tchekhovskoy (2015) Tchekhovskoy, A. 2015, in Astrophysics and Space Science Library, Vol. 414, The Formation and Disruption of Black Hole Jets, ed. I. Contopoulos, D. Gabuzda, & N. Kylafis, 45
  * Tchekhovskoy et al. (2011) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79
  * Thielemann et al. (2003) Thielemann, F. K., Argast, D., Brachwitz, F., et al. 2003, in From Twilight to Highlight: The Physics of Supernovae, ed. W. Hillebrandt & B. Leibundgut, 331
  * Thorne (1974) Thorne, K. S. 1974, ApJ, 191, 507
  * Tornatore et al. (2007) Tornatore, L., Borgani, S., Dolag, K., & Matteucci, F. 2007, MNRAS, 382, 1050
  * Valentini et al. (2020) Valentini, M., Murante, G., Borgani, S., et al. 2020, MNRAS, 491, 2779
  * Volonteri et al. (2005) Volonteri, M., Madau, P., Quataert, E., & Rees, M. J. 2005, ApJ, 620, 69
  * Volonteri & Rees (2005) Volonteri, M. & Rees, M. J. 2005, ApJ, 633, 624
  * Volonteri et al. (2007) Volonteri, M., Sikora, M., & Lasota, J.-P. 2007, ApJ, 667, 704
  * Wiersma et al. (2009) Wiersma, R. P. C., Schaye, J., Theuns, T., Dalla Vecchia, C., & Tornatore, L. 2009, MNRAS, 399, 574
  * Wilkins (1972) Wilkins, D. C. 1972, Phys. Rev. D, 5, 814
  * Zhang & Lu (2019) Zhang, X. & Lu, Y. 2019, ApJ, 873, 101