Introducing the BRAHMA simulation suite: Signatures of low mass black hole seeding models in cosmological simulations

The first"seeds"of supermassive black holes (BH) can range from $\sim10^2-10^6~M_{\odot}$. However, the lowest mass seeds ($\lesssim10^3 M_{\odot}$) are inaccessible to most cosmological simulations due to resolution limitations. We present our new BRAHMA suite of cosmological simulations that uses a novel flexible seeding approach to represent low mass seeds. Our suite consists of two types of boxes that model $\sim10^3~M_{\odot}$ seeds using two distinct but mutually consistent seeding prescriptions at different simulation resolutions. First, we have the highest resolution $[9~\mathrm{Mpc}]^3$ (BRAHMA-9-D3) boxes that directly resolve $\sim10^3~M_{\odot}$ seeds and place them within halos with dense and metal poor gas. Second, we have lower-resolution and larger-volume $[18~\mathrm{Mpc}]^3$ (BRAHMA-18-E4) and $\sim[36~\mathrm{Mpc}]^3$ (BRAHMA-36-E5) boxes that seed their smallest resolvable $\sim10^4~\&~10^5~\mathrm{M_{\odot}}$ BH descendants using new stochastic seeding prescriptions calibrated using the BRAHMA-9-D3 results. The three boxes together probe BHs between $\sim10^3-10^7 M_{\odot}$ at $z>7$ and we predict their key observables. The variation in the AGN luminosity functions is small (factors of $\sim2-3$) at the anticipated detection limits of potential future X-ray facilities ($\sim10^{43} \mathrm{ergs~s^{-1}}$ at $z\sim7$). Our simulations predict BHs $\sim10-100$ times heavier than expectations from local $M_*$ vs $M_{bh}$ relations, consistent with several JWST-detected AGN. For different seed models, our simulations merge BH binaries at $\sim1-15~\mathrm{kpc}$, with rates of $\sim200-2000$ per year for $\gtrsim10^3 M_{\odot}$ BHs, $\sim6-60$ per year for $\gtrsim10^4~M_{\odot}$ BHs, and up to $\sim10$ per year amongst $\gtrsim10^5 M_{\odot}$ BHs. These results suggest that the LISA mission has promising prospects for constraining seed models.

Prior to JWST, the observational space of high-z (z ∼ 6 − 7.5) BH populations was solely comprised of the most luminous quasars powered by ∼ 10 9 − 10 10 M⊙ BHs (Fan et al. 2001;Willott et al. 2010;Mortlock et al. 2011;Venemans et al. 2015;Jiang et al. 2016;Bañados et al. 2016;Reed et al. 2017;Matsuoka et al. 2018;Wang et al. 2018;Bañados et al. 2018;Matsuoka et al. 2019;Yang et al. 2019;Wang et al. 2021).Theoretical studies have shown that assembling these objects requires a combination of optimistic scenarios for BH seeding and growth.These include possible contributions from both light and heavy seeds growing via sustained Eddington or super-Eddington accretion (Sijacki et al. 2009;Bennett et al. 2023;Fragione & Pacucci 2023).Given that sustaining such accretion rates might be difficult, their assembly may also need to be aided by adequately efficient merger driven growth at the earliest times (z ≳ 10, Bhowmick et al. 2022b).
With JWST, we have now opened up a new frontier for high-z AGN science with spectroscopic confirmations of moderately bright (bolometric luminosities ∼ 10 45 erg s −1 ) AGN powered by ∼ 10 6 − 10 8 M⊙ BHs between z ∼ 4 − 11 (Larson et al. 2023;Harikane et al. 2023;Matthee et al. 2023;Goulding et al. 2023;Maiolino et al. 2023a,b).The mere existence of objects such as CEERS-1019 (Larson et al. 2023) and GN-z11 (Maiolino et al. 2023a) at redshifts as high as z ∼ 10 − 11, could pose a significant constraint on BH seeding and growth models.Furthermore, several of the z ∼ 4−11 JWST AGN are ∼ 10 − 100 times more massive compared to the expectations from local scaling relations (Pacucci et al. 2023;Goulding et al. 2023), which has been interpreted to be more favorable to heavy seed scenarios (Natarajan et al. 2023;Pacucci et al. 2023).Finally, since these recently discovered JWST AGN are likely powered by BHs much less massive than those hosted by the most luminous pre-JWST quasars, they may help resolve potential degeneracies between constraints on BH growth and BH seeding.
However, the strongest possible constraints on BH seeding can only be obtained via direct observations of the seed populations; i.e. before they start to grow via gas accretion and mergers.The vast majority of these postulated seeds are presumed to lie within the largely unexplored regime of "intermediate mass" black holes or IMBHs (∼ 10 2 − 10 6 M⊙).Detecting these IMBH populations is going to be crucial for constraining the seeding mechanisms (see review by Greene et al. 2020).While the more massive ∼ 10 5 M⊙ seeds lie at the very edge of the detection limits for JWST (Natarajan et al. 2017;Cann et al. 2018;Inayoshi et al. 2022), potential future X-ray missions using the proposed NASA X-ray probes for the Astrophysics Probe Explorer (APEX) programme, are likely going to detect them more readily 1 .
Detecting the lower mass end of the postulated seed population (≲ 10 4 M⊙) at z ≳ 7 is likely going to be infeasible even for the next generation of electromagnetic (EM) facilities.Here, gravitational waves (GW) offer a new window as they are able to probe merging BH binaries regardless of whether they are actively accreting.The Laser Interferometer Gravitational-Wave Observatory (LIGO; Abbott et al. 2009) took the first steps towards probing the elu-sive IMBH regime (∼ 10 2 − 10 5 M⊙) with the detection of GW190521 (Abbott et al. 2020), a merger which produced a ∼ 142 M⊙ BH remnant.We also note the recent remarkable achievement of the North American Nanohertz Observatory for Gravitational Waves (NANOGrav, Agazie et al. 2023a) that used pulsar timing arrays (PTA) to find evidence for the Hellings-Downs correlation expected from a stochastic GW background.This was further corroborated by other collaborations, namely the combined European and Indian PTA (EPTA + InPTA, Antoniadis et al. 2023), and the Chinese PTA (CPTA, Xu et al. 2023).This is likely the first observational evidence showing that SMBHs do merge (Agazie et al. 2023b).Notably, the measured GW background amplitude is somewhat higher than expectations from local scaling relations, further highlighting that they offer a brand new complementary probe to the EM observations of SMBH populations (Sato-Polito et al. 2023).While PTAs are most sensitive to high mass (∼ 10 9 M⊙) SMBH binaries, the upcoming Laser Interferometer Space Antenna (LISA; Baker et al. 2019) is expected to detect the bulk of the merging binaries in the IMBH regime, including BHs as small as ∼ 10 3 M⊙ merging at redshifts as high as z ∼ 15 (Amaro-Seoane et al.

2017).
Robust theoretical predictions for IMBH populations are needed to constrain seed models from the wealth of upcoming observational data that we expect over the next two decades.To date, empirical and semi-analytic models (SAMs) have been extensively used to the probe the impact of BH seeding on IMBH populations (e.g., Sesana et al. 2007;Volonteri & Natarajan 2009;Barausse 2012;Valiante et al. 2018;Ricarte & Natarajan 2018;Dayal et al. 2019;Sassano et al. 2021;Evans et al. 2023).With the advantage of low computational cost, these models have explored a wide range of seeding scenarios both in terms physical channels (Pop III, NSC and DCBH) as well as in their technical implementation.They generally predict strong signatures of seeding within IMBH populations, particularly within the merger rates measurable by LISA.However, these models do come with the drawback of not being able to directly probe the gas dynamics that is naturally expected to play a crucial role in BH seed formation; this is only possible in full cosmological hydrodynamic simulations In cosmological hydrodynamic simulations, modeling the formation of seeds and the resulting IMBH populations has been particularly challenging due to resolution and dynamic range limitations.Most uniform volume cosmological simulations (Di Matteo et al. 2012;Vogelsberger et al. 2014b;Sijacki et al. 2015;Khandai et al. 2015;Schaye et al. 2015;Volonteri et al. 2016;Dubois et al. 2016;Kaviraj et al. 2017;Tremmel et al. 2017;Nelson et al. 2019;Volonteri et al. 2020;Vogelsberger et al. 2020) have gas mass resolutions ranging between ∼ 10 5 − 10 7 M⊙.This naturally sets the minimum resolvable BH mass, which leaves the bulk of the IMBH populations (∼ 10 2 − 10 5 M⊙) unresolved.
Many cosmological simulations simply seed ∼ 10 5 −10 6 M⊙ BHs in sufficiently massive (≳ 10 10 M⊙) halos (e.g.Di Matteo et al. 2012;Vogelsberger et al. 2014b;Khandai et al. 2015;Nelson et al. 2019).This seeding prescription is consistent with our expection (based on observations) that all sufficiently massive galaxies should eventually contain a BH through some (currently unknown) seeding mechanism.These simulations do broadly reproduce the observed local supermassive BH populations (Habouzit et al. 2020) wherein the imprints of seeding are expected to get washed out.However, the "halo mass threshold" based approach to seeding cannot distinguish between the different theoretical seeding channels that have been postulated.This is a significant liability for probing IMBH as well as highest-z SMBH populations that are more likely to be sensitive to the specifics of their seeding origins.For simulations at the upper end of the typical range of gas mass resolutions i.e. ∼ 10 5 M⊙, it is still possible to emulate the heavy seed channels such as DCBHs.This is in fact done inTremmel et al. (2017) and Bellovary et al. (2019) wherein 10 5 M⊙ seeds are produced from gas cells that are dense, metal-poor (< 3 × 10 −4 Z⊙) and at temperatures of ∼ 9500 − 10000 K (ideal conditions for direct collapse).
However, simulating low mass seeds (≲ 10 3 M⊙) presents an additional challenge compared to their heavier counterparts like DCBHs.This is due to the fact that directly resolving them within our desired volumes would render the computations prohibitively expensive.The simulations that do attempt to model seed masses as low as ∼ 10 4 M⊙ (Ni et al. 2022) or ∼ 10 3 M⊙ (Taylor & Kobayashi 2014;Habouzit et al. 2017;Wang et al. 2019), do so without explicitly resolving the seed-forming gas to those masses.In our recent paper (Bhowmick et al. 2023, hereafter Paper I), we built a new sub-grid stochastic seed model that can faithfully represent the smallest resolvable (∼ 10 4 − 10 5 M⊙) descendants of ∼ 10 3 M⊙ seeds in lower resolution simulations.More specifically, we demonstrated that with a combination of stochastic seeding criteria based on galaxy total mass (including dark matter, stars and gas) and an environmental richness parameter (number of neighboring halos), we could represent seeds that are ∼ 10 − 100 times smaller than the masses of the individual gas cells, without having to directly resolve these seeds.
In this paper, we build on the results of Paper I to create a new suite of uniform volume cosmological hydrodynamic simulations: the BRAHMA simulations.Note that we choose not to use the exact seed model calibrations derived from the highly overdense zoom regions of Paper I in this work, due to cosmic variance.Instead, we re-calibrate our seed models using uniform volume simulations.Our primary simulations comprise of a series of boxes with volumes being successively increased (and mass resolutions decreased) by factors of 8 i.e. (9 Mpc) 3 , (18 Mpc) 3 and (36 Mpc) 3 .The smallest of these boxes directly resolve the ∼ 10 3 M⊙ seeds.In the remaining larger boxes, the higher mass descendants are seeded at the gas mass resolutions using the newly built stochastic seed model.We combine these three boxes to provide predictions for z ≥ 7 IMBH and SMBH populations for a range of seed models.The three boxes together provide predictions for BHs with masses ranging from ∼ 10 3 − 10 7 M⊙.We note here that we do not attempt to resolve the lightest Pop III seed masses (∼ 100 M⊙) as our underlying galaxy formation model is not built for simulating individual Pop III stars.Nevertheless, our target seed mass of ∼ 10 3 M⊙ is still at the lower end of the postulated seed mass spectrum.
Section 2 presents the basic simulation setup and the underlying galaxy formation model.Section 3 presents the details of the simulation suite.Section 4 describes the BH seed models that we use for the different boxes.Section 5 presents the detailed results followed by the key conclusions in section 6.

Identifying bound structures and substructures
We identify halos using the friends of friends (FOF) algorithm (Davis et al. 1985) with a linking length of 0.2 times the mean particle separation.Subhalos are computed using the SUBFIND (Springel et al. 2001) algorithm for each simulation snapshot.In addition to the standard halos and subhalos computed in most simulations, we also compute an additional catalog of FOFs using a factor of 3 shorter linking length compared to that used for halos.These are referred to as "best-friends-of-friends" or bFOFs.We studied their properties in relation to the FOF halos in Paper I. Here, we use them as seeding sites for our stochastic seed model as detailed in Section 4.2.
Note that both bFOFs and SUBFIND-based subhalos identify bound substructures within halos.The baryonic components of these substructures correspond to "galaxies".While it is more common to refer to SUBFIND-based subhalos as "galaxies", here we use bFOFs very extensively to describe our seed models.Therefore, similar to Paper I, we use the term "galaxies" to refer to bFOFs.In select few instances where we do refer to the SUBFIND-based subhalos, we shall use the term "subfind-subhalos" or "subfind-galaxies".

