Black hole spin evolution in warped accretion discs

Massive black holes (BHs) inhabiting galactic nuclei can be described by two parameters only, i.e. mass and spin, that change through cosmic time in response to accretion and merger events. While most numerical simulations accurately track the BH mass, spin evolution is rarely taken into account. In this work, we implement and validate a self-consistent sub-grid model for the evolution of the BH mass and spin via gas accretion in the hydrodynamics code GIZMO. The model assumes that accretion from resolved scales does not occur instantaneously, but is mediated by a sub-grid geometrically thin $\alpha$-disc. After validating our model semi-analytically, we test it in an idealized environment consisting of a circumnuclear disc, where gas accretion onto the accretion disc is consistently determined by GIZMO. In the absence of any accretion-related feedback, the spin evolution closely traces that observed in the semi-analytical models and depends on the free parameters of our implementation, such as the initial BH spin, angular momentum of the accretion disc, and the radius at which the gas inflow circularises. In GIZMO, we also couple our model with the biconical-outflow model presented in a companion paper, wherein the feedback axis is always aligned with the BH spin. In this last case, the evolution of the central BH differs significantly from the previous cases, since the feedback process modifies the gas dynamics and its inflow rates from resolved scales. Such an interaction cannot be modeled by simple semi-analytical models and should be treated using full N-body hydrodynamical simulations.


