This site uses cookies. By continuing to use this site you agree to our use of cookies. To find out more, see our Privacy and Cookies policy. Close this notification
NOTICE: Ukraine: Read IOP Publishing's statement.

The American Astronomical Society (AAS), established in 1899 and based in Washington, DC, is the major organization of professional astronomers in North America. Its membership of about 7,000 individuals also includes physicists, mathematicians, geologists, engineers, and others whose research and educational interests lie within the broad spectrum of subjects comprising contemporary astronomy. The mission of the AAS is to enhance and share humanity's scientific understanding of the universe.

The Institute of Physics (IOP) is a leading scientific society promoting physics and bringing physicists together for the benefit of all. It has a worldwide membership of around 50 000 comprising physicists from all sectors, as well as those with an interest in physics. It works to advance physics research, application and education; and engages with policy makers and the public to develop awareness and understanding of physics. Its publishing company, IOP Publishing, is a world leader in professional scientific communications.

Seeding Supermassive Black Holes with Self-interacting Dark Matter: A Unified Scenario with Baryons

, , and

Published 2021 June 16 © 2021. The American Astronomical Society. All rights reserved.
, , Citation Wei-Xiang Feng et al 2021 ApJL 914 L26

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.



Observations show that supermassive black holes (SMBHs) with a mass of ∼109 M exist when the universe is just 6% of its current age. We propose a scenario where a self-interacting dark matter halo experiences gravothermal instability and its central region collapses into a seed black hole. The presence of baryons in protogalaxies could significantly accelerate the gravothermal evolution of the halo and shorten collapse timescales. The central halo could dissipate its angular momentum remnant via viscosity induced by the self-interactions. The host halo must be on high tails of density fluctuations, implying that high-z SMBHs are expected to be rare in this scenario. We further derive conditions for triggering general relativistic instability of the collapsed region. Our results indicate that self-interacting dark matter can provide a unified explanation for diverse dark matter distributions in galaxies today and the origin of SMBHs at redshifts z ∼ 6–7.

Export citation and abstract BibTeX RIS

1. Introduction

Astrophysical observations of high-redshift quasars indicate that ∼109 M black holes exist when the universe is just 800 Myr old after the Big Bang (z ∼ 7); see Inayoshi et al. (2020) for a review. The origin of these supermassive black holes (SMBHs) is still a mystery. In particular, it is extremely puzzling how they could become so massive in such a short time. A popular idea is that there exist heavy seed black holes in the early universe and they grow massive by accreting baryons (Volonteri 2010; Natarajan 2011). Assuming Eddington accretion, we can relate the black hole mass (MBH) and its seed mass (Mseed) as (Salpeter 1964)

Equation (1)

where Δt is the elapse time and τ = (450/fEdd)[epsilonM /(1 − epsilonM )]Myr is the e-folding time. epsilonM is the radiative efficiency and commonly assumed to be 0.1 (Shakura & Sunyaev 1976), and fEdd is the Eddington ratio following fEdd = Lbol/LEdd, where Lbol is the observed bolometric luminosity and LEdd = 1.3 × 1038(MBH/M)erg s−1 is the Eddington luminosity. epsilonM measures the efficiency of conversion of mass energy to luminous energy by accretion, while fEdd characterizes the efficiency of accretion luminosity.

Consider J1007+2115, the most massive known quasar with MBH ≈ 1.5 × 109 M at z > 7.5 (Yang et al. 2020b). Taking fEdd ≃ 1, we estimate Mseed ∼ 104 M if it forms at z ∼ 30, i.e., Δt = 597 Myr to its observed z = 7.51. Such a seed is too massive to be produced from collapsed Population III stars (Inayoshi et al. 2020), but it could form through the direct collapse of pristine baryonic gas (Bromm & Loeb 2003; Begelman et al. 2006); see also Freese et al. (2010). The latter scenario predicts Mseed ∼ 105–106 M. However, observations show there are high-z SMBHs with fEdd much less than 1 (Matsuoka et al. 2019; Onoue et al. 2019). For example, J1205−0000 is observed at z = 6.7 with MBH = 2.2 × 109 M and fEdd = 0.16 (Onoue et al. 2019). The Eddington accretion then implies it grows from a seed with a mass of 2 × 108 M at z ∼ 30, too heavy to be produced via the direct collapse of gas.

There could be complications in those estimates. For example, the accretion luminosity may not be a constant over time, and SMBHs could experience a rapid accretion phase beyond the Eddington limit (Begelman 1979; Volonteri & Rees 2005; Alexander & Natarajan 2014). In addition, the radiative efficiency depends on black hole spin. For standard thin disks, which drive the hole to maximal spin, epsilonM ∼ 0.42, while for magnetohydrodynamic disks epsilonM ∼ 0.2 (Shapiro 2005) or even lower epsilonM ∼ 0.01 associated with super-Eddington accretion (McKinney et al. 2014; Sadowski et al. 2015). Furthermore, mergers could amplify the black hole mass. It's possible that the existence of SMBHs can be explained after taking into account these effects, but more work is needed to understand how they affect individual ones; see Inayoshi et al. (2020) for further discussion.