Star formation, chemical enrichment and feedback
The following sub-grid models are used to capture the cooling of gas to densities high enough to form stars, and the subsequent stellar evolution and chemical enrichment of gas.
• The radiative cooling processes include contributions from primodial species (H, H + , He, He + , He ++ based on Katz et al. 1996), as well as metals (with cooling rates interpolated from pre-calculated tables as in Smith et al. 2008) in the presence of a spatially uniform, time dependent UV background.
• Star formation occurs in gas cells with densities exceeding 0.1 cm −3 , with an associated time scale of 2.2 Gyr.These gas cells represent an unresolved multiphase interstellar medium described by an effective equation of state (Springel & Hernquist 2003;Vogelsberger et al. 2014a).
• The model for stellar evolution is adopted from Vogelsberger et al. ( 2013) with modifications for IllustrisTNG as described in Pillepich et al. (2018a).Star particles represent an unresolved single stellar population with fixed age and metallicity, with an underlying initial mass function adopted from Chabrier (2003).
• The chemical enrichment of gas occurs by following the evolution of seven species of metals (C, N, O, Ne, Mg, Si, Fe) in addition to H and He.Prior to enrichment caused by stellar evolution, the gas is assigned an initial metallicity of 7 × 10 −8 Z⊙.More details in Pillepich et al. (2018a).
• Stellar and Type Ia/II Supernova feedback are modelled as galactic scale winds (Pillepich et al. 2018b) that deposit mass, momentum and metals on to the gas surrounding the star particles.

BH accretion, feedback and dynamics
We now describe the sub-grid models implemented for all aspects of BH physics aside from their initial seeding • BH accretion in IllustrisTNG is modelled by the Eddington-limited Bondi-Hoyle formalism given by Ṁbh = min( ṀBondi , ṀEdd ) (1) where G is the gravitational constant, ρ is the local gas density, M bh is the BH mass, cs is the local sound speed, mp is the proton mass, and σT is the Thompson scattering cross section.Accreting BHs are assumed to radiate at bolometric luminosities given by with an assumed radiative efficiency of ϵr = 0.2.
• IllustrisTNG implements a 'thermal feedback' for Eddington ratios ( Ṁbh / ṀEdd ) higher than a critical value (given by ηcrit = min[0.002(MBH/10 8M⊙) 2 , 0.1]).Here, a fraction of the radiated luminosity is deposited to the neighboring gas at a rate of ϵ f,high ϵr ṀBHc 2 with ϵ f,high ϵr = 0.02 where ϵ f,high is the "high accretion state" coupling efficiency.For Eddington ratios lower than the critical value, feedback is in the form of kinetic energy that is injected onto the gas at irregular time intervals along a randomly chosen direction.The rate of injection is given by ϵ f,low ṀBHc 2 with ϵ f,low being the 'low accretion state' coupling efficiency (ϵ f,low ≲ 0.2).Readers interested in further details are encouraged to refer to Weinberger et al. (2017).
• The limited mass resolution prevents our simulations from fully resolving the BH dynamical friction force that is crucial for capturing the small-scale BH dynamics, particularly at their lowest masses.In fact, the DM particles are ∼ 10 times more massive than the BH seeds; this makes the seeds susceptible to encountering spuriously large kicks.To stabilize the dynamics, a BH is "re-positioned" to the nearest potential minimum within its "neighborhood" (defined by n gas ngb nearest neighboring gas cells).Two BHs are promptly merged when at least one of them is within the "neighbor search radius" (R Hsml ) of the other.The target neighbor counts (n gas ngb ) are listed in Table 1 for the different simulation boxes, and these are essentially set to have R Hsml values roughly similar to that in the Illustris-TNG suite.As we shall see in Section 5.5, the resulting R Hsml values range from ∼ 1 − 15 kpc at the time of their merger in the simulation.This essentially implies that BHs are promptly merged at separations ranging from ∼ 1 − 15 kpc.
As we build our seed models over the foundation of the success of the Illustris-TNG model in reproducing the observed low-z galaxy populations, we must also note that Illustris-TNG and its successors THESAN (Kannan et al. 2022) and Millennium-TNG (Pakmor et al. 2023) tend to underpredict the abundances of the z ≳ 12 galaxy populations observed by JWST (Kannan et al. 2023).This could be either due to missing physics (for example, weak stellar feedback at high-z, more top heavy initial mass function of high-z stellar populations, or lack of stochasticity in star formation due to an unresolved ISM) combined with insufficient resolution.In fact, the latter may significantly contribute to this tension given that lower resolutions delay the onset of star formation as we demonstrated in Figure 19 of Bhowmick et al. (2021).In any case, since BH seed formation is deeply intertwined with the star formation and stellar evolution physics, resolving this puzzle could have significant implications on our results, which we plan to explore in the future.

THE BRAHMA SIMULATION SUITE
Our primary simulation suite comprises of three boxes, namely BRAHMA-9-D3, BRAHMA-18-E4 and BRAHMA-36-E5.Each of these boxes has 1024 3 DM particles with initial conditions (IC) generated using MUSIC (Hahn & Abel 2011); as a result, their DM mass resolutions differ by factors of ∼ 8.We named our simulations as BRAHMA-*-D* / BRAHMA-*-E*, wherein the first '*' represents the box-length in comoving Mpc units and the second '*' represents the logarithm of the typical mass resolution of the gas cells and the BHs in M⊙ units.The letters 'D' and 'E' refer to whether the seeds are "direct gas based seeds (DGBs)" or "extended seed descendants (ESDs)", as defined in the next section.Accordingly, the box lengths of BRAHMA-9-D3, BRAHMA-18-E4 and BRAHMA-36-E5 are 9, 18 & 36 Mpc respectively, with gas mass resolutions of ∼ 10 3 , ∼ 10 4 and ∼ 10 5 M⊙ respectively.These primary set of boxes are run until z = 7. Figure 1 shows a visualization of the gas density field for the three boxes at z = 7.As illustrated in the visualizations, the gas mass resolution also sets the miminum (seed) BH mass that can be resolved in these simulations, which are chosen to be 2.3×10 3 , 1.8 × 10 4 and 1.5 × 10 5 M⊙ for BRAHMA-9-D3, BRAHMA-18-E4 and BRAHMA-36-E5 respectively (see more details in Section 4).The specific values for the seed masses were chosen by reducing the TNG seed mass of 8×10 5 M⊙/h (= 1.2×10 6 M⊙) by factors of 8.
We also have some lower resolution versions of the N dm and M dm are the number of DM particles in the simulation and the mass of each DM particle (3rd and 4th columns), Mgas is the typical mass of a gas cell (5th column: note that gas cells can refine and de-refine depending on the local density), and ϵ is the gravitational smoothing length (6th column).The 7th and 8th columns correspond to the seed mass and seed type used at the different resolutions.The 9th column shows the final redshift z final to which the simulations were run.The entire simulation suite is divided into two categories (primary and secondary), each with a unique purpose as stated in the Table.9 & 18 Mpc boxes with 512 3 DM particles.These are referred to as the BRAHMA-9-E4 and BRAHMA-18-E5 simulations, and have gas mass resolutions of ∼ 10 4 M⊙ and ∼ 10 5 M⊙ respectively.The lower count of DM particles and gas cells substantially reduces the computational expense of these boxes, which allows us to do two things: 1) Run a large number of boxes with different ICs to assess the cosmic variance in our predictions in Section 5.1.2) Run these boxes to low redshifts (z ∼ 0) and validate them against available observational constraints.Note that when we run the secondary boxes below z = 7, there is a caveat that the calibration of the BH seed models (discussed in the next section) in the lower resolution boxes is done based on the BRAHMA-9-D3 boxes, which are only run until z = 7.Another potential issue with running these relatively smaller volumes to z = 0 is the missing large scale modes in the overall structure growth, which are expected to become non-linear at these late times.Due to these reasons, we primarily focus on the z ≥ 7 BH populations in this paper, even though our simulations do also make predictions at lower redshifts.Nevertheless, it is still a valuable exercise to show that despite these caveats, all of our seed models do a good job at reproducing a variety of observed measurements for the low-z BH populations.In addition, we will also show that our predicted low-z BH populations are consistent with larger volume TNG simulations that do contain some of the large scale modes missing within the BRAHMA boxes.Lastly, we have some extremely small 4.5 Mpc boxes at the highest ∼ 10 3 M⊙ gas mass resolution, which are used to test the impact of removing the BH repositioning scheme on the BH merger rates (Appendix B).
Table 1 shows the full details of our simulation suite.Notably, these are only the first set of boxes that will be a small portion of the full BRAHMA simulation suite.At its completion, the full BRAHMA suite will include a much larger set of simulations that span a wider range of BH physics models including heavy DCBH seeds, alternate treatment of BH dynamics, as well as various accretion and feedback prescriptions.

BLACK HOLE SEED MODELS
As we mentioned, the key distinguishing feature of the BRAHMA simulation suite is our novel approach towards modeling BH formation, which we describe in this section.In this paper, we focus exclusively on making predictions for low mass seeds.Note that while the postulated masses of Pop III seeds can be as low as ∼ 100 M⊙, explicitly resolving ∼ 100 M⊙ objects would be extremely computationally expensive within our simulation volumes.Moreover, our underlying galaxy formation model was never designed to explicitly resolve the origin of individual Pop III stars and their remnant BHs.Due to all these considerations, we choose our target seed mass to be 2.3 × 10 3 M⊙, which we explicitly resolve in our smallest, highest resolution box i.e.BRAHMA-9-D3.These could potentially represent the most massive Pop III seeds, or relatively low mass NSC seeds.However, we remain agnostic about which specific theoretical seeding channel do our models represent.For the BRAHMA-18-E4 and BRAHMA-36-E5 boxes that cannot explicitly resolve our target 2.3 × 10 3 M⊙ seeds, we seed their smallest resolvable descendants, particularly at masses of 1.8 × 10 4 M⊙ and 1.5 × 10 5 M⊙ respectively.The details of the seeding criteria used in the different boxes follow.
4.1 When ∼ 10 3 M⊙ seeds are explicitly resolvable: Gas based seed model In the BRAHMA-9-D3 box, we place the explicitly resolved 2.3 × 10 3 M⊙ seeds in halos with pristine dense (star forming) gas, as illustrated in Figure 2. In order to distinguish them from the "unresolved" seeds representing their descendants within the larger volumes, we shall hereafter refer the 2.3×10 3 M⊙ seeds as "direct gas based" seeds or DGBs (with mass M DGB seed ).We refer to this seed model as the "gas-based seed model", which was developed, tested and systematically explored in our previous works (Bhowmick et al. 2021(Bhowmick et al. , 2022a,b) ,b) across a wide range of seed parameters and overdensities.In Bhowmick et al. (2021), we demonstrated that these seed models are reasonably well converged (to within ∼ 25% for the seed formation rates) at our adopted gas mass resolution of ∼ 10 3 M⊙.We seed the DGBs with the following two criteria • Dense & metal poor gas mass criterion: DGBs are placed in halos with a minimum threshold of dense (> 0.1 cm −3 i.e. the star formation threshold) & metal poor (Z < 10 −4 Z⊙) gas mass, denoted by Msfmp (in the units of M DGB seed ).Msfmp measures the efficiency with which pristine dense gas is converted to BH seeds.There are no empirical or theoretical constraints on Msfmp .In this work, we consider models with Msfmp = 5, 50, 150 & 1000.
• Halo mass criterion: We also ensure that DGBs can only form in those halos where the total mass exceeds a critical threshold, specified by Mh in the units of M DGB seed .Here we consider Mh values of 1000, 3000 & 10000.While DGB formation should only depend on the gas properties, we still adopt this criterion for two reasons.First, we want to avoid seeding in halos significantly below the atomic cooling threshold.This is because in our simulations, these (mini)halos are likely to have artificially suppressed star formation due to the absence of H2 cooling that is required to self-consistently capture the collapse of gas.Second, this criterion favors seed formation to occur within deeper gravitational potentials.This is consistent with the expectation that NSC seeds will grow more efficiently within sufficiently deep gravitational potential wells where runaway BH merger remnants are unable to easily escape the cluster.
The above two criteria dictate that the formation of DGBs is essentially driven by the interplay of three main processes namely, halo growth, dense gas formation, and metal enrichment.Generically speaking, at the highest redshifts when the gas is nearly pristine, DGB formation is primarily driven by either halo growth or dense gas formation.Over time, metal enrichment takes over and suppresses the formation of new DGBs.In Paper I, we discuss how the gas based seed parameters Mh and Msfmp determine which amongst the three processes dominates DGB formation (or suppression) at different epochs.We also demonstrated in Paper I that this interplay between the three processes is not just imprinted on the formation of DGBs, but also in the assembly of higher mass descendants.In the BRAHMA-18-E4 and BRAHMA-36-E5 simulations that can only resolve the higher mass descendants of the DGBs, we specifically try to capture this interplay using the seed models discussed in the next subsection.
4.2 When ∼ 10 3 M⊙ seeds cannot be explicitly resolved: Stochastic seed model The gas based seed models described in the previous subsection cannot be applied directly to the primary BRAHMA-18-E4 and BRAHMA-36-E5 boxes and the secondary BRAHMA-9-E4 and BRAHMA-18-E5 boxes, due to their inability to directly resolve the 2.3 × 10 3 M⊙ DGBs.As already mentioned, this is a longstanding limitation of larger volume cosmological simulations.Paper I was specifically aimed at addressing this, wherein we developed a new stochastic seed model that seeds the smallest resolvable descendants of the 2.3×10 3 M⊙ DGBs.
To distinguish these seeded descendants against their progenitor DGBs, we shall refer to them as "extrapolated seed descendants" or ESDs.We assume ESD masses (denoted by M ESD seed ) of 1.8×10 4 M⊙ and 1.5×10 5 M⊙ for BRAHMA-*-E4 and BRAHMA-*-E5 boxes respectively.These assumed ESD masses are chosen to be comparable to the gas mass resolutions of these larger volume boxes.In Paper I, we thoroughly discussed the motivation and construction of the various components of this model using a large suite of zoom simulations.Here, we summarize the key features: In this seed model, ESDs are stochastically placed in BRAHMA-18-E4 and BRAHMA-36-E5 galaxies (see Section 2.1) based on the galaxies wherein the descendants of the 2.3 × 10 3 M⊙ DGBs end up within the BRAHMA-9-D3 simulations.Note that unlike the gas-based seed model that places one  To obtain fits, we first assumed appropriate ztrans values for the different seed models via visual inspection.The fits were then performed to obtain the slopes at z < ztrans and z > ztrans using scipy.optimize.curve_fit.The final fitted parameters are shown in Table 2.These fits are used in our stochastic seed models that directly seeds the descendants at 1.8 × 10 4 M ⊙ and 1.5 × 10 5 M ⊙ masses within the BRAHMA-18-E4 and BRAHMA-36-E5 boxes respectively.To distinguish them from the DGBs, the seeds formed by our stochastic seed models are referred to as extrapolated seed descandants or ESDs.
DGB per halo, we do not place one ESD per halo.This is because the small halos that form the DGBs eventually merge together, such that the higher mass descendants occur in higher multiplicities within individual halos.To capture this, we seed one ESD per galaxy (recall that galaxies are bFOFs that trace halo substructures), which naturally leads to multiple ESDs per halo (see Paper I for further details).In the following subsections, we describe the two main stochastic seeding criteria that we use.