INTRODUCTION
According to the no-hair conjecture (Israel 1967(Israel , 1968Carter 1971;Hawking 1972;Robinson 1975), massive black holes (BHs) are thought to be completely characterized by three parameters: mass, charge, and spin. Since any electric charge would be quickly neutralized by charges in the surrounding plasma, astrophysical BHs can be described only in terms of their masses and spins.
Spins in particular have a number of fundamental consequences on the evolution of the BHs. The spin magnitude affects the position of the innermost stable circular orbit (ISCO; Bardeen 1970; Bardeen et al. 1972) and, as a consequence, the radiative efficiency of accretion processes 1 and the rate at which BHs can grow in time. Furthermore, ac-E-mail: e.cenci@campus.unimib.it 1 The dependence of the radiative efficiency on the spin is weaker for supercritically accreting BHs, see, e.g. Madau et al. (2014).
cording to the spin paradigm, high spins are responsible for the launching of the relativistic jets observed in active galactic nuclei (AGN) over the whole spectral range (Blandford & Znajek 1977).
The spin direction as well has been proposed to play a central role in the growth of the BHs and in the co-evolution with their host galaxies. As an example, feedback exerted by the accreting BH on to its environment has been suggested as a possible cause for the observed BH-host galaxy scaling relations (Kormendy & Ho 2013) and for the quenching of star formation in massive galaxies (Silk & Rees 1998;Fabian 1999Fabian , 2012. It has been argued that, to achieve such effects, the spin direction 2 must evolve in time (see, e.g. Nayakshin et al. 2012 for the galaxy-scale feedback from a wind and 2 The feedback reference direction is expected to be parallel to the BH spin both for a Blandford & Znajek (1977) jet as well as for a wind launched from the very central regions of an accretion disc (Bardeen & Petterson 1975). Cielo et al. 2018 for the feedback on to the intra-cluster medium mediated by a jet).
Finally, spin magnitudes and directions determine the recoil velocity that the remnant of a BH merger experiences due to anisotropic gravitational-wave emission (e.g. Koppitz et al. 2007), with such velocity being well above typical escape velocities from massive galaxies for some spin configurations (e.g. Campanelli et al. 2007;Baker et al. 2008;Herrmann et al. 2007;Schnittman & Buonanno 2007;Lousto & Zlochower 2011;Lousto et al. 2012).
The spin direction has a direct impact on the evolution of the spin magnitude. Prolonged accretion events with a fixed geometry result in maximally spinning BHs after a mass growth comparable to the initial BH mass (Bardeen 1970). On the contrary, episodic accretion events isotropically oriented and each having an accreted mass significantly smaller than the initial mass of the BH will on average decrease the BH spin, due to the larger size of the ISCO and therefore a larger magnitude of the angular momentum per unit of mass associated to retrograde accretion events (King & Pringle 2006). Dotti et al. (2013) showed that the two abovementioned accretion modes, often dubbed "coherent" and "chaotic", are the extremes of all the possible accretion configurations, and that, depending on the assumptions on the BH fueling geometry, the expected spin magnitudes for any given BH mass can seamlessly vary from 0 to ∼1, in agreement with the results of early N -body hydrodynamical simulations with the spin of the central BH being evolved in post-processing (Dotti et al. 2010;Maio et al. 2013) or onthe-fly (Dubois et al. 2014a,b). By implementing the model by Dotti et al. (2013) in a pre-existing semi-analytical galaxy formation model (Barausse 2012), and by assuming that the gas reservoir for the BH growth has the same dynamical properties of the host galaxy nuclei at 100 pc scales, 3 Sesana et al. (2014) managed for the first time to reproduce the observational constraints on BH spins available, without the introduction of any additional freely-tunable parameter.
The above-cited studies about accretion-driven spin evolution did not consider the possible effect of AGN feedback on to the dynamics of the gas fueling the accretion process. As an example, a prolonged accretion event whose accretion disc is initially misaligned with respect to the BH spin direction by more than π/2 would tend to re-align the BH spin with the gas angular momentum (Bardeen & Petterson 1975). If, however, the accretion process triggers a directional feedback aligned with the BH spin, the feedback could significantly alter the dynamics of the gas reservoir, modifying the following spin evolution in a non-linear fashion. This is a severe limitation of semi-analytical studies, that cannot follow in real time the impact of the spin evolution on the larger-scale gas dynamics. For this reason, we hereby present a new implementation for the coupled evolution of BH spins, unresolved accretion discs, and directional feedback in N -body, hydrodynamical simulations. Our implementation includes the spin evolution discussed in Fiacconi et al. (2018), that relaxes the small-warp approxima-tion (Scheuer & Feiler 1996;Martin et al. 2007;Perego et al. 2009;Dotti et al. 2013), the directional feedback presented in a companion paper (Sala et al. 2021), and a new sub-grid model for the self-consistent evolution of unresolved accretion discs around BHs. The model has been implemented in the publicly available code gizmo 4 (Hopkins 2015).
The manuscript is structured as follows: the new model for the unresolved accretion disc and spin evolution is described in Section 2, and its semi-analytical validation (in the absence of feedback) is presented in Section 3. The initial conditions of the idealized run used to test our full implementation are discussed in Section 4.1. The results of the tests and their discussion are presented in Sections 4 and 5, respectively.

THE SPIN EVOLUTION MODEL
In this section, we introduce our spin evolution model and validate it in a semi-analytic controlled environment. In our model, a BH particle represents a sub-grid system composed of a BH surrounded by a standard α-disc (Shakura & Sunyaev 1973).

Model description
The BH is described by its mass MBH = 10 6 MBH,6 M , angular momentum direction jBH = JBH/JBH, and spin parameter 5 aBH = cJBH/(GM 2 BH ), where JBH and JBH are the BH's angular momentum vector and magnitude, respectively, c is the speed of light in vacuum, and G is the gravitational constant.
To characterize the accretion disc, we need to specify its mass M disc = 10 4 M disc,4 M , total angular momentum J disc = j disc J disc (where j disc and J disc are the angular momentum direction and magnitude, respectively), and the accretion rate on to the BHṀacc−BH. The Eddington ratio can be expressed by f Edd =Ṁacc−BH/Ṁ Edd , wherė M Edd = 4πGMBHmp/(σTηc) is the Eddington accretion rate, mp is the proton mass, η = 0.1η0.1 is the radiative efficiency, and σT is the Thomson cross-section. In order to initialise the accretion disc properties, we set j disc = jgas, where jgas is the angular momentum direction of the gas surrounding the BH-disc system.
Following the Shakura & Sunyaev (1973) solution for the external regions of the disc, the radial viscosity ν1 scales with the cylindrical radius R as (Martin et al. 2007) where β = 3/4. The normalization is computed as (Frank et al. 2002;Perego et al. 2009 where α = 0.1α0.1 is related to the definition of ν1 (via ν1 ≡ αcsH), cs is the gas sound speed, and H is the disc scale-height. We compute the vertical viscosity ν2 under the approximation of a constant viscosity ratio ν2/ν1 = ξα −2 /2 (Lodato & Pringle 2007). In order to ensure that the adopted analytical prescription for viscosity is valid, we assume that vertical perturbations in the disc propagate diffusively, the latter translating in α > H/R (Pringle 1992). For our tests, we set α = 0.1 inside the disc, in agreement with the value inferred by available observational data for thin accretion discs (King et al. 2007), and ξ = 0.7 (Lodato & Pringle 2007;Perego et al. 2009). By equating the Lense-Thirring (Lense & Thirring 1918) precession timescale Ω −1 LT = c 2 R 3 /(2GJBH) to the characteristic time-scale τν 2 ∼ R 2 /ν2 of the vertical perturbation propagation in the disc out to a radius R (Lodato & Pringle 2006), we obtain the characteristic extension of the warp Rwarp, that divides the disc in two regions: the inner one (R Rwarp), where the disc and BH angular momenta are expected to be parallel, and the outer one (R Rwarp), where the disc maintains its original inclination (j disc ∼ jgas). For the adopted viscosity model, we obtain

Rwarp
Rg where Rg = GMBH/c 2 is the BH gravitational scale-radius. The radiative efficiency η in the disc is determined selfconsistently from the properties at the equatorial ISCO of a Novikov & Thorne (1973) disc (Bardeen 1970;Bardeen et al. 1972), under the assumption that the inner part of the disc is aligned with the BH equatorial plane due to the Bardeen-Petterson effect (Bardeen & Petterson 1975): where RISCO is the ISCO radius, computed as in Bardeen et al. (1972). Assigning a viscosity prescription allows to compute the steady-state solution for the disc mass and angular momentum density profiles (i.e. per unit surface). The disc total angular momentum is defined as the integrated angular momentum density out to an orbit with R = Rout that encompasses a mass equal to the disc total mass M disc . Assuming that the main contribution to the disc total angular momentum comes from the outer regions of the disc, we can consistently derive the total angular momentum magnitude as (Fiacconi et al. 2018 Once the BH-disc system is initialised, the time evolution is modelled according to the prescriptions of Fiacconi et al. (2018), consistently with the properties of the largescale resolved environment. At every time-step ∆t, we compute the Eddington ratio f Edd by means of Equation (5). To ensure sub-Eddington accretion, as required by the assumed disc model, we enforce f Edd ≤ 1. The standard thin-disc model is not supposed to hold in low-accretion regimes. In those cases, one should instead implement a more appropriate disc model (e.g. an advection-dominated accretion flow; Narayan & Yi 1994). However, we do not impose any lower limit to f Edd (as done in, e.g. Fiacconi et al. 2018), because the validity threshold is not well constrained (see, e.g. Chen 1995) and the BH spin evolution is negligible for very small accretion rates.
The BH and disc masses are updated at every timestep taking into account the disc draining rateṀacc−BH, radiative losses (through η), the mass accretion rate on to the discṀin given by the dynamics of the resolved-scales gas, and the mass outflow rateṀout: The BH angular momentum evolves in response to both the gas accretion at R = RISCO and the gravito-magnetic torque exerted by the material flowing through R ∼ Rwarp. While any process taking place at the ISCO only modifies the magnitude of the BH angular momentum, the Bardeen-Petterson torque, in steady-state warped conditions, is responsible for its change in direction (King et al. 2005;Fiacconi et al. 2018): Here ΛISCO is the specific angular momentum per unit mass carried by the gas flowing at the ISCO, and τgm is the characteristic time-scale over which the gravito-magnetic interaction significantly changes the BH spin direction (Martin et al. 2007;Lodato & Pringle 2006;Perego et al. 2009;Dotti et al. 2013 During the evolution of the BH angular momentum we cap aBH at the theoretical limit (in the presence of accreting matter) of 0.998 (Thorne 1974). Since the total angular momentum of the BH-disc system must be conserved, the disc angular momentum evolves aṡ whereJin is the angular momentum change due to the material flowing on to the accretion disc from resolved scales. As pointed out by Dotti et al. (2013), in some conditions Rwarp can exceed the total extent of the accretion disc Rout, i.e. the BH mass is larger than a critical value In this regime, the time-scale for alignment or counteralignment of the BH and disc angular momenta is drastically reduced and the disc cannot attain a steady warped state. Under these circumstances, we assume that JBH and J disc instantaneously reorient themselves towards the axis of the system's total angular momentum, Jtot. While JBH always tends to align itself to Jtot, J disc will end up aligned with Jtot only when (King et al. 2005) counter-aligning otherwise.
In order to accurately evolve the spin, we impose a time-step criterion to the sub-grid model, defined as a fraction (fixed at 10 per cent; as long as we properly resolve the investigated time-scales, the presented results are not significantly sensible to the choice of this fraction) of the minimum between the disc consumption time-scale, τ drain = M disc /Ṁacc−BH, and the gravito-magnetic interaction timescale τgm. Moreover, one of the underlying assumptions of our model is that the disc attains a steady-state warped profile, therefore introducing a further limitation on ∆t, which must be thus greater than the warp propagation time-scale τν 2 (Rwarp) (Martin et al. 2007):

Connecting the sub-grid model to simulations
In this section, we describe how we consistently couple our sub-grid model to the hydrodynamics code used for simulations. To this aim, at each time-step, we determine the boundary conditions for the sub-grid evolution from the average properties of the gas surrounding each BH particle. First, we determine the inflow rateṀin on to the disc as the accretion rate from resolved scales provided by the code.
To prevent the sub-grid disc from becoming selfgravitating, at every time-step we limitṀin to ensure M disc ≤ Msg, where Msg is computed as the disc mass enclosed within the self-gravitating radius R = Rsg, defined as the radius at which the Toomre parameter (Toomre 1964) Q is equal to unity: Every time the disc is depleted, we refill it in a selfconsistent way by taking the inflow rate from larger scales. The mass of the newly formed disc is instantaneously set to where M disc, seed is a user-defined seed-mass parameter for the disc, that we set to 10 5 M . Since the mass inflow from resolved scales could in principle be arbitrarily small, wheṅ Min∆t < M disc we refill the disc stochastically, with a probability q = ∆tṀin/M disc . More specifically, we randomly sample n ∈ (0, 1) from a uniform distribution: if n ≤ q, the disc is created with mass M disc ; otherwise, its mass is left to zero. We recreate the disc by choosing its angular momentum direction to be along jgas, and by setting f Edd = f Edd,0 , where f Edd,0 is a free parameter of the model.
We compute the angular momentum inflow rate on to the disc aṡ where Λin is the angular momentum per unit mass carried by the inflowing material. Unfortunately, when the angular momentum transfer is not properly resolved, Λin is too large to be supported by a self-gravitating disc, hence the gas inflowing from resolved scales cannot circularise and join the accretion disc. In these cases, we reduce Λin to account for any mechanism that would make the gas lose its angular momentum flowing from large scales down to the disc scale. Assuming that the gas circularises at a characteristic radius Rcirc, we limit Λin to J disc (Rcirc) /M disc (Rcirc), i.e. to the disc specific angular momentum per unit mass at Rcirc. In our model, Rcirc is specified as a parameter in units of the disc self-gravitating radius Rsg. Our treatment is different from that presented in Fiacconi et al. (2018), where they imposeṀin = 0 whenever Λin is too large. Their approach could excessively limit the inflow in simulations where the angular momentum transfer is not properly resolved. Moreover, at every time-step, we check if J disc /M disc > ΛISCO. If this condition is not satisfied, then the gas will not be able to settle into circular orbits and will fall on to the BH over time-scales much shorter than τν 2 (Rwarp). In these cases, we instantaneously add the disc mass and angular momentum to the BH and, immediately after, refill the disc as described above.

SEMI-ANALYTIC VALIDATION
To validate our prescriptions, we performed a first test of the BH evolution in a semi-analytically modelled environment, as in Dotti et al. (2013), wherein an initially non-rotating BH with MBH = 10 4 M grows in mass through subsequent accretion episodes. At the onset of each episode, we create a new accretion disc with mass M disc = min{Msg , M cloud }, where M cloud = 10 5 M is the maximum mass available for each episode. This procedure prevents the disc from becoming too massive for relatively large BH masses, since the cold gas fraction relative to stars is observed to decrease with increasing galaxy mass (e.g. di Serego Alighieri et al. 2007;Catinella et al. 2010). In order to compute the angular momentum of a newly formed disc, we must specify a value for f Edd , that determines the magnitude in Equation (5), and extract the direction by a Monte Carlo sampling of the distribution of misalignment angles θ between the angular momenta of the disc and the larger-scale gas reservoir. We control the degree of anisotropy in the fueling process by introducing a parameter F representing the fraction of discs forming with θ > π/2. We set the Eddington ratio to a fiducial average value f Edd = 0.1. However, we note that varying the choice for f Edd in [0.01, 1] does not qualitatively change the long-term BH spin evolution investigated with this model. For this validation test, we also assume that neither accretion nor outflows occur, i.e.Ṁin =Ṁout = 0.
Our results qualitatively reproduce those obtained by Dotti et al. (2013), although they exhibit a swifter alignment between the angular momenta, due to the stronger contribution of strongly misaligned configurations in our prescription for the gravito-magnetic torque (see Equation 8), whereas the model adopted by Dotti et al. (2013) holds only in the approximation of small misalignments (Perego et al. 2009).
In Figure 1, we present the evolution of the BH spin parameter aBH through several orders of magnitude in the growth of MBH (10 4 -10 9 M ). Independent of the initial BH-disc configuration, for relatively low BH masses, aBH rapidly grows up to its maximum value, because the angular momenta of the BH and disc align themselves on time-scales way shorter than the disc consumption time-scale. When Msg ≥ M cloud , the disc is refilled with the same mass at each episode, while J disc /JBH keeps decreasing as MBH grows. Choosing a different value for M cloud would only directly affect the value of MBH at which we observe this transition. In this regime, alignment becomes inefficient and the BH can spend more time growing in mass via retrograde accretion, thus significantly reducing aBH.
For large BH masses, JBH aligns itself with the average angular momentum of the gas reservoir, therefore aBH evolves towards an equilibrium value aeq set by the degree of anisotropy in the fueling process (Sesana et al. 2014). With respect to the prescription adopted by Dotti et al. (2013), our model predicts a swifter alignment with the reservoir, making prograde accretion events more frequent. Therefore, for large BH masses, aBH is biased toward higher values, the disc is recreated with smaller J disc /JBH at the onset of each accretion event, and we earlier enter the regime where Rwarp Rout. This results in values for aBH at large BH masses that are in better agreement with those computed analytically in the limit for J disc /JBH → 0.

NUMERICAL SIMULATIONS
In this work, we implement a physically motivated BH spin evolution sub-grid model in the publicly available N -body, mesh-less hydrodynamics code gizmo (Hopkins 2015), descendant of gadget2 (Springel 2005) and gad-get3 (Springel et al. 2008), although the model is flexible enough that it can be easily transported to other codes. In detail, our prescriptions follow the evolution of a BH-disc system associated to BH/sink particles in the code. Simulations were run on the CINECA cluster MARCONI100, with an typical usage of ∼ 10 CPU hours per Myr.

Numerical setup
We carried out simulations with a single BH particle in an idealised environment resembling a typical galactic nucleus, consisting of a gaseous circumnuclear disc (CND) embedded in a stellar bulge (Lupi et al. 2015). The spherical stellar component of our initial conditions follows a Hernquist (1990) density profile, where r is the radial spherical coordinate, M b = 5 × 10 8 M is the bulge total mass, and r b = 100 pc is the bulge scaleradius. The gas particles constitute a rotationally supported exponential disc in vertical hydrostatic equilibrium, with a surface density profile where MCND = 10 8 M is the disc total mass and RCND = 50 pc is the disc scale-radius. The vertical density profile and velocity field of the disc are calculated by means of the publicly available 6 code gd_basic (Lupi et al. 2015) taking into account the global potential of the bulge+disc+BH system. We initialised the gas component assuming an ideal equation of state with γ = 5/3 and uniform temperature T = 10 4 K. In order to reduce the pressure support of the disc and favour gas inflows towards the centre, the simulations are performed with a lower γ = 7/5, allowing us to mimic a mild radiative cooling without actually employing a dedicated sub-grid model ). We consider two initial BH masses: 10 7 and 5 × 10 7 M . The mass resolution is 10 3 M for both star and gas particles, translating into N b = 5 × 10 5 stellar particles and NCND = 10 5 gas particles. The spatial resolution is determined by the Plummerequivalent gravitational softening parameter ε, that is fixed at 0.16 and 1 pc for stellar and BH particles, respectively. For Table 1. Summary of the parameters assumed for our simulations. The subscript zero refers to quantities evaluated at initialization, i.e. at time t = 0. The runs' labels are chosen to recall the single parameter changed with respect to the Fiducial run: Rc stands for R circ ; Jd recalls a change in the initial disc angular momentum through a different choice for f Edd, 0 ; aBH stands for a BH, 0 ; Md and MBH stand for the initial disc and BH masses, respectively. The suffixes UL, VL, L, H, and VH refer to an initial configuration where a specific parameter is set to ultra-low, very-low, low, high, and very-high, respectively, relative to the Fiducial run. In run Fiducial+Feedback, we couple our model with the biconical-outflow model presented in Sala et al. (2021). gas particles/cells, we employ instead fully adaptive softening, i.e. the gravitational and hydrodynamic resolutions are both defined by the kernel size of each gas element, set to encompass an effective number of neighbours N ngb = 32. The minimum gravitational softening/kernel size, that also sets the maximum spatial resolution of the simulation for gas, is set to εgas = 0.16 pc. In order to couple our sub-grid model with the resolved scales in the hydrodynamics code, we determine the accretion rate on to the BH particle by means of the Bondi-Hoyle-Lyttleton (Bondi & Hoyle 1944;Bondi 1952;Hoyle & Lyttleton 1939) formula, 7 as implemented by Springel et al. (2005): where ρ is the density of the surrounding gas and v is the gas-BH relative velocity, both determined via mass-weighting over the nearest N ngb,BH neighbour particles of the BH. In our implementation, we set N ngb,BH = 3N ngb . Finally, αacc is a dimensionless parameter, typically employed to correct for the interstellar medium's dense gas not properly resolved in ∼kpc-scale simulations (Di Booth & Schaye 2009), which we set equal to one. In order to prevent the accretion properties from being un-physically affected by gas particles at large distances, we set a maximum accretion radius of 10 pc. Therefore, the BH particle kernel is defined as the region encompassing N ngb,BH neighbours, unless limited in size by the specified maximum accretion radius.
To compute the angular momentum inflow rate on to the disc,Jin, we assume that the angular momentum per unit mass carried by the inflowing material, Λin, is equal 7 We note, however, that any other prescription is equally valid.
to Jgas/Mgas, where Mgas is the total mass of the gas particles in the BH particle kernel, and Jgas is their total angular momentum.The direction of Jgas is also used to initialise j disc . We performed simulations initialising the BHdisc angular momenta misalignment angle to θ BH−disc = arccos (jBH · j disc ) = 5π/6, and varying the initial BH and disc masses and angular momentum ratio, as well as the parameter Rcirc/Rsg. As fiducial parameters for our model, adopted to initialise our Fiducial run, we took: MBH,0 = 10 7 M ; M disc,0 = 5 × 10 4 M ; f Edd,0 = 5 × 10 −3 ; 8 Rcirc/Rsg = 0.5; and aBH,0 = 0.5. Details on the choice of parameters for each performed simulation are reported in Table 1.

Results
In Figure 2, we show the effect of varying the parameter Rcirc/Rsg, which controls the angular momentum inflow on to the disc, on the misalignment angle θBH−gas = arccos (jBH · jgas) between the BH spin and the average angular momentum of the resolved gas reservoir (top panel), and on the Eddington ratio f Edd (bottom panel). We compare runs Rc-VL, Rc-L, Fiducial, Rc-H, and Rc-VH, for which we assume the same BH-disc initial configuration (see Table 1 for details). A smaller circularisation radius corresponds to a smaller upper limit for Λin, in turn related to a more compact disc (at fixed M disc ) and therefore to a larger accretion rate on to the BH (i.e. a larger f Edd ). This results in a swifter BH spin evolution for smaller values of Rcirc/Rsg and the difference in the evolution is more evident the smaller Rcirc/Rsg is.
In Figure 3, we show the impact of changing the initial 0 π/6 π/3 π/2 2π/3 5π/6 π θ BH−gas [rad] Rc-VL (R circ /R sg = 0.1) Rc-VH (R circ /R sg = 0.9) These runs share the same set of initial parameters, but for R circ /Rsg. A smaller R circ implies a smaller angular momentum inflow, resulting in a more compact disc with higher accretion rate and in a faster BH spin evolution.
Configurations with larger f Edd,0 correspond to lower initial J disc /JBH, resulting in more compact and denser discs with higher accretion rates, and leading to a faster evolution of the BH spin.
In runs Md-L and MBH-H, we changed the initial disc and BH mass, respectively, with respect to our Fiducial run. In run Md-L, we assume an initially less massive disc, with M disc, 0 = 10 4 M . In run MBH-H, we chose a heavier BH seed, with MBH, 0 = 5×10 7 M . The different time evolution of θ BH−disc , θBH−gas, and θ disc−gas = arccos (j disc · jgas), i.e. the misalignment angles between the angular momenta of the BH, disc, and gas reservoir in these runs, is presented in the lower panels of Figure 4. The upper panels show the time evolution of Jgas/Mgas (blue line) and Λin (black solid line), until J disc /M disc (magenta line) in run MBH-H drops below ΛISCO (dashed black line). All these runs share the same initial accretion rate f Edd, 0 = 5 × 10 −3 . This translates in the 10 4 M disc of run Md-L being initially more compact (less angular momentum, with equal f Edd ) with respect to runs Fiducial and MBH-H, hence in a faster BH evolution. In run MBH-H, the circularisation radius is closer to Rwarp These runs share the same set of initial parameters, but for f Edd, 0 . The smaller the initial Eddington ratio f Edd, 0 (i.e. larger initial total angular momentum of the disc), the slower is the BH spin evolution, because of a reduced accretion rate on to the BH from the less compact disc. than in the other runs, because of the more massive BH. The major contribution to the disc total angular momentum comes from regions closer to the BH, where we expect the disc to significantly modify its angular momentum direction. Both J disc and JBH tend to align/counter-align with Jtot. Therefore, a lower J disc /JBH, i.e. a reduced contribution of J disc to Jtot, results in a more evident evolution of J disc .
In run Md-L, In the Fiducial run, this effect is negligible, since the disc is more extended and less compact, whereas in run MBH-H, the evolution of J disc is much more pronounced. Moreover, J disc /JBH is small enough to expect the BH-disc system to evolve towards a configuration with counter-aligned angular momenta. Indeed, the angle θ BH−disc (black line in the bottom panels of Figure 4) eventually settles to ∼ π, whereas the angle θBH−gas (blue line in the bottom panel) approaches zero. Before J disc /M disc ΛISCO, the BH and disc angular momenta counter-align, with the BH spin close to being aligned with the average angular momentum direction of the reservoir, jgas. Therefore, J disc is almost counter-aligned with respect to jgas, and the coherent inflow of angular momentum directed as jgas acts to reduce J disc /M disc , which rapidly drops below the threshold of ΛISCO. In the evolution of J disc , our model assumes that the inflowing counter-rotating gas and the disc shock  With respect to our Fiducial run, run Md-L is initialised with a less massive disc (M disc, 0 = 10 4 M ), whereas run MBH-H presents a BH with a larger initial mass (M BH, 0 = 5 × 10 7 M ). In the top panels, we compare the evolution of J disc /M disc (magenta solid lines), Λ in (black solid lines), and Jgas/Mgas (blue solid lines). The black dashed lines represents Λ ISCO . In run Md-L, Λ ISCO exhibits a clear jump at ∼0.5 Myr, when the BH and disc switch from co-rotating to counter-rotating, because the ISCO radius differs in the two scenarios. This behaviour of Λ ISCO is seen in every run in which there is a switch between counter-rotating and co-rotating configurations (e.g. in the Fiducial run, but at later times than those presented in this figure). In the bottom panels, we show the evolution of θ BH−disc (black line), θ BH−gas (blue line), and θ disc−gas (magenta line). The more compact disc of run Md-L results in the BH and disc angular momenta ending up aligned over ∼1 Myr. In run MBH-H, after an initial slow evolution towards alignment, the BH and disc angular momenta rapidly counter-align.
instantaneously, settling to a more compact configuration with less angular momentum. This behaviour is not necessarily physical, but it is the best that can be done with a model that does not fully evolve self-consistently the whole angular momentum profile of the accretion disc. Such an improvement would significantly slow down any large scale simulation, and therefore is not considered in this study.
In Figure 5, we present the effect of changing the initial BH spin parameter aBH, 0 on the evolution of θBH−gas, by comparing runs aBH-L, Fiducial, and aBH-H, where we set aBH, 0 equal to 0.1, 0.5, and 0.8, respectively. A more rapidly spinning BH offers more resistance to changes in its angular momentum direction, resulting in a slower spin evolution (see top panel of Figure 5) and, after a transition period of a few Myr, higher accretion rates (bottom panel). In run aBH-H, accretion switches from retrograde to prograde at later times, allowing for f Edd to peak at larger values, up to an order of magnitude with respect to run aBH-L.
Finally, in our Fiducial+Feedback run, we coupled the spin evolution model to the new accretion/feedback model by Sala et al. (2021), assuming that biconical outflows are launched along the BH spin direction, and the outflow rateṀ out is self-consistently determined by the actual BH accretion rate provided by the accretion disc. The evolution of θBH−gas with time for this last run is shown in the top panel of Figure 5 as a blue dashed line, to facilitate the comparison with the no-feedback sibling run (Fiducial). While the initial evolution (up to ∼ 1 Myr) is similar, with a slightly slower alignment in the simulation with feedback due to the reduced mass accretion rate (accommodating for the generated outflow), the late evolution differs significantly. This is due to the increased relevance of the biconical feedback when the spin direction approaches the large-scale gaseous disc mid-plane and the feedback therefore impinges on to the dense gas distribution, strongly affecting its dynamics and further reducing the inflow rate on to the BH-disc system. These effects can be clearly seen in the bottom panel of Figure 5, where the accretion rate in Fiducial+Feedback is only a factor of two lower than in Fiducial during the first ∼Myr of the run, but then the difference increases up to a factor of ∼4-5.  Figure 5. Time evolution of the misalignment angle θ BH−gas (top panel) and the Eddington ratio f Edd (bottom panel), in runs aBH-L, aBH-H, Fiducial, and Fiducial+Feedback. These runs share the same set of initial parameters, but for a BH, 0 . Runs Fiducial+Feedback and Fiducial share the same a BH, 0 , but in the former we couple our model with the biconical-outflow model presented in Sala et al. (2021). A higher initial a BH results in a delayed alignment between the BH and disc angular momenta, and in a higher accretion rate, after a transition period of a few Myr.

CONCLUSIONS
We have introduced a novel sub-grid model for the evolution of the BH mass and spin and the surrounding accretion disc. Our model is a direct descendant of the model presented in Fiacconi et al. (2018), with some modifications concerning the limitation on the angular momentum inflow through Rcirc and the stochastic refill of the disc mass in the case of small mass inflow rates, as discussed in Section 2.2. Moreover, we do not impose any lower limit to f Edd , because the BH spin evolution can be neglected for very small accretion rates. In addition, the newly implemented model for the unresolved accretion disc accounts for the presence of self-consistently computed outflows (see Sala et al. 2021, for a detailed description of the implementation).
In the future, we are planning to further improve our model, in order to increase its accuracy and account for additional processes that can affect the BH spin. In particular, we will (i) employ a more sophisticated viscosity prescription, based on the numerical solutions for different warps (see, e.g. Ogilvie & Latter 2013;Tremaine & Davis 2014), (ii) account for the additional angular momentum coupling in BH binaries in gaseous environments (e.g. Gerosa et al. 2020) and the spin change after BH mergers, and (iii) consider a more appropriate coupling between resolved and sub-grid scales when the accretion disc is counter-aligned relative to the inflowing gas, that would relax the instantaneous shock assumption employed in this study.
We first assessed the validity of our model in a semianalytically modelled environment, qualitatively reproducing the results of Dotti et al. (2013). A crucial quantitative difference, however, is represented by the swifter alignment between the angular momenta in our model, that results from its validity also for strongly misaligned configurations.
We then interfaced our model with the publicly available code gizmo. We itemize the findings of our N -body, hydrodynamic tests below: • The radius at which the gas circularises, Rcirc, has profound effects on the BH evolution, with a swifter spin evolution and a higher accretion rate for smaller Rcirc.
• The initial total angular momentum of the disc, J disc , also plays an important role in the BH spin evolution, that becomes slower for larger values of J disc (or, equivalently, for smaller values of the initial Eddington ratio).
• Decreasing the initial accretion disc mass produces a disc with less angular momentum. The larger growth rate for M disc than for J disc translates into a larger f Edd at later times and, therefore, in a faster BH evolution.
• An initial larger (smaller) BH spin results in a slower (faster) evolution, due to increased (reduced) resistance to changes in its own angular momentum.
• The inclusion of BH feedback has the effect of altering the dynamics of the CND and reducing the gas inflow, in particular when the spin evolution leads to the feedback cone to cross the large-scale gas distribution. This, as a consequence, further slows down the evolution of the spin direction, when its relative angle respect to the CND angular momentum is about π/2.
In conclusion, our model is able to accurately follow the BH mass and spin evolution in hydrodynamic simulations, also when coupled with a sub-grid prescription for BH accretion and feedback, and can be easily applied to simulations on different scales, from galaxy mergers to cosmological simulations. In future works, we will employ it to investigate the evolution of BH spin in a cosmological environment, enabling us to make prediction for the Laser Interferometer Space Antenna (Amaro-Seoane et al. 2017; Barack et al. 2019) and pulsar timing arrays.

DATA AVAILABILITY STATEMENT
The data underlying this article will be shared on reasonable request to the corresponding author.