In this Letter, we study the scenario of gravothermal collapse of self-interacting dark matter (SIDM; Spergel & Steinhardt 2000; Kaplinghat et al. 2016; Tulin & Yu 2018) in explaining the origin of high-z SMBHs. Dark matter self-interactions can transport heat in the halo over cosmological timescales (Dave et al. 2001; Ahn & Shapiro 2005; Vogelsberger et al. 2012; Rocha et al. 2013). As a gravothermal system, the SIDM halo has negative heat capacity (Lynden-Bell & Wood 1968). The central region could become hot and collapse to a singular state with a finite mass at late stages of the evolution (Balberg & Shapiro 2002; Balberg et al. 2002). Thus SIDM has a natural mechanism in triggering dynamical instability, a necessary condition to form a black hole. Recent studies also show that SIDM is favored for explaining diverse dark matter distributions over a wide range of galactic systems; see Tulin & Yu (2018) for a review. It is intriguing to explore an SIDM scenario that may explain the origin of the high-z SMBHs and observations of galaxies at z ∼ 0.

We adopt a typical baryon mass profile for high-z protogalaxies, and show the collapse time can be shortened by a factor of 100, compared to the SIDM-only case. Even for the self-scattering cross section per unit mass σ/m ∼ 1 cm2 g−1, broadly consistent with the value used to explain galactic observations (Tulin & Yu 2018), the central halo could collapse sufficiently fast to form a seed for z ≳ 7. Depending on the halo mass, this scenario could explain the origin of high-z SMBHs with fEdd ∼ 1 and 0.1. It also has a built-in mechanism to dissipate angular momentum remanent of the central halo, i.e., viscosity induced by the self-interactions. We will further show when the 3D velocity dispersion of SIDM particles in the collapsed central region reaches 0.57c, the general relativistic (GR) instability can be triggered. We demonstrate a unified SIDM scenario that could explain observations of galaxies today and high-z SMBHs. In the appendices, we provide additional information.

2. Gravothermal Evolution

We use a conducting fluid model (Balberg et al. 2002; Koda & Shapiro 2011) to study the gravothermal evolution of an SIDM halo, as it yields high resolution for us to closely trace the collapse process. To capture the influence of baryons, we extend the original model with a baryonic component. The evolution of the halo can be described by the following equations:

Equation (2)

where Mχ , ρχ , νχ , and Lχ are dark matter mass, density, 1D velocity dispersion, and luminosity profiles, respectively, and they are dynamical variables and evolve with time; Mb is the baryon mass profile in the host galaxy; kB is the Boltzmann constant; G is the Newton constant; and Dt denotes the Lagrangian time derivative. Heat conductivity of the dark matter fluid κ can be expressed as $\kappa ={({\kappa }_{\mathrm{lmfp}}^{-1}+{\kappa }_{\mathrm{smfp}}^{-1})}^{-1}$, where ${\kappa }_{\mathrm{lmfp}}\approx 0.27C{\rho }_{\chi }{\nu }_{\chi }^{3}\sigma k/({\mathrm{Gm}}^{2})$ and κsmfp ≈ 2.1νχ k/σ denote conductivity in the long- and short-mean-free-path regimes, respectively, and we set C ≃ 0.75 based on calibrations with N-body simulations (Pollack et al. 2015; Essig et al. 2019). In the short-mean-free-path regime, heat conduction can be characterized by the self-interaction mean free path λ = m/ρχ σ and Kn = λ/H < 1, where $H={({\nu }_{\chi }^{2}/4\pi G{\rho }_{\chi })}^{1/2}$ is the scale height. In the long-mean-free-path regime, it is characterized by H and Kn > 1.

We assume the initial halo follows a Navarro–Frenk–White (NFW) profile (Navarro et al. 1996) with rs and ρs as its scale radius and density, respectively. The boundary conditions are Mχ = 0 at r = 0, Mχ = M200 and Lχ = 0 at r = r200, where M200 and r200 are the virial halo mass and radius, respectively; they are are equivalent to rs and ρs in specifying a halo for a given redshift z. We adopt the baryon mass profile ${M}_{b}{(r)\approx 0.1(4\pi {\rho }_{s}{r}_{s}^{3})(r/{r}_{s})}^{0.6}$, based on cosmological hydrodynamical simulations of protogalaxies at z ∼ 17 (Wise et al. 2008); see Appendix A. As an approximation, we assume the baryon mass profile is static and it does not evolve with time. We recast the fluid equations with dimensionless variables and solve them numerically using the method as in Balberg et al. (2002), Pollack et al. (2015), and Essig et al. (2019). The fiducial quantities relevant for later discussions are ${M}_{0}=4\pi {\rho }_{s}{r}_{s}^{3}$, ${t}_{0}=1/\sqrt{4\pi G{\rho }_{s}}$, and (σ/m)0 = 1/(rs ρs ); hence ${M}_{b}{(r)=0.1{M}_{0}(r/{r}_{s})}^{0.6}$. We then map dimensionless outputs from the simulations to physical ones assuming Planck cosmology, i.e., h = 0.67, Ωm = 0.315, and ΩΛ = 0.685 (Aghanim et al. 2020).

We further elaborate our assumptions made above. The initial halo is optically thin at its characteristic radius if (σ/m)(rs ρs ) < 1 (Pollack et al. 2015). Our numerical study takes (σ/m)(rs ρs ) = 0.2, and hence an NFW initial condition is well justified. As a concrete example, we take the simulated baryon profile from Wise et al. (2008), which is based on collisionless dark matter. In SIDM, the baryon profile could be diffuse because of halo core formation (Vogelsberger et al. 2014; Cruz et al. 2020). However, if baryon infall occurs early before a large core forms, the baryon distribution can be as compact as the one predicted in the collisionless limit (Robertson et al. 2018), or even more dense if SIDM collapse occurs (Sameie et al. 2021), albeit both simulations focus on systems at low redshifts.

