# Cosmic evolution of black hole spin and galaxy orientations: Clues from the NewHorizon and Galactica simulations

Sébastien Peirani<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> D, Yasushi Suto<fnref n="3" />,<fnref n="4" />,<fnref n="5" /> D, Ricarda S. Beckmann<fnref n="6" /> D, Marta Volonteri<fnref n="2" /> D, Yen-Ting Lin<fnref n="7" /> D, Yohan Dubois<fnref n="2" />, Sukyoung K. Yi<fnref n="8" />, Christophe Pichon<fnref n="2" />,<fnref n="9" />,<fnref n="10" />, Katarina Kraljic<fnref n="11" />, Minjung Park<fnref n="12" />, Julien Devriendt<fnref n="13" /> D, San Han<fnref n="8" /> D, and Wei-Huai Chen<fnref n="7" />,<fnref n="14" /> D

Received 25 December 2023 / Accepted 12 March 2024

## ABSTRACT

Black holes (BHs) are ubiquitous components of the center of most galaxies. In addition to their mass, the BH spin, through its amplitude and orientation, is a key factor in the galaxy formation process, as it controls the radiative efficiency of the accretion disk and relativistic jets. Using the recent cosmological high-resolution zoom-in simulations, NEWHORIZON and GALACTICA, in which the evolution of the BH spin is followed on the fly, we have tracked the cosmic history of a hundred BHs with a mass greater than $2 \times 10^4 M_\odot$. For each of them, we have studied the variations of the three-dimensional angle ($\Psi$) subtended between the BH spins and the angular momentum vectors of their host galaxies (estimated from the stellar component). The analysis of the individual evolution of the most massive BHs suggests that they are generally passing by three different regimes. First, for a short period after their birth, low-mass BHs ($M_{BH} < 3 \times 10^4 M_\odot$) are rapidly spun up by gas accretion and their spin tends to be aligned with their host galaxy spin. Then follows a second phase in which the accretion of gas onto low-mass BHs ($M_{BH} \le 10^5 M_\odot$) is quite chaotic and inefficient, reflecting the complex and disturbed morphologies of forming proto-galaxies at high redshifts. The variations of $\Psi$ are rather erratic during this phase and are mainly driven by the rapid changes of the direction of the galaxy angular momentum. Then, in a third and long phase, BHs are generally well settled in the center of galaxies around which the gas accretion becomes much more coherent ($M_{BH} > 10^5 M_\odot$). In this case, the BH spins tend to be well aligned with the angular momentum of their host galaxy and this configuration is generally stable even though BH merger episodes can temporally induce misalignment. We even find a few cases of BH-galaxy spin anti-alignment that lasts for a long time in which the gas component is counter-rotating with respect to the stellar component. We have also derived the distributions of cos($\Psi$) at different redshifts and found that BHs and galaxy spins are generally aligned. Our analysis suggests that the fraction of BH-galaxy pairs with low $\Psi$ values reaches maximum at $z \sim 4-3$, and then decreases until $z \sim 1.5$ due to the high BH-merger rate. Afterward, it remains almost constant probably due to the fact that BH mergers becomes rare, except for a slight increase at late times. Finally, based on a Monte Carlo method, we also predict statistics for the 2D projected spin-orbit angles $\lambda$. In particular, the distribution of $\lambda$ traces the alignment tendency well in the three-dimensional analysis. Such predictions provide an interesting background for future observational analyses.

**Key words.** methods: numerical – galaxies: evolution – galaxies: general – galaxies: interactions – galaxies: jets

## 1. Introduction

It is now well established that the evolution of supermassive black holes (BHs) and their host galaxy is intimately related (e.g., Silk & Rees 1998; Granato et al. 2004). The existence of strong scaling relations between BH masses and

galaxy properties, such as the velocity dispersion of their host bulge (Magorrian et al. 1998; Merritt & Ferrarese 2001a,b; Ferrarese & Merritt 2000; Hu 2008; McConnell & Ma 2013; van den Bosch 2016; Batiste et al. 2017; Baldassare et al. 2020), the stellar mass of their host bulge (Marconi & Hunt 2003; Häring & Rix 2004; Saglia et al. 2016; Sahu et al. 2019;

<footnote n="1">¹ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice Cedex 4, France
e-mail: speirani@oca.eu</footnote>

<footnote n="2">² Institut d'Astrophysique de Paris, CNRS and Sorbonne Université, UMR 7095, 98 Bis Boulevard Arago, 75014 Paris, France</footnote>

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

<footnote n="4">⁴ Research Center for the Early Universe, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan</footnote>

<footnote n="5">⁵ Research Institute, Kochi University of Technology, Tosa Yamada, Kochi 782-8502, Japan</footnote>

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

<footnote n="7">⁷ Institute of Astronomy and Astrophysics, Academia Sinica, No. 1, Section 4, Roosevelt Road, Taipei 10617, Taiwan</footnote>

<footnote n="8">⁸ Department of Astronomy and Yonsei University Observatory, Yonsei University, Seoul 03722, Republic of Korea</footnote>

<footnote n="9">⁹ IPhT, DRF-INP, UMR 3680, CEA, L'Orme des Merisiers, Bât 774, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="10">¹⁰ Korea Institute for Advanced Study, 85 Hoegi-ro, Dongdaemun-gu, Seoul 02455, Republic of Korea</footnote>

<footnote n="11">¹¹ Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS, UMR 7550, 67000 Strasbourg, France</footnote>

<footnote n="12">¹² Center for Astrophysics | Harvard & Smithsonian, Cambridge, MA, USA</footnote>

<footnote n="13">¹³ Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="14">¹⁴ Bienen School of Music, Northwestern University, 70 Arts Cir Dr, Evanston, IL 60208, USA</footnote>
Zhao et al. 2021), and the total stellar mass (Cisternas et al. 2011; Simmons et al. 2011; Reines & Volonteri 2015; Davis et al. 2019; Sahu et al. 2019; Ding et al. 2020; Bennert et al. 2021) does indeed suggest that BHs coevolve with their host galaxy.

Also, due to matter accretion, central BHs in galaxies often exhibit powerful jets and winds, commonly referred to as active galactic nuclei (AGNs). The resulting AGN feedback can heat and/or expel the surrounding gas from the center of galaxies. According to theoretical works, this mechanism is supposed to play an important role in reshaping the gas, stellar, and dark matter distributions in the central and outer parts of their host (e.g., Peirani et al. 2008, 2017, 2019; Duffy et al. 2010; Martizzi et al. 2013; Dubois et al. 2016; Ardila et al. 2021, and references therein) and in quenching the star formation in massive galaxies, producing more realistic galaxy population that are broadly consistent with observations (e.g., Di Matteo et al. 2005; Springel et al. 2005; Schawinski et al. 2006; Croton et al. 2006; Sijacki et al. 2007; Booth & Schaye 2009; Dubois et al. 2012; Choi et al. 2015; Kurinchi-Vendhan et al. 2023).

Black holes are generally characterized by two classical properties<fnref n="1" />: their mass $M_{BH}$ and their angular momentum vector $J_{BH}$. In addition, it is customary to introduce the dimensionless spin parameter $a$ assuming a Kerr metric for rotating bodies:

$$a \equiv \frac{cJ_{BH}}{GM_{BH}^2}, \qquad (1)$$

where $c$ is the speed of light and $G$ is the gravitational constant. Although one naturally expects that $a$ ranges from 0 to 1, it is not rare to find in the literature $a$ ranges from -1 to 1 which gives additional information on how the surrounding gas is accreted onto the BH. If the gas accretion disk does indeed have a prograde rotation, $a > 0$, while $a < 0$ for retrograde rotation.

From their initial values, both the BH the mass and spin evolve according to successive gas accretion phases and BH merger episodes. While the BH mass increases through cosmic time (though a small amount of mass is lost in the emission of gravitational waves during BH-BH mergers), the spin magnitude may decrease via different mechanisms such as gas accretion (Bardeen 1970), BH-mergers (Rezzolla et al. 2008; Barausse & Rezzolla 2009) and energy extraction by feedback (Blandford & Znajek 1977). However, the long-term evolution of the BH mass and spin is not independent but generally connected.

The amplitude and orientation of the BH spin control the efficiency of the conversion of the accreted gas into energy as

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

where $L_{BH}$ is the BH accretion luminosity, $\dot{M}_{BH}$ the BH mass accretion rate, $c$ then process. In addition, the spin of merging BHs (especially their relative orientation) has a strong influence on the gravitational wave emission and therefore on the expected recoil velocity of the merger remnant (Campanelli et al. 2007; González et al. 2007; Lousto & Zlochower 2011, 2013). Thus, modeling and characterizing the long-term evolution of the BH spin and its impact on the host galaxies are a valuable effort for many domains in modern astrophysics.

For this purpose, semi-analytical models and, above all, cosmological simulations are vital tools in characterizing BH mass growth and spin evolution. The task has however proven to be

difficult. The different physical processes that govern the BH evolution, in particular the gas accretion and the associated feedback processes, are still poorly understood and occur at physical scales (<pc) much lower than the resolution limit of the simulations. Therefore, one needs to employ empirical prescriptions and sub-grid physics.

While BHs are now commonly incorporated in hydrodynamical simulations generally as sink particles of a given mass, their spin is often neglected (e.g., Di Matteo et al. 2008; Booth & Schaye 2009; Dubois et al. 2014a; Vogelsberger et al. 2014; Schaye et al. 2015; Kannan et al. 2022). Nevertheless, increasing effort in the last decade has tried to fill the gap by proposing improved modeling including spin evolution in traditional cosmological simulations (e.g., Dotti et al. 2013; Dubois et al. 2014b, 2021; Fiacconi et al. 2018; Bustamante & Springel 2019; Talbot et al. 2021; Hopkins et al. 2023; Dong-Páez et al. 2023; Huško et al. 2024; Rennehan et al. 2023) as well as in general relativistic magneto-hydrodynamical simulations (e.g., Fedrigo et al. 2024; Cui et al. 2023; Koudmani et al. 2023, and references therein).

The bulk of the past effort on the cosmic evolution of black hole spin has focused almost exclusively on its amplitude and on its relative alignment with the gas disk using semi-analytical models (e.g., Volonteri et al. 2005, 2013; Shapiro 2005; King et al. 2008; Berti & Volonteri 2008; Fanidakis et al. 2011; Dotti et al. 2013; Sesana et al. 2014; Griffin et al. 2019; Izquierdo-Villalba et al. 2020) or cosmological simulations (e.g., Dubois et al. 2014b,c; Bustamante & Springel 2019; Beckmann et al. 2024; Sala et al. 2024). To date, there are only a few attempts that compute the relative orientation of the BH spin orientation with respect to that of their host galaxy.

A first detailed analysis of the evolution of the BH spin and orientation from numerical study in a cosmological context was presented in Dotti et al. (2013). They found the occurrence of two regimes. An early phase ($M_{BH} \le 10^7 M_\odot$) in which rapid alignment of the BH spin direction to the gas accretion disk angular momentum in each single episode leads to erratic changes in the BH spin orientation. For more massive BHs ($>10^7 M_\odot$), a single accretion episode does not modify significantly the BH spin direction, and the BH spin tends to align with the direction of the angular momentum of the accreting material.

Also, a pioneering work based on cosmological simulation was done by Dubois et al. (2014b). They found that BHs with mass $10^7 \le M_{BH} \le 10^8 M_\odot$ show a high level of alignment with their host galaxy, as gas accretion is mostly responsible for both the BH and stellar mass growth in this mass range. As BHs become more massive, spins are more randomly oriented with respect to their host galaxy angular momentum (see their Fig. 10). Moreover, mergers rapidly change the orientation of galaxies, in addition to coalescence for BHs, but the latter does not necessarily follow the same evolution while they merge. They also highlighted an early phase of misalignment when gas was turbulent and the galaxy disk not have coherent angular momentum structures. More recently, Beckmann et al. (2024) used the HORIZON-AGN simulation from which black hole spin evolution was post-processed following Dubois et al. (2014b). One of the main results is that merger-free galaxies tend to have higher BH spins which are preferentially aligned with their host galaxy's spin.

The aim of the present paper is to contribute to such theoretical efforts by presenting the most up-to-date detailed statistical study of the evolution of the BH spins relative to the galaxy angular momentum vectors. This is a complementary paper to Beckmann et al. (in prep.) which focuses on the

<footnote n="1">¹ A third one, the electric charge, is assumed to be negligible in most astrophysical settings (Gürlebeck 2015).</footnote>
evolution of the amplitude of the BH spin. Both analyses rely on the NEWHORIZON simulation (Dubois et al. 2021) that presents several advantages over HORIZON-AGN. First, it includes a sophisticated prescription of the BH model, wherein the BH spin is evolved for all BHs using accretion disk models, and the jets are launched along the direction of the BH spin vectors, with spin-dependent efficiencies. Second, it follows the evolution of a statistical number of galaxies and black holes. Third, it has a high spatial resolution (~34 pc) that accurately resolves a typical scale height of galactic disks (see, for instance, Park et al. 2021), and therefore allows for a reliable estimation of their stellar angular momentum. NEWHORIZON has also a sufficient resolution to capture the injection scale of gas turbulence, and has a multiphase interstellar medium. Furthermore, existing observations of outflows accelerated by active galactic nuclei suggest that some radio jets are inclined with respect to the galaxy disk (e.g., Morganti et al. 2015; Venturi et al. 2021), which strongly suggests misalignment between the central BH and the galaxy angular momentum. However, a recent observational analysis of a sample of 3682 radio AGNs (with reliable radio and optical position angle measurements) suggests a tendency of BH/galaxy spin alignment, especially for lower radio luminosity (Zheng et al. 2024). These observational trends give therefore additional motivation for this work.

It is interesting to note the similarity of the spin-orbit orientation between the BH-galaxy and star-planet systems. Planets are occasionally found to have misaligned or even retrograde orbits with respect to the direction of stellar rotation. Specifically, Kamiaka et al. (2019) found that a non-negligible fraction (~20%) of hot Jupiters exhibits the projected spin-orbit angle $\lambda > 30^\circ$ (see also, e.g., Ohta et al. 2005; Winn & Fabrycky 2015; Albrecht et al. 2022). The origin of those misalignment between the stellar spin and planetary orbits is not yet well-understood. Takaishi et al. (2020) performed a series of hydrodynamical simulations of the collapse of turbulent molecular cloud cores, and suggested the presence of two regimes: at the initial epoch when the protostar formed, $\Psi$ (i.e., the angle between the protostellar spin and the protoplanetary disk rotation axes) is very broadly distributed within ~130°. With the subsequent mass accretion from the proto-planetary disk to the protostar, however, the part of the angular momentum of the disk is transferred to the stellar spin, and $\Psi$ gradually decreases and tends to be aligned (≤20°).

