A GLOBAL THREE-DIMENSIONAL RADIATION MAGNETO-HYDRODYNAMIC SIMULATION OF SUPER-EDDINGTON ACCRETION DISKS

The history of the accretion rate $\dot{M}$ for this simulation is shown in Figure 1. Here the accretion rate is calculated as the net mass flux through the cylinder with radius 5r_s and height 5r_s. The first 5000t_s is the initial transient phase, when gas from the torus flows toward the black hole and forms the accretion disk. After that, $\dot{M}$ fluctuates around a mean value $-22.0\dot{M}_{\rm {Edd}}$ between 5000t_s and 1.25 × 10⁴t_s because of MRI turbulence, which shows that the accretion rate through the inner boundary has reached a steady state. As Keplerian rotation period in the Paczyński–Wiita potential at radius r is

$$\begin{eqnarray} &&t_k=755t_s\left(\frac{r}{20r_s}\right)^{1/2}\left(\frac{r/r_s-1}{19}\right), \end{eqnarray} \tag{ 7 }$$

the duration of the simulation is equivalent to 16.6 orbits at 20r_s and 49.4 orbits at 10r_s.

Zoom In Zoom Out Reset image size

Time evolution of the disk structures can also be studied with the space–time diagram as shown in Figure 2. At radii 10r_s and 20r_s, we calculate the azimuthally averaged vertical profiles of ρ, T and B_ϕ at each time. History of the averaged vertical profile is the space–time diagram for each radius. After the initial 5000t_s when the disk reaches the center, the very hot temperature and low density caused by the initial conditions go away. Density and temperature profiles are self-consistently determined by the MRI turbulence. Above a certain height at each radius, which is the location of the photosphere for effective absorption opacity (κ_esκ_ff)^1/2, gas temperature increases to 10⁸ ∼ 10⁹ K rapidly. This is caused by the dissipation of the buoyantly rising magnetic field from the disk mid-plane, which is consistent with the corona found in local shearing box simulations (Jiang et al. 2014b). However, we caution that the exact values of gas temperature in this region will change if Compton scattering is included. Surface density, as well as the height of effective absorption opacity photosphere, increases with radius, while the corona region shrinks with radius. This is also consistent with the change of corona properties with surface density as studied with local shearing box simulations (Jiang et al. 2014b). Starting near the disk mid-plane, B_ϕ reverses with time and magnetic field rises up buoyantly at each radius, which causes consistent fluctuations of density and temperature. As the local rotation period increases with radius, this process takes a longer time at larger radius. This so-called butterfly diagram has been observed widely in local shearing box (Stone et al. 1996; Miller & Stone 2000; Shi et al. 2010; Davis et al. 2010; Jiang et al. 2014b) and global simulations (O'Neill et al. 2011). This phenomenon is believed to be caused by a dynamo process in MRI-generated turbulence (Brandenburg et al. 1995; Gressel 2010; Blackman 2012). However, the period of the butterfly diagram in local shearing box simulations is usually found to be ∼10 orbits (O'Neill et al. 2011; Jiang et al. 2014b). As Keplerian rotation periods at 10 and 20r_s are only 252.9 and 755t_s, respectively (Equation (7)), Figure 2 shows that B_ϕ takes longer than 10 local orbital periods to flip in this simulation, especially around 10r_s. This suggests that properties of MRI turbulence are modified due to global effects compared with local shearing box simulations.

Zoom In Zoom Out Reset image size

Snapshots of 3D density ρ and radiation energy density E_r of the disk at time 1.13 × 10⁴t_s are shown in Figure 3. Density peaks at the disk mid-plane and decreases with height rapidly. Close to the central black hole, density drops quickly with decreasing radius. For a constant accretion rate, this implies that radial velocity increases rapidly in this region. The main body of the disk shows very turbulent structures with similar spatial distributions between E_r and ρ. However, ρ in the region near the rotation axis is very small, but relatively large E_r fills up this low-density region. A valley of E_r is formed between the rotation axis and disk mid-plane. Azimuthal variations of ρ and E_r can also be seen clearly in the full 3D simulation. Detailed structures of the disk will be studied quantitatively in the following sections.

Zoom In Zoom Out Reset image size

To see which part of the disk has reached inflow equilibrium, Figure 4 shows various mass fluxes through each radius defined as

$$\begin{eqnarray} \dot{M}_{\rm sum}&=&\oint \rho {\boldsymbol {v}}\cdot d{\boldsymbol S},\nonumber \\ \dot{M}_{\rm in}&=&\int _{-L_z/2}^{L_z/2}2\pi {\rm min}(v_r,0) r\rho dz, \nonumber \\ \dot{M}_{\rm out}&=&\int _{-L_z/2}^{L_z/2}2\pi {\rm max}(v_r,0) r\rho dz, \nonumber \\ \dot{M}_{\rm z}&=&\int _{0}^{r}2\pi v_z(z=\pm L_z/2) r\rho dr. \end{eqnarray} \tag{ 8 }$$

