Evolution of High-redshift Quasar Hosts and Promotion of Massive Black Hole Seed Formation

We construct DM merger trees based on the extended Press–Schechter formalism (Press & Schechter 1974; Lacey & Cole 1993; Cole et al. 2000) using the GALFORM semianalytic algorithm summarized in Parkinson et al. (2008). Our sample consists of 10⁴ merger tree realizations for the DM halos that end up as high-z quasar hosts with M_h = 10¹² M_⊙ at z = 6. For each tree, we adopt a minimum DM halo mass of ${M}_{{\rm{h}},\min }={10}^{5}\,{M}_{\odot }$. Halos smaller than this threshold do not significantly impact the gas evolution, because the critical virial temperature above which gas collapse can be induced by coolant H₂ is ∼10³ K (see Haiman et al. 1996; Tegmark et al. 1997), corresponding to halo mass higher than ${M}_{{\rm{h}},\min }$ (see also Figure 1). Reflecting the rarity of quasar host galaxies, the progenitor halos form in highly biased regions with ≳5 mass variance (Mo & White 2002). Note that the fraction of all matter in such rare halos is ≲10⁻⁷.

Zoom In Zoom Out Reset image size

Due to the photodissociation of H₂ exposed to LW radiation, we also consider the local LW intensity J_LW (at h ν = 12.4 eV, hereafter in units of 10⁻²¹ erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹) in order to follow the gas evolution in a given progenitor halo. Along each merger tree, we calculate the cumulative J_LW from neighboring star-forming galaxies (hereafter source halos). Based on the model developed by Dijkstra et al. (2014), the basic equations and assumptions we adopt are summarized as below.

We consider a DM halo with mass M_h (gas + DM) at redshift z, which is supposed to be the main progenitor in a merger tree. The average number of source halos (within mass range m ± dm/2) that populate a surrounding spherical shell (at a physical distance r with thickness dr) is calculated by

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{d}^{2}{ \mathcal N }(m,r)}{{dmdr}}{dmdr}=4\pi {r}^{2}{dr}{\left(1+z\right)}^{3}\,\displaystyle \frac{{\rm{d}}{n}_{\mathrm{ST}}(m,z)}{{\rm{d}}m}\\ \,\times \,[1+\xi ({M}_{{\rm{h}}},m,z,r)]\,{\rm{d}}m,\end{array}\end{eqnarray} \tag{ 1 }$$

where dn_ST/dm is the mass function of source halos (Sheth et al. 2001) and ξ denotes the nonlinear bias function (Iliev et al. 2003), which gives the deviation (from random) probability of finding a halo with mass m at distance r from the main progenitor. We set the minimum source halo mass to be ${m}_{\mathrm{ac},z}\simeq 6\times {10}^{6}{M}_{\odot }{\left({T}_{\mathrm{vir}}/{10}^{4}\,{\rm{K}}\right)}^{3/2}{\left[\left(1+z\right)/31\right]}^{-3/2}$, where the halo virial temperature is just above the hydrogen atomic-cooling threshold of T_vir = 10⁴ K, where radiative cooling by Lyα emission leads to star formation. In our model, we do not consider the production of LW radiation background by star formation activity in less massive DM halos. The maximum mass of source halos is determined so that the LW intensity converges toward the higher-mass bins, namely, in terms of averaged flux, contributions from the ${m}_{\max }$ halos vanish owing to their low abundance. The value of ${m}_{\max }$ ranges from ∼10⁶ to ∼10¹⁰ M_⊙ and is larger at lower z.

Following Dijkstra et al. (2014), we compute the average LW radiation flux that irradiates the target halo. The time-averaged production rate of LW photons (per unit stellar mass) emitted from a surrounding source galaxy is approximated by

$$\begin{eqnarray}&&\left\langle {Q}_{\mathrm{LW}}(t)\right\rangle ={Q}_{0}{\left[1+\left({t}_{6}/4\right)\right]}^{-3/2}{{\rm{e}}}^{-{t}_{6}/300},\end{eqnarray} \tag{ 2 }$$

where ${Q}_{0}={10}^{47}\,{{\rm{s}}}^{-1}{M}_{\odot }^{-1}$ and t ( = t₆ Myr) is the time after a single starburst in the star-forming halo. Thus, the specific LW luminosity from the halo is calculated by

$$\begin{eqnarray}&&{L}_{\mathrm{LW}}({m}_{\star },t)=\displaystyle \frac{h\langle \nu \rangle }{{\rm{\Delta }}\nu }\left\langle {Q}_{\mathrm{LW}}(t)\right\rangle {f}_{\mathrm{esc},\mathrm{LW}}\left(\displaystyle \frac{{m}_{\star }}{{M}_{\odot }}\right),\end{eqnarray} \tag{ 3 }$$

where the mean frequency and frequency width of the LW band (11.2 ≤ h ν/eV ≤ 13.6) are set to 〈ν〉 = 12.4 eV/h and Δν =2.4 eV/h, respectively. The total stellar mass is calculated by m_⋆ = f_⋆(Ω_b/Ω_m)m, assuming the star formation efficiency to be f_⋆ = 0.05. The escape fraction of LW photons from the halo is assumed to be unity (f_esc,LW = 1). This value tends to be lower for atomic-cooling halos with m ≳ 10⁷ M_⊙. As a reference, Schauer et al. (2015) calculated the LW escape fraction for a single Population III star in an atomic-cooling halo with 1D simulations and found f_esc,LW ≃ 0.7. However, this is considered to be a lower bound because the escape fraction would be higher for 3D calculations through directions with lower optical depths, and besides, a higher star formation rate (SFR) is expected in our case (rather than a single massive star). We estimate the LW luminosity at one freefall time after the burst of star formation: ${t}_{\mathrm{sf}}=\sqrt{3\pi /(32G{{\rm{\Delta }}}_{\mathrm{vir}}\bar{\rho })}\simeq 18\,\mathrm{Myr}\,{[(1+z)/31]}^{-3/2}$, where Δ_vir ≃ 18 π². Using Equations (1)–(3), we obtain the mean LW radiation intensity in the target halo as

$$\begin{eqnarray}&&{J}_{\mathrm{LW}}({M}_{{\rm{h}}},z)={\int }_{{m}_{\mathrm{ac},z}}^{{m}_{\max }}{\int }_{{r}_{\min }}^{{r}_{\max }}\displaystyle \frac{{d}^{2}{ \mathcal N }(m,r)}{{dmdr}}\cdot \displaystyle \frac{{L}_{\mathrm{LW}}}{16{\pi }^{2}{r}^{2}}\,{dmdr},\end{eqnarray} \tag{ 4 }$$

where ${r}_{\min }$ and ${r}_{\max }$ are the minimum and maximum distance, respectively, of the source halo from the target halo. In the absence of metal pollution, ${r}_{\min }$ can be safely set by adding the virial radii of the target and source halos. However, metal enrichment of the main progenitor is a main obstacle in the formation scenario of massive seed BHs, because efficient metal-line cooling (and possibly dust thermal emission) will likely lead to gas fragmentation during its gravitational collapse and thus suppress massive star formation. Generically, there are two types of enrichment processes: (1) genetic enrichment due to past star formation episodes in the progenitors, and (2) environmental enrichment owing to metal bubbles created by supernova (SN) explosions in nearby galaxies. In our model, we consider the environmental enrichment process by adopting the minimum distance to source halos as ${r}_{\min \ }=\max \{{r}_{\mathrm{vir}}({M}_{{\rm{h}}})+{r}_{\mathrm{vir}}(m),{r}_{{\rm{s}}}(m)\}$, where r_s is the size of the metal-polluted region surrounding the source halo

$$\begin{eqnarray}&&{r}_{{\rm{s}}}(m,t)={\left(\displaystyle \frac{{E}_{\mathrm{SN}}{m}_{\star }}{{m}_{0}{\rho }_{{\rm{s}}}}\right)}^{1/5}{t}^{2/5},\end{eqnarray} \tag{ 5 }$$

