# Rapid Growth of Seed Black Holes during Early Bulge Formation

Kohei Inayoshi  Kavli Institute for Astronomy and Astrophysics, Peking
University, Beijing 100871, China  Riouhei Nakatani  RIKEN Cluster for
Pioneering Research, 2-1 Hirosawa, Wako-shi, Saitama 351-0198, Japan  Daisuke
Toyouchi  Research Center for the Early Universe (RESCEU), The University of
Tokyo Hongo, 7-3-1, Bunkyo-ku Tokyo, 113-0033, Japan  Takashi Hosokawa
Department of Physics, Kyoto University, Kyoto 606-8502, Japan  Rolf Kuiper
Zentrum für Astronomie der Universität Heidelberg, Institut für Theoretische
Astrophysik, Albert-Ueberle-Straße 2, 69120 Heidelberg, Germany  Masafusa
Onoue  Kavli Institute for Astronomy and Astrophysics, Peking University,
Beijing 100871, China  Kavli Institute for the Physics and Mathematics of the
Universe (Kavli IPMU, WPI), The University of Tokyo, Chiba 277-8583, Japan
Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

###### Abstract

We study the early growth of massive seed black holes (BHs) via accretion in
protogalactic nuclei where the stellar bulge component is assembled,
performing axisymmetric two-dimensional radiation hydrodynamical simulations.
We find that when a seed BH with $M_{\bullet}\sim 10^{5}~{}M_{\odot}$ is
embedded in dense metal-poor gas ($Z=0.01~{}Z_{\odot}$) with a density of
$\gtrsim 100~{}{\rm cm}^{-3}$ and bulge stars with a total mass of
$M_{\star}\gtrsim 100~{}M_{\bullet}$, a massive gaseous disk feeds the BH
efficiently at rates of $\gtrsim 0.3-1~{}M_{\odot}~{}{\rm yr}^{-1}$ and the BH
mass increases nearly tenfold within $\sim 2$ Myr. This rapid accretion phase
lasts until a good fraction of the gas bounded within the bulge accretes onto
the BH, although the feeding rate is regulated owing to strong outflows driven
by ionizing radiation emitted from the accreting BH. The transient growing
mode can be triggered for seed BHs formed in massive dark-matter halos with
masses of $\gtrsim 10^{9}~{}M_{\odot}$ at $z\sim 15-20$ (the virial
temperature is $T_{\rm vir}\simeq 10^{5}~{}{\rm K}$). The host halos are
heavier and rarer than those of typical first galaxies, but are more likely to
end up in quasar hosts by $z\simeq 6$. This mechanism naturally yields a mass
ratio of $M_{\bullet}/M_{\star}>0.01$ higher than the value seen in the local
universe and the existence of such overmassive BHs provides us a unique
opportunity of detecting highly accreting seed BHs at $z\sim 15$ with AB
magnitude of $m_{\rm AB}\sim 26-29$ mag at $2~{}\mu{\rm m}$ (rest-frame 10 eV)
by the upcoming observations by the James Webb Space Telescope and Nancy Grace
Roman Space Telescope.

Supermassive black holes (1663); Quasars (1319); High-redshift galaxies (734)

##  1 Introduction

Supermassive black holes (SMBHs) are one of the most essential objects that
consist of galaxies and their coevolution with host galaxies are suggested by
the empirical relations between BH masses and galaxy properties seen in the
local universe (e.g., Magorrian et al., 1998; Ferrarese & Merritt, 2000;
Kormendy & Ho, 2013). Nevertheless, the origin of those SMBHs remains one of
the most intriguing and longest-standing unsolved puzzles in astrophysics. In
the past decades, the discoveries of very luminous quasars in the high-
redshift universe ($z>6$) have revealed that SMBHs with masses greater than
$\sim 10^{9}~{}M_{\odot}$ formed within the first billion years after the Big
Bang (e.g., Fan, 2006; Willott et al., 2010; Mortlock et al., 2011; Wu et al.,
2015; Bañados et al., 2018; Wang et al., 2021; Yang et al., 2021). Beyond the
rarest and most massive SMBHs that represent the tip of the iceberg of the
high-redshift BH population, a recent wide-field survey with Subaru Hyper
Suprime-Cam (HSC) brings the total number of $z>6$ quasars to nearly 200 and
enables us to construct their luminosity function down to the faint-end regime
(Matsuoka et al., 2016, 2018). The BH mass and host galaxy’s properties of
those low-luminosity quasars at $z>6$ have also been extensively investigated
in order to construct the BH mass distribution at $z>6$ and understand the
early coevolution between BHs and galaxies (e.g., Onoue et al., 2019; Izumi et
al., 2019).

The existence of those high-redshift quasars requires their quick assembly
mechanisms (Volonteri, 2012; Haiman, 2013; Inayoshi et al., 2020); for
instance rapid gas collapse into the nuclei of early protogalaxies (Volonteri
& Rees, 2005; Inayoshi et al., 2016; Toyouchi et al., 2021) and the formation
of massive heavy seed BHs through primordial star formation (Omukai, 2001;
Bromm & Loeb, 2003; Lodato & Natarajan, 2006; Shang et al., 2010; Latif et
al., 2013; Regan et al., 2014; Inayoshi et al., 2014; Becerra et al., 2015;
Chon et al., 2018; Wise et al., 2019; Regan et al., 2020; Sassano et al.,
2021) and runaway stellar mergers in dense regions (Devecchi & Volonteri,
2009; Sakurai et al., 2017; Tagawa et al., 2020; Chon & Omukai, 2020). In most
previous studies, those efficient BH assembly processes are considered to take
place in atomically-cooling dark-matter (DM) halos with virial temperature of
$T_{\rm vir}\gtrsim 10^{4}~{}{\rm K}$, (corresponding to DM halo masses of
$M_{\rm h}\gtrsim 10^{7}~{}M_{\odot}$ at $z\sim 10-15$), where typical first-
galaxies would form (Bromm & Yoshida, 2011). However, a majority of
cosmological simulations have suggested that the subsequent growth of seed BHs
formed in typical first galaxies with shallow gravitational potential would be
strongly quenched by various feedback processes including supernova (SN)
explosions of forming massive stars and radiative/mechanical outputs
associated with BH accretion (e.g., Dubois et al., 2013; Prieto & Escala,
2016; Anglés-Alcázar et al., 2017; Habouzit et al., 2017; Latif et al., 2018).

Recently, semi-analytical studies by Lupi et al. (2021) and Li et al. (2021b)
have proposed that BH seeding processes preferentially occur in highly-biased,
overdense regions of the universe where quasars form by $z\simeq 6$. The
quasar progenitor halos are substantially heavier than those hosting first
galaxies and their average mass exceeds $M_{\rm h}\sim 10^{9}~{}M_{\odot}$ by
$z\sim 20$ via quick assembly after the birth of BH seeds with
$M_{\bullet}\sim 10^{4-5}~{}M_{\odot}$. Large-scale cosmological simulations
demonstrate that cold gas streams feed the centers of such massive halos and
sustain high mass accretion rates until the mass of the galaxy reaches $\sim
10^{12}~{}M_{\odot}$ (Li et al. 2007; Sijacki et al. 2009; Di Matteo et al.
2012; Zhu et al. 2020; see also Dekel & Birnboim 2006). However, it is worth
emphasizing that due to numerical limitations, cosmological simulations with a
spatial resolution of $\sim O({\rm kpc})$ treat feedback effects with subgrid
models. Therefore, it is essential to directly resolve physical processes in
the nuclear region and investigate the growth phases of seed BHs hosted in
massive halos that will end up as high-redshift quasars.

The nature of mass accretion onto BHs in galactic nuclei has been explored by
analytical and numerical work. In the recent decade, radiation hydrodynamical
(RHD) simulations have shown that rapid accreting flows can feed the BH at
super-Eddington rates as long as a sufficient amount of gas already exists or
is efficiently supplied to the vicinity of the BH horizon scale of $r_{\rm
nuc}\simeq 10^{2-3}~{}r_{\rm Sch}$ (Ohsuga et al., 2005; Jiang et al., 2014;
Sa̧dowski et al., 2015; McKinney et al., 2015), where $r_{\rm Sch}\equiv
2GM_{\bullet}/c^{2}$ is the Schwarzschild radius of the BH, where $G$ is the
gravitational constant and $c$ is the speed of light. Note that the nuclear
size for a seed BH is as small as $r_{\rm nuc}\simeq 10^{-5}~{}{\rm
pc}~{}(M_{\bullet}/10^{5}~{}M_{\odot})$. However, ionizing radiation emitted
from the accreting BH heat the gas inflowing from the BH gravitational
influence radius of $\sim 1-10~{}{\rm pc}$, within which the BH gravitational
energy dominates the gas thermal energy. Thus, the mass accretion rate is
generally self-regulated below the Eddington value (Ciotti & Ostriker, 2001;
Milosavljević et al., 2009; Johnson et al., 2011; Park & Ricotti, 2011, 2012;
Jeon et al., 2012; Smith et al., 2017).

Accordingly, RHD simulations that cover the BH influence radii and resolve
sub-parsec scales showed that when the BH is embedded in a dense gas cloud so
that the gas inflowing rate from the BH influence radius substantially exceeds
the Eddington value, the global quasi-steady structure of rapid inflows is
maintained (Inayoshi, Haiman, & Ostriker, 2016). In this case, since the
emergent radiation flux is reduced by photon trapping and dust absorption in
the flow, the gas inflow is not prevented by radiative feedback but leads to
collapse of the ionized region surrounding the BH (see also Sakurai et al.,
2016; Yajima et al., 2017; Park et al., 2020). Moreover, anisotropic radiation
produced from the accreting BH toward the polar regions moderates the negative
feedback effect (Sugimura et al., 2017; Takeo et al., 2018) and mechanical
feedback associated with strong outflows completely evacuates the polar
regions but does not affect the gas dynamics of disk accretion (Takeo et al.,
2020). Recent three-dimensional RHD simulations have shown that massive gas
supply can be sustained when the dusty disk becomes sufficiently optically
thick to ionizing radiation (Toyouchi et al., 2021). Applying the conditions
to seed BHs in high-redshift protogalaxies, they found that BHs formed in
massive DM halos with $M_{\rm h}\gtrsim 10^{9}~{}M_{\odot}$ in the early
universe can experience a rapidly growing phase. We note that the halo
conditions nicely agree with those suggested by previous studies in the
context of seed formation (Lupi et al., 2021; Li et al., 2021b) and subsequent
BH growth (e.g., Di Matteo et al., 2012).

However, those RHD simulation studies that focus on the intermediate physical
scales within the gravitational influence radius of a seed BH ($\lesssim
1-10~{}{\rm pc}$) have not taken into account the existence of its host
galaxy, specifically stellar bulge components. A previous work by Park et al.
(2016) addressed the impact of bulge gravity on BH accretion, using
spherically symmetric one-dimensional RHD simulations that assume metal-free
gas with a uniform initial distribution. They found that when the bulge mass
is greater than a critical mass of $\sim 10^{6}~{}M_{\odot}$, the bulge
gravity leads to gas accumulation and increases the BH fueling rate.

In this paper, we investigate the dynamics of accreting and outflowing gas in
the intermediate missing region surrounding a massive seed BH in
protogalaxies, performing a series of axisymmetric two-dimensional (2D) RHD
simulations with a sufficiently large computation domain covering $\sim
0.1-100~{}{\rm pc}$ that enables us to capture the multi-scale physics
properly. We consider a massive metal-polluted cloud with metallicity of
$Z=0.01~{}Z_{\odot}$ concentrated inside a massive DM halo, adopting initial
conditions of gas motivated by cosmological simulations. We consider bulge
stars formed in the protogalaxies that cause two relevant effects on BH
feeding. First, since the bulge mass is considered to be heavier than the BH
mass (for instance, $M_{\star}/M_{\bullet}\simeq 10^{3}$ in the local
universe), a larger amount of gas can be accumulated from larger scales of
$\sim O({\rm kpc})$ due to additional stellar gravity and can efficiently be
delivered into the BH gravitational sphere of influence. This would lead to
rapid growth of BHs, as seen in the previous study. On the other hand, since
the bulge consists of a number of stars (presumably young stars in
protogalaxies), photoionization and heating caused by stellar radiation affect
the thermal properties of the surrounding gas and would prevent the gas from
feeding the BH. In our study, we quantify the critical conditions required for
rapid growth of BHs induced by the galactic bulges.

One of the most interesting questions is how and when the BH-galaxy
correlations have been established. To answer this, we apply our RHD
simulation results to seed formation in the cosmological framework of galaxy
formation and translate the conditions required for rapid BH accretion in
terms of the properties of the host protogalaxies and halos. We find that this
BH-growth mechanism naturally yields a mass ratio of
$M_{\bullet}/M_{\star}>0.01$ at $z\gtrsim 10$, which is significantly higher
than the canonical value seen in the local universe (Kormendy & Ho, 2013) but
is suggested from the observations of the brightest quasars at $z>6$ (Wang et
al., 2010, 2013; Pensabene et al., 2020). Furthermore, the highly accreting
seed BHs with $M_{\bullet}\simeq 10^{5-6}~{}M_{\odot}$ produce intense
radiation with luminosities of $\sim 10^{44-45}~{}{\rm erg~{}s}^{-1}$.
Therefore, the existence of such overmassive BHs provides us a unique
opportunity of detecting seed BHs in the very early universe with upcoming
deep near-infrared observations by e.g., James Webb Space Telescope (JWST) and
Nancy Grace Roman Space Telescope (RST). The result would shed the light on
the unresolved coevolutionary process between BHs and galaxies in the early
universe.

The rest of the paper is organized as follows. In Section 2, we first describe
the numerical method and settings of our 2D RHD simulations. In Section 3, we
present the simulation results and their dependence on physical parameters
that characterize the bulge and halo properties. Estimating the radiative
luminosity from the highly accreting BHs, we quantify the detectability of the
signature in upcoming observations by JWST and RST. In Section 4, based on the
simulation results, we provide a theoretical explanation for the conditions
required to promote rapid accretion onto BH seeds in protogalaxies. In Section
5, we discuss the conditions for rapid BH accretion in the framework of the
hierarchical structure formation, the $M_{\bullet}/M_{\star}$ ratio expected
in the early bulge formation stage, and implications for the measurements of
quasar lifetime. Finally, we summarize our findings in Section 6.

##  2 Method

###  2.1 Basic equations

We solve the axisymmetric two-dimensional hydrodynamical equations using the
open source code PLUTO (Mignone et al., 2007), which has been modified to
study massive star formation and the evolution of protoplanetary disks (e.g.,
Kuiper et al., 2010, 2011; Nakatani et al., 2018a, b). The detailed
description of the code and its applications to many other astrophysical
systems are shown in Kuiper et al. (2020). In particular, we update the
specific version of the code to investigate the physics of BH feeding and
feedback in the high-redshift protogalactic nuclei.

The basic equations we solve are following: the equation of continuity,

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mbox{\boldmath$v$})=0,$ (1)

and the equation of motion,

$\displaystyle\frac{\partial(\rho v_{r})}{\partial t}+\nabla\cdot(\rho
v_{r}\mbox{\boldmath$v$})=$ $\displaystyle-\frac{\partial p}{\partial
r}+\rho\frac{v_{\theta}^{2}+v_{\phi}^{2}}{r}-\rho\frac{\partial\Phi}{\partial
r}$ $\displaystyle+(\nabla\cdot\mbox{\boldmath$\sigma$})_{r}+f_{{\rm rad},r},$
(2) $\displaystyle\frac{\partial(\rho v_{\theta})}{\partial
t}+\nabla\cdot(\rho v_{\theta}\mbox{\boldmath$v$})=$
$\displaystyle-\frac{1}{r}\frac{\partial
p}{\partial\theta}-\rho\frac{v_{\theta}v_{\phi}}{r}+\rho\frac{v_{\phi}^{2}}{r}\cot\theta$
$\displaystyle+(\nabla\cdot\mbox{\boldmath$\sigma$})_{\theta}+f_{\rm
rad,\theta},$ (3) $\frac{\partial(\rho v_{\phi})}{\partial
t}+\nabla^{l}\cdot(\rho
v_{\phi}\mbox{\boldmath$v$})=(\nabla\cdot\mbox{\boldmath$\sigma$})_{\phi}$ (4)

where $\rho$ is the mass density, $v$ is the velocity, $p$ is the gas
pressure, $\Phi$ is the gravitational potential, $\sigma$ is the stress tensor
due to viscosity, and $f_{{\rm rad},r}$ and $f_{{\rm rad},\theta}$ are the
radial and polar components of the radiation force.

We solve the energy equation of

$\frac{\partial E}{\partial t}+\nabla\cdot(H\mbox{\boldmath$v$})=-\rho
v_{r}\frac{\partial\Phi}{\partial
r}+(\mbox{\boldmath$\sigma$}\cdot\nabla)\mbox{\boldmath$v$}+\Gamma-\Lambda,$
(5)

where $E=(1/2)\rho\mbox{\boldmath$v$}^{2}+\rho e$ and $H=E+p$ are the total
energy density and enthalpy per unit volume, and $e$ is the internal energy
per unit mass. The equation of state for ideal gas is assumed as

$p=(\gamma-1)\rho e,$ (6)

where the adiabatic exponent $\gamma$ depends on the chemical abundances and
gas temperature (e.g., Omukai & Nishi, 1998). The first and second terms on
the right-hand-side present work done by the gravitational force and viscous
heating. The last two terms are volumetric radiative cooling and heating (in
units of ${\rm erg~{}s^{-1}~{}cm^{-3}}$).

The gravitational potential is set to

$\Phi=-\frac{GM_{\bullet}}{r}+\Phi_{\star}(r)+\Phi_{\rm DM}(r),$ (7)

where $M_{\bullet}$ is the BH mass, $\Phi_{\star}$ and $\Phi_{\rm DM}$ are the
gravitational potential owing to bulge stars and a DM halo, respectively. We
here do not consider the self-gravity of gas. For the bulge stellar potential,
we assume that the stellar-mass density profile follows a Hernquist profile
and thus the potential is given by

$\Phi_{\star}(r)=-\frac{GM_{\star}}{r+r_{\rm c}},$ (8)

