Accelerated growth of seed black holes by dust in the early universe

We explore the effect of dust on the growth of seed black holes (BHs) in the early universe. Previous 1D radiation-hydrodynamic (RHD) simulations show that increased radiation pressure on dust further suppresses the accretion rate than the case for the chemically pristine gas. Using the Enzo+Moray code, we perform a suite of 3D RHD simulations of accreting BHs in a dusty interstellar medium (ISM). We use the modified Grackle cooling library to consider dust physics in its non-equilibrium chemistry. The BH goes through an early evolutionary phase, where ionizing BH radiation creates an oscillating HII region as it cycles between accretion and feedback. As the simulations proceed, dense cold gas accumulates outside the ionized region where inflow from the neutral medium meets the outflow driven by radiation pressure. In the late phase, high-density gas streams develop and break the quasi-spherical symmetry of the ionized region, rapidly boosting the accretion rate. The late phase is characterized by the coexistence of strong ionized outflows and fueling high-density gas inflows. The mean accretion rate increases with metallicity reaching a peak at Z$\sim$0.01-0.1$\,Z_\odot$, one order of magnitude higher than the one for pristine gas. However, as the metallicity approaches the solar abundance, the mean accretion rate drops as the radiation pressure becomes strong enough to drive out the high-density gas. Our results indicate that a dusty metal-poor ISM can accelerate the growth rate of BHs in the early universe, however, can stun its growth as the ISM is further enriched toward the solar abundance.