where m₀ = 100 M_⊙ is the stellar mass budget required to form an SN progenitor and ${E}_{\mathrm{SN}}={10}^{51}\,\mathrm{erg}$ is the explosion energy of an SN. The density ρ_s of gas surrounding the wind is considered to be ${\rm{\Delta }}{\bar{\rho }}_{{\rm{b}}}$, where ${\bar{\rho }}_{{\rm{b}}}$ is the intergalactic medium (IGM) baryon density, and Δ = 60 corresponding to the typical baryonic overdensity of halos at their virial radius for an NFW profile (Dijkstra et al. 2014). Similar to the production of LW radiation, we estimate the size of metal-enriched bubbles at t_sf. We note that metal enrichment through in situ star formation should be subdominant because intense LW radiation suppresses star formation in low-mass progenitors (see Section 4).

On the other hand, the maximum distance in the integration is given by ${r}_{\max }=\left({\lambda }_{\mathrm{LW},1}-{\lambda }_{\beta }\right)c/\left[{\lambda }_{\beta }H(z)\right]$, where the λ_LW,1 = 1110 Å and λ_β are wavelengths of the lowest LW energy and Lyβ line, respectively (see Haiman et al. 1997). We consider the redshift effect by cosmic expansion, where $H(z)={H}_{0}\left[{{\rm{\Omega }}}_{{\rm{m}}}\right(1+z{\left){}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}\right]}^{1/2}$ is the Hubble constant at redshift z and c is the light speed. LW photons emitted at $r\gt {r}_{\max }$ are redshifted into one of the Lyman series resonances and are converted into low-energy photons before reaching the target halo. The ${r}_{\max }$ is thus set as an absorbing screen, i.e., we exclude the contribution of J_LW from halos located at $r\gt {r}_{\max }$.

The main progenitor halo experiences vigorous halo mergers in the high-z universe. Successive merger events, in particular major mergers, inject energy into the gas in the parent halo. At early phase, energy loss through radiative cooling is inefficient, i.e., the cooling timescale is comparable to or longer than the Hubble timescale. Gas is heated through shock formation at the halo virial radius in an adiabatic manner. Subsequently, the energy is transported into the halo interior, leading to gas virialization with a nearly constant temperature profile (T_gas ∼ T_vir) across all radii (Wise & Abel 2007). Assuming that the virial equilibrium state is reached after a merger event, the virial theorem applies to the gas in the post-merger halo, where the internal and kinetic (turbulence) energy is balanced with the gravitational energy as

$$\begin{eqnarray}&&{e}_{\mathrm{tot}}={e}_{\mathrm{th}}+{e}_{{\rm{k}}}+{{\rm{\Phi }}}_{{R}_{\mathrm{vir}}}=\displaystyle \frac{1}{2}{{\rm{\Phi }}}_{{R}_{\mathrm{vir}}},\end{eqnarray} \tag{ 6 }$$

where e_tot, e_th, and e_k are the total, thermal, and kinetic energy per unit mass, respectively, and ${{\rm{\Phi }}}_{{R}_{\mathrm{vir}}}$ is the gravitational energy at the virial radius. In this work, we adopt the NFW potential for DM halos given by

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{R}_{\mathrm{vir}}}=-\displaystyle \frac{2{k}_{{\rm{B}}}{T}_{\mathrm{vir}}}{\mu {m}_{{\rm{p}}}}\cdot \displaystyle \frac{\mathrm{ln}(1+{c}_{\mathrm{vir}})}{\mathrm{ln}(1+{c}_{\mathrm{vir}})-{c}_{\mathrm{vir}}/(1+{c}_{\mathrm{vir}})},\end{eqnarray} \tag{ 7 }$$

where T_vir is the halo virial temperature, the concentration parameter of the DM density profile ${c}_{\mathrm{vir}}=1.9\,{({M}_{{\rm{h}}}/{10}^{7}\,{M}_{\odot })}^{-0.13}{[(1+z)/31]}^{-1}$ (Bullock et al. 2001), k_B is the Boltzmann constant, μ = 1.22 is the mean molecular weight, and m_p is the proton mass. Therefore, the total energy change owing to the halo evolution is given by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{mrg}}=-\displaystyle \frac{1}{2}{{\rm{\Phi }}}_{{R}_{\mathrm{vir}}}\left(\displaystyle \frac{2}{3}\displaystyle \frac{\dot{{M}_{{\rm{h}}}}}{{M}_{{\rm{h}}}}-\displaystyle \frac{1}{1+z}\displaystyle \frac{{dz}}{{dt}}\right),\end{eqnarray} \tag{ 8 }$$

where the first term on the right-hand side denotes the energy change associated with mass growth and the second term represents the cosmic expansion effect. In the generally turbulent virialized gas, the kinetic-to-thermal energy ratio is equal to 1 around the virial radius and decreases to 1/3 at the center (see Wise & Abel 2007). Adopting this branching ratio of the total injected energy, the thermal and kinetic heating rates associated with mergers are given by Γ_mrg,th = 3Γ_mrg/4 and Γ_mrg,kin = Γ_mrg/4, respectively. Combining Equations (6)–(8), the gas temperature follows the halo virial temperature as

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{T}}_{\mathrm{gas}}}{{\dot{T}}_{\mathrm{vir}}}=\displaystyle \frac{1}{2}\cdot \displaystyle \frac{\mathrm{ln}(1+{c}_{\mathrm{vir}})}{\mathrm{ln}(1+{c}_{\mathrm{vir}})-{c}_{\mathrm{vir}}/(1+{c}_{\mathrm{vir}})}.\end{eqnarray} \tag{ 9 }$$

This ratio is close to unity for a wide range of (M_h, z) halos of interest, e.g., ${\dot{T}}_{\mathrm{gas}}/{\dot{T}}_{\mathrm{vir}}\simeq 1.3$ and 0.81 for c_vir = 2 and 10. Note that our method is different from that adopted in previous papers (e.g., Yoshida et al. 2003; Lupi et al. 2021), where T_gas = T_vir is imposed. The treatment allows us to precisely calculate the radiative cooling rates and chemical reaction coefficients, which sensitively depend on the gas temperature.

The kinetic energy injected through mergers is stored in the halo as turbulence. During the viliarization process, turbulence plays an important role in massive star formation (e.g., McKee & Tan 2002). Namely, turbulence acts as a source of pressure, which stabilizes the gas against its self-gravity and delays the collapse until the cloud becomes massive enough to overcome the turbulent pressure. In addition to turbulence, the baryonic streaming motion relative to the DM produced in the epoch of cosmic recombination at z_rec ≃ 1100 also significantly delays gas collapse and star formation in protogalaxies. The streaming velocity is found to follow a Maxwell−Boltzmann distribution with the rms speed of σ = 30 km s⁻¹ at z = z_rec and decays as ${\tilde{v}}_{\mathrm{bsm}}={v}_{\mathrm{bsm}}(1+z)/(1+{z}_{\mathrm{rec}})$ (Tseliakhovich & Hirata 2010). We note that the volume fraction of the universe with streaming velocities of v_bsm ≥ A σ is estimated as ≃0.4, 8 × 10⁻³, and 5.9 × 10⁻⁶ for A = 1, 2, and 3, respectively.

Considering both the three-dimensional turbulence and coherent baryonic streaming velocity, we approximate the effective pressure by kinetic motion of gas as

$$\begin{eqnarray}&&{P}_{\mathrm{tur}}\approx \displaystyle \frac{1}{3}\rho {v}_{\mathrm{tur}}^{2}+\rho {\left[{\alpha }_{0}{\tilde{v}}_{\mathrm{bsm}}(z)\right]}^{2},\end{eqnarray} \tag{ 10 }$$