(Hernquist, 1990), where $M_{\star}$ is the total bulge mass and $r_{\rm c}$
is the core radius. Note that the bulge mass increases with time at a given
star formation rate (see §2.5). With the stellar mass distribution, the
velocity dispersion of stars at the half-mass radius of $R_{\rm
e}=(1+\sqrt{2})r_{\rm c}$ is analytically calculated as $\sigma_{\rm e}\simeq
0.23\sqrt{GM_{\star}/r_{\rm c}}$ (see Appendix A)111Park et al. (2016)
calculated the bulge size so that the average stellar mass density within
$r_{\rm c}$ is comparable to that of the Milky Way. Namely, they adopt $r_{\rm
c}=50~{}{\rm pc}(M_{\star}/10^{7}~{}M_{\odot})^{1/3}$, which yields
$M_{\star}\simeq 2.5\times 10^{11}~{}M_{\odot}(\sigma_{\rm e}/200~{}{\rm
km~{}s}^{-1})^{3}$. The bulge mass in this relation is significantly
overmassive compared to that observed in the local universe. Assuming that the
bulge stellar component in high-redshift protogalaxies obeys the
$M_{\star}$-$\sigma_{\rm e}$ relation (Kormendy & Ho, 2013), the core radius
is set to

$r_{\rm c}\simeq 6.5~{}{\rm
pc}\left(\frac{M_{\star}}{10^{7}~{}M_{\odot}}\right)^{0.47}.$ (9)

Note that the observed slope of the $R_{\rm e}$-$M_{\star}$ relation is known
to depend on the stellar mass, galaxy type, and redshift (e.g., Shen et al.,
2003; Mowla et al., 2019; Li et al., 2021a; Kawinwanichakij et al., 2021). For
the DM potential, we consider a Navarro, Frenk & White (NFW) potential
(Navarro et al., 1997; Mo et al., 1998), approximately written as

$\Phi_{\rm DM}(r)=-\frac{2k_{\rm B}T_{\rm vir}}{\mu m_{\rm
p}}~{}\frac{\ln(1+r/r_{\rm s})}{r/r_{\rm s}}~{}f(c_{\rm N}),$ (10)

where $k_{\rm B}$ is the Boltzmann constant, $\mu$ is the mean molecular
weight of gas, $T_{\rm vir}$ and $r_{\rm vir}$ are the virial temperature and
virial radius of the halo, $r_{\rm s}(\equiv r_{\rm vir}/c_{\rm N})$ is the
characteristic radius of the NFW density profile, $c_{\rm N}$ is the
concentration factor, and $f(c_{\rm N})=c_{\rm N}/[\ln(1+c_{\rm N})-c_{\rm
N}/(1+c_{\rm N})]$. We model the dependence of the mean concentration on halo
mass $M_{\rm h}$ and redshift $z$: $c_{\rm N}\simeq 1.56~{}(M_{\rm
h}/10^{9}~{}M_{\odot})^{-0.13}[(1+z)/21]^{-1}$ (Bullock et al., 2001).

In our axisymmetric simulations without magneto-hydrodynamical (MHD) effects,
angular momentum transport in the accreting flow is given by imposing explicit
viscosity. The viscous stress tensor is given by

$\sigma_{ij}=\rho\tilde{\nu}\left[\left(\frac{\partial v_{j}}{\partial
x_{i}}+\frac{\partial v_{i}}{\partial
x_{j}}\right)-\frac{2}{3}(\nabla\cdot\mbox{\boldmath$v$})\delta_{ij}\right],$
(11)

where $\tilde{\nu}$ is the shear viscosity and the bulk viscosity is
neglected. To mimic angular momentum transport associated with MHD turbulence
driven by the magneto-rotational instability (MRI) in a sufficiently ionized
disk (e.g., Balbus & Hawley, 1998; Stone & Pringle, 2001; McKinney & Gammie,
2004; Bai, 2011; Narayan et al., 2012), we assume the azimuthal components of
the shear tensor are non-zero and, in spherical polar coordinates, are given
by

$\sigma_{r\phi}=\rho{\tilde{\nu}}\frac{\partial}{\partial
r}\left(\frac{v_{\phi}}{r}\right),$ (12)
$\sigma_{\theta\phi}=\rho{\tilde{\nu}}\frac{\sin\theta}{r}\frac{\partial}{\partial\theta}\left(\frac{v_{\phi}}{\sin\theta}\right)$
(13)

(e.g., Stone et al., 1999; Fernández & Metzger, 2013; Inayoshi et al., 2019).
The strength of anomalous shear viscosity is calculated with the
$\alpha$-prescription (Shakura & Sunyaev, 1973),

${\tilde{\nu}}=\alpha\frac{c_{\rm s}^{2}}{\Omega_{\rm
K}}\cdot\exp\left(-\frac{|z|}{H}\right),$ (14)

where $\alpha$ is the viscous parameter, $c_{\rm s}$ is the sound speed,
$\Omega_{\rm K}\equiv(GM_{\bullet}/r^{3})^{1/2}$, and $H(\equiv c_{\rm
s}/\Omega_{\rm K})$ is the disk scale height. Note that the exponential factor
imposes that the viscous process is active near the midplane. The strength of
viscosity is set to

$\alpha=\alpha_{0}+\alpha_{\rm max}\exp\left[-\left(\frac{\rho_{\rm
crit}}{\rho}\right)^{2}\right].$ (15)

The first term corresponds to the strength of MRI turbulence and the value is
set to $\alpha_{0}=0.01$ (e.g., Zhu & Stone, 2018; Takasao et al., 2018). The
second term characterizes the torque caused by non-axisymmetric structure
(e.g., spiral arms) excited in a marginally unstable disk against its self-
gravity. The density threshold, above which viscosity turns on, is assumed to
be $\rho_{\rm crit}\equiv\Omega_{\rm K}^{2}/(\pi G)\simeq 2.2\times
10^{-21}(M_{\bullet}/10^{5}~{}M_{\odot})(r/10~{}{\rm pc})^{-3}~{}{\rm
g~{}cm}^{-3}$. This choice is motivated by the following reasons. The local
gravitational instability of a rotating disk is described by Toomre’s $Q$
parameter (Toomre, 1964) defined by

$Q\equiv\frac{c_{s}\kappa_{\Omega}}{\pi G\Sigma},$ (16)

where $\kappa_{\Omega}^{2}=4\Omega^{2}+d\Omega^{2}/d\ln r$ is the epicyclic
frequency and $\Sigma$ is the disk surface density. For a geometrically-thin
cold disk around a point mass, the $Q$-value is approximated to
$Q\simeq\Omega_{\rm K}^{2}/(2\pi G\rho)=0.5(\rho_{\rm crit}/\rho)$. Thus, the
second term of the right-hand-side of Eq. (15) is written as
$\propto\exp(-4Q^{2})$, which is a commonly used parametrization of the
effective viscosity adopted in semi-analytical models of a self-gravitating
disk; $\propto\exp(-Q^{A}/B)$ where $A>0$ and $B>0$ (Zhu et al. 2009;
Takahashi et al. 2013; see also Kratter & Lodato 2016). The value of
$\alpha_{\rm max}$ depends on the level of non-axisymmetric structures in a
disk. A previous 3D RHD simulation study of a dusty circum-nuclear disk around
a BH shows that the mass inflow velocity is as high as a substantial fraction
of the free-fall velocity due to strong torque caused by spiral-arms in the
disk, indicating $\alpha_{\rm max}\simeq O(1)$ (Toyouchi et al., 2021). We
here adopt $\alpha_{\rm max}=2$ (Hirano et al., 2014; Fukushima et al., 2020).
The critical $Q$-value for the onset of gravitational torque caused by spiral
arms in a disk is considered to be $Q\sim 1$, but the exact value depends on
various properties of the disk (e.g., cooling, heating, and disk irradiation).
We note that in our viscous model, the second term in Eq. (15) becomes larger
than $\alpha_{0}$ when $Q\leq 1.15$.

We consider cooling processes associated with chemical reactions and radiative
processes of metal-polluted gas, which is composed of 11 chemical species (H,
H2, e-, H+, H${}^{+}_{2}$, H-, He, He+, He++, C+, and O) and dust grains. We
solve the chemical reaction network of primordial gas among the 9 species
(hydrogen and helium) and adopt the same reaction rate coefficients compiled
in Li et al. (2021b). The cooling rates by collisional excitation and
ionization, radiative recombination of H, He, He+ atoms, their free-free
emission (Glover & Jappsen, 2007), and the fine-structure lines of C+ and O
(Hollenbach & McKee, 1989) are calculated in the optically-thin limit. The
cooling rate of H2 rovibrational transitions is included (Galli & Palla,
1998), but this effect is negligible under strong H2-dissociating photons
produced from bulge stars. We assume that the amount of heavy elements in gas-
phase and dust grains is proportional to the relative metallicity
$Z/Z_{\odot}$ and give the dust-to-gas mass ratio by $0.01~{}Z/Z_{\odot}$.

We implement photoionization of H, He, and He+ by intense BH radiation and
associated photoheating by solving radiation transfer. In addition, we
consider H2-dissociating radiation in the Lyman-Werner band and X-ray
radiation background produced from bulge stars surrounding the accreting BH
and its gravitational influence radius. As for X-ray irradiation, the
secondary ionization effect caused by energetic primary electrons are
considered (Inayoshi & Omukai, 2011, references therein). More detailed
treatments of BH and stellar radiation are described in §2.2.

###  2.2 Radiation transfer

We quantify the radiative heating rate $\Gamma$, ionization rate $\zeta$, and
radiation pressure force by solving the radiation transfer (RT) of both direct
and diffusion components. We adopt a hybrid RT scheme (e.g., Kuiper et al.,
2020): the direct component emitted from the central BH is solved by means of
the multi-frequency ray-tracing method, while the diffusive component owing to
thermal (re-)emission by dust grains is solved by a flux-limited diffusion
(FLD) approximation.

The direct component of the flux associated with BH accretion is the primary
source of radiative feedback that affects the properties of accreting flows.
The radiation flux can analytically be expressed as

$F_{\nu}(r,\theta,t)=\frac{\mathscr{L}_{\bullet,\nu}(\theta,t)}{4\pi
r^{2}}~{}e^{-(\tau_{\rm g,\nu}+\tau_{\rm d,\nu})},$ (17)

where $\nu$ is a frequency of radiation, $F_{\nu}(r,\theta,t)$ is the radial
component of the specific radiative flux at a distance of $r$ from the center
toward a direction of $\theta$ radian from the pole, and
$\mathscr{L}_{\bullet,\nu}/(4\pi r^{2})$ is the unattenuated specific
radiative flux that includes the anisotropic degree (see §2.3). The optical
depth of gas and dust to a photon at a frequency of $\nu$ are calculated as

$\displaystyle\tau_{\rm g,\nu}$ $\displaystyle=\int_{r_{\rm
in}}^{r}\sum_{i}n_{i}(\hat{r},\theta,t)\cdot\sigma_{i,\nu}~{}d\hat{r},$ (18)
$\displaystyle\tau_{\rm d,\nu}$ $\displaystyle=\int_{r_{\rm
in}}^{r}\rho(\hat{r},\theta,t)\cdot\kappa_{\rm d,\nu}~{}d\hat{r},$ (19)

respectively, where $n_{i}$ is the number density of the $i$-th component
($i=$H, He, and He+), $\sigma_{\nu,i}$ is the bound-free cross section (Verner
et al., 1996; Yan et al., 1998), and $\kappa_{\rm d,\nu}$ is the dust opacity.
In our simulations, we use the opacity table taken from Draine & Lee (1984) in
the RT calculation for non-ionizing radiation, while the attenuation level of
ionizing (EUV) radiation is calculated by using one single opacity value of
$\kappa_{\rm d,0}\simeq 2.8~{}[Z/(0.01~{}Z_{\odot})]~{}{\rm cm^{2}~{}g^{-1}}$
(per gas mass). The radiation flux absorbed by each component is used for
calculating the heating and ionization rates ($\Gamma$ and $\zeta$) owing to
EUV irradiation and the dust heating rate $S_{\bullet}$ caused by all the
photons emitted from the accreting BH. Those rates are estimated so that the
number of photons emitted along any line of sight equals the number of
photoionization in that direction over any time interval (e.g., Whalen &
Norman, 2006; Inayoshi et al., 2016; Nakatani et al., 2018a, b). We note that
the component of diffusive EUV radiation produced by radiative recombination
of gas is not considered because UV feedback onto accretion flows is dominated
by the direct component emitted from the central bright object (e.g., Hosokawa
et al., 2011). However, in order to calculate the size of an ionized region
precisely, we need to take into account that recombinations directly to the
ground state lead to immediate ionization of a hydrogen atom in the vicinity.
To treat this effect, we adopt the “on-the-spot” approximation using the
case-B recombination rate coefficient.

The diffusion component of infrared (IR) thermal radiation produced by dust
grains is treated with a grey-approximated FLD method by solving the moment
equation

$\frac{\partial E_{\rm rad}}{\partial t}+\nabla\cdot\mbox{\boldmath$F$}_{\rm
IR}=\rho\kappa_{\rm P}c(a_{R}T_{\rm d}^{4}-E_{\rm rad})+S_{\bullet},$ (20)

where $E_{\rm rad}$ is the radiation energy density, $T_{\rm d}$ is the dust
temperature, $\kappa_{\rm P}$ is the Planck-mean dust opacity, $a_{\rm R}$ is
the radiation constant, and $\mbox{\boldmath$F$}_{\rm IR}$ is the IR radiation
flux (re-)emitted from dust and is approximated as

$\mbox{\boldmath$F$}_{\rm IR}=-\frac{\lambda c}{\rho\kappa_{\rm R}}\nabla
E_{\rm rad},$ (21)

where $\kappa_{\rm R}$ is the Rosseland-mean dust opacity and $\lambda$ is the
flux limiter (Levermore & Pomraning, 1981). We solve the radiation transfer
equations with an implicit solver imposing a zero-gradient boundary condition
for the radiation energy. The dust temperature is calculated separately from
the radiation temperature by using the two-temperature linearization approach
(see more details in Kuiper et al., 2020).

With the direct component of radiation incident from the BH, we calculate the
outward radiation force through electron scattering, bound-free transitions,
and attenuation by dust grains. The diffusive component exerts the IR
radiation force to both the radial and polar directions.

###  2.3 BH radiation model

We inject photons from the unresolved central region (i.e., the sink cell) to
the computational domain, supposing that a compact accretion disk forms around
the nuclear BH. Following Sazonov et al. (2004), we consider a broad-band
radiation spectral energy distribution (SED) obtained from the averaged quasar
samples

$\displaystyle
L_{\bullet,\nu}=\begin{cases}L_{0}\left(\dfrac{\nu}{\nu_{0}}\right)^{-0.6}&(\nu_{\rm
min}\leq\nu<\nu_{0}),\vspace{10pt}\\\
L_{0}\left(\dfrac{\nu}{\nu_{0}}\right)^{-1.5}&(\nu_{0}\leq\nu\leq\nu_{\rm
max}),\end{cases}$ (22)

where $h\nu_{0}=10~{}{\rm eV}$, $h\nu_{\rm min}=1~{}{\rm eV}$, and $h\nu_{\rm
max}=1~{}{\rm keV}$. With this spectral shape, the bolometric (frequency
integrated) luminosity is calculated as $L_{\bullet}=\int_{\nu_{\rm
min}}^{\nu_{\rm max}}L_{\bullet,\nu}d\nu\simeq 3.3~{}\nu_{0}L_{0}$ and the
total EUV luminosity ($h\nu\geq h\nu_{\rm Ly}\equiv 13.6$ eV) is $L_{\rm
EUV}\simeq 1.5~{}\nu_{0}L_{0}$. The normalization of the luminosity is
determined as a function of the mass accretion rate onto the BH. We note that
the spectral power-law indices we adopt are consistent with those obtained
from the stacked UV (rest-frame wavelengths $600$-$2500~{}{\rm\AA}$) spectrum
of 53 luminous quasars at $z\simeq 2.4$ with a correction for the intervening
Lyman forest and Lyman continuum absorption (Lusso et al., 2015); namely, the
continuum spectrum is approximated by a broken power law with indices of
$\simeq-0.61\pm 0.01$ at $\nu<\nu_{0}$ and $\simeq-1.7\pm 0.61$ at
$\nu\geq\nu_{0}$, respectively.

We set a model for radiation luminosity emitted from the accreting BH as

$\displaystyle\frac{L_{\bullet}}{L_{\rm
Edd}}=\begin{cases}\eta_{0}\dot{m}&(\dot{m}<20),\vspace{10pt}\\\
2\left[1+\ln\left(\dfrac{\dot{m}}{20}\right)\right]&(\dot{m}\geq
20),\end{cases}$ (23)

where $L_{\rm Edd}$ is the Eddington luminosity, $\dot{m}$ is the BH mass
accretion rate normalized by the Eddington rate $\dot{M}_{\rm Edd}~{}[\equiv
L_{\rm Edd}/(\eta_{0}c^{2})]$, and $\eta_{0}=0.1$ is the radiative efficiency
in the sub-Eddington and mildly-Eddington regime (Shakura & Sunyaev, 1973).
The radiative luminosity at $\dot{m}>20$ is obtained from Watarai et al.
(2000) in a Slim-disk model (Abramowicz et al., 1988). We note that the two
expressions in Eq. (23) are broadly consistent with the results obtained with
general-relativistic, (magneto-) RHD simulations (Jiang et al. 2014; Sa̧dowski
et al. 2015; see also a model comparison in Inayoshi et al. 2020).

Furthermore, we assume isotropic/anisotropic radiation fields depending on the
bolometric luminosity emitted from the accreting BH. When the radiative
luminosity is lower than a critical value of $L_{\rm crit}$, we inject the
isotropic radiation flux from the center. On the other hand, when the BH is
rapidly growing and $L_{\bullet}>L_{\rm crit}$, the radiation flux with
$L_{\rm crit}$ is considered to be isotropic but the excess luminosity $\Delta
L_{\bullet}(\equiv L_{\bullet}-L_{\rm crit}$) is distributed into the
radiation flux anisotropically as $\Delta
L_{\bullet}(\mathcal{N}+1)|\cos^{\mathcal{N}}\theta|$, where $\mathcal{N}$
characterizes the anisotropic degree and is set to $\mathcal{N}=4$ in our
simulations (see also Ohsuga et al., 2005; Takeo et al., 2018). Therefore, the
anisotropic degree of the flux is characterized as
$\mathscr{L}_{\bullet}(\equiv\int\mathscr{L}_{\bullet,\nu}d\nu)=L_{\rm
crit}+\mathcal{F}(\theta)\Delta L_{\bullet}$ with a function of
$\mathcal{F}(\theta)$ defined by