Our approximation of a static baryon mass profile could be conservative in estimating the collapse time as the baryons would further contract when the collapse starts; see Sameie et al. (2021). We have also used a Hernquist profile (Hernquist 1990) to model the baryon distribution and obtained similar results if the baryons dominate the central potential as the simulated galaxies in Wise et al. (2008). Given these considerations, our model assumptions are justified. Nevertheless, it would be of interest to simulate high-z protogalaxies in SIDM and test our assumptions in a cosmological environment.

3. Roles of Baryons

Figure 1 shows the gravothermal evolution of the dark matter density versus enclosed mass (solid) in the presence of the baryons (dashed–dotted), where we fix (σ/m)(rs ρs ) = 0.2. The insert panel illustrates the average inner density versus evolution time with (solid) and without (dashed) including the baryon mass. The average inner density $\left\langle {\rho }_{\chi ,\mathrm{in}}\right\rangle $ is calculated within the central region where the enclosed mass equals to that of the seed black hole, as we will explain later. With the baryons, the halo does not form a large density core and it quickly evolves into the collapse phase (Sameie et al. 2018; Yang & Yu 2021). Its density keeps increasing and eventually becomes superexponential in the end. The collapse timescale is tc = 8.41t0, a factor of ∼100 shorter than the one predicted in the SIDM-only case with the same interaction strength. We also performed simulations with a Hernquist baryon profile and found a similar reduction factor in tc if the baryon distribution is as compact as that in Wise et al. (2008).

Figure 1.

Figure 1. Gravothermal evolution of the dark matter density versus enclosed mass in the presence of the baryonic potential (solid), as well as the fixed baryon profile (dashed–dotted). Each dark matter profile is labeled with its corresponding evolution time, and the vertical dotted line indicates the mass of the central halo that would eventually collapse into a seed black hole. The insert panel illustrates the evolution of the averaged dark matter density of the central halo with (solid) and without (dashed) including the baryons.

Standard image High-resolution image

We also see that as the central density increases for t ≳ 8.4t0, the enclosed mass for a central region remains almost a constant Min ≈ 1.8 × 10−3 M0. This is the region where the halo is in the short-mean-free-path regime. A similar phenomenon also occurs without including the baryons (Balberg et al. 2002). For the SIDM-only case we consider, the corresponding Min/M0 value is 4.2 × 10−2, which is larger than the one with the baryons. As the halo evolves further, the density continues increasing and the central halo (Kn ≲ 1) would eventually collapse into a singular state, a seed black hole. We assume the seed mass Mseed = Min, suggested by numerical studies of collapsed massive stars (Saijo et al. 2002).

4. Seeding Supermassive Black Holes

To explain the origin of high-z SMBHs, the initial halo must be sufficiently heavy and collapse fast enough. We first check the scaling relations MinM0M200 and ${t}_{{\rm{c}}}\propto {r}_{s}^{-1}{\rho }_{s}^{-3/2}\propto {M}_{200}^{-1/3}{c}_{200}^{-7/2}{(1+z)}^{-7/2}$ (Essig et al. 2019), where c200 = r200/rs is the halo concentration. Apparently, tc is very sensitive to c200. There is a tight correlation between c200 and M200 for halos at z ≲ 5, but the c200 distribution at higher redshifts is less known. There is a trend that c200 gradually becomes independent of M200 and its median asymptote to c200 ∼ 3 at z ∼ 5–10 (Dutton & Macciò 2014; Zhao et al. 2009). We fix c200 = 3, and leave two parameters M200 and z to vary.

Figure 2 shows benchmarks (red) that could explain the origin of the SMBHs J1205−0000 with the Eddington ratio fEdd = 0.16 (upper panel; Onoue et al. 2019) and J1007+2115 with fEdd = 1.06 (lower panel; Yang et al. 2020b). The black curves indicate their Eddington accretion history reconstructed using Equation (1). For reference, the gray ones denote the accretion history of other high-z SMBHs with fEdd ∼ 0.1 (upper; Matsuoka et al. 2019; Onoue et al. 2019) and those fEdd ∼ 1 (lower; Mortlock et al. 2011; Wu et al. 2015; Bañados et al. 2018; Wang et al. 2018; Yang et al. 2020b). We have checked that all of them (gray) could also be explained within our scenario. The direct collapse of pristine gas could provide massive enough seeds (magenta) for the SMBHs with fEdd ∼ 1, but not those with fEdd ∼ 0.1 (Onoue et al. 2019).

Figure 2.

Figure 2. SIDM benchmarks (red) that could explain the origin of the SMBHs J1205−0000 (labeled as "1"; upper panel) and J1007+2115 ("6"; lower panel) with an observed Eddington growth rate of fEdd = 0.16 (Onoue et al. 2019) and 1.06 (Yang et al. 2020b), respectively. The black curves indicate their Eddington accretion history. For each red arrow, the markers on the higher- and lower-z ends denote initial halo and seed masses, respectively, and the horizontal difference between the two ends indicates the timescale of gravothermal collapse. The blue shaded regions indicate the ratio of the critical density fluctuation to the halo mass variance. The magenta bands denote the mass range of the seed produced via the direct collapse of pristine gas. The gray curves are Eddington growth history of other high-z SMBHs with fEdd ∼ 0.1 (upper) and ∼1 (lower).