Here, $\dot{M}_{\rm sum}$ is the total mass flux through the cylinder with radius r, $\dot{M}_{\rm in}$, and $\dot{M}_{\rm out}$ are the inward and outward mass flux along the radial direction, respectively, $\dot{M}_{z}$ is the mass flux through the vertical direction. As the time-averaged value of $\dot{M}_{\rm sum}$ is almost a constant for different radii between time 10570t_s and 12080t_s up to ∼20r_s, this part of disk has reached inflow equilibrium and will be the focus of our analysis. Figure 4 also shows that starting from ∼4r_s, there is a significant outward mass flux along the radial and vertical directions. At 20r_s, $\dot{M}_{\rm in}=3.01\dot{M}_{\rm sum}$, $\dot{M}_{\rm out}=-1.72\dot{M}_{\rm sum}$ while $\dot{M}_{z}=-0.29\dot{M}_{\rm sum}$.

Zoom In Zoom Out Reset image size

Figure 5 shows the time and azimuthally averaged distribution of $\rho,\ {\boldsymbol {v}},\ E_r,\ {\boldsymbol {F}}_r$ in the r−z plane. Consistent with the snapshot shown in Figure 3, the disk clearly has two distinct components, namely the turbulent body of the disk and a strong outflow region within ∼45° from the rotation axis. Most of the mass is concentrated near the mid-plane of the disk, where accretion happens. The outflow starts from a place well inside the electron scattering photosphere and carries the lowest density gas in the disk. However, a significant amount of radiation energy is carried along with the outflow. The streamlines pointing toward the inner boundary are probably an artifact of the cylindrical coordinate we are using. The emerging flux from the photosphere at each radius is a composition of photons generated at different radii, which completely changes the radial profiles of the radiation flux compared with the classical one zone models where the radiation flux from photosphere at each radius is only determined by the photons generated locally.

Zoom In Zoom Out Reset image size

In order for the outward moving gas seen in the simulation to be truly astrophysical outflow, the gas has to be unbound from the gravitational potential. However, with radiative diffusion, the classical Bernoulli number is no longer a constant. One lower bound estimate is to treat the radiation acceleration as an effective reduction of the gravitational acceleration and we use the following quantity to determine whether the gas is bound or not:

$$\begin{eqnarray} &&E_t=\frac{1}{2}\rho v^2+\frac{\gamma P}{\gamma -1}-E_{\rm grav}+\frac{E_r}{3}, \end{eqnarray} \tag{ 9 }$$

where E_grav = −ρϕ. The first three terms in this equation are the classical Bernoulli constant, while the last term is to account for the balance of gravity due to radiation force. We azimuthally average E_t between 10,570t_s and 12,080t_s, which is shown in Figure 6. The outflow region seen in Figure 5 does have positive E_t while the turbulent part of the disk has negative E_t. Although Figure 5 shows that the gas with negative E_t beyond 30r_s can also move outward, this is just the dynamic motion of the torus and they cannot reach infinity. They will fall back at a larger radius, which is not captured by the simulation domain. We have done another simulation with similar setup but without radiation field. The gas can have similar large-scale outward motion but the Bernoulli constant is always negative.

Zoom In Zoom Out Reset image size

When both gas and radiation pressure gradients along the radial direction are negligible, gravitational force is balanced by the centrifugal force and the disk is in Keplerian rotation. This is what usually assumed in standard thin disk model. To check this, Figure 7 shows the radial profile of density-weighted rotation velocity for the bound gas, scaled with the Keplerian value $V_k\equiv \sqrt{r({\boldsymbol {\nabla }}\phi)_r}$. Beyond ∼8r_s, the disk is indeed Keplerian. Within ∼8r_s, V_ϕ stars to drop below V_k because a strong radiation pressure gradient is built up in this region. At r_ISCO, V_ϕ = 0.6V_k. However, we also caution that the exact number of V_ϕ/V_k within ∼6r_s, where the flow is effectively optically thin in this simulation, may change if Compton scattering is included and inner boundary is treated properly with general relativity.

Zoom In Zoom Out Reset image size

To see what is the dominant force to balance gravity and drive the outflow, vertical components of accelerations due to radiation force, gas pressure gradient and magnetic pressure gradient are calculated as