$\displaystyle\mathcal{F}(\theta)=\begin{cases}(\mathcal{N}+1)|\cos^{\mathcal{N}}\theta|&(\Delta
L_{\bullet}>0),\vspace{5pt}\\\ 1.0&({\rm otherwise}).\end{cases}$ (24)

It is worth mentioning that this anisotropic model injects radiation to the
equatorial region ($\theta\simeq\pi/2$) and the non-zero radiation flux
through the equator affects the thermal properties of an accretion disk via IR
re-emission from hot dust grains. In addition, RHD simulations performed by
Ohsuga et al. (2005) suggest that a super-Eddington accreting disk transports
a significant fraction of radiation energy through the equator via diffusion
and the luminosity is comparable to the Eddington value (although most
radiation is collimated toward the poles). Therefore, we adopt $L_{\rm
crit}=L_{\rm Edd}$.

###  2.4 Stellar irradiation

Bulge stars are not only sources of external gravitational potential but also
radiation sources associated with star formation activity. We here consider LW
radiation and X-rays, which play an important role in determining the thermal
and chemical properties of gas around the nuclear region. The LW and X-ray
radiation is mainly produced by short-lived massive stars and X-ray binaries,
and their luminosities are proportional to the star formation rate (SFR). We
assume that the bulge is composed of metal-poor stars with
$Z=0.02~{}Z_{\odot}$ whose mass $m_{\star}$ follows a Salpeter initial mass
function ($1\leq m_{\star}/M_{\odot}\leq 100$), yielding the number flux of LW
photons $Q_{\rm LW}\simeq 5.5\times 10^{53}~{}{\rm s}^{-1}~{}({\rm
SFR}/M_{\odot}~{}{\rm yr}^{-1})$. Note that the number flux of LW radiation is
lowered only by $\simeq 15\%$ over $5\times 10^{-4}\leq Z/Z_{\odot}\leq 0.2$
(Schaerer, 2003). Taking the mean energy of LW bands ($h\nu_{\rm
LW}=12.4~{}{\rm eV}$), the luminosity is calculated as
$L_{\rm\star,LW}=h\nu_{\rm LW}Q_{\rm LW}\simeq 1.1\times 10^{43}~{}{\rm
erg~{}s}^{-1}({\rm SFR}/M_{\odot}~{}{\rm yr}^{-1})$. We also assume that the
X-ray luminosity correlates with the SFR as seen in local star forming
galaxies: $L_{\rm\star,X}\simeq 6.7\times 10^{39}~{}{\rm erg~{}s}^{-1}({\rm
SFR}/M_{\odot}~{}{\rm yr}^{-1})$ in the $2-10$ keV band (e.g., Grimm et al.,
2003; Mineo et al., 2014).

We treat the stellar radiation field in the optically-thin limit (note that
the radiation flux produced from the central accreting BH dominates the
stellar radiation flux at the vicinity of the BH within $1~{}{\rm pc}$).
Assuming a spherically symmetric distribution of bulge stars, the radiation
energy density at a distance of $r$ from the center is given by

$\displaystyle
cE_{\star,i}(r)=\int_{0}^{r}\epsilon_{i}(\hat{r})\frac{\hat{r}^{2}}{r^{2}}d\hat{r}+\int_{r}^{\infty}\epsilon_{i}(\hat{r})\frac{\hat{r}}{r}\ln\sqrt{\frac{\hat{r}+r}{\hat{r}-r}}d\hat{r},$
(25)

where the first and second terms in the right hand side represent
contributions from inside and outside the distance of $r$, $\epsilon_{i}(r)$
is the radiation emissivity at each cell

$\displaystyle\epsilon_{i}(r)=\frac{L_{\star,i}}{2\pi}\cdot\frac{r_{\rm
c}}{r(r+r_{\rm c})^{3}}\hskip 10.0pt{\rm for}~{}i={\rm LW~{}and~{}X},$ (26)

where the functional form is given by the stellar mass distribution, and
$\rho_{\star}(r)=\Delta\Phi_{\star}(r)/(4\pi G)$. Integrating Eq. (25), we
obtain

$cE_{\star,i}(r)=\frac{L_{\star,i}}{4\pi r_{\rm c}^{2}}~{}\Xi(r/r_{\rm c}),$
(27)

where

$\displaystyle\Xi(x)\equiv\frac{1}{(1+x)^{2}}+\int_{x}^{\infty}\frac{\ln(\frac{x_{0}+x}{x_{0}-x})}{(1+x_{0})^{3}}\frac{dx_{0}}{x}.$
(28)

The function of $\Xi(x)$ is approximated as $O(1)$ at $x<1$ and decreases as
$\propto x^{-2}$ at $x\gg 1$.

In this paper, we consider H2 photodissociation by LW radiation (${\rm
H}_{2}+\gamma_{\rm LW}\rightarrow 2~{}{\rm H}$), H- photodetachment (${\rm
H}^{-}+\gamma_{\rm 1eV}\rightarrow{\rm H}+{\rm e}^{-}$) by photons with
$h\nu\simeq 1~{}{\rm eV}$, and ionization/heating of H and He by X-rays with
$h\nu=2$ keV. Those reaction rates are calculated in the same way as in
Inayoshi & Tanaka (2015). For the H- photodetachment rate, we set the minimum
rate coefficient at a range of the effective temperature of bulge stars
($10^{4}~{}{\rm K}\leq T_{\rm eff}\leq 10^{5}~{}{\rm K}$)222 While the H2
photodissociation rate hardly depends on $T_{\rm eff}$, the H- photodetachment
rate increases with lower $T_{\rm eff}$ because photons with $\simeq 1$ eV
($\ll h\nu_{\rm LW}$) contributes to the rate. With a realistic stellar
spectrum for metal-poor galaxies (Inoue, 2011), the H- photodetachment rate
boosts by a factor of $\sim 20$ from that for $T_{\rm eff}=10^{5}~{}{\rm K}$..
We also take into account secondary ionization by energetic electrons produced
by primary X-ray ionization (Shull & van Steenberg, 1985).

Table 1: Models considered

Model |  $M_{\star}$ ($M_{\odot}$) |  $T_{\rm vir}$ (${\rm K}$) | $\epsilon_{\star}$ |  $n_{\rm c}$ (${\rm cm}^{-3}$) | transition | References
---|---|---|---|---|---|---
Massive Halos${\dagger}$ |  |  |  |  |  |
B7T5N4 | $10^{7~{}~{}}$ | $10^{5}$ | 0.05 | $10^{4~{}~{}}$ | Y | §3.2.1
B7T5N3.5 | $10^{7~{}~{}}$ | $10^{5}$ | 0.05 | $10^{3.5}$ | Y | §3.2.1
B7T5N3 (fiducial) | $10^{7~{}~{}}$ | $10^{5}$ | 0.05 | $10^{3~{}~{}}$ | Y | §3.1
B7T5N2.5 | $10^{7~{}~{}}$ | $10^{5}$ | 0.05 | $10^{2.5}$ | Y | §3.2.1
B7T5N2 | $10^{7~{}~{}}$ | $10^{5}$ | 0.05 | $10^{2.5}$ | N | §3.2.1
B0T5N3 | $0~{}~{}~{}~{}$ | $10^{5}$ | 0.05 | $10^{3~{}~{}}$ | N | §3.2.2
B5T5N3 | $10^{5~{}~{}}$ | $10^{5}$ | 0.05 | $10^{3~{}~{}}$ | N | §3.2.2
B6T5N3 | $10^{6~{}~{}}$ | $10^{5}$ | 0.05 | $10^{3~{}~{}}$ | N | §3.2.2
B6.5T5N3 | $10^{6.5}$ | $10^{5}$ | 0.05 | $10^{3~{}~{}}$ | N | §3.2.2
B7T5N3-highSFE | $10^{7~{}~{}}$ | $10^{5}$ | 0.5 | $10^{3~{}~{}}$ | Y | §3.2.3
Normal Halos${\ddagger}$ |  |  |  |  |  |
B7T4N4 | $10^{7~{}~{}}$ | $10^{4}$ | 0.05 | $10^{4~{}~{}}$ | N | §3.2.4
B7T4N3 | $10^{7~{}~{}}$ | $10^{4}$ | 0.05 | $10^{3~{}~{}}$ | N | §3.2.4
B6T4N3 | $10^{7~{}~{}}$ | $10^{4}$ | 0.05 | $10^{3~{}~{}}$ | N | §3.2.4

Note. — Simulation runs and input parameters. Column (1) model ID, (2) initial
bulge mass, (3) virial temperature of the host DM halo, (4) star-formation
efficiency, (5) initial gas density at the center, and (6) Y(/N) indicates the
cases where a transition to rapid accretion does (not) occur by the simulation
termination at $t=4.0$ Myr (high-SFE case) and $t=3.5$ Myr (otherwise).

${\dagger}$: a massive dark-matter halo with $M_{\rm h}=2\times
10^{9}~{}M_{\odot}$ at $z=15$.

${\ddagger}$: a normal atomic-cooling halo with $M_{\rm h}=10^{8}~{}M_{\odot}$
at $z=10$.

###  2.5 Boundary and initial conditions

We set a computational domain of $r_{\rm min}\leq r\leq r_{\rm max}$ and
$\epsilon\leq\theta\leq\pi-\epsilon$, where $\epsilon=0.001$ radian is set to
avoid numerical singularity at poles. We adopt logarithmically spaced grids in
the radial direction and uniformly spaced grids in the polar direction. The
number of the grid points is set to ($N_{r},N_{\theta})=(250,240)$. For all
the simulations we conduct, the minimum and maximum radius are set to $r_{\rm
min}=0.1~{}{\rm pc}$ and $r_{\rm max}=100~{}{\rm pc}$.

We adopt the outflow boundary conditions at the innermost and outermost grid,
where zero gradients across the boundaries are imposed on physical quantities
in order to avoid spurious reflection of wave energy at the boundaries. We
also impose $v_{r}\leq 0$ at the inner boundary (i.e., inflowing gas from
ghost cells is prohibited). At the poles, the reflective condition is imposed
on the polar component of the velocity $v_{\theta}$. We set a density floor of
$n_{\rm min}=1~{}{\rm cm}^{-3}$ and a maximum radial velocity of $v_{\rm
max}=150~{}{\rm km~{}s}^{-1}$. This treatment allows us to avoid cells with a
high Mach number ($\gg 100$) in the outflowing regions near the inner and
polar boundary.

Figure 1: A schematic picture of accretion flows onto a massive BH embedded in
a protogalactic bulge with a mass of $M_{\star}$ and DM halo with a virial
temperature of $T_{\rm vir}$. There are four characteristic scales: the core
radius of the initial density distribution ($r_{0}$), the BH gravitational
influence radius ($r_{\rm B}$), the bulge core size ($r_{\rm c}$), and the
size of ionized regions created by the radiating BH ($r_{\rm ion}$). The
computational domain covers $r_{\rm min}\leq r\leq r_{\rm max}$, where those
characteristic radii are well resolved.

As initial conditions, we consider a neutral gas cloud with a temperature of
$T=10^{4}~{}{\rm K}$ and a spherically symmetric density profile with a flat
core with the central density of $n_{\rm c}$ and an envelope following
$\propto r^{-2}$,

$\displaystyle n(r)=\frac{n_{\rm c}}{1+(r/r_{0})^{2}},$ (29)

where the core radius $r_{0}$ is set so that the density profile becomes
consistent with that of high-redshift protogalaxies obtained by cosmological
simulations at $r\gg r_{0}$. Estimating the gas density at the outer boundary
of the halo by the minimum-energy truncated isothermal sphere model as $n_{\rm
vir}\simeq 0.07~{}{\rm cm}^{-3}[(1+z)/21]^{3}$ (Shapiro et al., 1999; Iliev &
Shapiro, 2001), the density profile outside the core is approximated as
$n(r)\simeq f_{n}n_{\rm vir}(r/r_{\rm vir})^{-2}$ or equivalently as

$\displaystyle n(r)$ $\displaystyle\simeq 2.5\times 10^{2}~{}f_{n}T_{\rm
vir,4}\left(\frac{r}{10~{}{\rm pc}}\right)^{-2}~{}{\rm cm}^{-3},$ (30)

and thus the core radius is approximated as

$\displaystyle r_{0}\simeq 5~{}f_{n}^{1/2}T_{\rm vir,4}^{1/2}n_{\rm
c,3}^{-1/2}~{}{\rm pc},$ (31)

where $T_{\rm vir,4}=T_{\rm vir}/(10^{4}~{}{\rm K})$, $n_{\rm c,3}=n_{\rm
c}/(10^{3}~{}{\rm cm}^{-3})$, and $f_{n}\sim O(1)$ is a numerical factor. We
note that the density profile with $f_{n}\simeq 4$ agrees with those of atomic
cooling haloes with $T_{\rm vir}\simeq 10^{4}~{}{\rm K}$ obtained by
cosmological simulations without star formation and stellar feedback (Wise et
al., 2008; Shang et al., 2010; Regan et al., 2014). After a seed BH forms at
the center of the halo, the density normalization becomes consistent with
$f_{n}\simeq 1$.

In addition, star formation in the halo is modeled by assuming a conversion
efficiency from gas into stars $\epsilon_{\star}$. In our fiducial case, we
adopt $\epsilon_{\star}=0.05$ (e.g., Visbal et al., 2015), which is motivated
by abundance matching and the observed UV luminosity function of galaxies at
$z\simeq 6$ (Bouwens et al., 2015). Note that the efficiency is calculated as
the average value over time and scales in galaxies. Alternatively, we suppose
that a star-formation episode with a high value of $\epsilon_{\star}$ lasts
within a few Myrs in the nuclear region before SN explosions of massive stars
begin to occur and regulate the star formation efficiency in the bulge.
According to numerical simulations that study star cluster formation from a
giant molecular cloud with a size of $10-100~{}{\rm pc}$, the star formation
efficiency (SFE) becomes as high as $\epsilon_{\star}\gtrsim 0.2-0.3$ when the
initial gas surface density is higher than $\sim 10^{3}~{}M_{\odot}{\rm
pc}^{-2}$ (Fukushima et al., 2020; Fukushima & Yajima, 2021). Referring to
those observational and theoretical studies, we assume
$f_{n}=f_{0}(1-\epsilon_{\star})$, where $f_{0}=1$ and $\epsilon_{\star}=0.05$
(fiducial case) and $f_{0}=4$ and $\epsilon_{\star}=0.5$ (high SFE case).
Following the definition of the star formation efficiency $\epsilon_{\star}$,
the SFR is approximated as

$\displaystyle{\rm SFR}$
$\displaystyle=\frac{4\pi\epsilon_{\star}}{1-\epsilon_{\star}}\rho(r)r^{2}V_{\rm
c},$ $\displaystyle\simeq 0.11~{}\epsilon_{\star}f_{0}~{}T_{\rm
vir,4}^{3/2}~{}M_{\odot}~{}{\rm yr}^{-1},$ (32)

where $V_{\rm c}$ is the halo circular velocity. This SFR is used to estimate
the emissivity of stellar irradiation (LW and X-rays) and to calculate the
bulge growth in mass. Note that in most cases, the gas mass within $r\simeq
r_{\rm c}$ for the given initial condition is lower than the total mass of
newly forming stars333Although we do not consider the self-gravity of gas in
our simulations, it would dominate the bulge gravity at the intermediate scale
of $\sim O(10~{}{\rm pc})$ only in the early stage of the bulge formation
where $M_{\star}\lesssim 10^{6}~{}M_{\odot}$. In this case, however, radiative
feedback associated with BH accretion blows the gas outward and then the gas
self-gravity becomes less important eventually (see §3.2.2).. Therefore, our
star-formation model implicitly assumes that the bulge growth is not led by
in-situ star formation but efficient migration of stars formed at larger
radii.

The initial velocity field is set to $(v_{r},v_{\theta},v_{\phi})=(0,0,j_{\rm
in}/R)$, where $R=r\sin\theta$ is the cylindrical radius. We assume that the
initial specific angular momentum $j_{\rm in}(R)$ is proportional to the
Keplerian velocity measured with the enclosed gas mass $\mathcal{M}$ within a
distance of $R$. Namely, we adopt $j_{\rm in}(R)=j_{0}\sqrt{G\mathcal{M}R}$
and set $j_{0}=0.3$. The model assumption is motivated by the fact that the
dynamics of a collapsing gas obeys a self-similar solution and the rotational
velocity in the accreting envelope where $\rho\propto r^{-2}$ is a good
fraction of the Keplerian velocity in agreement with hydrodynamical
simulations of collapsing gas in a DM halo (Abel et al., 2002; Yoshida et al.,
2008; Inayoshi et al., 2014).

###  2.6 Cases considered

In Table 1, we summarize the models we consider in this paper. We consider two
types of DM halo with different virial temperatures of $T_{\rm vir}\simeq
10^{5}~{}{\rm K}$ ($M_{\rm h}=2\times 10^{9}~{}M_{\odot}$ and $z=15$) and
$T_{\rm vir}\simeq 10^{4}~{}{\rm K}$ ($M_{\rm h}=10^{8}~{}M_{\odot}$ and
$z=10$). Note that the mean molecular weight is set to $\mu=0.5$ for
evaluating the virial temperature above $T_{\rm vir}=10^{4}~{}{\rm K}$. The
two cases are referred to “Massive Halo” model and “Normal Halo” model,
respectively. In the former case, we consider a relatively massive DM halo
because high-redshift quasar host galaxies form in rare, overdense regions of
the universe at $z>6$ and their progenitor halos would be likely as massive as
our choice at higher redshifts when seed BHs form (Li et al., 2021b). The
latter case corresponds to that of typical first galaxies where star formation
is trigged by hydrogen atomic cooling (Bromm & Yoshida, 2011).

We perform the 13 simulations shown in Table 1 with different bulge mass
($M_{\star}=0$, $10^{5}$, $10^{6}$, $10^{6.5}$, and $10^{7}~{}M_{\odot}$),
initial gas density at the center ($n_{\rm c}=10^{2}$, $10^{2.5}$, $10^{3}$,
$10^{3.5}$, and $10^{4}~{}{\rm cm}^{-3}$), and halo virial temperature
($T_{\rm vir}=10^{4}$ and $10^{5}~{}{\rm K}$). In all the simulations, we
adopt the initial BH mass to $M_{\bullet,0}=10^{5}~{}M_{\odot}$. The SFE and
the initial density normalization are assumed to be $\epsilon_{\star}=0.05$
and $f_{0}=1$, respectively, for most of the cases. To see the impact of their
parameter choice, we additionally consider a case (B7T5N3-highSFE), where a
higher SFE of $\epsilon_{\star}=0.5$ and $f_{0}=4$ are set.

