Systematic two-dimensional radiation-hydrodynamic simulations of super-Eddington accretion flow and outflow: Comparison with the slim disk model

Abstract

To what extent can the one-dimensional slim disk model reproduce the multi-dimensional results of global radiation-hydrodynamic simulations of super-Eddington accretion? With this question in mind, we perform a systematic simulation study of accretion flow onto a non-spinning black hole for a variety of black hole masses of (10–10⁷) M_⊙ and mass accretion rates of (1.4 × 10²–5.6 × 10³) L_Edd/c² (with L_Edd and c being the Eddington luminosity and the speed of light). In order to adequately resolve large-scale outflow structure, we extensively expand a simulation box to cover the space of 3000 r_S (with r_S being the Schwarzschild radius), larger than those in most previous studies, so that we can put relatively large angular momentum on the gas injected from the outer simulation boundary. The adopted Keplerian radius, at which the centrifugal force balances the gravitational force, is r_K = 300 r_S. The injected mass first falls and is accumulated at around this radius and then slowly accretes toward the central black hole via viscosity. We simulate such accretion processes, taking inverse and bulk Compton scattering into account. The simulated accretion flow is in a quasi-steady state inside r_qss ∼ 200 r_S. Within this radius the flow properties are, on the whole, in good agreement with those described by the slim disk model except that the radial density profile of the underlying disk is much flatter, ρ ∝ r^−0.73 (cf. ρ ∝ r^−3/2 in the slim disk model), due probably to efficient convection. We find very weak outflow from inside r ∼ 200 r_S, unlike the previous studies.

1 Introduction

It has long been believed that the Eddington luminosity is a classical limit to the luminosities of any accreting objects. However, we now know that this is no longer the case from both observational and theoretical grounds. In fact, there is growing observational evidence supporting the existence of super-Eddington accretion in several distinct classes of objects. In parallel with extensive observational studies, multi-dimensional simulations are being performed by a number of groups. Super-Eddington (or supercritical) accretion is attracting much attention among researchers.

Good candidates for super-Eddington accretors are ultra-luminous X-ray sources (ULXs), bright X-ray compact sources of X-ray luminosity of 10³⁹–10⁴¹ erg s⁻¹ (see Kaaret et al. 2017 for a recent review, and references therein). Quite a few ULXs have been discovered so far in off-nuclear regions of nearby galaxies, and the number is rapidly increasing. There are two main ideas to explain their high luminosities: sub-Eddington accretion onto intermediate-mass black holes (Makishima et al. 2000; Miller et al. 2004), and super-Eddington accretion onto stellar-mass black holes (Watarai et al. 2001; King et al. 2001). (Both scenarios require high accretion rates, but this is a separate issue and we do not go into the details here.) The discovery of ULX pulsars in which the central object is a neutron star supports the latter scenario (NGC 7793 P13: Fürst et al. 2016, Israel et al. 2017b; NGC 5907 ULX: Israel et al. 2017a; NGC 300 ULX-1: Kosec et al. 2018).

Other super-Eddington accretors are found in narrow-line Seyfert 1 galaxies (NLS1s), bright micro-quasars such as GRS 1915+105, ultra-soft X-ray sources (ULSs), etc. The NLS1s harbor less massive central black holes (with mass ≲ 10⁷ M_⊙; see Boller 2000) than broad-line Seyfert 1 galaxies (BLS1s). That is, NLS1s tend to have higher Eddington ratios than BLS1s with similar luminosities, and thus super-Eddington accretion is more feasible in the former (Wang & Zhou 1999; Mineshige et al. 2000). Jin et al. (2017), for example, analyzed RX J0439.6−5311 (NLS1) using the multi-wavelength spectrum and estimated the accretion rate to be ∼71 L_Edd/c⁻² in the outer disk (for the black hole mass ∼1 × 10⁷ M_⊙). The black hole binary GRS 1915+105 is also known to stay occasionally in the super-Eddington phase (see, e.g., Done et al. 2007; Vierdayanti et al. 2010).

One of the most prominent features of super-Eddington flow is photon trapping (Katz 1977; Begelman 1978). That is, photons generated in the accretion disk tend to be directly swallowed by a black hole before diffusing toward the surface of the disk. This occurs when accretion rate is high, since then the diffusion time t_diff becomes longer than the accretion time t_acc. The equation t_diff = t_acc leads to the photon-trapping radius r_trap (e.g., Ohsuga et al. 2002):

\begin{equation} r_{\rm trap}\equiv \frac{3}{2} \frac{\dot{M}}{L_{\rm Edd}/c^{2}} \frac{H}{r}r_{\rm S}. \end{equation}

(1)

Here, |$\dot{M}$| is the mass accretion rate and H is the half-thickness of the accretion disk. The diffused photons in the region r < r_trap accrete into the black hole with gases because photons cannot escape from the accretion disk. We wish to stress here that photon trapping is essentially a multi-dimensional effect (Ohsuga et al. 2002). This is a key issue to be discussed in the present paper.

The slim disk model is the one-dimensional accretion disk model including the photon-trapping effect as the advection of the photon entropy in the energy equation (Abramowicz et al. 1988; see chapter 10 of Kato et al. 2008 for a concise review). The basic equations for the radial structure of the slim disk are derived from the Navier–Stokes equations under the vertically one-zone approximation; that is, physical quantities are integrated in the |$z$|-direction (perpendicular to the equatorial plane). The slim disk model is a numerical model, but an approximate analytical expression is available, as obtained by Watarai (2006), which gives a simple formula describing the parameter dependence of mass density, temperature, and velocity on the equatorial plane of the slim disk model.

Outflow is another signature characterizing super-Eddington flow—see a pioneering discussion by Shakura and Sunyaev (1973). When the disk luminosity exceeds the Eddington luminosity L_Edd, this means that the radiation force is greater than the gravity by the definition of the Eddington luminosity. Part of the accretion flow gas is blown out from the disk, accelerated by radiation pressure, and then the outflow occurs at a high mass accretion rate. Note that such multi-dimensional motion as large-scale circulation (convection) and outflow are not explicitly considered in the slim disk model. We thus need to perform multi-dimensional radiation-hydrodynamic (RHD) simulations.

Ohsuga et al. (2005) pioneered two-dimensional RHD simulations for super-Eddington accretion flow—see also Kawashima et al. (2009, 2012), Hashizume et al. (2015), Ogawa et al. (2017), and Kitaki et al. (2017) for RHD simulations; Ohsuga et al. (2009), Ohsuga and Mineshige (2011), and Jiang, Stone, and Davis (2014) for radiation-magnetohydrodynamic (R-MHD) simulations; Sa̧dowski et al. (2014, 2015), McKinney et al. (2014), Fragile, Olejar, and Anninos (2014), and Takahashi et al. (2016) for general relativistic radiation-magnetohydrodynamic (GR-R-MHD) simulations. In these simulation studies the authors adopted somewhat unrealistic situations; that is, they commonly assumed relatively small angular momentum of accreting gas, which is either injected from the outer simulation boundary or provided from the initial gaseous torus. This was necessary for numerical reasons, since otherwise it would take an enormous computation time (corresponding to a long viscous timescale) to be completed within a few months. This leads to a quite narrow viscous accretion region in a quasi-steady state (within a few tens of the Schwarzschild radii in most cases), which makes it difficult to compare with the slim disk model.

In the present study, therefore, the initial angular momentum is set larger. We calculate the structure of the super-Eddington accretion flow and outflow by means of two-dimensional (2D) RHD simulations for a variety of black hole masses (M_BH) and mass accretion rates (⁠|$\dot{M}_{\rm BH}$|⁠), and compare the simulation results with those calculated based on the slim disk model. The main feature in this study compared with previous papers is a systematic study of the scaling relations produced by numerical simulations, which is due more to the volume of parameter space spanned (along with a larger radial extension) than to a higher spatial resolution (some of the studies presented in the introduction actually used better resolution). We also obtained the fitting formulas of the super-Eddington accretion disk for the first time.

The plan of this paper is as follows: We first explain our models and calculation methods in the next section. We then show our main results in section 3, and finish with a discussion in section 4.

2 Models and numerical methods

2.1 Radiation-hydrodynamic simulations