where ${v}_{\mathrm{tur}}^{2}=2\int {{\rm{\Gamma }}}_{\mathrm{mrg},\mathrm{kin}}{dt}$ is the kinetic specific energy accumulated through successive mergers and the coefficient of 1/3 is required to estimate the pressure due to isotropic turbulence (Chandrasekhar 1951a, 1951b). With pressure support from turbulence, gas collapse is delayed to different extents, with varying strengths of the streaming motion. In this work, we adopt α₀ = 4.7 in our fiducial model, in order to match the delay of collapse obtained from cosmological simulations (Hirano et al. 2018). The total gas pressure is therefore defined by P_tot = P_gas + P_tur.

With the energy injection processes defined above, in this section we describe our model for calculating the density evolution of a gas cloud concentrated in a DM halo that grows through successive merger episodes. There are three characteristic stages of the evolution: (1) initial adiabatic phase, (2) transition to isothermal gas due to radiative cooling, and (3) gravitationally collapsing phase in a runaway fashion. We model the gas dynamics in these stages based on a one-zone model (e.g., Omukai 2001), which is often used to follow the physical quantities at the center of a gravitationally collapsing cloud with a self-similar density profile ρ_gas ∝ r⁻². However, this profile does not apply to gas in hydrostatic equilibrium before the onset of gravitational collapse. Therefore, we construct a new method to model the three characteristic stages in a physically motivated way.

2.5.1. Adiabatic Stage

In the early stage, since the gas density is not high enough for radiative cooling to operate through collisionally excited transitions, the gas is adiabatically compressed in the DM halo as the underlying DM gravitational potential evolves. In the DM assembly history through mass accretion, the entropy profile K(r) of the adiabatic gas is characterized by a power-law outer profile of $K{(r)={K}_{\mathrm{vir}}(r/{R}_{\mathrm{vir}})}^{1.1}$ and a constant core with K₀ ≃ 0.1K_vir, where ${K}_{\mathrm{vir}}={k}_{{\rm{B}}}{T}_{\mathrm{vir}}/\left[(\mu {m}_{{\rm{p}}}){\bar{\rho }}_{{\rm{b}}}^{2/3}\right]$ is the gas entropy at the virial radius (Voit et al. 2003, 2005). This self-similar entropy profile is also found to be established inside high-z protogalaxies formed in DM halos more massive than 3 × 10⁶ M_⊙ at z = 10, while the core entropy for less massive halos is maintained at the IGM entropy when gas decouples from the cosmic microwave background (CMB; see more details in Visbal et al. 2014a). Motivated by both numerical simulations and galaxy cluster observations, we approximate the entropy profile as

$$\begin{eqnarray}&&K(r)\simeq {K}_{\mathrm{vir}}\left(\displaystyle \frac{r}{{R}_{\mathrm{vir}}}\right)+{K}_{0},\end{eqnarray} \tag{ 11 }$$

where ${K}_{0}=\max (0.1{K}_{\mathrm{vir}},{K}_{\mathrm{IGM}})$. Using the entropy profile and the equation of state given by ${P}_{\mathrm{gas}}=K(r){\rho }_{\mathrm{gas}}^{\gamma }$, where γ = 5/3, we calculate the density profile by solving the hydrostatic equation (the so-called Lane−Emden equation) for the cloud embedded in the DM potential:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{d}{{dr}}\left[\displaystyle \frac{{r}^{2}}{{\rho }_{\mathrm{gas}}}\displaystyle \frac{d(K{\rho }_{\mathrm{gas}}^{\gamma }+{P}_{\mathrm{tur}})}{{dr}}\right]=-4\pi G\left({\rho }_{\mathrm{gas}}+{\rho }_{\mathrm{DM}}\right).\end{eqnarray} \tag{ 12 }$$

Throughout this paper, we adopt the NFW density profile of DM halos of all masses characterized by a simple analytical form of

$$\begin{eqnarray}&&{\rho }_{\mathrm{DM}}(r)={\rho }_{{\rm{m}}}(z)\displaystyle \frac{{\delta }_{0}}{\left({c}_{\mathrm{vir}}r/{R}_{\mathrm{vir}}\right){\left(1+{c}_{\mathrm{vir}}r/{R}_{\mathrm{vir}}\right)}^{2}},\end{eqnarray} \tag{ 13 }$$

where ρ_m(z) is the mean matter density and

$$\begin{eqnarray}&&{\delta }_{0}=\displaystyle \frac{200}{3}\displaystyle \frac{{c}_{\mathrm{vir}}^{3}}{\mathrm{ln}(1+{c}_{\mathrm{vir}})-{c}_{\mathrm{vir}}/(1+{c}_{\mathrm{vir}})}\end{eqnarray} \tag{ 14 }$$

is the characteristic overdensity within the halo virial radius (Navarro et al. 1997).

We integrate this hydrostatic equation with respect to ρ_gas(r), imposing the regularity conditions at the center, i.e., ρ_gas = ρ₀ and d ρ_gas/dr = 0 at r = 0. Since the solution for adiabatic gas generally has the radius r₀ where ρ_gas(r₀) = 0, we determine the central density ρ₀ so that the enclosed gas mass at r ≤ r₀ satisfies M_gas = f_b M_h, where f_b = Ω_b/Ω_m is the baryonic fraction.

2.5.2. Isothermal Stage

As the gas temperature increases owing to gravitational compression and merger heating, radiative cooling processes begin to operate in the cloud and the adiabatic assumption no longer applies. When the radiative cooling timescale is shorter than the heating timescale, we solve the hydrostatic equation for the density profile assuming an isothermal equation of state:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{d}{{dr}}\left[{r}^{2}{c}_{\mathrm{eff}}^{2}\displaystyle \frac{d\mathrm{ln}{\rho }_{\mathrm{gas}}}{{dr}}\right]=-4\pi G\left({\rho }_{\mathrm{gas}}+{\rho }_{\mathrm{DM}}\right),\end{eqnarray} \tag{ 15 }$$

where ${c}_{\mathrm{eff}}\equiv \sqrt{{c}_{{\rm{s}}}^{2}+{v}_{\mathrm{tur}}^{2}/3+{\left({\alpha }_{0}{\tilde{v}}_{\mathrm{bsm}}\right)}^{2}}$ is the effective sound speed developed from the isothermal sound speed ${c}_{{\rm{s}}}\equiv \sqrt{{k}_{{\rm{B}}}{T}_{\mathrm{gas}}/(\mu {m}_{{\rm{p}}})}$. The solution of the isothermal Lane−Emden equation with the regularity condition does not have the radius where the density becomes zero but connects to the external medium with a density of ρ_ext = f_b ρ_DM. The central density is determined so that ρ_gas = ρ_ext at the virial radius.

From the analogy of the Bonnor−Ebert sphere, the isothermal gas cloud embedded in a DM halo potential has a critical mass for the onset of its gravitational collapse. Practically, for a given T_gas and ρ_DM(r), we construct the density profile with different values of the gas central density ρ₀ and thus obtain ρ_gas(R_vir) as a function of ρ₀. Since this function has a local maximum value and the value decreases with increasing halo mass, a hydrostatic equilibrium solution where ρ_gas(R_vir) = ρ_ext no longer exists for M_h ≥ M_h,crit (see the Appendix). In this case, the gas evolution is described by the freefall stage below.

2.5.3. Freefall Stage

Once the gas cloud becomes gravitationally unstable, the evolution of the gas density profile is well described by the Penston–Larson self-similar solution (Larson 1969; Penston 1969), which has a density profile with a flat core of the Jeans scale and an envelope with a power-law density distribution ρ_gas(r) ∝ r⁻². The central density increases over the freefall timescale as

$$\begin{eqnarray}&&\displaystyle \frac{d{\rho }_{\mathrm{gas}}}{{dt}}=\displaystyle \frac{{\rho }_{\mathrm{gas}}}{{t}_{\mathrm{ff}}},\end{eqnarray} \tag{ 16 }$$

where the freefall timescale is calculated with

$$\begin{eqnarray}&&{t}_{\mathrm{ff}}\equiv \sqrt{\displaystyle \frac{3\pi }{32G\left({\rho }_{\mathrm{gas}}+\langle {\rho }_{\mathrm{DM}}\rangle \right)}},\end{eqnarray} \tag{ 17 }$$