In Fig. 1, we summarize the characteristic physical scales that determine the
properties of accretion flows onto a seed BH embedded in a protogalaxy. In our
case, the typical values of the BH gravitational influence radius (the so-
called Bondi radius) for neutral gas, ($\mu=1.22$) defined by

$r_{\rm B}\equiv\frac{GM_{\bullet}}{c_{\rm s}^{2}}\simeq
6.4~{}\left(\frac{M_{\bullet}}{10^{5}~{}M_{\odot}}\right)\left(\frac{T}{10^{4}~{}{\rm
K}}\right)^{-1}~{}{\rm pc},$ (33)

is well resolved (i.e., $r_{\rm min}<r_{\rm B}<r_{\rm max}$). This indicates
the distance from the BH within which the gravitational energy dominates over
the thermal energy of the gas and thus gas accretion begins to occur unless BH
feedback plays an important role. In reality, however, when the BH grows via
mass accretion, the accreting flow releases its gravitational energy as
radiation, which heats the surrounding gas and forms an ionized bubble with a
size of $r_{\rm ion}$. In the intermediate region, stars form a galactic bulge
component with a total mass of $M_{\star}$ and a half-mass radius of $R_{\rm
e}=(1+\sqrt{2})r_{\rm c}\simeq 16~{}{\rm
pc}~{}(M_{\star}/10^{7}~{}M_{\odot})^{0.47}$. For both the “Massive” and
“Normal” halo cases, the halo virial radius is $r_{\rm vir}\simeq 1.6~{}{\rm
kpc}$ and $0.9~{}{\rm kpc}$, respectively. Therefore, our computational domain
covers the interior of the host DM halo ($r_{\rm max}\simeq 0.1~{}r_{\rm
vir}$).

##  3 Results

###  3.1 Fiducial case

Figure 2: Time evolution of the gas accretion rate onto a massive BH with an
initial mass of $M_{\bullet}=10^{5}~{}M_{\odot}$ embedded in a dense gas cloud
with the central density of $n_{\rm c}=10^{3}~{}{\rm cm}^{-3}$ hosted in a
stellar bulge with $M_{\star}=10^{7}~{}M_{\odot}$ and a massive DM halo with a
virial temperature of $T_{\rm vir}\simeq 10^{5}~{}M_{\odot}$ (model: B7T5N3).
The dashed curves show the Eddington accretion rate $\dot{M}_{\rm Edd}$ and
the critical rate for the onset of rapid mass accretion ($\sim
500~{}\dot{M}_{\rm Edd}$) obtained in previous RHD simulations (Inayoshi et
al., 2016). Open circles mark the six epochs at which we show the density and
temperature distributions in Fig. 3 and radial profiles of the physical
quantities in Fig. 4.  Figure 3: Distribution of the gas density for the
fiducial case at the six elapsed times denoted by open circles in Fig. 2. In
each panel, the location of the ionization front where the neutral fraction is
$x_{\rm HI}=0.95$ (thin contours) and the velocity vectors are overlaid. In
the first two panels of (a) and (b), the velocity vectors with
$|\mbox{\boldmath$v$}|\geq 30~{}{\rm km~{}s}^{-1}$ are shown. At the early
stage ($t<1.6$ Myr; phases $a$-$c$), the intense inflows of neutral gas lead
to collapse of the ionized region and form a dense gaseous accretion disk
within the BH influence radius. In the late stage ($t>1.6$ Myr; phases
$d$-$f$), the disk feeds the BH at rates of $\gg\dot{M}_{\rm Edd}$ and the
bipolar outflows driven by radiation decreases the gas supplying rate from
larger radii.

We first discuss the fiducial case where a BH with
$M_{\bullet}=10^{5}~{}M_{\odot}$ is embedded in a dense gas cloud with the
central density of $n_{\rm c}=10^{3}~{}{\rm cm}^{-3}$ hosted in a stellar
bulge with an initial mass of $M_{\star}=10^{7}~{}M_{\odot}$ hosted in a
massive DM halo with $M_{\rm h}=2\times 10^{9}~{}M_{\odot}$ at $z=15$ (the
corresponding virial temperature is $T_{\rm vir}\simeq 10^{5}~{}{\rm K}$). The
star formation efficiency is set to $\epsilon_{\star}=0.05$.

Fig. 2 shows the time evolution of the accretion rate $\dot{M}_{\bullet}$
(solid curve). As reference values, the Eddington accretion rate $\dot{M}_{\rm
Edd}$ and the critical rate for the onset of hyper-Eddington accretion
($\simeq 500~{}\dot{M}_{\rm Edd}$; Inayoshi et al. 2016) are shown with dashed
curves. At the early stage of $t<1.6$ Myr, the mass accretion proceeds
episodically. While the accretion rate exceeds the Eddington value during
those burst phases, the long-term average rate is $\simeq\dot{M}_{\rm Edd}$
due to radiative feedback associated with BH feeding. At $t\simeq 1.6$ Myr,
the mass accretion rate abruptly rises and reaches $\dot{M}_{\bullet}\simeq
10~{}M_{\odot}~{}{\rm yr}^{-1}$, which corresponds to $\simeq 2.5\times
10^{3}~{}\dot{M}_{\rm Edd}$. After the peak, the rate gradually decreases with
time but the super-Eddington accretion phase lasts until $t\simeq 3$ Myr. Note
that the BH mass is $\simeq 6.3\times 10^{5}~{}M_{\odot}$ at the end of the
simulation. The overall behavior of the accretion rate is consistent with
those found in previous RHD simulations where the BH is not embedded in the
external gravitational potential composed of stars and DM (Inayoshi et al.,
2016; Takeo et al., 2018, 2020; Park et al., 2020; Toyouchi et al., 2021).

Figure 4: Radial structure of the gas density (top), temperature (middle),
velocity (bottom) along the equator ($\theta=90^{\circ}$; left-hand panels)
and pole ($\theta=0^{\circ}$; right-hand panels). The bottom panels present
the rotational velocity along the equator and radial velocity along the polar
direction, respectively. In each panel, we show the profiles at different
epochs during and after the transition: $t=0.28~{}{\rm Myr}$ (phase $a$;
purple), $0.35~{}{\rm Myr}$ (phase $b$; green), $1.5~{}{\rm Myr}$ (phase $c$;
blue), $1.65~{}{\rm Myr}$ (phase $d$; yellow), and $1.8~{}{\rm Myr}$ (phase
$e$; black). In the top-left and bottom-left panels, we overlay the critical
density (see the viscous model below Eq. 15) and the Keplerian velocity for
$M_{\bullet}=6\times 10^{5}~{}M_{\odot}$.

Fig. 3 presents the distribution of the gas density in the domain of $0\leq
x\leq 8~{}{\rm pc}$ and $|z|\leq 6~{}{\rm pc}$ at six different elapsed times
(phases $a$-$f$ of Fig. 2) . In each panel, the location of the ionization
front where the neutral fraction is $x_{\rm HI}=0.95$ (thin contours) and the
velocity vectors are overlaid. Fig. 4 presents the radial profiles of the gas
density (top), temperature (middle), and velocity (bottom) along the equator
($\theta=90^{\circ}$; left panels) and the pole ($\theta=0^{\circ}$; right
panels), respectively. In the bottom panels, we show the rotational velocity
along the equator and radial velocity along the polar direction, respectively.
Note that the radial profiles at the phase ($f$) are not shown because they
are similar to those at the phase ($e$) except that the shock front of the
polar outflow moves forward. This indicates that the disk properties are in a
quasi-steady state although the net accretion rate onto the BH gradually
decreases.

At the early stage of $t<0.3$ Myr, the accreting BH emits nearly isotropic
radiation and creates an ionizing bubble. The prolate shape of the bubble is
caused by angular momentum of the inflowing gas with a lower density near the
pole and by anisotropic radiation produced when the BH feeding rate exceeds
the Eddington value in a short period. While the size of the ionizing bubble
is initially maximized at $r_{\rm ion}\sim 30~{}{\rm pc}$, it shrinks to $\sim
3~{}{\rm pc}$ by the phase ($a$) owing to efficient radiative recombination of
the gas surrounding the Eddington-limited radiating BH. As a result, the
equatorial inflow of neutral gas penetrates into the interior of the BH
gravitational influence radius of $r_{\rm B}\simeq 8~{}{\rm pc}$ for neutral
gas with $T\simeq 8000~{}{\rm K}$ (see the left panels in Fig. 4). The inflow
with strong ram pressure leads to fragmentation of the shell of the ionized
bubble. Previous studies also found that collapse of an ionized region
surrounding the accreting BH is an essential process to trigger the accretion
transition (Inayoshi et al., 2016; Sakurai et al., 2016; Toyouchi et al.,
2019, 2021). However, we note that gas rotation, anisotropic radiation, and
dust obscuration moderate the propagation of ionizing radiation to the
equatorial region and thus do promote the transition (see also Takeo et al.,
2018, 2020).

Figure 5: Time-evolution of the BH accretion rate (left) and the mass accreted
onto the BH (right) with different values of the central density for the
initial conditions; $n_{\rm c}=10^{4}~{}{\rm cm}^{-3}$ (cyan),
$10^{3.5}~{}{\rm cm}^{-3}$ (orange), $10^{3}~{}{\rm cm}^{-3}$ (purple;
fiducial case), and $10^{2.5}~{}{\rm cm}^{-3}$ (green). The transition of
rapid accretion tends to occur earlier with higher values of $n_{\rm c}$. In
spite of the stochastic nature of mass accretion, the mass accreted onto the
BH reaches similar values of $\Delta M_{\bullet}\simeq(3-6)\times
10^{5}~{}M_{\odot}$, which is substantially higher than that expected under
the Eddington growth (dashed curve).

At $t\simeq 0.35$ Myr (phase $b$), intense inflows of neutral gas feeds the
nuclear scale within $\sim 2~{}{\rm pc}$, but the BH feeding rate does not
increase at the same moment. This is primarily because the inflowing gas with
angular momentum forms a rotationally-supported disk and the rotational
velocity exceeds the Keplerian value at $r\sim 2~{}{\rm pc}$ (see the left-
bottom panel in Fig. 4), where the inflow speed through the disk slows down
significantly. In addition, the outward radiation force exerted through
electron scattering, absorption of UV radiation by atomic bound-free
transitions and dust grain, and re-emission of IR radiation from heated dust
prevent the inflowing gas from feeding the BH at rates of $\gtrsim\dot{M}_{\rm
Edd}$ (Toyouchi et al., 2019). Meanwhile, the ionized regions become as small
as $r_{\rm ion}\simeq 0.3~{}{\rm pc}$ and the inflow from larger radii
accumulates mass in the nuclear region (phase $c$; see Fig. 3 and the left
panels in Fig. 4).

At the late stages of $t\gtrsim 1.6$ Myr (phase $d$), the disk becomes opaque
to UV by dust, and the inner-edge of the dense dusty disk reaches the central
cell with $r_{\rm min}=0.1~{}{\rm pc}$. Meanwhile, the mid-plane density
approaches the profiles of $\rho_{\rm crit}=\Omega_{\rm K}^{2}/(\pi G)$ (the
dashed line in the top-left panel of Fig. 4), where the Toomre’s $Q$ parameter
is nearly set to unity within the disk. Since the disk feeds the BH at rates
exceeding the Eddington value, the radiation flux from the center blows the
gas away and launches strong outflows toward the polar regions at velocities
of $v_{\rm out}\simeq 150~{}{\rm km~{}s}^{-1}$, which is significantly faster
than the escape velocity from the system (see Fig. 3). In the phase ($e$), the
radiation-driven outflow collides with the inflowing gas and creates bipolar
low-density cavities surrounded by a dense shell (see the right panels in Fig.
4). As a result, the emergent UV radiation begins to heat the disk surface and
drives thermally-driven outflows. Although the radiation/mechanical feedback
effect suppresses the mass supply from larger radii, the dense accretion disk
coexists with outflows and feeds the central BH at super-Eddington rates
$\gtrsim\dot{M}_{\rm Edd}$ until $t\simeq 3.5$ Myr (phase $f$), when the
simulation terminates.

###  3.2 Parameter dependence

Next, we investigate the dependence of the result on the model parameters: (1)
the central density $n_{\rm c}$, (2) the bulge mass $M_{\star}$, (3) the star
formation efficiency $\epsilon_{\star}$, and (4) the halo virial temperature
$T_{\rm vir}$.

####  3.2.1 Central density $n_{\rm c}$

In Fig. 5, we present the time evolution of the BH accretion rate (left) and
the mass accreted onto the BH (right) with different values of the central
density for the initial conditions. The transition to rapid accretion tends to
occur earlier with higher values of $n_{\rm c}$. This is because the size of
the ionized region surrounding the BH becomes smaller due to efficient
recombination and thus collapse of the compact ionized region triggers the BH
feeding in a shorter dynamical timescale. Although this trend holds for
$n_{\rm c}\geq 10^{3}~{}{\rm cm}^{-3}$, the accretion transition for the
lowest density case ($n_{\rm c}=10^{2.5}~{}{\rm cm}^{-3}$) occurs earlier than
that for the fiducial case due to the stochastic nature of mass accretion
through filamentary structures and a compact disk at $\sim 1~{}{\rm pc}$ (see
the phase $b$ in Fig. 3). Moreover, the two cases ($n_{\rm c}=10^{2.5}$ and
$10^{3.5}~{}{\rm cm}^{-3}$) show accretion bursts twice by the end of the
simulations. In those cases, the vertical oscillation of the accretion disk
after its formation promotes mass loading into radiation-driven outflows and
the net accretion rate through the disk decreases quickly with time. However,
as a result of weak radiative feedback, continuous mass supply from large
radii triggers the second accretion burst within $\sim 1$ Myr. In spite of the
stochastic nature of mass accretion, the gas mass accreted onto the BH reaches
similar values of $\Delta M_{\bullet}\simeq(3-6)\times 10^{5}~{}M_{\odot}$,
which is substantially higher than that expected under the Eddington growth
(dashed curve). We note that the accreted mass is comparable to that contained
within the gaseous core in the initial state ($r\lesssim r_{0}$).

It is worth mentioning that all the cases in Fig. 5 experience the transitions
to rapid accretion exceeding the Eddington rate. As described in Fig. 4, the
transition is induced by continuous accretion of neutral gas whose inward ram
pressure pushes the ionization front (hereafter, I-front). In §4, we quantify
the conditions required for the onset of rapid mass accretion. On the
contrary, when the central density is lower than $n_{\rm c}\sim 10^{2}~{}{\rm
cm}^{-3}$ (model: B7T5N2), the ionized region quickly reaches the core radius
$r_{0}$. Since the density follows $\rho\propto r^{-2}$ outside the core, the
I-front further expands at speeds faster than the sound speed of ionized gas
without disturbing the density structure (the so-called R-type I-front) and
the entire cloud is ionized (e.g., Kitayama et al., 2004; Whalen & Norman,
2006). In this case, therefore the transition to rapid mass accretion does not
occur.

####  3.2.2 Bulge mass $M_{\star}$

Figure 6: Time-evolution of the mass accretion rate onto the BH with different
values of the bulge mass; $M_{\star}=10^{7}~{}M_{\odot}$ (purple),
$10^{6.5}~{}M_{\odot}$ (green), $10^{6}~{}M_{\odot}$ (cyan),
$10^{5}~{}M_{\odot}$ (orange), and no bulge (black). The transition to rapid
accretion phases occur when the bulge mass is as high as
$M_{\star}=10^{7}~{}M_{\odot}$, which corresponds to $\sim 100~{}M_{\bullet}$
(see also the discussion in §4).

Fig. 6 shows the time evolution of the mass accretion rate with four different
bulge masses at $0\leq M_{\star}/M_{\odot}\leq 10^{7}$. In Fig. 7, we also
present the density distribution for the two cases of $M_{\star}=0$ (no bulge)
and $10^{6.5}~{}M_{\odot}$ around the epochs when the simulations terminate,
respectively.

Without the bulge component, the accretion rate shows multiple bursts with a
period of $\simeq 0.6$ Myr and the time-averaged rate is significantly below
the Eddington value; namely $\langle\dot{M}_{\bullet}\rangle\simeq
0.2~{}\dot{M}_{\rm Edd}$. In this case, the I-front initially propagates up to
$\sim 25~{}{\rm pc}$ and heats the ambient gas outside the BH influence radius
$r_{\rm B}\simeq 5~{}{\rm pc}$. Within the hot ionized gas, the BH’s gravity
accelerates the inflow within a new sonic point in the hot region ($\simeq
0.4~{}{\rm pc}$ for $T\simeq 10^{5}~{}{\rm K}$), while the gas pressure (and
partially the radiation pressure force) pushes the gas outwards at $r\gtrsim
1~{}{\rm pc}$ and the mass accretion rate is suppressed. In the quiescent
phases, although the gas is in thermal-pressure equilibrium against the BH
gravity, the mass depletion owing to BH feeding reduces the outward pressure
gradient force. As a result, the dense shell surrounding the ionized gas
accretes to the BH and leads to burst-like accretion again. The overall
behavior of mass accretion is consistent with previous RHD simulations where
the BH is assumed to be embedded in a uniform density distribution (e.g.,
Ciotti & Ostriker, 2001; Milosavljević et al., 2009; Park & Ricotti, 2011,
2012; Inayoshi et al., 2016). In our case, where the gas density decreases
outward, the I-front continues to expand during the multiple episodes of
burst-like accretion and reaches $\gtrsim 40~{}{\rm pc}$ by $t\simeq 3.5$ Myr.

Figure 7: Distribution of the gas density for the case without bulge stars
(left) and with the bulge of $M_{\star}=10^{6.5}~{}M_{\odot}$ (right) at the
end of those simulations. The arrows present the velocity vectors. In both
cases, the accreting BH is surrounded by an expanding ionized region and thus
rapid inflows of neutral gas onto the central BH are not observed.

With the bulge mass increasing, the episodic behavior of mass accretion ceases
and the time-average rate increases because the gas density and mass inflow
rate are enhanced within the additional bulge gravitational potential. When
the bulge mass is lower than $M_{\star}<10^{7}~{}M_{\odot}$, however, the
ionized region continues to expand and heats the surrounding gas in the same
way as in the case without a bulge. Note that this statement is still valid if
the self-gravity of gas would be taken into account in our simulations, since
the total gas mass in the domain is lower than $10^{7}~{}M_{\odot}$. As a
result, the transition to rapid mass accretion does not occur by the end of
the simulation. Since the free-fall velocity at $r\simeq 50~{}{\rm pc}$ is
$\sim 20~{}{\rm km~{}s}^{-1}$ (the enclosed mass is dominated by the bulge),
the transition could not occur within another $\sim 2.5$ Myr even if all the
radiative output turned off.