In the present study, we consider super-Eddington accretion flow and outflow onto a black hole by injecting mass from the outer simulation boundary at a constant rate of |$\dot{M}_{\rm input}$| with a certain amount of angular momentum. (The parameter values will be specified in subsection 2.2.) The flux-limited diffusion approximation is adopted (Levermore & Pomraning 1981; Turner & Stone 2001). We also adopt the α viscosity prescription (Shakura & Sunyaev 1973) and assign α = 0.1 throughout the present study. General relativistic effects are incorporated by adopting the pseudo-Newtonian potential (Paczyńsky & Wiita 1980).

The basic equations and numerical methods are the same as in Kawashima et al. (2009, 2012), but upgraded to solve the energy equations with the implicit method. This 2D RHD code solves the axisymmetric two-dimensional radiation-hydrodynamic equations in the spherical coordinates (x, y, |$z$|⁠) = (r sin θcos ϕ, r sin θsin ϕ, r cos θ), where the aziumuthal angle ϕ is set to be constant. The continuity equation is given by

\begin{eqnarray} \frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \boldsymbol {v}\right)=0. \end{eqnarray}

(2)

Here, ρ is the gas mass density and |$\boldsymbol {v}=(v_{\rm r},v_{\theta },v_{\phi })$| is the gas velocity. The equations of motion are written as

\begin{eqnarray} &&\!\!\!\frac{\partial (\rho v_{\rm r})} {\partial t}+\nabla \cdot \left(\rho v_{\rm r}\boldsymbol {v}\right)\nonumber\\ &&\quad = -\frac{\partial p}{\partial r}+\rho \left[\frac{v_{\theta }^{2}}{r}+\frac{v_{\phi }^{2}}{r}-\frac{GM_{\rm BH}}{(r-r_{\rm S})^{2}}\right] +\frac{\chi }{c}F_{0\rm r}, \end{eqnarray}

(3)

\begin{eqnarray} \frac{\partial (\rho rv_{\theta })}{\partial t}+\nabla \cdot \left(\rho rv_{\theta }\boldsymbol {v}\right) & = & -\frac{\partial p}{\partial \theta }+\rho v_{\phi }^{2}\cot \theta +r\frac{\chi }{c}F_{0\theta }, \end{eqnarray}

(4)

\begin{eqnarray} \frac{\partial (\rho r\sin \theta v_{\phi })}{\partial t}+\nabla \cdot \left(\rho r\sin \theta v_{\phi }\boldsymbol {v}\right) & = & \frac{1}{r^{2}}\frac{\partial }{\partial r}\left(r^{3}\sin \theta t_{\rm r\phi }\right), \end{eqnarray}

(5)

where p is the gas pressure and r_S ≡ 2GM_BH/c² is the Schwarzschild radius, with G being the gravitational constant and c being the speed of light. χ = κ + ρσ_T/m_p is the total opacity, where κ is the free–free and free–bound absorption opacity (Rybicki & Lightman 1979), σ_T is the cross-section of Thomson scattering, and m_p is the proton mass. |$\boldsymbol {F}_{0}=(F_{0\rm r},F_{0\theta })$| is the radiative flux in the comoving frame, where the suffix 0 represents quantities in the comoving frame. Using the dynamical viscous coefficient η, t_rϕ is the viscous stress tensor described as

\begin{eqnarray} t_{\rm r\phi }&=&\eta r \frac{\partial }{\partial r}\left(\frac{v_{\phi }}{r} \right), \end{eqnarray}

(6)

\begin{eqnarray} \eta &=&\alpha \frac{p+\lambda E_{0}}{\Omega _{\rm K}}. \end{eqnarray}

(7)

Here, Ω_K is the Keplerian angular speed, E₀ is the radiation energy density, and λ represents the flux limiter of the flux-limited diffusion approximation (Levermore & Pomraning 1981; Turner & Stone 2001).

The energy equations of the gas and the radiation are given by

\begin{eqnarray} \frac{\partial e}{\partial t}{+}\nabla \cdot \left(e\boldsymbol {v}\right)&{=}&-p\nabla \cdot \boldsymbol {v}-4\pi \kappa B{+}c\kappa E_{0} +\Phi _{\rm vis}-\Gamma _{\rm Comp}, \end{eqnarray}

(8)

\begin{eqnarray} &&\!\!\!\frac{\partial E_{0}}{\partial t}+\nabla \cdot \left(E_{0}\boldsymbol {v}\right)\nonumber\\ &&\quad=-\nabla \cdot \boldsymbol {F}_{0}-\nabla \boldsymbol {v}:\boldsymbol {P}_{0}+4\pi \kappa B-c\kappa E_{0} +\Gamma _{\rm Comp}. \end{eqnarray}

(9)

Here, e is the internal energy density, which is linked to the thermal pressure by the ideal gas equation of state, p = (γ − 1)e = ρk_BT_gas/(μm_p) with γ = 5/3 being the specific heat ratio, k_B being the Boltzmann constant, μ = 0.5 is the mean molecular weight, and T_gas is the gas temperature. |$B=\sigma _{\rm SB}T_{\rm gas}^{4}/\pi$| is the blackbody intensity, where σ_SB is the Stefan–Boltzmann constant. |$\boldsymbol {P}_{0}$| is the radiation pressure tensor, and Φ_vis is the viscous dissipative function written as

\begin{eqnarray} \Phi _{\rm vis}=\eta \left[r\frac{\partial }{\partial r}\left(\frac{v_{\phi }}{r} \right)\right]^{2} . \end{eqnarray}

(10)

The Compton cooling/heating rate Γ_Comp is described as

\begin{eqnarray} \Gamma _{\rm Comp}&=&4\sigma _{\rm T}c\frac{k_{\rm B}\left(T_{\rm gas}-T_{\rm rad}\right)}{m_{\rm e}c^{2}}\left(\frac{\rho }{m_{\rm p}} \right)E_{0} . \end{eqnarray}

(11)

Here, m_e is the electron mass and T_rad ≡ (E₀/a)^1/4 is the radiation temperature with the radiation constant a = 4σ_SB/c.

The simulation settings are roughly the same as those in Kitaki et al. (2017), except for the initial conditions and the size of the simulation box. The computational box is set by r_in = 2r_S ≤ r ≤ r_out = 3000 r_S, and 0 ≤ θ ≤ π/2. Grid points are uniformly distributed in logarithm in the radial direction, ▵log₁₀r = (log₁₀r_out − log₁₀r_in)/N_r, and uniformly distributed in cos θ in the azimuthal direction, ▵cos θ = 1/N_θ, where the numbers of grid points are (N_r, N_θ) = (192, 192) throughout the present study.

2.2 Initial conditions and calculated models

All the parameter values of the calculated models are summarized in table 1. We start the calculations with an empty space around a black hole with mass M_BH, though we initially add a hot optically thin atmosphere with negligible mass for numerical reasons. Mass is injected continuously with a constant rate of |$\dot{M}_{\rm input}$| through the outer disk boundary at r = r_out and 0.45π ≤ θ ≤ 0.5π. The black hole mass and mass injection rate are free parameters and are set to be M_BH = 10, 10⁴, and 10⁷ (M_⊙), and |$\dot{M}_{\rm input} = 3\times 10^2$|⁠, 10³, 5 × 10³, 10⁴, and 10⁵ (L_Edd/c²). We set the injected mass to have an angular momentum corresponding to the Keplerian angular momentum at the Keplerian radius, r = r_K, which is a free parameter— that is, the initial specific angular momentum is |$\sqrt{GM_{\rm BH} r_{\rm K}}$|⁠. We thus expect that inflow material first falls toward the center and forms a rotating gaseous ring at around r ∼ r_K, from which the material slowly accretes inward via a viscous diffusion process. We allow mass to go out freely through the outer boundary at r = r_out and 0 ≤ θ ≤ 0.45π, and assume that mass at r = r_in is absorbed.

Table 1.

Calculated models.*

Model name	MBH	|$\dot{M}_{\rm input}$|	rK
	[M⊙]	[LEdd/c2]	[rS]
(a11)		300	100
(a12)		103	300
(a13)	101	5 × 103	300
(a14)		104	300
(a15)		105	300
(a41)		300	100
(a42)		103	300
(a43)	104	5 × 103	300
(a44)		104	300
(a45)		105	300
(a71)		300	100
(a72)		103	300
(a73)	107	5 × 103	300
(a74)		104	300
(a75)		105	300