where 〈ρ_DM〉 = ρ_m(z)δ₀ represents the averaged DM density. ¹

In the collapsing stage, compressional heating by the self-gravitating gas is taken into account and the rate is given by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{comp}}\equiv \displaystyle \frac{{P}_{\mathrm{gas}}+{P}_{\mathrm{tur}}}{{\rho }_{\mathrm{gas}}^{2}}\cdot \displaystyle \frac{d{\rho }_{\mathrm{gas}}}{{dt}}=\displaystyle \frac{{c}_{\mathrm{eff}}^{2}}{{t}_{\mathrm{ff}}}.\end{eqnarray} \tag{ 18 }$$

We note that the compressional heating rate is enhanced by turbulent pressure through the effective sound speed.

We consider the evolution of thermal and kinetic energy of the gas by solving the two energy equations:

$$\begin{eqnarray}&&\displaystyle \frac{{{de}}_{\mathrm{th}}}{{dt}}={{\rm{\Gamma }}}_{\mathrm{mrg},\mathrm{th}}+{{\rm{\Gamma }}}_{\mathrm{comp}}-{{ \mathcal L }}_{\mathrm{chem}}-{{ \mathcal L }}_{\mathrm{rad}},\end{eqnarray} \tag{ 19 }$$

where ${{ \mathcal L }}_{\mathrm{chem}}$ is the cooling/heating rate associated with chemical reactions and ${{ \mathcal L }}_{\mathrm{rad}}$ is the radiative cooling rate (note that all the rates are in units of erg s⁻¹ g⁻¹). While the compressional heating rate is included only in the collapse stage, the other effects are taken into account to calculate the gas temperature over the three evolutionary stages. The cooling term includes radiative cooling by H, He, He⁺, and He⁺⁺ (Glover & Jappsen 2007); H₂ (Glover & Abel 2008; Glover 2015a, 2015b); and cooling/heating associated with chemical reactions.

We solve the chemical reactions of primordial gas among the following nine species: H, H₂, e⁻, H⁺, ${{\rm{H}}}_{2}^{+}$, H⁻, He, He⁺, and He⁺⁺. In Table 1, we show the 35 reaction rate coefficients adopted in this work. In terms of photodissociation of H₂, H⁻, and ${{\rm{H}}}_{2}^{+}$ by external radiation emitted from nearby star-forming galaxies, the reaction rate is calculated by assuming the radiation spectral energy distribution (SED) to be a blackbody spectrum with T_rad = 2 × 10⁴ K. The SED model approximates more realistic spectra of observed metal-poor star-forming galaxies (Inoue 2011). The dissociation rates of H⁻ and ${{\rm{H}}}_{2}^{+}$ are calculated by a convolution with the cross section of the ith chemical species (i= H⁻ and ${{\rm{H}}}_{2}^{+}$) as

$$\begin{eqnarray}&&{k}_{{\rm{i}},\mathrm{pd}}={\int }_{0}^{\infty }\displaystyle \frac{4\pi J(\nu )}{h\nu }{\sigma }_{{\rm{i}}}(\nu )d\nu .\end{eqnarray} \tag{ 20 }$$

The cross sections we adopt are from references listed in Table 1.

Table 1. Chemical Reactions

Number	Reaction	Reference
	H Collisional Reactions
1	H + e− → H+ + 2e−	1
2	H+ + e− → H + γ	2*
3	H + e− → H− + γ	3
4	H− + H → H2 + e	4
5	H + H+ → ${{\rm{H}}}_{2}^{+}$ + γ	5
6	${{\rm{H}}}_{2}^{+}$+ H → H2 + H+	6
7	H2 + H → 3H	7
8	H2 + H+ → ${{\rm{H}}}_{2}^{+}$ + H	8
9	H2 + e− → 2H + e−	9
10	H− + e− → H + 2e−	10
11	H− + H+ → 2H	11
12	H− + H+ → ${{\rm{H}}}_{2}^{+}$ + e−	12
13	${{\rm{H}}}_{2}^{+}$ + e− → 2H	13
14	${{\rm{H}}}_{2}^{+}$ + H− → H2 + H	14
15	3H → H2 + H	15
16	2H + H2 → 2H2	16
17	2H2 → 2H + H2	17
18	H− + H → 2H + e−	18
19	H− + ${{\rm{H}}}_{2}^{+}$ → 3H	19
20	H2 + e− → H− + H	20
	Photodissociation and Detachment Reactions
21	H2 + γ → 2H	21
22	H− + γ → H + e−	22
23	${{\rm{H}}}_{2}^{+}$ + γ → H + H+	23
	He Reactions
24	He + e− → He+ + 2e−	24
25	He+ + e− → He + γ	25
26	He+ + e− → He++ + 2e−	26
27	He++ + e− → He+ + H+ +γ	27
28	H2 + He → 2H + He	28
29	H2 + He+ → He + H + H+	29
30	H2 + He+ → ${{\rm{H}}}_{2}^{+}$ + He	30
31	He+ + H → He + H+	31
32	He + H+ → He+ + H	32
33	He+ + H− → He + H	33
34	He + H− → He + H + e−	34
35	2H + He → H2 + He	35

References. (1) Abel et al. 1997. (2) Ferland et al. 1992, Case B. (3) McLaughlin et al. 2017. (4) Kreckel et al. 2010. (5) Coppola et al. 2011. (6) Karpas et al. 1979. (7) Mac Low & Shull 1986; Lepp & Shull 1983. (8) Savin et al. 2004; Coppola et al. 2011. (9) Trevisan & Tennyson 2002. (10) Janev et al. 1987. (11) Croft et al. 1999. (12) Poulaert et al. 1978. (13) Schneider et al.1994. (14) Dalgarno & Lepp 1987. (15) Abel et al. 2002; Orel 1987. (16) Jacobs et al. 1967. (17) Martin et al. 1998; Shapiro & Kang 1987. (18) Janev et al. 1987. (19) Dalgarno & Lepp 1987. (20) Schulz & Asundi 1967. (21) Wolcott-Green & Haiman 2011. (22) McLaughlin et al. 2017. (23) Stancil 1994. (24) Janev et al. 1987. (25) Hummer & Storey 1998. (26) Janev et al. 1987. (27) Ferland et al. 1992. (28) Dove et al. 1987. (29) Barlow 1984. (30) Barlow 1984. (31) Zygelman et al. 1989. (32) Kimura et al. 1993. (33) Peart & Hayton 1994. (34) Huq et al. 1982. (35) Glover & Abel 2008.

Download table as:  ASCIITypeset image

In Figure 1, we show the evolution of the main progenitors, i.e., the most massive halos at each epoch, for all the 10⁴ merger trees that grow to M_h = 10¹² M_⊙ at z = 6. In such overdense regions of the universe, the DM halo mass increases via rapid mergers. The median halo mass (dashed curve) reaches M_h ≃ 8 × 10¹⁰ M_⊙, 6 × 10⁸ M_⊙, 2 × 10⁷ M_⊙, and 8 × 10⁵ M_⊙ at z = 10, 20, 30, and 40, respectively, and the virial temperature exceeds the atomic-cooling threshold of T_vir ≃ 10⁴ K at z ≃ 34. Therefore, the gas cloud concentrated in the massive halo collapses at an epoch earlier than when typical first galaxies would form in atomic-cooling halos (M_h ≃ 10⁷ M_⊙ at z ≃ 10; see Bromm & Yoshida 2011), which are usually considered to be massive-seed-forming sites in most previous studies (e.g., Dijkstra et al. 2014). For illustration purposes, we highlight three merger trees: the blue (id 1, a less massive tree), orange (id 2, a tree comparable to the median evolution), and green curves (id 3, a more massive tree). In the following sections, we focus our analysis on these three representative cases.

