The ASTRID simulation: the evolution of Supermassive Black Holes

The formation mechanism of SMBH seeds is an active research topic  (see,e.g. Smith & Bromm, 2019; Volonteri et al., 2021, for recent review). Different seeding channels have been proposed, including remnants of the first generation (Population III) stars (Madau & Rees, 2001) creating $10\sim 10^{3}\,M_{\odot}$ BH seeds, direct collapse of pristine (metal-free) gas (Bromm & Loeb, 2003) and runaway collapse of nuclear star clusters (Davies et al., 2011) with a higher initial seed mass ($\sim 10^{4}-10^{6}\,M_{\odot}$). However, these processes cannot be modelled self-consistently at the resolution of current cosmological simulations.

In Astrid , we continue to follow the practice of seeding a BH after the formation of a sufficiently massive halo. We periodically run a Friends-of-Friends (FOF) group finder during the simulation, with a linking length of $0.2$ times the mean particle separation. Haloes with a total mass and stellar mass exceeding the seeding criteria { $M_{\rm halo,FOF}>M_{\rm halo,thr}$; $M_{\rm*,FOF}>M_{\rm*,thr}$} are eligible for seeding. We apply a mass threshold value $M_{\rm halo,thr}=5\times 10^{9}{h^{-1}M_{\odot}}$ and $M_{\rm*,thr}=2\times 10^{6}{h^{-1}M_{\odot}}$. The choice of $M_{\rm*,thr}$ is a conservative setting to make sure that BHs are seeded in halos with enough cold dense gas to form stars, and that there are at least some collisionless star particles in the BH neighbourhood to act as sources of dynamical friction. The value of $M_{\rm*,thr}$ is estimated from the global $M_{\rm halo}$-$M_{\rm*}$ relation. Most of the FOF halos with $M_{\rm halo}>M_{\rm halo,thr}$ satisfy the $M_{\rm*,thr}$ criteria.

Considering the complex astrophysical process likely to be involved in BH seed formation, we allow haloes with the same mass to have different SMBH seeds. Therefore, in Astrid , instead of applying a uniform seed mass for all the BHs, we probe a range of the BH seed mass $\,M_{\rm sd}$ drawn probabilistically from a power-law distribution

where $\mathcal{N}$ is the normalization factor. We set $M_{\rm sd,min}=3\times 10^{4}{h^{-1}M_{\odot}}$, $M_{\rm sd,max}=3\times 10^{5}{h^{-1}M_{\odot}}$, and a power-law index $n=-1$. In Section 4.4, we investigate the effect of $M_{\rm sd}$ on BH growth.

For each halo that satisfies the seeding criteria but does not already contain at least one SMBH particle, we convert the densest gas particle into a BH particle. The newborn BH particle inherits the position and velocity of its parent gas particle. The intrinsic mass of BH ($M_{\rm BH}$ hereafter) seeds are initially at $M_{\rm sd}$ (as a subgrid mass). However, we preserve in a separate mass variable the parent particle mass, representing a subgrid gas reservoir around the SMBH. Neighbouring gas particles are swallowed once (Eddington-limited) accretion has allowed $M_{\rm BH}$ to grow beyond the initial parent particle mass (i.e., once it depletes the gas reservoir of the parent gas particle).

A common approach applied in cosmological simulations is to anchor the BH at the halo potential minimum by forcefully repositioning the BH particle to the local (within the SPH kernel) potential minimum at every active timestep. While this method is effective at stabilizing the BH in the halo centre, it means that the BH trajectory is discontinuous and the BH velocity is ill-defined. In addition, should two galaxies pass close enough to each other, a BH merger can result immediately even though they are not gravitationally bound.

We instead use a sub-grid dynamical friction model (Chen et al., 2021a) to account for the unresolved frictional force exerted by the surrounding collisionless matter (stars and stellar-mass black holes) on SMBH dynamics. In this model, SMBHs gradually sink to the galaxy centre, as their relative velocity is dissipated by dynamical friction. Moreover, we can make a more physical prediction for BH mergers by calculating whether the (now well-defined) relative velocities of two SMBHs allow them to be gravitationally bound. The validation of our dynamical friction model in cosmological simulations was performed in Chen et al. (2021a). In this section, we briefly review the BH dynamics.

We model BH growth and AGN feedback as in the BlueTides simulation. BHs are allowed to grow by accreting mass from nearby gas particles. The gas accretion rate onto the BH is estimated via the Bondi-Hoyle-Lyttleton-like prescription applied to the smoothed properties of the $112$ gas particles within the SPH kernel of the BH:

$c_{s}$ and $\rho$ are the local sound speed and density of gas, $v_{\rm rel}$ is the relative velocity of the BH with respect to the nearby gas and $\alpha=100$ is a dimensionless fudge parameter to account for the underestimation of the accretion rate due to the unresolved cold and hot phase of the subgrid interstellar medium nearby. Note that hydrodynamically decoupled wind particles are not included in the density calculation of Eq. 8.

We allow for short periods of super-Eddington accretion in the simulation, but limit the accretion rate to $2$ times the Eddington accretion rate

$m_{p}$ is the proton mass, $\sigma_{T}$ the Thompson cross section, c the speed of light, and $\eta=0.1$ the radiative efficiency of the accretion flow onto the BH.

In summary, the total BH accretion rate in the Astrid simulation is:

Numerically, the physical BH mass, $\,M_{\rm BH}$, grows continuously with time, while we separately track the BH dynamic mass by stochastic accretion of nearby gas particles with a probability that statistically satisfies mass conservation. Given that $M_{\rm dyn}$ is temporarily boosted (and initially fixed) for the new-born BH particles for gravity and dynamical friction calculation, we use a separate mass tracer $M_{\rm trace}$ to perform gas accretion when $\,M_{\rm BH}<M_{\rm dyn,seed}$. This term is initialized as the parent particle mass that spawns the BH, and traces $\,M_{\rm BH}$ with stochastic gas accretion until $\,M_{\rm BH}>=M_{\rm dyn,seed}$. After that, $M_{\rm dyn}$ grows following $\,M_{\rm BH}$ by swallowing gas.

