Dusty Gas Accretion onto Massive Black Holes and Infrared Diagnosis of the Eddington Ratio

TL;DR

This study uses radiative-hydrodynamics simulations to explore dust's role in SMBH growth, revealing that dust influences accretion rates and spectral signatures, which can be used to estimate Eddington ratios via infrared observations.

Contribution

It provides new insights into dust's secondary role in black hole growth and proposes infrared diagnostics for Eddington ratios using JWST, ALMA, and Spitzer.

Findings

Dust causes gentle accretion with small fluctuations.

Eddington ratio ranges from 10^{-4} to 10^{-2} in dusty clouds.

Infrared flux ratios correlate with accretion activity.

Abstract

Evidence for dust around supermassive black holes (SMBHs) in the early Universe is strongly suggested by recent observations. However, the accretion mechanism of SMBHs in dusty gas is not well understood yet. We investigate the growth of intermediate-mass black-holes (IMBHs) of $\sim 10^{5}~M_{\odot}$ in dusty clouds by using one-dimensional radiative-hydrodynamics simulations. We find that the accretion of dusty gas onto IMBHs proceeds gently with small fluctuations of the accretion rate, whereas that of pristine gas causes more violent periodic bursts. At dust-to-gas mass ratios similar to the solar neighborhood, the time averaged luminosity becomes smaller than that for primordial gas by one order of magnitude and the time-averaged Eddington ratio ranges from $\sim 10^{-4}$ to $\sim 10^{-2}$ in clouds with initial gas densities of $n_{\rm H} = 10 - 1000~\rm cm^{-3}$. Our calculations show that the effect of dust opacity alone is secondary compared to the radiation pressure on dust in regulating the BH growth. We also derive spectral energy distributions at IR bands by calculating dust thermal emission and show that the flux ratio between $\lambda \lesssim 20~\rm \mu m$ and $\gtrsim 100~\rm \mu m$ is closely related to the Eddington ratio. Thermal emission from hot dust near the BH dominates only during the high accretion phase, producing higher flux density at $\lesssim 20~\rm \mu m$. Therefore, we suggest that the combinations of MIR observations by JWST and FIR observation by ALMA or Spitzer can be used to estimate the Eddington ratio of massive BHs.

Tables

Table 1. Table 1 Model parameters

Model	$M_{\mathrm{BH}}\left[M_{\odot }\right]$	$n_{\infty }\left[{\mathrm{cm}}^{-3}\right]$	Dust attenuation	$f_{\mathrm{rad}}^{\mathrm{dust}}$
M5-Z0	${10}^5$	100	$\times$	$\times$
M5-Z1	${10}^5$	100	$○$	$\times$
M5-Z1rad	${10}^5$	100	$○$	$○$

Equations

$$\begin{array}{rl}& \frac{\partial \rho }{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}(r^2\rho v)=0,\\ & \rho \left(\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial r}\right)=-\frac{\partial p}{\partial r}-\frac{G{M}_{\mathrm{BH}}\rho }{r^2}+f_{\mathrm{rad}},\\ & \rho \left(\frac{\partial e}{\partial t}+v\frac{\partial e}{\partial r}\right)=-p\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^{2}v\right)+\Gamma -\Lambda ,\end{array}$$

$$\begin{array}{r}d{\tau }_{\mathrm{d},\nu }=Q_{\nu }\pi a_{\mathrm{d}}^2\frac{D{m}_{\mathrm{H}}{n}_{\mathrm{H}}}{m_{\mathrm{d}}}dl=\frac{3Q_{\nu }D{m}_{\mathrm{H}}{n}_{\mathrm{H}}}{4{\rho }_{\mathrm{d}}{a}_{\mathrm{d}}}dl,\end{array}$$

$$a_{\mathrm{d},1}=\left(\frac{3-\alpha }{4-\alpha }\right)\left(\frac{a_{\mathrm{d},\max }^{4-\alpha }-a_{\mathrm{d},\min }^{4-\alpha }}{a_{\mathrm{d},\max }^{3-\alpha }-a_{\mathrm{d},\min }^{3-\alpha }}\right),$$

$$\begin{array}{rl}{f}_{\mathrm{d}}\equiv \frac{f_{\mathrm{rad}}^{\mathrm{dust}}}{f_{\mathrm{rad}}^{\mathrm{e}}}& =\frac{3Q_{\nu }D{m}_{\mathrm{H}}}{4{\rho }_{\mathrm{d}}{a}_{\mathrm{d}}{\sigma }_{\mathrm{T}}}\\ & =7.1\times 10^2{\left(\frac{a_{\mathrm{d}}}{0.1\mu \mathrm{m}}\right)}^{-1}\left(\frac{Z}{Z_{\odot }}\right).\end{array}$$

$$⟨f_{\mathrm{Edd}}⟩(Z=0)\approx 1\%\left(\frac{M_{\mathrm{BH}}{n}_{\mathrm{H}}}{10^7{\mathrm{M}}_{\odot }\mathrm{c}{\mathrm{m}}^{-3}}\right).$$