Following the method laid out in Section 2, in Figure 2 we present the redshift evolution of J_LW for the three representative trees and the median track. For all the cases, the LW background intensity gradually increases from higher redshifts, peaks at the intermediate redshifts, and decreases toward lower redshifts. This redshift dependence reflects the nature of the nonlinear bias function that boosts the abundance of halo pairs with comparable masses (Scannapieco & Barkana 2002). Namely, when the mass of the main progenitor is close to the atomic-cooling halo mass (m_ac,z ∼ 10⁷ M_⊙), a large number of source halos form nearby owing to the halo clustering effect, and thus the LW intensity is maximized. As the main progenitor grows, its mass difference from m_ac,z is larger, and thus the clustering effect of atomic-cooling sources becomes weaker so that their spatial distribution is approximated to be uniform (i.e., ξ ≪ 1). As a result, the LW intensity is dominated by the contribution from a large number of atomic-cooling source halos within the absorbing screen ($r\lesssim {r}_{\max }$) and begins to decline owing to the cosmic dilution effect at lower redshifts. For rapidly growing progenitor halos exceeding m_ac,z earlier, the LW intensity quickly rises at higher redshifts and the peak values become higher owing to stronger clustering at earlier epochs. Namely, the peak values of LW intensity in the overdense regions are J_LW ≃ 60 (id 1), J_LW ≃ 400 (id 2), J_LW ≃ 600 (median), and J_LW ≃ 6 × 10³ (id 3), which are significantly higher than the level of LW intensity irradiating typical atomic-cooling halos that are expected to form massive BH seeds (see Dijkstra et al. 2008; Agarwal et al. 2012; Johnson et al. 2013).

Zoom In Zoom Out Reset image size

In our semianalytical approach, we model metal pollution of the progenitor halos due to SN explosions that occur in source halos. Although we treat this effect by replacing the minimum distance between the target and source halos with r_s, there is no information on the time-dependent spatial distributions of DM halos in our framework. To examine the impact of the model assumptions, in Figure 2 we also show the case where the size of the metal-polluted bubbles (r_s) is doubled and the corresponding t_sf is comparable to the Hubble time at the redshift, or equivalent to setting Δ = 1 with the fiducial value of t_sf. In this case, the LW intensity is overall reduced at higher redshifts, indicating a significant contribution from nearby source halos with ≳m_ac,z to the LW radiation background. We note that our treatment simply removes the contribution from source halos within distances of r_s but does not address how likely it is for the main progenitor to be affected by environmental metal enrichment. Our argument nevertheless provides a conservative estimate of J_LW if the efficiency of environmental metal enrichment is low. As discussed in Section 4, the efficiency should be negligibly low because metal-polluted bubbles rarely penetrate the interior of the target halo (Chiaki et al. 2018).

In Figure 3, we present the histograms of the LW background intensity that irradiates the main progenitor halos for the 10⁴ trees at different redshifts. For the whole sample of the target halos in highly biased regions, the histogram resembles a probability distribution function (pdf) of J_LW, with the bar height in each bin (${\rm{\Delta }}\mathrm{log}{J}_{\mathrm{LW}}=0.3$) representing the number fraction of halos irradiated within $\mathrm{log}{J}_{\mathrm{LW}}\to \mathrm{log}{J}_{\mathrm{LW}}+{\rm{\Delta }}\mathrm{log}{J}_{\mathrm{LW}}$. From higher redshifts down to z ≃ 30, the mean value of J_LW in the pdf increases owing to a large number of clustered source halos with ≳m_ac,z and the J_LW distribution peaks around ≃270 at z ≃ 25. Toward lower redshifts, the target halo mass becomes higher than the typical mass of source halos. Therefore, the abundance of sources is hardly boosted by the clustering effect (Iliev et al. 2003). Moreover, the LW intensity is diluted by the cosmic expansion, lowering the mean value. While the dispersion of the pdf is larger at higher redshifts, reflecting the diversity of the progenitor mass, the pdf peaks at J_LW ≃ 60 by z = 10 when all the 10⁴ trees converge to the high-z quasar host. We note that our model does not consider LW radiation produced from DM minihalos with m < m_ac,z , where H₂ is the only coolant to induce star formation. However, strong LW background radiation in the overdense region likely suppresses its formation. Therefore, the histogram shown in Figure 3 counts the lower bound of the LW background intensity.

Zoom In Zoom Out Reset image size

In this section we focus our analysis on the gas properties in the main progenitors along the three representative merger trees. In Figure 4, we show the evolution of gas density (left panels) and temperature (right panels) at the central core as a function of redshift. In order to examine the impact of baryonic streaming motion, for each merger tree we assume two different v_bsm values, i.e., v_bsm = 1σ (top panels), and 2σ (bottom panels). Each curve corresponds to the representative case highlighted in Figure 1. Along with the three evolutionary tracks, we denote the epochs when the DM halo mass exceeds M_h = 10⁶ M_⊙, 10⁷ M_⊙, and 3 × 10⁷ M_⊙ in the left panels, and when the LW background intensity first crosses J_LW = 1, 10, 10², and 10³ in the right panels. In the following paragraphs, we first describe the gas properties with v_bsm = 1σ and then discuss the impact of the baryonic streaming motion on gas evolution in cases with v_bsm = 2σ.

Zoom In Zoom Out Reset image size

For the lowest-mass case (blue curve, tree id 1), the gas density gradually increases with the halo mass in the early stage (z > 30), where the gas cloud is supported by thermal and turbulent pressure against its self-gravity and DM gravitational force. After the halo mass reaches ≃10⁶ M_⊙, the cloud becomes gravitationally unstable owing to its low temperature and collapses over one freefall timescale at z ≃ 28. The gas temperature remains at T ≲ 10³ K owing to H₂ cooling, under a modest level of LW intensity (J_LW ∼ 1) at z > 35. In addition to LW radiation, the gas is heated by four major merger events around z ≃ 31–34, but the dynamical heating rate does not overcome the H₂ cooling rate in this case.

For the intermediate-mass case (orange curve, tree id 2), the evolution begins from a redshift higher than in the previous case. In this case, the gas temperature is substantially higher as a result of the combination of merger heating and intense LW irradiation with J_LW ≳ 1 in the early stage. As several episodes of halo mergers increase the halo mass to ∼10⁷ M_⊙ by z ≃ 30 (the corresponding halo virial temperature is T_vir ≃ 10⁴ K), the gas temperature reaches T ≃ 10⁴ K, where the atomic cooling via Lyα emission begins to operate. Although the LW intensity reaches J_LW ≳ 100 before the cloud gravitationally collapses, the level of LW intensity is not strong enough to suppress H₂ formation in the dense region (≳10² cm⁻³), where H₂ reforms owing to its self-shielding effect. As a result of efficient H₂ cooling, the gas temperature drops down to T ≃ 10³ K in the collapsing stage.

For the highest-mass case (green curve, tree id 3), the gas temperature quickly rises to T ≃ 10⁴ K owing to frequent mergers. Owing to the clustering effect of the massive parent halo, the LW intensity reaches J_LW ≳ 10³ at z ≃ 47, prominently higher than those seen in the less massive cases. Although the H₂ self-shielding becomes more effective as the central density increases up to ≳10⁴ cm⁻³, the gas collapses, keeping a nearly constant temperature of T ≃ 8000 K. Inside the dense and warm region, H₂ is collisionally dissociated and its radiative cooling does not alter the thermal evolution.

In cases where v_bsm = 2σ, the gas property evolution is shown in the bottom panels of Figure 4. Overall, the collapse of gas clouds is delayed owing to kinetic energy injection to the gas concentrated at the halo center. When the cloud begins to collapse, the corresponding halo masses reach M_h ≃ (3.5, 4.2, 5.9) × 10⁷ M_⊙. For comparison, the collapse halo masses are M_h ≃ (0.24, 2.1, 2.2) × 10⁷ M_⊙ for v_bsm = 1σ. The delay effect is more remarkable for the lower-mass cases because the halo circular velocity is lower than the effective sound speed boosted by injection of turbulence and streaming motion. As the gas collapse proceeds, H₂ forms efficiently in the modest J_LW environment, and eventually its cooling reduces the gas temperature in the low- and intermediate-mass cases.