To ensure stable accretion, the timestep of the BH particle is set not by the acceleration dynamics, but by being $2$ times longer than the smallest timestep found amongst the gas particles within its density kernel.

AGN thermal feedback is implemented as follows. The SMBH is assumed to radiate with a bolometric luminosity $L_{\rm Bol}$ proportional to the accretion rate $\dot{M}_{\rm BH}$:

with $\eta=0.1$ being the mass-to-light conversion efficiency in an accretion disk according to Shakura & Sunyaev (1973). 5% of the radiated energy is thermally coupled to the surrounding gas that resides within twice the radius of the SPH smoothing kernel of the BH particle.

The AGN feedback energy is isotropically imparted to the nearby gas particles, distributing the energy among them according to the SPH kernel. Gas particles that are heated by AGN feedback dissipate the energy in different channels depending on their thermal properties. Non star-forming gas cools as normal, while star-forming gas heated by the SMBH relax to the temperature corresponding to the effective equation of state for star-forming (c.f. Springel & Hernquist, 2003), on a time scale determined by the cooling time. Note that star-forming gas heated by other means cools to the effective equation of state on the longer relaxation timescale.

For every time step in the simulation, we log the full physical state of all the active BH particles, including the thermal states (density, entropy), dynamics (position, velocity, acceleration from gravity, dynamical friction, gas drag), growth state (mass, accretion rate, accretion of gas and BHs), neighbour information (number of surrounding particles, local velocity dispersion) etc, to a separate catalogue. This high time resolution BH catalogue allows us to follow the detailed growth history of all BHs in the simulation, build the merger trees of BHs, trace the trajectory of BH merger events down to $\,{\rm kpc}$ scale and allow post-processing of BH mergers at sub-resolution scale based on the orbits of BHs.

We start our analysis by showing the evolution of BH mass functions (BHMFs). In the left panel of Figure 2, the coloured solid lines show the BH mass function evolving from $z=7$ to $z=3$, illustrating how the BH mass gradually builds up with cosmic time.

In Astrid , the first $10^{8}\,M_{\odot}$, $10^{9}\,M_{\odot}$ and $10^{10}\,M_{\odot}$ BHs are in place at $z=6.6$, $z=5.3$ and $z=4.5$ respectively. With the 250 $h^{-1}{\rm Mpc}$ volume size of Astrid , we marginally being able to simulate the population of the first $\,M_{\rm BH}>10^{9}\,M_{\odot}$ QSOs at $z>6$, which are thought to power high-$z$ quasars. At our final redshift of $z=3$, the simulation contains a hundred $10^{9}\,M_{\odot}$ BHs and a handful ($\sim 3$) of ultramassive BHs with $\,M_{\rm BH}>10^{10}\,M_{\odot}$.

As the first sanity check, we show in the left panel of Figure 2 some observational inferences of the BHMF derived from high-$z$ AGN samples. The purple dashed line gives the estimate of the $z\sim 6$ BHMF derived in Willott et al. (2010) based on optical QSO samples at $z=6$. The brown dotted, dashed-dotted and dashed lines are the BHMF predictions at $z=5$, $z=4$ and $z=3$ from Ueda et al. (2014), based on X-ray observations of AGN samples (assuming mass-dependent radiative efficiency and averaged Eddington ratio $\bar{\lambda}_{\rm Edd}\sim 0.25$). At masses $\,M_{\rm BH}>10^{7}\,M_{\odot}$, Astrid produces a $\,M_{\rm BH}$ distribution in good agreement with these observational predictions from $z=6$ to $z=3$. A caveat is that there are large uncertainties (and inconsistencies with the simulation model) lying in the assumptions applied by these observational inferences of BHMFs. The match to within an order of magnitude across a wide range of $\,M_{\rm BH}$ and redshift still validates that the BH model in Astrid seems to be building up the $\,M_{\rm BH}$ population in a reasonable way.

BH populations at the small mass end are sensitive to the seeding prescription and to the level of growth achieved before reaching AGN feedback modulated self-regulation. In Astrid , the BH population with $\,M_{\rm BH}<10^{6}\,M_{\odot}$ is dominated by the seed BHs. The dashed vertical lines mark the range of the BH seed with $M_{\rm sd,min}=3\times 10^{4}{h^{-1}M_{\odot}}$ ($4.4\times 10^{4}\,M_{\odot}$), $M_{\rm sd,max}=3\times 10^{5}{h^{-1}M_{\odot}}$ ($4.4\times 10^{5}\,M_{\odot}$). Instead of spiking at a uniform $\,M_{\rm sd}$, the power-law seeding prescription applied in Astrid allows the BH population to extend down to a lower $\,M_{\rm sd}$ with a flatter slope, and allows a statistical diagnosis of the impact of $\,M_{\rm sd}$ on BH growth, which will be further discussed in Section 4.5. We will show there that $\,M_{\rm sd}$ has a limited impact on the higher mass BH population, for which growth is mostly dominated by the local environment, governing gas accretion.

Now we investigate the relationship between $\,M_{\rm BH}$ and luminosity across different redshifts. The AGN luminosity is calculated from the BH accretion rate using bolometric luminosity corrections. We calculate the bolometric luminosity of AGN via Eq. 11 with $L_{\rm Bol}\propto\dot{M}_{\rm BH}$ and the assumption of a uniform mass-to-light conversion efficiency of $\eta=0.1$. To compare with AGN observations in different bands such as the X-ray and UV, we apply the corresponding bolometric corrections.

In particular, the X-ray band luminosity is estimated by converting AGN $L_{\rm Bol}$ to the luminosity in the hard X-ray band [2-10] keV following the bolometric correction $L_{X}=L_{\rm Bol}/k$ from Hopkins et al. (2007) with

Here $L_{\rm Bol}$ is given by Eq. 11, and $L_{\odot}$ refers to the bolometric solar luminosity $L_{\odot}=3.9\times 10^{33}$ erg s^-1.

