# Sustained super-Eddington accretion in high-redshift quasars

Alessandro Lupi, Giada Quadri, Marta Volonteri, Monica Colpi, John A Regan

► To cite this version:

Alessandro Lupi, Giada Quadri, Marta Volonteri, Monica Colpi, John A Regan. Sustained super-Eddington accretion in high-redshift quasars. Astronomy & Astrophysics - A&A, 2024, 686, pp.A256. 10.1051/0004-6361/202348788 . hal-04372655

HAL Id: hal-04372655

https://hal.science/hal-04372655v1

Submitted on 31 Oct 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Sustained super-Eddington accretion in high-redshift quasars

Alessandro Lupi<fnref n="1" />,<fnref n="2" />,<fnref n="3" />✉, Giada Quadri<fnref n="3" />✉, Marta Volonteri<fnref n="4" />✉, Monica Colpi<fnref n="3" />,<fnref n="2" />, and John A. Regan<fnref n="5" />

Received 29 November 2023 / Accepted 19 April 2024

## ABSTRACT

Observations of $z \ge 6$ quasars provide information on the early evolution of the most massive black holes (MBHs) and galaxies. Current observations, able to trace both gas and stellar properties, reveal a population of MBHs that is significantly more massive than expected from the local MBH-stellar mass relation. The population lies on, but mostly above, the relation observed in the nearby Universe. This suggests that these objects grew very rapidly. To explain their presence when the Universe was less than 1 Gyr old and to assess the physical conditions for their rapid growth, we explored whether episodes of accretion above the Eddington limit can occur across cosmic epochs. By employing state-of-the-art high-resolution cosmological zoom-in simulations of a $z \sim 7$ quasar, where different accretion regimes are included consistently, together with their associated radiative and kinetic feedback, we show that super-Eddington phases can be sustained for relatively long timescales (tens of millions of years). This allows the MBH to rapidly grow by up to three orders of magnitude, depending on the strength of the kinetic feedback. We also show by means of a semianalytic calculation that the MBH spin remains moderate and does not take on extremely high values during the super-Eddington phases. This results in a lower feedback efficiency, which may allow the rapid growth required to explain over-massive high-redshift MBHs.

**Key words.** black hole physics – galaxies: evolution – galaxies: formation – galaxies: high-redshift – quasars: general – quasars: supermassive black holes

## 1. Introduction

Massive black holes (MBHs) are ubiquitous in the Universe and inhabit the centre of massive galaxies up to redshift $z \ge 6$ (e.g. Fan et al. 2006, 2023; Mortlock et al. 2011; Bañados et al. 2018; Maiolino et al. 2023), with masses up to $10^{9-10} M_\odot$. Observationally, they are commonly identified via gas accretion through the conversion of gravitational energy into radiation, which makes them shine as active galactic nuclei (AGN). They sometimes also produce powerful collimated jets. According to Soltan's argument (Soltan 1982), the evolution of the AGN luminosity function, and the local MBH-mass density (Marconi et al. 2004), MBHs gain most of their mass via radiatively efficient accretion. Hence, they should have formed from lower-mass black hole “seeds” (see, e.g. Inayoshi et al. 2020; Volonteri et al. 2021 for a review).

In this context, the observations of high-redshift MBHs help placing tight constraints on the minimum initial seed mass, which has to be about $M_{\text{seed}} \ge 10^4 M_\odot$ when we assume growth to occur via radiatively efficient accretion with $\eta_{\text{rad}} \sim 0.1$ at the Eddington limit. Several groups have studied the formation mechanisms of such “heavy” seeds that facilitate interpreting current observations. However, the formation of the most massive seeds requires rare peculiar conditions (e.g. Latif et al. 2015, 2022; Schauer et al. 2017; Lupi et al. 2021). An alternative possibility, which was also suggested by recent results (e.g. Davies et al. 2019; Yang et al. 2023), is that these MBHs grew rapidly because of intermittent phases of super-Eddington accretion

from lighter seeds (Madau et al. 2014; Volonteri et al. 2015; Lupi et al. 2016; Pezzulli et al. 2016).

The plausibility of super-Eddington accretion has been demonstrated to occur in the local Universe in tidal disruption events (Lin et al. 2017) and ultra-luminous X-ray sources (Bachetti et al. 2014) and was initially proposed in the context of non-spherical accretion flows within the “slim-disc” model (Abramowicz et al. 1988). Recent simulations of accretion discs showed that the slim-disc solution, in which the radiation is trapped in the innermost regions of the accretion disc and is advected inwards within the fluid, making the disc moderately luminous, is instead characterised by powerful radiatively driven outflows and jets that escape through a central funnel perpendicular to the disc itself (Sądowski & Narayan 2016; Sądowski et al. 2016; Jiang et al. 2019).

While most of the studies focussing on accretion above the Eddington rate were performed on the typical scales of the accretion disc, large-scale simulations have commonly neglected this regime and based their results on the common Eddington-limited accretion. Moreover, the inability to resolve the influence radius of the MBH and, in many cases, the low resolution preventing the dense gas from being properly resolved also required additional tuning parameters, such as the $\alpha$ boost factor on the accretion rate (Springel et al. 2005; Di Matteo et al. 2005) and the coupling efficiency parameter for the accretion-powered feedback (dumped in the form of thermal energy), which was commonly set to values in the range 0.05–0.15, depending on the numerical technique employed to reproduce

<footnote n="1">¹ Dipartimento di Scienza e Alta Tecnologia, Università degli Studi dell'Insubria, Via Valleggio 11, 22100 Como, Italy
e-mail: alessandro.lupi@uninsubria.it</footnote>

<footnote n="2">² INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy</footnote>

<footnote n="3">³ Dipartimento di Fisica “G. Occhialini”, Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy</footnote>

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

