Cosmological simulations of massive black hole seeds: predictions for next generation electromagnetic and gravitational wave observations
Colin DeGraf Debora Sijacki



2019

Time evolution of correlation functions in quantum many-body systems
    Luis Pedro García-Pintos
    Submitted to MNRAS
====================================================================


We study how statistical properties of supermassive black holes depend on the
frequency and 
conditions for massive seed formation in cosmological simulations of structure
formation. We develop a novel method to re-calculate detailed growth histories
and merger trees of black holes within the framework of the Illustris
simulation for several seed formation models, including a physically motivated
model where black hole seeds only form in progenitor galaxies that conform to
the conditions for direct collapse black hole formation. While all seed models
considered here are in a broad agreement with present observational
constraints on black hole populations from optical, UV and X-ray studies, we
find they lead to widely different black hole number densities and halo
occupation fractions which are currently observationally unconstrained. In
terms of future electromagnetic spectrum observations, the faint-end quasar
luminosity function and the low
mass-end black hole-host galaxy scaling relations are very sensitive to the
specific massive seed
prescription. Specifically, the direct collapse model exhibits a seeding efficiency which decreases rapidly with cosmic time and produces much fewer black holes in low mass galaxies, in contrast to the original Illustris simulation. We further find that the total black hole merger rate varies by more than one order of magnitude for different seed models, with the redshift evolution of the chirp mass changing as well. Supermassive black hole merger detections with LISA and International Pulsar Timing Array may hence provide the most direct means of constraining massive black hole seed formation in the early Universe.


quasars: general — galaxies: active — black hole physics
  — methods: numerical — galaxies: haloes






§ INTRODUCTION



It is widely understood that supermassive black holes (SMBHs) are found at the centre of massive galaxies <cit.>, and that properties of the host galaxy strongly correlate with black hole mass, suggesting a causal link between them <cit.>.  

Observations have also shown that high-redshift quasars are found at z ∼
6 <cit.> and even z > 7, with black hole masses
estimated at ∼ 10^9 M_⊙ <cit.>.  Reaching
such high masses at such an early cosmic time provides a challenge to the
theoretical models of supermassive black hole growth, and implies the
possibility of the formation of sufficiently massive seeds followed by growth
capable of reaching 10^9 M_⊙ in less than 10^9 yrs of cosmic time.
Cosmological simulations have shown that extended periods of Eddington growth
can be sustained by high density gas flows in very over-dense regions of the
Universe, providing an almost constant fuel source for the accreting black
hole, such that ∼ 10^9 M_⊙ black holes are produced by z ∼ 7
when starting with seeds of ∼ 10^5 M_⊙ <cit.>.  

The initial formation of supermassive black holes seeds, however, is an open
question, with many different mechanisms proposed which may lead to the
formation of light (stellar), intermediate and massive seeds <cit.>. Light seed formation has been proposed in the context of remnants of Population III (PopIII) stars <cit.>.  After the first generation of stars forms in pristine haloes at very high redshift, stars with masses from ∼ 25-140 M_⊙ and ≳ 260 M_⊙ are expected to collapse to form a black hole <cit.>, producing light seeds with M_BH,seed on the order of half the stellar mass <cit.>.  While seed formation from PopIII stellar remnants is potentially common (though with some uncertainty in the mass scales), the low-mass nature of these seeds would require substantially more growth and possibly super-Eddington accretion to reach ∼ 10^9  M_⊙ by z ∼ 7 making them potentially problematic as the progenitors of bright z ≳ 6 quasars.

At an intermediate seed mass scale, it is proposed that black holes may form
in dense nuclear star clusters (NSC), where runaway stellar collisions can
produce a very massive star which can then collapse to an intermediate mass
black hole, on the order of M_BH,seed∼ 10^3  M_⊙
<cit.>.  The higher initial
mass compared to PopIII remnants makes the prospects of reaching
bright quasar masses at z ∼ 6 more viable but not without a challenge of
sustained and very high growth rate.