To calculate the UV band AGN luminosity, we convert $L_{\rm Bol}$ to rest frame UV band absolute magnitude $M_{\rm UV}$ following the bolometric corrections from Fontanot et al. (2012)

where $f_{B}=10.2$, $\mu_{B}=6.7\times 10^{14}$Hz and $\Delta_{\mathrm{B,UV}}$=-0.48.

The right panel of Figure 2 shows the distribution of $\,M_{\rm BH}$ versus the AGN luminosity. The left $y$-axis is labelled in units of bolometric luminosity, while the right $y$-axis marks the corresponding $L_{\rm X}$ calculated via Eq. 12. The background shows the 2D histogram of the $L_{\rm BH}-\,M_{\rm BH}$ relation for the BH population at $z=3$, with the purple solid line giving the mean value of $L_{\rm BH}$ for each $\,M_{\rm BH}$ bin. We also give the average $L_{\rm BH}-\,M_{\rm BH}$ for higher redshifts from $z=7$ to $z=3$ (using coloured solid lines with the same colour scheme as the left panel) to show time evolution.

At the low mass end, $\,M_{\rm BH}\lesssim 10^{7}\,M_{\odot}$, $L_{\rm BH}$ exhibits a steep positive correlation with $\,M_{\rm BH}$ as the BH growth is dominated by the initial Bondi-Hoyle-like accretion. More massive BHs experience a higher accretion rate, which also gives on average a higher Eddington accretion ratio $\lambda_{\rm Edd}\equiv\dot{M}/\dot{M}_{\rm Edd}$. As shown by the $L_{\rm BH}-\,M_{\rm BH}$ distribution at $z=3$ in the background, some of the BHs with $\,M_{\rm BH}\sim[10^{6},10^{7}]\,M_{\odot}$ can temporarily reach Eddington-limited accretion (marked with a grey dashed-dotted line). At increased mass, $\,M_{\rm BH}\gtrsim 10^{7}\,M_{\odot}$, the $L_{\rm BH}-\,M_{\rm BH}$ curve flattens, making the overall Eddington accretion ratio lower for more massive BHs. This signifies that BH growth is self-regulated by AGN feedback, with massive BHs dumping thermal energy into the surrounding gas and the resulting hot and sparse gas slowing down the BH accretion. The overall $\lambda_{\rm Edd}$ tends to peak at $\,M_{\rm BH}\sim 10^{7}\,M_{\odot}$, a typical feature found in Bondi-Hoyle-like BH accretion models (see, e.g. DeGraf et al., 2012).

As shown by the coloured solid lines, the mean $L_{\rm BH}-\,M_{\rm BH}$ relation exhibits a clear time evolution across different redshifts, with the mean Eddington accretion ratio for BHs in the same $\,M_{\rm BH}$ regime decreasing when going to lower redshifts. For $\,M_{\rm BH}\sim 10^{8}\,M_{\odot}$, the average $\lambda_{\rm Edd}\sim 0.07$ at $z=3$, $\sim 0.12$ at $z=4$, and $\sim 0.2$ at $z=5\sim 6$. The decrease in $\lambda_{\rm Edd}$ with redshift is a consequence of persistent AGN feedback and lower cold gas fractions in BH environments at lower redshift (see discussions in, e.g., Anglés-Alcázar et al., 2015), and is mirrored in a similar decrease in star formation (see Bird et al., 2021). As we will show in Section 4.3.2, AGN at higher redshift is typically embedded in denser gas environments. Observational evidence also suggests an increasing $\lambda_{\rm Edd}$ at higher redshift (e.g. Lusso et al., 2012; Aird et al., 2015), supporting the trend found in the simulation.

In this section, we discuss the predictions for the AGN luminosity functions (in X-ray and UV bands) for the BH population in Astrid and compare them to current observational constraints at $z\sim 3-7$. The AGN luminosity functions (LFs: 1-point functions) provide basic statistics to gauge the accretion history of BHs as tracked in our simulation and hence are a fundamental way to validate our BH model. In addition, following Ni et al. (2020), we calculate the hydrogen column density, $N_{\rm H}$ due to the gas in the host galaxies and assess its role in AGN obscuration. Although we cannot resolve scales relevant for the obscuring AGN torus, galactic obscuration serves as a lower limit on the total amount of obscuration expected. We compare our derived galactic $N_{\rm H}$ as a function of X-ray luminosity to the current obscured fractions of AGN inferred from X-ray observations at $z>3$.

This section discusses the global evolution history of the BH mass density (BHMD) and BH accretion rate density (BHAD) from $z=10$ to $z=3$. The left panel of Figure 6 shows the BHMD as a function of redshift. The red to orange coloured lines denote the BHMD contributed by BHs in different ranges of $M_{\rm BH}$ (as indicated in the legend). The solid triangles indicate the upper limits derived from deep X-ray observations (Treister et al., 2013) or from the integrated X-ray background (Salvaterra et al., 2012). Note that those X-ray inferred upper limits do not include Compton-thick AGNs and might thus miss the contributions from a fraction of BHs. The black dotted line marked by empty circles was derived by Ueda et al. (2014) from observational constraints on the AGN XLF. The grey shaded area denotes the predicted range of BHMD for BHs with $\,M_{\rm BH}>10^{5}\,M_{\odot}$ reported by Volonteri & Reines (2016), which is derived from the galaxy mass function (Grazian et al., 2015) after making assumptions about the BH-stellar mass relationships. We can see that the Astrid BHMD is in overall good agreement with the observational constraints. BHs within the mass range [$10^{5}$, $10^{6}$] $\,M_{\odot}$ make the dominant contribution to BHMD.