<footnote n="5">⁵ Centre for Astrophysics and Space Sciences Maynooth, Department of Theoretical Physics, Maynooth University, Maynooth, Ireland</footnote>
local observations (Di Matteo et al. 2005; Dubois et al. 2014). Even though improved models for accretion and feedback were developed more recently, for instance correcting for the gas angular momentum (Tremmel et al. 2017) or dumping the MBH feedback as radiation-driven winds (Choi et al. 2012; Anglés-Alcázar et al. 2017), the physics of the accretion process above the Eddington limit was almost never considered.

In the last few years, however, several groups have started to account for the possibility of super-Eddington accretion, but often in a very simplistic way. For instance, some works only extended the standard accretion and feedback prescriptions above the Eddington limit(e.g. Zhu et al. 2022; Ni et al. 2022; Bhowmick et al. 2022), whereas others only considered the radiation trapping and the subsequent decrease in radiative efficiency (Lupi et al. 2016; Rennehan et al. 2023). A more physically motivated modelling that accounts for the impact of kinetic winds/jets during the super-Eddington phases on the interstellar medium of a galaxy host and on the MBH growth itself has instead only been considered in idealised setups such as in isolated galaxies (Massonneau et al. 2023a), circum-nuclear discs (Sassano et al. 2023), and atomic cooling haloes just after the formation of a heavy-seed MBH (Regan et al. 2019).

In this work, we move forward and study the MBH evolution in the quasar host originally studied in Lupi et al. (2019, 2022). We properly account for super-Eddington accretion phases in a full cosmological context. The simulation adopts the same initial conditions of the high-redshift quasar originally studied in Lupi et al. (2019, 2022, hereafter Papers I and II). This is Paper IV of a series of papers addressing properties of high-redshift quasar hosts and their MBHs. In Paper I, we presented and discussed the main evolution of the target galaxy and its central MBH, focussing on the stellar and gas tracers (total gas and [CII] emission), and found that super-Eddington phases were measured in the simulation, even though accretion was capped at the Eddington limit. In Paper II, we extended the analysis by focussing on the dynamics and morphology of the main galaxy as a function of redshift. In Paper III (Lupi et al., in prep.) we will discuss the evolution of the entire MBH population that forms during the simulation, and in Paper V (Quadri et al., in prep.) we will focus in detail on the impact of the super-Eddington regime on the galaxy host and on the properties of quasar outflows.

The paper is organised as follows. In Sect. 2, we recapitulate the setup of the simulation and describe the improvements relative to the previous works. In Sect. 3, we present our results. In Sect. 4, we draw our conclusions.

## 2. Numerical setup

The simulation followed the evolution of a massive halo ($M_{\text{halo}} \sim 3 \times 10^{12} M_{\odot}$ at $z = 6$) that is expected to represent a quasar host (Di Matteo et al. 2017). The initial conditions were accurately created to match the expected halo mass (Di Matteo et al. 2017; Tenneti et al. 2018) and the galaxy overdensity significance (Uchiyama et al. 2018; Mignoli et al. 2020) via MUSIC (Hahn & Abel 2013), adopting the Planck Collaboration XIII (2016) cosmological parameters, with $\Omega_{\text{m}} = 0.3089$, $\Omega_{\Lambda} = 0.6911$, $\Omega_{\text{b}} = 0.0489$, $\sigma_{8} = 0.8159$, $n_{s} = 0.9667$, and $H_{0} = 67.74 \text{ km s}^{-1} \text{ Mpc}^{-1}$, without a contribution from radiation and curvature. From a parent dark-matter-only simulation, we recursively zoomed-in on a Lagrangian volume extending up to 2.5 virial radii of the target halo, following the approach

by Fiacconi et al. (2017) to exclude any contamination by low-resolution dark matter particles within the virial radius.

The simulation was run with GIZMO (Hopkins 2015), a descendant of GADGET3 (Springel et al. 2008) and GADGET2 (Springel 2005) in its meshless finite-mass mode. The spatial resolution of the simulation was set to 40, 10, and $2.5 \text{ pc h}^{-1}$ for dark matter, stars, and MBHs, respectively, whereas fully adaptive softening was assumed for the gas component, down to a minimum of ~5 pc. The mass resolution was $\sim 10^4 M_{\odot}$ for baryons and $\sim 10^5 M_{\odot}$ for dark matter.

### 2.1. Baryonic physics

Our simulation was performed with state-of-the-art sub-grid prescriptions that allowed us to follow non-equilibrium chemistry of primordial species, star formation, and stellar feedback in detail, as well as MBH seeding, accretion, and feedback. Compared to the previous works of the series, we here slightly revised and improved many of the sub-resolution prescriptions, which we describe below.

- Gas thermodynamics and chemistry: We further extended our chemical network to include high-ionisation states of several important species that are commonly observed in quasar hosts, namely $C[I-IV]$, $O[I-VI]$, $N[I-V]$, and Fe[I-II], also accounting for their contribution to the low-temperature cooling of the gas. Since MBHs are commonly surrounded by a hot corona emitting in X-rays, we also incorporated detailed X-ray chemistry calculations in our network, accounting for the impact of Compton heating by the AGN (assuming $T_{\text{Compton}} = 3.23 \times 10^6 \text{ K}$ and $T_{\text{Compton}} = 8.41 \times 10^7 \text{ K}$ for soft and hard X-rays, respectively). This will be discussed in Lupi et al. (in prep.) in greater detail.

- Star formation: We slightly revised our estimate for the turbulent support of the gas as in Hopkins et al. (2013), accounting for the particle distribution inside the kernel, which gives $\sigma_{\text{turb}} = \|\nabla v\|/5$. Relative to Lupi et al. (2019, L19 hereon), we updated our star formation efficiency employing the Padoan et al. (2012) model, as in Lupi & Bovino (2020).