In summary, the size of the expanding I-front determines the nature of mass
accretion onto the BH, and the formation/acceleration of neutral-gas inflows
from larger radii is a key process to trigger the accretion transition. Based
on those findings, we give a simple analytical argument for the conditions
required for the transition in §4.

####  3.2.3 Star formation efficiency $\epsilon_{\star}$

Fig. 8 shows the time evolution of the mass accretion rate (solid), the
Eddington accretion rate (dotted), and the BH-to-bulge mass ratio of
$M_{\bullet}/M_{\star}$ (dashed) for the case with a high SFE
($\epsilon_{\star}=0.5$; B7T5N3-highSFE). The overall behavior of the
accretion rate is consistent with the cases with the lower SFE as shown in
Fig. 5. Compared to the fiducial case, the transition to rapid accretion
occurs earlier because a larger amount of gas is accumulated within the
gravitational influence radius of the fast growing bulge. With the high SFR,
intense stellar radiation heats the gas surrounding the BH but efficient
radiative cooling keeps the gas neutral. As a result, the attractive
gravitational force by the massive bulge is more important than the negative
feedback effect caused by the stellar irradiation. In addition, the BH-to-
bulge mass ratio evolves substantially within $\sim 3$ Myr and the ratio
reaches $\sim 0.04$ at the end of the simulation. This mass ratio is $\sim 10$
times higher than the BH-galaxy correlation seen in the local universe
(Kormendy & Ho, 2013) and is consistent with those of bright quasars at $z>6$
(see more discussion in §5.2).

####  3.2.4 Virial temperature of host DM halos $T_{\rm vir}$

We examine three cases of BH accretion in a typical atomic-cooling halo with a
virial temperature of $T_{\rm vir}=10^{4}~{}{\rm K}$ (see Table 1). For all
the cases, the I-front quickly expands and reaches the outer-most radius of
the computational domain. Since the I-front is R-type, the density profile
does not change from the initial distribution except within $1~{}{\rm pc}$,
where the BH gravitational force pulls the ionized gas toward the center and
increases the density as $\rho\propto r^{-3/2}$. However, the accretion rate
of the ionized gas is limited to the Eddington owing to the outward radiation
force via electron scattering. Due to the self-regulated nature of BH
radiative outputs, the gas ejection from the halo takes place in a relatively
slow manner, unlike stellar irradiation and supernova feedback with given
feedback energy (e.g., Kitayama et al., 2004; Kitayama & Yoshida, 2005).
Overall, the stunted growth of seed BHs in “normal” atomic-cooling halos is
consistent with previous studies that use cosmological hydrodynamical
simulations (Habouzit et al., 2017; Latif et al., 2018), although our
simulations focus on the early stage of the bulge formation without including
SN feedback.

###  3.3 Radiative luminosity of rapidly accreting seed BHs

Figure 8: Time-evolution of the mass accretion rate onto a BH (solid), the
Eddington accretion rate (dotted), and the BH-to-bulge mass ratio (dashed) for
the case with a high SFE ($\epsilon_{\star}=0.5$ and $f_{0}=4$). In this case,
both the BH and bulge evolve within the simulation time and the mass ratio
reaches $M_{\bullet}/M_{\star}\simeq 0.04$, which is consistent with those of
bright quasars at $z>6$ (see also §5.2).

Figure 9:  L͡eft: Radiative luminosity produced from a rapidly accreting BH
for the high-SFE case (purple), the fiducial case (green), and the high-
density case with $n_{\rm c}=10^{3.5}~{}{\rm cm}^{-3}$ (cyan). During the
growing phases from $M_{\bullet}=10^{5}~{}M_{\odot}$ to $\sim
10^{6}~{}M_{\odot}$, the luminosities becomes as high as
$L_{\bullet}\simeq(0.3-1.2)\times 10^{45}~{}{\rm erg~{}s}^{-1}$. Right:
Observed radiative flux at $\lambda_{\rm obs}=1.98~{}\mu{\rm m}$ emitted from
an accreting BH in a protogalaxy at $z=15$ predicted for the three models. The
viewing angle is set to $\theta_{\rm obs}=30^{\circ}$. The horizontal, dashed
line shows the $S/N=10$ detection limit of JWST/NIRCam imaging with the F200W
filter in a $10~{}{\rm ks}$ exposure time.

In the left panel of Fig. 9, we present the time evolution of the radiative
luminosity produced by a rapidly accreting seed BH for the high-SFE case
(purple), the fiducial case (green), and the high-density case with $n_{\rm
c}=10^{3.5}~{}{\rm cm}^{-3}$ (cyan). Here, the radiative luminosity is
calculated by using Eq. (23) as a function of the BH feeding rate. We note
that in general, this luminosity is not the one we observe but that injected
through the inner-most cells of the simulation. The two values are almost
identical when and after the accretion transition occurs because intense
radiation ionizes the surrounding gas and further blows the gas away, reducing
the amount of absorbers along the line of sight. However, in the early stage
before the accretion transition, the injected radiation flux is attenuated by
the surrounding neutral gas and dust to some degree. When the mass accretion
increases during the transition, the radiative luminosity also rises and
reaches $L_{\bullet}\simeq(0.3-1.2)\times 10^{45}~{}{\rm erg~{}s}^{-1}$,
corresponding to $L_{\bullet}/L_{\rm Edd}\sim 4-16$. We note that the
luminosity of the accreting BH for each case substantially dominates that of
the host galaxy; namely, the LW luminosity we assume as sources of stellar
feedback is $1.9\times 10^{42}~{}{\rm erg~{}s}^{-1}$ (fiducial case) and
$7.6\times 10^{43}~{}{\rm erg~{}s}^{-1}$ (high-SFE case), respectively.

With our spectral model for the injected radiation (see Eq. 22), $\sim 55\%$
of the total radiation energy is emitted as photons with
$h\nu<h\nu_{0}=10~{}{\rm eV}$ (the rest-frame wavelength is
$\lambda_{0}=0.124~{}\mu{\rm m}$), which are not absorbed significantly by the
intergalactic media. The specific radiation luminosity at $\nu=\nu_{0}$ is
calculated by $L_{\nu_{0}}=L_{0}=L_{\bullet}/(3.3\nu_{0})$. As an example, we
consider that a seed BH is accreting in a protogalaxy at $z\simeq 15$ (the
luminosity distance is $D_{\rm L}\simeq 167~{}{\rm Gpc}$). In this case, the
wavelength of the rest-frame $h\nu_{0}=10$ eV photons is redshifted to
$\lambda_{\rm obs}\simeq 1.98~{}\mu{\rm m}~{}[(1+z)/16]$, and the observed
radiative flux is calculated via

$F_{\nu_{\rm obs}}=\frac{(1+z)}{4\pi D_{\rm L}^{2}}\cdot\frac{L_{\rm
Edd}+\mathcal{F}(\theta_{\rm obs})\Delta L_{\bullet}}{3.3~{}\nu_{0}},$ (34)

where $\nu_{\rm obs}=\nu_{0}/(1+z)[=c/\lambda_{\rm obs}]$ is the observed
frequency, $\Delta L_{\bullet}=L_{\bullet}-L_{\rm Edd}$, and the anisotropic
degree of the flux is characterized by the function of
$\mathcal{F}(\theta_{\rm obs})$ (see Eq. 24). The right panel of Fig. 9
presents the radiative flux at $\lambda_{\rm obs}=1.98~{}\mu{\rm m}$ emitted
from the accreting BH in a protogalaxy at $z=15$, seen from a viewing angle of
$\theta_{\rm obs}=30^{\circ}$. We note that the x-axis is the time at the
source rest-frame (the redshift effect is not considered). After the
transition of mass accretion, the observed flux becomes as high as $\sim
0.01-0.2~{}{\mu\rm Jy}$, corresponding to $m_{\rm AB}\simeq 29-26$ mag. For
comparison, we overlay the $S/N=10$ detection limit of JWST/NIRCam imaging
with the F200W filter at $1.755\leq\lambda_{\rm obs}/\mu{\rm m}\leq 2.226$ in
a $10~{}{\rm ks}$ exposure time ($m_{\rm AB}\simeq 29~{}{\rm mag}$ denoted by
the dashed line; Rieke et al. 2019)444https://jwst-docs.stsci.edu/jwst-near-
infrared-camera/nircam-instrumentation/nircam-filters. For all the cases, the
radiation flux can be detectable (i.e., $m_{\rm AB}<29$ mag) in durations of
$\Delta t\gtrsim 2$ Myr at the source-rest frame. The corresponding duty cycle
is $f_{\rm duty}\sim 7\%$ of the cosmic time duration when the redshift
changes by $\Delta z=\pm 0.5$. Note that the host galaxy is as bright as $\sim
0.02~{}{\mu\rm Jy}$ ($m_{\rm AB}\gtrsim 28~{}{\rm mag}$) for the high-SFE case
and would not be a serious contamination for hunting seed BHs but rather be a
detectable extended source with JWST.

We finally estimate the number of rapidly accreting seed BHs detectable within
a JWST/NIRCam field of view (${\rm FoV}=9.7~{}{\rm arcmin}^{2}=0.0027~{}{\rm
deg}^{2}$) as

$\displaystyle N_{\rm det}^{\rm JWST}$ $\displaystyle=f_{\rm
duty}\mathcal{N}_{\rm h,seed}\Delta V_{\rm c}(z)\frac{\Delta\Omega_{\rm
obs}}{4\pi}$ $\displaystyle\simeq 1.1\left(\frac{f_{\rm
duty}}{0.1}\right)\left(\frac{\Delta\Omega_{\rm obs}}{10~{}{\rm
FoV}}\right)\left(\frac{\mathcal{N}_{\rm h,seed}}{10^{-4}~{}{\rm
cMpc}^{-3}}\right),$ (35)

where $\Delta\Omega_{\rm obs}$ is the observed solid angle and
$\mathcal{N}_{\rm h,seed}$ is the comoving number density of DM halos that
host growing seed BHs. Here, we consider only the BH luminosity, because the
seed BHs outshine their host galaxies. The number density of DM halos with
$M_{\rm h}\geq 2\times 10^{9}~{}M_{\odot}$ at $z=15$ is estimated as
$\mathcal{N}_{\rm h}\simeq 2.5\times 10^{-3}~{}{\rm Mpc}^{-3}$. Therefore, if
$\gtrsim 4\%$ of those DM halos at $z=15$ contain seed BHs, we are able to
detect one rapidly accreting seed BH by JWST in its 10 fields of view.

Deep imaging surveys by the RST will also enable us to hunt a larger number of
rapidly accreting BHs in the early universe. Assuming the survey area is as
large as $\Delta\Omega_{\rm obs}\simeq 40~{}{\rm deg}^{2}$ and the limiting
magnitude is $m_{\rm AB}\simeq 28-29$ with the F184 filter at
$1.68\leq\lambda_{\rm obs}/\mu{\rm m}\leq 2.00$ as planned in Akeson et al.
(2019), the expected number of detected seed BHs would be

$\displaystyle N_{\rm det}^{\rm RST}\simeq 1.6\times 10^{3}\left(\frac{f_{\rm
duty}}{0.1}\right)\left(\frac{\Delta\Omega_{\rm obs}}{40~{}{\rm
deg}^{2}}\right)\left(\frac{\mathcal{N}_{\rm h,seed}}{10^{-4}~{}{\rm
cMpc}^{-3}}\right).$ (36)

We note that if the limiting magnitude is $m_{\rm AB}<27~{}{\rm mag}$, the
observable duration becomes shorter and the number of detections is lowered.

The detailed modeling of their SED and the multi-band selection and
spectroscopic diagnosis of accreting seed BH candidates will be discussed in
future work. Note that some previous studies modeled the radiation spectrum of
an accreting seed BH using a spherically-symmetric one-dimensional flow
structure (Pacucci et al., 2015; Natarajan et al., 2017; Valiante et al.,
2018).

##  4 Analytical derivation of the conditions for rapid accretion

We here derive the conditions required for the onset of rapid mass accretion
onto a seed BH embedded within bulge stars. Based on the RHD simulation
results, the transition is triggered when (1) the I-front is confined within
the core region without propagating outward and (2) a sufficient amount of
neutral gas is supplied from larger radii without being suppressed by
radiative feedback. In what follows, we quantify the two conditions with
analytical expressions.

Let us consider that an accreting BH emits ionizing radiation at a rate of
$Q_{0}$ and the radiation propagates outward within a gas cloud with a
distribution given by Eq. (29). As a reference value, we estimate the size of
the I-front in a uniform density with $n_{\rm c}$ as

$\displaystyle r_{\rm ion}$ $\displaystyle=\left(\frac{3Q_{0}}{4\pi\alpha_{\rm
rec,B}n_{\rm c}^{2}}\right)^{1/3},$ $\displaystyle\simeq
25.5~{}\lambda^{1/3}M_{5}^{1/3}n_{\rm c,3}^{-2/3}\left(\frac{T_{\rm
HII}}{7\times 10^{4}~{}{\rm K}}\right)^{0.28}~{}{\rm pc}$ (37)

where $\alpha_{\rm rec,B}$ is the case-B radiative recombination rate
coefficient, $T_{\rm HII}$ is the temperature of the ionized gas, the number
flux of ionizing photons is given by $Q_{0}\simeq
0.173~{}L_{\bullet}/(h\nu_{\rm Ly})$ for the radiation spectral model (Eq.
22), and $\lambda(\equiv L_{\bullet}/L_{\rm Edd})$ is the Eddington ratio.
However, this estimate in Eq. (37) is no longer valid when the I-front size
reaches the radius where the density distribution decreases steeply. With the
gas remaining at rest, the critical power-law index of the density profile
($\rho\propto r^{-\beta}$) for continuous expansion of the I-front is given by
$\beta=3/2$ (Franco et al., 1990). Therefore, for the given initial density
profile, one defines the critical radius, outside which the power-law index is
steeper than the critical value and the I-front expansion is accelerated, by
$r_{\rm crit}=\sqrt{3}~{}r_{0}$ or

$r_{\rm crit}=27.4~{}T_{\rm vir,5}^{1/2}~{}n_{\rm c,3}^{-1/2}~{}{\rm pc}.$
(38)

Thus, the first condition for the transition is given by $r_{\rm crit}\gtrsim
r_{\rm ion}$ or

$M_{\bullet}\lesssim 2\times 10^{5}~{}T_{\rm vir,5}^{3/2}~{}n_{\rm
c,3}^{1/2}\left(\frac{\lambda}{0.5}\right)^{-1}~{}M_{\odot},$ (39)

where the Eddington ratio is set to $\lambda=0.5$, which is the typical value
before the transition to rapid mass accretion.

Next, we consider the condition where radiative feedback does not affect mass
inflows from larger radii. Namely, this requires the effective gravitational
influence radius to be larger than the size of the ionized region, i.e.,
$r_{\rm inf}\gtrsim r_{\rm ion}$. The effective gravitational influence radius
is calculated with the balance between the total (BH+bulge+DM) gravitational
force and the gas-pressure gradient force as

$\displaystyle\alpha\frac{c_{\rm s}^{2}}{r}\simeq$
$\displaystyle\frac{GM_{\bullet}}{r^{2}}+\frac{GM_{\star}}{(r+r_{\rm
c})^{2}}-\frac{2k_{\rm B}T_{\rm vir}f(c_{\rm N})}{\mu m_{\rm p}r_{\rm s}}g(x)$
(40)

where $x=r/r_{\rm s}$, $f(c_{\rm N})\simeq 5$ at $1\lesssim c_{\rm N}\lesssim
5$, and $g(x)=[x-(1+x)\ln(1+x)]/[x^{2}(1+x)]$. In the cases of interest, where
$r_{\rm ion}\lesssim r_{\rm crit}\ll r_{\rm s}$, the third term in the right-
hand-side can be approximated as $\approx 5xk_{\rm B}T_{\rm vir}/(\mu m_{\rm
p})$ and the first term is negligible compared to the second term. Therefore,
the force balance at $r_{\rm crit}\lesssim r\ll r_{\rm s}$, where
$\alpha\simeq 2$, is governed by the equation of

$2c_{\rm s}^{2}\simeq\frac{GM_{\star}r}{(r+r_{\rm c})^{2}}+\frac{5k_{\rm
B}T_{\rm vir}}{\mu m_{\rm p}}\left(\frac{r}{r_{\rm s}}\right).$ (41)

In addition, when $M_{\star}\ll 1.7\times 10^{8}~{}T_{\rm vir,5}n_{\rm
c,3}^{-1}~{}M_{\odot}$, the bulge size is sufficiently small ($r\gg r_{\rm
c}$), and the gravitational influence radius is expressed by555Another
solution of $r\gtrsim 0.4~{}(T_{\rm g}/T_{\rm vir})r_{\rm s}$ is not adopted..

$r_{\rm inf}\simeq\frac{GM_{\star}}{2c_{\rm s}^{2}}.$ (42)

Note that Eq. (42) is identical to the bulge-gravitational influence radius
that approximates the bulge to be a point gravitational source. Therefore, the
condition of $r_{\rm inf}\gtrsim r_{\rm ion}$ is rewritten as

$\frac{M_{\star}}{M_{\bullet}}\gtrsim 82~{}M_{5}^{-2/3}n_{\rm
c,3}^{-2/3}\left(\frac{\lambda}{0.5}\right)^{1/3}\left(\frac{T_{\rm
HII}}{7\times 10^{4}~{}{\rm K}}\right)^{1.28},$ (43)

where the sound speed is evaluated with the temperature in the ionized region
as $c_{\rm s}\simeq 29.6~{}{\rm km~{}s}^{-1}(T_{\rm HII}/7\times 10^{4}~{}{\rm
K})^{1/2}$. Indeed, for $M_{\bullet}=10^{5}~{}M_{\odot}$, the transient super-
Eddington accretion mode can be triggered when the bulge mass is as massive as
$M_{\star}\gtrsim 10^{7}~{}M_{\odot}$ (our fiducial case). Finally, combining
the two conditions given by Eqs. (39) and (43), we obtain the bulge-BH mass
relation required to feed the BH efficiently as