The right panel of Figure 6 shows the BHAD with different coloured lines showing contributions from different $L_{\rm X}$ bins. At $z>5$ the BHAD is dominated by AGN with relatively low luminosity in the range of $\log L_{\rm X}$ = [41, 42.5] erg s^-1, while at lower redshift ($z<5$) more luminous AGNs with $\log L_{\rm X}>42.5$ start to dominate the BH accretion rate. The black data points show the observational results from Vito et al. (2018), derived from samples of X-ray-detected AGN. As a comparison, we also show the simulation result from Illustris (Sijacki et al., 2015) (including all BHs) as the black dotted line. The results of Vito et al. (2018) effectively probe AGN down to the luminosity limit $\log L_{\rm X}>42.5$ for $z=3\sim 6$, but beyond this ($\log L_{\rm X}<42.5$) the sample is strongly affected by incompleteness. We can see that the dashed and dotted lines (corresponding to the $\log L_{\rm X}>42.5$ AGN population) show overall good agreement with Vito et al. (2018) at $z=6$ and at $z\sim 4$, but overshoot by about 0.5 dex at $z=5$.

The lower left panel shows the ratio of the stellar mass density (SMD) to the BHMD at each redshift and the lower right panel shows star formation rate density (SFRD) divided by the BHAD. Here SMD and SFRD are all calculated based on the galaxies with $M_{*}>10^{7}\,M_{\odot}$ in the simulation. Since the star formation process builds up galaxy mass faster than BH formation and accretion (as also indicated by the SFRD/BHAD ratio), the ratio between SMD and BHMD slowly increases with time, from $\sim 100$ at $z=10$ to $\sim 700$ at $z=3$. However, we note that the ratio of SFRD and BHAD does not show very strong time evolution, with SFRD/BHAD lying between $1000-3000$ for all redshifts. Some observations (e.g. Mullaney et al., 2012) and simulations (e.g. Ricarte et al., 2019; Thomas et al., 2019) also exhibit no strong redshift evolution of the relation between star formation rate and BH accretion rate, and treat this as evidence for or a consequence of BH-galaxy coevolution.

Since the Bondi-Hoyle prescription for BH accretion depends on the BH mass with $\dot{M}_{\rm BH}\propto\,M_{\rm BH}^{2}$ (c.f. Eq. 8), one would expect the early growth of BHs to be sensitive to the initial $\,M_{\rm sd}$. This is because in the early phase of BH growth the AGN feedback does not have a large impact on the local environment and the BH accretion is dominated by the $\,M_{\rm BH}^{2}$ term, where a larger $\,M_{\rm sd}$ would cause $\,M_{\rm BH}$ to build up more rapidly.

In Astrid , BHs are seeded with $\,M_{\rm sd}$ stochastically drawn from a power-law distribution within the mass range from $M_{\rm sd,min}=4.4\times 10^{4}\,M_{\odot}$ ($3\times 10^{4}{h^{-1}M_{\odot}}$) to $M_{\rm sd,max}=4.4\times 10^{5}\,M_{\odot}$, which is independent of the host environment at the seeding time. This allows us to statistically investigate how the growth of a BH is dependent on its seeding mass.

To this end, we select a batch of BHs that are seeded in a given time interval, trace them to a later time, and show how their final $M_{\rm BH}$ are dependent on their initial $\,M_{\rm sd}$. In order to allow the BHs enough time evolution, we select the BHs seeded at high redshift between $10<z<15$ (about $8.6\times 10^{4}$ BHs in total), and trace them to lower redshifts from $z=6$ to $z=3$. To take into account the fact that BHs also grow by merging with other BHs (with a different seed mass), we split the analysis into two parts. First, we exclude the BHs that have undergone merger events to cleanly diagnose how $M_{\rm BH}$ gained by the gas accretion is affected by original seed mass. Second, we include the BHs that have merged with other smaller mass BHs to address the possible bias introduced by selecting the no-merger BHs.

The upper panels of Figure 7 show the result of the $\,M_{\rm BH}-\,M_{\rm sd}$ relation for BHs that are seeded at $10<z<15$ and have not undergone a merger before the redshift inspected. The blue points in the left panel give the resultant $\,M_{\rm BH}$ at $z=3$ versus their original $\,M_{\rm sd}$, with the black solid and dashed-dotted line giving the running mean and median $\,M_{\rm BH}$ for each $\,M_{\rm sd}$ bin. The middle panel shows the statistics of the mean and median $\,M_{\rm BH}$/$\,M_{\rm sd}$ in each $\,M_{\rm sd}$ bin from $z=6$ to $z=3$ to show time evolution. The right panel shows the fraction of BHs in each $\,M_{\rm sd}$ bin that have grown beyond their initial $\,M_{\rm sd}$ at each redshift, with the solid and dashed lines representing $\,M_{\rm BH}>2\times\,M_{\rm sd}$ and $\,M_{\rm BH}>10\times\,M_{\rm sd}$ respectively.

We can see that the initial $\,M_{\rm sd}$ has a clear imprint on the final $\,M_{\rm BH}$, even after time evolution from $z=10$ to $z=3$. The median $\,M_{\rm BH}-\,M_{\rm sd}$ relation at $z=3$ displays a strong positive correlation especially at the lower $\,M_{\rm sd}$ end, indicating that for the BH population from the smallest $\,M_{\rm sd}$ bin, their BH accretion in general lags behind at the early stages and has not yet caught up even by $z=3$. The right panel tells that BHs with the smallest $\,M_{\rm sd}$ have a significantly lower chance to grow beyond the seed mass compared to BHs with more massive $\,M_{\rm sd}$: less than $50\%$ of BHs from the smallest $\,M_{\rm sd}$ bin have grown to $\,M_{\rm BH}>2\times\,M_{\rm sd}$ at $z=3$.

However, it is important to note that these results are not showing that $\,M_{\rm BH}$ is predominantly determined by the $\,M_{\rm sd}$. We can see from the left panel that the BHs from the lowest $\,M_{\rm sd}$ bin can still reach the highest $\,M_{\rm BH}$ range ($\,M_{\rm BH}>10^{9}\,M_{\odot}$) through pure gas accretion, indicating that the local environment specific to each BH plays a significant role in BH growth regardless of $\,M_{\rm sd}$.