Thus, they concluded that the isolated planetary disks formed in the turbulent cloud cores are unlikely to explain the observed misalignment, implying the importance of the interaction with the nearby stellar systems, or the gravitational planet-planet scattering after the gas disk dispersal.

Those results inspired us to ask the same question for the orientation between the central BH spin and the host galaxy angular momenta. Disk galaxies and planetary systems exhibit several similarities: they are both dominated by a central object (BH/star) and surrounded by either a galactic disk of stars and gas or a planetary disk of accreting gas. One may naturally ask whether the same physical mechanisms (inherent or external) operate in the formation process of the two astrophysical objects. It is definitely interesting to explore this possibility by exploring the connection of BH-galactic disk and central star-planetary disk as well as their mutual coevolution, to constrain the formation of planetary systems and galaxies, two major topics of contemporary astrophysics.

The paper is organized as follows. Section 2 briefly introduces the NEWHORIZON simulation and the numerical modeling used in this work (simulations and post-processing). Section 3 presents our main results on the 3D statistics of the BH-galaxy

spin angles while in Sect. 4 we derive statistics on the 2D projected angles to make predictions for future observational studies. We summarize our results and conclusions in Sect. 5.

## 2. Simulation data

Throughout this paper, we analyze the results of the NEWHORIZON<fnref n="2" /> simulation. The details of the simulation has been described in detail in Dubois et al. (2021), so we only summarize here its main features. Those who are interested in the BH spin results may skip this section and move to Sect. 3.

### 2.1. General

NEWHORIZON is a high-resolution zoom-in simulation from the HORIZON-AGN simulation (Dubois et al. 2014a), which extracts a spherical sub-volume with a radius of 10 comoving Mpc. A standard $\Lambda$CDM cosmology was adopted with the total matter density $\Omega_m = 0.272$, the dark energy density $\Omega_\Lambda = 0.728$, the baryon density $\Omega_b = 0.045$, the Hubble constant $H_0 = 70.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$, the amplitude of the matter power spectrum $\sigma_8 = 0.81$ and the power-law index of the primordial power spectrum $n_s = 0.967$, according to the WMAP-7 data (Komatsu et al. 2011). The initial conditions have been generated with MPGRAFIC (Prunet et al. 2008) at the resolution of 4096³ for NEWHORIZON in contrast to 1024³ for HORIZON-AGN. The dark matter mass resolution reaches $1.2 \times 10^6 M_\odot$ compared to $8 \times 10^7 M_\odot$ in HORIZON-AGN.

Both simulations were run with the RAMSES code (Teyssier 2002) in which the gas component is evolved using a second-order Godunov scheme and the approximate Harten-Lax-Van Leer-Contact (HLLC, Toro 1999) Riemann solver with linear interpolation of the cell-centered quantities at cell interfaces using a minmod total variation diminishing scheme. In NEWHORIZON, refinement is performed according to a quasi-Lagrangian scheme with the highest resolution of $\Delta x = 34 \text{ pc}$ at $z = 0$. The refinement is triggered in a quasi-Lagrangian manner, if the number of DM particles becomes greater than 8, or the total baryonic mass reaches 8 times the initial DM mass resolution in a cell. Extra levels of refinement are successively added at $z = 9, 4, 1.5$ and $0.25$ (i.e., $a = 0.1, 0.2, 0.4$ and $0.8$ respectively). The NEWHORIZON simulation is currently completed at redshift $z = 0.18$. Additionally, we have analyzed two other zoom simulations (nicknamed GALACTICA) focusing on isolated galaxies. For them, we have used exactly the same physics and mass resolution as NEWHORIZON but they are located in different regions of HORIZON-AGN (see, for instance, Park et al. 2021).

### 2.2. Gas and stellar physics

The gas follows an equation of state for an ideal monoatomic gas with an adiabatic index of $\gamma_{ad} = 5/3$. Gas cooling is modeled assuming equilibrium chemistry with rates tabulated by Sutherland & Dopita (1993) above $10^4 \text{ K}$ and by Rosen & Bregman (1995) below $10^4 \text{ K}$. Gas is also heated via a uniform ultraviolet radiation after the reionization epoch at $z = 10$ following Haardt & Madau (1996). Star formation (SF) is also included: stars can form out of gas cells with a hydrogen number density greater than $n_0 = 10 \text{ cm}^{-3}$ and a temperature lower than $2 \times 10^4 \text{ K}$, following a Schmidt relation:

$$d\rho_*/dt = \epsilon_* \rho_{\text{gas}}/t_{\text{ff}}, \quad (3)$$

<footnote n="2">² https://new.horizon-simulation.org/</footnote>
where $\rho_{\text{gas}}$ is the gas density, $t_{\text{ff}}$ is the free-fall time of the gas, and $\epsilon_{*}$ is the efficiency of star formation per free-fall time.

Contrary to HORIZON-AGN, NEWHORIZON adopts the efficiency of star formation that depends on the local turbulent Mach number and Jeans length (Kimm et al. 2017; Trebitsch et al. 2017, 2021). The initial mass function follows a Chabrier functional form (Chabrier 2005) with cutoffs at 0.1 and 150 $M_{\odot}$. Finally, a model of type II supernovae is based on the amount of linear momentum injected at the adiabatic or snow-plow phase (Kimm & Cen 2014; Kimm et al. 2015). The typical mass resolution for star particles in NEWHORIZON is ~$10^4 M_{\odot}$.

## 2.3. Black hole physics

Here we summarize the BH sub-grid physics which is implemented in NEWHORIZON, and particularly relevant to our current investigation of BH-galaxy spin relative orientation. It includes a model for BH mass growth and AGN feedback in alternating radio/quasar (jet/heating) mode (Dubois et al. 2012) coupled to a model of BH spin evolution (Dubois et al. 2014b).

### 2.3.1. Formation, mass growth and dynamics

First, supermassive BH seeds are allowed to form within any cells satisfying the following criteria:

(i) Both the stellar and gas densities exceed the threshold for star formation.

(ii) The local stellar velocity dispersion is larger than $20 \text{ km s}^{-1}$.

(iii) No preexisting BH can be found within a distance of 50 comoving kpc from the cell.

Then BHs form with an initial mass of $10^4 M_{\odot}$ and a spin parameter of $a = 0$. Then, their mass grows at a rate $\dot{M}_{\text{BH}}$ over time by accreting gas following an un-boosted Bondi-Hoyle-Lyttleton accretion rate $\dot{M}_{\text{Bondi}}$ and a spin-dependent radiative efficiency $\epsilon_r$:

$$ \dot{M}_{\text{BH}} = (1 - \epsilon_r) \dot{M}_{\text{Bondi}}, \quad (4) $$

$$ \epsilon_r = f_{\text{att}}(1 - e_{\text{isco}}) = f_{\text{att}} \left(1 - \sqrt{1 - 2/(3r_{\text{isco}})}\right), \quad (5) $$

$$ \dot{M}_{\text{Bondi}} = 4\pi\bar{\rho} \frac{(GM_{\text{BH}})^2}{(\bar{u}^2 + \bar{c}_s^2)^{3/2}}, \quad (6) $$

where $e_{\text{isco}}$ is the energy per unit rest mass energy of the innermost stable circular orbit (ISCO) of radius $r_{\text{isco}}$, $\bar{u}$ is the average BH-to-gas relative velocity, $\bar{c}_s$ the average gas sound speed and $\bar{\rho}$ the average gas density. Those quantities are calculated by averaging over a sphere of radius $4\Delta x$ (~150 pc) of the considered BH and using mass and kernel weighting (Dubois et al. 2012). $r_{\text{isco}}$ is in units of the half of the Schwarzschild radius and depends on the BH spin magnitude and its orientation with respect to the gas accretion disk angular momentum.

Finally, the Bondi-Hoyle-Lyttleton accretion rate is capped at the Eddington luminosity rate for the appropriate $\epsilon_r$:

$$ \dot{M}_{\text{Edd}} = \frac{4\pi G M_{\text{BH}} m_{\text{p}}}{\epsilon_r \sigma_{\text{T}} c}, \quad (7) $$

where $\sigma_{\text{T}}$ is the Thomson cross-section and $m_{\text{p}}$ the proton mass. According to Eq. (2), a fraction of the mass $\epsilon_r$ accreted is radiated away, while the rest of the mass $(1 - \epsilon_r)$ is accreted onto the BH and increases the BH mass. We introduce the Eddington ratio $\chi$:

$$ \chi = \begin{cases} \dot{M}_{\text{Bondi}} / \dot{M}_{\text{Edd}} & \text{if } \dot{M}_{\text{Bondi}} < \dot{M}_{\text{Edd}} \\ 1 & \text{if } \dot{M}_{\text{Bondi}} > \dot{M}_{\text{Edd}} \end{cases} \quad (8) $$

and distinguish the different AGN feedback at the threshold value of $\chi_{\text{trans}} = 0.01$ as discussed in Sect. 2.3.3. We note that in the case of the radio mode (see Sect. 2.3.3), $\epsilon_r$ used to estimate the effective growth of the BH is attenuated by a factor $f_{\text{att}} = \min(\chi/\chi_{\text{trans}}, 1)$ following Benson & Babul (2009). We also impose a maximum value of the BH spin at $a_{\text{max}} = 0.998$ due to the emitted photons by the accretion disk captured by the BH (Thorne 1974).

Furthermore, BHs merge when they get closer than $4\Delta x$ and the relative velocity of the pair is smaller than the escape velocity of the binary. The less massive BH of the binary is absorbed into the more massive one. It is then possible to know for each (existing) black hole, the percentage of its mass that has been gained through BH mergers at a specific redshift, $f_{\text{BH,merger}}$. $f_{\text{BH,merger}} = 0\%$ means no mass contribution from BH mergers while high values of $f_{\text{BH,merger}}$ indicate that a large percentage of mass has been gained through mergers over the BHs history. It is also worth mentioning that no recoil velocities have been applied to BHs during the merger process.

Finally, due to the finite force resolution effect, an explicit drag force is introduced for the gas onto the BH following Ostriker (1999), in order to avoid any spurious motion that can arise especially around high density gas regions. A detailed analysis of BH mergers in NEWHORIZON is presented in Volonteri et al. (2020), and an analysis of the population statistics of intermediate mass BHs in dwarf galaxies is found in Beckmann et al. (2023).

### 2.3.2. Spin evolution model

The BH spin is modeled on-the-fly in NEWHORIZON and updated according to the gas accretion and BH-BH mergers. This is a major improvement relative to HORIZON-AGN in which the evolution of the BH spin was not included and required some post-processing afterwards (Dubois et al. 2014b; Beckmann et al. 2024). The spin evolution follows specifically the model detailed in Dubois et al. (2014b,c) to which we refer the readers for a full description and technical details. The only difference is that a new model of spin evolution is adopted for low accretion rates $\chi < \chi_{\text{trans}}$. We summarize the main points here.

We stress that the accretion disk is not resolved in the simulation. Thus, at high accretion rates, a thin accretion disk solution by Shakura & Sunyaev (1973) is adopted. Then the evolution of the amplitude of the BH spin $a$ through gas accretion is estimated following Bardeen (1970):

$$ a_{n+1} = \frac{1}{3} r_{\text{isco}}^{1/2} \frac{M_{\text{BH},n}}{M_{\text{BH},n+1}} \left[ 4 - \left( 3r_{\text{isco}} \left( \frac{M_{\text{BH},n}}{M_{\text{BH},n+1}} \right)^2 - 2 \right)^{1/2} \right], \quad (9) $$

where $n$ refers to the value of the different variables at the $n$th timestep.

We note that Eq. (9) assumes that BH spin and accretion disk angular momentum are perfectly aligned or anti-aligned (i.e., only the BH spin amplitude changes but not its direction), which is not the case in general. A misaligned accretion disk experiences a torque due to the Lense-Thirring effect, and precesses the spin axis of the BH and warps the innermost parts of the disk (for large enough viscosity). The result of the Lense-Thirring effect is that the BH spin and the accretion disk angular momentum tend to align (or anti-align) with the total angular momentum. If one defines the total angular momentum of the system BH+disk at the $n$th timestep by $\mathcal{J}_{\text{tot},n} = \mathcal{J}_{\text{BH},n} + \mathcal{J}_{\text{d},n}$, the direction of the new BH spin is obtained by $\mathcal{J}_{\text{BH},n+1} = \mathcal{J}_{\text{tot},n} = \mathcal{J}_{\text{BH},n} + \mathcal{J}_{\text{d},n}$,
due to the conservation of total angular momentum. Its co- or counter-rotation with respect to the accretion disk is decided following criterion from King et al. (2005). In NEWHORIZON, the orientation of the BH spin is updated first, and then its amplitude (using Eq. (9)).

Finally, at lower accretion rates, $\chi < 0.01$, jets are assumed to be powered by energy extraction from black hole rotation (Blandford & Znajek 1977). Consequently, the black hole spin only decreases and the variations of $da/dt$ are obtained using the polynomial fits in McKinney et al. (2012). The orientation of the BH spin is updated following the same procedure as for high accretion rates.

The BH-BH mergers are also taken into account in the simulation. The spin of the remnant is calculated according to the spin of each BH prior the merger as well as the orbital angular momentum of the binary system. Specifically, we adopt the analytical fit of Rezzolla et al. (2008), assuming a random orientation between spins and the orbital angular momentum.

### 2.3.3. AGN feedback

The AGN feedback follows two different prescriptions, depending on the Eddington ratio $\chi$ defined in Eq. (8). Energy is released with an efficiency $\eta$ after each accretion episode in the form of

$$L_{\text{AGN,R,Q}} = \eta_{\text{R,Q}} \dot{M}_{\text{BH}} c^2, \quad (10)$$

where R and Q stand for the radio and quasar heating mode, respectively.