$\frac{M_{\star}}{M_{\bullet}}\gtrsim 1.03\times 10^{2}~{}\lambda~{}n_{\rm
c,3}^{-1}~{}T_{\rm vir,5}^{-1},$ (44)

where the dependence on $T_{\rm HII}$ is omitted. For our cases with
$M_{\bullet}=10^{5}~{}M_{\odot}$, the critical bulge mass for the onset of
rapid accretion is estimated from Eq. (44) as $M_{\star}\simeq
10^{7}~{}M_{\odot}$. This is consistent with the result of our parameter
studies discussed in §3.2.2.

In summary, for two given quantities to characterize the properties of gas in
a DM halo ($T_{\rm vir}$ and $n_{\rm c}$), we discuss whether seed BHs can
grow via mass accretion at rates exceeding the Eddington value. This transient
growing phase takes place when the ionizing radiation produced by the BH does
not break the gaseous core (see Eq. 39) and a massive bulge with
$M_{\star}\gtrsim 100~{}M_{\bullet}$ attracts a sufficient amount of gas
within its gravitational influence radius (see Eq. 44).

Finally, we note that the critical bulge-to-BH mass ratio in Eq. (44) is $\sim
10$ times higher than that obtained in the previous study by Park et al.
(2016), where the gravitational effect of bulge stars is considered in
spherically symmetric one-dimensional RHD simulations for BH accretion.
Compared to their work, we further take into account the multi-dimensional
effect (e.g., gas angular momentum), irradiation by bulge stars, and
metallicity and also adopt a different bulge model666 Park et al. (2016)
modeled it so that the mean stellar density within $r_{\rm c}$ is consistent
with that of the Milky Way. On the other hand, we model the bulge mass-size
relation so that the observed $M_{\star}-\sigma_{\rm e}$ relation is
reproduced (see §2), yielding a more concentrated bulge.. However, the
difference on the critical mass ratio is mainly caused by the initial density
distribution. Since Park et al. (2016) adopted a uniform density distribution
as their initial conditions, the I-front always has a maximum size without
continuous expansion. On the other hand, in our cases where the initial
density profile consists of the core and envelope that follows $n\propto
r^{-2}$. Therefore, the condition of Eq. (39) has to be required to trigger
rapid BH accretion. If we adopted a uniform density distribution, only the
condition of Eq. (43) would be considered. For $M_{\star}=10^{6}~{}M_{\odot}$,
the critical mass ratio estimated from Eq. (43) is
$M_{\star}/M_{\bullet}\simeq 18$, which is consistent with that in Park et al.
(2016).

##  5 Discussion

###  5.1 Feasible conditions of rapid BH accretion in the galaxy assembly

In this section, we discuss whether seed BHs formed in high-redshift
protogalaxies experience rapidly growing phases. Let us first consider a BH
with $M_{\bullet}=10^{5}~{}M_{\odot}$. From the condition of Eq. (39), the
host DM halo is required to be as massive as

$M_{\rm h}\gtrsim 2\times 10^{9}~{}M_{\bullet,5}~{}n_{\rm
c,3}^{-1/2}\left(\frac{1+z}{16}\right)^{-3/2}~{}M_{\odot},$ (45)

or equivalently, $T_{\rm vir}\gtrsim 10^{5}~{}M_{\bullet,5}^{2/3}n_{\rm
c,3}^{-1/3}~{}{\rm K}$. Now, we define the ratio of the galaxy stellar mass to
the DM halo mass as $f_{\star}(\equiv M_{\star}/M_{\rm h})$ and assume that
the bulge mass in the protogalaxy is comparable to the total stellar mass.
Then, using the condition of Eq. (44), we obtain

$f_{\star}\gtrsim 5.2\times 10^{-3}~{}M_{\bullet,5}^{-2/3}~{}n_{\rm
c,3}^{-1/6}\left(\frac{1+z}{16}\right)^{3/2}.$ (46)

The required value of $f_{\star}$ corresponds to $\gtrsim 3\%$ of the
conversion efficiency from gas into stars (i.e., the star formation efficiency
$\epsilon_{\star}$), assuming the cosmic mean baryon fraction $f_{\rm
b}=\Omega_{\rm b}/\Omega_{\rm m}\simeq 0.16$. This value is consistent with
those inferred by abundance matching and the observed UV luminosity function
of galaxies at $z\simeq 6$ (Bouwens et al., 2015). Moreover, this choice of
$\epsilon_{\star}\sim 0.03-0.05$ explains the cosmic reionization history
without violating the optical depth of the universe to electron scattering
measured by the Planck satellite (Visbal et al., 2015; Inayoshi et al.,
2021)777 An empirical model for linking galaxy star formation rates to the
properties of their host haloes (UniverseMachine; ) predicts the bulge-to-halo
mass ratio as $f_{\star}\sim 3\times 10^{-4}(M_{\rm
h}/10^{9}~{}M_{\odot})^{0.75}$ at $10<z<15$ over $10^{9}\lesssim M_{\rm
h}/~{}M_{\odot}\lesssim 10^{10}$. The error size with the $84\%$ confidence is
$\sim\pm(0.3-0.5)$ dex at $z\simeq 10$ and $\sim\pm(0.7-1.0)$ dex at $z\simeq
14$, respectively. Indeed, the upper envelope of the predicted values at the
halo mass range is consistent with $f_{\star}\sim 5\times 10^{-3}$.. Adopting
a single value of $\epsilon_{\star}=0.05$, therefore the condition of Eq. (46)
is rewritten as

$M_{\bullet}\gtrsim 5.2\times 10^{4}~{}n_{\rm
c,3}^{-1/4}\left(\frac{\epsilon_{\star}}{0.05}\right)^{-3/2}\left(\frac{1+z}{16}\right)^{9/4}~{}M_{\odot},$
(47)

which holds for $M_{\bullet}=10^{5}~{}M_{\odot}$ seed BHs at $z<20$.

Figure 10: Summary of the analytical argument in §5.1 on the redshift-halo
mass plane. In the region enclosed with the two conditions of Eqs. (45) and
(47), rapid mass accretion onto seed BHs with $M_{\bullet}=10^{5}~{}M_{\odot}$
can take place (red curve). The star symbols represent the mass and redshift
of the DM halos studied in this paper; Massive halo (red) and Normal halo
(blue), respectively (see Table 1). The gray curves show $10^{4}$ merger trees
of the high-redshift quasar main progenitors, corresponding to the $3-4\sigma$
mass variance, and the black dashed curve presents the median halo mass. For
comparison, the blue curve is the assembly history of the typical first
galaxies, corresponding to the $2\sigma$ mass variance. The dotted curves
indicate constant virial temperatures, the values of which are denoted by
numbers in the figure.  Figure 11: The $M_{\bullet}-M_{\star}$ relation for
the three simulations: B7T5N3-highSFE (red), B7T5N3 (magenta), and B6T5N3
(purple). Circle symbols show the $z>6$ quasar samples compiled by Izumi et
al. (2019, 2021); brighter ones with $M_{1450}<-25$ (blue) and fainter ones
with $M_{1450}>-25$ (cyan). The cross symbols are the observational samples in
the local universe provided by Kormendy & Ho (2013). Each line represents the
local relation of $M_{\bullet}/M_{\star}$ (black solid, with $\sim 1\sigma$
errors), and the best-fit relations for the brightest $z>6$ quasars by
Pensabene et al. (2020). The red dotted curve presents the evolutionary track
of the $M_{\bullet}/M_{\star}$ ratio after rapid accretion phases predicted by
a phenomenological model (see Appendix B). As a reference, the evolutionary
tracks of the $M_{\bullet}/M_{\star}$ ratio obtained by cosmological
simulations (Zhu et al., 2020) are overlaid (green curves).

To see the feasibility of rapid BH accretion in the hierarchical structure
formation paradigm, we consider two different assembly histories of DM halos.
One is the evolution of the high-redshift quasar main progenitors, i.e., the
most massive halos at each epoch among halos that end up in $M_{\rm
h}=10^{12}~{}M_{\odot}$ at $z=6$. In Fig. 10, we present $10^{4}$ merger trees
of the high-redshift quasar main progenitors (gray curves), corresponding to
the $3-4\sigma$ mass variance, and the median halo mass is shown by the dashed
curve. Along with the merger history of the high-redshift quasar hosts, there
is a parameter region where the two conditions given by Eqs. (45) and (47) are
satisfied and the seed BH can experience transient super-Eddington accretion
phases (red curve). The other one is the assembly history of the typical first
galaxies that form in DM halos with $M_{\rm h}\gtrsim 10^{7}~{}M_{\odot}$ at
$z\simeq 10$, corresponding to the $2\sigma$ mass variance (blue curve). Along
with this evolutionary track, seed BHs with $M_{\bullet}=10^{5}~{}M_{\odot}$
could grow via rapid accretion at $z<6$ (i.e., $T_{\rm vir}\gtrsim
10^{5}~{}{\rm K}$). It is also worth noting that strong suppression of BH
growth from $M_{\bullet}\sim 10^{5}~{}M_{\odot}$ in typical galaxies with
$T_{\rm vir}\simeq 10^{4}~{}{\rm K}$ is commonly seen in cosmological
hydrodynamical simulations. For instance, Latif et al. (2018) showed that the
early growth of a $10^{5}~{}M_{\odot}$ BH is totally quenched in an atomic-
cooling halo with the $2\sigma$ mass variance.

Next, we briefly discuss the cases with different BH masses. When the seed
mass is $M_{\bullet}=10^{4}~{}M_{\odot}$, the halo condition yields $T_{\rm
vir}\gtrsim 2\times 10^{4}~{}{\rm K}$, but the high value of $f_{\star}$ can
be achieved at lower redshift of $1+z\lesssim
7.8~{}(\epsilon_{\star}/0.05)^{2/3}$. On the other hand, when the seed mass is
$M_{\bullet}=10^{6}~{}M_{\odot}$, the halo condition is satisfied in massive
DM halos with $M_{\rm h}\gtrsim 10^{11}~{}M_{\odot}$ (or $T_{\rm vir}\gtrsim
4.6\times 10^{5}~{}{\rm K}$), where all the high-redshift quasar progenitor
halos are expected to have sufficiently massive bulges. However, the host
galaxies formed in those massive halos would be already polluted by heavy
elements. While our simulation results hold for lower metallicity environments
with $Z\simeq 0.01~{}Z_{\odot}$, the radiative feedback effect caused by the
accreting BH would be stronger in more metal-enriched environments. This would
quench their efficient growth of BHs and change the accretion mode to the
Eddington-limited one. To explore the critical metallicity to terminate super-
Eddington accretion and the impact of cosmological metal-enrichment processes
is left for future work.

###  5.2 Early coevolution of seed BHs with host galaxies

The empirical relation between the mass of SMBHs and the properties of their
host galaxies are considered to be one of the most important outcomes caused
by their coevolution over the cosmic timescale (e.g., Magorrian et al., 1998;
Ferrarese & Merritt, 2000; Kormendy & Ho, 2013). Theoretical models for
explaining the tight correlations have been proposed but the origin is still
unclear. To understand the nature of these correlations, it is critically
important to study them beyond the local universe, characterizing how and when
the relations have been established and evolved until now. So far, a large
number of observational studies have extensively investigated the redshift
dependence of the BH-to-bulge mass ratio of $M_{\bullet}/M_{\star}$ and
overall suggested its positive redshift dependence, i.e., the ratio increases
with redshift (Bennert et al., 2011; Schramm & Silverman, 2013; Ding et al.,
2020). Beyond $z\sim 6$, ALMA is a powerful tool to measure the dynamical mass
of gas in quasar host galaxies and allows us to explore the early stage of the
BH/galaxy correlation (e.g., Wang et al., 2010, 2013; Venemans et al., 2017).
In addition, observations with the Subaru HSC provide low-luminosity and less
massive BH samples, which are unique populations to determine the
$M_{\bullet}/M_{\star}$ ratio at $z>6$ (Izumi et al., 2019, 2021). Fig. 11
shows the distribution of $z>6$ quasars compiled in Izumi et al. (2021),
together with those in the local universe (Kormendy & Ho, 2013). First, the
brightest $z>6$ quasars with $M_{1450}<-25$ mag tend to have
$M_{\bullet}/M_{\star}$ ratios higher than those seen in the local universe.
Namely, the mass ratio for those brightest objects is boosted by a factor of
$\sim 10$ (blue dashed line; Pensabene et al. 2020). On the other hand, the
fainter quasars with $M_{1450}>-25$ mag appear to follow the local relation,
although those BHs are considered to grow at rates of $\gtrsim 0.05~{}{\rm
SFR}$ and will be overmassive at lower redshifts. We note that for all the
$z>6$ samples, the values of the x-axis are not the bulge mass of their host
galaxies but the dynamical mass measured by [CII] 158 $\mu{\rm m}$ lines. In
general, the dynamical mass is considered to be higher than the true bulge
mass. With a high-resolution ALMA observation, Izumi et al. (2021) found that
the gas dynamics of the core component of a low-luminosity quasar at $z=7.07$
(HSC J1243+0100) is governed by rotation associated with a compact bulge and
estimated its mass as $\sim 50\%$ of the [CII]-based dynamical mass.
Therefore, the correlation at $z>6$ might be shifted to the left if the
conversion factor from the dynamical mass to the bulge mass is taken into
account.

At the left-bottom conner of Fig. 11, we present the evolution tracks of the
$M_{\bullet}/M_{\star}$ ratio obtained from our simulations. When the bulge
mass is $M_{\star}<10^{7}~{}M_{\odot}$ (purple curve), the seed BH hardly
grows in mass owing to radiative feedback and the $M_{\bullet}/M_{\star}$
moves to the right, i.e., only the bulge mass grows. When the bulge is as
massive as $M_{\star}\geq 10^{7}~{}M_{\odot}$ (magenta and red curves), the
efficient mass accretion phases are triggered and thus the BH mass abruptly
increases within $\sim$ a few Myrs. As a result, the $M_{\bullet}/M_{\star}$
ratio becomes as high as $\sim 0.03-0.05$, which is consistent with those of
$z>6$ (see also the line obtained by Pensabene et al. 2020). In the high-SEF
model ($\epsilon_{\star}=0.5$), where the bulge mass reaches $M_{\star}\simeq
3\times 10^{7}~{}M_{\odot}$ within $\sim 3$ Myr, the BH-to-bulge mass ratio is
$M_{\bullet}/M_{\star}\simeq 0.03$ at the end of the simulation. After the
rapid BH growing phase terminates, both the BH and host galaxy evolve in mass
and reach the observed values as high-redshift quasars. As an example, we
present the subsequent BH/galaxy evolutionary track to a $z>6$ HSC quasar
(Izumi et al., 2019) predicted by a phenomenological model (see more details
in Appendix B). In summary, seed BHs formed in high-redshift quasar hosts can
be substantially overmassive during the early bulge assembly at $z\gtrsim 15$
and the mass ratio would approach the ratio of $M_{\bullet}/M_{\star}$ seen at
$z\sim 6$ via the subsequent growth.

It is worth noting that most cosmological simulations studying the BH-galaxy
coevolution concluded that seed BHs hardly grow in mass via gas accretion
because dense, cold gas is expelled by energetic SN feedback associated with
star formation (e.g., Sijacki et al., 2009; Habouzit et al., 2017; Lupi et
al., 2019). As a result of SN feedback, BH growth is strongly quenched until
the host galaxies become heavier than $M_{\star}\gtrsim 10^{9}~{}M_{\odot}$,
whose critical stellar mass depends on the sub-grid models for star formation,
SN feedback, and AGN feedback. There is a model parameter set for which the
most massive BHs in their simulations can reach $M_{\bullet}\simeq
10^{9}~{}M_{\odot}$ by $z\sim 6$. However, the predicted shape of the
$M_{\bullet}-M_{\star}$ relation shows that most of the BHs are undermassive,
compared to the local relation (e.g., Zhu et al., 2020; Valentini et al.,
2021). As a reference, the evolutionary tracks of the $M_{\bullet}/M_{\star}$
ratio obtained by Zhu et al. (2020)888 Zhu et al. (2020) have extensively
investigated the effect of feedback, BH seeding, and accretion models on the
BH growth. Among their simulation results, we show the cases where the initial
BH mass is $M_{\bullet}=10^{5}~{}M_{\odot}$ for comparison. are overlaid in
Fig. 11 (green curves).

As discussed in Inayoshi et al. (2020), there are important limitations of
current numerical simulations. First, most large-scale cosmological
simulations resolve the dynamics of DM/gas/stars on galactic scales at $\sim
O({\rm kpc})$, but do not resolve the BH gravitational influence radius.
Secondly, owing to simplified star formation models where gas particles denser
than a threshold are replaced with stars, dense clouds in the nuclear region
would be disrupted and thus the BH growth could be quenched as seen in many
cosmological simulation studies. In contrast, as shown in our simulations that
resolve sub-pc scales, a fraction of seed BHs that were born in highly biased
regions of the universe with mass variance of $\gtrsim 3-4~{}\sigma$ could be
fed through dense, cold accretion flows (see also Li et al., 2021b).
Additionally, the existence of such overmassive BHs in protogalaxies will
provide us a unique opportunity of detecting highly accreting seed BHs in the
very early universe at $z>10$, unlike the undermassive-BH scenarios supported
by cosmological simulations.

###  5.3 Young quasars with low radiative efficiencies

Our RHD simulations suggest the existence of high-redshift quasars that
accrete at super-Eddington rates and fade out in a timescale of $\sim{\rm
a~{}few}$ Myr. The duration of such rapid accretion is generally consistent
with lifetimes of $z\gtrsim 6$ quasars ($t_{\rm Q}\lesssim 1-10~{}{\rm Myr}$),
which are estimated by the measurement of the physical extents of hydrogen
Ly$\alpha$ proximity zones observed in the rest-frame UV spectra (e.g., Eilers
et al., 2018; Davies et al., 2019; Eilers et al., 2021). Since the inferred
quasar lifetimes are substantially shorter than the $e$-folding timescale
assuming Eddington-limited accretion, $t_{\rm Edd}~{}(\equiv
M_{\bullet}/\dot{M}_{\rm Edd})\simeq 45~{}{\rm Myr}$, some of those $z\gtrsim
6$ quasars are expected to undergo radiatively inefficient super-Eddington
accretion to grow up to $M_{\bullet}\sim 10^{9}~{}M_{\odot}$ in such a short
duration. For instance, Davies et al. (2019) proposed that two $z>7$ quasars
(ULAS J1120+0641 and ULAS J1342+0928) would have a small radiative efficiency
significantly below $\eta_{0}\simeq 0.1$. Based on the above, mupltiple
episodes of rapidly accretion would be an intriguing counterpart of such young
high-redshift quasars, although the observed quasars are hosted in
substantially heavier halos and exhibit star formation rates higher than we
consider. Moreover, the episodic nature of quasar activity produces a complex
structure of ionized and neutral gas surrounding the growing BH, unlike the
situation where continuous output of radiation from the BH is assumed (the so-
called “light-bulb” light curve model). The variable light-curve effect on the
formation of proximity zones needs to be explored in more detail (see e.g.,
Davies et al., 2020). To perform RT calculations of Ly$\alpha$ lines under a
realistic density profile obtained from simulations for BH accretion is left
for future investigations.