Moreover, the $\,M_{\rm BH}$ difference caused by different initial $\,M_{\rm sd}$ values becomes narrower as time evolves, as the $\,M_{\rm BH}$ becomes dominated by the accretion process determined by the local environment. The evolution of the $\,M_{\rm BH}$-$\,M_{\rm sd}$ relation shown in the right two panels becomes flatter, especially for $\,M_{\rm sd}>10^{5}\,M_{\odot}$ when moving to lower redshift. At $z=3$, the median $\,M_{\rm BH}/\,M_{\rm sd}$ ratio, as well as the fraction of $\,M_{\rm BH}$ growth is almost the same for BHs from $\,M_{\rm sd}\sim 10^{5}\,M_{\odot}$ and from the most massive $\,M_{\rm sd}$ bin. This indicates that for BHs seeded with $\,M_{\rm sd}>10^{5}\,M_{\odot}$ at $z>10$, the impact of $\,M_{\rm sd}$ on the $\,M_{\rm BH}$ at $z=3$ has been virtually erased by accretion and its associated exponential growth.

Selecting BHs that have not undergone any mergers might artificially eliminate the BH population residing in massive galaxies that are also likely to be gas-rich. Therefore in the lower panel of Figure 7 we repeat the same analysis including the growth of BHs due to mergers. In this case, the $\,M_{\rm sd}$ (in the lower panel of Figure 7) is that of the more massive BH progenitor in the merger history. The resultant $\,M_{\rm BH}$-$\,M_{\rm sd}$ relation displays the same qualitative features as in the upper panels, but with a higher $\,M_{\rm BH}$ than obtained in the upper panel.

Therefore, we can conclude that the initial $\,M_{\rm sd}$ has an impact on $\,M_{\rm BH}$ growth by gas accretion (and additionally by BH mergers, when they occur), whereby BHs from the smallest $\,M_{\rm sd}$ bin have a lower chance to grow beyond the seed mass even with the long period of time evolution from $z>10$ to $z=3$. However, the effect of the $\,M_{\rm sd}$ would be erased with time evolution, and the $\,M_{\rm BH}$ is not dominated by the initial $\,M_{\rm sd}$. The most massive BHs can grow from the smallest initial $\,M_{\rm sd}$, as the local environment for each specific BH plays a crucial role in their evolutionary history.

In Section 4.5, we have seen that the BHs seeded from the smallest $\,M_{\rm sd}$ bin can still grow to $>10^{9}\,M_{\odot}$ by gas accretion only. This indicates that gas accretion plays a predominant role in the mass assembly history of at least some BHs. It is interesting to analyse and compare the contribution of BH mergers and gas accretion to the overall $\,M_{\rm BH}$ distribution on a statistical basis and investigate its time evolution.

To this end, in Figure 8 we present the mass fraction of $\,M_{\rm BH}$ contributed by $\,M_{\rm sd}$ and by mergers as a function of $\,M_{\rm BH}$ at two different redshifts, $z=6$ and $z=3$. For each BH that has undergone mergers, we trace the more massive progenitors, use the $\,M_{\rm sd}$ from the massive progenitor and sum the mass of the merged BHs (calling this the mass contribution from mergers). The orange line in Figure 8 shows the mass fraction contributed by $\,M_{\rm sd}$. It lies between $M_{\rm sd,min}/\,M_{\rm BH}$ and $M_{\rm sd,max}/\,M_{\rm BH}$, as marked by the grey dashed lines. The purple and blue lines give the averaged mass fraction of $\,M_{\rm BH}$ obtained by BH accretion and by swallowing other BHs.

BHs in the mass range $\,M_{\rm BH}\lesssim 10^{6}\,M_{\odot}$ are dominated by the seed mass, with the $\,M_{\rm sd}$ fraction lying toward the upper limit of $M_{\rm sd,max}/\,M_{\rm BH}$. This is because the low mass BHs are more likely to trace the initial seed mass and BHs in the $10^{6}\,M_{\odot}$ regimes have a larger seed contribution compared to those with massive seeds. The $\,M_{\rm sd}$ fraction becomes less biased when moving to the more massive end because massive BHs gain most of their mass from the gas accretion which is determined by the local environment. Their mass, therefore, is less dependent on the initial $\,M_{\rm sd}$.

The averaged mass fraction of $\,M_{\rm BH}$ from BH mergers has an almost flat distribution at $z=3$, with no strong dependence on the BH mass. The strict BH merger criteria applied in Astrid allows BHs to have a wide range of merger histories and leads to a large scatter in the BH mass fraction contributed by the merger. The scatter of the BH merger component (blue shade) increases with $\,M_{\rm BH}$ as more massive BHs can involve BH mergers with larger mass contrast.

The BH merger mass fraction in the $\,M_{\rm BH}>10^{7}\,M_{\odot}$ regime at $z=3$ is slightly larger than that at $z=6$, as more merger events happen at lower redshift. However, the averaged contribution to $\,M_{\rm BH}$ from BH mergers is less than 10% for all mass bins and at both redshifts. For the population of the most massive BHs with $\,M_{\rm BH}>10^{7}\,M_{\odot}$, on average, more than 90% of the BH mass is obtained through the gas accretion, with only a few per cent from the BH mergers and negligible contribution from the seed mass. Therefore, we conclude that in Astrid , BH mergers are a subdominant factor in BH growth until at least $z=3$.

Simulations where BHs are centered using repositioning rather than dynamic friction merge BHs more aggressively and thus make different predictions. For example, Illustris-TNG (Weinberger et al., 2018) predicts that the SMBHs with $\,M_{\rm BH}>10^{8.5}\,M_{\odot}$ grow most of their mass via mergers with lower mass BHs (though this claim is for BH population at $z=0$). A more thorough comparison of the BH mergers with other simulations will be presented in our follow-up work.