- Stellar mechanical feedback: We redetermined the scalings in Martizzi et al. (2015) to improve the agreement with their results, that is, we properly accounted for the initial fraction of kinetic and thermal energy during the Sedov-Taylor phase $f_{\text{kin}} \sim 0.28$.

- Stellar radiative feedback: Instead of the cost-effective approximated radiation transport of L19, we included on-the-fly radiation transport as in Lupi et al. (2020), with the reduced speed of light set to $c_{\text{red}} = 1000 \text{ km s}^{-1}$, which is high enough compared to the gas motions to ensure consistent results. In addition, we followed X-ray chemistry, and we therefore now follow radiation in 11 photobins ranging from 0.7 keV up to 10 keV, where 2 bins are used to cover soft (0.2-2 keV) and hard (2-10 keV) X-rays.

### 2.2. MBH accretion/feedback and dynamics

In addition to the changes above, we devised a novel set of prescriptions for MBH growth and dynamics that we discuss in detail in the following. The MBH seeding is instead identical to that in L19 and occurs in galaxies with a stellar mass $> 10^8 M_{\odot}$ that do not yet host an MBH, that is, galaxies that are identified through an on-the-fly Friends-of-Friends algorithm (see L19 for details).
### 2.2.1. MBH dynamics

In most cosmological simulations, the mass and spatial resolution of seed MBHs is not sufficient to accurately resolve the dynamical friction bringing MBHs to the centre of galaxies, and often also the interaction with other particles, leading to spurious scattering of the MBHs. For this reason, most simulations include an ad hoc MBH pinning procedure (Di Matteo et al. 2005; Schaye et al. 2015; Barai et al. 2018) that moves the MBH to the potential minimum inside its kernel at every time step. Although effective, this procedure can produce unphysical behaviour such as superluminal motions or artificial suppression of the MBH wandering. A more physically motivated prescription artificially corrects the dynamics accounting for the unresolved dynamical friction effect (Dubois et al. 2013; Tremmel et al. 2015), as long as the mass ratio of the MBH relative to the other tracers (gas, stars, and dark matter) is high enough (Tremmel et al. 2015; Pfister et al. 2019). Lupi et al. (2019) ensured a reasonable dynamical evolution by seeding the MBH with an already high mass of $M_{BH} = 10^6 M_{\odot}$. In this work, we instead opted for decoupling the MBH mass into a physical mass (used for accretion) and a dynamical mass (for the dynamics; Anglés-Alcázar et al. 2017), which evolve together as soon as the first reaches the second. The physical mass was set to $M_{BH} = 10^5 M_{\odot}$, and the dynamical mass $M_{dyn,MBH}$ to a value ten times higher. This ensured a better dynamical evolution. From a physical point of view, this initially higher dynamical MBH mass can be considered as an unresolved stellar envelope (a nuclear stellar cluster) surrounding the MBH, which is commonly found in many galaxies in the local Universe (see Neumayer et al. 2020, for a review). From a numerical point of view, however, it was simply used to avoid the spurious scattering by other particles, which despite the very high resolution adopted in this simulation are still more massive than individual stars.

Unlike Tremmel et al. (2017) and Ma et al. (2021), the unresolved dynamical friction in this work was implemented following the more accurate model by Pfister et al. (2019), where the dark matter and stellar distributions are accounted for separately (thus allowing for different distribution functions), and the high-velocity part of the distribution function was also considered. Moreover, we also included gas-driven dynamical friction as described in Tanaka & Haiman (2009), following the prescription in Escala et al. (2004).

### 2.2.2. MBH accretion and feedback

One of the main novelties of this work is the inclusion of three accretion regimes covering the entire range of accretion rates. We also account for the different feedback mechanisms relevant in each of them<fnref n="1" />. We classify them in terms of the Eddington ratio $\lambda = \dot{M}_{BH}/\dot{M}_{Edd}$, where $\dot{M}_{BH}$ is the accretion rate, $\dot{M}_{Edd} = 16 L_{Edd}/c^2$ (Madau et al. 2014), $L_{Edd}$ is the Eddington luminosity, and c is the speed of light. In principle, the radiative efficiency depends on the MBH spin, which in the current simulation is not evolved over time. We set the spin magnitude to a constant value of a = 0.7, which gives the commonly adopted radiative efficiency in the standard sub-Eddington radiatively efficient accretion regime ($\eta_{rad} = 0.103$), and we further assumed that the spin is parallel to the angular momentum of the gas within the MBH kernel $l_{gas,MBH}$. While accretion is still mod-

elled using the Bondi-Hoyle-Lytteton prescription, the resulting MBH feedback was implemented in radiative and kinetic forms.

Radiative feedback was injected as in Lupi et al. (2020), assuming a composite black body plus X-ray corona spectrum for the MBH, with a bolometric luminosity defined as $L_{bol} = \eta_{rad} \dot{M}_{BH}c^2$, with $\eta_{rad}$ depending on the accretion regime (see below). The fraction of $L_{bol}$ associated with the X-ray corona was determined according to Duras et al. (2020) as

$$f_{2-10 \text{ keV}}^{-1} = \frac{L_{\text{bol}}}{L_{2-10 \text{ keV}}} = 12.76 \left[ 1 + \left( \frac{\log(L_{\text{bol}}/L_{\odot})}{12.15} \right)^{18.78} \right]. \quad (1)$$

From this expression, we determined the soft X-ray fraction $f_{0.2-2\text{keV}}$ assuming a power-law spectrum with a slope -1.7 for the corona (Regan et al. 2019), and we finally obtained the residual black-body component as

$$f_{\text{BB}} = 1 - f_{0.2-2\text{keV}} - f_{2-10\text{keV}}. \quad (2)$$

Kinetic feedback was instead implemented for the different regimes as<fnref n="2" />