*Here, M_BH is the black hole mass, |$\dot{M}_{\rm input}$| is the mass injection rate, and r_K is the Kepler radius (see text). The other parameters are the same among all models: the inner radius is r_in = 2r_S, the outer radius is r_out = 3000 r_S, the α viscosity is α = 0.1, and the metallicity is Z = Z_⊙.

3 Results

3.1 Overall flow structure

We first give an overview of the simulated flow structure plotted in figures 1 and 2. These are the two-dimensional color contours of the time-averaged mass density and gas temperature of all the calculated models, respectively. As is clearly seen in figure 1, the injected gas accumulates at around r_K (∼300 r_S, except for the panels at the left end) and forms a puffed-up structure inside ∼10³ r_S. The accumulated gas then slowly accretes toward the center. After a sufficient time over the viscous timescale, an accretion flow settles in a quasi-steady state (Ohsuga et al. 2005).

Fig. 1.

Time-averaged density contours of super-Eddington accretion flow onto black holes with masses of M_BH = 10¹, 10⁴, and 10⁷ (M_⊙) from the top to the bottom, for mass injection rates of |$\dot{M}_{\rm input} = 300$|⁠, 10³, 5 × 10³, 10⁴, and 10⁵(L_Edd/c²) from left to right. Note the different color scales for different black hole masses. Each panel looks similar if we adjust the color scale according to the density normalization relationship of |$\rho _0\propto M_{\rm BH}^{-1}$| (see Kitaki et al. 2017). (Color online)

Fig. 2.

As figure 1 but for the time-averaged gas temperature distributions. The color scales are the same in all the panels. (Color online)

If we consider the panels in figure 1 from left to right, we notice that the puffed-up region grows with the increase of the mass injection rate, although the Keplerian radius is kept the same (except in the panels at the left end). This is because the trapping radius increases with an increase in the mass accretion rate, obeying |$r_{\rm trap} \sim (\dot{M}_{\rm BH}c^2/L_{\rm Edd})r_{\rm S}$|⁠. If we consider the panels from top to bottom in figure 1, on the other hand, we see no big changes, as long as we change the color scale according to the relationship |$\rho _0 \propto M_{\rm BH}^{-1}$|⁠. This scaling of density with the black hole mass is linked to mass conservation (Kitaki et al. 2017).

Let us next compare the different models in terms of the temperature distribution displayed in figure 2. As we see from the left to the right panels, we understand that the high-temperature regions (with red color) shrink as |$\dot{M}_{\rm BH}$| increases, while the low-temperature regions expand. As we see from the upper to the bottom panels, the color is much bluer near the equatorial plane. This reflects the fact that the higher the black hole mass, the cooler the accretion disk becomes, because the flux is almost equal to the Eddington luminosity divided by surface area in the super-Eddington accretion disk; in other words, the disk temperature at a fixed r/r_S obeys |$\sigma _{\rm SB} T^{4}_{\rm disk} \propto L_{\rm Edd}/(2\pi r^{2})\propto M^{-1}_{\rm BH}$| (see the discussion in Kitaki et al. 2017).

The overall flow structure looks rather similar to those obtained by previous studies, e.g., Kawashima et al. (2012) and Ogawa et al. (2017). As was noted by Kitaki et al. (2017), we may distinguish three characteristic zones: accretion disk, funnel, and over-heated regions (see their figure 2). The first one is the accretion disk located at and around the equatorial plane, r ∼ r_in − 1000 r_S and θ ∼ 30°–90°—see Model (a12) in figure 1. This disk is puffed up by the radiation pressure, and gas falls toward the center by transporting the angular momentum outward via the viscous process.

The second zone is the funnel region located around the polar (rotational) axis, r ∼ r_in − r_out and θ ∼ 0°–30°—see Model (a12) in figure 2. The Thomson scattering optical depth of the funnel in the |$z$|-direction is τ_e ∼ 1. The funnel is characterized by high gas temperature, k_BT_funnel ≳ 10 keV (see figure 2), and by a very fast velocity, |$v$|_r ≳ 0.2 c. The funnels in the first two columns from the left, which show models with low mass injection rates of |$\dot{M}_{\rm input}=300$| and 10³ (× L_Edd/c²), are widely extended from the polar direction to the direction of θ ∼ 45°, whereas the funnels in the third and fourth columns from the left, which show models with large mass injection rates of |$\dot{M}_{\rm input}= 5\times 10^{3}$| and 10⁴ (× L_Edd/c²), are rather narrow around the polar axis (see figure 2). This is because the funnel is collimated by the thickness of the puffed-up accretion disk when the accretion rate is relatively high.

The third zone is the over-heated region near the black hole at r ∼ 5 r_S and θ ∼ 45° (see, e.g., figure 2 in Kitaki et al. 2017 for the details). The gas temperature is very high (≳ 10 keV) there.

3.2 Mass accretion rate, outflow rate, and net flow rate

One of the most advantageous points in the present study is that we take a relatively large simulation box in order to increase the angular momentum of the accretion material as much as possible, compared with those assigned in the previous studies. This is a great advantage in investigating the super-Eddington accretion flow and outflow, since we can achieve a quasi-steady state of inflow and outflow in a much wider spatial range than before.

We discuss here the properties of the simulated multi-dimensional accretion flow structure in quasi-steady state in comparison with that calculated by the one-dimensional slim disk model. For this purpose, it is essential to know to what extent a quasi-steady state is realized and thus to calculate mass inflow rate and outflow rate at each radius. In the present study, we calculate these rates according to

\begin{eqnarray} \dot{M}_{\rm in}(r)&=&\int _{4\pi }d\Omega r^{2}\rho (r,\theta ){\rm min}\lbrace v_{\rm r}(r,\theta ),0\rbrace , \end{eqnarray}

(12)

\begin{eqnarray} \dot{M}_{\rm out}(r)&=&\int _{4\pi }d\Omega r^{2}\rho (r,\theta ){\rm max}\lbrace v_{\rm r}(r,\theta ),0\rbrace . \end{eqnarray}

(13)

Here, |$\dot{M}_{\rm in} (<0)$| is the (time-averaged) mass inflow rate, and |$\dot{M}_{\rm out} (>0)$| is the (time-averaged) mass outflow rate. (Note that here we calculate the mass outflow rate, irrespective of the outflow speed; that is, we do not distinguish the outflow that can reach infinity from the one that cannot.) The net flow rate is then calculated by

\begin{eqnarray} \dot{M}_{\rm net}(r)&=&\dot{M}_{\rm out}(r)+\dot{M}_{\rm in}(r). \end{eqnarray}

(14)

By the quasi-steady flow we mean the flow in which |$\dot{M}_{\rm net}$| has no (or negligibly small) radial dependence.

Figure 3 shows typical examples of the radial distributions of the inflow rate, the outflow rate, and the net rate for Models (a12) and (a72). (We wish to note that the other models show similar results.) We see that |$\dot{M}_{\rm net}$| is nearly constant in the spatial range of r ∼ 2–200 r_S. We also notice that both the absolute value of the mass inflow rate and the mass outflow rate increase significantly beyond the radius of r ∼ 200 r_S. Such features can be understood since the injected gas from the outer boundary at r = r_out is accumulated around the Keplerian radius, r ∼ r_K, due to the angular momentum barrier. We thus conclude that a quasi-steady flow is achieved in the range of r ∼ 2–200 r_S.

Fig. 3.

Radial profiles of the (time-averaged) mass inflow rate, |$\dot{M}_{\rm in}$| (red), the (time-averaged) mass outflow rate, |$\dot{M}_{\rm out}$| (green), and the net flow rate, |$\dot{M}_{\rm net}\equiv \dot{M}_{\rm in} + \dot{M}_{\rm out}$| (blue) for Models (a12) and (a72) in the upper and lower panels, respectively. Note that the horizontal scale is logarithmic. The vertical black lines at r ∼ 200 r_S indicate the radius of r_qss, inside which the flow is in a quasi-steady state (see the text for the precise definition). The inset in each panel shows the same but in a narrower spatial range of r ≤ 200 r_S, in which the horizontal axis is now in a linear scale. (Color online)