As noted in Section 3.2 and Figure 4, depending on the main cooling processes inducing star formation, the evolutionary tracks of the gas clouds embedded in the main progenitors of high-z quasar hosts are classified into three cases: (i) H₂ cooling, (ii) initial H Lyα cooling followed by H₂ cooling after a short isothermal collapse, and (iii) H Lyα cooling when temperature is kept above 8000 K by compression along a wide density range. In Figure 5, we present the number count of merger trees for the three types with different baryonic streaming velocities, denoted as (i) H₂, (ii) H−H₂, and (iii) H−H. Without the streaming velocity, 74% of the trees experience gas collapse via H₂ cooling, while the rest (26%) form atomically cooling gas clouds (cases H−H₂ and H−H). With nonzero streaming motion (v_bsm ≠ 0), nearly all cases enter the atomic-cooling stage because the halo mass reaches m_ac,z via mergers owing to the significant delay effect. As the streaming velocity increases, the gas mass becomes higher at the onset of gravitational collapse, and thus the compressional heating rate during the collapse stage is higher owing to the accumulation of kinetic energy. Therefore, the number of trees where gas isothermally collapses with T ≃ 8000 K (case H−H) increases monotonically from 14% to 27% with increasing streaming velocity from v_bsm = 1σ to 3σ.

Zoom In Zoom Out Reset image size

In Figure 6, we show the distributions of the halo virial temperature (top panels) and LW background intensity (middle panels) for the three types of gas collapse. For each case, the values of T_vir and J_LW are measured at the epoch when the gas cloud first enters its unstable stage. In contrast to cases with v_bsm = 0, where gas collapse is led by H₂ cooling in less massive halos with T_vir ∼ 10³−10⁴ K, the streaming velocity delays the cloud collapse until after the halo grows across the atomic-cooling threshold of T_vir ≳ 10⁴ K. The virial temperature for the H−H cases is generally higher than that for the H−H₂ cases, and the mean value of T_vir for each case increases with larger streaming velocity. This trend is more clearly shown in the distributions of J_LW, namely, the mean LW background intensity for the H−H cases is 〈J_LW〉 ≳ 10³, which is ≃10 times higher than that for the H−H₂ cases. The higher value of J_LW is mainly caused by the delay of gas collapse until the halo mass becomes massive enough to be exposed by a larger number of LW source halos. In addition, compressional heating in collapsing clouds is stronger with larger v_bsm, and the minimum LW intensity required to keep isothermal collapse is extended to lower values.

Zoom In Zoom Out Reset image size

In the main progenitors of high-z quasar hosts, massive gas clouds form owing to the significant delay effect of cloud collapse by rapid halo mergers and intense LW irradiation from nearby star-forming galaxies. The mass accretion rate onto the central region of a gravitationally collapsing cloud is approximated as $\dot{M}\simeq {M}_{\mathrm{gas}}/{t}_{\mathrm{ff}}$, where M_gas and t_ff are the gas mass and freefall timescale at the onset of gravitational collapse, respectively. Since the cloud is supported by thermal and kinetic energy of the gas, the accretion rate can be written as ≃${c}_{\mathrm{eff}}^{3}/G$ (Larson 1969; Penston 1969, etc.), which depends only on the gas thermal and kinetic temperature (see below Equation (15)). In the bottom panels of Figure 6, we show the distributions of $\dot{M}\equiv {c}_{\mathrm{eff}}^{3}/G$, for which we adopt the minimum temperature value in the cloud collapse stage at n ≳ 10³ cm⁻³. The accretion rate is broadly distributed over $\dot{M}\simeq 3\times {10}^{-3}-5\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. The vertical line in the bottom panels indicates a reference value of 0.1 M_⊙ yr⁻¹, above which the outer envelope of an accreting protostar is bloated, due to rapid heat injection through mass accretion, and the emission of stellar ionizing photons is strongly suppressed. For v_bsm = 1σ, the majority of the H−H₂ cases yield $\dot{M}\gtrsim 0.1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. With v_bsm > 1σ, all the cases have sufficiently high accretion rates exceeding the reference value (see more discussion in Section 5).

Metal enrichment is considered to be a major obstacle in forming massive BH seeds through star formation because efficient radiative cooling via metal fine-structure lines will induce gas fragmentation and suppress the formation of massive stars. In order to quantify the critical metallicity, we calculate the cooling rate by C _II and O_I, assuming that the number fractions of carbon and oxygen nuclei in the gas phase with respect to hydrogen nuclei are x_C,gas = 0.927 × 10⁻⁴(Z/Z_⊙) and x_O,gas = 3.568 × 10⁻⁴(Z/Z_⊙) (Pollack et al. 1994), and all the carbon and oxygen are in the form of C _II and O_I, respectively. This treatment is justified for warm gas with T ≃ 8000 K (Omukai et al. 2008).

In Figure 7, we present the metal-line cooling rate (dashed) and heating rate associated with mergers and gravitational compression (solid) as a function of the density of gas embedded in the two representative progenitor halos (tree id 2 and 3) with v_bsm = 1σ. In order to examine the cooling effect by metal lines against heating, the H₂ cooling is turned off, and metal-line cooling is calculated but not included in the thermal evolution. The metallicity for each case is set so that the cooling rate is marginally balanced with the heating rate at least once during the collapse phase. Namely, the critical metallicity is estimated as Z_crit ≃ 1.9 × 10⁻³ Z_⊙ (tree id 2) and 2.5 × 10⁻³ Z_⊙ (tree id 3), respectively. These values are higher than the critical metallicity of Z_cirt ∼ 3 × 10⁻⁴ Z_⊙ (in the absence of dust) obtained by Omukai et al. (2008), where the effect of turbulence and merger heating is not included. Although the critical metallicity depends on the relative abundance of metals produced in SN ejecta, we use Z_crit = 2 × 10⁻³ Z_⊙ as a reference value in the following discussion.

Zoom In Zoom Out Reset image size

Throughout this paper, we do not consider the genetic pollution process through mergers of metal-rich minihalos, given that the star-forming efficiency is strongly suppressed by intense LW radiation in the overdense region. However, we note that this treatment is justified only when the "actual" LW intensity is as high as the average value shown in Figure 3. Otherwise, H₂ cooling induces star formation in weak LW radiation pockets. We do not quantify this effect that reduces the number of the main progenitors where gas is kept pristine. As a reference, Lupi et al. (2021) found ∼30% of the atomic-cooling halos in the overdense region to be polluted genetically. Since some of those polluted halos do not belong to the merger history of the final massive quasar host halo, more than 70% of our main progenitor samples should remain pristine (or sufficiently metal-poor). On the other hand, together with the metal enrichment effect, we also exclude the contribution of LW flux from such lower-mass halos, making our treatment conservative.

Next, we discuss the modeling of environmental pollution led by SN-driven bubbles from nearby star-forming halos. One important caveat is that the progenitor halo is assumed to be immediately enriched once the bubble front reaches the halo virial radius. However, the instantaneous enrichment process considered in many previous studies in the literature may not be realistic. In fact, metals in SN ejecta cannot penetrate into the halo center but pollute the halo superficially in the outer region with low densities of ≲10 cm⁻³ (Chen et al. 2017; Chiaki et al. 2018), leaving the gas in the halo interior unpolluted, even for low-mass halos. If more energetic pair-instability SNe occur in nearby source halos, the ejecta with stronger ram pressure deeply penetrate into the target halo and induce metal mixing at the shock front (Chen et al. 2017). To consider this uncertainty, we introduce the metal mixing efficiency f_mix, which is the fraction of metals mixed with the interior gas in the target halo and is treated as a free parameter below.

