Mocking Faint Black Holes during Reionization

To investigate the potential abundance and impact of nuclear black holes (BHs) during reionization, we generate a neural network that estimates their masses and accretion rates by training it on 23 properties of galaxies harbouring them at $z=6$ in the cosmological hydrodynamical simulation Massive-Black II. We then populate all galaxies in the simulation from $z=18$ to $z=5$ with BHs from this network. As the network allows to robustly extrapolate to BH masses below those of the BH seeds, we predict a population of faint BHs with a turnover-free luminosity function, while retaining the bright (and observed) BHs, and together they predict a Universe in which intergalactic hydrogen is $15\%$ ionized at $z=6$ for a clumping factor of 5. Faint BHs may play a stronger role in H reionization without violating any observational constraints. This is expected to have an impact also on pre-heating and -ionization, which is relevant to observations of the 21 cm line from neutral H. We also find that BHs grow more efficiently at higher $z$, but mainly follow a redshift-independent galaxy-BH relation. We provide a power law parametrisation of the hydrogen ionizing emissivity of BHs.


INTRODUCTION
Black holes (BHs) have been prime candidates for the ionization of the Universe (e.g. Setti 1970 andArons &McCray 1970) ever since the early days of the unavailing search for the intergalactic medium (IGM, e.g. Field 1959). With the detection of 22 faint BHs at > 4, Giallongo et al. (2015) revived the question of their role in this process. With the optical depth of the intergalactic free electrons being as low as = 0.054 ± 0.007 (Planck Collaboration et al. 2020), hydrogen reionization is expected to end late enough for stars to be its main driver (e.g. Robertson et al. 2015;Bouwens et al. 2015). Nevertheless, an observational picture where stars alone are responsible for reionization might need to rely on an escape fraction, esc , of ionizing photons from galaxies during the epoch of reionization (EoR) higher than what observed at lower redshifts (see e.g. Naidu et al. 2018, but note also the high individual esc found by e.g. Vanzella et al. 2016or Fletcher et al. 2019, as well as on a population of unobserved faint galaxies. There is thus room for BHs even in the picture of a stellar-dominated EoR. The question remains how large their contribution is. In the most extreme case, Madau & Haardt (2015) found a BH-only reionization scenario under the Giallongo et al. (2015) constraints to positively match the evolution of the volume filling factor of ionized hydrogen, H II . This scenario however fails at re-★ E-mail: eide@mpa-garching.mpg.de producing several other observations. BHs alone would yield IGM temperatures and heating that are too high (see e.g. the comparison of BH-only models to the compilation of IGM temperatures of Garaldi et al. 2019), which is followed by too early adiabatic cooling. Furthermore, in the BHs dominated model of Madau & Haardt (2015), He II reionization would be completed prematurely at ∼ 4.2, shortly after H I reionization, which is at odds with the observed extended He II reionization process (Worseck et al. 2016. The observational constraints on the ionizing ouptut from high-(e.g. Onoue et al. 2017;Parsa et al. 2018;Matsuoka et al. 2018;Kulkarni et al. 2019), as well as theoretical inferences (e.g. Finkelstein et al. 2019), indicate that BHs supply a significant contribution to the ionizing budget, albeit subdominant to that of stars. As discussed by D' Aloisio et al. (2017), BHs can provide an elegant explanation to a flat redshift evolution of the ionizing emissivity, justify the low optical depths in the He II Lyman forest, and importantly, explain the origin for the large variations in the opacity of H I Lyman forest along different sightlines (as investigated by e.g. Chardin et al. 2015).
The interplay between BHs and their host galaxies shapes them both (e.g. Di Matteo et al. 2005). Observations have revealed that massive BHs exist already by = 7.5 (Bañados et al. 2018;Fan et al. 2019), and simulations do not rule this out as unfeasible (e.g. Feng et al. 2015;Di Matteo et al. 2017). The growth of BHs can be captured well by simulations (e.g. Sĳacki et al. 2015;DeGraf et al. 2012;Weinberger et al. 2018;Huang et al. 2018), however, the question of their formation remains still open (e.g.; Regan & Haehnelt 2009;DeGraf et al. 2015a;Inayoshi et al. 2019, for recent reviews). A common numerical approach in large cosmological volume simulations Khandai et al. 2015;Sĳacki et al. 2015;Crain et al. 2015;Weinberger et al. 2018), is to seed galaxies above a mass threshold with a BH of mass close to the mass resolution (typically BH seeds of 10 4−5 M within halos of 10 10−11 M ). This approach leads to a population of BHs at = 0 that matches observations (e.g. Kormendy & Ho 2013, for a recent review), however it does not shed light on the abundance and properties of faint/small mass BHs at higher , a population which can be important during the initial stages of the EoR.
In this work we attempt a novel approach to model a highpopulation of small BHs and study their impact on the EoR. This is done by training a neural network with the properties of the BHs and host galaxies modelled in the cosmological hydrodynamical simulation MassiveBlack-II (MBII, Khandai et al. 2015). The network is then used to mock the BH population (down to halo and BH masses lower than what was assumed and seeded in MBII) at redshifts relevant for the EoR. The paper is structured as follows: in section 2 we introduce the simulations and methods employed to develop the neural network; in section 3 we present our results in terms of BH and galactic properties, as well as the impact on the EoR; in section 4 we discuss some caveats and advantages of our new approach and give our conclusions.

METHODS
In the following we will introduce the simulations and neural network adopted in our work.
MBII was post-processed with the multifrequency ionizing radiative transfer (RT) code CRASH (Ciardi et al. 2001;Maselli et al. 2009;Graziani et al. 2013Graziani et al. , 2018 between = 18 and = 6 (hereafter Eide2018 and Eide2020, respectively; Eide et al. 2018Eide et al. , 2020, to study the impact of various ionizing and heating sources on the physical properties of the IGM during the EoR. To this aim, stars, high and low mass X-ray binaries (XRBs), shock-heated interstellar gas (ISM) and black holes were identified in MBII and assigned spectra depending on their physical characteristics as mass, accretion rates, ages, metallicities or local star formation. The corresponding ionizing emissivities were evaluated and used as input for the radiative transfer calculations.
In this work we make use of the cosmological environment and galactic properties provided by MBII and use them to train a neural network to generate the mass and accretion rate of BHs hosted by such galaxies. In a companion paper we plan, instead, to use numerical simulations as those discussed in Eide2018 and Eide2020 to investigate more in detail the possible impact of the neural network generated BHs on the reionization process of hydrogen and helium.

Cosmological, Galactic and BH Properties
Here we present the 23 galactic and cosmological properties that we use as input to our neural network. From MBII we retrieve the stellar mass * (in 10 10 ℎ −1 M ), the mean stellar metallicity Z, the star formation rate SFR (in M yr −1 ), the mean stellar age (in yr), the dark matter halo mass ℎ (in 10 10 ℎ −1 M ), and the galactic gas mass gas (in 10 10 ℎ −1 M ). We also derive some geometrical and kinematic properties of the galaxies by doing a principal component analysis of the velocities and positions of the gas and stellar particles (see e.g. VanderPlas et al. 2012). We find the galactic gas number density gas (in cm −3 ), the mean velocity of the gas { Additionally, we consider some cosmological properties at the site of each galaxy. Using the cosmic gas number density , we calculate and grid onto 1024 3 regularly spaced cells the overdensity = /¯, where¯is the volume averaged number density. As Di Matteo et al. (2017) found that the tidal field plays a central role in the growth of BHs, we follow their prescription to calculate and grid it. We evaluate the strain tensor in Fourier space,ˆ= 2ˆ/ ( ) from the Fourier transform of the aforementioned gridded overdensity field (following Dalal et al. 2008), and find the tidal field as = − Tr /3. We calculate the eigenvalues of the tidal tensor, and retain the largest one, 1 . As we did for the overdensity, we read off 1 from the grid at the site of the galaxy.
We additionally need to evaluate the accretion rate, luminosity and ionizing emissivity of the BHs. For this, we follow the approach taken in Eide2018 and Eide2020. In line with Shakura & Sunyaev (1973) and the feedback model employed in MBII, we write the bolometric luminosity as = BH 2 (in erg s −1 ), where BH is the BH accretion rate, = 0.1 is an efficiency parameter and is the speed of light. The ionizing emissivity BH (in phots s −1 ) is derived by rescaling the integrated ionizing spectrum with the bolometric luminosity. The spectrum is determined observationally by Krawczyk et al. (2013) and it is essentially a broken power law at hydrogen-ionizing frequencies, with ( ) ∝ and = −1 for ℎ P > 0.2 keV, where is the frequency and ℎ P is the Planck constant. The integral of the ionizing spectrum gives the emissivity. From the rescaled spectrum we also derive the AB luminosity of the BHs, AB,BH (in erg s −1 Hz −1 ).

The Neural Network
We now describe how we construct and train the neural network which ultimately is used to predict the BH masses, BH , and accretion rates, BH . In essence, these are derived from the aforementioned 23 galactic properties, and the network is trained and validated on existing BHs at = 6.
We use TensorFlow 1 with the Keras 2 interface to construct the neural network, while we employ RMSProp as optimiser 3 . Starting from a single input layer (23) with the same number of units as we have learning parameters (23), we consecutively add layers with larger number of units to the network and test its accuracy after adding each new layer. We eventually arrive at a multilayered deep network where introduction of additional hidden layers lead to overfitting and modelling of the noise in the data because of the too many free parameters. We then introduce dropout layers, ( , ), which when enabled ( = 1) randomly remove a fraction of the connections to the preceding layer, helping to increase the versatility of the network and to prevent overfitting (Hinton et al. 2012). The maximally connected network can be described as one that takes an input vector x of our 23 learning (and prediction) parameters and forward feeds it through several hidden layers and before finally reaching an output layer which returns the predictions y NN , where NN 0 = NN BH and NN 1 = NN BH . In its most complex form, it has the following structure, where the arrows indicate that the outputs a −1 of the layer − 1 are used to compute the activation of the units in the next hidden layer , Here, W is the matrix of elements of the connection weights between unit of layer − 1 and unit of layer , ( , 1) = (1 − ) −1 D is a matrix where a fraction of the connections are dropped, b is a bias, and ReLU is the activation function. The final layer has a linear activation function, i.e. y = a = W a −1 + b . In the following we omit ' ' from the notation as we always assume the standard vale = 0.5 (Hinton et al. 2012).
In Fig. 1 we show how the combinations of the hidden layers affect the accuracy 1-MSE (mean-square error) of the predictions, with MSE = −1 =1 ( NN , − , ) 2 for = 0, 1 and where are the validation values. We test the network on the portion of the dataset that it has not been trained on. The figure shows combinations of = 1 . . . 5 hidden layers in addition to the output layer. A network where = 1 only has two ReLU layers, ( ) 3 . We also show combinations of the dropout layers . Networks where all the dropout layers are enabled are labeled as 'ddd' in the figure, whereas a combination such as e.g.
(1)( (0)) 2 3 = 4 or (1) (0) 2 = is labeled 'dnn'. The networks without any dropout layers, labeled 'nnn', usually have the smallest errors, but also the largest potential for overfitting the data. The networks with the smallest error (marked in the figure with a hatch) are 4 for the BH mass, and 5 for the BH accretion rate. The 'd**' networks, where a dropout is applied right after the input layer, generally present larger errors, particularly for the predicted BH accretion rates. Unsurprisingly, the predictions are better with dropout layers for the > 1 layered networks. The 'ndn' class of networks particularly sets itself apart with consistently good 3 For an overview of gradient descent optimization algorithms (including RMSProp) we refer the reader to https://ruder.io/optimizing-gradientdescent/index.html predictions. We find that the 2 network is the one that strikes the best balance between simplicity and predictive power, and is the one we apply in our work.
The training of all the networks was done at = 6, where MBII has a sizable population of 2, 734 BHs. As the distribution of accretion rates is not uniform in our sample of galaxies hosting BHs, we whitened the input data before training. The whitening is done by duplicating galaxies with rare accretion rates, where the properties of the duplicates are added gaussian noise ∼ N (0, 0.05 ) based on the variance of these properties. This extends the training sample and prevents the network from being biased towards only predicting the most common BH masses and accretion rates.
Furthermore, we did not train the networks on galaxies holding BHs with masses equal to those of the seeds, but rather restricted ourselves to BH > 1.1 BH,seed , as the BH properties just after seeding do not immediately reflect the properties of the host galaxy. This left 62% of the data available for training and verification. This also means that any prediction in the range 1 ≤ BH / BH,seed ≤ 1.1 can be used to evaluate the predictive power of the network for masses that it has not been trained for.
We used 7/8 of the full sample of galaxies in the whitened set at = 6 for training the network, before validating it on the remaining 1/8 of the set. By evaluating the MSE of the predictions of the network versus the validation data, we conclude that our ability to predict the BH masses and accretion rates with our chosen

Relation between Galactic and Black Holes Properties
We now turn to examine if any of the 23 galactic properties plays a dominant role in predicting the BH masses and accretion rates. We do this by generating (i) a network which is the same as the original one except that now one parameter is removed, and (ii) a network using solely this parameter. In both cases we estimate the MSE on the predicted BH masses and accretion rates.
In Fig. 2 we show the networks' accuracy, 1-MSE, in predicting BH masses and accretion rates. A low 1 − MSE for the models plotted in red reflects a poorer performance of the network without the component under consideration. Conversely, a high 1 − MSE for the models in blue means that the predictive power of this singleparameter network is better.
We first note that the multi-parameter networks have a higher accuracy (> 88%) compared to the single-parameter networks (< 88%), and are hence performing better. The best single parameter for predicting the BH accretion rate and mass is * , for which the networks recover BH and BH with an accuracy of 88% and 75%, respectively. As for the multi-parameter networks, they perform worst when removing * , yielding an accuracy of 99% and 89% for BH and BH , respectively.
While the relevance of * is clear both for single-and multiparameter networks, this is not the case for the other parameters. The three next-most important parameters for the determination of BH are the mean stellar age , the stellar ionizing emissivity and the component of the mean velocity of the stars * for the multi-parameter network. For the single-parameter network, instead, these are the dark matter halo mass ℎ , the galactic gas mass gas and the stellar AB luminosity AB . Similarly, for BH the most relevant quantities in the multi-parameter network are the component of the mean velocity of the stars * , the component of the stellar velocity dispersion * and the component of the mean velocity of the gas gas , indicating that the network captures the dependency between accretion and environmental kinematics; while AB , SFR and gas yield the highest accuracies in the singleparameter networks. We recover the same order of importance for the single-parameter networks by calculating the correlation coefficient between BH or BH and the parameter in question. Again, it should be noted that the predictions of these single-parameter networks are far less accurate (∼ 60%) than the multi-parameter networks (> 89%). As a further test, we create a network with only * , AB , gas and ℎ as input parameters. It recovers BH and BH with an accuracy of 97% and 85%, respectively. This highlights our need for the full network's complexity if the goal is to recover the accretion rates as precisely as possible.

Populating Galaxies with Black Holes
Our network was generated from = 6 BHs and their host galaxies. We now turn to examine how it performs at other redshifts, for galaxies it has not been trained for. This can also reveal any redshift evolution in BH properties.
To do so, we use the network to seed BH-hosting galaxies at various with BHs with predicted masses NN BH and accretion rates NN BH , and compare them to the BH and BH directly obtained from MBII. We show the deviation between the generated and true values at various redshifts in Fig. 3. In the BH,seed < BH < 10 6.5 ℎ −1 M mass range, the deviations vary from 20% larger to 30% smaller, with the largest ones at BH > 10 7 ℎ −1 M , where we have a poorer statistic of the training set. At BH < 1.1 BH,seed where the network has not been trained for (this mass limit is indicated by a vertical dashed line in the figure), the predictions are 2-4% larger than the true values at = 6, while at = 9 they are ∼ 15% lower. This indicates that our network is very powerful in predicting masses it has not been trained for. The predicted accretion rates deviates from being between ∼ 30% larger to ∼ 60% smaller for 10 5 < BH /(ℎ −1 M /(0.98 Gyr)) < 10 7.5 . Also in this case, the predictions are best within the most common range of accretion rates. At the high mass and accretion rate end, the network underpredicts the true values at > 6 and overpredicts them at < 6. This indicates that the BH formation efficiency declines with decreasing . We also see this effect within the central mass and accretion ranges, albeit in a much more moderate fashion-e.g. at = 9 ( = 5), BH = 10 6 ℎ −1 M BHs are on average predicted to be 16% less (4% more) massive, and the accretion rate of BH = 10 6 ℎ −1 M /(0.98 Gyr) BHs is predicted to be ∼ 50% (∼ 4%) lower-indicating that this is not merely an effect caused by lacking statistics of our training set.
Next, we populate all galaxies in the range = 5 − 18 with a BH using the neural network including all the 23 physical properties described in sections 2.2 and 2.3. We thereby create a much larger population of BHs than is present in MBII. In Fig. 4 we show the resulting mass function at various redshifts. While at = 18 the BH population is limited to the range 10 4 < BH /(ℎ −1 M ) < 10 5 , the peak of the mass function shifts towards higher values with decreasing redshift, and by ∼ 6 we have BHs with masses as high as all times. This is more than a magnitude lower than the seed mass BH,seed of MBII, and reflects that the predictions of the network are not restricted by the mass range it was trained on. Note that the generated BH mass function is not dissimilar to those for a range of physical BH seed models at = 15 − 18 (e.g.; Volonteri et al. 2008). Our generated BH population appears to exploit reasonably well the actual resolution of the simulation, introducing BHs at smaller masses and earlier time when they are indeed expected to form. We note here that seeding halos of mass smaller than the one used in the MBII prescription is not a mere extrapolation, but is made possible by the fact that, even if the mass falls outside of the range used for the training, all the other 22 properties are not restricted by any limit. Hence the robustness of our procedure.
In Fig. 5 we show the UV luminosity function (LF) of the BHs at = 6, and compare it to the LF of the MBII BHs, as well as to the observationally determined LFs of Giallongo et al. (2015) and Kulkarni et al. (2019). As our network slightly underpredicts the highest accretion rates, we have a small deficit of bright BHs compared to both the MBII-seeded BHs and the Giallongo et al. (2015) observations. It should be noted, though, that the bright end of the observed LF may be overestimated (Parsa et al. 2018), so that our conservative result might be more realistic. This is further corroborated by the recent compilation of Kulkarni et al. (2019), based on 66 QSOs at 5.5 < < 6.5, as our predicted LF matches their observations at all AB . The agreement of our LF with the original one from the MBII and the Giallongo et al. (2015) LF is extremely good in the range −17 < AB < −15. Our network also predicts a substantial population of faint BHs which are not present in MBII, and yields a LF with a knee at AB = −15 and no turnover at least down to AB = −5.
In Fig. 6 we plot the comoving volume averaged emissivity,¯, in comparison to values inferred from observations. The predicted emissivity increases exponentially from = 18, when¯= 7.6×10 41 phots s −1 cMpc −1 , to = 5, where¯= 1.2×10 52 phots s −1 cMpc −1 . This evolution can be parametrized as a power law, log¯( ) = −0.5097 + 53.86, using a least-square fit to the predictions. We find that the predicted emissivity is much higher than that inferred by Mason

Impact on the Reionization Process
The final question we address in our study is whether such a population of faint BHs could have a significant impact on the EoR. While we plan to run simulations as those presented in Eide2018 and Eide2020 including these faint BHs, here we limit the analysis to a simpler approach. We calculate the filling factor HII of ionized hydrogen (H II) as (Madau et al. 1999): where esc is the escape fraction of ionizing photons,¯is the volume averaged ionizing emissivity,¯H is the average cosmic hydrogen number density and¯r ec = (¯H ( )) −1 is the recombination time, for which we assume a clumping factor = 1, 5, 10 and a case-A recombination coefficient at = 10 4 K. We calculate HII for the MBII BHs, as well as for those seeded by our neural network, assuming esc = 1 for both. As a comparison, we also calculate HII for the stars of MBII, assuming esc = 0.15 as in Eide2018 and Eide2020. These are shown in Fig. 7. For = 5 we find that the population of mainly faint BHs seeded with our neural network results in a reionization history in which the BHs have a central, albeit not dominant, role, reaching HII > 0.15 (0.5) at = 6 (5). This is in stark contrast to the massive BHs of MBII, which reside only in the most massive galaxies and yield HII < 0.05 (0.2) at the same redshifts. As expected, the stars dominate the reionization process, producing HII ∼ 1 already at = 6. Finally, we should note again that the contribution from the network generated BHs should be regarded as an upper limit, as not every galaxy is in reality expected to host an active BH. We defer to future work a refinement of this approach.

DISCUSSION AND CONCLUSIONS
In the cosmological hydrodynamical simulation MassiveBlack-II (Khandai et al. 2015), galaxies with a halo mass in excess of ℎ,seed = 5 × 10 10 ℎ −1 M are populated with seed black holes (BHs) with BH,seed = 5.5 × 10 5 ℎ −1 M . While this prescription assures that the BH population has physical properties consistent with observations at 6, a different seeding procedure, with BHs hosted also in smaller galaxies, might have a strong impact, among others, on the role played by BHs in the reionization process of the intergalactic medium and the related 21 cm signal. To investigate this in more detail, we have trained a neural network using the properties of galaxies harboring BHs at = 6. This network allowed us to mock BHs in all galaxies down to the resolution limit of the simulations at all redshifts, corresponding to halos of mass ∼ 9 × 10 6 ℎ −1 M . By design and through training, the network replicates the properties of the pre-existing BHs in the simulation.
Our network predicts the BH masses and accretion rates of existing BHs with great precision (> 99% and > 93%, respectively). Interestingly, we mock BHs with masses below the MBII seed mass when applying the network to all galaxies, also those with halo masses below ℎ,seed . Although the seeding procedure is extrapolated to lower masses, our predictions of BH and BH are robust because they are constrained in 23 dimensions with a high accuracy (e.g. the predictions of the mass function in Fig. 4 where BHs are lighter at higher ). In fact, a galaxy with ℎ < ℎ,seed may still share up to 22 other parameters with galaxies hosting BHs in MBII, and thus be tightly constrained in these other dimensions. Additionally, as the networks have been trained on galaxies with BH > 1.1 BH,seed , the predictions in the range 1 ≤ BH ≤ BH,seed have been used to confirm the strong predic-tive power of the networks for masses that they had not been trained for.
We find that removal of one parameter, including ℎ , from our network did not lead to a significant deterioration of its predictions. Similarly, not a single one of the input parameters provides predictions as accurate as the full network. The exercise of removing parameters from the network, nevertheless, highlighted that the stellar mass of the galaxy, * , is the most important parameter. Alone, it can predict the BH mass with an accuracy of 0.88, while ℎ has an accuracy of ∼ 0.80. It is harder to infer the effect of the velocity dispersion. From the well-known BH - * relation (Ferrarese & Merritt 2000) we expect the velocity dispersion to be a dominant parameter, but we cannot directly infer its role as it is not a single input to our network, but it is rather decomposed along each coordinate axis , and .
Even though the formation efficiency of BHs is declining with decreasing redshift, this does not necessarily imply that their growth is decoupled from the stellar growth. On the contrary, our network has a strong dependence on stellar properties, such as stellar mass and age. Observations indicate that the SFR history is closely related to the BH accretion history (see e.g. the review by Madau & Dickinson 2014), but they are not identical. Furthermore, our powerlaw parametrisation of the BH emissivity with a slope of −0.5 is similar to that of −0.45 which has been found for the stellar UV density at > 9 Madau 2018), although our BH emissivity relation lacks the turnover to a slower growth at ≤ 9, which is seen instead for the stellar UV density. The strong dependence on the tidal field, overdensity, gas mass and halo mass also indicates that an environment that promotes stellar growth also positively influences BH growth. Such highly biased regions are in fact required to avoid quenching of the growth of the lightest BHs by SN-feedback (Inayoshi et al. 2019).
While we took great care in the training of the network, its performance is still somewhat limited by the size of the training sample, both in terms of number of objects and range of masses covered. Furthermore, our network was trained on = 6 galaxies hosting a BH, as only at that time does MBII produce a sizable population of BHs. This situation would improve by adopting larger and/or higher resolution simulations, such as B T  or Illustris TNG300 (Nelson et al. 2018), or employ simulations specifically designed for this task. Nevertheless, we found that the network predicts the properties of the majority of BHs at all redshifts with high accuracy, indicating also that there is no significant evolution in the relation between the environment and the BHs' properties, in line with Huang et al. (2018). However, our slight deficiency of brighter BHs at > 6 (and surplus at < 6) points to these being formed more efficiently at early times (see also e.g. DeGraf et al. 2012DeGraf et al. , 2015b. Our results suggest that a galaxy at > 6 with properties identical to those of one at = 6 is more likely to host a brighter BH. We also note that recent work which relaxes repositioning of the BHs (as done instead in MBII and most large scale cosmological simulations such as the previously mentioned Illustris) and uses additional dynamical friction (e.g. Tremmel et al. 2018;Pfister et al. 2019;Barausse et al. 2020) should provide more realistic predictions for the early BH populations and their BH merger rates, possibly leading to lower occupation fractions and central BH masses in the galaxies. In the future, different scale simulations (such as those mentioned above) could be used as additional training sets.
Our slight deficiency of the brightest BHs at lower in turn ensures a perfect match at = 6 to the recent LF of Kulkarni et al. (2019), and a perfect match at −15 > AB > −17 to the LF of Giallongo et al. (2015). The most interesting feature of our results is however the large population of faint, AB > −15, BHs. Such a population is entirely possible, as the pre-existing BHs (and the combined contributions from other energetic X-ray emitting sources) in MBII are unable to account for more than a few per cent of the unresolved X-ray background , leaving ample margin for a higher contribution at high redshift.
This predicted population of BHs is unable to drive EoR alone, but it may play an important role nevertheless. Our mocked BHs do not yield enough ionizing photons to fulfil the constraints on the ionizing budget calculated from observational constraints by Mason et al. (2019). However, our emissivites are an order of magnitude larger than those inferred from integrating the LF of the brighter QSOs of Kulkarni et al. (2019). Our BHs leave a significant imprint on the H II volume filling factor, which at = 5 ranges from = 0.41 with a clumping factor = 10 to = 0.83 with = 1. The existing BHs in MBII can at best yield = 0.23 with = 1, but while this population satisfy the bright end of the LF down to > 2, it does not include the fainter population that our network predicts. Our population of mocked BHs is neither negligible, nor is it as dominating as the one of Madau & Haardt (2015). Further work is needed to investigate whether they will induce an extended He II reionization epoch as observations imply (Worseck et al. 2016 without providing undue heating (see e.g. D' Aloisio et al. 2017;Garaldi et al. 2019). We plan to investigate this more in detail with numerical simulations following the work of Eide2018 and Eide2020.
A more prominent population of high-, small mass BHs could also have an important impact on the 21 cm signal from neutral hydrogen in the IGM, by partially ionizing and heating the gas prior to full reionization (e.g. Madau et al. 1997).
Our conclusions can be summarized as follows.
• We train a neural network on properties of BH hosting galaxies at = 6. For our training sample, this predicts the mass, BH , and accretion rate, BH of BHs with an accuracy > 99% and > 93%, respectively. These properties at other redshifts are also predicted with high precision.
• BH and BH are predicted with the most relevant single parameter, the stellar mass * , with an accuracy of 88% and 75%, respectively. Removing * degrades the network to accuracies of 98.6% and 88.8%. The predictions of our network are robust, even when single parameters are ill-defined.
• The neural network is slightly less effective at predicting the brightest and most massive BHs at > 6, and conversely predicts a population of slightly brighter BHs at < 6. This points to a decrease in BH formation efficiency with decreasing .
• Populating all galaxies with a nuclear BH, we predict a substantial population with mass below that of the seeds at all redshifts. This results in a LF at = 6 with a knee at AB = −15 and a lack of turnover at least down to AB = −5.
• Our predicted population of BHs can contribute significantly to H reionization, yielding a Universe in which H is ∼ 15% ionized by BHs at = 6 for a clumping factor of 5. The bright BHs alone, which are well reproduced by MBII, predict instead a Universe that is only ∼ 5% ionized at the same redshift. knowledges funding from NSF ACI-1614853, NSF AST-1616168, NASA ATP 19-ATP19-0084 and NASA ATP 80NSSC20K0519, ATP 80NSSC18K101.

DATA AVAILABILITY
No new data were generated or analysed in support of this research.