$$\begin{eqnarray} a_{\rm rad}&=&-(\boldsymbol { S_r}({\boldsymbol {P}}))_z/\rho,\nonumber \\ a_{\rm gas}&=&-({\boldsymbol {\nabla }}P)_z/\rho,\nonumber \\ a_{\rm mag}&=&-({\boldsymbol {\nabla }}(B^2/2))_z/\rho, \nonumber \\ a_{\rm grav}&=&-({\boldsymbol {\nabla }}\phi)_z. \end{eqnarray} \tag{ 10 }$$

Vertical profiles of these accelerations at radii 10 and 20r_s are shown in Figure 8, where they are scaled with the gravitational acceleration at the disk mid-plane of each radius r: a_r0 = GM_BH/(r − r_s)². Near the disk mid-plane when z is small, the gravitational acceleration is linearly proportional to z, as assumed in local shearing box simulations. This is also the turbulent body of the disk, where gravity is almost balanced by radiation force vertically. When the vertical height becomes comparable to the radius, vertical gravitational acceleration decreases with height. This is the outflow region, where radiative acceleration is much larger than gravitational acceleration. During the transition region between the turbulent disk and outflow, the radiation acceleration can actually point toward the disk mid-plane, where magnetic pressure becomes the dominant force. This is the valley of E_r as shown in Figure 3, which can also be seen clearly in the right panel of Figure 5. As the outflow carries a significant amount of radiation energy density generated from the inner region, which is actually larger than the value of E_r from the local disk part, the direction of net radiation flux reverses when the outflow touches the turbulent disk part. The fact that magnetic pressure is dominant in this transition region means that the outflow is collimated by magnetic field (Ohsuga & Mineshige 2011), although it is driven by radiation pressure.

Zoom In Zoom Out Reset image size

The turbulent state of MRI is found to have some empirical scaling relations based on local shearing box and global isothermal simulations (Blackman et al. 2008; Guan et al. 2009; Sorathia et al. 2012), which are then used as criteria to determine how MRI is resolved (Hawley et al. 2011, 2013). This is particularly important for our simulation as it is too expensive to double the resolution to check convergence. However, we caution that the dimensionless numbers from MRI turbulence in the radiation pressure dominated case may not be the same as in the gas pressure dominated regime (Jiang et al. 2013b). As MRI turbulence is only important for the bound gas, the analysis is restricted to the part with E_t < 0 and the time average is done between 10,570 and 12,080t_s. Spatial average at each radius is done along the azimuthal and vertical directions.

The first important dimensionless number is the ratio between stress and total pressure, which corresponds to the parameter α in the standard thin and slim disk models. In the simulation, angular momentum transfer is dominated by Maxwell and Reynolds stress from the MRI turbulence, which can be directly calculated as

$$\begin{eqnarray} &&W_{r\phi }=-B_rB_{\phi }+\rho v_r\delta v_{\phi }, \end{eqnarray} \tag{ 11 }$$

where δv_ϕ is the difference between v_ϕ and the azimuthally averaged v_ϕ for the same r and z. Then the spatially averaged W_rϕ at each radius r can be scaled with the spatially averaged total pressure P + E_r/3 to get an effective α at each radius. The radial profile of the time-averaged effective α is shown in the top panel of Figure 9, while the ratio between the spatially averaged radiation and gas pressure is shown at the bottom panel of the same Figure. As radiation pressure is ∼20 times the gas pressure around 20r_s and the ratio increases to 100 around 8r_s, the total pressure is truly dominated by radiation pressure. Around 20r_s, the effective α is ∼0.04, which is a little bit larger than the number found in radiation pressure dominated local shearing box simulations with zero net vertical flux (Jiang et al. 2013a). The effective alpha increases rapidly to 0.09 around 8r_s. Within 8r_s, due to strong pressure support and significant sub-Keplerian rotation (see Figure 7), the effective α drops to ∼0.02. A similar radial profile of α is also observed in global GRMHD simulation for the pure gas pressure dominate case (Penna et al. 2013; Noble et al. 2010), although the exact values of α and the location of the peak are different in our simulation. Penna et al. (2013) attribute the variation of α to the change of the shearing rate in the turbulent region and a mean magnetic field component in the laminar flow, which is consistent with the results of shearing box simulations (Pessah et al. 2008) that show a decline in α with a decreasing shear rate.

Zoom In Zoom Out Reset image size