- $\lambda < 2.5 \times 10^{-3}$ (ADAF regime, Yuan & Narayan 2014): In this regime, the disc is optically thin and geometrically thick due to inefficient cooling. Ions and electrons decouple, resulting in a two-temperature accretion flow where the radiative efficiency decays roughly as $\eta_{rad} \propto \lambda^{2/3}$ (Xie & Yuan 2012). In this regime, a jet is launched along $l_{gas,MBH}$ and forms a cylinder with the base defined by the MBH kernel with an efficiency determined as (Blandford & Znajek 1977; Tchekhovskoy & Giannios 2015)

$$\eta_{\text{jet}} \approx 2.5 \frac{a}{(1 + \sqrt{1 - a^2})^2} \left( \frac{\phi}{\phi_{\text{MAD}}} \right)^2 \approx 0.417 \Phi^2, \quad (3)$$

where $\phi$ is the maximum magnetic flux in the disc, $\phi_{MAD}$ is a critical value corresponding to a magnetically arrested disc (MAD, Narayan et al. 2003), and $\Phi = \phi/\phi_{MAD}$. The mass loading of the jet $\beta_{jet} \equiv \dot{M}_{jet}/\dot{M}_{BH}$ in our simulation was determined assuming energy conservation and a jet velocity $v_w = 0.1c = 3 \times 10^4 \text{ km s}^{-1}$, which gives

$$\beta_{\text{jet}} = \frac{\dot{M}_{\text{jet}}}{\dot{M}_{\text{BH}}} = 2\eta_{\text{jet}} \left(\frac{c}{v_w}\right)^2 \approx 83.4\Phi^2. \quad (4)$$

- $2.5 \times 10^{-3} < \lambda < 1$ (sub-Eddington radiatively efficient accretion): This regime occurs in typical AGN, where the disc is geometrically thin and optically thick, and it can be understood in terms of the Shakura and Sunyaev solution (Shakura & Sunyaev 1973). In this case, we assumed the MBH feedback to be in the form of radiation (again assuming a composite black-body plus corona spectrum) and bipolar line-driven winds whose mass loading is given by momentum conservation during the matter-radiation interaction (Choi et al. 2012; Anglés-Alcázar et al. 2017), that is,

$$\beta_w = \frac{L_{BH}}{\dot{M}_{BH} v_w c} = \eta_{\text{rad}} \frac{c}{v_w} \approx 1, \quad (5)$$

where $L_{BH} = \eta_{rad} \dot{M}_{BH}c^2$ is the accretion luminosity. The resulting energy coupling efficiency for the wind, which was launched assuming a semi-aperture of 45 deg (Sala et al. 2021), can be estimated in

$$\epsilon = \frac{\dot{M}_w v_w^2}{2 L_{BH}} = \frac{\dot{M}_w v_w^2}{2 \eta_{\text{rad}} \dot{M}_{BH} c^2} = \frac{\beta v_w^2}{2 \eta_{\text{rad}} c^2} = \frac{v_w}{2c} = 0.05. \quad (6)$$

<footnote n="1">¹ During the drafting of this work, Rennehan et al. (2023) also proposed inclusion of the three MBH accretion regimes in simulations, but not the impact of jets during the super-Eddington stage.</footnote>

<footnote n="2">² The final normalisations of the equations in this section correspond to the spin magnitude assumed in our simulation, $a = 0.7$.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Evolution of the MBH accretion rate in the HMAD (left panel) and MAD (right panel) simulations. The red dots correspond to the instantaneous accretion rate, i.e. at every time step of the simulations, whereas the solid green line corresponds to the accretion rate averaged over 1 Myr. The dashed blue line represents the instantaneous Eddington accretion rate at the MBH mass.</figcaption></figure>

– $\lambda > 1$ (super-Eddington accretion): In this regime, we assumed radiative plus kinetic MBH feedback, where the radiative component was determined according to the slim-disc solution, with the radiative efficiency obtained by Lupi et al. (2016),

$$ \eta_{\text{rad}} = \frac{A(a)}{16\lambda} \left[ \frac{0.985}{1/\lambda + B(a)} + \frac{0.015}{1/\lambda + C(a)} \right], \quad (7) $$

where $A(a)$, $B(a)$, and $C(a)$ are spin-dependent coefficients that in our case assumed values 1.915, 0.795, and 0.017 respectively. Although radiation is expected to be trapped in the innermost regions of the slim disc, we neglected any change in the spectral shape of the radiation spectrum, leaving this exploration for a future study. The kinetic efficiency in this case was computed according to Tchekhovskoy & Giannios (2015) and Sdowski et al. (2016) for a magnetically driven jet, that is,

$$ \eta_{\text{jet}} = 1.3a^2\Phi^2 \approx 0.637\Phi^2, \quad (8) $$

corresponding to a mass loading $\beta_{jet} \approx 127.4\Phi^2$.

A crucial parameter affecting the jet efficiency is the limiting magnetic field relative to the MAD limit. We considered two cases, one case at the MAD limit ($\Phi = 1$) and the other case at half the limit ($\Phi = 0.5$; HMAD hereon; see Sdowski et al. 2016, for a discussion). At every accretion event, we estimated the accretion rate on the MBH, $\dot{M}_{BH}$, from its 96 nearest gas neighbours, compared it with the Eddington limit, and then determined $\eta_{rad}$ and the value of $\beta$ for the corresponding regime.

An important aspect that must be considered is that for very high accretion rates, the large mass-loading factor might yield a total gas mass affected by the MBH accretion/feedback process over a time step $\Delta t$ that exceeds the available mass in the MBH kernel. In these cases, we followed Regan et al. (2019) and assumed that the estimated accretion only occurred for a fraction of time,

$$ f_{\text{acc}} = \min \left\{ 1, \frac{M_{\text{ngbs}}}{(1 + \beta)\Delta M_{\text{BH}}} \right\}, \quad (9) $$

where $\Delta M_{BH} = \dot{M}_{BH}\Delta t$ and $\beta$ is the mass-loading factor of the kinetic wind/jet. At this point, we randomly flagged enough gas

