Three-dimensional radiation hydrodynamics simulations of wandering intermediate-mass black holes considering the anisotropic radiation and dust sublimation

By performing three-dimensional radiation hydrodynamics simulations, we study Bondi-Hoyle-Lyttleton accretion onto 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 (~20 km/s) and gas density is ~10^4/cm^3, the gas mainly accretes from near the equatorial plane of the accretion disk at a time-averaged rate of 0.6% 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 ~10^-8cm/s^2 due to the gravity of the shell. For denser gas (~10^6/cm^3), the time-averaged accretion rate is also 0.6% of the Bondi-Hoyle-Lyttleton rate.However, the BH is decelerated at ~10^-7cm/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>several 10km/s without increasing mass in interstellar gas with density of ~10^4/cm^3, and slow down and grow into supermassive BHs in galaxies with the density of ~10^6/cm^3.


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 et al. 2010;Janka 2012;Wongwathanarat et al. 2013;Tauris et al. 2017;Kapil et al. 2023), to stars and BHs in globular clusters (Zwart & McMillan 2002;Gurkan et al. 2004;Kaaz et al. 2019;Shi et al. 2021;Gonzalez et al. 2021), to compact objects in X-ray binaries via stellar wind accretion, and to intermediate-mass BHs (IMBHs) floating in the remnant galactic disk (Matteo et al. 2005(Matteo et al. , 2012(Matteo et al. , 2017;;Dubois et al. 2012;Sĳacki et al. 2015;Latif et al. 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 et al. 1994Ruffert et al. , 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 et al. 2008;Park & Ricotti 2012; ★ E-mail: ogata@ccs.tsukuba.ac.jp 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 disk, however, the accretion rates can not decrease significantly because the gas accretes through the region near the disk plane, where the radiation effect is relatively weak.Although Bondi-Hoyle-Lyttleton accretion onto the compact objects under the anisotropic radiation field produced by the accretion disk has been investigated in previous studies (Fukue & Ioroi 1999;Hanamoto et al. 2001;Ogata et al. 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 onto 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 hydrody-namics 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  = 0.1 ⊙ as our target objects.The rest of this paper is organized as follows.We describe the method of our threedimensional 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.

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 et al. 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.

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  ∞ , positive direction velocity  ∞ , and temperature  ∞ = 180 K at the upstream boundary.We set the size of a calculation box as  out = 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 disk and set the sink radius as  sink = 2.7 × 10 −3 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  max = 10, and the minimum cell size is Δ 10 = 3.85 × 10 −4 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  max = 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 disk 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).

Basic Equations
We solve the following governing equations: the equation of continuity, the equation of motion, and the equation of energy, where  is total energy defined as Here, , , and  are the gas density, gas velocity, and gas pressure, respectively, and  and  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 eV -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 H, H 2 , H − , H + , H + 2 , e, CO, C II , O I , O II , and O III .The number density of the -th species is calculated as where   =   / H is the number ratio of -th species to hydrogen nuclei, and   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 × / ⊙ = 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 K, above which we ignore the contributions of dust grains on radiation, force attenuation, and thermal processes of gas components.

Radiative Source
In our simulations, the BH luminosity  is estimated with the mass accretion rate .We adapt the model of Watarai et al. (2000) as  where  E (≡ 4  BH / es ) and  E (≡  E / 2 ) are the Eddington luminosity and Eddington accretion rate.Here,  BH ,  es = 0.4 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   ∝  −  in the frequency range 11.2 eV < ℎ < 1 keV, where ℎ 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 et al. 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  sink  is the photon emissivity of the BH, and  and  are the cell numbers.The anisotopy factor F  (Θ, ) represents the direction dependence of the radiation as Here, the angle from the rotation axis of the disk to the center of each cell (Θ) is defined as = (cos , sin  cos , sin  sin ), (10) where ) is the distance of the origin,  is the angle between -axis and the disk rotation axis projected onto -plane, and  is the angle between -axis and the disk rotation axis, respectively.We set  = 0 in the present simulations.We also solve the diffuse IR photons produced as the dust thermal emission.

Cases examined
In our simulations, we fix the initial BH mass and the metallicity as  BH = 10 4  ⊙ and  = 0.1  ⊙ .The simulation parameters are summarized in Table 1.Hereafter, we label each model according to the anisotropy factor (F  ), the number density ( ∞ ), and gas velocity ( ∞ ) considered in each simulation.For example, 'isoN4V20' represents the simulation model with the isotropic radiation field,  ∞ = 10 4 cm −3 , and  ∞ = 20 km s −1 .
In the present study,  = /2 (edge-on) is mainly adopted.The gas inflow direction (− direction) at the upstream boundary is perpendicular to the rotation axis of the disk 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 -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 -axis.Thus, the gas falling into the central region would form the accretion disk with the rotation axis perpendicular to the -axis.Even if the rotation axis of the accretion disk is not perpendicular to the -axis at the beginning, it should eventually become perpendicular to the -axis since the gas in the initial disk is swallowed by the BHs, and the supplied gas becomes main component of the disk.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 disk, 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  ∞ = 10 4 cm −3 and velocity  ∞ = 20 km s −1 as the fiducial model.We additionally investigate the cases with high-velocity  ∞ = 100 km s −1 and highdensity  ∞ = 10 6 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 been reported the presence of relatively dense gas, ≳ 10 4 cm −3 .The high-velocity situation, The gray and cyan lines show the models with isotropic radiation ('isoN4V20') and the pole-on radiation ('poleN4V20').The magenta, green, and orange lines represent the edge-on fiducial ('edgeN4V20'), high-velocity ('edgeN4V100'), and high-density models ('edgeN6V20').The black dotted lines represent the Eddington accretion rate for dusty gas  E,dg .

Overview
Before examining the numerical results in more detail, we overview the flow structure and the accretion rates to BHs. Figure 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 (− 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 ( D <  ∞ <  R ), as also shown in Park & Ricotti (2013).Inside the shocked region, the gas slows to  D at the ionization front.If the velocity  ∞ is higher than the critical value for the R-type ionization front  R , 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.
Figure 2 shows the time evolution of the mass accretion rate.The black dashed line means the Eddington accretion rate for dusty gas  E,dg (≡  es  E /( es +  d )), where the dust opacity for UV light  d is given as  d = 2.8 × 10 2 (/ ⊙ ) cm 2 g −1 (e.g.Yajima et al. 2017).The time-averaged accretion rates  and luminosity  obtained in our simulations are summarized in Table 1.Here, we assume that the luminosity  is the function of the accretion rate  (see Eq.6).In the isotropic radiation model 'isoN4V20', the accretion rate oscillates periodically between 4.5 × 10 −2  E and 1.8 × 10 −3  E with a period of 5.5 × 10 −4 Myr.The time-averaged accretion rate is about 2 × 10 −2  E (4.3 × 10 −6  ⊙ yr −1 ).The accretion bursts (rapid increase in accretion rate) occur in the model 'edgeN6V20' with a period of ∼ 0.01 Myr.The time-averaged accretion rate is about 2  E (4.2 × 10 −4  ⊙ 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 × 10 −2  E (4 × 10 −6  ⊙ 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'.

Isotropic case
Figure 3 shows a flow structure of model 'isoN4V20' around the BH.We find that an ionized region ( H+ / H ∼ 1) extends up to  ∼ 5 pc, and a dense shock shell is formed at the upstream ionization front ( ∼ −0.5pc).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 × 10 −2  E , is three orders of magnitude smaller than the classical Bondi-Hoyle-Lyttleton accretion rate  BHL .We will show the comparison with previous work Toyouchi et al. (2020) that does not include the dust sublimation in Section 4.2.
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 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  E,dg , the gravity overcomes the radiation force in the entire region (Fig. 4-1).The gas moves towards the BH, and the When the accretion rate exceeds  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  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.

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).Figure 5 presents the flow structure in the meridian ( = 0) and equatorial ( = 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. 5-a).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 -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 disk.As a result, the time-averaged accretion rate in the quiescent phase is much smaller than Bondi-Hoyle-Lyttleton rate (∼ 0.3% of  BHL ), and nearly comparable to  dg .
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 -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 -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 ∼  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, ∼ 0.15Myr, 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 ( shell ) exceeds  ∞  2 shell  2 ∞ /  BH because the gravity,   BH  shell / 2 shell , becomes larger than the ram pressure (gas pressure),  ∞  2 ∞ .Since the distance of the dense shock shell from the BH,  shell , is about 0.5pc, the critical column number density obtained from the above relation is ∼ 10 23 cm −2 , which is roughly consistent with that the burst in the present simulation.The timescale on which the column number density of the dense shock shell reaches 10 23 cm −2 is approximately This timescale is evaluated as ∼ 1.1Myr and is longer than the free-fall time so the burst interval is determined by this timescale.

Higher velocity model
Figure 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  R .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  HII / ∞ ∼  2 ∞ 1 − √︃ 1 − 4 2 s,HII / 2 ∞ /2 2 s,HII under the condition  s,HII ≫  ∞ , where  HII ,  s,HII , and  ∞ are the gas density of ionized gas and sound speeds of ionized and neutral gas (Spitzer et al. 1978).If  ∞ ≫  s,HII , then we obtain  HII ∼  ∞ .We confirm in our simulations that the condition of  ∞ ≫  s,HII is satisfied.The somewhat low-density region appears over several pc downstream ( > 0,  ≲ 1 pc, dark-purple in Fig. 6-a).This is because the gas is accelerated by radiation force and gas pressure gradient force in the ionized region.
The radiation force does not affect the the motion of gas at distances within the Bondi-Hoyle-Lyttleton radius from the -axis so that the flow structure is similar to the classical Bondi Hoyle-Littleton accretion.Figure 6-(b) shows the flow structure in the 10 −2 pc scale.Around the disk 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 -axis.Thus, a quasi-steady flow similar to Bondi Hoyle-Littleton accretion appears.The accretion rate is roughly comparable to  BHL (see Table 1).

Dense model
Model 'edgeN6V20' differs from the fiducial model in that it produces more intense, shorter-interval accretion bursts.Figure 7 shows the time evolution of the gas density distribution in the meridian plane ( = 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. 7-a).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. 7-b).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, ∼ 0.01Myr, is much shorter than that of the fiducial model, ∼ 0.15Myr.This is because the  shell 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  shell ).
At the burst phase, the radius of the ionized region,  HII , is comparable to or slightly smaller than the Bondi-Hoyle-Lyttleton radius,  BHL , at around the equatorial plane where gas mainly accretes.On the other hand, we find  HII >  BHL in regions other than the equatorial plane.This means that the condition for super-Eddington accretion " HII <  BHL " (Inayoshi et al. 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.

Pole-on case
Figures 8-(a) and -(b) show the flow structure on a 1 pc and 10 −2 pc scales of the pole-on radiation model ('poleN4V20').The flow structure is quasi-steady and axisymmetric with respect to the -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. 8-(a)).As shown in Fig. 8-(b), this model also shows that gas accretion occurs mainly from near the equatorial plane ( = 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 -axis at the downstream region and then flows out in the + 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 − direction.

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. Figure 9 shows the acceleration in the -direction   () at the elapsed time  = 8.5 × 10 −2 Myr, which is calculated as where Ω and  sink are the solid angle and area vector.Here, a negative (positive) value of   () means that the BH accelerates in the upstream (downstream) direction and the velocity of BH relative to the interstellar gas increases (decreases).
It is found that the acceleration   ( out ) is about −10 −8 cm s −2 in the models with  ∞ = 20 km s −1 and  ∞ = 10 4 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  ∞ = 100 km/s ('edgeN4V100'), we find   ( out ) ∼ 10 −10 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 (  > 0) onto 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   ( out ) ∼ 10 −7 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, ∼ 10 −7 cm s −2 .
In the following, we discuss the evolution of the mass and velocity of the IMBHs ( BH ∼ 10 4  ⊙ ) floating in early galaxies ( ≳ 6) using our results of model 'edgeN4V20', 'edgeN4V100', and 'edgeN6V20'.Here, we use ( out ) at  ∼ 0.1Myr presented in Fig. 9 instead of the time-averaged value since the acceleration at larger distance ( ≳several pc) does not change with time as much.
In the case that the gas density is ∼ 10 4 cm −3 , the IMBHs continue to float without significant mass growth.For the model 'ed-geN4V100', the acceleration is ∼ 10 −10 cm s −2 so that the timescale of the slowdown is ∼ 3 Gyr.This is much longer than the age of the universe at  = 6 (∼ 900 Myr).This implies that the speed of IMBHs does not change so much.On the other hand, the acceleration timescale of model 'edgeN4V20', ∼ 6 Myr, is much shorter than 900 Myr.Thus, if the initial velocity is relatively small, ∼ 10 km s −1 , the IMBHs will continue to speed up until the velocity reaches several × 10 km s −1 where the acceleration timescale 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 timescale of mass growth derived from the accretion rate (∼ 4 × 10 −6  ⊙ yr −1 ) is ∼ 2 Gyr, which is much longer than 900 Myr.Thus, we conclude that the IMBHs continue to float at the velocity of ≳ several × 10 km s −1 without significant mass growth.
For higher gas number densities (∼ 10 6 cm −3 ), IMBHs would slow down and increase in mass.Although we have not performed the simulations of high-density, high-velocity model ( ∞ = 10 6 cm −3 and  ∞ = 100km s −1 ), the classical Bondi-Hoyle-Lyttleton type flow is expected to emerge because the flow structure in model 'ed-geN4V100' 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 ∼ 6 × 10 −4  ⊙ yr −1 , and the acceleration is evaluated to be −4 × 10 −8 cm s −2 based on the prediction of Ostriker (1999).The timescales of mass growth and slowdown of IMBHs are expected to be ∼ 20 Myr and ∼ 7 Myr, respectively.This means that the relative velocity decrease faster than mass growth if the initial velocity is ∼ 100 km s −1 .This behavior is thought to continue even after the relative velocity drops to a few × 10 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 ∼ 10 −7 cm s −2 and the timescale of the slowdown, ∼ 0.4 Myr, is much shorter than the high-velocity case.The accretion rate (timescale of the mass growth), for  ∞ = 20km s −1 and 100km 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 disks around BHs has been investigated by Sugimura et al. (2017), and the accretion rate has been reported to be 1% of the Bondi rate.Adopting  ∞ and  ∞ as number density and temperature of the gas, the timescale of the mass growth is 16 Myr, which is much shorter than the age of the universe.To sum up, the IMBHs floating at speeds of ∼ 10 − 100km s −1 rapidly slow down and grow.Eventually supermassive BHs may be formed.
In the present 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.

Comparison with previous works
Our present 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 the present 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% of  BHL in the model 'edgeN4V20', while it is 20 − 30% for the model of the previous work that have the same luminosity as that of the quiescent phase of model 'ed-geN4V20'.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.

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 et al. (1997Ruffert et al. ( , 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 onto 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 disks around the BHs.The rotation axis of the disk 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 disk 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 the present 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 disk would be multi-temperature blackbody radiation (Kato et al. 2008).Additionally, X-ray is produced due to Compton scattering by the high-temperature plasma around the disk (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 et al. 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 multi-temperature 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 et al. 2018).This is because photoionization for accretion disk 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 disks.The disk could produce much stronger angular dependence if the mass accretion rate is above the Eddington rate (Watarai et al. 2005;Ohsuga 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 disk undergoes precession motion.The multidimensional simulations of the accretion disk around the BHs are needed to obtain more realistic models of the radiation field (Machida et al. 2000;Hawley & Krolik 2001;Ohsuga et al. 2009;Ohsuga & Mineshige 2011).Twodimensional 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 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 et al. 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% of the Bondi-Hoyle-Lyttleton accretion rate (e.g.Li et al. 2020;Bosch-Ramon 2022).They indicated that acceleration rate is 40 − 80% smaller than without outflow.For further study on the impacts of outflow, we will include these effects in future works.

SUMMARY AND CONCLUSIONS
In the present work, we study the Bondi-Hoyle-Lyttleton accretion mechanism onto IMBHs (10 4  ⊙ ) by three-dimensional radiation hydrodynamics simulations, taking into account the anisotropic radiation field originating from the accretion disk.We consider the situation in which the gas with relatively low metallicity ( = 0.1 ⊙ ) flows in perpendicular to the rotation axis of the accretion disk (parallel to the disk 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 disk 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 (∼ 10 4 cm −3 ) and the relative velocity between IMBHs and the gas is low (∼ 20 km s −1 ), time-averaged accretion rate is 0.6 % of the Bondi-Hoyle-Lyttleton accretion rate (4 × 10 −6  ⊙ 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 disk works to prevent gas accretion so that the gas mainly accretes through the disk equatorial plane.The gravity of the dense shock shell accelerate the IMBHs at ∼ 10 −8 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 (∼ 100 km s −1 ), the accretion rate, 4 × 10 −6  ⊙ 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 disk, 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 (∼ 10 6 cm −3 ), the timeaveraged accretion rate is 0.6 % of the Bondi-Hoyle-Lyttleton accretion rate, ∼ 4 × 10 −4  ⊙ 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 ∼ 10 4 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 ( ≳ 6) at the velocity of ≳ several × 10 km s −1 without significantly mass growth, if the gas density of the galaxies is ∼ 10 4 cm −3 .On the other hand, if the density of interstellar medium is extremely high, ∼ 10 6 cm −3 , even the IMBHs with the initial velocity of about ∼ 100 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.

Figure 1 .
Figure 1.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 (light gray).The black-filled circle represents the sink region.The solid and dotted blue lines show streamlines that pass through or accrete to BH.The yellow and gray regions present ionized and dust sublimation regions, respectively.The shocked region at the ionization front is painted dark yellow.The shapes of ionized and dust sublimation regions are different in the models with isotropic and pole-on radiation fields.

Figure 2 .
Figure2.Time evolution of mass accretion rates.The vertical axis is normalized by the Eddington accretion rate  E estimated with the opacity of electron scattering.The gray and cyan lines show the models with isotropic radiation ('isoN4V20') and the pole-on radiation ('poleN4V20').The magenta, green, and orange lines represent the edge-on fiducial ('edgeN4V20'), high-velocity ('edgeN4V100'), and high-density models ('edgeN6V20').The black dotted lines represent the Eddington accretion rate for dusty gas  E,dg .

Figure 3 .
Figure 3. Gas flow structure on a 1 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 + direction.

Figure 4 .
Figure 4. 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 black solid lines and filled gray 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 black filled circles mark the epochs for each snapshot 1-4.

Figure 5 .
Figure 5. Same as Fig. 3 but for the edge-on fiducial model ('edgeN4V20').Each panel shows the flow structure on (a)1 pc, (b)10 −1 and (c)10 −2 pc scales at  = 0.2 Myr (0.04 Myr before the next accretion burst).Each column represents the flow structure in the meridian (xy-) and equatorial (xz-) planes from top to bottom.

Figure 6 .
Figure 6.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 black solid lines in panel (b) represent the dust sublimation regions.

Figure 7 .
Figure 7.The flow structure on 10 −1 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.

Figure 8 .
Figure 8. Same as Fig. 5 but for the pole-on model ('poleN4V20').Each panel shows the flow structure on (a) 1 pc and (b) 10 −2 pc scales.We only show the snapshot in the xy-plane due to the axisymmetry of flow.

Figure 9 .
Figure 9. Acceleration in the -direction of BHs at  ∼ 8.5 × 10 −2 Myr as a function of .The solid (dotted) lines mean that the acceleration is negative (positive) and the speed of the BHs increases (decreases).The line colors are the same as in Fig. 2. The star markers correspond to the positions of the ionization fronts.

Table 1 .
Model parameters with the results of the time-averaged accretion rate and the luminosity.