##  6 Summary

We study the early growth of massive seed BHs via accretion at the centers of
protogalaxies where the stellar bulge component is assembled, performing
axisymmetric two-dimensional RHD simulations. We find that when a seed BH with
$M_{\bullet}\sim 10^{5}~{}M_{\odot}$ is embedded in dense metal-poor gas
($Z=0.01~{}Z_{\odot}$) with a density of $n_{\rm c}\gtrsim 100~{}{\rm
cm}^{-3}$ and bulge stars with a total mass of $M_{\star}\gtrsim
100~{}M_{\bullet}$, a massive gaseous disk feeds the BH efficiently at super-
Eddington rates of $\gtrsim 0.3-1~{}M_{\odot}~{}{\rm yr}^{-1}$. The peak rate
of $\dot{M}_{\bullet}\sim 10~{}M_{\odot}~{}{\rm yr}^{-1}$ corresponds to $\sim
2.5\times 10^{3}~{}\dot{M}_{\rm Edd}$. This rapid accretion phase lasts until
a good fraction of the gas bounded within the bulge accretes onto the BH,
although the feeding rate is regulated owing to strong outflows driven by
ionizing radiation emitted from the accreting BH. As a result, the BH mass
increases to $\sim 10^{6}~{}M_{\odot}$ within $\sim 2$ Myr after the onset of
accretion bursts. In contrast, when the seed BH is surrounded by diffuse
gaseous media and/or less massive bulges, the BH accretion proceeds
episodically due to radiative feedback and thus the average accretion rate is
limited below the Eddington rate.

We provide the analytical expressions of the conditions required for the onset
of super-Eddington accreting phases of seed BHs embedded within growing bulge
stars. This transient growing phase takes place when (1) the ionizing
radiation produced by the BH does not break the gaseous core (see Eq. 39) and
(2) a massive bulge with $M_{\star}\gtrsim 100~{}M_{\bullet}$ attracts a
sufficient amount of gas within its gravitational influence radius (see Eq.
44). In the framework of the cosmological hierarchical structure formation,
those conditions can be achieved for seed BHs formed in massive DM halos with
masses of $\gtrsim 10^{9}~{}M_{\odot}$ at $z\sim 15-20$ (the virial
temperature is $T_{\rm vir}\simeq 10^{5}~{}{\rm K}$). The host halos are
heavier and rarer than those of typical first galaxies with $T_{\rm vir}\simeq
10^{4}~{}{\rm K}$, but are more likely to end up in quasar hosts by $z\simeq
6$.

The rapid accretion mechanism found with our simulations that resolve sub-pc
scales naturally yields a high BH-to-bulge mass ratio of
$M_{\bullet}/M_{\star}\simeq 0.03-0.05$. The ratio is significantly higher
than that seen in the local universe, but is consistent with those of the
brightest high-redshift quasars with $M_{1450}<-25$ mag. The rarity of those
overmassive BHs in the brightest quasars can also be explained by the fact
that the transient rapid growth would take place in highly biased regions of
the universe with mass variance of $3-4\sigma$. In contrast, most large-scale
cosmological simulations that resolve the dynamics of DM/gas/stars on galactic
scales at $\sim O({\rm kpc})$ predict that SMBHs (or their seeds) tend to be
undermassive at high redshifts compared to the local BH-galaxy correlation.
Both observational and numerical studies are further required to better
understand the early development of the BH-galaxy correlation.

Moreover, the existence of such overmassive BHs provides us a unique
opportunity of detecting highly accreting seed BHs by the upcoming
observations by JWST. In fact, seed BHs accreting at super-Eddington rates
produce radiative luminosities of $L_{\bullet}\simeq(0.3-1.2)\times
10^{45}~{}{\rm erg~{}s}^{-1}$, corresponding to $L_{\bullet}/L_{\rm Edd}\sim
4-16$. If those BHs are in $z\simeq 15$ protogalaxies, the radiation flux at
$\lambda_{\rm obs}\simeq 2~{}\mu{\rm m}$ ($h\nu=10$ eV in the rest frame) is
estimated as $\sim 0.01-0.2~{}{\rm\mu Jy}$, which corresponds to $m_{\rm
AB}=26-29$ mag and can be detectable with the $S/N=10$ detection limit of JWST
NIRCam imaging with the F200W filter in a 10 ks exposure time. If $\gtrsim
4\%$ of the DM halos where the transition conditions are satisfied contain
seed BHs, we are able to detect one rapidly accreting seed BH by JWST in its
10 fields of view.

We greatly thank Luis Ho, Takuma Izumi, and Mingyang Zhuang for constructive
discussions. K. I. acknowledges support from the National Natural Science
Foundation of China (12073003, 12003003, 11721303, 11991052, 11950410493), the
National Key R&D Program of China (2016YFA0400702). R. N. is supported by the
Special Postdoctoral Researcher (SPDR) Program at RIKEN and by a Grant-in-Aid
for Research Activity Start-up (19K23469). T. H. acknowledges support from the
JSPS KAKENHI grant Nos. 17H06360, 17H01102, and 19H01934. R. K. acknowledges
financial support via the Emmy Noether and Heisenberg Research Grants funded
by the German Research Foundation (DFG) under grant no. KU 2849/3 and 2849/9.
R. K. also acknowledges financial support via the JSPS Invitational Fellowship
for Research in Japan under the Fellowship ID S20156. The numerical
simulations were performed with the Cray XC50 at the Center for Computational
Astrophysics (CfCA) of the National Astronomical Observatory of Japan and with
the High-performance Computing Platform of Peking University.

##  Appendix A bulge size-mass relation

We consider the gravitational effect of bulge stars on the thermal dynamics of
gas. In our simulations, we adopt a Hernquist density profile (Hernquist,
1990)

$\rho(r)=\frac{M_{\star}}{2\pi}\frac{r_{\rm c}}{r(r+r_{\rm c})^{3}},$ (A1)

where $M_{\star}$ is the total bulge mass and $r_{\rm c}$ is the core radius.
Integrating the density profile, the enclosed mass is given by

$M(r)=M_{\star}\frac{r^{2}}{(r+r_{\rm c})^{2}},$ (A2)

and thus the half-mass radius is calculated by

$R_{\rm e}=(1+\sqrt{2})r_{\rm c}.$ (A3)

Assuming isotropic motion of stars, the velocity dispersion is analytically
calculated as

$\sigma^{2}(r)=\frac{GM_{\star}}{r_{\rm
c}}\left[xy^{3}\ln\left(\frac{y}{x}\right)-\frac{x}{y}\left(\frac{1}{4}+\frac{y}{3}+\frac{y^{2}}{2}+y^{3}\right)\right],$
(A4)

where $x=r/r_{\rm c}$ and $y=1+x$. The velocity dispersion at the half-mass
radius is

$\sigma_{e}^{2}\simeq 0.055~{}\frac{GM_{\star}}{r_{\rm c}}.$ (A5)

In our paper, we quantify the core radius of $r_{\rm c}$ so that the
correlation between $M_{\star}$ and $\sigma$ becomes consistent with that
observed in the local universe (Kormendy & Ho, 2013), as shown in Eq. (9).

##  Appendix B Phenomenological model for the $M_{\bullet}/M_{\star}$ ratio

We provide a phenomenological model for calculating the redshift-dependent
$M_{\bullet}/M_{\star}$ ratio. First, we characterize the mass growth of a BH
using an exponential function of time $t$,

$\frac{M_{\bullet}}{M_{\bullet i}}=\exp\left[\langle
A\rangle~{}\frac{t-t_{i}}{t_{\rm Edd}}\right],$ (B1)

where $t_{\rm Edd}=45$ Myr is the Salpeter timescale, $M_{\bullet i}$ is the
BH mass at $t=t_{i}$, $A$ is the product of the Eddington ratio and duty cycle
of the BH active phase, and $\langle\cdot\rangle$ means the time-averaged
value over the cosmic time duration of interest. Next, we express the mass
growth of the host galaxy with a functional form of $M_{\star}\propto
e^{-Bz}$, which is well-known to nicely reproduce the mass growth of DM halos
$M_{\rm h}$ (e.g., Wechsler et al., 2002; Neistein & Dekel, 2008; Fakhouri et
al., 2010). In fact, this redshift dependence leads to the halo-mass growth
rate of $d(\ln M_{\rm h})/dt\propto(1+z)^{5/2}$, which has been understood
based on the extended Press-Schechter formalism and also derived by a fit to
merger trees from cosmological N-body simulations (Dekel et al., 2013).
Despite the complex nature of galaxy formation processes, cosmological
hydrodynamical simulations have shown that the galaxy mass assembly history
appears to follow a similar function form (e.g., see Fig. 7 in Zhu et al.,
2020). Note that the value of $B$ for galaxy assembly is different form that
for halo assembly. Motivated by this fact, we express the mass growth of the
host galaxy as

$\frac{M_{\star}}{M_{\star i}}=\exp\left[\langle
B\rangle\left(1-\frac{1+z}{1+z_{i}}\right)\right],$ (B2)

where $M_{\star i}$ is the galaxy mass at $z=z_{i}$ and $\langle B\rangle$ is
a parameter to characterize the growth speed of the galaxy. Combining the two
equations, we calculate the mass ratio of $M_{\bullet}/M_{\star}$ as

$\displaystyle\frac{M_{\bullet}}{M_{\star}}$ $\displaystyle=\frac{M_{\bullet
i}}{M_{\star i}}~{}\exp\left[\langle A\rangle~{}\frac{t-t_{i}}{t_{\rm
Edd}}-\langle B\rangle\left(1-\frac{1+z}{1+z_{i}}\right)\right],$
$\displaystyle\simeq\frac{M_{\bullet i}}{M_{\star i}}~{}\exp\left[\langle
A\rangle(D^{3/2}-1)\frac{t_{i}}{t_{\rm Edd}}-\langle
B\rangle(1-D^{-1})\right],$ (B3)

where $D(z)\equiv(1+z_{i})/(1+z)$ and $t_{\rm H}(z_{i})=(2/3H_{0})[\Omega_{\rm
m}(1+z_{i})^{3}]^{-1/2}$ is the Hubble time at $z=z_{i}$, which is a good
approximation in the matter dominant universe ($z>1$). Finally, we determine
$\langle A\rangle$ and $\langle B\rangle$ by

$\langle A\rangle=\frac{t_{\rm Edd}}{t_{i}}\frac{\ln(M_{\bullet f}/M_{\bullet
i})}{D_{f}^{3/2}-1}~{}~{}~{}~{}~{}{\rm and}~{}~{}~{}~{}~{}\langle
B\rangle=\frac{\ln(M_{\star f}/M_{\star i})}{1-D_{f}^{-1}},$ (B4)

where $M_{\bullet f}$ and $M_{\star f}$ are the BH and galaxy mass at
$z=z_{f}$ and $D_{f}=D(z_{f})$.

As an example, we adopt a high-redshift HSC quasar (J1208-0200). For this
quasar, the BH mass, dynamical mass, and redshift are measured as
$M_{\bullet}=7.1\times 10^{8}~{}M_{\odot}$, $M_{\rm dyn}=1.3\times
10^{10}~{}M_{\odot}$, and $z=6.2$, respectively (Izumi et al., 2019). For
simplicity, the galaxy mass is assumed to be equal to the dynamical mass
(i.e., $M_{\star}=M_{\rm dyn}$). Let us consider the BH after the transient
rapid accretion phases ($M_{\bullet i}=1.3\times 10^{6}~{}M_{\odot}$,
$M_{\star i}=3\times 10^{7}~{}M_{\odot}$, and $z_{i}=15$) to be its seed (see
the red curve in Fig. 11). Following the method above, one obtain $\langle
A\rangle=0.4765$ and $\langle B\rangle=11.039$. The evolution track of the
$M_{\bullet}/M_{\star}$ ratio is shown in Fig. 11 (red dashed curve). The mass
ratio decreases initially, reaches a minimum value of $\sim 0.01$ at $z\simeq
10$, and increases to the observed value.