The ratio between Maxwell stress and magnetic pressure, α_m, is always found to be ∼0.4–0.5 in unstratified isothermal simulations (Guan et al. 2009; Sorathia et al. 2012), although a smaller number ∼0.3–0.4 is also reported for stratified isothermal simulations (Hawley et al. 2011). A similar number is also found for unstratified radiation pressure dominated shearing box simulations with net vertical flux. For the case with zero net vertical flux, α_m is ∼0.25–0.3 (Jiang et al. 2013a, 2013b). For our simulation, α_m at each radius is calculated as the ratio between the spatially averaged stress and magnetic pressure. The radial profile of the time-averaged α_m is shown in the top panel of Figure 10. Between 10 and 20r_s, α_m varies between 0.25 and 0.3. It increases rapidly to 0.75 around 4r_s, where large radial motion is formed and the rotation is significant sub-Keplerian.

Zoom In Zoom Out Reset image size

The radial profile of the ratio between time and spatially averaged radial and azimuthal components of magnetic pressure is shown at the bottom panel of Figure 10. Between ∼10–20r_s, $B_r^2/B_{\phi }^2$ varies between ∼0.1–0.15, which is also similar to the number in radiation pressure dominated local shearing box simulations (Jiang et al. 2013a), but smaller than the number found in isothermal runs (Hawley et al. 2013). Within 10r_s, the radial component of magnetic pressure increases toward the central black hole due to the rapid inflow motion. It becomes the dominant component within 4r_s.

Azimuthally and time averaged (between 10,570 and 12,080t_s) radial profiles of surface density are shown in the top panel of Figure 11, while the bottom panel of the same figure shows the radial profile of density-weighted inflow velocity, scaled with the average sound speed $\sqrt{\Pi /\Sigma }$ at each radius, where Π is the sum of vertically integrated gas pressure and zz component of the radiation pressure tensor. Within ∼8r_s, the difference between Paczyński–Wiita and Newtonian potential becomes significant. Density-weighted inflow velocity increases rapidly with decreasing radius and becomes comparable to the average sound speed near the inner boundary. At the same time, surface density decreases with decreasing radius almost linearly in this region. As the fiducial surface density Σ₀ corresponds to electron scattering optical depth 6.49 × 10³, the inner region of the disk with Σ = 0.05Σ₀ is still optically thick for electron scattering. However, it becomes effectively optically thin within ∼6r_s, which may have important implications for steep power-law state of luminous black hole X-ray binaries (Dexter & Blaes 2014). Note, however, that the radial profiles of these quantities may be modified once general relativity is used self-consistently and Compton scattering is included. Beyond 10r_s, the density-weighted inflow velocity is much smaller than the average sound speed, but both the surface density and inflow velocity change more rapidly with radius.

Zoom In Zoom Out Reset image size

The radial profile of the radiation flux F_{r, z} measured from the surfaces of the simulation box is shown in Figure 12. We find that F_{r, z} is roughly constant with radius, varying by less than a factor of ∼three throughout the inner 20r_s and peaking near ∼7r_s. This is notably flatter than what is found in the standard thin disk model (Shakura & Sunyaev 1973). The total radiation luminosity, which is the integrated vertical component of the radiation flux from the top and bottom boundaries within 20r_s is 10 L_Edd for an average accretion rate of ${\sim } 20\dot{M}_{\rm {Edd}}$ at this region, as shown in the middle panel of Figure 12. Therefore, the radiative efficiency η, which is defined as $\eta \equiv L_r/(\dot{M}c^2)$, is ∼4.5%. This is actually comparable to the efficiency of the standard thin disk model and much larger than the radiative efficiency of the slim disk model, which will be smaller by almost a factor of 10 for the same accretion rate. The outflow not only carries 29% of net accreted mass flux, but also carries significant mechanical energy. The total kinetic energy luminosity associated with the outflow measured from the top and bottom surfaces of the simulation box within 20r_s is ∼20% of the total radiation luminosity from the same region, as shown in the middle and bottom panels of Figure 12.

Zoom In Zoom Out Reset image size

The fact that the radiative efficiency from the simulation is much larger than the predicted value of slim disk model suggests that there are additional important energy transport mechanisms that are not included in the slim disk model. In this model, photons only leave the disk along the vertical directions via diffusive process, which has an average speed of c/τ (Mihalas & Mihalas 1984) for an electron scattering optical depth τ measured from the disk mid-plane to the surface of the disk. The surface density of the accretion disk in the super-Eddington regime is so large that τ varies from 4 × 10⁴ at 20r_s to 200 at the inner boundary. The inflow velocity in the slim disk model is much larger than c/τ within the photon trapping radius r_trap. Therefore photons do not have time to escape from the surfaces of the disk before they are advected toward the black hole, which is the origin of photon trapping effect in the slim disk model.