Galaxy mass criterion
In Paper I, we demonstrated that despite the 2.3 × 10 3 M⊙ DGBs forming in halos biased towards low metallicities, their higher mass descendants end up within unbiased galaxies that are fully characterized by their total mass.This is because the halos forming the DGBs are in a transient state of rapid metal enrichment, which can be clearly seen in the rightmost panel of Figure 2. Therefore, in BRAHMA-*-E4 and BRAHMA-*-E5, ESDs are seeded in galaxies that exceed a threshold total mass (including dark matter, stars and gas).The placement of these ESDs is meant to emulate the assembly of 1.8 × 10 4 M⊙ and 1.5 × 10 5 M⊙ descendants respectively.Unlike many cosmological simulations that use a mass threshold that is constant and uniform across the entire simulation volume (for e.g.Vogelsberger et al. 2014a;Khandai et al. 2015;Nelson et al. 2018;Ni et al. 2022), our mass thresholds are stochastically drawn from galaxy mass distributions predicted for the assembly of 1.8 × 10 4 M⊙ and 1.5 × 10 5 M⊙ descendants (denoted by M galaxy total ) within the BRAHMA-9-D3 boxes.We compute M galaxy total by tracing BH growth along galaxy merger trees (see Paper I for the detailed methodology).As we can see in Figure 3, these M galaxy total distributions exhibit significant scatter and have a redshift dependence that is distinct for the different gas based seed parameters.The redshift dependence of M galaxy total is driven by a complex interplay of various processes, namely halo growth, dense gas formation and metal enrichment.
To implement these distributions within our stochastic seed models, we model them as a log-normal function i.e where M th is the mass threshold.µ ≡ log 10 M galaxy total (z) is the redshift dependent mean of the distributions and σ quantifies the scatter.The redshift evolution of the mean galaxy mass threshold M galaxy total (z) is described by a double power law: wherein ztrans roughly corresponds to when there is a change in the driving physical process for DGB formation.At z > ztrans DGB formation is assumed to be driven by halo growth or dense gas formation; at z < ztrans, metal enrichment takes over as the primary driver to suppress DGB formation.Mtrans is the mean of the distribution at z = ztrans.
Overall, the galaxy mass criterion is characterised by the five parameters ztrans, Mtrans, α, β and σ.The first four of these parameters are determined by fitting the M galaxy total vs z relation to Eq. 5; these are shown as dashed lines in Figure 3.The fifth parameter σ is assumed to be redshift independent based on the shaded regions of Figure 3; its values are determined by computing the standard deviations of M galaxy total within each redshift bin, and marginalizing over all the bins.Table 2 shows the best fit values for ztrans, Mtrans, α, β, and σ for each set of gas based seed parameters Mh and Msfmp .These values are used for implementing the galaxy mass criterion.Finally, there may be seed models for which it is not possible to obtain robust fits for β because they only start producing enough descendants at later times (z ≳ 12) .We have one such case, namely the assembly of 1.5 × 10 5 M⊙ BHs for [ Mh , Msfmp ] = [3000,5].In this case, we simply assume β = 0.

Galaxy environment criterion
We also found in Paper I that applying only a galaxy mass criterion leads to an underestimation of the small-scale clustering of ESDs compared to the actual descendants of 2.3 × 10 3 M⊙ DGBs.This is essentially because the merger dominated growth of DGBs accelerates the assembly of descendants in galaxies living in rich environments with more extensive merger histories.To ensure that the seeded ESDs faithfully represent the actual descendants of 2.3 × 10 3 M⊙ DGBs, we apply a galaxy environment criterion that preferentially places these ESDs in rich environments.Here we define the environment of a galaxy as the number of neighboring halos (N ngb ) that exceed the mass of its host halo and are located within 5 times the virial radius (Rvir).
The galaxy environment criterion is implemented by adding a seed probability (P env seed ) less than unity to suppress ESD formation in galaxies with ≤ 1 neighboring halos.P env seed is assigned to linearly increase with the galaxy mass as (6) where p0 and p1 denote the seeding probability in galaxies with 0 and 1 neighbors respectively, at the mean M galaxy total of the total mass distributions applied within the galaxy mass criterion.γ > 0 ensures that for a given N ngb , the seed probability increases with galaxy mass.This is because a more massive galaxy, accompanied by N ngb massive neighboring halos, lives in a more extreme environment compared to a lower mass galaxy with an identical count of N ngb halos of lesser mass.
Our galaxy environment criterion contains three parameters namely p0 < 1, p1 < 1 and γ > 0. By applying suitable choices for these parameters, we can ensure that the ESDs preferentially form in rich environments, examples of which are shown in Figure 4.Note also that these ESD seeding sites are fully metal enriched (right panels of Figure 4), which is expected from the fact that their parent DGBs formed in halos undergoing rapid metal enrichment.In Paper I, we systematically explored a wide range of values for these parameters, and found that a choice of p0 ∼ 0.1 − 0.2 and p1 ∼ 0.3 − 0.4 with γ = 1.6 reproduces the correct small scale clustering as well as overall abundances for the ESDs.In this paper, we adopt p0 = 0.1, p1 = 0.3 & γ = 1.6.

Inclusion of unresolved minor mergers to the BH growth
In Bhowmick et al. (2021), we determined that at these highest redshifts and lowest BH masses, the Bondi accretion rates are very small and the BH growth is primarily driven by BH mergers.However, our stochastic seed models do not capture (by construction) the growth due to those mergers where the primary BH is resolved (M1 > M ESD seed ) but the secondary BH is unresolved (M DGB seed < M2 < M ESD seed ).We refer to them as "light minor mergers", and contrast it against the "heavy mergers" where both primary and secondary BHs are resolved (M1 > M2 > M ESD seed ).In Paper I, we demonstrated that these light minor mergers make a substantial contribution to the BH growth.Therefore, we explicitly add their contribution in the BRAHMA-*-E4 and BRAHMA-*-E5 boxes as an extra mass growth factor.Based on the results from BRAHMA-9-D3 (shown in Figure A1), we add a mass growth of 4 M ESD seed during every resolved heavy merger event for all these boxes using the stochastic seed model.For further details and validation for this choice, please refer to Appendix A.

Summarizing the construction of the stochastic seed models
Having discussed all the components of our stochastic seed models, we now summarize the key steps involved in their construction: (i) We run the BRAHMA-9-D3 boxes for the different combinations of the gas-based seed parameters Mh and Msfmp .These boxes directly resolve our target 2.3 × 10 3 M⊙ DGBs.
(ii) For each combination of Mh and Msfmp , we determine the total masses (stars, gas and DM) of galaxies where 1.8 × 10 4 M⊙ and 1.5 × 10 5 M⊙ descendants assemble in BRAHMA-9-D3 (shown in Figure 3).Based on the results, we calibrate the parameters for the galaxy mass criterion (ztrans, Mtrans, α, β, and σ).(iii) The parameters for the galaxy environment criterion (p0, p1 and γ) are neither very sensitive to the gasbased seed parameters nor the overdensity of the simulated region.Therefore, we simply assume the parameters values that were derived from the zoom simulations of Paper I i.e. p0 = 0.1, p1 = 0.3 & γ = 1.6.
(iv) To include the contributions from unresolved light minor mergers to the growth of our ESDs in the stochastic seed models, we first used the BRAHMA-9-D3 boxes to determine their contribution to the BH growth (shown in Figure A1).Based on the results, we assumed that 4 M ESD seed of additional mass growth is contributed for every resolved heavy merger.

Seed model nomenclature
Our gas based seed models applied to BRAHMA-9-D3 are characterised by three parameters, namely Msfmp , Mh and M DGB seed .For each of these seed parameters, we calibrate the parameters for the stochastic seed model that represents the underlying gas based seed model within the lower resolution BRAHMA-*-E4 and BRAHMA-*-E5 boxes.Therefore, for both In these boxes, we apply the stochastic seed model with M ESD seed = 1.8 × 10 4 M ⊙ .We find that for the 9 Mpc boxes, the number densities can vary up to factors of ≲ 2 at z ≲ 15 due to cosmic variance.We expect this variance to be smaller in the larger volume boxes.
gas based seed model as well as the corresponding calibrated stochastic seed model, we use a common label of SM*_FOF* where the '*'s correspond to the values of Msfmp and Mh respectively.As an example, Msfmp = 5 and Mh = 3000 will be labelled as SM5_FOF3000.

Stochastic seed model validation
We start by validating the capability of our stochastic seed model to adequately represent the descendants of 2.3 × 10 3 M⊙ DGBs at lower resolutions.In the left panel of Figure 5, the solid lines show the BRAHMA-9-D3 predictions for the evolution of the galaxy total masses wherein the 1.8 × 10 4 M⊙ (top) and 1.5 × 10 5 M⊙ (bottom) descendants assemble from 2.3 × 10 3 M⊙ DGBs.The dashed lines show the same quantity for the BRAHMA-18-E4 (top panel) and BRAHMA-36-E5 (bottom panel) boxes, wherein 1.8 × 10 4 M⊙ and 1.5 × 10 5 M⊙ ESDs are placed using the stochastic seed model with the calibrated parameters listed in Table 2.The reasonably good match between the solid versus dashed lines depicts the successful calibration of the galaxy mass criterion.In Paper I, we demonstrated that once we calibrate for the total galaxy mass, the ESDs also naturally get placed in galaxies with the correct baryonic properties such as stellar mass, star formation rates and metallicities.
Recall again that not all seed models can be represented in the BRAHMA-18-E4 and BRAHMA-36-E5 boxes, due to which there are missing dashed lines in Figure 5.For the least restrictive SM5_FOF1000 model, a substantial fraction of the 1.8 × 10 4 M⊙ descendants end up assembling in galaxies that are below the resolution limit of BRAHMA-18-E4.Due to this, we do not attempt to represent the SM5_FOF1000 model in BRAHMA-18-E4.At the other end, the most restrictive SM150_FOF3000 and SM5_FOF10000 produce too few 1.5× 10 5 M⊙ descendants in BRAHMA-9-D3, to perform a robust calibration of the galaxy mass criterion for BRAHMA-36-E5.As a result, we refrain from representing SM150_FOF3000 and SM5_FOF10000 models in BRAHMA-36-E5.
In the right panel of Figure 5, we compare the two point clustering predictions for the ≥ 1.8 × 10 4 M⊙ and ≥ 1.5 × 10 5 M⊙ BHs, between the different boxes.We can clearly see that with our choice of parameters for the galaxy environment criterion i.e. p0 = 0.1, p1 = 0.3 & γ = 1.6, ESDs placed within BRAHMA-18-E4 and BRAHMA-36-E5 are able to successfully capture the two-point clustering of the higher mass descendants predicted by BRAHMA-9-D3.Notably, for the ≥ 1.5 × 10 5 M⊙ BHs, there are some differences between BRAHMA-9-D3 and BRAHMA-36-E5 at the smallest ≲ 0.04 Mpc scales for the SM5_FOF3000 seed model; this is likely due to decreased robustness of the stochastic seed model calibration because of smaller statistics.
In Paper I, we demonstrated using zoom simulations that when our stochastic seed model is calibrated to reproduce the galaxy masses and the small scale clustering of the higher mass descendants of 2.3 × 10 3 M⊙ DGBs, it also produces the correct BH abundances.In Figure 6, we show this again in our uniform boxes for one of our seed models, namely SM5_FOF3000.More specifically, we compare the BRAHMA-9-D3 predictions for the number density of ≥ 1.8 × 10 4 M⊙ BHs (light green line), against that of the secondary BRAHMA-9-E4 boxes that seed 1.8 × 10 4 M⊙ ESDs using our calibrated stochastic seed models.The BRAHMA-9-E4 box that uses the same IC as BRAHMA-9-D3 (black line) produces a very good match.In addition, since we eventually apply the stochastic seed models within volumes larger than BRAHMA-9-D3, we find it imperative to estimate the cosmic variance by running the same seed model on five other different IC realizations of BRAHMA-9-E4 (grey lines of different thicknesses).Amongst these realizations, we see that the variations among the number densities are well within factors of ∼ 2 at z ≲ 15.At the highest redshifts, the variations are slightly higher likely because of higher Poisson error.Overall, Figures 5 and 6 together demonstrate that our stochastic seed models can indeed faithfully represent the descendants (particularly their host galaxy masses, star formation rates, metallicities and their spatial clustering) of 2.3 × 10 3 M⊙ DGBs within BRAHMA-18-E4 and BRAHMA-36-E5 boxes.