particles around the MBH to guarantee $\Delta M = (1 + \beta)\dot{M}_{BH}\Delta t$, where part of the mass is accreted and part is kicked away in a kinetic wind/jet. In the case of $M_{dyn,MBH} > M_{MBH}$, we simply selected particles to be kicked away in an outflow, in order to ensure mass conservation. In order to prevent the outflowing gas particles from propagating over large distances without interacting (over a single time step they might end up well outside the galaxy), we updated the time step $\Delta t$ of the kicked particles to $\Delta t' = \min\{\Delta t, C_{CFL}\Delta x/v_{w}\}$, where $\Delta x$ is the effective gas cell size (Hopkins 2015), and $C_{CFL}$ is the Courant factor.

## 3. Results

With the model just described, we ran the same initial conditions of Lupi et al. (2019) down to $z \sim 9.7$, when the MBH has already formed near the centre of its host, but the galaxy mass is still low enough for supernovae to stunt the MBH growth completely. At this point, we split the simulation into two equivalent runs, MAD and HMAD, which only differ by the magnetic flux parameter $\Phi$. The simulations were run down to $z \sim 7.5$ in order to follow the early growth of the MBH seeds in the galaxy host, whose evolution we discuss in this section.

In Fig. 1, we report the accretion rate on the MBH in the two cases, as directly obtained from the simulation (red dots) and averaged over a timescale of 1 Myr (solid green line). As a comparison, we also report as the dashed blue line the Eddington accretion rate at every step of the evolution. The instantaneous $\dot{M}_{BH}$ exhibits very large excursions, which are associated with the intermittency between large inflows and the MBH feedback self-regulation.

The evolution can be divided into three main stages on average, however (see also Lupi et al. 2019). (i) In the early stages after MBH seeding ($t \lesssim 530$ Myr), supernova feedback strongly perturbs the gas in the galaxy, preventing the MBH (which is also offset from the galaxy centre) from efficiently accreting. (ii) When the MBH settles in the centre of the galaxy and the potential well becomes deep enough to sustain large gas inflows to the centre ($M_{\star} \approx 10^{10} M_{\odot}$), $\dot{M}_{BH}$ rapidly grows, easily exceeding Eddington by a factor of a few (MAD) or a few tens (HMAD). At this stage, the HMAD simulation shows an almost unimpeded super-Eddington growth for about 60 Myr, followed by (iii) a decrease in the accretion rate to values between a fraction and
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Evolution of the MBH mass in the MAD (shown as a solid orange line) and HMAD (shown as a dashed blue line) simulations. For comparison, we also show two dot-dashed black lines corresponding to the MBH growth for a constant Eddington accretion. The thick line starts at *t* = 500 Myr (consistent with the HMAD case), and the thin line starts at *t* = 350 Myr, immediately after seeding.</figcaption></figure>

100% Eddington, corresponding to the self-regulation stage. In both cases, there are sudden drops in the accretion rates at around *z* = 8.3 and *z* = 9 for the HMAD and MAD cases, respectively. These rapid variations are associated with the potential misalignment of the gas in the galaxy nuclear region with respect to the galactic disc plane, which results in winds/jets that directly hit the disc instead of escaping perpendicular to it. In these phases, the gas is expelled from the galaxy centre and falls back on a timescale comparable with the free-fall time of the gas, that is, 5 Myr $\lesssim t_{ff} \lesssim 50$ Myr for a typical gas density between 1 and $100\ \text{cm}^{-3}$.

The evolution just discussed reflects the mass growth reported in Fig. 2 for both simulations. The three stages just discussed can be noticed also in this figure. The time at which the MBH starts to grow is slightly different between the two simulations because of the small differences that appeared when the simulations were restarted from the parent simulation, which built up over time and changed the MBH history. Except for this and the plateaus in the growth history when the accretion rate drops, the growth exhibits a similar behaviour, that is, an initial rapid growth well above the Eddington limit, which lasts longer in the HMAD case due to the lower MBH feedback efficiency, followed by a self-regulation phase when the MBH settles at the Eddington limit or at a fraction of it.

Despite the quite high efficiency of the jet feedback during super-Eddington phases, the HMAD MBH is able to grow by almost three orders of magnitude in less than 100 Myr, whereas the MAD MBH stays one order of magnitude below throughout the same time interval. Another important aspect to keep in mind here is that when the MBH in the HMAD case starts to self-regulate its growth, the MBH mass is so high that super-Eddington phases begin to require extremely high inflow rates in the galaxy nucleus, which are increasingly less likely as the galaxy evolves and the gas fraction diminishes. In the Eddington-limited case, an initially more massive MBH or an earlier start of accretion are needed to ensure a growth in mass that is comparable to the two cases we simulated here. In Fig. 3, we report the gas density profile from the two simulations (HMAD on the left and MAD on the right) computed in spherical

shells at four different times: *z* = 9.6 represents the initial growth phases, when the accretion rate in HMAD starts to rise towards the super-Eddington regime; *z* = 8.8 is a super-Eddington phase in both simulations; *z* = 8.2 represents a feedback-dominated phase, when the MBH in MAD slightly suppresses its growth, settling around the Eddington limit, whereas the MBH in HMAD launches a wind/jet through the disc that is able to carve a cavity of ~300 pc (as discussed above); and *z* = 7.8, when both simulation settle around/below the Eddington limit, starting a self-regulated growth phase. Although the two simulations are formally independent, the gas distribution is in general very similar, but for the phases in which MBH feedback succeeds at expelling gas from the centre. This is consistent with most recent simulations that suggested that the central MBH does not directly affect the entire galaxy host, whose evolution is mainly determined (except for a few transient phases) by the cosmological environment.