Scaling relations between central BH mass and host galaxy properties are of fundamental importance to studying BH and galaxy evolution through cosmic time. In Figure 9, we show the scaling relations between $\,M_{\rm BH}$ and $\,M_{*}$ from $z=6$ to $z=3$ in Astrid . Here we use the total stellar mass of the galaxy in the subhalo and do not apply bulge-disk decomposition to identify the galaxy bulge. We only consider the central galaxies in each FOF halo (the main subhalo in each FOF halo with the greatest stellar mass) and do not include satellite galaxies in this analysis. For each subhalo that contains multiple BHs, we select the most massive one. We have also repeated the analysis by selecting the most luminous BH for each subhalo and this does not lead to any significant difference. The grey dashed and the brown dotted line shows the best-fitting relations from Kormendy & Ho (2013); Reines & Volonteri (2015) based on the AGN population in the nearby universe. We also collect the data points from observations of high-$z$ QSOs. The squared data points are properties of three $z>7$ QSOs summarized in Izumi et al. (2019) and the circles with error bars show Subaru HSC observations of $z>6$ QSOs from a recent compilation by Izumi et al. (2021). The observational stellar masses presented here are the inferred dynamical masses of the host galaxies. The star and triangle-shaped data points are for $z>4$ and $2<z<4$ quasars collected from Kormendy & Ho (2013).

Although the box size of Astrid is relatively large, it is not large enough to produce any $10^{8}\,M_{\odot}$ BH at $z>7$ or $\sim 10^{9}\,M_{\odot}$ BHs at $z>6$, and so we cannot directly compare with $z>6$ QSO observations. The most massive BHs ($\,M_{\rm BH}>10^{8}\,M_{\odot}$) at $z=6$ reside in $\,M_{*}\sim 10^{11}\,M_{\odot}$ galaxies, and minimal extrapolation of these Astrid results leads to qualitatively good agreement with the HSC observations. The BH population with mass $\,M_{\rm BH}<10^{6}\,M_{\odot}$ is dominated by BHs which are at the inefficient initial growth phase close to the seed mass. There is a large variation with regard to the host galaxy mass, ranging from $\,M_{*}=10^{8}\sim 10^{10}\,M_{\odot}$. At $\,M_{\rm BH}>10^{6}\,M_{\odot}$, where BHs go beyond the seed mass, Astrid produces a clear positive correlation between the stellar mass of a galaxy and the mass of its central BH, in concordance with the expectation that galaxies and BHs grow commensurately in a globally averaged sense.

Compared with fits to AGN populations in the local universe from Kormendy & Ho (2013); Reines & Volonteri (2015), Astrid has $\,M_{\rm BH}$ averaging about 0.5 dex lower for each $\,M_{*}$ bin, with the discrepancy getting narrower at the more massive end with $\,M_{*}>10^{11}\,M_{\odot}$. The lower averaged $\,M_{\rm BH}$ predicted by the simulation can be partially attributed to the population of small BHs that are difficult to obtain using observations. Moreover, we note that this discrepancy is of the same order as likely systematic errors in the observations and that these systematics will be worse for low mass objects.

As a comparison with other simulations, we also plot using a pink dashed-dotted line the $\,M_{\rm BH}-\,M_{*}$ relation measured from the SIMBA simulation Davé et al. (2019). We can see that even with different implementations of the BH accretion and feedback, Astrid and SIMBA produce $\,M_{\rm BH}-\,M_{*}$ relations in good agreement with each other.

Since active BH growth is induced by the same cold dense gas that fuels star formation, one expects that the star formation rate and BH accretion rate should be positively correlated. This is supported by a range of AGN observations (e.g. Mullaney et al., 2012; Chen et al., 2013; Delvecchio et al., 2015; Lanzuisi et al., 2017) and is also established in many simulations (e.g. Thomas et al., 2019; Ricarte et al., 2019).

In Figure 10, we show the instantaneous $\dot{M}_{\rm BH}$ - SFR relationship from $z=6$ to $z=3$. The solid lines are for the median $\dot{M}_{\rm BH}$ for each SFR bin, with the error bars showing the range corresponding to the 16-84th percentile. The grey squares are observational samples collected from the work of Delvecchio et al. (2015) who combined IR band selected star-forming galaxies and X-ray detected AGNs over the redshift range covering $1.5<z<2.5$.

The SFR broadly trace the $\dot{M}_{\rm BH}$ for galaxies with $\mathrm{SFR}\gtrsim 1\,M_{\odot}/{\rm yr}$, and show reasonable agreement at the order of magnitude level compared to observational data for active star-forming galaxies (Delvecchio et al., 2015). The curves flatten at the low SFR end, as smaller galaxies with lower SFR host BHs near the seed mass, where $\dot{M}_{\rm BH}$ is dominated by $\,M_{\rm BH}$ and is not on average coupled to SFR.

The $\dot{M}_{\rm BH}$-SFR relation at different redshifts does not show obvious time evolution for active star-forming galaxies with $\mathrm{SFR}\gtrsim 100\,M_{\odot}/{\rm yr}$. For AGN hosts with SFR $<100\,M_{\odot}/{\rm yr}$, the $\dot{M}_{\rm BH}$-SFR relation exhibits substantial scatter, and the median $\dot{M}_{\rm BH}$ for the same SFR decreases when going to higher redshift. This is because, for smaller hosts, a larger fraction of BHs have not built up the mass to reach the self-regulated regime and $\dot{M}_{\rm BH}$ is dominated by the initial Bondi accretion that is sensitive to $\,M_{\rm BH}$, where $\,M_{\rm BH}$ is smaller at higher redshift. We note that the $\,M_{\rm BH}-\,M_{*}$ relation shows a similar trend with redshift for smaller galaxies.

As another probe of the connection or feedback between SMBH growth and galaxy build-up, observations are used to investigate the relationship between AGN gas obscuration and the host galaxy mass, with some finding a positive correlation between $N_{\rm H}$ and $M_{*}$ (e.g. Rodighiero et al., 2015; Buchner et al., 2017; Lanzuisi et al., 2017).