$-\chi < \chi_{\text{trans}}$ (“radio mode”) BHs power jets that continuously release mass, momentum and energy. Bipolar jets are assumed as a cylinder of size $\Delta x$ in radius and semi-height, centered on the BH (Dubois et al. 2010). The AGN axis is aligned to the spin of the BH direction (without any opening angle). The jets are launched with a speed of $10^4$ km s$^{-1}$. We note that $\eta_R$ is not a free parameter but computed from the BH spin following the results of McKinney et al. (2012; see also Dubois et al. 2021, for the interpolating function).

$-\chi > \chi_{\text{trans}}$ (“quasar mode”) BHs release only thermal energy into the gas (Teyssier et al. 2011) within a sphere of radius $\Delta x$ (isotropic and uniformly distributed). The efficiency of the feedback in quasar mode is given by $\eta_Q = \epsilon_c \eta_c$, where $\eta_c = 0.15$ is calibrated on the local $M_{\text{BH}}$-galaxy mass in lower resolution (~kpc) simulations (Dubois et al. 2012).

### 2.4. Galaxy and black hole catalogs

Dark matter halos and galaxies are identified using the ADAPTAHOP structure finder (Aubert et al. 2004; Tweed et al. 2009) at different redshifts using a local density threshold of 178 and 80 times the average DM and stellar densities, respectively. ADAPTAHOP allows to separate substructures from their host halos/galaxies. Since NEWHORIZON is a zoom simulation, low-mass resolution dark matter particles might “pollute” some halos, especially when they are located close to the boundary of the high resolution area. We remove those DM halos and their embedded galaxies in the following statistical analysis, except in Appendix A where we present individual evolution histories of several BHs in “contaminated” DM halos in which low resolution DM particles represent less than 0.1% of the total mass of the halos. Also, to give an order of magnitude, ~35% of galaxies with a stellar mass greater than $10^6 M_{\odot}$ at $z = 0.18$ lies in non contaminated halos. In the following, we define the galaxy mass by the value returned by ADAPTAHOP.

To link BHs to galaxies, we adopt the same methodology developed in previous studies using either HORIZON-AGN or NEWHORIZON (e.g., Volonteri et al. 2016; Smethurst et al. 2024; Beckmann et al. 2024). In the first step, we loop over all galaxies, from the most to the least massive ones, and identify and associate for each of them the most massive BH to be contained within 2 effective radii of the galaxy’s center. Such objects are then flagged as a primary BH, and removed from the list of non-allocated BHs. Then we repeat the second loop, and label all BHs within two effective radii as secondary BHs. All non-allocated BHs are finally removed from the sample since they are considered too far (“wandering”). In other words, a galaxy can contain multiple BHs, but any BH is associated uniquely to a single galaxy. In this scheme, we use a shrinking sphere approach (Power et al. 2003) to determine precisely the galaxy center. We also estimate the effective radius $R_e$ of each galaxy by taking the geometric mean of the half-mass radius of the projected stellar densities along each of the simulation’s Cartesian axes.

Initially, our sample extracted from NEWHORIZON contains 572 BHs (379 primary and 193 secondary BHs) at $z = 0.18$. This sample is completed with two additional primary BHs provided by the two GALACTICA zoom-in simulations. However, BHs with mass below $2 \times 10^4 M_{\odot}$ are discarded from our statistical analysis because they are too close to the initial seed BH mass and likely to suffer from mass resolution. Furthermore, we mainly focus on primary black holes. All of these constraints finally lead to a sample of 102 primary BHs at $z = 0.18$.

## 3. Three-dimensional statistics

### 3.1. General

In this section, we mainly focus on the cosmic evolution of the 3D angle (hereafter $\Psi$) between the BH spin vector, $J_{\text{BH}}$, and the angular momentum vector of the host galaxy stellar component, $J_{\text{star}}$, as well as its correlations to specific galaxy and BH properties. If not specified, the galaxy spin is estimated from the total angular momentum vector of all star particles associated with the galaxy within one effective radius $R_e$. The angle $\Psi_{\text{star,BH}}$ between the two vectors is then simply computed by

$$\cos(\Psi_{\text{star,BH}}) = \frac{J_{\text{star}} \cdot J_{\text{BH}}}{\|J_{\text{star}}\| \|J_{\text{BH}}\|}. \quad (11)$$

As such, $\Psi_{\text{star,BH}}$ can only lie between 0 and 180 degrees.

We similarly define the angles, $\Psi_{\text{gas,BH}}$ and $\Psi_{\text{star,gas}}$ between the angular momentum of the gas accretion disk ($J_{\text{gas}}$) and the BH spin and the stellar spin, respectively. We recall that the gas accretion disk properties are estimated within the four closest cells in radius from the BH position and using mass and kernel weighting. Figure 1 presents a schematic view of the different galaxy components as well as the relevant angles we are particularly interested in.

Let us begin the analysis by taking a look at the BH population properties as well as host galaxies at different redshift. In Fig. 2, we plot the primary BHs mass against their host galaxy mass at three different redshift ($z = 2, 1$ and $0.18$). We also use a color code to indicate the associated value of $\cos(\Psi_{\text{star,BH}})$. Since the volume covered by NEWHORIZON is relatively small, and due to our tight constraints on BH and galaxy selection, the total number of BHs in our analysis with $M_{\text{BH}} > 10^5 M_{\odot}$ is quite low: at $z = 2$, only 7 BHs has been identified, and this number goes up to 11 at lower redshifts. As already pointed out in Dubois et al. (2021), on average, central massive BHs in NEWHORIZON grow significantly only above a
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Schematic view of the different galaxy components and the relevant angles studied throughout the paper. We note that $\Psi_{\text{star,BH}}$ is the 3D angle between the BH spin and the angular momentum of the host galaxy. The latter was estimated by considering all star particles within a sphere of radius $R_e$. Furthermore, $\Psi_{\text{gas,BH}}$ is the 3D angle between the BH spin and the gas accretion disk gas ($J_{\text{gas}}$) estimated within the four closest cells (in radius) from the BH position and using mass and kernel weighting.</figcaption></figure>

stellar mass threshold of a few $10^{10} M_{\odot}$. For galaxies with stellar masses lower than $5 \times 10^9 M_{\odot}$, however, BHs growth is in general regulated and limited by supernovae feedback (see also Dubois et al. 2015; Habouzit et al. 2017; Trebitsch et al. 2017; Lapiner et al. 2021). Moreover, due to their low seed mass, BHs may not remain strongly attached to the centers of their host galaxies, especially in the dwarf galaxy regime as studied in Beckmann et al. (2023). Consequently, most of the BHs population grows little over the course of the simulation.

As far as the BH-galaxy spin orientation is concerned, two trends seem to emerge. First, at high redshift ($z = 2$), most of BHs tend to have their spin aligned with the host galaxy spin. Second, the most massive BHs ($>10^5 M_{\odot}$) at all redshifts also tend to have a spin aligned to $J_{\text{star}}$. In the next sections, we investigate the origin of these different trends.

Finally, since we estimate the angular momentum of each galaxy within one effective radius $R_e$, it is instructive to derive the variations of $R_e$ as a function of the total stellar mass. One major advantage of NEWHORIZON is its high resolution (~34 pc), enabling to describe even the size of low-mass galaxy ($10^6–10^7 M_{\odot}$) with several resolution elements. As reported previously in Dubois et al. (2021), NEWHORIZON reproduces fairly well the size-mass relation similar to observations at all redshifts (see their Fig. 12).

## 3.2. Cosmic evolution of $\Psi_{\text{star,BH}}$

### 3.2.1. Individual histories

In order to have a better view on the cosmic evolution of relevant BH properties and $\Psi_{\text{star,BH}}$, we explore further in Figs. 3–8 the

<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Variations of the primary BH mass with respect to their host (stellar) galaxy mass at three different redshifts. In each panel, the color code indicates the value of the corresponding angle between the BH spin and the angular momentum of the stellar component computed within one effective radius. The most massive BHs ($>10^5 M_{\odot}$) tend to have a spin well aligned to $J_{\text{star}}$ at all redshifts.</figcaption></figure>

evolution of four BHs that represent different histories. We also present in Appendix A the evolution of other BHs (see Table 1).

The first BH, BH-G9685, is extracted from a GALACTICA zoom. Its host galaxy is selected from the HORIZON-AGN volume, and exhibits a very quiet evolution, dominated by a smooth accretion and no major merger. Before $z = 1$, the galaxy has a clear disk with spiral arms but turns into an S0 galaxy at a lower redshift mainly due to the quenching of star formation. Its final stellar mass is $2.4 \times 10^{10} M_{\odot}$ at $z = 0$. The projected distribution
<table><caption>**Table 1.** Simulated BHs that are examined in detail throughout this paper.</caption><thead><tr><th>BH id</th><th>BH mass</th><th>z</th><th>Figures</th><th>Comments</th></tr></thead><tbody><tr><td>166</td><td>1.4 × 10<sup>8</sup> M<sub>☉</sub></td><td></td><td>5, 6</td><td>BH-mergers</td></tr><tr><td>455</td><td>6.9 × 10<sup>7</sup> M<sub>☉</sub></td><td></td><td>A.1</td><td>BH-mergers</td></tr><tr><td>1049</td><td>1.8 × 10<sup>7</sup> M<sub>☉</sub></td><td></td><td>7, 8, 10</td><td>Anti-aligned spin</td></tr><tr><td>936</td><td>3.3 × 10<sup>6</sup> M<sub>☉</sub></td><td></td><td>A.1</td><td>Smooth accretion</td></tr><tr><td>549</td><td>5.8 × 10<sup>5</sup> M<sub>☉</sub></td><td>0.18</td><td>9, A.1</td><td>Anti-aligned spin</td></tr><tr><td>132</td><td>1.5 × 10<sup>5</sup> M<sub>☉</sub></td><td></td><td>A.1</td><td>Off centered BH</td></tr><tr><td>146</td><td>7.8 × 10<sup>7</sup> M<sub>☉</sub></td><td></td><td>A.2</td><td>&lt;1% LR DM cont.</td></tr><tr><td>348</td><td>7.1 × 10<sup>7</sup> M<sub>☉</sub></td><td></td><td>A.2</td><td>&lt;1% LR DM cont.</td></tr><tr><td>541</td><td>3.5 × 10<sup>7</sup> M<sub>☉</sub></td><td></td><td>A.2, 10</td><td>&lt;1% LR DM cont.</td></tr><tr><td>G9685</td><td>4.8 × 10<sup>6</sup> M<sub>☉</sub></td><td>0</td><td>3, 4, 10</td><td>Smooth accretion</td></tr><tr><td>G648</td><td>1.6 × 10<sup>7</sup> M<sub>☉</sub></td><td>0.26</td><td>A.2</td><td>Smooth accretion</td></tr></tbody></table>

**Notes.** The upper part of the table indicates primary BHs extracted from NEWHORIZON at $z = 0.18$. The first five lines correspond to the most massive BHs of our fiducial sample while BH-549 and BH-132 present specific features namely an anti-aligned BH-galaxy spin configuration and a off-centered BH respectively. The last three lines gathers three other massive primary BHs but with a “contamination” of <0.1% in mass from lower resolution (LR) DM particles. For these latter, we just show their evolution in Appendix A while they are discarded from any statistical analysis. The lower part of the table is related to the additional two BHs from the GALACTICA zooms, G9685 and G648.

of the gas density and stars (using $u-g-r$ bands) at different redshifts is shown in Fig. 3. The associated BH is seeded at redshift $z = 7.5$ ($T_{\text{Universe}} = 0.7$ Gyr) and “survives” until $z = 0$ (i.e., it did not merge with a more massive BH). It reaches the final mass of $4.8 \times 10^6 M_{\odot}$. It exhibits very quiet evolution as well, as indicated in the upper panel of Fig. 4.

This BH underwent only two merger episodes with other BHs in its life; a major merger at $z \sim 4$ and a minor one at $z \sim 1$ (minor and major merger are delimited by a 1:4 mass ratio). The percentage of its mass gained through the BH mergers is $f_{\text{BH,merger}} = 1.53\%$ at $z = 0$. For a very short period right after its birth, the BH is rapidly spun up ($|a|$ reaches 0.865), driven by gas accretion. During this early phase, its spin tends to be aligned with its host ($\Psi_{\text{star,BH}} < 50^{\circ}$). Then, the evolution of the BH encompassed two long phases. The first one lasts until $T_{\text{Universe}} \sim 3.9$ Gyr ($z \sim 1.7$), and is characterized by accretion of gas that is rather chaotic, suggested by the evolution of the $\Psi_{\text{gas,BH}}$, and inefficient (i.e., low Eddington rates), mainly caused by the complex and disturbed morphology of the forming proto-galaxy at high redshifts (see for instance the upper panel of Fig. 3). Note also that the BH tends to be in radio mode and the amplitude of its spin generally decreases because the extracted spin energy powers the jets. Consequently, its mass does not grow significantly. During this phase, the BH is in general not settled at the center of its host (Bellovary et al. 2019; Pfister et al. 2019). Also, the time evolution of the polar angle $\theta_{\text{BH}}$ indicates that the orientation of the BH spin does not vary significantly but instead, is almost fixed for a long period. On the contrary, the orientation of the angular momentum of the star component, through the evolution of $\theta_{\text{star}}$, exhibits rapid variations which lead to the erratic variations seen in the evolution of $\Psi_{\text{star,BH}}$. This is consistent with Lodato & Pringle (2006) and Lodato & Gerosa (2013) who showed that for low accretion rates, if the Eddington ratio is sufficiently low, the alignment time between BH spin and disk becomes very long.

The second long phase starts once the BH is well settled in the center of the galaxy with a mass close to $10^5 M_{\odot}$, where the

<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Projected distribution of the gas density (left column) and stars ($u-g-r$ band images using SUNSET, right column) for BH-G9685 and its host galaxy at three different epochs. Vectors correspond to the spin of the stellar (cyan), gas accretion disk (red) and BH (black). ① At high redshift, a proto-galaxy is forming and is displaying a disturbed morphology. ② Once the galactic disk is formed, the three spins are well aligned. ③ This Galaxy has a very quiet evolution and turns to a S0 galaxy due to the quenching of star formation. The three spins remains aligned.</figcaption></figure>