Another important quantity is the total amount of metals carried into the target halo through multiple SN-driven bubbles. Let us consider a source halo m with a distance of r_s from the target halo with a size of r_vir(M_h). The mass of metals produced by multiple SNe in the source halo is given by ${m}_{\mathrm{met}}={N}_{\mathrm{sn}}{m}_{\mathrm{ej}}$, where ${N}_{\mathrm{sn}}\simeq {m}_{\star }/{m}_{0}$ is the number of SNe and m_ej is the average mass of metals produced by one SN. We here adopt m_ej = 0.746 M_⊙, which corresponds to the metal ejecta mass produced by a 13 M_⊙ stellar progenitor (Chiaki et al. 2018). Assuming that a fraction f_esc,m of the metals is launched isotropically by the SN bubble, the mass of the metals that reach the target halo is given by f_esc,m m_met (r_vir/r_s)²/4. Therefore, due to SN bubbles produced from one source halo, the gas metallicity in the target halo increases by

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}Z\simeq \displaystyle \frac{{m}_{\star }{m}_{\mathrm{ej}}}{{f}_{b}{M}_{{\rm{h}}}{m}_{0}}\cdot \displaystyle \frac{{f}_{\mathrm{esc},{\rm{m}}}{f}_{\mathrm{mix}}}{4}{\left(\displaystyle \frac{{r}_{\mathrm{vir}}}{{r}_{{\rm{s}}}}\right)}^{2}\\ \quad \simeq 9.3\times {10}^{-5}\,{Z}_{\odot }\,{f}_{\mathrm{mix}}\left(\displaystyle \frac{{f}_{\mathrm{esc},{\rm{m}}}}{0.5}\right)\left(\displaystyle \frac{m}{{M}_{{\rm{h}}}}\right){\left(\displaystyle \frac{5{r}_{\mathrm{vir}}}{{r}_{{\rm{s}}}}\right)}^{2},\end{array}\end{eqnarray} \tag{ 21 }$$

where f_esc,m ≃ 0.5 is motivated by 3D high-resolution hydrodynamical simulations of SN-driven galactic outflows (Li et al. 2017).

As discussed in Section 3.1, the LW intensity peaks when the target halo reaches the atomic-cooling threshold because (1) source halos with m_ac,z are the most abundant population in number and (2) two halos with comparable mass are strongly clustered. This circumstance will also maximize the efficiency of environmental enrichment. Assuming M_h = m_ac,z, we estimate the number of source halos with mass of m ≥ m_ac,z located within r_s (≃5r_vir typically) from the target halo for the three representative trees as N_s ≃ 0.4 (tree id 1), 6 (tree id 2), and 86 (tree id 3), respectively. As a result, the gas metallicity in the target halo is calculated as Z = N_sΔZ ≃ 9.3 × 10⁻⁵ Z_⊙ f_mix N_s. Therefore, we obtain the conditions where the environmental enrichment process affects the thermal evolution of gas in the target halo as Z > Z_crit, or equivalently

$$\begin{eqnarray}&&{N}_{{\rm{s}}}\gt 21.5{f}_{\mathrm{mix}}^{-1}.\end{eqnarray} \tag{ 22 }$$

Since f_mix ≤ 1, the gas evolution in the main progenitor surrounded by ≲21 nearby source halos within r_s is unlikely to be affected by metal-line cooling. On the other hand, if the mixing efficiency is as high as f_mix ≳ 0.25, metal enrichment will play an important role in changing the gas evolution in rapidly growing halos (tree id 3), reducing the number fraction of H−H collapse cases (see Figure 5).

Additionally, inhomogeneous density distributions inside the source halos and the nonsteady SFR that forms SNe in the earlier stage change the velocity and shape of expanding bubbles. Those effects could result in either an overestimation or underestimation of the bubble size. To discuss the efficiency of environmental enrichment precisely, we need to further study a variety of situations with different physical parameters, as well as the metal mixing efficiency f_mix. We leave this to future work.

We quantify the critical metallicity and discuss the impact of metal-line cooling on the thermal evolution of gas clouds. However, dynamical evolution of a collapsing cloud with Z ≳ Z_crit that is composed of a warm outer envelope (T ≃ 8000 K) and a cool central core has not been fully understood; especially, long-term behavior of the mass inflow rate onto the central newly formed protostar is still uncertain. Recent cosmological simulations suggest that rapid mass inflows may occur even with metal pollution above the critical metallicity in atomic-cooling halos (Regan et al. 2020a), but widespread star formation limits the final mass of the central star to ≲10⁴ M_⊙ (Regan et al. 2020b). On the other hand, when the metallicity is lower than the critical value, the collapsing gas cloud fragments only at the central region and forms a compact disk, in which a vast majority of the clumps merge with the central protostar via inward migration (Inayoshi & Haiman 2014; Chon & Omukai 2020). As a result, the stellar growth is not quenched by metal pollution. Future work is needed to investigate the star formation in the overdense regions where high-z quasars form to quantify the impact of metal pollution on the gas dynamics.

We apply the obtained mass accretion rate to estimate the final protostellar mass distribution at the end of star formation episodes. Due to the existing angular momentum at large scales, the rapidly accreted pristine gas settles into a disk, which becomes gravitationally unstable and thus results in fragmentation and clump formation (e.g., Oh & Haiman 2002). Most clumps can migrate inward and merge with the central protostar before forming stars (Inayoshi & Haiman 2014), yielding the accretion rate onto the stellar surface through the disk ${\dot{M}}_{\star }=\eta \dot{M}$, where η( < 1) denotes the conversion efficiency from the global accretion rate to that through the accretion disk. Hydrodynamical simulations find that mass accretion through the unstable disk proceeds episodically and the time-averaged value of the efficiency is η ≃ 0.6 (Sakurai et al. 2016; D. Toyouchi et al., in preparation).

When the time-averaged accretion rate is higher than a critical rate, ${\dot{M}}_{\star }\gtrsim {\dot{M}}_{\mathrm{crit}}$, the accreting star continues to expand its envelope with a lower surface temperature of T_eff ≃ 5000 K, from which UV radiation hardly emits. As a result, stellar radiative feedback does not prevent the central star growing via mass accretion (Omukai & Palla 2001; Hosokawa et al. 2012, 2013; Schleicher et al. 2013; Sakurai et al. 2015, 2020; Haemmerlé et al. 2018). Since the value of ${\dot{M}}_{\mathrm{crit}}$ ranges from 0.01 to 0.04 M_⊙ yr⁻¹, depending on the treatment of the stellar evolution calculations and their boundary conditions (Hosokawa et al. 2013; Haemmerlé et al. 2018), we adopt the highest ${\dot{M}}_{\mathrm{crit}}=0.04\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ as a reference value. This choice leads to a lower bound of the stellar/BH mass. With ${\dot{M}}_{\star }\gtrsim {\dot{M}}_{\mathrm{crit}}$, the final stellar mass is determined either by its finite lifetime or by the onset of stellar collapse triggered by the general relativistic (GR) instability (Shibata et al. 2016; see a review by Woods et al. 2019 and references therein). The final mass is also affected by fuel supply through mass accretion onto the star. Woods et al. (2017) have investigated the final mass of stars accreting at a constant rate over ≃0.01–10 M_⊙ yr⁻¹ (radiative feedback is neglected) and found that the final mass linearly increases with the accretion rate below ∼0.03 M_⊙ yr⁻¹ but is saturated around ∼a few × 10⁵ M_⊙ owing to the GR instability. The relation between the critical mass and accretion rate is fitted as