Let us define the quasi-steady radius, r_qss, as the radius inside which a quasi-steady state is realized; i.e., the net flow rate is constant, |$\dot{M}_{\rm net}(r)\propto r^{0}$|⁠. Technically, we evaluate r_qss in the following way:

1. We first use an estimated value of r_qss, say, r_qss = 100 r_S.

2. We search for the radial mesh index i_qss such that the following inequality holds: |$r_{i_{\rm qss}}\le r_{\rm qss} < r_{i_{\rm qss}+1}$|⁠.

3. The mean net flow rate, |$\langle \dot{M}_{\rm net}\rangle$|⁠, and its standard deviation, σ_net, are calculated by averaging the net flow rate over the range between r = r_in and |$r_{i_{\rm qss}+1}$|⁠; that is,

   \begin{eqnarray} \langle \dot{M}_{\rm net}\rangle &\equiv & \frac{1}{i_{\rm qss}}\sum _{i=1}^{i_{\rm qss}}\dot{M}_{\rm net}(r_{i}), \end{eqnarray}

   (15)

   \begin{eqnarray} \sigma _{\rm net}&\equiv & \sqrt{\frac{1}{i_{\rm qss}}\sum _{i=1}^{i_{\rm qss}} \left[\dot{M}_{\rm net}(r_{i})-\langle \dot{M}_{\rm net}\rangle \right]^{2}}. \end{eqnarray}

   (16)
4. If the relationship |$\left|\dot{M}_{\rm net}(r_{i_{\rm qss}+5})\right| \ge \left|\langle \dot{M}_{\rm net}\rangle \right|+1.5\sigma _{\rm net}$| holds for the first time, we define the radius r_qss as

   \begin{eqnarray} r_{\rm qss}&\equiv & (r_{i_{\rm qss}}+r_{i_{\rm qss}+1})/2 \end{eqnarray}

   (17)

   and end the loop. Otherwise, we repeat the same procedurefrom the first step but by adding 1 to i_qss.

The quasi-steady radius, r_qss, as indicated in figure 3, and the black hole accretion rate, |$\dot{M}_{\rm BH}\equiv |\dot{M}_{\rm in}(r=r_{\rm in})|$|⁠, are listed in table 2. We understand from table 2 that a quasi-steady state is realized inside (1–2) × 10² r_S, and that |$\dot{M}_{\rm BH}$| exceeds several tens of L_Edd/c², meaning that the super-Eddington accretion flow is actually occurring (see, e.g., Watarai et al. 2001).

Table 2.

Net mass accretion rates in quasi-steady regions.*

Model	|$\dot{M}_{\rm BH}$| [LEdd/c2]	rqss[rS]
(a11)	1.4 × 102	1.2 × 102
(a12)	2.8 × 102	2.0 × 102
(a13)	7.9 × 102	1.6 × 102
(a14)	8.3 × 102	1.6 × 102
(a15)	5.5 × 103	1.6 × 102
(a41)	1.4 × 102	9.9 × 101
(a42)	2.7 × 102	1.9 × 102
(a43)	7.5 × 102	1.6 × 102
(a44)	9.4 × 102	1.5 × 102
(a45)	5.5 × 103	1.6 × 102
(a71)	1.4 × 102	9.6 × 102
(a72)	2.7 × 102	2.0 × 102
(a73)	8.0 × 102	1.6 × 102
(a74)	9.3 × 102	1.5 × 102
(a75)	5.6 × 103	1.6 × 102

*The time-averaged mass accretion rate onto the black hole, |${\dot{M}_{\rm BH}}$|⁠, and quasi-steady radius, r_qss, inside which the quasi-steady state is realized (see figure 3). Note that |$\dot{M}_{\rm BH}\equiv |\dot{M}_{\rm in}(r=r_{\rm in})|$|⁠.

3.3 Scaling relations of the flow structure

The rather systematic variations found in the panels of figures 1 and 2 indicate the existence of simple scaling laws for the functional dependences of the density and temperature distributions on M_BH and |$\dot{M}_{\rm BH}$|⁠. To demonstrate that this is really the case, we plot in figure 4 the radial distributions of mass density ρ and gas temperature T_gas on the equatorial plane for several models.

Fig. 4.

[Top] Radial distributions of the mass density ρ (left panels) and temperature T_gas (right panels) on the equatorial plane for various black hole masses of M_BH = 10 (red), 10⁴ (green), and 10⁷ M|$_{\odot}$| (blue), but for a fixed |$\dot{M}_{\rm input}=10^{3}\, L_{\rm Edd}/c^{2}$|⁠. [Bottom] As the top panels for a variety of mass injection rates of |$\dot{M}_{\rm input}=300$| (red), 10³ (green), 5 × 10³ (blue), 10⁴ (yellow), and 10⁵ L_Edd/c² (purple) for a fixed M_BH = 10 M|$_{\odot}$|⁠. (Color online)

First of all, we notice that each line is nearly straight inside the quasi-steady radius, meaning that density and temperature are power-law functions of radius there. Next, we find a roughly constant interval between each line. We thus expect the following universal scaling relations to hold for any physical quantities, f:

\begin{eqnarray} f&=&A_{f}\times \left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{a}\left(\frac{\dot{M}_{\rm BH}}{L_{\rm Edd}/c^{2}}\right)^{b}\left(\frac{r}{r_{\rm S}}\right)^{c}, \end{eqnarray}

(18)

where A_f, a, b, and c are numerical constants that depend on the physical quantity f but are independent of M_BH, |$\dot{M}_{\rm BH}$|⁠, and r. The best-fit values on the equatorial plane are:

\begin{eqnarray} \rho &=& (9.08\pm 1.25) \times 10^{-6}\ [{\rm g \, cm^{-3}}] \nonumber \\ &&\times \left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{-1.00} \left(\frac{\dot{M}_{\rm BH}}{L_{\rm Edd}/c^{2}}\right)^{1.04} \left(\frac{r}{r_{\rm S}}\right)^{-0.73}, \end{eqnarray}

(19)

\begin{eqnarray} T_{\rm gas}&=& (3.85\pm 0.33) \times 10^{7}\ [{\rm K}] \nonumber \\ &&\times \left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{-0.24} \left(\frac{\dot{M}_{\rm BH}}{L_{\rm Edd}/c^{2}}\right)^{0.24} \left(\frac{r}{r_{\rm S}}\right)^{-0.54}, \end{eqnarray}

(20)

\begin{eqnarray} E_{0}&=&(2.36 \pm 0.14) \times 10^{15}\ [{\rm erg \, cm^{-3}}] \nonumber \\ &&\times \left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{-1.00} \left(\frac{\dot{M}_{\rm BH}}{L_{\rm Edd}/c^{2}}\right)^{1.02} \left(\frac{r}{r_{\rm S}}\right)^{-1.73}, \end{eqnarray}

(21)

\begin{eqnarray} v_{\rm r}&=& (-0.36\pm 0.01)\ [c]\nonumber \\ &&\times \left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{0.00} \left(\frac{\dot{M}_{\rm BH}}{L_{\rm Edd}/c^{2}}\right)^{0.02} \left(\frac{r}{r_{\rm S}}\right)^{-1.11}, \end{eqnarray}

(22)

\begin{eqnarray} v_{\phi }&=& (0.81\pm 0.02)\ [c]\nonumber \\ &&\times \left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{0.00} \left(\frac{\dot{M}_{\rm BH}}{L_{\rm Edd}/c^{2}}\right)^{0.01} \left(\frac{r}{r_{\rm S}}\right)^{-0.60}. \end{eqnarray}

(23)

We also calculate the standard deviations of a, b, and c, confirming that they are sufficiently small, much less than unity. In the next section we will compare these scaling laws with those for the slim disk model.

3.4 The effective temperature

In this subsection we investigate the functional dependence of a more directly observable quantity, i.e., the effective temperature, T_eff. To obtain the effective temperature distributions, we solve the gray radiative transfer equation with isotropic scattering:

\begin{eqnarray} \mu \frac{dI}{dz}&=&\frac{1}{4\pi }\epsilon _{\rm ff} - (\alpha _{\rm ff} + \rho \kappa _{\rm es}) I + \rho \kappa _{\rm es}J. \end{eqnarray}

(24)