$$\begin{array}{rl}{r}_{\mathrm{HII}}& ={\left(\frac{3f_{\mathrm{Edd}}{L}_{\mathrm{Edd}}}{4\pi n_{\mathrm{H}}^2{\alpha }_{\mathrm{B}}{\bar{e}}_{\mathrm{ion}}}\right)}^{\frac{1}{3}}\\ & =33\mathrm{pc}{\left(\frac{f_{\mathrm{Edd}}}{10^{-2}}\right)}^{\frac{1}{3}}{\left(\frac{n_{\mathrm{H}}}{10^2\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-\frac{2}{3}}{\left(\frac{M_{\mathrm{BH}}}{10^5{\mathrm{M}}_{\odot }}\right)}^{\frac{1}{3}},\end{array}$$

$$\tau =\frac{3\times 10^{-2}{Q}_{\nu }{m}_{\mathrm{H}}{n}_{\mathrm{H}}}{4{\rho }_{\mathrm{d}}{a}_{\mathrm{d}}}\left(\frac{Z_{\mathrm{crit}}}{Z_{\odot }}\right)r_{\mathrm{HII}}\approx 1.$$

$$\begin{array}{r}{Z}_{\mathrm{crit}}\approx 0.2Z_{\odot }{\left(\frac{f_{\mathrm{Edd}}}{10^{-2}}\right)}^{-\frac{1}{3}}{\left(\frac{M_{\mathrm{BH}}{n}_{\mathrm{H}}}{10^7{\mathrm{M}}_{\odot }\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-\frac{1}{3}}\\ \approx 0.2Z_{\odot }{\left(\frac{M_{\mathrm{BH}}{n}_{\mathrm{H}}}{10^7{\mathrm{M}}_{\odot }\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-\frac{2}{3}}.\end{array}$$

$$\int \frac{x_{\mathrm{HI}}{L}_{\nu }{\sigma }_{\nu }{e}^{-{\tau }_{\nu }}}{4\pi r^{2}h\nu }d\nu \sim {\alpha }_{\mathrm{B}}{n}_{\mathrm{H}}(1-x_{\mathrm{HI}}{)}^2.$$

$${\dot{M}}_{\mathrm{B},\mathrm{HII}}\sim 4\pi {\lambda }_{\mathrm{B}}{G}^2{M}_{\mathrm{BH}}^2{\rho }_{\mathrm{HII}}{c}_{\mathrm{s},\mathrm{HII}}^{-3},$$

$$\begin{array}{r}\dot{M}\sim 4\pi {\lambda }_{\mathrm{B}}{G}^2{\left(1-\frac{f_{\mathrm{d}}{\sigma }_{\mathrm{T}}L}{4\pi cG{M}_{\mathrm{BH}}{m}_{\mathrm{p}}}\right)}^2{M}_{\mathrm{BH}}^2{\rho }_{\mathrm{HII}}{c}_{\mathrm{s},\mathrm{HII}}^{-3},\\ ={\dot{M}}_{\mathrm{B},\mathrm{HII}}{\left[1-\frac{7.1\times 10^2{\sigma }_{\mathrm{T}}\eta c}{4\pi G{M}_{\mathrm{BH}}{m}_{\mathrm{p}}}\left(\frac{Z}{Z_{\odot }}\right)\dot{M}\right]}^2.\end{array}$$

$$\int \pi a_{\mathrm{d}}^2{Q}_{\nu }{F}_{\nu }d\nu =\int 4\pi a_{\mathrm{d}}^2{Q}_{\nu }\pi B_{\nu }(T_{\mathrm{d}})d\nu .$$

$$F\sim 4{\sigma }_{\mathrm{SB}}{T}_{\mathrm{d}}^4\bar{Q}(T_{\mathrm{d}}),$$

$$T_{\mathrm{d}}(r)={\left[\frac{L{e}^{-{\tau }_{\mathrm{d}}(r)}}{16\pi {\sigma }_{\mathrm{SB}}\bar{Q}(T_{\mathrm{d}})r^2}\right]}^{\frac{1}{4}}.$$

$$\begin{array}{r}{T}_{\mathrm{d}}(r)=4.7\times 10^2\mathrm{K}{\left(\frac{f_{\mathrm{Edd}}}{10^{-2}}\right)}^{\frac{1}{4}}{\left(\frac{M_{\mathrm{BH}}}{10^5{\mathrm{M}}_{\odot }}\right)}^{\frac{1}{4}}\\ \times {\left(\frac{\bar{Q}}{10^{-2}}\right)}^{-\frac{1}{4}}{\left(\frac{r}{0.1\mathrm{pc}}\right)}^{-\frac{1}{2}}.\end{array}$$

$$f_{\mathrm{Edd}}\sim 3\times 10^{-2}{\left(L_{\mathrm{bol},40}^{1/2}{f}_{14/140}\right)}^{2/3}.$$

$$M_{\mathrm{BH}}\sim 8\times 10^3{\mathrm{M}}_{\odot }{L}_{\mathrm{bol},40}\left(\frac{f_{\mathrm{Edd}}}{10^{-2}}\right).$$

$$\frac{dv}{dt}=-n_{\mathrm{H}}\frac{m_{\mathrm{H}}}{m_{\mathrm{d}}}\pi a_{\mathrm{d}}^2{v}^2.$$

$$\begin{array}{rl}\tau & \sim \frac{v}{|dv/dt|}=\frac{4}{3}\frac{{\rho }_{\mathrm{d}}{a}_{\mathrm{d}}}{n_{\mathrm{H}}{m}_{\mathrm{H}}v}\\ & =6.0\times 10^3\mathrm{yr}\left(\frac{a_{\mathrm{d}}}{0.1\mu \mathrm{m}}\right){\left(\frac{n_{\mathrm{H}}}{10^2\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1}{\left(\frac{v}{10\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{-1}\end{array}$$

$$\begin{array}{r}{r}_{\mathrm{sub}}=6.8\times 10^{-3}\mathrm{pc}{\left(\frac{f_{\mathrm{Edd}}}{10^{-2}}\right)}^{\frac{1}{2}}\\ \times {\left(\frac{M_{\mathrm{BH}}}{10^5{\mathrm{M}}_{\odot }}\right)}^{1/2}{\left[\frac{\bar{Q}(T_{\mathrm{sub}})}{0.2}\right]}^{-1/2}.\end{array}$$

$${\tau }_{\mathrm{sp}}\sim 1\times 10^5\mathrm{yr}\left[1+{\left(\frac{T}{10^6\mathrm{K}}\right)}^{-3}\right]\frac{(a_{\mathrm{d}}/0.1\mu \mathrm{m})}{(n_{\mathrm{H}}/1\mathrm{c}{\mathrm{m}}^3)}.$$

$$n_{\infty }\gtrsim 10^5\mathrm{c}{\mathrm{m}}^{-3}{\left(\frac{M_{\mathrm{BH}}}{10^5{\mathrm{M}}_{\odot }}\right)}^{-1}{\left(\frac{T}{10^4\mathrm{K}}\right)}^{3/2},$$

$$n_{\infty }\gtrsim 2\times 10^3\mathrm{c}{\mathrm{m}}^{-3}{\left(\frac{M_{\mathrm{BH}}}{10^5{\mathrm{M}}_{\odot }}\right)}^{-1}\left(\frac{T}{10^4\mathrm{K}}\right){\left(\frac{Z}{Z_{\odot }}\right)}^{-1}.$$

Full text

1. Introduction

Understanding the growth mechanism of massive black holes is one of the major challenges of modern astrophysics. It is known that the number density of supermassive black holes (SMBHs) steeply decreases in the early Universe (e.g., Richards et al., 2006). Yet, recent observations have detected SMBHs with masses $\gtrsim 10^{9}~{}{\rm M_{\odot}}$ even at cosmic times less than $\rm\sim 1~{}Gyr$ (Mortlock et al., 2011; Wu et al., 2015). However, the existence of SMBHs at high redshifts alone does not suggest clues about the initial seed mass and the subsequent growth history.

Various mechanisms have been suggested for the formation channel of SMBH seeds, that typically are assumed to have masses $M_{\rm BH}\sim 10^{5}~{}{\rm M_{\odot}}$: (1) growth from stellar mass BHs of Population III star remnants by gas accretion (Alvarez et al., 2009; Jeon et al., 2012; Park & Ricotti, 2013); (2) direct collapse of supermassive stars (Omukai, 2001; Begelman et al., 2006; Volonteri & Begelman, 2010; Agarwal et al., 2012; Latif et al., 2013; Inayoshi et al., 2014; Sugimura et al., 2014); (3) formation in dense star clusters via collisions among stars (Rees, 1978; Portegies Zwart & McMillan, 2002; Devecchi et al., 2012; Katz et al., 2015; Yajima & Khochfar, 2016). After the formation of SMBH seeds, it is widely believed that further growth proceeds though gas accretion. However, the accretion mechanism has not been understood yet because the accreting gas is subject to feedback from radiation emitted near the accretion disks around BHs.

One of the mechanisms limiting the gas accretion rate is the balance between gravitational attraction and radiation pressure on free-electron, the so-called Eddington limit. Cosmological simulations showed that massive black holes of $\sim 10^{5}~{}{\rm M_{\odot}}$ could grow up to supermassive ones of $\sim 10^{9}~{}{\rm M_{\odot}}$ by $z\sim 6$ using Eddington limited Bondi-Hoyle accretion prescription and a simple thermal feedback model (Li et al., 2007; Di Matteo et al., 2008, 2012; Sijacki et al., 2009). Di Matteo et al. (2012) carried out cosmological simulations in large cosmological volumes of $(0.75~{}\rm Gpc)^{3}$ in comoving unit and showed the accretion rate of SMBHs in rare massive galaxies was near the Eddington limit most of the time, resulting in SMBHs with $\sim 10^{9}~{}{\rm M_{\odot}}$ at $z\sim 7$. However, due to computational limitations, in such cosmological simulations gas dynamics at the Bondi radius (where the gravity of BHs is dominant) is not well resolved, thus requiring sub-grid feedback models. Therefore there is a large uncertainty in the estimation of the gas accretion rate.

Using high-resolution radiation hydrodynamics simulations resolving the Bondi radius, Milosavljević et al. (2009) showed that the neighboring gas is ionized by the radiation from a central BH and the thermal pressure of Hii regions pushed gas away from the BH against the gravity. As a result, the accretion rate was significantly suppressed even at lower luminosities than Eddington. Park & Ricotti (2011) showed that the gas accretion periodically changed and the time-averaged accretion rate was $\sim 1~{}\%$ of the Bondi rate regardless of some parameters, e.g., radiative efficiency, black hole mass, background density (see also, Park & Ricotti, 2012).

Inayoshi et al. (2016) suggested that gas accretion rate could exceed the Eddington limit when BHs accrete from extremely high-density gas clouds where the size of ionized bubble is smaller than Bondi radius (see also, Park et al., 2014a; Sakurai et al., 2016). In addition, assuming the anisotropic radiation feedback, Sugimura et al. (2016) showed that gas efficiently accretes onto a BH along the shadowed region and the accretion rate exceeds the Eddington limit. Thus, it is still unclear how much the growth of BHs is regulated by radiative feedback.

In local galaxies, it is well known that BH mass tightly correlates with bulge mass or velocity dispersion (Kormendy & Ho, 2013). This implies co-evolution of BHs with galaxies. Park et al. (2016) showed that the growth rate of massive BHs can be enhanced under the influence of gravitational potential of the bulge. In addition, as star formation proceeds, gas surrounding a BH is metal/dust enriched through type-I/II supernovae, and stellar winds. Observations of high-redshift quasars at $z\gtrsim 6$ suggest that a large mass in dust exists around SMBHs. Dust is detected via its thermal emission or dust extinction (Bertoldi et al., 2003; Priddey et al., 2003; Maiolino et al., 2004; Wang et al., 2013). In addition, recent discoveries of high-redshift sub-millimeter galaxies indicate that some galaxies can become dust rich at early times (Riechers et al., 2013; Watson et al., 2015). Theoretically, recent simulations show that the metallicity near the galactic centers could reach $\sim 0.01~{}Z_{\odot}$ even at $z\sim 10$ (e.g., Wise et al., 2012) or even higher depending on the halo mass (e.g., Ricotti & Gnedin, 2005; Ricotti et al., 2016).

The metallicity and dust amount of massive galaxies in an over-dense region could reach the level of solar-neighborhood even at $z\gtrsim 6$ (e.g., Yajima et al., 2015). Therefore BHs were likely to grow in a dusty medium even in the early Universe. If dust exists around a BH, radiation from the inner parts of an accretion disk around a BH can be obscured. This changes the observational properties of accreting BHs and the dynamics of accreting gas. In addition to dust opacity, the radiation force on dust can play a roll in determining the growth rate of BHs (Ciotti & Ostriker, 2007; Namekata et al., 2014; Hensley et al., 2014). However, the interplay between dust and the photo-ionization feedback, which is the main feedback mechanism suppressing the growth of stellar or intermediate-mass BHs (IMBHs), has not been studied in sufficient detail. In this work we investigate the impacts of dust on the growth of BHs by using one-dimensional radiation hydrodynamics simulations resolving the both the Bondi radius and ionized bubbles simultaneously. We also estimate self-consistently the thermal emission from dust in our time-dependent models and show that the emission from a hot dust component at $\sim 20\rm\mu m$ is prominent only during luminosity bursts, while a warmer dust component produces a flux at $100\rm\mu m$ that is roughly proportional to the mean accretion rate. We thus conclude that the flux ratio at $20\rm\mu m/100\rm\mu m$ is a good proxy for the Eddington ratio.

The paper is organized as follows. We describe our models in §2. In §3, we present simulation results that include time-average Eddington ratios, accretion histories with and without dust, dependences of the Eddington ratios on metallicity and BH mass, spectral energy distributions (SEDs) in the infrared (IR) considering dust thermal emission. We discuss the dust destruction processes and the condition for hyper-accretion in §5, and summarize our results in §6.

2. Model

We solve the dynamics of gas surrounding a BH under radiative feedback using one-dimensional radiation hydrodynamics simulations. In this work, we use a hydrodynamics code Zeus-MP (Stone & Norman, 1992; Hayes et al., 2006). Park & Ricotti (2011) incorporated the radiative transfer of X-ray and UV photons and chemical reactions of primordial gas into the Zeus-MP. Here we furthermore incorporate dust attenuation and radiation pressure on dust to the code. In spherical symmetric coordinates the basic equations, i.e., the conservations of mass, momentum and energy, are:

[TABLE]

where $f_{\rm rad}$ is radiation force, $\Gamma$ is the heating rate by photo-ionization of hydrogen and helium, and $\Lambda$ is the radiative cooling rate. In this work, we consider radiation force on dust ($f_{\rm rad}^{\rm dust}$), free electrons ($f_{\rm rad}^{\rm e}$), and neutral hydrogen ($f_{\rm rad}^{\rm HI}$). For single-sized grains of dust the optical depth is

[TABLE]

where $a_{\rm d}$ is dust radius, $m_{\rm d}$ is mass of a dust grain, $Q_{\rm\nu}$ is the absorption coefficient to geometrical cross section, $m_{\rm H}$ is hydrogen mass, $\rho_{\rm d}$ is mass density of a dust grain, and $D$ is dust-to-gas mass ratio. If the wavelength of the radiation is shorter than $\sim 2\pi a_{\rm d}$, then $Q_{\rm\nu}\sim 1$ (Draine & Lee, 1984). Here we consider compact spherical dust grains, i.e., $m_{\rm d}=4\pi a_{\rm d}^{3}\rho_{\rm d}/3$. Here we assume that the dust-to-gas mass ratio is proportional to metallicity as $D\equiv M_{\rm dust}/M_{\rm gas}=0.01\left({Z}/{Z_{\odot}}\right)$, motivated by observations of a nearly constant dust-to-metal mass ratio in local galaxies (Draine et al., 2007).

Observations of local galaxies indicate that dust has a continuum size distributions. For example, the dust in the Milky Way shows a power-law size distribution ${dn_{\rm d}}/{da_{\rm d}}\propto a_{\rm d}^{-3.5}$, often referred to as the MRN distribution (Mathis et al., 1977). However, for sake of simplicity here we use a single dust size model with a fiducial size $a_{\rm d}=0.1~{}\rm\mu m$. The choice of the fiducial size is motivated below. At wavelengths $\lambda<2\pi a_{\rm d}$, where $Q_{\rm\nu}=1$, the opacity of a power law distribution of dust radii ${dn}/{da}\propto a^{-\alpha}$, is the same of a single size dust model for

[TABLE]

where $a_{\rm d,1}$ is the dust size of the equivalent single-size model, $a_{\rm d}^{\rm min}$ and $a_{\rm d}^{\rm max}$ are minimum and maximum grain sizes in the grain distribution. Assuming realistic minimum and maximum dust sizes of $a_{\rm d}^{\rm min}=8.1\times 10^{-3}~{}\rm\mu m$ and $a_{\rm d}^{\rm max}=1.0~{}\rm\mu m$ in a MRN model ( i.e., $\alpha=3.5$), the equivalent single dust size model has $a_{\rm d,1}=0.1~{}\rm\mu m$.

Near active galactic nuclei (AGNs), the main components of dust are graphite and silicate because of their high sublimation temperature $T\gtrsim 1500~{}\rm K$, whereas icy dust can be easily sublimated. Here we consider a mix of graphite and silicate dust grains with the mass ratio of $1:1$. The difference of mass density between the graphite and silicate is small. In addition, when $Q_{\nu}=1$, the absorption cross section depends only on the surface area of a dust grain. Therefore, the simulation results are rather insensitive to the assumed mass ratio of the dust components.

Within fully ionized region, the radiation force on dust can be larger than that from Compton scattering on electrons by a factor

[TABLE]

UV photons absorbed by dust are reemitted as IR photons via thermal emission of dust. However, the IR photons can escape from neighbor gas without further interaction with dust due to the lower absorption cross-section to the IR wavelengths. If the hydrogen column density exceeds $10^{22}~{}\rm cm^{-2}$, the IR photons can be absorbed by dust and impart additional momentum to the gas. In this work, we focus on spatial scales up to an ionized bubble. On these scales the absorption of IR photons by dust is negligible. Therefore, we neglect dust opacity to IR light.

Aside from the addition of dust physics, the simulations have the same physics and initial conditions as the ones presented in Park & Ricotti (2012). For our fiducial simulations we investigate the gas dynamics around BHs of $10^{5}~{}{\rm M_{\odot}}$ embedded in gas clouds with uniform initial gas density. The radiative luminosity is $L=\eta\dot{m}c^{2}$ with a constant radiative efficiency $\eta=0.1$ and $\dot{m}$ estimated at the inner boundary in the logarithmically spaced radial grid. We assume a power-law spectrum $L_{\rm\nu}\propto\nu^{-1.5}$ with the frequency range $\nu\geq\nu_{\rm L}$, where $\nu_{\rm L}$ is the Lyman-limit frequency. However the results can be easily scaled to different mass BHs and different densities as the problem is basically scale-free. Park & Ricotti (2012) showed that simulations with a fixed $M_{\rm BH}n_{\rm H}$ show identical behavior once the length and time scales are re-normalized appropriately. Park & Ricotti (2012) used the condition of $M_{\rm BH}n_{\rm H}=10^{7}~{}\rm{\rm M_{\odot}}~{}cm^{-3}$ as a fiducial run. For this parameters combination, if gas accretes onto the BH at the Bondi rate, the luminosity is close to the Eddington luminosity. Here we use the same fiducial run with $M_{\rm BH}n_{\rm H}=10^{7}~{}\rm{\rm M_{\odot}}~{}cm^{-3}$ to compare the impact of dust attenuation and radiation force on the gas accretion onto BHs to the dust-free case. Table 1 summarizes three different models we ran to understand the impact of dust.: (1) dust-free case (M5-Z0), (2) dusty case with $Z=Z_{\odot}$ without radiation pressure on dust (M5-Z1), and (3) dusty case with $Z=Z_{\odot}$ including radiation pressure on dust (M5-Z1rad). In all simulations we include $f_{\rm rad}^{\rm e}$ and $f_{\rm rad}^{\rm HI}$. The inner and outer boundaries are $2\times 10^{-2}$ pc and $2\times 10^{3}~{}\rm pc$, respectively. The number of cells is $n_{\rm cell}=256$. We have checked that our results are convergent and do not change when increasing the number of cells from $128$ to $400$. The cell size changes with radial distance with the equatorial ratio. To check that indeed the results are scalable to different BH masses we run simulations with BH masses of $10^{3},10^{4}$ and $10^{6}~{}{\rm M_{\odot}}$ keeping constant $M_{\rm BH}n_{\rm H}=10^{7}~{}\rm{\rm M_{\odot}}~{}cm^{-3}$.

3. Results

3.1. Time evolution of luminosity under radiative feedback

Figure 1 shows time-averaged BH luminosity $\left<f_{\rm Edd}\right>\equiv{\left<L\right>}/{L_{\rm Edd}}$, in units of the Eddington luminosity $L_{\rm Edd}\cong 1.3\times 10^{43}~{}{\rm erg~{}s^{-1}}~{}\left({M_{\rm BH}}/{10^{5}~{}{\rm M_{\odot}}}\right)$, for the simulations in Table 1. We calculate $\left<f_{\rm Edd}\right>$ at various metallicities ($Z=0-1.0~{}Z_{\odot}$), and initial gas densities $n_{\rm\infty}=10-1000~{}\rm cm^{-3}$. Regardless of dust amount, $\left<f_{\rm Edd}\right>$ linearly increases with gas density. This means that the accretion rate normalized by the Bondi rate is constant as a function of the gas density. The Eddington ratio $\left<f_{\rm Edd}\right>$ for gas of primordial composition and $n_{\rm\infty}=100~{}\rm cm^{-3}$ (i.e., M5-Z0 run), is $\sim 1\%$ of the Bondi rate. Thus for the case with $Z=0~{}Z_{\odot}$, if $\left<f_{\rm Edd}\right><1$ we have:

[TABLE]

This value is consistent with previous studies (Milosavljević et al., 2009; Park & Ricotti, 2011, 2012). As the metallicity and dust content increases, $\left<f_{\rm Edd}\right>$ decreases. At $Z=10^{-2}~{}Z_{\odot}$, $\left<f_{\rm Edd}\right>$ is similar to the cases without dust. When the metallicity is higher than $\sim 0.1~{}Z_{\odot}$, $\left<f_{\rm Edd}\right>$ becomes smaller than that of primordial gas by a factor greater than $\sim 2$. The M5-Z1rad run with radiation pressure on dust and solar metallicity has $\left<f_{\rm Edd}\right>\sim 10^{-3}$, while the same model without radiation pressure on dust (M5-Z1) has $\left<f_{\rm Edd}\right>\sim 7\times 10^{-3}$ that is similar to the M5-Z0 run without dust. This suggests that the radiation force on dust significantly reduces the accretion rate, whereas the effect of dust opacity alone is secondary.

Figure 2 shows the time evolution of the luminosity for the simulations in Table 1. As shown in previous works (Park & Ricotti, 2011), the accretion of primordial gas periodically changes. During the bursts the luminosity is near the Eddington limit. When the luminosity is high, the Hii region expands far beyond the Bondi radius. At the ionizing front, the low-density hot ionized gas is almost in pressure equilibrium with the outside Hi gas, and suppress the gas inflow from the neutral region. This leads to a decrease of accretion rate. The minimum luminosity becomes smaller than the maximum one by about 5-6 orders of magnitude.

After each burst, the gas density within the Hii bubble decreases with time. In addition, the lowered luminosity allows recombination of hydrogen in Hii region with radiative cooling. Due to the recombination, the ionization degree and temperature in the Hii region decrease. These break the pressure equilibrium and allows gas inflow, resulting in the burst of luminosity again.

The mid panel shows the luminosity of M5-Z1 run (with dust). The burst cycle is shorter than that of M5-Z0 run (without dust). This is because the cycle time scale is proportional to the sound crossing time over ionized region $\sim r_{\rm HII}/c_{\rm s}$ (Park & Ricotti, 2012), and the size of Hii region ($r_{\rm HII}$) is decreased due to the dust attenuation. In order for the dust attenuation to work, the optical depth of dust in a Hii region should be higher than unity. We can roughly estimate the critical metallicity (and therefore dust abundance) that can have an effect at reducing the size of Hii region and therefore the period between bursts. First we estimate the size of Strömgren sphere:

[TABLE]

where $\alpha_{\rm B}$ is the case-B recombination coefficient, and $\bar{e}_{\rm ion}$ is the mean energy of ionizing photons. We here set $T=7\times 10^{4}~{}\rm K$ as the temperature of ionized region in the above estimation. For a power-law spectrum, $\bar{e}_{\rm ion}$ is estimated by $\frac{\alpha}{\alpha-1}\times 13.6~{}{\rm eV}$, where $\alpha$ is a slope of the power-law spectrum, and it is $40.8~{}\rm eV$ for $\alpha=1.5$. Using the estimated size of ionized region above, we can derive the critical metallicity for the dust attenuation as

[TABLE]

Thus, the critical metallicity is given by

[TABLE]

If the metallicity is higher than $Z_{\rm crit}$, the size of Hii region decreases and so does the period between bursts. Note that, even in the case with the dust attenuation (but no radiation pressure on dust), the size of ionized bubble is larger than the Bondi radius. Therefore, the gas density inside the Hii region is roughly the same as the run with no dust (M5-Z0), and so is the mean accretion rate. However, the peak luminosities are somewhat smaller than in the M5-Z0 run due to the smaller Hii bubbles. This might lead to the lower $\left<f_{\rm Edd}\right>$ than the primordial gas cases by a factor $\sim 2$.

The lower panel shows the luminosity for the run that includes the effect of both dust opacity and radiation pressure on dust (M5-Z1rad). Unlike M5-Z0 and M5-Z1 runs, the gas accretion of M5-Z1rad proceeds more gently and the luminosity varies between the burst and quiescent phases by one or two orders of magnitude (instead of 5 orders of magnitude as in the other runs). The maximum and minimum luminosities are $\sim 10^{-2}$ and $10^{-4}$ of $L_{\rm Edd}$, respectively. Due to the high radiation force on dust that, as shown in Eq. (4), at solar metallicity is nearly 700 times larger than Compton scattering on electrons, the gas inflow is significantly suppressed when $f_{\rm Edd}\gtrsim 10^{-3}$, resulting in the lower maximum luminosity than in the M5-Z0 run by one order of magnitude. In the case of M5-Z1rad run, the period between bursts is not determined by the sound crossing time of the Hii region ($r_{\rm HII}/c_{\rm s}$). As will be shown in Sec. 3.2, the Hii bubble in the M5-Z1rad run maintains almost constant size, whereas the one in the M5-Z0 run recombines and collapses when the luminosity drops to the minimum (see also, Park & Ricotti, 2012). In the case of primordial gas, the expansion of the Hii bubbles significantly reduces the gas density and accretion rate, and then the rapid gas inflow during the collapse of the ionized bubble allows the luminosity to reach nearly the Eddington limit. On the other hand, in the case of dusty gas, the dust attenuation regulates the expansion of the Hii bubble, and the radiation force suppresses the gas inflow significantly. Thus, the dust plays a role in diminishing the large periodic variability of the luminosity and lead to more gentle gas accretion onto the BH.

Figure 3 presents the density profiles at the high and low luminosity phases for runs M5-Z0 and M5-Z1rad. At the high luminosity phase, the density of M5-Z0 run steeply increases at $r\lesssim 1~{}\rm pc$, and it reaches $n_{\rm H}\sim 10^{4}~{}\rm cm^{-3}$ at the inner boundary. The inflow velocity exceeds sound speed at $\sim r_{\rm B}/2$, and hence the density profiles are roughly $\propto r^{-3/2}$ as in the Bondi profile (free falling gas).

On the other hand, the density of M5-Z1rad slowly increases as the radial distance decreases. In the case with dust, the net inward force is significantly reduced due to the radiation force. By introducing an effective gravitational constant as $G^{\prime}=G\left(1-f_{\rm d}f_{\rm Edd}\right)$, we estimate the transonic radius by $r_{\rm t}\sim G^{\prime}M_{\rm BH}/2c_{\rm s}^{2}$. When $f_{\rm Edd}\gtrsim 10^{-3}$, the transonic radius is smaller than the inner boundary of the calculation box. As a result, the gas inflow of M5-Z1rad is always sub-sonic in our calculation box. In addition, the temperature at $r\lesssim 0.1~{}\rm pc$ increases as the radial distance decreases due to compression heating. Therefore, the pressure gradient force also works on suppressing the gas inflow. Thus the gas density does not increase steeply.

As shown in the lower panel of the figure, the gas density is decreased due to the photo-ionization feedback during the quiescent low-luminosity phase. At the ionization front, the gas is almost in pressure equilibrium with the gas outside the Hi region, i.e., $2n_{\rm HII}T_{\rm HII}\sim n_{\rm HI}T_{\rm HI}$, where the factor $\sim 2$ is due to the increased particle number from ionization. Since the temperature of the Hii region reaches $\sim 7\times 10^{4}~{}\rm K$ via photo-ionization of hydrogen and helium (Park & Ricotti, 2011), the density of the Hii region is roughly lower than the ambient Hi gas by a factor $\sim 14$. The gas density is $\lesssim 100~{}\rm cm^{-3}$ even at the inner boundary in the simulation. Therefore the optical depth by dust inside the Hii region does not exceed unity as will be discussed in Sec. 3.2.

Figure 4 shows the gas outward acceleration due to radiation pressure normalized by the gravitational acceleration. Both $f_{\rm rad}^{\rm dust}$ and $f_{\rm rad}^{\rm e}$ are roughly constant inside the Hii region. The number ratio of dust grains to free electrons is almost constant inside the Hii region. Therefore, the $f_{\rm rad}^{\rm dust}$ is higher than $f_{\rm rad}^{\rm e}$ by a factor $\sim 700$ as shown by Eq. (4). During the burst of accretion the dust optical depth exceeds unity inside the Hii region, because the gas density increases near the BH as shown in Figure 3. Therefore the radiation pressure decreases as the distance increases due to the dust opacity. $f_{\rm rad}^{\rm e}$ sharply drops outside the ionizing front because the electron abundance decreases, while $f_{\rm rad}^{\rm dust}$ gradually decreases since dust exists both in the ionized and neutral gas. The radiation force on Hi, $f_{\rm rad}^{\rm HI}$, monotonically increases with the radial distance in the Hii region. Where $\tau<1$, $f_{\rm rad}^{\rm HI}/f_{\rm grav}$ is simply proportional to the neutral fraction $x_{\rm HI}$. Inside the Hii region, hydrogen is in ionization equilibrium:

[TABLE]

For $\tau\ll 1$ and $x_{\rm HI}\ll 1$, $x_{\rm HI}\propto n_{\rm H}r^{2}$. $f_{\rm rad}^{\rm HI}$ increases somewhat more slowly than $r^{2}$. This is due to the decrease of gas density as the radial distance increases at $r\lesssim 10~{}\rm pc$. Near the ionization front, $f_{\rm rad}^{\rm HI}$ steeply increases with increasing $x_{\rm HI}$, and decreases exponentially in the neutral region. 111 We notice that there is a dip in $f_{\rm rad}^{\rm HI}$ near the ionizing front. At this radius, the flux steeply decreases as ${\rm exp}(-N_{\rm HI}\sigma_{\rm\nu,HI})$, while the neutral fraction increases at $\sim 1-2$ cells behind, producing the dip. This is due to the finite resolution. Increasing the number of cells makes this dip smaller. Since the radiation force on Hi is negligible, this feature does not affect our results.

3.2. Size of ionized bubble

The time evolution of the position of the ionization front (defined where $x_{e}=50~{}\%$) is presented in Fig. 5. The size of ionized bubble $r_{\rm HII}$ changes with the luminosity in the case of the M5-Z0 run without dust. When the pressure equilibrium at the ionization front breaks due to the decreased density and luminosity, the neutral gas free falls into the BH. In this case, the inflowing gas remains neutral until reaching near the BH, resulting in the collapse of the ionized bubble just before the next bursts of luminosity.

On the other hand, in the dusty gas simulation (M5-Z1rad) the size of the ionized bubble remains nearly constant, with $r_{\rm HII}\lesssim 20~{}\rm pc$. Even during the high luminosity phases, the size is several times smaller than that of the dust-free case (M5-Z0). This is because dust opacity regulates the evolution of the size of the ionized bubble. Figure 6 shows the optical depth of dust during the high and low luminosity phases. When the luminosity is high, the optical depth of dusty gas reaches unity at $r\sim 10~{}\rm pc$ (M5-Z1rad). Therefore the ionizing front cannot propagate far beyond $10~{}\rm pc$, whereas in the dust-free case (M5-Z0) the ionizing front (I-front) reaches $\sim 40~{}\rm pc$. During the low luminosity phase, the optical depth does not exceed unity in the ionized bubble due to the lower gas density as shown in Fig. 3. In this case, the size of ionized bubble is determined by the ionizing luminosity. The gas density steeply increases over the transition from ionized to neutral regions. This density jump significantly increases the optical depth due to the high-density dust in the neutral region. Therefore, regardless of different luminosities, the optical depth exceeds unity just behind the ionizing front.

3.3. Metallicity dependence

As star formation proceeds, the interstellar medium in galaxies becomes metal/dust enriched by type-II supernovae (e.g., Wise et al., 2012). Recent observations indicated a large dispersion of dust amount in high-redshift galaxies (Watson et al., 2015). Here we study the dependence of the BH accretion luminosity on the gas metallicity and therefore the dust abundance (since we assume a constant dust-to-metal mass ratio). Figure 7 shows $\left<f_{\rm Edd}\right>$ as a function of metallicity. As the metallicity increases, $\left<f_{\rm Edd}\right>$ monotonically decreases, and ranges from $\left<f_{\rm Edd}\right>\sim 10^{-2}$ at $Z=10^{-2}~{}Z_{\odot}$ to $\left<f_{\rm Edd}\right>\sim 10^{-3}$ at $Z=Z_{\odot}$. We examine the metallicity dependence based on the assumption of steady spherical inflow of ionized gas ( i.e., Bondi accretion inside the ionized bubble). In the case without dust, we estimate the accretion rate of Bondi-like inflow in the ionized region ($\dot{M}_{\rm B,HII}$) as follows:

[TABLE]

where $\lambda_{\rm B}$ is the dimensionless mass accretion rate, $\rho_{\rm HII}$ and $c_{\rm s,HII}$ are the density and the sound speed of ionized gas. The value of $\lambda_{\rm B}$ depends on the polytropic index $\gamma$ of the equation of state, i.e., $P=K\rho^{\gamma}$, and ranges from $\sim 1.12$ for an isothermal gas to $1/4$ for the adiabatic case. When we assume that the gas inflow is isothermal and consists of ionized hydrogen and helium with the temperature $T_{\rm HII}\sim 7\times 10^{4}~{}\rm K$, the accretion rate from Eq. (10) is $\sim 770$ times smaller than the Bondi rate calculate from the density and temperature of the ambient medium. This reduction of the accretion rate is roughly similar to our simulation results and previous works (Milosavljević et al., 2009; Park & Ricotti, 2012; Sugimura et al., 2016). Note that, however, the inflow of primordial gas causes periodic bursts as shown in the Sec. 3.1. These bursts can increase the time averaged accretion rate with respect to that estimated above (see also, Park & Ricotti, 2012).

Next, we modify the above accretion rate by taking into account the radiation pressure on dust. Due to radiation pressure, the net inward force is reduced. In the optically thin regime ($\tau_{\rm d}<1$), the fractional reduction of the inward gravitational force is independent to the radial distance, and depends only on the BH luminosity. In this case, the accretion rate is

[TABLE]

By solving this second order equation, we derive an analytical estimate of the accretion rate and luminosity. The analytically estimated Eddington ratios are shown as dashed ($n_{\rm\infty}=100~{}\rm cm^{-3}$) and dotted ($n_{\rm\infty}=300~{}\rm cm^{-3}$) lines in the Figure 7. As shown in the figure, the modified Bondi rate considering the radiation force roughly explains the metallicity dependence of the accretion rate in the simulations with radiation pressure on dust. Note that, $f_{\rm Edd}$ in the high metallicity and high density cases differ from the value estimated above, whereas it in the case with $n_{\rm\infty}=100~{}\rm cm^{-3}$ our analytical estimate is good fit to the simulation results. When the density and metallicity are high, the size of Hii bubble is close to the Bondi radius (e.g., Park & Ricotti, 2012). Therefore, the gas density in the ionized region can become somewhat higher than that estimated assuming pressure equilibrium across the ionization front due to the BH gravitational potential, resulting in higher accretion rate.

3.4. Mass dependence

Here we study the BH mass dependence on the gas accretion rate for our fiducial runs with $M_{\rm BH}n_{\rm H}=10^{7}~{}\rm{\rm M_{\odot}}~{}cm^{-3}$. Park & Ricotti (2012) showed that the accretion rate normalized by the Bondi rate did not change for different BH masses when $M_{\rm BH}n_{\rm H}=10^{7}~{}\rm{\rm M_{\odot}}~{}cm^{-3}$. Figure 8 shows the time averaged Eddington ratios for different BH masses. As shown in Park & Ricotti (2012), $\left<f_{\rm Edd}\right>\sim 1\%$ independently of the BH mass in the case without dust. We find that the same result holds even when including dust. Thus, regardless of the dust and metallicity of the gas, we confirm that the results are scale-free under the condition $M_{\rm BH}n_{\rm H}=10^{7}~{}\rm{\rm M_{\odot}}~{}cm^{-3}$.

3.5. Thermal emission from dust

Dust releases the energy absorbed at UV wavelengths as thermal emission in the IR. Here we investigate the IR properties of the dusty gas accreting onto the BH. The dust temperature $T_{\rm d}$ can be estimated by assuming radiative equilibrium:

[TABLE]

This equation assumes equilibrium between dust absorption of the radiation emitted by the BH (left-hand side) and the thermal emission from dust (right-hand side). By solving this equation, we estimate the dust temperature of each gas shell. We assume $Q_{\nu}=1$ at UV wavelengths for the left-hand side of the equation, and use $Q_{\nu}$ estimated in Laor & Draine (1993) for the right-hand side. In order to estimate the thermal emission by dust we post-process the simulation results, calculating the absorbed flux by dust. In the estimation of the dust thermal emission, we neglect the absorption of ionizing photons by hydrogen and helium. This is because the absorbed energy is converted into recombination and cooling radiation, e.g., Ly$\alpha$ photons. These recombination UV photons are eventually absorbed by dust. In addition, by considering only dust absorption, the IR properties do not change significantly even if we consider different choices for the SED from the accretion disk, e.g., a SED extending to soft-UV range. Figure 9 shows the dust temperature as a function of distance from the BH. During the high-accretion phase, the dust temperature is $\sim 800~{}\rm K$ near the inner boundary of the simulation, and decreases as the distance increases due to the geometrical attenuation of the flux $\propto r^{-2}$. Due to the strong attenuation of UV flux outside the I-front, the dust temperature sharply drops at $r\gtrsim 10~{}\rm pc$. In regions where $T_{\rm d}\gtrsim 100~{}\rm K$, the dust temperature of silicate is somewhat lower than graphite. This is because $Q_{\rm\nu}$ of silicates is higher than graphite in that temperature range, thus the efficient thermal photon emissivity results in the lower temperature.

In order to better understand the radial profile of the dust temperature, we re-write Equation (12) as:

[TABLE]

where $F=\int F_{\nu}d\nu$ is radiation flux and $\bar{Q}(T_{\rm d})={\int B_{\nu}(T_{\rm d})Q_{\rm\nu}d\nu}/{\int B_{\nu}(T_{\rm d})d\nu}$ is the frequency-averaged dust absorption coefficient. Therefore, the dust temperature is

[TABLE]

For $\tau\ll 1$, this equation is scaled as follows:

[TABLE]

As shown by the figure, the dust temperature roughly decreases as $\propto r^{-1/2}$. However, note that, $\bar{Q}$ decreases as the dust temperature decreases. Therefore, the radial profile of the dust temperature is somewhat shallower than $r^{-1/2}$.

The lower panel of Fig. 9 shows the radial dependence of the fraction of energy absorbed by the dust per unit radial log bin, normalized by the total luminosity. During the quiescent accretion phase, most of energy is absorbed by dust near the ionizing front. This is because the optical depth at $r<10~{}\rm pc$ is much smaller than unity as shown in Fig. 6. Therefore, most of IR emission comes from cold dust with $T_{\rm dust}\sim 20~{}\rm K$ near the ionizing front. As the accretion rate increases, a fraction of the radiation is absorbed by dust near the inner boundary of calculation box due to higher gas and dust density. This results in a large contribution to the IR luminosity by the inner hot dust with $T_{\rm dust}\gtrsim 100~{}\rm K$.

By integrating the thermal emission of each spherical shell, we derive SEDs. Figure 10 shows the SEDs at $f_{\rm Edd}=1.1\times 10^{-2}$ (red line), $0.8\times 10^{-3}$ (green line), and $1.2\times 10^{-4}$ (blue line). The red and blue lines correspond to the ones in Fig. 9. We also show the frequency ranges of Spitzer, JWST, Hershel and ALMA. Due to the contribution to the IR luminosity by hot dust, the SED during the high accretion phase shows stronger IR emission at $\lambda\lesssim 50~{}\rm\mu m$. The specific luminosity at these wavelengths varies by several orders of magnitude between the quiescent phase with low Eddington ratio and the bursty phase at higher Eddington ratio. On the other hand, at $\lambda\gtrsim 100~{}\rm\mu m$ the specific luminosity remains nearly constant during the accretion cycle, because the cold dust with $T_{\rm d}\lesssim 50~{}\rm K$ near the ionization front, always absorbs and reprocesses a large fraction of the radiation. Recently, Shimizu et al. (2017) showed that the IR SEDs of AGNs selected from their hard X-ray luminosity using Swift BAT, correlates with the X-ray luminosity $L_{\rm X}$: with increasing $L_{\rm X}$, the flux at short wavelength $\lambda\sim 10~{}\rm\mu m$ also increases. If the BH masses of the observed AGNs in the sample are similar each other, these variations of the SEDs can be explained by the different phases of accretion and the Eddington ratios as explained above. An estimate of the masses and Eddington ratios of these hard X-ray selected AGNs is currently in progress (T. Shimizu, private communication). With the caveat of being able to isolate the AGN contribution from the stellar contribution in the IR SEDs, the results of this study can be used to test our model predictions.

We also show the contributions of graphite and silicate grains, separately. The stretching resonance of astronomical silicate grains produces the bump in the SED at $\lambda\sim 9.8~{}\rm\mu m$. At shorter wavelengths ($\lambda\lesssim 8~{}\rm\mu m$), the contribution from graphite grains is dominant. In this work, we have not considered re-absorption of the thermal emission by dust, because the absorption cross-section of dust to IR photons is quite small and we assume that the ambient gas is optically thin at IR wavelengths. However, if the hydrogen column density of the gas reservoir fueling the BH is much higher than $10^{22}~{}\rm cm^{-2}$, some IR radiation emitted by hot dust can be absorbed by the dust and re-processed to thermal emission at lower dust temperature. This can somewhat suppress the IR flux at short wavelength.

Fig. 11 shows the relation between the Eddington ratio, $f_{\rm Edd}$, and the flux density ratio between $14$ and $140~{}\rm\mu m$ ($f_{14/140}$) multiplied by the square root of bolometric luminosity, $L_{\rm bol,40}$ (in units of $10^{40}~{}\rm erg~{}s^{-1}$). As shown in Fig. 10, the flux at $\lambda\lesssim 10~{}\rm\mu m$ depends sensitively on the presence of a hot dust component. Therefore we choose the flux at $14~{}\rm\mu m$ as diagnostic of the high-accretion phase. The wavelengths $14$ and $140~{}\rm\mu m$ correspond to the peak emission of two black body spectra with temperatures of $207.0$ K and $20.7~{}\rm K$, respectively. We find that this flux ratio tightly correlates with $f_{\rm Edd}$. The correlation between $f_{\rm Edd}$ and $f_{14/140}$ varies with the BH mass, and dust size, mainly because the temperature of the warm and hot dust components depend on these physical parameters. However, $f_{\rm Edd}$ tightly correlates with $f_{14/140}$ regardless of the assumed initial background density. Yet, we can use our knowledge of $L_{\rm bol,40}$ to correct for these effects and recover $f_{\rm Edd}$ regardless of the different physical parameters of the SMBH, that are generally unknown observationally. The different panels in the figure illustrate the tight correlation between $f_{\rm Edd}$ and $L_{\rm bol,40}^{1/2}f_{14/140}$, that holds by varying different physical parameters by many orders of magnitude. We find that the correlation can be roughly fitted by the following relationship:

[TABLE]

Therefore, keeping into account that our model is still very simplistic, it may be possible to use IR observations of warm and hot dust obscuring SMBHs to estimate $L_{\rm bol,40}$ and $f_{14/140}$, and therefore derive $f_{\rm Edd}$ and the SMBH mass as

[TABLE]

Since $L_{\rm bol,40}\equiv f_{\rm Edd}L_{\rm Edd}\propto f_{\rm Edd}M_{\rm BH}$, we can use Eq. (16) to show that $f_{\rm Edd}\propto M_{\rm BH}^{1/2}f_{14/140}$. Therefore, for a fixed BH mass, we indeed find $f_{\rm Edd}\propto f_{14/140}$. Instead, keeping constant $M_{\rm BH}n_{\rm\infty}$, while varying the density and/or BH mass, we find $r_{\rm HII}\propto f_{\rm Edd}^{1/3}M_{\rm BH}$ (see Eq. 6). Using Eq. 15, the dust temperature at $r_{\rm HII}$ is $T_{\rm d}(r_{\rm HII})\propto M_{\rm BH}^{-1/4}$ (here we ignore a weak dependence on $f_{\rm Edd}$). Therefore, the weak mass dependence in the figure indicates $f_{14/140}$ decreases as $T_{\rm d}(r_{\rm HII})$ decrease with increasing the BH mass (and $L_{\rm bol,40}$).

We also study the SEDs for different dust size models, $a_{\rm d}=0.05$ and $0.02~{}\rm\mu m$. The right panel of Fig. 11 shows the correlations in the cases of different dust sizes. Unlike the tight correlations found for different BH masses and densities, models with the different dust sizes produce a systematic change in shape of the relationship. Since $Q_{\rm abs}$ at IR band decreases with dust size, smaller dust grains become hotter when irradiated by the same UV flux as shown in Eq. (15). Moreover, the size of the Hii region becomes smaller as the dust size decreases (see Eq. (2)). Thus, models with smaller dust tend to have higher $f_{14/140}$ for a fixed $f_{\rm Edd}$. In the same panel we also show the correlation for the low-metallicity (low-dust) case with $Z=0.1~{}Z_{\odot}$. As the metallicity decreases, the BH luminosity and the ionize bubble radius, $r_{\rm HII}$, increase. As a result, the dust temperature does not change significantly, producing the same correlation as in the fiducial case with $Z=1~{}Z_{\odot}$. Therefore we suggest the uncertainty of metallicity is not needed to be considered.

Can JWST or ALMA observe thermal emission from the dust obscuring massive BHs in nearby galaxies? For low-redshift AGNs, the flux from hot dust at rest-frame $\lambda\sim 14~{}\rm\mu m$ can be observed by JWST, while the warm dust emission at $\lambda\sim 140~{}\rm\mu m$ is in the wavelength range covered by Spitzer. Even if we measure the flux density ratio between $14$ and $450~{}\rm\mu m$ (instead of $140~{}\rm\mu m$) that is in the waveband covered by ALMA, as similar relationship between $f_{\rm Edd}$ and the flux ratio is found. Therefore ALMA can also be used to trace the cold dust at $z\sim 0$.

Using the specific luminosities derived above, we estimate the flux densities from dust around massive BHs at a specific redshift $z$ as $F_{\rm\nu}=(1+z)L_{\rm\nu 0}/4\pi D_{\rm L}^{2}$ where $D_{\rm L}$ is the luminosity distance, $\nu$ is the frequency in the observer rest frame, and $\nu_{0}=\nu(1+z)$ is the frequency in the galaxy’s rest frame. Here we artificially put BHs at $z=0.05-10$, and discuss their observability by JWST or ALMA. Figure 12 shows the flux densities as a function of redshift. The dotted lines show the sensitivities of JWST and ALMA. We show the flux densities at $14~{}\rm\mu m$ ($F_{14}$), $140~{}\rm\mu m$ ($F_{140}$) and $450~{}\rm\mu m$ ($F_{450}$) in the observer frame. Due to the negative K-correction, the flux density at $450~{}\rm\mu m$ decreases slowly with increasing redshift, whereas that at $14~{}\rm\mu m$ decreases more steeply at $z\gtrsim 1$. At higher redshift, the flux at $14~{}\rm\mu m$ in observer rest frame corresponds to a wavelength in the galaxy’s rest frame of $14~{}{\rm\mu m}/(1+z)$ which can be shorter than the peak wavelength of modified black body of even hot dust, resulting in the decreases of the flux. As a result, the flux density at $450~{}\rm\mu m$ becomes higher than that at $14~{}\rm\mu m$ at $z>1$. The relation between $F_{14}$ and $F_{450}$ changes depending on the BH mass and $f_{\rm Edd}$. As shown above, $F_{14}$ is lower than $F_{140}$ or $F_{450}$ when $f_{\rm Edd}$ is low. Even at high-accretion phase ($f_{\rm Edd}\sim 10^{-2}$), $F_{450}$ is higher than $F_{14}$ in the case of massive BH with $10^{6}~{}\rm{\rm M_{\odot}}$. Under the constraint $M_{\rm BH}n_{\rm H}={\rm const}$, we have $r_{\rm HII}\propto M_{\rm BH}$ (see Eq. 6). Therefore, the dust temperature at the ionization front is lower in massive BHs (see Eq. 15), leading to a higher $F_{450}$.

In summary, we find that only massive BH of mass $\gtrsim 10^{6}~{}\rm{\rm M_{\odot}}$ at $z\lesssim 0.1$ can be observed by ALMA. On the other hand, JWST will allow us to observe the dust thermal emission from more distant massive BHs. It will be able to probe hot dust obscuring BHs of $10^{6}~{}{\rm M_{\odot}}$ with $f_{\rm Edd}\sim 10^{-2}$ up to $z\sim 0.5$. When $f_{\rm Edd}\sim 10^{-2}$, even a $10^{4}~{}{\rm M_{\odot}}$ BH at $z\lesssim 0.1$ can be observed. Thus we suggest JWST and ALMA will be powerful tools to probe IMBHs via the observation of dust thermal emission. In general, the IR flux increases as the product $M_{\rm BH}n_{\rm\infty}$ increases. For instance, for the case of a $M_{\rm BH}=10^{6}~{}{\rm M_{\odot}}$ accreting from gas with density $n_{\rm\infty}=100~{}\rm cm^{-3}$ and metallicity $Z=0.01~{}Z_{\odot}$, ALMA and JWST can probe BHs up to $z\sim 0.5$. However, even upcoming new telescopes will be difficult to observe IMBHs at $z\gg 1$ in the IR.

In our simulations, we have investigated the accretion dynamics of spherically symmetric clouds with isotropic radiation field in which the accretion rate and luminosity can be suppressed significantly. Recently, Sugimura et al. (2016) suggested the gas efficiently accreted onto a BH if the radiation is anisotropic due to shadowing effects near the edge of the BH accretion disk or torus. In this case, $f_{\rm Edd}$ can become $\gtrsim 1$, and even relatively massive BHs at high-redshifts $z\gtrsim 1$ may be observable by JWST or ALMA, although the dust temperature and SEDs can differ from our current works.

4. Discussion

4.1. Dust decoupling

In this work we have assumed that the motion of dust is completely coupled to the gas through collisions between gas and dust particles. Here we simply estimate the coupling time scale. At first, we consider the momentum equation for the relative velocity $v$ between dust and gas:

[TABLE]

Therefore, the time scale for the dust coupling with gas is roughly

[TABLE]

This coupling timescale is much shorter than dynamical timescale and the period between bursts in our current simulations. Thus the assumption of the complete coupling is reasonable.

4.2. Dust destruction

The dust formation/destruction processes have not been considered in our simulations. Even at innermost cell ($=0.02~{}\rm pc$), the dust temperature does not exceed the dust sublimation temperature. For the sublimation temperatures of graphite grains, $T_{\rm sub}\sim 1800~{}\rm K$, we estimate the destruction radius as follows:

[TABLE]

Thus, the dust can survive sublimation by photoheating even inside the inner boundary of the calculation box.

On the other hand, the destruction by thermal sputtering process, i.e., collisions between dust and gas, is likely to be effective near the BH because of the high-density and temperature found at the innermost cells in our simulations. The destruction time scale by the thermal sputtering is estimated (Draine & Salpeter, 1979; Draine, 2011) as follows:

[TABLE]

During the quiescent phase, the density, temperature and inflow velocity at the innermost cell are $\sim 10^{3}~{}{\rm cm^{-3}},~{}3\times 10^{5}~{}\rm K$ and $\lesssim 10~{}\rm km~{}s^{-1}$, respectively. Therefore, during this phase in the duty cycle, the sputtering time scale is shorter than the inflow time scale, $t_{\rm dyn}\sim r/v_{\rm in}$ and the dust is likely to be destroyed inside the inner boundary of the calculation box. However, during the peaks of accretion, the temperature somewhat decreases and the velocity increases, resulting in the longer sputtering time scale. Therefore, during the accretion burst dust can survive even close to a BH, and can absorb UV radiation emitted during the burst. However, as shown in Figure 9, the contribution to the opacity by dust absorption at radii smaller than the inner boundary in our simulations is not significant. In addition, Compton heating increases the temperature of gas at such the scale (Park et al., 2014b), leading to enhancement of dust destruction. Yet, photo-ionization heating of metals is not included in this work. This heating can expand the destruction radius. These detailed heating processes will be considered in our future works.

In this work, we adopt a single-size dust model in the calculations. When considering a realistic dust size distribution, the destruction radius depends on the dust size. This is because small dust grains can have higher temperature, hence they reach the sublimation temperature at larger distance. In addition, the destruction time scale by thermal sputtering is proportional to the size. Thus, smaller dust grains are more easily destroyed near the BH compared to larger dust grains. This indicates that the dust size distribution changes with radial distance.

4.3. Revisited condition for hyper-Eddington accretion

In this work we have investigated the density dependence up to $n_{\rm H}=1000~{}\rm cm^{-3}$. This condition allows the ionizing front to propagate far from the Bondi radius in the case without dust. Recently, Inayoshi et al. (2016) suggested that BHs could grow at super-Eddington rate if they are embedded in very high-density gas clouds where the size of ionized bubble is smaller than the Bondi radius. This is because the ionized bubble is gradually shrunk as the gas density near the ionization front increases due to the gravitational force by the BH. Finally the ionized bubbles disappear, resulting in no radiation force on electrons. As a result the growth rate of BH becomes close to the Bondi rate which can be much higher than the Eddington limit, i.e., so-called the hyper-Eddington accretion. They estimated the critical gas density for the hyper-Eddington accretion by comparing the Strömgren radius to the Bondi radius. By calculating the Strömgren radius for the density inside the ionized region (which is smaller than $n_{\rm\infty}$ by a factor $\sim 2(T_{\rm HII}/T_{\rm HI})$), we derive the critical density:

[TABLE]

where we have assumed $T_{\rm HII}=7\times 10^{4}~{}\rm K$. On the other hand, in the case of dusty gas, the size of ionized region is regulated by the dust opacity if the metallicity is higher than a critical value ( i.e., if $Z>0.1Z_{\odot}$). In this regime, the size of ionization region is close to the photon mean free path, i.e., $r_{\rm HII}\sim(f_{\rm d}\sigma_{\rm T}n_{\rm HII})^{-1}$. Thus, the critical density for super-Eddington accretion in a dusty gas is:

[TABLE]

This critical density is much lower than that for the primordial gas case. However, the effect of dust in the regime of super-Eddington accretion, in addition to reducing the size of the ionized bubble, is to reduce the accretion rate due to the radiation pressure on dust. Radiation pressure is effective at radii larger than the size of the ionization region, as long as $\tau_{\rm d}\lesssim 1$. Therefore, if the ram pressure force can overcome the radiation force on dust at the radius of $\tau_{\rm d}\sim 1$, the dusty gas may be able to efficiently accrete onto the BH. The accretion dynamics of high-density dusty clouds will be investigated in future works with further improvements to the code.

5. Summary

In this paper we have studied the accretion of dusty gas onto an intermediate mass black hole (IMBH) by using one-dimensional radiation hydrodynamics simulations. Park & Ricotti (2011) showed that the growth of BH was significantly regulated due to the photo-ionization feedback in the case of primordial gas. As star formation proceeds, high-redshift galaxies are metal/dust enriched due to type-II supernovae. Recent observations have detected dust-rich galaxies even at $z\gtrsim 6$ (e.g., Riechers et al., 2013). This indicates some BHs grew in dusty gas with the feedback. In this paper, we investigate the effects of dust on the growth of BH and observable diagnostics due to the thermal emission from heated dust. Dust affects the accretion rate and period of the bursts of accretion mainly because of radiation pressure on dust but also because of the dust opacity that reduces the size of the ionization region.

By assuming as a fiducial model accretion onto a BHs of $10^{5}~{}{\rm M_{\odot}}$ embedded in a uniform density medium, we investigate the dependence of the BH growth rate on the gas density (in the range $n_{\rm\infty}=10$ to $1000~{}\rm cm^{-3}$), and on the metallicity (in the range $Z=0$ to $1~{}Z_{\odot}$). We find that the accretion of dusty gas onto IMBHs proceeds gently with small fluctuations of the accretion rate, whereas that of primordial gas causes periodic bursts. For dust-to-gas mass ratio similar to the solar neighborhood, the time averaged luminosity becomes smaller than that for primordial gas by one order of magnitude. The time averaged Eddington ratio is $\left<f_{\rm Edd}\right>\sim 10^{-3}$ for the initial gas density $n_{\rm\infty}=100~{}\rm cm^{-3}$. Our calculations show that the effect of dust opacity is secondary with respect to radiation pressure on dust. Neglecting radiation pressure on dust but including the effect of dust opacity, the growth rate of IMBHs and $\left<f_{\rm Edd}\right>$ are smaller than that of primordial gas by a factor $\lesssim 2$. For both primordial and dusty clouds, the Eddington ratio, $\left<f_{\rm Edd}\right>$, linearly increases with the initial gas density. In addition, assuming the constraint $(M_{\rm BH}/{\rm M_{\odot}})(n_{\rm\infty}/{\rm cm^{-3}})=10^{7}$, we study the dependence of the growth rate on the BH mass. We show that $\left<f_{\rm Edd}\right>$ is constant for the different BH masses in the both cases with and without dust.

Finally, we derive the SEDs at IR bands by calculating dust thermal emission. Our modeled SEDs show that the flux ratio at $\lambda\lesssim 20~{}\rm\mu m$ and $\gtrsim 100~{}\rm\mu m$ depends sensitively on the Eddington ratio, but is nearly independent of the other parameters in the problem. This is because at high Eddington ratios the thermal emission from hot dust near the BH produces a higher flux density at $\lesssim 20~{}\rm\mu m$. While the emission at $\gtrsim 100~{}\rm\mu m$, that is produced by warmer dust further out, near the ionization front, is nearly independent of $\left<f_{\rm Edd}\right>$. Therefore, we suggest that the combinations of MIR observations by JWST and FIR observation by Hershel or ALMA can be combined to provide a novel method to estimate the Eddington ratio of BHs throughout their duty cycle, including their short bursty phase and their longer quiescent phase.

Acknowledgments

The numerical simulations were performed on the computer cluster, Draco, at Frontier Research Institute for Interdisciplinary Sciences of Tohoku University. This work is supported in part by MEXT/JSPS KAKENHI Grant Number 15H06022 (HY) and 15J03873 (KS).