However, radiative diffusion is not the only vertical cooling mechanism in the very optically thick medium with MRI turbulence. In local shearing box simulations, fluctuations of magnetic pressure is found to be anti-correlated with the density fluctuation (Blaes et al. 2011; Jiang et al. 2013a), which is buoyant. Because photon diffusion time is much longer than the local dynamical timescale especially near the disk mid-plane, photons rise vertically with the gas and escape near the disk photosphere. The magnetic buoyancy is closely related to the butterfly diagram as shown in Figure 2. Then the average energy transport speed is determined by the turbulent motion of the fluid instead of optical depth. The advective energy transport becomes more important with increasing optical depth, when radiative diffusion is inefficient.

When radiation reaches the outflow region, which is still within the electron scattering photosphere, the photons advect outward with the outflow. Because the velocity of the outflow is so large (>0.1c), radiative diffusion is still not the dominant cooling mechanism in this region.

To assess the importance of magnetic buoyancy, we calculate the cross correlations between fluctuations of density, magnetic pressure and vertical motion of the gas along the azimuthal direction for each (r, z) at each time as

$$\begin{eqnarray} \sigma _{\rho,B^2}&=&\frac{\langle (\rho -\overline{\rho })(B^2-\overline{B^2})\rangle }{\sigma _{\rho }\sigma _{B^2}},\nonumber \\ \sigma _{V_z,B^2}&=&\frac{\langle (|v_z|-\overline{|v_z|})(B^2-\overline{B^2})\rangle }{\sigma _{v_z}\sigma _{B^2}}, \end{eqnarray} \tag{ 12 }$$

where $\sigma _{\rho },\ \sigma _{B^2}$, and $\sigma _{v_z}$ are the standard deviations while $\overline{\rho },\ \overline{B^2},\ \overline{|v_z|}$ are the mean values of ρ, B² and |v_z| along the azimuthal direction. The average 〈 · 〉 is also done along the azimuthal direction. The averaged cross correlations between times 10,570 and 12,080t_s are shown in Figure 13. It is clear that in the turbulent part of the disk, there is a strong anti-correlation between density and magnetic pressure fluctuations, while the fluctuation of vertical motion is strongly correlated with magnetic pressure fluctuation in the same region. In the outflow region, the cross correlations have opposite signs as in the turbulent part. However, energy transport in the outflow region is dominated by the mean motion of the flow, not the turbulent fluctuations.

Zoom In Zoom Out Reset image size

The average advective energy transport velocities along the vertical and radial directions at each radius in the turbulent part of the disk can be calculated as

$$\begin{eqnarray} \overline{V_{E,z}}&=&\frac{\langle v_z(E_r+P/(\gamma -1))\rangle }{\langle E_r+P/(\gamma -1)\rangle },\nonumber \\ \overline{V_{E,r}}&=&\frac{\langle v_r(E_r+P/(\gamma -1))\rangle }{\langle E_r+P/(\gamma -1)\rangle }, \end{eqnarray} \tag{ 13 }$$

where the average 〈 · 〉 is done along the vertical and azimuthal directions and only for the gas with E_t < 0 (Equation (9)). Radial profiles of $\overline{V_{E,z}}$, $\overline{V_{E,r}}$ as well as the average diffusive speed c/τ used in the slim disk model averaged between times 10,570 and 12,080t_s are shown in Figure 14. It is clear that within ∼5r_s, most of the energy is advected toward the black hole. However, beyond that, vertical energy advection speed is the dominant one. If r_trap is still defined as the radius within which photons are advected toward the black hole before they have time to escape, r_trap should be 5r_s, instead of ∼330r_s (Equation (1) of Ohsuga et al. 2002) as in the slim disk model for the accretion rate in this simulation.

Zoom In Zoom Out Reset image size

The importance of vertical advective energy transport can also be shown by considering the local heating and different cooling mechanisms as done in local shearing box simulations (Hirose et al. 2009b; Jiang et al. 2013a). The local heating rate Q₊ at each (r, z) can be approximated as

$$\begin{eqnarray} &&Q_+=1.5\Omega \overline{W_{r\phi }}, \end{eqnarray} \tag{ 14 }$$

where Ω = V_k/r is the Keplerian angular velocity and $\overline{W_{r\phi }}$ is the azimuthally averaged stress. The dominant cooling mechanisms we consider are azimuthally averaged radiation flux gradients along vertical and radial directions as

$$\begin{eqnarray} &&Q_-\equiv d({\boldsymbol {F}}_r)_z/dz+d({\boldsymbol {F}}_r)_r/dr. \end{eqnarray} \tag{ 15 }$$