spatial separation of its position and the center of its host galaxy tends to zero. The BH mass and spin start to increase again mainly through the accretion of gas that is now more coherent, characterized by high Eddington ratios. The BH spin is also well aligned with the spin of the galaxy and this configuration continues down to $z = 0$. The lower panel of Fig. 4 shows the evolution of the polar angle $\theta$ that describes the BH/gas/star spins direction relative to a fixed coordinate system of the simulation. This angle contains information on how much redirection the BH/gas/star spin vectors have experienced. We note that the 3D orientation of these three vectors are coherently changing over the time during the third phase. The correlation between $J_{\text{star}}$ and $J_{\text{gas}}$ can be indeed explained by the fact that, on one hand, $J_{\text{star}}$ is correlated with the angular momentum of the gas component at galactic scales (a few kpc), which is expected for instance in regular S0 or disk galaxies. On the other hand, we also found that the gas angular momentum vector at those galactic scales shows a significant degree of correlation with the gas angular momentum at
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Evolution of relevant properties of the central black hole in the GALACTICA run G9685. From top to bottom: (1) BH mass evolution. The red and green dashed lines indicate epochs of major and minor mergers respectively; (2) BH spin evolution $|a|$. Red colors indicate when the gas accretion disk is in counter-rotation (i.e., $a < 0$); (3) the Eddington ratio evolution ($\chi$). The red line delimits the quasar dominant mode ($\chi > 0.01$) and the radio dominant mode ($\chi < 0.01$); (4) the evolution of the radiative efficiency of the gas accretion disk ($\epsilon_r$); (5) the variation of the distance separation between the location of the BH and the center of its host galaxy ($d_{star,BH}$); (6) the evolution of the angle $\Psi_{star,BH}$ between the BH spin and the angular momentum of the stellar component (blue) estimated within one effective radius. We also plot the angle between the BH spin and the angular momentum of the gas accretion disk (green) as well as the angle between the stellar and the accreted gas component (red); (7) the evolution of the polar angle $\theta$ describing the orientation of the BH spin (blue), the stellar (red) and gas accretion disk (green) angular momenta relative to a fixed reference frame of the simulation. The fraction of mass gained through BH mergers is $f_{BH,merger} = 1.53\%$. The cosmic evolution of $\Psi_{star,BH}$ follows three different regimes.</figcaption></figure>

scales $4\Delta x \sim 136$ pc (i.e., within which the properties of the gas accretion disk are estimated in the simulation due to the spatial resolution limit), in agreement with Dubois et al. (2014c).

The second example, BH-166, is extracted from the NEWHORIZON simulation. The projected images of the stellar and gas distributions of its host galaxy at different redshifts are

displayed in Fig. 5. The evolution of its main properties is plotted in Fig. 6. The BH was born at $z = 9.4$ ($T_{Universe} = 0.53$ Gyr) and became quite massive, $M_{BH} = 1.4 \times 10^8 M_{\odot}$ at $z = 0.18$ (the last output of simulation). It is hosted by a galaxy with a final mass of $2.4 \times 10^{11} M_{\odot}$. Its evolution significantly differs from BH-G9685 due to the frequent mergers experienced (3 major and
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Same as Fig. 3 but for BH-166 and its host galaxy. We can notice the presence of other BHs (black dots) that are going to merge soon with the main one (associated with the spin vector). ① At high redshift, the proto-galaxy has again a disturbed morphology. The accretion of gas onto the BH is rather chaotic. ② Here also, once the galactic disk is formed, the three spins are well aligned. ③ Although the main BH experiences several mergers and its host galaxy one minor merger, the galactic disk is stable with all spin remaining aligned.</figcaption></figure>

9 minor mergers). As a result, it acquired a high fraction of the mass through the BH mergers, $f_{BH,merger} = 27.3\%$ at $z = 0.18$.

Nevertheless, the evolution of the main properties of BH-166 is quite similar to BH-G9685. Just after its birth, the amplitude of its spin shortly reaches the value $a = 0.42$. We note that from its birth to $T_{\text{Universe}} = 0.81$ Gyr, we were unable to estimate the variations of $\Psi_{\text{star,BH}}$ as this BH is not associated to a galaxy yet following our selection criteria (i.e., it is a “wandering” BH). BH-166 also experienced a long phase of an inefficient gas accretion until $T_{\text{Universe}} \sim 2.5$ Gyr. During this period, its mass grows merely through BH mergers, and the BH spin magnitude, along with its direction, does not evolve significantly. Yet, the variations of $\Psi_{\text{star,BH}}$ are erratic simply due to the rapid changes of the stellar component spin, as suggested by the evolution of $\theta_{\text{star}}$.

After its mass exceeds $10^5 M_{\odot}$ at $z \sim 1.9$ ($T_{\text{Universe}} \sim 3.5$ Gyr), the BH starts to be settled in the host galaxy as suggested by the evolution of $d_{\text{star,BH}}$, and evolves through a smooth gas accretion and minor merger episodes.

At $z = 0.27$ ($T_{\text{Universe}} \sim 10.6$ Gyr), BH-166 experienced a major merger with BH-796 (merger mass ratio: 1–1.4). This

resulted in the decrease of the spin amplitude and misalignment of the spin orientation with respect to the galaxy spin (up to 50°). However, this misalignment lasts only for a short period ($\Delta t \sim 100$ Myr), after which the BH spin is rapidly realigned to the spin of the host galaxy. As in the case of BH-G9685, the spin orientations of the BH, accreted gas, and the galaxy are coherently changing over time during this last phase.

The third example, BH-1049, is also a massive object with $M_{\text{BH}} = 1.8 \times 10^7 M_{\odot}$ at $z = 0.18$, hosted by a galaxy with a final mass of $5.3 \times 10^{10} M_{\odot}$. Images of the projected distribution of stars and gas at different redshifts are displayed in Fig. 7, and the evolution of its main properties is presented in Fig. 8. Its mass evolution is dominated by a smooth accretion with one major merger at $z \sim 2.6$ and four subsequent minor mergers, leading to a value of $f_{\text{BH,merger}} = 0.79\%$ at $z = 0.18$. Right after being seeded at $z = 7.5$ ($T_{\text{Universe}} = 0.7$ Gyr), its spin magnitude increases until it reaches a maximum value ($a = 0.998$) with $\Psi_{\text{star,BH}} < 50°$. Similarly to the previous two cases, BH-1049 crosses a regime dominated by chaotic gas accretion until it is completely settled at the center of its host galaxy at $T_{\text{Universe}} \sim 3.6$ Gyr ($M_{\text{BH}} > 10^5 M_{\odot}$). Since then, its spin is well aligned with $J_{\text{star}}$, and this configuration remains for about 5.6 Gyr until $T_{\text{Universe}} \sim 9.2$ Gyr.

Later, its evolution exhibits a much more complex behavior. One can notice from the time evolution of $\Psi_{\text{star,gas}}$ and $\Psi_{\text{star,BH}}$ at $T_{\text{Universe}} \sim 7.8$ Gyr, that the gas accretion disk does indeed start a retrograde rotation relative to the BH spin and stellar component. This state lasts until $T_{\text{Universe}} \sim 9.8$ Gyr. During this phase, the BH spin also starts to be misaligned with that of the galaxy, and eventually becomes fully anti-aligned. Then, it becomes gradually realigned with the galaxy spin down to $z = 0.18$. In fact, the host galaxy undergoes two successive phases of the cosmic gas accretion from larger scales. As a result, two gas disks with different rotations emerged and coexisted: the original one in the inner parts (and corotating with the stellar disk) and a second one in the outer parts which is counter-rotating. The original one is progressively replaced by the new counter-rotating disk until it totally disappears at $T_{\text{Universe}} \sim 7.8$ Gyr when the gas accretion disk spin flips. Meanwhile, such retrograde accretion decreases the spin magnitude toward $a = 0$ and then grows increasingly negative $a < 0$, which drives the BH spin orientation toward anti-alignment.

Once the BH spin and the angular momentum of the gas are re-aligned ($9.4 \lesssim T_{\text{Universe}} \lesssim 9.8$ Gyr), the BH spin, radiative efficiency and mass significantly increase. However, due to the second phase of cosmic accretion of gas, the same phenomenon happens (coexistence of two gas disks) and the situation reverses: once the inner disk has been replaced by the new gas disk ($T_{\text{Universe}} \sim 9.8$ Gyr), the angular momentum of the gas accretion disk becomes aligned with that of the stellar component (second spin flip) until the end of the simulation. Consequently, the BH spin and accretion disk angular momentum start to be anti-aligned again ($T_{\text{Universe}} \sim 9.8$ Gyr). In the final phase of the evolution of $\Psi_{\text{star,BH}}$, the accretion disk drives the orientation of the BH spin toward alignment. While this example may be an extreme case, it illustrates very well how the long phase of cosmic inflow of gas from larger scales affects the orientation of the BH spin. The detailed study of this galaxy will be presented in a separate paper (Peirani et al., in prep.) which focuses on the formation and evolution of counter-rotating gas and stellar disks.
<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Same as Fig. 4 but for the Cosmic evolution of BH-166. Its evolution differs from BH-G9685 mainly by a higher number of encountered major and minor mergers during its life, respectively 3 and 9. Hence, the mass gained through BH mergers is here much higher, with a value of $f_{BH,merger} = 27.3\%$ at $z = 0.18$. Note that from its birth to $T_{Universe} = 0.81$ Gyr, we were unable to estimate the variations of $\Psi_{star,BH}$ as this BH is not affected to a galaxy yet following our selection criteria (it is a “wandering” BH).</figcaption></figure>

The fourth and last example concerns the BH-549. Here again, Fig. 9 shows images of the projected distribution of stars and gas at different redshifts. The evolution of its main properties is plotted in Fig. A.1. BH-549 is an interesting case since the time evolution of $\Psi_{star,BH}$ indicates that its spin is anti-aligned with the galaxy spin during the last ~4.5 Gyr. We find that the anti-alignment is triggered by the merger with another galaxy at $T_{Universe} \sim 7.2$ Gyr that destroyed the galactic gas disk, as also demonstrated by Park et al. (2019) by studying other NEWHORIZON galaxies. Then, a new disk of gas is rebuilt. Due to the specific impact parameter of the merger, however, the new

gas disk is counter-rotating with respect to the original stellar component. This leads to a stable configuration in which the BH and galaxy spins are anti-aligned over ~4.5 Gyr. It is worth mentioning that Capelo & Dotti (2017) found similar trends using idealized simulation of galaxy mergers while highlighting the role of the ram pressure shocks in the modification of the angular momentum budget of the gas distribution. Appendix A presents the evolution of other BHs as well. Our analysis of typical cases identified the occurrence of three successive phases that an intermediate mass BH ($>10^6 M_{\odot}$) experiences during its lifetime. We summarize those phases in Fig. 10 that shows the coevolution of
<figure><img src="image_7.png" /></figure>

**Fig. 7.** Same as Fig. 3 but for BH-1049 and its host galaxy. This system has a much more complex evolution than BH-G9685 and BH-166 as the host galaxy is passing by two successive phases of cosmic gas accretion with a retrograde rotation with respect to the galactic disk. We illustrate here the first phase: ① Due to a cosmic inflow of gas, there is the coexistence of two disks of gas rotating in opposite direction. The central and original one (dark blue) is progressively disappearing and replacing by the new disk (light blue) from the outer parts. All the spins are still aligned. ② Once the original disk has completely disappeared, The spin of the accreted gas (red) flips to follow the rotation of the newly formed disk. This latter one is counter-rotating with respect to the stellar component. The BH is still aligned with the spin of the galaxy. ③ After some time, the BH aligns with the spin of the gas accretion disk and becomes anti-aligned with the galaxy spin. We note that a new episode of (retrograde) cosmic gas accretion is already in progress and will repeat the same scenario.

three BH-galaxy systems using the color code according to the value of $\cos(\Psi_{\text{star,BH}})$.

### 3.2.2. Statistical analysis

We now turn to examining the statistical evolution of $\Psi_{\text{star,BH}}$ using the whole samples of BHs at different redshifts. In particular, we examine if the distribution of $\cos(\Psi_{\text{star,BH}})$ exhibits a departure from the uniform distribution expected from the isotropic BH spin orientation.

Figure 11 plots the normalized distribution of $\cos(\Psi_{\text{star,BH}})$ for our sample of primary BHs at three different redshifts. The

left panels are derived from our fiducial samples of $M_{\text{BH}} > 2 \times 10^4 M_{\odot}$. The general tendency of alignment ($\cos(\Psi_{\text{star,BH}})$ close to 1) is clearly seen between the black hole spin and the angular momentum of the host galaxy at high ($z = 2$), intermediate ($z = 1$) and low redshifts ($z = 0.18$). In Sect. 3.3, we study in detail to what extent the values of $\Psi_{\text{star,BH}}$ depend on the system properties including the BH mass, the host galaxy mass, the BH spin magnitude, the fraction of mass gained through BH mergers, the distance between the BH position and the center of the galaxy, and the morphology of host galaxies.

We confirmed that changing the definition of the radius within which we estimate the angular momentum of the host galaxies (i.e., $R_e$) has no significant impact on the distribution. For instance, the same trends are also seen even if we compute the angular momentum of the galaxies within $5 R_e$ (black lines) and $0.2 R_e$ (red dotted lines). It is also the case when we take into account all star particles selected by ADAPTAHOP, as it was considered in Dubois et al. (2014b) and Beckmann et al. (2024). For clarity, we do not show the results in Fig. 11. We note that for the result with $0.2 R_e$ at high redshift, the alignment seems to be less pronounced than for those with $R_e$ and $5 R_e$, but it is likely due to the resolution limit required to estimate $J_{\text{star}}$ accurately for regions within $0.2 R_e$.