RESULTS
Having validated our seed models, we now discuss the key predictions of our simulations in terms of observable properties of z ≥ 7 BH and AGN populations detectable with the ongoing JWST mission as well as proposed facilities such as LISA.vs. redshift relations (Figure 3).ztrans, Mtrans, α, β is obtained by fitting the mean trends using the double power-law function shown in Equation 5. σ is the standard deviation.Columns 9 to 11 show the parameter values for the galaxy environment criterion; i.e. p 0 , p 1 and γ.These were obtained in Paper I by exploring a range of possible values to obtain the best match with the small-scale BH clustering and overall BH counts predicted by the gas-based seed model.These parameter values are used to represent SM5_FOF3000, SM150_FOF3000 and SM5_FOF10000 with the stochastic seed models, within the BRAHMA-18-E4 and BRAHMA-36-E5 boxes.
Figure 7. Evolution of BH number densities for the different seed models.Hereafter, unless otherwise stated, the dark and light green lines are the lenient SM5_FOF1000 and SM5_FOF3000 seed models respectively, which are applied to the BRAHMA-18-E5 boxes that seed 1.5 × 10 5 M ⊙ ESDs.The maroon and pink lines are the strict SM150_FOF3000 and SM5_FOF10000 seed models respectively, applied to the BRAHMA-9-E4 boxes that seed 1.8 × 10 4 M ⊙ ESDs.Likewise, lines of different thicknesses distinguish the BRAHMA-9-D3, BRAHMA-18-E4, BRAHMA-36-E5 boxes.The leftmost panels show BHs above 2.3 × 10 3 M ⊙ , which are only probed by the BRAHMA-9-D3 boxes.The middle panels show BHs ≥ 1.8 × 10 4 M ⊙ , which are probed by both BRAHMA-9-D3 and BRAHMA-18-E4 boxes.The rightmost panels show BHs ≥ 1.5 × 10 5 M ⊙ that can be probed by all three boxes; however, for the BRAHMA-36-E5 box, we could only calibrate the more lenient SM5_FOF1000 and SM5_FOF3000 seed model.This is because for the more restrictive SM150_FOF3000 and SM5_FOF10000 seed models, there are too few 1.5 × 10 5 M ⊙ BHs formed within BRAHMA-9-D3 to perform any robust calibration for BRAHMA-36-E5.The different seed models produce significantly different BH number density evolution.

BH number density evolution
The z ≥ 7 BH number density evolution for the different seed models is shown in Figure 7.The BRAHMA-18-E4 and BRAHMA-36-E5 simulations produce a reasonably good match with the BRAHMA-9-D3 predictions (compare lines of different thickness with the same color).Small differences (up to factors of ∼ 2) owing to cosmic variance are apparent (revisit Figure 6).The larger boxes naturally provide us with much more statistical power and smaller cosmic variance.Additionally, the larger boxes allow us to extend our predictions to higher redshifts and smaller number densities, than in BRAHMA-9-D3.
We first look at the number density evolution of ≥ 2.3 × 10 3 M⊙ BHs (1st panel of Figure 7) probed only by BRAHMA-9-D3, which start to form at z ∼ 25 − 30.For ease of language, we shall hereafter refer to SM5_FOF1000 and SM5_FOF3000 as "lenient" seed models, and SM150_FOF3000 and SM5_FOF10000 as "strict" seed models.This distinction is appropriate given that the lenient seed models produce ≳ 10 times more numerous DGBs compared to the strict seed models.The ranges hereafter quoted in the number densities and various other observables are intended to bracket the variation between these lenient and strict seed models.At z ∼ 25, the number densities of ≥ 2.3 × 10 3 M⊙ BHs range from ∼ 0.01 − 0.5 Mpc −3 .Starting from these highest redshifts, the number densities continue to increase owing to dense gas formation and halo growth until z ∼ 10 − 12 when they reach ∼ 3 − 30 Mpc −3 .After z ∼ 10 − 12 they flatten due to metal enrichment and start to decrease due to BH mergers.
Before looking at the number densities of higher mass BHs Figure 8. Fraction of black hole growth that is contributed from mergers, and its evolution from z = 11 to z = 0; this is plotted as a function of BH mass for all the BHs at a given snapshot.The top and bottom panels are predictions for secondary BRAHMA-9-E4 and BRAHMA-18-E5 boxes respectively.Here we only show the results for one of the seed models (SM5_FOF1000), but our main takeaways are true regardless of the seed model.At z ≳ 7, the BH growth fully dominated by mergers.As we go to lower redshifts, we start to see increased contribution from gas accretion particularly for ≳ 10 6 M ⊙ BHs.By z ≲ 3, the more massive BHs have a higher contribution of their BH growth from gas accretion, with ≳ 10 8 M ⊙ BHs having ∼ 90% of the BH mass acquired via gas accretion.
that grow out of the initial 2.3 × 10 3 M⊙ seeds, let us first look at the relative contributions of BH growth from mergers vs gas accretion.We show this in Figure 8 by taking one each amongst the secondary BRAHMA-9-E4 (top panel) and BRAHMA-18-E5 (bottom panel) boxes that are evolved all the way to z ∼ 0. Specifically, we plot the ratio of the BH mass assembled via mergers vs. accretion from z = 11 to z ∼ 0. At z ≳ 7, the contribution from gas accretion is very small.This is also seen in previous works (Taylor & Kobayashi 2014;Habouzit et al. 2017).For the lowest mass BHs, the small accretion rates may be contributed by the M 2 bh scaling of the Bondi-Hoyle accretion rates.However, at z ≥ 7, the accretion driven BH growth is small even for the most massive BHs formed in our simulations.This is likely due to the influence of stellar feedback that can efficiently expel gas from the centers of low mass halos at high-z, as shown in Habouzit et al. (2017).In fact, the fraction of BH mass accumulated via gas accretion becomes ≳ 50% only after z ∼ 3 for ≳ 10 6 M⊙ BHs.At z ≲ 3, gas accretion drives up growth rates of ≳ 10 6 M⊙ BHs.
Keeping in mind the merger dominated BH growth at z ≥ 7, we now focus on the number density evolution of higher mass BHs in Figure 7.For the ≥ 1.8 × 10 4 M⊙ BHs in the BRAHMA-18-E4 box, the SM5_FOF3000, SM5_FOF10000 and SM150_FOF3000 seed models are reasonably well probed by both BRAHMA-9-D3 and BRAHMA-18-E4 (2nd panels).The most lenient SM5_FOF1000 seed model could only be shown for BRAHMA-9-D3, since most of 1.8 × 10 4 M⊙ BHs assemble in galaxies that are too small to be resolved in BRAHMA-18-E4 (revisit Figure 3).The ≥ 1.8 × 10 4 M⊙ BHs start to assemble between z ∼ 20 − 25, and continue to increase even after the 2.3 × 10 3 M⊙ DGBs are suppressed at z ∼ 10.By z = 7, the number densities of ≥ 1.8 × 10 4 M⊙ BHs range between ∼ 0.2 − 1 Mpc −3 for the different seed models.
For the ≥ 1.5 × 10 5 M⊙ BHs (3rd panels), we calibrate the lenient seed models for the BRAHMA-36-E5 simulations, wherein the assembly of these BHs starts at z ∼ 20.For the strict seed models, we could not calibrate the BRAHMA-36-E5 simulations due to extremely poor statistics of ≥ 1.5×10 5 M⊙ BHs in BRAHMA-9-D3; however, we do obtain some reasonable statistics within the BRAHMA-18-E4 boxes, wherein they start assembling at z ∼ 13.The number densities of ≥ 1.5×10 5 M⊙ BHs range from ∼ 0.01 − 0.2 Mpc −3 at z ∼ 7. Finally, the number density evolution of ≥ 1.2 × 10 6 M⊙ BHs can be reasonably well probed only for the lenient seed models in the BRAHMA-36-E5 boxes, wherein they start assembling at z ∼ 15−16.For these two models, the predicted z ∼ 7 number densities of ≥ 1.2 × 10 6 M⊙ BHs are ∼ 0.007 − 0.01 Mpc −3 .
The variations in the number densities of ∼ 10 4 − 10 6 M⊙ BHs across different seeding models can be up to factors of ∼ 100, which demonstrates that the memory of the initial ∼ 10 3 M⊙ seeds is retained within the number densities of higher mass ∼ 10 4 − 10 6 M⊙ BHs.This is a natural consequence of the merger dominated BH growth at z ≥ 7, wherein the growth of higher mass BHs strongly depends on the avail-ability of earlier generations of lower mass progenitors, all the way down to the initial seed (DGB) population.In Paper I, we describe the details of how the interplay of halo growth, dense gas formation and metal enrichment determines the impact of seeding criteria on the formation of 2.3 × 10 3 M⊙ DGBs at different redshifts.To briefly summarize, the halo mass seeding threshold Mh is more consequential to seeding at the highest redshifts (z ≳ 15) where the lack of halo growth is the major impediment.On the other hand, the dense & metal poor gas mass threshold Msfmp is consequential at all redshifts because the lack of dense gas impedes seeding at the highest redshifts, and the lack of metal poor halos impedes seeding at lower redshifts (z ≲ 10).
We now compare our number density evolution predictions to existing observational constraints.Note that our number densities include both active and inactive BHs, for which we have observational constraints only at z ∼ 0. Therefore, we use the lower-resolution secondary boxes (revisit Table 1) to push our predictions all the way to z ∼ 0. In Figure 9 which shows the number density evolution of ≥ 1.2 × 10 6 M⊙ to z ∼ 0, we applied the lenient seed models on BRAHMA-18-E5 boxes resolving 1.5 × 10 5 M⊙ ESDs.Recall that for the more restrictive strict seed models, we could not calibrate them for 1.5 × 10 5 M⊙ ESDs.Therefore, we run them on BRAHMA-9-E4 boxes that resolve 1.8 × 10 4 M⊙ ESDs.The seed model variations in the number density predictions become smaller with time.In fact, for the ≥ 1.2×10 6 M⊙ BHs, all the seed models predict similar number densities at z ∼ 0, which do approach the observational constraints (black markers in Figure 9) derived by integrating the BH mass function (BHMF) measurements from Merloni & Heinz (2008), Cao (2010) and Shen et al. (2020) (hereafter M08,C10 and S20).In addition, we compare our results with those of TNG100 and TNG300 boxes, which seed 1.2 × 10 6 M⊙ BHs in halos above 7.3 × 10 10 M⊙.Our results are consistent with the TNG boxes at z ∼ 0. At higher redshifts (z ≳ 3), the TNG number densities are lower than the more lenient seed models, but higher than the strict seed models.

BH mass functions
Figure 10 shows the z ≥ 7 BH mass functions (BHMFs) for the different seed models.While the BRAHMA-9-D3 boxes can probe the BHMFs from ∼ 10 3 − 10 5 M⊙, the larger BRAHMA-18-E4 and BRAHMA-36-E5 boxes together extend the predictions up to ∼ 10 7 M⊙.Notably, at the most massive end probed by BRAHMA-9-D3 i.e. (∼ 10 5 M⊙), the BRAHMA-18-E4 and BRAHMA-36-E5 boxes do predict somewhat higher BH abundances compared to BRAHMA-9-D3.In Appendix A, we show that this difference is primarily due to cosmic variance.Importantly, we see that the differences in BHMF predictions between the different seed models are larger than the cosmic variance between different boxes for the same seed model.For more restrictive seed models, predicted BHMFs are steeper, leading to the largest variations (factors ∼ 10) predicted by different seed models at the highest BH masses of ∼ 10 6 − 10 7 M⊙ at z = 7.This is a natural consequence of the merger dominated growth of BHs at these epochs, which makes it harder to grow BHs when fewer seeds are formed.
However, it is possible that the above trends in the BHMFs may be different in the regime of the observed luminous ∼ 10 9 − 10 10 M⊙ quasars at z ≳ 7, with masses that are substantially higher than the BHs in our volumes.These quasars are believed to reside in regions that are much more overdense than the ones probed by our simulations.As seen in Bhowmick et al. (2022b) using "constrained" simulations, gas accretion starts to dominate over mergers much earlier (z ∼ 9) in these extreme regions.In the future, we aim to run much larger boxes to understand how the BHMF properties in this "hypermassive" regime compare to the predictions of our current work.
We can put our high-z ∼ 10 3 − 10 7 M⊙ BHMF predictions in the context of local measurements from ∼ 10 6 − 10 9 M⊙ BHs (shown as black data points in Figure 10).Notably, the abundances of ∼ 10 6 M⊙ BHs predicted by our lenient seed models at z ≥ 7, are similar to the z ∼ 0 measurements (this can also be seen in Figure 9).However, our high-z BHMF predictions have substantially steeper slopes compared to the local BHMFs.Therefore, the abundances of ≳ 10 7 M⊙ BHs at z ≳ 7 are ≳ 10 times lower than the observed local BH populations.In Figure 11, we take our secondary BRAHMA-18-E5 simulation boxes and trace the evolution of the BHMF to z = 0. Note that our simulations are still relatively small, such that they only probe BH masses below the 'knee' of the presumably Schechter-like shape of the overall BHMF.Therefore, we are essentially tracing the redshift evolution of the slope of the low mass end of the BHMFs.As we approach the local universe, the predicted slopes flatten significantly and become similar to that of the observed BHMFs.This is due to the continued growth of massive BHs driven primarily by gas accretion, particularly from z ∼ 3 to z ∼ 0 (revisit Figure 8).
We now compare in more detail the observational measurements of the z ∼ 0 BHMF to that of our predictions.While the measurements span ∼ 10 6 − 10 9 M⊙, our BRAHMA-18-E5 and BRAHMA-9-E4 boxes do not have enough volume to produce BH masses above ∼ 10 8 M⊙ and ∼ 10 7 M⊙ respectively.For the ∼ 10 6 − 10 7 M⊙ regime that is sufficiently well probed by the BRAHMA-18-E5 boxes, there are some differences amongst the various observed BHMF measurements (as also reflected in the overall number density measurements shown in Figure 9).Our BRAHMA-18-E5 boxes with the lenient models produce the best agreement with the most recent S20 results, as also seen in the number density evolution.The BRAHMA-9-E4 boxes with the strict models also produce z ∼ 0 BHMFs that approach the local observations at ∼ 10 6 − 10 7 M⊙, but here the uncertainties due to Poisson variance and cosmic variance are higher than the BRAHMA-18-E5 boxes.Additionally, we see that all the seed models converge to produce similar BHMFs by z ∼ 0, particularly for the most massive ∼ 10 6 − 10 7 M⊙ BHs probed by our secondary boxes.We also plot the TNG100 and TNG300 predictions and show that they are similar to our BRAHMA boxes at ∼ 10 6 − 10 7 M⊙ (recall that the TNG seed mass is 1.2 × 10 6 M⊙).This further reinforces that our seed models do not really have any consequences on the BHMFs of ≳ 10 6 M⊙ BHs at z ∼ 0. This is in stark contrast to the trends at z ≥ 7 wherein the simulated BHMFs have strong differences particularly at ≳ 10 6 M⊙.Again, this is because at z ≥ 7, all ∼ 10 3 − 10 7 M⊙ BHs are primarily undergoing merger driven BH growth that tends to retain the memory of the seeding origins (revisit Figure 8).On the other hand, by z ∼ 0, gas accretion erases the memory of the seeding origins of ≳ 10 6 M⊙ BHs.The blue circles show the very recent BHMF constraints at z ∼ 5, inferred from the recent faint AGN observations at z ∼ 5 (referred to as "little red dots", Matthee et al. 2023).As we approach z ∼ 0, we note the following: 1) the differences between the various seed model predictions become smaller, and 2) the BHMF slopes become flatter due to formation of more massive BHs contributed by gas accretion at z ≲ 3. The BHMFs approach the local constraints at z ∼ 0.