In Figure 11, we show the simulation predictions for the $N_{\rm H}$ - $\,M_{*}$ relation. We investigate the $L_{\rm X}>10^{42.5}$ erg s^-1 AGN populations and their host galaxies at $z=6$ and $z=3$. The black solid line surrounded by a shaded area shows the running median and 16-84th percentiles of the $N_{\rm H}$ distribution including all sightlines for AGN lying in each $\,M_{*}$ bin. For comparison, the red squares in Figure 11 show observational results from Rodighiero et al. (2015) for a sample of $z\sim 2$ AGN hosts in the COSMOS field. The purple shaded region is the linear fit of $N_{\rm H}$ to $\,M_{*}$ for the $2<z<4$ QSO population measured by Lanzuisi et al. (2017), based on a sample of X-ray detected AGN and their far-UV detected host galaxies.

Astrid predicts an overall higher normalization for $N_{\rm H}$ at $z=6$ compared to $z=3$, as the galactic obscuration at higher redshift is larger. Considering the large variations brought about by different AGN sightlines, Astrid predicts an $N_{\rm H}$-$\,M_{*}$ relation broadly consistent with observations. We remind the reader that $N_{\rm H}$ calculated in the simulation is the galactic scale gas obscuration, which does not include resolving the nuclear torus, and is therefore a lower limit to the obscuration that would be probed by AGN observations. The consistency with observations implies that a significant fraction of the obscuration observed in AGN is due to galaxy-scale gas in the host (as also claimed in, e.g. Buchner et al., 2017).

Astrid predicts a weak positive correlation between $N_{\rm H}$ and $\,M_{*}$ at both $z=6$ and $z=3$, implying that as galaxy mass increases there are higher chances to have an additional component to the amount of gas along the line of sight towards the AGN. However, we also note that the variation predicted by the simulation is large. At $z=3$, the $1\sigma$ variation for massive galaxies with $\,M_{*}>10^{11}\,M_{\odot}$ is about 1.5 dex, ranging from $N_{\rm H}=10^{22}\rm{cm}^{-2}$ to $10^{23.5}\rm{cm}^{-2}$. As discussed in our previous work, Ni et al. (2020), simulations predict a large angular variation of $N_{\rm H}$ values along different AGN sightlines, as the AGN feedback modulates the surrounding gas density and launches gas outflows that lead to windows of low obscuration within even very massive galaxy hosts. We therefore do not see a tight power-law relation between $N_{\rm H}$ and $\,M_{*}$ due to this large variation between different sightlines.

There is compelling evidence that almost every massive galaxy in the nearby universe contains an SMBH at its centre. However, this may not be true for fainter AGN in low mass galaxies and the simulation will help us determine what to expect and why. The fraction of galaxies hosting AGN as a function of galaxy mass and other properties also provide useful information about the growth and fueling of the first BHs. In addition, it is interesting to investigate the statistics of BH occupation numbers in galaxies using a simulation that includes physically-based modelling of BH mergers and how they sink to the centre of their hosts. In Astrid , BHs do not merge immediately following the mergers of host galaxies, but instead, their relative kinetic energies are dissipated by the force of dynamical friction and they become gravitationally bound. As a consequence, large galaxies in Astrid will in general host multiple BHs. In this section, we show statistics of the BH occupation in galaxies.

Figure 12 provides a visualization of the galaxy hosts for BHs. The upper panel illustrates how BHs reside in galaxies associated with different subgroups, all from the same large FOF group. We show some examples of large systems at $z=4$ and $z=3$. Different subgroups are identified by coloured circles that correspond to the virial radius $R_{200}$ of that subgroup. The dashed circles identify the central galaxy (associated with the most massive subgroup). The coloured crosses mark the positions of the BHs corresponding to different subgroups, with the cross size scaled by the BH mass. The lower panels of Figure 12 show face-on and edge-on zoom-in views of the galaxy hosts for a selection of massive BHs at $z=4$ and $z=3$. The colour of the galaxies is determined by the age of the stars, with redder colours representing older stellar populations.

The upper panel of Figure 12 illustrates how massive galaxies can host multiple BHs. In Figure 13, we present quantitative information on BH occupation as a function of the stellar mass in the host galaxy. The upper panel shows the fraction of galaxies that contain at least one BH (occupation fraction) in each stellar mass bin, from $z=6$ to $z=3$. The green line includes all BHs, while the red line includes only the BHs that have grown beyond the seed mass regime with $\,M_{\rm BH}>10^{6}\,M_{\odot}$.

We also apply luminosity thresholds to predict BH occupation for bright AGNs and also for fainter AGNs that could be detected by the future X-ray observatory Lynx. The Lynx flux limit of $1\times 10^{-19}$ ergs cm^-2 $s^{-1}$ in $2-10$ keV band (see, e.g. The Lynx Team, 2018) allows it to probe an AGN population with $L_{\rm X}\gtrsim 10^{40-41}$ erg s^-1 at redshifts $z=7\sim 3$. The blue and grey dashed lines represent the occupation number of BHs with $L_{\rm X}>10^{41}$ erg s^-1 and $L_{\rm X}>10^{43}$ erg s^-1 respectively.

The lower panel of Figure 13 displays the BH occupation number as a function of galaxy stellar mass, showing the predicted number of BHs for galaxies with given $\,M_{*}$. The coloured lines give the mean occupation number and the shaded area covers the 16-84th percentiles for all galaxies in each $\,M_{*}$ bin.

