Forming massive seed black holes in high-redshift quasar host progenitors
MMNRAS , 1–15 (2021) Preprint 11 February 2021 Compiled using MNRAS L A TEX style file v3.0
Forming massive seed black holes in high-redshift quasar hostprogenitors
Alessandro Lupi, ★ Zoltán Haiman, and Marta Volonteri Dipartimento di Fisica “G. Occhialini”, Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy Department of Astronomy, Columbia University, 550 West 120th Street, New York, NY, 10027, U.S.A. Sorbonne Universitès, UPMC Univ Paris 6 et CNRS, UMR 7095, Institut d’Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The presence of massive black holes (BHs) with masses of order 10 M (cid:12) , powering brightquasars when the Universe was less than 1 Gyr old, poses strong constraints on their formationmechanism. Several scenarios have been proposed to date to explain massive BH formation,from the low-mass seed BH remnants of the first generation of stars to the massive seed BHsresulting from the rapid collapse of massive gas clouds. However, the plausibility of some ofthese scenarios to occur within the progenitors of high-z quasars has not yet been thoroughlyexplored. In this work, we investigate, by combining dark-matter only N-body simulationswith a semi-analytic framework, whether the conditions for the formation of massive seedBHs from synchronised atomic-cooling halo pairs and/or dynamically-heated mini-haloes arefulfilled in the overdense regions where the progenitors of a typical high-redshift quasar hostform and evolve. Our analysis shows that the peculiar conditions in such regions, i.e. stronghalo clustering and high star formation rates, are crucial to produce a non-negligible number ofmassive seed BH host candidates: we find ≈ 𝑧 =
6. This demonstrates that the progenitors of high-redshift quasar host haloescan harbour early massive seed BHs. Our results further suggest that multiple massive seedBHs may form in or near the quasar host’s progenitors, potentially merging at lower redshiftsand yielding gravitational wave events.
Key words: stars: black holes – quasars: supermassive black holes – methods: numerical –cosmology: first stars
The observations of quasars at redshift 𝑧 > M (cid:12) (Fan et al. 2006; Mortlocket al. 2011; Bañados et al. 2018; Decarli et al. 2018) pose tight con-straints on the formation and growth of these objects, and representa challenge for theoretical models (see, e.g. Inayoshi et al. 2020, fora recent review). Models of the evolution of the quasar luminosityfunction over the range 0 < 𝑧 ∼ < seeds of MBHspowering these quasars were in place early on, and subsequentlygrew via accretion and mergers (e.g. Soltan 1982; Small & Bland-ford 1992; Kauffmann & Haehnelt 2000; Shankar et al. 2009; Shenet al. 2020). To date, several scenarios have been proposed to explainthe early formation of MBHs, either relying on established physi-cal processes or in some cases based on unconventional physics. Inthe so-called ‘light-seed’ scenario, BHs formed already at 𝑧 (cid:38) ★ E-mail: [email protected] (cid:12) (Abel et al. 2002; McKee & Tan 2008). More recent studiesfound that fragmentation in the stellar accretion disc and radiativefeedback reduce the typical masses of these stars, but with theirIMF still top-heavy, producing remnant seed BHs in the range 10–1000M (cid:12) (Stacy et al. 2012; Hirano et al. 2015; Hirano & Bromm2017; Hirano et al. 2018; Wollenberg et al. 2020; Kimura et al.2020). Although these seeds can form at very early times, sustained © a r X i v : . [ a s t r o - ph . GA ] F e b A. Lupi et al. accretion at the Eddington limit is required to reach the MBHmasses observed in high-redshift quasars. This condition is unlikelyto be realistic, especially at early times, when galaxies are small andsupernova feedback and/or stellar UV radiation can easily expel gasfrom the halo hosting the BH.A possible solution to this limitation is the occurrence of rela-tively short super-Eddington phases. Exceeding the Eddington rateeven by a relatively modest factor, as expected in so-called slim diskmodels, could allow initially small BHs to efficiently grow by ordersof magnitude in mass in a few Myr (Madau et al. 2014; Volonteriet al. 2015; Lupi et al. 2016; Pezzulli et al. 2016). Given a suffi-ciently large inflow rate from larger radii, the Eddington rate canbe exceeded by a much larger factor (several orders of magnitude),due to the trapping of the radiation in the inflowing plasma (Begel-man 1979), resulting in rapid growth via "hyper-Eddington" ac-cretion (Volonteri & Rees 2005; Tanaka & Haiman 2009; Pacucciet al. 2015; Inayoshi et al. 2016; Sakurai et al. 2016). Caveats tothis scenario include mechanical feedback from a jet, which cansignificantly affect the dynamics of the accreting gas (Regan et al.2019), and the need for the accreted gas to efficiently shed angularmomentum (which may be facilitated in the early stages of growthby the presence of a nuclear star cluster; Alexander & Natarajan2014).Another possible solution is the formation of intermediate-mass BHs with masses of about 10 M (cid:12) in mildly metal-enrichednuclear stellar clusters, either via collisions between stars (PortegiesZwart & McMillan 2002; Gürkan et al. 2004; Omukai et al. 2008;Devecchi & Volonteri 2009; Devecchi et al. 2010, 2012; Katz et al.2015; Boekholt et al. 2018; Reinoso et al. 2018; Das et al. 2020;Tagawa et al. 2020) or the runaway merger of stellar-mass BHs(Davies et al. 2011; Miller & Davies 2012; Lupi et al. 2014; Tagawaet al. 2015). This range of masses is also obtained in models invokingphysics beyond the standard model, for instance from the collapseof self-interacting (Balberg et al. 2002; Pollack et al. 2015) or otherforms of dissipative dark matter (D’Amico et al. 2018; Latif et al.2019). However, because of the relatively low initial mass of theseseeds, it is not clear yet whether they can represent the seeds of theMBHs in high-redshift quasars.Finally, in the ‘heavy-seed’ scenario, even more massive seedscan form at lower redshift ( 𝑧 (cid:46) formation and cooling(Omukai 2001; Dijkstra et al. 2008; Agarwal et al. 2012; Latif et al.2013; Visbal et al. 2014b; Latif et al. 2015), and likely proceedsthrough the intermediate stage of a supermassive star (Hosokawaet al. 2012, 2013; Woods et al. 2017; Haemmerlé et al. 2018), ora quasi-star (Begelman et al. 2008; Begelman 2010; Dotan et al.2011; Choi et al. 2013; Fiacconi & Rossi 2017). In these models,which have been collectively dubbed to produce a "direct collapseblack hole (DCBH)", a strong UV radiation field is thought to berequired to efficiently suppress H formation in the haloes, the maincooling channel in cold pristine gas, able to trigger fragmentationand inhibit the monolithic collapse of 10 − M (cid:12) .The large required UV intensity ( 𝐽 ∼ > − erg s − cm − Hz − ; e.g. Wolcott-Green & Haiman 2019) istypically realised in the case of synchronised pairs, i.e. pairs of The Eddington limit here refers to the accretion rate (cid:164) 𝑀 Edd ≡ 𝐿 Edd / 𝑐 where 𝐿 Edd is the Eddington luminosity and 𝑐 the speed of light. atomic cooling haloes (ACHs; separated by ∼ < formation (Dijkstra et al. 2008; Agarwalet al. 2012; Dijkstra et al. 2014; Latif et al. 2014; Visbal et al. 2014b;Chon et al. 2016; Regan et al. 2017; Chon et al. 2018).Another requirement in this scenario is that the gas avoidscooling and star-formation prior to reaching the ACH stage. H -freeACHs may be produced by a combination of a lower UV backgroundflux ( 𝐽 ∼ . −
1) and the unusually rapid assembly of a subsetof ACHs: the corresponding reduction in H cooling (Fernandezet al. 2014) and the accompanying increase in dynamical heating(Yoshida et al. 2003; Inayoshi et al. 2015, 2018) of gas in the mini-haloes (haloes with virial temperatures below the atomic coolingthreshold, i.e. 𝑇 vir (cid:46) K) could naturally result in star-free ACHs.While these halos would form H and cool once they reach the ACHthreshold, it has recently been suggested that they may neverthelesshave large central gas inflow rates and form relatively massive seedBHs (Wise et al. 2019; Sakurai et al. 2020) of at least 10 M (cid:12) (Reganet al. 2020b).In this work, we investigate the plausibility of the"synchronised-pair" and "dynamical-heating" scenarios (both sep-arately and in combination) for the formation of massive seed BHsas progenitors of 𝑧 > ∼ M (cid:12) )dark matter halo by 𝑧 ≈
6, i.e. typical of haloes hosting brightquasars at this redshift. This is achieved by combining semi-analyticprescriptions with N-body simulations resolving the merger historyof such a massive halo.The rest of this paper is organised as follows: in § 2 we sum-marise the underlying N-body simulations and present our improvedsemi-analytic framework, in § 3 we present our main results, in § 4we discuss the caveats of our work, and in § 5 and § 6 we discussand summarise our conclusions and the implications of this work.
In this work, we investigate whether the so-called "direct-collapse"scenario can explain the formation of MBHs in high-redshiftquasars. To this aim, we build a semi-analytic model on top ofthe high-resolution dark matter (DM) only equivalent of the hy-drodynamic simulation presented in Lupi et al. (2019), a zoom-insimulation of a 100 Mpc comoving box targeting a high-sigmapeak, i.e. a massive halo with 𝑀 vir ∼ × M (cid:12) at 𝑧 =
6. Thesimulation is set to guarantee a mass resolution of 10 M (cid:12) within2.5 𝑅 vir of the halo at 𝑧 =
6, which corresponds to an initial vol-ume of ∼ 𝑧 = Δ 𝑡 = Λ CDM Universe with the Planck Collaboration et al. (2016) cos-mology parameters, i.e. Ω m = . Ω Λ = . Ω b = . 𝐻 = .
74 km s − Mpc − , 𝜎 = . 𝑛 s = . 𝑚 DM ≈ × M (cid:12) , and the spatial resolu-tion (Plummer-equivalent gravitational softening) is set to 640 pccomoving at 𝑧 ≥
15 and switched to 40 proper pc at 6 ≤ 𝑧 < MNRAS , 1–15 (2021) assive seed black holes in quasar hosts trees (also checking the consistency of the identified haloes through-out cosmic history using consistent trees (Behroozi et al. 2013b).Then, we start from the highest redshift available, and follow thecosmic evolution of all haloes identified. This approach, which mir-rors that employed in semi-analytic models, allows us to consistentlyfollow the evolution of all haloes, and check at every step whetherthey fulfill a defined set of criteria for direct collapse that we nextdescribe. During the cosmic history, haloes grow via accretion from the cos-mic web and merger of smaller substructures, with baryons fol-lowing the dark matter and accumulating within the haloes. As thehalo grows in mass with redshift 𝑧 , gas heats up to the halo virialtemperature 𝑇 vir = . × (cid:16) 𝜇 . (cid:17) (cid:18) Ω m Ω 𝑧 m Δ vir 𝜋 (cid:19) / (cid:18) ℎ𝑀 vir M (cid:12) (cid:19) / + 𝑧 , (1)where 𝜇 = .
22 is the mean molecular weight for a fully neutralmedium, Ω 𝑧 m = Ω m ( + 𝑧 ) /[ Ω m ( + 𝑧 ) + Ω Λ ] , Δ vir = 𝜋 + 𝑑 − 𝑑 , 𝑑 = Ω 𝑧 m −
1, and 𝑀 vir is the halo virial mass in solarunits (Barkana & Loeb 2001).Given the primordial composition, baryons can cool only viahydrogen and helium atomic cooling (above 𝑇 ∼ K) and viamolecular hydrogen cooling at lower temperatures. If we define theH cooling time 𝜏 H = 𝛾 − 𝑛 gas 𝑘 B 𝑇 Λ H 𝑛 𝑓 H , (2)where 𝑛 gas = 𝜌 /( 𝜇𝑚 H ) , 𝜌 is the gas density, 𝑘 B is the Boltzmannconstant, 𝑚 H is the hydrogen mass, 𝑛 H = . 𝜌 / 𝑚 H is the hydrogennuclei density, 𝛾 = / 𝑓 H = 𝑛 H / 𝑛 H is themolecular hydrogen number fraction, and Λ H is the H cooling rate,which can be approximated as 10 − . ( 𝑇 vir / ) . erg s − cm for120 K ≤ 𝑇 ≤ 𝜏 H ≈ 𝐻 Ω / ( + 𝑧 ) − / , (3)we obtain an estimate for the minimum H fraction required for effi-cient cooling. If we assume the presence of a non-negligible Lyman-Werner background 𝐹 LW = 𝜋𝐽 − erg s − cm − Hz − , we canwrite the equilibrium H fraction able to form in the halo as 𝑓 H = 𝑘 H − 𝑛 e 𝑘 diss , (4)where 𝑘 H − = . × − 𝑇 . cm s − is the coefficient for thereaction H+e − → H − + 𝛾 , assumed to be the dominant process andthe bottleneck in H formation, 𝑛 𝑒 = 𝑥 𝑒 𝑛 H is the electron density,and 𝑘 diss = 𝐹 LW / × − s − = J / . × s − . For 𝑥 e , weassume the relic abundance 𝑥 e = . × − √ Ω m /( Ω b ℎ ) (Peebles1993).By equating this equilibrium fraction to the minimum requiredfor cooling, and assuming 𝑇 = 𝑇 vir , 𝜌 = 𝜂 gas 𝜌 max , with 𝜌 max ≡ 𝜇𝑚 H Ω b ℎ ( 𝑇 vir / ) . the maximum gas density reachableby adiabatic contraction (Fernandez et al. 2014) and 0 ≤ 𝜂 gas ≤ 𝑀 𝑏 vir , min ( M (cid:12) ) ≈ . 𝛼 . 𝜉 . 𝑛 − . , ( + 𝑧 ) . , (5) where 𝛼 vir = 𝑇 vir /( 𝑀 vir / M (cid:12) ) / , 𝜉 = 𝐽 ( ℎ √ Ω m )/( 𝜇𝑥 e ) , and 𝑛 H , = . 𝜂 gas ( 𝜌 max / 𝑚 H )( 𝑇 vir / ) − . .Recently, Schauer et al. (2019) and Kulkarni et al. (2020) in-ferred the minimum mass for H cooling from hydrodynamic sim-ulations (performed with different codes, setups and sub-grid pre-scriptions), finding significant discrepancies. Hence, it is importantto check how our Eq. (5), based on simple scaling relations and anapproximate form of the H cooling function, compares to thesenumerically derived results. In Fig. 1, we show a comparison of theminimum masses over the redshift range 𝑧 = −
30 for differentvalues of 𝐽 , i.e. 0.01 (solid lines), 1 (dashed lines), and 10 (dot-ted lines). The results by Kulkarni et al. (2020) are shown in blue,our Eq. (5) with 𝜂 gas = . 𝐽 =
0) as a green dot-dashedline. From the figure, it is clear that our results exhibit a decreasingtrend with redshift similar to Kulkarni et al. (2020), although with ahigher normalisation resulting from 𝜂 gas = .
2. Indeed, if 𝜂 gas = ∼
2, getting closer toKulkarni et al. (2020). Nevertheless, the difference always remainsto within a factor of ∼ three over most of the relevant redshift range.On the other hand, the results by (Schauer et al. 2019) always remainin between, but are higher for the case of no impinging LW flux.To conclude, our Eq. (5) represents a reasonable estimate, inrough agreement with state-of-the-art numerical results over the rel-evant ranges of redshift and 𝐽 . Nonetheless, the minimum masseswe find are typically a factor of ∼ − 𝑧 ∼ <
16 and above 𝐽 ∼ > 𝑧 ∼ 𝑧 =
15, it resulted in a few percent higher number of can-didates relative to our fiducial case. However, we caution that ourability to fully test this caveat is limited by the resolution of oursimulations (with a particle mass of 𝑀 dm ≈ × M (cid:12) , we areunable to resolve mini-haloes with masses much below 10 M (cid:12) ).Above this mass threshold, a plausible process that is able tocounteract radiative cooling is the compressional heating providedby the mass accretion onto the halo. This has been dubbed "dynam-ical heating" in the literature, and the corresponding heating rate Γ dyn can be approximated as 𝑛 Γ dyn = 𝑛 H . 𝜇 𝑘 B 𝛾 − (cid:18) 𝛼 vir 𝑀 − / d 𝑀 vir d 𝑡 + 𝑇 vir + 𝑧 d 𝑧 d 𝑡 (cid:19) , (6)where the first term on the right-hand side represents compressionalheating accompanying the mass growth of the halo (Yoshida et al.2003), and the second term accommodates the decrease of the char-acteristic densities and the virial temperature (at fixed halo mass)towards lower redshift, due to the expansion of the Universe. Numer-ically, Γ dyn is computed by finite difference between halo massesfrom consecutive snapshots.In order to prevent H cooling, we need 𝑛 Γ dyn ≥ 𝑓 H 𝑛 Λ H .Since, for short enough time-scales (order of Myr at the high red-shifts considered in this work), the second term on the right-handside of Eq. (6) is small, it can be neglected, translating the conditioninto a critical mass accretion rate (cid:18) d 𝑀 vir d 𝑡 (cid:19) crit = . 𝜇 𝑓 H 𝑛 H ( 𝛾 − ) 𝑘 B 𝛼 − Λ H 𝑀 / . (7) MNRAS000
15, it resulted in a few percent higher number of can-didates relative to our fiducial case. However, we caution that ourability to fully test this caveat is limited by the resolution of oursimulations (with a particle mass of 𝑀 dm ≈ × M (cid:12) , we areunable to resolve mini-haloes with masses much below 10 M (cid:12) ).Above this mass threshold, a plausible process that is able tocounteract radiative cooling is the compressional heating providedby the mass accretion onto the halo. This has been dubbed "dynam-ical heating" in the literature, and the corresponding heating rate Γ dyn can be approximated as 𝑛 Γ dyn = 𝑛 H . 𝜇 𝑘 B 𝛾 − (cid:18) 𝛼 vir 𝑀 − / d 𝑀 vir d 𝑡 + 𝑇 vir + 𝑧 d 𝑧 d 𝑡 (cid:19) , (6)where the first term on the right-hand side represents compressionalheating accompanying the mass growth of the halo (Yoshida et al.2003), and the second term accommodates the decrease of the char-acteristic densities and the virial temperature (at fixed halo mass)towards lower redshift, due to the expansion of the Universe. Numer-ically, Γ dyn is computed by finite difference between halo massesfrom consecutive snapshots.In order to prevent H cooling, we need 𝑛 Γ dyn ≥ 𝑓 H 𝑛 Λ H .Since, for short enough time-scales (order of Myr at the high red-shifts considered in this work), the second term on the right-handside of Eq. (6) is small, it can be neglected, translating the conditioninto a critical mass accretion rate (cid:18) d 𝑀 vir d 𝑡 (cid:19) crit = . 𝜇 𝑓 H 𝑛 H ( 𝛾 − ) 𝑘 B 𝛼 − Λ H 𝑀 / . (7) MNRAS000 , 1–15 (2021)
A. Lupi et al.
16 18 20 22 24 26 28 30 z M v i r , m i n ( M (cid:12) ) J = 0 . J = 1 J = 10 Kulkarni+20 ( J = 0 . Kulkarni+20 ( J = 1) Kulkarni+20 ( J = 10) Schauer+2019
Figure 1.
Minimum mass for efficient H cooling from Eq. (5) for differentvalues of 𝐽 (red lines), compared with the results by Kulkarni et al. (2020),shown as blue lines, and the redshift-independent result (for 𝐽 =
0) bySchauer et al. (2019).
We follow the evolution of all mini-haloes in the simulation,and track whether at any stage they exceed the minimum mass for H cooling. For any mini-halo that exceeds this threshold, we furthercheck whether they fulfill the dynamical-heating criterion. For anyhalo above the H cooling minimum mass which, during any stepof the evolution, fails the latter criterion, we assume that it remainsa plausible host of a direct-collapse BH only for one H coolingtime 𝜏 H (Fernandez et al. 2014), after which it is assumed to formstars and is excluded from the candidate list, unless the conditionfor dynamical heating was fulfilled again within the time 𝜏 H . Forexample, a halo could experience a strong accretion or a merger ata given time, events that result in a dynamical-heating rate muchhigher that the critical value, then the accretion drops below thethreshold for a short time, and then increases again above threshold.If the time during which the halo is below the critical threshold isshorter than the halo H cooling time 𝜏 H , we retain the halo as acandidate. While dynamical heating may prevent H cooling and star forma-tion in the mini-halo stage, it is unclear whether it can help avoidcatastrophic cooling (via H atoms and then H ) and the ensuing starformation once a halo crosses the atomic-cooling threshold (ACT),even if the halo is pristine (see more discussion of this below). Aplausible scenario proposed to achieve this, and keep even ACHsfree of H cooling, involves synchronised pairs of ACHs. In thisscenario, once the first halo crosses the atomic cooling threshold(Halo 1), it cools and forms stars in Δ 𝑡 SF = 𝑡 ff ≡ √︁ ( 𝜋 )/( 𝜌 ) ,subsequently providing a strong enough LW flux able to prevent theformation of H in the companion halo (Halo 2). This allows gas inHalo 2 to collapse nearly isothermally at 𝑇 ∼ K, leading to theconditions (i.e. large accretion rate) required for the formation of aDCBH.In order for this process to occur, in our model we require thetime difference Δ 𝑡 between the onset of SF in Halo 1, occurringat 𝑡 = 𝜏 ff , after the ACT crossing, and the ACT crossing in Halo 2 to be shorter than Δ 𝑡 SN = 𝜏 ff , < Δ 𝑡 <𝜏 ff , + Δ 𝑡 SN . In addition, to ensure that the LW flux is large enoughto prevent H formation, we also require the halo separation tobe 200 < Δ 𝑟 <
500 physical pc. Since, in our model, we alsoexplore the case of star formation occurring in mini-haloes that hadnot been sufficiently dynamically heated, in these cases we excludefrom the list of synchronised pair candidates all haloes that hadformed stars during their mini-halo phase. Nevertheless, below wealso compare such results with a simpler approach based on thenumber of progenitors, as discussed in Section 2.5.
An important role in determining the abundance of potentialDCBHs is played by the value of 𝐽 chosen for the analysis. Inthis work, we consider three different cases. In the first two, 𝐽 iskept constant at 0.01 and 1.0, respectively, which brackets the fidu-cial values in Trenti et al. (2009), and values expected for the meanLW background at the redshifts considered here. In the third, weself-consistently compute the local 𝐽 due to nearby star-forminghaloes, which can be much larger than the background, dependingon the halo’s location. We follow Dijkstra et al. (2014) and computethe local flux impinging on each halo as the sum of the fluxes overall star-forming haloes in the tracked region, with the contributionof each halo given by 𝐽 = ℎ (cid:104) 𝜈 (cid:105) Δ 𝜈 𝑓 esc , LW 𝑄 LW 𝜋𝑟 𝑀 ★ , (8)where 𝑀 ★ = 𝑓 ★ ( Ω b / Ω m ) 𝑀 vir . Here we use (cid:104) 𝜈 (cid:105) = . × Hzand Δ 𝜈 = . × Hz as the average frequency and the bandwidthof the LW band, extending from 11.2 to 13.6 eV, and 𝑓 esc , LW is theescape fraction of LW radiation from the halo, which we set to100 per cent. The latter assumption is reasonable in most cases, inparticular far away from the irradiating halo, and when the relativegas velocities are large, hence doppler-shifting the LW photonsout of resonance. On the other hand, the measured LW 𝑓 esc candrop significantly closer to the halo, if the stellar ionisation frontexpands slowly and/or if the stellar spectrum softens, allowing forH to form in the outskirts of the halo and absorb LW radiation(see, e.g., Schauer et al. 2017). Finally, 𝑟 is the distance from theLW source, and 𝑄 LW = 𝑄 [ + 𝑡 LW /( )] − / exp [− 𝑡 LW /(
300 Myr )] (9)is the LW photon production rate per solar mass formed in starsper single burst of SF, at a time 𝑡 LW after a SF burst, with 𝑄 thenormalisation, which we set to 10 s − M (cid:12)− for both PopIII andPopII stars. We adopt the photon yields from Schaerer (2003) andthe same initial mass function for both populations, in order to keepthe same functional form for the time dependence of the LW flux.The decaying trend of 𝑄 LW reflects the fact that massive stars, thatare the main contributors to the LW radiation field, live only for a fewMyr. Before discussing the redshift evolution, we need to point outa possible limitation of this treatment of the LW radiation. Eq. (8) isself-consistent only in the case of SF occurring in a single burst andthen instantaneously stopping. If instead SF occurs continuously, asone would expect, 𝐽 should be determined by the convolution of 𝑄 LW with the time-dependent star formation rate (SFR) in the halo.Unfortunately, such a treatment cannot be easily implemented in a MNRAS , 1–15 (2021) assive seed black holes in quasar hosts semi-analytic framework, hence we opt here for an approximationin which 𝑄 LW in each halo is evolved at every time-step Δ 𝑡 as 𝑄 (cid:48) LW ≈ 𝑄 LW [ + (cid:68) (cid:164) 𝑄 LW 𝑄 LW (cid:69) (cid:48) Δ 𝑡 ] 𝑀 ★ 𝑀 (cid:48) ★ , (10)where the primed quantities correspond to the updated value at 𝑡 LW + Δ 𝑡 , (cid:68) (cid:164) 𝑄 LW 𝑄 LW (cid:69) (cid:48) ≈ − − (cid:18) + 𝑡 (cid:48) LW − (cid:104) 𝜏 (cid:48) LW (cid:105) (cid:19) − (11)is the average decay rate of the LW flux of the entire population,and (cid:104) 𝜏 LW (cid:105) (cid:48) = (cid:104) 𝜏 LW (cid:105) 𝑀 ★ 𝑄 LW + 𝑡 (cid:48) LW ( 𝑀 (cid:48) ★ − 𝑀 ★ ) 𝑄 𝑀 ★ 𝑄 LW + ( 𝑀 (cid:48) ★ − 𝑀 ★ ) 𝑄 (12)is the flux-weighted average time of the SF bursts at 𝑡 (cid:48) . In or-der to accurately follow the variations of 𝑄 LW , we also needto impose a time-step limiter to the integration, which we de-fine as Δ 𝑡 = min { . |(cid:104) 𝑄 LW / (cid:164) 𝑄 LW (cid:105)| , Δ 𝑡 sim } . If the integrationtime-step is shorter than the time interval between the simula-tion snapshots, we sub-cycle Eqs. (10) to (12) assuming 𝑀 (cid:48) ★ = 𝑀 ★ + Δ 𝑀 ★, sim ( Δ 𝑡 / Δ 𝑡 sim ) , with Δ 𝑀 ★, sim = 𝑀 ★ [ 𝑡 + Δ 𝑡 sim − 𝑀 ★ ( 𝑡 )] .It can be easily seen that our approximation reduces to the exact so-lution for a single burst, in which 𝑀 (cid:48) ★ = 𝑀 ★ yields (cid:104) 𝜏 LW (cid:105) (cid:48) = (cid:104) 𝜏 LW (cid:105) ,and Eq. (10) reduces to an explicit time integration of the 𝑄 LW de-cay with time. On the other hand, in the case of a constant SFR with 𝑀 (cid:48) ★ (cid:29) 𝑀 ★ , Eq. (12) reduces to (cid:104) 𝜏 LW (cid:105) (cid:48) (cid:46) 𝑡 (cid:48) LW , that gives an almostexact solution in which 𝑄 LW slowly decays with time accordingto the mass ratio between young stars and old stars. Although notexplicitly shown, we verified the accuracy of our approximation indifferent cases, finding an extremely good agreement always within ∼ 𝐽 case, we always enforce a minimum 𝐽 , min = .
01 to be present in the box, aimed at mimicking a globalbackground from sources outside the box (Trenti et al. 2009). Sucha floor also avoids the divergence in the H fraction calculations,which should be changed accounting for collisional H dissociationin the absence of radiation. The criteria described above are applied to every output of thesimulation, starting from the highest redshift down to 𝑧 =
10, similarto the procedures in semi-analytic models. Although the results forsynchronised pairs and dynamically heated haloes are presented indifferent Figures, they are tracked simultaneously in the code, i.e.,they can co-exist in the same volume for a given set of parameters.In particular, our procedure is the following:(i) We read a new rockstar output containing the list of haloesto add and evolve.(ii) If a halo progenitor is not present (the current halo has beenidentified for the first time), we add it to the list of haloes to evolve.If its virial temperature is below 10 K, the halo is consideredto be a candidate for becoming a DCBH host, according to thecriteria outlined above. Otherwise (this only occurs when the haloforms at the boundary of the high-resolution region, hence it iscontaminated by massive particles and the mini-halo phase is notresolved) it is conservatively considered as star-forming. These star-forming haloes are excluded from the list of DCBH host candidates,but are a potential source of metal contamination (see below). (iii) If a progenitor is present, and is a mini-halo, we determinethe halo accretion rate during the last step and check the dynamicalheating criterion. If the criterion is not matched, the halo is removedfrom the candidate list, but evolved further as a possible source ofcontamination. If the target halo has just crossed the ACT, andthe dynamical heating criterion has always been matched, we flagit as dynamically heated and store it in a dedicated catalogue ofcandidates.(iv) If, on the other hand, the progenitor was already above theACT, the halo is simply updated, accounting for star formation andmetal enrichment, the latter computed according to Dijkstra et al.(2014).(v) among all updated haloes, we identify those that have justcrossed the ACT, and look for nearby haloes matching the synchro-nisation criteria described above. Because of tidal stripping duringhalo mergers, it can occur that some haloes keep oscillating aroundthe ACT for some time after their first crossing. In order to avoid ourprescription to spuriously identify pairs multiple times, we removethe pairs from the candidate list as soon as they have crossed theACT for the first time. Also for synchronised pairs, we store the pairproperties in a catalogue.(vi) We then repeat the above for all haloes, checking each halofor metal pollution inhibiting the formation of DCBHs (see nextsubsection);(vii) Finally, before moving to the next snapshot, we check formergers among the haloes, identifying those haloes sharing the samedescendant; then, we update the list of active haloes by creatinga unique virtual descendant for each merger, represented by itsmost massive progenitor; we then consistently update the numberof progenitors of this virtual descendant, and also take into accountwhether any of the progenitors had already crossed the ACT in thepast or was contaminated with metals.
An important process that can affect our results is the possiblecontamination by metals, via genetic pollution (when progenitorsalready had metals) and via environmental pollution (when themetal bubbles created after SN explosions in nearby galaxies reachthe target halo). In our model, we consider both of these processes.For genetic pollution, in this work we consider two differentmodels. A first approach consists of tracking only the number ofprogenitors 𝑁 prog of each candidate halo, and, assuming that eachprogenitor has a fixed 10% probability of becoming star-forming(SF) during its history (Dijkstra et al. 2014). The advantage of thissimple treatment is that SF in mini-haloes does not need to befollowed explicitly. Instead, we filter the candidate list by applyingtwo different criteria: in order to avoid metal-pollution, we requireeither 𝑁 prog < 𝑁 prog <
5, which correspond to less than 20% or50% probability of being star forming during the mini-halo phase,respectively. We obtained results from this simple approach as anacademic exercise, and in order to be able to compare results withprevious work by Dijkstra et al. (2014). In the second approach,instead, which we employ as our fiducial case, we self-consistentlyaccount for SF in mini-haloes that can cool via H , and excludethese from the list of candidate dynamically heated mini-haloes.For environmental pollution, we follow all haloes able to formstars in our cosmological volume, i.e. mini-haloes not dynamicallyheated and all ACHs, and compute the size of the metal bubblestarting from the time of the first SN explosion. Hence, the metalbubbles start to expand in our model Δ 𝑡 = 𝑡 ff + Δ 𝑡 SN after the halois first declared star-forming. MNRAS000
5, which correspond to less than 20% or50% probability of being star forming during the mini-halo phase,respectively. We obtained results from this simple approach as anacademic exercise, and in order to be able to compare results withprevious work by Dijkstra et al. (2014). In the second approach,instead, which we employ as our fiducial case, we self-consistentlyaccount for SF in mini-haloes that can cool via H , and excludethese from the list of candidate dynamically heated mini-haloes.For environmental pollution, we follow all haloes able to formstars in our cosmological volume, i.e. mini-haloes not dynamicallyheated and all ACHs, and compute the size of the metal bubblestarting from the time of the first SN explosion. Hence, the metalbubbles start to expand in our model Δ 𝑡 = 𝑡 ff + Δ 𝑡 SN after the halois first declared star-forming. MNRAS000 , 1–15 (2021)
A. Lupi et al.
To model SF, we assume that the stellar mass formed in ACHs is 𝑓 ★ = .
05 of the baryonic mass, translating to 𝑀 ★ = 𝑓 ★ Ω b / Ω m 𝑀 vir .For mini-haloes, instead, we consider a lower star-forming efficiency 𝑓 ★ = . 𝜌 gas = 𝜂 gas 𝜌 max ) has elapsed since (i)dynamical heating conditions have failed or (ii) the ACT has beenexceeded (if no previous SF in the mini-halo had already occurred),and then proceeds smoothly following the increase in stellar masswith time in the halo. Following Weaver et al. (1977) and Madauet al. (2001), we define the bubble radius as 𝑅 bubble ( 𝑡 ) = (cid:18) 𝜋 (cid:19) / (cid:18) 𝐿 SN m H (cid:19) / 𝑡 / cm , (13)where 𝐿 = const is the average SN luminosity associated to thestellar population in the halo, 𝑡 the time since the burst of SF,and 𝑛 ≡ Δ gas Ω b 𝜌 crit ( + 𝑧 ) / 𝑚 H is the gas density (in cm − ) thebubble expands through. Now, we assume for simplicity that 𝐿 SN ≈ 𝐸 SN 𝜂 SN 𝑀 ★ ( 𝑡 )/ 𝑡 , where 𝐸 SN = erg is the energy releasedper SN, 𝑀 ★ ( 𝑡 ) = 𝑓 ★ ( Ω b / Ω m ) 𝑀 vir ( 𝑡 ) is the stellar mass at time 𝑡 , and 𝜂 SN ≈ .
011 M (cid:12)− is the number of SNe per stellar massformed, and replace them in the bubble radius equation, we obtain,for 𝑡 > Δ 𝑡 SN , 𝑅 bubble ( 𝑡 ) ≈
23 pc (cid:34) 𝑓 ★ Ω b Ω m 𝑀 vir ( 𝑡 ) M (cid:12) Myr 𝑡 (cid:18) 𝑡 − Δ 𝑡 SN Myr (cid:19) (cid:35) / 𝑛 − / . (14)We notice that this equation reduces, for 𝑡 (cid:29) Δ 𝑡 SN , to the equationreported in Dijkstra et al. (2014), apart from the slightly lower nor-malisation. In Dijkstra et al. (2014), the baryon overdensity throughwhich the bubble expands is assumed to be Δ gas ≈
60, typical of theoverdensities found within haloes (60 (cid:46) Δ gas (cid:46) Δ gas ∼ −
10. This results in the bubble size being underesti-mated (Dijkstra et al. 2014; Habouzit et al. 2017). In our model,instead, we also account how the overdensity changes outside thehalo, by fitting the asymptotic case of an extremely rare halo with 𝑀 vir = M (cid:12) from Barkana (2004), and simply express Δ gas as Δ gas = (cid:40) 𝑅 ≤ . 𝑅 vir − . 𝑥 + . 𝑥 − . 𝑥 + . otherwise , (15)where 𝑥 ≡ log ( 𝑅 / 𝑅 vir ) , 𝑅 is the bubble radius, and the factor0.787 has been chosen to join the inner constant value and theextrapolation of the fitting function without discontinuity. To con-sistently find the bubble radius outside 0 . 𝑅 vir , when Δ gas alsodepends on 𝑅 bubble , we iterate until convergence. At every step,we flag all candidate haloes that are enriched by metal bubbles, bychecking whether the bubble overlaps with the halo, i.e. the haloseparation 𝑟 < 𝑅 bubble − 𝑅 vir , candidate . A possible caveat in thisevolution comes from the stellar mass scaling in the bubble radiusequation, which does not consistently keep track of the superposi-tion of metal bubbles from multiple SF events in the same halo, andmight result in either an underestimate or an overestimate of theactual bubble size, depending on how the bubbles interact. In this section, we assess the importance of the two mechanismsdiscussed above for the formation of DCBH host candidates. We
Table 1.
List of the models explored, with the differences in their assump-tions. The first column is the name of the model, the second the time delaybetween when a halo is declared star-forming and the formation of the firststellar population, and the third is the time delay before the first SNe ex-plode. The fourth column reports the density in units of the maximum densityachievable via adiabatic contraction and the fifth the assumed value of theLW flux. The first two models (top of the table) exclude SF in mini-haloes,while the other four models (bottom of the table) include SF in the subset ofmini-haloes which can cool, according to the criteria discussed in the text.Model name Δ 𝑡 SF Δ 𝑡 SN 𝜂 gas 𝐽 Basic models: no SF in mini-haloesVisbal 10 Myr 10 Myr - -NoMini 𝑡 ff ( 𝑛 gas ) 𝑡 ff ( 𝑛 gas ) 𝑡 ff ( 𝑛 gas ) 𝑡 ff ( 𝑛 gas ) 𝑡 ff ( 𝑛 gas ) first focus on the synchronised-pair scenario, and then discuss therole of dynamical heating, but we again stress that we follow bothat the same time in our models. A list of the cases we explored,with the differences in their assumptions, is summarised in Table 1.In particular, while in model ‘Visbal’ we assume the same SF andsynchronisation parameters as in Visbal et al. (2014b), all othermodels employ different SF-synchronisation criteria, i.e. we assumethat SF occurs one free-fall time after the halo has been declaredeligible for SF, with the free-fall time estimated at a density 𝜌 gas = 𝜂 gas 𝜌 max . For our fiducial model, we assume 𝜂 gas = .
2, consistentwith the results in Visbal et al. (2014a), and explore the case of 𝜂 gas = In order to disentangle the effect of the large-scale overdensity fromthe additional physical processes included in our model, we firstconsider a simplified model, which we dub ‘Visbal’, in which wefollow the prescriptions employed in Visbal et al. (2014b), and laterinclude the additional physical processes and criteria described in§ 2. To give an idea of the number of halo candidates in the high-resolution region of our reference simulation, we show in Fig. 2 theredshift distribution of haloes crossing the ACT down to 𝑧 =
10 (toppanel). The first few hundred haloes appear already above 𝑧 = 𝑧 = 𝑧 = ( ) . To appreciate the relevance of these numbers, they have tobe normalised by volume. In the middle panel we show the effectivevolume of the Universe sampled as a function of redshift, identifiedas a sphere with radius the mean maximum separation along 𝑥, 𝑦 ,and 𝑧 . Compared to full box simulations, our run is a zoom-insimulation of an overdense region, in which a Lagrangian volumecentred around the target halo is selected, and then evolved with We note that the simulations in Visbal et al. (2014b) employed 𝜎 = . 𝜎 = . , 1–15 (2021) assive seed black holes in quasar hosts time following the collapse of such region. This results in a high-resolution volume that shrinks significantly with time, from about ∼
150 comoving Mpc at 𝑧 =
100 down to 8 cMpc at 𝑧 =
6, thelatter corresponding to the volume of a sphere with radius 2.5 timesthe virial radius of the target halo ( 𝑅 ∼
69 kpc at 𝑧 = The comparison with our region is shown in the bottompanel, with our results as a black histogram and the average den-sity region’s as a red dashed line. Such a large density of ACHs,especially at very high redshift, is unique to the overdense regionswhere 𝑧 = ∼ 𝜎 peak, where our simulated region produces about a 100 times largernumber of ACHs with respect to an average density one. The underlying N-body simulation used in our study represents astrongly biased region, which evolves into an ultra-rare ∼ M (cid:12) halo by 𝑧 =
6. As a result, this region is highly overdense - by a factorof ≈
20 compared to the global average at 𝑧 =
6. In order to assessthe impact of this bias, we first compare our results with those ofVisbal et al. (2014b), who considered the formation of synchronised(sub)halo pairs in a typical region of the Universe in a (15 Mpc) box. For this analysis, we make the same assumptions as Visbal et al.(2014b), i.e. (i) star formation is suppressed in all mini-haloes, (ii) the first halo to cross the ACT starts to form stars 10 Myr after thecrossing, (iii) the second halo must cross the ACT within 10 Myrfrom the SF event, and (iv) the halo pair separation should fall withinthe ranges (a) [0.2-0.5] kpc, (b) [0.2-0.75] kpc, or (c) [0.2-1.0] kpc.The lower limit on the separation (0.2 kpc) corresponds to a distancewithin which ram-pressure and/or tidal stripping could prohibit gascollapse and DCBH formation in the target halo, while the upperlimit of 0.5-1 kpc represents the uncertainty in the distance overwhich the critical LW flux may extend.The results of this comparison are reported in Fig. 3. The toppanel of this figure shows the redshift distribution, the middle onethe cumulative number of synchronised pairs, and the bottom onethe number of pairs per unit volume per redshift bin we found.The corresponding results by Visbal et al. (2014b), obtained in avolume of (15 cMpc) are also shown as coloured crosses in thebottom panel. Cyan, orange, and green colours correspond to thethree separation intervals considered, respectively, as labelled in thelegend of the figure.The number of synchronised pairs in our simulation can reachup to a few hundred down to 𝑧 =
10 in the most optimistic case (greenhistogram), with a rate that varies between a few (blue histogram)and a few tens (green histogram), depending on the allowed rangeof separations (the smaller the maximum separation, the lower thenumber of candidates). Compared to the results in Visbal et al.(2014b), we observe the same trend as a function of the maximumseparation allowed, while the overall number of candidates in therange 𝑧 = [ − ] is 4, 11, and 29 for separations (a), (b), and (c),respectively. For comparison, Visbal et al. (2014b) identified 2, 5,and 17 pairs in their box in the interval Δ 𝑧 = .
25 around 𝑧 = Δ 𝑧 = ≈
45 cMpc at To compute the Halo mass function, we employed the public tool colos-sus (Diemer 2018). d N / d z V e ff ( c M p c ) . . . . . . . z − − d N / d z d V e ff ( c M p c − ) Figure 2.
Redshift distribution of the haloes crossing the ACT in our effectiveLagrangian volume (top panel). The effective volume is ≈
150 comovingMpc at z=100, but subsequently shrinks as the region collapses and as asmaller fraction is tracked at high resolution (middle panel). The bottompanel shows the formation rate of ACHs per unit volume in our simulation(black histogram, compared to the average ACH formation rate density inthe background Universe (obtained from Watson et al. 2013; red dashedcurve). Overall, the number density of ACHs in our volume is about a factorof 100 higher than in an average region of the Universe. 𝑧 =
10, these correspond to 0.11, 0.27, and 0.9 pairs for separations(a), (b), and (c), respectively. This demonstrates that the overdenseenvironment of our simulation favours the formation of such systemsby a factor of ∼ −
40. This comparison is based on the effectivevolume at 𝑧 =
10, and some of the difference is due to the factor of ∼ at 𝑧 =
100 instead, we find that our biased region has a factorof 6–14 more synchronised pairs compared to Visbal et al. (2014b).
Having checked how the large-scale cosmological environment af-fects the total number of synchronised pairs formed, we next con-sider the role of genetic and environmental pollution. In this case,we consider our NoMini model, in which star formation occurs onefree-fall time (estimated assuming 𝜂 max = .
2) after the ACT cross-
MNRAS000
MNRAS000 , 1–15 (2021)
A. Lupi et al. − d N p a i r s / d z [0.2-1.0][0.2-0.75][0.2-0.5] N p a i r s ( ≥ z )
10 15 20 25 z − − − d N p a i r s / d z d V ( c M p c − ) Figure 3.
Redshift distribution of the 73, 251, and 496 synchronised pairsfound in our biased region, following the semi-analytic prescriptions in Vis-bal et al. (2014b) for the three spatial separations between the pair membersconsidered: [0.2-0.5] (green histogram), [0.2-0.75] (orange histogram), and[0.2,1.0] kpc (blue histogram), respectively. The top panel shows the for-mation rate (per unit redshift) in the simulated volume, the middle one thecumulative number of pairs, and the bottom one the formation rate densityper comoving unit volume. The three crosses with the same colour schemeshow the corresponding results by Visbal et al. (2014b) in an average densityvolume of (15 cMpc) . ing, and a synchronisation time-scale of 5 Myr, corresponding tothe time when the first SNe explode. In addition, we require the haloseparation to fall in the range [0.2-0.5] kpc. Also in this case, weassume that no prior SF occurs in mini-haloes.The results are shown in Fig. 4, where we report the totalnumber of pairs in black (34) and the genetic-pollution filteredpairs in blue (34, requiring 𝑁 prog <
5) and in orange (19, with thestricter requirement of 𝑁 prog < ∼ 𝑡 ff ( 𝜂 gas = . ) ; typically slightly longerthan 10 Myr) here; the maximum delay time is also 10 Myr (in Visbalet al. 2014b) and 5 Myr in our case. In facts, this difference results d N p a i r s / d z All N prog < N prog < N prog < . . . . . . . z N p a i r s ( ≥ z ) Figure 4.
Redshift distribution of the 34 synchronised pairs in our NoMinimodel (black histogram), with the additional criteria imposed to excludegenetically (blue and orange histograms for 𝑁 prog < 𝑁 prog < 𝑁 prog < 𝑁 prog < in a shorter time window and a longer time separation between theACT crossing in the two haloes of the pair, both concurring to thereduction of the number of candidates.If we now consider the additional criteria, we observe that theblue histogram completely overlaps with the black one, demonstrat-ing that all candidate haloes have fewer than 5 progenitors. On theother hand, only ∼
50% remain if we require a single progenitor.In addition to genetic pollution, environmental pollution due to ex-panding metal bubbles from nearby haloes appears very important,reducing the number of synchronised pairs by about a factor of 3.Another important effect of these additional criteria can be observedin the redshift distribution of the candidates, which becomes sig-nificantly skewed towards higher redshift in the most stringent case(purple histogram), i.e. almost all the pairs form at 𝑧 >
We next consider all of the physical processes described in § 2,including SF in the mini-haloes. In this case, we have to considertwo additional parameters that play a crucial role in preventingSF in low-mass haloes, i.e. the maximum density reached by thegas via adiabatic compression (defined by the 𝜂 parameter) andthe LW radiation background, which sets the minimum halo massable to efficiently form H and cool (defined by the 𝐽 parameter;see, e.g. Machacek et al. 2001; Trenti et al. 2009; Schauer et al.2017; Kulkarni et al. 2020, for a discussion). In our parameterexploration, we assume 𝜂 = { . , . } and 𝐽 = { . , . } , plus MNRAS , 1–15 (2021) assive seed black holes in quasar hosts z p a i r s Fiducial LowFlux HighFlux HighDensity
Model N p a i r s ( z ≥ ) J = local − η ρ = 0 . J = 0 . − η ρ = 0 . J = 1 − η ρ = 0 . J = local − η ρ = 1 Figure 5.
Comparison among the different model parameters (each columncorresponding to a parameter combination, as reported in the legend). Thetop panel shows the redshift distribution of the synchronised pairs as grey-filled circles, whereas the coloured ones correspond to the haloes that havenot been environmentally polluted. In these four models, the genetic pollu-tion criterion is automatically verified by the lack of SF in one of the haloesof the pair. In the bottom panel, we show the total number of pairs down to 𝑧 =
10, with the grey histogram corresponding to the total number, and thecoloured one to chemically pristine ones only. an additional case where we allow 𝐽 to vary locally according toEq. 8. An important difference from the basic models is that here,the possibility of SF to occur in mini-haloes allows us to naturallyaccount for genetic pollution resulting from SF events in the haloprogenitors, leaving environmental pollution the only additionalcheck to be applied a posteriori.In Fig. 5, we compare our fiducial model with alternativechoices of the parameters. In the top panel, we report the redshiftdistribution of the synchronised pair candidates, shown as grey cir-cles, with the pristine haloes identified by the colours circles. Ourfiducial model is in red, the models with constant 𝐽 are in blue(LowFlux; 𝐽 = .
01) and orange (HighFlux; 𝐽 =
1) respec-tively, and the model with variable 𝐽 and 𝜂 𝜌 = 𝑧 =
10, with the totalnumber in grey and the coloured histograms corresponding to thehaloes that have not been environmentally polluted. While coveringthe redshift range 11 ∼ < 𝑧 ∼ <
20, with 27 candidate pairs identified,the formation of suitable pairs in our fiducial model is concentratedaround 𝑧 ∼
15, with the only pristine (hence most plausible) candi-dates appearing at 𝑧 ∼ . 𝑧 ∼ . 𝐽 is assumed, instead, almost all pairs are able to form stars (withno metal-free cases at all), because of the early formation of many stars in mini-haloes and the subsequent pollution by SNe (secondcolumn, LowFlux). This is also confirmed by the 𝐽 = 𝐽 becomes less ef-fective, with 19 pairs forming, (rightmost column, HighDensity).Moreover, only one of them is pristine, making this case slightlyless favourable.Note that all numbers above include haloes found anywhere inour simulation box. Some of the candidates we identified up to thispoint do not end up in the final massive quasar-host halo by 𝑧 = As discussed above, a possible common path to DCBH formationinvolves the suppression of H cooling (and subsequent SF) inmini-haloes via dynamical heating (Yoshida et al. 2003; Inayoshiet al. 2018). Furthermore, it has recently been suggested that thesedynamically-heated halos may form very massive stars and leavebehind massive BH seeds, even though they form H and cool oncethey reach the ACH stage (Wise et al. 2019; Sakurai et al. 2020) thusforming several stars that compete for accretion (Regan et al. 2020b),so that the seed mass is probably lower than in the synchronised paircase. Here, we consider this channel and report the results obtainedin our framework, with otherwise the same set of assumptions as inthe four different models discussed in the previous section.In Fig. 6, we show four panels, the top ones reporting theredshift distribution of dynamically heated mini-haloes, and thebottom ones the cumulative number down to 𝑧 =
10. In the left-hand panels, we show the total number, without any constrainton the metallicity within the mini-halo. In the right-hand panels,we report the effective numbers after removal of the metal-pollutedcandidates. Our fiducial model is again reported in red, the LowFluxconstant low- 𝐽 case in blue, the HighFlux constant 𝐽 = 𝐽 case with 𝜂 𝜌 = 𝐽 produced by the few star-forming mini-haloes and by thelarge population of ACHs (which form stars at a much higher ratethan mini-haloes), and the relatively compact size of the regionconsidered (with volume ∼ at 𝑧 =
6, enclosing a totalmass of ∼ . × M (cid:12) ). The number of candidates significantlydrops when we consider either a higher gas density (2,698; greenhistogram, HighDensity) or a constant moderately high LW back-ground (228; orange histogram, HighFlux). In the most pessimisticcase LowFlux of an extremely low LW background, instead, thenumber of dynamically heated haloes becomes extremely small(10; blue histogram), because of the low minimum halo mass ableto cool via H cooling and the relatively high dynamical heatingrate needed to overcome cooling. MNRAS000
6, enclosing a totalmass of ∼ . × M (cid:12) ). The number of candidates significantlydrops when we consider either a higher gas density (2,698; greenhistogram, HighDensity) or a constant moderately high LW back-ground (228; orange histogram, HighFlux). In the most pessimisticcase LowFlux of an extremely low LW background, instead, thenumber of dynamically heated haloes becomes extremely small(10; blue histogram), because of the low minimum halo mass ableto cool via H cooling and the relatively high dynamical heatingrate needed to overcome cooling. MNRAS000 , 1–15 (2021) A. Lupi et al. d N h e a t e d / d z Including metal polluted haloes
AC haloes J = local − η gas = 0 . J = 0 . − η gas = 0 . . . . . . . . z N h e a t e d ( z ≥ ) J = 1 . − η gas = 0 . J = local − η gas = 1 . d N h e a t e d / d z Excluding metal polluted haloes
AC haloes J = local − η gas = 0 . J = 0 . − η gas = 0 . . . . . . . . z N h e a t e d ( z ≥ ) J = 1 . − η gas = 0 . J = local − η gas = 1 . Figure 6.
Redshift distribution (top panels) and cumulative number (bottom panels) of dynamically heated haloes in our model, including (left panels) orexcluding (right panels) the metal-polluted haloes. The light grey histograms in the background correspond to the total population of ACHs formed in the run,and the different colours to a specific 𝐽 − 𝜂 gas combination, as reported in the legend. If we apply the metal pollution constraints, the constant 𝐽 cases see a net reduction in number (to about 27% in our fiducialmodel and about 10% in the other cases), which result in the lackof any dynamically heated haloes in the 𝐽 = .
01 case, and ashift to higher redshifts of the peak of the distribution. As for thevariable LW flux cases, instead, our reference case exhibits a oneorder of magnitude decrease in number, with the peak reachedaround 𝑧 ∼
15, whereas the higher 𝜂 𝜌 case is almost completelysuppressed, because of the typically shorter cooling timescales inthe haloes, resulting in higher SFRs and subsequent metal pollutionof their neighbourhood. In general, in most of the cases explored,the number of dynamically heated haloes is larger than that ofsynchronised pairs, suggesting that the conditions to dynamicallyheat mini-haloes are easier to realise than those guaranteeing halosynchronisation, especially in the overdense environment typical ofmassive galaxies at high-redshift.To assess the importance of the self-consistent 𝐽 estimationin our analysis, we report in Fig. 7 the 𝑥 − 𝑦 projection of the 𝐽 fieldin our cosmological volume at 𝑧 =
10. From the map, we can clearlynotice that, because of the large over-density in the region ( 𝛿 = . 𝑧 = 𝜎 peak), many galaxies areforming, and 𝐽 reaches quite high values compared to an averagedensity region of the Universe. In particular, the most massive haloesin the region, which form at the centre of the map, are able toproduce a 𝐽 as high as 500 or more, although we also expect thesame region to be strongly polluted by the metals produced by SNe, The overdensity 𝛿 in our simulation at 𝑧 =
100 has been derived bycreating a convex hull around our high-resolution region, defining the re-gion’s volume, and then deriving its average density relative to the cosmicaverage density. The peak height 𝜈 ≡ 𝛿 / 𝜎 is then obtained using the massvariance 𝜎 = . 𝑀 = . × M (cid:12) region assuming the PlanckCollaboration et al. (2016) cosmological parameters. hence not necessarily representing the most favourable location toform DCBHs (at 𝑧 = 𝐽 for the dynamically heated haloes in our reference model (redline), compared with the median value for the entire population ofACHs (black line), as a function of redshift. In the overdense regionwe consider, the value for the general population of ACHs is alreadyhigher than that of typical high-redshift star-forming regions (e.g.Dayal & Ferrara 2018), reaching median values of about 100 below 𝑧 =
15. While such enhancement is enough to guarantee that mostof the mini-haloes are dynamically heated below 𝑧 =
15, only thehaloes exposed to a 2–to-10 times higher than average LW flux canbe dynamically heated at higher redshift. This can be explainedwith these haloes forming in denser filamentary regions (see Fig. 7)where the mass accretion rate is higher, and the nearby SF is able toreduce H formation sufficiently to prevent molecular cooling. Ourresults are also in agreement with previous studies of high-redshiftmassive haloes by Valiante et al. (2016), where the 𝐽 distributionextended to higher values than usual because of the higher SFRin the region, but also because of the larger halo clustering, aneffect that cannot be tracked in semi-analytic models based on theextended Press–Schechter formalism (Lacey & Cole 1993) and itsvariants.Although the high LW flux seems to suggest that dynamicalheating is much more favoured compared to synchronised pairs, wehave to keep in mind that as soon as the halo exceeds the ACT,H starts to form more efficiently because of the large increase indensity facilitated by atomic cooling, and a 𝐽 that was sufficientto prevent cooling in the mini-halo phase is not guaranteed any-more to avoid H formation in the ACH (Omukai & Palla 2001;Oh & Haiman 2002). Additionally, any metal-pollution can resultin even more efficient cooling in the ACH stage. Nevertheless, theresults in Wise et al. (2019) suggest that large accretion rates to- MNRAS , 1–15 (2021) assive seed black holes in quasar hosts Figure 7.
Projection of the 𝐽 distribution in our high-resolution region at 𝑧 =
10. While the low-density regions exhibit a low 𝐽 , as expected, the mostmassive filaments and knots in the cosmic web, where most of the haloes form, show a significantly higher LW flux, reaching a peak of more than 500 in thecentre, where the most massive haloes (progenitors of the expected quasar host), i.e. the most rapidly star-forming sites, form. Such a high LW flux does notnecessarily imply a more likely formation of DCBHs in the central region, because of the correspondingly higher metal pollution due to the large number ofSNe occurring in the same region. wards the centre of the halo may persist in the ACH stage, and,as mentioned above, spherically symmetric zoom-in simulations ofthese candidates (Sakurai et al. 2020) find that radiation from thegrowing central protostar does not shut down this rapid flow, whichproduces a supermassive star. Regan et al. (2020a) also find thatrapid inflow may occur even in the face of some metal pollution,but Regan et al. (2020b) find that widespread star formation limitsthe final BH mass to < M (cid:12) .To summarise our results, we report in Table 2 the number ofcandidate DCBH haloes in the dynamical-heating (DH; second col-umn) and in the synchronised-pairs scenarios (SP; third column).In the last two columns we also report the number of pristine can-didates for these two cases, respectively. Although we have assessed the number of candidate progenitorhaloes expected to form in the overdense region that subsequentlyevolves to a high-redshift massive halo, our analysis so far has not
Table 2.
Number of candidate DCBH haloes in the dynamical-heating (sec-ond column) and synchronised-pair (third column) scenarios for our mod-els. The corresponding number of metal-free subsets of these candidates areshown in the fourth and fifth columns. Note that not all of these candidatesend up in the massive quasar-host halo by 𝑧 = considered the subsequent evolution of these haloes, i.e. if theyare going to merge with the main progenitor of our target halo by 𝑧 = MNRAS000
Number of candidate DCBH haloes in the dynamical-heating (sec-ond column) and synchronised-pair (third column) scenarios for our mod-els. The corresponding number of metal-free subsets of these candidates areshown in the fourth and fifth columns. Note that not all of these candidatesend up in the massive quasar-host halo by 𝑧 = considered the subsequent evolution of these haloes, i.e. if theyare going to merge with the main progenitor of our target halo by 𝑧 = MNRAS000 , 1–15 (2021) A. Lupi et al.
10 12 14 16 18 20 22 24 z h J i ACHsDyn. heated ACHs
Figure 8.
Median 𝐽 at the locations of the dynamically heated haloes (inred) and all the ACHs (in black) as a function of redshift for our fiducialmodel. The black line shows that the median 𝐽 in the region quickly rises(above 𝑧 =
20) and then settles around 100, a value much higher than thatfound in typical star forming regions at high redshift, but comparable tothe values found in previous studies around massive quasar hosts (Valianteet al. 2016). The dynamically heated population exhibits instead a factor of2–to–10 higher LW flux above 𝑧 =
18, consistent with their being in denserfilaments around the central halo, where the number of haloes and the SFRare higher, hence producing a stronger local LW flux, able to suppress H formation and cooling. and grow to ∼ M (cid:12) , as observed by 𝑧 ≈
6, we checked how manysynchronised pairs and dynamically-heated haloes belong to themerger history of our target halo. The results are reported in Table 3for the same four cases discussed above, with the main results of ouranalysis highlighted in bold. For both scenarios, about 30–to–70 percent of the candidates appear as progenitors of the target halo, i.e.9 synchronised pairs and 5,463 dynamically heated haloes in ourfiducial model.When accounting for metal enrichment, excluding pollutedhaloes from the candidate list, the numbers get reduced to aboutone third (or less) of the total, which translates in ‘only’ one syn-chronised pair and 1,389 dynamically heated haloes in our fiducialmodel. The large number of dynamically heated haloes can be eas-ily explained by the strong 𝐽 produced by the highly star-forminggalaxies in the simulated region. These galaxies create a large LW"bubble" (i.e. a region in which the local LW flux dominates theglobal background) around the target halo, able to more efficientlysuppress H formation in the mini-haloes around it, hence boostingthe number of dynamically-heated mini-haloes in its surrounding.This does not mean that the seed BH formation via dynamical heat-ing is less likely in the target halo (the number is typically higherthan that of synchronised pairs anyway), simply that the overdensityboost affects also nearby halos. The formation of many seeds withina relatively small volume is favourable to the eventual formationof MBH pairs as these halos merge with one another during theircosmic assembly. Despite the large parameter variation we explored here, there aresome caveats in our analysis that must be considered. First, the
Table 3.
Number of synchronised pairs (SP) and dynamically heated (DH)haloes belonging to the merger history of the target halo in our simulation(QSO prog) in the four models considered in the analysis (second and thirdcolumns). The number of pristine halos among these candidates is shown inthe fourth and fifth columns, respectively, and is highlighted in bold face asthe main result of this paper. The number in the parentheses corresponds tothe pristine fraction (percentage) of these progenitors.Model QSO prog Pristine QSO prog (%)DH SP DH SPFiducial 5463 9 (25.4%) (11.1%)LowFlux 7 1 ( 0.0%) ( 0.0%)HighFlux 90 8 ( 4.4%) (37.5%)HighDensity 1359 11 ( 0.5%) ( 9.1%) assumption of an 100% LW radiation escape fraction can resultis a more effective suppression of H formation, hence cooling,boosting the population of both dynamically heated haloes andsynchronised pairs (although the effect is more moderate for thelatter case). Indeed, while in low-mass haloes gas is strongly affectedby stellar feedback, likely leaving LW radiation to escape efficiently,more massive haloes could instead exhibit lower escape fractions,reducing the effective 𝐽 reaching the nearby mini-haloes (see, e.gSchauer et al. 2017).Secondly, the model we employed for metal bubble propaga-tion, despite the improvements we made by considering the lowerdensities outside the haloes, is still very approximate and does notconsider how the bubble expansion velocity and shape change dueto mixing (see Habouzit et al. 2016b, for a discussion). This couldresult in either an overestimation or underestimation of the metalpollution, which then affects gas cooling and the possible forma-tion of DCBHs. Moreover, the assumption of a constant SFR inEq. (14), equal to the average SFR since the first SF burst could alsoresult in an overprediction/underprediction of the bubble size, sinceyounger stars would create metal bubbles at later times in possiblydifferent density conditions, and would then need time to propagateoutwards.Third, although the mass resolution in our simulation is alreadyhigh, allowing us to properly resolve the ACT, it is still too low toresolve the lowest-mass minihalos, near the minimum halo mass forH cooling for low and moderate LW fluxes. This inability might inprinciple result in the overprediction of the DCBH candidates (seeSection 2). Our results highlight how a non-negligible number of synchronisedpairs can form in the overdense regions where quasar hosts areexpected to form, as long as the background LW radiation fieldis larger than 𝐽 = .
01. Although most of these pairs would bealready contaminated by metals in our models, employing simplemodels for metal propagation and mixing, we find that some ofthese haloes remain plausibly pristine DCBH host candidates. Inour fiducial case, we find that only two pairs match all the conditionsto be a plausible DCBH host candidate, of which only one belongsto the merger history of the quasar host. If we relax the spatialseparation condition for synchronised pairs from [0.25-0.5] kpcto the intervals explored by Visbal et al. (2014b), as shown inSection 3.1.1, the number of plausible candidates belonging to thequasar host merger history increases from 9 to 48 (up to 750 pc)
MNRAS , 1–15 (2021) assive seed black holes in quasar hosts and 94 (up to 1 kpc), of which 7 and 19 are pristine, respectively. Ingeneral, our findings confirm that the large overdensity in the regionwhere a massive quasar host later forms exhibits ideal conditionsfor the formation of close pairs. Because of the semi-analytic natureof our framework, a detailed analysis of the gas evolution within thehalo is not possible, and is deferred to a future study.As for dynamical heating, our results suggest that a very largenumber of dynamically heated ACHs could form in the overdenseregions where quasar hosts subsequently assemble (14,698 in ourfiducial model, of which 1,522 appear between 𝑧 =
16 to 𝑧 = ∼ × M (cid:12) halo at 𝑧 =
6. Wise et al. (2019) considered a somewhat less overdenseregion (a Lagrangian volume centred around an 3 × M (cid:12) halo at 𝑧 =
6) and found about 600 such haloes between 𝑧 =
16 and 𝑧 = 𝐽 (cid:38)
1) is present, able tosimply reduce the equilibrium fraction of H forming in the haloes,rather than completely suppressing it. This result is also in line withWise et al. (2019), where the dynamically heated haloes exhibited 𝑓 H ∼ − − − , much smaller than the values in the absenceof any LW flux (Haiman et al. 1996; Tegmark et al. 1997). On theother hand, if no LW radiation is present, the accretion rates neededto balance H cooling would be unrealistically large.Overall, we note that the fulfilment of the dynamical heatingcriterion does not guarantee that the ACH will collapse withoutfragmenting. This is because even if the heating rate matches theH cooling in the mini-halo stage, the subsequent evolution of theACH is more uncertain. Indeed, once the ACT is crossed, coolingby atomic H lines (and metals, if present), together with the rapidchange in the thermodynamic conditions (gas becomes ionised afterthe ACT has been exceeded, just before atomic cooling is activated),could still boost the production of H , hence favouring fragmenta-tion and reducing the accretion rate towards the centre, preventingthe formation of a DCBH. This is in line with Wise et al. (2019),who found that in the DCBH host candidates they identified intheir hydrodynamic simulations, dynamical heating is insufficientto prevent H formation during the ACH stage. Nevertheless, in afew such haloes, they also found that the accretion rates remainedhigh, plausibly leading to the formation of a supermassive star anda DCBH. Sakurai et al. (2020) recently followed the subsequentevolution of these candidates in spherically symmetric radiation-hydrodynamic simulations, and indeed found that the collapse isnot prevented by radiative feedback, and produces a supermassivestar, suggesting that dynamical heating alone may be a sufficientcondition to lead to the formation of a DCBH. In a more realisticcosmological environment, however, Regan et al. (2020b) find thataccretion is suppressed shortly after the formation of the first star bythe formation of several other stars that compete for accretion andby the motion of the star(s) in the messy non-spherical potential ofhigh-redshift galaxies, resulting in DCBHs of about 10 M (cid:12) . Whilesuch masses can be enough to explain billion solar mass MBHs at 𝑧 ∼
6, they still represent a potential problem to explain MBHsabove 10 M (cid:12) at 𝑧 >
7, even assuming an optimistic sustainedaccretion at the Eddington limit.Given these possible limitations on the accretion rate in dy-namically heated mini-haloes, pair synchronisation still appears tobe the best case to trigger massive DCBH formation. While we find J N h a l o e s DH candidatesPristine DH candidates
Figure 9. 𝐽 distribution in dynamically heated mini-haloes (at the time ofthe ACT crossing) from our fiducial model, including (blue histogram) orexcluding (orange histogram) metal polluted haloes. a large number of mini-haloes that are dynamically heated and avoidH cooling, the vast majority of these fail the synchronised-pair cri-terion (we find 9 synchronised pairs that end up in the quasar-hosthalo at 𝑧 = synchronised multiplets . Totest whether these conditions actually occur in our overdense region,we evaluated the 𝐽 distribution sampled by dynamically heatedmini-haloes in our fiducial model, excluding haloes already flaggedas synchronised pairs. The results are shown in Fig. 9, where theblue histogram corresponds to the full distribution (including metalpolluted haloes), and the orange one to the pristine haloes only.Although the full distribution extends up to 𝐽 ∼ , most ofthe haloes in the high-LW flux tail are metal-rich, as such high 𝐽 values can be reached only near the very massive (hence highlystar-forming) haloes present in the region (see also Habouzit et al.2016a). When we remove the metal polluted haloes, the distributionshifts to lower LW fluxes, but with 38 haloes still exhibiting 𝐽 > ∼ 𝐽 values considered for efficientDCBH formation (Dijkstra et al. 2014), the non-negligible numberof such haloes, and most importantly the fact that one of them alsoexhibits a value close to a more conservative 𝐽 = , hint thatthe clustering of LW sources in overdense regions might representa viable alternative for DCBH formation to pair synchronisation. In this work, we combined a semi-analytic framework with N-bodysimulations, in order to assess the formation of DCBHs in overdenseregions where high-redshift quasar hosts are expected to form. Theframework exploited the dark-matter-only version of the simulationpresented in Lupi et al. (2019), where the particle mass resolutionis high enough too resolve the last stages of the mini-halo phase andthe transition to the atomic cooling halo phase.
MNRAS000
MNRAS000 , 1–15 (2021) A. Lupi et al.
We considered two different mechanisms proposed to be im-plicated in DCBH formation, the "synchronised pairs" scenario (Di-jkstra et al. 2008; Agarwal et al. 2012; Dijkstra et al. 2014; Visbalet al. 2014b; Habouzit et al. 2017; Chon et al. 2018; Regan et al.2019), where the time and spatial proximity between two ACHs (ofwhich one has already formed stars) guarantees a strong enough LWflux to suppress H formation (Shang et al. 2010; Wolcott-Greenet al. 2011; Latif et al. 2014; Wolcott-Green & Haiman 2019), andthe "dynamical heating" scenario (Yoshida et al. 2003; Inayoshiet al. 2018; Wise et al. 2019), where H formation is suppressed bycompressional heating caused by rapid halo growth.The goal of this paper is to identify halos that meet the follow-ing criteria: (i) metal-free, (ii) belonging to the quasar host mergerhistory, in which H formation (and cooling) has been significantlysuppressed in (iiia) the mini-halo phase, because of the combinedeffect of dynamical heating and a moderate LW flux), and/or (iiib) in the atomic-cooling phase, because of the presence of a nearbystar-forming ACH (a synchronised pair).This work shows that the very dense environments that lead tothe formation of the massive quasar hosts at high redshift lead to ahigh number of ACHs in a relatively small volume. A large fraction(as high as 85%) of these is dynamically heated during their mini-halo phase, helped by the high local levels of LW radiation dueto a strong clustering of nearby sources. The number density ofsynchronised pairs is also higher, by approximately an order ofmagnitude, than in average density regions of the Universe.Metal pollution decreases by a factor of about five the numberof halos eligible to host DCBH formation. In our fiducial case wefound one pristine synchronised pair of halos and ∼ ∼ LISA ) and forthird-generation terrestrial instruments, which will be able to detectmergers between seed BHs up to the redshift of their formation withhigh signal-to-noise.
ACKNOWLEDGEMENTS
AL acknowledges support from MIUR under the grant PRIN 2017-MB8AEZ. ZH acknowledges support from NSF grants 2006176 and1715661, and NASA grant NNX17AL82G. The analysis reportedin this work has been performed on the Horizon cluster, hosted bythe Institut d’Astrophysique de Paris.
DATA AVAILABILITY STATEMENT
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93Agarwal B., Khochfar S., Johnson J. L., Neistein E., Dalla Vecchia C., LivioM., 2012, MNRAS, 425, 2854Alexander T., Natarajan P., 2014, Science, 345, 1330Bañados E., et al., 2018, Nature, 553, 473Balberg S., Shapiro S. L., Inagaki S., 2002, ApJ, 568, 475Barkana R., 2004, MNRAS, 347, 59Barkana R., Loeb A., 2001, Phys. Rep., 349, 125Begelman M. C., 1979, MNRAS, 187, 237Begelman M. C., 2010, MNRAS, 402, 673Begelman M. C., Volonteri M., Rees M. J., 2006, MNRAS, 370, 289Begelman M. C., Rossi E. M., Armitage P. J., 2008, MNRAS, 387, 1649Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013a, ApJ, 762, 109Behroozi P. S., Wechsler R. H., Wu H.-Y., Busha M. T., Klypin A. A.,Primack J. R., 2013b, ApJ, 763, 18Boekholt T. C. N., Schleicher D. R. G., Fellhauer M., Klessen R. S., ReinosoB., Stutz A. M., Haemmerlé L., 2018, MNRAS, 476, 366Choi J.-H., Shlosman I., Begelman M. C., 2013, ApJ, 774, 149Chon S., Hirano S., Hosokawa T., Yoshida N., 2016, ApJ, 832, 134Chon S., Hosokawa T., Yoshida N., 2018, MNRAS, 475, 4104D’Amico G., Panci P., Lupi A., Bovino S., Silk J., 2018, MNRAS, 473, 328Das A., Schleicher D. R. G., Leigh N. W. C., Boekholt T. C. N., 2020, arXive-prints, p. arXiv:2012.01456Davies M. B., Miller M. C., Bellovary J. M., 2011, ApJ, 740, L42Dayal P., Ferrara A., 2018, Phys. Rep., 780, 1Decarli R., et al., 2018, ApJ, 854, 97Devecchi B., Volonteri M., 2009, ApJ, 694, 302Devecchi B., Volonteri M., Colpi M., Haardt F., 2010, MNRAS, 409, 1057Devecchi B., Volonteri M., Rossi E. M., Colpi M., Portegies Zwart S., 2012,MNRAS, 421, 1465Diemer B., 2018, ApJS, 239, 35Dijkstra M., Haiman Z., Mesinger A., Wyithe J. S. B., 2008, MNRAS, 391,1961Dijkstra M., Ferrara A., Mesinger A., 2014, MNRAS, 442, 2036Dotan C., Rossi E. M., Shaviv N. J., 2011, MNRAS, 417, 3035Fan X., et al., 2006, AJ, 131, 1203Fernandez R., Bryan G. L., Haiman Z., Li M., 2014, MNRAS, 439, 3798Fiacconi D., Rossi E. M., 2017, MNRAS, 464, 2259Gürkan M. A., Freitag M., Rasio F. A., 2004, ApJ, 604, 632Habouzit M., et al., 2016a, MNRAS, 456, 1901Habouzit M., Volonteri M., Latif M., Dubois Y., Peirani S., 2016b, MNRAS,463, 529Habouzit M., Volonteri M., Dubois Y., 2017, MNRAS, 468, 3935Haemmerlé L., Woods T. E., Klessen R. S., Heger A., Whalen D. J., 2018,MNRAS, 474, 2757Haiman Z., Loeb A., 2001, ApJ, 552, 459Haiman Z., Thoul A. A., Loeb A., 1996, ApJ, 464, 523Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H., 2003,ApJ, 591, 288Hirano S., Bromm V., 2017, MNRAS, 470, 898Hirano S., Hosokawa T., Yoshida N., Omukai K., Yorke H. W., 2015, MN-RAS, 448, 568Hirano S., Yoshida N., Sakurai Y., Fujii M. S., 2018, ApJ, 855, 17Hopkins P. F., 2015, MNRAS, 450, 53Hosokawa T., Omukai K., Yorke H. W., 2012, ApJ, 756, 93Hosokawa T., Yorke H. W., Inayoshi K., Omukai K., Yoshida N., 2013, ApJ,778, 178Inayoshi K., Visbal E., Kashiyama K., 2015, MNRAS, 453, 1692Inayoshi K., Haiman Z., Ostriker J. P., 2016, MNRAS, 459, 3738Inayoshi K., Li M., Haiman Z., 2018, MNRAS, 479, 4017Inayoshi K., Visbal E., Haiman Z., 2020, ARA&A, 58, 27Katz H., Sijacki D., Haehnelt M. G., 2015, MNRAS, 451, 2352Kauffmann G., Haehnelt M., 2000, MNRAS, 311, 576Kimura K., Hosokawa T., Sugimura K., 2020, arXiv e-prints, p.arXiv:2012.01452Koushiappas S. M., Bullock J. S., Dekel A., 2004, MNRAS, 354, 292MNRAS , 1–15 (2021) assive seed black holes in quasar hosts Kulkarni M., Visbal E., Bryan G. L., 2020, arXiv e-prints, p.arXiv:2010.04169Lacey C., Cole S., 1993, MNRAS, 262, 627Latif M. A., Schleicher D. R. G., Schmidt W., Niemeyer J. C., 2013, MNRAS,436, 2989Latif M. A., Bovino S., Van Borm C., Grassi T., Schleicher D. R. G., SpaansM., 2014, MNRAS, 443, 1979Latif M. A., Bovino S., Grassi T., Schleicher D. R. G., Spaans M., 2015,MNRAS, 446, 3163Latif M. A., Lupi A., Schleicher D. R. G., D’Amico G., Panci P., Bovino S.,2019, MNRAS, 485, 3352Lodato G., Natarajan P., 2006, MNRAS, 371, 1813Lupi A., Colpi M., Devecchi B., Galanti G., Volonteri M., 2014, MonthlyNotices of the Royal Astronomical Society, 442, 3616Lupi A., Haardt F., Dotti M., Fiacconi D., Mayer L., Madau P., 2016, MN-RAS, 456, 2993Lupi A., Volonteri M., Decarli R., Bovino S., Silk J., Bergeron J., 2019,MNRAS, 488, 4004Machacek M. E., Bryan G. L., Abel T., 2001, ApJ, 548, 509Madau P., Rees M. J., 2001, ApJ, 551, L27Madau P., Ferrara A., Rees M. J., 2001, ApJ, 555, 92Madau P., Haardt F., Dotti M., 2014, ApJ, 784, L38Mayer L., Kazantzidis S., Escala A., Callegari S., 2010, Nature, 466, 1082McKee C. F., Tan J. C., 2008, ApJ, 681, 771Miller M. C., Davies M. B., 2012, ApJ, 755, 81Mortlock D. J., et al., 2011, Nature, 474, 616Oh S. P., Haiman Z., 2002, ApJ, 569, 558Omukai K., 2001, ApJ, 546, 635Omukai K., Palla F., 2001, ApJ, 561, L55Omukai K., Schneider R., Haiman Z., 2008, ApJ, 686, 801Pacucci F., Volonteri M., Ferrara A., 2015, MNRAS, 452, 1922Peebles P. J. E., 1993, Principles of Physical Cosmology. Princeton Univer-sity PressPezzulli E., Valiante R., Schneider R., 2016, MNRAS, 458, 3047Planck Collaboration et al., 2016, A&A, 594, A13Pollack J., Spergel D. N., Steinhardt P. J., 2015, ApJ, 804, 131Portegies Zwart S. F., McMillan S. L. W., 2002, ApJ, 576, 899Regan J. A., Visbal E., Wise J. H., Haiman Z., Johansson P. H., Bryan G. L.,2017, Nature Astronomy, 1, 0075Regan J. A., Downes T. P., Volonteri M., Beckmann R., Lupi A., TrebitschM., Dubois Y., 2019, MNRAS,Regan J. A., Haiman Z., Wise J. H., O’Shea B. W., Norman M. L., 2020a,The Open Journal of Astrophysics, 3, E9Regan J. A., Wise J. H., Woods T. E., Downes T. P., O’Shea B. W., NormanM. L., 2020b, The Open Journal of Astrophysics, 3, 15Reinoso B., Schleicher D. R. G., Fellhauer M., Klessen R. S., BoekholtT. C. N., 2018, A&A, 614, A14Sakurai Y., Inayoshi K., Haiman Z., 2016, MNRAS, 461, 4496Sakurai Y., Haiman Z., Inayoshi K., 2020, MNRAS, 499, 5960Schaerer D., 2003, A&A, 397, 527Schauer A. T. P., et al., 2017, MNRAS, 467, 2288Schauer A. T. P., Glover S. C. O., Klessen R. S., Ceverino D., 2019, MonthlyNotices of the Royal Astronomical Society, 484, 3510Schlaufman K. C., Thompson I. B., Casey A. R., 2018, ApJ, 867, 98Shang C., Bryan G. L., Haiman Z., 2010, MNRAS, 402, 1249Shankar F., Weinberg D. H., Miralda-Escudé J., 2009, ApJ, 690, 20Shen X., Hopkins P. F., Faucher-Giguère C.-A., Alexander D. M., RichardsG. T., Ross N. P., Hickox R. C., 2020, MNRAS, 495, 3252Small T. A., Blandford R. D., 1992, MNRAS, 259, 725Soltan A., 1982, MNRAS, 200, 115Springel V., 2005, MNRAS, 364, 1105Stacy A., Greif T. H., Bromm V., 2012, MNRAS, 422, 290Tagawa H., Umemura M., Gouda N., Yano T., Yamai Y., 2015, MNRAS,451, 2174Tagawa H., Haiman Z., Kocsis B., 2020, ApJ, 892, 36Tanaka T., Haiman Z., 2009, ApJ, 696, 1798Tegmark M., Silk J., Rees M. J., Blanchard A., Abel T., Palla F., 1997, ApJ,474, 1 Trenti M., Stiavelli M., Michael Shull J., 2009, ApJ, 700, 1672Valiante R., Schneider R., Volonteri M., Omukai K., 2016, MNRAS, 457,3356Visbal E., Haiman Z., Bryan G. L., 2014a, MNRAS, 442, L100Visbal E., Haiman Z., Bryan G. L., 2014b, MNRAS, 445, 1056Volonteri M., Rees M. J., 2005, ApJ, 633, 624Volonteri M., Haardt F., Madau P., 2003, ApJ, 582, 559Volonteri M., Silk J., Dubus G., 2015, ApJ, 804, 148Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R., Yepes G.,2013, MNRAS, 433, 1230Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377Wise J. H., Regan J. A., O’Shea B. W., Norman M. L., Downes T. P., Xu H.,2019, Nature, 566, 85Wolcott-Green J., Haiman Z., 2019, MNRAS, 484, 2467Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS, 418, 838Wollenberg K. M. J., Glover S. C. O., Clark P. C., Klessen R. S., 2020,MNRAS, 494, 1871Woods T. E., Heger A., Whalen D. J., Haemmerlé L., Klessen R. S., 2017,ApJ, 842, L6Yoshida N., Abel T., Hernquist L., Sugiyama N., 2003, ApJ, 592, 645This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000