Here, I is the specific intensity, ε_ff = 1.4 × 10⁻²⁷T^1/2(ρ/m_p)² [erg s⁻¹ cm⁻³] (m_p is the proton mass) is the emissivity, α_ff = 1.7 × 10⁻²⁵T^−7/2(ρ/m_p)² [cm⁻¹] is the absorption coefficient, κ_es = σ_T/m_p, σ_T is the Thomson cross-section, J = (1/4π)∫IdΩ = cE₀/(4π) is the mean intensity, and μ is the direction cosine. In the present study we fix μ ≡ 1 for simplicity.

We solve equation (24) numerically in the |$z$|-direction at a fixed cylindrical radius, R = rsin θ. Using the value E₀(∝J) calculated from 2D RHD code, the solution is

\begin{eqnarray} I(\tau )&=&\int _{\tau }^{\tau (-z_{\rm max})}S(\tau ^{\prime })e^{-(\tau ^{\prime }-\tau )}d\tau ^{\prime }, \end{eqnarray}

(25)

where

\begin{eqnarray} S(\tau )&\equiv &\frac{1}{4\pi } \frac{\epsilon _{\rm ff}+\rho \kappa _{\rm es}cE_{0}}{\alpha _{\rm ff}+\rho \kappa _{\rm es}} \end{eqnarray}

(26)

and

\begin{eqnarray} \tau (z)&\equiv &\int ^{z_{\rm max}}_{z}(\alpha _{\rm ff}+\rho \kappa _{\rm es})dz. \end{eqnarray}

(27)

We set |$z$|_max to be the outer boundary of the simulation box; that is,

\begin{equation} z_{\rm max}\equiv \sqrt{r_{\rm out}^{2}-R^{2}}. \end{equation}

(28)

Figure 5 shows the solution in equation (25) at radius R = 10 r_S for Models (a12) and (a13). Let us examine the case of Model (a12) first (see the top panel). The specific intensity I near the equatorial plane (⁠|$z$| ∼ 0) is equal to the blackbody |$B=\sigma _{\rm SB}T_{\rm gas}^{4}/\pi$| and the mean intensity J = cE₀/4π. This is because the optical depth is large within the accretion disk. The larger |$z$| is, the lower the gas temperature T_gas becomes, and so does the specific intensity |$I=B\propto T^{4}_{\rm gas}$| at around |$z$| ∼ 1r_S–10 r_S.

Fig. 5.

The specific intensity I [blue, equation (25)], the blackbody |$B=\sigma _{\rm SB}T_{\rm gas}^{4}/\pi$| (red), and the mean intensity J = cE₀/4π (green) at radius R = 10 r_S for Models (a12) and (a13) in the top and bottom panels, respectively. The vertical black dashed line indicates the position of τ(⁠|$z$|⁠) = 1. The optical depth at |$z$| = 100 r_S in the bottom panel is τ(⁠|$z$| = 100 r_S) ∼ 2. (Color online)

In the middle region of |$z$| ∼ 10 r_S–30 r_S, the layer is marginally optically thick (τ ≳ 1) and so the specific intensity I no longer matches the blackbody intensity B. This is because a decrease of intensity by scattering of photons out of the ray (−ρκ_esI) is dominant over an increase of intensity by scattering of photons into the ray (ρκ_esJ). We confirm that the other terms (ε_ff/4π and −α_ffI) are of minor importance in equation (24). Hence, the intensity I should become weaker and weaker with an increase of |$z$| until |$z$| ∼ 30 r_S, where τ ∼ 1 holds. We can say that this layer corresponds to a photosphere.

Above the photosphere the intensity I stays roughly constant, since radiation hardly interacts with gas there. We can thus approximate the specific intensity reaching a distant observer to be that at the outer simulation boundary, |$z$| = |$z$|_max [see equation (28)]. We can then calculate the effective temperature T_eff at radius R by

\begin{equation} T_{\rm eff}(R)\sim [\pi I(z=z_{\rm max})/\sigma _{\rm SB}]^{1/4}. \end{equation}

(29)

Let us next examine the case of Model (a13)—see the bottom panel in figure 5. The specific intensity I behaves in a similar way to that of Model (a12), but the |$z$|-dependence of the intensity is not exactly the same between them. The specific intensity I is roughly constant above the photosphere in Model (a12), whereas it still decreases slowly even above the photosphere at |$z$| ∼ 300 r_S in Model (a13). (Note that the higher the mass injection rate, the larger the scale height becomes.) We expect that I will stay nearly constant above ∼3000 r_S, although this is not numerically confirmed. In this paper, therefore, we calculate T_eff(R) by inserting the intensity at |$z$| = |$z$|_max into equation (29). We should note that this T_eff(R) is likely to be overestimated. We confirm that the shape of the intensity curve (i.e., I, J, B in figure 5) looks the same at different R in Model (a12). We also confirm that the intensity curve in all models presents the same behavior as in figure 5.

We calculate the effective temperature as a function of R, T_eff(R), for various models and plot them in figure 6. As one can see in the top panel of figure 6, the effective temperature is proportional to |$M_{\rm BH}^{-1/4}$| as long as |$\dot{M}$| is moderately large, |$\dot{M} \lesssim 10^{3}\, L_{\rm Edd}/c^2$|⁠. We confirm that this relation holds for other models with different M_BH. As |$\dot{M}$| increases, however, the effective temperature (T_eff) profile obviously becomes flatter and the innermost temperature significantly drops (see the bottom panel in figure 6).

Fig. 6.

Radial profiles of the effective temperature T_eff(R) for various black hole masses, M_BH (top panel), and for various mass injection rates, |$\dot{M}_{\rm input}$| (bottom panel). Note that the effective temperature distribution for high mass injection rates of |$\dot{M}_{\rm input}=5\times 10^{3}$|⁠, 10⁴, and 10⁵(L_Edd/c²) may not be so accurate due to the limited size of the computational box (see text). (Color online)

We wish to note again that the effective temperature calculated for models with |$\dot{M}_{\rm input}=5\times 10^{3}$|⁠, 10⁴, and 10⁵ L_Edd/c² is likely to be overestimated. This is because the photosphere is not in the numerical box size at larger distances for the high mass accretion rate model. We think that the flatter profile at high accretion rates is linked to the overestimate of T_eff due to the numerical box size, an effect which increases at larger distances from the black hole.

Why does the effective temperature decrease as the mass accretion rate increases? The multi-dimensional photon-trapping effect may be a reason (Ohsuga et al. 2002). To demonstrate that this is the case, we calculate three types of radiation advection, inward, outward, and net advection rates of radiation energy, as functions of r by

\begin{eqnarray} L_{\rm in}(r)&=&\int _{4\pi }d\Omega {\rm min}\lbrace 0, r^{2}\left[F_{r}+E_{0}v_{\rm r}\right] \rbrace <0, \end{eqnarray}

(30)

\begin{eqnarray} L_{\rm out}(r)&=&\int _{4\pi }d\Omega \, {\rm max}\lbrace 0, r^{2}\left[F_{r}+E_{0}v_{\rm r}\right] \rbrace \, >0, \end{eqnarray}

(31)

and

\begin{eqnarray} L_{\rm net}&=&L_{\rm in}+L_{\rm out}. \end{eqnarray}

(32)

Figure 7 shows |L_in|, |L_out|, and |L_net| as functions of the radius. This figure clearly shows that the inward advection of radiation energy is dominant over the outward advection near the black hole, and that the higher the mass accretion rate, the larger is |L_in|. We confirm that |L_in| is about 3.1 times larger in Model (a13) than in Model (a12) at r = 10 r_S. Thus, the radiative flux emerging from the innermost part is significantly reduced as |$\dot{M}$| increases. This is just a qualitative argument and its quantitative assessment is left as future work. Discussion regarding to what extent the boundary conditions, the spatial resolution, and the computational size affect the surface temperature will also be produced in the future.

Fig. 7.

Radial profile of the inward and outward advection rates of radiation—see equations (30)–(32). In both panels the inward advection of radiation (i.e., photon trapping) dominates inside a radius of several tens of r_S. (Color online)

Finally, we give here scaling laws of the effective temperature only for the cases with |$\dot{M}_{\rm input}=300$| and 10³(L_Edd/c²), that is,