Across all redshifts from $z=6$ to $z=3$, galaxies with $\,M_{*}=10^{9}-10^{10}\,M_{\odot}$ begin hosting BHs that have grown beyond the seed mass or are relatively active with luminosities $L_{\rm X}>10^{41}$ erg s^-1. This is qualitatively consistent with the findings of many other simulations such as Illustris, Illustris-TNG, Horizon-AGN and SIMBA, where $\,M_{*}\sim 10^{10}\,M_{\odot}$ galaxies have a probability of ¿ 80 per cent to host an efficient accretor at $z=3\sim 4$ (with threshold $L_{\rm bol}>10^{43}$ erg s^-1, see Habouzit et al. (2022) for a comprehensive review). On the other hand, most galaxies below $10^{9}\,M_{\odot}$ contain only BHs near the seed mass. The deficiency of BH growth in low mass galaxies could be explained by the lack of cold dense gas that fuels BH accretion in the shallower potential well of smaller galaxies. Moreover, multiple simulation studies have found that BH growth in low mass galaxies with $\,M_{*}<10^{9}\,M_{\odot}$ is regulated or stunted by supernovae (SN) feedback so that low mass galaxies generally have a lower probability of hosting a BH (see, e.g., Dubois et al., 2015; Habouzit et al., 2017).

The green line in the lower panel of Figure 13 shows that the mean occupation number of BHs quickly increases with stellar mass. A similar feature is found in other cosmological simulations where BHs are permitted to evolve dynamically without being fixed to halo centres (e.g. Volonteri et al., 2016; Tremmel et al., 2018; Ricarte et al., 2021). In particular, Ricarte et al. (2021) shows that the occupation number of wandering BHs scales with halo or galaxy mass as a power-law relation.

We note that the occupation number of more massive and more luminous BHs becomes flattened for large galaxies. At $z=3$, the most massive $\,M_{*}>10^{11}\,M_{\odot}$ galaxies on average host 2 BHs with $\,M_{\rm BH}>10^{6}\,M_{\odot}$, and about $1\sim 2$ actively accreting BHs with $L_{\rm X}>10^{41}$ erg s^-1. The slower increase in luminous AGN with galaxy mass can be explained by efficient AGN feedback in massive galaxies which self-regulates the BHs and so suppresses the number of rapid accretors. A similar feature is also seen in some other simulations like Illustris-TNG and SIMBA that apply very strong AGN feedback for massive BHs (see, e.g. Habouzit et al., 2021, 2022, for a review). In those cases, they find that massive galaxies with $\,M_{*}\sim 10^{11}\,M_{\odot}$ are even less likely to host AGN with rapid accretion at $z\leq 3$ compared to galaxies with $\,M_{*}\sim 10^{10}\,M_{\odot}$. Another reason of the flat feature of the BH occupation number in massive galaxies is likely caused by BH mergers. As we will show in Section 6, more massive galaxies, in general, have undergone more BH mergers. Also, more massive BHs on average have a higher merger rate compared to seed mass BHs. Therefore, even if massive galaxies contain several unmerged seed mass BHs, the average number of central actively accreting BHs remains $1\sim 2$.

We begin with some illustrations of the BH mergers in Astrid . Figure 14 shows three examples of BH merger events, happening at different redshifts and with different $\,M_{\rm BH}$. The top three panels plot the trajectories of the merging BHs together with the evolution of their host galaxies shown in the background. All the snapshots are centred on the position of the primary black hole (BH1, red point). The blue lines show the trajectories of the secondary black hole (BH2, the blue point) in the reference frame of BH1 over a time period with a redshift interval $dz=0.2$. The smooth trajectory of the merging BHs is brought about by the well-defined BH particle velocity.

The lower panels of Figure 14 show the evolutionary histories of the BH masses and luminosities for the merger examples, with the red and blue lines representing BH1 and BH2 separately. We give the X-ray band luminosity (c.f. Eq. 12) modelled from the BH accretion rate. Note that the merging BHs in these three examples are well above the Lynx detection limit of $L_{\rm X}\gtrsim 10^{40-41}$ erg s^-1.

In both examples 1 and 2, the BH merger takes place not long after the encounter or merger of their host galaxies (as shown in the first panel). For the first few close encounters, the BH binaries have large relative velocities that do not satisfy the merger criteria. With their kinetic energy gradually dissipated by dynamical friction, the binary systems become gravitationally bound and the BHs (numerically) merge at the last encounter. Example 1 shows a relatively eccentric trajectory before the merger, while in Example 2 the BH relative orbit is more circular. The masses of the BHs in the Example 2 binary are very close to each other with $q\equiv M_{2}/M_{1}=0.8$. They grow commensurately during the merger process, as also shown by their respective light curves $L_{\rm X}$ in the lower panel. In both examples 1 and 2, the luminosities of the binary BHs are similar at the time of the merger, as the two BHs are close in $\,M_{\rm BH}$ and are embedded in similar gas environments by the time they reach the merger criterion $|\bf{\Delta r}|$ $<3h^{-1}\,{\rm kpc}$.

The third panel shows an example where the merger happens a long time after the first encounter between the host galaxies. In this case, BH2 lingers in the host galaxy before eventually sinking to the galaxy centre and merging with BH1. In this case the two BHs have a larger mass difference, with $M_{2}\sim 10^{6}\,M_{\odot}$ and $M_{2}\sim 10^{8}\,M_{\odot}$ at $z=3.5$. BH2 accretes strongly, reaching a larger mass of $10^{7}\,M_{\odot}$ at $z=3.1$ during the process of sinking to the galaxy centre. It has slightly lower final luminosity compared to BH1, as $M_{2}$ is still an order of magnitude less than $M_{1}$.

In this section, we examine where the merging BHs reside with respect to the $\,M_{\rm BH}$ - $\,M_{*}$ relation. In the upper panel of Figure 18, we identify all the BHs that have undergone at least 1 merger event and plot their host galaxy mass in the $\,M_{\rm BH}$ - $\,M_{*}$ plane. The cross markers are coloured and scaled by the number of the BH progenitors. The black dashed-dotted lines are a repeat of those shown in Figure 9, giving the average $\,M_{\rm BH}$ - $\,M_{*}$ relation from $z=6$ to $z=3$ in Astrid . It is interesting to see that the cross markers evenly populate the area around the mean $\,M_{\rm BH}$ - $\,M_{*}$ relation, indicating that the systems hosting BH mergers do not have a biased distribution compared to the overall population. We note that some other simulations that perform BH mergers without modelling dynamical friction also find a similar feature (see DeGraf et al., prep).