Standard image High-resolution image

As the redshift of the initial halo increases, the favored halo mass becomes smaller, because the seed black hole has more time to grow. To explain the origin of the SMBHs with fEdd ∼ 1, the mass is in a range of M200 ∼ 109–1011 M for z ∼ 11–9. For those with fEdd ∼ 0.1, M200 needs to be relatively higher, ∼ 1011–1012 M, as their growth rate is much smaller and a heavier seed is required. We have checked that the overall trend holds for halos with z ≳ 11.

As an example, we take the case with (M200/M, z) = (6.8 × 1011, 8) that seeds J1205−0000, the most challenging SMBH, to demonstrate our derivation. For the halo, ρs ≈ 8.1 × 107 M kpc−3 and rs ≈ 10 kpc. Hence the fiducial parameters are t0 = 15 Myr, M0 = 1.1 × 1012 M, and (σ/m)0 = 5.8 cm2 g−1. The seed mass is Min = 1.8 × 10−3 M0 ≈ 1.9 × 109 M and the collapse time tc = 8.4t0 ≈ 124 Myr, and the self-scattering cross section σ/m ≈ 1.2 cm2 g−1. Since z = 8 is equivalent to t = 642 Myr after the Big Bang, the seed forms at 766 Myr (z = 7). For an SIDM-only halo with the same parameters, we find tc ≈ 103 t0 ≈ 15 Gyr, too long to form a seed.

To speed up the collapse process in the absence of the baryonic influence, one would need to take much larger σ/m and c200 (Pollack et al. 2015), or consider dissipative self-interactions (Choquette et al. 2019; Essig et al. 2019; Huo et al. 2020). For comparison, our scenario predicts σ/m ∼ 1 cm2 g−1 within a minimal elastic SIDM model that has been shown to explain dark matter distributions in the spirals (Ren et al. 2019), Milky Way satellites (Sameie et al. 2020), and dark-matter-deficient galaxies (Yang et al. 2020a). It is important to note dwarf galaxies at present that favor a large density core are those with diffuse baryon distributions (Kaplinghat et al. 2020). Thus their host SIDM halos would still be in the core-expansion phase and a shallow density profile is expected. In addition, in many well-motivated particle physics realizations of SIDM (see Tulin & Yu 2018), the cross section diminishes toward cluster scales. Thus, the stringent bounds on σ/m from galaxy clusters (Kaplinghat et al. 2016; Andrade et al. 2020) can be avoided.

5. Density Fluctuations

For the benchmark cases, the halo mass is in the range of 109–1012 M for z ∼ 11–8. We use the standard Press–Schechter formalism (Press & Schechter 1974) to examine conditions for realizing those halos in the early universe. The halo mass function scales as ${dn}(M,z)/{dM}\propto \exp [-{\delta }_{{\rm{c}}}^{2}(z)/2{\sigma }^{2}(M)]$ (Mo et al. 2010), where δc (z) is the critical density fluctuation at z and σ(M) the mass variance. We shaded the regions with various values of δc (z)/σ(M) in Figure 2 (blue). As expected, the halo mass increases as the density fluctuation increases and more massive halos form at later times.

The halos for seeding the SMBHs with a sub-Eddington (Eddington) accretion rate correspond to δc (z)/σ(M) ∼ 4–6 (3–5). In addition, the baryon concentration of host galaxies needs to be high as well such that the gravothermal collapse could occur fast enough. Thus our scenario predicts that high-z SMBHs should be rare. Indeed, observations show they are extremely rare in the universe. Quasar surveys indicate that the number density of luminous (LAGN ≳ 1046 erg s−1) high-z SMBHs with MBH ∼ 109 M is ≲ 10−7 Mpc−3 (Kulkarni et al. 2019; Inayoshi et al. 2020; Trakhtenbrot 2020). The ratio of their mass to the dynamical (gas+stars) mass is MBH/Mb ∼ 1/100–1/30 (Trakhtenbrot 2020). Taking Mb/M200 ∼ 0.2, we find MBH/M200 ∼ (2–7) × 10−3, broadly consistent with our prediction. We also note that baryon infall can occur for a halo heavier than 5 × 103[(1 + z)/10]1.5 M (Mo et al. 2010), and all of the benchmarks satisfy this condition easily.

6. Angular Momentum

The angular momentum remnant of the inner halo could counter gravity and even slow down the gravothermal collapse. Besides, there is an upper limit on the specific angular momentum of a black hole, JBH/MBH ≤ (G/c)MBH ≈ 1.4 × 10−4(MBH/107 M) kpc km s−1 (Kerr 1963). Consider the benchmark with (M200/M, z) = (6.8 × 1011, 8) again, dark matter particles within the radius rin = 0.063rs ≈ 0.63 kpc of the initial NFW halo would collapse to a seed, where the total enclosed mass is Min. We find Jχ /Min ≈ 8 kpc · km s−1 for the halo within rin, based on a universal fitting formula (Liao et al. 2017). This is a factor of 100 larger than the allowed value for a 109 M seed.