Since this rapid growth may have important implications for the MBH-galaxy correlations, we report in Fig. 4 the evolution of the MBH together with the galaxy host stellar mass for the two cases, using the same colour scheme as in Fig. 2. We also show for comparison the observations of low-redshift AGN by Reines & Volonteri (2015) as grey dots, elliptical galaxies by Kormendy & Ho (2013) as grey stars, and high-redshift (only at *z* ≥ 6) observations by ALMA as green squares (Neeleman et al. 2021) and JWST as red diamonds (Yue et al. 2024; Stone et al. 2023, 2024; Harikane et al. 2023; Maiolino et al. 2023). We also report as magenta crosses the “little red dots” observed by Greene et al. (2024), with the stellar masses estimated according to JWST photometry alone and including ALMA photometry by Labbe et al. (2023), where the magenta shaded area connects the two mass estimates for each source<fnref n="3" />. Finally, we also report as a cyan hexagon GN-z11 (Maiolino et al. 2024), the highest-redshift MBH to date, which is thought to be accreting at super-Eddington rates, and the results by L19 as a dot-dashed black line. The overlap between GN-z11 and the L19 simulation is simply due to the higher seeding mass ($10^6 M_{\odot}$) in L19 and not to the physical evolution in the simulation. In both cases, the MBH starts at the lower end of the distribution and grows rapidly, but only the HMAD case is able to reach the upper limit of the distribution because of a less effective MBH feedback. Then, when the MBH starts to self-regulate, even the growth of the HMAD MBH slows down, moving toward the region occupied by giant ellipticals (see L19 for a discussion).

Finally, we note that the super-Eddington phase was able to start only after the weakening of the impact of supernova feedback on the galaxy, in agreement with other simulations (Dubois et al. 2014; Anglés-Alcázar et al. 2017; Lupi et al. 2019). This prevented the simulated MBHs from reaching the region occupied by high-redshift over-massive MBHs, especially those observed in low-mass galaxies (e.g. Maiolino et al. 2023). This might be a potential issue in explaining the overmassive MBHs, unless an even weaker super-Eddington feedback or more favourable inflow conditions occur around these MBHs.

## 4. Discussion and conclusions

We have performed high-resolution cosmological zoom-in simulations of a massive halo at high redshift that is expected to be a quasar host. The simulations include for the first time in

<footnote n="3">³ The emission in this case is dominated by the AGN and dust, hence, the stellar mass estimates are only indicative, but likely lower than $10^{10} M_{\odot}$ because the little red dots are selected to be point sources.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Gas density profiles for HMAD (left panel) and MAD (right panel) in spherical shells at different times: $z = 9.6$ (sub-Eddington, solid blue line), $z = 8.8$ (super-Eddington, dashed orange line), $z = 8.2$ (not accreting in HMAD and mildly sub-Eddington in MAD, dotted green line), and $z = 7.8$ (self-regulated, dot-dashed red line).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Evolution of the MBH mass in the MAD (shown as a solid orange line) and HMAD (shown as a dashed blue line) simulations relative to the galaxy stellar mass compared with the local AGN observations by Reines & Volonteri (2015, grey dots), local ellipticals by Kormendy & Ho (2013, grey stars), and high-redshift ($z \ge 6$) observations by ALMA (Neeleman et al. 2021, green squares) and JWST (Yue et al. 2024; Stone et al. 2023, 2024; Harikane et al. 2023; Maiolino et al. 2023, red diamonds). We also show the “little red dots” reported by Greene et al. (2024) as magenta crosses, with the stellar masses estimated with or without ALMA photometry, respectively, connected by magenta shaded areas (see Labbe et al. 2023, for details), and the $z \sim 10$ MBH GN-z11 (Maiolino et al. 2024), thought to be super-Eddington, as a cyan hexagon. For completeness, we also show the results of the simulation in L19 as a dot-dashed black line.</figcaption></figure>

a cosmological simulation a sub-resolution model for the three main accretion regimes, from ADAF up to super-Eddington. The model accounts for MBH accretion-powered feedback in the form of radiation and kinetic winds/jets, with an efficiency that is determined consistently for the different cases. The results we found here will be complemented by a more thorough analysis of the MBH-galaxy interplay and the properties of the MBH outflows in a forthcoming paper (Quadri et al., in prep.). They show that super-Eddington phases might be sustained over timescales

of a few tens of millions of years in massive systems where large gas inflows are frequent, systems that are quite different from those investigated in previous works. However, some important caveats are listed below.

- The resolution of the simulation is not high enough to properly resolve the Bondi-Hoyle-Lyttleton radius of the MBH. This might affect the estimated accretion rate. Nonetheless, a higher accretion rate is reflected in a stronger feedback, which has the potential of clearing out the central region of the galaxy (as sometimes occurred in our simulations). This would hinder accretion for a few tens of millions of years.

- The launch direction of the kinetic feedback is instantaneously aligned with the gas within the MBH kernel, without taking into account the actual alignment of an unresolved accretion disc. This means that if the gas in the galaxy has settled into a disc-like configuration, the kinetic feedback escapes most of the time without interacting significantly with the galaxy, which means that accretion is not suppressed for a long time. In a more realistic case, the alignment of the accretion disc might occur on much longer timescales and might in this way affect the impact of the feedback.

- The results of this simulation provide a more optimistic outlook on the role of super-Eddington accretion in growing high-redshift black holes compared to previous numerical investigations in a galactic/cosmological context (Johnson et al. 2011; Regan et al. 2019; Massonneau et al. 2023a). In particular, even though Johnson et al. (2011) studied the early growth of MBH seeds in highly overdense regions, finding that radiation is able to regulate MBH accretion, their exploration was limited to the very first stages, that is, a pristine atomic cooling halo that grew little mass over the 1 Myr timescale they considered. In our case, the environment in which the MBH is embedded is completely different, as it is a metal-enriched massive galaxy with massive inflows. In metal-rich conditions and at high densities, gas cooling becomes very efficient and thus suppresses the effect of accretion-powered photoheating (in agreement with the ineffective super-Eddington growth found by Regan et al. 2019, in similar physical conditions). To assess the reasons for the difference with Regan et al. (2019) and Massonneau et al. (2023a), we measured the gas density and
<figure><img src="image_6.png" /><figcaption>Fig. 5. Evolution of the MBH spin (solid blue line) and MAD-ness (dashed orange line) in the HMAD (left panel) and MAD (right panel) simulations, as modelled in Ricarte et al. (2023).</figcaption></figure>