## References

  * Abel et al. (2002) Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93, doi: 10.1126/science.295.5552.93
  * Abramowicz et al. (1988) Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, ApJ, 332, 646, doi: 10.1086/166683
  * Akeson et al. (2019) Akeson, R., Armus, L., Bachelet, E., et al. 2019, arXiv e-prints, arXiv:1902.05569.  https://arxiv.org/abs/1902.05569
  * Anglés-Alcázar et al. (2017) Anglés-Alcázar, D., Faucher-Giguère, C.-A., Quataert, E., et al. 2017, MNRAS, 472, L109, doi: 10.1093/mnrasl/slx161
  * Bañados et al. (2018) Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473, doi: 10.1038/nature25180
  * Bai (2011) Bai, X.-N. 2011, ApJ, 739, 50, doi: 10.1088/0004-637X/739/1/50
  * Balbus & Hawley (1998) Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1, doi: 10.1103/RevModPhys.70.1
  * Becerra et al. (2015) Becerra, F., Greif, T. H., Springel, V., & Hernquist, L. E. 2015, MNRAS, 446, 2380, doi: 10.1093/mnras/stu2284
  * Bennert et al. (2011) Bennert, V. N., Auger, M. W., Treu, T., Woo, J.-H., & Malkan, M. A. 2011, ApJ, 742, 107, doi: 10.1088/0004-637X/742/2/107
  * Bouwens et al. (2015) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34, doi: 10.1088/0004-637X/803/1/34
  * Bromm & Loeb (2003) Bromm, V., & Loeb, A. 2003, ApJ, 596, 34, doi: 10.1086/377529
  * Bromm & Yoshida (2011) Bromm, V., & Yoshida, N. 2011, ARA&A, 49, 373, doi: 10.1146/annurev-astro-081710-102608
  * Bullock et al. (2001) Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559, doi: 10.1046/j.1365-8711.2001.04068.x
  * Chon et al. (2018) Chon, S., Hosokawa, T., & Yoshida, N. 2018, MNRAS, 475, 4104, doi: 10.1093/mnras/sty086
  * Chon & Omukai (2020) Chon, S., & Omukai, K. 2020, MNRAS, 494, 2851, doi: 10.1093/mnras/staa863
  * Ciotti & Ostriker (2001) Ciotti, L., & Ostriker, J. P. 2001, ApJ, 551, 131, doi: 10.1086/320053
  * Davies et al. (2019) Davies, F. B., Hennawi, J. F., & Eilers, A.-C. 2019, ApJ, 884, L19, doi: 10.3847/2041-8213/ab42e3
  * Davies et al. (2020) —. 2020, MNRAS, 493, 1330, doi: 10.1093/mnras/stz3303
  * Dekel & Birnboim (2006) Dekel, A., & Birnboim, Y. 2006, MNRAS, 368, 2, doi: 10.1111/j.1365-2966.2006.10145.x
  * Dekel et al. (2013) Dekel, A., Zolotov, A., Tweed, D., et al. 2013, MNRAS, 435, 999, doi: 10.1093/mnras/stt1338
  * Devecchi & Volonteri (2009) Devecchi, B., & Volonteri, M. 2009, ApJ, 694, 302, doi: 10.1088/0004-637X/694/1/302
  * Di Matteo et al. (2012) Di Matteo, T., Khandai, N., DeGraf, C., et al. 2012, ApJ, 745, L29, doi: 10.1088/2041-8205/745/2/L29
  * Ding et al. (2020) Ding, X., Silverman, J., Treu, T., et al. 2020, ApJ, 888, 37, doi: 10.3847/1538-4357/ab5b90
  * Draine & Lee (1984) Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89, doi: 10.1086/162480
  * Dubois et al. (2013) Dubois, Y., Pichon, C., Devriendt, J., et al. 2013, MNRAS, 428, 2885, doi: 10.1093/mnras/sts224
  * Eilers et al. (2018) Eilers, A.-C., Hennawi, J. F., & Davies, F. B. 2018, ApJ, 867, 30, doi: 10.3847/1538-4357/aae081
  * Eilers et al. (2021) Eilers, A.-C., Hennawi, J. F., Davies, F. B., & Simcoe, R. A. 2021, ApJ, 917, 38, doi: 10.3847/1538-4357/ac0a76
  * Fakhouri et al. (2010) Fakhouri, O., Ma, C.-P., & Boylan-Kolchin, M. 2010, MNRAS, 406, 2267, doi: 10.1111/j.1365-2966.2010.16859.x
  * Fan (2006) Fan, X. 2006, New A Rev., 50, 665, doi: 10.1016/j.newar.2006.06.077
  * Fernández & Metzger (2013) Fernández, R., & Metzger, B. D. 2013, ApJ, 763, 108, doi: 10.1088/0004-637X/763/2/108
  * Ferrarese & Merritt (2000) Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9, doi: 10.1086/312838
  * Franco et al. (1990) Franco, J., Tenorio-Tagle, G., & Bodenheimer, P. 1990, ApJ, 349, 126, doi: 10.1086/168300
  * Fukushima & Yajima (2021) Fukushima, H., & Yajima, H. 2021, MNRAS, 506, 5512, doi: 10.1093/mnras/stab2099
  * Fukushima et al. (2020) Fukushima, H., Yajima, H., Sugimura, K., et al. 2020, MNRAS, 497, 3830, doi: 10.1093/mnras/staa2062
  * Galli & Palla (1998) Galli, D., & Palla, F. 1998, A&A, 335, 403
  * Glover & Jappsen (2007) Glover, S. C. O., & Jappsen, A.-K. 2007, ApJ, 666, 1, doi: 10.1086/519445
  * Grimm et al. (2003) Grimm, H.-J., Gilfanov, M., & Sunyaev, R. 2003, MNRAS, 339, 793, doi: 10.1046/j.1365-8711.2003.06224.x
  * Habouzit et al. (2017) Habouzit, M., Volonteri, M., & Dubois, Y. 2017, MNRAS, 468, 3935, doi: 10.1093/mnras/stx666
  * Haiman (2013) Haiman, Z. 2013, in Astrophysics and Space Science Library, Vol. 396, Astrophysics and Space Science Library, ed. T. Wiklind, B. Mobasher, & V. Bromm, 293, doi: 10.1007/978-3-642-32362-1_6
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359, doi: 10.1086/168845
  * Hirano et al. (2014) Hirano, S., Hosokawa, T., Yoshida, N., et al. 2014, ApJ, 781, 60, doi: 10.1088/0004-637X/781/2/60
  * Hollenbach & McKee (1989) Hollenbach, D., & McKee, C. F. 1989, ApJ, 342, 306, doi: 10.1086/167595
  * Hosokawa et al. (2011) Hosokawa, T., Omukai, K., Yoshida, N., & Yorke, H. W. 2011, Science, 334, 1250, doi: 10.1126/science.1207433
  * Iliev & Shapiro (2001) Iliev, I. T., & Shapiro, P. R. 2001, MNRAS, 325, 468, doi: 10.1046/j.1365-8711.2001.04422.x
  * Inayoshi et al. (2016) Inayoshi, K., Haiman, Z., & Ostriker, J. P. 2016, MNRAS, 459, 3738, doi: 10.1093/mnras/stw836
  * Inayoshi et al. (2019) Inayoshi, K., Ichikawa, K., Ostriker, J. P., & Kuiper, R. 2019, MNRAS, 486, 5377, doi: 10.1093/mnras/stz1189
  * Inayoshi et al. (2021) Inayoshi, K., Kashiyama, K., Visbal, E., & Haiman, Z. 2021, arXiv e-prints, arXiv:2103.12755.  https://arxiv.org/abs/2103.12755
  * Inayoshi & Omukai (2011) Inayoshi, K., & Omukai, K. 2011, MNRAS, 416, 2748, doi: 10.1111/j.1365-2966.2011.19229.x
  * Inayoshi et al. (2014) Inayoshi, K., Omukai, K., & Tasker, E. 2014, MNRAS, 445, L109, doi: 10.1093/mnrasl/slu151
  * Inayoshi & Tanaka (2015) Inayoshi, K., & Tanaka, T. L. 2015, MNRAS, 450, 4350, doi: 10.1093/mnras/stv871
  * Inayoshi et al. (2020) Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27, doi: 10.1146/annurev-astro-120419-014455
  * Inoue (2011) Inoue, A. K. 2011, MNRAS, 415, 2920, doi: 10.1111/j.1365-2966.2011.18906.x
  * Izumi et al. (2019) Izumi, T., Onoue, M., Matsuoka, Y., et al. 2019, PASJ, 71, 111, doi: 10.1093/pasj/psz096
  * Izumi et al. (2021) Izumi, T., Matsuoka, Y., Fujimoto, S., et al. 2021, ApJ, 914, 36, doi: 10.3847/1538-4357/abf6dc
  * Jeon et al. (2012) Jeon, M., Pawlik, A. H., Greif, T. H., et al. 2012, ApJ, 754, 34, doi: 10.1088/0004-637X/754/1/34
  * Jiang et al. (2014) Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2014, ApJ, 796, 106, doi: 10.1088/0004-637X/796/2/106
  * Johnson et al. (2011) Johnson, J. L., Khochfar, S., Greif, T. H., & Durier, F. 2011, MNRAS, 410, 919, doi: 10.1111/j.1365-2966.2010.17491.x
  * Kawinwanichakij et al. (2021) Kawinwanichakij, L., Silverman, J. D., Ding, X., et al. 2021, arXiv e-prints, arXiv:2109.09766.  https://arxiv.org/abs/2109.09766
  * Kitayama & Yoshida (2005) Kitayama, T., & Yoshida, N. 2005, ApJ, 630, 675, doi: 10.1086/432114
  * Kitayama et al. (2004) Kitayama, T., Yoshida, N., Susa, H., & Umemura, M. 2004, ApJ, 613, 631, doi: 10.1086/423313
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Kratter & Lodato (2016) Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271, doi: 10.1146/annurev-astro-081915-023307
  * Kuiper et al. (2010) Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010, ApJ, 722, 1556, doi: 10.1088/0004-637X/722/2/1556
  * Kuiper et al. (2011) —. 2011, ApJ, 732, 20, doi: 10.1088/0004-637X/732/1/20
  * Kuiper et al. (2020) Kuiper, R., Yorke, H. W., & Mignone, A. 2020, ApJS, 250, 13, doi: 10.3847/1538-4365/ab9a36
  * Latif et al. (2013) Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013, MNRAS, 433, 1607, doi: 10.1093/mnras/stt834
  * Latif et al. (2018) Latif, M. A., Volonteri, M., & Wise, J. H. 2018, MNRAS, 476, 5016, doi: 10.1093/mnras/sty622
  * Levermore & Pomraning (1981) Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321, doi: 10.1086/159157
  * Li et al. (2021a) Li, J., Silverman, J. D., Ding, X., et al. 2021a, ApJ, 918, 22, doi: 10.3847/1538-4357/ac06a8
  * Li et al. (2021b) Li, W., Inayoshi, K., & Qiu, Y. 2021b, arXiv e-prints, arXiv:2105.12637.  https://arxiv.org/abs/2105.12637
  * Li et al. (2007) Li, Y., Hernquist, L., Robertson, B., et al. 2007, ApJ, 665, 187, doi: 10.1086/519297
  * Lodato & Natarajan (2006) Lodato, G., & Natarajan, P. 2006, MNRAS, 371, 1813, doi: 10.1111/j.1365-2966.2006.10801.x
  * Lupi et al. (2021) Lupi, A., Haiman, Z., & Volonteri, M. 2021, MNRAS, 503, 5046, doi: 10.1093/mnras/stab692
  * Lupi et al. (2019) Lupi, A., Volonteri, M., Decarli, R., et al. 2019, MNRAS, 488, 4004, doi: 10.1093/mnras/stz1959
  * Lusso et al. (2015) Lusso, E., Worseck, G., Hennawi, J. F., et al. 2015, MNRAS, 449, 4204, doi: 10.1093/mnras/stv516
  * Magorrian et al. (1998) Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285, doi: 10.1086/300353
  * Matsuoka et al. (2016) Matsuoka, Y., Onoue, M., Kashikawa, N., et al. 2016, ApJ, 828, 26, doi: 10.3847/0004-637X/828/1/26
  * Matsuoka et al. (2018) Matsuoka, Y., Strauss, M. A., Kashikawa, N., et al. 2018, ApJ, 869, 150, doi: 10.3847/1538-4357/aaee7a
  * McKinney et al. (2015) McKinney, J. C., Dai, L., & Avara, M. J. 2015, MNRAS, 454, L6, doi: 10.1093/mnrasl/slv115
  * McKinney & Gammie (2004) McKinney, J. C., & Gammie, C. F. 2004, ApJ, 611, 977, doi: 10.1086/422244
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228, doi: 10.1086/513316
  * Milosavljević et al. (2009) Milosavljević, M., Bromm, V., Couch, S. M., & Oh, S. P. 2009, ApJ, 698, 766, doi: 10.1088/0004-637X/698/1/766
  * Mineo et al. (2014) Mineo, S., Gilfanov, M., Lehmer, B. D., Morrison, G. E., & Sunyaev, R. 2014, MNRAS, 437, 1698, doi: 10.1093/mnras/stt1999
  * Mo et al. (1998) Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319, doi: 10.1046/j.1365-8711.1998.01227.x
  * Mortlock et al. (2011) Mortlock, D. J., Warren, S. J., Venemans, B. P., et al. 2011, Nature, 474, 616, doi: 10.1038/nature10159
  * Mowla et al. (2019) Mowla, L., van der Wel, A., van Dokkum, P., & Miller, T. B. 2019, ApJ, 872, L13, doi: 10.3847/2041-8213/ab0379
  * Nakatani et al. (2018a) Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., & Kuiper, R. 2018a, ApJ, 857, 57, doi: 10.3847/1538-4357/aab70b
  * Nakatani et al. (2018b) —. 2018b, ApJ, 865, 75, doi: 10.3847/1538-4357/aad9fd
  * Narayan et al. (2012) Narayan, R., Sa̧dowski, A., Penna, R. F., & Kulkarni, A. K. 2012, MNRAS, 426, 3241, doi: 10.1111/j.1365-2966.2012.22002.x
  * Natarajan et al. (2017) Natarajan, P., Pacucci, F., Ferrara, A., et al. 2017, ApJ, 838, 117, doi: 10.3847/1538-4357/aa6330
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Neistein & Dekel (2008) Neistein, E., & Dekel, A. 2008, MNRAS, 388, 1792, doi: 10.1111/j.1365-2966.2008.13525.x
  * Ohsuga et al. (2005) Ohsuga, K., Mori, M., Nakamoto, T., & Mineshige, S. 2005, ApJ, 628, 368, doi: 10.1086/430728
  * Omukai (2001) Omukai, K. 2001, ApJ, 546, 635, doi: 10.1086/318296
  * Omukai & Nishi (1998) Omukai, K., & Nishi, R. 1998, ApJ, 508, 141, doi: 10.1086/306395
  * Onoue et al. (2019) Onoue, M., Kashikawa, N., Matsuoka, Y., et al. 2019, ApJ, 880, 77, doi: 10.3847/1538-4357/ab29e9
  * Pacucci et al. (2015) Pacucci, F., Ferrara, A., Volonteri, M., & Dubus, G. 2015, MNRAS, 454, 3771, doi: 10.1093/mnras/stv2196
  * Park & Ricotti (2011) Park, K., & Ricotti, M. 2011, ApJ, 739, 2, doi: 10.1088/0004-637X/739/1/2
  * Park & Ricotti (2012) —. 2012, ApJ, 747, 9, doi: 10.1088/0004-637X/747/1/9
  * Park et al. (2016) Park, K., Ricotti, M., Natarajan, P., Bogdanović, T., & Wise, J. H. 2016, ApJ, 818, 184, doi: 10.3847/0004-637X/818/2/184
  * Park et al. (2020) Park, K., Wise, J. H., Bogdanović, T., & Ricotti, M. 2020, ApJ, 905, 92, doi: 10.3847/1538-4357/abc336
  * Pensabene et al. (2020) Pensabene, A., Carniani, S., Perna, M., et al. 2020, A&A, 637, A84, doi: 10.1051/0004-6361/201936634
  * Prieto & Escala (2016) Prieto, J., & Escala, A. 2016, MNRAS, 460, 4018, doi: 10.1093/mnras/stw1285
  * Regan et al. (2014) Regan, J. A., Johansson, P. H., & Haehnelt, M. G. 2014, MNRAS, 439, 1160, doi: 10.1093/mnras/stu068
  * Regan et al. (2020) Regan, J. A., Wise, J. H., Woods, T. E., et al. 2020, The Open Journal of Astrophysics, 3, 15, doi: 10.21105/astro.2008.08090
  * Rieke et al. (2019) Rieke, M., Arribas, S., Bunker, A., et al. 2019, BAAS, 51, 45
  * Sakurai et al. (2016) Sakurai, Y., Inayoshi, K., & Haiman, Z. 2016, MNRAS, 461, 4496, doi: 10.1093/mnras/stw1652
  * Sakurai et al. (2017) Sakurai, Y., Yoshida, N., Fujii, M. S., & Hirano, S. 2017, MNRAS, 472, 1677, doi: 10.1093/mnras/stx2044
  * Sassano et al. (2021) Sassano, F., Schneider, R., Valiante, R., et al. 2021, MNRAS, 506, 613, doi: 10.1093/mnras/stab1737
  * Sazonov et al. (2004) Sazonov, S. Y., Ostriker, J. P., & Sunyaev, R. A. 2004, MNRAS, 347, 144, doi: 10.1111/j.1365-2966.2004.07184.x
  * Sa̧dowski et al. (2015) Sa̧dowski, A., Narayan, R., Tchekhovskoy, A., et al. 2015, MNRAS, 447, 49, doi: 10.1093/mnras/stu2387
  * Schaerer (2003) Schaerer, D. 2003, A&A, 397, 527, doi: 10.1051/0004-6361:20021525
  * Schramm & Silverman (2013) Schramm, M., & Silverman, J. D. 2013, ApJ, 767, 13, doi: 10.1088/0004-637X/767/1/13
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Shang et al. (2010) Shang, C., Bryan, G. L., & Haiman, Z. 2010, MNRAS, 402, 1249, doi: 10.1111/j.1365-2966.2009.15960.x
  * Shapiro et al. (1999) Shapiro, P. R., Iliev, I. T., & Raga, A. C. 1999, MNRAS, 307, 203, doi: 10.1046/j.1365-8711.1999.02609.x
  * Shen et al. (2003) Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978, doi: 10.1046/j.1365-8711.2003.06740.x
  * Shull & van Steenberg (1985) Shull, J. M., & van Steenberg, M. E. 1985, ApJ, 298, 268, doi: 10.1086/163605
  * Sijacki et al. (2009) Sijacki, D., Springel, V., & Haehnelt, M. G. 2009, MNRAS, 400, 100, doi: 10.1111/j.1365-2966.2009.15452.x
  * Smith et al. (2017) Smith, A., Becerra, F., Bromm, V., & Hernquist, L. 2017, MNRAS, 472, 205, doi: 10.1093/mnras/stx1993
  * Stone & Pringle (2001) Stone, J. M., & Pringle, J. E. 2001, MNRAS, 322, 461, doi: 10.1046/j.1365-8711.2001.04138.x
  * Stone et al. (1999) Stone, J. M., Pringle, J. E., & Begelman, M. C. 1999, MNRAS, 310, 1002, doi: 10.1046/j.1365-8711.1999.03024.x
  * Sugimura et al. (2017) Sugimura, K., Hosokawa, T., Yajima, H., & Omukai, K. 2017, MNRAS, 469, 62, doi: 10.1093/mnras/stx769
  * Tagawa et al. (2020) Tagawa, H., Haiman, Z., & Kocsis, B. 2020, ApJ, 892, 36, doi: 10.3847/1538-4357/ab7922
  * Takahashi et al. (2013) Takahashi, S. Z., Inutsuka, S.-i., & Machida, M. N. 2013, ApJ, 770, 71, doi: 10.1088/0004-637X/770/1/71
  * Takasao et al. (2018) Takasao, S., Tomida, K., Iwasaki, K., & Suzuki, T. K. 2018, ApJ, 857, 4, doi: 10.3847/1538-4357/aab5b3
  * Takeo et al. (2020) Takeo, E., Inayoshi, K., & Mineshige, S. 2020, MNRAS, 497, 302, doi: 10.1093/mnras/staa1906
  * Takeo et al. (2018) Takeo, E., Inayoshi, K., Ohsuga, K., Takahashi, H. R., & Mineshige, S. 2018, MNRAS, 476, 673, doi: 10.1093/mnras/sty264
  * Toomre (1964) Toomre, A. 1964, ApJ, 139, 1217, doi: 10.1086/147861
  * Toyouchi et al. (2019) Toyouchi, D., Hosokawa, T., Sugimura, K., Nakatani, R., & Kuiper, R. 2019, MNRAS, 483, 2031, doi: 10.1093/mnras/sty3012
  * Toyouchi et al. (2021) Toyouchi, D., Inayoshi, K., Hosokawa, T., & Kuiper, R. 2021, ApJ, 907, 74, doi: 10.3847/1538-4357/abcfc2
  * Valentini et al. (2021) Valentini, M., Gallerani, S., & Ferrara, A. 2021, MNRAS, 507, 1, doi: 10.1093/mnras/stab1992
  * Valiante et al. (2018) Valiante, R., Schneider, R., Zappacosta, L., et al. 2018, MNRAS, 476, 407, doi: 10.1093/mnras/sty213
  * Venemans et al. (2017) Venemans, B. P., Walter, F., Decarli, R., et al. 2017, ApJ, 837, 146, doi: 10.3847/1538-4357/aa62ac
  * Verner et al. (1996) Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487, doi: 10.1086/177435
  * Visbal et al. (2015) Visbal, E., Haiman, Z., & Bryan, G. L. 2015, MNRAS, 453, 4456, doi: 10.1093/mnras/stv1941
  * Volonteri (2012) Volonteri, M. 2012, Science, 337, 544, doi: 10.1126/science.1220843
  * Volonteri & Rees (2005) Volonteri, M., & Rees, M. J. 2005, ApJ, 633, 624, doi: 10.1086/466521
  * Wang et al. (2021) Wang, F., Yang, J., Fan, X., et al. 2021, ApJ, 907, L1, doi: 10.3847/2041-8213/abd8c6
  * Wang et al. (2010) Wang, R., Carilli, C. L., Neri, R., et al. 2010, ApJ, 714, 699, doi: 10.1088/0004-637X/714/1/699
  * Wang et al. (2013) Wang, R., Wagg, J., Carilli, C. L., et al. 2013, ApJ, 773, 44, doi: 10.1088/0004-637X/773/1/44
  * Watarai et al. (2000) Watarai, K.-y., Fukue, J., Takeuchi, M., & Mineshige, S. 2000, PASJ, 52, 133, doi: 10.1093/pasj/52.1.133
  * Wechsler et al. (2002) Wechsler, R. H., Bullock, J. S., Primack, J. R., Kravtsov, A. V., & Dekel, A. 2002, ApJ, 568, 52, doi: 10.1086/338765
  * Whalen & Norman (2006) Whalen, D., & Norman, M. L. 2006, ApJS, 162, 281, doi: 10.1086/499072
  * Willott et al. (2010) Willott, C. J., Delorme, P., Reylé, C., et al. 2010, AJ, 139, 906, doi: 10.1088/0004-6256/139/3/906
  * Wise et al. (2019) Wise, J. H., Regan, J. A., O’Shea, B. W., et al. 2019, Nature, 566, 85, doi: 10.1038/s41586-019-0873-4
  * Wise et al. (2008) Wise, J. H., Turk, M. J., & Abel, T. 2008, ApJ, 682, 745, doi: 10.1086/588209
  * Wu et al. (2015) Wu, X.-B., Wang, F., Fan, X., et al. 2015, Nature, 518, 512, doi: 10.1038/nature14241
  * Yajima et al. (2017) Yajima, H., Ricotti, M., Park, K., & Sugimura, K. 2017, ApJ, 846, 3, doi: 10.3847/1538-4357/aa8269
  * Yan et al. (1998) Yan, M., Sadeghpour, H. R., & Dalgarno, A. 1998, ApJ, 496, 1044, doi: 10.1086/305420