INTRODUCTION
Bright quasars powered by supermassive black holes (SMBHs) shine at the edge of the observable universe at redshifts greater than six. Only after less than a billion years, these distant beacons already have masses above 10 9 M (e.g., Fan et al. 2001;Willott et al. 2003;Mortlock et al. 2011;Wu et al. 2015;Bañados et al. 2018), comparable to the most massive BHs in the local universe. Their rapid growth during the early universe raises fundamental astrophysical questions about their formation, evolution, and connection with galaxy formation (Kormendy & Ho 2013).
Numerical simulations suggest that intermediate-mass black holes (IMBHs) in the range of M BH = 10 2 -10 5 M tive feedback, where UV and X-ray photons heat and ionize the surrounding interstellar medium (ISM), creating an H ii region around an accreting BH (Milosavljević et al. 2009;Park & Ricotti 2011. This hot bubble filled with low-density gas regulates the gas supply to the BHs. Due to the cycles between accretion and feedback, the accretion behavior is highly oscillatory creating constantly fluctuating H ii regions. That introduces turbulence in the ISM, however Park, Wise, & Bogdanović (2017, hereafter, PWB17) find that the the thermal energy dominates the turbulent kinetic energy.
This feedback-dominated accretion occurs when average size of the H ii region R s exceeds the Bondi radius r B = GM BH /c 2 s . Here, M BH is the BH mass, and c s is the sound speed of the gas. By definition, the Bondi radius is the scale within which the gravitational energy of the BH dominates over the thermal energy of the gas, and R s is another scale that the BH radiation heats and ionizes. Note that this scale is by far larger than the accretion disk scale and regulates gas supply to the BH from galactic scales. The condition for feedback-limited regime assuming the gas temperature T ∞ = 10 4 K and the standard radiative efficiency η = 0.1, can be simplified as the product of the BH mass and the gas density as M BH n H,∞ ≤ 10 9 M cm −3 where n H,∞ is the hydrogen number density of the gas. The mean accretion rate can be expressed assuming the spectral index is α spec = 1.5 as the following (Park & Ricotti 2012) Ṁ BH ∼ 6 × 10 −6 M 2 BH,4 n H,∞,3 T ∞,4 M yr −1 , (1) where M BH,4 = M BH /(10 4 M ), n H,∞,3 = n H,∞ /(10 3 cm −3 ), and T ∞,4 = T ∞ /(10 4 K). Note that the mean accretion rate is proportional to the gas pressure for a given BH mass (i.e., Ṁ BH ∝ n H,∞ T ∞ ).
In more extreme conditions, the accretion flow makes a transition to feeding-dominated regime (Pacucci et al. 2015).
When R s ≤ r B which is translated as M BH n H,∞ ≥ 10 9 M cm −3 assuming the same T ∞ and η, the ionization front becomes unstable and is accepted as the condition for hyperaccretion (Park et al. 2014;Inayoshi et al. 2016). There have been several works to study the geometry of the hyperaccretion flow such as anisotropic radiation  or biconicaldominated accretion flow (BDAF) and decretion disk (Park et al. 2020).
Several scenarios have been suggested for exponential accretion (e.g., Alexander & Natarajan 2014). Among them, it might be worth noting that some works pointed out the possible connection with the early assembly of stellar bulge (Park et al. 2016;Inayoshi et al. 2021). Park et al. (2016) pioneered the possibility of bulge-driven growth of seed BHs. This work predicts that there exists a critical bulge mass of ∼ 10 6 M above which the rapid growth of seed BHs is triggered no matter of the central BH masses. This mechanism could favor DCBHs as the seeds of high-redshift quasars because immediately their formation, efficient star formation is triggered by positive radiative feedback (Barrow et al. 2018). This newly formed nuclear star cluster would then aid in subsequent growth as the young stellar bulge component increases the effective accretion radius supplying more gas from large scale to the seed BHs. Recently, Inayoshi et al. (2021) perform asymmetric 2D simulations, suggesting obese BHs (OBGs; Agarwal et al. 2013) as the progenitor of the high-redshift quasars, which are significantly overmassive compared to the observed bulge-to-BH mass relation (Kormendy & Ho 2013).
Another critical uncertainty in the growth of seed BHs, which is the main focus of the current study, is the role of dust. Observation discovers that the early universe is rapidly contaminated with dust by supernovae. For example, a star-forming galaxy A1689-zD1 at z > 7 that is heavily enriched in dust, and the dustto-gas ratio is close to that of the Milky Way (Watson et al. 2015). Atacama Large Millimeter Array (ALMA) observations also detected dust-contaminated quasars at 6.0 z < 6.7 with the estimated dust masses of M d = 10 7 -10 9 M . Cosmological simulations also predict that the early universe is rapidly contaminated by the ongoing star-formation. A single pair-instability supernova can enrich a halo of 10 7 M to a metallicity of 10 −3 Z , triggering a transition from Pop III to Pop II star formation (Bromm et al. 2003;Wise et al. 2012). Yajima et al. (2015) show that the amount of dust in an overdense region of massive galaxies can reach the solar neighborhood level even at z 6. Interestingly, the contamination of the dust is closely associated with the bulge-driven growth of the BHs discussed above since the assembly of the stellar bulge leads to metal enrichment of the ISM around the BHs. Recently, Ishibashi (2021) explores the parameter space for dust photon trapping (DPT) which might lead to the supercritical BH growth.
Previous work by Yajima et al. (2017) was performed using 1D RHD simulations for a relatively short period of time compared to the time scale in the current study. They discover that the radiation pressure on dust is the important parameter compared to the dust attenuation effect. The average accretion rate is lower than the case of the primordial gas. The spectral energy distribution (SED) is calculated from the dust thermal emission and they suggest that the flux ratio between λ 20µm and 100 µm shows a close relation to the Eddington ratio. Future multi-messenger observations of the cosmic dawn will reveal how the first luminous structures coevolve with massive BHs. Some examples include radiation from the first galaxies and BHs with JWST (James Webb Space Telescope) (e.g., Natarajan et al. 2017), X-ray observations of actively accreting BHs by ATHENA (Advanced Telescope for High-ENergy Astrophysics), and gravitational waves from merging BHs by LISA (Laser Interferometer Space Antenna).
This paper builds upon the previous 1D work by Yajima et al. (2017) to investigate the effect of the dust contamination onto BHs growth in the early universe using 3D RHD simulations with various metallicities. We emphasize that we investigate the case of the feedbackdominated regime with dust physics included, where the growth of the BH is highly suppressed for chemically pristine gas. We find that the dust cooling lowers the mean accretion rate due to lower thermal pressure in the early phase, however, boosts the accretion rates when the spherical symmetry is broken and the dense stream of gas develops an inflow channel. We also find that metallicity affects the mean accretion rates. In Section 2, we describe the numerical simulations. In Section 3 we present the simulation results. We discuss the results in Section 4 and summarize in Section 5.

Radiation hydrodynamic simulations
We perform 3D RHD simulations to study the role of dusty gas on the growth seed BHs under the influence of radiative feedback. We use the adaptive mesh refinement (AMR) code Enzo equipped with the Moray package to solve radiative transfer equations (Wise & Abel 2011;Bryan et al. 2014;Brummel-Smith et al. 2019).
We use the chemistry and cooling library Grackle O + 2 , O 2 , Mg, Si, SiO and SiO 2 . In this study, we also trace the evolution of 2 grain species abundances: silicate (enstatite, MgSiO 3 ) and graphite (C). We calculate the cooling by dust thermal emission, H 2 formation on grain surfaces, and continuum opacity of the silicate and graphite. We describe the detail of the dust physics in Section 2.4. Table 1 shows the list of simulations. We fix the BH mass M BH = 10 4 M and gas density n H,∞ = 10 3 cm −3 . We also use the same domain size of (80 pc) 3 with the resolution of ∆D min = 0.625 pc. The IDs in the first column of Table 1 show the metallicities as Zmn where n is Z/Z = 10 −n . The 3rd column is the time τ tran from an early oscillatory phase to a late stable phase with higher accretion rates, and the 4th and 5th columns are the mean accretion rates for each phase, respectively.

Accretion rates of dusty gas
In the simulation, we follow the accretion flows in a dusty ISM around a seed BH. To quantify their behavior and evolution, we compare our results to the standard Bondi solution. The thermal state of the accretion flow is affected by adiabatic gas heating, shocks, photoionization, recombination, and dust heating and sublimation. To calculate the bolometric accretion luminosity, we use the Eddington-limited Bondi recipe that sources the radiative transfer equation. The Bondi radius is normalized as where c ∞ is the sound speed of the gas outside the ionized region that is not affected by the BH radiation or gravity and T ∞ is the relevant gas temperature. We normalize the gas temperature of T ∞ = 10 3 K due to dust cooling and adopt c ∞ = 9.1 km s −1 (T ∞ /10 4 K) 1/2 assuming isothermal gas with a mean molecular weight of 1. The Bondi accretion rate is defined aṡ where the dimensionless accretion rate λ B depends on the equation of state and ρ ∞ is the gas density far from the BH. Since the gas is under the influence of the BH radiation, we use the gas density ρ HII and sound speed c HII of the cell that contains the BH, arriving aṫ where r B,HII = GM BH /c 2 HII is the local Bondi radius under the influence of radiation. By using the ratio between the temperatures of the ionized and neutral region Note-MBH = 10 4 M and nH,∞ = 10 3 cm −3 , D box =80 pc, and ∆Dmin=0.625 pc for all runs. We assume Z = 0.02. The upper and lower errors for each accretion rate represent the 5% and 95% (2σ) percentiles.
where the pressure equilibrium across the ionization front is applied and the ionization fraction in the H ii region is ∼ 1 as 2ρ HII T HII = ρ ∞ T ∞ . Then the BH accretion rate is whereṀ BH,d is 1.8 × 10 −6 times of the Bondi rateṀ B when T HII = 6 × 10 4 K and T ∞ ∼ 3 × 10 2 K due to dust cooling. For comparison, the mean accretion rate for the case of pristine gas is ∼ 1% of the Bondi rate because the temperature ratio is δT ∼ 6 with T ∞ ∼ 10 4 K (Park & Ricotti 2012).

Size of the ionized region
As the BH accretes and shines, an H ii region forms, and in this paper we track how it evolves. Below we describe the connection between the BH accretion luminosity and the size of the resulting H ii region. The BH luminosity is calculated from the accretion rate as where we adopt the radiative efficiency of η = 0.1 assuming a thin disk model (Shakura & Sunyaev 1973) and c is the speed of light. We simply adopt the Eddington luminosity as the cap for the luminosity since the accretion rates mostly stay below L Edd . However, it is worth noting that global 3D magneto-hydrodynamic simulations show that the accretion luminosity converges to in the order of 10 L Edd even whenṀ BH 10 3 L Edd /c 2 (Sakurai et al. 2016;Jiang et al. 2019). The Eddington luminosity is where m p is the proton mass and σ T is the Thomson cross-section for electrons. Therefore the mean size of an ionized sphere is calculated as where N ion is the total number of ionizing photons in the range of 13.6 eV ≤ E ≤ 100 keV with a spectral index of α spec = 1.5 for power-law energy distribution, and α rec is the case B recombination coefficient. In the case of pristine gas, comparing Equations (2) and (9), the parameter space that we explore with M BH = 10 4 M and n H,∞ = 10 3 cm −3 belongs to the feedback-limited regime as M BH n H,∞ = 10 7 M cm −3 that is less than the critical threshold of 10 9 M cm −3 (Park et al. 2014). Photo-ionization, photo-heating, and gas cooling are computed by the adaptive ray-tracing module Moray that couples these rates to the hydrodynamic equations (see section 2.2 of PWB17). Here, we use 4 energy bins of (28.4, 263.0, 2435.3, 22551.1) eV, each with a fractional luminosity of (0.6793, 0.2232, 0.0734, 0.0241), respectively (see section 2.3 in PWB17 for details).

Dust modeling
For a given metallicity Z, the initial carbon and iron abundances are set as A(C) = 8.43 and A(Fe) = 7.50 (Asplund et al. 2009) for the solar metallicity, where we use Z = 0.02 (Grevesse & Sauval 1998). 2 For compact, spherical dust grains with a radius a d , we calculate the dust temperature T d from the thermal balance between radiation heating, collisional heat transfer from gas to dust and dust thermal emission as where F ν is the flux of radiation from a BH. Q ν is the absorption coefficient, and κ P,d is the Planck-mean opacity of dust per unit dust mass taken from Nozawa et al. (2008). m d = 4π 3 d a 3 d is the mass of a dust particle, where d is the bulk density of a grain. We use d = 3.20 and 2.28 g cm −3 for silicate and graphite, respectively. The gas cooling rate, H 2 formation rate and grain opacity are calculated with the derived dust temperature (in more detail, see Chiaki et al. 2015). Silicate and graphite sublimate above the dust temperatures 1222 K and 1800 K, respectively (Pollack et al. 1994). We model that nuclei tied up in grains are released into the gas phase as where "s" denotes the solid phase.
We add the radiation pressure on dust grains of uniform size in addition to the radiation forces on free electrons and neutral hydrogen. The optical depth is expressed as where m H is hydrogen mass, ρ d is the dust mass density, and D is the dust-to-gas mass ratio. We use the value D ≡ M dust /M H from each cell. D is set to 0.005(Z/Z ) for both silicate and graphite. We assume a single-sized grain model with a d = 0.1 µm. Note that observations of dust in the Milky Way indicates that the size distribution is well approximated by the power-law (Mathis et al. 1977) The radiation pressure on dust is much larger than the case for Compton scattering on electrons by a factor of f dust rad f e rad ∼ 7.1 × 10 2 a d 0.1 µm which is proportional to the metallicity Z/Z . This ratio indicates that the radiation pressure on dust dominates over the case for electron for Z 10 −3 Z . UV photons are absorbed and re-emitted in the IR by dust grains, which can then escape from the system without much absorption. In this study, we do not consider the radiation pressure by the IR, however, it should be considered when the column density is above 10 22 cm −2 .

RESULTS
All the simulations listed in the Table. 1 with different metallicities which span 6 orders of magnitude, show a similar evolution pattern. They start from a highly oscillatory quasi-spherical mode that we refer to as the early phase and make a transition to high accretion mode with broken symmetry with smaller variations that we refer to as the late phase.

Accretion rates evolution
The left panel of Fig. 1 shows the early evolution of accretion rates compared to the previous simulations with zero-metallicity (grey line) from PWB17. The light red line is for Zm6 with Z = 10 −6 Z while dark red is for Zm0 with Z = Z . Bondi accretion rate for the temperature of T ∞ = 8 × 10 3 K (dot-dashed) and Eddington rate for the M BH = 10 4 M with η = 0.1 (dashed) are shown for references. All the simulations immediately start from decreasing accretion rates and the runs with dust drop to a much lower minimum rate ofṀ BH ∼ 10 −8Ṁ yr −1 at t ∼ 0.8 Myr. The accretion rates shows a repeated pattern of fast rising and exponential decaying since then.
Previous simulations with zero-metallicity find that the mean accretion rate is approximately 2 orders of magnitude lower than the classical Bondi accretion rate for T ∞ = 10 4 K. With dust included in the study, the mean accretion rates are approximately two orders of magnitude lower than the case for zero-metallicity gas. Note that the mean accretion rate for pristine gas is ∼ 7 × 10 −6 M yr −1 , whereas with dust it is ∼ 2 × 10 −7 M yr −1 as shown in the 4th column of Table 1. Using Eq. (6) to compare the current accretion rate relative to the zero-metallicity case, we find that it scales to whereṀ BH,8000K is the accretion rate for T ∞ = 8000 K. With the temperature contrast between the ionized and neutral region δT ∼ 6 and δT d ∼ 200, the ratio becomes Ṁ BH / Ṁ BH,8000K 0.02 which is consistent with the mean accretion rates in the simulations. Fig. 1 shows that the period of bursting accretion is similar to the case of the pristine gas despite the much lower accretion rates. This similarity is explained by the pressure equilibrium across the ionization front. It is affected by the change in temperature in the following quantities in Eq. (14) that originate from the Strömgren radius Eq. (9): the electron number density n e = (δT /δT d )n H,∞ (considering a completely ionized medium), the hydrogen number density n H = (δT /δT d )n H,∞ , the recombination rate α rec ∝ T −1/2 ∞ , and the ionizing photon luminosity N ion = 0.02N ion,8000 . The last relation comes from the decreased accretion rate as discussed in the previous subsection. The mean size of the dust-affected ionized region is R s d = 0.02 1 3 δT d δT 2 3 300 K 8000 K 1 6 R s 8000 1.6, (15) which expect 60% larger in radius compared to PWB17. This size is consistent with PWB17 when compared to the average size of R s . The behavior is also consistent in terms of the average period between the cycles as Zm0 and Zm6 runs show longer periods in the left panel of Fig. 1. This is due to the fact that the period is proportional to the R s assuming a similar gas depletion rate inside the H ii region (see Park & Ricotti 2011, for details). Fig. 2 shows the density, H ii fraction, temperature, and radial velocity snapshots at t = 5.5, 6.0, 6.5, 7.0, 7.5, and 8.0 Myr. Ionizing radiation from the BH creates a low-density, highly ionized, and hot bubble, similar to the cases of PWB17. However, the temperature structures show a clear difference. At scales larger than 20 pc, the gas temperature is not affected by the radiation, and a spherical density front forms from a fluctuation in the ionized region. The gas temperature drops to T ∞ ∼ 3×10 2 K from the initial condition of T ∞ = 10 4 K. The lower gas temperature and associated lower thermal pressure decrease the mean accretion rates (Park & Ricotti 2012). The radial velocity snapshots in the bottom row clearly show the fluctuation of the ionized region due to the cycles between accretion and feedback. The fluctuation drives density waves which are seen as enhanced density shells at r ∼ 20 pc. Also, note that high-energy X-ray photons selectively propagate out of ionized region while most of the low-energy UV photons are trapped inside the low-density region. Streaks of enhanced H ii fraction are observed in the H ii fraction snapshots at 10 < r < 20 pc. This creates the intermediate temperature layer of T int ∼ 10 3 K between the gas reservoir and the ionized region.

Evolution of thermal structure
There is another thin layer of low temperature gas just outside the ionized region. Gas inflow from reservoir and outflow due to radiation feedback converge just outside the H ii region and increase the gas density. Enhanced gas density promotes the radiative gas cooling further.

Stable and enhanced accretion rates
The right side of Fig. 1 shows the early and late phase accretion rates up to t = 23.5 Myr. The early phase of the accretion rates of all the simulations with different metallicities is characterized by the oscillatory Radial Velocity (km/s) Figure 2. Early evolution for Zm0 run. From top to bottom, the panels show slices of density, H ii fraction, temperature, and radial velocity at t = 5.5, 6.0, 6.5, 7.0, 7.5, and 8.0 Myr (left to right). The early phase is characterized by the oscillatory behavior of the ionized region due to the accretion and feedback cycles.
behavior while the late phase accretion rates are increased continual accretion rates. The accretion pattern makes a noticeable change at t ∼ 7.2 Myr for Zm0 while the transition of the other runs occurs earlier around t ∼ 4.7 Myr. All of the transition times τ tran are listed in the 3rd column of Table 1, generally occurring later for higher metallicities. The oscillatory accretion rates in the early phase vary less and the overall accretion rates increase. The run Zm0 shows a longer variation period, e.g., the accretion rate isṀ BH ∼ 3 × 10 −5 M yr −1 at t ∼ 17 Myr and drops toṀ BH ∼ 10 −7 M yr −1 at t ∼ 20 Myr. Other runs also show the minor variation of accretion rates, but clearly show a higher average than Zm0. Last column of Table 1 shows the mean accretion rates Ṁ late for 10.0 ≤ t ≤ 25.0 Myr. Note that the average accretion rates are ≥ 10 −5 M yr −1 which is higher than the Ṁ early by ∼ 2 orders of magnitude except for Zm0. Here the thin dense shell of gas becomes clumpy, allowing for the X-ray photons to leak preferentially through the porous medium, which heat and ionize the affected regions.

Asymmetric thermal structure
The radial velocity slices show that inflows and outflows develop in different directions. The small outflow structure grows in size at t = 15.0 Myr. The increasing BH luminosity can not drive the gas clumps away from the central BH. In the subsequent snapshots, the strong outflow and high-density inflow grow simultaneously. The ionized region contains strong outflows that are surrounded by a cold dense shell. The morphology of the ionized region at t ≥ 20 Myr becomes highly random depending on the direction of the outflows. Fig. 4 shows the phase diagrams in the early (t = 0.22 Myr) and late phase (t = 17.0 Myr) for the gas within r ≤ 40 pc for tor Zm0 run. In the early phase, most of the gas cools to T ∼ 3 × 10 2 K, and the temperature of the gas inside the ionized region increases to T = 2 -6 × 10 4 K. In the early phase, the dashed line    shows the line for the pressure equilibrium p early ∝ ρ gas T gas ∼ 1 × 10 −18 g cm −3 K .

Evolution of thermal pressure
Pressure equilibrium is maintained between the neutral and ionized region, and the gas above the line is the dense gas located outside the ionized region shown in Fig. 2. In the late phase, the gas temperature of the neutral region increases to T ∼ 2 × 10 3 K meaning that the thermal pressure also increases by ∼ 1 order of magnitude compared to the early phase. The dot-dashed line shown in the right phase diagram is The gas below the thermal equilibrium line is the transitional outflow which is supported by the momentum in addition to the thermal pressure. Most of the gas is generally in pressure equilibrium, however, momentumdriven outflows drive more gas above this pressure equilibrium line.

Thermal structure of gas and dust
The vorticity ω = ∇ × v provides a measure of turbulence in units of s −1 that is approximately the inverse eddy turnover time. Fig. 5 shows the magnitude of vorticity | ω| with velocity streams at t = 15, 20, 25, and 30 Myr for Zm0 run. Outflows dominate the interior of the ionized region while low-velocity inflows are dominant at larger scales. Just outside of the ionization front, the outflows and inflows converge, driving turbulence and causing the outflow structure and direction are highly variable. Fig. 6 shows slices of the gas temperature, H ii fraction, and MgSiO 3 density and temperature from left to right at t = 30 Myr for Zm0 run. The thermal structure of the ionized region is highly asymmetric depending on the stage of the evolution. As discussed earlier, the gas is characterized by three different temperature structures. The ionized region has the highest temperature at T HII ∼ 6 × 10 4 K that is surrounded by cooler dense gas at T ∼ 10 3 K. The H ii fraction shows that the high-density gas is partially ionized by the leaking X-ray photons. The MgSiO 3 density, shown as a fraction of the initial gas density, has a similar structure to the gas density because we use a fixed ratio of dust-to-gas. The dust density is depleted inside the ionized region with an average dust temperature of T MgSiO3 ∼ 20 K and a minor enhancement only nearby the BH. The dust temperature increases only in the high-density region and not in the ionized region because of its low density. Dust sublimation is not captured by our simulations because we do not resolve the sublimated regions that have T d > 1,500 K. Fig. 7 shows the mean accretion rate in physical units as a function of metallicity. Blue circles show the accretion rates at early phase Ṁ early while red squares show the mean rate at late phase Ṁ late . The values of Ṁ early for each metallicity is calculated by taking the mean rate at 2.0 ≤ t ≤ 4.0 Myr while the mean at 10.0 ≤ t ≤ 25.0 Myr is taken for Ṁ late . The Ṁ early and Ṁ late are listed in the 4th and 5th columns of Table 1. The upper and lower limits for the mean rates shown as error bars represent the 5%/95% (2σ) percentiles, respectively. The small symbols show the medians for each phase and the offset from the mean value indicates that the distribution is not Gaussian.

Mean accretion rates as a function of metallicity
In the early phase, the accretion rates are about an order of magnitude lower than the case of zero metallicity, shown as a dot-dashed line. It becomes even lower with high metallicity runs with Z ≥ 10 −1 Z , consistent with Yajima et al. (2017), who find that the mean accretion rate is one order of magnitude lower than the case for primordial gas for Z ≤ 0.01Z . The standard deviations for each metallicity are shown as error bars which do not vary much as a function of metallicity.
In the late phase, the mean accretion rates Ṁ late are in general higher than the zero-metallicity case, and it is the highest for Z = 10 −2 Z by about an order of magnitude. Ṁ late increases with metallicity for low metallicities until it peaks Z ∼ 10 −2 Z . The Z ∼ 10 −1 Z case has similarly high rates, however, the mean accretion rate drops below the primordial case when the metallicity reaches the solar value Z = Z .

DISCUSSION
We have built upon the previous work of Yajima et al. (2017), finding that the early phases of our simulations are in good agreement. By evolving the simulations significantly longer, we find that the accretion transitions to an increased and less variable rate. This late phase might last until the available gas is consumed by star formation or BH accretion. With increasing mass, the BHs will become more susceptible to hyperaccretion.
Cosmological simulations with star formation and feedback, the inclusion of dust, and accretion and feedback from seed BHs will provide insightful clues toward a more complete picture of BH accretion lifecycles in the early universe. Our simulations are idealized to focus on the dust physics on the growth of seed BHs, however, the dust and metals produced from supernovae and the radiative feedback from the BHs should be considered together to understand the role of dust contamination on BH growth. In particular, we found that in the early phases of accretion, pressure equilibrium  across the ionization front regulates the gas accretion rates, and accordingly, the accretion rate estimate given in Eq. (14) matches with the simulations. However, in the late phase, the spherical symmetry is broken, and thus the accretion rates from Eq. (1) only apply to the early phase. Furthermore, our simulations are capable of tracking the sublimation of the dust for graphite and silicates, however, the sublimation of dust might occur only at the small distance from the BHs, which is not resolved in our simulations. Dust temperatures increase near the BH, but at the levels of T d ∼ 30 K. Dust sublimation near the BHs seems to have a minor impact on the accretion rate due to the low fraction of the dust compared to hydrogen and helium gas. Thus, the accretion rates from our simulations should not be affected greatly in higher resolution simulations that better captures dust sublimation.
Our simulations include physical processes over a range of temperatures and photon energies, such as a non-equilibrium chemical network with dust species and the radiation transport of X-rays. The resulting physical structures will therefore provide solid basis for observational predictions. Existing radio telescopes, such as ALMA, and space-based IR telescopes, such as JWST, may detect the signatures of the rapidly growing BHs in the early universe. For example, Yajima et al. (2017) describe an IR observational signature, the flux ratio of 20 µm and 100 µm, that is closely related to the Eddington ratio. Radiative transfer calculations (e.g. Narayanan et al. 2021) compute detailed SEDs from the stars, black holes, and gas, ranging from X-rays to the sub-mm, using the simulation physical fields as input.
Our simulated systems are highly variable and turbulent and have anisotropic features in the temperatures and abundances of different dust species and will vary to some degree depending on viewing angle. The X-ray photons that leak through the porous medium could be detectable by future X-ray missions, such as ATHENA, will also provide clues about actively accreting seed BHs.
We will address such observational signatures in a later study.

SUMMARY
In this study, we run a suite of 3D RHD simulations using Enzo+Moray with dust physics to study the role of dust in enhancing the growth of seed BHs in the early universe. We find that a dusty ambient medium, compared to the zero-metallicity case, induces changes in the thermal structures surrounding the BH, leading to a change in accretion behavior. Our simulations show distinctive phases as a function of time (i.e., the early and late phases). These two phases are different in the accretion rates and thermal structure of the neighboring gas. The following list the main discoveries.
• The early phase: Radiative cooling by dust decreases the gas temperature below T ∼ 3 × 10 2 K, and the associated drop in thermal pressure reduces the mean accretion rate by approximately two orders of magnitude compared to the case of chemically pristine gas. Nevertheless, radiation from the accreting BH ionizes and heats the surrounding gas. The system cycles between accretion and feedback modes, causing oscillatory behavior in both the accretion rates and H ii region size.
• On average, thermal pressure equilibrium across the ionization front is maintained, and the density in the H ii region is low due to the temperature contrast across the ionization front. Low densities and high recombination rates make the size of the ionized region comparable to the pristine case despite the low accretion rate.
• Outflows within the ionized region and inflows from a larger scale converge at the ionization front and enhance the gas density. Hard X-rays leak into the neutral region, heating and ionizing the gas, that result in distinctive regions with varying temperatures.
• Transition: The high-density gas outside the ionization front becomes clumpy and develops outflow and inflow channels. As the high-density gas fuels the BHs, the accretion rate and thus the BH luminosity increase rapidly.
• The late phase: The quasi-spherical symmetry is broken, and a high accretion rate and strong feedback are maintained simultaneously. The ionized region is surrounded by cold dense gas and its shape evolves randomly depending on the dominant direction of the outflows.
• The mean accretion rate in the late phase increases as a function metallicity for Z ≤ 0.1Z . However, it drops to below the level for the case pristine gas when the metallicity approaches the solar abundance.
Our numerical simulations suggest that dust contamination by star formation can play a critical role in the growth of seed BHs in the early universe. When the ISM is metal-poor (Z ≤ 0.1Z ), the growth rate is enhanced by a dust ambient medium, suggesting that the growth of seeds is closely connected with co-eval star formation. In this sense, this finding is closely related to the previous work on the role of stellar bulges in boosting accretion rates (Park et al. 2016). However, our results also suggest that when the ISM is sufficiently enriched near the solar metallicity, BH growth might be suppressed by strong feedback, implying that there is an optimal metallicity range for the rapid growth of seed black holes.
KHP and JHW are supported by the National Science Foundation grant OAC-1835213 and NASA grants NNX17AG23G and 80NSSC20K0520. GC was supported by Overseas Research Fellowships of the Japan Society for the Promotion of Science (JSPS). Numerical simulations presented were performed using the opensource Enzo and the visualization package yt (Turk et al. 2011). This work used the HIVE cluster, which is supported by the National Science Foundation under grant number 1828187. All the simulations presented here were supported in part through research cyberinfrastructure resources and services provided by the Partnership for an Advanced Computing Environment (PACE) at the Georgia Institute of Technology, Atlanta, Georgia, USA.