temperature around the MBH, finding that super-Eddington bursts were followed by an increase in the gas temperature up to $10^6 < T/K < 10^9$ (above the peak of the cooling curve) due to the jet shocking with the gas. This is consistent with the results in Regan et al. (2019) and Massonneau et al. (2023a). Despite this strong heating, the density around the MBH did not change dramatically and always remained in the range $n_H \approx 10^{4-5} \text{ cm}^{-3}$ on average. This value is significantly higher than that in Massonneau et al. (2023a). This suggests that the combination of a deep potential well, the pressure of the inflowing gas from larger scales, and a short free-fall time (see Regan et al. 2019 for a discussion) prevented the gas from escaping from the galaxy nucleus. As a consequence, after each jet event, the accretion rate suppression we observed was short-lived and was followed by a rapid rise as soon as the gas cooled down. However, the limited resolution of our simulation might still underestimate the effect of the jet feedback by not properly resolving the early expansion of the shocked gas, and further investigations are needed.

- The magnitude of the spin in the simulations was kept fixed. This has profound implications for the evolution. First, since the jet launched during super-Eddington accretion phases extracts rotational energy from the MBH (Blandford & Znajek 1977), the spin should decrease rapidly, thus overcoming the spin-up due to coherent accretion. This would likely result in a lower spin during the first accretion phases and hence, in a lower kinetic efficiency for the jet, even in the MAD case, which would enable the MBH to grow almost unimpeded until the MBH spin and mass become high enough for feedback to start playing a role (Massonneau et al. 2023b). In order to test this idea, we computed the spin evolution semi-analytically by means of the model in Ricarte et al. (2023), employing the accretion rates and the MBH masses from our simulations. In particular, we employed the 1 Myr average quantities instead of the full data sample in order to limit the effect of strong fluctuations on our conclusions. The results are reported in Fig. 5, where we show the evolution of the MBH spin and of the disc MAD-ness as defined in Ricarte et al. (2023) for $f_{Edd} > 0.03$ (at lower Eddington ratios, a MAD disc is assumed),

$$ \frac{\phi}{\phi_{\text{MAD}}} = \frac{(f_{\text{Edd}}/f_c)^\alpha}{1 + (f_{\text{Edd}}/f_c)^\alpha} \quad (10) $$

where $f_c = 1.88$ and $\alpha = 1.29$.

We clearly observe that during the super-Eddington phases, the MBH spin decreases to values around 0.5–0.6, which correspond to jet efficiencies that are lower by 1.5–2 times than those employed in the simulations. Only during the later phases, when the accretion rate settles at sub-Eddington values, the spin is able to grow beyond the initial value. During the evolution, the MAD-ness of the disc also exhibits strong variations, with average values during the super-Eddington phases of about 0.6 (HMAD) and 0.8 (MAD), which is in between the values assumed in our simulations. Combining these effects, we obtain a net efficiency of the super-Eddington jets of $\eta_{\text{jet}} \propto a^2(\phi/\phi_{\text{MAD}})^2$, which is more in line with our HMAD simulation than with the MAD case.

To conclude, we have shown that despite the unfavourable conditions associated with assuming a relatively high spin and a consequently strong jet feedback, MBHs in gas-rich environments at high redshift can support long-lasting super-Eddington accretion phases and grow rapidly in mass within their host galaxies. Moreover, because of this rapid growth, these MBHs can easily move above the local correlations before they start to regulate themselves. This might explain the formation of overmassive systems.

*Acknowledgements.* A.L. and M.C. acknowledge support by the PRIN MUR “2022935STW”. J.R. acknowledges support from the Royal Society and Science Foundation Ireland under grant number URF/R1/191132. J.R. also acknowledges support from the Irish Research Council Laureate programme under grant number IRCLA/2022/1165.

## References

* Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, ApJ, 332, 646

* Anglés-Alcázar, D., Faucher-Giguère, C.-A., Quataert, E., et al. 2017, MNRAS, 472, L109

* Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473

* Bachetti, M., Harrison, F. A., Walton, D. J., et al. 2014, Nature, 514, 202

* Barai, P., Gallerani, S., Pallottini, A., et al. 2018, MNRAS, 473, 4003

* Bhowmick, A. K., Blecha, L., Ni, Y., et al. 2022, MNRAS, 516, 138

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

* Choi, E., Ostriker, J. P., Naab, T., & Johansson, P. H. 2012, ApJ, 754, 125

* Davies, F. B., Hennawi, J. F., & Eilers, A.-C. 2019, ApJ, 884, L19

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

* Di Matteo, T., Croft, R. A. C., Feng, Y., Waters, D., & Wilkins, S. 2017, MNRAS, 467, 4243

* Dubois, Y., Gavazzi, R., Peirani, S., & Silk, J. 2013, MNRAS, 433, 3297

* Dubois, Y., Pichon, C., Welker, C., et al. 2014, MNRAS, 444, 1453
Duras, F., Bongiorno, A., Ricci, F., et al. 2020, A&A, 636, A73

Escala, A., Larson, R. B., Coppi, P. S., & Mardones, D. 2004, ApJ, 607, 765

Fan, X., Strauss, M. A., Richards, G. T., et al. 2006, AJ, 131, 1203

Fan, X., Bañados, E., & Simcoe, R. A. 2023, ARA&A, 61, 373