AGN luminosity functions
In Figure 12, we show predictions for the z ≥ 7 AGN bolometric luminosity functions (LFs) for the different seed models.The BRAHMA-9-D3 simulations are able to probe only faint AGN up to ∼ 10 42 erg s −1 ; this is far below the detection limits of not only JWST, but also any potential future missions including NASA APEX.The BRAHMA-18-E4 and BRAHMA-36-E5 boxes help us extend our LF predictions up to ∼ 10 44 erg s −1 at z = 7 (left most panel), which could be within the reach of an APEX X-ray probe, and close to the JWST detection limit at z = 7.We do not show the BRAHMA-18-E4 and BRAHMA-36-E5 predictions at luminosities ≲ 10 41 erg s −1 , since they are artificially suppressed due to the resolution limits of the seed (ESD) masses.
For a given seed model, the BRAHMA-18-E4 and BRAHMA-36-E5 predictions are consistent with the BRAHMA-9-D3 results at luminosities of ∼ 10 41 − 10 42 erg s −1 where they overlap with each other.Additionally, when we combine the predictions for all three boxes, they align with each other to approximately form a single continuous LF that ranges from 10 38 − 10 44 erg s −1 .This is yet another signature of the success of our stochastic seed models developed in Paper I.
The AGN fractions in Figure 13 show that more BHs are active in the strict models compared to the lenient ones.In other words, even if we increase the abundances of massive (∼ 10 6 − 10 7 M⊙) BHs by making the seed model less restrictive, only a certain fraction of them end up in environments that support sufficient AGN activity.This limits the imprint of seed models on the AGN LFs despite the significant seed model variations in the BHMFs that include both active and inactive BHs.At the highest luminosities potentially accessible with an APEX X-ray probe, the seed model variations in the LFs are generally small (factors of ∼ 2 − 3).As a result, it may be challenging to constrain seed models using measurements of z ≥ 7 AGN LFs alone.
In the rightmost panel of Figure 12, we replot the z = 7 AGN LF predictions against observational measurements of pre-JWST z ∼ 6 quasars performed by S20.z ∼ 6 is the highest redshift at which bolometric LF measurements have been carried out to date.Low luminosity high-z AGNs detected by JWST have not yielded precise constraints of the Figure 12.Predictions of the z ≥ 7 bolometric luminosity functions (LFs) for the different seed models using our primary boxes.For a given seed model, we show the BRAHMA-9-D3, BRAHMA-18-E4 and BRAHMA-36-E5 as lines of different thickness.In the 4th panel, we replot the z = 7 predictions along with the observational constraints at z ∼ 6.The two thin solid lines are predictions from TNG100 and TNG300 at z = 6.The errorbars show 1σ Poisson errors.The dotted, dotdashed, dashed vertical lines show the detection limits of Athena, JWST, and a potential NASA APEX X-ray probe (with assumed limit similar to the proposed AXIS telescope: Reynolds et al. 2023) respectively, derived using bolometric corrections of Vasudevan & Fabian (2007).The seed model variations at the most luminous end probed by our simulations (i.e.∼ 10 42 − 10 44 erg s −1 ), are markedly smaller (by factors of ∼ 2 − 3) than the variations in the massive end of the BHMFs.Given these variations and the anticipated hurdles in constraining the faint end AGN luminosity functions such as detection limits and obscuration, constraining seed models using the AGN luminosity functions alone is likely going to be challenging.
bolometric LFs yet, likely because the sample size is still very small and the selection function (including the target selection completeness) with instruments like NIRSpec/MSA is not known (as stated in Harikane et al. 2023).Additionally, the bolometric luminosities for these AGNs are also uncertain as these AGNs are currenly only observed in rest frame UV.Our boxes are not large enough to probe the S20 quasar luminosities and number densities (these objects were probed in the constrained simulations of rare high-z overdense regions in Bhowmick et al. 2022b).Nevertheless, it is encouraging to see that our predictions do broadly align with straightforward extrapolations of the S20 constraints.We also see that the TNG100 and TNG300 predictions are similar to BRAHMA-36-E5 at ∼ 10 43 erg s −1 , and are also in reasonable agreement with the S20 constraints at ∼ 10 45 erg s −1 .
At lower redshifts (z ≲ 5), the observations are able to probe luminosities as low as ∼ 10 43 erg s −1 , that partially overlap with the simulated LFs for the secondary BRAHMA-18-E5 and BRAHMA-9-E4 boxes shown in Figure 14.However, the measured z ∼ 0 − 5 LFs suggest a flattening of the slope at ∼ 10 43 − 10 45 erg s −1 .Due to this, both the BRAHMA as well as the TNG boxes predict more numerous ∼ 10 43 −10 45 erg s −1 AGN than the observed LFs at z ∼ 1−5 even though there is reasonable agreement between observations and TNG at the bright end (≳ 10 45 erg s −1 ).In fact, several other cosmological simulations also produce a similar tension with the observed LFs at ∼ 10 43 − 10 45 erg s −1 (see Figure 5 of Habouzit et al. 2022).In addition, the fact that the simulated LFs are not very sensitive to differences in BH seed models at these luminosities and redshifts, suggests that the uncertain seeding origin is less likely to explain this discrepancy with observations (see also Figure 4 of Bhowmick et al. 2021).Possible reasons for these discrepancies have been discussed in several works (Weinberger et al. 2018;Bhowmick et al. 2021;Habouzit et al. 2022); these include 1) uncertainties in the modeling of AGN obscuration in observations, and 2) uncertainties in the modeling of BH accretion, feedback and dynamics in simulations.In future work, we will tackle this issue by exploring alternative models for BH accretion, feedback and dynamics, as well as by predicting the amount of obscuration as a function of AGN luminosity.

Stellar mass vs BH mass relations
In Figure 15, we show the stellar mass vs. black hole mass (M bh − M * ) relations at z = 7, 10 & 11 for the simulated galaxies (Subfind-subhalos) for the various seed models shown in different columns.The three boxes together probe stellar masses ranging from ∼ 10 4 − 10 8 M⊙.The BRAHMA-18-E4 and BRAHMA-36-E5 boxes do a reasonably good job in modeling the BH growth above 8 M ESD seed , which is demonstrated by the fact that the scaling relations predicted for the three boxes match reasonably well with each other for a given seed model2 .That being said, there is slightly higher BH growth in BRAHMA-18-E4 and BRAHMA-36-E5 compared to BRAHMA-9-D3, which is more readily seen in the BHMFs discussed in the previous subsection.As we demonstrate in Appendix A, this difference is primarily due to cosmic variance.Comparing the M bh − M * relations for the different seed models, we can clearly see that the strict seed models produce smaller BH masses at fixed stellar mass, compared to the lenient seed models.
For all the seed models, the z ≳ 7 BH masses are significantly above expectations (at fixed stellar mass) from simple extrapolations of the local observational constraints (shown by the orange color in Figure 15).However, the host galaxies do tend to grow faster than the BHs from z = 11 − 7, leading to a rightward shift in the predicted M bh − M * relations with time.In Figure 16, we use the BRAHMA-18-E5 and BRAHMA-9-E4 boxes to continue tracing the preferential build up of stellar mass to z ∼ 0. All of our predictions approach the local measurements by z ∼ 0, thereby serving as further validation for our seed models.
Recent JWST discoveries of high-redshift, moderateluminosity AGN are pushing our current frontiers for M bh − M * relation measurements at early times.We compare them to our predictions in Figure 15.The pre-JWST highest red- The fractions are exactly zero for the > 10 42 erg s −1 AGN in BRAHMA-9-D3 boxes (lower left panel), because there are no such AGNs formed due to limited volume.Generally, we see lower AGN fractions for the more lenient seed models.Therefore, even if we produce more BHs with more lenient seed models, only a certain fraction of them end up in environments that support enough gas accretion to become AGN.This explains why the seed model variations in the LFs are smaller than that in the BHMFs.
shift GN-z11 galaxy has now acquired a new title for hosting the highest redshift AGN (Maiolino et al. 2023a) at z ∼ 10.6.With an estimated BH mass of ∼ 10 6 M⊙ in a ∼ 10 9 M⊙ galaxy, its stellar to BH mass ratio is comparable to the local scaling relations (black square in Figure 15, top row).Therefore, its BH mass is significantly below simple extrapolations of our simulated z = 11 M bh − M * relations for all the seed models.This suggests that the seeding mechanism for the GN-z11 BH may not be encompassed by our explored set of seed parameters.On the other hand, the z = 10.1 AGN spectroscopically confirmed by Goulding et al. (2023), has a predicted mass of ∼ 10 7 − 10 8 M⊙ in a ∼ 1.4 × 10 8 M⊙ galaxy (black circle in Figure 15, middle row).This AGN has a significantly overmassive BH compared to the local relations, which is consistent with extrapolations of our predicted M bh − M * relations, particularly for our lenient seed models.
There are also several recent measurements of the M bh −M * relations for AGNs between z ∼ 4 − 7 (Maiolino et al. 2023b;Harikane et al. 2023;Pacucci et al. 2023) (grey diamonds and circles in Figure 15).We can compare our z = 7 predictions with three objects amongst these samples that have spectroscopically confirmed redshifts between z = 6 − 7 (black diamonds in Figure 15, bottom row).Two of these objects have BH masses close to the upper end of the scatter in the local scaling relations of Reines & Volonteri (2015), and one of them lies well within the scaling relations.All three objects lie below the M bh − M * relation predictions from our lenient seed models.However, these objects do align well with simple extrapolations of predictions from our strict seed models.
Notably, the vast majority of the z ∼ 4 − 6 AGN analyzed by Pacucci et al. (2023) and Maiolino et al. (2023b) are well above the local scaling relations.While the error bars in these measurements (not shown for clarity) are still very large for both BH masses (up to ∼ 1 dex) and the host stellar masses (typically ∼ 0.5−1 dex but can be up to ∼ 2 dex), the feasibility of assembling such overmassive AGN at high-z is still an interesting theoretical question with significant potential implications for BH seeding.In the top panels of Figure 16, we can see that these "overmassive" AGNs align well with BRAHMA-18-E5 predictions at z = 5 for the lenient seed models.The strict seed model predictions from BRAHMA-9-E4 tend to be more consistent with the BHs lying at the lower end of the observed scatter.By z = 3 (middle panels of Figure 16), the simulated M bh − M * relations are below the overmassive JWST AGN, as they continue their approach to become consistent with the local measuremnets at z ∼ 0. Overall, while overmassive black holes (BHs) measured by Goulding et al. (2023) and Pacucci et al. (2023) may indicate heavy DCBH seeds (as proposed by Natarajan et al. 2023;Pacucci et al. 2023), our predictions show that is it is not impossible to assemble them from lighter seeds if they can form and merge efficiently enough.