As a comparison, we also calculate the gradients of fluid frame radiation flux as $d({\boldsymbol {F}}_{r,0})_z/dz,\ d({\boldsymbol {F}}_{r,0})_r/dr$. Time-averaged vertical profiles of these heating and cooling mechanisms at 20r_s are shown in Figure 15. In radiation pressure dominated disks, if radiative diffusion is the only cooling mechanism, the heating rate should be Q_c = cΩ²/κ_es (Shakura & Sunyaev 1973). As shown in Figure 15, $d({\boldsymbol {F}}_{r,0})_z/dz$ roughly tracks Q_c, which is larger than $d({\boldsymbol {F}}_{r,0})_r/dr$. The actual heating and cooling rates are much larger than Q_c in this super-Eddington flow and the additional cooling is caused by advection. The total cooling rate Q₋ roughly follows Q₊ albeit large fluctuations, as the simulation duration only corresponds to roughly one thermal time at this location. The radial radiation flux gradient $d({\boldsymbol {F}}_r)_r/dr$ is dominated by the turbulent fluctuations at this radius, which can be large in localized regions. However, when vertically averaged, there is no net inward radial advective energy flux as shown in Figure 14. Therefore, vertically integrated $d({\boldsymbol {F}}_r)_z/dz$ is the most important cooling term, which is much larger than the diffusive cooling.

Zoom In Zoom Out Reset image size

Significant vertical advection also causes more mass to be concentrated toward the disk mid-plane compared with the slim disk model as shown in Figure 2 and 3. In radiation pressure dominated disks with electron scattering as the dominated opacity, vertical density profile is not constrained by the hydrostatic equilibrium. Instead, it is related to the vertical dissipation profile (Hirose et al. 2009b; Jiang et al. 2013a). When photons are able to move away from the disk mid-plane more easily compared with diffusion, the disk mid-plane is cooler and density scale height becomes smaller.

With a recently developed numerical algorithm to solve the time-dependent radiative transfer equation, we have performed a high-resolution global 3D radiation MHD simulation, maintaining a steady state accretion rate ∼220 L_Edd/c² with inflow equilibrium up to 20r_s. Surprisingly, the total radiative luminosity emitted from the disk photosphere is ∼10 L_Edd, yielding a radiative efficiency of 4.5%. This efficiency is significantly larger than the value predicted by the slim disk model, where photons are assumed to leave the disk vertically only via the diffusive process.

With MRI turbulence, vertical advective energy transport caused by magnetic buoyancy is found to be more important than pure diffusive process, especially near the disk mid-plane where the photon mean free path is much smaller than the typical size of turbulence eddies. Photons generated deep inside the disk are advected toward the photosphere and this process is not limited by the large optical depth. This effectively increases the height of dissipation and significantly reduces the cooling timescale. The disk is also thinner than what slim disk model assumes because of this additional cooling mechanism. This is consistent with the model proposed by Socrates & Davis (2006), except that most dissipation is still inside the electron scattering photosphere in the simulation. Consequently, inflow velocity is also reduced. Rapid inflow motion only exists inside ∼8r_s, where GR effects mimicked by the Paczyński–Wiita potential become significant.

The simulation also shows strong radiation-driven outflow near the rotation axis, which was not included in the original slim disk model (Abramowicz et al. 1988). However, radiative driven outflow from the slim disks has been expected for a long time (Watarai & Fukue 1999). Figure 8 shows that at each radius, when the height z becomes comparable to radius r, gravitational acceleration starts to drop with height and becomes systematically smaller than radiation acceleration. This is roughly consistent with the outflow region shown in the 2D plane of Figure 5. Photons generated in the turbulent body of the disk first enter the outflow region and leave the disk with the outflow. As the outflow starts from a place very close to the central black hole and it picks up photons from different radii, the radiation leaving the photosphere of the disk at each radius is a composition of photons generated at different locations.

The trapping radius in the slim disk model is defined in terms of the radiative diffusion timescale. However, if advective cooling is more important than radiative diffusion, the trapping radius defined in this traditional way is irrelevant. We should compare the radial and vertical energy advection speeds to consider photon trapping effect. As shown in Figure 14, the radial advection speed only becomes larger than the vertical advection speed within ∼5r_s in this simulation, which is well inside the simulation domain.

Unlike the radiation pressure dominated thin disks (Jiang et al. 2013a), radiation pressure dominated slim disks are expected to be thermally stable (Kato et al. 1998) because of the strong radial advection in the standard slim disk model. Beyond ∼8r_s in our simulation, the disk is still radiation pressure dominated but radial advection is much weaker than vertical advection, which can also stabilize the disk in principle. Because the duration of the simulation is only about one thermal time within 20r_s, this simulation cannot give a definite answer on the thermal stability of super-Eddington disks. However, this will be a focus of future work.