\begin{eqnarray} T_{\rm eff}=(2.93\pm 0.01) \times 10^{7}\ [{\rm K}]\left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{-0.25}\left(\frac{r}{r_{\rm S}}\right)^{-0.47} \end{eqnarray}

(33)

\begin{equation} {\rm for}\, \dot{M}_{\rm input} =300\, L_{\rm Edd}/c^{2}\, {\rm and} \end{equation}

(34)

\begin{eqnarray} T_{\rm eff}=(2.47\pm 0.01) \times 10^{7}\ [{\rm K}] \left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{-0.25}\left(\frac{r}{r_{\rm S}}\right)^{-0.44} \end{eqnarray}

(35)

for |$\dot{M}_{\rm input} =10^{3}\, L_{\rm Edd}/c^{2}$|⁠. As we will see in the next section, these functional dependences on M_BH and r are in good agreement with those of the slim disk model.

The reason why we do not show the scaling laws for models with higher mass injection rates is that the location of the photosphere is very close to the outer boundary of the calculation box so that the effective temperature calculations may not be so reliable. We need an even larger computational box size in a future study.

4 Discussion

4.1 Comparison with the slim disk model

We obtained the fitting formulas for a super-Eddington accretion disk for the first time by systematic study of the scaling relations produced by numerical simulations. It will be interesting to examine how well the one-dimensional slim disk model can reproduce our simulation results. This is a fundamental issue, but surprisingly it has been poorly investigated in past simulation studies due probably to the limited spatial resolution. We are not in a position to answer to this question.

For comparison purposes, we use the data on the equatorial plane from Watarai (2006):

\begin{eqnarray} \rho ^{\rm slim}(r)&=& 1.2 \times 10^{-5}\ [{\rm g \, cm^{-3}}]\nonumber \\ &&\times \left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{-1} \left(\frac{\dot{M}}{L_{\rm Edd}/c^{2}}\right) \left(\frac{r}{r_{\rm S}}\right)^{-3/2}, \end{eqnarray}

(36)

\begin{eqnarray} T_{\rm gas}^{\rm slim}(r)&=& 5.3 \times 10^{7}\ [{\rm K}]\nonumber \\ &&\times \left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{-1/4} \left(\frac{\dot{M}}{L_{\rm Edd}/c^{2}}\right)^{1/4} \left(\frac{r}{r_{\rm S}}\right)^{-5/8}, \end{eqnarray}

(37)

\begin{eqnarray} T_{\rm eff}^{\rm slim}(r)&=& 4.5 \times 10^{7}\ [{\rm K}] \left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{-1/4} \left(\frac{r}{r_{\rm S}}\right)^{-1/2}, \end{eqnarray}

(38)

\begin{eqnarray} v_{\rm r}^{\rm slim}(r)&=& -1.1 \times 10^{-1}\ [c] \, \left(\frac{r}{r_{\rm S}}\right)^{-1/2}, \end{eqnarray}

(39)

\begin{eqnarray} v_{\phi }^{\rm slim}(r)&=& 7.1 \times 10^{-1}\ [c] \, \left(\frac{r}{r_{\rm S}}\right)^{-1/2}. \end{eqnarray}

(40)

Here, α = 0.1 and γ = 4/3 in these equations.

Let us compare these functional dependences with those obtained by our RHD simulations—see equations (19)–(23). We soon notice that the dependences on the black hole mass, the mass accretion rate, and the radius of each physical quantity are in reasonable agreement between the two. Two important exceptions are the radial dependences of the mass density ρ and of the radial velocity |$v$|_r. In the former, in particular, we find ρ ∝ r^−1.5 according to the slim disk model while the density profile is much shallower: ρ ∝ r^−0.5. Why is it so different?

It is interesting to note in this respect that a similar discrepancy was found in a simulation study of RIAF (radiatively inefficient accretion flow). Igumenshchev and Abramowicz (1999) were the first to demonstrate by their hydrodynamic simulations that pure ADAF (advection-dominated accretion flow) appears when the α viscosity parameter is relatively large (α ∼ 0.1), leading to a steep density profile, ρ ∝ r^−1.5, while convection arises when α is small (α ∼ 0.01), giving rise to a much flatter density profile, ρ ∝ r^−0.5. The latter type of flow is sometimes called convection-dominated accretion flow (CDAF). Machida, Matsumoto, and Mineshige (2001) examined the radial density profile by performing 3D MHD simulations and confirmed the existence of large-scale circulation. The density profile is accordingly flatter: ρ ∝ r^−0.5, |$v$|_r ∝ r^−1.3. These previous studies considered the case of a low mass accretion rate (⁠|$\dot{M}_{\rm BH}\le L_{\rm Edd}/c^{2}$|⁠). Here, we stress that even if we set α = 0.1, the convection occurs in the present study with high mass accretion rate (⁠|$\dot{M}_{\rm BH}\ge L_{\rm Edd}/c^{2}$|⁠), and the radial profile of the mass density is ρ ∝ r^−0.73.

We thus checked the simulation movies of the RHD simulations, and confirmed the occurrence of large-scale circulation (or convection) within the accretion disk (see also Ohsuga et al. 2005). The two-dimensional velocity map in the R–|$z$| plane also supports the CDAF-type flow. We may thus tentatively conclude that the flatter density profile in our RHD simulation data could be the result of the convection, which is not properly considered in the slim disk model. We should then note that the density profile may depend on the adopted α value. This point needs to be checked in future radiation-MHD (R-MHD) simulations. Note that the |$v$|_r profile is determined by the quasi-steady condition, |$\dot{M} = -2\pi r v_{\rm r} \rho H \sim$| const., with H(∼ r) being the scale-height of the inflow disk.

The effective temperature profile in equations (33) and (35) also agrees well with that of the slim disk model—see equation (38). The effective temperature of the standard disk (Shakura & Sunyaev 1973) is proportional to |$T_{\rm eff}^{\rm standard}\propto r^{-3/4}$|⁠. But, when the mass accretion rate becomes higher (⁠|$\dot{M}\gtrsim L_{\rm Edd}/c^{2}$|⁠), the radial dependence of the effective temperature becomes flatter than that of the standard disk—see equation (38). This is because of the photon-trapping effects. We can also understand the behavior of the effective temperature by this relationship: |$T_{\rm eff}^{4}\propto F\propto L_{\rm Edd}/(2\pi r^{2})$| (Kato et al. 2008).

In our RHD simulations, |$T_{\rm gas} \sim T_{\rm rad}\propto E_{0}^{1/4}$| is established from equations (20) and (21), which means that the accretion disk is optically thick. This relation is consistent with one of the assumptions needed for constructing the slim disk model.

4.2 Comparison with previous simulations

In this subsection we compare our RHD simulation results with those of previous simulations to stress what is new in the present study.

Hashizume et al. (2015) performed RHD simulations using the same code used in Ohsuga et al. (2005), but the computational box was set to be larger (r_out = 5000 r_S). The mass injection rate in Hashizume et al. (2015) was |$\dot{M}_{\rm input}=10^{3}\, L_{\rm Edd}/c^{2}$|⁠, and the initial Keplerian radius of the injected gas was r_K = 100 r_S. The important difference between our study and Hashizume et al. (2015) lies in that Compton effects were not taken into account in their simulations. According to their figure 4, the net flow rate is roughly constant in radius at r ≲ 100 r_S. The outflow rate is negligible near the black hole (at r ≲ 60 r_S), while it is substantial in the outer region at r ≳ 60 r_S. Such a separation of the innermost region without outflow and the outer region with significant outflow is also observed in our simulation data (see figure 3), although the separating radius (i.e., quasi-steady radius, r_qss) is much less in their simulations. This is because the outflowing gas density becomes significantly lower when we include the effects of Compton cooling as shown in figures 1c and 1d of Kawashima et al. (2009).

Sa̧dowski et al. (2015) performed GR-R-MHD simulation of super-Eddington accretion flow onto a 10 M_⊙ black hole for various simulation parameters (black hole spin, initial magnetic field strength and configurations, etc.). According to their figure 6, a quasi-steady state is achieved inside 30 r_S, while outflow mainly emerges outside ∼10 r_S. Although the trend that the outflow hardly emerges from the black hole vicinity is consistent with our simulation results, the radial extent, in which outflow is negligible, is significantly narrower in their simulations. This difference seems to arise in the fact that they adopted a small radius for the centroid location of the initial torus (∼21 r_S).