Finally, we also investigate how the results depend on the lower limit of the BH mass. The right panels of Fig. 11 show the distributions for BHs with $M_{\text{BH}} > 4 \times 10^4 M_{\odot}$. They present an even stronger trend of the BH-galaxy spin alignment. Although our study is focused on primary BHs, it is also interesting to take a look at trends obtained from the sample of secondary BHs with $M_{\text{BH}} > 2 \times 10^4 M_{\odot}$ at $z = 0.18$. Only 30 BHs satisfy the criteria. Secondary BHs are indeed mainly off-centered from the galaxy center by definition, and therefore their mass growth is very modest due to an inefficient gas accretion. In this case, the distribution of $\cos(\Psi_{\text{star,BH}})$ (Fig. 12) tends to be more uniform (i.e., there is no preferential orientation between the BH and galaxy spins), in contrast to the primary BHs. This may be explained by the fact that off-centered BHs have generally low accretion rates and low Eddington ratios (see for instance the evolution of BH-132 in Fig. A.1). Consequently, and as mentioned before, the realignment between BH spin and the accretion disk takes a long time (Lodato & Pringle 2006; Lodato & Gerosa 2013). Thus, the direction of the BH spin remains virtually unchanged. On the contrary, the orientation of the stellar spin changes more rapidly and can explain the distribution in Fig. 12.

A complementary way to characterize the statistics of $\Psi_{\text{star,BH}}$ is the cumulative distribution of $\cos(\Psi_{\text{star,BH}})$, from which one can easily extract the percentage of the aligned BH-galaxy systems below a given angle. The upper panel of Fig. 13 plots those at $z = 3, 2, 1$, and $0.18$. It appears that the BH-galaxy systems tend to be more aligned at higher redshifts. About ~44% of our BH sample have $\Psi_{\text{star,BH}} \le 30^\circ$ at $z = 3$, while this number drops to ~24% at $z = 0.18$. Yet, we notice that the trend is reversed between $z = 1$ and $z = 0.18$. To see more clearly the variation of such quantities, we plot in the lower panel of Fig. 13, the time evolution of the fraction of BH-galaxy pairs that satisfy $\Psi_{\text{star,BH}} \le 45^\circ$ (green line), $\Psi_{\text{star,BH}} \le 30^\circ$ (blue line) and $\Psi_{\text{star,BH}} \le 15^\circ$ (red line). Two different regimes are clearly visible. First, until $z \sim 1.5$, the fraction of aligned BH-galaxy pairs is decreasing, independently of the choice of the $\Psi_{\text{star,BH}}$ maximum value. Then it seems to remain approximately constant until the end of the simulation. However, we note a slight increase from $T_{\text{Universe}}$ close to 11 Gyr. We subsequently see in the next section that these variations are strongly correlated with the BH merger rate. However, it is worth mentioning that we
<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Same as Fig. 4 but for the cosmic evolution of BH-1049. This example displays a phase in which the accretion of gas onto the BH has a retrograde accretion ($7.8 \le T_{\text{Universe}} \le 9.8$ Gyr) which drives the orientation of $J_{\text{BH}}$ toward anti-alignment with respect to the stellar component.</figcaption></figure>

cannot exclude that adding an extra refinement level in the simulation at $z \sim 1.5$ and $z \sim 0.25$ might cause the spurious artifacts in the evolution shown in the lower panel of Fig. 13. The refinement results indeed in a sudden better force resolution and an enhanced gas condensation, which may temporally impact the gas accretion in the BH.

### 3.2.3. Effect of BH mergers

The individual evolution of BHs presented in Figs. 4–8 and in Appendix A suggests that the BH mergers suddenly decreases the BH spin magnitude as well as the rapid change of its 3D orientation. We dedicate this section to study their impact statistically.

Let us first examine the correlation of $\cos(\Psi_{\text{star,BH}})$ with the BH mass and the BH spin parameter. It is shown in Fig. 14 for primary BHs with a mass greater than $2 \times 10^4 M_{\odot}$, where the color of the symbols indicates the percentage of mass gained through BH-merger ($f_{\text{BH,merger}}$). At $z = 3$, most BHs with $M_{\text{BH}} \sim 10^4 M_{\odot}$ (i.e., have doubled or tripled their seed mass) have in general $|a|$ close to unity; they are still in the process of spinning up through early gas accretion as suggested by $f_{\text{BH,merger}}$ close to 0%. This behavior was already seen in the individual evolution of black hole right after their birth in Figs. 4 and 8 (see also Fig. 23 of Dubois et al. 2021). During this phase, their spin is more likely to be aligned with the angular momentum of the stellar component. We also note the existence of a population of BHs with very low $f_{\text{BH,merger}}$ and high $\Psi_{\text{star,BH}}$ values.
<figure><img src="image_9.png" /></figure>

**Fig. 9.** Same as Fig. 3 but for BH-549 and its host galaxy. This system is displaying an interesting feature since from high redshift to $z \sim 0.7$, there is no stable disk and coherent gas accretion onto the BH. However, the host galaxy merges with another galaxy at this specific epoch as illustrated in the panel Ⓓ. After the merger, a gas of disk if forming but rotating in counter-rotation with respect to the stellar component. As a result, the BH spin aligns with the accreted material angular momentum vector and becomes anti-aligned with the galaxy spin as illustrated in Ⓒ.

This population of BHs most likely experienced the phase of chaotic gas accretion seen again in the individual evolution of BHs. We also note that other BHs have already undergone major mergers (red filled circles), which tend to decrease the magnitude of their spin while increasing the scatter in the distribution of $\cos(\Psi_{\text{star,BH}})$. At low redshifts, the same trend is observed in the $\cos(\Psi_{\text{star,BH}}) - |a|$ diagram; BH mergers tend to decrease the magnitude of the BH spin. Regarding the correlation between $\cos(\Psi_{\text{star,BH}})$ and $f_{\text{BH,merger}}$ no clear trend is visible. In fact, even BHs with high $f_{\text{BH,merger}}$ can have their spin well-aligned with their host. This is because their spin can be realigned rapidly as suggested by the evolution of BH-166 in Fig. 6 and several BHs in Appendix A. This is the case for the two most massive black holes of our sample, BH-166 and BH-455, which have a high $f_{\text{BH,merger}}$ value (27.3% and 23.8% respectively) while having low values of $\Psi_{\text{star,BH}}$ (7.7° and 4.9° respectively).

We also study the correlation of $f_{\text{BH,merger}}$ with $\cos(\Psi_{\text{star,BH}})$ in the first column of Fig. 18, at $z = 2, 1$ and $0.18$. A visual inspection indicates that BHs which have experienced a few

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

**Fig. 10.** Variations of the BH mass vs the galaxy mass with color coding according to the value of $\cos(\Psi_{\text{star,BH}})$. We show three typical evolutions of intermediate mass BHs to illustrate the three different regimes followed by the variations of $\Psi_{\text{star,BH}}$: (1) a tendency to BH-galaxy spin alignment right after their birth and for relatively short periods ($M_{\text{BH}} < 3 \times 10^4 M_{\odot}$ and $M_{\text{star}} \lesssim 2 \times 10^7 M_{\odot}$), (2) erratic variations ($M_{\text{BH}} \lesssim 10^5 M_{\odot}$), (3) new tendency to alignment at late time ($5 \times 10^9 M_{\odot} \lesssim M_{\text{star}}$). Note that in the case of BH-1049, the blue region at late time corresponds to the BH spin anti-alignment situation.

merger episodes in their life, for instance those with $f_{\text{BH,merger}} < 10\%$ (blue filled circles), are more likely to have their spin aligned with the angular momentum of their host galaxy. This trend is clearly visible at every redshift. On the contrary, BHs that have undergone more merger episodes and characterized by high $f_{\text{BH,merger}}$ values (e.g., $f_{\text{BH,merger}} > 10\%$, red filled circles), have their spins more uniformly distributed at $z = 2$ and $z = 1$. This is not the case at $z = 0.18$ because, as explained in the next paragraph, BH spins tend to be realigned to the galaxy spins since the subsequent BH mergers occur rarely.

Finally, Fig. 15 plots the BH merger rate of our primary BHs in the volume of the studied simulation. We note first an increase of the number of BH mergers prior to $z = 2$, which is associated with a period of BH formation. Between $z = 2$ and $z = 1$, on the other hand, the number of BH mergers decreases, and becomes nearly constant after $z = 1$. Such trends are clearly correlated with the global evolution presented in the lower panel of Fig. 13. At high redshift, the number of BHs that are aligned below a threshold angle decreases due to the high BH merger rate. At lower redshifts, BH mergers become rarer and have less impact on the variations of $\Psi_{\text{star,BH}}$. Thus, the BH spin becomes progressively realigned with the galaxy spin. This also explains the trend in the left panels of Fig. 18 at $z = 2$ and $z = 1$: BHs that have undergone more merger episodes (red filled circles) have generally exhibited the spin misalignment relative to the galaxy spin. Between $z = 1$ and $z = 0.18$, as BH mergers become less frequent, their impact of the BH-galaxy spin alignment is reduced. Thus, during the last 5 Gyr, BH spins become realigned with their host galaxy spin. This explains why at $z = 0.18$ even merger-dominated BHs can have their spin aligned with
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Normalized distributions of $\cos(\Psi_{\text{star,BH}})$, the angle between the BH spin, and the angular momentum of the stellar component of the host galaxy. The latter was computed within spheres of radius $R_e$ (blue histograms), $5 \times R_e$ (black lines), and $0.2 \times R_e$ (dotted red lines). The results are shown at three different redshifts ($z = 2, 1$ and $0.18$) as well as two lower BH mass limits, $2 \times 10^4 M_{\odot}$ (left column) and $4 \times 10^4 M_{\odot}$ (right column). In each panel, the number of primary black holes considered to derive the histograms is indicated.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Normalized distribution of $\cos(\Psi_{\text{star,BH}})$ for secondary BHs at $z = 0.18$. The sample consists of 30 BHs with a mass greater than $2 \times 10^4 M_{\odot}$. Contrary to primary BHs, the distribution tends to be uniform suggesting no privileged orientation between the BH spin and the angular momentum of their host galaxies.</figcaption></figure>

$J_{\text{star}}$. These results are in good agreement with Beckmann et al. (2024), who found from HORIZON-AGN that BHs dominated by non-merger growth ($f_{\text{BH,merger}} < 10\%$) are more likely to be aligned to their galaxy spins than BHs dominated by merger growth.

### 3.3. Dependence on BH and galaxy properties

Finally, we investigate the potential dependence of $\cos(\Psi_{\text{star,BH}})$ on the BH and galaxy properties. Let's consider first the correlation of $\cos(\Psi_{\text{star,BH}})$ and the BH mass. The normalized distribution of $\cos(\Psi_{\text{star,BH}})$ for our samples of primary BHs are shown in the left panels of Fig. 16. For clarity, we have divided our sample into two groups: BHs with mass lower (blue filled circles) or greater (red filled circles) than $10^5 M_{\odot}$. As already found in the last sections, BHs with $M_{\text{BH}} > 10^8 M_{\odot}$ preferentially show the spin alignment with $J_{\text{star}}$. This result is consistent with

Dubois et al. (2014b). However, they also found that for more massive BHs ($M_{\text{BH}} > 10^8 M_{\odot}$), spins seem to be more randomly oriented with respect to their host galaxy angular momentum. Their finding was not confirmed in the present analysis, mainly because our sample lacks in very massive black holes due to the limited volume size. The spins of less massive BHs are more likely to be aligned at high-$z$ since most of them are still in the process of spinning up through an early efficient gas accretion. This effect is less visible for intermediate mass BHs at low redshifts, as already indicated in Fig. 2.

Second, we examine the correlation between $\cos(\Psi_{\text{star,BH}})$ and the host galaxy mass; see middle panels of Fig. 16 for two separated samples of galaxies selected with a mass lower (blue circles) or greater (red circles) than $10^9 M_{\odot}$. In general, the two subsamples present a quite similar trend that BH-galaxy spins are more likely to be aligned in both samples. However, as already noticed in Fig. 2, the trend is more prominent for low-mass galaxies at high redshifts; see the blue histogram at $z = 2$.

Finally, the right panels of Fig. 16 plot the correlations between $\cos(\Psi_{\text{star,BH}})$ and the dimensionless BH spin parameter. At high redshift ($z > 1$), BHs with a higher spin magnitude ($|a| > 0.8$) tend to have a spin more aligned with the galaxy spin than those with $|a| < 0.8$. This is consistent since BHs spin up rapidly due to an efficient gas accretion (Volonteri et al. 2013). In this case, the angular momentum transferred to the black hole at the accretion phase is sufficient to keep the BH spin mostly aligned, in agreement with Bustamante & Springel (2019). At low redshifts, the trends are even more pronounced. This is due to the fact that minor/major mergers become less frequent. Consequently, BH spin magnitudes tend to increase through a coherent gas accretion, which efficiently aligns their spin to the angular momentum vector of the stellar component.

Additionally, we consider the effect of retrograde gas accretion disks. Figure 17 shows the correlation of $\cos(\Psi_{\text{star,BH}})$ and
<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Cumulative distributions of $cos(\Psi_{star,BH})$ at different redshifts (upper panel). The black line represents the cumulative distribution of a uniform distribution. In the lower panel, we show additionally show the evolution of the percentage of BHs or BH-galaxy systems that present an alignment angle $\Psi_{star,BH}$ lower than 45° (red line), 30° (blue line), or 15° (red line). The two vertical dashed lines indicate the epochs where an additional level of refinement was added in the simulation.</figcaption></figure>

a; the positive and negative signs of a indicate the prograde (red) and retrograde or counter-rotating (blue) gas disks, respectively. The two histograms in the upper panel indicate counter-rotating gas disks generally lead to BH-galaxy spins less aligned.

Middle panels of Fig. 18 plot the distances between the BH positions and the center of their host galaxy, $d_{star,BH}$, against $cos(\Psi_{star,BH})$. Those BHs located closer to the center of their host galaxy (<1 kpc) have generally their spin well-aligned with $J_{star}$.

On the contrary, for off-centered BHs, the values of $cos(\Psi_{star,BH})$ tend to be more uniformly distributed. This trend was already noticed when studying a sample of secondary BHs that are off-centered according to our definition.

Finally, we look into the dependence of $cos(\Psi_{star,BH})$ on the morphology type of the host galaxies. For each galaxy, we have computed the radial, tangential and vertical velocity components of stellar particles. In doing so, we choose the z-axis of the cylindrical coordinate as the total angular momentum vector direction of the stellar component. Then, we estimate the rotation velocity $V$ of the galaxy by averaging the tangential velocity component. We also compute the velocity dispersion $\sigma$ from the dispersion of the radial $\sigma_r$, the tangential $\sigma_{\theta}$ and the vertical velocity $\sigma_z$ components with respect to their averaged values; $\sigma^2 = (\sigma_r^2 + \sigma_{\theta}^2 + \sigma_z^2)/3$. In what follows, we use the quantity $V/\sigma$ as a proxy of galaxy morphology.