BH binary formation and merger rates
We now look at the formation rates of BH binaries at z ≥ 7 for the different seed models.These BH binaries are possible precursors to LISA mergers; however, the limited resolution of our simulations do not allow us to probe the hardening of these binaries all the way to the actual merger.As mentioned earlier, in our simulations, these binaries are promptly merged when one of the pairs comes within the neighbor search radius (R Hsml ) of its companion.R Hsml forms a sphere around the BH enclosing 256, 128 and 64 neighboring gas cells in BRAHMA-9-D3, BRAHMA-18-E4 and BRAHMA-36-E5 respectively.In Figure 17, we plot the distributions of R Hsml (the maximum distance between the primary and secondary BHs) for all the binaries at the time of their "merger" in the simulations.These distributions are essentially a measure of the minimum separations up to which we are able to probe these binaries.This minimum separation naturally depends on the overall resolution as seen in the third panel (compare lines of different thicknesses with the same color) where ≥ 1.5 × 10 5 M⊙ binaries merge at smaller separations in BRAHMA-9-D3 (∼ 1 − 5 kpc) compared to BRAHMA-18-E4 (∼ 2 − 10 kpc) and BRAHMA-36-E5 (∼ 4 − 15 kpc).At the same time, due to the refinement and derefinement of the gas cells, the minimum separation is also influenced by the local spatial resolution of the gas around the binary.In particular, the binaries residing within higher gas densities can be traced to smaller separations because the individual gas cells are smaller in size.This is essentially responsible for the differences in R Hsml distributions between the different seed models, as well as the different BH mass thresholds.As we make our seed model more restrictive, the BH binaries have smaller separations at the time of their merger (see lines of differ-Figure 14.Here we use the secondary boxes to show the luminosity functions predicted by our seed models to z ∼ 1.The dark and light green dashed lines show predictions for the lenient seed models using BRAHMA-18-E5.The maroon and pink dashed lines show the strict seed model predictions using BRAHMA-9-E4.The two thin solid lines are predictions from TNG100 and TNG300.The grey points show different observational measurements.The LF predictions for different seed models become very similar at z ≲ 3. Our simulations have higher abundances of ∼ 10 43 − 10 44 erg s −1 AGN compared to observations, as is the case with Illustris-TNG and several other simulations. ent colors in a given panel of Figure 17).This is because for more restrictive seed models, DGBs and their descendants assemble in increasingly massive halos that typically also have higher gas densities (higher spatial resolution) at their centers.For the same reason, we also see that higher mass BH binaries living in higher mass halos can be traced to smaller separations (compare lines of the same color and linestyle in different panels).In BRAHMA-9-D3, the ≥ 2.3 × 10 3 M⊙ (1st panel) binaries merge at ∼ 3 − 15 kpc separations, whereas the ≥ 1.8 × 10 4 M⊙ (2nd panel) and ≥ 1.5 × 10 5 M⊙ (third panel) binaries merge at ∼ 1 − 5 kpc separations.
Overall, the separations at which our binaries merge range from ∼ 1 − 15 kpc depending on the overall simulation resolution as well as the local gas environment.However, even at scales larger than these separations, we are not able to accurately probe the orbital dynamics of these binaries due to the BH repositioning scheme.In other words, because of the repositioning scheme, we are prone to underestimating the time scales for the assembly of these ∼ 1 − 15 kpc scale binaries (let alone their continued evolution to merger).While not desirable, it is a necessary measure to ensure the kinematic stability of our BHs.This is particularly crucial for our smallest BHs close to the DGB and ESD seed mass, since they are highly vulnerable to "artificial" numerical kicks from the ∼ 10 times heavier background DM particles.Nevertheless, in Appendix B, we test the impact of removing the repositioning scheme and find that it reduces the BH binary formation rates by a factor of ∼ 3 at z ≳ 7.While the impact is not negligible, it is still relatively small.This also implies that if we were to account for the missing dynamical friction using subgrid prescriptions such as those proposed by Tremmel et al. (2017) or Ma et al. (2023), the actual reduction in the BH binary formation rates at ∼ 1−15 kpc separations should be within a factor of three.In the future we plan to explore these sub-grid dynamical friction models, which will allow us to trace the merging BH binaries up to ∼ 0.1 kpc (related note: dynamical friction is expected to be a dominant binary hardening mechanism only at ≳ 0.1 kpc scales; Kelley et al. 2017;Sayeb et al. 2021).
Keeping in mind that we are promptly merging the ∼ 1 − 15 kpc scale BH binaries in our simulations, we compute the merger rates as: where d 2 N is the number of BH binaries that would be observed (if they promptly merged as done in the simulations) between z and z + dz, over a time-interval dt obs in the observer frame. 1 Vc dN dz is the comoving number density of merger events recorded in the simulation per unit redshift interval.4πr 2 c drc dz is the comoving volume per unit redshift interval at redshift z (rc is the comoving distance).Lastly, the extra 1 1+z factor is to convert the time interval dt in the rest frame of a merger event between z to z + dz, to the observer frame.For simplicity of language, we shall simply refer to the rates of these promptly merging BH binaries as "merger rates".
The merger rates are plotted in Figure 18 for different mass thresholds for each binary component.Let us first focus on the z ≥ 7 predictions made by our primary boxes (solid lines).We start with the leftmost panel, wherein both the primary and secondary BH masses (M1 and M2 respectively) are ≥ 2.3 × 10 3 M⊙.For the most lenient SM5_FOF1000 seed model, the earliest mergers occur at z ≳ 25; the merger rates peak at ∼ 2000 events per year at z ∼ 12.In contrast, for the more restrictive seed models like SM150_FOF3000, the first merger events occur at z ∼ 20 − 23, with maximum rates of ∼ 200 events per year at z ∼ 10.Generally, depending on the seed model, these rates can peak anywhere between z ∼ 12 − 10, and subsequently drop due to the slow down and eventual suppression of DGB formation due to metal enrichment.
For more massive BH binaries, the larger boxes that use the stochastic seed model make more statistically robust predictions.For the merger rates of ≥ 1.8 × 10 4 M⊙ BH binaries, we use the BRAHMA-18-E4 predictions.At z ∼ 15, we predict ∼ 8 − 10 mergers per year for the lenient seed models, and ≲ 1 mergers per year for the strict seed models.At z ∼ 7 − 8, we predict ∼ 6 − 7 mergers per year for the strict seed models and ∼ 20 − 60 mergers per year for the lenient seed models.
We plot the merger rates of ≥ 1.5 × 10 5 M⊙ BHs in the third panel of Figure 18.For the predictions at z ≥ 7, we make use of the BRAHMA-36-E5 boxes.For the lenient seed models, we predict ∼ 2 − 8 mergers per year.Finally, for the   Figure 18.The solid lines show predictions for the simulated z ≥ 7 merger rates of the BH binaries that were promptly merged at ∼ 1 − 15 kpc separations.These are upper limits to the true merger rates to be measured by LISA.The 1st panel from the left shows primary and secondary BHs above 2.3 × 10 3 M ⊙ , which are only probed by the BRAHMA-9-D3 boxes.The 2nd panel shows primary and secondary BHs above 1.8 × 10 4 M ⊙ , probed most effectively by the BRAHMA-18-E4 boxes.Note that the BRAHMA-18-E4 predictions are not present for the SM5_FOF1000 model since the calibration was not possible.The rightmost two panels show merging BHs above 1.5 × 10 5 M ⊙ and 1.2 × 10 6 M ⊙ , probed most effectively by the BRAHMA-36-E5 boxes.Here, we could only calibrate the SM5_FOF1000 and SM5_FOF3000 seed models.The dashed lines show the secondary BRAHMA-18-E5 and BRAHMA-9-E4 boxes, which extend the merger rate predictions to z ∼ 0. The thin solid lines with squares show predictions from Illustris and IllustrisTNG simulations.The black markers correspond to SAMs and empirical model predictions.At z ≳ 7, our seed models predict between ∼ 200 − 2000 mergers per year amongst ≥ 2.3 × 10 3 M ⊙ BHs depending on the seeding criteria.Likewise, we predict ∼ 6 − 60 mergers per year amongst ≥ 1.8 × 10 4 M ⊙ BHs, and up to ∼ 10 mergers per year amongst ≥ 1.5 × 10 5 M ⊙ BHs.Between the different seed models, the merger rate variations at a given redshift can be up to factors of ∼ 100.
≥ 1.2×10 6 M⊙ BHs, there are only a handful of z ≥ 7 mergers in the lenient seed models, from which we can roughly infer ∼ 0.02 − 0.1 mergers per year.
Next, we use the less expensive BRAHMA-18-E5 and BRAHMA-9-E4 boxes to look at what our seed models predict for the merger rates at z < 7. The first thing we note is that while the ≥ 2.3 × 10 3 M⊙ merger rates peak anywhere between z ∼ 10 to 12 depending on the seed model, the merger rates of the more massive BHs continue to increase below z = 7, owing to successive mergers and accretion episodes.The ≥ 1.5 × 10 5 M⊙ and ≥ 1.2 × 10 6 M⊙ BH merger rates peak at z ∼ 5 and z ∼ 2 respectively.The lenient seed models predict ∼ 0.7 mergers per year amongst ≥ 1.2 × 10 6 M⊙ BHs.

Comparison to the Illustris family of simulations
We can compare our merger rate predictions to two predecessors of BRAHMA; i.e the Illustris and IllustrisTNG simulations.Both these simulations were run using AREPO with the same BH repositioning scheme to stabilize the small scale BH dynamics.Illustris placed seeds of mass 1.5 × 10 5 M⊙ within 7.1 × 10 10 M⊙ halos.At z ≥ 7, for our lenient SM5_FOF1000 and SM5_FOF3000 seed models, BRAHMA predicts ≳ 100 times more mergers amongst ≥ 1.5 × 10 5 M⊙ BHs compared to Illustris.This is due to the stricter halo mass threshold adopted in Illustris, compared to the typical halos in which 1.5 × 10 5 M⊙ BHs assemble (∼ 10 9 − 10 10 M⊙) within BRAHMA.Our merger rates will therefore also be higher than other simulations that adopt similar halo mass seeding thresholds, such as MASSIVEBLACK II (Khandai et al. 2015), EAGLE (Salcido et al. 2016) and ASTRID (Chen et al. 2022;De-Graf et al. 2023).There is only a handful of ≥ 1.2 × 10 6 M⊙ BH mergers at z ≥ 7 in BRAHMA (for SM5_FOF1000) due to limited volume; but these events still imply higher merger rates in BRAHMA compared to Illustris at z ≥ 7 for our lenient seed models.
At z ≲ 1, the differences in the ≥ 1.5 × 10 5 M⊙ merger rates between BRAHMA and Illustris become much smaller.This is driven by the slow down of new 2.3 × 10 3 M⊙ DGB formation (and subsequent assembly of new 1.5 × 10 5 M⊙ BHs) due to metal enrichment in BRAHMA.Similarly, the ≥ 1.2 × 10 6 M⊙ BHs also merge at similar rates in BRAHMA and Illustris at z ≲ 1.
Because Illustris and IllustrisTNG have very similar merger rates for ≥ 1.2 × 10 6 M⊙ BHs, the same conclusions can be drawn for the comparison between BRAHMA vs IllustrisTNG.Notably, the ≥ 1.2 × 10 6 M⊙ BH merger rates are slightly higher in TNG300 than in TNG100, despite the two simulations producing very similar BH number densities.This is likely due to the larger volume of TNG300 that captures a higher number of extreme overdense regions where galaxies and BHs are expected to merge more frequently.

Comparison to SAMs
We can also compare our results to SAMs even though their seeding prescriptions, treatment of mergers, as well as the sample selection used are generally very different from our work.We compare our predictions against the results of Barausse (2012, B12), Ricarte & Natarajan (2018, RN18), and Dayal et al. (2019, D19).Based on their sample selection, we placed their predictions on different panels of Figure 18.Furthermore, to facilitate the most fair comparison, we only look at those SAM results which are made under assumptions similar to those of our prompt mergers.For example, B12 merges their BHs immediately after the satellite and central galaxies merge.For RN18, we make comparisons with their "optimistic scenario" wherein the merging probability of BHs (after their host galaxies merge) is unity.Similarly, for D19, we compare the BRAHMA merger rates with the "instantaneous (ins)" merging scenario where BHs and galaxies are assumed to merge as soon as the host halos merge.
In general, the overall merger rates in BRAHMA are ∼ 10−100 times higher than the SAM predictions for "light seeds" (1st and 2nd panels in Figure 18).The major contributor to this is the fact that these SAMs adopt significantly more stringent seeding criteria.For example, B12 allows seed formation only after z = 20, and RN18 also places both "light" and "heavy" seeds between z = 15−20.On the other hand, for most of our seed models the 2.3 × 10 3 M⊙ DGBs start forming at z ≳ 20.As for D19, their assumed halo mass threshold of ≥ 10 8 M⊙ is much higher (due to the minimum halo mass they resolve) than the ones adopted in BRAHMA (∼ 10 6 − 10 7 M⊙).
If we only consider the massive mergers amongst ≳ 10 5 M⊙ BHs (3rd panel in Figure 18), our low mass seed models in BRAHMA predict higher rates at high redshifts (z ≳ 10) than the B12 and R18 SAMs which model "heavy" ∼ 10 5 M⊙ seeds.This is because the seeding criteria in the B12 and R18 SAMs initialize ∼ 10 5 M⊙ seeds faster than the rates at which similar mass BHs are assembled from ∼ 10 3 M⊙ seeds in BRAHMA.Notably, at lower redshifts (z ≲ 5), the BRAHMA predictions do become similar to the B12 and R18 SAMs.
Lastly, B12 predicts much lower merger rates for ≳ 10 6 M⊙ BHs in their heavy seed model (4th panel in Figure 18), which they say is due to ejections via gravitational-wave recoil.As a result, at all redshifts, their predictions are significantly lower than BRAHMA, which does not include recoiling BHs.
The above comparisons further highlight the profound implications of seed models for LISA.They also underscore the importance of exploring our seed models coupled with other galaxy formation models with alternate treatments of star formation, metal enrichment and stellar feedback; we shall pursue this in the future.In addition, the merger rates will be influenced by the detailed BH dynamics on small scales that our simulations do not capture.Due to this, it is best to interpret the simulated merger rates as upper limits to the actual GW event rates detectable by LISA.In a follow-up paper, we shall use post processing models (Kelley et al. 2017;Sayeb et al. 2021) to not just account for hardening due to dynamical friction up to ∼ 0.1 kpc separations, but also the continued hardening of the BH binaries at ≲ 0.1 kpc separations due to processes that are presumed to occur at scales below the simulation resolution, such as stellar loss cone scattering, viscous dissipation due to circumbinary gas disk, and GW emission and recoil.