Let us next compare our results and theirs in terms of the velocity profiles. The azimuthal velocity |$v$|_ϕ is grossly sub-Keplerian (see their figure 13), which is consistent with our result. The mass-density-weighted and azimuthally averaged radial velocity (〈|$v$|_r〉_θ) approximately obeys the relationship of 〈|$v$|_r〉_θ ∝ r⁻² at r ≲ 30 r_S (see their figure 16), which is much steeper than our results: 〈|$v$|_r〉_θ ∝ r^−1.25. This radial dependence is very close to that on the equatorial plane: |$v$|_r ∝ r^−1.11—see equation (22). (This similar radial dependence is not so surprising, since mass density is at a maximum around the equatorial plane.) To conclude, the radial dependence of the accretion velocity on the equatorial plane in Sa̧dowski et al. (2015) is very different from our results (⁠|$v$|_r ∝ r^−1.11). This discrepancy could arise due to different treatments of disk viscosity (or magnetic processes). We adopted the α-viscosity model, whereas they solved the MHD processes in the axisymmetric geometry with a sub-grid magnetic dynamo. Again, full three-dimensional R-MHD simulations are necessary to settle this issue.

4.3 Why is outflow so weak from the innermost region?

From figure 3, we understand that the inflow rate is roughly constant in radius, while the outflow rate is very small near the black hole. This feature is in good agreement with the slim-disk formulation (since the mass flow rate is assumed to be constant) but does not quantitatively agree with the previous RHD simulation results. How can we understand this?

The small outflow rate could be due to (1) low density at the launching point of outflow (i.e., the inflow surface), (2) slow outflow velocity, or (3) a combination thereof— see equation (13). To examine which is the case, we plot the gas streamlines and the radial velocity contour lines |$v$|_r = 0, 10⁻¹, 10⁻², and 10⁻³[c] for Model (a12) in figure 8. We readily understand that outflow emerges even from near the black hole (r ≲ 100 r_S) and that the outflow speed is not small: ||$v$|_r| ≳ 0.1c. When we follow each streamline near the black hole, we see that the outflow is accelerated up to nearly the speed of light. Thus, we conclude that the mass density at the flow surface should be very small near the black hole to account for the small outflow rates.

Fig. 8.

Streamlines (red lines) and loci of constant radial velocities of |$v$|_r = 10⁻¹ (blue line), 10⁻² (yellow line), 10⁻³ (purple line), and 0 [c] (green line) for Model (a12). The black arrows indicate the direction of gas flow along the streamlines. (Color online)

The mass density ρ^surf at the flow surface for Model (a12) is plotted as a function of radius in figure 9. Here, by the surface we mean the places where the radial velocity vanishes: |$v$|_r = 0. (These are the places where outflow is launched.) The best-fit line (in the log–log plot) in the range of R = 10–100 r_S is ρ^surf ∝ R^0.4. That is, density is decreasing as matter accretes. This supports that the gas density at the outflow launching site is indeed very small.

Fig. 9.

Radial distribution of mass density (red) at the surface of the inflow region and its best-fit line (green) at R = 10–100 [r_S] for Model (a12). Here, by the “surface” we mean the location of vanishing radial velocity, |$v$|_r = 0; namely, |$v$|_r < 0 (or |$v$|_r > 0) below (above) the surface (see the green line in figure 8). (Color online)

To summarize, the high-speed (⁠|$v$|_r ≳ 0.1 [c]) outflow is driven even from the innermost region, but its gas density is negligibly small, leading to a very small outflow rate compared with inflow rate.

Finally, let us comment on the radial density profiles in other studies. According to the slim disk model, mass density on the equatorial plane is expressed as ρ^slim(r) ∝ R^−3/2 from equation (36). Since the scale-height of the slim disk is H ∼ R, we expect that density at the surface is roughly proportional to the density at the equatorial plane. That is, mass density at the outflow launching site should rapidly grow inward in the slim disk model. This does not agree with our simulation study, which shows a much flatter density profile. In the GR-R-MHD simulation, by contrast, density profile seems flat, since we find roughly Σ ∝ r (see figure 10 of Sa̧dowski et al. 2015).

Fig. 10.

Snapshot of the color contour map about the mass density with the arrows of the velocity in Model (a12). The length of the arrow corresponds to the amplitude of the velocity. There are circulating motions in the accretion flow. We remove the parts of the arrows which are outside the boundaries of the figure. (Color online)

The much flatter density profile in our results is due probably to the occurrence of radial convection. It is very plausible that this occurs, since entropy increases inward (in the direction of gravity), a condition for being convectively unstable (see Narayan & Yi 1994). Note that convection is not taken into account in the slim disk formulation. It is not yet clear why convection is not so efficient in the GR-R-MHD simulation. Careful simulation work is needed.

4.4 Convection in super-Eddington accretion flow

Figure 10 shows the convection (circulating motion) in the super-Eddington accretion disk. We estimate the timescale of the convection t_conv and the radiative diffusion t_diff. The convection timescale is calculated as t_conv = D_conv/|$v$| ∼ 0.54 [s]. Here, D_conv ∼ 2π × 15[r_S] is the typical circumference of the convection, |$v$| ∼ 5.2 × 10⁸ [cm s⁻¹] is the typical velocity of the convection, at around (R, |$z$|⁠) ∼ (70 r_S, 40 r_S) in figure 10. The radiative diffusion timescale is calculated as t_diff = 3τ_eH/c ∼ 15 [s]. Here, H ∼ 80[r_S] is the scale height and τ_e ∼ 6.3 × 10² is the optical depth in scattering at the radius R ∼ 70[r_S]. We thus understand that the convection occurs because of the relation t_conv ≲ t_diff.

Let us explain why convection occurs when this inequality holds. In a disk with high mass accretion rate (⁠|$\dot{M}c^{2}\ge L_{\rm Edd}$|⁠), photon-trapping effects occur as the advection of radiation entropy. Thus, entropy increases inward; i.e., in the direction of the gravity, a condition for convective instability (Narayan & Yi 1994).

From here, we consider the criterion of the convection in the slim disk model in analogy with the ADAF in Narayan and Yi (1994). In Watarai (2006), the pressure of the slim disk model on the equatorial plane is

\begin{eqnarray} p&=&7.8\times 10^{15}\ [{\rm dyn \, cm^{-2}}]\nonumber \\ &&\times f^{-1/2}\left(\frac{\alpha }{0.1}\right)^{-1}\left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{-1} \left(\frac{\dot{M}c^{2}}{L_{\rm Edd}}\right) \left(\frac{r}{r_{\rm S}}\right)^{-5/2}. \end{eqnarray}

(41)

Here, f depends on |$\dot{M}$| and r, and is given, from equation (24) in Watarai (2006), by

\begin{eqnarray} f&=&f(\hat{r})=\frac{1}{2}\left(\hat{r}^{2}+2-\hat{r}\sqrt{\hat{r}^{2}+4}\right), \end{eqnarray}

(42)

\begin{eqnarray} \hat{r}&\equiv &D\frac{r}{r_{\rm S}}\frac{L_{\rm Edd}}{\dot{M}c^{2}}. \end{eqnarray}

(43)

Here, D is the numerical coefficient (D ∼ 2.18), and we see that f approaches unity as |$\hat{r}$| vanishes (i.e., when |$\dot{M}c^{2}/L_{\rm Edd}$| is large and r/r_S is small). Actually, ρ also depends on f and α in such a way that ρ ∝ f^−3/2α⁻¹, but we used the approximation f ∼ 1 and set α = 0.1 in equations (36)–(40).

From equation (13) in Narayan and Yi (1994), in a rotating medium the condition for dynamical convective instability is

\begin{eqnarray} N^{2}_{\rm eff}&\equiv &-\frac{1}{\rho }\frac{dp}{dr}\frac{d\ln (p^{\frac{1}{\gamma }}/\rho )}{dr} + \left(\frac{v_{\phi }}{r}\right)^{2} <0 . \end{eqnarray}

(44)

Here, γ is the specific heat ratio. When we use the p, ρ, and |$v$|_ϕ with factors f and α of the slim disk model in equations (41), (36), and (40), |$N^{2}_{\rm eff}$| becomes