Super-Eddington accretion disks for non-spinning black holes have been studied extensively with 2D axisymmetric radiation hydro or MHD simulations (Ohsuga et al. 2005, 2009; Ohsuga & Mineshige 2011, 2013; Sa¸dowski et al. 2014). Approximate numerical algorithms for radiative transfer such as FLD and M1 closure are used in these calculations. These 2D simulations also find a strong radiation-driven outflow in a funnel region near the rotation axis, which is similar to the azimuthally averaged spatial structures of our 3D simulation shown in Figure 5.

However, our 3D simulation differs from previous 2D calculations in many important aspects. In contrast to our results, radiative efficiencies reported from these calculations are lower, consistent with slim disk model predictions. This is because radiative diffusion is still the dominant cooling mechanism in these 2D calculations. Since the vertical advective energy transport found in our simulation is driven by magnetic buoyancy, it is not surprising that this transport is absent in hydro simulations with a parameterized α viscosity. Even for the 2D MHD simulations, it is well known that a self-contained dynamo cannot operate in 2D because of the anti-dynamo theorem. The transient turbulence in 2D is dominated by the channel solution before it dies away, in contrast to the non-axisymmetric turbulence in 3D (Hawley et al. 1995). The butterfly diagram, as well as the associated magnetic buoyancy, is also different in 2D compared with 3D. In fact, the anti-correlation between density and magnetic pressure fluctuations shown in Figure 13 will be zero if axisymmetry is assumed. It is notable, however, that the radiative efficiency reported by a recent 3D GR radiation MHD simulation using M1 closure for a rapid spinning black hole with a similar accretion rate (McKinney et al. 2014) is also much smaller than what we find here. The discrepancy requires further investigation.

Because previous 2D and 3D simulations adopt either FLD or M1 as approximate radiative transfer algorithms, the results can also differ from ours where the full radiative transfer equation is solved, particularly when radiation forces dominate the dynamics. A small change of the direction of radiation flux can change the flow structures significantly. Since the radiation field in the outflow is composed of photons originating from different radii, M1 will merge these photons into a single beam near the photosphere (Sa¸dowski et al. 2014), which may affect the collimation of the radiation-driven outflow. As radiation flux in FLD points toward any gradient of radiation energy density, the valley of E_r shown in the right panel of Figure 5 will likely be changed with FLD.

Stratified local shearing box simulations (Turner 2004; Hirose et al. 2009a, 2009b; Blaes et al. 2011; Jiang et al. 2013a), which focus on the region near the disk mid-plane, share many similar properties with the turbulent body of the disk in our simulation. As discussed in Section 4.5, the dimensionless numbers, such as α, α_m, and $B^2_r/B^2_{\phi }$, measured in the bound gas are similar to the values reported by Hawley et al. (2013) and Jiang et al. (2013a). Vertical structures of the disk at each radius, such as the radiation pressure supported optically thick part and strong corona above the photosphere for effective absorption opacity, are also consistent with local shearing box simulations (Jiang et al. 2013a, 2014b). Significant advective energy transport along the vertical direction is also found in radiation pressure dominated local shearing box simulations (Jiang et al. 2013a).

Global effects that cannot be captured in local shearing box models show up when the height is comparable to the radius and at the inner region close to the central black hole. The vertical component of gravitational acceleration due to the central black hole is always assumed to be proportional to z in local shearing box models. This is a good approximation near the disk mid-plane but grossly overestimates the gravitational acceleration when z is comparable and larger than the radius r. Since this is the region where outflow forms, local shearing box models would significantly underestimate outflow rates. Within ∼8r_s, significant inflow motion develops and the rotation becomes sub-Keplerian, which also cannot be included in local shearing box models.

A common feature of spectra from ULXs is the turnover in the ∼2–10 keV range (Gladstone et al. 2009; Feng & Soria 2011; Kawashima et al. 2012), which suggests the existence of an emission component with similar temperature. Since Compton scattering has not been accounted for in this simulation, we have made no attempt to post-process our simulations to generate detailed spectral predictions. Nevertheless, the radiation flux shown in the top panel of Figure 12 can give a characteristic temperature of the radiation field from the simulation as T_f ≡ (F_{r, z}/(a_rc))^1/4. This effective temperature varies from 0.75 to 0.92 keV within 20r_s in this simulation. Therefore, spectra based on color-corrected blackbodies with this range of effective temperature should be broadly consistent with the observed spectra of ULXs.