$$\begin{eqnarray}&&{M}_{\star ,\mathrm{GR}}\simeq \left[0.83\mathrm{log}\left(\displaystyle \frac{{\dot{M}}_{\star }}{\,{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)+2.48\right]\times {10}^{5}\,{M}_{\odot }\end{eqnarray} \tag{ 23 }$$

at ${\dot{M}}_{\star }\geqslant 0.1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, which is used for our analysis.

On the other hand, when the stellar accretion rate is lower than the critical rate, ${\dot{M}}_{\star }\lesssim {\dot{M}}_{\mathrm{crit}}$, the star evolves to the main-sequence stage and begins to emit strong ionizing radiation, quenching the stellar growth. Here, we simply consider that ionizing radiation from the star heats the disk surface and thus photoevaporation suppresses the accretion rate (McKee & Tan 2008; Hosokawa et al. 2011; Tanaka et al. 2013). This process becomes important when the ionization front reaches the stellar gravitational influence radius for ionized gas with a temperature of 2 × 10⁴ K, defined by

$$\begin{eqnarray}&&{R}_{\inf ,\star }\equiv \displaystyle \frac{{{GM}}_{\star }}{{c}_{{\rm{s}},\mathrm{ion}}^{2}}\simeq 0.17\,\mathrm{pc}\left(\displaystyle \frac{{M}_{\star }}{{10}^{4}\,{M}_{\odot }}\right),\end{eqnarray} \tag{ 24 }$$

and the ionized gas breaks out through the neutral infalling gas. The photoevaporation rate can be expressed as

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{pe}}\simeq 2.1\times {10}^{-2}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}{\left(\displaystyle \frac{{{\rm{\Phi }}}_{\mathrm{ion}}}{{10}^{52}\,{{\rm{s}}}^{-1}}\right)}^{1/2}{\left(\displaystyle \frac{{R}_{\mathrm{disk}}}{0.1\,\mathrm{pc}}\right)}^{1/2},\end{eqnarray} \tag{ 25 }$$

where Φ_ion is the ionizing photon number flux and R_disk is the size of the accretion disk. The photon flux is approximated as Φ_ion ≃ 1.6 × 10⁵² s⁻¹(M_⋆/10⁴ M_⊙) in the range of 10³ ≲M_⋆/M_⊙ ≲ 10⁵ for main-sequence stars (Johnson et al. 2012). We evaluate the mass outflow rate owing to photoevaporation by setting ${R}_{\mathrm{disk}}\simeq {R}_{\inf ,\star }$ as

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{pe},\min }\simeq 3.5\times {10}^{-2}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\left(\displaystyle \frac{{M}_{\star }}{{10}^{4}\,{M}_{\odot }}\right),\end{eqnarray} \tag{ 26 }$$

which gives a lower bound for the rate because the outflow of ionized gas is mainly driven from larger radii (i.e., a lager surface area). Therefore, equating ${\dot{M}}_{\star }={\dot{M}}_{\mathrm{pe},\min }$, we obtain the feedback-regulated stellar mass as ${M}_{\star ,\mathrm{fb}}\simeq {\dot{M}}_{\star }{t}_{\mathrm{pe}}$ or

$$\begin{eqnarray}&&{\dot{M}}_{\star ,\mathrm{fb}}\simeq 2.9\times {10}^{3}\,{M}_{\odot }\left(\displaystyle \frac{{\dot{M}}_{\star }}{0.01\,{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)\end{eqnarray} \tag{ 27 }$$

at ${\dot{M}}_{\star }\leqslant {\dot{M}}_{\mathrm{crit}}$, where t_pe( ≃ 2.9 × 10⁵ yr) is the characteristic photoevaporation timescale (note that this expression is valid when the stellar lifetime is longer than t_pe). The final mass at the intermediate accretion rate (${\dot{M}}_{\mathrm{crit}}\leqslant {\dot{M}}_{\star }\leqslant 0.1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$) is estimated by performing logarithmic interpolation.

In Figure 8, we show the mass distribution of massive BH seeds formed in the high-z quasar progenitor halos, calculated with the method described above (see also the bottom panels in Figure 6). Note that the number fraction from the different types of gas evolution is stacked at the same mass bin. Without the streaming motion (v_bsm = 0; left panels), the BH mass is widely distributed from 500 M_⊙ (250 M_⊙) to ≳2 × 10⁵ M_⊙ for η = 0.6 (0.3, respectively) with a few peaks corresponding to the virial temperatures of halos when the BHs form by gas collapse. Overall, the cases with high accretion rates ${\dot{M}}_{\mathrm{in}}$ (H−H₂ and H−H cases) are responsible for high-mass BH formation beyond ∼10⁴ M_⊙, while the H₂ case with lower accretion rates yields less massive BHs with <10⁴ M_⊙. The number of BH seeds above 2 × 10⁵ M_⊙ is limited because the GR instability induces direct collapse of accreting supermassive stars. The shape of the mass distribution at 10⁴ ≲ M_•/M_⊙ ≲ 10⁵ depends on the accretion efficiency, namely, the smaller value of η(=0.3) yields a distribution skewed toward lower masses. With nonzero streaming motion (v_bsm = 1σ; right panels), the less massive population with <10⁴ M_⊙ decreases abruptly since nearly all the cases experience the atomic-cooling stage and thus the central stars accrete at high rates without strong radiative feedback. We note that the BH mass distribution for higher streaming velocities is similar to that for v_bsm = 1σ, but their contribution to the total BH mass distribution is less important because regions with v_bsm ≥ 2σ are rarer.

Zoom In Zoom Out Reset image size

As discussed in Section 4, the number fraction of the cases with highest mass accretion rates (H−H cases) would be reduced by the effect of line cooling via atomic carbon and oxygen, which are produced in nearby source halos through SNe and carried into the quasar main progenitor halos with interest. The level of reduction depends on the metal mixing efficiency in the main progenitor, namely, the enrichment effect could be neglected if the mixing efficiency is lower than ∼20%. Nevertheless, the overall shape of the BH mass distribution still holds.

How do those massive seed BHs formed in overdense regions grow to be SMBHs that are observed as high-z quasars at z ≃ 6–7? In previous studies in the literature, the subsequent growth of their BHs via gas accretion and/or mergers and the required conditions have been discussed (e.g., Tanaka & Haiman 2009; Valiante et al. 2016). Recently, large-scale cosmological simulations have been exploring the evolution of SMBHs and the coevolution of their host galaxies, including various feedback processes due to SNe and active galactic nucleus activity with subgrid models. These simulations have generally found that massive seed BHs formed in protogalaxies hardly grow via gas accretion because dense, cold gas is expelled by energetic SN feedback associated with star formation. However, it is worth noting that most simulations in which SN feedback quenches BH growth have focused on "typical" atomic-cooling halos that will grow to ∼10¹⁰–10¹¹ M_⊙ by z ≃ 6 (e.g., Habouzit et al. 2017; Latif et al. 2018; Smith et al. 2018).

On the other hand, as pointed out by Inayoshi et al. (2020), the progenitor halos of high-z quasar hosts with M_h ≃ 10¹² M_⊙ at z ≃ 6 form in rarer regions and have reached M_h ∼ 10⁸ M_⊙ with deeper gravitational potential by the time when star formation takes place (z ∼ 20–35). In such massive halos, a large amount of cold gas is supplied to the nuclear region through filamentary structures of the proto-cosmic web (Di Matteo et al. 2012), and the seed BHs can be fed at high rates significantly exceeding the Eddington limit when the metallicity of inflowing gas is as low as ≲0.01 Z_⊙ (Toyouchi et al. 2021; see also Inayoshi et al. 2016). The critical halo mass required for the onset of rapid mass accretion exceeding the Eddington rate is M_h ≃ 10⁹ M_⊙, almost independent of redshift. Most of the quasar progenitor halos of interest can reach this mass threshold after the birth of seed BHs in ≃20–50 Myr, within which intense starbursts would take place and form protogalaxies. Exploring the nature of BH growth embedded in such a protogalaxy is left for future investigations.

This process is a possible way to form IMBH populations. Observations of IMBHs in the local universe have the potential to constrain high-z BH (seed) formation (see the review by Greene et al. 2020). Furthermore, if those IMBHs form binaries through galaxy mergers and dynamical processes during the cosmic history, the seed-forming channel also provides a significant number of gravitational-wave events (e.g., Chon et al. 2018; Hartwig et al. 2018; Regan et al. 2020b), which will be detectable by space-based gravitational-wave detectors such as the Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017) and Tianqin (Luo et al. 2016), as well as third-generation terrestrial instruments.