SUMMARY AND DISCUSSION
Modeling of low mass seeds in cosmological simulations has been a longstanding challenge due to dynamic range limitations.We have developed a new suite of cosmological hydrodynamic simulations named BRAHMA, which adopts a novel flexible seeding framework that enables the representation of seeds that are ∼ 10 − 100 times below the simulation mass resolution.With the exception of the seed models, our simulations adopt the underlying galaxy formation model from Illustris-TNG.Our primary set of simulations is comprised of three simulation volumes that use two distinct seeding prescriptions depending on whether we can directly resolve our target seed mass.In BRAHMA-9-D3, the smallest of these volumes (9 Mpc box-length) that explicitly resolves our target seed mass of 2.3 × 10 3 M⊙, we use a gas-based seeding model.Here, seeds are placed inside halos exceeding critical thresholds of dense & metal-poor gas mass and total halo mass ( Msfmp and Mh in the units of the seed mass).In the larger-volume lower-resolution boxes, we do not resolve our target seed mass.Instead, we seed the higher mass descendants of the target seed masses using our new stochastic seed model developed in Paper I. The BRAHMA-18-E4 (∼ 18 Mpc box-length) and BRAHMA-36-E5 (∼ 36 Mpc box-length) simulations resolve 1.8 × 10 4 & 1.5 × 10 5 M⊙ descendants respectively.Calibrated based on where these descendants assemble in BRAHMA-9-D3, the stochastic seed model places them over a broad range of galaxy masses, with a stronger probability of seeding them in rich environments with at least one neighboring halo.Using this simulation suite, we predict the z ≥ 7 IMBH and SMBH populations with masses ranging from ∼ 10 3 − 10 7 M⊙, for multiple sets of gas-based seed parameters Mh and Msfmp .
The following are our key findings: • At z ≥ 7, the BH growth in our simulations is dominated by BH mergers.This is likely due in part to stellar feedback that makes gas unavailable to fuel gas accretion in galactic nuclei.In addition, the M 2 bh scaling within the Bondi accretion model also contributes to slower BH growth for low mass seeds.This merger-driven growth at high-z tends to retain the imprint of seeding, as the growth of higher mass BHs explicitly relies on the availability of enough seeds to fuel the merger driven BH growth.
• As we make the seed models more restrictive, the z ≥ 7 BHMFs become steeper, leading to significant suppression over the entire BH mass range probed by our simulations (10 3 − 10 5 M⊙).The largest variations (factors ∼ 10) due to different seeding choices are seen for the most massive ∼ 10 6 − 10 7 M⊙ BHs.
• The variations in the AGN LFs are smaller than those of the BHMFs.This is because regardless of the total number of BHs assembled in a given seeding model, there is only a limited set of environments that have enough gas to fuel these BHs to become AGN.At the luminosities potentially reachable with upcoming instruments such as the proposed NASA APEX X-ray probe, the seed model variations are typically within factors of ∼ 2 − 3.This suggests that it may be difficult to distinguish between these seed models using AGN luminosity functions alone.
• On the M * − M bh plane, our simulations generally predict BH masses ∼ 10 − 100 times higher than the local scaling relations, similar to many of the detected "overmassive" JWST AGN at z = 4−10 (Pacucci et al. 2023;Goulding et al. 2023).While these"overmassive" AGN may very well be indicative of heavy (≳ 10 5 M⊙) seeding origins, our simulations do suggest that it is not impossible to assemble these AGN from relatively lower mass seeds if they can form in abundant numbers and merge with one another efficiently enough.
• Our simulations trace the evolution of BH binaries to separations of ∼ 1 − 15 kpc, after which they are promptly merged.The resulting merger rates can therefore be interpreted as possible upper limits to the actual GW event rates.Depending on the seed model, our simulations produce ∼ 200 − 2000 mergers per year amongst ≳ 10 3 M⊙ BHs, ∼ 6 − 60 mergers per year amongst ≳ 10 4 M⊙ BHs, and up to ∼ 10 mergers per year amongst ≳ 10 5 M⊙ BHs at z ≳ 7.For our most lenient seed models, we also predict a handful of mergers between ≳ 10 6 M⊙ BHs over the course of a few years.The significant seed model variations in the merger rates and their redshift distributions highlight the promise of LISA in producing the strongest constraints for seeding.
• Our overall merger rates are generally ∼ 10 − 100 times higher than existing predictions from several "light seed model" SAMs under similar assumptions of prompt BH mergers following the host galaxy mergers.This difference largely originates from more stringent seeding criteria used in SAMs in terms of the allowed redshift range as well as halo mass thresholds for seed formation.
• Despite the significant seed model variations seen amongst the BH populations at z ≥ 7, they become negligible by z ∼ 0. This is because of the significant accretion driven BH growth at z ≲ 3, which erases the imprint of seeding by z ∼ 0. By z ∼ 0, the BH number densities, the BHMFs, and the M * − M bh relations approach the local observational constraints for ≳ 10 6 M⊙ BHs.
The inability to resolve BH dynamical friction in cosmological simulations is a major challenge in probing the dynamics of BH binaries on small scales.A number of studies that have explored dynamical friction either as subgrid prescriptions in zoom simulations (Pfister et al. 2019;Ma et al. 2021) or via direct of direct (unsoftened) two-body encounters in idealized galaxy simulations (Partmann et al. 2023), are finding that low mass seeds have difficulty in sinking to the centers and coalescing with other BHs particularly within high-z halos with shallow gravitational potentials.For the BH that do merge, GW radiation could kick the remnants out of the halo and prevent future mergers.This certainly poses a challenge to the feasibility of merger dominated BH growth that our simulations exhibit at high-z.At the same time, if the BHs cannot efficiently sink to the gas rich halo centers, it also greatly diminishes the prospect of accretion driven BH growth to support the Eddington or super-Eddington accretion rates that apparently need to be sustained by low mass seeds to assemble the highest-z quasars.Furthermore, the theoretical mechanisms for heavier DCBH seeds are likely too stringent to produce sufficient seeds to explain the abundance of SMBHs in the low-z universe.Given these considerations, we cannot yet rule out the avenue of merger-driven growth of low mass seeds.Several channels have been proposed to overcome the barrier posed by the "sinking problem".These include: 1) High seeding efficiency for low mass seeds, so that a small fraction of them can still end up merging.2) Low mass seeds could be embedded within dense nuclear star clusters (not resolvable in the simulations) that would enhance their effective gravitational mass and expedite the sinking as well as preventing the escape of recoiling BHs. 3) Hardening can be expedited by triple BH interactions.We will explore these scenarios in future simulations with alternative treatments for BH dynamics that will include the missing dynamical friction force on scales above the spatial resolution of the simulations.We also plan to use post-processing analytical models to trace the hardening of these BHs at sub-resolution scales to make predictions for LISA.This will not only include dynamical friction, but also processes that are presumably responsible for hardening of the BH binaries on sub-kiloparsec scales; i.e. drag due to circumbinary disks, stellar scattering, and gravitational waves.Despite the above challenges, given the high rates at which BH binaries assemble in our simulations at ∼ 1 − 15 Mpc separations, even if only a small fraction (≲ 10%) of them end up actually coallescing to produce GWs, LISA would still detect enough events to provide strong constraints for seed formation.
In the future, we will continue to expand the BRAHMA suite to include simulations with more seed models, particularly for heavier DCBH seeds.In addition, we will consider alternate BH accretion and AGN feedback models beyond the current Bondi scheme, and explore the prospect of accretion driven BH growth for low mass seeds at high-z.Finally, we will include simulations with alternative galaxy formation models to explore the impact of subgrid recipes of star formation and metal enrichment on BH seeding and growth.
The BRAHMA simulations introduced in this paper together address a major gap in the existing landscape of cosmological hydrodynamic simulations; i.e. the representation of the low mass seeds in larger simulation volumes.This allows us to make predictions for IMBH and SMBH populations produced by low mass seeding channels for current and upcoming observational facilities such as JWST, Athena, NASA APEX, and most importantly, LISA.Our results highlight the transformative role these observatories will play in unveiling the origins of the diverse populations of observed SMBHs and AGN in our Universe.M ⊙ for all four seed models we considered.The rightmost two panels show results for M ESD seed = 1.5 × 10 5 M ⊙ ; we show this only for the two lenient models for which the stochastic seed criteria could be robustly calibrated; i.e.SM5_FOF1000 and SM5_FOF3000.The top panels show the total merger rate.The middle panels show the mass growth rate due to mergers as a function of redshift, which is defined as the total mass of all merging secondary BHs per unit redshift.The light minor mergers dominate the contribution to the BH mass growth for the majority of the time, and they need to be explicitly included in the stochastic seed models.Notably, the relative dominance of the light minor mergers does tend to recede with time, as quantified in the bottom panels that show the mass growth (∆M light minor ) due to the light minor mergers between successive heavy mergers.Except for the most lenient SM5_FOF1000 seed model with M ESD seed = 1.8 × 10 4 M ⊙ , it is not possible to do a robust power-law fit for ∆M light minor due to the very large uncertainties (green dotted line in the bottom left panel).Therefore we decided to adopt a constant ∆M light minor = 4 M ESD seed (blue horizontal line) for the stochastic seed models, as it is broadly consistent with the ∆M light minor results for all the seed models.ZOOM_REGION_z5-D3 (green), ZOOM_REGION_z5-E4 (brown) and ZOOM_REGION_z5-E5 (grey) correspond to the same zoom region at the resolutions of BRAHMA-9-D3, BRAHMA-18-E4 and BRAHMA-36-E5 respectively.The three different resolutions use the gas-based seed model to seed 2.3 × 10 3 M ⊙ , 1.8 × 10 4 M ⊙ and 1.5 × 10 5 M ⊙ DGBs respectively.The BHMFs predicted by all three zooms are similar at the most massive end.Bottom panels: The thin green solid lines correspond to BHMFs of BRAHMA-9-D3.For the BRAHMA-18-E4 and BRAHMA-36-E5 boxes, we compare the BHMFs predicted by the stochastic seed model with 1.8 × 10 4 M ⊙ and 1.5 × 10 5 M ⊙ ESDs (thicker green lines), against that of the gas-based seed model with 1.8×10 4 M ⊙ and 1.5×10 5 M ⊙ DGBs (brown and grey lines respectively).Very importantly, the stochastic seed model included the minor mergers using the constant value assumption.Amongst the BRAHMA-18-E4 and BRAHMA-36-E5 boxes, the gas based seed model and stochastic seed model predictions are similar, particularly at the massive end of the BHMFs beyond what can be probed by BRAHMA-9-D3.This establishes that the differences amongst the predictions of BRAHMA-9-D3, BRAHMA-18-E4 and BRAHMA-36-E5, are largely due to cosmic variance.Therefore, the constant value assumption does not compromise our BHMF predictions, even for BH masses larger than what can be effectively probed by the BRAHMA-9-D3 boxes.
instructive to assess the impact of our repositioning scheme on the merger rates.
In Figure B1, we remove the BH repositioning scheme and re-run one of our seed models (SM5_FOF1000) on a BRAHMA-4.5-D3box that explicitly resolves the 2.3 × 10 3 M⊙ DGBs.We find that the removal of the repositioning scheme reduces the BH merger rates at z ≳ 10 only by factors of ∼ 3 (solid vs dotted lines in Figure B1).This is due to a delay in the occurrence of these merger events, which also causes an increase in the merger rate at lower redshifts; as a result, the peak in merger rates of ≥ 2.3 × 10 3 M⊙ BHs gets shifted from z ∼ 12 to z ∼ 7 when BH repositioning is turned off.Overall, the impact of BH repositioning is much smaller than the differences (factors ∼ 10 − 100) seen between our simulations vs. SAMs.We emphasize that in this test, we did not substitute the BH repositioning with an alternate model to correct for the missing dynamical friction.Inclusion of such a model could bring the predictions even closer to the simulations that use repositioning.Overall this suggests that our BH repositioning scheme is not the major contributor to the differences between the BRAHMA predictions vs. SAMs.

Figure 1 .
Figure 1.Schematic comparison of the three primary BRAHMA simulations i.e.BRAHMA-9-D3, BRAHMA-18-E4 and BRAHMA-36-E5.For each panel, we show the 2D projected gas density (left) and metallicity (right) profiles.Each of these boxes was run using 1024 3 DM particles.As a result they differ in mass resolutions by factors of 8, and spatial resolutions by factors of 2. The panel dimensions span the full box width and height divided into 1024 × 1024 pixels.The line of sight thickness of the pixels is chosen to be 0.8 % of the total box depth.The yellow crosses correspond to positions of BHs.More massive BHs are depicted by larger crosses.The BRAHMA-9-D3 boxes resolve down to 2.3 × 10 3 M ⊙ BHs seeded using the gas-based seed model.The BRAHMA-18-E4 and BRAHMA-36-E5 boxes resolve down to 1.8 × 10 4 & 1.5 × 10 5 M ⊙ BHs respectively, seeded using the stochastic seed model.

Figure 2 .
Figure 2.An example formation site of 2.3 × 10 3 M ⊙ seeds in the BRAHMA-9-D3 box using the gas-based seed model.These seeds are hereafter referred to as direct gas based seeds or DGBs.Left, middle and right panels show the 2D projected plots of the gas density, star formation rate density and gas metallicity, averaged over a 7.4 kpc slice.The yellow circle marks the DGB formation site comprised of dense & metal poor gas.