\begin{eqnarray} N^{2}_{\rm eff}&=&\frac{v_{\phi ,0}^{2}}{r_{\rm S}^{2}}\left(\frac{r}{r_{\rm S}}\right)^{-3}\left[1-g(\hat{r})\right], \end{eqnarray}

(45)

\begin{eqnarray} g(\hat{r})&=&\frac{p_{0}f(\hat{r})}{\rho _{0}v_{\phi ,0}^{2}}\frac{3\gamma -1}{2\gamma }\left(\frac{5}{2}-\frac{\hat{r}}{\sqrt{\hat{r}^{2}+4}}\right)\nonumber\\ &&\times\left(\frac{2\hat{r}}{\sqrt{\hat{r}^{2}+4}}-\frac{3\gamma -5}{3\gamma -1}\right). \end{eqnarray}

(46)

Here, ρ₀ ≡ 1.2 × 10⁻⁵ [g cm⁻³], p₀ ≡ 7.8 × 10¹⁵ [dyn cm⁻²], and |$v$|_{ϕ, 0} ≡ 2.1 × 10¹⁰ [cm s⁻¹], with γ = 4/3. We note that the α coefficient of the viscosity does not appear in |$N_{\rm eff}^{2}$| when we consider the dependence of α in ρ, p.

The criterion of convection instability is satisfied in the region |$g(\hat{r}) > 1$| (i.e., |$N_{\rm eff}^{2} < 0$|⁠). We obtain |$\hat{r}\lesssim 1.83$| by solving the relation |$g(\hat{r}) > 1$| with γ = 4/3. This range can be represented as

\begin{eqnarray} r \lesssim 0.84\left(\frac{\dot{M}}{L_{\rm Edd}/c^{2}}\right)r_{\rm S}. \end{eqnarray}

(47)

This result means that the slim disk is convectively unstable within the photon-trapping radius (r ≲ r_trap).

4.5 Future issues

There are a number of future issues to be discussed. Our simulations are restricted to Newtonian dynamics, but for discussing the Blandford–Znajek (BZ) processes that are very efficient when the central black hole is rapidly spinning we definitely need global GR-RHD simulations. In addition, we had better solve the MHD processes in purely three dimensions, since angular momentum transport by the MHD processes could be a key to examining the existence (or absence) of large-scale circulation (or convection motion) and thus to constrain the radial velocity and density profiles. Such simulation studies are extremely expensive and are impossible at present. Hence, we need careful treatments. For example, we may solve the innermost part by GR-R-MHD simulations to properly solve the gas flow dynamics near the black hole and evaluate the strengths and directions of the BZ flux, and solve the outflow dynamics in a rather large simulation box by Newtonian R-MHD simulations. The latter is essential to discuss spectral formation of high-luminosity objects, such as ULXs, since outflow material can produce Compton up-scattering of the radiation from the innermost region (Kawashima et al. 2009, 2012). Possible line emission needs to be studied (see, e.g., Pinto et al. 2016), since it could contain fruitful information from the outflow material.

Finally, we need to comment on the dependence of our results on the adopted value of the initial angular momentum (or r_K). If the quasi-steady radius increases further when we increase r_K no significant outflow is launched from any radius, in contradiction with the powerful jets from some ULXs (IC342 X-1 and Holmberg II X-1; Cseh et al. 2014), and baryonic jets from a ULS (M81 ULS1; Liu et al. 2015) and SS433. The existence of powerful outflow (or jet) is also indicated by observations of ULX nebulae (e.g., Pakull & Mirioni 2003; Grisé et al. 2006; Soria et al. 2010; Cseh et al. 2012). This issue also needs further investigation.

Acknowledgments

Numerical computations were mainly carried out on a Cray XC30 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. This work is supported in part by JSPS Grant-in-Aid for Scientific Research (C) (17K05383 S.M.; 15K05036 K.O.). This research was also supported by MEXT as “Priority Issue on Post-K computer” (Elucidation of the Fundamental Laws and Evolution of the Universe) and JICFuS.

Appendix 1. A Derivations of scaling relations

1.1 Black hole mass dependence, a

In this appendix we describe how to calculate the index a from the simulation data. Let us, for example, compare the results of model 1 and model 2 with the same |$\dot{M}$| and r_K but with the different masses M₁ and M₂. The power-law relation f(M) ∝ M^a leads to

\begin{equation} \frac{f(M_{1}; r,\theta )}{f(M_{2}; r,\theta )} = \frac{M_{1}^{a}}{M_{2}^{a}}. \end{equation}

(A1)

Here, f is any physical quantity, e.g., the mass density ρ. We can then calculate a(r, θ) by

\begin{equation} a(r,\theta )\equiv \Bigg[\log \frac{f(M_{1}; r,\theta )}{f(M_{2}; r,\theta )}\Bigg]\Bigg/\Bigg[\log \left(\frac{M_{1}}{M_{2}}\right)\Bigg] \end{equation}

(A2)

at each grid point of (r, θ) by comparing the values of model 1 and model 2. Using equation (A2), we can calculate the mass dependence in the following way:

1. We adopt three values for the black hole mass: M_BH = 10, 10⁴, and 10⁷ (M_⊙). That is, there are ₃C₂ = 3 combinations of models for a fixed mass injection rate (⁠|$\dot{M}_{\rm input}$|⁠). We thus calculate three indices, a_i(r, θ) (i = 1, 2, 3) for each |$\dot{M}_{\rm input}$|⁠.

2. We average the indices a_i(r, θ) over the three combinations for each |$\dot{M}_{\rm input}$|⁠; i.e.,

   \begin{eqnarray} \langle a(r,\theta )\rangle &\equiv &\frac{1}{3}\sum _{\rm i=1}^{3}a_{\rm i}(r,\theta ). \end{eqnarray}

   (A3)

3. The index 〈a〉 for each |$\dot{M}_{\rm input}$| is calculated by averaging 〈a(r, θ = π/2)〉 over the spatial range between r = 5r_S and r_{qss, 2}. The inner boundary 5r_S is chosen to remove the effects of the inner boundary r = r_in, while the outer boundary r_{qss, 2} ≡ 50 r_S (⁠|$\dot{M}_{\rm input}=300\, L_{\rm Edd}/c^{2}$|⁠), 100 r_S (⁠|$\dot{M}_{\rm input}\ge 10^{3}\, L_{\rm Edd}/c^{2}$|⁠) is chosen to remove the outflow effects around r ∼ r_qss where the mass outflow rate becomes large; in other words, we consider that the mass inflow rate is almost independent of the radius in r ≲ r_{qss, 2}.

4. We have confirmed that the derived values of 〈a〉 for each |$\dot{M}_{\rm input}$| are rather insensitive to the |$\dot{M}_{\rm input}$| values, so we simply averaged them.

1.2 Accretion-rate dependence, b

The derivation method for the accretion-rate dependence, b, is the same as that of a but we replace M_BH by |$\dot{M}$| and a by b. Here, the results of the models with |$\dot{M}_{\rm input}=300\, L_{\rm Edd}/c^{2}$| are not used, because the initial angular momentum is different among other models. The number of combinations of models with different |$\dot{M}$| is ₄C₂ = 6.

1.3 Radial dependence, c, and coefficients, A_f

The index c is calculated by fitting to the radial profile of each physical quantity by a power-law relation, f = B_f(r/r_S)^c, with B_f being a constant. The spatial range of fitting is the same as before; namely, between r = 5 r_S and r_{qss, 2}.

The coefficients |$A_{\rho },\ldots , A_{v_{\phi }}$| in equations (19)–(23) are calculated in the following way. The coefficient B_f introduced above includes an M_BH-dependent part and an |$\dot{M}_{\rm BH}$|-dependent part. To remove such dependences, we convert B_f to |$\tilde{A}_{f}$| by

\begin{equation} \tilde{A}_{f} \equiv B_{f}\left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{-a} \left(\frac{\dot{M}_{\rm BH}}{L_{\rm Edd}/c^{2}}\right)^{-b}. \end{equation}

(A4)

Here, we use the indices a and b obtained above. Then, we calculate the coefficients |$A_{\rho },\ldots , A_{v_{\phi }}$| by averaging all |$\tilde{A}_{f}$| values.