Fiacconi, D., Mayer, L., Madau, P., et al. 2017, MNRAS, 467, 4080

Greene, J. E., Labbe, I., Goulding, A. D., et al. 2024, ApJ, 964, 39

Hahn, O., & Abel, T. 2013, Astrophysics Source Code Library [record ascl:1311.011]

Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023, ApJ, 959, 39

Hopkins, P. F. 2015, MNRAS, 450, 53

Hopkins, P. F., Cox, T. J., Hernquist, L., et al. 2013, MNRAS, 430, 1901

Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27

Jiang, Y.-F., Blaes, O., Stone, J. M., & Davis, S. W. 2019, ApJ, 885, 144

Johnson, J. L., Khochfar, S., Greif, T. H., & Durier, F. 2011, MNRAS, 410, 919

Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511

Labbe, I., Greene, J. E., Bezanson, R., et al. 2023, ApJ, submitted [arXiv:2306.07320]

Latif, M. A., Bovino, S., Grassi, T., Schleicher, D. R. G., & Spaans, M. 2015, MNRAS, 446, 3163

Latif, M. A., Whalen, D. J., Khochfar, S., Herrington, N. P., & Woods, T. E. 2022, Nature, 607, 48

Lin, D., Guillochon, J., Komossa, S., et al. 2017, Nat. Astron., 1, 0033

Lupi, A., & Bovino, S. 2020, MNRAS, 492, 2818

Lupi, A., Haardt, F., Dotti, M., et al. 2016, MNRAS, 456, 2993

Lupi, A., Volonteri, M., Decarli, R., et al. 2019, MNRAS, 488, 4004

Lupi, A., Pallottini, A., Ferrara, A., et al. 2020, MNRAS, 496, 5160

Lupi, A., Haiman, Z., & Volonteri, M. 2021, MNRAS, 503, 5046

Lupi, A., Volonteri, M., Decarli, R., Bovino, S., & Silk, J. 2022, MNRAS, 510, 5760

Ma, L., Hopkins, P. F., Ma, X., et al. 2021, MNRAS, 508, 1973

Madau, P., Haardt, F., & Dotti, M. 2014, ApJ, 784, L38

Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2023, A&A, submitted

Maiolino, R., Scholtz, J., Witstok, J., et al. 2024, Nature, 627, 59

Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS, 351, 169

Martizzi, D., Faucher-Giguère, C.-A., & Quataert, E. 2015, MNRAS, 450, 504

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

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

Mignoli, M., Gilli, R., Decarli, R., et al. 2020, A&A, 642, L1

Mortlock, D. J., Warren, S. J., Venemans, B. P., et al. 2011, Nature, 474, 616

Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69

Neeleman, M., Novak, M., Venemans, B. P., et al. 2021, ApJ, 911, 141

Neumayer, N., Seth, A., & Böker, T. 2020, A&ARv, 28, 4

Ni, Y., Di Matteo, T., Bird, S., et al. 2022, MNRAS, 513, 670

Padoan, P., Haugbølle, T., & Nordlund, Å. 2012, ApJ, 759, L27

Pezzulli, E., Valiante, R., & Schneider, R. 2016, MNRAS, 458, 3047

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

Planck Collaboration XIII. 2016, A&A, 594, A13

Regan, J. A., Downes, T. P., Volonteri, M., et al. 2019, MNRAS, 486, 3892

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

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

Ricarte, A., Narayan, R., & Curd, B. 2023, ApJ, 954, L22

Sala, L., Cenci, E., Capelo, P. R., Lupi, A., & Dotti, M. 2021, MNRAS, 500, 4788

Sassano, F., Capelo, P. R., Mayer, L., Schneider, R., & Valiante, R. 2023, MNRAS, 519, 1837

Sdowski, A., Lasota, J. P., Abramowicz, M. A., & Narayan, R. 2016, MNRAS, 456, 3915

Schauer, A. T. P., Regan, J., Glover, S. C. O., & Klessen, R. S. 2017, MNRAS, 471, 4878

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

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

Sądowski, A., & Narayan, R. 2016, MNRAS, 456, 3929

Sądowski, A., Lasota, J.-P., Abramowicz, M. A., & Narayan, R. 2016, MNRAS, 456, 3915

Soltan, A. 1982, MNRAS, 200, 115

Springel, V. 2005, MNRAS, 364, 1105

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

Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS, 391, 1685

Stone, M. A., Lyu, J., Rieke, G. H., & Alberts, S. 2023, ApJ, 953, 180

Stone, M. A., Lyu, J., Rieke, G. H., Alberts, S., & Hainline, K. N. 2024, ApJ, 964, 90

Tanaka, T., & Haiman, Z. 2009, ApJ, 696, 1798

Tchekhovskoy, A., & Giannios, D. 2015, MNRAS, 447, 327

Tenneti, A., Di Matteo, T., Croft, R., Garcia, T., & Feng, Y. 2018, MNRAS, 474, 597

Tremmel, M., Governato, F., Volonteri, M., & Quinn, T. R. 2015, MNRAS, 451, 1868

Tremmel, M., Karcher, M., Governato, F., et al. 2017, MNRAS, 470, 1121

Uchiyama, H., Toshikawa, J., Kashikawa, N., et al. 2018, PASJ, 70, S32

Volonteri, M., Silk, J., & Dubus, G. 2015, ApJ, 804, 148

Volonteri, M., Habouzit, M., & Colpi, M. 2021, Nat. Rev. Phys., 3, 732

Xie, F.-G., & Yuan, F. 2012, MNRAS, 427, 1580

Yang, J., Wang, F., Fan, X., et al. 2023, ApJ, 951, L5

Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529

Yue, M., Eilers, A. C., Simcoe, R. A., et al. 2024, ApJ, 966, 176

Zhu, Q., Li, Y., Li, Y., et al. 2022, MNRAS, 514, 5583