The right panels in Fig. 18 show the correlation between $cos(\Psi_{star,BH})$ and $V/\sigma$. In those panels, we select more rotation-dominated galaxies with $V/\sigma > 0.6$, and more dispersion-dominated galaxies with $V/\sigma < 0.6$. In general, disk galaxies (i.e., objects with high $V/\sigma$ values) clearly have their spin well-aligned with that of the associated BH spin. The trend is however less pronounced for more spheroidal galaxies, especially at high redshifts. This can be explained by the fact that BH spins tend to align with the angular momentum of the accreted material, which is in general related to the large scale kinematics of the galaxy. To this regard, Sesana et al. (2014) found that different galaxy morphologies result in different spin distributions, in particular that BHs hosted in ellipticals tend to have lower spins than those hosted in spirals. Another possible explanation for this discrepancy is that mergers contribute more significantly to the build-up of such galaxies at high redshifts (e.g., Boylan-Kolchin et al. 2006; Trujillo et al. 2011; Volonteri & Ciotti 2013). It could also be a selection effect since galaxies with low $V/\sigma$ values tend to be less massive and characterized by a more turbulent gas. Consequently, BHs lie in a gas environment owing a more random angular momentum which affects their spin direction. Yet, at $z = 0.18$, the trend looks quite similar for the two subsamples. This is consistent with Beckmann et al. (2024) who have derived the distribution of $cos(\Psi_{star,BH})$ from the HORIZON-AGN simulation at $z = 0$ by classifying the galaxy morphology according to the bulge-to-total ratio (B/T). They found that bulgeless systems (with $B/T < 0.1$) exhibit no statistical difference in the BH-galaxy spin alignment relative to the rest of the sample ($B/T > 0.1$).

# 4. Projected angles statistics

From an observational point of view, the accurate measurement of the “true” 3D angle $\Psi_{star,BH}$ may prove to be a difficult task since the stellar component as well as spin vectors are “seen” in 2D projection on the sky. Instead, it might be relevant to focus on the so-called 2D projected angles $\lambda$ (e.g., Ohta et al. 2005; Benomar et al. 2014), which is the angle between the projected BH spin and the galaxy spin on the sky as schematically shown in Fig. 19. As for $\Psi$, $\lambda$ lies between 0 and 180°.

We therefore would like to present relevant statistics of $\lambda$ which can be directly derived from our 3D statistical analysis. This could potentially help the interpretation of future observations of BH-galaxy spin alignment through the measurement of $\lambda$, for instance by combining the projected velocity field (to estimate the direction of $J_{star}$) and the direction of jet to infer the direction of the BH spin.
<figure><img src="image_14.png" /><figcaption>**Fig. 14.** Correlation of $\cos(\Psi_{\text{star,BH}})$ against the BH mass (for the primary BH with a mass greater than $2 \times 10^4 M_{\odot}$; left panels) and the BH spin magnitude (right panels), with color coding according to the percentage of mass gain through BH-merger $f_{\text{BH,merger}}$ at each specific redshift. At $z=3$, most BHs whose mass is close to their seed mass ($\sim[1-3] \times 10^4 M_{\odot}$) are still in the process of spinning up through early gas accretion (i.e., blue filled circles). Their dimensionless spin parameters are close to unity, and tend to be aligned with the stellar component angular momentum. Other BHs have already undergone major mergers (red and green filled circles) which tend to decrease the magnitude of their spin while increasing the scatter in the distribution of $\cos(\Psi_{\text{star,BH}})$. At low redshift, the same trend is observed in the $\cos(\Psi_{\text{star,BH}})-|a|$ diagram. The corresponding correlations between $|a|$ and the BH mass can be directly seen in Fig. 23 of Dubois et al. (2021).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 15.** Evolution of the merger rates of primary black holes with a mass greater than $2 \times 10^4 M_{\odot}$ in the simulated volume (green line). The red and blue lines correspond respectively to major and minor mergers delimited by a 1:4 mass ratio. We have slightly smoothed the curves to emphasize the general trends.</figcaption></figure>

We recall that for a 3D isotropic distribution of BH spins, the distribution of $\cos(\Psi_{\text{star,BH}})$ is uniform. In 2D, it is the distribution of $\lambda$ which is uniform and for this reason, we study in the

following the probability density function $P(\lambda|\Psi, i_o)$, where $i_o$ is the galaxy disk inclination relative to the observer.

### 4.1. 3D and 2D projected angle relation

We use a simple Monte-Carlo method to derive the projected angle from a given BH-galaxy pair of our samples with a random orientation in space as well as a random $i_o$. We describe the main steps here. First, for a given studied BH-galaxy pair, we have rotated both the angular momentum of the stellar component (galaxy spin, $J_{\text{star}}$) and the associated BH spin from their original 3D orientation in space in the simulation so that $J_{\text{star}}$ coincides now exactly with the z-axis of the frame presented in Fig. 19. Then we have randomly drawn an angle between 0 and $2\pi$ and rotate accordingly the BH spin with respect to the z-axis to generate a new configuration. This operation does not modify the original value of the 3D angle $\Psi_{\text{star,BH}}$. The final step consist in applying an orbital inclination $i_o$. To do so, we randomly draw an angle between 0 and $\pi/2$ and rotate the BH spin with respect to the y-axis. The projected angle is then the angle between the z-axis and the projection of the 3D BH spin onto the $x=0$ plane.

The relation between the 3D ($\Psi$) and 2D projected ($\lambda$) angles has already been studied in works related to the stellar spin-orbit misalignment angle in planetary systems (Fabrycky & Winn 2009; Crida & Batygin 2014). But in the majority of cases, $i_o$
<figure><img src="image_16.png" /><figcaption>**Fig. 16.** Variations of $\cos(\Psi_{\text{star,BH}})$ from our sample of primary BHs with a mass $M_{\text{BH}} > 2 \times 10^4 M_{\odot}$ with respect to the BH mass (first column), the host galaxy mass (second column), and BH spin (third column) at three different redshifts ($z = 2, 1$, and $0.18$). In each panel, we divided the BH sample into two subsamples, characterized by the blue and red colors, and show the corresponding distributions of $\Psi_{\text{star,BH}}$ in the upper part of the panels. Here, the angular momentum of each associated galaxy is calculated within spheres of radius $R_e$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 17.** Variations of $\cos(\Psi_{\text{star,BH}})$ with respect to the BH spin magnitude ($a$) at $z = 0.18$. Negative values of $a$ is associated with a retrograde gas accretion.</figcaption></figure>

is close to $\pi/2$ since the observer is almost exactly in the orbital plane of the transiting planets. In that case, a simple analytical expression of the probability density function of $\lambda$ for fixed $\Psi$, $P(\lambda|\Psi, i_o = \pi/2)$, can be derived: Eq. (4) in Crida & Batygin (2014) or Eq. (19) in Fabrycky & Winn (2009).

To test the validity of our sampling, we compute $P(\lambda|\Psi, i_o = \pi/2)$ for two BH-galaxy pairs from our sample at $z = 0.18$ for which $\Psi_{\text{star,BH}} = 7.7^\circ$ (BH-166) and $\Psi_{\text{star,BH}} = 153.2^\circ$ (BH-549), respectively. For each of them, we have generated 10 000 random realizations using the steps described above, and computed the corresponding projected angles. The probability density functions are shown in Fig. 20. The good agreement with the analytical expressions derived by Crida & Batygin (2014) ensures that our sampling is reliable. Moreover, it is instructive to see that for $i_o = \pi/2$, the projected angle can take many values ranging from $[0-\Psi]$ (if $\Psi < \pi/2$) or from $[\Psi-\pi]$ (if $\Psi > \pi/2$). We have also checked that our sampling is in good agreement with analytical solutions for more complex cases ($i_o \neq \pi/2$). For instance we have correctly reproduced the probability functions $P(\lambda|\Psi, i_o = 80^\circ)$ presented in Fig. 3 of Fabrycky & Winn (2009).

## 4.2. 2D projected angles statistics

Now, since $\cos(\Psi_{\text{star,BH}})$ and therefore $\Psi_{\text{star,BH}}$ have their own distribution at $z = 0.18$ (the lower left panel of Fig. 11 for $\cos(\Psi_{\text{star,BH}})$), what is the corresponding distribution of $\lambda$?

To answer this question, we have applied the Monte-Carlo method to the 102 BH-galaxy pairs of our sample at $z = 0.18$. For each original BH-galaxy pair, we have generated 500 random configurations and computed the corresponding projected
<figure><img src="image_18.png" /><figcaption>**Fig. 18.** Same as Fig. 16 but for the variations of $\cos(\Psi_{\text{star,BH}})$ with respect to the percentage of BH mass gained through BH mergers ($f_{\text{BH,merger}}$, first column), the distance separation between the BH and the center of its host galaxy (second column), and $V/\sigma$ as a proxy for galaxy morphology (third column).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 19.** Schematic view of the projected angle. The x-axis is pointing toward an observer. The angle $i_o$ is the inclination which defines the angle through which the galactic plan is seen by the observer. The projected angle $\lambda$ of the 3D angle $\Psi_{\text{star,BH}}$ is contained in the y-z plane (with $x=0$).</figcaption></figure>

angles. The resulting distribution of $\lambda$ is shown in the upper left panel of Fig. 21. We also show the spline interpolated function to describe it. The tendency of spin alignment seen in the 3D analysis is still observable in the 2D case. In the lower left panel, we show the different distributions obtained using spline functions for different values of $i_o$ (every 10°). Since observational measurements might be preferentially obtained for high $i_o$ values (i.e., close to $\pi/2$), this result could reveal some selection effects. We notice that the tendency of alignment is less and less clear as $i_o$ decreases (more face-on views). Note that we have assumed in our model a uniform distribution of $i_o$. However, it is now well known that galaxies do not form everywhere but in filaments and nodes implying some preferred directions for the orientation of the spins within the cosmic web, as suggested by numerical simulations (e.g., Codis et al. 2018; Kraljic et al. 2020; Zhang et al. 2023) or observational analysis (e.g., Welker et al. 2020; Kraljic et al. 2021; Desai & Ryden 2022). This effect could be taken into account for even more realistic modeling.

If the orientation of jets is used as a proxy to the direction of the projection of any BH spin vector, the sense of its rotation actually cannot be known. Thus, to be more consistent with future observational measurements based on jet orientation, we plot in the right side of Fig. 21 the same statistics as the left part by redefining $\lambda$ in the range between 0 and $\pi/2$. In that case, among the two possible projected angles $\lambda$ and
<figure><img src="image_20.png" /><figcaption>**Fig. 20.** Probability density function $P(\lambda|\Psi, i_o = \pi/2)$ of the 2D projected angle ($\lambda$) for two fixed values of $\Psi_{star,BH}$, 7.7° (top panel) and 153.2° (lower panel). In the two cases, we have considered an orbital inclination of $i_o = \pi/2$. The histograms are the results obtained from a Monte Carlo method (using 10 000 random lines-of-sight) while the red lines are analytical solutions.</figcaption></figure>