Fortunately, dark matter self-interactions that lead to heat conductivity also induce viscosity, which dissipates the angular momentum remnant. In the long-mean-free-path regime, we find decays as

Equation (3)

where ${J}_{\chi }^{i}$ and ${J}_{\chi }^{f}$ are the initial and final angular momenta of the central halo within rin, respectively, and k = 2(β + 3)/3(β + 5) for a power law of ρχ rβ ; see Appendix C. For an NFW (cored) profile, β = −1 (0) and hence k = 1/3 (2/5). Consider the benchmark, we have rin = 0.063rs , Mχ = Min = 1.8 × 10−3 M0, ρχ (rin) = 14ρs , σ/m = 0.2/(rs ρs ), and ${\nu }_{\chi }({r}_{\mathrm{in}})\,\approx 0.48\sqrt{4\pi G{\rho }_{s}}{r}_{s}$. Taking k = 1/3, we find the timescale for achieving ${J}_{\chi }^{f}\sim {10}^{-2}{J}_{\chi }^{i}$ is Δt ≈ 0.17t0, much shorter than that of gravothermal collapse tc = 8.4t0. We have checked that the other five benchmarks in Figure 2 satisfy the dissipation condition. In SIDM, viscosity and conductivity share the same microscopic nature, and both effects are critical for seeding the SMBHs in our scenario.

7. Relativistic Instability

As the central density increases, the velocity dispersion of the collapsed central region increases as well (Balberg et al. 2002), and it would eventually approach the relativistic limit. To see the fate of the central halo where Kn ≲ 1, we examine conditions for reaching GR instability. Motivated by early studies on globular cluster systems (King 1966; Merafina & Ruffini 1989), we assume that the number density of SIDM particles in the central halo at late stages follow a truncated Maxwell–Boltzmann distribution

Equation (4)

where T, epsilon, and p are temperature, energy, and momentum, respectively, and epsilonc is the cutoff energy, above which the particle escapes to the boundary. Given the distribution in Equation (4), we use the method in Merafina & Ruffini (1989) and solve the Tolman–Oppenheimer–Volkov equation to find density and pressure profiles for the collapsed central region, where we impose the boundary condition kB T = 0.1 mc2. For each configuration, we follow Chandrasekhar's criterion (Chandrasekhar 1964), and calculate the critical adiabatic index γcr and the pressure-averaged adiabatic index 〈γ〉. The sufficient condition for the system to collapse into a black hole is 〈γ〉 < γcr. We will discuss technical details in a companion paper (W.-X. Feng et al. 2021, in preparation).

Figure 3 shows the averaged adiabatic index $\left\langle \gamma \right\rangle $ (red) and the critical index γcr (black) versus 3D velocity dispersion at the center for each configuration denoted as a dot. As the velocity dispersion increases, its averaged index $\left\langle \gamma \right\rangle $ gradually decreases from its nonrelativistic limit for monatomic ideal gas 5/3 toward the ultrarelativistic limit 4/3. In contrast, the critical index γcr increases from the Newtonian limit 4/3 (Shapiro & Teukolsky 1983). This is because as the pressure starts to dominate the energy density toward the GR limit, it destabilizes the system instead. The relativistic instability occurs when the 3D central velocity dispersion approaches 0.57c and $\left\langle \gamma \right\rangle ={\gamma }_{\mathrm{cr}}\approx 1.62$, at which the corresponding fractional binding energy is 0.033 (W.-X. Feng et al. 2021, in preparation).

Figure 3.

Figure 3. The pressure-averaged adiabatic index $\left\langle \gamma \right\rangle $ (red) and the critical index γcr (black) versus the central 3D velocity dispersion for each GR configuration (dot). When $\left\langle \gamma \right\rangle \lt {\gamma }_{\mathrm{cr}}$, the system triggers the GR instability. In the Newtonian limit, $\left\langle \gamma \right\rangle =5/3$ for a monatomic ideal gas, and the instability condition is $\left\langle \gamma \right\rangle \lt 4/3$.

Standard image High-resolution image

8. Conclusions

We have presented a scenario that could explain the origin of high-z SMBHs with Eddington and sub-Eddington accretion rates. The presence of baryons in protogalaxies could deepen the gravitational potential and expedite the gravothermal collapse of an SIDM halo. The favored self-scattering cross section is broadly consistent with the one used to explain diverse dark matter distributions of galaxies. In this scenario, dark matter self-interactions induce viscosity that dissipates the angular momentum remnant of the central halo. The initial halo must be on high tails of density fluctuations, which may explain why high-z SMBHs are extremely rare in observations. We also checked that the GR instability condition can be satisfied. The upcoming and future facilities are expected to search for quasars with a wide range of luminosities (Trakhtenbrot 2020). The observations would provide a more complete picture of populations of high-z SMBHs and further test our scenario.

