https://orcid.org/0000-0002-0547-3208
ABSTRACT
By performing three-dimensional radiation hydrodynamics simulations, we study Bondi–Hoyle–Lyttleton accretion on to intermediate-mass black holes (BHs) wandering in the dusty gas. Here, we take into account the anisotropic radiation feedback and the sublimation of dust grains. Our simulations show that when the relative velocity between the BH and the gas is small (|$\sim 20\rm\, km~s^{-1}$|) and gas density is |$\sim 10^4 \rm cm^{-3}$|, the gas mainly accretes from near the equatorial plane of the accretion disc at a time-averaged rate of 0.6 per cent of the Bondi–Hoyle–Lyttleton rate. An ionized region like two spheres glued together at the equatorial plane is formed, and the dense shock shell appears near the ionization front. The BH is accelerated at |$\sim 10^{-8}\, \rm cm~s^{-2}$| due to the gravity of the shell. For denser gas (|$\sim 10^6 \rm cm^{-3}$|), the time-averaged accretion rate is also 0.6 per cent of the Bondi–Hoyle–Lyttleton rate. However, the BH is decelerated at |$\sim 10^{-7}\, \rm cm~s^{-2}$| due to gravity of the dense downstream gas although the dense shock shell appears upstream. Our simulations imply that intermediate-mass BHs in the early universe keep floating at |$\gtrsim {\rm several}\, 10\, \rm km~s^{-1}$| without increasing mass in interstellar gas with density of |$\sim 10^4\, \rm cm^{-3}$|, and slow down and grow into supermassive BHs in galaxies with the density of |$\sim 10^6\, \rm cm^{-3}$|.
1 INTRODUCTION
Bondi–Hoyle–Lyttleton accretion is the accretion phenomenon that occurs under the condition where there is relative motion between the central object and the surrounding gas, in which the gas is attracted by gravity and accretes to the central object (for recent reviews, see Edgar 2004). It occurs when the gas accretes to neutron stars and stellar-mass black holes (BHs) which are kicked by supernova explosions (Hobbs et al. 2005; Wong, Willems & Kalogera 2010; Janka 2012; Wongwathanarat, Janka & Müller 2013; Tauris et al. 2017; Kapil et al. 2023), to stars and BHs in globular clusters (Zwart & McMillan 2002; Gurkan, Freitag & Rasio 2004; Kaaz, Antoni & Ramirez-Ruiz 2019; Gonzalez et al. 2021; Shi, Grudić & Hopkins 2021), to compact objects in X-ray binaries via stellar wind accretion, and to intermediate-mass BHs (IMBHs) floating in the remnant galactic disc (Matteo, Springel & Ilernquist 2005; Dubois et al. 2012; Matteo et al. 2012, 2017; Sijacki et al. 2015; Latif, Volonteri & Wise 2018).
Bondi–Hoyle–Lyttleton accretion flow was investigated through numerous hydrodynamics simulations from the 1970s to the late 1990s (Hunt 1971; Shima et al. 1985; Fryxell et al. 1988; Ho et al. 1989; Matsuda et al. 1991; Ruffert 1994; 1996), which showed that numerically obtained the accretion rates are consistent with the analytical estimates of Hoyle & Lyttleton (1939) and Bondi & Hoyle (1944). Subsequently, numerical simulations considering the radiation from the central object were performed (Blondin et al. 1990; Taam et al. 1991; Milosavljevic, Couch & Bromm 2008; Park & Ricotti 2012; Park & Bogdanović 2017; Sugimura & Ricotti 2020; Toyouchi et al. 2020), showing that the accretion rates are significantly reduced by radiation feedback e.g. the radiation force and ionization heating. In the case in which the central object is compact object with an accretion disc, however, the accretion rates cannot decrease significantly because the gas accretes through the region near the disc plane, where the radiation effect is relatively weak. Although Bondi–Hoyle–Lyttleton accretion on to the compact objects under the anisotropic radiation field produced by the accretion disc has been investigated in previous studies (Fukue & Ioroi 1999; Hanamoto, Ioroi & Fukue 2001; Ogata, Ohsuga & Yajima 2021), the effect of fluid dynamics such as gas pressure gradient force has not been taken into account. Therefore, radiation hydrodynamics simulation of Bondi–Hoyle–Lyttleton accretion including an anisotropic radiation field is required.
The increase in the absorption coefficient due to dust also decreases the accretion rate because the radiation force becomes more effective. In fact, Yajima et al. (2017) and Toyouchi et al. (2020) showed that accretion on to the central BHs is suppressed by the strong radiation force acting on dusty gas through radiation hydrodynamics simulations. Toyouchi et al. (2020) revealed that the accretion rate decreases with increasing metallicity (dust amount) and is comparable to the Eddington accretion rate for dusty gas. In these previous studies, dust sublimation was not considered, and isotropic radiation fields were adopted. More realistically, dust sublimation should also be taken into account, because the radiative force is reduced in dust sublimation regions.
In this study, we perform three-dimensional radiation hydrodynamics simulations to investigate the effects of radiation anisotropy and dust sublimation on the Bondi–Hoyle–Lyttleton accretion mechanism. We focus on the growth process of supermassive BHs observed in the early universe and select IMBHs floating in dusty gas with metalicity Z = 0.1Z⊙ as our target objects. The rest of this paper is organized as follows: We describe the method of our three-dimensional radiation hydrodynamics simulations in Section 2. The simulation results are given in Section 3, and we discuss the evolution of IMBHs in Section 4. Finally, we provide a summary and conclusion in Section 5.
2 SIMULATION METHODS
In this study, we perform three-dimensional radiation hydrodynamics simulations with sfumato-m1 (Fukushima & Yajima 2021). This code is based on SFUMATO-RT (Sugimura et al. 2020), which is the modified version of sfumato (Matsumoto 2007; Matsumoto, Dobashi & Shimoikura 2015). In sfumato-m1, the module with the M1-closure technique is adapted instead of the adoptive ray-tracing solver developed in Sugimura et al. (2020). We further include the anisotropic radiation field and dust sublimation for the current purpose.
2.1 Basic setup
We perform the three-dimensional radiation hydrodynamics simulations with the Cartesian coordinate and simulate the gas flow relative to the BH. The gas has constant number density n∞, positive x-direction velocity v∞, and temperature T∞ = 180 K at the upstream boundary. We set the size of a calculation box as Rout = 12.6 pc, which is much larger than Bondi–Hoyle–Lyttleton radius and the size of ionized region around a BH in our simulations.
We adopt the sink method to mask the accretion disc and set the sink radius as |$R_{\rm sink}=2.7\times 10^{-3}~\rm pc$|, which is small enough to resolve dust sublimation region. We use the nested grid method to resolve the gas structure inside ionized region and dust sublimation region efficiently. In our simulations, the maximum refinement level is lmax = 10, and the minimum cell size is |$\Delta _{10}=3.85\times 10^{-4}~\rm pc$|. The BH is fixed at the origin and grows due to mass accretion. However, the increase in the BH mass is negligibly small in the elapsed time of our simulations. Also, we neglect acceleration of the BH, self-gravity of gas, and magnetic field for simplicity.
In the model with the anisotropic radiation, we perform the convergence check tests with the higher refinement level lmax = 11. As a result, we confirm good convergence in terms of the flow structure and the accretion rate. We also perform the simulations with the same sink size and resolution of Toyouchi et al. (2020), where they investigated Bondi–Hoyle–Lyttleton accretion of the dusty gas around luminous object by three-dimensional radiation hydrodynamics simulations. We confirm that our simulation results are consistent with their study.
The M1-closure method employed in this study is known to be more diffusive than the ray-tracing method (e.g. Asahina et al. 2020). By the test calculation of an ionized bubble produced by the anisotropic radiation, it is found that the distance from the radiation source to the ionization front using the M1-method is 1.2–1.5 times larger (smaller) than the analytical estimation around the disc equatorial plane (rotation axis). This discrepancy may be due not only to differences in the calculation methods for radiative transfer but also to differences in the treatment of photons emitted by recombination. Our method also takes into account the recombination photons, but the Case-B recombination is assumed in the analytical solution. The expansion of the ionized region due to recombination photons is more pronounced near the equatorial plane because many recombination photons enter from the ionized regions sandwiching the equatorial plane. Simulations using more accurate radiation transfer methods are left as an important future work (e.g. Sugimura & Ricotti 2020).
2.2 Basic equations
We solve the following governing equations: the equation of continuity,
the equation of motion,
and the equation of energy,
where E is total energy defined as
Here, ρ, |$\boldsymbol{v}$|, and P are the gas density, gas velocity, and gas pressure, respectively, and |$\boldsymbol{g}$| and |$\boldsymbol{f}$| are the gravity and radiation force, and Γ and Λ are the specific heating and cooling rates. We estimate the adiabatic exponent γ as in Omukai & Nishi (1998). In this simulations, we also solve the transfer of far-ultraviolet (FUV; 11.2– 13.6 eV) and extreme-ultraviolet (EUV; 13.6 eV–1 keV) photons emitted from the sink using the M1-closure approximation method. We also compute the transfer of infrared (IR) photons mainly emitted from dust grains as thermal emission. Further details on our radiative transfer method are described in Fukushima & Yajima (2021).
We calculate the non-equilibrium chemical reactions for the eleven species of |$\rm H, H_2, H^-, H^+, H_2^+, e, CO, C_{II}, O_I, O_{II}$|, and |$\rm O_{III}$|. The number density of the i-th species is calculated as
where yi = ni/nH is the number ratio of i-th species to hydrogen nuclei, and Ri is the corresponding chemical reaction rate. We also consider the dust grains in the gas, assuming the dust-to-gas mass ratio of 0.01 × Z/Z⊙ = 10−3. The temperature of dust grains is estimated from the energy balance between the absorption/emission of radiation and energy transfer with the gas. With the updated chemical abundances, we compute Γ and Λ for solving the energy equation. The details of them are shown in Fukushima & Yajima (2021). We also include dust sublimation. We set the temperature of dust sublimation as |$10^3~\rm K$|, above which we ignore the contributions of dust grains on radiation, force attenuation, and thermal processes of gas components.
2.3 Radiative source
In our simulations, the BH luminosity L is estimated with the mass accretion rate |$\dot{M}$|. We adapt the model of Watarai et al. (2000) as
where LE (≡ 4πcGMBH/κes) and |$\dot{M}_{\rm E}~(\equiv L_{\rm E}/\eta c^2)$| are the Eddington luminosity and Eddington accretion rate. Here, MBH, |$\kappa _{\rm es}=0.4~\rm cm^2~g^{-1}$| and η = 0.1 are the BH mass, opacity of electron scattering, and the radiative efficiency. We assume that the radiation is emitted from the sink region with a single power-law spectrum Lν ∝ ν−α in the frequency range 11.2 eV < hν < 1 keV, where h and ν are planck constant and photon frequency, respectively. Here, we adopt the power-law index of α = 1.5, corresponding to the spectral energy distribution of active galactic nuclei in the high-accretion state (Sazonov, Ostriker & Sunyaev 2004).
We inject the extreme EUV and FUV photons at the sink region with the directional dependence as described below. The number of photons injected per unit time in a single cell is given as
where |$\dot{N}_{\nu }^{\rm sink}$| is the photon emissivity of the BH, and k and j are the cell numbers. The anisotopy factor |$\mathcal {F}_j(\Theta , R)$| represents the direction dependence of the radiation as
Here, the angle from the rotation axis of the disc to the centre of each cell (Θ) is defined as:
where |$R(=\sqrt{x^2+y^2+z^2})$| is the distance of the origin, φ is the angle between y-axis and the disc rotation axis projected on to yz-plane, and ψ is the angle between x-axis and the disc rotation axis, respectively. We set φ = 0 in the present simulations. We also solve the diffuse IR photons produced as the dust thermal emission.
2.4 Cases examined
In our simulations, we fix the initial BH mass and the metallicity as MBH = 104 M⊙ and Z = 0.1 Z⊙. The simulation parameters are summarized in Table 1. Hereafter, we label each model according to the anisotropy factor (|$\mathcal {F}_j$|), the number density (n∞), and gas velocity (v∞) considered in each simulation. For example, ‘isoN4V20’ represents the simulation model with the isotropic radiation field, |$n_{\infty }=10^4\rm ~cm^{-3}$|, and |$v_{\infty }=20\rm ~km~s^{-1}$|.
Model parameters with the results of the time-averaged accretion rate and the luminosity.
| Model | |$n_{\infty }(\rm cm^{-3}$|) | |$v_{\infty }(\rm km~s^{-1})$| | Radiation field | |$\overline{{\dot{M}}}~({\rm M}_{\odot }\,\rm yr^{-1})$| | |$\overline{L}~(\mathrm{ L}_{\odot })$| | |$\overline{{\dot{M}}}/\dot{M}_{\rm BHL}$| |
|---|---|---|---|---|---|---|
| isoN4V20 | 104 | 20 | isotropic | 4.3 × 10−6 | 6.3 × 106 | 6.2 × 10−3 |
| edgeN4V20 | 104 | 20 | edge-on | 4.4 × 10−6 | 6.5 × 106 | 6.4 × 10−3 |
| edgeN6V20 | 106 | 20 | edge-on | 4.2 × 10−4 | 4.6 × 108 | 6.1 × 10−3 |
| edgeN4V100 | 104 | 100 | edge-on | 4.4 × 10−6 | 6.5 × 106 | 5.6 × 10−1 |
| poleN4V20 | 104 | 20 | pole-on | 4.2 × 10−6 | 6.2 × 106 | 6.1 × 10−3 |
| Model | |$n_{\infty }(\rm cm^{-3}$|) | |$v_{\infty }(\rm km~s^{-1})$| | Radiation field | |$\overline{{\dot{M}}}~({\rm M}_{\odot }\,\rm yr^{-1})$| | |$\overline{L}~(\mathrm{ L}_{\odot })$| | |$\overline{{\dot{M}}}/\dot{M}_{\rm BHL}$| |
|---|---|---|---|---|---|---|
| isoN4V20 | 104 | 20 | isotropic | 4.3 × 10−6 | 6.3 × 106 | 6.2 × 10−3 |
| edgeN4V20 | 104 | 20 | edge-on | 4.4 × 10−6 | 6.5 × 106 | 6.4 × 10−3 |
| edgeN6V20 | 106 | 20 | edge-on | 4.2 × 10−4 | 4.6 × 108 | 6.1 × 10−3 |
| edgeN4V100 | 104 | 100 | edge-on | 4.4 × 10−6 | 6.5 × 106 | 5.6 × 10−1 |
| poleN4V20 | 104 | 20 | pole-on | 4.2 × 10−6 | 6.2 × 106 | 6.1 × 10−3 |
Model parameters with the results of the time-averaged accretion rate and the luminosity.
| Model | |$n_{\infty }(\rm cm^{-3}$|) | |$v_{\infty }(\rm km~s^{-1})$| | Radiation field | |$\overline{{\dot{M}}}~({\rm M}_{\odot }\,\rm yr^{-1})$| | |$\overline{L}~(\mathrm{ L}_{\odot })$| | |$\overline{{\dot{M}}}/\dot{M}_{\rm BHL}$| |
|---|---|---|---|---|---|---|
| isoN4V20 | 104 | 20 | isotropic | 4.3 × 10−6 | 6.3 × 106 | 6.2 × 10−3 |
| edgeN4V20 | 104 | 20 | edge-on | 4.4 × 10−6 | 6.5 × 106 | 6.4 × 10−3 |
| edgeN6V20 | 106 | 20 | edge-on | 4.2 × 10−4 | 4.6 × 108 | 6.1 × 10−3 |
| edgeN4V100 | 104 | 100 | edge-on | 4.4 × 10−6 | 6.5 × 106 | 5.6 × 10−1 |
| poleN4V20 | 104 | 20 | pole-on | 4.2 × 10−6 | 6.2 × 106 | 6.1 × 10−3 |
| Model | |$n_{\infty }(\rm cm^{-3}$|) | |$v_{\infty }(\rm km~s^{-1})$| | Radiation field | |$\overline{{\dot{M}}}~({\rm M}_{\odot }\,\rm yr^{-1})$| | |$\overline{L}~(\mathrm{ L}_{\odot })$| | |$\overline{{\dot{M}}}/\dot{M}_{\rm BHL}$| |
|---|---|---|---|---|---|---|
| isoN4V20 | 104 | 20 | isotropic | 4.3 × 10−6 | 6.3 × 106 | 6.2 × 10−3 |
| edgeN4V20 | 104 | 20 | edge-on | 4.4 × 10−6 | 6.5 × 106 | 6.4 × 10−3 |
| edgeN6V20 | 106 | 20 | edge-on | 4.2 × 10−4 | 4.6 × 108 | 6.1 × 10−3 |
| edgeN4V100 | 104 | 100 | edge-on | 4.4 × 10−6 | 6.5 × 106 | 5.6 × 10−1 |
| poleN4V20 | 104 | 20 | pole-on | 4.2 × 10−6 | 6.2 × 106 | 6.1 × 10−3 |
In this study, ψ = π/2 (edge-on) is mainly adopted. The gas inflow direction (−x direction) at the upstream boundary is perpendicular to the rotation axis of the disc in this model. Such a situation appears naturally. The angular momentum vector of the gas injected from the upstream boundary relative to the origin (BH) is perpendicular to the x-axis. Therefore, although the gas density at the upstream boundary is assumed to be uniform in our simulations, if the density is non-uniform, the total angular momentum vector of the gas will be perpendicular to the x-axis. Thus, the gas falling into the central region would form the accretion disc with the rotation axis perpendicular to the x-axis. Even if the rotation axis of the accretion disc is not perpendicular to the x-axis at the beginning, it should eventually become perpendicular to the x-axis since the gas in the initial disc is swallowed by the BHs, and the supplied gas becomes main component of the disc. The pole-on situation (ψ = 0) is realized only at the moment if the black hole happens to enter the gas cloud from the direction of the rotation axis of the disc, but it is adopted in order to compare with the edge-on models. In addition, for comparison with the model of anisotropic radiation, we perform the simulation with isotropic radiation.
We adopt the model with the gas number density |$n_{\infty }=10^4\rm ~cm^{-3}$| and velocity |$v_{\infty }=20\rm ~km~s^{-1}$| as the fiducial model. We additionally investigate the cases with high-velocity |$v_{\infty }=100~\rm km~s^{-1}$| and high-density |$n_{\infty }=10^6~\rm cm^{-3}$| only in the models with edge-on anisotropic radiation. Such high-density environments would appear, for example, in merger remnant galaxies (e.g. Fiacconi et al. 2013; Lima et al. 2017). Also, Katz et al. (2023) has reported the presence of relatively dense gas, |$\gtrsim 10^4 \rm cm^{-3}$|. The high-velocity situation, |$v_\infty \sim 100~\rm km~s^{-1}$|, might be realized when galaxy mergers occur (e.g. Mayer et al. 2007).
3 RESULTS
3.1 Overview
Before examining the numerical results in more detail, we overview the flow structure and the accretion rates to BHs. Fig. 1 schematically shows a flow structure around a BH. Ionizing photons emitted from the accreting BH create an ionized region over several pc. The shock structure is formed upstream (−x direction) if the relative velocity of the ambient gas to the ionization front is between the critical values of the D-type and R-type ionization fronts (vD < v∞ < vR), as also shown in Park & Ricotti (2013). Inside the shocked region, the gas slows to vD at the ionization front. If the velocity v∞ is higher than the critical value for the R-type ionization front vR, the shock does not appear. Dust grains are heated and sublimate close to the BH. In this region, radiation force cannot push out the gas, and the gas tends to accrete to the BH. In the outer region, where the dust grains are not sublimated, the gas is accelerated by radiation force and thermal pressure. Then, it tends to pass through the ionized region without falling into the BH.
A schematic picture of a gas flow around a BH. We show a case with an edge-on radiation field. The BH floats in the interstellar medium. The filled circle represents the sink region. The solid and dotted lines show streamlines that pass through or accrete to BH. The shapes of ionized and dust sublimation regions are different in the models with isotropic and pole-on radiation fields.
Fig. 2 shows the time evolution of the mass accretion rate. The black dashed line means the Eddington accretion rate for dusty gas |$\dot{M}_{\rm E,dg}~(\equiv \kappa _{\rm es}\dot{M}_{\rm E}/(\kappa _{\rm es}+\kappa _{\rm d}))$|, where the dust opacity for UV light κd is given as |$\kappa _{\rm d}=2.8\times 10^2 (\mathit{ Z}/Z_{\odot })~\rm cm^2~g^{-1}$| (e.g. Yajima et al. 2017). The time-averaged accretion rates |$\overline{{\dot{M}}}$| and luminosity |$\overline{L}$| obtained in our simulations are summarized in Table 1. Here, we assume that the luminosity L is the function of the accretion rate |$\dot{M}$| (see equation 6). In the isotropic radiation model ‘isoN4V20’, the accretion rate oscillates periodically between |$4.5\times 10^{-2}\dot{M}_{\rm E}$| and |$1.8\times 10^{-3}\dot{M}_{\rm E}$| with a period of |$5.5\times 10^{-4}~\rm Myr$|. The time-averaged accretion rate is about |$2\times 10^{-2}\dot{M}_{\rm E}$| (4.3 × 10−6 M⊙ yr−1). The accretion bursts (rapid increase in accretion rate) occur in the model ‘edgeN6V20’ with a period of ∼0.01 M yr. The time-averaged accretion rate is about |$2\dot{M}_{\rm E}$| (4.2 × 10−4 M⊙ yr−1). Although it does not appear in the figure, model ‘edgeN4V20’ also exhibits accretion bursts at the interval with ∼0.15 Myr (we will discuss later). The time-averaged value in this model, |$2 \times 10^{-2} \dot{M}_{\rm E}$| (4 × 10−6 M⊙ yr−1, see Table 1), is slightly larger than the accretion rate shown in Fig. 2 because the bursts increase the time-averaged rate. The quasi-steady state is achieved in the models ‘edgeN4V100’ and ‘poleN4V20’, and the time-averaged rate is almost the same as in model ‘edgeN4V20’.
Time evolution of mass accretion rates. The vertical axis is normalized by the Eddington accretion rate |$\dot{M}_{\rm E}$| estimated with the opacity of electron scattering. The dotted lines represent the Eddington accretion rate for dusty gas |$\dot{M}_{\rm E,dg}$|.
3.2 Isotropic case
Fig. 3 shows a flow structure of model ‘isoN4V20’ around the BH. We find that an ionized region (nH +/nH ∼ 1) extends up to x ∼ 5 pc, and a dense shock shell is formed at the upstream ionization front (x ∼ −0.5 pc). The gas velocity gradually increases to ∼50 km s−1 in the ionized region after passing through the dense shock shell. Such acceleration is caused by the thermal pressure increased by ionization heating. Increasing the gas velocity and increasing the gas temperature reduce the accretion rate since the Bondi–Hoyle–Lyttleton accretion rate decreases when the velocity and temperature are high. The outward radiation force also woks to reduce the accretion rate. Thus, the time-averaged accretion rate, |$2\times 10^{-2}\dot{M}_{\rm E}$|, is three orders of magnitude smaller than the classical Bondi–Hoyle–Lyttleton accretion rate |$\dot{M}_{\rm BHL}$|. We will show the comparison with previous work (Toyouchi et al. 2020) that does not include the dust sublimation in Section 4.2.
Gas flow structure on a |$1~\rm pc$| scale for the isotropic radiation model (‘isoN4V20’). Each panel shows distributions of the gas number density (upper left), velocity (upper right), and ionization degree (lower right). In the lower-left panel, we show the radiation force normalized by gravity in each cell. The white arrows represent streamlines to the +x direction.
The reason why the periodic oscillation of the accretion rate occurs can be understood from Fig. 4. This figure shows the time evolution of flow structure and the radiation force normalized by gravity on the |$10^{-2}~\rm pc$| scale in one period of oscillation. The time evolution of the accretion rate is also plotted. When the gas accretion rate is less than |$\dot{M}_{\rm E,dg}$|, the gravity overcomes the radiation force in the entire region (Fig. 4-1). The gas moves towards the BH, and the flow structure becomes similar to Bondi–Hoyle–Lyttleton accretion. When the accretion rate exceeds |$\dot{M}_{\rm E,dg}$| (Fig. 4-2), the radiation force reduces the accretion rate. The reason why gas accretion does not stop completely is because the radiation force is not effective in the region where dust sublimes near the BH (inside the black solid lines). When the accretion rate decreases to |$\dot{M}_{\rm E,dg}$|, gravity overcomes radiation force. Thus, the accretion rate increases again (Fig. 4-3) and then begins to decrease (Fig. 4-4) due to the strong radiation force.
Time evolution of the gas flow in one oscillation period for the isotropic radiation model (‘isoN4V20’). Each snapshot is the same as the bottom left panel in Fig. 3. The solid lines and filled circles represent the dust sublimation and sink regions. The upper right panel shows the time evolution of the accretion rate in one period of oscillation. The filled circles mark the epochs for each snapshot 1–4.
3.3 Edge-on case
3.3.1 Fiducial model
In the case of ‘edgeN4V20’ (fiducial model), the accretion rate is almost constant for most of the time (quiescent phase), but accretion bursts occur periodically (burst phase). Fig. 5 presents the flow structure in the meridian (z = 0) and equatorial (y = 0) planes in the quiescent phase. Although the radiation field is anisotropic, the profiles of density and velocity on the ∼1 pc scale are almost similar to that of the isotropic radiation model (see Fig. 5a). Also in this model, an ionized region extends up to several pc. Around the surface of the dense shock shell, the gas pressure gradient force by the high-density and high-temperature gas in the shell moves the gas away from the BH. Indeed, we find that some gas moves away from the x-axis along the shell (see streamlines). The effect of anisotropic radiation can be seen in Fig. 5(b). It is found that the shape of the ionization region is like two spheres glued together at the equatorial plane. The dense shock shell appears near the ionization front. The velocity and temperature of the gas that passes through the shell and enters the ionized region rapidly increases by gas pressure gradient force and the ionization heating. These increases in velocity and temperature work to reduce the rate of accretion by the Bondi–Hoyle–Lyttleton mechanism. We find in Fig. 5(c) that the gas accretes mainly from the equatorial plane since the radiation force is ineffective because BH irradiation is considerably weakened by absorption. In contrast, the strong radiation force works to prevent gas accretion around the rotation axis of the disc. As a result, the time-averaged accretion rate in the quiescent phase is much smaller than Bondi–Hoyle–Lyttleton rate (|$\sim 0.3~{{\ \rm per\ cent}}$| of |$\dot{M}_{\rm BHL}$|), and nearly comparable to |$\dot{M}_{\rm dg}$|.
Same as Fig. 3 but for the edge-on fiducial model (‘edgeN4V20’). Each panel shows the flow structure on (a) |$1~\rm pc$|, (b) 10−1, and (c) |$10^{-2}~\rm pc$| scales at |$t=0.2~\rm Myr$| (|$0.04~\rm Myr$| before the next accretion burst). Each column represents the flow structure in the meridian (x y-) and equatorial (x z-) planes from top to bottom.
The accretion bursts only increase the time-averaged accretion rate by about a factor of 2, and also hardly affect acceleration at all (we will discuss later). The bursts occur when part of the dense shock shell (near the x-axis) flows into the sink. The dense shock shell is almost never moved due to the balance between the ram pressure gradient force and the gas pressure gradient force in the quiescent phase. However, part of the dense shock shell around the x-axis gradually moves inward due to the gravity of the BH and then flows into the sink leading to the burst phase. The accretion rate of the burst phase is |$\sim \dot{M}_{\rm E}$| and about half of the total accreting gas accretes during the burst phase. Due to the bursts, the time-averaged accretion rate (Table 1) is approximately twice the accretion rate during the quiescent phase (see Fig. 2).
The burst interval, |$\sim 0.15\,\rm Myr$|, can be understood approximately as follows. In the situation where gas pressure and ram pressure are approximately balanced at the dense shock shell, the gravity begins to pull the shell when the column number density of the shell (Nshell) exceeds |$n_{\infty } r_{\rm shell}^2 v_{\infty }^2/GM_{\rm BH}$| because the gravity, |$GM_{\rm BH}N_{\rm shell}/r_{\rm shell}^2$|, becomes larger than the ram pressure (gas pressure), |$n_{\infty } v_{\infty }^2$|. Since the distance of the dense shock shell from the BH, rshell, is about |$0.5\rm\, pc$|, the critical column number density obtained from the above relation is |$\sim 10^{23}\, \rm cm^{-2}$|, which is roughly consistent with that the burst in the present simulation. The time-scale on which the column number density of the dense shock shell reaches |$10^{23}\, \rm cm^{-2}$| is approximately |$\left(n_{\infty } r_{\rm shell}^2 v_{\infty }^2/GM_{\rm BH}\right)/(n_{\infty } v_{\infty })\sim v_{\infty } r_{\rm shell}^2/GM_{\rm BH}$|. This time-scale is evaluated as |$\sim 1.1\, \rm Myr$| and is longer than the free-fall time so the burst interval is determined by this time-scale.
3.3.2 Higher velocity model
Fig. 6 shows the flow structure in the model ‘edgeN4V100’. The most noticeable difference from the fiducial model is the absence of shock waves. This is because the relative velocity is largar than vR. We also confirm that the density does not change so much at the ionization front. The density change at the ionization front can be estimated analytically with the mass and momentum conservation laws as |${\rho _{\rm HII}}/{\rho _{\infty }} \sim {v_{\infty }^2} \left(1- \sqrt{1-{4c_{\rm s,HII}^2}/{v_{\infty }^2}} \right)/{2c_{\rm s,HII}^2}$| under the condition cs,HII ≫ c∞, where ρHII, cs,HII, and c∞ are the gas density of ionized gas and sound speeds of ionized and neutral gas (Spitzer 1978). If v∞ ≫ cs,HII, then we obtain ρHII ∼ ρ∞. We confirm in our simulations that the condition of v∞ ≫ cs,HII is satisfied. The somewhat low-density region appears over several pc downstream (x > 0, y ≲ 1 pc, dark-purple in Fig. 6a). This is because the gas is accelerated by radiation force and gas pressure gradient force in the ionized region.
Same as Fig. 5 but for the edge-on high-velocity model (‘edgeN4V100’). Each panel shows the flow structure on (a) 1 pc and (b) 10−2 pc scales. The solid lines in panel (b) represent the dust sublimation regions.
The radiation force does not affect the the motion of gas at distances within the Bondi–Hoyle–Lyttleton radius from the x-axis so that the flow structure is similar to the classical Bondi–Hoyle–Littleton accretion. Fig. 6(b) shows the flow structure in the 10−2 pc scale. Around the disc rotation axis, the radiation force is basically stronger than gravity, but only near the BH is it weaker than the gravity due to sublimation of the dust grains. Since the Bondi–Hoyle–Lyttleton radius is around 9 × 10−3 pc, and since the dust sublimation radius around the rotation axis is ∼7 × 10−3 pc, the radiation force does not affect the motion of gas at distances within the Bondi–Hoyle–Lyttleton radius from the x-axis. Thus, a quasi-steady flow similar to Bondi–Hoyle–Littleton accretion appears. The accretion rate is roughly comparable to |$\dot{M}_{\rm BHL}$| (see Table 1).
3.3.3 Dense model
Model ‘edgeN6V20’ differs from the fiducial model in that it produces more intense, shorter-interval accretion bursts. Fig. 7 shows the time evolution of the gas density distribution in the meridian plane (z = 0) for the model ‘edgeN6V20’. As shown in Fig. 2, the accretion bursts occur periodically in this model. The burst mechanism is similar to that of fiducial model. In the quiescent phase, the dense shock shell exists near the ionization front, and it gradually becomes denser. After a while, the gravity causes part of the shell to reach the sink, and the burst phase begins (see Fig. 7a). In the burst phase, the accretion rate exceeds the Eddington limit. However, the radiation force is ineffective around the equatorial plane due to the attenuation of radiation and the thermal pressure of inoized bubbles cannot push outward these dense gas. After the significant part of the dense shock shell falls into the BH through the vicinity of the equatorial plane, it transitions to the quiescent phase (see Fig. 7b). After that, the part of the burst shell falls into the BH, and the burst accretion occurs again. The above process is repeated, resulting in periodic bursts. The interval of the accretion burst, |$\sim 0.01\rm\, Myr$|, is much shorter than that of the fiducial model, |$\sim 0.15\rm\, Myr$|. This is because the rshell in this model, ∼0.1 pc, is smaller than that of the fiducial model (it was shown in Section 3.3.1 that the burst interval depends on the rshell).
The flow structure on |$10^{-1}~\rm pc$| scale for the edge-on high-density model (‘edgeN6V20’). Each panel shows the number density distribution in the equatorial (xy-) plane at (a) the burst and (b) quiescent phases. We also show the mass accretion rates and the elapsed times in each panel.
At the burst phase, the radius of the ionized region, RHII, is comparable to or slightly smaller than the Bondi–Hoyle–Lyttleton radius, RBHL, at around the equatorial plane where gas mainly accretes. On the other hand, we find RHII > RBHL in regions other than the equatorial plane. This means that the condition for super-Eddington accretion ‘RHII < RBHL’ (Inayoshi, Haiman & Ostriker 2016) is satisfied locally. The effect of dust, which increases the radiation force and narrows the ionized region, on the condition for super-Eddington accretion should be investigated in detail in the future.
3.4 Pole-on case
Figs 8(a) and (b) show the flow structure on a 1 and 10−2 pc scales of the pole-on radiation model (‘poleN4V20’). The flow structure is quasi-steady and axisymmetric with respect to the x-axis. The structure on pc scale is almost the same as the fiducial model. The ionized region extends to several pcs, and the dense shock shell is formed at the ionization front (see Fig. 8a). As shown in Fig. 8(b), this model also shows that gas accretion occurs mainly from near the equatorial plane (x = 0 plane), where the gravity is stronger than the radiation force. In addition, the gas that is attracted by the gravity but does not fall into the BH reaches near the x-axis at the downstream region and then flows out in the +x direction. Unlike the fiducial model, the vortex flows appear in the upstream region. This is due to the radiation force pushing the gas in the −x direction.
Same as Fig. 5 but for the pole-on model (‘poleN4V20’). Each panel shows the flow structure on (a) |$1~\rm pc$| and (b) |$10^{-2}~\rm pc$| scales. We only show the snapshot in the xy-plane due to the axisymmetry of flow.
4 DISCUSSION
4.1 Evolution of IMBHs
Gravity by non-uniform gas distribution around the BHs and momentum transport to the BHs by gas accretion induce acceleration of BHs. Fig. 9 shows the acceleration in the x-direction ax(R) at the elapsed time |$t = 8.5\times 10^{-2}~\rm Myr$|, which is calculated as
where dΩ and |$\mathrm{ d}\boldsymbol{S}_{\rm sink}$| are the solid angle and area vector. Here, a negative (positive) value of ax(R) means that the BH accelerates in the upstream (downstream) direction and the velocity of BH relative to the interstellar gas increases (decreases).
Acceleration in the x-direction of BHs at |$t\sim 8.5\times 10^{-2}~\rm Myr$| as a function of R. The solid (dotted) lines mean that the acceleration is negative (positive) and the speed of the BHs increases (decreases). The line colours are the same as in Fig. 2. The star markers correspond to the positions of the ionization fronts.
It is found that the acceleration ax(Rout) is about |$-10^{-8}\rm\, cm~s^{-2}$| in the models with v∞ = 20 km s−1 and |$n_{\infty }=10^4~\rm cm^{-3}$| (‘isoN4V20’, ‘poleN4V20’, and ‘edgeN4V20’). This is mainly due to gravity from the dense shock shell on the ionization front. This can be understood from the fact that the acceleration suddenly increases at around the star marks (the positions of the ionization front) and becomes nearly constant outside of them. In the model with the high velocity v∞ = 100 km s−1 (‘edgeN4V100’), we find |$a_x (R_{\rm out}) \sim 10^{-10}\rm\, cm~s^{-2}$|, which is consistent with the result of Ostriker (1999), in which the acceleration of the central object of Bondi–Hoyle–Lyttleton flow has been evaluated. Positive acceleration is due to the accretion of gas with positive momentum (ρvx > 0) on to the BH. Although the gas density outside the sink in the upstream region is higher than that in the downstream (see Section 3.3.2), the negative acceleration due to this density difference is smaller than the positive acceleration via the gas accretion. In the high-density case (’edgeN6V20’), the acceleration in the burst phase is about |$a_x (R_{\rm out})\sim 10^{-7}\rm cm~s^{-2}$|, and we find the gravity of the dense downstream gas induces the positive acceleration. Although the acceleration becomes negative due to the gravity of the dense shock shell in the quiescent phase, time-averaged acceleration is positive, |$\sim 10^{-7}\rm cm~s^{-2}$|.
In the following, we discuss the evolution of the mass and velocity of the IMBHs (MBH ∼ 104 M⊙) floating in early galaxies (z ≳ 6) using our results of model ‘edgeN4V20’, ‘edgeN4V100’, and ‘edgeN6V20’. Here, we use a(Rout) at |$t\sim 0.1\rm Myr$| presented in Fig. 9 instead of the time-averaged value since the acceleration at larger distance (R ≳ several pc) does not change with time as much.
In the case that the gas density is |$\sim 10^4\rm ~cm^{-3}$|, the IMBHs continue to float without significant mass growth. For the model ‘edgeN4V100’, the acceleration is |$\sim 10^{-10}\rm ~cm~s^{-2}$| so that the time-scale of the slowdown is |$\sim 3~\rm Gyr$|. This is much longer than the age of the universe at z = 6 (|$\sim 900~\rm Myr$|). This implies that the speed of IMBHs does not change so much. On the other hand, the acceleration time-scale of model ‘edgeN4V20’, |$\sim 6~\rm Myr$|, is much shorter than |$900~\rm Myr$|. Thus, if the initial velocity is relatively small, |$\sim 10 ~ \rm km~s^{-1}$|, the IMBHs will continue to speed up until the velocity reaches |${\rm several} \times 10 ~ \rm km~s^{-1}$| where the acceleration time-scale becomes comparable to the age of the universe. The above discussion does not take into account the change in the mass of the IMBHs, but this is appropriate. In both model ‘edgeN4V100’ and ‘edgeN4V20’, the time-scale of mass growth derived from the accretion rate (|$\sim 4\times 10^{-6}~{\rm M}_{\odot }~\rm yr^{-1}$|) is |$\sim 2~\rm Gyr$|, which is much longer than |$900~\rm Myr$|. Thus, we conclude that the IMBHs continue to float at the velocity of |$\gtrsim {\rm several} \times 10 ~ \rm km~s^{-1}$| without significant mass growth.
For higher gas number densities (|$\sim 10^6~\rm cm^{-3}$|), IMBHs would slow down and increase in mass. Although we have not performed the simulations of high-density, high-velocity model (|$n_\infty = 10^6~\rm cm^{-3}$| and |$v_\infty =100\, \rm km~s^{-1}$|), the classical Bondi–Hoyle–Lyttleton type flow is expected to emerge because the flow structure in model ‘edgeN4V100’ is similar to the Bondi–Hoyle–Lyttleton accretion flow due to the small impact of radiation, and because the radiation impact is thought to be more ineffective as the density is increased. Hence, the accretion rate can be inferred from the classical Bondi–Hoyle–Lyttleton accretion theory as |$\sim 6\times 10^{-4}\,{\rm M}_{\odot }\rm\, yr^{-1}$|, and the acceleration is evaluated to be |$-4\times 10^{-8}~\rm cm~s^{-2}$| based on the prediction of Ostriker (1999). The time-scales of mass growth and slowdown of IMBHs are expected to be |$\sim 20~\rm Myr$| and |$\sim 7 ~\rm Myr$|, respectively. This means that the relative velocity decrease faster than mass growth if the initial velocity is |$\sim 100~\rm km~s^{-1}$|. This behaviour is thought to continue even after the relative velocity drops to |${\rm a\, \, few} \times 10~\rm km~s^{-1}$|. This is because that the slowdown is expected to be more pronounced due to the decrese of the velocity, while the mass accretion rate remains the same. Indeed, in the model ‘edgeN6V20’, the deceleration rate is |$\sim 10^{-7}\rm ~cm~s^{-2}$| and the time-scale of the slowdown, |$\sim 0.4 ~\rm Myr$|, is much shorter than the high-velocity case. The accretion rate (time-scale of the mass growth), for |$v_\infty =20 \rm$| and |$100\, \rm km~s^{-1}$| is almost the same. Subsequently, a shift from Bondi–Hoyle–Lyttleton accretion to Bondi accretion might occur. Bondi accretion in an anisotropic radiation from accretion discs around BHs has been investigated by Sugimura et al. (2017), and the accretion rate has been reported to be 1 per cent of the Bondi rate. Adopting n∞ and T∞ as number density and temperature of the gas, the time-scale of the mass growth is |$16~\rm Myr$|, which is much shorter than the age of the universe. To sum up, the IMBHs floating at speeds of |$\sim 10-100\rm\, km~s^{-1}$| rapidly slow down and grow. Eventually supermassive BHs may be formed.
In this discussion, we only consider the interaction between a single BH and the surrounding gas. To more accurately investigate the evolution of IMBHs, it should be necessary to consider their interactions with stars and other BHs.
4.2 Comparison with previous works
This work demonstrates that the dust absorption significantly changes the flow structure around moving BHs. In the low-velocity cases (model ‘edgeN4V20’ and ‘edgeN6V20’), the accretion rate is several times smaller than the analytical solution of Park & Ricotti (2013), in which dust-free situation is supposed. One reason for the such a small accretion rate may be due to the radiation force acting on the dust (see Toyouchi et al. 2020). However, the presence of dust does not significantly affect the gas flow in the high-velocity case. For the model ‘edgeN4V100’, the accretion rate is almost consistent with the analytical solution of Park & Ricotti (2013) and classical Bondi–Hoyle–Lyttleton rate. This is because the radiation force does not change the flow structure, and a classical Bondi–Hoyle–Littleton type flow emerges.
It is also important to take into account the dust sublimation. The time-averaged accretion rate is not much different between our simulation ‘isoN4V20’ and the simulations by Toyouchi et al. (2020), but oscillation of the accretion rate occurs only in our simulation ‘isoN4V20’. Since the situational setting is almost the same, the cause of this difference is thought to be the treatment of dust sublimation. In addition, the result that the moving IMBH is accelerated by the gravity when the dense shock shell is formed is consistent with previous studies (Park & Bogdanović 2017; Toyouchi et al. 2020), but the acceleration obtained with model ‘isoN4V20’ is slightly smaller than that of Toyouchi et al. (2020). This is because, in our simulations, the dense shock shell is formed somewhat further away from the IMBH. The reason for this is expected to be the small sink size and the implementation of dust sublimation, but the details are to be worked out in the future.
Here, we compare this study with our previous work that investigated the Bondi–Hoyle–Littleton accretion of dusty gas in an anisotropic radiation field without solving the hydrodynamics equations. The accretion rate is 0.6 per cent of |$\dot{M}_{\rm BHL}$| in the model ‘edgeN4V20’, while it is 20 per cent to 30 per cent for the model of the previous work that have the same luminosity as that of the quiescent phase of model ‘edgeN4V20’. In the model ‘edgeN4V20’, the gas is pushed out due to thermal pressure on the ionization front and a large amount of gas passes through without being swallowed by the BH. Thus, the accretion rate is drastically reduced. This means that hydrodynamics calculations are necessary, at least when shock is formed. Otherwise, the accretion rate would be overestimated.
4.3 Future work
In this study, we assume the uniform density distribution of the ambient gas. However, the interstellar medium is inhomogeneous, which could alter the acceleration rate and acceleration of BHs. Ruffert (1997, 1999) have reported that non-uniformity leads to flows orbiting around the sink, which reduces the accretion rate. Recently, Lescaudron et al. (2022) suggested that turbulent medium decelerates BHs using the three-dimensional magneto-hydrodynamics simulations. In their simulations, the deceleration rate increases by several tens of times that of the analytical solution with a uniform density (Ostriker 1999). However, the moving BHs could be accelerated if the radiation is taken into account. Although Sugimura et al. (2018) has calculated the accretion of gas with angular momentum on to a static BH, the simulations of moving BHs in the inhomogeneous gas, considering the radiation feedback, have not yet been performed. These are important future work.
There would also be room to modify the handling of accretion discs around the BHs. The rotation axis of the disc is likely to be perpendicular to the direction of gas motion (edge-on) as described in Section 2.4. However, the gas accreting to the central region could have random angular momentum if the ambient gas is inhomogeneous. In such a case, the edge-on disc is not maintained. In addition, immediately after the BHs encounter the gas clouds, the situation is likely to be different from edge-on case. In these cases, the shape and location of the dense shock shell and ionization front could be different from the results of this work. To be realistic, the improvements mentioned above are needed.
In our simulations, the gas flowing into the sink is assumed to accrete to the BH and change the luminosity immediately. However, more precisely, the luminous intensity should increase with a delay of about the viscous time. This time lag would not be negligible in the case that the angular momentum of the accreting gas is large. Clarifying these points is an important future work.
Spectral energy distribution (SED) of BHs depend on mass accretion rates. In this study, we assume a power-law SED (see also Section 2). If the mass accretion rate is near or above the Eddington rate, the SED of the accretion disc would be multitemperature blackbody radiation (Kato, Fukue & Mineshige 2008). Additionally, X-ray is produced due to Compton scattering by the high-temperature plasma around the disc (Kawashima et al. 2012). On the other hand, if the mass accretion rate is much lower than the Eddington rate, the spectral energy widely distributes in the range between radio and gamma-rays (Narayan & Yi 1995; Manmoto et al. 1996; Yuan, Quataert & Narayan 2003). Since the ionization and heating rate depend on the SED, the different SEDs could vary in flow structure and accretion rate. Two-dimensional radiation hydrodynamics simulations on Bondi scale, in which the multitemperature blackbody radiation is taken into consideration, showed that the critical gas density, above which a transition to super-Eddington accretion occurs, slightly decreases compared to the case of power-law spectra (Takeo et al. 2018). This is because photoionization for accretion disc spectra is less efficient than that for single power-law spectra. We need to study the Bondi–Hoyle–Lyttelton accretion considering more realistic SED.
In this study, we assume the cosine function for the angular distribution of the radiation field caused by accretion discs. The disc could produce much stronger angular dependence if the mass accretion rate is above the Eddington rate (Ohsuga et al. 2005; Watarai et al. 2005). Conversely, the angular dependence of the radiation field is weak if the mass accretion rate is much lower than the Eddington rate. Also, the angular distribution varies with time as the accretion disc undergoes precession motion. The multidimensional simulations of the accretion disc around the BHs are needed to obtain more realistic models of the radiation field (Machida, Hayashi & Matsumoto 2000; Hawley & Krolik 2001; Ohsuga et al. 2009; Ohsuga & Mineshige 2011). Two-dimensional radiation hydrodynamics simulations around a static BH have been performed with parameters such as the size of the shadow region and angular distribution of the anisotropic radiation field (e.g. Sugimura et al. 2017; Takeo et al. 2018). They suggested that the accretion rate becomes much higher in the case of a large shadow size and strong anisotropy. Investigation of the dependence of shadow size and anisotropy in the case of moving BHs is left as future work.
Outflow affects the gas flow around BHs. A strong bipolar outflow has been observed around compact objects with mass accretion rates above the Eddington rate (e.g. Fabrika 2004). Also, the jet has been detected in many sources (e.g. Walker et al. 2018). The previous studies showed that outflow pushes out the ambient gas, and accretion rate is reduced to approximately |$20-30~{{\ \rm per\ cent}}$| of the Bondi–Hoyle–Lyttleton accretion rate (e.g. Li et al. 2020; Bosch-Ramon 2022). They indicated that acceleration rate is |$40-80~{{\ \rm per\ cent}}$| smaller than without outflow. For further study on the impacts of outflow, we will include these effects in future works.
5 SUMMARY AND CONCLUSIONS
In this work, we study the Bondi–Hoyle–Lyttleton accretion mechanism on to IMBHs (104 M⊙) by three-dimensional radiation hydrodynamics simulations, taking into account the anisotropic radiation field originating from the accretion disc. We consider the situation in which the gas with relatively low-metallicity (Z = 0.1 Z⊙) flows in perpendicular to the rotation axis of the accretion disc (parallel to the disc plane). We take into account the radiation force acting on the dusty gas and decrease in the absorption coefficient due to the dust sublimation. The luminosity of the accretion disc is supposed to increase with the mass accretion rate. Our major findings are summarized as follows:
If the density of dusty gas is relatively high (|$\sim 10^4~\rm cm^{-3}$|) and the relative velocity between IMBHs and the gas is low (|$\sim 20~\rm km~s^{-1}$|), time-averaged accretion rate is 0.6 per cent of the Bondi–Hoyle–Lyttleton accretion rate (|$4 \times 10^{-6}\,{\rm M}_{\odot }\rm\, yr^{-1}$|). It is found that an ionized region like two spheres glued together at the equatorial plane appears around the IMBH and the dense shock shell is formed nearby the ionization front. The radiation force around the rotation axis of the disc works to prevent gas accretion so that the gas mainly accretes through the disc equatorial plane. The gravity of the dense shock shell accelerate the IMBHs at |$\sim 10^{-8}\rm\, cm~s^{-2}$|. Note that the accretion rate periodically increases, but does not significantly affect the time-averaged accretion rate and acceleration.
In the case of high relative velocity (|$\sim 100~\rm km~s^{-1}$|), the accretion rate, |$4 \times 10^{-6}\,{\rm M}_{\odot }\rm\, yr^{-1}$|, is approximately equal to the Bondi–Hoyle–Lyttleton accretion rate. This is because the radiation force hardly prevents the gas accretion. Even around the rotation axis of the disc, the radiation force is significantly small due to the dust sublimation. The dense shock shell is not formed, and the IMBHs are decelerated due to the gas accretion.
When the gas density is extremely high (|$\sim 10^6~\rm cm^{-3}$|), the time-averaged accretion rate is 0.6 per cent of the Bondi–Hoyle–Lyttleton accretion rate, |$\sim 4\times 10^{-4}\,{\rm M}_{\odot }\rm\, yr^{-1}$|. Due to the gravity by dense downstream gas, the IMBHs are decelerated although the dense shock shell is formed at the vicinity of the ionization front. The accretion rate oscillates periodically as in the case of the density of |$\sim 10^4~\rm cm^{-3}$|. This is because when the surface density of the shell increases to a certain degree, part of the shell falls due to the gravity of the IMBHs. Periodic variations in accretion rate occur even when the radiation is isotropic. This result differs from previous work, possibly because we have taken into account dust sublimation.
Based on the present results, the IMBHs are expected to continue floating in the early galaxies (z ≳ 6) at the velocity of |$\gtrsim {\rm several} \times 10 \, {\rm km~s^{-1}}$| without significantly mass growth, if the gas density of the galaxies is |$\sim 10^4~\rm cm^{-3}$|. On the other hand, if the density of interstellar medium is extremely high, |$\sim 10^6~\rm cm^{-3}$|, even the IMBHs with the initial velocity of about |$\sim 100 \, {\rm km~s^{-1}}$| will slow down due to momentum transport caused by the mass accretion. If the Bondi accretion begins with decreasing the velocity, the IMBHs could grow rapidly. This may lead to the evolution from IMBHs to SMBHs.
ACKNOWLEDGEMENTS
We thank an anonymous referee for fruitful comments. This work was supported by JSPS KAKENHI Grant Numbers JP22KJ0435 (EO), JP21H04488, JP18K03710 (KO), JP20H04724, JP21H04489 (HY), and JP23K13139 (HF), and Japan Society for the Promotion of Science and JST FOREST Program, Grant Number JPMJFR202Z (HY). A part of this research was funded by the MEXT as ‘Program for Promoting Researches on the Supercomputer Fugaku’ (Toward an unified view of the universe: from large-scale structure to planets, Grant Number JPMXP1020200109) (KO, HF, HY), by MEXT as ‘Program for Promoting Researches on the Supercomputer Fugaku’ (Structure and Evolution of the Universe Unravelled by Fusion of Simulation and AI; Grant Number JPMXP1020230406) and used computational resources of supercomputer Fugaku provided by the RIKEN Center for Computational Science (Project ID: hp230204, KO), and by Joint Institute for Computational Fundamental Science (JICFuS, KO). Numerical computations were performed with computational resources provided by the Multidisciplinary Cooperative Research Program in the Center for Computational Sciences, University of Tsukuba, Oakforest-PACS operated by the Joint Center for Advanced High-Performance Computing (JCAHPC), Cray XC 50 at the Center for Computational Astrophysics (CfCA) of the National Astronomical Observatory of Japan (NAOJ), the FUJITSU Supercomputer PRIMEHPC FX1000 and FUJITSU Server PRIMERGY GX2570 (Wisteria/BDEC-01) at the Information Technology Center, The University of Tokyo.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author.
References
APPENDIX A: OSCILLATION OF ACCRETION RATE FOR ISOTROPIC MODEL
Here, we explain why the second of the two peaks appearing in the burst has a larger accretion rate (luminosity). Fig. A1 represents the number density (nH), ratio of radial component of total radiation force to gravity (fR/|g|), ratio of radial component of radiation force due to EUV to gravity (|$f^R_{\rm EUV}/|g|$|), and ratio of radial component of radiation force due to FUV to gravity (|$f^R_{\rm FUV}/|g|$|) as a function of the distance from the IMBH (R). The light-pink and light-blue means the results at phase 4–2 (|$t\sim 0.21715~\rm Myr$|) and phase 4–2* (|$t\sim 0.2173~\rm Myr$|), respectively (see Fig. 4). Phase 4–2 roughly corresponds to the first peak, while phase 4–2* is a little before the second peak. The accretion rate in phase 4–2* is approximately equal to the accretion rate in Phase 4–2.
It is found that the density is larger in the phase 4–2* than in phase 4–2 at slightly outside the dust sublimation region. This difference in density produces a difference in radiation flux (radiation force) through attenuation by absorption. Indeed, we find that the radiation force in the phase 4–2* is slightly weaker than that in the phase 4–2. In phase 4–2, the accretion rate starts to decrease due to radiation force, but in phase 4–2*, the accretion rate continues to increase further due to the weaker radiation force. Thus, the accretion rate (luminosity) is larger in second peak. Note that the main component of radiation force is due to FUV. The gas density is high in phase 4–2* because relatively large amount of gas that is not absorbed in the first peak remains around the BH. Conversely, the gas density in phase 4–2 is low because a large amount of gas is swallowed in the second peak.