Moderate to large degrees of geometric beaming have previously been invoked to explain ULX luminosities (e.g., King et al. 2001; King 2009). The angular distribution of the directly calculated specific intensities indicate how much beaming is present in our simulation. Our calculation neglects general relativistic effects, assuming that the specific intensities with the same angle ${\boldsymbol {n}}$ from different locations of the disk will be parallel rays at infinity. We spatially average I from the top and bottom boundaries of the disk within 20r_s for each angle. We average intensities with different azimuthal angles to give the distribution of I with respect to the polar angles θ, which is shown in Figure 16. Although we only have four polar angles, Figure 16 shows that the intensity peaks around θ = 40°, but is fairly isotropic when compared to the large beaming factors typically invoked to explain ULXs. It is possible that the obscuration may be underestimated for higher inclinations because the true photosphere is outside our domain for r ≳ 20r_s, but the opening angle of our funnel is about ∼30° so the beaming is very unlikely to be larger than a factor of a few, even if additional obscuration would be present. These results are consistent with constraints from emission line nebulae, which suggest emission is approximately isotropic in many ULXs (e.g., Pakull & Mirioni 2002; Moon et al. 2011).

Zoom In Zoom Out Reset image size

If super-Eddington accretion disks based on a standard slim disk model are used to explain ULXs without adopting a large beaming factor, a very large mass accretion rate $\dot{M}$ is required to generate the observed super-Eddington luminosity because of the low radiative efficiency. However, even when thin disk efficiencies are assumed, mass transfer in stellar-mass black hole binaries is only barely able to explain the observed ULX populations when sophisticated modeling of the binary evolution is carried out (Rappaport et al. 2005). This strongly suggests that one needs both high radiative efficiency and super-Eddington luminosities to explain the ULX population. Vertical radiative advection naturally provides the necessary high efficiencies. We note that an alternative model that would provide super-Eddington fluxes and high efficiencies was proposed by Begelman (2002). This model relies on non-linear evolution of photon bubbles to produce low density channels that allow more rapid photon diffusion. However, we do not see any evidence for such structures or enhancement of photon diffusion in our simulation. This is perhaps not surprising, as one might imagine that such structures have difficulty forming in turbulent MHD flows. Even in higher resolution local shearing box simulations, no evidence of photon bubbles has ever been found (Turner et al. 2005; Blaes et al. 2007; Tao & Blaes 2011). Hence, there is no evidence from numerical simulations that photon bubbles play a significant role in super-Eddington accretion flows.

If mass estimates based on emission line widths are correct, the majority of observed quasars are accreting below or, at most, slightly above the Eddington limit (e.g., Kollmeier et al. 2006). Nevertheless, our results may still be relevant to earlier stages of black hole growth, particularly at high redshift, where sustained super-Eddington accretion may be necessary to explain the large masses inferred in the early universe (Madau et al. 2014; Volonteri & Silk 2014).

Although the total radiative luminosity is as large as 10 L_Edd for an accretion rate ∼220 L_Edd/c², the radiation does not halt accretion. The outflowing gas and photons leave the system through the low density funnel near the rotation axis while mass is accreted in the plane of the disk. Our simulation implies that the super-Eddington accretion disks onto a non-spinning black hole can convert ∼4%–5% of the rest mass energy to radiative and mechanical energies with a ratio ∼5:1. This not only allows for rapid black hole growth, but also implies that there may be significant levels of both radiative and mechanical AGN feedback. Such feedback may play a critical role in the growth of black holes in the early universe and the evolution of the host galaxies (Ciotti & Ostriker 2007; Ciotti et al. 2010; Choi et al. 2012).

The simulation can be improved in many aspects. General relativity effects are approximated with the Paczyński–Wiita potential in this simulation. Although this pesudo-Newtonian potential captures several features of non-spinning black holes, studying the effects of black hole spin and jet formation require GR radiation MHD simulations. Extending our radiative transfer algorithm to GR and repeating this simulation for a spinning black hole are the next step.

The Compton process is not included in the current simulation, even though it is crucial to determine the gas temperature in the corona region. The spectrum of the radiation field from the disk will also be significantly affected by the Compton process (Kawashima et al. 2012; Schnittman et al. 2013). Compton cooling will not change the region of the disk near the disk mid-plane, where gas and radiation are in thermal equilibrium, but it may significantly alter the dynamics and thermodynamics of the coronal and funnel regions of the flow. Adding Compton process to the time-dependent radiative transfer equation for specific intensities is not as straightforward as doing this for radiation moment equations, but we are making progress on schemes to include Compton scattering in future simulations.

The simulation only reaches inflow equilibrium up to 20r_s due to limited computer power. The small dynamic range makes it hard to see the radial profile of radiation flux, which is an important observable quantity. As this simulation already shows, the existence of radiation-driven outflow makes the radiation flux at each radius be quite different from predictions of one zone models. With improved computer power and more efficient code, extending this simulation to a larger radial range will be done.