Massive seed formation is generally proposed to occur in haloes with
T_vir > 10^4 K containing (almost) pristine gas which can potentially
collapse directly to a black hole without significant fragmentation. Such a
collapse requires a massive gas cloud without significant rotation (which
would collapse to a rotationally supported disk rather than a massive black
hole), and without significant cooling from metals or molecular hydrogen, either of which would cause the cloud to fragment rather than collapse to a massive black hole.  
The gas content in primordial galaxies at sufficiently high redshift tends to
be pristine, thus lacking the metals needed for metal-line cooling, and a
nearby source of photo-dissociating Lyman-Werner radiation can disrupt and
prevent formation of molecular hydrogen <cit.>. The net result is a massive cloud capable of collapsing
directly to a black hole seed (referred to as a `Direct Collapse Black Hole',
or DCBH) with mass M_BH,seed∼ 10^4-10^6  M_⊙
<cit.>. Moreover, high baryonic streaming
velocities with respect to the dark matter may promote the growth of DCBHs
<cit.> with alternatively, or additionally
sufficiently high gas turbulent velocities preventing excessive fragmentation as well.  While DCBH seeds can be very massive (and
thus easier to grow to the high masses observed at z ∼ 6-7), the strict
criteria for formation suggests they may be relatively rare to form (see
e.g.  ,  
and recent reviews by   and
 ).   

Given the uncertainty in formation mechanisms and resolution limitations,
cosmological simulations often rely upon a simple prescription to insert a
black hole seed of ∼ 10^5  M_⊙ into sufficiently well-resolved
haloes <cit.>. This method
may be considered to broadly encompass each scenario described above: it is
comparable to a recently formed DCBH, or could be a PopIII or NSC remnant
which formed earlier and grew by several orders of magnitude before being
added to the simulation.  However the number densities of black hole and host
halo occupation fractions could differ largely. Furthermore, the final masses
of the largest black holes have been shown not to be strongly sensitive to the
seed mass provided it is rather large (∼ 10^5 M_⊙ or above), and this
model 
has been shown to well reproduce observed populations of black holes
<cit.>.  However, by construction these
models are unable to directly investigate early black hole formation and
growth. Furthermore, they treat each black hole's seed equivalently, and so
cannot address how sensitive individual black hole histories or time-dependent
statistical properties of the whole population are on the frequency or method
by which supermassive black hole seeds form. In particular, they typically
rely upon two hand-selected mass scales: the threshold halo mass above which
black hole seeds are inserted, and the mass of those seeds.  These simple
assumptions may not be realistic, and suggest that more advanced models for
seed formation should be explored, and there have been several recent analyses investigating this, using both simulations and semi-analytic models <cit.>. 
This is especially timely in view of the
flood of next generation observational data, from, for example, JWST, WFIRST,
eROSITA, Athena, LSST, IPTA, EPTA, NANOGrav, PPTA, and LISA which will push into the presently
inaccessible regime of lower
mass black holes and higher redshifts, together with entirely new constrains
from 
gravitational wave observations. Hence, in this paper we address these
questions by studying how changing the seed frequency and criteria for seed
formation can impact individual black hole growth and the overall populations
of black holes across cosmic time. 

To avoid the prohibitive computational expense of repeatedly re-running a
large cosmological simulation with many varying parameters for black hole
seeding, we developed a post-processing method of re-calculating the growth of
black holes, which we apply to the Illustris simulation <cit.>.
Using the full output of black holes and their surrounding gas properties at
every timestep, we are able to calculate the efficiency with which black holes
should grow, even if the frequency and conditions for seed formation are
changed.  We apply this method to a variety of prescriptions for seed
formation, including a less frequent stochastic seed model, a gas-spin model,
and a physically motivated model based on the mass, metallicity, and gas spin
of progenitor galaxies based on the DCBH mechanism.  For each of these models,
we study which galaxies we expect to form black holes, and how the growth
history of each individual black hole from the original simulation is changed
by each of these models.  Using the full set of re-computed black hole
histories, we characterize the impact on the black hole mass function and
quasar luminosity function, late-time accretion efficiencies, black hole-host
scaling relations, and black hole merger rates, quantifying how the frequency
and mass ratios of black hole mergers are impacted by each seed model.

The outline of the paper is as follows.  In Section <ref> we
discuss the Illustris simulation (Section <ref>) and
associated construction of galaxy merger trees (Section <ref>), our post-processing method of calculating black hole
evolution (Section <ref>), and the seeding models used in our
analysis (Section <ref>).  In Section <ref> we
discuss the results of our post-processing analysis.  Section <ref> shows how the seed model affects seed formation and
the galaxies in which seeds form; Section <ref> shows the
impact on global black hole statistics over cosmic time; and
Section <ref> shows the impact on black hole mergers. Finally, we
summarize our conclusions in Section <ref>. 

 


§ METHOD




 §.§ Illustris Simulation



In this work we analyse the
Illustris[http://www.illustris-project.org; <cit.>.]
suite of simulations with a cosmological box 106.5  Mpc on a side,
performed with 
the moving mesh code AREPO <cit.>. We only consider the highest
resolution simulation incorporating the full physics model <cit.> which has target
gas cell mass m_ gas = 1.26 × 10^6  M_⊙ and dark matter
particle mass m_ DM=6.26 × 10^6  M_⊙.  This run uses a standard ΛCDM cosmology, with Ω_m,0=0.2726, Ω_Λ,0=0.7274, Ω_b,0=0.0456, σ_8=0.809, n_s=0.963, H_0=70.4  km s^-1 Mpc^-1 <cit.>.

For the purpose of this study the modelling of black hole seeding, growth and
feedback is most relevant, so we briefly summarize it here. Black holes are
treated as collisionless sink particles. In Illustris, their seed mass is 10^5 h^-1  M_⊙ and they are seeded in all haloes with mass
above 5 × 10^10 h^-1  M_⊙ if void of a black hole
particle. This seeding model is loosely motivated by the direct collapse scenario
<cit.>, but may remain broadly consistent with lighter seed
formation models 
provided that their occupation fraction is similar and that they grow
efficiently enough in the early Universe. After seeding, black holes may grow
in mass either through gas accretion, which is modelled assuming a
Bondi-Hoyle-like rate <cit.> capped at the
Eddington limit, or by merging with other black holes which happen to be within
each others' smoothing lengths. In the absence of feedback such a black hole
growth prescription would result in unrealistically large black hole
masses. It is hence necessary to incorporate feedback prescriptions
which in Illustris is accomplished by invoking three separate modes:
“quasar”, “radio” and “radiative”, which depend on the accretion
efficiency <cit.>.

We note that while this black hole model leads to a wide range of black hole
properties in good agreement with observational constraints
<cit.>, the seeding prescription is rather simplistic. Thus, the
aim here is to explore different assumptions regarding seeding to pin down how
much they can affect the properties of a representative cosmological black hole
population, and in turn how this can inform us about the likely seeding
scenario. 





 §.§ Constructing galaxy merger trees



The Illustris simulation identifies dark matter haloes using a
friends-of-friends (FOF) halo-finder <cit.> with linking-length of
0.2 and baryonic matter assigned to the same group as the nearest dark
matter particle.  Following this, the SUBFIND algorithm <cit.> is run to produce a catalog of gravitationally self-bound substructures within each halo, which we use to determine galaxy properties.  Furthermore, to track galaxies over time (in particular to identify the progenitors of a given galaxy), we use the SubLink catalog produced by <cit.>.  This catalog identifies descendant subhaloes across snapshots by summing a merit function (based on particle binding energies) for each particle common to a given pair of subhaloes in different snapshots.  The descendant is then defined as the subhalo with the highest summed merit score.  In this manner, we have access to a complete tree of all resolved subhaloes, such that we can identify galaxy mergers and the progenitors of any given galaxy, a key point when investigating if a given galaxy ever satisfied the conditions for SMBH formation at any earlier time.  





 §.§ Post-processing black hole growth



                          
One of the challenges in investigating black hole seeding mechanisms with
cosmological simulations is the large computational expense to run such
simulations.  Ideally, these simulations must have high resolution to model
the regions in which black hole seeding occurs, large volume so as to include
rare objects such as massive high-redshift quasars <cit.>, and be run to low-redshift to provide information on
the SMBH population in the local Universe where observational constraints are strongest (both for predictive purposes, and also to confirm the accuracy of the simulation).  Furthermore, the parameter space to investigate is large, spanning seeding environment (where and when black hole seeds form) and the seed black hole itself[Note that these considerations assume a fixed black hole accretion and feedback model, and that relaxing this would further enlarge the parameter space for exploration.].  
Hence rather than re-running the entire Illustris simulation for each
variation in these parameters <cit.>, which would be numerically prohibitive, we take a
post-processing approach in which we 
re-calculate the entire growth history of the full black hole population,
using the high temporal resolution output from the original run, which
contains local gas properties of each black hole at every timestep.

 Black hole accretion within the simulation is given by

    Ṁ_BH = 4 πα G^2 M_BH^2 ρ/c_s^3

<cit.>, with an imposed upper limit of the Eddington rate [Ṁ_Edd=(4 π G M_BH m_p) / (ϵ_r σ_T c)].
Thus the accretion rate depends upon black hole mass (M_BH) and the
density (ρ) and sound speed (c_s) of the nearby gas, with dimensionless
α = 100, radiative efficiency ϵ_r = 0.2, gravitational
constant G, speed of light c, proton mass m_p and Thompson
cross-section σ_T.  Since the simulation saves each of these properties every time
the black hole properties are calculated (with typical timesteps on the order
of 2 × 10^4 yrs), all the information necessary to re-calculate
accretion rates is saved. 

We perform this growth recalculation by taking the necessary properties (M_BH, ρ, c_s) to calculate the expected accretion rate (Equation <ref>, with imposed Eddington limit).  We repeat this calculation at each timestep to find the expected mass at the following timestep (i.e. M_BH,i+1 = M_BH,i + Ṁ_BH,iΔ t).
In this way we are able to reproduce the full accretion history of the entire black hole population, but with the freedom to change certain parameters (e.g. seeding a subset of the black hole population or changing the initial seed mass).
We also incorporate black hole mergers in our new calculation.  Because we are
using the gas properties local to black holes in the original simulation, our black holes are
necessarily subsets of the original Illustris black holes (though often with different
masses).  Thus the mergers for our new black hole populations are also
subsamples of the mergers occurring in the original simulation (and for each
merger, the simulation saves the merger time and the ID and mass of each black
hole).  We determine mergers in our recalculated history by extracting mergers in the original simulation which involve only black holes seeded in our new model, at which point we combine the masses of the merging black hole pair (using our newly-calculated mass histories).  

We note that black holes in the Illustris simulation merge as soon as they are
within the black hole smoothing length of one another, without considering the
hardening/merging timescale, which may delay when a given black hole pair
should merge.  In this analysis we use the merging times of the original
simulation, but note that an extended hardening times may delay black hole
mergers and thus impact our predicted black hole merger rate.  However, an initial investigation
finds that adding a delay of 1Gyr has minimal impact on the overall merger rates, except at the highest
redshifts where the first black hole mergers occur, consistent with results from <cit.>, who found minimal impact for delay timescales of 5Gyr in gas-poor galaxies. 

We note that the primary limitation of our model is the inability to
incorporate self-consistent feedback and self-regulation into post-processing
calculations.  In the original simulation, AGN feedback is deposited into the
surrounding gas as the black hole grows; in our post-processing analysis we
are able to modify the accretion history of the black hole, but the feedback
energy deposited remains unchanged.  At low black hole masses we expect this
to have a minimal impact, as the black holes are generally too small to have a
strong impact on the gas supply of their host halo.  At high black hole mass,
we know that quasar feedback plays an important role, and acts as a catalyst
for self-regulation of further growth <cit.>.  Thus our recalculated growth histories may
underestimate the final black hole mass, but this will only be relevant for black holes whose recalculated mass is substantially below the original mass at the time when feedback-induced self-regulation begins (we further discuss this issue and potential implications in Section <ref>).  

We further note that the Illustris simulation output is incomplete for black holes below z ∼ 0.38. This missing data necessarily prevents recalculating growth and mergers at lower redshifts, so we limit our analysis to z ≥ 0.5.  This limitation has a minimal impact on our results, however, as it is after the majority of black hole growth.  We do not expect any qualitative changes at lower redshifts.





 §.§ Seeding models





  §.§.§ Illustris seeding



Within our framework, the simplest seeding prescription is to re-use the data from seeding black holes in the original Illustris simulation, i.e. when a halo crosses a specified threshold of 5 × 10^10 h^-1 M_⊙.  By only seeding black holes that were originally seeded in Illustris, we explicitly have the local gas properties necessary for re-calculating growth, and the black hole mergers in this new case will simply match the black hole merger trees from Illustris.  If using the same seed mass as Illustris (10^5  h^-1 M_⊙), this will mean that each black hole should be `re-grown' in a way that matches the data from the original simulation, which we have explicitly confirmed. 



  §.§.§ Stochastic seeding



An obvious variation on using the Illustris seed criterion is to simply seed a
fixed fraction (f_seed) of Illustris-seeded black holes, accomplished
using a random number generator to determine whether a given Illustris seed
should be considered `seedable' in our re-calculation.  As this model
continues to only seed subset of black holes which have been seeded by
Illustris, we again have the complete set of local gas properties.
Furthermore, since all black holes are chosen from the original Illustris
catalog, all black hole mergers can be extracted from the original Illustris
merger trees.






We note that although we seed using a fixed probability, the total number of black holes is not the same fraction of original black holes (i.e. N_BH,f_seed f_seed× N_BH,Illustris), due to black holes which involve multiple black hole seeds that eventually merge together.  Furthermore, this is a mass-dependent effect, since high-M_BH,Illustris black holes tend to have undergone the largest number of mergers and thus have the largest number of chances to have at least one seed in our new calculation (e.g. for f_seed=0.1, at z = 0.5 18% of Illustris black holes above 10^6  M_⊙ have been reseeded in our new model, compared to 58% of Illustris black holes above 10^8  M_⊙. This is discussed further in Section <ref> and especially Figure <ref>). This different impact on high- and low-mass black holes means that the impact on overall black hole populations, such as the black hole mass function and luminosity function as functions of cosmic time, is more complicated than a simple decrease in total black hole number (for further details, see Section <ref>).



  §.§.§ Spin- and metallicity-based seeding



In Figure <ref> we investigate the galaxies seeded by black holes
in the Illustris simulation, showing the distribution in both spin-redshift
space and 
metallicity-redshift space (both important criteria in direct-collapse
formation models), where λ and Z are computed for the gas
within the gas half-mass radius of the galaxy.  This distribution shows that
gas spin in the center of newly black hole-seeded galaxies is essentially
redshift independent, and metallicity only shows a slight increase at low
redshift (below z ∼ 1.5).  However, this is explicitly the case for
haloes which have just recently crossed the 5 × 10^10 h^-1  M_⊙ threshold; we consider the impact on other galaxies in more detail in Section <ref>. 

As a variant on the stochastic seeding model, here we seed a subset of
Illustris black holes based on the dimensionless spin parameter of gas within
the gas half-mass radius of the host galaxy at the time when the Illustris
black hole is seeded <cit.>.  Specifically, we
impose a cut on spin such that only black holes in galaxies with λ <
λ_crit are seeded, where we select λ_crit to
correspond to a given fraction of galaxies (i.e. to be comparable to a given
value of f_seed).  Since λ is approximately
redshift-independent, seeding based on gas spin in this manner is
qualitatively similar to the random seed selection model (note that the 50th-,
25th-, and 10th-percentiles are marked in Figure <ref>,
corresponding to f_seed = 0.5, 0.25, 0.1 in Section <ref>).  As a further test, we imposed an additional criterion based on gas metallicity (see also Section <ref>).  The metallicity distribution is also roughly redshift-independent for z > 2, so this additional criterion tends to only affect a small number of low-redshift, low-mass black holes, and otherwise produces qualitatively similar results.







  §.§.§ Progenitor-based direct collapse seeding



The models described in Sections <ref> and
<ref> are based solely upon instantaneous properties when
a black hole particle was first inserted into the Illustris simulation.
However, we note that the insertion of a sink particle in Illustris (as in
other cosmological simulations) is not necessarily intended to represent the
actual formation of a seed black hole, but rather a point at which we would expect a seed to have previously formed and grown to a “resolvable” mass.
A more physically meaningful approach is to consider which black holes seeded
in Illustris have a progenitor galaxy which satisfies the necessary criteria
for the formation of a black hole seed, in this case by direct collapse.  We
are able to do this by considering the full history of all Illustris galaxies
(and the merger trees involving them; see Section <ref>), checking
for a set of conditions necessary for DCBHs where we impose three criteria:


  
  * Galaxy mass: We impose that M_subhalo > 3 × 10^9  M_⊙ (where M_subhalo is the total mass contained within the SUBFIND-defined subhalo, as discussed in Section <ref>).  This is largely a resolution-based criterion (corresponding to ∼ 500 dark matter particles) such that we have sufficient gas particles to reasonably resolve the remaining criteria.



  
  * Gas spin: Sufficiently high gas spin implies a more gravitationally
    stable structure which may not lead to the further collapse necessary to form a DCBH seed. A rough approximation for this criterion is

    λ < m_d^2  Q_c  T_ vir/8  j_d  T_gas,

where m_d is the fraction of gas involved in infall, Q_c is the Toomre stability parameter, j_d is the angular momentum in the disc, T_vir is the virial temperature of the halo, and T_gas is the gas temperature <cit.>. We calculate T_vir and T_gas directly from the simulation, noting that T_gas is typically on the order of 1-3 × 10^4K.  The remaining parameters (m_d, j_d, and Q_c) are much more sensitive to the structure of the disk which is less well-resolved in our simulation, especially in the low-mass galaxies at high-z relevant for the DCBH formation.  As such, we follow <cit.> and use m_d = j_d = 0.05 and Q_c = 2 rather than use the resolution-sensitive values from the simulation.  



We show an example of these constraints using 1D histograms in
Figure <ref>.  The solid (dashed) black histogram shows the
distribution of λ (λ_max) of galaxies newly-seeded in
the original Illustris simulation, where there is a minor overlap (i.e. λ < λ_max).  For comparison, the green histograms show the distribution for λ and λ_max for all galaxies above 3 × 10^9 M_⊙ (our resolution-based mass threshold for DCBH seeding) at z = 4.  Here we see that both λ and λ_max are higher than for the newly-seeded haloes in Illustris, but again there is a region of overlap where in principal DCBH formation can occur.

  
  * Gas metallicity: In addition to rotational support, one of the main
    obstacles for massive DCBH seed formation is fragmentation of the
    collapsing gas cloud.  Specifically, metals within the gas can facilitate
    cooling, which leads to fragmentation of the cloud prior to collapse into
    a single supermassive black hole seed.  To account for this, we impose a
    metallicity threshold of Z < 10^-5 Z_⊙ for the galactic gas
    <cit.>.  We note that our
    results are not sensitive to the exact threshold, since Illustris galaxies
    rarely have non-zero metallicity near this value (e.g., at z > 5,
    <0.05% of galaxies with non-zero metallicity have metallicity below
    10^-3 Z_⊙). This is also shown in Figure <ref>,
    where the green 1D histogram shows the distribution of metallicities for
    galaxies above 3 × 10^9  M_⊙ at z = 4: the peak of
    non-zero metallicity galaxies occurs at Z ∼ 10^-3.5  Z_⊙, and
    nearly one-third of galaxies have zero metallicity. 




We note
      however that this is sensitive to the metallicity enrichment model of
      the simulation. In the Illustris simulation, star particles are treated
      as single-age stellar populations which deposit metals into the
      neighboring gas cells based on the stellar evolution model, which
      includes supernovae (including core collapse, see
        and
       , and Type Ia, see  ) and AGB winds <cit.> [note that the simulation directly models nine distinct elements, though we only consider total metallicity in this analysis].  The metals are also distributed by stellar feedback via metal-loaded winds, and metal mixing is modeled purely by advection <cit.>.  

Our results may be sensitive to this model as a less efficient metal
mixing model could lead to more haloes satisfying the metallicity criteria;
conversely, a full incorporation of PopIII star formation may contribute more
high-redshift metallicity, thus decreasing the haloes satisfying our metal
criterion as there is no high redshift metallicity floor in Illustris.


These criteria are used to determine which Illustris black holes have progenitor galaxies with the necessary conditions for DCBH seed formation.  Specifically, for the host subhalo of every newly seeded black hole in the Illustris simulation we check every progenitor galaxy in the SubLink galaxy merger tree <cit.>.  Any Illustris black hole which has at least one progenitor galaxy simultaneously satisfying all three criteria is considered a `seedable' black hole according to our DCBH model. We note here that our method of post-processing black hole growth relies on the output data provided for each Illustris black hole. As such, we are unable to directly calculate growth for black holes which do not exist in the original simulation. Thus our seeding model here does not seed black holes as soon as a galaxy satisfies the DCBH criteria, but rather when the original simulation seeds a black hole which had a progenitor galaxy which satisfies the criteria. 





To illustrate our DCBH seeding method, in Figure <ref>
we show the full galaxy progenitor trees of three z = 5 galaxies with
newly-seeded black holes in Illustris, colour-coded by gas spin and symbols
distinguishing between galaxies with M_* = 0 and M_* > 0 (vertical lines
represent individual galaxy evolution with time; horizontal lines represent
galaxy mergers).  The first is a representative moderate-sized galaxy tree
with no seedable progenitors; the second is a very large merger tree with a
single seedable progenitor galaxy (circled in red), and the third shows a
representative small merger tree containing a single seedable progenitor
galaxy (circled in red).


As seen in Figure <ref>, the conditions for DCBH seed formation typically occur significantly earlier than the Illustris seeding.  As such, although we apply our post-processing analysis to a given black hole with an initial seed mass of 10^5 M_⊙ starting when Illustris seeds, there is a substantial period of time since actual seeding (when DCBH conditions are satisfied) during which the seed could grow. For high-redshift seeds, black holes tend to have several hundred Myr in which to grow, introducing two additional parameters to characterize this early growth: the seed mass when the initial DCBH forms, and the typical growth rate from that time until Illustris seeds the BH.  Among z > 5 seeds, we find that if the initial seed mass is 10^4 M_⊙, 96% of DCBHs would be able to reach 10^5 M_⊙ by the time they are included in Illustris if limited to sub-Eddington rates, and 83% would be able to do so even if initially formed at only 10^3 M_⊙.  As such, we find that the seed mass used in Illustris is fully consistent with a wide range of lower formation masses, followed by a reasonable period of sub-Eddington accretion. Recent work by <cit.> investigated this early growth directly, using GADGET-3 simulations which seed at M_BH,seed=10^3  M_⊙ based on local gas properties, and found that the mean mass for black holes in haloes at the threshold for seeding in Illustris is ∼ 10^5.26  h^-1 M_⊙, confirming that a seed mass of ∼ 10^5  M_⊙ is reasonable, and that we can expect our DCBH seeds to reach this mass by the time the Illustris black holes are included.


We further note that an additional constraint for DCBH formation is a nearby
source of photo-dissociating Lyman Werner radiation, capable of dissociating
molecular hydrogen and thus prevent fragmentation resulting from molecular
cooling.  However, the Lyman-Werner intensity threshold needed remains very
poorly constrained <cit.>. Furthermore, recently <cit.> proposed
that lower Lyman-Werner fluxes can still be sufficient, with the primary
driver for DCBH formation being high gas inflow rates onto the halo.  Hence
estimating for the incoming Lyman-Werner flux in the Illustris simulation would
be highly uncertain, given the resolution limitations and the lack of
PopIII stars.  Given these caveats, we do not include any
Lyman-Werner constraint in our analysis here.  Note that we expect our seed
model to over-estimate the DCBH formation rate compared to a simulation which
fully incorporates this additional criterion, and thus may be considered an
upper limit on the formation efficiency.




This DCBH seeding model is the most physically meaningful of all the models we consider, since we only take into account black holes in galaxies which have, at some point in the past, satisfied the necessary criteria for DCBH formation.  Because our post-processing growth calculations rely upon the local gas properties which are only saved around Illustris black holes, this model by construction only includes black holes which are also seeded in the original Illustris simulation.  As such, we again have the complete history of any mergers involving any combination of black holes seeded according to this new model.  Note that this approach ignores galaxies which satisfy the criteria for DCBH seed formation but which never have a black hole form in Illustris.  However, the only cases in which this can occur are galaxies which satisfy the above criteria but 
whose host haloes do not cross the 5 × 10^10 h^-1 M_⊙ mass threshold for seeding.  These black holes would likely not grow into the regime in which we consider our results well-resolved (M_BH > 10^6  M_⊙), and so will not largely impact our results.  






We note that we are only able to check progenitor galaxy conditions in the saved snapshots, which sets the time-resolution for this model. This timescale (∼ 40 Myr for z ∼ 5-10) is longer than the typical timescale for DCBH formation, such that it is possible to miss a galaxy which satisfies the criteria during a short time period between snapshots.  
The typical variation of λ is relatively slow compared to the snapshot timescales (see the generally gradual colour variation in individual vertical lines in Figure <ref>), so we would expect this to be a sub-dominant effect.  The metallicity evolution can be fairly rapid, however, suggesting more frequent snapshots might show additional galaxies which satisfy the conditions for DCBH formation.  By sub-sampling the available snapshots and extrapolating to shorter inter-snapshot timescales, if we assume that the collapse time is on the order of ∼ 5 Myr, it is possible we might underestimate the seed frequency by up to 50%.  We expect this is a conservative upper limit based on significant extrapolation (and the unresolved metal enrichment); nonetheless, we would expect that even an increase at this upper limit should not qualitatively change the results presented here.  A more precise treatment would require a full re-simulation which checks for seed conditions on-the-fly, which is beyond the scope of this analysis.  








  §.§.§ Individual black hole growth





In Figure <ref> we illustrate our post-processing
technique, showing the full history of a SMBH and all its progenitors using
the four different seeding models described above. Each panel shows the mass
history for the most massive black hole from the original Illustris simulation
as well as all black holes it merged with, having re-calculated the history
according to each seed model.  Over-plotted in colour on each panel is the most
massive progenitor history for each of the four models, determined by tracking
the final black hole back in time, and at each black hole merger selecting the
more massive progenitor black hole at the time of the merger (and as such the
coloured lines are identical in all four panels). There are significant
differences in the total number of seeded black holes and the masses such
black holes reach at both early and late times.     


The top panel shows the history for the f_seed=1 model, i.e. the
seeding conditions as in the original Illustris simulation.  This seed model
involves a 
large number of black hole mergers, with black holes being seeded continually throughout cosmic time.  The second panel shows the f_seed=0.1 model.  By construction this model exhibits an equivalent seeding pattern as the original simulation, merely more rare.  This shows that forming only a subset of the original seed population can have a significant impact on the overall black hole growth history, affecting not only the seed number but also how efficiently the most massive progenitor grows (comparing the dark green and blue curves), and even the final mass (in this case resulting in a black hole ∼0.5 dex less massive).  In the third panel, we show seeding using a cut on the instantaneous gas spin and metallicity when Illustris first seeds black holes, as described in Section <ref>.  Again we see a moderate impact on the final black hole mass, as well as a significant impact on the growth histories reaching that point.

Finally, in the bottom panel of Figure <ref> we show
the analogous history according to the DCBH seeding model.  Again we see a
moderate impact on the final M_BH (though weaker than seeding with
∼ 10% efficiency). We furthermore note a dramatic impact on when DCBH
seeds form: rather than forming continuously throughout the simulation, the
seeding occurs preferentially at high redshifts (z > 4) when gas clouds are
more likely to consist of pristine gas capable of collapsing directly to a
massive black hole without fragmenting.  This example growth history
demonstrates the impact the seed formation can have on black hole populations,
which we investigate quantitatively throughout the remainder of the paper.




§ RESULTS





 §.§ Black hole seeds and their host galaxies



Having performed the full recalculated growth histories of black holes
according to the models described in Section <ref>, we
consider the galaxies in which seeds form in each model.  In
Figure <ref> we show the distribution of dark matter mass
(M_DM), gas mass (M_gas), stellar mass (M_*), and baryon
fraction (f_baryon) of galaxies with newly-seeded black holes as a
function of redshift, for the original Illustris seeding model (right panels)
and our DCBH-based seeding (left panels). The most obvious differences between
the seeding models is that the original Illustris seeding tends to take place
in galaxies with higher masses and at lower redshifts than our DCBH-based
seeding.  This is to be expected, since the idea behind Illustris seeding (and
that of most cosmological simulations) is to insert a black hole roughly
consistent with the black hole - host galaxy scaling relation possibly well
after the actual seed should have formed.  In contrast, our new DCBH model
finds galaxies in which the conditions for DCBH formation are satisfied, which
tend to occur at earlier times and in lower-mass galaxies <cit.>. 










Furthermore, we note that the masses and baryon fractions of galaxies hosting
Illustris seeds tend to be relatively redshift independent (though M_* does
increase slightly with time), since Illustris always seeds haloes as they
cross a fixed halo mass threshold (M_halo > 5 × 10^10  h^-1 M_⊙). In contrast, in our DCBH-based seeding we find
that at low redshift (below z ∼ 3), only galaxies with relatively low
baryon fractions are seeded, as they are more likely to have the
low-metallicity necessary to prevent excessive fragmentation of gas.  

The top left-hand panel of Figure <ref> shows the total number of galaxies and haloes above
several mass thresholds as a function of redshift (galaxies - solid lines; haloes - dot-dashed lines), showing that
the number of high-mass galaxies grows monotonically with time, and the number
of low-mass galaxies peaks at moderate redshift (below which we find a slight
decrease due to increasing number of galaxy mergers).  In particular, we note
that the dot-dashed blue line (> 5 × 10^10 h^-1 M_⊙) is
the number of haloes above the mass threshold for black hole seeding in the
original Illustris simulation, and thus corresponds roughly to the number of
black holes seeded in Illustris (neglecting mergers).  We compare
these total halo numbers to the number of galaxies satisfying our criteria
for DCBH seed formation at redshift z (dotted line), and the cumulative
number of galaxies satisfying the same seedability criterion at any earlier
time (dashed line).  This shows that the total number of seeded black holes by
z = 0.5 is comparable between the original Illustris simulation (blue line)
and our DCBH-based model (dashed black line); however, where the original
simulation has a monotonically increasing seeded number, our DCBH seed model
peaks at z ∼ 8 where the conditions for direct collapse are most common.
Furthermore, we note that our total number of DCBH seeds is actually higher than the original simulation.  This increase
is due to the lower host galaxy mass (see also Figure <ref>), and thus represents a population of black holes in lower-mass galaxies than included in the Illustris black hole population.  Given how little the scaling relation is affected, this would represent a correspondingly low-M_BH population and would thus have minimal impact on the population of resolved black holes we consider in this analysis (see also discussion in Section <ref>).

From this figure, we find a very large number of galaxies which could
potentially form DCBH seeds.  Based on our model criteria (see
Section <ref>), we predict ∼ 10^-2  cMpc^-3 possible DCBH seed candidates by z ∼ 6, far above the
number density of high-redshift quasars <cit.>. However, we note that the majority of these objects will not
grow significantly.  In fact, at z ∼ 5, we find ∼ 9 × 10^-4cMpc^-3 possible DCBH seeds whose host halo will have reached the
mass threshold for seeding in Illustris, and our regrowth calculations suggest
only ∼ 10^-5cMpc^-3 will have reached M_BH = 10^7
M_⊙. This is still well above the observed quasar number
density, so incorporating additional seed criteria such as a Lyman-Werner
radiation threshold or a gas inflow cut <cit.> would be a
worthwhile exercise, bearing in mind a possibility that only a fraction of
DCBH seeds grows sufficiently to power high redshift quasars.

In the lower left-hand panel of Figure <ref> we plot the
fraction of 
galaxies at redshift z (with M_gal > 3 × 10^9  M_⊙) which would have formed a DCBH seed at some earlier time,
according to our DCBH model.  We find that the seeded fraction reaches a peak
of ∼ 25% at z ∼ 8, while at lower redshifts the fraction drops,
suggesting only ∼ 10% of local galaxies above 3 × 10^9 M_⊙ should have a black hole formed
through direct collapse. We also note a slight upturn in the seeded fraction
below z ∼ 3.  However, this increase is not due to more direct collapse
seeds; rather it is a result of the decrease in N_gal (M > 3 ×
10^9  M_⊙) caused by low redshift galaxy mergers (seen in the
red curve of the upper panel). 




In the right-hand panel of Figure <ref> we plot the fraction
of original Illustris black holes which are seeded according to several
alternate seeding models (i.e. N_BH,new/N_BH,Illustris).  In
particular, we consider seeding based on gas spin only (see
Section <ref>), having selected a maximum gas spin value
corresponding to approximately 
50, 25, and 10 per cent of Illustris black hole seeds
(λ_g,50, λ_g,25, and λ_g,10, respectively),
and for our DCBH seed model.  We also
show the total number of black holes seeded in Illustris for reference.  Consistent with Figure <ref> which showed λ had
minimal z-dependence, seeding according to instantaneous galaxy spin alone is
roughly redshift-independent, similar to the simpler
f_seed-based seeding (not plotted, but would correspond to horizontal
lines by construction).  In contrast to this, DCBH seeding grows with time
until z ∼ 4.  Below this the fraction of seeded black holes decreases as
few new DCBH seeds form and existing black holes proceed to merge together.  


This peak at z ∼ 4 is qualitatively different from the left-hand panel of
Figure <ref>, which demonstrated that the peak in DCBH seeding
occurs at z ∼ 8.  This difference is in the definition of what is being
plotted: left-hand panel shows the fraction of galaxies which satisfy the
criteria for seeding a black hole at redshift z, while the right-hand panel shows the fraction of Illustris black holes which have at least one progenitor which satisfied our seeding criteria at some earlier time.  At lower-redshift, a given Illustris black hole has a larger number of progenitor galaxies than at higher redshift, and is thus more likely to have at least one progenitor satisfying the necessary criteria, hence the increase in seedable fraction.  At lower redshifts (z ≲ 4), however, most of the progenitor galaxies
capable of forming a DCBH have already merged into a large enough halo for
Illustris to insert a black hole; thus most black holes inserted into Illustris at at z ≲ 4 tend not to have any DCBH-seedable progenitor galaxies making low-z DCBH formation very rare. 


In Figure <ref> we show the distribution of black hole mass at z = 0.5 as a function of redshift at which the black hole is seeded (z_seed) for the original Illustris simulation (top) and our DCBH seed model (middle).  We define z_seed as the seed redshift of the most massive progenitor within the given black hole merger tree, rather than the earliest seed redshift within the tree, noting that the most massive progenitor is often not the first to be seeded (e.g. see Figure <ref>), but is the most relevant quantity to consider here.

We find a direct correlation between M_BH(z=0.5) and z_seed,
consistent between the original Illustris simulation, our DCBH seeding, and our f_seed-based models (not plotted).  
As previously discussed, there are few low-z_seed black holes in the DCBH model.  This is further quantified in the bottom panel of Figure <ref> where we plot the fraction of original Illustris black holes which should be seeded according to our DCBH model.  This shows a rapid decline in the fraction of Illustris black holes with z_seed≲ 4 which should form via our DCBH model.  



 §.§ Global black hole properties



Having considered the galaxies in which black holes should form, we use the
full recalculated growth history of all black holes to investigate the impact
the seeding model has on measurements of the full black hole population. In
particular, in Figure <ref> we show the black hole mass function
(BHMF) and
quasar luminosity function (QLF) from the original Illustris simulation, for our
fractional seed models and for our DCBH-based
seeding model. We also show observational data from <cit.> and <cit.>, and the
sensitivity limits of the Athena <cit.> and eROSITA
<cit.> surveys.  For the <cit.> data, the Athena and
eROSITA limits, we convert from X-ray to bolometric luminosities using the SED
of <cit.>.

The BHMF for the fractional seed models shows generally expected behavior: as
f_seed decreases, the BHMF decreases.  At high redshift, this shift
in BHMF is proportional to the decrease in f_seed, as the lower
f_seed decreases the number of seeded black holes without affecting
the growth of those black holes which are seeded.  At low redshift, however,
black hole mergers become significant for high-mass black holes, and thus the
post-seed growth behaviour depends on seeding probability as well.  Since
massive, low-redshift black holes tend to have a large number of progenitor
black holes, most of the high-M_BH black holes are seeded regardless
of f_seed. 

We show this explicitly in Figure <ref>, the upper panel of which
shows the correlation between the mass of each black hole at z = 0.5,
M_BH(z = 0.5), and the total number of progenitor seeds it had,
N_prog, while the lower panel shows the fraction of black holes which
have at least 2, 4, or 10 progenitors (corresponding to f_seed = 0.5,
0.25, 0.1, respectively). This shows that <N_prog >
grows with M_BH, and the mass scale above which we would expect most
black holes to have at least one progenitor seed would be M_BH(z =
0.5) ∼ 1.5 × 10^7  M_⊙, 6 × 10^7  M_⊙, and 4 × 10^8  M_⊙ for f_seed
= 0.5, 0.25, and 0.1, respectively.  We further note that these scales are
redshift-dependent, as higher redshift black holes tend to have fewer
progenitors. In addition to the original Illustris data, we show the number of
progenitor black holes in our DCBH seeding model using dashed lines.  As
expected, N_prog is lower in this model, as fewer black holes are
seeded. 

The decrease in M_BH among high-mass black holes is due to the
decreased mass gained via mergers (since fewer progenitors were seeded).  This is further compounded by a lower accretion rate, both because M_BH∝
M_BH^2 (Equation <ref>), and the higher accretion rate in the original simulation produced stronger feedback and thus lower gas density and higher temperature than we should have in the new seed model, which each contribute to a lower Ṁ. Furthermore, the most massive
black holes require a relatively early seed in a host with sufficient gas
density to fuel an extended phase of near-Eddington growth (e.g. the most
massive progenitor curves in Figure <ref>).  When the
seed fraction is decreased, there are fewer chances and thus a generally
delayed onset of Eddington growth, producing a lower black hole mass.  The net
effect is that the number of high-mass black holes tends to be roughly
unchanged but with slightly lower mass (as seen in
Figure <ref>), producing a slightly steeper slope in
the BHMF. We further note that this also means that although the high-end mass function is lower in our DCBH model, the black hole occupation fraction is still high in massive galaxies (rapidly approaching 1 for total galaxy masses ∼few× 10^12 M_⊙). Conversely, low mass black holes with single progenitors are less likely to be seeded, but those that do have unchanged masses. Thus the decrease in the low-end BHMF does represent a comparable decrease in the black hole occupation fraction among low-mass galaxies as has also been found to be the case when seeding at lower masses based on local gas properties <cit.>.  

The BHMF for the DCBH model is roughly consistent with the f_seed∼
0.25 model, except for the low-redshift, low-M_BH end, where the
DCBH seed model falls off drastically (e.g. below M_BH∼ 10^7 M_⊙ at z = 2) due to the difference in redshifts at which black hole
seeds form.  In particular, most seeds in the DCBH model form at an early time
(see Figure <ref>), in contrast to f_seed based
models which continue seeding black holes at all redshifts as is the case in
the original Illustris model.  Since the DCBH
model seeds very few black holes at low redshift and M_BH(z = 0.5) is
directly correlated with z_seed, nearly all black holes have had
sufficient time to reach at least ∼ 10^6.5 M_⊙ (given
our accretion prescription), and
thus the mass function decreases rapidly below this scale. Thus our model
predicts very few low redshift, low mass black holes which formed via direct
collapse route. However, we note that this does not necessarily mean that there
will be a turnover in the overall observed BHMF. Our model only considers
DCBHs, but we also expect black holes may form via NSCs or from PopIII
remnants. Black holes formed through these mechanisms may contribute
significantly to the mass function, providing the potential to partially or
even completely overcome the turnover found in the DCBH
population. Alternatively, or additionally, it is possible that some of
the direct collapse seeds will not grow by accretion as efficiently as in our
model which may lead to a less significant decrease in the low mass-end BHMF
shape. Furthermore, our analysis does not include black holes which may have been seeded by direct-collapse but whose host halo never reaches the threshold for seeding in the original simulation (5 × 10^10  h^-1 M_⊙).  As discussed in Section <ref>, this will represent a low-M_BH population, and thus should have minimal impact on the results shown here.  In particular, we note that the mean mass for black holes in haloes below 10^11 M_⊙ in the original simulation is only 2.6 × 10^5  M_⊙ at z = 0.5 and less than 3% have M_BH > 10^6  M_⊙, suggesting these black holes will have minimal impact on the results considered here, where we generally limit our analysis to M_BH > 10^6  M_⊙.  This is supported by results of <cit.>, who found that seeding black holes at earlier times with M_BH,seed=10^3   h^-1 M_⊙ resulted in comparable mass black holes when reaching the Illustris-seed threshold (for haloes at M_DM=10^10  h^-1 M_⊙, ⟨ log(M_BH) ⟩∼ 5.18 with standard deviation ∼ 0.55).  <cit.> also looked at low-M_BH,seed simulations, and found that the black holes would not grow to 10^6  M_⊙ until the host galaxy reached a stellar mass of a few × 10^8  M_⊙ (and a few × 10^9 in simulations which used a delayed cooling model for SN feedback), broadly consistent with Illustris' halo mass threshold.



Similar to the BHMF, the QLF tends to exhibit a simple normalization shift
equivalent to f_seed at either high redshift or low-L_BH,
the regimes in which most black holes have undergone few or no mergers and
thus modifying the seed prescription changes which black holes form but
does not affect their subsequent growth.  As with the high-mass BHMF, the low
redshift, high-luminosity black holes tend to have undergone numerous mergers,
and thus the number is generally unaffected by the seed model.  However, we
find that the low-redshift, high-end of the QLF is somewhat less affected by
the seed 
prescription than the BHMF is, such that current observational data cannot
differentiate between our seed formation models. This is broadly consistent with
the semi-analytic model of <cit.>, who also
found the mass function to be more sensitive to seed model than the luminosity
function.


We note that we find good agreement with observational data across a wide
range of luminosities and redshifts.  At low redshift and high luminosities,
all of our seed models reproduce the observed data well (though noting that
the simulations are
volume-limited, and thus do not produce the brightest quasars). At the
lowest luminosities the rarer seed models under-predict the luminosity
function, but additional seed models which provide lower-mass and
lower-redshift seeds would be expected to fill this gap in.  At moderate
redshifts (z ∼ 2) we find that the rarest seed models (including the DCBH
model) match the faint end most effectively. Although the brightest objects
are underestimated, we expect this to be at least partially due to the limited
boxsize and somewhat overestimating the impact of feedback in our DCBH seed
model (see the discussion of the scaling relation in
Figure <ref>, below).  At the highest redshifts (z ∼
4), all our models remain broadly consistent with observations, though this
is largely due to the uncertainty in the observed high-redshift faint-end
luminosity function. We note that the upcoming eROSITA data will not help
differentiate between the different black hole seeding models (see redshift
dependent eROSITA limits), which will only become possible (both at high and
low redshifts) with the Wide Field Imager of the Athena mission.



In Figure <ref> we show the total comoving number density
of AGN binned by bolometric luminosity to investigate cosmic downsizing,
comparing the original Illustris simulation (dashed lines) to the DCBH model
(solid lines).  We also show observational constraints in the form of shaded
regions, showing the areas spanned by the data from <cit.> and a
range of models fitted to observational data <cit.>.  We
find that Illustris overpredictss low-luminosity AGN
(10^43-10^44 erg s^-1, blue) at high-redshift, whereas the DCBH
model more closely matches the observed data.  Similarly, at intermediate
luminosities (10^44-10^45 erg s^-1, green) the DCBH model
matches observations very closely, where the original simulation over-predicted
the number density, especially at high-redshift.  In the brightest sample
(10^45-10^46 erg s^-1, red) the DCBH under-predicts the number
density at all but the highest redshifts.  However, we note that the brightest AGN are least
well-constrained in our model, as they are the objects where the AGN feedback
(which is unaffected by our DCBH model) plays the strongest role, which we can better understand by looking at the relation between black holes and their host galaxies.
Furthermore, we emphasize that the comoving number densities from <cit.> are based on analytic fits to the observed data.  The observed data includes significant scatter (see Figure <ref>) and can vary significantly depending on the fitting model used, especially when extrapolating toward faint luminosities.  


In Figure <ref> we show the impact the DCBH seeding
method has on the scaling relation between the stellar mass of a galaxy
(M_b,*; defined as the stellar mass within the stellar half-mass
radius of the subhalo) and the corresponding M_BH (the mass of the
central black hole within the galaxy) for the original Illustris simulation
(contours) and our new DCBH model (colourmaps). At high redshift we see
excellent agreement between the two, with the exception of some
low-M_BH outliers in the DCBH model.  These outliers are caused when
the black hole which undergoes the fastest growth in the original simulation
is not seeded. Instead there is a lower-mass black hole that grows up toward
the main relation, but will not reach it until a lower redshift.


At lower redshift, we note that the DCBH model only fills the high-mass end of
the scaling relation, without any low-mass objects <cit.>.  This
lack of a low-mass end is due to the low-redshift, low-mass galaxies not satisfying
the conditions for direct-collapse seed formation while those that are seeded
did so at sufficiently high-redshift, providing enough time for both the black
hole and the galaxy to grow to at least a moderate mass scale[We
  stress again this result may also be sensitive to the accretion model.  If a
  less-efficient accretion model were adopted, the slower growth of black
  holes seeded at high-redshift would also lead to a larger population of
  low-mass black holes at low redshift.]. We emphasize that this suggests that
low-redshift, low-mass galaxies rarely host massive black holes formed
via direct collapse; however they may still host SMBHs formed through
other channels, e.g. from NSCs or PopIII stellar remnants.  We also note that
the threshold for seeding in these low-mass galaxies may be sensitive to the
metal enrichment model in Illustris <cit.>. Thus the precise mass scale below which we do not
find DCBHs may be model-dependent, but the qualitative result should hold.





In addition, we note that below z ∼ 2 the high-mass end in the DCBH model
lies below the relation from the original simulation.  As discussed earlier in
this section, black holes in the DCBH model often have lower masses than in
the original simulation (see e.g. the mass function in
Figure <ref>). However, we note that one caveat of our model is
that the feedback from the original simulation remains unchanged, and thus the
growth of massive black holes in the DCBH model may be underestimated, as the
gas supply around black holes which reach the self-regulated regime in the
original simulation has been impacted by a stronger feedback than it would be
with the lower-mass DCBH.  As such, while we do expect the DCBH model to have
a lower high-M_* end, our results provide a lower bound rather than a direct
prediction.



We also consider the global black hole mass density as a function of redshift
for a range of M_BH-bins in the upper panel of Figure <ref>.
As noted above, high-M_BH are likely somewhat under-massive due to
overestimation of the effect of feedback. In this figure we hence compare the
original Illustris mass densities (dashed lines) to a range of possible values
in our DCBH model (shaded region). The shaded areas show the region bounded
by two constraints: the explicit calculation of M_BH from our DCBH
model, and an adjusted value to account for this overestimated feedback.  For the
adjusted mass, we compare the best-fitting relation for the scaling relation
(see Figure <ref>) of the original simulation and the
DCBH model, and at the high-mass end where the DCBH model is lower, we
increase M_BH of the DCBH by the difference between the two
relations.  The lower-M_BH bound is thus expected to be a lower limit
as the black hole growth has been underestimated. Similarly, the
higher-M_BH bound is likely the upper limit, as we expect the DCBHs to have a somewhat lower mass than the original simulation due to the less
frequent seeding.  As such, the shaded regions represent a well-defined
constraint for the mass densities. The width of the shaded regions is clearly
mass dependent, with only the most massive black holes above 10^9  
  M_⊙ being affected significantly by the overestimated feedback in our
DCBH model. 


In the lower panel of Figure <ref> we show the ratio of the DCBH mass
density to that of the original Illustris simulation. Here we see that the
density of low-mass black holes (<10^7  M_⊙) tends to drop off below
z∼3-4, as new seeds rarely form and thus black holes tend to grow out of
the lowest mass bin. For higher mass black holes, we see that at high-redshift
we have much lower mass densities in the DCBH model, but as time passes the
more efficient growth starts to increase the mass density up toward that of
the original Illustris simulation.




From Figure <ref> several interesting features can be deduced. First,
consistent with our findings above, black hole mass density of the DCBH model
is systematically lower with respect to the Illustris prediction. 
Second, by splitting the contribution of total black hole mass density into different mass
bins, we see that at high redshift the density is dominated by low-mass black holes.
Hence, the discrepancy with regard to current high-z observational estimates for black hole
mass density is likely due to underestimating this low-mass population, requiring either a slightly higher DCBH seed efficiency (but see our discussion above that our DCBH seeding model is likely rather conservatively high), or other channels of black hole formation contributing to the total mass density (e.g. SMBHs from PopIII or NSC seeds).  Furthermore, at low redshift the overall density is dominated by high mass black holes, suggesting that, if the growth of SMBHs from low mass seeds is negligible, matching current observational densities would require more efficient growth than our post-processing analysis finds, possibly up to including some super-Eddington phases.  It will be interesting to see if future JWST, Athena, IPTA, and LISA observations will sufficiently inform us about these growth channels to break these degeneracies in theoretical models <cit.>.
 










 §.§ Black hole merger rates




In addition to the impact on the overall black hole population and its cosmic
growth, changing the seeding prescription strongly modifies the rate of black hole mergers, as well
as the mass ratio of the merging black holes.  

Naively one might expect that the number of mergers N_merger scales roughly with f_seed^2, 
if f_seed were the probability of any given black hole having formed.
However, this ignores the full history of the merger trees: more massive black
holes tend to have undergone a large number of mergers (see
Figure <ref> and the discussion in
Section <ref>), and thus have a much higher probability
of at least one of those progenitors having been seeded.  For this reason,
mergers involving more massive black holes will be less sensitive to the
seeding prescription with the number of mergers in the DCBH model becoming
significantly closer to the Illustris prediction for higher mass black holes.




In Figure <ref> we plot the 2-dimensional distribution
of the number density of mergers, binned by both M_1 and M_2 (the masses
of the more massive and less massive black holes involved in the merger,
respectively).  For fixed f_seed models, we find that the merger
frequency peaks at M_1 ∼ 10^6.3 M_⊙ and M_2 ∼ 10^6 M_⊙.  The overall M_1 - M_2 distribution appears roughly comparable (barring an overall normalization shift), but we do note that the lower-mass mergers are affected more strongly by a lower f_seed, so lower f_seed results in a distribution weighted slightly more toward high-mass mergers. 

The DCBH based seeding model, however, produces a qualitatively different
merger-mass distribution which peaks at M_1 ∼ 10^7.75 M_⊙,
M_2 ∼ 10^7 M_⊙ and generally has M_2 within ∼ 1 dex
of M_1.  We also see this quite clearly in the 1D distributions of M_1 and
M_2 in Figure <ref>, which show that mergers in the
DCBH model peak at both higher M_1 and M_2 (see also Figure <ref> and <ref>).

We also note that an earlier analysis of the EAGLE simulation found merger rates dominated by low-mass black holes, but much more strongly peaked toward seed mass black holes.  They also found that decreasing the black hole seed mass to M_BH,seed∼ 10^4  M_⊙ (but maintaining the halo threshold for seed formation) dramatically affected the merger distribution, eliminating nearly all mergers except those with M_2 ∼ M_BH,seed and a small population of M_1, M_2 ∼ 10^8  M_⊙, reinforcing the importance of improved seed models when looking at merger rates <cit.>.


In Figure <ref> we plot 1-dimensional histograms of
M_2 for various bins of M_1. There is some slight M_2-dependence, with
high-M_2 values tending to be less sensitive to f_seed due to the
increase in progenitor black holes. Other than this small effect, the
distribution is qualitatively similar for all f_seed models beyond an
overall normalization shift.  However, as in
Figure <ref>, the DCBH-based seeding model shows a
dramatically different distribution, dropping off rapidly at low M_2.  This
dropoff is consistent with the dropoff in mass function seen in
Figure <ref>, which shows dramatically fewer low-mass black holes
in the DCBH model at low redshifts.  Because the DCBH seeding model has most
black holes form at relatively high redshift, at lower redshift there are far
fewer low-M_BH black holes available to undergo mergers, and thus we
find fewer low-M_2 mergers.  This also means that the low-M_2 dropoff is
partially a redshift-dependent factor, and so high redshift mergers exhibit a
less extreme dropoff (see the left- vs. right- panels of
Figure <ref>).  




In Figure <ref> we show the number density of mergers for
original Illustris and our DCBH model as a function of primary mass (M_1,
the more massive black hole involved in the merger; left panel) and secondary
mass (M_2, the less massive black hole; center panel).  In addition to a significant decrease in merger number, we again see that the DCBH model tends to have much higher secondary black hole masses, while the primary masses tend to be relatively similar (though without reaching as high- or low- masses).  We note that the Illustris simulation predicts that M_1 should peak near ∼ 10^7  M_⊙, broadly consistent with the findings from the EAGLE simulation, though EAGLE produced more mergers with M_2 ∼ M_BH,seed and relatively fewer intermediate-M_2 mergers <cit.>. In addition, we show (right panel) the merger number as a function of chirp mass, a measure of the effective mass of a binary system affecting the gravitational wave signal of the merger, defined as

    M_chirp=(M_1 M_2)^3/5/(M_1 + M_2)^1/5 .

We see that the DCBH model peaks at a higher M_chirp (∼
10^7-10^8  M_⊙, compared to ∼ 10^6-10^7  M_⊙ from the original
Illustris simulation), and does not extend to as low values (ending just below
10^6  M_⊙), suggesting DCBH-DCBH mergers are most common at
slightly higher masses then the original Illustris simulation predicts. It is
worth emphasizing that the merger
density of our DCBH model remains consistent with predicted constraints from
Pulsar Timing Array <cit.>, though the current PTA
constraints remain very broad, hence future tighter constraints would help
differentiating the models.


In Figure <ref> we look at the redshift
distribution of mergers, showing the cumulative number of mergers before
redshift z, binned by M_1 (left panels) and M_2 (right panels), for
f_seed=1 (top), f_seed=0.1 (middle), and the DCBH model
(bottom).  Again we find the f_seed models are qualitatively similar
to one another, but the DCBH seed model is significantly different.  Mergers
with low-M_1 in the DCBH seed model tend to occur exclusively at high
redshifts. Most DCBH mergers with M_1 < 10^7  M_⊙ occur by z
∼ 3.5, whereas in f_seed-based models N_mergers grows
with time throughout the simulation. Similarly, there are few M_2 < 10^6  M_⊙ DCBH mergers below z ∼ 3.5, though 10^6  M_⊙
< M_2 < 10^7  M_⊙ do continue to lower redshifts.  

In Figure <ref> we show the relation between chirp mass and
merger redshift.  In the original Illustris simulation (orange contours), we see that
the typical M_chirp remains roughly constant with time. In these models, black holes are seeded all the way to z = 0, so at all redshifts we find mergers spanning the full range of masses (see also Figure <ref>, which shows comparable increases in N_mergers for all M_BH). In contrast, however, the DCBH mergers (blue contours) show a strong increase in M_chirp with decreasing redshift, as expected from our earlier analysis: at low-redshift, low mass black holes tend to be much rarer (Figure <ref>), resulting in far fewer low-mass mergers (the flattening of the N_mergers for low-M_BH in Figure <ref>).  Thus the redshift distribution of M_chirp detections can provide a strong means of discriminating between seed models.

Finally, in Figure <ref> we consider the rate at which mergers
occur throughout the observable universe, both at a given chirp mass
(dN/dlogMdt, black curves) and the cumulative number above a given chirp mass
(dN/dt, red curves).  Here we see that the merger rate peaks for a chirp
mass of ∼ 3-4 × 10^6  M_⊙ (original Illustris) and ∼
4 × 10^6  M_⊙ (DCBH model).  For the Illustris simulation,
this peak is comparable to the peak seen in Figure <ref>, as we
would expect; although the observable volume is redshift-dependent,
M_chirp is approximately independent of redshift (see
Figure <ref>), so the merger rate is comparable to the merger
density.  However, for the DCBH model, M_chirp evolves strongly with
redshift, shifting the merger rate toward a lower chirp mass. From the
cumulative rate of mergers above a given chirp mass (red curves), we see that
the DCBH model predicts ∼ 1   dex fewer mergers, such that we would only
expect to detect ∼ 0.1 supermassive black hole merger per year.  However, we note that this rate is only for
mergers involving black holes which are well resolved by the simulation, so
detectors sensitive to lower-mass ranges can be expected to make more
detections than our simulation predicts. 



For example, we note that LISA <cit.> has primary sensitivity
which extends down to ∼ 10^4  M_⊙, well below our limit of ∼
10^6  M_⊙.  Given the steepness of the halo mass function (and the
associated increase of low mass galaxies, see Figure <ref>),
we note that the actual number of detectable mergers will likely be higher as
a result of the lower-mass mergers unresolved by our simulation 
This is consistent with previous analyses <cit.>, which have found that the distribution of merger rates peaks at low black hole masses, tending toward the black hole seed mass used in the model.
However, at higher
masses where we resolve well the black holes, we expect our estimates to be much closer to complete.

We also note that this predicted rate assumes that black holes form only via direct collapse,
and the incorporation of black holes formed via alternative seed pathways (e.g. PopIII remnants or in NSCs)
would be expected to significantly increase the overall detection rate.  In
particular, a recent paper by <cit.> used a semi-analytic model to
predict high redshift merger rates, and found that for z > 5 mergers, LISA
detections of mergers between a DCBH and a stellar mass seed would outnumber
the DCBH-DCBH merger detections by two orders of magnitude or
more. Furthermore, similar to <cit.>, <cit.> found
that the merger rate would peak at very low masses, with LISA-detectable
mergers (based on SNR > 7) peaking at M_BH(1+z) ∼ 10^4-10^5
 M_⊙, and <cit.> found that black hole mergers in dwarf galaxies (below the masses considered in this work) will also be detectable by LISA suggesting that the predictions we make for the rate of high-mass DCBH-DCBH mergers represents only a small subset of all LISA detections. 

As such, our predicted merger rate can be treated as a lower-limit for
detections. However, this is dependent on the treatment of merging timescales.
In particular, the black hole model in Illustris repositions black holes
toward the potential minimum of its host halo to prevent N-body scattering
with dark matter or stellar particles, which we note could result in
underestimating the migration time for a secondary black hole brought into a
halo by a galaxy merger.  Furthermore, the simulation merges black holes as
soon as they are within a smoothing length of one another <cit.>, although the actual hardening
timescale may be quite significant <cit.>. 

Adding long hardening timescales would substantially delay the merger of any given black hole pair, decreasing the rate of high-redshift mergers which may be detected.  However, it would be expected to produce a corresponding increase in late-time merger rates, and the lower-redshift at which the merger occurs would make them more easily detectable. Thus we expect that a more realistic hardening timescale would not substantially impact the detection rate of black hole mergers, only that it might lead to a distribution weighted more toward low-redshift (and we further note that an initial investigation suggests that incorporating a delay of ∼1 Gyr would not significantly impact the overall merger rates, consistent with findings of  ).  






§ CONCLUSIONS



In this work, we have developed a novel method to investigate the impact black
hole seeding models have on the cosmic evolution of the entire black hole
population. Once a given seed model is assumed, our method recalculates entire
black hole growth and merging histories taking advantage of the Illustris cosmological simulation of galaxy
formation. Hence our method presents a more self-consistent and comprehensive
approach than the traditional semi-analytic models while at the same time
allowing us to test a large number of seeding models which would be
computationally prohibitively expensive with full hydrodynamical simulations
of structure formation.  Recall the caveat to this approach is that the impact of AGN feedback on the gas supply around the black hole is unchanged during the post-processing approach; thus for black holes with lower mass compared to the original simulation our model may over-estimate the effect of AGN feedback and thus decrease the predicted accretion rate.  However, this effect will generally be significant only for a small fraction of black holes, at the highest mass end where feedback is strongest (for details, see Section <ref> and Figure <ref>).  

In our analysis we test stochastic and physically-motivated seed models based on the direct collapse black hole (DCBH) formation scenario, find the impact on low- and high-redshift black hole populations and merger rates, and determine which observables provide a means for distinguishing between simulated models for seed formation.  From this, we draw the following conclusions:



  * Our novel tool produces post-processed black hole growth histories based on the simulation data available at each timestep. These histories match the original simulation when using the same seeding parameters, but alternate conditions for seeding can also be easily incorporated at a negligible fraction of the original simulation's computational expense.  In particular, we implement simple stochastic seeding models (our f_seed-based models) and a physically motivated seeding model which uses galaxy merger trees to find progenitor galaxies in atomically cooling haloes with low metallicity and low central gas spin in which may be the forming sites of black hole seeds via direct collapse (our DCBH model).


  * Galaxy properties of newly-seeded black holes in the original Illustris simulation are largely redshift independent (though with some increase in metallicity below z ∼ 2), and seeding continues throughout cosmic time.  In contrast, our DCBH seed model is biased toward higher redshifts.  Furthermore, the low redshift seeds (below z ∼ 3.5) tend to be found in galaxies with relatively low baryon fractions, which have a tendency for low-metallicity while still having enough gas to form DCBH seeds.


  * About 10-20% of all simulated galaxies with total mass above 3 × 10^9  M_⊙ have a progenitor which satisfies all our criteria for direct collapse formation, suggesting only ∼ 20% of Illustris black holes may have formed in this manner.


  * All of our seed models (down to f_seed=0.1) produce quasar
  luminosity functions (QLFs) consistent with current observational
  constraints.  Except at very high redshift (where the observed QLF is poorly
  constrained), the bright end of the QLF (L_BH > 10^43.5 erg s^-1) is
  similar for all seed models.  The faint-end QLF shows significant dependence
  on seeding, however, particularly at low redshift, suggesting that upcoming
  surveys such as Athena should be able to put constraints on the seeding
  model. We note, however, that this is based on a single-formation mechanism
  via direct collapse. Additional formation pathways, such as PopIII remnants
  and NSCs, can provide additional black hole populations, particularly at
  lower masses.


  * The scaling relation between black hole and host galaxy masses is
  strongly affected by the seeding model assumed, with the DCBH seed model
  resulting in very few low black hole mass objects at low redshift, in contrast to the original Illustris simulation.  Hence, future observational campaigns that aim to probe the low mass end of the scaling relations may provide much needed insight into the likely nature of black hole seeding.



  * Total black hole mass density is sensitive to the seeding model adopted,
  and our results 
  imply that either DCBH seed need to undergo more efficient growth to match current observational constraints (possibly including a phase of super-Eddington
  accretion) or that the other seed
  formation channels lead to non-neglible cosmic black hole mass
  growth.  Future JWST, Athena, LISA, and PTA observations may be able to shed light on this degeneracy <cit.>.


  * Black hole - black hole mergers are dramatically affected by the seeding
  model.  The total number of mergers scales more strongly than linearly with
  f_seed, with the DCBH model expecting approximately one order of
  magnitude fewer mergers than the original Illustris simulation.
  Furthermore, the expected merger frequency is mass-dependent, with mergers
  involving the lowest masses depending most strongly on the seeding
  frequency. Thus high mass-ratio mergers are most reduced within our DCBH seed model: black holes tend to reach the highest masses only at low-redshift, and the lack of low redshift seeding means there are very few low-mass DCBHs at late times with which the large black hole can merge.


  * The redshift evolution of merger chirp mass is strongly dependent on the
  seed formation model as well.  The continued seeding across all redshifts in
  the original simulation results in a roughly redshift-independent
  distribution of chirp masses.  In contrast, the lack of low-redshift seed
  formation in our DCBH model results in a typical chirp mass that strongly
  increases with time, reaching ∼ 10^8 M_⊙ at low-z, compared to ∼
  10^6.5  M_⊙ in the original Illustris simulation.





In summary, this work demonstrates that seed formation mechanisms can have significant impact on black hole populations across all mass and redshift ranges.
We have shown that differing seed models, including a more physically-motivated mechanism for seeding black holes in cosmological simulations based on the direct collapse scenario, are capable of broadly matching current observations across all redshifts. However, different seeding models have 
distinct predictions for higher redshift and lower mass black holes, which upcoming surveys such as JWST and Athena may be able to constrain. Perhaps most significantly, we have shown that seeding models have a strong effect on black hole merger rates, particularly at low redshift and high mass ratios, such that upcoming gravitational wave detections [e.g. LISA <cit.>, International Pulsar Timing Array <cit.>, European Pulsar Timing Array <cit.>, Parkes Pulsar Timing Array <cit.>, and NANOGrav <cit.>] may be able to differentiate between seed formation pathways for supermassive black holes.




§ ACKNOWLEDGMENTS


We would like to thank Martin Haehnelt, Alberto Sesana, Marta Volonteri and
Priyamvada Natarajan for useful discussions and comments on this work and the referee for a helpful report. CD and DS acknowledge ERC
starting grant 638707 and support from the STFC. This research 
used: The Cambridge Service for Data Driven Discovery
(CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC
HPC Facility (www.dirac.ac.uk). The DiRAC component of
CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. DiRAC is part of the National e-
Infrastructure. Simulations were
run on the Harvard Odyssey and CfA/ITC clusters, the Ranger and Stampede supercomputers at the Texas Advanced Computing Center as part of XSEDE, the Kraken
supercomputer at Oak Ridge National Laboratory as part
of XSEDE, the CURIE supercomputer at CEA/France as
part of PRACE project RA0844, and the SuperMUC computer at the Leibniz Computing Centre, Germany, as part
of project pr85je. 



 mn2e       
        