We thank Stuart L. Shapiro for helpful and friendly discussions, and Masafusa Onoue for clarifying the Eddington growth rate of J2216–0016. W.-X.F. acknowledges the Institute of Physics, Academia Sinica, for the hospitality during the completion of this work. This work was supported by U.S. Department of Energy under grant No. de-sc0008541 (H.-B.Y.), NASA grant 80NSSC20K0566 (H.-B.Y.), and in part by the Kavli Institute for Cosmological Physics at the University of Chicago through an endowment from the Kavli Foundation and its founder Fred Kavli (Y.-M.Z.). This project was made possible through the support of a grant from the John Templeton Foundation (H.-B.Y., ID# 61884). The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation.

Appendix A: The Gas Density Profile

To model the gas distribution of protogalaxies, we adopt simulation results in Wise et al. (2008) (simulation B). Their simulated gas and dark matter distributions are fitted with a single power law of ρb r−2.4 and an NFW profile, respectively. We find the following ansatz works well for the gas:

Equation (A1)

where ρb,s is the scale density of the gas and rs is the scale radius of the simulated halo. The corresponding mass profile is

Equation (A2)

We use simulation data shown in Figure 4 (right, panel b) of Wise et al. (2008) to fix the model parameters, rs = 73 pc, ρs = 2.6 M pc−3, and ρb,s = 0.19 M pc−3; see Figure 4 for comparison. Since 1.67 × ρb,s /ρs ≈ 0.1, we take ${M}_{b}{(r)=0.1{M}_{0}(r/{r}_{s})}^{0.6}$ for the static baryon distribution in our semianalytical simulations, shown as the dashed–dotted line in the left panel of Figure 5. Note that the results from Wise et al. (2008) have high enough resolutions for setting initial conditions in our simulations, where we trace the collapse process with the conducting fluid model.

Figure 4.

Figure 4. Dark matter (red) and gas (black) density profiles after fitting to the simulated ones in Wise et al. (2008); see their Figure 4 (right, panel b).

Standard image High-resolution image
Figure 5.

Figure 5. Left: evolution of dark matter mass profiles (solid) with (σ/m)rs ρs = 0.2, together with the fixed baryon mass profile (dashed–dotted). Each dark matter profile is labeled with its corresponding evolution time. The dashed line indicates the mass of the central halo with Kn < 1. Right: corresponding Kn value versus enclosed mass. The dotted horizontal line indicates Kn = 1, the boundary between short- and long-mean-free-path regimes, where Kn < 1 and > 1, respectively.

Standard image High-resolution image

Appendix B: The Numerical Procedure

The procedure of our semianalytical simulations is largely based on the treatment given in Balberg et al. (2002), Pollack et al. (2015), and Essig et al. (2019). We first translate a relevant physical quantity x to a dimensionless one $\hat{x}$ as $\hat{x}=x/{x}_{0}$, where x0 is its corresponding fiducial value built from the halo parameters ρs and rs , as shown in Table 1.

Table 1. Fiducial Quantities Used in Our Numerical Simulations

${M}_{0}=4\pi {\rho }_{s}{r}_{s}^{3}$ ${(\sigma /m)}_{0}={\left({r}_{s}{\rho }_{s}\right)}^{-1}$
${\nu }_{0}={\left(4\pi G{\rho }_{s}\right)}^{1/2}{r}_{s}$ ${L}_{0}={\left(4\pi \right)}^{5/2}{G}^{3/2}{\rho }_{s}^{5/2}{r}_{s}^{5}$
${t}_{0}={\left(4\pi G{\rho }_{s}\right)}^{-1/2}$  

Download table as:  ASCIITypeset image

The self-gravitating halo is segmented to N = 182 evenly log-spaced concentric shells in radius $\{{\hat{r}}_{1},{\hat{r}}_{2},\cdots ,{\hat{r}}_{N}\}$ with ${\hat{r}}_{1}={10}^{-4}$ and ${\hat{r}}_{N}=100$. The halo is assumed to be in a quasi-hydrostatic equilibrium, and each shell is assumed to be in its local thermal equilibrium. The values of extensive quantities (${\hat{M}}_{i}$, ${\hat{L}}_{i}$) and intensive quantities (${\hat{\rho }}_{i}$, ${\hat{\nu }}_{i}$) are taken as the value at ${\hat{r}}_{i}$ and the average between values at ${\hat{r}}_{i}$ and ${\hat{r}}_{i-1}$, respectively. We fix the baryon mass profile ${\hat{M}}_{b,i}$ as

Equation (B1)

Consequently, we only use one set of Lagrangian zone radius for the halo through the simulations and dynamically update the enclosed baryon mass according to Equation (B1). The workflow is as follows:

  • 1.  
    Compute the initial 1D velocity dispersion profile ${\hat{\nu }}_{\chi ,i}$ based on the input ${\hat{r}}_{i}$, ${\hat{\rho }}_{\chi ,i}$, and ${\hat{\rho }}_{b,i}$ under the hydrostatic equilibrium condition,
    Equation (B2)
  • 2.  
    Compute the luminosity profile ${\hat{L}}_{\chi ,i}$ based on ${\hat{r}}_{i}$, ${\hat{\rho }}_{\chi ,i}$, ${\hat{\nu }}_{\chi ,i}$, and $\hat{\sigma }$ according to Equation (2) of the main text.
  • 3.  
    Allow a small passage of time ${\rm{\Delta }}\hat{t}$ and compute the specific energy change ${\rm{\Delta }}{\hat{u}}_{\chi ,i}$, ${\hat{u}}_{\chi ,i}\equiv 3{\hat{\nu }}_{\chi ,i}^{2}/2$, due to heat conduction,
    Equation (B3)
    where the dark matter density is fixed. We then update ${\hat{u}}_{\chi ,i}$ with ${\hat{u}}_{\chi ,i}+{\rm{\Delta }}{\hat{u}}_{\chi ,i}$. The time step ${\rm{\Delta }}\hat{t}$ is sufficiently small, i.e., $| {\rm{\Delta }}{\hat{u}}_{\chi ,i}/{\hat{u}}_{\chi ,i}| \lt {10}^{-4}$. such that the linear approximations used in step 4 below are valid.
  • 4.  
    Upon updating ${\hat{u}}_{\chi ,i}$, the ith dark matter halo shell is no longer virialized. To return to hydrostatic equilibrium, we perturb ${\hat{r}}_{i}$, ${\hat{\rho }}_{\chi ,i}$, and ${\hat{\nu }}_{\chi ,i}$, while keeping the mass ${\hat{M}}_{\chi ,i}$ and specific entropy ${\hat{s}}_{\chi ,i}=\mathrm{ln}({\hat{\nu }}_{\chi ,i}^{3}/{\hat{\rho }}_{\chi ,i})$ of the shell fixed. We treat mass conservation, specific entropy conservation, kinetic energy conservation, and hydrostatic equilibrium relations, shown in the main text, at the linear order and solve them for all shells simultaneously. For the hydrostatic equilibrium relation, we take the sum of ${\hat{M}}_{\chi ,i}$ and ${\hat{M}}_{b,i}={\hat{M}}_{b}({\hat{r}}_{i})$ to compute the gravitational potential. For numerical accuracy, we iteratively perform the perturbation 10 times until hydrostatic equilibrium is established everywhere.
  • 5.  
    Reestablishing hydrostatic equilibrium gives new values for ${\hat{r}}_{i}$, ${\hat{\rho }}_{\chi ,i}$, and ${\hat{\nu }}_{\chi ,i}$. We return to step 1 and update the luminosity ${\hat{L}}_{i}$.
  • 6.  
    Track the Knudsen number Knλ/H for the innermost shell. The evolution is terminated when Kn drops below 10−4.

The above procedure is coded in C++ with the eigen 3.2 library for linear algebra (Guennebaud et al. 2010).

In Figure 5, we show evolution of dark matter mass profile (left panel) and the corresponding Kn value versus enclosed mass (right panel). These results are complementary to those presented in Figure 1 of the main text.

Appendix C: Angular Momentum Dissipation

Dark matter self-interactions provide an important avenue to transport angular momentum. To estimate this effect, we keep track of the collapsing central halo in a Lagrangian zone manner, i.e., the number of particles in each mass sphere is conserved, while allowing its corresponding radius to change over the evolution. This is consistent with the method we use for solving the conduction fluid equations. We further assume the mass distribution is spherically symmetric, which is particularly motivated in SIDM (Dave et al. 2001; Peter et al. 2013), and write the moment of inertia as ${I}_{\chi }={{kM}}_{\mathrm{in}}{r}_{\mathrm{in}}^{2}$, where k = 2(β + 3)/3(β + 5) for a power law of ρχ rβ . For an NFW profile, β = − 1 and k = 1/3. For a density core, we have β = 0 and k = 2/5. The angular momentum is given by

Equation (C1)

where ω is the rotational frequency of the inner region and rin is its boundary that changes with time. This leads to

Equation (C2)

The bulk velocity is ${v}_{\phi }={r}_{\mathrm{in}}\omega \sin \theta $ along the ϕ-direction of the rotational axis. θ is the polar angle. The bulk velocity increases (decreases) dvϕ /dt > 0 (<0) as the boundary of the inner region shrinks (expands) drin/dt < 0 (>0). An increase in vϕ will drag the ambient regions just outside the boundary and exert a shear pressure on the boundary bulk surface,

Equation (C3)

where m is the dark matter particle mass, ${A}_{\mathrm{in}}=4\pi {r}_{\mathrm{in}}^{2}$ is the surface area of the inner region, ${\eta }_{{r}_{\mathrm{in}}}\equiv \eta ({r}_{\mathrm{in}})$ is the viscosity of the SIDM fluid, and ${N}_{{r}_{\mathrm{in}}}$ is number of particles, on the bulk surface. The subscript + / − indicates if the quantity increases/decreases. Given the Lagrangian zone setup, ${N}_{{r}_{\mathrm{in}}}$ is a constant through the evolution, and from Equation (C3) we can show

Equation (C4)

The bulk momentum can be transported out through the shear pressure due to viscosity. Combining Equations (C2) and (C4), we obtain the rate of momentum transport to the surroundings:

As the total angular momentum is conserved, the loss of the angular momentum of the inner region (shrinking case) is given by

Equation (C5)

We have

Equation (C6)

where ${J}_{\chi ,\mathrm{in}}^{i}$ is the initial angular momentum of the inner region. As for conductivity (Balberg et al. 2002), the viscosity for both long-mean-free-path and short-mean-free-path regimes can be combined into a single expression,

Equation (C7)

where we have used the gravitational scale height $H=\sqrt{{\nu }^{2}/4\pi G\rho }$, the mean free path λ = 1/n σ, and the relaxation time tr = 1/(α n ν σ) with number density n, cross section σ, and α = (16/π)1/2 ≈ 2.26 for hard spheres.

We evaluate η at the boundary rin and take $\bar{v}\simeq \sqrt{3}\nu $, and obtain

Equation (C8)

in the long-mean-free-path limit, and

Equation (C9)

in the short-mean-free-path limit, where we have used the density ${\rho }_{{r}_{\mathrm{in}}}=\rho ({r}_{\mathrm{in}})={mn}({r}_{\mathrm{in}})$ and ${\nu }_{{r}_{\mathrm{in}}}=\nu ({r}_{\mathrm{in}})$ the 1D velocity dispersion at the boundary.

For the benchmark discussed in Section 6,we estimate the characteristic timescale to dissipate the angular momentum remnant Δt ∼ 0.2t0 with fixed Min = 1.8 × 10−3 M0 and its corresponding rin = 0.063rs , i.e.,the radius at which the enclosed mass of the initial NFW halo is Min. We do not take into account the change in the radius of the mass sphere and the velocity dispersion, as they are negligible for Δt ∼ 0.2t0; see the left panel of Figure 5. The estimation is based on Equation (C8), and this is well justified as the collapsed region is in the long-mean-free-path regime at rrin for Δt ∼ 0.2t0, as shown in the right panel of Figure 5. For a sphere with smaller Min more toward the center, the timescale is even shorter, as we can see from the scaling relation ${\rm{\Delta }}t\sim {{GM}}_{\chi }/({\rho }_{\chi }{\nu }_{\chi }^{3}r\sigma /m)$. For an NFW halo, we have ρχ r−1, Mχ r2, νχ r0.2, and hence Δtr1.4. For a cored halo, ρχ r0, Mχ r3, νχ r0, and Δtr2. For both cases, Δt decreases as r reduces. Thus it is reasonable to collectively treat particles in the mass sphere Min = 1.8 × 10−3 M0, which would collapse to a seed, and estimate the overall timescale for dissipating angular momentum.

Appendix D: The Adiabatic Index

We consider a perfect fluid with energy density ρ(r)c2 and pressure p(r) in a Schwarzschild metric (Chandrasekhar 1964)

where ${e}^{2{\rm{\Phi }}}=\exp [2{\int }_{r}^{\infty }({dp}/{dr}^{\prime} )/(p+\rho {c}^{2}){dr}^{\prime} ]$ and ${e}^{2{\rm{\Lambda }}}\,={[1-2{GM}(r)/{{rc}}^{2}]}^{-1}$. The critical adiabatic index is

Equation (D1)

and the pressure-averaged adiabatic index is 〈γ〉 ≡ ∫e3Φ+Λ γ(r)pr2 dr/(∫e3Φ+Λ pr2 dr). These expressions are fully relativistic, and we will provide their derivations in a companion paper (W.-X. Feng et al. 2021, in preparation).

Appendix E: The Sample of High-z SMBHs

In Table 2, we list high-z SMBHs shown in Figure 2 of the main text, in the order of their labeling number in the figure. The Eddington ratio is calculated as fEdd = Lbol/LEdd, where Lbol is the observed bolometric luminosity and LEdd = 1.3 × 1038 (MBH/M) erg s−1 is the Eddington luminosity.

Table 2. The Sample of High-z SMBHs Shown in Figure 2 of the Main Text

LabelName MBH (109M) z fEdd Ref.
1J1205−0000 ${2.2}_{-0.6}^{+0.2}$ ${6.699}_{-0.001}^{+0.007}$ ${0.16}_{-0.02}^{+0.04}$ Onoue et al. (2019)
2J1243+0100 ${0.33}_{-0.2}^{+0.2}$ ${7.07}_{-0.01}^{+0.01}$ ${0.34}_{-0.02}^{+0.02}$ Matsuoka et al. (2019)
3J2239+0207 ${1.1}_{-0.2}^{+0.3}$ ${6.245}_{-0.007}^{+0.008}$ ${0.17}_{-0.05}^{+0.04}$ Onoue et al. (2019)
4J2216−0016 ${0.7}_{-0.23}^{+0.14}$ ${6.109}_{-0.008}^{+0.007}$ ${0.15}_{-0.03}^{+0.05}$ Onoue et al. (2019)
5J1208−0200 ${0.71}_{-0.52}^{+0.24}$ ${6.144}_{-0.010}^{+0.008}$ ${0.24}_{-0.08}^{+0.18}$ Onoue et al. (2019)
6J1007+2115 ${1.5}_{-0.2}^{+0.2}$ ${7.5149}_{-0.0004}^{+0.0004}$ ${1.06}_{-0.2}^{+0.2}$ Yang et al. (2020b)
7J1342+0928 ${0.78}_{-0.19}^{+0.33}$ ${7.5413}_{-0.0007}^{+0.0007}$ ${1.5}_{-0.4}^{+0.5}$ Bañados et al. (2018)
8J1120+0641 ${2.0}_{-0.7}^{+1.5}$ ${7.085}_{-0.003}^{+0.003}$ ${1.2}_{-0.5}^{+0.6}$ Mortlock et al. (2011)
9J0038−1527 ${1.33}_{-0.25}^{+0.25}$ ${7.021}_{-0.005}^{+0.005}$ ${1.25}_{-0.19}^{+0.19}$ Wang et al. (2018)
10J0100+2802 ${12.4}_{-1.9}^{+1.9}$ ${6.30}_{-0.01}^{+0.01}$ ${0.99}_{-0.22}^{+0.22}$ Wu et al. (2015)

Download table as:  ASCIITypeset image

Please wait… references are loading.