$\pi-\lambda$, we always choose the smallest angle which ranges from 0 to $\pi/2$. We note that the tendency of spin alignment is still present but less pronounced than when taking into account the sense of the BH spins. Note also that, for a higher consistency, only BHs with a low Eddington ratio at the studied redshift (or in a redshift interval) should be considered. For instance, at $z = 0.18$, there are only 5 and 8 BHs that have a Eddington ratio greater that 0.01 and 0.001 respectively. We have checked that the main trends displayed in the right part of Fig. 21 remain similar even if we remove these BHs from the sample. Additionally, we have considered the whole sample of BHs (102 in total) within a small redshift interval (for instance $0.18 < z < 0.28$ which corresponds to the last 1 Gyr), and focused on those BHs when their Eddington ratio was low (i.e., radio mode). Even in those cases, we made sure that our trends remain the same. Also, it is worth mentioning that the orientation of radio jets might not be necessarily aligned with BH spins because the propagation of jets from the accretion disk up to galactic scale can be strongly perturbed by the multiphase structure of the interstellar medium (e.g., Mukherjee et al. 2018; Cielo et al. 2018; Junor et al. 1999), or because of the jet precession (Dunn et al. 2006; Krause et al. 2019; Ubertosi et al. 2024), or also due to relative motions between the jets and the hot atmosphere in intra-galaxy clusters (e.g., Heinz et al. 2006; Morsony et al. 2010,

2013; O'Dea & Baum 2023, and references therein). All of these effects might affect the observed distribution of $\lambda$.

Additional statistics on the 2D projected angles according to the galaxy morphology are derived in Appendix B.

## 5. Conclusion and discussion

The BH spin is a key parameter in black hole physics as well as in the formation process of the host galaxies. Characterizing the evolution of its amplitude and orientation through cosmic time helps distinguish among different scenarios for BH-galaxy coevolution. Using the NEWHORIZON and GALACTICA simulations, we have investigated the long-term evolution of the 3D angle between BH spins and the angular momentum vectors of their host galaxies. These simulations have sufficient resolution to capture the injection scale of gas turbulence, and has a multiphase interstellar medium. Catalogs of BH-galaxy pairs have been produced by selecting primary BHs with mass greater than $2 \times 10^4 M_\odot$ and the angular momentum of each galaxy has been estimated using all star particles within one effective radius. Our main conclusions can be summarized as follows:

1. The analysis of typical individual evolution of the most massive BHs of our catalogs suggests the occurrence of three phases in the evolution of $\Psi_{star,BH}$:

- Right after being seeded with an initial mass of $10^4 M_\odot$ and no spin, most of BHs are rapidly spun up to a high or a maximum spin value as they double or triple their mass through gas accretion (Fig. 14), in agreement with Dubois et al. (2021), Bustamante & Springel (2019), and Sala et al. (2024). During this relatively short period, BH spins are likely to be well aligned with the spin of their host galaxies. And as a consequence, low-mass BHs at high redshift present in general low values of $\cos(\Psi_{star,BH})$ (Figs. 2, 10, and 14);

- A second phase ($M_{BH} \lesssim 10^5 M_\odot$) is dominated by a rather chaotic and inefficient gas accretion, reflecting the patchy morphology of forming proto-galaxies. In this phase, $\Psi_{star,BH}$ presents erratic variations mainly driven by the rapid change of the stellar angular momentum orientation, rather than that of the BH spin direction (see the evolution of $\theta_{BH}$ and $\theta_{star}$ in the individual evolution shown in Figs. 4–8, A.1, and A.2). This phase is indeed characterized by low accretion rates and low Eddington ratios, which favor a very long alignment time between BH spin and the accreted material (Lodato & Pringle 2006; Lodato & Gerosa 2013).

- The last phase follows, in which the gas accretion becomes much more coherent and in which BHs are generally well settled in the center of their host ($10^5 M_\odot \lesssim M_{BH}$). From the study of the most massive BHs of our samples, the transition happens generally when the host galaxy has a mass ranging from $5 \times 10^9$ to $5 \times 10^{10} M_\odot$. The BH spins are likely to be well aligned with the angular momentum of the galaxy, and this configuration can last over Gyr time periods, even though BH merger episodes can temporally misalign their spin.

2. The majority of the most massive BHs of our samples ($M_{BH} > 10^5 M_\odot$) have BH spins highly aligned with $J_{star}$, even at low redshift $z = 0.18$. This result is consistent with those from Dubois et al. (2014b). However, in that paper they also found that for more massive BHs ($M_{BH} > 10^8 M_\odot$), spins are more randomly oriented with respect to their host galaxy angular momentum and tend toward $\cos(\Psi_{star,BH}) = 0$. In these galaxies, BHs indeed accrete, at very low levels, from the turbulent hot intracluster gas that bears no connection to the angular momentum of the stellar component of the galaxies. Such trend could not
<figure><img src="image_21.png" /><figcaption>**Fig. 21.** Left column: probability density function $P(\lambda|\Psi, i_o)$ of the 2D projected angle ($\lambda$) obtained from the distribution of $\Psi_{star,BH}$ at $z = 0.18$. We also use a spline function to describe the obtained PDF. In the lower left panel, we show the different distributions obtained through spline function description for different values of $i_o$ (every 10°). Right column: we generated the plots using the same statistics but without specifying the sense of the BH vector to maintain consistency with observations.</figcaption></figure>

be confirmed in the present analysis since our sample of very massive black holes is too limited.

3. Statistically, the distribution of $\cos(\Psi_{star,BH})$ at a given redshift indicates that BH spins are more likely to be aligned with their host galaxy spin, in agreement with Beckmann et al. (2024) and with a recent observational analysis (Zheng et al. 2024). This tendency to alignment is more pronounced at high redshift ($z = 3–2$), while it decreases until $z \sim 1.5$ due to BH-merger events, and becomes nearly constant. We also note that at late time, BH spin starts again to be more aligned with the galaxy spin since the BH mergers become less frequent.

4. BH merger events is one of the main drivers of 3D angle misalignment. Our analysis suggests that BHs which have undergone more merger episodes in their history, and characterized by high $f_{BH,merger}$ values, are more likely to have their spin misaligned than those dominated by a smooth gas accretion (Fig. 14 and first column in Fig. 18). This is expected since the orbital angular momentum is redistributed in BH mergers. Our investigation of individual BH history also suggests that after a BH-merger episode, the BH spin tends to realign to the galaxy spin in a relatively short period ($\Delta t \sim 100–200$ Myr). A similar timescale is suggested in Fig. 9 of Bustamante & Springel (2019) regarding the average variations of the $\Psi_{gas,BH}$ after a merger.

5. Long phase of retrograde accretion of gas induced by either cosmic inflows or galaxy mergers can produce strong

effects, by anti-aligning the BH-spin in a long period, as suggested by the evolution of BH-1049 and BH-549. From a statistical point of view, counter-rotating gas disk with respect to the stellar component are more likely to lead to BH-galaxy spin misalignment (Fig. 17). Dubois et al. (2014b) found that about 10–30% of BHs have spins counter-aligned with their host galaxy. This occurs in galaxies where the central regions have been deprived of cold gas or because BH coalescences have flipped the spin direction. Bustamante & Springel (2019) confirmed this statement and found that 20% of BH merger events may lead to momentarily counter-rotating gas disk. The individual evolution of BH-1049 and BH-549 present different scenarios as the retrograde gas accretion onto the BH is not induced by a punctual effect (such as a BH merger). Instead, it is due to the continuous accretion supplied by the galactic gas disk. For this reason, a BH-galaxy spin anti-alignment configuration could persist for an extended period, 0.5 Gyr and >4 Gyr for BH-1049 and BH-547 respectively. The present work therefore predicts that galaxies presenting strong star-gas spin angle misalignment (e.g., from the ManGA survey: Beom et al. 2022; Xu et al. 2022; Katkov et al. 2024) are likely to host BHs owning to a clear misaligned or anti-alignment with respect to the stellar component. We will present in a forthcoming paper, a detailed analysis about the formation of counter-rotating gas/stellar disks (Peirani et al., in prep.).

6. The analysis of the 2D misalignment angle ($\lambda$) revealed that the distribution tends to keep record of the alignment
tendency seen in the 3-dimensional predictions. It might then
be possible to identify such signal from observational analysis
measuring 2D misalignment angles (Lin et al., in prep.).

This paper is complementary to Beckmann et al. (in prep.)
which focuses on the cosmic evolution of BH spin magni-
tude using NEWHORIZON and the same BHs catalogs. They
also found that the evolution of BH spins follow three dis-
tinct phases while highlighting the role of BH mergers during
the second phase which tend to increase the scattering of |q|.
This is consistent as the magnitude and orientation of spins are
connected: high/low values of spin magnitude are generally
associated with BH-galaxy spin alignment/misalignment, as sug-
gested by Figs. 14 and 16.

The theoretical predictions presented in this paper heavily
rely on sub-grid models, which may not completely follow all
relevant physical processes, even in the currently best simula-
tions. In particular, sub-grid modeling of BH physics should be
interpreted with caution. For instance, the gas accretion disk
falling onto the BH is not spatially resolved in our simulation.
Thus, we have assumed that the gas accretion disk around the
BH retains the same angular momentum of gas component esti-
mated at 4Δx ~ 136 pc from the location of each BH. While
this scale is much larger than the physical scale of the accre-
tion disk, it is expected to be a reasonably good approximation
to assume that its angular momentum, especially its orientation,
is conserved during the accretion process of the gas toward the
disk. This is consistent with the result of Maio et al. (2013) using
high resolution simulations down to 1 pc, and even in the pres-
ence of strong star formation feedback. Yet, we note that there
are different suggestions in the literature. Levine et al. (2010)
found that the angular momentum vector on scales ~100 pc may
vary substantially from the direction of angular momentum on
kiloparsec scales between z = 4 and z = 3. This is mainly due to
the interaction of clump of gas moving toward the center of the
disk, which results in dramatic change of the nuclear gas angular
momentum. Hopkins et al. (2012) also predicted a weak correla-
tion between the nuclear axis and the large-scale disk axis, using
high-resolution simulations of gas inflows from galaxy to parsec
scales around AGN. Sudden misalignments may be caused either
by single massive clumps falling into the center slightly off-axis,
or due to gravitational instabilities. The two latter works, how-
ever, focused on gas rich galaxies at high redshifts where the gas
turbulence may still play a strong role in the gas angular momen-
tum alignment.

Another important assumption in the present work is that two
BHs instantly merge when they come close less than 4Δx. In
reality, however, they would first form a binary for a while. If
there is a sufficient amount of gas supply, the BH spins may
evolve during the merger process. If the circumbinary disk is
misaligned with respect to the BH binary orbital plane, it leads
to tearing or warping of the disk structure, which may affect the
spin alignment (or not) between the two BHs (see, for instance,
Nixon et al. 2011, 2013; Gerosa et al. 2015; Moody et al. 2019;
Nealon et al. 2022; Bourne et al. 2023). However, our analysis
suggests that BH spins tend to rapidly align after the BH merger
especially during long and coherent phases of gas accretion.
Thus, we expect that these missing processes should not signifi-
cantly affect our predictions, at least statistically.

We should also note that we have not included any analyt-
ically model for the unresolved dynamical friction experienced
by BHs from stars and DM particles. This has been shown to
play a crucial role in sinking BHs closer to the galaxy cen-
ter, and also in enhancing the BH merger rate (Pfister et al.
2019; Chen et al. 2022; Ma et al. 2023). The dynamical fric-

tion from gas is expected, however, to have a lower impact
especially at increased resolution due to instabilities in the
wake (Beckmann et al. 2018) and the turbulent nature of the
gas (Lescaudron et al. 2023). The absence of dynamical friction
in our model might at some point delay the alignment of the
BH spin since off-centered BHs are more likely to have mis-
aligned spins. On the contrary, a higher BH merger rate should
statistically reduce the tendency of alignment. Thus, more quan-
titative and reliable predictions need the refinement of the sub-
grid physics.

Finally, we have assumed that all BH at birth have an initial
mass of $10^4 M_{\odot}$ with no spin. Depending on the initial BH spin,
however, accretion onto the BH after birth can alter the spin mag-
nitude and/or torque the BH spin alignment (McKernan & Ford
2023). Furthermore, Sala et al. (2024) recently employed a very
similar numerical scheme to compute the BH spin evolution
in idealized or cosmological simulations. They found the rapid
increase of the spin parameter for $5 \times 10^5 M_{\odot} < M_{BH} < 10^6 M_{\odot}$
(see their Fig. 13), which is consistent with ours qualitatively (or
see Fig. 23 of Dubois et al. 2021, for a more direct comparison).
Since the mass range is very close to their initial BH seed mass
of $5 \times 10^5 M_{\odot}$, the result may be suspected to be biased due to
the mass resolution. In fact, our BH seed mass is $10^4 M_{\odot}$, 50
times smaller than theirs, and we found the similar trend around
$10^4 M_{\odot} < M_{BH} < 2 \times 10^4 M_{\odot}$. Even though we selected BHs
of $M_{BH} > 2 \times 10^4 M_{\odot}$, it is likely that they still suffer from
the numerical artefact to some extent. A similar statement can
be drawn from Fig. 7 of Bustamante & Springel (2019) as they
used a seed mass close to ~$10^6 M_{\odot}$. This is why we examined
the behavior of four examples in detail (BH-G9685, BH-166,
BH-1049 and BH-541), which are supposed to be free from such
possible numerical effects, at least in the late epochs ($z < 2$). In
any case, our primary focus is the orientation of the BH spin and
galaxy orbit, and thus the present paper is very complementary
to those previous works that are more focused on the cosmic evo-
lution of the amplitude of the BH spin (e.g., Dubois et al. 2014b,
2021; Bustamante & Springel 2019; Sala et al. 2024).

**Acknowledgements.** We warmly thank the referee for an insightful review that considerably improved the quality of the original manuscript. S.P. acknowledges the support from the JSPS (Japan Society for the Promotion of Science) long-term invitation program, and is grateful for the hospitality at Department of Physics, the University of Tokyo. This research is partly supported by the JSPS KAKENHI grant Nos. JP18H01247 and 23H01212 (Y.S.). S.K.Y. acknowledges support from the Korean National Research Foundation (2020R1A2C3003769, 2022R1A6A1A03053472). This work was granted access to the HPC resources of CINES under the allocations c2016047637, A0020407637, and A0070402192 by Genci, KSC-2017-G2-0003 by KISTI, and as a “Grand Challenge” project granted by GENCI on the AMD Rome extension of the Joliot Curie supercomputer at TGCC. The large data transfer was supported by KREONET which is managed and operated by KISTI. Y.-T.L. is grateful for support from the National Science and Technology Council of Taiwan under grants MOST 111-2112-M-001-043 and NSTC 112-2112-M-001-061. C.P. is partially supported by the grant Segal (https://www.secular-evolution.org/) ANR-19-CE31-0017 of the French Agence Nationale de la Recherche and by the National Science Foundation under Grant No. NSF PHY-1748958. This work was carried within the framework of the Horizon project (http://www.projet-horizon.fr). Most of the numerical modeling presented here was done on the Horizon cluster at Institut d'Astrophysique de Paris (IAP).

References

Albrecht, S. H., Dawson, R. I., & Winn, J. N. 2022, PASP, 134, 082001
Ardila, F., Huang, S., Leauthaud, A., et al. 2021, MNRAS, 500, 432
Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376
Baldassare, V. F., Dickey, C., Geha, M., & Reines, A. E. 2020, ApJ, 898, L3
Barausse, E., & Rezzolla, L. 2009, ApJ, 704, L40
Bardeen, J. M. 1970, Nature, 226, 64
Batiste, M., Bentz, M. C., Raimundo, S. I., Vestergaard, M., & Onken, C. A. 2017, ApJ, 838, L10

Beckmann, R. S., Slyz, A., & Devriendt, J. 2018, MNRAS, 478, 995

Beckmann, R. S., Dubois, Y., Volonteri, M., et al. 2023, MNRAS, 523, 5610

Beckmann, R. S., Smethurst, R. J., Simmons, B. D., et al. 2024, MNRAS, 527, 10867

Bellovary, J. M., Cleary, C. E., Munshi, F., et al. 2019, MNRAS, 482, 2913

Bennert, V. N., Treu, T., Ding, X., et al. 2021, ApJ, 921, 36

Benomar, O., Masuda, K., Shibahashi, H., & Suto, Y. 2014, PASJ, 66, 94

Benson, A. J., & Babul, A. 2009, MNRAS, 400, 2208

Beom, M., Bizyaev, D., Walterbos, R. A. M., & Chen, Y. 2022, MNRAS, 516, 3175

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

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

Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53

Bourne, M. A., Fiacconi, D., Sijacki, D., Piotrowska, J. M., & Koudmani, S. 2023, MNRAS, submitted [arXiv:2311.17144]

Boylan-Kolchin, M., Ma, C.-P., & Quataert, E. 2006, MNRAS, 369, 1081

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

Campanelli, M., Lousto, C., Zlochower, Y., & Merritt, D. 2007, ApJ, 659, L5

Capelo, P. R., & Dotti, M. 2017, MNRAS, 465, 2643

Chabrier, G. 2005, in The Initial Mass Function 50 Years Later, eds. E. Corbelli, F. Palla, & H. Zinnecker, Astrophys. Space Sci. Lib., 327, 41

Chen, N., Ni, Y., Tremmel, M., et al. 2022, MNRAS, 510, 531

Choi, E., Ostriker, J. P., Naab, T., Oser, L., & Moster, B. P. 2015, MNRAS, 449, 4105

Cielo, S., Babul, A., Antonuccio-Dealogu, V., Silk, J., & Volonteri, M. 2018, A&A, 617, A58

Cisternas, M., Jahnke, K., Bongiorno, A., et al. 2011, ApJ, 741, L11

Codis, S., Jindal, A., Chisari, N. E., et al. 2018, MNRAS, 481, 4753

Crida, A., & Batygin, K. 2014, A&A, 567, A42

Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11

Cui, Y., Hada, K., Kawashima, T., et al. 2023, Nature, 621, 711

Davis, B. L., Graham, A. W., & Cameron, E. 2019, ApJ, 873, 85

Desai, D. D., & Ryden, B. S. 2022, ApJ, 936, 25

Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604

Di Matteo, T., Colberg, J., Springel, V., Hernquist, L., & Sijacki, D. 2008, ApJ, 676, 33

Ding, X., Silverman, J., Treu, T., et al. 2020, ApJ, 888, 37

Dong-Páez, C. A., Volonteri, M., Beckmann, R. S., et al. 2023, A&A, 673, A120

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

Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2010, MNRAS, 409, 985

Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662

Dubois, Y., Pichon, C., Welker, C., et al. 2014a, MNRAS, 444, 1453

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

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

Dubois, Y., Volonteri, M., Silk, J., et al. 2015, MNRAS, 452, 1502

Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948

Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109

Duffy, A. R., Schaye, J., Kay, S. T., et al. 2010, MNRAS, 405, 2161

Dunn, R. J. H., Fabian, A. C., & Sanders, J. S. 2006, MNRAS, 366, 758

Fabrycky, D. C., & Winn, J. N. 2009, ApJ, 696, 1230

Fanidakis, N., Baugh, C. M., Benson, A. J., et al. 2011, MNRAS, 410, 53

Fedrigo, G., Cattorini, F., Giacomazzo, B., & Colpi, M. 2024, Phys. Rev. D, 109, 103024

Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9

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

Gerosa, D., Veronesi, B., Lodato, G., & Rosotti, G. 2015, MNRAS, 451, 3941

González, J. A., Hannam, M., Sperhake, U., Brügmann, B., & Husa, S. 2007, Phys. Rev. Lett., 98, 231101

Granato, G. L., De Zotti, G., Silva, L., Bressan, A., & Danese, L. 2004, ApJ, 600, 580

Griffin, A. J., Lacey, C. G., Gonzalez-Perez, V., et al. 2019, MNRAS, 487, 198

Gürlebeck, N. 2015, Phys. Rev. Lett., 114, 151102

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

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

Häring, N., & Rix, H.-W. 2004, ApJ, 604, L89

Heinz, S., Brüggen, M., Young, A., & Levesque, E. 2006, MNRAS, 373, L65

Hopkins, P. F., Hernquist, L., Hayward, C. C., & Narayanan, D. 2012, MNRAS, 425, 1121

Hopkins, P. F., Wetzel, A., Wheeler, C., et al. 2023, MNRAS, 519, 3154

Hu, J. 2008, MNRAS, 386, 2242

Huško, F., Lacey, C. G., Schaye, J., Nobels, F. S. J., & Schaller, M. 2024, MNRAS, 527, 5988

Izquierdo-Villalba, D., Bonoli, S., Dotti, M., et al. 2020, MNRAS, 495, 4681

Junor, W., Biretta, J. A., & Livio, M. 1999, Nature, 401, 891

Kamiaka, S., Benomar, O., Suto, Y., et al. 2019, AJ, 157, 137

Kannan, R., Garaldi, E., Smith, A., et al. 2022, MNRAS, 511, 4005

Katkov, I., Gasymov, D., Kniazev, A., et al. 2024, ApJ, 962, 27

Kimm, T., & Cen, R. 2014, ApJ, 788, 121

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

Kimm, T., Katz, H., Haehnelt, M., et al. 2017, MNRAS, 466, 4826

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

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

Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18

Koudmani, S., Somerville, R. S., Sijacki, D., et al. 2023, MNRAS, submitted [arXiv:2312.08428]

Kraljic, K., Davé, R., & Pichon, C. 2020, MNRAS, 493, 362

Kraljic, K., Duckworth, C., Tojeiro, R., et al. 2021, MNRAS, 504, 4626

Krause, M. G. H., Shabala, S. S., Hardcastle, M. J., et al. 2019, MNRAS, 482, 240

Kurinchi-Vendhan, S., Farcy, M., Hirschmann, M., & Valentino, F. 2023, MNRAS, submitted [arXiv:2310.03083]

Lapiner, S., Dekel, A., & Dubois, Y. 2021, MNRAS, 505, 172

Lescaudron, S., Dubois, Y., Beckmann, R. S., & Volonteri, M. 2023, A&A, 674, A217

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

Lodato, G., & Gerosa, D. 2013, MNRAS, 429, L30

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

Lousto, C. O., & Zlochower, Y. 2011, Phys. Rev. Lett., 107, 231102

Lousto, C. O., & Zlochower, Y. 2013, Phys. Rev. D, 87, 084027

Ma, L., Hopkins, P. F., Kelley, L. Z., & Faucher-Giguère, C.-A. 2023, MNRAS, 519, 5543

Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285

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

Marconi, A., & Hunt, L. K. 2003, ApJ, 589, L21

Martizzi, D., Teyssier, R., & Moore, B. 2013, MNRAS, 432, 1947

McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184

McKernan, B., & Ford, K. E. S. 2023, MNRAS, submitted [arXiv:2309.15213]

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

Merritt, D., & Ferrarese, L. 2001a, MNRAS, 320, L30

Merritt, D., & Ferrarese, L. 2001b, ApJ, 547, 140

Moody, M. S. L., Shi, J.-M., & Stone, J. M. 2019, ApJ, 875, 66

Morganti, R., Oosterloo, T., Oonk, J. B. R., Frieswijk, W., & Tadhunter, C. 2015, A&A, 580, A1

Morsony, B. J., Heinz, S., Brüggen, M., & Ruszkowski, M. 2010, MNRAS, 407, 1277

Morsony, B. J., Miller, J. J., Heinz, S., et al. 2013, MNRAS, 431, 781

Mukherjee, D., Bicknell, G. V., Wagner, A. Y., Sutherland, R. S., & Silk, J. 2018, MNRAS, 479, 5544

Nealon, R., Ragusa, E., Gerosa, D., Rosotti, G., & Barbieri, R. 2022, MNRAS, 509, 5608

Nixon, C. J., King, A. R., & Pringle, J. E. 2011, MNRAS, 417, L66

Nixon, C., King, A., & Price, D. 2013, MNRAS, 434, 1946

O'Dea, C. P., & Baum, S. A. 2023, Galaxies, 11, 67

Ohta, Y., Taruya, A., & Suto, Y. 2005, ApJ, 622, 1118

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

Park, M.-J., Yi, S. K., Dubois, Y., et al. 2019, ApJ, 883, 25

Park, M. J., Yi, S. K., Peirani, S., et al. 2021, ApJS, 254, 2

Peirani, S., Kay, S., & Silk, J. 2008, A&A, 479, 123

Peirani, S., Dubois, Y., Volonteri, M., et al. 2017, MNRAS, 472, 2153

Peirani, S., Sonnenfeld, A., Gavazzi, R., et al. 2019, MNRAS, 483, 4615

Pfister, H., Volonteri, M., Dubois, Y., Dotti, M., & Colpi, M. 2019, MNRAS, 486, 101

Power, C., Navarro, J. F., Jenkins, A., et al. 2003, MNRAS, 338, 14

Prunet, S., Pichon, C., Aubert, D., et al. 2008, ApJS, 178, 179

Reines, A. E., & Volonteri, M. 2015, ApJ, 813, 82

Rennehan, D., Babul, A., Moa, B., & Davé, R. 2023, arXiv e-prints [arXiv:2309.15898]

Rezzolla, L., Barausse, E., Dorband, E. N., et al. 2008, Phys. Rev. D, 78, 044002

Rosen, A., & Bregman, J. N. 1995, ApJ, 440, 634

Saglia, R. P., Opitsch, M., Erwin, P., et al. 2016, ApJ, 818, 47

Sahu, N., Graham, A. W., & Davis, B. L. 2019, ApJ, 876, 155

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

Schawinski, K., Khochfar, S., Kaviraj, S., et al. 2006, Nature, 442, 888

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

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

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

Shapiro, S. L. 2005, ApJ, 620, 59
Sijacki, D., Springel, V., Di Matteo, T., & Hernquist, L. 2007, MNRAS, 380, 877

Silk, J., & Rees, M. J. 1998, A&A, 331, L1

Simmons, B. D., Van Duyne, J., Urry, C. M., et al. 2011, ApJ, 734, 121

Smethurst, R. J., Beckmann, R. S., Simmons, B. D., et al. 2024, MNRAS, 527, 10855

Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776

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

Takaishi, D., Tsukamoto, Y., & Suto, Y. 2020, MNRAS, 492, 5641

Talbot, R. Y., Bourne, M. A., & Sijacki, D. 2021, MNRAS, 504, 3619

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

Teyssier, R., Moore, B., Martizzi, D., Dubois, Y., & Mayer, L. 2011, MNRAS, 414, 195

Thorne, K. S. 1974, ApJ, 191, 507

Toro, E. 1999, Riemann Solvers and Numerical Methods for Fluid Dynamics (Berlin, Heidelberg: Springer-Verlag)

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

Trebitsch, M., Dubois, Y., Volonteri, M., et al. 2021, A&A, 653, A154

Trujillo, I., Ferreras, I., & de La Rosa, I. G. 2011, MNRAS, 415, 3903

Tweed, D., Devriendt, J., Blaizot, J., Colombi, S., & Slyz, A. 2009, A&A, 506, 647

Ubertosi, F., Schellenberger, G., O'Sullivan, E., et al. 2024, ApJ, 961, 134

van den Bosch, R. C. E. 2016, ApJ, 831, 134

Venturi, G., Cresci, G., Marconi, A., et al. 2021, A&A, 648, A17

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518

Volonteri, M., & Ciotti, L. 2013, ApJ, 768, 29

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

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

Volonteri, M., Dubois, Y., Pichon, C., & Devriendt, J. 2016, MNRAS, 460, 2979

Volonteri, M., Pfister, H., Beckmann, R. S., et al. 2020, MNRAS, 498, 2219

Welker, C., Bland-Hawthorn, J., van de Sande, J., et al. 2020, MNRAS, 491, 2864

Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409

Xu, H., Chen, Y., Shi, Y., et al. 2022, MNRAS, 511, 4685

Zhang, B., Lee, K.-G., Krolewski, A., et al. 2023, ApJ, 954, 49

Zhao, Y., Ho, L. C., Shangguan, J., et al. 2021, ApJ, 911, 94

Zheng, X., Zhang, Y., & Röttgering, H. 2024, A&A, 686, A169
## Appendix A: Other individual evolutions

We present in this appendix additional individual evolution of BHs summarized in Table 1. Most of them are displaying a typical $\Psi_{star,BH}$ evolution with different regimes described in Section 3.2.1. This is the case for the BH-146, BH-348 and BH-541 owing mass greater than $10^7 M_{\odot}$ but are contaminated with less than 0.1% of their mass with low resolution DM particles. We reasonably believe that this should not affect the results, in

particular the estimation of $J_{star}$, though such objects are discarded from the statistical study conducted in Section 3. On the contrary, BH-549 and BH-132 show different evolutions. BH-549 is displaying a long phase where its spin is anti-aligned with the stellar angular momentum (see Fig. 9 and the end Section 3.2.1 for more details). As far as BH-132 is concerned, it is off-centered from the galaxy center. The efficiency of both the gas accretion and the Eddington ratio are generally low and the BH spin is most of the time mis-aligned with $J_{star}$.

<figure><img src="image_22.png" /><figcaption>**Fig. A.1.** Cosmic evolution of BH-455 (upper left panel), BH-936 (upper right), BH-549 (lower left panel), and BH-132 (lower right panel). Same as Fig. 4 but for BH-455 and BH-936 display typical $\Psi_{star,BH}$ evolution. BH-549 is displaying a long phase where its spin is anti-aligned with the stellar angular momentum. In this example, the gas accretion disk is in counter-rotation with respect to the stellar component. As far as BH-132 is concerned, it is off-centered from the galaxy center. In this latter case, the efficiency of gas accretion and Eddington ratio are generally low and BH spin is mis-aligned with $J_{star}$.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Fig. A.2.** Cosmic evolution of BH-G648 (upper left panel), BH-146 (upper right), BH-348 (lower left panel), and BH-541 (lower right panel). BH-G648 is a black hole from a GALACTICA zoom. The other BHs are extracted from NEWHORIZON and are hosted by DM halos that have less than 0.1% of their mass composed by low resolution particles (contamination).</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Fig. A.3. Projected distribution of gas and starts (u-g-r band images) at z=0.18 of the different BH-galaxy pairs studied in this Appendix.</figcaption></figure>
**Appendix B: 2D projected angles statistics and galaxy morphologies**

In order to help the comparison with observational analysis, we
derive some statistics on the 2D projected angles depending on
the galaxy morphology. We select from our BH-galaxy sample
at z = 0.18, either galaxies with V/σ < 0.4 (i.e., spheroidal-

dominated galaxies, see the definition at the end of section 3.3)
or V/σ > 0.6 (i.e., disk-dominated galaxies).

We then repeat our Monte-Carlo method and derive the new
distributions of λ with and without specifying the sense of the
BH vector in Figs. B.1 and B.2. Nevertheless, the results are
quite similar to those of presented Fig. 21 without any constraint
on galaxy morphology.

<figure><img src="image_25.png" /><figcaption>**Fig. B.1.** Same as Fig. 21, but selecting galaxies with V/σ > 0.6, i.e., preferentially disk-dominated galaxies (and discard spheroidal-dominated galaxies).</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. B.2. Same as Fig. 21, but selecting galaxies with V/σ < 0.4, i.e., preferentially spheroidal-dominated galaxies.</figcaption></figure>