Figure 3 .
Figure3.Relation between redshift and the galaxy total masses (M galaxy total ) at which 1.8 × 10 4 M ⊙ and 1.5 × 10 5 M ⊙ BHs respectively assemble from 2.3 × 10 3 M ⊙ DGBs in our BRAHMA-9-D3 suite.Each data point corresponds to a recorded instance of assembly when the BH growth is traced along the galaxy merger tree.The left, middle and right panels show the different gas-based seeding models namely SM5_FOF3000, SM150_FOF3000 and SM5_FOF10000 respectively.For each panel, different colors correspond to different values of M assembly = 1.8 × 10 4 & 1.5 × 10 5 M ⊙ .Solid lines show the mean trend and the shaded regions show ±1σ standard deviations.As motivated in Paper I, we fit the mean trends by a double power-law (dashed lines).To obtain fits, we first assumed appropriate ztrans values for the different seed models via visual inspection.The fits were then performed to obtain the slopes at z < ztrans and z > ztrans using scipy.optimize.curve_fit.The final fitted parameters are shown in Table2.These fits are used in our stochastic seed models that directly seeds the descendants at 1.8 × 10 4 M ⊙ and 1.5 × 10 5 M ⊙ masses within the BRAHMA-18-E4 and BRAHMA-36-E5 boxes respectively.To distinguish them from the DGBs, the seeds formed by our stochastic seed models are referred to as extrapolated seed descandants or ESDs.

Figure 4 .
Figure 4. Similar to Figure2, but here we show an example formation site of ESDs within the BRAHMA-18-E4 box using the stochastic seed model.These ESDs are meant to faithfully represent the higher mass descendants of unresolvable 2.3 × 10 3 M ⊙ DGBs.The quantities are averaged over a slice of thickness 22 kpc.The yellow circle marks the ESD formation site, which is a sufficiently massive galaxy that lives in a rich environment (i.e. with ≥ 2 neighboring halos).

Figure 5 .
Figure 5. Left panels: Solid lines show the mean trends in the evolution of the total masses of galaxies wherein 1.8 × 10 4 M ⊙ (top) and 1.5 × 10 5 M ⊙ (bottom) BHs assemble from 2.3 × 10 3 M ⊙ DGBs within the BRAHMA-9-D3 boxes.The error bars in the solid lines show the standard deviation.The dashed lines show the same for ESDs seeded within the BRAHMA-18-E4 and BRAHMA-36-E5 respectively, using the stochastic seed models calibrated based on the BRAHMA-9-D3 results.The close agreement between the solid and dashed lines indicates the success of our galaxy mass criterion.Right panels: The two point clustering of ≥ 1.8 × 10 4 M ⊙ (top) and ≥ 1.5 × 10 5 M ⊙ (bottom) BHs for the BRAHMA-9-D3 (solid lines), BRAHMA-18-E4 and BRAHMA-36-E5 (dashed lines) simulations.Here, the close agreement between the solid vs dashed lines indicate the success of the galaxy environment criterion.

Figure 6 .
Figure 6.Here we assess the cosmic variance in the BH number densities.The light green dashed line corresponds to a BRAHMA-9-D3 box that applies the SM5_FOF3000 model using gas-based seeding with M DGB seed = 2.3 × 10 3 M ⊙ .The remaining solid lines correspond to a set of the secondary BRAHMA-9-E4 boxes with six different ICs.Only the black line applies the same IC as the BRAHMA-9-D3 box.In these boxes, we apply the stochastic seed model with M ESD seed = 1.8 × 10 4 M ⊙ .We find that for the 9 Mpc boxes, the number densities can vary up to factors of ≲ 2 at z ≲ 15 due to cosmic variance.We expect this variance to be smaller in the larger volume boxes.
The model parameters for the stochastic seed model, calibrated for each of the gas-based seeding parameters (SM5_FOF1000, SM5_FOF3000, SM150_FOF3000 and SM5_FOF10000) applied to the BRAHMA-9-D3 suite.Columns 1 to 3 show the gas-based seeding parameters Mh and Msfmp and their corresponding labels.For each set of Mh and Msfmp values, the remaining columns list the parameters of the stochastic seed model.Columns 4 to 8 show the parameter values used for the galaxy mass criterion, which are derived from gas-based seed model predictions of the M galaxy total

Figure 9 .
Figure9.Evolution of number densities of ≥ 1.2×10 6 M ⊙ traced to z ∼ 0 by the secondary boxes (dashed lines).The lenient seed models (light and dark green color) are run on BRAHMA-18-E5.The strict seed models (pink and maroon dashed color) are run on BRAHMA-9-E4.The black markers show observational constraints at z = 0 from C10, M08 and S20.The two thin solid lines are predictions from TNG100 and TNG300.The different seed model predictions show significantly smaller differences at z ∼ 0, compared to z ≳ 7.All of our seed model predictions approach the local measurements at z ∼ 0.

Figure 10 .
Figure10.Black hole mass functions (BHMFs) from our primary set of boxes i.e.BRAHMA-9-D3, BRAHMA-18-E4 and BRAHMA-36-E5 at z ≥ 7. The grey points show observational constraints at z ∼ 0. The three boxes together probe the BHMFs between ∼ 10 3 − 10 7 M ⊙ at z ≥ 7, wherein the BHMFs become steeper for more restrictive seed models.As a result, more massive BHs exhibit stronger differences in their abundances between the different seed models.Lastly, the z ≥ 7 BHMFs are substantially steeper than the local constraints for all the seed models.

Figure 11 .
Figure 11.Here we use the secondary boxes to continue tracing the evolution of the BHMFs to z ∼ 0. The lenient models applied to BRAHMA-18-E5 are shown by the dark and light green dashed lines.The strict models applied to BRAHMA-9-E4 are shown by the maroon and pink green dashed lines.The two thin solid lines are predictions from TNG100 and TNG300.The grey points show the local measurements.The blue circles show the very recent BHMF constraints at z ∼ 5, inferred from the recent faint AGN observations at z ∼ 5 (referred to as "little red dots",Matthee et al. 2023).As we approach z ∼ 0, we note the following: 1) the differences between the various seed model predictions become smaller, and 2) the BHMF slopes become flatter due to formation of more massive BHs contributed by gas accretion at z ≲ 3. The BHMFs approach the local constraints at z ∼ 0.

Figure 13 .
Figure 13.Predictions for the AGN fractions at z = 7 as a function of BH mass for the different seed models from the BRAHMA-9-D3 (left) and BRAHMA-18-E4 (right) simulations.The different rows correspond to different bolometric luminosities.The fractions are exactly zero for the > 10 42 erg s −1 AGN in BRAHMA-9-D3 boxes (lower left panel), because there are no such AGNs formed due to limited volume.Generally, we see lower AGN fractions for the more lenient seed models.Therefore, even if we produce more BHs with more lenient seed models, only a certain fraction of them end up in environments that support enough gas accretion to become AGN.This explains why the seed model variations in the LFs are smaller than that in the BHMFs.

Figure 15 .
Figure 15.Stellar mass vs. black hole mass (M bh − M * ) relations at z ≥ 7 for the BRAHMA-9-D3, BRAHMA-18-E4 and BRAHMA-36-E5 simulation boxes.The black dashed lines correspond to M * = M bh .More specifically, we show the total stellar mass vs. total BH mass within subfind-galaxies.The different columns the show different seed models.The black points correspond to spectroscopically confirmed JWST AGN at their respective redshifts.Grey points in the bottom panel indicate all spectroscopically confirmed JWST AGN at z ∼ 4−7.The solid orange lines show z ∼ 0 measurements.Solid blue line roughly outlines the scatter of Reines & Volonteri (2015), and the dashed and dotted blue lines are mean trends of Kormendy & Ho (2013) and Terrazas et al. (2016) respectively.All the seed models broadly predict overmassive BHs at z ≥ 7 compared to the local relations.However, at fixed stellar mass, BH masses are smaller for more restrictive seed models.Finally, the M bh − M * relations shift between z = 10 to z = 7, showing that galaxy growth is faster than BH growth across this redshift range.

Figure 16 .
Figure 16.Similar to the previous figure, but here we trace the evolution of the M bh − M * relations to z ∼ 0 using the secondary boxes.The left two panels show the BRAHMA-18-E5 predictions for the lenient seed models.The right two panels show the BRAHMA-9-E4 predictions for the strict seed models.At z = 5, the largest BHs in BRAHMA-18-E5 have masses similar to the JWST z = 4 − 7 AGNs (grey points).Our simulations approach the local measurements by z ∼ 0.

Figure 17 .
Figure17.Distributions of the neighbor search radii (R Hsml ) of the BH binaries at the time when the simulations merge them.We show the larger R Hsml amongst the two merging BHs.These distributions essentially probe the typical separations up to which the BH binary evolution can be tracked in our simulations.Left to right panels show increasing BH mass thresholds for the merging BHs, similar to the previous figure.The different colors correspond to different seed models.For a given seed model, lines of different thicknesses distinguish the BRAHMA-9-D3, BRAHMA-18-E4, BRAHMA-36-E5 boxes.Typical separations at which the BHs merge range from ∼ 1 − 15 kpc and essentially depend on the overall resolution and "local" resolution at the merging sites.Therefore, the merger rates in our simulations essentially probe the formation rates of BH binaries at ∼ 1 − 15 kpc separations.

Figure A1 .
Figure A1.Here we compare the contributions of heavy mergers versus light minor mergers (as defined in Appendix A) to the merger driven BH growth within the BRAHMA-9-D3 boxes.The green lines show heavy mergers where the masses of both primary and secondary BHs are ≥ M ESD seed .The orange lines show the light minor mergers where the secondary BH mass is < M ESD seed but the primary BH mass is ≥ M ESD seed .The olive lines show the total contribution from both types of mergers; i.e. all mergers with primary BHs ≥ M ESD seed .The first four columns from the left show results for M ESD seed = 1.8 × 10 4 M ⊙ for all four seed models we considered.The rightmost two panels show results for M ESD seed = 1.5 × 10 5 M ⊙ ; we show this only for the two lenient models for which the stochastic seed criteria could be robustly calibrated; i.e.SM5_FOF1000 and SM5_FOF3000.The top panels show the total merger rate.The middle panels show the mass growth rate due to mergers as a function of redshift, which is defined as the total mass of all merging secondary BHs per unit redshift.The light minor mergers dominate the contribution to the BH mass growth for the majority of the time, and they need to be explicitly included in the stochastic seed models.Notably, the relative dominance of the light minor mergers does tend to recede with time, as quantified in the bottom panels that show the mass growth (∆M light minor ) due to the light minor mergers between successive heavy mergers.Except for the most lenient SM5_FOF1000 seed model with M ESD seed = 1.8 × 10 4 M ⊙ , it is not possible to do a robust power-law fit for ∆M light minor due to the very large uncertainties (green dotted line in the bottom left panel).Therefore we decided to adopt a constant ∆M light minor = 4 M ESD seed (blue horizontal line) for the stochastic seed models, as it is broadly consistent with the ∆M light minor results for all the seed models.

Figure A2 .
Figure A2.Black points show BHMFs for the SM5_FOF1000 model using the BRAHMA-9-D3 suite.We compare the BHMFs produced by the constant value assumption (green) vs the best fit power-law assumption (maroon) for the modeling of the light minor mergers within the stochastic seed models applied to BRAHMA-9-E4 boxes.The best-fit power-law assumption uses ∆M light minor described by the the dashed green line in the left-most bottom panel of FigureA1.The constant value assumption uses ∆M light minor = 4 M ESD seed .While the best-fit power law assumption naturally produces BHMFs close to BRAHMA-9-D3, the constant value assumption only marginally overestimates (≲ 2 factor) the BHMFs.Therefore this assumption does reasonably well given its simplicity.

Figure A3 .
Figure A3.Top panels: BHMFs predicted by the SM5_FOF3000 model when applied to the zoom region used in Paper I.ZOOM_REGION_z5-D3 (green), ZOOM_REGION_z5-E4 (brown) and ZOOM_REGION_z5-E5 (grey) correspond to the same zoom region at the resolutions of BRAHMA-9-D3, BRAHMA-18-E4 and BRAHMA-36-E5 respectively.The three different resolutions use the gas-based seed model to seed 2.3 × 10 3 M ⊙ , 1.8 × 10 4 M ⊙ and 1.5 × 10 5 M ⊙ DGBs respectively.The BHMFs predicted by all three zooms are similar at the most massive end.Bottom panels: The thin green solid lines correspond to BHMFs of BRAHMA-9-D3.For the BRAHMA-18-E4 and BRAHMA-36-E5 boxes, we compare the BHMFs predicted by the stochastic seed model with 1.8 × 10 4 M ⊙ and 1.5 × 10 5 M ⊙ ESDs (thicker green lines), against that of the gas-based seed model with 1.8×10 4 M ⊙ and 1.5×10 5 M ⊙ DGBs (brown and grey lines respectively).Very importantly, the stochastic seed model included the minor mergers using the constant value assumption.Amongst the BRAHMA-18-E4 and BRAHMA-36-E5 boxes, the gas based seed model and stochastic seed model predictions are similar, particularly at the massive end of the BHMFs beyond what can be probed by BRAHMA-9-D3.This establishes that the differences amongst the predictions of BRAHMA-9-D3, BRAHMA-18-E4 and BRAHMA-36-E5, are largely due to cosmic variance.Therefore, the constant value assumption does not compromise our BHMF predictions, even for BH masses larger than what can be effectively probed by the BRAHMA-9-D3 boxes.

Figure B1 .
Figure B1.Here we test the impact of removing the BH repositioning scheme on the merger rates.We consider two BRAHMA-4.5-D3boxes with the SM5_FOF1000 seed model.The solid lines correspond to the box where BH repositioning is applied.Dashed lines correspond to the other box where repositioning is removed and the BHs simply follow the natural dynamics (limited by the simulation resolution).Different colors show the merger rates of BHs above different mass thresholds.The removal of BH repositioning delays the mergers, such that it is reduced by factor of ∼ 3 at z ≳ 10.

Table 1 .
Full simulation suite: L box is the box length in Mpc (2nd column).